ElastoHydrodynamic Lubrication

Chapter 6 ElastoHydrodynamicLubrication In the final chapterof this bookthe differentcharacteristic elementsof the solversfor the two dimensionalPoissonproblem(Chapter3), the HydrodynamicLubricationproblem discussedin Chapter 4, and the Dry Contact problemdiscussedin Chapter 5 are all broughttogctherin a Multigrid solverfor the EHL circularcontactproblem.

6.1

Introduction

In the previouschapterthe dry contactbetween elastic bodieswas considered.From the resultsit may be clearthat even if sufficientlubricantis suppliedto form a film separating tile surfaces,tim shapeand thicknessof this film will be stronglyinfluencedby the deformationof the contactingelements.This situationis typical for lubricatedconcentrated contacts,where concentrated refersto the conditionthat the radii of curvatureof the two surfacesare of oppositesign (contraform)or of the samesign but with a large differencein value. Practicalexamplesare the contactsbetweenrollingelementand the inner and outerracewayin rollingelementbearings,tile contactsbetweengear teeth, and the contactsi)etweena cam and followerin a cam tappet mechanism. In the secondhalf of the twentieth centurythe stud)' of the film formationand the performance of thesecontactshas evolvedinto a separatedisciplinereferredto as ElastoHydrodynamicLubrication,abbreviatedto EHL. Characteristic for EHL contactsare the extremeconditionsof a very thin lubricantfilm, high pressure(up to 3.0 GPa for contacts betweensteelsurfaces)and high shearrates. Consequently, the elasticdeformationis not the only characteristicelement in EHL studies. Under these extreme conditions the lubricantbehaviourcan no longerbe characterized with a constantviscosityand density. Pressuredependenceof the viscosityand the densityand also non-Newtonianlubricant behaviourplay an importantrole in EHL contacts. For a detailed account of the historyof EHL the readeris referredto Dowson[35], [37] and Gohar [45]. Only in asymptoticand much simplified casesan estimateof a characteristic value suchas the centralor minimumfilm thicknesscould be obtainedanalytically. To obtainthe completepressureprofileand film shape,a numericalapproachwas needed; simultaneously solvingpressureand film thicknessfrom the Reynoldsequation,the film equationincludingelasticdeformation,and the force balancecondition. Classicalwork in this senseis the work of Dowsonand Higginson[34] (incorporated in [35]) for the line 179

180

CHAPTER 6. ELASTOHYDROD}WA.~IICLUBRICATION

contactproblemand of Hamrock andDowson[55] and [56] for the point contactproblem. The film thicknessformulaspresentedin thoseworks,which were obtainedfrom the numericalresults,are still widely usedtoday. Unfortunatelythe solutionalgorithmsfor the EHL problemturned out to be rather time consuming.Cpu timestendedto be O ( u 3) if N is the nun~berof nodeson the grid. This seriouslylimited the numberof nodestl~at could be usedand as a resultvery dense gridscould not be lltilised. A major step forward was formed by tl~e iI~troductioz~ of nmltigrid techniquesby Lubrecht[71, 72, 73] as an alternativemethodfor the nulnericalsolutionof EHL problems. Their furtherdevelopment is recordedin Venner[99] and Wijnant [110]. Theseworks have resultedin algorithmsthat enablethe solutionof line and point contactproblemsin an amountof work of O ( N l n ( N ) ) if -V is the total nunlberof grid points. This low work countfacilitatesthe useof densegrids,usingwork stationsand even personalcomputers. Moreover,the algorithmsare stablefor a wide rangeof loadconditions.The efficiencywith which tile basicequationscoul(l now t)e solvedcreatedroom for more complexproblems to be studied. The numericalsinmlationsfor theseproblemshave led to a significantly increasedunderstanding of mechanismsand phenomenathat can play a role in EHL contacts.Some of theseaspectswill be discussedin tl~eadvancedtechniquessectionat the end of this chapter. This chapterdescribeshow an efficientsolverfor the steadystateEHL circularcontact problemcan be obtainedby combiningthe lessonsof the previouschapters.The resulting algorithmformsan excellentbasisfor the studyof more complexEIIL problems.

6.2

Equations

For completeness the equationsdescribingthe EHL point contactproblemas given in Chapter1 are repeatedhere. The Reynoldsequationfor the steadystate casereads:

0 pha Op

0 pha Op

)+

)-

O(ph) = 0 "'"

(6.1)

ox

with p - 0 on the boundariesand the cavitationconditionp >_ 0 everywherein the domain. The viscosityrl is given by either the Barnsequation[5]" = ,j0 e x p ( o

(6.:)

where7/0 is the atmospheric viscosityand a is the pressureviscositycoefficient,or by tile empiricalrelationproposedby Roelands[93] T](p) - 7]oexp((ln(,]o)+ 9.67)(-1 + (1 + P)~)) Po

(6.3)

wherez is the pressu,'eviscosityindex, typicallyz - 0.6. and Po - 1.96. l0 s is a constant. The densityis either assumedto be constant,p - po, or it is given by the Dowsonand Higginson,'elation[35]:

p(p)-Po

5.9- 108 + 1.34p 5.9. l0 s + p

(6.4)

6.3. D I M E N S I O N L E S S E Q U A T I O N S

181

where p0 is the atmosl)heric densityand p is given in Pa. The equationfor the gapheight or film thicknesswas given as:

x2 h(x, y ) = h0 + + +

y2

2 ,+~r E'

[.+o0 /-

p(x', y') dx' dg'

where Rx and Ry are the reducedradii of curvatureof the surfacesin x and y direction respectively,and E' is the reducedmodulusof elasticity.For the circularcontactstudied in this chapterRv = Rx. The constanth0 is determinedby the force balancecondition: e+or e+r'~

p(x',y')dx'dy'

w= o --OO

(6.6)

o --O()

In physicalterms - h 0 representsthe mutual approachof two remotepointsin the solids. By remoteit is meant that the elasticdeformationat thesepointsis negligible. From theseequationsone can solvethe pressureand film thicknessas a functionof the operatingconditions.In particular,one is interestedin characteristic solutionparameters such as the minimumfilm thickness:

h,,, - rain h(x, y)

(6.7)

or tim central film thicknesswhich is defined as the film thicknessin the centre of the contact,i.e. at (x, y ) = (0, 0)"

= h(x = 0, v = o)

(6.8)

Alternatively,it can be defined as the film thicknessat the locationwhere (Op/Ox) = 0 anti (Op/Oy) = O.

h~ -- h(x, g)l(op/o~)=(op/ov)=o

(6.9)

For moderatelyto highly loadedcasesthereis no differencebetweenthesetwo definitions. In thosecasesthe pressureresembles the Hertziandry contactpressurewith the maximum occurringat the centreof the contact. However.for low load casesthe maximumof the pressureprofile occurson the centreline but in front of the point x - 0.

6.3

Dimensionless Equations

Assumingthat the Barus viscosity-pressure equationis used and that the lubricantbehaves incompressible, tile dimensionalequationsgiven above feature 8 parameters,urn, ~10,p0, (x, Rx. Ru, E' and w. This itnpliesthat the minimunifilm thicknessor the central filln thicknessarc'. afunctionof these8 parameters.However,theseparametersare not all independent.To reducethe numberof independentparameters,dimensionless variables can t)e intro(tuced.An ol)timal choiceis the set of dimensionless variablesthat reduces the numberof parametersto a minimum. This is for examplethe caseif the equationsare madedimensionless usingthe Hertziandry contactparameters, and the viscosityand density at ambientpressure.For the circularcontact(R, = By) the dimensionless variables are definedas:

182

C H A P T E R 6. E L A S T O H 3 " D R O D Y N A A I I CL U B R I C A T I O N

.X

~"

~x/o

P=

H-

P/Ph =

-

,1/,7o

v/a

hR~/a2 =

(6.10)

plao

Substitutionin Equation(6.1) gives: 0" OP" 0 " OP" ~-~ + ~-: ~ - = 0 OX ~ j ~ 0,~' /

O(fiH) OX

(6.11)

for X e [X~,Xb] and Y e [ - ~ , ~ ; ] where: tSHa with ~ _ 12umr/0R~ = (l~ aaph The boundaryconditionsare P(Xa, )') - P(Xb, ~') -- P ( X , - Y ~ ) = P ( X , ~ ) = 0. The cavitationconditionimposesP ( X , Y) >_ O. For an incompressible lubricant/~= 1. For a compressible lubricantone obtains: 5.9 10 9 s + 1.34phP 5.9. l0 s + p h P The Barusviscosity-pressure equation(6.2) in dimensionless form reads: p-

(6.12)

(6.13)

0 - exp(6P) with O = c~ph. The Roelandsequationin dinlensionless variablesis given by:

0 - exp((ln(,/0)+ 9.67)(-1 + (1 + Php)~)) (6.14) P0 Substitutionof the dimensionless variables in Equation(6.5) gives the followingdimensionlessfilm thicknessequation:

H ( X , I " ) - Ho + T + T +

(6.15)

-~

-~

~ v / ( X - x')~ + (~'- v,)~

with H0 deternlinedby the dimensiol~less force balancecolldition: Foo F~ ~

P(X, Y)dXd]"

2~r

(6.16)

3

If the Barus viscosity-pressure equationis used and tile lubricantis assumedto behave incompressible only two parametersremain.5 and ~. This impliesthat the dimensionless minimum film thicknessHm and the dimensionless central film thicknessHc (and any other characteristicsolutionvallle) will be a function of ~ and A only. This greatly reducesthe nulnberof computations to be I)erforlnediI~ parameterstudies.

183

6.4. DIMENSIONLESSP A R A M E T E R S

6.4

Dimensionless Parameters

In l l~e past,(lifl'erez~tsets()f (limelmionless l)arametershave been proposedto characterize EHL contacts. It does not really matter much which set is used, as long as it usesthe mini~numnumberand no redundantparametersare introduced. To presentthe resultsin designchartsor surveydiagrams,a convenientset of dimensionlessvariablesfor EHL problemswas proposedby Blok and co-workers,see [35]. For the line contactthey were pointedout by Moes [83] and for the point contactby Moes and Bosma[84]. Their derivationcan be foun(l in [86]. They are now often referredto as the ,~Ioesdimensionless parametersand for the circularcontactthey are definedas:

w

(71oU,) -3/1

(6.17)

and ( 710u~)1/4

L = a.E' E'R~

(6.18)

and the associated dimensionless film thicknessparameteris:

H^t = h ( -~

qoU~) -1/2 E'R~

(6.19)

Design charts are createdgiving the dimensionless minimum or central film thickness defined in this way as a functionof M and L. For EHL contactsdifferent regimesof asymptoticbehaviourcan be distinguished.In termsof ~I and L theseregimesare listed in Table 6.1. !regime

l ?~I [ L

rigid isoviscous small 0 rigid piezoviscous small large elasticisoviscous large O elasticpiezoviscous largc large Table 6.1" Regimesof EHL solutions. Tile parameters31. L anti H At are a convenientway to presenttile resultsin designcharts. However, from a numerical solution point of view the dimensionless equationsas given above arebetter suitedas all varial)lesre~nainrougillyill the range [0, 1]. Therefore,M and L will be usedin this work only to characterizea given load situationand to present designgraphs. Computedpressureprofilesand film shapeswill all be given in termsof the dimensionless variablesP and H as definedin tile previoussection. For convenience, the relationsbetween6, ,~ and H as the5'appearin the dimensionless equationsand 3I, L, and H M are given"

184

C H A P T E R 6. E L A S T O H Y D R O D Y N A M I C L U B R I C A T I O N

i' 1287r3'~1/3 =

3M 4 /

71"

HM-

H ~

(6.20)

At this pointit is noted that three parametersare often usedin the literature,to characterize a given contact.For the point contactproblemthey are referredto as the Hamrock and Dowsonparameters[55] definedas:

IV = U= G-

w 2 E IR x

rl~ E'Rz aE'

(6.21)

and the associateddimensionless film thicknessis defineda s H H D - h / R x . The disadvantageof usingtheseparametersis that it is not the minimumset of parameters.One reasonmentionedin favourof theseparameters is that they representparameters that one generallyvaries independentlyin experiments,i.e. load, speedand lubricant. However, this advantagedoes not outweighthe extra work neededin parametricstudies. Also, in termsof theseparametersthe differentasymptoticregimesare not so clearlydefined. Finally, the situationchangesif the lubricantis consideredcompressible and/or the Roelandspressureviscosityrelationis used. In additionto M, L or O, ~ two parameters of the group c~, z and 710have to be given whereas thethird parameteris determinedby the relation: c~po = In(r/o) + 9.67

(6.22)

z

whereP0 = 1.96.108 as mentionedbefore.

6.5

Discrete Equations

The equationsare discretizedon a uniformgrid with meshsize h in X and Y direction on a rectangulardomain Xa _ X <_ Xb, -)~,_< Y < ]'~. A secondorder accurate discretization of the Poiseuilletermsin Reynoldsequationis obtainedin the sameway as describedin Section4.2 for the hydrodynamicproblem. This givesthe followingdiscrete Reynoldsequationfor the point (i, j) of the grid:

185

6.5. D I S C R E T E E Q U A T I O N S

~~ih _

l/2,j i+lj +

p.h.

1/2 ',J-1 --

nh

h2

(~i,j-l/2 -[- ~ i , j + l / 2 )

ph.

~

',J -['- ~ i , j + l / 2

ph

i,3+1

h2 -o

(6.23)

with p.a ,,j - 0 for points on the boundary,and the cavitationcondition ph. > O. The coefficients~i~l/2,j and ~i,j+l/2 h are defined by:

h

(6.24)

with"

h )3 (,~,j = P(P;J)(H',i O(Rh. -

h

(6.25)

(pH)~ in Equation (6.23) denotesthe discrete'wedge'term If 9 Hh. ,o were a given function independentof the pressureany second order discretizationwould be suitable, e.g. a secondorder central discretizationsuch as'

(pH)~ " ~I~+'J H~+ 1,j - pa ',: H~_ ,,i (6.26) 2h where t~,• -/~(P/a+ld). However, in EHL this discretizationcan only be applied for low loads. As the centralterm for the t)oint (i, j) doesnot appearin the stencilwith increasing load, and thus dependenceof H on P, numerousiterative processeswill be unstable when applied to this discretization.Therefore,an upstream(one-sided)discretizationis preferred. A secondorder upstremndiscretizationis given by: z,h H h,-2j j + 0. 5 ~',-2,j (pH) .x ~ 1.opi,j'-hHh,,j - 2p~_~,jH~_l, h

(6.27)

This discretizationcan be used for all t)oints i > 2. At the first line near the boundary (i = 1) it can l)e replacedby a first order (liscretization" -h Hh - p~_~ .Hh ,o 'J ' - ' J1 (6.28) h Alternativelyone can computethe film thicknessin one extra point (i = - 1 ) and use the secondorder discretizationthroughout. As was shown in Section5.2, (5.14), tim film thicknessequationcan be discretizedas: ( -pH)x=Pi,J h-

Hh,,j -

H 0 +X~ V + V~]+ Z Z i'

K -hh i,i',)a' i'

(6.29)

186

CHAPTER

6.

ELASTOHYDRODYNAMIC LUBRICATION

where, for a secondorder discretizationbasedon an approximationof the pressureby a piecewiseconstantfunctionon a squareh x h around at X[,, 1~', the coefficientsIs are given by: hh

_

i,i',j,j, ~

2

.~j+h/2

dX'd1"'

",\'i+h/2

-#o

(6.30)

- v')'

v/(X, - +

and can be computedanalytically: !~ h h 2

/ ~i "p"~

~,,,,j,j, = 75{IX~,larcsi~h ~ +

() -[Xv[arcsinh '"' "-

-[.\'.,[a,csi,,h ~I -; "

(Xp~

Ilplarcsinh )~ I)p ]arcsinh I"

(v)

\

p/

II;.[arcsinh Xv

kXv/

{ '~m .

Xm

+lXm[arcsinh -~.,) +11 mlarcsinh( ~ ) }

(6.31)

with:

.Yp

=

X i - Xi, + h / 2

X,,,

=

Xi - Xr - h/2

);,,

-

}'j - I"j, - h / 2

(6.32)

On a uniform grid the coefficientsK,h.~I.j4, are a function of l i - i'[ and IJ - J ' ] only. Thus they can be pre-computedand storedin a single array of n~ x n u points. In the ~.t will be use(t, designatingK i,i',j,j' hh with [ i - i ' [ - k followingsections therefore also h "hh and {j - j'[ = I. Finally the discretizedforce balallce equationis derived, it is the stone equationas obtainedfor the dry contactproblem,(5.16): h2 ~ E ij

6.6

ph _ 27r i,j-

3

(6.33)

Model Problem

In order to designa stablerelaxationschemethat efficientlyreduceshigh frequencyerror components,it is necessaryto understandthe nature of the equations.Due to the exponential viscosity-pressure relation the coefficient~ in Eqllation (6.11) varies many orders of magnitudeover the domain. In the inlet region~ > > 1 whereasin the Hertzianconta(:t region~ < < 1. With the changeof { the characterof the problemchanges.This can be

187

6.6. MODEL PROBLEAI

illustratedby a simplemodel problemthat is obtainedwhen the force balanceequation is neglected,the coefficient{ is taken constant,,and f i - 1 is assumed: "02p

,{ ~ + ~

02p '~ Oy2t

OH =o OX

(6.34)

with: "+~ H ( X , Y ) - Ho + -X~2 + - ~~.2+ - ~ 2 f_+~ ~ o-~

P(X', Y') d X ' d Y '

~(X - X') 2 + ( Y - y,)2

(6.35)

for an arbitraryvalue of/4o. P = 0 is usedas boundarycondition.This model problem haslittle physicalmeaning,as is illustratedby the fact that the valueof H0 is irrelevantas it, cancelsin the equation.However,the problemis characteristic for the local behaviour of the EHL problem. For largevaluesof { the differentialaspectsas represented by the secondorder derivativesof the pressllredeterminethe behaviour.For smallvaluesof the OH/OX term dominatesthe behaviour,and, becausethe film thicknessis given by an integralequation,the problembehavesas an integralproblem. Understanding how to solve this model problemfor large and small valuesof ~ respectively,forms the key to understanding how to constructan efficient solver for the completeEHL problem. Consistentwith the discretizationof the fifll problem thediscretizedmodel problemis given by:

P/h_l,~ -2Ph. ,,: + P/h+l,~+ ~ pa. ,o-1 -2ph++Pa ij+l_ h2

h2

1. 5 u

i,j

_ 2u,

- l ,j

+ 0.su,

- 2,j

h

=0

(6.36)

with: Hh

H0 + X2

]7

l(hh i'

p.n

(6.37)

j'

*hh where the coefficientsI'~,,i,j,y are given by (6.31).

6.6.1 Relaxation

for Large ( Values

For large valuesof { Equation (6.34) is dominatedby the first two terms. Hence, it behavesessentiallythe sameas the Poisson2d problem.For this problemit was already shownthat Gauss-Seidel one point relaxationis an efficientsmootherwhen H h is given in all points. For tile presentproblmnthis relaxationcan t)e describedas follows. Let, /5./'~ denotethe m~rrentapproximationto ph, o an
---

p i.3h.

"-}-

w6~j

(6.38)

188

CHAPTER 6. ELASTOHYDRODYhL4MICLUBRICATION

to where w is the relaxationfactor and 6zhis the changethat shouldbe appliedto/3h. ,3 ~,3 solvethe equationat point (i, j) and can be written as (see Chapter2)" '4

,o

\

0P/' ,,3

(6.39)

]

with r.Z,3 h. the dynamicresidualdefinedas: + 7..h

--

,,j --

--~ h 2 _( z,j-1

~,3&1

1.Sfl?j ,-

+

h2

-" j + 0.SH _ ,j 2Hi_l, h

(6.40)

and O(Lhph),,j/OP~j the derivativeof the discreteoperatorworkingon ph at the point {i, j) with respectto p Z,3 h..

O(Lhph)iO_ = , OP~j ,,j

4(

1-5Koh,~ -- 2Khhl,O + 0 .0I~2,0 - -hh

h2

h

(6.41)

Next given the new approximationp!,. /~.h. z ,./ in each point, a new approximation z,3 to Hh. z,3 is computedin eachpoint usingEquation(6.37). Notice that it is not a full Gauss-Seidel relaxation.When relaxingthe pressure,tile changesalreadymade in pointsrelaxedpreviouslyare taken into accountin tile first two termsof the equation but not in the H!'. l, 3 terms, i.e. the discretewedge term izl the dynamicresidualdefined by Equation (6.40) only uses/}~ijvalueswhich were computedwith the old values/5.h.,o. Updating/7/.h,,jfor the changesappliedto points(Xi,))) l)reviouslyrelaxedis feasiblebut tendsto be computationallyexpensive,unlessit would besufficientto updateonly a few neighbouring points. Besides,as tile Multilevel Multi-integrationdescribedin the previouschapterprovidesa fast and efficientway to obtain all valueso f / } h aboveis z,3 at once, a schemeas described preferred. The actionof relaxationupon a Fouriercomponentof the solutioncan be analysedby meansof local mode analysis,see Section2.8. The behaviourdependson the parameter ( / h 2. For large valuesof ( / h 2 the error amplificationgraph is the same as given in Figure 2.9. The asymptotic smoothingfactoris fi = 1/2 for w - 1. However,with decreasing ~/h 2 the relaxationgraduallybecomesunstable. This is illustratedin Figure 6.1 which givesp(01,02) for ~/h 2 = 0.5. Tile figureshowsthat certainlow frequencycomponents are amplified. For this particularcase apl)lying tile changeswith an un(terrelaxation factor w = 0.5 will stabilizethe schemeat the expenseof a poor smoothingof high frequency components.For smallervaluesof ~/h 2 the amplificationof low frequencycomponents becomesstrongerand an increasinglysmallervalue of w is neededto stabilizethe scheme. Alternatively,the relaxationcan be applied in line relaxationmanner. Instead of scanningthe grid I)oint by point it is scannedline by line solving simultaneously the equationsof one line. Assumingline relaxationin the X directionfor a given line j the

189

6.6. M O D E L PROBLEI~I

1.5

-

1.2 0.9 0.6 0.3

5 1.0

o.

-0.5

0.6

02/n o.5

O1/n

Figure 6.1" Error amplification factor tL(Ol, 02) for Gauss-Seidelpointwiserelaxation applied to the model problemfor ( / h ~ - 0.5.

changesto be applied to/Sh, for 0 < i < nx of that line are simultaneouslysolved from a t,3 systemof equations,see Section 2.10.2: .4J~_~= r~

(6.42)

__r~a vector with the current residualsr z.h. Both where _6~ is a vector of changes6h,,j and ~h ,3 " i k are defined by: are vectorsof nx elements. The matrix coefficientsAi,

A~i,k= O ( L h p---h)i'J

(6.43)

OP~j

for 0 < k <

n~ and 0 < i <

n~. This gives: "hh

A~,k = _ for [ i - k ] >

"hh

"hh

l ' 5 h l i - k l ' ~ 2hli-k-ll'~ + 0"5hli-k-21'~ h

(6.44)

1 and

A~, i

A~,i_

1

--

-

--

A~,i+ 1 - -

2u i + o.5u h2 h2 h2

h 1 . 5 I ( ~, - 2/( 0hh , 0 + 0.5K~ho , h, 1.5Khh -hh + 0"5h3,o -hh t , o - 21(2,o h

if i > 1 if i < n ~ - 1

(6.45)

190

C H A P T E R 6. ELASTOH~"DROD'I%~AAIlC L U B R I C A T I O N

A~,~_1 and Ai,i+ J 1 are zero for i - 1 and i - n ~ - 1 respectively. Next the systemof Equations(6.42) shouldbe solved. For the two dimensionalPoissonproblemdiscussed in Chapter2 the systemof equationsto be solvedfor eachline was a tridiagonalsystemwhich can be solvedquickly. This would also apply here if H did not dependon P. However, due to this dependencethe matrix A j is a full matrix and its solutionwould normally requireO(n a) operations.Tile work count of one relaxationsweepwhich involvesn v - 1 lines would be O(nvna) and assumingn , - ny = v/N one relaxationwould cost O(N 2) operationsand thusbe very expensive.However, to obtain the line relaxationefficiency it, is generallynot neededto solve tile systemexactly. In practiceit is often sufficientto take into accountonly the terms in the summatioI~srelated to the direct neighboursof .-h0 decreaseswith increasingdistance.In particular a point i. This is justified because I~k, it is sufficientto solve a hexadiagonalsystemdefining all ,4~,k = 0 for k < i - 3 and k > i + 2. This hexadiagonalsystemcan then be solved quickly in O(n~) operations usingGaussianelimination.In that casetile work count of the relaxationis O(N) as for a pointwisescheme. After the systenlof equationsfor tile line j is solved tilechanges5~,j are sinmltaneously appliedto all pointsof the line with somerelaxationfactor wg~"

/sh,,j_/5h,,3+ wg~5~j

(6.46)

Subsequently, the systemof equationsfor the next line is constructedand solved. This processis repeatedfor each line of tile grid. Finally using the new approximation/3h a new approximation/}h to H h c a n be computed. The behaviourof the relaxationcan againbe analysedwith Local Nlode Analysis.Tile conclusionis the same as for the pointwiserelaxation. For large ~/h 2 it is stable with excellentsmoothing:# = 0.44. \Vith decreasing~/h 2 the processbecomesunstabledue to the amplificationof low frequencycomponents.The processcan be stabilizedusing underrelaxation but not for the small valuesof ~/h 2 that may occurin EHL contacts.In fact, the small valuesof co neededwould completelyremovethe smoothingpropertiesof the scheme. Summarizing,in this sectionGauss-Seidelpoint- and line-relaxationapplied to the model problemwere explained.In both casestile stabilityof tile relaxationis determined by the ratio ~/h,2. For large valuesof ~/h 2 both schen~esare stable and provide good smoothing.Hence, a nmltigridsolverllsing eitller of theseschemeswill yield an efficient solver. Note that this is only true if { / h 2 remainslargeon all the gridsthat are used. With decreasing{ / h 2 both schemesbecomeunstabledlle to the amplificationof smoottl error components.The amplificationof low frequencycomponentsis related to the changes occurringin H/h,j as a result of the cllangesapplied to all pressurevalues. As Hh.z,j is defined by a summationover all pressures,it will accumulatethe changesmade to the pressuresin the entire field. One can accountfor the changesby using updatedvalues of/7/h in the dynalnicresidual i.e lll)(tated for changesw9s5h. ,o applied to points (i j) already relaxed. This does hell) to stabilizetile scheme,however, the improvementis only marginal. Even conlbinedwith llI,(hurelaxationit does not result in a relaxation schemethat is stable and has acceptai)lestnoothingbehaviourfor tile small valuesof { / h 2 anticipatedto occurin EHL contacts.

191

6.6. klODEL PROBLEM 6.6.2 Relaxation for Small

s rV a l u e s

In the previmlssectionit was explainedthat tile Galiss-Sei(lelpoint- and line-relaxation are unstablefor small valuesof {/h 2. This instabilityis dim to the accumulation of changes in the summationinvolved in the definitionof H h To l)e precise: a change61,3 .h- applied -hh in point (i, j) will causea changeto /:/h. that is l)roportional t o /(i,i',j,j' As K o( 1 / r t,,3~ with r being the distancebetweenthe points (i,j) and (i',j'), the changein/7/O,j, will be proportionalto 6h,,a/r"The changein/7/0,a, as a resultof the entire relaxationsweepwill be tile summationof all such changes.It is this accumulationoccurringwhen re-computing tile film thicknesses, that causesthe instabilityof the entire process.Chapter5 showsthat the cure for this kind of instabilityis distributiverelaxation.When relaxinga point (i,j) changesare not only appliedto tile point (i, j) but also to a few of its neighboursin such a way that tile local Equation(6.34) is solved,as usual, but that the resultingchangesin ,o) are essentiallylocal. Using such a distributiverelaxation the integral (andthus in H.h may indeed takecare of the stabilityproblem. However,it still does not provideefficient smoothing.For very small ( valuesEquation(6.34) reducesto: 9

OH

~. 0 (6.47) 0X which is an equationin X directiononly. The only couplingbetweenequationsin tile Y direction is the coupling in the elastic deforlllationintegrals. Consequently,as for tlw anisotropic2D Poissollproi)leni (Secti(~ll 2.10.2). errors that are smooth in the X directionand oscillatoryin the t" directioncan notl)e efficientlyreducedby the relaxation. However,good smoothingcan be obtainedif lille relaxationis applied. A distributiveJacobiline relaxationis describedas follows. Given an approximation /3h and the associatedvaluesof/7/h for eachline j changes6-h. to be applieddistributively, are solvedfrom:

.4J~_~= r~9~

(6.48)

Both are vectors whereO"h _ is a vectorof challges6h,4and Ei .h. a vectorwith the resi(luaisr l,-h. 3 9 .j of n r - 1 elements.The nlatrix coefficients.4,. kare definedby:

Ai'k=

Ophj

1 lO(Lhp___h)i,j + O(Lhph)io OPkh- ,+, , 0Ph+~,,

4~

O(LhPh)i'J+ O(Lhph)io~ (6.49) OR" OP~j k,,+~ ,-1 J

for 0 < k < nx and 0 < i < nz. For details regardingthe systemof equationsthe reader is retbrredto Apt)endix C. This systemof eqllationsis constructedand solved for each line j. Subsequently,these changesare applied distributivelyin the same way as for the l)ointwisedistril)utiverelaxational)t~lied to the dry contactproblemdiscussedin the previouschapter The new a l)proxinlalioliat /Sh 1,3 is given I)v: " 9

p h - ~h

-h -h + a,,,_~ + ~h - ( 0 , + , j + 0,_,,, ,,,+,)/4);

(6.50)

With respectto the actual implementationof the relaxationit is noted that again the matrix A j is a fllll matrix but, as is shown ill Appendix C, it can be truncatedto a hexadiagonalsystemwhich call be solvedquickly as mentionedbefore.

C H A P T E R 6. E L A S T O H Y D R O D Y N A h H CL U B R I C A T I O N

192

The action of this line relaxationon the error can be analysedusing Local Mode Analysis.Figure 6.2 showsthe error amplification factor I~(01,02). The figure showsthat for small 0t, 02 the relaxationis stable. Moreover, it efficientlyreduceshigh frequency components.The asymptoticsnmothingrate p = 0.41. The relaxationis also stablefor large valuesof ~/h 2, but it is mucll lessefficientthan the Gauss-Seidel line relaxationin termsof smoothingl)ehaviour.

g _

0.80.60.4

~

1.0

0.2

-o.s 01/~

0.6

5

o.s

1.0

021n

.o

Figure 6.2: Error amplification factor p(Ol, 02) for distributiveJacobiline relaxationapplied to the modelproblemfor ( / h 2 = 0.0.

6.6.3

Relaxationfor Varying ~ Values

Neitherof the relaxationprocesses describedin the previoussectionis suitedto efficiently solvethe modelprol)lemfor all valuesof ~. The Gauss-Seidel lille relaxationis an excellent smootherfor large~/h 2 but is unstablefor small( / h 2. Oil the other hand tile distributive line relaxationis a good smootherfor very small ~/h 2 bllt rapidly losesefficiencywith increasing( / h 2. The objectiveis to obtain an efficientmultigridsolver. In that caseeven though~ is a constantthe ratio ( / h 2 will dependon the grid that is used. Therefore,a good choiceto obtain an efficientcoarsegrid correctioncycleseen~sto be to use different relaxationson differentgrids. 9On grids where{/h 2 > {ti,,,it the Gauss-Seidel line relaxationis used. 9Oil grids where{ / h 2 < ~limit tile .Jacobir

line relaxationis used.

193

6.7. R E L A X A T I O N OF THE EHL PROBLEI"~I

Where ~limit is a switch parameterto be defined. From practicaltestsit was found that an efficientcycle can be obtainedusing~limit~0.3 and someunderrelaxation in both the processes.Typicallywi~ =0.6 and wg~=0.8. In the EHL problemthe coefficient~ variesover the grid. The next step towardsa relaxationis thereforeto considerthe followingmodel problem:

0I

OP"

0I

OP"

o-~ ~b-~ +b-V" ,~N =,

OH

ox =0

(6.5x)

with P = 0 on the boundariesand H given by Equation (6.5) for different functions ((X, Y). Note that the constantH0 is still irrelevantfor this problem. To approximate the behaviourfor the full circularcontactproblem~(X, l') can be chosenas: if X 2 + ~-2 > f(x, )) - { (x +,))~,0 otherwise.

1;

(6.52)

When discretizedon a uniformgrid one obtainsexactly thediscrete Reynolds equation (6.23) if/5 - 1 is substituted.A suitablerelaxationprocessfor this modelproblemis now obtainedby combining thetwo relaxationson a singlegrid. This is justifiedby the fact that relaxation,as fat"as its effect on the error is concerned,is a local process. Tile combinedrelaxationschemecan be describedas follows. For a given approximation /Sh and associatedapproximation/~/h a new approximation/Sh is obtainedby scanningtile grid line by line, for each linej constructing a systemof equationsof the form (6.48). The matrixcoefficientsfor row i are definedby Equation(6.43) for the GaussSei(lel line relaxationschemeif the local value of ( h / h 2 > (u,m' and by Equation(6.49) z,3 for the distrii)utiverelaxationotherwise. For further details the reader is referred to AppendixC. The resultingsystemof equationsis truncatedto a hexadiagonal systemand solved. Sul)sequentlv, the change~ih is addedto the currentapproximation if at the point (i, j) the coefficients~/h 2 > ~ti,-,,itand it is also distributedto its four direct neighboursif ( / h 2 <_ ~,,,,it. Each type of changecan be (and will be) appliedwith a certain underrelaxation factor. Usingthe samerelaxationfa('torsas for tile problemwith constant~ this type of hybrid relaxationyielr an etfi('ient mllltilevelsolver fi)rEqllation (6.51) with (6.52) and (6.15).

6.7

R e l a x a t i o n of t h e E H L P r o b l e m

In the previoussectionsa stablerelaxationprocesswith good smoothingpropertiesfor Equation(6.23) with t5 = 1 was developed.In this equationthe variableH is definedby Equation(6.29) and ~ is a knownfunctionof X al~(tI'. The next stepstowardsan efficient r('laxati 0 and to include the relaxationof the force balance equation. In ad(lition OH/OX sl~o~ll(ll)e replacedby O(DH)/OX. Tl~is last changeis trivial and requiresonly someminorchangesin the coefficientsof the systemof equations that is solvedfor each line. see AppendixC.

194

CHAPTER 6. ELASTOH]'DROD}'NAMICLUBRICATION

The non-linearityis simply introducedby replacingEquation(6.52) with (6.25). However, when computingthe changeto be applied at each point the variation of { is not taken into accountin the relaxation. As a result the changeshave to be applied using s o m e underrelaxation. Cavitationwas treated in the hydrodynamicprot)lem (Chapter 4) as follows; after a change6h. is applied thenew pressure/~h is checked. If its value is smaller than zero it z ,3 1,3 is set to zero. In this way when relaxing the nextneighbour,cavitationof the previous point is taken into account. However. for a line relaxationprocess,the situationis a bit, more complicatedas all changesfor a line are solved simultaneously.The simplestway to incorporatecavit.ationwould l)e to apply all changesto a given line and set negative valuesto zero when needed. However.as the changesof the line are solvedsimultaneously this implies that within a single line the changefor a point i is computedas if a point i - 1 is not cavitatedat. all. As a result convergence of the relaxationmay stall due to the fact that a few pointsaroundthe cavitationboundaryconti~mouslyswitch back and forth betweencavitatedand pressurized.The alternativeis to create the systemof equations for a line taking into accountcavitation. Wheneverone of the neighboursof a point i has zero pressure,only the central term of the local system is used. As a result the change for this point is computedas it would be done for a pointwiseJacobi relaxation. Next the tbrce balanceequationis addedto the system.How to relax the force balance equationwas explainedin Section5.-1. After a ll,lml~er of relaxationson a singlegrid, the presentvalue of the variable/7/ois updatedaccordingto: /]0 = t7/0 + CUno(l~'fh -h 2 Y~./Shi,j)

(6.53)

i,)

where Wuo is some small constallt allr Iv f h is the right hand side of the force balance equation. For a single grid ~.fh = 2rr/3. The better the first, approximationto H0 the quicker the convergence.However, the nature of tile conw;rgenceof the force balance equationis determinedby the constantwho too. If this value is too large and the changes occur too often the solution can not digest the cl~angeof H0 and the residual of tt, e force balanceequationwill oscillate,whicll will also hamperconvergenceof the Reynolds equation. If Who is too small, convergencewill be too slow. In general the convergence behaviourwill be a dampedoscillationaround zero. The importanceof the force balanceequationshould not be underestimated.In the model problem with constant~, the integrationconstantH0 cancelsfrom the problem and therebythe force balanceequation. This implies that for the EHL problemthe val,le of H0 is in fact determinedby the boundarylayer that connectsthe region of large and small ~. Hence, the value of H0 is stronglycoupledto the coefficients~ in this region. The developedrelaxationcan llow be al~t)lier to the solutionof tim EHL problenlon a singlegrid. As an illustratiot~its convergenceI~el~avio~rfor two load conditionsis illustrated below. Table 6.2 specifiestl~e load casesin terms of differentsets of dimensionless parameters.Table 6.3 gives valuesof the parametersfor a contactbetweensteel surfac~.'s that would result in theseconditions.Also given are the maximum Hertzian pressureand the Hertzian contactradius. Figure 6.3 shows the error reductionas a fi~nction of the nund)er of relaxationsfor case 1. Resultsare presentedfor four grids with decreasingmeshsize coveringthe domain

6.7. R E L A X A T I O N OF THE EHL PROBLE,~I

195

I parameter!i case 1 lease 2 31 20 200 L 10 10 c3 9.89 21.31 ,\ 0.20 9 . 3 5 . 1 0 -2 IV 1.73.10 -7 1.73.10 -6 U 8.85-10 -12 8.85" 10 -12 G 4972 4972 Table 6.2: l hlues of the differentsetsof dimensionless parameterscharacterizingthe two load casesconsidered,and of the Hertzian dry contact parameters.

I parameterI a 7/0 Rx E' tts w Ph a

V~h,~ I Dimensi~ ] case 1 [ case 2 2 . 2 . 1 0 -8 [ N - l m 2] 4 0 . 0 . 1 0 -3 [Nm-2s] 1.6.10 -2 [m] 2.26" 1011 [N-Ira 2] 1.6 [ms-1] 10.0 100.0 [N] 0.45" 109 0.97" 109 [Nm -2] 1.02-10 -4 2.20.10 -4 [m]

Table 6.3: Values of the parametersfor the contact betweensteel surfacesfor the two load cases.

196

C H A P T E R 6. E L A S T O H ~ ' D R O D Y N A M I CL U B R I C A T I O N

X E [-4.5, 1.5] and t" C [-3.0.3.0]. The coarsestgrid for which resultsare presentedhas (16 + 1) x (16 + 1) points. Tim llext finer grid (32 + 1) x (32 + 1) points etc. For ea(:h casethe Hertzian I)ressureprofile was taken as a first approximationto ph and the first approximation to H h was obtainedfrom Equation(6.29) usingthe first approximationto H0 that is given as input. For the presentcase H0 = - 0 . 5 3 was used. The value of H0 was changedaccordingto Equation(6.53) each tenth relaxationwith a value WHo = 0.1.

I

10-2

I

1

!

10 -2

10-4

10-4

rn

r

10-6

10-6 "

"%

10-8 -

-

"%

,, 0

200

16'16

"~ 400

600 S

~"

32*32 10-8 64*64 128"128 " " " 800

1000

0

",,,, N,,

,,

200 400 600 800

32.*32 64.*64-128"128"'"

-

1000

S

Figure 6.3" Residualnorm (rn) of the Reynoldsequation (left) and absolutevalue of tile residual (r) of the force balance equation (right) on different gridsas a function of the number of relaxations(s). :'ll -- 20, L - 10, incompressible lubricant, Barus viscositypressureequation. Ttle left graph in Figure 6.3 showsthe evolutionof tile residualnorm of the Reynolds equation.Tile graph oil the right showstile evolutionof the absolutevalue of the residual of the force balancecondition.The figureshowsthat the relaxationindeedconvergesbut with decreasing meshsize the convergence speeddecreases.This is exactly the behaviour one would expectand it is causedby the fact that the low frequencyerror components are slow to converge.Notice the ripple on the Reynoldsresidualmlrve. This ripple is caused by the changes appliedto H0 when relaxing the force balanceequation. The residual of this equationdiminishesin a similar way as the residualof the Reynoldsequation. This impliesthat both equationshave rollghly the sameconvergespeed. In the first few relaxationsthe decrease is i~otnmnotonouslv an(t the residualof the forcebalance equation may increasein th(, first nmnberof iterations.This behaviollris causedby the choiceof the Hertzian pressuredistributionas an initial pressuredistribution.Tt~is apl)roximation exactlysatisfiesthe force balalmecondition. As an illustrationthe Figures6.4 and 6..5 show the computedapproximations to the dimensionless pressureph and film thicknessH h on a grid with n ~ - n y - 32. Although the resolutionis very crudethe sol~ltioncontainsall characteristic aspectsof the solution

197

6.7. R E L A X A T I O N O F T H E E H L PROBLEI~I

1

f

0.5 -

/~

2~

.~

~~~ ~

F i g u r e 6.4: C o m p u t e d dimensionless pressureph as a function o f X a n d ~t" on a grid

~vith 7~ --- nu -- 32. l~l - 20, L -- 10, incompressiblelubricant, B a r u s viscosity-pressure equation.

9

,,L

-

0.5

' -0~

5

-

2

15

1

F i g u r e 6.5: C o m p u t e ddimensionless film thicknessH j' as a function o f X and ~" oil a grid

with nx - n y -- 32. ~ I - 20, L -- 10, incompressible lubricant, Barus viscosity-pressure equation.

198

C H A P T E R 6. E L A S T O H Y D R O D Y N A M I CL U B R I C A T I O N

of an EHL point contact problem. The pressureresemblesthe Hertzian semi-elliptical pressurebut the singularityat X 2 + ~-2= 1 is replacedby a gradual pressureincreasein the inlet region and a rapid drop to the cavitationpressurein the outlet. Just beforethe cavitatedregion a 'pressurespike' occurs.This spike is clue to the exponentialviscositypressureequation: it does not occur for the isoviscouscase L - 0. Furthermore,in solutionsfor small values of L it only appearsas a local maximum in the pressure. As will be shown later, the spike actually forms a 'shield' wrapped around the solution in the outlet region. However, this can only be seen on a sufficientlyfine grid. This 'spike region' is the two dimensionalequivalentof the pressurespike first shown by Petrusevich for the line contact problem. Figure 6.5 showsthe computedfilm thickness. For convenienceof interpretationthe film thicknesshas been plotted upsidedown. In this way it can be seen as the shapeof the ball for the caseof a contactbetweena ball and a rigid fiat surface. Due to the elastic deformationthe surfaceis flattened in the central region. This effect decreasestowards the outside of the contact as a result of wl~ich a horseshoeshape occurs. The overall minimum film thicknessoccursin the 'side.-lobes'of the horse-shoe.The flatteningin the central region is in accordancewith the asymptoticbehaviourof the Reynoldsequation. For small ~ this equationreducesto:

O(pH) -

OX

~0

(6.54)

with the solution ~H = c(l'). This behaviourwill be more clearly visible in a contour plot of the solutionon a much finer grid, which will be discussedlater. The characteristiceffectsdisplayedin thesefigureswill becomestrongerwith increasing load as is shown by the results obtained for case 2. For this case the domain of computationis chosensmaller: X E [-2.5, 1.5] and I~" E [-2, 2]. Becauseof the higher load a smaller domain will be sufficientto accuratelypredict the film thicknessvalues. Figure 6.6 showstile evolutionof t l~e residualof the Reynoldsequationand the absolute value of the residualof the force balanceequationa.s a function of the number of relaxations. The coarsestgrid for which resultsare shown has (32 + 1) • (32 + 1) points. The resultsfor the grid with (16 + 1) • (16 + 1) points have been omitted as this grid does not yield a physicallyacceptablesolution. The convergedresult yields a minimum fihn thicknesswhich is smaller than zero. From a mathematicaland numericalpoint of view this is not a problem. However, such a solutionhas no physicallymeaning. This bizarre result indicatesthat on the (16 + 1) • (16 + 1) grid, the discretization error is locally larger than the solution. Figure 6.6 showsthat also for the higher load case, where the distributiverelaxation is certainly used in the contact centre, the relaxationconvergesin the same way as for the previouscase. The residual of the force balanceequationconvergences along with the residualof the Reynoldsequationexcept for the first few iterations. The Figures6.7 and 6.8 show the approximatesolutionsto the pressureP and the film thicknessH as a function of X and Y on a grid with n~ = ny = 32. The effects shown in the Figures 6.4-6.5 have indeed becomestronger. The pressuredistributionmore closely resembles the semi-ellipticHertzian pressureprofile with the singularityat X 2 + y2 = 1 replacedby a continuouspressureincreasein a small boundarylayer. The pressurespike appearsto

6.7. R E L A X A T I O N OF THE EHL PROBLEAI

I

I

I

199

I

10-2

I

I

i

\ )i

I

I

10-2 10-4

10-4 rll

r

10-6

10 -6

10-8 ,,i 0

200

\

10-8

64*64 128"128~

400 600 800

1000

0

200

S

64*64 128"128

400

600

800

1000

S

F i g u r e 6.6: Residual norm(rn)of the Reynoldsequation (left) and absolutevalue of the residual (r) of the force balance equation (right) on differentgrids as"a function of the numberof relaxations(sl. ~I -= 200. L - 10, incompressible lubricant, B a r u s viscositypress~weequation.

I-

'//.-..

X

15

-2

F i g u r e 6.7: Computeddimensionless pressureph as a function of X and Y on a grid with nx = l~y - 32. 1~I - 200, L -- 10, incompressible lubricant, B a r u s viscosity-pressure

equation.

200

C H A P T E R 6. E L A S T O H Y D R O D Y N A M I CL U B R I C A T I O N

H 05 \ 1-

15

~\/ -2.5 -2 -1~

l!

.,,2~1 152

~

,~.2 15 Figure 6.8: Computeddimensionless fihn thicknessH h as a fi~nctionof X and }" on a grid with nx = n u = 32. ~I - 200 L - 10, incompressible lubricant, Barus viscosity-pressure equation.

be absentnow, however,this is not true. With increasingload, the width of the pressure spike expressedin the dimensionless coordinateX decreases and thereforeit will only be visible on a suffici(;ntlyfine grid. The film thicknessresemblesthe dry contact result. In that case H would be zero in the unit disc X 2 + I .2 < 1. For the EHL case the asymptoticsolutionsatisfiesEquation (6.54) which givesa small film thicknessthat is constantas a functionof X but still variesas a functionof Y. Finally, it is noted that the side lobeswhere the minimumfilln thickness occurs have decreasedin (relative)size, and have movedto the outsideof the contactdisc X 2 + ~-2 < 1. Consequently, to resolve the minimumfilm thicknessacc~lratelyfor high loads,rather fine grids are needed. Summarizing,it is concludedthat a suitablerelaxationprocesshas been developed from the point of viev,,of stability. However, its speedof convergence is far too low to efficientlyobtain solutionson densegrids. This problemwill be overcomeby the use of Multigrid techniquesto be describedin the followingsections.

6.8

Coarse Grid Correction Cycle

Having developeda stable relaxationprocessfor the non-linearEHL problem,the next step is to introducecoarsergrids and ~Iultigrid techniques to overcomethe slow convergence offine grid relaxation. The constructionof a coarsegrid correctiollcycle for the EHL problemcan be done in the sameway as for the dry contactproblem,see Section5.4. The FAS schemeshould be applied to all equationsas was (loire for the dry contact t)roblem ill the previous chapter, and once more, specialattentioll sllouhl be given to the intergrid transfersin the cavitatedregion and the treatmentof the force balanceequation. To obtain good

6.8. C O A R S E G R I D C O R R E C T I O N C ) ' C L E

201

asvlnptoticconvergencethe transferof the residualsin the region near the cavitation boundaryis very important.The best approachis to use full weightingand for all points 9 j, (l', ) :/: ( i , j ) appearingin the stencildefine the residualto be zero if the point has zero pressure.This is a slightlymore advancedapproachthan simply injectingthe residual. The developedcoarsegrid correctioncycle provesto be an efficient tool to solve the EHL problemfor a wide range of load conditions.However. its performancewill depend on the choiceof the relaxationfactorsand on the accuracyof the first approximationof H0. The relaxationfactorsshould be small enoughto avoid instability. For most cases good performancecan be obtainedtaking ~timit= 0.3, 0.2 <_ wj, <_ 0.6, 0.4 < wg~ < 0.8, and WHo ~ 0.1. Also, in many casesW-cycles give better performancethan V-cycles. The performanceof the developedcoarsegrid correctioncycle for the two load conditions consideredin the previoussectionis outlinedbelow. W-cycle 1 2 3 4 5 6 7 8 9 l0 20

level 3 64 x 64

level 4 128 x 128

level 5 256 x 256

6 . 4 3 . 1 0-i 1.39.10-1 2 . 7 0 . 1 0-2 4 . 0 3 . 1 0-a 7 . 4 4 . 1 0-~ 1.48. l0 -4 3.17. l0 -5 7 . 5 2 . 1 0-6 1.99.10-6 5 . 6 1 . 1 0-7 6 . 3 7 . 1 0-19

7 . 3 0 . 1 0-l 2 . 1 8 . 1 0-I 7 . 0 1 . 1 0-2 2 . 2 4 . 1 0-2 6 . 2 1 . 1 0-a 8 . 8 7 . 1 0-4 1.87-10-4 4 . 0 8 . 1 0-5 9 . 6 8 . 1 0-6 2 . 3 4 . 1 0-6 6 . 2 2 . 1 0-12

9 . 0 2 . 1 0-1 3 . 9 0 . 1 0-1 2 . 1 3 . 1 0-I 1 . 7 9 . 1 0-l 1 . 0 6 . 1 0-I 7.88. l0 -2 3 . 2 8 . 1 0-2 1 . 1 5 . 1 0-2 2 . 3 6 . 1 0-3 4 . 9 7 . 1 0-4 6 . 6 2 . 1 0 -l~

Table 6.4" R e s i d u a ln o r m o f the R e y n o l d sequation on three differentgrids as a function o f the n u m b e r o f W(2, 1) cvcles. 31 - 20. L - 10, incompressiblelubricant, Barus viscosity-pressure equation. R e l a x a t i o nparameters:wia - 0.4, Wg~ - 0.6 and WHo -- 0.1. The Tables 6.4 and 6.5 show the residual of the Reynoldsequation and tile force balanceequationas a function of tile number of II'(2, 1) cycles for three grids for the case 31 = 20 and L - 10. Tile domain was chosena s - 4 . 5 < X < 1.5 and - 3 . 0 < ~'" _< 3.0. and the coarsestgrid (level 1) used in the cycle has (16 + 1) • (16 + 1) points. The Hertzian pressureprofile was taken as an initial approximationfor ph. The first approximationfor H h was obtainedby straightforwardcomputationaccordingto Equation(6.29) with an initial approximation to H0 (H0 = -0.53). Throughoutthe cycle ~Iultilevel ~illlti-integrationis used for tile fast evaluationof the elastic deformation integrals. The algorithmas describedill the previouschapteris generallyused with 6 th order transferswhich ensuresaccurateevaluationfor all casespresentedin this chapter. The results were further obtained using underrelaxationfactors wia = 0.4 for the Jacobichangesand Was = 0.6 for the Gmlss-Sei(lelchanges. The relaxationfactor used in the relaxationof the force balanceconditionwas chosenWHo - 0.1. These valuesare r suitablefor low and moderatelyloaded cases. However, using W-cycles with increasing ~

C H A P T E R 6. E L A S T O H ] ' D R O D ] ' N A M I CL U B R I C A T I O N

202

W-cycle 1 2 3 4 5 6 7 8 9 10 20

level 3 64 • 64

level 4 128 x 128

level 5 256 • 256

2 . 6 2 . 1 0 -2 3 . 0 3 . 1 0 -3 2 . 5 4 . 1 0 -4 1.78. 10 -l 9.11. 10 -5 2 . 6 7 - 1 0 -`5 8.58. 10 - 6 2.83 910-6 9 . 3 0 - 1 0 -7 2 . 8 9 . 1 0 -~ 1.23-10 -12

2 . 7 6 . 1 0 -2 7 . 2 9 . 1 0 -3 1.89-10 -3 2.27- 10 -4 3 . 1 2 - 1 0 -5 1 . 0 9 . 1 0 -7 3.83. 10 -7 1.31 910-7 1 . 3 9 . 1 0 -7 3 . 9 6 . 1 0 -8 7 . 8 7 . 1 0 -13

6 . 9 1 . 1 0 -3 1 . 0 2 . 1 0 -2 7.82. 10 - 4 1.65 10 9 -3 2 . 9 5 . 1 0 -'i 1 . 3 3 . 1 0 -4 8 . 1 3 . 1 0 -5 2 . 0 7 . 1 0 -5 1.25-10 - 6 8 . 2 2 . 1 0 -~ 2 . 3 9 . 1 0 -12

Table 6.5" Residualnorm of the force balanceon three differentgrids as a functionof the numberof W(2, 1) cycles. J~l = 20. L = 10. incompressible lubricant, Barus viscositypressureequation. Relaxationparameterswi~ = 0.4, w~ --- 0.6, ~,'~Io= 0.1.

numberof gridlevelsthe numberof visits to the coarsergrids also increases.Consequently, on very fine grids, a smallervalue of Wtfo may give better results. This problemdoes not occur when V-cycles are used, but in that case the asymptoticperformanceis less good. Table 6.4 and 6.5 show that the cycle convergesand asymptoticallythe error reduction per cycle is roughly a factor of 3. This appliesboth to the residualof the Reynolds equationand to the residualof the force balanceequation. Note t h a t this is the asymptotic factor. Performancewill often be better in tile first few cvclesand on coarsergrids. Also larger underrelaxationfactors can be used. which, on the coarsergrids gives better performance. However, this is not necessarilytrue on very fine grids and for many cycles. If the underrelaxationfactors are chosentoo large this may lead to a situation where the cavitationboundaryis movingback and forth which resultsin poor asymptotic convergence. As the conditionsof the problemvary significantlywith varyingp a r a m e t e rvalues,ttle relaxationparameterswhich give optimal perfornlancevary as well. The generaltendency is that. with increasingload more under relaxationis neededand WHo should be smaller. A safe choice for many casesis a~'j, - 0.2 and Wgs = 0.4. T h e force balance constant shouldbe taken WUo < 0.05 on very fine grids fi)r W-cycleswhereasit, can be taken O(0.1) for V-cycles. This conservativechoice leads to a less than optimal performanceon the coarsergrids where a larger factor could be used. Nevertheless,as will be shown in tile next section,this performanceis sufficientlygood to obtain a solutionwith an error t h a t is below the discretizationerror in only a few cycles. The resultsfor load casenumber2 with theserelaxationfactorsare given in Tables6.6 and 6.7. These tables show that botl~ ill terms of the residualof the Reynolds equation and in terms of the residualof the force balanceequationthese parametersgive a stable convergence with an asymptoticerror reductionof approximately2 per cycle.

6.9. FULL M U L T I G R I D W-cycle 1 2 3 4 5 6 7 8 9 10 20

203 level 3 64•

level 4 128• 128

level 5 2 5 6 x 256

3 . 9 7 . 1 0-l 1.43-10-l 4 . 9 3 . 1 0-2 1.71.10-3 6 . 0 0 . 1 0-3 2 . 1 5 . 1 0-'3 8 . 0 2 . 1 0-4 3 . 1 8 . 1 0-4 1.38.10-4 6 . 7 7 . 1 0-5 6 . 2 8 . 1 0-7

4 . 3 1 . 1 0-l 1.64.10-l 6.23-10-2 2 . 2 4 . 1 0-2 8 . 1 7 . 1 0-a 3 . 1 0 . 1 0-a 1.21.10-3 4 . 8 9 . 1 0-4 2 . 0 6 . 1 0-4 9 . 0 6 . 1 0-5 1.93.10-7

4 . 7 6 - 1 0-l 2 . 1 0 . 1 0-l 1.03-10-1 6 . 0 3 . 1 0-2 2 . 9 0 . 1 0-2 1 . 4 6 . 1 0-2 6 . 0 0 . 1 0-3 2 . 3 3 . 1 0-a 9 . 5 1 . 1 0-4 4 . 0 1 . 1 0-4 2 . 4 7 . 1 0-7

Table 6.6" Residualnorm of the Reynoldsequationon three differentgrids as a function of the number of IV(2,1) cycles. 31 = 200, L - 10, incompressible lubricant, Barus viscosity-pressure equation. Relaxation parameters wj~ = 0.2, wg~ = 0.4, and WHo = 0.05. -2.5
6.9

Full M u l t i G r i d

Tile algorithmfor tile EHL problemis completedwith the extensionto a Full MultiGrid (FMG) algorithm. This is a very straightforward step. Insteadof only using the coarser gridsfor convergence acceleration,they flllfill an additionaltask in generatingan accurate first approximationon the grid. With respectto the solutionof the problem for high load conditions(large M, large L) it is important to separatethese two functionsof tile coarsegrid. With increasingload tile discretizationerror on a given grid increases. If the grid is ratller coarsethe solutionof the discreteproblem on this grid will be a poor approximationto the solutionof tile differentialproblem. Consequently,it will be a bad initial approximationto tile solutionon the next filler grid. Characteristicfor situationswhere tile discretization error is large is tile occurrenceof a significantnegative film thicknessvalues in tile solution. In combinationwith tile non-linearityand the cavitationconditionthesenegativevaluesof the film thicknesscan seriouslyfrustratethe convergence of the relaxationprocessand even causeit to diverge. Hence, the first grid used in the FMG algorithmmust have a sufficientlysmall mesh size. In very extreme cases tilecoarsestgrid in the cycle shouldalso be finer, but this is not often the case. In tile correctioncycle the coarsergrids serve to generatea smoothcorrectionto the fine grid problem,basedon residualinformationfrom the more accuratefine grids. Summarizing,for high loads tile FMG processshould not be started with simple relaxationson the coarsestgrid. Rather one shouldstart with a few correctioncycleson a sufficientlyfine grid. Subsequently, the FMG processcan be continuedas usual. In Chapters4 and 5 it was shownthat usinga FMG algorithmwith one or two cycles per level yieldeda solutionwith an errorsmall comparedto the discretization error. Below it will be shownthat this also appliesto the EHL problem. Table 6.8 showsthe residualof tile Reynoldsequationfor four grids as a functionof

204

C H A P T E R 6. E L A S T O H 3 " D R O D Y N A M I C L U B R I C A T I O N W-cycle 1 2 3 4 5 6 7 8 9 10 20

level 3 64 x 64

level 4 128 x 128

level 5 256 x 256

1 . 0 5 - 1 0 .2 8 . 2 8 . 1 0 .3 5 . 6 1 - 1 0 -a 3 . 1 7 . 1 0 .3 1.68. 10 - 3 8 . 6 9 . 1 0 -4 4 . 3 8 - 1 0 -4 2 . 1 9 . 1 0 -4 1 . 0 8 . 1 0 -4 5 . 2 6 . 1 0 -5 3 . 1 5 . 1 0 -8

2 . 2 9 - 1 0 -2 2 . 7 6 . 1 0 .3 4 . 4 2 . 1 0 -4 3 . 5 9 . 1 0 -5 5.00" 10 -5 2 . 5 9 . 1 0 .5 8.59. 10 - 6 1 . 6 9 . 1 0 -6 5 . 8 2 . 1 0 -7 1 . 0 7 . 1 0 -6 1 . 5 4 . 1 0 -s

1 . 0 6 . 1 0 -2 1 . 3 0 - 1 0 .3 1.54. 10 - 3 4 . 7 2 . 1 0 -4 1.25" 10 -4 4 . 1 3 - 1 0 -5 1 . 5 5 . 1 0 -'5 6 . 0 4 - 1 0 -6 1 . 9 0 . 1 0 -6 4 . 2 1 - 1 0 -7 6 . 7 9 - 1 0 -9

Table 6.7: Residual norm of the force balance equation on three different grids as a function o f the n u m b e r of II'(2, 1) cycles. ~I - 200, L = 10, incompressiblelubricant, B a r u s viscosity-pressure equation. Relaxation parameters: wja - 0.2, wg.~ - 0.4, and W H o - 0 . 0 5 , - 2 . 5 <_ X _< 1 . 5 , - 2 < V _< 2.

the n u m b e r of W ( 2 , 1) cycles carried out in a Full l~lultigrid a l g o r i t h m for load case 1. From this table it can be seen t h a t when using a FXIG l)rocesswith a single IV(2, 1) cycle per level, the n o r m of the residualof the Reynoldsequationfor the final solutionon level 5 is a l r e a d y smallert h a n the residualnorm for the solutiont h a t was o b t a i n e dapplying6 cycless t a r t i n gwith the Hertzian a p p r o x i m a t i o nas an initial solution,see Table 6.4. T h e work investedto obtain the FI~IG solutionis only slightlylarger t h a n the work of a single cycle. Consequently,using the convergedsolutionof a coarsergrid as the initial solution on a finer grid, really pays off. To d e t e r m i n ehow m a n y cycles are neede(t to obtain a solution with an error of the m a g n i t u d eof the discretizationerror Tables 6.9-6.10 can be used. These tables list the central and the m i n i m u m film thicknessfor the solutionspresentedin Table 6.8. Tim W-cycle 1 2 3 4 5 10

level 3 64 x 64

level 4 128 x 128

level 5 level 6 256x256 512x512

2 . 3 9 . 1 0 -2 5 . 3 1 - 1 0 -3 1.32 10 9 -3 3.67.10-t 1.04 1() 9 -t 1.44. lO -7

1.45 10 -2 1.48 10 -3

7 . 4 5 . 1 0 -3 7 . 6 8 - 1 0 -4 2 . 5 3 - 1 0 -4 7.38" 10 -5 2 . 5 0 . 1 0 -5 1 . 5 2 - 1 0 -7

2.82 10 -'l 5.89 10-5 1.28

10 -5

1.62

10 -~

5 . 3 1 . 1 0 -3 1 . 0 2 . 1 0-3 6.37" 10 -4 3.14" 10 - 4 2 . 2 8 . 1 0 -4 4 . 1 6 . 1 0 -5

Table 6.8: Residualnorm of the Reynoldsequation on four differentgrids as a function of the n u m b e r of IV(2, 1) cycles used in a F•IG algorithm. 11I - 20, L -- 10, incompressible lubricant, Barus viscosity-pressure equation. Relaxationparameters:wja - 0.4, wgs - 0.6,

and WHo -- 0.1.

205

6.9. FULL MULTIGRID W-cycle 1 2 3 4 5 l0

level 3 64 • 64

level 4 128 x 128

level 5 256 x 256

4.7940.10-1 4.9336.10-1 4.9644-10-1 4.8185.10-1 4.9358.10-1 4.9655-10-l 4.8190.10-~ 4.9357.10-l 4.9652.10-1 4.8188.10-l 4.9357.10-l 4.9653-10-l 4.8187.10-1 4.9357.10-l 4.9653.10-1 4.8187.10-l 4.9357.10-1 4.9653.10-l

level 6 512 x 512 4.9739.10-1 4.9724-10-1 4 . 9 7 2 9 . 1 0-1 4 . 9 7 2 7 . 1 0-1 4 . 9 7 2 8 . 1 0-l 4.9727-10-1

Table 6.9: Central tilm thicknessdetinedas H h at (X, Y) = (0, 0) on four differentgrids as a functionof the numberof II'(2, 1) cyclesin a FMG algorithm. M = 20, L = 10, incompressible lubricant,Barus viscosity-pressure equation. W-cycle 1 2 3 4 5 10

level 3 64 x 64

level 4 128 • 128

level 5 256 • 256

level 6 512 x 512

2.9531.10-1 2.9826.10-1 2.9863.10-l 2.9874-10-~ 2.9877.10-1 2.9878.10-~

3.0399-10-1 3.0509.10-t 3.0515.10-1 3.0516.10-l 3.0516.10-1 3.0516.10-1

3.0605. l0 -1 3.0631.10-i 3.0629.10-1 3.0630.10-1 3.0630.10-l 3.0630.10-1

3 . 0 6 6 1 . 1 0- 1 3 . 0 6 5 2 . 1 0-1 3.0657-10-1 3 . 0 6 5 5 . 1 0-1 3 . 0 6 5 5 . 1 0-1 3 . 0 6 5 5 . 1 0-1

Table 6.10: I~linimumfihn thicknesson three differentlevelsas a functionof the number of IV(2, 1) cyclesin a FMG algorithm. M = 20, L - 10, incompressible lubricant,Barus

viscosity-pressure equation. resultsshow that alreadyusingtwo |t"(2, l) cyclesper gridlevelan approximationto the minimumand centralfilm thicknessis obtainedwith an error that is very small compared to the discretizationerror. This can be seenfrom a comparisonof the differencebetween the valueafter 2 cyclesand 10 cycleswith the differencebetweenthe valueson the different grids. Finally, the order of the discretizationerror itself can be checkedcomparingthe valuesobtainedon the differentgrids. Each time the mesh size is refined the changein the centralfilm thicknessbecomesapproximatel)" 4 timessmallerwhich is consistentwith the fact that a secondorder accuratediscrizationhas been used. To concludethe resultsfor this load case,Figures6.9 and 6.10 show the dimensionless pressureand film thicknessfor the level 6 solutionas a functionof X and Y. Comparing Figure 6.9 with Figure 6.4 showsthat the higher resolutionresultsin a much more pronouncedpressurespike. Furthermore,this spike forms a shieldof high pressurevalues precedingthe drop to the cavitationpressurep h _ O. Figure 6.11 (left) showsa contour plot of the film thickness.This figureclearlyshowsthe horse-shoe shapeof the film thicknessand the flat region in the centre. In accordancewith Equation (6.54) for constantp the contourlines tend to be straightlines. An alternativeway to presentthe film thickness resultsis by means of a pseudointerferometry graph. Insteadof drawingcontourlines, all intensityZ is definedaccording to"

206

CHAPTER 6. E L A S T O H Y D R O D Y N A M I CLUBRICATION

F i g u r e 6.11" Contourplot (left, & H = 0.01) and pseudo-interferometry plot (right, A H = 0.05) of the computeddimensionless film thicknessH h as a functionof X a n d ~'. M = 20, lubricant,B a r u s viscosity-pressure equation. L = 10, incompressible

207

6.9. FULL M U L T I G R I D

Z(X, Y) - Z0 + 20 cos

AH

/

(6.55)

where Z = 1 gives a squareh~ • h~ of white color, Z = 0 gives a black square,and squareswith grey shadesare obtainedfor 0 < Z < 1. Figure 6.11 (right) shows sucha pseudo-interferometry graph of the film thicknessfor load case 1. This type of graph is interestingsinceit allowsa direct comparisonwith the film thicknessmapsobtainedwith a popularmeasurement techniquereferredto as 'Optical Interferometry'. This technique was developedin the 1960's at Imperial College, London, by Foord et al. [44] and has presentlymaturedto a level where the film thicknessmapscan be obtainedfor transient problems,and also very thin fihns, e.g. see Cann [29] and Kaneta [67, 68]. In an optical interferometryrig, a steel ball is run againsta glassdisk that is coveredwith a semi reflectivelayer (and a spacerlayer for the thin film version). Light is shoneinto the contact. Part of it is reflectedby the layer, and part of it goes through the film and is reflectedlater by the ball. The result is an interferencepattern representinga film thicknessmap of the contact. Tables6.11-6.13 show the resultsobtainedusinga FMG algorithmfor load case2. In termsof the residualnorm of the Reynoldsequation,startingwith the convergedsolution of a coarsergrid on level 5 savesthe work of about 5 cycles. Furthermore,the valuesof the centraland the minimumfilm thicknessgivenin Table 6.12 and 6.13 can be compared. A comparisonof the valuesafter many cycleswith thoseobtainedafter 2 cycles shows that alreadyafter 2 cyclesthe error in the centraland minimumfilm thicknessis small comparedto the discretizationerror. Finally, these valuesconvergein a secondorder mannerwith decreasing meshsize. The differencebetweenthe valuesof the centraland minimumfilm thicknessobtainedon two consequentgrids decreaseswith a factor four, when decreasing the meshsize. Finally, the Figures6.12 and 6.13 showthe computeddimensionless pressure and film thicknessfor the level 6 solutionas a functionof X and Y. ComparingFigure 6.12 with Figure 6.7 showsthat becauseof the higher resolutionthe pressurespike is now clearly visiblein the solution.Figure6.14 showsa contourplotand a pseudo-interferometry graph of the dimensionless film thickness.ComparingFigure 6.14 with Figure 6.11 showsthat for the higher load case the almost flat region in the center with straightcontourlines coversalmostthe entire Hertzian dry contactregion. Consequentlythe side-lobeswhere the minimumfilm thicknessoccurshave becomevery small. The Full Multigrid algorithmcompletesa fast and efficient solver for the EHL circular contact problem. The computingtime requiredby the solver is proportionalto O ( N In(N)) if N is the numberof nodeson the grid. With the algorithmpresented, the circularcontactproblemcan be solvedefficientlyfor a wide rangeof load conditions.The programin its standardform as given in AppendixH will run many load caseswithout any problem. However, asthe EHL problemis a non-linearproblemin which different effectsplay a strongeror weaker role with variationof the parametervalues,it is very difficult to give one set of relaxationparametersand domainsizesthat ensuresgoodconvergeneefor all cases.The programcan be an efficientsolverfor many casesif it is used with somecare. This care appliesto the choiceof the size of the domainwhich must be sufficientlylarge on the inlet side, in particularfor easeswith small M values. On the

208

C H A P T E R 6. E L A S T O H Y D R O D Y N A M I CLUBRICATION

Figure 6.14: Contourplot (left, A H - 2.5. 10-3) and pseudo-interferometry plot (right, A H = 2 . 5 . 1 0 -2) of the computeddimensionless film thicknessH h as a functionof X and Y. M = 200, L -- 10, incompressible lubricant,Barus viscosity-pressure equation.

209

6.9. FULL MULTIGRID W-cycle 1 2 3 4 5 10

level 3 64x64

level 4 128• 128

level 5 level 6 256x256 512x512

9 . 7 2 . 1 0 -2 2 . 1 8 . 1 0 -2 9 . 7 8 . 1 0 -a 4 . 9 0 . 1 0 -3 2 . 5 9 . 1 0 -a 1.22-10 -4

5 . 8 9 . 1 0 -2 1 . 0 7 . 1 0 -2 3 . 2 3 . 1 0 -a 1 . 2 4 . 1 0 -a 5 . 1 3 . 1 0 -4 1 . 4 2 . 1 0 -s

3 . 6 3 . 1 0 -2 3 . 4 3 . 1 0 -3 9 . 8 9 . 1 0 -4 4 . 0 3 - 1 0 -4 1.86-10 -4 9 . 3 3 . 1 0 -6

3 . 3 0 . 1 0 -2 1.35.10-3 5 . 4 8 . 1 0 -4 2 . 1 4 . 1 0 -4 1.08- 10-4 6 . 1 2 . 1 0 -6

Table 6.11" Residualnorm of the Reynoldson four differentgridsas a functionof the number of IV(2, 1) cyclesin a FMG algorithm. M = 200, L = 10, incompressible lubricant,

Barus viscosity-pressure equation.

W-cycle 1 2 3 4 5 10

level 3 64 • 64

level 4 128 x 128

level 5 256 x 256

level 6 512 x 512

7.4704-10-2 8.3089-10-2 8 . 5 6 6 5 . 1 0-2 8 . 6 3 2 1 . 1 0-2 8 . 6 5 1 2 . 1 0-2 8 . 6 6 8 2 . 1 0-2

9 . 4 5 9 5 . 1 0-2 9 . 6 1 9 0 . 1 0-2 9 . 6 2 8 6 . 1 0-2 9 . 6 3 5 1 . 1 0-2 9.6377-lO -2 9 . 6 3 8 7 . 1 0-2

9.8932- 10-2 9 . 8 8 7 9 . 1 0-2 9.8895- 10-2 9 . 8 9 0 1 . 1 0-2 9.8902-10-2 9 . 8 9 0 2 . 1 0-2

9 . 9 4 7 3 . 1 0-2 9 . 9 5 4 2 . 1 0-2 9 . 9 5 3 8 . 1 0-2 9 . 9 5 3 9 . 1 0-2 9 . 9 5 3 8 - 1 0-2 9 . 9 5 3 8 . 1 0-2

Table 6.12: Central film thicknessdefinedas H h at (X, Y) = (0, 0) on four differentgrids as a functionof the numberof W(2, 1) cyclesin a FMG algorithm. M - 200, L - 10,

incompressible lubricant,Barus viscosity-pressure equation.

other hand for large M values, i.e. M > 500 care shouldbe taken t h a t the F M G process is startedon a sufficientlyfine grid, see Figure 5.10. Finally, with increasingvalue of the L parameter,the problembecomes more difficult to solve. In that case the piezoviscouseffects becomevery strong andthe film becomes very stiff. This implies t h a t the processis very sensitive tothe changesof H0 and to the quality of the first approximationon a given grid. In fact it becomesvery sensitiveto any change. Changesin the positionof the cavitationb o u n d a r ycondition,or interpolation errors in the coarsegrid correctionmay render the algorithm unstablefor the extreme casesof L > 20. The use of a less extreme viscosity-pressure relation may make these problemseasierto solve. Generally,if stability problemsoccur for a particularload case variousmeasurescan be taken. For example, one should ensuret h a t tile first a p p r o x i m a t i o nto H0 is a good one. Secondly,the grid oll which the F M G processis started should be sufficientlyfine. Next, the coarsestgrid in the cycle shouldbe sufficientlyfine. This can often be ensured by choosinga smallerdomain for higher load cases(large M). And finally, the relaxation parameterscan be reduced.

210

C H A P T E R 6. E L . 4 S T O H Y D R O D Y N A M I C L U B R I C A T I O N W-cycle 1 2 3 4 5 10

level 3 64 • 64

level 4 128 • 128

level 5 256 x 256

level 6 512 • 512

3.0780.10-2 3.3918.10-2 3.4848.10-2 3.5373.10-2 3.5579-10-2 3.5701.10-2

3.8239.10-2 3.9675.10-2 4.0179.10-2 4.0302-10-2 4.0331.10-2 4.0342.10-2

4.0236-10-2 4.1386.10-2 4.1440- 10-2 4.1448.10-2 4.1450.10-2 4.1452.10-2

4.1414.10-2 4.1702.10-2 4.1695.10-2 4.1696.10-2 4.1695.10-2 4.1696-10-2

Table 6.13: Minimum film thicknesson four differentgrids as a functionof the numberof W(2, 1) cyclesin a F M G algorithm. M = 200, L = 10, incompressible lubricant, Barus viscosity-pressure equation.

6.10

Design Graphs

The programas describedcan be used to computethe film thicknessand pressuredistribution as a function of the contactoperatingconditions.Consequently,the evolutionof specific aspectsof the solution (pressureor film thickness) canbe analysedusing parametric studies. Obviously,the film thicknessitself is important, as, comparedwith an estimateof the roughnesson the surfaces,it can yield an indicationof the lubrication condition of the contact in real applications. The questionarises what film thickness value should be used, the minimum film thicknessbecauseit is the smallestvalue, or the central film thickness because it extendsover a much larger region? As an example a designgraph for the central film thicknessis given in Figure 6.15. This figure shows the central film thicknessdefined as the film thicknessat the point (0, 0) as a functionof M and L. Note that the dimensionless film thicknessaccordingto Equation (6.19) is used. In the chart the four asymptoticregimesmentionedin Table 6.1 are indicated. The markersindicatethe resultsof numericalsolutions.Also shownin the figure are the predictionsgiven by two film thicknessformulas. The dashedlines indicate the predictionsof the formula for the central film thicknesspresentedby Hamrock and Dowson [56]. This formula was obtainedby curve-fittingusing a set of numericalresults as input. As can be seen from tile figure, the Hamrock and Dowson equation yields accurateresultsfor M > 10 and 2.5 < L < 25. However,for small valuesof M and L the predictionyields valueswhich are much too low. An alternativeway to constructa film thicknessformulais to createa functionfit using solutionsfor the asymptoticregimesas buildingblocks. This appraochwas introducedby Moes, e.g. see [86, 87, 99]. It has the advantagethat it yields a formula that is applicable in the entire parameterrange,includingthe asymptoticregimes.For the line contactcase the solutionsfor the asymptoticregimescan be derivedanalytically.For the point contact problemapproximatesolutionsfor the asymptoticregimescan be usedas buildingblocks, see [87]. A formula to predict the central film thicknessfor ttle circular contact can be obtainedusing the followingapproximateasymptoticsolutions:

6.10. DESIGN GRAPHS .

.

211 .

.

.

.

.

i

.

.

.

.

.

.

.

.

!

+~o~?cou~ .

I0

HM

.

.

.

.

.

.

.

.

.

1

E~-p~~

k. 25

9 .~,~- ~ ............... '. . . . . . - . . . . + _ . _ _ : I 0 : ,so,,~o,.,s.~__~::::: ........... , ........... ,. 5

ci~ 1

2.5 1. Elastic4soviscous

0.1

.

.

1

.

.

.

.

.

1

10

.

.

.

.

.

.

.

,

.

.

.

.

.

1 o0 M

.

.

*0

J

1000

.

.

.

.

.

.

.

.

10000

Figure 6.15" Centralfilm thickness chart Hc Mas a functionof M and L. The drawn lines indicatethe predictionsof (6.61). The dashedlines indicatethe predictionsof Hamrock and Dowson [56].

9Rigid lsoviscous: H~{ = 41.4M -2

(6.56)

HrA~ = 0.91L 2/3

(6.57)

He~t= 2.42M -2/t5

(6.58)

9Rigid Piezoviscous:

9Elastic Isoviscous:

A more accurateequation based on a fit through numericallyobtained values for this case, see [99], is:

H~[= 1.961~1-1/9

(6.59)

H ~ = 1.25M-l/12L3/4

(6.60)

9Elastic Piezoviscous:

212

C H A P T E R 6. E L A S T O H u 1 6 9

LUBRICATION

Introducinga parameters to ensurea smoothtransitionbetweenthe asymptoticregimes, the followingfilm thicknessformulais proposed:

+ 0.1]-3/8] 2~/3-~-[(H~M)-8-'~-(HeMp)_8]_s/8 } 1/s (6.61) where: 3I6

s=-1+ 2

exp

H~ \

"

5H~z.

The drawn lines in Figure 6.15 have been obtained~sing Equation (6.61) with (6.56), (6.57), (6.59), and (6.60). As can be seenfrom this graph the formulaaccuratelypredicts the computedcentralfilm thicknessover the entire rangeof parametervalues.

6.11

Conclusion

In this chapterthe variousaspectsof the previoussolvershave been combinedinto an efficientsolver for the circular EHL problem. These includethe solversfor the Poisson problem,for the HydrodynamicLubricationproblem,for the Dry Contactproblemand in particularthe Multilevel Multi-integration.As eachstepin the EHL algorithmrequiresat most O ( N l n ( N ) ) o p e r a t i o n s , the total cpu time neededby the solveris also O(NIn(N)). This enablesthe solutionof the EHL problemon densegridsand thusmakes thealgorithm well suitedfor extensiveparameterstudies. In this chapteronly a few pressureand film thicknessresultswere presented. Once the pressureand film thicknessare known, the pressureand film thicknessdistributioncan be used as input to computemany other quantities. For exampleusing the multilevel multi-integrationthe sub-surfacestresses can be computedefficiently. The algorithm(:an be extendedin many ways to cater for more complex(realistic) EHL problems. For exanlpledifferent viscosity-pressure equationsand densitypressure equationscan be used. The solver can be extendedto elliptic contactsas described in [110]. Non-Newtonianlubricant behaviourcan also be incorporated,by means of effectiveviscosities,involvingonly minor changesto the program. A related extension involvestile 'slip-at-the-wall'as describedby Ehret et al. [38, 39]. More involved are the addition of fractionalfilm contentas a variableto study the effectsof starvationand variablelubricantsupply,see Chevalieret al. [31, 32]. The algorithmhas servedas a basis for a solver of the time dependentproblemof roughness movingthroughthe contact,seeVennerand Lubrecht[100]-[106]. Alternatively a transientanalysiswas used to obtain insight in the dynamicsof the contactsuch as stiffnessand damping,see Wijnant and Venner [109]-[112]. In generaltheseefficientnumericaltools allow one to generatehuge amountsof data in relativelyshort times. As an example: the transient2d EHL solver easily produces 2 • 104 • 512 • 512 - O(101~ data points. As the calculationsare performedin double precisionjust storingthe pressureand film thicknessdistributionat eachtime step,would requireO(100) GBytes. It is clearthat the old adagio"calculationis not aboutgenerating

6.12. ADVANCED TOPICS

213

numbersbut all aboutgeneratinginsight"is particularlyvalid for thesesolvers.In many casesthe vast amountof numbershas been reducedto simpleapproximativeequations that can be interpretedin physicaltermsand solvedon a pocketcalculator.

6.12

Advanced Topics

As mentionedabove,the low complexityof the algorithmcreatesroom for more complex EHL problemsto be studied.Someof theseproblemsare briefly describedbelow. 6.12.1 Roelands Equation, Compressible Lubricant The resultspresentedin this chapterwereobtainedassumingan incompressible lubricant and usingthe Barusviscosity-pressure equation.This wasdoneon purposeas it represents the simplestsituationand hencethe problemcan be describedby only two parameters. However, as was explainedin Chapter 1 the Barus equationis valid only for a limited range of pressures. With increasingpressuresthe predictedviscosityis generallytoo high. An improvedequationis the Roelandsequation(6.3). In addition,due to the high pressures occurringin EHL contacts,it is no longerjustifiedto assume the lubricantto be incompressible. The relationoften usedto modelthe increaseof the densitywith pressure is the Dowsonand Higginsonequation(6.4). The programgiven in Appendix H is fully written in terms of a variable density. It can be used to solve the pressureand film thicknessdistributionusing either of the viscosityequationsand/or assumingan incompressible or compressible lubricant. This sectionpresentssomeresultsobtainedfor a compressible lubricantand usingthe Roelands viscosity-pressure equation.From a numericalpointof view, addingtheseequationsmakes ttle problemeasierto solve. First, becausethe Roelandsequationpredictsa viscosity that increaseslessrapidlywith pressurethan accordingto the Barusequation.Secondly, becausethe compressibility addssomettexibilityto the system.As a result,the problem becomesless sensitiveto changesin the value of H0 when relaxing the force balance equation. The resultspresentedapplyto the sameload casesconsidered throughoutthis chapter. However, asmentionedin Section6.4 usingRoelands'equationand assuminga compressible lubricanttile problemis no longergovernedby two parameters.In additionto M and L, two other parametershave to t)e given. Tile resultspresentedhere have beenobtained usingc~ - 2.2.10 -s [Pa -1] and '10 = 4 0 . 1 0 -3 [Pas]. UsingEquation(6.22) this gives z = 0.67. Tables6.14 and 6.15 give the dimensionless minimumand centralfilm thicknessas a functionof the gridlevelobtainedusinga FMG algorithmwith 3 W(2,1) cyclesper level for the two load cases. The Figures6.16 and 6.18 show the film thicknesscontourplots and pseudo-interferometry graphs. Finally, Figures6.17 and 6.19 show the computed pressureat the centerlineY = 0 as a functionof X and at the line X = 0 as a function of Y. For referencethe resultsobtainedusingthe Barusviscosity-pressure equationand assumingan incompressible lubricantare shownin the latter two figures. From Tables6.14 and 6.15 it can be concludedthat the convergence behaviourof the

214

CHAPTER 6. E L A S T O H Y D R O D Y N A M I CLUBRICATION level n. • ny Hc 2 32• 3.8174.10-1 3 64 x 64 4.1904. 10 -1 4 128x 128 4.2872-10-l 5 2 5 6 x 256 4.3116.10-I 6 5 1 2 x 512 4.3177-10-l

Hm 2.6291-10-1 2.8622.10-1 2.9094.10-1 2.9218-10-1 2.9237.10-1

Table 6.14: Central film thicknessdefinedas H h at (X, Y) = (0, 0) and minimumfilm

thicknessas a functionof the gridlevelcomputedwith a FMG algorithmwith 3 W(2,1) cycles.M - 20, L = 10, compressible lubricant,Roelandsviscosity-pressure equation. level nx x ny 2 32x32 3 64 x 64 4 128 x 128 5 256 x 256 6 512x512

H~ 4.0557.10-2 7.0686- 10-2 7.8872. 10 -2 8.0935.10-2 8.1447.10-~

H~ 1.6921" 10-2 3.3080" 10-2 3 . 7 1 2 0 " 10 -2

3.8480.10-2 3.8760.10-2

Table 6.15: Central film thicknessdefinedas H h at ( X , ) ' ) - (0, 0) and minimumfilm thicknessas a functionof the gridlevelcomputedwith a FMG algorithmwith 3 W(2,1) cycles.M - 200, L = 10, compressible lubricant,Roelandsviscosity-pressure equation. minimum and centralfilm thicknesswith decreasingmeshsize displaystile secondorder discretization.Each time the mesh size is refined the changein the calculatedvalues decreasesby roughly a factor of 4. Comparingthe valuesof the central and minimum film thicknessgiven in Table 6.14 for load case 1 with thosegiven in Tables 6.9 and 6.10 it can be concludedthat the minimum film thicknesshas hardly changedwhereasthe centralfilm thicknessis now smallerby about 13%. The sameappliesto the secondload caseas can be seenfrom a comparisonof the valuesgiven in Table 6.15 with thosegiven in Tables6.12 and 6.13. For this load casea reductionof 18% of the centralfilm thickness is obtained. However, not only the value of the centralfilm thicknesshas changed. As can be seenfrom a comparisonof Figures6.16 and 6.18 with the Figures6.11 and 6.14 the shapeof the film in the centralregionhas also changed.This is more clearly observed in Figures 6.17 and 6.19. The originof the differentfilm thicknessshapescan be understood from the equations. As has been explainedbefore, due to the high viscosityand small film thickness,the Reynoldsequationin the centralregionof the contactreducesto:

O(pH) .~ 0 OX

(6.62)

with the solution for the incompressible case H = c(Y). The function c(Y) will be determinedby the flow into the contact. In principle it will depend on the viscositypressureequationthat is used. However,as in the inlet regionthe pressures are small,the

6.12. A D V A N C E D TOPICS

215

F i g u r e 6.16" Contourplot (left, A H = 0.01) and pseudo-interferometry plot (right, A H 0.05) of the computeddimensionless film thickness H a as a functionof X and Y. M - 20, lubricant,Roelandsviscosity-pressure equation. L --- 10, compressible

2

2

1.5-

1.5

I-

0.5

-2

-1.5

-1

~.5 X

0

0.5

"0 1

1.5

,

,

,

,

,

1.5-

IHPI

0.5 -

O' -2.5

2

2

-1.5

.

0.5

0j -1.5

.

.

.

.

...

/.-

-1

" ..................

-1

0.5

".......................

~.5

' 0

0.5

1

0 1.5

Y

Figure 6.17: Dimensionless pressureP and film thicknessH as a functionof X at the line Y = 0 (left) and as a functionof Y at the line X = 0 (right) for M = 20, L = 10. Drawn line: incompressible lubricant, B a r u s viscosity-pressure equation. Dashed line: compressible lubricant,Roelandsviscosity-pressure equation.

216

C H A P T E R 6. E L A S T O H Y D R O D Y N A M I CL U B R I C A T I O N

F i g u r e 6.18: Contourplot (left A H - 2.5- 10 -3) and pseudo-interferometry plot (right A H - 2 . 5 . 1 0 -2) of the computeddimensionless film thicknessH h as a functionof X a n d lubricant,Roelandsviscosity-pressure equation. Y. M - 200, L - 10, compressible

1.5

1.25

0.5

-

1,5

1 25

,

,

,

,

,

-

0.4

0.4

1-

1-

0.3

0.3

P

0.75 - H P

0.75 -

H

0.2 0.5

0.5

0.5-

-

0.1

0.25 -

0' -2.5

0,2

-2

-1.5

-1

-0.5

0

05

1

0 1.5

0.1

0 25 -

0 -1 5

-1 ..0.5 0

X

0.5

I

0 1.5

Y

F i g u r e 6.19: Dimensionless pressureP and film thicknessH as a functionof X at the line Y = 0 (left) and as a functionof Y at the line X = 0 (right) for M = 200, L = 10. lubricant. B a r u s viscosity-pressure equation. Dashed line: D r a w n line: Incompressible

compressible lubricant,Roelandsviscosity-pressure equation.

6.12. ADVANCED TOPICS

217

dimensionless viscositywill roughlybe unity. After all, for small valuesof the pressure the differencein the viscositypredictedby these equationsis small. It can therefore be concludedthat c()') is the same when using either the Roelands equationor the Barus equation. If c(I') is the same for both casesthis implies that the level of the film thicknesswill be the same for both casesand thus neither the minimum nor the central film thicknesswill changemuch if the Roelandsequationis used insteadof the Barus equation. Consequently,the changes observed in the central film thickness,and the centralregion of the contactmust be due to the compressibility.This is indeedthe case.

For the compressible case Equation(6.62) imposes/~H= ~(Y) with ~(Y) also determined by the inlet region of the contact. However, as the pressuresare low in the inlet region,compressibility effectsonly play a minor role as/~ ~ 1. Hence ~(Y) .~ c(Y). Thus the level of the product/~Hfor the compressible casewill roughlybe the sameas the level of H for the incompressible case. The consequences for the minimum and central film thicknesscan be derivedstraightforwardly.The minimumfilm thicknessoccurscloseto the contactperipheryX 2 + }.-2 ~ 1 wherethe pressures are small. Hence, p = 1 and thus the minimumfilm thicknesswill roughlybe the sameas in the incompressible case. For the central film thicknessthe situationis different. If ~H is the same for an incompressible and compressible lubricant,H can no longer be the same in the central region. In this regionthe pressures are high and/~ > 1. Along a line in the X direction,the relationpH = c(}') dictatesthat the productof the densityand the film thicknessremains constant. Consequently,changesin the densityhave to be compensated by changesin the film thickness.With a semi-ellipticalpressurevariationalong thisline, this implies that the film thicknessprofile will no longer be a constant,but will display a curved line with a local minimum at X - 0. This behaviourwill be similar for every line in the X direction,however,each line will have a differentconstant.Moreover,the density increasealongeach of theselines is different,as the pressuredecreases goingtowardsthe peripheryof the contact. Thesechangesin the film thicknessimply that the contourplot will no longershow straightlines in the centerof the contact. This is illustratedby the Figures 6.16 and 6.18 which showtile contourplots and pseudo-interferometry graphsof the film thicknessfor the solutionsobtainedass~lminga compressible lubricantand using tile Roelandspressureviscosityequation. It shouldbe noted that the changeof the centralfilm thickness due to the compressibility of the lubricantcan be predictedquite accurately.Becausethe conditionsin the inlet of the contactwill be the samefor an incompressible and compressible lubricantthe functionc(Y) will also be the same. As a first order approximationone can use:

H~=

1-

ff(P~) /

H~

(6.63)

if H~ is the central film thickness obtained for the compressiblecase,//~ its value for the incompressiblecase, and P~ the pressure in the centre of the contact. For the present cases Pc ~ 1 which gives f = 1.15 for load case 1 (Ph = 0.45 [GPa]) and/~ = 1.21 for load case 2 (Ph = 0.97 [GPa]). Consequently, the central film thickness obtained for the compressible case should be 15% lower than for the incompressiblecase, for load case 1. For load case 2 it should be about 20% lower. These values are indeed close to the variations obtained

218

CHAPTER 6. ELASTOHYDRODYNAMICLUBRICATION

numerically.This approachcan be usedto predict theeffectof compressibility on the film thicknessfor any densitypressureequation,e.g. see [103] whereresultsobtainedusingthe densitypressureequationproposedby Jacobson and Vinet [62] are presented. Finally, the pressuredistributionchangestoo. For load case 1 the compressibility has a significanteffect on the pressureprofile. This resultsin lower pressuresin the central regionand in particularin a much lower pressurespike. However,as the compressibility predictedby the Dowsonand Higginsonequationis limited to 34%, the differencesin the pressureprofile decreasewith increasingload. For load case 2, the only visibleeffect is the decreaseof the pressurespike which becomeslesspronounced.

6.12.2 Time Dependent Problems The EHL problemis generallya transientproblem. Only in the idealizedsituationof perfectlysmoothsurfacesand stationaryoperatingconditionsdoes the problemreduce to a steadystate problem. If the effect of surfaceroughnessor micro geometryare to be studied,the problembecomesinherentlytime dependent.Thus tile pressureand film thicknesshaveto be solvedas a functionof time with the microgeometrymovingthrough the contact.The variabletime addsanotherdimensionto the problem andaugmentsthe time requiredfor a numericalsolution.The stationarysolverpresentedin this chaptercan be extendedto includetransienteffects. For the time dependentproblemthe Reynolds equationis given by:

0 ph3Op)

0 ph3tOp

O(ph) O(ph) =0 (6.64) Ox Ot with p = 0 on the boundaryand tile cavitationconditionp(x, y, t) > O. The film thickness equationfor the generalcaseincludingsurfaceroughness can be written as:

x2 y2 h(x, y) = ho(t)+ - ~ + ~Ry

r(:r, y, t ) + ~

2 /§ /§ o0

p(x', y') dx' dy'

~ ~/(x - x') 2 + ( y - y,)2

(6.65)

where:

~(x,y,t) - ~ ( ~ - ~,,t,y) + ~ ( ~ - ~ t , y )

(6.66)

standsfor the undeformed roughness of the two surfaces,and U 1 and u2 are the velocitiesof the two surfaces.The integrationfuntionho(t) is in generaldeterminedby an equationof motion,see [110]. If the accelerationforcesare neglectedand a constantload is assumed, this equationreducesto the usualconditionof force balance.When both load, speedand geometryvary as a functionof time the sameapproachcan be used,see Messd [82].

w(t) =

oojf+oo

p(x', y', t) dx' dV

(6.67)

Using the dimensionless variablesdefinedby (6.10) and definingthe dimensionless time

as:

219

6.12. ADVANCED TOPICS

(6.68)

T : tum/a resultsin the followingdimensionless Reynoldsequation"

0I

OP ~ 0 i +

OX ~ - ~ j

OP ~

O(pH)

~ OYj -

OX

a(pH) = 0 OT

(6.69)

where: - pH3= with ,X = 12u,,,rloR~

FIA

a3ph

The boundaryconditionsare P(Xa, Y ) = P(Xb, Y ) = P ( X , - Y , , ) = P(X, Y , ) = O. The cavitationconditionimposesP(X, ]q T) > O. The dimensionless film thicknessequation is given by: X2

].-2

H(X, l, T) - Ho(T) + --~ + ~ -

P(X', Y', T) dX'dY'

2 f _ +~= / _ oo ~ ~ ( x , ~, 7) + -~ v/(X _ X,)2 + ( y _ y,)2 (6.70)

with" TO(X, Y, T) = 7r

- u~/umT,Y) + TC2(X - u2/umT,Y)

(6.71)

and H0 is determinedby the dimensionless force balancecondition:

f _ + ~ f + ~ P ( X }, T) dX dY = -21r O0

O0

~

(6.72)

3

The introductionof the roughnessmovingthroughthe contactadds one additionalparameterto the problem,i.e. only the parameteru~/um needsto be specified.If ul/u~ is given then by definitionu2/um followsfrom ul + u2 = 2urn. The objectiveis now to solveP(X, }, T) and H(X, Y, T) as a functionof time. In the sameway as was done for the steadystateproblem,the natureof the transientproblem can be analysedby lookingat the limit for very highviscosities.For the transientproblem one obtains:

_O(pH) _ O(pH) = 0 (6.73) OX OT whilst for the steadystate problemt~H = c(Y) is obtainedas explainedbefore. For the caseof surfaceroughness the stationaryconditionwill imply that insidethe high viscosity regionall roughnesswill be flattened. For the transientproblemthis equationimposes pH = p H ( X - T ) . Thus, the productpH will tend to be a functionof X - T i n the high viscosityregion. This is a naturalconsequence of demandingmassconservation in a flow betweenflexiblewalls wherethe viscoustermscan be neglected.Note that this will happenregardless of the dimensionless surfacevelocitiesUl/Um and u2/um. This implies that any changein the film thicknessinducedby roughness uponenteringa high viscosity regionwill be propagatedthroughthe high viscosityregionat the dimensionless speedof

220

C H A P T E R 6. E L A S T O H Y D R O D Y N A M I C L U B R I C A T I O N

unity. In physicalterms this is the averagevelocityof the surfaces.This behaviourleads to very interestingand sometimesconfusingaspectsin transientsolutionsof the problem which at any time will consistof a combinationof components 9Only for pure rolling the situationreducesto a simplesituationas in this case the velocityof the surfacesequals the averagevelocity. For more detailsregardingthe transientaspectsof the problemthe readeris referredto publicationssuch as [76, 77], [99]-[112] and the referencestherein. The problemcan be discretizedin exactly thesameway as the steadystate problem. In additionto the grid in spacea 'mesh'in time can be assumedwith meshsize (timestep) hT. Subsequently, the discreteReynoldsequationfor each timestepcan be derivedas:

-

z+l/2,j,k

h2 ~,j phh ~h ph -1/~,k i,j-l,k -- (~i,j-~/2,k + ,,:+~/2,k) h2

+

h ph. ~,j,k+ ~i,j+l/2,k ,,j+~,k (pg)x-h

(pH)~

-0

(6.74)

with the coefficientsdefinedby:

•

-

(~,,:,k+ ~i,j+l,k)/2

(6.75)

The discrete'wedge'term can be taken as -h h ( p H ) h : 1 5pij,kHi,i, k - 2 f h_' ,j,k H ih- I 9

,j,k -'[- 0 . 5 ~ h 2,j,k H~-2,j,k

(6.76)

-h H h -h H h Hh (pH)~._.l'5pi,:,k i,j,k--2p,,j,k-I i,j,k-1+0"5fi~,j,k-2 i,j,k-2

(6.77)

x--

h

and the discrete'squeeze'term a s

hy

At this point it is noted that an accurateapproximationof these terms is very important. The wedgeand squeezeterm togetherform a so-calledadvectionoperator,i.e. they describea propagationmechanismin the solutionalong a characteristicX - T. If the discretizationerror is too large this will lead to artificial amplitudedecay effectsin the solution. One should at least choosea discreteapproximationto the combinedwedge and squeezeterm such that it has a zero discretization error for the characteristic direction. If the meshsize and timestepare chosenequal the standardsecondorder upstreaIn discretizationgiven abovesatisfiesthis requirement.A more advancedoption is to use a combined discretization of the wedgeand squeezeterm as is commonfor advectiveterms in ComputationalFluid Dynamics,see [107]. Having obtainedthe discretesystemof equationsthe next step is to solve the equations. For this purposethe multigrid solver as presentedin this chaptercan be used. The extensioninvolvesincorporationof the discretewedgeterm in the discrete Reynolds equationand in the definitionof the systemsof equationsto be solvedfor each line in

6.12. ADVANCED TOPICS

221

the line relaxation. Subsequently, the coarsegrid correctioncycle can be used for each timestep,usingthe solutionof a previoustimestepas a first approximation.Alternatively the more advancedF-cycleexplainedat the end of Chapter4 can be used. 6.12.3 Starved Lubrication In the previouschapterson lubricationit was alwaysassumedthat pressureis generated wheneverthe gap between the two bodies narrows in the direction of motion. This supposesthat the gap between the two bodies is completelyfilled with oil, and thus that sufficientamountsof oil are presentin the inlet. Wheneverthe quantity of oil is insufficient,the generatedoil film thicknesswill be smallerthan predictedby the fully flooded theory. In order to correctlydescribethis problemthe Reynolds equation (6.1) has to be expandedto includethe parameter0 which describesthe ratio of the oil film thicknessto the gap height. As a consequence 0 is dimensionless, and 0 < 0 <_ 1. An extendedReynoldsequationwas proposedby Elrod et al. [40, 41] describingthe flow in the completeand incompletezone:

0 ph3~_~,_..~ 0 phacgp

0--~(1-~7/ ) +

~y (1-~ 0--~y) - u~

O(pOh)o:r

=0

(6.78)

Please note that this equationreducesto tile classicalequation in the completezone p ~_ 0 and 0 = 1. In the incompletezone no pressureis generated"p = 0 and 0 < 0 _~ 1. This complete equation is valid for x~ < x < Xb and -y~ < y < Ya with the boundary conditionsp(xa, y)=p(xb,y) =p(x, -Ya) =p(x, y~) = 0. These boundaryconditionshave to be extendedto include a boundaryconditionfor 0: O(xa,y) = Oo(y). When the amountof oil availablein tile inlet is reduced,the pressurebuild-upwill be delayed,and the generatedfilm thicknesswill be reduced. Chevalieret al. [31, 32] showedthat it is advantageous to expressthe solutionin terms of the amount of lubricantpresenton the surfacesin front of the contact, Hoit. Furthermore,tile starvedfilm thicknesswas expressedas a ratio, by dividingit by the fully floodedfilm thicknessHell.

Figure 6.20: Film thickness graphsusingpseudo interferometric techniques showingthe influenceof the inletoil fihnHoit,M = 100, L = 10.

222

C H A P T E R 6. E L A S T O H ~ ' D R O D Y N A M I C L U B R I C A T I O N

Figure 6.21" Film thicknessgraphs H ( X = O, }') showingthe influenceof the inlet oil film Holt, M = 100, L = 10.

Figure 6.20 shows the film thicknessevolutionas a function of the inlet parameter Hoiz. The drawn line representsthe inlet and outlet meniscus. When the amount of availablelubricantis reduced,the film thicknessdiminishes,but also the film thickness variationreduces.The film thicknessprofile becomesmore and more Hertzian. This can be observedeven better in Figure 6.21 which showsa crosssectionof the film thickness distributionfor x - 0. In order to reduce the complexity,one can study the evolutionof the central film thicknessH~. The central film thicknessdependson three parameters:M, L and Holt, the amountof oil presenton the surfaces.However,it was found that the ratio H~/Hc H is virtually independentof ~1I and L. Consequently, H~/H~ H can be expressedas a function of Hoil/Hcl! only. Figure 6.22 showsthis relationfor variousoperatingconditions,indeedshowingonly a small dependenceon M and L. The followingequationaccuratelydescribesthe film amountof oil available thicknessratio 7~ - H~/H~H , as a functionof the dimensionless on the surfacesr - Hoit/HcI I"

"R -

F

~/1 +r'Y

(6.79)

The value of 3' for EHL contactslies between2 and 5. Ill this equationand in Figure6.22, two asymptotescan be observed: r ~ 0 and r ~ cc which give 7~ ~ r and 7~ ~ 1 respectively. Of these two, the asymptotefor thin films is the most interesting,as it showsthat for very thin oil films, the EHL film thicknessis the same as the film in front of the contact. In other words: the contactbecomesvery efficient in building up an oil film when very little oil is availablein tile inlet. Thus for r ~ 0 and Hc ~ Hoil or in dimensionalterms: when ho,t/lzc// -+ 0 then hc --+ ho,l.

223

6.12. A D E 4 N C E D TOPICS 1

!/

I

r-lit,-.7

" i El I ,e~),,

-

!I

M = 10, L = I O M = 100, L - I O

0 +

M=IO00. L=tO r-i M= 10, L= 5 x M= 100, L= 5 z~ M= 10, L= 2 M= 100, L= 2 0

0.8 -

06

0.4

0.2 - i 0IIIIt1III 0 0 5 1 1.5 2

2.5 3 Hod/Hcff

3.5

4

4.5

5

Figure 6.22" Numericalcalculationsof the central fihn thicknessHc as a function of Hoit for varioussetsof operatingconditions,and a circular contact.