B.17 Propriétés du matériau. Loi de comportement élastoplastique cyclique de Hujeux en non-saturé

Utilisé par les éléments 2D et 3D (NPAR(1) = 2 ou 3).
On a pour ce modèle MODEL=NPAR(15)=2 $\longrightarrow$ loi élastoplastique de Hujeux pour les sols non-saturés.
Le modèle 17 est un modèle nonlinéaire.
Le programme prend NCON=44 et JDETAT=38 en 2D et 42 en 3D.

8 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_{4ela}$)
	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)	(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)	(PREF) Pression de référence par défaut mise à $-1.\times10^6.$
	41-50	PROP(37,n)	non utilisé
	51-60	PROP(38,n)	(PDESAT) pression de désaturation dans un chemin de drainage (pression d?entrée d?air)
	61-70	PROP(39,n)	(PRESAT) pression de resaturation dans un chemin d'humidification (pression d?entrée d?air)
	71-80	PROP(40,n)	(CCAP) coefficient $k'$ reliant $p_c^{NS}$ à la contrainte capillaire
ligne 6:
	1-10	PROP(41,n)	(DDIS) paramètre $D_{10}$ de la courbe granulométrique
	11-20	PROP(42,n)	(PCL) pression de préconsolidation limite; non utilisé dans la version actuelle
	21-30	PROP(43,n)	(PENTD) pente d?évolution de la pression de désaturation; non utilisé dans la version actuelle
	31-40	PROP(44,n)	(P0) pression interstitielle à laquelle le matériau est saturé à 100%
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. Loi de comportement de Hujeux généralisée aux sols non-saturés.
   De même que pour les sols saturés, les déformations plastiques pour les matériaux partiellement saturés seront générées par les variations de la contrainte effective qui, dans ce cas, a pour expression:
   $\underline{\underline{\sigma}}' = \underline{\underline{\sigma}} - \pi_c \underline{\underline{I}}$
   où $\pi_c$ représente la contrainte capillaire. $\pi_c$ est une fonction de la pression interstitielle négative ainsi que de la densité et la granulométrie du matériau. Pour une densité et une granulométrie données, les variations de la pression capillaire en fonction de la pression négative doivent avoir les caractéristiques suivantes:
   • une continuité entre les domaines saturés et non saturés,
   • un palier de pression capillaire maximale pour les valeurs élevées de pression négative.
   On se donne une expression de la fonction $\pi_c(-p_{int})$ qui assure ces conditions:
   $\pi_c = \pi_{cmax} \tanh \frac{-p_{int}}{\pi_{cmax}}$
   où $\pi_{cmax}$ caractérise le palier de pression capillaire maximale. $\tanh$ représente la fonction tangente hyperbolique.
   L'influence de la densité et de la granulométrie sur la pression capillaire est prise en compte dans le paramètre $\pi_{cmax}$ qui s'écrit sous la forme suivante déduite du modèle microstructural:
   $\pi_{cmax} =$
   où

   $K(e)$ est une fonction de l'indice des vides $e$, déduite des valeurs calculées par le modèle pour les 4 arrangements de base par une régression parabolique:

   $K(e) = 0.32 e^2 + 4.06 e + 0.11$

   $r$ est un paramètre granulométrique du matériau. Des calculs avec le modèle microstructural montrent que $D_{10}$ est le plus représentatif du paramètre $r$.

   $T$ est la tension superficielle à l'interface eau-air.

   A part la généralisation du concept des contraintes effectives, cette loi de comportement basée sur le modèle multimécanisme de Hujeux est une loi de comportement élastoplastique 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^{NS}} \right) r_k$
        avec:
        • $ppp=p_k$ si IECOUL>0
        • $ppp=p'$ si IECOUL<0
        et: $\partial_t r_k = \lambda_k \frac{\left( 1-r_k\right)^2}{a_k} \mbox{ (écrouissage en déviateur)}$
        où: $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^{NS} r_4$
          avec $p'$ la contrainte moyenne effective qui s'écrit:
          $p' = p - \pi_c$
          où p est la contrainte moyenne totale et:
          $p_c^{NS} = p_{ci} \exp \left( \beta \varepsilon_v^p \right) - k' \pi_c \mbox{ (écrouissage en densité)}$
        • La loi d'écrouissage est définie par : $\partial_t r_4 = \lambda_c \frac{\left( 1 - r_4 \right)^2}{c}$

   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 effective généralisée 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. On suppose être non saturé dans les cas suivants:
   • $P_{int} \leq \texttt{PDESAT}$ sur un chemin de drainage ( $\Delta P_{int} < 0$)
   • $P_{int} > \texttt{PRESAT}$ sur un chemin d?humidification ( $\Delta P_{int} > 0$)
   Dans ce cas le degré de saturation est une fonction de la pression intersticielle.
   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 }$
   • $K(S) = K_r(n) \left( \frac{S-S_r}{1-S_r} \right)^3$
   où

   $S_r$ est la saturation résiduelle

   $\alpha_s$ caractérise la courbure de la fonction $p(S)$.

   $p_0$ est la pression interstitielle à laquelle le materiau se sature

   La perméabilité relative $K_r$ est fonction de la porosité $n$ (POROS) suivant la loi suivante :
   $K_r = K_{r0} \frac{n^3}{\left( 1-n \right)^2}$
   où $K_{r0}$ est la perméabilité à la porosité initiale $n_0$; 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.