r7.01.28 Loi de Mohr-Coulomb

Résumé:

Ce document présente la méthode de résolution de la loi de Mohr-Coulomb dans Code_Aster .

Formulation en termes de contraintes principales

Cette formulation n’est valable que sous l’hypothèse d’une isotropie du matériau [, ]. En effet, cette condition est nécessaire pour garantir que la méthode de retour radial préserve les directions principales. Son intérêt réside dans le fait qu’elle simplifie l’écriture des équations et autorise de ce fait des méthodes de résolution très performantes (car quasi-analytiques).

Le comportement élastique est purement linéaire.

La surface de charge est caractérisée par six plans dans l’espace des contraintes principales \(\left({\sigma}_{1},{\sigma}_{2},{\sigma}_{3}\right)\) . Chacun de ces plans est caractérisé par une équation du type:

(3526)\[{F}_{13}\left({\sigma}_{1},{\sigma}_{3}\right)={\sigma}_{1}-{\sigma}_{3}+\left({\sigma}_{1}+{\sigma}_{3}\right)\sin\varphi -\mathrm{2c}\cos\varphi =0\]

Où \(\phi\) et \(c\) sont des données matériau et caractérisent respectivement l’angle de frottement interne et la cohésion du matériau.

La loi est non associée et le potentiel d’écoulement plastique \({G}_{13}\) associé à la surface de charge \({F}_{13}\) s’écrit de la même façon:

(3527)\[{G}_{13}\left({\sigma}_{1},{\sigma}_{3}\right)={\sigma}_{1}-{\sigma}_{3}+\left({\sigma}_{1}+{\sigma}_{3}\right)\sin\psi -\mathrm{2c}\cos\psi\]

Où \(\psi\) est une donnée matériau et caractérise l’angle de dilatance du matériau.

Lorsque \(\psi =\varphi\) , la loi d’écoulement plastique devient associée.

Une représentation graphique de la surface de charge de Mohr-Coulomb dans l’espace des contraintes principales se trouve sur la figure et dans le plan \(\pi\) sur la figure .

On observe que les six plans s’intersectent deux à deux suivant une arête anguleuse, et se rejoignant au sommet du cône caractérisé par l’équation:

(3528)\[p=c\mathrm{cot}\left(\varphi \right)\]

Ces arêtes, au nombre de six, ainsi que le sommet du cône, forment des singularités qui posent des problèmes d’intégration numérique, car les dérivées de la surface de charge ne sont pas définies à ces endroits. On discutera plus loin des méthodes permettant de résoudre cette difficulté.

Figure 2-1: Représentation de la surface de charge de Mohr-Coulomb dans l’espace tridimensionnelle des contraintes principales

Figure 2-2: Représentation de la surface de charge de Mohr-Coulomb dans le plan \(\pi\) des déviateurs des contraintes (tout vecteur représenté dans ce plan correspond à une contrainte déviatorique).

Intégration locale de la loi de Mohr-Coulomb

Le taux de déformation plastique est donné à l’aide de la formule de Koiter:

(3529)\[d{\epsilon}^{p}=\sum_{j=1}^{m}d{\lambda}_{j}\frac{\partial {G}_{j}}{\partial \sigma }=\sum_{j=1}^{m}d{\lambda}_{j}{{n}_{G}}_{j}\]

Où \(m\) caractérise le nombre de mécanismes actifs, égal à un, deux ou six suivant les situations suivantes:

• la contrainte finale se situe à l’intérieur de la surface de charge, le point est régulier et \(m=1\) ;

• la contrainte finale se situe sur une arête du cône, le point est singulier et \(m=2\) ;

• la contrainte finale ne se situe ni à l’intérieur de la surface de charge ni sur une arête. Elle est alors projetée au sommet du cône, le point est singulier et \(m=6\) ;

La contrainte finale \({\sigma}^{+}\) est calculée à partir d’une prédiction élastique notée \({\sigma}^{e}\) et d’une correction \({\sigma}^{c}=C.d{\epsilon}^{p}\) de sorte que:

(3530)\[{\sigma}^{+}={\sigma}^{e}-d{\sigma}^{c}={\sigma}^{e}-C.d{\epsilon}^{p}={\sigma}^{e}-\sum_{j=1}^{m}d{\lambda}_{j}C.{{n}_{G}}_{j}\]

Les multiplicateurs plastiques \(d{\lambda}_{j}\) sont calculés en injectant l’équation () dans l’équation (), ce qui donne:

(3531)\[\sum_{j=1}^{m}d{\lambda}_{j}\left[{(C.{{n}_{G}}_{j})}_{1}-{(C.{{n}_{G}}_{j})}_{3}+({(C.{{n}_{G}}_{j})}_{1}+{(C.{{n}_{G}}_{j})}_{3})s\right]={\sigma}_{1}^{e}-{\sigma}_{3}^{e}+({\sigma}_{1}^{e}+{\sigma}_{3}^{e})s-\mathrm{2c}\cos\phi\]

Dans ce qui suit, on détaille les expressions correspondant aux différentes situations mentionnées plus haut. La procédure de résolution est rappelée dans le synoptique de la figure .

Cas où un seul mécanisme est actif

On détecte que la prédiction active un critère quand:

(3532)\[{F}_{13}\left({\sigma}^{e}\right)={\sigma}_{1}^{e}-{\sigma}_{3}^{e}+\left({\sigma}_{1}^{e}+{\sigma}_{3}^{e}\right)s-\mathrm{2c}\cos\varphi \ge 0\]

On a l’expression suivante pour la normale:

(3533)\[{n}_{G}=\langle \begin{array}{ccc}t+1& 0& t-1\end{array}\rangle\]

Et donc:

(3534)\[C.{n}_{G}=2\langle \begin{array}{ccc}\left(K+\frac{G}{3}\right)t+G& \left(K-\frac{2}{3}G\right)t& \left(K+\frac{G}{3}\right)t-G\end{array}\rangle\]

Que l’on introduit dans l’équation (), qui devient:

(3535)\[4d\lambda \left[G+(K+\frac{G}{3})ts\right]={\sigma}_{1}^{e}-{\sigma}_{3}^{e}+({\sigma}_{1}^{e}+{\sigma}_{3}^{e})s-\mathrm{2c}\cos\phi\]

D’où on déduit le multiplicateur plastique:

(3536)\[d\lambda =\frac{{\langle {F}_{13}({\sigma}^{e})\rangle }_{+}}{4\left[G+(K+\frac{G}{3})ts\right]}\]

Où \({\langle \rangle }_{+}\) désigne la partie positive d’une grandeur. On obtient finalement:

(3537)\[\begin{split}\left\lbrace \begin{array}{c}{\sigma}_{1}^{+}\\ {\sigma}_{2}^{+}\\ {\sigma}_{3}^{+}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}{\sigma}_{1}^{\text{e}}\\ {\sigma}_{2}^{\text{e}}\\ {\sigma}_{3}^{\text{e}}\end{array}\right\rbrace -\frac{{\langle {F}_{13}\left({\sigma}^{e}\right)\rangle }_{+}}{2G+2\left(K+\frac{G}{3}\right)ts}\left\lbrace \begin{array}{c}\left(K+\frac{G}{3}\right)t+G\\ \left(K-\frac{2}{3}G\right)t\\ \left(K+\frac{G}{3}\right)t-G\end{array}\right\rbrace\end{split}\]

Si on décompose le taux de déformation plastique en une partie déviatorique \(d{\varepsilon}_{v}^{p}=\mathit{trace}(d{\varepsilon}^{p})\) et une partie sphérique \(d{e}^{p}=d{\epsilon}^{p}-\frac{d{\epsilon}_{v}^{p}}{3}1\) , on a:

(3538)\[\begin{split}\lbrace \begin{array}{c}d{\epsilon}_{v}^{p}=\mathrm{2t}d\lambda \\ d{e}^{p}=d\lambda \langle \begin{array}{ccc}\frac{t}{3}+1& \frac{-\mathrm{2t}}{3}& \frac{t}{3}-1\end{array}\rangle \end{array}\end{split}\]

Soit l’expression suivante de l’incrément de déformation équivalente:

(3539)\[d{e}^{p}=d\lambda \sqrt{{t}^{2}+3}\]

Cas où deux mécanismes sont actifs

Formulation de la solution

On vérifie que \({\sigma}^{+}\) obtenue à partir de l’étape précédente () vérifie toujours:

(3540)\[{\sigma}_{1}^{+}\ge {\sigma}_{2}^{+}\ge {\sigma}_{3}^{+}\]

Si ce n’est pas le cas, alors il convient d’activer deux mécanismes de la façon suivante:

(3541)\[\begin{split}\lbrace \begin{array}{cccc}{\sigma}_{2}^{+}\ge {\sigma}_{1}^{+}\ge {\sigma}_{3}^{+}& \Rightarrow & {F}_{13}\mathit{et}{F}_{23}\mathit{actifs}& \textcolor[rgb]{1,0,0}{\left[\text{LEFT}\right]}\\ {\sigma}_{1}^{+}\ge {\sigma}_{3}^{+}\ge {\sigma}_{2}^{+}& \Rightarrow & {F}_{13}\mathit{et}{F}_{12}\mathit{actifs}& \textcolor[rgb]{1,0,0}{\left[\text{RIGHT}\right]}\end{array}\end{split}\]

La définition des mécanismes LEFT et RIGHT est purement conventionnelle, et obéit à la logique géométrique représentée dans la figure .

Figure 3.2.1-1: Définition des mécanismes LEFT et RIGHT

On a les expressions suivantes pour les normales:

(3542)\[{n}_{G}^{13}=\langle \begin{array}{ccc}t+1& 0& t-1\end{array}\rangle {n}_{G}^{12}=\langle \begin{array}{ccc}t+1& t-1& 0\end{array}\rangle\]

Donc:

(3543)\[\begin{split}\lbrace \begin{array}{c}C.{n}_{G}^{13}=2\langle \begin{array}{ccc}\left(K+\frac{G}{3}\right)t+G& \left(K-\frac{2}{3}G\right)t& \left(K+\frac{G}{3}\right)t-G\end{array}\rangle \\ C.{n}_{G}^{12}=2\langle \begin{array}{ccc}\left(K+\frac{G}{3}\right)t+G& \left(K+\frac{G}{3}\right)t-G& \left(K-\frac{2}{3}G\right)t\end{array}\rangle \\ C.{n}_{G}^{23}=2\langle \begin{array}{ccc}\left(K-\frac{2}{3}G\right)t& \left(K+\frac{G}{3}\right)t+G& \left(K+\frac{G}{3}\right)t-G\end{array}\rangle \end{array}\end{split}\]

On introduit ces expressions dans l’équation (). Pour le mécanisme LEFT, on obtient le système suivant à résoudre:

(3544)\[\begin{split}\lbrace \begin{array}{c}4d{\lambda}^{13}\left[G+\left(K+\frac{G}{3}\right)ts\right]+2d{\lambda}^{23}\left[G\left(1-t-s\right)+\left(\mathrm{2K}-\frac{G}{3}\right)ts\right]={F}_{13}\left({\sigma}^{e}\right)\\ 2d{\lambda}^{13}\left[G\left(1-t-s\right)+\left(\mathrm{2K}-\frac{G}{3}\right)ts\right]+4d{\lambda}^{23}\left[G+\left(K+\frac{G}{3}\right)ts\right]={F}_{23}\left({\sigma}^{e}\right)\end{array}\end{split}\]

Pour le mécanisme RIGHT, on obtient:

(3545)\[\begin{split}\lbrace \begin{array}{c}4d{\lambda}^{13}\left[G+\left(K+\frac{G}{3}\right)ts\right]+2d{\lambda}^{12}\left[G\left(1+t+s\right)+\left(\mathrm{2K}-\frac{G}{3}\right)ts\right]={F}_{13}\left({\sigma}^{e}\right)\\ 2d{\lambda}^{13}\left[G\left(1+t+s\right)+\left(\mathrm{2K}-\frac{G}{3}\right)ts\right]+4d{\lambda}^{12}\left[G+\left(K+\frac{G}{3}\right)ts\right]={F}_{12}\left({\sigma}^{e}\right)\end{array}\end{split}\]

Pour la mise en œuvre algorithmique, il est plus agréable de réécrire () et () comme suit:

(3546)\[\begin{split}\lbrace \begin{array}{c}Ad{\lambda}^{1}+{B}^{\mathit{side}}d{\lambda}^{2}={F}_{1}\left({\sigma}^{e}\right)\\ {B}^{\mathit{side}}d{\lambda}^{1}+Ad{\lambda}^{2}={F}_{2}\left({\sigma}^{e}\right)\end{array}\end{split}\]

Avec les multiplicateurs plastiques suivants:

(3547)\[d{\lambda}^{1}=d{\lambda}^{13}\]

Pour les surfaces de charge:

(3548)\[{F}_{1}={F}_{13}\]

L’expression de \(A\) :

(3549)\[A=4\left[G+\left(K+\frac{G}{3}\right)ts\right]\]

Et l’expression de \({B}^{\mathit{side}}\) :

(3550)\[\begin{split}{B}^{\mathit{side}}=\lbrace \begin{array}{c}2\left[G\left(1-t-s\right)+\left(\mathrm{2K}-\frac{G}{3}\right)ts\right]\textcolor[rgb]{1,0,0}{\left[\text{LEFT}\right]}\\ 2\left[G\left(1+t+s\right)+\left(\mathrm{2K}-\frac{G}{3}\right)ts\right]\textcolor[rgb]{1,0,0}{\left[\text{RIGHT}\right]}\end{array}\end{split}\]

La solution du système d’équations () existe si et seulement si son déterminant est non nul, soit:

(3551)\[\begin{split}det=∣\begin{array}{cc}A& {B}^{\mathit{side}}\\ {B}^{\mathit{side}}& A\end{array}∣={A}^{2}-{\left({B}^{\mathit{side}}\right)}^{2}\ne 0\end{split}\]

On peut montrer que cela n’arrive jamais pour des valeurs «physiques» des paramètres \(K\) , \(G\) , \(\varphi\) et \(\psi\) . On a alors comme solutions:

(3552)\[\begin{split}\lbrace \begin{array}{c}d{\lambda}^{1}=\frac{A{F}_{1}\left({\sigma}^{e}\right)-{B}^{\mathit{side}}{F}_{2}\left({\sigma}^{e}\right)}{det}\\ d{\lambda}^{2}=\frac{A{F}_{2}\left({\sigma}^{e}\right)-{B}^{\mathit{side}}{F}_{1}\left({\sigma}^{e}\right)}{det}\end{array}\end{split}\]

Pour le mécanisme LEFT, on obtient finalement:

(3553)\[\begin{split}\left\lbrace \begin{array}{c}{\sigma}_{1}^{+}\\ {\sigma}_{2}^{+}\\ {\sigma}_{3}^{+}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}{\sigma}_{1}^{\text{e}}\\ {\sigma}_{2}^{\text{e}}\\ {\sigma}_{3}^{\text{e}}\end{array}\right\rbrace -2d{\lambda}^{13}\left\lbrace \begin{array}{c}\left(K+\frac{G}{3}\right)t+G\\ \left(K-\frac{2}{3}G\right)t\\ \left(K+\frac{G}{3}\right)t-G\end{array}\right\rbrace -2d{\lambda}^{23}\left\lbrace \begin{array}{c}\left(K-\frac{2}{3}G\right)t\\ \left(K+\frac{G}{3}\right)t+G\\ \left(K+\frac{G}{3}\right)t-G\end{array}\right\rbrace\end{split}\]

Et pour le mécanisme RIGHT:

(3554)\[\begin{split}\left\lbrace \begin{array}{c}{\sigma}_{1}^{+}\\ {\sigma}_{2}^{+}\\ {\sigma}_{3}^{+}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}{\sigma}_{1}^{\text{e}}\\ {\sigma}_{2}^{\text{e}}\\ {\sigma}_{3}^{\text{e}}\end{array}\right\rbrace -2d{\lambda}^{13}\left\lbrace \begin{array}{c}\left(K+\frac{G}{3}\right)t+G\\ \left(K-\frac{2}{3}G\right)t\\ \left(K+\frac{G}{3}\right)t-G\end{array}\right\rbrace -2d{\lambda}^{12}\left\lbrace \begin{array}{c}\left(K+\frac{G}{3}\right)t+G\\ \left(K+\frac{G}{3}\right)t-G\\ \left(K-\frac{2}{3}G\right)t\end{array}\right\rbrace\end{split}\]

Ou encore, en simplifiant l’écriture:

(3555)\[{\sigma}^{+}={\sigma}^{e}-2d{\lambda}^{1}{v}_{t}^{1}-2d{\lambda}^{2}{v}_{t}^{2}\]

Avec les multiplicateurs plastiques suivants:

(3556)\[d{\lambda}^{1}=d{\lambda}^{13}\]

Et les vecteurs:

(3557)\[\begin{split}{v}_{t}^{1}=\left\lbrace \begin{array}{c}\left(K+\frac{G}{3}\right)t+G\\ \left(K-\frac{2}{3}G\right)t\\ \left(K+\frac{G}{3}\right)t-G\end{array}\right\rbrace\end{split}\]

De la même façon, le taux de déformation plastique s’écrit:

(3558)\[\begin{split}\lbrace \begin{array}{c}d{\epsilon}_{v}^{p}=\mathrm{2t}\left(d{\lambda}^{1}+d{\lambda}^{2}\right)\\ d{e}^{p}=d{\lambda}^{1}\left(\begin{array}{ccc}\frac{t}{3}+1& \frac{-\mathrm{2t}}{3}& \frac{t}{3}-1\end{array}\right)+d{\lambda}^{2}\end{array}\lbrace \begin{array}{c}\left(\begin{array}{ccc}\frac{-\mathrm{2t}}{3}& \frac{t}{3}+1& \frac{t}{3}-1\end{array}\right)\textcolor[rgb]{1,0,0}{\left[\text{LEFT}\right]}\\ \left(\begin{array}{ccc}\frac{t}{3}+1& \frac{t}{3}-1& \frac{-\mathrm{2t}}{3}\end{array}\right)\textcolor[rgb]{1,0,0}{\left[\text{RIGHT}\right]}\end{array}\end{split}\]

Soit l’expression suivante de l’incrément de déformation équivalente:

(3559)\[d{e}^{p}=\sqrt{\left({\left(d{\lambda}^{1}\right)}^{2}+{\left(d{\lambda}^{2}\right)}^{2}\right)\left({t}^{2}+3\right)+(-{t}^{2}+\mathit{sign}\times 6t+3)d{\lambda}^{1}d{\lambda}^{2}}\]

Avec \(\mathit{sign}=\lbrace \begin{array}{c}-1\phantom{\rule{2em}{0ex}}\textcolor[rgb]{1,0,0}{\left[\text{LEFT}\right]}\\ +1\phantom{\rule{2em}{0ex}}\textcolor[rgb]{1,0,0}{\left[\text{RIGHT}\right]}\end{array}\) .

Choix du deuxième mécanisme

On propose un critère simple et purement géométrique pour déterminer à coup sûr si le deuxième mécanisme, au cas où il serait activé, se trouve à gauche (LEFT) ou à droite (RIGHT). On définit le vecteur \({t}_{G}\) perpendiculaire à la direction d’écoulement comme représenté dans la figure , de la manière suivante:

(3560)\[\begin{split}\left\lbrace \begin{array}{c}t-1\\ -2\\ t+1\end{array}\right\rbrace ={t}_{G}\perp {n}_{G}={n}_{G}^{13}=\left\lbrace \begin{array}{c}t+1\\ 0\\ t-1\end{array}\right\rbrace\end{split}\]

La droite passant par \(0\) et parallèle à \({n}_{G}\) , représentée en pointillés, caractérise les états de contraintes tels que \(\sigma .{t}_{G}=0\) . On a donc le choix suivant:

• Si \(\sigma .{t}_{G}>0\) et qu’un deuxième mécanisme doit être actif, ce mécanisme est nécessairement situé à droite, soit le mécanisme RIGHT;

• Si \(\sigma .{t}_{G}<0\) et qu’un deuxième mécanisme doit être actif, ce mécanisme est nécessairement situé à gauche, soit le mécanisme LEFT;

Figure 3.2.2-1: Critère de sélection des mécanismes LEFT et RIGHT

Cas de la projection au sommet du cône

On vérifie que \({\sigma}^{+}\) obtenue à partir de l’étape précédente () ou () vérifie toujours:

(3561)\[{\sigma}_{1}^{+}\ge {\sigma}_{2}^{+}\ge {\sigma}_{3}^{+}\]

Si ce n’est pas le cas, alors il convient d’effectuer une projection au sommet du cône d’équation:

(3562)\[{p}_{c}=c\mathrm{cot}\varphi\]

On impose donc:

(3563)\[\begin{split}\lbrace \begin{array}{c}{p}^{+}={p}^{e}-Kd{\epsilon}_{v}^{p}:=c\mathrm{cot}\varphi \\ {\sigma}^{+}:={p}^{+}1\end{array}\end{split}\]

De la même façon, le taux de déformation plastique s’écrit:

(3564)\[d{\epsilon}^{P}=d{\lambda}^{1}{n}_{G}^{13}+d{\lambda}^{2}{n}_{G}^{12}+d{\lambda}^{3}{n}_{G}^{23}\]

Soit:

(3565)\[\begin{split}\lbrace \begin{array}{c}d{\epsilon}_{v}^{p}=2t\left(d{\lambda}^{1}+d{\lambda}^{2}+d{\lambda}^{3}\right)=\frac{1}{K}\left({p}^{\text{e}}-c\mathrm{cot}\varphi \right)\\ d{e}^{p}=d{\lambda}^{1}\left(\begin{array}{ccc}\frac{t}{3}+1& \frac{-\mathrm{2t}}{3}& \frac{t}{3}-1\end{array}\right)+d{\lambda}^{2}\left(\begin{array}{ccc}\frac{t}{3}+1& \frac{t}{3}-1& \frac{-\mathrm{2t}}{3}\end{array}\right)+d{\lambda}^{3}\left(\begin{array}{ccc}\frac{-\mathrm{2t}}{3}& \frac{t}{3}+1& \frac{t}{3}-1\end{array}\right)\end{array}\end{split}\]

Figure 3.3-1: Synoptique de résolution de la loi de Mohr-Coulomb

Expression de la matrice tangente consistante dans la base principale

Cas où un seul mécanisme est actif

La matrice tangente cohérente \(T\) dans la base principale est obtenue en dérivant l’équation (), ce qui donne:

(3566)\[\begin{split}d{\sigma}^{+}=C.d\epsilon -\frac{2d{F}_{13}\left({\sigma}^{e}\right)}{A}\left\lbrace \begin{array}{c}\left(K+\frac{G}{3}\right)t+G\\ \left(K-\frac{2}{3}G\right)t\\ \left(K+\frac{G}{3}\right)t-G\end{array}\right\rbrace =C.d\epsilon -\frac{2d{F}_{13}\left({\sigma}^{e}\right)}{A}{v}_{t}\end{split}\]

Sachant que:

(3567)\[\begin{split}d{F}_{13}\left({\sigma}^{e}\right)=d{\sigma}_{1}^{e}-d{\sigma}_{3}^{e}+\left(d{\sigma}_{1}^{e}+d{\sigma}_{3}^{e}\right)s=2\left\lbrace \begin{array}{c}\left(K+\frac{G}{3}\right)s+G\\ \left(K-\frac{2}{3}G\right)s\\ \left(K+\frac{G}{3}\right)s-G\end{array}\right\rbrace .d\epsilon =2{v}_{s}.d\epsilon\end{split}\]

On obtient:

(3568)\[d{\sigma}^{+}=\underset{T}{\underset{⏟}{\left(C-\frac{4}{A}\underset{D}{\underset{⏟}{{v}_{t}\otimes {v}_{s}}}\right)}}.d\epsilon\]

Où \(A\) est donné par l’équation ().

Cas où deux mécanismes sont actifs

De la même façon, la matrice tangente cohérente \(T\) est obtenue en dérivant les équations () et (), ce qui donne:

(3569)\[d{\sigma}^{+}=\underset{T}{\underset{⏟}{\left(C-\frac{4}{det}\underset{D}{\underset{⏟}{\left({v}_{t}^{1}\otimes {v}_{s}^{1}+{v}_{t}^{2}\otimes {v}_{s}^{2}-\frac{{B}^{\mathit{side}}}{A}\left({v}_{t}^{1}\otimes {v}_{s}^{2}+{v}_{t}^{2}\otimes {v}_{s}^{1}\right)\right)}}\right)}}.d\epsilon\]

Où \(A\) est donné par l’équation (), \({B}^{\mathit{side}}\) par (), \(det\) par () et \(v\) par ().

Cas de la projection au sommet du cône

D’après l’équation (), on a trivialement \(T=0\) .

Expression de la matrice tangente consistante dans la base globale

Le paragraphe § 4 permet de construire la matrice tangente consistante dans la base principale, notée \(T\) . Il convient désormais de ramener cette matrice dans la base globale (cartésienne), que l’on notera \(\stackrel{̄}{T}\) .

Remarque importante:

Il est à noter que la construction de cette matrice tangente consistante est une étape cruciale à la fois pour la robustesse et la performance de l’algorithme:

• Premièrement, il est parfaitement connu qu’une telle matrice permet un taux de convergence quadratique pour le processus de Newton;

• Deuxièmement, cette matrice rend compte de la rotation des directions principales au cours d’un incrément. Sans elle, la formulation de la loi de Mohr-Coulomb en termes de contraintes principales décrite au paragraphe § 2 ne serait pas complète, puisque les contraintes principales, maintenues fixes au cours de l’intégration locale de la loi ( § 3 ), ne pourraient pas tourner au niveau global de la structure.

Dans ce paragraphe, on décrit en détail la méthode permettant de construire \(\stackrel{ˉ}{T}\) à partir de \(T\) .

Quelques résultats sur les tenseurs symétriques isotropes d’ordre deux

On définit par \({S}^{3}\) l’espace des tenseurs symétriques d’ordre deux dans l’espace vectoriel de dimension \(n=3\) , et les tenseurs \(Y\in {S}^{3}\) et \(X\in {S}^{3}\) tels que:

(3570)\[Y(X):{S}^{3}\to {S}^{3}\]

La fonction tensorielle \(Y\left(X\right)\) est dite isotrope si:

(3571)\[R.Y(X).{R}^{t}=Y(R.X.{R}^{t})\]

Quelle que soit la rotation \(R\) . L’hypothèse d’isotropie implique que \(Y\) et \(X\) sont coaxiaux, c’est-à-dire qu’ils possèdent les mêmes directions principales \({d}_{\alpha =1,2,3}\) . On note:

(3572)\[\begin{split}\begin{array}{c}X=\sum_{\alpha =1}^{3}{x}_{\alpha}\underset{{E}_{\alpha}}{\underset{⏟}{\left({d}_{\alpha}\otimes {d}_{\alpha}\right)}}=\sum_{\alpha =1}^{3}{x}_{\alpha}{E}_{\alpha}\\ Y\left(X\right)=\sum_{\alpha =1}^{3}{y}_{\alpha}\left({d}_{\alpha}\otimes {d}_{\alpha}\right)=\sum_{\alpha =1}^{3}{y}_{\alpha}{E}_{\alpha}\end{array}\end{split}\]

Où \({y}_{\alpha}={y}_{\alpha}({x}_{1},{x}_{2},{x}_{3})\) et \({x}_{\alpha}\) représentent les valeurs propres de \(Y\) et \(X\) , respectivement.

Dérivée d’une fonction tensorielle isotrope d’ordre deux

On suppose que la fonction tensorielle isotrope \(Y\left(X\right)\) est différentiable par rapport à \(X\) , et on définit sa dérivée \(D\) telle que:

(3573)\[D\left(X\right)\stackrel{\scriptscriptstyle\mathrm{def}}{=}\frac{dY\left(X\right)}{dX}\]

Appliquée à l’équation (), on obtient l’expression suivante:

(3574)\[D(X)=\sum_{\alpha =1}^{3}({E}_{\alpha}\otimes \frac{d{y}_{\alpha}}{dX}+{y}_{\alpha}\frac{d{E}_{\alpha}}{dX})=\sum_{\alpha =1}^{3}({y}_{\alpha}\frac{d{E}_{\alpha}}{dX}+\sum_{\beta =1}^{3}\frac{\partial {y}_{\alpha}}{\partial {x}_{\beta}}{E}_{\alpha}\otimes \frac{d{x}_{\beta}}{dX})\]

Cas bidimensionnel de type contraintes planes (C_PLAN)

En dimension deux (cas C_PLAN), l’équation caractéristique \(det\left(X-{x}_{\alpha}I\right)=0\) donne une équation quadratique des valeurs propres \({x}_{\alpha =1,2}\) de \(X\) du type suivant:

(3575)\[{x}_{\alpha}^{2}-{I}_{1}{x}_{\alpha}+{I}_{2}=0\]

Avec:

(3576)\[\begin{split}\lbrace \begin{array}{c}{I}_{1}=\mathit{trace}\left(X\right)={X}_{11}+{X}_{22}\\ {I}_{2}=det\left(X\right)={X}_{11}{X}_{22}-{X}_{12}{X}_{21}\end{array}\end{split}\]

La résolution du problème spectral donne aisément les solutions suivantes pour les valeurs propres:

(3577)\[\begin{split}\lbrace \begin{array}{c}{x}_{1}=\frac{{I}_{1}+\sqrt{{I}_{1}^{2}-4{I}_{2}}}{2}\\ {x}_{2}=\frac{{I}_{1}-\sqrt{{I}_{1}^{2}-4{I}_{2}}}{2}\end{array}\end{split}\]

Et les vecteurs propres, tenant compte de la multiplicité des valeurs propres:

(3578)\[\begin{split}\lbrace \begin{array}{cc}{E}_{\alpha}=\frac{X+({x}_{\alpha}-{I}_{1})I}{2{x}_{\alpha}-{I}_{1}}& \mathit{si}{x}_{1}\mathrm{\ne }{x}_{2}\\ {E}_{1}=I& \mathit{si}{x}_{1}={x}_{2}\end{array}\end{split}\]

En particulier, Carlson et Hoger montrent que si \({x}_{1}\mathrm{\ne }{x}_{2}\) , on a:

(3579)\[\frac{d{x}_{\alpha}}{dX}={E}_{\alpha}\]

En utilisant les équations (), () et () dans (), on obtient l’expression de la dérivée \(D(X)\) , tenant compte de la multiplicité des valeurs propres:

(3580)\[\begin{split}D\left(X\right)=\lbrace \begin{array}{cc}\frac{{y}_{1}-{y}_{2}}{{x}_{1}-{x}_{2}}\left[{I}_{S}-{E}_{1}\otimes {E}_{1}-{E}_{2}\otimes {E}_{2}\right]+\sum_{\alpha =1}^{2}\sum_{\beta =1}^{2}\frac{\partial {y}_{\alpha}}{\partial {x}_{\beta}}{E}_{\alpha}\otimes {E}_{\beta}& \mathit{si}{x}_{1}\ne {x}_{2}\\ \left(\frac{\partial {y}_{1}}{\partial {x}_{1}}-\frac{\partial {y}_{1}}{\partial {x}_{2}}\right){I}_{S}+\frac{\partial {y}_{1}}{\partial {x}_{2}}I\otimes I& \mathit{si}{x}_{1}={x}_{2}\end{array}\end{split}\]

Avec la matrice identité \(I\) :

(3581)\[{(I)}_{ijkl}={\delta}_{\mathit{ik}}{\delta}_{\mathit{jl}}\]

La matrice de transposition \({({I}_{t})}_{ijkl}={\delta}_{\mathit{il}}{\delta}_{\mathit{jk}}\) et la matrice de symétrisation \({I}_{S}\) , telles que:

(3582)\[{({I}_{S})}_{ijkl}=\frac{1}{2}(I+{I}_{t})=\frac{1}{2}({\delta}_{\mathit{ik}}{\delta}_{\mathit{jl}}+{\delta}_{\mathit{il}}{\delta}_{\mathit{jk}})\]

Remarque:

On remarque que le terme \(\frac{{y}_{1}-{y}_{2}}{{x}_{1}-{x}_{2}}\left[{I}_{S}-{E}_{1}\otimes {E}_{1}-{E}_{2}\otimes {E}_{2}\right]\) dans la dérivée \(D(X)\) de la première équation de () exprime la rotation des directions principales dans le plan.

Cas bidimensionnel de type déformations planes (D_PLAN) et axisymétrique (AXIS)

La direction hors-plan \(\alpha =3\) étant fixe, l’expression de la dérivée \(D(X)\) s’obtient à partir du cas précédant. En effet, en isolant le terme \(\alpha =3\) dans l’équation (), on a l’expression suivante:

(3583)\[D\left(X\right)=\underset{{D}_{\mathrm{2D}}\left(X\right)}{\underset{⏟}{\sum_{\alpha =1}^{2}\left({y}_{\alpha}\frac{d{E}_{\alpha}}{dX}+\sum_{\beta =1}^{2}\frac{\partial {y}_{\alpha}}{\partial {x}_{\beta}}{E}_{\alpha}\otimes \frac{d{x}_{\beta}}{dX}\right)}}+\underset{{D}_{3}\left(X\right)}{\underset{⏟}{\sum_{\alpha =1}^{2}\frac{\partial {y}_{\alpha}}{\partial {x}_{3}}{E}_{\alpha}\otimes \frac{d{x}_{3}}{dX}+\sum_{\beta =1}^{3}\frac{\partial {y}_{3}}{\partial {x}_{\beta}}{E}_{3}\otimes \frac{d{x}_{\beta}}{dX}}}\]

Où \({D}_{\mathrm{2D}}\left(X\right)\) est donnée par l’équation (). Le terme complémentaire \({D}_{3}\left(X\right)\) s’écrit, en tenant compte de la multiplicité des valeurs propres, de la façon suivante:

(3584)\[\begin{split}{D}_{3}\left(X\right)=\lbrace \begin{array}{cc}\sum_{\alpha =1}^{2}\left(\frac{\partial {y}_{\alpha}}{\partial {x}_{3}}{E}_{\alpha}\otimes {E}_{3}+\frac{\partial {y}_{3}}{\partial {x}_{\alpha}}{E}_{3}\otimes {E}_{\alpha}\right)+\frac{\partial {y}_{3}}{\partial {x}_{3}}{E}_{3}\otimes {E}_{3}& \mathit{si}{x}_{1}\ne {x}_{2}\\ \frac{\partial {y}_{1}}{\partial {x}_{3}}{I}_{p}\otimes {E}_{3}+\frac{\partial {y}_{3}}{\partial {x}_{1}}{E}_{3}\otimes {I}_{p}+\frac{\partial {y}_{3}}{\partial {x}_{3}}{E}_{3}\otimes {E}_{3}& \mathit{si}{x}_{1}={x}_{2}\end{array}\end{split}\]

Où \({I}_{p}\) est la matrice de la projection orthogonale de \({I}_{S}\) dans le plan \(({e}_{x},{e}_{y})\) :

(3585)\[\begin{split}{I}_{p}=\lbrace \begin{array}{cc}\frac{1}{2}\left({\delta}_{\mathit{ik}}{\delta}_{\mathit{jl}}+{\delta}_{\mathit{il}}{\delta}_{\mathit{jk}}\right)& \text{si}i,j,k,l\in \lbrace 1,2\rbrace \\ 0& \text{sinon}\end{array}\end{split}\]

Cas tridimensionnel

En dimension trois, l’équation caractéristique \(det\left(X-{x}_{\alpha}I\right)=0\) donne une équation cubique des valeurs propres \({x}_{\alpha =1,2,3}\) de \(X\) du type suivant:

(3586)\[{x}_{\alpha}^{3}-{I}_{1}{x}_{\alpha}^{2}+{I}_{2}{x}_{\alpha}-{I}_{3}=0\]

Avec:

(3587)\[\begin{split}\lbrace \begin{array}{c}{I}_{1}=\mathit{trace}\left(X\right)\\ {I}_{2}=\frac{1}{2}\left[\mathit{trace}{\left(X\right)}^{2}-\mathit{trace}\left({X}^{2}\right)\right]\\ {I}_{3}=det\left(X\right)\end{array}\end{split}\]

La résolution du problème spectral donne aisément les solutions suivantes pour les valeurs propres:

(3588)\[\begin{split}\lbrace \begin{array}{c}{x}_{1}=-2\sqrt{Q}\cos\left(\frac{\theta}{3}\right)+\frac{{I}_{1}}{3}\\ {x}_{2}=-2\sqrt{Q}\cos\left(\frac{\theta +2\pi }{3}\right)+\frac{{I}_{1}}{3}\\ {x}_{3}=-2\sqrt{Q}\cos\left(\frac{\theta -2\pi }{3}\right)+\frac{{I}_{1}}{3}\end{array}\end{split}\]

Où \(Q\) et \(\theta\) sont donnés par:

(3589)\[Q=\frac{{I}_{1}^{2}-3{I}_{2}}{9}\]

Avec:

(3590)\[R=\frac{-2{I}_{1}^{3}+9{I}_{1}{I}_{2}-27{I}_{3}}{54}\]

Et les vecteurs propres, en tenant compte de la multiplicité des valeurs propres:

(3591)\[\begin{split}\lbrace \begin{array}{cc}{E}_{\alpha}=\frac{{x}_{\alpha}}{2{x}_{\alpha}^{3}-{I}_{1}{x}_{\alpha}^{2}+{I}_{3}}\left[{X}^{2}+\left({x}_{\alpha}-{I}_{1}\right)X+\frac{{I}_{3}}{{x}_{\alpha}}I\right]& \mathit{si}{x}_{1}\ne {x}_{2}\ne {x}_{3}\\ {E}_{\beta}=I-{E}_{\alpha}& \mathit{si}{x}_{\alpha}\ne {x}_{\beta}\\ {E}_{1}=I& \mathit{si}{x}_{1}={x}_{2}={x}_{3}\end{array}\end{split}\]

Dans la deuxième équation de (), \({E}_{\alpha}\) est calculé à l’aide la première équation. Sans donner les étapes intermédiaires de calcul, la dérivée \(D\left(X\right)\) , en tenant compte de la multiplicité des valeurs propres, s’écrit finalement:

(3592)\[\begin{split}D\left(X\right)=\lbrace \begin{array}{cc}\begin{array}{c}\sum_{\alpha =1}^{3}\frac{{y}_{\alpha}}{\left({x}_{\alpha}-{x}_{\beta}\right)\left({x}_{\alpha}-{x}_{\gamma}\right)}[\frac{d{X}^{2}}{dX}-\left({x}_{\beta}+{x}_{\gamma}\right){I}_{S}\\ -\left(2{x}_{\alpha}-{x}_{\beta}-{x}_{\gamma}\right){E}_{\alpha}\otimes {E}_{\alpha}-\left({x}_{\beta}-{x}_{\gamma}\right)\left({E}_{\beta}\otimes {E}_{\beta}-{E}_{\gamma}\otimes {E}_{\gamma}\right)]\\ +\sum_{a=1}^{3}\sum_{b=1}^{3}\frac{\partial {y}_{a}}{\partial {x}_{b}}{E}_{a}\otimes {E}_{b}\end{array}& \mathit{si}{x}_{1}\ne {x}_{2}\ne {x}_{3}\\ {s}_{1}\frac{d{X}^{2}}{dX}-{s}_{2}{I}_{S}-{s}_{3}X\otimes X+{s}_{4}X\otimes I+{s}_{5}I\otimes X-{s}_{6}I\otimes I& \mathit{si}{x}_{\alpha}\ne {x}_{\beta}={x}_{\gamma}\\ \left(\frac{\partial {y}_{1}}{\partial {x}_{1}}-\frac{\partial {y}_{1}}{\partial {x}_{2}}\right){I}_{S}+\frac{\partial {y}_{1}}{\partial {x}_{2}}I\otimes I& \mathit{si}{x}_{1}={x}_{2}={x}_{3}\end{array}\end{split}\]

Où \(\left(\alpha ,\beta ,\gamma \right)\) correspond à une permutation cyclique de \(\left(1,2,3\right)\) . \(I\) et \({I}_{S}\) sont donnés par les équations () et (), respectivement. En remarquant que \(X\) est un tenseur symétrique, il faut prendre garde à appliquer l’opérateur de dérivation symétrique pour l’évaluation de \(\frac{d{X}^{2}}{dX}\) , ce qui donne la forme suivante:

(3593)\[\begin{split}\begin{array}{cc}{\left(\frac{d{X}^{2}}{dX}\right)}_{ijkl}& =\frac{d\left({X}_{im}{X}_{\mathit{mj}}\right)}{d{X}_{kl}}=\frac{1}{2}\left({\delta}_{\mathit{ik}}{\delta}_{\mathit{lm}}+{\delta}_{\mathit{il}}{\delta}_{\mathit{km}}\right){X}_{\mathit{mj}}+\frac{{X}_{im}}{2}\left({\delta}_{\mathit{mk}}{\delta}_{\mathit{jl}}+{\delta}_{\mathit{ml}}{\delta}_{\mathit{kj}}\right)\\ & =\frac{1}{2}\left({\delta}_{\mathit{ik}}{X}_{\mathit{lj}}+{\delta}_{\mathit{il}}{X}_{\mathit{kj}}+{\delta}_{\mathit{jl}}{X}_{\mathit{ik}}+{\delta}_{\mathit{kj}}{X}_{\mathit{il}}\right)\end{array}\end{split}\]

Enfin, les expressions de \({s}_{i=1,6}\) sont les suivantes:

(3594)\[{s}_{1}=\frac{{y}_{\alpha}-{y}_{\gamma}}{{\left({x}_{\alpha}-{x}_{\gamma}\right)}^{2}}+\frac{1}{{x}_{\alpha}-{x}_{\gamma}}\left(\frac{\partial {y}_{\gamma}}{\partial {x}_{\beta}}-\frac{\partial {y}_{\gamma}}{\partial {x}_{\gamma}}\right) {s}_{2}=2{x}_{\gamma}\frac{{y}_{\alpha}-{y}_{\gamma}}{{\left({x}_{\alpha}-{x}_{\gamma}\right)}^{2}}+\frac{{x}_{\alpha}+{x}_{\gamma}}{{x}_{\alpha}-{x}_{\gamma}}\left(\frac{\partial {y}_{\gamma}}{\partial {x}_{\beta}}-\frac{\partial {y}_{\gamma}}{\partial {x}_{\gamma}}\right) {s}_{3}=2\frac{{y}_{\alpha}-{y}_{\gamma}}{{\left({x}_{\alpha}-{x}_{\gamma}\right)}^{3}}+\frac{1}{{\left({x}_{\alpha}-{x}_{\gamma}\right)}^{2}}\left(\frac{\partial {y}_{\alpha}}{\partial {x}_{\gamma}}+\frac{\partial {y}_{\gamma}}{\partial {x}_{\alpha}}-\frac{\partial {y}_{\alpha}}{\partial {x}_{\alpha}}-\frac{\partial {y}_{\gamma}}{\partial {x}_{\gamma}}\right) {s}_{4}=2{x}_{\gamma}\frac{{y}_{\alpha}-{y}_{\gamma}}{{\left({x}_{\alpha}-{x}_{\gamma}\right)}^{3}}+\frac{1}{{x}_{\alpha}-{x}_{\gamma}}\left(\frac{\partial {y}_{\alpha}}{\partial {x}_{\gamma}}-\frac{\partial {y}_{\gamma}}{\partial {x}_{\beta}}\right)+\frac{{x}_{\gamma}}{{\left({x}_{\alpha}-{x}_{\gamma}\right)}^{2}}\left(\frac{\partial {y}_{\alpha}}{\partial {x}_{\gamma}}+\frac{\partial {y}_{\gamma}}{\partial {x}_{\alpha}}-\frac{\partial {y}_{\alpha}}{\partial {x}_{\alpha}}-\frac{\partial {y}_{\gamma}}{\partial {x}_{\gamma}}\right)\]

Où \((\alpha ,\beta ,\gamma )\) correspond à une permutation cyclique de \((1,2,3)\) .

Remarque:

On remarque que le terme suivant (tout sauf trivial):

(3595)\[\begin{split}\begin{array}{c}\sum_{\alpha =1}^{3}\frac{{y}_{\alpha}}{\left({x}_{\alpha}-{x}_{\beta}\right)\left({x}_{\alpha}-{x}_{\gamma}\right)}[\frac{d{X}^{2}}{dX}-\left({x}_{\beta}+{x}_{\gamma}\right){I}_{S}-\left(2{x}_{\alpha}-{x}_{\beta}-{x}_{\gamma}\right){E}_{\alpha}\otimes {E}_{\alpha}\\ -\left({x}_{\beta}-{x}_{\gamma}\right)\left({E}_{\beta}\otimes {E}_{\beta}-{E}_{\gamma}\otimes {E}_{\gamma}\right)]\end{array}\end{split}\]

Qui apparaît dans la dérivée \(D(X)\) de la première équation de () exprime la rotation des directions principales dans l’espace tridimensionnelle.

Application au cas de Mohr-Coulomb

La transposition des formules précédentes à la mise en œuvre numérique mérite quelques précisions. On a tout d’abord les correspondances suivantes:

• \(X={\stackrel{̃}{\epsilon}}^{\mathit{pred}}\) et \({x}_{\alpha}={\epsilon}_{\alpha}^{\mathit{pred}}\) ;

• \(Y={\stackrel{̃}{\sigma}}^{+}\) et \({y}_{\alpha}={\sigma}_{\alpha}^{+}\) ;

• \({E}_{\alpha}={\stackrel{̃}{d}}_{\alpha}^{\mathit{pred}}\otimes {\stackrel{̃}{d}}_{\alpha}^{\mathit{pred}}\) ;

• \({\left(T\right)}_{\alpha \beta }=\frac{\partial {y}_{\alpha}}{\partial {x}_{\beta}}\) est la matrice tangente consistante dans la base principale calculée au paragraphe § 4 ;

La notation \({}^{\mathit{pred}}\) indique que l’on travaille avec des grandeurs «prédites» données en entrée par le processus de Newton, la notation \({}^{+}\) avec des grandeurs issues de la résolution locale de la loi de comportement, et la notation \(\tilde{\rbrace\) avec la base de Voigt. On notera que les directions principales prédites \({\tilde{d}}_{\alpha}^{\mathit{pred}}\) sont figées au cours de la résolution locale, ce qui est cohérent avec l’hypothèse d’isotropie adoptée (voir les explications du paragraphe § 5.1 ).

Disposant de toutes ces informations à l’issue de la résolution locale de la loi de comportement, on en déduit la matrice tangente consistante \(\stackrel{ˉ}{T}=\stackrel{ˉ}{D}\) exprimée dans la base de projection \(\stackrel{ˉ}{b}\) définie au paragraphe § 1.2 :

• L’équation () dans le cas bidimensionnel en contraintes planes (C_PLAN);

• L’équation () dans le cas bidimensionnel en déformation plane (D_PLAN) ou axisymétrique (AXIS);

• L’équation () dans le cas tridimensionnel (3D);

La seconde information importante concerne la convention d’écriture des différents tenseurs. En effet, par souci de généralité, la notation utilisée pour les tenseurs dans tout le paragraphe § 5 est la notation classique, faisant apparaître des tenseurs jusqu’à l’ordre quatre. Cette écriture est impropre à la résolution numérique, où l’on préfère utiliser des notations condensées rendues possibles par le fait que l’on travaille avec des tenseurs symétriques d’ordre deux (contraintes et déformations le sont toujours). On distingue deux formes de notations condensées correspondant à deux bases de projection (voir § 1.2 ):

• La base orthonormée \(\stackrel{ˉ}{b}\) des tenseurs symétriques d’ordre deux. C’est dans cette base que sont données les contraintes et les déformations à l’entrée et à la sortie de la résolution locale du comportement;

• La base dite de Voigt \(\tilde{b}\) , beaucoup plus commode à utiliser lors de la résolution numérique locale du comportement car elle évite d’avoir à manipuler des coefficients en \(\sqrt{2}\) lors des opérations matricielles;

Le schéma de résolution est résumé dans la figure .

Figure 5.2.4-1: Processus de résolution de la loi de Mohr-Coulomb. Mise en évidence de l’écriture dans les différentes bases de projection.

Bibliographie

[78] Jiang Hua , Xie Y. « A note on the Mohr-Coulomb and Drucker-Prager strenght criteria », Mechanics Research Communications, 38 , pp, 309-314, 2011

[79] Fotios E. Karaoulanis « Nonsmooth multisurface plasticity in principal stress space », 6th GRACM International Congress on Computational Mechanics, Thessaloniki, 2008

[80] Pankaj , Bićanić N. « Detection of multiple active yield conditions for Mohr-Coulomb elasto- plasticity », Computers and Structures, 62 , pp, 51-61, 1997

[81] S. W. Sloan , Booker J.R. « Removal of singularities in Tresca and Mohr-Coulomb yield functions », Communications in Applied Numerical Methods, 2 , pp, 173-179, 1986

[82] Ronaldo I. Borja , Sama K.M., Sanz P.F. « On the numerical integration of three-invariant elastoplastic constitutive models », Comput. Methods Appl. Mech. Engrg., 192 , pp, 1227- 1258, 2003

[83] X. Wang , Wang L.B., Xu L.M. « Formulation of the return mapping algorithm for elastoplastic soil models », Computers and Geotechnics, 31 , pp, 315-338, 2004