Recent developments in the theory of squeeze film bearings

Recent developments in the theory of squeeze film bearings F.P. da Silva and E.M. Almas*

Squeeze film air bearings have not received much attention for practical applications due mainly to load capacity limitations. However, it is possible that a better knowledge of several critical parameters might improve the load capacity. In this work the aim is to define some of these critical parameters and discuss the applicability of some methods-of analysis. Keywords: squeeze film bearings, load capacity, finite element method The first report on gas squeeze film theory is due to Tipei ~ . Since then, several contributions referring to theoretical and experimental work have been produced by Langlois 2 , Michael 3 , Salbu4 and Pan s'6 . In the 1960s squeeze film bearings had their main applications in guidance gyroscopes for ballistic missiles, and IBM research teams also tried to develop applications for computer systems. As a result of difficulties in obtaining significant values of load capacity the interest for this type of bearing has decreased. Improvements in prediction might, however, be obtained using techniques that have progressed significantly in the last decade and become common place in dynamic analysis, eg the finite element method (FEM).

There are actually two, or possibly three, independent components in these bearings - the shell containing the fluid (usually air) vibrating at high frequency, the source of this vibration, for instance a piezoelectric crystal, and the fluid itself. The behaviour of the fluid is well defined by the Reynolds equation including terms for the time-dependent squeeze film, and its solution does not present special difficulties. Michael 3 discussed extensively the possible ways of obtaining an approximate solution of the Reynolds differential equation by its substitution for a system of finite difference equations solved by numerical methods. One of the methods, by Crank and Nicolson 7 is presented in this paper. The work of Huxley a and, more recently, Pina da Silva9 examines the influence of the supporting surface

*CEMUL, Centro de Mecanica e Materials da Universidade Tdcnica de Lisboa, Av. Rovisco Pais, 1096 Lisboa, Codex, Portugal

List of symbols

Inl IKI

(P) (p) fnEx

fnMC fnVE 17

?, eo t gzz

ri/

re/ ri l

R/ At T

Global mass matrix Global stiffness matrix Global force vector Displacement vector Acceleration vector Experimental value for first natural frequency First natural frequency from Marcus-Goldberg method First natural frequency from finite element method Structural damping loss factor Equivalent viscous damping coefficient Voltage applied to transducer electrodes Transducer thickness Piezoelectric constant Internal radius of circle f External radius of circle / Force in circle/" of the interface transducer/ supporting member Position coordinate of F/ Resulting force in the interior node of circle j Resulting force in the exterior node of circle j Time increment First natural period of vibration Excursion ratio

TRIBOLOGY international

A ho cr /2 w B Pa x h p t U X T P H qz u;

Amplitude of oscillation Mean film thickness Squeeze number Fluid dynamic viscosity Frequency of oscillation Film length in the cone generator direction Ambient pressure Space coordinate Film thickness Film pressure Time coordinate Bearing surfaces relative tangential velocity Dimensionless space coordinate Dimensionless time coordinate Dimensionless film pressure Dimensionless film thickness Product PH Function representing product PH at node

(JX, nT) 0, 0'

Nt Winst W a

Finite difference formulation parameters Time divisions for one period of oscillation Instantaneous load capacity Bearing load capacity Cone half apex angle

0301-679X/86/010079-08 $03.00 © 1986 Butterworth & Co (Publishers) Ltd

79

da Silva and Almas - theory o f squeeze film bearings

dynamic behaviour on the performance of this type of bearing. This last reference considers discs and conical shapes. For the conical bearing, the theoretical prediction of the natural frequencies and the modes of vibration uses the Marcus-Goldberg 1° method and for the same experimental amplitudes of vibration, the modal shapes are plotted along the radial axis. Because the vibrational behaviour in such an important factor in the bearing output the load capacity - an investigation has been carried out for a conical squeeze air film bearing with geometrical and material parameters similar to those already used in previous work 9 . For this particular bearing (Fig 1), a complete theoretical study of its load capacity, including the evaluation of the amplitudes of vibration of the supporting member has been performed and is now presented.

L

R=25.5mm

I-

L ~ =

A discretization of the support member into finite elements enabled the continuum to be approximated by a dynamic system with a finite number of degrees of freedom, the solution then being obtained by a numerical method. The proposed theoretical model includes the forces generated by the piezoelectric transducer. The amplitudes of vibration of the finite element nodes at the resonant frequency are then obtained. Actually only the values corresponding to the interior surface of the support member contribute to the determination of the pressure distribution in the air film and the bearing load capacity. Assuming that the supported member of the bearing is stationary, the air film thickness between the bearing surfaces will be a function of the mean film thickness and the amplitude of vibration of the support member. This member is flexible and so has a non-uniform amplitude of vibration. The amplitude distribution is given by a set of values as described before and the film thickness will then be a time-dependent variable, non-uniform along the cone generator and described by a set of discrete values. The pressure distribution for a compressible fluid is governed by a nonlinear partial differential Reynolds equation. For the discretization of the film thickness variable, the finite difference method is used. With this method the Reynolds equation is substituted by a finite system of algebraic equations with the derivatives replaced by difference quotients. The solution of this system of equations is an approximation to the exact solution at a finite set of points. The bearing load capacity calculation is obtained by double integration of the pressure along the cone generator and for one period of vibration.

!

=

I

!

Fig 2 Support element geometry

Analysisproceaure The performance of this bearing, namely its working frequency and load capacity, is mainly a function of the geometry and material properties of the support member. Fig 2 presents the geometry that has been studied. The solution for systems with such geometry and boundary conditions can only be obtained using numerical analysis. In the finite element method, the continuum is divided into a discrete number of elements interconnected by the nodal points located at their boundaries. The supporting member has axial symmetry and also an axisymmetric loading, so the three-dimensional analysis is reduced to a bidimensional one, yielding a semi-analytical finite element formulation al . The finite elements are considered as solids of revolution with quadrilateral section. Once all element stiffness and mass matrices, and also the element force vectors are established, they are assembled to form the global stiffness matrix, the global mass matrix and the global force vector. The dynamic equilibrium equations are constructed using the Hamilton variational principle. IMI { ~ } + l K l { p ) =

{P}

(1)

The numerical integration of these equations can be performed either directly or using the modal analysis.

Fig I Squeeze film bearing model

80

The direct integration is a step-by-step procedure present in the equations without resorting to evaluation of the eigenvalues. It is very effective when the response for a few time steps is required. However, as in the present case, if the frequency of the external forces coincides with the

April 86 Vol 19 No 2

da Silva a n d A l m a s - t h e o r y o f squeeze f i l m bearings

system's first natural frequency and the other natural frequencies are not in the same range of the first, the response will be represented by the superposition of few natural modes of vibration and requires a long period of time to become steady. Furthermore, the use of techniques such as subspace iteration reduces the eigenproblem to the determination of the few required eigenvalues and eigenvectors. The mode superposition is, for these reasons, the most appropriate to the response analysis of the bearing support member. The finite element analysis requires some care in the mesh generation. Several hypotheses have been studied 12 . If the first natural frequency is compared with the experimental value fnEx = 18550 Hz, and with a theoretical value from the Ref 9, °ffnnG = 17700 Hz, it is possible to have a realistic image of the finite element performance. In fact, the best results are obtained by dividing the global geometry represented in Fig 2 into sections, each containing finite elements with uniform dimension and geometry. The section with greatest gradient of amplitudes of vibration, with respect to the cone, must have the smallest elements in the mesh. Also, the best results obtained correspond to those elements approaching the square shape or, when not possible, with the smallest dimension in the axis of the steepest gradient. Before presenting the best mesh configurations, the importance of the node numbering should be discussed. In fact, the matrix storage efficiency in the computer memory depends on the global matrix band width and the band width is a function of the maximum value computed for the difference between the highest and lowest node numbers in each finite element. Hence, proper numbering of the nodes reduces memory and computation time. The most simplified model with the mesh configuration represented in Fig 3 produced the result fnFE which is compared with fnEx and fnn G in Table 1.

Fig 3 Simplified mesh geometry simulated in this model. The amount of structural damping, associated with the oscillation of the bearing support member, may be represented by the loss factor and is estimated applying the half-power points technique to the experimental data 9 , in the neighbourhood of the first resonant frequencies. The result indicates the existence of light damping (r/='0.003), validating the procedure used. It also allows an equivalent viscous damping treatment, provided that the excitation is, as in this case, of the harmonic type. In these conditions the equivalent viscous damping coefficient is: 1 X = ~- 77 = 0.0015

I m p r o v e d models

The mesh in Fig 3 does not represent the whole geometry of the bearing member and it is necessary to extend the analysis to include the base and to ring in the outer edge of the cone as well as the piezoelectric ceramic element. The corresponding mesh is presented in Fig 4 with the first mode of vibration also plotted. The computed first natural frequency is then fnFE = 17404 Hz with a relative error of -6.2% about the experimental value. In order to obtain the dynamic response of the bearing member, a damping factor has been introduced, with the excitation created by the piezoelectric transducer also

Table 1 Comparison of relative error for various first natural frequencies

To obtain the value of the excitation force, one takes into consideration that the stress in the transducer/support member interface is a function of the voltage eo applied to the transducer electrodes, the transducer thickness t and the piezoelectric constant gzz. As shown in Fig 5, the interface area is then divided by the six nodes ( 6 1 - 7 5 ) into one solid circle ( l ) and five hollow circles ( 2 - 6 ) with internal radius rij and external radius rej. Assuming the stress as constant along the interface area, the amplitude of the resultant force F j and its position coordinate lj, are derived for each circle j as follows:

eo 2zr frreJi/ rdr

Fj-

gzz t

eo n(r2e/-r~.)

=

First natural frequency, Hz

Relative error about fnE x

fnF E = 18664

+ 0.6%

Fj lj

eo 2zr

f re/r 2 dr t rij 2e°rr(r3ej-r;')

- --

gzz

=

TRIBOLOGY international

(2)

gzz t

fn EX = 18550 fnn G = 17700

'

- 4.58%

(3)

3gzz t

81

da Silva and Almas - theory of squeeze film bearings

4

3

/9 2 15

y

18

.~4

39

Jllh

~ m 52 56

60

57

59

73

66

80

83

63

64

65

67

74

81

84

82

85

=l

54

61

68

70

71

72

75

55

62

69

76

77

78

79

0

Fig 4 Complete mesh configuration and first vibration mode (not to scale) R j = F/

therefore 3

3

2

2

2 (re/- ri/) l[--

3

(4)

(re/- ri/)

The forces F/are transferred to the nodes in the boundaries of the respective circle ], resulting in a set of forces R/and S~ (Fig 6). From the static equilibrium equations:

82

r e j - lj

re/- ri/

Si= Fj

(5)

lj - rij re/.-

ri]

Therefore, the resultant force in each node is the sum of the Sj component from the internal circle and the R/+ 1

April 86 Vol 19 No 2

da Silva and Almas - theory o f squeeze f i l m bearings

1

61

2

68 3 70 4

71

5

72

3mm 13.0104m~ mmJ mm / 4.3735mm im-i~

~i ~

~,~

~--~

6

75

4.3735mm ~ ~

i

computation time for several ratios At/T. These results refer to just four cycles of vibration and a standard value of 100 is attributed to the amplitude and time consumed with At~T= 1/16. The results corresponding to other At/T ratios are expressed as a percentage of this standard value. From these results and considerations it is clear that the time interval A t/T = 1/120 shows the best compromise. The amplitudes of vibration at the nodes of the finite element mesh are obtained and relevant results are presented in Fig 7 (nodes 10, 19, 25 and 31). Load capacity evaluation

Fig 5 Support/piezoelectric element interface

The load capacity in this type of bearing is highly dependent on the degree of compressibility of the gas between the two surfaces. With the simplification of the male cone surface being assumed as stationary, the film thickness will also change sinusoidally with time. The two major variables governing the performance of these bearings have been identified 9 as the excursion ratio e and the squeeze number o defined as:

o = 121~B2co/(Pah2o )

S

ri

re

(6)

At low squeeze numbers, when the frequency of oscillation co is low or there is a large mean clearance ho, gas flow will occur throughout the film. Only two forces oppose the flow: the forces due to viscosity which are small, and the inertial effects which are negligible. Under these conditions, the pressure and the force generated are proportional to the squeeze velocity rather than the displacement and their profiles are nearly symmetrical. This means that the film action will be that of viscous damping and no significant load capacity will be generated. The other limiting situation occurs for high squeeze numbers (high frequency and large excursion ratio). The viscous forces become important and flow resistance will be so high as to introduce compressibility effects. Flow will occur only in a narrow edge region and in the limit, the process 10

//~

Fig 6 System of equivalent forces applied to interface component from the external one. This equivalent system of harmonic forces is assumed in this model as concentrated at the interface nodes. As referred to before, the modal analysis was selected for the oscillation amplitude evaluation and only the first four natural eigenvalues are considered. Actually, the others have a negligible contribution to the mode superposition. The numerical integration requires the selection of the time coordinate increment At, as a function of two factors. On one hand, it must be small enough to allow a good representation of the period corresponding to the fourth natural frequency, and also should take into account that the system is forced to vibrate in a lightly damped resonant condition. This means a long transient period of many cycles of oscillation before the steady state is achieved, and short At values imply a large computing time. This is undesirable mainly when dealing with the large first natural period. In Table 2, a comparison is made between amplitudes and

TR I BO L O G Y international

l/IT

E E

//~11

I

0

/ ~It

=- 1

"o

o. E
F

~ ~~xxk \\\

\

,,V,x,.,. \

l_ I

~x

~x~ \

/

///X

//

,,

..///// x /

~. \ \ ,X,

Node 10 .

"\.Xx,,""..19 " x . . X

/

0 16.6

\

~25

""X-..

31

"X'..

I

I

I

17

17.5

18

Frequency, (f) kHz

Fig 7Amplitudes of vibration

83

da Silva a n d A l m a s - t h e o r y o f squeeze f i l m bearings

becomes non-dissipative. The compressibility forces are proportional to displacement rather than velocity and the resulting force becomes more unsymmetrical. The air film now acts as a nonlinear spring. This state corresponds to the ideal performance of a squeeze film bearing with maximum load capacity and no energy dissipation in the film.

F=(X,T, ae, ga ~- ~ ) - Ha~

For cases of uniform excursion, the asymptotic theory provides an analytical solution to the Reynolds equation. However, when the excursion has to be considered as non-uniform along the bearing surface, the analytical extensions of this theory proved to be very limited on the allowable variation for this parameter. Numerical methods are an alternative procedure having the advantage of obtaining increased accuracy as a function of time step size and the most important advantage of being able to consider any case of modal shape of the bearing surface.

The working frequency of real bearings is imposed by operating limitations. Then, for most cases, the air film action is represented by a nonlinear spring with some viscous damping added. The edge effect will be responsible for a transient period of increased pressure until the interior pressure balances the edge flow forces and equilibrium is achieved.

A study of numerical procedures based upon finite differences was performed by Michael 3 . Explicit, semi-explicit and implicit difference schemes were examined from the point of view of truncation error, stability and computational efficiency. His results indicate the superiority of the Crank-Nicolson formula for the present application of a gas bearing problem in one space variable. Assume, therefore, the interval 10,1 I of the x-axis subdivided into N+I equal sub-intervals AX and the T-axis into intervals AT. Equation (11) can be approximated by a family of difference equations by setting:

The pressure in the film can be predicted by the wellknown Reynolds equation by making the following assumptions: • • • •

(12)

the flow is laminar and isothermal the inertia effects are negligible the fluid is a perfect gas with Newtonian behaviour there is no pressure change across the film thickness

For a compressible fluid, this equation can be written as: a

L (~P) = L/n (u)

(h3p- ~x

x =j ,~x

a a =6/~ {2 fit (ph) + U -~-~i (ph ) )

(7)

T = (n + O) A T

Two more assumptions are considered for mathematical simplification. The conical axisymmetry reduces the two space variables to one space variable (x coordinate axis along the cone generator direction), and for a pure squeeze film condition, there will be no tangential velocity of one surface relative to the other, ie U = 0. Equation (7) then reduces to:

3 (h3p_ ap a O~x )= 12~ ~ ( p h )

Ox-

*=0u;+l+(1-0)4 3.

aT a~

(13)

_ (u/n+l-u~)IAT

[0(U;+~ 1 -Uj_l n + l l + ( l - 0 l ( u T + l -U/_l) ]

ax

(2/,x)

O2q,

(8)

aX 2

In dimensionless form Eq (8) becomes:

- [0' (UT+l 1- 2"; +1

U/-ln+l"I)

+ (1 - 0') (u:+ 1 - 2u~ + U~_l) ]/(2 ~XX)2 a aP aX (Hap ~ ) =°

3 (PH) ~-~

(9) where u f is a function, representing the product PH, defined at the nodes:

where:

X = x/B

P = P/Pa

T = cot

H = h/ho

(jAX, nAT)(j=O,I

..... N+l;n=0,1 .... )

and B is equal to the film length along the cone generator.

with u~, UN+ n 1 replaced by appropriate boundary values

Using the usual transformation qs = PH in Eq (9):

and u? as initial values.

a

[,I,(H aq.,

q, a H

a~

The amplitudes of oscillation were computed for fifteen finite element nodes ( 4 - 4 6 ) dividing the bearing surface into fourteen equal space intervals (Fig 4). Then:

(10)

Eq (10) is a quasi-linear differential equation better represented as:

~I' - [F , ( X, T, ,,It, -a~~) L ( ~ ) = 3m-f 0'4'

+ F z ( X , T , qs, ~ )

0 2 ~I'

AX = 1/14 From Michael's evaluations, the propagated error in the Crank-Nicolson method is proportional to [(AX) = + ( A T ) 2 I, and it seems reasonable to equalize the AXand AT contributions, ie making AX = AT. Then:

+

~-X~-- ] = 0

(11)

AT-

where:

3~

Fl ( X , T , * , ~ - ) =

84

1 IH(O~ 2 O~ H _ ~ 2 o ~-) _a., ~ - OX

02H ~-[

27r

1

Nt

14

corresponding to 88 equal divisions Nt for one period of oscillation.

April 86 Vol 19 No 2

da Silva a n d A l m a s - t h e o r y o f squeeze f i l m bearings

With AX and A T defined, the u/n + 1 are determined from the u n] by solving a system of algebraic equations (Eqs (1 I) -(13)). The equations are nonlinear and must be solved by iteration using the Newton-Raphson method. Algorithms .for this method 3'9'12 will enable the calculation of the pressure distribution for each desired instant time T. To obtain the load capacity, a double integration for the space coordinate interval and for one period of oscillation, must be performed. Numerical integration was used and Simpson's rule was the selected method.

60 V 90 V

11 100

,,.is

'.,, '.,. A Z

\

10

The instantaneous load capacity Winst is given by:

\

'., ',,,. \

Winst = 2 n sin 2 aBZpaf2 (P - 1) X d X

(14)

\

or taking a numerical formulation: Winnst = 2~" sin 2 o~B2pa [ ~

\

\~, "\ ,\\

,\ 10' \

[ 4 ( P n - 1 ) &X

\

+2(P n-1)zAX+4(P

n-1)3AX+...

+ 4(pn3 - 1) 13 AX] ]

n = 0, 1. . . .

(15)

where a is the cone half-apex angle (a = 45 o). 0

I

The time average of Winst for one cycle of oscillation represents the bearing load capacity W:

I

lill

I

I

I

10

I

\

-

\ ~e20 I ~ t , l "\ 100

h 0 (P)

1

2~r

[4' = 2--~ fo

Winst d T

(16)

Fig 9 Load capacity for 20 study cases

(17)

Winst are presented in Fig 8, corresponding to cases 1,2, 4 and 7, clearly showing the non-symmetrical behaviour of Winst in respect of the time. A transient period of increasing load capacity, due to the edge effect, is also evident in the family of curves converging to the steady state.

using a numerical formulation: 1 (0) . (1) u/(2) 44W~3)t+ W= 3 × 88 (Winst+4Winst + 2 " i n s t ..

+...

. (87) + 4 Winst +

w(88), ,, inst )

Twenty cases were considered for the present conical bearing application: cases 1 - 1 0 (e = 60 V); cases 11-20 (e = 90 V), e (at finite element node 4) = (1,0.9, 0.8 . . . . . 0.1). Some examples of computed instantaneous load capacity

6OO

6OO

40O

400

200

20O

A

t~

2PI J

0

a

Conclusions 2PI

0

b 600

A

Z U.

400

200

200

fl~

2PI

c

,~,,,~,,,~..._

2PI

o]

d

Fig 8 Instantaneous load capacity." (a) case 1; (b) case 2; (c) case 4; (d) case 7

TR I B O L O G Y international

The bearing load capacity for the study cases referred to is shown in Fig 9. It can be seen that load points taken for the same excursion ratio and characterized by an index number differing by 10 (eg 1 and 11, 5 and 15, etc) have identical values. This shows that the excursion ratio has an effect that overshadows other parameters in influence. Regarding the input voltage to the ceramic element eo, it clearly also has an important effect on the bearing performance. In fact, an increase in eo will cause a similar increase in the amplitudes of the bearing surface, the final result being a higher mean film thickness ho for the same load capacity w (as shown in Fig 9).

The prediction of the resonant frequency of the bearing unit is of primary importance because the maximum load capacity will occur at or very near to this frequency. Furthermore, the investigation of the flexibility of the bearing surface is required for the prediction of the film thickness shape. The FEM enabled both objectives to be achieved, producing a better and more complete representation of the bearing than other numerical procedures previously used. A new way of improving squeeze film bearing analysis is now possible with this numerical method. Ahnost any design parameter can be changed and its effect studied. The final result will be an improvement in bearing performance through an optimization procedure. There are, however, some modelling deficiencies which in most cases derive from the need to simplify the analysis.

85

da Silva and Almas - theory o f squeeze f i l m bearings

Some have been referred to in the text but mention is due to some less obvious conditions that have been assumed in this study. The axial stiffness of the securing member screws was taken as infinite with the interfaces between crystal and base and the crystal and the centre electrode assumed as of uniform and permanent contact. This electrode is an experimental set-up between two similarly shaped bearings and is intended to provide the cancellation of axial reactions at the centre plane. Because true symmetric arrangements are difficult to ensure in real bearings, it is hard to assess its influence on the dynamic behaviour of the bearing, namely the frequencies of resonance. For a dynamic system, undergoing a forced vibration and at resonance the damping factor is critical. To determine this factor the best procedure is experimental. In this case a mixed technique was used, but no a posteriori confirmation of results was possible. To obtain actual amplitude values and because of the reduced size of the bearing compared to available accelerometers, optical sensors are currently used. The main drawback of this technique is that it requires the removal of the supported member. For pressure measurements, a similar situation arises and special devices are being studied that would cause minimum bearing performance disturbance. Summarizing the present situation for these bearings it is possible to say that both theoretical and experimental inaccuracies have limited application of these bearings in the past and prevented their successful industrial use. Their potential due to the method's working simplicity would otherwise seem promising in particular applications.

86

References 1.

Tipei N. Equatiile lubrificatiei cu gaze. Communicarile Acad. R.P. Romine, 4, 1954, 699

2.

LangloisW.E. Isothermal squeeze films. Appl. Math. X X (2), 1962 MichaelW.A. Approximate methods for time dependent gas film lubrication problems. IBM Res. Rep. R J-205, 1962 and J. AppL Mech. 30, 1963, 509

3.

4. 5. 6.

Salbu E.O.J. Compressible squeeze films and squeeze bearings. Tran~ ASME, (D], J. Basic Eng., 1964, 355-364 Pan C.H.T. Analysis, design and prototype development of squeeze film bearings for AB-5 gyro. MIT, Rep. 6 4 - T R - 6 6 , 1964 Pan C.H.T. On asymptotic analysis of gaseous squeeze film bearings. Trans. ASME, (F), Z Lubr. Technol., 89, 1967, 245- 253

7,

Crank J. and Nicolson P. A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type. Proc. Cambridge Philos. Soc., 43, 1947, 5 0 - 6 7

8.

Huxley A.S. Progress in squeeze film technology. Pap. 16, 5th Joint ONR/UK Gas Bearing Meet., 19 72 Pina da Silva F. Squeeze film air bearings with flexible supports. Ph.D. Thesis, Imperial College, 1980

9.

10. Goldberg J.E., Bogdanoff J.l. and Marcus L. On the calculation of the axisymmetric modes and frequencies of conical shells. Z Acoust. Soc. Am,, 32 (G), 1960, 738-742 11. Zienkiewicz O.C. The finite element method. McGraw-Hill, New York, 1977 12. Matos Almas E.J.A. Dynamic analysis of a conical squeeze film bearing. M, Sc. Thesis, Instituto Superior T~cnico, Lisbon, 1984

April 86 Vo119 No 2