High order discontinuous Galerkin method for elastohydrodynamic lubrication line contact pr

HighOrderDiscontinuousGalerkinMethodfor

ElastohydrodynamicLubricationLineContactProblems

H.Lu,M.Berzins,C.E.GoodyerandP.K.Jimack∗

SchoolofComputing,UniversityofLeeds,LS29JT,U.K.

SUMMARY

InthispaperahighorderDiscontinuousGalerkinmethodisusedtosolvesteady-stateisothermallinecontactelastohydrodynamiclubricationproblems.Thismethodisfoundtobestableacrossawiderangeofloadsandisshowntopermitaccuratesolutionsusingjustasmallnumberofdegreesoffreedomprovidedsuitablegridsareused.Acomparisonismadebetweenresultsobtainedusingthisproposedmethodandthosefromaverylargefinitedifferencecalculationinordertodemonstrate

c2000JohnWiley&Sons,excellentaccuracyforatypicalhighlyloadedtestproblem.Copyright󰀁

Ltd.keywords:

elastohydrodynamiclubrication;discontinuousGalerkin;finiteelementmethod

1.INTRODUCTION

Inmostcomplexmechanicalsystemslubricantsareusedtoreducefrictionandprotectmovingpartsagainstwear.Inmanycases,thegeometryofthecontactingelementsdeterminestheshapeofthelubricantfilm.However,inthecasesofconcentratedcontactssuchasjournalbearingsandgears,thecontactingelementsdeformelasticallydefiningthefilmthickness:thisisknownaselastohydrodynamiclubrication(EHL).NumericalsolutionsofEHLareveryimportantforindustryinthedesignandevaluationofbothlubricantsandcomponents.Thesecasesrequirethesimultaneoussolutionofbothpressureandfilmthickness,alongwiththebehaviourofthelubricantundertheseextremeconditions.

ThekeyfeaturesofanEHLsolutionarethelowpressureinletregion;arapidriseinpressurethroughthecentreofthecontact,typicallyreachingthegiga-Pascalrange;acavitationboundaryintheoutflow;and,asharppressurespikepassedthecentreofthecontact,towardstheoutflow.Ithasbeenshownthatifthisspikeisnotresolvedsufficientlywellthenthecalculatedfrictioncanbeinaccurate[1].

Thefinitedifference(FD)methodisperhapsthemostwidelyusednumericalsolutionmethod.StabilitythroughwiderangesofoperatingconditionshasbeenattainedthroughuseofthedifferentFDtechniquesinthelowandhighpressureregions,asshownbyVenner[2]

2

H.LU,M.BERZINS,C.E.GOODYER,ANDP.K.JIMACK

andNurgatetal.[3].InordertoaccelerateconvergencethemultigridmethodwasfirstemployedbyLubrechtin1986[4].ForthefastcalculationoftheelasticdeformationBrandtandLubrecht[5]developedamultilevelmulti-integrationalgorithmwhichsignificantlyreducesthecomputationalcomplexityinapproximatingdeformationsateachpointinthecontact.Basedontheaboveachievements,FDcombinedwithmultigrid[2,6,7]hasbecomethemostpopularmethodforEHLproblemssinceitisbothefficientandstable.

OthertechniquesforEHLarealsoused.Theinversemethodof[8],whichismainlysuitableforhighlyloadedcases,updatesthepressureprofilefromfilmthicknessratherthantheotherwayaround.Thecoupledmethodof[9,10]calculatesthepressureandfilmthicknesssimultaneouslyinsteadofiteratingbetweenthemandcanbeappliedwitheitherFDorFE(finiteelement)discretizations.TheFEmethodhasalsobeensuccessfullyusedtosolveincompressibleEHLcontactproblemsbyanumberofauthors,e.g.[11,12,13].ForcompressibleEHLcontactproblems,upwindingisnecessaryinordertoensurestability.Forexample,in[9],oscillationsareobservedinthecentralregioninhighlyloadedcaseswhenusingatraditionalquadraticfiniteelementmethodtodiscretizetheproblem

InrecentyearstheDiscontinuousGalerkin(DG)methodhasbecomeapopularchoiceforsolvingconvection-dominatedpartialdifferentialequations[14,15,16].Sincediscontinuityisallowedoverelementinterfaces,wegettheopportunitytostablizethehigh-ordermethodbydefininganappropriatenumericalflux.DGcaneasilyhandleirregularmeshesandthedegreeoftheapproximatingpolynomialcanbeeasilychangedfromoneelementtoanother.Inthispaper,high-orderDGisusedtosolvesmoothEHLlinecontactproblems.Ournumericalexperimentsshowthatthisapproachcangiveveryaccuratesolutionswithonlyasmallnumberofunknowns:forexampleresultsobtainedusinglessthan200unknownsareshowntobecomparabletoFDresultsobtainedusingoverhalfamillionequallyspacedpoints.

2.GOVERNINGEQUATIONS

ThemathematicalmodeloflinecontactEHLproblemstypicallyconsistsofthreeequations,shownhereusingtheusualnon-dimensionalization,describedfullyin[6].TheReynoldsequationreads󰀂

dd(ρH)

−dX,PandHaretheunknownpressureandfilmthickness,ρandηarethedensity

andviscosity,andλisadimensionlessspeedparameter.(Thelubricantrheologyishighlynon-linearinpressure:inthisworkwehaveusedtheviscosity-pressurerelationshipofRoelands[17]anddensitymodelofDowsonandHiggison[18].)Theelasticityisincludedthroughthefilmthicknessequationwhichdefinesthecontactgeometryforagivenpressuresolution:

󰀄∞

X2

ln|X−X󰀈|P(X󰀈)dX󰀈.(2)H(X)=H00+

π−∞

ηλ

Theforcebalanceequationisaconservationlawfortheappliedload,givenby:

󰀄∞

π

P(X)dX−

−∞

HIGHORDERDISCONTINUOUSGALERKINMETHODFOREHLLINECONTACTPROBLEMS

3

Forphysicalreasons,allpressuresshouldbelargerthanorequaltothevapourpressureofthelubricant(zero).ThisisnotaccountedforintheReynoldsequation,hence,intheoutletregion,thecalculatedsolutionmayhavenegativepressures.ConsequentlytheReynoldsequationisonlyvalidinthepressurizedregionand,thecavitationposition,Xcav,isthereforetreatedasafreeboundary.Theboundaryconditionsare:

dP

P(Xin)=0,P(Xcav)=0and

(7)(v|∂Ωe∩Γef+v|∂Ωf∩Γef).

2

Followingtheapproachof[14]and[15],adiscreteformoftheReynoldsequationbecomes:

a(P,v)=l(P,v),

wherea(P,v)=

󰀃󰀁󰀄

󰀂

󰀄

(8)

󰀁v·󰀣󰀁Pdx+

Ωe

([P]󰀉(󰀣󰀁v)·n󰀊−[v]󰀉(󰀣󰀁P)·n󰀊)ds

Γint

Ωe∈Ph

+

andl(P,v)=

󰀃󰀁󰀄

󰀂

󰀄

󰀄

(P(󰀣󰀁v)·n−v(󰀣󰀁P)·n)ds,(9)

ΓD

(󰀁v·β)ρHdx−

Ωe

v(ρ(P−)H)(β·ne)ds

∂Ωe\\Γ−

Ωe∈Ph

−

󰀄

vρ(g)H(β·n)ds+

Γ−

󰀄

(󰀣󰀁v)·ngds.(10)

ΓD

c2000JohnWiley&Sons,Ltd.Copyright󰀂

Preparedusingcnmauth.cls

Commun.Numer.Meth.Engng2000;00:1–6

4

H.LU,M.BERZINS,C.E.GOODYER,ANDP.K.JIMACK

Intheaboveequations,

P−=limP(x−δβ),forx∈Γint

δ→0

(11)

andnistheouterunitnormalandn=neonΓint.

Ineachelement,Pisexpressedinthefollowingform:

Pe=

ep+1󰀃

e

ueiNi,

(12)

i=1

wherepeistheorderoftheapproximatingpolynomial,ueiaretheunknowncoefficientsand

e

Niarethelocalbasisfunctions.Thediscontinuousbasisfunctionsusedherearedescribedin[19].Notethatin(10)upwindingiseasilyimplementedbypickingtheleftvalueρ(P−)whencalculatingthenumericalflux(seethethirdtermin(10)).3.2.TheFilmThicknessEquation

Foragivenpressuredistributionthefilmthicknessmaybecalculatedasfollows:

󰀄Xcav

󰀂󰀂󰀂X2

H(X)=H00+ln|X−X|P(X)dX

πXin

N󰀄󰀃󰀂󰀂󰀂X2

=H00+ln|X−X|P(X)dX

πe=1Ωe

=H00+=H00+=H00+

X

2

(13a)(13b)(13c)(13d)(13e)

πX2

πX

2

N󰀄󰀃e=1

ln|X−X|

Ωe

󰀂

+1󰀄Np󰀃󰀃

e

e

p+1󰀃

e

ueiNi(X)dX

󰀂󰀂

i=1

󰀂

e=1i=1

e

Ωe

ln|X−X|Nie(X)dXuei

󰀂󰀂

π󰀄

+1Np󰀃󰀃e=1i=1

e

Ki(X)uei,

e

wheretheKi(X)aredefinedby:

e

Ki(X)

=

Ωe

ln|X−X|Nie(X)dX.

󰀂󰀂󰀂

(14)

e

HereKi(X)iscalculatedusingnumericalintegration:forX∈esingularquadratureis

󰀂󰀂

employedsinceln|X−X|hasaweaksingularityatX=X;elsewhereGaussianquadratureissatisfactory.

3.3.TheForceBalanceEquation

Theforcebalanceequationisdiscretizedaccordingto:

N󰀄󰀃e=1

e

p+1󰀃

Ωei=1

e

ueiNi(X)dX−

π

HIGHORDERDISCONTINUOUSGALERKINMETHODFOREHLLINECONTACTPROBLEMS

5

4.SOLUTIONPROCEDURE

Thissectiondescribesasimplesolutionprocedurethatmaybeusedonagivenmeshandwithagivenchoiceofthepolynomialdegreeoneachelement.

TheunderlyingrelaxationschemeassumesthatifXcavisgiventhenthediscreteformof(8),alongwiththefirsttwoconditionsin(4),maybewritteninthegeneralform:

L(U)=A(U)U−b(U)=0,

(16)

1NN

whereU=(u11,...,up1+1;...;u1,...,upN+1)aretheunknownpressurecoefficients.NotethatbothA(U)andb(U)dependonU.Thepressureisrelaxedaccordingto:

󰀁∂L

U←U+

∂U

isapproximatedby:

(18)

∂L

∂U

whichisafullmatrixsinceH(X)hasaglobalpressureintegral.

Thecavitationposition,Xcav,mustbechosensoastoimposethethirdconditionin(4).ThisisachievedbymovingtheentiregridaccordingtothecurrentvalueofdP

dP

dX>0thegridismovedleftandif

dX

(Xcav),

(19)

whereδisasmallpositiveconstant.

Theaboveingredientsareusedwithinthefollowingoverallsolutionprocedure.(i)(ii)(iii)(iv)(v)

StartwithaninitialchoiceofXcavandaninitialpressuredistribution.RelaxH00usingtheerrorinequation(15),asexplainedin[6]forexample.ComputethethicknessprofileH(X)from(13e).

Updatethepressuredistributionbyrepeating(18)untilconvergence.UpdatethecavitationpositionaccordingtothevalueofdP

dX(Xcav)=0isobtained.

5.NUMERICALRESULTS

InthissectionwedemonstratethepotentialoftheDGmethodbycomparingsolutionsforahighloadtestproblemagainstthoseobtainedusingastandardmulti-level,multi-integrationFDalgorithm.Itisshownin[1]thatinordertofullyresolvethepressurespikeuptohalfamillionFDgrid-pointsmayberequired.HerewecompareourDGresultsagainstincreasingresolutionsofFDgrids.FortheDGsolution16elementsareused(notofequalsize)andthepolynomialdegreeiseither10(inthepressurespikeregion)or8(elsewhere).

Figure1showsthepressureprofilecomputedforatypicalhighlyloadedcase.Theentirecontactisshownintheleftgraphwhilstadetailedviewofthepositionofthepressurespike

c2000JohnWiley&Sons,Ltd.Copyright󰀂

Preparedusingcnmauth.cls

Commun.Numer.Meth.Engng2000;00:1–6

6

1 0.9 0.8Non-dimensional pressure, P 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0-3

-2.5

H.LU,M.BERZINS,C.E.GOODYER,ANDP.K.JIMACK

FD 524289DG ???dof

1 0.9 0.8Non-dimensional pressure, P 0.7 0.6 0.5 0.4 0.3 0.2

-2-1.5-1-0.5 0 0.5Non-dimensional distance through contact, X

1

1.5

0.1 0.9

FD 524289FD 262145FD 131073FD 65537FD 16385FD 8193FD 4097FD 2049FD 1025FD 513FD 257DG ???dof

0.905 0.91 0.915Non-dimensional distance through contact, X

0.92

Figure1.PressuredistributionsobtainedusingDGandFDmethodsacrosstheentirecontact,left,

andaroundthepressurespike,right

Method

409781931638565537131073262145524289

DGFE

PeakPressure

0.90690.90840.90920.90660.90970.90970.9097

0.9169

FreeBoundaryPosition

1.0707

TableI.ComparsionofPressurePeakPositionandFreeBoundaryValues

isshownontheright.Thekeyfeaturesofinterestarethepeakvalueofpressure,itspositionandthepointatwhichthefreeboundaryoccurs.ThesevaluesareshowninTable1.ItisclearfromtheseresultsthattheDGsolutionisclosesttothehalfmillionpointFDresults.

Notethatsolutionefficiencyisnotconsideredinthiscommunication.TheFDcalculationswouldobviouslybenefitfromtheuseoflocallyrefinedmeshes,suchasusedin[1,20],whilsttheDGmethodcouldcertainlybeimplementedwithamoreefficientsolutionalgorithm(baseduponp-multigrid,[21]forexample).WhathasbeendemonstratedhoweveristhattheDGapproachcandeliververyhighaccuracyusingonlyasmallnumberofdegreesoffreedom.

6.DISCUSSION

InthispaperahighorderDGalgorithmisusedtosolvetheEHLlinecontactproblemforthefirsttimeand,unlikefortraditional(continuous)FEmethods,stabilityiseasilyobtainedthroughtheuseofasimpleupwindingstep.Ournumericalresultsshowthat,whensuitablegridsareused,highlyaccuratesolutionsmaybeobtainedwithveryfewdegreesoffreedom,thusillustratingtheexcitingpotentialofthisapproachforthesolutionofEHLproblems.Otherresultsofsimilarqualityhavebeenobtainedunderdifferentloadingconditionsandusingdifferentsolutionalgorithms(e.g.thepenaltymethoddescribedin[22]).

c2000JohnWiley&Sons,Ltd.Copyright󰀂

Preparedusingcnmauth.cls

Commun.Numer.Meth.Engng2000;00:1–6

HIGHORDERDISCONTINUOUSGALERKINMETHODFOREHLLINECONTACTPROBLEMS

7

Fortheresultsobtainedinthispaperaverycrudeadaptivestrategywasusedtoselectthemeshandthelocalpolynomialdegree.Currentlyweareinvestigatingmoresophisticatedh-p-refinementstrategiesforthismethod.Wearealsointheprocessofdevelopingamultigridversionofthesolutionalgorithmbaseduponthep-multigridmethod,describedin[21]forexample,andwillalsoconsider2-dpointcontactcases.

REFERENCES

1.GoodyerCE,FairlieR,HartDE,BerzinsM,andScalesLE.CalculationofFrictioninSteady-StateandTransientEHLSimulationsInA.Lubrechtetal.,editor,TransientProcessesinTribology:Proceedingsofthe30thLeeds-LyonSymposiumonTribology.Elsevier,2004.

2.VennerCH.MultilevelSolutionoftheEHLLineandPointContactProblems.PhDthesis,UniversityofTwente,Netherland,1991.

3.NurgatE,BerzinsMandScalesLESolvingEHLproblemsusingiterative,multigridandhomotopymethods.J.Tribol-TASME,121:28–34,1999.

4.LubrechtAA,TenNapelWE,BosmaR.Multigrid,analternativemethodforcalculatingfilmthicknessandpressureprofileinelastohydrodynamicallylubricatedlinecontacts.J.Tribol-TASME,108:28–34,1986

5.BrandtA,LubrechtAA.Multilevelmatrixmultiplicationandfastsolutionofintegralequations.J.Comput.Phys.,90:348–370,1990

6.VennerCH,LubrechtAA.MultilevelMethodsinLubrication.Elsevier,2000.

7.JingWang,ShiyueQu,PeiranYang.Simplifiedmultigridtechniqueforthenumericalsolutiontothesteady-stateandtransientEHLlinecontactsandthearbitraryentrainmentEHLpointcontacts.TribologyInternational,34:191–202,2001.

8.KwehCC,EvansHP,SnidleRW.Elastohydrodynamclubricationofheavilyloadedcircularcontacts.P.I.Mech.Eng.C-JMEC,203:133–148,1989.

9.HughesTG,ElcoateCD,EvansHP.Anovelmethodforintegratingfirst-andsecond-orderdifferentialequationsinelastohydrodynamiclubricationforthesolutionofsmoothisothermal,linecontactproblems.Int.J.Numer.Meth.Eng.,44:1099–1113,1999.

10.ElcoateCD,EvansHP,HughesTG.Onthecouplingoftheelastohydrodynamicproblem.P.I.Mech.

Eng.C-JMEC,212:307–318,1998.

11.TaylorC,O’CallaghanJF.Anumericalsolutionoftheelastohydrodynamiclubricationproblemusing

finiteelements.J.Mech.Eng.Science.14:229–237,1972.

12.RohdeSM,OhKP.Aunifiedtreatmentofthickandthinfileelastohydrodynamicproblemsbyusing

higherorderelementmethods.Proc.R.Soc.Lond.A.343:315–331,1975.

13.WuSR,OdenJT.Anoteonapplicationofadaptivefiniteelementstoelastohydrodynamiclubrication

problems.Commun.AplliedNumer.Meth.3:485–494,1987.

14.OdenJT,BabuskaIandBaumannCE.Adiscontinuoushpfiniteelementmethodfordiffusionproblems.

J.Comput.Phys.146:491–519,1998.

15.BaumannCE,OdenJT.Adiscontinuoushpfiniteelementmethodfordiffusionproblems.Comput.

Methods.Appl.Mech.Engrg.175:311–341,1998.

16.CockburnB,KarniadakisGE,ShuCW.DiscontinuousGalerkinMethods:Theory,Computationand

Applications.Springer,1999.

17.RoelandsCJA,VlugterJC,WatermannHI.Theviscositytemperaturepressurerelationshipoflubricating

oilsanditscorrelationwithchemicalconstitution.ASMEJ.BasicEngng.85:601,1963.

18.DowsonD,HigginsonGR.Elasto-hydrodynamicLubrication,TheFundamentalsofRollerandGear

Lubrication.PermagonPress,Oxford,GreatBritain,1966.

19.SzaboB,BabuskaI.IntroductiontoFiniteElementAnalysis.JohnWileyandSons,Inc,1989.

20.GoodyerCE,FairlieR,BerzinsM,andScalesLE.Adaptivemeshmethodsforelastohydrodynamic

lubrication.InECCOMASCFD2001:ComputationalFluidDynamicsConferenceProceedings.InstituteofMathematicsanditsApplications.ISBN0-905-091-12-4,2001.

21.FidkowskiKJ,DarmofalDL.Developmentofahigher-ordersolverforaerodynamicapplications.42nd

AIAAAerospaceSciencesMeetingandExhibit,AIAAPaper2004-0436,2004.

22.WuSR.ApenaltyformulationandnumericalapproximationoftheReynolds-Hertzproblemof

elastohydrodynamiclubrication.Int.J.Engng.Science.24:1001–1013,1986.

c2000JohnWiley&Sons,Ltd.Copyright󰀂

Preparedusingcnmauth.cls

Commun.Numer.Meth.Engng2000;00:1–6