# Nouveaux développements pour la propagation de fissures en fatigue de contact

Rachelle Ribeaucourt, Marie-Christine Baietto, Anthony Gravouil

► To cite this version:

Rachelle Ribeaucourt, Marie-Christine Baietto, Anthony Gravouil. Nouveaux développements pour la propagation de fissures en fatigue de contact. 8e Colloque national en calcul des structures, CSMA, May 2007, Giens, France. hal-01488183

HAL Id: hal-01488183

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

Submitted on 13 Mar 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
# Nouveaux développements pour la propagation de fissures en fatigue de contact

Rachelle Ribeaucourt* – Marie-Christine Baietto-Dubourg** –, Anthony Gravouil**

* Laboratoire LMT (ENS - UMR 8535)
Secteur Structure
61 Avenue du Président Wilson
94235 Cachan

ribeaucourt@lmt.ens-cachan.fr

** Laboratoire LaMCoS (INSA-UMR 5514)
Bât Jean d'Alembert
18-20 rue des Sciences
69621 Villeurbanne Cedex

{marie-christine.baietto, anthony.gravouil}@insa-lyon.fr

RÉSUMÉ. Cet article expose de nouveaux développements effectués dans le cadre du modèle de propagation de fissures en fatigue de contact, les prémices du modèle étant été présenté au 7ème CNCS (Giens-2005). Les nouveaux développements portent sur l'étape locale de contact intégrant l'histoire du chargement, un critère de convergence adapté aux échelles considérées, le post-traitement intégrant les spécificités du contact avec frottement le long des faces de la fissure et l'intégration de critères adaptés pour la propagation de fissures (SIF, J) en fatigue de contact. Quelques exemples de simulation sont illustrés.

ABSTRACT. This article gives new developments carried out within the model of contact fatigue crack propagation, the first steps of the model being presented at the 7th CNCS (Giens-2005). The new developments relate to the local stage of contact integrating the history of the loading, a convergence criterion adapted on scales considered, postprocessing integrating specificities of the contact with friction along the crack faces and the implementation of criteria adapted for the propagation of cracks (SIF, J) in contact fatigue. Some examples of simulation are illustrated.

MOTS-CLÉS : contact, frottement, fissure, fatigue, roulement, X-FEM, propagation.

KEYWORDS: contact, frottement, fissure, fatigue, roulement, X-FEM, propagation.
# 1. Introduction

Dans le cas de chargement de type tribologique (fatigue de roulement et fretting), les champs de contraintes et de déformations induits sont fortement localisés, multiaxiaux, non proportionnels avec de sévères gradients. Les fissures de fatigue se développant sous ces chargements sont notamment soumises à des conditions cycliques complexes alternant des phases de compression, traction, cisaillement. Leur comportement et leur propagation au cours de cycles de chargement sont donc fortement dépendantes (i) des conditions de contact et de frottement à leur interface ainsi que (ii) de l'évolution de ces conditions.

La démarche considérée pour simuler le comportement d'une fissure subissant un chargement de fatigue de contact, variable au cours du temps, a été abordée au précédent Colloque National en Calcul de Structures (7ème CNCS-2005). L'article associé exposait notamment les prémices du modèle numérique, basé sur la mécanique linéaire élastique de la rupture, prenant en compte le contact avec frottement le long des faces de la fissure où la discrétisation élément fini utilisée est la méthode X-FEM. Le processus itératif pour la détermination des conditions de contact était issu de la méthode LATIN et la loi de frottement interfaciale utilisée était une loi incrémentale de type Coulomb.

De nouveaux développements nous ont permis d'améliorer la qualité de la solution, par comparaison avec un modèle semi-analytique de propagation de fissures bidimensionnelles de fatigue de contact [DUB 02]. Cette amélioration est une conséquence (i) d'une étape locale appropriée au problème considéré, utilisant une loi incrémentale de type Coulomb, (ii) de l'implémentation d'un nouveau critère d'arrêt, (iii) ainsi que la prise en compte du contact avec frottement dans les techniques de calcul du post-traitement, permettant notamment l'obtention des facteurs d'intensité des contraintes (FIC). Les quantités de comparaison sont les champs de déplacements et les FIC. Au niveau des champs de déplacements, la comparaison permet une évaluation de l'erreur de discrétisation du modèle développé, et au niveau des FIC, la comparaison permet d'apprécier de la prise en compte de l'histoire du chargement et de son implication sur le cycle de fatigue de contact. Afin d'aborder ces trois aspects, nous allons exposer le modèle numérique développé en présentant le problème de référence, le modèle approché issu de la modélisation, suivi du modèle discret.

# 2. Description détaillée du modèle numérique

## 2.1. Le problème de référence

Le problème que nous avons à résoudre est celui du contact avec frottement le long des faces d'une fissure appartenant à un corps soumis à des sollicitations de contact cycliques. Nous supposons un comportement élastique linéaire, homogène et isotrope pour le solide, et l'hypothèse des petites perturbations. D'autre part, nous négligeons les forces d'inertie et de gravitation par rapport aux autres actions en présence. Fi-
nalement, le problème de contact surfacique est découplé du problème de structure
fissurée ce qui conduit à considérer les sollicitations de contact comme connues au
cours du temps.

Pour ce qui est de la fissure, elle est soumise à un chargement de contact tel que
l'aire de contact surfacique est de l'ordre du mm au cm, avec une pression normale
et des contraintes sous le contact de l'ordre du MPa au GPa. Ce chargement induit
de forts gradients en sous surface et dans une zone localisée, avec la présence de
champs de contraintes et déformations multiaxiaux et non proportionnels (rapport des
contraintes principales non constant au cours du temps). Soumise à des sollicitations
cycliques de traction/compression couplées à du cisaillement, des zones d'ouverture et
de fermeture apparaissent alors le long des faces de la fissure. Nous introduisons alors
μ_C, le coefficient de frottement le long des faces de la fissure. Nous avons supposé
que μ_C est constant en temps et en espace. Concernant l'aspect dissipatif du phénom-
mène de frottement, il nécessite de prendre en compte l'histoire du chargement le long
des faces de la fissure, conduisant à une stratégie incrémentale et à une modélisation
quasi-statique du problème au niveau global. Le phénomène de fatigue de contact est
alors représenté par une répétition de cycles de chargement, et chaque cycle est défini
comme une succession de pas de chargement. Par conséquent, nous utiliserons la no-
tion d'incrément du saut de déplacement tangentiel pour les conditions de frottement.

La mise en équation du problème mécanique pour un pas de chargement, $t_n$, donné et sous les hypothèses établies, considère un solide noté $\Omega$, soumis à des forces extérieures $\bar{F}$ imposées sur $\Gamma_F$ et à des déplacements imposés $\bar{U}_d$ sur $\Gamma_U$. On définit $\Gamma_C$ la géométrie de la fissure, définie par les deux surfaces de discontinuité $\Gamma_C^+$ et $\Gamma_C^-$, et $\vec{n}$ la normale extérieure à $\Omega$. Les notations $u$ et $\sigma$ seront utilisées pour qualifier les champs vectoriels du déplacement et des contraintes dans le massif, et $\omega$ et $t$ ceux de la fissure. Les équations continues et les conditions aux limites du problème de structure sont définies par les expressions suivantes :

$$ \operatorname{div} \sigma(x, t) = 0 \qquad \forall (x, t) \in \Omega \times [0, T], \quad (1) $$

$$ \sigma(x, t) = K\epsilon(x, t) \qquad \forall (x, t) \in \Omega \times [0, T], \quad (2) $$

$$ u(x, t) = \bar{U}(x) \qquad \forall (x, t) \in \Gamma_U \times [0, T], \quad (3) $$

$$ \sigma(x, t) \cdot \vec{n} = \bar{F}(x, t) \qquad \forall (x, t) \in \Gamma_F \times [0, T], \quad (4) $$

$K$ étant le tenseur élastique de Hooke. Sur $\Gamma_C^+$ et $\Gamma_C^-$ nous définissons le saut de déplacement normal (N) et tangentiel (T) par l'écriture

$$ [\omega_N \text{ ou } T] = \bar{\omega}_N \text{ ou } T - \omega_N^+ \text{ ou } T, \quad (5) $$

et l'incrément du saut de déplacement tangentiel (ou incrément de glissement) par
l'écriture $\Delta [\omega_T]$ entre le pas de chargement $t_n$ considéré et le précédent :

$$ \Delta [\omega_T] = [\omega_T]^n - [\omega_T]^{n-1} = (\omega_T^{-n} - \omega_T^{+n}) - (\omega_T^{-(n-1)} - \omega_T^{+(n-1)}). \quad (6) $$
Le contact interfacial sur $\Gamma_C$ est alors modélisé par
- les conditions de contact unilatéral

$$ [w_N] \quad := \quad \omega_N^{-} - \omega_N^{+} \quad \geq 0, \qquad (7) $$

$$ F_N \quad := \quad t_N^{+} = -t_N^{-} \quad \leq 0, \qquad (8) $$

$$ F_T \quad := \quad t_T^{+} = -t_T^{-}, \qquad (9) $$

$$ [w_N] \cdot F_N \quad = \quad 0, \qquad (10) $$

- les conditions de frottement

$$ \| F_T \| < \mu_C |F_N| \Rightarrow \quad \Delta[w_T] = \quad 0, \qquad (11) $$

$$ \| F_T \| = \mu_C |F_N| \Rightarrow \quad \exists \lambda \geq 0, \Delta[w_T] = \quad \lambda F_T. \qquad (12) $$

Les notations $F_N$ et $F_T$ désignent respectivement les composantes normale et tangentielle des forces de contact. A partir de ces équations, nous définissons trois problèmes. Tout d'abord le problème global défini par les équations (1) - (4) et qui fournit une solution au problème d'équilibre. Ensuite le problème local normal qui permet de déterminer les zones d'ouverture le long des faces de la fissure (équations (7) - (10)). Puis le problème local tangentiel qui fournit le lieu des discontinuités du champ de déplacement (i.e. le glissement) (équations (11) - (12)). Pour résoudre ces problèmes, nous faisons appel à une méthode de discrétisation de type élément fini.

## 2.2. Le problème discret

Le modèle développé est basé sur la méthode de discrétisation des éléments finis étendus, connue sous le sigle X-FEM (pour eXtended Finite Element Method). Cette méthode est appliquée aux problèmes de mécanique de la rupture depuis 1999 [MOË 99]. C'est une extension de la méthode des éléments finis qui permet de modéliser la présence d'un défaut (fissure ou autre) sans le mailler explicitement. Seul le maillage de la structure sans défaut (non fissurée) est nécessaire. Ensuite, la méthode de partition de l'unité introduite par Melenk et Babuska [MEL 96] est utilisée afin d'enrichir l'approximation élément fini standard par des fonctions additionnelles. Un choix pertinent de l'enrichissement permet de libérer le maillage de son obligation de se conformer aux discontinuités. Cette approche permet de représenter un champ de déplacement discontinu de part et d'autre d'une fissure et de prendre en compte l'aspect asymptotique du champ des contraintes à la pointe.

Cette méthode a plusieurs avantages. D'abord, elle est générale et non liée à un comportement particulier du matériau (élastique, plastique, ...). Ensuite, l'enrichissement reste valable quelle que soit la dimension du problème et quel que soit le type des éléments utilisés. Enfin, la croissance de fissure ne nécessite aucun remaillage, seule la géométrie de la fissure est nécessaire pour définir les éléments à enrichir. En quelque sorte, X-FEM étend les possibilités de la méthode éléments finis sans perdre les avantages.
<figure><img src="image_2.png" /><figcaption>Figure 1. (a) Description des ensembles $E_I$, $E_J$, $E_K$ de l'approximation du champ de déplacement élément fini enrichi (b) Description de la fonction $H$</figcaption></figure>

Concrètement, la présence de fissure dans un corps est associée à : (i) un champ de déplacement discontinu le long de la fissure et (ii) un champ de contraintes présentant une singularité en $\frac{1}{\sqrt{r}}$ en pointe de fissure. Par conséquent, la méthode X-FEM enrichit l'approximation standard élément fini par des fonctions additionnelles, (i) dis-continues pour modéliser les discontinuités du champ de déplacement, et (ii) d'autres pour caractériser le champ asymptotique en pointe de fissure. Le champ de déplacement est donc considéré comme la somme de deux champs : un champ continu et un discontinu, où l'enrichissement X-FEM se fait localement près de la fissure. Dans le cadre d'une fissure débouchante (à un front), l'approximation X-FEM du champ de déplacement est la suivante :

$$ u^h(x) = \underbrace{\sum_{i \in E_I} u_i \phi_i(x)}_{1} + \underbrace{\sum_{j \in E_J} b_j \phi_j(x) H(x)}_{2} + \underbrace{\sum_{k \in E_K} [\phi_k(x) (\sum_{l=1}^4 c_k^l F_l(r, \theta))]}_{3} \quad (13) $$

Le terme 1 correspond à l'approximation élément fini standard avec $E_I$ l'ensemble des noeuds du maillage, $u_i$ et $\phi_i$ respectivement les degrés de liberté et la fonction de forme associés au noeud $i$. Le terme 2 est un enrichissement lié à la discontinuité de la fissure, avec $E_J$ un sous-ensemble de noeuds de $E_I$, comme indiqué sur la figure 1(a), $b_j$ les degrés de liberté supplémentaires associés à l'utilisation de la fonction discontinue $H(x)$, fonction Heaviside généralisée. Les noeuds de $E_J$ sont tels que leur support intersecte la fissure mais sans contenir le front. En considérant $\Gamma_C$ comme une courbe paramétrique, pour chaque point $x$ du domaine nous pouvons trouver sur $\Gamma_C$ le point le plus proche $x^*$ de $x$. En ce point, les vecteurs normal et tangentiel à la courbe sont construits, $e_t$ and $e_n$, pris tels que $e_t \times e_n = e_z$, $e_z$ pointant vers l'extérieur de la page comme indiqué sur la figure 1(b). La fonction $H(x)$ est alors donnée par l'expression suivante

$$ H(x) = \begin{cases} 1 & \text{pour } (x - x^*) \cdot e_n > 0 \\ -1 & \text{pour } (x - x^*) \cdot e_n < 0 \end{cases} \quad (14) $$
Le terme 3 traduit la singularité des contraintes le long du front. $E_K$ est le sous-ensemble de noeuds de $E_I$ dont le support contient ce front, comme indiqué sur la figure 1(a), et $c_k^l$, $l,k=1,...,4$ les degrés de liberté correspondants. Les fonctions $F_l(r, \theta)$, $l=1..4$, sont des fonctions de forme supplémentaires avec pour expression

$$ \{F_l(r, \theta)\} \equiv \sqrt{r} \{\sin \frac{\theta}{2}, \cos \frac{\theta}{2}, \sin \frac{\theta}{2} \sin \theta, \cos \frac{\theta}{2} \sin \theta\} \quad (15) $$

où $(r, \theta)$ sont les coordonnées polaires dans un repère local lié à la pointe de la fissure.

Dans le cadre de cette méthode, la fissure n'étant pas maillée explicitement, des éléments dits "éléments d'interface" ont été introduits afin de gérer et de résoudre le contact le long des faces de la fissure. Pour tout élément fini coupé par la fissure, deux éléments d'interface en vis-à-vis $\Gamma^+$ et $\Gamma^-$ sont définis avec deux points de Gauss $(G_i^+, G_i^-)_{i=1,2}$ (voir la figure 2). Les inconnues en ces points, $\omega_i^+$ et $t_i^+$ pour $G_i^+$ ($\omega_i^-$ et $t_i^-$ pour $G_i^-$), sont exprimées dans le repère local de l'élément d'interface $\Gamma^+$ et les conditions de contact sont exprimées selon une loi d'interface pour deux points de Gauss en vis-à-vis. Cette loi, notée $\mathcal{G}$, reflète le type de contact envisagé (avec ou sans frottement, loi de frottement envisagée). Elle s'exprime par la relation

$$ \mathcal{G}(\omega_i^+, \omega_i^-, t_i^+, t_i^-) = 0 \quad \text{sur} \quad \Gamma^+ \cup \Gamma^- \quad (16) $$

Des développements nouveaux seront présentés concernant cette loi d'interface. La présence de ces "éléments d'interface" introduit des conditions aux limites supplémentaires assurant la continuité le long de la fissure :

$$ (\sigma \cdot \vec{n}_C)|_{\Gamma^+} = t^+, \quad (\sigma \cdot \vec{n}_C)|_{\Gamma^-} = t^-, \quad u|_{\Gamma^+} = \omega^+, \quad u|_{\Gamma^-} = \omega^- \quad (17) $$

<figure><img src="image_3.png" /><figcaption>Figure 2. Description des éléments d'interface $\Gamma^+$ et $\Gamma^-$, des points de Gauss $G_{1,2}^+$, $G_{1,2}^-$ et du vecteur $\vec{n}_C$</figcaption></figure>
## 2.3. La résolution du modèle discret : le schéma itératif

La stratégie de résolution utilise un schéma itératif qui est une situation dégénérée de la méthode LATIN conformément aux travaux de Champaney [CHA 96]. Cette stratégie (voir figure 3) consiste à :

- premièrement, diviser les équations du problème en deux sous-ensembles :

- un groupe d'équations locales, généralement non linéaires, définissant l'ensemble $L$, et

- un groupe d'équations linéaires, généralement globales, définissant l'ensemble $G$,

- deuxièmement, trouver la solution approchée en suivant une succession d'itérations dont le processus se fait en deux étapes. En effet, pour chaque iteration $i + 1$, la solution (indicée $s_{i+\frac{1}{2}}^L$) de l'ensemble $L$ est déterminée durant l'étape locale, puis en-suite, la solution (indicée $s_{i+1}^G$) de l'ensemble $G$ est déterminée durant l'étape globale.

Par souci de clarté nous omettrons les exposants $L$ et $G$, sachant que les quantités de chaque espace sont identifiables par les indices $i$ ou $i+\frac{1}{2}$. D'après le bilan des équations du problème et la décomposition que nous en avons fait, $L$ est constitué des équations (7)-(12) tandis que $G$ est constitué des équations (1)-(4).

<figure><img src="image_4.png" /><figcaption>Figure 3. Illustration géométrique du schéma itératif issu de la méthode LATIN</figcaption></figure>

Le nouvel itéré $s_{i+\frac{1}{2}}$ est construit à partir de $s_i$. Il appartient à la variété $L$ et est défini par la direction de montée $E^+$. Quant à l'itéré $s_{i+1}$, il est construit à partir de $s_{i+\frac{1}{2}}$ et appartient à la variété $G$. Il est obtenu par la direction de descente $E^-$.
Les espaces $E^+$ et $E^-$ sont caractérisables par des opérateurs, $H^+$ et $H^-$ respectivement. Pour que la première étape soit locale, l'opérateur $H^+$ doit l'être aussi. D'autre part, l'opérateur $H^-$ est choisi linéaire afin que l'étape globale soit linéaire. Par suite, la direction de recherche $E^+$ adoptée pour la résolution est associée à un opérateur constant linéaire $K_0$. Ainsi, les équations additionnelles définies par les opérateurs $E^+$ et $E^-$ n'impliquent que les quantités locales, $\omega$ et $t$, par les les équations de recherche :

$$ 
\begin{aligned}
t_{i+\frac{1}{2}} - t_i &= K_0(\omega_{i+\frac{1}{2}} - \omega_i) \\
t_{i+1} - t_{i+\frac{1}{2}} &= -K_0(\omega_{i+1} - \omega_{i+\frac{1}{2}})
\end{aligned} 
\qquad (18) 
$$

Pour assurer la convergence du schéma itératif [LAD p], l'opérateur constant linéaire $K_0$ doit être symétrique défini positif, et la loi de comportement doit être maximum monotone. Ces conditions étant respectées, l'algorithme converge vers la solution. L'opérateur $H^+$ utilisé par Dolbow et al. [DOL 01] est $K_0 = k_0 Id$, avec $Id$ la matrice identité et $k_0$ un réel positif. Ces auteurs ont étudié l'influence du paramètre $k_0$ sur la vitesse de convergence pour une situation de contact parfait, c'est à dire une situation où les quantités $\omega$ et $t$ sont telles que : $\omega^+ = \omega^-$, $t^+ = -t^-$. Ils ont également montré que la qualité de la solution n'était pas affectée par la variation de ce paramètre, mais que seul le nombre d'itérations à convergence était modifié. D'autre part, d'après les résultats de convergence énoncés dans [LAD p], la méthode converge toujours vers la solution du problème discrétisé. Pour un Lagrangien augmenté, la valeur de la pénalité influe sur la vitesse de convergence mais pas sur la convergence. De même pour la LATIN, la convergence est assurée par le paramètre $k_0$, et seule la vitesse de convergence est influencée par sa valeur. L'utilisation de la méthode LATIN pour traiter des problèmes de décomposition de domaine ou de contact statique en fatigue uniaxiale, a abouti à de bons résultats en choisissant le paramètre $k_0$ comme le rapport de la dimension de la rigidité sur la longueur caractéristique du problème. Nous utilisons donc comme opérateur $H^+$ la matrice $K_0 = k_0 Id$ avec $k_0 = \frac{E}{b}$, $E$ étant le module de Young et $b$ la longueur de la fissure choisie comme longueur caractéristique dans ce cas. $H^-$ est associé à l'opérateur $-K_0$.

Concernant la résolution de l'étape globale, nous cherchons $u_{i+1}$, $\sigma_{i+1}$, $\omega_{i+1}$, $t_{i+1}$, vérifiant

$$ \operatorname{div} \sigma_{i+1} = 0 \quad \text{dans } \Omega, \qquad (19) $$

$$ \sigma_{i+1} = \kappa \epsilon(u_{i+1}) \quad \text{dans } \Omega, \qquad (20) $$

$$ \sigma_{i+1} \cdot \vec{n} = \bar{F} \quad \text{sur } \Gamma_F, \qquad (21) $$

$$ u_{i+1} = \bar{U} \quad \text{sur } \Gamma_U, \qquad (22) $$

sous la contrainte de continuité le long de $\Gamma_C$

$$ u_{i+1}|_{\Gamma_C} = \omega_{i+1} \qquad (23) $$
et dans la direction (de recherche)

$$t_{i+1} - t_{i+\frac{1}{2}} = -\mathbf{K}_0(\omega_{i+1} - \omega_{i+\frac{1}{2}}). \quad (24)$$

A partir de la relation d'équilibre, de l'écriture de $\sigma$ en terme de déplacement, en tenant compte des équations de recherche de $E^-$, et sous le formalisme des multiplicateurs de Lagrange par l'introduction des efforts de contact $t_{i+1}$ (pour prendre en compte la contrainte de continuité dûe à la présence des éléments d'interface), la formulation variationnelle de ce problème est la suivante

$$ \left| \begin{array}{l} \text{trouver } u_{i+1}, \omega_{i+1}, t_{i+1} \text{ tel que } \forall \delta u, \delta \omega, \delta t \\ \\ a(u_{i+1}, \delta u) + \int_{\Gamma_C} \mathbf{K}_0 \omega_{i+1} \cdot \delta \omega \, d\Gamma \\ \\ - \int_{\Gamma_C} \delta t (u_{i+1} - \omega_{i+1}) \, d\Gamma - \int_{\Gamma_C} t_{i+1} \cdot \delta u \, d\Gamma + \int_{\Gamma_C} t_{i+1} \cdot \delta \omega \, d\Gamma \\ \\ = (f, \delta u) + (\bar{F}, \delta u) + \int_{\Gamma_C} (t_{i+\frac{1}{2}} + \mathbf{K}_0 \omega_{i+\frac{1}{2}}) \cdot \delta \omega \, d\Gamma \end{array} \right. $$

où $(u_{i+1}, \omega_{i+1})$ et $(\delta u, \delta \omega)$ appartiennent à l'espace fonctionnel $U_a$ définit par
$$U_a = \{(u, \omega) \in U_{ad}^2 / u|_{\Gamma_C^+} = \omega^+, u|_{\Gamma_C^-} = \omega^-\}, \text{ avec } U_{ad} = \{v \in \text{ régularité} / v = \bar{U} \text{ sur } \Gamma_U\}, \text{ C'est une formulation mixte qui fait intervenir le champ de déplacement } u \text{ et les quantités locales } \omega \text{ et } t.$$

La discrétisation éléments finis introduit l'écriture discrète des déplacements $\{U\} = \{N\}^T\{u\}$, avec $\{N\}$ la matrice d'interpolation et $\{U\}$ le vecteur des variables noyales sur tout le maillage par rapport aux fonctions de forme. Nous définissons les matrices $\mathbf{K}_{u\omega}$, $\mathbf{K}_{ut}$, $\mathbf{K}_{\omega\omega}$, $\mathbf{K}_{\omega t}$, permettant de lier les degrés de liberté du maillage (respectivement $u$ et $\omega$, $u$ et $t$, $\omega$ et $\omega$ et $t$). Ces écritures nous permettent d'aboutir à l'équation discrétisée du problème sous sa forme matricielle :

$$ \begin{bmatrix} \mathbf{K} & 0 & -\mathbf{K}_{ut} \\ 0 & \mathbf{K}_{\omega\omega} & \mathbf{K}_{\omega t} \\ -\mathbf{K}_{ut}^T & \mathbf{K}_{\omega t}^T & 0 \end{bmatrix} \begin{pmatrix} \{u_{i+1}\} \\ \{\omega_{i+1}\} \\ \{t_{i+1}\} \end{pmatrix} = \begin{pmatrix} \{F\} \\ \mathbf{K}_{\omega t}\{t_{i+\frac{1}{2}}\} + \mathbf{K}_{\omega\omega}\{\omega_{i+\frac{1}{2}}\} \\ 0 \end{pmatrix} \quad (25) $$

avec $\{F\}$ le vecteur des forces extérieures discrétisées, $\mathbf{K}$ matrice de rigidité du problème, $\{u_{i+1}\}$ les degrés de liberté associés à l'approximation du déplacement, c'est-à-dire les $u_i$, $b_j$ et $c_k^l$ de l'équation (13).

La stratégie employée pour résoudre le problème de contact interfacial soumis à un chargement en fatigue de contact se résume donc par la décomposition du cycle en une succession de pas de chargement. Pour chaque pas de chargement, le processus itératif résout successivement une étape locale suivie d'une étape globale où le système linéaire $\mathbb{K} \cdot V_n^{(i+1)} = F_{ext,n}^{(i)}$ est établi par l'équation (25). La matrice $\mathbb{K}$ est une matrice constante incluant la matrice de rigidité et des termes additionnels de couplage entre les variables locales et globales. Cette matrice est assemblée et factorisée uniquement à la première itération de chaque pas de chargement. Le membre de
droite $F_{ext,n}^{(i)}$ est un vecteur contenant les forces extérieures, les forces de contact et un terme nul (Les forces de contact impliquent l'état de la fissure à l'étape locale.). A la fin du processus itératif, nous tenons compte de l'état de contact interfacial pour le pas de chargement suivant. A la fin d'un cycle entier de chargement, nous disposons donc d'une information sur le comportement de la fissure dépendant du chargement appliqué. Cette information nous sera très utile pour l'estimation de la propagation en fatigue, laquelle se fera à la fin de chaque cycle de chargement.

L'initialisation du schéma pour le premier pas de chargement est construite en calculant une étape globale élastique avec les conditions locales $t_0 = 0 = \omega_0$. L'initialisation du schéma pour les autres pas de chargement calcule aussi une étape globale élastique mais en tenant compte de l'état de contact de la fissure au pas précédent.

A présent, abordons les nouveaux développements qui ont permis d'améliorer la qualité de la solution.

### 3. Les nouveaux développements apportés

Les nouveaux développements qui ont permis d'améliorer la qualité de la solution, concernent l'étape locale, l'indicateur de convergence ainsi que la prise en compte du contact avec frottement dans les techniques de calcul du post-traitement.

#### 3.1. L'étape locale

L'étape locale correspond à la résolution du problème de contact unilatéral avec frottement, c'est-à-dire la détermination des zones d'ouverture, de contact adhérent ou glissant. Le comportement de la fissure dépendant du chargement appliqué, nous tenons compte de l'état de contact interfacial au cours du cycle c'est-à-dire entre chaque pas de chargement. La nouvelle formulation tient compte de l'hystéresis liée au frottement (l'effet d'histoire) en définissant une nouvelle loi d'interface. La démarche développée se résume par l'évaluation de deux opérateurs, $\tilde{C}_I$ et $\tilde{G}_I$, fonctions des données locales,

$$ \tilde{C}_I = \frac{1}{2}(\omega_{N,i}^{-} - \omega_{N,i}^{+}) - \frac{1}{2k_0}(t_{N,i}^{-} - t_{N,i}^{+}), \quad (26) $$

$$ \tilde{G}_I = \frac{1}{2}(t_{T,i}^{+} - t_{T,i}^{-}) + \frac{k_0}{2}(\Delta\omega_{T,i}^{+} - \Delta\omega_{T,i}^{-}), \quad (27) $$

et la mise à jour des données locales (cf. table 1) avec $g = \mu_C |t_{N,i+\frac{1}{2}}^+|$.

Pour le premier pas de chargement, l'étape locale est résolue par l'approche développée par Dolbow et al. [DOL 01], où $\tilde{G}_I$ diffère par l'écriture suivante

$$ G_I = \frac{1}{2} \left\{ k_0 (\omega_{T,i}^{-} - \omega_{T,i}^{+}) - (t_{T,i}^{-} - t_{T,i}^{+}) \right\} \quad (28) $$
<table><caption>**Tableau 1.** Etat du point, composantes du déplacement et des forces de contact en fonction de la valeur des indicateurs $\tilde{C}_I$ et $\tilde{G}_I$ pour les pas de chargement avec un effet d'histoire.</caption><tr><td rowspan="2">C&#771;<sub>I</sub> &gt; 0</td><td colspan="2">C&#771;<sub>I</sub> &le; 0</td></tr><tr><td>|G&#771;<sub>I</sub>| &lt; g</td><td>|G&#771;<sub>I</sub>| &ge; g</td></tr><tr><td>Ouverture</td><td>Adhérence</td><td>Glissement</td></tr><tr><td>t<sup>-</sup><sub>N,i+1/2</sub> = -t<sup>+</sup><sub>N,i+1/2</sub> = 0</td><td>t<sup>+</sup><sub>N,i+1/2</sub> = -t<sup>-</sup><sub>N,i+1/2</sub> = k<sub>0</sub>C&#771;<sub>I</sub></td><td>t<sup>+</sup><sub>T,i+1/2</sub> = -t<sup>-</sup><sub>T,i+1/2</sub> = g G&#771;<sub>I</sub>/|G&#771;<sub>I</sub>|</td></tr><tr><td>&omega;<sup>+</sup><sub>N,i+1/2</sub> = &omega;<sup>+</sup><sub>N,i</sub> - 1/k<sub>0</sub>t<sup>+</sup><sub>N,i</sub><br>&omega;<sup>-</sup><sub>N,i+1/2</sub> = &omega;<sup>-</sup><sub>N,i</sub> - 1/k<sub>0</sub>t<sup>-</sup><sub>N,i</sub></td><td>&omega;<sup>+</sup><sub>N,i+1/2</sub> = &omega;<sup>-</sup><sub>N,i+1/2</sub> = 1/2(&omega;<sup>+</sup><sub>N,i</sub> + &omega;<sup>-</sup><sub>N,i</sub>) - 1/(2k<sub>0</sub>)(t<sup>+</sup><sub>N,i</sub> + t<sup>-</sup><sub>N,i</sub>)</td><td></td></tr><tr><td>&omega;<sup>+</sup><sub>T,i+1/2</sub> = &omega;<sup>+</sup><sub>T,i</sub> - 1/k<sub>0</sub>t<sup>+</sup><sub>T,i</sub><br>&omega;<sup>-</sup><sub>T,i+1/2</sub> = &omega;<sup>-</sup><sub>T,i</sub> - 1/k<sub>0</sub>t<sup>-</sup><sub>T,i</sub></td><td>&Delta;&omega;<sup>+</sup><sub>T,i+1/2</sub> = &Delta;&omega;<sup>-</sup><sub>T,i+1/2</sub><br>= 1/2(&Delta;&omega;<sup>+</sup><sub>T,i</sub> + &Delta;&omega;<sup>-</sup><sub>T,i</sub>) - 1/(2k<sub>0</sub>)(t<sup>+</sup><sub>T,i</sub> + t<sup>-</sup><sub>T,i</sub>)</td><td>&Delta;&omega;<sup>+</sup><sub>T,i+1/2</sub> = &Delta;&omega;<sup>+</sup><sub>T,i</sub> - 1/k<sub>0</sub>(t<sup>+</sup><sub>T,i+1/2</sub> - t<sup>+</sup><sub>T,i</sub>)<br>&Delta;&omega;<sup>-</sup><sub>T,i+1/2</sub> = &Delta;&omega;<sup>-</sup><sub>T,i</sub> - 1/k<sub>0</sub>(t<sup>-</sup><sub>T,i+1/2</sub> - t<sup>-</sup><sub>T,i</sub>)</td></tr></table>

L'opérateur $\tilde{C}_I$ est homogène à des déplacements et $G_I$ et $\tilde{G}_I$ à des forces de contact, ce qui est représentatif de chaque phénomène étudié (déplacement pour réprésenter ouverture/contact, force tangentielle pour adhérence/glissement).

Nous avons opté pour une résolution explicite du problème de contact, en résolvant d'abord le problème de contact unilatéral, puis ensuite le problème tangentiel. Nous avons vérifié numériquement et théoriquement que faire quelques itérations de point fixe pour garder le caractère implicite du problème n'apportait aucune amélioration. En effet, réinjecter les valeurs en sortie de l'étape locale à nouveau dans cette étape ne change en rien l'expression des quantités locales. Ceci est propre à la méthode LATIN, la solution locale appartenant déjà à l'ensemble L, aucune autre solution ne
peut-être trouvée à l'intersection de $E^+$ et $L$. La vérification théorique peut être trouvée en annexe du manuscript [RIB 06].

A partir des quantités locales, un indicateur de convergence est utilisé pour stopper le processus itératif. Les quantités locales issues du problème de contact avec frottement nécessite le développement d'un indicateur de convergence spécifique. C'est l'objet du paragraphe suivant.

## 3.2. Le nouvel indicateur de convergence

L'indicateur de convergence, noté $\eta$, calculé après chaque résolution globale, considère les quantités locales et se définit comme étant la distance entre les approximations globale $s_{i+1}$ et locale $s_{i+\frac{1}{2}}$ par l'écriture :

$$ \eta^2 = \frac{\| s_{i+1} - s_{i+\frac{1}{2}} \|^2}{\| s_{i+1} \|^2 + \| s_{i+\frac{1}{2}} \|^2} \quad (29) $$

Dans sa définition précédente [CHA 96, DOL 01], les champs locaux sont introduits sans faire de distinction entre les composantes normale et tangentielle du saut du champ de déplacement le long des faces de la fissure, et la norme utilisée était définie par :

$$ \| s \|^2 = \int_{\Gamma_C} t \mathbf{K}_0^{(-1)} t d\Gamma + \int_{\Gamma_C} \omega \mathbf{K}_0 \omega d\Gamma \quad (30) $$

Le problème de contact avec frottement nécessite cependant le développement d'un critère d'arrêt spécifique car l'ordre de grandeur de ces deux quantités peut être très différent. Comme cela est illustré figures 5(a) à 5(c), en présence d'une zone de contact, la composante normale est théoriquement nulle (cf. courbe de référence), et son approximation numérique devient donc très faible. L'évaluation de la composante tangentielle peut, quant à elle, atteindre des ordres de grandeur non négligeables (cf. figure 5(c)). Dans un tel cas, l'évaluation de la norme de la quantité totale est représentative de celle de la quantité tangentielle. L'arrêt du processus itératif risque donc d'être dépendant uniquement du problème tangentiel, sans représenter la contribution du problème normal. Afin de considérer les phénomènes séparément, une amélioration a été apportée en définissant un nouvel indicateur

$$ \eta^2 = \max(\eta_N) + \max(\eta_T) \quad (31) $$

où les quantités $\eta_N$ et $\eta_T$ représentent chaque problème (normal/tangentiel), et formulées par les relations suivantes

$$ \eta_N = \frac{\| s_{N,i+1} - s_{N,i+\frac{1}{2}} \|_{\infty}^2}{\| s_{N,i+1} \|_{\infty}^2 + \| s_{N,i+\frac{1}{2}} \|_{\infty}^2}, \quad (32) $$

$$ \eta_T = \frac{\| s_{T,i+1} - s_{T,i+\frac{1}{2}} \|_{\infty}^2}{\| s_{T,i+1} \|_{\infty}^2 + \| s_{T,i+\frac{1}{2}} \|_{\infty}^2}. \quad (33) $$
La norme utilisée ici est la norme infinie définie par

$$ \| s \|_{\infty}^{2} = \max_{\Gamma_C} \left( \frac{1}{k_0} t^2 + k_0 \omega^2 \right) \qquad (34) $$

La quantité $\eta_N$ est définie avec les composantes normales des quantités locales, et de même pour $\eta_T$ seules les composantes tangentielles sont utilisées. Ainsi, l'indicateur de convergence est représentatif de la position du point de Gauss le long de l'interface où le calcul de $\eta$ est maximal. Le processus itératif ne se contente plus d'être globalement correct, mais au contraire, la résolution du contact est atteinte à la précision demandée sur toute la fissure. Le critère est ainsi plus "restrictif" en ce sens où la solution doit être localement satisfaisante pour tous les points de Gauss. Ceci peut entraîner pour un pas de chargement donné un plus grand nombre d'itérations. Néanmoins ceci conduit généralement à un algorithme incrémental plus robuste et un nombre global d'itérations moindres (du fait d'un cumul d'erreur plus faible au court du temps).

Pour illustrer ces propos, considérons un massif de longueur 150 mm et de largeur 100 mm, contenant une fissure de longueur 3 mm perpendiculaire par rapport à l'horizontale et positionnée à l'origine du repère comme indiqué dans la figure 4(a). Nous allons mettre en évidence par un exemple simplifié de chargement de contact que le nouvel indicateur de convergence aboutit à des résultats satisfaisants en terme de champs locaux, notamment les sauts de déplacements le long de la fissure. Le co-

<figure><img src="image_5.png" /><figcaption>Figure 4. (a) Données du cas test pour la comparaison de deux indicateurs de convergence (b) Variations de $K_{II}$ pour le cas test figure 4(a)</figcaption></figure>

efficient de frottement le long des faces de la fissure est $\mu_C = 0.1$, et un chargement normal hertzien est appliqué en surface du massif, de la gauche vers la droite, en 5 pas de chargement, avec $P_0 = 49$ MPa, $p(y) = P_0 \sqrt{1 - \frac{(y-y_C)^2}{a^2}}$, $a$ étant la zone de contact du chargement, telle que $2a = 2$ mm. Le centre du chargement $y_C$ se déplace de la gauche vers la droite, sur l'axe horizontal (OY), en partant de la valeur $y_C = -2$ mm et en incrémentant de 1 mm jusqu'à $y_C = 2$ mm. Le maillage contient 110 éléments d'interface et le système linéaire à résoudre contient 10 534 degrés de liberté. La précision demandée est de $10^{-3}$.
Les faces de la fissure sont en contact, ce qui induit $K_I$ nul sur tout le cycle. Concernant les variations en $K_{II}$, elles sont rapportées figure 4(b) pour les deux indicateurs de convergence. L'analyse des FIC ne montre aucune différence significative, ce qui n'est pas le cas des champs de déplacement normaux. En effet les figures 5(a), 5(b) illustrent des champs de déplacement normaux obtenus au quatrième pas de chargement, pour des données identiques (même maillage, même conditions de chargement et géométriques), mais avec des indicateurs de convergence différents, en comparaison avec ceux obtenus avec le modèle 2D de référence. L'axe de la fissure est représenté en abscisse, avec l'origine correspondant à la surface. En ordonnée, le saut de déplacement est considéré en mm. Le premier calcul représente les résultats obtenus pour l'indicateur de convergence "global", tandis que le second calcul représente ceux obtenus pour l'indicateur de convergence "local". Des irrégularités de contact apparaissent sur les sauts de déplacement normaux le long de la face supérieure de la fissure, caractérisant l'aspect global de l'ancien indicateur. Figure 5(c), sont illustrés les champs de déplacements tangentiels pour l'indicateur "local" et le modèle de référence; sur ces champs, aucune différence notable n'a été observée entre l'utilisation des deux indicateurs.

<figure><img src="image_6.png" /><figcaption>Figure 5. Sauts de déplacement le long de la fissure, pour le quatrième pas d'un chargement de fatigue de contact : (a) normal avec l'ancien, (b) normal avec le nouvel (c) tangentiel avec le nouvel indicateur de convergence</figcaption></figure>

L'analyse des champs montre donc des imprécisions sur le problème normal, non représentées par les variations des FIC. Par conséquent, l'introduction du nouvel indicateur de convergence permet de s'assurer localement d'une bonne qualité de la solution en terme de champ de déplacement. Ceci permet aussi d'avoir une information sur l'estimation de l'erreur de discrétisation du modèle numérique basé sur X-FEM, par la comparaison des champs avec le modèle de référence. Des études comparatives avec le modèle 2D de référence ont été menées pour différentes configurations de fissures, de chargement et discrétisation. Elles peuvent être trouvées dans [RIB 06].

L'obtention des champs de déplacement permet une analyse des FIC pour la prévision de la propagation. Dans ce but, de nouveaux développements de post-traitement nous ont permis d'améliorer l'estimation des FIC en prenant soin d'intégrer le contact
interfacial dans l'extraction des quantités d'intérêt. Ces différents aspects vont être abordés dans la section suivante.

### 3.3. Le post-traitement

Dans le cadre X-FEM, nous utilisons des techniques intégrales indépendantes du contour afin d'évaluer les FIC, ces techniques étant beaucoup plus robustes et plus précises que les techniques basées sur les définitions locales des modes de rupture. Les nouveaux développements théoriques concernent l'utilisation des intégrales de contour en présence de contact avec frottement, et les méthodes de découplage des modes de fissuration. L'aspect novateur de ce travail réside dans la prise en compte des termes liés à la présence de contact avec frottement à toutes les étapes du calcul des FIC, ainsi que dans le choix stratégique du calcul de K2 en présence de contact en pointe.

La stratégie choisie pour le calcul des FIC est une approche globale, la méthode de l'intégrale $J$, conformément aux travaux de Palmer et Rice [PAL 73]. Initialement cette intégrale est définie pour une fissure ouverte en statique sous l'hypothèse de l'élasticité, et sur un contour ouvert $\Gamma$ (cf. figure 6), ayant la propriété d'indépendante du contour. Elle s'écrit en 2D sous la forme

$$J = \int_{\Gamma} (W n_1 - \sigma_{jk} n_k \frac{\partial u_j}{\partial x_1}) ds, \qquad (35)$$

où $W$ est la densité d'énergie de déformation, $n_i$ les composantes du vecteur extérieur unitaire de $\Gamma$, $\sigma_{jk}$ et $u_j$ les composantes des champs de contraintes et de déplacement respectivement. En présence de contact interfacial avec frottement, l'intégrale $J$ n'est plus indépendante du contour et dépend de la position des points $A^+$ et $A^-$ (les extrémités du contour $\Gamma$ situées sur chaque face de la fissure au point $A$) par la transmission des valeurs des contraintes de cisaillement en ces points [PAL 73, LEE 92]. Dans ce cas, nous avons identifié une expression indépendante du contour ouvert $\Gamma$ [PAL 73] ayant pour expression

$$J_{\Gamma} - t_{1}^{+}(A)(u_{1}^{-} - u_{1}^{+})(A). \qquad (36)$$

Ici, l'indice 1 indique la composante suivant l'axe de la fissure noté $x_1$ ; $t_1$ symbolise donc la contrainte de cisaillement s'exerçant le long des faces de la fissure ($t_1^+ = \sigma_{1k}^+ n_k = \sigma_{12}^+ n_2$) et $u_1^-$ et $u_1^+$ le déplacement tangentiel sur chaque faces. La relation entre l'intégrale $J$ et les FIC $K_I$ et $K_{II}$ en mode I et II respectivement, est :

$$J = G = \frac{k}{E}(K_{I}^{2} + K_{II}^{2}), \qquad (37)$$

où $k = 1 - \nu^2$ en déformation plane, $E$ est le module de Young et $\nu$ le module de Poisson.

Tant que la fissure est fermée, $u_2$ et $K_I$ sont nuls ; la valeur de $K_{II}$ peut donc être obtenu directement à partir de l'équation (37). Finalement si le contact s'établit
<figure><img src="image_7.png" /><figcaption>Figure 6. Définition du contour Γ utilisé pour l’intégrale J</figcaption></figure>

près du front de fissure, $K_I = 0$ et $K_{II}$ est calculé à partir de $J$. Pour les autres cas, la quantité $J$ masquant le caractère local du mode de sollicitation en ne délivrant qu'une information globale du comportement, nous faisons appel à une méthode de découplage des modes de rupture pour acquérir la connaissance locale des modes de sollicitation en pointe de fissure. Nous utilisons l'intégrale d'interaction $M^{(1,2)}$ introduite par Yau *et al.* en 1980, pour des problèmes 2D [SHI 88, YAU 80], mais aussi pour des fissures planes 3D [GOS 98] ou non-planes [NAK 89, NAK 91]. Le concept sous-jacent est basé sur l'utilisation de deux états indépendants, (1) et (2), d'un corps fissuré élastique, homogène et isotrope. Les variables associées sont décrites respectivement par les exposants (1) et (2). L'état (1) correspond au champ réel qui satisfait le problème de conditions aux limites, tandis que l'état (2) est un champ auxiliaire arbitraire. La superposition des deux états conduit à un nouveau décrit par l'exposant (0). Dans le cas d'une fissure rectiligne et libre d'effort (sans frottement ni contact, hypothèses de l'intégrale $J$), l'intégrale $J$ de l'état (0) est

$$J^{(0)} = \int_{\Gamma} (W^{(0)} \delta_{1j} - (\sigma_{ij}^{(1)} + \sigma_{ij}^{(2)}) \frac{\partial(u_i^{(1)} + u_i^{(2)})}{\partial x_1}) n_j \, ds, \quad (38)$$

où $\Gamma$ correspond au contour défini figure 6 et $W^{(0)} = \frac{1}{2}(\sigma_{ij}^{(1)} + \sigma_{ij}^{(2)})(\epsilon_{ij}^{(1)} + \epsilon_{ij}^{(2)})$. Après développements et réarrangement des termes, nous obtenons

$$J^{(0)} = J^{(1)} + J^{(2)} + M^{(1,2)} \quad (39)$$

avec $J^{(1)}$ et $J^{(2)}$ les intégrales $J$ des états (1) et (2) respectivement, et $M^{(1,2)}$, appelée l'intégrale d'interaction pour les états (1) et (2), est donnée par :

$$M^{(1,2)} = \int_{\Gamma} \left[ W^{(1,2)} \delta_{1j} - \sigma_{ij}^{(1)} \frac{\partial u_i^{(2)}}{\partial x_1} - \sigma_{ij}^{(2)} \frac{\partial u_i^{(1)}}{\partial x_1} \right] n_j \, ds, \quad (40)$$
où $W^{(1,2)}$ est l'énergie de déformation d'interaction définie par $W^{(1,2)} = \frac{1}{2}(\sigma_{ij}^{(1)} \epsilon_{ij}^{(2)} + \sigma_{ij}^{(2)} \epsilon_{ij}^{(1)})$. En utilisant l'équation (37) pour l'état (0), et après développement et simplification, nous aboutissons à l'écriture de $J^{(0)}$ sous la forme :

$$J^{(0)} = J^{(1)} + J^{(2)} + \frac{2k}{E}(K_{I}^{(1)} K_{I}^{(2)} + K_{II}^{(1)} K_{II}^{(2)}). \quad (41)$$

La comparaison des deux équations conduit à la relation :

$$M^{(1,2)} = \frac{2k}{E}(K_{I}^{(1)} K_{I}^{(2)} + K_{II}^{(1)} K_{II}^{(2)}). \quad (42)$$

Par conséquent, les FIC de l'état réel (1) peuvent être déterminés en évaluant l'intégrale d'interaction $M^{(1,2)}$ pour un choix judicieux des champs asymptotiques de l'état auxiliaire (2). Nous optons pour les trois états des modes purs asymptotiques définis par Westergaard, de sorte que :

- en mode I pur, $K_{I}^{(2)} = 1$, $K_{II}^{(2)} = 0$, $K_{III}^{(2)} = 0$, détermine $K_{I}^{(1)}$ en terme d'intégrale d'interaction par la relation :

$$K_{I}^{(1)} = \frac{E}{2k} M^{(1),\text{modeI}}, \quad (43)$$

- en mode II pur, $K_{I}^{(2)} = 0$, $K_{II}^{(2)} = 1$, $K_{III}^{(2)} = 0$, détermine $K_{II}^{(1)}$ en terme d'intégrale d'interaction par la relation :

$$K_{II}^{(1)} = \frac{E}{2k} M^{(1),\text{modeII}}, \quad (44)$$

L'utilisation de cette méthode pour des problèmes de contact interfacial a été effectuée par Dolbow et al [25], dans le cadre d'un problème de contact avec frottement en statique, mais en prennant l'intégrale J sous sa forme classique. Nous avons reformulé cette utilisation en considérant l'équation (36) à la place de J, conduisant ainsi aux équations suivantes :

$$J^{(0)} = J^{(1)} + J^{(2)} + M^{(1,2)} - t_{1}^{(1)}(A^{+})[u_{1}^{(2)}](A) - t_{1}^{(2)}(A^{+})[u_{1}^{(1)}](A) \quad (45)$$

$$J^{(0)} = J^{(1)} + J^{(2)} + \frac{2k}{E}(K_{I}^{(1)} K_{I}^{(2)} + K_{II}^{(1)} K_{II}^{(2)}) - t_{1}^{(1)}(A^{+})[u_{1}^{(2)}](A) - t_{1}^{(2)}(A^{+})[u_{1}^{(1)}](A).$$

La comparaison des deux équations conduit à la même constatation que précédemment (cf. (42)).

La méthode de l'intégrale de domaine introduite par Moran et Shih en 1987 [SHI 88, YAU 80, DOL 01] est généralement utilisé en élément fini afin de disposer d'une forme d'intégrale sur le domaine intérieur du contour $\Gamma$ plutôt qu'une intégrale sur le contour lui-même. Cela permet de lisser les oscillations parasites qui affectent la
solution numérique et d'augmenter la précision de la solution estimée. Nous avons
donc utilisé ce concept dans le cadre de contact entre les faces de fissure. Après refor-
mulation en tenant compte de l'équation (36), nous avons obtenu

$$
\begin{equation}
\begin{split}
J_{\Gamma} - t_{1}(A^{+})(u_{1}(A^{-}) - u_{1}(A^{+})) ={}& \\
& - \int_{A} \left[ W \delta_{1j} - \sigma_{ij} \frac{\partial u_{i}}{\partial x_{1}} \right] \frac{\partial q}{\partial x_{j}} \, dS + t_{1}(B^{+})(u_{1}(B^{-}) - u_{1}(B^{+})).
\end{split}
\tag{46}
\end{equation}
$$

Pour l’évaluation numérique de cette intégrale, le domaine A est obtenu comme une
collection d’éléments autour du front de fissure.

A présent nous allons exposer un exemple de simulation de propagation de fissure.

4. Quelques résultats intéressants

La complexité des conditions de contact unilatérales avec frottement à l'interface de fissures est ici illustrée sur un exemple de propagation de fissure soumise à une fatigue de contact de roulement. Le comportement de la fissure en fatigue de contact est illustré figures 7(a) et 7(b). Les courbes rouge et verte correspondent aux pas de chargement initial et final respectivement. Les sauts de déplacement sur un le cycle de chargement sont issus du modèle elfe_3d intégrant les nouveaux développements, pour le déplacement d'une charge normale hertzienne incrémentée avec $\delta l = 50 \ \mu m$, et pour une fissure inclinée $\theta = 15^\circ$, de longueur 3 mm.

<figure><img src="image_8.png" /><figcaption>Figure 7. Exemple de comportement de fissure</figcaption></figure>

L'analyse des sauts de déplacement normal le long des faces de la fissure révèle la présence d'une "bulle" d'ouverture. Cette bulle se déplace de la gauche vers la droite, c'est-à-dire de la surface libre vers la pointe de la fissure, pour yC allant de -0.9 à 0.3 mm. Nous observons une différence sur la valeur du FIC en mode I à partir de l'apparition de cette "bulle" jusqu'à sa disparition. D'autre part, après l'apparition de cette
"bulle" d'ouverture, la fissure s'ouvre progressivement en surface. Ce comportement engendre le déplacement le long de la fis- sure d'un point M qui passe progressive- ment d'un état en contact à celui d'ouverture. Ce comportement peut introduire une singularité d'un autre ordre, qui doit certainement voir une influence sur celle de la pointe de la fissure. Des différences entre les deux modèles, portant sur les valeurs des FIC, ont été identifié suite à l'apparition de cette "bulle" d'ouverture, les différences de valeur des FIC provenant du post-traitement qui suppose d'autres hypothèses. En effet, le modèle de référence utilise une approche de calcul des FIC en terme de dépla- cement, alors que le modèle elfe_3d utilise une formulation énergétique impliquant les contraintes et les déplacements. Ces résultats soulèvent des questions qui nécessitent des développements et des approfondissement théoriques.

A la fin du cycle de chargement, une analyse des grandeurs caractéristiques sur tout
le cycle permet d'obtenir la direction de fissuration. Des critères de propagation spéci-
fiques basés sur les travaux d'Hourlier [HOU 82] sont utilisés pour nos problèmes de
fatigue de contact, ainsi que des lois de propagation fournies par des essais en fatigue
de roulement. A titre d'exemple d'utilisation de ces critères, la figure 7(c) illustre la
propagation d'une fissure inclinée (zoom en pointe) soumise à un chargement de fa-
tigue de contact en roulement. Nous observons une bifurcation de la fissure dans la
direction perpendiculaire au chargement, sur plusieurs pas de chargement et en forme
de zigzag. Cette forme en zigzag est due à la présence de plusieurs maxima.

Des études paramétriques sur la longueur de la fissure, la valeur du coefficient de frottement, l'intensité et la direction de l'effort tangentiel appliqué en surface du solide fissuré ont été menées, et peuvent être retrouvées dans [RIB 06].

**5. Conclusion**

Les nouveaux développements bidimensionnels effectués pour la propagation de fissures en fatigue de contact ont permis d'intégrer une gestion du contact avec frottement entre les faces de fissures. Une description incrémentale du problème non linéaire de frottement engendré par les sollicitations multi-axiales et non proportionnelles a ainsi été implémentée. Un nouvel indicateur de convergence a été proposé. Il prend en compte séparément les phénomènes ouverture/contact et glissement/adhérence permettant ainsi une précision accrue en terme de champs locaux le long des faces de la fissure. Plusieurs étapes de validation ont été effectuées entre des résultats analytiques/numériques et ceux issus du modèle (cf. [RIB 06]). L'ensemble de ces valuations a montré la robustesse du modèle et la capacité de la stratégie du modèle à restituer des comportements de fissures variés et complexes. La simulations de propagation bidimensionnelle présentée dans le cas d'un contact tribologique de roulement utilise les nouvelles implémentations touchant aux critères de bifurcation ainsi qu'aux lois de fissuration appropriées. La validation du modèle a permis de montrer l'aptitude du code elfe_3d à reproduire le comportement de fissures de fatigue allant de configurations simples à des configurations représentatives de la fatigue tribologique, certaines de ces configurations ayant été illustrées.
## 6. Bibliographie

[CHA 96] CHAMPANEY L., *Une nouvelle approche modulaire pour l'analyse d'assemblages de structures tridimensionnelles*, Cachan : Ecole Normale Supérieure de Cachan, Thèse de doctorat, 1996.

[DOL 01] DOLBOW N., MOËS N., BELYTSCHKO T., « An extended finite element method for modeling crack growth with frictional contact », *Computer Methods in Applied Mechanics and Engineering*, vol. 190, 2001, p. 6825-6846.

[DUB 02] DUBOURG M.-C., LAMACK V., « A predictive rolling contact fatigue crack growth model : onset of branching, direction, and growth- role of dry and lubricated conditions on crack patterns. », *ASME Journal of Tribology*, vol. 124, n° 4, 2002, p. 680-688.

[GOS 98] GOSZ M., DOLBOW J., MORAN B., « Domain integral formulation for stress intensity factor computation along curved three-dimensional interface », *International Journal of Solids and Structures*, vol. 35, n° 15, 1998, p. 1763-1783.

[HOU 82] HOURLIER F., D'HONDT H., TRUCHON M., PINEAU A., « Fatigue crack path behavior under polymodal fatigue », *Institut de recherche de la sidérurgie française - IRSID*, n° RE 958, 1982.

[LAD p] LADEVÈZE P., *Nonlinear computational structural mechanics*, New York : Springer, 1998, 232 p.

[LEE 92] LEE K., JINN J., EARMME Y., « Finite element analysis of a subsurface crack on the interface of a coated material under a moving compressive load », *Wear*, vol. 155, 1992, p. 117-136.

[MEL 96] MELENK J., BABUSKA I., « The partition of unity finite element method : basic theory and applications », *Computers Methods in Applied Mechanics and Engineering*, vol. 139, 1996, p. 289-314.

[MOË 99] MOËS N., DOLBOW J., BELYTSCHKO T., « A finite element method for crack growth without remeshing », *International Journal for Numerical Methods in Engineering*, vol. 46, 1999, p. 131-150.

[NAK 89] NAKAMURA T., PARKS D., « Antisymmetrical 3-d stress field near the crack front of a thin elastic plate », *International Journal of Solids and Structures*, vol. 25, n° 12, 1989, p. 1411-1426.

[NAK 91] NAKAMURA T., « Three-dimensional stress fields of elastic interface cracks », *Journal of Applied Mechanics*, vol. 58, 1991, p. 939-946.

[PAL 73] PALMER A., RICE J., « The growth of slip surfaces in the progressive failure of over consolidated clay », *Proceedings of the Royal Society of London. Series A Mathem*, vol. 332, 1973, p. 527-548.

[RIB 06] RIBEAUCOURT R., *Gestion du contact avec frottement le long des faces de fissures dans le cadre de la méthode X-FEM. Application à la fatigue tribologique*, Lyon : Institut National des Sciences Appliquées, Thèse de doctorat, 2006.

[SHI 88] SHIH C., ASARO R., « Elastic-plastic analysis of cracks on bimaterial interfaces : part I - Small scale yielding », *Journal of Applied Mechanics*, vol. 55, 1988, p. 299-316.

[YAU 80] YAU J., WANG S., CORTEN H., « A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity », *Journal of Applied Mechanics*, vol. 47, 1980, p. 335-341.