Design procedures based on numerical methods for hydrodynamic lubrication

237

PaperVI II(iv)

Design procedures basedon numerical methodsfor hydrodynamic lubrication J. 0. Medwell

T h i s p a p e r p r e s e n t s a r e v i e w of numerical schemes for predicting the performance of cylindrical bore journal bearings which have been developed under the initial guidance of F. T. Barwell at Swansea. The schemes are based on the finite element and finite difference methods which are used t o solve the lubricant momentum and energy equations. T h e effects of recirculatory flow, bush conduction and deformation, shaft expansion and misalignment have been considered.

INTRODUCTION B C

P

=

-

C

DP

P r R Re T U

U V W X,Y,Z

8

V'

al

B

e fl Y

Y

X

Y Y

z

8 l . U P Q,

@ 0

si

bush thickness radical clearance specific heat journal diameter Young's Modulus local film thickness thermal conductivity bearing length mesh coordinate point local pressure radial coordinate journal radius WRC film Reynolds number [= 7 1 temper at u r e v e l o c i t y c o m p o n e n t in x direct ion velocity vector ui + U] + wk v e l o c i t y c o m p o n e n t in y d i rec t i on v e l o c i t y c o m p o n e n t in z direct ion Cartesian coordinates see Figure l(a) or circumferential coordinate the Laplacian operator see Figure l(a) see Figure l(a) eccentricity ratio molecular viscosity kinematic viscosity turbulent viscosity turbulent viscosity Poisson's ratio parametric coordinate directions lubricant density see Figure l(a) disslpatlon function rotational speed

The theory of operation of hydrodynamic j o u r n a l b e a r i n g s w a s we1 1 established by Reynolds over a century ago. H o w e v e r , a s p o i n t e d o u t b y B a r w e l l et a 1 [ I ] the application of this theory to design w a s not simple because I t is not possible to o b t a i n a general analytical solution to the Reynolds equation i t i s d i f f i c u l t to provide the correct value of viscosity for use in the e q u a t i o n because of its temperature dependence. Strictly, s i n c e large differences in temperature can arise the consequent variation of viscosity can only be determined by i n c l u d i n g the energy equation, c o n t a i n i n g t h e a s s o c i at ed dissipative terms, in the design procedure. This usually involves u s i n g i t e r a t i v e t e c h n i q u e s in order that the Reynolds and the e n e r g y e q u a t i o n s can be solved simultaneously.

t h e a v a i l a b i l i t y of modern day numerical and computational methods, this does not p o s e a n y particular problem. The numerical techniques which have been used most o f t e n a r e f i n i t e d i f f e r e n c e methods a n d , latterly, in the last decade or so, finite element methods. The former has been used in the generation of design charts and also forms t h e b a s i s of s o f t w a r e s e r v i c e s such as provided by ESDU. However, although these may be u s e d w i t h confidence in normal journal bearing d e s i g n c a l c u l a t i o n s t h e y c a n b e suspect w h e n o f f design conditions are met. Typical of these are: recirculatory flows that can occur in the highly eccentric lobes

of profile bore bearings, shaft misalignment and t h e r m o - e l a s t i c distortion in bearing bushes and shafts. T h i s p a p e r will d'iscuss h o w t h e s e conditions may be incorporated into standard numerical design procedures that have been developed, initiated by the late Professor F T Barwell. at Swansea. In this respect most of the b o u n d a r y conditions imposed in the solution procedures relate t o t h e e x t e n s i v e e x p e r i m e n t a l p r o g r a m m e s that h a v e been carried out in Swansea. These have b e e n c o n c e r n e d w i t h a s s e s s i n g t h e p e r f o r m a n c e of s t a n d a r d cylindrical bearings which w e r e fed by two diametrically opposed axial grooves at 90' to the load line as shown in Figures ](a) and I(b). Wide ranges of speeds (both laminar and turbulent lubricant films) clearance ratios and load capacities were investigated (see Ref [2]). However, in the interests of brevity, detailed boundary conditions have been omitted and attention focused on the governing equations and solution strategies. Design Procedures The equations of motion and energy of an incompressible lubricant flow are based on the conservation l a w s of m a s s , m o m e n t u m and e n e r g y . Generally, the three dimensional momentum equations may be written as

[31 and continuity 0.u

= 0

[41

Similarly, the energy equation may be written as pCp W.T

=

kV2T +

U@

In t h e above equations the lubricant film has been assumed to be laminar in motion and therefore the transport quantities (u and k) and the viscous dissipation term @ reflect t h e s t r e a m l i n e nature of the flow. I t is possible with extensive modification to cast the equations in a form suitable for turbulent motion. This can be accomplished by replacing the molecular properties w i t h corresponding t u r b u l e n t v a l u e s , and t h e n i n t r o d u c i n g a p p r o p r i a t e models to describe them. The additional e q u a t i o n s r e q u i r e d g e n e r a l l y i n c r e a s e s the difficulties in obtaining a solution to such problems. However, as mentioned earlier, in many practical applications the above equations can be s i m p l i f i e d c o n s i d e r a b l y by making appropriate assumptions. Typical of these are the neglect of cross stream velocities, all

fluid advection terms, all conduction effects in the lubricant film. Thus the problem is reduced to solving the Reynolds equation and an energy equation which simply balances the generated heat with that convected through the f i l m . A d d i t i o n a l l y , t u r b u l e n c e c a n be accounted for by r e p l a c i n g t h e m o l e c u l a r v i s c o s i t y by e x p r e s s i o n s that m o d i f y i t depending on the position considered in the film [3]. The above e q u a t i o n s m u s t be solved s u b j e c t t o t h e v a r i o u s i m p o s e d boundary conditions which will depend on the type of equations to be solved. For example, in their i n v e s t i g a t i o n of t h e r o l e o f i n e r t i a in lubricant films, Medwell et a1 [4] used the finite element method w h i c h , because of the elliptic nature of the equations of motion, required t h e s p e c i f i c a t i o n s i x b o u n d a r y values. These were either prescribed or cast in the form of gradients a n d , in some cases, updated during the solution procedure since t h e y c o u l d not b e s t a t e d a p r i o r i . A s mentioned earlier, the solution procedure is necessarily iterative in form because of the dependence of the lubricant film profile on the bearing attitude a n g l e and t h e r e f o r e r e q u i r e s c o n t i n u o u s u p d a t i n g . In t h e s e investigations the three dimensional films profile w a s generated conveniently by scaling a unit cube (suitably discretized) to produce the a p p r o p r i a t e film geometry through the relationships

w h e r e the overbar denotes a fixed coordinate. T h i s technique w a s first proposed by Gethin [ 2 ] and has been s u b s e q u e n t l y u s e d successfully in analysing other performance aspects of journal bearing operation [5], [6] and [7]. In all analyses i t is c o n v e n i e n t t o unwind the loaded part of the lubricant film because the radius of curvature of the film is v e r y l a r g e in c o m p a r i s o n w i t h t h e f i l m thickness Figures l(a) and l(b) s h o w s t h e journal bearing geometry and the unwrapped film discretized in a form suitable for the f i n i t e e l e m e n t m e t h o d of solution, using eight-node isoparametric elements. If heat conduction in the journal bush is important it is necessary t o include the appropriate f o r m of the energy equation and boundary conditions into the solution viz.

V'T

=

0

[81

wlth the inclusion of conduction terms in the fluid film also. A suitable discretization of the bush is also depicted on Figure l(c).

239 So far, the equations that have b e e n d i s c u s s e d a r e f u l l y three dimensional in nature and reference has only been made to the finite element method of solution. This i s because finite difference methods [ 8 ] and [ 9 ] solve parabolized form of the equations of motion w i t h solution stability being maintained using upwind techniques. Where recirculatory flows occur the v a l i d i t y o f these procedures becomes questionable. The importance of recirculatory flows has been demonstrated for isothermal lubricant films [ 4 ] by solving the equations of motion in which only fluid advertion across the film w a s ne.glected. A l t h o u g h t h e a n a l y s i s c o n f i r m e d that lubricant inertia does not affect bearing performance significantly i t does indicate that large areas of recirculation occur at h i g h e c c e n t r i c i t y ratios as shown on Figure 2. As pointed out earlier, this phenomena can be expected to appear in the lobes of profile bore bearings which can be operating in a highly eccentric c o n d i t i o n even though the overall bearing eccentricity is modest. T h e p r o b l e m in the type of rigorous analysis mentioned above is that i t can be e x p e n s i v e to compute and which can become prohibitive if thermal and elastic behaviour are included. Obviously, it i s important to the designer to use methods which will compare favourably with both the rigorous analyses and experimentally determined performances and be e c o n o m i c c o m p u t a t i o n a l l y . A n excellent exposition of this has been given by Gethin [5] w h o presented two models. One was based o n the s o l u t i o n of t h e f u l l e q u a t i o n s including bush conduction (i.e. equations 1 4 and 8) and in the other the equations of motion w e r e reduced to a generalized Reynolds equation, viz.

a

[Gg] +a z[ G g ]

= O R zdh iji-

dF dx

r91

where G

=

1:

(y-F)dy; F

=

FI FO

where

A consequence of embodying the laws of c o n s e r v a t i o n of m a s s and momentum into a Reynolds equation is the simplification of the boundary conditions which were now in the form of specified pressrues or gradients o f pressure. The imposition of several different thermal boundary conditions w a s assessed and discussed fully while the viscosity dependence on temperature was 'inbuilt' into the solution procedure in the form of the Walther equation. The strategy for the solution of both s e t s o f e q u a t i o n s w a s b r o a d l y the s a m e consisting of the follawing basic steps: (i)

the bearing geometry and boundary conditions for the hydrodynamic e q u a t i o n s w e r e s e t up for a n assumed value of attitude angle

which w a s used to generate a film thickness profile for a specified eccentricity ratio. The film thickness equation used was

h = C (1 +

E COS

(0 +

01

- fb)}

1101

(11) the load w a s calculated and the n e w attitude angle w a s evaluated which was then used to update the film geometry. This procedure was repeated unt i 1 convergence w a s achieved and any parameters such as velocity g r a d i e n t s and/or pressure gradients that may be required for s u b s i t u t i o n into the energy equation were calculated. (iii) the boundary conditions for the energy equation w e r e set u p and this w a s solved to give temperature distributions. (iv) the nodal values of viscosity were updated throughout the film and t h e s t e p s ( i ) , (ii), ( i i i ) and ( i v ) repeated until e i t h e r t h e velocity o r pressure fields converge to a p r e s c r i b e d tolerance. Some of the results from the two models are displayed on Figures 3 and 4 together with some experimental results. I t can be seen that generally the more sophisticated model does not p r e d i c t r e s u l t s w h i c h a r e s i g n i f i c a n t l y more accurate than those predicted by the method based o n R e y n o l d s e q u a t i o n . However, in terms of computing requirements the latter model used just 10% of the CPU time of the more rigorous analysis - a very significant saving. T h e s i m u l t a n e o u s solution of the Reynolds and energy equations using the above s t r a t e g y f o r m s t h e b a s i s f o r m a n y other journal bearing performance predictions using b o t h f i n i t e element and finite difference techniques. It w a s first used [ l o ] to assess the influence of heat dissipation in high speed b e a r i n g s w h e r e t h e lubricant f i l m w a s in turbulent motion. Thus was done by using very crude turbulence models [ 3 ] , viz. uX = and

Y

(1 +

0.000175 Re

1.06

uZ = u ( 1 + 0.000375 Re

)

1.08

)

to replace molecular viscosity terms in the R e y n o l d s e q u a t i o n s . T h e energy equation employed neglected crossfilm and streamwise cnduction effects (the adiabatic condition) which r e s u l t s in t h e t e m p e r a t u r e o f t h e lubricant film at any point being represented by a mean value. Both equations were cast in conventional finite d i f f e r e n c e f o r m b a s e d o n c e n t r a l differences. The Reynolds equation was solved

240

\A

90

L

E l

la)

Solution of full equation Solution based on Reynolds equation Experiment

(b)

--0

Journal bearing geometry Shaft surface

Unwrapped lubricant film and finite element mesh f o r fluid and solid FIGURE 1 12 l6

1I 0

/

Isothermal

02

01

04

03

05

07

06

08

E FIGURE 2

RECIRCULATING FLOW AT BEARING INLET

Solution of f u l l equation Solution based on Reynolds equation

FIGURE 5

THEORETICAL VARIATION OF LOAD WITH ECCENTRICITY RATIO FOR ISOTHERMAL AND ADIABATIC FLOW SPEED = 30,000 rev/min

--_

10

-

5

Isothermal

-

0

03

05

07

09

Eccentricity r a t i o FIGURE 3

VARIATION OF LOAD WITH ECCENTRICITY RATIO Bearing speed = 10,000 rev/min.

$

0.002 ,

‘/R = 1.0

0 1 0

I 10

I

Speed FIGURE 6

I

30

20 (x

I 40

I 50

I

60

lo3 rev/min.)

THEORETICAL VARIATION OF LOAD WITH SPEED FOR ISOTHERMAL AND ADIABATIC FLOW. E = 0.7

24 1

using a Gauss Seidel iterative technique while the simpler form of the energy equation allowed a straight through march to be used in obtaining a solution.

In polar coordinates the radical and circumferential components of strain w e r e written in terms of a displacement matrix as

Some of the findings of the investigation are displayed on Figures 5 and 6 .

1

The above method was extended by Medwell and Gethin [ll] to account for off-design bearing operation such as misalignment. Identical Reynolds and energy equations were solved using the same models of turbulent viscosity and employing a modified film thickness equation which included an extra term to account for the misalignment. This i s measured as the inclination of the shaft axis to that of the bush and is measured in the plane in which the shaft i s loaded. The main consequence of misalignment i s shown in Figure 7 where maximum and bulk temperatures of the lubricant film are plotted as a function speed for an eccentricity of 0 . 6 5 and a misalignment angle of only 0.12'. I t can be seen that although the bulk temperatures exhibit modest rises over the inlet lubricant temperature the maximum temperatures generated are dangerously high and would certaily be a limiting parameter in such a bearing operation and which confirms the importance of misalignment effects in journal bearing applications.

Generally, when both finite element and finite difference methods are applied to similar problems as above, it is found that the difference methods are less demanding on computational requirements. However, as pointed out earlier, when recirculatory flow occur,s the governing equations lose their parabolic form and the vality of upwind techniques introduced to achieve solution stability in the difference methods becomes questionable. A further disadvantage of the finite difference method is its inability to cope efficiently with coupled problems such as those where bush conduction and thermo-elastic deformations are important. I t is in these types of situations that the finite element method is particularly effective. This derives from the capacity of the method to cope which property changes in adjacent elements either side of the fluidsolid interface that makes i t extremely attractive to designers. An excellent example of this has been given by Gethin [12] in his theoretical investigation into the effect of thermoplastic bush deformation and shaft thermal e x p a n s i o n o n j o u r n a l bearing performance. He accounted for the distortion of the bush due to the temperature and pressure fields by introducing the relevant equations of solid mechanics. However, based on experimental evidence, computational economy w a s achieved by assuming that the distortion occurring on the bearing centre line only could be applied over the entire width of the bearing.

I

I

1

II

l a

The stresses were related to strain by the constitutive equations

or {o} = [D] {el}

[I21

where e l ,as and 7 r 6 arise from the applied pressure and thermal loading. Applying an energy minimization principle and equations ( 1 1 ) and (12) an expression of the following form was obtained In[BIT[D]

[B] (6)dQ

=

I I

pdQ - qdr

[I31

or [Kl (6)= {F} which can be solved for {a}, the displacement of the domain. In equation ( 1 3 ) the integral over the solid domain includes the thermal loading while the surface integral over r contains contributions arising from the imposed pressure that is at an elemental level.

where f(T) i s the temperature field over the element edge. Similarly

where f(p) is the pressure field over the element edge. The applied pressure obtained from solution of the equations while the thermal determined from solution of equation in the bush.

loading w a s hydrodynamic loading was the e n e r g y

A further a s s u m p t i o n m a d e in the analysis was that the bush was constrained in a rigid housing s o that the computational domain was confined to the load carrying zone of the bearing. The solution procedure follows that outlined earlier with the additional steps requried to calculate the distortion of the bush and shaft expansion which were then included in the calculation of the new film thickness. Gethin presented a large amount of data to illustrate the e f f e c t o f i n c l u d i n g distortion for a wide range of speeds (he assumed that the film w a s turbulent) and

242

which highlight the success of his method is the maximum film temperature variation with shaft speed as shown on Figure 8.

CONCLUSIQNS This paper has been concerned with the u s e of finite element and finite difference methods in the prediction of journal bearing performances. Several aspects of journal bearing operation have been considered and i t is concluded that the finite element method is a very powerful design tool which can cope w i t h many a.spects of bearing behaviour hitherto neglected.

REFERENCES

10. Bowen, E.R. a n d M e d w e l l J . O . ; "A therrnohydrodynamic analysis of journal bearings operating under turbulent conditions." Wear, Vo1.51, 345, 1977. 11. Medwell, J . O . a n d G e t h i n , D . T . ; "Synthesis of thermal effects in misaligned hydrodynamic j o u r n a l bearings". Int. J. Num. Methods in Fluids, Vo1.6, 445, 1986. 12. Gethin, D.T.; "An investigation into plain journal bearing b e h a v i o u r inlcuding thermo-elastic deformation of the bush." Proc. Instrn. Mech. Engrs. Vo1.199, 215, 1985.

1. Barwell, F.T. and Medwell, J.O.; "Journal B e a r i n g C a 1 cul a t ions" , Euromech

Maximum temperature

Colloquium No. 1 2 4 , T u r i n , Italy, 1979.

150

Gethin, D.T.; "Superlaminar f l o w i n j o u r n a l bearings" PhD Thesis, University of Wales, 1983.

U

3. Ng, C.N. and Pan, C.H.T.; "A linearised turbulent lubrication theory." Trans. A S M E , J . Basic Engineering. Vol. 87.675.1965.

n

2.

4. Medwell, J.O., Gethin D.T. and Taylor, C.; "A finite element analysis of the Navier-Stokes equations applied to high speed thin film lubrication." Trans. ASME J. Tribology Vol.109, 71, 1987.

e E!

+

2

100

Q

20

30

FIGURE 7

50

40

Shaft speed rev/min.

x

lom4

VARIATION OF MAXIMUM AND BULK TEMPERATURES WITH BEARING SPEED Inlet temperature = 40'C Shaft misalignment = 0.12" Eccentricity ratio = 0.65

5. Gethin, D.T.; "A finite element approach to ana 1 ys i ng thermohyd r o d y h a m i c lubrication in journal bearings. Tribology International, Vo1.21, 67, 1988 6. . Gethin, D.T. and Medwell, J.O.; "Analysis of high speed bearings operating with incomplete films." T r i b o l o g y International Vo1.18, 340, 1985. 7. Medwell, J.O. and Gethin D.T.; "A finite element analysis of journal bearing lubrication." Proc.lOth Leeds-Lyon Symposium on Tribology (Eds. D Dowson et al), Butterworth Scientific, UK, 1983.

90

-e

-

80 -

U

$

70

No deformation 0 Bush deformation 0 Shaft bush A

4-

8. Launder, B.E. and Leschziner, M.A.; "An efficient numerical scheme for the prediction of turbulent flow in thrust b e a r i n g s " Proc. 2nd Leeds-Lyon Symposium on Tribology (Eds. D Dowson et al) Mech Eng Publications, UK, 1975.

m

aI L

E60

-

e

50

-

0

9. King, K . F . a n d T a y l o r , C . M . : " A n estimation of the effect of f l u i d inertia on the performance of the plane inclined slider thrust bearing with particular regard to turbulent lubrication" Trans. ASME, J . Lub. Technology, Vol.99, 129, 1979.

10

I

I

20

30

Shaft speed rev/min x103 FIGURE 8

VARIATION OF MAXIMUM FILM TEMPERATURE WITH ROTATIONAL SPEED .

'/o

'/o

0.002, = 0.5, B = lOmm Eccentricity ratio = 0.5