# Simulation multi-échelle de fissures 3D à front courbe avec estimation directe des coefficients asymptotiques

Clément Roux-Langlois, Anthony Gravouil, Marie-Christine Baietto, Julien Réthoré

► To cite this version:

Clément Roux-Langlois, Anthony Gravouil, Marie-Christine Baietto, Julien Réthoré. Simulation multi-échelle de fissures 3D à front courbe avec estimation directe des coefficients asymptotiques. 12e Colloque national en calcul des structures, CSMA, May 2015, Giens, France. hal-01516622

HAL Id: hal-01516622

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

Submitted on 2 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
# Simulation multi-échelles de fissures 3D à front courbe avec estimation directe des coefficients asymptotiques

Clément Roux-Langlois <fnref n="1" />,<fnref n="2" />, Anthony Gravouil <fnref n="1" />, Marie-Christine Baietto <fnref n="1" />, Julien Réthoré <fnref n="1" />

**Résumé** — À partir des hypothèses de la mécanique de la rupture, la méthode proposée identifie les coefficients asymptotiques le long d'un front courbe, lors de la résolution du problème considéré. Parmi ces coefficients on retrouve notamment les facteurs d'intensité des contraintes et le T-stress. Une portion du domaine autour du front est discrétisée avec les séries asymptotiques. Elle est raccordée au sens faible au reste du domaine qui est traité avec une méthode X-FEM standard. Enfin, cette méthode est insérée dans un algorithme multigrilles localisées pour améliorer son efficacité numérique.

**Mots clés** — Mécanique de la rupture, X-FEM, FIC, multigrilles localisées, fissures courbes 3D.

## 1 Introduction

Sous les hypothèses de la mécanique linéaire élastique de la rupture, Williams [10] a proposé un développement en série asymptotique des champs mécaniques au voisinage d'une fissure 2D. Les facteurs d'intensité des contraintes (FIC) associés aux termes singuliers (ordre $n = 1$) de cette série sont les paramètres dominant la propagation dans de nombreuses applications (e.g. matériaux fragiles et fatigue). Des coefficients d'ordre supérieur de cette série jouent également un rôle dans la propagation comme le T-stress (ordre $n = 2$). Il est donc intéressant de développer des méthodes de simulation qui permettent d'identifier ces coefficients pour prédire la duré de vie et la rupture d'une structure.

L'approche la plus courante consiste à calculer ces coefficients asymptotiques (FIC et T) par post-traitement des champs mécaniques obtenus par des méthodes numériques classiques. Une approche énergétique aboutit à des intégrales de ces champs mécaniques indépendantes du contour qui permettent de calculer les FIC [2] et le T-stress. Une méthode plus simple consiste à projeter les champs calculés sur les premiers ordres des séries de Williams [3]. Cependant, la singularité diminue le taux de convergence des méthodes numériques classiques et l'opération de post-traitement est longue et les erreurs se cumulent. Pour retrouver un taux de convergence satisfaisant, il faut introduire la singularité dans la discrétisation [4]. La méthode des éléments finis étendus (X-FEM) de Moës et al. [5] enrichit une discrétisation élémentaire avec la singularité et la discontinuité et permet d'utiliser des maillages non-compatibles. Néanmoins elle ne permet pas d'obtenir directement les FIC comme inconnues de la simulation.

En 2D, en utilisant les premiers ordres de la série de Williams comme discrétisation au voisinage du front il est possible d'identifier directement les coefficients asymptotiques. Différentes approches ont été développées, des super-éléments, des enrichissements avec la méthode de la partition de l'unité [11]. Finalement Réthoré et al. [7] ont montré qu'une approche basée sur la décomposition de domaine donnait de meilleurs résultats. C'est une de ces méthodes, introduite en 2D dans [6] avec un raccord intégrale que nous avons étendu en 3D [8] qui est présentée dans ce résumé. La méthode est brièvement introduite et illustrée. Elle peut être couplée avec un algorithme multigrilles localisées pour être plus efficace et permettre le raffinement nécessaire à une identification précise des coefficients asymptotiques. Ainsi, la méthode permet de prendre en compte l'échelle de la structure, celle de la fissure et celle de la singularité.

## 2 La méthode DEK-FEM en 3D

La méthode est appelée DEK-FEM pour *Direct Estimation of the generalized stress intensity factors* $K_i$ ($i = I$ ou $II$) *Finite Element Method* [6]. Deux échelles sont à prendre en compte, celle de la

<footnote n="1">¹ PRES Université de Lyon, INSA-Lyon, CNRS : UMR5259, LaMCoS, France, (anthony.gravouil,julien.rethore,marie-christine.baietto)@insa-lyon.fr</footnote>

<footnote n="2">² École polytechnique fédérale de Lausanne (EPFL), Civil Engineering Institute & Institute of Materials, Suisse, clement.roux@epfl.ch</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Description de la décomposition du domaine Ω du problème aux limites considéré.</figcaption></figure>

structure et celle de la singularité. Comme représenté sur la Figure 1, le domaine Ω est décomposé en deux sous domaines : Ω<sub>W</sub> localisé autour du front (efforts imposés négligés) et le reste du domaine Ω<sub>X</sub>. Les deux domaines ne se recouvrent pas et partagent uniquement une interface commune Γ<sub>W</sub>. Pour traiter les échelles propres à ces domaines, des modèles spécifiques y sont utilisés. Le domaine autour du front est traité avec une approximation 3D des séries asymptotiques, et le domaine Ω<sub>X</sub> avec une méthode X-FEM. Ils sont raccordés sur leur interface.

## 2.1 Domaine X-FEM Ω<sub>X</sub>

Le domaine Ω<sub>X</sub> ne contenant pas le voisinage du front est traité avec la X-FEM [1]. Le front de la fissure n'appartenant pas à ce domaine, seul l'enrichissement de la discontinuité est utilisé

$$u_X^h(x) = \sum_{k \in N} u_k \varphi_k(x) + \sum_{k \in N^d} u_{ck} \varphi_k(x) H(x) = \sum_{k \in N \cup N^d} u_k \varphi_k^g(x). \quad (1)$$

Où $\varphi_k^g(x)$ est une fonction de forme généralisée. C'est la fonction de forme éléments finis $\varphi_k(x)$ pour les degrés de liberté standards $k \in N$, et le produit de ces fonctions de forme avec la fonction d'enrichissement $H(x)$ pour les degrés de liberté enrichis $N^d$. À partir de cette discrétisation, $K_X$ la matrice de rigidité du domaine X-FEM et $F$ les efforts généralisés sont définis par

$$K_{Xkj} = \int_{\Omega_X} \nabla^s \varphi_k^g : C : \nabla^s \varphi_j^g \, dV \quad \forall (j,k) \in [1; N+N^d]^2 \quad (2)$$

$$F_j = \int_{\Omega_X} f_d \cdot \varphi_j^g + \int_{\Omega_X} F_d \cdot \varphi_j^g \, dV \quad \forall j \in [1; N+N^d] \quad (3)$$

## 2.2 Patch asymptotique

Les séries de Williams sont obtenues pour une fissure plane dans des états plan et anti-plan. L'étude du problème plan fait apparaître deux modes, le mode *I* symétrique et le mode *II* antisymétrique, auquel le problème anti-plan permet d'ajouter un troisième mode. Le champ de déplacement au voisinage du front est alors décrit par la série $U(r, \theta) = \sum_{n=0}^{\infty} \sum_{i=I,II,III} b_i^n g_i^n(r, \theta)$. Les fonctions $g_i^n(r, \theta)$ sont connues en 2D mais pas pour une géométrie quelconque en 3D.

Nous proposons d'utiliser une approximation du développement en séries 3D basée sur celui connu en 2D. Nous utilisons les premiers ordres de la série en déplacement dans un repère local autour du front et supposons une évolution continue le long du front (dont l'abscisse curviligne est notée $s$). À ce champ de déplacement nous ajoutons la troisième rotation de corps rigide

$$u_{W+}^h(r, \theta, s) = \sum_{k=1}^{N_s} \left[ \sum_{n=0}^{n_{DEK}} \sum_{i=I,II,III} b_{ink} r^{n/2} g_i^n(\theta) + b'_k r \sin(\theta) e_s \right] \varphi_k(s) = \sum_{k=1}^{N_s(3n_{DEK}+1)} b_k G_k(r, \theta, s). \quad (4)$$

Où $G_k(r, \theta, s)$ est la fonction de forme du mode et de l'ordre correspondant. À partir de cette définition locale du champ de déplacement, le gradient est calculé dans la base globale. Des déformations et des contraintes compatibles sont donc calculées. La rigidité du patch asymptotique est alors simplement

$$K_{Wkj} = \int_{\Omega_X} \nabla^s G_k : C : \nabla^s G_j \, dV \quad \forall (j,k) \in [1; N_s \cdot n_{DEK} \cdot 3]^2 \quad (5)$$
<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Représentation et résultats de l’exemple numérique présenté</figcaption></figure>

Une discrétisation éléments finis est utilisée le long de l'abscisse curviligne. Une approximation des coefficients de la série de Williams est donc identifiée en chaque nœud. De plus, l'évolution de ces coefficients est directement connue, par exemple pour les FIC $K_i(s) = \sqrt{2\pi}\sum_{s=1}^{N_s} b_{i1s}\phi_k(s)$. Ce résultat est également vrai pour les coefficients d'ordre supérieur comme le T-stress.

## 2.3 Couplage entre les domaines

Les approximations des champs mécaniques pour les deux domaines sont incompatibles sur leur frontière. Il convient néanmoins que les déplacements à l'interface soient les plus continus possible. Cette contrainte permet de coupler les deux domaines pour traiter l'intégralité du domaine $\Omega$. Les déplacements sont raccordés en moyenne sur l'interface par l'intermédiaire d'un troisième champ $\lambda$

$$ \int_{\Gamma_W} \lambda(x) (u_X(x) - u_W(x)) = 0 \quad \forall \lambda \in \mathcal{L}. \qquad (6) $$

Il s'agit du champ de multiplicateurs de Lagrange qui sont les efforts à l'interface entre le deux domaines. L'espace $\mathcal{L}$ est approché par $N_L$ fonctions $L_k(x) : \lambda^h(x) = \sum_{k=1}^{N_L} \lambda_k L_k(x)$. Les opérateurs de couplage sont donc définis comme

$$ C_{Wkj} = - \int_{\Gamma_W} L_k G_j \, dS \quad \forall (j,k) \in [1; N_s \cdot n_{DEK} \cdot 3] \times [1; N_l] \qquad (7) $$

$$ C_{Xkj} = \int_{\Gamma_W} L_k \varphi_j^g \, dS \quad \forall (j,k) \in [1; N+\mathcal{N}^d] \times [1; N_l] \qquad (8) $$

En pratique, cet espace est discrétisé comme le champ de déplacement du patch asymptotique (i.e. $N_l = N_s$ et $L_k = G_k$). Ce choix est simple d'application et donne de bons résultats.

## 3 Exemple numérique

Plusieurs cas test ont été simulés pour valider la méthode. Comme illustration, une fissure plane à front circulaire est considérée sous un chargement de tension incliné par rapport à la normal du plan de la fissure. Cette configuration est représentée sur la Figure 2a. Les trois modes de rupture sont activés et deux évoluent le long du front, pour un milieu infini ils sont connus analytiquement.

Pour simuler ce problème, un maillage de la structure (cubique avec $l = 10a$) raffiné localement autour de la fissure est considéré. Il est uniforme dans la région de la fissure et est assez grossier ($h \simeq a/7$). À partir de ce maillage, les éléments totalement ou partiellement dans un rayon $r_W = 2h$ autour du front sont désactivés du domaine X-FEM et leur support constituent le patch asymptotique. Ces éléments sont néanmoins utilisés pour l'intégration avec une quadrature de Gauss d'ordre 7. Une troncature à l'ordre
<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Représentation de la discrétisation multigrilles localisées avec le patch en rouge.</figcaption></figure>

$n_{DEK} = 7$ est considérée, 16 éléments sont considérés le long du front. Pour palier à des problèmes de couplage, une discrétisation d'ordre 1 est utilisée pour les modes plans et une discrétisation d'ordre 0 pour le mode anti-plan. Bien que grossière, cette discrétisation permet une bonne évaluation de l'évolution des facteurs d'intensité des contraintes (voir Figure 2b). Le *t*-stress mode *I* est également obtenu et représenté.

## 4 Conclusion

Cette approche permet l'évaluation directe des évolutions des coefficients asymptotiques (dont FIC et T) le long de fronts courbes en 3D. Elle repose sur une approximation des développement asymptotiques dont une troncature est utilisée comme discrétisation autour du front. Pour être précise, cette approximation nécessite une forte discrétisation autour du front et le patch asymptotique doit être de dimension réduite. Pour réunir ces conditions en 3D, la méthode a été couplée avec une approche multigrilles localisées où le patch est introduit dans la grille la plus fine (comme représenté Figure 3). Cette approche et le patch asymptotique permettent de prendre en compte les différentes échelles du problème considéré.

## Références

[1] T.-P. Fries, T. Belytschko. *The extended/generalized finite element method: An overview of the method and its applications*, Int. J. for Numerical Methods in Engineering, **84**(3), 253-304, 2010.

[2] M. Gosz, B. Moran. *An interaction energy integral method for computation of mixed-mode stress intensity factors along non-planar crack fronts in three dimensions*, Engineering Fracture Mechanics, **69**(3), 299-319, 2002.

[3] J. Lachambre. *Développement d'une Méthode de Caractérisation 3D des Fissures de Fatigue à l'aide de la Corrélation d'Images Numériques obtenues par Tomographie X*, Thèse de doctorat, INSA de Lyon, 2014.

[4] P. Laborde, J. Pommier, Y. Renard, M. Salaün. *High-order extended finite element method for cracked domains*, Int. J. for Numerical Methods in Engineering, vol. **64**(3), 354-381, 2005.

[5] N. Moës, J. Dolbow, T. Belytschko. *A finite element method for crack growth without remeshing*, Int. J. for Numerical Methods in Engineering, **46**(1), 131-150, 1999.

[6] J.-C. Passieux, A. Gravouil, J. Réthoré, M.-C. Baietto. *Direct estimation of generalized stress intensity factors using a three-scale concurrent multigrid X-FEM*, Int. J. for Numerical Methods in Engineering, **85**(13), 1648-1666, 2011.

[7] J. Réthoré, J., S. Roux, F. Hild. *Hybrid analytical and extended finite element method (HAX-FEM) : A new enrichment procedure for cracked solids*, Int. J. for Numerical Methods in Engineering, **81**(3), 269-285, 2010.

[8] C. Roux, Langlois, A. Gravouil, M.-C. Baietto, J. Réthoré. *Three-dimensional simulation of crack with curved front with direct estimation of stress intensity factors*, Int. J. for Numerical Methods in Engineering, In Press 2014 (10.1002/nme.4811).

[9] S. Roux, F. Hild. *Stress intensity factor measurements from digital image correlation : post-processing and integrated approaches*, Int. J. of Fracture, **140**(1-4), 141-157, 2006.

[10] M. Williams. *On the stress distribution at the base of a stationary crack*, ASME J. of Applied Mechanics, **24**, 109-114, 1957.

[11] Q. Z. Xiao, B. L. Karihaloo. *Direct evaluation of accurate coefficients of the linear elastic crack tip asymptotic field*, Fatigue and Fracture of Engineering Materials and Structures, **26**(8), 719-729, 2003.