r5.03.36 Loi de comportement KICHENIN_NL combinant élastoplasticité etviscoélasticité non linéaire#

Résumé:

Ce document décrit le modèle de comportement KICHENIN_NL, destiné à décrire la réponse mécanique du polyéthylène. Il combine une branche élastoplastique à écrouissage cinématique linéaire avec une branche viscoélastique de Maxwell non linéaire. Une modélisation des grandes déformations est accessible via la cinématique GDEF_LOG.

En préambule, précisons que cette documentation s’appuie sur la note [Lorentz, 2022]. Elle n’en retient que les éléments principaux concernant la formulation et l’intégration numérique du modèle. On renvoie le lecteur intéressé à cette note et ses références pour plus de détails le cas échéant.

Modèle continu#

Conformément à la présentation rhéologique, la contrainte \(\sigma\) est la somme de deux contributions: celle issue de la branche élastoplastique \({\sigma}^{p}\) et celle issue de la branche viscoélastique \({\sigma}^{v}\) :

(2172)#\[\sigma ={\sigma}^{p}+{\sigma}^{v}\]

Branche élastoplastique#

Il s’agit d’un modèle de plasticité à seuil de vonMises et écrouissage cinématique linéaire classique. La déformation (mécanique) \(\epsilon\) se décompose de manière additive en une partie plastique \({\epsilon}^{p}\) et une partie élastique \(\epsilon -{\epsilon}^{p}\) , cette dernière étant liée à la contrainte de la branche élastoplastique par la relation élastique suivante:

(2173)#\[{\sigma}^{p}={C}^{p}:(\epsilon -{\epsilon}^{p})=3{K}^{p}\mathit{sph}(\epsilon -{\epsilon}^{p})S+2{\mu}^{p}\mathit{dev}(\epsilon -{\epsilon}^{p})\text{};\text{}3{K}^{p}=\frac{{E}^{p}}{1-2{\nu}^{p}}\text{}2{\mu}^{p}=\frac{{E}^{p}}{1+{\nu}^{p}}\]

où \(\mathit{sph}()\) et :math:`` désignent respectivement la partie sphérique normée et le déviateur d’un tenseur d’ordre deux, tandis que \(S\) est le tenseur identité normé d’ordre deux. Pour un tenseur \(a\) quelconque d’ordre deux, cela s’exprime par:

(2174)#\[S\equiv \frac{\mathit{Id}}{\sqrt{3}};\mathit{sph}(a)\equiv S:a=\frac{\mathit{tr}(a)}{\sqrt{3}};\mathit{dev}(a)\equiv a-(S:a)S\]

A noter que l’utilisation de \(S\) et \(\mathit{sph}()\) plutôt que \(\mathit{Id}\) et \(\mathit{tr}()\) permet d’éviter la présence des \(3\) et \(1/3\) dans les différentes expressions grâce à leur caractère normé.

On introduit maintenant la fonction seuil, qui fait intervenir une contrainte de rappel \(X=C{\epsilon}^{p}\) :

(2175)#\[F(\sigma ,{\epsilon}^{p})={\Vert \sigma -C{\epsilon}^{p}\Vert }_{\mathit{eq}}-{\sigma}^{c}\text{};\text{}{\Vert a\Vert }_{\mathit{eq}}\equiv \sqrt{\frac{3}{2}\mathit{dev}(a):\mathit{dev}(a)}\]

L’évolution de la déformation plastique est gouvernée par l’équation d’écoulement. La direction d’écoulement est normale à la surface seuil (plasticité associée):

(2176)#\[{\dot{\epsilon}}^{p}=\frac{3}{2}\dot{\kappa}\frac{N}{{\Vert N\Vert }_{\mathit{eq}}}\text{};\text{}N=\mathit{dev}(\sigma -C{\epsilon}^{p})\]

Enfin, l’évolution du multiplicateur plastique \(\dot{\kappa}\) est fixée par la condition de cohérence:

(2177)#\[\dot{\kappa}\ge 0;F(\sigma ,{\epsilon}^{p})\le 0;\dot{\kappa}F(\sigma ,{\epsilon}^{p})=0\]

Branche viscoélastique#

Il s’agit d’un modèle de Maxwell dont l’amortisseur est non linéaire de type Norton. La déformation (mécanique) \(\epsilon\) se décompose de manière additive en une partie visqueuse \({\epsilon}^{v}\) et une partie élastique \(\epsilon -{\epsilon}^{v}\) , cette dernière étant liée à la contrainte de la branche viscoélastique par la relation élastique suivante:

(2178)#\[{\sigma}^{v}={C}^{v}:(\epsilon -{\epsilon}^{v})=3{K}^{v}\mathit{sph}(\epsilon -{\epsilon}^{v})S+2{\mu}^{v}\mathit{dev}(\epsilon -{\epsilon}^{v})\text{};\text{}3{K}^{v}=\frac{{E}^{v}}{1-2{\nu}^{v}}\text{}2{\mu}^{v}=\frac{{E}^{v}}{1+{\nu}^{v}}\]

L’évolution de la déformation visqueuse est gouvernée par la loi d’évolution de Norton:

(2179)#\[{\dot{\epsilon}}^{v}=\frac{1}{{\eta}^{\gamma}}{\Vert {\sigma}^{v}\Vert }_{v}^{\gamma -1}V:{\sigma}^{v};\gamma \equiv \frac{1}{n};{\Vert a\Vert }_{V}\equiv \sqrt{a:V:a};V:a\equiv (1+{\nu}^{d})\mathit{dev}(a)+(1-2{\nu}^{d})\mathit{sph}(a)S\]

où on a introduit \(\gamma\) l’exposant inverse de \(n\) pour alléger les notations et \(V\) un tenseur isotrope d’ordre 4, sans unité, qui permet de «désaligner» \({\dot{\epsilon}}^{v}\) par rapport à \({\sigma}^{v}\) .

Variables en post-traitement#

Deux variables peuvent être combinées dans certains formulations en post-traitement pour prédire la ruine du matériau. Il s’agit de la déformation plastique cumulée et de la déformation visqueuse cumulée, définies respectivement comme suit:

(2180)#\[{\dot{\epsilon}}_{\mathit{cum}}^{p}=\sqrt{\frac{2}{3}{\dot{\epsilon}}^{p}:{\dot{\epsilon}}^{p}};{\dot{\epsilon}}_{\mathit{cum}}^{v}=\sqrt{\frac{2}{3}{\dot{\epsilon}}^{v}:{\dot{\epsilon}}^{v}}\]

Pour des raisons pratiques, elles sont calculées à l’issue de la phase d’intégration du comportement et stockées parmi les variables internes.

Intégration numérique#

Équations discrétisées#

“’intégration de la loi de comportement consiste à calculerles variables internes et lacontrainte à histoire de déformation (mécanique) donnée. Compte tenu du caractère additif de la contrainte (), on peut procéder indépendamment pour chaque branche. L’opérateur tangent du système sera lui aussi la somme des opérateurs tangents propre à chaque branche.

La discrétisation en temps des équations de comportement s’appuie sur un schéma d’Euler implicite (y compris pour la branche viscoélastique), c’est-à-dire que les différentes variables du problème sont exprimées à l’instant final du pas de temps considéré. On note \({Q}_{n}\) la valeur d’une quantité \(Q\) au début du pas de temps, \(\Delta Q\) son incrément pendant le pas de temps et (simplement) \(Q\) sa valeur en fin de pas de temps. L’état mécanique au début du pas de temps \(({\epsilon}_{n},{{\epsilon}^{p}}_{n},{{\epsilon}^{v}}_{n})\) est supposé connu ainsi que l’incrément de déformation \(\Delta \epsilon\) (et donc également la déformation \(\epsilon\) ). Il s’agit alors de calculer les incréments de variables internes \(\Delta {\epsilon}^{p}\) et \(\Delta {\epsilon}^{v}\) ainsi que la contrainte à la fin du pas de temps \(\sigma\) .

Remarque: la loi étant écrite par rapport aux vraies variables thermodynamiques, elle ne prend pas en compte un état initial en contraintes défini par le mot-clef facteur ETAT_INIT. Un état initial non-nul doit donc être défini par les variables internes de déformation.

Branche élastoplastique#

La discrétisation de la relation contrainte – déformation () s’écrit comme suit:

(2181)#\[{\sigma}^{p}=3{K}^{p}\mathit{sph}(\epsilon -{{\epsilon}^{p}}_{n}-\Delta {\epsilon}^{p})S+2{\mu}^{p}\mathit{dev}(\epsilon -{{\epsilon}^{p}}_{n}-\Delta {\epsilon}^{p})\]

D’après l’équation d’écoulement (), on constate que la trace de l’incrément de déformation plastique est nul, si bien qu’il en va de même pour la déformation plastique à tout instant (déformation plastique isochore). La relation () se simplifie et on peut y faire apparaître une contrainte élastique (ou contrainte d’essai élastique) notée \({\sigma}^{e}\) et une direction d’écoulement (d’essai) élastique notée \({N}^{e}\) , qui sont connues:

(2182)#\[{\sigma}^{p}={\sigma}^{e}-2{\mu}^{p}\Delta {\epsilon}^{p};{\sigma}^{e}=3{K}^{p}\mathit{sph}(\epsilon )S+2{\mu}^{p}\mathit{dev}(\epsilon -{{\epsilon}^{p}}_{n});{N}^{e}=\mathit{dev}({\sigma}^{e})-C{\epsilon}_{n}^{p}\]

Si l’essai élastique reste à l’intérieur du domaine d’élasticité, l’incrément de déformation plastique est nul (évolution élastique):

(2183)#\[{\Vert {N}^{e}\Vert }_{\mathit{eq}}\le {\sigma}^{c}\text{}\Rightarrow \text{}\Delta {\epsilon}^{p}=0\text{et}{\sigma}^{p}={\sigma}^{e}\]

Dans le cas contraire, il faut procéder à une correction plastique (évolution plastique). L’équation d’écoulement () et l’équation de cohérence () se discrétisent alors comme suit:

(2184)#\[{\Delta \epsilon }^{p}=\frac{3}{2}\frac{\Delta \kappa }{{\sigma}^{c}}\left[{N}^{e}-(2{\mu}^{p}+C)\Delta {\epsilon}^{p}\right];{\Vert {N}^{e}-(2{\mu}^{p}+C)\Delta {\epsilon}^{p}\Vert }_{\mathit{eq}}={\sigma}^{c}\]

Après quelques manipulations algébriques, on en déduit:

(2185)#\[N={\sigma}^{c}\frac{{N}^{e}}{{\Vert {N}^{e}\Vert }_{\mathit{eq}}};\Delta {\epsilon}^{p}=\frac{{N}^{e}-N}{2{\mu}^{p}+C};{\sigma}^{p}={\sigma}^{e}-2{\mu}^{p}\Delta {\epsilon}^{p}\]

On note maintenant \({H}^{p}\) l’opérateur tangent de la branche élastoplastique, tel que \(\delta {\sigma}^{p}={H}^{p}:\delta \epsilon\) . Dans le cas d’une évolution élastique, on a simplement d’après ():

(2186)#\[{H}^{p}=3{K}^{p}S\otimes S+2{\mu}^{p}{P}^{\mathit{dev}};{P}^{\mathit{dev}}\equiv {I}_{4}-S\otimes S\]

où on a introduit \({I}_{4}\) le tenseur identité d’ordre 4 et \({P}^{\mathit{dev}}\) le projecteur sur l’espace des déviateurs. Et dans le cas d’une évolution plastique, on procède par dérivation des équations dans (). On obtient après quelques calculs:

(2187)#\[{H}^{p}=3{K}^{p}S\otimes S+\left(\frac{2{\mu}^{p}}{2{\mu}^{p}+C}\right)\left(\frac{{\sigma}^{c}}{{\Vert {N}^{e}\Vert }_{\mathit{eq}}}2{\mu}^{p}+C\right){P}^{\mathit{dev}}-3{\mu}^{p}\left(\frac{2{\mu}^{p}}{2{\mu}^{p}+C}\right)\left(\frac{{\sigma}^{c}}{{\Vert {N}^{e}\Vert }_{\mathit{eq}}}\right)\left(\frac{{N}^{e}}{{\Vert {N}^{e}\Vert }_{\mathit{eq}}}\otimes \frac{{N}^{e}}{{\Vert {N}^{e}\Vert }_{\mathit{eq}}}\right)\]

L’opérateur tangent n’est pas continu à la transition entre les régimes élastique et plastique. A priori, la connaissance du régime d’évolution et donc le choix de l’opérateur tangent dépendent de la réalisation ou non de l’inégalité \(\Vert {N}^{e}\Vert \le {\sigma}^{c}\) . Toutefois, en phase de prédiction et si le pas précédent correspondait à une évolution plastique, on peut se trouver dans le cas où \(\Vert {N}^{e}\Vert \approx {\sigma}^{c}\) , à la précision de résolution près tandis qu’on s’attend à un nouveau pas de temps en régime plastique. C’est pourquoi on choisit «d’élargir» un peu le régime plastique pour le choix de l’opérateur tangent en prédiction: si \(\Vert {N}^{e}\Vert >{\sigma}^{c}(1-\text{prec})\) , avec prec fixé à \({10}^{-6}\) , alors on choisit l’opérateur tangent plastique.

Branche viscoélastique#

La discrétisation temporelle des équations de comportement () et () conduit respectivement à:

(2188)#\[{\sigma}^{v}={C}^{v}:(\epsilon -{\epsilon}^{v})=3{K}^{v}\mathit{sph}(\epsilon -{\epsilon}^{v})S+2{\mu}^{v}\mathit{dev}(\epsilon -{\epsilon}^{v})\text{};\text{}3{K}^{v}=\frac{{E}^{v}}{1-2{\nu}^{v}}\text{}2{\mu}^{v}=\frac{{E}^{v}}{1+{\nu}^{v}}\]
(2189)#\[\frac{{\Delta \epsilon }^{v}}{\Delta t}=\frac{1}{{\eta}^{\gamma}}{\Vert {\sigma}^{v}\Vert }_{v}^{\gamma -1}V:{\sigma}^{v}\]

On introduit la contrainte instantanée \({\sigma}^{i}\) qui est une quantité connue:

(2190)#\[{\sigma}^{i}={C}^{v}:(\epsilon -{\epsilon}_{n}^{v});{\sigma}^{v}={\sigma}^{i}-{C}^{v}:\Delta {\epsilon}^{v}\]

En substituant () dans (), on obtient une équation portant sur \({\sigma}^{v}\) :

(2191)#\[\left({I}_{4}+\frac{\Delta t}{{\eta}^{\gamma}}{\Vert {\sigma}^{v}\Vert }_{v}^{\gamma -1}{C}^{v}:V\right):{\sigma}^{v}={\sigma}^{i}\]

La méthode de résolution de cette équation tensorielle s’appuie sur le caractère isotrope des tenseurs d’ordre 4 qui y apparaissent. On commence par introduire l’inconnue \(x\ge 0\) telle que:

(2192)#\[{\Vert {\sigma}^{v}\Vert }_{v}=x{\Vert {\sigma}^{i}\Vert }_{v}\]

En séparant les parties sphérique et déviatorique de (), on obtient alors respectivement:

(2193)#\[\mathit{sph}({\sigma}^{v})=\frac{\mathit{sph}({\sigma}^{i})}{1+{b}_{1}{x}^{\gamma -1}};\mathit{dev}({\sigma}^{v})=\frac{\mathit{dev}({\sigma}^{i})}{1+{b}_{2}{x}^{\gamma -1}}\]

où on a introduit les coefficients:

(2194)#\[{b}_{0}=\frac{\Delta t}{{\eta}^{\gamma}}{\Vert {\sigma}^{i}\Vert }_{v}^{\gamma -1};{b}_{1}=3{K}^{v}(1-2{\nu}^{d}){b}_{0};{b}_{2}=2{\mu}^{v}(1+{\nu}^{d}){b}_{0}\]

La connaissance de l’inconnue \(x\) permet ainsi de calculer tout le tenseur de contrainte \({\sigma}^{v}\) . En substituant ce dernier par son expression en fonction de \(x\) dans (), il vient:

(2195)#\[\sum_{k=1}^{2}\frac{{a}_{k}}{{(1+{b}_{k}{x}^{\gamma -1})}^{2}}={x}^{2};{a}_{1}=(1-2{\nu}^{d}){\left(\frac{\mathit{sph}({\sigma}^{i})}{{\Vert {\sigma}^{i}\Vert }_{v}}\right)}^{2};{a}_{2}=(1+{\nu}^{d}){\left(\frac{\mathit{dev}({\sigma}^{i})}{{\Vert {\sigma}^{i}\Vert }_{v}}\right)}^{2}\]

Il s’agit là d’une équation scalaire en \(x\) qui admet une solution unique. Elle est résolue par une méthode de Newton. Comme \(n\le 1\) et donc \(\gamma \ge 1\) , on peut déterminer des bornes:

(2196)#\[\sqrt{\sum_{k=1}^{2}\frac{{a}_{k}}{{(1+{b}_{k})}^{2}}}\le x\le 1\]

On procède aux itérations de la méthode de Newton jusqu’à obtenir un encadrement de la solution \({x}_{m}\le x\le {x}_{p}\) tel que:

(2197)#\[{x}_{p}-{x}_{m}\le \text{RESI\_INTE\_RELA}\]

On contrôle ainsi la qualité de la solution \({\sigma}^{v}\) en relatif par rapport à \(\Vert {\sigma}^{i}\Vert\) .

L’opérateur tangent \({H}^{v}\) tel que \(\delta {\sigma}^{v}={H}^{v}:\delta \epsilon\) est obtenu en différentiant successivement les équations précédentes. On en donne directement l’expression:

(2198)#\[{H}^{v}=\frac{3{K}^{v}}{1+{b}_{1}{x}^{\gamma -1}}S\otimes S+\frac{2{\mu}^{v}}{1+{b}_{2}{x}^{\gamma -1}}{P}^{\mathit{dev}}-{x}^{\gamma -1}h\otimes {B}_{0}-(\gamma -1){x}^{\gamma -1}{b}_{0}h\otimes X\]

où on a introduit les différentes quantités intermédiaires suivantes:

\({A}_{1}=2(1-2{\nu}^{d})\left[\left(\frac{3{K}^{v}\mathit{sph}({\sigma}^{i})}{{\Vert {\sigma}^{i}\Vert }_{v}}S\right)-{\left(\frac{\mathit{sph}({\sigma}^{i})}{{\Vert {\sigma}^{i}\Vert }_{v}}\right)}^{2}\left(\frac{{C}^{v}:V:{\sigma}^{i}}{{\Vert {\sigma}^{i}\Vert }_{v}}\right)\right]\)

\({A}_{2}=2(1+{\nu}^{d})\left[\left(\frac{2{\mu}^{v}\mathit{dev}({\sigma}^{i})}{{\Vert {\sigma}^{i}\Vert }_{v}}\right)-{\left(\frac{\mathit{dev}({\sigma}^{i})}{{\Vert {\sigma}^{i}\Vert }_{v}}\right)}^{2}\left(\frac{{C}^{v}:V:{\sigma}^{i}}{{\Vert {\sigma}^{i}\Vert }_{v}}\right)\right]\)

\({B}_{0}=\frac{(\gamma -1)\Delta t}{{\eta}^{\gamma}}{{\Vert {\sigma}^{i}\Vert }_{v}}^{\gamma -1}\left(\frac{{C}^{v}:V:{\sigma}^{i}}{{\Vert {\sigma}^{i}\Vert }_{v}}\right);{B}_{1}=3{K}^{v}(1-2{\nu}^{d}){B}_{0};{B}_{2}=2{\mu}^{v}(1+{\nu}^{d}){B}_{0}\)

\(X={\left[\sum_{k=1}^{2}\frac{2{a}_{k}(1+{b}_{k}\gamma {x}^{\gamma -1})}{{(1+{b}_{k}{x}^{\gamma -1})}^{3}}\right]}^{-1}\left[\sum_{k=1}^{2}\frac{(1+{b}_{k}{x}^{\gamma -1}){A}_{k}-2{a}_{k}{x}^{\gamma -1}{B}_{k}}{{(1+{b}_{k}{x}^{\gamma -1})}^{3}}\right]\)

\(h=\frac{3{K}^{v}(1-2{\nu}^{d})}{{(1+{b}_{1}{x}^{\gamma -1})}^{2}}\left(\frac{\mathit{sph}({\sigma}^{i})}{{\Vert {\sigma}^{i}\Vert }_{v}}S\right)+\frac{2{\mu}^{v}(1+{\nu}^{d})}{{(1+{b}_{2}{x}^{\gamma -1})}^{2}}\left(\frac{\mathit{dev}({\sigma}^{i})}{{\Vert {\sigma}^{i}\Vert }_{v}}\right)\)

Variables en post-traitement#

Comme mentionné précédemment, la déformation plastique cumulée et la déformation visqueuse cumulée sont calculées à l’issue de la phase d’intégration du comportement. Après discrétisation en temps, l’équation () conduit ainsi à:

(2199)#\[{\Delta \epsilon }_{\mathit{cum}}^{p}=\sqrt{\frac{2}{3}{\Delta \epsilon }^{p}:{\Delta \epsilon }^{p}};{\Delta \epsilon }_{\mathit{cum}}^{v}=\sqrt{\frac{2}{3}{\Delta \epsilon }^{v}:{\Delta \epsilon }^{v}}\]

Bibliographie#

  • Lorentz E. (2022) TUYAUTERIE 3 – Relation de comportement robuste dédiée aux calculs de structures en polyéthylène. Note interne EDF R&D 6125-1723-2022-00150.

  • Kichenin J. (1992) Comportement thermomécanique du polyéthylène. Application aux structures gazières. Rapport de thèse, Ecole Polytechnique.

    1. Kichenin, K. Dang Van, K. Boytard (1996) Finite-element simulation of a new two-dissipative mechanism model for bulk medium-density polyethylene. Journal of Materials Science 31, 1653-1661.