r5.03.06 Modèle de Rousselier en grandes déformations

Résumé

On présente ici une variante du modèle de Rousselier qui permet de décrire la croissance plastique de cavités dans un acier. La relation de comportement est élastoplastique avec écrouissage isotrope, permet les changements de volume plastique et est écrite en grandes déformations. Ces dernières s’appuient sur la théorie proposée par Simo et Miehe, modifiée pour faciliter l’intégration numérique de la loi de comportement et replacer le modèle dans le cadre des matériaux standard généralisés.

Ce modèle est disponible dans la commande STAT_NON_LINE par l’intermédiaire du mot-clé RELATION:’ROUSSELIER’ sous le mot-clé facteur COMPORTEMENT et avec le mot-clé DEFORMATION:’SIMO_MIEHE’.

Ce modèle est implanté pour les modélisations tridimensionnelles (3D), axisymétrique (AXIS) et en déformations planes (D_PLAN).

On présente l’écriture et le traitement numérique de ce modèle.

Table des matières

Notations

On notera par:

\(\mathrm{Id}\)	matrice identité
\(\text{tr}A\)	trace du tenseur A
\({A}^{T}\)	transposé du tenseur A
\(\detA\)	déterminant de A
\(\tilde{A}\)	partie déviatorique du tenseur A définie par \(\tilde{A}=A-(\frac{1}{3}\text{tr}A)\text{Id}\)
\({A}_{H}\)	partie hydrostatique du tenseur A définie par \({A}_{H}=\frac{\text{tr}A}{3}\)
:	produit doublement contracté: \(A:B=\sum_{i,j}{A}_{ij}{B}_{ij}=\text{tr}({\text{AB}}^{T})\)
\(\otimes\)	produit tensoriel: \((A\otimes B{)}_{ijkl}={A}_{ij}{B}_{kl}\)
\({A}_{\text{eq}}\)	valeur équivalente de Von Mises définie par \({A}_{\text{eq}}=\sqrt{\frac{3}{2}\tilde{A}:\tilde{A}}\)
\({\nabla}_{X}A\)	gradient: \({\nabla}_{X}A=\frac{\partial A}{\partial X}\)
\(\lambda ,\mu ,E,\nu ,K\)	coefficients de l’élasticité isotrope
\({\sigma}_{y}\)	limite d’élasticité
\(\alpha\)	coefficient de dilatation thermique
\(T\)	température
\({T}_{\text{ref}}\)	température de référence

Par ailleurs, dans le cadre d’une discrétisation en temps, toutes les quantités évaluées à l’instant précédent sont indicées par \({}^{-}\) , les quantités évaluées à l’instant \(t+\Delta t\) ne sont pas indicées et les incréments sont désignés par \(\Delta\) . On a ainsi :

\(\Delta Q=Q-{Q}^{-}\)

Théorie de Simo et Miehe

Introduction

Nous rappelons ici les spécificités de la formulation proposée par SIMO J.C. et MIEHE C. [bib3] pour traiter les grandes déformations. Cette formulation a déjà été utilisée pour des modèles de comportement thermo-élasto-(visco)-plastique avec écrouissage isotrope et critère de Von Mises, [R5.03.21] pour un modèle sans effet des transformations métallurgiques et [R4.04.03] pour un modèle avec effet des transformations métallurgiques.

Les choix cinématiques permettent de traiter des grands déplacements et des grandes déformations mais également des grandes rotations de manière exacte.

Les spécificités de ces modèles sont les suivantes:

• tout comme en petites déformations, on suppose l’existence d’une configuration relâchée, c’est-à-dire localement libre de contrainte, qui permet de décomposer la déformation totale en une partie thermoélastique et une partie plastique,

• la décomposition de cette déformation en des parties thermoélastique et plastique n’est plus additive comme en petites déformations (ou pour les modèles grandes déformations écrits en taux de déformation avec par exemple une dérivée de Jaumann) mais multiplicative,

• les déformations élastiques sont mesurées dans la configuration actuelle (déformée) tandis que les déformations plastiques sont mesurées dans la configuration initiale,

• comme en petites déformations, les contraintes dépendent uniquement des déformations thermo-élastiques,

• si le critère de plasticité ne dépend que de la contrainte déviatorique, alors les déformations plastiques se font à volume constant. La variation de volume est alors uniquement due aux déformations thermo-élastiques,

• ce modèle conduit lors de son intégration numérique à un modèle incrémentalement objectif (cf. [§3.2.3]) ce qui permet d’obtenir la solution exacte en présence de grandes rotations.

Par la suite, on rappelle brièvement quelques notions de mécanique en grandes déformations.

Généralités sur les grandes déformations

Cinématique

Considérons un solide soumis à des grandes déformations. Soit \({\Omega}_{0}\) le domaine occupé par le solide avant déformation et \(\Omega (t)\) le domaine occupé à l’instant t par le solide déformé.

Figure 3.2.1-a: Représentation de la configuration initiale et déformée

Dans la configuration initiale \({\Omega}_{0}\) , la position de toute particule du solide est désignée par \(X\) (description lagrangienne). Après déformation, la position à l’instant \(t\) de la particule qui occupait la position \(X\) avant déformation est donnée par la variable \(x\) (description eulérienne).

Le mouvement global du solide est défini, avec \(u\) le déplacement, par:

\(x=\stackrel{ˆ}{x}(X,t)=X+u\) éq 3.2.1-1

Pour définir le changement de métrique au voisinage d’un point, on introduit le tenseur gradient de la transformation \(F\) :

\(F=\frac{\partial \stackrel{ˆ}{x}}{\partial X}=\text{Id}+{\nabla}_{X}u\) éq 3.2.1-2

Les transformations de l’élément de volume et de la masse volumique valent:

\(d\Omega =\text{Jd}{\Omega}_{o}\) avec \(J=\detF=\frac{{\rho}_{o}}{\rho}\) éq 3.2.1-3

où \({\rho}_{o}\) et \(\rho\) sont respectivement la masse volumique dans les configurations initiale et actuelle.

Différents tenseurs de déformations peuvent être obtenus en éliminant la rotation dans la transformation locale. Par exemple, en calculant directement les variations de longueur et d’angle (variation du produit scalaire), on obtient:

\(E=\frac{1}{2}(C-\text{Id})\) avec \(C={F}^{T}F\) éq 3.2.1-4

\(A=\frac{1}{2}(\text{Id}-{b}^{-1})\) avec \(b={\text{FF}}^{T}\) éq 3.2.1-5

\(E\) et \(A\) sont respectivement les tenseurs de déformation de Green-Lagrange et d’Euler-Almansi et \(C\) et \(b\) , les tenseurs de Cauchy-Green droit et gauche respectivement.

En description lagrangienne, on décrira la déformation par les tenseurs \(C\) ou \(E\) car ce sont des quantités définies sur \({\Omega}_{0}\) , et en description eulérienne par les tenseurs \(b\) ou \(A\) (définis sur \(\Omega\) ).

Contraintes

Le tenseur des contraintes utilisé dans la théorie de Simo et Miehe est le tenseur de Kirchhoff \(\tau\) défini par:

\(J\sigma =\tau\) éq 3.2.2-1

où \(\sigma\) est le tenseur eulérien de Cauchy. Le tenseur \(\tau\) résulte donc d’une «mise à l’échelle» par la variation de volume du tenseur de Cauchy \(\sigma\) .

Objectivité

Lorsqu’on écrit une loi de comportement en grandes déformations, on doit vérifier que cette loi est objective, c’est-à-dire invariante par tout changement de référentiel spatial de la forme :

\({x}^{\text{*}}=c(t)+Q(t)x\) éq 3.2.3-1

où \(Q\) est un tenseur orthogonal qui traduit la rotation du référentiel et \(c\) un vecteur qui traduit la translation.

Plus concrètement, si on réalise un essai de traction dans la direction \({e}_{1}\) , par exemple, suivi d’une rotation de 90° autour de \({e}_{3}\) , ce qui revient à effectuer un essai de traction selon \({e}_{2}\) , alors le danger avec une loi de comportement non objective est de ne pas retrouver un tenseur des contraintes uniaxial dans la direction \({e}_{2}\) (ce qui est notamment le cas avec la cinématique PETIT_REAC).

Formulation de Simo et Miehe

Par la suite, on notera \(F\) le tenseur gradient qui fait passer de la configuration initiale \({\Omega}_{0}\) à la configuration actuelle \(\Omega (t)\) , \({F}^{p}\) le tenseur gradient qui fait passer de la configuration \({\Omega}_{0}\) à la configuration relâchée \({\Omega}^{r}\) , et \({F}^{e}\) de la configuration \({\Omega}^{r}\) à \(\Omega (t)\) . L’indice \(p\) se réfère à la partie plastique, l’indice \(e\) à la partie élastique.

Figure 3.3-a: Décomposition du tenseur gradient \(F\) en une partie élastique \({F}^{e}\) et plastique \({F}^{p}\)

Par composition des mouvements, on obtient la décomposition multiplicative suivante:

\(F={F}^{e}{F}^{p}\) éq 3.3-1

Les déformations élastiques sont mesurées dans la configuration actuelle avec le tenseur eulérien de Cauchy-Green gauche \({b}^{e}\) et les déformations plastiques dans la configuration initiale par le tenseur \({G}^{p}\) (description lagrangienne). Ces deux tenseurs sont définis par:

\({b}^{e}={F}^{e}{F}^{\text{eT}}\) , \({G}^{p}=({F}^{\text{pT}}{F}^{p}{)}^{-1}\) d’où \({b}^{e}={\text{FG}}^{p}{F}^{T}\) éq 3.3-2

Toutefois, on emploiera alternativement une autre mesure des déformations élastiques \(e\) qui coïncide avec l’opposé des déformations linéarisées lorsque les déformations élastiques sont petites:

\(e=\frac{1}{2}(\text{Id}-{b}^{e})\) éq 3.3-3

Dans le cas d’un matériau isotrope, on peut montrer que le potentiel énergie libre ne dépend que du tenseur de Cauchy-Green gauche \({b}^{e}\) (où dans notre cas du tenseur \(e\) ) et en plasticité de la variable \(p\) liée à l’écrouissage isotrope. De plus, on suppose que l’énergie libre volumique se décompose, tout comme en petites déformations, en une partie hyperélastique qui ne dépend que de la déformation élastique et une autre liée au mécanisme d’écrouissage:

\(\Phi (e,p)={\Phi}^{\text{el}}(e)+{\Phi}^{\text{bl}}(p)\) éq 3.3-4

Si au lieu d’utiliser la contrainte de Cauchy \(\sigma\) , on utilise la contrainte de Kirchhoff \(\tau\) , l’inégalité de Clausius-Duhem s’écrit (on oublie la partie thermique):

\(\tau :D-\dot{\Phi}\ge 0\) éq 3.3-5

expression dans laquelle \(D\) représente le taux de déformation eulérien.

Sous les hypothèses précédentes, la dissipation s’écrit encore:

\((\tau +\frac{\partial \Phi }{\partial e}{b}^{e}):D+\frac{1}{2}\frac{\partial \Phi }{\partial e}:(F{\dot{G}}^{p}{F}^{T})-\frac{\partial \Phi }{\partial p}\dot{p}\ge 0\) éq 3.3-6

Le second principe de la thermodynamique requiert alors l’expression suivante pour la relation contrainte-déformation:

\(\tau =-\frac{\partial \Phi }{\partial e}{b}^{e}\) éq 3.3-7

On définit enfin les forces thermodynamiques associées à la déformation élastique et à la déformation plastique cumulée conformément au cadre des matériaux standard généralisés:

\(s=-\frac{\partial \Phi }{\partial e}\text{}\text{soit}\text{}\tau =s{b}^{e}\) éq 3.3-8

\(A=-\frac{\partial \Phi }{\partial p}\) éq 3.3-9

où la force thermodynamique A est l’opposé de la variable d’écrouissage isotrope \(R\) .

Il reste alors pour la dissipation:

\(\tau :(-\frac{1}{2}F{\dot{G}}^{p}{F}^{T}{b}^{e-1})+A\dot{p}=s:(-\frac{1}{2}F{\dot{G}}^{p}{F}^{T})+A\dot{p}\ge 0\) éq 3.3-10

Formulation originelle

Le principe de dissipation maximale appliqué à partir du seuil d’élasticité \(F\) , fonction de la contrainte de Kirchhoff \(\tau\) et de la force thermodynamique \(A\) permet d’en déduire les lois d’évolution de la déformation plastique et de la déformation plastique cumulée, soit:

\(-\frac{1}{2}F{\dot{G}}^{p}{F}^{T}{b}^{e\text{-1}}=\dot{\lambda}\frac{\partial F}{\partial \tau }\) éq 3.3.1-1

\(\dot{p}=\dot{\lambda}\frac{\partial F}{\partial A}\) éq 3.3.1-2

\(\dot{\lambda}\ge 0\text{}F\le 0\text{}F\dot{\lambda}=0\) éq 3.3.1-3

Formulation modifiée

L’approximation introduite ici sur la formulation originelle de Simo et Miehe porte sur l’expression de la loi d’écoulement, approximation d’autant plus réduite que les déformations élastiques sont petites, puisque \(\tau =s{b}^{e}\) . En effet, on exprime dorénavant le seuil d’élasticité comme une fonction des forces thermodynamiques et non plus des contraintes \(F(s,A)\le 0\) , et c’est par rapport à ces variables qu’on applique le principe de dissipation maximale, ce qui conduit aux lois d’écoulement suivantes:

\(-\frac{1}{2}F{\dot{G}}^{p}{F}^{T}=\dot{\lambda}\frac{\partial F}{\partial s}\) éq 3.3.2-1

\(\dot{p}=\dot{\lambda}\frac{\partial F}{\partial A}\) éq 3.3.2-2

\(\dot{\lambda}\ge 0\text{}F\le 0\text{}F\dot{\lambda}=0\) éq 3.3.2-3

Conséquences de l’approximation

En remplaçant la contrainte \(\tau\) par la force thermodynamique \(s\) associée à la déformation élastique dans l’expression du critère de plasticité, on introduit en fait une perturbation de la frontière du domaine de réversibilité de l’ordre de grandeur de \(2\parallel e\parallel\) . Par rapport à la formulation initiale, il en résulte évidemment une influence sur la limite d’élasticité observée mais aussi sur la direction d’écoulement: en particulier, la dérivée par rapport au temps de la variation de volume plastique s’écrit alors:

\({\dot{J}}^{p}=\dot{\lambda}{J}^{p}{b}^{e\text{-1}}:\frac{\partial F}{\partial s}\) éq 3.3.3-1

si bien que dans le cas où le critère \(F\) ne dépend que du déviateur du tenseur des contraintes \(s\) , on ne retrouve pas \({J}^{p}=1\) : le caractère isochore de la déformation plastique n’est plus parfaitement préservé. Nous serons alors amené à introduire une correction de volume a posteriori.

Dans la mesure où les déformations élastiques restent petites, les résultats obtenus avec ce modèle modifié ne s’écartent pas significativement de ceux obtenus avec l’ancienne formulation (cf. [bib4]), tandis que l’intégration numérique en sera simplifiée. En effet, on verra par la suite que ce modèle suit le même schéma d’intégration que celui des modèles écrits en petites déformations.

Remarque:

Cette nouvelle formulation des grandes déformations permet de replacer la théorie de Simo et Miehe dans le cadre des matériaux standard généralisés. D’un point de vue numérique, ceci a pour conséquence d’exprimer la résolution de la loi de comportement comme un problème de minimisation par rapport aux incréments de variables internes.

En effet, on rappelle que dans le cadre des matériaux standard généralisés, la donnée des deux potentiels l’énergie libre \(\Phi (\varepsilon ,a)\) et le potentiel de dissipation \(D(\dot{a})\) , fonction du tenseur de déformation \(\varepsilon\) et d’un certain nombre de variables internes \(a\) , permet de définir complètement la loi de comportement (on se place dans le cas des matériaux indépendant du temps).

\(\sigma =\frac{\partial \Phi }{\partial \varepsilon }\) , \(A=-\frac{\partial \Phi }{\partial a}\in \partial D(\dot{a})\) éq 3.3.3-2

où \(\partial D(\dot{a})\) est le sous différentiel du potentiel de dissipation \(D\) .

Les lois de comportement de type standard généralisé qui ne dépendent pas du temps sont caractérisées par un potentiel de dissipation positivement homogène de degré 1, qui se traduit par la propriété suivante:

\(\forall \dot{a}\text{}\forall \lambda >0\text{}D(\lambda \dot{a})=\lambda D(\dot{a})\text{}\Rightarrow \text{}\partial D(\lambda \dot{a})=\partial D(\dot{a})\) éq 3.3.3-3

Maintenant si on écrit le problème [éq 3.3.3-2] sous forme discrétisée en temps et si on utilise la propriété des sous différentiels [éq 3.3.3-3], on obtient le problème discrétisé suivant:

\(\sigma =\frac{\partial \Phi }{\partial \varepsilon }\) , \(A=-\frac{\partial \Phi }{\partial a}\in \partial D(\Delta a)\) éq 3.3.3-4

On peut montrer que l’équation [éq 3.3.3-4] est équivalente (cf. [bib5]) à résoudre le problème de minimisation par rapport aux incréments de variables internes \(\Delta a\) suivant:

\(-\frac{\partial \Phi }{\partial a}\in \partial D(\Delta a)\iff \Delta a=\text{Arg}\underset{\Delta {a}^{\text{*}}}{\text{Min}}\left[\Phi ({a}^{-}+{\mathrm{\Delta a}}^{\text{*}})+D(\Delta {a}^{\text{*}})\right]\) éq 3.3.3-5

L’application de l’équation [éq 3.3.3-5] au modèle de Rousselier en grandes déformations modifiées s’écrit:

\(\underset{\text{énergie continue}}{\Phi (e,p)\text{et}D({D}^{p},\dot{p})}\text{}\underset{\text{discrétisation}}{\text{=>}}\text{}\underset{\text{énergie discrétisée}}{\Phi ({e}^{\text{Tr}}+\Delta e,{p}^{-}+\Delta p)\text{et}D(\Delta e,\Delta p)}\) éq 3.3.3-6

\(\begin{array}{}A=-\frac{\partial \Phi }{\partial a}=\lbrace \begin{array}{}s=-\frac{\partial \Phi }{\partial e}\\ -R=-\frac{\partial \Phi }{\partial p}\end{array}\text{}\in \partial D(\Delta e,\Delta p)\\ \iff \underset{\Delta e,\Delta p}{\text{Min}}\left[\Phi ({e}^{\text{Tr}}+\Delta e,{p}^{-}+\Delta p)+D(\Delta e,\Delta p)\right]\text{}\end{array}\) éq. 3.3.3-7

On trouvera dans le paragraphe [§4], la relation qui lie le taux de déformation plastique \({D}^{p}\) une fois discrétisée et l’incrément de déformation élastique \(\Delta e\) , ainsi que la définition de \({e}^{\text{Tr}}\) .

On voit bien ici que si on prend la formulation initiale de Simo et Miehe, on ne peut plus écrire le problème de minimisation [éq 3.3.3-7] avec la contrainte de Kirchhoff \(\tau\) à cause du terme en \({b}^{e}\) dans l’expression:

\(\tau =-\frac{\partial \Phi }{\partial e}{b}^{e}\) éq 3.3.3-8

Modèle de Rousselier

Nous décrivons maintenant l’application des grandes déformations au modèle de Rousselier présenté en introduction.

Equations du modèle

Pour décrire un modèle thermoélastoplastique à écrouissage isotrope (l’équivalent en petites déformations au modèle à écrouissage isotrope et critère de Von Mises), Simo et Miehe proposent un potentiel élastique polyconvexe. Par raison de simplicité, on choisit ici le potentiel de Saint Venant qui s’écrit:

\(\Phi (e,p)=\frac{1}{2}\left[K{(\text{tr}e)}^{2}+2\mu \tilde{e}:\tilde{e}+\mathrm{6K}\alpha \Delta T\text{tr}e\right]+\underset{0}{\overset{p}{\int}}R(u)\text{du}\) éq 4.1-1

Conformément aux équations [éq 3.3-8] et [éq 3.3-9], les lois d’état qui dérivent du potentiel élastique ci-dessus s’écrivent alors:

\(s=-\left[K\text{tr}e\text{Id}+2\mu \tilde{e}+\mathrm{3K}\alpha \Delta T\text{Id}\right]\) éq 4.1-2

\(A=-R(p)\) éq 4.1-3

Le seuil d’élasticité est donné par:

\(F(s,R)={s}_{\text{eq}}+{\sigma}_{1}\text{Df}\exp(\frac{{s}_{H}}{{\sigma}_{1}})-R-{\sigma}_{y}\) éq 4.1-4

D’après les équations [éq 3.3.2-1] et [éq 3.3.2-2], les lois d’écoulement sont définies par:

\(-\frac{1}{2}F{\dot{G}}^{p}{F}^{T}=\dot{\lambda}\left[\frac{3\tilde{s}}{2{s}_{\text{eq}}}+\frac{\text{Df}}{3}\exp(\frac{{s}_{H}}{{\sigma}_{1}})\text{Id}\right]\) éq 4.1-5

\(\dot{p}=\dot{\lambda}\) éq 4.1-6

\(\dot{\lambda}\ge 0\text{}F\le 0\text{}F\dot{\lambda}=0\) éq 4.1-7

Traitement des points singuliers

En fait, l’équation d’écoulement [éq 4.1-5] traduit l’appartenance de la direction d’écoulement au cône normal à la surface du domaine d’élasticité. Elle n’est valide qu’aux points réguliers, caractérisés par:

\({s}_{\text{eq}}\ne 0\) éq 4.2-1

Il reste donc à caractériser le cône normal aux points singuliers, c’est-à-dire vérifiant:

\(\tilde{s}=0\text{et}{\sigma}_{1}Df\exp(\frac{{s}_{H}}{{\sigma}_{1}})-R={\sigma}_{y}\) éq 4.2-2

Le cône normal au convexe d’élasticité en un tel point est l’ensemble des directions d’écoulement qui réalisent le problème de maximisation suivant:

\({\Delta}^{\text{*}}(s\text{,R})=\underset{{D}^{p},\dot{p}}{\sup}\left[s:{D}^{p}-R\dot{p}-\Delta ({D}^{p},\dot{p})\right]\) éq 4.2-3

où \({\Delta}^{\text{*}}\) est la fonction indicatrice du convexe \(F\) et \(\Delta ({D}^{p},\dot{p})\) le potentiel de dissipation obtenu par transformée de Legendre-Fenchel de la fonction indicatrice de \(F\) :

\(\Delta ({D}^{p},\dot{p})=\underset{\begin{array}{}s,R\\ F(s,R)\le 0\end{array}}{\text{Sup}}\left[s:{D}^{p}-R\dot{p}\right]\) éq 4.2-4

Après quelques calculs, on obtient:

\(\Delta ({D}^{p},\dot{p})={\sigma}_{y}\dot{p}+{\sigma}_{1}\text{tr}{D}^{p}(\ln\frac{\text{tr}{D}^{p}}{Df\dot{p}}-1)+{I}_{{ℝ}^{+}}(\text{tr}{D}^{p})+{I}_{{ℝ}^{+}}(\dot{p}-\frac{2}{3}{D}_{\text{eq}}^{p})\) éq 4.2-5

avec

\({I}_{{ℝ}^{+}}(x)=\lbrace \begin{array}{}0\text{}\text{si}\text{}x\ge 0\\ +\infty \text{sinon}\end{array}\) éq 4.2-6

Pour \(\tilde{s}=0\) , \({\Delta}^{\text{*}}\) vaut:

\({\Delta}^{\text{*}}(s,R)=\underset{\begin{array}{}{D}^{p},\dot{p}\\ \text{tr}{D}^{p}\ge 0\\ \dot{p}-\frac{2}{3}{D}_{\mathrm{eq}}^{p}\ge 0\end{array}}{\text{Sup}}\left[\underset{G(\text{tr}{D}^{p})}{\underset{\underbrace{}}{{s}_{H}\text{tr}{D}^{p}-{\sigma}_{1}\text{tr}{D}^{p}(\ln\frac{\text{tr}{D}^{p}}{Df\dot{p}}-1)}}-R\dot{p}-{\sigma}_{y}\dot{p}\right]\) éq 4.2-7

En remarquant que pour \(\text{tr}{D}^{p}\ge 0\) , la fonction \(G(\text{tr}{D}^{p})\) est concave, le suprémum par rapport à la trace du taux de déformation plastique \({D}^{p}\) est obtenu pour:

\({G}^{'}(\text{tr}{D}^{p})=0\text{}\text{d'où}\text{}\text{tr}{D}^{p}=\mathrm{Df}\dot{p}\exp(\frac{{s}_{H}}{{\sigma}_{1}})\) éq 4.2-8

Remarque:

On retrouve bien alors pour la fonction indicatrice du seuil d’élasticité \(F\) .

\({D}^{\text{*}}(s\text{,R})=\underset{\begin{array}{}\dot{p}\\ \dot{p}\ge \frac{2}{3}{D}_{\text{eq}}^{p}\ge 0\end{array}}{\text{Sup}}\left[F\dot{p}\right]=\lbrace \begin{array}{cc}0& \text{si}F\le 0\\ +\infty & \text{sinon}\end{array}\) éq 4.2-9

En un point singulier, le cône normal, ensemble des directions d’écoulement admissibles, se caractérise donc par:

\(\text{tr}{D}^{p}=Df\dot{p}\exp(\frac{{s}_{H}}{{\sigma}_{1}})\) éq 4.2-10

\(\dot{p}\ge \frac{2}{3}{D}_{\text{eq}}^{p}\ge 0\) éq 4.2-11

Expression de la porosité

On a vu en introduction que l’inspiration microscopique du modèle de Rousselier se fonde sur une microstructure constituée d’une cavité et d’une matrice rigide plastique, donc isochore. Dans ce cas, la porosité est directement reliée à la déformation macroscopique parla relation eq. 1-1.

Dans cette expression, le changement de volume d’origine élastique est négligé. Sans précaution particulière, cette approximation peut s’avérer pénalisante en présence de compression élastique, même raisonnable, car elle conduit à une porosité éventuellement négative.

On lui préfère donc l’expression équivalente suivante, assortie d’une minoration explicite par la porosité initiale:

\(f=\max({f}_{0,}1-\frac{1-{f}_{0}}{J})\) éq 4.3-1

Rousselier propose quant à lui d’exprimer la porosité en se basant sur le taux de déformation plastique \({D}^{p}\) . La relation est écrite sous forme incrémentale:

\(\dot{f}=(1-f)\text{tr}{D}^{p}\) éq 4.3-2

Cela signifie que la variable porosité employée pour paramétrer le critère de plasticité \(F\) ne dépend que de la déformation plastique. En fait, le taux de déformation plastique est une quantité évaluée dans la configuration relâchée. Son transport dans la configuration actuelle (comme \(D\) ) s’exprime encore:

\({F}^{e}{D}^{p}{F}^{{e}^{T}}=-\frac{1}{2}F{\dot{G}}^{p}{F}^{T}\) éq 4.3-3

Dans ce cas, la loi d’évolution de la porosité s’exprime:

\(\dot{f}=(1-f)\text{tr}(-\frac{1}{2}F{\dot{G}}^{p}{F}^{T})\) éq 4.3-4

Pour éviter que l’intégration de la porosité n’interfère avec celle de la plasticité (puisque les deux variables sont couplées), il est nécessaire de séparer l’intégration de la loi de comportement en deux temps: intégration de la plasticité à porosité fixée à sa valeur au début du pas de temps, puis intégration de la porosité au moyen de l’équation 4.3-4 où l’évolution plastique est celle calculée à la phase précédente.

Remarque importante:

On a ainsi deux versions possibles du modèle (PORO_TYPE = 1 ou 2,cf. U4.43.01), selon que l’on choisisse respectivement la déformation totale ou la déformation plastique dans l’évolution de la porosité. Il a été remarqué que pour une porosité initiale f0 faible, le comportement au début d’évolution change fortement en fonction de ce paramètre. En effet, le choix semble avoir un impact déterminant sur la réponse à l’échelle de la structure, puisqu’il privilégie ou non les bifurcations de la zone où les déformations se localisent. Ainsi, cette forte sensibilité doit conduire l’utilisateur à une grande prudence dans l’utilisation de ce modèle. Des recherches sont en cours pour comprendre cette sensibilité et discriminer les deux variantes pour l’évolution de la porosité.

Relation ‘ROUSSELIER‘

Cette relation de comportement est disponible via l’argument ‘ROUSSELIER‘ du mot clé COMPORTEMENT sous l’opérateur STAT_NON_LINE, avec l’argument ‘SIMO_MIEHE‘ du mot-clé facteur DEFORMATION.

L’ensemble des paramètres du modèle est fourni sous les mots clés facteurs ‘ROUSSELIER‘ ou ‘ROUSSELIER_FO‘ et ‘TRACTION‘ (pour définir la courbe de traction) de la commande DEFI_MATERIAU ([U4.43.01]).

Remarque :

L’utilisateur doit bien s’assurer que la courbe de traction «expérimentale» utilisée, soit directement, soit pour en déduire la pente d’écrouissage est bien donnée dans le plan contrainte rationnelle \(\sigma =F/S\) - déformation logarithmique \(\ln(1+\Delta l/{l}_{0})\) où \({l}_{0}\) est la longueur initiale de la partie utile de l’éprouvette, \(\Delta l\) la variation de longueur après déformation,Fla force appliquée et Sla surface actuelle .

Contraintes et variables internes

Les contraintes sont les contraintes de Cauchy , \(\sigma\) calculées donc sur la configuration actuelle (six composantes en 3D, quatre en 2D).

Les variables internes produites dans le Code_Aster sont:

• V1, la déformation plastique cumulée \(p\) ,

• V2, la porosité \(f\) ,

• V3, l’indicateur de plasticité (0 si le dernier incrément calculé est élastique, 1 si solution plastique régulière, 2 si solution plastique singulière),

• V4 à V9, le tenseur de déformation élastique \(e\) .

Remarque :

Si l’utilisateur veut récupérer éventuellement des déformations en post-traitement de son calcul, il faut tracer les déformations de Green-Lagrange \(E\) , qui représente une mesure des déformations en grandes déformations (option EPSG_ELGA ou EPSG_ELNO de CALC_CHAMP). Les déformations linéarisées \(\varepsilon\) classiques mesurent des déformations sous l’hypothèse des petites déformations et n’ont pas de sens en grandes déformations.

Formulation numérique

Pour la formulation variationnelle, il s’agit de la même que celle donnée dans la note [R5.03.21] et qui se reporte à la loi de comportement avec écrouissage isotrope et critère de Von Mises en grandes déformations. Nous rappelons uniquement qu’il s’agit d’une formulation eulérienne, avec réactualisation de la géométrie à chaque incrément et à chaque itération, et que l’on tient compte de la rigidité de comportement et de la rigidité géométrique.

Nous présentons maintenant l’intégration numérique de la loi de comportement et donnons l’expression de la matrice tangente (options FULL_MECA et RIGI_MECA_TANG).

Expression du modèle discrétisé

Connaissant la contrainte \({\sigma}^{-}\) , la déformation plastique cumulée \({p}^{-}\) , les déformations élastiques \({e}^{-}\) , les déplacements \({u}^{-}\) et \(\Delta u\) , on cherche à déterminer \((\sigma ,p,e)\) .

Les déplacements étant connus, les gradients de la transformation de \({\Omega}_{0}\) à \({\Omega}^{-}\) , noté \({F}^{-}\) , et de \({\Omega}^{-}\) à \(\Omega\) , noté \(\Delta F\) , sont connus.

Pour intégrer ce modèle de comportement, on emploie un schéma d’Euler implicite, la porosité étant une fonction explicite de la déformation via l’équation 4.3-1, donc connue lors de l’intégration du comportement.

Une fois discrétisé, on obtient alors le système suivant:

\(F=\Delta F{F}^{-}\) éq 5.1-1

\(J=\detF\) éq 5.1-2

\(J\sigma =\tau\) éq 5.1-3

\(\tau =sb{}^{e}\text{}\) éq 5.1-4

\({b}^{e}=\text{Id}-2e\) éq 5.1-5

• Equations d’état:

\(s=-\left[2\mu \tilde{e}+K\text{tr}e\text{Id}+\mathrm{3K}\alpha \Delta T\text{Id}\right]\) éq. 5.1-6

\(A=-R(p)\) éq 5.1-7

Par la suite, on exprime les lois d’écoulement et le critère de plasticité directement en fonction du tenseur des déformations élastiques \(e\) .

• Lois d’écoulement

\(\begin{array}{}{D}^{p}\simeq -\frac{1}{2}F{\dot{G}}^{p}{F}^{T}=-\frac{1}{2\Delta t}\left[\underset{{b}^{e}}{\underset{\underbrace{}}{{\text{FG}}^{p}{F}^{T}}}-\Delta F\underset{{b}^{e-}}{\underset{\underbrace{}}{{F}^{-}{G}^{\text{p-}}{F}^{-T}}}\Delta {F}^{T}\right]\\ \text{}=-\frac{1}{2\Delta t}\left[\text{Id}-2e-\Delta F\left\lbrace \text{Id}-2{e}^{-}\right\rbrace \Delta {F}^{T}\right]\\ \text{}=(e-\underset{{e}^{\mathrm{Tr}}}{\underset{\underbrace{}}{\frac{1}{2}\left[\text{Id}-\Delta F\left\lbrace \text{Id}-2{e}^{-}\right\rbrace \Delta {F}^{T}\right]}})/\Delta t=(e-{e}^{\text{Tr}})/\Delta t\end{array}\) éq 5.1-8

En prenant les parties traces et déviatoriques de l’équation [éq 4.1-5], on obtient:

\(\text{tr}e-\text{tr}{e}^{\mathrm{Tr}}=\Delta pDf\exp(\frac{-3K\alpha \Delta T}{{\sigma}_{1}})\exp(\frac{-K\text{tr}e}{{\sigma}_{1}})\) éq 5.1-9

\(\tilde{e}=\lbrace \begin{array}{}\tilde{{e}^{\text{Tr}}}-\frac{3}{2}\Delta p\frac{\tilde{e}}{{e}_{\text{eq}}}\text{}\text{si solution régulière}\\ 0\text{et}\mathrm{\Delta p}\ge \frac{2}{3}(\tilde{e}-{\tilde{e}}^{\text{Tr}}{)}_{\text{eq}}\text{}\text{si solution singulière}\end{array}\) éq 5.1-10

• Conditions de cohérence

\(\begin{array}{}F=\lbrace \begin{array}{cc}2\mu {e}_{\mathrm{eq}}+{\sigma}_{1}Df\exp(\frac{-3K\alpha \Delta T}{{\sigma}_{1}})\exp(\frac{-K\mathrm{tr}e}{{\sigma}_{1}})-R-{\sigma}_{y}& \text{si solution régulière}\\ {\sigma}_{1}Df\exp(\frac{-3K\alpha \Delta T}{{\sigma}_{1}})\exp(\frac{-K\mathrm{tr}e}{{\sigma}_{1}})-R-{\sigma}_{y}& \text{si solution singulière}\end{array}\\ \text{avec}F\le 0\Delta p\ge 0F\Delta p=0\end{array}\) éq 5.1-11

Résolution du système non linéaire

L’intégration de la loi de comportement se résume donc à résoudre le système [éq 5.1-9], [éq 5.1-10] et [éq 5.1-11]. Nous verrons que cette résolution se ramène à celle d’une seule équation scalaire, dont l’inconnue \(x\) est l’incrément de la trace des déformations élastiques:

\(x=\text{tr}e-\text{tr}{e}^{\mathrm{Tr}}\) éq 5.2-1

Grâce à ce choix, que la solution soit élastique ou plastique, atteinte en un point singulier ou non, l’équation [éq 5.1-9] portant sur la trace de l’incrément élastique est toujours valide et permet d’exprimer l’incrément de déformation plastique cumulée:

\(\begin{array}{}\text{tr}e-\text{tr}{e}^{\mathrm{Tr}}=\Delta p\underset{G}{\underset{\underbrace{}}{Df\exp(\frac{-K\text{tr}{e}^{\mathrm{Tr}}}{{\sigma}_{1}})\exp(\frac{-3K\alpha \Delta T}{{\sigma}_{1}})}}\exp(\frac{-K(\text{tr}e-\text{tr}{e}^{\mathrm{Tr}})}{{\sigma}_{1}})\\ \Rightarrow \Delta p(x)=\frac{1}{G}x\exp\frac{Kx}{{\sigma}_{1}}\end{array}\) éq 5.2-2

Cette équation nous montre qu’on peut chercher \(x\ge 0\) pour garantir une déformation plastique cumulée positive et que la solution élastique est obtenue pour \(x=0\) . On remarque également que l’incrément de déformation plastique cumulée est une fonction continue et strictement croissante de \(x\) . Moyennant ces remarques, si on note par S le terme [éq 5.2-3] dans le critère de plasticité, il s’agit alors, là aussi, d’une fonction continue et strictement croissante de \(x\) :

\(F=2\mu {e}_{\text{eq}}-S(x)\text{}\text{avec}\text{}S(x)=-{\sigma}_{1}\text{G exp}(-\frac{\text{Kx}}{{\sigma}_{1}})+R(p(x))+{\sigma}_{y}\) éq 5.2-3

A ce stade, la démarche de résolution se décompose en deux temps.

Examen des points singuliers

Un tel point singulier est caractérisé par [éq 5.1-10] (bas) et [éq 5.1-11] (bas), donc en particulier par \(S(x)=0\) . Du fait des propriétés de \(S\) , cette équation admet au plus une solution positive, disons \({x}^{S}\) qui existe si et seulement si \(S(0)\le 0\) . La connaissance de \({x}^{S}\) permet d’en déduire le tenseur de déformations élastiques \(e\) , la déformation plastique cumulée \(p\) ainsi que les forces thermodynamiques \(s\) et \(R\) .

Finalement, ce point singulier sera solution si l’inégalité dans [éq 5.1-11] (bas) est vérifiée, c’est-à-dire si:

\({\mathrm{\Delta p}}^{s}\ge \frac{2}{3}({\tilde{e}}^{s}-{\tilde{e}}^{\text{Tr}}{)}_{\text{eq}}\) éq 5.2.1-1

Solution régulière

L’équation d’écoulement [éq 5.1-10] (haut) qui détermine la partie déviatorique du tenseur de déformation élastique permet d’en déduire une équation scalaire fonction de l’incrément de déformation plastique cumulée:

\(\tilde{e}-{\tilde{e}}^{\text{Tr}}=-\frac{3}{2}\mathrm{\Delta p}\frac{\tilde{e}}{{e}_{\text{eq}}}\text{}\Rightarrow \text{}\lbrace \begin{array}{}{e}_{\text{eq}}={e}_{\text{eq}}^{\text{Tr}}-\frac{3}{2}\Delta p\\ \tilde{e}=\frac{{e}_{\text{eq}}}{{e}_{\text{eq}}^{\text{Tr}}}{\tilde{e}}^{\text{Tr}}\end{array}\) éq 5.2.2-1

On constate qu’en raison de la positivité de \({e}_{\text{eq}}\) , la valeur licite de \(\Delta p\) est bornée:

\(\Delta p\le \frac{2}{3}{e}_{\text{eq}}^{\text{Tr}}\) éq 5.2.2-2

La condition de cohérence détermine maintenant \(x\) :

\(F=2\mu {e}_{\text{eq}}^{\text{Tr}}-S(x)-3\mu \Delta p\le 0\) éq 5.2.2-3

Etant donné cette expression, la majoration de la valeur licite de \(\Delta p\) se réduit à la seule condition \(S(x)\ge 0\) ou, de manière équivalente, à \(x\ge {x}^{S}\) .

La solution élastique est obtenue pour \(x\) = 0. C’est la solution du problème si et seulement si:

\(F(0)=2\mu {e}_{\text{eq}}^{\text{Tr}}-S(0)<0\) éq 5.2.2-4

Dans le cas contraire, on doit alors résoudre:

\(F(x)=2\mu {e}_{\text{eq}}^{\text{Tr}}-S(x)-\frac{3\mu }{G}x\exp(\frac{\text{Kx}}{{\sigma}_{1}})=0\text{}\text{avec}\text{}\lbrace \begin{array}{}x>{x}^{s}\text{}\text{si}{x}^{s}\text{existe}\\ x>0\text{}\text{sinon}\end{array}\) éq 5.2.2-5

Cette fonction est continue et strictement décroissante et tend vers \(-\infty\) avec \(x\) . Elle admet donc au plus une solution. La démonstration de l’existence de cette solution est immédiate. En effet, il suffit de prouver que \(F\) est positive sur la borne inférieure de l’intervalle de recherche.

Lorsque \({x}^{S}\) n’existe pas, \(F(0)>0\) puisque la solution n’est pas élastique.

Lorsque \({x}^{S}\) existe, la fonction vaut:

\(F({x}^{s})=2\mu {e}_{\text{eq}}^{\text{Tr}}-3\mu \Delta {p}^{s}>0\iff \Delta {p}^{s}<\frac{2}{3}{e}_{\text{eq}}^{\text{Tr}}\) éq 5.2.2-6

Cette condition est vérifiée puisqu’on a rejeté la solution singulière.

Déroulement du calcul

La démarche pour résoudre l’ensemble des équations du modèle est la suivante:

1. On recherche si la solution est élastique

• calcul de \(F(0)\)

• si \(F(0)<0\) , la solution du problème est la solution élastique \({x}^{\text{Sol}}=0\)

• sinon on passe en 2)

  • • • • • • 2. Si \(S(0)>0\) , la solution est plastique et régulière

• on passe en 4)

3. Si \(S(0)<0\) , on cherche si la solution est singulière

• on résout \(S({x}^{s})=0\)

• si \({x}^{s}\) vérifie l’inégalité \({\mathrm{\Delta p}}^{s}\ge \frac{2}{3}({\tilde{e}}^{s}-{\tilde{e}}^{\text{Tr}}{)}_{\text{eq}}\) , alors la solution est singulière \({x}^{\text{Sol}}={x}^{s}\)

• sinon, \({x}^{s}\) est une borne inférieure pour résoudre \(F(x)=0\) , on passe en 4)

4. La solution est plastique et régulière

• on résout \(F(x)=0\)

Résolution

Pour résoudre les deux équations \(S(x)=0\) et \(F(x)=0\) , on emploie une méthode de Newton avec bornes contrôlées couplée à de la dichothomie lorsque Newton donne une solution en dehors de l’intervalle des deux bornes. On présente maintenant la détermination des bornes pour chacun des cas précédents (points 2) 3) et 4) du paragraphe précédent).

Bornes supérieure et inférieure dans le cas où la fonction S est strictement positive à l’origine

On résout:

\(\lbrace \begin{array}{}F(x)=0\\ F(0)>0\end{array}\iff \lbrace \begin{array}{}\underset{{f}_{1}}{\underset{\underbrace{}}{2\mu {e}_{\text{eq}}^{\text{Tr}}-S(x)}}=\underset{3\mu \Delta p}{\underset{\underbrace{}}{\frac{3\mu }{G}x\exp(\frac{\text{Kx}}{{\sigma}_{1}})}}\\ {f}_{1}(0)>0\end{array}\) éq 5.4.1-1

où la fonction \(\mathrm{\Delta p}(x)\) est continue, strictement croissante et nulle à l’origine et la fonction \({f}_{1}(x)\) est continue, strictement décroissante et strictement positive à l’origine (voir [Figure 5.4.1-a]).

On pose:

\({f}_{1}=\underset{{f}_{2}}{\underset{\underbrace{}}{2\mu {e}_{\text{eq}}^{\text{Tr}}-R(x)-{\sigma}_{y}}}+{\sigma}_{1}\text{Gexp}(-\frac{\text{Kx}}{{\sigma}_{1}})\text{}\text{alors}\text{}{f}_{2}(x)<{f}_{1}(x)\text{}\forall x\ge 0\) éq 5.4.1-2

où la fonction \({f}_{2}(x)\) est continue, strictement décroissante. Dans ce cas, la résolution des équations:

\({f}_{2}(\Delta {p}^{\inf})=3\mu \Delta {p}^{\inf}\) et \({x}^{\inf}\exp(\frac{{\text{Kx}}^{\inf}}{{\sigma}_{1}})=G\Delta {p}^{\inf}\) éq 5.4.1-3

pour en déduire successivement \(\mathrm{\Delta p}\) puis \(x\) donne une borne inférieure \({x}^{\text{Inf}}\) qui correspond à la solution du modèle à écrouissage isotrope et critère de Von Mises. Si \({f}_{2}(0)<0\) , la borne inférieure est prise égale à zéro: \({x}^{\inf}=0\) .

La borne supérieure \({x}^{\text{Sup}}\) est telle que:

\({x}^{\text{Sup}}\exp(\frac{{\text{Kx}}^{\text{Sup}}}{{\sigma}_{1}})=\frac{G}{3\mu }{f}_{1}({x}^{\text{Inf}})\) éq 5.4.1-4

L’équation du type \(x\exp(\frac{\text{Kx}}{{\sigma}_{1}})=\text{constante}\) est résolue par une méthode de Newton.

Figure 5.4.1-a: représentation graphique des bornes supérieure et inférieure

Bornes supérieure et inférieure dans le cas où la fonction S est négative ou nulle à l’origine

Le système à résoudre est le suivant:

\(\lbrace \begin{array}{}S(x)=0\\ S(0)<0\end{array}\text{}\iff \text{}\lbrace \begin{array}{}R({p}^{-}+\frac{x}{G}\exp(\frac{\text{Kx}}{{\sigma}_{1}}))+{\sigma}_{y}={\sigma}_{1}\text{Gexp}(-\frac{\text{Kx}}{{\sigma}_{1}})\\ R({p}^{-})+{\sigma}_{y}<{\sigma}_{1}G\end{array}\) éq 5.4.2-1

La partie de gauche est une fonction continue, strictement croissante de \(x\) et strictement positive à l’origine, la partie de droite est une fonction continue, strictement décroissante de \(x\) et strictement positive à l’origine. Utilisant les propriétés de ces deux fonctions, une représentation graphique (cf.[Figure 5.4.2-a]) de ces fonctions montre que la borne supérieure \({x}^{\text{Sup}}\) est telle que:

\({\sigma}_{1}\text{G exp}(-\frac{{\text{Kx}}^{\text{Sup}}}{{\sigma}_{1}})=R({p}^{-})+{\sigma}_{y}\text{}\iff \text{}{x}^{\text{Sup}}=\frac{{\sigma}_{1}}{K}\log(\frac{{\sigma}_{1}G}{R({p}^{-})+{\sigma}_{y}})\) éq 5.4.2-2

La borne inférieure \({x}^{\text{Inf}}\) est telle que:

\(\begin{array}{}{\sigma}_{1}\text{G exp}(-\frac{{\text{Kx}}^{\text{Inf}}}{{\sigma}_{1}})=R({p}^{-}+\frac{{x}^{\text{Sup}}}{G}\exp(\frac{{\text{Kx}}^{\text{Sup}}}{{\sigma}_{1}}))+{\sigma}_{y}\\ \iff {x}^{\text{Inf}}={\langle \frac{{\sigma}_{1}}{K}\log(\frac{{\sigma}_{1}G}{R({p}^{-}+\frac{{x}^{\text{Sup}}}{G}\exp(\frac{{\text{Kx}}^{\text{Sup}}}{{\sigma}_{1}}))+{\sigma}_{y}})\rangle }^{+}\end{array}\) éq 5.4.2-3

Figure 5.4.2-a: représentation graphique des bornes supérieure et inférieure

Bornes supérieure et inférieure dans le cas où la fonction :math:`S`est strictement négative à l’origine et :math:`{x}^{s}`non solution

On résout le système suivant:

\(\lbrace \begin{array}{}F(x)=0\\ S(0)<0\\ S({x}^{s})=0\end{array}\text{}\iff \text{}\lbrace \begin{array}{}\underset{{f}_{1}}{\underset{\underbrace{}}{2\mu {e}_{\text{eq}}^{\text{Tr}}-S(x)}}=\underset{3\mu \Delta p}{\underset{\underbrace{}}{\frac{3\mu }{G}x\exp(\frac{\text{Kx}}{{\sigma}_{1}})}}\\ {f}_{1}(0)>0\\ 2\mu {e}_{\text{eq}}^{\text{Tr}}=\frac{3\mu }{G}{x}^{s}\exp(\frac{{\text{Kx}}^{s}}{{\sigma}_{1}})\end{array}\) éq 5.4.3-1

La solution \({x}^{\text{Sol}}\) est telle que \(S({x}^{\text{Sol}})>0\) .

Pour la borne inférieure, on prend \({x}^{\text{Inf}}={x}^{s}\) . Etant données les propriétés des fonctions \({f}_{1}\) (strictement décroissante) et \(3\mu \Delta p(x)\) (strictement croissante), la borne supérieure \({x}^{\text{Sup}}\) est telle que (cf. [Figure 5.4.3-a]):

\({x}^{\text{Sup}}\exp(\frac{{\text{Kx}}^{\text{Sup}}}{{\sigma}_{1}})=\frac{\mathrm{2G}}{3}{e}_{\text{eq}}^{\text{Tr}}\) éq 5.4.3-2

Cette équation est résolue par une méthode de Newton.

Figure 5.4.3-a: représentation graphique des bornes supérieure et inférieure

Correction de volume a posteriori

Dans le cas du modèle de Rousselier, le changement de volume joue un rôle crucial, si bien que l’erreur commise par l’approximation \(s=\tau\) dans les équations d’écoulement peut conduire à une évolution peut précise de la porosité, cf. [Bib2]. En suivant la proposition de cette référence, on va corriger a posteriori (i.e. après avoir calculé toutes les quantités) la trace de la déformation élastique.

En effet, en l’absence d’approximation, la partie hydrostatique de l’équation d’écoulement conduit à:

\(\frac{d}{\mathrm{dt}}(\ln{J}^{p})=\dot{p}Df\exp(\frac{\text{tr}\tau }{3{\sigma}_{1}})\) (éq. 5.5-1)

Après discrétisation en temps, on obtient alors une proposition corrigée pour les changements de volume plastique et élastique:

\({J}_{\mathrm{corr}}^{p}={J}^{p-}\exp(\text{tr}e-\text{tr}{e}^{\mathrm{Tr}});{J}_{\mathrm{corr}}^{e}=\frac{J}{{J}_{\mathrm{corr}}^{p}}\) (éq. 5.5-2)

On va alors chercher une nouvelle trace de déformation élastique telle que le changement de volume élastique corresponde à la valeur corrigée ci-dessus:

\({e}_{\mathrm{corr}}=\tilde{e}+t\mathrm{Id}\text{où}t\text{tel que :}\det(\mathrm{Id}-2{e}_{\mathrm{corr}})={{J}_{\mathrm{corr}}^{e}}^{2}\) (éq. 5.5-3)

Cela conduit à un polynôme de degré 3 en \(t\) , dont on choisira la solution la plus proche de \(e\) .

Expression de la matrice tangente du comportement

On donne ici l’expression de la matrice tangente (option FULL_MECA au cours des itérations de Newton, option RIGI_MECA_TANG pour la première itération).

Pour l’option FULL_MECA, celle-ci est obtenue en linéarisant le système d’équations qui régit la loi de comportement. Nous donnons ci-après les grandes lignes de cette linéarisation.

Pour l’option RIGI_MECA_TANG, il s’agit des mêmes expressions que celles données pour FULL_MECA mais avec \(\Delta p=0\) . En particulier, on aura \(\Delta F=\text{Id}\) .

La loi de comportement peut se mettre sous la forme générale suivante:

\(\sigma =\sigma (\tau ,\Delta F)\) éq 5.6-1

\(\tau =\tau (e)\) éq 5.6-2

\(e=e({e}^{\mathrm{Tr}},f)\) éq 5.6-3

\({e}^{\text{Tr}}={e}^{\text{Tr}}(\Delta F)\) éq 5.6-4

La linéarisation de ce système donne:

\(\delta \sigma =(\frac{\partial \sigma }{\partial \tau }:\frac{\partial \tau }{\partial e}:(\frac{\partial e}{\partial {e}^{\mathrm{Tr}}}:\frac{\partial {e}^{\mathrm{Tr}}}{\partial \Delta F}+\frac{\partial e}{\partial f}\frac{\partial f}{\partial J}\otimes \frac{\partial J}{\partial \Delta F})+\frac{\partial \sigma }{\partial \Delta F}):\delta \Delta F=H:\delta \Delta F\) éq 5.6-5

où \(H\) est la matrice tangente. Par la suite, on détermine séparément les cinq termes de l’équation précédente.

Dans la linéarisation du système, on utilisera souvent le tenseur \(C\) défini ci-dessous et les deux équations suivantes:

\({\mathrm{\delta a}}_{ij}=\frac{1}{2}({\delta}_{\text{ik}}{\delta}_{\text{jl}}+{\delta}_{\text{jk}}{\delta}_{\text{il}}){\mathrm{\delta a}}_{kl}\) éq 5.6-6

\(\delta {a}_{\text{pp}}=\delta {}_{kl}\delta {a}_{kl}\) éq 5.5-7

\({C}_{ijkl}=\frac{1}{2}({\delta}_{\text{ik}}{\delta}_{\text{jl}}+{\delta}_{\text{jk}}{\delta}_{\text{il}})\) éq 5.6-8

• Calcul de \(\frac{\partial \sigma }{\partial \tau }\) et de \(\frac{\partial \sigma }{\partial \Delta F}\)

La linéarisation de la relation qui lie la contrainte de Cauchy \(\sigma\) et la contrainte de Kirchhoff \(\tau\) donne:

\(J\sigma =\tau \text{}\iff \text{}\delta \sigma =\frac{1}{J}\delta \tau -(\frac{\sigma}{J}\otimes \frac{\partial J}{\partial \Delta F}):\delta \Delta F\) éq 5.6-9

En utilisant la relation [éq 5.6-6], on obtient pour \(\frac{\partial \sigma }{\partial \tau }\) :

\(\frac{\partial \sigma }{\partial \tau }=C\) éq 5.6-10

et pour \(\frac{\partial \sigma }{\partial \Delta F}\) :

\(\frac{\partial \sigma }{\partial \Delta F}=-\frac{\sigma}{J}\otimes \frac{\partial J}{\partial \Delta F}\) éq 5.6-11

avec

\(\begin{array}{}\frac{\partial J}{\partial {F}_{11}}={\mathrm{\Delta F}}_{22}{\mathrm{\Delta F}}_{33}-{\mathrm{\Delta F}}_{23}{\mathrm{\Delta F}}_{32}\\ \frac{\partial J}{\partial {F}_{22}}={\mathrm{\Delta F}}_{11}{\mathrm{\Delta F}}_{33}-{\mathrm{\Delta F}}_{13}{\mathrm{\Delta F}}_{31}\\ \frac{\partial J}{\partial {F}_{33}}={\mathrm{\Delta F}}_{11}{\mathrm{\Delta F}}_{22}-{\mathrm{\Delta F}}_{12}{\mathrm{\Delta F}}_{21}\\ \frac{\partial J}{\partial {F}_{12}}={\mathrm{\Delta F}}_{31}{\mathrm{\Delta F}}_{23}-{\mathrm{\Delta F}}_{33}{\mathrm{\Delta F}}_{21}\text{}\frac{\partial J}{\partial {F}_{21}}={\mathrm{\Delta F}}_{13}{\mathrm{\Delta F}}_{32}-{\mathrm{\Delta F}}_{33}{\mathrm{\Delta F}}_{12}\\ \frac{\partial J}{\partial {F}_{13}}={\mathrm{\Delta F}}_{21}{\mathrm{\Delta F}}_{32}-{\mathrm{\Delta F}}_{22}{\mathrm{\Delta F}}_{31}\text{}\frac{\partial J}{\partial {F}_{31}}={\mathrm{\Delta F}}_{12}{\mathrm{\Delta F}}_{23}-{\mathrm{\Delta F}}_{22}{\mathrm{\Delta F}}_{13}\\ \frac{\partial J}{\partial {F}_{23}}={\mathrm{\Delta F}}_{31}{\mathrm{\Delta F}}_{12}-{\mathrm{\Delta F}}_{11}{\mathrm{\Delta F}}_{32}\text{}\frac{\partial J}{\partial {F}_{32}}={\mathrm{\Delta F}}_{13}{\mathrm{\Delta F}}_{21}-{\mathrm{\Delta F}}_{11}{\mathrm{\Delta F}}_{23}\end{array}\) éq 5.6-12

• Calcul de \(\frac{\partial \tau }{\partial e}\)

La relation qui lie la contrainte de Kirchhoff \(\tau\) et le tenseur de déformation élastique \(e\) est donnée par:

\(\begin{array}{}\tau =s{b}^{e}=-2\mu e-\lambda \text{Tr}e\text{Id}+4\mu ee+2\lambda (\text{tr}e)e\\ -\mathrm{3K}\alpha \Delta T\text{Id}+\mathrm{6K}\alpha \Delta Te\end{array}\) éq 5.6-13

On obtient après linéarisation:

\(\delta \tau =2(\lambda \text{tr}e-\mu +\mathrm{3K}\alpha \Delta T)\delta e+\lambda (2e-\text{Id})\text{Tr}\delta e+4\mu (\delta ee+e\delta e)\) éq 5.6-14

d’où

\(\frac{\partial {\tau}_{ij}}{\partial {e}_{kl}}=2(\lambda \text{tr}e-\mu +\mathrm{3K}\alpha \Delta T){C}_{ijkl}+\lambda (2{e}_{ij}-{\delta}_{ij}){\delta}_{kl}+2\mu ({\delta}_{\text{ik}}{e}_{\text{lj}}+{\delta}_{\text{il}}{e}_{\text{kj}}+{e}_{\text{il}}{\delta}_{\text{kj}}+{e}_{\text{ik}}{\delta}_{\text{jl}})\) éq 5.6-15

• Calcul de \(\frac{\partial {e}^{\text{Tr}}}{\partial \Delta F}\)

La relation entre le tenseur de déformation élastique \({e}^{\text{Tr}}\) et l’incrément du gradient de la transformation \(\Delta F\) s’écrit:

\({e}^{\text{Tr}}=\frac{1}{2}(\text{Id}-\Delta F{b}^{e-}\Delta {F}^{T})\) éq 5.6-16

Sa linéarisation donne:

\(\frac{\partial {e}_{ij}^{\text{Tr}}}{\partial \Delta {F}_{kl}}=-\frac{1}{2}({\delta}_{\text{ik}}\Delta {F}_{\text{jp}}{b}_{\text{pl}}^{\text{e-}}+\Delta {F}_{\text{ip}}{b}_{\text{pl}}^{\text{e-}}{\delta}_{\text{jk}})\) éq 5.6-17

• Calcul de \(\frac{\partial e}{\partial {e}^{\text{Tr}}}\)

Cas élastique

Dans le cas élastique, le calcul de \(\frac{\partial e}{\partial {e}^{\text{Tr}}}\) est immédiat puisque \(\delta e=\delta {e}^{\text{Tr}}\) d’où

\(\frac{\partial e}{\partial {e}^{\text{Tr}}}=C\) éq 5.6-18

Cas plastique – Solution régulière

Pour déterminer \(\frac{\partial e}{\partial {e}^{\text{Tr}}}\) , on opère en deux étapes. Par la loi d’écoulement discrétisée, on calcule en premier \(\delta e\) en fonction de \(\delta {e}^{\text{Tr}}\) et \(\delta \Delta p\) . Puis la condition de cohérence permet d’en déduire \(\delta \Delta p\) en fonction de \(\delta {e}^{\text{Tr}}\) . Ces deux étapes sont détaillées par la suite.

La partie déviatorique de la loi d’écoulement discrétisée s’écrit:

\(\tilde{e}-{\tilde{e}}^{\text{Tr}}=-\frac{3}{2}\Delta p\frac{\tilde{e}}{{e}_{\mathrm{eq}}}\) éq 5.6-19

On obtient après linéarisation:

\(\underset{1/\alpha }{\underset{\underbrace{}}{(1+\frac{3}{2}\frac{\Delta p}{{e}_{\text{eq}}})}}\delta \tilde{e}=\delta {\tilde{e}}^{\text{Tr}}-\frac{3}{2}\frac{\tilde{e}}{{e}_{\text{eq}}}\delta \Delta p+\frac{9}{4}\Delta p\frac{\tilde{e}}{{e}_{\text{eq}}^{3}}(\tilde{e}:\delta \tilde{e})\) éq 5.6-20

Pour déterminer \(\tilde{e}:\delta \tilde{e}\) , on contracte l’équation [éq 5.6-20] avec \(\tilde{e}\) ce qui donne:

\(\tilde{e}:\delta \tilde{e}=\tilde{e}:\delta {\tilde{e}}^{\text{Tr}}-{e}_{\text{eq}}\delta \Delta p\) éq 5.6-21

d’où

\(\delta \tilde{e}=\underset{{A}_{1}}{\underset{\underbrace{}}{\alpha \left[\frac{9\Delta p}{4{e}_{\text{eq}}^{3}}\tilde{e}\otimes \tilde{e}+C\right]}}:\delta {\tilde{e}}^{\text{Tr}}\underset{{A}_{2}}{\underset{\underbrace{}}{-\frac{3}{2}\frac{\tilde{e}}{{e}_{\text{eq}}}}}\mathrm{\delta \Delta p}\) éq 5.6-22

Pour la partie trace de la loi d’écoulement discrétisée, on a:

\(\text{Tr}e-\text{Tr}{e}^{\text{Tr}}=Df\Delta p\exp(\frac{3K\alpha \Delta T}{{\sigma}_{1}})\exp(-\frac{K}{{\sigma}_{1}}\text{Tr}e)\) éq 5.6-23

ce qui donne, en posant \({k}_{1}=1+\frac{DfK\Delta p}{{\sigma}_{1}}\exp(\frac{3K\alpha \Delta T}{{\sigma}_{1}})\exp(-\frac{K}{{\sigma}_{1}}\text{Tr}e)\) :

\(\text{tr}\delta e=\begin{array}{}\underset{{\alpha}_{1}}{\underset{\underbrace{}}{\frac{1}{{k}_{1}}}}\text{tr}\delta {e}^{\mathrm{Tr}}\\ +\underset{{\alpha}_{2}}{\underset{\underbrace{}}{\frac{Df\exp(\frac{-K}{{\sigma}_{1}}\text{tr}e)\exp(\frac{3K\alpha \Delta T}{{\sigma}_{1}})}{{k}_{1}}}}\delta \Delta p\\ +\underset{{\beta}_{1}}{\underset{\underbrace{}}{\frac{\Delta p\exp(\frac{-K}{{\sigma}_{1}}\text{tr}e)\exp(\frac{3K\alpha \Delta T}{{\sigma}_{1}})}{{k}_{1}}}}D\delta f\end{array}\) éq 5.6-24

Dans le cas plastique, la condition de cohérence vaut:

\(2\mu {e}_{\text{eq}}+Df{\sigma}_{1}\exp(-\frac{3K\alpha \Delta T}{{\sigma}_{1}})\exp(-\frac{K}{{\sigma}_{1}}\text{Tr}e)-R-{\sigma}_{y}=0\) éq 5.6-25

d’où

\(\begin{array}{}\frac{3\mu }{{e}_{\mathrm{eq}}}(\tilde{e}:\delta \tilde{e})+D{\sigma}_{1}\exp(\frac{-3K\alpha \Delta T}{{\sigma}_{1}})\exp(\frac{-K\text{tr}e}{{\sigma}_{1}})\delta f\\ -DfK\exp(\frac{-3K\alpha \Delta T}{{\sigma}_{1}})\exp(\frac{-K\text{tr}e}{{\sigma}_{1}})\text{tr}\delta e-h\delta \Delta p=0\end{array}\) éq 5.6-26

En injectant la relation [éq 5.6-21] dans l’équation ci-dessus, on obtient alors, en posant

\({k}_{2}=3\mu +h+{\alpha}_{2}DfK\exp(\frac{-3K\alpha \Delta T}{{\sigma}_{1}})\exp(-\frac{K}{{\sigma}_{1}}\text{Tr}e)\) :

\(\delta \Delta p=\begin{array}{}\frac{3\mu }{{e}_{\mathrm{eq}}}\underset{{\alpha}_{3}}{\underset{\underbrace{}}{\frac{1}{{k}_{2}}}}\tilde{e}:\delta {\tilde{e}}^{\mathrm{Tr}}-\underset{{\alpha}_{4}}{\underset{\underbrace{}}{\frac{{\alpha}_{1}DfK\exp(\frac{-3K\alpha \Delta T}{{\sigma}_{1}})\exp(-\frac{K}{{\sigma}_{1}}\text{Tr}e)}{{k}_{2}}}}\text{tr}\delta {e}^{\mathrm{Tr}}\\ \underset{{\beta}_{2}}{\underset{\underbrace{}}{+\frac{\exp(\frac{-3K\alpha \Delta T}{{\sigma}_{1}})\exp(\frac{-K}{{\sigma}_{1}}\text{Tr}e)({\sigma}_{1}-{\beta}_{1}DfK)}{{k}_{2}}}}D\delta f\end{array}\) éq 5.6-27

En remplaçant \(\delta \Delta p\) par sa valeur obtenue ci-dessus dans les équations [éq 5.6–22] et [éq 5.6-24], on obtient:

\(\begin{array}{}\delta e=\underset{\mathrm{ddvetr}}{\underset{\underbrace{}}{\left[{A}_{1}+({A}_{2}+\frac{1}{3}{\alpha}_{2}\mathrm{Id})\otimes (\frac{3\mu {\alpha}_{3}}{{e}_{\mathrm{eq}}}\tilde{e})\right]}}:\delta \tilde{{e}^{\mathrm{Tr}}}+\underset{\mathrm{dtretr}}{\underset{\underbrace{}}{\left[\frac{1}{3}{\alpha}_{1}\mathrm{Id}+{\alpha}_{4}({A}_{2}+\frac{1}{3}{\alpha}_{2}\mathrm{Id})\right]}}\text{tr}\delta {e}^{\mathrm{Tr}}\\ +\underset{\mathrm{dedf}}{\underset{\underbrace{}}{\left[\frac{1}{3}{\beta}_{1}\mathrm{Id}+{\beta}_{2}({A}_{2}+\frac{1}{3}{\alpha}_{2}\mathrm{Id})\right]D}}\delta f\end{array}\) éq 5.6-28

d’où

\(\frac{\partial e}{\partial {e}^{\text{Tr}}}=\text{ddvetr}+(\text{dtretr}-\frac{1}{3}\text{ddvetr}:\text{Id})\otimes \text{Id}\) éq 5.6-29

Cas plastique – Solution singulière

La démarche est identique à celle utilisée précédemment.

On obtient pour la loi d’écoulement discrétisée:

\(\tilde{e}=0\text{}\iff \text{}\delta \tilde{e}=0\) éq 5.6-30

pour la partie déviatorique et pour la partie trace, la relation est identique à celle trouvée pour la solution régulière.

\(\text{tr}\delta e={\alpha}_{1}\text{tr}\delta {e}^{\mathrm{Tr}}+{\alpha}_{2}\delta \Delta p+{\beta}_{1}D\delta f\) éq 5.6-31

où \({\alpha}_{1}\) , \({\alpha}_{2}\) et \({\beta}_{1}\) ont les mêmes définitions qu’au paragraphe précédent.

La condition de cohérence permet alors de trouver \(\delta \Delta p\) en fonction de \(\partial {e}^{\text{Tr}}\) .

\(Df{\sigma}_{1}\exp(\frac{\mathrm{3K}\alpha \Delta T}{{\sigma}_{1}})\exp(-\frac{K}{{\sigma}_{1}}\text{Tr}e)-R-{\sigma}_{y}=0\) éq 5.6-32

soit après linéarisation:

\(\delta \Delta p={\alpha}_{4}\text{tr}\delta {e}^{\mathrm{Tr}}+{\beta}_{2}D\delta f\) éq 5.6-33

soit finalement:

\(\delta e=\underset{\mathrm{dtretr}}{\underset{\underbrace{}}{\frac{1}{3}\left[{\alpha}_{1}+{\alpha}_{4}{\alpha}_{2}\right]\mathrm{Id}\text{tr}\delta {e}^{\mathrm{Tr}}}}+\underset{\mathrm{dedf}}{\underset{\underbrace{}}{\frac{1}{3}\left[{\beta}_{1}+{\beta}_{2}{\alpha}_{2}\right]\mathrm{Id}D\delta f}}\) éq 5.6-34

d’où

\(\frac{\partial e}{\partial {e}^{\text{Tr}}}=\text{dtretr}\otimes \text{Id}\) éq 5.6-35

• Calcul de \(\frac{\partial f}{\partial J}\)

Compte tenu de la relation 4.3-1, la dérivée s’écrit simplement:

\(\lbrace \begin{array}{cc}\frac{\partial f}{\partial J}=D\frac{1-{f}_{0}}{{J}^{2}}& \text{si}f>{f}_{0}\\ \frac{\partial f}{\partial J}=0& \text{si}f={f}_{0}\end{array}\) éq 5.6-36

Remarque:

La matrice tangente n’est pas altérée par la correction de volume car celle-ci est effectuée a posteriori, c’est-à-dire après le calcul des contraintes.

Bibliographie

1. ROUSSELIER G. : “Finite deformation constitutive relations including ductile fracture damage, In three dimensional constitutive relations and ductile fracture”, Ed. Nemat-Nasser, North Holland publishing company, pp. 331-355.

2. LORENTZ E., BESSON J., CANO V., “Numerical simulation of ductile fracture: an efficient and robust implementation of the Rousselier constitutive law”, Comp. Meth. Appl. Mech. Engrg. 197, pp. 1965-1982, 2008.

3. SIMO J.C., MIEHE C., « Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation », Comp. Meth. Appl. Mech. Eng., 98, pp. 41-104, North Holland, 1992.SIDOROFF F., « Cours sur les grandes déformations », Rapport Greco n51, 1982.

4. LORENTZ E., CANO V.: “A minimisation principle for finite strain plasticity: incremental objectivity and immediate implementation”, Comm. Num. Meth. Eng. 18, pp. 851-859, 2002.

5. MIALON P.: “Eléments d’analyse et de résolution numérique des relations de l’élasto‑plasticité“, EDF, bulletin de la DER, série C mathématiques informatique, 3, pp. 57-89, 1986

Description des versions du document

Version Aster	Auteur(s) ou contributeur(s), organisme	Description des modifications
6.0	V.Cano, E.Lorentz EDF/R&D/AMA	Texte initial
9.4	E.Lorentz EDF/R&D/SINETICS	Modification de l”expression de la porosité, fonction de du changement de volume J (ce qui rend le modèle entièrement implicite) et introduction d’une correction sur le changement de volume élastique a posteriori.
12.1	J.M.Proix EDF/R&D/AMA	Fiche 19650: modification de l’ordre des variables internes