r7.02.03 Taux de restitution de l’énergie en thermo-élasticité non-linéaire#

Résumé:

On présente le calcul du taux de restitution de l’énergie par la méthode thêta en 2D ou en 3D pour un problème thermo-élastique non-linéaire. La relation de comportement élastique non linéaire est décrite en [R5.03.20].

Calcul du taux de restitution de l’énergie par la méthode thêta en grandes transformations#

On étend la relation de comportement du [§1] à de grands déplacements et de grandes rotations, dans la mesure où elle dérive d’un potentiel (loi hyperélastique). Cette fonctionnalité est déclenchée par le mot-clé DEFORMATION=”GROT_GDEP” dans la commande CALC_G.

Relation de comportement#

On désigne par:

  • \(\mathrm{E}\) le tenseur de déformations de Green-Lagrange,

  • \(\mathrm{S}\) le tenseur des contraintes de Piola‑Kirchoff 2,

  • \(\Psi (E)\) la densité d’énergie interne.

Le comportement du solide est supposé hyperélastique, à savoir que:

  • \(E\) est relié au champ de déplacement \(u\) mesuré par rapport à la configuration de référence \({\Omega}_{0}\) par:

\({E}_{ij}(u)=\frac{1}{2}({u}_{i,j}+{u}_{j,i}+{u}_{k,i}{u}_{k,j})\)

  • \(S\) est relié au tenseur des contraintes de cauchy

    ../../../../_images/Object_683.svg

par:

\({S}_{ij}=\det(F){F}_{\mathrm{ik}}^{\text{-1}}{T}_{kl}{F}_{\mathrm{jl}}^{\text{-1}}\)

\(F\) étant le gradient de la transformation qui fait passer de la configuration de référence \({\Omega}_{0}\) à la configuration actuelle \(\Omega\) , relié au déplacement par:

\({F}_{ij}=({\delta}_{ij}+{u}_{i,j})\)

La relation de comportement d’un matériau hyperélastique s’écrit sous la forme:

\({S}_{ij}=\frac{\partial \Psi }{\partial {E}_{ij}}=\frac{\partial \Psi }{\partial {E}_{jj}}={S}_{ji}\)

Cette relation décrit un comportement élastique non-linéaire, analogue à celle du [§1.1]. Elle offre la possibilité de traiter les problèmes de mécanique de la rupture sans y intégrer la plasticité. Et dans le cas d’un chargement radial monotone, elle permet d’obtenir des déformations et des contraintes de la structure semblables à celles que l’on obtiendrait si le matériau présentait un écrouissage isotrope. Le matériau hyperélastique a un comportement mécanique réversible, c’est-à-dire que tout cycle de chargement n’engendre aucune dissipation.

Ce modèle est choisi dans la commande CALC_G [U4.82.03] par l’intermédiaire du mot-clé:

RELATION: “ELAS”

pour une relation élastique « linéaire », c’est-à-dire que la relation entre les déformations et les contraintes considérées est linéaire,

RELATION: “ELAS_VMIS_LINE” ou “ELAS_VMIS_TRAC” ou “ELAS_VMIS_PUIS”

pour une relation de comportement élastique « non linéaire » (loi de Hencky-von Mises à écrouissage isotrope).

Une telle relation de comportement permet en toute rigueur de prendre en compte de grandes déformations et de grandes rotations. Toutefois, on se cantonne à de petites déformations pour assurer l’existence d’une solution et pour être identique à un comportement élastoplastique sous un chargement radial monotone [R5.03.20§2.1].

Énergie potentielle et relations d’équilibre#

Le chargement considéré se réduit à une densité surfacique non suiveuse \(R\) appliquée sur une partie \({\Gamma}_{0}\) du bord de \({\Omega}_{0}\) (hypothèse des charges mortes [R5.03.20 §2.2]).

On définit un espace des champs cinématiquement admissible \(V\) :

\(V=\lbrace v\text{admissibles,}v=0\text{sur}{\Gamma}_{0}\rbrace\)

Les relations d’équilibre en formulation faible sont:

\(\underset{{\Omega}_{0}}{\int}{F}_{\mathrm{ik}}{S}_{\mathrm{kj}}{v}_{i,j}d{\Omega}_{0}=\underset{\Gamma}{\int}{R}_{i}{v}_{i}d\Gamma\)

Elles peuvent être obtenues en minimisant l’énergie potentielle globale du système:

\(W(v)=\underset{{\Omega}_{0}}{\int}\Psi (E(v))d\Omega -\underset{\Gamma}{\int}{R}_{i}{v}_{i}d\Gamma\)

En effet, si cette fonctionnelle est minimale pour le champ de déplacement \(u\) , alors:

\(\begin{array}{cc}\delta W& =\underset{{\Omega}_{0}}{\int}\frac{\partial \Psi }{\partial {E}_{ij}}\delta {E}_{ij}d\Omega -\underset{\Gamma}{\int}{R}_{i}\delta {v}_{i}d\Gamma =0,\forall \delta v\in V\\ & =\underset{{\Omega}_{0}}{\int}{S}_{ij}\frac{1}{2}(\delta {v}_{i,j}+\delta {v}_{j,i}+\delta {v}_{p,i}{u}_{p,j}+{u}_{p,i}\delta {v}_{p,j})d\Omega -\underset{\Gamma}{\int}{R}_{i}\delta {v}_{i}d\Gamma \\ & =\underset{{\Omega}_{0}}{\int}{S}_{ij}({\delta}_{\mathrm{ip}}+{u}_{p,i})\delta {v}_{p,i}d\Omega -\underset{\Gamma}{\int}{R}_{i}\delta {v}_{i}d\Gamma \\ & =\underset{{\Omega}_{0}}{\int}{F}_{pi}{S}_{ij}\delta {v}_{p,j}d\Omega -\underset{\Gamma}{\int}{R}_{i}\delta {v}_{i}d\Gamma \\ & =\underset{{\Omega}_{0}}{\int}{F}_{\mathrm{ik}}{S}_{\mathrm{kj}}\delta {v}_{i,j}d\Omega -\underset{\Gamma}{\int}{R}_{i}\delta {v}_{i}d\Gamma =0\end{array}\)

Nous retrouvons donc les équations d’équilibre et la relation de comportement en ayant posé:

\({S}_{ij}=\frac{\partial \Psi }{\partial {E}_{ij}}\)

Expression Lagrangienne du taux de restitution de l’énergie en thermo-élasticité non-linéaire et en grandes transformations#

Par définition, le taux de restitution d’énergie \(G\) est défini par l’opposé de la dérivée de l’énergie potentielle à l’équilibre par rapport au domaine \(\Omega\) [bib1]. Il est calculé par la méthode thêta, qui est une méthode lagrangienne de dérivation de l’énergie potentielle [bib4] et [bib2]. On considère des transformations \({F}^{\eta}:M\to M+\eta \theta (M)\) du domaine en \({\Omega}_{0}\) un domaine \({\Omega}_{0}\) qui correspondent à des propagations de la fissure. A ces familles de configuration de référence ainsi définies \({\Omega}_{\eta}\) correspondent des familles de configurations déformées où la fissure s’est propagée. Le taux de restitution de l’énergie \(G\) est alors l’opposé de la dérivée de l’énergie potentielle \(W(u(\eta ))\) à l’équilibre par rapport à l’évolution initiale du fond de fissure : \(\eta\)

\(G=-{(\frac{\text{d}W(u(\eta ))}{\text{d}\eta })}_{\eta =0}\)

On note comme dans [bib4] par . la dérivation lagrangienne dans une propagation virtuelle de fissure de vitesse \(\theta\) . Soit \(\varphi (\eta ,M)\) un champ quelconque (\(\eta\) réel positif et \(M\) appartenant au domaine \({\Omega}_{0}\) ), nous noterons:

\(\stackrel{ˉ}{\varphi}(\eta ,M)=\stackrel{ˉ}{\varphi}(\eta ,{F}^{\eta}M)\text{et}\dot{\varphi}={(\frac{\partial \stackrel{ˉ}{\varphi}}{\partial \eta })}_{\eta =0}\)

L’énergie potentielle définie sur \({\Omega}_{\eta}\) est ramenée sur \({\Omega}_{0}\) ,

../../../../_images/Object_1023.svg

est supposé indépendant de \(\eta\) , la dérivation par rapport au paramètre de propagation \(\eta\) est alors aisée et le taux de restitution de l’énergie dans cette propagation est solution de l’équation variationnelle:

\({\int}_{{\Gamma}_{o}}G\theta \cdot m\text{dS}=G(\theta ),\forall \theta \in \Theta\)

avec:

\(-G(\theta )={\int}_{{\Omega}_{o}}\stackrel{.}{\stackrel{}{(\Psi (E,T))}}+\Psi (E,T){\theta}_{k,k}\text{d}\Omega -{\int}_{\Gamma}{R}_{i}\dot{{u}_{i}}+{R}_{i,k}{\theta}_{k}{u}_{i}+{R}_{i}{u}_{i}({\theta}_{k,k}-\frac{\partial \theta }{\partial {n}_{k}}{n}_{k})\text{d}\Gamma\)

Or:

\(\stackrel{.}{\stackrel{}{(\Psi (E,T))}}=\frac{\partial \Psi }{\partial {E}_{ij}}{\dot{E}}_{ij}+\frac{\partial \Psi }{\partial T}\dot{T}\)

Par la suite, nous ne considérerons que le terme \(\frac{\partial \Psi }{\partial {E}_{ij}}{\dot{E}}_{ij}\) , le terme thermique étant traité de la même façon qu’en petits déplacement [R7.02.01].

Et d’après la proposition 2 de [bib4]:

\(\begin{array}{c}{\dot{\mathrm{E}}}_{ij}=\frac{1}{2}({\dot{\mathrm{u}}}_{i,j}+{\dot{\mathrm{u}}}_{j,i}+{\dot{\mathrm{u}}}_{k,i}{\mathrm{u}}_{k,j}+{\mathrm{u}}_{k,i}{\dot{\mathrm{u}}}_{k,j})\\ \phantom{{\dot{\mathrm{E}}}_{ij}}-\frac{1}{2}({\mathrm{u}}_{i,\mathrm{p}}{\mathrm{q}}_{\mathrm{p},j}+{\mathrm{u}}_{j,\mathrm{p}}{\mathrm{q}}_{\mathrm{p},i}+{\mathrm{u}}_{k,\mathrm{p}}{\mathrm{q}}_{\mathrm{p},i}{\mathrm{u}}_{k,j}+{\mathrm{u}}_{k,i}{\mathrm{u}}_{k,\mathrm{p}}{\mathrm{q}}_{\mathrm{p},j})\end{array}\)

On peut éliminer \(\dot{u}\) de l’expression de

../../../../_images/Object_1114.svg

comme en petites déformations en remarquant que \(\dot{u}\) est cinématiquement admissible (cf. [bib3] pour les problèmes de régularité) et en utilisant l’équation d’équilibre:

\({\int}_{{\Omega}_{o}}\stackrel{.}{\stackrel{}{(\Psi (E))}}\text{d}\Omega -{\int}_{\Gamma}{R}_{i}\dot{{u}_{i}}\text{d}\Gamma =\)

\(\begin{array}{c}{\int}_{{\Omega}_{o}}\frac{\partial \Psi }{\partial {\mathrm{E}}_{ij}}\frac{1}{2}({\dot{\mathrm{u}}}_{i,j}+{\dot{\mathrm{u}}}_{j,i}+{\dot{\mathrm{u}}}_{k,i}{\mathrm{u}}_{k,j}+{\mathrm{u}}_{k,i}{\dot{\mathrm{u}}}_{k,j})d\Omega \\ -{\int}_{\Gamma}{\mathrm{R}}_{i}{\dot{\mathrm{u}}}_{i}d\Gamma -{\int}_{{\Omega}_{o}}\frac{\partial \Psi }{\partial {\mathrm{E}}_{ij}}\frac{1}{2}({\mathrm{u}}_{i,\mathrm{p}}{\theta}_{\mathrm{p},j}+{\mathrm{u}}_{j,\mathrm{p}}{\theta}_{\mathrm{p},i}+{\mathrm{u}}_{k,\mathrm{p}}{\theta}_{\mathrm{p},i}{\mathrm{u}}_{k,j}+{\mathrm{u}}_{k,i}{\mathrm{u}}_{k,\mathrm{p}}{\theta}_{\mathrm{p},j})d\Omega \\ =-{\int}_{{\Omega}_{o}}{\mathrm{S}}_{ij}({\mathrm{u}}_{i,\mathrm{p}}{\theta}_{\mathrm{p},j}+{\mathrm{u}}_{k,i}{\mathrm{u}}_{k,\mathrm{p}}{\theta}_{\mathrm{p},j})d\Omega \\ =-{\int}_{{\Omega}_{o}}{\mathrm{S}}_{ij}({\mathrm{d}}_{\text{ki}}+{\mathrm{u}}_{k,i}){\mathrm{u}}_{k,\mathrm{p}}{\theta}_{\mathrm{p},j}d\Omega \\ =-{\int}_{{\Omega}_{o}}{\mathrm{S}}_{ij}{\mathrm{F}}_{\text{ki}}{\mathrm{u}}_{k,\mathrm{p}}{\mathrm{q}}_{\mathrm{p},j}d\Omega \\ =-{\int}_{{\Omega}_{o}}{\mathrm{F}}_{\text{ik}}{\mathrm{S}}_{\text{kj}}{\mathrm{u}}_{i,\mathrm{p}}{\mathrm{q}}_{\mathrm{p},j}d\Omega \end{array}\)

Finalement, on obtient:

\(G(\theta )=\underset{{\Omega}_{0}}{\int}{\mathrm{F}}_{\mathrm{ik}}{\mathrm{S}}_{\mathrm{kj}}({\mathrm{u}}_{i,\mathrm{p}}{\theta}_{\mathrm{p},j})-\Psi (\mathrm{E}){\theta}_{k,k}\text{d}\Omega +\underset{\Gamma}{\int}{\mathrm{R}}_{i,k}{\theta}_{k}{\mathrm{u}}_{i}+{\mathrm{R}}_{i}{\mathrm{u}}_{i}({\theta}_{k,k}-\frac{\partial \theta }{\partial {\mathrm{n}}_{k}}{\mathrm{n}}_{k})\text{d}\Gamma\)

L’expression complète pour les chargements suivants:

  • densité surfacique non suiveuse \(R\) appliquée sur une partie \(\Gamma\) du bord de \({\Omega}_{0}\) ,

  • densité volumique non suiveuse \(\mathrm{f}\) appliquée sur le domaine \(\Omega\) ,

et en tenant compte de la thermique:

\(\begin{array}{cc}G(\theta )=& \underset{{\Omega}_{0}}{\int}{\mathrm{F}}_{\mathrm{ik}}{\mathrm{S}}_{\mathrm{kj}}({\mathrm{u}}_{i,\mathrm{p}}{\theta}_{\mathrm{p},j})-\Psi (\mathrm{E}){\theta}_{k,k}-\frac{\partial \Psi }{\partial \mathrm{T}}{\mathrm{T}}_{,k}{\theta}_{k}\text{d}\Omega \\ & +\underset{{\Omega}_{0}}{\int}{\mathrm{f}}_{i}{\mathrm{u}}_{i}{\theta}_{k,k}+{\mathrm{f}}_{i,k}{\theta}_{k}{\mathrm{u}}_{i}\text{d}\Omega \\ & +\underset{\Gamma}{\int}{\mathrm{R}}_{i,k}{\theta}_{k}{\mathrm{u}}_{i}+{\mathrm{R}}_{i}{\mathrm{u}}_{i}({\theta}_{k,k}-\frac{\partial \theta }{\partial {\mathrm{n}}_{k}}{\mathrm{n}}_{k})\text{d}\Gamma \end{array}\)

Implantation#

La comparaison des formules de \(G(\theta )\) du [§1.3] et du [§1.4] montre que les termes de \(G(\theta )\) sont très proches. L’introduction des grandes transformations nécessite peu de modification en post-traitement.

L’utilisation de la commande CALC_G en grandes transformations ne nécessite pas de syntaxe particulière. La présence du mot-clé DEFORMATION=”GREEN_LAGRANGE” sous le mot-clé facteur COMPORTEMENT de la commande STAT_NON_LINE indique ensuite à la commande CALC_G qu’il est nécessaire de récupérer le tenseur des contraintes de Piola-Kirchoff \(S\) et le gradient de la transformation \(F\) (routines NMGEOM et NMELNL).

Les types d’éléments finis sont les mêmes qu’en élasticité linéaire [R7.02.01 §2.4]. Ce sont les éléments isoparamétriques 2D et 3D.

Les chargements supportés sont ceux supportés en élasticité linéaire à condition que ce soient des charges mortes: typiquement une force imposée est une charge morte tandis que la pression est un chargement suiveur puisqu’il dépend de l’orientation de la surface, donc de la transformation.Cette restriction s’applique également aux conditions aux limites en déplacements.

Restriction#

Avec la relation de comportement précisée au §2, on a une formulation de \(G\) valable pour de grandes déformations pour des matériaux hyper-élastiques, mais… si l’on souhaite une cohérence avec le matériau réel qui, rappelons-le, est élasto-plastique, il est impératif de se cantonner à des déformations petites, les déplacements et les rotations pouvant être grandes.

Les conditions de chargements proportionnels et monotones, indispensables pour assurer la cohérence du modèle avec le matériau réel, conduisent à des restrictions importantes du champ des problèmes à même d’être traités par cette méthode (le thermique en particulier peut conduire à des décharges locales). Il ne peut donc s’agir que d’une solution palliative avant d’être en mesure de donner un sens au taux de restitution d’énergie dans le cadre de comportements plastiques.

Bibliographie#

  1. BUI H.D., Mécanique de la rupture fragile, Masson, 1977.

  2. DESTUYNDER Ph., DJAOUA M., Sur une interprétation de l’intégrale de Rice en théorie de la rupture fragile, Mathematics Methods in the Applied Sciences, Vol.3, pp.70-87, 1981.

  3. GRISVARD P., « Problèmes aux limites dans les polygones », Mode d’emploi- EDF- Bulletin de la Direction des Etudes et Recherches, SérieC, 1, 1986 pp.21-59.

  4. MIALON P. , « Calcul de la dérivée d’une grandeur par rapport à un fond de fissure par la méthode théta », EDF- Bulletin de la Direction des Etudes et Recherches, SérieC, n∞3, 1988, pp.1-28.

  5. MIALON P., Etude du taux de restitution de l’énergie dans une direction marquant un angle avec une fissure, note interne EDF, HI/4740-07, 1984.

  6. SIDOROFF F., Cours sur les grandes déformations, Ecole d’été, Sophia-Antipolis, 8 au 10septembre 1982.

  7. LORENTZ E., Une relation de comportement hyperélastique non linéaire, Note interne EDF DER HI-74/95/011/0, 1995.