# Méthode d'intégration robuste pour les fissures modélisées par X-FEM

Hans Minnebo, Nicolas Moes, Éric Béchet, Bertrand Burgardt

► To cite this version:

Hans Minnebo, Nicolas Moes, Éric Béchet, Bertrand Burgardt. Méthode d'intégration robuste pour les fissures modélisées par X-FEM. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01813060

HAL Id: hal-01813060

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

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éthode d'intégration robuste pour les fissures modélisées par X-FEM

H. Minnebo* — É. Béchet* — N. Moës* — B. Burgardt**

* GeM - Institut de Recherche en Génie Civil et Mécanique
École Centrale de Nantes - Université de Nantes - CNRS UMR 6183
1, rue de la Noë, BP 92101 - 44321 Nantes Cedex 3 - FRANCE

** Snecma Moteurs, Centre de Villaroche, 77550 Moissy-Cramayel

RÉSUMÉ. La prédiction de la durée de vie en fatigue de pièces fissurées n'est pas chose aisée avec la méthode des éléments finis classique. En effet, cette méthode nécessite de remailler la fissure à chaque étape de propagation. Ces actions de remaillage sont complexes en 3D et coûteuses en de temps pour les calculs industriels de fatigue. La méthode X-FEM (eXtended Finite Element Method), avec des fonctions d'enrichissement adaptées et combinées avec des level sets, va nous permettre d'introduire une fissure au sein d'un maillage sans modifier celui-ci. Cependant, l'introduction des fonctions d'enrichissement spécifiques à la pointe de fissure amène des termes singuliers à calculer pour la matrice de rigidité. Par conséquent, un schéma de Gauss-Legendre sera insuffisant pour avoir une grande précision, nous présentons donc un schéma d'intégration modifié pour les cas 2D et 3D et spécifiquement adapté aux fonctions introduites.

ABSTRACT. Fatigue lifetime prediction of cracked parts is a difficult problem with the standard finite element method. It is necessary to remesh the crack at each propagation step. This remeshing procedure is complex in 3D and take much time for industrial cases. The X-FEM, using adapted enrichment functions, together with the level set method, will enable use to embed a crack in a mesh without modifying it. Nevertheless, this new approach compels us to integrate properly singular functions. Obviously, a Gauss-Legendre scheme will not be sufficient for such a computation. We present a new modified integration scheme in 2D and 3D, which is specifically dedicated to the singular functions in X-FEM.

MOTS-CLÉS : Fissuration, singularité, partition de l'unité X-FEM

KEYWORDS: Fracture, singularity, partition of unity method, X-FEM
# 1. Introduction

L'analyse complète et précise de la propagation de fissure dans des pièces en rotation est un enjeu majeur dans l'industrie. En effet, une meilleure estimation de la durée de vie de ce type de pièces devrait permettre de réduire le besoin d'inspection qui fait perdre du temps d'utilisation. Une fissure est une entité évolutive dans le temps. Du fait de cette propriété, l'utilisation de la méthode éléments finis classique est très lourde à mettre en oeuvre. En effet, à chaque étape de propagation, un nouveau maillage, conforme à la nouvelle géométrie de la fissure est nécessaire. Nous allons nous placer dans le cadre de la mécanique linéaire élastique de la rupture. Les fonctions d'enrichissement utilisées avec la méthode X-FEM pour la fissuration présentent le désavantage d'avoir un gradient singulier, ce qui pose un problème pour l'intégration avec un schéma de Gauss. En quelque sorte, c'est l'allure des fonctions à intégrer qui va nous obliger à utiliser un schéma d'intégration non standard, que nous verrons en deuxième partie pour le cas 2D. Ensuite, pour le problème 3D, nous allons expliciter la méthode plus complexe basée sur la nature cylindrique du problème.

# 2. X-FEM et origine de la singularité

Comparée à la méthode classique, la méthode X-FEM correspond à une extension de l'espace fonctionnel. Pour la fissuration, les champs de déplacements s'écrivent de la même manière, seulement avec plus de fonctions d'interpolation. Si on appelle R l'ensemble des noeuds classiques, T (tip) l'ensemble des noeuds enrichis avec les quatre fonctions pointe de fissure et H (heaviside) l'ensemble des noeuds enrichis avec seulement la fonction discontinue.

$$ \mathbf{u}^h = \sum_{i \in R} \mathbf{N}_i \alpha_i + \sum_{i \in T} \sum_{j \in [1,4]} \mathbf{N}_i g_j \beta_{ij} + \sum_{i \in H} \mathbf{N}_i h \gamma_i \quad (1) $$

où les fonctions $g_i(r, \theta)$ sont définies de la manière suivante (cf [1]):

$$ g_i(r, \theta) = \{\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\} \quad (2) $$

Avec $r$ la distance au front et $\theta$ l'angle par rapport au plan tangent à la fissure, du côté de la matière. Ces fonctions sont obtenues à partir du développement de la solution asymptotique en pointe de fissure, dans le cas d'un solide infini coupé par un demi-plan infini, dans le cas élastique. Nous allons nous intéresser au cas 2D, en 3D, la méthode serait similaire.

Les fonctions $g_i(r, \theta)$ sont donc de la forme $\sqrt{r}f(\theta)$, où les fonctions $f$ sont harmoniques en $\theta$. Il apparaît intéressant se baser sur les coordonnées cylindriques, ce qui sera confirmé dans la partie suivante. Ainsi, les fonctions d'interpolation classiques en 2D (s, t, 1-s-t) peuvent se réécrire en fonction de ces coordonnées, pour simplifier, on prendra une fissure alignée sur un des axes.

$$ s = s_{\text{singulier}} + r\cos(\theta); \quad t = t_{\text{singulier}} + r\sin(\theta) \quad (3) $$
Les fonctions d'interpolation enrichies sont donc de la forme : $\sqrt{r}g(\theta) + r\sqrt{r}k(\theta)$.
En prenant le gradient, en coordonnées cylindriques (i.e. $d/dr$ et $\frac{1}{r}d/d\theta$), on obtient
des termes de la forme $\frac{1}{\sqrt{r}}g_1(\theta) + \sqrt{r}k_1(\theta)$. Rappelons que le gradient des fonctions
classiques donne des fonctions constantes.

Lors du calcul de la matrice raideur, on fait le produit des différents termes possibles,
on obtient donc une somme de termes de la forme :

$$ \frac{1}{r}f_1(\theta) ; \frac{1}{\sqrt{r}}f_2(\theta) ; f_3(\theta) ; \sqrt{r}f_4(\theta) ; rf_5(\theta) \quad (4) $$

où les fonctions $f_i$ sont harmoniques en $\theta$.

On a donc obtenu des termes singuliers en r. En utilisant les coordonnées cylindriques, cette singularité disparaît, mais les termes à intégrer ne sont pas des polynômes, il s'agit de fonctions harmoniques en $\theta$, faciles à calculer par une intégration classique, multipliées par des puissances de $\sqrt{r}$ pour lesquelles une méthode classique n'est pas du tout adaptées. C'est pour cela que nous avons développé un schéma d'intégration spécifique, pouvant servir aussi pour d'autres types de problèmes.

## 3. 2D, singularité et mapping

### 3.1. Triangle dans le plan, singularité en un de ses noeuds

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Triangle dans le plan, singularité en un de ses noeuds</figcaption></figure>

Dans le cas de la figure 1, la fissure s'arrête exactement sur un noeud de l'élément sans découper ce dernier. Il s'agit du cas le plus simple qui va nous permettre de débuer notre étude. On veut intégrer des fonctions singulières en $r^p, p = -1, -1/2, 0, 1/2, 1$

$$ I = \int_{\Omega_e} g(x, y)dS \quad (5) $$
En utilisant les coordonnées polaires :

$$x = x_0 + r\cos(\alpha + \beta_0); \quad y = y_0 + r\sin(\alpha + \beta_0) \qquad (6)$$

Ce que l'on peut écrire sous forme matricielle (ce qui sera intéressant pour généraliser par la suite).

$$\begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} x_0 \\ y_0 \end{bmatrix} + \begin{bmatrix} \cos(\beta_0) & -\sin(\beta_0) \\ \sin(\beta_0) & \cos(\beta_0) \end{bmatrix} \begin{bmatrix} \cos(\alpha) \\ \sin(\alpha) \end{bmatrix} \qquad (7)$$

Ce qui revient bien à :

$$R_{\beta_0} = \begin{bmatrix} \cos(\beta_0) & -\sin(\beta_0) \\ \sin(\beta_0) & \cos(\beta_0) \end{bmatrix} = \begin{bmatrix} \vec{x} \cdot \vec{n} & \vec{x} \cdot \vec{u} \\ \vec{y} \cdot \vec{n} & \vec{y} \cdot \vec{u} \end{bmatrix} = [P(\vec{n}, \vec{u}) \rightarrow (\vec{x}, \vec{y})] \quad (8)$$

On obtient pour le terme à intégrer :

$$I = \int_{\Omega_e} g(x(r, \alpha), y(r, \alpha)) r dr d\alpha \qquad (9)$$

avec $r \in [0; \frac{r_0}{\cos(\alpha)}]$; $\alpha \in [\alpha_1; \alpha_2]$ La fonction g(...).r est donc régulière. Cette opération revient, pour l'intégration numérique, à placer des points d'intégration sur un carré qui sera ensuite transformé en triangle, par dégênerescence d'un côté en un point. Les fonctions à intégrer sont donc maintenant en $r^p, p = 0, 1/2, 1, 3/2, 2$. Ces fonctions sont régulières, mais non polynomiales en r. Un changement de variable est donc recommandé sur cette variable. Nous proposons :

$$r = r_0 \frac{(1+u)^2}{4\cos(\alpha)}; \quad dr = r_0 \frac{(1+u)}{2\cos(\alpha)} du \qquad (10)$$

Ceci transforme l'intégrale en :

$$I = \int_{[-1;1] \times [\alpha_1; \alpha_2]} r_0 \frac{(1+u)}{2\cos(\alpha)} g(x(r(u, \alpha), \alpha), y(r(u, \alpha), \alpha)) r(u, \alpha) du d\alpha \quad (11)$$

Pour éliminer le terme en $\frac{1}{\cos(\alpha)}$, on procède au changement de variable (cf [2]) :

$$dt = \frac{d\alpha}{\cos(\alpha)}; \quad t = \frac{1}{2} \ln(\frac{1 + \sin(\alpha)}{1 - \sin(\alpha)}) \qquad (12)$$

Finalement, avec une dernière modification pour se placer sur un carré, on peut réécrire notre intégrale :

$$I = \int_{[-1;1]^2} r_0 \frac{(1+u)}{t_2 - t_1} g(x(r(u, \alpha(v)), \alpha(v)), y(r(u, \alpha(v)), \alpha(v))) r(u, \alpha(v)) du dv \qquad (13)$$
### 3.2. Triangle dans l'espace et singularité en un de ses noeuds.

Comme il a été démontré dans le cas précédent, on va pouvoir écrire la position d'un point du triangle de la même manière en 3D. Dans le repère local du triangle, les coordonnées sont :

$$ \begin{bmatrix} x_{loc} \\ y_{loc} \\ z_{loc} \end{bmatrix} = \begin{bmatrix} r\cos(\alpha) \\ r\sin(\alpha) \\ 0 \end{bmatrix} \qquad (14) $$

Pour écrire dans le système global, il suffit d'écrire à nouveau :

$$ \begin{bmatrix} x \\ y \\ z \end{bmatrix} = \begin{bmatrix} x_0 \\ y_0 \\ z_0 \end{bmatrix} + \left[ P(\vec{n}, \vec{u}, \vec{v}) - (\vec{x}, \vec{y}, \vec{z}) \right] \begin{bmatrix} r\cos(\alpha) \\ r\sin(\alpha) \\ 0 \end{bmatrix} \qquad (15) $$

où

$$ [P(\vec{n}, \vec{u}, \vec{v}) - (\vec{x}, \vec{y}, \vec{z})] = \begin{bmatrix} \vec{x} \cdot \vec{n} & \vec{x} \cdot \vec{u} & \vec{x} \cdot \vec{v} \\ \vec{y} \cdot \vec{n} & \vec{y} \cdot \vec{u} & \vec{y} \cdot \vec{v} \\ \vec{z} \cdot \vec{n} & \vec{z} \cdot \vec{u} & \vec{z} \cdot \vec{v} \end{bmatrix} \qquad (16) $$

Ensuite, il suffit de réaliser le même traitement que précédemment sur les coordonnées locales.

### 3.3. Triangle dans l'espace, singularité en un point du plan du triangle.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Triangle dans l'espace, singularité en un point du plan du triangle.</figcaption></figure>

Dans ce cas, la méthode consiste à procéder à une superposition en créant trois triangles fictifs, chacun des côtés avec le noeud singulier, sur lesquels sont intégrées les fonctions du vrai triangle.

## 4. Problème 3D, tétraèdre.

Les fonctions d'enrichissement sont les mêmes en 2D et 3D. En effet, le problème 3D est cylindrique. La singularité est toujours uniquement radiale. Les points d'intégration seront donc placés dans des positions préférentielles en rapport à la géométrie radiale du problème. Toutefois, plutôt que se référer à un tétraèdre, on va choisir un
<figure><img src="image_4.png" /><figcaption>**Figure 3.** *Tétraèdre coupé par un plan orthogonal au front de fissure.*</figcaption></figure>

cube : une des directions correspond à la direction de la fissure, ce qui revient à découper le tétraèdre par des plans orthogonaux au front ; ensuite, l’intégration (pondérée) sera faite sur les intersections de la même manière qu’en 2D, pour tenir compte de la singularité.

**5. Résultats et conclusion.**

<figure><img src="image_5.png" /><figcaption>**Figure 4.** *Convergence comparée pour les fonctions sur un triangle de référence.*</figcaption></figure>

La convergence est nettement accrue avec le nouveau schéma, on atteint très rapidement la précision machine. Cette convergence peut être optimisée en utilisant des nombres de points différents selon les directions. En effet, trois points seraient nécessaires radialement, environ dix le long du front de fissure, et orthoradialement, le nombre dépend de l'ouverture angulaire.

**6. Bibliographie**

1. N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, Int. J. Num. Meth. Eng., 46, (1999), 131-150.

2. E. Béchet, H. Minnebo, N. Moës, B. Burgardt, Improved implementation and robustness study of the X-FEM method for stress analysis around cracks, Int. J. Num. Meth. Eng. (2004), accepté.