# Une amélioration dans X-FEM du raccord entre enrichissement et éléments finis classiques

Elie Chahine, Patrick Laborde, Yves Renard

► To cite this version:

Elie Chahine, Patrick Laborde, Yves Renard. Une amélioration dans X-FEM du raccord entre enrichissement et éléments finis classiques. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01412116

HAL Id: hal-01412116

https://hal.science/hal-01412116

Submitted on 8 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 Attribution 4.0 International License
# Une amélioration dans XFEM du raccord entre enrichissement et éléments finis classiques

E. Chahine¹, P. Laborde², Y. Renard³

¹ Laboratory for Nuclear Materials, Nuclear Energy and Safety Research Department, Paul Scherrer Institute OVGA/14, CH-5232 Villigen PSI, Switzerland
elie.chahine@psi.ch

² Institut de Mathématiques de Toulouse, UMR CNRS 5215, Univ. Paul Sabatier (Univ. de Toulouse)
118 route de Narbonne, 31062 Toulouse Cedex 4, France
patrick.laborde@math.ups-tlse.fr

³ Institut Camille Jordan, CNRS UMR 5208, INSA de Lyon, Université de Lyon, 20 rue Albert Einstein, 69621 Villeurbanne Cedex, France
yves.renard@insa-lyon.fr

**Résumé** — On s'intéresse au calcul des déformations d'un corps élastique fissuré par la méthode des éléments finis étendue. On sait que la qualité de l'approximation est perturbée par la couche des éléments finis de transition entre zone d'enrichissement autour du fond de fissure et le reste du maillage. On propose de remplacer la couche de transition par une condition intégrale, de type "raccord mortar", à l'interface entre ces deux zones. Nous montrerons comment cette variante de XFEM permet d'améliorer sensiblement la précision.

**Mots clés** — Fissure, XFEM, approximation non-conforme.

## 1 Introduction

On s'intéresse à la simulation numérique des déformations d'un corps élastique linéaire iso-trope fissuré en utilisant un maillage éléments finis du domaine non fissuré. Sur ce maillage in-dépendant de la fissure, la méthode XFEM consiste à "enrichir" la base éléments finis à l'aide de fonctions qui prennent en compte la connaissance a priori que l'on a de la solution [10]. Par exemple, on dispose de la forme asymptotique du déplacement en pointe de fissure, ce qui nous fournit un enrichissement "singulier" que l'on peut appliquer à tous les nœuds d'une surface don-née autour du fond de fissure.

On sait que la précision se dégrade sur la couche des éléments situés à la transition entre la surface d'enrichissement singulier et le reste du domaine (voir [6] dans un autre contexte et [8]). De plus, dans le but de diminuer le nombre de degrés de liberté supplémentaires et d'améliorer le conditionnement, il est intéressant d'utiliser une fonction cut-off pour "globaliser" l'enrichis-sement sur la surface fixée [4]. On constate alors que la difficulté demeure au passage entre zone enrichie et zone non enrichie par les fonctions singulières. Il s'agit de remédier à cet inconvénient en s'affranchissant de la zone de transition.
<figure><img src="image_2.png" /><figcaption>Figure 1 – Décomposition du domaine fissuré.</figcaption></figure>

## 2 Une méthode d'approximation non-conforme

On considère un domaine fissuré $\Omega$ du plan, $\Gamma_C$ représentant la fissure. Pour simplifier, $\Omega$ est supposé polygonal et $\Gamma_C$ rectiligne. On note par $\Omega_1$ et $\Omega_2$ une partition du domaine non fissuré $\bar{\Omega}$ de telle sorte que la pointe de fissure $x^*$ appartienne à $\Omega_2$ et que la frontière de $\Omega_2$, notée $\partial\Omega_2$, soit polygonale (Figure 1).

Soit une triangulation $\mathcal{T}_h$ du domaine non fissuré $\bar{\Omega}$ telle que l'interface $\partial\Omega_2$ entre $\Omega_1$ et $\Omega_2$ coïncide avec des arêtes d'éléments de $\mathcal{T}_h$. Le maillage éléments finis est donc indépendant de la fissure, $h$ désignant le paramètre de discrétisation. Sur $\mathcal{T}_h$ on formule une méthode d'éléments finis $P_1$ dont les fonctions de base scalaires seront notées $\{\varphi_i\}_{i\in I}$. La fonction de forme $\varphi_i$ est associée au noeud $x_i$. On appelle $I_H$ l'ensemble des numéros de noeuds "enrichis" par la fonction de type Heaviside :

$$H(x) = \begin{cases} +1 & \text{if } (x-x^*) \cdot n \ge 0 \\ -1 & \text{ailleurs,} \end{cases} \tag{1}$$

où $n$ est un vecteur normal à $\Gamma_C$. Il s'agit des noeuds dont la fonction de forme associée a son support entièrement coupé par la fissure (voir [10]). Par ailleurs, $I(\Omega_k)$ (resp. $I_H(\Omega_k)$) est le sous-ensemble d'indices $i \in I$ (resp. $i \in I_H$) tels que $x_i \in \Omega_k$, $k \in \{1, 2\}$.

On désigne par $\mathcal{V}_1^h$ l'espace de fonctions éléments finis $P_1$ vectorielles associé à la triangulation induite sur $\Omega_1$ et enrichi par la fonction saut $H$ :

$$\mathcal{V}_1^h = \left\{ v_1^h : v_1^h = \sum_{i \in I(\Omega_1)} a_i \varphi_i + \sum_{i \in I_H(\Omega_1)} b_i H \varphi_i; a_i, b_i \in \mathbb{R}^2 \right\}. \tag{2}$$

L'espace $\mathcal{V}_2^h$ est constitué des fonctions éléments finis vectorielles sur $\Omega_2$ enrichi, non seulement à l'aide de la fonction $H$, mais aussi par les fonctions "singulières" de fond de fissure :

$$\mathcal{V}_2^h = \left\{ v_2^h : v_2^h = \sum_{i \in I(\Omega_2)} a_i \varphi_i + \sum_{i \in I_H(\Omega_2)} b_i H \varphi_i + \sum_{j=1}^4 c_j F_j; a_i, b_i, c_j \in \mathbb{R}^2 \right\}, \tag{3}$$

où les fonctions d'enrichissement $F_j$ en coordonnées polaires de fond de fissure sont données par (voir [10]) :

$$\{F_j(x)\}_{1 \le j \le 4} = \left\{ \sqrt{r \sin \frac{\theta}{2}}, \sqrt{r \cos \frac{\theta}{2}}, \sqrt{r \sin \frac{\theta}{2} \sin \theta}, \sqrt{r \cos \frac{\theta}{2} \sin \theta} \right\}. \tag{4}$$

Dans ce qui suit, le champ de déplacement solution défini dans le milieu fissuré $\Omega$ sera cherché dans l'espace $\mathcal{V}^h$ constitué des fonctions $v^h$ définies sur $\Omega$ telles que $(k \in \{1, 2\})$ :

$$v^h = v_k^h \quad \text{dans } \Omega_k, \quad \text{où } v_k^h \in \mathcal{V}_k^h. \tag{5}$$
L'espace d'approximation global $V^h$ peut être identifié à l'espace produit $V_1^h \times V_2^h$.

La surface d'enrichissement sur laquelle agissent les fonctions singulières de fond de fissure est le sous-domaine $\Omega_2$. Soulignons que les fonctions d'enrichissement introduites (dans la définition de $V_2^h$) opèrent "globalement" sur $\Omega_2$ et, de ce fait, ne rajoutent que 8 degrés de liberté supplémentaires à la méthode d'éléments finis.

Naturellement, la méthode d'approximation envisagée sera non-conforme, le déplacement solution dans $V^h$ n'étant pas continu à travers l'interface entre $\Omega_1$ et $\Omega_2$. Dans la section suivante, nous définissons une condition de raccord pour récupérer d'une certaine façon la propriété de continuité de la solution à travers l'interface.

### 3 Une condition de raccord inégal

Le saut du champ des déplacements à l'interface s'écrit :

$$|u| = (u_2 - u_1)|_\Gamma \qquad (6)$$

en appelant $\Gamma = \partial\Omega_2 \cap \Omega$ l'interface "coupée" par la fissure. La condition de continuité du déplacement à la traversée de $\Gamma$ peut s'écrire sous forme intégrale :

$$b(u, \mu) = 0 \quad \text{pour tout } \mu \text{ dans } \mathcal{W}, \qquad (7)$$

où $\mathcal{W}$ est un espace convenable de multiplicateurs définis sur $\Gamma$ et

$$b(v, \mu) = - \int_{\Gamma} \mu \cdot [v] \, d\Gamma. \qquad (8)$$

Afin de transcrire cette condition de continuité dans l'espace des fonctions d'approximation discontinues $V^h$, on introduit les multiplicateurs approchés suivants. La triangulation $\mathcal{T}_h$ choisie en Section 2 définit une subdivision $S_h$ en segments élémentaires de l'interface non coupée $\bar{\Gamma} = \partial\Omega_2$. L'espace $W^h$ des multiplicateurs discrets s'écrit

$$W^h = \{ \mu^h \in C^0(\bar{\Gamma})^2 : \mu_i^h|_{S_i} \in P_i, \forall S_i \in S_h, i \in \{1, 2\} \}. \qquad (9)$$

Remarquons que cette définition ne tient pas compte de la discontinuité à travers la fissure. On considère alors le problème hybride (voir [2]) :

$$\begin{align}
& \text{trouver } u^h \in V^h = V_1^h \times V_2^h, \lambda^h \in W^h \quad \text{tel que} \nonumber \\
& a(u^h, v^h) + b(v^h, \lambda^h) = L(v^h) \quad \text{pour tout } v^h \in V^h, \tag{10a} \\
& b(u^h, \mu^h) = 0 \quad \text{pour tout } \mu^h \in W^h. \tag{10b}
\end{align}$$

En notant $u = (u_1, u_2)$, la forme bilinéaire "brisée" $a$ est définie par :

$$a(u, v) = \sum_{k=1}^{2} a_k(u_k, v_k) = \sum_{k=1}^{2} \int_{\Omega_k} D\varepsilon(u_k) : \varepsilon(v_k) \, dx, \qquad (11)$$

où $D$ est l'opérateur des rigidités élastiques. Par ailleurs,

$$L(v) = \sum_{k=1}^{2} \int_{\Omega_k} g \cdot v_k \, dx + \int_{\Gamma_N} f \cdot v_1 \, d\Gamma, \qquad (12)$$

avec les notations $f$ et $g$ pour les forces connues s'exerçant respectivement sur la portion $\Gamma_F$ de la frontière et à l'intérieur du corps $\Omega$. On a supposé, pour simplifier, que le contact entre les lèvres de la fissure n'est pas activé au cours des déformations.
<figure><img src="image_3.png" /><figcaption>Figure 2 – Maillage éléments finis du domaine non fissuré et surface d'enrichissement singulier.</figcaption></figure>

L'espace d'approximation $V^h$ ne peut être identifié à un sous-espace de l'espace de Sobolev $H^1(\Omega, \mathbb{R}^2)$ classique en élasticité. Dans cette méthode non-conforme, la continuité à travers l'interface $\Gamma$ du champ de déplacement approché dans $V^h$ n'est imposée que sous une forme faible discrète ; c'est la condition de raccord intégral (10b). La condition de raccord ponctuelle :

$$[u] = 0 \text{ en chaque noeud sur } \Gamma, \qquad (13)$$

fournit un autre type d'approximation de la condition de continuité à la transition.

On peut établir une d'estimation d'erreur optimale quand le pas du maillage $h$ tend vers zéro. Plus précisément, on montre que l'erreur en norme de l'énergie est majorée par $Ch$ où $C$ est une constante indépendante de $h$. La preuve s'appuie sur une condition inf-sup qui traduit une compatibilité convenable entre l'espace des déplacements $V^h$ et l'espace des multiplicateurs $W^h$. La condition de raccord discrète envisagée permet d'assurer la propriété de coercivité de la forme bilinéaire brisée $a$ nécessaire pour appliquer les résultats mathématiques généraux sur les problèmes hybrides (voir [7]).

## 4 Simulations numériques

Les tests suivants sont réalisés sur le domaine fissuré :

$$\Omega = [-0.5; 0.5[\times] - 0.5; 0.5[\setminus \Gamma_C],$$

où $\Gamma_C$ désigne la fissure définie par :

$$\Gamma_C = [-0.5; 0] \times \{0\}.$$

La condition aux limites de déplacements imposés est fournie par le mode d'ouverture $u_I$ exact (voir [9]).

On choisit une méthode d'éléments finis $P_I$ sur un maillage structuré du milieu non fissuré $\bar{\Omega}$, comme indiqué sur la Figure 2. La surface d'enrichissement singulier (Section 2) est constituée de tous les éléments inclus dans la boule $B(x^*, 0.2)$ centrée sur la pointe de fissure $x^*$. La méthode XFEM a été implémentée en utilisant GETFEM++, une librairie éléments finis orientée objet en C++ développée dans notre équipe (voir [11]).
<figure><img src="image_4.png" /><figcaption>Figure 3 – Surfaces de niveau de la contrainte de Von Mises pour un problème en mode I résolu en P₁ par la méthode XFEM avec raccord intégral.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4 – Erreur relative en norme de l'énergie par rapport au nombre d'éléments dans chaque direction pour un problème en mode I (échelles logarithmiques).</figcaption></figure>

La Figure 3 présente les surfaces de niveaux de la contrainte de Von Mises sur la structure déformée. Autour du fond de fissure et dans la zone de transition au voisinage de la frontière de la surface d'enrichissement singulier, on observe que le champ reste très régulier, à comparer avec la méthode d'approximation conforme (voir [8]).

La Figure 4 compare les courbes d'erreur de convergence, quand le pas de maillage *h* tend vers zéro, obtenues par la méthode XFEM avec surface d'enrichissement fixe conforme (voir [1], [8]), par XFEM avec fonction cut-off (voir [4]) et par XFEM avec raccord intégral présentée ici. Ces différentes erreurs sont calculées en norme de l'énergie en fonction du nombre *N* = 1/h d'éléments dans chaque direction).

Les résultats numériques sont analogues pour un maillage non structuré ou bien pour un problème de fissuration en mode mixte.

Enfin, la Figure 5 présente les valeurs du premier facteur d'intensité de contrainte *K*<sub>I</sub> calculées par la méthode de l'intégrale *J* (voir [10] pour le calcul de l'intégrale de contour). La solution exacte est une combinaison linéaire donnée entre les deux modes *u*<sub>I</sub> et *u*<sub>II</sub>, égale à *P(x)* + 3*u*<sub>I</sub> + 5*u*<sub>II</sub>, où *P(x)* est une solution régulière fixée sur le domaine non fissuré. Les courbes de convergence pour la méthode d'éléments finis classique (FEM), pour XFEM avec fonction cut-off function et pour XFEM avec surface d'enrichissement fixe sont régulières, mais surestiment la valeur exact du FIC. La méthode XFEM avec raccord intégral s'approche mieux de cette valeur, même si elle
<figure><img src="image_6.png" /><figcaption>Figure 5 – *Erreur relative dans le calcul du facteur d'intensité de contrainte K*₁.</figcaption></figure>

présente quelques oscillations.

En conclusion, nous avons illustré comment le méthode XFEM avec raccord intégral introduite augmente la précision numérique en améliorant la qualité de l'approximation autour du fond de fissure et sur la zone de transition entre le surface d'enrichissement et le reste du domaine.

Ce travail correspond à un chapitre de la thèse [3] et à l'article [5].

Références

[1] E. Béchet, H. Minnebo, N. Moës and B. Burgardt. Improved implementation and robustness study of the X-FEM for stress analysis around cracks. *Int. J. Numer. Meth. Engng.*, 64 :1033–1056, 2005.

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

[3] E. Chahine. Etude mathématique et numérique de métodes d'éléments finis étendues pour le calcul en domaines fissurés. *Thèse, INSA Toulouse, 2008*.

[4] E. Chahine, P. Laborde, and Y. Renard. Crack tip enrichment in the XFEM method using a cut-off function. *Int. J. Numer. Meth. Engng.*, 64 :354–381, 2005.

[5] E. Chahine, P. Laborde, and Y. Renard. Integral matching XFEM : a non-conformal eXtended Finite Element approach. *A paraitre*.

[6] J. Chessa, H. Wang and T. Belytschko. On the construction of blending elements for local partition of unity enriched finite elements. *Int. J. Numer. Meth. Engng.*, 57 : 1015–1038, 2003.

[7] A. Ern, J.L. Guermond. *Eléments finis : théorie, applications, mise en oeuvre*, Springer, 2002.

[8] P. Laborde, Y. Renard, J. Pommier, and M. Salaün. High order extended finite element method for cracked domains. *Int. J. Numer. Meth. Engng.*, 64 :354–381, 2005.

[9] J. Lemaitre, J.L. Chaboche. *Mechanics of Solid Materials*, Cambridge University Press, 1994.

[10] N. Moës, J. Dolbow, and T. Belytschko. A finite element method for crack growth without remeshing. *Int. J. Numer. Meth. Engng.*, 46 :131–150, 1999.

[11] Y. Renard, J. Pommier. Getfem++. *An open source generic C++ library for finite element methods*, http://home.gna.org/getfem.