r7.01.16 Intégration des comportements mécaniques élasto-plastiques de Drucker-Prager, associé (DRUCK_PRAGER)et non-associé (DRUCK_PRAG_N_A) et post-traitements#
Résumé:
Ce document décrit les principes de plusieurs développements concernant la loi de comportement élastoplastique de Drucker-Prager en version associée (DRUCK_PRAGER)et non-associée (DRUCK_PRAG_N_A).
On s’intéresse d’abord à l’intégration proprement dite de la loi puis, cette loi étant adoucissante, à un indicateur de localisation de Rice et enfin au calcul de sensibilité par différentiation directe pour cette loi.
Pour l’intégration de la loi, on utilise un schéma implicite.
Table des matières
Intégration de la loi de comportement de Drucker-Prager#
Notations#
Les contraintes mécaniques sont comptées positives en traction, les déformations positives en extension.
\(u\) |
déplacements du squelette de composantes \({u}_{x},{u}_{y},{u}_{z}\) |
\(\varepsilon =\frac{1}{2}(\nabla u+{\nabla}^{T}u)\) |
tenseur des déformations linéarisées |
\(e=\varepsilon -\frac{\text{Tr}(\varepsilon )}{3}I\) |
déviateur des déformations |
\({\varepsilon}_{v}=\text{Tr}(\varepsilon )\) |
trace des déformations: variation de volume |
\({\varepsilon}^{p}\) |
Tenseur des déformations plastiques, |
\({\varepsilon}_{v}^{p}=\text{Tr}({\varepsilon}^{p})\) |
variation de volume plastique. |
\({e}^{p}\) |
déviateur des déformations plastiques |
\(p\) |
déformation plastique cumulée |
\(\sigma\) |
Tenseur des contraintes |
\(s=\sigma -\frac{\text{Tr}\left(\sigma \right)}{3}I\) |
déviateur des contraintes |
\({\sigma}_{\text{eq}}=\sqrt{\frac{3}{2}s:s}\) |
Contrainte équivalente de Von Mises |
\({I}_{1}=\text{Tr}\left(\sigma \right)\) |
trace des contraintes |
\({E}_{0}\) |
Module de Young |
\({\nu}_{0}\) |
coefficient de Poisson |
\(\varphi\) |
Angle de frottement |
\(c\) |
Cohésion |
\({\psi}^{0}\) |
Angle de dilatance initial |
On pose \(2\mu =\frac{{E}_{0}}{1+{\nu}_{0}}\) et \(K=\frac{{E}_{0}}{3(1-2{\nu}_{0})}\) .
Formulation en version associée#
Expression du comportement#
\(\sigma\) est le tenseur des contraintes, qui ne dépend que de \(\varepsilon\) et son histoire. On considère le critère de type Drücker-Prager :
où \(A\) est un coefficient donné et \(R\) est une fonction de la déformation plastique cumulée \(p\) (fonction d’écrouissage), de type linéaire ou parabolique :
écrouissage linéaire
Les coefficients \(h\) , \({p}_{\mathit{ultm}}\) et \({\sigma}_{Y}\) sont donnés.
écrouissage parabolique
\(R(p)={\sigma}_{Y}(1-(1-\sqrt{\frac{{\sigma}_{\text{Yultm}}}{{\sigma}_{Y}}})\frac{p}{{p}_{\text{ultm}}}{)}^{2}\) |
si \(p\in [0,{p}_{\mathit{ultm}}]\) |
\(R(p)={\sigma}_{{Y}_{\text{ultm}}}\) |
si \(p>{p}_{\mathit{ultm}}\) |
Les coefficients \({\sigma}_{\mathit{Yultm}}\) , \({p}_{\mathit{ultm}}\) et \({\sigma}_{Y}\) sont donnés.
Remarque 1 :
On peut se donner à la place de \(A\) et \({\sigma}_{Y}\) le coefficient de cohésion \(c\) et l’angle de frottement \(\varphi\) :
\(A=\frac{2\sin\varphi }{3-\sin\varphi }\)
\({\sigma}_{Y}=\frac{\mathrm{6c}\cos\varphi }{3-\sin\varphi }\)
Remarque 2 :
On a choisi dans ce document de privilégier la variable \(p\) *. La déformation plastique cumulée de cisaillement* \({\gamma}^{p}=p\sqrt{3/2}\) est également très utilisée en mécanique des sols.
En considérant une version associée on suppose que le potentiel de dissipation suit la même expression que celle de la surface de charge \(F\) . L’écoulement plastique se résume alors à :
avec :
La loi de normalité par rapport à la force généralisée \(R\) donne l’égalité entre l’incrément de déformation plastique cumulée et l’incrément du multiplicateur \(\lambda\) :
Résolution analytique de la formulation mécanique#
On se place dans ce chapitre dans le cadre d’accroissement fini. L’intégration de la loi suit un schéma implicite pur, et la résolution est analytique. L’incrément fini de déformation \(\Delta \epsilon\) est connu et fourni par l’itération de Newton global. On utilise par convention les notations suivantes : un indice – pour désigner une composante en début de pas de chargement, aucun indice pour une composante en fin de pas de chargement, et l’opérateur \(\Delta\) pour désigner l’accroissement d’une composante. Les équations traduisant le comportement élastique s’écrivent alors :
Les équations (4896) et (4898), compte tenu de la (4894), donnent :
D’où :
Au cas où l’incrément \(\Delta {e}^{p}\) soit non nul, l’incrément de déformation plastique cumulée peut s’écrire aussi :
En combinant les équations (4899) et (3684) on trouve :
d’où :
ce qui conduit à :
En combinant respectivement les équations (3686) et (3687), et les équations (3681) et (3683), on obtient :
En réinjectant l’équation sur \({I}_{1}\) et la relation \({\sigma}_{\text{eq}}={\sigma}_{\text{eq}}^{e}-3\mu \cdot \Delta p\) dans la formulation du seuil, on obtient l’équation scalaire en \(\Delta p\) :
On suppose que : \(F({\sigma}^{e},{p}^{-})>0\) .
Pour poursuivre la résolution, on doit maintenant distinguer plusieurs cas :
Cas où \({p}^{-}>{p}_{\mathit{ultm}}\)
On a : \(R({p}^{-}+\Delta p)=R({p}^{-})\)
L’équation scalaire devient donc : \(F({\sigma}^{e},{p}^{-})-\Delta p(\mathrm{3\mu }+9{\text{KA}}^{2})=0\)
On trouve :
Cas où \({p}^{-}\le {p}_{\mathrm{ultm}}\)
2a) Écrouissage linéaire
On a : \(R({p}^{-}+\Delta p)=R({p}^{-})+h\Delta p\)
L’équation scalaire devient donc : \(F\left({\sigma}^{e},{p}^{-}\right)-\Delta p\left(3\mu +9{\text{KA}}^{2}+h\right)=0\)
On trouve :
2b) Écrouissage parabolique
En exprimant de la même façon \(R({p}^{-}+\Delta p)\) en fonction de \(R({p}^{-})\) et de \(\Delta p\) , on trouve que l’équation scalaire s’écrit :
\(F\left({\sigma}^{e},{p}^{-}\right)+B\Delta p+G\Delta {p}^{2}=0\)
avec :
\(\left\{ \begin{array}{}G=-\frac{{\sigma}_{Y}}{{p}_{\text{ultm}}^{2}}{(1-\sqrt{\frac{{\sigma}_{\text{Yultm}}}{{\sigma}_{Y}}})}^{2}\\ B=-3\mu -9{\text{KA}}^{2}+\frac{2{\sigma}_{Y}}{{p}_{\text{ultm}}}(1-\sqrt{\frac{{\sigma}_{\text{Yultm}}}{{\sigma}_{Y}}})(1-(1-\sqrt{\frac{{\sigma}_{\text{Yultm}}}{{\sigma}_{Y}}})\frac{{p}^{-}}{{p}_{\text{ultm}}})\end{array} \right.\)
La seule racine positive du polynôme est :
2c) Vérification finale : Cas où \(({p}^{-}+\Delta p)>{p}_{\mathrm{ultm}}\)
Dans les deux cas précédents, une fois \(\Delta p\) calculé, il faut vérifier que \({p}^{-}+\Delta p\le {p}_{\mathit{ultm}}\) . Si cette inégalité n’est pas satisfaite, on a alors :
\(R({p}^{-}+\Delta p)=R({p}_{\mathrm{ultm}})\)
L’équation scalaire devient donc :
\(F({\sigma}^{e},{p}_{\text{ultm}})-\Delta p(3\mu +9{\text{KA}}^{2})=0\)
Le principe de la résolution analytique présentée ci-dessus est équivalente à déterminer le point \(({I}_{1},s)\) comme la projection du point \(({I}_{{}_{1}}^{e},{s}^{e})\) sur le critère (prédiction élastique-correction plastique). Cette méthode vient donc de la loi d’écoulement approximée sur un incrément fini, et peut être représentée par le graphique suivant :
Fig. 207 Projection sur le critère#
Projection au sommet du cône
L’intégration de la loi sur un incrément \(\Delta t\) fini peut être compliqué quand l’état de contrainte est proche du sommet du cône (voir la Figure Fig. 207), à cause du caractère non lisse de la surface critère. Il y a alors deux cas :
cas d’un état hydrostatique pur,
cas de projection dans un domaine non-admissible.
Dans le cas particulier d’un état hydrostatique pur , la dérivée de la contrainte de von Mises \({\sigma}_{\text{eq}}\) par rapport à \(\sigma\) n’est pas définie. La loi d’écoulement (3682) est indéterminée (il y a en effet un cône de normales possibles au critère), et les équations (3684), (3686), (3687), (3688) ne peuvent pas être écrites. Il reste la définition de \(\Delta p\) sur la trace des contraintes (équation (3683)). Comme dans le cas plus général, on doit distinguer plusieurs cas :
Cas où \({p}^{-}>{p}_{\mathit{ultm}}\) : \(R({p}^{-}+\Delta p)=R({p}^{-})\)
L’équation scalaire avec \({\sigma}_{\mathrm{eq}}=0\) devient : \(A{I}_{1}^{e}-\Delta p\cdot 9{\mathrm{KA}}^{2}=F({\sigma}^{e},{p}^{-})-\Delta p\cdot 9{\mathrm{KA}}^{2}=0\)
On trouve :
Cas où \({p}^{-}\le {p}_{\mathrm{ultm}}\)
2a) Écrouissage linéaire
On a : \(R({p}^{-}+\Delta p)=R({p}^{-})+h\Delta p\)
L’équation scalaire avec \({\sigma}_{\mathrm{eq}}=0\) devient :
\(A{I}_{1}^{e}-\Delta p\cdot 9{\mathrm{KA}}^{2}-R({p}^{-})+h\Delta p=F({\sigma}^{e},{p}^{-})-\Delta p\cdot 9{\mathrm{KA}}^{2}=0\)
On trouve alors :
2b) Écrouissage parabolique
En exprimant \(R({p}^{-}+\Delta p)\) en fonction de \(R({p}^{-})\) et de \(\Delta p\) , on trouve encore la solution (3693) :
avec la valeur de \(B\) modifiée par rapport au cas précédent :
\(B=-9{\text{KA}}^{2}+\frac{2{\sigma}_{Y}}{{p}_{\text{ultm}}}(1-\sqrt{\frac{{\sigma}_{\text{Yultm}}}{{\sigma}_{Y}}})(1-(1-\sqrt{\frac{{\sigma}_{\text{Yultm}}}{{\sigma}_{Y}}})\frac{{p}^{-}}{{p}_{\text{ultm}}})\)
2c) Vérification finale : Cas où \(({p}^{-}+\Delta p)>{p}_{\mathrm{ultm}}\)
Dans les cas 2a) et 2b), si l’inégalité \({p}^{-}+\Delta p\le {p}_{\mathit{ultm}}\) n’est pas satisfaite, on a :
\(R({p}^{-}+\Delta p)=R({p}_{\mathrm{ultm}})\)
L’incrément \(\Delta p\) est donné par l’équation (3695).
A cause de la résolution incrémentale, il se peut que la solution trouvée ne soit pas admissible, avec \({\sigma}_{\text{eq}}<0\) . Cela peut arriver quand l’état de contrainte à l’instant \({t}^{-}\) est proche du sommet du cône.
On choisit alors de projeter l’état de contrainte trouvé par prédiction élastique sur le sommet du cône, soit de se reporter à un état de contrainte purement hydrostatique. On fait un contrôle à posteriori de l’admissibilité de la solution \(\left({I}_{1},\mathrm{s}\right)\) , et on apporte éventuellement la correction.
Dans les détails :
On actualise l’état de contrainte par le biais des équations (3691), (3692), (3693), (3694).
On contrôle que la solution \(\left({I}_{1},\mathrm{s}\right)\) trouvée soit admissible, soit que \({\sigma}_{\text{eq}}<0\) où, de façon équivalente, que \({I}_{1}\) soit à l’intérieur de la surface critère :
\({I}_{1}\le \frac{R(p)}{A}\)
Calcul de l’opérateur tangent#
Calcul global de l’opérateur tangent#
On cherche à calculer la matrice cohérente : \(\frac{\partial \sigma }{\partial \varepsilon }=\frac{\partial \mathrm{s}}{\partial \varepsilon }+\frac{1}{3}\mathrm{I}\otimes \frac{\partial {I}_{1}}{\partial \varepsilon }\)
En dérivant le système d’équations (3686), on obtient :
Expression de \(\frac{\partial {\mathrm{s}}^{\mathrm{e}}}{\partial \varepsilon }\)
\(\frac{\partial {s}_{ij}^{e}}{\partial {\varepsilon}_{\text{pq}}}=2\mu ({\delta}_{\text{ip}}{\delta}_{\text{jq}}-\frac{1}{3}{\delta}_{ij}{\delta}_{\text{pq}})\)
Expression de \(\frac{\partial {I}_{1}^{e}}{\partial \varepsilon }\)
\(\frac{\partial {I}_{1}^{e}}{\partial {\varepsilon}_{\text{pq}}}=\mathrm{3K}{\delta}_{\text{pq}}\)
Calcul de \(\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \varepsilon }\)
\(\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial {\varepsilon}_{\text{pq}}}=\frac{3\mu }{{\sigma}_{\text{eq}}^{e}}{s}_{\text{pq}}^{e}\)
Calcul de \(\frac{\partial \Delta p}{\partial \varepsilon }\)
\(\frac{\partial \Delta p}{\partial {\varepsilon}_{\text{pq}}}=-\frac{1}{T(\Delta p)}.(\frac{3\mu }{{\sigma}_{\text{eq}}^{e}}{s}_{\text{pq}}^{e}+3\text{AK}{\delta}_{\text{pq}})\)
avec:
\(T(\Delta p)=\lbrace \begin{array}{c}\begin{array}{}-(3\mu +9{\text{KA}}^{2})\text{}\text{dans le cas}{p}^{-}+\Delta p\ge {p}_{\text{ultm}}\text{}(\text{écrouissage linéaire ou parabolique})\\ -(3\mu +9{\text{KA}}^{2}+h)\text{}\text{dans le cas}{p}^{-}+\Delta p<{p}_{\text{ultm}}\text{}(\text{écrouissage linéaire})\end{array}\\ B+\text{2G}\Delta \text{p}\text{dans le cas}{p}^{-}+\Delta p<{p}_{\text{ultm}}\text{}(\text{écrouissage parabolique})\text{}\end{array}\)
où \(B\) et \(G\) ont la même expression qu’au paragraphe [r7.01.16-formulation-version-associee].
Expression complète
Calcul initial de l’opérateur tangent#
On cherche a exprimer \(\frac{\partial {\sigma}^{-}}{\partial {\varepsilon}^{-}}\) . Pour cela on va chercher à calculer l’opérateur tangent par un calcul en vitesse: \(\frac{\partial \dot{\sigma}}{\partial \dot{\varepsilon}}\) .
En partant de l’expression: \(\stackrel{}{F}=\frac{\partial F}{\partial \stackrel{}{\sigma}}+\frac{\partial F}{\partial p}\stackrel{}{p}=0\) on démontre que:
\(\dot{p}=\frac{\mathrm{3\mu }}{{\sigma}_{\text{eq}}D}s.\dot{\epsilon}+\frac{3\text{AK}}{D}\dot{{\epsilon}_{v}}\) avec \(D=\mathrm{3\mu }+9{\text{KA}}^{2}+\frac{\partial R}{\partial p}\)
Des expressions: \(\dot{\sigma}=\mathrm{H}\left(\dot{\epsilon}-{\dot{\epsilon}}^{p}\right)\) et \({\stackrel{}{\epsilon}}^{p}=\stackrel{}{p}\frac{\partial F}{\partial \sigma }\) on démontre ensuite que:
\(\frac{\partial \dot{\sigma}}{\partial \dot{\varepsilon}}=H-(\frac{3\mu }{{\sigma}_{\mathrm{eq}}}s+3\mathrm{AK}I)\frac{\partial \Delta p}{\partial \epsilon }\)
qui n’est autre que l’expression de la matrice cohérente du système global du paragraphe précédent où \(\mathrm{\Delta p}=0\) .
Formulation en version non-associée#
La version non-associée de la loi Drucker-Prager n’a pas pour prétention de modéliser finement un comportement physique réaliste. Le but est de représenter le plus simplement possible une physique grossièrement réaliste, notamment dans le cas de la mécanique des sols pour laquelle l’angle de dilatance varie avec la déformation plastique.
Dans cette formulation, le critère s’écrit de la même manière que la version associée
\(F(\sigma ,p)={\sigma}_{\text{eq}}+{\text{AI}}_{1}-R(p)\le 0\)
La fonction \(R(p)\) peut prendre plusieurs formes, en fonction du type d’écrouissage (linéaire, parabolique, exponentiel):
écrouissage linéaire (identique à la version associée de la loi de comportement)
\(R(p)={\sigma}_{Y}+h.p\) |
si \(p\in [0,{p}_{\mathit{ultm}}]\) |
\(R(p)={\sigma}_{Y}+h.{p}_{\text{ultm}}\) |
si \(p>{p}_{\mathit{ultm}}\) |
Les coefficients \(h\) , \({p}_{\mathit{ultm}}\) et \({\sigma}_{Y}\) sont donnés.
écrouissage parabolique (différent de celui de la loi associée)
\(R(p)=({\sigma}_{Y}-{\sigma}_{\mathit{ult}}){(\frac{p}{{p}_{\mathit{ult}}})}^{2}+2({\sigma}_{Y}-{\sigma}_{\mathit{ult}})\) |
si \(p\in [0,{p}_{\mathit{ultm}}]\) |
\(R(p)={\sigma}_{\mathit{ult}}\) |
si \(p>{p}_{\mathit{ultm}}\) |
Les coefficients \({\sigma}_{\mathit{ult}}\) , \({p}_{\mathit{ultm}}\) et \({\sigma}_{Y}\) sont donnés.
écrouissage exponentiel (spécifique à la loi non associée)
\(R(p)=({\sigma}_{Y}-{\sigma}_{\mathit{ult}}){(\frac{p}{{p}_{\mathit{ult}}})}^{2}+2({\sigma}_{Y}-{\sigma}_{\mathit{ult}})\) |
si \(p\in [0,{p}_{\mathit{ultm}}]\) |
\(R(p)={\sigma}_{\mathit{ult}}\) |
si \(p>{p}_{\mathit{ultm}}\) |
Les coefficients \({\sigma}_{\mathit{ult}}\) , \({p}_{\mathit{ultm}}\) et \({\sigma}_{Y}\) sont donnés.
Le potentiel plastique se différencie donc de la surface de charge dans cette nouvelle formulation. Le potentiel plastique est le suivant: \(G(\sigma ,p)={\sigma}_{\text{eq}}+\beta (p){I}_{1}\)
où \(\beta \left(p\right)\) est une fonction qui décroît avec l’évolution de la déformation plastique. Dans le cas d’un écrouissage linéaire ou parabolique, cette fonction \(\beta \left(p\right)\) est linéaire. Elle s’écrit:
\(\beta (p)=\left\{ \begin{array}{cc}\text{}\beta ({\psi}^{0})(1-\frac{p}{{p}_{\text{ult}}})& \text{si}p\in \text{}\left[0,{p}_{\text{ult}}\right]\\ \begin{array}{}\\ 0\end{array}& \begin{array}{}\\ \text{si}p>{p}_{\text{ult}}\end{array}\end{array} \right.\)
où \({\psi}^{0}\) désigne l’angle de dilatance initial et \(\beta ({\psi}^{0})=\frac{2\sin({\psi}^{0})}{3-\sin({\psi}^{0})}\) .
Dans le cas d’un écrouissage exponentiel, cette fonction \(\beta (p)\) s’écrit:
\(\beta (p)=({\beta}_{0}-{\beta}_{\mathit{ult}})\exp(\frac{-p}{{p}_{c}})+{\beta}_{\mathit{ult}}\)
Résolution analytique#
La méthode d’intégration de la loi de comportement dans sa version non associée diffère de celle utilisée pour la version associée.
Interprétation de la loi d’écoulement#
La loi d’écoulement s’écrit:
\(\dot{\epsilon}=\lambda M(\sigma )\mathit{avec}\lambda \ge 0,F\ge 0\mathit{et}\lambda F=0\)
\(\lambda ` désigne ici le multiplicateur plastique, et :math:`M(\sigma )\) la direction d’écoulement. Cette direction d’écoulement appartient au sous-différentiel de \(\sigma ` , appelé :math:`N(\sigma )\) . Ce sous différentiel est défini de la manière suivante:
\(N(\sigma )=\lbrace M;\forall \sigma '\mathit{tel}G(\sigma ')\le G(\sigma ),(\sigma '-\sigma ):M\le 0\rbrace\)
Dans le cas où le potentiel d’écoulement est dérivable en \(\sigma ` , c’est à dire pour :math:`{\sigma}_{\mathit{eq}}\ne 0\) , la direction d’écoulement s’écrit:
\(M=\nabla G=\frac{3}{2}\frac{{\sigma}^{D}}{{\sigma}_{\mathit{eq}}}+\beta (p)I\)
Dans ce cas, on peut aussi montrer que \(\lambda =\dot{p}\) :
\(\dot{{\epsilon}_{p}}=\lambda \left(\frac{3}{2}\frac{{\sigma}^{D}}{{\sigma}_{\mathit{eq}}}+\beta (p)I\right)\)
\(\dot{{\epsilon}_{p}^{D}}=\lambda \frac{3}{2}\frac{{\sigma}^{D}}{{\sigma}_{\mathit{eq}}}\)
\(\dot{{\epsilon}_{p}^{D}}:\dot{{\epsilon}_{p}^{D}}=\frac{3}{2}{\lambda}^{2}\)
\(\lambda =\dot{p}\)
Dans le cas où le potentiel d’écoulement n’est pas dérivable en \(\sigma ` , c’est à dire pour :math:`{\sigma}_{\mathit{eq}}=0\) ,
:math:`` si et seulement si:
\(max(\sigma '-\sigma ):M\le 0\)
De là, en utilisant le fait que \(G(\sigma ')\le G(\sigma )\) , on aboutit à:
\(2\beta (p){M}_{\mathit{eq}}-\mathit{tr}(M)\le 0\)
Cette inégalité est une première caractérisation de la direction d’écoulement \(M\) dans le cas singulier.
Intégration implicite du comportement#
On définit le prédicteur élastique à l’incrément n+1:
\({\sigma}^{e}=E:(\epsilon -{\epsilon}_{p}^{n})\)
avec \(E\) la matrice de Hooke.
Cas élastique:
Si \(F({\sigma}^{e},p)\le 0\) , alors \(\lambda =0\) est solution.
On a \(\Delta {\epsilon}_{p}={\epsilon}_{p}^{n+1}-{\epsilon}_{p}^{n}=0\) et \(\Delta p={p}^{n+1}-{p}^{n}=0\)
Cas plastique singulier:
On utilise \({\sigma}_{\mathit{eq}}=0\) :
\({\sigma}^{D}=0\Rightarrow 2\mu ({\epsilon}^{D}-{\epsilon}_{p}^{D})=0\)
On a donc \({\epsilon}^{D}={\epsilon}_{p}^{D}\) . D’où: \(\Delta {\epsilon}_{p}=\frac{1}{2\mu }{\sigma}_{D}^{e}\) . On a finalement:
\(\Delta p=\sqrt{\frac{1}{3\mu }{\sigma}_{D}^{e}:{\sigma}_{D}^{e}}\)
Or, on a nécessairement \(F(\sigma )=0\) . On a donc:
\(\mathit{tr}(\sigma )=\frac{R(p)}{\alpha}\)
On a aussi \(\mathit{tr}(\sigma )=\mathit{tr}({\sigma}^{e})-3K\Delta \mathit{tr}({\epsilon}_{p})\) . D’où:
\(\Delta \mathit{tr}({\epsilon}_{p})=\frac{1}{3K}(\mathit{tr}({\sigma}_{e})-\frac{R(p)}{\alpha})\)
On connaît donc les incréments de déformations plastiques volumique et déviatoriques. Au vu de la condition obtenue précédemment sur la direction d’écoulement dans le cas singulier, cette solution est valable si:
\(2\beta (p)\Delta {\epsilon}_{\mathit{eq}}-\Delta \mathit{tr}({\epsilon}_{p})\le 0\)
\(\frac{\beta (p){\sigma}_{\mathit{eq}}^{e}}{3\mu }\le \frac{1}{9K}(\mathit{tr}({\sigma}^{e})-\frac{R(p)}{\alpha})\)
Cas plastique régulier :
Dans le cas où \({\sigma}_{\mathit{eq}}\ne 0\) , on a:
\(\sigma ={\sigma}^{e}-K\Delta \mathit{tr}({\epsilon}_{p})I-2\mu \Delta {\epsilon}_{p}^{D}\)
En utilisant l’expression de la direction d’écoulement dans le cas régulier, il vient:
\(\Delta \mathit{tr}({\epsilon}_{p})=3\lambda \beta (p)\) et \(\Delta {\epsilon}_{p}^{D}=\frac{3}{2}\lambda \frac{{\sigma}^{D}}{{\sigma}_{\mathit{eq}}}\)
On peut donc exprimer les composantes sphérique et déviatorique du tenseur \(\sigma\) en fonction du prédicteur élastique et du multiplicateur plastique:
\(\mathit{tr}(\sigma )=\mathit{tr}({\sigma}^{e})-3Kb(p)\lambda\)
\({\sigma}^{D}={\sigma}_{D}^{e}-3\mu \frac{{\sigma}^{D}}{{\sigma}_{\mathit{eq}}}\lambda\)
En particulier, cette dernière équation permet d’exprimer \({\sigma}_{\mathit{eq}}\) :
\({\sigma}_{\mathit{eq}}={\sigma}_{\mathit{eq}}^{e}-3\mu \lambda\)
Or on sait que \(\Delta p=\lambda ` . En insérant ces expressions dans l'équation :math:`F=0\) , il vient:
\({\sigma}_{\mathit{eq}}^{e}-3\mu \lambda +3\alpha (\mathit{tr}({\sigma}^{e})-9K\lambda b({p}^{n}+\lambda ))-R({p}^{n}+\lambda )=0\)
La résolution de cette dernière équation permet de déterminer l’inconnue \(\lambda\) , et donc de déterminer toutes les inconnues du problème incrémental.
Par ailleurs on peut montrer que le problème incrémental admet une unique solution si la fonction définie par \(f(\lambda )={\sigma}_{\mathit{eq}}^{e}-3\mu \lambda +3\alpha (\mathit{tr}({\sigma}^{e})-9K\lambda b({p}^{n}+\lambda ))-R({p}^{n}+\lambda )\) est décroissante. Dans l’intégration numérique de la loi de comportement, on s’assure que cette condition est vérifiée à chaque incrément.
Calcul de l’opérateur tangent#
Les équations sont très similaires à celles de la loi associée, elles sont simplement écrites ici dans un cadre plus général.
On cherche à calculer la matrice cohérente: \(\frac{\partial \sigma }{\partial \varepsilon }=\frac{\partial \mathrm{s}}{\partial \varepsilon }+\frac{1}{3}\mathrm{I}\otimes \frac{\partial {I}_{1}}{\partial \varepsilon }\)
Il vient alors:
\(\frac{\partial \sigma }{\partial \epsilon }=\frac{1}{3}\mathrm{I}\otimes \frac{\partial {I}_{1}}{\partial \epsilon }+(1-q)\frac{\partial {s}^{e}}{\partial \epsilon }-{s}^{e}\otimes \frac{\partial q}{\partial \epsilon }\)
avec \(q=3\mu \frac{\lambda}{{\sigma}_{\mathit{eq}}^{e}}\)
Variation de la contrainte#
\(\frac{\partial q}{\partial \epsilon }=\frac{\partial {I}_{1}^{e}}{\partial \epsilon }\frac{\partial q}{\partial {I}_{1}^{e}}+\frac{\partial q}{\partial {I}_{1}^{e}}\frac{\partial {\sigma}_{\mathit{eq}}^{e}}{\partial \epsilon }\)
\(\frac{\partial {I}_{1}}{\partial \epsilon }=\frac{\partial {I}_{1}^{e}}{\partial \epsilon }\frac{\partial {I}_{1}}{\partial {I}_{1}^{e}}+\frac{\partial {\sigma}_{\mathit{eq}}^{e}}{\partial \epsilon }\frac{\partial {I}_{1}}{\partial {\sigma}_{\mathit{eq}}^{e}}\)
\(\frac{\partial {s}^{e}}{\partial \epsilon }=2\mu (\underline{\underline{\underline{\underline{\underline{I}}}}}-\frac{1}{3}I\otimes I)\)
Variation des invariants élastiques#
\(\frac{\partial {\sigma}_{\mathit{eq}}^{e}}{\partial \epsilon }=3\mu \frac{{s}^{e}}{{\sigma}_{\mathit{eq}}^{e}}\)
\(\frac{\partial {I}_{1}^{e}}{\partial \epsilon }=3KI\)
Variation des invariants de la contrainte#
\(\frac{\partial q}{\partial {\sigma}_{\mathit{eq}}^{e}}=3\mu (\frac{\partial p}{\partial {\sigma}_{\mathit{eq}}^{e}}\frac{1}{{\sigma}_{\mathit{eq}}^{e}}-\frac{\lambda}{{({\sigma}_{\mathit{eq}}^{e})}^{2}})\)
\(\frac{\partial q}{\partial {I}_{1}^{e}}=3\mu \frac{\partial p}{\partial {I}_{1}^{e}}\frac{1}{{({\sigma}_{\mathit{eq}}^{e})}^{2}}\)
\(\frac{\partial {I}_{1}}{\partial {\sigma}_{\mathit{eq}}^{e}}=-3K\frac{\partial p}{\partial {\sigma}_{\mathit{eq}}^{e}}(\frac{\partial \beta }{\partial p}\lambda +\beta (p))\)
\(\frac{\partial {I}_{1}}{\partial {I}_{1}^{e}}=1-3K\frac{\partial p}{\partial {I}_{1}^{e}}(\frac{\partial \beta }{\partial p}\lambda +\beta (p))\)
Variation de p#
\(\frac{\partial f}{\partial p}=\frac{-\partial R}{\partial p}-3\mu -9\alpha K(\frac{\partial \beta }{\partial p}\lambda +\beta (p))\)
\(\frac{\partial p}{\partial {I}_{1}^{e}}=\frac{-1}{\frac{\partial f}{\partial p}}.3\alpha\)
\(\frac{\partial p}{\partial {\sigma}_{\mathit{eq}}^{e}}=\frac{-1}{\frac{\partial f}{\partial p}}\)
Variables internes des lois Drucker-Prager associée et non associée#
Le modèle associé comporte trois variables internes:
\(\mathit{V1}\) est la déformation déviatorique plastique cumulée \(p\)
\(\mathit{V2}\) est la déformation volumique plastique cumulée \(\sum\Delta {\varepsilon}_{V}^{p}\)
\(\mathit{V3}\) est l’indicateur d’état (1 si \(\Delta p>0\) , 0 dans le cas contraire).
Le modèle non associé comporte neuf variables internes:
\(\mathit{V1}\) est la déformation déviatorique plastique cumulée \(p\)
\(\mathit{V2}\) est la déformation volumique plastique cumulée \(\sum\Delta {\varepsilon}_{V}^{p}\)
\(\mathit{V3}\) est l’indicateur d’état (2 si \(\Delta p>0\) et \({\sigma}_{\mathit{eq}}=0\) , 1 si \(\Delta p>0\) et \({\sigma}_{\mathit{eq}}\ne 0\) , 0 dans le cas contraire).
\(V4\) À \(V9\) sont les déformations anélastiques.
Dans sa version non associée, la loi de comportement est écrite en déformations anélastiques (et non pas en contrainte incrémentale). Pour appliquer une condition initiale en contrainte, il est donc nécessaire d’utiliser ces variables internes \(V4\) à \(V9\) . Les cas-tests ssnv168e et ssnv168f en sont un exemple.
Indicateur de localisation de Rice pour la loi Drucker-Prager#
On définit l’indicateur de localisation du critère de Rice dans le cadre de la loi de comportement Drucker-Prager. Mais la définition d’un indicateur de localisation peut-être utilisé, de façon plus générale, dans des études en mécanique de la rupture, mécanique de l’endommagement, théorie de la bifurcation, mécanique des sols et mécanique des roches (et plus globalement dans le cadre des matériaux à loi de comportement adoucissante).
Cet indicateur défini un état à partir duquel l’évolution du système mécanique étudié (équations, d’équilibre, loi de comportement) peut perdre son caractère d’unicité. Cette théorie permet, en d’autres termes:
le calcul de l’état d’initiation possible de la localisation qui est perçu comme la limite de validité des calculs par éléments finis classiques;
la détermination «qualitative» des angles d’orientation des zones de localisation.
Le critère de localisation constitue une limite de fiabilité des calculs par éléments finis «classiques».
Les différentes façons d’étudier la localisation#
Dans le cadre des études menées en mécanique des sols, on a constaté une forte dépendance de la solution numérique en fonction de la discrétisation par éléments finis. Il apparaît une concentration de valeurs élevées des déformations plastiques cumulées au niveau des éléments finis et on note que cette «zone de localisation» change brutalement avec le raffinement du maillage. Ce phénomène de localisation est source de problèmes numériques et engendre des problèmes de convergences au sens des éléments finis.
La localisation peut être interprétée comme un phénomène instable, précurseur de mécanisme de rupture, caractérisant certains types de matériaux sollicités dans le domaine inélastique. Pour étudier les instabilités liées à la localisation on distingue, d’une part, les classes de matériaux à comportement dépendant du temps et d’autre part, celles ne dépendant pas du temps. Pour les matériaux à comportement indépendant du temps, l’approche communément utilisée est la méthode dite par bifurcation (c’est à cette méthode que l’on s’intéresse dans cette note). Elle consiste à analyser les pertes d’unicité du problème en vitesses. Pour les matériaux à comportement dépendant du temps, l’unicité du problème en vitesses est souvent garantie et ceci n’empêche pas l’observation des instabilités lors de leur déformation. Pour ces matériaux, on doit alors avoir recours à d’autres approches. La plus couramment utilisée est l’approche par perturbation. Cette approche ne sera pas traitée dans cette note, mais pour plus d’informations consulter les notes [bib1], [bib2].
Rudnicki et Rice [bib3] ont démontré que l’étude de la localisation des déformations en mécanique des roches s’inscrit dans le cadre de la théorie de la bifurcation. Celle-ci est basée sur la notion d’équilibre instable. Rice [bib4] considère que le point de bifurcation marque la fin du régime stable. Le début de la localisation est associé à une instabilité rhéologique du système et cette instabilité correspond localement à la perte d’ellipticité des équations qui gouvernent l’équilibre incrémental continu en vitesses. Rice propose ainsi un critère dit de «bifurcation par localisation» qui permet de détecter l’état à partir duquel, la solution des équations mathématiques qui gouvernent le problème aux limites considérées et l’évolution du système mécanique étudié (équations, d’équilibre, loi de comportement) perdent leur caractère d’unicité. Cette théorie permet le calcul de l’état d’initiation de la localisation qui est perçu comme la limite de validité des calculs par éléments finis classiques.
Approche théorique#
Ecriture du problème en vitesse#
On considère une structure occupant, à un instant \(t\) , l’ouvert \(\Omega\) de \({\mathrm{\Re }}^{3}\) . Le problème en vitesse consiste à trouver le champ des vitesses de déplacements \(v\) lorsque la structure est soumise aux vitesses de forces volumiques \({\dot{f}}_{d}\) , aux vitesses de déplacements imposés \({v}_{d}\) sur une partie \({\partial}_{1}\Omega\) de la frontière et aux vitesses d’efforts surfaciques \({\dot{F}}_{d}\) sur la partie complémentaire \({\partial}_{2}\Omega\) .
Dans l’écriture locale du problème, le champ des vitesses de déplacements \(v\) doit ainsi vérifier le problème:
\(v\) suffisamment régulier et \(v={v}_{d}\) sur \({\partial}_{1}\Omega\)
Les équations d’équilibre:
\(div[\mathrm{L}:\epsilon (v)]+{\dot{\mathrm{f}}}_{d}=0\) sur \(\Omega\)
\(\mathrm{L}:\varepsilon (v).\mathrm{n}={\dot{\mathrm{F}}}_{d}\) sur \({\partial}_{2}\Omega\)
\(\mathrm{n}\) étant la normale unitaire sortante à \({\partial}_{2}\Omega\) .
Les conditions de compatibilité(on se limite ici aux petites perturbations):
\(\epsilon (v)=\frac{1}{2}\left[\nabla v+(\nabla v{)}^{T}\right]\)
où l’opérateur \(L\) est défini de façon générale pour les lois de comportement écrites sous forme incrémentale par la relation:
\(\dot{\sigma}=L(\epsilon ,V):\dot{\epsilon}\)
avec:
\(L=\left\{ \begin{array}{ccc}E& & \text{si }F<0\text{ ou }F=0\text{ et }\frac{b:E:\dot{\epsilon}}{h}\le 0\\ H=& E-\frac{(E:a)\otimes (b:E)}{h}& \text{si }F=0\text{ et }\frac{b:E:\dot{\epsilon}}{h}>0\end{array} \right.\)
où \(\sigma\) est la contrainte, \(\epsilon\) la déformation totale, \(\mathrm{V}\) un ensemble de variables internes et \(F\) la surface seuil de plasticité. Les expressions de \(a,b,E\) et h dépendent de la formulation de la loi de comportement.
Résultats d’existence et d’unicité, perte d’ellipticité#
Nous donnons dans ce chapitre quelques résultats sans démonstrations. La référence pour ces démonstrations est toutefois spécifiée.
Une condition suffisante d’existence et d’unicité du problème précédent est: \(\dot{\sigma}:\dot{\epsilon}>0\) . Cette inégalité peut être interprétée comme une définition, dans le cas tridimensionnel, du non-adoucissement. La démonstration est faite par Hill [bib5] pour les matériaux standards et par Benallal [bib1] pour les matériaux non-standards.
La perte d’ellipticité correspond à l’instant pour lequel l’opérateur \(N.H.N\) devient singulier pour une direction \(N\) en un point de la structure. Cette condition est équivalente à la condition: \(\det\left(\mathrm{N}.\mathrm{H}.\mathrm{N}\right)=0\) . Il s’agit de la condition de «bifurcation continue» [1]
au sens de Rice aussi appelé tenseur acoustique. Rice et Rudnicki [bib3] montrent que cette condition de perte d’ellipticité du problème local en vitesse est une condition nécessaire à la bifurcation «continue ou discontinue» [2]
pour le solide. Les conditions aux limites ne jouent aucun rôle, seule la loi de comportement définit les conditions de localisation (seuil de localisation et orientation de la surface de localisation.
Les bifurcations continues fournissent ainsi la limite inférieure de la gamme de déformation pour laquelle les bifurcations discontinues peuvent se produire.
Résolution analytique pour le cas bidimensionnel.#
On pose \(N=({N}_{1},{N}_{2},0)\) avec \({N}_{1}^{2}+{N}_{2}^{2}=1\)
On a alors: \(N.H.N=\left[\begin{array}{ccc}{A}_{11}& {A}_{12}& 0\\ {A}_{21}& {A}_{22}& 0\\ 0& 0& C\end{array}\right]\) où Ortiz [bib6] montre que:
\(C={N}_{1}^{2}{H}_{1313}+{N}_{2}^{2}{H}_{2323}>0\)
\({A}_{11}={N}_{1}^{2}{H}_{1111}+{N}_{1}{N}_{2}({H}_{1112}+{H}_{1211})+{N}_{2}^{2}{H}_{1212}\)
\({A}_{22}={N}_{1}^{2}{H}_{1212}+{N}_{1}{N}_{2}({H}_{1222}+{H}_{2212})+{N}_{2}^{2}{H}_{2222}\)
\({A}_{12}={N}_{1}^{2}{H}_{1112}+{N}_{1}{N}_{2}({H}_{1122}+{H}_{1212})+{N}_{2}^{2}{H}_{1222}\)
\({A}_{21}={N}_{1}^{2}{H}_{1211}+{N}_{1}{N}_{2}({H}_{1212}+{H}_{2211})+{N}_{2}^{2}{H}_{2212}\)
Il suffit ainsi d’étudier le signe de \(\det(A)\) comme précisé par Doghri [bib7]:
\(\det(A)={a}_{0}{N}_{1}^{4}+{a}_{1}{N}_{1}^{3}{N}_{2}+{a}_{2}{N}_{1}^{2}{N}_{2}^{2}+{a}_{3}{N}_{1}{N}_{2}^{3}+{a}_{4}{N}_{2}^{4}\)
avec:
\({a}_{0}={H}_{1111}{H}_{1212}-{H}_{1112}{H}_{1211}\)
\({a}_{1}={H}_{1111}({H}_{1222}+{H}_{2212})-{H}_{1112}{H}_{2211}-{H}_{1122}{H}_{1211}\)
\({a}_{2}={H}_{1111}{H}_{2222}+{H}_{1112}{H}_{1222}+{H}_{1211}{H}_{2212}-{H}_{1122}{H}_{1212}-{H}_{1122}{H}_{2211}-{H}_{1212}{H}_{2211}\)
\({a}_{3}={H}_{2222}({H}_{1112}+{H}_{1211})-{H}_{1122}{H}_{2212}-{H}_{1222}{H}_{2211}\)
\({a}_{4}={H}_{1212}{H}_{2222}-{H}_{1222}{H}_{2212}\)
On pose alors \({N}_{1}=\cos\theta\) et \({N}_{2}=\sin\theta\) avec \(\theta \in ]\begin{array}{cc}-\frac{\pi}{2};& +\frac{\pi}{2}\end{array}]\) . On distingue alors deux cas:
si \(\theta \text{=+}\frac{\pi}{2}\) alors \(\det(A)=0\) si \({a}_{4}=0\) ;
si \(\theta \ne +\frac{\pi}{2}\) alors \(\det(A)=0\) si \(f(x)={a}_{4}{x}^{4}+{a}_{3}{x}^{3}+{a}_{2}{x}^{2}+{a}_{1}x+{a}_{0}=0\) avec \(x=\tan\theta\) .
Calcul des racines#
Pour résoudre un polynôme de degré n (comme celui défini ci-dessus, où n=4) on propose d’utiliser la méthode dite «Companion Matrix Polynomial». Le principe de cette méthode consiste à chercher les valeurs propres de la matrice (de type Hessenberg) d’ordre n associée au polynôme. Si l’on considère le polynôme \(P(x)={x}^{n}+{a}_{n-1}{x}^{n-1}+...+{a}_{k}{x}^{k}+...+{a}_{1}x+{a}_{0}\) . Chercher les racines de ce polynôme revient à chercher les valeurs propres de la matrice:
\(\left[\begin{array}{cccccc}0& 0& 0& 0& 0& -{a}_{0}\\ 1& 0& 0& 0& 0& -{a}_{1}\\ 0& 1& 0& 0& 0& ...\\ 0& 0& 1& 0& 0& -{a}_{k}\\ 0& 0& 0& 1& 0& ...\\ 0& 0& 0& 0& 1& -{a}_{n-1}\end{array}\right]\)
Cet indicateur est calculé par l’option INDL_ELGA de CALC_CHAMP [u4.81.04]. Il produit en chaque point d’intégration 5 composantes: la première est l’indicateur de localisation valant 0 si \(det(\mathit{N.H.N})>0\) (pas de localisation), et valant 1 sinon, ce qui correspond a une possibilité de localisation. Les autres composantes fournissent les directions de localisation.
Calculs de sensibilité#
L’analyse de sensibilité porte uniquement sur la version associée de la formulation décrite au chapitre r7.01.16-formulation-version-associee .
Sensibilité aux données matériaux#
Le problème direct#
Nous nous plaçons dans cette partie dans le cadre de la résolution de calculs non-linéaires.
Le calcul statique non-linéaire est résolu incrémentalement. Il nécessite donc à chaque pas de charge \(i\in \lbrace 1,I\rbrace\) la résolution du système d’équation non-linéaire :
avec
\({\mathrm{w}}_{k}\) est la fonction de forme du \(k\) ième degré de liberté de la structure modélisée,
\((\mathrm{R}({\mathrm{u}}_{i},{t}_{i}))\) est le vecteur des forces nodales.
La résolution de ce système se fait par la méthode de Newton-Raphson:
où \({K}_{i}^{n}=\frac{\partial R}{\partial u}{\mid}_{({u}_{i}^{n},{t}_{i})}\) est la matrice tangente au pas de charge \(i\) et à l’itération de Newton \(n\) .
La solution est donc donnée par:
avec \(N\) , le nombre d’itérations de Newton qui a été nécessaire pour atteindre la convergence.
Le calcul dérivé#
Préliminaires#
Dans le cadre du calcul de sensibilité, il est nécessaire d’insister sur les dépendances d’une grandeur par rapport aux autres. Nous allons ainsi expliciter que les résultats du calcul précédent dépendent d’un paramètre \(\Phi\) donné (module d’Young, limite d’élasticité, masse volumique, …) et cela de la manière suivante:
\({u}_{i}={u}_{i}(\Phi )\) , \({\lambda}_{i}={\lambda}_{i}(\Phi )\) .
Mais cela n’est pas suffisant. Aussi nous plaçons-nous dans le cadre d’un calcul incrémental avec loi de comportement de type Drucker-Prager.. Si l’on considère les inter-dépendances des paramètres à un niveau algorithmique, on peut écrire:
\(R=R({\sigma}_{i-1}(\Phi ),{p}_{i-1}(\Phi ),\Delta u(\Phi ))\)
\({\sigma}_{i}={\sigma}_{i-1}(\Phi )+\Delta \sigma ({\sigma}_{i-1}(\Phi ),{p}_{i-1}(\Phi ),\Delta u(\Phi ),\Phi )\)
\({p}_{i}={p}_{i-1}(\Phi )+\mathrm{\Delta p}({\sigma}_{i-1}(\Phi ),{p}_{i-1}(\Phi ),\Delta u(\Phi ),\Phi )\)
Où \(\Delta u\) est l’incrément de déplacement à convergence au pas de charge \(i\) .
Précisons le sens des notations que nous utiliserons pour les dérivées:
\(\frac{\partial X}{\partial Y}\) désigne la dérivée partielle explicite de \(X\) par rapport à \(Y\) ,
\(X{,}_{Y}\) désigne la variation totale de \(X\) par rapport à \(Y\) .
Dérivation de l’équilibre#
Compte tenu des remarques précédentes, exprimons la variation totale par rapport à \(\Phi\) :
Remarquons qu’ici \(\frac{\partial R}{\partial \Phi }=0\) : \(R\) ne dépend pas explicitement de \(\Phi\) mais implicitement comme nous le verrons en détail dans la suite.
Soit:
Où
\({K}_{i}^{N}\) est la dernière matrice tangente utilisée pour atteindre la convergence dans les itérations de Newton,
\(R{,}_{\Phi}{\mid}_{\mathrm{\Delta u}\ne \mathrm{\Delta u}(\Phi )}\) est la variation totale de \(R\) , sans tenir compte de la dépendance de \(\mathrm{\Delta u}\) par rapport à \(\Phi\) .
Le problème réside maintenant dans le calcul de \(R{,}_{\Phi}{\mid}_{\mathrm{\Delta u}\ne \mathrm{\Delta u}(\Phi )}\) .
Remarque:
Dans (3852), on a utilisé le fait que \({K}_{i}^{N}=\frac{\partial R({u}_{i},{t}_{i})}{\partial \Delta u}\) alors que dans (4583) on l’a défini par \({K}_{i}^{N}=\frac{\partial R({u}_{i},{t}_{i})}{\partial {u}_{i}^{N}}\) *. On a bien équivalence de ces deux définitions dans la* mesure où \({u}_{i}={u}_{i-1}+\mathrm{\Delta u}\) et que \(R\) dépend effectivement de \(\Delta u\) (et aussi bien sûr de \({\sigma}_{i-1}\) et \({p}_{i-1}\) ).
Remarque:
Si on dérive par rapport à \(\Phi\) directement (4583), on trouve \({K}^{n}=\frac{\partial {u}^{n+1}}{\partial \Phi }+{B}^{t}\lambda ,\Phi =-R{,}_{\Phi /\mathrm{\Delta u}\ne \mathrm{\Delta u}/\Phi }-{K}^{n}{,}_{\Phi}{\mathrm{\delta u}}^{n+1}\) *. Ce qui est la même chose à convergence et fait apparaître que l’erreur sur* \(\frac{\partial u}{\partial \Phi }\) dépend de \({K}^{-1}K{,}_{\Phi}\) .
Calcul de la dérivée de la loi de comportement#
Dans la suite, par souci de clarté, nous abandonnerons les indices \(i-1\) .
D’après (4582), on peut réécrire \(R{,}_{\Phi}{\mid}_{\mathrm{\Delta u}\ne \mathrm{\Delta u}(\Phi )}\) sous la forme:
On doit donc calculer \(\Delta \sigma {,}_{\Phi}{\mid}_{\Delta u\ne \Delta u(\Phi )}\) . Pour ce faire, nous allons utiliser les expressions qui interviennent dans l’intégration numérique de la loi de comportement.
Cas de l’élasticité linéaire#
Dans le cadre de l’élasticité linéaire, la loi de comportement s’exprime par:
\(\left\{ \begin{array}{}\Delta \tilde{\sigma}=\mathrm{2\mu }.\tilde{\epsilon}(\Delta u)\\ \text{Tr}(\Delta \sigma )=\mathrm{3K}.\text{Tr}(\epsilon (\Delta u))\end{array}\right.\)
ou bien:
où \(\text{Id}\) est le tenseur identité d’ordre 2.
Alors, en calculant la variation totale de (3854) par rapport à \(\Phi\) on obtient:
Soit:
Cas de l’élastoplasticité de type Drucker-Prager#
La loi de comportement de type Drucker-Prager s’écrit:
où \(S\) est le tenseur des souplesses élastiques et \(R\) est le critère de plasticité défini par:
dans le cas d’un écrouissage linéaire:
\(\begin{array}{}R(p)=h\cdot p+{\sigma}^{y}\text{pour}0\le p\le {p}_{\text{ultm}}\\ R(p)=h\cdot {p}_{\text{ultm}}\text{pour}p\ge {p}_{\text{ultm}}\end{array}\)
dans le cas d’un écrouissage parabolique:
\(\begin{array}{}R(p)={\sigma}^{y}\cdot (\text{1-}(\text{1-}\sqrt{\frac{{\sigma}_{\text{ultm}}^{y}}{{\sigma}^{y}}})\cdot \frac{p}{{p}_{\text{ultm}}}{)}^{2}\text{pour}0\le p\le {p}_{\text{ultm}}\\ R(p)={\sigma}_{\text{ultm}}^{y}\text{pour}p\ge {p}_{\text{ultm}}\end{array}\)
En termes numériques, cette loi de comportement est intégrée à l’aide d’un algorithme de retour radial: on fait une prédiction élastique (notée \({\sigma}^{e}\) ) que l’on corrige si le seuil est violé. On écrit donc:
Nous allons distinguer deux cas.
1er cas : \(\mathrm{\Delta p}=0\)
Ce qui revient à dire que lors du présent pas de charge, le point de Gauss considéré n’a pas vu d’accroissement de sa plastification. On se retrouve alors dans le cas de l’élasticité linéaire:
2eme cas : \(\mathrm{\Delta p}>0\)
Compte tenu des dépendances entre variables dans (3857), on peut écrire:
En outre, en accord avec l’intégration algorithmique de la loi, nous allons séparer parties déviatorique et hydrostatique.
Et donc, on calcule:
\(\frac{\partial \Delta \sigma }{\partial \Phi }\)
\(\frac{\partial \Delta \tilde{\sigma}}{\partial \Phi }=\frac{\partial \mathrm{2\mu }}{\partial \Phi }\cdot \tilde{\epsilon}(\mathrm{\Delta u})-\frac{\partial \mathrm{3\mu }}{\partial \Phi }\cdot \frac{\mathrm{\Delta p}}{{\sigma}_{\text{eq}}^{e}}\cdot {\tilde{\sigma}}^{e}-\mathrm{3\mu }\cdot \frac{\frac{\partial \mathrm{\Delta p}}{\partial \Phi }}{{\sigma}_{\text{eq}}^{e}}\cdot {\tilde{\sigma}}^{e}+\mathrm{3\mu }\cdot \frac{\mathrm{\Delta p}\cdot \frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \Phi }}{{\sigma}_{\text{eq}}^{{e}^{2}}}\cdot {\tilde{\sigma}}^{e}-\mathrm{3\mu }\cdot \frac{\mathrm{\Delta p}}{{\sigma}_{\text{eq}}^{e}}\cdot \frac{\partial {\tilde{\sigma}}^{e}}{\partial \Phi }\)
\(\frac{\partial \text{Tr}(\Delta \sigma )}{\partial \Phi }=\frac{\partial \mathrm{3K}}{\partial \Phi }\cdot \text{Tr}(\epsilon (\Delta u))-\frac{\partial \mathrm{9K}}{\partial \Phi }\cdot {\rm A}\cdot \Delta p-\mathrm{9K}\cdot \frac{\partial {\rm A}}{\partial \Phi }\cdot \mathrm{\Delta p}-\mathrm{9K}\cdot {\rm A}\cdot \frac{\partial \mathrm{\Delta p}}{\partial \Phi }\)
\(\frac{\partial \Delta \sigma }{\partial \sigma }\)
\(\frac{\partial \Delta \tilde{\sigma}}{\partial \sigma }=-\mathrm{3\mu }\cdot \frac{\frac{\partial \mathrm{\Delta p}}{\partial \sigma }}{{\sigma}_{\text{eq}}^{e}}\otimes {\tilde{\sigma}}^{e}+\mathrm{3\mu }\cdot \frac{\mathrm{\Delta p}}{{\sigma}_{\text{eq}}^{{e}^{2}}}\cdot \frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \sigma }\otimes {\tilde{\sigma}}^{e}-\mathrm{3\mu }\cdot \frac{\mathrm{\Delta p}}{{\sigma}_{\text{eq}}^{e}}\cdot J\)
où \(J\) est l’opérateur déviatorique défini par: \(J:\sigma =\tilde{\sigma}\)
\(\frac{\partial \text{Tr}(\Delta \sigma )}{\partial \sigma }=-\mathrm{9K}\cdot {\rm A}\cdot \frac{\partial \mathrm{\Delta p}}{\partial \sigma }\)
\(\frac{\partial \mathrm{\Delta \sigma }}{\partial p}\)
\(\frac{\partial \Delta \tilde{\sigma}}{\partial p}=-\frac{\mathrm{3\mu }}{{\sigma}_{\text{eq}}^{e}}\cdot \frac{\partial \mathrm{\Delta p}}{\partial p}\cdot {\tilde{\sigma}}^{e}\)
\(\frac{\partial \text{Tr}(\mathrm{\Delta \sigma })}{\partial p}=-9\cdot K\cdot {\rm A}\cdot \frac{\partial \mathrm{\Delta p}}{\partial p}\)
\(\mathrm{\Delta p}{,}_{\Phi}\)
On utilise le fait que: \((\sigma +\mathit{\Delta \sigma }{)}_{\text{eq}}=(\sigma +\mathit{\Delta \sigma }{)}_{\text{eq}}^{e}-3\mu \cdot \mathit{\Delta p}\)
\({\mathrm{\Delta p}}_{,\Phi }=\frac{1}{\mathrm{3\mu }}\cdot ((\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}^{{e}_{,\Phi }}-(\sigma +\mathrm{\Delta \sigma }{)}_{{\text{eq}}_{,\Phi }}-\frac{\partial \mathrm{3\mu }}{\partial \Phi }\cdot \mathrm{\Delta p})\)
Remarque:
Dans ces calculs ont été ou doivent être utilisés les résultats suivants: |
|
\(\frac{\partial {\tilde{\sigma}}^{e}}{\partial \Phi }=\frac{\partial \mathrm{2\mu }}{\partial \Phi }\cdot \tilde{\epsilon}(\mathrm{\Delta u})\) |
\(\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \Phi }=\frac{3}{2}\cdot \frac{(\frac{\partial \mathrm{2\mu }}{\partial \Phi }\cdot \tilde{\epsilon}(\mathrm{\Delta u})):(\tilde{\sigma}+\mathrm{2\mu }.\tilde{\epsilon}(\mathrm{\Delta u}))}{{\sigma}_{\text{eq}}^{e}}\) |
Tenseur d’ordre 2 |
Scalaire |
\(\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \sigma }=\frac{3}{2}\cdot \frac{{\tilde{\sigma}}^{e}}{{\sigma}_{\text{eq}}^{e}}\) |
\(\frac{\partial {\tilde{\sigma}}^{e}}{\partial \sigma }=J\) |
Tenseur d’ordre 2 |
Tenseur d’ordre 4 |
\(\frac{\partial \text{Tr}({\sigma}^{e})}{\partial \sigma }=\text{Id}\) |
\(\frac{\partial \text{Tr}({\sigma}^{e})}{\partial \Phi }=\frac{\partial \mathrm{3K}}{\partial \Phi }\cdot \text{Tr}(\epsilon (\mathrm{\Delta u}))\) |
Tenseur d’ordre 2 |
Scalaire |
On doit également calculer les dérivées partielles de l’incrément de déformation plastique cumulée par rapport aux paramètres matériaux, aux contraintes et à la déformation plastique cumulée (cf Annexe)
Celles-ci sont obtenues en dérivant l’équation résolue pour calculer l’incrément de déformation plastique cumulée lors du calcul direct.
Calcul de la dérivée du déplacement#
Une fois \(\Delta \sigma {,}_{\Phi}{|}_{\Delta \mathrm{u}\ne \Delta \mathrm{u}(\Phi )}\) calculé, on peut constituer le second membre \(R{,}_{\Phi}{\mid}_{\Delta u\ne \Delta u(\Phi )}\) en utilisant (3853). On résout alors le système (3852) et l’on obtient l’incrément de déplacement dérivé par rapport à \(\Phi\) .
Calcul de la dérivée des autres grandeurs#
Maintenant que l’on dispose de \(\Delta u{,}_{\Phi}\) , on doit calculer la dérivée des autres grandeurs. On sépare encore deux cas:
Élasticité linéaire
D’après (3857), on calcule comme suit la dérivée de l’incrément de contrainte:
\(\Delta \sigma {,}_{\Phi}=\Delta \sigma {,}_{\Phi}{\mid}_{\Delta u\ne \Delta u(\Phi )}+\mathrm{2\mu }.\tilde{\epsilon}(\Delta u{,}_{\Phi})+K.\text{Tr}(\epsilon (\Delta u{,}_{\Phi})).\text{Id}\)
L’incrément de déformation plastique cumulée, quant à lui, ne voit pas d’évolution:
\(\mathrm{\Delta p}{,}_{\Phi}=0\)
Elastoplasticité de type Drucker-Prager
Si \(\mathrm{\Delta p}=0\) , on retrouve le cas précédent.
Sinon, on obtient d’après (3858):
\(\Delta \sigma {,}_{\Phi}=\Delta \sigma {,}_{\Phi}{\mid}_{\Delta u\ne \Delta u(\Phi )}+\frac{\partial \Delta \sigma }{\partial \epsilon (\Delta u)}:\epsilon (\Delta u{,}_{\Phi})\)
Et pour la déformation plastique cumulée, on utilise la relation suivante:
\((\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}=(\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}^{e}-\mathrm{3\mu }\cdot \mathrm{\Delta p}\)
Celle-ci nous permet d’écrire que:
\({\mathrm{\Delta p}}_{,\Phi }=\frac{1}{\mathrm{3\mu }}\cdot ((\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}^{{e}_{,\Phi }}-(\sigma +\mathrm{\Delta \sigma }{)}_{{\text{eq}}_{,\Phi }}-\frac{\partial \mathrm{3\mu }}{\partial \Phi }\cdot \mathrm{\Delta p})\)
Les contraintes équivalentes sensibles sont calculées ainsi:
\((\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}^{{e}_{,\Phi }}=\frac{3}{2(\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}^{e}}\cdot ({\tilde{\sigma}}_{,\Phi }+\frac{\partial \mathrm{2\mu }}{\partial \Phi }\cdot \tilde{\epsilon}(\mathrm{\Delta u})+\mathrm{2\mu }\cdot \tilde{\epsilon}({\mathrm{\Delta u}}_{,\Phi })):(\tilde{\sigma}+\mathrm{2\mu }\cdot \tilde{\epsilon}(\mathrm{\Delta u}))\)
\((\sigma +\mathrm{\Delta \sigma }{)}_{{\text{eq}}_{,\Phi }}=\frac{3}{2(\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}}\cdot ({\tilde{\sigma}}_{,\Phi }+\Delta {\tilde{\sigma}}_{,\Phi }):(\tilde{\sigma}+\Delta \tilde{\sigma})\)
Une fois que tous ces calculs sont terminés, on réactualise toutes les grandeurs dérivées et on passe au pas de charge suivant.
Synthèse#
Pour résumer les paragraphes précédents, on représente les différentes étapes du calcul par le diagramme suivant:
Résolution du système (3852)
Sensibilité au chargement#
La démarche est ici assez proche de celle du paragraphe précédent. Nous la développons néanmoins entièrement dans un souci de clarté, afin que ce paragraphe puisse être lu de façon indépendante.
Le problème direct : expression du chargement#
Jusqu’à maintenant nous avons exprimé le problème direct sous la forme:
Les chargements sont rassemblés au second membre et comprennent les forces imposées \({L}_{i}\) et les déplacements imposés \({u}_{i}^{d}\) .
Supposons que le chargement en force imposée \({L}_{i}\) dépende d’un paramètre scalaire \(\alpha\) de la manière suivante:
Où
\({L}_{i}^{1}\) est un vecteur indépendant de \(\alpha\) ,
\({L}_{i}^{2}\) dépend linéairement de \(\alpha\) .
On désire calculer la sensibilité des résultats du calcul direct à une variation du paramètre \(\alpha\) .
Le problème dérivé#
Dérivation de l’équilibre#
Comme dans le chapitre précédent, en tenant compte des dépendances entre les différents champs, on dérive l’équilibre (3862) par rapport \(\alpha\) :
On a utilisé le fait que \({L}_{i}^{2}\) dépend linéairement de \(\alpha\) .
Soit:
Où
\({K}_{i}^{N}\) est la dernière matrice tangente utilisée pour atteindre la convergence dans les itérations de Newton,
\(R{,}_{\alpha}{\mid}_{\Delta u\ne \Delta u(\alpha )}\) est la variation totale de \(R\) , sans tenir compte de la dépendance de \(\Delta u\) par rapport à \(\alpha\) .
Le problème réside comme précédemment dans le calcul de \(R{,}_{\alpha}{\mid}_{\Delta u\ne \Delta u(\alpha )}\) .
Calcul de la dérivée de la loi de comportement#
D’après (4582), on peut réécrire \(R{,}_{\alpha}{\mid}_{\Delta u\ne \Delta u(\alpha )}\) sous la forme:
Pour ce faire, nous allons utiliser les expressions qui interviennent dans l’intégration numérique de la loi de comportement pour calculer \(\Delta \sigma {,}_{\alpha}{\mid}_{\Delta u\ne \Delta u(\alpha )}\) .
Cas de l’élasticité linéaire#
Dans le cadre de l’élasticité linéaire, la loi de comportement s’exprime par:
où \(\text{Id}\) est le tenseur identité d’ordre 2.
Alors, en calculant la variation totale de (3867) par rapport à \(\alpha\) on obtient:
Soit:
\(\Delta \sigma {,}_{\alpha}{\mid}_{\Delta u\ne \Delta u(\alpha )}=0.\)
Cas de l’élastoplasticité de type Drucker-Prager#
Comme précédemment, nous allons distinguer deux cas.
1er cas : \(\mathrm{\Delta p}=0\)
Ce qui revient à dire que lors du présent pas de charge, le point de Gauss considéré n’a pas vu d’accroissement de sa plastification. On se retrouve alors dans le cas de l’élasticité linéaire:
\(\Delta \sigma {,}_{\alpha}{\mid}_{\Delta u\ne \Delta u(\alpha )}=0.\)
2eme cas : \(\mathrm{\Delta p}>0\)
Compte tenu des dépendances entre variables, on peut écrire:
\(\left\{ \begin{array}{ccc}\mathit{\Delta \sigma }{,}_{\alpha}& =& \frac{\partial \mathit{\Delta \sigma }}{\partial \alpha }+\frac{\partial \mathit{\Delta \sigma }}{\partial \sigma }\cdot \sigma {,}_{\alpha}+\frac{\partial \mathit{\Delta \sigma }}{\partial p}\cdot p{,}_{\alpha}+\frac{\partial \mathit{\Delta \sigma }}{\partial \epsilon (\mathit{\Delta u})}\cdot \epsilon (\mathit{\Delta u}{,}_{\alpha})\\ \mathit{\Delta p}{,}_{\alpha}& =& \frac{\partial \mathit{\Delta p}}{\partial \alpha }+\frac{\partial \mathit{\Delta p}}{\partial \sigma }\cdot \sigma {,}_{\alpha}+\frac{\partial \mathit{\Delta p}}{\partial p}\cdot p{,}_{\alpha}+\frac{\partial \mathit{\Delta p}}{\partial \epsilon (\mathit{\Delta u})}\cdot \epsilon (\mathit{\Delta u}){,}_{\alpha}\end{array}\right.\)
En outre, en accord avec l’intégration algorithmique de la loi, nous allons séparer parties déviatorique et hydrostatique.
\(\left\{ \begin{array}{ccc}\mathit{\Delta \sigma }{,}_{\alpha}{\mid}_{\mathit{\Delta u}\mathrm{\ne }\mathit{\Delta u}(\alpha )}& =& \frac{\partial \Delta \tilde{\sigma}}{\partial \alpha }+\frac{1}{3}\cdot \frac{\partial \text{Tr}(\mathit{\Delta \sigma })}{\partial \alpha }\cdot \text{Id}\\ & +& \frac{\partial \Delta \tilde{\sigma}}{\partial \sigma }\cdot \sigma {,}_{\alpha}+\frac{1}{3}\cdot \frac{\partial \text{Tr}(\mathit{\Delta \sigma })}{\partial \sigma }\cdot \text{Id}\cdot \sigma {,}_{\alpha}\\ & +& \frac{\partial \Delta \tilde{\sigma}}{\partial p}\cdot p{,}_{\alpha}+\frac{1}{3}\cdot \frac{\partial \text{Tr}(\mathit{\Delta \sigma })}{\partial p}\cdot \text{Id}\cdot p{,}_{\alpha}\\ \mathit{\Delta p}{,}_{\alpha}{\mid}_{\mathit{\Delta u}\mathrm{\ne }\mathit{\Delta u}(\alpha )}& =& \frac{\partial \mathit{\Delta p}}{\partial \alpha }+\frac{\partial \mathit{\Delta p}}{\partial \sigma }\cdot \sigma {,}_{\alpha}+\frac{\partial \mathit{\Delta p}}{\partial p}\cdot p{,}_{\alpha}\end{array}\right.\)
Et donc, on calcule:
\(\frac{\partial \Delta \sigma }{\partial \alpha }\)
Dans la mesure où l’on n’a pas de dépendance explicite de \(\Delta \sigma ` par rapport à :math:\)alpha` , on obtient:
\(\frac{\partial \Delta \tilde{\sigma}}{\partial \alpha }=0.\)
\(\frac{\partial \text{Tr}(\Delta \sigma )}{\partial \alpha }=0.\)
\(\frac{\partial \Delta \sigma }{\partial \sigma }\)
\(\frac{\partial \Delta \tilde{\sigma}}{\partial \sigma }=-\mathrm{3\mu }\cdot \frac{\frac{\partial \mathit{\Delta p}}{\partial \sigma }}{{\sigma}_{\text{eq}}^{e}}\otimes {\tilde{\sigma}}^{e}+\mathrm{3\mu }\cdot \frac{\mathit{\Delta p}}{{\sigma}_{\text{eq}}^{{e}^{2}}}\cdot \frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \sigma }\otimes {\tilde{\sigma}}^{e}-\mathrm{3\mu }\cdot \frac{\mathit{\Delta p}}{{\sigma}_{\text{eq}}^{e}}\cdot J\)
où \(\mathrm{J}\) est l’opérateur déviatorique.
\(\frac{\partial \text{Tr}(\mathit{\Delta \sigma })}{\partial \sigma }=-\mathrm{9K}\cdot {\rm A}\cdot \frac{\partial \mathit{\Delta p}}{\partial \sigma }\)
\(\frac{\partial \Delta \sigma }{\partial p}\)
\(\frac{\partial \Delta \tilde{\sigma}}{\partial p}=-\frac{\mathrm{3\mu }}{{\sigma}_{\text{eq}}^{e}}\cdot \frac{\partial \mathit{\Delta p}}{\partial p}\cdot {\tilde{\sigma}}^{e}\)
\(\frac{\partial \text{Tr}(\mathit{\Delta \sigma })}{\partial p}=-9\cdot K\cdot {\rm A}\cdot \frac{\partial \mathit{\Delta p}}{\partial p}\)
\(\mathit{\Delta p}{,}_{\alpha}\)
On utilise le fait que: \((\sigma +\mathit{\Delta \sigma }{)}_{\text{eq}}=(\sigma +\mathit{\Delta \sigma }{)}_{\text{eq}}^{e}-\mathrm{3\mu }\cdot \mathit{\Delta p}\)
\({\mathit{\Delta p}}_{,\alpha }=\frac{1}{\mathrm{3\mu }}\cdot ((\sigma +\mathit{\Delta \sigma }{)}_{\text{eq}}^{{e}_{,\alpha }}-(\sigma +\mathit{\Delta \sigma }{)}_{{\text{eq}}_{,\alpha }}-\frac{\partial \mathrm{3\mu }}{\partial \alpha }\cdot \mathit{\Delta p})\)
On se réfèrera de nouveau à la remarque à la fin du [r7.01.16-cas-elastoplasticite-drucker-prager] pour les grandeurs dont le calcul n’a pas été détaillé ici.
Calcul de la dérivée du déplacement#
Une fois \(\Delta \sigma {,}_{\alpha}{\mid}_{\Delta \mathrm{u}\mathrm{\ne }\Delta \mathrm{u}(\alpha )}\) calculé, on peut constituer le second membre \(R{,}_{\alpha}{\mid}_{\Delta u\ne \Delta u(\alpha )}\) . On résout alors le système (3864) et l’on obtient l’incrément de déplacement dérivé par rapport à \(\alpha\) .
Calcul de la dérivée des autres grandeurs#
Maintenant que l’on dispose de \(\Delta u{,}_{\alpha}\) , on doit calculer la dérivée des autres grandeurs. On sépare encore deux cas:
Élasticité linéaire
D’après (3867), on calcule comme suit la dérivée de l’incrément de contrainte:
\(\Delta \sigma {,}_{\alpha}=0.+\mathrm{2\mu }.\tilde{\epsilon}(\Delta u{,}_{\alpha})+K.\text{Tr}(\epsilon (\Delta u{,}_{\alpha})).\text{Id}\)
L’incrément de déformation plastique cumulée, quant à lui, ne voit pas d’évolution:
\(\mathrm{\Delta p}{,}_{\alpha}=0\)
É lastoplasticité de type Drucker Prager
Si \(\mathrm{\Delta p}=0\) , on retrouve le cas précédent.
Sinon, on obtient:
\(\mathrm{\Delta \sigma }{,}_{\alpha}=\mathrm{\Delta \sigma }{,}_{\alpha}{\mid}_{\mathrm{\Delta u}\ne \mathrm{\Delta u}(\alpha )}+\frac{\partial \mathrm{\Delta \sigma }}{\partial \epsilon (\mathrm{\Delta u})}:\epsilon (\mathrm{\Delta u}{,}_{\alpha})\)
Et pour la déformation plastique cumulée:
\({\mathrm{\Delta p}}_{,\alpha }=\frac{1}{\mathrm{3\mu }}\cdot ((\sigma +\mathrm{\Delta \sigma }{)}_{\text{eq}}^{{e}_{,\alpha }}-(\sigma +\mathrm{\Delta \sigma }{)}_{{\text{eq}}_{,\alpha }}-\frac{\partial \mathrm{3\mu }}{\partial \alpha }\cdot \mathrm{\Delta p})\)
Une fois que tous ces calculs sont terminés, on réactualise toutes les grandeurs dérivées et on passe au pas de charge suivant.
Synthèse#
Pour résumer les paragraphes précédents, on représente les différentes étapes du calcul par le diagramme suivant:
Assemblage de \(R{,}_{\alpha}{\mid}_{\Delta u\ne \Delta u(\alpha )}\)
Calcul des termes \(\Delta \sigma {,}_{\alpha}{\mid}_{\Delta u\ne \Delta u(\alpha )}\)
Convergence du pas de charge \(n\) du calcul direct
Passage au pas
de charge \(n+1\)
Calcul de
\(\Delta \sigma {,}_{\alpha}\) et \(\mathrm{\Delta p}{,}_{\alpha}\)
Résolution du système (3864)
\(\Rightarrow \Delta u{,}_{\alpha}\)
Fonctionnalités et vérification#
La loi de comportement peut être définie par les mots-clés DRUCK_PRAG et DRUCK_PRAG_N_A pour la version non-associée (commande STAT_NON_LINE, mot clé facteur COMPORTEMENT). Elles sont associées aux matériaux DRUCK_PRAG et DRUCK_PRAG_FO (commande DEFI_MATERIAU).
La loi HOEK_BROWN est vérifiée par les cas tests suivants:
SSND104 |
Validation du comportement DRUCK_PRAG_N_A |
|
SSNP124 |
Essai biaxial drainé avec un comportement DRUCK_PRAGERadoucissant |
|
SSNP125 |
Documentation inexistante |
Validation de l’option INDL_ELGApour le comportement DRUCK_PRAGER |
SSNV168 |
Essai triaxial drainé avec un comportement DRUCK_PRAGERadoucissant |
|
WTNA101 |
Essai triaxial non-drainé avec un comportement DRUCK_PRAGERadoucissant |
|
WTNP114 |
Cas test de référence pour le calcul des déformations mécaniques |
|
SENSM12 |
[v1.01.190] |
Plaque sous pression en déformations planes (plasticité de DRUCK_PRAGER) |
SENSM13 |
[v1.01.192] |
Essai tri-axial avec le modèle de type 3D |
SENSM14 |
[v1.01.193] |
Cavité 2D calcul de sensibilité (Loi DRUCK_PRAGER) |
Bibliographie#
BENALLAL A. et COMI C.: The role of deviatoric and volumetric non-associativities in strain localisation (1993).
CANO V : Instabilités et rupture dans les solides élasto-visco-plastiques (1996).
RICE JR et RUDNICKI JW: A note on some features of the theory of localization of deformation (1980).
RICE JR : The localization of plastic deformations, in Theoretical and Applied Mechanics (1976).
HILL R : A general theory of uniqueness and stability in elastic-plastic solids (1958).
ORTIZ M : An analytical study of the localized failure modes of concrete (1987).
DOGHRI I : Etude de la localisation de l’endommagement (1989).
Annexe 1 Calcul des dérivées partielles de \(\Delta p\)#
Calcul de la dérivée partielle de l’incrément de déformation plastique dans le cas d’un écrouissage linéaire#
\(R(p)=h\cdot p+{\sigma}^{y}\) pour \(0\le p<{p}_{\text{ultm}}\)
\(\mathrm{\Delta p}=\frac{{\sigma}_{\text{eq}}^{e}+{\rm A}\cdot \text{Tr}({\sigma}^{e})-h\cdot {p}^{-}-{\sigma}^{Y}}{\mathrm{9K}\cdot {{\rm A}}^{2}+\mathrm{3\mu }+h}\)
donc:
\(\begin{array}{}\frac{\partial \mathrm{\Delta p}}{\partial \Phi }=\frac{1}{\mathrm{9K}.{{\rm A}}^{2}+\mathrm{3\mu }+h}\cdot (\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \Phi }+\frac{\partial {\rm A}}{\partial \Phi }\cdot \text{Tr}({\sigma}^{e})+{\rm A}\cdot \frac{\partial \text{Tr}({\sigma}^{e})}{\partial \Phi }-\frac{\partial h}{\partial \Phi }\cdot {p}^{-}-\frac{\partial {\sigma}^{y}}{\partial \Phi }\\ -\mathrm{\Delta p}\cdot (9\cdot \frac{\partial K}{\partial \Phi }\cdot {{\rm A}}^{2}+18\cdot K\cdot A\cdot \frac{\partial {\rm A}}{\partial \Phi }+\frac{\partial \mathrm{3\mu }}{\partial \Phi }+\frac{\partial h}{\partial \Phi }))\end{array}\)
\(\frac{\partial \mathrm{\Delta p}}{\partial \sigma }=\frac{1}{\mathrm{3\mu }+\mathrm{9K}\cdot {{\rm A}}^{2}+h}\cdot (A\cdot \frac{\partial \text{Tr}({\sigma}^{e})}{\partial \sigma }+\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \sigma })\)
\(\frac{\partial \mathrm{\Delta p}}{\partial p}=-h\cdot \frac{1}{\mathrm{3\mu }+\mathrm{9K}\cdot {{\rm A}}^{2}+h}\)
\(R(p)=h\cdot {p}_{\text{ultm}}+{\sigma}^{y}\) pour \(p>{p}_{\text{ultm}}\)
\(\mathrm{\Delta p}=\frac{{\sigma}_{\text{eq}}^{e}+{\rm A}\cdot \text{Tr}({\sigma}^{e})-h\cdot {p}_{\text{ultm}}-{\sigma}^{Y}}{\mathrm{9K}\cdot {{\rm A}}^{2}+\mathrm{3\mu }}\)
donc:
\(\begin{array}{}\frac{\partial \mathrm{\Delta p}}{\partial \Phi }=\frac{1}{\mathrm{9K}\cdot {{\rm A}}^{2}+\mathrm{3\mu }}\cdot (\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \Phi }+\frac{\partial {\rm A}}{\partial \Phi }\cdot \text{Tr}({\sigma}^{e})+{\rm A}\cdot \frac{\partial \text{Tr}({\sigma}^{e})}{\partial \Phi }-\frac{\partial h}{\partial \Phi }\cdot {p}_{\text{ultm}}-h\cdot \frac{\partial {p}_{\text{ultm}}}{\partial \Phi }-\frac{\partial {\sigma}^{y}}{\partial \Phi }\\ -\mathrm{\Delta p}\cdot (9\cdot \frac{\partial K}{\partial \Phi }\cdot {{\rm A}}^{2}+18\cdot K\cdot A\cdot \frac{\partial {\rm A}}{\partial \Phi }+\frac{\partial \mathrm{3\mu }}{\partial \Phi }))\end{array}\)
\(\frac{\partial \mathrm{\Delta p}}{\partial \sigma }=\frac{1}{\mathrm{3\mu }+\mathrm{9K}\cdot {{\rm A}}^{2}}(A\cdot \frac{\partial \text{Tr}({\sigma}^{e})}{\partial \sigma }+\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \sigma })\)
\(\frac{\partial \mathrm{\Delta p}}{\partial p}=0\)
Calcul de la dérivée partielle de l’incrément de déformation plastique dans le cas d’un écrouissage parabolique#
\(R(p)={\sigma}^{y}\cdot (1-(1-\sqrt{\frac{{\sigma}_{\text{ultm}}^{y}}{{\sigma}^{y}}})\cdot \frac{p}{{p}_{\text{ultm}}}{)}^{2}\) pour \(0\le p<{p}_{\text{ultm}}\)
\(\begin{array}{}\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \Phi }-(\frac{\partial \mathrm{3\mu }}{\partial \Phi }+{\mathrm{9{\rm A}}}^{2}.\frac{\partial K}{\partial \Phi }+18K.{\rm A}.\frac{\partial {\rm A}}{\partial \Phi }).\mathrm{\Delta p}-(\mathrm{3\mu }+\mathrm{9K}.{{\rm A}}^{2}).\frac{\partial \mathrm{\Delta p}}{\partial \Phi }+\frac{\partial {\rm A}}{\partial \Phi }.\text{Tr}({\sigma}^{e})+{\rm A}.\frac{\partial \text{Tr}({\sigma}^{e})}{\partial \Phi }\\ -\frac{\partial {\sigma}^{y}}{\partial \Phi }.(1-(1-\sqrt{\frac{{\sigma}_{\text{ultm}}^{y}}{{\sigma}^{y}}}).\frac{{p}^{-}+\mathrm{\Delta p}}{{p}_{\text{ultm}}}{)}^{2}\\ -{\mathrm{2\sigma }}^{y}.(1-(1-\sqrt{\frac{{\sigma}_{\text{ultm}}^{y}}{{\sigma}^{y}}}).\frac{{p}^{-}+\mathrm{\Delta p}}{{p}_{\text{ultm}}}).\\ (\frac{\partial {\sigma}_{\text{ultm}}^{y}}{\partial \Phi }.\frac{{p}^{-}+\mathrm{\Delta p}}{{\mathrm{2p}}_{\text{ultm}}.\sqrt{{\sigma}_{\text{ultm}}^{y}.{\sigma}^{y}}}-\frac{\partial {\sigma}^{y}}{\partial \Phi }.\frac{{p}^{-}+\mathrm{\Delta p}}{{\mathrm{2p}}_{\text{ultm}}.{\sigma}^{y}}.\sqrt{\frac{{\sigma}_{\text{ultm}}^{y}}{{\sigma}^{y}}}+\frac{\partial {p}_{\text{ultm}}}{\partial \Phi }.(1-\sqrt{\frac{{\sigma}_{\text{ultm}}^{y}}{{\sigma}^{y}}}).\frac{{p}^{-}+\mathrm{\Delta p}}{{p}_{{\text{ultm}}^{2}}}-\frac{\partial \mathrm{\Delta p}}{\partial \Phi }\frac{1-\sqrt{\frac{{\sigma}_{\text{ultm}}^{y}}{{\sigma}^{y}}}}{{p}_{\text{ultm}}})\\ 0\end{array}\)
\(\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \sigma }-(\mathrm{3\mu }+\mathrm{9K}\cdot {{\rm A}}^{2})\cdot \frac{\partial \mathrm{\Delta p}}{\partial \sigma }+A\cdot \frac{\partial \text{Tr}({\sigma}^{e})}{\partial \sigma }+{\mathrm{2\sigma }}^{y}\cdot (1-(1-\sqrt{\frac{{\sigma}_{\text{ultm}}^{y}}{{\sigma}^{y}}})\cdot \frac{{p}^{-}+\mathrm{\Delta p}}{{p}_{\text{ultm}}})\cdot \frac{1-\sqrt{\frac{{\sigma}_{\text{ultm}}^{y}}{{\sigma}^{y}}}}{{p}_{\text{ultm}}}\cdot \frac{\partial \mathrm{\Delta p}}{\partial \sigma }=0\)
\(-(\mathrm{3\mu }+\mathrm{9K}\cdot {{\rm A}}^{2})\cdot \frac{\partial \mathrm{\Delta p}}{\partial p}+{\mathrm{2\sigma }}^{y}\cdot (1-(1-\sqrt{\frac{{\sigma}_{\text{ultm}}^{y}}{{\sigma}^{y}}})\frac{{p}^{-}+\mathrm{\Delta p}}{{p}_{\text{ultm}}})\cdot \frac{1-\sqrt{\frac{{\sigma}_{\text{ultm}}^{y}}{{\sigma}^{y}}}}{{p}_{\text{ultm}}}\cdot (1+\frac{\partial \mathrm{\Delta p}}{\partial p})=0\)
\(R(p)={\sigma}_{\text{ultm}}^{y}\) pour \(p>{p}_{\text{ultm}}\)
\(\frac{\partial \mathrm{\Delta p}}{\partial \Phi }=\frac{1}{\mathrm{3\mu }+\mathrm{9K}\cdot {{\rm A}}^{2}}(\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \Phi }-(\frac{\partial \mathrm{3\mu }}{\partial \Phi }+\frac{\partial \mathrm{9K}}{\partial \Phi }\cdot {{\rm A}}^{2}+18K\cdot \frac{\partial {\rm A}}{\partial \Phi }\cdot A)\cdot \mathrm{\Delta p}+\frac{\partial {\rm A}}{\partial \Phi }\cdot \text{Tr}({\sigma}^{e})+{\rm A}\cdot \frac{\partial \text{Tr}({\sigma}^{e})}{\partial \Phi }-\frac{\partial {\sigma}_{\text{ultm}}^{y}}{\partial \Phi })\)
\(\frac{\partial \mathrm{\Delta p}}{\partial \sigma }=\frac{1}{\mathrm{3\mu }+\mathrm{9K}\cdot {{\rm A}}^{2}}(\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \sigma }+A\cdot \frac{\partial \text{Tr}({\sigma}^{e})}{\partial \sigma })\)
\(\frac{\partial \mathrm{\Delta p}}{\partial p}=0\)
Cas de la projection au sommet du cône#
Le principe de la résolution analytique consiste à déterminer les contraintes effectives comme la projection des contraintes élastiques sur le critère.
Il se peut qu’il n’y ait pas de solution.
Si la condition \(\mathrm{\Delta p}\le \frac{{\sigma}_{\text{eq}}^{e}}{\mathrm{3\mu }}\) n’est pas respectée, il faut trouver les contraintes effectives par projection au sommet du cône \(\mathrm{\Delta p}=\frac{{\sigma}_{\text{eq}}^{e}}{\mathrm{3\mu }}\) .
Dans ce cas, on obtient:
\(\frac{\partial \mathrm{\Delta p}}{\partial \Phi }=\frac{1}{\mathrm{3\mu }}\cdot (\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \Phi }-\mathrm{\Delta p}\cdot \frac{\partial \mathrm{3\mu }}{\partial \Phi })\)
\(\frac{\partial \mathrm{\Delta p}}{\partial \sigma }=\frac{1}{\mathrm{3\mu }}\cdot \frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \sigma }\)
\(\frac{\partial \mathrm{\Delta p}}{\partial p}=0\)