B.4 Propriétés du matériau. Loi de comportement viscoplastique cyclique VPMCC

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

6 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
	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)	(SATR) saturation résiduelle ($S_r$) en statique et pour un problème couplé non saturé, ou (IRTNU)numéro du niveau d'eau qui traverse les éléments du groupe en statique et en mécanique pure.
4	71-80	PROP(32,n)	(ALFAS) coefficient de la loi $S(p)$ en statique et pour un problème couplé non saturé.
ligne 5:
4	1-10	PROP(33,n)	(DSDPMX) déterminé par le programme maximum de $\frac{d(sat)}{dp}$
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)	(VISDEV) coefficient de viscosité du mécanisme en déviateur ($\mu_{dev}$)
	41-50	PROP(37,n)	(VISISO)coefficient de viscosité du mécanisme isotrope ($\mu_{iso}$)
	51-60	PROP(38,n)	(VISTRC)coefficient de viscosité du mécanisme de traction ($\mu_{trc}$)
	61-70	PROP(39,n)	(XNVIS) exposant de la loi d'écoulement ($n_{vis}$)
	71-80	PROP(40,n)	(PREF) Pression de référence, par défaut mise à $-1.\times10^6$
ligne 6:
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 de VPMCC est une loi viscoplastique à 4 mécanismes pour les sols, basée sur les concepts suivants:
   1. obéit au principe d'état critique modifié pour Mohr-Coulomb
   2. la dissipation plastique est divisée en 4 mécanismes:
      • 3 mé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.
      • 1 mé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}}$
   3. 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.
   4. 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 variation du degré de saturation ($S$) et de la perméabilité ($K$) en fonction de la porosité ($n$) et de la pression interstitielle ($p_{int}$) est définie par les lois d'évolution suivantes:
   • $S = S_r + \frac{1 - S_r}{1+\left( \alpha_s \frac{p_{int}}{ \frac{\rho_w}{q_z} } \right)^2 }$
     où $S_r$ est la saturation résiduelle et $\alpha_s$ caractérise la courbure de la fonction $p(S)$.
   • $K(S,n) = K_{sat}(n) k_r(S)$
     où$k_r(S)$ est une fonction dépendant du degré de saturation et qui varie entre $0$ et $1$ selon la relation donnée par van Genuchten:
     $k_r(S) = \left( \frac{S-S_r}{1-S_r} \right)^3$
   La perméabilité saturée $K_{sat}$ est fonction de la porosité $n$ selon la loi suivante :
   $K_{sat}(n) = K_{sat}(n_0) \frac{n^3}{n_0^3} \frac{\left(1-n_0 \right)^2}{\left( 1-n \right)^2}$
   où $K_{sat}(n_0)$ (AKS(i)) est la perméabilité saturée à la porosité initiale $n_0$ (POROS0); elle est fonction de la viscosité du fluide.
   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 ou 12)
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.