# Méthode hybride analytique / éléments finis étendus pour l'estimation directe des facteurs d'intensité des contraintes

Julien Réthoré, Stéphane Roux, François Hild

► To cite this version:

Julien Réthoré, Stéphane Roux, François Hild. Méthode hybride analytique / éléments finis étendus pour l'estimation directe des facteurs d'intensité des contraintes. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01422164

HAL Id: hal-01422164

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

Submitted on 23 Dec 2016

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 CC BY 4.0 - Attribution - International License
# Méthode hybride analytique / éléments finis étendus pour l'estimation directe des facteurs d'intensité des contraintes

J. Réthoré<fnref n="1" />, S. Roux<fnref n="2" />, F. Hild<fnref n="2" />

**Résumé** — On propose une approche alternative aux enrichissements X-FEM singuliers usuels qui permette une évaluation directe des facteurs d'intensité des contraintes avec une convergence quasi optimale. Il s'agit de désactiver la partie régulière de la solution dans la zone enrichie, de la remplacer par une solution analytique et de procéder par couplage de modèles.

**Mots clés** — Facteur d'intensité des contraintes, éléments finis étendus, couplage énergétique.

## 1 Introduction

Depuis une dizaine d'années, la méthode des éléments finis étendus [4] a été développée bien au delà du cadre de l'élasticité bidimensionnelle pour l'analyse de la propagation de fissures. Bien que les schémas de discrétisation soient aujourd'hui éprouvés, peu permettent une évaluation directe des facteurs d'intensité des contraintes. Nicaise et al. [5] montrent qu'utiliser un schéma de discrétisation introduisant, dans le contexte de la partition de l'unité [1], les facteurs d'intensité des contraintes $K_I$ et $K_{II}$ directement comme des inconnues du problème n'est pas satisfaisant en termes de convergence. En effet la partition de l'unité sous-jacente qui est maintenue active dans la zone enrichie réduit de manière importante le conditionnement du système.

On propose, comme solution alternative, d'introduire au voisinage de la pointe de fissure un enrichissement analytique qui permette de désactiver la partie régulière de la discrétisation par éléments finis. Ces deux modèles sont ensuite couplés en utilisant la méthode Arlequin [2]. Cette technique alternative à la méthodes des éléments finis étendus (X-FEM), est Hybride en couplant un champ Analytique à une technique X-FEM (elle est appelée HAX-FEM). Ses performances sont illustrées à l'aide d'un exemple après que nous en ayons décrit le principe.

## 2 Champ de déplacement au voisinage de la pointe d'une fissure

Si l'on considère un modèle 2D d'un corps élastique, le champ de déplacement $u_2$ peut s'écrire en notation complexe $u_2 = f(r)e^{i\theta}$ comme un développement d'une double série de Williams [6] :

$$u_2(r, \theta) = \sum_{i=I,II} \sum_{n} c_i^n \Phi_i^n(r, \theta) \quad (1)$$

<footnote n="1">¹ LaMCoS (INSA Lyon/CNRS/Université de Lyon)
20 Avenue Albert Einstein, 69621 Villeurbanne
julien.rethore@insa-lyon.fr</footnote>

<footnote n="2">² LMT-Cachan (ENS Cachan/CNRS/UPMC/PRES UniverSud Paris)
61 Avenue du Président Wilson, 94235 Cachan
{stephane.roux,francois.hild}@lmt.ens-cachan.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1 – Schématisation des domaines $\Omega_1$, $\Omega_2$, $\Omega_{12}$ et fonction de pondération $\alpha_1$ calculée sur le maillage éléments finis du modèle 1.</figcaption></figure>

où les $c_i^n$ sont les paramètre d'influence et $\Phi_i^n$ la famille des solutions élastiques de mode $i$ satisfaisant aux conditions de bords libres le long des lèvres de la fissure :

$$ \Phi_I^n(r, \theta) = r^{n/2} \left( \kappa e^{in\theta/2} - \frac{n}{2} e^{i(4-n)\theta/2} + \left(\frac{n}{2} + (-1)^n\right) e^{-in\theta/2} \right) \quad (2) $$

et

$$ \Phi_{II}^n(r, \theta) = ir^{n/2} \left( \kappa e^{in\theta/2} + \frac{n}{2} e^{i(4-n)\theta/2} - \left( \frac{n}{2} + (-1)^n \right) e^{-in\theta/2} \right) \quad (3) $$

avec $\kappa$ la constante de Kolossov, $\kappa = (3 - v)/(1 + v)$ en contrainte plane ou $\kappa = (3 - 4v)$ en déformation, $v$ étant le coefficient de Poisson. Dans ce qui suit, seuls les champs à énergie élastique finie ($n \ge 0$) seront considérés. Les champs pour $n$ pair sont continus au passage de la ligne fis-surée, alors que ceux obtenus pour $n$ impair sont discontinus. On obtient pour $n = 0$ les translations suivant les deux axes, pour $n = 2$ la composante de contrainte T et la rotation autour du front de la fissure, et pour $n = 1$ le champ de déplacement en $\sqrt{r}$ permettant de déterminer les facteurs d'intensité des contraintes. Les champs définis pour $n > 2$ permettent de rendre compte d'effets liés aux conditions aux limites.

### 3 Champ de déplacement lointain

Loin de la pointe de la fissure et de la singularité de contraintes qu'elle introduit, on utilise une description avec des éléments finis étendus du champ de déplacement :

$$ \mathbf{u}_1(\mathbf{x}) = \sum_{i \in \mathcal{N}_I} \mathbf{N}_i(\mathbf{x}) d_i + \sum_{i \in \mathcal{N}_{cut}} \mathbf{N}_i(\mathbf{x}) \mathcal{H}(\mathbf{x}) b_i \quad (4) $$

où $\mathbf{N}_i$ sont les fonctions de forme associées à des éléments finis standards supportés par les noeuds $\mathcal{N}_I$, $d_i$ les degrés de liberté correspondants. $\mathcal{H}(\mathbf{x})$ est la fonction de Heaviside généralisée qui permet grâce aux degrés de liberté associés $b_i$, supportés par les noeuds $\mathcal{N}_{cut}$ dont le support est entièrement coupé par la fissure, de décrire la discontinuité du déplacement. Ceci constitue une discrétisation par élément fini étendu usuelle utilisant les propriétés de partition de l'unité des fonctions de forme $\mathbf{N}_i$ [1].

### 4 Formulation et couplage

Le domaine d'étude $\Omega$ est découpé en 2 sous-domaines $\Omega_1$ et $\Omega_2$ (autour de la pointe de la fissure) avec une zone de recouvrement $\Omega_{12} = \Omega_1 \cap \Omega_2$ (Figure 1). Dans $\Omega_1$, respectivement $\Omega_2$, le champ de déplacement est décrit par $\mathbf{u}_1$, respectivement $\mathbf{u}_2$.
<figure><img src="image_3.png" /><figcaption>Figure 2 – Composante horizontale (a) et verticale (b) du champ de déplacement (en µm) pour la plaque fissurée pour $r_{int} = l_{rec} = 1$.</figcaption></figure>

A partir de la forme forte des équations d’équilibre, on écrit des formes faibles pondérées par une partition de l’unité :

$$ \alpha_1(\mathbf{x}) + \alpha_2(\mathbf{x}) = 1 \qquad (5) $$

$\alpha_1$ et $\alpha_2$ pondérant les énergies des domaines $\Omega_1$ et $\Omega_2$, respectivement. On obtient :

$$ a_i(\mathbf{u}_i, \mathbf{v}_i^*) = l_i(\mathbf{v}_i^*) \qquad (6) $$

où $a_i$ est la forme bilinéaire standard en élasticité pondérée par $\alpha_i$ :

$$ a_i(\mathbf{u}_i, \mathbf{v}_i^*) = \int_{\Omega_i} \alpha_i \boldsymbol{\varepsilon}(\mathbf{u}_i) : \mathbb{C} : \boldsymbol{\varepsilon}(\mathbf{v}_i^*) d\Omega \qquad (7) $$

et $l_i$ la forme linéaire :

$$ l_i(\mathbf{v}_i^*) = \int_{\partial \Omega_i} \alpha_i \mathbf{t}_d \cdot \mathbf{v}_i^* dS \qquad (8) $$

$\mathbb{C}$ étant l’opérateur de Hooke, $\boldsymbol{\varepsilon}$ le tenseur des déformations infinitésimales et $\mathbf{v}_i^*$ les champs de déplacement admissible pour chacun des modèles.

On impose en plus la compatibilité cinématique des deux descriptions via des multiplicateurs de Lagrange $\boldsymbol{\lambda}$. Le couplage s’écrit pour tout champ de multiplicateur admissible $\boldsymbol{\mu}^*$ :

$$ \Pi(\mathbf{u}_1 - \mathbf{u}_2, \boldsymbol{\mu}^*) = \int_{\Omega_{12}} \boldsymbol{\mu}^* \cdot (\mathbf{u}_1 - \mathbf{u}_2) d\Omega = 0 \qquad (9) $$

Il en résulte un système couplé qui, une fois discrétisé, s’écrit sous forme matricielle de la manière suivante :

$$ \begin{bmatrix} \mathbf{K}_1 & \mathbf{0} & \mathbf{C}_1^T \\ \mathbf{0} & \mathbf{K}_2 & \mathbf{C}_2^T \\ \mathbf{C}_1 & \mathbf{C}_2 & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{U}_1 \\ \mathbf{U}_2 \\ \mathbf{\Lambda} \end{bmatrix} = \begin{bmatrix} \mathbf{F}_1 \\ \mathbf{F}_2 \\ \mathbf{0} \end{bmatrix} \qquad (10) $$

où $\mathbf{K}_i$ sont les matrices de raideur pondérée associées à chacune des descriptions, $\mathbf{C}_i$ les termes de couplage et $\mathbf{F}_1$ les efforts extérieurs.
<figure><img src="image_4.png" /><figcaption>Figure 3 – Facteur d’intensité des contraintes $K_I$ normé en fonction de la taille du domaine couvert par le modèle analytique pour différentes valeurs de $n_{max}$. On compare la méthode proposée (a) à une approche par éléments finis étendus usuelle (b) où $r_{enr}$ désigne le rayon de la zone dans laquelle les champs sont enrichis à l’aide de composantes singulières.</figcaption></figure>

## 5 Exemple

On s’intéresse à une plaque fissurée de largeur $w$ et de longueur $L$. Une fissure de longueur $a$ au milieu de la plaque crée une intensification $K_I$ du champ de contrainte. On trouve, pour cette géométrie, une valeur théorique $K_I = C\sigma\sqrt{\pi a}$ où $C$ ne dépend que du rapport $a/w$ [3] :

$$C = 1.12 - 0.231 \left(\frac{a}{w}\right) + 10.55 \left(\frac{a}{w}\right)^2 - 21.72 \left(\frac{a}{w}\right)^3 + 30.39 \left(\frac{a}{w}\right)^4 \quad (11)$$

Dans le modèle numérique, $w$ vaut $7 \times 10^{-3}$ m, $a = 0.5w$ et $L = 17 \times 10^{-3}$ m. Le module de Young vaut 200 GPa et le coefficient de Poisson 0.3. On applique une contrainte uniaxiale $\sigma$ de 10 MPa, qui conduit à $K_{Io} = 2.98 \text{ MPa}\sqrt{m}$, valeur de référence à laquelle seront comparées différentes estimations de $K_I$. Le maillage par éléments finis est constitué de 49 × 119 quadrangles. La Figure 2 montre le champ de déplacement pour $r_{int} = \ell_{rec} = 1$ ($r_{int}$ paramétrant en taille d’élément la dimension de la zone dans laquelle le modèle par éléments finis est désactivé, et $\ell_{rec}$ l’épaisseur de la zone de couplage).

Pour différentes valeurs de $n_{max}$, ordre maximum des termes considérés dans le modèle analytique, on trace la valeur normée de $K_I$ en fonction de $r_{int}$ (Figure 3). On observe une forte dépendance à $r_{int}$ des résultats pour $n_{max} = 2$ (base de six champs décrivant les translations ($n = 0$), la contrainte T et la rotation ($n = 2$), et le champ d’intensification des contraintes ($n = 1$) qui permet l’identification des facteurs d’intensité des contraintes). En effet, à cause d’effets de taille finie, des composantes sous-singulières ($n > 2$) sont nécessaires. Ainsi, quand on augmente le nombre de termes, les résultats deviennent quasiment indépendants de la taille de la zone couverte par le modèle analytique. Pour $n_{max} = 7$ par exemple, la valeur normée de $K_I$ reste quasiment constante, variant entre 0.999 et 1.0012 (Figure 3a). Sur la Figure 3b, on reporte les résultats obtenus en utilisant une approche par éléments finis étendus classique. Quelle que soit le valeur de $n_{max}$, on observe une forte dépendance des résultats à la taille de la zone enrichie $r_{enr}$. En effet, quand celle-ci augmente, la zone d’influence de l’enrichissement augmente et les résultats s’améliorent d’abord (0.983 pour $n_{max} = 5$) puis se dégradent car le conditionnement du problème diminue fortement.

La Figure 4 montre le résultat d’une étude de convergence pour différentes valeurs de $n_{max}$ et pour les deux techniques. L’ordre de convergence de l’approche par éléments finis étendus
classique est très faible (< 0.5), ce qui est confirmé par Nicaise *et al.* [5]. Par contre, dès que $n_{max} = 3$, l'approche proposée permet d'atteindre un ordre de convergence quasi optimal ($\approx 2$).

<figure><img src="image_5.png" /><figcaption>Figure 4 – Erreur normée sur $K_I$ en fonction de la taille des éléments finis pour différentes valeurs de $n_{max}$ par la méthode proposée (a) et une approche par éléments finis étendus usuelle (b).</figcaption></figure>

## 6 Conclusion

La partition de l'unité [1] s'est avérée un outil théorique approprié pour enrichir une approximation par éléments finis. Néanmoins, si l'on attribue aux degrés de liberté additionnels une signification mécanique et que l'on requiert une grande précision sur leur détermination, des problèmes de conditionnement diminuent l'attrait de la méthode. Pour contourner ces problèmes, on propose de désactiver la partie régulière de la solution par éléments finis dans la zone enrichie en pointe de fissure. Moyennant l'adoption d'un enrichissement approprié écrit sous la forme d'une série de Williams, on propose une stratégie de résolution à deux modèles couplés grâce à la méthode Arlequin. L'approche hybride proposée s'avère extrêmement robuste et l'on obtient une convergence presque optimale.

## Remerciements

Ce travail a été effectué dans le cadre du projet PROPAVANFIS (“Méthodes avancées pour l'étude expérimentale et numérique de propagations de fissures sous chargements complexes”) soutenu par la fondation CETIM.

## Références

[1] Babuska I, Melenk J. The partition of unity method. *International Journal for Numerical Methods in Engineering* 1997 ; **40** :727-758.

[2] Ben Dhia H, Rateau G. The arlequin method as a flexible engineering design tool. *International Journal for Numerical Methods in Engineering* 2005 ; **62** :1442-1462.

[3] Bui H. *Mécanique de la rupture fragile*. Masson, 1978.
[4] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. *International Journal for Numerical Methods in Engineering* 1999; **46**(1):133–150.
[5] Nicaise S, Renard Y, Chahine E. Optimal convergence analysis for the extended finite element method. *Soumis* 2009.
[6] Williams M. On the stress distribution at the base of a stationary crack. *ASME Journal Applied Mechanics* 1957; **24**:109–114.