r5.03.02 Intégration des relations de comportement élasto-plastique de Von Mises

Résumé :

Ce document décrit les quantités calculées par l’opérateur STAT_NON_LINE nécessaires à la mise en œuvre de l’algorithme non linéaire quasi statique décrit en [R5.03.01] dans le cas des comportements élasto-plastiques. Ces quantités sont calculées par les mêmes sous-programmes dans l’opérateur DYNA_NON_LINE dans le cas d’une sollicitation dynamique [R5.05.05].

Cette description est présentée suivant les différents mots clés qui permettent à l’utilisateur de choisir la relation de comportement souhaitée. Les relations de comportement traitées ici sont:

• le comportement de Von Mises à écrouissage isotrope (linéaire ou non linéaire)

• le comportement de Von Mises à écrouissage cinématique linéaire (modèle de Prager)

La méthode d’intégration utilisée se base sur une formulation implicite directe. A partir de l’état initial, ou à partir de l’instant de calcul précédent, on calcule le champ de contraintes résultant d’un incrément de déformation. On calcule également l’opérateur tangent.

Notations générales et hypothèses sur les déformations

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.

Les incréments sont désignés par \(\Delta\) . On a ainsi:

\(\mathrm{Q}=\mathrm{Q}(t+\Delta t)=\mathrm{Q}(t)+\Delta \mathrm{Q}={\mathrm{Q}}^{-}+\Delta \mathrm{Q}\) .

Pour le calcul des dérivées , on notera: \(\dot{Q}\) dérivée de \(Q\) par rapport au temps

\(\sigma\)	tenseur des contraintes.
\(\text{~}\)	opérateur déviatoire : \({\tilde{\sigma}}_{ij}=\sigma -\frac{1}{3}{\sigma}_{kk}{\delta}_{ij}\) .
\({()}_{\mathrm{eq}}\)	valeur équivalente de Von Mises: \({\sigma}_{\mathit{eq}}=\sqrt{\frac{3}{2}{\stackrel{~}{\sigma}}_{ij}{\stackrel{~}{\sigma}}_{ij}}\)
\(\Delta \varepsilon\)	incrément de déformation.
\(A\)	tenseur d’élasticité.
\(\lambda ,\mu ,E,\nu ,K\)	coefficients de l’élasticité isotrope, respectivement: coefficients de Lamé, module d’Young, coefficient de Poisson et module de compressibilité.
\(3K=3\lambda +2\mu\)	module de compressibilité
\(\alpha\)	coefficient de dilatation thermique moyen.
\(t\)	temps.
\(T\)	température.
\({()}_{+}\)	partie positive.

Pour calculer les opérateurs tangents, on adoptera la convention d’écriture des tenseurs symétriques d’ordre deux sous forme de vecteurs à six composantes. Ainsi, pour un tenseur \(a\) :

\(\vec{a}={}^{t}\text{}\left[\begin{array}{ccc}{a}_{xx}& {a}_{yy}& {a}_{zz}\end{array}\begin{array}{ccc}\sqrt{2}{a}_{xy}& \sqrt{2}{a}_{xz}& \sqrt{2}{a}_{yz}\end{array}\right]\)

On introduit le vecteur hydrostatique \(\vec{1}\) et la matrice de projection déviatorique \(\mathrm{P}\) :

\(\vec{1}={}^{t}\text{}\left[\begin{array}{ccc}1& 1& 1\end{array}\begin{array}{ccc}0& 0& 0\end{array}\right]\)

\(\mathrm{P}=\text{Id}-\frac{1}{3}\vec{1}\otimes \vec{1}\)

Partition des déformations (petites déformations)

On écrit pour tout instant:

\(\epsilon (t)={\epsilon}^{e}(t)+{\epsilon}^{\mathit{th}}(t)+{\epsilon}^{p}(t)\)

avec la déformation élastique

\({\epsilon}^{e}(t)={\mathrm{A}}^{-1}(T(t))\sigma (t)\)

avec la déformation thermique

\({\epsilon}^{\mathit{th}}(t)=\alpha (T(t))\left(T(t)-{T}_{\mathit{ref}}\right)\mathrm{Id}\)

ou de façon plus générale:

\(\begin{array}{cccc}& {\epsilon}^{\text{th}}\left(T\right)& =& \alpha \left(T\right)\left(T-{T}_{\text{def}}\right)-\alpha \left({T}_{\text{ref}}\right)\left({T}_{\text{ref}}-{T}_{\text{def}}\right)\\ & & =& \widehat{\alpha}\left(T\right)\left(T-{T}_{\text{ref}}\right)\\ \text{et}& {\epsilon}^{\text{th}}\left({T}_{\text{ref}}\right)& =& 0\end{array}\)

La matrice d’élasticité \(A\) dépend de l’instant \(t\) par l’intermédiaire de la température. Le coefficient de dilatation thermique \(\alpha (T(t))\) est un coefficient de dilatation moyen qui peut dépendre de la température \(T\) . La température \({T}_{\mathrm{ref}}\) est la température de référence, c’est-à-dire celle pour laquelle la dilatation thermique est supposée nulle si le coefficient de dilatation moyen n’est pas connu par rapport à \({T}_{\mathrm{ref}}\) , on peut utiliser une température de définition du coefficient de dilatation moyen \({T}_{\text{def}}\) (définie par le mot-clé TEMP_DEF_ALPHA de DEFI_MATERIAU) différente de la température de référence [R4.08.01].

Ce qui conduit à : \(\dot{\epsilon}(t)=\stackrel{.}{\stackrel{⏞}{{\mathrm{A}}^{-1}(T(t))\sigma (t)}}+{\dot{\epsilon}}^{\mathit{th}}(t)+{\dot{\epsilon}}^{p}(t)\)

Ce choix est fait par souci de cohérence avec l’élasticité : il faut pouvoir trouver la même solution en élasticité (opérateur MECA_STATIQUE) et en élasto-plasticité (opérateur STAT_NON_LINE) lorsque les caractéristiques du matériau restent élastiques. Ce choix mène à la discrétisation:

\(\Delta \varepsilon =\Delta {\varepsilon}^{p}+\Delta ({A}^{-1}\sigma )+\Delta {\varepsilon}^{\mathrm{th}}\)

avec:

\(\Delta ({A}^{-1}\sigma )={A}^{-1}({t}^{-}+\Delta t)({\sigma}^{-}+\Delta \sigma )-{A}^{-1}({t}^{-}){\sigma}^{-}\)

et

\(\Delta {\varepsilon}^{\mathrm{th}}=(\alpha ({t}^{-}+\Delta t)(T-{T}_{\mathrm{ref}})-\alpha ({t}^{-})({T}^{-}-{T}_{\mathrm{ref}}))\mathrm{Id}\)

Réactualisation

Dans STAT_NON_LINE, sous le mot clé facteur COMPORTEMENT, plusieurs modes de calcul des déformations sont possibles :

• “PETIT”

• ‘SIMO_MIEHE’ [R5.03.21] (qui effectue le calcul en grandes déformations pour un écrouissage isotrope)

• ‘GDEF_LOG’ [R5.03.24] qui effectue le calcul en grandes déformations, en utilisant une formulation de type logarithmique

• ‘GROT_GDEP’ [R5.03.22] (qui effectue le calcul en grands déplacements et grandes rotations, mais en petites déformations)

• “PETIT_REAC” (qui est un substitut au calcul en grandes déformations, valable pour de petits incréments de charge, et pour de petites rotations [bib2]).

Cette dernière possibilité consiste à réactualiser la géométrie avant de calculer \(\Delta \epsilon\) :

On écrit \(x={x}_{0}+{u}_{i-1}+\Delta {u}_{i}^{n}\) , le calcul des gradients de \(\Delta {u}_{i}^{n}\) est donc fait avec la géométrie

au lieu de la géométrie initiale \({x}_{0}\) .

Conditions initiales

Elles sont prises en compte par l’intermédiaire de \({\sigma}^{-},{p}^{-},{u}^{-}\) .

En cas de poursuite ou reprise d’un calcul précédent, on a directement l’état initial \({\sigma}^{-},{p}^{-},{u}^{-}\) en partant de \(\sigma ,p,u\) du calcul précédent à l’instant spécifié.

Relation de von Mises à écrouissage isotrope

Expression des relations de comportement

Ces relations sont obtenues par les mots clés VMIS_ISOT_LINE, VMIS_ISOT_TRAC et VMIS_ISOT_PUIS.

On décrit ici ces relations en petites déformations (DEFORMATION=”PETIT”):

\(\begin{array}{c}\lbrace \begin{array}{c}{\dot{\epsilon}}^{p}=\frac{3}{2}\dot{p}.\frac{\stackrel{~}{\sigma}}{{\sigma}_{\text{eq}}}=\dot{\epsilon}-\stackrel{·}{\stackrel{⏞}{{A}^{-1}\sigma }}-{\dot{\epsilon}}^{\text{th}}\\ {\sigma}_{\text{eq}}-R\left(p\right)\le 0\\ (\begin{array}{c}\dot{p}=0\text{si}{\sigma}_{\text{eq}}-R\left(p\right)<0\\ \dot{p}\ge 0\text{si}{\sigma}_{\text{eq}}-R\left(p\right)=0\end{array}\end{array}\\ \begin{array}{cc}{\dot{\epsilon}}^{p}:& \text{vitesse de déformation plastique,}\\ p:& \text{déformation plastique cumulée}:p=\sqrt{\frac{2}{3}{\dot{\epsilon}}^{p}:{\dot{\epsilon}}^{p}}\\ {\epsilon}^{\text{th}}:& \text{déformation d'origine thermique}:{\epsilon}^{\text{th}}=\alpha \left(T-{T}_{\text{ref}}\right)\text{Id}\end{array}\end{array}\)

La fonction d’écrouissage \(R(p)\) est déduite d’un essai de traction simple monotone et isotherme.

Pour la relation de von Mises à écrouissage, la déformation plastique équivalente \({\epsilon}_{\mathit{eq}}^{p}=\sqrt{\frac{3}{2}{\stackrel{~}{\epsilon}}_{ij}^{p}{\stackrel{~}{\epsilon}}_{ij}^{p}}\) , la déformation plastique cumulée \(p\) et la variable d’écrouissage \(R(p)\) , coïncident. code_aster sauvegarde la même valeur de ces trois termes dans le variable interne «EPSPEQ» (déformation plastique équivalente).

L’utilisateur peut choisir un écrouissage linéaire (relation VMIS_ISOT_LINE) ou une courbe de traction donnée soit points par points (relation VMIS_ISOT_TRAC), soit par une expression analytique (relation

VMIS_ISOT_PUIS).

La loi de comportement VMIS_JOHN_COOK diffère des précédentes dans le sens où la fonction d’écrouissage dépend de la vitesse de la déformation plastique cumulée et de la température.

On décrit ici ces relations en petites déformations (DEFORMATION=”PETIT”):

\(\begin{array}{c}\lbrace \begin{array}{c}{\dot{\epsilon}}^{p}=\frac{3}{2}\dot{p}\frac{\stackrel{~}{\sigma}}{{\sigma}_{\text{eq}}}=\dot{\epsilon}-\stackrel{·}{\stackrel{⏞}{{A}^{-1}\sigma }}-{\dot{\epsilon}}^{\text{th}}\\ {\sigma}_{\text{eq}}-R\left(p,\dot{p},T\right)\le 0\\ \lbrace \begin{array}{c}\dot{p}=0\text{si}{\sigma}_{\text{eq}}-R\left(p,\dot{p},T\right)<0\\ \dot{p}\ge 0\text{si}{\sigma}_{\text{eq}}-R\left(p,\dot{p},T\right)=0\end{array}\end{array}\\ \begin{array}{cc}{\dot{\epsilon}}^{p}:& \text{vitesse de déformation plastique,}\\ p:& \text{déformation plastique cumulée,}\\ \dot{p}:& \text{vitesse de déformation plastique cumulée,}\\ {\epsilon}^{\text{th}}:& \text{déformation d'origine thermique}:{\epsilon}^{\text{th}}=\alpha \left(T-{T}_{\text{ref}}\right)\text{Id}\end{array}\end{array}\)

La fonction d’écrouissage \(R(p,\dot{p},T)\) est déduite d’une série d’essais de traction simple monotone à différentes vitesses de déformation et à différentes températures.

Relation VMIS_ISOT_LINE

Les données des caractéristiques de matériaux sont celles fournies sous le mot clé facteur ECRO_LINE ou ECRO_LINE_FO de l’opérateur DEFI_MATERIAU [U4.43.01].

/ ECRO_LINE =_F ( D_SIGM_EPSI = \({E}_{T}\) , SY = \({\sigma}_{y}\) )

/ ECRO_LINE_FO =_F ( D_SIGM_EPSI = \({E}_{T}\) , SY = \({\sigma}_{y}\) )

ECRO_LINE_FO correspond au cas où \({E}_{T}\) et \({\sigma}_{y}\) dépendent de la température et sont alors calculés pour la température du point de Gauss courant.

Le module d’Young \(E\) et le coefficient de Poisson \(\nu\) sont ceux fournis sous les mots clés facteurs ELAS ou ELAS_FO.

Dans ce cas la courbe de traction est la suivante:

C’est-à-dire:

\((\begin{array}{c}{\sigma}_{L}=E{\varepsilon}_{L}\text{si}{\varepsilon}_{L}\le \frac{{\sigma}_{y}}{E}\\ {\sigma}_{L}={\sigma}_{y}+{E}_{T}({\varepsilon}_{L}-\frac{{\sigma}_{y}}{E})\text{si}{\varepsilon}_{L}\ge \frac{{\sigma}_{y}}{E}\end{array}\)

Remarque:

\({\sigma}_{y}\) est la limite d’élasticité (le choix de \({\sigma}_{y}\) incombe à l’utilisateur: elle peut correspondre à la fin de linéarité de la courbe de traction réelle, ou bien à une limite d’élasticité réglementaire ou conventionnelle. Quoi qu’il en soit, on utilise ici l’unique valeur définie sous ECRO_LINE).

Lorsque le critère est atteint on a: \({\sigma}_{\mathrm{eq}}-R(p)=0\) . Pour identifier \(R(p)\) , on utilise les propriétés de l’état de contrainte uniaxial :

\(\sigma =(\begin{array}{ccc}{\sigma}_{L}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array})\) . donc \(\begin{array}{c}{\sigma}_{\text{eq}}={\sigma}_{L}\\ p={\varepsilon}_{L}^{P}={\varepsilon}_{L}-\frac{{\sigma}_{L}}{E}\end{array}\) et le critère s’écrit : \({\sigma}_{L}-R(p)=0\)

\({\sigma}_{L}-{\sigma}_{y}={E}_{T}({\varepsilon}_{L}-\frac{{\sigma}_{y}}{E})={E}_{T}(\frac{{\sigma}_{L}}{E}+p-\frac{{\sigma}_{y}}{E})\)

donc

\(({\sigma}_{L}-{\sigma}_{y})(1-\frac{{E}_{T}}{E})={E}_{T}p\) soit \(({\sigma}_{L}-{\sigma}_{y})=\frac{{E}_{T}.E}{(E-{E}_{T})}p\)

d’où la fonction d’écrouissage linéaire : \(R(p)=\frac{{E}_{T}E}{E-{E}_{T}}p+{\sigma}_{y}\)

Relation VMIS_ISOT_TRAC

Les données du matériau sont celles fournies sous le mot clé facteur TRACTION =_F (SIGM = f), de l’opérateur DEFI_MATERIAU.

f est une fonction à une ou deux variables représentant les courbes de traction simple. La première variable est obligatoirement la déformation, la deuxième si elle existe est la température (paramètre d’une nappe). Pour chaque température, la courbe de traction doit être telle que:

• les abscisses (déformations vraies) sont strictement croissantes,

• le premier point de la courbe correspond à \(({ϵ}_{1,}{\sigma}_{y})\) où \({ϵ}_{1}\ne 0\) et \({\sigma}_{y}\) est la limite d’élasticité,

• la pente entre 2 points successifs est inférieure à la pente élastique entre 0 et le premier point de la courbe \((E={\sigma}_{y}/{ϵ}_{1})\) .

Pour interpoler par rapport à la température, Code_Aster transforme d’abord toutes les courbes \(\sigma =f(ϵ)\) données par l’utilisateur en courbes \(\sigma (p)=R(p)\) de la façon suivante: si \(E\) est la pente élastique entre 0 et le premier point de la courbe \(({ϵ}_{1,}{\sigma}_{y})\) (où \({ϵ}_{1}\ne 0\) et \({\sigma}_{y}\) est la limite d’élasticité), tout point de la courbe \(({ϵ}_{i},{\sigma}_{i})\) devient le point \(({p}_{i},{\sigma}_{i})\) avec \({p}_{i}={ϵ}_{i}-{\sigma}_{i}/E\) (d’où \({p}_{1}=0\) ). Soit \(\theta\) la température considérée, s’il existe \(k\) tel que \(\theta \in [:ref:`{\theta}_{k},{\theta}_{k+1} <{\theta}_{k},{\theta}_{k+1}>\)]` où \(k\) désigne l’indice des courbes de traction contenues dans la nappe, on construit point par point la courbe \(R(p,\theta )\) en interpolant à partir de \(R(p,{\theta}_{k})\) et \(R(p,{\theta}_{k+1})\) pour toutes les valeurs de \(p\) de la réunion des valeurs des abscisses des courbes \(k\) et \(k+1\) (dans le cas où ces deux courbes sont prolongées linéairement ou par une fonction constante):

Si \({n}_{k}\) et \({n}_{k+1}\) sont les nombres de points des courbes \(k\) et \(k+1\) , le nombre de points \(n\) de la courbe \(R(p,\theta )\) vaut dans le cas général \({n}_{k}+{n}_{k+1}-1\) (cas où toutes les abscisses non nulles sont distinctes).

Si \(\theta\) est en dehors des intervalles de définition des courbes de traction, on extrapole conformément aux prolongements spécifiés par l’utilisateur dans DEFI_NAPPE [U4.31.03] et selon le principe précédent.

Remarque:

Pour éviter de générer des erreurs d’approximation importantes ou même d’obtenir par extrapolation de mauvaises courbes de traction, il vaut mieux ne pas utiliser de prolongement linéaire dans DEFI_NAPPE.

Dans le cas où le prolongement de la courbe la plus courte est «EXCLU», on arrête l’interpolation à cet endroit et le prolongement de la courbe interpolée est également «EXCLU».

On obtient donc dans tous les cas une fonction d’écrouissage linéaire par morceaux:

\(R(p,q)={s}_{i}+\frac{{s}_{i+1}-{s}_{i}}{{p}_{i+1}-{p}_{i}}(p-{p}_{i})\) pour \(p\in [:ref:`{p}_{i},{p}_{i+1} <{p}_{i},{p}_{i+1}>\)]` pour \(i+1\le n\) , avec \({p}_{1}=0\)

Le module d’Young correspondant à la température \(\theta\) est calculé de la façon suivante:

\(E={E}_{k}+\frac{q-{q}_{k}}{{q}_{k+1}-{q}_{k}}({E}_{k+1}-{E}_{k})\)

où , pour \(i=k\) ou \(i=k+1\) , \({E}_{i}\) est la pente élastique entre \(0\) et le premier point de la courbe \(\sigma =f(\epsilon )\) correspondant à la température \({\theta}_{i}\) .

C’est ce module de Young qui est utilisé dans l’intégration de la relation de comportement.

La limite d’élasticité à la température \(\theta\) vaut:

\({\sigma}_{y}=R(0,\theta )={\sigma}_{1}\)

L’utilisateur doit aussi donner le coefficient de Poisson \(\nu\) et un module d’Young fictif (qui ne sert que pour calculer la matrice de rigidité élastique si le mot clé NEWTON=_F(MATRICE=”ELASTIQUE”) est présent dans STAT_NON_LINE) par les mots clés:

/ ELAS =_F (NU = \(\nu\) , E = \(E\) )

/ ELAS_FO =_F (NU = \(\nu\) , E = \(E\) )

Relation VMIS_ISOT_PUIS

Les données du matériau sont celles fournies sous le mot clé facteur ECRO_PUIS ou ECRO_PUIS_FO de l’opérateur DEFI_MATERIAU [U4.43.01].

ECRO_PUIS=_F( SY= \({\sigma}_{y}\) , A_PUIS =a, N_PUIS =n)

La courbe d’écrouissage est déduite de la courbe uniaxiale reliant les déformations aux contraintes, dont l’expression est:

\(\epsilon =\frac{\sigma}{E}+a\frac{{\sigma}_{y}}{E}{\left(\frac{\sigma -{\sigma}_{y}}{{\sigma}_{y}}\right)}^{n}\) pour \(\sigma >{\sigma}_{y}\)

ce qui donne la courbe d’écrouissage:

\(R\left(p\right)={\sigma}_{y}+{\sigma}_{y}{\left(\frac{E}{a{\sigma}_{y}}p\right)}^{\frac{1}{n}}\)

La courbe représentative d’une telle fonction a l’allure suivante, pour différentes valeurs de n:

Les données du matériau sont celles fournies sous le mot clé facteur ECRO_COOK ou ECRO_COOK_FO de l’opérateur DEFI_MATERIAU [U4.43.01].

Attention: cette loi d’écrouissage n’est pas du tout identique à la loi de Ramberg-Osgood (souvent utilisée comme hypothèse dans les analyses simplifiées en mécanique de la rupture). La relation contraintes-déformations de la loi de Ramberg-Osgood est:

\(\epsilon =\frac{\sigma}{E}+\alpha \frac{\sigma}{E}{\left(\frac{\sigma}{{\sigma}_{y}}\right)}^{n-1}\) pour tout :math:`sigma ` .

Relation VMIS_JOHN_COOK

Les données du matériau sont celles fournies sous le mot clé facteur ECRO_COOK ou ECRO_COOK_FO de l’opérateur DEFI_MATERIAU [U4.43.01].

ECRO_COOK =_F(A=A,B=B,C=C,N_PUIS=n,M_PUIS=m,EPSP0=epsp0,TROOM=troom,TMELT=tmelt,)

La courbe d’écrouissage s’écrit de la manière suivante:

\(R\left(p,\dot{p},T\right)=\left(A+B{p}^{n}\right)\left(1+C\ln\left(\frac{\dot{p}}{\dot{{p}_{0}}}\right)\right)\left(1-{\left(\frac{T-{T}_{\mathit{room}}}{{T}_{\mathit{melt}}-{T}_{\mathit{room}}}\right)}^{m}\right)\)

ou de manière plus concise:

\(R\left(p,\dot{p},T\right)=\left(A+B{p}^{n}\right)\left(1+C{\dot{p}}^{\ast }\right)\left(1-{{T}^{\ast }}^{m}\right)\)

avec \({\dot{p}}^{\ast }=\lbrace \begin{array}{ccc}\frac{\dot{p}}{\dot{{p}_{0}}}& \mathrm{si}& \dot{p}\ge \dot{{p}_{0}}\\ 1& \mathrm{si}& \dot{p}\le \dot{{p}_{0}}\end{array}\) et \({T}^{\ast }=\lbrace \begin{array}{ccc}\frac{T-{T}_{\mathrm{room}}}{{T}_{\mathrm{melt}}-{T}_{\mathrm{room}}}& \mathrm{si}& T\ge {T}_{\mathrm{room}}\\ 0& \mathrm{si}& T\le {T}_{\mathrm{room}}\end{array}\)

Opérateur tangent. Option RIGI_MECA_TANG

Le but de ce paragraphe est de calculer l’opérateur tangent \({K}_{i-1}\) (option de calcul RIGI_MECA_TANG appelée à la première itération d’un nouvel incrément de charge) à partir des résultats connus à l’instant précédent \({t}_{i-1}\) .

Pour cela, si le tenseur des contraintes à \({t}_{i-1}\) est sur la frontière du domaine d’élasticité, on écrit la condition:

\(\dot{f}=0\)

qui doit être vérifiée (pour le problème continu en temps) conjointement à la condition:

\(f=0\)

avec:

\(f(\sigma ,p)={\sigma}_{\mathit{eq}}-R(p)\)

Si par contre le tenseur des contraintes à \({t}_{i-1}\) est à l’intérieur du domaine, \(f<0\) , alors l’opérateur tangent est l’opérateur d’élasticité.

Les quantités intervenant dans cette expression sont calculées à l’instant précédent \({t}_{i-1}\) , qui sont les seules connues au moment de la phase de prédiction. En calculant la différentielle:

\(\dot{f}=\frac{\partial f}{\partial \sigma }\dot{\sigma}+\frac{\partial f}{\partial p}\dot{p}=\frac{\partial f}{\partial \sigma }\stackrel{~}{\dot{\sigma}}+\frac{\partial f}{\partial p}\dot{p}\)

car \(\frac{\partial f}{\partial \sigma }\) est déviateur. Et donc

\(\dot{f}=\frac{\partial f}{\partial \sigma }\left(2\mu \stackrel{~}{\dot{\epsilon}}-2\mu {\dot{\epsilon}}^{P}\right)+\frac{\partial f}{\partial p}\dot{p}=\frac{\partial f}{\partial \sigma }\left(2\mu \dot{\epsilon}-2\mu {\dot{\epsilon}}^{P}\right)+\frac{\partial f}{\partial p}\dot{p}=0\)

Avec \(\sigma ={\sigma}^{-}=\sigma \left({t}_{\text{i-}1}\right)\) , \(\epsilon ={\epsilon}^{-}=\epsilon ({t}_{i-1})\) , \({\epsilon}^{p}={\epsilon}^{{p}^{-}}={\epsilon}^{p}\left({t}_{\text{i-}1}\right)\) et \(p{\text{=p}}^{-}=p\left({t}_{\text{i-}1}\right)\)

Remarque:

On ne tient pas compte dans cette expression de la variation des coefficients élastiques avec la température. C’est une approximation, sans conséquence importante, puisque cet opérateur ne sert qu’à initialiser les itérations de Newton. Par contre, la dépendance de l’opérateur tangent par rapport aux déformations thermiques est bien pris en compte au niveau de l’algorithme global [R5.03.01].

Quand on différencie par rapport à \(p\) la fonction-seuil \(f\) , on a

\(\frac{\partial f}{\partial p}=-\frac{dR}{dp}=-{R}^{'}\left(p\right)\)

On a aussi \({\dot{\epsilon}}^{p}=\frac{3}{2}\dot{p}\frac{\stackrel{~}{\sigma}}{{\sigma}_{\text{eq}}}\) et \(\frac{\partial f}{\partial \sigma }=\frac{3}{2}\frac{\stackrel{~}{\sigma}}{{\sigma}_{\text{eq}}}\) , l’expression de \(\dot{f}=0\) se ré-écrit donc

\(\frac{3}{2}\frac{\stackrel{~}{\sigma}}{{\sigma}_{\text{eq}}}\left(2\mu \dot{\epsilon}-2\mu \dot{p}\frac{3}{2}\frac{\stackrel{~}{\sigma}}{{\sigma}_{\text{eq}}}\right)-{R}^{'}\left(p\right)\dot{p}=0\)

Ce qui conduit à l’expression suivante du multiplicateur plastique:

\(\dot{p}=\frac{\left(3\mu \right)}{{\sigma}_{\text{eq}}}\frac{\left(\stackrel{~}{\sigma}.\dot{\stackrel{~}{\epsilon}}\right)}{3\mu +{R}^{'}\left(p\right)}\)

donc

\({\dot{\epsilon}}^{p}=\lbrace \begin{array}{c}\frac{9\mu }{2}\frac{\left(\stackrel{~}{\sigma}.\dot{\stackrel{~}{\epsilon}}\right)}{3\mu +{R}^{'}\left(p\right)}\frac{\stackrel{~}{\sigma}}{{\sigma}_{\text{eq}}^{2}}\text{si}f(\sigma ,p)=0\text{(plasticité)}\\ 0\text{si}f<0\text{(élasticité)}\end{array}\)

la contrainte s’écrit

\({\dot{\sigma}}_{ij}=3K{\dot{\epsilon}}_{kk}{\delta}_{ij}+2\mu \left({\dot{\stackrel{~}{\epsilon}}}_{ij}-{\dot{\epsilon}}_{ij}^{p}\right)\)

on rappelle que \(2\mu =\frac{E}{1+\nu }\) et \(3K=\frac{E}{1-2\nu }\) .

Remarque :

L’information \({\sigma}_{\text{eq}}^{-}-R\left({p}^{-}\right)=0\) est conservée sous forme d’une variable interne \(\xi\) qui vaut 1 dans ce cas et 0 si \({\sigma}_{\text{eq}}^{-}<R\left({p}^{-}\right)\) .

L’opérateur tangent lie le vecteur de déformations virtuelles \({\epsilon}^{\ast }\) à un vecteur de contraintes virtuelles \({\sigma}^{\ast }\) .

La matrice de rigidité tangente s’écrit pour un comportement élastique:

\({\sigma}^{\ast }=\left(3K\vec{1}\otimes \vec{1}+2\mu \mathrm{P}\right){\epsilon}^{\ast }\)

et pour un comportement plastique:

\({\sigma}^{\ast }=\left(3K\vec{1}\otimes \vec{1}+2\mu \mathrm{P}-{C}_{p}\mathrm{s}\otimes \mathrm{s}\right){\epsilon}^{\ast }\)

avec \(s\) le vecteur des contraintes déviatoriques associé à \({\sigma}^{-}\) défini par:

\({s}^{T}=\left({\stackrel{~}{s}}_{11}^{-},{\stackrel{~}{s}}_{22}^{-},{\stackrel{~}{s}}_{33}^{-},\sqrt{2}{\stackrel{~}{s}}_{12}^{-},\sqrt{2}{\stackrel{~}{s}}_{23}^{-},\sqrt{2}{\stackrel{~}{s}}_{31}^{-}\right)\)

et:

\({C}_{p}=\xi \frac{{\left(3\mu \right)}^{2}}{{\left({\sigma}_{\text{eq}}\right)}^{2}}\frac{1}{3\mu +{R}^{'}}\)

\(\xi =\lbrace \begin{array}{c}1\text{si}{\sigma}_{\text{eq}}^{-}-R\left({p}^{-}\right)=0\\ 0\text{sinon}\end{array}\)

Dans le cas du premier incrément de chargement, donc si l’état à l’instant précédent correspond à un état initial non contraint, l’opérateur tangent est identique à l’opérateur d’élasticité.

Calcul des contraintes et des variables internes

La décomposition des déformations permet d’écrire:

\(\Delta \varepsilon =\Delta {\varepsilon}^{p}+\Delta ({A}^{-1}\sigma )+\Delta {\varepsilon}^{\text{th}}\)

Soit, en prenant les parties sphériques et déviatoriques:

\(D\tilde{e}={\mathrm{De}}^{p}+D(\frac{\tilde{s}}{\mathrm{2m}})\) car \(\Delta {\tilde{\varepsilon}}^{\text{th}}=0.\)

\(\text{tr}\Delta \varepsilon =\Delta (\frac{\text{tr}\sigma }{\mathrm{3K}})+\text{tr}\Delta {\varepsilon}^{\text{th}}\) car \(\text{tr}\Delta {\varepsilon}^{p}=0.\)

Par discrétisation implicite directe des relations de comportement pour l’écrouissage isotrope, on obtientalors:

\(2\mu \Delta \tilde{\varepsilon}-({\tilde{\sigma}}^{-}+\Delta \tilde{\sigma})=\frac{3}{2}2\mu \Delta p\frac{{\tilde{\sigma}}^{-}+\Delta \tilde{\sigma}}{{({\tilde{\sigma}}^{-}+\Delta \tilde{\sigma})}_{\mathrm{eq}}}-2\frac{\mu}{2{\mu}^{-}}{\tilde{\sigma}}^{-}\)

\(\mathrm{tr}\sigma =\frac{\mathrm{3K}}{{\mathrm{3K}}^{-}}\mathrm{tr}{\sigma}^{-}+3K\mathrm{tr}\Delta \varepsilon -3K\mathrm{tr}{\Delta \varepsilon }^{\mathrm{th}}\)

\(\begin{array}{c}{({\sigma}^{-}+\Delta \sigma )}_{\text{eq}}-R({p}^{-}+\Delta p )\le 0\\ \Delta p =0\text{si}{({\sigma}^{-}+\Delta \sigma )}_{\text{eq}}<R({p}^{-}+\Delta p )\\ \Delta p \ge 0\text{si}{({\sigma}^{-}+\Delta \sigma )}_{\text{eq}}=R({p}^{-}+\Delta p )\end{array}\)

On définit, pour simplifier les notations, le tenseur \({\sigma}^{e}\) tel que :

\({\tilde{\sigma}}^{e}=\frac{2\mu }{2{\mu}^{-}}{\tilde{\sigma}}^{-}+2\mu \Delta \tilde{\varepsilon}\) et \(\text{tr}{\sigma}^{e}=\text{tr}\sigma\) .

Deux cas se présentent:

• \({({\sigma}^{-}+\Delta \sigma )}_{\text{eq}}<R({p}^{-}+\Delta p )\)

dans ce cas : \(\Delta p =0\text{soit}\tilde{\sigma}={\tilde{\sigma}}^{-}+\Delta \tilde{\sigma}={\tilde{\sigma}}^{e}\) , donc : \({({\tilde{\sigma}}^{e})}_{\text{eq}}<R({p}^{-})\)

• \({({\sigma}^{-}+\Delta \sigma )}_{\text{eq}}=R({p}^{-}+\Delta p )\)

dans ce cas : \(\Delta p \ge 0\) donc : \({({\tilde{\sigma}}^{e})}_{\text{eq}}\ge R({p}^{-})\)

On en déduit l’algorithme de résolution:

• si \({\tilde{\sigma}}_{\text{eq}}^{e}\le R({p}^{-})\) alors \(\Delta p=0\text{soit}\tilde{\sigma}={\tilde{\sigma}}^{-}+\Delta \tilde{\sigma}={\tilde{\sigma}}^{e}\)

• si \({\tilde{\sigma}}_{\text{eq}}^{e}>R({p}^{-})\)

alors il faut résoudre: \(\tilde{{\sigma}^{e}}=\tilde{{\sigma}^{-}}+\Delta \tilde{\sigma}+\frac{3}{2}2\mu \Delta p\frac{{\tilde{\sigma}}^{-}+\Delta \tilde{\sigma}}{{({\sigma}^{-}+\Delta \sigma )}_{\mathrm{eq}}}\)

doncen factorisant \({\tilde{\sigma}}^{-}+\Delta \tilde{\sigma}\) et en prenant la valeur équivalente de Von Mises :

\({\sigma}_{\text{eq}}^{e}=(1+\frac{3}{2}\frac{2\mu \Delta p }{{({\sigma}^{-}+\Delta \sigma )}_{\text{eq}}}){({\sigma}^{-}+\Delta \sigma )}_{\text{eq}}\)

soit:

\({\sigma}_{\text{eq}}^{e}=R({p}^{-}+\Delta p )+3\mu \Delta p\)

car: \({\sigma}_{\text{eq}}={({\sigma}^{-}+\Delta \sigma )}_{\text{eq}}=R(\overline{p}+\Delta p )\)

C’est une équation scalaire en \(\Delta p\) , linéaire ou non suivant \(R(p)\) . \(\Delta p\) est calculé de la façon suivante:

• dans le cas où l’écrouissage est linéaire (relation VMIS_ISOT_LINE), on obtient directement:

\(\Delta p =\frac{{\sigma}_{\mathrm{eq}}^{e}-{\sigma}_{y}-R'{p}^{-}}{R'+3\mu }\) \(R'=\frac{E{E}_{T}}{E-{E}_{T}}\) avec

• dans le cas où l’écrouissage est donné par une courbe de traction affine par morceaux, (relation VMIS_ISOT_TRAC), on tire parti de la linéarité par morceaux pour déterminer exactement \(\Delta p\) (voir § Annexe 1 );

• dans le cas d’un écrouissage défini par une loi en puissance (relation VMIS_ISOT_PUIS), \(\Delta p\) est solution de l’équation non linéaire: \(R({p}^{-}+\Delta p )+3\mu \Delta p -{\sigma}_{\text{eq}}^{e}=0\) . Cette équation est résolue par une méthode itérative (algorithme de type sécante). Au voisinage de l’origine, on linéarise \(R(p)\) , car la dérivée \({R}^{'}=\frac{E}{\text{an}}{(\frac{E}{a{\sigma}_{y}}p)}^{\frac{1}{n}-1}\) est infinie en \(p=0\) . Donc si \(p<{p}_{0}\) , on remplace \(R(p)\) par \({R}^{\text{lin}}(p)={\sigma}_{y}+\frac{p}{{p}_{0}}(R({p}_{0})-{\sigma}_{y})\) ,ce qui évite la recherche d’une solution numériquement presque nulle. En pratique, on choisit \({p}_{0}={10}^{-10}\) .

Une fois \(\Delta p\) déterminé, on calcule \(\sigma\) par:

\({\tilde{\sigma}}^{-}+\Delta \tilde{\sigma}=\frac{{\sigma}_{\text{eq}}^{e}-3\mu \Delta p }{{\sigma}_{\text{eq}}^{e}}.{\tilde{\sigma}}^{e}\)

et

\(\text{tr}({\sigma}^{-}+\Delta \sigma )=\text{tr}{\sigma}^{e}\) .

Les options RAPH_MECA et FULL_MECA effectuent toutes deux le calcul précédent, qui explicite le calcul de \(R({u}_{i}^{n})\) . On remarque qu’en réalité, \(R({u}_{i}^{n})={Q}^{T}{\sigma}_{i}^{n}\) où \({\sigma}_{i}^{n}\) est calculé non en fonction de \({u}_{i}^{n}\) , mais de \({\sigma}_{i-1}\text{et}{\Delta u }_{i}^{n}\) .

Remarque:

Cas particulier des contraintes planes.

Le modèle de Von Mises à écrouissage isotrope (VMIS_ISOT_LINE, VMIS_ISOT_PUIS ou VMIS_ISOT_TRAC) est également disponible en contraintes planes, c’est à dire pour les modélisations C_PLAN, DKT, COQUE_3D, COQUE_AXIS, COQUE_D_PLAN, COQUE_C_PLAN, TUYAU, TUYAU_6M..

Dans ce cas, le système à résoudre comporte une équation supplémentaire. Ce calcul est détaillé en annexe 2.

Opérateur tangent. Option FULL_MECA

L’option FULL_MECA permet de calculer la matrice tangente \({K}_{i}^{n}\) à chaque itération. L’opérateur tangent qui sert à la construire est calculé directement sur le système discrétisé précédent (on note pour simplifier : \(\tilde{\sigma}=\tilde{{\sigma}^{-}}+\Delta \tilde{\sigma},p={p}^{-}+\Delta p\) ) et on écrit les expressions seulement dans le cas isotherme.

• Si le tenseur des contraintes est sur la frontière du domaine, \(f=0\) alors on a, en différentiant l’expression de la loi de normalité en \(\tilde{\sigma}=\tilde{{\sigma}^{-}}+\Delta \tilde{\sigma}\) :

\(2\mu {\delta \varepsilon }^{p}=2\mu \delta \tilde{\varepsilon}-\delta \tilde{\sigma}=\frac{3}{2}2\mu (\deltap \frac{\tilde{\sigma}}{{\sigma}_{\text{eq}}}+\Delta p \frac{\delta \tilde{\sigma}}{{\sigma}_{\text{eq}}}-\frac{3}{2}\Delta p \frac{\tilde{\sigma}:d\tilde{\sigma}}{{\sigma}_{{\text{eq}}^{3}}}.\tilde{\sigma})\)

où \({\delta \varepsilon }^{p},\delta \tilde{\varepsilon},\delta \tilde{\sigma}\) représentent des accroissements infinitésimaux autour de la solution du problème élastoplastique incrémental obtenue précédemment.

Comme:

\(\frac{3}{2}\frac{\tilde{\sigma}:d\tilde{\sigma}}{{\sigma}_{\text{eq}}}={R}^{'}(p)\mathrm{dp}\)

en effectuant le produit tensoriel de la première équation par \(\tilde{\sigma}\) on a:

\(2\mu \tilde{\sigma}:\delta \tilde{\varepsilon}-\tilde{\sigma}:\delta \tilde{\sigma}=2\mu \deltap .{\sigma}_{\text{eq}}\)

en éliminant \(\deltap\) des deux dernières équations:

\(\tilde{\sigma}:\delta \tilde{\sigma}=\frac{2\mu \tilde{\sigma}:\delta \tilde{\varepsilon}}{1+\frac{3\mu }{{R}^{'}(p)}}\) .

• Si par contre si le tenseur des contraintes est à l’intérieur du domaine, \(f<0\) , alors l’opérateur tangent est l’opérateur d’élasticité.

En exprimant \(\deltap\) et \(\tilde{\sigma}:\delta \tilde{\sigma}\) dans la première équation, on obtient:

\(2\mu \delta \tilde{\varepsilon}-\delta \tilde{\sigma}=\frac{3\mu \Delta p }{{\sigma}_{\text{eq}}}\delta \tilde{\sigma}+{C}_{p}.{(\tilde{\sigma}:\delta \tilde{\varepsilon})}_{+}\tilde{\sigma},\)

avec :

\({C}_{p}=\frac{9{\mu}^{2}}{{\sigma}_{\mathrm{eq}}^{2}}(1-\frac{{R}^{'}(p)\Delta p}{{\sigma}_{\text{eq}}})\frac{1}{{R}^{'}(p)+3\mu }\)

La partie positive \({(\tilde{\sigma}:\delta \tilde{\varepsilon})}_{+}\) permet de regrouper en une seule équation les deux conditions:

• soit \(f<0\) , ce qui implique \(\Delta p =0\)

• soit \(f=0\)

On obtient alors:

\(\delta \tilde{\sigma}=\frac{2\mu }{a}\delta \tilde{\varepsilon}-\frac{{C}_{p}}{a}{(\tilde{\sigma}:\delta \tilde{\varepsilon})}_{+}\tilde{\sigma}\)

en posant:

\(a=1+\frac{3\mu \Delta p }{R(p+\Delta p )}\)

L’opérateur tangent lie le vecteur de déformations virtuelles \({\varepsilon}^{\ast }\) à un vecteur de contraintes virtuelles \({\sigma}^{\ast }\) .

La matrice de rigidité tangente s’écrit pour un comportement élastique:

\({\sigma}^{\ast }=(K\overrightarrow{1}\otimes \overrightarrow{1}+2\mup ){\varepsilon}^{\ast }\)

et pour un comportement plastique:

\({\sigma}^{\ast }=(K\overrightarrow{1}\otimes \overrightarrow{1}+\frac{2\mu }{a}P-\frac{{C}_{p}}{a}s\otimes s){\varepsilon}^{\ast }\)

avec \(s\) le vecteur des contraintes déviatoriques associé à \({\sigma}^{-}\) défini par:

\({s}^{T}=({\tilde{\sigma}}_{11}^{-},{\tilde{\sigma}}_{22}^{-},{\tilde{\sigma}}_{33}^{-},\sqrt{2}{\tilde{\sigma}}_{12}^{-},\sqrt{2}{\tilde{\sigma}}_{23}^{-},\sqrt{2}{\tilde{\sigma}}_{31}^{-})\)

et:

\(\xi =\lbrace \begin{array}{c}1\text{si}\Delta \varepsilon \text{conduit à une plastification}\text{et}\tilde{\sigma}.\Delta \tilde{\varepsilon}\ge 0\\ 0\text{sinon}\end{array}\)

On constate que l’opérateur tangent au système issu de la discrétisation implicite diffère de l’opérateur tangent au problème en vitesse (RIGI_MECA_TANG). On le retrouve en faisant : \(\Delta p =0\) dans les expressions de \({C}_{p}\) et \(a\) .

Variables internes des comportements VMIS_ISOT_LINE, VMIS_ISOT_PUIS, VMIS_ISOT_TRAC et VMIS_JOHN_COOK

Les relations de comportement VMIS_ISOT_LINE, VMIS_ISOT_PUIS et VMIS_ISOT_TRAC produisent deux variables internes :

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

• et \(\chi\) indicateur de plasticité à l’instant considéré (utile pour le calcul de l’opérateur tangent).

VMIS_JOHN_COOK utilise deux variables internes en plus des deux précédentes:

• \(\dot{{p}^{-}}\) vitesse de déformation plastique équivalente cumulée à l’instant moins,

• et \(\Delta {t}^{-}\) l’incrément de pas de temps à l’instant moins.

Relation de Von Mises à écrouissage cinématique linéaire

Expression de la relation de comportement, cas général

Cette relation est obtenue par le mot clé VMIS_CINE_LINE du mot clé facteur COMPORTEMENT.

Elle s’écrit(toujours en petites déformations):

\(\lbrace \begin{array}{c}{\dot{\varepsilon}}^{p}=\frac{3}{2}\dot{p}\frac{\tilde{\sigma}-\tilde{X}}{{(\sigma -X)}_{\text{eq}}}=\frac{3}{2}\dot{p}\frac{\tilde{\sigma}-X}{{(\sigma -X)}_{\text{eq}}}=\dot{\varepsilon}-\stackrel{·}{\stackrel{}{{A}^{-1}\sigma }}-{\dot{\varepsilon}}^{\text{th}}\\ X=C{\varepsilon}^{p},{\varepsilon}^{\text{th}}=\alpha (T-{T}_{\text{ref}})\mathrm{Id}\\ {(\sigma -X)}_{\text{eq}}-{\sigma}_{y}\le 0\\ \lbrace \begin{array}{c}\dot{p}=0\text{si}{(\sigma -X)}_{\text{eq}}-{\sigma}_{y}\le 0\\ \dot{p}\ge 0\text{si}{(s-X)}_{\text{eq}}-{\sigma}_{y}=0\end{array}\end{array}\) éq 4.1-1

\({\sigma}_{y}\) est la limite d’élasticité (le choix de \({\sigma}_{y}\) incombe à l’utilisateur: elle peut correspondre à la fin de linéarité de la courbe de traction réelle, ou bien à une limite d’élasticité réglementaire ou conventionnelle… Quoi qu’il en soit, on utilise ici l’unique valeur définie sous ECRO_LINE).

\(C\) est le coefficient d’écrouissage déduit des données par un essai de traction simple.

Dans ce cas(tenseur de contraintes uniaxial, tenseur de déformations plastiques isochore et orthotrope):

\(\sigma =(\begin{array}{ccc}{\sigma}_{L}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array})\) \(X=(\begin{array}{ccc}{X}_{L}& 0& 0\\ 0& -\frac{{X}_{L}}{2}& 0\\ 0& 0& -\frac{{X}_{L}}{2}\end{array})\)

\({(\sigma -X)}_{\mathrm{eq}}={\sigma}_{L}-\frac{3}{2}{X}_{L}\) et \({X}_{L}=C{\varepsilon}_{L}^{p}=C({\varepsilon}_{L}-\frac{{\sigma}_{L}}{E})\)

Les données matériau sont celles fournies sous le mot clé facteur ECRO_LINE ou ECRO_LINE_FO de l’opérateur DEFI_MATERIAU :

/ ECRO_LINE =_F( D_SIGM_EPSI = \({E}_{T}\) , SY = \({\sigma}_{y}\) )

/ ECRO_LINE_FO =_F ( D_SIGM_EPSI = \({E}_{T}\) , SY = \({\sigma}_{y}\) )

ECRO_LINE_FO correspond au cas où \({E}_{T}\) et \({\sigma}_{y}\) dépendent de la température et sont alors calculés pour la température du point de Gauss courant.

Le module d’Young \(E\) et le coefficient de Poisson sont ceux fournis sous les mots clés facteurs ELAS ou ELAS_FO.

Pour \({\varepsilon}_{L}>\frac{{\sigma}_{y}}{E}\) \({\sigma}_{L}={\sigma}_{y}+{E}_{T}({\varepsilon}_{L}-\frac{{\sigma}_{y}}{E})\) ,

mais on a également:

\(\lbrace \begin{array}{c}{\sigma}_{L}-\frac{3}{2}{X}_{L}={\sigma}_{y}\\ {X}_{L}=C({\varepsilon}_{L}-\frac{{\sigma}_{L}}{E})\end{array}\)

d’où, en éliminant \({X}_{L}\) et en identifiant:

\(C=\frac{2}{3}\frac{E{E}_{T}}{E-{E}_{T}}\) .

Expression de la relation de comportement en 1D

Pour des raisons de performances la relation est également écrite en 1D pour une utilisation avec des éléments finis de type poutre multi-fibres. Les équations précédentes sont identiques, les grandeurs \({\sigma}_{L}\) , \({X}_{L}\) et \({\varepsilon}_{L}\) sont des scalaires.

Les données matériau sont celles fournies sous le mot clé facteur ECRO_LINE de l’opérateur DEFI_MATERIAU[U4.43.01] :

/ECRO_LINE=_F(

♦D_SIGM_EPSI = \({E}_{T}\) [Réel]

♦SY = \({\sigma}_{y}\) [Réel]

◊ SIGM_ELS = sigmels [ R éel]

◊EPSI_ELU =epsielu [Réel]

)

Les opérandes sigm_els et espi_elu permettent de définir les bornes qui correspondent aux états limites de service et ultime, classiquement utilisées lors d’étude en génie civil.

◊SIGM_ELS=sgels

Définition de la contrainte limite de service.

◊EPSI_ELU=epelu

Définition de la déformation limite ultime.

Ces bornes sont obligatoires lorsque l’on utilise le comportement ecro_cine_1D (Cf. [:ref:` U4 . 5 1. 11 < U4 . 5 1. 11 >`] Comportements non-linéaires , [U4.42.07] DEFI_MATER_GC) . Dans les autres cas elles ne sont pas prises en compte.

La modélisation supportée est 1D, le nombres de variables internes est de 6.

• \(\mathrm{V1}\) :Critère ELS : CRITELS. Cette variable donne des informations par rapport à l’état limite de service. Cette variable représente la valeur absolue de la contrainte divisée par la contrainte limite à l’ELS du matériau. Si cette variable est dans \([0,1]\) le matériau respecte l’ELS.

• \(\mathrm{V2}\) :Critère ELU : CRITELU. Cette variable donne des informations par rapport à l’état limite ultime. Cette variable représente la valeur absolue de la déformation totale divisée par la déformation limite à l’ELU du matériau. Si cette variable est dans \([0,1]\) le matériau respecte l’ELU.

• \(\mathrm{V3}\) :Écrouissage cinématique : XCINXX. En 1D seul un scalaire est nécessaire.

• \(\mathrm{V4}\) :Indicateur plastique : INDIPLAS. Indique si le matériau a dépassé le critère élastique.

• \(\mathit{V5}\) :dissipation non récupérable : DISSIP. Lors de calculs sismiques il peut être utile à l’utilisateur de connaître l’énergie dissipée non récupérable. La variable dissip représente le cumul d’énergie non récupérable. L’incrément d’énergie non récupérable s’écrit sous la forme :

\(\Delta \mathit{Eg}=\frac{1}{2}({E}^{+}\Delta \varepsilon –({\sigma}^{+}-{\sigma}^{-}))\Delta \varepsilon\)

• \(\mathit{V6}\) :dissipation thermodynamique : DISSTHER . L’incrément de dissipation thermodynamique s’écrit sous la forme : \(\Delta \mathit{Eg}={\sigma}_{y}\dot{p}\) .

Opérateur tangent. Option RIGI_MECA_TANG

Le but de ce paragraphe est de calculer l’opérateur tangent \({K}_{i-1}\) (option de calcul RIGI_MECA_TANG appelée à la première itération d’un nouvel incrément de charge) à partir des résultats connus à l’instant précédent \({t}_{i-1}\) .

Pour cela, si le tenseur des contraintes à \({t}_{i-1}\) est sur la frontière du domaine d’élasticité, on écrit la condition:

\(\dot{f}=0\)

qui doit être vérifiée (pour le problème continu en temps) conjointement à la condition:

\(f=0\)

avec

\(f=f({\sigma}^{-},{X}^{-})={({\sigma}^{-}-{X}^{-})}_{\text{eq}}-{\sigma}_{y}\)

Si par contre le tenseur des contraintes à \({t}_{i-1}\) est à l’intérieur du domaine, \(f<0\) , alors l’opérateur tangent est l’opérateur d’élasticité.

On pose:

\({\sigma}^{\text{dev}}={\tilde{\sigma}}^{-}-{X}^{-}\text{et}\gamma =\lbrace \begin{array}{c}1\text{si}{({\sigma}^{-}-{X}^{-})}_{\text{eq}}-{\sigma}_{y}=0(\text{variable interne}\chi )\\ 0\text{sinon}\end{array}\)

Le problème en vitesse s’écrit dans ce cas:

\(\lbrace \begin{array}{c}{\dot{\varepsilon}}^{p}=\lbrace \begin{array}{c}\frac{1}{2\mu }\frac{3}{2}{(\frac{2\mu }{{\sigma}_{y}})}^{2}\frac{((\tilde{\sigma}-X).\dot{\tilde{\varepsilon}})(\tilde{\sigma}-X)}{C+2\mu }\text{si}(\sigma -X)-{\sigma}_{y}=0\\ 0\text{si}{(s-X)}_{\mathrm{eq}}-{\sigma}_{y}<0\end{array}\\ {\dot{\sigma}}_{ij}=K{\dot{\varepsilon}}_{kk}{\delta}_{ij}+2\mu ({\dot{\tilde{\varepsilon}}}_{ij}-{\dot{\varepsilon}}_{ij}^{p})\end{array}\)

L’opérateur tangent lie le vecteur de déformations virtuelles \({\varepsilon}^{\ast }\) à un vecteur de contraintes virtuelles \({\sigma}^{\ast }\) .

La matrice de rigidité tangente s’écrit pour un comportement élastique:

\({\sigma}^{\ast }=(K\overrightarrow{1}\otimes \overrightarrow{1}+2\mup ){\varepsilon}^{\ast }\)

et pour un comportement plastique:

\({\sigma}^{\ast }=(K\overrightarrow{1}\otimes \overrightarrow{1}+2\mup -{C}_{p}s\otimes s){\varepsilon}^{\ast }\)

avec \(s\) le vecteur des contraintes déviatoriques associé à \({\sigma}^{\mathrm{dev}}\) défini par:

\({s}^{T}=({\sigma}_{11}^{\mathrm{dev}},{\sigma}_{22}^{\mathrm{dev}},{\sigma}_{33}^{\mathrm{dev}},\sqrt{2}{\sigma}_{12}^{\mathrm{dev}},\sqrt{2}{\sigma}_{23}^{\mathrm{dev}},\sqrt{2}{\sigma}_{31}^{\mathrm{dev}})\)

et:

\({C}_{p}=\gamma \frac{3}{2}{(\frac{2\mu }{{\sigma}_{y}})}^{2}\frac{1}{2\mu +C}\)

Dans le cas du premier incrément de chargement, donc si l’état à l’instant précédent correspond à un état initial non contraint, l’opérateur tangent est identique à l’opérateur d’élasticité.

Calcul des contraintes et variables internes

La discrétisation implicite directe des relations continues conduit à résoudre:

\(\lbrace \begin{array}{c}2\mu \Delta {\varepsilon}^{p}=2\mu (\Delta \tilde{\varepsilon}+\frac{{\tilde{\sigma}}^{-}}{2{\mu}^{-}}-\frac{\tilde{\sigma}}{2\mu })=\frac{3}{2}2\mu \Delta p\frac{\tilde{\sigma}-X}{{\sigma}_{y}}\\ X=\frac{C}{{C}^{-}}{X}^{-}+C\Delta {\varepsilon}^{p}\\ {(\sigma -X)}_{\text{eq}}\le {\sigma}_{y}\\ \Delta p=0\text{si}{(\sigma -X)}_{\text{eq}}<{\sigma}_{y}\\ \Delta p\ge 0\text{sinon}\\ \text{tr}({\sigma}^{-}+\Delta \sigma )=\frac{3K}{3{K}^{-}}\text{tr}{\sigma}^{-}+3K\text{tr}\Delta \varepsilon -3K\text{tr}\Delta {\varepsilon}^{\text{th}}\end{array}\)

On pose encore:

\({\tilde{\sigma}}^{e}=\frac{2\mu }{2{\mu}^{-}}{\tilde{\sigma}}^{-}+2\mu \Delta \tilde{\varepsilon}-\frac{C}{{C}^{-}}{X}^{-}\) .

La première équation s’écrit aussi:

\((2\mu \Delta \tilde{\varepsilon}+\frac{2\mu }{2{\mu}^{-}}{\tilde{\sigma}}^{-})=\tilde{\sigma}+\frac{3}{2}2\mu \Delta p \frac{\tilde{\sigma}-X}{{\sigma}_{y}}\)

en retranchant \(X=\frac{C}{{C}^{-}}{X}^{-}+C\Delta {\varepsilon}^{p}\) a chaque terme, on obtient:

\(2\mu \Delta \tilde{\varepsilon}+\frac{2\mu }{2{\mu}^{-}}{\tilde{\sigma}}^{-}-\frac{C}{{C}^{-}}{X}^{-}=\tilde{\sigma}-X+\frac{3}{2}2\mu \Delta p \frac{\tilde{\sigma}-X}{{\sigma}_{y}}+C\Delta {\varepsilon}^{P}\)

ou encore, en utilisant la loi d’écoulement:

\({\tilde{\sigma}}^{e}=(\tilde{\sigma}-X)(\text{1+}\frac{3}{2}(2\mu +C)\frac{\Delta p }{{\sigma}_{y}})\)

On obtient encore une équation scalaire en \(\Delta p\) en prenant les valeurs équivalentes de Von Mises:

\({\sigma}_{\text{eq}}^{e}={\sigma}_{y}+\frac{3}{2}(2\mu +C)\Delta p\)

ce qui donne directement:

\(\Delta p =\frac{{\sigma}_{\text{eq}}^{e}-{\sigma}_{y}}{\frac{3}{2}(2\mu +C)}\)

Et \(\sigma\) est obtenu par: \(\tilde{\sigma}=\frac{2\mu }{2{\mu}^{-}}{\tilde{\sigma}}^{-}+2\mu \Delta \tilde{\varepsilon}+2\mu \Delta {\varepsilon}^{p}\)

En remarquant que: \(\stackrel{}{\Delta {\varepsilon}^{p}}=\frac{3}{2}\stackrel{}{\Delta p }\frac{\tilde{\sigma}-X}{{\sigma}_{y}}=\frac{3}{2}\stackrel{}{\Delta p }\frac{{\tilde{\sigma}}^{e}}{{\sigma}_{\text{eq}}^{e}}\) car : \(\frac{\tilde{\sigma}-X}{{\sigma}_{y}}=\frac{{\tilde{\sigma}}^{e}}{{\sigma}_{\text{eq}}^{e}}\)

on a donc:

\(\tilde{\sigma}=\frac{2\mu }{2{\mu}^{-}}{\tilde{\sigma}}^{-}+2\mu \stackrel{}{\Delta \tilde{\varepsilon}}-\frac{2\mu }{2\mu +C}\frac{{({s}_{\text{eq}}^{e}-{\sigma}_{y})}_{+}}{{\sigma}_{\text{eq}}^{e}}.{\tilde{\sigma}}^{e}\)

Les variables internes \(X\) sont calculées par:

\(X=\frac{C}{{C}^{-}}{X}^{-}+C\Delta {\varepsilon}^{p}=\frac{C}{{C}^{-}}{X}^{-}+\frac{3}{2}C\Delta p\frac{{\tilde{\sigma}}^{e}}{{\sigma}_{\text{eq}}^{e}}\)

Remarque: Cas particulier des contraintes planes.

La prise en compte directe de l’hypothèse des contraintes planes dans l’intégration du modèle de Von Mises à écrouissage cinématique linéaire n’a pas été faite dans Code_Aster . Pour prendre en compte cette hypothèse, c’est à dire pour utiliser un comportement élastoplastique de Von Mises avec un écrouissage cinématique linéaire (loi de Prager) avec les modélisations C_PLAN, DKT, COQUE_3D, COQUE_AXIS, COQUE_D_PLAN, COQUE_C_PLAN, TUYAU, TUYAU_6M, on peut :

• soit utiliser la méthode de condensation statique (due à R. de Borst [R5.03.03]) qui permet d’obtenir un état plan de contraintes à convergence des itérations globales de l’algorithme de Newton ;

• soit utiliser la comportement VMIS_ECMI_LINE (cf. [R5.03.16]).

Opérateur tangent. Option FULL_MECA

L’option FULL_MECA permet de calculer la matrice tangente \({K}_{i}^{n}\) à chaque itération. L’opérateur tangent qui sert à la construire est calculé directement sur le système discrétisé précédent (on note pour simplifier : \(\tilde{\sigma}=\tilde{{\sigma}^{-}}+\Delta \tilde{\sigma},p={p}^{-}+\Delta p\) ) et on écrit les expressions seulement dans le cas isotherme.

On pose \({\sigma}^{\mathrm{dev}}=\tilde{\sigma}-X\) et \(\gamma =\lbrace \begin{array}{c}1\text{si}\Delta p>0\text{et}(\tilde{\sigma}-X).\Delta \tilde{\varepsilon}\ge 0\\ 0\text{sinon}\end{array}\)

L’opérateur tangent lie le vecteur de déformations virtuelles \({\epsilon}^{\ast }\) au vecteur de contraintes virtuelles \({\sigma}^{\ast }\) . Alors la matrice de rigidité tangente s’écrit:

\({\sigma}^{\ast }=\left(K\vec{1}\otimes \vec{1}+2\mu {a}_{2}P-{C}_{p}\mathrm{s}\otimes \mathrm{s}\right){\epsilon}^{\ast }\)

avec \(s\) le vecteur de contraintes associé à \({\sigma}^{\mathrm{dev}}\) par:

\({s}^{T}=({\sigma}_{11}^{\mathrm{dev}},{\sigma}_{22}^{\mathrm{dev}},{\sigma}_{33}^{\mathrm{dev}},\sqrt{2}{\sigma}_{12}^{\mathrm{dev}},\sqrt{2}{\sigma}_{23}^{\mathrm{dev}},\sqrt{2}{\sigma}_{31}^{\mathrm{dev}})\)

et:

\({C}_{p}=\gamma \frac{3}{2}{(\frac{2\mu }{{\sigma}_{y}})}^{2}\frac{1}{2\mu +C}{a}_{1}\) avec \({a}_{1}=\frac{1}{1+\frac{3}{2}\frac{(2\mu +C)\Delta p}{{\sigma}_{y}}}\) et \({a}_{2}={a}_{1}(1+\frac{3}{2}C\frac{\Delta p }{{\sigma}_{y}})\)

Variables internes du modèle VMIS_CINE_LINE

Les variables internes sont au nombre de 7:

• le tenseur :math:`` stocké sur 6 composantes,

• la variable scalaire \(\khi\) .

Bibliographie

[bib1]

16. MIALON, Eléments d’analyse et de résolution numérique des relations de l’élastoplasticité. EDF - Bulletin de la Direction des Etudes et Recherches - Série C - N° 3 1986, p. 57 - 89.

[bib2]

E.LORENTZ, J.M.PROIX, I.VAUTIER, F.VOLDOIRE, F.WAECKEL «Initiation à la thermo‑plasticité dans le Code_Aster », Note EDF/DER/HI‑74/96/013

Description des versions du document

Relation VMIS_ISOT_TRAC: compléments sur l’intégration

La discrétisation implicite de la relation de comportement conduit à résoudre une équation en \(\Delta p\) (voir § 3.3 ):

\({\sigma}_{\mathrm{eq}}^{e}-3\mu \Delta p -R({p}^{-}+\Delta p )=0\)

On résout exactement l’équation en tirant parti de la linéarité par morceaux.

On examine d’abord si la solution pourrait être en dehors des bornes des points de discrétisation de la courbe \(R(p)\) , c’est-à-dire, si \(p\ge {p}_{n}\) est une solution possible. Pour cela:

• 

si \({\sigma}_{\mathrm{eq}}^{e}+3\mu ({p}^{-}-{p}_{n})-{\sigma}_{n}\ge 0\) , alors on est dans la situation suivante:

• si le prolongement à droite est linéaire alors: \(\Delta p=\frac{{\sigma}_{\text{eq}}^{e}-{H}_{n-1}}{{\alpha}_{n-1}+3\mu }\)

avec: \({\alpha}_{n-1}=\frac{{\sigma}_{n}-{\sigma}_{n-1}}{{p}_{n}-{p}_{n-1}}\text{}{H}_{n-1}={\sigma}_{n-1}+{\alpha}_{n-1}\left({p}^{-}-{p}_{n-1}\right)\)

• si le prolongement est constant: \(\Delta p=\frac{{\sigma}_{\mathit{eq}}^{e}-{\sigma}_{n}}{3\mu }\)

• sinon, la solution \(p\) est à chercher dans l’intervalle \([{p}_{i},{p}_{i+1}]\) tel que:

\({\sigma}_{i+1}>{\sigma}_{\text{eq}}^{e}+3\mu \left({p}^{-}-{p}_{i+1}\right)\) et \({\sigma}_{i}\le {\sigma}_{\text{eq}}^{e}+3\mu \left({p}^{-}-{p}_{i}\right)\)

alors la solution est : \(\Delta p=\frac{{\sigma}_{\text{eq}}^{e}-{H}_{i}}{{\alpha}_{i}+3\mu }\text{et}{p}^{-}+\Delta p\in \left[{p}_{i},{p}_{i+1}\right]\)

avec:

\({\alpha}_{i}=\frac{{\sigma}_{i+1}-{\sigma}_{i}}{{p}_{i+1}-{p}_{i}};{H}_{i}={\sigma}_{i}+{\alpha}_{i}\left({p}^{-}-{p}_{i}\right)\text{pour}i=1\text{à}n-1\)

Écrouissage isotrope en contraintes planes

Dans ce cas, le système à résoudre comporte une équation de plus: \(\Delta {\sigma}_{33}=0\) .On obtient alors le système suivant:

\(\lbrace \begin{array}{c}2\mu \Delta \tilde{\varepsilon}-\Delta \tilde{\sigma}=\frac{3}{2}2\mu \Delta p\frac{{\tilde{\sigma}}^{-}+\Delta \tilde{\sigma}}{{({\sigma}^{-}+\Delta \sigma )}_{\text{eq}}}\\ \text{tr}\Delta \sigma =3K\text{tr}\Delta \varepsilon \\ {({\sigma}^{-}+\Delta \sigma )}_{\text{eq}}-R({p}^{-}+\Delta p)\le 0\\ \Delta p=0\text{si}{({\sigma}^{-}+\Delta \sigma )}_{\text{eq}}<R({p}^{-}+\Delta p)\\ \begin{array}{}\Delta p\ge 0\text{si}{({\sigma}^{-}+\Delta \sigma )}_{\text{eq}}=R({p}^{-}+\Delta p)\\ \Delta {\sigma}_{33}=0\end{array}\end{array}\)

Avec cette hypothèse, \(\Delta \varepsilon\) n’est pas entièrement connu: \(\Delta {\varepsilon}_{33}\) ne peut être calculé uniquement à partir de \(\Delta {u}_{i}^{n}\) .

Remarque :

Dans le cas des modélisations autres que C_PLAN, donc par exemple pour les modélisations de coques (DKT, COQUE_3D), les hypothèses sur les termes de cisaillement transverses \(\Delta {\sigma}_{13}\) et \(\Delta {\sigma}_{23}\) sont définies par ces modélisations (en général, le comportement lié au cisaillement transverse est linéaire, élastique et découplé des équations ci-dessus). Ces termes n’entrent donc pas en ligne de compte ici.

On pose \(\Delta \varepsilon =\Delta {\varepsilon}^{q}+\Delta {\varepsilon}^{y}\) avec \(\Delta {\varepsilon}^{q}\) entièrement connu à partir de \(\Delta {u}_{i}^{n}\) et de l’élasticité, donc \(\Delta {\varepsilon}_{33}^{q}=-\frac{\nu}{1-\nu }(\Delta {\varepsilon}_{11}^{q}+\Delta {\varepsilon}_{22}^{q})\text{et}\Delta {\varepsilon}^{y}=(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& \Delta y\end{array})\) est inconnu.

Par rapport au système précédent, il y a une inconnue supplémentaire, \(\Delta y\) .

• Si \({({\tilde{\sigma}}^{-}+\Delta \tilde{\sigma})}_{\text{eq}}<R({p}^{-}+\Delta p)\text{alors}\Delta p=0\) donc \(2\mu \Delta \tilde{\varepsilon}=\Delta \tilde{\sigma},\) c’est-à-dire \(\Delta y =0.\)

• Sinon, la technique de résolution consiste à exprimer \(\Delta y\) en fonction de \(\Delta p\) . On obtient alors une équation scalaire non linéaire en \(\Delta p\) .

On pose:” \({\tilde{\sigma}}^{e}=\frac{2\mu }{2{\mu}^{-}}{\tilde{\sigma}}^{-}+2\mu \Delta {\tilde{\varepsilon}}^{q}\) . De la même façon que pour l’intégration hors contraintes planes, on obtient:

\({\tilde{\sigma}}_{e}+2\mu {\stackrel{}{\Delta \tilde{\varepsilon}}}^{y}=({\tilde{\sigma}}^{-}+\stackrel{}{\Delta \tilde{\sigma}})(1+\frac{3\mu \Delta p}{R(p+\Delta p)})\) .

Mais cette expression fait intervenir une inconnue supplémentaire \(\Delta y\) : En particulier :

\({\tilde{\sigma}}_{33}+2\mu \Delta {\tilde{\varepsilon}}_{33}^{y}=({\tilde{\sigma}}_{33}^{-}+\Delta {\tilde{\sigma}}_{33})(1+\frac{3\mu \Delta p}{R(p+\Delta p)})\)

or

\(\Delta {\tilde{\varepsilon}}_{33}^{y}=\frac{2}{3}\Delta y\)

et \(\mathrm{tr}({\sigma}^{-}+\Delta \sigma )=3K\mathrm{tr}{\Delta \varepsilon }^{q}+3{K}^{+}\Delta y+\frac{3{K}^{+}}{3{K}^{-}}\mathrm{tr}{\sigma}^{-}-3{K}^{+}\Delta {\varepsilon}^{\mathrm{th}}\)

Comme:

\({\tilde{\sigma}}_{33}^{e}+\Delta {\tilde{\sigma}}_{33}={\sigma}_{33}^{e}+\Delta {\sigma}_{33}=\frac{\mathrm{tr}({\sigma}^{-}+\Delta \sigma )}{3}=0-\frac{\mathrm{tr}({\sigma}^{-}+\Delta \sigma )}{3}\)

On obtient une équation liant \(\Delta y\) et \(\Delta p\) :

\({\tilde{\sigma}}_{33}+2\frac{2}{3}\Delta y=(1+\frac{3\mu \Delta p}{R(p+\Delta p)})(\frac{-\mathrm{tr}{\sigma}_{e}-3K\Delta y}{3})\)

avec:

\(\mathrm{tr}{\sigma}_{e}=\frac{3K}{3{K}^{-}}\mathrm{tr}{\sigma}^{-}+3K\mathrm{tr}\Delta {\varepsilon}^{q}-3K\Delta {\varepsilon}^{\mathrm{th}}\)

Soit:

\(\Delta y(\frac{4\mu }{3}+K(1+\frac{3\mu \Delta p}{R({p}^{-}+\Delta p)}))={\tilde{\sigma}}_{33}^{e}-\frac{\text{tr}{\sigma}_{e}}{3}(1+\frac{3\mu \Delta p}{R({p}^{-}+\Delta p)})\)

en remarquant que:

\({\tilde{\sigma}}_{33}^{e}={\sigma}_{33}^{e}-\frac{\text{tr}{\sigma}_{e}}{3}=0-\frac{\text{tr}{\sigma}_{e}}{3}\)

et en explicitant \(\mu ,K\) , on obtient:

\(\Delta y =\frac{3(1-2\nu )\Delta p }{E\Delta p +2(1-\nu )R(p+\Delta p )}{\tilde{\sigma}}_{33}^{e}\)

à reporter dans l’équation en \(\Delta p\) (identique aux cas précédents):

\({(\tilde{{\sigma}^{e}}+2\mu {\Delta \tilde{\varepsilon}}^{y})}_{\mathrm{eq}}-3\mu \Delta p-R({p}^{-}+\Delta p)=0\)

où \(\Delta y\) s’exprime en fonction de \(\Delta p\) puisque : \({\stackrel{}{\Delta \tilde{\varepsilon}}}^{y}=\frac{\Delta y }{3}(\begin{array}{ccc}-1& & \\ & -1& \\ & & 2\end{array})\)

L’équation scalaire en \(\Delta p\) ainsi obtenue est toujours non linéaire. Cette équation est résolue par une méthode de recherche de zéros de fonctions, basée sur un algorithme de sécante. Une fois la solution \(\Delta p\) connue on calcule \(\Delta y\) puis \(\sigma\) .