r5.03.08 Intégration des relations de comportement viscoélastoplastiques dans les opérateurs demécanique non linéaire

Résumé

Ce document décrit dans le cas de comportements viscoélastiques non linéaires et viscoélastoplastiques avec seuil, construites sur la base du modèle viscoélastique non linéaire de Lemaître, les ingrédients nécessaires à la mise en œuvre de l’algorithme de mécanique non linéaire, comme dansSTAT_NON_LINE décrit en [R5.03.01]. Les données d’entrée de toutes ces relations de comportement viscoélastiques non linéaires intégrées dans Code_Aster ont de manière générale la même forme. Seule la façon d’introduire la donnée principale (la fonction vitesse de déformation visqueuse) varie: elle est présentée suivant les différents mot-clés qui permettent à l’utilisateur de choisir la relation de comportement souhaitée.

Ces quantités sont calculées par une méthode d’intégration semi-implicite ou implicite. 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.

Relation continue

On se place dans l’hypothèse des petites perturbations et on scinde le tenseur des déformations en une partie élastique, une partie thermique, une partie anélastique (connue) et une partie visqueuse. Les équations sont alors:

\(\begin{array}{c}{\epsilon}_{\text{tot}}={\epsilon}_{e}+{\epsilon}_{\text{th}}+{\epsilon}_{a}+{\epsilon}_{v}\\ \sigma =\mathrm{A}\left(T\right){\epsilon}_{e}\\ {\dot{\epsilon}}_{v}=g\left({\sigma}_{\text{eq}},\lambda ,T\right)\frac{3}{2}\frac{\stackrel{~}{\sigma}}{{\sigma}_{\text{eq}}}\end{array}\)

avec:

\(\lambda\) : déformation visqueuse cumulée \(\dot{\lambda}=\sqrt{\frac{2}{3}{\dot{\epsilon}}_{v}:{\dot{\epsilon}}_{v}}\)

\(\stackrel{~}{\sigma}\) : déviateur des contraintes \(\stackrel{~}{\sigma}=\sigma -\frac{1}{3}\text{Tr}\left(\sigma \right)\mathrm{Id}\)

\({\sigma}_{\text{eq}}\) : contrainte équivalente \({\sigma}_{\text{eq}}=\sqrt{\frac{3}{2}\stackrel{~}{\sigma}:\stackrel{~}{\sigma}}\)

\(\mathrm{A}\left(T\right)\) : tenseur d’élasticité

La densité de puissance associée à la dissipation d’énergie viscoélastique s’écrit:

\({D}_{\mathit{intr}}=\stackrel{~}{\sigma}\cdot {\dot{\epsilon}}_{v}=\dot{\lambda}.{\sigma}_{\text{eq}}\)

Nature de la fonction g pour chacune des relations de comportement

Relation viscoélastique LEMAITRE

Dans ce cas, \(g\) s’exprime de façon explicite (\(\sigma\) est ici un scalaire):

\(g\left(\sigma ,\lambda ,T\right)={\left(\frac{1}{K}\frac{\sigma}{{\lambda}^{1/m}}\right)}^{n}\) avec \(\frac{1}{K}\ge 0,\frac{1}{m}\ge 0,n>0\)

Les données des caractéristiques de matériaux sont celles fournies sous les mots clés facteurs LEMAITRE ou LEMAITRE_FO de l’opérateur DEFI_MATERIAU.

LEMAITRE= _F(N=n, UN_SUR_K= \(\frac{1}{K}\) , UN_SUR_M = \(\frac{1}{m}\) )

Le module de Young \(E\) et le coefficient de Poisson \(\nu\) sont ceux fournis sous les mots clés facteurs ELAS ou ELAS_FO. Dans le cas \(1/m=0\) , cette relation de comportement tridimensionnel est de type NORTON. Si de plus \(n=1\) , cette relation de comportement se ramène à la loi viscoélastique linéaire de MAXWELL

Relation de LEMAITRE dépendant de la fluence (LEMAITRE_IRRA)

Ce paragraphe décrit la dépendance vis à vis de la fluence (et son traitement) de la relation de comportement viscoplastique de J.Lemaître, introduite pour la modélisation des assemblages combustibles et applicable aux éléments 2D et 3D massifs et aux éléments TUYAU, sous le mot-clé LEMAITRE_IRRA..

Formulation du modèle

Les équations sont les suivantes:

\(\lbrace \begin{array}{c}{\dot{\epsilon}}_{v}=\frac{3}{2}\dot{p}\frac{{\sigma}^{\mathrm{D}}}{{\sigma}_{\text{eq}}}\begin{array}{cccccc}& & & & & \end{array}\\ \dot{p}={\left[\frac{{\sigma}_{\text{eq}}}{{p}^{1/m}}\right]}^{n}{\left(\frac{1}{K}\frac{\varphi}{{\varphi}_{0}}+L\right)}^{\beta}{e}^{-\frac{Q}{R(T+{T}_{0})}}\\ \stackrel{\phantom{\rule{2em}{0ex}}\cdot \phantom{\rule{2em}{0ex}}}{\stackrel{⏞}{({\mathrm{A}}^{-1}(T).\sigma )}}={\dot{\epsilon}}_{\text{tot}}-{\dot{\epsilon}}_{v}-{\dot{\epsilon}}_{g}-{\dot{\epsilon}}_{\text{th}}\end{array}\)

avec:

\({T}_{0}=273,15°C\)

\(n>0\), \(1/K\ge 0\), \(1/m\ge 0\)

\({\phi }_{0}>0\), \(Q/R\ge 0\), \(L\ge 0\), \(\beta\) réel quelconque

Les coefficients sont fournis sous les mots-clés LEMAITRE_IRRA et ELAS de DEFI_MATERIAU et \(\phi\) est le flux neutronique (quotient de l’incrément de fluence, définie par le mot-clé AFFE_VARC de AFFE_MATERIAU, par l’incrément de temps).

La loi de grandissement est: \({\epsilon}_{g}(t)=f(T,{\Phi}_{t}){\epsilon}_{g}^{0}\) avec \({\epsilon}_{g}^{0}\) déformation uniaxiale unité dans un repère \({R}_{1}\) donné par l’utilisateur à l’aide du mot-clé MASSIF (voir [U4.42.01] et [U4.43.01]) et \(f(T,{\Phi}_{t})\) .est également une fonction définie par l’utilisateur dans DEFI_MATERIAU ([U4.43.01]).

Remarques:

La fluence, le temps et le flux \({\Phi}_{0}\) doivent être exprimés en des unités telles que le rapport \(\Phi /{\Phi}_{0}\) soit sans dimension. \(Q/R\) est en Kelvin. \(T\) est en \(°C\) .

Traitement de la dépendance vis-à-vis de la fluence

Le modèle décrit ci-dessus correspond en fait à une loi viscoélastique de Lemaître ordinaire, définie par les trois coefficients \(n\) , \(1/K’\) et \(1/m\) avec:

\(\frac{1}{K'}={\left(\frac{1}{K}\frac{\Phi}{{\Phi}_{0}}+L\right)}^{\beta /n}{e}^{-\frac{Q}{\text{nR}(T+{T}_{0})}}\)

Il suffit donc de calculer \(1/K’\) et de le fournir comme donnée au calcul à la place de \(1/K\) .

Par ailleurs, dans le calcul de la contrainte élastique, on ajoute aux déformations anélastiques (nulles par défaut) la déformation de grandissement exprimée ci-dessus, après avoir effectué le changement de repère entre le repère local et le repère \({R}_{1}\) .

Relation LEMA_SEUIL

Pour ce comportement, \(g\) s’exprime aussi de façon explicite (puisqu’il s’agit d’un cas particulier de la relation de LEMAITRE présenté ci dessus):

Si \(D\le 1\) alors \(g\left(\sigma ,\lambda ,T\right)=0\) (comportement purement élastique)

Si \(D>1\) alors \(g\left(\sigma ,\lambda ,T\right)=A\phantom{\rule{0.5em}{0ex}}\left(\frac{2}{\sqrt{3}}\sigma \right)\Phi ` avec :math:`A\ge 0\) , \(\Phi \ge 0\)

Avec: \(D=\frac{1}{S}\underset{0}{\overset{t}{\int}}{\sigma}_{\text{eq}}(u)\text{du}\)

Les données matériaux à renseigner par l’utilisateur doivent être fournies sous le mot clé LEMA_SEUIL ou LEMA_SEUIL_FO de l’opérateur DEFI_MATERIAU:

LEMA_SEUIL = _F(A=A, S=S)

Quant au paramètre \(\Phi\) , il s’agit du flux de neutron qui bombarde le matériau (quotient de l’incrément de fluence, définie par le mot-clé AFFE_VARC de AFFE_MATERIAU, par l’incrément de temps).

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

Relation VISC_IRRA_LOG

Pour cette relations, \(g\) ne s’exprime pas de façon explicite. Le comportement viscoélastoplastique est représenté par un essai de fluage unidimensionnel, à contrainte constante, qui fait intervenir le temps écoulé depuis l’instant où l’on applique la contrainte. La relation de comportement est ici définie par la donnée de quatre fonctions \({f}_{1,}{g}_{1,}{f}_{2,}{g}_{2}\) décrivant l’évolution de la déformation visqueuse au cours du temps:

\({\varepsilon}_{v}=\lambda ={f}_{1}(t){g}_{1}(\sigma ,T)+{f}_{2}(t){g}_{2}(\sigma ,T)\) éq 3.4-1

La fonction \(g\) se calcule alors numériquement en éliminant le temps \(t\) de la manière suivante:

1. pour un triplet donné \((\sigma ,\lambda ,T)\) , on résout en \(t\) l’équation [éq 3.2-1] par la méthode de Newton (voir [bib2]). On trouve une approximation de la solution \(t(\sigma ,\lambda ,T)\) ,

2. on obtient la valeur de la fonction \(g\) en \((\sigma ,\lambda ,T)\) en dérivant par rapport au temps l’équation [éq 3.2-1] (voir [bib1]):

\({\dot{\varepsilon}}_{v}=\dot{\lambda}=g(\sigma ,\lambda ,T)={f}_{1}^{'}(t){g}_{1}(\sigma ,T)+{f}_{2}^{'}(t){g}_{2}(\sigma ,T)\)

et en substituant dans cette nouvelle équation la valeur de \(t(\sigma ,\lambda ,T)\) trouvée précédemment. On trouve la formulation uniaxiale suivante:

\({\dot{\varepsilon}}_{v}=\dot{\lambda}=g(\sigma ,\lambda ,T)={f}_{1}^{'}(t(\sigma ,\lambda ,T)){g}_{1}(\sigma ,T)+{f}_{2}^{'}(t(\sigma ,\lambda ,T)){g}_{2}(\sigma ,T)\)

La forme des quatre fonctions \({f}_{1,}{g}_{1,}{f}_{2,}{g}_{2}\) est prédéfinie et l’utilisateur n’introduit que quelques paramètres dans le fichier de commande.

Pour VISC_IRRA_LOG, on a:

\({f}_{1}\left(t\right)=\ln(1+\omega .\Phi .t)\)

\({g}_{1}(\sigma ,T)=A.\sigma .{e}^{-\frac{Q}{T}}\)

\({f}_{2}(t)=\Phi .t\)

\({g}_{2}(\sigma ,T)=B.\sigma .{e}^{-\frac{Q}{T}}\)

Le paramètre \(\Phi\) , est le flux de neutrons. Il est doit être renseigné sous le mot-clé facteur AFFE_VARC via la variable de commande IRRA.

Les paramètres \(A\) , \(B\) , \(\omega\) et \(Q\) sont ceux fournis sous le mot clé facteur VISC_IRRA_LOG de l’opérateur DEFI_MATERIAU.

On notera que, pour toutes les fonctions:

\(t\) s’exprime en heures

\(T\) s’exprime en \(°C\)

\(\sigma\) s’exprime en \(\mathit{MPa}\)

Il faut donc rentrer des données cohérentes avec ces unités dans le fichier de commande et dans le fichier de maillage.

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

Intégration de la relation de comportement

Établissement de l’équation scalaire pour le schéma implicite et avec des coefficients élastiques constants

On désigne par \({\epsilon}_{\text{tot}}\) la déformation totale à l’instant \(t+\Delta t\) et par \(\Delta {\epsilon}_{\text{tot}}\) la variation de déformation totale au cours du pas de temps courant.

On appelle \({\epsilon}_{\mathrm{o}}\) la déformation à l’instant \(t+\Delta t\) résultant de la dilatation thermique et des déformations anélastiques (parmi lesquelles figurent éventuellement les déformations de grandissement, cf [§3.2]). On a donc:

\(\Delta {\epsilon}_{\mathrm{o}}=\left[\alpha \left(t+\Delta t\right)\left(T\left(t+\Delta t\right)-{T}_{\text{ref}}\right)-\alpha \left(t\right)\left(T\left(t\right)-{T}_{\text{ref}}\right)\right]\mathrm{Id}+{\epsilon}_{\mathrm{a}}\left(t+\Delta t\right)-{\epsilon}_{\mathrm{a}}\left(t\right)\)

où \(\mathrm{Id}\) est le tenseur identité d’ordre 2 en dimension 3.

On pose \(\Delta \epsilon =\Delta {\epsilon}_{\text{tot}}-\Delta {\epsilon}_{\mathrm{o}}\)

Comme on suppose ici que \(\mu\) est constant, on a la relation suivante entre les déviateurs de \(\Delta \sigma\) et \(\Delta \varepsilon\) :

\(\Delta \stackrel{~}{\sigma}=2\mu \left(\Delta \stackrel{~}{\epsilon}-\Delta {\epsilon}_{\mathrm{v}}\right)\) éq 4.1-1

Or, la loi d’écoulement s’écrit, de façon implicite:

\(\frac{\Delta {\epsilon}_{\mathrm{v}}}{\Delta t}=\frac{3}{2}g\left({\sigma}_{\text{eq}},{\lambda}^{-}+{\left(\Delta {\epsilon}_{v}\right)}_{\text{eq}},T\right)\frac{\stackrel{~}{\sigma}}{{\sigma}_{\text{eq}}}\) éq 4.1-2

On a donc, en éliminant \(\Delta {\epsilon}_{\mathrm{v}}\) entre [éq 4.1-1] et [éq 4.1-2]:

\(\begin{array}{c}2\mu \Delta \stackrel{~}{\epsilon}=\Delta \stackrel{~}{\sigma}+3\mu \Delta t.g\left({\sigma}_{\text{eq}},{\lambda}^{-}+{\left(\Delta {\epsilon}_{v}\right)}_{\text{eq}},T\right)\frac{\stackrel{~}{\sigma}}{{\sigma}_{\text{eq}}}\\ \left({\stackrel{~}{\sigma}}^{-}+2\mu \Delta \stackrel{~}{\epsilon}\right)=\left(1+3\mu \Delta t\frac{g\left({\sigma}_{\text{eq}},{\lambda}^{-}+{\left(\Delta {\epsilon}_{v}\right)}_{\text{eq}},T\right)}{{\sigma}_{\text{eq}}}\right)\stackrel{~}{\sigma}\end{array}\) éq 4.1-3

En posant \({\stackrel{~}{\sigma}}^{e}={\stackrel{~}{\sigma}}^{-}+2\mu \Delta \stackrel{~}{\epsilon}\) , on a donc:

\({\sigma}_{\text{eq}}^{e}=3\mu \Delta t.g\left({\sigma}_{\text{eq}},{\lambda}^{-}+{\left(\Delta {\epsilon}_{v}\right)}_{\text{eq}},T\right)+{\sigma}_{\text{eq}}\) éq 4.1-4

Or, on a d’après [éq 4.1-2]:

\({\left(\Delta {\epsilon}_{v}\right)}_{\text{eq}}=\Delta t.g\left({\sigma}_{\text{eq}},{\lambda}^{-}+{\left(\Delta {\epsilon}_{v}\right)}_{\text{eq}},T\right)\)

D’où:

\(\begin{array}{c}{\sigma}_{\text{eq}}^{e}=3\mu {\left(\Delta {\epsilon}_{v}\right)}_{\text{eq}}+{\sigma}_{\text{eq}}\\ {\left(\Delta {\epsilon}_{v}\right)}_{\text{eq}}=\frac{1}{3\mu }\left({\sigma}_{\text{eq}}^{e}-{\sigma}_{\text{eq}}\right)\end{array}\)

En substituant cette dernière expression dans [éq 4.1-4], on a:

\({\sigma}_{\text{eq}}^{e}=3\mu \Delta t.g\left({\sigma}_{\text{eq}},{\lambda}^{-}+\frac{1}{3\mu }\left({\sigma}_{\text{eq}}^{e}-{\sigma}_{\text{eq}}\right),T\right)+{\sigma}_{\text{eq}}\)

Si l’on pose, \({\sigma}_{\text{eq}}^{e},{\lambda}^{-},T\) et \(\Delta t\) étant connus:

\(f\left(x\right)=3\mu \Delta t.g\left(x,{\lambda}^{-}+\frac{1}{3\mu }\left({\sigma}_{\text{eq}}^{e}-x\right),T\right)+x-{\sigma}_{\text{eq}}^{e}\)

on peut alors calculer la quantité \({\sigma}_{\text{eq}}={\left({\sigma}^{-}+\Delta \sigma \right)}_{\text{eq}}\) comme étant la solution de l’équation scalaire: \(f\left(x\right)=0\) où \(x={\sigma}_{\text{eq}}\) , convention adoptée pour les paragraphes suivants.

Dans le cas d’une loi de Lemaître avec seuil, LEMA_IRRA_SEUIL, les équations précédentes ne sont utiles qu’une fois le seuil franchi. En effet en-deçà du seuil le comportement est élastique.

On discrétise le seuil implicitement:

\(D({\sigma}^{-}+\Delta \sigma )=\frac{1}{S}\underset{0}{\overset{t}{\int}}{\left({\sigma}^{-}+\Delta \sigma \right)}_{\text{eq}}(u)\text{du}\)

De la même façon que pour l’intégration des lois de comportement élasto-plastique de Von-Mises, on distingue alors deux cas.

\(D({\sigma}^{-}+\Delta \sigma )\le 1\) alors \(\phantom{\rule{0.5em}{0ex}}g\left({\sigma}^{-}+\Delta \sigma ,\lambda ,T\right)=0\)

\(D({\sigma}^{-}+\Delta \sigma )>1\) alors \(\phantom{\rule{0.5em}{0ex}}g\left({\sigma}^{-}+\Delta \sigma ,\lambda ,T\right)=A.\left(\frac{2}{\sqrt{3}}{\sigma}_{\text{eq}}\right)\phi\)

Il en résulte à partir de l’équation ci-dessus:

\(\phantom{\rule{0.5em}{0ex}}g\left({\sigma}^{-}+\Delta \sigma ,\lambda ,T\right)\ne A.\left(\frac{2}{\sqrt{3}}{\sigma}_{\text{eq}}\right)\phi\) implique \(D({\sigma}^{-}+\Delta \sigma )\le 1\)

Or \(g\) ne peut prendre que la valeur \(0\) ou \(A.(\frac{2}{\sqrt{3}}{\sigma}_{\text{eq}})\varphi\) donc

\(g({\sigma}^{-}+\Delta \sigma ,\lambda ,T)\ne A.(\frac{2}{\sqrt{3}}{\sigma}_{\text{eq}})\varphi\) implique \(D({\sigma}^{-}+\Delta \sigma )\le 1\) et \(\Delta \sigma =A.\Delta \varepsilon\)

soit \(g({\sigma}^{-}+\Delta \sigma ,\lambda ,T)=A.(\frac{2}{\sqrt{3}}{\sigma}_{\text{eq}})\varphi\) alors \(D({\sigma}^{-}+A.\Delta \varepsilon )>1\)

Le critère de franchissement du seuil peut donc s’écrire \(D({\sigma}^{-}+A.\Delta \varepsilon )>1\)

Or \(D({\sigma}^{-}+A.\Delta \varepsilon )=\frac{1}{S}\underset{0}{\overset{t}{\int}}{({\sigma}^{-}+A.\Delta \varepsilon )}_{\text{eq}}(u)\text{du}\)

En discrétisant le temps, on a alors:

\(D({\sigma}^{-}+A.\Delta \varepsilon )=\frac{1}{S}\underset{0}{\overset{{t}^{-}}{\int}}{({\sigma}^{-}+A.\Delta \varepsilon )}_{\text{eq}}(u)\text{du}+\frac{1}{S}\underset{{t}^{-}}{\overset{t}{\int}}{({\sigma}^{-}+A.\Delta \varepsilon )}_{\text{eq}}(u)\text{du}\)

\(D({\sigma}^{-}+A.\Delta \varepsilon )=\frac{1}{S}({D}^{-}S+\frac{{\sigma}_{\text{eq}}^{-}+{({\sigma}^{-}+A.\Delta \varepsilon )}_{\text{eq}}}{2}(t-{t}^{-}))\)

Résolution de l’équation scalaire : principe de la routine ZEROF2

On montre facilement que, si les conditions requises au paragraphe [§3] sur les caractéristiques des matériaux sont vérifiées, la fonction \(f\) est strictement croissante et l’équation \(f(x)=0\) admet une solution unique.

Si \({\sigma}_{\text{eq}}^{e}=0\) , alors la solution est \(x=0\) . Sinon, on a: \(f(0)=-{\sigma}_{\text{eq}}^{e}<0\)

Le problème consiste donc à trouver pour une fonction \(f\) quelconque la solution de l’équation \(f(x)=0\) sachant que cette solution existe, que \(f(0)<0\) et que \(f\) est strictement croissante.

L’algorithme adopté dans ZEROF2 est le suivant:

• on part de \({a}_{0}=0\) et \({b}_{0}={x}_{\mathrm{ap}}\) où \({x}_{\mathrm{ap}}\) est une approximation de la solution. Si c’est nécessaire (c’est-à-dire si \(f({b}_{0})<0\) ), on se ramène par la méthode des sécantes (\({z}_{n}=\frac{{a}_{n}f({b}_{n})-{b}_{n}f({a}_{n})}{f({b}_{n})-f({a}_{n})}\) puis \({a}_{n+1}={b}_{n}\) et \({b}_{n+1}={z}_{n}\) ) en une ou plusieurs itérations au cas où \(f(a)<0\) et \(f(b)>0\) :

(Dans le cas de la figure ci-dessus, cette première phrase s’est faite en une itération: \({a}_{1}={b}_{0}\) et \(f({b}_{1})>0\) ).

• on calcule \({N}_{d}\) = partie entière \((\sqrt{{N}_{\max}})\) où \({N}_{\max}\) est le nombre maximum d’itérations que l’on s’est donné. On résout alors l’équation par la méthode des sécantes en utilisant toutefois la méthode de dichotomie à chaque fois que \(n\) est multiple de \({N}_{d}\) :

Si \({N}_{d}\) divise \(n\)

\({z}_{n}=\frac{{a}_{n}+{b}_{n}}{2}\)

sinon

\({z}_{n}=\frac{{a}_{n}f({b}_{n})-{b}_{n}f({a}_{n})}{f({b}_{n})-f({a}_{n})}\)

finsi

\(n=n+1\)

si \(\mid f(z)\mid >\varepsilon\)

si \(f(z)<0\)

\({a}_{n+1}={z}_{n}\) \({b}_{n+1}={b}_{n}\)

sinon

\({a}_{n+1}={a}_{n}\) \({b}_{n+1}={z}_{n}\)

finsi

aller en 1)

sinon

La solution est: \(x={z}_{n}\to \text{FIN}\)

finsi

Cette deuxième partie de l’algorithme permet de traiter en un nombre d’itérations raisonnable les cas où \(f\) est très fortement non-linéaire, alors que la méthode des sécantes aurait convergé trop lentement. Ces cas de forte non-linéarité se rencontrent notamment avec la loi de LEMAITRE, pour des valeurs de \(\frac{n}{m}\) grandes.

Calcul de la contrainte à la fin du pas de temps courant

D’après [éq 4.1-3], si \(x\) est la solution de l’équation scalaire, en posant:

\(b(x,{\sigma}_{\text{eq}}^{e})=\frac{1}{1+3\mu \Delta t\frac{g(x,{\lambda}^{-}+\frac{1}{3\mu }({\sigma}_{\text{eq}}^{e}-x),T)}{x}}=\frac{x}{{\sigma}_{\text{eq}}^{e}}\)

on a:

\(\stackrel{~}{\sigma}=b\left(x,{\sigma}_{\text{eq}}^{e}\right){\stackrel{~}{\sigma}}^{e}\) éq 4.3-1

Dans le cas où \({\sigma}_{\text{eq}}^{e}=0\) , ce qui équivaut d’après l’équation scalaire à \(x=0\) , on prolonge \(b\) par continuité. Pour cela, on pose \(y(x)={\lambda}^{-}+\frac{1}{3\mu }({\sigma}_{\text{eq}}^{e}-x)\) et \(G(x)=g(x\text{,y}(x),T)\) . La dérivée de \(G\) s’exprime en fonction des dérivées partielles de \(g\) au point \((x\text{,y}(x),T)\) :

\(\text{G'}\left(x\right)=\frac{\partial g}{\partial x}\left(x\text{,y}\left(x\right),T\right)-\frac{1}{3\mu }\frac{\partial g}{\partial y}\left(x\text{,y}\left(x\right),T\right)\)

Le prolongement de \(b\) par continuité donne alors:

\(b\left(0,0\right)=\frac{1}{1+3\mu \Delta t\text{G'}\left(0\right)}\)

et on a, toujours dans le cas où \({\sigma}_{\text{eq}}^{e}=0\) , \(\tilde{\sigma}=0\) .

Une fois que l’on a calculé \(\stackrel{~}{\sigma}\) , on obtient \(\sigma ` par la relation (:math:`K\) est ici supposé constant):

\(\sigma ={\sigma}^{-}+\Delta \sigma =\stackrel{~}{\sigma}+\left(\frac{1}{3}\text{Tr}\left({\sigma}^{-}\right)+K.\text{Tr}\left(\Delta \epsilon \right)\right).\mathrm{Id}\) éq 4.3-2

Les variables internes sont les suivantes: \(V1\) : déformation plastique cumulée, \(V2\) : densité de dissipation d’énergie viscoélastique cumulée.

Schéma semi-implicite

Avec un schéma numérique implicite [éq 4.1-2], dans le cas, par exemple, où \(g\) ne dépend pas de \(\lambda\) , seule intervient par le calcul de \(\Delta {\epsilon}_{\mathrm{v}}\) la valeur de la contrainte en fin de pas de temps. Il peut en résulter des erreurs numériques importantes si la contrainte varie fortement au cours du temps (voir [bib2]).

Pour remédier à cela et améliorer la résolution, on discrétise la loi d’écoulement de façon semi‑implicite:

\(\frac{\Delta {\epsilon}_{\mathrm{v}}}{\Delta t}=\frac{3}{2}g\left({\left({\sigma}^{-}+\frac{\Delta \sigma }{2}\right)}_{\text{eq}},{\lambda}^{-}+\frac{{\left(\Delta {\epsilon}_{v}\right)}_{\text{eq}}}{2},{T}^{-}+\frac{\Delta T}{2}\right)\frac{\left({\stackrel{~}{\sigma}}^{-}+\frac{\Delta \stackrel{~}{\sigma}}{2}\right)}{{\left({\sigma}^{-}+\frac{\Delta \sigma }{2}\right)}_{\text{eq}}}\) éq 4.4-1

Pour transformer de la façon la plus économique ce qui a été programmé précédemment (en suivant la formulation implicite [éq 4.1-2]), il suffit de diviser chaque membre de l’équation [éq 4.4-1] par 2:

\(\frac{\left(\Delta {\epsilon}_{\mathrm{v}}/2\right)}{\Delta t}=\frac{3}{2}\frac{g\left({\left({\sigma}^{-}+\frac{\Delta \sigma }{2}\right)}_{\text{eq}},{\lambda}^{-}+\frac{{\left(\Delta {\epsilon}_{v}\right)}_{\text{eq}}}{2},{T}^{-}+\frac{\Delta T}{2}\right)}{2}\frac{\left({\stackrel{~}{\sigma}}^{-}+\frac{\Delta \stackrel{~}{\sigma}}{2}\right)}{{\left({\sigma}^{-}+\frac{\Delta \sigma }{2}\right)}_{\text{eq}}}\)

et de faire la même chose avec la relation [éq 4.1-1]:

\(\frac{\Delta \stackrel{~}{\sigma}}{2}=2\mu \left(\frac{\Delta \stackrel{~}{\epsilon}}{2}-\frac{\Delta {\epsilon}_{\mathrm{v}}}{2}\right)\)

On constate que ce système est de la même forme que celui constitué par les équations [éq 4.1-1] et [éq 4.1-2], la donnée étant \(\frac{\Delta \epsilon }{2}\) au lieu de \(\Delta \epsilon ` , les inconnues étant respectivement :math:\)frac{Delta sigma }{2}` et \(\frac{\Delta {\epsilon}_{\mathrm{v}}}{2}\) au lieu de \(\Delta \sigma ` et :math:\)Delta {epsilon}_{mathrm{v}}` et la fonction \(\frac{g}{2}\) remplaçant la fonction \(g\) .

On peut donc utiliser la résolution des paragraphes [§4.1] à [§4.3] ainsi que l’algorithme correspondant en introduisant \(\frac{\Delta \epsilon }{2}\) et en divisant la fonction \(g\) par \(2\) . Il reste alors à multiplier les résultats \(\frac{\Delta \sigma }{2}\) et \(\frac{\Delta {\epsilon}_{\mathrm{v}}}{2}\) par \(2\) pour obtenir les incréments de contrainte et de déformation visqueuse calculés par le schéma semi-implicite (le \(\Delta \sigma ` et le :math:\)Delta {epsilon}_{mathrm{v}}` de l’équation [éq 4.4-1]).

On remarquera que le calcul de l’opérateur tangent n’est pas affecté par cette modification du schéma numérique. En effet, on a évidemment:

\(\frac{\partial \Delta \sigma }{\partial \Delta \epsilon }=\frac{\partial \left(\frac{\Delta \sigma }{2}\right)}{\partial \left(\frac{\Delta \epsilon }{2}\right)}\)

Prise en compte de la variation des coefficients élastiques en fonction de la température

On a, si \(\mathrm{A}\) est le tenseur d’élasticité:

\(\Delta \epsilon =\Delta {\epsilon}_{\mathrm{v}}+\Delta \left({\mathrm{A}}^{-1}.\sigma \right)\)

avec:

\(\Delta \left({\mathrm{A}}^{-1}.\sigma \right)={\mathrm{A}}^{-1}\left({T}^{-}+\Delta T\right).\left({\sigma}^{-}+\Delta \sigma \right)-{\mathrm{A}}^{-1}\left({T}^{-}\right).{\sigma}^{-}\)

Ceci se traduit dans les équations du [§4.4] par:

\(2\mu \left(\frac{\Delta \stackrel{~}{\epsilon}}{2}\right)-\left({\stackrel{~}{\sigma}}^{-}+\frac{\Delta \stackrel{~}{\sigma}}{2}\right)=3\mu \Delta t\frac{g\left({\left({\sigma}^{-}+\frac{\Delta \sigma }{2}\right)}_{\text{eq}},{\lambda}^{-}+\frac{{\left(\Delta {\epsilon}_{v}\right)}_{\text{eq}}}{2},{T}^{-}+\frac{\Delta T}{2}\right)}{2}\frac{\left({\stackrel{~}{\sigma}}^{-}+\frac{\Delta \sigma }{2}\right)}{{\left({\sigma}^{-}+\frac{\Delta \sigma }{2}\right)}_{\text{eq}}}-{\stackrel{~}{\sigma}}^{-}.\left(\frac{2{\mu}^{-}+2\mu }{4{\mu}^{-}}\right)\)

En posant:

\({\stackrel{~}{\sigma}}^{e}=\left(\frac{2{\mu}^{-}+2\mu }{4{\mu}^{-}}\right){\stackrel{~}{\sigma}}^{-}+2\mu \left(\frac{\Delta \stackrel{~}{\epsilon}}{2}\right)\)

et

\(\text{Tr}\left({\sigma}^{\mathrm{e}}\right)=\left(\frac{3{K}^{-}+3K}{6{K}^{-}}\right)\text{Tr}\left({\sigma}^{-}\right)+3K.\text{Tr}\left(\frac{\Delta \epsilon }{2}\right)\)

on se ramène exactement au cas précédent [§4.4].

Calcul de l’opérateur tangent

Dans le cas où \({\sigma}_{\mathrm{eq}}^{e}=0\) et \(x=0\) , on prend le tenseur d’élasticité comme opérateur tangent.

Sinon, on obtient cet opérateur en dérivant l’équation [éq 4.3-1] par rapport à :math:`Delta epsilon ` :

\(\frac{\partial \stackrel{~}{\sigma}}{\partial \Delta \epsilon }=\frac{\partial \Delta \stackrel{~}{\sigma}}{\partial \Delta \epsilon }=\frac{\partial b\left(x,{\sigma}_{\text{eq}}^{e}\right)}{\partial \Delta \epsilon }{\stackrel{~}{\sigma}}^{e}+b\left(x,{\sigma}_{\text{eq}}^{e}\right)\frac{\partial {\stackrel{~}{\sigma}}^{e}}{\partial \Delta \epsilon }\)

puis en dérivant aussi [éq 4.3-2] par rapport à :math:`Delta epsilon ` :

\(\frac{\partial \Delta \sigma }{\partial \Delta \epsilon }=\frac{\partial \Delta \stackrel{~}{\sigma}}{\partial \Delta \epsilon }+K{\mathrm{I}}_{3}\frac{\partial \text{Tr}\left(\Delta \epsilon \right)}{\partial \Delta \epsilon }=\frac{\partial \Delta \stackrel{~}{\sigma}}{\partial \Delta \epsilon }+K.{\mathrm{I}}_{3}^{\mathrm{t}}{\mathrm{I}}_{3}\)

On notera que, dans ces équations, les tenseurs d’ordre 2 et d’ordre 4 sont respectivement assimilés à des vecteurs et à des matrices. \({\mathrm{I}}_{3}\) est ici un tenseur d’ordre 2, assimilé à un vecteur:

\({}^{t}\text{}{\mathrm{I}}_{3}=\left(\mathrm{1,1,1,0,0,0}\right)\)

On a de plus:

\(\frac{\partial b(x,{\sigma}_{\text{eq}}^{e})}{\partial \Delta \varepsilon }=\frac{\partial b}{\partial x}(x,{\sigma}_{\text{eq}}^{e})\frac{\partial x}{\partial \Delta \varepsilon }+\frac{\partial b}{\partial {\sigma}_{\text{eq}}^{e}}(x,{\sigma}_{\text{eq}}^{e})\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \Delta \varepsilon }\)

Il faut donc calculer \(\frac{\partial x}{\partial \Delta \varepsilon }\) . Pour cela, on dérive implicitement l’équation scalaire par rapport à \(\Delta \varepsilon\) .

Pour simplifier, on omettra par la suite dans l’écriture de \(g\) et de ses dérivées le paramètre \(T\) .

On a alors:

\(\left[3\mu \Delta t{G}^{'}(x)+1\right]\frac{\partial x}{\partial \Delta \varepsilon }+\Delta t\frac{\partial g}{\partial y}(x,y)\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \Delta \varepsilon }=\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \Delta \varepsilon }\)

D’où:

\(\frac{\partial x}{\partial \Delta \varepsilon }=\frac{1-\Delta t\frac{\partial g}{\partial y}(x,y)}{1+3\mu \Delta t{G}^{'}(x)}\frac{\partial {\sigma}_{\text{eq}}^{e}}{\partial \Delta \varepsilon }\)

\(\frac{\partial x}{\partial \Delta \varepsilon }=\frac{1-\Delta t\frac{\partial g}{\partial y}(x,y)}{1+3\mu \Delta t{G}^{'}(x)}\frac{3\mu }{{\sigma}_{\mathrm{eq}}^{e}}{\tilde{\sigma}}^{e}\)

avec l’expression de \(\text{G'}(x)\) obtenue au [§4.3].

On obtient finalement l’expression suivante de l’opérateur tangent:

\(\frac{\partial \Delta \sigma }{\partial \Delta \varepsilon }=K{I}_{3}^{t}{I}_{3}+2\mu \left[\gamma {\tilde{\sigma}}^{e}{}^{t}\text{}{\tilde{\sigma}}^{e}+b(x,{\sigma}_{\text{eq}}^{e})A\right]\)

avec

\(A=\frac{\partial \Delta \tilde{\varepsilon}}{\partial \Delta \varepsilon }={J}_{6}-\frac{1}{3}{I}_{3}{}^{t}\text{}{I}_{3}\) où \({J}_{6}\) est la matrice identité de rang 6.

\(\gamma =\frac{3}{2{({\sigma}_{\mathrm{eq}}^{e})}^{3}}\left[{\sigma}_{\mathrm{eq}}^{e}\frac{1-\Delta t\frac{\partial g}{\partial y}(x,y)}{1+3\mu \Delta t{G}^{'}(x)}-x\right]\)

Remarque:

Dans le cas de la loi VISC_IRRA_LOG, on vérifie facilement que:

\(\begin{array}{}{G}^{'}(x)=\frac{1}{{f}_{1}^{'}{g}_{1}+{f}_{2}^{'}{g}_{2}}[{g}_{1}{g}_{1}^{'}({f}_{1}^{{'}^{2}}-{f}_{1}{f}_{1}^{\text{''}})+{g}_{2}{g}_{2}^{'}({f}_{2}^{{'}^{2}}-{f}_{2}{f}_{2}^{\text{''}})+{g}_{1}{g}_{2}^{'}({f}_{1}^{'}{f}_{2}^{'}-{f}_{1}^{\text{''}}{f}_{2})\\ +{g}_{2}{g}_{1}^{'}({f}_{1}^{'}{f}_{2}^{'}-{f}_{1}{f}_{2}^{\text{''}})-\frac{1}{\mathrm{3m}}({f}_{1}^{\text{''}}{g}_{1}+{f}_{2}^{\text{''}}{g}_{2})]\\ \frac{\partial g}{\partial \lambda }(x,y,T)=\frac{{f}_{1}^{\text{''}}{g}_{1}+{f}_{2}^{\text{''}}{g}_{2}}{{f}_{1}^{'}{g}_{1}+{f}_{2}^{'}{g}_{2}}\end{array}\)

où \({f}_{1,}f{'}_{1,}f'{'}_{1,}{f}_{2,}f{'}_{2,}f'{'}_{2}\) désignent les valeurs de \({f}_{1}\) et \({f}_{2}\) et de leurs dérivées au point \(t(x,y,T)\) et où \({g}_{1,}g{'}_{1,}{g}_{2,}g{'}_{2}\) désignent les valeurs de \({g}_{1}\) et \({g}_{2}\) et de leur dérivée par rapport à \(s\) au point \((x,T)\) (voir [bib1]).

Bibliographie

1. de BONNIERES P.: Ecriture sous forme standard généralisée des lois de comportement viscoplastiques du Zircaloy, note EDF-DER HI-71/7940-Indice A, 1992

2. de BONNIERES P., ZIDI M.: Introduction de la viscoplasticité dans le module de thermomécanique de Cyrano3: principe, description et validation, note EDF-DER HI-71/8334, 1993.

Description des versions du document

Version Aster	Auteur(s) Organisme(s)	Description des modifications
7.4	de BONNIERES	Texte initial
8.4	LECLERCQ, R.MASSON	Lois GATT-MONERIEet LEMA_SEUIL
9.3	de BONNIERES	Suppression ZIRC_CYRA2, ZIRC_EPRI.
16	F.VOLDOIRE	Ajout de la dissipation d’énergie visqueuse