# Étude de la stabilité de modèles d'endommagement à gradient : cylindre sous torsion

Kim Pham, Jean-Jacques Marigo, Corrado Maurini

► To cite this version:

Kim Pham, Jean-Jacques Marigo, Corrado Maurini. Étude de la stabilité de modèles d'endommagement à gradient : cylindre sous torsion. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01418169

HAL Id: hal-01418169

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

Submitted on 16 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
# Etude de la stabilité de modèles d'endommagement à gradient: cylindre sous torsion

K. Pham¹, J-J. Marigo¹, C. Maurini¹

¹ Institut Jean Le Rond d'Alembert
Université Paris 6, 4 Place Jussieu, Case 162, 75252 Paris Cedex 05
pham@lmm.jussieu.fr, marigo@lmm.jussieu.fr, maurini@lmm.jussieu.fr

**Résumé** — Nous étudions un cylindre mince que l'on soumet à une torsion croissante. Celui-ci est composé d'un bimatériaux adoucissant-incassable. Le matériau adoucissant est modélisé à l'aide d'un modèle d'endommagement à gradient. En suivant une approche variationnelle, on s'intéresse alors à l'évolution et à la stabilité de la structure au cours du chargement. Des expériences numériques sont alors comparées aux résultats théoriques.

**Mots clés** — endommagement, régularisation, modèle à gradient, stabilité, matériau adoucissant

## Introduction

L'usage de matériaux fragiles adoucissants comme le béton ou la céramique est courante dans de nombreuses constructions. Cependant, leur étude reste encore difficile pour de multiples raisons : d'un point de vue théorique, le problème d'évolution est mal posé [1] et l'utilisation de modèles standards d'endommagement ne peut conduire qu'à l'apparition de solutions instables. Ceci se traduit d'un point de vue numérique par de la "localisation" i.e. une sérieuse dépendance au maillage du champs d'endommagement. L'introduction de modèles non-locaux, vus comme des régularisés des modèles standards, a pour but de pallier à ces problèmes [2,4]. Ceux-ci peuvent être vus comme approximation de la rupture fragile [3]. Dans la formulation adoptée dans [2], la régularisation intervient directement dans la forme postulée de l'énergie de déformation W

$$W(\epsilon, \alpha, \nabla\alpha) = \left( \frac{1}{2} E \ell^2 \nabla\alpha \cdot \nabla\alpha + \frac{1}{2} A(\alpha) \epsilon \cdot \epsilon + w(\alpha) \right) \quad (1)$$

où $\alpha$ et $\epsilon$ désignent l'endommagement et la déformation tandis que $E$, $A(\alpha)$ et $w$ représentent respectivement le module de Young du matériau sain, la loi de comportement et l'énergie volumique dissipée dans le matériau élastique endommageable au cours d'un processus d'endommagement homogène ($\nabla\alpha = 0$). Le paramètre $\ell$ joue le rôle d'une longueur interne au modèle et mesure l'influence de la régularisation. On se propose d'étudier dans le cadre de cette classe de modèles d'endommagement un cylindre mince sous torsion. La première partie est consacrée à la mise en place du problème. Dans la deuxième partie, nous étudions l'évolution et la stabilité du problème au cours du temps. Enfin dans une dernière partie, nous simulons le problème numériquement et comparons les résultats issus de la simulation à la théorie.

## 1 Mise en place du problème

On considère un cylindre mince $\Omega$ de hauteur L, de rayon extérieur $R$, et d'épaisseur $e$. On suppose que $e$ est petite devant les autres dimensions de la structure, ainsi que devant la longueur
<figure><img src="image_2.png" /><figcaption>Figure 1 – Cylindre mince sous torsion</figcaption></figure>

interne $\ell$ (1) au modèle non-local. Le cylindre est constitué de deux matériaux $M1$ (élastique adoucissant) et $M2$ (élastique incassable) qui occupent respectivement les domaines $\Omega_1$ (d'épaisseur $\eta_1 e$) et $\Omega_2$ (d'épaisseur $\eta_2 e$) avec $\eta_1 + \eta_2 = 1$. On suppose en outre que les deux matériaux sont élastiques isotropes, leurs lois de comportement pouvant alors s'écrire localement

$$A_1(\alpha)\varepsilon = \lambda_1(\alpha)\mathrm{Tr}(\varepsilon)\mathbf{1} + 2\mu_1(\alpha)\varepsilon \quad A_2\varepsilon = \lambda_2\mathrm{Tr}(\varepsilon)\mathbf{1} + 2\mu_2\varepsilon$$

où $\lambda_1(\alpha)$ et $\mu_1(\alpha)$ sont des fonctions qui dépendent localement de l'endommagement. On appelle $w_1(\alpha)$ l'énergie volumique dissipée dans le matériau 1 au cours d'un processus d'endommagement homogène au sein du matériau $M1$. Enfin la non localité du modèle d'endommagement se fait en rajoutant dans l'énergie de déformation un terme en gradient d'endommagement. L'énergie volumique de déformation s'écrit alors

$$W_1(\varepsilon, \alpha, \nabla\alpha) = \frac{1}{2} E_1 \ell^2 \nabla\alpha \cdot \nabla\alpha + \frac{1}{2} A_1(\alpha) \varepsilon \cdot \varepsilon + w_1(\alpha) \quad \text{dans } M1 \quad (2)$$

$$W_2(\varepsilon, \alpha, \nabla\alpha) = \frac{1}{2} A_2 \varepsilon \cdot \varepsilon \quad \text{dans } M2 \quad (3)$$

où $\ell$ désigne la longueur interne au modèle adoucissant régularisé. On applique alors une torsion croissante au cylindre via les conditions de chargement en déplacement suivantes (cf (Fig.1)) :

$$- \text{en } z = 0 : u_{\theta} = u_{r} = 0 \text{ et } (\sigma \cdot e_z)_z = 0$$

$$- \text{en } z = L : u_r = 0, u_{\theta} = rt \text{ et } (\sigma \cdot e_z)_z = 0$$

$$- \text{en } r = R : \text{bord libre, } \sigma \cdot e_r = 0$$

où $t$ est un paramètre de chargement croissant de 0 à $+\infty$.

## 2 Problème d'évolution

On définit les variables d'endommagement et de déplacement les espaces d'admissibilité suivants

$$\mathcal{D} = \{ \beta | \beta \in H^1(\Omega_1, \mathbb{R}) \text{ avec } \beta \ge 0 \text{ p.p.} \} \quad (4)$$

$$C_t = \{ u | u \in H^1(\Omega, \mathbb{R}^3) \text{ avec } u_{\theta}(r, \theta, L) = tr \} \quad (5)$$

Pour un couple $(u, \alpha)$ pris dans $\mathcal{D} \times C_t$, l'énergie potentielle du système s'écrit donc

$$\mathcal{E}(u, \alpha) = \int_{\Omega_1} W_1(\varepsilon(u), \alpha, \nabla\alpha) d\Omega + \int_{\Omega_2} W_2(\varepsilon(u), \alpha, \nabla\alpha) d\Omega \quad (6)$$

L'évolution du problème se formule de manière variationelle à partir de l'énergie potentielle en écrivant la condition d'optimalité du premier ordre. Il faut cependant veiller à tenir compte ici de
la condition d'irreversibilité $\dot{\alpha}_t \ge 0$. L'équation d'évolution au temps $t$ s'écrit alors formellement

$$ \text{Trouver } (\alpha_t, u_t) \in \mathcal{D} \times \mathcal{C}_t \text{ tel que} \\ d\mathcal{E}(u_t, \alpha_t)(v - \dot{u}_t, \beta - \dot{\alpha}_t) \ge 0 \quad \forall (v, \beta) \in \dot{\mathcal{C}}_t \times \mathcal{D} \qquad (7) $$

En utilisant les arguments classiques des calculs de variation (voir [2] pour plus de détails), (7) peut se séparer en deux problèmes implicitement liés. Le premier est un problème classique d'élasticité à endommagement fixé et s'obtient en faisant $\beta = \dot{\alpha}_t$ dans (7)

$$ \int_{\Omega_1} A_1(\alpha_t) \varepsilon(u_t) \varepsilon(v - \dot{u}_t) d\Omega + \int_{\Omega_2} A_2 \varepsilon(u_t) \varepsilon(v - \dot{u}_t) d\Omega = 0 \quad \forall v \in \dot{\mathcal{C}}_t \qquad (8) $$

En faisant $v = \dot{u}_t$ dans (7), on obtient un problème en endommagement, qui après intégration par parties s'écrit localement :

* Dans $\Omega_1$ :

- Irreversibilité : $\dot{\alpha}_t \ge 0$

- Critère d'endommagement : $\frac{1}{2}A'(\alpha_t)\varepsilon(u_t)^2 + w'(\alpha_t) - E_0\ell^2\Delta\alpha_t \ge 0$

- Condition de Kuhn-Tucker : $\dot{\alpha}_t \cdot (\frac{1}{2}A'(\alpha_t)\varepsilon(u_t)^2 + w'(\alpha_t) - E_0\ell^2\Delta\alpha_t) = 0$

* Sur $\partial\Omega_1$ : Conditions Limites de type Neumann $\frac{\partial\alpha_t}{\partial n} = 0$

Avec la montée en chargement, plusieurs phases d'évolution vont se succeder. Nous nous proposons dans la suite de faire une étude analytique du problème en étudiant chacune de ces phases en terme d'évolution et de stabilité.

## 2.1 Phase élastique

Durant la phase élastique, l'endommagement est nul dans le cylindre. On connait alors explícitement la forme du déplacement associé pour chaque pas de chargement. Celui-ci s'écrit sous la forme

$$ u_{\text{elas}} = \frac{trz}{L} e_{\theta} \qquad (9) $$

En réinjectant (9) dans le critère d'endommagement () et en utilisant la nullité de l'endommagement dans la structure (pour la phase élastique), on voit que le critère d'endommagement n'est pas satisfait (inégalité stricte) tant que $t < t_{\text{elas}}$ où $t_{\text{elas}}$ est donné par

$$ t_{\text{elas}} = \sqrt{\frac{-2w'_{1}(0)L^2}{\mu'_{1}(0)R^2}} \qquad (10) $$

Pour $t = t_{\text{elas}}$, le critère en endommagement () est nul et l'endommagement apparaît pour la première fois dans la structure. D'autre part, d'après (10), pour $w'_{1}(0) = 0$, on remarque cette phase élastique n'a pas lieu (i.e. $t_{\text{elas}} = 0$), le cylindre pouvant s'endommager dès le premier instant de mise en chargement. Par exemple, si $w_1$ est une loi proportionnelle à $\alpha^2$, on trouve que $t_{\text{elas}} = 0$ et le cylindre s'endommagera donc dès le premier instant de mise en chargement. En revanche, pour des énergies volumiques dissipées $w_1$ proportionnelle à $\alpha$, le temps $t_{\text{elas}}$ sera non nul et on observera donc bien une phase élastique.

## 2.2 Phase endommagée pour $e \ll 1$

### 2.2.1 Evolution de la branche axisymétrique

Une fois que le critère d'endommagement est atteint, l'endommagement va pouvoir croître dans la structure. Dans cette partie nous nous interessons à des endommagements solution du problème d'évolution de type axisymétrique i.e. qui ne dépendent que de la variable $r$. En étudiant l'énergie dans le cas où l'on se restreint à des endommagements axisymétriques, on peut montrer
qu'il existe à chaque instant un seul et unique champ d'endommagement axisymétrique $\alpha_e \in \mathcal{D}$ qui soit solution du problème d'évolution. Lorsque l'épaisseur faible, $\alpha_e$ tend vers une valeur homogène dans toute la barre qui ne dépend que tu temps et qui est donnée par

$$ \alpha_e = \alpha_0 + o(1) \\ \text{avec } -w'_1(\alpha_0)/\mu'_1(\alpha_0) = t^2 R^2 / 2L^2 \quad (11) $$

En prenant comme loi d'endommagement

$$ \mu_1 = \mu_1^0 (1 - \alpha)^2 \quad \text{et} \quad w_1(\alpha) = \alpha, \quad \alpha \in [0, 1] \quad (12) $$

on trouve

$$ \alpha_0 = 1 - \frac{1}{t^2}, \quad t \ge 1 \quad (13) $$

### 2.2.2 Stabilité de la branche endommagée axisymétrique

Suivant le temps de chargement, la taille de la structure ou les paramètres matériaux, il est possible que la solution axisymétrique ne soit plus un chemin stable d'évolution. D'un point de vue énergétique, cela signifierait que cet état ne serait plus un minimum local pour l'énergie potentielle mais seulement un point selle. Le système cherchera alors à évoluer naturellement vers un autre état d'énergie moindre. C'est à cette question que nous allons tenter de répondre dans cette partie en étudiant l'énergie potentielle autour du point $(\alpha_e, u_e)$, où $u_e$ représente le champ de déplacement optimal associé à $\alpha_e$ (8). L'état $(\alpha_e, u_e)$ sera dit stable pour un chargement $t$ donné si et seulement s'il existe un voisinage détats admissibles dans lequel tout autre accessible a une énergie potentielle supérieure. Pour ce faire, nous allons chercher des bifurcations particulières en déplacement de la forme $u_e + hv(r, \theta)e_z$ (entraînant donc un gauchissement du cylindre) et en endommagement de la forme $\alpha_e + h\beta(r, \theta)$ avec $h$ et $\beta$ positifs. La positivité de $\beta$ est fondamentale car on ne peut aller tester pour la stabilité des états voisins qui violeraient la condition d'irréversibilité. En développant l'énergie potentielle au deuxième ordre autour du point critique axisymétrique, nous obtenons

$$ \mathcal{E}(u_e + hv e_z, \alpha_e + h\beta) = \mathcal{E}(\alpha_e, u_e) + h\mathcal{E}'(\alpha_e, u_e)(\beta, v) + \frac{h^2}{2}\mathcal{E}''(\alpha_e, u_e)(\beta, v) + o(h^2) \quad (14) $$

Durant la phase élastique, le critère en endommagement étant strictement positif, l'énergie première est strictement positive pour des $\beta$ non nuls. La phase élastique est donc stable. En revanche, lorsque l'endommagement apparaît dans la structure ($t \ge t_{\text{elas}}$), l'énergie première ne suffit plus pour conclure. Il faut alors s'interesser à l'énergie seconde. Celle-ci peut se réeccrire après calculs sous la forme condensée suivante

$$ \mathcal{E}''(u_e, \alpha_e)(\beta, v) = C(v, \beta) \cdot (\mathcal{R}_e^t(v, \beta) - 1) \quad \text{avec } C(v, \beta) > 0 \quad (15) $$

La stabilité de l'état $(u_e, \alpha_e)$ équivaut donc à étudier la première valeur propre du quotient de Rayleigh

$$ \alpha_e \text{ est stable ssi } \lambda_e^t = \inf_{\beta, v} \mathcal{R}_e^t > 1 \quad (16) $$

La valeur asymptotique de la première valeur propre du quotient de Rayleigh peut se calculer de façon explicite. Pour le choix (12) de loi d'endommagement, celle-ci s'exprime en fonction de 3 paramètres adimensionnés $\eta_2\mu_2/\eta_1\mu_1^0$, $t/t_{\text{elas}}$ et $E_1\ell/(2R\sigma_1)$

$$ \lambda_0^t = \lim_{e \to 0} \lambda_e^t = \frac{4}{3} \cdot \min \left\{ 1, 1 - \frac{1}{1 + \frac{\eta_2 \mu_2}{\eta_1 \mu_1^0} \cdot \left(\frac{t}{t_{\text{elas}}}\right)^4} + \left( \frac{E_1 \ell / (2R\sigma_1)}{1 + \frac{\eta_2 \mu_2}{\eta_1 \mu_1^0} \cdot \left(\frac{t}{t_{\text{elas}}}\right)^4} \cdot \frac{t_{\text{elas}}}{t} \right)^{2/3} \right\} \quad (17) $$
La présence d'un matériau incassable (i.e. $\frac{\eta_2\mu_2}{\eta_1\mu_1^0} > 0$) rend donc nécessairement la solution homogène stable pour des temps assez longs puisque $\lambda_0' \to \frac{4}{3}$ lorsque $t \to \infty$. En revanche, si celui-ci est absent (i.e. $\frac{\eta_2\mu_2}{\eta_1\mu_1^0} = 0$), on trouve que $\lambda_0' \to 0$ lorsque $t \to \infty$. Ceci confirme le résultat attendu lorsque qu'il n'y a aucune d'élasticité résiduelle dans la structure : la solution homogène ne peut être stable puisque la structure localise l'endommagement sur une zone caractéristique fixée de taille $\ell$. Nous traçons sur (Fig.2) à présent $\lambda_0'$ pour $\ell = 0.01$ et $\frac{\eta_2\mu_2}{\eta_1\mu_1^0} = 1$ (les autres constantes étant prises égales à 1). Pour ce choix de fraction, on constate que l'état est instable juste après la fin du

<figure><img src="image_3.png" /><figcaption>Figure 2 – $\ell = 0.01$ et $\frac{\eta_2\mu_2}{\eta_1\mu_1^0} = 1$</figcaption></figure>

régime élastique. Il va y avoir alors bifurcation dans le chemin d'évolution avec une localisation de l'endommagement. Après un certain laps de temps, on remarque sur (Fig.2) que l'état homogène redevient stable. Néanmoins, cela ne signifie pas que le système va naturellement revenir vers cet état uniforme d'endommagement. En effet, celui-ci peut se retrouver fort éloigné et il est possible qu'il soit nécessaire de franchir de hautes barrières d'énergie avant de l'atteindre de nouveau. La partie numérique va nous permettre de mieux comprendre les bifurcations et les phénomènes non linéaires d'instabilité dans la structure.

## 3 Etude numérique comparative

Le modèle d'endommagement à gradient (1) sur lequel se base notre étude nous offre l'avantage de pouvoir travailler avec des éléments finis classiques. En effet, l'endommagement se calcule de la même manière que le déplacement, c'est à dire à chaque noeud du maillage et non pas après post-traitement aux points de Gauss comme dans une formulation locale du problème d'endommagement. Dans le cas du modèle adoucissant (1), l'énergie n'étant plus convexe par rapport au couple déplacement-endommagement, on n'est plus assuré que le minimum global de l'énergie soit l'unique point critique. De ce fait, la résolution de la formulation variationelle, équivalente à la condition d'optimalité d'ordre 1, nous assure seulement d'avoir trouvé un point critique qui peut être au choix un point selle ou un minimum local. Pour obtenir à chaque pas de temps un point critique de l'énergie, on utilise un algorithme de minimisation alternée qui s'appuie sur les propriétés mécaniques du système. En effet, on sait que l'énergie totale est convexe par rapport à chacune des variables isolément (l'autre étant fixée), on est donc assuré de converger vers un point critique à l'issue de ce processus. D'autre part, bien que le problème de départ sur le cylindre mince soit tridimensionnel, l'hypothèse de faible épaisseur nous permet de travailler sur la forme asymptotique de l'énergie définit sur $[0, 2\pi]$ qui s'écrit après homogénéisation dans l'épaisseur

$$ \int_{0}^{2\pi} \frac{1}{2} \eta_1 \frac{E_1 \ell^2}{R^2} \alpha_{,\theta}^2 + \frac{1}{2} (\eta_1 \mu_1(\alpha) + \eta_2 \mu_2) \left( \frac{tR}{L} + \frac{v_{,\theta}}{R} \right)^2 + w(\alpha) d\theta \quad (18) $$
En prenant la loi quadratique (12) pour l'endommagement, on peut alors utiliser des éléments finis 1D pour la simulation. L'algortihme consiste alors à minimiser alternativement dans chaque direction, l'autre direction étant fixée. On prend le choix suivant de coefficients,

$$E_1 = 1, \frac{\eta_2 \mu_2}{\eta_1 \mu_1^0} = 1, \ell = 5.10^{-3}, L = 1, R = 1, \sigma_1 = 1 \text{ d'où } t_{\text{elas}} = 1 \quad (19)$$

(Fig.3) montre les différentes étapes de la formation d'endommagement au cours du temps (Fig.3-

<figure><img src="image_4.png" /><figcaption>Figure 3 – $\alpha$ sur $[0, 2\pi]$ pour $t_A = 1, t_B = 1.1, t_C = 1.3, t_D = 1.5, t_E = 1.6, t_F = 1.7$</figcaption></figure>

A à Fig.3-F) qu'on a reporté sur (Fig.2). La courbe verte continue représente la courbe homogène qui est tracée d'après les valeurs analytiques (13). La courbe bleue pointillée est la courbe obtenue numériquement par l'algorithme de minimisation alternée. Tant qu'on est dans la phase élastique (Fig.3-A), il n'y a pas d'endommagement dans la structure. Une fois atteint le temps $t_{\text{elas}}$, il y a comme prévu (Fig.3-B) bifurcation dans le chemin d'évolution. Il se crée un "pustule" qui va croître au cours du temps en longueur à endommagement maximal fixé. (Fig.3-B) et (Fig.3-C) sont des snapshots de la simulation pris dans le domaine d'instabilité de la solution axisymétrique (cf. (Fig.2)). De ce fait, la solution localisée bleue ne peut bifurquer vers la solution homogène verte car cette dernière est instable. En revanche (Fig.3-D) a été prise dans la zone de stabilité de la solution axisymétrique. Cependant la solution localisée ne peut toujours pas "sauter" vers cette dernière pour des raisons d'irreversiblité. En effet comme le montre (Fig.3-D), la solution localisée bleue présente des zones plus endommagées que la solution homogène verte. Bien que la solution homogène soit donc stable de nouveau, il est donc impossible de bifurquer à nouveau vers celle-ci. En revanche, (Fig.3-E) et (Fig.3-F) montrent que une fois que les deux solutions se sont rejointes, on ne constate plus de bifurcations dans l'évolution du problème. Ceci vient donc confirmer l'analyse de stabilitée de la solution homogène.

## Références

[1] A. Benallal, R. Billardon, and G. Geymonat. Bifurcation and localization in rate independent materials. In Q. Nguyen, editor, C.S.I.M Lecture Notes on Bifurcation and Stability of Dissipative Systems. Springer-Verlag.

[2] A. Benallal and J.-J. Marigo. Bifurcation and stability issues in gradient theories with softening. Modelling Simul. Mater. Sci. Eng., 15:283-295, 2007.

[3] B. Bourdin, G. Francfort, and J.-J. Marigo. Numerical experiments in revisited brittle fracture. J. Mech. Phys. Solids, 48:797-826, 2000.

[4] E. Lorentz and S. Andrieux. Analysis of non-local models through energetic formulations. Int. J. Solids Struct., 40:2905-2936, 2003.