# Simulation d'un milieu fissuré tridimensionnel avec estimation directe des facteurs d'intensité des contraintes

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

► To cite this version:

Clément Roux, Anthony Gravouil, Marie-Christine Baietto, Julien Réthoré. Simulation d'un milieu fissuré tridimensionnel avec estimation directe des facteurs d'intensité des contraintes. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717828

HAL Id: hal-01717828

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

Submitted on 26 Feb 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
# Simulation d'un milieu fissuré tridimensionnel avec estimation directe des facteurs d'intensité des contraintes.

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

**Résumé** — L'objet du travail présenté est le développement d'un simulation d'un milieu fragile fissuré afin d'obtenir directement et précisément les facteurs d'intensité des contraintes (FIC). Cette méthode généralise la stratégie bidimensionnelle DEK-FEM, développée précédemment par J.C. Passieux et al. [15], aux problèmes tridimensionnels. Le principe est de décomposer le domaine étudié en deux, une partie analytique autour du front et le reste du domaine simulé avec des éléments finis étendus. En pointe de fissure, le développement asymptotique du champ de déplacement est utilisé comme discrétisation, aussi les degrés de liberté associés sont directement proportionnels aux coefficients asymptotiques. Les champs asymptotiques étant bidimensionnels, l'hypothèse d'une évolution continue le long du front est introduite pour son utilisation en trois dimensions. Les deux domaines ont des discrétisations incompatibles sur leur interface, ils sont raccordés en moyenne avec une méthode de Mortar [4]. Enfin, la singularité en pointe de fissure étant localisée, la stratégie est rendu compatible avec une méthode multi-grille localisée. Cette approche permet de ne pas recourir à un post-traitement des champs mécaniques, tel que l'intégrale d'interaction [17], pour obtenir les facteurs d'intensité des contraintes. Le mode anti-plan est intégré aux champs analytiques, et les autres termes de la séries comme la contrainte T sont également obtenus directement.

**Mots clés** — Mécanique de la rupture, XFEM, FIC, multi-grille localisée

# 1 Introduction

L'étude présentée porte sur la fissuration de milieux fragiles pour lesquels les phénomènes non-linéaires sont localisés en pointe de fissure, ce qui est le cas pour de nombreux métaux à basses températures. Elle a donné lieu à la mécanique linéaire élastique de la rupture [8] qui est le cadre de ce travail. Cette théorie est utilisée pour simuler le comportement de milieux continus fissurés, et les caractériser. Un développement asymptotique en pointe de fissure en est extrait [12], et est utilisé comme discrétisation autour du front [15] en complément d'une approche X-FEM [10]. La simulation de fissure fragile nécessite de concentrer l'effort de discrétisation autour du front, ce qui est réalisé par une approche multi-grille localisée [16].

La mécanique de la rupture, ou l'étude de milieux fissurés a débuté par une approche énergétique (griffith et Irwin) proposant un critère de propagation basé sur le taux de restitution d'énergie [8]. Il existe plusieurs méthodes pour évaluer ce taux, il peut notamment être relié au facteur d'intensité des contraintes. Ces facteurs caractérisent la singularité en pointe de fissure et sont issus de l'étude d'un milieu continu présentant une discontinuité [12]. Cette étude mène à l'expression des champs mécaniques en pointe de fissure sous la forme d'une série infinie, c'est à partir de cette série que la discrétisation autour du front est obtenue.

La simulation des milieux fissurés présente des difficultés dues à la prise en compte de la singularité et aux phénomènes non-linéaires qui interviennent autour du front. Cette singularité dégrade la précision et le taux de convergence de la méthode des éléments finis. L'introduction d'éléments singuliers [3] permet d'en améliorer la précision, notamment pour l'évaluation des facteurs d'intensités des contraintes. Enfin, la méthode des éléments finis étendus [10] améliore le taux de convergence. Cette méthode est basée sur la propriété de partition de l'unité de la méthode éléments finis [2], et permet de prendre une fissure en compte sans qu'elle soit compatible avec le maillage.

<footnote n="1">¹ PRES Université de Lyon, INSA-Lyon, CNRS : UMR5259, France, LaMCoS, {clement.roux,anthony.gravouil,julien.rethore,marie-christine.baietto}@insa-lyon.fr</footnote>
À partir des champs simulés le calcul, de l’intégrale d’interaction permet d’obtenir les facteurs d’intensité des contraintes [11]. Cependant l’utilisation des champs asymptotiques dans la discrétisation permet d’obtenir directement ces facteurs, dans un premier temps sur un maillage conforme [6], puis en utilisant la partition de l’unité [9]. La méthode utilisée ici permet d’obtenir ces coefficients et les autres termes de la série plus précisément en introduisant un domaine analytique en pointe de fissure, raccordé sur son interface [15]. Néanmoins la précision de l’identification est dépendante de la finesse de la discrétisation en pointe de fissure. Afin de limiter le coût du calcul en conciliant l’échelle structurelle et l’échelle de la singularité, la méthode a été couplée avec une approche multi-grille localisée [16].

## 2 Principe de la méthode

L’étude asymptotique des champs mécanique en pointe de fissure, section 2.1, montre que ceux-ci sont singuliers localement. Cette partie du domaine est donc traitée séparément, avec des champs spécifiques permettant d’en capter la singularité. Le domaine est donc décomposé en deux parties traitées séparément. Une approche X-FEM est utilisée pour la partie avec les lèvres de la fissure qui sont simplement considérées comme une discontinuité du champ de déplacement. Les deux parties sont raccordées en moyenne par une méthode de Mortar, sur une interface sans recouvrement.

### 2.1 Champ asymptotique au voisinage du front

Pour des chargements et des géométries de milieux fissurés simples, il est possible de déterminer les champs mécanique autour du front de la fissure à partir des équations de la mécanique. Ce problème est présenté brièvement ci-dessous pour une fissure rectiligne semi-infinie pour un chargement lointain dans un état d’élasticité plane et dans le cas d’un chargement antiplan. Ce raisonnement a été mené par [12, 18], et est présenté dans [5]. Dans une zone suffisamment proche du front, la géométrie s’approche de ces conditions. Cependant si la fissure est courbe, on peut être amené à devoir considérer un domaine analytique petit pour que la fissure puisse être considérée comme rectiligne.

#### 2.1.1 Développement asymptotique 2D

Au voisinage d'une fissure rectiligne, portée par $x$ de normale $y$, l'équilibre mécanique $\text{div}(\boldsymbol{\sigma}) = 0$ et les équations de compatibilité des déformations dans un état de déformations planes impliquent qu'il existe un potentiel $\phi$ appelé fonction d'Airy vérifiant $\triangle \triangle \phi = 0$. En cherchant une solution de la forme $\phi(r, \theta) = r^{\alpha+2}\psi(\theta)$ à cette équation, $r$ et $\theta$ sont les coordonnées polaires en pointe de fissure, le problème devient d'une équation différentielle d'ordre quatre linéaire homogène en $\psi$. La résolution de cette équation sous la condition que les lèvres soient libres d'efforts, permet d'exprimer les contraintes sous la forme d'une double série infinie. La relation de comportement élastique linéaire permet d'obtenir le champs de déplacement, noté dans le plan complexe $U = u_x + iu_y$, sous la forme d'une double série :

$$ U = \sum_{n=-\infty}^{\infty} \sum_{i=I,II} a_i^n \varphi_i^n(r, \theta), \text{ avec } \left\{ \begin{array}{l} \varphi_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], \\ \varphi_{II}^n(r, \theta) = i 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], \end{array} \right. \quad (1) $$

avec $\kappa = 3 - 4\nu$ en déformations planes. Le mode I est un mode symétrique, par rapport à $x$ tandis que le mode II est un mode antisymétrique.

Les termes de ces séries pour $n < 0$ induisent une énergie de déformation élastique infinie et sont donc exclus de la solution. Les termes d'ordre 0 (i.e. $n = 0$) correspondent aux deux translations de solide rigide, et le terme $II$ d'ordre 2 à la rotation autour du front. Les termes d'ordre 1 présentent des contraintes proportionnels à $1/\sqrt{r}$, qui tend vers l'infini à la pointe de la fissure (singulier), ils caractérisent la singularité et les coefficients $a_i^1$ associés sont proportionnels aux facteurs d'intensité des contraintes. Cette approche donne de bon résultats tant que la plasticité reste confinée devant la zone d'influence de ces termes. Enfin le terme symétrique d'ordre 2 donne accès à la contrainte T (T-stress).

#### 2.1.2 Mode III antiplan

Une approche complémentaire consiste à s'intéresser au problème antiplan, pour lequel $u = u_3(x_1, x_2)z$. Cette approche idéalisée, associée aux équations d'équilibre induit que $u_3$ doit être une fonction har-
monique (i.e., il existe une fonction holomorphe $f(z)$ tel que $\mu u_3 = \text{Re}[f(z)]$). Un retour vers les contraintes impose $\sigma_{13} - i\sigma_{23} = f'(z)$. Si l'on considère des solutions de la forme $f(z) = (A_n + iB_n)r^{n/2}e^{in\theta/2}$, l'hypothèse que les lèvres de la fissures sont libres d'efforts permet de mettre le champ de déplacement sous la forme d'une série infinie :

$$ u_3 z = \sum_{n=0}^{\infty} a_{III}^n \varphi_{III}^n(r, \theta), \text{ avec } \begin{cases} \varphi_{III}^n(r, \theta) = r^{n/2} \cos(n\theta)z & \text{si } n \text{ paire,} \\ \varphi_{III}^n(r, \theta) = r^{n/2} \sin(n\theta)z & \text{si } n \text{ impaire.} \end{cases} \quad (2) $$

Le terme d'ordre 0 correspond à la translation de solide rigide hors plan $(x, y)$. Le terme d'ordre 2 correspond à un rotation antiplane autour de la normale à la fissure passant par le front. Comme précédemment, le terme d'ordre 1 correspond à une singularité et est proportionnel au troisième facteur d'intensité des contraintes.

### 2.1.3 Discrétisation le long du front

Les développements, présentés ci-dessus ne sont valables que pour des géométries bidimensionnelles. Dans le cas général, il convient de les aménager. Des travaux [1] sont en cours pour mettre en place des développements asymptotiques en trois dimensions, mais ils n'ont pour le moment pas abouti. Les développements asymptotiques en situation plane et antiplane sont considérés afin de pouvoir caractériser les trois modes de ruptures. De plus, la troisième dimension impose une évolution des coefficients asymptotiques le long du front, cette évolution est supposée continue.

Afin de pouvoir traiter cette évolution numériquement, l'évolution suivant l'abscisse curviligne du front est discrétisée au sens des éléments finis, avec $\varphi^s(s)$ pour fonctions de forme. Le champs de déplacement, réunissant les modes I, II et III présentés équation (1) et (2) prend donc la forme :

$$ u_W(\mathbf{X}) = u_W(r, \theta, s) = \sum_{s} \sum_{n=0}^{n_{max}=\infty} \sum_{i=I,II,III} (a_i^n)^s \varphi_i^n(r, \theta) \varphi^s(s). \quad (3) $$

L'abscisse curviligne est explicitée au moyen d'une troisième fonction de niveau, confère partie 2.2. Les trois développements en série comprennent trois translations de corps rigide et deux rotations pour s'adapter aux déplacements de structure, néanmoins il manque une dernière rotation qui est ajoutée. Le domaine analytique se comporte de manière analogue à une poutre, l'information étant concentrée sur les nœuds FEM de la discrétisation en s. Le passage à l'information volumique est découplé et se fait grâce aux fonctions du champ asymptotique. En pratique, la série utilisée est une troncature de celle présentée équation (3), pour une bonne identification des termes singuliers les termes conservés vont jusqu'à $n_{max} = 7$, confère [15].

## 2.2 Décomposition du domaine

La singularité due à la fissure est traitée spécialement avec les champs analytiques présentés précédemment. Les coefficients $a_i^n$ associés à ces champs s'approchent des facteurs d'intensité des contraintes et des autres paramètres de mécanique de la rupture. Ils sont utilisés comme degrés de liberté, et sont donc directement obtenus lors de la résolution. Cependant, au voisinage du front, lorsque le développement asymptotique est utilisé en parallèle de la discrétisation éléments finis [19], les degrés de libertés associés aux fonctions asymptotiques ne sont pas précisément les coefficients asymptotiques. Aussi, le développement asymptotique est considéré seul en pointe de fissure pour une identification précise [15].

Le domaine considéré $\Omega$, est donc divisé en une partie $\Omega_W$ au voisinage du front et son complémentaire $\Omega_X$, comme illustré figure 1. La décomposition étant arbitraire, il convient de raccorder les déplacements $u_X$ (de $\Omega_X$) à ceux $u_W$ (de $\Omega_W$). Le raccordement se fait sur leur interface commune $\Gamma_W$, les champs de déplacements n'étant pas compatibles le raccord est fait en moyenne au sens de Mortar. Le raccord est pris en compte dans la rigidité par des multiplicateurs de Lagrange.

Le domaine $\Omega_W$ est discrétisée avec le champ asymptotique tronqué, tandis que le domaine $\Omega_X$ est traité dans un cadre X-FEM. Le front n'étant pas présent dans la partie $\Omega_X$, seul des enrichissements discontinus sont utilisés. Le support de la discontinuité est stocké sous la forme d'une fonction de niveau [14] sur un maillage structuré annexe. Cette fonction est proportionnelle à la distance avec la fissure. Afin de calculer $r$ et $\theta$ dans la zone analytique, une seconde fonction permet de déterminer la position
du front, elle est orthogonale à la première et proportionnelle à la distance au front. Enfin l'utilisation de l'abscisse curviligne dans $\Omega_W$ nous amène à introduire une troisième fonction de niveau orthogonale aux deux autres.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Représentation de la décomposition de domaine entre domaine analytique et zone X-FEM.</figcaption></figure>

Cette décomposition introduit un paramètre géométrique : le rayon du domaine analytique $r_W$. Dans la situation bidimensionnelle, il a été montré que $r_W = 2$ éléments X-FEM permet d'obtenir précisément les facteurs d'intensité des contraintes [15]. Pour traiter le cas général d'une fissure courbe avec le champ asymptotique considéré, il faut que la zone de rayon deux éléments autour du front soit suffisamment petite pour que la fissure puisse y être considérée comme rectiligne. Pour cette raison et pour augmenter la précision de l'identification, le maillage de la zone X-FEM peut avoir besoin d'être raffiné autour de la pointe de fissure. Ce raffinement localisé peut être obtenu par une approche multi-grille localisée, cette approche est présentée à la section 4.

# 3 Illustration de la méthode

## 3.1 Barreau sous chargement symétrique

Dans un premier temps la méthode est appliquée à l'étude d'un problème symétrique (par rapport à la fissure). Ce problème n'active que le mode I. L'étude porte sur un barreau de $7 \times 18 \times 7 \text{ mm}^3$ soumis à une contrainte de traction $\sigma \cdot n = 10 \text{ MPa}$, le matériau considéré est un acier de module d'Young $E = 200 \text{ GPa}$ et de coefficient de Poisson $\nu = 0.3$. Ce barreau, représenté figure 2(a), est fissuré dans son plan de symétrie, et le front de la fissure est inclinée (i.e. $a_i = 0,7 \text{ mm}$ et $a_f = 3,4 \text{ mm}$) afin que les paramètres de rupture évoluent le long du front.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – (a) Représentation du problème traité : géométrie et CL, (b) maillage déformé avec champs de déplacement et (c) évolutions de $K_I$ et $T$ normalisés pour deux discrétisations de $\Omega_W$.</figcaption></figure>

Le domaine $\Omega_X$ est discrétisé sur un maillage d'éléments cubiques à huit nœuds d'arête $dh = 1/3 \text{ mm}$, cette discrétisation apparaît figure 2(b). Pour le domaine $\Omega_W$, l'évolution le long du front est dis-
crétisée avec une longueur élémentaire $ds$. La simulation est menée en parallèle pour plusieurs de ces longueurs.

Comme illustré figure 2(c), pour deux discrétisations du front $ds = 0.38 \text{ mm}$ et $ds = 0.8 \text{ mm}$, l'évaluation des coefficients asymptotiques est peu dépendante de la discrétisation du front. Cette propriété est vérifiée pour les longueurs $ds > dh_{X-FEM}$. Cette simulation permet de capturer l'évolution de $K_I$ et $T$ le long du front. Comme illustré à la figure 2(c), $K_I$ augmente avec l'allongement de la fissure, le mode I est bien prépondérant : $\frac{K_{II}}{K_I} < 1\%$ et $\frac{K_{III}}{K_I} < 1\%$. Malgré les difficultés conceptuelles liées au développement asymptotique en pointe de fissure au niveau d'un bord libre, on constate que la chute de $K_I$ et/ou l'augmentation de la triaxialité des contraintes (représentée par la contrainte $T$) sont cohérentes avec une propagation plus lente au niveau des bords libres.

## 3.2 Barreau en flexion trois points

Afin de solliciter tous les modes singuliers, un barreau d'acier ($E = 200 \text{ GPA}$, $v = 0.3$) soumis à de la flexion trois points est considéré. La géométrie d'ensemble est la même que précédemment (i.e. $7 \times 18 \times 7 \text{ mm}^3$) mais la fissure est de longueur $a = 2.33 \text{ mm}$ constante et le front est incliné d'un angle $\gamma = 30^\circ$, comme représenté figure 3(a). L'effort appliqué est $F_x = 500 \text{ N}$, il est repris sur deux appuis distants de $H_u = 16.7 \text{ mm}$. La discrétisation est également des éléments à huit nœuds d'arête $dh = 1/3 \text{ mm}$.

<figure><img src="image_4.png" /><figcaption>Fig. 3 – (a) Représentation du problème traité : géométrie et CL, (b) maillage déformé avec champs de déplacement et (c) évolutions de $K_I$, $K_{II}$ et $K_{III}$ et $T$ normalisés pour deux discrétisations de $\Omega_W$, $ds$ en mm.</figcaption></figure>

Dans cette configuration, 3(c) l'estimation des coefficients asymptotiques est également peu dépendante de la discrétisation le long du front. $K_I$ présente une évolution faible puisque la longueur de fissure est constante, il diminue près des bords où le chargement est moins critique. $K_{II}$ est bien antisymétrique et présente une évolution en accord avec [7]. De même $K_{III}$ évolue d'une manière similaire à celle constatée dans [7], évolution faible puisque l'inclinaison de la fissure est constante.

## 4 Couplage avec une approche multi-grille localisée

La méthode est couplée à une approche multi-grille localisée afin de pouvoir localement augmenter la précision de la discrétisation autour du front, et pour le cas des fissures quelconques évoqué précédemment. Ce couplage consiste à mettre en place une méthode multi-grille localisée [16], et à introduire le domaine $\Omega_W$ sur la grille la plus fine. Cette approche donne une précision du niveau de celle de la grille la plus fine, mais le solveur itératif, capturant la solution d'ensemble à moindre coût (résolution sur la grille grossière), converge plus rapidement.

La méthode multi-grille consiste à approcher le problème sur la grille grossière puis à affiner la solution sur les grilles plus fines avant de remonter ces informations vers la grille la plus grossière. Il est donc nécessaire de mettre en place des opérateurs de passage entre les grilles. Ces opérateurs sont obtenus par une interpolation directe termes à termes, qui a été montrée efficace dans [16]. Au niveau de la grille la plus fine décomposée en $\Omega_X$ et $\Omega_W$, ces opérateurs sont obtenus par une approche en moyenne de type Mortar, en effet les champs de déplacement ne sont pas compatibles.
## 4.1 Illustration de l'approche multi-grille

Pour illustrer cette approche, un barreau en acier (E = 200 GPA, ν = 0.3) de dimension 7*18*7 mm³ sous un chargement de traction (σ.n = 10 MPa) est étudié. La fissure est dans le plan de symétrie de longueur a = 1.75 mm, le problème a donc deux plans de symétrie. La configuration ne sollicite que le mode I et son évolution doit être symétrique par rapport au milieu de l'éprouvette. La discrétisation retenue est une grille grossière d'éléments cubiques à 8 nœuds d'arête $dh_g = 0.5$ mm. A cette grille s'ajoutent trois grilles. La grille la plus fine est de dimension : $dh_f = dh_g/2^3 = 62.5$ µm, ces grilles apparaissent figure 4(b).

<figure><img src="image_5.png" /><figcaption>Fig. 4 – (a) Représentation du problème multi-grille traité : géométrie et CL, (b) maillage déformé avec champs de déplacement autour de la fissure.</figcaption></figure>

En utilisant l'erreur de stagnation des multi-grilles : $\eta = |\frac{\Delta uK\Delta u}{uKu}|$, où $\Delta u$ représente l'incrément de déplacement de la grille fine. Figure 5(a), il apparaît que la méthode converge, et que toutes les grilles convergent à la même vitesse. Cette erreur de stagnation représente le ratio de l'énergie de déformation de la correction d'une grille fine $\Delta uK\Delta u$ par rapport à l'énergie de déformation de sa grille grossière $uKu$.

<figure><img src="image_6.png" /><figcaption>Fig. 5 – (a) Erreur de stagnation des multi-grille $\eta$, (b) évolutions de $K_I$ et $T$ normalisés comparées aux valeurs obtenues en déformation planes.</figcaption></figure>

Les évaluations des coefficients mode I obtenues pour cette simulation sont représentées figure 5(b), les résultats sont bien symétriques. Le facteur d'intensité des contraintes mode I comme le T-stress ont tendance à freiner la propagation au niveau des bords libres. En effet, $K_I$ diminue, traduisant une singularité moins marquée, et la contrainte T s'accentue, traduisant l'augmentation de la triaxialité des contraintes [13]. A cœur, ces valeurs ont une amplitude comparable à celles obtenues en déformation
plane. Ici la fissure est relativement courte (i.e.a/W = 1/4). Comme dans [13] la contrainte T obtenue en trois dimension est légèrement plus faible que celle obtenue en déformation plane. On retrouve également que mode I est prépondérant : $\frac{K_{II}}{K_I} < 1\%$ et $\frac{K_{III}}{K_I} < 1\%$.

# 5 Conclusion

La méthode DEK-FEM présentée permet d'obtenir directement les coefficients asymptotiques de la mécanique de la rupture lors de la simulation d'un milieu fissuré. Les facteurs d'intensité des contraintes et la contrainte T sont obtenus. L'approche élastique permet de traiter des cas de fissuration fragile. Comme l'approche X-FEM, cette stratégie ne nécessite pas de remaillage, aussi elle parait à même de traiter des problèmes de propagation en trois dimensions, avec des lois de propagations multi-paramètres. De plus, la distribution des coefficients asymptotiques est estimée naturellement à partir de données surfaciques sur l'interface $\Gamma_W$ et de la vérification de l'équilibre dans le domaine $\Omega_W$. Enfin le couplage avec le cadre multi-grille est prometteur quand à l'aptitude de la méthode à traiter efficacement des problèmes complexes, notamment présentant des fissures au rayon de courbure variable. Pour ce type de fissure, la discrétisation du front et la finesse du maillage pourront évoluer conjointement avec la courbure pour évaluer au mieux la singularité introduite par la fissure.

Ce travail a été effectué dans le cadre du projet ANR RUPXCUBE (ANR-09-BLAN-0009-01) « Identification tridimensionnelle des lois locales de propagation des fissures par micro-tomographie X, mesures de champ étendus et simulations par éléments finis étendus ». Dans le cadre de ce projet, des essais de fissuration de fatigue sont suivis par tomographie X et également traités par une approche asymptotique. La confrontation de ces essais à la modélisation présentée ici nous permettra de valider cette approche asymptotique.

# Références

[1] T. Apel, D. Leguillon, C. Pester, and Y. Yosibash, *Edge singularities and structure of the 3-D williams expansion*, C. R. Mecanique **336** (2008), 629-635.

[2] I. Babuška and J. Melenk, *The partition of unity method*, International Journal for Numerical Methods in Engineering **40** (1997), 727-758.

[3] R. S. Barsoum, *Application of quadratic isoparametric finite elements in linear fracture mechanics*, International Journal of Fracture **10** (1974), 603-605.

[4] F. B. Belgacem, *The mortar finite element method with lagrange multipliers*, Numerische Mathematik **84** (1999), 173-197.

[5] H. D. Bui, *Mécanique de la rupture fragile*, Masson, 1978.

[6] L. N. Jr. Gifford and P. D. Hilton, *Stress intensity factors by enriched finite elements*, Engineering Fracture Mechanics **10** (1978), no. 3, 485-496.

[7] V. Lazarus, F.-G. Buchholz, and M. Fulland, *Comparison of predictions by mode II or mode III criteria on crack front twisting in three or four point bending experiments*, International Journal of Fracture **153** (2008), 141-151.

[8] J.B. Leblond, *Mécanique de la rupture fragile et ductile*, Hermes, 2003.

[9] X. Y. Liu, Q. Z. Xiao, and B. L. Karihaloo, *X-FEM for direct evaluation of mixed mode SIFs in homogeneous and bi-materials*, International Journal for Numerical Methods in Engineering **59** (2004), no. 8, 1103-1118.

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

[11] B. Moran and C.F. Shih, *Crack tip and associated domain integrals from momentum and energy balance*, Engineering Fracture Mechanics **27** (1987), no. 6, 615-642.

[12] N.I. Muskhelishvili, *Some basic problems of the mathematical theory of elasticity*, Journal of the Franklin Institute **258** (1954), 170-, Second edition, translated from the Russian by J.R.M. Radok, 704 pages. Groningen (Holland), P. Noordhoff Ltd., 1953.

[13] T. Nakamura and D.M. Parks, *Determination of elastic T-stress along 3-D crack fronts using an interaction integral*, International Journal of Solids and Structures **29** (1992), 1597-1611.

[14] S. Osher and J. Sethian, *Fronts propagating with curvature dependent speed: algorithms based on the hamilton-jacobi formulations*, J. Comput. Phys (1988), no. 79.
[15] J. C. Passieux, A. Gravouil, J. Réthoré, and M. C. Baietto, *Direct estimation of generalized stress intensity factors using a three-scale concurrent multigrid X-FEM*, International Journal for Numerical Methods in Engineering **85** (2010), no. 13, 1648–1666.

[16] J. Rannou, *Prise en compte d'effets d'échelle en mécanique de la rupture tridimensionnelle par une approche X-FEM multigrille localisée non-linéaire*, Ph.D. thesis, Institut National des Sciences Appliquées de Lyon, 2008.

[17] M. Stern, E. B. Becker, and R. S. Dunham, *A contour integral computation of mixed-mode stress intensity factors*, International Journal of Fracture **12** (1976), 359–368.

[18] M. Williams, *On the stress distribution at the base of a stationary crack*, ASME Journal of Applied Mechanics **24** (1957), 109–114.

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