r7.01.44 Modèle de comportement CSSM

Résumé :

Ce document présente la formulation et l’intégration numérique du modèle de comportement élastoplastique CSSM (« Critical State Soil Model »). Le modèle CSSM est une combinaison de deux modèles visant à représenter respectivement le comportement monotone et cyclique drainé des sols. Le premier composant s’appuie sur le modèle de Cam-Clay modifié [RoBu68], qui décrit les phénomènes de dilatance ou contractance ainsi que des effets de confinement sous des chargements monotones et drainés. Le second composant est un modèle multimécanisme de type Iwan [Iwan66], conçu pour représenter le comportement non-linéaire cyclique des sols à faibles niveaux de déformation.

Le modèle de comportement CSSM comporte un total de onze paramètres : deux pour décrire le comportement réversible (élastique), cinq pour le comportement irréversible lié au composant Cam-Clay modifié, et trois pour le comportement irréversible lié au composant Iwan. Un dernier paramètre pondère les contributions des deux composants sur le comportement au cisaillement du modèle CSSM.

Les équations de comportement sont résolues en utilisant une méthode d’intégration implicite directe en temps, implémentée dans l’outil MFront [HMPS15].

Formulation du modèle

Le modèle CSSM s’appuie sur l’association des équations de comportement similaires aux modèles de Cam-Clay modifié et d’Iwan. Ceux-ci sont respectivement détaillés dans les documentations [r7.01.38] et [r7.01.48].

Une description rhéologique du modèle CSSM consiste à superposer ces deux précédents modèles comme deux composants, le choix se portant sur :

• Une association en parallèle du comportement déviatorique.

• Une association en série du comportement volumique.

Ces associations sont très schématiquement représentées sur la Fig. 266.

Fig. 266 Représentation rhéologique des composants mis en parallèle (comportement déviatorique) et en série (comportement volumique). Les composants 1 et 2 voient la même contrainte moyenne \(\sigma_m\) mais pas le même déviateur des contraintes \(\boldsymbol{\sigma}_d\).

Énergie et lois d’état

Les variables d’état du modèle CSSM sont le tenseur des déformations totales \(\boldsymbol{\varepsilon}\) et une collection de variables internes propres à ses deux composants :

• Pour le premier composant : trois variables internes, respectivement tensorielle et scalaires, notées \(\boldsymbol{\varepsilon}^p,\xi\) et \(\gamma\).

• Pour le second composant : \(N\) variables internes tensorielles, notées \(\left(\boldsymbol{\alpha}^i\right)_{1\leq i\leq N}\).

Le potentiel d’énergie libre \(\psi\) est défini comme suit :

(4103)\[\psi(\boldsymbol{\varepsilon},\boldsymbol{\varepsilon}^p,\boldsymbol{\alpha}^1,\dots,\boldsymbol{\alpha}^N, \xi,\gamma) = \psi_{e,v}\left(\varepsilon_v-\varepsilon_v^p-\sum_{i=1}^N\alpha^i_v\right)+\left[\rho\psi_{e,d}(\boldsymbol{\varepsilon}_d-\boldsymbol{\varepsilon}_d^p) +(1-\rho)\psi_{e,d}\left(\boldsymbol{\varepsilon}_d-\sum_{i=1}^N\boldsymbol{\alpha}^i_d\right)\right] +\psi_{h,1}(\xi,\gamma) + \psi_{h,2}(\boldsymbol{\alpha}^1,\dots,\boldsymbol{\alpha}^N)\]

Dans ce potentiel, on note que :

• Le premier terme \(\psi_{e,v}\) tient lieu d’énergie élastique volumique.

• Le second terme, entre crochets, représente l’énergie élastique de cisaillement, celle-ci étant pondérée entre les composants 1 et 2, \(\rho\psi_{e,d}\) et \((1-\rho)\psi_{e,d}\), au prorata d’un coefficient \(\rho\).

• Les deux derniers termes \(\psi_{h,1}\) et \(\psi_{h,2}\) sont les énergies stockées par écrouissage.

Les sections suivantes précisent les expressions des fonctions \(\psi_{e,v},\psi_{e,d},\psi_{h,1}\) et \(\psi_{h,2}\).

Énergie élastique

Les énergies élastique volumique et de cisaillement sont quadratiques :

(4104)\[\psi_{e,v}(x)= \frac{K}{2}x^2,\quad \psi_{e,d}(\boldsymbol{x})= \mu\boldsymbol{x}:\boldsymbol{x}\]

\(K\) est le module de compressibilité et \(\mu\) est le module de cisaillement.

Énergie stockée par écrouissage du premier composant

L’énergie stockée par écrouissage \(\psi_{h,1}\), pilotée par les variables d’écrouissage \(\xi\) et \(\gamma\), est définie comme :

(4105)\[\psi_{h,1}(\xi,\gamma) = p_{c0}\left(\frac{\exp\left(-\beta\xi\right)-1}{\beta}+\eta\frac{\exp\left(-\omega\gamma\right)-1}{\omega}\right)\]

Les paramètres \(\beta\) (indice d’incompressibilité plastique) et \(p_{c0}\) (pression critique initiale) sont ceux invoqués dans l’expression de la pression critique du modèle Cam-Clay modifié [r7.01.48].

Les indices \(\eta\) et \(\omega\), qui portent sur l’écrouissage piloté par la variable \(\gamma\), ont notamment pour effet de moduler les réponses du modèle lors de chargements de compression isotrope. Ceci est illustré dans la r7.01.44-cas_tests.

Énergie stockée par écrouissage du second composant

L’énergie stockée par écrouissage \(\psi_{h,2}\), pilotée par les variables d’écrouissage \(\left(\boldsymbol{\alpha}^i\right)_{1\leq i\leq N}\), est définie comme :

(4106)\[\psi_{h,2}(\boldsymbol{\alpha}^1,\dots,\boldsymbol{\alpha}^N) = \sum_{i=1}^N \left(\frac{H^i_d}{2}\boldsymbol{\alpha}^i_d:\boldsymbol{\alpha}^i_d+\frac{H^i_v}{2}(\alpha^i_v)^2\right)\]

où \(H_d^i\) et \(H_v^i\) sont respectivement les modules d’écrouissage cinématique déviatorique et volumique associés à \(\boldsymbol{\alpha}^i\).

Lois d’état

Ayant établi l’expression du potentiel d’énergie via ses composantes, on peut écrire la densité volumique de dissipation intrinsèque \(D\) comme :

(4107)\[D=\boldsymbol{\sigma}:\dot{\boldsymbol{\varepsilon}}-\dot{\psi}=\boldsymbol{X}:\dot{\boldsymbol{\varepsilon}}^p+ \boldsymbol{A}^i:\dot{\boldsymbol{\alpha}}^i+p_c\dot{\xi}+S\dot{\gamma}\]

en définissant le tenseur des contraintes \(\boldsymbol{\sigma}\) et les forces irréversibles \(\boldsymbol{X},\boldsymbol{A}^i,p_c\) et \(S\) conjuguées aux variables internes par les lois d’état :

(4108)\[\boldsymbol{\sigma} =\frac{\partial \psi}{\partial \boldsymbol{\varepsilon}},\quad\boldsymbol{X} =-\frac{\partial \psi}{\partial \boldsymbol{\varepsilon}^p},\quad \boldsymbol{A}^i =-\frac{\partial \psi}{\partial \boldsymbol{\alpha}^i},\quad p_c = -\frac{\partial \psi}{\partial \xi},\quad S =-\frac{\partial \psi}{\partial \gamma}\]

Ces lois peuvent être précisées à partir de (4225), (4228), (4105) et (4106) :

(4109)\[\begin{split}\begin{align} \boldsymbol{\sigma} &= K\left(\varepsilon_v-\varepsilon_v^p-\sum_{j=1}^N\alpha^j_v\right)\boldsymbol{I} + 2\mu\left(\rho\left(\boldsymbol{\varepsilon}_d-\boldsymbol{\varepsilon}_d^p\right)+(1-\rho)\left(\boldsymbol{\varepsilon}_d-\sum_{j=1}^N\boldsymbol{\alpha}^j_d\right)\right)\vphantom{\Big{(}}\\ \boldsymbol{X} &= K\left(\varepsilon_v-\varepsilon_v^p-\sum_{j=1}^N\alpha^j_v\right)\boldsymbol{I} + 2\mu\rho\left(\boldsymbol{\varepsilon}_d-\boldsymbol{\varepsilon}_d^p\right)\vphantom{\Big{(}}\\ \boldsymbol{A}^i &= K\left(\varepsilon_v-\varepsilon_v^p-\sum_{j=1}^N\alpha^j_v\right)\boldsymbol{I} + 2\mu(1-\rho)\left(\boldsymbol{\varepsilon}_d-\sum_{j=1}^N\boldsymbol{\alpha}^j_d\right)-H_d^i\boldsymbol{\alpha}^i_d-H_v^i\alpha^i_v\boldsymbol{I}\vphantom{\Big{(}}\\ p_c &= p_{c0}\exp\left(-\beta\xi\right)\vphantom{\Big{(}}\\ S &= \eta p_{c0}\exp\left(-\omega\gamma\right)\vphantom{\Big{(}} \end{align}\end{split}\]

Remarque :

À des fins de commodité pour la résolution des équations de comportement avec MFront dans la r7.01.44-integration_numerique, nous écrirons avantageusement les trois premières lois d’état (4109)-1, (4109)-2 et (4109)-3 sous la forme :

(4110)\[\boldsymbol{\sigma} = \boldsymbol{\sigma}(\boldsymbol{\varepsilon}^e)=K\varepsilon_v^e\boldsymbol{I} + 2\mu\boldsymbol{\varepsilon}_d^e,\quad \boldsymbol{X} = \boldsymbol{X}(\boldsymbol{\varepsilon}^x) = K\varepsilon_v^x\boldsymbol{I} + 2\mu\rho\boldsymbol{\varepsilon}_d^x,\quad \boldsymbol{A}^i = \boldsymbol{\sigma}-\boldsymbol{X}_d-H_d^i\boldsymbol{\alpha}^i_d-H_v^i\alpha^i_v\boldsymbol{I}\]

où sont introduites deux déformations élastiques, \(\boldsymbol{\varepsilon}^e\) et \(\boldsymbol{\varepsilon}^x\), respectivement en bijection à \(\boldsymbol{\sigma}\) et \(\boldsymbol{X}\). Celles-ci sont facilement exprimables en fonction du tenseur des déformations totales et des variables internes par les relations :

(4111)\[\boldsymbol{\varepsilon}^e = \boldsymbol{\varepsilon} - \left(\mathbb{J}+\rho\mathbb{K}\right):\boldsymbol{\varepsilon}^p -\left(\mathbb{J}+(1-\rho)\mathbb{K}\right):\sum_{i=1}^N\boldsymbol{\alpha}^i,\quad \boldsymbol{\varepsilon}^x = \boldsymbol{\varepsilon}^e - \left(1-\rho\right)\mathbb{K}:\left(\boldsymbol{\varepsilon}^p-\sum_{i=1}^N\boldsymbol{\alpha}^i\right)\]

Lois d’évolution

Les lois d’évolution sont déduites des critères de plasticité propres aux deux composants du modèle CSSM ci-dessous définis.

Critère de plasticité et règle d’écoulement du premier composant

L’expression du critère de plasticité est l’équation d’une ellipse dans le plan \((X_m,X_{eq})\) :

(4112)\[f(\boldsymbol{X},p_c,S;\xi) = \sqrt{\left(\frac{X_{eq}}{M}\right)^2 + \left(X_m+p_c-S\right)^2}+S-R(\xi),\quad \text{avec} \quad R(\xi) = p_{c0}\exp(-\beta\xi)\]

Le paramètre \(M\), appelé pente d’état critique, porte la même interprétation que dans le modèle Cam-Clay modifié [r7.01.48]. Géométriquement, \(p_c-S\) positionne le centre de l’ellipse suivant l’axe hydrostatique, et \(S-R(\xi)\) précise sa taille.

Remarque :

Suivant l’expression du critère de plasticité (4112), il faut 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. Par ailleurs, le critère de plasticité 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 premier composant du modèle CSSM présente, tout comme pour le modèle Cam-Clay modifié [r7.01.48], un écrouissage combiné cinématique-isotrope piloté par \(\xi\). À cet écrouissage s’ajoute celui piloté par la variable d’écrouissage \(\gamma\) du fait de l’intervention de sa force conjuguée \(S\) dans le critère de plasticité (4112).

Le domaine de réversibilité initial associé au critère de plasticité (4112) s’explicite à l’aide des lois d’état (4109) par :

(4113)\[f_{0} = \sqrt{\left(\frac{\rho\sigma_{eq}}{M}\right)^2 + \left(\sigma_m+p_{c0}(1-\eta)\right)^2}-p_{c0}(1-\eta)\]

Lorsque \(\rho=0\), le critère de plasticité initial ne dépend donc plus de la contrainte équivalente de von Mises \(\sigma_{eq}\). Néanmoins, quelque soit la valeur de \(\rho\), la limite en traction isotrope est zéro, et vaut \(2p_{c0}(1-\eta)\) en compression isotrope. La Fig. 267 représente l’effet de \(\eta\) sur le domaine de réversibilité initial. Plus celui-ci s’approche de un, plus la taille du domaine tend vers zéro par réduction homothétique autour de l’origine.

Fig. 267 Domaines de réversibilité définis par le critère de plasticité (4113) pour trois valeurs de \(\eta\).

L’évolution des trois variables internes \(\boldsymbol{\varepsilon}^p,\xi\) et \(\gamma\) respecte la règle de normalité au critère de plasticité (4112) :

(4114)\[\begin{split}\begin{align} \dot{\boldsymbol{\varepsilon}}^p &= \dot{\lambda}\frac{\partial f}{\partial \boldsymbol{X}} = \dot{\lambda}\frac{\cfrac{3}{2M^2}\boldsymbol{X}_d + \left(X_m+p_c-S\right)\cfrac{\boldsymbol{I}}{3}}{\sqrt{\left(\cfrac{X_{eq}}{M}\right)^2 + \left(X_m+p_c-S\right)^2}}\\ \dot{\xi} &= \dot{\lambda}\frac{\partial f}{\partial p_c} = \dot{\lambda}\frac{X_m+p_c-S}{\sqrt{\left(\cfrac{X_{eq}}{M}\right)^2 + \left(X_m+p_c-S\right)^2}} \\ \dot{\gamma} &= \dot{\lambda}\frac{\partial f}{\partial S} = \dot{\lambda}-\dot{\lambda}\frac{X_m+p_c-S}{\sqrt{\left(\cfrac{X_{eq}}{M}\right)^2 + \left(X_m+p_c-S\right)^2}} \end{align}\end{split}\]

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

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

On obtient donc à la suite de (4114) les relations :

(4116)\[\dot{\xi}=\dot{\varepsilon}_v^p, \quad \dot{\gamma}=\dot{\lambda}-\dot{\varepsilon}_v^p, \quad \dot{\lambda}=\sqrt{(\dot{\varepsilon}_v^p)^2+(M\dot{\varepsilon}_{eq}^p)^2}\]

La première variable d’écrouissage \(\xi\) est la déformation plastique volumique. La seconde \(\gamma\) est la différence entre la déformation plastique cumulée \(\lambda\) et celle volumique. Étant donné que (4116) montre que \(\dot{\gamma}\geq 0\), l’écrouissage porté par sa force conguguée \(S\) est monotone décroissant et évanescent d’après (4109)-5 (\(\dot{S}\leq 0\) et \(\lim_{\gamma\rightarrow\infty} S = 0\)).

Critère de plasticité et règle d’écoulement du second composant

La loi d’évolution des variables internes \(\left(\boldsymbol{\alpha}^i\right)_{1\leq i\leq N}\) provient de critères de plasticité se différenciant uniquement par leur taille. Ces critères sont hémi-elliptiques dans les plans \((A^i_m,A^i_{eq})\) :

(4117)\[F_i(\boldsymbol{A}^i) = \sqrt{\left(A^i_{eq}\right)^2 + \left(\frac{R_i}{C}\left\langle A^i_m+C\right\rangle\right)^2} -R_i,\quad \text{avec}\quad \langle x \rangle = \frac{x+|x|}{2}\]

La constante \(C\) positionne la transition entre les ellipses et les droites des critères de plasticité, cette valeur étant supposée commune aux \(N\) critères. Les constantes \(R_i\) sont les seuils maximaux atteignables par les normes équivalentes de von Mises \(A^i_{eq}\). La Fig. 268 illustre trois des critères (4117) pour des valeurs croissantes \(R_1,R_2\) et \(R_3\).

Fig. 268 Domaines de réversibilité définis par les critères de plasticité (4117) pour trois valeurs croissantes de \(R_1,R_2\) et \(R_3\).

Remarque :

Lorsque \(C\rightarrow 0\), les critères de plasticité (4117) s’approchent de ceux utilisés dans le modèle d’Iwan [r7.01.38] (critères de von Mises), avec dans le cas présent une coupure interdisant le demi-espace de traction (\(A^i_m>0\)).

Le domaine de réversibilité initial associé aux critères de plasticité (4117) s’explicite à l’aide des lois d’état (4109) par :

(4118)\[F_{1,0} = \sqrt{\left((1-\rho)\sigma_{eq}\right)^2 + \left(\frac{R_1}{C}\left\langle \sigma_m+C\right\rangle\right)^2} -R_1\]

ayant supposé \(R_1\) la plus petite valeur des \(R_i\). Lorsque \(\rho=1\), le critère de plasticité initial ne dépend donc plus de la contrainte équivalente de von Mises \(\sigma_{eq}\).

En combinant (4113) et (4118), le domaine d’élasticité initial du modèle CSSM est donc défini par l’intersection \(\left\lbrace f_{0}\leq 0 \cap F_{1,0}\leq 0\right\rbrace\).

L’évolution des variables internes \(\left(\boldsymbol{\alpha}^i\right)_{1\leq i\leq N}\) respecte la règle de normalité aux critères de plasticité (4117) :

(4119)\[\begin{split}\dot{\boldsymbol{\alpha}}^i &= \dot{\lambda}_i\frac{\partial F_i}{\partial \boldsymbol{A}^i} = \dot{\lambda}_i\frac{\cfrac{3}{2}\boldsymbol{A}^i_d + \left(\cfrac{R_i}{C}\right)^2\left\langle A^i_m+C\right\rangle\cfrac{\boldsymbol{I}}{3}}{\sqrt{\left(A^i_{eq}\right)^2 + \left(\cfrac{R_i}{C}\left\langle A^i_m+C\right\rangle\right)^2}}\\\end{split}\]

où \(\dot{\lambda}_i\) sont les multiplicateurs plastiques satisfaisant les conditions de cohérence suivantes :

(4120)\[\dot{\lambda}_i\geq 0, \quad F_i\leq 0, \quad \dot{\lambda}_iF_i=0\]

Détermination des modules d’écrouissage cinématique et des seuils associés

Dans cette partie, on explicite comment la détermination des modules modules d’écrouissage cinématique \(H_d^i,H_v^i\) et des seuils \(R_i\) peut être simplifiée en s’appuyant sur deux paramètres alternatifs : \(\gamma_{\mathrm{hyp}}\), \(n_\mathrm{hyp}\).

Calcul des modules d’écrouissage déviatorique et des seuils

Le calcul des modules d’écrouissage déviatorique et des seuils suit la même méthodologie que celle présentée dans la documentation [r7.01.38] en supposant vérifiées les conditions suivantes :

• La surface de plasticité du premier composant n’est pas atteinte (\(f<0\)).

• Les déformations volumiques de chaque \(i<N\)-ième mécanisme du second composant sont constantes (\(\dot{\alpha}_v^i=0\)).

En quelques mots, la seconde condition est toujours vérifiée si le \(i<N\)-ième mécanisme plastifie suffisamment de sorte que \(\alpha_v^i\) puisse se stabiliser par un écrouissage cinématique suivant l’axe hydrostatique, où bien plus simplement lorsque \(-\sigma_m\geq C\).

Le calcul des modules d’écrouissage déviatorique et des seuils se ramène alors à s’intéresser à un chemin monotone en chargement de cisaillement à pression constante (\(-\dot{\sigma}_m=0\)), par exemple dans le plan \((x,y)\). Dans ce chargement, la relation entre les taux de déformation \(\dot{\varepsilon}_{xy}\) et de contrainte \(\dot{\sigma}_{xy}\) s’écrit :

(4121)\[\dot{\varepsilon}_{xy} = \left(2\rho\mu + \frac{1}{\cfrac{1}{2(1-\rho)\mu}+\sum_{i=1}^N\cfrac{\mathbf{H}(\sqrt{3}|\sigma_{xy}|-R_i)}{H_d^i}}\right)^{-1}\dot{\sigma}_{xy}\]

avec \(\mathbf{H}\) la fonction échelon unité. En supposant la réponse connue en \(N\) points \((\varepsilon^{(i)}_{xy},\sigma^{(i)}_{xy})_{1\leq i\leq N}\), on peut en déduire :

(4122)\[\begin{split}R_i = \sqrt{3}|\sigma^{(i)}_{xy}|,\quad H_i^d = \begin{cases} \left\langle\left(\cfrac{\sigma^{(i+1)}_{xy}-\sigma^{(i)}_{xy}}{\varepsilon^{(i+1)}_{xy}-\varepsilon^{(i)}_{xy}}-2\rho\mu\right)^{-1} - \cfrac{1}{2(1-\rho)\mu} - \sum_{j=1}^{i-1} \cfrac{1}{H_d^j}\right\rangle^{-1} &\quad \text{si} \quad i<N\\ 0 &\quad \text{si} \quad i=N \end{cases}\end{split}\]

Notons alors que (4122)-2 entraîne \(H_d^{m}=0\Longrightarrow\left(H_d^i\right)_{m+1\leq i\leq N}=0\). \(m-1<N\) modules sont donc nuls avec systématiquement \(H_d^N=0\), ce qui permet au modèle de prédire un état critique conformément au système que nous établissons (4241).

Le calcul des paramètres d’écrouissage \(H_d^i\) est réalisé en fixant \(N=10\) points d’interpolation sur une échelle logarithmique de base dix, entre \(2\varepsilon^{(1)}_{xy}=10^{-5}\) et \(2\varepsilon^{(10)}_{xy}=10^{-2}\). Les points \((\varepsilon^{(i)}_{xy},\sigma^{(i)}_{xy})_{1\leq i\leq 10}\) satisfont une relation « hyperbolique modifiée », couramment utilisée en dynamique des sols, qui s’écrit :

(4123)\[\mu_{\mathrm{secant}} = \frac{\mu}{1+\left(\cfrac{2\varepsilon_{xy}}{\gamma_{\mathrm{hyp}}}\right)^{n_{\mathrm{hyp}}}},\quad\text{avec}\quad \mu_{\mathrm{secant}}=\frac{\sigma_{xy}}{2\varepsilon_{xy}}\]

\(\gamma_{\mathrm{hyp}}\) et \(n_\mathrm{hyp}\) sont donc deux paramètres du modèle CSSM à renseigner, lesquels ont une interprétation commune aux paramètres d’écrouissage documentés dans [r7.01.38].

Expression des modules d’écrouissage volumique

L’obtention des modules d’écrouissage volumique \(H_v^i\) est justifiée par la volonté d’intégrer efficacement les équations d’évolution présentes dans (4117) et (4120), sans nécessiter la résolution d’une équation non-linéaire pour les multiplicateurs plastiques. En effet, on prouve que cela est possible lorsque la condition suivante est satisfaite :

(4124)\[H_v^i = \frac{3}{2}\left(\frac{C}{R_i}\right)^2H_d^i\]

Sous cette condition, l’intégration numérique du modèle CSSM permet alors de réduire de le nombre d’équations à résoudre. Cette procédure est détaillée dans la r7.01.44-integration_numerique.

Équations de l’état critique

Dans ce paragraphe, nous cherchons à établir les équations du modèle CSSM correspondant à l’état critique. Cet état est défini par des contraintes et une déformation volumique constantes, asympotiquement atteint dans un chargement monotone. Nous supposons que cet état est atteint pour chacun des composants 1 et 2 du modèle.

Pour le premier composant, si les contraintes restent constantes, il n’y a plus d’écrouissage dû à \(\xi\) ni à \(\gamma\). Dans cette situation, on prouve que les conditions suivantes sont vérifiées :

(4125)\[X_{eq}+M X_m=0,\quad X_m+p_c=0\]

Or, les lois d’état (4109)-1 et (4109)-2 montrent que \(X_m=\sigma_m\). (4125) est donc équivalent à :

(4126)\[X_{eq}+M \sigma_m=0,\quad \sigma_m+p_c=0\]

En ce qui concerne le second composant, suite à (4122)-2 et (4124), \(H_d^{m}=H_v^{m}=0=0\) pour un certain \(m\leq N\). Les lois d’état (4109)-1 et (4109)-3 montrent donc que \(A^m_m=\sigma_m\). La vérification du critère (4117) conduit alors à :

(4127)\[A^m_{eq} = R_m\sqrt{1-\left\langle \sigma_m/C+1\right\rangle^2}\]

De plus, on rappelle (4109)-1, (4109)-2 et (4109)-3 qui établissent que \(\boldsymbol{\sigma}_d=\boldsymbol{X}_d+\boldsymbol{A}^m_d\). En admettant que ces deux derniers tenseurs sont colinéaires et de même sens, les équations (4126)-1 et (4127) impliquent que :

(4128)\[\sigma_{eq}+M\sigma_m=R_m\sqrt{1-\left\langle \sigma_m/C+1\right\rangle^2}\]

Enfin, on montre que la déformation volumique totale reste constante seulement si \(\sigma_m+C<0\). Dans le cas contraire, celle-ci évolue uniquement comme \(\dot{\varepsilon}_v=\dot{\alpha}_v^m>0\).

En bref, on retiendra que l’état critique prédit par le modèle CSSM peut être défini par le jeu d’équations suivant :

(4129)\[\begin{split}\begin{align} \dot{\boldsymbol{\sigma}}=\boldsymbol{0}\\ \dot{\varepsilon}_v=0\\ \sigma_m+p_c = 0\\ \sigma_m+C<0 \\ \sigma_{eq}+M\sigma_m=R_m\sqrt{1-\left\langle \sigma_m/C+1\right\rangle^2} \end{align}\end{split}\]

La Fig. 269 représente trois lignes d’état critique (4241)-5 selon des valeurs de \(R_m\) dans le plan méridien des contraintes. Pour \(\sigma_m+C<0\), on y indique que le modèle CSSM prédit un état de dilatance (\(\dot{\varepsilon}_v=\dot{\alpha}_v^m>0\)).

Fig. 269 Exemples de lignes d’état critique (4241)-5 pour trois valeurs de \(R_m\) dans le plan méridien des contraintes.

Suppression de l’écrouissage négatif du modèle

Le modèle CSSM hérite du même écrouissage que celui du modèle Cam-Clay modifié, piloté par la variable \(\xi\). Or nous avons montré à la suite du critère de plasticité (4112) que cette variable produit deux écrouissages :

• L’un cinématique : ce mécanisme est associé à sa force conjuguée congugée \(p_c\) par la loi d’état (4109)-4, jouant le rôle de force de rappel hydrostatique. Cet écrouissage déplace la surface de plasticité dans le plan méridien des contraintes sans modifier sa taille ou sa forme.

• L’autre isotrope : ce mécanisme est lié à la paramétrisation du critère à travers la fonction \(R(\xi)\). Cet écrouissage modifie la taille de la surface de plasticité.

L’écrouissage cinématique seul ne peut pas produire un écrouissage négatif, en raison de la convexité du potentiel à partir duquel \(p_c\) dérive. En revanche, un écrouissage isotrope négatif peut survenir si \(R(\xi)\) diminue, ce qui potentiellement pose des problèmes de stabilité numérique dans un calcul éléments finis. Pour éviter cette situation, il est proposé de remplacer la fonction \(R(\xi)\), après discrétisation incrémentale des instants de chargement, par :

(4130)\[R_{n+1} = R_n\max\left\lbrace 1,\exp\left(-\beta\Delta\xi\right)\right\rbrace , \quad \text{avec}\quad R_0=p_{c0},\quad \Delta\xi = \xi_{n+1}-\xi_n\]

\(\left(R_n\right)_{n\geq 0}\) est une suite croissante. On ne tient alors compte que d’un écrouissage positif isotrope cumulé qui aurait été prédit par la fonction initiale \(R(\xi)\) (4112)-2.

Remarque :

Si des applications du modèle CSSM devaient nécessiter un intérêt à prédire un écrouissage négatif, on pourra au besoin lever la restriction imposée par (4130) dans des développements ultérieurs.

Paramètres

Classification

Onze paramètres sont nécesssaires à définir le modèle CSSM. Ceux-ci sont classifiés dans le Tableau 62.

Tableau 62 Paramètres.

Intervention	Appellation	Définition	Unité	Intervalle	Symbole
Élasticité	BulkModulus ShearModulus ShearModulusRatio	Module de compressibilité Module de cisaillement Rapport du module de cisaillement du composant 1 sur le module de cisaillement	[Pa] [Pa] [-]	\(\mathbb{R}^{+*}\) \(\mathbb{R}^{+*}\) \([0;1]\)	\(K\) \(\mu\) \(\rho\)
Composant 1	CritStateSlope InitCritPress IncoPlastIndex IsoHardRatio IsoHardIndex	Pente d’état critique Pression critique initiale Indice d’incompressibilité plastique Rapport de réduction homothétique du domaine d’élasticité initial Indice d’écrouissage par agrandissement homothétique du domaine d’élasticité initial	[-] [Pa] [-] [-] [-]	\(\mathbb{R}^{+*}\) \(\mathbb{R}^{+*}\) \(\mathbb{R}^{+}\) \([0;1[\) \(\mathbb{R}^{+}\)	\(M\) \(p_{c0}\) \(\beta\) \(\eta\) \(\omega\)
Composant 2	HypDistortion HypExponent MinCritPress	Distorsion de référence de la relation « hyperbolique modifiée » (4123) Paramètre de courbure de la relation « hyperbolique modifiée » (4123) Pression minimale pour laquelle l’état critique (4241) est atteignable	[-] [-] [Pa]	\(\mathbb{R}^{+*}\) \([0;1[\) \(\mathbb{R}^{+*}\)	\(\gamma_{\mathrm{hyp}}\) \(n_{\mathrm{hyp}}\) \(C\)

Cas particuliers

De par sa conception, le modèle CSSM est un modèle hybride qui combine deux modèles de type Cam-Clay modifié (composant 1) et d’Iwan (composant 2). Cette combinaison est réalisée par la pondération du comportement en cisaillement via le paramètre \(\rho\in [0;1]\). Ainsi les cas particuliers peuvent être déduits.

• Si \(\rho=1\) : le domaine d’élasticité initial du modèle CSSM est caractérisé par l’intersection \(\left\lbrace f_{0}\leq 0 \cap F_{1,0}\leq 0\right\rbrace\) où d’après (4113) et (4118) :

  (4131)\[f_{0} = \sqrt{\left(\frac{\sigma_{eq}}{M}\right)^2 + \left(\sigma_m+p_{c0}(1-\eta)\right)^2}-p_{c0}(1-\eta), \quad F_{1,0} = \left|\frac{R_1}{C}\left\langle \sigma_m+C\right\rangle\right|-R_1\]

  Il est facile de voir que \(\left\lbrace f_{0}\leq 0 \cap F_{1,0}\leq 0\right\rbrace =\left\lbrace f_{0}\leq 0\right\rbrace\) et que seul le composant 1 du modèle CSSM agit.

• Si \(\rho=0\) : le domaine d’élasticité initial du modèle CSSM est à l’opposé caractérisé par \(\left\lbrace f_{0}\leq 0 \cap F_{1,0}\leq 0\right\rbrace\) avec :

  (4132)\[f_{0} = \left|\sigma_m+p_{c0}(1-\eta)\right|-p_{c0}(1-\eta), \quad F_{1,0} = \sqrt{\left(\sigma_{eq}\right)^2 + \left(\frac{R_1}{C}\left\langle \sigma_m+C\right\rangle\right)^2} -R_1\]

  Dans ce cas tant qu’on ne s’intéresse qu’au comportement en cisaillement, seul le composant 2 du modèle CSSM agit.

Nous illustrons dans la r7.01.44-cas_tests que sélectionner \(\rho\in\lbrace 0,1\rbrace\) ne permet pas de prédire de manière satisfaisante à la fois le comportement d’un matériau sol sous chargements monotones et cycliques, allant des très faibles aux niveaux moyens de déformation. C’est donc la combinaison des deux composants qui améliore les prédictions, et justifie en fait l’intérêt du modèle CSSM.

Intégration numérique

Résolution implicite des équations de comportement dans MFront

Le modèle de comportement CSSM est intégré de manière implicite via l’outil MFront [HMPS15]. Cette intégration intervient sur l’incrément du tenseur des déformations totales \(\Delta \boldsymbol{\varepsilon}\) à l’instant actuel \(t_{n+1}\) en chaque point d’intégration, sachant l’état des variables internes et du tenseur des contraintes à l’instant précédent \(t_n\). Le vecteur des inconnues numériques est noté \(\boldsymbol{x}\). Ses composantes sont :

(4133)\[\begin{split}\begin{Bmatrix} \boldsymbol{x}_1 \\ \boldsymbol{x}_2 \\ x_3 \\ x_4 \\ x_5 \end{Bmatrix}= \begin{Bmatrix} \Delta\boldsymbol{\varepsilon}^e \\ \Delta\boldsymbol{\varepsilon}^x \\ \Delta\lambda \\ \Delta\lambda_m \\ \Delta\xi \end{Bmatrix},\quad \text{avec}\quad \Delta x_i = x_{i,n+1}-x_{i,n}\end{split}\]

On rappelle que les deux tenseurs des déformations élastiques \(\Delta\boldsymbol{\varepsilon}^e\) et \(\Delta\boldsymbol{\varepsilon}^x\) ont été définis (4111).

De plus, il convient d’observer ci-dessus que les incréments des multiplicateurs plastiques \(\left(\Delta\lambda_i\right)_{1\leq i\leq m-1}\) n’interviennent pas. En effet, grâce à la condition (4124), on établit que :

(4134)\[\forall 1\leq i\leq m-1,\quad \Delta\lambda_i = \frac{2}{3}\frac{\left \langle F_i\left(\boldsymbol{\sigma}_{n+1}-\boldsymbol{X}_{d,n+1} -H_d^i\boldsymbol{\alpha}_{d,n}^i-H_v^i\alpha_{v,n}^i\boldsymbol{I}\right)\right \rangle}{H_d^i}\]

et :

(4135)\[\frac{\partial F_i}{\partial \boldsymbol{A}^i}\Bigg{|}_{n+1} = \frac{\partial F_i}{\partial \boldsymbol{A}^i}\left(\boldsymbol{\sigma}_{n+1}-\boldsymbol{X}_{d,n+1} -H_d^i\boldsymbol{\alpha}_{d,n}^i-H_v^i\alpha_{v,n}^i\boldsymbol{I}\right)\]

Les incréments \(\left(\Delta\boldsymbol{\alpha}^i=\Delta\lambda_i\cfrac{\partial F_i}{\partial \boldsymbol{A}^i}\Bigg{|}_{n+1}\right)_{1\leq i\leq m-1}\) peuvent donc à chaque instant être exprimés en fonction des inconnues \(\Delta\boldsymbol{\varepsilon}^e\) et \(\Delta\boldsymbol{\varepsilon}^x\).

Remarque :

Le vecteur des inconnues numériques \(\boldsymbol{x}\) (4133) ne comprend pas la variable d’écrouissage scalaire \(\gamma\) car, suite à (4116)-1 et (4116)-2, celle-ci est déduite de \(\lambda\) et \(\xi\).

Système à résoudre dans le cas d’un chargement élastique

La prédiction élastique suppose l’incrément du tenseur des déformations totales élastique, ce qui siginifie \(\Delta\boldsymbol{\varepsilon}^e=\Delta\boldsymbol{\varepsilon}^x=\Delta\boldsymbol{\varepsilon}\). Si cette prédiction respecte les deux critères de plasticité \(f\) et \(F_m\), le résidu du système d’équations à résoudre, noté \(\boldsymbol{r}\), est :

(4136)\[\begin{split}\begin{Bmatrix} \boldsymbol{r}_1 \\ \boldsymbol{r}_2 \\ r_3 \\ r_4 \\ r_5 \end{Bmatrix} = \begin{Bmatrix} \Delta\boldsymbol{\varepsilon}^e - \Delta\boldsymbol{\varepsilon} + \left(\mathbb{J}+(1-\rho)\mathbb{K}\right):\sum_{i=1}^{m-1}\Delta\boldsymbol{\alpha}^i \\ \Delta\boldsymbol{\varepsilon}^x - \Delta\boldsymbol{\varepsilon}^e -\left(1-\rho\right)\mathbb{K}:\sum_{i=1}^{m-1}\Delta\boldsymbol{\alpha}^i\\ \Delta\lambda \\ \Delta\lambda_m \\ \Delta\xi\\ \end{Bmatrix}\end{split}\]

Ci-dessus, les deux premières lignes rendent compte des décompositions (4111) dans le cas où les incréments \(\Delta\boldsymbol{\varepsilon}^p\) et \(\Delta\boldsymbol{\alpha}^m\) tous tous deux nuls. De même, les trois dernières lignes sont cohérentes à la prédiction élastique si celle-ci est vérifiée (pas d’incrément des variables internes dont l’évolution est relative aux deux critères de plasticité \(f\) et \(F_m\)).

Système à résoudre dans le cas d’un chargement plastique

Si ni le critère de plasticité \(f\) ni \(F_m\) n’est vérifié après la prédiction élastique, le résidu du système d’équations à résoudre lors de l’intégration est donné par :

(4137)\[\begin{split}\begin{Bmatrix} \boldsymbol{r}_1 \\ \boldsymbol{r}_2 \\ r_3 \\ r_4 \\ r_5 \end{Bmatrix} = \begin{Bmatrix} \Delta\boldsymbol{\varepsilon}^e - \Delta\boldsymbol{\varepsilon} + \left(\mathbb{J}+\rho\mathbb{K}\right):\Delta\lambda\cfrac{\partial f}{\partial\boldsymbol{X}}\Big{|}_{n+1}+ \left(\mathbb{J}+(1-\rho)\mathbb{K}\right):\left(\Delta\lambda_m\cfrac{\partial F_m}{\partial \boldsymbol{A}^m}\Big{|}_{n+1} +\sum_{i=1}^{m-1}\Delta\boldsymbol{\alpha}^i\right) \\ \Delta\boldsymbol{\varepsilon}^x - \Delta\boldsymbol{\varepsilon}^e + \left(1-\rho\right)\mathbb{K}:\left(\Delta\lambda\cfrac{\partial f}{\partial\boldsymbol{X}}\Big{|}_{n+1} -\Delta\lambda_m\cfrac{\partial F_m}{\partial \boldsymbol{A}^m}\Big{|}_{n+1} -\sum_{i=1}^{m-1}\Delta\boldsymbol{\alpha}^i\right)\\ \cfrac{f_{n+1}}{K} \\ \cfrac{F_{m,n+1}}{K} \\ \Delta\xi - \Delta\lambda\cfrac{\partial f}{\partial p_c}\Big{|}_{n+1} \\ \end{Bmatrix}\end{split}\]

Dans le précédent système, les deux premières lignes décrivent les décompositions (4111) dans le cas où les incréments \(\Delta\varepsilon^p\) et \(\Delta\alpha^n\) sont a priori non nuls. Les deux suivantes expriment que les surfaces relatives aux deux critères de plasticité \(f\) et \(F_m\) sont atteintes (les équations sont adimensionnées par le module de compressibilité). La dernière rend compte de l’évolution de la variable d’écrouissage scalaire \(\xi\).

Remarque :

Le cas où la prédiction élastique respecte exclusivement le critère de plasticité \(f\) ou \(F_m\) peut être déduit aisément à partir des deux résidus (4136) et (4137).

La méthode itérative de Newton-Raphson est utilisée pour résoudre le système d’équations \(\boldsymbol{r}=\boldsymbol{0}\). 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 :

(4138)\[\begin{split}J = \begin{Bmatrix} \cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\boldsymbol{\varepsilon}^e} & \cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\boldsymbol{\varepsilon}^x} & \cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\lambda} & \cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\lambda_m} & \cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\xi} \\ \cfrac{\partial \boldsymbol{r}_2}{\partial \Delta\boldsymbol{\varepsilon}^e} & \cfrac{\partial \boldsymbol{r}_2}{\partial \Delta\boldsymbol{\varepsilon}^x} & \cfrac{\partial \boldsymbol{r}_2}{\partial \Delta\lambda} & \cfrac{\partial \boldsymbol{r}_2}{\partial \Delta\lambda_m} & \cfrac{\partial \boldsymbol{r}_2}{\partial \Delta\xi} \\ \cfrac{\partial r_3}{\partial \Delta\boldsymbol{\varepsilon}^e} & \cfrac{\partial r_3}{\partial \Delta\boldsymbol{\varepsilon}^x} & \cfrac{\partial r_3}{\partial \Delta\lambda} & \cfrac{\partial r_3}{\partial \Delta\lambda_m} & \cfrac{\partial r_3}{\partial \Delta\xi} \\ \cfrac{\partial r_4}{\partial \Delta\boldsymbol{\varepsilon}^e} & \cfrac{\partial r_4}{\partial \Delta\boldsymbol{\varepsilon}^x} & \cfrac{\partial r_4}{\partial \Delta\lambda} & \cfrac{\partial r_4}{\partial \Delta\lambda_m} & \cfrac{\partial r_4}{\partial \Delta\xi} \\ \cfrac{\partial r_5}{\partial \Delta\boldsymbol{\varepsilon}^e} & \cfrac{\partial r_5}{\partial \Delta\boldsymbol{\varepsilon}^x} & \cfrac{\partial r_5}{\partial \Delta\lambda} & \cfrac{\partial r_5}{\partial \Delta\lambda_m} & \cfrac{\partial r_5}{\partial \Delta\xi} \end{Bmatrix}\end{split}\]

Chaque bloc est obtenu par des calculs analytiques détaillés dans la r7.01.44-jacobien . La valeur du critère d’arrêt du Newton-Raphson est prise à \(10^{-14}\) (valeur 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, s’obtient par une règle de dérivées en chaîne :

\[\mathbb{D} = \frac{\partial \Delta \boldsymbol{\sigma}}{\partial \Delta\boldsymbol{\varepsilon}}=\mathbb{C} : \frac{\partial \Delta\boldsymbol{\varepsilon}^e}{\partial \Delta\boldsymbol{\varepsilon}}\]

Le premier terme \(\mathbb{C}=3K\mathbb{J}+2\mu\mathbb{K}\) désigne le tenseur d’élasticité linéaire isotrope. Le second terme \(\cfrac{\partial \Delta\boldsymbol{\varepsilon}^e}{\partial \Delta\boldsymbol{\varepsilon}}\) s’écrit quant à lui :

\[\cfrac{\partial \Delta\boldsymbol{\varepsilon}^e}{\partial \Delta\boldsymbol{\varepsilon}}= 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], [Helf20].

Variables internes

Les variables internes sont résumées dans le Tableau 63.

Tableau 63 Variables internes du modèle.

Appellation	Définition	Symbole	Composantes 2D	Composantes 3D
ElasticStrain	Tenseur des déformations élastiques	\(\boldsymbol{\varepsilon}^e\)	V1-V4	V1-V6
ElasticStrain1	Tenseur des déformations élastiques dans le composant 1	\(\boldsymbol{\varepsilon}^x\)	V5-V8	V7-V12
CumPlastStrain1	Déformation plastique cumulée dans le composant 1	\(\lambda\)	V9	V13
CumPlastStrain2	Déformation plastique cumulée dans le composant 2 relatif au critère \(F_m\)	\(\lambda_m\)	V10	V14
VolPlastStrain1	Déformation plastique volumique dans le composant 1	\(\xi\)	V11	V15

Les variables auxiliaires supplémentaires dans le champ des variables internes au sens de Code_Aster sont présentées dans le Tableau 64.

Tableau 64 Variables auxiliaires du modèle.

Appellation	Définition	Symbole	Composantes 2D	Composantes 3D
OldBkStress	Tenseurs d’écrouissage cinématique à l’instant \(t_n\)	\(\left(H_d^i\boldsymbol{\alpha}_{d,n}^i+H_v^i\alpha_{v,n}^i\boldsymbol{I}\right)_{1\leq i \leq N-1}\)	V12-47	V16-V69
OldYieldSize	Paramétrisation à l’instant \(t_n\) du critère de plasticité (4112) vérifiant la relation (4130)	\(R_n\)	V48	V70
PlasticIndicator	0 si \(f<0\) et \(F_m<0\), 1 si \(f=0\) et \(F_m<0\), 2 si \(f<0\) et \(F_m=0\), 3 si \(f=0\) et \(F_m=0\)	\(I_p\)	V49	V71

Vérification et validation

Cas-tests

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

Tableau 65 Cas-tests et documentations associées.

Référence cas-test	Référence documentation	Description
ssnv160f	[v6.04.160]	Essai de compression isotrope
ssnv205c	[v6.04.205]	Essai de cisaillement cyclique
ssnv207c	[v6.04.207]	Essai de cisaillement cyclique avec micro-décharge
comp012i	[v6.07.112]	Test de compatibilité avec la commande CALC_GEO_MECA
wtnv122e	[v7.31.122]	Essai de compression triaxiale non drainée

Exemples de réponses au point matériel

Les Fig. 270, Fig. 271 et Fig. 272 présentent quelques réponses obtenues à l’aide de l’exécutable MTest (réponses au point matériel) :

• Essai de compression isotrope.

• Essais de compression triaxiale pour plusieurs contraintes de confinement.

• Essais de cisaillement cyclique pour plusieurs amplitudes de déformations.

Les paramètres utilisés dans ces simulations sont regroupés dans le Tableau 66.

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

Intervention	Appellation	Définition	Symbole	Valeur
Élasticité	BulkModulus ShearModulus ShearModulusRatio	Module de compressibilité total Module de cisaillement total Rapport du module de cisaillement du composant 1 sur le module de cisaillement total	\(K\) \(\mu\) \(\rho\)	516 MPa 238 MPa 0.1
Composant 1	CritStateSlope InitCritPress IncoPlastIndex IsoHardRatio IsoHardIndex	Pente d’état critique Pression critique initiale Indice d’incompressibilité plastique Rapport de réduction homothétique du domaine d’élasticité initial Indice d’écrouissage par agrandissement homothéTique du domaine d’élasticité initial	\(M\) \(p_{c0}\) \(\beta\) \(\eta\) \(\omega\)	\(1.38\) \(100\) kPa \(30\) \(0.99\) \(32\)
Composant 2	HypDistortion HypExponent MinCritPress	Distorsion de référence de la relation « hyperbolique modifiée » (4123) Paramètre de courbure de la relation « hyperbolique modifiée » (4123) Pression minimale pour laquelle l’état critique (4241) est atteignable	\(\gamma_{\mathrm{hyp}}\) \(n_{\mathrm{hyp}}\) \(C\)	\(2.10^{-4}\) \(0.78\) \(448\) kPa

La Fig. 270 présente en ligne continue la réponse de l’essai de compression isotrope pour les valeurs du Tableau 66. Dans la phase de charge, la réponse est élastique jusqu’à atteindre la pression \(-\sigma_m=2p_{c0}(1-\eta)=0.2\) kPa (\(I_p=0\)) puis devient irréversible, avec activation uniquement de la plasticité du premier composant (\(I_p=1\)). La décharge est élastique. En comparaison, la ligne discontinue montre la réponse dans le cas où \(\eta=0\), avec les autres paramètres inchangés. Dans cette situation, la réponse reste élastique jusqu’à une pression de \(-\sigma_m=200\) kPa. Ainsi, l’utilisation du paramètre \(\eta\) permet d’introduire une compaction initiale plus forte pour des faibles niveaux de pression, où la déformation volumique \(\varepsilon_v\) est plus élevée en compression. Cette possibilité est intéressante pour modéliser le comportement fortement contractant des sols à porosité élevée.

Fig. 270 Essai de compression isotrope piloté par la pression \(-\sigma_m\).

La Fig. 271 présente les réponses d’essais de compression triaxiale pour trois niveaux de confinement \(p_{\mathrm{conf}}\). On observe que plus le confinement est élevé, plus la pression \(-\sigma_m\) atteinte asympotiquement est importante. De plus, lorsque cette pression dépasse \(C=448\) kPa, la déformation volumique \(\varepsilon_v\) se stabilise. En revanche, lorsque cette pression est inférieure, comme on le voit clairement pour \(p_{\mathrm{conf}}=100\) kPa, la déformation volumique continue d’augmenter, indiquant un comportement dilatant. Cette observation est cohérente avec les équations de l’état critique référencées dans (4241)-4.

Fig. 271 Essais de compression triaxiale à trois contraintes de confinement \(-\sigma_{xx}=-\sigma_{yy}=p_{\mathrm{conf}}\).

La Fig. 272 présente, en ligne continue, la réponse des essais de cisaillement à pression constante \(-\sigma_{m}=200\) kPa pour les valeurs du Tableau 66, avec plusieurs amplitudes de distorsion \(\gamma=2\varepsilon_{xy}\in[10^{-6};10^{-2}]\). Pour chaque valeur de distorsion, un cycle est réalisé selon la cronologie \(\gamma\rightarrow-\gamma\rightarrow\gamma\). La contrainte de cisaillement \(\sigma_{xy}\) présente des boucles d’hystérésis, qui ne sont pas fermées en raison du fait qu’un seul cycle ne permet pas de stabiliser l’écrouissage associé aux composantes volumiques des variables internes. Les réponses sont comparées à celles prédites avec \(\rho=0\) et \(\rho=1\), les autres paramètres restant inchangés. Les prédictions du modèle CSSM rejoignent alors celles du second composant seul (cas \(\rho=0\)), et s’écartent plus significativement de celles du premier composant seul (cas \(\rho=1\)) sur les évolutions du module de cisaillement sécant normalisé \(\mu_{\mathrm{secant}}/\mu\) et de l’amortissement réduit.

Fig. 272 Essais de cisaillement cyclique à pression constante \(-\sigma_{m}=200\) kPa.

Exemple de comparaison expérimentale

Les Fig. 273, Fig. 274 et Fig. 275 confrontent les prédictions du modèle CSSM avec des données expérimentales sur des matériaux en enrochement lâche provenant de deux essais de compression triaxiale aux confinements \(p_{\mathrm{conf}}=100\) kPa et \(p_{\mathrm{conf}}=200\) kPa, ainsi que d’essais de cisaillement cyclique.

Sur les essais de compression triaxiale, le modèle CSSM a été calibré « au mieux » en ajustant les cinq paramètres \(p_{c0},M,\beta,\omega\) et \(C\), en supposant les autres paramètres du Tableau 66 et en considérant trois cas distincts ayant déjà fait l’objet de précédentes comparaisons :

• \(\rho=0\leftrightarrow\) effet du composant 2 du modèle seulement sur le comportement en cisaillement.

• \(\rho=1\leftrightarrow\) effet du composant 1 du modèle seulement.

• \(\rho=0.1\leftrightarrow\) effet des deux composants du modèle.

Le dernier cas correspond à considérer que la rigidité de cisaillement dans le composant 1 \(\rho\mu\) (rigidité « statique ») est dix fois plus faible que la rigidité de cisaillement totale \(\mu\) (rigidité « dynamique »). Ce rapport de dix n’est pas déraisonnable, comparé aux rapports observés entre les modules dynamique et statique pour les matériaux sols.

La Fig. 273 confronte les prédictions pour le cas où \(\rho=0\) avec les données expérimentales. Il est clair que le composant 2 du modèle à lui seul ne parvient pas à reproduire de manière réaliste les évolutions de la contrainte équivalente de von Mises \(\sigma_{eq}\) et de la déformation volumique \(\varepsilon_v\) observées lors des essais de compression triaxiale. En revanche, le composant 2 offre une représentation assez précise du module de cisaillement sécant normalisé \(\mu_{\mathrm{secant}}/\mu\) lors d’essais de cisaillement. Quant à l’amortissement réduit, la prédiction surévalue progressivement l’évolution expérimentale avec l’augmentation du niveau de distorsion \(\gamma\). Cette prédiction est très similaire à celle du modèle d’Iwan mais également de Hujeux [r7.01.23].

Fig. 273 Comparaison des prédictions dans le cas \(\rho=0\) aux données expérimentales.

La Fig. 274 confronte les prédictions pour le cas où \(\rho=1\) avec les données expérimentales. Dans ce cas, le composant 1 du modèle améliore sensiblement l’évolution de la contrainte équivalente de von Mises lors des essais de compression triaxiale. Toutefois, seule la déformation volumique pour le confinement \(p_{\mathrm{conf}}=200\) kPa est correctement reproduite. Pour le confinement \(p_{\mathrm{conf}}=100\) kPa, la prédiction manque de dilatance. En ce qui concerne les essais de cisaillement, les prédictions du module de cisaillement sécant normalisé \(\mu_{\mathrm{secant}}/\mu\) et de l’amortissement réduit sont largement détériorées par rapport à celles obtenues pour le cas \(\rho=0\).

Fig. 274 Comparaison des prédictions dans le cas \(\rho=1\) aux données expérimentales.

La Fig. 275 présente enfin les prédictions dans le cas où \(\rho=0.1\). Cette combinaison des composants 1 et 2 du modèle permet :

• D’améliorer l’évolution observée de la dilatance pour le confinement \(p_{\mathrm{conf}}=100\) kPa lors des essais de compression triaxiale, tout en conservant une représentation correcte de la contrainte équivalente de von Mises. Cette représentation est similaire à celle obtenue avec le composant 1 seul dans le cas où \(\rho=0\).

• D’obtenir des prédictions proches de celles du composant 2 pour les essais de cisaillement lorsque \(\rho=1\).

En résumé, la combinaison des deux composants du modèle permet de tirer parti des meilleures caractéristiques de chaque composant individuel, offrant ainsi une représentation améliorée pour les deux types d’essais monotone et cyclique.

Fig. 275 Comparaison des prédictions dans le cas \(\rho=0.1\) aux données expérimentales.

Remarque :

Les comparaisons précédentes valident la pertinence du modèle CSSM sur des chargements avec cisaillement, avec une calibration des paramètres \(p_{c0},M,\beta,\omega\) et \(C\) sur les essais de compression triaxiale. Néanmoins, pour évaluer plus globalement,et rigoureusement, les capacités du modèle, il sera nécessaire de le comparer à des essais de compression isotrope. Cette étape permettra d’estimer plus précisément le paramètre \(\eta\), ici pris proche de un dans le Tableau 66, et affiner la valeur de \(\omega\), qui influencent significativement la compaction initiale, comme nous l’avons aperçu sur la Fig. 270 avec \(\eta\in\lbrace 0,0.99\rbrace\).

Annexe : expression de la matrice jacobienne

Dans le cas où ni le critère de plasticité \(f\) ni \(F_m\) n’est vérifié après la prédiction élastique, on rappelle l’expression du système d’équations non-linéaires (4137) et la définition de la matrice jacobienne (4249) :

(4139)\[\begin{split}\begin{Bmatrix} \boldsymbol{r}_1 \\ \boldsymbol{r}_2 \\ r_3 \\ r_4 \\ r_5 \end{Bmatrix} = \begin{Bmatrix} \Delta\boldsymbol{\varepsilon}^e - \Delta\boldsymbol{\varepsilon} + \left(\mathbb{J}+\rho\mathbb{K}\right):\Delta\lambda\cfrac{\partial f}{\partial\boldsymbol{X}}\Big{|}_{n+1}+ \left(\mathbb{J}+(1-\rho)\mathbb{K}\right):\left(\Delta\lambda_m\cfrac{\partial F_m}{\partial \boldsymbol{A}^m}\Big{|}_{n+1} +\sum_{i=1}^{m-1}\Delta\boldsymbol{\alpha}^i\right) \\ \Delta\boldsymbol{\varepsilon}^x - \Delta\boldsymbol{\varepsilon}^e + \left(1-\rho\right)\mathbb{K}:\left(\Delta\lambda\cfrac{\partial f}{\partial\boldsymbol{X}}\Big{|}_{n+1} -\Delta\lambda_m\cfrac{\partial F_m}{\partial \boldsymbol{A}^m}\Big{|}_{n+1} -\sum_{i=1}^{m-1}\Delta\boldsymbol{\alpha}^i\right)\\ \cfrac{f_{n+1}}{K} \\ \cfrac{F_{m,n+1}}{K} \\ \Delta\xi - \Delta\lambda\cfrac{\partial f}{\partial p_c}\Big{|}_{n+1} \\ \end{Bmatrix}\end{split}\]

(4140)\[\begin{split}J = \begin{Bmatrix} \cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\boldsymbol{\varepsilon}^e} & \cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\boldsymbol{\varepsilon}^x} & \cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\lambda} & \cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\lambda_m} & \cfrac{\partial \boldsymbol{r}_1}{\partial \Delta\xi} \\ \cfrac{\partial \boldsymbol{r}_2}{\partial \Delta\boldsymbol{\varepsilon}^e} & \cfrac{\partial \boldsymbol{r}_2}{\partial \Delta\boldsymbol{\varepsilon}^x} & \cfrac{\partial \boldsymbol{r}_2}{\partial \Delta\lambda} & \cfrac{\partial \boldsymbol{r}_2}{\partial \Delta\lambda_m} & \cfrac{\partial \boldsymbol{r}_2}{\partial \Delta\xi} \\ \cfrac{\partial r_3}{\partial \Delta\boldsymbol{\varepsilon}^e} & \cfrac{\partial r_3}{\partial \Delta\boldsymbol{\varepsilon}^x} & \cfrac{\partial r_3}{\partial \Delta\lambda} & \cfrac{\partial r_3}{\partial \Delta\lambda_m} & \cfrac{\partial r_3}{\partial \Delta\xi} \\ \cfrac{\partial r_4}{\partial \Delta\boldsymbol{\varepsilon}^e} & \cfrac{\partial r_4}{\partial \Delta\boldsymbol{\varepsilon}^x} & \cfrac{\partial r_4}{\partial \Delta\lambda} & \cfrac{\partial r_4}{\partial \Delta\lambda_m} & \cfrac{\partial r_4}{\partial \Delta\xi} \\ \cfrac{\partial r_5}{\partial \Delta\boldsymbol{\varepsilon}^e} & \cfrac{\partial r_5}{\partial \Delta\boldsymbol{\varepsilon}^x} & \cfrac{\partial r_5}{\partial \Delta\lambda} & \cfrac{\partial r_5}{\partial \Delta\lambda_m} & \cfrac{\partial r_5}{\partial \Delta\xi} \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 ci-dessous.

On notera pour cela :

• Les directions d’écoulement de \(\Delta\boldsymbol{\varepsilon}^p,\left(\Delta\boldsymbol{\alpha}^i\right)_{1\leq i\leq m}\) et \(\Delta\xi\) :

(4141)\[\begin{split}\begin{align} &\boldsymbol{N}_{n+1} = \frac{\partial f}{\partial \boldsymbol{X}}\Big{|}_{n+1} = \frac{\cfrac{3}{2M^2}\boldsymbol{Y}_{d,n+1}+Y_{m,n+1}\cfrac{\boldsymbol{I}}{3}}{T_{eq}(\boldsymbol{Y}_{n+1})} , \quad \text{avec} \quad \boldsymbol{Y}=\boldsymbol{X}+\left(p_c-S\right)\boldsymbol{I},\quad T_{eq}(\boldsymbol{Y}) = \sqrt{\left(\frac{Y_{eq}}{M}\right)^2+Y_m^2}\\ &\boldsymbol{M}^i_{n+1} = \frac{\partial F_i}{\partial \boldsymbol{A}^i}\Big{|}_{n+1} = \frac{\cfrac{3}{2}\boldsymbol{A}^i_{d,n+1}+\left(\cfrac{R_i}{C}\right)^2\left\langle A^i_{m,n+1}+C\right\rangle\cfrac{\boldsymbol{I}}{3}}{T^i_{eq}(\boldsymbol{A}^i_{n+1})} , \quad \text{avec} \quad T^i_{eq}(\boldsymbol{A}^i) = \sqrt{\left(A^i_{eq}\right)^2+\left(\cfrac{R_i}{C}\left\langle A^i_{m}+C\right\rangle\right)^2}\\ &N_{v,n+1} = \frac{\partial f}{\partial p_c}\Big{|}_{n+1} = \boldsymbol{I}:\boldsymbol{N}_{n+1} \\ \end{align}\end{split}\]

• Leurs trois dérivées suivantes :

(4142)\[\begin{split}\begin{align} &\boldsymbol{N}_{n+1,\boldsymbol{Y}_{n+1}} = \frac{\partial \boldsymbol{N}}{\partial \boldsymbol{Y}}\Big{|}_{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{Y}_{n+1})}\\ &\boldsymbol{M}^i_{n+1,\boldsymbol{A}^i_{n+1}} = \frac{\partial \boldsymbol{M}^i}{\partial \boldsymbol{A}^i}\Big{|}_{n+1} = \frac{\cfrac{3}{2}\mathbb{K}+\left(\cfrac{R_i}{C}\right)^2\mathbf{H}\left(A^i_{m,n+1}+C\right)\cfrac{\mathbb{J}}{3}-\boldsymbol{M}^i_{n+1}\otimes\boldsymbol{M}^i_{n+1}}{T^i_{eq}(\boldsymbol{A}^i_{n+1})}\\ &N_{v,n+1,\boldsymbol{Y}_{n+1}} = \frac{\partial N_{v}}{\partial \boldsymbol{Y}}\Big{|}_{n+1} = \frac{\cfrac{\boldsymbol{I}}{3}-N_{v,n+1}\boldsymbol{N}_{n+1}}{T_{eq}(\boldsymbol{Y}_{n+1})} \end{align}\end{split}\]

• Les dérivées des forces \(\boldsymbol{Y}=\boldsymbol{Y}(\boldsymbol{\varepsilon}^x,\lambda,\xi)\) et \(\left(\boldsymbol{A}^i\right)_{1\leq i \leq m}=\left(\boldsymbol{A}^i(\boldsymbol{\varepsilon}^e,\boldsymbol{\varepsilon}^x)\right)_{1\leq i \leq m}\) :

(4143)\[\begin{split}\begin{align} &\boldsymbol{Y}_{n+1,\boldsymbol{\varepsilon}^x} = \frac{\partial \boldsymbol{Y}}{\partial\Delta\boldsymbol{\varepsilon}^x}\Big{|}_{n+1} = \left(\mathbb{J}+\rho\mathbb{K}\right):\mathbb{C}\\ &\boldsymbol{Y}_{n+1,\lambda} = \frac{\partial \boldsymbol{Y}}{\partial\Delta\lambda}\Big{|}_{n+1} = \omega S_{n+1}\boldsymbol{I}\\ &\boldsymbol{Y}_{n+1,\xi} = \frac{\partial \boldsymbol{Y}}{\partial\Delta\xi}\Big{|}_{n+1} = -\left(\omega S_{n+1}+\beta p_{c,n+1}\right)\boldsymbol{I}\\ &\boldsymbol{A}^i_{n+1,\boldsymbol{\varepsilon}^e} = \frac{\partial \boldsymbol{A}^i}{\partial\Delta\boldsymbol{\varepsilon}^e}\Big{|}_{n+1} = \mathbb{C} \\ &\boldsymbol{A}^i_{n+1,\boldsymbol{\varepsilon}^x} = \frac{\partial \boldsymbol{A}^i}{\partial\Delta\boldsymbol{\varepsilon}^x}\Big{|}_{n+1} = -\rho\mathbb{K}:\mathbb{C} \end{align}\end{split}\]

• Les dérivées de tous les \(\left(\Delta\boldsymbol{\alpha}^i=\Delta\lambda_i\boldsymbol{M}_{n+1}^i\right)_{1\leq i\leq m-1}\) :

(4144)\[\begin{split}\begin{align} \frac{\partial \Delta\boldsymbol{\alpha}^i}{\partial \Delta\boldsymbol{\varepsilon}^e}=\chi_i\left(\Delta\lambda_i\boldsymbol{M}^i_{n+1,\boldsymbol{A}^i_{n+1}}+\frac{2}{3H_d^i}\boldsymbol{M}^i_{n+1}\otimes \boldsymbol{M}^i_{n+1}\right):\boldsymbol{A}^i_{n+1,\boldsymbol{\varepsilon}^e}\\ \frac{\partial \Delta\boldsymbol{\alpha}^i}{\partial \Delta\boldsymbol{\varepsilon}^x}=\chi_i\left(\Delta\lambda_i\boldsymbol{M}^i_{n+1,\boldsymbol{A}^i_{n+1}}+\frac{2}{3H_d^i}\boldsymbol{M}^i_{n+1}\otimes \boldsymbol{M}^i_{n+1}\right):\boldsymbol{A}^i_{n+1,\boldsymbol{\varepsilon}^x} \end{align}\end{split}\]

où \(\chi_i=1\) si \(\left \langle F_i\left(\boldsymbol{\sigma}_{n+1}-\boldsymbol{X}_{d,n+1} -H_d^i\boldsymbol{\alpha}_{d,n}^i-H_v^i\alpha_{v,n}^i\boldsymbol{I}\right)\right \rangle>0\), zéro sinon.

Première ligne

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

(4145)\[\begin{split}\begin{align} &\frac{\partial \boldsymbol{r}_1}{\partial \Delta\boldsymbol{\varepsilon}^e} = \mathbb{I}+\left(\mathbb{J}+(1-\rho)\mathbb{K}\right):\left(\Delta\lambda_m\boldsymbol{M}^m_{n+1,\boldsymbol{A}^m_{n+1}}:\boldsymbol{A}^m_{n+1,\boldsymbol{\varepsilon}^e} +\sum_{i=1}^{m-1}\frac{\partial \Delta\boldsymbol{\alpha}^i}{\partial \Delta\boldsymbol{\varepsilon}^e}\right)\\ &\frac{\partial \boldsymbol{r}_1}{\partial \Delta\boldsymbol{\varepsilon}^x} = \left(\mathbb{J}+\rho\mathbb{K}\right):\Delta\lambda\boldsymbol{N}_{n+1,\boldsymbol{Y}_{n+1}}:\boldsymbol{Y}_{n+1,\boldsymbol{\varepsilon}^x} +\left(\mathbb{J}+(1-\rho)\mathbb{K}\right):\left(\Delta\lambda_m\boldsymbol{M}^m_{n+1,\boldsymbol{A}^m_{n+1}}:\boldsymbol{A}^m_{n+1,\boldsymbol{\varepsilon}^x} +\sum_{i=1}^{m-1}\frac{\partial \Delta\boldsymbol{\alpha}^i}{\partial \Delta\boldsymbol{\varepsilon}^x}\right)\\ &\frac{\partial \boldsymbol{r}_1}{\partial \Delta\lambda}=\left(\mathbb{J}+\rho\mathbb{K}\right): \left(\boldsymbol{N}_{n+1}+\Delta\lambda\boldsymbol{N}_{n+1,\boldsymbol{Y}_{n+1}}:\boldsymbol{Y}_{n+1,\lambda}\right)\\ &\frac{\partial \boldsymbol{r}_1}{\partial \Delta\lambda_m}=\left(\mathbb{J}+(1-\rho)\mathbb{K}\right):\boldsymbol{M}^m_{n+1} \\ &\frac{\partial \boldsymbol{r}_1}{\partial \Delta\xi}=\left(\mathbb{J}+\rho\mathbb{K}\right): \Delta\lambda\boldsymbol{N}_{n+1,\boldsymbol{Y}_{n+1}}:\boldsymbol{Y}_{n+1,\xi} \end{align}\end{split}\]

Deuxième ligne

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

(4146)\[\begin{split}\begin{align} &\frac{\partial r_2}{\partial \Delta\boldsymbol{\varepsilon}^e} = -\mathbb{I}-(1-\rho)\mathbb{K}: \left(\Delta\lambda_m\boldsymbol{M}^m_{n+1,\boldsymbol{A}^m_{n+1}}:\boldsymbol{A}^m_{n+1,\boldsymbol{\varepsilon}^e}+\sum_{i=1}^{m-1}\frac{\partial \Delta\boldsymbol{\alpha}^i}{\partial \Delta\boldsymbol{\varepsilon}^e}\right)\\ &\frac{\partial r_2}{\partial \Delta\boldsymbol{\varepsilon}^x} = \mathbb{I}+(1-\rho)\mathbb{K}:\left( \Delta\lambda \boldsymbol{N}_{n+1,\boldsymbol{Y}_{n+1}}:\boldsymbol{Y}_{n+1,\boldsymbol{\varepsilon}^x}- \Delta\lambda_m\boldsymbol{M}^m_{n+1,\boldsymbol{A}^m_{n+1}}:\boldsymbol{A}^m_{n+1,\boldsymbol{\varepsilon}^x}-\sum_{i=1}^{m-1}\frac{\partial \Delta\boldsymbol{\alpha}^i}{\partial \Delta\boldsymbol{\varepsilon}^x}\right)\\ &\frac{\partial r_2}{\partial \Delta\lambda} = (1-\rho)\mathbb{K}:\left(\boldsymbol{N}_{n+1}+\Delta\lambda\boldsymbol{N}_{n+1,\boldsymbol{Y}_{n+1}}:\boldsymbol{Y}_{n+1,\lambda}\right)\\ &\frac{\partial r_2}{\partial \Delta\lambda_m} =-(1-\rho)\mathbb{K}:\boldsymbol{M}^m_{n+1} \\ &\frac{\partial r_2}{\partial \Delta\xi}=(1-\rho)\mathbb{K}:\Delta\lambda\boldsymbol{N}_{n+1,\boldsymbol{Y}_{n+1}}:\boldsymbol{Y}_{n+1,\xi} \end{align}\end{split}\]

Troisième ligne

La troisième ligne du système est :

(4147)\[\begin{split}\begin{align} &\frac{\partial r_3}{\partial \Delta\boldsymbol{\varepsilon}^e} = \boldsymbol{0}\\ &\frac{\partial r_3}{\partial \Delta\boldsymbol{\varepsilon}^x} = \cfrac{\boldsymbol{N}_{n+1}:\boldsymbol{Y}_{n+1,\boldsymbol{\varepsilon}^x}}{K}\\ &\frac{\partial r_3}{\partial \Delta\lambda}=\cfrac{\boldsymbol{N}_{n+1}:\boldsymbol{Y}_{n+1,\lambda}-\omega S_{n+1}}{K}\\ &\frac{\partial r_3}{\partial \Delta\lambda_m}=0\\ &\frac{\partial r_3}{\partial \Delta\xi}=\cfrac{\boldsymbol{N}_{n+1}:\boldsymbol{Y}_{n+1,\xi}+\omega S_{n+1}+\beta R_{n+1}\mathbf{H}(\Delta\xi)}{K}\\ \end{align}\end{split}\]

Quatrième ligne

La quatrième ligne :

(4148)\[\begin{split}\begin{align} &\frac{\partial r_4}{\partial \Delta\boldsymbol{\varepsilon}^e} = \cfrac{\boldsymbol{M}^m_{n+1}:\boldsymbol{A}^m_{n+1,\boldsymbol{\varepsilon}^e}}{K}\\ &\frac{\partial r_4}{\partial \Delta\boldsymbol{\varepsilon}^x} = \cfrac{\boldsymbol{M}^m_{n+1}:\boldsymbol{A}^m_{n+1,\boldsymbol{\varepsilon}^x}}{K}\\ &\frac{\partial r_4}{\partial \Delta\lambda}=0\\ &\frac{\partial r_4}{\partial \Delta\lambda_m}=0\\ &\frac{\partial r_4}{\partial \Delta\xi}=0\\ \end{align}\end{split}\]

Cinquième ligne

La cinquième et dernière ligne :

(4149)\[\begin{split}\begin{align} &\frac{\partial r_5}{\partial \Delta\boldsymbol{\varepsilon}^e} = \boldsymbol{0}\\ &\frac{\partial r_5}{\partial \Delta\boldsymbol{\varepsilon}^x} = -\Delta\lambda N_{v,n+1,\boldsymbol{Y}_{n+1}}\boldsymbol{Y}_{n+1,\boldsymbol{\varepsilon}^x}\\ &\frac{\partial r_5}{\partial \Delta\lambda}=-N_{v,n+1}-\Delta\lambda N_{v,n+1,\boldsymbol{Y}_{n+1}}\boldsymbol{Y}_{n+1,\lambda}\\ &\frac{\partial r_5}{\partial \Delta\lambda_m}=0\\ &\frac{\partial r_5}{\partial \Delta\xi}=1-\Delta\lambda N_{v,n+1,\boldsymbol{Y}_{n+1}}\boldsymbol{Y}_{n+1,\xi}\\ \end{align}\end{split}\]

Bibliographie

[Iwan66]

Iwan, W.D. A Distributed-Element Model for Hysteresis and Its Steady-State Dynamic Response. Journal of Applied Mechanics, 1966, Volume 33, Issue 4, pp.893-900.

[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.

[HMPS15] (1,2,3)

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.

[Helf20]

Helfer, T. Assisted computation of the consistent tangent operator of behaviours integrated using an implicit scheme. Theory and implementation in MFront. Technical report, 2020.