# Implémentation robuste pour maîtriser le conditionnement et la précision des modélisations X-FEM

Marcel Ndeffo, Patrick Massin, Nicolas Moes

► To cite this version:

Marcel Ndeffo, Patrick Massin, Nicolas Moes. Implémentation robuste pour maîtriser le conditionnement et la précision des modélisations X-FEM. 12e Colloque national en calcul des structures, CSMA, May 2015, Giens, France. hal-01517283

HAL Id: hal-01517283

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

Submitted on 3 May 2017

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
# Implémentation robuste pour maîtriser le conditionnement et la précision des modélisations X-FEM

M. Ndeffo<fnref n="1" />, P. Massin<fnref n="1" />, N. Moës<fnref n="2" />

**Résumé** — Dans le cas général d'une discontinuité positionnée arbitrairement dans un maillage, la méthode X-FEM peut présenter une convergence sous-optimale. En effet, lorsque l'interface X-FEM se rapproche des nœuds du maillage, le conditionnement augmente exponentiellement, tandis qu'une erreur numérique notable, se propage dans toute la structure. Ces limitations sont intrinsèques à la formulation X-FEM. Nous discuterons de différents remèdes à cet état de fait.

**Mots clés** — X-FEM, conditionnement, robustesse, précision.

## 1. Introduction

L'augmentation du conditionnement et la diminution de la précision, étudiés dans [2], ne sont que des symptômes de la colinéarité entre les fonctions d'enrichissement et les fonctions de forme de la partition de l'unité [4]. La colinéarité introduit des difficultés numériques lors de l'assemblage, puis lors de la résolution du système linéaire, difficultés quantifiées par le conditionnement et la précision; ces difficultés sont intrinsèques à la formulation X-FEM.

Nous envisageons une nouvelle formulation des fonctions de discontinuité, dont le domaine de définition est « pratiquement » orthogonal au domaine des fonctions de forme. En fait, les ingrédients de la nouvelle formulation se retrouvent dans les formulations [6] [5] [1], mais la nouvelle formulation introduit une autre interprétation de la discontinuité, qui a un impact notable sur le comportement de la fonction « saut » lors des calculs élémentaires.

Durant l'exposé, nous allons montrer sur quelques exemples, l'importance du gain en précision avec la nouvelle formulation. Mais la phénoménologie décrite ci-dessus est tout-à-fait générale: une étude de robustesse sur des maillages non structurés ou une interface positionnée arbitrairement, mènerait aux mêmes observations.

## 2. Limitations de l'enrichissement « saut » de déplacement X-FEM

Pour décrire une interface, la formulation standard [1] du « saut » déplacement X-FEM, s'appuie sur une fonction constante par morceaux telle que l'approximation en déplacement s'écrit :

$$u^h = \sum_{i \in I} a_i \Phi_i + \sum_{j \in I_H} b_j H \Phi_j \quad (1)$$

où **I** représente l'ensemble des nœuds du maillage ; **I_H** représente l'ensemble des nœuds enrichis par la fonction « saut » **H**, c'est-à-dire l'ensemble des nœuds dont le support est coupé par l'interface notée **Γ**. La fonction « saut » s'exprime de la manière suivante:

<footnote n="1">¹ EDF R&D, UMR EDF-CNRS-CEA 8193, {marcel.ndeffo,patrick.massin}@edf.fr</footnote>

<footnote n="2">² GeM, Ecole Centrale Nantes, {nicolas.moes}@ec-nantes.fr</footnote>
$$ H : x \rightarrow \begin{cases} +1 & \text{si } x \in \Omega^{+} \equiv \{u | lsn(u) \geq 0\} \\ -1 & \text{si } x \in \Omega^{-} \equiv \{u | lsn(u) < 0\} \end{cases} \qquad (2) $$

où « lsn » représente la distance signée à l'interface ou level-set normale.

Toutefois, lorsque l'interface Γ se rapproche des nœuds du maillage les fonctions de discontinuité $H\Phi_i$ se rapprochent des fonctions de formes $\Phi_i$ de la partition de l'unité (cf. figure 2, à gauche). Il existe un certain nombre de critères « estimatifs » pour quantifier ce phénomène et éliminer les degrés de liberté associés aux fonctions de « saut ».

## 2.1. Les critères estimatifs

Ils permettent d'estimer l'information relative entre l'intersection du support du nœud avec le domaine $\Omega^+$ et l'intersection du support du nœud avec le domaine $\Omega^-$. Par exemple, si l'information contenue dans le domaine $\Omega^+$ est très supérieure à l'information du côté $\Omega^-$, le support du nœud peut-être considéré comme appartenant totalement au domaine $\Omega^+$. Par conséquent, le nœud ne « voit » pas la discontinuité et les degrés de liberté de saut doivent-être éliminés (cf. figure 2).

Il existe diverses procédures d'estimation de cette information relative :

- une estimation géométrique ou fit-to-vertex, s'appuyant sur la distance signée représentant implicitement l'interface [3] [2],

- une estimation volumique, s'appuyant sur la comparaison directe de la mesure respective des domaines $\Omega^+$ et $\Omega^-$ [7],

- une estimation de la rigidité, s'appuyant sur la comparaison de la mesure respective des domaines $\Omega^+$ et $\Omega^-$, pondérée par la dérivée des fonctions de forme [2],

- et bien d'autres variantes [7]...

## 2.2. Les limitations des critères estimatifs

Toutefois les critiques suivantes, peuvent être faites des critères estimatifs :

- d'une part, ils ne garantissent pas le contrôle du conditionnement juste avant l'élimination. Tant que l'élimination n'est pas effectuée, le conditionnement augmente exponentiellement (cf. figure 1), ce qui peut conduire à l'arrêt du solveur direct, de manière arbitraire. Ce phénomène est exacerbé avec les éléments quadratiques (cf. tableau 1) ;

- d'autre part, une analyse de ces différents critères est faite dans [2]. Elle montre que les éliminations conduisent à une absence de convergence de l'erreur en énergie sur un exemple simple de compression homogène d'un cube traversé par une interface inclinée.

Pour contrôler le conditionnement, il est possible de mettre en place un pré-conditionneur efficace [8]. Mais cela ne permet pas de recouvrir la précision [2].

Du fait de la diminution exponentielle de l'information relative entre $\Omega^+$ et $\Omega^-$ dans la matrice de rigidité, les fonctions $H\Phi_i$ et $\Phi_i$ deviennent très rapidement colinéaires : si l'interface coupe une arête à 10% de sa longueur, jusqu'à 9 chiffres significatifs peuvent être perdus. Dans la pratique, même pour des modélisations simples, on observe une perte d'au moins 12 chiffres significatifs [2]. De plus, cette erreur se répand dans toute la structure et ne peut pas être considérée comme une erreur localisée.

Cette perte de précision est inextricablement liée à la formulation de la fonction saut (équation 2). Cette formulation implique des sommations hétérogènes au moment de l'assemblage, entre l'information du domaine $\Omega^+$ et celle du domaine $\Omega^-$.
<table><caption>Table 1 – La distance au nœud sommet le plus proche normalisée par la longueur du côté valant 10⁻⁷, on montre la dépendance du nombre de conditionnement 10⁶, par rapport au paramètre γ pour les éléments linéaires en partie supérieure et pour les éléments quadratiques en partie inférieure. En d'autres mots, l'information sur le petit support diminue dans la matrice de rigidité en 10⁻⁷. [7]</caption><thead><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>δ = γ</td><td>δ = 2γ</td><td>δ = γ</td><td>δ = 4γ</td><td>δ = 3γ</td><td>δ = γ</td><td>δ = 4γ</td><td>δ = 7γ</td></tr><tr><td>δ = 3γ</td><td>δ = 4γ</td><td>δ = 3γ</td><td>δ = 6γ</td><td>δ = 5γ</td><td>δ = 3γ</td><td>δ = 6γ</td><td>δ = 9γ</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Figure 1 – analyse du conditionnement de la matrice de rigidité, en faisant varier la distance entre l’interface et une rangée de nœuds du maillage, sur l’exemple d’un cube en compression homogène [2]. Le conditionnement est anarchique (courbe rouge) : les critères estimatifs n’agissent que sur les pics de conditionnement, mais pas avant. Le pré-conditionneur permet de contrôler le conditionnement [8] (courbe bleue).</figcaption></figure>

## 2.3. Nouvelle formulation du « saut » de déplacement

Pour représenter la cinématique de discontinuité, il est possible de s’appuyer sur un autre choix de fonction constante par morceaux $H'$ tel que [5]:

$$H': x \rightarrow \begin{cases} 0 & \text{si } x \in \Omega^{+} \equiv \{u / lsn(u) \geq 0\} \\ 1 & \text{si } x \in \Omega^{-} \equiv \{u / lsn(u) < 0\} \end{cases} \qquad (3)$$

En d'autres mots, $H'$ est la fonction caractéristique du domaine $\Omega$, que nous noterons aussi $\chi_{\Omega}(x)$. Pour les nœuds dont la position géométrique $x_i$ appartient à $\Omega^+$, la fonction de saut $H'\Phi_i$ transporte rigoureusement l'information sur le saut de déplacement, même lorsque la mesure de l'intersection du support du nœud avec le domaine $\Omega$ tend vers zéro. En effet, la fonction $H'\Phi_i$ conserve l'information sur la discontinuité et empêche les sommations hétérogènes du paragraphe précédent : $H'\Phi_i$ tend à être orthogonale à la fonction nodale $\Phi_i$, lorsque l'interface est « rasante » (cf. figure 2, à droite). De même, pour les nœuds dont la position géométrique $x_i$ appartient à $\Omega$, on envisage la fonction caractéristique du domaine $\Omega^+$, notée $\chi_{\Omega^+}(x)$, pour décrire le saut de déplacement. Cette construction de la fonction « saut » suit alors une règle de **complémentarité** (de domaines) qui peut-
être étendue aux branchements [9]. A l'aide de fonctions caractéristiques $\chi$, nous proposons de généraliser cette approche et introduire une nouvelle formulation de l'approximation en déplacement:

$$u^h = \sum_{i \in I} a_i \Phi_i + \sum_{j \in \{\Omega^+ \cap I_H\}} b_j \chi_{\Omega^+} \Phi_j + \sum_{j \in \{\Omega^- \cap I_H\}} b_j \chi_{\Omega^-} \Phi_j \quad (4)$$

$$\Delta u^h(x) = u^h(x^+) - u^h(x^-) = \sum_{j \in \{\Omega^+ \cap I_H\}} b_j \Phi_j - \sum_{j \in \{\Omega^- \cap I_H\}} b_j \Phi_j \quad (5)$$

<figure><img src="image_3.png" /><figcaption>Figure 2 – Comportement des fonctions de « saut » lorsque l'interface se rapproche des nœuds du maillage.</figcaption></figure>

### 3. Références bibliographiques

[1] N. Moës, J. Dolbow, T. Belytschko. A finite element method for crack growth without remeshing, Int. J. Numer. Meth. Engng, Volume 46, 135-150, 1999.

[2] M. Siavelis, M.L.E. Guiton, P. Massin, N. Moës. Large sliding contact along branched discontinuities with X-FEM, publié on line à Computational Mechanics, Volume 52, 201-219, 2012.

[3] K. Dréau, N. Chevaugeon, N. Moës. Studied X-FEM enrichment to handle material interfaces with higher order finite element, Computer Methods in Applied Mechanics and Engineering, Volume 199, 1922-1936, 2010.

[4] M. Ndeffo, P. Massin, Z. Qian, N. Moës. Implantation d'éléments X-FEM 3D quadratiques avec ou sans contact, CSMA, 2013.

[5] K. W. Cheng, T.P. Fries. Higher-order XFEM for curved strong and weak discontinuities, Int. J. Numer. Meth. Engng, Volume 82, 564-590, 2010.

[6] A. Hansbo, P. Hansbo. Improved implementation and robustness study of the x-fem method for stress analysis around cracks, Comput. Methods Appl. Mech. Engrg., Volume 193, 3523-3540, 2004.

[7] Documentation du Code_Aster. R7.02.12, Version 12, www.code-aster.org.

[8] E. Bechet, H. Minnebo, N. Moes, B. Burgardt. Improved implementation and robustness study of the x-fem method for stress analysis around cracks, International Journal for Numerical Methods in Engineering, Volume 64, 1033-1056, 2005.

[9] C. Daux, N. Moës, J. Dolbow et al. Arbitrary branched and intersecting cracks with the X-FEM, International Journal for Numerical Methods in Engineering, Volume 48, 1747-1760, 2000.