r7.01.48 Modèle de comportement MCC

Résumé :

Ce document présente la formulation et l’intégration numérique du modèle de comportement élastoplastique Cam-Clay modifié dans le cadre des Matériaux Standard Généralisés. Ce modèle est dédié à représenter de façon qualitative la phénoménologie des sols granulaires, tenant compte des effets de dilatance/contractance, d’adoucissement/durcissement, ainsi que de l’état critique. Les paramètres du modèle de comportement sont au nombre total de sept. Trois définissent le comportement élastique non-linéaire isotrope ; trois représentent le domaine de réversibilité initial. Un écrouissage combiné cinématique-isotrope porté par la déformation plastique volumique est piloté par un dernier paramètre.

Les équations de comportement sont résolues suivant un principe variationnel pour lequel les conditions d’optimalité sont équivalentes à un système d’équations non-linéaires régissant les équations d’évolutions incrémentales des variables internes. Celles-ci sont résolues en employant une méthode d’intégration implicite directe en temps dans MFront.

Formulation du modèle

Le modèle Cam-Clay modifié est formulé selon la construction usuelle aux Matériaux Standard Généralisés, à partir d’un potentiel d’énergie et d’un critère de plasticité (ou d’un potentiel de dissipation). Le cas des déformations linéarisées et des conditions isothermes est supposé.

Lois d’état

Les variables d’état sont le tenseur des déformations totales \(\boldsymbol{\epsilon}\), le tenseur des déformations plastiques \(\boldsymbol{\epsilon}^p\), et une variable d’écrouissage scalaire \(\xi\). Son inteprétation sera précisée lors de l’établissement des lois d’évolution.

Le potentiel d’énergie libre \(\psi\) (densité volumique) s’écrit sous la forme :

(4225)\[\psi(\boldsymbol{\epsilon},\boldsymbol{\epsilon}^p,\xi) = \psi_e(\boldsymbol{\epsilon}-\boldsymbol{\epsilon}^p) + \psi_h(\xi)\]

où \(\psi_e\) désigne l’énergie restituée par décharge élastique et \(\psi_h\) représente l’énergie stockée par écrouissage. Leurs expressions sont indiquées ultérieurement.

L’expression de la densité volumique de dissipation intrinsèque \(D\) s’obtient comme :

(4226)\[D=\mathbf{\boldsymbol{\sigma}}:\dot{\boldsymbol{\epsilon}}-\dot{\psi} = \mathbf{\boldsymbol{\sigma}}:\dot{\boldsymbol{\epsilon}} - \left(\frac{\partial \psi}{\partial \boldsymbol{\epsilon}}: \dot{\boldsymbol{\epsilon}}+\frac{\partial \psi}{\partial \boldsymbol{\epsilon}^p}:\dot{\boldsymbol{\epsilon}}^p +\frac{\partial \psi}{\partial \xi}\dot{\xi}\right)= \boldsymbol{\sigma}:\dot{\boldsymbol{\epsilon}}^p + p_c\dot{\xi}\]

cette dernière égalité étant due au caractère non-dissipatif de la déformation totale \(\boldsymbol{\epsilon}\) (la plasticité est sans effet visqueux).

Ci-dessus, on aura également défini \(\boldsymbol{\sigma}\) le tenseur des contraintes et \(p_c\), celle-ci dénommée la pression critique. Ces deux forces thermodynamiques sont naturellement dérivées par les lois d’état suivantes :

(4227)\[\begin{split}\begin{align} \boldsymbol{\sigma} &= \frac{\partial \psi}{\partial \boldsymbol{\epsilon}} = \frac{\partial \psi_e}{\partial \boldsymbol{\epsilon}} =-\frac{\partial \psi_e}{\partial \boldsymbol{\epsilon}^p} \\ p_c &=-\frac{\partial \psi}{\partial \xi} = -\frac{\partial \psi_h}{\partial\xi} \end{align}\end{split}\]

Énergie élastique

Le potentiel d’énergie élastique \(\psi_e\), établi dans [Bour97], et mis à profit ici, prend l’expression suivante :

(4228)\[\psi_e(\boldsymbol{\epsilon}^e) = \frac{K}{\kappa^2}\left(\Upsilon(\boldsymbol{\epsilon}^e)+\kappa\mathrm{tr}(\boldsymbol{\epsilon}^e)-1\right)\qquad\text{avec}\qquad \Upsilon(\boldsymbol{\epsilon}^e) = \exp\left(-\kappa\mathrm{tr}(\boldsymbol{\epsilon}^e)+\frac{\mu}{K}(\kappa\boldsymbol{\epsilon}^e_d):(\kappa\boldsymbol{\epsilon}^e_d)\right)\]

où \(\boldsymbol{\epsilon}^e=\boldsymbol{\epsilon}-\boldsymbol{\epsilon}^p\) désigne le tenseur des déformations élastiques, ayant noté \(\boldsymbol{\epsilon}^e_d=\boldsymbol{\epsilon}^e-\cfrac{\mathrm{tr}(\boldsymbol{\epsilon}^e)}{3}\boldsymbol{I}\) son déviateur (\(\boldsymbol{I}\) tenseur identité d’ordre deux).

L’expression du tenseur des contraintes se déduit de la loi d’état (4227)-1 :

(4229)\[\boldsymbol{\sigma} = \Upsilon(\boldsymbol{\epsilon}^e)\left(-\frac{K}{\kappa}\boldsymbol{I}+2\mu\boldsymbol{\epsilon}^e_d\right)+\frac{K}{\kappa}\boldsymbol{I}\]

Dans (4227) et (4229), \(K\) et \(\mu\) sont respectivement les modules de compressibilité et de cisaillement initiaux. Le coefficient \(\kappa\geq 0\) introduit une non-linéarité par l’effet de la contrainte moyenne \(\sigma_m\) sur les modules de compressibilité et de cisaillement tangents. En particulier, cette élasticité non-linéaire généralise celle plus classique du modèle Cam-Clay modifié, uniquement sur le comportement volumique (voir [r7.01.14]_), au comportement déviatorique. Cet effet se trouve détaillé [Bacq23] sur des chargements isotropes et de cisaillement.

Remarque :

Le potentiel d’énergie élastique (4228) est bien défini pour \(\kappa=0\) puisqu’à l’aide d’un développement limité à l’ordre deux en \(\kappa\boldsymbol{\epsilon}^e\), on a :

\[\psi_e(\boldsymbol{\epsilon}^e) \underset{\|\kappa\boldsymbol{\epsilon}^e\|\ll 1}{=} \frac{K}{2}\mathrm{tr}(\boldsymbol{\epsilon}^e)^2 + \mu\boldsymbol{\epsilon}^e_d:\boldsymbol{\epsilon}^e_d + \mathcal{O}(\|\kappa\boldsymbol{\epsilon}^e\|^3)\]

Donc pour \(\kappa=0\), le potentiel d’énergie élastique (4228) prédit un comportement linéaire élastique isotrope.

Énergie stockée par écrouissage

Le potentiel d’énergie stockée par écrouissage \(\psi_h\) s’écrit :

(4230)\[\psi_h(\xi) = \frac{p_{c0}}{\beta}\left(\exp(-\beta\xi)-1\right)\]

On en déduit l’expression de la pression critique d’après la loi d’état (4227)-2 :

(4231)\[p_c = p_{c0}\exp(-\beta\xi)\]

\(\beta\geq 0\) est appelé indice d’incompressibilité plastique et \(p_{c0}>0\) est la pression critique initiale.

Lois d’évolution

L’expression du critère de plasticité est l’équation d’une ellipse dans le plan méridien des contraintes :

(4232)\[f(\boldsymbol{\sigma},p_c; \xi) = \sqrt{\left(\frac{\sigma_{eq}}{M}\right)^2 + \left(\sigma_m+p_c\right)^2} - R(\xi)\]

avec \(\sigma_m=\cfrac{\mathrm{tr}(\boldsymbol{\sigma})}{3}\) la contrainte hydrostatique et \(\sigma_{eq}=\sqrt{\cfrac{3}{2}\boldsymbol{\sigma}_d:\boldsymbol{\sigma}_d}\) la contrainte équivalente de von Mises, ayant noté \(\boldsymbol{\sigma}_d=\boldsymbol{\sigma}-\sigma_m \boldsymbol{I}\) le déviateur du tenseur des contraintes.

Dans (4232), le paramètre \(M\), appelé pente d’état critique, module le rapport des axes du domaine de réversibilité dans le plan méridien des contraintes. Géométriquement, la pression critique \(p_c\) positionne le centre de l’ellipse suivant l’axe hydrostatique, et enfin la fonction \(R(\xi)\) précise sa taille.

Remarque :

Suivant l’expression du critère (4232), il convient de noter que \(p_c\) joue le rôle de force de rappel hydrostatique. Sa variation prédit donc un écrouissage de type cinématique suivant l’axe hydrostatique. De plus, le critère est paramétré, comme le précise la notation \(f( ; \xi)\), via la fonction \(R(\xi)\). Son évolution conduit donc à prédire un écrouissage de type isotrope. Ainsi, le modèle Cam-Clay modifié est à écrouissage combiné cinématique-isotrope. Ces deux mécanismes sont tous deux pilotés par la variable d’écrouissage scalaire \(\xi\).

La règle d’écoulement du modèle Cam-Clay modifié respecte la loi de normalité, ainsi le tenseur des déformations plastiques \(\boldsymbol{\epsilon}^p\) et la variable \(\xi\) évoluent comme :

(4233)\[\begin{split}\begin{align} \dot{\boldsymbol{\epsilon}}^p &= \dot{\lambda}\frac{\partial f}{\partial \boldsymbol{\sigma}} = \frac{\cfrac{3}{2M^2}\boldsymbol{\sigma}_d + \left(\sigma_m+p_c\right)\cfrac{\boldsymbol{I}}{3}}{\sqrt{\left(\cfrac{\sigma_{eq}}{M}\right)^2 + \left(\sigma_m+p_c\right)^2}}\\ \dot{\xi} &= \dot{\lambda}\frac{\partial f}{\partial p_c} = \frac{\sigma_m+p_c}{\sqrt{\left(\cfrac{\sigma_{eq}}{M}\right)^2 + \left(\sigma_m+p_c\right)^2}} \end{align}\end{split}\]

où \(\dot{\lambda}\) est le multiplicateur plastique donné par la condition de cohérence suivante :

(4234)\[\dot{\lambda}\geq 0, \quad f\leq 0, \quad \dot{\lambda}f=0\]

À partir de la loi d’écoulement (4233), on en déduit que \(\dot{\xi}=\mathrm{tr}(\dot{\boldsymbol{\epsilon}}^p)\). La variable d’écrouissage scalaire du modèle Cam-Clay modifié est donc la déformation plastique volumique.

Remarque :

On se propose d’établir l’expression du potentiel de dissipation \(\phi\) du modèle, en procédant par la transformation de Legendre-Fenchel de la fonction indicatrice du domaine de réversibilité défini par le critère de plasticité dans (4232). On obtient :

(4235)\[\phi(\dot{\boldsymbol{\epsilon}}^p,\dot{\xi};\xi) = \sup_{f(\boldsymbol{\sigma},p_c; \xi)\leq 0}\lbrace \boldsymbol{\sigma}:\dot{\boldsymbol{\epsilon}^p}+p_c\dot{\xi}\rbrace =R(\xi)\sqrt{\left(M\dot{\epsilon}_{eq}^p\right)^2+\mathrm{tr}(\dot{\boldsymbol{\epsilon}}^p)^2} + I_{\lbrace 0\rbrace }\left(\mathrm{tr}(\dot{\boldsymbol{\epsilon}}^p)-\dot{\xi}\right)\]

avec la fonction indicatrice \(I_{\lbrace 0\rbrace }(x)=0\) si \(x=0\), \(I_{\lbrace 0\rbrace }(x)=+\infty\) sinon, et \(\dot{\epsilon}_{eq}^p=\sqrt{\cfrac{2}{3}\dot{\boldsymbol{\epsilon}}^p_d: \boldsymbol{\epsilon^p_d}}\) le taux de déformation plastique équivalent de von Mises, ayant noté \(\dot{\boldsymbol{\epsilon}}^p_d=\dot{\boldsymbol{\epsilon}}^p-\cfrac{\mathrm{tr}(\dot{\boldsymbol{\epsilon}}^p)}{3}\boldsymbol{I}\) son déviateur.

Expression de la paramétrisation à l’état

La formulation du modèle Cam-Clay modifié s’achève en précisant la définition de la fonction \(R(\xi)\), celle-ci étant responsable d’un écrouissage de type isotrope avec la déformation plastique volumique \(\xi=\mathrm{tr}(\boldsymbol{\epsilon}^p)\). Rappelons que cet écrouissage provient d’une paramétrisation à l’état dans le critère de plasticité (4232) (et dans le potentiel de dissipation (4235)). Son expression est déduite du maintien constant d’une limite en traction notée \(\sigma_0\geq 0\) quel que soit la valeur de \(\xi\). Cette condition s’exprime donc comme suit :

(4236)\[\forall \xi, \quad f(\boldsymbol{\sigma}=\sigma_0\boldsymbol{I}, p_c;\xi) = 0\]

À partir de l’expression du critère de plasticité (4232) ainsi que de l’expression de la pression critique (4231), l’égalité (4236) est vérifiée si :

(4237)\[R(\xi) = \sigma_0+p_c=\sigma_0+p_{c0}\exp(-\beta\xi)\]

Géométriquement, le domaine de réversibilité du modèle Cam-Clay modifié est donc délimité suivant l’axe hydrostatique par le segment \([-p_c-R(\xi),-p_c+R(\xi)]=[-\sigma_0-2p_c,\sigma_0]\) dans le plan méridien des contraintes. La figure suivante le représente pour différentes valeurs de déformation plastique volumique \(\xi\).

Fig. 296 Domaine de réversibilité du modèle Cam-Clay modifié dans le plan méridien des contraintes pour plusieurs valeurs de déformation plastique volumique \(\xi\). Une diminution de \(\xi\) (contractance) conduit à agrandir le domaine de réversibilité. Une augmentation de \(\xi\) (dilatance) le réduit. Tout du long, la limite en traction isotrope reste constante égale à \(\sigma_0\) du fait de vérifier (4236).

Équations de l’état critique

Dans ce paragraphe, on établit les équations de l’état critique du modèle Cam-Clay modifié. Cet état correspond à un chargement monotone en déformation pour lequel l’état des contraintes et la déformation volumique restent constants. Pour cela, le taux de déformation plastique volumique \(\dot{\xi}=\mathrm{tr}(\dot{\boldsymbol{\epsilon}}^p)\) s’annule. À partir de la loi d’écoulement (4233)-2, cette condition survient lorsque :

(4238)\[\sigma_m+p_c=0\]

La pression moyenne \(-\sigma_m\) est ainsi égale à la pression critique \(p_c\). De plus, à l’aide de l’expression du critère de plasticité (4232), on en déduit que la contrainte équivalente \(\sigma_{eq}\) est dans cette situation :

(4239)\[\sigma_{eq} = MR(\xi)\]

L’expression de la fonction \(R(\xi)\) dans (4237) permet de conclure à l’aide de (4238) :

(4240)\[\sigma_{eq} = M(\sigma_0+p_{c}) \Longrightarrow \sigma_{eq} + M(\sigma_m-\sigma_0)=0\]

En résumé, l’état critique prédit par le modèle Cam-Clay modifié est défini par le jeu d’équations suivant :

(4241)\[\begin{split}\begin{align} \dot{\boldsymbol{\sigma}}=\boldsymbol{0}\\ \mathrm{tr}(\dot{\boldsymbol{\epsilon}})=0\\ \sigma_m+p_c = 0\\ \sigma_{eq}+M(\sigma_m-\sigma_0)=0 \end{align}\end{split}\]

Il s’agit d’un état de plasticité isochore (sans variation de volume) et parfait (sans variation des contraintes). Géométriquement, l’état des contraintes à l’état critique est positionné sur la demi-droite de pente \(M\) sectionnant la surface du domaine de réversibilité de forme elliptique à son sommet d’abscisse \(-\sigma_m=p_c\).

Formulation incrémentale

La formulation du modèle Cam-Clay modifié est désormais discrétisée en temps afin de mettre en place une approche de résolution numérique s’appuyant sur un schéma d’intégration implicite des équations de comportement. Celles-ci sont obtenues suivant un principe variationnel incrémental ([Mial86], [OrSt99], etc.) expliqué ci-dessous.

Principe de minimisation incrémentale

Pour mettre en œuvre une résolution incrémentale des équations de comportement du modèle Cam-Clay modifié, on considère une discrétisation d’instants \(t_0, t_1, \dots, t_{n+1}=t_n+\Delta t\). Sur l’intervalle \([t_n,t_{n+1}]\), l’évolution des variables internes est approchée par un schéma Euler-implicite :

(4242)\[\dot{\boldsymbol{\epsilon}}^p \approx \frac{\boldsymbol{\epsilon}^p-\boldsymbol{\epsilon}^p_{n}}{\Delta t}, \quad \dot{\xi} \approx \frac{\xi-\xi_n}{\Delta t}\]

Les valeurs des deux variables internes à l’instant \(t_{n+1}\), qui sont \(\boldsymbol{\epsilon}^p_{n+1}\) et \(\xi_{n+1}\), sont recherchées en résolvant le principe variationnel local suivant :

(4243)\[(\boldsymbol{\epsilon}_{n+1}^p,\xi_{n+1}) = \underset{\boldsymbol{\epsilon}^p,\xi}{\mathrm{argmin}}\left\lbrace \psi(\boldsymbol{\epsilon}_{n+1},\boldsymbol{\epsilon}^p,\xi) + \Delta t\left\langle \phi\left(\frac{\boldsymbol{\epsilon}^p-\boldsymbol{\epsilon}^p_n}{\Delta t},\frac{\xi-\xi_n}{\Delta t};\xi(\tau)\right)\right\rangle\right\rbrace\]

Le dernier terme au membre de droite, entre crochets, désigne une approximation cohérente du potentiel de dissipation \(\phi\) moyenné sur l’intervalle \([t_n,t_{n+1}]\). Suivant les développements présentés dans [Bacq23] (voir aussi [BRAV24]), une expression licite pour cette approximation s’écrit comme :

(4244)\[\Delta t\left\langle \phi\left(\frac{\boldsymbol{\epsilon}^p-\boldsymbol{\epsilon}^p_n}{\Delta t},\frac{\xi-\xi_n}{\Delta t};\xi(\tau)\right)\right\rangle = \phi\left(\boldsymbol{\epsilon}^p-\boldsymbol{\epsilon}^p_n,\xi-\xi_n;\xi_n\right) + \frac{H_n}{2} \left(\left(M\|\boldsymbol{\epsilon}^p-\boldsymbol{\epsilon}^p_{n}\|_{eq}\right)^2 +\mathrm{tr}\left(\boldsymbol{\epsilon}^p-\boldsymbol{\epsilon}^p_n\right)^2\right)\]

En quelques mots, le premier terme au membre de droite ci-dessus s’interprète comme le potentiel de dissipation intégré sur l’intervalle \([t_n,t_{n+1}]\) en figeant la dépendance à l’état à l’instant \(t_n\), c’est-à-dire à \(\xi(\tau) = \xi_{n}\). Le second terme approche les conséquences de l’évolution de cette dépendance sur le pas d’intégration, dans le cas d’un écrouissage positif, pour lequel le module \(H_n\) est une quantité positive ou nulle estimée par la solution établie à l’instant \(t_n\). Son expression est la suivante :

(4245)\[H_n = \max\left\lbrace \frac{R'(\xi_n)}{R(\xi_n)}\left(\sigma_{m,n}+p_c(\xi_n)\right),0\right\rbrace\]

On notera que ce terme possiblement non-nul n’existe qu’en conséquence de la paramétrisation à l’état du modèle Cam-Clay modifié à partir de \(R'(\xi_n)\neq 0\), cette dépendance étant, rappelons-le, la source de l’écrouissage isotrope du modèle.

Équations d’optimalité

Le potentiel d’énergie totale incrémentale exprimé par (4243), à partir de (4244), peut se réécrire de la manière suivante :

(4246)\[(\boldsymbol{\epsilon}_{n+1}^p,\xi_{n+1}) = \underset{\boldsymbol{\epsilon}^p,\xi}{\mathrm{argmin}}\left\lbrace \psi(\boldsymbol{\epsilon}_{n+1},\boldsymbol{\epsilon}^p,\xi) + \phi\left(\boldsymbol{\epsilon}^p-\boldsymbol{\epsilon}^p_n,\xi-\xi_n;\xi_n\right) + \frac{H_n}{2} \left(\left(M\|\boldsymbol{\epsilon}^p-\boldsymbol{\epsilon}^p_{n}\|_{eq}\right)^2 +\mathrm{tr}\left(\boldsymbol{\epsilon}^p-\boldsymbol{\epsilon}^p_n\right)^2\right)\right\rbrace\]

La résolution de ce problème implique alors l’obtention des équations d’optimalité du premier ordre. Elles sont données comme tel :

(4247)\[\begin{split}\begin{align} &\boldsymbol{\sigma}_{n+1} = \frac{\partial \psi_e}{\partial \boldsymbol{\epsilon}_{n+1}}, \quad p_{c,{n+1}} =-\frac{\partial \psi_h}{\partial \xi_{n+1}} & (\text{lois d'état})\\ &f_{n+1} = \sqrt{\left(\frac{\sigma_{eq,{n+1}}}{M}\right)^2+\left(\sigma_{m,n+1}+p_{c,n+1}\right)^2}-\left(R(\xi_n)+H_n\Delta\lambda\right)& (\text{critère de plasticité})\\ &\Delta \boldsymbol{\epsilon}^p = \Delta\lambda \frac{\partial f_{n+1}}{\partial \boldsymbol{\sigma}_{n+1}},\quad \Delta \xi = \Delta\lambda \frac{\partial f_{n+1}}{\partial p_{c,{n+1}}}& (\text{lois d'écoulement})\\ &\Delta \lambda\geq 0,\quad f_{n+1}\leq 0, \quad \Delta\lambda f_{n+1}=0& (\text{condition de cohérence}) \end{align}\end{split}\]

en ayant noté par concision \(\Delta\boldsymbol{\epsilon}^p=\boldsymbol{\epsilon}^p_{n+1}-\boldsymbol{\epsilon}^p_n\), etc., les incréments de \(\boldsymbol{\epsilon}^p\), \(\xi\) et \(\lambda\) sur \([t_n,t_{n+1}]\).

Remarque :

Sur la ligne (4247)-2 portant sur la définition du critère de plasticité, le rayon du domaine de réversibilité s’exprime comme \(R(\xi_n)+H_n\Delta\lambda\), étant dans le cas général légèrement différent de \(R(\xi_{n+1})\) pour un pas de temps non nul. Cette différence résulte de la cohérence des équations de comportement comme conditions d’optimalité à la résolution du problème de minimisation incrémentale dans (4246). À l’opposé, insistons sur le fait que vouloir résoudre la deuxième ligne du système (4247), en insérant l’expression du rayon \(R(\xi_{n+1})\) en lieu de \(R(\xi_n)+H_n\Delta\lambda\), n’est en général par une condition d’optimalité d’un problème de minimisation (voir par exemple [BoLe90] à ce sujet).

Paramètres

Classification

Sept paramètres sont nécesssaires à définir le modèle Cam-Clay modifié, ceux-ci étant classifiés ci-dessous.

Tableau 70 Paramètres.

Effet	Appellation	Définition	Unité	Intervalle	Symbole
Élasticité non-linéaire isotrope	BulkModulus ShearModulus SwellingIndex	Module de compressibilité Module de cisaillement Indice de non-linéarité élastique	[Pa] [Pa] [-]	\(\mathbb{R}^{+*}\) \(\mathbb{R}^{+*}\) \(\mathbb{R}^{+}\)	\(K\) \(\mu\) \(\kappa\)
Domaine d’élasticité initial	InitCritPress CritStateSlope TensileYieldStress	Pression critique initiale Pente d’état critique Limite d’élasticité en traction isotrope	[Pa] [-] [Pa]	\(\mathbb{R}^{+*}\) \(\mathbb{R}^{+*}\) \(\mathbb{R}^{+}\)	\(p_{c0}\) \(M\) \(\sigma_0\)
Écrouissage cinématique-isotrope	IncoPlastIndex	Indice d’incompressibilité plastique	[Pa]	\(\mathbb{R}^{+}\)	\(\beta\)

Équivalence aux paramètres [r7.01.14]

On indique dans le tableau suivant l’équivalence des sept paramètres du modèle Cam-Clay modifié ici présenté à ceux nécessaires à définir le modèle CAM_CLAY [r7.01.14] [r7.01.14]_.

Remarque :

Le tableau suivant propose une identification du comportement élastique du modèle Cam-Clay modifié ici présenté au modèle CAM_CLAY [r7.01.14]_. Néanmoins, celle-ci n’est pas exacte. Dans cette dernière, l’élasticité adoptée est linéaire sur le comportement déviatorique et suppose un état de contrainte hydrostatique non nul (\(\boldsymbol{\sigma}=-\widetilde{P}_0\boldsymbol{I}\)) à déformation élastique nulle. Le lecteur est invité à parcourir [r7.01.14]_ pour obtenir plus de détails sur l’élasticité y étant formulée équation 3.5-6.

Tableau 71 Équivalence aux paramètres du modèle CAM_CLAY [r7.01.14]. Pour lever toute confusion, ces derniers sont sur-indexés par le symbole \(\widetilde{\cdot}\) pour les différencier de ceux définis dans ce document.

Relation aux paramètres [r7.01.14]	Référence [r7.01.14]
\(K\approx\widetilde{k}_0\widetilde{P}_0+\widetilde{K}_{cam}\)	équation 3.3-3
\(\mu\approx\widetilde{\mu}\)	équation 3.3-1
\(\kappa=\widetilde{k}_0\)	équation 3.3-3
\(p_{c0}=\widetilde{P}_{cr}^0\)	équation 3.3-4
\(M=\widetilde{M}\)	équation 3.2-1
\(\sigma_0=-\widetilde{P}_{trac}\)	équation 3.2-1
\(\beta=\widetilde{k}\)	équation 3.3-4

Intégration numérique

Résolution implicite des équations de comportement dans MFront

Le modèle de comportement Cam-Clay modifié est intégré de manière implicite (directive \(\texttt{@DSL Implicit}\)) via l’outil MFront [HMPS15]. Cette intégration intervient sur l’incrément de déformation totale imposée \(\Delta \boldsymbol{\epsilon}\) à l’instant actuel \(t_{n+1}\) (avec \(\boldsymbol{\epsilon}_{n+1}=\boldsymbol{\epsilon}_n+\Delta \boldsymbol{\epsilon}\)) en chaque point d’intégration, sachant l’état des variables internes et du tenseur des contraintes à l’instant précédent \(t_n\). La prédiction élastique suppose l’incrément du tenseur des déformations totales élastique (bloc \(\texttt{@Predictor}\)). Si le critère de plasticité \(f_{n+1}\) présenté (4247)-2 n’est pas vérifié, la correction plastique doit assurer une solution positionnée sur la surface du domaine de réversibilité (bloc \(\texttt{@Integrator}\)). Dans cette situation, le résidu du système d’équations à résoudre lors de l’intégration est noté \(\boldsymbol{r}\) et le vecteur des inconnues est noté \(\boldsymbol{x}\). Ils sont écrits comme :

(4248)\[\begin{split}\begin{Bmatrix} \boldsymbol{r}_1 \\ r_2 \\ r_3 \end{Bmatrix} =\begin{Bmatrix} -\Delta\boldsymbol{\epsilon} + \Delta\boldsymbol{\epsilon}^e + \Delta\lambda\cfrac{\partial f_{n+1}}{\partial \boldsymbol{\sigma}_{n+1}} \\-\Delta\xi+\Delta\lambda\cfrac{\partial f_{n+1}}{\partial p_{c,n+1}}\\f_{n+1}/K=0\end{Bmatrix},\qquad \begin{Bmatrix} \boldsymbol{x}_1 \\ x_2 \\ x_3 \end{Bmatrix} =\begin{Bmatrix} \Delta\boldsymbol{\epsilon}^e \\ \Delta\xi \\ \Delta\lambda \end{Bmatrix}\end{split}\]

Dans le précédent système, \(\boldsymbol{r}_1=\boldsymbol{0}\) satisfait la décomposition additive du tenseur des déformations totales et la loi d’écoulement normal sur le tenseur des déformations plastiques. La deuxième \(r_2 =0\) traduit le respect de la loi de normalité sur la variable d’écrouissage scalaire. L’équation \(r_3 =0\) exprime que la solution se situe sur la surface du domaine de réversibilité à l’instant \(t_{n+1}\) (l’équation est adimensionnée par le module de compressibilité \(K\)). La méthode itérative de Newton-Raphson est utilisée pour résoudre ce système d’équations \(\boldsymbol{r}=\boldsymbol{0}\), celui-ci étant non-linéaire. Pour cela, on recourt au calcul de la matrice jacobienne \(J\) à chaque itération du schéma, exprimée comme la matrice par bloc suivante :

(4249)\[\begin{split}J = \begin{Bmatrix}\cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\boldsymbol{\epsilon}^e} &\cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\xi} & \cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\lambda} \\\cfrac{\partial r_2}{\partial \Delta\boldsymbol{\epsilon}^e} &\cfrac{\partial r_2}{\partial \Delta\xi} &\cfrac{\partial r_2}{\partial \Delta\lambda} \\\cfrac{\partial r_3}{\partial \Delta\boldsymbol{\epsilon}^e} &\cfrac{\partial r_3}{\partial \Delta\xi} &\cfrac{\partial r_3}{\partial \Delta\lambda}\end{Bmatrix}\end{split}\]

L’obtention de chaque bloc est détaillée en annexe par des calculs analytiques. La valeur du critère d’arrêt du Newton-Raphson est prise à \(10^{-14}\) (cette valeur est renseignée dans la directive \(\texttt{@Epsilon}\)).

Calcul de l’opérateur tangent algorithmique

L’opérateur tangent algorithmique \(\mathbb{D}\), calculé une fois les équations de comportement résolues (bloc \(\texttt{@TangentOperator}\)), s’obtient par une règle de dérivées en chaîne :

\[\mathbb{D} = \frac{\partial \Delta \boldsymbol{\sigma}}{\partial \Delta\boldsymbol{\epsilon}}=\mathbb{C}_{n+1} : \frac{\partial \Delta\boldsymbol{\epsilon}^e}{\partial \Delta\boldsymbol{\epsilon}}\]

Le premier terme \(\mathbb{C}_{n+1}\) désigne le tenseur d’ordre quatre d’élasticité. Celui-ci résulte de l’expression du potentiel d’énergie élastique (4228) comme :

\[\mathbb{C}_{n+1} = \Upsilon(\boldsymbol{\epsilon}^e_{n+1}) \left(3K\mathbb{J}+2\mu\mathbb{K} - 2\kappa\mu\left(\boldsymbol{I}\otimes\boldsymbol{\epsilon}_{d,n+1}^e+\boldsymbol{\epsilon}_{d,n+1}^e\otimes\boldsymbol{I}\right)+\frac{\left(2\kappa\mu\right)^2}{K}\boldsymbol{\epsilon}_{d,n+1}^e\otimes \boldsymbol{\epsilon}_{d,n+1}^e\right)\]

avec \(\mathbb{J}\) et \(\mathbb{K}\) les projecteurs orthogonaux d’ordre quatre, tels que pour tout tenseur d’ordre deux symétrique \(\boldsymbol{a}\) :

\[\mathbb{J}:\boldsymbol{a} = \frac{\mathrm{tr}(\boldsymbol{a})}{3}\boldsymbol{I} \qquad \text{et} \qquad \mathbb{K}:\boldsymbol{a} = \boldsymbol{a}-\frac{\mathrm{tr}(\boldsymbol{a})}{3}\boldsymbol{I}\]

On note que \(\mathbb{C}_{n+1}\) n’est pas constant en raison de la non-linéarité de l’élasticité lorsque \(\kappa>0\).

Le second terme \(\cfrac{\partial \Delta\boldsymbol{\epsilon}^e}{\partial \Delta\boldsymbol{\epsilon}}\) s’écrit quant à lui :

\[\cfrac{\partial \Delta\boldsymbol{\epsilon}^e}{\partial \Delta\boldsymbol{\epsilon}}= J^{-1}_{eel}\]

où \(J^{-1}_{eel}\) représente le premier bloc supérieur gauche de l’inverse de la matrice jacobienne dans (4249). On trouvera la preuve de ce calcul en référence [HMPS15].

Variables internes

Les variables internes, dont les incréments constituent les inconnues du système non-linéaire (4248) à résoudre (directive \(\texttt{@StateVariable}\)), sont résumées dans le tableau suivant.

Tableau 72 Variables internes du modèle.

Appellation	Définition	Symbole	Composantes 2D	Composantes 3D
eel	Tenseur des déformations élastiques	\(\boldsymbol{\epsilon}^e\)	V1-V4	V1-V6
lam	Déformation plastique cumulée	\(\lambda\)	V5	V7
epv	Déformation plastique volumique	\(\xi\)	V6	V8

Quelques variables auxiliaires supplémentaires (directive \(\texttt{@AuxiliaryStateVariable}\)) apparaissent dans le champ des variables internes au sens de Code_Aster. Elles sont résumées dans le tableau ci-dessous.

Tableau 73 Variables auxiliaires du modèle.

Appellation	Définition	Symbole	Composantes 2D	Composantes 3D
OldYieldSize	Rayon à l’instant précédent \(t_n\) du domaine de réversibilité	\(R(\xi_n)\)	V7	V9
YieldSizeGrowth	Accroissement positif du domaine de réversibilité sur l’intervalle \([t_n,t_{n+1}]\)	\(H_n \approx \frac{\langle R(\xi_{n+1})-R(\xi_n)\rangle_+}{\Delta\lambda}\)	V8	V10
PlasticIndicator	0 si élastique, 1 si plastique	\(Ip\)	V9	V11
Dissipation	Densité volumique d’énergie dissipée	\(\overline{D}=\int_0^{t_{n+1}}\left(\boldsymbol{\sigma}:\dot{\boldsymbol{\epsilon}}-\dot{\psi}\right)\mathrm{d}\tau\)	V10	V12

Vérification et validation

Cas-tests

Les références des cas-tests et documentations associées sont données dans le tableau suivant.

Tableau 74 Cas-tests et documentations associées.

Référence cas-test	Référence documentation	Description
ssnv160g	[v6.04.160]	Essai de compression isotrope
ssnv202c	[v6.04.202]	Essai de compression oedométrique
wtnv122d	[v7.31.122]	Essai de compression triaxiale non drainée
comp012h	[v6.07.112]	Test de compatibilité avec la commande CALC_GEO_MECA

Exemples de réponses au point matériel

Les trois figures suivantes présentent quelques réponses obtenues à l’aide de l’exécutable MTest (réponses au point matériel) :

• Essai de compression isotrope ;

• Essai de traction isotrope ;

• Essais de compression triaxiale pour plusieurs contraintres de confinement.

Les paramètres du modèle utilisées dans ces simulations sont regroupés dans le tableau ci-dessous.

Tableau 75 Paramètres du modèle simulé.

Effet	Appellation	Définition	Symbole	Valeur
Élasticité non-linéaire isotrope	BulkModulus ShearModulus SwellingIndex	Module de compressibilité Module de cisaillement Indice de non-linéarité élastique	\(K\) \(\mu\) \(\kappa\)	160 MPa 100 MPa 50
Domaine d’élasticité initial	InitCritPress CritStateSlope TensileYieldStress	Pression critique initiale Pente d’état critique Limite d’élasticité en traction isotrope	\(p_{c0}\) \(M\) \(\sigma_0\)	1 MPa 1 10 kPa
Écrouissage cinématique-isotrope	IncoPlastIndex	Indice d’incompressibilité plastique	\(\beta\)	50

Fig. 297 Essai de compression isotrope piloté par la contrainte moyenne \(\sigma_m\). Dans la phase de charge, la réponse est élastique non-linéaire (\(Ip=0\)) jusque \(-\sigma_m=2p_{c0}+\sigma_0=2.01\) MPa, puis irréversible (\(Ip=1\)). L’écrouissage est alors positif. La décharge depuis \(-\sigma_m=10\) MPa jusque zéro est élastique non-linéaire (\(Ip=0\)).

Fig. 298 Essai de traction isotrope piloté par la déformation verticale \(\epsilon_{zz}\). Après la phase de charge élastique, la contrainte sature à \(\sigma_m=\sigma_0=10\) kPa, sans montrer d’écrouissage ni positif ni négatif. On note que la non-linéarité de l’élasticité est imperceptible. L’expression de la contrainte (4229) l’explique par la condition \(\sigma_m\ll K/\kappa\), situation rencontrée dans ce chargement avec les paramètres du modèle simulé.

Fig. 299 Essais de compression triaxiale à plusieurs contraintes de confinement \(-\sigma_{xx}=-\sigma_{yy}\). La contrainte équivalente \(\sigma_{eq}\) présente un écrouissage négatif aux plus faibles confinements et un écrouissage positif aux plus forts. En parallèle, l’évolution de la déformation plastique volumique \(\xi=\mathrm{tr}(\boldsymbol{\epsilon}^p)\) passe d’un régime de dilatance (\(\dot{\xi}>0\)) à contractance (\(\dot{\xi}<0\)). On notera que la première phase de chargement à contrainte équivalente \(\sigma_{eq}=0\) correspond à la mise sous confinement avec \(-\sigma_{xx}=-\sigma_{yy}=-\sigma_{zz}\).

Bibliographie

[RoBu68]

Roscoe, K.H. and Burland, J.B. On the Generalized Stress-Strain Behavior of Wet Clay. Engineering Plasticity, Cambridge University Press, Cambridge, 1968, pp.535-609.

[HaNg75]

Halphen, B. and Nguyen, Q.S. Sur les matériaux standard généralisés. Journal de Mécanique, 1975, Volume 14, pp.39-63.

[Mial86]

Mialon, P. Éléments d’analyse et de résolution numérique des relations de l’élasto-plasticité. EDF, Bulletin de la Direction des Etudes et Recherches, 1986, Série C, Numéro 3, pp. 57-89.

[BoLe90]

Borja, R.I. and Lee, S.R. Cam-Clay plasticity, Part 1: Implicit integration of elasto-plastic constitutive relations. Computer Methods in Applied Mechanics and Engineering, 1990, Volume 78, Issue 1, pp.49-72.

[Bour97]

Bourgeois, E. Mécanique des milieux poreux en transformation finie : Formulation des problèmes et méthodes de résolution. Thèse de doctorat, ENPC, 1997.

[OrSt99]

Ortiz, M. and Stainier, L. The variational formulation of viscoplastic constitutive updates. Computer Methods in Applied Mechanics and Engineering, 1999, Volume 171, Issues 3-4, pp.419-444.

[HMPS15] (1,2)

Helfer, T., Michel, B., Proix, J.M., Salvo, M., Sercombe, J. and Casella, M. Introducing the open-source mfront code generator: Application to mechanical behaviours and material knowledge management within the PLEIADES fuel element modelling platform. Computers & Mathematics with Applications, 2015, Volume 70, Issue 5, pp.994-1023.

[Bacq23] (1,2)

Bacquaert, G. Comportement des géomatériaux pour la modélisation et la sûreté des ouvrages géotechniques. Thèse de doctorat, Sorbonne Université, 2023.

[BRAV24]

Bacquaert, G., Raude S., Aves-Fernandes, V., Voldoire F. and Kondo, D. A standard thermodynamic-based extension of the Modified Cam-Clay soil model and its applications. European Journal of Mechanics - A/Solids, 2024, Volume 103, 105122.

Annexe : expression de la matrice jacobienne

On rappelle l’expression du systèmes d’équations non-linéaires (4248) et la définition de sa matrice jacobienne (4249) :

(4250)\[\begin{split}\boldsymbol{r}= \begin{Bmatrix} \boldsymbol{r}_1 \\ r_2 \\ r_3 \end{Bmatrix} =\begin{Bmatrix} -\Delta\boldsymbol{\epsilon} + \Delta\boldsymbol{\epsilon}^e + \Delta\lambda\cfrac{\partial f_{n+1}}{\partial \boldsymbol{\sigma}_{n+1}} \\- \Delta\xi+\Delta\lambda\cfrac{\partial f_{n+1}}{\partial p_{c,n+1}}\\f_{n+1}/K=0\end{Bmatrix},\qquad J = \begin{Bmatrix} \cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\boldsymbol{\epsilon}^e} &\cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\xi} & \cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\lambda} \\\cfrac{\partial r_2}{\partial \Delta\boldsymbol{\epsilon}^e} &\cfrac{\partial r_2}{\partial \Delta\xi} &\cfrac{\partial r_2}{\partial \Delta\lambda} \\\cfrac{\partial r_3}{\partial \Delta\boldsymbol{\epsilon}^e} &\cfrac{\partial r_3}{\partial \Delta\xi} &\cfrac{\partial r_3}{\partial \Delta\lambda}\end{Bmatrix}\end{split}\]

Dans MFront, cette matrice \(J\) peut être obtenue par perturbation numérique ou analytiquement, comme c’est le cas présenté par la suite. Ses composantes sont détaillées dans le cas d’un incrément avec correction plastique. Par commodité, on notera :

• Les directions d’écoulement de \(\Delta\boldsymbol{\epsilon}^p\) et \(\Delta\xi\) :

(4251)\[\begin{split}\begin{align} &\boldsymbol{N}_{n+1} = \frac{\partial f_{n+1}}{\partial \boldsymbol{\sigma}_{n+1}} = \frac{\cfrac{3}{2M^2}\boldsymbol{A}_{d,n+1}+A_{m,n+1}\cfrac{\boldsymbol{I}}{3}}{T_{eq}(\boldsymbol{A}_{n+1})} \\ &N_{v,n+1} = \frac{\partial f_{n+1}}{\partial p_{c,n+1}} = \mathrm{tr}\left(\boldsymbol{N}_{n+1}\right) \qquad \text{avec} \qquad \boldsymbol{A}=\boldsymbol{\sigma}+p_c\boldsymbol{I},\quad T_{eq}(\boldsymbol{A}) = \sqrt{\left(\frac{A_{eq}}{M}\right)^2+A_m^2} \end{align}\end{split}\]

• Les deux dérivées suivantes :

(4252)\[\begin{split}\begin{align} &\boldsymbol{N}_{n+1,\boldsymbol{A}_{n+1}} = \frac{\partial \boldsymbol{N}_{n+1}}{\partial \boldsymbol{A}_{n+1}} = \frac{\cfrac{3}{2M^2}\mathbb{K}+\cfrac{\mathbb{J}}{3}-\boldsymbol{N}_{n+1}\otimes\boldsymbol{N}_{n+1}}{T_{eq}(\boldsymbol{A}_{n+1})}\\ &N_{v,n+1,\boldsymbol{A}_{n+1}} = \frac{\partial N_{v,n+1}}{\partial \boldsymbol{A}_{n+1}} = \frac{\cfrac{\boldsymbol{I}}{3}-N_{v,n+1}\boldsymbol{N}_{n+1}}{T_{eq}(\boldsymbol{A}_{n+1})} \end{align}\end{split}\]

Première ligne

La dérivation de chaque terme de la première ligne du système présenté (4250) fournit :

(4253)\[\begin{split}\begin{align} &\frac{\partial \boldsymbol{r}_1}{\partial \Delta\boldsymbol{\epsilon}^e} = \mathbb{I} + \Delta\lambda\boldsymbol{N}_{n+1,\boldsymbol{A}_{n+1}}:\mathbb{C}_{n+1}\\ &\frac{\partial \boldsymbol{r}_1}{\partial \Delta\lambda}=\boldsymbol{N}_{n+1}\\ &\frac{\partial \boldsymbol{r}_1}{\partial \Delta\xi}=\Delta\lambda\boldsymbol{N}_{n+1,\boldsymbol{A}_{n+1}}:\frac{\partial \boldsymbol{A}_{n+1}}{\partial \Delta\xi} \end{align}\end{split}\]

où \(\frac{\partial \boldsymbol{A}_{n+1}}{\partial \Delta\xi} = p_{c,n+1}'\boldsymbol{I} = -\beta p_{c,n+1}\boldsymbol{I}\) conformément à l’expression de la pression critique définie (4231) et \(\mathbb{I}\) désigne le tenseur identité d’ordre quatre.

Deuxième ligne

La deuxième ligne du système s’écrit :

(4254)\[\begin{split}\begin{align} &\frac{\partial r_2}{\partial \Delta\boldsymbol{\epsilon}^e} = \Delta\lambda N_{v,n+1,\boldsymbol{A}_{n+1}}:\mathbb{C}_{n+1}\\ &\frac{\partial r_2}{\partial \Delta\lambda}=N_{v,n+1}\\ &\frac{\partial r_2}{\partial \Delta\xi}=-1 + \Delta\lambda N_{v,n+1,\boldsymbol{A}_{n+1}}:\frac{\partial \boldsymbol{A}_{n+1}}{\partial \Delta\xi} \end{align}\end{split}\]

Troisième ligne

Enfin, la troisième ligne du système est :

(4255)\[\begin{split}\begin{align} &\frac{\partial r_3}{\partial \Delta\boldsymbol{\epsilon}^e} = \boldsymbol{N}_{n+1}:\mathbb{C}_{n+1}/K\\ &\frac{\partial r_3}{\partial \Delta\lambda}=-H_n/K\\ &\frac{\partial r_3}{\partial \Delta\xi}=\boldsymbol{N}_{n+1}:\frac{\partial \boldsymbol{A}_{n+1}}{\partial \Delta\xi}/K \end{align}\end{split}\]