Rough surface elastohydrodynamic lubrication and ... - DiVA portal

:

L ICENTIATE T H E S I S

Rough Surface Elastohydrodynamic Lubrication and Contact Mechanics

Andreas Almqvist

Luleå University of Technology Department of Applied Physics and Mechanical Engineering, Division of Machine Elements :|: -|: - -- ⁄ -- 

2004:35

ROUGH S URFACE E LASTOHYDRODYNAMIC L UBRICATION AND C ONTACT M ECHANICS

A NDREAS A LMQVIST

Luleå University of Technology Department of Applied Physics and Mechanical Engineering, Division of Machine Elements

2004 : 35 | ISSN : 1402 − 1757 |ISRN:LTU-LIC--04/35--SE

Cover figure:

A modeled surface topography pressed against a rigid plane, assuming linear elastic surface material. The theory describing the contact mechanics tool used to produce this result is given in Chapter 3

Title page figure:

Elementary surface features passing each other inside the EHD lubricated conjunction, see Fig. 6.5 for details.

ROUGH S URFACE E LASTOHYDRODYNAMIC L UBRICATION AND C ONTACT M ECHANICS c Andreas Almqvist (2004). This document is freely available at Copyright http://epubl.ltu.se/1402-1757/2004/35 or by contacting Andreas Almqvist, [email protected] The document may be freely distributed in its original form including the current author’s name. None of the content may be changed or excluded without permissions from the author.

ISSN: ISRN:

1402-1757 LTU-LIC--04/35--SE

This document was typeset in LATEX 2ε .

Abstract In the field of tribology, there are numerous theoretical models that may be described mathematically in the form of integro-differential systems of equations. Some of these systems of equations are sufficiently well posed to allow for numerical solutions to be carried out resulting in accurate predictions. This work has focused on the contact between rough surfaces with or without a separating lubricant film. The objective was to investigate how surface topography influences contact conditions. For this purpose two different numerical methods were developed and used. For the lubricated contact between rough surfaces the Reynolds equation were used as a basis. This equation is derived under the assumptions of thin fluid film and creeping flow. In highly loaded, lubricated, non-conformal contacts of surfaces after running-in, the load concentration no longer results in plastic deformations, however large elastic deformations will be apparent. It is the interaction between the hydrodynamic action of the lubricant and the elastic deformations of the surfaces that, in certain applications, enable the lubricant film to fully separate the surfaces. This is commonly referred to as full film elastohydrodynamic (EHD) lubrication. Typical machine elements that operate in the full film EHD lubrication (FL) regime include rolling element bearings, cams and gears. Unfortunately, a cost effective way of machining engineering surfaces seldom results in a surface topography that influence contact conditions in the same way as a surface after running-in. Such topographies may prevent the lubricant from fully separating the surfaces because of deteriorated hydrodynamic action. In this case the applied load is carried in part by the lubricant and in part by surface asperities and/or surface active lubricant additives. This could also be the case in lubricant starved contacts, which is a common situation in not only grease lubricated contacts but also in many liquid lubricated contacts, such as high speed operating rolling element bearings. The load sharing between the highly compressed lubricant and the surface and/or surface active lubricant additives is the reason why this lubrication regime is most commonly referred to as mixed EHD lubrication (ML). Machine elements that while running operate in the FL regime may experience a transition into the ML regime at stops or due to altered operating conditions. It is not possible to simulate direct contact between the surfaces using a numerical method based on Reynolds equation. A parameter study, of elementary surface features passing each other inside the EHD lubricated conjunction, was performed. The results obtained, even though no direct contact could be simulated, does indicate that a transition from the FL to the ML regime would occur for certain combinations of the varied parameters. At start-ups, the contact in a rolling element bearing could be both starved and drained from lubricant. In this case the hydrodynamic action becomes negligible in terms of load carrying capacity. The load is carried exclusively by surface asperities and/or surface aci

tive lubricant additives. This regime is referred to as boundary lubrication (BL). Operation conditions could also make both FL and ML impossible to achieve, for example, in the case in a low rpm operating rolling element bearing. The BL regime is in this work modeled as the unlubricated frictionless contact between rough surfaces, i.e., a dry contact approach. A variational principle was used in which the real area of contact and contact pressure distribution are those which minimize the total complementary energy. A linear elastic-perfectly plastic deformation model in which energy dissipation due to plastic deformation is accounted for was used. The dry contact method was applied to the contact between four different profiles and a plane. The variation in the real area of contact, the plasticity index and some surface roughness parameters due to applied load were investigated. The surface roughness parameters of the profiles differed significantly.

Contents I

The Thesis

1 Introduction 1.1 Elastohydrodynamic lubrication 1.2 Lubrication regimes . . . . . . . 1.3 Surface topography . . . . . . . 1.4 Objectives . . . . . . . . . . . . 1.5 Outline of this thesis . . . . . .

1 . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

3 3 4 5 5 7

2 Full Film EHD Lubrication 2.1 Deterministic roughness models . . . . . . . . . 2.2 Governing equations . . . . . . . . . . . . . . . 2.3 The Block-Jacobi method . . . . . . . . . . . . . 2.3.1 Ordinary Jacobi . . . . . . . . . . . . . . 2.3.2 Coupled systems of equations . . . . . . 2.4 A brief overview of the multilevel technique . . . 2.4.1 Grid levels . . . . . . . . . . . . . . . . 2.4.2 Intergrid transfer operators . . . . . . . . 2.4.3 Two level solver for the Poisson equation 2.4.4 Two-dimensional functional operators . . 2.5 Dimensionless formulation . . . . . . . . . . . . 2.6 Discrete formulation . . . . . . . . . . . . . . . 2.7 Solution method . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

9 9 11 15 15 16 17 17 18 19 20 21 23 24

3 Dry Elasto-Plastic Contact 3.1 Statistical roughness models . . 3.2 Deterministic roughness models 3.3 Numerical solution techniques . 3.4 Governing equations . . . . . . 3.5 Spectral analysis . . . . . . . . 3.6 Dimensionless formulation . . . 3.7 Discrete formulation . . . . . . 3.8 Solution method . . . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

27 27 28 28 29 30 31 32 33

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

iii

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . .

. . . . . . . .

. . . . . . . .

4 Surface Characterization 4.1 Theoretical model . . . 4.2 Measured topographies 4.3 Parameter study . . . . 4.4 Conclusions . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

35 35 35 36 39

5 Reynolds vs. CFD 5.1 The CFD approach . . . . . . . . . 5.2 CFD - Governing equations . . . . . 5.3 The model problem . . . . . . . . . 5.3.1 Interpolation of solution data 5.3.2 Error estimation . . . . . . 5.4 The results of the comparison . . . . 5.5 Discussion and concluding remarks

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

41 41 42 42 43 44 46 46

6 Simulations of Rough FL 6.1 The different overtaking situations . . . . . . . . . . . . . . . . . . . . . 6.2 The Dent-Ridge overtaking . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51 51 54 58

7 Simulations of the dry contact 7.1 Varying the applied load . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61 62 64

8 Concluding Remarks

65

9 Future Work

67

II Appended Papers

69

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

A A.1 Introduction . . . . . . . . . . . . . . . . . . . . . . A.2 Governing equations . . . . . . . . . . . . . . . . . A.2.1 Boundary conditions and cavitation treatment A.3 Numerics . . . . . . . . . . . . . . . . . . . . . . . A.3.1 The numerics for the CFD approach . . . . . A.3.2 The numerics for the Reynolds approach . . A.3.3 Error estimation . . . . . . . . . . . . . . . A.3.4 Interpolation of solution data . . . . . . . . . A.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . A.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . A.6 Conclusions . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

71 74 74 76 76 77 77 78 79 80 81 84

B.1 Introduction . . . B.2 Theory . . . . . . B.2.1 Equations B.2.2 Numerics

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

87 90 91 91 92

B . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

B.2.3 Error estimation . . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . B.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C.1 C.2 C.3 C.4 C.5

Introduction . . . . . . . Theory . . . . . . . . . . Surface characterization . Results . . . . . . . . . . Conclusions . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

93 94 98 103 106 107 108 109 112

Preface This licentiate thesis comprises the results from numerical simulations of both lubricated and unlubricated contacts, specifically on the influence of surface topography on contact conditions. The work has been carried out at the Department of Applied Physics and Mechanical Engineering, the Division of Machine Elements, Luleå University of Technology and the thesis is based on the three papers, A, B and C found in part II. Graduate studies are time consuming which means little time for other activities. This makes it hard work being the father of three wonderful children. My greatest gratitude is therefore given to my wife Ulrika who by being a such fantastic mother and taking care of our children has made these studies possible. For all those inspiring moments with my children I am also very grateful. My parents, grand parents and other close relatives and friends also deserve my gratitude for all their support. I would like to thank my supervisor Dr. Roland Larsson for introducing me to the subject of tribology and especially the subject of EHL and for the inspiring and encouraging meetings we have had and hopefully will continue to have. I would also like to thank all my colleagues at the Division of Machine Elements for their support, their co-operation and for introducing new concepts in tribology to me. I must also say that I am grateful to be a part of this particular division because of the friendly and relaxed atmosphere. Special thanks goes to my friend and former colleague Dr. Torbjörn Almqvist for his co-operation, enthusiasm and all the stimulating discussions we have had. Another special thanks goes to Mr. Fredrik Sahlin for his co-operation, and his support and help in scripting especially in Perl and LATEX 2ε . From the Department of Mathematics I would like to thank my associate supervisor Dr. Inge Söderqvist for his support in the field of Scientific Computing. I would also like to thank Mr. Reynold Näslund for stimulating and encouraging discussions and his help with mathematics in general. Finally, I would like to thank my sponsors; Fortum, Indexator, SKF Statoil, Volvo Car Corp.and the national research programme HiMeC and the national graduate school in scientific computning, NGSSC, both financed by the Swedish Foundation for Strategic Research (SSF).

vii

Nomenclature Pa −1

α

Pressure-viscosity coefficient

A

Functional space

µ¯

Dimensionless dynamic viscosity

ρ¯

Dimensionless density

∆T

Dimensionless step size in time

s

∆t

Step size in time

s

∆X

Dimensionless step size in x

m

∆x

Step size in spatial coordinate x

m

∆z

Step size in spatial coordinate z

m

γ˙

Shear rate

E

Error

ε

Reynolds equation coefficient

s−1

εspatial Discretization error in space εtime

Discretization error in time

εx

Discretization error in space

εy

Discretization error in space

Total potential complementary energy

I

Intergrid transfer operator

Λ

Wavelength

L

Functional operator

m

m s−1

u (x, y, z) Velocity field x (x, y, z) Spatial coordinates

m ix

µ

Dynamic viscosity

Pa s

µ0

Dynamic viscosity at ambient pressure

Pa s

ν

Poissont’s ratio

Ω

Integration domain 2D / 3D contact

ωx

Spatial frequency

1/m

ωy

Spatial frequency

1/m

O (n)

Mathematical order of the number n

Φ

Solution variable

ψi

Topography of surface i

ψavg

Mean surface height

ρ

Density

kg m−3

ρ0

Density at ambient pressure

kg m −3

Σ

Dimensionless stress function

σ

Total stress tensor (A)

Pa

τ

Lubricant shear stress

Pa

τ0

Eyring stress

Pa

τ1

Lubricant shear stress at surface 1

Pa

τm

Lubricant midplane shear stress

Pa

τxz

Lubricant shear stress

Pa

f

Fourier transformation of f

1 ∇

Dimensionless reduced wavelength

ξ

Film thinning effect measure

A

Dent / Ridge Amplitude

ai

Matrix coefficient

aξ

Contact semi-width in the rolling direction

Ai j

Matrices

b

Hertzian half-width b =

bi

Matrix coefficient

bξ

Contact semi-width in the transversial direction

m / m2

m

m

(8wRx ) / (πE )

m

m

m

C1

Constant = 5.9 10 8

C2

Constant = 1.34

ci

Matrix coefficient

d

Elastic deformation

m

de

Elastic deformation

m

di

Matrix coefficient

dp

Plastic deformation

m

E

Modulus of elasticity

Pa

E

Effective mod. of elast. 2/E = (1 − ν21 )/E1 + (1 − ν22)/E2

Pa

F

Discrete right hand side (Reynolds equation)

fi

Matrix coefficient

G

Discrete right hand side (Film th. equation)

G

Shear modulus of elasticity of the lubricant

g

Determinant of the metric tensor

gi

Matrix coefficient

gs

The gap between the undeformed surfaces

H

Dimensionless film thickness

h

Film thickness

m

h0

Integration constant

m

Hs

Hardness of the softer material

Pa

H00

Dimensionless integration constant

hc

Central film thickness

m

hmin

Minimum film thickness

m

K

Elastic deformation integral kernel

L

Moes dimensionless speed parameter

li

Amplitude variation of dent/ridge

Lx

Length of the Fourier window in the x-dir.

m

Ly

Length of the Fourier window in the y-dir.

m

L1

Discrete operator (Reynolds equation)

Pa

m

L2

Discrete operator (Film thickness equation)

M

Mass flux per unit width

M

Moes dimensionless load parameter

mi

Wavelength variation of dent/ridge

P

Dimensionless pressure

p

Pressure

Pa

P0 ph

Constant in the viscosity expression Hertzian pressure p h = (wE ) / (2πRx )

Pa

pmax

Maximum pressure

Pa

Ra

Average roughness

µm

Rd

Radius of the lower surface (A)

m

Ri

Radius of the surface i

m

Rk

Kurtosis

Rq

Root mean square (RMS) roughness

Ru

Radius of the upper surface (A)

m

Rx

Reduced radius of curvature in x-dir. 1/R x = 1/R11 + 1/R22

m

Rz

Average maximum height

Rsk

Skewness

S

Non-Newtionian slip factor

s

Slide-to-roll ratio 2(u 1 − u2)/(u1 + u2 )

T

Dimensionless time

t

Time

u

Solution to the one-dimensional Poisson equation

u1

Surface velocity

m s −1

u2

Surface velocity

m s −1

us

Sum of velocities u s = u1 + u2

m s−1

W

Applied load, 2D / 3D contact

Pa m −1 / Pa

w

Applied load

X

Dimensionless spatial coordinate

kg m −1 s−1

µm

µm

s

Pa m −1 m

x

Spatial coordinate in the rolling direction

m

xc

Centre of dent / ridge x c (t) = xs − u t

m

xs

Initial placement of dent / ridge

m

z

Spatial coordinate across the film

m

z


m

z1


m

zvisc

Pressure-viscosity index

Part I

The Thesis

1

Chapter 1

Introduction Machines consists of machine elements and their safe and efficient operation relies on carefully designed interfaces between the machine elements. The functional design of interfaces covers geometry, materials, lubrication and surface topography. Incorrect choice of these design parameters may lead to both lowered efficiency and shortened service life. Poor geometry could lead to large stress concentrations which in turn may lead to rapid fatigue. Large stress concentrations also implicitly imply a temperature rise because of energy dissipation due to plastic deformations. The choice of materials is also of great importance, for example, electrolytic corrosion may drastically reduce service life. Contact fatigue contact to low ductility would not only lower the service life but could lead to third body abrasion due to spalling which in turn could end up lowering the service life of other components. A lubricant serves several crucial objectives. When its main objective is to lower friction, the action of additives is of concern. If the interface is subjected to excessive wear, the lubricant’s ability to from a separating film is even more crucial. In this case the bulk properties of the lubricant have to be carefully chosen. All real surfaces are rough at some scale and the surface topography certainly influences the contact between machine elements. The design parameters mentioned are also mutually dependent. That is, they affect the way the others influence the operation of the system. For example, a change in geometry could require to another choice of materials. The new choice of material may change the objectives of the lubricant and force the operation into another lubrication regime. All these four design parameters are of great importance, however in this thesis work it is surface topography which is the main focus.

1.1 Elastohydrodynamic lubrication Elastohydrodynamic (EHD) lubrication is the type of hydrodynamic lubrication in which the elastic deformations of the contacting surfaces cannot be neglected. This is often the case in non-conformal (concentrated) contacts. For example, the contact between the roller and the raceway in a typical roller bearing as shown in Fig. 1.1 will operate in the elastohydrodynamic lubrication regime. The actual contact zone for a rolling bearing is, in general, elliptic in shape. Depending on the design parameters previously mentioned and the actual running conditions, the shape of the ellipse will change as shown in Fig. 1.2. 3

CHAPTER 1. INTRODUCTION

4

Figure 1.1: A typical rolling element bearing In any case, the contact region is small and the concentrated load implies high pressure which will lead to large elastic deformations and could also lead to plastic deformation. For a bearing in operation, it is the large elastic deformations that causes fatigue which in turn can lead to shortened service life due to for example spalling. When the contact is starved of lubricant, or when running conditions do not allow for a hydrodynamic action that fully separates the surfaces, the risk for plastic deformation increases. x

y

Roller width

Figure 1.2: Elliptic contacts (x/aξ )2 + (y/bξ)2 = 1 If bξ exceeds the minimum width of the raceway and the roller, then the contact will be truncated which could lead to increased stresses in the material; at least for the unlubricated contact. In the case of a contact with b ξ /aξ > 4 but with bξ still less than the minimum width of the raceway and the roller, the centerline in the rolling direction can be approximated to a line contact, Evans et al. [1].

1.2 Lubrication regimes The lubrication regimes in EHD lubrication are commonly divided into: Boundary Lubrication (BL), Mixed Lubrication (ML) and Full Film Lubrication (FL). In the BL regime the lubricant’s hydrodynamic action is negligible and the load is carried directly by surface asperities or by surface active additives. In the ML regime the load is carried by the

1.3. SURFACE TOPOGRAPHY

5

lubricant’s hydrodynamical action, the surface active additives and/or directly by surface asperities. When the hydrodynamic action of the lubricant fully separates the surfaces and the load is carried totally by the lubricant film the contact enters the FL regime. In terms of traction the modified Stribeck curve shown in Fig. 1.3 gives a good visualization of the different regimes.

ML

FL

ξ

BL

ln (µu/p)

Figure 1.3: A modified Stribeck curve for EHD lubrication

1.3 Surface topography Depending on the lubrication regime, the surface topography will have a different influence on operation. In the BL regime the topography is preferably chosen so that friction is reduced but without increasing the rate of wear. In the ML regime the objective of the topography is to support the hydrodynamic action of the lubricant, to enable bonding of the surface active additives and to minimize friction in the contact spots without increasing wear. In the FL regime, traction may be reduced by carefully chosen topographies and, even though there is no direct contact, the topography must also prevent fatigue which lead to excessive wear in from of spalling.

1.4 Objectives This work concerns the contact between rough surfaces; with or without a separating lubricant in between. The main objective was to investigate how the surface topography influences different contact conditions. The long term goal is to develop engineering tools, to predict different contact conditions, that takes into account the surface topography. For this purpose two different numerical methods were developed. A method of surface characterization was also developed, to be used to facilitate the interpretation of the data acquired when using deterministic roughness models. The FL regime was initially set in focus and a theoretical model based on Reynolds equation for line contacts was considered. The model was to account for transient effects,


6 i.e., a time dependent model, and include: • Linear elastic surface materials • Surface roughness, measured or modeled • Non-Newtonian rheology

The Reynolds equation derived by Conry et al. [2] was chosen to handle the non-Newtonian rheology. A coupled solver was developed in which pressure and film thickness are solved simultaneously. The numerical method uses a multilevel technique to facilitate the solution process. A schematic picture of the theoretical model is given in Fig. 1.4.

Load R1 E1 ,ν1 µ0 ,α

hy

rap

g po

To

E2 ,ν2 R2

Figure 1.4: A schematic picture of the theoretical FL model used in this work The objective set for this thesis was to develop, verify and use this tool for systematic investigations of the local effects caused by elementary surface features passing each other inside the lubricated conjunction. To complement this work, a contact mechanics model was considered in order to be able to investigate the influence of surface topography for contacts operating in the BL regime. The theoretical model is schematically depicted in Fig. 1.5. The material model used was initially a linear elastic one. However, to be able to simulate running-in the material model was extended to allow for linear elastic-perfectly plastic materials. The specification of this model can be summarized as follows:

1.5. OUTLINE OF THIS THESIS

7

Load

H1

R1 E1 ,ν1

y

ph

gra

po To

E2 ,ν2

H2 R2

Figure 1.5: A schematic picture of the theoretical BL model used in this work

• 2D and 3D contacts • linear elastic-perfectly plastic surface material • Surface roughness, measured or modeled A model accounting for the energy dissipation due to plastic deformation was also considered, see Tian and Bhushan [3]. The numerical solution technique is, however, based on the model of Stanley and Kato [4] with some modification necessary to account for plastic deformations. The BL regime is here modeled as the unlubricated, contact mechanics problem, accounting for plastic deformations, and include the possibility of simulating the running-in process of real surfaces.

1.5 Outline of this thesis The thesis is based on the three papers A, B and C found in part II. In Part I the content of these papers is extended with more a complete theoretical description and derivations and the inclusion of material not addressed in the papers. Part I is written in the form of a monograph.

8


Chapter 2, of this thesis introduces the topic of FL including the theory of the Reynolds approach used in this work. It includes a shortened derivation of the non-Newtonian Reynolds equation according to Conry et al. [2]. The full non-linear integro-differential system of equations including expressions for compressibility, viscosity, elastic deflection, film thickness and force balance are also described. The numerical ‘Block-Jacobi’ method used to solve the coupled system consisting of the Reynolds equation and the film thickness equation may also be found here. A multilevel technique used to facilitate the numerical solution process of the one dimensional Poisson problem and a possible modification to facilitate the numerical solution of coupled systems are then briefly described. The line contact problem is restated in dimensionless form and the discrete analog is derived. Finally, a description of the numerical solution process is given. The theory that constitutes the basis for the dry elasto-plastic contact method is found in Chapter 3 together with an introduction to the contact mechanics problem in tribology. Different models of the elastic contact are then discussed in connection to the present model. A description of how spectral analysis can be used in order to determine the elastic deformation integral is then given. The dimensionless and discrete formulation and the numerical solution process, for the dry elasto-plastic method, may be also found in this chapter. The interpretation of the data acquired when using a deterministic model, simulating either FL or the dry contact, may be facilitated by making use of the surface characterization technique developed in Chapter 4. This technique is based on truncated Fourier series and is applied to a sample topography. In Chapter 5, the Reynolds based approach described in Chapter 2 is compared with a Computational Fluid Dynamics (CFD) based approach. This comparison consist of simulations of one stationary and two different transient problems. The results are encouraging from several viewpoints: verification of the codes, the possibilities to further develop the CFD approach given by Almqvist and Larsson [5], and the justification of using a Reynolds approach under the running conditions chosen. In this comparison, the roughness has the form of a single ridge passing through the EHD lubricated conjunction. The topic of roughness in the FL regime more thoroughly addressed in Chapter 6. In a sliding contact, where the surface topography of both surfaces is, a continuously changing effective surface roughness occurs. For a line contact, the surface topography is superimposed by elementary surface features, i.e., dents and ridges. When passing each other inside the conjunction, these surface features cause local effects. In this chapter the influence of these local effects on the film thickness and the pressure is the subject of investigation. In Chapter 7 the method for modeling the dry elasto-plastic contact between rough surfaces is used to simulate the contact conditions of four different topographies. The simulations are restricted to two-dimensions and profiles taken from the measurements of four surfaces form the basis of the load parameter study performed. The content of the second part of this thesis, part II, is simply the three published papers that form the basis of the thesis. These are as published except for some additional corrections and proofreading.

Chapter 2

Full Film EHD Lubrication Due to increasing demands on performance the lubricant film thickness in EHD lubricated contacts has decreased over the years, and will probably continue to do so. However, the surface machining processes are financially constrained which slows down the decrease in surface roughness and sometimes even causes an increase in commercially available machine elements. As a result, the ratio of film thickness to surface roughness will continue to decrease, which in turn will affect the performance of the EHD lubricated contact. Performance increases must also be matched with improved reliability hence a detailed understanding of EHD lubrication of rough surfaces deserves attention. In this chapter, the theoretical basis of the EHD lubrication model is described. First the adopted deterministic roughness model is described along with previous work done using such a model. A derivation of the one-dimensional non-Newtonian Reynolds-Eyring equation according to Conry et al. [2] is also given. The integro-differential problem for the isothermal EHD lubricated line contact consist of the Reynolds equation derived, an equation for the film thickness including the elasticity of the rollers, an equation for force balance and two semi-empirical equations for the fluid density and the fluid viscosity. These equations are restated in dimensionless form and then discretized. The coupled solution method, which allows for simultaneous solution of pressure and film thickness, referred to as the ‘Block-Jacobi’ method, is explained. A brief overview of the multilevel technique used to facilitate the solution method is also given.

2.1 Deterministic roughness models There are many ways to model surface topographies. Mathematically, the models are either statistic or deterministic. Throughout this work the latter approach is used. In many theoretical studies of EHD lubrication, using deterministic models, one surface is considered smooth and the other as being rough. This is a suitable approximation when modeling a rolling contact or a contact in which the roughness of one of the surfaces is of minor importance. Because of the assumption of parabolic surfaces (Hertzian theory) it is possible to simply sum the roughness of the two contact surfaces and model the contact of the effective roughness and a perfectly smooth surface. This approach will be referred to as ‘one-sided roughness’. 9

CHAPTER 2. FULL FILM EHD LUBRICATION

10

Many extensive, systematic studies have been carried out using the one-sided roughness assumption. Based on the knowledge that the surface and hence roughness deform inside the EHD lubricated conjunction, Lubrecht et al. [6] investigated how sinusoidal Fourier components deform as a function of wavelength and the contact operating conditions, including slide-to-roll ratio. They combined the results from the numerical simulations of the Newtonian line contact in a single formula,

where

Ad 1 = ˜ ˜2 Ai 1 + 0.125∇1 + 0.04∇ 1

(2.1)

3/4 1 = Λ M √1 , ∇ b L1/2 s

(2.2)

that describes the relationship between the reduced amplitude A d /Ai and the initial wavelength Λ, the Hertzian half width b, Moes parameters M and L and the slide-to-roll ratio s. The relation described by Eq. (2.1) is visualized in Fig. 2.1. According to Eq. (2.2), 1

Ad /Ai

0.8

0.6

0.4

0.2

0 0.1

1

10

100

˜1 ∇

1 Figure 2.1: A d /Ai as a function of ∇ a shorter wavelength (Λ) or larger slide-to-roll ratio (s) will lead to a shorter reduced di 1 and thus a less reduced amplitude A d /Ai . In such cases the mensionless wavelength ∇ risk of film breakdown is high. The part of Eq. (2.2) describing the influence of: • applied load (w), • surface mean velocity ((u 1 + u2 )/2), • reduced radius of curvature of the contacting bodies (R x ), • the effective modulus of elasticity (E ), • bulk viscosity (µ 0 ), • the pressure-viscosity coefficient (α)

2.2. GOVERNING EQUATIONS

11

concealed in M, L and b will also affect the reduction in roughness amplitude. If the 1 is expressed in the physical input parameters above: parameter ∇ 1 = ∇

pi1/2 Λ w1/4 3/4

2 Rx

1/2

E 1/4 µ0

1/2

us

α1/2 s1/2

,

(2.3)

it can be seen that a decrease of w, or an increase in any of the parameters u s , Rx , E , µ0 , α will lead to a decreased reduction in roughness amplitude. In an EHD lubricated contact, that is subjected to sliding, the surface topographies of both the surfaces are of significance and a continuously changing effective surface roughness occurs. A number of investigations have been carried out on this subject, for example Evans et al., Tau et al., Chang, Venner and Morales-Espejel [1, 7, 8, 9]. The one-sided roughness approach is not valid in such situations and a ‘two-sided’ roughness treatment is needed. The topic of two-sided roughness is addressed in this work and results of an extensive parameter study may be found in Chapter 6.

2.2 Governing equations A modified Reynolds equation, based on on the Eyring theory of non-Newtonian flow, is derived in one dimension. Johnson and Tevaarwerk [10] proposed the nonlinear constitutive equation for a lubricant under isothermal conditions given by Eq. (2.4). τ 1 dτ τ0 ˙γ = + sinh , (2.4) G dt µ τ0 where G is the shear modulus of elasticity of the fluid, τ is fluid shear stress, µ is the dynamic viscosity, and τ 0 is the Eyring shear stress. This is a Maxwell rheological model where the total shear strain rate is the sum of an elastic term and a nonlinear viscous term based on the Eyring’s theory of viscosity. The modified Reynolds equation is derived from on the Eyring equation, the nonlinear viscous portion of Eq. (2.4), and under the assumption of plain strain rate. Fig. 2.2 shows a fluid element in the thin lubricating film between two solids. (It should be noted that the velocity of the solids at the surface is approximated by the velocity u in the x-direction.) In this case the equation of equilibrium is the x-direction takes the form: ∂p ∂τxz = , (2.5) ∂z ∂x where p is fluid pressure. If the film thickness is denoted by h (x), then according to Fig. 2.2, 0 ≤ z ≤ h (x) and z = z − z1 . Assuming p = p (x), integration of Eq. (2.5) with respect to z yields: dp (2.6) τxz = τ1 + z , dx where τ1 is the shear stress at surface 1. Substituting Eq. (2.6) into the constitutive equation (˙γ = (τ0 /µ) sinh(τxz /τ0 )) yields: τ1 + z dp du τ0 dx = sinh . (2.7) γ˙ = dz µ τ0


12 z

u2

z

z

z1 u1 x

x

Figure 2.2: A description of the line contact region Assuming the velocity of the lower surface (z = 0) is u1 and that the viscosity does not vary across the film ,i.e., µ = µ (x), integration of Eq. (2.7) with respect to z gives the following expression for the velocity profile in the x-direction:

z τ1 + s dp τ0 dx sinh ds, (2.8) u x, z = u1 + τ0 0 µ which after evaluation of the integral expression becomes: τ1 + z dp τ20 τ1 dx − cosh . u x, z = u1 + dp cosh τ0 τ0 µ

(2.9)

dx

Introduction of the midplane shear stress: τm = τ1 +

h dp 2 dx

(2.10)

and a dimensionless function, defined as: Σ=

h dp 2τ0 dx

makes it possible to rewrite Eq. (2.9) as: τ0 h 1 2z τm τm u x, z = u1 + −Σ 1− −Σ . cosh − cosh 2µ Σ τ0 h τ0

(2.11)

(2.12)

Application of the boundary condition u (x, h) = u 2 and utilizing hyperbolic relations gives: τm τm τ0 h u2 = u1 + + Σ − cosh −Σ cosh 2µΣ τ0 τ0 (2.13) τm τ0 h = u1 + sinh sinh Σ. µΣ τ0


13

Eq. (2.13) can be rearranged to allow for determination of the midplane shear stress as follows: τm µ (u2 − u1) Σ . (2.14) sinh = τ0 τ0 h sinh Σ The mass flux per unit width, M (x), is defined as: M (x) =

h 0

ρu x, z dz .

(2.15)

After substitution of Eq. (2.12) into Eq. (2.15) and integration, the expression for mass flux per unit width becomes: τm τm τ0 h h + Σ − sinh −Σ − M (x) = ρu1 h + ρ sinh 2µΣ 2Σ τ0 τ0 (2.16) τm −Σ . h cosh τ0 Expanding and rearranging the terms in the brackets of Eq. (2.16) gives τm τ0h2 sinh Σ M (x) = ρu1 h + ρ sinh + 2µ Σ τ0 τm sinhΣ − Σ coshΣ cosh . Σ2 τ0

(2.17)

Substitution of Eq. (2.14) into Eq. (2.17) yields: M (x) =

=

(u1 + u2 ) ρh + 2 ρτ0 h2 sinhΣ − Σ coshΣ τm cosh 2 2µ Σ τ0 (u1 + u2 ) ρh + 2

=

(2.18)

τm ρτ0 sinhΣ − Σ coshΣ Σ cosh 2µ Σ3 τ0 (u1 + u2 ) ρh + 2 τm dp ρh3 3 (sinhΣ − Σ coshΣ) cosh 12µ Σ3 τ0 dx h2

The dimensionless quantity inside the brackets of Eq. (2.18) is the non-Newtonian slip factor, which after using the hyperbolic relation cosh 2 (x) − sinh2 (x) = 1, and substitution of Eq. (2.14) is given by:

2 µ (u2 − u1) Σ 3 (sinh Σ − Σ coshΣ) 1+ . (2.19) S (x) = Σ3 τ0 h sinh Σ


14

The equation of continuity takes the following form: d dx

h(x) 0

ρu x, z dz = 0.

(2.20)

or in terms of the mass flux per unit width: d M (x) = 0, dx

(2.21)

thus requiring constant mass flux along the x-direction. Substitution of Eq. (2.18) into Eq. (2.21) using the expression of the non-Newtonian slip factor given by Eq. (2.19) finally yields the stationary one-dimensional ReynoldsEyring equation: d ρh3 dp (u1 + u2) d S (ρh) . (2.22) = dx 12µ dx 2 dx The representative lubricant stress τ 0 , characterizes the transition from Newtonian to non-Newtonian fluid behavior. When using the Eyring model an infinitely large τ 0 characterizes a Newtonian fluid, and by using L’Hospital rule it can be shown that S approaches unity. In this case Eq. (2.22) reduces to the conventional Reynolds equation. The value of S is always equal to or greater than one, which means that the effective viscosity µ/S is always less or equal to µ. It is possible to restate Eq. (2.22) in transient form, incorporating squeeze effects. The one-dimensional transient Reynolds-Eyring equation reads: ∂ ∂ ρh3 ∂p (u1 + u2) ∂ S (ρh) + (ρh). (2.23) = ∂x 12µ ∂x 2 ∂x ∂t Eq. (2.23) combined with an equation for the film thickness, an equation of force balance, and two semi-empiric equations for the density and the viscosity respectively, constitutes the basis of the EHD lubricated line contact problem studied in this work. The equation for film thickness is given by: h (x,t) = h0 (t) +

x2 + d (x,t) + ψ1 (x,t) + ψ2 (x,t) , 2Rx

(2.24)

where h0 is an integration constant, R x is the reduced radius of curvature in the x-direction given by 1/R x = 1/Rx1 + 1/Rx2 , d is the elastic deformation of the contacting solids and ψ1 and ψ2 the topography of surfaces 1 and 2 respectively. The elastic deformation in one-dimension is given by: d (x,t) = −

4 πE

∞ −∞

ln |x − x|p x ,t dx ,

(2.25)

where E is the effective modulus of elasticity given by 2/E = (1 − ν21 )/E1 + (1 − ν22 )/E2 . There are several different ways to model viscosity. In this thesis, the Roelands expression is used throughout and is given by: p zvisc αP0 −1 + 1 + . (2.26) µ (p) = µ0 exp zvisc P0

2.3. THE BLOCK-JACOBI METHOD

15

where α, P0 and zvisc are mutually dependent: zvisc (ln (µ0 ) + 9.67), α

P0 =

(2.27)

The Dowson and Higginson equation is the semi-empiric expression used to model the density C1 + C2 p ρ (p) = ρ0 . (2.28) C1 + p At all times, the force balance condition (2.29)

I

p (x,t) dx = w (t)

(2.29)

determines the integration constant h 0 (t). Also, the cavitation condition (2.30) is used to ensure that all negative pressures obtained during the solution process are removed. p ≥ 0.

(2.30)

2.3 The Block-Jacobi method In this section a generalized form of the Jacobi method, taking into account the coupling of the integro-differential system consisting of the Reynolds- and the film thickness- equation, is derived. To increase comprehension, the ordinary Jacobi method is revisited first.

2.3.1 Ordinary Jacobi Let L (x) be general non-linear discrete operator L : R n → Rn . One way to solve the problem: ⎤ ⎡ ⎡ ⎤ 0 L1 (x) ⎥ ⎢ .. ⎥ ⎢ L (x) = 0 ⇔ ⎣ ... (2.31) ⎦ = ⎣ . ⎦, Ln (x)

0

is to apply the Newton-Raphson algorithm. In this case the determining system for the changes, δx, is ⎤⎡ ⎡ ∂L ⎤ ⎤ ⎡ 1 1 · · · ∂L δx1 L1 (x) ∂x1 ∂xn ⎥ ⎢ . ⎥ ⎥ ⎢ . . . . .. ⎥⎢ (2.32) Jδx = −L (x) ⇔ ⎢ ⎦ = − ⎣ .. ⎦. . . ⎦ ⎣ .. ⎣ .. ∂Ln ∂Ln δxn Ln (x) ··· ∂x1

∂xn

Unfortunately, the need to compute the full Jacobian matrix J of the operator L and of solving the full linear system for highly resolved problems leads to an algorithm of high complexity. In the case of rough EHD lubrication, high resolution is needed and the Newton-Raphson algorithm is therefore not suitable. However, if the operator is either linear by nature or in a linearized form, i.e., L (x) ≡ Ax − b

(2.33)


16

where A is a n × n matrix and x and b are n × 1 matrices, then the ordinary Jacobi method could be applied. If ⎤ ⎡ a11 · · · a1n ⎥ ⎢ . . .. A = ⎣ ... (2.34) ⎦ . . an1

···

⎡

and

a11 ⎢ .. diag(A) = ⎣ . 0

ann ··· .. . ···

⎤

0 .. .

⎥ ⎦,

(2.35)

ann

then the changes δx are obtained by solving the following diagonal matrix system diag(A) δx = −L (x) .

(2.36)

Since only the diagonal is used, the complexity of the Jacobi method is far less than that of the Newton-Raphson method.

2.3.2 Coupled systems of equations Let us assume that a non-linear operator L is the discrete representation of a coupled system of two analytical equations, i.e., L1 (v1 , v2 ) , (2.37) L (v1 , v2 ) ≡ L2 (v1 , v2 )

where

L1 (v1 , v2 ) L2 (v1 , v2 )

=

f1 f2

.

(2.38)

and (v1 , v2 ) represents the two dependent variables of the system. If this system is linearized at a point ( vˇ 1 , vˇ 2 ) it is possible to write it in general matrix form as: f A11 (ˇv1 , vˇ 2 ) A12 (ˇv1 , vˇ 2 ) v1 L (v1 , v2 ) ≡ − 1 . (2.39) A21 (ˇv1 , vˇ 2 ) A22 (ˇv1 , vˇ 2 ) v2 f2 If the ordinary Jacobi method is applied to the system L (v 1 , v2 ) = 0, the solution would be obtained by solving Eq. (2.36) with v1 , (2.40) x= v2 A=

A11 (ˇv1 , vˇ 2 ) A12 (ˇv1 , vˇ 2 ) A21 (ˇv1 , vˇ 2 ) A22 (ˇv1 , vˇ 2 )

and b=

f1 f2

(2.41)

.

(2.42)

In this case the coupling between the two equations would not be properly addressed and the solution method would most probably not converge.

2.4. A BRIEF OVERVIEW OF THE MULTILEVEL TECHNIQUE

17

One way to include the coupling between the equations in the solution process is to improve the Jacobi method by introducing a slightly more expensive step δx = [δv 1 , δv2 ]T . Consider the system Eq. (2.39), where A l m are n × n matrices, and v l and fl are n × 1 matrices. A step that address the coupling may be solved from the system J˜ 11 (ˇv1 , vˇ 2 ) J˜ 12 (ˇv1 , vˇ 2 ) L1 δv1 − , (2.43) δv2 L2 J˜ 21 (ˇv1 , vˇ 2 ) J˜ 22 (ˇv1 , vˇ 2 ) where J˜ l m = diag(Al m ) , l = 1, 2 m = 1, 2. Because of the block matrix structure, this method will be referred to as the ‘Block-Jacobi’ method. This makes a simple and not so expensive solution of Eq. (2.43) easy to implement. However, the most important feature of the the ‘Block-Jacobi’ method is that the coupling between the equations is considered.

2.4 A brief overview of the multilevel technique Numerical linear system solvers like the Gauss-Seidel relaxation process and different types of Jacobi iterative processes are often not affordable due to their relatively slow rate of convergence. A multilevel technique applied to solve differential equations is based on the use of a conventional smoother, such as the Gauss-Seidel or Jacobi methods. The smoother itself is conventionally used to solve the problem with specified tolerance and resolution, i.e. using one discrete set of points referred to as a grid for the relaxation process. A decrease in the computational time for such a smoother can be achieved by either accepting a larger error or a lower resolution i.e. a coarser grid. Using multilevel methods it is actually possible to obtain a solution with both specified error and resolution without carrying out all the computational work on the fine grid level corresponding to the requested resolution. Instead the computations are restricted and performed on coarser grids and then prolonged back to the specified resolution. A brief overview of the multilevel technique used in the numerical solution is presented here. Multilevel techniques both for solving differential equations and to facilitate the solution procedure for different integral equations are described in detail by Venner and Lubrecht [11]. The brief description given here is based on a model problem consisting of the one-dimensional Poisson problem given in Eq. (2.46). A two level solution procedure is visualized for the one-dimensional functional operator defined for the Poisson problem. A description of how to modify this solution procedure for a two-dimensional functional operator in order to solve the coupled problem consisting of the Reynolds equation and the film thickness equation is then given.

2.4.1 Grid levels A grid level is defined as a grid with number of nodes determined by the particular level. In this thesis work, the number of grid nodes n ∆x k on the grid level k with mesh size ∆x are obtained as follows: k n∆x (2.44) k = 2 + 1.


18 The mesh size ∆x may be chosen as:

∆x = (b − a)/ n∆x k −1

(2.45)

where the solution domain is given by the interval [a, b]. This means that if the grid level k corresponds to the mesh size ∆x then the grid level k − 1 corresponds to the mesh size 2∆x, i.e., a coarser grid.

2.4.2 Intergrid transfer operators The intergrid transfer operators transfers the functions between different grid levels. The transfer to a coarser grid level is referred to as restriction and the transfer to a finer grid level is referred to as prolongation. These operators can be defined as follows: Definition 1 (Restriction operator) Let def .

∆x v∆x i = v(xi ),

x∆x i ∈ [a = x1 , x1 + ∆x, . . . , x1 + i∆x, . . . , xn = b], ∆x = (b − a)/ n∆x k −1 , where n∆x k is the number of nodes and ∆x is the corresponding mesh size. Then the restriction operator I2∆x ∆x transfers discrete functions v∆x to the next coarser grid with mesh size 2∆x =

b−a n2∆x k−1 −1

cording to ∆x v2∆x = I2∆x ∆x v

Definition 2 (Prolongation operator) Let def .

v2∆x = v(x2∆x i i ), ∈ [a = x1 , x1 + 2∆x, . . ., x1 + i2∆x, . . ., xn = b] , x2∆x i 2∆x =

b−a n2∆x k−1 − 1

,

where n2∆x k−1 is the number of nodes and 2∆x is the corresponding mesh size. Then the prolongation operator I∆x 2∆x transfers discrete functions v2∆x to the next finer grid with mesh size ∆x according to 2∆x v∆x = I∆x 2∆x v

ac-

2.4. A BRIEF OVERVIEW OF THE MULTILEVEL TECHNIQUE

19

2.4.3 Two level solver for the Poisson equation The one-dimensional Poisson problem is given by Eq. (2.46). d2 v = f (x) , v (0) = 0, v (1) = 0. dx2

(2.46)

It is possible to restate the problem in operator form, i.e., L (v)

=

d2 v , dx2

=

f (x) , v (0) = 0, v (1) = 0.

def .

(2.47) L (v)

Using second order finite differences, the problem defined in Eq. (2.47) is discretized into: Li (v)

=

vi−1 − 2vi + vi+1 (∆x)2

, (2.48)

Li (v)

=

fi , i ∈ {2, . . . , n − 1}, v1 = 0, vn = 0.

The iterative solution process of this problem can be facilitated by a multilevel technique based on only two grid levels. For this purpose, the Gauss-Seidel relaxation process is assumed to be the smoother. Based on an initial guess on the finest grid level (v ∆x ) the solution can be obtained according to: • Use a few iterations of the Gauss-Seidel method to obtain an approximate solution v˜i ∆x and then compute the residuals r i∆x = Li∆x v˜∆x − fi∆x , where Li∆x (v) =

∆x ∆x v∆x i−1 − 2vi + vi+1

(∆x)2

and f i∆x is equal to the discrete right hand side function f in Eq. (2.46). 2∆x = ∆x • Restrict the solution and the residuals to the coarse grid, v 2∆x = I2∆x I ∆x v˜i , and rI 2∆x ∆x I∆x ri , where I is the coarse grid index. • Calculate the coarse left hand side L i2∆x v2∆x of the Poisson equation to obtain the corresponding coarse right hand side, f I2∆x = LI2∆x v2∆x − rI2∆x .

• Use a few iterations the Gauss-Seidel method to obtain the approximate solution, v˜I 2∆x , to the system LI2∆x v2∆x = fI2∆x . • Calculate the coarse grid error, E I2∆x = v˜I 2∆x − VI , where VI = v˜∆x 2I . 2∆x and correct the solution, vˆ ∆x = v˜ ∆x + • Prolong the coarse grid error, E i∆x = I∆x i i 2∆x EI ∆x Ei . ∆x • Use a few Gauss-Seidel iterations to obtain the approximate solution, vˆi , on the fine grid.

20


Note the FAS right hand side, given by f I2∆x = LI2∆x v2∆x − rI2∆x on levels j where def . j < k and k is the finest level, is not equal to the discretized right hand side f I2∆x = f x2∆x i of Eq. (2.48). This two level multilevel technique can be extended to a more general multilevel technique, i.e., the multilevel V-cycle containing m levels. That is, the restriction process continues until the coarsest level is reached and then the prolongation process corrects and transfers the solution back to the finest level. (Venner and Lubrecht [11]).

2.4.4 Two-dimensional functional operators The coupled integro-differential problem that consists of the Reynolds equation given by Eq. (2.23), and the film thickness equation given by Eq. (2.24) viz. ∂ ρh3 ∂p ∂ (u1 + u2) ∂ S (ρh) − (ρh) = 0, − ∂x 12µ ∂x 2 ∂x ∂t (2.49) 2 x h (x,t) − h0 (t) − d (x,t) = + ψ1 (x,t) + ψ2 (x,t) , 2Rx is two-dimensional in the functional space A 2 . For example, the two coordinates x and y can be thought of as one point (tuple) in R 2 , i.e, (x1 , x2 ) ∈ R2 . In the same way the pressure p and the film thickness h is a tuple in A 2 , i.e, (p, h) ∈ A 2 . This means that if Eq. (2.49) is restated in operator form, the operators map functions from A 2 to R, i.e., L : A 2 → R. In comparison, the one-dimensional functional operator previously defined for the onedimensional linear Poisson problem map functions from A to R, i.e., L : A → R. The operator system equivalent to Eq. (2.49) can be stated as L1 (v1 , v2 )

=

f1 (x,t) , (2.50)

L2 (v1 , v2 )

=

f2 (x,t) ,

where (v1 , v2 ) represents (p, h), the operators L 1 and L2 are defined as the left hand x2 + ψ1 (x,t) + ψ2 (x,t). Note that the sides of Eq. (2.49) and f 1 (x,t) ≡ 0 and f 2 (x,t) = 2R x deformation, (d (x,t)), is a part of the operator, (L 2 ), because of the explicit dependence of pressure, here represented by v 1 . The multilevel technique applied to the solution of the one-dimensional Poisson problem can be modified to allow for two-dimensional functional operators. In this case, the smoother (section 2.3) is applied to the system Eq. (2.50) which after a few iterations produce the approximative solution (v˜ 1 i , v˜2 i ). The residuals rl i are then calculated according to: ∆x − fl∆x rl∆xi = Ll i v∆x 1 , v2 i , (l = 1, 2), to allow for the determination of the FAS right hand sides 2∆x 2∆x v2∆x − rl2∆x fl2∆x 1 , v2 I = Ll i I on the next coarser grid level. The coarse grid errors (E l I ) are then calculated by 2∆x El2∆x I = v˜l I − Vl I ,

2.5. DIMENSIONLESS FORMULATION

21

where Vl I = v˜∆x l 2I . These errors are prolonged in order to correct the solution (v˜ 1 , v˜2 ). Together with the description for the two-level technique for the one-dimensional Poisson problem, this description outlines the multilevel technique used to facilitate the numerical solution of Eq. (2.49).

2.5 Dimensionless formulation The smooth Newtonian EHL line contact is dependent on the six input parameters w, E , us = (u1 + u2 ) , Rx α, η0 . For a smooth Reynolds-Eyring approach two extra input parameters are used; (u1 − u2 ) which could be introduced as a slide-to-roll ratio s = 2 (u 1 − u2 ) / (u1 + u2) and τ0 for the Eyring shear stress. For a rough Reynolds-Eyring approach the additional input parameters needed depend on how the surface topographies are modeled, e.g., for a single ridge/dent modeled by an amplitude and a wavelength only there are 8 + 2 = 10 input parameters. In order to restate the EHL line contact problem in dimensionless form the following set of dimensionless parameters where introduced, (Hamrock [12]): T = t/ (2b/us) ,

X = x/b,

H = h/ b2 /Rx , P = p/ph , ρ¯ = ρ/ρ0, where

b=

(2.51)

µ¯ = µ/µ0,

2w 8wRx , ph = = πE πb

wE 2πRx

(2.52)

If this set of dimensionless parameters is used to restate the EHD lubricated line contact problem in dimensionless form, it becomes clear that the smooth incompressible problem is, mathematically, a two-parameter problem, see Moes [13], i.e., the solution is dependent on the parameters M and L, viz.

2b3 ph 21/2 w π = 1/2 . (2.53) M= 1/2 1/2 4Rx µ0 us Rx E 1/2 µ us 0

2π E 3/4 µ0 us α = . (2.54) 1/4 M 21/4 Rx only. In this case, the Moes parameters given by Eq. (2.53) and Eq. (2.54) enables systematic parameter studies. That is, there is a one-to-one mapping between a given Moes pair (M, L) and parameters such as p max , hmin , hc , etc. It is also possible to use these parameters as determining factors for any two other independent parameters. Using the transformation given by Eq. (2.51) also prevents the introduction of truncation errors when solving the problem numerically. That is, because of the particular choice of transformation the values of the dimensionless parameters are restricted so as not to assume extremely large nor extremely small values. For example, the Hertzian pressure in a rolling element bearing is typically O 109 Pa whereas the film thickness is typically O 10−9 − 10−6 m and in terms of dimensionless quantities these parameters are O 100 . L = αph

1/4 1/4


22

In this work, the compressible solution is sought. This introduces an expression for the pressure-density relation which affect the solution. The density at ambient pressure, (ρ 0 ) does, however, not affect the solution The determining parameters are here represented by λ and α¯ which in terms of input parameters is given by: λ=

3π2 µ0 us Rx E 6us R2x µ0 = ph b3 w2

and α¯ = αph =

αwE . 2πRx

(2.55)

(2.56)

Using the dimensionless parameters given by Eq. (2.51), Eq. (2.52) and the expressions Eq. (2.55) and Eq. (2.56) the resulting dimensionless Reynolds-Eyring equation (Eq. (2.23)) is transformed into: 3 ¯ ∂P ∂ ∂ ∂ ρH ¯ + ¯ S (ρH) (ρH) . (2.57) =λ ∂X µ¯ ∂X ∂X ∂T which is, explicitly dependent on the first determining parameter λ. Applying the same transformation, the dimensionless form of the film thickness equation (Eq. (2.24)) reads: H (X, T ) =

H0 (T ) +

X2 1 − 2 π

∞ −∞

ln |X − X |P X , T dX (2.58)

+

¯ 2 (x,t) ¯ 1 (x,t) + ψ ψ

which, if ψi ≡ 0, i = 1 ∧ 2 is an equation totally independent of input parameters. Note that the dimensionless parameter H 0 (T ) has, in addition to the dimensionless form, a contribution from the dimensionless restated elastic deformation, i.e., 1 ∞ ln (b) P (X, T ) dX. H0 (T ) = h0 (T ) / b2 /Rx − π −∞

(2.59)

The force-balance equation Eq. (2.29) transforms into:

I

P (X, T ) dX =

π . 2

(2.60)

which is also an equation totally independent of input parameters when restated in dimensionless form. The Roelands expression for the viscosity depends on the determining parameter α¯ and the parameter p h when transformed into dimensionless form, viz. ¯ 0 αP Pph zvisc µ¯ (P) = exp −1 + 1 + . (2.61) ph zvisc P0 where

ph zvisc (2.62) (ln(µ0 ) + 9.67), α¯ or P0 is a determining parameter depending on which one of them P0 =

Note that, either zvisc that is specified.

2.6. DISCRETE FORMULATION

23

The dimensionless form of the Dowson-Higginson expression for the density is given by ρ¯ (P) =

C1 + C2 Pph , C1 + Pph

(2.63)

where the constants C1 and C2 is chosen so that ρ −→ 1.34 as p −→ ∞. Together with the cavitation condition P ≥ 0, (2.64) this completes the set of dimensionless equations used in a Reynolds-Eyring approach for the EHD lubricated line contact problem.

2.6 Discrete formulation The set of equations yielding the EHD lubricated line contact problem given in section 2.5 are here discretized. The discrete form of the time-dependent Reynolds Eq. (2.57) used in this work can be written as: 1 k k k k k k k εi− 1 Pi−1 − εi− 1 + εi+ 1 Pi + εi+ 1 Pi+1 − 2 (∆X)

2

2

2

2

k k λ k k ¯ ¯k ¯k 2∆X 3ρi Hi − 4ρi−1 Hi−1 + ρi−2 Hi−2 − λ ∆T

where εki± 1 2

(2.65)

k−1 = 0, H ρ¯ ki Hik − ρ¯ k−1 i i

3 εki + εki±1 ρ¯ ki Hik Sik k , εi = = . 2 µ¯ ki

(2.66)

Which is of the second order in space (index i) and of the first order in time (index k). The following discrete form of the film thickness equation Eq. (2.58) is used: Hik − H0k − where Ki j

=

(Xi )2 1 + 2 π

∑ Ki j Pkj + ψ¯ ki = 0,

(2.67)

j

∆X ∆X Xi − X j + ln Xi − X j + −1 − 2 2 ∆X ∆X ln Xi − X j − − 1 . Xi − X j − 2 2

(2.68)

The discretization technique used for the force-balance equation (2.60) is the trapetzodial rule. In the Roelands equation (2.61) and the Dowson-Higginson equation (2.63), used for determining the viscosity and density respectively, the dimensionless pressure P is simply replaced by its discrete representation Pi .


24

2.7 Solution method Equations (2.65) and (2.67) can be written in operator form (Eq. (2.38)) if the k − 1 term in Eq. (2.65) is moved to the right hand side of the equation. If the coefficients in Eq. (2.38) are linearized using the previously determined values of pˇ k and hˇ k obtained from any relaxation process, then it would be possible to rewrite the system in the form of Eq. (2.39), viz. k k−1 A11 pˇ k , hˇ k A12 pˇ k , hˇ k p f1 = . (2.69) hk A22 pˇ k , hˇ k A21 pˇ k , hˇ k fk2 where the time dependence (k) is written out explicitly. As an example, for second order discretization in both the X- and the T - directions, the matrix A 11 is a tri-diagonal matrix with one sub and one super diagonal, A 12 is a tri-diagonal matrix with two sub diagonals. The matrix A21 depends on the discretization order of the integral expression for the elastic deformation and is a full matrix. The diagonal matrix A 22 is however, independent of the discretization order. If the Block-Jacobi method is applied to solve Eq. (2.69) it is straightforward to rewrite the system in indexes form, i.e., in the form of a 2 × 2 matrix system ai bi pi f1 i = , (2.70) ci d i hi f2 i with coefficients: ai

=

bi

=

k k ε + ε 1 1 i+ 2 (∆X)2 i− 2 3 1 + −λ ρ¯ ki 2∆X ∆X 1

−

ci

=

∆X ln

bi

=

1

∆X 2

(2.71)

and f1 i

=

− −

f2 i

=

1 2

(∆X) λ

k k εki− 1 Pi−1 + εki+ 1 Pi+1 2

2

λ k−1 k−1 k k − 4ρ¯ ki−1 Hi−1 − ρ¯ ki−2Hi−2 ρ¯ Hi 2∆X ∆T i

H0k +

Xi2 1 − 2 π

(2.72)

∑ Ki j Pkj − ψ¯1ki − ψ¯2ki . j =i

The solution of this system can be vectorized and thus very easily implemented in any higher level programming language. Unfortunately, the nature of the EHD lubricated line contact problem introduces some difficulties that have to be dealt with. The solution process consist of an inner relaxation loop for Eq. (2.70) and an outer loop to converge the

2.7. SOLUTION METHOD

25

force-balance criterion. When the solution process is to be facilitated by the multilevel technique described in section 2.4.4 the force-balance condition is relaxed only at the coarsest grid level. The convergence of the force-balance equation is very important in order to obtain sufficient accuracy of H 0 within each time step. This can be seen from the pressure solution of the coupled system Eq. (2.73), viz. Pi =

di f i − bi gi , a i d i − b i ci

(2.73)

where the nominator contains the time dependence of the Reynolds equation. The nominator of Eq. (2.73) expresses numerically the time derivative, in dimensionless form, of the film thickness and thus the time derivative of H 0 (T ) . That is, a small ∆T will enhance the influence of an error in H 0 (T ) and the physical explanation of this behavior is the squeeze mechanism represented by ∂/∂T in Eq. (2.57). Moreover, this is not only important when using a coupled solution method; a serial/segregated solution method will be affected in the same way. At the outlet of the EHD lubricated conjunction the cavitation condition Eq. (2.64) is imposed. Thus, the solution of the Reynolds equation, which admits a negative pressure, is not valid here and the resulting non-Newtonian parameter S will drop to zero. Setting S = 1, for ∂P/∂X = 0 in the cavitated region, leads to an effective viscosity µ/S which equals the ambient viscosity µ 0 . This gives an almost smooth extension of S in the region of zero pressure. The importance of such an extension relates to the Block-Jacobi relaxation of Eq. (2.70), since Eq. (2.64) implies P = 0 in the cavitation zone and that the whole domain is considered at all iterations. In the solution process, this correction successively moves the starting point of the cavitation zone to a certain, converged, point.

Chapter 3

Dry Elasto-Plastic Contact Machining processes aim to produce surface topographies to a given cost, that have a positive influence on contact conditions. For some applications product life may be of greatest importance in which case a machining process must be chosen so that wear is minimized. In other applications, low friction may be of primary importance and one of the tasks of the surface topography might consequently become to minimize the real area of contact. However, this may significantly reduce service life since a small contact area causes a very high pressure and thus an increased risk of wear. This chapter includes theory for both 2D and 3D dry elasto-plastic contacts. A contact mechanics problem can be governed by the minimum potential complementary energy theory which is an approach based on a deterministic roughness model. That is, pressures and displacements can be obtained by solving a minimum value of an integral equation, i.e., a variational problem, see Kalker [14], for rough surface contact mechanics. Combined with an equation of force-balance, the equation for potential complementary energy including constraints provides the basis for the solution of the dry elasto-plastic contact. It should be mentioned that it is possible to use pressure and gap constraints together with the Boussinesq equation, i.e., Eq. (3.5) to solve the elastic problem, Lubrecht and Ioannides [15]. For both the models described above it can also be argued that, if the pressure constraint includes an limiting value p ≤ H s , where Hs is the hardness of the softer surface material, and that the load is such that only a very small fraction of the resulting pressure profile is truncated by this limiting value, then this is one way of modeling the dry contact assuming elastic perfectly plastic surface material. This way of modeling elastic perfectly plastic surfaces does not, however, allow for computation of the resulting plastically deformed surface since it does not include in the energy dissipation due to plastic deformation. Using the variational approach, however, it is possible to account for energy dissipation and thus allow computation of the resulting plastically deformed surface. This model is meant to be used for prediction of parameters determining different contact conditions, such as the real area of contact, and not for determining the sub-surface stress field in detail.

3.1 Statistical roughness models Modeling topographies of real surfaces is a difficult task because of their random structures. Various approaches to modeling the dry contact between surfaces have been re27

28

CHAPTER 3. DRY ELASTO-PLASTIC CONTACT

ported in the literature. Greenwood and Williamson [16] (GW), considered topographies consisting of hemispherically shaped asperities of uniform radius. They assumed that the asperities had a Gaussian distribution, in height, about a mean plane. Greenwood and Tripp [17] extended this approach to handle the contact between two rough surfaces. Further improvements to the GW approach were added by Whitehouse and Archard [18], Nayak [19], Onions and Archard [20], Bush et al. [21, 22] and Whitehouse and Phillips [23, 24]. The input needed for these types of models is deduced from surface measurements. More precisely, all of these models are based on asperity curvature. This parameter depends on the measurement resolution (Poon and Bhushan [25]) and a mean value of the asperity curvature is often a poor approximation since it can vary considerably in between individual asperities. It is, therefore, a difficult task to find the correct input value for the model. Moreover, use of only a few parameters to describe a surface generates a one-tomany mapping, i.e., the same set of parameters can be deduced for surfaces obtained by completely different machining processes.

3.2 Deterministic roughness models In line with improvements in affordable computing power, deterministic approaches to the dry contact problem have been developed. These are becoming more and more realistic as computational speed increases. However, simplified models of the material and/or the topography are still needed. This is necessary to minimize the computing time and to be able to study different effects independently. For example, the Hertzian contact assumes linear elastic, frictionless materials and parabolic surface profiles for which an analytical solution is available. That is, the contact between a single asperity, approximated by a parabolic profile, and a rigid plane may be studied analytically. Westergaard [26] showed that it is possible to solve analytically the contact between a single frequency sinusoidal surface (mathematically: sin (2πnx)) in contact with a rigid plane. On a micro scale these analytical solutions represent quite realistic topographies. However, the coupling to asperities at other scales is neglected and a more complete representation is therefore needed. In this work the total surface content, i.e., all data of the measured or modeled surface topography is used to determine different contact conditions such as the real area of contact. It should be mentioned that a deterministic roughness model benefits from surface characterization (Chapter 4) which can be used tp deduce the data and possibly facilitates the interpretation of the results obtained.

3.3 Numerical solution techniques In addition to the different ways of modeling the mating materials and topography, there are numerous deterministic numerical techniques that can be applied to solve contact problems between rough surfaces. Lubrecht and Ioannides [15] applied multilevel techniques to facilitate the numerical solution of the elastic contact problem. Ren and Lee [27] applied a moving grid method to reduce storage of the influence matrix when the conventional matrix inversion approach is used to solve the contact problem of linear elastic bodies with rough surfaces. Björklund and Andersson [28] extended the conventional matrix inversion approach by incorporating friction induced deformations. This was done by using the assumption of linear elastic material. Ju and Farris [29] introduced an FFT-based method to


29

solve the elastic contact problem. Stanley and Kato [4] combined this FFT-based method with a variational principle to solve both the 2D and the 3D contact problem of rough surfaces. It should be noted that the contact between realistic topographies, under relatively small loads, leads to plastic deformations. See for example, Tian and Bhushan [3] who based their theoretical model on a variational principle for both linear elastic and linear elastic perfectly plastic materials. In contrast to the model given by Stanley and Kato, the present model incorporates linear elastic-perfectly plastic materials and the energy dissipation due to plastic deformation is accounted for. In this way, not only the in-contact topography and the corresponding pressure distribution but also the unloaded plastically deformed topography is obtained.

3.4 Governing equations Assuming frictionless linear elastic contact, the variational problem including constraints to be solved ([14]) is given by (3.1). A general formulation is introduced here that admits both 2D and 3D problems.

1 p de dΩ − p ψ dΩ , min () = min p≥0 p≥0 2 Ω Ω (3.1)

S

p dS = W,

where p is the pressure, d e is the normal, elastic deformation, ψ is the gap between the undeformed surfaces and W is the applied load. The mathematical model given by Eq. (3.1) has one equivalent in the system p = 0,

h > 0,

p > 0,

h = 0,

S

(3.2)

p dS = W,

where h = de + ψ + const, see Lubrecht and Ioannides [15]. Restated, this model may be expressed as: h = de + ψ + const, S

p dS = W,

(3.3)

p ≥ 0. However, some modification of Eq. (3.1) is necessary to solve the contact problem assuming linear elastic - perfectly plastic surface materials, model depicted in Fig. 3.1.


30

σ

σs

Elastic

Perfectly plastic

ε

Figure 3.1: Material model for linear elastic - perfectly plastic surfaces To account for the energy dissipation due to plastic deformations, the variational statement given by Eq. (3.1) has to be restated (see Tian and Bhushan [3]) and the variational approach reads:

1 min () = min p de dΩ − p (ψ − d p) dΩ . 0≤p≤Hs 0≤p≤Hs 2 Ω Ω (3.4)

p dS S

= W,

where Hs is the hardness of the softer material and d p is the normal, plastic deformation. For two elastic half spaces, assuming plain strain, the total normal surface displacement de (x) for a given pressure distribution can be obtained from Eq. (3.5). de (x) =

∞

−∞

K (x − s) p (s) ds + const.

(3.5)

where bold symbols indicate vector valued parameters, e.g., x = x in 2D and x = (x, y) in 3D and the integral kernel K is given by K (x − s) = − and K (x − s, y − t) =

4 ln|x − s| , πE

2 1 πE (x − s)2 + (y − t)2

(3.6)

(3.7)

for 2D and 3D contacts respectively.

3.5 Spectral analysis In the case of a sinusoidal pressure p (x) of one frequency only it is possible to find the exact, closed form solution for the normal surface displacement d e (x), see for example

3.6. DIMENSIONLESS FORMULATION

31

Johnson [30]. In 2D, the relationship between pressure and elastic deformation reads: p (x) = α cos (2πnx/Lx ) ⇐⇒ de (x) =

4 p (x) + const. nE

(3.8)

Thus, the displacement d e (x) is a sinusoid of the same frequency as the pressure but with frequency dependent amplitude. The 3D analog is given by: p (x, y) = de (x, y) =

α cos (2πnx/Lx ) cos (2πmy/Ly ) (3.9) 1 4 √ p (x, y) + const. 2 E n + m2

According to Spectral/Fourier theory it is possible to rewrite the relations given by Eq. (3.8) and Eq. (3.9) in a more general form. The Fourier transformation f of f may be defined for the 1D case as:

1 ∞ f(ωx ) = f (x) γ (ωx , x) dx (3.10) 2π −∞ and the 2D analog as: f(ωx , ωy ) =

1 (2π)2

∞ ∞ −∞

−∞

f (x, y) γ (ωx , x) dx γ (ωy , y) dy,

(3.11)

where γ (ωz , z) = exp(−iωz z) and i is the imaginary unit. The relation re-stated in general form then becomes de (ω) = χ (ω) p(ω) , (3.12) where de and p denotes either the 1D or 2D Fourier transformation of the elastic deformation and pressure respectively and χ denotes the transfer function given by: 4 , nE

(3.13)

2 1 √ , 2 E n + m2

(3.14)

χ= in 1D and χ= in 2D.

3.6 Dimensionless formulation The mathematical model given by Eq. (3.1) is equivalent to the model given by Eq. (3.3). This equivalence may be used to explain the simplified similarity analysis used to form the dimensionless factors for the dry contact problem. The analysis is based on the linear elastic model but the results may be used also for the linear elastic - perfectly plastic surface material model. If the dimensionless for the 2D contact problem forms is stated as X = x/xr , P = p/pr , H = h/hr ,

(3.15)


32

the corresponding dimensionless Eq. (3.3) model is transformed into: 4pr xr H = − πE h r

∞

−∞ ln |X

¯ + const, − S| P (S) dS + h1r ψ

∞

−∞ P (X)

dX =

w xr pr ,

(3.16)

P ≥ 0, where the dimensionless form of d e (for 2D contacts) is written out explicitly. It is quite easy to perform an optimum similarity analysis for this system. That is, since there are only two groups of distinctly different dimensionless functions: 4pr xr w = const1 , = const2 . πE hr xr p r

(3.17)

If the solution of the nominally flat 2D contact problem with rough surfaces is sought it is convenient to use L x which is the length in the x-direction of the modeled roughness as the reference length x r . An arbitrary choice of the constants (const 1 and const2 ) in Eq. (3.17) can be made, for example xr = Lx , pr =

w 4w , hr = , Lx E

(3.18)

so that const1 = 1/π and const2 = 1. In terms of these reference quantities the system Eq. (3.16) transforms into: H = − π1

∞

−∞ ln|X

E ¯ − S| P (S) dS + 4w ψ + const,

∞

−∞ P (X)

dX = 1,

(3.19)

P ≥ 0. This is a one-parameter system determined by the fraction E / (4w) only. That is, all possible solutions, in terms of physical quantities, can be deduced from a parameter study ¯ is held of this particular fraction as long as the dimensionless form of the topography ( ψ) constant.

3.7 Discrete formulation The solution of Eq. (3.4) and Eq. (3.1) for 2D and 3D contacts respectively is described in Section 3.8. Included in this description is the procedure of minimizing , i.e., finding the zero gradient of with respect to the pressure p. In this section the discrete analog to the variational problem, Eq. (3.4) and other necessary discrete equations used in the solution process are given. The elastic deformation Eq. (3.5) is discretized into dei =

∑ Ki j p j , j

(3.20)

3.8. SOLUTION METHOD

33

where Ki j is an arbitrary discretization of K. In the solution process however, the relation Eq. (3.12) is used to facilitate the computation of the elastic deformation. That is, the transfer functions Eq. (3.13) and Eq. (3.14) are discretized and used to calculate the elastic deformation in the frequency domain. The discretization of the transfer function is dependent on the particular discrete Fourier transformation which is being computed by the fast Fourier transformation (FFT) routine used in the solution process. The most straight forward way of determining the transfer function is by simply checking how the sinusoidal functions A1 sin (2πnx) + A2 cos (2πnx), A11 sin (2πnx)sin (2πmy) + A12 sin (2πnx)cos (2πmy) +

(3.21)

A21 cos (2πnx)sin (2πmy) + A 22 cos (2πnx)cos (2πmy), in 1D, 2D respectively, are being transformed into the frequency domain. For example, in 1Dthe Fourier transformation fˆ (omegax) of f (x) = sin (2πnx/Lx ) will be pure imaginary (ℜ fˆ (omegax ) ≡ 0) and different from zero only at ω x = n and ωx = N − n. It is also possible to apply a multilevel technique, i.e. multilevel multi integration/summation to facilitate the computation of the elastic deformation, see [15]. Numerically, the equation defining is a quadratic form, because of the specific relation between the elastic deformation and the pressure, viz. i = ∑ pi dei − ∑ pi ψi + d pi i

=

1 2

i

∑ pi ∑ Ki j p j − ∑ pi ψi + d pi , i

j

(3.22)

i

which restated in matrix form is given by 1 = pKp − p (ψ − dp ) . 2

(3.23)

de f .

Thus the gradient ∇ = ∂/∂p of is given by ∇ =

Kp − (ψ − dp ) (3.24)

=

de + dp − ψ.

3.8 Solution method The relationship Eq. (3.12), which can be interpreted as the complete contact of sinusoidal elastic half spaces, assuming plain strain, together with Eq. (3.4), for the 2D contact, forms the basis of the solution technique used in this work. The solution method used for the 3D contact is based on Eq. (3.12) together with Eq. (3.1) and provides the pure elastic solution only. However, the elasto-plastic model is also valid in the case of a 3D contact.

34


For the linear elastic contact, the numerical solution procedure proposed by Stanley and Kato [4] is adopted. Some modification of this procedure is necessary to solve the contact problem assuming linear elastic - perfectly plastic materials. That is, ∇ has to be modified according to ∇ = d e + dp − ψ to account for the energy dissipated due to plastic deformations. The procedure is then iterated until d e , dp and p satisfies Eq. (3.4), which is the method used here for the 2D contact. Note that the 3D contact could also be solved using the elastic-plastic solution procedure presented here.

Chapter 4

Surface Characterization Investigations where a deterministic surface roughness model is adopted for the simulation of rough FL or dry contacts has the advantage of working with the total surface content compared to working with models having parameters deduced from the topography as input. However, this might also be a major drawback of the model since it involves working with a great amount of data. Surface characterization may be used to deduce significant properties of the surface topography and thus reduce the amount of data used. This is important when working with measurements of surface topographies but might also be useful when working with modeled topographies.

4.1 Theoretical model The limitations of the surface measurement apparatus used may introduce problems when characterizing surfaces. These may lead to different characterizations of the same surface, which in turn may lead to contradictory conclusions as to how the topography influences the contact conditions. However, depending on the parameters sought and the application (rough FL or dry contact), it might be possible to overcome the problem of measurement resolution. For a 2D profile, it is possible to filter the topography through its discrete Fourier series, as shown by ψn (x) =

n

∑ ak sin (2πkx/L) + bk cos (2πkx/L) ,

(4.1)

k=1

where the number of frequencies n ≤ N/2 − 1 and N is the number of sampling points.

4.2 Measured topographies In this work, the form of truncated Fourier series Eq. (4.1) has been used for surface characterization of measured topographies used in dry contact simulations. To give some insight to this type of surface characterization technique, a profile taken from a measured surface was used, Lord and Larsson [31]. The plastically deformed surface profile, loaded with 100 kN and hence a mean pressure of approximately 250 MPa is shown in Fig. 4.1. The upper edge of the filled red region, which displays the plastic deformation, is the orig35

CHAPTER 4. SURFACE CHARACTERIZATION

36 1.5 1

z (µ m)

0.5 0 -0.5 -1 -1.5 0

0.1

0.2 x (mm)

0.3

0.4

Figure 4.1: Plastically deformed surface profile in the case of 100 kN load. inal profile and the black line shows the plastically deformed profile. It can be observed, under these conditions, that the maximum depth of plastic deformation is approximately 0.5 µm.

4.3 Parameter study In this section a parameter study is performed using the profile shown in Fig. 4.1. Fig. 4.2 shows how the real area of contact for both linear elastic- and linear elastic perfectly plastic- materials varies with the number of frequencies (n) used to represent the topography. The number sampling points (N) used were 512 which leads to N/2 − 1 = 255 frequencies. It is clear that the limitations of the apparatus used does not affect the percentage 100 Elastic Elasto-Plastic

Ar /An (%)

80 60 40 20 0

50

100

150

200

250

n

Figure 4.2: Percent real area of contact as a function of the number of frequencies n under 100 kN load. of real area of contact A r /An since a constant value is attained for n ≈ 75 and n ≈ 50 for the elastic and the elasto-plastic surface models respectively, which is significantly smaller

4.3. PARAMETER STUDY

37

than N/2 − 1 = 255. The lowest value of n for which A r /An can be considered as constant, depends on the applied load. However, for the theoretical model to be valid, the ratio Ar /An needs to be small; maybe even smaller than for the case presented in Fig. 4.2. The surface roughness parameters R a , Rz , the skewness (Rsk ) and the kurtosis (R k ) were also considered to be functions of n, in the same way as described for A r . In Fig. 4.3 the average roughness (R a ) values variation with n is shown. It can be seen 0.4

Ra (µ m)

0.3

0.2

0.1

0

Elastic Elasto-Plastic Undeformed 50

100

150

200

250

n

Figure 4.3: The R a as a function of the number of frequencies n under 100 kN load. that the values of R a as a function of n for the undeformed profile and the elastic model are approximately the same and that the deviations between those and the elasto-plastic model are visible but still quite small. The R a -values converges with n approximately as did Ar /An . The extreme parameter R z displayed in Fig. 4.4 shows more spread between the undeformed, the elastic and the elasto-plastic results than the R a -value. The trend, however, is 2.5

Rz (µ m)

2 1.5 1 Elastic Elasto-Plastic Undeformed

0.5 0

50

100

150

200

250

n

Figure 4.4: The R z as a function of the number of frequencies n under 100 kN load. the same but the value of n needed to converge the R z -value is larger.

CHAPTER 4. SURFACE CHARACTERIZATION

38

The skewness (Rsk ) and the kurtosis (R k ) variation with n are shown in Fig. 4.5 and Fig. 4.6 respectively. These values do like A r /An and converge to a constant value for 0.5

Rsk (-)

0

-0.5

-1

-1.5

Elastic Elasto-Plastic Undeformed 50

100

150

200

250

n

Figure 4.5: The skewness (R sk ) as a function of the number of frequencies n under 100 kN load.

4.5 4

Rk (-)

3.5 3 2.5 2 Elastic Elasto-Plastic Undeformed

1.5 1

50

100

150

200

250

n

Figure 4.6: The skewness (R sk ) as a function of the number of frequencies n under 100 kN load. n significantly lower than N/2 − 1 = 255. The differences in spread are explained by the functions determining the skewness k = 3 and the kurtosis k = 4 given in Eq. (4.2). 1 NRkq

N

∑ (ψi − ψavg)k ,

(4.2)

i=1

where N

ψavg = 1/(N) ∑ ψi i=1

(4.3)

4.4. CONCLUSIONS

39

is the mean surface height. The results of this parameter study shows that all parameters investigated are accurately predicted using a surface profile representation with a value of n < N/2 − 1.

4.4 Conclusions The numerical method described in Chapter 3, for the 2D contact between linear elastic perfectly plastic rough surfaces, has been combined with surface characterization. Based on simulations of a surface profile originating from a measured topography the following conclusions can be drawn: • The high frequency content of the measured profile has a negligible influence on the real area of contact. • The average roughness value (R a ) converges to a constant value for n < N/2 − 1. However the value of n is larger than for the real area of contact study. • The average maximum height value (R z ) follows the convergence of the R a value but the spread between the undeformed surface and the resulting R z from the pure elastic and the elasto-plastic models is larger. • The rate of convergence for the skewness (R sk ) and the kurtosis (R k ) is approximately the same as for the real area of contact when considering the parameters as functions of the number of frequencies (n). Due to the functions determining these values, the spread is larger for the skewness values (undeformed, elastic and elasto-plastic).

Chapter 5

Reynolds vs. CFD, A Comparison of Different FL Approaches When simulating full film elastohydrodynamic lubrication (FL), the Reynolds equation is the dominat partial differential equation used for prediction of the fluid flow. Very few attempts have been carried out using the full Navier-Stokes equations. In this chapter the Reynolds approach described in chapter 2 is compared with a CFD approach for simulations of the FL regime. The CFD approach is implemented as an extension to the commercial software CFX 4 in which the full Navier-Stokes equations are used as basis for simulating the FL regime, Almqvist and Larsson [5, 32]. This approach enables numerical simulations outside the contact regions, as well as in the thin film region. The numerical experiments show that, under the running conditions chosen, only small deviations between the two approaches can be observed. The results are encouraging from several viewpoints: verification of the codes, the possibilities of further developments of the CFD approach and the justification of using a Reynolds approach under the running conditions chosen.

5.1 The CFD approach A CFD approach certainly has advantages compared with a Reynolds approach when simulating EHD lubricated line contacts. It is not assumed that the inertia inside the fluid is negligible and the thin film assumption is not applied, which in some cases determine the accuracy of the predicted pressure- and film thickness- profiles. By not utilizing these assumptions on the fluid behavior, it is possible to extend the computational domain to include the flow outside the EHD lubricated conjunction. A disadvantage of a CFD compared to a Reynolds approach is the reduction in efficiency, mainly due to the enlarged size of the problem. For 3D cases, four equations have to be solved in order to obtain the fluid pressure in comparison with the single Reynolds equation. Another disadvantage is that it is not straightforward to modify a commercial CFD-code to solve FL problems. 41

CHAPTER 5. REYNOLDS VS. CFD

42

5.2 CFD - Governing equations Transient simulations of rough surfaces needs special attention. That is, the topographies will in this case travel through the computational domain and the grid must be updated in time. History is of importance and must be accounted for in the equations of momentum and continuity. According to Hawkins and Wilkes [33], the transient and advective, or convective terms need to be corrected to include the volume changes and grid velocities. This may be done by inclusion of the determinant of the metric tensor g, and the advective term to allow for a grid velocity ∂x/∂t, (CFX4.4 [34]). The corrected equation of momentum is given by: √ ∂x 1 ∂ ρu g +∇· ρ u− ⊗ u = ∇ · σ. (5.1) √ g ∂t ∂t √ where g relates a volume in the physical space to a unit volume in the computational space. The grid velocity is included in the advective term ∂x/∂t. The symbol σ denotes the total stress tensor and the symbol ⊗ denotes the vector product u i × u j . The body forces are assumed to be vanishingly small compared to the other terms and can therefore be neglected. This also applies to the equation of continuity given by: √ ∂x 1 ∂ ρ g + ∇ · ρu − = 0. (5.2) √ g ∂t ∂t The equations described above form the basis of the isothermal CFD-approach used in Paper A. In [5] thermal effects are also considered. In this case the energy equation is, also modified to account for rough surface simulations. An Eyring rheological model, developed by Johnson and Tevaarwerk [10], where the elastic term is neglected, was used for the dynamic part of the stress tensor (σ). In this case the constitutive equation takes the form: σ =

−pδ −

2µ ∇ · uδ + µ˙γ, 3

γ˙ =

τ ∇u + (∇u)T = τ0 sinh τ0

(5.3)

.

In the equations given by 5.3, the viscosity (µ) is modeled by the Roelands equation, Eq. (2.26). The Dowson-Higginson expression for the density of the lubricant [35] given by Eq. (2.28) is used to model the density.

5.3 The model problem The choice of model problem was set to be transient, non-Newtonian and to include sliding in order to validate the use of a Reynolds approach under such conditions. A relatively lightly loaded line contact with a Hertzian pressure of 0.6 GPa where a single ridge defined by ψ given in Eq. (5.4) passes through the conjunction under sliding was chosen for the comparison. Isothermal contact conditions were assumed. An Eyring model of the nonNewtonian behavior of the fluid was adopted; having a reduced shear stress to shear rate

5.3. THE MODEL PROBLEM

43

dependence when the pressure and shear rates are high. ψ (x,t) = A 10

−10

x−xc (t) 2 λ

2π (x − xc (t)) cos λ

(5.4)

Three different cases were investigated, one with perfectly smooth surfaces, one with a ridge attached to the fast moving surface and the one with a ridge attached to the slower moving surface, Case 1, Case 2 and Case 3 respectively. Input data for the three different cases can be found in Table 5.1. A relatively large slide-to-roll ratio (s = ±5/3) was

Table 5.1: Input data used in the simulations Parameter R1 R2 µ0 zvisc Load E ν τ0 ρ0 λ A

Value 0.01 2.0 0.14 0.6 105 206 0.3 3.3 870 0.1 0.2

Dimension m m Pa s N m−1 GPa MPa kg m −3 mm µm

chosen in order to enhance the non-Newtonian behavior of the lubricant. The velocity of the ridge uridge equals 1.1 m/s for s = +5/3 and 0.1 m/s for s = −5/3.

5.3.1 Interpolation of solution data The numerical set of data representing the solutions obtained from the different approaches is not defined on exactly the same grid points nor at exactly the same points in time. In order to compare the solutions it is necessary to evaluate them on the same grid points. This was done using linear interpolation of the Reynolds solutions on to the CFD solution grid points, both spatially and in time. This is actually a subject of discussion since the form of interpolation used, e.g., linear, second order, third order, spline interpolation of different orders, etc., affects the results. Linear interpolation was chosen since the others produced too large negative pressures in zones of large pressure gradients. These pressures were of minor significance for the comparison, since they were placed in the inlet and the outlet zone of the EHD lubricated conjunction. However, the ‘boundary zones’ of the traveling ridge does also show large pressure gradients. Thus, the interpolated data by itself induce deviations between the two approaches. The possible true deviations because of the different mathematical models might, because of the interpolation of solution data, become partly concealed.

44


5.3.2 Error estimation To determine possible deviations between the two approaches, due to the different mathematical models, the confidence intervals of errors for each approach must be well known. Conclusions can be only drawn about possible deviations if these intervals do not overlap. Therefore, an extensive error analysis was carried out to validate the results. The strategy used was to check the maximum error in the whole domain (both spatially and in time). For the Reynolds approach, the defined error functions are two-dimensional. In addition, the CFD-approach has one more dimension, i.e. the error is measured across the film as well (z-direction). All errors measured were obtained by halving the step size. The time error measure is defined as: ∆t (5.5) εtime (∆t) = max|Φi, j,2n − Φi, j,n (∆t) |. i, j,n 2 That is, the solution profile Φ i, j,k , which in this comparison is pressure profile, is considered as a function of the step size (∆t) in time where i, j and k are indexes for x, z and t respectively. However, this error will of course include the error in the spatial domain, e.g., the finer the spatial step size (∆x and/or ∆z) the smaller the error ε time . The same principle was adopted to define the spatial errors by considering the solution profiles as functions of the spatial step size. ∆x 1 (5.6) εx (∆x) = max|Φ2i, j,n − Φi, j,n (∆x) |. 3 i, j,n 2 ∆z 1 εz (∆z) = max|Φi,2 j,n (5.7) − Φi, j,n (∆z) |. 3 i, j,n 2 These errors (Eq. (5.6) and Eq. (5.7)), include the time errors in the same way as Eq. (5.5) includes the spatial errors. In addition, both the time and spatial error measures include a convergence error since the solution processes in both the Reynolds and CFD approach are not numerically exact. In the Reynolds approach, the convergence error is the maximum residual value of the Reynolds equation at each time step. This is sufficient since the residual value, and thus the errors of the film thickness equations, are small in comparison. The CFD approach controls convergence by summarizing the errors of each control volume. It was found that, for both the approaches, a clear relationship between the step size in time (∆t) and the spatial step size (∆x) exists. That is, for the running conditions specified in Table 5.1, u s ∆t/2 must be less than or equal to the step size ∆x. Applying this condition yields the appropriate behavior of the time error Eq. (5.5), e.g., εtime ∆t2 1 (5.8) εtime (∆t) 2 However, in the case of a slowly moving ridge u ridge ≤ u/4 the transient effects are well resolved even with values of u s ∆t/2 larger than ∆x. The time steps and meshes for the three different cases with input data given by Table 5.1 are shown in Table 5.2 whilst Table 5.3 summarizes the discretization errors for the three cases. The error denoted by ε spatial equals the maximum of the errors ε x and εz . Case 1 is the stationary solution for s = ±5/3 used as initial guess for the transient Case 2 and Case 3. The largest time errors are, as expected, the ones for Case 3 because of the enhanced transient effects due to the fast moving ridge.

5.3. THE MODEL PROBLEM

45

Table 5.2: Meshes and time steps used in the simulations Parameter Mesh CFD Mesh Reynolds ∆x CFD ∆x Reynolds ∆t CFD ∆t Reynolds

Case 1 400x20x1 1024x0x1 1.5 0.568 -

Case 2 400x20x1 1024x0x1 1.5 0.568 60.0 0.948

Case 3 400x20x1 1024x0x1 1.5 0.568 1.8 0.474

Dimension µm µm µs µs

In terms of equations, the higher velocity of the ridge enhance the influence of the squeeze term in Reynolds equation, viz. ∂ρ ∂h ∂ (ρh) = h +ρ . ∂t ∂t ∂t

(5.9)

The second term in Eq. (5.9) may be expanded to yield: ρ

∂h ∂ ∂ψ = ρ (h0 + δ) + ρ , ∂t ∂t ∂t

(5.10)

If the second term in Eq. (5.10) is expanded, the expression given by Eq. (5.11) for the influence of the velocity of the ridge on the squeeze effect is obtained: −20 ln(10) u ridge ∂ψ x − xc (t) ρ = ψ (x,t) − ∂t λ λ (5.11) 2 x−x (t) c 2π uridge 2π (x − xc (t)) −10 λ sin A 10 , λ λ since xc (t) = xs − uridge t. From Eq. (5.11) it is clear that the velocity of the ridge linearly enhance squeeze effects. In order to resolve these enhanced transients numerically, a smaller step size in time (∆t) is needed. In Case 2, (u ridge ∆t/2)/∆x 0.08 and (uridge ∆t/2)/∆x = 2.0 for the Reynolds solution and the CFD solution respectively. In Case 3, (uridge ∆t/2)/∆x 0.46 and (u ridge ∆t/2)/∆x = 0.66 for the Reynolds solution and the CFD solution respectively. As indicated by Table 5.3, the error will generally be larger in the case of positive slide-to-roll ratio, even if the relation (u ridge ∆t/2)/∆x 1 applies. Table 5.3: Discretization errors (%) Parameter s εtime CFD εtime Reynolds εspatial CFD εspatial Reynolds

Case 1 ±5/3 0.34 0.021

Case 2 −5/3 0.49 0.057 1.1 0.12

Case 3 +5/3 3.61 1.04 1.1 0.17


46

5.4 The results of the comparison

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2

0

-0.2

-0.1 x (mm)

0

0.1

P (GPa)

H (µm)

In this section, results from the numerical simulations are given. In Figures 5.1, 5.2 and 5.3 corresponding to Case 1, 2 and 3 respectively, the CFD solution is denoted by small circles and the Reynolds equation solution by solid lines. The deviations in the case of a perfectly smooth contact, visualized in Fig. 5.1, were extremely small. In this case the sign of the slide-to-roll ratio does not affect the result because of the isothermal assumption. It should also be noted that the pressure spike is not present and that the exit constriction has a smooth geometrical curvature due the strong non-Newtonian behavior of the lubricant. In Fig. 5.2, Case 2 is visualized by some

0

Figure 5.1: Pressure- and film thickness- profiles for a smooth contact, s = ±5/3. CFD solution (o), Reynolds (-). intermediate time steps of the passage of the slower moving ridge (u ridge = 0.1 m/s, s = −5/3) through the EHD lubricated conjunction. Once again the deviations are very small. The vertical line shown in the figure indicates the center of the undeformed ridge. Under the running conditions presented in Table 5.1, no phase difference between the ridge and the film thickness distortion can be observed, see for example Venner and Lubrecht [11]. This is attributed to the moderate pressure and the strong non-Newtonian behavior. Finally, Fig. 5.3 also shows some intermediate time steps for the passage of a ridge through the EHD lubricated conjunction. The ridge is attached to the faster moving surface and moving with velocity 1.1 m/s. The speed of the slow moving surface is again 0.1 m/s. The maximum deviation between the two approaches is approximately 1%. This is, however, smaller than the accumulated errors in the numerical solutions. So, if finer time and/or spatial resolutions are used together with stronger demands of the convergence error, this deviation will probably vanish.

5.5 Discussion and concluding remarks The numerical experiments performed in this comparison show that the deviations between the two approaches different approaches are small. One reason is that the thin film approximation is accurate because the shape of the traveling ridge have a small film thickness to

1 0.8

1.2

0.6

0.8

0.4

0.4

0.2 0

2

1 0.8

1.2

0.6

0.8

0.4

0.4

0.2 -0.2

-0.1

0

0.1

0

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2

0

-0.2

-0.1

0

0.1

0

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2

0

-0.2

-0.1

0

0.1

0

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2

0

-0.2

-0.1 x (mm)

0

0.1

P (GPa)

0.1

P (GPa)

0

P (GPa)

H (µm) H (µm)

-0.1

1.6

0

H (µm)

-0.2

P (GPa)

2 1.6

0

H (µm)

47

P (GPa)

H (µm)

5.5. DISCUSSION AND CONCLUDING REMARKS

0

Figure 5.2: Pressure and film thickness profiles at five intermediate time steps when the ridge is attached to the slower moving surface, u ridge = 0.1 m/s, s = −5/3. CFD solution (o), Reynolds (-).

wavelength ratio. Larger deviations would be expected if the ratio increased. A Newtonian fluid would probably also increase the deviations. The choice of running conditions might also be a reason for the good agreement between the two approaches. From a CFD viewpoint, the results are encouraging because of the possibilities of extending the simulations of the EHD lubricated conjunction to include more of the inlet and outlet zones. If the numerically complicated flow problem inside the conjunction can be simulated, it should also be possible to simulate the full inlet and outlet flow. Such

1 0.8

1.2

0.6

0.8

0.4

0.4

0.2 0 1 0.8

1.2

0.6

0.8

0.4

0.4

0.2 -0.2

-0.1

0

0.1

0

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2 -0.2

-0.1

0

0.1

0

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2

0

-0.2

-0.1

0

0.1

0

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2

0

-0.2

-0.1 x (mm)

0

0.1

P (GPa)

0.1

P (GPa)

H (µm)

0

2

0

H (µm)

-0.1

1.6

0

H (µm)

-0.2

P (GPa)

H (µm)

0

P (GPa)

H (µm)

2 1.6

P (GPa)


48

0

Figure 5.3: Pressure and film thickness profiles at five intermediate time steps when the ridge is attached to the slower moving surface, u ridge = 1.1 m/s, s = +5/3. CFD solution (o), Reynolds (-).

simulations would add more clarity to the reformation process, which is a very important topic of investigation. CFD would also be a suitable tool for simulating starved contacts by means of two-phase flows. Another use of CFD would be for simulating wear debris traveling through the conjunction (particle transport simulation). The comparison also enhances the reliability of the Reynolds approach which, for the running conditions chosen, including large non-Newtonian effects predicts approximately the same pressure and film thicknesses as the CFD solution. In particular, this comparison

5.5. DISCUSSION AND CONCLUDING REMARKS

49

has also shown that it is possible to use the Block-Jacobi relaxation method as a smoother for the multilevel FAS algorithm. Moreover, the Reynolds approach described in chapter 2 has been verified. The choice of step size, both spatially and in time, is crucial when it comes to numerical convergence in both approaches. The passage of a single ridge with small amplitude to wavelength ratio has been simulated. If the wavelength is kept constant and the amplitude is enlarged it would almost certainly be necessary to consider the speed of the asperity instead of the surface mean speed when choosing the step size in time, e.g., u ridge ∆t ≤ ∆x. On top of this, asperities with large amplitude to wavelength ratios have to be well resolved, meaning small spatial step sizes and thus even smaller step sizes in time are required. This also applies to situations where one of the surfaces can be considered as perfectly smooth. When dealing with real applications the roughness of both surfaces are often of comparable size. Combined with sliding, this leads to a continuously changing effective surface roughness inside the EHD lubricated conjunction. In such a case, the squeeze effects are underestimated if the surface roughness is imposed on one of the contacting surfaces. Numerically, using either of the two approaches, it is necessary to consider the amplitude to wavelength ratio as being time dependent and hence the spatial step size should be chosen so that the largest amplitude to wavelength ratio in time is resolved. The rheological model used in this work is the Eyring model which has a reduced shear stress to shear rate dependence when the pressure and shear rates are high. In order to enable CFD simulations at high contact pressures, some form of ‘damped’ relationship between shear stress and shear rate is necessary in order to avoid a singularity in the momentum equation Eq. (5.1). The numerical simulations performed in this comparison are isothermal. It is, however, necessary to include thermal computations in order to have an appropriate fluid flow model. Such a project is ongoing and will hopefully add more understanding to transient EHD lubrication. Based on the comparison the following conclusions can be drawn: • As far as the numerical results are considered, deviations between the two approaches are small. It must be emphasized, however, that the amplitude to wavelength ratios used in this work are small. • The numerics of both the approaches seem stable. However, the case when the ridge is attached to the faster moving surface is somewhat more difficult to simulate due to the greater influence of transient effects. • The deviation in pressure between the two approaches is approximately 10-20 times higher than the deviations in film thickness. This is also the case for numerical errors measured for each approach separately.

Chapter 6

Simulations of Rough FL In the case of two-sided roughness, the hydrodynamic action caused by asperities passing each other inside the conjunction gives rise to local effects that may lead to plastic deformation and/or collapse of the lubricant film. For an elliptic contact, there are numerous possible situations which would give rise to such local effects. For example, a ridge with a certain lay passing an other ridge with a different lay. However, for a 2D line contact there are precisely four different ‘asperity overtaking’ situations, shown in Fig. 6.1, assuming that the infinitely wide ridges and dents can be approximated by a single wave with certain wavelength and amplitude. The dents/ridges shown in Fig. 6.1 are mathematically described by Eq. (5.4).

6.1 The different overtaking situations A number of transient simulations were carried out to investigate the influence of the elementary surface features, shown in Fig. 6.1, on both the film thickness and the pressure. The theory describing the mathematical model and the solution procedure is found chapter 2 and it should be pointed out that the effects of interest are enhanced by assuming a non-Newtonian fluid behavior. The focus is on the film thinning effect caused by a dent-ridge overtaking event. The other three possible overtaking situations have also been studied briefly and are described below. In the ‘Ridge-Ridge’ case an infinitely wide ridge attached to the faster moving surface, with speed u2 , passes an identical ridge attached to the slower moving surface, with speed u1 , inside the EHD lubricated conjunction. This is the case where the highest pressure and the thinnest film are to be expected. In Fig. 6.2 the dimensionless film thickness H is plotted as a function of both dimensionless time T and dimensionless space X. The figure shows the overtaking situation and it is clear in that the global minimum occurs as a result of this. In the ‘Ridge-Dent’ situation, with the height of the ridge equal to the depth of the dent and s = 4/3, shown in Fig. 6.3, the global minimum film thickness occurs in the outlet region as the faster moving ridge exits the conjunction. That is, it is not the result of the overtaking process. The possibility of internal cavitation is highest in the ‘Dent-Dent’ case which for s = 4/3 is shown in Fig. 6.4, as the overtaking creates a relatively low viscosity region which 51

CHAPTER 6. SIMULATIONS OF ROUGH FL

52

Ridge-Ridge

u1

u2 Ridge-Dent

u1

u2 Dent-Dent

u1

u2 Dent-Ridge

u1

u2

Figure 6.1: The four different asperity overtaking situations in an EHD lubricated line contact

H (X,T )

0.8 0.7

0.5

u1 0.6

1 -1

0.5

u2

-0.5 0.4

0 0.5 X

1 11

13

12

14

15

0.3

T

Figure 6.2: Dimensionless film thickness in the Ridge-Ridge overtaking situation with s = 4/3. causes a relaxation of the elastic deformation. Moreover, at the end of the event the contact geometry forms a rapidly convergent gap that produce a large pressure spike. The pressure

H (X,T )

6.1. THE DIFFERENT OVERTAKING SITUATIONS

53

0.7 0.5

u1 0.6

1 -1

0.5

u2

-0.5 0.4

0 0.5 X

1 11

13

12

14

15

0.3

T

Figure 6.3: Dimensionless film thickness in the Ridge-Dent overtaking situation with s = 4/3. spike in this case is of comparable magnitude to the pressure spike in the Ridge-Ridge event with ridges of half the amplitude. There is a local minimum film thickness formed due to the overtaking, in connection with the pressure build up, but the global minimum occurs as the faster moving dent exits the conjunction.

H (X,T )

0.9 0.8

0.5

u1 0.7

1 -1

0.6

u2

-0.5 0.5

0 0.5 X

1 11

13

12

14

15

0.4

T

Figure 6.4: Dimensionless film thickness in the Dent-Dent overtaking situation with s = 4/3. The ‘Dent-Ridge’ case, like the Ridge-Dent case, would seem intuitively to be quite harmless. However, for certain choices of the point of overtaking this case has a significant film thinning effect. This effect is shown in Fig. 6.5 as the global minimum film thickness occurs in the beginning of the overtaking process. This is because, when the dent enters the contact zone the lubricant it holds will have lower viscosity than the fluid downstream has and in the beginning of the overtaking process this low viscosity fluid inside the dent will reduce the pressure and thus relax the elastic deformation upstream the ridge. In doing so, the ridge will regain some of its undeformed height and expands toward the opposite


54

H (X,T )

0.9 0.5

0.8

u1

0.7

1 -1

0.6

u2

-0.5 0.5

0 0.5 X

1 11

13

12

14

15

0.4

T

Figure 6.5: Dimensionless film thickness in the Dent-Dent overtaking situation with s = 4/3. surface. Because of this film thinning effect the case of a Dent-Ridge overtaking can be seen as a potential breakdown mechanism which is the reason for the interest in it in this work.

6.2 The Dent-Ridge overtaking For the simulation of the Dent-Ridge overtaking situation, five different slide-to-roll ratios (s) were initially selected; 1/3, 2/3, 3/3 , 4/3 5/3. The amplitude A and the wavelength Λ, in Eq. (5.4), were modeled as A = 2 · 10 −7 li and Λ = 10−4 mi (in dimensions), where l and m are integers. A positive l i , where i is the surface index, sets a ridge amplitude, whereas a negative l i sets the depth of the dent. Fig. 6.6 shows intermediate time steps of the Dent-Ridge overtaking event for dent parameters l2 = −5, m2 = 1, ridge parameters l 1 = 1, m1 = 1 and s = 4/3. It is clear that the Dent-Ridge overtaking causes a thinner film than the reference simulation of a single ridge, i.e., a ridge attached to the slower moving surface with parameters l 1 = 1 and m1 = 1 contacting a perfectly smooth surface, passing the conjunction. This is due to the relaxation of the elastic deformation of the ridge at the beginning of the event. In order to evaluate the influence of different parameters, i.e., s, w, (l 1 , m1 ) and (l2 , m2 ), on the film thinning effect occurring as a result of the Dent-Ridge overtaking a measure was defined by: max

ξ (l, m) =

−1.3≤X≤0.0

|H l,m (X, T ) − H re f (X, T ) | min|H l,m (X, 0) |

(6.1)

X

Note that this measure assumes a large value when the actual film thickness is small. Simulations of a ridge passage only for each of the slide-to-roll ratios chosen were used as reference for the Dent-Ridge simulations. The exact position of the overtaking was set as −0.3253 · b in order to have a suitable ‘overtaking window’ for all slide-to-roll

6.2. THE DENT-RIDGE OVERTAKING

55 H Hre f P

1.5 1 0.5 0 -2

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5

-1

-0.5

0

0.5

1

1.5

1.5 1 0.5 0 -2 1.5 1 0.5 0 -2 1.5 1 0.5 0 -2

X

Figure 6.6: Dimensionless film thickness (red), reference film thickness (blue) and pressure (black) at four different time steps showing the dent caused thinning effect for s = 4/3. ratios tested. The ridge parameters for the reference film thickness were held constant at l1 = 1 and m1 = 1. The set of simulations where both dent amplitude (l 2 ) and dent wavelength (m 2 ) were varied are shown as a contour map of the measure ξ in Fig. 6.7. The dent parameters were varied to give a total of forty-two combinations in the ranges l 2 = −1, −2, . . ., −6 and mi = 1, 2, . . . , 7. The figure shows that the larger the amplitude of the dent, the larger the film thinning effect and thus the thinner the film thickness. As far as the dent wavelength is concerned, an optimal and worst case wavelength in the range of wavelengths simulated becomes clear from Fig. 6.7. This is interesting since in many tribological applications dents are thought of as lubricant depots and these simulations show that there is an optimal dent wavelength to maintain the greatest film thickness possible, at least when considering an EHD lubricated line contact. When studying the different slide-to-roll ratios, the dent amplitude (l 2 ) was varied as in Fig. 6.7 but here the dent wavelength (m 2 ) as well as the ridge parameters were held constant, i.e., m2 = 1, l1 = 1 and m1 = 1. As Fig. 6.8 shows, the amount of film thinning is clearly dependent of the slide-to-roll ratio and that the thinning increases with increasing s. The film thinning is both greater in magnitude and increase more rapidly with l 2 for


56 -1

0.45 0.4

-2

0.35 -3 l2

0.3 0.25

-4

0.2 -5 0.15 -6 1

2

3

4 m2

5

7

6

Figure 6.7: Contour map of ξ for s = 5/3, y-axis corresponds to dent amplitude, x-axis to dent wavelength. l 1 = 1 and m1 = 1. 0.5

0.4

s = 5/3 s = 4/3 s=1 s = 2/3 s = 1/3

ξ

0.3

0.2

0.1

0 -1

-2

-3

-4

-5

-6

l2

Figure 6.8: ξ as a function of the slide-to-roll ratio (s). x-axis corresponds to dent amplitude. higher slide-to-roll ratios. The overall maximum of ξ is approximately 0.4, i.e., the dent causes a film thickness decrease that corresponds to 40% of the smooth surface minimum film thickness value. In the case of s = 5/3, the effect of ridge amplitude, i.e., the cases of l 1 = 1 and 2, were also studied. In both cases m 2 was varied in the range 1, 2, . . . , 7 and the dent parameter l2 were taken as −l1 . The results are shown in Fig. 6.9. These results confirms that the measure ξ assumes both a local maximum and a local minimum. That is, a dent wavelength (m2 ) somewhere between 1 and 3 yields the thinnest film and a value somewhere between 4 and 6 the thickest film. Fig. 6.9 also shows that the film thinning effect is greater in the case of doubled ridge amplitude (l 1 = 2). Doubled ridge amplitude does not increase the film thinning effect more than a few percent in the parameter range chosen. However, the difference seems to increase with increasing wavelength. Together with the results

6.2. THE DENT-RIDGE OVERTAKING 0.3

57

l1 = 1 l1 = 2

0.28

ξ

0.26 0.24 0.22 0.2 0.18 1

2

3

4 m2

5

6

7

Figure 6.9: ξ as a function of ridge amplitude for s = 5/3. x-axis corresponds to dent wavelength. shown in Fig. 6.7, this indicates a tendency toward film breakdown if one or a combination of the parameters, ridge amplitude (l 1 ), dent amplitude (l 2 ) and dent wavelength (m 2 ) are combined in an unfortunate way. In real applications, such combinations will occur frequently. In order to investigate the effect of load (w), another set of simulations with s = 5/3 using the same set of dent and ridge parameters as those in the parameter study of the slide-to-roll ratio, shown in Fig. 6.8, were carried out. The results, shown in Fig. 6.10, indicate that a doubling of the load gives rise to thinner films. It can also be seen that the minimum film thickness corresponds to approximately 45% of the smooth film thickness and that the slope of ξ is greater in the case of the doubled load. 0.5

w = 10 kN w = 20 kN

ξ

0.4

0.3

0.2

0.1 -1

-2

-3

-4

-5

-6

l2

Figure 6.10: ξ as a function of load for s = 5/3. x-axis corresponds to dent amplitude. In all these simulations, several significant input parameters were held constant, some


58

of these are shown in Table 6.1. It is not possible to draw general conclusions from the Table 6.1: Fixed parameters Parameter R1 R2 E1 E2 ν1 ν2 µ0 α p0

Value 0.01 2.0 206 · 109 206 · 109 0.3 0.3 0.14 2.10 · 10−8 1.98 · 108

Dimension m m Pa Pa Pa s 1/Pa Pa

investigation presented here, but a trend toward film thinning when an asperity is free to recover its shape in the valleys or dents between asperities on the opposing surface. Clearly, at a certain combination of lubricant and roughness parameters film breakdown will occur. Such a breakdown will be due to the low viscosity regions, leading to local low pressure or even internal cavitation, in the valleys or dents rather than to very high contact pressures as might intuitively be thought to be the cause of film breakdown. In the cases presented here, no actual film breakdown was observed. However, in a line contact of finite length the film thinning will be greater. This is due to side leakage from the valleys which will lead to even lower pressures allowing asperities to expand more rapidly and to a greater extent. This will eventually lead to film breakdown. Further investigations into internal cavitation phenomena and side leakage effects are therefore needed to obtain knowledge about the film thinning effects studied here. It is, of course, also necessary to introduce more realistic surface roughness features. The smaller the wavelengths become, the smaller the amount of sliding is required to cause this type of local film thinning. However, as indicated in Almqvist and Larsson [36], the validity of a Reynolds approach becomes questionable as the ratio of asperity amplitude to asperity wavelength increase. Therefore, the full Navier-Stokes equations may well form the theoretical basis for future research in the area of EHL lubrication. Using more realistic surface topographies will also increase the need for surface characterization when evaluating results and in deriving tools more directly applicable for engineering purposes. The technique described in chapter 4 might in this sense become an important tool.

6.3 Conclusions Based on the asperity overtaking investigation performed the following conclusions can be drawn: • Certain effects in a rolling/sliding EHD lubricated contact cannot be studied without using a two-sided roughness treatment. • A Dent-Dent overtaking event gives rise to a pressure spike equivalent to a RidgeRidge overtaking event where the ridges have half the amplitude of the dent. Bearing

6.3. CONCLUSIONS

59

capacity is therefore not significantly affected by this type of overtaking event. • A Ridge-Ridge overtaking event produces the thinnest films and largest pressure spikes compared to the other three possible overtaking situations for a given ridge amplitude. However, film breakdown could not be simulated in this work because of the model used. • Local pressure effects resulting from a ridge-dent overtaking event are of comparable magnitude to the Dent-Dent case. Because of the involvement of a ridge, the film becomes slight thinner causing this case which implies less bearing capacity. • In a sliding contact, a dent overtaking a ridge can cause film breakdown. This might be the result of a too large a ridge, a too wide or deep dent, too high a load, or a combination of these parameters. • There is a dent wavelength in a Dent-Ridge overtaking event that minimize the effects of film thinning and thus optimizes the film forming capability. • As the slide-to-roll ratio increases, the film thickness decreases as a result of the Dent-Ridge overtaking event. • If the depth of the dent equals the height of the ridge, then a larger ridge implies a thinner film. Doubling the height , however, does not have a devastating effect. • The contact load influences the film thinning effect. Load related effects are somewhat larger than those for increasing ridge amplitude.

Chapter 7

Simulations of the dry contact Two-dimensional profiles from four different steel surfaces (Lord and Larsson [31]) were studied, each originating from different machining processes. These were hardened and for the numerical simulations their hardness is set to H s = 4.0 GPa. The modulus of elasticity (E) and the Poisson ratio (ν) were assumed to be 210 GPa and 0.3 respectively for all four surfaces. The surfaces were all initially unidirectional ground and the profiles are taken perpendicular to this direction. However, as Table 7.1 shows, their 2D profiles differ because of additional machining. Due to these different additional machining processes, the surfaces will generate different contact conditions. These contact conditions are here represented by various statistical surface parameters, e.g., R a , Rz , Rsk , Rk , but also the percentage of real area of contact A r /An and the percentage of plastically deformed area A p /An . A parameter study varying the load was preformed and the contact conditions described above the outcome of the study. Care was taken to ensure that the metrological limitations of the optical topometer, used to measure the topographies, do not significantly influence the resulting contact conditions. That is, because of the size of the measured area, some low frequency content will be lost and because of limitations in resolving the topography, some high frequency content will be lost. The surface characterization technique proposed in Chapter 4 was used. The investigation showed that the measurement resolution was sufficient and the accuracy of the outcome from a parameter study varying the load could be verified. Table 7.1: Data for the four undeformed surface profiles. # A B C D

Machining Unidirectional ground + Phosphated and dephosphated Unidirectional ground Unidirectional ground + 3 x chemically deburred Unidirectional ground + shot-peened

R a (µm) 0.93

Rz (µm) 5.50

Rsk -1.00

Rk 4.40

0.30 0.14

1.69 0.78

-0.25 -0.02

3.20 3.50

0.35

2.22

-0.07

3.10

To enrich the data given in Table 7.1 the Abbott-Firestone bearing ratio curve and the 61

CHAPTER 7. SIMULATIONS OF THE DRY CONTACT

62

probability density distribution (height distribution) of the four different topographies are visualized in Fig. 7.1. It is clear, both from Table 7.1 and Fig. 7.1, that surface profile A is much rougher than the others. This surface was phosphated and then dephosphated. The surface will inherit some special properties from this, but in the numerical simulations the hardness of this surface is set to 4.0 GPa as for the other surfaces. Surfaces B and D differs in that D was shot-peened. This can be seen in the table as increased R-parameter values for surface profile D, and in the figure both as an overall increase but also as a small peak on the profiles at the upper part of the surface. Surface C differs from B because it has been chemically deburred; a process which has smoothed the profile significantly. There is a 7

7

A B C D

6

5 z (µ m)

z (µ m)

5 4 3

4 3

2

2

1

1

0 0

20

40

60

80

A B C D

6

100

%

0

P (z)

Figure 7.1: Abbot-Firestone bearing ratio (left) and the height distribution (right). close relationship between the Abbott-Firestone bearing ratio and the height distribution P (z), i.e., P (z) can be obtained by differentiating the Abbot-Firestone bearing ratio with respect to z. Because of this form of relationship P (z) shows an amplification of certain properties of the Abbot-Firestone bearing ratio which makes it a good complement for the characterization of surfaces.

7.1 Varying the applied load The real area of contact is an important quantity in terms of surface performance under different running conditions. A small real area of contact may reduce friction but could lead to unwanted plastic deformation. The percentage of real area of contact A r /An for the four surfaces is plotted as a function of load in the left of Fig. 7.2. It can be seen that the real area of contact increases almost linearly with load. On the right of Fig. 7.2 the plasticity index, here defined as A p /Ar , is plotted for the four surfaces. It should be noted that the relationship between A p /Ar and applied load is quite linear in the high load region. It is also clear that surface profile A is subjected to the largest amount of plastic deformation, which might be explained by its R-parameter values. Even though the R-parameter values of surface profile A are significantly larger than those for the other surface profiles they do not differ significantly in terms of the percentage of real area of contact. For surface profile C, the plasticity index is low and the real area of contact high which is also suggested by its comparably small R-parameter values. Surface profiles B and D have quite similar Rparameter values, and the higher A p /Ar ratio for D is probably explained by its slightly lager R-parameter values. In Fig. 7.3, the R a - and the Rz -values for the four surfaces are

7.1. VARYING THE APPLIED LOAD

63

10

0.8 A p /Ar (-)

Ar /An (%)

8

1 A B C D

6 4 2 0 0

0.6 A B C D

0.4 0.2

20

40 60 Load (kN)

80

0 0

100

20

40 60 Load (kN)

80

100

Figure 7.2: Percent real area of contact (left) and plasticity index (right). plotted as functions of load in the left and right sub-figure. The variation with load is small and the decrease is linear. In terms of magnitude of the R-parameters, Table 7.1 suggests 6

1

5

A B C D

0.6

Rz (µ m)

Ra (µ m)

0.8

0.4

3 2

0.2 0 0

A B C D

4

1 20

40 60 Load (kN)

80

0 0

100

20

40 60 Load (kN)

80

100

Figure 7.3: R a (left) and Rz (right). a ranking of the surface profiles, i.e.: A, D, B, C. From Fig. 7.3 it is clear that this ranking is preserved when increasing the load. In the same way, Table 7.1 suggests a ranking in terms of the skewness R sk . As shown to the left of Fig. 7.4. this ranking is also preserved when increasing the load. However, 0

5

A B C D

-0.2

4.5 Rk (-)

Rsk (-)

-0.4 -0.6 -0.8

A B C D

4 3.5

-1 3

-1.2 -1.4 0

20

40 60 Load (kN)

80

100

2.5 0

20

40 60 Load (kN)

80

100

Figure 7.4: Skewness (left) and kurtosis (right). the Rsk -parameter values seem not to posses the strong linear behavior seen for the Rparameters, except for surface profile A. The kurtosis, R k , is visualized to the right in

CHAPTER 7. SIMULATIONS OF THE DRY CONTACT

64

Fig. 7.4. Considering the ranking suggested by Table 7.1, for the kurtosis, a variation can be observed as the R k -parameter values for C rapidly decrease in the low load region when increasing the applied load. From this figure it can also be concluded that load alters the kurtosis of C the most and D the least.

7.2 Conclusions Results from the numerical method, described in Chapter 3, for the 2D contact between linear elastic perfectly plastic rough surfaces is presented. The model takes into account the energy dissipation due to plastic deformation. Based on simulations of four different surface profiles originating from measured topographies the following conclusions can be drawn: • The four surface profiles show similar behavior in terms of the real area of contact, plasticity index and R-parameters, with increasing the applied load. • If the R-parameter values are used to characterize the roughness of the four examined surface profiles, the real area of contact can be considered almost independent of the surface roughness. However, as the load increase the deviations in the real area of contact increase and surface profile C, which in terms of its R-parameters is the smoothest, attains the largest real area of contact. Profile A which in terms of its R-parameters is the roughest surface, does not show a significant difference compared to the profiles B and D. • The plasticity index is highly dependent of surface roughness at least if the Rparameter values are being used to characterize the roughness. Here the ranking of R-parameter made up by Table 7.1 is preserved. It should be noticed that these results are in part also due to the assumption that all surfaces have the same hardness, i.e., 4.0 GPa.

Chapter 8

Concluding Remarks The Reynolds based numerical tool developed here allows simulation of the transient effects occurring in line contacts that operates in the FL regime. Until now the tool has been used for simulation of elementary surface features passing each other inside the lubricated conjunction. The purpose of this work was to investigate whether any local effects that may arise may lead to a transition from the FL to the ML regime. This was done and these local effects were found to significantly affect the film formation. The contact mechanics based tool has been implemented and tested. The accuracy of the results predicted by simulating the dry contact between surface profiles originating from measured topographies were verified by the surface characterization technique also developed.

65

Chapter 9

Future Work The long term goal is to develop engineering tools, to predict different contact conditions, that takes into account the surface topography. Suggested future work is therefore to perform systematic investigations of rough surfaces, in both the lubricated an unlubricated situations, using the developed numerical tools. It should also be noted that the developed tools can not handle the real contact case, but simplified forms. Improvements of the numerical tools are therefore needed, for example to account for other physical properties of the contact between rough surfaces as thermodynamical effects and effects due to lubricant and surface chemistry.

67

68

CHAPTER 9. FUTURE WORK

Part II

Appended Papers

69

Paper A

71

73 Tribology International 37 (2004) 61–69

A comparison between computational fluid dynamic and Reynolds approaches for simulating transient EHL line contacts T. Almqvist 1 , A. Almqvist, R. Larsson Luleå University of Technology, Department of Mechanical Engineering, Division of Machine Elements, Luleå, SE-971 87 Sweden Received 9 April 2003; received in revised form 12 June 2003; accepted 20 June 2003

Abstract When simulating elastohydrodynamic lubrication (EHL), the Reynolds equation is the predominating partial differential equation for prediction of the fluid flow. Also very few attempts have been carried out using the full momentum and continuity equations separately. The aim of this investigation is to compare two different approaches for simulation of EHL line contacts where a single ridge travels through an EHL conjunction. One of the approaches is based on the Reynolds equation, addressing the coupling between the pressure and the film thickness. The solver uses the advantages of multilevel techniques to speed up the convergence rate. The other approach is based on commercial CFD software. The software uses the momentum and continuity equations in their basic form, enabling numerical simulations outside the contact regions, as well as in the thin film region to be carried out. The numerical experiments show that, under the running conditions chosen, only small deviations between the two approaches can be observed. The results are encouraging from several viewpoints: verification of the codes, the possibilities of further developments of the CFD approach and the justification of using a Reynolds approach under the running conditions chosen. Keywords: Elastohydrodynamic lubrication, Surface roughness, Fluid mechanics

1 Corresponding author. E-mail address: [email protected]

PAPER A.

74

A.1

Introduction

Theoretical investigations in the field of elastohydrodynamically lubricated (EHL) contacts are almost universally based on the Reynolds equation, derived more than 100 years ago by Osborne Reynolds [37]. The equation is derived from the Navier-Stokes equations and only a few attempts have so far been carried out to solve EHL problems numerically using the momentum and continuity equations as a basis. Almqvist and Larsson [5, 32] used the Navier-Stokes equations to simulate steady state EHL line contacts under both isothermal and thermal conditions. Schäfer et al. [38] used the Stokes equations to simulate isothermal EHL line contacts and investigated the pressure variation across the lubricant film. The use of CFD (computational fluid dynamics) software when simulating EHL line contacts certainly has advantages compared with the Reynolds equation based numerical tools. The assumption of small inertia inside the fluid and the thin film approximation are not applied, which in some cases determine the accuracy of the predicted pressure and film thickness profiles. By not utilizing these assumptions about the fluid behavior, it is possible to extend the computational domain to include the flow outside the EHL conjunction. Additionally, the greater flexibility for modifying the rheology has advantages compared to a Reynolds approach, i.e. it is more straightforward to implement complex rheological models. A disadvantage of a CFD compared to a Reynolds approach is the reduction in efficiency, mainly due to the enlarged size of the problem. For 3D cases, four equations have to be solved in order to obtain the fluid pressure in comparison with the single Reynolds equation. Another disadvantage is that it is not straightforward to modify a CFD code to solve EHL problems. The validity of the Reynolds equation has recently been discussed by Odyck and Venner [39]. They found large differences between the Reynolds and Stokes equations in predicted load when the ratio ω is O (0.1). The subject of this paper is the investigation of two different approaches of simulating transient isothermal EHL line contacts numerically. That is, a Reynolds approach is compared to a CFD approach where the momentum and continuity equations are used to simulate the fluid flow. The Eyring rheological model is adopted in both these approaches and no elastic contributions are considered. The comparison is based on transient simulations of a relatively lightly loaded EHL line contact, 0.6 GPa maximum Hertzian pressure, where a single ridge is passing through the conjunction under sliding. A relatively large slide-toroll ratio was chosen in order to enhance the non-Newtonian behavior of the lubricant. An extensive error analysis is carried out to validate the results. That is, in order to determine the deviations between the two approaches, due to the different mathematical models, the confidence intervals of errors in each approach must be well known. If these intervals do not overlap, conclusions can be drawn about the different approaches. It will be shown that there is almost no difference between the solutions in the case of slide-to-roll ratio ±5/3 when smooth surfaces are considered. The transient effects caused by the ridge are very similar in both slide-to-roll ratio cases.

A.2

Governing equations

The governing equations for an isothermal fluid flow are those for momentum and continuity. When solving the fluid flow equations for moving grids, the equations in the CFD code are changed to account for volume changes in the transient terms by inclusion of the determinant of the metric tensor g, and the advective term to allow for a grid velocity

A.2. GOVERNING EQUATIONS

75

∂x/∂t, (CFX4.4 [34]). The equation of continuity reads: √ 1 ∂ ρ g ∂x + ∇ · ρu − =0 √ g ∂t ∂t

(A.1)

The density and time are here denoted by ρ and t, respectively, and u is the velocity field. The equation of momentum reads: √ 1 ∂ ρu g ∂x +∇· ρ u− ⊗ u = ∇ · σ. (A.2) √ g ∂t ∂t The symbol σ denotes the total stress tensor and the symbol ⊗ denotes the vector product ui × u j . The body forces are assumed to be vanishing small compared to the other terms, and have therefore been neglected. An Eyring rheological model is used for the dynamic part of the stress tensor. Here the total stress tensor attains the following form: 2µ ∇ · uδ + µ˙γ, 3

σ =

−pδ −

γ˙ =

∇u + (∇u)T

=

τ0 sinh

τ τ0

(A.3)

.

The pressure is denoted by p, and µ is the dynamic viscosity. The Eyring stress is denoted by τ0 and γ˙ˇ is the rate of deformation tensor. When modeling the fluid flow by the Reynolds equation, the form used by Conry et al. [2] is adopted here with the squeeze term incorporated: ∂ ρh3 ∂p us ∂ (ρh) ∂ (ρh) S − = 0. (A.4) − ∂x 12µ ∂x 2 ∂x ∂t The film thickness is denoted here by h. u s is the sum of the surface velocities u 1 and u2 . The non-Newtonian factor S is as according to Johnson and Tevaarwerk [10]:

µ2 (u2 − u1)2 Σ2 3 (Σ cosh Σ − sinhΣ) 1+ S = , 3 Σ sinh2 Σ (τ0 h)2 (A.5) h p . Σ = 2τ0 x The equations for film thickness, density and viscosity are common to the two approaches where the film thickness equation reads: h (x,t) = h0 (t) +

x2 + d (x,t) + ψ (x,t) . 2Rx

(A.6)

The reduced radius of curvature is here given by 1/R x = 1/Ru + 1/Rd and the symbol h0 (t) is an integration constant, which is determined by the force balance equation:

I

p (x,t) dx = w (t) .

(A.7)

PAPER A.

76

where I represents the whole solution domain. The elastic deformation is denoted d (x,t) and is defined as:

2 (1 − ν) ∞ d (x,t) = p x ,t ln |x − x|dx , (A.8) πE −∞ where Poissont’s ratio is denoted by ν and E is the modulus of elasticity. In the CFD approach, the elastic deformation equation is solved for the upper and lower boundaries separately. In the Reynolds approach, the deformations of the surrounding solids are composed into a single equation. The equation for the ridge ψ is defined as (see [40]): ψ (x,t) = A 10

−10

x−xc (t) 2 λ

2π (x − xc (t)) cos , λ

(A.9)

where A is the amplitude, x c (t) = xs + uridget is the position of the ridge and λ is the ridge wavelength. The equation for the viscosity is the Roelands equation, [41], and reads: p zvisc µ (p) = µ0 exp (ln(µ0 ) + 9.67) −1 + 1 + , (A.10) P0 where µ0 denotes the viscosity at an ambient pressure, P0 is a constant and zvisc is the pressure-viscosity index. The Dowson-Higginson density expression [35] reads: ρ (p) = ρ0

C1 + C2 p , C1 + p

(A.11)

where ρ0 denotes the density at ambient pressure and C 1 , C2 are constants.

A.2.1 Boundary conditions and cavitation treatment Specified pressure boundary conditions are used in the hydrodynamic equations at the flow boundaries in the CFD code, i.e. where the fluid enters or leaves the computational domain. The value of the pressure is set to atmospheric. When the pressure is specified, the boundary conditions used on the momentum equation are of Neumann type, or zero velocity gradients, CFX4.4 [34], since a Neumann boundary condition implies fully developed flow. In the Reynolds equation approach, the pressure at the boundaries is set to zero, but at the outlet the pressure gradient is also forced to zero by the cavitation condition P ≥ 0. In the CFD-approach this is not possible. A pressure correction must be calculated to force the mass flow to fulfill the continuity criterion. A way to achieve this, in the CFD approach is to modify the density in regions of negative pressure [5]. The advantage of this is that continuity is fulfilled and that internal cavitation are handled without any further modifications to the solution algorithm.

A.3

Numerics

This section explains the numerics of the two approaches. There are some similarities such as the order of the numerical schemes, first order in time and second in the spatial dimensions. However, the numerics differ in many ways because of the differences in

A.3. NUMERICS

77

theoretical models. The solution process for the Reynolds approach is based on a local relaxation method, e.g. it is based on the Block-Jacobi method. The CFD approach uses the Stonet’s method to converge the discretized system. The two approaches also differ in discretization methods. For example, the CFD approach uses the finite volume method whereas the Reynolds approach is based on finite differences.

A.3.1 The numerics for the CFD approach The CFD approach is based on commercial CFD software, CFX4.4 [34]. The software has the facilities of pre-processing, front-end, solver and post-processing modules. The code is based on a pre-processing module where geometry building and meshing is performed and a front-end module, where specifications of the fluid flow can be set through a command language. In the solver module, the discrete representation of the equations are solved and, finally, a post-processing module for the graphical analysis of the solution. In the CFD approach, the equations of momentum and continuity are solved numerically. The velocity field is obtained from the momentum equations, leaving the continuity equation for calculation of the pressure. Actually, the code combines the equations of momentum and continuity to obtain a pressure correction equation where the pressure correction is added to the velocity field so that continuity is enforced [42]. Therefore, if mass continuity is a demand in the solution domain the exit boundary condition used in the Reynolds equation is not appropriate. The solution algorithm used in this work is the ‘simple’ algorithm, CFX4.4 [34]. The code has a number of utility routines where it is possible to modify the code to account for more specific user requirements. To implement the expressions for effective viscosity and density, the user routines USRVIS and USRDEN are used. The elastic deformation in the contact is implemented through the user routine USRGRD and the velocity boundary conditions, through the routine USRBCS, see CFX4.4 [34]. The code uses a finite volume discretization of second order accuracy in the diffusive terms, and hybrid discretization for the advective terms (when the local Reynolds number exceeds two, the discretization for the advective terms are replaced by an upwind discretization). Because of the dominance of the diffusive terms in the contact region in EHL conjunctions, the scheme is assumed to be of second order accuracy. The order of the scheme is of importance when the discretization errors are calculated. The relaxation algorithms for the momentum and pressure correction equations are Stonet’s method. The CFD solver works in a segregated manner which means the equation for each variable has to be solved in turn, while the other variables are assumed to be constant. In the outer iterations, the matrix coefficients are updated. Such approach easily leads to instabilities and long computational times. The meshes used in the simulations were structured and uniform. The convergence criterion was set on the mass source residual (mass continuity error) and force balance.

A.3.2 The numerics for the Reynolds approach In the Reynolds approach used in this work, the Reynolds equation and the equation for the film thickness are coupled so that pressure and film thickness are solved simultaneously [1]. Multilevel techniques [11] have been used to accelerate the convergence of Reynolds approaches based on serial solution methods. For example, the pressure is the only unknown parameter. Here the multilevel technique is extended to accelerate the solution process of the coupled approach. The full multilevel algorithm is based on a specially

78

PAPER A.

developed iterative relaxation process. This relaxation process is a generalized form of the Jacobi method that addresses the coupling between the Reynolds and the film thickness equation. A block matrix system by the discrete forms of the two equations can be composed and therefore, the smoother will be referred to as the ‘Block-Jacobi’ method. If the ordinary Jacobi method is applied, only the main diagonal of the determining blockmatrix would be used to give pressure and film thickness increments. This kind of iterative method will not converge. In order to account for the coupling between pressure and film thickness, the Block-Jacobi method makes use of the main diagonals in each of the four blocks (Appendix). The Block-Jacobi method, which works as a smoother for the full multilevel algorithm, is stable but has the disadvantage of slow convergence. When the Block-Jacobi method is applied to solve the 2n − by − 2n (n is the number of spatial nodes) block matrix system, it is reduced to a 2 − by− 2 matrix system Eq. (A.12). Pi fi ai bi = , (A.12) ci d i Hi gi with coefficients given in the Appendix. The nature of the EHL line contact problem introduces some difficulties that have to be dealt with. The convergence of the force-balance equation is very important in order to obtain sufficient accuracy of h 0 within each time step. This can be seen from the pressure-part solution of the coupled system Eq. (A.13), viz. di f i − bi gi Pi = , (A.13) a i d i − b i ci where the nominator contains the time dependence of the Reynolds equation. The nominator of Eq. (A.13) expresses numerically the time derivative, in dimensionless form, of the film thickness and thus the time derivative of h 0 (t) . That is, a small ∆t will enhance the influence of an error in h 0 (t) and the physical explanation of this behavior is the squeeze mechanism represented by ∂/∂t in Eq. (A.4). Moreover, this is not only important when using a coupled solution method; a serial solution method will be affected in the same way.

A.3.3 Error estimation In order to compare the solutions between the CFD and Reynolds approaches, it is important to estimate the errors. The strategy used in this work is to check the maximum error in the whole domain (both spatially and in time). For the Reynolds approach, the defined error functions are two-dimensional. In addition, the CFD approach has one more dimension, i.e. the error is measured across the film as well. All errors measured were obtained by halving the step size. The time error measure is defined as: ∆t (A.14) εtime (∆t) = max|Φi, j,2n − Φi, j,n (∆t) |. i, j,n 2 That is the solution profiles Φ i, j,k (pressure) is considered as a function of the step size (∆t) in time where i, j and k is suffices for x, z and t respectively. However, this error will of course include the error in the spatial domain, e.g., the finer the spatial step size (∆x and/or ∆z) the smaller the error ε time . The same principle was adopted to define the spatial errors by considering the solution profiles as functions of the spatial step size. ∆x 1 εx (∆x) = max|Φ2i, j,n − Φi, j,n (∆x) |. (A.15) i, j,n 3 2

A.3. NUMERICS

79 1 εz (∆z) = max|Φi,2 j,n 3 i, j,n

∆z − Φi, j,n (∆z) |. 2

(A.16)

These errors Eq. (A.15) and Eq. (A.16), include the time errors in the same way as Eq. (A.14) includes the spatial errors. In addition, both the time- and the spatial- error measures include a convergence error since the solution processes in both the Reynolds and the CFD approach are not numerically exact. In the Reynolds approach, the convergence error is the maximum residual value of the Reynolds equation at each time step. This is sufficient since the residual value, and thus the errors of the film thickness equations, are small in comparison. The CFD approach controls convergence by summarizing the errors at each control volume. It was found that in both approaches there exists a pronounced relationship between the step size in time and the spatial step size ∆x. That is, for the running conditions specified in Table A.1 u s ∆t/2 has to be less or equal to the step size ∆x. Applying this condition yields the appropriate behavior of the time error Eq. (A.14), i.e., εtime ∆t2 1 . (A.17) εtime (∆t) 2 However, in the case of a slowly moving ridge u ridge ≤ u/4 the transient effects are well resolved even with values of u s ∆t/2 larger than ∆x.

A.3.4 Interpolation of solution data The numerical set of data representing the solutions obtained from the different approaches is not defined on exactly the same grid points nor in exactly the same points in time. In order to compare the solutions it is necessary to evaluate them on the same grid points. Here, linear interpolation of the Reynolds solutions on to the CFD solution grid points, both spatially and in time, is performed.

Table A.1: Input data used in the simulations Parameter Ru Rd µ0 zvisc Load E ν τ0 ρ0 λ A

Value 0.01 2.0 0.14 0.6 105 206 0.3 3.3 870 0.1 0.2

Dimension m m Pa s N m−1 GPa MPa kg m −3 mm µm

PAPER A.

80

A.4

Results

In this section results from the numerical simulations are given. In the figures shown, the CFD solution is denoted by circles and the Reynolds equation solution by solid lines. In total, three running conditions are compared: • A smooth line contact under sliding • A ridge attached to the faster surface • A ridge attached to the slower surface

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2

0

-0.2

-0.1 x (mm)

0

0.1

P (GPa)

H (µm)

All three take place under the same amount of sliding (|s| = 5/3). In Fig. A.1, a smooth contact comparison is shown. The upper surface moves with a velocity of 0.1 m/s and the lower surface with 1.1 m/s. As can be observed, the two solutions correspond well to each other. According to the figure, the pressure spike is not present and the exit constriction at the outlet has a smooth geometrical curvature due to the strong non-Newtonian behavior of the lubricant. Fig. A.2 shows intermediate time steps of the passage of a ridge through the

0

Figure A.1: Pressure- and film thickness- profiles for a smooth contact, s = ±5/3. CFD solution (o), Reynolds (-). EHL conjunction. The upper surface with the ridge moves with velocity 0.1 m/s, and the lower surface at 1.1 m/s. Once again the deviations are very small. The vertical line shown in the figure indicates the center of the undeformed ridge and, under running conditions as presented in Table A.1, no phase difference between the ridge and film disturbance can be observed [40]. This is attributed to the moderate pressure and the strong non-Newtonian behavior. Finally, A.3 again shows intermediate time steps for passage of a ridge through the EHL conjunction. The ridge is attached to the upper surface and moves with velocity 1.1 m/s. The lower surface moves at 0.1 m/s. The maximum deviation between the two approaches is approximately 1%. This is, however, smaller than the accumulated errors in the numerical solutions. So, if finer time and/or spatial resolutions are used in comparison with stronger demands of the convergence error, this deviation will probably vanish. The full set of data used for the simulations is contained in Table A.1, for time steps and meshes in Table A.3 and in Table A.2 for the discretization errors. The discretization error is now denoted by ε spatial and is the maximum of the errors ε x and εz .

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2 0 1 0.8

1.2

0.6

0.8

0.4

0.4

0.2 -0.2

-0.1

0

0.1

0

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2 -0.2

-0.1

0

0.1

0

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2

0

-0.2

-0.1

0

0.1

0

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2

0

-0.2

-0.1 x (mm)

0

0.1

P (GPa)

0.1

P (GPa)

H (µm)

0

2

0

H (µm)

-0.1

1.6

0

H (µm)

-0.2

P (GPa)

H (µm)

0

P (GPa)

81

P (GPa)

H (µm)

A.5. DISCUSSION

0

Figure A.2: Pressure and film thickness profiles at five intermediate time steps when the ridge is attached to the slower moving surface, u ridge = 0.1 m/s, s = −5/3. CFD solution (o), Reynolds (-).

A.5

Discussion

The numerical experiments performed so far show very small deviations between the two solution approaches. The geometries used in this work have a small film thickness to wavelength ratio, which means that the thin film approximation is still accurate. If this ratio is enlarged, one may expect larger deviations. From a CFD viewpoint, the results are encouraging due to the extending possibilities of simulating EHL conjunctions. If the

PAPER A. 1 0.8

1.2

0.6

0.8

0.4

0.4

0.2 0 1 0.8

1.2

0.6

0.8

0.4

0.4

0.2 -0.2

-0.1

0

0.1

0

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2 -0.2

-0.1

0

0.1

0

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2

0

-0.2

-0.1

0

0.1

0

2

1

1.6

0.8

1.2

0.6

0.8

0.4

0.4

0.2

0

-0.2

-0.1 x (mm)

0

0.1

P (GPa)

0.1

P (GPa)

H (µm)

0

2

0

H (µm)

-0.1

1.6

0

H (µm)

-0.2

P (GPa)

H (µm)

0

P (GPa)

2 1.6

P (GPa)

H (µm)

82

0

Figure A.3: Pressure and film thickness profiles at five intermediate time steps when the ridge is attached to the slower moving surface, u ridge = 1.1 m/s, s = +5/3. CFD solution (o), Reynolds (-).

numerically complex flow problems inside the contact regions can be simulated, the full inlet and outlet flow should also be possible to simulate. Such simulation will give more clarity in the reformation process or when starvation may be expected. How wear debris travels through the EHL conjunctions (particle transport simulation) is another subject for study. The comparison also enhances the reliability of the Reynolds approach, that even in a case with large non-Newtonian effects and predicts approximately the same pressure and

A.5. DISCUSSION

83

Table A.2: Discretization errors (%) Parameter s εtime CFD εtime Reynolds εspatial CFD εspatial Reynolds

Case 1 ±5/3 0.34 0.021

Case 2 −5/3 0.49 0.057 1.1 0.12

Case 3 +5/3 3.61 1.04 1.1 0.17

Table A.3: Meshes and time steps used in the simulations Parameter Mesh CFD Mesh Reynolds ∆x CFD ∆x Reynolds ∆t CFD ∆t Reynolds

Case 1 400x20x1 1024x0x1 1.5 0.568 -

Case 2 400x20x1 1024x0x1 1.5 0.568 60.0 0.948

Case 3 400x20x1 1024x0x1 1.5 0.568 1.8 0.474

Dimension µm µm µs µs

film thicknesses as the CFD solution. In particular, this work has also shown that it is ˝ possible to use the Block UJacobi relaxation method as a smoother for the multilevel FAS algorithm. The choice of step size both spatially and in time, is crucial when it comes to numerical convergence in both approaches. In this work the passage of a single ridge with small amplitude to wavelength ratio has been simulated. If the wavelength is kept constant and the amplitude is enlarged one would most certainly have to consider the speed of the asperity instead of the surface mean speed when choosing the step size in time, for example, uridge ∆t ≤ ∆x. On top of this, asperities with enlarged amplitude to wavelength ratios have to be well resolved, meaning small spatial step sizes, and thus even smaller step sizes in time are required. This applies to situations where one of the surfaces can be considered as perfectly smooth. Other situations arise when dealing with real applications. In real applications both surfaces are rough. Combined with sliding, this leads to a continuously changing effective surface roughness inside the EHL conjunctions. The squeeze effects are underestimated if the surface roughness is imposed on one of the contacting surfaces. Numerically, using either of the two approaches, it is necessary to consider the amplitude to wavelength ratio as time dependent and hence the spatial step size should be chosen so that the largest amplitude to wavelength ratio in time is resolved. The rheological model used in this work is the Eyring model which has a reduced shear stress to shear rate dependence when the pressure and shear rates are high. In order to enable CFD simulations at high contact pressures, some form of damped relationship between shear stress and shear rate is necessary in order to avoid a singularity in the momentum equation Eq. (A.2). However, in the CFD approach there are enhanced flexibilities to allow for modification of the rheology to include more advanced models. (This is because the momentum equation is used in its basic form.) The numerical simulations performed in this work are isothermal. It is, however, nec-

PAPER A.

84

essary to include thermal computations in order to have an appropriate fluid flow model. Such a project is ongoing and will hopefully add more understanding to transient EHL.

A.6

Conclusions

• Considering numerical results, deviations between the two approaches are small. It must be emphasized here that the amplitude to wavelength ratios in this work are small. • The numerics of the two approaches seem stable. However the case when the ridge is attached to the faster moving surface is somewhat more difficult to simulate due to the greater influence of the transient effects. ˝ times • The deviation in pressure between the two approaches is approximately 10 U20 higher compared to the deviations in film thickness. This is also the case for numerical errors measured for each approach separately.

Acknowledgment The authors gratefully acknowledge the financial support from the Swedish Research Council (VR), the National Graduate School in Scientific Computing (NGSSC) and the national research programme HiMeC financed by the Foundation of Strategic Research and Fortum, Indexator, SKF, Statoil, Volvo AB and Volvo Car Corp.

Appendix The dimensionless discrete form of the time-dependent Reynolds Eq. (A.4) can be written as: 1 k Pk − εk k k + εk Pk + ε ε P − 2 i i− 1 i−1 i− 1 i+ 1 i+ 1 i+1 (∆X)

2

λ 2∆X

2

2

2

k k k +ρ k ¯ ki−2Hi−2 3ρ¯ i Hi − 4ρ¯ ki−1Hi−1 − λ ∆T

where εki± 1 2

(A.18)

Hik−1 = 0, ρ¯ ki Hik − ρ¯ k−1 i


(A.19)

Sik is the non-Newtonian factor after Johnson and Tevaarwerk [10]. For the film thickness the following dimensionless form is used: Hik − H0k −

(Xi )2 1 + 2 π

∑ Ki j Pkj + ψ¯ ki = 0, j

(A.20)

A.6. CONCLUSIONS

85

where Ki j

∆X ∆X − 1 − Xi − X j + ln Xi − X j + 2 2

=

∆X ∆X − 1 . Xi − X j − ln Xi − X j − 2 2

(A.21)

The equations (A.18) and (A.20) can be written in operator form if the k − 1 term in Eq. (A.18) is moved to the right hand side, viz. Lk1 Pk , H k = Lk2

P ,H k

k

F k−1 (A.22)

=

G. k

If the coefficients in Eq. (A.22) are linearized by the previously determined values of P and H obtained from any relaxation process, then it would be possible to rewrite Eq. (A.22) as a block matrix system, i.e.,

A11 A21

A12 A22

k Pnew k Hnew

=

F k−1 Gk

.

(A.23)

If the Jacobi method is applied to solve the above block-matrix system then the updated, ‘new’ pressure and film thickness would be obtained from the following equations: k diag(A11 ) Pnew

k k = F k−1 − A∗11 Pold − A∗12 Hold

k diag(A22 ) Hnew

k k = Gk − A∗21 Pold − A∗22 Hold .

(A.24)

where diag(Alm ) are diagonal matrices with the main diagonal of A lm . A∗lm is the same as Alm except that the main diagonal is zero. The system identified in Eq. (A.24) does not address the coupling between pressure and film thickness. Thus, a Jacobi relaxation of the ˝ blockUmatrix system Eq. (A.23) will not converge. However, the system can be modified to address the coupling, e.g., k k diag(A11 ) Pnew + diag(A12 ) Hnew


k k + diag(A22 ) Hnew diag(A22 ) Hnew

= G

(A.25) k

k − A∗21 Pold

k − A∗22 Hold .

The Block-Jacobi relaxation method is based on solving the block-matrix system Eq. (A.23) for the ‘new’ pressure and film-thickness profiles. It is therefore straightforward to rewrite the system to indices form, i.e. to the two-by-two matrix system Eq. (A.12), viz.

ai ci

bi di

Pi Hi

=

fi gi

.

PAPER A.

86 with coefficients: ai

=

bi

=

k k ε + ε 1 i+ 21 (∆X)2 i− 2 3 1 + −λ ρ¯ ki 2∆X ∆X 1

−

ci

=

∆X ∆X ln 2

bi

=

1

1

(A.26)

and fi

=

− −

gi

=

(∆X)2 λ

k k εki− 1 Pi−1 + εki+ 1 Pi+1 2

2


H0k +

Xi2 1 − 2 π

∑ Ki j Pkj − ψ¯ ki . j =i

(A.27)

Paper B

87

89 Proceedings of the 30th Leeds-Lyon Symposium on Tribology

The effect of two-sided roughness in rolling/sliding ehl line contacts A. Almqvist 1 and R. Larsson Luleå University of Technology, Department of Mechanical Engineering, Division of Machine Elements, Luleå, SE-971 87 Sweden

Abstract In most theoretical studies carried out to date on the effect of surface roughness in elastohydrodynamic lubrication (EHL) one surface is considered smooth and one as being rough. In real tribological contacts however, both surfaces normally have similar roughness heights. When modeling a rolling contact it is possible to simply sum the roughness of the two contact surfaces but in a sliding EHL contact, a continuously changing effective surface roughness occurs. The aim of this work was to investigate the influence of elementary surface features such as dents and ridges on the film thickness and pressure. This was done numerically using transient non-Newtonian simulations of an EHL line contact using a coupled smoother combined with a multilevel technique. Four different ‘overtaking’ phenomena were investigated; ridge-ridge, dent-ridge, ridge-dent, and dent-dent. It was shown that the minimum film-thickness produced by a ridge is further reduced in a dent-ridge overtaking event. The squeeze effect seen in the ridge-ridge case resulted in large deformations and film thickness heights comparable to the corresponding smooth case just before the overtaking event occurred. These local effects arising from simulating two-sided roughness were compared to simulations using a traditional ‘one-sided rough surface contacting a perfectly smooth surface’. Keywords: Elastohydrodynamic lubrication, Transient effects, Coupled method, Non-Newtonian, Surface roughness

1 Corresponding author. E-mail address: [email protected].

90

PAPER B.

B.1 Introduction A machined surface is never perfectly smooth. The lubrication of such surfaces in elastohydrodynamic lubrication (EHL) contacts is normally a highly transient process. Two different types of transient effects occur in an EHL line contact subjected to a rolling/sliding motion. In the first case, pure rolling, the film thickness will change transiently at each position along the contact as a result of the shifted roughness profile. It is possible to simplify the theoretical model of the pure rolling case by assuming all roughness to be located at one surface while letting the other surface be smooth. The roughness profile used in this model is simply the sum of roughness profile heights of both surfaces. In the second case, a combined rolling and sliding motion, the roughness profiles of each surface shift relative to each other resulting in asperity ‘overtaking’ phenomena occurring within the contact. Most simulations of EHL with rough surfaces assume the rolling case, i.e. that the roughness is located on only one of the surfaces [43, 44, 45, 11]. The reason for this is that it becomes much more difficult to define a simple model if the roughness and relative movement of both surfaces has to be considered. Clearly, a realistic analysis of a contact experiencing a sliding motion is not possible if the second transient effect is not considered. Changes in the film thickness in cases with roughness on only one surface are very much controlled by the local inlet conditions. One example is the distortion phenomena where the passage of an asperity causes a distortion in the film thickness which travels through the contact at the speed of the fluid and not the speed of the asperity through the contact. However, in the case where both surfaces are assumed to be rough and have a relative motion, local effects may occur inside the EHL conjunction that cause film breakdown and/or extreme pressure peaks. Such effects would be due to the constantly changing sum of the surface roughness when distortions from asperity overtaking phenomena inside the contact area are considered. This would in turn cause changing volumes between the asperities which will influence film pressure and also bearing capacity. A number of investigations taking into account roughness on both surfaces have been presented, for example, [1, 7, 8]. However no systematic study of the two-sided roughness problem (such as has been done by Venner [45] for the one-sided roughness case) has been attempted. This paper presents our first attempt to make such a systematic study of two-sided roughness. In order to describe the pressure and film thickness between sliding rough surfaces in a highly pressurized hydrodynamically lubricated conjunction, the numerical solver has to be able to handle transient effects. Moreover, if the simulation aims to model the function of a real application, the solver needs to include a non-Newtonian description of the fluid, internal cavitation, thermal effects, plastic deformations, etc. There are still much to be done in this area, but the character of many of these effects makes it possible to model them in simplified way whilst still obtaining valuable information. The most commonly used numerical model of an EHL line contact is based on Reynolds equation coupled with an equation for the film thickness, empirical laws for viscosity and density, and a forcebalance equation. See, for example, Venner and Lubrecht [11] for a detailed view on multi-level techniques for the numerical solution of the equations. The transient squeeze effect that occurs when a dent and ridge overtake each other, almost causing film breakdown, is the main focus of this work. This event was chosen since it would intuitively not be considered harmful. However, simulations show that a

B.2. THEORY

91

film breakdown can occur as a result of such an overtaking event. Five different slide-toroll ratios, a test with a doubled load and a test with doubled ridge amplitude were carried out. It was shown that, for all amplitude cases, there exists a wavelength for which the dent ridge overtaking produces a minimum film thickness.

B.2 Theory This work focuses on one of the four different overtaking phenomena that may occur in an EHL line contact. Pressure build-up and film formation related to these phenomena could have been studied without taking into account non-Newtonian and thermal effects in order to simplify the analysis. However, Non-Newtonian effects are included according to Conry et al. [2]).

B.2.1

Equations

The line contact problem to be solved comprises the time dependent non-Newtonian Reynolds equation (B.1), ∂ ρh3 ∂p us ∂ (ρh) ∂ (ρh) S − =0 (B.1) − ∂x 12µ ∂x 2 ∂x ∂t where S =

3 (Σ cosh Σ − sinhΣ) Σ3

Σ =

h p , 2τ0 x

1+

µ2 (u2 − u1)2 (τ0 h)

2

Σ2 , sinh2 Σ

the film thickness equation (B.2), h (x,t) = h0 (t) +

x2 + d (x,t) + ψ1 (x,t) + ψ2 (x,t) , 2Rx

(B.2)

where the elastic deformation is given by d (x,t) =

4 πE

∞ −∞

p x ,t ln |x − x|dx .

Roelands equations for the viscosity (B.3) αP0 p zvisc µ (p) = µ0 exp −1 + 1 + . zvisc P0

(B.3)

where α, P0 and zvisc are mutually dependent: P0 =

zvisc (ln (µ0 ) + 9.67), α

and Dowson and Higginson equation for the density (B.4) ρ (p) = ρ0

C1 + C2 p , C1 + p

(B.4)

PAPER B.

92 At all times, the force balance condition (B.5)

I

p (x,t) dx = w (t) .

(B.5)

determines the integration constant h 0 (t). Also, the cavitation condition (B.6) is used to ensure that all negative pressures obtained during the solution process are removed. P ≥ 0.

(B.6)

The surface dents/ridges were modeled by Eq. (B.7) −7

ψi (x,t) = 2 · 10 li 10

−10

x−xc (t) 10−4 mi

2

2π (x − xc (t)) cos , 10−4mi

(B.7)

where the dent/ridge central position is given by x ci (t) = xsi + uit, and xsi is the initial position. To simplify varying the amplitudes and wavelengths of the surface dents/ridges, two help parameters where introduced; l i and mi . A variation of l i corresponds to a variation in the amplitude of the dent/ridge whilst a variation of m i corresponds to a variation in the wavelength of the dent/ridge. The line contact equations (B.1)-(B.6) were transformed into dimensionless form by the set of dimensionless variables: X = x/b H = h/ b2 /Rx P = p/ph ρ¯ = ρ/ρ0 , T = t/ (2b/us) µ¯ = µ/µ0

where b=

B.2.2

2w 8wRx , ph = πE πb

Numerics

In the Reynolds approach used in this work, the Reynolds equation and the equation for the film thickness are coupled so that pressure and film thickness are solved simultaneously (Evans et al. [1], Holmes [46]). Multilevel technique (Venner, Venner and Lubrecht [44, 45, 11]) have been used to accelerate the convergence of the Reynolds approaches based on serial solution methods, i.e., the pressure is the only unknown parameter. Here the multilevel technique is extended to accelerate the solution process of the coupled approach. The full multilevel algorithm is based on a specially developed iterative relaxation process. This relaxation process is a generalized form of the Jacobi method that addresses the coupling between the Reynolds and the film thickness equations. A block matrix system of the discrete forms of the two equations can be composed and hence the smoother will be referred to as the ‘Block-Jacobi’ method. If the ordinary Jacobi method is applied, only the main diagonal of the determining block-matrix would be used to give pressure and film thickness increments. This kind of iterative method would not converge. In order to account for the coupling between pressure and film thickness, the Block-Jacobi method makes use of the main diagonals in each of the four blocks (see Appendix). The Block-Jacobi method is stable but has the disadvantage of slow convergence. Convergence can, however, be accelerated by using

B.2. THEORY

93

the applied multilevel technique. When the Block-Jacobi method is applied to solve the 2n-by-2n (n is the number of spatial nodes) block matrix system, it is reduced to a 2-by-2 matrix system Eq. (B.8). Pi fi ai bi = , (B.8) ci d i Hi gi with coefficients given in the Appendix. The nature of the EHL line contact problem introduces some difficulties that have to be dealt with. The convergence of the force-balance equation is very important in order to obtain sufficient accuracy of h 0 within each time step. This can be seen from the pressure solution of the coupled system Eq. (B.9), viz. Pi =

di f i − bi gi , a i d i − b i ci

(B.9)

where the nominator contains the time dependence of the Reynolds equation. The nominator of Eq. (B.9) expresses numerically the time derivative, in dimensionless form, of the film thickness and thus the time derivative of h 0 (t) . That is, a small ∆T will enhance the influence of an error in h 0 (t) and the physical explanation of this behavior is the squeeze mechanism represented by ∂/∂t in Eq. (B.1). Moreover, this is not only important when using a coupled solution method; a serial solution method will be affected in the same way. In the EHL conjunction outlet the cavitation condition Eq. (B.6) is imposed. The solution of the Reynolds equation, which admits a negative pressure, is not valid here and the resulting non-Newtonian parameter S will therefore drop to zero. Correction of S by forcing S = 1 when ∂p/∂x = 0 leads to an effective viscosity µ/S which equals the ambient viscosity µ0 . This gives an almost smooth extension of S in the region of zero pressure. The importance of such an extension relates to the Block-Jacobi relaxation of Eq. (B.8), since Eq. (B.6) implies P = 0 in the cavitation zone and the whole domain is considered at all iterations. In the solution process, this correction successively moves the starting point of the cavitation zone to a certain, converged, point.

B.2.3

Error estimation

In order to validate the results it is important to estimate the errors. The strategy used in this work was to check the maximum error in the whole domain (both spatially and temporally). In the case of a line contact, the defined error functions are two-dimensional. The time error measure is defined by Eq. (B.10) since all errors measured were obtained by halving the step size. ∆t (B.10) − Φi,k (∆t) |. εtime (∆t) = max|Φi,2k i,k 2 That is the solution profiles Φ i, j,k (pressure) is considered as a function of the step size (∆t) in time. The indices i and k represents the spatial and the time coordinate respectively. However, this error will of course include the error in the spatial domain, e.g., the finer the spatial step size ∆x the smaller the error εtime . The same principle was adopted to define the spatial errors by considering the solution profiles as functions of the spatial step size. ∆x − Φi,k (∆x) |. (B.11) εx (∆x) = max|Φ2i,k i,k 2

94

PAPER B.

This treatment of numerical errors is further discussed in [47]. The spatial error includes the time errors in the same way as the time error includes the spatial error. In addition, both the time and the spatial error measures include a convergence error since the solution processes is not numerically exact. The convergence error is the maximum residual value of the Reynolds equation at each time step. This measure of the convergence error is sufficient since the residual value and thus the error in the film thickness equation are small in comparison. It was found that a relationship existed between the step size in time and the spatial step size. That is, for the running conditions in this work, u s ∆t/2 has to be less or equal to the step size ∆x (∆T ≤ ∆X). Applying this condition yields the appropriate behaviour of the time error Eq. (B.10), i.e., εtime ∆t2 1 . εtime (∆t) 2 However, in the case of a slowly moving asperity u asperity ≤ u/4 the transient effects are well resolved even with values of u s ∆t/2 larger than ∆x. In this work the maximum time error in the film thickness was less than 5% of the corresponding smooth minimum film thickness. This maximum time error was found in the simulations with s = 1/3 and with the largest dent amplitude to dent wavelength ratio. All the results shown were simulated with at least four different time steps and in some cases, spatial discretisation was checked by a set of three different spatial mesh sizes. It was found that the time error in all cases was larger than the spatial error even when using the smallest number of nodes at the finest level, e.g., 257. Therefore, all simulations performed have a spatial resolution of 257 nodes since the computational time in this case is decreased by roughly a factor of four compared to a resolution of 1025 nodes.

B.3 Results and discussion The focus of this section is on the film thickness thinning effect caused by a dent-ridge overtaking event. The other three possible cases of surface features overtaking each other have also been briefly studied and are described below. In the case of ridge-ridge overtaking the pressure spike produced is by far the largest, when comparing surface features of the same amplitude and wavelength, and this type of overtaking also produce the thinnest film thickness. In the beginning of a dent-dent overtaking event, a pressure decrease due to the diverging surfaces occurs and at the end of the event the contact geometry causes a rapidly convergent gap that produce a large pressure spike. The pressure spike produced in this case is of comparable magnitude to the pressure spike in the ridge-ridge event with ridges of half the amplitude. The case of the ridge-dent overtaking event was found to have less extreme pressure spikes and minimum film thicknesses. However, it is not possible to investigate such a case without a two-sided surface roughness treatment. Thus, it is possible that the transient effects may be more significant in a contact with roughness on both surfaces. The case studied in more detail, a dent-ridge overtaking event, is most interesting because of the potential film breakdown which may occur as an elastically deformed asperity is allowed to recover to its original height. Because of the lubricant film thinning effect in a dent-ridge overtaking event, several different parameter studies were performed. Five different slide-to-roll ratios (s) were initially selected, starting from 5/3, decreasing to 1/3 in 1/3 steps. Sets of simulations with fixed dent wavelength and varying dent

B.3. RESULTS AND DISCUSSION

95

amplitude were carried out, keeping the ridge amplitude and wavelength constant. Also, for the case of s = 5/3, both dent amplitude and wavelength were varied giving a total of 42 combinations. A set of simulations with doubled load and varying the dent amplitude only, were also carried out. In addition, again for the case of s = 5/3, the dent amplitude was kept constant and simulations with doubled ridge amplitude varying only the dent wavelength were carried out. Fig. B.1 shows intermediate time steps of the dent-ridge overtaking event in the case of dent parameter l 2 = −5, m2 = 1 and s = 4/3. It is clear that the dent causes a thinning effect when overtaking the ridge. This is due to the pressure release, experienced by the already compressed ridge, under the time of the overtaking causing a thinning of the film thickness to occur. To compare the different slide-to-roll ratios, a measure of the minimum film thickness caused by the dent was defined, see Eq. (B.12). max

ξ (l, m) =

−1.3≤X≤0.0

|H l,m (X, T ) − H re f (X, T ) | min|H l,m (X, 0) |

(B.12)

X

Simulations of the ridge passage only for each of the slide-to-roll ratios tested were used as a reference for the dent-ridge simulations. The exact position of the overtaking was set as −0.3253 · b, in order to have a suitable ‘window’ for all slide-to-roll ratios tested. It is important to remember that this measure assumes a large value when the actual film thickness is small. The set of simulations where both dent amplitude (l 2 ) and dent wavelength (m2 ) were varied are shown as a contour map of the measure ξ in Fig. B.2. Values of l 2 and m2 for the dent were varied to give a total of forty-two combinations in the ranges l2 = −1 . . . − 6 and m2 = 1 . . . 7 respectively. Ridge parameters were held constant, l 1 = 1 and m1 = 1. The figure shows that the larger the amplitude of the dent, the larger were the effects of the dent and thus the thinner the film thickness. As far as the dent wavelength is concerned, an optimal and worst case wavelength in the range of wavelengths simulated were found. This is interesting since in many tribological applications dents are thought of as lubricant depots and these simulations show that there is an optimal dent wavelength to maintain the greatest film thickness possible, at least when considering an EHL line contact. When comparing the different slide-to-roll ratios, the dent amplitude (l 2 ) was varied as in Fig. B.2 but here the dent wavelength (m 2 ) was held constant, i.e., m 2 = 1. Ridge parameters were also kept constant using the same values as previously. As Fig. B.3 shows, the film thinning effect is clearly dependent of the slide-to-roll ratio and that the effects increase with increasing s. The effects caused by the dent are of both greater magnitude and increase more rapidly with dent l 2 with a higher slide-to-roll ratio. The overall maximum of ξ is approximately 0.4, i.e., the dent causes a film thickness decrease that corresponds to 40% of the smooth surface minimum film thickness value. In the case of s = 5/3, the effect of ridge amplitude, i.e., the cases of l 1 = 1 and 2, where also studied. In both cases m 2 was varied in the range 1, 2, . . . , 7 and the dent parameter l2 were taken as −l1 . The results are shown in Fig. B.4. These results confirms that the measure ξ assumes both a local maximum and a local minimum. That is, a dent wavelength (m2 ) somewhere between 1 and 3 yield the thinnest film and a value somewhere between 4 and 6 yield the thickest film. Fig. 6.9 also shows that the film thinning effect is greater in the case of doubled ridge amplitude (l 1 = 2). Doubled ridge amplitude does not increase the film thinning effect more than a few percent in the parameter range chosen. However,

PAPER B.

96 H Hre f P

1.5 1 0.5 0 -2

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5

-1

-0.5

0

0.5

1

1.5

1.5 1 0.5 0 -2 1.5 1 0.5 0 -2 1.5 1 0.5 0 -2

X

Figure B.1: Dimensionless film thickness (−), reference film thickness (−·) and pressure (−;) at four different time steps showing the dent caused thinning effect for s = 4/3.

the difference seem to increase with increasing wavelength. Together with the results shown in Fig. 6.7, this indicates a tendency toward film breakdown if one or a combination of the parameters, ridge amplitude (l 1 ), dent amplitude (l 2 ) and dent wavelength (m 2 ) are combined in a unfortunate way. In real applications, such combinations will occur frequently. In order to investigate the effect of load (w), another set of simulations with s = 5/3 using the same set of dent and ridge parameters as those in the parameter study of the slideto-roll ratio, shown in Fig. 6.8, were carried out. The results, shown in Fig. 6.10, indicate that a doubling the load gives rise to thinner films. It can also be seen that the minimum film thickness corresponds to approximately 45% of the smooth film thickness and that the slope of ξ is greater in the case of the doubled load. In all these simulations, several significant input parameters where held constant, some of these are shown in Table B.1. It is not possible to draw general conclusions from the investigation presented here, but a trend toward film thinning when an asperity is free to recover its shape in the valleys or dents between asperities on the opposing surface. Clearly, at a certain combination of lubricant and roughness parameters film breakdown will occur. Such a breakdown will be due to the low viscosity regions, leading to local low pressure or even internal cavitation,

B.3. RESULTS AND DISCUSSION -1

0.15

0.2

0.2

5

0.25

-2

-4

0.3 5

l2

0. 3

0.2 -3

97

0.2 0.25

0.3

0.25

0.35

0.3

0.4

-6 1

0.3

0.3

5

0.4

-5

0.15

0.15

0.35

0.4

0.45 2

3

4 m2

5

6

7

Figure B.2: Contour map of ξ for s = 5/3, y-axis corresponds to dent amplitude, x-axis to dent wavelength. l 1 = 1 and m1 = 1. 0.5

0.4

s = 5/3 s = 4/3 s=1 s = 2/3 s = 1/3

ξ

0.3

0.2

0.1

0 -1

-2

-3

-4

-5

-6

l2

Figure B.3: ξ as a function of the slide-to-roll ratio (s). x-axis corresponds to dent amplitude. in the valleys or dents rather than to very high contact pressures as might intuitively be thought to be the cause of a film breakdown. In the different cases presented here, no actual film breakdown was observed. However, in a line contact of finite length the film thinning will be greater. This is due to side leakage from the valleys which lead to even lower pressures allowing asperities to expand more rapidly and to a greater extent. This will eventually lead to film breakdown. Further investigations on internal cavitation phenomena and side leakage effects is therefore requested to achieve knowledge about the film thinning effect studied here. It is, of course, also necessary to introduce more realistic surface roughness features. The smaller the wavelengths become, the smaller the amount of sliding is required to cause this type of local film thinning.

PAPER B.

98 0.3

l1 = 1 l1 = 2

0.28

ξ

0.26 0.24 0.22 0.2 0.18 1

2

3

4 m2

5

6

7

Figure B.4: ξ as a function of ridge amplitude for s = 5/3. x-axis corresponds to dent wavelength. 0.5

w = 10 kN w = 20 kN

ξ

0.4

0.3

0.2

0.1 -1

-2

-3

-4

-5

-6

l2

Figure B.5: ξ as a function of load for s = 5/3. x-axis corresponds to dent amplitude.

B.4 Conclusions • Certain effects in a rolling/sliding EHL contact cannot be studied without using a two-sided roughness treatment. • A dent-dent overtaking event gives rise to a pressure spike equivalent to a ridgeridge overtaking event where the ridges have half the amplitude of the dent. Bearing capacity is therefore not significantly affected by this type of overtaking event. • A ridge-ridge overtaking event produces the thinnest films and largest pressure spikes compared to the other three possible overtaking phenomena for a given asperity amplitude. However, film breakdown could not be simulated in this work because of the model used.

B.4. CONCLUSIONS

99

Table B.1: Fixed parameters Parameter R1 R2 E1 E2 ν1 ν2 µ0 α p0

Value 0.01 2.0 206 · 109 206 · 109 0.3 0.3 0.14 2.10 · 10−8 1.98 · 108

Dimension m m Pa Pa Pa s 1/Pa Pa

• Local pressure effects resulting from a ridge-dent overtaking event are of comparable magnitude to the dent-dent case. Because of the involvement of a ridge, the film becomes slight thinner causing this case to have a somewhat greater negatively effect on bearing capacity. • In a sliding contact, a dent overtaking a ridge can cause film breakdown. This could occur with too large a ridge, a too wide or deep dent, too high a load, or a combination of these parameters. • There is a dent wavelength in a dent-ridge overtaking event that optimizes film thickness. • As the slide-to-roll ratio increases, the film thickness decreases in a dent-ridge overtaking event. • A ridge with larger amplitude will cause only slightly thinner films. • The contact load influences the film thinning effect. Load related effects are somewhat larger than those for increasing ridge amplitude.

Appendix The dimensionless discrete form of the time-dependent Reynolds Eq. (B.1) can be written as: 1 k Pk − εk k k + εk Pk + ε ε P − 2 i i− 1 i−1 i− 1 i+ 1 i+ 1 i+1 (∆X)

2

λ 2∆X

2

2

2

k k k +ρ k ¯ ki−2Hi−2 3ρ¯ i Hi − 4ρ¯ ki−1Hi−1 − λ ∆T

where εki± 1 2

(B.13)

Hik−1 = 0, ρ¯ ki Hik − ρ¯ k−1 i


(B.14)

PAPER B.

100

Sik is the non-Newtonian factor after Johnson and Tevaarwerk [10]. For the film thickness the following dimensionless form is used: Hik − H0k − where Ki j

=

(Xi )2 1 + 2 π

∑ Ki j Pkj + ψ¯1ki + ψ¯2ki = 0,

(B.15)

j

∆X ∆X −1 − Xi − X j + ln Xi − X j + 2 2 ∆X ∆X −1 . Xi − X j − ln Xi − X j − 2 2

(B.16)

The equations (B.13) and (B.15) can be written in operator form if the k − 1 term in Eq. (B.13) is moved to the right hand side of the equation and the terms not depending on the pressure and/or the film thickness is moved to the left hand side of B.13 viz. Lk1 Pk , H k = F k−1 Lk2

P ,H k

k

(B.17) =

G. k

If the coefficients in Eq. (B.17) are linearized by the previously determined values of P and H obtained from any relaxation process, then it would be possible to rewrite Eq. (B.17) as a blockmatrix system, i.e., k k−1 Pnew F A11 A12 = . (B.18) k A21 A22 Hnew Gk If the Jacobi method is applied to solve the above block-matrix system then the updated, ‘new’ pressure and film thickness would be obtained from the following equations: k diag(A11 ) Pnew


k diag(A22 ) Hnew

= G

(B.19) k

k − A∗21 Pold

k − A∗22 Hold .

where diag(Alm ) are diagonal matrices with the main diagonal of A lm . A∗lm is the same as Alm except that the main diagonal is zero. The system identified in Eq. (B.19) does not address the coupling between pressure and film thickness. Thus, a Jacobi relaxation of the ˝ blockUmatrix system Eq. (B.18) will not converge. However, the system can be modified to address the coupling, e.g., k k + diag(A12 ) Hnew diag(A11 ) Pnew


k k + diag(A22 ) Hnew diag(A22 ) Hnew

= G

(B.20) k

k − A∗21 Pold

k − A∗22 Hold .

The Block-Jacobi relaxation method is based on solving the block-matrix system Eq. (B.18) for the ‘new’ pressure and film-thickness profiles. It is therefore straightforward to rewrite

B.4. CONCLUSIONS

101

the system to indices form, i.e. to the two-by-two matrix system Eq. (B.8), viz. ai bi Pi fi = . ci d i Hi gi with coefficients: ai

=

bi

=

k k ε + ε 1 1 i+ 2 (∆X)2 i− 2 3 1 + −λ ρ¯ ki 2∆X ∆X 1

−

ci

=

∆X ln

bi

=

1

∆X 2

(B.21)

and fi

=

− −

gi

=

1 2

(∆X) λ

k k εki− 1 Pi−1 + εki+ 1 Pi+1 2

2


H0k +

Xi2 1 − 2 π

∑ Ki j Pkj − ψ¯1ki − ψ¯2ki . j =i

(B.22)

Paper C

103

105 Proceedings of the 11th Nordtrib Conference

On The Dry Elasto-Plastic Contact of Nominally Flat Surfaces A. Almqvist 1 , F. Sahlin, R. Larsson, S. Glavatskih Luleå University of Technology, Department of Mechanical Engineering, Division of Machine Elements, Luleå, SE-971 87 Sweden

Abstract A model to be used for numerical simulation of the contact of linear elastic perfectly plastic rough surfaces was developed. Energy dissipation due to plastic deformation is taken into account. Spectral theory and an FFT-technique are used to facilitiate the numerical solution process. Results of simulations using 4 two-dimensional profiles with different topographies in contact with a rigid plane for a number loads are reported. From the results it is clear that the real area of contact (A p ) changes almost linearly with load and is only slightly affected by the difference in topography. A plasticity index is defined as the ratio of plastically deformed area (A p ) and Ar . Plastic deformation occurs even at low loads and there is a significant difference in plasticity index between the surface profiles considered. An investigation on how the spectral content of the surface profile influences the results presented is also performed. This is to ensure that the metrological limitations of the optical profiler used to measure the surfaces do not have a significant influence. It is concluded that the highest frequencies of the measured profile have a negligible influence on the real area of contact. Keywords: Dry contact, FFT, Elasto-plastic, Surface roughness

1 Corresponding author. E-mail address: [email protected].

PAPER C.

106

C.1

Introduction

Surface machining has the purpose to produce topographies, in a cost effective way, that has a positive influence on contact conditions. For some applications the product life may be of highest importance and a machining process must be chosen so that wear is minimized. In other applications, low friction is of primary importance and one of the tasks of the surface topography might consequently become to minimize the real area of contact. However, this may significantly reduce service life since a small contact area causes a very high pressure and thus an increased risk of wear. Modeling topographies of real surfaces is a difficult task because of their random structures. Various approaches to model the contact between surfaces have been reported in the literature. Greenwood and Williamson [16] (GW), considered topographies consisting of hemispherically shaped asperities of uniform radius. They assumed that the asperities had a Gaussian distribution, in height, about a mean plane. Greenwood and Tripp [17] extended this approach to handle the contact between two rough surfaces. More improvements to the GW approach were added by: Whitehouse and Archard [18], Nayak [19], Onions and Archard [20], Bush et al. [21, 22] and Whitehouse and Phillips [23, 24]. However, the input needed in these types of models is deduced from surface measurements. More precisely, all of them are based on asperity curvature. This parameter depends on the measurement resolution (Poon and Bhushan [25]) and does also vary in between individual asperities. It is, therefore, a difficult task to find the correct input value for the model. Moreover, use of only a few parameters to describe a surface generates a oneto-many mapping, i.e., the same set of parameters can be deduced for surfaces obtained by completely different machining processes. Along with the accelerating computer development, deterministic approaches to the contact problem have been developed. These are getting more and more realistic as the computational speed increase. However, simplified models of the material and/or the topography are still needed. This is required to minimize the computing time and to be able to study effects independently. For example, the Hertzian contact assumes linear elastic, frictionless materials and parabolic surface profiles for which an analytical solution is available. Westergaard [26] showed that it is possible to solve the contact between a single frequency sinusoidal surface in contact with a rigid plane analytically. On a micro scale these analytical solutions represent quite realistic topographies. However, the coupling to asperities at other scales is neglected and a more complete representation is therefore needed. Except different ways to model the mating materials and the topography, there are numerous deterministic numerical techniques that can be applied to solve contact problems between rough surfaces. Lubrecht and Ioannides [15] applied multilevel techniques to facilitate the numerical solution of the elastic contact problem. Ren and Lee [27] applied a moving grid method to reduce storage of the influence matrix when the conventional matrix inversion approach is used to solve the contact problem of linear elastic bodies with rough surfaces. Björklund and Andersson [28] extended the conventional matrix inversion approach by incorporating friction induced deformations. This was done by using the assumption of linear elastic material. Ju and Farris [29] introduced an FFTbased method to solve the elastic contact problem. Stanley and Kato [4] combined this FFT-based method with a variational principle to solve both, the 2D and the 3D, contact problem of rough surfaces. It should be noted that, the contact between realistic topographies, under relatively

C.2. THEORY

107

small loads, leads to plastic deformations. See for example, Tian and Bhushan [3] who based the theoretical model on a variational principle for both linear elastic and linear elastic perfectly plastic materials. In contrast to the model of Stanley and Kato, the present model incorporate linear elastic perfectly plastic materials. Energy dissipation due to plastic deformation is also accounted for. In this way not only the in-contact topography and the corresponding pressure distribution but also the unloaded plastically deformed topography is obtained. Care is taken to ensure that the metrological limitations of the optical topometer used to measure the topographies do not have a significant influence. That is, because of the measurement size, some low frequency content will be lost and because of limitations in resolving the topography, some high frequency content will be lost. An investigation on how the spectral content of the surface profile, originating from measured topography, influences the results presented is also performed. Both under the assumption of pure elastic material and elasto-plastic material. The elasto-plastic contact between four surface profiles and a rigid plane was also investigated, while applied load was varied.

C.2

Theory

A contact mechanics problem can be governed by the minimum potential energy theory. Pressures and displacements can be obtained by solving a minimum value of an integral equation, i.e., a variational problem, see Kalker [14]. For the variational approach, assuming frictionless linear elastic contact, the equations to be solved are given as

1 pudS + pgdS , min ( f ) = min p≥0 p≥0 2 S S (C.1)

S

pdS = W,

where p is the pressure, u is the normal displacement, g is the gap between the undeformed surfaces and W is the applied load. For two elastic half spaces, assuming plain strain, the total normal surface displacement u (x) in 2D for a given pressure distribution can be obtained from the equation u (x) = −

2 πE ∗

∞ −∞

ln |x − s| p (s) ds + const,

(C.2)

where 1/E ∗ = 1 − v21 /E1 + 1 − v22 /E2 . In the case of a sinusoidal pressure profile p (x) of one frequency only, it is possible to find the exact, closed form solution for the normal surface displacement u (x), see for example [30]. The relationship p (x) = α cos (ωx) ⇐⇒ u (x) =

2α cos (ωx) + const, ωE ∗

where ω = 2πk/L, can be interpreted as the complete contact of two sinusoidal, elastic half spaces, assuming plain strain. It is possible to rewrite this relationship in a more general form U (ω) = H (ω) P (ω) (C.3) where U (ω), P (ω) are the Fourier transformations of normal displacement and pressure respectively and H (ω) = (2α) / (ωE ∗ ) is the transfer function. Note that the constant

PAPER C.

108

determining the absolute value of elastic deformation in Eq. (C.2) is omitted. Eq. (C.3) forms together with Eq. (C.1) the basis of the solution technique to be used in this paper. For the linear elastic contact, the numerical solution procedure proposed by [4] is adopted here. However, some modification of this procedure is necessary to solve the contact problem assuming linear elastic - perfectly plastic materials. To account for the energy dissipation due to plastic deformations, the gradient of f has to be modified according to ∇ f = u e + u p + g, where ue and u p is the elastic and the plastic displacement respectively. The procedure is then to be iterated until u e , u p and p satisfy Eq. (C.1) under the constriction 0 ≤ p ≤ Hs , where Hs is the hardness of the softer material. In this way not only the linear elastic perfectly plastic loaded topography and the corresponding pressure distribution, but also the unloaded plastically deformed profile are obtained. A full description of the method and numerical technique will be presented in a separate paper.

C.3

Surface characterization

The limitations of the surface measurement apparatus may introduce difficulties when characterizing surfaces. These may lead to different characterizations of the same surface, which in turn may lead to contradictive conclusions on how the topography influences the contact conditions. However, depending on the sought parameters, it is possible to overcome the problem of measurement resolution. For a 2D profile, it is possible to filter the topography through its discrete fourier series, as shown by gn (x) =

n

∑ ak sin (2πkx/L) + bk cos(2πkx/L) .

k=1

The number of frequencies n ≤ N/2 − 1, where N is the number of sampling points. In this work N = 512. Fig. C.1 shows how the real area of contact both in the case of linear elasticand linear elastic perfectly plastic- materials, depends on the number of frequencies (n) used to represent the topography. It is clear that the limitations of the apparatus used do 100

Ar /An [%]

80

Linear elastic perfectly plastic Linear elastic

60 40 20 0 0

100 150 50 Number of frequencies

Figure C.1: Percent real area of contact as a function of the number of frequencies n for surface profile D loaded with 100kN. not affect the percentage of real area of contact A r /An since a constant value is attained for n significantly smaller than N/2 − 1. The lowest value of n, for which A r /An can be

C.4. RESULTS

109

considered as constant, depends on the applied load. However, for the theoretical model to be valid, the ratio, A r /An need to be small, maybe even smaller than for the case of 100kN load presented in Fig. C.1. The surface roughness parameters R a , Rq , skewness (Rsk ), kurtosis (Rk ) and the surface mean separation were also considered to be functions of n, in the same way as described for A r . The results of that investigation showed that all parameters could accurately be predicted by using a surface profile representation with a value of n < N/2 − 1. Hence, for the results presented in this work the measurement resolution can be considered sufficient.

C.4

Results

Two-dimensional profiles from four different steel surfaces [31] are studied, each originating from different machining processes. All of these are hardened and their hardness is for the numerical simulations set to Hs = 4.0 GPa. The modulus of elasticity (E) and the Poisson ratio (ν) were assumed to be 210 GPa and 0.3 respectively for all four surfaces. The surfaces are initially unidirectional ground and the profiles are taken perpendicular to this direction. But as Table C.1 shows, their 2D profiles differ because of additional machining of the surfaces. Due to these different additional machining processes, the surfaces will generate different contact conditions. These contact conditions are here represented by various statistical surface parameters, e.g., R a , Rq , Rsk , Rk , but also the percentage of real area of contact A r /An , the percentage of plastically deformed area A p /An and the mean separation of the surfaces. Mean separation is defined as the distance between the rigid plane and the mean height value of the loaded surface profile. To enrich the data given by Table C.1: Data for the four undeformed surface profiles. # A B C D

Machining Unidirectional ground + Phosphated and dephosphated Unidirectional ground Unidirectional ground + 3 x chemically deburred Unidirectional ground + shot-peened

R a (µm) 0.93

Rq (µm) 1.2

Rsk -1.0

Rk 4.4

0.3 0.14

0.38 0.18

-0.25 -0.022

3.2 3.5

0.35

0.45

-0.069

3.1

Table C.1 the Abbott-Firestone bearing ratio curve and the probability density distribution of the four different topographies are visualized in Fig. C.2. It is clear, both from Table C.1 and Fig. C.2, that surface profile A is much rougher than the other ones. This surface is phosphated and then dephosphated. It inherits some special properties from that, but in the numerical simulations the hardness of this surface is set to 4.0GPa as for the other surfaces. Surfaces B and D differs because of the fact that D is shot-peened. This can be seen in the table as increased R-parameter values for surface profile D, and in the figure both as an overall increase but also as a small peak on the profiles at the upper part of the surface. Surface C differs from B because it is chemically deburred which has smoothed the profile significantly. There is a close relationship between the Abbott-Firestone bearing ratio and the height distribution P (z), i.e., P (z) can be obtained by differentiating the

PAPER C.

110

8

×10−6

8 A B C D

6 z [m]

z [m]

6 4

A B C D

2

2 0 0

4

×10−6

20

40

60

80

100

0 P (z)

%

Figure C.2: Abbot-Firestone bearing ratio (left) and the height distribution (right). Abbot-Firestone bearing ratio with respect to z. Because of this form of relationship P (z) shows an amplification of certain properties of the Abbot-Firestone bearing ratio which makes it a good complement in characterization of surfaces. In Fig. C.3 surface profile D and the corresponding plastic deformation (shaded area) after unloading of 100kN is shown to give some additional insight.

Profile height [m]

1.5

×10−6

1 0.5 0 −0.5 −1 −1.5

1.5

2

2.5

3 x[m]

3.5

4

4.5

5 ×10−4

Figure C.3: Undeformed profile D and the corresponding plastic deformation (shaded area) after unloading of 100 kN. The real area of contact is an important quantity in terms of surface performance under different running conditions. A small real area of contact may reduce friction but could lead to unwanted plastic deformation. The percentage of real area of contact A r /An for the four surfaces is plotted as a function of load in the left of Fig. C.4. It can be seen that the real area of contact increases almost linearly with load. In the right of Fig. C.4 the plasticity index here defined as A p /Ar is plotted for the four surfaces. It should be noted that the relationship between A p /Ar and applied load is quite linear in the high load region. It is also clear that surface profile A is subjected to the largest amount of plastic deformation, which might be explained by its R-parameter values. Even though the R-parameters of surface profile A are significantly larger than those for the other surface profiles they do not differ significantly in terms of the percentage of real area of contact. For surface profile C, the plasticity index is low and the real area of contact is high, which also is suggested by its comparably small R-parameter values. Surface profiles B and D have quite similar R- parameter values, and the reason of the higher A p /Ar ratio for D is probably explained by its slightly lager R- parameter values. In Fig. C.5, the R a - and the Rq -values for the

C.4. RESULTS 10

1

A B C D

8

0.8

6

A p /Ar

Ar /An [%]

111

4 2

0.6 0.4

A B C D

0.2

0 0

2

8

4 6 Load [N]

0 0

10 ×104

2

4 6 Load (N)

8

10 ×104

Figure C.4: Percent real area of contact (left) and plasticity index (right).

four surfaces are plotted as functions of load in the left and right sub-figure. The variation with load is small and the decrease is linear. In terms of magnitude of the R-parameters,

1.2

×10−6

1.2

1

0.6

1 A B C D

0.8 Rq [m]

Ra [m]

0.8

0.4

0.2

0.2 5

Load [N]

10

15 ×104

A B C D

0.6

0.4

0 0

×10−6

0 0

5

Load [N]

10

15 ×104

Figure C.5: Ra (left) and Rq (right).

Table C.1 suggests a ranking of the surface profiles, i.e.: A, D, B, C. From Fig. C.5 it is clear that this ranking is preserved when increasing the load. In the same way, Table C.1 suggests a ranking in terms of the skewness R sk . As shown to the left in Fig. C.6. this ranking is also preserved when increasing the load. However, the Rsk -parameter values seem not to posses the strong linear behavior, as for the R-parameters, except for surface profile A. The kurtosis R k is visualized to the right in Fig. C.6. Considering the ranking suggested by Table C.1, for the kurtosis, a variation can be observed as the R k -parameter values for C rapidly decrease in the low load region when increasing the applied load. From this figure it can also be concluded that load alters the kurtosis of C the most and D the least. Finally, the surface mean separation was considered as a function of load in Fig. C.7. Also here, the ranking in terms of the R-parameter values is kept when varying the applied load.

PAPER C.

112 5

0 A B C D

4.5 Rk

Rsk

−0.5

−1

−1.5 0

4

A B C D

3.5 3 5

Load [N]

10

15 ×104

0

5

Load (N)

10

15 ×104

Figure C.6: Skewness (left) and kurtosis (right).

Mean separation [m]

2

×10−6 A B C D

1.5 1 0.5 0 0

2

4 6 Load [N]

8

10 ×104

Figure C.7: Mean separation as a function of load.

C.5

Conclusions

Results from a numerical method for the contact between linear elastic perfectly plastic rough surfaces is presented. The model takes into account the energy dissipation due to plastic deformation. Based on simulations of four different surface profiles originating from measured topographies the following conclusions can be drawn: • The highest frequencies of the measured profile have a negligible influence on the real area of contact. • The four surface profiles show similar behavior in terms of the real area of contact, the plasticity index, the R-parameters and mean separation, when increasing the applied load. • The real area of contact is almost independent of the R-parameter values. • The plasticity index is highly dependent of the R-parameter values.

Bibliography [1] T. G. Huges R. W. Snidle H. P. Evans, C. D. Elcoate. Transient elastohydrodynamic analysis of rough surfaces using a novel coupled differential deflection method. Proc Instn Mech Engrs, 215, 2001. [2] C. A. Cusano T. F. Conry, S. Wang. A reynolds-eyring equations for elastohydrodynamic lubrication in line contacts. Transactions of ASME, 109:648–658, 1987. [3] Bharat Bhushan and Xuefeng Tian. A numerical three-dimensional model for the contact of rough surfaces by variational principle. ASME Journal of Tribology, 118:33–42, 1996. [4] H. M. Stanley and T. Kato. An fft-based method for rough surface contact. ASME Journal of Tribology, 119:481–485, 1997. [5] T. Almqvist and R. Larsson. The navier-stokes approach for thermal ehl line contact solutions. Tribology International, 35:163–170, 2002. [6] A. A. Lubrecht J. A. Greenwood C. H. Venner, F. Couhier. Amplitude reduction of waviness in transient ehl line contacts. Proceedings of the 23rd Leeds-Lyon Symposium on Tribology, pages 103–112, 1996. [7] H. P. Evans R. W. Snidle J. Tau, T. G. Huges. Elastohydrodynamic response of transverse ground gear teeth. Proceedings of the 28th Leeds-Lyon Symposium on Tribology, 1995. [8] L. Chang. A deterministic model for line contact partial elastohydrodynamic lubrication. ASME Tribology International, 28:75–84, 1995. [9] G. E. Morales-Espejel C. H. Venner. Amplitude reduction of small-amplitude waviness in transient elastohydrodynamically lubricated line contacts. Proc Instn Mech Engrs, 213, 1999. [10] J. L Tevaarwerk K. L. Johnson. The shear behaviour of elastohydrodynamic oil films. Proc R Soc Lond A, 356:215–236, 1977. [11] A. A. Lubrecht C. H. Venner. Multilevel Methods in Lubrication. Tribology Series 37, 2000. [12] B. J. Hamrock. Fundamentals of fluid film lubrication. McGraw-Hill series in mechanical engineering, 1994. 113

114

BIBLIOGRAPHY

[13] H. Moes. Optimum similarity analysis with applications to elastohydrodynamic lubrication. Wear, 159:57–66, 1992. [14] J. J. Kalker. Variational principles in contact mechanics. J. Inst. Math. Appl., 20:199, 1977. [15] A.A. Lubrecht and E. Ioannides. A fast solution of the dry contact problem and the associated sub-surface stress field. AMSE Journal of Tribology, 113:128–133, 1989. [16] J. A. Greenwood and J. B. P. Williamson. Contact of nominally flat surfaces. Proc. Roy. Soc. (London), Series A, 295:300–319, 1966. [17] J. A. Greenwood and J. H. Tripp. The contact of two nominally flat rough surfaces. Proc. Inst. Mech. Eng., 185:625–633, 1971. [18] D. J. Whitehouse and J. F. Archard. The properties of random surface of significance in their contact. Proc. Roy. Soc. (London), Series A, 316:97–121, 1970. [19] P. R. Nayak. Random process model of rough surfaces. ASME J. Lubr. Technol., 93:398–407, 1971. [20] R. A. Onion and J. F. Archard. The ccontact of surfaces having a random structure. Journal of Physics D: Applied Physics, 1973. [21] A. W. Bush, R. D. Gibson, and T. R. Thomas. The elastic contact of rough surface. Wear, 35:87–111, 1975. [22] A. W. Bush, R. D. Gibson, and G. P. Keogh. Strongly anisotropic rough surfaces. ASME J. Lubr. Technol., 101:15–20, 1979. [23] D. J. Whitehouse and M. J. Phillips. Discrete properties of random surfaces. Filos. Trans. R. Soc. (London) Series A, 290:267–298, 1978. [24] D. J. Whitehouse and M. J. Phillips. Two-dimensional discrete properties os random surfaces. Filos. Trans. R. Soc. (London) Series A, 305:441–468, 1982. [25] B. Bhushan and C. Y. Poon. Comparison of surface roughness measurements by stylus profiler, afm and non-contact optical profiler. Wear, 190:76–78, 1995. [26] H. M. Westergaard. Bearing pressure and cracks. ASME Journal of Applied Mechanics, 6:49–53, 1939. [27] N. Ren and S. Lee. Contact simulation of three-dimensional rough surfaces using moving grid method. ASME Journal of Tribology, 115:597–601, 1993. [28] Stefan Björklund and Sören Andersson. A numerical method for real elastic contacts subjected to normal and tangetial loading. Wear, 179:117–122, 1994. [29] Y. Ju and T. N. Farris. Spectral analysis of two-dimensional contact problems. ASME Journal of Tribology, 118:320–328, 1996. [30] K. L. Johnson. Contact Mechanics. The Press Syndicate of the University of Cambridge, 1985.

BIBLIOGRAPHY

115

[31] J. Lord and R. Larsson. Film forming capability in rough surface ehl investigated using contact resistance. Submitted for publication. [32] T. Almqvist and R. Larsson. Comparison of reynolds and navier-stokes approaches for solving isothermal ehl line contacts. Proceedings of the WTC Conference, Vienna, Austria, 2001. [33] N. S. Wilkes I. R. Hawkins. Moving grid in harwell flow3d. AEA-InTec-0608, 1991. [34] CFX 4.4. User guide. AEA Technology, UK, 1999. [35] G. R. Higginson D. Dowson. Elasto-hydrodynamic lubrication: The fundamentals of roller or gear lubrication. Pergamon Press, Oxford, 1966. [36] T. Almqvist and R. Larsson. Some remarks on the validity of the reynolds equation for the modeling of lubricant flow on the roughness scale. Submitted for publication, ASME Journal of Tribology, 2004. [37] O. Reynolds. On the theory of lubrication and its application to mr. beauchamps towerŠs experiments, including an experimental determination of the viscosity of olive oil. Philos Trans R Soc, 177:157–234, 1886. [38] W. B. Rowe N. H. Woolley C. T. Schäfer, P. Giese. Elastohydrodynamically lubricated line contact based on the navier-stokes equations. Proceedings of the 26th Leeds-Lyon Symposium on Tribology, Leeds, pages 57–69, 1999. [39] C. H. Venner D. E. A. Odyck. Stokes flow in thin films. ASME Journal of Tribology, 125:121–134, 2003. [40] A. A. Lubrecht C. H. Venner. Transient analysis of surface features in an ehl line contact in the case of sliding. Transactions of the ASME, 116:86–93, 1994. [41] C. J. A. Roelands. Correlational aspects of the viscosity-temperature-pressure relationships of lubricating oils. Ph.D. Thesis, Technische Hogescool Delft, The Netherlands, (V. R. B.), Groningen, The Netherlands, 1966. [42] M. Peric J. H. Ferziger. Computational methods in fluid dynamics. Springer, Berlin, 1997. [43] D. Dowson P. Ehret and C. M. Taylor. Time-dependent solutions with waviness and asperities in ehl point contacts. Proceedings of the 29th Leeds-Lyon Symposium on Tribology, 1996. [44] C. H. Venner. Multilevel solution of the EHL line and point contact problems. Ph.D. Thesis, University of Twente, Enschede, The Netherlands, 1991. [45] C. H. Venner. Transient analysis of surface features in an ehl line contact in the case of sliding. ASME Journal of Tribology, 116:186–193, 1994. [46] M. A. J. Holmes. Transient analysis of the point contact elastohydrodynamic lubrication problem using coupled solution methods. Ph.D. Thesis, Cardiff University, year = 2002,.

116

BIBLIOGRAPHY

[47] R. Larsson T. Almqvist, A. Almqvist. A comparison between computational fluid dynamic and reynolds approaches for simulating transient ehl line contacts. ASME Tribology International, 37:61–69, 2004.