# Méthodes éléments finis robustes pour gérer l'incompressibilité en grandes déformations

Dina Al Akhrass, Sébastien Fayolle, Sylvain Drapier, Julien Bruchon

► To cite this version:

Dina Al Akhrass, Sébastien Fayolle, Sylvain Drapier, Julien Bruchon. Méthodes éléments finis robustes pour gérer l'incompressibilité en grandes déformations. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722340

HAL Id: hal-01722340

https://hal.science/hal-01722340v1

Submitted on 3 Mar 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Méthodes éléments finis robustes pour gérer l'incompressibilité en grandes déformations

Dina AL AKHRASS <fnref n="1" />,<fnref n="2" />, Sébastien FAYOLLE <fnref n="2" />, Sylvain DRAPIER <fnref n="1" />, Julien BRUCHON <fnref n="1" />

**Résumé** — Le présent travail traite de la gestion de l'incompressibilité en grandes déformations, dans le cadre de la mécanique des solides. Nous proposons pour cela de travailler avec un modèle de grandes déformations basé sur des déformations logarithmiques, adapté pour une formulation mixte à trois et deux champs. Les performances des différents éléments obtenus sont ensuite comparées, puis validées sur un cas industriel.

**Mots clés** — Incompressibilité, grandes déformations, éléments finis mixtes, stabilisation

## 1 Introduction

Dans le cadre de la mécanique des solides, la simulation d'applications faisant intervenir des non-linéarités comme les grandes déformations ou des lois de comportement non-linéaires présente des difficultés numériques. En effet, lorsque l'on fait des calculs impliquant une loi de comportement élastoplastique, et que l'on atteint un niveau de déformation plastique important, cela peut engendrer le phénomène d'incompressibilité. Il est connu que les éléments finis classiques, utilisant comme seule inconnue les déplacements, ne sont alors plus adéquats pour rendre compte correctement de ce phénomène, donnant lieu par exemple à des effets de blocage. Depuis les années 70, différentes stratégies ont été proposées pour réduire ou éliminer le blocage volumique et les oscillations de pression dans les solutions éléments finis. On peut citer entre autres les éléments sous-intégrés, les méthodes *Enhanced Assumed Strain* (EAS) [6], B-bar et Fbar [9], ou les formulations mixtes. La méthode des éléments finis mixtes [1] est une méthode robuste, générique et qui donne des résultats précis. Elle fournit un cadre théorique qui permet de développer une formulation mixte à trois champs d'inconnues, en déplacement, gonflement et pression par exemple [11], dans laquelle une relation renforçant la contrainte d'incompressibilité et liant gonflement et pression est imposée via un multiplicateur de Lagrange. Lorsqu'il existe une relation bi-univoque entre la pression et le gonflement, comme dans le cas de la plasticité J2, ce dernier peut être eliminé pour obtenir une formulation à deux champs (déplacement et pression) [10]. La méthode des éléments finis mixtes est un moyen efficace de gérer l'incompressibilité, à condition que les interpolations des différentes inconnues vérifient la condition de stabilité de Ladyjenskaia-Brezzi-Babuska (LBB). Ainsi, pour la formulation à deux champs, l'élément P2/P1 (quadratique en déplacement et linéaire en pression) vérifie bien cela pour le problème mécanique. Il est cependant possible de travailler avec des éléments qui ne vérifient pas la condition, comme l'élément P1/P1 (linéaire en déplacement et en pression) mais il est alors nécessaire de le stabiliser en modifiant la formulation variationnelle du problème mécanique. De nombreuses méthodes de stabilisation ont été développées ces dernières décennies. L'une des plus connues conduit au mini-élément (P1+/P1), avec la prise en compte d'un terme "bulle" en déplacement dans chaque élément, qui, une fois condensé, peut être considéré comme un terme de stabilisation [4] ; plus récemment, les méthodes VMS (Variational Multi-Scale) comme la méthode OSGS (Orthogonal Sub-Grid Scale) [3] ont étendu et formalisé cette démarche, en considérant les inconnues du problème comme la somme d'une partie résoluble (la solution éléments finis) et d'une partie que l'on doit approcher qui est liée à ce que le maillage ne permet pas de capturer. Les méthodes évoquées sont génériques.

Pour ce qui est du formalisme de grandes déformations, les premiers ont été introduits dans les années 70, et sont basés sur des concepts hypoélastiques [8]. Une des faiblesses de ce type de modèle est

<footnote n="1">¹ Ecole des Mines de Saint-Etienne, Centre SMS et LGF UMR CNRS 5307, alakhrass@emse.fr, drapier@emse.fr, bruchon@emse.fr</footnote>

<footnote n="2">² EDR R&D, Departement AMA et LaMSID UMR EDF-CNRS-CEA 8193, dina.al-akhrass@edf.fr, sebastien.fayolle@edf.fr</footnote>
de dissiper de l'énergie lors d'un cycle de charge/décharge élastique. Dans les années 80, des formalismes basés sur des concepts hyperélastiques, et ne présentant pas ce problème, ont émergé. Parmi ces derniers, il existe une catégorie utilisant le tenseur des déformations logarithmiques, qui ont l'avantage de permettre la décomposition additive des déformations en partie élastique et plastique. Celui que l'on considère ici fait partie de cette catégorie [7]. Il présente de nombreux avantages, dont l'un des principaux est l'extension naturelle des algorithmes d'intégration des lois de comportement, du cadre des petites vers les grandes déformations.

Dans ce papier, on se propose d'adapter ce formalisme de grandes déformations à une formulation mixte à trois, puis deux champs. Le formalisme de grandes déformations est d'abord décrit dans la section 2, puis la formulation mixte à trois et deux champs sont présentées dans la section 3. Enfin dans la section 4, des applications, dont une industrielle, sont présentées.

## 2 Description du formalisme de grandes déformations

Afin de mettre en place le cadre des grandes déformations, un modèle de type hyperélastique développé par Miehe, Apel et Lambrecht est considéré. La mise en place de ce formalisme se fait en trois étapes.

### Prétraitement géométrique

Cette étape consiste à définir le tenseur des déformations logarithmiques $\underline{E}$. Un point clé de la mise en place du cadre de la plasticité finie est la définition de la mesure de déformations élastiques $\underline{E}^e$. On suppose qu'elle est fonction du tenseur des déformations de Cauchy-Green droit $\underline{C}$ et du tenseur de métrique plastique $\underline{G}^p$ :

$$ \underline{E}^e = \underline{E}^e(\underline{C}, \underline{G}^p) \qquad (1) $$

On utilise la décomposition additive du tenseur des déformations

$$ \underline{E}^e = \underline{E} - \underline{E}^p \qquad (2) $$

où le tenseur des déformations logarithmiques $\underline{E}$ est tel que

$$ \underline{E} = \frac{1}{2} \ln(\underline{F}^T \underline{F}) = \frac{1}{2} \ln(\underline{C}) \qquad (3) $$

avec $\underline{F}$ le gradient de déformation.

Le tenseur des déformations logarithmiques plastiques est tel que

$$ \underline{E}^p = \frac{1}{2} \ln(\underline{G}^p) \qquad (4) $$

Ainsi, le tenseur des déformations logarithmiques permet de passer des caractéristiques multiplicatives de l'élasto-plasticité en grandes déformations, à la structure additive de la théorie de la géométrie linéaire. Le tenseur des contraintes $\underline{T}$ qui est conjugué au tenseur $\underline{E}$ est relié au second tenseur des contraintes de Piola-Kirchhoff $\underline{S}$ par la relation

$$ \underline{T} = \underline{S} : \underline{P}^{-1} \qquad (5) $$

avec

$$ \underline{P} \equiv_L = 2 \frac{\partial \underline{E}}{\partial \underline{C}} \qquad (6) $$
## Modèle constitutif dans l'espace des déformations logarithmiques

Comme nous l'avons vu ci-dessus, le choix de $\underline{E}$ comme tenseur des déformations nous permet de faire une décomposition additive des déformations en partie élastique et plastique. La loi de comportement va nous fournir, à partir du tenseur des déformations logarithmiques $\underline{E}$, des déformations plastiques $\underline{E}^p$, et des variables internes $\alpha$, le tenseur des contraintes $\underline{T}$ et le module tangent élasto-plastique associé $\underline{\mathbb{E}}^{ep}$.

$$ \{\underline{E}, \underline{E}^p, \alpha\} \Rightarrow \text{Loi de comportement} \Rightarrow \{\underline{T}, \underline{\mathbb{E}}^{ep}\} \quad (7) $$

Le module tangent détermine le taux de contrainte en respectant le taux de déformation logarithmique

$$ \dot{\underline{T}} = \underline{\mathbb{E}}^{ep} : \dot{\underline{E}} \quad (8) $$

## Post-traitement géométrique

Cette étape consiste à retrouver les tenseurs standards et leur module tangent associé, à l'aide des transformations géométriques introduites au dessus. Ainsi, le second tenseur de Piola-Kirchhoff $\underline{S}$ est obtenu par la relation (5), et le module tangent $\underline{\mathbb{C}}^{ep}$ par

$$ \underline{\mathbb{C}}^{ep}_{\equiv L} = 2 \frac{\partial \underline{S}}{\partial \underline{C}} = \left( 2 \frac{\partial \underline{E}}{\partial \underline{C}} \right)^T : \underline{\mathbb{E}}^{ep} : \left( 2 \frac{\partial \underline{E}}{\partial \underline{C}} \right) + \underline{T} : \left( 4 \frac{\partial^2 \underline{E}}{\partial \underline{C}^2} \right) \quad (9) $$

# 3 Formulation éléments finis mixtes

On se propose ici d'adapter le formalisme de grandes déformations présenté ci-dessus à une formulation mixte à trois champs puis à deux champs.

## 3.1 Formulation mixte à trois champs

Afin de mettre en place la formulation mixte à trois champs, le gradient de déformation $\underline{F}$ est enrichi par l'introduction de la variable de gonflement notée $\theta$, qui permet de mesurer la variation de volume et qui est reliée de façon faible au Jacobien $J$. Plusieurs relations sont possibles, parmi lesquelles nous avons choisi

$$ \ln(J) = \theta \quad (10) $$

et on définit alors le gradient de déformation enrichi comme

$$ \underline{\tilde{F}} = \left( \frac{\exp(\theta)}{J} \right)^{1/3} \underline{F} \quad (11) $$

Nous avons fait ce choix de fonction, car en plus de permettre la simplification de l'expression des équations, il est facile de montrer qu'avec cette forme de gradient de déformation enrichi, $\theta$ correspond à la partie sphérique de la déformation logarithmique enrichie [2].

L'équilibre de la structure est obtenu par le minimum de l'énergie potentielle $E(\mathbf{u}, \theta)$

$$ E(\mathbf{u}, \theta) = \int_{\Omega_0} \Phi(\underline{\tilde{C}}') d\Omega_0 - W_{\text{ext}}(\mathbf{u}) \quad (12) $$

avec $W_{\text{ext}}$ le potentiel des efforts extérieurs.

On introduit alors un multiplicateur de Lagrange $p$ qui nous permet d'imposer de façon faible, la relation entre $J$ et $\theta$ (10). On recherche alors le point selle du Lagrangien $Q$, qui est tel que

$$ Q(\mathbf{u}, \theta, p) = E(\mathbf{u}, \theta) + \int_{\Omega_0} p(\ln(J) - \theta) d\Omega_0 \quad (13) $$

Le point selle est caractérisé par la condition d'optimalité
$$ \forall \delta \mathbf{u}, \delta p, \delta \theta, \quad \delta Q = 0 \qquad (14) $$

La première variation de Q s'écrit

$$ \delta Q = \int_{\Omega_0} \left[ \frac{1}{2} \tilde{S} : \delta \tilde{C} + p \left( \frac{\delta J}{J} - \delta \theta \right) + \delta p (\ln(J) - \theta) \right] d\Omega_0 - \delta W_{\text{ext}} \quad (15) $$

On introduit la variation du gradient eulérien des déplacements $\delta \underline{L}$ qui est tel que

$$ \delta \underline{L} = \frac{\partial \delta \mathbf{u}}{\partial \mathbf{x}} = \delta \underline{F} \underline{F}^{-1} \qquad (16) $$

Finalement, la variation de Q peut s'écrire

$$ \delta Q = \int_{\Omega_0} \left[ (\tilde{\tau}^D + p \underline{11}) : \delta \underline{L} + \delta \theta \left( \frac{1}{3} \mathrm{tr}(\tilde{\tau}) - p \right) + \delta p (\ln(J) - \theta) \right] d\Omega_0 - \delta W_{\text{ext}} \quad (17) $$

avec $\tilde{\underline{\tau}} = \underline{\tilde{F}} \underline{\tilde{S}} \underline{\tilde{F}}^T$ et $\tilde{\underline{\tau}}^D$ sa partie déviatorique.

Ce qui nous donne les équations d'équilibre suivantes

$$ \left\{ \begin{array}{l} \delta Q^u = \displaystyle\int_{\Omega_0} (\tilde{\underline{\tau}}^D + p \underline{11}) : \delta \underline{L} d\Omega_0 - \delta W_{\text{ext}} = 0 \\ \\ \delta Q^\theta = \displaystyle\int_{\Omega_0} \delta \theta \left( \frac{1}{3} \mathrm{tr}(\tilde{\underline{\tau}}) - p \right) d\Omega_0 = 0 \\ \\ \delta Q^p = \displaystyle\int_{\Omega_0} \delta p (\ln(J) - \theta) d\Omega_0 = 0 \end{array} \right. \qquad (18) $$

Afin de résoudre ces équations, on utilise la méthode des éléments finis. Les espaces continus sont approchés par des espaces discrétisés appropriés, dans lesquels se trouvent les solutions approchées ($\mathbf{u}_h, \theta_h, p_h$). Le domaine $\Omega_0$ est discrétisé par le maillage $T_h(\Omega_0)$. Les trois champs sont discrétisés de sorte que

$$ \mathbf{u}_h = \sum_{k=1}^{nbnodes} \sum_{i=1}^{dim} N_{ki}^u U_{ki} \qquad (19) $$

$$ \theta_h = \sum_{k=1}^{nbnodes} N_k^\theta \Theta_k \qquad (20) $$

$$ p_h = \sum_{k=1}^{nbnodes} N_k^p P_k \qquad (21) $$

avec $\mathbf{N}^u, \mathbf{N}^\theta$ et $\mathbf{N}^p$ les fonctions de forme pour respectivement le champ des déplacements, de gonflement et de pression, et $U, \Theta$ et $P$ les variables nodales. Ainsi, pour chaque élément, les résidus sont tels que

$$ \left\{ \begin{array}{l} \delta Q_h^u = \displaystyle\int_{\Omega_e} (\tilde{\underline{\tau}}^D + P \underline{11}) : \underline{B} d\Omega_e - \delta W_{\text{ext}} = 0 \\ \\ \delta Q_h^\theta = \displaystyle\int_{\Omega_e} \mathbf{N}^\theta \left( \frac{1}{3} \mathrm{tr}(\tilde{\underline{\tau}}) - P \mathbf{N}^p \right) d\Omega_e = 0 \\ \\ \delta Q_h^p = \displaystyle\int_{\Omega_e} \mathbf{N}^p (\ln(J) - \Theta \mathbf{N}^\theta) d\Omega_e = 0 \end{array} \right. \qquad (22) $$

avec $\underline{B}$ le gradient des fonctions de forme tel que

$$ \delta L_{ij} = \sum_{k=1}^{nbnoeuds} \delta U_{ki} B_{ij}^k, \quad B_{ij}^k = \frac{\partial N_{ki}^u}{\partial x_j} \qquad (23) $$
Afin de résoudre l'équation le système (22), nous utilisons la méthode de Newton-Raphson. Dans cette méthode, le système d'équations linéaires à résoudre à la (i + 1)ième itération d'équilibre du (n + 1)ième pas de temps est

$$ \left( \begin{array}{ccc} \underline{K}^{uu} & \underline{K}^{u\theta} & \underline{K}^{up} \\ \underline{K}^{\theta u} & \underline{K}^{\theta\theta} & \underline{K}^{\theta p} \\ \underline{K}^{pu} & \underline{K}^{p\theta} & \underline{K}^{pp} \end{array} \right)^{(n+1,i)} \left( \begin{array}{c} \delta U \\ \delta \Theta \\ \delta P \end{array} \right)^{(n+1,i+1)} = - \left( \begin{array}{c} \delta Q_h^u \\ \delta Q_h^\theta \\ \delta Q_h^p \end{array} \right)^{(n+1,i)} \quad (24) $$

avec les termes de la matrice tangente qui sont tels que

$$ \begin{aligned} K_{inlm}^{uu} = & \int_{\Omega_0} [-B_{nj}B_{mj'}(\tilde{\tau}_{ij'}^D + P\delta_{ij'})\delta_{lj}] d\Omega_0 \\ & + \int_{\Omega_0} B_{nj}B_{mj'} \left[ \tilde{\tau}_{lj}\delta_{j'i} + \tilde{\tau}_{jj'}\delta_{li} - \frac{2}{3}\tilde{\tau}_{lj'}\delta_{ij} - \frac{2}{3}\tilde{\tau}_{ij}\delta_{lj'} + \frac{2}{9}tr(\tilde{\tau})\delta_{lj'}\delta_{ij} \right] d\Omega_0 \end{aligned} \quad (25) \\ & + 2 \int_{\Omega_0} B_{nj} (D_{ijlj'}^D)^D B_{mj'} d\Omega_0 \quad (26) $$

$$ K_{inr}^{u\theta} = \int_{\Omega_0} \frac{2}{3} N_r^\theta (B_{nj} \tilde{\tau}_{ij}^D + B_{nj} D_{ijkl}^D \delta_{kl}) d\Omega_0 = (K_{inr}^{\theta u})^T \quad (26) $$

$$ K_{ins}^{up} = \int_{\Omega_0} N_s^p B_{nj} \delta_{ij} d\Omega_0 = (K_{ins}^{pu})^T \quad (27) $$

$$ K_{rs}^{\theta p} = - \int_{\Omega_0} N_r^\theta N_s^p d\Omega_0 = (K_{sr}^{p\theta})^T \quad (28) $$

$$ K_{rt}^{\theta\theta} = \int_{\Omega_0} N_r^\theta N_t^\theta (tr(\underline{\tilde{\tau}}) + \delta_{ij} D_{ijkl} \delta_{kl}) d\Omega_0 \quad (29) $$

$$ K_{sv}^{pp} = 0 \quad (30) $$

où

$$ \underline{D} = \frac{1}{2} \underline{\tilde{F}} \underline{\tilde{F}} \underline{\underline{C}}^{ep} \underline{\underline{F}}^{T} \underline{\underline{F}}^{T} $$

## 3.2 Formulation mixte à deux champs

Il existe de nombreuses lois de comportement pour lesquelles l'incompressibilité plastique est prise en compte de manière naturelle. Ce sont les lois pour lesquelles le critère de plasticité ne dépend pas de la pression hydrostatique, ce qui est par exemple le cas pour le critère de Von Mises ou encore Tresca. Il existe alors une relation bi-univoque entre la pression et le gonflement et ce dernier peut alors être éliminé du champ des inconnus, de sorte que l'on se retrouve alors avec une formulation mixte à deux champs [10]. On définit alors un gradient de déformation enrichi qui est tel que

$$ \underline{\tilde{F}} = J^{-1/3} \underline{\underline{F}} \quad (31) $$

L'équilibre de la structure est ainsi régi par le minimum de l'énergie potentielle $E(\mathbf{u})$

$$ E(\mathbf{u}) = \int_{\Omega_0} \Phi(\underline{\underline{C}}) - W_{ext} \quad (32) $$

Afin de mettre en place la formulation à deux champs, on introduit un multiplicateur de Lagrange $p$ qui renforce la contrainte d'incompressibilité, qui dans ce cas s'écrit

$$ \ln(J) = \frac{p}{2\kappa} \quad (33) $$

avec $\kappa$ le module de compressibilité.
$$Q_L(\mathbf{u}, p) = E(\mathbf{u}) + \int_{\Omega_0} p(\ln(J) - \frac{p}{2\kappa})d\Omega_0 \quad (34)$$

Le point selle est caractérisé par les conditions d'optimalité

$$\delta Q = 0, \quad \forall \delta \mathbf{u}, \delta p \qquad (35)$$

La première variation de Q s'écrit

$$\delta Q = \int_{\Omega_0} \frac{1}{2} \tilde{S} : \delta \underline{C} d\Omega_0 + \int_{\Omega_0} ptr(\delta \underline{L}) d\Omega_0 + \int_{\Omega_0} \delta p \left( \ln(J) - \frac{p}{\kappa} \right) d\Omega_0 - \delta W_{ext} = 0 \quad (36)$$

Il est facile de montrer que l'on peut l'écrire sous la forme

$$\delta Q = \int_{\Omega_0} (\tilde{\underline{\tau}}^D + p \underline{11}) : \delta \underline{L} d\Omega_0 + \int_{\Omega_0} \delta p \left( \ln(J) - \frac{p}{\kappa} \right) d\Omega_0 - \delta W_{ext} = 0 \quad (37)$$

Ceci nous donne les équations d'équilibre suivantes :

$$\left\{ \begin{array}{l} \delta Q^u = \displaystyle\int_{\Omega_0} (\tilde{\underline{\tau}}^D + p \underline{11}) : \delta \underline{L} d\Omega_0 - \delta W_{ext} = 0 \\ \\ \delta Q^p = \displaystyle\int_{\Omega_0} \delta p \left( \ln(J) - \frac{p}{\kappa} \right) d\Omega_0 = 0 \end{array} \right. \quad (38)$$

Comme pour la formulation mixte à trois champs, les deux champs u et p sont alors discrétisés (19) et (21)), de même que les résidus. De la même façon, la méthode de Newton-Raphson est utilisée pour résoudre le problème (35), de sorte qu'à la $(i+1)$ième itération d'équilibre du $(n+1)$ième pas de temps, il faut résoudre le système

$$\left( \frac{\underline{K}^{uu}}{\underline{K}^{pu}} \quad \frac{\underline{K}^{up}}{\underline{K}^{pp}} \right)^{(n+1,i)} \left( \begin{array}{c} \delta U \\ \delta P \end{array} \right)^{(n+1,i+1)} = - \left( \begin{array}{c} \delta Q^u \\ \delta Q^p \end{array} \right)^{(n+1,i)} \quad (39)$$

avec les termes de la matrice tangente qui sont tels que

$$\begin{aligned} K_{inlm}^{uu} = & \int_{\Omega_0} [-B_{nj}B_{mj'}(\tilde{\tau}_{ij'}^D + P\delta_{ij'})\delta_{lj}] d\Omega_0 \\ & + \int_{\Omega_0} B_{nj}B_{mj'} \left[ \tilde{\tau}_{lj}\delta_{j'i} + \tilde{\tau}_{jj'}\delta_{li} - \frac{2}{3}\tilde{\tau}_{lj'}\delta_{ij} - \frac{2}{3}\tilde{\tau}_{ij}\delta_{lj'} + \frac{2}{9}tr(\tilde{\tau})\delta_{lj'}\delta_{ij} \right] d\Omega_0 \end{aligned} \quad (40) \\ & + 2 \int_{\Omega_0} B_{nj} (D_{ijlj'}^D)^D B_{mj'} d\Omega_0 \end{aligned}$$

$$K_{ins}^{up} = \int_{\Omega_0} N_s^p B_{nj} \delta_{ij} d\Omega_0 = (K_{ins}^{pu})^T \quad (41)$$

$$K_{rs}^{pp} = - \int_{\Omega_0} \frac{N_r^p N_s^p}{\kappa} d\Omega_0 \quad (42)$$

# 4 Résultats numériques

## 4.1 Cube déformé

Le premier cas test considéré ici consiste en un cube de 10mm de côté ayant un comportement élasto-plastique parfaitement plastique (figure 1). Il est soumis à une pression de 1000Pa exercée sur un quart de sa face supérieure. Le module de Young est de $E = 200GPa$ et le coefficient de Poisson de $\nu = 0.3$. Le critère de plasticité considéré est celui de Von Mises et l'écrouissage est isotrope linéaire avec une limite d'élasticité $\sigma_Y = 150MPa$. La figure 1 représente la répartition de la trace de la contrainte aux points de Gauss avec des éléments linéaires, quadratiques, P2/P1/P1 (issus de la formulation à trois
<figure><img src="image_2.png" /><figcaption>Fig. 1 – Répartition de la trace de la contrainte aux points de Gauss (MPa) avec des éléments : (a) P1 (b) P2 (c) P2/P1/P1 (d) P2/P1</figcaption></figure>

champs, quadratique en déplacements, linéaire en pression et gonflement) et P2/P1 (pour celle à deux champs).

Ainsi, on voit apparaitre d'importantes oscillations avec les éléments linéaires (on a une répartition en damier). Avec les éléments quadratiques, on observe également des oscillations, un peu moins importantes, en particulier au niveau de la zone où la pression est exercée. Cependant avec des éléments P2/P1/P2 et P2/P1 la répartition est régulière et identique.

## 4.2 Cas industriel : Evaluation des contraintes au niveau de la zone d'affouillement d'un organe d'admission

Il s'agit d'estimer la durée de vie en fluage d'un organe d'admission pour une centrale thermique en acier 15CDV910, soumis à une pression interne constante et une température variant entre 475 et 565 degrés Celsius. Cet organe présente une fissure et afin d'éviter sa propagation, un affouillement est réalisé au niveau de celle-ci. Le matériau est élasto-visco-plastique et suit une loi de Lemaitre [5]. La figure 2 représente la géométrie de l'organe d'admission, avec la zone affouillée entourée d'un cercle.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Géométrie de l'organe d'admission</figcaption></figure>

Au niveau de la zone d'affouillement, il y a apparition d'oscillations sur les valeurs de pression dues à l'incompressibilité, même avec des éléments sous-intégrés, comme le montre la figure 3(a). La figure 3(b) montre la répartition de la pression obtenue avec des éléments P2/P1/P1 ; il apparaît clairement qu'il n'y a plus d'oscillations et que la répartition est régulière. L'utilisation d'éléments sous-intégrés n'est pas suffisante pour obtenir des résultats fiables et précis.
<figure><img src="image_4.png" /><figcaption>Fig. 3 – Repartition de la pression hydrostatique aux points de Gauss au niveau de la zone d'affouillement avec des éléments : (a) sous-intégrés (b) P2/P1/P1</figcaption></figure>

**5 Conclusions**

Dans ce papier, un formalisme de grandes déformations basé sur les déformations logarithmiques, est adapté à une formulation mixte à trois champs, puis deux champs. Les exemples numériques montrent que les résultats obtenus ne présentent pas d'oscillations sur les valeurs de pression et de blocage volumétrique et que le phénomène d'incompressibilité est bien géré. L'adaptation du formalisme aux méthodes de stabilisation (mini-élément et OSGS) est en cours.

**Références**

[1] F. Brezzi, M. Fortin. *Mixed and Hybrid Finite Element*, Springer-Verlag, 1991.

[2] M. Brunig. *Formulation and numerical treatment of incompressibility constraints in large strain elastic-plastic analysis*, International Journal for Numerical Methods in Engineering, 1047-1068, 1999.

[3] M. Chiumenti, Q. Valverde, C. Agelet de Saracibar, M. Cervera. *A stabilized formulation for incompressible plasticity using linear triangles and tetrahedra*, International Journal of Plasticity, 1487-1504, 2004.

[4] T.J.R Hughes. *The Finite Element Method: Linear Static and Dynamic Finite Element Analysis*, Prentice-Hall, 1987.

[5] E. Galenne. *Comportement de l'acier soumis à la corrosion*, Documentation Code_Aster, 2009.

[6] E.P Kasper, R.L Taylor. *A mixed-enhanced strain method Part I : Geometrically linear problems*, Computers and structures, 237-250, 2000.

[7] C. Miehe, N. Apel, M. Lambrecht. *Anisotropic additive plasticity in the logarithm strain space modular kinematic formulation and implementation based on incremental minimization principles for standard materials*, Computational Methods in Applied Mechanics and Engineering, 5383-5425, 2002.

[8] J.C Natgtegaal, J.E. de Jong. *Some computational aspects of elastic-plastic large strain analysis*, International Journal of Numerical Methods in Engineering, 15-41, 1981.

[9] D.P. Recio, R.M. Natal Jorge, L.M.S. Dinis. *Locking and hourglass phenomena in an element-free Galerkin context : the B-bar method with stabilization and an enhanced strain method*, International Journal for Numerical Methods in Engineering, 1329-1357, 2006.

[10] T. Sussman, K.J Bathe. *A finite element formulation for nonlinear incompressible elastic and inelastic analysis*, Computers & Structures, 357-409, 1987.

[11] R.L Taylor. *A mixed-enhanced formulation for tetrahedral finite elements*, International Journal for Numerical Methods in Engineering, Wiley, 205-227, 2000.