# Approche Isogéométrique avec XFEM pour la mécanique de la rupture

Emmanuel de Luycker, David Benson, Ted Belytschko, Yuri Bazilevs, Ming-Chen Hsu

► To cite this version:

Emmanuel de Luycker, David Benson, Ted Belytschko, Yuri Bazilevs, Ming-Chen Hsu. Approche Isogéométrique avec XFEM pour la mécanique de la rupture. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717072

HAL Id: hal-01717072

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

Submitted on 25 Feb 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
# Approche Isogéométrique avec XFEM pour la mécanique de la rupture

Emmanuel DE LUYCKER<fnref n="1" />, David J. BENSON<fnref n="2" />, Ted BELYTSCHKO<fnref n="3" />, Yuri BAZILEVS<fnref n="2" />, Ming-Chen HSU<fnref n="2" />

**Résumé** — Ce travail propose d'investiguer l'extension des méthodes XFEM [1][2] à l'approche Isogéométrique [3][4]. Après présentation du cas test utilisé et discussion des méthodes permettant d'imposer convenablement les conditions aux limites, nous présentons les conditions permettant une convergence optimale quel que soit le degré des fonctions NURBS utilisées pour le calcul. Nous conclurons ensuite sur l'efficacité d'une telle approche en comparaison aux éléments fins classiques enrichis.

**Mots clés** — XFEM, Analyse Isogéométrique, NURBS, convergence optimale.

## 1 Introduction

Les méthodes XFEM ont prouvé leur efficacité pour la modélisation sans remaillage de discontinuités. Elles trouvent un domaine de prédilection dans celui de la mécanique de la rupture. Et permettent de garantir une convergence optimale même en présence d'une singularité de type pointe de fissure. L'objet de ce travail [5] est d'évaluer la compatibilité des méthodes XFEM avec une approche Isogéométrique. Cette dernière, permet de construire une approximation de plus haut degré tout en gardant un nombre raisonnable de degrés de liberté. Cette capacité provient du support élevé des fonctions de base lorsqu'elle sont utilisées avec une continuité élevée.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Fonctions de base pour 5 éléments en dimension 1.</figcaption></figure>

La Fig. 1 illustre cette propriété : Fig. 1(a) en degré 1 et Fig. 1(c) en degré 2 sont représentés des maillages typiques de la méthode des éléments finis (même si des polynômes de Lagrange sont plus classiques). Le tableau 1 permet de constater les différences induites par les stratégies de remaillage, celons que l'on augmente le degré avant ou après raffinement, la continuité est affectée ainsi que le nombre de points de contrôle. De manière générale, on essaie de raffiner le maillage avant d'élever le degré pour garantir une continuité maximale (sur l'ensemble du domaine) et un nombre minimal de points de contrôle, c'est le raffinement de type *k* [4].

<footnote n="1">¹ Institut Clément Ader, Université Toulouse 3, Emmanuel.De-Luycker@iut-tarbes.fr</footnote>

<footnote n="2">² University of California San Diego, Structural Engineering, dbenson@ucsd.edu</footnote>

<footnote n="3">³ Northwestern University, Mechanical Engineering</footnote>
<table><caption>Tableau 1 – Raffinement de type h, p et k [3].</caption><thead><tr><th>Degré</th><th>Multiplicité</th><th>Éléments</th><th>Continuité</th><th>Points de contrôle</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>1</td><td>C<sup>0</sup></td><td>2</td></tr><tr><td>2</td><td>2</td><td>1</td><td>C<sup>1</sup></td><td>3</td></tr><tr><td>1</td><td>1</td><td>5</td><td>C<sup>0</sup></td><td>6</td></tr><tr><td>2</td><td>1</td><td>5</td><td>C<sup>1</sup></td><td>7</td></tr><tr><td>2</td><td>2</td><td>5</td><td>C<sup>0</sup></td><td>11</td></tr></tbody></table>

L'approche Isogéométrique permet surtout de mailler exactement des géométries telles que des cercles, des tubes (Cf. Fig. 2) ou des courbes avec très peu d'éléments alors qu'un niveau élevé de raffinement est nécessaire avec des éléments finis classiques (linéaires par exemple).

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Maillage d'un tube par 12 éléments linéaires suivant le rayon et quadratiques suivant les deux autres directions (orthoradiale et longitudinale) ; champ de température grossier.</figcaption></figure>

## 2 Cas test

Nous considérons une géométrie simple fendue sur sa moitié telle que représentée Fig. 3(a).

Le champ de déplacement théorique en mode I pur donné équation (1) est imposé sur le bord du domaine.

$$ \underline{u}_{K_1}^{\infty}(r, \theta) = \frac{\sqrt{r}}{2\sqrt{2\pi}\mu} \left[ \begin{array}{l} (-\frac{1}{2} + k) \cos \frac{\theta}{2} - \cos \frac{3\theta}{2} \\ (\frac{1}{2} + k) \sin \frac{\theta}{2} - \sin \frac{3\theta}{2} \end{array} \right] \quad (1) $$

La précision du calcul est évaluée en comparant le facteur d'intensité des contraintes imposé et celui identifié à l'intérieur du domaine. Selon le type d'enrichissement, $K_I$ représente directement le degré de liberté associé à l'enrichissement en pointe de fissure ou alors, celui-ci est identifié en minimisant $J(K_I)$ définie équation (2) sur $\Omega_{\Gamma_I}$, un domaine entourant la pointe de fissure.

$$ J(K_I) = \int_{\Omega_{\Gamma_I}} (\underline{u}^h - \underline{u}_{K_I}^{\infty} K_I)^2 d\Omega \quad (2) $$

L'introduction de la méthode de Gauss pour calculer de manière discrète cette intégrale conduit à l'équation (3) dans laquelle $w_g$ et $\underline{x}_g$ sont respectivement les poids et positions des points de gauss.

$$ K_I = \frac{\sum_{g \in G_{\Gamma_I}} w_g \underline{u}^h(\underline{x}_g) \cdot \underline{u}_{K_I}^{\infty}(\underline{x}_g)}{\sum_{g \in G_{\Gamma_I}} w_g \underline{u}_{K_I}^{\infty}(\underline{x}_g) \cdot \underline{u}_{K_I}^{\infty}(\underline{x}_g)}. \quad (3) $$

Pour pouvoir intégrer correctement les champs asymptotiques en pointe de fissure, une méthode reposant sur un paramétrage quasi-polaire des sous-éléments [6] a été utilisée telle qu'illustrée Fig. 3(b).
<figure><img src="image_4.png" /><figcaption>Fig. 3 – Cas test : géométrie fissurée en mode I pur et stratégie d’intégration.</figcaption></figure>

### 3 Imposition des conditions aux limites

Les fonctions NURBS ne satisfaisant pas la propriété du delta de Kronecker, un soin particulier doit être porté à l'imposition des conditions aux limites. Une minimisation des moindres carrés est utilisée ici. Pour cela un ensemble de points $\mathcal{P}$ est considéré affin d'évaluer le champ de déplacement. La fonctionnelle $J$ est alors créée équation (4).

$$J = \frac{1}{2} \sum_{p \in \mathcal{P}} \left\| \sum_{A \in \mathcal{A}^p} N_A(s^p) q_A - u_d(x^p) \right\|^2 \quad (4)$$

C'est la somme des carrés des écarts sur le déplacement théorique (en mode I pur), $\mathcal{A}^p$ représente alors l'ensemble des indices des fonctions de base non nulles autours du point $p$.

En écrivant la stationnarité de $J$ vis à vis des paramètres $q_A$, recherchés, on obtient le système donné équation (5).

$$\sum_{p \in \mathcal{P}} \begin{pmatrix} N_1^p N_1^p & \cdots & N_1^p N_n^p \\ \vdots & \ddots & \vdots \\ N_n^p N_1^p & \cdots & N_n^p N_n^p \end{pmatrix} \begin{pmatrix} q_x^1 & q_y^1 & q_z^1 \\ \vdots & \vdots & \vdots \\ q_x^n & q_y^n & q_z^n \end{pmatrix} = \sum_{p \in \mathcal{P}} \begin{pmatrix} u_x^p N_1^p & u_y^p N_1^p & u_z^p N_1^p \\ \vdots & \vdots & \vdots \\ u_x^p N_n^p & u_y^p N_n^p & u_z^p N_n^p \end{pmatrix} \quad (5)$$

Finalement, la valeur du déplacement des points de contrôle est ensuite obtenue par inversion du système.

Une étude de convergence à permis de vérifier la précision de notre méthode d'imposition des déplacements aux bords. La Fig. 4 illustre bien la convergence optimale observée sur la norme $L_2$ après raffinement et ce pour un degré variant de 1 à 4.

### 4 Enrichissements et compatibilité

Nous cherchons ici à utiliser les capacités de XFEM pour prendre en compte la fissure sans remaillage. Nous avons donc commencé par introduire un enrichissement de type Heaviside. Ensuite, pour vérifier l'optimalité de la méthode, nous avons également introduit un enrichissement permettant de prendre en compte les champs asymptotiques en pointe de fissure. Dans notre cas, nous avons simplement ajouté le déplacement théorique en mode I pur défini équation (1) à notre approximation.

La compatibilité entre les différentes zones est un point capital qui nécessite quelques précautions dans le cas d'une approche Isogéométrique.
<figure><img src="image_5.png" /><figcaption>Fig. 4 – Convergence de l’erreur sur la norme L₂ de l’imposition des déplacements aux bords.</figcaption></figure>

Beaucoup de méthodes développées dans XFEM pour assurer la compatibilité sont basées sur l'hy-
pothèse de nœuds interpolant aux coins des élément, hypothèse qui est fondée pour des éléments fins
possédant des fonctions de base satisfaisant la propriété du delta de Kronecker. Dans le cas des NURBS
les points de contrôle jouent le rôle des nœuds élément fins mais de manière plus générale. Le shifting
par exemple qui consiste à soustraire la valeur de l'enrichissement au nœud n'est pas applicable dans
ce cas. Nous nous focalisons donc sur une technique de blending, multipliant notre enrichissement par
une fonction évanescente sur le bord du domaine enrichi. Cette approche est schématisée Fig. 5 ou B<sub>H</sub>
et B<sub>T</sub> sont des pondérations respectivement des enrichissements Heaviside et pointe de fissure. Elle per-
met d'assurer une bonne compatibilité entre la pointe de fissure et le reste du domaine (fissuré ou non
enrichi).

<figure><img src="image_6.png" /><figcaption>Fig. 5 – Compatibilité assurée par blending, représentation des fonctions B<sub>H</sub> et B<sub>T</sub>.</figcaption></figure>

Cette approche par *blending* est introduite dans l'approximation équation (6). En effet, pour garantir une convergence optimale, la zone enrichie en pointe de fissure doit rester comparable après raffinement du maillage. Il est donc déconseillé de n'enrichir que l'élément contenant la pointe de fissure qui voit sa taille réduire avec le degré de raffinement. L'allure de $B_T$ doit donc rester la même au cours du raffinement.

$$
\underline{u} = \sum_{A \in \mathcal{N}} N_A \underline{q}_A + B_H H \sum_{A' \in \mathcal{H}} N_{A'} \tilde{q}_{A'} + B_T \underline{u}_{K_1}^{\infty} \tilde{q} \quad (6)
$$

Même si la zone enrichie garde une taille constant avec cette approche, nous avons pu constater
numériquement une convergence dégradée avec les fonctions B<sub>H</sub> et B<sub>T</sub> telles que représentées Fig. 5 qui
ne sont que C<sup>0</sup> sur le domaine. L'explication vient de la différence de continuité entre la partie classique
et les enrichissements.
Pour palier ce problème, il convient de garder un même support pour les enrichissements et pour la partie continue de la solution. Cette contrainte provient du caractère hautement continu de la solution NURBS qui, contrairement à une solution élément fini classique, peut être $C^{p-1}$ (ou $p$ est le degré) sur l'ensemble du domaine. Ici, la continuité élevée de l'approximation NURBS impose une forme particulière à la fonction support de nos enrichissement. De plus, pour vérifier une convergence, le support de l'enrichissement doit rester invariant quel que soit le degré de raffinement du maillage.

$$ \underline{u} = \sum_{A \in \mathcal{N}} N_A q_A + H \sum_{A' \in \mathcal{H}} B_{A'}^{\mathcal{H}} N_{A'} \tilde{q}_{A'} + \underline{u}_{K_1}^{\infty} \sum_{A'' \in \mathcal{T}} B_{A''}^{\mathcal{T}} N_{A''} \bar{q}_{A''} \quad (7) $$

Nous proposons donc la formulation de l'équation (7) pour laquelle $N_A$ sont les fonctions d'interpolation NURBS, $q_A$, les positions des points de contrôle, $H$ est la fonction Heaviside, $\underline{u}_{K_1}^{\infty}$ est l'enrichissement en pointe de fissure, $B_{A'}^{\mathcal{H}}$ et $B_{A''}^{\mathcal{T}}$ permettent de construire le support des enrichissements en s'appuyant sur les fonctions utilisées pour le calcul. Enfin, $\tilde{q}_{A'}$ et $\bar{q}_{A''}$ sont les degrés de liberté introduits pour l'enrichissement. Ici, $\mathcal{N}$ est l'ensemble des indices des fonctions de base, $\mathcal{H}$ est l'ensemble des indices des fonctions de base utilisées comme support pour l'enrichissement de type Heaviside et $\mathcal{T}$ est l'ensemble des indices des fonctions de base utilisées comme support pour l'enrichissement en pointe de fissure.

# 5 Convergence optimale

Cette formulation permet d'atteindre une convergence optimale, testée jusqu'au degré 4 en terme de taille d'élément.

<figure><img src="image_7.png" /><figcaption>Fig. 6 – Convergence sur la norme $H_1$ pour la formulation de l'équation (7).</figcaption></figure>

La Fig. 6(a) atteste de la vitesse de convergence de la méthode. Asymptotiquement, les tendances NURBS tangentent par le haut leurs homologues éléments finis de même degré. Une conclusion hâtive qualifierait les éléments finis de plus efficaces.

En pratique, c'est le coût effectif du calcul qu'il convient de comparer. Celui-ci est très fortement lié au nombre de degrés de liberté dans le maillage, et c'est la que les NURBS présentent un très grand avantage. En terme de degrés de liberté, l'efficacité des NURBS sur les éléments fins classiques ($C^0$ entre les éléments) est clairement illustrée Fig. 6(b). Pour atteindre une erreur relative de $10^{-8}$ sur le facteur d'intensité des contraintes, avec des fonctions de degré 4, il faut presque 10 fois plus de degrés de liberté aux éléments finis ($C^0$ entre les éléments) qu'aux éléments NURBS (ayant une multiplicité de 1 et donc une continuité $C^3$ sur l'ensemble du domaine).

Une telle approche permet donc de tirer pleinement partie des deux méthodes ; l'approche Isogéométrique qui, grâce à la continuité accrue de l'approximation construite, s'avère être très efficace pour représenter
les phénomènes classiques de la mécanique des milieux continus et l'approche XFEM qui étend cette efficacité aux phénomènes discontinus ou singuliers.

Références

[1] T. Belytschko, T. Black. *Elastic crack growth in finite elements with minimal remeshing*, International Journal for Numerical Methods in Engineering, 45(5), p601-p620,1999.

[2] N. Moës, J. Dolbow, T. Belytschko. *A finite element method for crack growth without remeshing*, International Journal for Numerical Methods in Engineering, 46, p131-p150, 1999.

[3] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs. *Isogeometric analysis : CAD, finite elements, NURBS, exact geometry and mesh refinement*, Computer Methods in Applied Mechanics and Engineering, 194(39-41), p4135-p4195, 2005.

[4] J. A. Cottrell, T. J. R. Hughes, Y. Bazilevs. *Isogeometric analysis : toward integration of CAD and FEA*, John Wiley & Son, 2009.

[5] E. De Luycker, D. J. Benson, T. Belytschko, Y. Bazilevs, M.-C. Hsu. *X-FEM in isogeometric analysis for linear fracture mechanics*, International Journal for Numerical Methods in Engineering, 87, p541-p565, 2011.

[6] P. Laborde, J. Pommier, Y. Renard, and M. Salaün. *High-order extended finite element method for cracked domains*, International Journal for Numerical Methods in Engineering, 64, p354-p381, 2005.