Computation methods for hydrodynamic problems
(Reynolds' equation)
G C Singhal This paper deals with the application o f numerical methods to finding approximate solutions for a two
The basic equation of hydrodynamic lubrication developed by Reynolds is a differential equation giving the relationship between the pressure gradient, the viscosity of lubricant, form of oil film, clearance and speed. There are two possible methods for the approximate solution of Reynolds' equation for a finite journal bearing: an analytic method, where the differential equation may be solved by a truncated series solution; and a numerical method, where the Reynolds' equation is replaced by a finite difference equation. The advantage of the series method is that the solution may be obtained quickly. The advantage of the numerical method is its versatility: the method can be modified easily to incorporate alternative boundary conditions.
NOTATION bs
de/dr
bW
1 - 2d¢/dt
2N' the squeeze film effect 2/V~ the wedge film effect
E
e/c, the eccentricity ratio
e
Eccentricity of the journal relative to the bearing, in cm
H
h/c, the nondimensional oil film thickness
K
Number of iterations
L
Length of the bearing in the axial direction (z-direction), in cm
m,n
Number of mesh points in the circumferential and axial directions respectively
P*
p H 3/2 nondimensional height parameter
p
Pressuregenerated in the oil film, in kg/cm 2 p ( c / R ) 2 / N ' nondimensional pressure
R
Radius of the journal
# ~2
Absolute viscosity of the oil, in kg-s/cm 2 Over-relaxation factor
FINITE DIFFERENCE FORM OF REYNOLDS' EQUATION To solve the two-dimensional Reynolds' equation by a finite difference technique, the pressure is evaluated at the gridpoints of a rectangular mesh which covers the developed surface of the bearing. A series of m points in the circumferential direction and n points in the axial direction, gives an (m x n) mesh. The position of any point on this mesh is described by two integer coordinates (i,j). These represent the appropriate mesh line number in the x and z directions, counting respectively from the line of maximum film thickness and from one end of bearing (Figure I). Although the pressure height parameter (P*) curves are smoother than the pressure ~ ) curves in the circumferential direction, at high eccentricity ratios, the higher derivatives are greatest around the point of minimum film thickness. Therefore, to economize on computing time and storage, a mesh graded in the circumferential direction is used. Let the mesh length in the circumferential direction be a function of its position designated 2uQ, where A x i = x i - x i _ I . The general Reynolds' equation in terms of the pressure height parameter (P*) is
c
Radial clearance in the bearing, in cm
a2p
D
Journal diameter, in cm
a
+ (2/rR)2L__aE-a2P * ~-
3 7r~
e e[~- sin 2 ~r2 -2cos 2 1 r ~ ] P *
48 7r3 - H3a
Mahindra Owen Limited, 148 Bombay-Poona Road, Pimpri, Poona 411 0~8, India
volume 13 number 3 may 1 9 8 1
[-bwesin21rx-+2bscos2~-E ]
(I~
(See Singha[ ~.)
0010-4485/81/030151-04 $02.00 © 1981 IPC Business Press
151
i.l,m~ _..~ Y" )t'k"~X..~&~,
Substituting equations (4) and (5) in equation (I) and neglecting third and higher order derivatives gives: 1 1 2~TR 2 2 { -Pi*,j [A~/(A2/+I )+ (A~)~--~ (~---)
l
j=l
[kk\\\\\\\\\' ¢,\\\\\\\\\\\~
" ~
3 e +--n2e( • sin2 2~E-2 cos 2 ~.E)] 2H ,~
[\\\\\\\\\\\{\\\\\\\\\\\,]
1
1 ~-p* _ P~-I ,j ~.~i(A~i+A.~i+ 1 ) i+1 ,l Axi+ 1 (AEi+A~i+I) 2~R 2 1 +(Pi*,j-1 +Pi*,/+l) ( ~ ) ~2 (A2) ~ +
z n
--- I-- ~ - - + - - + - - f l - - - L - - - F _ +-_ F - _ ] _ . q -~, ~--Pd-- d--~d--fl-+f-1
48 ~3 - Ha a [ - b w e s i n 2 n ' Y + 2 b s C O S 2 7 r E
- -~- + - ~ - - + - f- 4 ' ~ ' + - b -4- - t - 4 J --r--~-P g,_~- ,,-~j +,+,~I- H -- +- -I --i-- --I--F u- #--F " ~-4-- I-~ ---~----~. I 2
I
--~, I
I
--r~--i--4 I I ',
I
i
3
21~.R
I~
]
(6)
This simplifies to: X
,,, = e;-1,/ [ ~/(~% P.*.
m
I
+ ~%+~) / DEN]
"
1
Figure 1. G e o m e t r y of nodal points in finite difference
+P*i+1 ,j [~'~i+1 (AEi + AE/+I ) / DEN]
egpression
2rrR 2
+ (Pi,%-1 + P i * , j + l ) ( - - ~ )
I
" 2~ I(A-E) DEh'
By Taylor's theorem: 12
aP, P~+l ,/'= Pi*,j+AXi+l (~--)i,j 1
a2P *
+ H31~ [ b w e sin 2 n 2 - 2 b s cos 2 n 3]
1
where
a3P *
+~(~x;. )~(~T) +~(Z~x;+~)~%--~--). i,j a4p *
1
1
1 2nR 2 DEN = (AEi" AEi+I )+ (A~)----~ • (--L--)
I,j
+ __(Axi+ 1 )4.( ~ _ ) . + higher order terms 24 %J and aP* P*i-l,j : P,,j . * . - A x/. (a~ - - ) ..
(2)
I,J
1
+- (zXx;) ~
2
1
a2P *
(a-U-)u-6
(zXx;) ~ (~-U-)
. + higher order terms
(3)
/,J
aZ p ,
i,j 1
-P*. .
-',-~
+
1
P*
~+
[ ~ J '~i+1 (AEi+AEi+I) AEi asP *
P'i+1 ,j
AEi+I
-P* .
1
(4)
1
(-a~-) : 2 l " //+ (P* .+P*. _)] ;,/ ,(,~)2 2(/,,2)2 /,i-~ i,~+~ 1 a4P * - --12 (Az)2 (a~ -) + higher order terms
152
(8)
1
At/: DE~ [ a% (a% + ~%++)] I A2i = DE~-N [
a+P *
+ (AEi+l)2] (a_~__)
Equation (I) is transformed into its equivalent finite difference form using equations (4) and (5). This gives the value of pressure height parameter (P*) in terms of the value at the neighbouring points.
1
)]
The terms containing third and higher order derivatives can be ignored (See Smalley and Lloyd 2.) Taking the mesh length in the z direction as constant, the form of the second derivative with respect to z is: O2P *
e
(.-.sin 2 2 a E - 2 c o s 2 n V ) H
Comparing equations (7) and (8) gives:
+-3 (~ %- a x~+~- ) (a-~-) I - I-2 [ (A2/)2-A-EiAxi+I
3he
P.*.=A i,j li P~_Ij+A2i P*i+I,j+A3i (P.*..+P.*.+.)+A4i i,j-i i,j
Eliminating the first derivatives and rearranging (2) and (3) gives:
(a-~-) =2
+
Influence coefficients
aaP *
a4P *
1
+ --24 (A xi)4 ( ~ ) .
(7)
A3i-DEN I
I
~'~i+I (A~i + ~XEi+I ) 2nR
[(-L;--)
2
(10)
I
~2(A7) ]
(11)
24
A 1 i - DEN H 3/2 lra ( b w t s i n 2 n E - 2 b s C O S 2 n Y )
(12)
All the 'influence coefficients' except A4i are independent of the righthand side of equation (1). For a steadily loaded bearing 1
A4i-
(5)
(9)
DEN
24 H 3/2
tr3 tsin 2 h E
(13)
and the other coefficients remain the same.
computer-aided design
NUMERICAL
METHODS
There are two methods available for solving the finite difference equivalent of Reynolds' equation: the direct and the iterative methods. The direct method of solution yields the exact answer in a finite number of steps (ignoring roundoff error). The boundary conditions must be known on a rectangle with sides parallel to the axes of the variables of separation. It is therefore not suited for the solutions with internal boundaries of the form p = aP/ax = o (This gives a curved boundary for the end of the pressure curve for a journal bearing). Hence an iterative method of solution is used for the numerical work.
ITERATIVE
Jacobi method
Each new pressure evaluated from the finite difference equation implied a change from the previous value at that point. To increase the rate of convergence the overrelaxation method increases the change at each point. An intermediate value of the dependent variable is calculated as in the Gauss--Seidel method and the difference between this intermediate value and the current value is found and multiplied by a factor, usually between 1 and 2, before adding it to the current value to give the new value. This may be represented symbolically as:
,i
p, K i - l , j +A~ p . K +A ( _*K +peK I .'i i + l , j "~3i'Pi, j-1 -gj+l /
+ A 4i
(14)
The procedure consists of assuming convenient initial pressure values for all the nodal points and evaluating the pressure value at each internal point in turn according to (14). The iterative procedure starts from the point (1) and proceeds in a systematic manner by varying i from I to m - I, keepingj constant. This procedure is repeated f o r j varying from I to n - I . The procedure is continued till it produces significant improvement in values. Initial values of the dependent variable at all the nodal points can be assumed arbitrarily or may be started with a coarse net as illustrated by Scarborough 3. The final values do not depend on the initial values chosen, as the system is self-correcting. The number of iterations for a chosen convergence will be higher if the initial approximation of values is poor. If the new values at the adjacent nodes are used the rate of convergence is faster.
Gauss-Seidel method If, in each iteration, updated values are used the procedure is known as Gauss-Seidel iteration. The new value of the variable at each node is computed from the most recent value available at the adjacent nodes. This procedure converges at twice the rate of the Jacobi iteration 2. The scheme may be represented symbolically as p.K+I_A p . K + I + ,~ if -"li" i-l,j -"2i
p.K + l ij =p;K+~ K+I+
+ A 3i (P/:j-1
[Ali
p . K + l +A2~ p . K i-l,/ i i+l,j
p.K
i,j+1)+A4i
~
.
(16)
P'*j]
= P.*, u (1-~) + ~Q where D is the value from the Gauss-Seidel method
O p t i m u m over-relaxation factor
The pressure P,. i is a function of constants and four surround"rig pressure values, so the Jacobi scheme may be represented symbolically as:
U'"
Successive over-relaxation method
METHODS
These are the effective methods for numerical solution of the finite difference form of Reynolds' equation. These methods will be discussed and the importance of the Gauss-Seidel method with optimum over-relaxation factor outlined.
p,K+I=A"
becomes excessive. In many cases the time required can be greatly reduced by the use of successive over relaxation.
p.K
,"p.K + I i+1,f~zI3i ( i,j-1
+ .K PT,j * I ) + A 4 i
(1.5)
The convergence rate of the Gauss-Seidel scheme depends on the number of mesh points: if the number of mesh points is large, the time required to obtain a desired accuracy
volume 13 number 3 may 1981
The over-relaxation factor g2 is defined as the increased change of the value of the dependent variable between two successive iterations "divided by the original change given by the simple Gauss-Seidel method. In general, over-relaxation speeds up the convergence, but there exists a unique over-relaxation factor for any given particular system which causes the fastest rate of convergence. One of the most useful methods of determining the over-relaxation factor is that of Carre4, which calculates the over-relaxation factor during the iteration and proceeds towards the optimum factor. The optimum value for the over-relaxation can be obtained from:
1 - (1 - # 1 2 ) 1/2 ~opt = 2 [
(17)
]
/112
[4 +
]
#1 =1 . . . . 2 [m 2 + ( 1 r D n / L ) 2
The optimum over-relaxation factor is thus given as a function of length and radius ratio and the number of mesh points in the circumferential and axial directions. As the number of mesh points tends to infinity, '~opt tends to 2.
Table 1. Optimum over-relaxation factors usingapproximate formulae and Carre's method m
n
L/R
Carre's
Equation(1 7)
0.5 0.5 0.5 0.5 0.5 0.5 1.5 1.5 1.5
1.42 1.48 1.56 1.65 1.65 1.62 1.47 1.62 1.74
1.45 1.48 1.56
24 48 24 48 96
8 8 8 8 16 16 8 8 8
1.67 1.68 1.48 1.60 1.75
192 24 48
8 16 16
1.5 1.5 1.5
1.76 1.65 1.69
1.86 1.66 1.70
24 48 96
192
1.69
153
Table 1 provides information on the choice of the optimum over-relaxation factors. The value obtained by Carre's method shows that it is more influenced by the number of mesh points. If successive over-relaxation is used with the optimum over-relaxation factor, the saving in computer time compared with the Gauss-Seidel method is approximately a factor of 3.
CONVERGENCE CRITERION The iterations are continued until the following condition is satisfied for all nodal points (i,j)in a partivular iteration: I
i
P is some suitable small fraction known as the convergence factor, and the superscripts K and K - 1 refer to the Kth and (K-1)th iterations. This criterion is not used at points on the lines ~ = 0 and ~ = 1, as the pressures here are zero.
CONCLUSIONS N.umerical methods are preferred over analytical methods for finding approximate solutions of the hydrodynamic equation because of the versatility of these methods in adapting to different boundary conditions. Among the various numerical methods the Gauss-Seidel method with optimum over-relaxation factor yields the fastest rate of convergence and consumes about 1/6th of the computer time used by the Jacobi method. For practical purposes, Carre's method of finding the optimum over-relaxation factor is one of the most useful. As the number of mesh points tends to infinity, ~oot tends to 2.
154
For this work, the CDC 3600 computer installed at the Tata Institute of Fundamental Research is used.
REFERENCES 1 Singhal, G C 'Analysis of hydrodynamic journal bearings with axial and circumferential grooves' Dissertation Report (1973) pp 44-47 2 Smalley, A J and Lloyd, T 'Steadily loaded journal bearings' Conference on lubrication wear V-180, Proc. Instn Mech.Engrs Pt 3K (1966) 3 Scarborough, J B On numerical mathematical analysis Oxford and IHN Publishing Company (1966) 4 Carre, B A 'The determination of the optimum accelerating factor for successive over-relaxation' Comput. J. Vol 4 (1962)
BI BLIOG RAPHY Cameron, A Principles of lubrication Longmans, London (1966) Coleman, R L and Lathan, N Y 'A brief comparison of accuracy of time dependent integration schemes, for Reynolds' equation' Mech. Technol. Int. TR59 (1971) p42 Lloyd, T and Mccallion, H 'Recent development in fluid film lubrication' Conference on Lubrication Wear, Vo1182, Proc.lnstn. Mech.Engrs Pt 3A pp 36-50 (IME) Steva, G P 'Finite element analysis of hydrodynamic bearing' Aerotech Magazine (May 1974) p 41 Wada, S and Hayashi, H 'Application of finite element method to hydrodynamic lubrication problems' ISME Bulletin (1977) pp 1234-1244
computer-aided design