427
Paper XIV(i)
A new numerical technique for the analysis of lubricating films. Part I: Incompressible, isoviscous lubricant C.H.T. Pan, A. Perlman and W. Li
This paper describes a new algorithm, which produces high resolution results, for the numerical computation of pressure and flux profiles in a lubricating fluid-film. The new method avoids the use of elementary functions to approximate the pressure field. Natural functions are generated by performing a directional integration of the pressure-flux relationship along a mesh line in the discretized domain. Pressure profiles which accurately depict details can be produced with a coarse computational mesh. Illustrative examples are given for a flat slider and a journal bearing. 1
INTRODUCTION
A high resolution of film pressure is often
desired in engineering studies. Critical issues may be related to the determination of load capacity, the estimation of lubricant heating, the calculation of minimum surface separation, the analysis of elastohydrodynamic deformation, and the description of film rupture. Adaptations of the conventional Finite Difference Method (F.D.M.) and the Finite Element Method (F.E.M.) are common practices in studies of hydrodynamic lubrication (1,Z). Use of a fine computational mesh is usually mandatory if good numerical accuracy is desired. The proposed new method is capable of producing very accurate numerical results without being dependent on the use of fine computational meshes. By performing directional integration of the relationship between pressure gradient and flux locally, a high resolution algorithm is derived and given the name Local Partial Differential Equation Method (L.P.D.E.M.).
E
eccentricity ratio of journal bearing
e
azimuthal coordinate
Q
dimensionless film flux vector
'PX
x-component of dimensionless film flux
z
dimensionless time
-
2
THEORY
Conventional F.D.M. and F.E.M. are dependent on the use of elementary functions which may not be most suitable to describe the details of the desired fields. In the derivation of the new algorithm, emphasis is placed on a practical procedure to compute the numerical relationship between the pressure field and the flux components according to the theory of lubrication, which is presented as: - aH Divergence Statement: div Q t - = 0 (1) az
1.1 Notation
Flux Formula:
5
=
-H3 grad P t % ( 2 )
A, B
constants of local flux field
CsD
constants of continuous flux field
H
dimensionless film thickness
,I
mth moment integral of the nth reciprocal power of H
LID
length/diameter ratio of journal bearing
2 . 1 Directional Integration of the Local P.D.E.
P
dimensionless film pressure
V
unit vector along the sliding direction
"X
directional cosine between x and sliding
X
local coordinate in arbitrary direction; also, abscissa of Cartesian system
Y
ordinate of Cartesian system
2
axial coordinate
It is presumed that a resonable mesh size can be selected according to the knowledge of the field property of the squeeze term in Equation (1) and the geometrical parameter of the bearing (e.g. the lengthldiameter ratio of a journal bearing), so that the flux components are smoothly varying in the discretized domain. The primary approximation to be imposed is that the flux components are adequately represented by truncated Taylor expansions. Consider an arbitrarily oriented mesh line, on which the distance x is measured from a local
-
v
Where is the unit vector aimed along the sliding direction. The flux formula used here is retricted to an isoviscous film. If temperature variation is to be considered, then factors to account for velocity distortion due to viscosity variation would have to be added (3).
428
origin. The linear Taylor expansion of flux is: cpx=A+Bx
(3)
where A and B are local constants. The proposed pressure profile generating equation in the x-direction is obtained by eliminating 'px between Equations ( 2 ) and (3): dP = Vx(1/H2) - A (l/H3) - B (x/H3) (4) dx Since A and B are local constants, Equation ( 4 ) can be integrated in a closed form:
-
Subscript "k" refers to the local origin situated at the k-th mesh point and the functions ,I are the exact integral from the local origin: Im(x)
=
(xm/Hn) dx ik
(6)
Accurate numerical treatment of these integrals is an essential part of the new algorithm in its practical implementation. For the theoretical film profile of an eccentric journal bearing, it is possible to apply the method of Sommerfeld to calculate their precise values (4). However, to make available a more versatile procedure, which may be applied to any sectionally smooth film profile, it is proposed that the mesh interval be divided into four subintervals, and the approximation of a polygonal periphery to span the subintervals be assumed. These integrals can be thus replaced by simple algebraic expressions as shown in the APPENDIX. This approximation is consistent with the discretization scheme. The local constants are linked to the discretized P-field through the "exact finite integrals".
The local constants, A and B, are determined by the increments of film pressure from the local origin to either adjacent mesh node. Thus,
- +
I13
I13-*13 - +
I03
I03-IO3 -I03
3.1 Pressure and Fluxes at Mesh Nodes Consider a local rectangular domain centered at the intersection of two mutually perpendicular mesh lines. A central approximation of the divergence operator may be derived for the uniform rectangular mesh by using spatial derivatives of the flux fields given by Equation (3). Thus, in coordinates (x,y)
where "c" denotes values at the central point or the local origin. Upon substitution of Equation (7) for Bx and By, a discretized approximation of Equation (1) is established, involving the local origin and its four adjacent points. The five point algorithm presented above can be written for all internal mesh nodes to form a complete, linear, non-homogeneous, algebraic system to define the discretized P-field. The system is blockwise tridiagonal, and can be readily solved by standard efficient matrix methods. Upon solving the discretized P-field, the local constants, A and B, can be calculated according to Equations (7) at all internal mesh nodes. At a boundary node, Equation (7) can only be used to find B along the boundary. However, Equations (1) and (8) can then be used to find B across. Afterwards, Equation (5) is applied toward the interior to obtain an equation, which is used for the calculation of A across. The set of values of A at all mesh nodes is the discretized cpx-field. 3.2 Computation of Profiles
2.2 Calculation of Local Constants
+[:
Computation of flux and pressure profiles between adjacent mesh nodes. Algebraic details of the first step depend on the discretization scheme. In the following, the orthogonal rectangular mesh system will be used to obtain illustrative examples. (4)
Equation (3) does not assure continuity between Ak and Ak+l-Bk+l(fbX). Consequently, there is a corresponding amguity in the pressure profile as given by Equation (5) and two distinct values of dP/dx exist at every internal mesh node. In order to construct a smooth pressure profile, it is necessary to replace Equation ( 3 ) by
A
'pcorr=
r
1
where superscript "+" or "-" defines the upper limit of integration of Equation (6) to be k+l. 3 COMPUTATION PROCEDURE The required computation procedure of the lubricating film involves the following steps: ( 1 ) Computation of the pressure field at discretized mesh nodes. (2) Computation of flux components at internal mesh nodes. ( 3 ) Computation of the cross flux at boundary nodes.
+
C x
+ D.x2
(9)
The coefficients C and D are selected to ensure continuity of P and of 'pcorr. That is,
+
102vx
+ +
- 103Ak
-
113ck
+
- 123Dk
(11)
123 is to be computed by the general procedure previously described for other.,I The pressure at a point between mesh nodes can now be calculated as
(12)
429
4 EXAMPLES 4.1 Flat Slider
A preliminary trial of the new algorithm was ap-
plied to a flat inclined slider. This served as as a convenient model problem since its geometry naturally suggests the use of the formulas given in the APPENDIX. The flat slider considered has the following parameters: Length/Width Ratio = 1.0 Inlet/Exit Gap Ratio = 9.0 z Calculations were made with a coarse mesh ( 6 ~ 6 ) ~ and also with a fine mesh (18x18). An encouraging experience from this example Figure 1. Film Pressure of Finite Journal Brg is the relative insensitivity of the result on 16x8 L.P.D.E.M., L/D=l, ~=0.9 the mesh size. Since all nodes of the coarse mesh are repeated in the fine mesh, one can compare the computed pressure values at these nodes directly. The ratios of the computed film pressure using the coarse mesh to that using the fine mesh at all common mesh nodes are shown in Table 1. Indices "i" and "j" respectively mark nodes along and across sliding. The entrance of of the slider film is at (i=O) and the centerline is at (j=3). The biggest discrepancy is a modest 5.9%. It is of interest to note that the coarse mesh calculation onsistently yields a lower value of film pressure than that from the fine mesh calculation. The residual inaccuracy in film pressure exhibits a second order trend. Figure 2. Azimuthal Flux of Finite Journal Brg 16x8 L.P.D.E.M. , L/D=l, ~ 0 . 9 1.000 1.000 1.000 1.000 Table 1.
1.000 0.972 0.974 0.974
1.000 0.973 0.972 0.971
1.000 0.970 0.969 0.969
1.000 0.960 0.965 0.967
1.000 0.941 0.962 0.968
1.000 1.000 1.000 1.000
Mesh Size Sensitivity of L.P.D.E.M. Square Flat Slider, 9:l Gap Ratio
4.2 Finite Length Journal Bearing The journal bearing of finite length was used as another model to evaluate the new method for two dimensional calculations. Specific examples are for ~ 0 . 9 0 and L/D=l. Mesh sizes ranging from 8x8 to 16x8 were used. Results calculated with the 16x8 mesh using the L.P.D.E.M. are shown in Figures 1, 2, and 3. Results from coarser mesh calculations were substantially similar. Direct comparison of the centerline pressure profile is shown in Figure 4 between 8x8 and 16x8 computations. The profile of the 8x8 computation is remarkably accurate, not only for the level of the peak, but also for its overall shape. It is especially noteworthy that the highest mesh node film pressure of the 8x8 calcuation is only onethird of the peak pressure. The journal bearing was also calculated by the standard F.D.M. using various meshes. There is much more dependence on mesh size of both the level and the profile shape of film pressure. An interesting tell-tale symptom of computational pathology is the rather bizzare shape of the azimutahl flux as shown in Figure 5. An abrupt dip is seen on either side of the location of the minimum gap. Some improvement is gained by mesh refinement, but the anomalous tendency remains evident even in the 16x8 calculation. There is no trace of similar difficulty in L.P.D.E.M. results.
Figure 3. Axial Flux of Finite Journal Bearing 16x8 L.P.D.E.M., L/D=l, ~=0.9
'A
a 4
a
16x8
R I 0
t
I -4
-
Figure 4. Centerline Pressure, Finite Journal Brg, L.P.D.E.M., L/D 1, E = 0.9 Thus, numerical evidence strongly suggests that a coarse mesh calculation by L.P.D.E.M. can achieve a level of precision which is very close to the corresponding asymptotic limit. Using a conventional F.D.M., similar accuracy can be attained only by a combination of mesh refinement and extrapolation.
430
7 ACKNOWLEDGEMENT Partial support from the Digital Equipment Corporation is gratefully acknowledged. References /
z
(a) 8x8 Mesh
(b)
16x8 Mesh
Figure 5. Azimuthal Flux of Finite Journal Brg Calculated by F.D.M., L/D=l, ~ 0 . 9 5 CONCLUSIONS A new method for the numerical analysis of an incompressible lubricating film has been demonstrated. Based on the illustrative example of a flat slider of finite width, and a comparative study against the F.D.M., using a finite length journal bearing as the model problem, L.P.D.E.M. is seen to possess three important relative advantages. The same relative advantages are believed to hold against the F.E.M. since the same fundamental issues are involved. They are: (1) Accurate results can be obtained with a very coarse computational mesh. (2) Accurate flux fields are obtained. ( 3 ) Intermesh profiles can be generated.
(1) Raimondi, A.A. and Boyd, J. 'A Solution for for the Finite Journal Bearing and Its Application to Analysis and Design', Trans. ASLE, 1958, 1, 159-209. (2) Reddi, M.M. 'Finite Element Solution of the Incompressible Lubrication Problem', Trans. ASME. ser. F, 1969, 9 l , 524-533. ( 3 ) Dowson, D. and Hudson, J.D. 'Thermo-hydrodynamic analysis of the infinite sliderI. The plane-inclined sliderbearing, bearing', Instn. Mech. Engrs., Proc. of the Lubr. and Wear Group, 1964, 34-44. (4) Somerfeld, A. 'Zur hydrodynamischen Theorie der Schmiermittelreibung', Z. Math. Phys., 1904, 50, 97-155.
APPENDIX There are four combinations of (m,n) needed of the integral ,I as defined by Equation (6). The closed form algebraic expressions of them based on the subdivided secant approximation of the film thickness profile are derived as follows. The secant approximation of the film thickness profiles allows one to write dx = (6x/6H) dH
6x and 6H are increments of the local coordinate and the film thickness of the secant segments. It can be shown that, with the lower limit of integration placed at the vertex x of the secant segment, integrated formulas of Eqaution (6) are
6 PROSPECTS OF FURTHER DEVELOPMENTS Although present examples are based on the field equations which are restricted to films of uniform viscosity, generalization to treat the effects of cross-film viscosity variation is a matter of rounding out details. The same can be said about extending L.P.D.E.M. to be used with non-rectangular mesh setups. The high resolution capability of the new algorithm opens up the possibility of performing fluid film analysis on small computers. In all likelihood, a 10x10 mesh would be very adequate for most problems. Since accurate pressure gradient is inherent in the high resolution result, L.P.D.E.M. will enhance theoretical studies of film rupture in its various forms. Because the flux fields can be accurately calculated with ease by L.P.D.E.M., it will be a powerful adjunct to numerical studies of thermohydrodynamic problems, where convective heat transfer is one of the dominating factor.
(A.1)
where