Thesis 3
Thesis 3
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 - GÉNIE MÉCANIQUE - GÉNIE CIVIL
par
Koffi Espoir KOUMI
*ScSo : Histoire, Géographie, Aménagement, Urbanisme, Archéologie, Science politique, Sociologie, Anthropologie
The present PhD thesis deals with contact problems between heterogeneous ma-
terials. Nowadays heterogeneous materials are extensively used in several industrial
domains (automotive, aeronautics, aerospace, . . . ). Heterogeneous materials involve
porous materials, aluminum alloys, composites materials (woven composites, inter-
locks 3D, interlocks 2D), metallic or ceramics materials containing impurities (poro-
sities/precipitates). In this work, a contact model based on semi-analytical method
is proposed. A conjugate gradient algorithm is used for a fast resolution of contact
equations. The model can account for one or more isotropic/anisotropic inhomoge-
neities. An approach taking inspiration from the Eshelby equivalent inclusion me-
thod is used in the contact solver to account for the effect of inhomogeneities. 2D
and 3D Fast Fourier Transforms (FFT) are used to speed up the computation. A
numerical method is implemented in order to take into account interactions bet-
ween many heterogeneities. The semi-infinite space/ matrix can be either elastic or
visco-elastic. The model developed in the present PhD thesis can solve indentation,
rolling/sliding or fretting contact problems between heterogeneous elastic materials,
homogeneous or heterogeneous visco-elactic materials. In the case of visco-elastic
materials, the model permits to get the solution in terms of contact pressure distri-
bution, subsurface stresses, apparent friction coefficient, both in the transient and
then steady-state regimes. The model has been validated by performing a compari-
son with the results of a finite element model. The CPU time and memory necessary
are greatly reduced in comparison with the classical finite element method. The mo-
del developed is fast, robust and extremely easy to use. An original experimental
approach was proposed in order to measure the displacement fields at interface of two
contacting bodies. A good agreement between experimental results and numerical
simulations is obtained. Finally, the model is applied on some industrial applica-
tions. A coupling between a finite element model and the semi-analytical code allow
to take into account the effects of structure on contact problem.
Introduction 1
1 Etat de l’art 5
1.1 Contexte industriel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1.1 Présentations des turboréacteurs . . . . . . . . . . . . . . . . 7
1.1.2 Objectif de l’étude . . . . . . . . . . . . . . . . . . . . . . . . 9
1.1.3 Impuretés matériaux, inclusions, défauts . . . . . . . . . . . . 9
1.1.4 Contact aube/disque. . . . . . . . . . . . . . . . . . . . . . . . 10
1.1.5 Modes de défaillances du contact aube/disque sous charge-
ment de fretting. . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2 Matériaux composites. . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.2.1 Nature de la matrice . . . . . . . . . . . . . . . . . . . . . . . 21
1.2.2 Forme des renforts . . . . . . . . . . . . . . . . . . . . . . . . 23
1.2.3 Structure des pièces composites. . . . . . . . . . . . . . . . . . 24
1.3 Homogénéisation et méthodes numériques multi-échelles . . . . . . . . 27
1.3.1 Théorie de l’homogénéisation classique . . . . . . . . . . . . . 27
1.3.2 Quelques approches d’homogénéisation . . . . . . . . . . . . . 31
1.3.3 Méthode FE2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
1.3.4 Calcul direct . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
1.4 Modélisation en mécanique des contacts . . . . . . . . . . . . . . . . 38
1.4.1 Approches analytiques . . . . . . . . . . . . . . . . . . . . . . 38
1.4.2 Approches numériques . . . . . . . . . . . . . . . . . . . . . . 39
1.4.3 Cas spécifique du contact aube/disque, aube métallique . . . . 41
1.4.4 Cas spécifique du contact aube/disque pour une aube composite 43
1.5 Synthèse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
ii
iii
iv
Bibliographie 297
vi
Présentation du manuscrit
Le premier chapitre introduit le contexte industriel notamment la problématique
du contact dans les moteurs aéronautiques. Une présentation succincte des tur-
boréacteurs sera faite afin de faciliter la compréhension du contexte industriel de
l’étude. Un zoom sera fait sur le contact aube disque ainsi que sur les différents
modes de défaillances. Une classification des matériaux composites suivant plusieurs
critères sera réalisée. Les différentes méthodes multi-échelles d’estimation du com-
portement ou de modélisation des matériaux multi-échelles seront présentées. Les
différentes approches de résolution des problèmes de contact seront ensuite détaillées.
Un accent particulier sera mis sur la difficulté inhérente à la modélisation du contact
entre matériaux hétérogènes.
Le deuxième chapitre présente les méthodes semi-analytiques appliquées à la
mécanique du contact. Ces approches sont utilisées depuis plus d’une dizaine
d’années maintenant par l’équipe du Professeur Daniel Nélias au LaMCos de l’INSA
de Lyon. Ces méthodes sont basées sur des solutions analytiques et utilisent des
méthodes numériques spécifiques. Ainsi donc, dans ce chapitre seront présentés suc-
cessivement les méthodes de FFT et DC-FFT ainsi que les algorithmes de Gradient-
conjugué. La méthode utilisée pour la discrétisation du massif et des équations du
contact sera également détaillée.
Dans le troisième chapitre, le modèle permettant de prendre en compte des
hétérogénéités isotropes/anisotropes dans une matrice élastique isotrope sera
présenté. Les notions d’inclusion équivalente d’Eshelby et les fonctions potentiels
seront reprises en se basant sur les travaux d’Eshelby [ESH 57, ESH 59], Moschovo-
dis et Mura [MOS 75]. Le modèle sera validé par une comparaison avec la méthode
des éléments finis. Une méthode numérique a été mise en œuvre pour la prise en
compte de l’interaction entre plusieurs hétérogénéités. Des cas d’applications sur
massif revêtu, avec revêtement élastique seront présentés. Beaucoup d’améliorations
numériques ont été apportées afin de rendre le code plus robuste et plus rapide.
Une parallélisation OpenMP/MPI et une nouvelle méthode de décomposition ont
été introduites. Les gains en terme de temps de calcul obtenus peuvent aller jusqu’à
11 sur un PC portable classique comparé au cas initial.
Le quatrième chapitre s’intéressera à la théorie du contact viscoélastique. Une ap-
proche originale basée sur les méthodes semi-analytiques pour résoudre le problème
de contact entre matériaux viscoélastiques, homogènes ou hétérogènes est présentée.
Des applications au cas d’indentation, de roulement et de fretting sont présentées.
Le modèle permet entre autres de quantifier le coefficient de frottement apparent
induit par l’effet de la viscoélasticité dans le cas du roulement, glissement ou du
fretting. Le modèle fournit des solutions en termes de champs de pression mais aussi
de champs de contraintes aussi bien en régime permanent que transitoire. Il s’agit
des résultats tout à fait nouveaux comparés à ceux existants dans la littérature.
Dans le chapitre 5, il s’agira principalement de présenter le dispositif expérimental
développé au LAMCOS pour mesurer des champs de déplacements à l’interface des
corps en contact. Ce dispositif a fait l’objet de très nombreuses modifications durant
cette thèse afin d’arriver à des mesures quantitatives et qualitatives. Le dispositif
ainsi que le protocole expérimental seront présentés. La tomographie sera utilisée afin
de connaitre la microstructure des éprouvettes d’essais en contact. Les techniques
de corrélation d’images 2D permettront de remonter aux champs de déplacements à
l’interface des deux corps en contact. Les résultats obtenus seront comparés à ceux
obtenus au travers d’une simulation numérique ou de calculs analytiques. Le dispo-
sitif expérimental a fait l’objet d’un dépôt de brevet SNECMA/LAMCOS.
Dans le dernier chapitre, le modèle sera appliqué sur trois cas industriels différents.
Une méthodologie multi-échelle de calculs sur pièces sera présentée. Un prototype
de couplage entre un modèle éléments finis et le code semi-analytique sera réalisé
pour résoudre le contact aube/disque. Ensuite l’étude de la nocivité de défauts dans
les billes céramiques sera effectuée et enfin une application dans le cas d’un matériau
composite idéalisé. Les résultats dans le cas de tissages 3D ne pouvant être présentés
pour des raisons de confidentialité.
Enfin, une conclusion générale sera formulée ainsi qu’un ensemble de perspectives.
Etat de l’art
Sommaire
1.1 Contexte industriel . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1.1 Présentations des turboréacteurs . . . . . . . . . . . . . . . . 7
1.1.2 Objectif de l’étude . . . . . . . . . . . . . . . . . . . . . . . . 9
1.1.3 Impuretés matériaux, inclusions, défauts . . . . . . . . . . . . 9
1.1.4 Contact aube/disque. . . . . . . . . . . . . . . . . . . . . . . 10
1.1.5 Modes de défaillances du contact aube/disque sous charge-
ment de fretting. . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2 Matériaux composites. . . . . . . . . . . . . . . . . . . . . . . 21
1.2.1 Nature de la matrice . . . . . . . . . . . . . . . . . . . . . . . 21
1.2.2 Forme des renforts . . . . . . . . . . . . . . . . . . . . . . . . 23
1.2.3 Structure des pièces composites. . . . . . . . . . . . . . . . . 24
1.3 Homogénéisation et méthodes numériques multi-échelles . 27
1.3.1 Théorie de l’homogénéisation classique . . . . . . . . . . . . . 27
1.3.2 Quelques approches d’homogénéisation . . . . . . . . . . . . . 31
1.3.3 Méthode FE2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
que celles des étages d’après. La soufflante comprime le flux d’air qui serait
divisé en un flux primaire et secondaire. Le flux secondaire fournit l’essentiel
de la poussée (80%).
– le compresseur Basse Pression (BP) et Haute Pression (HP). Le flux d’air passe
dans les compresseurs BP et HP où il est accéléré avant d’arriver à la chambre
de combustion.
– la chambre de combustion. C’est l’endroit où est injecté le carburant. L’ajout
de carburant à l’air comprimé permet la combustion du mélange ce qui aug-
mente l’énergie et la température du flux.
– les turbines Haute Pression et Basse Pression. Une partie de l’énergie cinétique
des gaz de combustion est utilisée pour actionner le compresseur Haute Pres-
sion, l’ensemble du compresseur Basse Pression et la soufflante. Les gaz de
combustion sont ensuite éjectés à grande vitesse à l’arrière du moteur et leur
énergie cinétique restante assure la propulsion par réaction.
Au sein du moteur, l’ensemble des pièces en contact peuvent être soumis à des
sollicitations de matage, de fretting, de type roulement (avec ou sans glissement).
Suivant la nature des sollicitations, différents types endommagement peuvent être
rencontrés : rayure, écaillage, usure, oxydation, grippage, fissuration, formation de
débris . . . Dans le moteur, les endommagements liés à des problématiques de fretting,
peuvent être observés sur des pièces aussi diverses que les aubes, disques, paliers,
cannelures, roulements, douilles, tambours, joints d’étanchéité, . . . Les conséquences
peuvent alors aller de la perte de cote, à la perte de fonctionnalité, voire à la rup-
ture des pièces en fatigue. Pour les motoristes, il est essentiel d’acquérir une solide
Figure 1.2 – Exemple de défauts de type carbures et inclusions dans les roulements
hybrides ([AZE 15]).
10
Figure 1.3 – Papillon de fatigue dans le 100cr6 autour d’une inclusion d’alumine
[SAN 93].
dans un futur proche par une nouvelle génération de moteurs plus économiques,
plus performants, plus respectueux de l’environnement : les moteurs LEAP avec des
aubes en matériaux composites.
11
(a) (b)
12
Usure
Il existe plusieurs types d’usure : l’usure liée au frottement, l’usure par érosion
(enlèvement de matière par un fluide chargé de particules en contact avec la surface
du matériau), l’usure par cavitation (induite par des ondes de chocs créées par
l’implosion de bulles de vapeur dans les liquides) . . . Dans notre cas de figure, l’usure
observée sur l’aube est une usure liée au frottement (fig.1.9) . Plusieurs approches
de prédiction de l’usure de frottement existent dans la littérature. Elles peuvent être
13
Sollicitations polycycliques
Vol
Sollicitation
Sollicitations oligocycliques
Décollage Atterrissage
Temps
(a) (b)
Figure 1.7 – Sollicitations mécaniques en pieds d’aubes (a) durant les différentes
phases de vols (b).
14
15
Fissuration
16
17
√
σ πa θ θ 3θ
σxx = √ × cos 1 − sin sin
2πr 2 2 2
√
σ πa θ θ 3θ
σyy = √ × cos 1 + sin sin
2πr 2 2 2
√
σ πa θ θ 3θ
σxy = √ × cos sin cos
2πr 2 2 2
(1.2)
Les contraintes tendent donc vers l’infini quand (r → 0). En réalité les
contraintes en pointe de fissure sont finies du fait de la plasticité. Mais en
élasticité pure, les contraintes tendent vers l’infini. Il est alors difficile de tra-
vailler avec ces champs de contraintes. Les équations précédentes peuvent se
mettre sous la forme :
K
σi,j = √ × fi,j (θ) (1.3)
2πr
où K est une variable finie représentant le facteur d’intensité de contraintes
du mode considéré et s’exprime en MPa.m1/2 .
Dans le cas général K s’exprime sous la forme :
√
K = β.σ. πa (1.4)
où β représente un facteur de correction géométrique.
L’amorçage et la propagation des fissures dans le cas du contact aube/disque ont
été et continuent d’être l’objet de nombreuses études [MER 11, DIC 06c, SUN 12].
18
19
Figure 1.13 – Mécanismes de dégradation observés sur un tissé 3D CMO [ELH 06].
20
21
que les thermodurcissables (polyester, époxy). Par contre, compte tenu du fait qu’ils
sont plus ductiles, les thermoplastiques résistent mieux à la fissuration. Ils existent
également des matrices thermostables pour des applications hautes températures.
(a) Volets froids en composites thermostructu- (b) Prototype de mélangeur en CMC sur moteur
raux C/SiC (Avion de combat Rafale). CFM56-5C.
22
23
Figure 1.17 – renforts en forme de (a) particules (b) fibres courtes (c) fibres longues.
24
Taffetas
Satin
Tissus
Sergé
Tressé circulaire
Hybride Adanur et Liao (1998)
Circulaire
Préformes Textiles
Tressé Figuré
Biaxial
Plat Triaxial
Stratifié cousu (non tissé)
Drapier et Wisnom (1999)
Couche à couche
Tissé Interlock 3D Tissage angle Interlock 3D
Interlock orthogonal Adanur et Liao (1998)
Solide 2-pas
Tressé Cartésien 4-pas
Multi-pas
Tubulaire
Tressé 3D
Chaîne de tricotage multiaxial Dexter (1998)
Tricoté Sandwich tricoté
Tricoté 3D multiaxial
Adanur et Liao (1998)
25
matériau composite considéré. Entre l’échelle micro et macro une échelle méso peut
être rajoutée.
26
x1
Figure 1.21 – Image virtuel d’un VER du composite tissé tridimensionnel CMO
orienté dans le repère du contact.
27
[Link] Représentation
Avant de construire le VER, il faudrait vérifier l’hypothèse de séparabilité des
échelles.
Représentation du VER
28
l
Figure 1.22 – Exemple de microstructure périodique.
[Link] Localisation
29
Les moyennes des champs microscopiques doivent être égales aux champs macrosco-
piques.
< ε >= E, et < σ >= Σ (1.6)
L’égalité entre la moyenne des champs locaux et les champs macroscopiques imposés
est vérifiée pour certains types de conditions limites données.
Le tenseur de localisation
Cette étape consiste à établir des relations explicites (ou implicites) entre
les champs microscopiques et macroscopiques. Les tenseurs de localisation des
déformations A et des contraintes B peuvent se mettre sous la forme fonctionnelle
générique suivante :
30
Or, les contraintes locales sont déterminées par la relation de localisation (eq.1.11).
D’où :
E =< s(x) : B(x) : σ 0 >=< s(x) : B(x) >: σ 0 (1.15)
L’expression des modules de souplesses homogénéisés S Hom se met sous la forme :
31
Σ Comportement global E
Homogénéisation
Homogénéisation
Localisation
Localisation
Comportement local
= ( ,…)
= ( ,…)
32
uniquement de déterminer les champs mécaniques moyens par phase. Les champs
locaux microscopiques sont donc moyennés. Il n’y alors aucune information sur la
variation des champs locaux, ni sur leurs répartitions spatiales.
33
Toutes les approches présentées dans cette partie du manuscrit sont difficilement
transposables au cas d’un contact entre matériaux hétérogènes. L’une des limitations
majeures étant l’hypothèse de séparabilité des échelles.
u(x) = ε0 .x + v (1.17)
divσ = 0
σ = c : ε, surV
σ et ε periodiques
< σ >= Σ et < ε >= E (1.18)
34
35
Ω3
Γ23
Ω2
Ω Γ12
Ω1
Il s’agit d’une approche très originale. Cette méthode a été développée par Ben
Dhia [BEN 98, BEN 99] pour la résolution numérique de problèmes mécaniques pour
lesquels une représentation macro (grossière) est suffisante, sauf en des zones loca-
lisées (zones d’intérêts). Dans ces zones, des représentations micro sont nécessaires.
La méthode peut aussi être utilisée pour coupler des modèles différents comme par
exemple des modèles discrets/continus. Cette méthode est fondée sur trois consti-
tuants clés : la superposition de modèles, la partition des énergies et les couplages
énergétiques des modèles superposés.
Considérons le cas d’un solide élastique, occupant un domaine Ω0 (Fig.1.28).
Ce domaine contient par exemple une fissure qui occupe un sous-domaine Ω1 , dans
lequel la modélisation héritée de celle mise en place dans Ω0 , n’est pas assez précise.
Le sous-domaine est donc notre zone d’intérêt micro dans lequel on a besoin d’une
modélisation très fine. Le concept de la méthode Arlequin est de concevoir le modèle
local complet souhaité et de le superposer au modèle global dans la zone d’intérêt.
Dans la zone d’intérêt, il existe deux représentations mécaniques différentes donc
deux états mécaniques différents. Une partition énergétique est ensuite réalisée dans
la zone d’intérêt par l’introduction de fonctions poids qui forment une partition
de l’unité. Le dernier constituant clé de la méthode est le raccord entre les deux
domaines. La liaison entre les deux domaines est faite via l’introduction d’une zone
de recouvrement Ω1g (Fig.1.29). La méthode Arlequin peut s’appliquer à plusieurs
problèmes industriels : pièces tournantes de turboréacteur en structures composites
36
(a) (b)
stratifiées [TOU 12] ; au couplage stochastique-déterministe [PRU 08, PRU 09] ; aux
problèmes d’impact [BEN 04]. Elle a également été utilisée pour la modélisation du
fretting wear [BEN 11] dans le cas d’un contact entre matériaux homogènes. On
peut très bien imaginer une extension de l’approche dans le cas d’un contact entre
matériaux hétérogènes. La difficulté de la méthode Arlequin réside dans le choix
des multiplicateurs de Lagrange recollant les domaines dans leur zone de transition.
Cette difficulté est principalement numérique. Pour certains types de problème cette
méthode est beaucoup plus lourde à mettre en place, avec des temps de calculs plus
élevés que le calcul direct. C’est le cas du contact hétérogène avec une, deux ou une
dizaine d’hétérogénéités par exemple.
Ces dernières solutions même si elles peuvent fournir des résultats assez précis,
sont difficilement envisageables dans un contexte industriel du fait du temps de calcul
et du coût humain que cela peut engendrer. Au problème hétérogène se superpose
la problématique assez complexe de la modélisation du contact.
37
Ω0
Ω1
Ω1f
Ω1g
Ω0
Figure 1.29 – Représentation des modèles macroscopique Ω0 et microscopique Ω1
avec une zone de recouvrement Ω1g .
38
39
finis classiques. Il faut également rappeler que les temps de calculs sont assez élevés,
et le sont d’autant plus que l’existence de forts gradients de contraintes oblige à
utiliser un maillage assez fin. Il s’agit cependant de la méthode la plus aboutie en
mécanique numérique.
40
le cas des contacts bidimensionnels. Les approches de Kalker d’origine utilisent des
algorithmes de type Newton-Raphson. La méthode de Newton Raphson sera rem-
placée plus tard par celle de Gauss-siedel dans l’article de Jaeger [JAE 04]. Plusieurs
techniques d’accélération de temps de calcul seront ensuite introduites. C’est le cas
par exemple des méthodes multi-grilles [BRA 90, LUB 91], les transformées de Fou-
rier rapides [JU 96a, NOG 97, POL 00, LIU 01]. La finesse des discrétisations ren-
due possible par ces méthodes les rendent quasiment incontournables [AI 99] dans
l’étude des contacts rugueux. L’équipe du Professeur Daniel Nélias développe depuis
quelques années un code de contact robuste basé sur les méthodes semi-analytiques.
Les premiers travaux remontent à ceux de Christophe Jacq [JAC 02a, JAC 02b] qui
utilise les approches semi-analytiques pour résoudre un contact élastoplastique 3D.
Le code s’est développé depuis et a été appliqué avec succès à des problématiques :
– de contacts thermo-élasto-plastiques [BOU 05, BOU 07]
– de plasticité et d’accumulation de déformation plastique en mécanique du
contact [BOU 07]
– de roulement élasto-plastique [NEL 07b, CHA 11c]
– d’usure de contacts mécaniques [GAL 06, GAL 07b, GAL 10a, GAL 10b]
– de contact entre matériaux anisotropes [BAG 13, BAG 12b, GAO 15]
– de contact entre matériaux hétérogènes [LER 10, LER 11, KOU 14b]
– du grenaillage [CHA 12a, CHA 11c, CHA 11d]
– de prise en compte des effets dynamiques dans la résolution du contact
– de prise en compte des contacts viscoélastiques hétérogènes, du roulement
viscoélastique sur massif homogène/hétérogène [KOU 14a, KOU 15]
– de fretting contact en présence de massif viscoélastique
Les méthodes semi-analytiques, comparées aux autres méthodes numériques,
présentent l’avantage d’être très rapides. Elles peuvent prendre en compte des
géométries complexes 3D, des phénomènes physiques complexes comme la plasti-
cité, la viscoélasticité . . . Les approches développées dans cette thèse se baseront sur
ces méthodes semi-analytiques.
41
42
Tous les travaux mentionnés ci-dessus se sont limités au cas des matériaux ho-
mogènes isotropes. A notre connaissance, aucune étude complète n’existe de nos
jours sur le contact aube/disque dans le cas des matériaux composites tissés 3D.
43
Figure 1.32 – Modélisation du contact aube-disque par zoom strucural [SIN 02]
1.5 Synthèse
Ce premier chapitre a permis de mettre en exergue la problématique du contact
entre matériaux hétérogènes dans le cas du contact aube composite/disque, aube
revêtu/disque ou des roulements hybrides, . . . Il s’agit d’une problématique très peu
étudiée dans la littérature et qui nécessite de revoir les outils de modélisation et
de résolution des problèmes de contact. Les différentes catégories de composites ont
été ensuite présentées et un accent particulier a été mis sur la spécificité des com-
posites tissés 3D CMO et CMC utilisés par SNECMA. Les différentes méthodes
de modélisation multi-échelles, utilisées dans la littérature pour la modélisation
des matériaux hétérogènes, ont été ensuite détaillées. Il a été montré que les
méthodes d’homogénéisation ne peuvent s’appliquer au contact hétérogène du fait
du caractère très local d’un problème de contact. Les méthodes semi-analytiques
44
ainsi que les autres approches de résolution des problèmes de contact ont été
présentées. A la complexité des problèmes de contact s’ajoute la complexité inhérente
à la modélisation des matériaux hétérogènes. Finalement, ce sont les méthodes
semi-analytiques qui seront retenues et utilisées dans ce manuscrit. Ces méthodes
présentent beaucoup d’avantages. Outre le fait que le temps de calcul et l’espace
mémoire sont considérablement réduits par rapport à la méthode des éléments fi-
nis, ces méthodes peuvent également prendre en compte simultanément plusieurs
géométries et phénomènes complexes comme la plasticité, la viscoélasticité, . . .
Le chapitre prochain présentera les bases de la méthode semi-analytique afin de
faciliter la compréhension du reste du manuscrit.
45
46
Méthodes semi-analytiques en
mécanique du contact
Sommaire
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.2 Mise en équation du problème . . . . . . . . . . . . . . . . . 50
2.2.1 Cinématique du contact . . . . . . . . . . . . . . . . . . . . . 50
2.3 Discrétisation numérique . . . . . . . . . . . . . . . . . . . . . 52
2.3.1 Equations du problème de contact . . . . . . . . . . . . . . . 53
2.3.2 Contact normal . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.3.3 Contact tangentiel . . . . . . . . . . . . . . . . . . . . . . . . 54
2.3.4 Couplage du problème normal et tangentiel . . . . . . . . . . 55
47
48
2.1 Introduction
Les méthodes semi-analytiques sont de plus en plus utilisées dans la résolution
des problèmes de contact complexes. En plus d’être simples d’utilisation en compa-
raison aux méthodes éléments finis, elles permettent d’avoir des gains considérables
en termes de temps de calcul. L’équipe de Nélias développe depuis quelques années
des modèles de contact basés sur ces méthodes. Les premiers travaux remontent
à ceux de Christophe Jacq [JAC 01, JAC 02a] qui utilise pour la première fois
les méthodes semi-analytiques pour la résolution d’un contact élasto-plastique en
présence d’indents et de défauts de surface. Les travaux de Jacq [JAC 01, JAC 02a]
sont basés sur les solutions intégrales de Chiu [CHI 78a, CHI 77a]. Les solutions
analytiques élémentaires développées par Jacq [JAC 01, JAC 02a] permettent
ainsi de déterminer les contraintes résiduelles induites par un cuboı̈de élémentaire
de déformation plastique noyé dans un massif semi-infini. Les solutions analy-
tiques ont été intégrées à un algorithme de résolution d’un contact tridimension-
nel afin de prendre en compte l’effet des déformations plastiques induites en sous-
couche. Les développements de Jacq [JAC 01, JAC 02a] ouvrirent la porte à d’autres
problématiques de contacts élastoplastiques et de contact entre surfaces fractales. Le
code s’est développé depuis et a permis de prendre en compte des phénomènes phy-
siques assez complexes. En se basant sur les travaux de Jacq, Boucly [BOU 05] prend
en compte les aspects d’échauffement thermique dans la résolution du problème de
contact élastoplastique. En parallèle, le modèle de plasticité est enrichi avec un al-
gorithme de retour radial initialement développé par Simo et Taylor [SIM 85]. La fi-
nesse de discrétisation rendue possible par cette méthode fait qu’elle a été également
utilisée pour l’étude des contacts rugueux. Antaluca et al [ANT 05, ANT 08] étudia
l’effet sur, des surfaces dentées, d’un chargement tangentiel (hypothèse de glisse-
ment total). Il étudia également le cas d’une charge roulante sur un massif ayant
un comportement élasto-plastique pour des applications sur roulement à rouleaux
cylindriques [NEL 07a, CHA 11c]. Le cas d’un contact roulant/glissant en présence
d’aspérités en surface a été également étudié par [BOU 07]. En se basant sur l’ap-
proche simplifiée de Antaluca en glissement total [ANT 05], Gallego [GAL 07a] pro-
pose un modèle élastique de contact sous condition de stick-slip. Le modèle ainsi
développé a été appliqué avec succès à des problématiques de fretting wear en glis-
sement total [GAL 06], glissement partiel [GAL 07c]. Il est dès lors possible de si-
muler les 3 modes de fretting [GAL 10a, GAL 10c]. Dans le modèle de Gallego, les
problèmes de contact normal et tangentiel sont résolus alternativement. Une exten-
sion du modèle de Gallego au cas de massifs élastiques revêtus a été introduite par
[WAN 10]. Chaise [CHA 11a, CHA 12a, CHA 11b, CHA 12b] a étendu l’usage des
modèles semi-analytiques aux problématiques d’impacts. Son modèle a été utilisé
pour étudier la mise en compression des surfaces par des procédés de grenaillage.
Le modèle a été amélioré afin de prendre en compte les effets d’inertie dans les
problématiques d’impact. Fulleringer [FUL 11] durant sa thèse a déterminé des so-
lutions analytiques permettant de prendre en compte l’effet de la plasticité dans la
49
Suite aux diverses thèses, il est de nos jours très difficile de faire ressortir un
algorithme unique pour le code semi-analytique. Plusieurs phénomènes sont pris en
compte de diverses manières. Pour faciliter la compréhension du reste du manus-
crit, nous allons présenter la méthode juste dans le cas d’un problème de contact
normal/tangentiel élastique avec calculs des contraintes en sous couche. Ce cas per-
mettra de se familiariser avec la méthode tant d’un point de vue mise en équation
que modélisation numérique. Cette étape est très importante et facilitera amplement
la compréhension du document. Plusieurs algorithmes numériques avancés seront
présentés et seront intensivement utilisés dans les chapitres 3 et 4.
(1) (2)
x3 = f (1) (x1 , x2 ), x3 = f (2) (x1 , x2 ), (2.1)
50
W
M2
O
hi x1
x2
x3
M1
On supposera que les pentes de f (1) et f (2) sont assez faibles pour approcher ces
surfaces par le plan X1 − X2 . Cette hypothèse est très importante pour la suite des
travaux présentés dans cette thèse.
La force normale W est transmise au contact par l’intermédiaire des contraintes
surfaciques (p (x1 , x2 )) dans la zone de contact ΓC . La distance entre les deux corps
en contact h(x1 , x2 ), qui permet de définir la condition de contact, est construite à
partir de la séparation initiale des corps hi (x1 , x2 ), le déplacement de corps rigide
(1) (2)
δ = δ (1) + δ (2) et le déplacement normal élastique des deux surfaces u3 = u3 + u3 ,
Le contact tangentiel fait intervenir, des glissements à l’interface des deux corps en
contact. Le terme glissement sera défini par le vecteur s à deux composantes sx1
et sx2 . st sera le vecteur glissement à l’instant t. L’opérateur ˙ = dtd sera l’opérateur
dérivé par rapport au temps. Soit xt(1) et xt(2) , les coordonnées
! de deux points ! en
t(1) t(2)
t(1) u1 t(2) u1
contact des corps 1 et 2 à l’instant t. Soit uτ = t(1) et uτ = t(2) les
u2 u2
déplacements élastiques tangentiels des surfaces en contact. Quand le système passe
d’un instant t à t0 , la variation de chaque terme sera supposée linéaire. Les points
51
0 0
en contact subissent un déplacement de corps rigide et se retrouve en xt (1) et xt (2) .
Le glissement s’écrit
n 0 0
o
0
t t
ut(1) ut(2) utτ (1) utτ (2) − ẋt(1) − ẋt(2) (t − t0 ) (2.6)
s = ṡ (t − t ) = τ − τ − −
˙ − xt(2)
Les vecteurs st et ṡt sont colinéaires. Le terme xt(1) ˙ est la contribution
des déplacements de corps rigide et est fonction de la variation de δx1 , δx2 . Il ressort
que :
˙ − xt(2)
xt(1) ˙ 0
(t − t0 ) = xt(1) − xt (1) − xt(2) − xt (2) ,
0
(2.7)
t(1) 0 0
− δx1 t (1) + δx1 t(2) − δx1 t (2)
˙ ˙
0 δx1
x −x
t(1) t(2) (t − t ) = 0 0 . (2.8)
δx2 t(1) − δx2 t (1) + δx2 t(2) − δx2 t (2)
Pour des raisons de simplicité, ∆ sera utilisée pour indiquer la différence d’une
grandeur entre les instants t et t0 . Ce qui revient à écrire :
t(1) t(2)
˙ ˙
0 ∆δx1 + ∆δ x1
xt(1) − xt(2) (t − t ) = . (2.9)
∆δx2 t(1) + ∆δx2 t(2)
∆ut1 − ∆δx1 t
t
s = . (2.10)
∆ut2 − ∆δx2 t
52
x1
x2
Figure 2.2 – Discrétisation de la surface de contact.
où u3 (x1 , x2 ) représente le déplacement élastique total des deux corps au point
(1) (2)
(x1 , x2 ), i.e. u3 (x1 , x2 ) = u3 (x1 , x2 ) + u3 (x1 , x2 ). Les variables p (x1 , x2 ) et
u3 (x1 , x2 ) sont inconnues. L’équation comprend donc deux inconnues. Il est donc
nécessaire de trouver une relation entre les champs de pression et le déplacement
normal u3 (x1 , x2 ). Dans le cas d’un contact entre matériaux élastiques, le
déplacement normal et les champs de pression sont reliés entre eux par la relation
suivante :
Sous forme continue
(1 − ν 2 ) +∞ +∞
Z Z
p (ξ, η) dξdη
u3 (x1 , x2 ) = q (2.16)
πE −∞ −∞ (ξ − x1 )2 + (η − x2 )2
N1 X
X N2
u3 (i, j) = up3 (i, j) = K3p (i − k, j − l) p (k, l) (2.17)
k=1 l=1
53
Ces résultats sont basés sur les solutions de Love [LOV 20] pour un massif élastique.
Les expressions analytiques de p (x1 , x2 ) et de u3 (x1 , x2 ) sont à priori inconnues.
Et pour des géométries de contact complexes, il est quasiment impossible de faire
une résolution analytique du problème de contact. C’est ici qu’intervient l’approche
numérique (de la partie semi-analytique), on utilisera la discrétisation présentée
à la figure 2.2 ainsi que la méthodologie du gradient conjugué, la méthode FFT
ainsi que d’autres méthodes numériques avancées pour résoudre efficacement les
équations du problème. Il faut toutefois rappeler que dans la plupart des travaux
de thèses présentés au paragraphe 2.1, les solutions analytiques élémentaires ne
sont pas connues à priori. Il faudrait d’abord déterminer les solutions analytiques
élémentaires avant de dérouler la méthode semi-analytique (cas dynamique, . . . ).
54
Il arrive assez souvent que le couplage normal/tangentiel ne soit pas pris en compte
pour des raisons de simplification. A chaque itération de calcul on résolvera de
manière séquentielle une fois le problème normal puis une fois le problème tangentiel.
55
56
convolution linéaire y(t) entre deux fonctions x(t) et h(t) est définie par :
Z +∞
y(t) = x(t)h(t − τ )dτ ≡ x(t) ∗ h(t) (2.26)
−∞
Il est possible de faire une analogie entre l’équation ci-dessus (eq.2.26) et l’équation
donnant le déplacement en fonction du champ de pression, équation qui sera rappelée
ci-dessous :
(1 − ν 2 ) +∞ +∞
Z Z
p (ξ, η) dξdη
u3 (x1 , x2 ) = q (2.27)
πE −∞ −∞ 2 2
(ξ − x1 ) + (η − x2 )
Dans les problèmes de contact, la variable y(t) est inconnue, x(t) et h(t) sont sup-
posées connues. Soit eh(ω) et x
e(ω) les transformées de Fourier respectives des fonc-
tions x(t) et h(t). La transformée de Fourier du produit de convolution de l’équation
2.26 conduit à :
ye(ω) = x
e(ω)e
h(ω) (2.28)
Il s’agit dès lors d’une convolution continue. Il est donc possible d’avoir la variable
y(t) en utilisant une transformée de Fourier inverse.
Z +∞
y(t) = 1/2π ye(ω)eiωt dt (2.29)
−∞
57
(
= 0 if x < 0,
H(x) (2.34)
= 1 if x ≥ 0
La convolution discrète se fait sur une longueur infinie à partir des valeurs de
longueurs finies. Ici x et h ont le même nombre de termes j ∈ [0, N − 1]. La fonction
Heaviside H est active lorsque j − r < 0 ce qui permet d’éviter un indice négatif
pour h et se traduit par le remplacement de celui-ci par j − r + N . On a donc une
sommation circulaire, et ainsi introduit une périodicité de l’échantillon (périodicité
de longueur L0 ).
[Link] Source de l’erreur des problèmes de contact résolus avec les FFT
L’utilisation de la convolution cyclique induit des erreurs compte tenu de la non-
périodicité spatiale des variables de contact. Cette convolution induit un phénomène
de recouvrement. Ce phénomène est illustré à la figue 2.3 sur une convolution discrète
unidimensionnelle d’un échantillon de pression p et de coefficients d’influence K. La
convolution cyclique suppose une périodicité de la pression. Au début et en fin du
processus de convolution, les pressions ajoutées à gauche et à droite interfèrent avec
les coefficients d’influence. Le résultat sera donc biaisé à cause de ce phénomène
de recouvrement. Liu et al. [LIU 00a] ont trouvé une solution à ce problème de
recouvrement en introduisant une technique de zero-padding (Fig.2.4) sur la source
que sont les champs de pression puis une technique de wrap-around (Figs.2.5 et 2.6)
sur les coefficients d’influence que sont les coefficients de Green. Cette méthode est
58
connue sous le nom DC-FFT (Discrete Convolution and Fast Fourier Transform) et
sera présentée ci-dessous.
59
K
j
K−j
wrap−around et zero padding
KN=0
0
0 N 2N−1
Kj
K−j
wrap−around et zero padding
KN=0
0 N 2N−1
60
2.4.3 Cas 2D
La transformée de Fourier bidimensionnelle (2D-FFT) consiste à utiliser une
transformée de Fourier suivant une première direction, suivie d’une deuxième suivant
l’autre direction. L’algorithme est utilisé dans le cas d’une source surfacique. L’image
peut quant à elle être :
61
– surfacique
N1 X
X N2
u3 (i, j) = up3 (i, j) = Kcp (i − k, j − l) p (k, l) (2.35)
k=1 l=1
– volumique
N1 X
X N2
p
σIJ (i, j, k) = p(l, m)(CIJ (i − l, j − m, k, ν2 )) (2.36)
l=1 m=1
2.4.4 Cas 3D
La transformée de Fourier tridimensionnelle consiste à utiliser la transformée
de Fourier rapide suivant les trois directions x1 , x2 , x3 . Cet algorithme est prin-
cipalement utilisé dans le cas d’une source volumique et d’une image volumique
N
X 3 −1 N
X2 −1 N
X 1 −1
σij (x1 , x2 , x3 ) = Bijkl (x1 − xI1 , x2 − xI2 , x3 − xI3 )ε∗kl (xI1 , xI2 , xI3 )
xI3 =0 xI2 =0 xI1 =0
(2.37)
C’est le plus rapide dans le cas des sources volumiques. Le domaine d’étude
N1 ×N2 ×N3 est d’abord multiplié par 2 dans les trois directions afin d’obtenir
un produit de convolution cyclique. La transformée de Fourier 3D requiert
alors O(2N1 × 2N2 × 2N3 × log(2N1 × 2N2 × 2N3 )) opérations. La transformée
inverse requiert également le même nombre. Compte tenu du fait que la 3D-
DC-FFT requiert 3 transformées directes, le nombre d’opérations pour la triple
sommation est de l’ordre de O(2N1 × 2N2 × 2N3 + 3 × 2N1 × 2N2 × 2N3 ×
log(2N1 × 2N2 × 2N3 )).
62
Ax = b (2.39)
end while
63
64
Incrément de déplacement
Contraintes résiduelles résiduel dures
(3D-FFT)
(2D-FFT)
oui
Fin
2.7 Synthèse
Ce chapitre a passé en revue les bases de la méthode semi-analytique, méthode
qui sera utilisée dans cette thèse pour la modélisation du contact. L’ex-
posé s’est limité au cas du contact élastique normal/tangentiel afin de fa-
ciliter sa compréhension. De manière générale ; qu’il s’agisse d’un problème
d’impact, d’un contact élastique, d’un contact élastoplastique, . . . ; la mise
65
66
Sommaire
3.1 Etat de l’art . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.2 Formulation générale du problème de contact . . . . 71
3.3 Prise en compte des hétérogénéités . . . . . . . . . . . 71
3.3.1 Solution en espace infini . . . . . . . . . . . . . . . . . . 71
3.3.2 Solution en espace semi-infini . . . . . . . . . . . . . . . 74
3.3.3 Déplacement normal dû aux hétérogénéités . . . . . . . 76
3.3.4 Prise en compte de l’inclinaison de l’inclusion . . . . . . 77
3.3.5 Prise en compte de l’effet des hétérogénéités dans l’al-
gorithme de contact . . . . . . . . . . . . . . . . . . . . 78
3.4 Validation du modèle développé . . . . . . . . . . . . 81
67
68
69
a été résolue par Nélias et co-auteurs [JAC 02a, FUL 11, LER 10, LER 11,
LER 13a], puis ensuite par [ZHO 1a]. Ces modèles sont basés sur les méthodes
semi-analytiques développées par [JAC 02b].
Dans ce chapitre, le cas des hétérogénéités isotropes inclinées et anisotropes
a été résolu. Le modèle a été validé par la méthode des éléments finis. Une
méthode numérique a été mise en œuvre pour la prise en compte de l’inter-
action entre plusieurs hétérogénéités. La méthode a été validée dans le cas
d’un revêtement élastique. Un accent particulier a été mis sur la réduction
du temps de calcul. Une nouvelle méthode de décomposition ainsi qu’une pa-
rallélisation du code de calcul ont été mises en œuvre. En combinant la nouvelle
méthode de décomposition et la parallélisation, il est possible d’avoir un gain
en terme de temps de calculs d’un facteur 11 sur (4 cœurs). Tous les ingrédients
sont désormais réunis pour résoudre le problème de contact entre matériaux
hétérogènes ou encore entre matériaux composites. Les temps de calculs sont
très bas comparés à des modèles éléments finis classiques.
70
h(x1 , x2 ) > 0
contact : hi (x1 , x2 ) = 0 et p(x1 , x2 ) > 0
séparation : hi (x1 , x2 ) > 0 et p(x1 , x2 ) = 0 (3.3)
Toute la difficulté de la résolution du contact entre matériaux hétérogènes
réside dans la détermination du champ de déplacement u∗3 et du champ de
contraintes σ induit par la présence des hétérogénéités. Dans la suite de ce
chapitre, nous tâcherons de détailler l’approche utilisée pour calculer σ et u∗3
ainsi que l’algorithme permettant de coupler le problème de contact et la prise
en compte des hétérogénéités.
71
σij,j = 0 (3.4)
et σij = 0 à l’infini.
Si on se place dans le cadre de l’élasticité linéaire, les composantes du tenseur
des contraintes sont données par la loi de Hooke :
M
σij0 + σij = Cijkl M
(u0k,l + uk,l ) = Cijkl (ε0kl + εkl ) dans D − Ω (3.5)
La méthode de l’inclusion équivalente d’Eshelby consiste à remplacer
l’hétérogénéité par une inclusion ayant les même propriétés matériaux que
la matrice mais soumise à une déformation imaginaire appelée déformation
libre ou eigenstrain ε∗ donnée par la relation :
I
Cijkl M
(ε0kl + εkl ) = Cijkl (ε0kl + εkl − ε∗kl ) dans Ω (3.6)
L’équation d’équivalence (Eq.3.6) est une équation nécessaire et suffisante au
problème d’équivalence entre inclusion et hétérogénéité. Si la déformation ap-
pliquée est uniforme alors l’eigenstrain ε∗ij et la déformation induite εij sont
uniformes et reliées entre elles par la relation suivante, valide uniquement dans
Ω.
I M
∆Cijkl = Cijkl − Cijkl
Quand l’hétérogénéité contient une déformation inélastique (thermique, . . . )
initiale (inclusion hétérogène), les équations précédentes (Eq.3.6 - Eq.3.8) de-
viennent :
σij0 + σij = Cijkl
I
(ε0kl + εkl − εine
kl ) dans Ω
72
I ∗
Cijkl (ε0kl + εkl − εine ) = Cijkl
M
(ε0kl + εkl − εine
kl − εkl ) (3.10)
Ainsi,
I
σij = Cijkl (Sklmn ε∗∗
mn − ε
ine M
) = Cijkl (Sklmn ε∗∗ ∗∗
mn − εkl ) (3.11)
où
ε∗∗ = ε∗ + εine
1
Dijkl = [Ψ,ijkl − 2νδkl φ,ij − (1 − ν)(δkl φil + δki φ,jl + δjl φ,ik + δli φ,jk )]
8π(1 − ν)
Z (3.15)
Ψ (x) = |x − x0 |dx0
Ω
Z
1
φ(x) = dx0
Ω |x − x0 |
73
Z θ0
0 1
F (θ , k) = dw (3.16)
0 (1 − k )
2 1/2
1/2
a23
0 −1
θ = sin 1− 2
a1
3(a21 − a22 )
k= (3.17)
(a21 − a23 )
a1 , a2 , a3 étant les demi-axes de l’inclusion ellipsoı̈dale. Le tenseur d’Eshelby
Sijkl est obtenu à partir de l’équation Eq. (3.15) :
74
𝜀∗
O 𝑥 ,𝑥
O 𝑥 ,𝑥 O 𝑥 ,𝑥
= + - O 𝑥 ,𝑥
𝜀∗
𝑥 ,𝑥 𝑥 ,𝑥
𝑥 ,𝑥 𝑥 ,𝑥
3 −1 n
nX X2 −1 n
X1 −1
σij (x1 , x2 , x3 ) = Bijkl (x1 − xI1 , x2 − xI2 , x3 − xI3 )ε∗kl (xI1 , xI2 , xI3 )
xI3 =0 xI2 =0 xI1 =0
3 −1 n
nX X2 −1 n
X1 −1
M
Bijkl (x) = Cijmn Dmnkl (x) ∀ x ∈ D-Ω (3.20)
M
Bijkl (x) = Cijmn (Dmnkl (x) − Imnkl ) ∀ x ∈ Ω (3.21)
où : Iijkl = 12 (δil δjk + δik δjl ), est le tenseur identité d’ordre 4.
Dans le cas d’une inclusion unique centrée en (xI1 , xI2 , 0) dans l’espace semi-
infini, la contrainte surfacique σ n en un point (x1 , x2 , x3 ) de la surface est
donnée par :
σ n (x01 , x02 , 0) = − B33kl (x01 − xI1 , x02 − xI2 , −xI3 , )ε∗kl (xI1 , xI2 , xI3 )
(3.22)
− B33kl (x01 − xI1 , x02 − xI2 , xI3 , )ε∗skl (xI1 , xI2 , −xI3 )
75
Dans l’équation 3.19, chaque composante du tenseur Mij () est obtenue par
une double intégration de Fij () sur 2∆x1 × 2∆x2 centrée en (xI1 , xI2 , 0).
Z xI1 +∆x1 Z xI2 +∆x2
Mij (x1 − xI1 , x2 − xI2 , x3 ) = Fij (x1 − x01 , x2 − x02 , x3 )dx01 x02
xI1 −∆x1 xI2 −∆x2
(3.23)
La 3D-FFT est alors utilisée pour accélérer la résolution des sous-problèmes
(1) et (2) et la 2D-FFT pour le (3). On utilise les techniques de wrap around
et de zero-padding pour s’affranchir de l’erreur de périodicité induite par la
FFT [LIU 00b] (cf. chapitre2).
L’effet d’une pression uniforme sur une surface rectangulaire a été ana-
lysé en détail par Love [LOV 52] et Johnson [JOH 85b]. Les coefficients K3p
représentent les coefficients d’influence reliant σ n en un point (x01 , x02 , 0) au
déplacement normal au point (x1 , x2 , 0).
76
3 −1 n
nX X2 −1 n
X 1 −1
σij (x1 , x2 , x3 ) = Bijkl (x1 − xI1 , x2 − xI2 , x3 − xI3 )ε∗kl (xI1 , xI2 , xI3 )
xI3 =0 xI2 =0 xI1 =0
3 −1 n
nX X2 −1 n
X1 −1
+ S
Bijkl (x1 − xI1 , x2 − xI2 , x3 + xI3 )ε∗skl (xI1 , xI2 , −xI3 )
xI3 =0 xI2 =0 xI1 =0
2 −1 n
nX X 1 −1
77
W W
∗
2
o x o x
y dz
α
2 2
z z
(a) (b)
78
O ,
O , O ,
, = , + , - O ,
,
∗
,
, ,
,
79
Méthode de l’inclusion
équivalente
Contraintes élastiques en sous-
couches (3D-FFT)
Prise en compte hétérogénéité
miroir
Détermination de l’eigenstrain
ε*
Détermination de la surcontrainte
due à ε* (3D-FFT) Contraintes en surface
Prise en
compte
Somme des contraintes puis des espace
déplacements
semi-infini
80
Sphere
Subsurface layer
Solid Inclusion
𝑥1
𝑜 30o
𝑥1
30o
𝑥2 = 𝑥2,
𝑥3 𝑥3,
81
Région Géométrie(mm)
Sphère R=31
Solide L1 = L2 = L3 = 60
Géométrie Position
a1 = 0.4a∗ , a2 = a3 = 0.1a∗ dx3 = 0.4a∗
Le modèle éléments finis est un peu lourd, il s’agit d’un modèle à 1.5 × 106 de
ddls. Le choix d’un modèle pareil est justifié principalement pour 2 raisons :
– Pour éviter l’influence des conditions limites sur le problème de contact,
il faut que les dimensions des corps en contact soient au moins 20 fois
82
1 1 1
83
E
k= (3.26)
3(1 − 2ν)
E
µ= (3.27)
2(1 + ν)
La pression équivalente de Hertz correspondant à cette configuration est :
P0 = 4750M P a et la demi-largeur de contact a∗ = 1mm. On étudie donc
l’effet de l’inclusion sur la pression de contact et les contraintes en sous-couche.
L’étude sera faite ici sur une hétérogénéité unique de forme ellipsoı̈dale de
demi-axes a1 = 0.4a∗ , a2 = 0.1a∗ , a3 = 0.1a∗ et située à une profondeur
dx3 = 0.3a∗ .
84
I EI νI EI
Cijkl = δij δkl + (δik δjl + δil δjk ) (3.28)
(1 + ν I )(1 − 2ν I ) 2 × (1 + ν I )
ν I = νM et E I = γE M (3.29)
I M
Cijkl = γCijkl
I M M
∆C = Cijkl − Cijkl = (γ − 1)Cijkl (3.30)
L’équation (3.10) devient :
M
(γ − 1)Cijkl M ∗
Sklmn ε∗mn + Cijkl M 0
εkl = (1 − γ)Cijkl εkl (3.31)
5.5
1.5 γ →0
γ →0 γ = 0.5
γ = 0.5 5
γ =1
γ =1 γ =2
1.25 γ =2 4.5 γ →∞
γ →∞
4
1 3.5
P /P0 P/P0 3
0.75
2.5
2
0.5
1.5
0.25 1
0.5
0 0
−1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5
x1 /a∗ x1 /a∗
(a) (b)
85
4
θ = 00
θ = 300
3.5 θ = 600
θ = 900
θ = 1050
3
2.5
P /P0 2
1.5
0.5
0
−1.5 −1 −0.5 0 0.5 1 1.5
x1 /a∗
86
J :J =J
K:K=K
J :K=K:J =0
J : I2 = I2 : J = I2
K : I2 = I2 : K = 0 (3.33)
µI = µM , k I = δk M (3.34)
I M
∆Cijkl = Cijkl − Cijkl = 3(k I − k M )Jijkl = 3(δ − 1)k M Jijkl (3.35)
L’équation (3.10) devient :
87
M ∗
(δ − 1)k M δij Skkmn ε∗mn + Cijkl εkl = −(δ − 1)k M δij ε0kk (3.37)
où ε0kk représente la partie sphérique de ε0kl .
Les figures 3.16a et 3.16b présentent la pression de contact adimensionnée pour
différentes valeurs du ratio δ.
On observe une augmentation locale de la distribution du champ de pression
quand le matériau tend à devenir incompressible (i.e. quand k I ou δ → ∞).
Le pic de pression est moins élevé dans le cas d’une inclusion parallèle à la
surface (θ = 0, Fig. 3.16a) en comparaison avec le cas (θ = 45◦ , Fig. 3.16b).
Une surpression tout à fait non négligeable peut être générée du seul fait de
l’augmentation du module de compressibilité.
I
Cijkl = 3k I Jijkl + 2µI Kijkl (3.38)
k I = k M , µI = ηµM (3.39)
I M
∆Cijkl = Cijkl − Cijkl = 2(η − 1)µM Kijkl (3.40)
L’équation (3.10) devient :
88
89
1.25
dx3 /a∗ = 0.25
dx3 /a∗ = 0.3
dx3 /a∗ = 0.4
dx3 /a∗ = 0.5
1 dx3 /a∗ = 0.6
dx3 /a∗ = 0.7
P/P0 0.75
0.5
0.25
0
−2 −1 0 1 2
x1 /a∗
(a)
1.25
dx3 /a∗ = 0.25
dx3 /a∗ = 0.3
dx3 /a∗ = 0.4
dx3 /a∗ = 0.5
1 dx3 /a∗ = 0.6
dx3 /a∗ = 0.7
P/P0 0.75
0.5
0.25
0
−2 −1 0 1 2
x1 /a∗
(b)
90
(a) (b)
(c)
Figure 3.12 – Contrainte σ33 /P0 dans le plan x2 = 0 pour différentes valeurs
de γ = E I /E M (avec ν I = ν M = 0.3) pour une inclusion ellipsoı̈dale isotrope
(a1 = 0.4a∗ , a2 = a3 = 0.1a∗ , dx3 /a∗ = 0.4, θ = 30◦ ) ; (a) γ = 0.25, (b) γ = 4, et (c)
γ → ∞.
91
(a) (b)
(c)
Figure 3.13 – Contrainte σ13 /P0 dans le plan x2 = 0 pour différentes valeurs
de γ = E I /E M (avec ν I = ν M = 0.3) pour une inclusion ellipsoı̈dale isotrope
(a1 = 0.4a∗ , a2 = a3 = 0.1a∗ , dx3 /a∗ = 0.4, θ = 30◦ ) ; (a) γ = 0.25, (b) γ = 4, et (c)
γ → ∞.
92
0.5 1
0
0.5
−0.5
−1
σn /P0 τ /P0 0
−1.5
−2
−0.5
Figure 3.14 – Contrainte normale et cisaillement (σn /P0 et τ /P0 ) dans le plan
x2 = 0 pour différentes valeurs de γ = E I /E M , pour une inclusion ellipsoı̈dale
isotrope (a1 = 0.4a∗ , a2 = a3 = 0.1a∗ , dx3 /a∗ = 0.4, θ = 30◦ ) ; (a) σn /P0 et (b)
τ /P0 .
0.5 1.5
0 1
−0.5
0.5
−1
σn /P0 τ /P0 0
−1.5
−0.5
−2
−1
−2.5 γ = 0.25 γ = 0.25
γ=4 γ=4
γ→∞ γ→∞
−3 −1.5
0 60 120 180 240 300 360 0 60 120 180 240 300 360
ξ(degre) ξ(degre)
(a) (b)
Figure 3.15 – Contrainte normale et cisaillement (σn /P0 et τ /P0 ) dans le plan
x2 = 0 pour différentes valeurs de γ = E I /E M , pour une inclusion sphérique (a1 =
a2 = a3 = 0.2a∗ , dx3 /a∗ = 0.3) ; (a) σn /P0 et (b) τ /P0 .
93
1.25 1.5
δ →0 δ →0
δ = 1/2 δ = 1/2
δ =1 δ =1
δ =2 1.25 δ =2
1 δ →∞ δ →∞
1
0.75
P /P0 P /P0
0.75
0.5
0.5
0.25
0.25
0 0
−1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5
x1 /a∗ x1 /a∗
(a) (b)
1.5 5
η →0 η →0
η = 1/2 4.5 η = 1/2
η =1 η =1
1.25 η =2 η =2
η →∞ 4 η →∞
3.5
1
3
0.75 2.5
P /P0 P /P0 2
0.5
1.5
1
0.25
0.5
0 0
−1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5
x1 /a∗ x1 /a∗
(a) (b)
Figure 3.17 – Étude de l’influence du rapport η = µI /µM sur les champs de pression
pour une hétérogénéité ellipsoı̈dale isotrope (a1 = 0.4a∗ , a2 = a3 = 0.1a∗ , dx3 =
0.3a∗ ) ; (a) θ=0 et (b) θ=45◦ .
94
1.5 1.5
Hertz Hertz
Cas 1 Cas 1
Cas 2 Cas 2
1.25 Cas 3 1.25 Cas 3
Cas 4 Cas 4
1 1
P /P0 P /P0
0.75 0.75
0.5 0.5
0.25 0.25
0 0
−1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5
x2 /a∗ x1 /a∗
(a) (b)
95
8
X
φ (x) = E (cn ) (3.42)
n=1
X8
ψ (x) = F (cn ) (3.43)
n=1
avec
et
p
R= c1 2 + c2 2 + c3 2 (3.49)
96
x3
2a3
x1
2a2
x2
2a1
c1 = (x1 − a1 , x2 − a2 , x3 − a3 )
c2 = (x1 + a1 , x2 − a2 , x3 − a3 )
c3 = (x1 + a1 , x2 + a2 , x3 − a3 )
c4 = (x1 − a1 , x2 + a2 , x3 − a3 ) (3.50)
c5 = (x1 − a1 , x2 + a2 , x3 + a3 )
c6 = (x1 − a1 , x2 − a2 , x3 + a3 )
c7 = (x1 + a1 , x2 − a2 , x3 + a3 )
c8 = (x1 + a1 , x2 + a2 , x3 + a3 )
Ces fonctions potentiels ont été implémentées dans le code de calcul par Julien
Leroux [LER 13a] dans sa thèse.
Soit un contact entre une sphère rigide de rayon R = 31mm, soumis à une
force normale FN = 10000N , et un solide élastique semi-infini contenant une
hétérogénéité de forme cubique. Le module de Young et le coefficient de Poisson
de l’espace semi-infini sont respectivement : E M = 210GP a et ν M = 0.3.
Les propriétés matériaux du cube sont anisotropes cubiques, E I = 420GP a,
µI = 14GP a, ν I = 0.45, E I 6= 2 × µI × (1 + ν I )
Les champs de pression et de contraintes sont présentés respectivement dans
les figures 3.21 et 3.22. La particularité des hétérogénéités cubiques réside
dans le fait que les coins du cube agissent parfois comme des concentrateurs
de contraintes.
97
1.4
1.2
0.8
P /P0
0.6
0.4
0.2
0
−1.5 −1 −0.5 0 0.5 1 1.5
x1 /a∗
∆sτ (x1 , x2 )
qτ (x1 , x2 ) = −µ · p(x1 , x2 ) · ∀(x1 , x2 ) ∈ Γsl
k ∆sτ (x1 , x2 ) k
tot ∗ ∗
∆uτ (x1 , x2 ) + u1 (x1 , x2 ) + u2 (x1 , x2 ) − ∆δτ = ∆sτ (x1 , x2 ) ∀(x1 , x2 ) ∈ Γsl
k qτ (x1 , x2 ) k< µ · p(x1 , x2 ) ∀(x1 , x2 ) ∈ Γst
∆uτ (x1 , x2 ) + u∗1 (x1 , x2 ) + u∗2 (x1 , x2 ) − ∆δτ = 0
tot
∀(x1 , x2 ) ∈ Γst
P
q(x1 , x2 )S = Qτ
Γp
où le domaine Γst définit la zone en adhérence (stick) et Γsl la zone de glisse-
ment (slip).
Les champs de déplacements tangentiels u∗i sont à la fois causés par les eigens-
trains et par le champ de pression σn . Les champs u∗i sont donnés par :
XXX 0 0
0 0 0
u∗i (x1 , x2 ) = u
Dikl x1 − x1 , x2 − x2 , x3 ε∗kl x1 , x2 , x3
x3 x2 x1
XXX 0 0
0 0 0
us
+ Dikl x1 − x1 , x2 − x2 , x3 εs∗
kl x1 , x2 , x3
x3 x2 x1
XX 0 0
0 0
+ Kip x1 − x1 , x2 − x2 σ n
x1 , x2
x2 x1
(i = 1, 2) (3.52)
98
(a) (b)
(c)
u
Le tenseur d’ordre 3 Dikl est donné par :
u 1
Dikl = (ψ,kli − 2νδkl φ,i − 2 (1 − νm ) [φ,k δij + φ,j δik ]) (3.53)
8π (1 − ν)
Les fonctions potentiels φ et ψ sont données soit par [MUR 87a] dans le cas
des hétérogénéités de formes ellipsoı̈dales, soit par [MAC 58] dans le cas des
hétérogénéités parallélépipédiques rectangles.
L’algorithme de la figure 3.5 peut toujours s’utiliser, pour faire un couplage
entre hétérogénéités, contact normal et contact tangentiel.
99
∗
2
o x1
x3
100
ψ ∗ ψ −1 ∗
Cijkl (εkl −εine ine
kl ) = Cijkl (εkl −εkl −εkl ) ⇒ σij = Cijkl ((Cklmq ) σmq −εkl ) (3.57)
ψ ψ
σij∗ −Cijkl ((Cklmq )−1 σmq
∗
−ε∗kl ) = Cijkl (Cklmq )−1 (σmq
ine 0
+σmq )−σijine −σij0 , (ψ = 1, 2, . . . , n)
(3.59)
L’équation 3.59 est nécessaire mais pas suffisante pour la résolution du
problème d’hétérogénéités multiples.
Elle n’est pas résolvable en l’état. Il s’agit là d’une équation à deux inconnues.
Il faudrait pour cela trouver une relation entre σ ∗ et ε∗ ainsi qu’entre σ ine et
εine .
Pour ce faire, une méthodologie numérique proposée initialement par [ZHO 11]
puis reprise ensuite par [LER 13a] sera utilisée :
−1 ∗ −1p
I
(Cα,β,γ CM − I d )σα,β,γ I
+ Cα,β,γ ε∗α,β,γ = (I d − Cα,β,γ
I
C M )(σα,β,γ
0
+ σα,β,γ ),
(0 ≤ α ≤ N1 − 1, 0 ≤ β ≤ N2 − 1, 0 ≤ γ ≤ N3 − 1) (3.60)
I
où Cα,β,γ est le tenseur matériau élastique de l’inclusion hétérogène [α, β, γ]
centrée en xα1 , xβ2 , xγ3 .
101
De façon analogue,
N
X 3 −1 N
X 2 −1 N
X1 −1
p
σα,β,γ = Bα−ξ,β−ζ,γ−ϕ εpξ,ζ,ϕ
ϕ=0 ζ=0 ξ=0
(0 ≤ α ≤ N1 − 1, 0 ≤ β ≤ N2 − 1, 0 ≤ γ ≤ N3 − 1) (3.62)
3 −1 N
NX 2 −1 N1 −1
−1
X X
I
Cα,β,γ CM −I d
Bα−ξ,β−ζ,γ−ϕ ε∗ξ,ζ,ϕ + Cα,β,γ
I
ε∗α,β,γ
ϕ=0 ζ=0 ξ=0
3 −1 N2 −1 N 1 −1
N
!
−1
X X X
= I d − Cα,β,γ
I
CM Bα−ξ,β−ζ,γ−ϕ εpξ,ζ,ϕ + σα,β,γ
0
,
ϕ=0 ζ=0 ξ=0
(0 ≤ α ≤ N1 − 1, 0 ≤ β ≤ N2 − 1, 0 ≤ γ ≤ N3 − 1) (3.63)
Ax = b (3.64)
102
g0 = Aε∗0 − b
w 0 = g0 (3.65)
v = Awp−1
ρp−1 = −(gp−1 .wp−1 )/([Link]−1 )
ε∗p = ε∗p−1 + ρp−1 wp−1 (3.66)
gp = gp−1 + ρp−1 v
if (gp .gp )/(b.b) < 2 then
F in
end if
γp−1 = −(gp .v)/([Link]−1 )
wp = gp + γp−1 wp−1
103
g0 = Ax0 − b
d0 = g0
g̃ 0 = g0 (3.67)
r0 = g0
v = Adp−1
(rp−1 .g̃ 0 )
ρp−1 = −
(v.g0 )
qp−1 = rp−1 + ρp−1 v
w = Aqp−1
([Link]−1 )
ωp−1 = −
kwk2
rp = qp−1 + ωp−1 w
∗ ∗
εp = εp−1 + ρp−1 dp−1 + ωp−1 qp−1
if (rp .rp )/(b.b) < 2then (3.68)
F in
end if
1 rp .g̃ 0
γp−1 = −
ωp−1 v.g̃ 0
dp = rp + γp−1 (dp−1 + ωp−1 .v)
Algorithme Orthodir
g0 = Ax0 − b
w = g0
v = Aw (3.69)
w
w1 =
kvk
v
Aw1 =
kvk
104
w = Awp−1
v = Aw
for i = 1 to p − 1
αi = ([Link] )
w = w − αi .vi
v = v − αi .Awi
end for
1
wp = .w
kvk
1
Awp = .v (3.70)
kvk
ρp = −(gp−1 .Awp )
ε∗p = ε∗p−1 + ρp wp
gp = gp−1 + ρp Awp
if (gp .gp )/(b.b) < 2 then
F in
end if
Il s’agit bien sûr d’une solution approchée, compte tenu du fait que l’eigens-
train est supposée constante dans chaque hétérogénéité. Cependant une bonne
précision peut être obtenue si la discrétisation est assez fine et les dimensions
des hétérogénéités pas trop grandes (≤ 0.4a∗ ).
105
O O
(a) (b)
Figure 3.24 – Contact entre une sphère rigide et un massif revêtu.
revêtement dur. Le rayon R de l’indenteur est fixé égal à 0.21mm. Les pro-
priétés matériaux du substrat sont choisies telles que Es = 210000M P a et
νs = 0.3. La pression équivalente de Hertz correspondant à la configura-
tion homogène sans revêtement est P0 = 2106M P a et la demi-largeur de
contact a∗ = 0.003mm. L’épaisseur ec du revêtement est fixé à 0.25a∗ . Dans
la modélisation semi-analytique, le revêtement sera représenté par une série
d’hétérogénéités de formes cubiques (cf. Fig.3.24) isotropes de propriétés Ec , νc
s’étendant sur une profondeur ec . Ici, le formalisme d’Eshelby sera donc ap-
pliqué aux hétérogénéités de formes cubiques. La méthodologie présentée dans
les sections (3.7, 3.7.2,3.7.3) sera utilisée pour prendre en compte les influences
mutuelles entre les différents cuboı̈des.
La comparaison sera faite ici entre le code semi-analytique et un code multi-
grille développé par Hugo Boffy durant sa thèse [BOF 12b, BOF 14]. L’avan-
tage, par rapport à la solution classique de O’Sullivan, est qu’il est possible
à la fois de comparer les champs de pression et les champs de contraintes
en sous-couche. La figure 3.25 présente la comparaison du champ de pression
entre le code semi-analytique et la solution basée sur les méthodes multi-grilles
[BOF 12b, BOF 14]. La comparaison des champs de contraintes le long des
axes X1 et X3 est présentée aux figures 3.26 et 3.27. Une bonne corrélation
existe entre les deux solutions.
Les nouveaux outils ainsi introduits permettent de prendre en compte des ma-
trices non symétriques. Il devient alors possible de résoudre des problèmes
de fretting (coefficient de frottement non nul) en présence de plusieurs
hétérogénéités. Au-delà du cas d’un matériau revêtu, ce formalisme per-
met de traiter le cas des géométries d’hétérogénéités complexes. Il est pos-
sible de décrire les fibres d’un matériau composite tissé par une succession
d’hétérogénéités de formes cubiques (représentation schématique à la figure
3.28a) ou comme une succession d’hétérogénéités cubiques et ellipsoı̈dales
106
1.4
Code SA, sans influences mutuelles
Code SA, avec influences mutuelles
1.2 multigrilles, Boffy et al.(2012)
0.8
P /P0
0.6
0.4
0.2
0
−1.5 −1 −0.5 0 0.5 1 1.5
x1 /a∗
107
1.5 1.5
Code SA, sans influences mutuelles Code SA, sans influences mutuelles
Code SA, avec influences mutuelles Code SA, avec influences mutuelles
1 multigrilles, Boffy et al.(2012) 1 multigrilles, Boffy et al.(2012)
0.5 0.5
σ11/P0 σ22/P0
0 0
−0.5 −0.5
−1 −1
−1.5 −1.5
−1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5
x1 /a∗ x1 /a∗
(a) (b)
1.5 1.5
Code SA, sans influences mutuelles Code SA, sans influences mutuelles
Code SA, avec influences mutuelles Code SA, avec influences mutuelles
1 multigrilles, Boffy et al.(2012) 1 multigrilles, Boffy et al.(2012)
0.5 0.5
−0.5 −0.5
−1 −1
−1.5 −1.5
−1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5
x1 /a∗ x1 /a∗
(c) (d)
108
1.5 1.5
Code SA, sans influences mutuelles Code SA, sans influences mutuelles
Code SA, avec influences mutuelles Code SA, avec influences mutuelles
1 multigrilles, Boffy et al.(2012) 1 multigrilles, Boffy et al.(2012)
0.5 0.5
σ11/P0 σ22/P0
0 0
−0.5 −0.5
−1 −1
−1.5 −1.5
0 0.5 1 0 0.5 1
x3 /a∗ x3 /a∗
(a) (b)
1.5 1.5
Code SA, sans influences mutuelles Code SA, sans influences mutuelles
Code SA, avec influences mutuelles Code SA, avec influences mutuelles
1 multigrilles, Boffy et al.(2012) 1 multigrilles, Boffy et al.(2012)
0.5 0.5
−0.5 −0.5
−1 −1
−1.5 −1.5
0 0.5 1 0 0.5 1
x3 /a∗ x3 /a∗
(c) (d)
109
(a) (b)
110
111
3 −1 n
nX X2 −1 n
X 1 −1
σij (x1 , x2 , x3 ) = Bijkl (x1 − xI1 , x2 − xI2 , x3 − xI3 )ε∗kl (xI1 , xI2 , xI3 )
xI3 =0 xI2 =0 xI1 =0
2 −1 n
nX X1 −1
112
O ,
O ,
, = O
, - ,
,
∗
, , ,
6
T2 /T1
T3 /T1
5
ratio 3
0
0 25 50 75 100 125 150
N
113
114
Mémoire Centrale
Réseau de communication
Déformation Déformation Déformation
macroscopique liée au macroscopique liée au macroscopique liée au
problème de contact problème de contact problème de contact
Mémoire locale 2 :
Mémoire locale1 :
(informations relatives Mémoire locale 3 :
(coefficients élastiques,
aux familles d’ (informations relatives
données d’entrées du
hétérogénéités 1 à N/2) aux familles d’ hé-
contact, …)
térogénéités N/2 à N)
Réseau
Figure 3.32 – Exemple de communication
d’architecture à mémoire distribuée.
115
Tp
Tnp = Ts + (3.72)
np
Bien évidemment il est difficile voire impossible d’atteindre Tnp . Cela suppo-
serait que la part parallélisable du code soit parfaitement répartie entre les
différents processeurs sans que l’exécution parallèle n’entraine de surcoût.
La partie séquentielle de la résolution de contact sera par exemple : la lec-
ture et l’écriture de fichiers, la résolution du problème normal et tangentiel.
Cependant les algorithmes de calcul de contraintes et de prise en compte
des hétérogénéités sont fortement parallélisables.
– Équilibrage des tâches. Un autre critère de performance est la régularité en
termes de temps d’exécution des différentes tâches de la partie parallélisable
du code. En effet, si une tâche requiert un temps d’exécution plus élevé que
les autres, c’est elle qui va donc déterminer le temps minimal d’exécution en
parallèle. Les divers algorithmes de la partie parallélisable du code doivent
être parallélisés de manière à avoir un bon équilibrage des tâches sur les
différents processeurs.
– Granularité. La granularité d’une tâche parallèle peut se définir comme le
nombre d’opérations arithmétiques que la tâche réalise et le nombre de
données que la tâche doit recevoir. Concrètement il s’agit du rapport entre le
temps de calculs et celle des échanges entre processus. En effet, l’exécution en
mode parallèle entraı̂nera des coûts supplémentaires qui vont réduire l’effica-
116
cité. Sur des machines à mémoire partagée, ces surcoûts sont faibles, liés aux
fonctions du système d’exploitation qui activent les différents processeurs et
gèrent la synchronisation entre processus. Sur un système à mémoire dis-
tribuée, les surcoûts sont principalement liés aux transferts de données entre
processeurs qui induisent des délais supplémentaires lors de l’exécution en
mode parallèle. De manière simplifiée, ne seront parallélisées que les boucles
parallélisables (bien sûr) ayant un temps d’exécution élevé.
– Extensibilité. Un algorithme ou une méthode seront dits extensibles si leur
efficacité ne décroı̂t pas avec le nombre de processeurs. Dans le cas d’un en-
vironnement OpenMP, on peut établir l’extensibilité en exécutant le même
calcul sur un nombre croissant de processeurs et mesurer la décroissance
du temps de calcul. Aucun algorithme ne sera extensible à l’infini, ce qui
compte c’est d’avoir une bonne extensibilité dans la gamme de processeurs
disponibles. Dans le cas de la parallélisation dans l’environnement à mémoire
partagée, OpenMP, l’extensibilité est bonne entre 2 et 16 processeurs. Ce
qui est très satisfaisant, puisque le code semi-analytique se veut un code
de calcul Bureau d’études (BE) léger devant tourner sur des ordina-
teurs de bureaux ayant entre 4 et 8 Cœurs. Par contre dans le cas MPI,
l’extensibilité dépend fortement du type de problème à résoudre (couplage
normal/tangentiel, nombre d’hétérogénéités, . . . ).
Le code de contact fait intervenir des calculs répétitifs sur de grands nombres
de données (calcul du tenseur Bijkl ). Ces grands nombres de données appa-
raissent sous forme de boucles de programme portant sur des tableaux.
Lorsque que le code de calcul effectue par exemple des boucles de programme,
les différents processeurs peuvent se partager le travail, en partitionnant l’en-
semble ou une partie des itérations. Évidemment toutes les boucles d’un pro-
gramme ne sont pas parallélisables. Il y a une série de règles assez strictes à
respecter.
– Dépendance
L’exécution en parallèle signifie que les différentes itérations seront effectuées
dans un ordre qui est fonction du nombre de processeurs, du mode de
découpage, et qui de fait est très aléatoire. La parallélisation ne modifie
pas le résultat du calcul si et seulement si il n’existe pas de dépendance
entre deux ou une série d’instructions. Il existe deux types de dépendances :
la dépendance de données et la dépendance de sorties.
la dépendance de données : Deux instructions ou une série d’instructions
InstA et InstB présentent une dépendance de données, lorsque :
117
σ
bij (s1 , s2 , s3 ) = B ε∗kl (s1 , s2 , s3 ) + B
bijkl (s1 , s2 , s3 )b ε∗skl (s1 , s2 , s3 )
bijkl (s1 , s2 , s3 )b
+M cij (s1 , s2 , s3 )b σ n (s1 , s2 )
(3.77)
118
σ
bij (s1 , s2 , s3 ) = B ε∗kl (s1 , s2 , s3 ) + M
bijkl (s1 , s2 , s3 )b cij (s1 , s2 , s3 )b
τzz (s1 , s2 )
(3.78)
+Mczx (s1 , s2 , s3 )b τzx (s1 , s2 ) + M czy (s1 , s2 , s3 )bτzy (s1 , s2 )
ij ij
Les tenseurs Bijkl et ε∗kl font intervenir des transformées de Fourier tridimen-
sionnelles (3D-FFT). Les 3D-FFT sont très coûteuses en temps de calcul com-
parées aux 2D-FFT. Il convient de paralléliser en priorité les routines 3D-FFT.
Le tenseur Bijkl possède la symétrie en ij, en kl mais pas la symétrie majeure.
Il existe alors 36 termes distincts pour le tenseur Bijkl . Il est alors possible de
réaliser simultanément la transformée de Fourier de ces 36 (Bijkl ) + 6(ε∗kl )
termes. Il s’agit alors (section 3.9.2) d’une portion de code fortement pa-
rallélisable. Ainsi quand l’on est présence de 42 processeurs, la loi d’Amdhal
appliquée à cette portion de code permet de dire :
Transformée inverse
σij (x1 , x2 , x3 ) = F F T −1 (b
σij (s1 , s2 , s3 )) (3.80)
119
calcul)
3.5
2.5
de
2
Gain(T emps
1.5
0.5
0
0 1 2 3 4
Cœurs
120
121
Réseau de communication
Déformation Déformation Déformation
macroscopique liée au macroscopique liée au macroscopique liée au
problème de contact problème de contact problème de contact
Mémoire locale 2 :
Mémoire locale1 :
(informations relatives Mémoire locale 3 :
(coefficients élastiques,
aux familles d’ (informations relatives
données d’entrées du
hétérogénéités 1 à N/2) aux familles d’ hé-
contact, …)
térogénéités N/2 à N)
Réseau
Figure 3.34 – Mise en oeuvre de communication
la parallélisation MPI code semi-analytique.
3.10 Synthèse
Ce chapitre s’inscrit dans la continuité des thèses de Julien Leroux [LER 13a]
et de Benjamin Fulleringer [FUL 11] qui se sont intéressés au cas des
hétérogénéités élastiques isotropes noyées dans un massif semi-infini élastique
isotrope. Dans ce chapitre la prise en compte des hétérogénéités anisotropes ou
isotropes inclinées dans le cas du contact entre corps élastiques isotropes a été
présentée. Le modèle a été validé par comparaison avec un modèle éléments-
finis réalisé avec le logiciel commercial Abaqus v6.11. Une étude fine, de l’effet
du module de compressibilité/cisaillement de l’hétérogénéité sur le problème
de contact, a été réalisée grâce à une décomposition du tenseur d’élasticité
C sur la base des tenseurs isotropes J et K. Les résultats de cette ana-
lyse ont été utilisés en bureau d’études (BE) pour étudier la nocivité de
défauts dans les billes céramiques de roulements hybrides. La méthode per-
met d’accéder aux champs de contraintes normales/tangentielles à l’interface
hétérogénéité/matrice. Ces variables s’avèrent utiles dans la mise en place des
lois de décohésions hétérogénéité/matrice. Dans le cas des hétérogénéités ani-
sotropes en particulier, il est indispensable de prendre en compte la bonne
orientation du repère d’anisotropie par rapport au repère du contact. Outre
les hétérogénéités de forme ellipsoı̈dale, le modèle pourrait également prendre
en compte des hétérogénéités de formes parallélépipédiques.
Une méthode numérique permettant de prendre en compte les interactions
mutuelles entre plusieurs hétérogénéités a été ensuite présentée. De nouveaux
122
123
124
Théorie du contact
viscoélastique
Sommaire
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . 127
4.2 Généralités sur les matériaux viscoélastiques . . . . 130
4.3 Contact normal sur massif homogène viscoélastique 134
4.3.1 Equations du modèle . . . . . . . . . . . . . . . . . . . . 134
4.3.2 Mise en œuvre numérique . . . . . . . . . . . . . . . . . 137
4.3.3 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . 138
4.3.4 Quelques résultats . . . . . . . . . . . . . . . . . . . . . 141
4.4 Contact normal sur massif viscoélastique hétérogène 146
4.4.1 Equations du modèle . . . . . . . . . . . . . . . . . . . . 146
4.4.2 Validation. . . . . . . . . . . . . . . . . . . . . . . . . . 151
4.4.3 Résultats . . . . . . . . . . . . . . . . . . . . . . . . . . 154
125
126
4.1 Introduction
Le contact viscoélastique, contrairement au cas élastique, est très peu étudié
dans la littérature. Cependant, de nombreuses applications et solutions indus-
trielles mettent en œuvre la problématique du contact entre deux matériaux
dont l’un au moins est viscoélastique. Les applications industrielles les plus
répandues sont le contact pneu/chaussée, le phénomène de fretting dans le
cas des matériaux composites à matrice organique, la biomécanique, l’étude
de revêtements viscoélastiques pour des applications pétrolières, automo-
biles, aéronautiques et aérospatiales. Dans le cas du contact pneu-chaussée,
les industriels du pneumatique cherchent à mieux comprendre et à mieux
modéliser l’adhérence (frottement, usure) et ses facteurs influents (texture,
vitesse,. . . ). En effet, la viscoélasticité de la gomme du pneu introduit un
différentiel de chargement entre l’avant et l’arrière de la zone de contact, in-
duisant un couple résistant. La détermination de cette résistance au roule-
ment est une problématique qui intéresse fortement les industriels du pneu-
matique. Les concepteurs de chaussée sont quant à eux plutôt intéressés par
la compréhension des mécanismes de dégradations des chaussées, donc aux
champs de contraintes en surface et en sous-couche. Cela passe bien évidement
par une modélisation fine du contact entre matériaux viscoélastiques. Certaines
entreprises comme Hexcel, leader dans la fabrication des matériaux compo-
sites, sont fortement intéressées par les problématiques de fretting dans le cas
de ces matériaux. L’un des objectifs de ces études, est de comprendre l’effet
de la matrice viscoélastique sur le comportement à long terme en fretting de
divers matériaux composites [TER 11, TER 09]. Dans la biomécanique, cer-
tains industriels sont intéressés par la tenue et le confort des implants ou des
prothèses. Dans le cas des implants orthopédiques [GER 05] par exemple, le
matériau utilisé est un alliage T i − 6Al − 4V qui est en contact avec des tis-
sus mous dont le comportement est supposé viscoélastique linéaire. Enfin dans
le monde aéronautique, les revêtements viscoélastiques sont particulièrement
appréciés au niveau des structures froides soumises à des sollicitations de fret-
ting.
La plupart des modèles de la littérature qui se sont intéressés au contact
viscoélastique peuvent être subdivisés en 3 grandes parties : l’indentation, le
roulement/glissement et le fretting. Contrairement au contact élastique dont
les premiers travaux remontent à 1885, le premier modèle traitant du contact
viscoélastique a été publié en 1960 par Lee et Radok [LEE 60]. Ils ont donné
les solutions en terme de champs de pression ainsi que l’évolution de la zone de
contact entre une sphère rigide et un massif viscoélastique linéaire. La solution
fournie par [LEE 60] n’est valable que pour une augmentation monotone de la
zone de contact. Une pression négative est obtenue lorsque la zone de contact
décroit, résultat qui n’est pas réaliste. Hunter [HUN 60] et Graham [GRA 67]
ont proposé une amélioration du modèle d’origine de Lee et Radok [LEE 60]
127
128
129
mentionner que l’approche développée dans ce chapitre peut traiter le cas des
matériaux viscoélastiques linéaires avec un ou plusieurs temps de relaxation,
peut prendre en compte l’effet de la rugosité des surfaces sur le problème de
contact, ceci quel que soit le chargement appliqué : force/déplacement im-
posé, chargement de type rampe, Heavy-side, triangulaire, quelconques, . . . Le
modèle fournit des solutions en termes de champs de pression, de cisaille-
ments, de coefficient de frottement mais surtout permet d’obtenir des champs
de contraintes en sous-couche. Ces champs de contraintes représentent des
variables indispensables dans l’étude de l’endommagement en sous-couche.
Contrairement à une grande majorité des modèles existants dans la littérature,
l’approche développée ici fournit des résultats aussi bien en régime transitoire
qu’en régime permanent. Le chapitre sera subdivisé en 3 grandes parties. La
première partie s’intéressera au problème d’indentation entre deux massifs
viscoélastiques homogènes ou hétérogènes. L’effet des hétérogénéités sur le
problème de contact sera étudié. La deuxième partie s’intéressera au problème
de roulement/glissement. La dernière partie s’intéressera à des problématiques
de fretting dans le cas du contact viscoélastique. Des comparaisons entre
les résultats obtenus par le modèle et la méthode des éléments finis seront
présentées régulièrement tout au long de ce chapitre. Le gain en terme de
temps de calculs entre le modèle et la méthode des éléments finis peut at-
teindre 110, dans le cas viscoélastique homogène. Il s’agit d’un modèle très
rapide.
Les développements présentés dans ce manuscrit se limiteront au cadre de la
viscoélasticité linéaire.
130
viscoélasticité linéaire
131
132
" n
#
X
R(t) = µ0 + µi exp(−t/τi ) H(t) (4.6)
i=0
Z t
dσ(ξ)
ε(t) = J(t − ξ) dξ (4.8)
0 dξ
133
W(t)
( )( ),
X2 o X1
( )( ),
X3
Figure 4.5 – Illustration du contact normal viscoélastique dans le cas d’une confi-
guration sphère/plan.
134
135
(1) (2)
u3 (x1 , x2 , t) = u3 (x1 , x2 , t) + u3 (x1 , x2 , t) (4.13)
∞ ∞ t
(1 − ν1 )p(x01 , x02 , ξ)
Z Z Z
u3 (x1 , x2 , t) = ( p 0 2 0 2
J (1) (t − ξ)
−∞ −∞ 0 π (x1 − x1 ) + (x2 − x2 )
0 0
(1 − ν2 )p(x1 , x2 , ξ) (2) ∂p(x01 , x02 , ξ) 0 0
+ p 0 0
J (t − ξ)) dx1 x2 dξ (4.14)
π (x1 − x1 )2 + (x2 − x2 )2 ∂ξ
Z ∞ Z ∞ Z t
0 0
u3 (x1 , x2 , t) = (G(1) (x1 − x1 , x2 − x2 , t − ξ) +
−∞ −∞ 0
0 0 ∂p(x01 , x02 , ξ) 0 0
G(2) (x1 − x1 , x2 − x2 , t − ξ)) dx1 x2 dξ (4.15)
∂ξ
où, G(1,2) (x1 , x2 , t) est le tenseur de Green viscoélastique, représentant le
déplacement normal de la surface viscoélastique si celle-ci était soumise à une
force unité. G(1,2) (x1 , x2 , t) peut être décomposé en une partie spatiale et une
partie temporelle.
(1 − ν1,2 ) (1,2)
G(1,2) (x1 , x2 , t) = p 2 2
J (t) = G(x1 , x2 )×(1−ν1,2 )×J (1,2) (t) (4.16)
π x1 + x 2
R∞
En discrétisant la forme intégrale −∞ ()dξ et en utilisant la méthode des
différences finies pour évaluer la dérivée ∂p ∂ξ
, l’équation 4.15 peut se réécrire
sous la forme :
Z ∞Z ∞
u3 (x1 , x2 , α∆t) = G(x1 − x01 , x2 − x02 )dx01 x02 ×
−∞ −∞
α
X
(1 − ν1 )J (1) [(α − k)∆t] + (1 − ν2 )J (2) [(α − k)∆t]
k=0
×[p(x01 , x02 , k) − p(x01 , x02 , k − 1)] (4.17)
où ∆t est le pas de la discrétisation temporelle.
Avec cette expression générale de u3 (x1 , x2 , t), il est possible de remonter aux
différents cas limites : contact élastique/élastique, élastique/viscoélastique, ri-
gide/viscoélastique (cas le plus étudié dans la littérature).
A l’aide des équations développées dans cette partie, on peut résoudre
complètement le problème de contact normal entre deux matériaux
viscoélastiques et ainsi obtenir les champs de pression à chaque pas de la
discrétisation temporelle.
136
N1 X
X N2
0(2)
σij (x1 , x2 , x3 , α∆t) = p(x01 , x02 , α∆t)(Cijp (x1 − x01 , x2 − x02 , x3 , ν2 ))
x01 =1 x02 =1
(4.18)
137
+(1 − ν2 )J (2)
[(α − k)∆t]) × p(i0 , j 0 , k − 1) (4.20)
4.3.3 Validation
Pour la validation, on considérera un contact normal sphère/plan. Les deux
corps en contact sont supposés homogènes. Le diamètre de la sphère est D =
62mm et une force normale d’intensité W (t) = 10000N ×H(t) est appliquée. a∗
et P0 sont respectivement le rayon de contact et la pression maximum de Hertz
correspondant au cas homogène élastique quand le module de cisaillement
instantané R(0) = µ est utilisé (cf. Eq. 4.22). Autrement dit le rayon de contact
138
1 2
0.9
1.75
0.8
0.7
1.5
µ × J (t)
0.6
R(t)/µ
0.5 1.25
0.4
1
0.3
0.2
0.75
0.1
0 0.5
0 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1
t/τ t/τ
(a) (b)
1 t
J(t) = + (4.21)
µ η
et :
R(t) = µ exp(−t/τ ) (4.22)
où µ est le module purement élastique, η la viscosité, et τ = η/µ le temps
de relaxation. Le coefficient de Poisson du massif viscoélastique est supposé
constant ν M =0.3. Le module de cisaillement instantané est R(0) = µ =
80.77GP a et le temps de relaxation τ = 25 secondes.
Trois calculs seront réalisés :
1. Contact sphère viscoélastique/massif viscoélastique. La sphère et le massif
seront considérés comme viscoélastiques. En considérant cette configura-
tion, la demi-largeur de contact et la pression de Hertz équivalentes à
t = 0 sont respectivement, a∗ (t = 0) = a∗ = 1.263mm et P max(t = 0) =
P0 = 2992.894M P a
2. Contact sphère élastique/massif viscoélastique. Le massif sera considéré
viscoélastique et la sphère élastique. Ceci revient simplement à faire
tendre uniquement τ → ∞. Pour cette configuration a∗ (t = 0) = a∗ =
1.263mm et P max(t = 0) = P0 = 2992.894M P a
3. Contact sphère rigide/massif viscoélastique. La sphère sera considérée
rigide. Ceci revient dans nos formulations à faire tendre τ → ∞ et µ →
∞. a∗ (t = 0) = a∗ = 1.00mm et P max(t = 0) = P0 = 4750.923M P a
139
Les figures 4.7a, 4.7b, 4.7c présentent respectivement les comparaisons des
champs de pression dans le cas du contact sphère viscoélastique/massif
viscoélastique, sphère élastique/massif viscoélastique, sphère rigide/massif
viscoélastique.
Une bonne corrélation entre le code semi-analytique et le code éléments finis
commercial (Abaqus v6.12) peut être observée.
Méthode EF 3D Code SA
Temps CPU 76500s(∼ 21h15min) 672s(∼ 12min)
140
2000
Mesures experimentales
Decomposition en serie de Pronys
1800
R(t)
1600
1400
1200
0 200 400 600 800 1000
t(s)
t t
RP M M A (t) = 1429.71 + 184.62 exp(− ) + 191.06 exp(− ) (4.23)
8.93 117.96
On peut donc en déduire la fonction de fluage en utilisant la relation 4.5 :
141
1.2 1.6
1.4
1
1.2
W (mN )
0.8
W (t)/Wmax
0.6 0.8
0.6
0.4
0.4
0.2
0.2
0 0
0 0.5 1 1.5 2 0 500 1000 1500 2000
t/T t(sec)
(a) (b)
142
100
Simulation(Code SA)
Cas experimental
80
60
W (N )
40
20
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07
δ(mm)
0.6
0.5
0.4
0.3
0.2
0.1
0
0 0.2 0.4 0.6 0.8 1 1.2
x1 /a∗
143
100 λ = 0.1
λ = 0.2
λ=1
λ = 10
80 λ = 100
60
W (N )
40
20
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07
δ(mm)
1.5 ω = 0.0092s−1
ω = 0.3s−1
1.4
1.3
W (t)/Wlim
1.2
1.1
1
0.9
0.8
0.7
0.6
0.5
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4
δ(t)/δlim
Figure 4.13 – Courbe d’indentation dans le cas d’un chargement sinusoı̈dal pour
différentes pulsations d’excitation.
144
1.5 1.5
Chargement Chargement
Reponse du contact Reponse du contact
1.5 1.5
W (t)/Wlim
δ(t)/δlim
W (t)/Wlim
δ(t)/δlim
1 1
1 1
plus approfondie montre qu’il existe une pulsation de résonance pour laquelle
cette énergie est maximale.
L’analyse des champs de pression (Fig.4.14) montre qu’il y a un déphasage
entre le chargement et la réponse du contact. L’effet visqueux entraı̂ne donc
un retard de la réponse du contact sur le chargement.
145
146
A-A
W(t) W(t)
( )( ),
( )( ),
∗
2
X2 o X1
o x1
A A dx
M Ɵ
( )( ),
2 2
X3
x3 ( )( ),
(a) (b)
Figure 4.15 – Illustration du contact normal dans le cas d’une sphère et d’un massif
viscoélastique hétérogène. (a) Vue 3D et (b) Coupe dans le plan X2 = 0 2D.
M
Considérons une matrice infinie D avec des propriétés élastiques Cijkl (α∆t)
contenant une inclusion de forme ellipsoı̈dale, cubique, cylindrique, . . . Ω de ten-
I
seur élastique Cijkl . L’ensemble est soumis à un champ uniforme ε0 à l’infini. La
présence du domaine inclusionnaire modifie grandement les champs de déformations.
En appliquant le formalisme d’Eshelby à un instant t = α∆t, on peut remplacer
l’hétérogénéité ellipsoı̈dale par une inclusion ayant les mêmes propriétés élastiques
que la matrice à l’instant t = α∆t mais sujette à une déformation imaginaire libre
de contrainte, eigenstrain ε∗ donnée par :
I
Cijkl (ε0kl + εkl ) = Cijkl
M
(α∆t) × (ε0kl + εkl − ε∗kl ) dans Ω (4.30)
L’équation 4.30 est une équation nécessaire au problème d’équivalence entre in-
clusion et hétérogénéité. La relation entre la déformation compatible εij et l’eigens-
train ε∗ij est donnée par :
εij = Sijkl (α∆t) × ε∗kl (4.31)
où Sijkl (α∆t) est le tenseur d’Eshelby à l’instant t = α∆t. Le tenseur d’Eshelby
classique ne dépend que du coefficient de Poisson de la matrice. Compte tenu du
fait que dans nos approches, le coefficient de Poisson est supposé constant, on montre
147
Données d’entrée :
• Propriétés géométriques et
matériaux
• Chargement
• Simulation pour ∈ [0, ]
Discrétisation temporelle en
pas de temps ∆ = ; ∈
0…….
=0
Détermination eigenstrain ε*
Contraintes en surface
Détermination de l’eigenstress
(3D-FFT)
Déplacements ∗ et
= +1 Déplacements induits par la
contraintes induits par les
présence des hétérogénéités
contraintes en surfaces.
U ∗ , U ∗ (éventuellement)
Prise en
compte
Somme des contraintes puis des espace
déplacements
semi-infini
= ?
148
O ,
O ,
, = O
, - ,
,
∗
, , ,
que :
où :
I M
∆Cijkl (α∆t) = Cijkl − Cijkl (α∆t)
Comme cela a été mentionné au chapitre précédent, les solutions d’Eshelby clas-
siques ne sont valables qu’en espace infini. Le problème de contact par contre fait
intervenir des espaces semi-infinis. On peut s’inspirer de la décomposition proposée
par Chiu [CHI 78a] (cf. chapitre 3) ou de la nouvelle méthode de décomposition
présentée à la section 3.8 du chapitre précédent. La nouveauté ici, réside dans le
fait que la méthode de décomposition sera utilisée à chaque pas de temps de la
discrétisation temporelle α∆t (Fig.4.17). Pour une simulation avec plusieurs pas de
temps, le gain en terme de temps de calcul devient très vite conséquent lorsque la
nouvelle méthode de décomposition est utilisée.
Le champ de contrainte en tout point de l’espace maillé par n1 × n2 × n3 pa-
rallélépipèdes est donné par :
149
3 −1 n
nX X2 −1 n
X1 −1
σij (x1 , x2 , x3 , α∆t) = Bijkl (x1 − xI1 , x2 − xI2 , x3 − xI3 , α∆t)ε∗kl (xI1 , xI2 , xI3 , α∆t)
xI3 =0 xI2 =0 xI1 =0
2 −1 n
nX X 1 −1
M
Bijkl (x, α∆t) = Cijmn (α∆t)(Dmnkl (x, α∆t) − Imnkl ) ∀ x ∈ Ω (4.36)
où Iijkl = 12 (δil δjk + δik δjl) est le tenseur identité d’ordre 4. Le tenseur d’ordre 2
Mij (voir Annexe A) et le tenseur d’ordre 4 Dijkl dépendent uniquement du coeffi-
cient de Poisson qui est supposé constant. Ainsi :
et
Mij (x, α∆t) = Mij (x) (4.38)
1
Dijkl = [Ψ,ijkl − 2νδkl φ,ij − (1 − ν)(δkl φil + δki φ,jl + δjl φ,ik + δli φ,jk )]
8π(1 − ν)
(4.41)
Ψ (x) et φ(x) sont les potentiels harmoniques et bi-harmoniques (cf. chapitre 3) :
150
Pour une inclusion unique centrée en (xI1 , xI2 , xI3 ) dans l’espace semi-infini, la
contrainte normale σ n au point (x,1 , x,2 , 0) est obtenue par :
σ n (x01 , x02 , 0, α∆t) = − B33kl (x01 − xI1 , x02 − xI2 , −xI3 , α∆t)ε∗kl (xI1 , xI2 , xI3 , α∆t)
(4.42)
− B33kl (x01 − xI1 , x02 − xI2 , xI3 , α∆t)ε∗skl (xI1 , xI2 , −xI3 , α∆t)
Déplacement normal de la surface viscoélastique
Les déplacements normaux u∗3 dus aux hétérogénéités sont uniquement générés
par le champ de pression σ n permettant de créer une surface libre (massif semi-
infini). A la différence du cas élastique, ces déplacements s’écrivent sous la forme :
2 −1 n
nX 1 −1
X p
u∗3 (x1 , x2 , α∆t) = K 3 (x1 − x01 , x2 − x02 ) ×
x02 =0 x01 =0
α
X
(1 − ν2 )J (2) [(α − k)∆t] × [σ n (x01 , x02 , k∆t) − σ n (x01 , x02 , (k − 1)∆t)] (4.43)
k=0
4.4.2 Validation.
Pour la validation, on considérera un contact normal entre une sphère rigide et
un matériau viscoélastique contenant une ou plusieurs hétérogénéités. Le diamètre
de la sphère rigide est D = 62mm et une force normale W (t) = 10000N × H(t) est
appliquée. a∗ et P0 sont respectivement le rayon de contact et la pression maximum
de Hertz correspondant au cas homogène élastique quand le module de cisaillement
instantané est utilisé R(0) = µ.
Propriétés matériaux
Afin de valider notre modèle, une comparaison avec la méthode des éléments
finis est réalisée en utilisant la version v6.11-2 du code EF commercial Abaqus. Le
contact entre une sphère rigide et un massif viscoélastique hétérogène contenant une
hétérogénéité ellipsoı̈dale élastique isotrope est simulé. Une technique de partition
particulière a été utilisée pour réaliser l’interface hétérogénéité/matrice. Les pro-
priétés géométriques ainsi que la position du centre de l’inclusion sont normalisées
par la demi-largeur de contact a∗ (Tableau 4.2).
151
Géométrie Position
a1 = 0.4a∗ , a2 = a3 = 0.1a∗ dx3 = 0.4a∗
152
1.2 t=0
1
t=τ /2
0.8
P /P0
0.6
t=τ
0.4
0.2
SAM
FEM
0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x1 /a∗
153
4.4.3 Résultats
Les résultats présentés ici se limiteront à un cas académique. Le contact entre
une sphère rigide et un massif viscoélastique hétérogène est considéré. Les propriétés
viscoélastiques retenues pour le massif sont identiques à celles utilisées dans la par-
tie validation. Le but de cette partie sera d’étudier l’effet de la présence d’une ou
plusieurs hétérogénéités sur la solution du problème de contact viscoélastique. La
simulation dans le domaine temporel varie entre [0, τ ] et est divisée en 40 pas de
temps. a∗ (t) représente le rayon de contact instantané à un instant t donné corres-
pondant au cas homogène élastique c’est-à-dire sans hétérogénéité. Compte tenu du
fait que le matériau s’assouplit avec le temps, il est évident que a∗ (t) augmentera
avec le temps t. On peut introduire une variable γ qui est le ratio entre le mo-
dule de Young de l’inclusion élastique et le module de Young instantané du massif
viscoélastique (Eq. 4.44).
γ = E I /E M (0) (4.45)
La distribution de pression est présentée aux figures 4.19, 4.20, 4.21 et 4.22 pour
différentes valeurs du module γ = 1, 2, 0.5 et 0.2. Une hétérogénéité inclinée de demi-
axes a1 = 0.4a∗ , a2 = a3 = 0.1a∗ située à la profondeur dx3 = 0.4a∗ avec un angle
θ = 30◦ est considérée. Les résultats sont tracés dans le plan x2 = 0. Dans le premier
cas (Fig. 4.19) c’est-à-dire pour γ = 1 on observe que pour t = 0, l’hétérogénéité
n’a aucun effet sur la distribution de la pression de contact (ce qui est normal).
Cependant à t = τ /2 et t = τ la surpression due à l’hétérogénéité passe de 6% à 18%
comparé au cas homogène viscoélastique. Ce résultat reste cohérent compte tenu du
fait que la matrice est de moins en moins raide alors que l’hétérogénéité quant à elle
demeure élastique. Le cas d’une hétérogénéité dure (γ = 2) est maintenant considéré
Fig. 4.20. Ce cas entraine une surpression de 6%, 14% et 24% respectivement pour
t = 0, t = τ /2, t = τ . Pour γ = 0.5 (Fig. 4.21), on remarque que la surpression
induite par l’hétérogénéité élastique est négligeable pour t = τ /2. Alors qu’à t = 0
une réduction du maximum du champ de pression de 9.5% est obtenue ainsi qu’une
augmentation du maximum du champ de pression de 10.5% pour t = τ . Pour γ = 0.2
(Fig. 4.22) le champ de pression a diminué de 17.15% à t = 0 et de 12% pour t = τ /2.
Pour t = τ , l’effet de l’hétérogénéité élastique sur la distribution de pression est
négligeable.
154
1.2 1.2
t=0 t=0
1 1
t=τ t=τ
0.4 0.4
0.2 0.2
homogene homogene
heterogene heterogene
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x2 /a∗ x1 /a∗
(a) (b)
1.2 1.2
1 t=0 1 t=0
t=τ t=τ
0.4 0.4
0.2 0.2
homogene homogene
heterogene heterogene
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x2 /a∗ x1 /a∗
(a) (b)
155
1.2 1.2
1 t=0 1 t=0
t=τ t=τ
0.4 0.4
0.2 0.2
homogene homogene
heterogene heterogene
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x2 /a∗ x1 /a∗
(a) (b)
1.2 1.2
1 t=0 1 t=0
t=τ t=τ
0.4 0.4
0.2 0.2
homogene homogene
heterogene heterogene
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x2 /a∗ x1 /a∗
(a) (b)
156
1.2 1.2
1 t=0 1 t=0
t=τ t=τ
0.4 0.4
0.2 0.2
homogene homogene
heterogene heterogene
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x2 /a∗ x1 /a∗
(a) (b)
157
GPa, E3I = 50 GPa, µI12 = 83 GPa, µI13 = 400 GPa, µI23 = 20 GPa, ν12 I
= 0.15,
I I
ν13 = 0.26, ν23 = 0.40.
L’orientation du repère d’orthotropie par rapport au repère de contact est donnée
par la convention ZXZ des angles d’Euler : ϕ = 45◦ , θ = 60◦ , ψ = 30◦ . L’on peut re-
marquer qu’à l’instant initial, l’hétérogénéité se comporte comme une hétérogénéité
dure dans le plan x2 = 0 et comme une hétérogénéité souple dans le plan x1 = 0.
Ceci est principalement dû aux propriétés anisotropes et à l’orientation de l’inclu-
sion. Dans le plan x2 = 0 la surpression due à l’hétérogénéité varie de 2.1% (à t = 0)
à 49% (à t = τ ).
1.2 1.2
1 t=0 1 t=0
t=τ t=τ
0.4 0.4
0.2 0.2
homogene homogene
heterogene heterogene
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x2 /a∗ x1 /a∗
(a) (b)
158
1.2 1.2
1 t=0 1 t=0
t=τ t=τ
0.4 0.4
0.2 0.2
homogene homogene
heterogene heterogene
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x2 /a∗ x1 /a∗ t
(a) (b)
159
ξ
0.25 0.5
0
0.25
−0.25
−0.5
−0.25
−0.75
t=0 t=0
t = τ/2 t = τ/2
t=τ t=τ
−1 −0.5
0 60 120 180 240 300 360 0 60 120 180 240 300 360
ξ(degre) ξ(degre)
(a) (b)
160
(a) (b)
(c)
Figure 4.27 – Valeur adimensionnée de la contrainte de von Mises (σV M /P0 ) dans
le plan x2 = 0 en présence d’une inclusion unique isotrope sphérique (a1 = a2 =
a3 = 0.15a∗ , dx3 /a∗ = 0.4 ; (a) t = 0, (b) t = τ /2 et (c) t = τ .
161
162
Un contact sans frottement sera considéré dans les parties 4.6 et 4.7. Soit ω
la vitesse angulaire de la sphère dans le cas d’un mouvement de roulement pur
(Fig.4.28), et v la vitesse linéaire de glissement dans le cas d’un glissement pur.
Remarquons que :
v =R×ω (4.46)
v étant la vitesse de glissement et R le rayon de la sphère. Contrairement aux
autres approches existantes dans la littérature, le modèle présenté dans le cadre
de cette thèse peut prendre en compte des mouvements de roulement/glissement
à vitesse non constante. Une convention sur les repères sera introduite et restera
valable tout au long de ce chapitre. R0 (O0 , X10 , X20 , X30 ) correspond au repère du
contact à l’instant initial à t = 0. Pour t 6= 0 le repère courant du contact sera défini
par R(O, X1 , X2 , X3 ), O étant la projection du centre de la sphère sur la surface de
contact. R(O, X1 , X2 , X3 ) suit donc le mouvement. Les équations du contact seront
résolues dans le repère R(O, X1 , X2 , X3 ). L’incrément de déplacement ∆→ −
r de la
sphère à chaque pas de temps ∆t est de :
∆→
−
r = ∆t.→
−
v (4.47)
−→0 −→0 −→0
X1 , X2 , X3 sont respectivement les vecteurs unitaires suivant les axes (OX10 ),
(OX20 ), (OX30 ).
163
FN
v=ωR
ω
ܺଵ ≡ ܺଵ FN
O0 ܺଵ
Oo≡O v=ωR
ܺଶ ߙ. ∆ݎԦ
ω
ܺଶ ≡ܺଶ
ܺଷ ≡ܺଷ O
ܺଵ
ܺଶ
ܺଷ
ܺଷ
Z ∞ Z ∞
u3 (x1 , x2 , α∆t) = G(x1 − x01 , x2 − x02 )dx01 x02 ×
−∞ −∞
α
X
(1 − ν1 )J (1) [(α − k)∆t] + (1 − ν2 )J (2) [(α − k)∆t]
k=0
−→ 0 −→
×[p(x10 − α.∆→
−r .X10 , x20 − α.∆→
−
0
r .X20 , k) −
−→ 0 −→
p(x 0 − α.∆→−
r .X 0 , x 0 − α.∆→−
0
1 1r .X 0 , k − 1)]
2 2 (4.50)
164
Z ∞ Z ∞
u3 (x1 , x2 , α∆t) = G(x1 − x01 , x2 − x02 )dx01 x02 ×
−∞ −∞
α
X
((1 − ν1 )J (1) [(α − k)∆t] × [p(x01 , x02 , k) − p(x01 , x02 , k − 1)]
k=0
−→ 0 −→
+(1 − ν2 )J (2) [(α − k)∆t] × [p(x10 − α.∆→ −
r .X10 , x20 − α.∆→ −
0
r .X20 , k) −
−→ 0 −→
p(x 0 − α.∆→−r .X 0 , x 0 − α.∆→−
0
1 r .X 0 , k − 1)])
1 2 2 (4.51)
165
La méthode des rectangles peut être utilisée pour évaluer la valeur de l’intégrale
M (t). D’un point de vue macroscopique, les variables M (t) et µapp (t) /µv (t)
sont reliées entre elles par la relation suivante :
4.6.2 Validations
Le modèle semi-analytique développé sera validé par deux approches différentes.
Dans un premier temps, par une comparaison avec la méthode éléments finis en
utilisant le logiciel commercial Abaqus v6.11. Ensuite en utilisant l’approche basée
sur les éléments de frontières présentée par [CAR 13].
– Validation en utilisant la méthode éléments finis
La simulation éléments finis se fera entre [0, 2τ ]. Cet intervalle de temps est
divisé en 100 pas de temps fixes. Ainsi :
Le calcul a été réalisé en trois steps . Le step 1 consiste à faire une trans-
lation de la sphère rigide par rapport au centre du contact O d’une distance
égale à −v × τ . Dans le step 2, la charge normale est appliquée de manière
instantanée. Le matériau étant considéré comme élastique. Le dernier step,
i.e. le step 3 combine à la fois une translation de l’indenteur rigide à chaque
incrément de temps d’une distance ∆XF EM = v × ∆tF EM et la résolution du
contact viscoélastique.
– Validation en utilisant la méthode basée sur les éléments de frontière [CAR 13]
Carbone et al. [CAR 13], ont introduit récemment une nouvelle méthode de
prédiction du coefficient de frottement apparent dans le cadre du roulement
entre une sphère rigide et un matériau viscoélastique. Les résultats obtenus
ne sont valides qu’en régime permanent. Le cas transitoire n’étant pas pris en
compte. Une revue bibliographique étendue montrait, qu’il s’agissait de l’un
des rares travaux qui fournissait des résultats exacts en régime permanent. Le
début du mouvement de roulement/glissement (régime transitoire) n’étant pas
modélisé.
Afin de faciliter la comparaison entre les 3 différentes méthodes, les propriétés
matériaux, géométries et chargements considérés seront identiques à ceux utilisés
dans [CAR 13]. Ainsi pour la validation, une sphère rigide de rayon R = 10mm rou-
lant à vitesse constante sur un massif viscoélastique sera considérée. Les propriétés
matériaux du massif viscoélastique seront données par les équations 4.57 et 4.58.
Fonction de fluage :
1 1 t 1 1 t
J(t) = − exp(− ) = + (1 − exp(− )) (4.57)
µ0 µ1 τ µ∞ µ1 τ
166
Fonction de relaxation :
" #
t
R(t) = µ0 + (µ∞ − µ0 ) exp(− µ0 ) (4.58)
µ∞
τ
10
1
9
0.9
8
0.8
7
0.7
µ∞ × J (t)
6
R(t)/µ∞
0.6
0.5 5
0.4 4
0.3 3
0.2 2
0.1 1
0 0
0 0.5 1 1.5 2 0 1 2 3 4 5 6 7
t/τ t/τ
(a) (b)
167
0.2
Carbone,2013
FEM
SAM
0.15
P/P0
0.1
0.05
0
−1.5 −1 −0.5 0 0.5 1 1.5
x1 /a∗
0.4
Carbone,2013
FEM
SAM
0.3
P/P0
0.2
0.1
0
−2 −1 0 1 2
x1 /a∗
168
0.04
0.035
0.03
0.025
µapp
0.02
0.015
0.01
0.005 SAM
Carbone,2013
0
0 0.25 0.5 0.75 1 1.25 1.5 1.75 2
t/τ
de [CAR 13]. Les résultats de cette comparaison sont présentés à la figure 4.32
pour une valeur adimensionnée de la vitesse de glissement de vτ /a∗ = 0.6. Il
est quand même important de mentionner qu’ici le régime permanent est at-
teint pour t > τ /2. La méthode semi-analytique se présente alors comme une
méthode de détermination instantanée du coefficient de frottement apparent.
4.6.3 Résultats
[Link] Cas sphère rigide/massif viscoélastique
Dans ce cas de figure, il n’y a aucune différence entre glissement et roulement,
seule compte la vitesse de translation du centre de la sphère.
1. Chargement à vitesse constante
Il s’agit du cas le plus répandu dans la littérature. La vitesse de roule-
ment/glissement est supposée constante. Le déplacement de corps rigide sui-
vant l’axe de glissement/roulement X1 sera une rampe (cf. Fig. 4.33).
L’effet de la nature du chargement (normal) sur la réponse transitoire du
contact roulant/glissant sera étudié.
Chargement en déplacement normal imposé : Un chargement en
déplacement normal imposé est fixé à δ = 0.1mm. Les simulations seront
faites pour trois valeurs différentes de la vitesse adimensionnée vτ /a∗. Le cas
correspondant à t = 0 n’est pas présenté, car ce cas correspond à la solution
classique de Hertz. Les figures 4.34 et 4.35 présentent respectivement la distri-
bution de pression, le coefficient de frottement apparent et le moment résultant
pour différentes valeurs de t. Le régime permanent est ici atteint pour t = τ /2.
169
1.8
1.6
1.4
1.2
δx1 /a∗
1
0.4
0.2
0
0 0.5 1 1.5 2
t/τ
170
t = τ /10 t = τ /10
t = τ /5 t = τ /5
0.5 t = τ /4 0.5 t = τ /4
t=τ t=τ
t = 2τ t = 2τ
0.4 0.4
P/P0
P/P0
0.3 0.3
0.2 0.2
0.1 0.1
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x1 /a∗ x1 /a∗
(a) (b)
t = τ /10
t = τ /5
0.5 t = τ /4
t=τ
t = 2τ
0.4
P/P0
0.3
0.2
0.1
0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
x1 /a∗
(c)
171
0.05 0.1
0.045 0.09
0.04 0.08
0.035 0.07
0.03 0.06
µapp
M
0.025 0.05
0.02 0.04
0.015 0.03
0.01
∗
vτ /a =0.4 0.02 vτ /a∗ =0.4
0.005
vτ /a∗ =0.8 0.01
vτ /a∗ =0.8
vτ /a∗ =1.2 vτ /a∗ =1.2
0 0
0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2
t/τ t/τ
(a) (b)
1.5
1 t=0
0.5 sphere
t = τ /2
u/δ 0 t=τ
−0.5 t = 3τ /2
−1 t = 2τ
−1.5
−4 −3 −2 −1 0 1 2 3 4
x1 /a
Figure 4.36 – Déformation de la surface viscoélastique pour différents pas de temps
dans le cas d’un chargement en déplacement imposé de δ = 0.1mm, pour une vitesse
vτ /a∗ = 1.2.
172
0.6 0.6
t = τ /4 t = τ /4
t = τ /2 t = τ /2
0.5 t=τ 0.5 t=τ
t = 2τ t = 2τ
t = 3τ t = 3τ
t = 4τ t = 4τ
0.4 0.4
P/P0
P/P0
0.3 0.3
0.2 0.2
0.1 0.1
0 0
−4 −3 −2 −1 0 1 2 3 4 −4 −3 −2 −1 0 1 2 3 4
x1 /a∗ x1 /a∗
(a) (b)
0.6
t = τ /4
t = τ /2
0.5 t=τ
t = 2τ
t = 3τ
t = 4τ
0.4
P/P0
0.3
0.2
0.1
0
−4 −3 −2 −1 0 1 2 3 4
x1 /a∗
(c)
173
0.08 1.1
1
0.07
0.9
0.06 0.8
0.05 0.7
0.6
µapp
0.04
M
0.5
0.03 0.4
0.3
0.02
vτ /a∗ =0.4 0.2 vτ /a∗ =0.4
0.01 vτ /a∗ =0.8 vτ /a∗ =0.8
0.1
vτ /a∗ =1.2 vτ /a∗ =1.2
0 0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
t/τ t/τ
(a) (b)
δx1 = K × t2 (4.59)
v =2×K ×t (4.60)
K = accélération/2 (mm/s/s)
La position et la valeur du pic du coefficient de frottement est fonction de
l’accélération (Fig.4.40). Cette méthode est intéressante pour essayer d’esti-
mer la vitesse linéaire qui produit le frottement apparent maximal. Si K est
très petit la solution obtenue tend vers la solution en régime permanent. Les
courbes tendent donc vers une valeur asymptotique. Au lieu de faire plusieurs
simulations ou plusieurs essais pour différentes valeurs de vitesses, on peut faire
une seule simulation avec une très faible valeur de K permettant de couvrir
toute la gamme de vitesses utiles [CAR 13].
174
2 25
Kt2 /a∗ τ 2 (K = 2) 2Ktτ /a∗ (K = 2)
1.8 Kt2 /a∗ τ 2 (K = 1) 2Ktτ /a∗ (K = 1)
Kt2 /a∗ τ 2 (K = 1/2) 2Ktτ /a∗ (K = 1/2)
1.6 20
Kt2 /a∗ τ 2 (K = 1/4) 2Ktτ /a∗ (K = 1/4)
1.4 Kt2 /a∗ τ 2 (K = 1/8) 2Ktτ /a∗ (K = 1/8)
1.2 15
δx1 /a∗ τ 2 vτ /a∗
1
0.8 10
0.6
0.4 5
0.2
0 0
0 0.5 1 0 0.1 0.2 0.3 0.4 0.5
t/τ t/τ
(a) (b)
Figure 4.39 – (a) Déplacement de corps rigide suivant X1 (b) vitesse de glissement
suivant X1 dans le cas d’un roulement/glissement à accélération constante.
0.06
K=2
K=1
K=1/2
K=1/4
K=1/8
0.04
µapp
0.02
0
0 0.1 0.2 0.3 0.4 0.5
t/τ
175
3. Chargement cyclique
Des chargements cycliques seront considérés (Fig.4.41). Le coefficient de frot-
tement de Coulomb ne sera pas pris en compte. L’on se limitera unique-
ment à l’étude du coefficient de frottement apparent induit par la viscosité
du matériau, parfois appelé frottement visqueux. La vitesse, le déplacement
maximal et la fréquence de sollicitation sont reliés par la relation suivante :
δmax
vmax = (4.61)
4×λ
L’étude de l’effet de la fréquence λ peut se faire à déplacement tangentiel
constant δmax ou à vitesse constante vmax .
Ici l’étude sera faite à vitesse constante afin d’avoir la même valeur limite pour
le coefficient de frottement apparent. La figure 4.42 montre le coefficient de
frottement pour différentes fréquences de glissement pour un chargement en
force normale imposée FN = 1.48N . Pour λ = 1, la solution en régime perma-
nent est atteinte quasiment après le premier cycle. Ensuite plus la fréquence
augmente, plus il faut de cycles pour atteindre la solution en régime perma-
nent.
Dans le cas d’un chargement en déplacement imposé (Fig.4.43), il faut encore
moins de cycles (à fréquence identique) pour tendre vers la solution en régime
permanent. Ces constatations restent cohérentes avec les résultats présentés
plus haut dans la partie [Link].
Il est important de noter l’augmentation de la taille de la zone de contact ainsi
qu’une réduction du pic de pression avec la fréquence de sollicitation λ que
l’on soit en force imposée (Fig.4.44) ou en déplacement imposé (Fig.4.45). Les
distributions de pression ne sont présentées que pour l’instant t = 2τ .
Les variations de la zone de contact et des pics de pressions, en fonction de
la fréquence de sollicitations λ, sont résumées dans le tableau 4.5 pour un
chargement en force imposée et dans le tableau 4.6 pour un chargement en
déplacement imposé.
176
0.8
λ=1
λ=2
0.6 λ=4
λ=8
0.4
0.2
δx1 /a∗
0
−0.2
−0.4
−0.6
−0.8
0 0.5 1 1.5 2
t/τ
0.08 0.08
Regime Transitoire Regime Transitoire
0.06 Regime Permanent 0.06 Regime Permanent
0.04 0.04
0.02 0.02
µapp
µapp
0 0
−0.02 −0.02
−0.04 −0.04
−0.06 −0.06
−0.08 −0.08
0 0.5 1 1.5 2 0 0.5 1 1.5 2
t/τ t/τ
(a) (b)
0.08 0.08
Regime Transitoire Regime Transitoire
0.06 Regime Permanent 0.06 Regime Permanent
0.04 0.04
0.02 0.02
µapp
µapp
0 0
−0.02 −0.02
−0.04 −0.04
−0.06 −0.06
−0.08 −0.08
0 0.5 1 1.5 2 0 0.5 1 1.5 2
t/τ t/τ
(c) (d)
177
0.05 0.05
Regime Transitoire Regime Transitoire
0.04 Regime Permanent 0.04 Regime Permanent
0.03 0.03
0.02 0.02
0.01 0.01
µapp
µapp
0 0
−0.01 −0.01
−0.02 −0.02
−0.03 −0.03
−0.04 −0.04
−0.05 −0.05
0 0.5 1 1.5 2 0 0.5 1 1.5 2
t/τ t/τ
(a) (b)
0.05 0.05
Regime Transitoire Regime Transitoire
0.04 Regime Permanent 0.04 Regime Permanent
0.03 0.03
0.02 0.02
0.01 0.01
µapp
µapp
0 0
−0.01 −0.01
−0.02 −0.02
−0.03 −0.03
−0.04 −0.04
−0.05 −0.05
0 0.5 1 1.5 2 0 0.5 1 1.5 2
t/τ t/τ
(c) (d)
1.2
λ=1
λ=2
1 λ=4
λ=8
0.8
P /P0
0.6
0.4
0.2
−0.2
−2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5
x1 /a∗
178
1.2
λ=1
λ=2
1 λ=4
λ=8
0.8
P/P0
0.6
0.4
0.2
−0.2
−1.5 −1 −0.5 0 0.5 1 1.5
x1 /a∗
179
4. Chargement quelconque
Tous les cas présentés précédemment sont plus ou moins idéalisés, ces cas
sont plutôt rencontrés lors d’essais de laboratoires bien instrumentés. Dans
cette partie, il sera également montré que le modèle développé, est capable de
prendre en compte le cas d’un chargement quelconque (vitesse constante ou
accélération constante ou chargement cyclique). Dans les applications indus-
trielles nous sommes plutôt confrontés à des trajets de chargement quelconques
(Fig.4.46), qui peuvent certes dans certaines conditions se ramener à des cas
académiques (vitesse constante, accélération constante, . . . ). L’évolution du co-
efficient de frottement apparent sera présentée dans le cas de deux chargements
tangentiels quelconques. Une modélisation précise du roulement/glissement
viscoélastique évitera toute surestimation de µapp et donc de la force tangen-
tielle FV , et donc de mieux prédire la dégradation de ces types de matériaux.
0.06
Regime Transitoire
1.4
0.04 Regime Permanent
1.2
1 0.02
µapp
0.2 −0.04
0
−0.06
−0.2 0 0.5 1 1.5 2
0 0.5 1 1.5 2
t/τ t/τ
(a) (b)
Dans le cas où la sphère est rigide, et bien sûr lorsque le frottement de Cou-
lomb est négligé, il n’y a aucune différence entre un mouvement de glissement et
de roulement. Par contre, si les deux matériaux sont viscoélastiques, il n’y a plus
d’équivalence possible entre roulement et glissement.
180
181
0.05 0.06
τ1 = 0.25τ2
τ1 = 0.5τ2
τ1 = 2τ2
0.04 0.05 τ1 = 4τ2
τ1 = 1000τ2
τ1 → ∞
0.04
0.03
µapp µapp
0.03
0.02
0.02
τ1 = 0.25τ2
τ1 = 0.5τ2
0.01 τ1 = 2τ2 0.01
τ1 = 4τ2
τ1 = 1000τ2
τ1 → ∞
0 0
0 1 2 3 4 0 1 2 3 4 5 6 7 8
t/τ2 t/τ2
(a) (b)
182
0.04 0.06
Glissement τ1 = 0.25τ2 Glissement τ1 = 2τ2
Roulement τ1 = 0.25τ2 Roulement τ1 = 2τ2
0.05
0.03
0.04
µapp µapp
0.02 0.03
0.02
0.01
0.01
0 0
0 1 2 3 4 0 1 2 3 4
t/τ2 t/τ2
(a) (b)
0.06 0.06
Glissement τ1 = 4τ2 Glissement τ1 = 1000τ2
Roulement τ1 = 4τ2 Roulement τ1 = 1000τ2
0.05 0.05
0.04 0.04
µapp µapp
0.03 0.03
0.02 0.02
0.01 0.01
0 0
0 1 2 3 4 0 1 2 3 4
t/τ2 t/τ2
(c) (d)
183
0.05
0.045
0.04
0.035
µapp
0.03
0.025
0.02
0.015
Roulement vτ /a∗ =1.58
0.01 Roulement vτ /a∗ =0.79
0.005 Glissement vτ /a∗ =1.58
Glissement vτ /a∗ =0.79
0 −1 0 1 2 3 4
10 10 10 10 10 10
τ1 /τ2
Figure 4.49 – Valeur maximale du coefficient de frottement apparent, max µapp (t),
t∈[0,4τ2 ]
en fonction du ratio τ1 /τ2 , dans le cas d’un chargement en déplacement normal
imposé δ = 0.1mm.
184
γ = E I /E M (0) (4.64)
où E I est le module de Young de l’hétérogénéité élastique isotrope.
ν I = ν M = 0.3 (4.65)
185
FN FN
( )( ), ( )( ),
Oo≡O ≡ o
A A
≡ A A
( )( ), ( )( ),
≡
FN
( )( ),
=0 =2
O0 o
A A
( )( ),
=
Figure 4.50 – Vue 3D du mouvement de roulement/glissement d’une sphère sur un
massif viscoélastique hétérogène contenant une hétérogénéité ellipsoı̈dale, le mouve-
ment ayant lieu suivant l’axe X1 .
A-A
FN FN
A-A
( )( ), ( )( ),
≡
o
O ≡O O 0
o
dx
M Ɵ M
Ɵ
≡
2 2 A-A 2 2
( )( ), FN x3 ( )( ),
( )( ),
=0 =2
O0 o
M Ɵ
2 2
( )( ),
186
Données d’entrée :
• Propriétés géométriques
et matériaux
• Chargement
• Simulation pour ∈
[0, ]
• Roulement/glissement
• Profile du déplacement
Discrétisation temporelle en
pas de temps ∆ = ; ∈
0…….
=0
Détermination eigenstrain ε*
Contraintes en surface
Détermination de la surcontrainte
due à ε* (3D-FFT)
Déplacements ∗ et
= +1 Déplacements induits par la
contraintes induits par les
présence des hétérogénéités
contraintes en surfaces.
U ∗ , U ∗ (éventuellement)
Prise en
compte
Somme des contraintes puis des espace
déplacements
semi-infini
= ?
187
4.7.2 Résultats
188
0.8 0.8
P/P0
P/P0
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3
x1 /a∗ x01 /a∗
(a) (b)
0.04
0.035
0.03
0.025
µapp
0.02
0.015
0.01
0.005 Homogene
Heterogene
0
0 0.25 0.5 0.75 1 1.25 1.5 1.75 2
t/τ
189
190
0.8 0.8
P/P0
P/P0
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3
x1 /a∗ x01 /a∗
(a) (b)
0.8 0.8
P/P0
P/P0
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3
x1 /a∗ x01 /a∗
(c) (d)
0.8 0.8
P/P0
P/P0
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3
x1 /a∗ x01 /a∗
(e) (f)
0.8 0.8
P/P0
P/P0
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3
x1 /a∗ x01 /a∗
(g) (h)
191
0.04
0.035
0.03
0.025
µapp
0.02
0.015
0.01
Homogene
0.005 Heterogene θ = 300
Heterogene θ = 600
0
0 0.25 0.5 0.75 1 1.25 1.5 1.75 2
t/τ
Figure 4.56 – Coefficient de frottement apparent pour différentes orientations an-
gulaires θ de l’hétérogénéité ellipsoı̈dale γ = E I /E M (0) = 1 (a1 = 0.4a∗ , a2 =
0.5a∗ , a3 = 0.2a∗ , dx3 = 0.4a∗ ) ; pour un déplacement normal imposé de δ = 0.1mm
et pour une vitesse vτ /a∗ = 0.6.
l’approche.
192
0.8 0.8
P/P0
P/P0
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3
x1 /a∗ x01 /a∗
(a) (b)
0.8 0.8
P/P0
P/P0
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3
x1 /a∗ x01 /a∗
(c) (d)
0.8 0.8
P/P0
P/P0
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3
x1 /a∗ x01 /a∗
(e) (f)
193
0.04
0.035
0.03
0.025
µapp
0.02
0.015
0.01 Homogene
γ = 0.1(E I = E 0 )
0.005 γ = 1(E I = E ∞ )
γ=4
0
0 0.25 0.5 0.75 1 1.25 1.5 1.75 2
t/τ
Figure 4.58 – Coefficient de frottement apparent pour différentes valeurs de γ quand
le massif viscoélastique contient une hétérogénéité ellipsoı̈dale (a1 = 0.4a∗ , a2 =
0.5a∗ , a3 = 0.2a∗ , dx3 = 0.4a∗ , θ = 30◦ ) ; pour un déplacement normal imposé de
δ = 0.1mm et pour une vitesse vτ /a∗ = 0.6.
dans le cas des hétérogénéités ellipsoı̈dales cet effet est beaucoup plus marqué en
régime transitoire. Il faut toutefois noter que certaines distributions d’hétérogénéités
sphériques peuvent avoir des effets beaucoup plus importants sur le coefficient de
frottement apparent.
Il est possible d’architecturer des matériaux viscoélastiques hétérogènes de
manière à augmenter ou à réduire le coefficient de frottement apparent. Les mêmes
simulations peuvent également être réalisées sur un matériau composite réel.
194
0.04 120
γ = 4, θ = 80◦
100 γ = 4, θ = −80◦
0.035
80
γ = 4, θ = 45◦
γ = 4, θ = −45◦
0.03 = 0.01, θ = 80◦
hom
40 γ = 0.01, θ = 45◦
0.025
20 γ = 0.01, θ = −45◦
µapp
hom
0.02 0
γ = 4, θ = 80◦
γ = 4, θ = −80◦ −20
0.015 γ = 4, θ = 45◦
(µheter
−40
γ = 4, θ = −45◦
0.01 Homogene −60
γ = 0.01, θ = 80◦
γ = 0.01, θ = −80◦ −80
0.005
γ = 0.01, θ = 45◦ −100
γ = 0.01, θ = −45◦
0 −120
0 0.25 0.5 0.75 0 0.25 0.5 0.75
t/τ t/τ
(a) (b)
195
0.04
20
0.035
15
app − µapp )/µapp (%)
0.03 10
hom
0.025 5
µapp
hom
0.02 0
0.015 −5
(µheter
0.01 −10
γ=4
0.005 Homogene −15 γ=4
γ = 0.01 γ = 0.01
0 −20
0 0.25 0.5 0.75 0 0.25 0.5 0.75
t/τ t/τ
(a) (b)
196
(a)
(b)
(c)
(d)
(e)
Figure 4.61 – Contrainte de von Mises pour différents pas de temps dans le cas
d’une hétérogénéité sphérique de rayon r = 0.15a∗ située à dxI3 = 0.3a∗ (γ = 0.4) ;
pour un déplacement normal imposé de δ = 0.1mm et pour une vitesse vτ /a∗ = 0.6 ;
(a) t = 0, (b) t = τ /2, (c) t = τ , (d) t = 3τ /2, (e) t = 2τ .
197
198
Contact Normal
Contact Tangentiel
∆sτ (x1 , x2 , t)
qτ (x1 , x2 , t) = −µcoul · p (x1 , x2 , t) · ∀ (x1 , x2 ) ∈ Γsl (t) (4.69)
k ∆sτ (x1 , x2 , t) k
tot
∆uτ (x1 , x2 , t) − ∆δτ (t) = ∆sτ (x1 , x2 , t) ∀ (x1 , x2 ) ∈ Γsl (t) (4.70)
k qτ (x1 , x2 , t) k< µcoul · p (x1 , x2 , t) ∀ (x1 , x2 ) ∈ Γst (t) (4.71)
∆utot
τ (x1 , x2 , t) − ∆δτ (t) = 0 ∀ (x1 , x2 ) ∈ Γst (t) (4.72)
X
q (x1 , x2 , t) S = Qτ (t) (4.73)
Γp
199
Les champs de déplacements présentés dans les équations 4.67 et 4.74 sont donnés
par :
Z ∞ Z ∞ Z t
q (1) 0 0
uq1,2,3
x1
(x1 , x2 , t) = x1
(T1,2,3 (x1 − x1 , x2 − x2 , t − ξ) +
−∞ −∞ 0
q (2) 0 0 ∂qx1 (x01 , x02 , ξ) 0 0
x1
T1,2,3 (x1 − x1 , x2 − x2 , t − ξ)) dx1 x2 dξ (4.76)
∂ξ
Z ∞ Z ∞ Z t
qx2 q (1) 0 0
u1,2,3 (x1 , x2 , t) = x2
(T1,2,3 (x1 − x1 , x2 − x2 , t − ξ) +
−∞ −∞ 0
q (2) 0 0 ∂qx2 (x01 , x02 , ξ) 0 0
x2
T1,2,3 (x1 − x1 , x2 − x2 , t − ξ)) dx1 x2 dξ (4.77)
∂ξ
Z ∞ Z ∞ Z t
p(1) 0 0
up1,2 (x1 , x2 , t) = (T1,2 (x1 − x1 , x2 − x2 , t − ξ) +
−∞ −∞ 0
p(2) 0 ∂p(x01 , x02 , ξ) 0 0
0
T1,2 (x1 − x1 , x2 − x2 , t − ξ)) dx1 x2 dξ (4.78)
∂ξ
qx1 qx2
Les tenseurs T1,2,3 , T1,2,3 contiennent des propriétés spatiales et matériaux (la
fonction de fluage des matériaux viscoélastiques en contact). Ces tenseurs sont
présentés dans l’annexe A.
Une approche numérique identique à celle présentée en 4.3.2 sera utilisée pour
la résolution du problème.
4.9.2 Résultats
Le couplage normal/tangentiel ne présente d’intérêt qu’en présence de coefficient
de frottement de Coulomb µcoul . Ici, un problème de contact avec un indenteur
sphérique et un solide élastique semi-infini viscoélastique homogène sera considéré.
Le coefficient de frottement de Coulomb est fixé à µcoul = 0.7. L’indenteur de rayon
R = 10mm sera supposé rigide.
Les fonctions de fluage et de relaxation du massif semi-infini viscoélastique sont
données par :
Fonction de fluage :
1 1 t 1 1 t
J(t) = − exp(− ) = + (1 − exp(− )) (4.79)
µ0 µ1 τ µ∞ µ1 τ
200
0.8
0.6
0.4
0.2
δx1 /a∗
0
−0.2
−0.4
−0.6
−0.8
0 0.5 1 1.5 2
t/τ
δx1
Figure 4.62 – Chargement tangentiel a∗
∈ [−0.62, +0.62].
Fonction de relaxation :
" #
t
R(t) = µ0 + (µ∞ − µ0 ) exp(− µ0 ) (4.80)
µ∞
τ
201
0.6 1.6
1.4
0.5
1.2
0.4
1
δ(mm) FN (N )
0.3 0.8
0.6
0.2
0.4
0.1
0.2
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
t/τ t/τ
(a) (b)
202
1.2 1.2
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−0.2 −0.2
−0.4 −0.4
−0.6 −0.6
−0.8 −0.8
P/P0 P/P0
−1 −1
Qx1 /µcoul × P0 Qx1 /µcoul × P0
−1.2 −1.2
−1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5
x1 /a∗ x1 /a∗
(a) (b)
1.2 1.2
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−0.2 −0.2
−0.4 −0.4
−0.6 −0.6
−0.8 −0.8
P/P0 P/P0
−1 −1
Qx1 /µcoul × P0 Qx1 /µcoul × P0
−1.2 −1.2
−1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5
x1 /a∗ x1 /a∗
(c) (d)
1.2 1.2
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−0.2 −0.2
−0.4 −0.4
−0.6 −0.6
−0.8 P/P0 −0.8
Qx1 /µcoul × P0 P/P0
−1 −1
−P/P0 Qx1 /µcoul × P0
−1.2 −1.2
−1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5
x1 /a∗ x1 /a∗
(e) (f)
203
1.2 1.2
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−0.2 −0.2
−0.4 −0.4
−0.6 −0.6
−0.8 −0.8
P/P0 P/P0
−1 −1
Qx1 /µcoul × P0 Qx1 /µcoul × P0
−1.2 −1.2
−2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5
x1 /a∗ x1 /a∗
(a) (b)
1.2 1.2
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−0.2 −0.2
−0.4 −0.4
−0.6 −0.6
−0.8 −0.8
P/P0 P/P0
−1 −1
Qx1 /µcoul × P0 Qx1 /µcoul × P0
−1.2 −1.2
−2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5
x1 /a∗ x1 /a∗
(c) (d)
1.2 1.2
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−0.2 −0.2
−0.4 −0.4
−0.6 −0.6
−0.8 P/P0 −0.8
Qx1 /µcoul × P0 P/P0
−1 −1
−P/P0 Qx1 /µcoul × P0
−1.2 −1.2
−2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5
x1 /a∗ x1 /a∗
(e) (f)
204
0.15
δx1 /a∗ ∈ [-0.08,0.08]; λ/τ = 8
0.1
0.05
δx1 /a∗
0
−0.05
−0.1
−0.15
0 0.5 1 1.5 2 2.5 3 3.5 4
t/τ
δx1
Figure 4.66 – Chargement en déplacement tangentiel imposé a∗
∈ [−0.08, +0.08],
λ/τ = 8.
Le chargement en déplacement tangentiel imposé est tel que δax1∗ ∈ [−0.08, +0.08].
32 cycles de fretting ont été réalisés entre [0, 4τ ], ce qui fait un ratio λ/τ = 8,
Fig.4.66. Les simulations sont réalisées pour différentes valeurs du coefficient de
frottement de Coulomb µcoul = 0.4, 0.7, 1.0.
Lorsque µcoul = 0.4 (Fig.4.67b), il n’y a pas une grande différence entre le cas
élastique et le cas viscoélastique. Le contact se retrouve très vite en glissement total.
Le coefficient de frottement apparent s’obtient facilement en faisant la somme µapp =
µcoul + µv . Pour µcoul = 0.7 (Fig.4.67d), il est important de noter que le cas élastique
n’est pas en stick-slip alors que dans le cas viscoélastique, le contact est glissement
total. L’énergie dissipée par fretting est quasiment multipliée par 3. Pour µcoul = 1.0
(Fig.4.67f), la boucle élastique reste quasiment inchangée comparé au cas µcoul = 0.7.
Par contre, le cas viscoélastique atteint toujours le régime de glissement total. A
noter que le régime de glissement total, dans le cas viscoélastique, ne disparait que
pour µcoul = 1.5 (Fig.4.68b). Il est important de remarquer que µv ne varie quasiment
pas avec le coefficient de frottement de Coulomb (Figs. 4.67a,4.67c,4.67e,4.68a).
205
0.08 1
Massif Elastique
0.8 Massif viscoElastique Qx1 /FN
0.06 Massif viscoElastique (Qx1 + FV )/FN
0.6
0.04
0.4
µv = FV /FN
0.02
0.2
Q/FN
0 0
−0.2
−0.02
−0.4
−0.04
−0.6
−0.06
−0.8 µcoul=0.4
−0.08 −1
−0.1 −0.05 0 0.05 0.1 −0.15 −0.1 −0.05 0 0.05 0.1 0.15
δx1 /a∗ δx1 /a∗
(a) (b)
0.08 1.2
Massif Elastique
1 Massif viscoElastique Qx1 /FN
0.06 Massif viscoElastique (Qx1 + FV )/FN
0.8
0.04 0.6
0.4
µv = FV /FN
0.02
0.2
Q/FN
0 0
−0.2
−0.02
−0.4
−0.04 −0.6
−0.8
−0.06
−1 µcoul=0.7
−0.08 −1.2
−0.1 −0.05 0 0.05 0.1 −0.15 −0.1 −0.05 0 0.05 0.1 0.15
δx1 /a∗ δx1 /a∗
(c) (d)
0.08 1.6
Massif Elastique
1.4 Massif viscoElastique Qx1 /FN
0.06 1.2 Massif viscoElastique (Qx1 + FV )/FN
1
0.04
0.8
µv = FV /FN
0.02 0.6
0.4
0 0.2
Q/FN
0
−0.02 −0.2
−0.4
−0.04
−0.6
−0.06 −0.8
−1 µcoul=1.0
−0.08 −1.2
−0.1 −0.05 0 0.05 0.1 −0.15 −0.1 −0.05 0 0.05 0.1 0.15
δx1 /a∗ δx1 /a∗
(e) (f)
206
0.08 2.4
Massif Elastique
2 Massif viscoElastique Qx1 /FN
0.06 Massif viscoElastique (Qx1 + FV )/FN
1.6
0.04 1.2
0.8
µv = FV /FN
0.02
0.4
Q/FN
0 0
−0.4
−0.02
−0.8
−0.04 −1.2
−1.6
−0.06
−2 µ =1.5
coul
−0.08 −2.4
−0.1 −0.05 0 0.05 0.1 −0.15 −0.1 −0.05 0 0.05 0.1 0.15
δx1 /a∗ δx1 /a∗
(a) (b)
Figure 4.68 – Boucles de fretting de µv pour µcoul = 1.5 (a) / Boucles de fretting
de µapp pour µcoul = 1.5 (b).
La valeur du coefficient de frottement est toujours fixée à µcoul = 0.7. Cette fois-
ci la fréquence du déplacement tangentiel est identique pour les deux simulations
λ/τ = 8. Par contre, l’amplitude varie de la sollicitation varie : δax1∗ ∈ [−0.32, +0.32]
pour la première simulation et δax1∗ ∈ [−0.08, +0.08] pour la seconde, Fig.4.71.
Lorsque l’amplitude est élevée, le contact se retrouve très vite en glissement total
(Figs. 4.72c et 4.72d), même dans le cas élastique. Les boucles obtenues dans le cas
viscoélastique deviennent très proches de celles obtenues dans le cas élastique. Dans
ce cas de figure, le coefficient de frottement µv , induit par le seul fait de la dissymétrie
des champs de pression, est plus élevé. Ici µv est proche de 0.1 (Figs. 4.72a et 4.72b),
ce qui représente près de 14% du coefficient de frottement de Coulomb. Ce qui peut
207
0.05
δx1 /a∗
0
−0.05
−0.1
208
0.08 0.08
0.06 0.06
0.04 0.04
µv = FV /FN
µv = FV /FN
0.02 0.02
0 0
−0.02 −0.02
−0.04 −0.04
−0.06 −0.06
−0.08 −0.08
−0.1 −0.05 0 0.05 0.1 −0.1 −0.05 0 0.05 0.1
δx1 /a∗ δx1 /a∗
(a) (b)
1.2 1.2
Massif Elastique Massif Elastique
1 Massif viscoElastique Qx1 /FN 1 Massif viscoElastique Qx1 /FN
Massif viscoElastique (Qx1 + FV )/FN Massif viscoElastique (Qx1 + FV )/FN
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
Q/FN
Q/FN
0 0
−0.2 −0.2
−0.4 −0.4
−0.6 −0.6
−0.8 −0.8
−1 µcoul=0.7 −1 µcoul=0.7
−1.2 −1.2
−0.15 −0.1 −0.05 0 0.05 0.1 0.15 −0.15 −0.1 −0.05 0 0.05 0.1 0.15
δx1 /a ∗
δx1 /a∗
(c) (d)
Figure 4.70 – Boucles de fretting de µv pour différentes fréquences (a) λ/τ = 1 (b)
λ/τ = 4 / Boucles de fretting de µapp pour différentes fréquences (c) λ/τ = 1 (d)
λ/τ = 4.
209
0.2
δx1 /a∗
0
−0.2
−0.4
210
0.12 0.08
0.1
0.06
0.08
0.06 0.04
µv = FV /FN
0.04
µv = FV /FN
0.02
0.02
0 0
−0.02
−0.02
−0.04
−0.06 −0.04
−0.08
−0.06
−0.1
−0.12 −0.08
−0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 −0.1 −0.05 0 0.05 0.1
∗
δx1 /a δx1 /a∗
(a) (b)
1.4 1.2
Massif Elastique Massif Elastique
1.2 Massif viscoElastique Qx1 /FN 1 Massif viscoElastique Qx1 /FN
Massif viscoElastique (Qx1 + FV )/FN Massif viscoElastique (Qx1 + FV )/FN
1 0.8
0.8 0.6
0.6
0.4
0.4
0.2
Q/FN
0.2
Q/FN
0
0
−0.2
−0.2
−0.4
−0.4
−0.6 −0.6
−0.8 −0.8
−1 µcoul=0.7 −1 µcoul=0.7
−1.2 −1.2
−0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 −0.15 −0.1 −0.05 0 0.05 0.1 0.15
δx1 /a ∗
δx1 /a∗
(c) (d)
211
4.11 Synthèse
Dans ce chapitre, la résolution du contact entre matériaux viscoélastiques
homogènes/hétérogènes a été présentée. A noter que le contact viscoélastique
est très peu étudié dans la littérature. Dans le modèle, les hétérogénéités sont
élastiques (isotropes/anisotropes) et seul l’un des corps en contact peut contenir
des hétérogénéités. Le chapitre est subdivisé en trois parties.
La première partie a été consacrée à un problème d’indentation entre matériaux
viscoélastiques homogènes et hétérogènes. Dans le cas de l’indentation, seul le
problème de contact normal est résolu. L’approche développée a été validée avec un
modèle éléments-finis. Plusieurs types de chargement peuvent être pris en compte.
Des cas d’applications à l’indentation du PMMA (PolyMethylMethAcrylate) ont été
réalisés. Un déphasage entre la réponse du contact et le chargement imposé peut être
observé lorsqu’on est présence d’un chargement sinusoı̈dal. Dans le cas du contact
viscoélastique hétérogène, les résultats obtenus au chapitre 3 dans le cas de la matrice
élastique ne sont plus valides. L’effet de l’hétérogénéité sur le problème de contact
varie avec le temps et nécessite donc une modélisation fine. Le modèle peut prendre
en compte une ou plusieurs hétérogénéités et peut traiter le cas des matériaux com-
posites constitués d’une matrice viscoélastique et des renforts élastiques.
La deuxième partie s’est intéressée à des problématiques de roulement/glissement
entre deux corps viscoélastiques ou entre un corps rigide et un massif viscoélastique
[CAR 13]. L’approche a été validée avec une solution éléments finis et la solution
de référence de Carbone [CAR 13] dans le cas des massifs homogènes. Le roule-
ment/glissement entre massif viscoélastique hétérogène a été également étudié. Il
a été montré par exemple qu’une hétérogénéité unique sphérique ne modifie pas le
coefficient de frottement apparent. Le modèle, contrairement à la plupart des autres
méthodes de la littérature, permet de décrire l’évolution des contraintes en sous-
couche et notamment au voisinage de l’hétérogénéité.
212
La dernière partie est consacrée à l’étude du fretting dans le cadre des matériaux
viscoélastiques. Le contact tangentiel est résolu et le couplage problème normal
et tangentiel a été pris en compte. L’analyse s’est limitée au cas purement ho-
mogène pour des raisons de précision et de concision. La transition stick-slip varie
considérablement selon que l’on soit en force normale imposée ou en déplacement
normal imposé. L’analyse des boucles de fretting montre que le coefficient de frot-
tement apparent augmente considérablement. L’énergie dissipée (aire sous la boucle
de fretting) peut être 2 à 5 fois supérieure à celle dissipée dans le cas d’un matériau
élastique et varie grandement en fonction de l’amplitude et la fréquence du charge-
ment. Il est dès lors possible de choisir ou de concevoir au mieux des revêtements
viscoélastiques en fonction des applications considérées réduisant ainsi le nombre
d’essais.
L’étude réalisée dans ce chapitre a permis d’avoir une compréhension approfondie du
contact viscoélastique, car il s’agit d’un problème très peu étudié dans la littérature.
Le modèle peut prendre en compte des matériaux viscoélastiques avec plusieurs
hétérogénéités et donc un matériau composite constitué d’une matrice viscoélastique
et des renforts élastiques. L’avantage ici est qu’il est possible de résoudre un problème
de contact en prenant en compte toute la microstructure (simulation dite full field),
tout en gardant des temps de calculs raisonnables. Jusqu’à présent, l’approche
développée a été validée de manière numérique. Il convient maintenant de propo-
ser un dispositif expérimental permettant une validation expérimentale de notre
approche.
213
214
Sommaire
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 217
5.2 Dispositif expérimental . . . . . . . . . . . . . . . . . . . . . 217
5.2.1 Cahier de charges . . . . . . . . . . . . . . . . . . . . . . . . . 217
5.2.2 Principe de fonctionnement . . . . . . . . . . . . . . . . . . . 218
5.3 Eprouvette . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223
5.3.1 Eprouvette homogène . . . . . . . . . . . . . . . . . . . . . . 224
5.3.2 Eprouvette hétérogène . . . . . . . . . . . . . . . . . . . . . . 225
5.3.3 Eprouvette composites . . . . . . . . . . . . . . . . . . . . . . 225
5.4 Réalisation des essais . . . . . . . . . . . . . . . . . . . . . . . 226
5.4.1 Propriétés matériaux des corps en contact . . . . . . . . . . . 226
215
216
5.1 Introduction
Un dispositif expérimental développé au LAMCOS et considérablement modifié
tout au long de cette thèse, a été utilisé pour réaliser des essais de fretting. L’origina-
lité de l’approche réside dans le fait qu’on peut mesurer les champs de déplacements
à l’interface des surfaces en contact. Le dispositif a mis beaucoup de temps à se
mettre en place suite à de très nombreux essais infructueux.
Pour les matériaux homogènes, les solutions analytiques ainsi que les essais clas-
siques prédisent assez bien ce qui se passe à l’interface des matériaux en contact.
Par contre, avec l’avènement des matériaux composites, notamment dans le cas du
contact aube/disque, il devient indispensable de connaı̂tre finement ce qui se passe
à l’interface afin de valider les modèles numériques, de développer des lois d’endom-
magements phénoménologiques, . . . La plupart des tribomètres existants se limitent
à la mesure du coefficient de frottement de Coulomb et ne permettent pas de mesu-
rer les champs de déplacements ou de déformations réelles à l’interface des corps en
contact pour les trois modes de fretting. Le but de notre dispositif est donc de four-
nir une solution économique, simple et efficace permettant de répondre à ce besoin.
Le dispositif expérimental a fait l’objet d’un dépôt de brevet SNECMA/LAMCOS.
Le chapitre se structurera en quatre grandes parties. Dans un premier temps la ver-
sion finale fonctionnelle du dispositif expérimental sera présentée. Ensuite un zoom
sera fait sur les différents types d’éprouvettes d’essais utilisées. Le protocole d’essais
sera détaillé dans la troisième partie. Enfin, quelques résultats d’essais sur diverses
éprouvettes seront présentés. Une comparaison essai/calcul sera effectuée dans le cas
de l’éprouvette homogène.
217
Le but de cette démarche expérimentale est donc de réussir à obtenir les champs de
déplacements et de déformations lors des essais des trois modes de fretting, tels que
définis par [MOH 95] (Fig.5.1).
218
Le dispositif est composé d’un bâti rigide permettant d’éviter des chargements
parasites et d’avoir des mesures assez précises. En effet, les déplacements à mesurer
étant de l’ordre du micron, la moindre déformation du tribomètre pourrait entrainer
d’importantes erreurs de mesure. Une photo du dispositif expérimental (tribomètre+
caméra) ainsi que tous les composants électroniques qui le composent (contrôleur,
. . . ) sont présentés à la figure 5.2a. La figure 5.2b présente la configuration CAO.
Le tribomètre comprend plusieurs composants
– un bâti rigide.
– une plaque de verre ayant de très bonne qualité optique. Les caractéristiques
de la vitre utilisée sont données dans l’annexe B.
– une éprouvette de forme cylindrique dont le bout est surmonté d’une ca-
lotte sphérique. L’éprouvette pouvant être homogène, hétérogène ou composite
([Link] 5.3). L’éprouvette doit présenter un mouchetis naturel pour
pouvoir utiliser la corrélation d’images.
– Un capteur dénommé MC3A (de la société PM Instrumentation) mesurant
trois composantes de forces FX , FY , FZ , et trois composantes de moments
MX , MY , MZ .
– L’éprouvette est donc montée sur le capteur, l’ensemble (capteur, éprouvette)
est mis en contact avec un verre rigide (Fig.5.3) ayant de très bonne qualité
optique.
– Trois actionneurs permettent de reproduire les trois modes de fretting. Le
premier actionneur est un moteur pas à pas configuré pour appliquer une
force dans une direction verticale (Z), perpendiculaire à l’interface de contact.
Le second actionneur est également un moteur pas à pas électrique configuré
pour appliquer un déplacement alternatif suivant une direction parallèle Y , co-
linéaire à la surface de contact. Un troisième actionneur constitué d’un moteur
pas à pas permettant d’appliquer un moment autour de l’axe Z, MZ .
– Un système asservi permettant de contrôler les différents types de chargements
appliqués.
– La surface ou interface de contact est éclairée à l’aide d’une lampe
pieuvre de part et d’autre de l’échantillon (cf. Fig 5.4). Il s’agit d’un éclairage
à incidence rasante de manière à éviter toute surbrillance de la surface à
étudier. Il est possible d’envisager l’utilisation d’une lampe annulaire à LED
en forme de parabole, de type CCS LDR2 afin d’améliorer l’homogénéité de
l’éclairage sur l’interface de contact.
– Deux capteurs LVDT sont utilisés pour mesurer le déplacement de corps ri-
gide suivant Z. Le premier se trouve avant le capteur de force, le second au
niveau de l’échantillon. Le dispositif peut être aussi utilisé pour faire des essais
d0 indentation, même si ce n’est pas sa vocation première.
– Une caméra CCD (de 29 Mégapixels) et des moyens informatiques configurés
pour déterminer les champs de déplacements sur la surface de contact et
ainsi remonter aux champs de déformations. La caméra haute résolution (cf.
Fig.5.5) est dotée d’un objectif télé-centrique (caractérisé par un grandisse-
219
(a)
(b)
220
(a)
(b)
221
ment constant quel que soit l’éloignement de l’objet par rapport à l’objectif)
est placée en dessous de la vitre transparente.
– A partir des images enregistrées par la caméra, il est possible de remonter aux
champs de déplacements et donc de déformations en utilisant une technique
de corrélation d’images numériques.
En utilisant ce tribomètre, les conditions de contact ne sont pas perturbées lors
de la mesure, ce qui très difficile à réaliser. En résumé, le dispositif expérimental
décrit plus haut permet :
– de soumettre l’interface de contact (éprouvette/vitre) à des sollicitations de
fretting en mode I, II ou III.
– d’identifier et de suivre les positions des points en contact au fur et à mesure
du chargement.
222
– de faire une capture des images à l’aide d’une caméra et des logiciels d’acqui-
sition d’images.
– d’analyser les images enregistrées en utilisant des techniques de corrélations
d’images de façon à déterminer les champs de déplacements de ladite interface.
Le principe tel que présenté paraı̂t plus ou moins simple. Par contre, plusieurs
difficultés restent à surmonter avant d’aboutir à un dispositif fonctionnel.
– Il était par exemple très difficile de faire coı̈ncider le point d’application de la
force et le centre de la zone de contact.
– d’éviter de perturber la surface de contact.
– Le capteur multi-composants à jauge utilisé initialement se déformait
énormément (cf Fig. 5.6). Il en résultait donc un important déplacement de
corps rigide suivant Y lors des sollicitations en mode I. Il fallait donc avoir
un capteur le plus rigide possible. Ce principe est contraire au caractère in-
trinsèque même des capteurs à jauge, dont les mesures sont basées sur une
déformation d’un corps d’épreuve. L’utilisation des capteurs piézoélectriques
éviterait ce type de problème. Cependant les capteurs piézoélectriques se
révèleront, après analyse, trop lourds. De plus, ces types de capteurs réalisent
une une mesure relative et non une mesure absolue . Finalement, un
capteur à jauge MC3A (cf. Annexe B) à très haute performance, avec de très
faibles déflections, sera retenu.
5.3 Eprouvette
La géométrie de l’éprouvette peut être décomposée en deux parties : un corps
cylindrique et une extrémité constituée d’une calotte sphérique (Fig.5.7a).
La détermination des champs de déplacements de la surface de contact est basée
sur le principe de la corrélation d’images. Cela nécessite la présence d’un mouchetis
sur la surface de contact. Deux possibilités s’offrent à nous en ce qui concerne la
conception et la fabrication des éprouvettes :
223
(a) (b)
mesure.
– Concevoir des éprouvettes dont la microstructure peut servir de mouchetis. En
jouant sur la microstructure, il est alors possible de diminuer ou d’augmenter le
niveau de gris. C’est cette deuxième option qui a été retenue pour la réalisation
des essais (Fig.5.7b).
Le processus de fabrication des éprouvettes homogènes et hétérogènes sera présenté
dans un premier temps.
224
225
(a) (b)
Figure 5.9 – Tomographie d’une éprouvette hétérogène contenant (a) une ou (b)
deux billes céramiques.
226
227
où u(x) décrit le champ de déplacement d’une image déformée à partir d’une
image de référence. Deux approches différentes seront mises en œuvre pour le post-
traitement des résultats par corrélation d’images.
228
(a) (b)
229
les éléments de référence (éléments finis) sont des quadrilatères. La méthode globale
sera privilégiée pour le post-traitement des essais sur éprouvettes composites. Les
éprouvettes composites présentent de fortes variations de microstructures et donc
de fortes variations de mouchetis. La méthode globale sera donc la seule à même de
donner des résultats fiables.
(a) (b)
230
231
232
la courbe d’indentation, n’est pas très bonne. Ceci est lié principalement à la
résolution du capteur LVDT utilisé pour la mesure du déplacement normal. La
précision du capteur est de l’ordre de 10 microns. Il serait possible d’améliorer
considérablement les résultats d’indentation en choisissant un capteur avec une
meilleure précision, et surtout en utilisant une table anti-vibration.
2100
1800
1500
1200
FN (N )
900
600
300
Experimental
Code SA
0
0 50 100 150 200 250
δ(µm)
233
234
1.5
FN /F0
1
0.5
0
0 50 100
δ/δ0
235
La figure 5.20 présente une distribution des champs de pression dans le cas de
l’éprouvette homogène. La gamme de force est limitée entre [0, 1000N ]. Il est possible
de mesurer la demi-largeur de contact et tracer la courbe a∗ = ff onc (F ), a∗ étant la
demi-largeur de contact. Une bonne corrélation essai-calcul est obtenue (Fig.5.21).
La même démarche peut être appliquée aux éprouvettes composites.
10
7
a∗ (mm)
1 Experimental
Analytique/Numerique
0
0 200 400 600 800 1000
FN (N )
236
5.7 Synthèse
Dans ce chapitre, un dispositif expérimental original a été mis en place afin de
valider les méthodes numériques développées dans le cas du contact entre matériaux
hétérogènes ou composites. Le dispositif expérimental a fait l’objet d’un brevet
SNECMA/LAMCOS.
Les différentes composantes du dispositif ont été présentées dans un premier temps.
Quelques difficultés liées la mise en place du dispositif ont été énumérées succincte-
ment.
Dans la deuxième partie, diverses éprouvettes d’essais ont été présentées :
éprouvettes homogènes, hétérogènes et composites. Un accent particulier a été
mis sur le processus de fabrication des éprouvettes polymères. Des compositions
spécifiques de résines permettent en effet d’avoir un très bon contraste pour la
corrélation d’images. Dans le cas des éprouvettes polymères hétérogènes, des ana-
lyses par tomographie ont été réalisées afin de localiser la position des hétérogénéités
par rapport à la surface de contact.
Le protocole expérimental a ensuite été présenté dans le cas des trois modes de fret-
ting. Plusieurs tests ont été réalisés afin de s’assurer de la reproductibilité des essais.
Une comparaison essais/calculs a été réalisée dans le cas des essais en mode II. La
comparaison s’est limitée au cas des éprouvettes polymères. Une bonne corrélation
essais/calculs a été obtenue dans le cas de l’éprouvette homogène. Une comparai-
son n’était pas cependant envisageable dans le cas des éprouvettes composites car
les propriétés matériaux ainsi que la mésostructure ne sont pas connues. Dans la
dernière partie un film spécifique a été utilisé afin d’évaluer la taille de la zone de
contact.
237
238
Sommaire
6.1 Méthodologie multi-échelle de calculs sur pièces . . . . . . . 240
6.1.1 Description du modèle éléments finis . . . . . . . . . . . . . . 240
6.1.2 Cycles de chargement . . . . . . . . . . . . . . . . . . . . . . 241
6.1.3 Résultats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 241
6.2 Etude de la nocivité des défauts dans les billes céramiques. 247
6.3 Applications sur matériaux composites idéalisés . . . . . . . 250
6.4 Synthèse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 254
239
x2
x3
x1
(a) (b)
Figure 6.1 – Modèles éléments finis des portées extrados de l’aube et du disque (a)
aube ; (b) disque.
240
2.5
1.5
0.5
−0.5
−1
−1.5
Force normale (normalisee par F0 )
−2 Force tangentielle X (normalisee F0 )
Force tangentielle Y (normalisee F0 )
−2.5
0 3 6 9 12 15
Increments Chargement
Dans la suite, les résultats ne seront présentés que pour l’incrément de charge-
ment N ◦ 6 (cf. Fig.6.2).
6.1.3 Résultats
Les surfaces de la portée seront reconstruites pour leur utilisation dans le code
SA. Elles seront définies par des plans et des rayons en sortie de portée. Le maillage
du code semi-analytique utilisera 256 points dans la direction x (direction radiale
du moteur) et 512 points dans la direction y (direction axiale du moteur). La taille
de maille sera de 31 microns suivant x et de 282 microns suivant y.
– Champs de pression
241
242
1.2 1.2
1 1
0.8 0.8
P /P0 P /P0
0.6 0.6
0.4 0.4
0.2 0.2
0 0
−2 −1 0 1 2 −40 −30 −20 −10 0 10 20 30 40
x/a∗ y/a∗
(a) (b)
relativement raisonnables. Les mêmes calculs d’usure sur des milliers de cycles avec
la méthode éléments finis s’avèrent extrêmement coûteux en temps de calcul. Les
lois d’usure utilisées sont obtenues sur des essais bien instrumentés et peuvent être
facilement intégrées au code de calcul semi-analytique. La méthode semi-analytique
peut être vue comme une méthode de zoom permettant d’accélérer les calculs de
contact. Ainsi il est possible de simuler 60000 cycles d’usure avec bouclage (bouclage
structure (EF)/contact (SA)) tous les 100 cycles. Le bouclage consiste donc à in-
troduire les géométries usées, obtenues avec le code semi-analytique, dans le modèle
structurel éléments finis afin d’avoir les nouveaux torseurs d’efforts.
243
244
0.7
0.6 0.6
0.5
0.4 0.4
qx /µP0 qx /µP0
0.3
0.2 0.2
0.1
0 0
−2 −1 0 1 2 −40 −30 −20 −10 0 10 20 30 40
x/a∗ y/a∗
(a) (b)
0.15 0.15
0.1 0.1
qy /µP0 qy /µP0
0.05 0.05
0 0
−2 −1 0 1 2 −40 −30 −20 −10 0 10 20 30 40
x/a∗ y/a∗
(c) (d)
(a) (b)
245
(a) (b)
246
247
Piste Piste
∗ ∗
2 x1 2 x1
− −
− −
défaut
défaut
x3 x3
Bille Bille
(a) (b)
Figure 6.11 – Étude de défauts dans les billes céramiques (a) Cas d’une
hétérogénéité sphérique (b) Cas d’une hétérogénéité ellipsoı̈dale.
248
1 1
F0 F0
0.8 2.01 ∗ F0 0.8 2.01 ∗ F0
5.6 ∗ F0 5.6 ∗ F0
0.6 77.5 ∗ F0 0.6 77.5 ∗ F0
0.4 0.4
0.2 0.2
σ11 /Re σ22 /Re
0 0
−0.2 −0.2
−0.4 −0.4
−0.6 −0.6
−0.8 −0.8
−1 −1
0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4
x3 /a∗ x3 /a∗
(a) (b)
1 0
F0
0.8 2.01 ∗ F0
0.6 5.6 ∗ F0
77.5 ∗ F0 −0.2
0.4
0.2
σ33 /Re −0.4
0
−0.2 σV Ms /Re
−0.4 −0.6
−0.6
−0.8
−0.8 F0
−1 2.01 ∗ F0
−1.2 5.6 ∗ F0
77.5 ∗ F0
−1.4 −1
0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4
x3 /a∗ x3 /a∗
(c) (d)
Figure 6.12 – Valeur adimensionnée des contraintes suivant l’axe Z dans le cas
d’un massif homogène (sans défauts).
249
1 1
0.5 0.5
β11 0
β22 0
−0.5 −0.5
−1 −1
0 0.5 1 1.5 2 0 0.5 1 1.5 2
(dx3 − r)/a∗ (dx3 − r)/a∗
(a) (b)
1 1
0.5 0.5
β33 0
βV M 0
−0.5 −0.5
−1 −1
0 0.5 1 1.5 2 0 0.5 1 1.5 2
(dx3 − r)/a∗ (dx3 − r)/a∗
(c) (d)
250
(a) (b)
251
1.2
0.8
P/P0
0.6
0.4
0.2
t=0
t = τ /5
0
−1.5 −1 −0.5 0 0.5 1 1.5
x1 /a∗
(a) (b)
(a) (b)
252
(a) (b)
(a) (b)
253
6.4 Synthèse
Le modèle semi-analytique a été appliqué sur 3 cas industriels différents. Un
couplage entre un modèle éléments-finis et le code semi-analytique a tout d’abord été
réalisé pour résoudre le contact aube/disque. Le modèle éléments finis est utilisé pour
simuler une phase de décollage et atterrissage moteur. Les torseurs des efforts sont
récupérés du modèle éléments finis et introduits en entrée du code semi-analytique
avec un maillage plus fin, ce qui permet de mieux décrire les champs de pression, de
cisaillements et les contraintes en sous-couche. Le code semi-analytique a été ensuite
utilisé pour réaliser une étude qualitative en ce qui concerne la nocivité des défauts
dans les billes céramiques de roulements hybrides. Un premier critère de contrôle des
billes a été proposé. Cette étude a permis de réduire considérablement le volume de
la zone à contrôler. Le modèle a enfin été appliqué sur un matériau composite idéalisé
constitué de renforts sphériques noyés dans une matrice viscoélastique semi-infini.
254
Synthèse
La problématique du contact entre matériaux hétérogènes est de plus en
plus présente dans de nombreuses applications industrielles : roulements hy-
brides pour l’aéronautique et l’aérospatiale, contact entre matériaux composites
. . . L’introduction des aubes de soufflante en matériaux composites dans les nou-
veaux moteurs LEAP, développés par la société SNECMA du groupe SAFRAN et
GE Aviation, a permis par exemple de réduire de près de 500 Kg le poids de l’avion. Il
s’agit donc d’une véritable rupture technologique qui doit être suivie d’une évolution
des modèles de résolution des problèmes de contact, plus précisément de contact
entre matériaux hétérogènes. Dans ce manuscrit, une méthodologie de résolution de
contact entre matériaux hétérogènes est présentée. Les méthodes semi-analytiques
sont utilisées pour la résolution du problème contact. Les modèles proposés sont
systématiquement validés par comparaison avec la méthode des éléments finis. Le
temps de calcul ainsi que l’espace mémoire nécessaire sont considérablement réduits
par rapport aux éléments finis.
Dans le premier chapitre, le contexte industriel de la thèse a été introduit ainsi
que la problématique du contact entre matériaux hétérogènes dans le cas du contact
aube/disque, roulements hybrides, . . . . Les différents modes d’endommagements du
contact aube/disque ont été présentés. Une classification des matériaux composites
suivant plusieurs critères a été réalisée afin de mieux cerner la spécificité des com-
posites CMO et CMC utilisés par SNECMA. Les différentes méthodologies multi-
échelles utilisées dans la modélisation des matériaux hétérogènes ont été présentées.
Les méthodes de résolution des problèmes de contact ont été ensuite détaillées. Un
accent particulier a été mis sur la difficulté inhérente au contact entre matériaux
hétérogènes.
Dans le second chapitre, les bases de la méthode semi-analytique ont été passées
en revue afin de faciliter la compréhension du reste du manuscrit. Les algorithmes de
gradient conjugué CGM et les Transformées de Fourrier Rapides furent présentés.
L’utilisation de ces transformées de Fourier permet de réduire considérablement les
temps de calculs.
Le troisième chapitre s’intéresse à la prise en compte d’hétérogénéités iso-
tropes/anisotropes inclinées dans la résolution des problèmes de contact. Le mas-
255
sif est supposé élastique. Une étude fine, de l’effet du module de compressibi-
lité/cisaillement de l’hétérogénéité sur le problème de contact a été réalisée. Une
méthode numérique permettant de prendre en compte l’interaction entre plusieurs
hétérogénéités a été ensuite présentée. Un accent particulier a été mis sur la réduction
du temps de calcul et la robustesse du code. La mise en place d’une nouvelle méthode
de décomposition permet de diviser le temps de calcul par 2 ou 4 selon les configura-
tions. Une parallélisation OpenMP a été également introduite permettant d’avoir un
gain de 1.87 sur 2 coeurs et de 3.55 sur 4 coeurs. L’algorithme d’une parallélisation
MPI a été également présenté.
Le chapitre 4 s’intéresse à la théorie du contact viscoélastique. Une approche
originale basée sur les méthodes semi-analytiques a été mise en place pour résoudre
le problème de contact entre matériaux viscoélastiques, homogènes ou hétérogènes.
Le modèle peut traiter à la fois des cas d’indentation, de roulement/glissement
et de fretting. Le modèle fournit des solutions en termes de champs de pression,
de cisaillements, de coefficients de frottement apparent mais aussi et surtout les
champs de contraintes en sous-couche. Les résultats obtenus sont valables aussi bien
en régime permanent que transitoire.
Le chapitre 5 présente essentiellement le dispositif expérimental développé en
laboratoire pour mesurer les champs de déplacements à l’interface des corps en
contact. La technique de corrélations d’images a été utilisée. Le protocole d’essais
ainsi que les éprouvettes utilisées ont été présentés. Les résultats obtenus seront
comparés à ceux obtenus au travers de la simulation numérique. Le dispositif a fait
l’objet d’un dépôt de brevet SNECMA/LAMCOS.
Dans le dernier chapitre, le modèle sera appliqué sur trois cas industriels
différents. Un prototype de couplage entre un modèle éléments finis et le code semi-
analytique sera réalisé pour résoudre le contact aube/disque. Il s’agit d’une stratégie
multi-échelle de calculs sur pièce. Le modèle a été ensuite utilisé pour étudier la
nocivité des défauts dans les billes céramiques. Une application dans le cas d’un
matériau composite idéalisé sera présentée. Les résultats dans le cas de tissages 3D
SNECMA/Safran Composites ne pouvant être présentés pour des raisons de confi-
dentialités.
Perspectives
Les perspectives pouvant être apportées à ce travail de thèse sont nombreuses.
Un premier chantier concerne la mise en place d’une chaı̂ne automatisée, partant de
la description de la microstructure des composites tissés 3D via le logiciel WiseTex
jusqu’à la résolution du problème de contact avec la méthode semi-analytique. Les
fibres anisotropes du composite tissé peuvent être décrites par une série de cubes
ou d’ellipsoı̈des. Cette description se limitera à l’échelle mésoscopique. Le défi scien-
tifique n’est pas énorme, cependant ce travail est indispensable d’un point de vue
industriel.
256
257
258
Coefficients d’influences
x23 + x22 − x2 ρ 1 − 2ν ρ − x2 + x3
ν x3 x1 x2
Sxp1 x1 (x1 , x2 , x3 , E, ν) = arctan + arctan + ,
π x3 x1 π x1 2π (x21 + x23 )ρ
(A.2)
x23 + x22 − x2 ρ 1 − 2ν ρ − x1 + x3
ν x3 x1 x2
Sxp2 x2 (x1 , x2 , x3 , E, ν) = arctan + arctan + ,
π x3 x1 π x2 2π (x22 + x23 )ρ
(A.3)
259
x23 + x22 − x2 ρ
1 x3 x1 x2 1 1
Sxp3 x3 (x1 , x2 , x3 , E, ν) = arctan − + ,
2π x3 x1 2π ρ x21 + x23 x22 + x23
(A.4)
x3 1 1 − 2ν
Sxp1 x2 (x1 , x2 , x3 , E, ν) = − − ln (ρ + x3 ) , (A.5)
2π ρ 2π
x23 x1
Sxp2 x3 (x1 , x2 , x3 , E, ν) = , (A.6)
2π (x22 + x23 )ρ
x23 x2
Sxp1 x3 (x1 , x2 , x3 , E, ν) = . (A.7)
2π (x1 + x23 )ρ
2
La contribution des cisaillements qx1 sur les contraintes du massif semi-infini est
donnée par
−x21 + x3 x2
x3 1 ν x2 1
Sxqx1
1 x1
(x1 , x2 , x3 , E, ν) = − 1+ + − ln (ρ − x2 ) ,
2π ρ (ρ + x3 )(ρ − x2 ) π ρ + x3 π
(A.9)
x3 x2 ν x2
Sxqx1
2 x2
(x1 , x2 , x3 , E, ν) =− − + ln (ρ − x2 ) , (A.10)
2π ρ(ρ + x3 ) π ρ + x3
x23 x2
Sxqx1
3 x3
(x1 , x2 , x3 , E, ν) = , (A.11)
2π ρ(x1 + x23 )
2
x3 x1 ν x1 1
Sxqx1
1 x2
(x1 , x2 , x3 , E, ν) = − − − ln (ρ − x1 ) , (A.12)
2π ρ(ρ + x3 ) π ρ + x3 2π
x3 1
Sxqx1
2 x3
(x1 , x2 , x3 , E, ν) = − , (A.13)
2π ρ
260
x23 + x22 − x2 ρ
x3 x1 x2 1
Sxqx1
1 x3
(x1 , x2 , x3 , E, ν) = 2 2 + arctan . (A.14)
2π ρ(x1 + x3 ) 2π x3 x1
La contribution des cisaillements qx2 sur les contraintes du massif semi-infini est
donnée par
x3 x1 ν x1
Sxqx2
1 x1
(x1 , x2 , x3 , E, ν) =− − + ln (ρ − x1 ) , (A.16)
2π ρ(ρ + x3 ) π ρ + x3
−x22 + x3 x1
x3 1 ν x1 1
Sxqx2
2 x2
(x1 , x2 , x3 , E, ν) = 1+ + − ln (ρ − x1 ) ,
2π ρ (ρ + x3 )(ρ − x1 ) π ρ + x3 π
(A.17)
x23 x1
Sxqx2
3 x3
(x1 , x2 , x3 , E, ν) = (A.18)
2π ρ(x2 + x23 ),
2
x3 x2 ν x2 1
Sxqx2
1 x2
(x1 , x2 , x3 , E, ν) = − − − ln (ρ − x2 ) , (A.19)
2π ρ(ρ + x3 ) π ρ + x3 2π
x23 + x21 − x1 ρ
x3 x2 x1 1
Sxqx2
2 x3
(x1 , x2 , x3 , E, ν) = 2 2 + arctan , (A.20)
2π ρ(x2 + x3 ) 2π x3 x2
x3 1
Sxqx2
1 x3
(x1 , x2 , x3 , E, ν) = − . (A.21)
2π ρ
Identiquement aux p contraintes, on peut définir les déplacements élastiques en
surface. On note ρ̄ = x21 + x22 . La contribution des pressions p sur les déplacements
élastiques en surface du massif semi-infini est donnée par
261
avec
(1 + ν)(1 − 2ν) ρ̄ − x2
Uxp1 (x1 , x2 , E, ν) = − 2x1 arctan − x2 ln ρ̄ , (A.23)
2πE x1
(1 + ν)(1 − 2ν) ρ̄ − x1
Uxp2 (x1 , x2 , E, ν) =− 2x2 arctan − x1 ln ρ̄ , (A.24)
2πE x2
(1 − ν 2 )
Uxp3 (x1 , x2 , E, ν) = − (x2 ln (ρ̄ − x1 ) + x1 ln (ρ̄ − x2 )) . (A.25)
πE
La contribution des pressions qx1 sur les déplacements élastiques en surface du
massif semi-infini est donnée par
1 − ν2 1+ν
Uxqx1 (x1 , x2 , E, ν) = − x1 ln (ρ̄ − x2 ) − ln (ρ̄ − x1 ) , (A.27)
πE πE
ν(1 + ν)
Uxqx1 2 (x1 , x2 , E, ν) = − ρ̄, (A.28)
πE
(1 + ν)(1 − 2ν) ρ̄ − x2
Uxq3x1 (x1 , x2 , E, ν) = −2x1 arctan + x2 ln ρ̄ . (A.29)
2πE x1
La contribution des pressions qx2 sur les déplacements élastiques en surface du
massif semi-infini est donnée par
ν(1 + ν)
Uxq1x2 (x1 , x2 , E, ν) = − ρ̄, (A.31)
πE
262
1 − ν2 1+ν
Uxq2x2 (x1 , x2 , E, ν) = − x2 ln (ρ̄ − x1 ) − ln (ρ̄ − x2 ) , (A.32)
πE πE
(1 + ν)(1 − 2ν) ρ̄ − x1
Uxq3x2 (x1 , x2 , E, ν) = −2x2 arctan + x1 ln ρ̄ . (A.33)
2πE x2
1
Dijkl = {ψ,klij − 2νδkl φ,ij
8π (1 − ν)
− (1 − ν) [φ,kj δil + φ,ki δjl + φ,lj δik + φ,li δjk ]} (A.34)
Inclusion ellipsoı̈dale
Eshelby fut un des premiers à montrer que les coefficients d’influence Dijkl à
l’intérieur de l’inclusion sont constants. Considérons une ellipsoı̈de de demi-axe a1 ,
a2 , a3 suivant x1 ,x2 , et x3 respectivement. Les fonctions potentielles peuvent alors
s’exprimer comme,
Z ∞
U ds
φ = πa1 a2 a3 (A.35)
λ ∆
1 ∞ ∂ s2 U 2 1 U 2s
Z
ψ = πa1 a2 a3 − ds (A.36)
2 λ ∂s ∆ 2 ∆
Où,
xi xi
U (s) = 1 − 2 (A.37)
a +s
q i
∆ (s) = (a21 + s) + (a22 + s) + (a23 + s) (A.38)
263
φ=V (A.42)
a2i Vi
ψ,i = xi V − (A.43)
Avec,
Z ∞
U (s) 1
V (x) = πa1 a2 a3 ds = [I (λ) − xr xr Ir (λ)] (A.44)
λ ∆ (s) 2
Z ∞
U (s) 1
Vi (x) = πa1 a2 a3 2
ds = [Ii (λ) − xr xr Iri (λ)] (A.45)
λ (ai + s) ∆ (s) 2
Où,
s
a21 − a23 π
0 < Θ = arcsin 2
< (A.48)
a1 + λ 2
s
a21 − a22
K= <1 (A.49)
a21 − a23
s
a21 − a23
w = arcsin (A.50)
a21 + s
264
265
Où,
s
a21 + λ
bp = (A.63)
a23 + λ
s
a21 − a23
dp = (A.64)
a23 + λ
– Sphère (a1 = a2 = a3 = a)
4πa3
I (λ) = √ , (A.65)
a2 + λ
4πa3
I1 (λ) = I2 (λ) = I3 (λ) = 3 (A.66)
3 (a2 + λ) 2
266
Dispositif expérimental
B.1 Objectifs.
Cet annexe est fortement complémentaire du chapitre 5. Le chapitre 5 s’est li-
mité, pour des raisons de concisions, à une description très succincte du dispositif
expérimental.
267
Précises et compactes, ces capteurs remplacent les montages complexes utilisant plusieurs
capteurs mono axe. Ils s'installent facilement sur les machines-outils et les bancs de test.
Les gammes standard des capteurs 6 axes sont de 440 N à 4400 N.
Les capteurs d'effort multiaxe AMTI sont de très grandes performances et ont une conception
qui permet d’avoir un couplage entre voies très faible.
Cette génération de capteurs de force multiaxe présente de très bonnes sensibilité et précision
de mesure dès le premier Newton quelques que soit la gamme et le modèle choisi.
Nous proposons des capteurs multiaxes IP 67 ou IP68 pour une immersion sous 10 m d’eau.
268
Caractéristiques et dimensions
Caractéristiques / Model 100 250 500 1000
Etendue de mesure Fz (N) 440 1100 2200 4400
Etendue de mesure Fx,Fy (N) 220 560 1100 2200
Etendue de mesure Mz (Nm) 5.6 14 28 56
Etendue de mesure Mx,My (Nm) 11 28 56 110
Sensibilité Fz (µV/V/N) 1.35 0.67 0.34 0.17
Sensibilité Fx,Fy (µV/V/N) 5.4 2.7 1.35 0.67
Sensibilité Mz (µV/V/Nm) 121 106 53 26
Sensibilité Mx,My (µV/V/Nm) 266 138 71 35
Rigidité Fz (x10^7 N/m) 2.8 7.5 15.0 30.0
Rigidité Fxy (x10^7 N/m) 0.2 0.5 1.0 2.0
Rigidité Mz (x10^4 Nm/rad) 0.2 0.5 1.1 2.2
Fréquence de résonnance Mxy (Hz) 300 500 700 1000
269
270
271
272
VN-8MC / VN-11MC
VN-16MC / VN-29MC
NANO STAGE PIXEL SHIFT CAMERA FOR EXTENDED RESOLUTIONS
VN Series pixel shift cameras are designed for applications where the object is stationary and extremely
high resolution is required. Vieworks advanced pixel shift technology based on a precise piezoelectric stage
allows image captures as high as 99 million pixels using the VN-11MC cameras. Even higher resolutions, up
to 260 million pixels can be obtained using the VN-29MC camera. These cameras are ideal for applications
such as FPD inspection, document/film scanning, research and scientific imaging.
W W W. S T E M M E R- I M A G I N G . C O M
[Link]
273
VN-8MC / VN-11MC
VN-16MC / VN-29MC
NANO STAGE PIXEL SHIFT CAMERA FOR EXTENDED RESOLUTIONS
Specifications
Model VN-8M VN-11M VN-16M VN-29M
×1 Mode 3296 × 2472, 8.1M 4008 × 2672, 10.7M 4872 × 3248, 15.8M 6576 × 4384, 28.8M
Resolution ×4 Mode 6592 × 4944, 32.6M 8016 × 5344, 42.8M 9744 × 6496, 63.3M 13152 × 8768, 115.3M
(H × V)
×9 Mode 9888 × 7416, 73.3M 12024 × 8016, 96.4M 14616 × 9744, 142.4M 19728 × 13152, 259.5M
Sensor
KAI-08050 KAI-11002 KAI-16000 KAI-29050
(Truesense Imaging)
Sensor Size
4/3㶆 㪊㪌㩷䋘 㪊㪌㩷䋘 㪊㪌㩷䋘
(Optical Format)
Sensor Type Progressive Scan Interline Transfer CCD
Pixel Size 5.5 ঙ × 5.5 ঙ 9.0 ঙ × 9.0 ঙ 7.4 ঙ × 7.4 ঙ 5.5 ঙ × 5.5 ঙ
8.1M Mode : 16.3 fps 10.7M mode : 6.4 fps 15.8M Mode : 4.2 fps 28.8M Mode: 5 fps
Max. Frame Rate
32.6M Mode : 4.1 fps 42.8M mode : 1.6 fps 63.3M Mode : 1.1 fps 115.3M Mode: 1.3 fps
(40 ঐ)
73.3M Mode : 1.8 fps 96.4M mode : 0.7 fps 142.4M Mode : 0.5 fps 259.5M Mode: 0.6 fps
Electronic Shutter Global Shutter
Pixel Data Format 8/10/12 bit
Interface Base Camera Link
Data Output
40/80 ঐ 30/40 ঐ 30/40 ঐ 40/80 ঐ
Pixel Clock Speed
Free-Run , Overlap, Fast, Double
Trigger Mode
Programmable Exposure Time and Trigger Polarity
Dynamic Range 62
Shift Range 0 ~ 15 ঙ, 1 ঘ step
Shift Resolution 0.001 ঙ
Shift Control Manual Mode or Sequence Mode (4/9 Shot Mono, 4/16/36 Shot Color)
Shift Latency <8
Dimension / Weight 90 䋘 × 90 䋘 × 123.5 䋘㩷/ 1.2 䋗
Temperature Operating : 10㷄 ~ 40㷄, Storage : -30㷄 ~ 65㷄
Lens Mount F-mount, Custom mount available upon request
10 ~ 14 V DC 10 ~ 14 V DC 10 ~ 14 V DC 10 ~ 14 V DC
Power
MAX. 6 W MAX. 10 W MAX. 8 W MAX. 12 W
Configuration
Configurator
Software
W W W. S T E M M E R- I M A G I N G . C O M
274
W W W. S T E M M E R- I M A G I N G . C O M
275
Mechanical Dimensions
Unit : mm
forerrors
acceptedfor
liabilityisisaccepted whichmay
errorswhich containedininthis
becontained
maybe document.
thisdocument.
For more information please contact local distributor or
visit our website at [Link]
Noliability
Reproduction in whole or in part without written permission is prohibited. Vieworks Co., Ltd. is not responsible for any technical or
typographical errors and reserves the right to make changes to products, specifications and documentation without prior notice.
[Link]
withoutnotice.
RA14-11A-016
Subjecttototechnical
K-SONO28-11/2012∙ ∙Subject
K-VIEW5-09/2013 changewithout
technicalchange
W W W. S T E M M E R- I M A G I N G . C O M I M A G I N G I S O U R PA S S I O NISO-9001, ISO-13485
GERMANY UNITED KINGDOM FRANCE SW I T Z E R L A N DCorporate Headquarters
THE NETHERLANDS
AU ST R I A IRELAND L I E C H T E N ST E I N
Vieworks Co., Ltd. # 604 Suntechcity 㪠㪠, 307-2 Sangdaewon-dong, Jungwon-gu,B ESeongnam-si,
LG I U M L U X E M B O U R G
Gyeonggi-do, 462-806 South Korea
Phone: +49 89 80902-0 Phone: +44 1252 780000 Phone: +33 1 45069560 Phone: +41 55 415 90 90 Phone: +31 575 495159
info@[Link] tel +82-70-7011-6161 info@[Link]
info@[Link] fax +82-31-737-4936 info@[Link]
e-mail sales@[Link] website [Link]
info@[Link]
276
Optical Components
Ask About Our Build-to-Print and Custom Capabilities
Windows and
Optical Flats
OEM
Retroreflectors
Prisms and
Spherical
Lenses
Clear Aperture
$ All CVI Melles Griot low loss, high energy AR coatings available
$ CVI Melles Griot beamsplitter coatings available
$ CVI Melles Griot partially reflective coatings available
Laser Grade Square Windows
Mirrors
N-BK7
12.7 6.4 l/10 SQW-0525-C
Window dimensions layout 25.4 2.0 l/10 SQW-1008-C
25.4 6.4 l/4 SQW-1025-C
25.4 9.5 l/10 SQW-1037-C
Components
Polarization
277
278
279
B.3 Eprouvette
B.3.1 Billes céramiques.
280
B.3.4 Dimensions.
281
Gamme d'usinage
51,150
2
50
0.1
50
282
B.4 Tomographie
B.4.1 Phoenix v—tome—x s, système informatisé utilisé
pour la tomographie.
283
284
285
286
287
288
289
B.5 Essais
B.5.1 Réglages acquisitions images
Le réglage de la ”mise au point” se fait facilement avec une molette de réglage
assez précise. Le réglage de la luminosité peut se faire de trois manières différentes :
– Ouverture du diaphragme (la netteté étant meilleure quand l’ouverture est
maximum) ;
– Réglage de la luminosité de la lampe (autour de 80 − 90%) ;
– Réglage du temps d’exposition (de l’ordre de 0.1 ms suivant les paramètres
précédents) ;
La figure B.15 ci-dessous montre une image avec un bon réglage et son histogramme
de répartition des pixels par niveau de gris. La fonction saturation permet de montrer
les zones saturée des deux extrémités (endroit o le niveau de gris est maximum ou
minimum). Ces zones sont à minimiser pour permettre une bonne corrélation des
images.
290
Il est important de :
291
– choisir une zone d’étude pas trop large avec la fonction polygone, le minimum
nécessaire à l’exploitation des résultats (la surface étant bombée, les bords sont
plus flous)
– Déplacer le point de référence hors de la zone de contact
– Choisir un Subset d’environ 90 +/- 10
– Choisir un Step ≤ 25
292
Pressure -micro ®
®
Common Applications
GASKETED INTERFACE
Packaging & Converting:
QLSUROOHULPSUHVVLRQVKHDWVHDOLQJODPLQDWLRQ
Automotive:
EUDNHSDGFODPSLQJFOXWFKEDWWHU\ IXHOFHOOLPSDFW
JDVNHWEROWHGMRLQWODPLQDWLRQZHOGLQJZLSHU
Electronics:
KHDWVLQN%*$FRQQHFWRUODPLQDWLRQ/&'
ERQGLQJZDIHUERQGLQJSROLVKLQJVRODUFHOOV
TIRE TREAD
Aerospace:
FRPSRVLWHOD\XSIXHOFHOOODPLQDWLRQLPSDFW
EROWHGLQWHUIDFH
Ergonomics:
ELRPHFKDQLFVERG\PDSSLQJ
293
STEP 1 STEP 2
Place Pressurex-micro®*UHHQ¿OPEHWZHHQ 7KHVHQVRU¿OPFDSWXUHVDSHUPDQHQWLPDJH
any two contacting surfaces of relative pressure distribution
Zo
om High
pressure
ed
vie
w
Heat sealer
Low
Pressurex-micro® Green pressure
placed here
How It Works
ROLL DIMENSION PRESSURE RANGE
6LPSO\ SODFH 3UHVVXUH[PLFUR® *UHHQ SUHVVXUH LQGLFDWLQJ FILM TYPE
UHYHDOHG3UHVVXUH[PLFUR®*UHHQLVH[WUHPHO\WKLQDQGÀH[LEOH PMG-Custom* ~ ~ ~
PP ZKLFKHQDEOHVLWWRFRQIRUPWRFXUYHGVXUIDFHV,WLV *&XVWRP¿OPW\SHVDYDLODEOHXSRQUHTXHVW
LGHDOIRULQYDVLYHLQWROHUDQWHQYLURQPHQWVDQGWLJKWVSDFHVQRW
DFFHVVLEOHE\FRQYHQWLRQDOHOHFWURQLFWUDQVGXFHUV
The Science
0LFURVFRSLF SLJPHQWHG SDUWLFOHV DGKHUHG WR WKH GRQRU VXE SPECIFICATIONS
VWUDWHDUHDWWUDFWHGWRWKHFKHPLFDOO\VXUIDFHWUHDWHGUHFHLYHU OPERATING TEMPERATURE &WR&
VKHHW 7KH UHFHLYHU VKHHW FRDWLQJ KDV EHHQ IRUPXODWHG
WR FUHDWH D QRQFRYDOHQW ERQG ZLWK WKH GRQRU FRDWLQJ DV HUMIDITY RANGE 10% to 90% RH
294
Figure B.19 – Eprouvette hétérogéne contenant une seule bille céramique pour des
essais en mode II, efforts normaux de 500N , 1000N , 1500N et 2000N : εxx , εyy et
εxy .
295
296
297
298
299
300
[CHA 12b] Chaise T., Li J., Nelias D., Kubler R., Taheri S., Douchet G.,
Robin V., Gilles P.
Modelling of multiple impacts for the prediction of distortions and residual stresses
induced by Ultrasonic Shot Peening (USP). Journal of Materials Processing Tech-
nology, vol. 212, 2012, p. 2080-2090, Elsevier.
[CHE 08a] Chen P., Xiong J., Shen Z.
Thickness effect on the contact behavior of a composite laminate indented by a
rigid sphere. Mechanics of Materials, vol. 40, no 4, 2008, p. 183–194, Elsevier.
[CHE 08b] Chen W., Liu S., Wang Q.
FFT-based numerical methods for elasto-plastic contacts of nominally flat sur-
faces. Journal of Applied Mechanics, vol. 75, 2008, p. 011022-1-11.
[CHE 11] Chen W., , Wang Q., Huan Z., Luo X.
Semi-analytical viscoelastic contact modeling of polymer-based materials. Journal
of Tribology, vol. 133, 2011, page 041404.
[CHI 77a] Chiu Y.
On the stress field due to initial strains in a cuboid surrounded by an infinite
elastic space. Journal of Applied Mechanics, vol. 44, 1977, page 587.
[CHI 77b] Chiu Y.
On the stress field due to initial strains in a cuboid surrounded by an infinite
elastic space. Journal of Applied Mechanics, vol. 44, 1977, p. 587-590.
[CHI 78a] Chiu Y.
On the stress field and surface deformation in a half space with a cuboidal zone
in which initial strains are uniform. Journal of Applied Mechanics, vol. 45, 1978,
p. 302-306.
[CHI 78b] Chiu Y.
On the stress field and surface deformation in a half space with a cuboidal zone
in which initial strains are uniform. Journal of Applied Mechanics, vol. 45, 1978,
page 302.
[COR 99] Cormier N. G., Smallwood B. S., Sinclair G. B., Meda G.
Aggressive submodelling of stress concentrations. International Journal for Nu-
merical Methods in Engineering, vol. 46, no 6, 1999, p. 889–909.
[COU 05] Courbon J., Lormand G., Dudragne G., Daguier P., Vincent
A.
Influence of inclusion pairs, clusters and stringers on the lower bound of the endu-
rance limit of bearing steels. Tribology International, vol. 36, 2005, p. 921-928.
[COU 08] Couégnat G.
Approche multiéchelle du comportement mécanique de matériaux composites à
renfort tissé. Thèse de doctorat, Université Sciences et Technologies-Bordeaux I,
2008.
[CRO 56] Crossland B.
Effect of large hydrostatic pressures on the torsional fatigue strength of an alloy
301
steel. Proc. Int. Conf. on Fatigue of Metals, vol. 138 Institution of Mechanical
Engineers London, 1956.
[DAN 93] Dang-Van K.
Macro-micro approach in high-cycle multiaxial fatigue. ASTM Special Technical
Publication, vol. 1191, 1993, p. 120–120, AMERICAN TECHNICAL PUBLI-
SHERS LTD.
[DIC 06a] Dick T., Cailletaud G.
Fretting modelling with a crystal plasticity model of Ti6Al4V. Computational
Materials Science, vol. 38, no 1, 2006, p. 113–125.
[DIC 06b] Dick T., Paulin C., Cailletaud G., Fouvry S.
Experimental and numerical analysis of local and global plastic behaviour in fret-
ting wear. Tribology International, vol. 39, no 10, 2006, p. 1036–1044.
[DIC 06c] Dick T.
Modélisation multiéchelle du phénomene de fretting dans le contact aube-disque.
Thèse de doctorat, École Nationale Supérieure des Mines de Paris, 2006.
[DIN 04] Dini D., Nowell D.
Flat and rounded fretting contact problems incorporating elastic layers. Interna-
tional Journal Of Mechanical Sciences, vol. 46, 2004, page 16351657.
[DRA 01] Dragon-Louiset M.
On a predictive macroscopic contact-sliding wear model based on micromechanical
considerations. International journal of solids and structures, vol. 38, no 9, 2001,
p. 1625–1639, Elsevier.
[DUB 03] Dubourg M., Chateauminois A., Villechaise B.
In situ analysis and modeling of crack initiation and propagation within model
fretting contacts using polymer materials. Tribology International, vol. 36, no 2,
2003, p. 109–119, Elsevier.
[DVO 92a] Dvorak G. J.
Transformation field analysis of inelastic composite materials. Proceedings of the
Royal Society of London A : Mathematical, Physical and Engineering Sciences,
vol. 437 The Royal Society, 1992, p. 311–327.
[DVO 92b] DVORAK G. J., BENVENISTE Y.
On transformation strains and uniform fields in multiphase elastic media. Pro-
ceedings of the Royal Society of London A, vol. 437, 1992, page 291310.
[DVO 94] Dvorak G .-J., Bahei-el din Y., Wafa A.
The modelling of inelastic composite materials with the transformation field ana-
lysis. Model. Simulat. Mater. Sci. Eng, vol. 2, 1994, page 571586.
[ELH 06] El Hage C.
Modélisation du comportement élastique endommageable de matériaux compo-
sites à renfort tridimensionnel. Thèse de doctorat, 2006.
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
RÉSUMÉ :
Cette thèse s’intéresse à la problématique du contact entre matériaux hétérogènes. L’industrie (automobile,
aéronautique, spatiale, . . . ) s’intéresse de plus en plus à ces types de matériaux. Il s’agira par exemple des alliages
métalliques, des matériaux poreux, matériaux composites (composites tissés, interlocks 3D, interlocks 2D), des billes
céramiques contenant des impuretés (porosités/précipités),. . . Dans ce manuscrit, un modèle de contact basé sur les
méthodes semi-analytiques a été développé. Un algorithme de gradient conjugué est utilisé afin de résoudre rapidement
le problème de contact. Le modèle permet de prendre en compte la présence d’une ou de plusieurs hétérogénéités
isotropes/anisotropes dans le problème de contact. Une approche inspirée de la méthode de l’inclusion équivalente
proposée par Eshelby est utilisée dans le solveur de contact pour prendre en compte l’effet de ces hétérogénéités. Les
méthodes de transformées de Fourier rapides (FFT) permettent d’accélérer les calculs. Une méthode numérique a
été mise en œuvre afin de prendre en compte l’interaction entre plusieurs hétérogénéités. Le massif peut e élastique
ou viscoélastique. L’approche développée dans la thèse peut résoudre à la fois les problèmes d’indentation, de
roulement/glissement ou de fretting en présence de matériaux élastiques hétérogènes, viscoélastiques homogènes ou
hétérogènes. Les solutions sont données en termes de champs de pressions, de cisaillements et de contraintes. Dans
le cas des matériaux viscoélastiques le code de calcul est capable de fournir le coefficient de frottement apparent
ainsi que toutes les variables de contact aussi bien en régime permanent que transitoire. Le modèle a été validé
par comparaison avec la méthode des éléments finis classiques en utilisant le logiciel commercial Abaqus v6.11.
Le temps de calcul ainsi que l’espace mémoire nécessaire sont considérablement réduits par rapport à la méthode
éléments finis. La parallélisation a été introduite dans le code de contact afin de réduire toujours plus le temps de
calcul. Il s’agit d’un code robuste, rapide et facilement utilisable en Bureau d’Etudes. Une approche expérimentale
originale a été mise en place afin de mesurer les champs de déplacements à l’interface des corps en contact. De bonnes
corrélations essais/calculs ont été obtenues. Enfin quelques applications industrielles ont été présentées. Un cou-
plage entre un code éléments finis structurel et le code semi-analytique de résolution de contact a été également réalisé.
MOTS-CLéS : contact, fretting, semi-analytique, matériaux hétérogènes, composites tissés, viscoélasticité, anisotropie,
approche multi-échelle, couplage de code, parallélisation