2. Formulation variationnelle mixte du problème

2.1. Formulation dans le cadre des petites déformations

Soit un solide \(\Omega\) soumis à:

• un champ de déplacement imposé \(u={u}_{0}\) sur \({\Gamma}_{u}\)

• un champ de contrainte imposée \(t=\sigma .n={t}_{0}\) sur \({\Gamma}_{t}\)

• un champ d’effort volumique \(f\) sur \(\Omega\)

Dans le cas classique des éléments finis en déplacement (modélisation 3D ou D_PLAN ou AXIS dans code_aster ), lorsque le problème dérive d’une énergie, le problème résolu est le suivant:

trouver \(u\in V\) avec \(\sigma\) vérifiant la relation de comportement, qui minimise l’énergie potentielle:

\[\Pi (u)=\frac{1}{2}\underset{\Omega}{\int}\sigma :\varepsilon d\Omega -\underset{\Omega}{\int}f.ud\Omega -\underset{{\Gamma}_{t}}{\int}t.ud\Gamma\]

Comme nous l’avons expliqué dans la section [Difficultés liées au traitement de l’incompressibilité], cette formulation ne convient pas lorsqu’on cherche à se rapprocher de la solution incompressible, c’est-à-dire de la condition \(div(u)=0\) ou \(\text{tr}(\varepsilon )=0\) . Pour contourner cette difficulté, une solution est de traiter séparément la partie sphérique du tenseur des déformations (la partie qui pose des problèmes numériques) et sa partie déviatorique. On aura donc:

(2.1)\[\varepsilon (u,g)={\varepsilon}^{D}(u)+\frac{g}{3}I\]

Le problème précédent se ramène donc à la résolution d’un problème à 2 variables, \(u\) et \(g\) , sous la contrainte \(g=\text{tr}(\varepsilon (u))\) . Il peut être ramené à la résolution d’un problème sans contrainte en introduisant un multiplicateur de Lagrange \(p\) ; il s’écrit:

trouver \(u\in V\) , \(p\) et \(g\) solutions du problème de point-selle pour le lagrangien:

(2.2)\[ℒ(u,p,g)=\underset{\Omega}{\int}\left[\sigma :({\varepsilon}^{D}(u)+\frac{g}{3}I)+p(div(u)-g)\right]d\Omega -\underset{\Omega}{\int}f.ud\Omega -\underset{{\Gamma}_{t}}{\int}t.ud\Gamma\]

Ce problème peut être résolu, en écrivant les conditions d’optimalité :

(2.3)\[\begin{split}\left\lbrace \begin{array}{c}\frac{\partial \mathrm{ℒ}}{\partial u}=\underset{\Omega}{\int}({\sigma}^{D}+pI):\delta \varepsilon d\Omega -\underset{\Omega}{\int}f.\delta ud\Omega -\underset{{\Gamma}_{t}}{\int}t.\delta ud\Gamma =0\\ \frac{\partial \mathrm{ℒ}}{\partial p}=\underset{\Omega}{\int}(div(u)-g)\delta pd\Omega =0\\ \frac{\partial \mathrm{ℒ}}{\partial g}=\underset{\Omega}{\int}(\frac{1}{3}\text{tr}(\sigma )-p)\delta gd\Omega =0\end{array}\right.\end{split}\]

• la première équation correspond à l’équation d’équilibre,

• la deuxième équation traduit la relation cinématique liant \(g\) à \(u\),

• la troisième équation donne l’expression du multiplicateur de Lagrange \(p\) ,

lorsque le problème ne dérive pas d’une énergie, on peut directement utiliser le système d’équations (2.3).

Dans le cas où il existe une relation biunivoque entre la pression et le gonflement comme par exemple pour un matériau élastoplastique avec un critère de plasticité de type von Mises, il est possible d’expliciter le gonflement et donc de supprimer la troisième équation du système (2.3). On obtient alors le système de deux équations à deux inconnues qui suit:

(2.4)\[\begin{split}\left\lbrace \begin{array}{}\underset{\Omega}{\int}({\sigma}^{D}+pI):\delta \varepsilon d\Omega -\underset{\Omega}{\int}f.\delta ud\Omega -\underset{{\Gamma}_{t}}{\int}t.\delta ud\Gamma =0\\ \underset{\Omega}{\int}(div(u)-\frac{p}{\kappa})\delta pd\Omega =0\end{array}\right.\end{split}\]

Où \(\kappa\) est le module de compressibilité.

2.2. Formulation en grandes déformations

Comme pour les petites déformations, il est possible de proposer une formulation variationnelle valable pour les grandes déformations. Le principe est identique, mais on s’appuie dans ce cas sur la décomposition du tenseur gradient de la transformation \(F=I+\frac{\partial u}{\partial X}\) proposée par Flory [bib3] :

\[F={F}^{s}\stackrel{ˉ}{F} \ \text{avec} \ {F}^{s}={J}^{1/3}I \ \ \text{et} \ \ \stackrel{ˉ}{F}={J}^{-1/3}F \ \ \text{et} \ \ J=\det(F)\]

L’idée là encore, est d’enrichir la cinématique par le biais d’une variable de gonflement \(g\) , a priori indépendante des déplacements, et liée faiblement à la variation de volume par une relation faible:

\[B(J)\underset{\mathit{faible}}{\approx}B\mathrm{°}A(g)=B(A(g))\]

Plusieurs relations ont été testées:

\[\begin{split}\left\lbrace \begin{array}{c}J=1+g\\ {J}^{2}=1+g\\ \ln(J)=g\\ J=\exp(g)\end{array} \right.\end{split}\]

Pour certaines simulations, de petites différences ont été observées. Pour les éléments INCO_UPG avec déformation SIMO_MIEHE, c’est finalement la relation linéaire qui a été implantée dans le code: donc \(B(J)=J\) et \(A(g)=1+g\) . Pour les éléments INCO_UPG avec déformation GDEF_LOG, c’est la relation en logarithme qui a été retenue.

Néanmoins, ce choix n’étant pas forcément définitif, on propose d’écrire le problème dans le cas général. On introduit donc un gradient de déformation enrichi:

(2.5)\[\tilde{F}={(\frac{A(g)}{J})}^{\frac{1}{3}}F\]

La formulation faible du problème s’appuie sur la recherche du point selle du Lagrangien \(ℒ\) , dans lequel le multiplicateur de Lagrange \(p\) et un troisième champ \(g\) , indépendant des deux autres, assurant de façon faible que la relation entre \(J\) et \(g\) est vérifiée:

(2.6)\[\mathrm{ℒ}(u,g,p)=\underset{{\Omega}_{0}}{\int}\psi (\tilde{F})d{\Omega}_{0}-{W}_{\mathit{ext}}(u)+\underset{{\Omega}_{0}}{\int}p\left[B(J)-B\mathrm{°}A(g)\right]d{\Omega}_{0}\]

où \({W}_{\mathrm{ext}}\) représente le potentiel des efforts extérieurs et \(\psi (\tilde{F})\) l’énergie de déformation.

Ce problème peut être résolu comme en petites déformations en écrivant les conditions d’optimalité. La variation du Lagrangien s’écrit:

(2.7)\[\delta ℒ=\underset{{\Omega}_{0}}{\int}\left[P:\delta \tilde{F}+p(\frac{\partial B(J)}{\partial J}J\frac{\delta J}{J}-\frac{\partial B°A(g)}{\partial g}\delta g)+\delta p(B(J)-B°A(g))\right]d{\Omega}_{0}-\delta {W}_{\mathrm{ext}}(u)\]

avec \(P\) le premier tenseur des contraintes de Piola-Kirchhoff:

En injectant la variation de la transformation enrichie et l’expression de la contrainte de Kirchhoff \(\tau =P\tilde{{F}^{T}}\) [1], on obtient la forme suivante pour la variation du Lagrangien:

(2.8)\[\begin{split}\begin{array}{ccc}\delta \mathrm{ℒ}& =& \underset{{\Omega}_{0}}{\int}({\tau}^{d}+p\frac{\partial B(J)}{\partial J}JI):\delta Ld{\Omega}_{0}\\ & +& \underset{{\Omega}_{0}}{\int}(\frac{\text{tr}(\tau )}{3}\frac{\partial A(g)/\partial g}{A(g)}-p\frac{\partial B\mathrm{°}A(g)}{\partial g})\delta gd{\Omega}_{0}+\underset{{\Omega}_{0}}{\int}(B(J)-B\mathrm{°}A(g))\delta pd{\Omega}_{0}\\ & -& \delta {W}_{\mathit{ext}}(u)\end{array}\end{split}\]

où l’on a introduit le gradient eulérien du déplacement (\(x\) représente le vecteur position à la fin de l’incrément): \(\delta L=\frac{\partial \delta u}{\partial x}=\delta F.{F}^{-1}\)

En résumé, le système à résoudre est le suivant :

(2.9)\[\begin{split}\begin{array}{c}\frac{\partial \mathrm{ℒ}}{\partial u}=\underset{{\Omega}_{0}}{\int}\delta L:({\tau}^{d}+p\frac{\partial B(J)}{\partial J}JI)d{\Omega}_{0}-\delta {W}_{\mathit{ext}}=0\\ \frac{\partial \mathrm{ℒ}}{\partial g}=\underset{{\Omega}_{0}}{\int}\delta g(\frac{\text{tr}(\tau )}{3}\frac{\partial A(g)/\partial g}{A(g)}-p\frac{\partial B\mathrm{°}A(g)}{\partial g})d{\Omega}_{0}=0\\ \frac{\partial \mathrm{ℒ}}{\partial p}=\underset{{\Omega}_{0}}{\int}\delta p(B(J)-B\mathrm{°}A(g))d{\Omega}_{0}=0\end{array}\end{split}\]

Remarque : La contrainte de Kirchhoff issue de la loi de comportement, s’écrit donc \(\tau ={\tau}^{D}+p\frac{\partial B(J)}{\partial J}JI\)

En ce qui concerne l’obtention de la matrice tangente, elle demande bien sûr un peu plus de calcul qu’en petites déformations, et a la particularité de ne pas être symétrique dans le cas général. Elle se trouve dans le code sous la forme suivante :

\[\begin{split}K=\left[\begin{array}{ccc}{K}_{\text{UU}}& {K}_{\text{UG}}& {K}_{\text{UP}}\\ {K}_{\text{GU}}& {K}_{\text{GG}}& {K}_{\text{GP}}\\ {K}_{\text{PU}}& {K}_{\text{PG}}& {K}_{\text{PP}}\end{array}\right]\end{split}\]

Les calculs ne sont pas détaillés ici. Le lecteur pourra se reporter à la lecture de [bib8].

Remarque: Cette formulation permet de régulariser à moindre coût, les modèles d’endommagement ductile où la variable d’endommagement est directement reliée à la variation de volume. En effet, pour contrôler la localisation de l’endommagement et de la déformation, l’idée est de pénaliser les forts gradients d’endommagement. Comme dans cette formulation à 3 champs, le gonflement local est traité comme une variable nodale, son gradient est facilement accessible numériquement (sous réserve d’une interpolation au moins linéaire).

Dans l’esprit des formulations à second gradient de déplacement ([bib9], [bib10]), elle est enrichie par un terme quadratique en gradient de gonflement. La variation du lagrangien s’écrit alors:

(2.10)\[\begin{split}\begin{array}{ccc}\delta ℒ& =& \underset{{\Omega}_{0}}{\int}\left({\tau}^{d}+p\frac{\partial B(J)}{\partial J}JI\right):\delta L+\left(\frac{\text{tr}(\tau )}{3}\frac{\partial A(g)/\partial g}{A(g)}-p\frac{\partial B\mathrm{°}A(g)}{\partial g}\right)\delta gd{\Omega}_{0}\\ & +& \underset{{\Omega}_{0}}{\int}\left[\left(B(J)-B\mathrm{°} A(g)\right)\delta p+c\nabla g\cdot \nabla \delta g\right]d{\Omega}_{0}-\delta {W}_{\mathit{ext}}(u)\end{array}\end{split}\]

\(c\) est un paramètre à déterminer et homogène à une force. Ce paramètre introduit en quelque sorte une longueur interne de couplage entre les points matériaux. Le terme ajouté ici est isotrope: on considère que la longueur interne à introduire est identique dans toutes les directions. Pour l’application à l’endommagement ductile des aciers, cette hypothèse semble tout à fait recevable. Cette formulation est utilisable pour le modèle de Rousselier, [R5.03.07], moyennant la définition du mot-clé C_CARA sous l’opérande NON_LOCAL de DEFI_MATERIAU (voir test ssnp122a)

L’obtention de la formulation à deux champs en grandes déformations suit le même principe que En ce qui concerne l’obtention de la matrice tangente, elle demande bien sûr un peu plus de calcul qu’en petites déformations, et a la particularité de ne pas être symétrique dans le cas général. Elle se trouve dans le code sous la forme suivante :

[1]

Il s’agit là d’un choix et l’on aurait tout aussi bien pu transporter \(P\) par la déformation compatible \(F\) au lieu de la déformation enrichie \(\tilde{F}\). Mais cette dernière a l’avantage de préserver la symétrie de la formulation dans le cas élastique. En outre, elle se prête mieux à l’architecture d’un code de calcul dans lequel l’intégration des lois de comportement est bien différenciée du calcul des termes propres aux éléments finis (les routines du comportement n’ont pas besoin de connaître l’existence des deux mesures de déformation) mais seulement la déformation enrichie).