# Préservation de l'orientation en hyperélasticité compressible

François Peyraut, Zhi-Qiang Feng, Nadia Labed

► To cite this version:

François Peyraut, Zhi-Qiang Feng, Nadia Labed. Préservation de l'orientation en hyperélasticité compressible. 8e Colloque national en calcul des structures, CSMA, May 2007, Giens, France. hal-01495784

HAL Id: hal-01495784

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

Submitted on 26 Mar 2017

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
# Préservation de l'orientation en hyperélasticité compressible

François Peyraut * — Zhi Qiang Feng ** — Nadia Labed ***

* Laboratoire de Recherche et d'Etude sur les Matériaux, les Procédés et les Surfaces (LERMPS), Université de Technologie de Belfort-Montbéliard
90010 Belfort cedex, France

{francois.peyraut, nadia.labed}@utbm.fr

** Laboratoire de Mécanique d'Evry, Université d'Evry-Val d'Essonne
40 rue du Pelvoux, 91020 Evry, France

feng@iup.univ-evry.fr

*** Laboratoire M3M, Université de Technologie de Belfort-Montbéliard
90010 Belfort cedex, France

RÉSUMÉ. Cette communication met l'accent sur le développement d'un algorithme destiné à garantir numériquement la préservation de l'orientation dans le cadre de l'analyse éléments finis des grandes déformations hyperélastiques. Cet algorithme consiste à empêcher l'apparition de valeurs propres négatives de la matrice gradient des déformations en contrôlant le pas de charge de la méthode de Newton-Raphson. Plusieurs exemples numériques sont présentés pour illustrer l'efficacité de cet algorithme. Ces exemples impliquent des conditions aux limites et des lois de comportement variées (lois d'Ogden, de Gent, de Blatz-Ko). Un cas tridimensionnel impliquant un maillage de plusieurs milliers d'éléments conduisant à des déformations non homogènes est également étudié.

ABSTRACT. This paper is focused on the orientation preserving condition in the case of large hyperelastic deformations. A numerical algorithm ensuring this condition by the finite element method is presented. This algorithm prevents the occurrence of negative eigenvalues of the deformation gradient matrix by selecting the appropriate load step number for the Newton-Raphson method. Several numerical examples illustrate the efficiency of the method. Different boundary conditions and energy densities (Ogden, Gent, and Blatz-Ko models) are considered. A three dimensional case, involving thousands of elements with non-homogeneous deformation, is also studied.

MOTS-CLÉS : hyperélasticité, condition de préservation de l'orientation, méthode des éléments finis.

KEYWORDS: hyperelastic material, orientation preserving condition, finite element method.
# 1. Introduction

Cette communication se propose de faire le point sur les problèmes de préservation de l'orientation dans le cadre de l'hyperélasticité compressible. Dans certaines conditions de calcul, il a en effet été établi que cette condition n'est pas satisfaite, ce qui conduit à la divergence de la méthode de Newton-Raphson (Peyraut et al., 2001). Un algorithme basé sur l'analyse préalable des valeurs propres de la matrice gradient des déformations a été proposé pour pallier ce problème (Peyraut et al., 2004). Comme les valeurs propres sont évaluées par les formules de Cardan sur des matrices de petite taille (3x3), le coût de cet algorithme est marginal en comparaison du temps CPU requis pour réaliser le calcul éléments finis. Il a été implémenté dans le code de calcul aux éléments finis FER, développé par le Laboratoire de Mécanique de l'Université d'Evry (http://gmfe16.cemif.univ-evry.fr:8080/~feng/FerSystem.html). Son efficacité a été testée avec des exemples académiques utilisant différente combinaison de géométries et de conditions aux limites (Peyraut, 2003 ; Peyraut et al., 2003). Plusieurs lois de comportement ont été utilisées : densités de Blatz-Ko (Blatz et al., 1962), d'Ogden (Ogden, 1972) et de Gent (Gent, 1996), modèle « 3é » (Lainé et al., 1999). Un exemple tridimensionnel impliquant des déformations non homogènes a également été testé. Enfin, des problèmes de contact et d'impact en grandes déformations hyperélastiques ont été traités avec les modèles de Blatz-Ko (Feng et al., 2003), d'Ogden (Feng et al., 2006a) et avec le modèle « 3é » (Feng et al., 2006b).

# 2. Hyperélasticité compressible

L'hyperélasticité postule l'existence d'une densité d'énergie de déformation W fonction du tenseur des déformations de Green-Lagrange E [1a]. Les plus courantes concernent les matériaux isotropes pour lesquels W ne dépend que des trois invariants de la matrice de dilatation $C=F^T F$ ($I_1$=Trace($C$) ; $I_2$=(I₁²-Trace($C$²))/2 ; $I_3$=det($C$)) ou de ses trois valeurs propres $λ_1, λ_2$ et $λ_3$. La matrice gradient des déformations F, qui s'exprime en fonction des coordonnées lagrangienne a et du déplacement u, décrit la transformation géométrique entre les configurations non déformée et déformée [1b]. Son déterminant [1c] représente le rapport entre deux volumes élémentaires $dV_0$ et $dV$ pris chacun dans ces deux configurations. La positivité de ce déterminant traduit le principe de préservation de l'orientation. Trois densités différentes sont utilisées dans cet article pour les applications numériques : modèles de Blatz-Ko [2], d'Ogden [3] et de Gent [4].

$$ E = \frac{C - I}{2} \quad \text{a) } F(a) = I + \frac{\partial u}{\partial a} \quad \text{b) } \det(F) = \frac{dV}{dV_0} > 0 \quad \text{c)} \qquad [1] $$

$$ W = \frac{\mu}{2} \left\{ \frac{I_2}{I_3} + 2I_3^{1/2} - 5 \right\} \qquad [2] $$
$$
W = \sum_{i=1}^{N} \frac{\mu_i}{\alpha_i} (\lambda_1^{\alpha_i/2} + \lambda_2^{\alpha_i/2} + \lambda_3^{\alpha_i/2} - 3) + \sum_{i=1}^{N} \frac{\mu_i}{\alpha_i \beta_i} (I_3^{-\alpha_i} \beta_i / 2 - 1) \quad [3]
$$

$$
W(I_1^*, J) = -\frac{\mu}{2} J_m \ln \left( 1 - \frac{I_1^* - 3}{J_m} \right) + \frac{1}{d} \left( \frac{J^2 - 1}{2} - \ln J \right) \quad [4]
$$

Où μ représente le module de cisaillement. N, αᵢ et βᵢ sont des paramètres du modèle d'Ogden. Jₘ est une constante et d un paramètre matériau caractérisant l'incompressibilité du modèle de Gent. Enfin, dans le modèle de Gent, l'invariant I₁ est remplacé par I₁* = I₃⁻¹/³ I₁ de manière à isoler la partie isochore des déformations.

3. Condition de préservation de l'orientation

Pour traiter numériquement les problèmes de préservation de l'orientation, un
algorithme basé sur l'analyse des valeurs propres de F a été développé (Peyraut et
al., 2004). Le principe consiste à choisir le paramètre de chargement α de la
méthode de Newton-Raphson incrémental de manière à prévenir l'apparition de
valeurs propres négatives. Pour cela, on considère le déplacement **u**_α obtenu à la
première itération par la résolution du système linéaire tangent [5a] :

$$
K_T u_\alpha = \alpha f \quad \text{a) ;} \quad \lambda_\alpha^i = (\lambda_1^i - 1)\alpha + 1, \quad i=1,2,3 \quad \text{b)} \tag{5}
$$

Où K_T représente la matrice de raideur tangente et f les forces extérieures. Il est facile d'établir, en utilisant la linéarité de [5a] par rapport à α, que les valeurs propres λ^i_α (i=1,2,3) de F_α = Id + ∇(u_α) dépendent linéairement de α [5b]. Cette dépendance linéaire est illustrée sur la figure 1. Dans le cas où une valeur propre de F_1 est négative, par exemple λ^i_1 < 0, la droite définie par [5b] décroit de 1 à λ^i_1 et intersecte l'axe horizontal au point d'abscisse α_0 définie par [6a]. Comme λ^i_1 est négatif, α_0 est strictement compris entre 0 et 1. L'intervalle ]0,1] se divise alors en deux parties. Dans la première, définie par ]0, α_0[, λ^i_α est positif et l'orientation est préservée. Dans la seconde, définie par [α_0,1], λ^i_α est négatif et l'orientation n'est pas préservée. Le paramètre de chargement α doit donc être choisi dans l'intervalle ]0, α_0[. Il doit par ailleurs être le plus proche possible de α_0, sans l'atteindre, pour réduire le temps de calcul. On utilise alors un nombre de pas de charge entier n [6b], plutôt qu'un paramètre de chargement réel α□, ce qui correspond à la pratique usuelle dans les codes de calcul aux éléments finis. Ce nombre de pas de charge répond au double objectif de préserver l'orientation et de minimiser le temps de calcul.
$$ \alpha_0 = \frac{1}{1 - \lambda_1^i} \quad \text{a) } n = \text{partie entière de } \left(\frac{1}{\alpha_0}\right) + 1 \quad \text{b) } \qquad [6] $$

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Paramètre optimal $\alpha_0$</figcaption></figure>

## 4. Application numérique

### 4.1. Déformation homogène

Cette application concerne un cube de 0.5 m de côté chargé par une pression uniforme $p=1$ MPa sur sa face supérieure (fig. 2). Sa base est encastrée et ses quatre faces latérales simplement supportées. Ces conditions aux limites induisant des déformations homogènes, un seul élément brique à 8 nœuds a été utilisé pour réaliser le maillage. Trois lois de comportement différentes ont été considérées. Les paramètres matériaux de chacune sont donnés dans le tableau 1. Un excellent accord a été observé entre le nombre de pas optimal $n$ prévu par [6b] et des calculs éléments finis réalisés avec FER (tableau 2). FER converge en effet si $n$ est utilisé alors qu'il diverge pour $n-1$.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Maillage</figcaption></figure>

<table><caption>**Tableau 1.** Paramètres matériaux</caption><thead><tr><th>Blatz-Ko [2]</th><th>Ogden [3]</th><th>Gent [4]</th></tr></thead><tbody><tr><td>μ=220711.2 Pa</td><td>N=2 (modèle à 2 paramètre)</td><td>μ=3E5 Pa</td></tr><tr><td></td><td>μ<sub>1</sub>=1.85E3 Pa ; μ<sub>2</sub>=-9.2E3 Pa</td><td>J<sub>m</sub>=97</td></tr><tr><td></td><td>α<sub>1</sub>=4.5 ; α<sub>2</sub>=-4.5 ; β<sub>1</sub>=β<sub>2</sub>=1</td><td>d=0.001 Pa<sup>-1</sup></td></tr></tbody></table>
<table><caption>**Tableau 2.** Prédiction FER – comparaison avec des calculs FER et ANSYS</caption><thead><tr><th>Modèles</th><th>Blatz-Ko</th><th>Ogden</th><th>Gent</th></tr></thead><tbody><tr><td>n optimal [6b]</td><td>2</td><td>11</td><td>3</td></tr><tr><td>Calcul FER avec n-1</td><td>divergence</td><td>divergence</td><td>divergence</td></tr><tr><td>Calcul FER avec n</td><td>convergence</td><td>convergence</td><td>convergence</td></tr></tbody></table>

## 4.2. Déformation non homogène

Cet exemple concerne un cas tridimensionnel impliquant des déformations non homogènes. Il s'agit d'un cube de 1 m de côté, chargé sur sa face supérieure par une pression uniforme de 1 MPa et maillé avec 8000 éléments brique à 8 nœuds (fig. 3). La base du cube est encastrée et sa face supérieure est supposée indéformable. Le modèle de Blatz-Ko est considéré pour cet exemple avec les propriétés matérielles du tableau 1. L'équation [6b] indique que le nombre de pas de charge *n* doit être égal à 3. Cette prédiction est conforme aux résultats des calculs éléments finis (fig. 4). En effet, dans le cas où *n* est égal à 1, le cube traverse l'encastrement situé à z=0, tous les éléments déformés du maillage ayant un volume négatif. L'orientation n'est pas préservée. Si *n* est égal à 2, les éléments situés aux coins de la face inférieure du cube ont un volume négatif. A nouveau, l'orientation n'est pas préservée. Seuls quelques éléments sont concernés par un volume négatif car on est proche de la situation où l'orientation est préservée (*n*=3).

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Maillage</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Déformée à la première itération du premier pas de charge</figcaption></figure>
**5. Conclusions**

La problématique de la préservation de l'orientation dans le cadre de l'hyperélasticité compressible a été abordée. Un algorithme basé sur l'analyse des valeurs propres de la matrice gradient des déformations a été proposé. Il consiste à prévenir l'apparition de valeurs propres négatives en contrôlant le pas de charge de la méthode de Newton-Raphson. Deux exemples numériques, avec des déformations homogènes et non homogènes, illustrent son efficacité.

**6. Bibliographie**

Blatz P.J., Ko W.L., « Application of Finite Elastic Theory to the Deformation of Rubbery Materials », *Transactions of the Society of Rheology*, Vol. 6, 1962, p. 223-251.

Feng Z.-Q., Peyraut F., Labed N., « Solution of large deformation contact problems with friction between Blatz-Ko hyperelastic bodies », *International Journal of Engineering Science*, Vol. 41, 2003, p. 2213-2225.

Feng Z.-Q., Peyraut F., He Q.-C., « Finite deformations of Ogden's materials under impact loading », *International Journal of Non-Linear Mechanics*, Vol. 41, Issue 4, 2006a, p. 575-585.

Feng Z.-Q., Vallée C., Fortuné D., Peyraut F., « The 3é hyperelastic model applied to the modeling of 3D impact problems », *Finite Elements in Analysis and Design*, Vol. 43, Issue 1, 2006b, p. 51-58.

Gent A.N., « A new constitutive relation for rubber », *Rubber Chem. Technol.*, 1996, 69:59-61.

Lainé E., Vallée C., Fortuné D., « Nonlinear isotropic constitutive laws: choice of the three invariants, convex potentials and constitutive inequalities », *Int. J. Eng. Sci.*, 37, 1999, p. 1927-1941.

Ogden R.W., « Large deformation isotropic elasticity: on the correlation of theory and experiment for compressible rubberlike solids », *Proc. R. Soc. Lond. A.* 328, 1972, p. 567-583.

Peyraut F., Labed N., « Préservation de l'orientation et convergence de Newton-Raphson avec le modèle hyperélastique compressible de Blatz-Ko », *Revue Européenne des Eléments Finis*, 10, 2001, p. 595-605.

Peyraut F., « Orientation preservation and Newton-Raphson convergence in the case of an hyperelastic sphere subjected to an hydrostatic pressure », *Computer Methods in Applied Mechanics and Engineering*, Vol. 192, Nos 9-10, 2003, p. 1107-1117.

Peyraut F., Labed N., « Condition de préservation de l'orientation - application en hyperélasticité compressible », *Revue Européenne des Eléments Finis*, Vol. 12, N°1, 2003, p. 99-116.

Peyraut F., Feng Z.-Q., « Implémentation éléments finis d'une condition optimale de préservation de l'orientation. Cas du contact en grandes déformations hyperélastiques », *Revue Européenne des Eléments Finis*, Vol. 13, N°8, 2004, p. 899-914.