# Implantation d'éléments X-FEM 3D quadratiques avec ou sans contact

Marcel Ndeffo, Patrick Massin, Zhengzhou Qian, Nicolas Moes

► To cite this version:

Marcel Ndeffo, Patrick Massin, Zhengzhou Qian, Nicolas Moes. Implantation d'éléments X-FEM 3D quadratiques avec ou sans contact. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722079

HAL Id: hal-01722079

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

Submitted on 2 Mar 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
# Implantation d'éléments X-FEM 3D quadratiques avec ou sans contact

Marcel NDEFFO ¹,² *, Patrick MASSIN ¹, Zhengzhou QIAN ¹, Nicolas MOES ²

1 LaMSID, UMR EDF-CNRS-CEA 8193, 1 Avenue du Général de Gaulle, 92141 Clamart CEDEX, {marcel.ndeffo*, zengzhou.qian, patrick.massin}@edf.fr

2 GeM, Ecole Centrale de Nantes, UMR CNRS 6183, 1 rue de la Noé, 44321 Nantes Cedex 3, nicolas moes@ec-nantes.fr

* Auteur correspondant

**Résumé** — L'approche X-FEM permet d'insérer dans un maillage éléments finis une surface de discontinuité sans changer le maillage et ce quelle que soit la position de la fissure. Afin d'avoir des résultats plus précis et une meilleure approximation éléments finis de la géométrie des fissures, on souhaite utiliser des éléments quadratiques. Après l'étape de construction géométrique, les enrichissements des éléments X-FEM pour prendre en compte des discontinuités ou des singularités affectent particulièrement le conditionnement de la matrice de rigidité. On propose des solutions pour capturer l'information sur la solution, malgré des petites valeurs propres dans le spectre de la matrice.

**Mots clés** — X-FEM, interpolation quadratique, géométrie quadratique, conditionnement

## 1. Introduction

L'approche X-FEM, introduite par Moës et al. [1] permet d'insérer dans un maillage éléments finis une surface de discontinuité sans changer le maillage et ce quelle que soit la position de la fissure. Ceci est rendu possible par un enrichissement local de l'approximation éléments finis utilisant la technique de la partition de l'unité qui conserve le caractère creux et défini positif de la matrice de raideur. Le problème majeur est alors l'intégration de la rigidité. On découpe donc les éléments traversés par l'interface en sous-éléments afin de trouver des intégrandes continues de part et d'autre de l'interface de discontinuité.

Afin d'avoir des résultats plus précis et une meilleure approximation par les éléments finis de la géométrie des fissures, on peut penser utiliser des éléments quadratiques [2]. Pour augmenter l'ordre de convergence avec des éléments quadratiques et retrouver un ordre de convergence optimal, il est nécessaire de monter en précision sur l'approximation des interfaces [3]. En effet, ces interfaces sont décrites au moyen de la méthode level set [4] et correspondent à l'iso-zéro d'une courbe de niveau sur le maillage. En chaque nœud on stocke sa distance par rapport à l'interface. A l'intérieur des éléments finis, la fonction level set est interpolée à l'aide des fonctions de forme. La position de l'interface dans ces éléments est donc seulement une approximation de la position exacte de cette dernière et la montée en ordre sur l'approximation est intrinsèquement liée à l'utilisation de fonctions de forme d'ordre supérieur dans l'interpolation. Une meilleure approximation de l'interface conduit à utiliser un sous découpage respectant la courbure de l'interface.

On propose ici une méthode spécifique pour construire les modèles quadratiques, différente de celles présentées jusqu'à présent [2,3]. En 2D, on fait le choix de sous-découper de part et d'autre de l'interface en sous triangles à 6 nœuds [5]. Le cadre théorique 3D est développé et on présente une construction de la surface 3D à partir des fonctions de forme élémentaire. Cette méthode repose sur la construction d'un nouveau maillage 3D quadratique pour l'intégration de la rigidité et des efforts au voisinage de l'interface afin de retrouver des ordres de convergence optimaux.
On montre qu'après l'étape de construction géométrique, les enrichissements des éléments X-FEM pour prendre en compte des discontinuités ou des singularités affectent particulièrement le conditionnement de la matrice de rigidité [6,7]. Pour de nombreux problèmes mécaniques dans Code_Aster, le conditionnement des matrices de rigidité XFEM quadratiques, est élevé, en général supérieur à 10¹⁶. Dans certains cas d'étude, on constate que les solveurs directs, robustes et éprouvés tel que MUMPS, donnent des résultats faux lors de la résolution du système matriciel issu de la discrétisation élément fini [8]. Afin de résoudre cette difficulté nous proposons une stratégie à plusieurs échelles: un pré-conditionnement local pour chaque type d'enrichissement et un pré-conditionnement global pour le système matriciel.

Par ailleurs, des multiplicateurs de Lagrange sont utilisés sur les éléments d'interface afin de pouvoir prendre en compte des conditions de contact-frottement au niveau de l'interface. En formulation mixte, une attention particulière doit être portée sur le choix de l'espace discret des inconnues représentant les pressions de contact par rapport à l'espace des champs de déplacement afin de satisfaire la condition LBB [9]. On propose dans [12] un choix particulier des dimensions de ces espaces pour les éléments quadratiques précédemment construits qui permet de satisfaire des conditions de convergence optimales pour des problèmes de contact-frottement sur interfaces courbes. Ce point ne sera pas développé dans le cadre de cette présentation mais les lecteurs intéressés trouveront des renseignements complets sur ce point dans [12].

## 2. Génération des éléments finis quadratiques XFEM

Afin de pouvoir prendre en compte des discontinuités ou des singularités en pointe de fissure, la discrétisation éléments finis du champ de déplacements devient :

$$u^h = \sum_{i \in \text{FEM}} \varphi_i a_i + \sum_{i \in \text{Heav}} H \varphi_i b_i + \sum_{i \in \text{CT}} \sum_{j=1}^{4} \varphi_i F_j c_{i,j} \quad (1)$$

Le premier terme correspond à la discrétisation standard du champ de déplacement. Le second terme prend en compte les discontinuités d'interface. Le dernier terme enfin, sur une zone localisée en pointe de fissure prend en compte la singularité des champs en mécanique de la rupture en élasticité. Quatre fonctions sont alors classiquement utilisées pour représenter cette singularité [1].

Une attention particulière doit être portée lors de l'intégration numérique des termes de rigidité et de second membre d'un élément traversé par la fissure. En effet, sur un élément traversé par une fissure, les gradients des déplacements peuvent être discontinus, et dans ce cas l'intégration numérique de Gauss-Legendre sur la totalité de l'élément n'est pas applicable. Afin de se replacer dans des conditions de régularité classiques, il convient de procéder à une intégration sur des domaines où l'intégrande est au moins continue. Pour un élément traversé par une fissure, il faut donc intégrer séparément de part et d'autre de la fissure (ceci apparaît pour la première fois dans [1] pour le 2D). Plusieurs procédures sont envisageables, et facilement mises en œuvre en 2D. Les difficultés apparaissent avec le 3D, notamment pour les éléments quadratiques.

Afin de réaliser l'intégration de part et d'autre de l'interface lorsque celle-ci est approximée de manière quadratique via une discrétisation nodale du champ des level set qui permet de repositionner numériquement l'interface, des sous-éléments quadratiques utilisés uniquement pour l'intégration numérique vont être générés. On s'appuiera pour les générer sur les intersections entre l'interface générée par le champ de level set et les éléments parents quadratiques initiaux. Les éléments parents sont décomposés dans un premier temps en sous-éléments virtuels tétraédriques quadratiques comme illustré ci-dessous. Ce sont ces sous-éléments tétraédriques qui sont à leur tour sous-découpés en tétraèdres quadratiques de part et d'autre de l'interface.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Exemple de sous découpage d'un élément héxaédrique en sous-éléments tétraédriques en 3D, puis intersection de l'interface avec les sous-éléments tétraédriques.</figcaption></figure>

Pour générer les sous-éléments quadratiques tétraédriques de part et d'autre de l'interface, on utilise la famille des points d'intersection entre l'interface et les arêtes des sous-tétraèdres issus de la découpe des éléments finis initiaux (ici des hexaèdres). On trouve ensuite les milieux des segments constitués de ces points d'intersection et des nœuds sommet, et les milieux des segments joignant les points d'intersection entre eux. Enfin les quadrangles sur les faces des tétraèdres coupées par l'interface sont découpés en sous-triangles. L'ensemble de ces points permet de générer les sous-tétraèdres de découpage comme illustré en partie inférieure de Fig.1.

# 3. Conditionnement

## 3.1. Conditionnement local

Deux types d'enrichissement pour prendre en compte des discontinuités ou des singularités affectent particulièrement le conditionnement de la matrice de rigidité [6,7]. Chaque type d'enrichissement génère un problème de conditionnement local. Nous déployons donc deux stratégies « locales » pour traiter ces deux problèmes de conditionnement XFEM.

Lorsque l'interface passe près d'un nœud, le conditionnement de la matrice de rigidité XFEM se dégrade. La fonction Heaviside devient « presque» colinéaire aux fonctions de forme, cf. Fig. 2. Par exemple, pour un seg3, on estime la matrice de rigidité locale du nœud milieu et son conditionnement, dans l'élément de référence. La matrice de rigidité locale est proportionnelle (au module d'Young près) à:

$$ B = \begin{pmatrix} \int_{-1}^{1} (\varphi_{i,x})^2 \, d\mathbf{x} & \int_{-1}^{1} \varphi_{i,x} \varphi_{i,x} H \, d\mathbf{x} \\ \int_{-1}^{1} \varphi_{i,x} \varphi_{i,x} H \, d\mathbf{x} & \int_{-1}^{1} (\varphi_{i,x} H)^2 \, d\mathbf{x} \end{pmatrix} \quad (2) $$
La valeur de la fonction de forme étant donnée par : $\varphi(x) = (1+x)(1-x)$, on en déduit que :

$$ B = \frac{8}{3} \begin{pmatrix} 1 & (1-\varepsilon)^3 \\ (1-\varepsilon)^3 & 1 \end{pmatrix} \quad \text{et} \quad \text{cond}(B) = \frac{1+(1-\varepsilon)^3}{1-(1-\varepsilon)^3} = \theta\left(\frac{1}{\varepsilon}\right) \quad (3) $$

où $\varepsilon$ est la distance entre l'interface et le nœud N2, cf. Fig.2.

Le conditionnement tend ainsi vers l'infini lorsque la distance entre l'interface et le nœud N2 tend vers zéro ce qui n'est pas acceptable. Cela correspond au fait que les deux fonctions de forme de la Fig. 2 sont quasi identiques à la limite.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Elément SEG3 coupé par une interface. Effet de la fonction H discontinue sur la fonction de forme de l'élément. Les deux fonctions de forme sont quasi identiques si l'interface se rapproche du nœud N2.</figcaption></figure>

Pour améliorer le conditionnement, on envisage une réécriture des ddls Heavisides. Les nouvelles fonctions d'enrichissement forment alors localement une base orthogonale. En particulier, l'information sur le petit support n'est pas écrasée par l'information sur le reste du support.

<table><thead><tr><th>Écriture n°1</th><th>Écriture n°2</th></tr></thead><tbody><tr><td>$\varphi(x) \rightarrow \frac{\varphi(x) + H\varphi(x)}{2}$</td><td>$\varphi(x) \rightarrow \varphi(x)$</td></tr><tr><td>$H\varphi(x) \rightarrow \frac{\varphi(x) - H\varphi(x)}{2}$</td><td>$H\varphi(x) \rightarrow \frac{\varphi(x) - H\varphi(x)}{2}$</td></tr></tbody></table>
$$B \approx \frac{4}{3} \begin{pmatrix} 2 - 3\epsilon & 0 \\ 0 & 3\epsilon \end{pmatrix}$$

$$B \approx \frac{4}{3} \begin{pmatrix} 2 & -3\epsilon \\ -3\epsilon & 3\epsilon \end{pmatrix}$$

<caption>Fig. 3. Réécriture des fonctions de forme.</caption>

Il suffit alors de mettre à l'échelle la matrice. Sur l'exemple du seg3, si on utilise le pré-conditionneur trivial suivant : $D_c \approx \sqrt{\frac{4}{3}} \begin{pmatrix} 1 & 0 \\ 0 & \sqrt{3\epsilon} \end{pmatrix}$, on aboutit à des formes du type :

<table><thead><tr><th>Écriture n°1</th><th>Écriture n°2</th></tr></thead><tbody><tr><td>$B = D_c B' D_c \text{ et } B' = \begin{pmatrix} 2 - 3\epsilon & 0 \\ 0 & 1 \end{pmatrix}$</td><td>$B = D_c B' D_c \text{ et } B' = \begin{pmatrix} 2 & -\sqrt{3\epsilon} \\ -\sqrt{3\epsilon} & 1 \end{pmatrix}$</td></tr></tbody></table>

où le conditionnement de la matrice locale pré-conditionnée $B'$, est de l'ordre de 1, avec les deux écritures proposées. Néanmoins la deuxième écriture divise le nombre d'opérations par deux, puisque les fonctions de forme habituelles ne sont pas modifiées et c'est donc celle que nous implémentons dans Code_Aster.

Pour les fonctions d'enrichissement crack-tip (décrivant la singularité en pointe de fissure), le mauvais conditionnement est principalement causé par l'enrichissement géométrique [6]. Pour améliorer le conditionnement, nous envisageons de remplacer l'enrichissement crack-tip scalaire (8 fonctions scalaires en 2D et 12 en 3D), par un enrichissement vectoriel (3 fonctions vectorielles seulement). Il s'agit d'injecter dans la modélisation XFEM, des fonctions singulières définies dans leur base locale (en fond de fissure). Néanmoins, cette stratégie d'enrichissement n'est pas encore opérationnelle dans la version 11 du Code_Aster, nous renvoyons donc aux articles [7,10] pour de plus amples explications.

## 3.2. Conditionnement global

Nous proposons ici une stratégie plus globale, pour compléter les 2 approches locales précédentes. Cette stratégie est analogue à l'algorithme de Strouboulis [11] proposé dans le cadre de GFEM. Il s'agit de changer le comportement spectral du solveur direct: le nouveau solveur récupère toute l'information des modes associés aux plus grandes valeurs propres et atténue celles venant des quasi-modes rigides. D'un point de vue mathématique, on souhaite projeter les solutions orthogonalement au noyau du problème initial. Nous allons montrer qu'on peut construire itérativement cette projection, en s'appuyant sur les propriétés spectrales du solveur direct, sans construire tacitement un projecteur.

### 3.2.1. Principe de conditionnement global en précision machine infinie

Pour résoudre l'équation linéaire $Kx = b$, avec $K$ singulier, on résout $(K + \mu I)x = b + \mu x$. On recherche de manière équivalente la solution du **problème de point fixe** suivant : $(K + \mu I)x_{n+1} = b + \mu x_n$.

Le problème de point fixe est bien posé dès que l'application $f: x \to (K + \mu I)^{-1}(b + \mu x)$ est contractante, lorsque le rayon spectral de $\nabla f(x)$ est inférieur à 1. De plus, le problème converge dans l'orthogonal du noyau, si le premier terme de la suite $x_n$ appartient à l'orthogonal du noyau. Par la suite on considère une matrice symétrique $K$ mais la démonstration peut s'étendre aux matrices non symétriques, en introduisant des valeurs propres complexes.
## Preuve

Puisque $K$ est symétrique donc diagonalisable: $K = PDP^{-1}$ et $\nabla f(x) = P((D + \mu I)^{-1}\mu I)P^{-1}$ avec $D = \text{diag}(\lambda_1, \dots, \lambda_k, \dots, \lambda_{\text{dim}(K)})$. Si le rayon spectral est inférieur à 1, alors $|\frac{\mu}{\lambda_k + \mu}| < 1$ pour toutes les valeurs propres.

En choisissant $\mu > 0$, la condition devient $\mu < \lambda_k + \mu$ soit $\lambda_k > 0$ ou $\lambda_k < -2\mu$.

* Si $K$ est définie positive, $\lambda_k > 0$ et le problème de point fixe est bien posé.

* Si $K$ est seulement positive, $\exists \lambda_k$ tel que $\lambda_k = 0$ et la recherche de point fixe échouera. Cependant l'obtention d'une telle matrice est exceptionnelle car les erreurs d'arrondis sur les coefficients de la matrice, suffisent à transformer une matrice avec une valeur propre nulle, en une matrice avec une valeur propre presque nulle.

* Si $K$ n'est pas positive, il suffit que les valeurs propres négatives soient assez grandes en valeur absolue alors $|\lambda_k| > 2\mu$.

Enfin, $\forall y \in \ker(K), Ky = 0_{\mathbb{R}^n}$, et on a le résultat suivant $y^t x_{n+1} = y^t x_n$. En effet $y^t(K + \mu I)x_{n+1} = \mu y^t x_{n+1}$ vu que $K$ est symétrique. Donc $y^t(b + \mu x_n) = \mu y^t x_n$ car $b$ appartient forcément à l'orthogonal du noyau (sinon l'alternative de Fredholm est brisée et le système n'a pas de solution). Ainsi, pour garantir que la suite reste dans l'orthogonal du noyau, il suffit que le premier terme soit dans l'orthogonal du noyau ($x_0 = 0_{\mathbb{R}^n}$ convient).

Quand $f$ est contractante, la suite $x_0, f(x_0), f \circ f(x_0), \dots, f^n(x_0)$ converge vers une unique limite $x^*$. Comme chaque terme de la suite appartient à l'orthogonal du noyau, alors la limite appartient à l'orthogonal du noyau. Puisque l'orthogonal du noyau est un espace complet, la suite $(x_n)$ converge, donc elle de Cauchy, donc elle converge dans l'orthogonal du noyau. Par unicité $x^*$ est le projeté orthogonal recherché.

### 3.2.2. Vitesse de convergence en précision machine finie

Dans le paragraphe précédent, nous avons montré dans quelles circonstances, la suite admet une limite. La relation de récurrence de la suite s'écrit alors $(K + \mu I)X_{n+1} = b + (\mu I)X_n$. Si $K$ est une matrice réelle donc diagonalisable dans $\mathbb{C}: K = PDP^{-1}D = \text{diag}(\lambda_1, \dots, \lambda_{\text{dim}(K)})$ avec $\lambda_i \in \mathbb{C}$, alors $(PDP^{-1} + \mu I)X_{n+1} = b + \mu X_n$.

Par changement de variables: $\tilde{x}_n = P^{-1}X_n$ et $\tilde{b} = P^{-1}b$, il vient: $(D + \mu I)\tilde{x}_{n+1} = \tilde{b} + \mu \tilde{x}_n$. Si $\tilde{x}_k$ désigne la k-ième composante du vecteur $\tilde{X}$ dans la base de vecteurs propres, alors, $(\lambda_k + \mu)\tilde{x}_{n+1,k} = \tilde{b}_k + \mu \tilde{x}_{n,k}$. C'est une suite (scalaire) qu'on résout par la méthode classique. Si $\lambda_k = 0$, c'est une suite arithmétique, si $\lambda_k \neq 0$, c'est une suite arithmétique-géométrique.

Finalement la solution obtenue devient,

$$ \lambda_k = 0, \qquad \tilde{x}_{n,k} = n \frac{\tilde{b}_k}{\mu} + \tilde{x}_{0,k} \qquad (4) $$

$$ \lambda_k \neq 0, \qquad \tilde{x}_{n,k} = \left(\frac{\mu}{\lambda_k + \mu}\right)^n \tilde{x}_{0,k} - \frac{\tilde{b}_k}{\lambda_k} + \frac{\tilde{b}_k}{\lambda_k} $$

## Remarques

* Le cas $\lambda_k = 0$ n'est pas singulier malgré une expression différente car l'application $\lambda_k \to \frac{\mu}{\lambda_k + \mu} \tilde{x}_{0,k} - \frac{\tilde{b}_k}{\lambda_k} + \frac{\tilde{b}_k}{\lambda_k}$ est prolongeable par continuité en 0.
* L'appartenance à l'orthogonal du noyau est moins triviale. En fait, la démonstration théorique ci-dessus ne tient pas compte de la précision machine: si $\mu$ est trop faible, $\mu I$ est négligeable et n'influence pas la matrice $K$. On estime que $\mu \gg \lambda_{min}$ pour recouvrir la plus petite valeur propre ce qui permet d'éliminer les composantes associées aux plus petites valeurs propres puisque : $\mu \gg \lambda_k \Rightarrow \tilde{x}_{n,k} = (\frac{\mu}{\lambda_k+\mu})^n(\tilde{x}_{0,k} - \frac{\tilde{b}_k}{\lambda_k}) + \frac{\tilde{b}_k}{\lambda_k} \approx (\tilde{x}_{0,k} - \frac{\tilde{b}_k}{\lambda_k}) + \frac{\tilde{b}_k}{\lambda_k} \approx \tilde{x}_{0,k}$. Par ailleurs, nous avons vu précédemment que $X_0 \in \text{Ker}(K)^\perp$ soit encore $\tilde{x}_{0,k} \approx 0$. On retrouve donc une autre conclusion théorique: si $X_0 \in \text{Ker}(K)^\perp$ la suite reste dans $\text{Ker}(K)^\perp$.

Si on fixe $\mu$ à une valeur intermédiaire $\mu = 10^{-8}$, et que l'on trace en fonction itérations, l'évolution du facteur d'atténuation $f_n(\lambda) = (\frac{10^{-8}}{\lambda+10^{-8}})^n$, on obtient la fonction de transfert de la Fig.4.

<figure><img src="image_4.png" /><figcaption>Fig. 4. Fonction de transfert de l'agorithme de résolution.</figcaption></figure>

L'algorithme se comporte donc comme un filtre passe haut dès la première itération et capte approximativement le projeté sur l'orthogonal du noyau de la solution. La « fréquence de coupure », est donnée par le paramètre $\mu$. L'information suivant les valeurs propres supérieures à $\mu$ est entièrement capturée, l'information suivant les plus petites valeurs propres est éliminée. Par contre, pour les valeurs propres intermédiaires (au voisinage de $\mu$), l'information n'est pas parfaitement captée, ni parfaitement éliminée: c'est la limitation de l'algorithme. Derrière cette difficulté, se cache la problématique de la définition numérique du noyau. Quand la queue du spectre est bruitée, la définition du noyau devient imprécise d'un point de vue spectral. Dans ce cas là, on risque de filtrer de l'information hors noyau, ce qui perturbe la solution. Dans la pratique, cette perturbation reste faible pour des choix de $\mu = 10^{-8}$.

On constate de plus une légère dépendance de la solution au paramètre $\mu$, compte tenu des composantes résiduelles suivant les petites valeurs propres. Ces composantes résiduelles, le « noyau numérique » dépendent de $\mu$ suivant l'équation $\tilde{x}_{n,k} = n \frac{\tilde{b}_k}{\mu}$ quand $\lambda_k \to 0$. On peut augmenter $\mu$ pour atténuer ces composantes, mais on court le risque de filtrer des modes qui ont une information spectrale importante. Du coup le choix de $\mu$ dépendrait fortement de la nature du spectre, ce que l'on cherche à éviter.

On propose dans ce cas d'utiliser les dérivées de (4) par rapport à $\mu$ pour corriger cette dépendance. La dérivée peut en outre être utilisée comme un outil de « monitoring » car elle transporte de
l'information sur les modes n'ayant pas convergés et les degrés de liberté impactés par ces modes. En reliant les degrés de liberté aux nœuds, on peut estimer le type des modes : modes rigides liés à l'insuffisance globale de conditions aux limites, modes locaux quasi-rigides XFEM. Comme les modes rigides impactent l'ensemble de la structure, la dérivée ne s'annulera pas sur un grand nombre de nœuds. A l'inverse, les modes XFEM étant naturellement localisés sur les nœuds enrichis, la dérivée sera nulle sur la majeure partie de la structure, mais non nulle sur les degrés de liberté XFEM.

### 3. Conclusion

Nous avons proposé pour nos éléments XFEM quadratiques la construction d'un pré-conditionnement Heaviside peu onéreux, s'appuyant sur des opérations algébriques élémentaires au niveau des degrés de liberté Heaviside. En plus de cette solution, nous formulons une amélioration du solveur direct peu coûteuse et robuste (le coût est de quelques descentes-remontées) pour résoudre les problèmes mal conditionnés, avec deux objectifs :

* fixer le conditionnement, donc les petits pivots de la matrice de rigidité, de telle sorte que le conditionnement reste fini,

* filtrer les modes quasi-rigides contribuant à la propagation d'erreurs d'arrondis dans la solution.

### Références

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

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

[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* **199**:1922–1936, 2010.

[4] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations, *J. Comput. Phys.* **79**: 12–49, 1988.

[5] A. Caron, P. Massin and N. Moës. Eléments X-FEM quadratiques en contact frottant. *Actes du 10eme Colloque National en calcul des Structures*, Giens 9-13 mai 2011.

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

[7] 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*, 64:1033–1056, 2005.

[8] 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*, DOI: 10.1007/s00466-012-0807-6.

[9] F. Brezzi., M. Fortin. Mixed and hybrid finite element methods, Springer-Verlag, 1991.

[10] S.Aluri, A.Kobayashi, M.Nakagaki. An assumed displacement hybrid finite element model for linear fracture mechanics, *Int. J. Fracture*, 11(2):257-271, 1975.

[11] T. Strouboulis, I. Babuska, K. Copps. The design and analysis of Generalized Finite Element Method, *Comp. Meth. Appl. Mech. Engng.*, 2000.

[12] G. Ferté, P. Massin, N. Moës. Higher-order interface problems with X-FEM. Part II: weak discontinuity problems, submitted.