B.5 Propriétés du matériau. Modèle VEEP - BB

Utilisé par les éléments 2D et 3D (NPAR(1) = 2 ou 3).
On a pour ce modèle MODEL=NPAR(15)=4 $\longrightarrow$ loi viscoplastique-élastoplastique "béton bitumineux"
Le modèle 4 est un modèle nonlinéaire.
le modèle 4 est une loi destinée à la recherche. L'utilisateur devra vérifier qu'elle correspond à ses besoins avant tout usage.
Le programme prend NCON=42 et JDETAT=43.

7 lignes de données

lecture dans EFNL2 et EFNL3
	READ (LECG,1010) (PROP(j,n), j=1,NCON)
1010	FORMAT(8F10.0)
	READ (LECG,1070) (INIPEL(j), j=1,4)
1070	FORMAT(4I5)

note	colonnes	variable	définition
ligne 1:
1	1-10	PROP(1,n)	(XKI) module de compressibilité isotrope sous la pression de référence PREF(=PROP(36)) ($K_i$)
	11-20	PROP(2,n)	(XGI) module de cisaillement sous la pression de référence PREF(=PROP(36)) ($G_i$)
	21-30	PROP(3,n)	(XN) exposant de la loi élastique ($n$)
	31-40	PROP(4,n)	(PHI) angle de friction plastique parfait ($\phi$) (en degrés)
	41-50	PROP(5,n)	(PSI) angle de dilatance (en degrés) ($\psi$)
	51-60	PROP(6,n)	(BETA) coefficient de la loi d'état critique, qui caractérise l'écrouissage en densité. ($\beta$)
	61-70	PROP(7,n)	(PCI) pression critique initiale ($p_{ci}$)
	71-80	PROP(8,n)	(A) coefficient de la loi d'écrouissage des mécanismes déviatoitres ($a_m$)
ligne 2:
2	1-10	PROP(9,n)	(B) coefficient d'enchevetrement utilisé dans l'expression des seuils en déviatoire ($b$)
	11-20	PROP(10,n)	(ACYC) coefficient $a$ cyclique ($a_{cyc}$)
	21-30	PROP(11,n)	(ALFA) coefficient pour le domaine de comportement(dilatance) ($\alpha$)
	31-40	PROP(12,n)	(RAYELA) rayon du domaine élastique ($r_{ela}$) $0 < r_{ela} < r_{hys} < r_{mbl} \leq 1$
	41-50	PROP(13,n)	(RAYHYS) rayon du domaine hystérétique ($r_{hys}$ ) $0 < r_{hys} \leq 1$
	51-60	PROP(14,n)	(RAYMBL) rayon du domaine mobilisé ($r_{mbl}$) $0 < r_{mbl} \leq 1$
	61-70	PROP(15,n)	(C) coefficient de la loi d'écrouissage du mécanisme isotrope ($c$)
	71-80	PROP(16,n)	(D) coefficient du mécanisme isotrope reliant la droite d'état critique et la droite de consolidation isotrope ($d$)
ligne 3:
	1-10	PROP(17,n)	(CCYC) coefficient $c$ en cyclique ($c_{cyc}$)
	11-20	PROP(18,n)	(DLTELA) rayon du domaine élastique du mécanisme isotrope ($r_4$)
	21-30	PROP(19,n)	(XKIMIN) coefficient de la loi d'écrouissage isotrope valeur usuelle = 0
	31-40	PROP(20,n)	(XM) coefficient du domaine de comportement ($m$)
	41-50	PROP(21,n)	(FACINC) facteur pour la taille des sous-incréments valeur usuelle = 0.2
3	51-60	PROP(22,n)	(XKAUX) module de compressibilité isotrope pour l'opérateur auxiliaire $\texttt{XKAUX} \geq \texttt{XKI}$
3	61-70	PROP(23,n)	(XGAUX) module de cisaillement pour l'opérateur auxiliaire $\texttt{XGAUX} \geq \texttt{XGI}$
	71-80	PROP(24,n)	(IECOUL) indicateur du type de loi d'écoulement valeurs possibles :$\pm 1$,$\pm3$ valeur usuelle = 1
ligne 4:
	1-10	PROP(25,n)	(INCMAX) nombre maximum de sous-incréments autorisés
	11-20	PROP(26,n)	(AKS(1)) perméabilité saturée à la porosité POROS0 dans la direction $o_x$ (inutile en 2D)
	21-30	PROP(27,n)	(AKS(2)) perméabilité saturée à la porosité POROS0 dans la direction $o_y$
	31-40	PROP(28,n)	(AKS(3)) perméabilité saturée a la porosité POROS0 dans la direction $o_z$
	41-50	PROP(29,n)	(AK0) coefficient de poussée des terres initial ($K_0$)
	51-60	PROP(30,n)	= 0 ; PROP(30) sert à stocker la valeur compactée de l'initialisation des mécanismes, elle est calculée en fonction des valeurs de INIPEL.
4	61-70	PROP(31,n)	(IRTNU) numéro du niveau d'eau qui traverse les éléments du groupe en statique et en mécanique pure.
	71-80	PROP(32,n)	=0 (non utilisé actuellement)
ligne 5:
	1-10	PROP(33,n)	=0 (non utilisé actuellement)
5	11-20	PROP(34,n)	(FACNAP) coefficient d'une éventuelle nappe phréatique initiale
5	21-30	PROP(35,n)	(ZNAPPE) altitude de la nappe phréatique
	31-40	PROP(36,n)	(VIS) coefficient de viscosité du bitume pour un taux de cisaillement de $1 \, \frac{m}{t}$ (où $t$ est l'unité de temps utilisée)
	41-50	PROP(37,n)	(EN) exposant de la loi de comportement du fluide ($n_v$)
	51-60	PROP(38,n)	(XKF) module de compressibilité isotrope du fluide ($K_f$)
	61-70	PROP(39,n)	(XGF) module de cisaillement du fluide ($G_f$)
	71-80	PROP(40,n)	(VFFLU) pourcentage de vide remplis par le fluide
ligne 6:
	1-10	PROP(41,n)	(DSATF) degré de saturation du fluide visqueux
	11-20	PROP(42,n)	(PREF) pression de référence par défaut mise à $-1.\times10^{6}$
ligne 7:
6	1-5	INIPEL(1)	initial. mec. déviateur 1, plan $yz$
	6-10	INIPEL(2)	initial. mec. déviateur 2, plan $zx$
	11-15	INIPEL(3)	initial. mec. déviateur 3, plan $xy$
	16-20	INIPEL(4)	initialisation du mec. isotrope

notes

1. La loi VEEP-BB est une loi viscoélastique-élastoplastique cyclique à 4 mécanismes pour les matériaux constitués d'une phase solide (élastoplastique) et d'une phase fluide (viscoélastique, i.e. béton bitumineux), basée sur les concepts suivants:
   

   Pour la phase solide :

   (a)

   obéit au principe d'état critique modifié pour Mohr-Coulomb

   (b)

   la dissipation plastique est divisée en 4 mécanismes:

    $$

   3mécanismes en déformations planes dans 3 plans orthotropes coïncidant avec les axes du repère. La rotation de ces plans due à l'écrouissage n'est pas prévue dans cette loi de comportement.
   La surface de charge pour chaque mécanisme déviatoire est définie par:

   $f_k=q_k-\sin \phi p_k \left( 1-b \log \frac{ppp}{p_c} \right) r_k$

   avec:

   • $p_c=p_{ci} \exp b \varepsilon_v^p$ (écrouissage en densité)
   • $ppp=p_k$ si IECOUL>0
   • $ppp=p$ si IECOUL<0

   et: $\partial_t r_k = \lambda \frac{\left( 1-r_k\right)^2}{a_k}$
   où: $\lambda = \frac{\left( \frac{f_k}{f_0} \right)^{n_{vis}}}{\mu_{dev}}$
   $a_k = a_{cyc} + \left( a_m - a_{cyc} \right) \alpha_k$
   $\alpha_k = \texttt{PROP(11)} \left( \frac{r_k-r_{hys}}{r_{mbl}-r_{hys}} \right)^m$
   Le choix de PROP(11)=ALFA est donc très important au niveau de la dilatance.

   • PROP(11)=0, dilatance inhibée
   • PROP(11)>0, dilatance

   La loi de dilatance peut être choisie parmi les 3 hypothèses suivantes avec la donnée de JECOUL (=ABS(IECOUL) ) :

   • si JECOUL=1, $\frac{\partial_t \varepsilon_v^{p_k}}{\partial_t \varepsilon^{p_k}} = \sin \psi - \frac{q_k}{p_k} \cos \theta$
     où $\theta$ est l'angle des contraintes et des incréments de déformation plastique (=0 en monotone);
   • si JECOUL=2, $\frac{\partial_t \varepsilon_v^{p_k}}{\partial_t \varepsilon^{p_k}} = \sin \psi - \frac{q_k}{p_k}$ (Hypothèse de Roscoe)
   • si JECOUL=3, $\frac{\partial_t \varepsilon_v^{p_k}}{\partial_t \varepsilon^{p_k}} = \sin \psi$ (Hypothèse de Mohr-Coulomb)

   Les 2 premières hypothèses ne diffèrent qu'en cyclique. La première est proposée par Hujeux, la seconde par Lassoudière.

    $$

   1mécanisme en compression isotrope:

   • La surface de charge est définie par: $f_c=p-d p_c r_4$
   • La loi d'écrouissage est définie par : $\partial_t r_4 = \lambda_c \frac{\left( 1 - r_4 \right)^2}{c}$
     avec $\lambda_c = \left( \frac{f_c}{f_0} \right)^{\mu_{iso}}$

   Pour la phase fluide:
   La phase fluide obéit à la loi de Maxwell Norton:

   $\partial_t \varepsilon^p = \partial_t \frac{q}{3G_f} +A_f \left( \frac{q}{p_a} \right)^n v$

   qui est résolue par la méthode de Runge-Kutta.

   (c)

   Le chargement cyclique est pris en compte en appliquant, pour les 3 mécanismes en déviateur et pour le mécanisme isotrope, le concept de module plastique aux seuls modules d'écrouissage intrinsèques de chaque mécanisme, ainsi que la notion de domaine de comportement pour les mécanismes déviatoires.

   (d)

   La partie élastique a un comportement non linéaire isotrope, de la forme:

   • $K=K_i \left( \frac{p}{p_{ref}} \right)^n$ (module de compressibilité isotrope)
   • $G = G_i \left( \frac{p}{p_{ref}} \right)^n$ (module de cisaillement)

   où $p_{ref}$ est la pression de référence et $p$ est la contrainte isotrope

2. Les valeurs limites de $b$ permettent d'obtenir les seuils Cam-Clay et Mohr-Coulomb:
   • $b = 0$, seuil de plasticité Mohr-Coulomb
   • $b = 1$, seuil de plasticité Cam-Clay

3. XKAUX et XGAUX sont les paramètres utilisés pour calculer la matrice d'élasticité auxiliaire, qui sert à résoudre le système non linéaire.
   De manière générale, donner à XKAUX et XGAUX des valeurs égales ou plus grandes que les valeurs de XKI et XGI. L'augmentation de ces valeurs peut conduire à une convergence mais peut augmenter le temps CPU.

4. La prise en compte de la non saturation en eau n'est pas disponible dans cette version de la loi de comportement.
   Remarque: Dans un modèle purement mécanique, on peut vouloir modéliser la présence d'un niveau d'eau, qui traverse le matériau. La masse volumique du matériau est alors:
   • $(1-n) \rho_s$ au dessus du niveau d'eau
   • $(1-n)(\rho_s-\rho_w)$ sous le niveau d'eau (i.e. la masse déjaugée)
   Le niveau d'eau, préalablement défini dans le paragraphe "Chargements" (§4.5), est donné par son numéro IRTNU dans le tableau PROP à la place de la saturation résiduelle inutile en mécanique pure.
   Cette possibilité ne peut pas être utilisée en dynamique (cf. note (2) §8.2.2.1).

5. La prise en compte de la présence d'une nappe phréatique n'est utile que pour initialiser les pressions interstitielles PINT aux points d'intégration des éléments du groupe, par:
   PINT = FACNAP * DENW * ABS(GZ) * (ZIPT-ZNAPPE)
   où

   ZIPT est la cote du point d'intégration

   DENW est la masse volumique du fluide

   GZ est la gravité dans la direction $z$.

   Ceci suppose que les inconnues nodales aient été préalablement initialisées par ailleurs (avec ICON = 10, 12 ou 13)

6. Pour l'initialisation des mécanismes, il existe actuellement deux possibilités :
   1. INIPEL(i) = 1: le point de contrainte normalisé est sur le seuil, le programme calcule le rayon $r_i$ en monotone pour être plastiquement admissible , on ajuste donc le seuil sur les valeurs des contraintes sauf si on trouve $r < r_{ela}$ . On a donc :
      • Si $r_i \geq r_{ela}$, $\texttt{RAY(i,1)} = r_i$
      • Si $r_i < r_{ela}$, $\texttt{RAY(i,1)} = r_{ela}$ et INIPEL(i) = -1
      on est dans un domaine élastique, il ne peut y avoir de changements de sens (mécanisme cyclique enclenché) ou de chargement plastique tant que l'on reste dans ce domaine.
   2. INIPEL(i) = -1: Le programme calcule le rayon $r_i$ en monotone correspondant à l'état de contrainte:
      • Si $r_i < r_{ela}$, $\texttt{RAY(i,1)} = r_{ela}$
      • Si $r_i \geq r_{ela}$ le programme s'arrête
      Dans tous les cas où INIPEL(i)=-1, on ne pourra donc pas avoir de changement de sens initial ou de décharge. Il faut d'abord agrandir la surface monotone par une charge (i.e. il faut passer par un état INIPEL(i).
   Remarque: De manière générale, donner INIPEL(i) = 1.