r3.06.08 Éléments finis traitant la quasi-incompressibilité

Résumé:

Dans certaines situations, le comportement mécanique du matériau impose que la dilatation volumique reste nulle, autrement dit que la déformation se fasse à volume constant : élasticité isotrope avec coefficient de Poisson égal à \(0.5\), écoulements plastiques parfaits en analyse limite …

On propose ici de traiter cette condition d” «incompressibilité» ou de «quasi-incompressibilité» en utilisant une formulation valable aussi bien dans le cas compressible que dans le cas quasi-incompressible. Pour cela, on utilise une formulation variationnelle à 3 champs où les inconnues sont le déplacement, la déformation volumique et le multiplicateur de Lagrange associé (qui correspondrait à la pression dans le cas incompressible). On propose deux versions de cette formulation: l’une pour les petites déformations, l’autre valable en présence de grandes déformations. Dans la situation d’une relation biunivoque entre la pression et le gonflement, cas de la plasticité de Von Mises, il est possible de venir éliminer l’inconnue de gonflement. On a alors une formulation à deux champs déplacement/pression.

Après quelques rappels sur les difficultés que posent la résolution des problèmes incompressibles, on décrit les éléments finis mixtes implantés (en 3D et en 2D, plan et axisymétrique en petites et en grandes déformations), et on présente également les grandes lignes de l’intégration dans code_aster (modélisations INCO_UP, INCO_UPG, INCO_UPO).

Cette modélisation est nécessaire pour pratiquer les analyses limites et pour modéliser des comportements élastiques pour des coefficients de Poisson proche de \(0.5\) . Elle peut aussi être utile dans le cas de modélisations engendrant de fortes déformations plastiques et pour lesquelles les modélisations traditionnelles peuvent être insuffisantes et engendrer des oscillations de contraintes.

Formulation variationnelle mixte du problème

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 [r3.06.08-difficulte-inco], 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:

(357)\[\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:

(358)\[ℒ(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é :

(359)\[\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 (359).

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 (359). On obtient alors le système de deux équations à deux inconnues qui suit:

(360)\[\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é.

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:

(361)\[\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:

(362)\[\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:

(363)\[\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:

(364)\[\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 :

(365)\[\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:

(366)\[\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).

Discrétisation par éléments finis mixtes

Choix de la discrétisation

Lorsqu’on utilise une formulation mixte, il est nécessaire de discrétiser à la fois l’espace des déplacements, du multiplicateur de Lagrange p et du «gonflement» \(g\) . L’expérience acquise sur les éléments mixtes, notamment à 2 champs pour les éléments incompressibles, permet de savoir que la discrétisation de ces champs ne peut être quelconque, sous peine d’obtenir des phénomènes d’oscillations (notamment au niveau des pressions) ou des phénomènes de blocage (éléments ne pouvant pas se déformer ou trop rigides). Ainsi il est nécessaire d’avoir un nombre de points de Gauss de pression suffisamment important pour vérifier la condition d’incompressibilité presque partout et un nombre de points de Gauss de pression suffisamment faible pour avoir plus de degrés de liberté à calculer que de contraintes à vérifier. Une des conditions nécessaires pour obtenir des résultats satisfaisants est la vérification par l’élément fini considéré de la condition LBB (LADYJENSKAIA, BREZZI, BABUSKA). On peut trouver dans [bib5] et [bib6] des exemples d’éléments satisfaisant la condition LBB.

Ici le problème est un peu différent lorsque la formulation 3 champs est retenue.

En l’état actuel, les discrétisations utilisées ne sont pas les mêmes dans la version HPP et dans la version grandes déformations.

Petites déformations

Pour les petites déformations, nous nous sommes inspirés des usages classiques des formulations mixtes (ex. [bib7]), en utilisant un élément de type \(\mathit{P2}/\mathit{P1}/\mathit{P1}\) pour la formulation à 3 champs. Autrement dit, le déplacement est quadratique, la pression et le gonflement sont tous les deux linéaires. Les éléments finis utilisés pour la formulation à 3 champs sont donc les suivants:

en 2D:	\(u\)	triangle à 6 nœuds	quadrilatère à 8 nœuds
	\(p,g\)	triangle à 3 nœuds	quadrilatère à 4 nœuds
en 3D:	\(u\)	tétraèdre à 10 nœuds	héxaèdre à 20 nœuds	pentaèdre à 15 nœuds
	\(p,g\)	tétraèdre à 4 nœuds	héxaèdre à 8 nœuds	pentaèdre à 6 nœuds

Pour chaque type d’élément, on utilise une seule famille de points de Gauss:

• 3 points pour les triangles

• 4 points pour les quadrilatères

• 4 points pour les tétraèdres

• 8 points pour les héxaèdres

• 21 points pour les pentaèdres

Pour la formulation à deux champs, un élément du type \(\mathit{P2}/\mathit{P1}\) a été introduit. Le déplacement a donc une interpolation quadratique tandis que la pression est interpolée linéairement. Dans le cas de l’utilisation d’une discrétisation en triangles ou tétraèdres linéaires, deux méthodes de stabilisation ont été introduite. La première correspond à l’élément fini stabilisé \(\mathit{P1}\text{+/}\mathit{P1}\) . Le \(+\) correspond à l’introduction d’un degré de liberté supplémentaire au centre de l’élément dans l’interpolation des déplacements. Ce degré supplémentaire est appelé communément «bulle». Cette méthode de stabilisation ne fonctionne que sur des éléments simplexes (triangle en 2D et tétraèdre en 3D). Elle a l’avantage d’utiliser très peu de degrés de liberté. La deuxième méthode de stabilisation correspond à la méthode Orthogonal Sub-Grid Scale (OSGS) [bib11]. L’avantage de cette méthode est de fonctionner pour toutes les topologies d’éléments. Son principal inconvénient est d’introduire une troisième inconnue (et donc des degrés de liberté supplémentaires) correspondant au champ de pression projeté \(\pi\) sur l’espace orthogonal aux champs de déplacement.

Les éléments finis utilisés pour la formulation à 2 champs sont donc les suivants:

Interpolation		P1+/P1	P1/P1 OSGS	P2/P1	P1/P1 OSGS	P2/P1	P1/P1 OSGS	P2/P1
en 2D:	\(u\)	triangle à 3 nœuds	triangle à 3 nœuds	triangle à 6 nœuds	quadrilatère à 4 nœuds	quadrilatère à 8 nœuds
	\(p\)	triangle à 3 nœuds	triangle à 3 nœuds	triangle à 3 nœuds	quadrilatère à 4 nœuds	quadrilatère à 4 nœuds
	\(\pi\)		triangle à 3 nœuds		quadrilatère à 4 nœuds
en 3D:	\(u\)	tétraèdre à 4 nœuds + bulle	tétraèdre à 4 nœuds	tétraèdre à 10 nœuds	cube à 8 nœuds	cube à 20 nœuds	pentaèdre à 6 nœuds	pentaèdre à 15 nœuds
	\(p\)	tétraèdre à 4 nœuds	tétraèdre à 4 nœuds	tétraèdre à 4 nœuds	cube à 8 nœuds	cube à 8 nœuds	pentaèdre à 6 nœuds	pentaèdre à 6 nœuds
	\(\pi\)		tétraèdre à 4 nœuds		cube à 8 nœuds		pentaèdre à 6 nœuds

Les familles de points de Gauss utilisées sont les mêmes que celles de la formulation à 3 champs. On notera que pour les éléments \(\mathit{P1}\text{+/}\mathit{P1}\) , on n’utilise qu’un seul point de Gauss pour l’intégration.

Grandes déformations

A partir de la version 11, le choix des interpolations en grandes déformations est identique à celui des petites déformations. Les éléments sont du type \(\mathit{P2}/\mathit{P1}/\mathit{P1}\) pour les formulations à 3 champs et \(\mathit{P2}/\mathit{P1}\) pour la formulation à 2 champs. Les éléments finis utilisés pour la formulation à 3 champs sont donc les suivants:

en 2D:	\(u\)	triangle à 6 nœuds	quadrilatère à 8 nœuds
	\(p,g\)	triangle à 3 nœuds	quadrilatère à 4 nœuds
en 3D:	\(u\)	tétraèdre à 10 nœuds	héxaèdre à 20 nœuds	pentaèdre à 15 nœuds
	\(p,g\)	tétraèdre à 4 nœuds	héxaèdre à 8 nœuds	pentaèdre à 6 nœuds

Ce sont les mêmes familles de points de Gauss que celles des petites déformations qui ont été utilisées.

Écriture du problème discret

On aborde ici dans un premier temps, l’écriture du problème discret dans le cadre de la formulation à 3 champs. Soit \({u}^{e}\) , \({p}^{e}\) et \({g}^{e}\) , les vecteurs des inconnues nodales élémentaires (respectivement déplacement, pression et gonflement). Si \({N}^{u}\) , \({N}^{p}\) et \({N}^{g}\) sont les fonctions de formes (respectivement des interpolations des déplacement, de la pression et du gonflement) associées à l’élément fini considéré:

\[\begin{split}\begin{array}{c} u={N}^{u}{u}^{e} \\ p={N}^{p}{p}^{e} \\ g={N}^{g}{g}^{e} \end{array}\end{split}\]

Écriture en petites déformations

\(B\) est la matrice de dérivation classique permettant de passer de \({u}^{e}\) à \(\varepsilon\) :

\[\varepsilon =B{u}^{e}\]

Dans la formulation, on distingue \({e}_{\text{dev}}\) et \({e}_{\text{dil}}\) , ce qui nous amène à définir les opérateurs \({B}_{\text{dev}}\) et \({B}_{\text{dil}}\) tels que: \({\epsilon}^{D}={B}_{\text{dev}}{U}^{e}\) et \(\frac{\text{tr}\epsilon }{3}={B}_{\text{dil}}{U}^{e}\)

La forme discrétisée des équations du problème à 3 champs [éq ] s’écrit:

\[\begin{split}\begin{array}{l} {F}_{u}=\underset{\Omega}{\int}{B}^{T}({\sigma}^{D}+pI)d\Omega ={F}_{\mathrm{ext}}\\ {F}_{p}=\underset{\Omega}{\int}{({N}^{p})}^{T}({B}_{\mathrm{dil}}u-g)d\Omega =0\\ {F}_{g}=\underset{\Omega}{\int}{({N}^{g})}^{T}(\frac{1}{3}\text{tr}(\sigma )-p)d\Omega =0\end{array}\end{split}\]

La matrice tangente du problème est symétrique et s’appuie sur les termes suivants:

\[\begin{split}\begin{array}{l} {K}_{\mathrm{uu}}=\frac{\partial {F}_{u}}{\partial {u}^{e}}=\underset{\Omega}{\int}{B}_{\mathrm{dev}}^{T}D{B}_{\mathrm{dev}}d\Omega \\ {K}_{\mathrm{up}}=\frac{\partial {F}_{u}}{\partial {p}^{e}}=\underset{\Omega}{\int}{B}_{\mathrm{dil}}^{T}{N}^{p}d\Omega \\ {K}_{\mathrm{ug}}=\frac{\partial {F}_{u}}{\partial {g}^{e}}=\frac{1}{3}\underset{\Omega}{\int}\text{tr}({B}_{\mathrm{dev}}^{T}D){N}^{g}d\Omega \\ {K}_{\mathrm{pp}}=\frac{\partial {F}_{p}}{\partial {p}^{e}}=0 \\ {K}_{\mathrm{pg}}=\frac{\partial {F}_{p}}{\partial {g}^{e}}=-\underset{\Omega}{\int}{({N}^{p})}^{T}{N}^{g}d\Omega \\ {K}_{gg}=\frac{\partial {F}_{g}}{\partial {g}^{e}}=\frac{1}{9}\underset{\Omega}{\int}{({N}^{g})}^{T}\text{tr}(D){N}^{g}d\Omega\end{array}\end{split}\]

En ce qui concerne la formulation à 2 champs, elle se déduit aisément de la précédente. Les formes discrétisée des équations nous donnent:

\[\begin{split}\begin{array}{l} {F}_{u}=\underset{\Omega}{\int}{B}^{T}({\sigma}^{D}+pI)d\Omega ={F}_{\mathit{ext}}\\ {F}_{p}=\underset{\Omega}{\int}{({N}^{p})}^{T}({B}_{\mathrm{dil}}u-\frac{{N}^{p}}{\kappa})d\Omega =0\end{array}\end{split}\]

La matrice tangente du problème est symétrique et s’appuie sur les termes suivants:

\[\begin{split}\begin{array}{l} {K}_{\mathrm{uu}}=\frac{\partial {F}_{u}}{\partial {u}^{e}}=\underset{\Omega}{\int}{B}_{\mathrm{dev}}^{T}D{B}_{\mathrm{dev}}d\Omega\\ {K}_{\mathrm{up}}=\frac{\partial {F}_{u}}{\partial {p}^{e}}=\underset{\Omega}{\int}{B}_{\mathrm{dil}}^{T}{N}^{p}d\Omega\\ {K}_{\mathit{pp}}=\frac{\partial {F}_{p}}{\partial {p}^{e}}=-\frac{1}{\kappa}\underset{\Omega}{\int}{({N}^{p})}^{T}{N}^{p}d\Omega\end{array}\end{split}\]

Écriture en grandes transformations

L’écriture étant un peu fastidieuse, le lecteur pourra se reporter à la lecture de [bib8] pour avoir plus d’informations.

Intégration dans Code_Aster des éléments finis incompressibles

Présentation générale de l’élément incompressible en petites déformations

Les éléments finis sont intégrés dans Code_Aster en 2D déformations planes, en 2D axisymétrique et en 3D. Les 3 modélisations sont accessibles en utilisant les options suivantes pour AFFE_MODELE:

• 3D_INCO_UPG, 3D_INCO_UP ou 3D_INCO_UPO pour le 3D et respectivement pour la formulation à 3, à 2 champs et à 2 champs stabilisés avec la méthode OSGS,

• D_PLAN_INCO_UPG, D_PLAN_INCO_UP ou D_PLAN_INCO_UPO pour le 2D en déformations planes et respectivement pour la formulation à 3, à 2 champs et à 2 champs stabilisés avec la méthode OSGS,

• AXIS_INCO_UPG, AXIS_INCO_UP ou AXIS_INCO_UPO pour le 2D axisymétrique et respectivement pour la formulation à 3, à 2 champs et à 2 champs stabilisés avec la méthode OSGS.

Dans le catalogue des éléments, les éléments incompressibles peuvent s’appliquer sur les mailles:

Mailles	Formulation	Nombre de nœuds en déplacements	Nombre de nœuds en pression ou gonflement	Nombre de nœuds de gradient de pression projeté
TRIA3	2 champs	3	3
TRIA3	2 champs OSGS	3	3	3
TRIA6	2 et 3 champs	6	3
QUAD4	2 champs OSGS	4	4	4
QUAD8	2 et 3 champs	8	4
HEXA20	2 et 3 champs	20	8
TETRA4	2 champs	4	4
TETRA4	2 champs OSGS	4	4	4
TETRA10	2 et 3 champs	10	4
PENTA6	2 champs OSGS	6	6	6
PENTA15	2 et 3 champs	15	6

Dans les routines d’initialisations des éléments incompressibles, on définit:

• 1 seule famille de points de Gauss (cf § Choix de la discrétisation ),

• 2 familles de fonctions de formes associées respectivement aux déplacements (fonctions de formes de degré 2) et aux termes de pression et de gonflement (de degré 1) si on est en formulation 3 champs .

Prenons comme exemple l’élément tétraédrique à 10 nœuds: les degrés de liberté en déplacement sont portés par tous les nœuds, en revanche, seuls les 4 nœuds sommets possèdent les degrés de liberté \(p\) et \(g\) .

Les composantes accessibles pour le champ DEPL sont donc:

• les déplacements: DX, DY et DZ en 3D à tous les nœuds,

• la pression: PRES pour les nœuds sommets,

• le gonflement (formulation à 3 champs): GONF pour les nœuds sommets.

Présentation générale de l’élément incompressible en grandes déformations

Les éléments finis sont intégrés dans code_aster en 2D déformations planes, en 2D axisymétrique et en 3D. Les 3 modélisations s’appuyant sur une formulation à 3 champs sont accessibles en utilisant les options suivantes pour AFFE_MODELE:

• 3D_INCO_UPG pour le 3D,

• D_PLAN_INCO_UPG pour le 2D en déformations planes,

• AXIS_INCO_UPG pour le 2D axisymétrique.

Les 3 modélisations utilisant le formalisme de grandes déformations de GDEF_LOG et s’appuyant sur une formulation à 2 champs sont accessibles en utilisant les options suivantes pour AFFE_MODELE:

• 3D_INCO_UP pour le 3D,

• D_PLAN_INCO_UP pour le 2D en déformations planes,

• AXIS_INCO_UP pour le 2D axisymétrique.

Dans le catalogue des éléments, les éléments incompressibles peuvent s’appliquer sur les mailles:

Mailles	Nombre de nœuds en déplacements	Nombre de nœuds en pression (et gonflement)
TRIA6	6	3
QUAD8	8	4
HEXA20	20	8
TETRA10	10	4
PENTA15	15	6

Dans les routines d’initialisation des éléments incompressibles, on définit:

• 1 seule famille de points de Gauss (cf § Choix de la discrétisation ),

• 2 familles de fonctions de formes associées respectivement aux déplacements et à la pression (fonctions de formes de degré 2) et aux termes de gonflement (fonctions de formes de degré 1).

Prenons comme exemple l’élément tétraédrique à 10 nœuds: les degrés de liberté en déplacement et en pression sont portés par tous les nœuds, en revanche, seuls les 4 nœuds sommets possèdent les degrés de liberté de gonflement.

Les composantes accessibles pour le champ DEPL sont donc

• les déplacements: DX, DY et DZ en 3D à tous les nœuds,

• la pression: PRES pour les nœuds sommet,

• le gonflement: GONF pour les nœuds sommet.

Emploi de la modélisation

Les modélisations INCO_UP, INCO_UPG et INCO_UPO peuvent être utilisées avec les opérateurs de mécaniques non-linéaires STAT_NON_LINE et DYNA_NON_LINE. Il est aussi possible d’utiliser l’opérateur de mécanique linéaire MECA_STATIQUE cependant ceci est fortement déconseillé car les résultats obtenus peuvent être de qualité médiocre. La version petites déformations est accessible en utilisant DEFORMATION=’PETIT’ sous COMPORTEMENT, la version grandes déformations en utilisant DEFORMATION=’SIMO_MIEHE’ ou DEFORMATION=’GDEF_LOG’. Les relations de comportement utilisables sont celles disponibles respectivement en petites déformations et en grandes déformations SIMO_MIEHE ou GDEF_LOG pour les modélisation INCO_UPG. Les modélisations INCO_UP et INCO_UPO sont actuellement limitées aux relations ELAS et VMIS_ISOT_XXX.

Il n’est donc pas possible d’utiliser les modélisations avec les commandes: CALC_MATR_ELEM/CALC_VECT_ELEM/ASSE_MATRICE/ASSE_VECTEUR/RESOUDRE

Étant donnée la forme de la matrice tangente pour les formulations à 3 champs (INCO_UPG), il est souvent nécessaire d’utiliser le solveur MUMPS pour résoudre les systèmes linéaires.

Il est conseillé d’utiliser le critère de convergence par contrainte de référence RESI_REFE_RELA sauf pour les éléments INCO_UPO car l’option n’est pas disponible.

Formulation des termes élémentaires du second membre

Les charges peuvent être la pesanteur, des forces surfaciques réparties, des pressions. Les termes élémentaires sont calculés de façon classique pour les degrés de liberté de déplacement et on affecte la valeur nulle pour les degrés de liberté de pression et de gonflement.

Calcul des déformations et des contraintes

Dans cette formulation, il convient de distinguer le champ de contrainte issu de la loi de comportement \({\sigma}_{\mathit{ldc}}\), du champ de contrainte qui vérifie l’équilibre et qui est défini par la relation \(\sigma ={\sigma}_{\mathit{ldc}}^{D}+pI\) . C’est ce dernier champ qui est stocké dans SIEF_ELGA ainsi que la relation liant le multiplicateur \(p\) et \({\sigma}_{\mathit{ldc}}\).

En petites déformations, les composantes de SIEF_ELGA sont:

• SIXX, SIYY, SIZZ, SIXY en 2D ainsi que SIXZ et SIYZ en 3D: composantes du tenseur \(\sigma ={\sigma}_{\mathit{ldc}}^{D}+pI\) ,

• SIP qui est égal à \(\left(\frac{1}{3}\text{tr}({\sigma}_{\mathit{ldc}})-p\right)\) ,

En grandes déformations, les composantes de SIEF_ELGA sont:

• SIXX, SIYY, SIZZ, SIXY en 2D ainsi que SIXZ et SIYZ en 3D: composantes du tenseur \(\sigma =\left(\frac{{\tau}^{d}}{J}+p\frac{\partial B(J)}{\partial J}I\right)\) ,

• SIP qui est égal à \(\left(\frac{\text{tr}(\tau )}{\mathrm{3J}}-p\frac{\partial B(J)}{\partial J}\right)\) ,

Il est également possible de recalculer EPSI_ELGA, qui est le champ de déformation au sens classique.

On peut également réaliser un calcul de charge limite avec POST_ELEM.

Validation

Cas élastique incompressible

Le test SSLV130 (cf. [V3.04.130]) permet de vérifier la validité de la modélisation dans le cas d’un cylindre élastique incompressible soumis à une pression interne. Son équivalent en grandes déformations existe également: test SSNV112 (cf. [V6.04.112]).

Cas élasto-plastique

Le but de cet exemple est d’illustrer l’apport de la modélisation incompressible dans le cas où les déformations plastiques sont importantes par rapport aux déformations élastiques. On étudie pour cela une éprouvette entaillée en axisymétrique, soumis à un déplacement imposé. La géométrie et le chargement sont représentés sur la figure ci-dessous. Le maillage est constitué de 548 TRI6.

Fig. 134 Géométrie et conditions aux limites

Le comportement du matériau est de type élastoplastique à écrouissage isotrope linéaire(VMIS_ISOT_LINE). Les paramètres sont les suivants:

• \(E=200000\mathit{MPa}\)

• \(\nu =0.3\)

• \({\sigma}_{y}=200\mathit{MPa}\)

• \({E}_{T}=1000\mathit{MPa}\)

Sur la figure suivante, on compare la contrainte \({\sigma}_{yy}\) obtenue sur le chemin \(\mathit{FC}\) (cf. [Géométrie et conditions aux limites]) avec la modélisation classique AXIS et la modélisation AXIS_INCO_UPG.

Fig. 135 \({\sigma}_{yy}\) le long de la ligne \(\mathit{FC}\)

On voit très clairement que la solution obtenue avec la formulation incompressible permet de s’affranchir des oscillations parasites.

Bibliographie

[bib1]

J.R. Hugues: The finite element method, Prentice-Hall, Inc. Englewood Cliffs, N-J. 07632, 1987.

[bib2]

J.C. Simo, M.S. Rifai : A class of mixed assumed strain methods and the method of incompatible modes, Int. Jnal Num. Meth. Engg, Vol. 29, pp1595-1638, 1990.

[bib3]

R.J. Flory : Thermodynamic relations of high elastic materials. Trans. Faraday Soc., vol.57, 1961, pp. 829-838.

[bib4]

V. Cano, E. Lorentz: Introduction dans le Code_Aster d’un module de comportement en grandes déformations élastoplastique avec écrouissage isotrope. Note EDF/DER HI‑74/98/006/0 du 26/08/1998

[bib5]

M. Girault, P. Raviart: Finite element methods for Navier-Stokes equations. Theory and algorithms. Springer Verlag, 1986.

[bib6]

P. Mialon, B. Thomas: Incompressibilité en plasticité: sous-intégration et autres techniques numériques. Note EDF/DER HI-72/6404 du 19/01/1990.

[bib7]

A.G.K. Jinka, M. Bellet, L. Fourment : A new three-dimensional finite element model for the simulation of powder forging processes : application to hot forming of P/M connecting road. Int. Jnal Num. Meth. Engg, Vol. 40, pp3955-3978, 1997.

[bib8]

S. Michel-Ponnelle, R. Bargellini: Validation d’éléments finis à 3 champs pour traiter les grandes déformations plastiques et régulariser le modèle de Rousselier, Note EDF R&D H-B30-2009-00094, 2009

[bib9]

P. Germain, 1973, La méthode des puissances virtuelles en mécanique des milieux continus; Première partie : Théorie du second gradient, J. Méca. 12, 235-274

[bib10]

R. Chambon, D. Caillerie, T. Matsuchima,2001, Plastic continuum with microstructure, local second gradient theories for geomaterials: Localization studies, Int. J. Solids Structures 38, 8503-8527

[bib11]

C. Agelet De Saracibar, M. Chiumenti, Q. Valverde, M. Cervera, On the orthogonal subgrid scale pressure stabilization of finite deformation J2 plasticity, Computer Methods in Applied Mechanics and Engineering, Vol. 195, n°9-12, pp. 1224-1251, 2006.

Description des versions

Version Aster	Auteur(s) Organisme(s)	Description des modifications
6.2	S.MICHEL-PONNELLE, E.LORENTZ EDF R&D AMA	Version initiale
6.4	S.MICHEL-PONNELLE, E.LORENTZ EDF R&D AMA	Mise à jour légère pour la version 6.4
7.2	S.MICHEL-PONNELLE, E.LORENTZ EDF R&D AMA	Ajout de la formulation en grandes transformations
7.4	S.MICHEL-PONNELLE, E.LORENTZ EDF R&D AMA	Ajout des éléments pentaèdres
9.4	S.MICHEL-PONNELLE EDF R&D AMA E.LORENTZ EDF R&D SINETICS	Nouvelle formulation éléments quasi-incompressibles en grandes transformations
10.3	S.FAYOLLE EDF R&D AMA	Nouvelle formulation à 2 champs en petite déformation \(U-P\)
11.2	S.FAYOLLE EDF R&D AMA	Nouvelles formulations OSGS et INCO_LOG
11.3	S.FAYOLLE EDF R&D AMA	Nouvelles formulations INCO_LUP
12.1	S.FAYOLLE EDF R&D AMA	Renommage des formulations (fiches 16002, 21945 et 21921)