# Méthodes X-FEM d'ordre supérieur en mécanique de la rupture

Patrick Laborde, Julien Pommier, Yves Renard, Michel Salaün

► To cite this version:

Patrick Laborde, Julien Pommier, Yves Renard, Michel Salaün. Méthodes X-FEM d'ordre supérieur en mécanique de la rupture. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01812961

HAL Id: hal-01812961

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

Submitted on 12 Jun 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.

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

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Méthodes XFEM d'ordre supérieur en mécanique de la rupture

Patrick Laborde * — Julien Pommier **

Yves Renard ** — Michel Salaün ***

* UPS - Toulouse 3, Laboratoire MIP (CNRS UMR 5640)
118 route de Narbonne, F-31062 Toulouse cedex 4
laborde@mip.ups-tlse.fr

** INSA Toulouse, Laboratoire MIP (CNRS UMR 5640)
Complexe scientifique de Rangueil, F-31077 Toulouse
{pommier, renard}@insa-toulouse.fr

*** ENSICA, 1 pl. Emile Blouin, F-31056 Toulouse cedex 5
msalaun@ensica.fr

**RÉSUMÉ.** Le calcul par éléments finis en mécanique de la rupture bénéficie du progrès considérable apporté par la méthode XFEM, ces cinq dernières années. Sur le problème non-régulier de la propagation de fissures, la précision de cette nouvelle méthode est sensiblement meilleure qu'avec une méthode d'éléments finis classique de même degré. On peut toutefois observer que la vitesse de convergence (quand le pas de maillage diminue) reste inchangée. Nous proposons diverses améliorations de XFEM : enrichissement d'une zone élargie autour du fond de fissure (fixée indépendamment du pas), condition de liaison entre les degrés de liberté sur cette zone enrichie, traitement de l'interface entre la zone d'enrichissement et le reste du domaine. Notre approche est validée sur un cas test où la précision optimale attendue est atteinte.

**ABSTRACT.** The XFEM method is a significant advance in computational fracture mechanics. For the nonsmooth problem of the fracture propagation in an elastic body, the accuracy of this new method is better than with a classical finite element method of the same degree. But it may be observed that the rate of convergence (when the mesh parameter decreases) remains unchanged. We present some improvements of the original XFEM method: enrichment on a enlarged zone around the crack tip (fixed independently of the mesh parameter), bonding condition between the singular enrichment degrees of freedom, treatment of the transition layer at the boundary of the enriched zone. Numerical simulations show that the expected optimal rate of convergence is achieved.

**MOTS-CLÉS :** Fissure, éléments finis, XFEM, vitesse de convergence optimale, raccord ponctuel

**KEYWORDS:** Fracture, finite elements, XFEM, optimal rate of convergence, pointwise matching.
# 1. Introduction

Dans beaucoup d'applications industrielles en calcul de structures (aéronautique, matériaux composites, assemblages collés, structures en béton..), les ingénieurs doivent de nos jours simuler de plus en plus finement la propagation de fissures. On connaît bien les difficultés de la prise en compte de problèmes de mécanique de la rupture dans une modélisation éléments finis : nécessité d'adapter le maillage à la ligne de fracture, raffinement drastique indispensable du maillage en fond de fissure, remaillage pour suivre la progression incrémentale de la fissure. En permettant le découplage du maillage éléments finis et de la géométrie de la fissure, la méthode XFEM (eXten-ded Finite Element Method) permet de remédier à ces inconvénients. Cette méthode introduite par [MOE 99] a montré son efficacité et a suscité un certain nombre de développements en 2D ou en 3D, ces cinq dernières années.

Une première étude a montré que la vitesse de convergence de la méthode XFEM (quand le paramètre de maillage $h$ diminue) n'était pas optimale [STA 03]. L'objet de la présente communication est de proposer des alternatives pour surmonter cette difficulté.

Dans la suite, nous commençons par rappeler le principe d'enrichissement de la base éléments finis qu'introduit la méthode XFEM classique (section 2). Puis nous présentons une première variante, dans l'objectif d'augmenter la vitesse de convergence de l'approximation. La zone d'enrichissement autour de la pointe de fissure est élargie, en fait on la fixe indépendamment du paramètre de maillage $h$ qui diminue (section 3).

On introduit ensuite une condition de liaison entre les degrés de liberté relatifs à cet enrichissement afin notamment d'améliorer le conditionnement (section 4). La zone de transition, c'est-à-dire la couche d'éléments entre la zone d'enrichissement fixée et le reste du domaine, contribue à dégrader la qualité de l'approximation. On propose alors une dernière variante dans le cadre d'une méthode non conforme : supprimer la transition dans l'approximation et la remplacer par une condition de continuité du déplacement aux noeuds à l'interface (section 5). De plus amples détails sur ce travail peuvent être trouvés dans [LAB 04].

# 2. Position du problème

On considère dans le plan un milieu fissuré $\Omega$, la ligne de fracture étant notée $\Gamma_C$ et la pointe de la fissure $x_0$. Soit $u$ le champ de déplacement solution du problème d'élasticité linéarisée isotrope correspondant. L'expression asymptotique des composantes de $u$ au voisinage du fond de fissure est une combinaison linéaire de fonctions singulières $F_1, ..., F_4$. Celles-ci sont définies dans un voisinage de $x_0$ en coordonnées polaires $(r, \theta)$ par :

$$F_1 = \sqrt{r} \sin \frac{\theta}{2}, \quad F_2 = \sqrt{r} \cos \frac{\theta}{2}, \quad F_3 = \sqrt{r} \sin \frac{\theta}{2} \cos \theta, \quad F_4 = \sqrt{r} \cos \frac{\theta}{2} \cos \theta.$$
On effectue une triangulation du domaine non fissuré $\bar{\Omega}$ paramétrée par $h > 0$. On note $\varphi_1, ..., \varphi_N$ la base vectorielle des fonctions polynômiales par morceaux de degré $k$ définies sur le maillage. Soit par ailleurs $H$ une fonction de type Heaviside définie dans un voisinage de la fissure par la valeur $+1$ d'un côté de $\Gamma_C$ et $-1$ de l'autre.

L'espace d'approximation fourni par la méthode XFEM classique s'écrit alors :

$$V_h = \left\{ v_h = \sum_{i=1}^{N} a_i \varphi_i + \sum_{i \in I_H} b_i H \psi_i + \sum_{i \in I_F} \sum_{j=1}^{4} c_{ij} F_j \psi_i \text{ avec } a_i, b_i, c_{ij} \in \mathbb{R} \right\}$$

Dans cette définition, $\psi_1, \psi_2, ...$ désigne la base éléments finis vectorielle des fonctions polynômiales par morceaux de degré 1. Les fonctions $\psi_1, \psi_2, ...$ constituent localement une partition de l'unité. Les ensembles d'indices $I_H$ et $I_F$ numérotent les degrés de liberté enrichis par la fonction saut $H$ ou les fonctions singulières $F_j$ respectivement. On met dans la base d'approximation la fonction $H\psi_i$ si le support de $\psi_i$ est entièrement coupé par $\Gamma_C$ et l'on rajoute $F_j\psi_i$ seulement pour une fonction de forme $\psi_i$ associée à un sommet du triangle contenant $x_0$ [MOE 99].

La Figure 1 fournit une comparaison entre la méthode XFEM et la méthode des éléments finis classique pour une approximation polynomiale de degré $k = 1, 2, 3$. Compte-tenu de la singularité de la solution, on sait que la méthode des éléments finis admet une vitesse de convergence en $\sqrt{h}$, et cela quel que soit le degré $k$. Les simulations ont été réalisées sur un problème test : fissure rectiligne sur un carré en mode I d'ouverture. Même si la comparaison est à l'avantage de XFEM (niveau d'erreur plus faible à maillage fixé), on peut noter que la vitesse de convergence reste toujours en $\sqrt{h}$ (lorsque le paramètre de discrétisation $h$ tend vers 0). Cette observation, qui avait déjà été faite par [STA 03], montre que l'on ne retrouve pas l'optimalité d'une méthode d'éléments finis classique, et cela bien que la singularité soit prise en compte dans l'approximation.

### 3. Elargissement de la zone d'enrichissement

Une première cause de non-optimalité provient de ce que l'enrichissement par les fonctions singulières ne concerne que le triangle contenant le fond de fissure. Dans l'étude de vitesse de convergence, la taille de cet élément tend vers zéro comme $h$ et l'influence de l'enrichissement va diminuant. Il semble donc naturel d'élargir la zone de l'enrichissement à une région autour du fond de fissure indépendante de la taille du maillage. Ainsi, on choisit d'enrichir avec les $F_j$ tous les noeuds contenus dans un disque $B(x_0, R)$ centré au fond de fissure et de rayon $R$ indépendant de $h$. Dans l'expression de l'approximation XFEM précédente, l'ensemble d'indices $I_F$ est donc remplacé par un ensemble $I_F(R)$ convenable. Pour le terme d'enrichissement en $H$ dans la définition de l'espace d'approximation, on choisit dorénavant la partition de l'unité fournie par les fonctions de forme $\varphi_i$ afin de rester cohérent avec l'approximation éléments finis de degré $k$ dans le domaine.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Vitesse de convergence de la méthode des éléments finis classique (FEM) et de la méthode XFEM pour une approximation polynômiale de degré 1, 2 ou 3. Echelles log., abscisses : $1/h$, ordonnées : erreur relative en norme de l'énergie.</figcaption></figure>

Reprenons le même problème test qu'en Figure 1 avec (i) une partition de l'unité de degré $k$ pour l'enrichissement par la fonction saut $H$, (ii) une zone d'enrichissement $B(x_0, R)$ dont le rayon est de l'ordre du dixième du côté du domaine. L'erreur en norme de l'énergie est alors en $h^\alpha$, où la vitesse de convergence $\alpha$ vaut 0.9, 1.8 ou 2.6 selon le degré $k = 1, 2$ ou $3$ (au lieu de $\alpha = 0.5$ auparavant). On obtient ainsi une précision presque optimale. L'idée d'élargir la zone d'enrichissement singulier a été développée indépendamment par [BEC 04] et [LAB 04].

Le nombre de degrés de liberté va naturellement augmenter quelque peu, mais surtout on peut observer une dégradation du conditionnement quand on monte en degré $k$ ou lorsque $h$ diminue. Ce phénomène peut notamment s'expliquer par le fait que les fonctions d'enrichissement $F_j\varphi_i$ ne sont pas linéairement indépendantes au niveau élémentaire. Ainsi, pour une partition de l'unité $P$, on a les deux relations de dépendance linéaire :

$$ \phi_2(F_1 - F_4) + \phi_1 F_3 = 0, \quad \phi_2(F_3 - F_2) + \phi_1 F_4 = 0 $$

$$ \text{avec } \phi_1(x) = x_1, \quad \phi_2(x) = x_2 \text{ et } \phi_3(x) = 1 - x_1 - x_2. $$

Si on utilise une partition de l'unité $P_2$, on trouve six relations de cette nature. Ceci justifie l'utilisation que l'on fait habituellement d'une partition de l'unité $P_1$ pour les fonctions singulières [MOE 99] [STA 03]. Même dans ce cas-là, l'élargissement de zone d'enrichissement va multiplier ces situations de non-unisolvance locale, ce qui explique en partie la dégradation du conditionnement.
### 4. Globalisation des degrés de liberté enrichis

Pour pallier ces inconvénients, on introduit une condition de liaison entre les degrés de liberté correspondant à l'enrichissement sur $B(x_0, R)$ en fond de fissure. Précisément, pour chaque composante du vecteur déplacement, on impose à tous les degrés de liberté relatifs à l'enrichissement singulier d'être égaux entre eux. En notant $T(R, x)$ la fonction continue et affine par morceaux qui vaut 1 aux sommets dans $B(x_0, R)$ et 0 aux autres sommets, l'espace d'approximation s'écrira ici

$$ \left\{ v_h = \sum_{i=1}^{N} a_i \varphi_i + \sum_{i \in I_H} b_i H \varphi_i + \sum_{k=1}^{2} \sum_{j=1}^{4} c_{jk} F_j T(R, x) e_k : a_i, b_i, c_{jk} \in \mathbb{R} \right\} $$

où $e_1, e_2$ représente la base canonique dans le plan. Le conditionnement s'en trouve sensiblement amélioré, mais la vitesse de convergence est moins bonne que précédemment. Sur le cas test, celle-ci est maintenant de 0.5, 1.5, 2.6 pour des polynômes de degré 1, 2, 3 respectivement.

### 5. Raccord ponctuel au bord de la zone d'enrichissement

L'étude de l'erreur d'interpolation (qui gouverne l'erreur de convergence) montre que la perte de la précision dans le cas précédent se concentre sur la zone de transition, c'est-à-dire sur la couche d'éléments situés à la frontière de la zone d'enrichissement.

L'idée est alors de revenir à la méthode originale de la partition de l'unité en approximation : la théorie assure une précision optimale [MEL 96]. On peut ici envisager une partition de l'unité à deux termes, correspondant à un recouvrement du domaine par deux sous-domaines $\Omega_1$ et $\Omega_2$, le fond de fissure étant disons à l'intérieur de $\Omega_1$ sans appartenir à $\Omega_2$. L'étude de l'erreur d'interpolation en norme $L^2(\Omega)$ montre que l'ordre de l'erreur ne dépend pas en fait de la largeur de la zone de transition $\Omega_1 \cap \Omega_2$.

D'où l'idée de faire tendre cette dimension vers zéro et d'introduire une condition de raccord (ponctuel ou intégral) entre la zone d'enrichissement $\Omega_1$ et son complémentaire $\Omega_2$. On choisit une condition de continuité du déplacement aux noeuds sur l'interface entre $\Omega_1$ et $\Omega_2$. Naturellement l'approximation sera non conforme, au sens où la continuité du champ de déplacement éléments finis enrichi n'est pas assurée partout sur l'interface.

On reprend le cas test du carré fissuré en mode I. On choisit pour $\Omega_1$ la réunion des triangles du maillage contenus dans le disque $B(x_0, R)$. Les courbes de convergence (Figure 2) montrent que l'on atteint l'optimalité. En fait, on observe même une légère super-convergence, et cela pour un coût quelque peu inférieur et un conditionnement du même ordre qu'avec XFEM classique.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Vitesse de convergence pour XFEM avec raccord ponctuel (XFEM - point-wise matching)</figcaption></figure>

Remerciements

Les auteurs tiennent à remercier J. P. Lambelin, J. J. Pesqué et B. Toson du CEA/CESTA pour de motivantes discussions, et sont reconnaissants au CEA/CESTA pour son soutien.

**6. Bibliographie**

[BEC 04] BÉCHET E., MINNEBO H., MOËS N., BURGARDT B., « Improved implementation and robustness study of the X-FEM method for stress analysis around cracks », soumis à *Int. J. Appl. Numer. Meth. Engrg.*

[LAB 04] LABORDE P., POMMIER J., RENARD Y., SALAÜN M., « High order extended finite element method for cracked domains », rapport de recherche n° 04-27, décembre 2004, laboratoire MIP, soumis à *Int. J. Appl. Numer. Meth. Engrg.*

[MEL 96] MELENK J.M., BABUŠKA I., « The partition of unity finite element method : Basic theory and applications. », *Comput. Meths. Appl. Mech. Engrg.*, vol. 139, 1996, p. 289-314.

[MOE 99] MOËS N., DOLBOW J., BELYTSCHKO T., « A finite element method for crack growth without remeshing. », *Internat. J. Numer. Meth. Engr.*, vol. 46, 1999, p. 131-150.

[STA 03] STAZI F.L., BUDYN E., CHESSA J., BELYTSCHKO T., « An extended finite element method with higher-order elements for curved cracks. », *Computational Mechanics*, vol. 31, 2003, p. 38-48.