These
These
THÈSE
Présentée devant
l’Institut National des Sciences Appliquées de Lyon
pour obtenir
le GRADE DE DOCTEUR
École doctorale :
Mécanique, Énergétique, Génie Civil, Acoustique
Spécialité :
MÉCANIQUE
par
Emilien PIERRES
...
Résumé
M OTS CLÉS : Fretting, fatigue, rupture, éléments finis étendus, contact frottement,
propagation des fissures, étude expérimentale, contact deux corps
Table des matières
Introduction 1
1 Synthèse bibliographique 7
1.1 Problèmes de rupture sous chargement de fretting fatigue . . . . . . . . . 10
1.1.1 Caractérisation des problèmes de fretting . . . . . . . . . . . . . 10
1.1.2 Amorçage des fissures sous chargement de fretting-fatigue . . . . 14
1.1.3 Propagation des fissures de fretting . . . . . . . . . . . . . . . . 15
1.2 Mécanique de la rupture en présence de non linéarités de contact avec
frottement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.2.1 Problème de référence . . . . . . . . . . . . . . . . . . . . . . . 17
1.2.2 Modes de fissuration . . . . . . . . . . . . . . . . . . . . . . . . 18
1.2.3 Analyse asymptotique du problème de fissuration . . . . . . . . . 18
1.2.4 Approches énergétiques . . . . . . . . . . . . . . . . . . . . . . 21
1.2.5 Calcul des facteurs d’intensité des contraintes . . . . . . . . . . . 26
1.2.6 Critères de propagation en fatigue . . . . . . . . . . . . . . . . . 29
1.3 Simulation de la propagation des fissures en fatigue par la méthode des
éléments finis étendus . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
1.3.1 La méthode des éléments finis étendus X-FEM . . . . . . . . . . 36
1.3.2 Représentation géométrique des fissures par fonctions de niveau
dans le cadre de la X-FEM . . . . . . . . . . . . . . . . . . . . . 38
1.4 État de l’art de la modélisation des problèmes de contact dans le cadre de
la méthode des éléments finis étendus . . . . . . . . . . . . . . . . . . . 41
i
Table des matières
ii
Table des matières
Bibliographie 165
iii
Table des matières
iv
Table des figures
v
Table des figures
vi
Table des figures
vii
Table des figures
2.28 Facteurs d’intensité des contraintes en mode II et III calculés le long des
deux fronts de fissures pour µΓ = 1, 0.5 et 0. . . . . . . . . . . . . . . . . 78
2.29 Géométrie, conditions aux limites et chargement du problème. . . . . . . 80
2.30 Maillages de la structure : (a)-pour la configuration A ; (b)-pour les confi-
gurations B et C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
2.31 Discrétisations de la fissure : (a)-pour la configuration A ; (b)-pour la
configuration B ; (c)-pour la configuration C (modèle multi-échelle). . . . 81
2.32 (a) Représentation amplifiée de la solution du problème en déplacement
global u et (b) représentation du champ des efforts locaux T à l’interface
de la fissure pour le cas de référence A. . . . . . . . . . . . . . . . . . . 82
2.33 Représentation des efforts d’interface sur le plan de fissure et tracé de
la frontière entre la zone de contact et la zone d’ouverture pour chaque
configuration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
2.34 Géométrie, chargement et conditions aux limites du problème. . . . . . . 85
2.35 (a) Maillage X-FEM en volume (3072 tétraèdres) et (b) discrétisation raf-
finée de l’interface (832 points d’intégration). . . . . . . . . . . . . . . . 85
2.36 Résultats X-FEM : (a)-Représentation amplifiée du champ de déplace-
ment U(× 120) ; (b) Efforts d’interface T. . . . . . . . . . . . . . . . . . 86
2.37 Résultats ANSYS : (a)-Représentation amplifiée du champ de déplace-
ment U(× 120) ; (b) Efforts d’interface T. . . . . . . . . . . . . . . . . . 86
2.38 Exemple de de chargement sinusoïdal P de valeur moyenne σm . Chaque
cycle [Tm ; Tm+1 ] est discrétisé en n piquets de temps tn . . . . . . . . . . . 88
2.39 LEGEND. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
viii
Table des figures
3.10 Champ des efforts T entre les faces de la fissure 2 au pas de temps 1 (Q =
Qmax . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
3.11 Maillage X-FEM localement raffiné au voisinage du front de la fissure
2 et domaine d’intégration D pour le calcul des facteurs d’intensité des
contraintes : Le champ d’extension virtuel q vaut 1 en pointe de fissure et
0 aux bord du domaine D (cf. section 1.2.5). . . . . . . . . . . . . . . . . 104
3.12 Variations du facteur d’intensité de contraintes en mode I au cours du
cycle de fretting calculées avec le modèle X-FEM et le modèle de réfé-
rence pour chaque fissure. . . . . . . . . . . . . . . . . . . . . . . . . . . 106
3.13 Variation du facteur d’intensité de contraintes en mode II en pointe de
fissure 1 au cours du cycle de fretting calculée avec le modèle X-FEM et
le modèle de référence. . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
3.14 Variation du facteur d’intensité de contraintes en mode II en pointe de
fissure 2 au cours du cycle de fretting calculée avec le modèle X-FEM et
le modèle de référence. . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
3.15 Variation de l’angle de phase ΨM = arctan KKIII au cours du cycle de
fretting. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
3.16 Notations pour la simulation numérique de la propagation des fissures en
2D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
3.17 Variation du facteur d’intensité de contraintes en mode II en pointe de
fissure 1 au cours du cycle de fretting calculée avec le modèle X-FEM
pour 6 valeurs différentes de µΓ : 0.2 ; 0.4 ; 0.6 ; 0.8 ; 1 et 1.2. . . . . . . . 111
3.18 Variation du facteur d’intensité de contraintes en mode II en pointe de
fissure 2 au cours du cycle de fretting calculée avec le modèle X-FEM
pour 6 valeurs différentes de µΓ : 0.2 ; 0.4 ; 0.6 ; 0.8 ; 1 et 1.2. . . . . . . . 112
3.19 Variation du facteur d’intensité de contraintes en mode I en pointe de
fissure 1 calculé avec le modèle X-FEM au cours du cycle de fretting
pour 4 valeurs différentes de σt : -115 MPa ; 0 MPa ; 57.5 MPa et 115 MPa.113
3.20 Variation du facteur d’intensité de contraintes en mode II en pointe de
fissure 1 calculé avec le modèle X-FEM au cours du cycle de fretting
pour 4 valeurs différentes de σt : -115 MPa ; 0 MPa ; 57.5 MPa et 115 MPa.113
3.21 Variation de l’angle de phase ΨM des sollicitations en pointe de fissure 1
au cours du cycle de fretting pour 4 valeurs différentes de σt : -115 MPa ;
0 MPa ; 57.5 MPa et 115 MPa. . . . . . . . . . . . . . . . . . . . . . . . 114
3.22 Parcours de propagation de la fissure 2 pour la configuration A : µΓ = 0.2
et σt = -115 MPa (10 pas de propagation de 2 µm). . . . . . . . . . . . . 116
3.23 Parcours de propagation de la fissure 2 pour la configuration B : µΓ = 1.2
et σt = -115 MPa (10 pas de propagation de 2 µm). . . . . . . . . . . . . 116
3.24 Parcours de propagation de la fissure 2 pour la configuration C : µΓ = 1.2
et σt = 0 (20 pas de propagation de 2 µm). . . . . . . . . . . . . . . . . . 117
3.25 Parcours de propagation de la fissure 2 pour la configuration D : µΓ = 1.2
et σt = +115 MPa (20 pas de propagation de 2 µm). . . . . . . . . . . . . 118
3.26 Représentation schématique du dispositif d’essai de fretting à deux vérins. 119
ix
Table des figures
x
Table des figures
4.3 Traces de fretting pour les essais sphère/plan conduits à (a) N = 30 000
cycles, (b) N = 50 000 cycles et (c) N = 70 000 cycles. . . . . . . . . . . 137
4.4 Coupes transversales post-mortem de l’éprouvette fissurée et reconstruc-
tion de la géométrie tridimensionnelle des fissures pour un essai conduit
à 70 000 cycles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
4.5 Géométries reconstruites des fissures tridimensionnelles pour les essais
conduits à (a) N = 30 000 cycles, (b) N = 50 000 cycles et (c) N = 70 000
cycles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
4.6 Boucle d’hystérésis de fretting en glissement partiel : L’énergie dissipée
Wd par glissement correspond à l’aire de la boucle et l’énergie totale Wt =
Wd + Wl correspond à l’aire du rectangle qui la contient [FOU 95]. . . . . 139
4.7 Évolution du rapport f aux cours des cycles des essais de fretting permet-
tant de calculer le coefficient de frottement local µ, le rapport de glisse-
ment c/a et la quantité u. . . . . . . . . . . . . . . . . . . . . . . . . . . 140
4.8 Évolution de l’effort tangentiel Q mesurée au cours d’un cycle de fretting
et discrétisée en 25 pas de temps tn . . . . . . . . . . . . . . . . . . . . . 140
4.9 (a) Distributions des pressions normales p(x,y) et (b) des contraintes tan-
gentielles q(x,y,tn = 1) sur la surface de contact sphère/plan interpolées
sur le maillage surfacique X-FEM de 30 × 30 éléments carrés de 12µm de
côté. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
4.10 Distribution du risque d’amorçage d calculé à partir du critère de Dang
Van en volume au voisinage de la surface de contact sphère/plan. . . . . . 141
4.11 Maillage de l’éprouvette (46 266 tétraèdres) raffiné au voisinage de la
surface de contact. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
4.12 (a) Définition des fonctions de niveaux pour la modélisation des fissures
de fretting et (b) dicrétisation de l’interface de ces fissures : 902 paires de
points d’intégration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
4.13 Représentation amplifiée du champ de déplacement global U, des efforts
de contact T, du glissement −[WT ] et de l’ouverture [WN ] entre les faces
de la fissure pour les étapes 1, 4 et 7 du cycle de fretting. . . . . . . . . . 144
4.14 Représentation amplifiée du champ de déplacement global U, des efforts
de contact T, du glissement −[WT ] et de l’ouverture [WN ] entre les faces
de la fissure pour les étapes 10, 13 et 16 du cycle de fretting. . . . . . . . 145
4.15 Représentation amplifiée du champ de déplacement global U, des efforts
de contact T, du glissement −[WT ] et de l’ouverture [WN ] entre les faces
de la fissure pour les étapes 19, 22 et 25 du cycle de fretting. . . . . . . . 146
4.16 Facteurs d’intensité des contraintes KI , KII et KIII calculés aux pas de
temps t1 (Q = -Qmax ) et t13 (Q = Qmax ) le long des fronts des fissures
à l’aide de l’intégrale d’interaction I en 3D (les valeurs positives sont
orientées dans le sens des x positifs). . . . . . . . . . . . . . . . . . . . . 148
4.17 Traces de fretting pour des essais sphère/plan conduits à : a) 100 000, b)
250 000, c) 500 000 et d) 750 000 cycles. . . . . . . . . . . . . . . . . . 150
xi
Table des figures
4.18 Section de l’éprouvette fissurée dans le plan médian du contact pour l’es-
sai à 250 000 cycles (représentation inversée). . . . . . . . . . . . . . . . 151
4.19 Reconstitution tridimensionnelle de la géométrie des fissures à partir de
coupes métallographiques pour les essais de fretting sphère/plan conduits
à : a) 100 000, b) 250 000, c) 500 000 et d) 750 000 cycles avec σt = 140
MPa (représentations inversées : surface vers le bas). . . . . . . . . . . . 152
4.20 (a) Distributions des pressions normales p(x,y) et (b) des contraintes tan-
gentielles q(x,y,tn = 1) sur la surface de contact sphère/plan. La zone d’in-
térêt est discrétisée en 256×256 cellules carrées de 20µm de côté. . . . . 153
4.21 Distribution du risque d’amorçage d calculé à partir du critère de Dang
Van : a) en volume et b) sur la surface de contact sphère/plan. . . . . . . . 154
4.22 KI : facteurs d’intensité des contraintes en mode I calculés au pas de temps
t1 (Q = Qmax ) le long des fronts de fissure à l’aide de l’intégrale d’inter-
action I en 3D pour les essais conduits à : a) 100 000, b) 250 000, c) 500
000 et d) 750 000 cycles (Représentation inversée ; les valeurs positives
sont orientées dans le sens des x positifs). . . . . . . . . . . . . . . . . . 155
4.23 KII : facteurs d’intensité des contraintes en mode II calculés au pas de
temps t1 (Q = Qmax ) le long des fronts de fissure à l’aide de l’intégrale
d’interaction I en 3D pour les essais conduits à : a) 100 000, b) 250 000,
c) 500 000 et d) 750 000 cycles (Représentation inversée ; les valeurs
positives sont orientées dans le sens des x positifs). . . . . . . . . . . . . 156
4.24 KIII : facteurs d’intensité des contraintes en mode III calculés au pas de
temps t1 (Q = Qmax ) le long des fronts de fissure à l’aide de l’intégrale
d’interaction I en 3D pour les essais conduits à : a) 100 000, b) 250 000,
c) 500 000 et d) 750 000 cycles (Représentation inversée ; les valeurs
positives sont orientées dans le sens des x positifs). . . . . . . . . . . . . 157
4.25 Évolution du facteur d’intensité des contraintes KI en fonction de la taille
de fissure 2 (longueur dans le plan médian : l f iss ) pour le problème considéré.158
xii
Liste des tableaux
xiii
Liste des tableaux
xiv
Introduction
1
Introduction
La prédiction de la durée de vie en fatigue sous chargement de fretting peut être divisée
en 3 étapes : (1) le calcul des champs de contraintes et de déformations multiaxiales au
voisinage de l’interface entre les composants en contact au cours du cycle de chargement ;
(2) la prédiction du risque d’amorçage, des sites et de l’orientation initiale des fissures de
fretting ; (3) la modélisation du problème de fissuration et la simulation de la propagation
des fissures en fatigue.
La résolution de chaque étape peut être basée sur une technique différente. L’étape
1 est facilement réalisable à l’aide de la méthode des éléments finis en 2D ou en 3D.
Cependant, la précision élevée requise pour le calcul de l’étape 2 implique d’importants
coûts de calculs. Les méthodes semi-analytiques sont plus largement utilisées. La mé-
thode employée ici a comme point de départ les travaux de Kalker [KAL 90]. L’hypothèse
de massif semi infini est adoptée, hypothèse classique pour la modélisation du contact
sphère/plan ou cylindre/plan et par extension pour les configurations pour lesquelles les
dimensions de l’aire de contact sont faibles par rapport aux dimensions des corps et aux
rayons de courbure des surfaces en contact. Les résolutions des problèmes normal et tan-
gentiel sont découplées. Cette hypothèse est justifiée tant que les propriétés élastiques des
matériaux en contact sont peu différentes. Elle est exacte quand les matériaux sont iden-
tiques. Les temps de calcul sont très courts et une solution précise du cycle de charge à
l’interface des deux corps est obtenue à l’échelle de ce contact.
L’étape 2 est basée sur des critères de fatigue multiaxiale basés sur le concept de
plan critique tels que les critères de Ruiz [RUI 84] ou de Dang Van [DAN 93] ou sur les
concepts de contrainte, de déformation, ou d’énergie équivalente.
L’analyse du problème de fatigue (étape 3) est généralement basée sur la mécanique
élastique linéaire de la rupture. De nombreuses lois de propagation sont fondées sur des
2
Introduction
critères énergétiques globaux tels que les facteurs d’intensité des contraintes. Des cri-
tères de propagation en fatigue multiaxiale non proportionnelle permettent également de
déterminer la direction de propagation des fissures.
Les modèles basés sur la méthode des distributions continues de dislocations per-
mettent une résolution très précise du problème de contact avec frottement entre les faces
de la fissure sous sollicitations cycliques. Celui proposé par Dudoug [DUB 89, DUB 92a]
est parfaitement adapté à la problèmatique de la fatigue tribologique – pour laquelle l’hy-
pothèse de massif semi-infini, comme nous l’avons vu, est justifiée – qui demande une
description très fine des cycles de charge complexes et des gradients sévères de contrainte
et de déformation générés dans des zones confinées sous le contact, et permet de détermi-
ner précisément l’évolution des non linéarités de contact avec frottement à l’interface des
fissures. L’état de contact avec frottement interfacial est calculé à l’échelle locale. Une
formulation incrémentale permet de prendre en compte les effets d’hystérésis liés à l’his-
toire des conditions de contact avec frottement interfacial. Les temps de calcul très réduits.
La méthode permet de modéliser des problèmes cycliques complexes (i.e. de fretting, de
roulement, des séquences d’ouverture-fermeture-adhérence-glissement) sous chargement
multiaxial non proportionnel et l’interaction de plusieurs fissures [DUB 92b]. Ce modèle
semi-analytique développé dans le cadre de la mécanique élastique linéaire de la rupture
en deux dimensions est cependant difficilement extensibles aux problèmes tridimension-
nels. Toutefois, comme nous le verrons, il peut être utilisé pour valider les résultats en 2D
issus d’autre modèles.
La méthode des éléments finis (MEF) est largement utilisée en mécanique de la rup-
ture. Sa robustesse, son aptitude à représenter des géométries complexes en 3D et la pos-
sibilité de pouvoir y intégrer des lois de comportement matériaux complexes en font un
outil privilégié pour traiter une grande variété de problèmes. Elle est cependant diffici-
lement applicable aux problèmes de discontinuités mobiles : interface fluide-structure,
changement de phase ou fissuration. En effet, un maillage conforme à la géométrie de
la surface de discontinuité est requis. Dans le cas des fissures, il est donc nécessaire de
reconstruire un nouveau maillage à chaque pas de propagation. Malgré les améliorations
apportées aux mailleurs automatiques, le coût de cette opération peut rapidement deve-
nir rédhibitoire en 3D pour des géométries de fissures complexes. Au-delà de ces aspects
topologiques, le transfert d’informations lors de la procédure de remaillage fait appel à
des techniques de projection de champs qui soulèvent encore aujourd’hui de nombreuses
questions liées notamment à la vérification des équations de conservation.
De nouvelles techniques numériques autorisant une description des discontinuités in-
dépendante de la discrétisation en volume ont été récemment développées. Dans ce travail,
nous nous intéresserons tout particulièrement à la méthode des éléments finis étendus (X-
FEM) [MOë 99] qui repose sur le principe de partition de l’unité développé par Melenk et
Babuska [MEL 96]. Cette méthode permet de modéliser les discontinuités géométriques
telles que les fissures de façon implicite, c’est-à-dire indépendamment du maillage élé-
ments finis. Elle est donc bien adaptée à la simulation de la propagation des fissures. En
effet, une simple mise à jour des fonctions d’enrichissement sans remaillage du domaine
ni projection de champs est requise lorsque la fissure progresse. De plus, cette méthode
3
Introduction
possède la plupart des avantages de la MEF tels que la possibilité de modéliser des struc-
tures complexes tridimensionnelles et d’intégrer des lois de comportement non linéaires.
La méthode X-FEM associée à la technique des fonctions de niveau [MOë 02] permet
en outre de modéliser précisément des géométries de fissures complexes éventuellement
issues d’observations expérimentales, i.e. imagerie 3D par microtomographie, imagerie
par laser confocal, coupes métallographiques.
Le développement de modèles X-FEM avec contact et frottement entre les faces de
la fissure posent toutefois quelques difficultés. La discrétisation de l’interface des fis-
sures étant dépendante du maillage éléments finis, une description précise du problème
de contact requiert un raffinement local du maillage en volume. Le rapport des échelles
mises en jeux dans les problèmes de fretting fatigue est très élevé. L’opération de maillage
du domaine peut alors s’avérer très délicate pour les problèmes tridimensionnels. De plus,
au regard du problème de la structure, le maillage n’est pas optimal et entraine un sur-
coût numérique. Ce procédé est en contradiction avec le concept d’indépendance entre le
maillage et les fissures inhérent à X-FEM. Par ailleurs, plusieurs travaux ont mis en évi-
dence des problèmes d’instabilité faisant apparaître des oscillations numériques dans la
solution de contact. La stabilité des modèles proposés dans la littérature dépend fortement
de la formulation du problème de contact et de la méthode de discrétisation du problème
utilisée.
4
Introduction
5
Introduction
Cette thèse a été réalisée au Laboratoire de Mécanique des Contacts et des Structures
(LaMCoS) de l’INSA Lyon) dans le cadre d’un projet de recherche en collaboration avec
l’entreprise SKF (thèse CIFRE). Les développements du modèle X-FEM avec contact et
frottement interfacial ont été réalisés au sein du code de calcul prototype ELFE_3D du
laboratoire LaMCoS.
6
Chapitre 1
Synthèse bibliographique
Après une brève introduction sur le phénomène de fissuration par fatigue, on s’intéresse
dans la première partie de ce chapitre aux problèmes de rupture sous chargement de
fretting fatigue. Le problème de contact entre deux corps et les principaux résultats
expérimentaux et théoriques sur l’amorçage et le propagation des fissures de fretting
sont présentés. Puis, les différentes approches de la mécanique de la rupture en présence
de non-linéarité de contact avec frottement sont décrites. Les principaux critères de
propagations en fatigue sont donnés. Ensuite, la simulation numérique de la propagation
des fissures par la méthode des éléments finis étendus (X-FEM) couplée à la méthode des
fonctions de niveau est présentée. Enfin, la dernière section de ce chapitre est consacrée
à l’état de l’art de la modélisation des problèmes de contact dans le cadre de la X-FEM.
Sommaire
1.1 Problèmes de rupture sous chargement de fretting fatigue . . . . . . . . 10
1.1.1 Caractérisation des problèmes de fretting . . . . . . . . . . . . . . 10
1.1.2 Amorçage des fissures sous chargement de fretting-fatigue . . . . . 14
1.1.3 Propagation des fissures de fretting . . . . . . . . . . . . . . . . . 15
1.2 Mécanique de la rupture en présence de non linéarités de contact avec
frottement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.2.1 Problème de référence . . . . . . . . . . . . . . . . . . . . . . . . 17
1.2.2 Modes de fissuration . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.2.3 Analyse asymptotique du problème de fissuration . . . . . . . . . . 18
1.2.4 Approches énergétiques . . . . . . . . . . . . . . . . . . . . . . . 21
1.2.5 Calcul des facteurs d’intensité des contraintes . . . . . . . . . . . . 26
1.2.6 Critères de propagation en fatigue . . . . . . . . . . . . . . . . . . 29
1.3 Simulation de la propagation des fissures en fatigue par la méthode des
éléments finis étendus . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
1.3.1 La méthode des éléments finis étendus X-FEM . . . . . . . . . . . 36
1.3.2 Représentation géométrique des fissures par fonctions de niveau
dans le cadre de la X-FEM . . . . . . . . . . . . . . . . . . . . . . 38
7
1.4 État de l’art de la modélisation des problèmes de contact dans le cadre
de la méthode des éléments finis étendus . . . . . . . . . . . . . . . . . . 41
8
Introduction :
On observe expérimentalement plusieurs mécanismes de rupture des pièces mécaniques.
La rupture fragile est caractérisée par la quasi-absence de déformation plastique macrosco-
pique et par la propagation très rapide des fissures avec une faible consommation d’énergie. La
rupture suit alors des plans cristallographiques. On parle de rupture par clivage. On l’observe gé-
néralement dans le cas des céramiques, de certains métaux et polymères à basse température.
À l’inverse, dans le cas de la rupture ductile, la déformation plastique macroscopique est im-
portante. On observe généralement une grande déformation plastique avant que la fissure ne se pro-
page, liée généralement à la présence d’inclusions ou de précipités. La majeure partie de l’énergie
est alors dissipée par plasticité.
La rupture par fatigue se distingue des autres types d’endommagement car elle se produit
pour des variations de contraintes répétitives dans le temps et généralement inférieures à la limite
d’élasticité du matériau. Elle se caractérise par la propagation relativement lente des fissures sous
l’influence des sollicitations cycliques, en présence de plasticité confinée en pointe de fissure.
On constate que la rupture des composants de structure en fonctionnement normal est le plus
souvent due à la fatigue. La durée de vie des pièces sollicitées en fatigue est généralement me-
surée en nombre de cycles N. On appelle endurance la capacité de résistance à la fatigue de ces
pièces. Expérimentalement, on caractérise généralement l’endurance d’une pièce mécanique ou
d’un matériau en construisant sa courbe de Wöhler (Contrainte σ - Nombre de cycles N), par une
série d’essais de fatigue d’amplitude de chargement variable. On trace alors le nombre de cycles
à rupture en fonction de l’amplitude de la contrainte appliquée pour chaque cas. La figure 1.1
représente schématiquement l’allure de la courbe obtenue. On distingue trois domaines :
i. une zone de fatigue oligocyclique où la rupture sous forte contrainte survient après un petit
nombre de cycles (< 104 ), précédée d’une déformation plastique locale importante ;
ii. une zone de fatigue limitée pour laquelle le nombre de cycles à rupture Nr croît lorsque la
contrainte σ décroît ;
iii. une zone d’endurance illimitée sous faible contrainte, inférieure à la limite d’endurance σD ,
pour laquelle la rupture ne survient pas avant 107 cycles.
9
1. Synthèse bibliographique
Remarque 1 : Dans certains cas comme la fatigue en milieu corrosif, il n’existe pas d’asymptote
horizontale.
Dans le cas des problèmes d’endurance limitée (zone ii), le phénomène de fatigue inclut plu-
sieurs étapes allant de la modification microstructurelle du matériau à la propagation de fissures,
conduisant après un grand nombre de cyles à la rupture de la structure. L’importance relative de
chaque étape du phénomène de rupture dépend de la géométrie, du chargement et des propriétés
du matériau [LEM 04, GER 80]. La présence de défauts dans la matière ou de concentrations de
contraintes localisées en surface favorise l’amorçage précoce de fissures. La propagation de ma-
crofissures est alors le stade critique et la prédiction de la durée de vie de la structure passe par
l’étude du comportement de ces fissures. Dans ce mémoire, nous nous intéressons plus particuliè-
rement aux problèmes de rupture sous chargement de fretting fatigue.
Le contact réel entre les composants d’une structure industrielle est généralement mal dé-
fini et difficilement reproductible. Les études expérimentales et la modélisation des problèmes de
fretting sont basées sur des géométries simplifiées : contact plan/plan, contact cylindre/plan ou
contact sphère/plan. Les configurations plan/plan présentent cependant une discontinuité dans la
distribution des pressions de contact qu’il est difficile d’évaluer. Les configurations cylindre/plan
ou sphère/plan sont donc souvent préférées.
Considérons un contact statique sphère/plan chargé sous l’effet d’une force normale P. La
théorie de Hertz [HER 81] permet de déterminer la taille de la zone de contact (disque de diamètre
10
Problèmes de rupture sous chargement de fretting fatigue
2a), la distribution de pressions normales p(x, y) en surface et le mouvement de corps rigide relatif.
Hertz introduit la simplification selon laquelle les corps en contact se comportent comme des
massifs semi infinis. Cette hypothèse est justifiée lorque les dimensions de la zone de contact sont
faibles par rapport aux dimensions de chaque corps (l’influence des conditions aux limites est
négligeable) et aux rayons de courbures des surfaces en contact.
Expérimentalement, on distingue deux types de comportement en fonction de valeur de la
force tangentielle Q. Le coefficient de frottement local à l’interface des corps en contact est noté
µ. L’application d’une force tangentielle Q < µ.P va se traduire par un glissement annulaire en
11
1. Synthèse bibliographique
périphérie du contact, alors que l’aire de contact centrale est adhérente (cf. figure 1.4.a). Cattanéo
[CAT 38] puis Mindlin [MIN 49] sont les premiers à interpréter ce comportement. Ils développent
un formalisme analytique afin d’obtenir la solution en cisaillement et la répartition de zones de
glissement et d’adhérence (a et c). Le problème est résolu sans frottement par la théorie de Hertz
puis les contraintes tangentielles q(x, y) obtenues à l’aide de la loi de Coulomb sont superposées
à la solution. Ce modèle est basé sur les équations de Boussinesq-Cerruti qui fournissent une
relation intégrale entre les déplacements et les contraintes en surface. Ils négligent ainsi l’interac-
tion des problèmes normal et tangantiel dans l’hypothèse des matériaux élastiques aux propriétés
semblables définies par le constante de Dundurs β. Ils mettent ainsi en avant la différence entre
le coefficient de frottement µ et l’évolution du rapport Q/P. L’accroissement de la force tangen-
tielle favorise l’expansion de la zone de glissement et réduit la zone d’ahérence (cf. figure 1.4.b).
Lorsque la valeur de la force tangentielle est supérieure à µ.P, la zone collée disparait et conduit à
un déplacement relatif entre les deux massifs (cf. figure 1.4.c).
L’analyse du contact deux corps basée sur les méthodes analytiques est difficile dans la confi-
guration sphère/plan car elle implique une descritpion tridimensionnelle. Elle est largement dé-
taillée dans les ouvrages de Hills et al. [HIL 94] et de Johnson [JOH 85].
- la condition de grippage pour laquelle le cycle est fermé (cf. figure 1.5.a). Le déplacement
tangentiel est accommodé par les déformations de la structure. Cette condition est non dis-
sipative et peut engendrer une destruction du matériau lors de l’ouverture du contact.
- la condition de glissement partiel pour lequel le cycle a une forme elliptique (cf. figure
1.5.b). Il se caractérise par le présence d’une zone adhérente centrale et d’une zone de
micro-glissement en périphérie.
12
Problèmes de rupture sous chargement de fretting fatigue
- la condition de glissement total pour laquelle le cycle est très ouvert et présente des paliers
car la force tangentielle Q est constante et égale à µ.P lors du glissement (cf. figure 1.5.c).
L’énergie dissipée lors du glissement est représentée par l’aire du cycle. La trace de fretting
est assimilée à une zone de glissement marquée par une zone de recouvrement.
13
1. Synthèse bibliographique
F IGURE 1.6: (a) Carte de régimes de glissement (Running Conditions Fretting Map :
RCFM) et (b) carte de réponse au fretting du matériau (Material Response Fretting Map :
MRFM).
du type d’endommagement prédominant. Dans le même but, les cartes de réponse au fretting du
matériau (Material Response Fretting Maps : MRFM) pour les problèmes de fretting ont aussi été
proposées (cf. figure 1.6.b). Celles-ci sont construites à partir des traces de fretting expérimentales
[FOU 96].
On distingue donc deux types principaux de problèmes de fretting. Les problèmes de fretting
wear pour lesquels l’usure est l’endommagement prédominant et les problèmes de fretting fatigue
auxquels nous nous intéressons ici et pour lesquels la condition de glissement partiel à l’interface
des corps en contact provoque l’amorçage et la propagation de fissures.
14
Problèmes de rupture sous chargement de fretting fatigue
F IGURE 1.7: Résultats expérimentaux d’un essai de fretting cylindre/plan [MUñ 06] :(a)
Trace de fretting ; (b) Coupe transversale de l’éprouvette fissurée.
L’analyse de l’amorçage des fissures de fretting est basée sur les cartes de réponse au fretting
du matériau construites expérimentalement et les critères d’amorçages analytiques. Ces derniers
sont basés sur les critères en fatigue multiaxiale et sont adaptés aux problèmes de fatigue de contact
tels que les problèmes de fretting. On trouve notamment les critères de plan critique suivants :
- le critère de McDiarmid [MCD 91] pour lequel l’amorçage des fissures a lieu dans le plan
critique où les contraintes de cisaillement sont maximales.
- le critère de Dang Van [DAN 93, FOU 96, FOU 02, BAI 10, PIE 10b] qui établit une condi-
tion d’amorçage à grand nombre de cycles. Il permet de prédire le site et la direction d’amor-
çage des fissures à partir du champ mésoscopique des contraintes sous le contact (partie
alternée du cisaillement et pression hydrostatique) et des propriétés en fatigue du matériau.
Ce critère a notamment été utilisé dans [PIE 10c, LAM 97b].
- Le critère de Liu [LIU 05] pour lequel l’amorçage des fissures a lieu dans le plan formant
un angle α (dépendant du matériau) avec le plan critique où l’amplitude des contraintes
normales est maximale.
- le critère SWT (Smith-Watson-Topper) [SMI 70] a été développé pour des nombres de
cycles compris entre 2.104 et 106 . Ce critère de fatigue stipule que l’amorçage des fissures
se produit dans le plan qui maximise le produit de l’amplitude de déformation normale par
la contrainte normale maximale au cours d’un cycle. Ces critères ont notamment permis de
prédire précisément l’amorçage de fissures de fretting dans [SZO 98, FRI 05, PRO 06].
15
1. Synthèse bibliographique
non proportionnel et d’amplitude variable. Des essais de fretting permettent d’établir la courbe
de Wöhler du matériau et de déterminer l’abattement de la limite d’endurance σD en présence
de fretting [SZO 98, SHI 06]. Cette courbe passe par un minimum. En effet, des contraintes de
fretting très importantes conduisent à un état de glissement total et l’usure de la surface de contact
empêche les fissures de se développer [WAT 81, NAI 05].
De plus, les déformations locales induisent un état compressif variable le long des fissures,
ce qui génère des séquences d’ouverture-fermeture-glissement-adhérence [LAM 96, LAM 97a].
L’état de contact entre les faces de la fissure influence l’état de mixité et conditionne le bran-
chement et la vitesse de propagation des fissures. La modélisation numérique des problèmes de
fretting fatigue requiert donc une modélisation fine du contact et la prise en compte du frottement
entre les lèvres de la fissure. Cependant, la valeur du coefficient de frottement est difficile à évaluer
expérimentalement.
Endo et Goto (1976) [END 76], puis Alic et Hawley (1979) [ALI 79] observent expérimenta-
lement l’apparition d’un réseaux de micro fissures obliques (≈ 45◦ ) sous la zone de micro glisse-
ment de la surface de contact. Au cours des essais, une fissure principale se propage alors que les
autres fissures s’arrêtent et disparaissent avec l’usure de la surface. Ils étudient expérimentalement
la propagation en 2D de fissures planes perpendiculaires à la surface, soumises à un chargement
de fretting fatigue couplé à un chargement de fatigue conventionnelle. Il ont ainsi relié les vitesses
de propagation de ces fissures aux valeurs du facteur d’intensité des contraintes KI calculées ana-
lytiquement (cf. section 1.2.3). La comparaison de ces résultats avec ceux d’une fissure soumise
à un chargement de fatigue uniaxial a montré que les vitesses de fissuration sous chargement de
fretting fatigue dans la zone proche de la surface de contact sont très supérieures à celles observées
sous chargement de fatigue pure. Ces vitesses sont équivalentes à partir d’une certaine profondeur
variable en fonction des paramètres de chargement surfacique. Par ailleurs, ces essais ont montré
que la fissure branche perpendiculairement à la surface au delà de cette profondeur. L’influence du
chargement de fretting fatigue est alors négligeable.
Des essais sur les polymères ont permis d’observer les géométries complexes des fissures tri-
dimensionnelles de fretting fatigue [BAI 03, CHA 06c]. L’expérience montre que la durée de vie
des pièces mécaniques sollicitées en fretting est divisée en quatre stades [LAM 97b, CHA 06b,
DUB 00] (cf. figure 1.8). Le stade I correspond à l’amorçage de fissures sous l’influence du
contact. Le stade II est la croissance de certains de ces défauts en fissures macroscopiques dans
une zone confinée sous le contact par l’action combinée du chargement local de fretting et de la
contrainte volumique. Le stade III correspond à la propagation en volume d’une fissure principale
sous l’influence majoritaire d’une contrainte volumique statique ou cyclique. Enfin, le stade IV
correspond à la rupture de la pièce.
L’amorçage des fissures fretting peut se produire pour un faible nombre de cycles. La majeure
partie de la vie des pièces industrielles est donc dépendante de la propagation de ces fissures.
Le comportement de ces fissures est encore mal maîtrisé. Malgré l’importance des enjeux indus-
triels liés aux problèmes de rupture par fretting, les techniques de développement actuelles ne per-
mettent pas de prendre en compte la complexité du phénomène. Il est donc nécessaire développer
des méthodes de prédiction de durée de vie globales adaptées aux spécificités de ces problèmes,
permettant de prendre en compte les phénomènes d’amorçage et d’usure des surfaces et de si-
muler la propagation des fissures de fretting. La modélisation de leur comportement requiert une
description multi-échelle des phénomènes locaux tels que la plasticité en pointe ou le contact avec
16
Mécanique de la rupture en présence de non linéarités de contact avec frottement
F IGURE 1.8: Représentation schématique des différents stades de propagation des fis-
sures de fretting [CHA 06b].
frottement entre les faces de la fissure. Elle passe par l’étude des mécanismes de propagation en
fatigue en mode mixte et l’utilisation de critères de propagation multiaxiale dédiés aux chargement
non proportionnels. Dans la section suivante, les outils de la mécanique de la rupture permettant
cette analyse sont présentés.
17
1. Synthèse bibliographique
Γ = Γ+ ∪ Γ−
Notations :
∂Ω = ∂1 Ω ∪ ∂2 Ω ∪ Γ
Remarque 2 : Une des particularité du problème de rupture est que la frontière du domaine peut
évoluer dans le temps lorsque la fissure se propage. Toutes les variables dépendent donc du temps
ainsi que Ω, ∂Ω, ∂1 Ω, ∂2 Ω, Γ+ , Γ− et F .
Les inconnues de ce problème sont le champ de déplacement u(t), le tenseur des contraintes
volumiques σ(t), et le vecteur des efforts surfaciques de contact t(t) sur les lèvres de la fissure.
Dans la suite de ce chapitre, différentes approches sont présentées afin d’analyser le comportement
des fissures et de prédire leur propagation à partir de la solution de ce problème.
Remarque 3 : La discontinuité normale [u2 ] n’est jamais négative car les faces de la fissure ne
peuvent pas s’interpénétrer.
18
Mécanique de la rupture en présence de non linéarités de contact avec frottement
(1.1) pour les milieux bidimensionnels relativement simples et sans considérer les phénomènes
plastiques en pointe de fissure :
où les fonctions fi j sont les fonctions de Westergaard [WES 39]. Il définit ainsi les facteurs d’in-
tensité des contraintes KI , KII et KIII associés à chacun des modes de sollicitation, qui quantifient
√
l’intensité de la singularité des contraintes. Ils s’expriment en MPa m.
Au voisinage immédiat du font de fissure, les contraintes présentent une singularité en √1r .
Elles sont infinies en pointe de fissure. On peut donc difficilement établir des lois de propagation
précises basées sur leur valeur. C’est pourquoi on utilise généralement les Facteurs d’Intensités
des Contraintes (FIC). Deux écritures des FIC sont alors définies [BUI 78] :
- les facteurs d’intensité des contraintes statiques. Ils s’expriment en fonction du champ des
19
1. Synthèse bibliographique
√
KI = lim 2πrσ22 (θ = 0)
r→0
√
KII = lim 2πrσ21 (θ = 0) (1.2)
r→0
√
KIII = lim 2πrσ23 (θ = 0)
r→0
- les facteurs d’intensité des contraintes cinématiques, qui sont proportionnels à la disconti-
nuité du champ de déplacement au voisinage du front de fissure :
q
KI = lim µ 2π
[u2 (θ = π)]
r→0 k+1 r
q
µ 2π (1.3)
KII = lim k+1 r [u1 (θ = π)]
r→0
q
KIII = lim 4µ 2π
r [u3 (θ = π)]
r→0
où (r, θ) sont les coordonnées cylindriques du repère centré sur le front de fissure (cf. figure
1.10), [ui (θ = π)] = ui (θ = +π) − ui (θ = −π) est le saut de déplacement lié à la disconti-
nuité, µ = E/2(1 + ν) et k est la constante de Kolosov :
3 − 4ν en déformations planes
k= 3−ν
1+ν en contraintes planes
On peut alors introduire ces facteurs dans les solutions analytiques de Westergaard :
1 r
p
u1 (r, θ) = 2µ 2π KI cos θ2 (k − cosθ) + KII sin θ2 (k + 2 + cosθ)
1 r (1.4)
p
u2 (r, θ) = 2µ 2π KI sin θ2 (k − cosθ) − KII cos θ2 (k − 2 + cosθ)
2 r
p θ
u3 (r, θ) = µ 2π KIII sin 2
A partir de ces résultats, Irwin établit une relation entre le taux de restitution d’énergie G (sous
la forme d’une intégrale de contour) et les facteurs d’intensité des contraintes pour une fissure non
débouchante tridimensionnelle :
1 − ν2 2 1+ν 2
G= KI + KII2 + KIII (1.5)
E E
Remarque 4 : Dans le cas d’une fissure débouchante, la nature de la singularité en pointe peut
dépendre de l’angle de la fissure avec la surface libre ainsi que du matériau. De même, en présence
de contact avec frottement entre les faces de la fissure, cette singularité évolue. Nous reviendrons
sur ce point ultérieurement.
Remarque 5 : Les valeurs des contraintes tendent vers l’infini au front de fissure, indépendam-
ment de la charge, de la géométrie de la fissure ou du matériau. Une zone plastique existe donc
20
Mécanique de la rupture en présence de non linéarités de contact avec frottement
autour du front et le problème ne satisfait pas l’hypothèse d’élasticité linéaire. En particulier dans
le cadre des problèmes de fissuration par fretting, on peut faire l’hypothèse de plasticité confinée.
On suppose alors que la redistribution des contraintes liée à la plasticité ne modifie pas les champs
élastiques. Les équations (1.2) et (1.3) sont alors valides lorsque la taille de la zone plastique r p
est suffisamment petite par rapport à la taille de la zone de dominance de la singularité élastique
(zone de K-dominance), liée à la taille de la fissure a. On admet généralement cette hypothèse
pour r p /a < 3%.
En théorie, on peut donc obtenir les valeurs des facteurs d’intensité des contraintes à partir
des champs de contrainte ou de déplacement à partir des expressions (1.2) ou (1.3). Cependant,
dans la majorité des cas, les différentes méthodes numériques ou expérimentales ne permettent
pas d’obtenir une approximation précise de ces champs suffisamment proche du front de fissure.
L’identification des facteurs d’intensité des contraintes par cette approche peut alors s’avérer dif-
ficile. En revanche, les approches énergétiques sont basées sur les champs plus lointains du front
pour lesquels on a une meilleure estimation ou sur la variation d’énergie globale du problème.
Les critères de propagation énergétiques sont beaucoup plus précis et faciles à calculer. Ils sont
présentées dans la section suivante.
21
1. Synthèse bibliographique
où Wl représente l’énergie élastique. Par ailleurs, la puissance des efforts extérieurs s’écrit :
Z Z
∂u ∂u
Pext = F dS + f dV (1.9)
∂2 Ω ∂t Ω ∂t
Z Z
!
∂K ∂u ∂u ∂Welast ∂A
= F dS + f dV − − 2γS (1.11)
∂t ∂2 Ω ∂A Ω ∂A ∂A ∂t
Z Z
∂u ∂u ∂Welast
G= F dS + f dV − (1.12)
∂2 Ω ∂A Ω ∂A ∂A
La théorie de Griffith constitue à elle seule un critère de propagation de fissure. Elle peut être
utilisée notamment pour étudier le risque de propagation en mode I de fissures dans les matériaux
fragiles. Cependant, cette approche globale ne permet pas de prendre en compte les effets locaux
et de déterminer la direction de propagation des fissures dans les problèmes tridimensionnels en
mode mixte. En revanche, les intégrales de contour sont calculées à l’aide de champs locaux au
voisinage étendu du front de fissure et permettent de découpler les modes de sollicitation.
22
Mécanique de la rupture en présence de non linéarités de contact avec frottement
Intégrale J :
En 1968, à partir des travaux d’Eshelby (1956), Rice [RIC 68] et Cherepanov [CHE 67] intro-
duisent indépendamment le concept d’intégrales indépendantes du contour. Rice définit l’intégrale
J comme le taux de diminution d’énergie potentielle stockée dans un domaine fermé autour de la
pointe de la fissure lorsque cette dernière se propage :
Z
∂ui σi j εi j
J= Wl δ1 j − σi j n j ds avec Wl = (1.13)
C ∂x1 2
où C est un contour orienté autour de la pointe de la fissure, n sa normale unitaire sortante (cf.
figure 1.13) et Wl l’énergie de déformation pour un matériau élastique linéaire.
Dans [RIC 68, BUI 78], on montre que cette intégrale est indépendante du contour C si :
i. les faces de la fissure sont libres d’effort et les forces volumiques sont nulles ;
ii. la fissure est rectiligne (en 2D) ou plane (en 3D) et dans le prolongement de l’axe e2 entre
les points A et la pointe P.
L’intégrale J est la première composante J1 du vecteur des forces configurationnelles
[MAU 92]. En prenant le contour C comme étant un cercle centré sur la pointe de fissure et dont le
rayon tend vers 0, on peut montrer que J est identique au taux de restitution de l’énergie : J ≡ G.
On peut ainsi relier l’intégrale J aux facteurs d’intensité des contraintes par l’intermédiaire de la
formule d’Irwin (1.5) [RIC 68].
La présence de contact avec frottement est à l’origine d’une dissipation d’énergie entre les
faces de la fissure. L’hypothèse i n’est plus valide et l’intégrale J n’est plus indépendante du
23
1. Synthèse bibliographique
contour [PAL 73]. À partir de la formule (1.13), si l’on considère le contour fermé C1 ∪ [A+ B+ ] ∪
C2 ∪ [B− A− ] (cf. figure1.14), on peut écrire [GOS 02, DOL 01, RIB 07a] :
Z Z
∂ui ∂ui
JC1 + Wl δ1 j − σi j (−n j ) ds − JC2 + Wl δ1 j − σi j n j ds = 0 (1.14)
A+ B+ ∂x1 B− A− ∂x1
où σ12 est la contrainte de cisaillement due au frottement entre les faces de la fissure. On peut
alors écrire une intégrale indépendante du contour C (cf. figure 1.13) en présence de contact avec
frottement entre les faces de la fissure :
Z
∂ui
˜ n j ds − σ12 (A) u− +
J= Wl δ1 j − σi j 1 (A) − u1 (A) (1.18)
C ∂x1
Intégrale d’interaction :
24
Mécanique de la rupture en présence de non linéarités de contact avec frottement
Le concept d’intégrale d’interaction pour le calcul des facteurs d’intensité des contraintes
est introduit par Yau et al. en 1980 [YAU 80]. Il est fréquemment utilisé dans la littérature en
élasticité linéaire 2D [MOë 99], en 3D [NAK 92, GOS 98, RAJ 00, GRA 02, GOS 02, RAN 08b]
et en dynamique [RéT 05].
Deux états d’équilibre indépendants d’un corps élastique homogène et isotrope sont consi-
dérés : un état réel (ℜ) qui satisfait les conditions aux limites et un état auxiliaire fictif (aux)
pour lequel les lèvres de la fissure sont libres de tout effort. La superposition de ces deux états
d’équilibre conduit à un nouvel état d’équilibre pour lequel est écrite l’intégrale J˜ pour un contour
C:
∂ uℜ +uaux
(i i )
J˜ℜ+aux = aux
R ℜ+aux
C Wl δ1 j − σi j + σi j
ℜ
∂x1 n j ds
h i (1.19)
aux (A)
−,ℜ −,aux +,ℜ +,aux
− σℜ
12 (A) + σ12 u1 (A) + u1 (A) − u1 (A) + u 1 (A)
où
1 ℜ
σi j + σaux
Wlℜ+aux = ij εℜ
ij + εaux
ij
2
Le champ des contraintes auxiliaires est nul sur les lèvres de la fissure : σaux
12 (A) = 0. En
réorganisant les termes (1.19) devient :
J˜ℜ+aux = J˜ℜ+ J˜aux
∂uaux
h i
∂uℜ ℜ (A) u−,aux (A) − u+,aux (A)
Wlℜ,aux δ1 j − σℜ − σaux
R
+ C ij
i
∂x1
i
i j ∂x1 n j ds + σ12 1 1
(1.20)
où
1 ℜ aux
Wlℜ,aux = σi j εi j + σaux εℜ
ij ij
2
L’intégrale d’interaction est alors définie par :
∂uaux
Z ℜ
i aux ∂ui aux
I ℜ,aux = Wlℜ,aux δ1 j − σℜ
ij − σ ij n j ds + σℜ
12 (A) [u1 (A)] (1.21)
C ∂x 1 ∂x1
E
KIaux = 1 , KIIaux = 0 , KIII
aux = 0 =⇒ K ℜ =
I 2(1−ν2 )
I ℜ,modeI
E
KIaux = 0 , KIIaux = 0 , KIII
aux = 1 =⇒ K ℜ =
III 2(1+ν) I
ℜ,modeIII
Les méthode énergétiques présentées dans cette section permettent de calculer les grandeurs
caractéristiques de la mécanique élastique linéaire de la rupture en prenant en compte les non li-
néarités de contact avec frottement entre les faces de la fissure. L’intégrale d’interaction permet de
25
1. Synthèse bibliographique
découpler mes modes de rupture I, II et III. Dans le cadre des méthodes éléments finis, il est diffi-
cile de calculer avec suffisamment de précision les intégrales de contour. Elles sont généralement
transformées en intégrales de domaine, i.e. de surface en 2D ou de volume en 3D, afin d’obte-
nir une précision maximale. Le calcul numérique des intégrales de domaine est présenté dans la
section suivante.
où le contour C du plan (P, e1 , e2 ) contient le front de fissure (cf. figure 1.13). n sa normale exté-
rieure. Le vecteur δa(s) = δl(s)e1 est appelé l’extension virtuelle de la fissure. Si δl(s) = 1 et si
la fissure est rectiligne dans le plan considéré, on retrouve les expressions de l’intégrale J˜ et de
l’intégrale d’interaction en fonction du choix du tenseur A :
Ai j = W δi j − σl j ∂ul ˜
⇒ L(s) = J(s)
∂xi
1
aux ℜ
aux + σaux εℜ δ − σℜ ∂ul − σaux ∂ul
(1.25)
Ai j = 2 σℜ ε
kl kl kl kl i j l j ∂xi l j ∂xi ⇒ L(s) = I(s)
Lorsque la fissure est non plane près du front, on fera tendre le contour d’intégration vers
le point P afin de se rapprocher des conditions de fissure rectiligne pour lesquelles ces formules
intégrales ont été développées.
Dans le cadre des méthodes éléments finis, les intégrales de contour sont difficiles à calculer
avec suffisamment de précision. Il est alors possible de les transformer en intégrales de domaine
(de surface en 2D ou de volume en 3D) [MOR 87, SUO 92, MOë 02]. Celles-ci lissent les éven-
tuelles oscillations numériques de la solution discrétisée et permettent d’obtenir à moindre coût
une précision beaucoup plus élevée. Dans le cas tridimensionnel, une nouvelle intégrale curviligne
L̄ est définie sur une partie du front F(s) :
Z
L̄(s) = L(s) ds (1.26)
F(s)
26
Mécanique de la rupture en présence de non linéarités de contact avec frottement
et Se à l’extérieur, par Sg et Sd aux extrémités et par les lèvres de la fissure Γ+ et Γ− (cf. figure
1.16). Une nouvelle intégrale est définie sur ce domaine D :
Z
H= Al j, j ql dV (1.28)
D
27
1. Synthèse bibliographique
où le vecteur normal n à la surface Si est orienté vers le front de fissure. À partir de 1.28 et 1.30
on écrit :
Z Z Z Z
Al j, j ql dV = Al j ql n j dS + Al2 ql nΓ dS − Al j ql, j dV (1.32)
D Si Γ+ ∪Γ− D
Lorsque Si tend vers ce front, le volume D est appelé D0 et une nouvelle expression de l’inté-
grale L̄ est obtenue :
Z Z Z
L̄ = Al j, j ql dV + Al j ql, j dV − Al2 ql nΓ dS (1.33)
D0 D0 Γ+ ∪Γ−
De plus, en suivant les préconisations de [MOë 02] pour le calcul de l’intégrale d’interaction,
la construction du champ d’extension virtuel s’appuie sur une base locale conforme à la géométrie
de la fissure au voisinage du front de fissure : q = q1 · e1 (cf. figure 1.15). Ce choix conduit à :
Z Z Z
L̄ = A1 j, j q1 dV + A1 j q1, j dV − A12 q1 nΓ dS (1.34)
D0 D0 Γ+ ∪Γ−
Dans [SUO 92], on montre en 2D que cette intégrale de domaine est exactement égale à l’ex-
pression sur le contour donnée par l’équation (1.24) et qu’elle est indépendante de ce contour. En
3D, l’hypothèse d’invariance de L(s) le long de la partie du front F(s) est faite. La qualité de l’ap-
proximation numérique de la solution au voisinage du front à une forte influence sur le résultat.
Il existe donc une taille du domaine d’intégration optimale dépendante du rayon de courbure du
front de fissure (en 3D) et de la finesse locale du maillage des éléments finis [RAN 08a].
28
Mécanique de la rupture en présence de non linéarités de contact avec frottement
da
= C (∆K)m (1.35)
dN
où C et m sont des constantes intrinsèques au matériau déterminées expérimentalement.
- le régime C de propagation instable pour lequel la vitesse de propagation est supérieure à
celle prévue dans le regime de Paris. ∆K tend alors vers la ténacité du matériau K1C pour
laquelle on observe expérimentalement la rupture de la pièce.
La loi de Paris rend compte uniquement du régime de propagation stable B. Forman et al.
[FOR 67] ont alors proposé une relation plus complète permettant de tenir compte des deux parties
asymptotiques qui définissement la rupture (régime C : ∆K ≥ K1C propagation instable) et le seuil
de non fissuration (régime A : ∆K ' KS ). On trouve aussi dans la littérature un certain nombre de
lois de propagations purement phénoménologiques et dérivées de la loi de Paris prenant en compte
l’influence du rapport de charge R = σσmax
min
, de la sollicitation moyenne et de la fréquence du cycle
de chargement.
En outre, la propagation en fatigue est très fortement liée aux phénomènes non-linéaires tels
que la plasticité confinée en pointe de fissure ou le contact avec frottement entre les lèvres de
la fissure. Elber [ELB 70] a constaté que, même avec des rapports de charge positifs, la fissure
est fermée pendant une partie du cycle. Elle ne commence à s’ouvrir qu’à partir d’une valeur de
contrainte σouv > σmin appelée contrainte d’ouverture. Ainsi, l’expérience montre que la vitesse de
propagation d’une fissure n’est pas toujours proportionnelle à l’amplitude du facteur d’intensité
des contraintes appliqué ∆K mais plutôt à un facteur d’intensité des contraintes effectif ∆Ke f f =
Kmax − Kouv permettant de prendre en compte les effets de refermeture en pointe de fissure liés la
29
1. Synthèse bibliographique
F IGURE 1.17: Allure générale du taux d’accroissement d’une fissure sollicitée en mode
I en fonction de l’écart de variation du facteur d’intensité des contraintes ∆K lors d’un
cycle de chargement.
plasticité confinée, la rugosité des faces, l’oxydation ou une transformation de phase induite par
le champ mécanique [ELB 71, SUR 98]. La propagation de la fissure dépend également des effets
d’histoire tels que les phénomènes de surcharge. Par exemple, une augmentation ponctuelle de la
contrainte maximum σmax provoque une diminution temporaire de la vitesse de propagation.
30
Mécanique de la rupture en présence de non linéarités de contact avec frottement
d’un chargement multiaxial non proportionnel pour lequel ces rapports varient. Premièrement, on
s’intéresse aux différents critères de direction de propagation en mode mixte proportionnel utilisés
dans la littérature :
Critère MTS (Maximum Tangential Stress) : Erdogan et Sih proposent en 1963 le critère de
rupture en mode mixte selon lequel, dans le cas d’un problème plan, la rupture d’un matériau
fragile se produit dans la direction θ pour laquelle la contrainte circonférencielle σθθ est maximale
et la contrainte de cisaillement σrθ est nulle [ERD 63]. Ce qui revient à résoudre :
∂σθθ ∂2 σθθ
= 0 et <0 (1.37)
∂θ ∂θ2
À partir des expressions asymptotiques des contraintes en pointe de fissure [WES 39], l’angle
de propagation en mode mixte proportionnel pour les problèmes 2D est :
s
2
1 KI KI
θ0 = 2 · atan · − sign (KII ) · + 8 (1.38)
4 KII KII
Critère MSS (Maximum Shear Stress) : Ce critère en cisaillement a été proposé par Otsuka et
al. [OTS 75] en complément du critère d’Erdogan et Sih pour les problèmes de propagation de
fissures en mode II pur :
∂σrθ
=0 (1.39)
∂θ
Des essais bidimensionnels en mode mixte (I+II) sur des éprouvettes en aluminium ont validé
ces critères.
Critère de Sih de densité minimale d’énergie de déformation : En 1974, Sih propose ce cri-
tère pour lequel la fissure se propage dans la direction où la densité d’énergie de déformation ou
d’énergie élastique W (εi j ) est minimale [SIH 74]. Il définit alors le facteur de densité d’énergie de
déformation S tel que :
dW
S=r (1.40)
dV
avec
dW 1 ν 1+ν 2
σ2x + σ2y + σ2z − (σx σy + σy σz + σz σx ) + σxy + σ2yz + σ2zx
= (1.41)
dV 2E E E
où r est la distance depuis le front de fissure. Sih exprime ensuite le facteur de densité d’énergie
de déformation S en fonction des facteurs d’intensité des contraintes :
31
1. Synthèse bibliographique
avec
En comparaison avec les critères dits explicites présentés ci-dessus, les critères suivants sont
qualifiés d’implicites car ils sont basés sur des quantités calculées à l’extrémité d’une extension
virtuelle de la fissure de longueur s et de direction θ (cf. figure 1.18) :
Critère de taux de restitution d’énergie maximale Gθmax : C’est une extension du critère de
Griffith [LEM 76]. La fissure se propage suivant la direction θ pour laquelle, à l’extrémité de l’ex-
tension de longueur donnée s et orientée par θ, le taux de restitution d’énergie G(s) est maximal.
Ce critère donne des résultats similaires à ceux du critère MTS [LEM 04].
32
Mécanique de la rupture en présence de non linéarités de contact avec frottement
l’extension virtuelle de la fissure [AME 79]. Les premiers termes des développements limités de
ces facteurs d’intensité des contraintes sont obtenus grâce aux potentiels complexes de Muskheli-
shvili et des techniques de représentation conformes [RIC 68, BUI 78]. Ils s’expriment en fonction
des facteurs d’intensité des contraintes calculés en pointe de fissure :
k10 (s, θ)
K11 (θ) K12 (θ) KI
= (1.45)
k20 (s, θ) K21 (θ) K22 (θ) KII
On note k1∗ (s, θ) et k2∗ (s, θ) les valeurs de ces FIC pour une longueur de segment infinitésimale :
De nombreux essais ont été nécessaires pour prendre en compte la diversité des mécanismes
de propagation mis en jeu en mode mixte et pour proposer des critères de propagation fiables. Ils
permettent d’obtenir des prévisions en accord avec les essais sous chargement uniaxial et mul-
tiaxial proportionnel dans la plupart des cas. Les critères principaux de direction de propagation
en mode mixte non proportionnel utilisés dans la littérature sont présentés dans la section suivante.
Critères de Hourlier : À partir des travaux d’Amestoy [AME 79, TRU 81, AME 87], Hourlier et
al. [HOU 82] proposent trois critères en espace et en temps basés sur la valeur des FIC à l’extrémité
d’une extension infinitésimale de la fissure (s → 0) :
i. Critère 1 : k1∗ (θ,t)max . La fissure se propage dans la direction θ pour laquelle k1∗ atteint son
maximum absolu en espace et en temps ;
ii. Critère 2 : ∆k1∗ (θ)max . La fissure se propage dans la direction θ pour laquelle l’écart de
variation de k1∗ (θ) est maximal pendant le cycle ;
da
iii. Critère 3 : dN (θ)max . La fissure se propage dans la direction θ pour laquelle le taux d’ac-
da
croissement de la fissure dN (θ) = f (k1∗ (θ,t), k2∗ (θ,t)) est maximal.
33
1. Synthèse bibliographique
F IGURE 1.19: Exemple de variation du critère direction k1∗ (θ)max au cours d’un cycle de
chargement proportionnel et non proportionnel.
Les directions obtenues à partir de ces critères ont été comparées à des résultats expérimentaux
[DHO 82]. Le critère 2 a fourni une estimation en accord avec les essais sous chargement non
proportionnel (mode I cyclique et mode II statique). Le critère 3 a donné des résultats acceptables.
Par contre, le critère 1 n’a pas permis de rendre compte du comportement de la fissure dans ce
cas. D’après [DHO 82], ce critère est sensible au phénomène de refermeture. Il donne cependant
de bons résultats pour un rapport de charge R élevé.
Pour des essais sous chargement de fatigue tribologique, ces critères donnent une bonne ap-
proximation de la direction de branchement [LAM 96, LAM 97a, BAI 03, CHA 06c, RIB 07a].
On pourra donc les utiliser pour la simulation de problèmes de rupture sous chargement multiaxial
non proportionnel tels que les problèmes de fretting fatigue.
Remarque 6 : L’effet de refermeture et le contact avec frottement entre les faces de la fissure
ont une influence importante sur la direction de propagation [DHO 82]. La modélisation de ces
phénomènes est donc indispensable afin de simuler précisément la propagation des fissures.
Des critères de propagation en fatigue multiaxiale sous chargement non proportionnels adaptés
aux problèmes de fretting ont été proposés. Ils permettent de calculer la direction de propagation
à partir des variables de la mécanique de la rupture présentées dans la section 1.2. Les lois de
propagation en vitesse dépendent du mode de sollicitation et du matériau du problème. Elles sont
généralement déterminées empiriquement en couplant l’observation expérimentale de la propaga-
tion des fissures et le calcul des facteurs d’intensité de contraintes. Cette analyse puis la simulation
numérique de la propagation des fissures requiert une méthode numérique permettant de modéliser
le problèmes de rupture bidimensionnels et tridimensionnels, les géométries des fissures de fretting
relativement complexes en 3D et les phénomènes multi-échelles tels que le contact et le frottement
entre le faces des fissures. Dans la section suivante, la méthode des éléments finis est présentée.
Elle s’inscrit comme la méthode la plus adaptée à la simulation de problèmes d’interfaces mobiles
dans le temps.
34
Simulation de la propagation des fissures en fatigue par la méthode des éléments finis étendus
Dans la littérature, on trouve un certain nombre de méthodes numériques utilisées pour la mo-
délisation des problèmes de rupture (étape (i)). La méthode la plus répandue est la méthode des
éléments finis (MEF) [DHO 98, CAR 00]. Sa robustesse, son aptitude à représenter des géomé-
tries complexes en 3D et la possibilité de prendre en compte des lois de comportement matériaux
complexes en font un outil privilégié pour traiter une grande variété de problèmes. Cependant, la
modélisation de problèmes de rupture par la MEF présente deux difficultés majeures. D’une part,
le maillage des éléments finis doit être conforme à la géométrie de la fissure. De plus, l’utilisa-
tion d’éléments finis standards ne permet pas de capturer précisément la singularité du champ de
contrainte en pointe de fissure. On utilise généralement un maillage très fin d’éléments singuliers
dits de Barsoum (taille des éléments en pointe de fissure de l’ordre d’un centième de la taille de la
fissure) [BAR 74], ce qui accroît les temps de calcul. De plus, le maillage de la structure doit être
reconstruit lorsque la fissure évolue. Malgré les améliorations apportées aux mailleurs automa-
tiques, le coût de cette opération peut rapidement devenir rédhibitoire en 3D pour des géométries
de fissures complexes. Il est aussi nécessaire de projeter des champs de l’ancien vers le nouveau
maillage pour les problèmes d’évolution [RéT 05]. Cette opération est non seulement très coûteuse
35
1. Synthèse bibliographique
numériquement mais aussi très délicate car il faut respecter les équations de conservation. Des dé-
veloppements récents tels que les techniques de morphing de maillage couplées à des méthodes
de réduction de modèle permettent d’éviter ces opérations dans certains cas encore relativement
simples [GAL 10].
D’autre méthodes comme la méthode des éléments de frontière (BEM) [POT 92, MI 94,
BON 95, TUH 97, LEI 00, NIK 01, KOL 05] ou les méthodes sans maillage [BEL 94, DUF 04,
KRY 99, FLE 97, RAO 04] ont été appliquées aux problèmes de fissuration.
En particulier, les développements de modèles semi-analytiques à partir de la méthode des dis-
tributions continues de dislocations permettent une résolution très précise des problèmes de rup-
ture avec contact et frottement et mettent en jeu des coûts numériques et des temps de calcul très
faibles [COM 77, DUB 92a, DUB 92b, DUB 92c, HIL 96, KAL 90]. Ils ont notamment été validés
pour des chargements de fatigue de contact (fretting, roulement) [BAI 02b, BAI 02a, LAM 97a].
Cependant, ils ne sont valides que pour des géométries simples et des matériaux homogènes sup-
posés élastiques isotropes et sont difficilement applicables aux problèmes tridimensionnels. Dans
la suite de ce travail, cette méthode sera utilisée comme référence pour la validation du calcul des
facteurs d’intensité des contraintes en 2D dans le chapitre 3 de ce manuscrit.
La méthodes des éléments finis étendus (X-FEM) est une extension de la MEF. Elle est parti-
culièrement adaptée à la modélisation de problèmes de rupture et à la simulation de propagation
de fissures. En effet, elle possède les principaux avantages de la MEF mais il n’est pas nécessaire
de prendre en compte les fissures lors du maillage de la structure. Cette méthode connaît un essor
important dans de nombreux secteurs de recherche industriels. Elle est présentée en détails dans
la section suivante.
36
Simulation de la propagation des fissures en fatigue par la méthode des éléments finis étendus
où r et θ sont les coordonnées polaires dans le repère centré sur le front de fissure (cf. figure
√
1.21). Le terme r permet d’introduire la singularité dans l’approximation du champ de dépla-
cement. Ces fonctions forment une base de la solution asymptotique en pointe de fissure (1.4).
La discontinuité en arrière du front de fissure est assurée par la fonction B2 discontinue lorsque
θ = ±π. De nombreux travaux portant sur ces fonctions singulières ont permis d’améliorer la
précision de la méthode [BOR 07, WYA 08, BOR 06, XIA 06, CHA 06a], notamment en 3D
[SUK 00, MOë 02, GRA 02, ARE 05]
L’approximation du champ de déplacement enrichi s’écrit :
4
u(x,t) ' ∑ ui (t)Φi (x) + H(x) ∑ a j (t)Φ j (x) + ∑ Bl ∑ blk (t)Φk (x) (1.48)
i∈N j∈N f iss l=1 k∈N f ront
où ui sont les degrés de liberté standards éléments finis, Φi sont les fonctions de forme des élé-
ments finis et a j et blk sont les degrés de liberté supplémentaires liés à la fonction saut H et aux
fonctions d’enrichissements singuliers Bl respectivement [MOë 99]. Grâce aux enrichissements
asymptotiques en pointe de fissure, la solution du problème à précision donnée est obtenue pour
des maillages environ 10 fois plus grossiers que ceux utilisées avec la MEF [MOë 02, GRA 02].
Afin de réduire le nombre de degrés de liberté pour des raisons de coûts numériques, tous les
noeuds du maillage d’éléments finis ne sont pas enrichis. Les noeuds des éléments traversés par la
37
1. Synthèse bibliographique
fissure sont enrichis avec la fonction saut et les éléments contenant le front de fissure sont enrichis
avec les fonctions singulières. On parle d’enrichissement topologique local (cf. figure 1.22).
Une autre stratégie, dite d’enrichissement géométrique, propose d’étendre la zone d’en-
richissement singulière à une région de taille fixée indépendamment de la discrétisation. La
zone d’enrichissement discontinue reste identique. Cette technique permet d’améliorer l’ordre
de convergence de la méthode et la qualité de la solution au voisinage du front de fissure
[LAB 05, BéC 05, CHA 06a, NIC 09]. Cependant, le conditionnement de la matrice de raideur est
fortement détérioré. Un pré-conditionneur numériquement coûteux est alors requis. Par ailleurs,
d’autres développements de la méthode ont été proposées pour améliorer la qualité de la solution :
traitement des blending éléments (éléments qui ne vérifient pas la partition de l’unité), nouveaux
enrichissement discontinus [HAN 04, MOë 06].
Dans le cadre de ce travail, une version standard (enrichissement topologique) de la X-FEM
est utilisée. Cette technique permet de limiter au maximum le nombre de degrés de liberté in-
troduits par les enrichissements. En effet, la précision obtenue localement est suffisante pour les
applications traitées dans ce mémoire.
38
Simulation de la propagation des fissures en fatigue par la méthode des éléments finis étendus
information géométrique. À chaque pas de propagation un nouveau segment est ajouté [BEL 99].
Il est possible d’étendre cette technique aux problèmes tri-dimensionnels. On effectue alors un
triangulation de la surface de la fissure. L’étape de propagation peut toutefois s’avérer complexe.
La modélisation précise des fissures en trois dimensions reste délicate, notamment dans le cas
de fissures non planes. Les techniques de représentation implicites telles que la Fast Marching Me-
thod (FMM) [CHO 03, SUK 03] ou la méthode des fonctions de niveau (Level Set Method : LSM)
fournissent une représentation précise des fissures en 3D et sont très bien adaptées à la représenta-
tion de surfaces mobiles au cours du temps. Ces méthodes ont été développées initialement pour
la mécanique des fluides, l’imagerie ou la CAO. C’est la méthode des fonctions de niveau qui
est adoptée dans ce travail. Cette méthode, initialement proposée par Sethian et Osher [OSH 88],
consiste en une représentation implicite d’une interface par un champ scalaire de distances signées
défini sur l’ensemble du domaine (cf. figure 1.23).
Les fonctions de niveau présentent l’intérêt d’être régies par des lois d’évolution de type eu-
lérien qui sont tout à fait adaptées aux problèmes d’interface mobile. On peut décrire mathémati-
quement une fonction de niveau φ(x) associée à une surface Γ par :
où x désigne le point courant de l’espace, x∗ est le point de Γ le plus proche de x et nΓ est le vecteur
normal unitaire à la surface (suivant l’application, son orientation peut être arbitraire). La surface
Γ est localisée par l’iso-0 de φ. La fonction de niveau doit alors respecter les deux conditions
suivantes :
39
1. Synthèse bibliographique
L’expression (1.51) signifie que la fonction de niveau est un champ de distances, on dit alors
qu’elle est initialisée à la distance signée. On obtient la loi d’évolution en temps de la fonction de
niveau en dérivant la relation (1.50) :
∂φ
+Vφ k∇φk = 0 (1.52)
∂t
où Vφ est un champ de vitesse défini sur tout le volume.
Dans la cadre de la mécanique de la rupture, Moës et Gravouil [MOë 02, GRA 02] ont proposé
d’utiliser un jeu de deux fonctions de niveau pour modéliser la fissure. Une première, notée ϕ, est
utilisée pour représenter la surface de discontinuité. La seconde, notée ψ, permet de modéliser le
front de fissure localisé à l’intersection des surfaces iso-0 de ces fonctions de niveaux (cf. figure
1.24).
On peut également exploiter ces deux fonctions de niveau pour définir le repère local associé
au front de fissure. Il est alors nécessaire qu’elles soient localement orthogonales. On obtient
40
État de l’art de la modélisation des problèmes de contact dans le cadre de la méthode des
éléments finis étendus
aisément les coordonnées d’un point x dans la base locale en fonction des valeurs des fonctions
de niveaux ϕ(x) et ψ(x) en ce point. De plus, ce formalisme permet de construire facilement les
fonctions d’enrichissement X-FEM au front de fissure. La définition des fonctions de niveau a
donc une incidence importante sur la qualité du calcul des facteurs d’intensité des contraintes. Les
aspects numériques de ce calcul sont décrits dans la section suivante.
La méthode des éléments finis étendus est une méthode tout à fait adaptée à la simulation
numérique de la propagation des fissures. Elle permet de modéliser des interfaces mobiles dans le
temps et d’économiser les coûts numérique associés aux étapes de remaillages et de projection de
champs normalement requises avec la MEF. De plus, la technique des fonctions de niveau permet
de modéliser des géométries de fissures complexes en 3D et facilite les étapes de propagation géo-
métrique et numérique (mise à jour des enrichissements X-FEM) des fissures [GRA 02, DUF 07].
Depuis une dizaine d’années, de nombreux travaux autour de la méthode des éléments finis éten-
dus ont été publiés, notamment en mécanique de la rupture. Des applications ayant trait à la rup-
ture par fatigue [DUF 04, ELG 07b, RIB 07b, ARE 05, RAN 08b, WYA 08, GIN 09a, RAN 10],
à la rupture fragile élasto-dynamique [RéT 05, MEN 06, MEN 08, COM 08, GRé 07, PRA 08,
ELG 09, GRA 09, ROZ 08] ou encore à la déchirure ductile [BOR 06] ont clairement mon-
tré l’intérêt de la méthode X-FEM. Des développements récents ont été proposés afin d’amé-
liorer la qualité de la solution dans la zone enrichie et d’atteindre des précisions élevées
[DRé 10, LEG 08, LAB 05, BéC 05, CHA 06a, NIC 09]. La X-FEM a atteint un niveau de ma-
turité qui lui vaut désormais d’être implantée dans un nombre croissant de codes commerciaux et
industriels. L’atout principal de cette technique est qu’elle peut s’intégrer relativement facilement
dans des codes éléments finis existants.
Cependant, la méthode standard ne prévoit pas la modélisation du contact avec frottement
entre les faces de la fissure. De nombreux travaux sur les problèmes de contact dans le cadre de
la méthode des éléments finis étendus ont mis en évidence un certain nombre de difficultés. L’ob-
jectif de la section suivante est d’identifier les problèmes associées à la modélisation du contact
avec frottement à l’interface avec X-FEM et d’évaluer la robustesse des différentes méthodes em-
ployées afin de proposer un modèle optimal pour la simulation des problèmes de fissuration sous
chargement de fretting.
Le contact avec ou sans frottement interfacial a été implémenté dans de nombreux modèles
41
1. Synthèse bibliographique
X-FEM en 2D [DOL 01, VIT 08, GIN 08, BéC 09, GUI 08] et en 3D [GéN 06, GéN 07, SIA 10,
PIE 10a, LIU 10]. Plusieurs modèles 2D quasi-statiques ont notamment été proposés pour simu-
ler la propagation de fissures en fatigue avec contact et frottement interfacial sous chargement
multiaxial non proportionnel [RIB 07b, LIU 08, GIN 09b, GIN 09a, BAI 10] et avec contact uni-
latéral et plasticité confinée en pointe de fissure [ELG 07b]. La modélisation du contact avec/sans
frottement interfacial pose toutefois deux difficultés majeures.
Premièrement, les problèmes de rupture avec contact et frottement font intervenir plusieurs
phénomènes à différentes échelles : celle de la structure, celle de la fissure et celle des non-
linéarités à l’interface. Ces échelles peuvent différer les unes des autres de plusieurs ordres de
grandeur. Dans le cadre de la X-FEM, des modèles multi-échelles basés sur des techniques multi-
grilles [PAS 10, RAN 08b, FIS 05] ou de décomposition de domaine [GUI 08] ont été proposés
afin de pallier cette difficulté.
Une solution de contact précise requiert une discrétisation spatiale de l’interface adaptée à
la dimension caractéristique du problème local [RIB 07b]. La plupart des modèles présentés ci-
dessus ne permettent pas une description multi-échelle du problème de contact avec frottement. En
effet, la discrétisation du problème d’interface est lié au maillage de la structure. On procède alors
à un raffinement local du maillage de la structure dans la zone d’intérêt afin d’adapter la discréti-
sation de l’interface (cf. figure 1.25). Dans le cas des problèmes de fretting fatigue par exemple, le
rapport des échelles mises en jeux est très élevé. On observe expérimentalement des variations de
la taille de la zone de contact adhérent ou glissant entre les faces des fissures de l’ordre du micron
alors que les dimensions de la structure peuvent être de l’ordre du mètre. L’opération de maillage
du domaine peut alors s’avérer délicate, voire impossible pour les problèmes tridimensionnels. De
plus, au regard du problème de la structure, le maillage n’est pas optimal et entraine un surcoût
numérique. Par ailleurs, dans le cadre de ce travail, ce procédé est en contradiction avec le concept
de la méthode des éléments finis étendus. En effet, on cherche à modéliser une interface (telle
qu’une fissure) indépendamment du maillage en volume.
42
État de l’art de la modélisation des problèmes de contact dans le cadre de la méthode des
éléments finis étendus
Deuxièmement, plusieurs travaux ont mis en évidence des problèmes d’instabilité venant per-
turber la solution de contact (avec ou sans frottement). Les méthodes concernées font apparaitre
des oscillations numériques (cf. figure 1.26) dans le champ discrétisé des efforts d’interface. Ce
problème est lié à un non respect de la condition inf-sup LBB (Ladyzhenskaya-Babuška-Brezzi).
Il dépend fortement de la discrétisation du problème et de la formulation utilisée. On distingue
alors trois types de formulation du problème de rupture avec contact et/sans frottement entre les
faces de la fissure induisant des comportement différents.
i. Les formulations primales [LIU 08] pour lesquelles le problème de contact est formulé en
déplacement (cf. figure 1.27.a). Une méthode d’enrichissement X-FEM spécifique permet
de modéliser la discontinuité géométrique liée à la présence de la fissure et d’intégrer la loi
de comportement d’interface (contact avec frottement). Le champ des efforts d’interface t
est considéré comme une composante des efforts externes. Au cours de la résolution itéra-
tive du problème global, il est re-calculé à partir de la solution en déplacement u. Le modèle
proposé par Liu et al. [LIU 08] est stable. Il ne présente pas de problèmes d’oscillations. Par
ailleurs, il est montré que la convergence de la méthode de résolution non linéaire utilisée
est très élevée en comparaison avec les modèles proposés dans la littérature.
Note : Dans ce modèle, la discrétisation du problème de contact avec frottement d’interface
est fortement dépendante du maillage de la structure. Une solution de contact fine requiert
un raffinement local du maillage en volume de la structure.
ii. Les formulation duales (cf. figure 1.27.b) sont généralement basées sur la méthode des
multiplicateurs de Lagrange. Le problème de contact est formulé entre le champ discrétisé
des déplacements u et le champ discrétisé des multiplicateurs de Lagrange λ équivalent aux
efforts de contact. La majorité de ces modèles sont concernés par les problèmes d’instabilité
numériques : oscillations, pics de contraintes et convergence faible). L’expérience montre
que le choix des espaces de multiplicateurs de Lagrange pour les problèmes normal et
tangentiel à l’interface influence fortement la stabilité de la solution. Cette question est
largement détaillée dans les références [HAN 04, BEC 09, NIS 08, LAD 02]. Moës et al.
ont proposé une méthode de réduction de cet espace [MOë 06, BéC 09, GéN 07, SIA 10]
dans le but de stabiliser la solution de contact. Pour ce modèle, la discrétisation du problème
d’interface est également fortement dépendante du maillage de la structure. Une description
43
1. Synthèse bibliographique
44
État de l’art de la modélisation des problèmes de contact dans le cadre de la méthode des
éléments finis étendus
45
1. Synthèse bibliographique
46
Chapitre 2
Sommaire
2.1 Formulation faible mixte du problème de structure avec interface in-
terne en présence de contact et frottement . . . . . . . . . . . . . . . . . 49
2.2 Discrétisation X-FEM multi-échelle dédiée aux problèmes de contact
et frottement interfacial . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
2.3 Solution non linéaire du problème de rupture avec contact et frotte-
ment interfacial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
2.3.1 La méthode du Lagrangien Augmenté . . . . . . . . . . . . . . . . 56
2.3.2 Résolution non-linéaire basée sur la méthode LATIN . . . . . . . . 60
2.3.3 Stabilisation de la méthode LATIN . . . . . . . . . . . . . . . . . . 65
47
2.4 Exemples tridimensionnels et validation du modèle . . . . . . . . . . . . 71
2.4.1 Problème tridimensionnel de contact avec frottement : stabilité de
la méthode numérique . . . . . . . . . . . . . . . . . . . . . . . . 71
2.4.2 Problème de contact unilatéral : intérêt de l’approche multi-échelle 80
2.4.3 Problème de contact avec frottement : comparaison avec une solu-
tion numérique de référence . . . . . . . . . . . . . . . . . . . . . 85
2.5 Modélisation quasi-statique et simulation multi-échelle en temps de la
propagation des fissures . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
48
Formulation faible mixte du problème de structure avec interface interne en présence de contact
et frottement
Comme nous l’avons vu, la modélisation dans le cadre de la X-FEM du contact et du frotte-
ment entre les faces d’une fissure ou entre deux corps demeure délicate. Un modèle basé sur une
formulation mixte issue des méthodes de décomposition de domaine est proposé dans le but de
s’affranchir des difficultés posées. Cette formulation présente a priori de bonnes propriétés de sta-
bilité [ELG 07a, RIB 07b]. Elle est présentée dans la première section de ce chapitre. L’interface
est considérée comme une structure autonome avec ses propres variables primales et duales et sa
propre discrétisation [LAD 85]. Cette considération permet une description multi-échelle précise
et optimale du problème local de contact à l’échelle de la fissure indépendamment du problème
global à l’échelle de la structure [PIE 10a].
On considère l’interface Γ (d’une fissure, entre deux corps ou deux milieux) comme une struc-
ture autonome possédant ses propres variables primales et duales et sa propre discrétisation. Au
temps t, le problème est alors divisé en deux sous problèmes [LAD 99] :
- le problème local de contact à l’échelle de la fissure. Les inconnues de ce problème sont
le champ des déplacements surfaciques w(t) et le champ des efforts surfaciques de contact
t(t). La solution de ce problème obéit aux lois de comportement d’interface normale et
tangentielle (loi de contact unilatéral et loi de coulomb par exemple) :
49
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
Afin de construire une formulation faible de ce problème, l’idée consiste à utiliser la méthode
des multiplicateurs de Lagrange à partir du principe des puissances virtuelles (PPV) [GER 80].
On définit alors les puissances virtuelles associées au problème global en volume et au problème
local à l’interface :
Z Z
∗ ∗ ∗
Pglobal = Tr [σ(t)ε(u )] − f(t) · u dV − F(t) · u∗ dS (2.8)
Ω ∂2 Ω
où les espaces fonctionnels des champs de déplacement réel u et virtuel u∗ sont définis par :
u ∈ Ud = u/u|∂1 Ω = ud + régularité
(2.9)
u∗ ∈ U0∗ = u∗ /u∗ |∂1 Ω = 0 + régularité
et Z
∗
Plocal =− t(t) · w∗ dS (2.10)
Γ
où w∗ est le champ de déplacement virtuel à l’interface tel que :
50
Formulation faible mixte du problème de structure avec interface interne en présence de contact
et frottement
On peut alors construire une formulation faible associée à l’ensemble du problème à partir du
PPV modifié :
∗ ∗ ∗
Pglobal + Plocal + Pcouplage = 0 ∀u∗ ∈ U0∗ , ∀w∗ ∈ W ∗ , ∀λ∗ ∈ Λ∗ , ∀t ∈ [0; T ] (2.14)
équivalent à :
∗ )] ∗ ∗ ∗
R R R R
0 = − Ω Tr [σ(t)ε(u dV + ∂2 Ω F(t) · u dS + Ω f(t) · u dV + Γ λ(t) · u dS
∗
R
+ Γ (t(t) − λ(t)) · w dS
(2.15)
∗
R
+ Γ (u(t) − w(t)) · λ dS
51
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
Remarque 7 : Dans ce travail, on fait l’hypothèse des petits déplacements et des petites déforma-
tions. En conséquence, le calcul de l’état local de contact est toujours réalisé dans la configuration
non déformée.
La formulation faible mixte à trois champs (2.15) est la base du modèle X-FEM à deux
échelles développé dans ce mémoire. On considère les problèmes local et global comme deux
sous-problèmes "autonomes". Ils sont reliés au sens faible par une condition de continuité des
déplacements et des efforts. Cette approche permet l’utilisation de discrétisations indépendantes
pour chaque sous-problème, et donc une modélisation multi-échelle du problème général. Dans
la section suivante, une méthode est proposée afin d’appliquer cette stratégie dans cadre de la
méthode des éléments finis étendus.
52
Discrétisation X-FEM multi-échelle dédiée aux problèmes de contact et frottement interfacial
en plusieurs triangles lorsqu’au moins une de ses arrêtes est plus longue qu’une dimension ca-
ractéristique lc définie par l’utilisateur. Ce critère permet d’uniformiser la distribution des points
d’intégration sur l’interface et de l’adapter à l’échelle caractéristique du problème de contact.
De plus, un critère de forme est utilisé afin d’optimiser le sous-découpage des triangles. Pour
cela, les rapports des longueurs entre des différentes arrêtes des éléments d’interface sont calcu-
lés. Ce critère permet de détecter les éléments allongés. En fonction de leur degré de distorsion,
on sélectionne le mode de subdivision optimal afin d’obtenir des triangles plus réguliers. Ce déve-
loppement a pour but d’améliorer la qualité de l’intégration numérique du problème local. Chaque
élément d’interface est testé puis sous-découpé si nécessaire jusqu’à l’obtention d’une distribution
de points d’intégration uniforme et adaptée à l’échelle requise. Ce processus a été automatisé et
est très peu coûteux numériquement. Il est illustré dans la figure 2.3.
Remarque 8 : Il ny a pas de dégrés de liberté associés aux noeuds des éléments d’interface. Le
"maillage" des éléments d’interface n’est donc pas nécessairement conforme. La loi de contact
locale s’applique pour chaque paire de points d’intégration en vis à vis indépendamment des
autres paires.
Un premier exemple tridimensionnel simple est proposé pour illustrer l’efficacité cette stra-
tégie. On considère un domaine Ω parallélépipédique de dimensions (120mm×120mm×120mm)
et une fissure débouchante Γ non plane à front courbe (cf. figure 2.4). On considère le maillage
élément finis structuré relativement grossier représenté dans la figure 2.5. La figure 2.6 représente
les discrétisations de l’interface de la fissure par la méthode (a) (cf. figure 2.2) et avec la méthode
de raffinement (b) proposée pour des valeurs de lc = 5mm et lc = 3mm.
La méthode de raffinement permet d’adapter la discrétisation de l’interface à l’échelle du pro-
blème de contact indépendamment du maillage de la structure. Cette modélisation à deux échelles
est possible grâce au couplage au sens faible des problèmes local et global décrit par l’équation
53
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
F IGURE 2.4: Géométrie du domaine (en mm) et de la fissure définie par deux fonctions
de niveau φ et ψ.
(2.15). Les gains numériques liés à cette stratégies est illustré dans la section 2.4. Dans la section
suivante, nous proposons deux méthodes de résolution non linéaires pour ce formalisme.
54
Solution non linéaire du problème de rupture avec contact et frottement interfacial
55
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
(2.20) et (2.4).
∗
dS + Ω fn+1 · u∗ dV
R R
+ ∂2 Ω Fn+1 · u
R (i) (i)
R (i+1) (i+1)
+ Γ tn+1 + αwn+1 w∗ dS − Γ λn+1 + αwn+1 w∗ dS (2.21)
R (i+1) (i+1)
+ Γ un+1 − wn+1 λ∗ dS
Cette formulation est valable quelle que soit la méthode de discrétisation en espace choisie
(X-FEM, MEF, BEM,...). Dans le cadre du modèle X-FEM à deux échelles présenté dans la sec-
tion précédente, on exprime le champ de déplacement global à l’aide de l’équation (1.48). Les
champs locaux s’expriment en fonction des fonctions de formes primales et duales associées à la
discrétisation de l’interface :
3
w(x,t) ' ∑ wm (t)Ψm (x) (2.23)
m=1
3
t(x,t) ' ∑ tm (t)Ψ0m (x) (2.24)
m=1
3
λ(x,t) ' ∑ λm (t)Ψ0m (x) (2.25)
m=1
56
Solution non linéaire du problème de rupture avec contact et frottement interfacial
où Ψ et Ψ0 sont respectivement les bases des fonctions de forme primales et duales linéaires
de éléments d’interface d’ordre 0 ou 1 avec trois degrés de liberté par point d’intégration pour
les déplacements locaux (wx , wy , wz ) et les efforts locaux (tx , ty , tz ). On fait une approximation
similaire pour les champs virtuels w∗ , t∗ et λ∗ .
On définit alors les champs discrétisés du problème par :
ui Fui fui
U ≡ aj F≡ 0 f≡ 0 (2.26)
blk n 0 n
0 n
W ≡ wi T ≡ ti Λ ≡ λi (2.27)
p p p
où
Ui+1 = ∆Ui+1 + Ui
Wi+1 = ∆Wi+1 + Wi (2.29)
Λi+1 = ∆Λi+1 + Λi
K est la matrice de raideur globale du problème de structure enrichi. Kww est l’opérateur
associé au coefficient de pénalité α :
Z
Kww = α Ψi Ψ j dS (2.30)
Γ
57
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
Kuλ et Kwλ sont équivalents à des opérateurs mortier (ou mortar) [MCD 00, COM 03,
DUR 06, KIM 07, RAN 08a]. Ils autorisent le couplage au sens faible entre les problèmes local et
global et l’utilisation de discrétisations non compatibles :
R R R
Kuλ = Γ Φi Ψm dS Γ HΦ j Ψm dS Γ Bl Φk Ψm dS (2.31)
Z
Kwλ = Ψi Ψ0j dS (2.32)
Γ
À chaque itération le système linéaire (2.28) est résolu, puis on calcule le champ des efforts
d’interface Ti+1 à partir de Wi+1 et des équations de la loi de comportement d’interface (2.20).
On considère qu’on obtient une solution convergée au sens d’un indicateur d’erreur pertinent η tel
que :
(i)
lim Xn+1 ≡ Xn+1 (2.33)
η→0
Un critère de convergence local spécifique aux problèmes de contact a été développé par Ri-
beaucourt et al. [RIB 07a]. Il est défini comme la distance entre deux approximations locale et
globale successives, en terme de déplacements et d’efforts à l’interface en distinguant les gran-
deurs relatives au problème normal et au problème tangentiel :
η = max(ηN ; ηT ) (2.34)
avec
(i+1) (i+ 21 ) (i+1) (i+ 12 )
k XN − XN k2∞ k XT − XT k2∞
ηN = , ηT = (2.35)
(i+1) (i+ 12 ) (i+1) (i+ 12 )
k XN k2∞ + k XN k2∞ k XT k2∞ + k XT k2∞
où la norme est définie par :
1
k X k2∞ = max(αT2 + W2 ) (2.36)
α
On peut également utiliser le critère d’erreur global basé sur la norme du second membre du
système linéaire (2.28) [ELG 07a]. Ce dernier permet de contrôler le respect de la condition de
couplage entres les problèmes local et global.
Le premier problème est proposé. Il a pour objectif de démontrer les propriétés de conver-
gence de cette méthode, appelée AL (Augmented Lagrangian) dans la suite de ce manuscrit. On
considère un domaine 2D carré Ω de dimensions (0.1m×0.1m) soumis à une contrainte de com-
pression verticale P=1 MPa. Une interface Γ, perpendiculaire à la direction du chargement, divise
Ω en deux parties égales. On impose une condition de déplacement nul sur la surface inférieure
(uX = uY = 0) et une condition de déplacement horizontal nul sur la surface supérieure (uX =0). Le
matériau est homogène, élastique et isotrope de module de Young E=206 GPa et de coefficient de
Poisson ν=0.3. Le coefficient de frottement de la loi de Coulomb à l’interface est µΓ = 0.2.
Deux maillages structurés sont considérés. Le premier est constitué de 722 éléments finis
linéaires triangulaires et le second de 361 éléments finis linéaires quadrangulaires. La solution de
contact est relativement uniforme sur toute l’interface. Il n’est donc pas nécessaire de procéder à
un raffinement spécifique de la discrétisation du problème de contact. Les bases des fonctions de
forme primale et duale des éléments d’interface sont d’ordre 0 (un point d’intégration pour chaque
face par élément d’interface).
58
Solution non linéaire du problème de rupture avec contact et frottement interfacial
Le problème est résolu à une précision de 10−4 en considérant le critère d’erreur local η
(2.34). Différentes valeurs du coefficient de pénalité α ont été testées afin d’obtenir un taux de
convergence optimal dans chaque configuration.
Dans les deux cas, un taux de convergence optimal est obtenu pour une valeur du coefficient
de pénalité proche de α = 2.10+11 Pa.m−1 . Pour cette valeur, le problème est résolu en 35 et 36
itérations avec le maillage de triangles et le maillage de quadrangles respectivement. L’expérience
montre que cette valeur optimale est proche de la valeur du module de Young pour des problèmes
élastiques linéaires. La méthode AL permet d’obtenir une solution globale très satisfaisante d’un
point de vue qualitatif (cf. figures 2.8.a et 2.9.a). De plus, la solution de contact est parfaitement
stable et régulière (cf. figures 2.8.b et 2.9.b) et ne présente pas d’oscillations numériques.
59
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
Remarque 12 : On obtient des résultats quasi-identiques avec des bases des fonctions de forme
primale et duale d’ordre 1 pour les éléments d’interface (deux points d’intégration pour chaque
face par élément d’interface linéaire).
60
Solution non linéaire du problème de rupture avec contact et frottement interfacial
axial non proportionnel) [RIB 07b]. Par ailleurs, la méthode LATIN a récemment été adaptée
par Guidault et al. dans le cadre d’un modèle multi-échelle X-FEM couplé à une technique de
décomposition de domaines pour la simulation de la propagation de fissures [GUI 08]. Ce modèle
a été étendu aux problèmes de contact et frottement interfacial dans [GUI 10].
Le concept de la méthode LATIN repose sur la division des équations du problème en deux
sous-ensembles :
- un sous-ensemble L d’équations locales éventuellement non linéaires. Ici, cet ensemble
correspond aux équations de la loi de comportement d’interface (2.20) ;
- un sous-ensemble G d’équations linéaires, généralement globales. Dans notre approche, cet
ensemble est composé de l’équation établissant la formulation faible mixte (2.15), de la loi
de comportement en volume (2.17) et des conditions aux limites (2.4).
(0)
À partir d’une solution initiale Xn+1 issue d’un calcul élastique linéaire, une approximation
de la solution du problème non linéaire Xn+1 = (un+1 , σn+1 , wn+1 , tn+1 ) ≡ G ∩ L à un instant tn+1
donné est obtenue par un processus itératif en deux étapes illustré schématiquement dans la figure
2.11. Des conditions supplémentaires sont définies entre les solutions des sous-espaces L et G,
appelées directions de recherches locales E+ et globales E− . Elles s’expriment par :
(i+ 21 ) 1
1
(i)
Xn+1 − Xn+1 ∈ E+ =⇒ t(i+ 2 ) − t(i) = kl w(i+ 2 ) − w(i) (2.37)
(i+ 12 ) 1
1
(i+1)
Xn+1 − Xn+1 ∈ E− =⇒ t(i+1) − t(i+ 2 ) = −kg w(i+1) − w(i+ 2 ) (2.38)
où kl et kg sont les opérateurs strictement positifs définis par l’utilisateur et exprimés en Pa.m−1
tels que :
kl,N 0 kg,N 0
kl = et kg = (2.39)
0 kl,T 0 kg,T
61
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
(i+ 1 )
i. l’étape locale. On calcule la solution Xn+12 du problème local définie par l’ensemble
d’équations L (i.e. le problème de contact avec frottement entre les faces de la fissure)
et l’équation 2.37. À partir de la solution précédente (w(i) ,t(i) ), on utilise un indicateur de
contact local proposé par Champaney [CHA 96, CHA 99] afin de déterminer les zones de
contact et d’ouverture et l’état de contact tangentiel local : glissement ou adhérence. On
applique alors la loi de comportement d’interface (2.20) localement et l’équation (2.37).
1 1
On obtient alors une nouvelle solution du problème local de contact (w(i+ 2 ) ,t(i+ 2 ) ).
(i+1)
ii. l’étape globale. On cherche une solution Xn+1 du problème global défini par l’ensemble
d’équations G. On introduit alors l’équation (2.38) dans la formulation faible (2.15) et on
obtient :
h i
(i+1) R (i+1)
Kε(un+1 )ε(u∗ ) dV + Γ λn+1 · u∗ dS
R
0 = − Ω Tr
∗ fn+1 · u∗
R R
+ Ω Fn+1 · u
∂2 dS + Ω dV
(i+ 21 ) (i+ 12 ) R (i+1)
(i+1)
· w∗ dS − · w∗ dS
R
+ Γ tn+1 + kg wn+1 Γ λn+1 + kg wn+1
(2.40)
R (i+1) (i+1)
+ Γ un+1 − wn+1 · λ∗ dS
La convergence de la méthode LATIN a été démontrée dans [LAD 99] sous certaines condi-
tions. Si la loi de comportement est monotone et si les interfaces représentent des conditions aux
limites, des liaisons parfaites ou des surfaces de contact sans frottement, il est montré que le pro-
cessus converge lorsque les directions de recherche locale et globale sont conjuguées et positives,
c’est-à-dire kl = kg . Afin d’assurer la convergence de la méthode pour une classe plus large de
62
Solution non linéaire du problème de rupture avec contact et frottement interfacial
comportements, i.e. contact avec frottement, il est classiquement introduit une étape de relaxa-
tion :
Remarque 13 : On peut faire l’analogie entre les formulations (2.21) et (2.40) construites à
partir de la même formulation faible mixte (2.15). La valeur de la direction de recherche k = kl =
kg influence uniquement le taux de convergence de la méthode. Elle ne modifie pas la solution du
problème. En pratique, il existe une valeur quasi optimale kopt liée à la rigidité du matériau et à
la taille de l’interface (longueur en 2D et aire en 3D) [PIE 10a].
À partir des champs discrétisés définis par (2.26) et (2.27), on construit le système linéaire
suivant :
K 0 −Kuλ Ui+1 F
0 Kww Kwλ Wi+1 = Kwλ · Ti+ 12 + Kww · Wi+ 12 (2.42)
−Kuλ KTwλ
T
0 Λi+1 0
où les opérateurs Kuλ , Kww et Kwλ sont identiques à ceux de la méthode AL.
(i+1)
À chaque itération, on calcule le nouveau champ des efforts de contact Tn+1 à partir de la
solution de (2.42) et de l’équation (2.38). Ce processus en deux étapes locale et globale est répété
jusqu’à convergence de la solution. On considère le critère d’erreur locale η dédié aux problèmes
de contact (2.34). On remarque alors que Λn+1 ∼ = Tn+1 à convergence pour une précision donnée.
Remarque 14 : La méthode LATIN, appelée SLA pour "Standard LAtin" dans la suite de ce ma-
nuscrit, est particulièrement bien adaptée aux problèmes multi-échelles. En effet, elle permet de
découpler la résolution du problème local et celle du problème global.
63
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
d’interface et de la forme du système linéaire du solveur non linéaire. En effet, pour le problème
simple traité dans ce début de chapitre, la méthode du Lagrangien Augmenté ne présente pas
d’oscillations numériques. Cependant, pour des problème de contact avec un glissement important,
les expérimentations numériques montrent des oscillations apparaissent quelle que soit le solveur
utilisé : la méthode LATIN ou la méthode du Lagrangien Augmenté. Dans la section suivante, le
problème de la condition LBB est étudié afin de proposer une méthode générale de stabilisation
dédié au modèle X-FEM avec contact et frottement proposé dans ce chapitre.
64
Solution non linéaire du problème de rupture avec contact et frottement interfacial
A BT
Y F
= (2.43)
B 0 Z G
où C est une constante dépendant du problème. Cette propriété garantit la solvabilité, c’est-à-dire
que Y et Z restent bornés quels que soient F et G. En effet, il s’agit d’un autre moyen d’écrire
l’unicité de la solution [AUR 04].
Quelques détails sont donnés ici sur la condition de stabilité à respecter pour la formulation
faible à trois champs présentée dans le paragraphe 2.1 afin de proposer une stratégie de stabilisation
générale des systèmes linéaires associés. On se place dans le cadre de la méthode LATIN. Le
65
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
système linéaire discrétisé (2.42) peut être écrit sous la forme matricielle par blocs suivante :
A BT
Y F
= (2.45)
B 0 Z 0
avec
K 0 T −Kuλ
A≡ B ≡
0 Kww Kwλ
! (2.46)
U F
Y≡ Z ≡ (Λ) F≡
W Kwλ · Ti+ 1 + Kww · Wi+ 1
2 2
En procédant à une élimination de Gauss par blocs, on obtient le complément de Schur associé
àZ:
A BT
Y F
= (2.47)
0 CS Z FS
avec
CS = B A−1 BT et FS = B A−1 F (2.48)
On cherche à savoir si CS est inversible. En mécanique des solides, A−1 est habituellement
inversible si l’on prend en compte les conditions aux limites du problème, et donc, le problème
revient à savoir si B BT est inversible. Pour que B BT soit inversible, il est nécessaire que :
La traduction de cette inégalité algébrique dans le contexte de la méthodes des éléments finis est :
YTh BTh Zh
max >0 (2.51)
Yh ∈Yh kYh kY · kZh kZ
h h
où h est une longueur caractéristique du maillage éléments finis, i.e. celle du plus petit élément par
exemple.
A ce niveau, il est très important de remarquer que des techniques de discrétisation qui dé-
pendent de la résolution spatiale sont considérées. Autrement dit, on considère des familles de
matrices et de solutions qui dépendent de h. Il vient alors :
Ah BTh
Yh Fh
= (2.52)
Bh 0 Zh 0
Par conséquent, il faut étudier la condition (2.51) lorsque h tend vers 0 et vérifier si ces familles dé-
génèrent vers un système singulier. La condition stricte de Ladyzhenskaya-Babuška-Brezzi (LBB)
assure que B BT ne dégénère pas vers 0 lorsque h décroit :
YTh BTh Zh
inf sup ≥β>0 (2.53)
Zh ∈Zh Yh ∈Yh kYh kY · kZh kZ
h h
66
Solution non linéaire du problème de rupture avec contact et frottement interfacial
Propriété :
La condition LBB est vérifiée si et seulement si :
où k et l sont les degrés des polynômes des espaces éléments finis utilisés pour Y et Z respec-
tivement. CY et CZ sont des constantes qui dépendent du problèmes mais pas de h. On retrouve
l’analogue de la condition générale sur la solvabilité des méthodes numériques (2.44).
Pour les formulations mixtes à trois champs, il existe dans la littérature plusieurs méthodes de
stabilisation qui permettent de vérifier la condition (2.54) [AUR 04]. Ici, la méthode de stabilisa-
tion par pénalité est adoptée. On introduit un terme supplémentaire dans la matrice et le second
membre du système linéaire (2.45) :
A BT
Y F
= (2.55)
B −εD Z −εd
où D est un opérateur symétrique défini positif, d est homogène à DZ et ε est un paramètre positif
de pénalité. Le complément de Schur en Z est calculé (pivot de Gauss par blocs) :
A BT
Y F
= (2.56)
0 CS0 Z FS0
avec
CS0 = B A−1 BT + εD et FS0 = B A−1 F + εd (2.57)
D est non singulier et B A−1 BT est positif. B A−1 BT + εD est donc non singulier et la
matrice condensée sur Z est inversible. Vérifions que la matrice condensée sur Y est également
inversible. Il vient :
1 T
A + B D B Y = F + D−1 d
(2.58)
ε
→ A + 1ε BT D B est non singulier car A est non singulier et BT D B est positif.
A BT
Y F
= (2.59)
B −εD Z −εd
α kYk2Y ≤ YT A Y ∀Y ∈ Y (2.60)
67
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
- Condition Inf-Sup : il existe une constante positive β indépendante du maillage (h) telle
que :
YT BT Z
inf sup ≥β>0 (2.61)
Z∈Z \0 Y∈Y \0 kYkY · kZkZ
Bilan :
L’étude de la stabilité du système linéaire (2.45) (ε = 0) se résume à étudier les conditions de
stabilité (2.63) et (2.64) lorsque β tend vers zéro, c’est-à-dire lorsque la condition LBB (2.61)
n’est plus vérifiée. Les second membres des inégalités (2.63) et (2.64) tendent alors vers +∞, ce
qui ne garantit plus la stabilité. En revanche, pour les systèmes linéaires (2.55) avec ε > 0, les
propriétés (2.63) et (2.64) restent bornées, même lorsque β tend vers zéro.
Ici, on généralise cette méthode à une formulation faible à trois champs et au système linéaire
associé à la méthode LATIN. On introduit un opérateur de stabilisation sur la condition de cou-
plage faible entre les champs de déplacements global u et local w à la fois au premier et au second
membre du système linéarisé [PIE 10b]. L’équation (2.42) devient alors :
K 0 −Kuλ Ui+1 F
0 Kww Kwλ Wi+1 = Kwλ · Ti+ 12 + Kww · Wi+ 12 (2.65)
T T
−Kuλ Kwλ Kλλ Λi+1 Kλλ · Λi
68
Solution non linéaire du problème de rupture avec contact et frottement interfacial
69
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
représentent la variation du nombre d’itérations nécessaires pour atteindre cette précision pour
différentes valeurs de la direction de recherche k et du coefficient de stabilisation ε :
- Dans le cas du maillage de triangles, une direction de recherche k = 2.10+13 Pa.m−1 et un
coefficient de stabilisation ε = 1.10+14 Pa.m−1 donnent un taux de convergence optimal (61
itérations).
- Dans le cas du maillage de quadrangles, une direction de recherche k = 1.10+13 Pa.m−1 et
un coefficient de stabilisation ε = 2.10+14 Pa.m−1 donnent un taux de convergence optimal
(28 itérations).
Les solutions locales et globales sont stables dans les deux cas considérés. Elles sont représen-
tées dans les figures 2.17 et 2.18. Ce résultat montre les bonnes propriétés de stabilité de ce solveur
70
Exemples tridimensionnels et validation du modèle
non linéaire appelé NLLA (Non Locking LAtin) pour les problèmes de contact adhérent en 2D.
L’introduction d’un opérateur de stabilisation supplémentaire permet de supprimer les éventuelles
oscillations numériques dues à un non respect de la condition LBB et d’optimiser la convergence
de la méthode. La figure 2.19 représente l’évolution de l’erreur η pour chaque configuration dans le
cas le plus favorable pour les trois méthodes de résolutions utilisées dans ce mémoire. La méthode
NLLA permet d’améliorer les temps de calcul de 94% et de 70% par rapport à la méthode SLA
pour les maillages de triangles et de quadrangles respectivement. Les performances des solveurs
LA et NLLA sont semblables.
Bilan :
Une formulation faible mixte du problème de rupture avec contact et frottement interfacial a été
proposée. Elle permet une description intrinsèque de la fissure dans le cadre de la méthode des
éléments finis étendus. Un solveur basé sur la méthode du Lagrangien Augmenté et un solveur
basé sur la méthode LATIN adaptés de cette formulation sont utilisées. Des problèmes d’instabilité
numériques de la solution locale de contact ont été mis en évidence pour différentes discrétisations
de l’interface et du maillage en volume. Une méthode de stabilisation est alors proposée. Elle
permet de s’affranchir des problèmes d’oscillations numériques et d’augmenter fortement le taux
de convergence de la résolution non linéaire. On peut donc conclure que les propriétés de stabilité
du modèle dépendent principalement du solveur utilisé et de la discrétisation des champs de la
formulation faible. Les expérimentations numériques montrent que les instabilités numériques
sont plus sévères dans le cas des problèmes de contact glissant. De plus, la modélisation multi-
échelle des problèmes tridimensionnels conduit à des discrétisations en volume et à l’interface
généralement plus complexes et plus favorables aux instabilité numériques. Ces configurations
sont susceptibles de conduire à des problèmes de convergence. Afin de tester le modèle 3D X-
FEM à deux échelles proposé et de démontrer les performances de la formulation faible mixte et
du solveur non linéaire stabilisé NLLA, trois problèmes tridimensionnels sont proposés dans la
section suivante.
71
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
bases des fonctions de forme primale et duale des éléments d’interface sont d’ordre 0 avec un
point d’intégration par élément d’interface. Ce qui conduit à une distribution de 1670 paires de
points d’intégrations réparties uniformément sur les deux faces de la fissure (cf. figure 2.21.b). À
l’interface, on considère la loi de contact unilatéral et la loi de frottement de coulomb. Le pro-
blème est résolu pour trois différentes valeurs du coefficient de frottement local µΓ = 1 ; 0.5 ; ou 0,
72
Exemples tridimensionnels et validation du modèle
F IGURE 2.21: Maillage en volume et discrétisation raffinée de l’interface (lc = 1.8 mm).
conduisant à différents états de contact : adhérence, glissement partiel ou glissement total. L’ap-
proximation de la solution est calculée à une précision de 10−3 avec l’indicateur d’erreur (2.34).
On utilise tout d’abord le solveur SLA. L’objectif est de montrer les éventuels problèmes d’in-
stabilité numériques associés à chaque cas. Un taux de convergence maximum est obtenu pour
une direction de recherche k = 1.10+13 Pa.m−1 . Puis, on utilise le solveur NLLA afin d’illustrer
les propriétés de stabilité et les performances de cette nouvelle approche. Un taux de convergence
optimal est obtenu pour une direction de recherche k = 1.10+13 Pa.m−1 et un coefficient de stabi-
lisation ε = 1.10+11 Pa.m−1 . Les efforts de contact T projetés dans le plan (X ;Z) et le glissement
−
local [WT ] = (W+ T − WT ) sont représentés pour chaque valeur de µΓ dans les figures 2.22, 2.23 et
2.24.
Dans le cas adhérent (µΓ = 1) ou le cas du glissement partiel (µΓ = 0.5), les approximations de
la solution du problème calculées avec les méthodes SLA et NLLA sont stables et très semblables.
On n’observe pas d’oscillations numériques de la solution de contact. On remarque cependant que
la solution locale SLA est perturbée dans le voisinage des fronts de la fissure alors que la solution
NLLA est parfaitement régulière. Par ailleurs la méthode NLLA permet d’obtenir des temps de
calcul inférieurs de 70% par rapport ceux obtenus avec la méthode SLA (cf. Table 2.1).
Dans le cas du glissement total (µΓ = 0), la solution de contact SLA est fortement oscillante
(cf. figure 2.24.a). Ces instabilités numériques introduisent des pics de contraintes sur les faces
de la fissure et viennent perturber la solution globale du problème. En conséquence, la méthode
converge très lentement vers une solution erronée du problème. Par contre, la méthode NLLA
fournit une solution de contact régulière et une approximation précise de la solution du problème
(cf. figure 2.24.b). Elle permet par ailleurs d’augmenter considérablement le taux de convergence
(cf. tableau 2.1). L’Évolution de l’erreur local η au cours des itérations pour les trois cas considérés
avec les modèles SLA-X-FEM et NLLA-X-FEM est représentée dans la figure 2.25. On remarque
73
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
Modèle X-FEM-NLLA 87 86 86
TABLE 2.1: Nombre d’itérations à convergence et gain de temps CPU procuré par le
méthode NLLA par rapport à la méthode SLA.
clairement un palier de convergence pour le modèle SLA-X-FEM qui ne permet pas d’atteindre
des précisions élevées. En revanche, la convergence du modèle NLLA-X-FEM est rapide et très
régulière. Elle permet d’atteindre des précisions élevées à moindre coût numérique. Avec cette
méthode stable, on remarque que le nombre d’itérations à convergence est indépendant de la va-
leur du coefficient de frottement, c’est-à-dire de la nature des conditions de contact à l’interface
(adhérence, glissement partiel ou glissement total).
Ce premier cas-test tridimensionnel a permis de démontrer les bonnes propriétés de stabilité du
modèle multi-échelle X-FEM associé au solveur NLLA et son aptitude à "capturer" précisément
74
Exemples tridimensionnels et validation du modèle
différentes types de solutions de contact avec frottement à l’interface. Ceci est possible grâce à
l’utilisation d’une formulation faible mixte couplée à une description multi-échelle du problème.
Une approximation précise de la solution de contact entre les faces de la fissure est un pré-requis
essentiel pour le calcul précis des facteurs d’intensité des contraintes le long du front. En effet,
les champs des efforts sont pris en compte dans le calcul de l’intégrale d’interaction (cf. section
[Link]). La méthode d’intégration numérique en volume et sur la surface de contact utilisée dans
le cadre du modèle X-FEM proposé dans ce mémoire et le calcul des FIC est présentée ci-après.
75
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
F IGURE 2.25: Évolution de l’erreur local η au cours des itérations de la résolution non
linéaire pour les trois cas considérés avec les modèles SLA-X-FEM et NLLA-X-FEM.
76
Exemples tridimensionnels et validation du modèle
Les facteurs d’intensité de contraintes (FIC) KI , KII et KIII sont calculés le long du front de fis-
sure à l’aide de l’intégrale d’interaction I en 3D présentée dans la section [Link]. La discrétisation
du front de fissure est basée sur le maillage de la structure : 1 point de calcul Pi par élément fini
traversé par le front. Le domaine d’intégration D est un parallélépipède de dimensions r1 × r2 × r3
(suivant les directions e1 , e2 , e3 du repère local associé à la fissure) centré sur le point Pi et doté de
8 × 216 points d’intégration (cf. figure 2.26). L’intégrale d’interaction sur le volume D et sur les
faces de la fissure s’écrit :
Z Z Z
aux
I ℜ,aux = A1 j, j q1 dV + A1 j q1, j dV + tℜ
2 u1,1 q1 dS (2.67)
D D Γ+ ∪Γ−
avec
aux aux ℜ
A1 j = Wlℜ,aux δ1 j − σℜ u
1 j 1,1 − σ1j u1,1 (2.68)
où le champ d’extension virtuel q est colinéaire à l’orientation de la fissure et vérifie les conditions
(1.29). L’extension virtuelle et le champ d’extension virtuelle suivent une loi en cos2 (cf. figure
2.27.a) :
δl(s) = cos2 π2 xr33
(2.69)
q = cos2 π2 xr11 cos2 π2 xr22 cos2 π2 xr33 e1
où x1 , x2 et x3 sont les coordonnées associées au repère local du front de fissure (cf. figure 1.10)
issues du formalisme des fonctions de niveau. Les conditions (1.29) sont respectées lorsque le
parallélépipède est entièrement contenu dans le volume de la structure. Lorsque D est traversé
par un bord de la structure, δl(s) et q sont multipliés par une fonction rampe valant 0 hors de la
structure et atteignant linéairement la valeur 1 au centre du domaine d’intégration [RAN 08a] (cf.
figure 2.27.b).
77
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
F IGURE 2.28: Facteurs d’intensité des contraintes en mode II et III calculés le long des
deux fronts de fissures pour µΓ = 1, 0.5 et 0.
Pour µΓ = 1 et 0.5, le glissement relatif entre les faces de la fissure est nul ou très faible. Les
valeurs de KII sont très proches. Pour µΓ = 0, le glissement à l’interface est beaucoup plus élevé.
Les valeurs de KII sont donc beaucoup plus importantes. Pour ce problème, les valeurs de KIII
calculés le long des fronts sont très faibles en comparaison avec les valeurs de KII et indépendantes
de la valeur du coefficient de frottement µΓ .
78
Exemples tridimensionnels et validation du modèle
79
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
80
Exemples tridimensionnels et validation du modèle
C - On utilise ici le modèle X-FEM à deux échelles. Un maillage relativement grossier de 3456
tétraèdres (R = 5.5) identique au cas B est utilisé (cf. figure 2.30.b). La discrétisation de la
fissure est raffinée indépendamment du maillage en volume (cf. figure 2.2.b). On choisit
81
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
une taille critique pour les arrêtes des éléments d’interface lc = 1.8 mm, ce qui conduit à
une distribution de 2076 paires de points de d’intégrations à l’interface (cf. figure 2.31.c).
82
Exemples tridimensionnels et validation du modèle
Dans le cas B, la discrétisation de l’interface est trop grossière et ne permet pas de décrire avec
suffisamment de précison le problème de contact entre les faces de la fissure. En conséquence, la
frontière entre la zone d’ouverture et la zone de contact n’est pas localisée précisément (cf. figure
2.33.b) : x=0.0167 mm (erreur ∼ = 10%). En effet, le degré de précision est restreint par la finesse
de la discrétisation d’éléments d’interface dépendante du maillage de la structure.
On peut donc conclure que le modèle X-FEM à deux échelles combiné au solveur NLLA
permet de résoudre précisément et à moindre coût numérique les problèmes de rupture multi-
échelles bidimensionnels et tridimensionnels avec contact et frottement entre les faces de la fissure.
On peut noter par ailleurs que les solutions de contact ne présentent pas d’oscillations numériques
quelles que soient les discrétisations en volume et à l’interface utilisées. Ceci vient confirmer les
bonnes propriétés de stabilité du modèle X-FEM à deux échelles basé sur une formulation faible
à trois champs et combiné au solveur NLLA. Il fournit une solution stable et permet de capturer
précisément les conditions de contact indépendamment du problème de structure. Dans la section
suivante, un cas-test de validation est proposé. L’objectif est de comparer les résultats du modèle
X-FEM à deux échelles avec ceux d’un logiciel éléments finis industriel.
83
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
-0.0182 -0.0 136 -0.00909 -0.00455 0.0 0.00455 0.00909 0.0136 0.0182 0.0227 0.0271 0.0318
y
~x
-== = ====:: : : : : = =
J .79t+il0 7 7.59e+007
-0.0 182 -0.0136 -0.00909 -0.00455 0.0 o.o04.s.s o.oo9o9 o.otJ6 [Link] 0.0221 [Link] [Link]
y
3.86t+007 7.72e+007 [Lx
-0.0182 ·0.0 136 -0.00909 -0.00455 0.0 o.o04.s.s o.oo9o9 o.on6 o.OI82 o.o221 [Link] o.OJI8
y
3.8 h +00 7 7.6 l e+007 [Lx
F IGURE 2.33: Représentation des efforts d’interface sur le plan de fissure et tracé de la
frontière entre la zone de contact et la zone d’ouverture pour chaque configuration.
84
Exemples tridimensionnels et validation du modèle
F IGURE 2.35: (a) Maillage X-FEM en volume (3072 tétraèdres) et (b) discrétisation raf-
finée de l’interface (832 points d’intégration).
Premièrement, le modèle X-FEM à deux échelles est utilisé pour modéliser ce problème. Un
maillage en volume structuré et relativement grossier de 3072 tétraèdres est utilisé (cf. figure
85
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
2.21.a). La taille caractéristique des tétraèdres est de 5 mm. La discrétisation de l’interface est
raffinée à l’aide de la stratégie proposée dans la section 2.2. Une distribution de 832 points d’inté-
gration est obtenue (cf. figure 2.21.b).
Le problème est résolu à une précision de η = 10−4 en considérant le critère d’erreur locale
η (2.34). Le champ de déplacement global U amplifié (×120) et le champ des efforts d’interface
sont représentés figure 2.36. La fissure est complètement fermée. La solution locale de contact est
régulière et stable. On n’observe pas d’oscillations numériques.
Le problème est ensuite modélisé avec le logiciel éléments finis ANSYS. Ce modèle est pris
comme référence. Il faut noter l’importante différence entre les maillages utilisés. Le maillage
ANSYS est composé de 41 000 éléments finis. Un raffinement conséquent est appliqué le long de
la fissure et plus particulièrement au voisinage du front (cf. figures 2.37). La taille caractéristique
des éléments finis varie du millimètre au micron.
86
Exemples tridimensionnels et validation du modèle
2.36 et 2.37 respectivement. Ces deux résultats sont très proches : X-FEM : kUkmax = 16.1 µm et
kUkmax = 49.7 MPa ; ANSYS : kUkmax = 16.6 µm et kUkmax = 45.9 MPa. Les surfaces des iso-
valeurs de la norme du déplacement global sont quasi-identiques (cf. figures 2.36.a et 2.37.a). Cette
comparaison démontre la précision élevée de la solution obtenue avec le modèle X-FEM à deux
échelles pour un maillage de la structure relativement grossier et une description fine et intrinsèque
de l’interface. En effet, à précision donnée, la solution du problème requiert un maillage en volume
environ dix fois plus fin avec la MEF. Le modèle proposé permet de capturer les non linéarités
locales de contact avec frottement à une échelle pertinente et de réduire significativement les temps
de calcul.
Dans la section suivante, la stratégie en temps qui sera employée pour la simulation de la
propagation des fissures sous chargement cyclique est présentée.
87
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
Comme nous l’avons vu dans ce chapitre, la méthode de résolution non linéaire NLLA est
basée sur une version incrémentale en temps de méthode LATIN initialement proposée par Cham-
paney [CHA 99]. Le cycle de chargement borné par les piquets de temps Tm et Tm+1 est discrétisé
en n piquets de temps tn et autant de pas de chargement (cf. figure 2.38). On définit ainsi l’échelle
temporelle fine, appelée échelle nř1, et l’échelle temporelle grossière , appelée échelle nř2.
L’histoire des conditions de contact entre les faces de la fissure à une influence très importante
sur le comportement de la fissure. En effet, les problèmes non linéaires dépendent du temps. Une
modélisation rigoureuse des problèmes de rupture avec contact et frottement entre les faces de la
fissure requiert donc une formulation quasi-statique. En conséquence, au pas de temps tn+1 , les
champs locaux de contact W et T sont initialisés avec les valeurs calculées au pas de temps tn . Au
premier pas de temps, on considère que ces champs sont nuls.
Pour chaque pas de chargement, on calcule la solution Xn et on calcule les facteurs d’intensité
des contraintes Ki (tn ). À la fin de chaque cycle, c’est-à-dire au pas de temps Tm+1 de l’échelle
da
temporelle nř2, on calcule alors l’accroissement par cycle dN et la direction de propagation θ de
chaque fissure à l’aide de la loi et du critère de propagation adaptés au matériau et à la nature de
la sollicitation du problème. On définit alors la nouvelle géométrie des fissures pour un nombre
de cycle donné. Les différents pas de propagation de la fissure constituent la troisième échelle de
temps du problème, appelée échelle nř3.
88
Modélisation quasi-statique et simulation multi-échelle en temps de la propagation des fissures
Bilan :
Le modèle X-FEM multi-échelle avec contact et frottement interfacial proposé dans ce chapitre
satisfait les pré-requis énoncés dans le chapitre 1. Le problème de rupture est modélisé par la
méthode des éléments finis étendus. L’utilisation d’une formulation faible mixte permet une des-
cription multi-échelle du problème de contact avec frottement entre les faces de la fissure. En effet,
la discrétisation de l’interface est indépendante du maillage de la structure et adaptée à l’échelle
des non linéarités locales. Plusieurs méthodes de résolution du problème non linéaire ont été pro-
posées. Des exemples ont mis en évidence les problèmes d’instabilité de la solution de contact.
La méthode de stabilisation du solveur basé sur la méthode LATIN permet de s’affranchir des
éventuelles oscillations numériques de la solution de contact et améliore grandement le taux de
convergence de la méthode itérative. Plusieurs exemples bidimensionnels et tridimensionnels ont
montré les performances du modèle proposé en ce qui concerne la finesse et la stabilité de la solu-
tion de contact. La comparaison avec un modèle de référence a permis de valider quantitativement
les résultats obtenus.
L’étape suivante consiste à intégrer ce modèle dans une stratégie expérimentale et numérique
multi-modèle de prédiction de durée de vie des composants de structures mécaniques soumis aux
problèmes de fretting fatigue.
89
2. Modèle multi-échelle dédié aux problèmes de rupture avec contact et frottement entre
les faces de la fissure
90
Chapitre 3
Sommaire
91
3.1 Stratégie de prédiction de durée de vie dédiée aux problèmes de fretting
fatigue . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
3.1.1 Problème expérimental de fretting fatigue de référence . . . . . . . 95
3.1.2 Résolution du problème de contact deux-corps . . . . . . . . . . . 96
3.1.3 Prédiction du risque d’amorçage de fissures de fretting . . . . . . . 100
3.1.4 Modélisation X-FEM avec contact et frottement interfacial et vali-
dation du calcul des facteurs d’intensité des contraintes en 2D . . . 101
3.1.5 Étude paramétrique de la propagation des fissures sous chargement
de fretting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
3.2 Étude expérimentale de la propagation des fissures de fretting en confi-
guration cylindre/plan . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
3.3 Simulation numérique X-FEM 2D de la propagation de fissures expé-
rimentale sous chargement de fretting . . . . . . . . . . . . . . . . . . . 124
3.3.1 Analyse du problème expérimental : Résolution du contact cy-
lindre/plan et prédiction du risque d’amorçage des fissures . . . . . 124
3.3.2 Construction d’une loi de propagation expérimentale dédiée aux
problèmes de fretting fatigue . . . . . . . . . . . . . . . . . . . . . 125
3.3.3 Simulation numérique X-FEM de la propagation des fissures sous
chargement de fretting fatigue . . . . . . . . . . . . . . . . . . . . 127
92
Stratégie de prédiction de durée de vie dédiée aux problèmes de fretting fatigue
Les problèmes de fretting sont généralement observés à la jonction entre les composants de
structures soumis à des vibrations ou à un chargement cyclique. Les déplacements relatifs ré-
sultants localisés à l’interface provoquent un endommagement par usure des surfaces ou par fis-
suration. La nature de cet endommagement et sa sévérité dépendent d’un très grand nombre de
paramètres (géométries macroscopiques, rugosité des surfaces, matériaux, chargement, physico-
chimie,...) qui conditionnent les sollicitations locales cycliques et leur évolution temporelle. La
simulation des problèmes de fretting requiert donc une connaissance précise des conditions de
contact entre les deux solides et par conséquent, une estimation du coefficient de frottement à
l’interface et de ses variations.
L’objectif de ce chapitre est de proposer une stratégie globale intégrant les avantages de plu-
sieurs modèles afin de prédire la durée de vie d’éprouvettes d’essais soumises à un chargement
de fretting fatigue. Une stratégie a été mise en place, combinant les avantages de trois modèles,
chacun étant dédié à l’un de ces objectifs :
i. La modélisation du contact deux-corps et le calcul des conditions de contact à l’interface.
ii. La prédiction du risque d’amorçage de fissures à partir du critère de Dang Van.
iii. La simulation X-FEM de la propagation des fissures de fatigue sous chargement de fretting.
Cette stratégie et les modèles qui la composent sont présentés dans la première partie de ce
chapitre. On s’intéressera en premier lieu à un cas expérimental issu de la littérature. La modéli-
sation X-FEM du problème de fissuration est validée en comparant les résultats obtenus avec ceux
d’un modèle de référence basé sur la méthode des distributions continues des dislocations. Puis,
une étude paramétrique est réalisée dans le but de caractériser la propagation des fissures de fret-
ting et d’évaluer l’influence du coefficient de frottement entre les faces de la fissure (difficilement
mesurable expérimentalement) et d’une contrainte statique appliquée à l’éprouvette. Cette stra-
tégie complète est ensuite appliquée à la simulation et à l’analyse d’essais cylindre/plan réalisés
dans le cadre de ce travail.
93
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
1. Les données des essais (boucles de fretting) sont analysées afin de déterminer les para-
mètres de chargement expérimentaux. Cette première étape permet de déterminer le mode
d’endommagement à partir du régime de glissement observé.
2. Le problème contact normal et tangentiel quasi-statique à l’interface des deux solides est ré-
solu afin de quantifier l’aire de contact, les zones d’adhérence et de glissement et de calculer
les distributions des pressions normales et des contraintes de cisaillement surfaciques.
3. Le risque d’amorçage des fissures est ensuite calculé à l’aide du critère de fatigue multiaxial
de Dang Van [DAN 93]. Ce critère de plan critique est basé sur le champ des contraintes
dans la zone sous le contact calculé à partir du résultat de l’étape 2. Il permet également de
déterminer les sites et les directions d’amorçage des fissures de fretting [ROB 06]. La géo-
métrie des fissures peut éventuellement être extraite à partir d’observations expérimentales
(étape 3bis).
4. À partir des données de chargement en surface (étape 2) et de la géométrie initiale des
fissures (étape 3), le problème de fatigue sous chargement de fretting fatigue est modélisé à
l’aide du modèle X-FEM avec contact et frottement interfacial proposé dans ce mémoire.
5. Les facteurs d’intensité des contraintes sont calculés en pointe de fissure au cours du cycle
de chargement à partir de la solution du problème de rupture.
6. La propagation des fissures de fretting est simulée à l’aide des critères de propagation en
fatigue multiaxiale dédiés aux chargement non proportionnels adaptés au problème de fret-
ting.
Dans les paragraphes suivants, un essai issu de la littérature est considéré. L’objectif est de
présenter les différents modèles et techniques utilisées pour accomplir les différentes étapes de la
94
Stratégie de prédiction de durée de vie dédiée aux problèmes de fretting fatigue
simulation de cet essai, puis de valider la stratégie proposée. Cette démarche est conduite en 2D.
Dureté Hv 160
Expérimentalement, on observe l’amorçage d’un réseau de fissures courtes dans la zone an-
nulaire de micro-glissement en périphérie de la zone d’adhérence. Ces fissures forment un angle
variable de 15◦ à 35◦ en volume avec la surface de contact en fonction de leur position par rapport
au centre du contact (cf. figure 3.2). La propagation de la majorité des fissures stoppe et seules
deux fissures principales se propagent selon une trajectoire semi elliptique en surface et avec un
angle de 29◦ par rapport à la surface de contact. Elles bifurquent ensuite progressivement pour
former un angle de 67◦ avec la surface.
95
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
F IGURE 3.2: Essai de fretting fatigue réalisé par Dubourg et al. [DUB 00] : sites d’amor-
çage et profile des fissures observées.
Cet essai est modélisé par une approche 2D. Son analyse est donc conduite dans le plan médian
(y=0).
n n n
uz,i = ∑ ai j · p j + ∑ bi j · q j uz,i = ∑ ai j · p j
j=1 j=1 j=1
n n =⇒ simpli f ication =⇒ n (3.1)
ux,i = ∑ ci j · q j + ∑ di j · p j ux,i = ∑ ci j · q j
j=1 j=1 j=1
96
Stratégie de prédiction de durée de vie dédiée aux problèmes de fretting fatigue
F IGURE 3.3: Zone d’intérêt potentielle pour la résolution du contact deux-corps discréti-
sée en N cellules de longueur ∆x et de coordonnée xi .
Une analyse 2D des essais sphère/plan est conduite en première approche. Un problème équi-
valent cylindre/plan est défini. Les travaux [NOW 90, HIL 94] ont montré que l’amorçage et la
propagation des fissures sont étroitement liés à la dimension du contact. La pression maximum de
Hertz P0 = 150 MPa et la dimension de l’aire de contact a = 1.785 mm sont identiques en 2D et en
97
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
3D. Les paramètres de charge linéiques (P et Q) et le rayon du cylindre RC sont alors déterminés
à partir cette contrainte. Le rapport c/a entre les rayons des zones circulaires d’adhérence et de
contact est égal à 0.408. Les nouveaux paramètres 2D pour la simulation de l’essai de fretting
fatigue sont donnés dans le tableau 3.2.
RC P Qmax F µ a c/a
98
Stratégie de prédiction de durée de vie dédiée aux problèmes de fretting fatigue
F IGURE 3.5: Représentation des distributions des contraintes normales p(x) et de cisaille-
ment q(x,tn ) (MPa) sur l’interface de contact cylindre/plan pour une cycle de fretting com-
plet (9 pas de temps sur 80 sont représentés).
Le chargement est considéré ici comme monotone. On fait l’hypothèse que le régime de glis-
sement partiel engendre une usure de surface très faible ne venant pas perturber la géométrie du
contact, ni la valeur du coefficient de frottement, ni le champ des contraintes normales et tangen-
tielles. Une solution précise du cycle de charge à l’interface des deux corps est obtenue à l’échelle
de ce contact. C’est un des points clé de la stratégie de prédiction de durée de vie proposée. En
99
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
avec
6t−1 − 3 f−1
α= β = 2t−1 (3.7)
f−1
où σi et σ j sont les contraintes principales et f−1 et t−1 sont respectivement les limites de fatigue
en traction et en torsion alternée pour un rapport R = -1.
Il est possible de tracer la distribution du risque d’amorçage en volume. L’amorçage de fissures
de fatigue a lieu pour les couples (tan(t), σ(t)) tels que d > 1. Le plan correspondant indique le
lieu et la direction initiales des fissures. Ce critère a été appliqué avec succès pour prédire le risque
d’amorçage de fissures sous sollicitations cycliques multiaxiales [DAN 93], notamment dans le
cas du fretting [FOU 96, FOU 02, BAI 10].
À partir du champ des contraintes en volume sous le contact calculé dans la section 3.1.2,
le risque d’amorçage de fissures est déterminé à l’aide du critère de Dang Van. La figure 3.7
représente la distribution du risque d’amorçage d et l’angle d’amorçage θ dans le domaine d’intérêt
en fonction de la position (x/a, z/a). x/a = ±1 correspond aux extrémités de la zone de contact. Le
risque d est supérieur à 1 au voisinage des deux zones de micro-glissement située de x/a = -1.12 à
-0.7 et x/a = 0.7 to 1.12. Il diminue progressivement vers la frontière de la zone d’adhérence. Les
valeurs maximales sont situées aux positions x1 /a = -0.92 et x2 /a = 0.92 pour un angle θ = ±29◦
avec la surface. Ce résultat est en bon accord avec les observations expérimentales [DUB 00].
Remarque 18 :La distribution de l’angle θ en profondeur donne une idée du parcours des fissures
en fatigue. Cette indication est cependant limitée à une faible profondeur. En effet, la présence
des fissures après l’amorçage vient modifier le champ des contraintes cycliques et le parcours
initialement prédit.
100
Stratégie de prédiction de durée de vie dédiée aux problèmes de fretting fatigue
Dans le cadre de la stratégie proposée, le critère de Dang Van permet de déterminer la géomé-
trie initiale des fissures de fretting. Ce résultat ainsi que les distributions des contraintes en surface
déterminées dans le section 3.1.2 sont utilisés dans la section suivante comme données d’entrée
pour la modélisation X-FEM du problème de fissuration avec contact et frottement.
101
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
102
Stratégie de prédiction de durée de vie dédiée aux problèmes de fretting fatigue
F IGURE 3.10: Champ des efforts T entre les faces de la fissure 2 au pas de temps 1 (Q =
Qmax .
de l’intégrale d’interaction I présentée dans la section [Link]. Le domaine D utilisé est un carré
de coté 5 µm indépendant du maillage et centré sur la pointe de la fissure (cf. figure 3.11).
L’objectif ici est de valider le calcul des FIC Les résultats obtenus avec le modèle X-FEM sont
103
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
104
Stratégie de prédiction de durée de vie dédiée aux problèmes de fretting fatigue
Modèle de référence :
Il s’agit du modèle développé par Dubourg [DUB 89] basé sur la méthode des distributions conti-
nues de dislocations pour la modélisation des zones de discontinuités de déplacement à l’inter-
face d’une fissure et la résolution du contact unilatéral avec frottement. La formulation initiale
[DUN 64] utilisée par Comninou [COM 77] et et de Hills et al. [HIL 85] a été modifiée et géné-
ralisée pour être intégrée dans un algorithme de résolution du contact unilatéral avec frottement
[DUB 89, DUB 92a]. L’état de contact avec frottement interfacial est calculé à l’échelle locale.
Une formulation incrémentale permet de prendre en compte les effets d’hystérésis liés à l’histoire
des conditions de contact avec frottement interfacial.
Ce modèle est très précis pour des temps de calcul très réduits. Il est parfaitement adapté
à la problèmatique de la fatigue tribologique – pour laquelle l’hypothèse de massif semi-infini,
comme nous l’avons vu, est justifiée – qui demande une description très fine des cycles de charge
complexes et des gradients sévères de contrainte et de déformation générés dans des zones confi-
nées sous le contact, et permet de déterminer précisément l’évolution des non linéarités de contact
avec frottement à l’interface des fissures. La méthode permet de modéliser des problèmes cy-
cliques complexes (i.e. de fretting, de roulement, des séquences d’ouverture-fermeture-adhérence-
glissement) sous chargement multiaxial non proportionnel et l’interaction de plusieurs fissures
[DUB 92b]. Ce modèle a été validé par comparaison avec des simulations expérimentales pour
des problèmes de fretting [LAM 96, BAI 02a] et de roulement [BAI 02b]. Il permet donc de vali-
der la modélisation X-FEM d’un essai de fretting fatigue.
Chaque fissure du problème est discrétisée ici avec 77 points d’intégration. Le cycle de fretting
est discrétisé en 201 pas de temps. Les FIC KI et KII sont calculés à chaque pas de temps tn
directement à partir du champ de déplacement, i.e. ouverture et glissement entre les faces de la
fissure [KRE 75] (cf. équation 1.3). Les figures 3.12, 3.13 et 3.14 représentent l’évolution des FIC
KI et KII au cours du cycle de fretting en fonction de l’effort tangentiel Q calculés à la pointe des
fissures 1 et 2 avec le modèle X-FEM et le modèle de référence.
Les résultats obtenus avec le modèle X-FEM avec contact et frottement interfacial concordent
très bien avec ceux du modèle de référence. Bien que la solution du problème de contact pour
chaque modèle soit identique, les résultats du calcul des FIC est légèrement différent aux pas
de temps pour lesquels on observe l’adhérence des faces des fissures, i.e. au début des phase
de décharge (Q = +Qmax ) et de charge (Q = -Qmax ). En effet, le calcul des FIC avec le modèle
de référence repose sur une approche en déplacements. Lorsqu’il y a adhérence, l’évolution du
glissement relatif entre les faces de la fissure est nulle et les variations de KII et ∆KII sont nulles.
En revanche, le calcul des FIC avec le modèle X-FEM repose sur une approche énergétique en
contraintes et en déplacements (intégrale d’interaction) qui capture les variations des contraintes
lorsque le déplacement relatif est nul.
Bilan :
Le modèle X-FEM avec contact et frottement proposé dans ce mémoire et le calcul des facteurs
d’intensité des contraintes en mode mixte à l’aide de l’intégrale d’interaction en 2D sont validés.
Le solveur stabilisé NLLA et le critère d’erreur local (2.34) permettent d’optimiser la convergence
des problèmes normal et tangentiel en déplacements (w) et en contraintes (t) sur l’intégralité de
l’interface de la fissure. L’évolution des non linéarités de contact avec frottement à l’interface des
fissures et les FIC sont calculés précisément.
105
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
Ce modèle s’insère dans la stratégie multi-modèles pour la simulation numérique des pro-
blèmes de fretting fatigue proposée. Les avantages des différents modèles sont couplés dans une
approche expérimentale et numérique. À partir de résultats d’essais de fretting, le problème de
contact entre les corps en contact a été résolu à l’aide d’un modèle semi-analytique pour le calcul
de l’évolution des distributions de contraintes normales et tangentielles dans la zone de contact au
cours du cycle. Le risque d’amorçage, les sites et l’orientation initiale des fissures de fretting ont
été déterminés en utilisant le critère de Dang Van Ces résultats sont ensuite utilisés comme don-
nées d’entrée la modélisation X-FEM du problème de fissuration. Le modèle X-FEM avec contact
et frottement proposé dans ce mémoire et le calcul des facteurs d’intensité des contraintes en mode
mixte à l’aide de l’intégrale d’interaction en 2D ont été validés par comparaison avec un modèle
de référence basé sur la méthode des distributions continues de dislocations.
La simulation de la propagation des fissures requiert une loi de propagation en vitesse et l’uti-
lisation d’un critère de direction de propagation en fatigue multiaxiale non proportionnelle adapté
au problèmes de fretting. Dans la section suivante, les directions de propagation calculées avec
différents critères sont comparés avec celles mesurées expérimentalement. Nous étudierons en-
suite l’influence du coefficient de frottement interfacial et d’une contrainte de traction volumique
additionnelle sur le comportement des fissures.
106
Stratégie de prédiction de durée de vie dédiée aux problèmes de fretting fatigue
107
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
108
Stratégie de prédiction de durée de vie dédiée aux problèmes de fretting fatigue
KII
F IGURE 3.15: Variation de l’angle de phase ΨM = arctan KI au cours du cycle de
fretting.
Les figures 3.12, 3.13 et 3.14 montrent que la fissure 1 est sollicitée en mode mixte I et II,
puis en mode II pur et de nouveau en mode mixte. La fissure 2 est sollicité en mode II pur, puis
en mode mixte et de nouveau en mode II pur. Il est possible de calculer l’angle de phase ΨM (3.8)
qui caractérise la mixité de la sollicitation au front [SUR 98] et de tracer son évolution au cours
du cycle (cf. figure 3.15).
KII
ΨM = arctan (3.8)
KI
ΨM = 0 indique une sollicitation en mode I pur et ΨM = π/2 caractérise une sollicitation en
mode II pur. La variation de ΨM au cours du cycle de fretting indique que les deux fissures sont
sollicitées en mode mixte non proportionnel. Les critères de direction de propagation en fatigue
multiaxiale non proportionnelle présentés dans la section [Link] sont donc utilisés dans le cadre
de la stratégie de simulation numérique proposée dans ce mémoire de thèse. Ces critères sont
basés sur l’évolution de l’amplitude des facteurs d’intensité des contraintes et de la mixité des
sollicitations au cours du cycle. L’objectif du paragraphe suivant est de montrer l’influence du
choix de ce critère sur le trajet de propagation calculé.
109
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
1 et 2 est de 88 µm. D’après [BRO 89], une extension de 1% est préconisée. Un saut de 2 µm est
considéré ici, soit légèrement plus de 2%. L’angle de propagation δθ entre la fissure initiale et
l’extension calculée (positif dans le sens trigonométrique) après trois pas de propagation de 2 µm
est répertorié dans le tableau 3.3 en fonction du critère utilisé.
TABLE 3.3: Direction de bifurcation δθ des fissures de fretting par rapport à la direc-
tion initiale (θ = 29◦ ) en fonction du critère de propagation en fatigue multiaxiale non
proportionnelle utilisé.
Des résultats différents sont obtenus en fonction du critère utilisé. Le critère de Hourlier 2
(∆k1∗ (θ)max ) fournit les résultats les plus proches des observations expérimentales [BAI 02a]. En
effet, il est adapté aux problèmes de fatigue multiaxiale avec un fort rapport de mixité.
110
Stratégie de prédiction de durée de vie dédiée aux problèmes de fretting fatigue
frottement µΓ allant de 0.2 à 1.2. Les figures 3.17 et 3.18 représentent les variations de KII au cours
du cycles de fretting pour 6 valeurs différentes de µΓ .
Les variations des facteurs d’intensité des contraintes en mode I (KI ) sont indépendantes de la
valeur de µΓ . En revanche, il influence fortement les valeurs des FIC en mode II. Plus le coefficient
de frottement est élevé, plus l’amplitude ∆KII = KIImax − KIImin est faible. En effet, le frottement
réduit le glissement entre les faces de la fissure. Pour la fissure 1, KIImax est indépendant de la va-
leur de µΓ alors que KIImin augmente lorsque µΓ augmente. En effet, au début du cycle : Q = +Qmax ,
la fissure 1 est complètement ouverte et le frottement n’a aucune influence sur la sollicitation en
mode II et la valeur de KII . Pendant la phase de décharge (t0 à t40 ), les faces de la fissure 1 sont
en contact et le frottement restreint le glissement. Les même tendances sont observées pendant
la phase de charge (t40 à t80 ). Pour la fissure 2, KIImin est constant et KIImax diminue lorsque le
frottement augmente.
La mixité des sollicitations est indépendante du coefficient de frottement. En effet, les valeurs
de KII sont quasi constantes lorsque KI est non nul. Cependant, la prédominance du mode II
décroît. Le rapport ∆KI /∆KII varie de 0.42 à 0.89 lorsque µΓ varie de 0.2 à 1.2. Les critères de
propagation en fatigue multiaxiale non proportionnelle tels que le critère ∆k1∗ (θ)max dépendent
fortement de ce rapport. En conséquence, l’augmentation du coefficient de frottement entre les
faces de la fissure accroît l’angle de bifurcation de la fissure 2 de 45◦ à 55◦ par rapport à la
direction initiale (θ = 29◦ ) (cf. tableau 3.4).
111
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
Pour les fissures 1 et 2, l’amplitude de variation ∆KI = KImax − KImin du facteur d’intensité des
112
Stratégie de prédiction de durée de vie dédiée aux problèmes de fretting fatigue
contraintes en mode I augmente lorsque σt augmente (cf. figure 3.19). Les valeurs de KII au cours
113
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
du cycle augmentent lorsque σt augmente (cf. figure 3.20). En effet, une contrainte de traction
(σt >0) favorise l’ouverture des fissures et réduit le frottement à l’interface. Le glissement local
est donc plus élevé. En revanche, l’amplitude ∆KII = KIImax − KIImin est indépendante de la valeur
de σt . En conséquence, le rapport de mixité ∆KI /∆KII varie de 0,52 à 1,22 lorsque σt varie de
-115 MPa (-D/2) à 115 MPa (+D/2). La prédominance du mode II décroît lorsque σt augmente
(cf. figure 3.21). Cependant, cette contrainte de traction accroît l’angle de bifurcation de la fissure
2 : δθ varie de 49◦ pour σt = -115 MPa à 60◦ pour σt = +115 MPa (cf. figure 3.5). En effet, elle
favorise la propagation des fissures dans la direction perpendiculaire à la direction de la charge.
Le chargement en surface et la contrainte volumique ont une influence différente sur le com-
portement des fissures :
- Le chargement de fretting constitue la sollicitation en fatigue responsable de l’amorçage et
de la propagation des fissures dans une zone confinée sous le contact. Son action décroît
avec la profondeur. La croissance de la fissure peut éventuellement stopper.
114
Stratégie de prédiction de durée de vie dédiée aux problèmes de fretting fatigue
Dans la suite de ce chapitre, des essais de fretting cylindre/plan réalisés dans le cadre de ce
travail sont présentés. L’objectif est d’observer l’évolution des fissures (longueur, orientation) en
115
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
116
Stratégie de prédiction de durée de vie dédiée aux problèmes de fretting fatigue
117
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
118
Étude expérimentale de la propagation des fissures de fretting en configuration cylindre/plan
Le dispositif d’essai utilisé est schématisé dans la figure 3.26. Il comporte une machine de
traction servo-hydraulique MTS à deux vérins. Une contrainte statique de traction σt est appliquée
à l’éprouvette à l’aide du premier vérin. La section de l’éprouvette dans la zone d’intérêt est 8 mm
×10 mm (cf. figure 3.27). Un poinçon cylindrique de rayon R = 80 mm et de largeur l = 8 mm est
mis en contact avec le plan rectifié de l’éprouvette par une force P = 8 000 N. Un roulement est
mis en opposition du poinçon afin de limiter les sollicitations de flexion sur l’éprouvette. L’état
des surfaces en contact est rectifié (Ra = 0.4). Le matériau étudié qui constitue l’éprouvette plane
est un acier 35NCD16 traité thermiquement pour atteindre une limite à rupture σR = 1270MPa. Le
matériau du poinçon est un acier 100C6 durci par traitement thermique pour atteindre une limite
élastique σe0.2 = 1500MPa. Les propriétés mécaniques de ces deux matériaux sont données dans
les tableaux 3.6 et 3.7 respectivement.
L’ensemble poinçon-roulement est mis en mouvement cyclique par le second vérin à une fré-
quence F = 13 Hz. L’amplitude du déplacement relatif δ entre l’éprouvette et le poinçon est mesuré
à l’aide d’un capteur laser. L’effort tangentiel Q est déterminé au cours de l’essai par une mesure
différentielle entre deux cellules de force liées à chaque extrémité de l’éprouvette. Les paramètres
de charge sont choisis afin de se placer en régime de glissement partiel. La valeur de δ est asservie
telle que Q = 4000N pendant toute la durée de l’essai. Les paramètres d’essais sont répertoriés
dans le tableau 3.8.
119
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
200 GPa 0.3 1270 MPa 1127 MPa 400 460 MPa 590 MPa 400 MPa
E ν σe0.2 Dureté Hv
Le mode d’endommagement principal est la fissuration avec une usure très faible de la surface
de contact. Les traces de fretting observées sont composées d’une zone d’adhérence centrale de
largeur 2c et d’une zone de micro-glissement en périphérie de la zone de contact de largeur 2a
(cf. figure 3.28). Les valeurs de a = 1.02 mm et de c = 0.62 mm sont quasi constantes quelle que
soit la contrainte de traction σt appliquée à l’éprouvette. Deux fissures principales s’amorcent en
périphérie symétriquement par rapport au centre de la surface de contact.
120
Étude expérimentale de la propagation des fissures de fretting en configuration cylindre/plan
RC P Qmax δ F
F IGURE 3.28: Trace de fretting obtenue après un essai de 200 000 cycles avec une
contrainte de traction statique σt = 180 MPa. La zone d’adhérence (c ∼ = 0.62 mm) se
distingue de la zone de micro-glissement caractéristique du régime de glissement partiel.
L’usure de la surface de contact est très faible et la taille de la zone de contact (a ∼
= 1.02
mm) est quasi constante quels que soient la contrainte de traction appliquée et le nombre
de cycles de l’essai.
L’observation in-situ de la propagation de ces fissures dans les matériaux métalliques par une
technique non destructive telle que la micro-tomographie demeure complexe. En effet, cette tech-
nique requiert notamment de réaliser les essais à l’intérieur du dispositif d’observation et d’effec-
tuer des rotations complètes de l’échantillon [BUF 05]. De plus, elle n’est pas applicable à tous les
matériaux. En effet, la présence de précipités ou d’inclusions dans la matière permet, par une tech-
nique de corrélation d’images, de capturer les champs de déplacement des configurations initiale
et déformée et de reconstituer la géométrie précise des fissures en 3D [RAN 10]. En l’absence de
tels marqueurs, il est difficile d’observer les fissures de reconstituer précisément le scénario. Il est
donc choisi ici d’analyser les essais conduits à différents nombres de cycles afin de reconstruire
l’évolution des fissures au cours du temps. Après chaque essai, l’éprouvette post-mortem est dé-
coupée suivant un plan perpendiculaire à la trace de fretting (cf. figure 3.29). L’observation de
plusieurs coupes d’un même échantillon donne des résultats similaires et permet de valider cette
analyse. Une technique de numérisation d’image permet d’identifier précisément la géométrie de
la fissure, i.e. longueur et angle de propagation.
Des essais ont été réalisés pour six valeurs de la contrainte de traction statique σt allant de 0 à
280 MPa. Les autres paramètres d’essais (R, P, Q et F) sont constants. On obtient alors six courbes
121
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
de propagation, i.e. longueur des fissures a en fonction du nombre de cycles N, représentées figure
3.30. Ces résultats permettent d’étudier le taux de croissance des fissures en fonction des para-
mètres de charge et de la profondeur sous la surface. Pour des longueurs de fissure a inférieures
à 1.5 mm, plus la contrainte de traction statique est importante, plus le taux de croissance des
fissures da/dN est élevé. Il varie de 8.5e-03 µm/cycle pour σt = 0 MPa à 9.8e-03 µm/cycle pour
σt = 280 MPa. Le comportement des fissures est principalement influencé par le chargement de
fretting en surface. Il correspond au stade II de la propagation des fissures (cf. figure 1.8 et section
[Link]). La contrainte de traction additionnelle intensifie les sollicitations de fatigue.
122
Étude expérimentale de la propagation des fissures de fretting en configuration cylindre/plan
Cette étude expérimentale a permis d’observer les géométries des fissures de fretting à diffé-
rents stades de leur propagation. Elle fournit une estimation de la longueur et de l’orientation des
fissures en fonction du nombre de cycles de charge et de la contrainte de traction statique appli-
quée à l’éprouvette. Ces résultats peuvent donc être mis en relation avec les résultats de simulation
numérique de ces essais afin de construire une loi de propagation.
Dans le cadre de la stratégie de prédiction de durée de vie proposée dans ce mémoire, ces
résultats fournissent également la géométrie des fissures utilisée comme donnée d’entrée pour la
simulation X-FEM. Les données de chargement du problème sont déterminées dans la section
suivante par la résolution du problème de contact deux-corps.
123
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
124
Simulation numérique X-FEM 2D de la propagation de fissures expérimentale sous chargement
de fretting
Les champs de contrainte et de déformation sont calculés dans la zone d’intérêt sous le contact
à partir du chargement des contraintes cycliques surfaciques. Le risque d’amorçage des fissures est
déterminé en utilisant le critère de Dan Van (cf. section 3.1.3). Le risque d’amorçage d et l’angle
d’amorçage θ dans le domaine d’intérêt pour une contrainte statique σt = 0 sont présentés figure
3.32. Le risque d obtenu est supérieur à 1 au voisinage des deux zones de micro-glissement situées
de x/a = -0.98 à -0.72 et x/a = 0.72 à 0.98. Les deux fissures principales s’amorcent en périphérie
de la zone de micro-glissement aux positions x1 = -0.95 mm et x2 = +0.95 pour lesquelles la valeur
du critère est maximale. L’angle d’amorçage calculé à ces positons est de ±76◦ par rapport à la
surface de contact. Ce résultat est en bon accord avec les observations expérimentales.
Les essais sont modélisés en 2D (cf. section 3.1.4) en appliquant la stratégie de simulation nu-
mérique proposée dans ce mémoire. Le chargement de fretting (i.e. les distributions de pressions
normales p(x) et de contraintes de cisaillement q(x,tn ) sur la surface de contact) déterminé dans
125
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
la section précédente est utilisé comme donnée d’entrée du modèle X-FEM à deux échelles avec
contact et frottement interfacial. Les géométries des fissures expérimentales obtenues par obser-
vation métallographique des coupes transversales des échantillons sont considérées. Le cycle de
fretting est discrétisé ici en 40 pas de temps. On considère la loi de contact unilatéral et un frotte-
ment de Coulomb µΓ = 0.3 entre les faces des fissures. Les facteurs d’intensité des contraintes KI
et KII sont calculés pour un cycle de chargement complet.
Remarque 21 : Ici, la valeur du coefficient de frottement µΓ entre les faces des fissures est in-
connue a priori. Elle est définie arbitrairement. Il est effectivement très difficile de l’évaluer ex-
périmentalement car les fissures se situent à l’interface de contact cylindre/plan et l’amplitude
du glissement interfacial ne peut être mesurée in-situ. Néanmoins, il a été montré dans la section
[Link] que la valeur de µΓ n’a pas d’influence sur ∆KI et ne fait varier que ∆KII pour un pro-
blème de fissuration sous chargement de fretting donné. La valeur du coefficient b de la loi de
propagation en vitesse a donc un effet correctif sur l’approximation de µΓ .
126
Simulation numérique X-FEM 2D de la propagation de fissures expérimentale sous chargement
de fretting
Les coefficients C = 1.3 · 1012 , m = 3.4 et b = 0.79 fournissent une bonne approximation nu-
mérique du taux d’accroissement expérimental des fissures da/dN quelle que soit la contrainte
de traction statique σt appliquée à l’éprouvette. La figure 3.33 q représente d’une part les points
expérimentaux obtenus à partir de la valeur numérique de ∆K = ∆KI2 + b · ∆KII2 et du taux d’ac-
croissement da/dN calculé à partir des résultats des essais, et d’autre part la courbe de la loi de
propagation (3.10) construite à partir des résultats de la simulation numérique X-FEM des essais.
La courbe de la loi de Paris expérimentale de l’acier 35NCD16 pour un chargement de fatigue
conventionnelle en mode I pur (R = 0.1) en milieu ambiant [BOU 89] est également tracée.
La loi de propagation en mode mixte proposée est en bon accord avec la majorité des points
expérimentaux. Cependant, pour les fissures dont la taille est supérieur à 1 500 µm, les points ex-
périmentaux sont plus proches de la courbe de la loi de propagation en fatigue conventionnelle en
mode I. Comme nous l’avons vu, l’influence du chargement de fatigue de contact décroît avec la
profondeur. Au delà d’une certaine profondeur, les fissures se propagent en mode I perpendiculai-
rement à la surface de contact sous l’influence majoritaire de la contrainte de traction σt , corres-
pondant au stade III de la propagation des fissures sous chargement de fretting (cf. figure 1.8). Le
chargement cyclique en surface influence peu la mixité des sollicitations en pointe de fissures à
ce stade. Le comportement des fissures se rapproche alors de celui en fatigue conventionnelle en
mode I.
Une loi expérimentale de propagation en vitesse des fissures sous chargement de fretting a
été construite à partir des résultats de la modélisation X-FEM 2D avec contact et frottement des
essais de fretting fatigue. Elle est en accord avec les lois issues de la littérature pour la fatigue
conventionnelle (uniaxiale). Deux comportements différents ont été identifiés en fonction de la
longueur des fissures. Le comportement des fissures de petite taille sous la zone de contact (stade
II, cf. figure 1.8) peut être décrit à l’aide de la loi de propagation en mode mixte sous chargement
de fretting proposée (3.10). Au delà de la zone d’influence des sollicitations tribologiques (stade
III), la fissure se propage selon une loi en fatigue uniaxiale perpendiculairement à la surface.
Ces données essentielles pour la simulation numérique X-FEM de la propagation des fissures de
fatigue sous chargement de fretting. Un exemple est présenté dans la section suivante. L’objectif
est de retrouver numériquement le trajet de propagation expérimental des fissures obtenu dans la
section 3.2.
127
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
numérique des essais. L’objectif ici est de comparer les résultats de la simulation numérique de la
propagation des fissures sous chargement de fretting fatigue et ceux des essais correspondants.
Le problème est résolu pour chacune des 40 étapes de chargement avec une précision de
10−4 en considérant l’indicateur d’erreur η (2.34). On considère la loi de contact unilatéral et un
frottement de Coulomb µΓ = 0.3 entre les faces des fissures. Le calcul des facteurs d’intensité
des contraintes en pointe de fissure indique que le mode I est prédominant pour ce problème.
On choisit ici d’utiliser le critère de propagation en fatigue multiaxiale non proportionnelle de
Hourlier 1 (k1∗ (θ,t)max ) (cf. paragraphe [Link]). La simulation de la propagation des fissures est
conduite en plusieurs étapes. Pour réduire les temps de calcul, le pas de propagation da varie en
fonction de la longueur de la fissure. On réalise 12 pas de 5 µm, 10 pas de 8 µm, 7 pas de 18
µm et 15 pas de 20 µm. La figure 3.35 représente les trajets de la fissure 2 (de droite) issus de la
simulation numérique de la propagation des fissures et des mesures expérimentales entre 100 000
et 500 000 cycles.
De 100 000 à 250 000 cycles, la longueur de la fissure calculée numériquement (a250Knum =
242 µm) est supérieure celle mesurée expérimentalement (a250Kexp = 146 µm). En revanche, la
longueur de la fissure calculée numériquement (a500Knum = 647 µm) entre 250 000 et 500 000
cycles est légèrement inférieur à la longueur mesurée (a250Kexp = 700 µm). L’erreur commise par
l’approximation linéaire par morceaux de l’évolution du taux de croissance expérimental (3.9) et
l’utilisation de pas de propagation supérieurs à 1% de la longueur de la fissure sont à l’origine de
cet écart. Le critère (k1∗ (θ,t)max ) utilisé donne une approximation correcte de la direction de pro-
pagation des fissures. Globalement, la simulation de la propagation des fissures est en bon accord
avec les mesures expérimentales. Ces résultats démontrent la validité de la stratégie numérique
proposée dans ce mémoire pour simuler la propagation des fissures sous chargement de fretting.
128
Simulation numérique X-FEM 2D de la propagation de fissures expérimentale sous chargement
de fretting
F IGURE 3.34: Coupes métallographiques des éprouvettes des essais conduits à 100 000,
250 000, 500 000 et 750 000 cycles pour une contrainte de traction σt = 70 MPa. La
fissure 2 (de droite) est représentée.
129
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
130
Simulation numérique X-FEM 2D de la propagation de fissures expérimentale sous chargement
de fretting
Bilan :
131
3. Études expérimentale et numérique bidimensionnelles couplées de la propagation des
fissures sous chargement de fretting-fatigue
132
Chapitre 4
Deux séries d’essais de fretting fatigue en configuration sphère/plan avec différents pa-
ramètres de charge et différents matériaux ont été réalisés dans le cadre de ce projet.
Les géométries tridimensionnelles spécifiques des fissures de fretting à différents stades
de propagation sont reconstituées à partir des coupes métallographiques post-mortem
des éprouvettes. La stratégie de simulation numérique dédiée aux problèmes de fretting
proposée dans le chapitre 3 est étendue ici aux problèmes tridimensionnels. À partir
des données d’essais enregistrées (Q-δ) et en quantifiant l’énergie dissipée à l’interface
dans l’anneau de glissement [FOU 95], l’évolution du frottement local µ est détermi-
née. L’aire de contact, les distributions de pression normale et de cisaillement cyclique
ainsi que la répartition glissement-adhérence sont déterminées. Le risque d’amorçage
des fissures en surface et en volume est analysé à partir des champs des contraintes
et des déformations calculés. Ces conditions de chargement et les faciès de fissuration
expérimentaux 3D sont utilisés comme données d’entrée pour la simulation numérique
X-FEM multi-échelle de ces essais. Les géométries 3D complexes des fissures expéri-
mentales sont décrites précisément par la méthode des fonctions de niveau. L’objectif
est de démontrer les performances du modèle et sa capacité à capturer précisément la
solution du problème de contact entre les faces des fissures en trois dimensions. Les fac-
teurs d’intensité des contraintes sont calculés le long des fronts des fissures au cours du
cycle afin d’étudier numériquement la propagation des fissures de fatigue sous charge-
ment multiaxial non proportionnel.
Sommaire
4.1 Simulation X-FEM 3D d’un essai de fretting fatigue . . . . . . . . . . . 135
133
4.1.1 Essais de fretting fatigue sphère/plan . . . . . . . . . . . . . . . . 135
4.1.2 Analyse des essais : Résolution du contact sphère/plan et prédiction
du risque d’amorçage des fissures . . . . . . . . . . . . . . . . . . 137
4.1.3 Modélisation 3D X-FEM multi-échelle du problème expérimental
de fretting fatigue . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
4.1.4 Calcul des facteurs d’intensité des contraintes . . . . . . . . . . . . 147
4.2 Étude expérimentale et numérique X-FEM de la propagation tridi-
mensionnelle des fissures sous chargement de fretting fatigue . . . . . . 149
4.2.1 Essais de fretting fatigue sphère/plan précontraints : étude expéri-
mentale de la propagation des fissures . . . . . . . . . . . . . . . . 149
4.2.2 Analyse des essais : Résolution du contact sphère/plan et prédiction
du risque d’amorçage des fissures. . . . . . . . . . . . . . . . . . . 153
4.2.3 Modélisation 3D X-FEM multi-échelle du problème expérimental
de fretting fatigue . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
134
Simulation X-FEM 3D d’un essai de fretting fatigue
F IGURE 4.1: Géométrie, chargement et conditons aux limites des essais de fretting
sphère/plan ERC-SKF.
Ces essais a été réalisés sur un banc expérimental dédié à l’analyse des problèmes de fretting
au Centre de Recherche et d’Ingénierie SKF. Une bille de rayon R = 6.4 mm est mise en contact
avec une éprouvette parallélépipédique de dimensions (24mm×16mm×4mm) avec une force P =
120 N (cf. figure 4.1). Un déplacement cyclique sinusoïdal d’amplitude δ = 9 µm et de fréquence
F = 20 Hz est imposé à la bille. Dans le but d’accélérer le phénomène d’amorçage des fissures,
une contrainte de flexion statique σ f = 440 MPa est également imposée à l’éprouvette plane. Le
matériau étudié est un acier de roulement de module de Young E = 210 GPa et de coefficient de
Poisson ν = 0.3. Les surfaces en contact sont polies et nettoyées.
La force tangentielle Q générée par le déplacement relatif des deux corps en contact est me-
surée pendant la durée des essais à l’aide d’une cellule de force liée à l’éprouvette. Cette mesure
permet de calculer le rapport entre les charges tangentielle et normale : f = Q/P. f est égal au
coefficient de frottement local µ en condition de glissement total et inférieur à µ en condition de
glissement partiel (en présence d’une zone d’adhérence). Les paramètres d’essais (P et δ) ont été
sélectionnés afin de se placer en élasticité linéaire et en régime de glissement partiel à l’interface
135
4. Couplage expérimental et numérique pour l’analyse de la propagation des fissures
tridimensionnelles sous chargement de fretting-fatigue
F IGURE 4.2: Boucle d’hystérésis (Q-δ) de fretting obtenue à partir de la mesure expéri-
mentale de l’effort tangent et du déplacement relatif cyclique.
entre les deux corps pendant toute la durée des essais. En effet, la forme allongée des boucles
d’hystérésis (Q-δ) tracées à partir des données expérimentales est représentative d’une condition
de glissement partiel (cf. figure 4.2). L’amorçage des fissures est effectivement l’endommagement
majoritaire. L’usure de la surface de contact est très faible. L’effort tangentiel maximum enregistré
est Qmax = 95 N. Les données de ces essais sont répertoriées dans le tableau 4.1.
E ν P Qmax σf F δ
TABLE 4.1: Propriétés mécaniques de l’acier de roulement étudié et paramètres des essais
de fretting précontraints.
Les traces des essais de fretting conduits à 30 000 cycles, 50 000 cycles et 70 000 cycles sont
présentées dans la figure 4.3. La zone d’adhérence au centre se distingue de la zone de micro-
glissement annulaire. Deux fissures principales s’amorcent en périphérie de la zone de contact
symétriquement par rapport au centre de la trace. Elles se propagent en surface selon une trajectoire
semi-elliptique conformément aux résultats observés lors d’essais similaires [CAD 05, DUB 03].
Les géométries des fissures tridimensionnelles en volume sont extraites par une technique de
reconstitution 3D à partir de coupes métallographiques. Des essais sont conduits pour différents
nombres de cycles. Les éprouvettes post-mortem sont ensuite découpées selon des plans parallèles
au plan médian (X,Z) de coordonnée yi . Puis les géométries des fissures en 3D sont reconstituées
à partir des images métallographiques de ces coupes transversales (cf. figure 4.4). Contrairement
aux essais cylindre/plan pour lesquels la géométrie des fissures est quasi constante le long de la
136
Simulation X-FEM 3D d’un essai de fretting fatigue
F IGURE 4.3: Traces de fretting pour les essais sphère/plan conduits à (a) N = 30 000
cycles, (b) N = 50 000 cycles et (c) N = 70 000 cycles.
Les géométries tridimensionnelles des fissures sont donc extraites précisément à différents
stades de propagation à partir de l’observation de coupes métallographique transversales. Ces géo-
métries seront utilisées dans la section 4.1.3 comme données d’entrée pour le modèle X-FEM via
la technique de représentation par fonctions de niveau. Dans le cadre de la stratégie de prédiction
de durée de vie proposée, la première étape est la détermination du chargement surfacique en trois
dimensions associé au contact sphère/plan à partir des données des essais.
137
4. Couplage expérimental et numérique pour l’analyse de la propagation des fissures
tridimensionnelles sous chargement de fretting-fatigue
F IGURE 4.5: Géométries reconstruites des fissures tridimensionnelles pour les essais
conduits à (a) N = 30 000 cycles, (b) N = 50 000 cycles et (c) N = 70 000 cycles.
138
Simulation X-FEM 3D d’un essai de fretting fatigue
à partir des boucles d’hystérésis (Q-δ) enregistrées au cours de l’essai. Les relations liant P, Q, δ,
c et a et un critère de glissement énergétique développé par Fouvry et al. [FOU 95], basé sur la
quantité Wd , permet de déduire le coefficient de frottement local µ dans la zone de glissement
annulaire à l’aide des équations (4.1) et (4.2) (cf. figure 4.7).
6 (1 + (1 − u)5/3 ) − 65 u(1 + (1 − u)2/3 )
A= · (4.1)
5 u(1 + (1 − u)2/3 )
Q f
u= = (4.2)
µP µ
où f = Q/P et A est le rapport de l”énergie de glissement calculée par intégration de la boucle
de fretting sur l’énergie totale (cf. figure 4.6). Après environ 2 500 cycles, les premières macro-
fissures expérimentales apparaissent en surface et le régime de glissement est stable. Le coefficient
de frottement local µ est alors égal à 0.92 (cf. figure 4.7).
139
4. Couplage expérimental et numérique pour l’analyse de la propagation des fissures
tridimensionnelles sous chargement de fretting-fatigue
F IGURE 4.7: Évolution du rapport f aux cours des cycles des essais de fretting permettant
de calculer le coefficient de frottement local µ, le rapport de glissement c/a et la quantité
u.
F IGURE 4.8: Évolution de l’effort tangentiel Q mesurée au cours d’un cycle de fretting
et discrétisée en 25 pas de temps tn .
140
Simulation X-FEM 3D d’un essai de fretting fatigue
F IGURE 4.9: (a) Distributions des pressions normales p(x,y) et (b) des contraintes tan-
gentielles q(x,y,tn = 1) sur la surface de contact sphère/plan interpolées sur le maillage
surfacique X-FEM de 30 × 30 éléments carrés de 12µm de côté.
est constitué de 46 266 tétraèdres. Il est localement structuré et raffiné au voisinage de la zone
d’intérêt afin de conserver la finesse de la solution de contact sphère/plan (cf. figure 4.11). La
taille caractéristique des éléments finis dans cette zone est 12 µm.
Les fissures extraites de l’essai de fretting conduit à 70 000 cycles sont considérées. Les géo-
métries complexes de ces fissures sont décrites à l’aide de fonctions de niveau dans le cadre du
modèle X-FEM à deux échelles (cf. figure 4.12.b). Le front de fissure est défini par l’intersection
des surfaces iso-0 des fonctions de niveau φ et ψ (cf. figure 1.24). Une seule paire de fonction de
niveau est utilisée pour décrire les deux fissures. Le signe de la fonction de niveau φ est donc né-
141
4. Couplage expérimental et numérique pour l’analyse de la propagation des fissures
tridimensionnelles sous chargement de fretting-fatigue
F IGURE 4.12: (a) Définition des fonctions de niveaux pour la modélisation des fissures
de fretting et (b) dicrétisation de l’interface de ces fissures : 902 paires de points d’inté-
gration.
gatif entre les deux fissures et positif en dehors (cf. figure4.12.a). Le vecteur normal nΓ du repère
local associé à la fissure est orienté automatiquement vers l’extérieur en accord avec le signe de
la fonction de niveau. La fissure 1 est située à la position x = -170 µm dans le plan médian (y=0).
Sa longueur en surface est 691 µm et sa profondeur maximale de 246 µm (≈ 25 éléments finis).
La fissure 2 est en x = +170 µm dans le plan médian (y=0). Sa longueur en surface est 603 µm
et sa profondeur maximale de 96 µm (≈ 10 éléments finis). La discrétisation de l’interface des
fissures est raffinée par subdivision indépendamment du maillage éléments finis (cf. section 2.2).
Une taille critique lc = 7 µm pour les arrêtes des éléments d’interface est définie. Une distribution
de 902 paires de points d’intégration répartie sur les deux fissures représentée dans la figure 4.12.b
est obtenue.
Les propriétés mécaniques sont celles du matériau 100C6 des essais. Le module de Young et
le coefficient de Poisson du matériau valent respectivement E = 210 GPa et ν = 0.3. L’hypothèse
142
Simulation X-FEM 3D d’un essai de fretting fatigue
de comportement élastique linéaire en volume est retenue. La loi de contact unilatérale et la loi de
frottement de Coulomb avec un coefficient de frottement µΓ = 0.5 sont considérées à l’interface.
Le problème est résolu à une précision de 10−4 avec l’indicateur d’erreur locale η (2.34). Le
solveur stabilisé NLLA est utilisé. Les valeurs de la direction de recherche et du coefficient de
stabilisation sont respectivement k = 1+14 Pa.m−1 et ε = 1+12 Pa.m−1 . Les solutions locale et glo-
bale du problème, i.e. le champ de déplacement global U, les efforts de contact T, le glissement
−[WT ] = W− + + −
T − WT et l’ouverture [WN ] = WN − WN entre les faces de la fissure, sont représen-
tées figures 4.13, 4.14 et 4.14 pour 9 pas de temps parmi les 25 considérés.
Note : La partie centrale de l’aire de contact est enfoncée dans la matière sous l’influence du char-
−
gement normal P. Étant donnée le signe de la fonction de niveau φ, les vecteurs [WT ] = W+ T − WT
sont dirigés vers le haut. Par souci de clarté, les vecteurs −[WT ] sont tracés.
Les fissures sont soumises à des séquences de conditions de contact avec frottement complexes
pendant le cycle de fretting. Au temps t1 , l’effort tangentiel est maximal : Q = +95 N. La fissure
2 est ouverte et la fissure 1 est fermée. Les solutions en déplacement global U et local W sont
concordantes. La pression normale de contact P induit un glissement maximal k[WT ]k = 1.4 µm
entre les faces de la fissure 2 et un glissement quasi nul entre les faces de la fissure 1 en contact
avec frottement. Les efforts de contact entre les faces de la fissures 1 sont réguliers et de valeur
maximale kTk = 0.87 GPa.
L’effort tangentiel Q diminue jusqu’au pas de temps t7 pour lequel il est nul. Les deux fissures
sont alors fermées ([WN ] ∼= 0). Les efforts de contact entre leurs faces sont très faibles et le glisse-
ment maximal est égal à 2 µm et 1.6 µm pour les fissures 1 et 2 respectivement. Puis Q augmente
ensuite progressivement jusqu’à atteindre une valeur Q = -95 N au pas de temps t13 . La fissure 1
est alors ouverte : k[WN ]k ∼= 1.15 µm. Les faces de la fissure 2 sont adhérentes : [WT ] = 0. Elles le
restent jusqu’au pas de temps t18 : ∆[WT ] = 0. Les efforts de contact entre les faces de la fissures
2 sont réguliers et atteignent une valeur maximale kTk = 1.47 GPa.
Au pas de temps t19 , l’effort tangentiel Q est de nouveau nul, le glissement entre les faces de
la fissure 1 atteint une valeur maximale stable k[WT ]k = 3.12 µm.
Au pas de temps t25 , le chargement est identique à celui du temps t1 . Cependant, les solu-
tions de contact sont différentes. En effet, la formulation quasi statique du modèle X-FEM à deux
échelles prend en compte les effets d’histoire. Les conditions de contact avec frottement pendant
le cycle complet ont une forte influence sur la solution du problème. De plus, l’hypothèse initiale
pour laquelle les champs locaux d’interface W et T sont nuls au temps t0 peut influencer la so-
lution de contact. Cette influence diminue au cours du cycle. En effet, les effets d’hystérésis se
limitent aux portions de cycles pour lesquels une fissure reste fermée, l’ouverture "réinitialisant"
l’état de contact. sont négligeables à partir du moment où les fissures sont ouvertes. Un deuxième
cycle est alors simulé afin de se placer en régime stabilisé.
Le modèle X-FEM à deux échelles permet de capturer les conditions non linéaires de contact
entre les faces des fissures 3D non planes à une échelle pertinente. De très petits déplacements
relatifs : k[W]k ≤ 3.12 µm, i.e. 200 fois plus petit que la longueur de la fissure (691 µm) et 7700 fois
plus petit que la longueur de la structure (24 mm) sont capturés. Les solutions au cours du cycle
sont stables (pas d’oscillations) et les effets d’hystérésis sont pris en compte. Ceci est un point clé
de la modélisation des problèmes de rupture avec contact et frottement tels que les problèmes de
fretting. En effet, le calcul des facteurs d’intensité des contraintes en 3D présenté dans la section
suivante en dépend (cf. section [Link]).
143
4. Couplage expérimental et numérique pour l’analyse de la propagation des fissures
tridimensionnelles sous chargement de fretting-fatigue
\0
a
C?
N "
r<i
a
\0
0
0;
M
N
'!!
-... 0
144
Simulation X-FEM 3D d’un essai de fretting fatigue
}. }. }. 8
ID
0
,;,
<n
,: ::
E
·::'j:~i:t.~Ji~
~
".\ .
.'
..,., ,,
E
(::
xl!< 8
ID
0
,;,
----- M
N
'!!
E
::;J
145
4. Couplage expérimental et numérique pour l’analyse de la propagation des fissures
tridimensionnelles sous chargement de fretting-fatigue
·.~.
'~~,;:;~
0
jf'
E
"''i'
N "
,.;
Ê
._,
0:::
~
146
Simulation X-FEM 3D d’un essai de fretting fatigue
Bilan :
Dans la première partie de ce chapitre, nous avons montré les performances du modèle X-FEM à
deux échelles. Il permet de capturer précisément les conditions de contact avec frottement entre
les faces de la fissure. Il s’insère dans la stratégie multi-modèles proposée pour la simulation de la
propagation des fissures sous chargement de fretting en trois dimensions. Les facteurs d’intensité
des contraintes sont calculés le long des fronts de fissures non plans à l’aide de l’intégrale d’inter-
action en 3D. Il est donc possible d’étudier numériquement la propagation des fissures en reliant
ces grandeurs caractéristiques de la mécanique de la rupture au taux de croissance des fissures.
Dans ce but, une deuxième série d’essais a été réalisée dans le cadre de ce travail afin d’éva-
luer précisément l’évolution des fissures, i.e. longueur, orientation, et géométrie 3D, au cours des
cycles de charge.
147
4. Couplage expérimental et numérique pour l’analyse de la propagation des fissures
tridimensionnelles sous chargement de fretting-fatigue
Pa.m1ll
Pa.m1 n.
5.66e•006
4.8e+006
_ _ _
-
Km
===Ïï[Link]-
Km
A~A+nnn
0.000171
0.000171
===:::::~24e•n~
?
Pa.m111
[Link].
1.02e•007
1.7e•007
===-- - --
==== =-- - --
8.51e+006
0.000171
0.000171
Kn
Kn
">1A+nM
Pa.m1n.
_ 127e•007 Pa.m n.
1
7.53e+006
====----
==:lili___
377e+006
F: ~"iA+MR
0.000171
0.000171
KI
KI
=
0.000122
0.000122
t13
tl
F IGURE 4.16: Facteurs d’intensité des contraintes KI , KII et KIII calculés aux pas de
temps t1 (Q = -Qmax ) et t13 (Q = Qmax ) le long des fronts des fissures à l’aide de l’intégrale
d’interaction I en 3D (les valeurs positives sont orientées dans le sens des x positifs).
148
Étude expérimentale et numérique X-FEM de la propagation tridimensionnelle des fissures sous
chargement de fretting fatigue
Une expertise des fissures en volume est réalisée par reconstitution des images métallogra-
phiques des coupes transversales des éprouvettes fissurées. Le protocole est similaire à celui pré-
senté dans la section 4.1.1. Les étapes de polissage et de digitalisation des fissures ont toutefois
été automatisées pour augmenter la précision des résultats. La section de l’éprouvette dans le plan
médian du contact pour l’essai à 250 000 cycles est présentée dans la figure 4.18. Les deux fissures
sont parfaitement symétriques en volume. Ce résultat révèle la précision élevé du protocole des
essais et des techniques de découpe et de polissage des échantillons. Les géométries 3D recons-
truites pour les essais conduits à 100 000, 250 000, 500 000 et 750 000 cycles sont représentés
dans la figure 4.19.
Note : la représentation des fissures est inversée (surface en bas) par souci de clarté.
Ces essais de fretting sphère/plan fournissent une représentation précise de l’évolution de la
géométrie des fissures en volume et en surface au cours Ils sont modélisés en trois dimensions en
149
4. Couplage expérimental et numérique pour l’analyse de la propagation des fissures
tridimensionnelles sous chargement de fretting-fatigue
F IGURE 4.17: Traces de fretting pour des essais sphère/plan conduits à : a) 100 000, b)
250 000, c) 500 000 et d) 750 000 cycles.
appliquant la stratégie de simulation numérique proposée dans ce mémoire. Premièrement, les ré-
sultats expérimentaux sont analysés. Le contact sphère/plan est résolu dans la section suivante afin
de déterminer le chargement cyclique en surface. L’objectif est ensuite de modéliser précisément
ces essais afin de caractériser numériquement la propagation des fissures en 3D sous chargement
de fretting.
150
Étude expérimentale et numérique X-FEM de la propagation tridimensionnelle des fissures sous
chargement de fretting fatigue
F IGURE 4.18: Section de l’éprouvette fissurée dans le plan médian du contact pour l’essai
à 250 000 cycles (représentation inversée).
151
4. Couplage expérimental et numérique pour l’analyse de la propagation des fissures
tridimensionnelles sous chargement de fretting-fatigue
152
Étude expérimentale et numérique X-FEM de la propagation tridimensionnelle des fissures sous
chargement de fretting fatigue
F IGURE 4.20: (a) Distributions des pressions normales p(x,y) et (b) des contraintes tan-
gentielles q(x,y,tn = 1) sur la surface de contact sphère/plan. La zone d’intérêt est discré-
tisée en 256×256 cellules carrées de 20µm de côté.
153
4. Couplage expérimental et numérique pour l’analyse de la propagation des fissures
tridimensionnelles sous chargement de fretting-fatigue
Les champs des contraintes et des déformations en volume sont calculés à partir du charge-
ment cyclique en surface. 64 cellules sont considérées en profondeur avec un pas de 5 µm. Le
risque d’amorçage des fissures de fretting en 3D est déterminé en utilisant le critère de Dang Van
(cf. section 3.1.3). Les figures 4.21.a et4.21.b représentent le risque d’amorçage d en volume (sur-
face iso-valeur d = 0.579) et sur la surface de contact respectivement. Deux zones opposées sont
identifiées au voisinage du bord de la zone de contact. d atteint une valeur maximum de 0.724
pour les positions x = 1.86 mm et x = -1.86 mm sur le plan médian (y=0). Ces résultats sont en bon
accord avec les géométries des fissures expérimentales observées à la fois en surface et en volume.
Dans la section suivante, les distributions des contraintes cycliques en surface déterminées ici
et les géométries des fissures issues des résultats expérimentaux présentées dans la section 4.2.1
sont utilisés comme données d’entrée pour la modélisation X-FEM du problème de fissuration
avec contact et frottement.
154
Étude expérimentale et numérique X-FEM de la propagation tridimensionnelle des fissures sous
chargement de fretting fatigue
F IGURE 4.22: KI : facteurs d’intensité des contraintes en mode I calculés au pas de temps
t1 (Q = Qmax ) le long des fronts de fissure à l’aide de l’intégrale d’interaction I en 3D
pour les essais conduits à : a) 100 000, b) 250 000, c) 500 000 et d) 750 000 cycles
(Représentation inversée ; les valeurs positives sont orientées dans le sens des x positifs).
Les valeurs de K varient le long des fronts. Les valeurs maximum sont obtenues au point le plus
bas de ces fronts pour KI et KII . Pour KIII , le comportement est différent. Les valeurs maximum
de KI pour la fissure 2 (des valeurs similaires sont obtenues pour la fissure 1) augmentent avec le
nombre de cycles et la taille de la fissure, puis diminuent. Ainsi (cf. figure 4.25), KI passe de 9.6
√
à 15.2, 22.2 puis 17.1 MPa m. Le comportement est caractéristique d’un chargement de fatigue
tribologique, différent de celui obtenu en fatigue conventionnelle où K augmente avec la taille de
la fissure. En effet, lors de l’amorçage et du début de la propagation, les fissures sont situées dans
155
4. Couplage expérimental et numérique pour l’analyse de la propagation des fissures
tridimensionnelles sous chargement de fretting-fatigue
F IGURE 4.23: KII : facteurs d’intensité des contraintes en mode II calculés au pas de
temps t1 (Q = Qmax ) le long des fronts de fissure à l’aide de l’intégrale d’interaction I en
3D pour les essais conduits à : a) 100 000, b) 250 000, c) 500 000 et d) 750 000 cycles
(Représentation inversée ; les valeurs positives sont orientées dans le sens des x positifs).
la zone en proche surface du contact deux-corps. Les contraintes dans cette zone confinée sont très
élevées et activent les mécanismes moteurs de la propagation des fissures. Lors de l’extension, les
fronts de fissures s’enfoncent en profondeur et s’éloignent de cette zone. Les valeurs de KI (ainsi
que celle de KII ) diminuent et la propagation ralentit. L’arrêt de ces fissures est possible. Cette
transition stade II/stade III (cf. figure 1.8) se situe ici à une profondeur d’environ 1 500 µm.
Bilan :
Des essais de fretting sphère/plan précontraints ont été réalisées au Laboratoire de Tribologie
et Dynamique des Systèmes (LTDS) de l’École Centrale de Lyon avec pour objectif d’étudier la
propagation des fissures tridimensionnelles sous chargement de fretting fatigue. Les géométries
des fissures en 3D ont été reconstituées à partir de coupes métallographiques des éprouvettes d’es-
156
Étude expérimentale et numérique X-FEM de la propagation tridimensionnelle des fissures sous
chargement de fretting fatigue
F IGURE 4.24: KIII : facteurs d’intensité des contraintes en mode III calculés au pas de
temps t1 (Q = Qmax ) le long des fronts de fissure à l’aide de l’intégrale d’interaction I en
3D pour les essais conduits à : a) 100 000, b) 250 000, c) 500 000 et d) 750 000 cycles
(Représentation inversée ; les valeurs positives sont orientées dans le sens des x positifs).
sais conduits à 100 000, 250 000, 500 000 et 750 000 cycles. Dans le cadre de la stratégie de
simulation en trois dimensions de la propagation des fissures en fatigue multiaxiale non propor-
tionnelle, le problème de contact sphère/plan a été résolu incrémentalement pour un cycle complet
par une méthode semi-analytique visant à déterminer les distributions des contraintes normales et
tangentielles cycliques en surface. Le risque d’amorçage des fissures a été déterminé en utilisant
le critère de Dang Van. Ces résultats sont utilisés comme données d’entrée pour la modélisation
X-FEM 3D à deux échelles avec contact et frottement interfacial du problème de fissuration. Les
géométries 3D des fissures sont décrites à l’aide de fonctions de niveau. Les facteurs d’intensité
des contraintes KI , KII et KIII ont été calculés le long des fronts de fissures à l’aide d’une inté-
grale d’interaction en volume (2.67) pour les quatre géométries de fissures correspondant à des
nombres de cycles de charges différents. L’étude des variations des FIC a permis d’identifier la
transition entre les stades II et III de la propagation des fissures en fretting (cf. figure 1.8). Cette
première approche qualitative de la propagation des fissures en trois dimensions permet d’analyser
157
4. Couplage expérimental et numérique pour l’analyse de la propagation des fissures
tridimensionnelles sous chargement de fretting-fatigue
158
Conclusions et perspectives
159
Conclusions et perspectives
Mécanique des Contacts et des Structures (LaMCoS) en collaboration avec l’entreprise SKF. Il
est focalisé sur la prédiction de l’amorçage et du comportement des fissures sous chargement de
fatigue multiaxiale non proportionnelle avec comme application l’analyse de la fissuration des
roulements à billes sollicités en fretting fatigue. Dans le chapitre 1, une synthèse bibliographique
des résultats théoriques et expérimentaux portant sur les problèmes de fretting et notamment sur
la fissuration sous sollicitations de fretting fatigue a été présentée. La maîtrise des mécanismes de
rupture requiert une analyse pluri-disciplinaire des différents phénomènes physiques couplés. Les
besoins mis en évidence pour la simulation numérique de ces problèmes sont :
i la prise en compte du chargement cyclique multiaxial non proportionnel à l’interface entre
les composants de structures en contact.
ii la détermination du risque d’amorçage et de la géométrie initiale des fissure.
iii la prise en compte de plusieurs fissures et de leurs interactions.
iv la prise en compte du contact avec frottement entre les faces des fissures.
v la modélisation de fissures non planes en trois dimensions. En effet, la mixité des sollicita-
tions et les phénomènes non linéaires couplés tels que le contact avec frottement interfacial
génèrent des géométries de fissures complexes en 3D.
vi la simulation de la propagation des fissures en utilisant une loi de fatigue adaptée au ma-
tériau étudié et aux sollicitations, et des critères de direction de propagation en fatigue
multiaxiale non proportionnelle.
Les points iii, v et vi sont aujourd’hui assez bien maitrisés dans le cadre de la X-FEM. Ce-
pendant la modélisation du contact avec frottement entre les faces des fissures pose un certains
nombre de difficultés. Un bref état de l’art de la modélisation des problèmes de contact dans le
cadre de la méthode des éléments finis étendus a été établi. Un certain nombre de travaux ont mis
en évidence deux difficultés majeures liées principalement au type de la formulation du problème
de contact : primale, duale ou mixte.
D’une part, les différents phénomènes physiques mis en jeu font intervenir différentes échelles
spatiales, celle de la structure, celle de la fissure et celle des non-linéarités à l’interface, pouvant
différer les unes des autres de plusieurs ordres de grandeur. Une solution de contact précise re-
quiert une discrétisation spatiale de l’interface adaptée à la dimension caractéristique du problème
local [RIB 07b]. La discrétisation de l’interface des fissures étant généralement dépendante du
maillage éléments finis sous-jacent, une description précise du problème de contact requiert un
raffinement local du maillage en volume. Cette étape peut s’avérer délicate, voire prohibitive pour
les problèmes tridimensionnels. De plus, au regard du problème de la structure, le maillage n’est
pas optimal et entraine un surcoût numérique.
D’autre part, plusieurs travaux ont mis en évidence des problèmes d’instabilité de la solution
de contact liés à un non respect de la condition inf-sup LBB (Ladyzhenskaya-Babuška-Brezzi).
La solution locale de contact (avec ou sans frottement) est perturbée par des oscillations numé-
riques. Ces instabilités ont un impact sur la solution locale du problème (notamment sur les efforts
d’interface) et sur le calcul des facteurs d’intensité des contraintes en pointe de fissure.
À partir de l’étude des différentes formulations utilisées, un modèle X-FEM multi-échelle
avec contact et frottement interfacial a été proposé dans le chapitre 2 afin de répondre aux exi-
gences listées ci-dessus. Il est basé sur une formulation faible mixte à trois champs permettant
une description multi-échelle du problème de contact avec frottement entre les faces de la fissure
160
Conclusions et perspectives
[PIE 10a]. L’interface est considérée comme une structure autonome avec ses propres variables
et sa propre discrétisation indépendante du maillage de la structure et adaptée à l’échelle des non
linéarités locales. Les solveurs non linéaires du Lagrangien Augmenté et de la méthode LATIN
ont été implémentés. De plus, une méthode de stabilisation de la formulation du problème a été
proposée [PIE 10b]. L’étude de la condition LBB et les expérimentations numériques ont montré
que cette méthode est stable. Le solveur stabilisé NLLA permet de s’affranchir des éventuelles os-
cillations numériques de la solution de contact et améliore grandement le taux de convergence de
la méthode itérative. Plusieurs exemples bidimensionnels et tridimensionnels ont mis en évidence
les bonnes propriétés de convergence et la précision élevée du modèle proposé. Le gain en temps
CPU peut atteindre 75%. L’ensemble de ces outils numériques a été mis en oeuvre au sein du code
de calcul ELFE_3D développé au LaMCoS. La comparaison avec le code industriel ANSYS a
permis de valider quantitativement les résultats en 3D.
Des essais de fretting fatigue cylindre/plan réalisés au LTDS de l’École Centrale de Lyon dans
le cadre de ce travail de thèse ont permis de caractériser expérimentalement la propagation des
fissures sous chargement de fretting. Les géométries des fissures ont été expertisées par coupes
métallographiques post-mortem des éprouvettes à différents stades de propagation. Ces résultats
ont permis de tracer une carte de l’évolution des fissures en fonction des paramètres d’essais et du
nombre de cycles de charge. Le modèle X-FEM avec contact et frottement interfacial proposé dans
ce mémoire permet de simuler numériquement les essais de fissuration sous chargement de fret-
ting. Une loi de propagation en vitesse dédiée au problèmes de fretting et au matériau étudié (acier
35NCD16) a été déterminée à partir des résultats de la modélisation X-FEM du problème de fissu-
ration et de l’analyse expérimentale des vitesses de propagation. Puis, la simulation numérique des
essais a été conduite en 2D à l’aide de cette loi et un critère de propagation en fatigue multiaxiale
dédiés aux chargements non proportionnels. Les comparaison des résultats numériques avec les
mesures expérimentales montrent une bonne adéquation. Une quantité plus importante de données
expérimentales permettrait toutefois d’améliorer la précision des résultats. La stratégie proposée
constitue un outil efficace pour la prédiction de la durée de vie des composants de structures sous
sollicitations de fatigue multiaxiale non proportionnelle. Elle permet de répondre aux besoins i, ii,
iii, iv, et vi listés ci dessus.
161
Conclusions et perspectives
162
Conclusions et perspectives
163
Conclusions et perspectives
164
Bibliographie
165
Bibliographie
166
Bibliographie
[BUF 05] B UFFIÈRE J., P ROUDHON H., F ERRIE E., L UDWIG W., M AIRE E., C LOETENS P.
Three dimensional imaging of damage in structural materials using high resolution micro-
tomography. Nuclear Instruments and Methods in Physics Research B, vol. 238, 2005, p. 75-
82.
[BUI 78] B UI H.
Mécanique de la rupture fragile. Masson, Paris, France, 1978, p. 215.
[CAD 05] C ADARIO A., A LFREDSSON B.
Fretting fatigue crack growth for a spherical indenter with constant and cyclic bulk load. Engi-
neering fracture mechanics, vol. 72(11), 2005, p. 1664-1690.
[CAR 00] C ARTER B., WAWRZYNEK P., I NGRAFFEA A.
Automated 3-d crack growth simulation. International Journal for Numerical Methods in En-
gineering, vol. 47, 2000, p. 229-253.
[CAT 38] C ATTANÉO C.
Sul contatto di due corpi elastici : distribuzione locale degli sforzi. Academia dei Lincei, Ren-
dicotti, vol. 27(6), 1938, p. 342-384 ; 434-436 ; 474-478.
[CHA 96] C HAMPANEY L.
Modular analysis of assemblages of three-dimensional structures with unilateral contact condi-
tions. Thèse de Doctorat, École Normale Supérieure de Cachan,France, 1996, p. 138 p.
[CHA 99] C HAMPANEY L.
Modular analysis of assemblages of three-dimensional structures with unilateral contact condi-
tions. Computers and Structures, vol. 73, 1999, p. 249-266.
[CHA 06a] C HAHINE E., L ABORDE P., R ENARD Y.
Crack tip enrichment in the XFEM method using a cut-off function. International Journal for
Numerical Methods in Engineering, vol. 342(7), 2006, p. 527-532.
[CHA 06b] C HAMBON L., J OURNET B.
Modelling of fretting fatigue in a fracture-mechanics framework. Tribology International,
vol. 39, 2006, p. 1220-1226.
[CHA 06c] C HATEAUMINOIS A., BAIETTO -D UBOURG M., G AUTHIER C., S CHIRRER R.
In situ analysis of the fragmentation of polystyrene films within sliding contacts. Tribology
International, vol. 38, 2006, p. 931-942.
[CHE 67] C HEREPANOV G.
The propagation of cracks in a continuous medium. Journal of Applied Mathematics and Me-
chanics, vol. 31(3), 1967, p. 503-512.
[CHO 03] C HOPP D., S UKUMAR N.
Fatigue crack propagation of multiple coplanar cracks with the coupled extended finite ele-
ment/fast marching method. International Journal of Engineering Science, vol. 41(8), 2003,
p. 845-869.
[COL 86] C OLOMBIÉ C.
Usure induite par petits débattements. Thèse Ecole Centrale de Lyon, France, 1986.
[COM 77] C OMNINOU M.
The interface crack with friction in the contact zone. Transactions of the ASME, Journal of
Applied Mechanics, vol. 44, 1977, p. 780-781.
167
Bibliographie
168
Bibliographie
169
Bibliographie
170
Bibliographie
171
Bibliographie
172
Bibliographie
173
Bibliographie
174
Bibliographie
175
Bibliographie
176
Bibliographie
177
Bibliographie
178
Bibliographie
179
FOLIO ADMINISTRATIF
THÈSE SOUTENUE DEVANT L’INSTITUT NATIONAL DES SCIENCES APPLIQUÉES DE LYON
RÉSUMÉ :
La prévision du comportement et de la propagation de fissures de fatigue est un domaine de recherche complexe, dont les ob-
jectifs peuvent se résumer simplement par "augmenter la durée de vie des structures, tout en évitant leur surdimensionnement".
Le fretting est attribué à l’interaction de deux composants en contact soumis à des déplacements relatifs cycliques de faible am-
plitude. Les sollicitations tribologiques à l’interface des composants des structures mettent en jeu des amplitudes de contraintes
très élevées dans des volumes de matière très localisé[Link] régime de glissement locale, i.e. l’évolution des conditions de glisse-
ment et d’adhérence, détermine la nature de l’endommagement initial : fissuration ou usure. Le travail présenté dans ce mémoire
est focalisé sur la prédiction de l’amorçage et du comportement des fissures de fretting fatigue. Ces fissures tridimensionnelles
sont sollicitées en mode mixte selon des cycles complexes et leur propagation en fatigue est conditionnée par les phénomènes
non-linéaires tels que la plasticité localisée en front de fissure et le contact avec frottement entre les lèvres. Ainsi, la maîtrise des
mécanismes de fissuration sous chargement de fretting fatigue requiert une analyse pluri-disciplinaire des différents phénomènes
physiques couplés afin de prendre en compte les effets globaux et locaux, les sollicitations multiaxiales non proportionnelles, les
lois de comportement et les conditions de contact avec frottement interfacial.
Un modèle tridimensionnel éléments finis étendus multi-échelles dédié au contact avec frottement entre les faces de la fissure est
proposé. Une formulation faible mixte à trois champs permet une définition intrinsèque de la fissure avec sa propre discrétisation
indépendante du maillage de la structure. Un solveur stabilisé adapté de la méthode LATIN est implémenté. Les propriétés de
stabilité et les performances du modèle sont illustrées dans plusieurs exemples bidimensionnels et tridimensionnels. Le modèle
est validé par comparaison avec un code éléments finis industriel.
Une stratégie multi-modèle globale basée sur l’analyse expérimentale et la simulation numérique de la propagation des fissures
est développée afin de prédire la durée de vie de composants en fretting fatigue. Des essais de fretting fatigue sont réalisés afin
d’analyser l’amorçage et la propagation de fissures. Les sollicitations tribologiques au cours du cycle sont déterminées par la
résolution du contact deux-corps. La prédiction du risque d’amorçage des fissures est conduite. Ces résultats sont utilisés comme
données d’entrée pour la modélisation X-FEM des essais de fretting bidimensionnels et tridimensionnels. La simulation de la
propagation des fissures est réalisée à l’aide de critères de fatigue multiaxiale non proportionnels et d’une loi de propagation
expérimentale.
MOTS-CLÉS : : Fretting, fatigue, rupture, éléments finis étendus, contact frottement, propagation des fissures, étude expérimen-
tale, contact deux corps.