Friction prediction for rough surfaces in an

elastohydrodynamically lubricated contact

Yuanyuan Zhang

To cite this version:

Yuanyuan Zhang. Friction prediction for rough surfaces in an elastohydrodynamically lubricated
contact. Materials and structures in mechanics [[Link]-ph]. Université de Lyon, 2019. English.
�NNT : 2019LYSEI063�. �tel-02956182�

HAL Id: tel-02956182

[Link]
Submitted on 2 Oct 2020

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
 N°d’ordre NNT : 2019LYSEI063

THESE de DOCTORAT DE L’UNIVERSITE DE LYON

opérée au sein de
l’INSA de Lyon

Ecole Doctorale N° EDA162

(Mécanique, Energétique, Génie Civil, Acoustique)

Spécialité de doctorat : Génie Mécanique

Soutenue publiquement le 05 septembre 2019, par :

Yuanyuan ZHANG

Friction Prediction for Rough Surfaces

in an Elastohydrodynamically
Lubricated Contact

Devant le jury composé de :

CAYER-BARRIOZ Juliette Directrice de recherche CNRS ECL Présidente

EVANS Pwt Professeur Cardiff University Rapporteur
KŘUPKA Ivan Professeur Brno University of Technology Rapporteur
VENNER Cornelis. H Professeur University of Twente Examinateur
BIBOULET Nans Maître de Conférences INSA Lyon Examinateur
LUBRECHT Antonius. A Professeur INSA Lyon Directeur de thèse

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
Cette thèse est accessible à l'adresse : [Link]
© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Département FEDORA – INSA Lyon - Ecoles Doctorales – Quinquennal 2016-2020

SIGLE ECOLE DOCTORALE NOM ET COORDONNEES DU RESPONSABLE

CHIMIE CHIMIE DE LYON M. Stéphane DANIELE

Institut de recherches sur la catalyse et l’environnement de Lyon
[Link] IRCELYON-UMR 5256
Sec. : Renée EL MELHEM Équipe CDFA
Bât. Blaise PASCAL, 3e étage 2 Avenue Albert EINSTEIN
secretariat@[Link] 69 626 Villeurbanne CEDEX
INSA : R. GOURDON
directeur@[Link]

E.E.A. ÉLECTRONIQUE, M. Gérard SCORLETTI

ÉLECTROTECHNIQUE, École Centrale de Lyon
AUTOMATIQUE 36 Avenue Guy DE COLLONGUE
69 134 Écully
[Link] Tél : [Link].97 Fax [Link].17
Sec. : M.C. HAVGOUDOUKIAN [Link]@[Link]
[Link]@[Link]

E2M2 ÉVOLUTION, ÉCOSYSTÈME, M. Philippe NORMAND

MICROBIOLOGIE, MODÉLISATION UMR 5557 Lab. d’Ecologie Microbienne
Université Claude Bernard Lyon 1
[Link] Bâtiment Mendel
Sec. : Sylvie ROBERJOT 43, boulevard du 11 Novembre 1918
Bât. Atrium, UCB Lyon 1 69 622 Villeurbanne CEDEX
Tél : [Link].62 [Link]@[Link]
INSA : H. CHARLES
secretariat.e2m2@[Link]

EDISS INTERDISCIPLINAIRE Mme Emmanuelle CANET-SOULAS

SCIENCES-SANTÉ INSERM U1060, CarMeN lab, Univ. Lyon 1
Bâtiment IMBL
[Link] 11 Avenue Jean CAPELLE INSA de Lyon
Sec. : Sylvie ROBERJOT 69 621 Villeurbanne
Bât. Atrium, UCB Lyon 1 Tél : [Link].09 Fax : [Link].16
Tél : [Link].62 [Link]@[Link]
INSA : M. LAGARDE
[Link]@[Link]

INFOMATHS INFORMATIQUE ET M. Luca ZAMBONI

MATHÉMATIQUES Bât. Braconnier
43 Boulevard du 11 novembre 1918
[Link] 69 622 Villeurbanne CEDEX
Sec. : Renée EL MELHEM Tél : [Link].52
Bât. Blaise PASCAL, 3e étage zamboni@[Link]
Tél : [Link].46
infomaths@[Link]

MATÉRIAUX DE LYON M. Jean-Yves BUFFIÈRE

Matériaux INSA de Lyon
[Link] MATEIS - Bât. Saint-Exupéry
Sec. : Stéphanie CAUVIN 7 Avenue Jean CAPELLE
Tél : [Link].70 69 621 Villeurbanne CEDEX
Bât. Direction Tél : [Link].70 Fax : [Link].28
[Link]@[Link] [Link]@[Link]

MEGA MÉCANIQUE, ÉNERGÉTIQUE, M. Jocelyn BONJOUR

GÉNIE CIVIL, ACOUSTIQUE INSA de Lyon
Laboratoire CETHIL
[Link] Bâtiment Sadi-Carnot
Sec. : Stéphanie CAUVIN 9, rue de la Physique
Tél : [Link].70 69 621 Villeurbanne CEDEX
Bât. Direction [Link]@[Link]
mega@[Link]

ScSo ScSo* M. Christian MONTES

Université Lyon 2
[Link] 86 Rue Pasteur
Sec. : Véronique GUICHARD 69 365 Lyon CEDEX 07
INSA : J.Y. TOUSSAINT [Link]@[Link]
Tél : [Link].76
[Link]@[Link]

Cette thèse est accessible à l'adresse : [Link]

*ScSo : Histoire, Géographie, Aménagement, Urbanisme, Archéologie, Science politique, Sociologie, Anthropologie
© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
Cette thèse est accessible à l'adresse : [Link]
© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Abstract

The friction of interfacial surfaces greatly influences the performance of mechanical ele-
ments. Friction has been investigated experimentally in most studies. In this work, the friction
is predicted by means of numerical simulation under an elastohydrodynamic lubrication (EHL)
rough contact condition.
The classical Multigrid technique performs well in limiting computing time and memory
requirements. However, the coarse grid choice has an important influence on code robustness
and code efficiency to solve the rough problem. In the first part of this work, a coarse grid con-
struction method proposed by Alcouffe et al. is implemented in the current time-independent
EHL Multi-Grid code. Then this modified solver is extended to transient cases to solve the
rough contact problem.
The friction curve is usually depicted as a function of “Λ ratio”, the ratio of oil film thickness
to root-mean-square of the surface roughness. However this parameter is less suitable to plot
friction variations under high pressure conditions (piezoviscous elastic regime). In the second
part of this work, the friction coefficient is computed using the modified EHL code for many op-
erating conditions as well as surface waviness parameters. Simulation results show that there
is no single friction curve when the old parameter "Λ ratio" used. Based on the Amplitude
Reduction Theory, a new scaling parameter depends on operating condition and waviness pa-
rameters is found, which can give a unified friction curve for high pressure situation.
For more complex rough surfaces, a power spectral density (PSD) based method is pro-
posed to predict friction variations in the third part of this work. The artificial surface rough-
ness is employed to test the rapid prediction method firstly. Good agreement is found between
the full numerical simulation and this rapid prediction. Then the rapid prediction method is
applied to analyze the friction variation of measured surface roughness. A comparison is also
made between predictions and experiments.
Both the new scaling parameter and the friction increase predicted by the PSD method
show good engineering accuracy for practical use.

Keywords: Elastohydrodynamic lubrication, Numerical simulation, Piezoviscous elastic regime,

Amplitude Reduction Theory, Friction variation

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
Cette thèse est accessible à l'adresse : [Link]
© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Résumé

Le frottement à l’interface des surfaces influence les performances des éléments mécaniques.
Le frottement a été étudié expérimentalement dans la plupart des études. Dans ce travail,
le frottement est prédit à l’aide d’une simulation numérique dans des conditions de contact
rugueux avec une lubrification élastohydrodynamique (EHL).
La technique classique Multigrille fonctionne bien pour limiter le temps de calcul et les
besoins en mémoire. Cependant, le choix de la grille grossière a une influence importante
sur la robustesse du code et son efficacité pour résoudre le problème brut. Dans la première
partie de ce travail, une méthode de construction de grille grossière proposée par Alcouffe et al.
est implémenté dans le code EHL Multigrille indépendamment du temps. Ensuite ce solveur
modifié est étendu aux cas transitoires pour résoudre le problème de contact avec rugosité.
La courbe de frottement est généralement représentée en fonction du « Λ ratio », le rapport
entre l’épaisseur du film d’huile et la valeur moyenne quadratique de la rugosité de la surface.
Cependant, ce paramètre est moins approprié pour tracer les variations de frottement dans
des conditions de haute pression (régime élasto piézo-visqueux). Dans la deuxième partie de
ce travail, le coefficient de frottement est calculé à l’aide du code EHL modifié pour de nom-
breuses conditions de fonctionnement ainsi que pour les paramètres d’ondulation de surface.
Les résultats de la simulation montrent qu’il n’y a pas de courbe de frottement unique lorsque
l’ancien paramètre « Λ ratio » est utilisé. En se basant sur la théorie de la réduction d’amplitude,
un nouveau paramètre de dimensionnement qui dépend des conditions de fonctionnement et
des paramètres d’ondulation est trouvé, ce qui peut donner une courbe de frottement unique
pour les situations de haute pression.
Pour les surfaces rugueuses plus complexes, une méthode basée sur la densité spectrale de
puissance (PSD) est proposée pour prédire les variations de frottement dans la troisième partie
de ce travail. La rugosité artificielle de la surface est utilisée pour tester d’abord la méthode
de prédiction rapide. Un bon accord est trouvé entre la simulation numérique complète et
cette prédiction rapide. La méthode de prédiction rapide est ensuite appliquée pour analyser
la variation de frottement de la rugosité de surface mesurée. Une comparaison est également
faite entre les prédictions et les expériences.
Le nouveau paramètre d’échelle et l’augmentation du frottement prédite par la méthode
PSD montrent une bonne précision technique pour une utilisation pratique.

Mots clés: Lubrification élastohydrodynamique, Simulation numérique, Régime élasto piézo-

visqueux, Théorie de réduction d’amplitude, Variation de frottement

ii

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
Cette thèse est accessible à l'adresse : [Link]
© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Contents

1 Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.1 Methods to solve the rough contact problem . . . . . . . . . . . . . . . . . . 3
1.2.2 Friction in rough EHL contact problem . . . . . . . . . . . . . . . . . . . . . 6
1.3 Research aims and Outlines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3.1 Research aims . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3.2 Outlines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Numerical model 10
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 Transient EHL model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.2 Dimensionless equations and parameters . . . . . . . . . . . . . . . . . . . 11
2.3 The finite difference scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.4 Transfer operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4.1 Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4.2 Injection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.5 Coarse grid operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.6 Relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.7 Implementation of the Multi-Grid method . . . . . . . . . . . . . . . . . . . . . . . 22
2.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3 Friction influence of harmonic surface waviness 25

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 Lubricant rheological models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.1 Relative friction coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.2 Numerical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.4 Time-dependent solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.5 Effect of operating conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.6 Effect of surface anisotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.6.1 Longitudinal and transverse wavy cases . . . . . . . . . . . . . . . . . . . . . 38
3.6.2 Purely longitudinal wavy case . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4 Friction of complex rough surfaces 44

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.2 Power spectral density friction method . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.2.1 PSD friction model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.2.2 Model validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.3 The artificial surface roughness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.3.1 Surface roughness power spectrum . . . . . . . . . . . . . . . . . . . . . . . 47
4.3.2 Friction increase prediction of a rough surface . . . . . . . . . . . . . . . . . 48
4.3.3 Comparison between the EHL simulation and the PSD prediction . . . . . 50
4.4 Measured surface roughness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

iii

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CONTENTS

4.4.1 Friction prediction under a specific operating condition . . . . . . . . . . . 54

4.4.2 Operating condition effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.4.3 Friction curves for measured surface roughness . . . . . . . . . . . . . . . . 59
4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

5 Conclusion and perspective 64

5.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5.2 Perspective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

Appendix A Construction of the coarse grid operator 66

Appendix B Derivation of matrix A j for line relaxation 71

B.0.1 Gauss-Seidel line relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
B.0.2 Jacobi distributive line relaxation . . . . . . . . . . . . . . . . . . . . . . . . 72

Appendix C Derivation of the scaling parameter θ2 74

Appendix D The relation between the elastic deformation and corresponding pressure
for 2D wavy surfaces 75

Bibliography 78

iv

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 List of Figures

1.1 World primary energy consumption (red column: Non-OECD, blue column: OECD).
(Source: IEEJ Outlook 2019 and Scenario) . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 The variation of mean global surface temperature relative to 1880-2017. (Source:
NASA/GISS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Total consumption by End-Use Sector, from 2000 to 2017. (Source: Data from the
U.S. Energy Information Administration) . . . . . . . . . . . . . . . . . . . . . . . . 2
1.4 Pressure flow factors. (Source: Reference [21]) . . . . . . . . . . . . . . . . . . . . . 4
1.5 Relative amplitude as a function of ∇2 under pure rolling, where ∇2 is dimension-
less wavelength parameter, A i and A d are amplitude of surface roughness and
deformed surface roughness respectively. (Source: Reference [58] ) . . . . . . . . 5
1.6 Results obtained from measurements compared with theoretical attenuation curve
defined by [56]. (Source: Reference [61] ) . . . . . . . . . . . . . . . . . . . . . . . . 5
1.7 Friction coefficient versus speed for different loads. (Source: Ref. [83]) . . . . . . 7

2.1 Mesh point (x i , y j ) and it’s related mesh region r i , j . . . . . . . . . . . . . . . . . . 13

2.2 Mesh point (i , j ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3 Interpolation process (green points: coarse grid points, black dots: fine grid points,
blue dots: middle points on the fine grid, red point: central point on the fine grid). 17
2.4 Weighting factors for the interpolation (blue points: coarse grid points, black
dots: fine grid points). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.5 Flow chart of the hybrid relaxation process . . . . . . . . . . . . . . . . . . . . . . . 22
2.6 Implementation of the Multi-Grid method with a two level "V" cycle. . . . . . . . 23
2.7 The time-dependent "V" cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.1 Shear stress-shear rate relationship for the EHL contact. . . . . . . . . . . . . . . . 26

3.2 Comparison of the relative deformed amplitude (Ad /Ai ) as a function of f (r )∇2
for the current model (blue squares) with those on Reference [53] (solid line) . . 29
3.3 Top view of the surface waviness with λ/a h = 0.5 and Ai = 0.5H c: (a) the isotropic
surface waviness r = 1, (b) the longitudinal surface waviness r = 2, (c) the trans-
verse surface waviness r = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.4 Central line pressure P (X , 0) (black lines) and central line film thickness H (X , 0)
(blue lines) of isotropic surface waviness (r = 1) for M = 1000, L = 10, λ/a h = 0.5
and Ai = 0.5H c during a time period. The central pressure line (red line) for the
smooth case is plotted as a reference. . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.5 Central line pressure P (X , 0) (black lines) and central line film thickness H (X , 0)
(blue lines) of longitudinal surface waviness (r = 2) for M = 1000, L = 10, λ/a h =
0.5 and Ai = 0.5H c during a time period. The central pressure line (red line) for
the smooth case is plotted as a reference. . . . . . . . . . . . . . . . . . . . . . . . . 32
3.6 Central line pressure P (X , 0) (black lines) and central line film thickness H (X , 0)
(blue lines) of transverse surface waviness (r = 0.5) for M = 1000, L = 10, λ/a h =
0.5 and Ai = 0.5H c during a time period. The central pressure line (red line) for
the smooth case is plotted as a reference. . . . . . . . . . . . . . . . . . . . . . . . . 33

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 LIST OF FIGURES

3.7 The dimensionless central film thickness H c rough as a function of the dimension-
less time T for: M = 1000, L = 10, λ/a h = 0.5 and Ai = 0.5H c: (a) the isotropic
surface wavy case, (b) the longitudinal surface wavy case r = 2, (c) the transverse
surface wavy case r = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.8 The relative friction coefficient µr /µs as a function of the dimensionless time T
for: M = 1000, L = 10, λ/a h = 0.5 and Ai = 0.5H c: (a) isotropic surface wavy case,
(b) longitudinal surface wavy case r = 2, (c) transverse surface wavy case r = 0.5.
(Blue dotted line: average value of the relative friction coefficient.) . . . . . . . . . 34
3.9 Relative friction coefficient as a function of H c/Ai for a specific operating condition 35
3.10 Effect of the load parameter M on the relative friction coefficient for L = 10 and
λ/a h = 0.5: (a) relative friction coefficient as a function of H c/Ai , (b) relative fric-
tion coefficient as a function of M 0.33 · (H c/Ai ) . . . . . . . . . . . . . . . . . . . . 35
3.11 Effect of material parameter L on the relative friction coefficient for M = 2000
and λ/a h = 0.5: (a) relative friction coefficient as a function of H c/Ai , (b) relative
friction coefficient as a function of L −1.1 · (H c/Ai ) . . . . . . . . . . . . . . . . . . . 36
3.12 Effect of wavelength λ/a h on the relative friction coefficient for M = 1000 and
L = 10: (a) relative friction coefficient as a function of H c/Ai , (b) relative friction
coefficient as a function of (λ/a h )0.67 · (H c/Ai ) . . . . . . . . . . . . . . . . . . . . . 36
3.13 Relative friction coefficient as a function of the classical parameter "lambda ratio"
i.e. H c/Ai for a large range of operating conditions. . . . . . . . . . . . . . . . . . 37
3.14 Relative friction coefficient as a function of the new parameter θ2 , simulation re-
sults: black circles; fitted curve: the black dashed line. . . . . . . . . . . . . . . . . 38
3.15 Relative friction coefficient (µr /µs ) as a function of H c/Ai for different r (1 ≤ r ≤
32) values for: M = 1000, L = 10 and λ y /a h = 0.5 (left), zoom from 2.3 − 2.7 (right). 38
3.16 Relative friction coefficient (µr /µs ) as a function of ff(r ) · (H c/Ai ) for different r
(1 ≤ r ≤ 32) values for: M = 1000, L = 10 and λ y /a h = 0.5. . . . . . . . . . . . . . . . 39
3.17 Relative friction coefficient (µr /µs ) as a function of (H c/Ai ) for different r (0 ≤
r ≤ 1) values for: M = 1000, L = 10 and λx /a h = 0.5(upper), zoom from 2.4 − 2.6
(lower). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.18 Relative friction coefficient (µr /µs ) as a function of ff(r ) × (H c/Ai ) for different r
(0 ≤ r ≤ 1) values for: M = 1000, L = 10 and λx /a h = 0.5. . . . . . . . . . . . . . . . 40
3.19 ff (r ) as a function of r . Numerical results: red squares. Fitted curve: solid lines. . 41
3.20 Comparison between the transient relative friction coefficient and that of the sta-
tionary case. Transient results: black line. Stationary results: magenta dash-
dotted line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.21 Relative friction coefficients for the purely longitudinal wavy case: (a) relative fric-
tion coefficient as a function of original "lambda ratio" H c/Ai parameter, (b) rel-
ative friction coefficient as a function of the new parameter θ2∗ . . . . . . . . . . . 42

4.1 Flow chart for the relative friction coefficient prediction . . . . . . . . . . . . . . . 46

4.2 Surface roughness (a) and its 2D power spectral density (b) . . . . . . . . . . . . . 47
4.3 Power spectral density C i so of the self-affine surface(Figure 4.2(a)) with H = 0.8. . 48
4.4 The selected artificial surface roughness (a), amplitude distribution of this sur-
face roughness (b) and its power spectral density (c). . . . . . . . . . . . . . . . . . 48
4.5 The ratio of the deformed amplitude and the initial amplitude fitted as Equation
(3.17) (a) and the deformed surface roughness in frequency domain (b). . . . . . 49
4.6 Comparison between the initial surface roughness (a) and the deformed surface
roughness (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.7 Pressure increase distribution in frequency (a) and space (b) domains, respectively. 50

vi

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 LIST OF FIGURES

4.8 Shear stress distributions for the smooth case (a) and for the rough case (b). . . . 50
4.9 The generated surface roughness patch (a), the roughness patch in the high pres-
sure zone (b) and the periodical roughness pattern for full the numerical simula-
tion (c). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.10 Top view of the deformed surface roughness for a specific time step for a full nu-
merical simulation (a) and for a PSD prediction (b). Central line r r d (x, 0) of the
deformed surface roughness for the full numerical simulation (c) and for the PSD
prediction (d). Central line p(x, 0) of the pressure distribution for the full numer-
ical simulation (e) and for the PSD prediction (f). . . . . . . . . . . . . . . . . . . . 51
4.11 The relative friction as a function of dimensionless time employing the surface
roughness pattern in Figure 4.9 (c) for the full numerical simulation method. . . 52
4.12 Top view of the twenty generated artificial random rough surfaces from N◦ 1 to
N◦ 20 with a same standard divation value σ = 0.05µm and a same set of operating
conditions listed in Table 4.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.13 Measured surface roughness dART: (a) corrected surface roughness and (b) raw
surface roughness. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.14 Effective prediction areas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.15 An extracted surface patch of dART: (a) surface roughness height of this surface
patch; (b) deformed surface patch; (c) pressure fluctuation of the surface patch. . 55
4.16 Extracted 529 surface patch (left) and their relative friction coefficient values (right). 56
4.17 Corrected relative friction coefficient values for 528 surface patches (left) and its
histogram (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.18 Relative friction coefficient as a function of total number of surface patches under
the operating condition in Table 4.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.19 A extracted surface patch (left) and its average initial amplitude as a function as q. 57
4.20 The deformed surface patch shown in Figure 4.19 (left) and its pressure variations
for case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.21 The deformed surface patch shown in Figure 4.19 (left) and its pressure variations
for case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.22 The deformed surface patch shown in Figure 4.19 (left) and its pressure variations
for case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.23 The deformed surface patch shown in Figure 4.19 (left) and its pressure variations
for case 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.24 The deformed surface patch shown in Figure 4.19 (left) and its pressure variations
for case 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.25 Moes parameters M and L as a function of u r for the roughness dART. . . . . . . 59
4.26 The relative friction coefficient as a function of u r for the roughness dART. . . . . 60
4.27 The relative friction coefficient as a function of "Λ ratio" for the roughness dART. 60
4.28 Surface roughness dARTEb (left) and the top view of this roughness (right). . . . 60
4.29 Moes parameters M and L as a function of the rolling speed for the surface rough-
ness dARTEb. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.30 The relative friction coefficient as a function of u r for the surface roughness dAR-
TEb. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.31 The relative friction coefficient as a function of the "Λ ratio" for the surface rough-
ness dARTEb. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.32 The relative friction coefficient as a function of the rolling speed u r for the surface
roughness dART and dARTEb, respectively. . . . . . . . . . . . . . . . . . . . . . . . 62

vii

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 LIST OF FIGURES

A.1 Influences of the coarse grid operator L h on central point, east point and north-
east point. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
A.2 Influences of the injection operator J hH on nine coincidental points (blue points). 69

D.1 Pressure distribution and the corresponding elastic deformation. . . . . . . . . . 75

D.2 Amplitude of the elastic deformation A D d
as a function of initial pressure ampli-
p
tude A i for the following cases: (a) isotropic, (b) purely transverse, (c) purely lon-
gitudinal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
D.3 Amplitude of the elastic deformation A D d
as a function of wavelength λ (λ = min(λx , λ y ))
for the following cases: (a) isotropic, (b) purely transverse, (c) purely longitudinal. 76
D.4 Amplitude of the elastic deformation A D d
as a function of anisotropy parameter
r = λx /λ y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

viii

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 List of Tables

1.1 U.S. CO2 emissions from end-use sectors, 2008-2017. (Source: U.S. Energy Infor-
mation Administration, August 2018 Monthly Energy Review) . . . . . . . . . . . . 3

3.1 Relative friction coefficient versus the number of mesh points for: M =1000, L=10,
λx /a h = 0.5, r=0.4 and H c/Ai = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.2 Relative friction coefficient versus different surface anisotropy parameters for:
M = 1000, L = 10, λ y /a h = 0.5 and H c/Ai = 2. . . . . . . . . . . . . . . . . . . . . . 42

4.1 Operating condition parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.2 Relative friction coefficients as a function of the mesh points for two prediction
schemes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.3 The relative friction coefficient obtained by EHL simulation and PSD prediction
for 20 artificial random rough isotropic surfaces. . . . . . . . . . . . . . . . . . . . 53
4.4 Measured operating condition and lubricant parameters. . . . . . . . . . . . . . . 55
4.5 Operating conditions of selected cases. . . . . . . . . . . . . . . . . . . . . . . . . . 57

C.1 Range of the exponent for each parameter. . . . . . . . . . . . . . . . . . . . . . . . 74

ix

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
Cette thèse est accessible à l'adresse : [Link]
© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Nomenclature

Ad dimensionless deformed amplitude in the center of the contact

Ai dimensionless initial amplitude

p
3
ah the radius of the contact area a h = 3wR x /(2E 0 ) m

C 2D 2d power spectral density m4

C i so 2d power spectral density of an isotropic surface m4

E 1 , E 2 elastic moduli of the two contact bodies Pa

E0 reduced modulus of elasticity 2/E 0 = (1 − v 12 )/E 1 + (1 − v 22 )/E 2 Pa

f the friction force induced by the shearing of the lubricant N

F the dimensionless friction force

G dimensionless material parameter G = αE 0

Ge elastic shear modulus Pa

G∞ the limiting elastic shear modulus Pa

h film thickness m

H dimensionless film thickness H = hR x /a h2

h0 mutual approach m

H0 dimensionless mutual approach

hc central film thickness m

Hc dimensionless central film thickness for a smooth case H c = h c R x /a h2

H c rough dimensionless central film thickness for a rough case

hx, h y dimensionless mesh sizes on the fine grid

H x, H y dimensionless mesh sizes on the coarse grid

h
IH interpolation operator

I hH restriction operator

L x , L y lengths of final topography m

Ll coarse grid operator on the l th level

L dimensionless material parameter (Moes) L = G(2U )0.25

M 2d dimensionless load parameter (Moes) M = W2 (2U )−0.75

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 LIST OF TABLES

p pressure Pa

ps pressure for smooth cases Pa

ph Hertzian pressure Pa

δp pressure fluctuations Pa

∆P dimensionless pressure fluctuations ∆P = δp/p h

q x , q y wavenumbers in x and y directions respectively 1/m

qr roll-off wavenumber 1/m

Rq root mean square of the surface roughness m

Rx reduced radius of curvature in x : 1/R x = 1/R 1x + 1/R 2x m

Ry reduced radius of curvature in y : R y = R x m

r wavelength ratio used to describe the surface anisotropy r = λx /λ y

rr undeformed surface roughness m

rrd deformed surface roughness m

RR dimensionless surface roughness RR = r r · R x /a h2

SRR slide to roll ratio

t time s

T dimensionless time T = t ū/a h

∆T dimensionless time step

u 1 , u 2 velocities of lower surface and upper surface respectively m/s

ur mean velocity of contact surfaces u r = 0.5 × (u 1 + u 2 ) m/s

U dimensionless speed parameter U = (η 0 u r )/(E 0 R x )

∆U slide-to-roll ratio ∆U = δu/u r = (u 1 − u 2 )/u r

Urat slip parameter Urat = u 1 /u r

v 1 , v 2 Poisson ratios of the two contact bodies

w normal load N

W2 2d dimensionless load parameter W2 = w/(E 0 R x2 )

x coordinate in the rolling direction m

X dimensionless coordinate X = x/a h

y coordinate perpendicular to x m

Y dimensionless coordinate Y = y/a h

xi

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 LIST OF TABLES

z pressure viscosity index

α pressure viscosity index 1/Pa

ᾱ dimensionless viscosity index ᾱ = αp h

τ shear stress induced by the shearing of the lubricant Pa

τ0 the Eyring stress Pa

τL the limiting shear stress Pa

p
∇2 dimensionless wavelength parameter ∇2 = (λ/a h ) M /L

λ̄ dimensionless speed parameter

λt the time constant for the fluid

λx , λ y wavelength in x, y direction m

η viscosity Pa · s

η0 the atmospheric viscosity Pa · s

η̄ dimensionless viscosity η̄ = η/η 0

ρ density of the lubricant Kg · m−3

ρ0 atmospheric density of the lubricant Kg · m−3

ρ̄ dimensionless density ρ̄ = ρ/ρ 0

θ2 dimensionless new lambda ratio parameter

θ2∗ dimensionless new lambda ratio parameter for purely longitudinal rough surfaces

γ̇ strain rate 1/s

µ friction coefficient

σ standard deviation of surface roughness m

Abbreviations

ART Amplitude Reduction Theory

EHL Elastohydrodynamic Lubrication

FFT Fast Fourier Transform

PSD Power Spectral Density

Superscripts

d deformed

h the fine grid

H the coarse grid

xii

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 LIST OF TABLES

i so isotropic

l the l th level

Subscripts

a, b inlet, outlet

i,d initial, deformed

r, s rough, smooth

st start

x, y space domain

q x , q y frequency domain

xiii

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Chapter 1

Introduction

Contents
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.1 Methods to solve the rough contact problem . . . . . . . . . . . . . . . . . 3
1.2.2 Friction in rough EHL contact problem . . . . . . . . . . . . . . . . . . . . 6
1.3 Research aims and Outlines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3.1 Research aims . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3.2 Outlines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

1.1 Background
Energy is the most important of all resources, which is needed to support economic and so-
cial progress and build a better quality of life. According to the report of IEEJ Outlook 2019
(Institute of Energy Economics, Japan), the world primary energy consumption will continue
growing from 2018 to 2050. Most of this growth comes from non-OECD (non-Organization for
Economic Cooperation and Development) countries, where demand is driven by strong eco-
nomic growth (shown in Figure 1.1). Fossil energy like oil, coal and natural gas are still the
largest energy source of the world, and the reserves of fossil energy are limited. Meanwhile,
the over-consumption of fossil fuel leads to the over-release of carbon dioxide (CO2 ), methane,
oxynitride (NOx ) and particulates into the air, which disturb the natural balance of the atmo-
sphere. The rapid rise in the carbon dioxide contributes to the serious global warming problem.
NASA (National Aeronautics and Space Administration) reported that the global surface tem-
perature has been persistently increasing since the late 19th century (shown in Figure 1.2). An
investigation from the U.S. Energy Information Administration (EIA) shows that the industrial
and transportation sectors have consumed the most energy (shown in Figure 1.3) as well as
produced the most CO2 emissions (show as Table 1.1).
However, a substantial amount of energy is not put to useful purposes. Researches [1–4]
show that a considerable amount of energy in industrial and transportation is consumed to
overcome friction. For instance, energy consumed to overcome friction over the total energy
consumption in heavy-duty vehicles is 33%, in paper machines is 32%, in passenger cars is
33% and in mineral mining industry is 40%. Recently, the increasing environment awareness
requires efforts to improve energy efficiency and reduce CO2 production. Correct lubrication
between engineering provides sufficient separation of the roughness present on the contact
surfaces, which contributes to reducing friction losses. Better understanding and control of
friction in mechanical components has the potential to offset large energy savings and CO2
emission reduction [5]. Studies [6–8] estimated that with the implementing advanced tribolog-
ical technologies, such as using new contact surface, materials, lubricants, energy losses due to
friction and wear could potentially be reduced by 40% in the long term (15 years) and by 18% in
the short term (8 years )and CO2 emissions globally can also reduced by 1,460 MtCO2 (million
tonnes CO2 ) in the short term and by 3,140 MtCO2 in the long term.

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 1. INTRODUCTION

Figure 1.1: World primary energy consumption (red column: Non-OECD, blue column: OECD).
(Source: IEEJ Outlook 2019 and Scenario)

Figure 1.2: The variation of mean global surface temperature relative to 1880-2017. (Source:
NASA/GISS)

Figure 1.3: Total consumption by End-Use Sector, from 2000 to 2017. (Source: Data from the
U.S. Energy Information Administration)

Hence, understanding the mechanisms of friction and improving frictional behavior be-
tween engineering contact surfaces still remains an important issue in today’s research, not

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 1. INTRODUCTION

only for meeting the increasing requirements of energy efficiency and CO2 emissions reduc-
tion but also for providing a theoretical tool in element design and optimization.

Table 1.1: U.S. CO2 emissions from end-use sectors, 2008-2017. (Source: U.S. Energy Informa-
tion Administration, August 2018 Monthly Energy Review)
Units: million metric tons
year
end-use sector 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017
transportation 1898 1832 1848 1817 1779 1805 1823 1848 1886 1902
industrial 1608 1400 1508 1498 1489 1508 1511 1456 1428 1409
residential 1234 1127 1210 1149 1043 1100 1115 1037 982 956
commercial 1075 1007 1025 990 932 958 970 932 894 875

1.2 Literature review

Elastohydrodynamic lubrication (EHL) is the type of lubrication for frictional pairs having elas-
tic contact under very high pressure and forming lubricant film in non-conformal contacts,
such as rolling bearings, gears, human synovial joints and so on [9]. Lubricant film and sur-
face roughness play an important role for improving reliability and effectiveness of mechanical
parts as well as reducing friction losses. The majority of the published work on the influence
of surface roughness on friction has been experimental, the minority of theoretical work has
been done on friction prediction. This section represents the literature review on the methods
to solve the rough EHL contact problem as well as friction in rough EHL contacts.

1.2.1 Methods to solve the rough contact problem

Typically, the EHL model consists of five equations [10], in which the Reynolds equation Equa-
tion 1.1 is a partial differential equation and the film thickness equation Equation 1.2 contains
an integral term, both equations are required to be solved simultaneously, making these equa-
tions very complex. There are many approaches that can be used to solve this EHL model: the
inverse method [11], the Newton-Raphson method [12], the homotopy method [13], the finite
element method [14], the Multigrid method [15,16] and the Navier-Stokes approach [17]. So far,
the Multigrid algorithm has been considered as one of the most efficient methods and applied
frequently to EHL problems.
∂ ρh 3 ∂p ∂ ρh 3 ∂p ∂(ρh)
( )+ ( ) − ur =0 (1.1)
∂x 12η ∂x ∂y 12η ∂y ∂x
x2 y 2 2 +∞ +∞ P (x 0 , y 0 )
Z Z
h(x, y) = h 0 + + − r r (x, y) + 2 d x 0d y 0 (1.2)
π −∞ −∞
p
2 2 0 2
(x − x ) + (y − y )0 2

Where p represents the pressure, h denotes the film thickness, h 0 is the mutual approach and
u r = (u 2 + u 2 )/2 is the mean velocity of two contact surfaces. ρ and η are the density and vis-
cosity of the lubricant, respectively. The x axis is aligned with the direction of the mean velocity
ū, and the y axis is perpendicular to the x direction.
In engineering, no surface is perfectly smooth, the order of magnitude of the surface rough-
ness is often the same as or greater than that of the film thickness predicted by smooth contact
conditions [18]. Therefore, the surface roughness should be considered. Generally, there are
two methods to treat the rough lubrication problem numerically.

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 1. INTRODUCTION

One approach is called the "stochastic" method. Early work was conducted on the hydrody-
namic lubrication (HL) problem. Theoretical analysis of the implementation of the stochastic
theory on rough HL contact problem was described by Christensen [19, 20]. Then Patir and
Cheng [21, 22] proposed an average flow model to determine the effects of surface roughness
on rough-lubricated contacts. In this model the Reynolds equation is simplified as an average
Reynolds equation using a independent flow factor (shown in Figure 1.4). Since the pioneering
studies on the "average flow model", a number of authors have extended and generalized this
work. Tripp [23] re-computed the flow factors using a perturbation expansion of the pressure in
a nominal parallel film. When small roughness amplitude is considered, the flow factors calcu-
lated in Reference [23] agree well with that of Patir and Cheng [21]. Hu and Zheng [24] studied
the influence of boundary conditions, grid systems and statistics of rough surfaces on the flow
factors. Lunde and Tonder [25] calculated the flow factors for an isotropic rough bearing and
found that the boundary conditions of the selected bearing part can not affect the flow pass-
ing through. Subsequently, Zhu and Cheng [26] extended the flow factors method in the point
EHL contact problem. Some authors [27,28] applied the flow factors to deal with the cavitation
problem. Letalleur et al. [29] studied the flow factors for two rough cases: smooth-rough sta-
tionary case and rough-rough unstationary case. Sahlin et al. [30] developed a novel method
using a homogenization technique to compute the flow factors. However, in stochastic model,
roughness asperities are mainly treated as rigid.

Figure 1.4: Pressure flow factors. (Source: Reference [21])

Another way is to incorporate the surface roughness term r r (x, y, t ) in the film thickness
equation (shown as Equation 1.2) and to solve the system of equations directly. Due to the
limitation of computation of speed and storage space, preliminary research [31–35] studied
the steady state line rough contact problem, where the surface roughness is time-independent
and one-dimensional model was considered. Later on, the stationary two dimensional rough
contact problem [36–39] were carried out. With the increasing development of computational
technique, transient cases were studied by many authors [40–51]. Based on the previous stud-
ies on transient rough contact problem, Venner and Lubrecht et al. [52–60] published a series
of papers on the "Amplitude Reduction Theory" describing the relation between the surface
roughness deformation and operating conditions. They found that under very high pressure
situations (piezoviscous elastic regime), the surface roughness will deform, and this deforma-

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 1. INTRODUCTION

tion depends on operating conditions. A master curve (shown in Figure 1.5) describes this
relation quantitatively. Then Šperka et al. [61] verified the "Amplitude Reduction Theory" by
measuring the deformed surface roughness on an optical test rig, the comparison of the mea-
sured results and the predicted results can be found in Figure 1.6. Recently, some extension
work on the rough contact problem were addressed [62–65]. One of the most important studies
about the influence of surface roughness on friction will be represented in the next subsection.

Figure 1.5: Relative amplitude as a function of ∇2 under pure rolling, where ∇2 is dimensionless
wavelength parameter, A i and A d are amplitude of surface roughness and deformed surface
roughness respectively. (Source: Reference [58] )

Figure 1.6: Results obtained from measurements compared with theoretical attenuation curve
defined by [56]. (Source: Reference [61] )

The surface roughness term r r (x, y), which is often considered to be the same order of mag-
nitude as the oil film thickness, incorporated in Equation 1.2 makes the coefficient (ρh 3 )/(12η)
in the Reynolds equation jump orders of magnitude, which leads to a significant variation in
EHL equations. From a mathematical point of view, the coefficient (ρh 3 )/(12η) is continu-
ous, while in the numerical simulation, the coefficient causes a strong discontinuity in discrete

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 1. INTRODUCTION

Reynolds equation. Work by Alcouffe et al [66] proposed an efficient way to overcome the above
discontinuity through constructing the coarse grid in a Multi-grid code. The present work em-
ploys the Multigrid techniques [10] to solve the rough contact problem. In addition, the idea
provided by Alcouffe et al [66] is also applied to improve the code robustness and code effi-
ciency.

1.2.2 Friction in rough EHL contact problem

Most heavily loaded machine elements are working under elastohydrodynamically lubricated
conditions. Understanding the frictional behavior in such contacts plays an important role for
reducing friction, preventing wear as well as improving service life. To reveal the relation be-
tween fluid properties and friction, research has been conducted as follows: Crook [67] used a
disc machine for measuring frictional traction. It was shown that the rolling friction (the trac-
tion due to rolling) is independent of load and simply proportional to the film thickness in the
elasto-hydrodynamic regime. Johnson and Cameron [68] measured traction in a rolling contact
disc machine and results showed that the traction first increases and then decreases when the
sliding speed increases. Johnson and Roberts [69] observed the visco-elastic behavior of film
thickness through measuring shear forces on a rolling-contact test rig. Evans and Johnson [70]
constructed traction maps depending upon pressure, temperature as well as shear rate for dif-
ferent fluids, where different areas represent different traction behavior. Zhang et al. [71] stud-
ied the elliptical contact between rib face and roller end in tapered roller bearings by means
of a full numerical simulation. They found that the elastic deformation has a non-negligible
influence on the friction coefficient. Yu and Medley [72] studied the influence of lubricant ad-
ditives on friction via a side-slip disc machine. They concluded that the limiting shear stress,
which is a useful parameter for predicting friction, is affected by the lubricant additives. Jacod
et al. [73] predicted the coefficient of friction over a wide range of operating conditions and
obtained a single generalized friction curve based on a full numerical simulation for a non-
Newtonian EHL contact model. Vicente et al. [74] explored friction in rolling-sliding, soft-EHL
contacts numerically and experimentally. Numerical calculations of the Couette friction are in
good agreement with measured results. Very recently, Liu et al [75] calculated the friction co-
efficient in a gear contact interface numerically, based on a thermal starved EHL model. They
found that the maximum friction coefficient appears at the engaging-in point where a consid-
erable slide-to-roll ratio exists. Björling et al. [76] measured the friction under EHL conditions
on a ball-on-disc test rig for aged and fresh oils. Results showed that there is no difference in
friction. In addition, Zhang [77] measured the EHL friction for a wide range of base fluids and
compared the friction values for five different operating conditions. The study underlined the
importance of molecular structure of the base fluid in determining the EHL friction.
Studies [78–80] showed that surface roughness has a significant impact on the friction be-
havior of lubricated surfaces. A useful tool to investigate the frictional behavior between rough
surfaces is the classical Stribeck curve, showing the friction coefficient is a function of a ratio of
the averaged oil film thickness to the combined surface roughness. The original research about
the Stribeck curve dates back to the 19th century. In 1879, Thurston gave precise values of the
friction coefficient and he was probably the first person to report that the friction coefficient
passed through a minimum as the load increased [81, 82]. Twenty years later, Stribeck [83, 84]
systematically published results of a carefully conducted and wide-ranging series of experi-
ments on journal bearings, which are frequently referred to as ‘the Stribeck curve’ (shown in
Figure 1.7). Gümbel [85] organised Stribeck’s experimental results in a single curve by plotting
the friction against the parameter ηω/p̄, where η is the lubricant viscosity, ω is the angular ve-
locity of the shaft and p̄ is the load per unit length. At the same time, Hersey [86] conducted

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 1. INTRODUCTION

experiments on journal bearings and plotted the friction coefficient against the load, speed,
temperature, viscosity and rate of oil supply. He showed that hydrodynamic friction should be
a function of ηn/p in which n is the rotational speed and p is the pressure. Many years later,
Wilson and Barnard [87] replotted the Stribeck curve by introducing a new variable i.e. zn/p,
where the lower-case z stands for the lubricant viscosity. Subsequently, McKee [88] provided a
similar dimensionless group Z N /P . Vogelpohl et al. [89] incorporated the boundary and fluid
friction coefficient and showed a transition from the hydrodynamic lubrication regime to the
mixed lubrication regime. All of the work mentioned above is performed under low pressure
conditions, in the isoviscous rigid regime [90].

Figure 1.7: Friction coefficient versus speed for different loads. (Source: Ref. [83])

The situation for non-conforming contacts, such as those occurring in rolling element bear-
ings, gears and cams, is somewhat different [91]. Shotter [92] experimentally showed that the
friction increases with the surface roughness. Tallian and his co-workers [93, 94] proposed a
ratio |ξ0 | between the elastohydrodynamic film thickness and the composite root mean square
roughness to represent the mixed elastodydrodynamic regime ( 1 < |ξ0 | < 4 ). Poon [95] was
concerned with the transition from the boundary to the mixed regime with a dimensionless
parameter 1 É ξ É 2 and the transition from mixed to full EHL region with 2 < ξ É 2.4 by using
electrical-conductivity measurements. Bair and Winer [96] plotted the reduced traction coef-
ficient as a function of a lambda ratio by performing sliding-rolling experiments. They found
that when the lambda ratios is less than 2 the contact moves into the mixed regime. In gen-
eral, the Stribeck curve can be divided into three regimes [97]: λ > 3 represents the full-film
regime, 1 É λ É 3 is the mixed EHL regime and λ < 1 indicates the boundary regime. However,
study [98] shows that this lambda ratio is not a suitable parameter to determine lubrication
states especially when some aspects such as non-Newtonian, thermal and transient effects are
considered. Transition locations from mixed to boundary lubrication regime or from full-film
to mixed lubrication regime are still ambiguous. Therefore, an appropriate grouping including
the speed, film thickness and roughness is required. Schipper [99] suggested a so-called Lu-
brication number L 0 , which takes viscosity, speed and pressure into consideration, to detect
the variation of the friction coefficient. Recently, Gelinck [100] extended Johnson’s model [101]
to calculate the coefficient of friction for the whole mixed EHL regime. Lu and Khonsari [102]

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 1. INTRODUCTION

examined the behavior of the Stribeck curve theoretically and experimentally on a journal bear-
ing and found a good agreement. Wang et al. [103] presented a numerical approach developed
on the basis of deterministic solutions of mixed lubrication to evaluate sliding friction. Mean-
while, they measured the sliding friction on a commercial test rig. Both results were plotted
against sliding velocities and also showed good agreement. Kalin [104] investigated changes
of the Stribeck curve when one or two surfaces in the contact are non-fully wetted. After-
wards, Kalin [105] tested the variations of the friction coefficient with diamond-like carbon
coatings (DLC). Zhang [106] developed a numerical approach assuming the asperity interac-
tion friction is proportional to the contact area to predict the mixed EHL friction coefficient.
Bonaventure [107] and his co-authors conducted rolling-sliding experiments with random sur-
face roughness, they found that the onset of ML occurs at a higher entrainment product η 0 u e
(in which η 0 is inlet viscosity and u e is entrainment speed) and a relevant roughness scalar
parameter was obtained to predict the onset position.
Most of the work on Stribeck curve was done by experiments. Current study employs the
Amplitude Reduction Theory [53] to study the frictional behavior in piezoviscous elastic regime
[108] by means of numerical simulation.

1.3 Research aims and Outlines

1.3.1 Research aims

Long term successful operation of mechanical devices greatly depends on correct lubrication
of the mechanical elements to provide sufficient separation of the roughness present on the
contact surfaces. However, lubrication provides another important role, reducing friction be-
tween rough contact surfaces.
The objective of the present research project is to develop an efficient and robust Multi-
Grid-based algorithm to study the frictional behavior between rough contact surfaces. Current
MultiGrid codes show the required efficiency, but are not sufficiently robust to treat the rough
surface problem in a general way. Difficulties may lie in the following aspects:

(i) The efficient construction of the coarse grid of EHL Multi-Grid model to guarantee the
code robustness and code efficiency of impact, rough surface EHL contact problems.

(ii) Tests of the increased robustness of the new EHL Multi-Grid solver.

(iii) Implementation applied to test the code robustness and code efficiency of rough surface
EHL contact problems.

(iv) The extension of the developed Multi-Grid lubrication code to transient contact prob-
lems.

(v) The computation of the friction coefficient of rough contact surfaces.

(vi) The unification of friction curves that differ according to operating conditions.

(vii) The extension of the lambda ratio parameter predicting the transition from mixed to full-
film regimes.

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 1. INTRODUCTION

1.3.2 Outlines
According to the research aims listed in the previous sub-section. The layout of this thesis is as
follows:
Chapter 1: This chapter first emphasizes the important role of friction played in energy
consumption and environmental issues. Subsequently, a literature review on the methods to
solve rough contact problems and friction in rough contact problems is illuminated. The ob-
jective and structure of the present thesis are given in the last section.
Chapter 2: This chapter represents the numerical model and algorithm for solving the tran-
sient rough EHL contact problem. The governing equations for transient EHL model are intro-
duced first. Then the method proposed by Alcouffe et al [66] is employed to construct transfer
operators as well as coarse grid operator. Finally the Multi-Grid method [10] is implemented.
Chapter 3: In this chapter, lubricant rheological models are illustrated in the first place. The
relative friction coefficient, an indicator for the full-film-mixed lubrication regime transition,
is proposed in methodology section. Then the relative friction coefficient is calculated numer-
ically for isotropic as well as anisotropic harmonic surface roughness respectively. Finally, a
single friction curve is obtained using a new "lambda ratio" parameter.
Chapter 4: A rapid analytical prediction method using the power spectral density [109] is
proposed to study a more complex surface topography in this chapter, firstly. Then an artificial
surface roughness is employed to test this rapid prediction method. Finally, the prediction
method is applied to predict friction for measured rough surfaces.
Chapter 5: The main results of current work are summarized and recommendations of
future work are also made.

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Chapter 2

Numerical model

Contents
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 Transient EHL model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.2 Dimensionless equations and parameters . . . . . . . . . . . . . . . . . . 11
2.3 The finite difference scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.4 Transfer operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4.1 Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4.2 Injection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.5 Coarse grid operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.6 Relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.7 Implementation of the Multi-Grid method . . . . . . . . . . . . . . . . . . . . . 22
2.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.1 Introduction
Multi-grid methods have been used successfully to treat Elastohydrodynamic lubrication (EHL)
problems in the past [110–112]. However, when taking the surface roughness into account, film
thickness and viscosity jump violently, both of them are strongly discontinuous parameters in
discrete equations and will influence code robustness and code efficiency. The paper by R.
Alcouffe [66] proposed an efficient way to solve this problem through constructing the coarse
grid in a Multi-grid code. In this chapter, the Multigrid method is applied to solve the transient
EHL model, and the algorithm outlined in Reference [66] is also implemented.

2.2 Transient EHL model

2.2.1 Governing equations

The lubrication of rough surfaces in EHL contacts is inherently a highly transient process.
Study [40] shows that the surface roughness induced by the transient effect has a remarkable
influence on the pressure and film thickness profiles. For the time-dependent problem [10],
the Reynolds equation is given as:

∂ ρh 3 ∂p ∂ ρh 3 ∂p ∂(ρh) ∂(ρh)
( )+ ( ) − ur − =0 (2.1)
∂x 12η ∂x ∂y 12η ∂y ∂x ∂t
| {z } | {z } | {z }
poiseuille couette transient

with p = 0 on the boundaries and the cavitation condition p Ê 0 everywhere. Where p is the
pressure, h is the film thickness whose expression is shown as Equation 2.2 and u r = (u 1 +u 2 )/2
is the mean velocity (u 1 and u 2 are the velocities of two contact surfaces respectively). The

10

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 2. NUMERICAL MODEL

direction of the x axis is as same as that of the mean velocity u r , the y axis is perpendicular to
x axis and t is time.
The equation used to describe the gap between the two contact bodies is the film thickness
equation:

x2 y2 2 +∞ Z +∞ p(x 0 , y 0 , t )
Z
h(x, y, t ) = h 0 (t ) + + − r r (x, y, t ) + d x 0 d y 0 (2.2)
πE 0
p
2R x 2R y −∞ −∞ (x − x 0 )2 + (y − y 0 )2
| {z }
elastic deformation

in which r r (x, y, t ) stands for surface roughness. R x and R y represent the reduced radius of
curvature in x and y direction respectively. h 0 denotes the rigid body approach. E 0 is called the
reduced elastic modulus and its expression can be found below. The elastic deformation term
is calculated with the approach named multilevel multi-integration [10, 113].

2 1 − v 12 1 − v 22
= +
E0 E1 E2

and E 1 and E 2 are the elastic moduli of the two contact bodies. v 1 and u 2 are the Poisson ratios.
In order to have a load balance. The integral of the pressure distribution should be equal to
the applied load w. Z Z +∞ +∞
p(x 0 , y 0 , t )d x 0 d y 0 = w(t ) (2.3)
−∞ −∞
In the Reynolds equation (2.1), ρ is the density and η is the viscosity of the lubricant. Both
of them are functions of pressure. A simply density pressure relation is given by Dowson and
Higginson [114]:
5.9 × 108 + 1.34p
ρ(p) = ρ 0 (2.4)
5.9 × 108 + p
where ρ 0 is the atmospheric density. The simplest viscosity pressure relation is proposed by
Barus [115]:
η(p) = η 0 exp(αp) (2.5)
in which η 0 is the atmospheric viscosity and α is the pressure viscosity coefficient. However,
this exponential Barus equation usually predicts a higher viscosity value when the pressure is
very large. A more realistic relation is derived by Roelands [116]:
p z
η(p) = η 0 exp[(ln(η 0 ) + 9.67)(−1 + (1 + ) )] (2.6)
p0

where η 0 is the atmospheric viscosity and z is the pressure viscosity index, typically z = 0.6 and
p 0 = 1.98 × 108 Pa.

2.2.2 Dimensionless equations and parameters

To simplify the equation system and generalize the EHL model, the equations described above
are made dimensionless using dimensionless variables based on the Hertzian dry contact solu-
tion [117]. For the dry point contact case, the pressure distribution profile required for contact
deformation reads:

p h 1 − (x/a h )2 − (y/a h )2 if x 2 + y 2 ≤ a h2
( p
p(x, y) = (2.7)
0 otherwise

11

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 2. NUMERICAL MODEL

with a h the radius of the contact area:

r
3 3wR x
ah = (2.8)
2E 0
and p h is referred to as the Hertzian pressure:

3w
ph = . (2.9)
2πa h2

Then the dimensionless variables are introduced to simplify the EHL model:

X = x/a h Y = y/a h
P = p/p h H = hR x /a h2
(2.10)
η̄ = η/η 0 ρ̄ = ρ/ρ 0
T = u r t /a h ᾱ = αp h

Substituting the dimensionless variables in Equation 2.1 yields:

∂ ρ̄H 3 ∂P ∂ ρ̄H 3 ∂P ∂(ρ̄H ) ∂(ρ̄H )

( )+ ( )− − =0 (2.11)
∂X η̄λ̄ ∂X ∂Y η̄λ̄ ∂Y ∂X ∂T

with X ∈ [X a , X b ] and Y ∈ [−Y a , Y a ]. Where λ̄ = (12u r η 0 R x2 )/(a 3 p h ). And the boundary condi-
tions are P (X a , Y a ) = P (X a , Y ) = P (X b , Y ) = P (X , Y a ) = P (X , −Y a ) = 0. The cavication condition
is P (X , Y , T ) ≥ 0.
The dimensionless film thickness equation becomes:

X2 Y 2 2 P (X 0 , Y 0 , T )
+∞ Z +∞
Z
H (X , Y , T ) = H0 (T ) + + − RR(X , Y , T ) + 2 d X 0d Y 0
π
p
2 2 −∞ −∞ 0 2
(X − X ) + (Y − Y ) 0 2
(2.12)
where H0 (T ) is determined by the dimensionless force balance equation:
Z +∞ Z +∞
2π
P (X 0 , Y 0 , T )d X 0 d Y 0 = (2.13)
−∞ −∞ 3

The dimensionless density equation for a compressible lubricant reads:

5.9 × 108 + 1.34p h P

ρ̄(P ) = (2.14)
5.9 × 108 + p h P

The dimensionless forms of viscosity equations are:

Barus: η̄ = exp(ᾱP ) (2.15)

and
ph z
Roelands: η̄ = exp((ln(η 0 ) + 9.67)(−1 + (1 + P ) )). (2.16)
p0
Beside the dimensionless variables mentioned in Equation 2.10, two dimensionless numbers
are often used to reduce the number of parameters, they are referred as Moes dimensionless
parameters [118, 119]. For point contact they are defined as [119]:

w 2η 0 u r −3/4
M= ( ) (2.17)
E 0 R x2 E 0Rx

12

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 2. NUMERICAL MODEL

and
2η 0 u r 1/4
L = αE 0 ( ) . (2.18)
E 0Rx
For convenience, the Moes parameters can be used to re-write the parameters ᾱ and λ̄:

128π3 1/3
λ̄ = ( )
3M 4
(2.19)
L 3M 1/3
ᾱ = ( )
π 2
Hamrock and Dowson [120] introduced three parameters to simplify the study of film thick-
ness. For point contact they are written as:
w
W=
E 0 R x2
η0 ur (2.20)
U=
E 0Rx

G = αE 0

2.3 The finite difference scheme

The second-order self-adjoint elliptic partial differential equation considered by Alcouffe [66]
is
− ∇ · (D(x, y, t )∇U (x, y, t )) + σ(x, y, t )U (x, y, t ) = f (x, y, t ) (x, y) ∈ Ω (2.21)
Compared to this equation, the Reynolds equation is of the same type with U = P , σ = 0, D =
−(ρ̄H 3 )/(η̄λ̄) and f = ∂(ρ̄H )/∂X + ∂(ρ̄H )/∂T . Rearranging the Reynolds equation (2.11) yields:

∂ ∂P ∂ ∂P ∂(ρ̄H ) ∂(ρ̄H )
− (D )− (D )= + (X , Y ) ∈ Ω (2.22)
∂X ∂X ∂Y ∂Y ∂X ∂T
In the present work, the calculation domain Ω is a rectangle [X a , X b ] × [−Y a , Y a ]. This do-
main is covered with a uniform grids with a system of straight lines parallel to the coordinate
axes. The mesh size in the two directions is hx = (X b − X a )/N x and h y = 2×Y a /N y , in which N x
and N y are the number of mesh points in both directions.

Figure 2.1: Mesh point (x i , y j ) and it’s related mesh region r i , j

13

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 2. NUMERICAL MODEL

To derive the difference scheme, a mesh region r i , j (shown as Figure 2.1) defined by the
lines x = x i − h x /2 , x = x i + h x /2 , y = y j − h y /2 and y = y j + h y /2 for each mesh point (x i , y j ) is
selected. In terms of each mesh point (x i , y j ), P (x i , y j , t k ) := P i , j ,k is unknown, now integrating
Equation 2.22 over the corresponding mesh region r i , j :

∂ ∂P ∂ ∂P ∂(ρ̄H ) ∂(ρ̄H )
Z Z Z Z
− [ (D )+ (D )]d xd y = [ + ]d xd y (2.23)
ri , j ∂X ∂X ∂Y ∂Y ri , j ∂X ∂T

According to Green’s Theorem [121], Equation 2.23 can be expressed as:

∂P ∂P ∂(ρ̄H ) ∂(ρ̄H )
Z Z Z
− [(D )d y − (D )d x] = [ + ]d xd y (2.24)
ci , j ∂X ∂Y ri , j ∂X ∂T

where c i , j is the boundary of r i , j and the integration path along this boundary is anticlockwise.
Supposing f (x i , y j , t k ) := f i , j ,k , the double integrals of the right hand side of Equation 2.24 can
be simply approximated by means of
Z Z
.
f (x, y, t )d xd y = f i , j ,k · a i , j . (2.25)
ri , j

where a i , j = hx· h y is the area of the rectangle region r i , j shown in Figure 2.2.

Figure 2.2: Mesh point (i , j ).

Referring again to Figure 2.2, the line integal of Equation 2.24 over the four boundaries of
ri , j is approximated by means of central differences as:

∂P ∂P P i , j ,k − P i +1, j ,k
Z
.
− [(D )d y − (D )d x] = (h y)[D i +1/2, j ,k ( )
ci , j ∂X ∂Y hx
P i , j ,k − P i −1, j ,k
+ D i −1/2, j ,k ( )]
hx
(2.26)
P i , j ,k − P i , j +1,k
+ (hx)[D i , j +1/2,k ( )
hy
P i , j ,k − P i , j −1,k
+ D i , j −1/2,k ( )]
hy

14

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 2. NUMERICAL MODEL

Rewriting Equation 2.26 and combining Equation 2.25 gives:

A i , j ,k (P i , j +1,k − P i , j ,k )+A i , j −1,k (P i , j −1,k − P i , j ,k )+

(2.27)
B i , j ,k (P i +1, j ,k − P i , j ,k ) + B i , j −1,k (P i −1, j ,k − P i , j ,k ) = F i , j ,k

where
1 hx
A i , j ,k = −( )( )(D i , j ,k + D i , j +1,k )
2 hy

1 hy
B i , j ,k = −( )( )(D i , j ,k + D i +1, j ,k )
2 hx

F i , j ,k = (hx · h y) f i , j ,k
In terms of f i , j ,k , the same discrete schemes used in Reference [10] is adopted. At this point,
the right hand side of Equation 2.27 can be taken as:
.
F i , j ,k = h y(1.5ρ̄ i , j ,k Hi , j ,k − 2ρ̄ i −1, j ,k Hi −1, j ,k + 0.5ρ̄ i −2, j ,k Hi −2, j ,k )
hx · h y (2.28)
+ (1.5ρ̄ i , j ,k Hi , j ,k − 2ρ̄ i , j ,k−1 Hi , j ,k−1 + 0.5ρ̄ i , j ,k−2 Hi , j ,k−2 )
ht
where ht is the mesh size in time domain. A more detailed derivation of the above difference
scheme can be found in Reference [122].

2.4 Transfer operators

Intergrid transfers are used for connecting the fine grid with the coarse grid. After a number
of relaxations the error on the fine grid is smooth enough to be approximate on the coarse
h
grid. Hence a restriction operator I H is needed to transfer the approximated solution Pe and
h
the residual r h . When the low frequency errors have been eliminated on the coarse grid, it
h
is necessary to define a new error υh (υh = P h − Pe ) on the fine grid to correct the fine grid
h
approximate solution Pe . The classical bi-linear interpolation works quite well for most load
cases. However when D jumps by orders of magnitude, Alcouffe [66] proposed a more efficient
interpolation operator and this type of operator allows D∇P to be continuous over the whole
calculation domain and gives a more reasonable physical representation on the coarse grid
[123].

2.4.1 Interpolation
Having defined the coefficients A and B in Equation 2.27, it is time to define the interpolation
operator. In matrix form, the interpolation can be represented as:

υh = I H
h H
υ (2.29)

where υh and υH are the fine grid and coarse grid error vectors respectively. I Hh
is the interpo-
lation operator and the superscripts h and H stand for the fine grid and the coarse grid respec-
tively. The new coarse grid construction method proposed by Alcouffe et al. [66] is used here,
the interpolation process will be illustrated as follows:
The first step is to interpolate the fine grid points (black points shown in Figure 2.3 (b)) coin-
ciding with the coarse grid points (green points shown in Figure 2.3 (a)):

υhiF, j F,k = υiC

H
, jC ,k (2.30)

15

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 2. NUMERICAL MODEL

where υh is the error on the fine grid, and υH is the error on the coarse grid. Subscripts (i F, j F, k)
and (iC , jC , k) are applied for illustrate the mesh points on the fine grid and on the coarse grid
at the k th time step respectively.
The second step is to obtain the middle points, represented as blue dots in Figure 2.3 (c), on
the fine grid. Along horizontal lines, the expression for middle points is:

(B ihF, j F,k υiC

H
, jC ,k
+ B ihF +1, j F,k υiC
H
+1, jC ,k
)
υhiF +1, j F,k = (2.31)
(B ihF, j F,k + B ihF +1, j F,k )

A similar expression can be derived for vertical lines, which reads:

(A hiF, j F,k υiC

H
, jC ,k
+ A hiF, j F +1,k υiC
H
+1, jC ,k
)
υhiF, j F +1,k = (2.32)
(A hiF, j F,k + A hiF, j F +1,k )

Finally, the central point represented as a red point in Figure 2.3 (d) on the fine grid, which is
obtained as:

υhiF +1, j F +1,k = (A hiF +1, j F +1,k υhiF +1, j F +2,k + A hiF +1, j F,k υhiF +1, j F,k +
B ihF, j F +1,k υhiF, j F +1,k + B ihF +1, j F +1,k υhiF +1, j F +1,k )/ (2.33)
(A hiF +1, j F +1,k + A hiF +1, j F,k + B ihF, j F +1,k + B ihF +1, j F +1,k )

The above pointwise description (from Equation 2.30 to Equation 2.33) can be replaced by
the matrix expression Equation 2.29, in which the matrix is large and complex. A simply way to
describe this matrix is by using a stencil notation. As was shown in Figure 2.4, in the interpo-
lation process, the stencil provides weighting factors for dividing the coarse grid value in point
(iC , jC , k) to the coinciding fine grid point (i F, j F, k) as well as its 8 adjacent points. Observing
those pointwise expressions, the contribution of the coarse grid point to the 9 corresponding
h
fine grids can be written as a stencil I H in Equation 2.34.

NWihF, j F,k NihF, j F,k N E ihF, j F,k

 
 
 
h
h C ihF, j F,k E ihF, j F,k
 
IH  Wi F, j F,k
= (2.34)


 
 
SWihF, j F,k S ihF, j F,k SE ihF, j F,k

where

C ihF, j F,k = 1,

A hiF, j F,k B ihF, j F,k

NihF, j F,k = , E ihF, j F,k = ,
A hiF, j F,k + A hiF, j F +1,k B ihF, j F,k + B ihF +1, j F,k

A hiF, j F −1,k B ihF −1, j F,k

S ihF, j F,k = , WihF, j F,k = ,
A hiF, j F −2,k + A hiF, j F −1,k B ihF −1, j F,k + B ihF −2, j F,k

16

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 2. NUMERICAL MODEL

Figure 2.3: Interpolation process (green points: coarse grid points, black dots: fine grid points,
blue dots: middle points on the fine grid, red point: central point on the fine grid).

Figure 2.4: Weighting factors for the interpolation (blue points: coarse grid points, black dots:
fine grid points).

17

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 2. NUMERICAL MODEL

[B ihF, j F +1,k A hiF, j F,k /(A hiF, j F,k + A hiF, j F +1,k )]

N E ihF, j F,k =
A hiF +1, j F +1,k + A hiF +1, j F,k + B ihF +1, j F +1,k + B ihF, j F +1,k
[B ihF, j F,k A hiF +1, j F,k /(B ihF, j F,k + B ihF +1, j F,k )]
+ ,
A hiF +1, j F +1,k + A hiF +1, j F,k + B ihF +1, j F +1,k + B ihF, j F +1,k

[B ihF −1, j F +1,k A hiF, j F,k /(A hiF, j F,k + A hiF, j F +1,k )]
NWihF, j F,k =
A hiF −1, j F +1,k + A hiF −1, j F,k + B ihF −2, j F +1,k + B ihF −1, j F +1,k
[B ihF −1, j F,k A hiF −1, j F,k /(B ihF −1, j F,k + B ihF −2, j F,k )]
+ ,
A hiF −1, j F +1,k + A hiF −1, j F,k + B ihF −2, j F +1,k + B ihF −1, j F +1,k

[B ihF, j F −1,k A hiF, j F −1,k /(A hiF, j F −1,k + A hiF, j F −2,k )]

SE ihF, j F,k =
A hiF +1, j F −1,k + A hiF +1, j F −2,k + B ihF, j F −1,k + B ihF +1, j F −1,k
[B ihF, j F,k A hiF +1, j F −1,k /(B ihF, j F,k + B ihF +1, j F,k )]
+ ,
A hiF +1, j F −1,k + A hiF +1, j F −2,k + B ihF, j F −1,k + B ihF +1, j F −1,k

[B ihF −1, j F −1,k A hiF, j F −1,k /(A hiF, j F −1,k + A hiF, j F −2,k )]
SWihF, j F,k =
A hiF −1, j F −1,k + A hiF −1, j F −2,k + B ihF −2, j F −1,k + B ihF −1, j F −1,k
[B ihF −1, j F,k A hiF −1, j F −1,k /(B ihF −1, j F,k + B ihF −2, j F,k )]
+ .
A hiF −1, j F −1,k + A hiF −1, j F −2,k + B ihF −2, j F −1,k + B ihF −1, j F −1,k

2.4.2 Injection
In general, the restriction operator matrix is the transposed matrix of the interpolation operator
[10]:
I hH = (I H
h T
) (2.35)
In order to derive the restriction operator, a basis function e il , j whose value is 1 at the point (i , j )
on the l th grid is employed.
h H h
X X 〈I H e iC , jC , e i F, j F 〉
W H = (I hH W h )(iC , jC ) = [(I H
h T
) W h ](iC , jC ) = H
WihF, j F (2.36)
iF jF 〈e iC ,eH
, jC iC , jC
〉

in which W H and W h are vectors of unknowns on the coarse and fine grid respectively. In terms
H
of the orthogonal basis e il , j , the dot products are 〈e iC , e H 〉 = HxHy and 〈e ihF, j F , e ihF, j F 〉 =
, jC iC , jC
hxhy.
Supposing the vectors of unknowns on coarse grid and on fine grid are marked as W H =
H
[u iC , jC ,k
] and W h = [u ihF −1, j F +1,k , u ihF −1, j F,k , u ihF −1, j F −1,k , u ihF, j F +1,k , u ihF, j F,k ,
u ihF, j F −1,k , u ihF +1, j F +1,k , u ihF +1, j F,k , u ihF +1, j F −1,k ]. According to Equation 2.36, one can obtain:

H hx · h y
u iC , jC ,k = (C h uh + NWihF, j F,k u ihF −1, j F +1,k + SWihF, j F,k u ihF −1, j F −1,k
H x · H y i F, j F,k i F, j F,k
+ SE ihF, j F,k u ihF +1, j F −1,k + N E ihF, j F,k u ihF +1, j F +1,k + NihF, j F,k u ihF, j F +1,k (2.37)

+ S ihF, j F,k u ihF, j F −1,k + WihF, j F,k u ihF −1, j F,k + E ihF, j F,k u ihF +1, j F,k ).

18

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 2. NUMERICAL MODEL

Hence, the stencil of the restriction operator I hH is:

nw ihF, j F,k n ihF, j F,k ne ihF, j F,k

 
 
 
hx · h y 
 wh
I hH c ihF, j F,k e ihF, j F,k

= i F, j F,k
 (2.38)
Hx · Hy 


 
sw ihF, j F,k s ihF, j F,k se ihF, j F,k

where

c ihF, j F,k = 1,

A hiF, j F,k B ihF, j F,k

n ihF, j F,k = , e ihF, j F,k = ,
A hiF, j F,k + A hiF, j F +1,k B ihF, j F,k + B ihF +1, j F,k

A hiF, j F −1,k B ihF −1, j F,k

s ihF, j F,k = , w ihF, j F,k = ,
A hiF, j F −2,k + A hiF, j F −1,k B ihF −1, j F,k + B ihF −2, j F,k

[B ihF, j F +1,k A hiF, j F,k /(A hiF, j F,k + A hiF, j F +1,k )]

ne ihF, j F,k =
A hiF +1, j F +1,k + A hiF +1, j F,k + B ihF +1, j F +1,k + B ihF, j F +1,k
[B ihF, j F,k A hiF +1, j F,k /(B ihF, j F,k + B ihF +1, j F,k )]
+ ,
A hiF +1, j F +1,k + A hiF +1, j F,k + B ihF +1, j F +1,k + B ihF, j F +1,k

[B ihF −1, j F +1,k A hiF, j F,k /(A hiF, j F,k + A hiF, j F +1,k )]
nw ihF, j F,k =
A hiF −1, j F +1,k + A hiF −1, j F,k + B ihF −2, j F +1,k + B ihF −1, j F +1,k
[B ihF −1, j F,k A hiF −1, j F,k /(B ihF −1, j F,k + B ihF −2, j F,k )]
+ ,
A hiF −1, j F +1,k + A hiF −1, j F,k + B ihF −2, j F +1,k + B ihF −1, j F +1,k

[B ihF, j F −1,k A hiF, j F −1,k /(A hiF, j F −1,k + A hiF, j F −2,k )]

se ihF, j F,k =
A hiF +1, j F −1,k + A hiF +1, j F −2,k + B ihF, j F −1,k + B ihF +1, j F −1,k
[B ihF, j F,k A hiF +1, j F −1,k /(B ihF, j F,k + B ihF +1, j F,k )]
+ ,
A hiF +1, j F −1,k + A hiF +1, j F −2,k + B ihF, j F −1,k + B ihF +1, j F −1,k

[B ihF −1, j F −1,k A hiF, j F −1,k /(A hiF, j F −1,k + A hiF, j F −2,k )]
sw ihF, j F,k =
A hiF −1, j F −1,k + A hiF −1, j F −2,k + B ihF −2, j F −1,k + B ihF −1, j F −1,k
[B ihF −1, j F,k A hiF −1, j F −1,k /(B ihF −1, j F,k + B ihF −2, j F,k )]
+ .
A hiF −1, j F −1,k + A hiF −1, j F −2,k + B ihF −2, j F −1,k + B ihF −1, j F −1,k

19

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 2. NUMERICAL MODEL

Subsequently, the solution P H and the residual r H on the coarse grid are restricted as:
H H h
P iC , jC ,k = I h P i F, j F,k
hx · h y h
= [c Ph + n ihF, j F,k P ihF, j F +1,k + s ihF, j F,k P ihF, j F −1,k
H x · H y i F, j F,k i F, j F,k (2.39)
+ w ihF, j F,k P ihF −1, j F,k + e ihF, j F,k P ihF +1, j F,k + nw ihF, j F,k P ihF −1, j F +1,k
+ sw ihF, j F,k P ihF −1, j F −1,k + se ihF, j F,k P ihF +1, j F −1,k + ne ihF, j F,k P ihF +1, j F +1,k ].

Special attention is needed when restricting the residual r h . The factor (hx · h y)/(H x · H y) in
Equation 2.38 is used to restrict functions mathematically, which is not used to restrict flows
physically. For Reynolds equation, the restriction operator I hH must be replaced by J hH when
restricting the right hand side of the discrete Reynolds equation [124] .

nw ihF, j F,k n ihF, j F,k ne ihF, j F,k

 
 
 
Hx · Hy H  h
J hH = c ihF, j F,k e ihF, j F,k

I  w i F, j F,k
= (2.40)

hx · h y h 


 
sw ihF, j F,k s ihF, j F,k se ihF, j F,k
yielding:
H H h
r iC , jC ,k = J h r i F, j F,k . (2.41)

2.5 Coarse grid operator

After defining the transfer operators, the coarse grid operator will be described here. The coarse
grid operator is formed by a restriction equation:
J hH (L h υh ) = J hH r h (2.42)
h h H
In which υ = IH υ and equation 2.42 can be represented as:

J hH L h I H
h H
(υ ) = J hH r h (2.43)
Thus, the coarse grid operator can be defined as:
L H = J hH L h I H
h
(2.44)
Assuming the stencil of L l at point (i , j ) on level l and at k t h time step is:
 l
L nw (i , j , k) L ln (i , j , k) L lne (i , j , k)

 
 
L l =  L lw (i , j , k) L lc (i , j , k) L le (i , j , k)  (2.45)
 
 
 
L lsw (i , j , k) L ls (i , j , k) L lse (i , j , k)
The construction of the coarse grid operator can be seen in Appendix A.
The coarse grid operator defined here involves 9 points, while as equation 2.27 shows an
operator referring to 5 points on the finest grid, reads:
A hiF, j F,k
 
0 0
 
 
h h h h h h h
 
L (i F, j F, k) = B i F −1, j F,k −A i F, j F,k − A i F, j F −1,k − B i F −1, j F,k − B i F, j F,k B i F, j F,k 

 (2.46)
 
 
0 A hiF, j F −1,k 0

20

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 2. NUMERICAL MODEL

2.6 Relaxation

Employing the grid operator L l defined before, the discrete Reynolds equation for point (i , j , k)
of the grid l at the k t h time step is shown as follows:

L lc P il , j ,k + L ln P il , j +1,k + L lw P il −1, j ,k + L ls P il , j −1,k + L le P il +1, j ,k +

L lnw P il −1, j +1,k + L lsw P il −1, j −1,k + L lse P il +1, j −1,k + L lne P il +1, j +1,k = (2.47)
(h l · h l )(ρ̄H )lx + (h l · h l )(ρ̄H )lt

with P il , j ,k = 0 for points on the boundary and P il , j ,k > 0 for the cavitation condition. In Equa-
tion 2.47, h l stands for the mesh size on the l t h grid. In terms of the couette and the transient
term, their discrete forms are:

ρ̄ H l −ρ̄ l
 l
Hl
 i , j ,k i , j ,k il −1, j ,k i −1, j ,k

if i = 1
h
(ρ̄H )lx = (2.48)
 1.5ρ̄ i , j ,k Hi , j ,k −2.0ρ̄ i −1, j ,k Hi −1, j ,k +0.5ρ̄ i −2, j ,k Hi −2, j ,k
l l l l l l

hl
if i Ê 2

and 

 0 if nt = 0

 ρ̄ l

H l −ρ̄ li −1, j ,k Hil−1, j ,k
(ρ̄H )lt =
i , j ,k i , j ,k
ht if nt = 1 (2.49)

ρ̄ l l
−2.0ρ̄ li −1, j ,k Hil−1, j ,k +0.5ρ̄ li −2, j ,k Hil−2, j ,k


 1.5 H
i , j ,k i , j ,k

ht if nt Ê 2
with
X i2,k Y j2,k
Hil, j ,k K il,il 0 , j , j 0 ,k P il 0 , j 0 ,k
XX
= H0 (k) + + + (2.50)
2 2 i0 j0

For the discrete Equation 2.47, Reference [10] shows that an iterative approach of combining
the Gauss-Seidel line relaxation and the Jacobi distribution line relaxation is pretty stable and
efficient. Hence, in this work the same method is applied. Figure 2.5 shows the relaxation
process represented in Reference [10], where ξl is the local coefficient and its definition can
been seen as Equation 2.51.
ρ̄(P il , j ,k )(Hil, j ,k )3
l
ξi , j ,k = (2.51)
η̄(P il , j ,k )λ̄

with ξli ±1/2, j ,k = (ξi , j ,k + ξi ±1, j ,k )/2 and ξli , j ±1/2,k = (ξi , j ,k + ξi , j ±1,k )/2.
Whether the Gauss-Seidel or the Jacobi distribution line relaxation is used, for each grid
point (i , j , k) a new approximation P̄ il , j ,k to P il , j ,k is computed by:

P̄ il , j ,k = P̃ il , j ,k + ωδli , j ,k (2.52)

with
l
ωg s δi , j ,k

Guass-Seidel
ωδli , j ,k = (δli +1, j ,k +δli −1, j ,k +δli , j +1,k +δli , j −1,k )
ω j a [δli , j ,k − ] Jacobi distribution

4

where ωg s and ω j a are the relaxation factors for the Guass-Seidel line relaxation and the Ja-
cobi distribution line relaxation respectively. The changes δli , j ,k for the line relaxation in the X
direction can be obtained from:
A j δlj = r lj (2.53)

21

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 2. NUMERICAL MODEL

in which δlj is a vector of changes δli , j ,k and r lj is the residual vector r il, j ,k , both are n x element
j
vectors. A j is a matrix of coefficients A i ,m , whose derivation can be found in Appendix B.

Figure 2.5: Flow chart of the hybrid relaxation process

2.7 Implementation of the Multi-Grid method

So far, pre-preparations for implementing the Multi-Grid method to solve the transient EHL
equation system have already been prepared. Final step is to organize all steps together. Figure
2.6 shows a simple two-level "V" cycle which is used to illustrate the implementation of the
Multi-Grid method. The corresponding steps are as follows:
h
• Several relaxations on grid h to obtain an approximate solution P̃ by:

h
L h 〈P̃ 〉 = F h , (2.54)

compute residual on grid h below:

h
r h = F h − L h 〈P̃ 〉. (2.55)

• Compute Alcouffe’s coefficients A, B and D, prepare the transfer operators I hH and I H

h
,
H
construct the coarse grid operator L .

22

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 2. NUMERICAL MODEL

Figure 2.6: Implementation of the Multi-Grid method with a two level "V" cycle.

h
• Coarsen pressure P̃ and residual r h to grid H :
H h
P̃ = I hH 〈P̃ 〉
(2.56)
r H = J hH 〈r h 〉

with the right hand side term F H expressed as:

h
F H = L H 〈P̃ 〉 + r H . (2.57)

H
• Relax coarse grid solution P̃ on grid H through:
H
L H 〈P̃ 〉 = F H . (2.58)

if grid H is the coarsest grid, the mutual approach H0 should be updated.

h
• Interpolate and correct P̃ from grid H to h using:
h h h H h
P̄ = P̃ + I H 〈P̃ − I hH P̃ 〉. (2.59)

• Relaxations on grid h.

23

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 2. NUMERICAL MODEL

For simplicity, the current code uses "V" cycles to solve the EHL equation system. In terms of
the transient EHL model, the transient term ∂ρ̄H /∂X should be arranged in calculation cycles.
Here use the equation below to denote the transient Reynolds equation with the first order
backward discretization of transient term:

ρ̄ li , j ,k Hil, j ,k − ρ̄ li −1, j ,k Hil−1, j ,k

l l l l
L 〈P 〉 − (h · h )(ρ̄H )lx l
= (h · h )( l
). (2.60)
ht
In order to construct the transient "V" cycle, Equation 2.60 will be rewritten as:

ρ̄ li , j ,k Hil, j ,k (−ρ̄ li −1, j ,k Hil−1, j ,k )

L l 〈P l 〉 − (h l · h l )(ρ̄H )lx − (h l · h l )( ) = (h l · h l ) . (2.61)
ht ht
Observing Equation 2.61, one finds that the right hand side term of the Reynolds equation is
(−ρ̄ li −1, j ,k Hil−1, j ,k )
(h l ·h l ) ht . As shown in Figure 2.7, the right hand side term of the Reynolds equation
(−ρ̄ li −1, j ,k Hil−1, j ,k )
for the (nt )t h time step is the value of (h l · h l ) ht for the (nt − 1)t h time step.

Figure 2.7: The time-dependent "V" cycles

2.8 Conclusion
In this chapter, the coarse grid construction method proposed by Alcouffe et al. [66] is imple-
mented in the current EHL Multi-Grid code. Instead of the bi-linear interpolation stencil used
in the class code [10], the new interpolation stencil applying Alcouffe’s method is derived. Then
the new restriction stencil is obtained through the Equation 2.35. A nine point coarse grid op-
erator is also constructed by the Galerkin method in Appendix A. Finally, the implementation
of the Multi-Grid method is introduced briefly.

24

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Chapter 3

Friction influence of harmonic surface waviness

Contents
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 Lubricant rheological models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.1 Relative friction coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.2 Numerical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.4 Time-dependent solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.5 Effect of operating conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.6 Effect of surface anisotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.6.1 Longitudinal and transverse wavy cases . . . . . . . . . . . . . . . . . . . . 38
3.6.2 Purely longitudinal wavy case . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

3.1 Introduction
As mentioned in Chapter 1, the friction of interfacial surfaces greatly influences the perfor-
mance of mechanical elements. Studies show that the surface roughness amplitude, wave-
length and anisotropy affect friction. Thus this chapter investigates the effect of surface wavi-
ness on the friction evolution, in which the surface waviness is defined by Equation (3.1). In this
chapter, a relative friction coefficient is proposed to indicate the transition from the full-film to
the mixed lubrication regime. The definition of the relative friction coefficient is illustrated in
Section 3.3. Subsequently, the influence of operating conditions and surface waviness topog-
raphy on friction are studied. Finally, a single friction curve is obtained depending on a new
"Lambda ratio" parameter including operating conditions as well as surface waviness parame-
ters. The waviness is defined as:

−10[max(0, λXx −/aX̄ )2 ] (X − X̄ ) Y

RR(X , Y , T ) = Ai × 10 h cos(2π )cos(2π ) (3.1)
λx /a h λ y /a h

where X̄ = X st +Urat × T with Urat = u 1 /u r , Ai is the initial amplitude of the surface waviness,
λx and λ y are the wavelengths in x and y direction respectively. At the same time, the param-
eter r (r = λx /λ y ) is used to determine the surface waviness anisotropy. For convenience, the
parameter λ is defined as λ = min(λx , λ y ). The exponential term is used to avoid discontinuous
derivatives when the waviness moves into the calculation domain.

3.2 Lubricant rheological models

The lubricant rheology greatly determines the friction in a tribological contact [125]. Once
the rheology of the lubricant is given, the friction force can be obtained by integration of the

25

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

shear stress over the contact area [126]. Figure 3.1 shows a typical shear stress-shear rate re-
lationship for an EHL contact. In the low shear strain rate range, the lubricant is Newtonian
and the shear stress increases linearly with the shear strain rate. This linear model is valid for
predicting hydrodynamic friction [127]. However, in the EHL regime of lubrication, the Newto-
nian assumption is insufficient, especially when large slide-to-roll ratios are employed. Stud-
ies [128–130] show that the lubricant in the inlet zone shows shear-thinning behavior (i.e. the
fluid viscosity decreases with increasing shear stress rate) which influences the film thickness.
Some frequently-used non-Newtonian fluid models are listed as follows:

Figure 3.1: Shear stress-shear rate relationship for the EHL contact.

• Ree-Eyring model (sinh-law model)

In 1977, Johnson and Tevaarwerk [131] proposed a simple constitutive equation for an
isothermal lubricant based on a "nonlinear Maxwell" model, which reads:
1 d τ τ0 τ
γ̇ = + sinh( ) (3.2)
Ge d t η τ0
in which G e is the elastic shear modulus, τ is the stress, τ0 is referred to as the Eyring
stress and γ̇ is the shear strain rate. The strain rate shown as Equation 3.2 consists of
two components where the first elastic term is always neglected [132]. Thus Equation 3.2
reduces to:
τ0 τ
γ̇ = sinh( ). (3.3)
η τ0
Observing Equation 3.3, when the Eyring stress τ0 approaches infinity, the limit of Equa-
tion 3.3 becomes the Newtonian constitutive equation, i.e.:
τ0 τ τ
lim sinh( ) = . (3.4)
τ0 →∞ η τ0 η

• Bair-Winer model

For larger strain rates the lubricant exhibits a limiting shear stress τL which is the thresh-
old of the shear stress. Bair and Winer [133] modified the "Maxwell" model applying this
limiting shear stress, yields:
1 d τ τL τ
γ̇ = + ln(1 − ) (3.5)
G∞ d t η τL

26

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

where G ∞ is the limiting elastic shear modulus.

• Carreau-Yasuda model (power-law model)

Carreau and Yasuda [134] offered a more general model whose equation is:

η = η 0 [1 + (λt γ̇)2 ](n−1)/2 (3.6)

where λt is the time constant for the fluid and n − 1 is the power-law slope.

• Circular fluid model

In order to overcome the difficulty lies in incorporating the Bair-Winer model into the
Reynolds equation, Lee and Hamrock [135] suggested an appropriate lubricant rheolog-
ical circular model, it reads:
τ τ
γ̇ = [1 − ( )2 ]−1/2 (3.7)
η τL

• Actual Ree-Eyring model

Bair [136] proposed an actual Ree-Eyring model for shear-shinning lubricants [137]:

N
τ= x i τi sinh−1 (λi γ̇)
X
(3.8)
i =1
PN
in which λi = η/τi is a characteristic time of the fluid and x i is a weighting factor ( i =1 =
1).

3.3 Methodology

3.3.1 Relative friction coefficient

Even though the previous section described several non-Newtonian rheological model, in this
section, we will use a Newtonian model to describe friction variations resulting from pressure
variations. Two arguments can be used.

• When the shear stress variations are small, a linearised model can be used. See the fol-
lowing example sinh.

• A linear model allows an FFT based sum over all wavelengths which a non-linear model
does not.

Example:
τ0 τ
γ̇ = sinh
η τ0
τ0 τ/τ0
= (e − e −τ/τ0 )
2η
τ0 τ τ2 τ τ2
= (1 + + 2 + ... − 1 + − 2 + ...)
2η τ0 2τ0 τ0 2τ0
τ
≈
η

27

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

For Newtonian lubricant, the above equation can be rewritten as:

du
τ=η (3.9)
dz
where u is the velocity in rolling direction, which is defined as [138, 139]:

1 ∂p 2 z
u= (z − zh) + (u 2 − u 1 ) + u 1 . (3.10)
2η ∂x h

Substituting Equation 3.10 into Equation 3.9, yields:

1 ∂p u2 − u1
τ= (2z − h) + η (3.11)
2 ∂x h
where two parts contribute to the shear stress, the first term is a parabolic part due to the
Poiseuille flow and the second term is a linear part due to the Couette flow. Thus the viscous
shear force on the lower surface (z = 0) is obtained by integrating Equation 3.11 [140]:

∂p h u2 − u1
Ï Ï
f =− ( · )d xd y + (η · )d xd y. (3.12)
∂x 2 h

In the full-film EHL regime, the friction force is dominated by the viscous shear force, the di-
mensionless friction force is:
∂P (X , Y , T ) H (X , Y , T ) SRR
Ï Ï
F friction = − ( · )d X d Y + (η̄ · )d X d Y (3.13)
∂X 2 H (X , Y , T )

in which SRR = (u 2 − u 1 )/ū is the slide-to-roll ratio.

Experimental results [141] also found that the Poiseuille force is lower than the Couette
force. Hence the dimensionless friction force is simplified as:

SRR
Ï
F (X , Y , T ) = η̄[P (X , Y , T )] · d X dY . (3.14)
H (X , Y , T )

Different lubricant oils give different curves of friction coefficient. For the sake of simplicity, a
relative friction coefficient is proposed to detect friction variations, which is defined as:

µr Fr F s Fr
(T ) = (T )/ = (T ) (3.15)
µs w w Fs

where subscripts r and s are used to distinguish the rough and smooth case.

3.3.2 Numerical solution

The numerical simulation uses the domain −2.5 É X É 1.5 and −2.0 É Y É 2.0 with 513 × 513
equal-spaced points. The time step is selected equal to the spatial mesh size on the finest grid,
i.e. with ∆T = hx = h y = 0.0078125. Meanwhile, the calculation starts with X st = −2.5 and the
surface topography moves into the high pressure zone with the velocity of the rough surface u 1 .
The monitoring time should be long enough so that ‘steady oscillations’ of the results occur.
The present work considers small-amplitude roughness and a small slip parameter is selected
i.e. Ur at = 1.01(SRR = 0.02). This small slip assumption and and small amplitude allow us to
use the Amplitude Reduction Theory [53] for pure rolling, as shown in [61].

28

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

The numerical code used in the present work is validated by comparing the current results
with results from Reference [53]. Figure 3.2 shows the relative deformed amplitude (Ad /Ai ) as
a function of f (r )∇2 . This figure shows that results from the current model with SRR = 0 (blue
squares) fall onto the master curve on Reference [53].

Figure 3.2: Comparison of the relative deformed amplitude (Ad /Ai ) as a function of f (r )∇2 for
the current model (blue squares) with those on Reference [53] (solid line)

The choice of the mesh size influences computing time as well as precision. A large mesh
size leads to large discretization errors and a small mesh size causes long computation times.
This work applies an intermediate mesh size (513 × 513 points). Table 3.1 suggests that it yields
a friction precision better than 1%.

Table 3.1: Relative friction coefficient versus the number of mesh points for: M =1000, L=10,
λx /a h = 0.5, r=0.4 and H c/Ai = 2.

Mesh points Relative friction coefficient

129 × 129 1.489
257 × 257 1.457
513 × 513 1.462
1025 × 1025 1.465

3.4 Time-dependent solution

Defined by Equation 3.1, Figure 3.3 shows the top view of the surface waviness for three cases:
the isotropic surface waviness of r = 1 and λ/a h = 0.5 (a), the longitudinal surface waviness of
r = 2 and λ/a h = 0.5 (b) and the transverse surface waviness of r = 0.5 and λ/a h = 0.5 (c).

29

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

Figure 3.3: Top view of the surface waviness with λ/a h = 0.5 and Ai = 0.5H c: (a) the isotropic
surface waviness r = 1, (b) the longitudinal surface waviness r = 2, (c) the transverse surface
waviness r = 0.5.

For the transient case, in the high pressure zone the Poiseuille term vanishes, the Reynolds
equation reduces to a transport equation:

∂(ρ̄H ) ∂(ρ̄H )
− − ≈ 0. (3.16)
∂X ∂T
In this case, the final film thickness and pressure profiles depend on the profiles of steady state
case as well as on the inlet disturbances [46]. Figures 3.4-3.6 show dimensionless pressure dis-
tribution P (X , 0) and dimensionless film thickness distribution H (X , 0) for three surface wavi-
ness cases, in which H c is the central film thickness for the smooth case. For the operating
condition of M = 1000 and L = 10, the dimensionless central film thickness for the smooth
case is: H c = 0.0243. In this work, we only focus on small surface waviness amplitude values
i.e. Ai /H c ∈ [1, 10]. Figures 3.4-3.6 show that the propagation speed of the pressure increase
is the same as that of the film thickness variation. This is due to the small slide-to-roll ratio,
the velocity of the rough surface is nearly the same as the entrainment velocity ū. In addition,
the wavelength of the pressure increase and the film thickness variation are the same as the
wavelength of corresponding initial surface waviness.
Figure 3.7 shows the dimensionless central film thickness H c rough (0, 0, T ) as a function of
time T for the above three rough cases. The initial amplitude of the surface waviness is Ai =
0.01215, while from those figures, one can observe that the amplitude of the central film thick-
ness is smaller. This means that the surface roughness is deformed under current operating
conditions. For those cases in Figure 3.7, the deformed amplitudes are: Ad = 0.0057 for r = 1,
Ad = 0.00375 for r = 2 and Ad = 0.0057 for r = 0.5. Reference [53] studied this deformation
over a large range of operating conditions, and a single equation is obtained as Equation 3.17.

Ad 1
= (3.17)
Ai 1 + 0.15 f (r )∇2 + 0.015( f (r )∇2 )2
where (
e 1−1/r , if r > 1
f (r ) =
1, otherwise
p
r = λx /λ y and ∇2 = (λ/a h ) M /L with λ = min(λx , λ y ).

30

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

Figure 3.4: Central line pressure P (X , 0) (black lines) and central line film thickness H (X , 0)
(blue lines) of isotropic surface waviness (r = 1) for M = 1000, L = 10, λ/a h = 0.5 and Ai = 0.5H c
during a time period. The central pressure line (red line) for the smooth case is plotted as a
reference.

31

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

Figure 3.5: Central line pressure P (X , 0) (black lines) and central line film thickness H (X , 0)
(blue lines) of longitudinal surface waviness (r = 2) for M = 1000, L = 10, λ/a h = 0.5 and Ai =
0.5H c during a time period. The central pressure line (red line) for the smooth case is plotted
as a reference.

32

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

Figure 3.6: Central line pressure P (X , 0) (black lines) and central line film thickness H (X , 0)
(blue lines) of transverse surface waviness (r = 0.5) for M = 1000, L = 10, λ/a h = 0.5 and Ai =
0.5H c during a time period. The central pressure line (red line) for the smooth case is plotted
as a reference.

33

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

Figure 3.7: The dimensionless central film thickness H c rough as a function of the dimensionless
time T for: M = 1000, L = 10, λ/a h = 0.5 and Ai = 0.5H c: (a) the isotropic surface wavy case,
(b) the longitudinal surface wavy case r = 2, (c) the transverse surface wavy case r = 0.5.

Equation 3.15 shows that the relative friction coefficient is time-dependent. Figure 3.8 de-
picts the relative friction coefficient as a function of dimensionless time (T ) for the isotropic,
longitudinal and transverse cases with the same operating condition: M=1000, L=10, Ai =
0.5H c and dimensionless waviness: λ = 0.5. Because the surface topography RR(T ) is periodi-
cal, the relative friction coefficient (µr /µs )(T ) is also periodical. For each operating condition,
the relative friction coefficient is defined as its average value (shown as blue dotted lines in Fig-
Pi =N µ
ure 3.8) i.e. i =0 T ( µrs )i /NT in which NT is the number of the relative friction coefficient values
during one time period. Although the operating conditions for the three surfaces are the same,
their averaged friction values are very different i.e. µr /µs = 1.79 for r = 1, µr /µs = 1.57 for r = 2
and µr /µs = 1.48 for r = 0.5.

Figure 3.8: The relative friction coefficient µr /µs as a function of the dimensionless time T for:
M = 1000, L = 10, λ/a h = 0.5 and Ai = 0.5H c: (a) isotropic surface wavy case, (b) longitudinal
surface wavy case r = 2, (c) transverse surface wavy case r = 0.5. (Blue dotted line: average
value of the relative friction coefficient.)

3.5 Effect of operating conditions

The relative friction coefficient can be plotted as a function of the classical parameter "lambda
ratio" i.e. H c/Ai (Ai is varying and H c is fixed) for a specific operating condition: Figure 3.9.
This figure shows that as H c/Ai increases, the relative coefficient decreases monotonically.

34

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

Early work on conformal condition shows that with the "Λ ratio" increasing, the friction coef-
ficient decreases to a minimum value (the transition position from the full-film to the mixed
lubrication regime) first then increases.

Figure 3.9: Relative friction coefficient as a function of H c/Ai for a specific operating condition

The following work will detect the influence of the operating conditions M , L and surface
wavelength on the relative friction coefficient for the isotropic case.

Figure 3.10: Effect of the load parameter M on the relative friction coefficient for L = 10 and
λ/a h = 0.5: (a) relative friction coefficient as a function of H c/Ai , (b) relative friction coefficient
as a function of M 0.33 · (H c/Ai )

Figure 3.10 displays the influence of the load parameter M on the relative friction coeffi-
cient for L = 10 and λ/a h = 0.5. From Figure 3.10 (a), one can observe that for each M value,
there is a single curve. As M increases, the relative friction coefficient decreases, meanwhile
each single curve has a small left-shift. It seems that there is a scaling factor can be used to
combine all curves together. Figure 3.10 (b) shows that M 0.33 is an appropriate factor to scale
all curves.

35

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

Figure 3.11 shows the influence of the material parameter L on the relative friction coeffi-
cient for M = 1000 and λ/a h = 0.5. From Figure 3.11 (a), one can observe that the parameter L
affects the relative friction coefficient. As L increases, the relative friction coefficient increases,
each single curve has a small right-shift. Once again, there exists a scaling parameter L −1.1
which can be applied to shift all curves together (shown in Figure 3.11 (b)).

Figure 3.11: Effect of material parameter L on the relative friction coefficient for M = 2000 and
λ/a h = 0.5: (a) relative friction coefficient as a function of H c/Ai , (b) relative friction coefficient
as a function of L −1.1 · (H c/Ai )

Figure 3.12: Effect of wavelength λ/a h on the relative friction coefficient for M = 1000 and
L = 10: (a) relative friction coefficient as a function of H c/Ai , (b) relative friction coefficient as
a function of (λ/a h )0.67 · (H c/Ai )

The relative friction coefficient is affected not only by the operating conditions M and L,
but also by the wavelength of the surface waviness. Figure 3.12 describes the effect of the wave-
length λ/a h on the relative friction coefficient for M = 1000 and L = 10. As λ/a h increases, the
relative friction coefficient decreases, each single curve has a small left-shift (shown in Figure

36

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

3.12 (a)). The parameter λ/a h0.67 is a suitable factor to unify those curves (shown in Figure 3.12
(b)).
Figure 3.13 presents the relative friction coefficient as a function of the classical "lambda
ratio" H c/Ai for many different operating conditions. This figure shows that for each operating
condition, a very smooth curve is obtained, however, one does not obtain a single curve as
for the low pressure case using the parameter H c/Ai . According to the Amplitude Reduction
Theory [53, 58], for very high pressures, the surface roughness will deform. Instead of using
this simple parameter Ai or a measured surface roughness parameter σ, it is better to use the
deformed parameter Ad .

Figure 3.13: Relative friction coefficient as a function of the classical parameter "lambda ratio"
i.e. H c/Ai for a large range of operating conditions.

Employing the scaling factors of Figures 3.10-3.12, it is possible to combine all results ob-
tained for different values of λ/a h , M , L as well as H c/Ai into a single curve using a dimen-
sionless parameter θ2 . Figure (3.14) shows the relative friction coefficient as a function of the
new parameter θ2 for 500 ≤ M ≤ 2000, 5 ≤ L ≤ 15, 0.25 ≤ λ/a h ≤ 1.0 and 1.0 ≤ H c/Ai ≤ 10. After
curve-fitting, the single curve can be described by the following equation:
µr
= 1 + 0.56 θ2−2 + 0.23 θ2−4 (3.18)
µs

where θ2 = L −1.1 M 0.33 (λ/a h )0.67 (H c/Ai ). The physical justification of this scaling parameter
can be understood from a simplified analysis given in Appendix C.

37

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

Figure 3.14: Relative friction coefficient as a function of the new parameter θ2 , simulation re-
sults: black circles; fitted curve: the black dashed line.

3.6 Effect of surface anisotropy

The previous section demonstrated the influence of the operating conditions and the wave-
length on the relative friction coefficient for the isotropic case. Figure 3.8 shows that the surface
anisotropy affects the relative friction coefficient under the same operating conditions. This
section describes the influence of the surface anisotropy on the relative friction coefficient in
more detail.

3.6.1 Longitudinal and transverse wavy cases

Figure 3.15 shows the relative friction coefficient (µr /µs ) as a function of H c/Ai for different r
(1 ≤ r ≤ 32) values. It can be observed that the results of the longitudinal case show the same
trend as those of the isotropic wavy case (r = 1). That is, a decreasing trend of the relative fric-
tion coefficient when H c/Ai increases. For each r value a single curve exists but it shifts to the
left compared to the curve with r = 1. For the case considered, the relative friction coefficient
monotonically decreases as r increases (shown in the right small figure in Figure 3.15).

Figure 3.15: Relative friction coefficient (µr /µs ) as a function of H c/Ai for different r (1 ≤ r ≤
32) values for: M = 1000, L = 10 and λ y /a h = 0.5 (left), zoom from 2.3 − 2.7 (right).

38

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

From Figure 3.15, one can conclude that the results for r = 16 and for r = 32 are almost
identical and the left shift is limited. It is assumed that a function can be used to shift all results
onto a single curve. The single curve is obtained as shown in Figure 3.16 when the horizontal
coordinate H c/Ai in Figure 3.15 is multiplied by the following function:

1
ff(r ) = 1.23 − 0.23( )1.4 (3.19)
r
in order to distinguish the surface roughness anisotropy function in Equation 3.17, in this work,
the anisotropy function is marked as ff(r ).

Figure 3.16: Relative friction coefficient (µr /µs ) as a function of ff(r ) · (H c/Ai ) for different r
(1 ≤ r ≤ 32) values for: M = 1000, L = 10 and λ y /a h = 0.5.

The transverse case is studied for M = 1000, L = 10 and λx /a h = 0.5. Figure 3.17 shows the
relative friction coefficient (µr /µs ) as a function of H c/Ai for different r (0 ≤ r ≤ 1) values. Re-
sults for the transverse case are more complicated. It can be observed from Figure 3.17 that as
r increases (r varies from 0 to 0.33), the relative friction coefficient decreases gradually (shown
as blue lines) while for r varying from 0.4 to 1, the relative friction coefficient increases (shown
as black lines).
Once again, a function can be found to scale all results and form a single curve as shown in
Figure 3.18.
2.7r 2 + 0.1
ff (r ) = 3 (3.20)
r + 1.7r 2 + 0.1
Combining all anisotropy factors, the anisotropy function is expressed as:
(
(2.7r 2 + 0.1)/(r 3 + 1.7r 2 + 0.1), if r < 1
ff (r ) = (3.21)
1.23 − 0.23(1/r )1.4 , otherwise

where r = λx /λ y . Figure 3.19 shows the comparison of the simulation results of the anisotropy
factors and the final fitted function (shown as Equation 3.21).

39

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

Figure 3.17: Relative friction coefficient (µr /µs ) as a function of (H c/Ai ) for different r (0 ≤ r ≤
1) values for: M = 1000, L = 10 and λx /a h = 0.5(upper), zoom from 2.4 − 2.6 (lower).

Figure 3.18: Relative friction coefficient (µr /µs ) as a function of ff(r ) × (H c/Ai ) for different r
(0 ≤ r ≤ 1) values for: M = 1000, L = 10 and λx /a h = 0.5.

40

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

Figure 3.19: ff (r ) as a function of r . Numerical results: red squares. Fitted curve: solid lines.

3.6.2 Purely longitudinal wavy case

Study [51] showed that the transient purely longitudinal case remains a stationary problem.
This phenomenon can also be seen in the current study. Figure 3.20 presents the compari-
son between the transient results and the stationary results for M = 1000, L = 10, H c/Ai = 2,
λ y /a h = 1.0 and λx = ∞. In order to compare these results, the stationary result is also plotted
as a function of the dimensionless time T in Figure 3.20. Figure 3.20 shows that the value of
the relative friction coefficient is the same as the value of the stationary case, i.e. µr /µs = 1.673.
This is due to the waviness term (Equation 3.1) reduces to RR(X , Y ) = Ai × cos(2π λyY/ah ), which
is a time-independent term.

Figure 3.20: Comparison between the transient relative friction coefficient and that of the sta-
tionary case. Transient results: black line. Stationary results: magenta dash-dotted line.

41

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

Figure 3.21: Relative friction coefficients for the purely longitudinal wavy case: (a) relative fric-
tion coefficient as a function of original "lambda ratio" H c/Ai parameter, (b) relative friction
coefficient as a function of the new parameter θ2∗

The influence of the waviness anisotropy on friction can be found in Table 3.2. It is reported
that with r increasing from 1 to 32 the relative friction coefficient (µr /µs ) decreases, while the
case r = ∞ predicts the highest relative friction coefficient value.

Table 3.2: Relative friction coefficient versus different surface anisotropy parameters for: M =
1000, L = 10, λ y /a h = 0.5 and H c/Ai = 2.

r 1 2 4 8 16 32 ∞
µr /µs 1.789 1.566 1.508 1.491 1.489 1.487 2.041

Figure 3.21 gives the results of the relative friction coefficient for the purely longitudinal
wavy case for a relative large range of operating conditions: M ∈ [500, 2000], L ∈ [5, 15] and
λ y /a h ∈ [0.25, 0.8]. As usual, in Figure 3.21 (a), a very smooth curve as a function of H c/Ai for
each operating condition. While, the old scaling parameter θ2 can not scale all cases together.
For the purely longitudinal cases considered, there is indeed a new parameter θ2∗ expressed
as Equation (3.22), and shown in Figure 3.21 (b) to unify all cases onto a single curve. The
operating conditions are: 500 ≤ M ≤ 2000, 5 ≤ L ≤ 15 and 0.25 ≤ λ/a h ≤ 0.8.

θ2∗ = M 0.4 L −1.2 [−1.1(λ y /a h )2 + 2.1(λ y /a h ) + 0.2](H c/Ai ). (3.22)

42

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

3.7 Conclusion
In this chapter, the relative friction coefficient is predicted by means of numerical simula-
tion for a large range of transient EHL operating condition and surface waviness anisotropy
(r = λx /λ y ). A relation is derived for the relative friction coefficient for isotropic and anisotropic
harmonic surfaces. For the isotropic wavy case, all relative friction coefficient values can be
unified onto a single curve applying a dimensionless parameter θ2 . For anisotropic wavy cases,
a roughness anisotropy function ff (r ) is used to scale all results together, when the same op-
erating conditions are considered. Finally, results of all cases can be combined into a single
equation:
µr
= 1 + 0.56[ff (r )θ2 ]−2 + 0.23[ff (r )θ2 ]−4 (3.23)
µs
where (
(2.7r 2 + 0.1)/(r 3 + 1.7r 2 + 0.1), if r < 1
ff (r ) =
1.23 − 0.23(1/r )1.4 , otherwise

and θ2 = M 0.33 L −1.1 (λ/a h )0.67 (H c/Ai ) with λ = min(λx , λ y ).

However, from a lubrication point of view, the purely longitudinal wavy case (r = ∞) has
a very different frictional behavior. A separate scaling parameter θ2∗ is obtained for the purely
longitudinal wavy case.
Early work on conformal contacts showed that the onset of mixed lubrication regime was
roughly around a "lambda ratio" equal to 3, in which the "lambda ratio" is defined as the ratio of
oil film thickness to the combined surface roughness. The current work reveals that besides the
classical "lambda ratio" parameter, the operating conditions as well as the surface topography
play an important role on the onset of the mixed lubrication regime.

43

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Chapter 4

Friction of complex rough surfaces

Contents
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.2 Power spectral density friction method . . . . . . . . . . . . . . . . . . . . . . . 44
4.2.1 PSD friction model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.2.2 Model validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.3 The artificial surface roughness . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.3.1 Surface roughness power spectrum . . . . . . . . . . . . . . . . . . . . . . 47
4.3.2 Friction increase prediction of a rough surface . . . . . . . . . . . . . . . . 48
4.3.3 Comparison between the EHL simulation and the PSD prediction . . . . 50
4.4 Measured surface roughness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.4.1 Friction prediction under a specific operating condition . . . . . . . . . . 54
4.4.2 Operating condition effects . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.4.3 Friction curves for measured surface roughness . . . . . . . . . . . . . . . 59
4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

4.1 Introduction
The previous chapter has already predicted the friction increase for simple harmonic surface
waviness, while the range of wavelengths is limited. In reality, the real surface roughness is
more complicated. It contains very different wavelengths and amplitudes. On the other hand,
the full numerical simulation is time-consuming for real rough surfaces. In this chapter a rapid
prediction method based on the roughness power spectral density (PSD) is provided to pre-
dict the friction increase due to the roughness. Section 4.2 will describe this method in the
first place. Section 4.3 tests this method for an artificial fractal surface roughness. Finally, sec-
tion 4.4 predicts the transition of the full-film to the mixed lubrication regime by applying this
method for real measured surface roughness.

4.2 Power spectral density friction method

4.2.1 PSD friction model

The power spectral density (PSD) is a mathematical tool that decomposes a rough surface into
harmonic components of different frequencies [109], which enables the pressure increase to
be calculated analytically for each frequency component. Subsequently, the friction variations
for the whole rough surface can be obtained. At last, the relative friction coefficient is obtained.
The calculation process is as follows:

44

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

A rough surface topography r r x,y can be expressed in the frequency domain by means of
the Fourier transform:
4
) (r r x,y )e −i (q x x+q y y)
X
r r q x ,q y = ( (4.1)
N x N y x,y
where r r x,y is the discrete form of the surface roughness r r (x, y), q x and q y are the wavenum-
bers in the x and y direction respectively. In general, Equation 4.1 is computed by the fast
Fourier transform (FFT) algorithm.
Combine Equation 4.1 and Equation 3.17, the deformed surface roughness r r qdx ,q y in the
frequency domain is:
Ad
r r qdx ,q y = ( )q ,q · |r r q x ,q y |. (4.2)
Ai x y
According to the relation between the pressure and the elastic deformation of the waviness
given in Appendix D, the pressure increase in the frequency domain follows the expression
below
πE 0
δp q x ,q y = (r r q x ,q y − r r qdx ,q y ) (4.3)
2k(r )λ
where λ is defined as λ = min(2π/q x , 2π/q y ). With the inverse discrete Fourier transform, the
pressure increase in the space domain is obtained:

4
δp x,y = ( δp q x ,q y · e i (q x x+q y y) .
X
) (4.4)
N x N y q x ,q y

According to the friction force Equation 3.14, the ratio of the shear stress τr /τs can be derived
as:
τr (x, y) η r (x, y) h s (x, y) η r (x, y) h s (x, y)
= · = · (4.5)
τs (x, y) η s (x, y) h r (x, y) η s (x, y) h s (x, y) − a d (x, y)
It is easy to obtain the shear stress distribution τs (x, y) for the smooth surface case, where the
pressure distribution for the smooth surface case can be replaced by a semi-elliptical pressure
distribution: ( p
p h 1 − (x/a h )2 − (y/a h )2 , if x 2 + y 2 ≤ a h2
p s (x, y) = (4.6)
0, otherwise.
Afterwards, the pressure distribution for roughness cases is computed by p s + δp. The shear
stress distribution τr (x, y) for a rough surface case is obtained as:

η r (x, y) h s (x, y)
τr (x, y) = · · τs (x, y). (4.7)
η s (x, y) h s (x, y) − a d (x, y)

Finally, the shear forces for both the smooth case and the rough case are computed by integrat-
ing the shear stress τs (x, y) and τr (x, y), respectively. The relative friction coefficient is then
calculated according to Equation 3.15. A detailed description of the prediction process of the
relative friction coefficient is shown in Figure 4.1.

45

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

Figure 4.1: Flow chart for the relative friction coefficient prediction

4.2.2 Model validation

To validate the model described in subsection 4.2.1, the relative friction coefficient evaluated
from a full numerical simulation is compared with that predicted by PSD under the same op-
erating conditions. The same calculation domain, mesh size and slide-to-roll ratio as in sub-
section 3.3.2 are used here. An artificial fractal rough surface is chosen to validate the model
mentioned in the previous subsection and the operating condition parameters are listed in Ta-
ble 4.1.

Table 4.1: Operating condition parameters.

Parameter Value Units

w 600 N
ur 0.84 m/s
Rx 0.018 m
E0 2.26 × 1011 Pa
α 2.2 × 10−8 Pa−1
η0 4 × 10−2 Pa · s
hc 0.233 µm
σ 5 × 10−8 m
Lx = L y 8.29 × 10−4 m
qr 0 m−1
Hurst exponent 0.8

46

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

Table 4.2 presents the friction ratio as a function of the number of mesh points for the two
methods. The results predicted by the two schemes are basically identical. The ratio of the
friction coefficients predicted by the PSD method changes slightly (< 0.08%) with decreasing
mesh size. However, in the full numerical simulation, a large mesh size leads to a relative large
error. This is because some high frequency components of the rough surface can not be cor-
rectly represented on such a large mesh size. In this chapter, the precision of the numerical
results simulated by a 513 × 513 points is considered acceptable.

Table 4.2: Relative friction coefficients as a function of the mesh points for two prediction
schemes.
Mesh points N x × N y hx = h y µr /µs (PSD) µr /µs (EHL)
257 × 257 1/64 1.534 1.487
513 × 513 1/128 1.534 1.518
1025 × 1025 1/256 1.534 1.510

4.3 The artificial surface roughness

In this section, an artificial fractal surface roughness is selected to test the rapid friction predic-
tion method.

4.3.1 Surface roughness power spectrum

Many parameters like the root-mean square (R q ), standard deviation (σ), skewness (Sk), kurto-
sis (K ) et al. are employed to describe a measured surface roughness [142]. The power spectral
density (PSD) is perhaps the most used method for the surface description. Reference [143]
gives the calculation of the PSD for a surface roughness:

C q2D
x ,q y
= (L x L y )−1 |r r q x ,q y |2 (4.8)

where L x and L y are the length of the measured surface roughness profiles in the x and y di-
rections. A 2D surface roughness is shown in Figure 4.2 (a) and it can be represented as a 2D
PSD C 2D in frequency space. When the surface is isotropic, its power spectral density is radially
symmetric (reported in Figure 4.2 (b)).

Figure 4.2: Surface roughness (a) and its 2D power spectral density (b)

47

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

Many real rough surfaces show self-affine characteristic. A surface is self-similar when it
has the same statistical properties as a magnified version of itself [144].The power spectral den-
sity of a self-affine surface often follows the following equation [109, 145] and shown in Figure
4.3: 
C 0 , if q L ≤ q < q r


C i so (q) = C 0 q −2−2H , if q r ≤ q < q s (4.9)


0, otherwise
q
in which q = q x2 + q y2 , q L is the long-wavelength cut-off wave vector, q r is the long-wavelength
roll-off wave vector and q s is the short-wavelength cut-off wave vector. H is the Hurst exponent,
which is related to the fractal dimension D f of the surface [144] : H = 3 − D f .
Once the power spectral density of a rough surface is obtained, other parameters like the
00
0
RMS roughness h rms , RMS slope h rms and RMS curvature h rms can be derived [146].

Figure 4.3: Power spectral density C i so of the self-affine surface(Figure 4.2(a)) with H = 0.8.

4.3.2 Friction increase prediction of a rough surface

The friction increase is determined by the non-linear viscosity variations from pressure vari-
ations caused by roughness deformation. In this subsection, these variables are plotted by
employing the artificial surface roughness with q r = 50000 m−1 shown in Figure 4.4.

Figure 4.4: The selected artificial surface roughness (a), amplitude distribution of this surface
roughness (b) and its power spectral density (c).

48

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

The first step is to compute the deformed surface roughness based on the Amplitude Re-
duction Theory [53] with the operating condition in Table 4.1. Figure 4.5 (a) shows the ratio
of the deformed amplitude and the initial amplitude described by Equation 3.17. Combining
the initial amplitude distribution in Figure 4.4, the deformed surface roughness distribution is
depicted in Figure 4.5 (b).

Figure 4.5: The ratio of the deformed amplitude and the initial amplitude fitted as Equation
(3.17) (a) and the deformed surface roughness in frequency domain (b).

Figure 4.6 compares the initial surface roughness with the deformed surface roughness. It
can be observed that the initial surface roughness is deformed.

Figure 4.6: Comparison between the initial surface roughness (a) and the deformed surface
roughness (b).

Subsequently, the pressure increase can be computed by Equation 4.3 and Equation 4.4.
Figure 4.7 gives the pressure increase distribution in frequency and space domains, respec-
tively. Then the pressure distribution as well as film thickness distribution for the rough case
can be obtained as p s + δp and h s − r r d respectively. Shear stress for both cases are predicted
by η(u 2 − u 1 )/h, which are shown in Figure 4.8. Those values are unrealistically high, and more
complex rheological models are required for an absolute shear stress prediction, but here we
are only interested in the relative shear. The shear stress for rough case is higher than that of
smooth case, this is because pressure variations caused by roughness deformations make a big
contribution to the viscosity. Finally, the friction force is calculated by Equation 3.14 and the
relative friction coefficient is given by Equation 3.15. For this rough surface, friction forces for

49

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

the smooth case and the rough case are F s = 1.0606×106 N and F r = 1.3624×106 N, respectively.
The relative friction coefficient is µr /µs = 1.2846.

Figure 4.7: Pressure increase distribution in frequency (a) and space (b) domains, respectively.

Figure 4.8: Shear stress distributions for the smooth case (a) and for the rough case (b).

4.3.3 Comparison between the EHL simulation and the PSD prediction
The chosen artificial surface topography is generated by means of fractals without a roll-off
region. The friction prediction only occurs in the high pressure zone (X 2 + Y 2 ≤ 1.0), hence
in the full numerical EHL simulation this generated surface roughness is located in the high
pressure zone (shown in Figure 4.9 (b)). In order to make results keep periodical, the "patch"
also needs to be periodical shown in Figure 4.9 (c).

Figure 4.9: The generated surface roughness patch (a), the roughness patch in the high pressure
zone (b) and the periodical roughness pattern for full the numerical simulation (c).

50

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

Figure 4.10: Top view of the deformed surface roughness for a specific time step for a full nu-
merical simulation (a) and for a PSD prediction (b). Central line r r d (x, 0) of the deformed sur-
face roughness for the full numerical simulation (c) and for the PSD prediction (d). Central
line p(x, 0) of the pressure distribution for the full numerical simulation (e) and for the PSD
prediction (f).

As shown in the previous subsection, one can find that the deformed amplitude is a fun-
damental intermediate variable. It is necessary to compare this parameter in the first place.
The resulting deformed micro-geometry, of which the original surface topography is shown in
Figure 4.9 (a) for a full numerical EHL simulation and a PSD prediction, are presented in Figure
4.10. In terms of numerical simulation results, the deformed micro-geometry r r d is obtained
by h s −h r and removing data outside the high-pressure zone (X 2 +Y 2 ≤ 1). Once again, for this
specific surface, the operating conditions are the ones given in Table 4.1 where M = 1000 and

51

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

L = 10. Both the deformed surface topographies r r d are shown in the same high pressure re-
gion. The maximum Hertzian pressure reaches 1.66 GPa and the maximum surface roughness
height deformed significantly from 1 × 10−7 m to 3.5 × 10−8 m. In addition, it is shown that the
height distribution of the deformed surface roughness from the EHL simulation and the PSD
prediction are very similar. Magnitudes of the central line r r d (x, 0) and p(x, 0) for both meth-
ods are also similar. However, the full numerical simulation shows a more "smooth" results in
terms of the deformed surface roughness, this is because small wavelength components can
not be represented correctly in full numerical simualtion. On the other hand, the Amplitude
Reduction Theory is valid for small roughness amplitude values, this is the reason for some
components shown in subfigure(d) whose maginitude is a little higher. In addition, the pres-
sure distribution predicted by the PSD method is the combination of a semi-elliptical pressure
distribution and the pressure deflection i.e. p s + δp x,y .

Figure 4.11: The relative friction as a function of dimensionless time employing the surface
roughness pattern in Figure 4.9 (c) for the full numerical simulation method.

Figure 4.11 shows the periodical relative friction coefficient variation when the surface rough-
ness pattern in Figure 4.9 (c) is used. As mentioned before, the relative friction coefficient value
should be the averaged value of the relative friction coefficient in one time period. The PSD
method only gives the value of the relative friction coefficient for one time step, hence one
should predict this value for each time step within one time period.
Twenty artificial random rough surfaces are generated (shown in Figure 4.12) with the same
input parameters i.e. the standard deviation σ = 5 × 10−8 m, lengths of final topography L x =
L y = 8.29×10−4 m, roll–off wave number q r = 0 m−1 and Hurst exponent=0.8. These generated
roughness are used to compare results from the full numerical simulation as well as the PSD
prediction.
The averaged relative friction coefficients for these artificial random rough surfaces are
given in Table 4.3, showing that the two different prediction methods give close results. The
average deviation is around 8%. It seems that the averaged relative friction values simulated
by the full numerical simulation are higher than those predicted by PSD. This is beacuse in
PSD prediction pressure spike effects are not taken into consideration. Here we need to notice
that when the number of mesh points is selected as 513 × 513, the calculation time of the full
numerical simulation is almost 3 days, this is because more than 2000 time steps are needed.
Meanwhile the calculation time of the PSD prediction is only 15 minutes.

52

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

Figure 4.12: Top view of the twenty generated artificial random rough surfaces from N◦ 1 to
N◦ 20 with a same standard divation value σ = 0.05µm and a same set of operating conditions
listed in Table 4.1.

Table 4.3: The relative friction coefficient obtained by EHL simulation and PSD prediction for
20 artificial random rough isotropic surfaces.

Surface number µr /µs (EHL) µr /µs (PSD) Deviation(%)

1 1.40 1.30 7
2 1.45 1.36 6
3 1.12 0.96 14
4 1.15 1.03 11
5 1.42 1.34 6
6 1.23 1.11 9
7 1.37 1.29 6
8 1.14 1.00 11
(to be continue)

53

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

Surface number µr /µs (EHL) µr /µs (PSD) Deviation(%)

9 1.19 1.07 10
10 1.37 1.29 6
11 1.21 1.09 9
12 1.42 1.37 4
13 1.34 1.36 1
14 1.42 1.36 4
15 1.29 1.17 9
16 1.09 0.93 15
17 1.34 1.24 7
18 1.25 1.12 10
19 1.29 1.19 8
20 1.33 1.22 8

4.4 Measured surface roughness

The surface roughness used in this section is the surface roughness measured from discs em-
ployed in friction experiments in Reference [140].

4.4.1 Friction prediction under a specific operating condition

In this subsection, the measured surface roughness dART is collected within an area of 1 cm2
and having a sampling interval of hx = h y = 3.653 µm [140]. This surface roughness is mea-
sured from a finished disc made of AISI 52100 steel with a 60 mm diameter and an 8 mm
thickness. The surface roughness is then corrected by removing large scales (through a high
pass-filtering operation), this corrected surface roughness and the original measured surface
roughness are shown in Figure 4.13. The root-mean-square of this surface roughness is R q =
0.038µm.
Since the area of high pressure happens in an area of Hertz contact zone, small square win-
dows with the length of 2a h are used to extract effective prediction areas (shown in Figure 4.14).
For a measured surface roughness, more effective prediction areas are obtained.

Figure 4.13: Measured surface roughness dART: (a) corrected surface roughness and (b) raw
surface roughness.

54

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

Figure 4.14: Effective prediction areas.

Table 4.4: Measured operating condition and lubricant parameters.

Parameter Value Units Parameter Value Units

w 10 N η0 0.7590 Pa · s
◦
ur 0.29 m/s Tr oom 22 C
Rx 9.525 mm SRR 25 %
E0 210 GPa ah 87.95 µm
α 2.0 × 10−8 Pa−1 ph 617.23 MPa

Figure 4.15: An extracted surface patch of dART: (a) surface roughness height of this surface
patch; (b) deformed surface patch; (c) pressure fluctuation of the surface patch.

55

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

For each extracted sub-surface, the friction prediction process described in subsection 4.2.1
is applied to predict its relative friction coefficient under a specific operating condition. Mea-
sured operating conditions and lubricant parameters are listed in Table 4.4. Under this specific
operating condition, Moes parameters M = 9.09 and L = 16.23, Figure 4.15 shows a deformed
surface geometry and pressure increase distribution of a surface patch, meanwhile its relative
friction coefficient is 1.03.

Figure 4.16: Extracted 529 surface patch (left) and their relative friction coefficient values
(right).

Figure 4.16 shows the extracted 529 surface patches (left) and their values of the relative
friction coefficient (right). For this measured surface roughness dART under this specific op-
erating condition, the relative friction coefficient is defined as the averaged value of that of
sub-surfaces. From this figure, one can observe that all relative friction values are smaller than
2 except for a point whose value reaches 4, hence this high value point will be removed. Figure
4.17 shows the corrected relative friction coefficients. At this time, the mean relative friction
coefficient is 1.09 and the corresponding standard deviation is 0.071.

Figure 4.17: Corrected relative friction coefficient values for 528 surface patches (left) and its
histogram (right).

Figure 4.18 shows the mean relative friction coefficient as a function of total number of
extracted surface patches with operating condition in Table 4.4. It can be observed that the
variation of mean relative friction for different the number of surface patch are small. For the
sake of saving computation time, the number of surface patch for the rest study is selected as

56

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

529.

Figure 4.18: Relative friction coefficient as a function of total number of surface patches under
the operating condition in Table 4.4.

4.4.2 Operating condition effects

The method to predict friction variations for measured rough surfaces has already described
in the previous subsection. In this subsection, effects of operating conditions on friction will
be investigated. Except for the rolling speed u r which is varying, the values of all the other
parameters are fixed as shown in Table 4.4. Equation C.1 shows that the pressure fluctuation
δp makes a significant contribution to the relative friction coefficient. Equation D.1 shows that
the amplitude deformation r r d of surface roughness is also an important variable. Therefore,
it is necessary to investigate δp and r r d . A surface patch shown extracted from the surface
roughness dART in figure 4.19 is extracted from the rough surface dART to study the variations
of δp and r r d when different operating conditions are employed (shown in table 4.5).

Table 4.5: Operating conditions of selected cases.

Cases u r (m/s) M L µr /µs

1 0.010 114.79 6.97 5.21
2 0.046 36.40 10.22 2.02
3 0.100 20.35 12.41 1.54
4 0.154 14.73 13.82 1.38
5 0.195 12.38 16.65 1.31

Figure 4.19: A extracted surface patch (left) and its average initial amplitude as a function as q.

57

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

Figure 4.20: The deformed surface patch shown in Figure 4.19 (left) and its pressure variations
for case 1.

Figure 4.21: The deformed surface patch shown in Figure 4.19 (left) and its pressure variations
for case 2.

Figure 4.22: The deformed surface patch shown in Figure 4.19 (left) and its pressure variations
for case 3.

Figure 4.23: The deformed surface patch shown in Figure 4.19 (left) and its pressure variations
for case 4.

58

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

Figure 4.24: The deformed surface patch shown in Figure 4.19 (left) and its pressure variations
for case 5.

Figures 4.20-4.24 show the deformed surface patches and their pressure variations as func-
tions of the wave vector q for all the operating conditions in table 4.5. The values at q = 0
in these figures are removed, because there are no deformations of sub-surfaces when q = 0.
According to observation, when the rolling speed u r increases, the value of the deformed sub-
surface increases while the value of the pressure increase decreases.

4.4.3 Friction curves for measured surface roughness

Friction variations for two measured surface roughness are plotted in this subsection. For the
roughness dART (shown in Figure 4.13), the variations of the Moes load parameter M and ma-
terial parameter L using geometry and lubricant parameter in Table 4.4 are shown in Figure
4.25.

Figure 4.25: Moes parameters M and L as a function of u r for the roughness dART.

Figure 4.26 shows variations of the relative friction coefficient employing the operating con-
ditions shown in Figure 4.25. It can be observed from this figure that as rolling speed u r in-
creases, the relative friction coefficient decreases. This decrease trend can also be found in
Reference [140] Figure 6.1 (b). According to Hamrock-Dowson central film thickness equa-
tion [120], h c is a function of rolling speed u r , hence Figure 4.26 can be re-plotted as a function
of the "Λ ratio" (shown in Figure 4.27).

59

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

Figure 4.26: The relative friction coefficient as a function of u r for the roughness dART.

Figure 4.27: The relative friction coefficient as a function of "Λ ratio" for the roughness dART.

Figure 4.28: Surface roughness dARTEb (left) and the top view of this roughness (right).

60

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

Another surface roughness dARTEb (shown in Figure 4.28) is collected from a finished and
roughening steel disc, with its root-mean-square R q = 0.132µm. The size of this disc is the same
as the one for surface roughness dART.

Figure 4.29: Moes parameters M and L as a function of the rolling speed for the surface rough-
ness dARTEb.

Figure 4.30: The relative friction coefficient as a function of u r for the surface roughness dAR-
TEb.

For the surface roughness dARTEb, operating conditions can be seen in Figure 4.29. And
friction variation curves for this roughness is shown in Figure 4.30 and Figure 4.31. The classical
Stribeck curve has a minimum friction point, which is often regarded as the transition point
from the mixed lubrication to hydrodynamic lubrication regime. However, for friction curves

61

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

shown in these two figures, minimum friction points do not occur. Reference [147] defined the
transition position as the intersection between two tangent lines to the friction curve (shown
as blue crosses in Figure 4.27 and 4.31). For conformal contact, the transition from the mixed
lubrication to full-film lubrication regime occurs around Λ = 3.0. For non-conformal contact
cases considered in this subsection, the onset of the "mixed regime" at hc/R q varying from 4 to
5. At the same time, two different surface roughness show different friction behaviour(shown
in Figure 4.32).

Figure 4.31: The relative friction coefficient as a function of the "Λ ratio" for the surface rough-
ness dARTEb.

Figure 4.32: The relative friction coefficient as a function of the rolling speed u r for the surface
roughness dART and dARTEb, respectively.

62

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES

4.5 Conclusion
This chapter proposes a rapid analytical method using the power spectral density to predict
friction increase which is due to non-linear viscosity variations from pressure variations caused
by surface roughness deformation. This method is validated for an artificial surface rough-
ness by comparing the results with those of a full numerical simulation. A good agreement is
found between the full numerical simulation and the PSD prediction. Then the PSD predic-
tion method is also employed to analyse friction variations for a measured surface roughness.
Prediction results show that as the "Λ ratio" increases, the friction decreases. The minimum
friction point describing the transition from the mixed lubrication to the full-film lubrication
regime does not occur for high pressure cases. The same result has also been found in experi-
mental work [140, 147].

63

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Chapter 5

Conclusion and perspective

Contents
5.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5.2 Perspective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5.1 Conclusion
Proper lubrication and low friction play essential roles in energy conservation. The Stribeck
curve is a good tool to describe the friction variation throughout the entire lubrication regime.
However, most of the work to obtain this friction curve was related to low pressure condition
(conformal contact), and the work was mainly experimental. The motivation for conduct-
ing this work is to determine the frictional behaviour under high pressure conditions (non-
conformal contact) when the surface roughness is taken into consideration. The main conclu-
sions are listed below:

(i) Current MultiGrid codes show good efficiency to solve lubricated contact problems, but
are not sufficiently robust to treat the rough surface problem in a general way. This is be-
cause for very rough surfaces, large variations of the coefficient ρh 3 /η in Reynolds equa-
tion occur on a small scale. Alcouffe et al [66] proposed an efficient way to restore the
performance by constructing the coarse grid operator and the intergrid transfers. An new
transient EHL code is modified by implementing the coarse grid construction method
into the existing EHL MultiGrid code. The new code shows better performance in solving
the rough contact problem.

(ii) The Stribeck curve shows that the "Λ ratio" is a suitable parameter to plot friction vari-
ation. According to the Amplitude Reduction Theory [53, 58], under very high pressure,
the surface roughness will deform and this deformation depends on the operating con-
ditions as well as on the surface roughness parameters. This means that the old "Λ ratio"
is not a proper parameter. In this work, using harmonic surface waviness, an elaborate
scaling parameter θ2 ff(r ) including the old "Λ ratio", operating conditions and surface
anisotropy parameters is found. Using this new parameters, all simulation results can be
unified into a single curve and a curve-fitting equation is obtained.

(iii) A rapid analytical prediction method using the power spectral density is proposed in this
work to predict friction increase for a complex surface roughness. This friction increase is
due to non-linear viscosity variations from pressure variations caused by surface rough-
ness deformation. This method is validated by the comparison between predictions and
full numerical simulations when an artificial surface roughness is employed. Then this
rapid prediction method is applied to analyse measured surface roughness.

64

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 CHAPTER 5. CONCLUSION AND PERSPECTIVE

5.2 Perspective
The work carried out in this thesis offers a useful tool to detect friction variations under high
pressure conditions. However, what we presented in this work is related to Newtonian lubri-
cant, sufficient lubrication and isothermal model. Future work can be suggested as follows:

• In current work, a small slide-to-roll ratio is used to prediction frictional force and a
unique curve is obtained. Big slide-to-roll ratio value should be applied to detect fric-
tion variations in a further study.

• The EHL model used in this work is isothermal and Newtonian. Thermal effects and
non-Newtonian lubricant can be considered in future study.

• The work studied in this thesis is numerical, more relevant experiments should be con-
ducted to give comparable results.

65

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
Cette thèse est accessible à l'adresse : [Link]
© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Appendix A

Construction of the coarse grid operator

The coarse grid operator is:

L H = J hH L h I H
h
(A.1)
In order to get the coarse grid operator, orthogonal basis vectors are needed. The first step is to
h H
compute I H 〈e iC , jC ,k
〉:

h H h h h h h h
IH 〈e iC , jC ,k 〉 = C i F, j F,k e i F, j F,k + NWi F, j F,k e i F −1, j F +1,k + SWi F, j F,k e i F −1, j F −1,k

+ SE ihF, j F,k e ihF +1, j F −1,k + N E ihF, j F,k e ihF +1, j F +1,k
(A.2)
+ NihF, j F,k e ihF, j F +1,k + S ihF, j F,k e ihF, j F −1,k
+ WihF, j F,k e ihF −1, j F,k + E ihF, j F,k e ihF +1, j F,k

Then, substituting Equation A.2 in L h 〈I H

h H
e iC , jC ,k 〉 reads:

L h 〈I H
h H
e iC , jC ,k 〉 = L h [C ihF, j F,k e ihF, j F,k + NWihF, j F,k e ihF −1, j F +1,k + SWihF, j F,k e ihF −1, j F −1,k
+ SE ihF, j F,k e ihF +1, j F −1,k + N E ihF, j F,k e ihF +1, j F +1,k
(A.3)
+ NihF, j F,k e ihF, j F +1,k + S ihF, j F,k e ihF, j F −1,k
+ WihF, j F,k e ihF −1, j F,k + E ihF, j F,k e ihF +1, j F,k ]

in which the influence of L h need to act on every term, as illustrated in Figure A.1 for 3 terms:

• The influence of L h on the central point (i , j ) (shown as Figure A.1 (a)):

L h 〈C ihF, j F,k e ihF, j F,k 〉 = C ihF, j F,k L h 〈e ihF, j F,k 〉

= C ihF, j F,k [L hc (i F, j F, k)e ihF, j F,k
+ L hnw (i F, j F, k)e ihF −1, j F +1,k + L hsw (i F, j F, k)e ihF −1, j F −1,k
+ L hse (i F, j F, k)e ihF +1, j F −1,k + L hne (i F, j F, k)e ihF +1, j F +1,k
+ L hn (i F, j F, k)e ihF, j F +1,k + L hw (i F, j F, k)e ihF −1, j F,k
+ L hs (i F, j F, k)e ihF, j F −1,k + L he (i F, j F, k)e ihF +1, j F,k ]

• The influence of L h on the east point (i + 1, j ) (shown as Figure A.1 (b)):

L h 〈E ihF, j F,k e ihF +1, j F,k 〉 = E ihF, j F,k L h 〈e ihF +1, j F,k 〉
= E ihF, j F,k [L hc (i F + 1, j F, k)e ihF +1, j F,k
+ L hnw (i F + 1, j F, k)e ihF, j F +1,k + L hsw (i F + 1, j F, k)e ihF, j F −1,k
+ L hse (i F + 1, j F, k)e ihF +2, j F −1,k + L hne (i F + 1, j F, k)e ihF +2, j F +1,k
+ L hn (i F + 1, j F, k)e ihF +1, j F +1,k + L hw (i F + 1, j F, k)e ihF, j F,k
+ L hs (i F + 1, j F, k)e ihF +1, j F −1,k + L he (i F + 1, j F, k)e ihF +2, j F,k ]

66

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 APPENDIX A. CONSTRUCTION OF THE COARSE GRID OPERATOR

• The influence of L h on the north-east point (i + 1, j + 1) (shown as Figure A.1 (c)):

L h 〈N E ihF, j F,k e ihF +1, j F +1,k 〉 = N E ihF, j F,k L h 〈e ihF +1, j F +1,k 〉
= N E ihF, j F,k [L hc (i F + 1, j F + 1, k)e ihF +1, j F +1,k
+ L hnw (i F + 1, j F + 1, k)e ihF, j F +2,k
+ L hsw (i F + 1, j F + 1, k)e ihF, j F,k
+ L hse (i F + 1, j F + 1, k)e ihF +2, j F,k
+ L hne (i F + 1, j F + 1, k)e ihF +2, j F +2,k
+ L hn (i F + 1, j F + 1, k)e ihF +1, j F +2,k
+ L hw (i F + 1, j F + 1, k)e ihF, j F +1,k
+ L hs (i F + 1, j F + 1, k)e ihF +1, j F,k
+ L he (i F + 1, j F + 1, k)e ihF +2, j F +1,k ]

Figure A.1: Influences of the coarse grid operator L h on central point, east point and north-east
point.

From the three equations above one observes that the operator L h influences fifteen points
around point (i , j ) in total. By that analogy, with respect to nine terms in Equation A.3 twenty
five points around point (i , j ) are affected by the operator L h . Coefficients of 25 points around
(i , j ) is denoted as:
IC (i F − 2, j F + 2, k) IC (i F − 1, j F + 2, k) IC (i F, j F + 2, k) IC (i F + 1, j F + 2, k) IC (i F + 2, j F + 2, k)
 
 
 
 IC (i F − 2, j F + 1, k) IC (i F − 1, j F + 1, k) IC (i F, j F + 1, k) IC (i F + 1, j F + 1, k) IC (i F + 2, j F + 1, k)
 
 
 
 IC (i F − 2, j F, k) IC (i F − 1, j F, k) IC (i F, j F, k) IC (i F + 1, j F, k) IC (i F + 2, j F, k)  (A.4)
 
 
 
 
 IC (i F − 2, j F − 1, k) IC (i F − 1, j F − 1, k) IC (i F, j F − 1, k) IC (i F + 1, j F − 1, k) IC (i F + 2, j F − 1, k)
 
 
IC (i F − 2, j F − 2, k) IC (i F − 1, j F − 2, k) IC (i F, j F − 2, k) IC (i F + 1, j F − 2, k) IC (i F + 2, j F − 2, k)

where:

IC (i F − 2, j F + 2, k) = NWihF, j F,k L hnw (i F − 1, j F + 1, k)

67

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 APPENDIX A. CONSTRUCTION OF THE COARSE GRID OPERATOR

IC (i F − 2, j F + 1, k) = NWihF, j F,k L hw (i F − 1, j F + 1, k)

IC (i F − 2, j F, k) = NWihF, j F,k L hsw (i F − 1, j F + 1, k)

+ SWihF, j F,k L hnw (i F − 1, j F − 1, k) + WihF, j F,k L hsw (i F − 1, j F, k)

IC (i F − 2, j F − 1, k) = SWihF, j F,k L hw (i F − 1, j F − 1, k) + WihF, j F,k L hsw (i F − 1, j F, k)

IC (i F − 2, j F − 2, k) = SWihF, j F,k L hsw (i F − 1, j F − 1, k)

IC (i F − 1, j F + 2, k) = NWihF, j F,k L hn (i F − 1, j F + 1, k) + NihF, j F,k L hnw (i F, j F + 1, k)

IC (i F − 1, j F + 1, k) = C ihF, j F,k L hnw (i F, j F, k) + NWihF, j F,k L hc (i F − 1, j F + 1, k)

+ NihF, j F,k L hw (i F, j F + 1, k) + WihF, j F,k L hn (i F − 1, j F, k)

IC (i F − 1, j F, k) = C ihF, j F,k L hw (i F, j F, k) + NWihF, j F,k L hs (i F − 1, j F + 1, k)

+ SWihF, j F,k L hn (i F − 1, j F − 1, k) + NihF, j F,k L hsw (i F, j F + 1, k)
+ S ihF, j F,k L hnw (i F, j F − 1, k) + WihF, j F,k L hc (i F − 1, j F, k)

IC (i F − 1, j F − 1, k) = C ihF, j F,k L hsw (i F, j F, k) + SWihF, j F,k L hc (i F − 1, j F − 1, k)

+ S ihF, j F,k L hw (i F, j F − 1, k) + WihF, j F,k L hs (i F − 1, j F, k)

IC (i F − 1, j F − 2, k) = SWihF, j F,k L hs (i F − 1, j F − 1, k) + S ihF, j F,k L hsw (i F, j F − 1, k)

IC (i F, j F + 2, k) = NWihF, j F,k L hne (i F − 1, j F + 1, k) + N E ihF, j F,k L hw (i F, j F − 1, k)

+ NihF, j F,k L hn (i F, j F + 1, k)

IC (i F, j F + 1, k) = C ihF, j F,k L hn (i F, j F, k) + NWihF, j F,k L he (i F − 1, j F + 1, k)

+ N E ihF, j F,k L hw (i F + 1, j F + 1, k) + NihF, j F,k L hc (i F, j F + 1, k)
+ WihF, j F,k L hne (i F − 1, j F, k) + E ihF, j F,k L hnw (i F + 1, j F, k)

IC (i F, j F, k) = C ihF, j F,k L hc (i F, j F, k) + NWihF, j F,k L hse (i F − 1, j F + 1, k)

+ SWihF, j F,k L hne (i F − 1, j F − 1, k)
+ SE ihF, j F,k L hnw (i F + 1, j F − 1, k) + N E ihF, j F,k L hsw (i F + 1, j F + 1, k)
+ NihF, j F,k L hs (i F, j F + 1, k) + S ihF, j F,k L hn (i F, j F − 1, k)
+ WihF, j F,k L he (i F − 1, j F, k) + E ihF, j F,k L hw (i F + 1, j F, k)

68

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 APPENDIX A. CONSTRUCTION OF THE COARSE GRID OPERATOR

The last step is to compute J hH 〈L h I H

h H
e iC , jC ,k 〉. The way to extend the restriction coefficients
on 25 points in Equation A.4 is the same as that of computing L h 〈I H
h H
e iC , jC ,k 〉. There is no need
to do a restriction for 25 points. As was shown in Figure A.2, besides coincidental points (blue
points) themselves, 3 points contribute to each corner point (shown as Figure A.2 (a)), 5 points
influence each middle point (shown as Figure A.2 (b)) and 8 points impact on central point
(shown as Figure A.2 (c)).

Figure A.2: Influences of the injection operator J hH on nine coincidental points (blue points).

Hence, the coarse grid operator whose stencil is:

 H
L nw (iC , jC , k) L nH (iC , jC , k) L ne
H

(iC , jC , k)
 
 
H  H H H
L =  L w (iC , jC , k) L c (iC , jC , k) L e (iC , jC , k)  (A.5)

 
 
H H H
L sw (iC , jC , k) L s (iC , jC , k) L se (i c, jC , k)

in which:
H
L nw (iC , jC , k) = IC (i F − 2, j F + 2, k)c ihF −2, j F +2,k + IC (i F − 1, j F + 2, k)w ihF −1, j F +2,k
+ IC (i F − 2, j F + 1, k)n ihF −2, j F +1,k + IC (i F − 1, j F + 1, k)nw ihF −1, j F +1,k

H
L sw (iC , jC , k) = IC (i F − 2, j F − 2, k)c ihF −2, j F −2,k + IC (i F − 2, j F − 1, k)s ihF −2, j F −1,k
+ IC (i F − 1, j F − 2, k)w ihF −1, j F −2,k + IC (i F − 1, j F − 1, k)sw ihF −1, j F −1,k

H
L se (iC , jC , k) = IC (i F + 2, j F − 2, k)c ihF +2, j F −2,k + IC (i F + 1, j F − 2, k)e ihF +1, j F −2,k
+ IC (i F + 2, j F − 1, k)s ihF +2, j F −1,k + IC (i F + 1, j F − 1, k)se ihF +1, j F −1,k

H
L ne (iC , jC , k) = IC (i F + 2, j F + 2, k)c ihF +2, j F +2,k + IC (i F + 1, j F + 2, k)e ihF +1, j F +2,k
+ IC (i F + 2, j F + 1, k)n ihF +2, j F +1,k + IC (i F + 1, j F + 1, k)ne ihF +1, j F +1,k

69

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 APPENDIX A. CONSTRUCTION OF THE COARSE GRID OPERATOR

L nH (iC , jC , k) = IC (i F, j F + 2, k)c ihF, j F +2,k + IC (i F − 1, j F + 2, k)e ihF −1, j F +2,k

+ IC (i F + 1, j F + 2, k)w ihF +1, j F +2,k + IC (i F, j F + 1, k)n ihF, j F +1,k
+ IC (i F − 1, j F + 1, k)ne ihF −1, j F +1,k + IC (i F + 1, j F + 1, k)nw ihF +1, j F +1,k

L sH (iC , jC , k) = IC (i F, j F − 2, k)c ihF, j F −2,k + IC (i F − 1, j F − 2, k)e ihF −1, j F −2,k

+ IC (i F + 1, j F − 2, k)w ihF +1, j F −2,k + IC (i F, j F − 1, k)s ihF, j F −1,k
+ IC (i F − 1, j F − 1, k)se ihF −1, j F −1,k + IC (i F + 1, j F − 1, k)sw ihF +1, j F −1,k

H
Lw (iC , jC , k) = IC (i F − 2, j F, k)c ihF −2, j F,k + IC (i F − 2, j F + 1, k)s ihF −2, j F +1,k
+ IC (i F − 2, j F − 1, k)n ihF −2, j F −1,k + IC (i F − 1, j F, k)w ihF −1, j F,k
+ IC (i F − 1, j F + 1, k)sw ihF −1, j F +1,k + IC (i F − 1, j F − 1, k)nw ihF −1, j F −1,k

L eH (iC , jC , k) = IC (i F + 2, j F, k)c ihF +2, j F,k + IC (i F + 2, j F + 1, k)s ihF +2, j F +1,k

+ IC (i F + 2, j F − 1, k)n ihF +2, j F −1,k + IC (i F + 1, j F, k)e ihF +1, j F,k
+ IC (i F + 1, j F + 1, k)se ihF +1, j F +1,k + IC (i F + 1, j F − 1, k)ne ihF +1, j F −1,k

L cH (iC , jC , k) = IC (i F, j F, k)c ihF, j F,k + IC (i F, j F + 1, k)s ihF, j F +1,k + IC (i F, j F − 1, k)n ihF, j F −1,k
+ IC (i F − 1, j F, k)e ihF −1, j F,k + IC (i F + 1, j F, k)w ihF +1, j F,k
+ IC (i F − 1, j F + 1, k)se ihF −1, j F +1,k + IC (i F − 1, j F − 1, k)ne ihF −1, j F −1,k
+ IC (i F + 1, j F − 1, k)nw ihF +1, j F −1,k + IC (i F + 1, j F + 1, k)sw ihF +1, j F +1,k

70

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Appendix B

Derivation of matrix A j for line relaxation

According to the discrete Reynolds equation in Chapter 2, the dynamic residual r il, j ,k is
defined as:
r il, j ,k = p f il, j ,k −(L lc P il , j ,k + L ln P il , j +1,k + L lw P il −1, j ,k + L ls P il , j −1,k + L le P il +1, j ,k
+ L lnw P il −1, j +1,k + L lsw P il −1, j −1,k + L lse P il +1, j −1,k + L lne P il +1, j +1,k )
+ h l (1.5ρ̄ li , j ,k Hil, j ,k − 2.0ρ̄ li −1, j ,k Hil−1, j ,k + 0.5ρ̄ li −2, j ,k Hil−2, j ,k ) (B.1)

(h l )2
+ (1.5ρ̄ i , j ,k Hi , j ,k )
ht
and the left hand side equation is:

L li , j ,k 〈P l 〉 =(L lc P il , j ,k + L ln P il , j +1,k + L lw P il −1, j ,k + L ls P il , j −1,k + L le P il +1, j ,k

+ L lnw P il −1, j +1,k + L lsw P il −1, j −1,k + L lse P il +1, j −1,k + L lne P il +1, j +1,k )
− h l (1.5ρ̄ li , j ,k Hil, j ,k − 2.0ρ̄ li −1, j ,k Hil−1, j ,k + 0.5ρ̄ li −2, j ,k Hil−2, j ,k ) (B.2)

(h l )2
− (1.5ρ̄ i , j ,k Hi , j ,k ).
ht
Reference [10] recommend that the switch parameter ξllimit = 0.3 is a good choice to have an
efficient smooth performance.

B.0.1 Gauss-Seidel line relaxation

When the local coefficient ξl satisfies the below condition, the Gauss-Seidel line relaxation is
applied.
ξli ±1/2, j ,k ξli , j ±1/2,k
l
> ξlimit and > ξllimit . (B.3)
(h l )2 (h l )2
j
Then the matrix A i ,m is given by:

∂L li , j ,k 〈P l 〉
j
A i ,m = ( l
)P h =P̃ h (B.4)
∂P m, j ,k

for 0 < m < n x and 0 < i < n x .

j
The matrix A i ,m have different expressions for different conditions:
• if |i − m| > 1:
j
A i ,m = − h l (1.5ρ̄ li , j ,k K |il l−m|,0 − 2.0ρ̄ li −1, j ,k K |il l−m−1|,0 + 0.5ρ̄ li −2, j ,k K |il l−m−2|,0 )
(h l )2 (B.5)
− (1.5ρ̄ i , j ,k K |il l−m|,0 )
• if i = m: ht
j
A i ,i = L lc − h l (1.5ρ̄ li , j ,k K 0,0
ll
− 2.0ρ̄ li −1, j ,k K 1,0
ll
+ 0.5ρ̄ li −2, j ,k K 2,0
ll
)
(h l )2 (B.6)
ll
− (1.5ρ̄ i , j ,k K 0,0 )
ht
71

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 APPENDIX B. DERIVATION OF MATRIX A J FOR LINE RELAXATION

• if i > 1:
j
A i ,i −1 = L lw − h l (1.5ρ̄ li , j ,k K 1,0
ll
− 2.0ρ̄ li −1, j ,k K 0,0
ll
+ 0.5ρ̄ li −2, j ,k K 1,0
ll
)
(h l )2 (B.7)
ll
− (1.5ρ̄ i , j ,k K 1,0 )
ht
• if i < n x − 1:
j
A i ,i +1 = L le − h l (1.5ρ̄ li , j ,k K 1,0
ll
− 2.0ρ̄ li −1, j ,k K 2,0
ll
+ 0.5ρ̄ li −2, j ,k K 3,0
ll
)
(h l )2 (B.8)
ll
− (1.5ρ̄ i , j ,k K 1,0 ).
ht

B.0.2 Jacobi distributive line relaxation

When the local coefficient ξl can not satisfy the condition Equation B.3, the Jacobi distributive
line relaxation is used.
j
Subsequently, the matrix A i ,m is:

∂L li , j ,k 〈P l 〉 l l l l l l l l
1 ∂L i , j ,k 〈P 〉 ∂L i , j ,k 〈P 〉 ∂L i , j ,k 〈P 〉 ∂L i , j ,k 〈P 〉
j
A i ,m =[ l
− ( + + + )]P h =P̃ h (B.9)
∂P m, 4 ∂P l l
∂P m−1, l
∂P m, l
∂P m,
j ,k m+1, j ,k j ,k j +1,k j −1,k

for 0 < m < n x and 0 < i < n x . For convenience, introducing a parameter ∆K m, j to simply the
j
expression of A i ,m :

ll ll 1 ll ll ll ll
∆K m,n = K m,n − (K m−1,n + K m+1,n + K m,n+1 + K m,n−1 ). (B.10)
4
j
For different i values, the matrix of A i ,m are represented as:

• for |i − m| > 2:
j
A i ,m = − h l (1.5ρ̄ li , j ,k ∆K |il l−m|,0 − 2.0ρ̄ li −1, j ,k ∆K |il l−m−1|,0 + 0.5ρ̄ li −2, j ,k ∆K |il l−m−2|,0 )
(h l )2 (B.11)
− (1.5ρ̄ i , j ,k ∆K |il l−m|,0 )
ht
• for i = m:
j 1
A i ,i = L lc − (L le + L lw + L ln + L ls )
4
− h l (1.5ρ̄ li , j ,k ∆K 0,0
ll
− 2.0ρ̄ li −1, j ,k ∆K 1,0
ll
+ 0.5ρ̄ li −2, j ,k ∆K 2,0
ll
) (B.12)
(h l )2 ll
− (1.5ρ̄ i , j ,k ∆K 0,0 )
ht
• for i > 2:
j 1
A i ,i −2 = − L lw − h l (1.5ρ̄ li , j ,k ∆K 2,0
ll
− 2.0ρ̄ li −1, j ,k ∆K 1,0
ll
+ 0.5ρ̄ li −2, j ,k ∆K 0,0
ll
)
4
(B.13)
(h l )2 ll
− (1.5ρ̄ i , j ,k ∆K 2,0 )
ht

• for i > 1:
j 1
A i ,i −1 = L lw − (L lc + L ln w + L ls w)
4
− h l (1.5ρ̄ li , j ,k ∆K 1,0
ll
− 2.0ρ̄ li −1, j ,k ∆K 0,0
ll
+ 0.5ρ̄ li −2, j ,k ∆K 1,0
ll
) (B.14)
l 2
(h ) ll
− (1.5ρ̄ i , j ,k ∆K 1,0 )
ht

72

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 APPENDIX B. DERIVATION OF MATRIX A J FOR LINE RELAXATION

• for i < n x − 1:

j 1
A i ,i +1 = L le − (L lc + L ln e + L ls e)
4
− h l (1.5ρ̄ li , j ,k ∆K 1,0
ll
− 2.0ρ̄ li −1, j ,k ∆K 2,0
ll
+ 0.5ρ̄ li −2, j ,k ∆K 3,0
ll
) (B.15)
(h l )2 ll
− (1.5ρ̄ i , j ,k ∆K 1,0 )
ht

• for i < n x − 2:

j 1
A i ,i +2 = − L le − h l (1.5ρ̄ li , j ,k ∆K 2,0
ll
− 2.0ρ̄ li −1, j ,k ∆K 3,0
ll
+ 0.5ρ̄ li −2, j ,k ∆K 4,0
ll
)
4
(B.16)
(h l )2 ll
− (1.5ρ̄ i , j ,k ∆K 2,0 ).
ht

73

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Appendix C

Derivation of the scaling parameter θ2

According to the Barus [115] viscosity-pressure equation, the shear stress ratio can be ap-
proximated as:

τr (ᾱ · ∆P )2 (ᾱ · ∆P )3 (ᾱ · ∆P )4

≈ e ᾱ·∆P = 1 + ᾱ · ∆P + + + + ... (C.1)
τs 2! 3! 4!
where the pressure increase ∆P is the following function of the deformation [117], i.e.:

π2 Ai Ad
∆P = (1 − ).
2λ/a h Ai

Using a first order approximation of the dimensionless pressure increase ∆P , ᾱ·∆P reduces to:

π2 Ai L 3M 1/3 π2 H c H c −1
ᾱ · ∆P ≈ ᾱ ≈[ ( ) ] ( ) (C.2)
2λ/a h π 2 2λ/a h Ai

where ᾱ is expressed as ᾱ = (L/π)(3M /2)1/3 . Defining HcD the dimensionless film thickness film
thickness value using the well-known Hamrock-Dowson Equation [120], HcD = 1.69G 0.53U 0.67W2−0.067 (1−
0.61exp(−0.73k)) with k = 1.03 for circle contact (i.e. HcD = 1.2L 0.53U 0.49 M −0.067 ). Now the di-
mensionless central film thickness Hc can be rewritten as:

R x2
Hc = · HcD (C.3)
a h2

in which R x2 /a h2 is expressed as R x2 /a h2 = (3/2)−2/3 M −2/3U −1/2 . Substituting Equation C.3 into

Equation C.2 gives:
ᾱ · ∆P ≈ 1.6467[L −1.03 M −0.4 (λ/a h )1 (H c/Ai )]−1 . (C.4)
Applying a second order approximation of ∆P , ᾱ∆P yields:

π2 Ai Ad
ᾱ · ∆P ≈ (1 − ) ≈ 0.24[L −1.03 M −0.1 (λ/a h )0 (H c/Ai )]−1 (C.5)
2λ/a h Ai

where Ad /Ai ≈ 1 − 0.15∇2 ≈ 1 − 0.15(λ/a h )(M /L)0.5 .

Observing Equation C.4 and Equation C.5, the exponent for the parameter M , L and λ/a h
are summarized in Table C.1. Hence the expression of the θ2 parameter using M 0.33 , L −1.1 and
(λ/a h )0.67 obtained by curve-fitting, has coefficients that fall in the range outlined above.

Table C.1: Range of the exponent for each parameter.

Parameter Range
M −0.10 ∼ 0.40
L −1.03 ∼ −1.53
λ/a h 0∼1
H c/Ai 1

74

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Appendix D

The relation between the elastic deformation and corre-

sponding pressure for 2D wavy surfaces

Reference [117] gives the relation between the elastic deformation and pressure for 1D wavy
surfaces is:
πE 0 ∆
δp = (D.1)
λ
and for 2D wavy surfaces, which consists of two separate components:
πE 0 ∆x πE 0 ∆ y
δp = + . (D.2)
λx λy
For the surface waviness expressed as Eq.3.1, there is no explicit equation to describe this rela-
tion. Hence, it is necessary to give such relation. Assuming a pressure distribution:
p 2πx 2πy
p(x, y) = A i cos( )cos( ) (D.3)
λx λy
p
with x ∈ [−4, 4] and y ∈ [−4, 4]. Where A i is the initial amplitude of the pressure distribution.
The equation to compute the corresponding elastic deformation is:
2 p(x 0 , y 0 )
Ï
D(x, y) = d x 0d y 0. (D.4)
πE 0
p
0 2 0
(x − x ) + (y − y ) 2

Figure D.1: Pressure distribution and the corresponding elastic deformation.

75

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 APPENDIX D. THE RELATION BETWEEN THE ELASTIC DEFORMATION AND
CORRESPONDING PRESSURE FOR 2D WAVY SURFACES
Figure D.1 shows that the elastic deformation due to the pressure shown as Figure D.1 (a)
has the same harmonic waviness. Thus the elastic deformation Equation D.4 can be re-written
as:
2πx 2πy
D(x, y) = A D
d cos( )cos( ) (D.5)
λx λy

in which A D
d
is the amplitude of the elastic deformation.

Figure D.2: Amplitude of the elastic deformation A D d

as a function of initial pressure amplitude
p
A i for the following cases: (a) isotropic, (b) purely transverse, (c) purely longitudinal.

Figure D.2 shows the relation between the amplitude of the elastic deformation A Dd
and the
p
initial pressure amplitude A i for three r values. From this figure, it is can be found that the
p
deformation amplitude A D d
is proportional to the initial pressure amplitude A i :

p
AD 0
d πE ∝ A i . (D.6)

Figure D.3: Amplitude of the elastic deformation A D d

as a function of wavelength λ (λ =
min(λx , λ y )) for the following cases: (a) isotropic, (b) purely transverse, (c) purely longitudi-
nal.

76

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 APPENDIX D. THE RELATION BETWEEN THE ELASTIC DEFORMATION AND
CORRESPONDING PRESSURE FOR 2D WAVY SURFACES
Figure D.3 shows the relation between the amplitude of the elastic deformation A D d
and
the wavelength λ for three r values. From this figure, it is can be found that the deformation
amplitude A D
d
is proportional to the the wavelength λ:

AD 0
d πE ∝ λ. (D.7)

Figure D.4 shows the relation between the amplitude of the elastic deformation A Dd
and the
anisotropy parameter r . From this figure, it is can be found that the deformation amplitude A D
d
is proportional to a curve-fitting function k(r ):

AD 0
d πE ∝ k(r ) (D.8)

where (
−0.3r 1.7 + 1.0, if 0 É r É 1
k(r ) = −1.7
(D.9)
−0.3r + 1.0, r > 1
with r = λx /λ y .

Figure D.4: Amplitude of the elastic deformation A D

d
as a function of anisotropy parameter
r = λx /λ y .

Combing Equation D.6, Equation D.7, Equation D.8 and Equation D.9, the relation between
p
the initial amplitude of pressure A i and the corresponding elastic deformation A D
d
is:

p πE 0
Ai = AD (D.10)
2k(r )λ d

with λ = min(λx , λ y ) and

(
−0.3r 1.7 + 1.0, if 0 É r É 1
k(r ) = −1.7
−0.3r + 1.0, r > 1
where r = λx /λ y .

77

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Acknowledgments

To all those people who have supported me in carrying out my PhD study.
Firstly, I would like to express my sincere gratitude to my supervisor Professor Antonius Lu-
brecht for his guidance, encouragement and support of this research work. I would also like to
thank my supervisor for opportunities I have been given to participate in scientific conferences
as well as two months academic exchanges. Meanwhile, I am very grateful to my supervisor for
caring my life in Lyon. During these three years, my supervisor has taught me more than this
specific subject, he has taught me how to think more logically and how to solve problems more
efficiently. Without his persistent help, this work would not have materialized.
I would also like to extend my thanks to Professor Cornelis Venner and his group for their
warmly welcome and help in the MultiGrid method studying and the first manuscript writing
during the period I have visited the University of Twente. I also want to thank my friends:
Mr. Haichao Liu, Mr. Balan Ramani and Mrs. Xiaozhen He for their help when I was in the
University of Twente. Special thanks to Mr. Binbin Zhang for taking care of me in Enschede.
I would like to show my appreciations to my colleagues: Dr. Romain Bugnicourt, Mr. Martin
Denni, Ms. Simona Dahdah and Ms. Marjolaine Gonon-Caux. With them, we built a pleasant
and friendly working environment. I would like to offer my special thanks to Dr. Marie-Pierre
Noutary for teaching me the coarse grid construction method. I also thankful to Dr. Nans
Biboulet for his time and advice of my work. I would like to thank the members of the secretary
group of LaMCoS for their help over the last three years.
Thanks to my following dear friends with whom I spent a happy time: Ms. Jing Lu, Ms. Pei
Niu, Mr. Minhu Lv, Ms. Jingya Yuan, Mr. Minyang Lou et al. I would especially like to thank Dr.
Hanfeng Gu and Dr. Teng Zhang, for their selfless help when I first came to Lyon.
This research work is financially supported by the China Scholarship Council (CSC).
Finally, my deepest appreciation goes to my parents: Mr. Ming Zhang and Mrs. Shulan
Zeng. Meanwhile, I also thank my boyfriend for his support and waiting.
Thank you for all of you again.

78

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
Cette thèse est accessible à l'adresse : [Link]
© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 Bibliography

[1] K. Holmberg, P. Andersson, and A. Erdemir, “Global energy consumption due to friction
in passenger cars,” Tribology International, vol. 47, pp. 221–234, 2012.

[2] K. Holmberg, R. Siilasto, T. Laitinen, P. Andersson, and A. Jäsberg, “Global energy con-
sumption due to friction in paper machines,” Tribology International, vol. 62, pp. 58–77,
2013.

[3] K. Holmberg, P. Andersson, N.-O. Nylund, K. Mäkelä, and A. Erdemir, “Global energy con-
sumption due to friction in trucks and buses,” Tribology International, vol. 78, pp. 94–
114, 2014.

[4] K. Holmberg, P. Kivikytö-Reponen, P. Härkisaari, K. Valtonen, and A. Erdemir, “Global

energy consumption due to friction and wear in the mining industry,” Tribology Interna-
tional, vol. 115, pp. 116–139, 2017.

[5] P. M. Lee and R. Carpick, “Tribological opportunities for enhancing america’s energy ef-
ficiency.,” Tech. Rep. DE-AR00282, A Report to the Advanced Research Project Agency-
Energy (ARPA-E) at the U.S. Department of Energy, 2017.

[6] K. Holmberg and A. Erdemir, “Global impact of friction on energy consumption, econ-
omy and environment,” FME Trans, vol. 43, no. 3, pp. 181–5, 2015.

[7] K. Holmberg and A. Erdemir, “Influence of tribology on global energy consumption, costs
and emissions,” Friction, vol. 5, no. 3, pp. 263–284, 2017.

[8] V. W. Wong and S. C. Tung, “Overview of automotive engine friction and reduction
trends–effects of surface, material, and lubricant-additive technologies,” Friction, vol. 4,
no. 1, pp. 1–28, 2016.

[9] P. M. Lugt and G. E. Morales-Espejel, “A review of elasto-hydrodynamic lubrication the-

ory,” Tribology Transactions, vol. 54, no. 3, pp. 470–496, 2011.

[10] C. H. Venner and A. A. Lubrecht, Multi-Level Methods in Lubrication. Elsevier, 2000.

[11] D. Dowson and G. R. Higginson, “A numerical solution to the elasto-hydrodynamic prob-

lem,” Journal of Mechanical Engineering Science, vol. 1, no. 1, pp. 6–15, 1959.

[12] R. Wolff and A. Kubo, “The application of newton-raphson method to thermal elastohy-
drodynamic lubrication of line contacts,” Journal of tribology, vol. 116, no. 4, pp. 733–740,
1994.

[13] E. Nurgat, M. Berzins, and L. Scales, “Solving ehl problems using iterative, multigrid, and
homotopy methods,” Journal of tribology, vol. 121, no. 1, pp. 28–33, 1999.

[14] K. P. Oh and S. M. Rohde, “Numerical solution of the point contact problem using the
finite element method,” International Journal for Numerical Methods in Engineering,
vol. 11, no. 10, pp. 1507–1518, 1977.

[15] A. A. Lubrecht, W. E. Ten Napel, and R. Bosma, “Multigrid, an alternative method for cal-
culating film thickness and pressure profiles in elastohydrodynamically lubricated line
contacts,” Journal of tribology, vol. 108, no. 4, pp. 551–556, 1986.

79

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 BIBLIOGRAPHY

[16] A. A. Lubrecht, W. E. Ten Napel, and R. Bosma, “Multigrid, an alternative method of solu-
tion for two-dimensional elastohydrodynamically lubricated point contact calculations,”
Journal of tribology, vol. 109, no. 3, pp. 437–443, 1987.

[17] T. Almqvist and R. Larsson, “The navier–stokes approach for thermal ehl line contact
solutions,” Tribology International, vol. 35, no. 3, pp. 163–170, 2002.

[18] Q. J. Wang and Y.-W. Chung, eds., History of EHL Development, pp. 1684–1684. Boston,
MA: Springer US, 2013.

[19] H. Christensen and K. Tonder, “The hydrodynamic lubrication of rough bearing surfaces
of finite width,” Journal of Lubrication Technology, vol. 93, no. 3, pp. 324–329, 1971.

[20] H. Christensen, “Some aspects of the functional influence of surface roughness in lubri-
cation,” Wear, vol. 17, no. 2, pp. 149–162, 1971.

[21] N. Patir and H. Cheng, “An average flow model for determining effects of three-
dimensional roughness on partial hydrodynamic lubrication,” Journal of lubrication
Technology, vol. 100, no. 1, pp. 12–17, 1978.

[22] N. Patir and H. Cheng, “Application of average flow model to lubrication between rough
sliding surfaces,” Journal of Lubrication Technology, vol. 101, no. 2, pp. 220–229, 1979.

[23] J. Tripp, “Surface roughness effects in hydrodynamic lubrication: the flow factor
method,” Journal of lubrication technology, vol. 105, no. 3, pp. 458–463, 1983.

[24] Y. Hu and L. Zheng, “Some aspects of determining the flow factors,” Journal of tribology,
vol. 111, no. 3, pp. 525–531, 1989.

[25] L. Lunde and K. Tonder, “Pressure and shear flow in a rough hydrodynamic bearing, flow
factor calculation,” Journal of tribology, vol. 119, no. 3, pp. 549–555, 1997.

[26] D. Zhu and H. Cheng, “Effect of surface roughness on the point contact ehl,” Journal of
tribology, vol. 110, no. 1, pp. 32–37, 1988.

[27] S. R. Harp and R. F. Salant, “An average flow model of rough surface lubrication with
inter-asperity cavitation,” Journal of Tribology, vol. 123, no. 1, pp. 134–143, 2001.

[28] G. Bayada, V. Carlos, et al., “An average flow model of the reynolds roughness including a
mass-flow preserving cavitation model,” Journal of Tribology, vol. 127, no. 4, pp. 793–802,
2005.

[29] N. Letalleur, F. Plouraboueé, and M. Prat, “Average flow model of rough surface lubrica-
tion: flow factors for sinusoidal surfaces,” Journal of tribology, vol. 124, no. 3, pp. 539–546,
2002.

[30] F. Sahlin, A. Almqvist, R. Larsson, and S. Glavatskih, “Rough surface flow factors in full
film lubrication based on a homogenization technique,” Tribology International, vol. 40,
no. 7, pp. 1025–1034, 2007.

[31] P. Goglia, T. Conry, and C. Cusano, “The effects of surface irregularities on the elastohy-
drodynamic lubrication of sliding line contacts. part i—single irregularities,” Journal of
Tribology, vol. 106, no. 1, pp. 104–112, 1984.

80

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 BIBLIOGRAPHY

[32] P. Goglia, C. Cusano, and T. Conry, “The effects of surface irregularities on the elastohy-
drodynamic lubrication of sliding line contacts. part ii—wavy surfaces,” Journal of Tri-
bology, vol. 106, no. 1, pp. 113–119, 1984.

[33] L. G. Houpert and B. J. Hamrock, “Elastohydrodynamic lubrication calculations used as

a tool to study scuffing,” Proceedings of the 12th Leeds-Lyon Symposium on Tribology,
Tribology Series, pp. 146–155, 1985.

[34] C. H. Venner and W. E. Ten Napel, “Surface roughness effects in an ehl line contact,”
Journal of tribology, vol. 114, no. 3, pp. 616–622, 1992.

[35] C. Hooke, “The behaviour of low-amplitude surface roughness under line contacts,” Pro-
ceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribol-
ogy, vol. 213, no. 4, pp. 275–285, 1999.

[36] A. A. Lubrecht, “The numerical solution of the elastohydrodynamically lubricated line-

and point contact problem using multigrid techniqes,” PhD. Thesis, University of Twente,
1987.

[37] A. A. Lubrecht, W. Ten Napel, and R. Bosma, “The influence of longitudinal and trans-
verse roughness on the elastohydrodynamic lubrication of circular contacts,” Journal of
tribology, vol. 110, no. 3, pp. 421–426, 1988.

[38] C. Kweh, H. Evans, and R. Snidle, “Micro-elastohydrodynamic lubrication of an elliptical

contact with transverse and three-dimensional sinusoidal roughness,” Journal of Tribol-
ogy, vol. 111, no. 4, pp. 577–584, 1989.

[39] C. C. Kweh, M. Patching, H. Evans, and R. Snidle, “Simulation of elastohydrodynamic

contacts between rough surfaces,” Journal of tribology, vol. 114, no. 3, pp. 412–419, 1992.

[40] X. Ai and H. S. Cheng, “A transient ehl analysis for line contacts with measured surface
roughness using multigrid technique,” Journal of tribology, vol. 116, no. 3, pp. 549–556,
1994.

[41] L. Chang, M. Webster, and A. Jackson, “On the pressure rippling and roughness deforma-
tion in elastohydrodynamic lubrication of rough surfaces,” Journal of tribology, vol. 115,
no. 3, pp. 439–444, 1993.

[42] L. Chang and M. Webster, “A study of elastohydrodynamic lubrication of rough surfaces,”

Journal of Tribology, vol. 113, no. 1, pp. 110–115, 1991.

[43] J. Greenwood and K. Johnson, “The behaviour of transverse roughness in sliding elasto-
hydrodynamically lubricated contacts,” Wear, vol. 153, no. 1, pp. 107–117, 1992.

[44] J. Greenwood and G. Morales-Espejel, “The behaviour of transverse roughness in ehl

contacts,” Proceedings of the Institution of Mechanical Engineers, Part J: Journal of En-
gineering Tribology, vol. 208, no. 2, pp. 121–132, 1994.

[45] A. A. Lubrecht and C. H. Venner, “Aspects of two-sided surface waviness in an ehl line
contact,” Tribology Series, vol. 25, pp. 205–205, 1993.

[46] G. E. Morales Espejel, Elastohydrodynamic lubrication of smooth and rough surfaces.

PhD thesis, University of Cambridge, 1993.

81

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 BIBLIOGRAPHY

[47] K. F. Osborn and F. Sadeghi, “Time dependent line ehd lubrication using the multi-
grid/multilevel technique,” Journal of tribology, vol. 114, no. 1, pp. 68–74, 1992.

[48] C. H. Venner, A. A. Lubrecht, and W. E. Ten Napel, “Numerical simulation of the over-
rolling of a surface feature in an ehl line contact,” Journal of tribology, vol. 113, no. 4,
pp. 777–783, 1991.

[49] C. H. Venner, Multilevel solution of the EHL line and point contact problems. PhD thesis,
Netherlands, 1991.

[50] C. H. Venner and A. A. Lubrecht, “Numerical simulation of a transverse ridge in a circular

ehl contact under rolling/sliding,” Journal of tribology, vol. 116, no. 4, pp. 751–761, 1994.

[51] C. H. Venner and A. A. Lubrecht, “Numerical simulation of waviness in a circular ehl

contact, under rolling/sliding,” in Tribology Series, vol. 30, pp. 259–272, Elsevier, 1995.

[52] C. H. Venner, F. Couhier, A. A. Lubrecht, and J. A. Greenwood, “Amplitude reduction of

waviness in transient ehl line contacts,” Tribology Series, vol. 32, pp. 103–112, 1997.

[53] C. H. Venner and A. A. Lubrecht, “Amplitude reduction of non-isotropic harmonic pat-

terns in circular ehl contacts, under pure rolling,” in Tribology Series, vol. 36, pp. 151–162,
Elsevier, 1999.

[54] A. A. Lubrecht, D. Graille, C. H. Venner, and J. A. Greenwood, “Waviness amplitude re-

duction in ehl line contacts under rolling-sliding,” Journal of tribology, vol. 120, no. 4,
pp. 705–709, 1998.

[55] C. H. Venner and G. E. Morales-Espejel, “Amplitude reduction of small-amplitude wavi-

ness in transient elastohydrodynamically lubricated line contacts,” Proceedings of the In-
stitution of Mechanical Engineers, Part J: Journal of Engineering Tribology, vol. 213, no. 6,
pp. 487–504, 1999.

[56] A. A. Lubrecht and C. H. Venner, “Elastohydrodynamic lubrication of rough surfaces,”

Proceedings of the institution of mechanical engineers, Part J: Journal of Engineering Tri-
bology, vol. 213, no. 5, pp. 397–404, 1999.

[57] B. Jacod, P. M. Lugt, M. L. Dumont, J. H. Tripp, and C. H. Venner, “Amplitude reduction of

waviness in elastohydrodynamic lubrication using an eyring fluid model,” Proceedings of
the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, vol. 214,
no. 4, pp. 343–350, 2000.

[58] C. H. Venner and A. A. Lubrecht, “An engineering tool for the quantitative prediction of
general roughness deformation in ehl contacts based on harmonic waviness attenua-
tion,” Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineer-
ing Tribology, vol. 219, no. 5, pp. 303–312, 2005.

[59] A. D. Chapkov, C. H. Venner, and A. A. Lubrecht, “Roughness amplitude reduction under

non-newtonian ehd lubrication conditions,” Journal of tribology, vol. 128, no. 4, pp. 753–
760, 2006.

[60] J. Wang, C. H. Venner, and A. A. Lubrecht, “Amplitude reduction in ehl line contacts under
rolling sliding conditions,” Tribology international, vol. 44, no. 12, pp. 1997–2001, 2011.

[61] P. Šperka, I. Křupka, and M. Hartl, “Experimental study of real roughness attenuation in
concentrated contacts,” Tribology International, vol. 43, no. 10, pp. 1893–1901, 2010.

82

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 BIBLIOGRAPHY

[62] M. Holmes, H. P. Evans, T. G. Hughes, and R. W. Snidle, “Transient elastohydrodynamic

point contact analysis using a new coupled differential deflection method part 1: theory
and validation,” Proceedings of the Institution of Mechanical Engineers, Part J: Journal of
Engineering Tribology, vol. 217, no. 4, pp. 289–304, 2003.

[63] H. P. Evans, R. W. Snidle, and K. J. Sharif, “Deterministic mixed lubrication modelling

using roughness measurements in gear applications,” Tribology International, vol. 42,
no. 10, pp. 1406–1417, 2009.

[64] Y. Zhang, H. Liu, C. Zhu, C. Song, and Z. Li, “Influence of lubrication starvation and sur-
face waviness on the oil film stiffness of elastohydrodynamic lubrication line contact,”
Journal of Vibration and Control, vol. 24, no. 5, pp. 924–936, 2018.

[65] Z. Gu, C. Zhu, H. Liu, and X. Du, “A comparative study of tribological performance of
helical gear pair with various types of tooth surface finishing,” Industrial Lubrication and
Tribology, 2018.

[66] R. E. Alcouffe, A. Brandt, J. E. Dendy, Jr, and J. W. Painter, “The multi-grid method for the
diffusion equation with strongly discontinuous coefficients,” SIAM Journal on Scientific
and Statistical Computing, vol. 2, no. 4, pp. 430–454, 1981.

[67] A. Crook, “The lubrication of rollers iv. measurements of friction and effective viscosity,”
Phil. Trans. R. Soc. Lond. A, vol. 255, no. 1056, pp. 281–312, 1963.

[68] K. L. Johnson and R. Cameron, “Fourth paper: shear behaviour of elastohydrodynamic

oil films at high rolling contact pressures,” Proceedings of the Institution of Mechanical
Engineers, vol. 182, no. 1, pp. 307–330, 1967.

[69] K. L. Johnson and A. D. Roberts, “Observations of viscoelastic behaviour of an elastohy-

drodynamic lubricant film,” Proc. R. Soc. Lond. A, vol. 337, no. 1609, pp. 217–242, 1974.

[70] C. R. Evans and K. L. Johnson, “Regimes of traction in elastohydrodynamic lubrication,”

Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical En-
gineering Science, vol. 200, no. 5, pp. 313–324, 1986.

[71] Z. Zhang, X. Qiu, and Y. Hong, “Ehl analysis of rib-roller end contact in tapered roller
bearings,” Tribology transactions, vol. 31, no. 4, pp. 461–467, 1988.

[72] H. Yu and J. B. Medley, “Influence of lubricant additives on friction in a disc machine,”

Tribology Series, vol. 32, pp. 475–486, 1997.

[73] B. Jacod, C. H. Venner, and P. M. Lugt, “A generalized traction curve for ehl contacts,”
Journal of tribology, vol. 123, no. 2, pp. 248–253, 2001.

[74] J. De Vicente, J. R. Stokes, and H. A. Spikes, “The frictional properties of newtonian fluids
in rolling–sliding soft-ehl contact,” Tribology Letters, vol. 20, no. 3-4, pp. 273–286, 2005.

[75] H. Liu, C. Zhu, Z. Sun, and C. Song, “Starved lubrication of a spur gear pair,” Tribology
International, vol. 94, pp. 52–60, 2016.

[76] M. Björling, K. Berglund, A. Spencer, and R. Larsson, “The effect of ageing on elastohy-
drodynamic friction in heavy-duty diesel engine oils,” Proceedings of the Institution of
Mechanical Engineers, Part J: Journal of Engineering Tribology, vol. 231, no. 6, pp. 708–
715, 2017.

83

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 BIBLIOGRAPHY

[77] J. Zhang, A. Tan, and H. Spikes, “Effect of base oil structure on elastohydrodynamic fric-
tion,” Tribology Letters, vol. 65, no. 1, p. 13, 2017.

[78] W. C. Emmens, “The influence of surface roughness on friction,” in International Confer-

ence on Controlling Sheet Metal Forming Processes, Michigan, pp. 63–70, 1988.

[79] S. Wu and H. S. Cheng, “A friction model of partial-ehl contacts and its application to
power loss in spur gears,” Tribology Transactions, vol. 34, no. 3, pp. 398–407, 1991.

[80] B. H. Lee, Y. T. Keum, and R. H. Wagoner, “Modeling of the friction caused by lubrication
and surface roughness in sheet metal forming,” Journal of materials processing technol-
ogy, vol. 130, pp. 60–63, 2002.

[81] R. H. Thurston, Friction and Lubrication: Determinations of the Laws and Coëfficients of
Friction by New Methods and with New Apparatus. Railroad gazette, 1879.

[82] D. Dowson, History of tribology. Addison-Wesley Longman Limited, 1979.

[83] R. Stribeck, “Die wesentlichen eigenschaften der gleit-und rollenlager,” Zeitschrift des
Vereines Deutscher Ingenieure, vol. 46, pp. 1341–1348, 1902.

[84] R. Stribeck and M. Schröter, Die wesentlichen Eigenschaften der Gleit-und Rollenlager:
Untersuchung einer Tandem-Verbundmaschine von 1000 PS. Springer, 1903.

[85] L. Gümbel, “Das problem der lagerreibung,” Mbl. Berlin. Bez. Ver. dtsch. Ing, vol. 5, 1914.

[86] M. D. Hersey, “The laws of lubrication of horizontal journal bearings,” Journal of the
Washington Academy of Sciences, vol. 4, no. 19, pp. 542–552, 1914.

[87] R. E. Wilson and D. P. Barnard, “The mechanism of lubrication,” in SAE Technical Paper,
SAE International, 1922.

[88] S. A. McKee, “The effect of running-in on journal bearing performance,” Mech. Eng,
vol. 49, pp. 1335–1340, 1927.

[89] G. Vogelpohl, “Die stribeck-kurve als kennzeichen des allgemeinen reibungsverhaltens

geschmierter gleitflächen,” Z. VDI, vol. 96, no. 9, pp. 261–68, 1954.

[90] K. L. Johnson, “Regimes of elastohydrodynamic lubrication,” Journal of Mechanical En-

gineering Science, vol. 12, no. 1, pp. 9–16, 1970.

[91] H. A. Spikes, “Mixed lubrication—an overview,” Lubrication Science, vol. 9, no. 3, pp. 221–
253, 1997.

[92] B. A. Shotter, “Experiments with a disc machine to determine the possible influence of
surface finish on gear tooth performance,” in Proc. Int. Conf. Gearing, vol. 120, Instn
Mech. Engrs London, 1958.

[93] T. E. Tallian, Y. P. Chiu, D. F. Huttenlocher, J. A. Kamenshine, L. B. Sibley, and N. E. Sin-

dlinger, “Lubricant films in rolling contact of rough surfaces,” ASLE TRANSACTIONS,
vol. 7, no. 2, pp. 109–126, 1964.

[94] T. E. Tallian, J. I. McCool, and L. B. Sibley, “Paper 14: Partial elastohydrodynamic lubri-
cation in rolling contact,” in Proceedings of the Institution of Mechanical Engineers, Con-
ference Proceedings, vol. 180, pp. 169–186, SAGE Publications Sage UK: London, England,
1965.

84

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 BIBLIOGRAPHY

[95] S. Y. Poon and D. J. Haines, “Third paper: frictional behaviour of lubricated rolling-
contact elements,” Proceedings of the Institution of Mechanical Engineers, vol. 181, no. 1,
pp. 363–389, 1966.

[96] S. Bair and W. O. Winer, “Regimes of traction in concentrated contact lubrication,” Jour-
nal of Lubrication Technology, vol. 104, no. 3, pp. 382–386, 1982.

[97] G. Stachowiak and A. W. Batchelor, Engineering tribology. Butterworth-Heinemann,

2013.

[98] P. M. Cann, E. Ioannides, B. Jacobson, and A. Lubrecht, “The lambda ratio—a critical re-
examination,” Wear, vol. 175, no. 1-2, pp. 177–188, 1994.

[99] D. J. Schipper and A. De Gee, “On the transitions in the lubrication of concentrated con-
tacts,” Journal of tribology, vol. 117, no. 2, pp. 250–254, 1995.

[100] E. Gelinck and D. J. Schipper, “Calculation of stribeck curves for line contacts,” Tribology
International, vol. 33, no. 3-4, pp. 175–181, 2000.

[101] K. L. Johnson, J. A. Greenwood, and S. Y. Poon, “A simple theory of asperity contact in

elastohydro-dynamic lubrication,” Wear, vol. 19, no. 1, pp. 91–108, 1972.

[102] X. Lu, M. M. Khonsari, and E. R. Gelinck, “The stribeck curve: experimental results and
theoretical prediction,” Journal of tribology, vol. 128, no. 4, pp. 789–794, 2006.

[103] W. Wang, S. Wang, F. Shi, Y. Wang, H. Chen, H. Wang, and Y. Hu, “Simulations and
measurements of sliding friction between rough surfaces in point contacts: From ehl
to boundary lubrication,” Journal of Tribology, vol. 129, no. 3, pp. 495–501, 2007.

[104] M. Kalin, I. Velkavrh, and J. Vižintin, “The stribeck curve and lubrication design for non-
fully wetted surfaces,” Wear, vol. 267, no. 5-8, pp. 1232–1240, 2009.

[105] M. Kalin and I. Velkavrh, “Non-conventional inverse-stribeck-curve behaviour and other

characteristics of dlc coatings in all lubrication regimes,” Wear, vol. 297, no. 1-2, pp. 911–
918, 2013.

[106] X. Zhang, Z. Li, and J. Wang, “Friction prediction of rolling-sliding contact in mixed ehl,”
Measurement, vol. 100, pp. 262–269, 2017.

[107] J. Bonaventure, J. Cayer-Barrioz, and D. Mazuyer, “Transition between mixed lubrication

and elastohydrodynamic lubrication with randomly rough surfaces,” Tribology Letters,
vol. 64, no. 3, p. 44, 2016.

[108] D. Dowson, J. F. Dunn, and C. M. Taylor, “The piezo-viscous fluid, rigid solid regime of
lubrication,” Proceedings of the Institution of Mechanical Engineers, Part C: Journal of
Mechanical Engineering Science, vol. 197, no. 1, pp. 43–52, 1983.

[109] T. D. Jacobs, T. Junge, and L. Pastewka, “Quantitative characterization of surface topogra-

phy using spectral analysis,” Surface Topography: Metrology and Properties, vol. 5, no. 1,
p. 013001, 2017.

[110] C. H. Venner, W. E. t. ten Napel, and R. Bosma, “Advanced multilevel solution of the ehl
line contact problem,” Journal of Tribology, vol. 112, no. 3, pp. 426–431, 1990.

85

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 BIBLIOGRAPHY

[111] C. H. Venner, “Higher-order multilevel solvers for the ehl line and point contact prob-
lem,” Journal of tribology, vol. 116, no. 4, pp. 741–750, 1994.

[112] C. H. Venner and A. A. Lubrecht, “Multigrid techniques: a fast and efficient method
for the numerical simulation of elastohydrodynamically lubricated point contact prob-
lems,” Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineer-
ing Tribology, vol. 214, no. 1, pp. 43–62, 2000.

[113] A. A. Lubrecht and E. Ioannides, “A fast solution of the dry contact problem and the asso-
ciated sub-surface stress field, using multilevel techniques,” Journal of tribology, vol. 113,
no. 1, pp. 128–133, 1991.

[114] D. Dowson and G. R. Higginson, “Elasto-hydrodynamic lubrication: the fundamentals of

roller and gear lubrication. vol. 23,” 1966.

[115] C. Barus, “Isothermals, isopiestics and isometrics relative to viscosity,” American Journal
of Science (1880-1910), vol. 45, no. 266, p. 87, 1893.

[116] C. J. A. Roelands, “Correlational aspects of the viscosity-temperature-pressure relation-

ship of lubricating oils,” Doctoral thesis, Technische Hogeschool te Delt, 1966.

[117] K. L. Johnson, Contact Mechanics. Cambridge university press, 1987.

[118] H. Moes, “Discussion on a contribution by K. Jakobsen and H. Christensen,” Proc. Inst.

Mech. Eng., Part J: J. Eng. Tribol, vol. 183, no. 3, pp. 205–206, 1969.

[119] H. Moes and R. Bosma, “Design charts for optimum bearing configurations: 1—the full
journal bearing,” Journal of Lubrication Technology, vol. 93, no. 2, pp. 302–305, 1971.

[120] B. J. Hamrock and D. Dowson, “Isothermal elastohydrodynamic lubrication of point con-

tacts: Part 1—theoretical formulation,” Journal of Lubrication Technology, vol. 98, no. 2,
pp. 223–228, 1976.

[121] K. F. Riley, M. P. Hobson, and S. J. Bence, “Mathematical methods for physics and engi-
neering,” 1999.

[122] W. McGuire, R. H. Gallagher, and R. D. Ziemian, Matrix structural analysis. 2000.

[123] M.-P. Noutary, N. Biboulet, and A. A. Lubrecht, “A robust piston ring lubrication solver:
Influence of liner groove shape, depth and density,” Tribology International, vol. 100,
pp. 35–40, 2016.

[124] M.-P. Noutary, A robust Reynolds solver for textured surfaces in the piston ring cylinder
liner contact. Theses, Université de Lyon, Nov. 2017.

[125] A. Ponjavic, Local rheology of lubricants in the elastohydrodynamic regime. PhD thesis,
Imperial College London, 2014.

[126] A. V. Olver and H. A. Spikes, “Prediction of traction in elastohydrodynamic lubrication,”

Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tri-
bology, vol. 212, no. 5, pp. 321–332, 1998.

[127] S. Bair and W. O. Winer, “The high pressure high shear stress rheology of liquid lubri-
cants,” Journal of tribology, vol. 114, no. 1, pp. 1–9, 1992.

86

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 BIBLIOGRAPHY

[128] S. Bair, “Shear thinning correction for rolling/sliding elastohydrodynamic film thick-
ness,” Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineer-
ing Tribology, vol. 219, no. 1, pp. 69–74, 2005.

[129] P. Kumar, M. M. Khonsari, and S. Bair, “Full ehl simulations using the actual ree-eyring
model for shear-thinning lubricants,” Journal of Tribology, vol. 131, no. 1, p. 011802, 2009.

[130] S. Bair, “The rheological assumptions of classical ehl: What went wrong?,” Tribology In-
ternational, vol. 131, pp. 45–50, 2019.

[131] K. L. Johnson and J. L. Tevaarwerk, “Shear behaviour of elastohydrodynamic oil films,”

Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences,
vol. 356, no. 1685, pp. 215–236, 1977.

[132] S. Huge and J. Zhang, “History, origins and prediction of elastohydrodynamic friction,”
Tribology Letters, vol. 56, no. 1, pp. 1–25, 2014.

[133] S. Bair and W. O. Winer, “A rheological model for elastohydrodynamic contacts based on
primary laboratory data,” Journal of Lubrication Technology, vol. 101, no. 3, pp. 258–264,
1979.

[134] K. Yasuda, R. C. Armstrong, and R. E. Cohen, “Shear flow properties of concentrated so-
lutions of linear and star branched polystyrenes,” Rheologica Acta, vol. 20, pp. 163–178,
Mar 1981.

[135] R.-T. Lee and B. Hamrock, “A circular non-newtonian fluid model: Part i—used in elasto-
hydrodynamic lubrication,” Journal of Tribology, vol. 112, no. 3, pp. 486–495, 1990.

[136] S. Bair, “The actual eyring models for thixotropy and shear-thinning: Experimental vali-
dation and application to ehd,” in ASME/STLE 2004 International Joint Tribology Confer-
ence, pp. 439–445, American Society of Mechanical Engineers, 2004.

[137] P. Kumar and M. M. Khonsari, “On the role of lubricant rheology and piezo-viscous prop-
erties in line and point contact ehl,” Tribology International, vol. 42, no. 11-12, pp. 1522–
1530, 2009.

[138] H. Rahnejat, Tribology and dynamics of engine and powertrain: fundamentals, applica-
tions and future trends. Elsevier, 2010.

[139] S. Wen and P. Huang, Principles of tribology. John Wiley & Sons, 2012.

[140] J. Bonaventure, Influence of random surface roughness on friction in elastohydrodynamic,

mixed and boundary lubrication. PhD thesis, Lyon, 2017.

[141] C. R. Evans, Measurement and mapping of the rheological properties of elastohydrody-

namic lubricants. PhD thesis, University of Cambridge, 1984.

[142] B. Bhushan, Modern tribology handbook, two volume set. CRC press, 2000.

[143] J. C. Stover, Optical scattering: measurement and analysis, vol. 2. SPIE optical engineering
press Bellingham, 1995.

[144] A.-L. Barabási and H. E. Stanley, Fractal concepts in surface growth. Cambridge university
press, 1995.

87

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 BIBLIOGRAPHY

[145] B. Persson, O. Albohr, U. Tartaglino, A. Volokitin, and E. Tosatti, “On the nature of surface
roughness with application to contact mechanics, sealing, rubber friction and adhesion,”
Journal of physics: Condensed matter, vol. 17, no. 1, p. R1, 2004.

[146] P. R. Nayak, “Random process model of rough surfaces,” Journal of Lubrication Technol-
ogy, vol. 93, no. 3, pp. 398–407, 1971.

[147] D. J. Schipper, “Transitions in the lubrication of concentrated contacts, 1988,” University

of Twente: Enschede, The Netherlands, p. 194, 1988.

88

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 FOLIO ADMINISTRATIF

THESE DE L’UNIVERSITE DE LYON OPEREE AU SEIN DE L’INSA LYON

NOM : ZHANG DATE de SOUTENANCE : 05/09/2019

Prénoms : Yuanyuan

TITRE : La prédiction du frottement des surfaces rugueuses dans un contact avec lubrification
élastohydrodynamique

NATURE : Doctorat Numéro d'ordre : 2019LYSEI063

Ecole doctorale : Mécanique, Energétique, Génie Civil, Acoustique (MEGA)

Spécialité : Génie Mécanique

RESUME :

Le frottement à l’interface des surfaces influence les performances des éléments mécaniques. Le frottement a été
étudié expérimentalement dans la plupart des études. Dans ce travail, le frottement est prédit à l'aide d'une
simulation numérique dans des conditions de contact rugueux avec une lubrification élastohydrodynamique (EHL).
La technique classique Multigrille fonctionne bien pour limiter le temps de calcul et les besoins en mémoire.
Cependant, le choix de la grille grossière a une influence importante sur la robustesse du code et son efficacité
pour résoudre le problème brut. Dans la première partie de ce travail, une méthode de construction de grille
grossière proposée par Alcouffe et al. est implémenté dans le code EHL Multigrille indépendamment du temps.
Ensuite ce solveur modifié est étendu aux cas transitoires pour résoudre le problème de contact avec rugosité.
La courbe de frottement est généralement représentée en fonction du « λ ratio », le rapport entre l'épaisseur du
film d'huile et la valeur moyenne quadratique de la rugosité de la surface. Cependant, ce paramètre est moins
approprié pour tracer les variations de frottement dans des conditions de haute pression (régime élasto piézo-
visqueux). Dans la deuxième partie de ce travail, le coefficient de frottement est calculé à l'aide du code EHL
modifié pour de nombreuses conditions de fonctionnement ainsi que pour les paramètres d'ondulation de surface.
Les résultats de la simulation montrent qu'il n'y a pas de courbe de frottement unique lorsque l'ancien paramètre «
λ ratio » est utilisé. En se basant sur la théorie de la réduction d'amplitude, un nouveau paramètre de
dimensionnement qui dépend des conditions de fonctionnement et des paramètres d'ondulation est trouvé, ce qui
peut donner une courbe de frottement unique pour les situations de haute pression.
Le nouveau paramètre d’échelle et l’augmentation du frottement prédite par la méthode PSD montrent une bonne
précision technique pour une utilisation pratique.

MOTS-CLÉS : Lubrification élastohydrodynamique, Simulation numérique, Régime élasto piézo-visqueux, Théorie

de réduction d’amplitude, Variation de frottement

Laboratoire (s) de recherche : Laboratoire de Mécanique des Contacts et des Structures

UMR CNRS 5259 - INSA de Lyon
Bâtiment Sophie Germain
Avenue Jean Capelle
69621 Villeurbanne Cedex FRANCE

Directeur de thèse: LUBRECHT Antonius. A

Président de jury : CAYER-BARRIOZ Juliette

Composition du jury : EVANS Pwt KŘUPKA Ivan

VENNER Cornelis. H CAYER-BARRIOZ Juliette
BIBOULET Nans LUBRECHT Antonius. A

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés
 BIBLIOGRAPHY

90

Cette thèse est accessible à l'adresse : [Link]

© [Y. Zhang], [2019], INSA Lyon, tous droits réservés