NUMERICALMETHODS IN HYDRODYNAMIC LUBRICATION
J. H. Vohr
INTRODUCTION
Analysis of the performance of fluid-film bearings involves the solution of Reynolds
equation, derived previously in the chapter on Hydrodynamic Lubrication (Volume II). For
an incompressible lubricant, the steady state form of this equation is
(1)
While there are many analytical solutions to this equation for relatively simple geometries,
the usual method of solution is by numerical analysis, employing either finite difference or
finite element techniques.
1-3
The following provides a brief description of these procedures.
FINITE DIFFERENCE SOLUTIONS
Assuming constant lubricant viscosity, Equation 1 can be rewritten
(2)
The purpose of finite difference approximations is to reduce a continuous differential
equation, such as Equation 2, to a series of algebraic equations that can be solved for the
pressure distribution P(x,y) for a known distribution of h(x,y) and µ(x,y).
First, the distribution of all functions in all coordinates concerned are replaced by samplings
at several discrete points. Usually the spacing between points along each coordinate axis is
chosen to be constant, but uneven spacing may be more suitable for some special applications.
The governing differential equations are then written at each point of the coordinate grid
with suitable approximations for the functions and their derivatives. Most common of such
approximations are central difference formulas applicable to a uniform grid spacing as shown
in Figure 1.
(3)
(4)
(5)
(6)
(7)

imposed, Equation 14 may be rewritten
P(2,2) = 1/4 [P(3,2) + P(2,3) ] + 3/8 10/0.03 P(3,2) −6/4 1/(0.03)
3
(15)
Similarly, at the i = 2, j = 3 grid point
P(2,3) = 1/4 [P(3,3) + P(2,2) ] + 3/8 10/0.03 P(3,3) −6/4 1/(0.03)
3
(16)
and so forth for i = 3, j = 2 and i = 3, j = 3.
Starting with some assumed initial distribution for the interior pressures P(i,j) [P(i,j) =
0 everywhere is an adequate choice], one can solve Equation 15 for P(2,2). This “updated”
value for P(2,2) may now be used in Equation 16, together with the initially assumed zero
for the rest of the pressures, to obtain an “updated” value for P(2,3). One proceeds in this
fashion up the i = 3 grid columns, solving for P(3,2) and P(3,3), always using the most
recently updated values of surrounding pressures when solving for P(i,j). Successive sweeps
through the finite difference grid continue until successive interior pressures P(i,j) converge
to a final distribution. Typically, 30 sweeps may be required for convergence.
Several improved procedures of applying relaxation methods exist such as “overrelaxa-
tion”, which strives to accelerate the convergence process by multiplying the change in
pressure from one sweep to the next by a factor such as 1.2. A drawback of relaxation
procedures is their proneness to numerical instability. This phenomenon sets in depending
on the value of the coefficients, and of the first guess distribution.
As seen in the example above, the equation written for point (i,j) is solved for P(i,j).
Therefore, the value of P(i,j) is corrected at each iteration by the influence of its immediate
neighbors. As a consequence, the influence of each point propagates one grid interval each
iteration and the number of iterations involved in achieving steady-state solutions is pro-
portional to the number of points on the longest side of the grid. In addition, the number
of operations is proportional to the total number of grid points N × M. Therefore, the total
time needed for converging to a steady-state solution is proportional to N
2

q
2
+ … a
2L
q
L
= r
2
. . . . . . .
a
L1
q
1
+ a
L2
q
2
+ … a
LL
q
L
= r
L
(17)
where
96 CRC Handbook of Lubrication
93_104 4/11/06 12:09 PM Page 96
Copyright © 1983 CRC Press LLC
Obviously, most of the “a’s” are zero. The system (Equations 17) may be written as in
matrix notation as

and <P(i)>
j+1
. These equations can
be written in matrix form as
[ A ]
j
{ P }
j
+ [ B ]
j
{ P }
j–1
+ [ C ]
j
{ P }
j+1
= { R }
j
(20)
where [A]
j
, [B]
j
, and [C]
j
are matrices which are M × M in size, M being the total number
of grid points in the i direction, i.e., the “height” of each column vector <P> j. Matrices
[A]
j
, [B]

2
= [ E ]
3
{ P }
3
+ { F }
3
(22)
where
(23)
If next written at j = 3, Equation 20 will involve unknown column pressure <P>
2
,
Volume II 97
93_104 4/11/06 12:09 PM Page 97
Copyright © 1983 CRC Press LLC
<P>
3
, and <P>
4
. By using Equation 22 to eliminate <P>
2
from this equation, one can
obtain an equation containing only the unknown column pressures <P>
3
and <P>
4
. We
can write this equation formally like Equation 22, i.e.,
{ P }

+[ B ]
j
[E]
j
) { P }
j
= [ C ]
j
{ P }
j+1
+{ R }
j
− [ B ]
j
{ F }
j
(26)
Comparing Equation 26 with Equation 25 we can deduce the recursion relationship for
determining [E]
j
and <F>
j
.
[E]
j+1
= –[T]
j
[C]
j
{F}

a
> , use the recursion
relationships (Equation 27) to determine [E]
j
and <F>
j
for all values of j from 2 to
N.
2. At the end of the bearing (j = N), set <P>
N
equal to the known boundary pressure
<P
a
> .
3. Use Equation 25 successively from j = N to j = 2 to determine all column pressures
<P>
j
.
While the matrices which must be inverted here are of the order M × M, in the direct
inversion method the matrix is of the order M × N. If M is chosen to be significantly
smaller than N, computation cost of the columnwise method will be significantly less than
that of the direct method.
A more detailed description of the columnwise solution method, including a discussion
of cyclical boundary conditions, is given in Reference 1.
Calculation of Flow and Power Loss
Once the discretized pressure field is obtained, quantities such as bearing flow and power
loss may be obtained from approximate finite difference expressions. For example, lubricant
flow per unit width in the direction of sliding at the i,j grid point would be evaluated by
the following central difference expression
(28)

simplest form by solving for the pressure in a one-dimensional, step slider bearing following
the sequence presented by Reddi.
4
Reddi presents a proof that if a functional φ(P) is given by
(31)
integrated over the region R, comprising a lubricant film, then
δ
φ
(P) = 0(32)
if and only if Pis a solution of the incompressible lubrication problem. In Equation 31, µ,
h, and Pare the lubricant film viscosity, thickness, and pressure, respectively, is the
vector velocity of the bearing surface, and q represents the outflow flux of lubricant from
the region R in in.
3
/sec/in. (m
2
/sec/m) across the boundary C
2
.
In Equations 32 and 33 pressure Pis described by a set of functions which are continuous
over R and which satisfy prescribed boundary conditions. In the finite element method,
region R is broken up into small elements and the distribution of pressure Pover each
element is described by convenient interpolation functions particular to the element. For
example, triangular elements are convenient for two-dimensional solutions as shown in Figure
3a.
The most convenient interpolation function is a linear one, i.e., the pressure P
m
over the
m
th

1
and ᐉ
2
, respectively. Since P is a function of x only, the finite elements
are linear rather than triangular. The pressures of the nodes of element 1 are P
A
and P
B
while the pressures at the nodes of element 2 are P
B
and P
C
. Using a linear interpolation
function for the pressure distributions over elements 1 and 2
P
1
(x) = P
A
+ x/l
1
(P
B
– P
A
) (34)
P
2
(x′) = P
B
+ x′/l

2
= { 0 (1 −x′/l
2
)x′/l
2
} (38)
By differentiating Equation 36, one obtains the gradient of P
m
(x). For two-dimensional
problems, the gradient would be a column vector, i.e.,
(39)
The matrix formed from differentiating the row vector T
m
first with respect to x and then
with respect to y is denoted as R
m
. In our case, ∂T
m
/
∂
y is zero so R
m
reduces to a row
vector, where
(40)
Since we now have expressions for P and P over each element m, the integral over R in
Equation 31 can be written as a summation of integrals over each element m. Using the
general notation of Reference 4, Equation 31 becomes
(41)
In our simple one-dimensional case, the element property term C

102 CRC Handbook of Lubrication
93_104 4/11/06 12:09 PM Page 102
Copyright © 1983 CRC Press LLC
Substituting these into the expression for V, integrating over each element, and summing
(48)
Finally, for Q, flows q
1
and q
2
leaving the bearing per unit width at the leading and trailing
edge boundaries, respectively, are
(49)
Using expressions (Equation 38) for T
1
and T
2
, x = 0 at the leading edge boundary and
x′ = ᐉ
2
at the trailing edge:
(50)
Substituting these expressions into Equation 45 for Q, integrating along the leading and
trailing edge boundaries, and summing
(51)
Now, Equation 32 states that if the pressure P is a solution of the incompressible lubrication
problem, then δφ(P) = 0 over the region R
(52)
where P
i
are the unknown pressures at the nodes of the elements. Since δP

associated with a finite element solution. As such, it is intended to provide only a “flavor”
of how the method works. In this simple case, the answer could have been obtained more
easily by a less involved procedure, so our example does not illustrate the power of the
method. Numerical advantages the finite element method has over finite difference techniques
are discussed in Reference 4.
The numerical techniques described in the preceding pages may also be applied to solution
of the energy equation derived in the chapter on Hydrodynamic Lubrication (Volume II).
Application of the columnwise influence coefficient finite difference method to solving both
the pressure and energy equations for a tilting pad thrust bearing are given in Reference 8.
Application of the finite element method to solving for both the static and dynamic char-
acteristics of tilting pad and step journal bearings is provided in References 9 and 10.
REFERENCES
1. Castelli, V. and Shapiro, W., Improved method for numerical solutions of the general incompressible
fluid film lubrication problem, J. Lubr. Technol., Trans. ASME, 89(2), 211, April 1969.
2. Castelli, V. and Pirvics, J., Review of numerical methods in gas bearing film analysis, J. Lubr. Technol.,
Trans. ASME, 90, October 1968.
3. Castelli, V., Numerical methods, in MTI Gas Bearing Design Manual, Wilcock, D. F., Ed., Mechanical
Technology Inc., Latham, N.Y., 1972.
4. Reddi, M. M., Finite element solution of the incompressible lubrication problem, J. Lubr. Technol., Trans.
ASME, Ser. F, 91(3), 524, July 1969.
5. Booker, J. F. and Huebner, K. H., Application of finite elements to lubrication; an engineering approach,
J. Lubr. Technol., Trans. ASME, Ser. F, 24(4), 313, October 1972.
6. Allaire, P. E., Nicholas, J. C., and Gunter, E. J., Systems of finite elements for finite bearings, J.
Lubr. Technol., Trans. ASME, Ser. F, 99(2), 187, April 1977.
7. Hays, D. F., A variational approach to lubrication problems and the solution of the finite journal bearing,
J. Basic Eng., Trans. ASME, 81(1), 13, March 1959.
8. Castelli, V. and Malanoski, S. B., Method for solution of lubrication problems with temperature and
elasticity effects: application to sector, tilting-pad bearings, J. Lubr. Technol., Trans. ASME, 91, 634,
October 1969.
9. Nicholas, J. C., Gunter, E. J., and Allaire, P. E., Stiffness and damping coefficients for the five-pad

in Denver in which an entire grandstand with 21,000 seats slides back on hydrostatic water
films to convert the stadium from football to baseball.
BASIC THEORY
Almost any fluid may be used in a hydrostatic bearing but this chapter will discuss only
incompressible fluids (liquids) in laminar flow. Fluid compressibility and turbulence intro-
duce special problems in gas bearings and high speed applications. The literature on these
topics is abundant.
Load
Load is the product of pressure and bearing area:
w = k(P
p
– P
o
) A(1)
The constant k is less than 1 and is a function only of the bearing geometry. One major
method of analysis
2
uses a “load coefficient” (a
f
) for the same purpose. The constant may
be thought of as the efficiency of the design in using the recess pressure to support the load.
For the geometry in Figure 1 for instance, k is approximately equal to (a – b)
2
/a
2
and A
= a
2
. The load as a function of supply pressure for this bearing is therefore:
Volume II 105

be important. Occasionally it becomes a problem because most of the pump power is
converted to heat when the fluid pressure is expended in the bearing. This can be serious
in high-precision applications such as astronomical telescopes and machine tools.
In the hydrostatic steadyrest of a large lathe, for example, a bearing with a geometry
resembling Figure 1 supports a heavy steel forging near midspan. The dimensions are a =
100 mm (3.94 in.) and b = 25 mm (0.98 in.). Oil is supplied at a rate of 0.025 l/sec (0.4
gpm) and a supply pressure of 56 MPa (8000 psi) by a fixed-displacement pump with an
efficiency of 0.5. The oil has a specific heat of 890 J/l(3.2 Btu/gal). From Equation 4, the
input power is calculated in SI units as follows:
H = 0.025 (56 × 10
6
)/1000 (0.5) = 2800 watts
In customary U.S. units, the power calculation would be:
H = 0.4(231)(8000)/60(12)(550)(0.5) = 3.73 hp
If all this power is converted to heat in the small oil flow, the resulting temperature rise
of the fluid would be 69°C (124°F). With a reservoir temperature of 38°C (100°F) this would
result in an oil temperature of 107°C (224°F) leaving the bearing. This local flow of hot oil
can cause a “hot spot” on the object being supported. In the lathe application, much of the
heat is absorbed by the rotating forging and is distributed by conduction and radiation. In
stationary applications local heating can cause distortions which may interfere with supporting
the load.
It is possible to minimize the power consumption with a rotating load because friction
(shear) losses come down and pump losses go up as the film thickness increases. The simple
case of a flat thrust bearing is given in Reference 6. Amore general discussion of optimizing
for minimum loss is given in Reference 3, based upon the parameter S
H
developed by Rowe
et al.
4
In applications sensitive to bearing heating, either of the above analytical methods

and has been verified in experiments.
8
Integration of the pressure distri-
bution in this bearing gives the load as a function of pocket pressure:
(8)
Equation 3 written for this geometry gives the flow.
(9)
Arestrictor in the supply line forces the bearing flow to assume the characteristic of the
restrictor. Equations 5 to 7 result in the flow-vs load relations illustrated in Figure 4 since
the pocket pressure in these equations is a linear function of load (Equation 8). By equating
each restrictor flow equation with the flow equation for this bearing (Equation 9), load
capacity expressions can be derived as below.
Orifice Restriction
(10)
Constant Flow
(11)
Volume II 109
Copyright © 1983 CRC Press LLC
FIGURE 5. Effect of restrictor type on stiffness of a thrust bearing (Reference 8).
FIGURE 6. Double-acting thrust bearing.
Volume II 111
Copyright © 1983 CRC Press LLC
FIGURE 7. Bearing pad coefficients for annular thrust pad bearing. Annular recess is
centrally located within bearing width (R
1
+ R
4
= R
2
+ R

1
,/R
4
1/2 1/4
(R
3
– R
2
)/R
4
− R
1
) 0.60 0.73
Copyright © 1983 CRC Press LLC
Values of β(Equation 21 of Reference 2) are calculated over a range of loads for each
bearing. Performance coefficients for film thickness, pump power, stiffness and flow are
then taken from Figure 8. In a constant-flow bearing, pump power is proportional to load
so k
H
equals β. Film thickness coefficients over the range of loads in this example are as
follows:
Film thickness is calculated from the following relation from Reference 2:
The load capacity for each bearing is an inverse function of the film thickness cubed. When
the load increases to the value at which the pocket pressure reaches the supply pressure, the
calculated film thickness becomes indeterminate. This happens at 60.3 kN on bearing A
and 82 kN on bearing B. Pump power (neglecting efficiency) is calculated from:
To find the total power required, the calculated values must be divided by the pump efficiency
(of the order of 0.5 to 0.9 for most gear, vane, and piston-type pumps).
Combining the individual bearing performances curves gives the results shown in Figure
9 for the overall assembly. Figure 9a shows that the total usable load capacity is about +70

p
P
s
16.58 × 10
−6
W 12.2 × 10
−6
W
Bearing A Bearing B
W, N
β
k
h
β
k
h
82,000 — — 1.00 1.00
70,000 — — .85 1.06
60,300 1.00 1.00 .74 1.11
50,000 .83 1.07 .61 1.19
30,000 .50 1.26 .37 1.41
10,000 .17 1.80 .12 2.05
Copyright © 1983 CRC Press LLC