NASA Technical Reports Server (NTRS) 19920021478: Fast methods to numerically integrate the Reynolds equation for gas fluid films

28 Pages·0.82 MB·English

by  NASA Technical Reports Server (NTRS)

NASA Technical Memorandum 105415
/
Fast Methods to Numerically Integrate the
Reynolds Equation for Gas Fluid Films
Florin Dimofte
Lewis Research Center
Cleveland, Ohio
Prepared for the
STLE-ASME Joint Tribology Conference
St. Louis, Missouri, October 13-16, 1991
(_ASA-T',I-IO_41 b) FAST '._F_THOOS TO
NLIMFt_,ICALLY INTFGRAI£ IH_Z RFYNGLL.)S
63/37 0115601

FAST METHODS TO NUMERICALLY INTEGRATE THE REYNOLDS EQUATION FOR
GAS FLUID FILMS
Jt
Florin Dimofte
National Aeronautics and Space Administration
Lewis Research Center
Cleveland, Ohio 44135
SUMMARY
The alternating direction implicit (ADI) method is adopted, modified, and applied to the Reynolds
equation for thin, gas fluid films. An efficient code is developed to predict both the steady-state and
dynamic performance of an aerodynamic journal bearing. An alternative approach is shown for hybrid
journal gas bearings by using Liebmann's iterative solution (LIS) for elliptic, partial differential
equations.
The results are compared with known design criteria from experimental data. The developed
methods show good accuracy and very short computer running time in comparison with methods based
on an inverting of a matrix. The computer codes need a small amount of memory and can be run on
either personal computers or on mainframe systems.
INTRODUCTION
The Reynolds partial differential pressure equation for a compressible, thin, fluid film is nonlinear
and involves a great deal of numerical effort to solve. From the early sixties to the present time several
methods were developed for numerically solving this equation. Major contributions were made by
Castelli and Pirvics, Booy, Coleman, Castelli in Gross' book, and others (refs. 1 to 4). The columnwise
influence coefficients method introduced by Castelli was one of the best at that time. The present work is
focused basically on two methods (alternating direction implicit method (ADI) and Liebmann's iterative
solution (LIS)), which add important advantages to the numerical solution technique of the compressible
Reynolds equation.
The Reynolds equation in dimensionless form for a compressible-fluid-film journal bearing subject
to isothermal conditions is (ref. 5)
Op2 0 0p2 O(ph) O(ph) (z)
a (h3 ___) + (ha __..) = 2A O(opeh) + 2Az + i4fA
Both steady-stateand dynamic analysestocalculatethestaticand dynamic performanceofan
aerodynamic journalbearingareshown intheappendix A. Startingfrom equation(1),thefollowing
partialdifferentiaelquationaredevelopedinappendixA:
National Research Council-NASA Research Associate at Lewis Research Center.
(_)
![+_o
N
_Q__!+ ho[
(_)
I i_Q°cosO + O_s_inoO + 2QocosO + _0c29o osO
_h o As OQo 2fA /
+2--- +i Q1= - _(----_ _ _
_2 2Qo2/3 & _]
(4)
1 [°_Qo inO t_Q° _0._
02ho A, OQ0 2fA I
+2_ - +i
Q2 = - __ho[_s_2 - _cosO_ + 2Qo sinO + _sin2Ooz
_Z2 2Qo2/30z
Equation (2)isusedto perform thesteady-statecharacteristicwsh;ileequations(3)and (4)areused to
compute thedynamic performance. The numericalintegrationtechniquesfortheseequationsare
developedand discussedinherein.
SYMBOLS
aij, bij, cij,
coefficientosfgridnode eq.(6)
dip eij, rij
:a supplyorificeradius,m
B¢ symbolicdamping coefficienutsedinstabilitaynalysis(seeeq.(A-27))
dynamic damping coefficientN,.m-l.s -1
Bij
C journalbearingradialclearance,m
CD orificesupplydischargecoefficient
D journalbearingdiameter,m
d orificseupplyholediameter,m
2
e eccentricity (see fig. 1), m
maximumiterationsteperror(seeeq.(12))
emax
F load capacity (the resulting force of the pressure distribution), N
F F/(PaLD) dimensionless load capacity (eq. (18))
Fr,F t component of F along and perpendicular to the centerline_ respectively (see fig. 1)
Fro,Fro steady-state component of F
Fx,Fy dynamic component of F along x and y directions (see fig. 1)
f whirl frequency ratio, u/f_
fo unstable whirl frequency ratio, Uo/fl
Gm dimensionless supply flow part that is function of Pm/Ps (see eq. (B-11))
h dimensionless film thickness, _'/C
ti" film thickness, m
ho steady-state component of h
hm film thickness at the supply orifice (see fig. 13)
i -11/2 the imaginary unit
ii the unit vector along xi direction (see fig. 14)
K stiffness, N.m -1
K dimensionless stiffness (see eq. (21))
Kij dynamic stiffness coefficient
Kc,Kco symbolic stiffness coefficient used in stability analysis (see eqs. (A-27) and (A-29))
L . bearing length, m
l length vector part of the Cm boundary of integral (B-7) (see also fig. 14)
M rotor mass allocated to one bearing; for a symmetric rotor M is half of the rotor mass, kg
corresponding rotor mass, allocated to one bearing, required to make the bearing unstable, kg
M,N number of grid points in i and j directions
n unit vector outward, perpendicular to the Cm boundary (see nlj in fig. 14)
total number of supply orifices
nO
P pressure, Pa
P dimensionless pressure
Pa ambient pressure, Pa
boundary pressure at bearing edges
Pla' P2a
Pm pressure downstream of the supply system
PO steady-state component of the pressure p (see eq. (A-6))
Ps supply pressure, Pa
Pl' P2 perturbation components of pressure p (see eq. (A-6))
Q dimensionless pressure variable, (ph) z
Q0, QI, Qu steady-state and perturbation components of Q (see eq. (A-7))
dimensionless flow through the pocketed orifice supply restrictor system
qm
R journal bearing radius, m
P gas constant (see eq. (B-9)), J.kg-l.K -1
r,t coordinates (fig. 1)
ST threshold speed, f/(M c C/W) 1/2
T absolutetemperature(seeeq. (B-9)),K
t time,s
V, axialvelocity(shaftspeedalongzdirection)m,.s-I
W bearingload,N
x,y whirlamplitudeofthejournalcenterinbearing(seeeqs.(.4_-20a)nd (A-25))
fluid film coordinates (see eq. (A-l) and fig. 14)
Xl,3
Zij impedance for translatory motion
g axial coordinate parallel to rotor axis
-/ adiabatic exponent
Az increment in z direction
AO increment in 0 direction
pocketed orifice compensation factor, a2/dC
eccentricity ratio, e//R
eccentricity ratio under static load
%
dimensionless tangential whirl amplitude
%¢1
dimensionless radial whirl amplitude
0 angular coordinate from centerline (see fig. 1)
A bearing number (see eq. (A-2) and (B-5))
Az equivalent bearing number in z direction (see eq. (A-2))
As restrictor coefficient (see eq. (B-9))
A relaxation coefficient (see eq. (10))
# dynamic viscosity, N-s.m -2
v whirl frequency, rad.s -1
unstable whirl frequency
v0
T dimensionless time (see eq. (A- 1))
attitude angle, deg
attitude angle under static load, deg
n rotation frequency, rad-s -1
nl shaft rotation frequency, rad.s -1
f_2 bearing housing rotation frequency, rad.s -1
Subscripts:
grid index
r centerline direction (see fig. 1)
t perpendicular to center-line direction (see fig. 1)
x x-direction (direction of the static load on bearing, see fig. 1)
y y-direction perpendicular to x (see fig. 1)
Superscripts:
n old iteration step
n+l new iteration step
(n+ 1)/2 half step iteration between n and n + 1
NUMERICAL COMPUTATION
In order to integrate equations (2) to (4), a finite-difference numerical method is adopted. This
method can handle a plain, cylindrical, journal bearing of finite length (fig. 1) with a high degree of
accuracy. Using a rectangular grid with constant spacing between nodes in both directions, the partial
derivatives become
(5)
"_ ij Qld+12A
(e)
(7)
_Q) _- qij+l - 2Qid + Qi,j-1
(8)
-o:2zql : qi+l,j - 2qi,j -k qi-l,j
j
where
(9)
AO, Az = node spacing
6
At node (i_) of the grid, equation (2) has the general form
aidQ01_l, j + bidQ0i+t, j + cidQ0ij_l + didQ0ij+l + eijQ0ij = ri.i (10)
Thus, a linear system of M × N equations results if the grid has M and N points in the z and 0 direc-
tions, respectively. Unfortunately, the coefficients of equation (10) depend on the value of the variable Q
so that an iterative method has to be used to compute the steady-state component of Q.
First, the system of equations (10) must be solved. Then, since
Q0 _ p0 (11)
and all the steady-state characteristics, the equations (3) and (4) can be solved, using the same finite
difference method, to calculate the dynamic performance.
ALTERNATING DIRECTION IMPLICIT SCHEME
The linear system (10) can be efficiently solved using the alternating direction implicit (ADI) scheme
(ref. 6), which changes the linear system of equation (10) into a tri-diagonal linear system using a proper
iterative manner. The tri-diagonal system can be solved with Thomas' algorithm. In order to use the
ADI method, equation (10) was written as follows for the first and second half step of the iteration n + 1:
aij (Qi_lj) (n+l)/2 W ei_i (Qij) (n+1)/2 + bij (Qi+lj) (n+l)/2 __ rij __ cij (Qij-1) n - dij (iQj+l) n
(12)
and then
cij (_Q ij-lJ _n+l -t- el.i (Qij) n÷l _ di.i (_f-'__i5+1/_n+l = ri,j _ aij (Qi-lj) (n+1)/2 - bij (Qi+lS) (n+1)/2
(13)
This same scheme was applied to equations (3) and (4) to obtain the dynamic characteristics.
Using the ADI method, the computer running time for one case Cone eccentricity ratio) is reduced
from 1 hr to 1 to 5 min (using a 20×9 grid and a 386, 20-MHz personal computer). To decrease the total
time to reach solution convergence, equation ill) can be used alone.
The convergence is difficult for equation (2), (3), or (4) for bearing numbers (eq. (A-2)) greater than
10. Also, the convergence of the numerical solution is sometimes difficult when the stability of the
bearing is analyzed, especially at high bearing numbers. One possibility to improve the convergence is to
use a relaxation technique:
(14)
qn-t-1 __+ _Qn+l _{_ ( I -- )_) Qn
where
A=0--+ 1
A second possibility is to convert the dimensions of the grid from a large one (e.g., from
Q
NxM = 40x33 or 30x25 or 20x17), used to compute the steady-state component Q0 of the variable
(eq. (A-7)) and the derivatives of Q0 and h0, to a small one (e.g., 10x9) and then to compute the
dynamic components Q1 and Q2 of the variable Q.
LIEBMANN'S ITERATIVE SOLUTION
The discretization of the Reynolds equation results in a large linear system of equations that needs
to be solved. For example, a 20 by 9 grid involves 180 linear algebraic equations like equation (10).
There are a maximum of five unknown terms per line in equation (10). For a large-sized grid, many of
the terms will be zero. For such a sparse system, the elimination methods waste significant computer
memory to store zeros. To avoid these difficulties, the LIS method (ref. 7) can be used as an alternative
to ADI. In this technique equation (10) is expressed as
Rij - Aij Qi-lj - Bij Qi+lj - Cij Qij-1 - Dij Qij+l
(15)
Qij : EiJ
and solved iteratively for j = 1 to N and i = 1 to M. Because equation (10) is diagonally dominant
for most of the bearing running conditions, this procedure will eventually converge to a stable solution.
Note that the coefficients C, D, and E depend on the value of the variable Q from the previous
iteration. As with the ADI method, a relaxation technique can be used to improve the convergence of
this method. The relaxation coefficient, A, can be variable and under control of the maximum error for
each iterationstepcalculatedas
(16)
¢hnax-Max
for
(17)
i= 1--,M, j = 1-* N
One of the most difficult cases regarding the convergence is the hybrid, compressible-fluid, film
bearing. This case balances the flow through the supply restriction system with the flow through the
fluid film to calculate the downstream pressure in the supply port. This requires an additional iteration
loop. Also, the equation of the flow through a restriction system, like a pocketed orifice, is a combination
of four relations that is a function of the type of flow (choked or unchoked) and the direction of the flow
(in or out). The analysis of a hybrid, compressible fluid bearing is outlined in appendix B.