# Prise en compte d'une limite d'extensibilité en grandes transformations

Nicolas Bonfils, Nicolas Chevaugeon, Gilles Marckmann, Nicolas Moes

► To cite this version:

Nicolas Bonfils, Nicolas Chevaugeon, Gilles Marckmann, Nicolas Moes. Prise en compte d'une limite d'extensibilité en grandes transformations. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01422167

HAL Id: hal-01422167

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

Submitted on 23 Dec 2016

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 CC BY 4.0 - Attribution - International License
# Prise en compte d'une limite d'extensibilité en grandes transformations

N. Bonfils¹, N. Chevaugeon¹, G. Marckmann¹, N. Moës¹

¹ Institut de Recherche en Génie Civil et Mécanique (GeM) UMR CNRS 6183
Ecole Centrale de Nantes
1 rue de la Noë - BP.92101 - 44321 NANTES CEDEX 3
{nicolas.bonfils,nicolas.chevaugeon,gilles.marckmann,nicolas.moes}@ec-nantes.fr

**Résumé** — Une bonne modélisation du comportement des matériaux élastomères est essentielle pour aborder des phénomènes comme la fissuration ou l'endommagement. De nombreux modèles hyperélastiques existent, mais la plupart d'entre eux sont complexes et dépendent de nombreux paramètres. Dans cette étude, nous proposons d'enrichir un modèle classique et simple (néo-hookéen) en imposant une contrainte matérielle de limite d'extensibilité. Nous aboutirons donc à un problème de minimisation sous contraintes d'inégalités que nous résoudrons par une méthode type Lagrangien Augmenté.

**Mots clés** — limite d'extensibilité, grandes transformations, élastomères.

## 1 Introduction

La propagation de fissures dans les matériaux élastomères est un phénomène encore assez mal compris. La nature complexe du comportement semble être à l'origine de phénomène de bifurcation de fissure ("knotty tearing") non pris en compte dans les modèles numériques [9]. Des études récentes [12] ont montré l'apparition d'une zone rigide (cristalline) en fond de fissure dans les matériaux élastomères. Le lien entre la cristallisation et la bifurcation de fissure a également été mis en évidence [10, 11]. Nous avons voulu prendre en compte ce phénomène en imposant une limite d'extensibilité des chaînes. Deux méthodes sont possibles alors :

- La première consiste à prendre en compte cette limite dans la forme de la densité d'énergie (Modèle de Gent [7], ou le modèle 8 chaînes [1])

- La deuxième consiste à imposer la limite comme une contrainte matérielle.

Le différence essentielle entre les deux modèles est que le deuxième introduit une liaison interne cinématique dans le matériau (ou contrainte matérielle). L'introduction d'une telle contrainte (au même titre que l'incompressibilité) restreint l'espace admissible pour le tenseur de déformation. C'est cette méthode que nous avons choisi pour la représentation de la limite d'extensibilité.

## 2 Gestion de la contrainte matérielle

### 2.1 Généralités

De manière très générale, une contrainte matérielle simple s'écrit à l'aide d'une fonction à valeur scalaire $\gamma$ telle que $\gamma(\underline{F}, \underline{X}) = 0$ [13]. Cette relation interdit le tenseur gradient de la transformation $\underline{F}$ d'être arbitraire. Cette contrainte doit être invariante pour des mouvements de solides
rigides. On l'exprime donc préférablement en fonction de $\underline{C} = \underline{F}^T \underline{F} : \gamma(\underline{C}, \underline{X}) = 0$. Pour exemple, citons seulement quelques contraintes matérielles classiques :

- L'incompressibilité (déjà largement traitée dans la littérature) : $\gamma(\underline{C}, \underline{X}) = \det(\underline{C}) - 1 = 0$.

- L'inextensibilité dans une direction donnée $\underline{a} : \gamma(\underline{C}, \underline{X}) = \underline{a}\underline{C}\underline{a} - \underline{1} = 0$

- Le comportement des cristaux [5] : $\gamma(\underline{C}, \underline{X}) = \text{tr}(\underline{C}) - 3$.

## 2.2 Hyperélasticité contrainte

Le principe de déterminisme de Truesdell [13] autorise une décomposition a priori du tenseur PK-2 (Piola-Kirchhoff 2) en deux parties pour les matériaux contraints :

$\underline{S} = \underline{S}^r + \underline{S}_R$, avec $\underline{S}^r$ la partie active ou régulière des contraintes et $\underline{S}_R$ la partie réactive. La partie active dépend directement de l'énergie libre spécifique (déterminée constitutivement) et la partie réactive dépend des contraintes internes. Si on explicite ces dépendances :

$$ \underline{S}_R = \xi \frac{\partial \gamma(\underline{C})}{\partial \underline{C}} (\underline{\xi} \in \mathbb{R}) \text{ et } \underline{S}^r = 2 \frac{\partial \psi^r}{\partial \underline{C}} $$

## 3 Modèle hyperélastique considéré

Il existe dans la littérature un certain nombre de modèles hyperélastiques. Notre choix s'est rapidement orienté sur le modèle néo-hookéen compressible pour sa simplicité. Le potentiel néo-hookéen est une fonction de deux invariants tensoriels $I_1 = \text{tr}(\underline{C})$ et $I_3 = \det(\underline{C}) = J^2$ :

$$ \psi^r(\underline{C}) = \frac{\mu}{2}(I_1 - 3) - \mu \ln(J) + \frac{\lambda}{2}(\ln(J))^2 $$

où $\psi^r$ est le potentiel d'énergie libre régulier. En adoptant le formalisme Lagrangien, la partie régulière du tenseur des contraintes de Piola-Kirchhoff 2 s'écrit :

$$ \underline{S}^r = 2 \frac{\partial \psi^r}{\partial \underline{C}} = \mu(I - \underline{C}^{-1}) + \lambda \ln(J) \underline{C}^{-1} $$

Enfin la partie régulière du tenseur d'élasticité Lagrangien d'ordre 4 (issue de la linéarisation de la relation de comportement) :

$$ \mathbb{C}^r = 2 \frac{\partial \underline{S}^r}{\partial \underline{C}} = \lambda \underline{C}^{-1} \otimes \underline{C}^{-1} + 2(\mu - \lambda \ln(J)) \mathbb{J} \text{ avec } \mathbb{J} = - \frac{\partial \underline{C}^{-1}}{\underline{C}} $$

## 4 Contrainte de limite d'extensibilité et conditions de Kuhn-Tucker

Le cadre classique étant posé, il reste à introduire la contrainte de limite d'extensibilité des chaînes. Elle peut s'imposer de deux manières différentes :

1. Le restriction cinématique peut être représentée en terme d'élongations principales [8] : $\gamma_2(\underline{C}) = \max(\lambda_1, \lambda_2, \lambda_3) - \lambda_m \le 0$, où $\lambda_i$ sont les valeurs propres (élongations principales) du gradient de la transformation $\underline{F}$ et $\lambda_m$ la limite d'extensibilité.

2. Mais elle peut également s'imposer sur le premier invariant du tenseur des dilatations $\underline{C}$. On note le premier invariant $I_1 = \text{tr}(\underline{C}) = \lambda_i^2$. Si $I_m$ est la valeur limite pour le premier invariant, la contrainte cinématique peut s'écrire $\gamma_1(\underline{C}) = I_1 - I_m \le 0$ [8].
Remarque : Le potentiel total se note, en présence d'une telle contrainte :

$\psi(\underline{C}) = \psi^r(\underline{C}) + I_-(\gamma(\underline{C}, \underline{X}))$ avec $I_-(x)$, la fonction indicatrice de $x$, telle que :

$$I_-(x) = \begin{cases} 0 & \text{si } x \le 0 \\ +\infty & \text{sinon} \end{cases}$$

C'est $\gamma_2$ que nous retiendrons comme contrainte pour modéliser la limite d'extensibilité. Le tenseur des contraintes se décompose donc en deux parties :

$$\underline{S} = \underline{S}^r + q \frac{\partial \gamma_2(\underline{C})}{\partial \underline{C}} = \underline{S}^r + q \underline{I}$$

Avec les conditions de Kuhn-Tucker suivantes :

$$q \ge 0, \quad \mathrm{tr}(\underline{C}) - I_m \le 0, \quad q \mathrm{tr}(\underline{C}) = 0$$

on transforme maintenant ces conditions en une équation [2, 6]

$$q - \chi(g)g = 0$$

avec :

$$g = q + \rho_n(\mathrm{tr}(\underline{C}) - I_m), \quad (\rho_n \ge 0)$$

et :

$$\chi(x) = \begin{cases} 0 & \text{si } x < 0 \\ 1 & \text{si } x \ge 0 \end{cases}$$

À l'aide de cette nouvelle écriture on peut exprimer le tenseur des contraintes de la manière suivante :

$$\underline{S} = \underline{S}^r + \chi(g)g\underline{I}$$

Le tenseur d'élasticité total se décompose également en deux parties distinctes :

$$\underline{C} = 2 \frac{\partial \underline{S}^r}{\partial \underline{C}} + 2 \frac{\partial (\chi(g)g\underline{I})}{\partial \underline{C}}$$

# 5 Formulation forte

On considère un corps **B** de surface $\Omega$ et de frontière $\Gamma$ en déformation plane. Nous notons $C_O$ sa configuration non-déformée. Ce corps, soumis à des sollicitations planes, se déforme. Nous noterons alors $\omega$ la nouvelle surface qu'il occupe et $\gamma$ sa frontière. Il occupe alors la configuration déformée **C**. Nous faisons un partage de la frontière en deux. Une partie notée $\Gamma_u$ sur laquelle s'appliquent les conditions aux limites de type Dirichlet. Une autre partie notée $\Gamma_T$ sur laquelle s'appliquent les conditions aux limites de type Neumann.

Compte tenu de ces notations, nous pouvons exprimer la formulation forte du problème en description Lagrangienne.

$$\underline{\mathrm{Div}}(\underline{P}) + \underline{B} = 0 \quad \text{sur } \Omega$$

$$\underline{P} \cdot \underline{N} = \underline{T}_d \quad \text{sur } \Gamma_T$$

$$\underline{u} = \underline{U}_d \quad \text{sur } \Gamma_u$$

$$\underline{P} = J^{-1} \underline{F} \cdot \underline{S} = J^{-1} \underline{F} (\underline{S}^r + \chi(g)g\underline{I}) \quad \text{sur } \Omega$$

$$\underline{C} = \underline{F}^T \underline{F} \quad \text{sur } \Omega$$

$$q - \chi(g)g = 0 \quad \text{sur } \Omega$$

$$g = q + \rho_n(\mathrm{tr}(\underline{C} - I_m)) \quad \text{sur } \Omega$$
<figure><img src="image_2.png" /><figcaption>Figure 1 – Algorithme des contraintes actives</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2 – Cas test</figcaption></figure>

## 6 Formulation variationnelle mixte

Nous noterons $V^0$ l'ensemble des champs de déplacements cinématiquement admissibles et $Q$ l'ensemble des champs de multiplicateurs de Lagrange. La formulation faible à deux champs s'écrit :

Trouver le couple $(\underline{u}, q) \in V^0 \times Q$ tel que $\forall \delta\underline{v} \in V^0$ et $\forall \delta q \in Q$

$$ \int_{\Omega} \underline{S}(\underline{u}) : \delta \dot{\underline{E}} d\Omega = \int_{\Omega} \underline{B} \cdot \delta \underline{v} d\Omega + \int_{\Gamma_T} \underline{T}_d \cdot \delta \underline{v} dA = \delta W(\underline{\phi}, \delta \underline{v}) = 0 $$

$$ \int_{\Omega} (q - \chi(g)g) \delta q d\Omega = \delta W(\underline{\phi}, \delta q) = 0 $$

En introduisant la loi de comportement :

$$ \int_{\Omega} (\underline{S}^r + \chi(g)g \underline{I}) : \delta \dot{\underline{E}} d\Omega = \int_{\Omega} \underline{B} \cdot \delta \underline{v} d\Omega + \int_{\Gamma_T} \underline{T}_d \cdot \delta \underline{v} dA = \delta W(\underline{\phi}, \delta \underline{v}) = 0 $$

$$ \int_{\Omega} (q - \chi(g)g) \delta q d\Omega = \delta W(\underline{\phi}, \delta q) = 0 $$

Nous avons introduit trois types de non-linéarités. Une non-linéarité géométrique, une matérielle, et une dernière avec la fonction $\chi$. Les deux premières se linéarisent classiquement avec un processus itératif type Newton. La dernière, par le processus otératif introduit appelé algorithme des contraintes actives[2, 4, 6].

## 7 Algorithme et linéarisation

Procédons tout d'abord à la mise en place de l'algorithme des contraintes actives. On commence par inititaliser la fonction $\chi$. La fonction est alors mise à jour après convergence de l'algorithme de Newton. L'algorithme complet est détaillée Figure 1.
Nous précisons maintenant le détail de la linéarisation effectuée dans l’algorithme de Newton. Nous partons de la formulation faible exprimée plus haut, à la différence près que nous connaissions entièrement la fonction $\chi$, grâce à la mise en place de l’algorithme des contraintes actives.

Si nous notons le travail virtuel interne :

$$ 
\begin{align*}
\delta W_{\text{int}}(\phi, \delta\underline{v}) &= \int_{\Omega} \underline{S} : \delta \dot{\underline{E}} \\
D\delta W_{\text{int}}(\phi, \delta\underline{v})[\underline{u}] &= \int_{\Omega} D\underline{S}[\underline{u}] : \delta \dot{\underline{E}} d\Omega + \int_{\Omega} \underline{S} : D\delta \dot{\underline{E}}[\underline{u}] d\Omega \\
D\underline{S}[\underline{u}] &= \mathbb{C} : D\underline{E}[\underline{u}] \\
\mathbb{C} &= 2 \frac{\partial \underline{S}^r}{\partial \underline{C}} + 2 \frac{\partial(\chi(g)g\underline{I})}{\partial \underline{C}}
\end{align*}
$$

Maintenant que $\chi$ est complètement connu, nous pouvons écrire :

$$ 
\begin{align*}
\mathbb{C} &= \lambda \underline{C}^{-1} \otimes \underline{C}^{-1} + 2(\mu - \lambda \ln(J))J + 2\chi \rho_n \mathbb{I} \\
D\delta W_{\text{int}}(\phi, \delta\underline{v})[\underline{u}] &= \int_{\Omega} \delta \dot{\underline{E}} : (\lambda \underline{C}^{-1} \otimes \underline{C}^{-1} + 2(\mu - \lambda \ln(J))J + 2\chi \rho_n \mathbb{I}) : D\underline{E}[\underline{u}] \\
&\quad + \int_{\Omega} (\underline{S}^r + \chi(q + \rho_n(tr(\underline{C}) - I_m)\underline{I})) : ((\underline{\nabla}_0 \underline{u}^T)\underline{\nabla}_0 \delta\underline{v}) d\Omega \\
D\delta W_{\text{int}}(\phi, \delta\underline{v})[\underline{\Delta} p] &= \int_{\Omega} D\underline{S}[\underline{\Delta} p] : \delta \dot{\underline{E}} d\Omega = \int_{\Omega} \chi \Delta p \operatorname{tr}(\delta \dot{\underline{E}}) d\Omega
\end{align*}
$$

Le terme $D\delta W(\phi, \delta q)[\Delta p]$ se linéarise de la manière suivante :

$$ D\delta W(\phi, \delta q)[\Delta p] = \int_{\Omega} (1 - \chi) \Delta p \delta q d\Omega $$

Enfin, pour le terme $D\delta W(\phi, \delta q)[\underline{u}]$, nous obtenons :

$$ D\delta W(\phi, \delta q)[\underline{u}] = -2\rho_n \int_{\Omega} \chi \operatorname{tr}(D\underline{E}[\underline{u}]) \delta q d\Omega $$

# 8 Discrétisation et résultats numériques

Du fait du caractère mixte de la formulation établie, les espaces discrets ne peuvent être choisis au hasard. Les éléments utilisés doivent non seulement vérifier la condition inf-sup (ou LBB [3]), mais aussi donner de bons résultats quant à la convergence de l’algorithme des contraintes actives. Or, il nous est très vite apparu que certains éléments, bien vérifiant la LBB, n’assuraient pas la convergence de l’algorithme. La raison principale de ce souci de convergence vient de la nature irrégulière du champ $\chi$. Après analyse de plusieurs éléments ainsi qu’une étude sur un cas linéaire monodimensionnel, nous avons extrait un élément satisfaisant toutes les conditions requises (l’élément retenu est l’élément $P_2^+/P_{-1}$, interpolation quadratique enrichie avec une fonction bulle pour le déplacement et interpration linéaire discontinue pour la pression).

Le cas test présenté Figure 2 est le cas d'une plaque rectangulaire trouée en déformations planes sollicitée en traction sur sa face supérieure. Le matériau est de type néo-hookéen avec $\lambda = 0.5$ et $\mu = 0.4225$. Le maillage est non-structuré et compte 600 éléments. Deux simulations sont effectuées : l'une en prenant en compte la contrainte matérielle ($I_m = 10$), et l'autre sans aucune contrainte. La Figure 3 montre une cartographie mettant en relief $tr(\underline{C})$ sur la plaque dans le cas contraint et non contraint. On remarque bien, dans le cas contraint, une limitation de $tr(\underline{C})$ à la valeur $I_m$.
<figure><img src="image_4.png" /><figcaption>Figure 3 – Visualisation de $I_1 = tr(\underline{C})$. Cas non contraint (gauche) et contraint (droite)</figcaption></figure>

Références

[1] E. Arruda and M.C. Boyce. A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials. *Rubber chemistry and technology*, Vol.41, page389–page412, 1993.

[2] H. BenDhia and I. Vautier and M. Zarroug. Problèmes de contact frottant en grandes transformations : du continu au discret. *Revue Europénne des Éléments Finis*, Vol.9, page243–page261, 2000.

[3] F. Brezzi and M. Fortin. Mixed and Hybrid Finite Element Methods. Springer, 1991.

[4] G. Dumont. Algorithme des contraintes actives et contact unilatéral sans frottement. *Revue Europénne des Éléments Finis*, Vol.4, page55–page73, 1995.

[5] J.L. Ecricksen. Constitutive theory for some constrained elastic crystals. *International journal of solids and structures*, Vol.22, page951–page964, 1986.

[6] Samuel Geniaut. *Approche X-FEM pour la fissuration sous contact des structures industrielles*, Ecole Centrale de Nantes, 2006.

[7] A.N. Gent. A new constitutive relation for rubber. *Rubber. Chem. Technol.*, Vol.69, page59–page61, 1996.

[8] C.O. Horgan and G. Saccomandi. Constitutive modeling of rubber-like and biological materials with limiting chain extensibility. *Math. Mech. Solids*, Vol.7, page353–page371, 2002.

[9] Grégory Legrain. *Extension de l'approche X-FEM aux grandes transformations pour la fissuration des milieux hyperélastiques*, Ecole Centrale de Nantes, 2006.

[10] Y. Marco and I. Moro and N. Ait Hocine and V. Lesaux and P. Charrier. Cristallisation induite sous sollicitations cycliques dans les élastomères. 18ème congrès de mécanique, 2007.

[11] Nicolas Saintier. *Fatigue multiaxiale dans un élastomère de type NR chargé : mécanismes d'endommagement et critère local d'amorçage de fissure*, Ecole des Mines de Paris, 2001.

[12] S. Trabelsi and P.- Albouy and J. Rault. Stress-induced crystallization properties of natural and synthetic cis-polyisoprene. *Rubber chemistry and technology*, 2004.

[13] C. Truesdell and W. Noll. The non-linear field theories of mechanics. Springer-Verlag, 1965.