# Développement d'une méthode éléments finis enrichis adaptée à la modélisation de structures multifissurées

Emilien Durif, Julien Réthoré, Alain Combescure

► To cite this version:

Emilien Durif, Julien Réthoré, Alain Combescure. Développement d'une méthode éléments finis en-
richis adaptée à la modélisation de structures multifissurées. 9e colloque national en calcul des struc-
tures, CSMA, May 2009, Giens, France. hal-01419925

HAL Id: hal-01419925

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

Submitted on 20 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
# Développement d'une méthode Eléments Finis enrichis adaptée à la modélisation de structures multifissurées

E. Durif, J. Réthoré, A. Combescure

INSA / LAMCOS, UMR/CNRS 5559

Bât. J. d'Alembert, 20 avenue Albert Einstein, 69621 Villeurbanne Cedex
{emilien.durif}@insa-lyon.fr

**Résumé** — Cette méthode se base sur la méthode des éléments finis étendus et les propriétés de partition de l'unité qui autorisent l'ajout de fonctions d'enrichissement. Celles-ci correspondent aux champs de déplacements analytiques au voisinage d'une fissure. Les degrés de liberté supplémentaires correspondent directement aux facteurs d'intensité des contraintes. Cette méthode a été validée en mode 1 et en mode mixte. Une étude paramétrique permet de déterminer les paramètres influents sur la précision des résultats. Enfin une étude permet d'étudier l'interaction de deux fissures en écran.

**Mots clés** — méthode des Eléments Finis Etendus, enrichissement analytique, interaction de fissures.

## 1 Introduction

La méthode des Eléments Finis Etendus (eXtended Finite Element Method) permet de créer des maillages indépendants des fissures ([4] et [1]). La propriété de partition de l'unité ([3]) autorise l'ajout de fonctions d'enrichissement dans la formulation en déplacement qui permettent de représenter les discontinuités et les singularités en pointes de fissures. Dans HX-FEM [4], on utilise des fonctions d'enrichissement de type Heavisde et des fonctions singulières qui décrivent la base de la solution analytique asymptotique en pointe de fissure de Weestergaard. Cette formulation entraîne l'ajout de degrés de liberté au problème. Ce nombre devient assez conséquent dans le cas d'une structure présentant un grand nombre de fissures. De plus dans ce type de problème, la détermination des facteurs d'intensité des contraintes est utile pour déterminer les fissures critiques dans la structure. Ceci nécessite un post traitement utilisant les intégrales d'interaction ou la minimisation par les moindres carrés, qui peuvent s'avérer compliquées lorsque des fissures interagissent. L'enrichissement proposé s'appuyant sur le champ de déplacement analytique au voisinage d'une fissure permet de déterminer directement les facteurs d'intensité des contraintes et de diminuer le nombre de degrés de liberté supplémentaires du problème.

## 2 Formulation de l'enrichissement

L'enrichissement développé s'appuie sur la solution analytique en déplacement de Griffith-Inglis ([5]). Ce modèle représente les fissures comme une somme continue de dislocations. On obtient le champ de déplacement dû à la présence d'une fissure ($Ud(x,y)$) en mode 1 et 2 au voisinage de celle-ci qui dépend des facteurs d'intensité des contraintes et des paramètres matériaux (1).
$$ \underline{U}d(x,y) = \sum_{i=1}^{2} K_i \underline{F}_i(x,y) \qquad (1) $$

Le milieu est élastique linéaire et isotrope et le champ de déplacement global est obtenu par le principe de superposition (2).

$$ \underline{U}(x,y) = \underline{U}c(x,y) + \underline{U}d(x,y) \qquad (2) $$

La partie continue, $\underline{U}c(x,y)$ est décomposée sur une base des éléments finis standards alors que la partie discontinue, $\underline{U}d(x,y)$ sur la base des fonctions d'enrichissement ($F_j(x,y)$) en utilisant les propriétés de la partition locale de l'unité [3] (équation 3).

$$ \underline{U}^h(x,y) = \sum_{i=1}^{N} N_i(x,y) \underline{u}_i + \sum_{j=1}^{2} \sum_{i=1}^{N_H} K_j^i \cdot N_i(x,y) \underline{F}_j(x,y) \qquad (3) $$

$N_i(x,y)$ sont les fonctions de forme éléments finis standards. N est le nombre de noeuds total de la structure et $N_H$ est un nombre défini de noeuds autour de la fissure. Pour chaque noeud $N_H$, il y a deux degrés de liberté supplémentaires ($K_1^i$ et $K_2^i$). Afin de limiter ce nombre, on impose pour tous les noeuds $i \in N_H$, $K_j^i = K_j$. On n'obtient alors plus que deux degrés de liberté supplémentaires par fissure qui correspondent directement aux facteurs d'intensité des contraintes. La nouvelle approximation éléments finis étendus est donnée par (4). La fonction $\Phi(x,y)$ est appelée fonction de troncature. Elle est égale à 1 sur tous les éléments qui contiennent des noeuds $N_H$ et décroît linéairement vers 0 sur les éléments adjacents (4).

$$ \underline{U}^h(x,y) = \sum_{i=1}^{N} N_i(x,y) \underline{u}_i + \sum_{j=1}^{2} K_j \underline{F}_j(x,y) \Phi(x,y) \qquad (4) $$

Cette formulation a l'avantage de faire correspondre directement les degrés de liberté supplémentaires aux facteurs d'intensité des contraintes $K_j$ et qui sont alors directement calculés lors de la résolution du système linéaire sans nécessiter un post traitement particulier.

## 3 Etude de validité de la formulation Eléments Finis Etendus analytique en mode 1 et en mode mixte

On choisit de vérifier la validité du modèle en sollicitant en traction une éprouvette comportant une fissure inclinée d'un angle $\beta$. En mode 1, $\beta$ est fixé à 0 et on fait varier la taille relative de la fissure par rapport la largeur de la structure avec le paramètre $\lambda = \frac{a}{W}$. a et W sont respectivement, la demi-longueur de fissure, et la largeur du massif sollicité en traction. [2] donne la dépendance de $K_1$ en fonction de $\lambda$ (équation 5). On fait varier $\lambda$ de 0,02 à 0,25. En mode mixte les facteurs d'intensité des contraintes en mode 1 et 2 dépendent de l'inclinaison de la fissure, $\beta$ et sont donnés par (6). Dans ce cas, on garde une longueur de fissure fixe, $\lambda = 0,1$ et on fait varier $\beta$ de 0 à $\pi$. Le matériau choisi a pour module d'Young, $E = 210\,\text{GPa}$ et un coefficient de Poisson, $\nu = 0,3$. On applique une contrainte $\sigma_A = 1\,\text{MPa}$. On utilise un maillage uniforme de 45 éléments quadrangles dans la largeur. Les dimensions de la fonction de troncature sont constantes pour les valeurs de $\lambda$: on prend une forme circulaire avec un rayon $r_b = 0,4W$. La figure 2 compare les résultats donnés par le calcul éléments finis enrichis et les formules analytiques. On remarque un bon comportement pour la sollicitation en mode mixte. On vérifie bien que $K_1$ est proche de 1 et $K_2$ nul, lorsque $\beta = 0$. Pour $\beta = \frac{\pi}{4}$, $K_2$ est maximal et égale à $K_1$. On trace aussi les valeurs des facteurs d'intensité des contraintes normés $K_1/\sigma_A\sqrt{\pi a}$ en fonction de $\lambda$ et on les compare avec celles donnés par (5). Globalement l'enrichissement éléments finis étendus se comporte comme la loi analytique.
Néanmoins une plus grande différence est notable entre les deux courbes lorsque $\lambda$ est grand. On rappelle que les différents calculs sont menés à $r_B$ constant donc à rapport $a/r_B$ variable. Ainsi lorsque la taille de la fissure se rapproche de la taille de la fonction de troncature, la précision des calculs diminue. La précision des calculs fluctue lorsque l'on prend des tailles relatives de la fonction de troncature ou de longueur de fissure non constantes. Ainsi les dimensions relatives de la fissure et de la fonction de troncature semblent avoir un impact sur la précision des résultats.

$$K_1 \approx \sigma_A \sqrt{\pi a} [1 - 0,025\lambda^2 + 0,06\lambda^4] \left[\cos\left(\frac{\pi}{2}\lambda\right)\right]^{-1/2} \quad (5)$$

$$\left\{ \begin{array}{l} K_1 = \sigma_A \sqrt{\pi a} \cos^2(\beta) \\ K_2 = \sigma_A \sqrt{\pi a} \cos(\beta) \sin(\beta) \end{array} \right. \quad (6)$$

<figure><img src="image_2.png" /><figcaption>Figure 1 – Fissure centrée et inclinée dans une sollicitation de traction.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2 – Evolution du facteur d’intensité des contraintes en mode 1 et en mode mixte en fonction de la taille relative de la fissure ($\lambda$) et de l’inclinaison de la fissure ($\beta$) : comparaison avec les modèles analytiques.</figcaption></figure>
# 4 Etude paramétrique de l'enrichissement

Les paramètres étudiés sont référencés dans le tableau 1. La convergence du modèle en fonction du nombre de degrés de liberté est étudiée avec une taille de fissure et une taille de fonction de troncature fixe. Les éléments utilisés sont des quadrangles à quatre noeuds qui forment un maillage uniforme de la structure. On étudie la convergence en imposant sur les bords de la structure les déplacements analytiques de Griffith-Inglis. On détermine par inversion de la matrice de rigidité les facteurs d'intensité des contraintes. L'erreur relative ($ε_r$) par rapport à la solution analytique est tracée en figure 3. On obtient moins de 10% d'erreur sur la détermination des facteurs d'intensité des contraintes avec au minimum dix éléments dans la largeur, sachant que le maillage est uniforme donc sans raffinement particulier au voisinage des fissures, ce qui montre la bonne précision de la méthode.

<table><caption>Tableau 1 – Liste des paramètres à étudier</caption><thead><tr><th>Paramètres</th><th>λ</th><th>γ</th><th>Nddl</th></tr></thead><tbody><tr><td>Relations</td><td>λ = <sup>a</sup>/<sub>W</sub></td><td>γ = <sup>R<sub>b</sub></sup>/<sub>W</sub></td><td>Nddl = 2 · N + 2</td></tr><tr><td>Intervalles</td><td>λ ∈ [0, 05; 0, 2]</td><td>γ ∈ [0, 2; 0, 5]</td><td>N<sub>elem</sub> ∈ [244; 51844]</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 3 – Erreur relative $ε_r$ en fonction du nombre de degrés de liberté pour l'enrichissement analytique.</figcaption></figure>

La validation du modèle dans une sollicition de traction a mis en évidence l'existence d'effets de bord. On étudie alors la variation simultanée de la taille relative de la fissure ($λ = a/W$) et de la taille relative de la fonction de troncature ($γ = r_B/W$) avec les valeurs répertoriées dans le tableau 1. Globalement lorsque $γ$ augmente et $λ$ diminue, l'erreur relative sur la détermination des facteurs d'intensité des contraintes diminue (figure 4). A $λ$ fixé, donc à longueur de fissure fixe, plus le domaine d'influence de l'enrichissement est grand, meilleure est la précision. En revanche à $r_B$ fixe, plus la fissure est grande (plus le rapport entre la taille de la fissure et le rayon de la fonction de troncature augmente) moins la précision est bonne. La taille de la fonction de troncature doit donc couvrir suffisamment les alentours de la fissure pour garantir une bonne précision. Le paramètre $γ/λ$ est donc le paramètre à contrôler pour s'assurer d'une bonne précision.
<figure><img src="image_5.png" /><figcaption>Figure 4 – Erreur relative du calcul de $K_1$ en fonction des paramètres $\gamma$ et $\lambda$</figcaption></figure>

# 5 Multifissuration

On étudie ici un massif comprenant plusieurs fissures. On reste dans le cas d'un matériau homogène, élastique linéaire et isotrope. On raisonne alors par superposition. On considère donc le problème élastique sans fissure auquel on ajoute le champ de déplacement discontinu de chaque fissure. On recherche le déplacement approximé sous la forme de (7).

$$ \underline{U}^h(x,y) = \sum_{i=1}^{N} N_i(x,y)u_i + \sum_{k=1}^{N_{\text{fiss}}} \sum_{j=1}^{2} K_{kj} \underline{F}_j^k(x,y)\Phi_k(x,y) \quad (7) $$

Le nombre de dégrés de liberté du problème comprend ceux de la partie éléments finis standard auquel on ajoute deux degrés de liberté par fissure ($K_{kj}$). Chaque fissure possède un domaine d'influence propre, défini par la fonction de troncature $\Phi_k(x, y)$. Les domaines d'influence peuvent être disjoints ou non, et à ce moment là il n'y a pas de couplage entre les termes enrichis des différentes fissures. Pour mettre en évidence le couplage entre deux fissures, on traite le cas de l'écrantage (figure 5). Le paramètre $\eta_y = \Delta_y/W$ contrôle l'écartement relatif.

La figure 6 donne l'évolution de $K_1$ en fonction de $\eta_y$. Sur cette même figure on trace l'évolution de $K_1$ dans le cas d'une fissure seule ne subissant pas l'influence d'une autre. Dans ce cas $K_1$ n'est pas constant par rapport à $\eta_y$. Ceci vient des effets de bords. Lorsque la fissure est éloignée des bords, $K_1$ devient quasiment constant. On remarque que lorsque $\eta_y$ est grand, les deux courbes sont très proches. En revanche lorsque $\eta_y$ devient petit, elles sont très différentes. Dans notre cas les domaines d'influence se coupent pour $\eta_y$ inférieur à 0,2. On remarque sur le graphe 6 qu'à partir de ce moment, $K_1$ chute de façon plus significative. Ceci met en évidence la prise en compte de l'interaction numérique des deux fissures.

# 6 Conclusion et discussion

L'enrichissement analytique proposé a été validé pour le cas d'un massif comprenant une seule fissure. Les résultats montrent que la méthode développée donnent des résultats prometteurs en terme de détermination des facteurs d'intensité des contraintes. Ils montrent aussi qu'il faut contrôler les paramètres de calculs. Les fissures doivent restées de faibles dimensions devant celles de la structure. Tous les calculs ont été réalisés avec un maillage de quadrangles uniforme. Pour encore améliorer la qualité des résultats on pourrait essayer d'affiner le maillage au voisinage des fissures. La méthode permet la détermination directe des facteurs d'intensité des contraintes ce qui
<figure><img src="image_6.png" /><figcaption>Figure 5 – Illustration du problème d’écrantage de deux fissures</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6 – Détermination de $K_1$ dans le cas de l’écrantage ou d’une fissure seule</figcaption></figure>

permet d’envisager d’utiliser cette méthode afin de gérer l’intéraction de plusieurs fissures. Cette étude montre que la méthode des Eléments Finis Etendus n’est pas seulement efficace pour traiter le cas de propagation de fissures mais elle s’adapte aussi très bien aux cas où on veut traiter le problème de multifissuration par l’utilisation de nouvelles fonctions d’enrichissement.

Références

[1] T. Belytschko, T. Black, Elastic crack growth in finite elements with minimal remeshing *Int. J. Numer. Meth. Engng.*, Éditeur, page601–page620, 1999.

[2] J. Lemaitre, J.L. Chaboche *Mécanique des matériaux solides*, Dunod, 1985.

[3] J.M. Melenk, I. Babuska. The Partition of Unity Finite Element Method Basic Theory and Applications. *Computer Methods in Applied Mechanics and Engineering*, Éditeur, page289–page314, 1996.

[4] N. Moës, J. Dolbow, T. Belytschko, A Finite Element Method for Crack without remeshing *Int. J. Numer. Meth. Engng.*, Éditeur, page131–page150, 1999.

[5] J. Weertman. *Dislocation Based Fracture Mechanics*, World Scientific, Singapore, 1996.