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

HAL
open science

# Implémentation robuste de lois cohésives non régulières en formulation X-FEM

Guilhem Ferte, Patrick Massin, Nicolas Moes

► To cite this version:

Guilhem Ferte, Patrick Massin, Nicolas Moes. Implémentation robuste de lois cohésives non régulières en formulation X-FEM. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722099

HAL Id: hal-01722099

https://hal.science/hal-01722099

Submitted on 2 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.

Public Domain
# Implémentation robuste de lois cohésives non régulières en formulation X-FEM

Guilhem FERTE <fnref n="1" />,<fnref n="2" />*, Patrick MASSIN <fnref n="1" />, Nicolas MOES <fnref n="2" />

**Résumé** — Une méthodologie est proposée pour introduire des lois cohésives non régulières, c'est à dire dont la dérivée est discontinue d'un régime de fonctionnement à l'autre, dans le formalisme des éléments finis étendus (X-FEM). Nous mettons en évidence des problèmes de discrétisation de la variable interne, sur laquelle repose la détermination du régime de fonctionnement, avec une formulation standard. Nous adaptons à X-FEM deux formulations à plusieurs champs de la littérature afin d'y remédier tout en gardant une implémentation simple: une pour les éléments linéaires, l'autre pour les quadratiques. La méthode est illustrée sur des calculs quasi-statiques faisant appel à des techniques de pilotage du chargement.

**Mots clés** — zone cohésive, X-FEM, pilotage du chargement

## 1. Introduction

Afin de réduire l'effort de remaillage lors de calcul de propagation de fissures ou d'évolution d'interface, une extension de la méthode des éléments finis, dénommée X-FEM, a été introduite par [1]. Celle-ci permet l'introduction de surfaces de discontinuités dans un modèle sans qu'il soit besoin de les mailler, par un enrichissement local de la base de fonctions tests avec des fonctions discontinues. Pour avoir une description physique des champs au voisinage de la pointe, Barenblatt [2] a proposé d'introduire des forces cohésives s'exerçant mutuellement entre les lèvres de la fissure, forces obéissant à une loi adoucissante à partir d'un certain seuil afin de rendre compte de la fissuration.

<figure><img src="image_2.png" /><figcaption>Fig. 1: modèle de zone cohésive.</figcaption></figure>

La fissure est alors décrite par (voir fig.1):

- une zone fermée,

- une zone partiellement ouverte, mais dont les bords sont soumis à des forces cohésives, de longueur caractéristique $l_c$ et d'ouverture caractéristique $\alpha_c$,

- une zone ouverte.

<footnote n="1">1 Laboratoire de mécanique des structures industrielles durables
EDF R&D, 1 avenue du général de Gaulle, 92141 Clamart, France.</footnote>

<footnote n="2">2 GeM, Ecole Centrale de Nantes,
1 rue de la Noë, BP 92101, 44321 Nantes, France</footnote>

<footnote>* Auteur correspondant</footnote>
Par cet article, nous mettons en évidence les différents problèmes numériques se posant lorsque l'on implante les éléments cohésifs en formulation X-FEM. L'objectif est d'apporter des réponses simples à ces différentes problématiques, pour aboutir à une formulation robuste tant pour les éléments finis linéaires que pour les quadratiques.

## 2. Une formulation explicite

### 2.1. Formulation générale d'un problème de fissuration cohésive

Nous considérons un solide occupant une position $\Omega$ dans sa configuration de référence, et traversé par une interface $\Gamma$, qui représente un trajet potentiel de fissuration, et sépare $\Omega$ en deux parties $\Omega_i, i \in \{1, 2\}$. Le bord $\partial\Omega$ se décompose en deux parties $\Gamma_u$ et $\Gamma_g$, sur lesquelles sont imposées un déplacement $\mathbf{u}_0$ et une force surfacique $\mathbf{g}$, respectivement (voir fig.2). On note $\mathbf{u}_i$ le champ de déplacement dans $\Omega_i$, de sorte que son saut est calculé par $[[\mathbf{u}](\mathbf{x})=\mathbf{u}_2(\mathbf{x})-\mathbf{u}_1(\mathbf{x})]$. Si on appelle $\mathbf{n}$ la normale à $\Gamma$ sortante à $\Omega_2$, la composante normale de ce saut est négative en décollement et positive en interpénétration. Enfin, si $\Gamma_i$ est la lèvre de la fissure sur le côté $i$, on appelle $\mathbf{t}_c$ le vecteur contrainte cohésif exercé par $\Gamma_1$ sur $\Gamma_2$ (fig.2).

<figure><img src="image_3.png" /><figcaption>Fig. 2: notations du problème.</figcaption></figure>

Le vecteur $\mathbf{t}_c$ est relié à $[[\mathbf{u}]]$ par la loi cohésive, qui est supposée présenter les caractéristiques suivantes (voir fig.3):

* l'adhésion initiale est parfaite ou quasi-parfaite (la pente initiale est très forte),

* il y a des changements de régime brutaux, se traduisant par une discontinuité de la dérivée.

Pour les éléments linéaires, nous considérerons par exemple la loi cohésive pénalisée suivante [3]:

$$ \mathbf{t}_c = \frac{\sigma_c}{\alpha} \exp\left(-\frac{\sigma_c}{G_c} \alpha\right) [[\mathbf{u}]] + \kappa \langle [[\mathbf{u}]_n]_+ \mathbf{n} \rangle \quad (1) $$

où $\alpha$ est la variable interne définie par $\alpha|_{t} = \max\{\alpha_0, \max_{s \in [0,t]} [[u(s)]]_{\text{eq}}\}$ avec $[[u]]_{\text{eq}} = \sqrt{\langle [[u]]_n \rangle^2 + [[\mathbf{u}]]_\tau^2}$. Avec ces définitions, une formulation faible explicite du problème est:

$$ \int_{\Omega} \sigma(\mathbf{u}):\epsilon(\mathbf{u}^*) d\Omega - \int_{\Gamma} \mathbf{t}_c([[\mathbf{u}]]) \cdot [[\mathbf{u}^*]] d\Gamma = \int_{\Gamma_g} g \mathbf{u}^* d\Gamma_g \quad (2) $$
<figure><img src="image_4.png" /><figcaption>Fig. 3: loi cohésive régularisée</figcaption></figure>

## 2.2. Discrétisation du déplacement et inconvénient d'une formulation explicite

La méthode des éléments finis étendus est couramment utilisée de pair avec une fonction distance signée $\phi$ dite « level-set », dont l'iso-zéro définit la position de l'interface. L'interpolation des éléments finis classiques peut alors être enrichie avec une fonction Heaviside afin de décrire le saut de déplacement à travers l'interface [1] :

$$ \mathbf{u}_h(\mathbf{x}) = \sum_{i \in N_h(x)} \mathbf{a}_i N_i(\mathbf{x}) + \sum_{j \in N_h(x) \cap K_h} \mathbf{b}_j N_j(\mathbf{x}) H(\phi(\mathbf{x})) \quad (3) $$

où $N_h$ est la fonction de forme du noeud $i$, linéaire ou quadratique; $N_h(x)$ est l'ensemble des nœuds dont le support contient $\mathbf{x}$ ; $K_h=N_h(\Gamma)$ est l'ensemble des nœuds dont le support est strictement intersecté par $\Gamma$ (voir fig.6); et $H$ est la fonction Heaviside utilisée pour représenter la discontinuité $H(\phi)=\begin{cases} -1 & \text{si } \phi < 0 \\ 1 & \text{si } \phi > 0 \end{cases}$.

Avec une telle discrétisation et la formulation faible (2), nous considérons un test de déchirement d'une bande (peeling test), comme montré sur la fig.4. Afin de tester la robustesse de la formulation, nous choisissons les paramètres matériau de sorte à ce que la longueur cohésive, estimée par la formule $l_c \sim \frac{E}{1-v^2} \frac{G_c}{\sigma_c^2}$, s'étale sur au moins une dizaine d'éléments, et un paramètre de pénalisation grand devant l'inverse du saut de déplacement caractéristique $\alpha_c = \frac{G_c}{\sigma_c}$ de la loi cohésive: $\kappa=10^6 \frac{G_c}{\sigma_c}$.

<figure><img src="image_5.png" /><figcaption>Fig. 4: peeling test.</figcaption></figure>
Si nous traçons alors la contrainte cohésive explicite déduite du saut de déplacement, nous observons des oscillations très importantes. Le problème, connu sous le nom de *surface locking* est dû au fait que la contrainte à l'interface est discrétisée de façon trop riche.

<figure><img src="image_6.png" /><figcaption>Fig. 5: oscillations de la contrainte cohésive avec une formulation explicite.</figcaption></figure>

# 3. Formulations avec champs duaux

## 3.1. Formulation avec champs duaux pour les éléments finis linéaires

Pour y remédier, il est nécessaire de restreindre l'espace $L_h$ sur lequel la contrainte cohésive est définie, afin que la condition inf-sup suivante soit vérifiée:

$$ \exists \alpha > 0 \text{ indépendant de } h, \forall \mu_h \in L_h, \sup_{v_h \in V_h} \frac{\int_{\Gamma} \mu_h[\|v_h\|] d\Gamma}{h^{1/2} \|v_h\|_{1,\Omega}} \geq \alpha \|\mu_h\|_{0,\Gamma} \quad (4) $$

où $\|\cdot\|_{1,\Omega}$ désigne la norme de $H^1(\Omega)$ et $\|\cdot\|_{0,\Gamma}$ celle de $L^2(\Gamma)$.

<figure><img src="image_7.png" /><figcaption>Fig. 6: notations pour la discrétisation des champs primaux et duaux.</figcaption></figure>

Béchet et Moës ([4] et [5]) abordèrent le problème dans le cadre du contact, et proposèrent l'adjonction d'un champ « dual » défini sur $L_h$ pour résoudre le problème. Les degrés de libertés duaux sont alors portés par $M_h$, qui sont les noeuds sommets de $K_h$. Un champ sur les éléments intersectés (colorés sur la fig.6) est déduit par interpolation linéaire, dont on prend la trace sur $\Gamma$ pour définir le champ dual (fig.6). Cependant, l'espace étant trop riche en l'état, on le restreint
en imposant des relations d'égalité entre les DDL originaux portés par $M_h$: un DDL réellement indépendant $I$ sera alors partagé par un groupe de plusieurs noeuds, dont la somme des fonctions de forme définit une fonction de forme « étendue » $\Psi_I$ pour le DDL $I$ (voir fig.7).

En un mot, le principe de l'algorithme de restriction [5] est de prendre un ensemble le plus grand possible d'arêtes indépendantes, c'est à dire non connectées entre elles, dont chacune définira un DDL dual, et d'attribuer aux noeuds n'appartenant pas à une arête indépendante l'un des DDL des arêtes voisines (fig.7). L'espace restreint ainsi construit vérifie (4), comme démontré dans [5].

<figure><img src="image_8.png" /><figcaption>Fig. 7: algorithme de restriction pour le champ dual.</figcaption></figure>

Les auteurs écrivent alors le problème consistant dans le cas du contact, qui revient à résoudre un problème du type [4]:

$$ \left[ \begin{matrix} K_u^1 & 0 & 0 \\ 0 & K_u^2 & A^T \\ 0 & A & \frac{1}{\kappa} C \end{matrix} \right] \left[ \begin{matrix} \{a\} \\ \{b\} \\ \{\lambda\} \end{matrix} \right] = \left[ \begin{matrix} \{L_{\text{méca}}^1\} \\ \{L_{\text{méca}}^2\} \\ 0 \end{matrix} \right] \quad (5) $$

En éliminant $\{\lambda\}$ de (5), on s'aperçoit que ceci revient à utiliser un déplacement consistant $C^{-1}A\{b\}$ pour écrire la pénalisation à l'interface. En écrivant les expressions matricielles, on retrouve la définition d'un projecteur: le saut de déplacement consistant est en fait la projection $L_2$ du saut primal sur $L_h$. Avec le même paradigme pour notre loi cohésive, nous introduisons un saut de déplacement consistant $w \in L_h$ qui devra vérifier $w+[[u]] \in L_h^{\perp}$ (nous inversons les conventions de signe par rapport à la fig.2). Nous souhaitons alors définir la contrainte cohésive non comme fonction de $[[u]]$ mais comme fonction de $w$. Compte-tenu de ces considérations, nous écrivons le Lagrangien du problème comme:

$$ L(u, \lambda, w) = \int_{\Omega} \Phi_{\epsilon l}(\epsilon(u)) d\Omega - \int_{\Gamma} g \cdot u d\Gamma_g - \int_{\Gamma} \lambda \cdot (w + [[u]]) d\Gamma + \int_{\Gamma} \Pi(w) d\Gamma \quad (6) $$

où $\Pi(w)$ est la densité d'énergie cohésive dissipée dans l'ouverture $\Pi(w) = \int_{0}^{w} t_c(w') d w'$. On déduit alors les équations du problème en écrivant les conditions d'optimalité du lagrangien:

$$ \forall u^* \in V_h, \int_{\Omega} \sigma(u):\epsilon(u^*) d\Omega + \int_{\Gamma_g} g \cdot u^* d\Gamma_g - \int_{\Gamma} \lambda \cdot [[u^*]] d\Gamma = 0 \quad (7) $$

$$ \forall \lambda^* \in L_h, -\int_{\Gamma} ([[u]] + w) \cdot \lambda^* d\Gamma = 0 \quad (8) $$
$$ \forall \mathbf{w}^* \in L_h, \int_{\Gamma} (\mathbf{t}_c(\mathbf{w}) - \lambda) \cdot \mathbf{w}^* d\Gamma = 0 \qquad (9) $$

Nous retrouvons bien la projection $L_2$ avec (8), tandis que (9) permet d'écrire la contrainte cohésive en fonction du saut consistant. Notons qu'une telle formulation n'est pas sans analogies avec celle proposée par Gravouil dans le cadre du couplage multi-échelles [6]. En reprenant le test précédent, nous observons alors une nette amélioration de la contrainte cohésive à l'interface, à présent complètement régulière.

<figure><img src="image_9.png" /><figcaption>Fig. 8: contrainte cohésive régulière pour la formulation avec champs duaux.</figcaption></figure>

## 3.2. Discrétisation de la variable interne

Cependant, le calcul peine toujours à converger sur certains pas de temps. Ceci s'explique par la discrétisation de la variable interne aux points de Gauss -adoptée en première intention en tant que méthode la plus simple et courante pour les lois non linéaires. Il y a beaucoup plus d'indicateurs de régime que de DDL d'interface (voir fig.9), donc ces indicateurs oscillent d'une itération de la méthode de Newton à une autre, ce qui ralentit la convergence.

<figure><img src="image_10.png" /><figcaption>Fig 9: discrétisations possibles et loi de complémentarité.</figcaption></figure>

Si nous considérons, pour nous aider à raisonner, un problème de complémentarité linéaire (contact unilatéral), nous voyons que la condition de complémentarité est violée pour les points de collocation
(les points de Gauss) dans la zone de transition (fig.9). En revanche, une discrétisation à raison d'un indicateur par degré de liberté de $L_h$ permet de respecter cette condition (fig.9).

Ainsi, il n'est plus question d'évaluer une contrainte cohésive au point de Gauss dans les expressions des matrices et vecteurs élémentaires, mais nous décomposons à la place $t_c(x_g) = \sum_j \psi_j(x_g) t_c(\alpha_j, w_j)$. Le problème converge plus vite avec cette formulation, puisque le temps CPU pour notre calcul passe à 700 secondes, contre 1695 secondes pour la formulation explicite initiale.

## 3.2. Formulation pour les éléments finis quadratiques

Reprenons la formulation à trois champs de (6), que nous augmentons d'un terme $\int_{\Gamma} \frac{r}{2} (w + \|u\|^2) d\Gamma$, sans présupposer d'une discrétisation particulière:

$$ L(u, \lambda, w) = \int_{\Omega} \Phi_{\dot{\epsilon}l}(\epsilon(u)) d\Omega - \int_{\Gamma} g \cdot u d\Gamma_g - \int_{\Gamma} \lambda \cdot (w + \|u\|^2) d\Gamma + \int_{\Gamma} \Pi(w) d\Gamma + \int_{\Gamma} \frac{r}{2} (w + \|u\|^2) d\Gamma \quad (10) $$

Les conditions d'une discrétisation consistante des conditions d'optimalité de ce Lagrangien telles qu'exposées dans [8] sont les suivantes:

- le respect de la condition inf-sup déjà évoquée,

- la surjectivité de l'application trace de $V_h$ dans l'espace de $w$.

Avec un déplacement linéaire, ceci impose un enrichissement de $V_h$ avec des fonctions bulles élaborées, mais avec un déplacement quadratique, il existe un set de discrétisation consistant simple pour le problème (voir [7] et [8]):

- le déplacement est quadratique,

- le champ $w$ est discrétisé par collocation aux points de Gauss de l'interface,

- le champ $w$ est linéaire sur l'interface.

Avec ce set de discrétisation, comme démontré dans [7], la condition d'optimalité déduite de la dérivation par rapport à $w$ s'écrit comme une loi de comportement locale: on parle de condensation statique.

La généralisation de ces résultats à X-FEM nécessite de trouver une construction de champ P1 sur l'interface qui respecte la condition inf-sup avec le déplacement P2. On peut réutiliser l'espace défini précédemment, mais celui-ci devient trop pauvre et génère des ordres de convergence sous-optimaux: on perd ainsi le gain de précision des éléments P2. Nous avons montré que l'on peut trouver une construction optimale de la façon suivante:

- on impose une relation d'égalité sur les arêtes intersectées près d'un noeud (voir fig.11),

- on liste les composantes du graphe des arêtes intersectées restantes. Pour chaque composante pour laquelle aucun noeud n'est encore soumis à une relation, on impose une relation d'égalité sur une des arête de la composante, choisie au hasard (voir fig.11).

Lors d'une rupture brutale, la réponse de la structure peut présenter des instabilités. Ainsi, la grandeur de contrôle peut diminuer au cours du chargement. On parle alors de *snap-back* ou retour en arrière. Pour réaliser un calcul, il devient nécessaire passer en dynamique ou d'utiliser une méthode de pilotage du chargement. Nous reprenons le pilotage par prédiction élastique proposé dans [9] et étendu à XFEM dans un acte de la session 2011 de ce congrès [10]. Pour des structures critiques, le calcul parvient alors à son terme tant avec les formulations à éléments quadratiques que linéaires.
<figure><img src="image_11.png" /><figcaption>Fig. 10: espace dual stable et optimal pour un déplacement quadratique.</figcaption></figure>

# 5. Conclusion et perspectives

La poursuite de notre travail consiste à réaliser des calculs de propagation cohésive sur trajet inconnu. Pour ce faire, à partir du résultat d'un calcul cohésif à un instant i, nous cherchons à détecter le front de décohésion. A partir de ce dernier, nous déterminerons une direction de propagation à partir d'un critère de la littérature, soit déduit du champ local de contraintes au voisinage de la pointe, soit déduit du calcul de facteurs d'intensité de contraintes équivalent, soit déduit d'une minimisation énergétique globale. Nous actualiserons alors la surface potentielle de fissuration en exploitant la littérature disponible sur les actualisations de level-sets, avant d'effectuer le calcul pour l'instant i+1.

# Références

[1] N. Moës, J. Dolbow, T. Belytschko. A finite element method for crack growth without remeshing. *International Journal for Numerical methods in engineering*, **46**: 135-150, 1999.

[2] J.I. Barenblatt. The mathematical theory of equilibrium cracks in brittle fracture. *Adv. Appl. Mech.*, **7**: 55-125, 1962.

[3] G. Meschke, P. Dumstorff. Energy-based modelling of cohesive and cohesionless cracks via X-FEM. *Computational Methods in Applied Mechanics Engineering*, **196**: 2338-2357, 2007.

[4] N. Moës, E. Béchet, M. Tourbier. Imposing Dirichlet boundary conditions in the extended finite element method. *International Journal for numerical methods in engineering*, **67**: 1641-1669, 2006.

[5] E. Béchet, M. Moës, B. Wohlmuth. A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method. *Int. J. Num. Meth. Eng.*, **78**: 931-954, 2009.

[6] E. Pierres, M.C. Baietto, A. Gravouil. A two-scale extended finite element method for modelling 3D crack growth with interfacial contact. *Comput. Meth. Appl. Mech. Eng.*, **199**: 1165-1178, 2010.

[7] E. Lorentz. A mixed interface finite element for cohesive zone models. *Computer methods in applied mechanics and engineering*, **198**: 317-320, 2008.

[8] D. Doyen, A. Ern, S. Piperno. A three-field augmented Lagrangian formulation of unilateral contact problems with cohesive forces. *ESIAM, Math. Model. Numer. Anal.*, **44**: 323-346, 2010.

[9] E. Lorentz, P. Badel. A new path-following constraint for strain-softening finite element simulations. *International Journal for numerical methods in engineering*, **60**: 499-526, 2004.

[10] P. Massin, G. Ferté, A. Caron, N. Moës. Pilotage du chargement en formulation X-FEM: application aux lois cohésives, *CSMA 2011*, Giens.