A Simplified Mass Conserving Algorithm for Journal Bearing under Large Dynamic Loads *

A simplified mass conserving solution approach, consisting of analytical and numerical methods, for performance evaluation of dynamically loaded journal bearings is presented. The analytical technique is used to determine position and velocity ofjournal center for given force components. Subsequently a finite difference formulation of universal Reynolds equation is used to calculate realistic oil flow. The proposed formulation is applied for analysis of an engine main bearing. The entry and exit flow and maximum pressure in the bearing are determined over complete cycle and matched with published results obtained by numerical scheme. The suggested hybrid computational method typically takes 55 s on ICL DRS 6000, and 29 s on 150 MHz Pentium-Pro computer.


INTRODUCTION
Journal bearings, used in reciprocating and rotary compressors, and internal combustion engines are subjected to large fluctuating loads.For example, a rotor in a rotary compressor for air conditioner whirls largely by unbalanced forces due to the eccentric rotation of the rolling piston, and by gas forces induced by the difference in pressure between compression and suction gases.Under large dynamic forces the rotor vibrates severely and the amplitude of vibration approaches the bearing radial clearance.The rotor-bearing system under such severe conditions becomes non-linear due to large variation in the stiffness and damping coeffi- cients of the oil film.This system can be studied effectively only by accounting for the non-linear forces produced by fluid film journal bearings.Such non-linear dynamics problem must be ana- lyzed sequentially by integration of the Reynolds equation at various discrete time steps.
A comprehensive range of solution schemes has been developed to evaluate the performance of bearing under large dynamic loads.These solutions can be categorized as rigorous and rapid methods.The rigorous methods involve oil feed features, transient evaluation of cavitation boundaries in lubricant film, etc. to predict minimum film thick- ness, maximum pressure, oil flow and power loss accurately if only massive computational facilities are made available.On the other hand, rapid solu- tion techniques are easy to use and predict fairly accurate results of minimum film thickness and maximum pressure.Since these parameters have profound effects on lubricant flow and power dis- sipation, a solution scheme is required which takes less computation time on desktop computer and predicts bearing performance accurately.
Rapid design tools for analysis of dynamically loaded journal bearings are established analytically (Booker, 1971;Ritchie, 1975;Hirani, 1998), graphi- cally (Booker, 1965) and curve fit forms (Moes  and Bosma, 1981; Goenka, 1984).Solutions can be obtained within a few seconds from an interac- tive computer terminal, by entering the initial data describing the application.Most of these methods are based on Gumbel and Reynolds cavitation boundary conditions.The load capacity and maxi- mum pressure can be predicted with reasonable accuracy using these boundary conditions.However, Gumbel boundary condition does not satisfy the continuity of flow at film rupture as well as reformation, so cannot predict the start and end of pressure curve accurately.Similarly Reynolds cavitation boundary condition represents rupture interface in a reasonable way but does not conserve when the film is re-established.Therefore the lub- ricant flow and to some extent power loss cannot be accurately determined using Reynolds boundary condition.
The collective work of Jakobsson-Floberg and Olsson, often referred as JFO theory (Dowson and  Taylor, 1974), has provided boundary conditions which conserve the mass at the interface between the full film and the cavitation zones at rupture and reformation boundary.This condition is consis- tent with conservation of mass and provides better insight of the phenomenon (oil-streamers in cavita- tion region).Therefore the performance of bearing can be predicted more precisely using JFO theory.The disadvantage of this condition is difficulty to accommodate in computer codes because complication arises in numerical treatment of the non- linear boundary conditions.Elrod (1981) developed a finite difference algo- rithm that incorporate JFO theory in a very simple manner.This computational scheme, known as Elrod algorithm or mass conserving algorithm is valid for complete geometry of the clearance space in shaft-bearing system.This algorithm automati- cally handles cavitation regions and operates satis- factorily with or without cavitation being present.If cavitation pressure is reached in a given application, the algorithm automatically introduces a cavity conforming to the requirements of mass con- tinuity, i.e.JFO boundary conditions.In Fig. the oil flow predicted using Gumbel, Reynolds and JFO theory is compared with the experimental flows for three arrangements.This indicates neces- sity of the application of mass conserving algorithm for dynamically loaded journal bearings for real- istic prediction of oil flow.
The analysis of a dynamically loaded journal bearing follows an altogether different approach, and is more complex and extensive compared to that for the case of steadily loaded journal bearing.In the analysis of steadily loaded fluid bearing the film geometry is generally specified or assumed.The fluid pressure consequent to that geometry is then integrated to find the load capacity of bearing.For the analysis of dynamically loaded bearing the calculation must proceed in the reverse direction, i.e. for a known load pattern, and the film geometry (its rate of change) is to be evaluated, which makes it an 'inverse problem'.Implementation of mass conserving algorithm for indirect problem requires a number of iterations at every time step.Moreover, the algorithm often diverge from true solution, and requires a time step splitting to remove the instability.This further increases the computational time significantly, and often requires the extensive computer facilities such as the super computer.Therefore due to high computational cost for inverse problem, the mass conserving algorithm is not very popular with design engineers or engineers in the field.
To determine the oil flow in dynamically loaded bearing within a few seconds various researchers (Martin and Lee, 1983; Martin and Stanojevic,  1990; Conway-Jones et al., 1991; Boedo and Booker, 1990) have suggested different curve fit methods.But there is lack of agreement between these curve fit equations and are deficient to determine oil flow accurately.Therefore a solution scheme is required which is rapid on the desktop computer and also predicts bearing performance accurately.Kumar and Booker (1991) showed that mass conservation algorithm is unimportant in determin- ing journal motion except probably when the load applied to the bearing is supported by an incom- plete film.Vincent et al. (1996) worked on mass conserving boundary conditions and predicted sim- ilar trajectory to the trajectory obtained with the Reynolds boundary condition.
To overcome the drawbacks pointed out, a hybrid solution approach for performance evalu- ation of dynamically loaded journal bearings is represented in this paper.This contains the ana- lytical and numerical approaches.The eccentricity ratio and journal center velocity for known force components are calculated using an analytical approach (Hirani, 1998).For known eccentricity ratio and velocity an Alternating Direction Implicit (ADI) fast-converging numerical scheme is adopted to solve universal Reynolds equation.
This formulation is applied to analyze a main crankshaft bearing.The proposed hybrid scheme predicts results comparable to those obtained by complete implicit numerical method (Paranjpe and Goenka, 1990).
The positive portion of pressure (P) lies between angle 01 and 02 that can be determined as (2) The journal center velocity can be evaluated at any time step as (Hirani, 1998) I3 cos 05 +/2 sin 0 { }_F(C/R) The position of the journal center at next time step is obtained as ct+At } t t (4)

NUMERICAL FORMULATION
The pressure profile of the fluid film within the clearance is governed by Reynolds equation.The two-dimensional transient form of Reynolds equa- tion (Vijayaraghavan and Keith, 1989) This is an expression of the conservation of mass: the first term on the left is the squeeze term; the second is the net mass flow rate in the x (flow) direction; and third is the net mass flow rate in the z (axial) direction.The flow rate in the x direction consists ofcontributions due to shear (Couette flow) and pressure gradient (Poiseuille flow); whereas the flow rate in the z direction is solely due to the pressure.Elrod and Adams (1974) assume compres- sible flow and modified by introducing a non- dimensional density (a= p/pc) term in the above equation to modify for compressible (due to the presence of air/gas in liquid) fluid as where Pc is the density of lubricant at cavitation pressure.It should be noted that c actually has a dual meaning because it may be interpreted as the ratio of densities in the full film region, while it denotes the fractional film content in the cavitation region, such as P in full film region (c >_ 1), pc ch in cavitated region (c < 1).
In addition, Elrod and Adams (1974) employ a switch function, which is defined as in full film region, gin cavitation region.
The switch function permits pressure terms to be retained (in full film) or neglected (in cavitation region).The pressure, non-dimensional density and the switch function are related through bulk mod- ulus as OP OP g p --p O .Oc Integration of above equation yields P Pc + g(log c). (6) Within the cavitation region the pressure is assumed to be uniform and equal to Pc. Since/5 is large (-500 MPa), c is slightly greater than unity in the full film zone, therefore Eq. ( 6) may be approximated as P Pc + g(a-1). (7) Introducing the switch function (g) and replacing the pressure in terms of a, Eq. ( 5) can be rewritten as 0 (oh)+ 0 oz U g/3h 0 x/h +zz 12# -0.
Equation ( 8) is an elliptical partial differential equation.In the cavitated region (g-0), Eq. ( 8 In the full film region, pressure increases gradually to a maximum value and then reduces.This behavior is governed by Eq. ( 5) for g l.There is continuity in the film geometry and all neighboring points influence parameter at any node point.This requires a central finite differencing scheme.On the other hand, in cavitation region there is disconti- nuity in flow, abrupt change in film geometry and fluid is unable to signal same to the upstream flow.Therefore one-sided or upwind differencing should be used.Before switching over to numerical method it is preferable to non-dimensionalize the partial differential equation (Eq.( 5)).This makes the range solution parameters very limited and consequently convergence will be faster.Let 5) in non-dimensional form is (10) The numerical scheme is based on the modified Elrod algorithm (Vijayaraghavan and Brewe, 1992).
The convective (shear) flow term, which exists both in the full film and cavitated region compared to pressure induced flow terms that exist only in full film portion, is treated separately.

Shear Flow
For convenience let E 0.5cH.The shear flow term can be written as 00 00 Using finite differencing SE) aEi+ + bEi + cEi-1 2 To present flow correctly, it is necessary to intro- duce central difference form of Eq. ( 8) in the full film portion of the bearing, and a one-sided backward difference form in the cavitated region.
Vijayaraghavan and Keith (1989) suggest the fol- lowing scheme to switch over from central differen- cing to upwind differencing (Forsythe and Wason, 1964) on change of full film to cavitation region: SE) --{ (g/+l --gi)Ei+l + [4 (gi+l -+-2gi-+-gi-l)]Ei [4 (gi q-gi-1 )]Ei-I }. ( If all points are within the full film (g 1), Eq. ( 12) which is the desired central difference form.On the other hand, if all grid points are within cavitation region then SE) Ei-Ei-1 The above formulation is based on the implicit assumption that the surface velocity (U) is positive.
In dynamically loaded bearings, velocity may be positive at some instants and negative at others.
If the velocity is negative (reversed), then convective flow must be forward differenced relative to the original backward (for cavitation region) scheme, such as (Press et al., 1992) <0.
In such cases, modified formulation (Vijayaraghavan and Brewe, 1992) should be used Values of g on the half step may be computed as gizk /2 gi + gil
The expression for pressure-induced flow in the z direction can be written in a similar fashion.

Alternating Direction Implicit Scheme
To evaluate the c, hence pressure distribution, a time march solution of Eq. ( 8) is required.In present study, an ADI, which is based on the concept of time splitting is used.The idea is to divide each time step into two steps.In the first half, changes in the x (flow) direction are calculated based on the previous time step c values in the z (axial) direction.In the second sub-step of time, the recently calculated c distribution is used to evaluate the flow direction gradients while the axial changes are updated.The advantage of this ADI method is that each sub-step requires only the solution of a simple tridiagonal system of equations.Tridiagonal has non-zero elements only the diagonal plus and minus one column; therefore one does not require storing complete matrix, but only non-zero compo- nents as three vectors.A complete solution scheme at every time step includes an iteration cycle that consists of a sweep in the circumferential direction in the first half time step, followed by a sweep in the axial direction.An Euler implicit time differencing scheme is used for time-march solution.These are expressed as Euler Implicit Scheme O(H) (cH)/--/xi/2 (cH) -/x- for first half; O( H) for second half. nV2 Circumferential sweep fl H3g--z q-(Hot) -A-t-.
This equation can be rewritten at each node in a general form as --A-/2 --A{/2 I-A?-/2 di aioi__l,j q-bioQ,j -}-cioi+l,j where Hi- IuI bi [(2 gi+l/2,j gi-l/2,j)H,] +A-i/2 These equations can be represented by a cyclic tri- diagonal system (Press et al., 1992) and solved by tridiagonal solver using Sherman-Morrison for- mula (Press et al., 1992).In the above equation, mainly three parameters film thickness, cx, and g are unknown.The film thickness can be determined using analytical approach.The values of c are determined for assumed g which typically come from the previous time step, and then iterating g to adjust them according to new time step.
Axial Sweep During the next half time step, the axial direction is differenced giving (Hc,) Generally the journal and bearing are aligned, and the pressure is symmetric about the center plan of bearing.Therefore, to reduce computational efforts the calculations are only made over half the axial length of bearing.

COMPUTATIONAL PROCEDURE
A computer code was written to evaluate the vari- ous important parameters.A computational grid of 37 nodal points in the circumferential direction and 11 nodal points in the axial direction is selected.This computational grid size is chosen as a com- promise between time and accuracy.For a given load cycle (force vs. crank-angle), angular velocity and bearing geometrical parameters the necessary steps required for analytical and numerical anal- ysis of dynamically loaded journal bearing are as follows: (1) Assign initial values for c, q5,8, and cij, gi#.
Initially, all nodes are assumed to be in the full film, and all c (except at groove locations) and g are equal to 1.
(5) Determine Ii at various nodes using ADI finite difference scheme.Repeat steps 2-5, until last time step is reached.
If convergence is satisfied, then move to step 7 else assign el-CNDATA, 1-ONDATA and c.
To illustrate the efficiency of the proposed sim- plified mass conserving algorithm a main crank- shaft bearing (Paranjpe and Goenka, 1990) is analyzed.A comparison for the distribution of inflow (Fig. 2), out-flow (Fig. 3) and power loss (Fig. 4) is made with results obtained by complete numerical solution (COMJOB) of inverse prob- lem (Paranjpe and Goenka, 1990).The results of finite element method using Reynolds boundary condition are also plotted for comparison.The proposed simplified method takes 55s on ICL DRS 6000, and 29s on 150MHz Pentium-Pro computer.

CONCLUSIONS
A hybrid solution scheme for performance evalua- tion of bearings under large dyanmic loads is pre- sented.A main crankshaft bearing is analyzed.The results of in-and out-flow, and maximum pressure over load cycle are compared with published plots.
Average oil flow into the bearing is found to be equal to the average oil flow out of the bearing.Very short computational time and close match- ing of results with those obtained by numerical methods indicates the rapidness and accuracy of the ---o--FEJOB (Paranipe and Goenka,1990)   ---CaM JOB ............  proposed hybrid scheme.Therefore one can say that the suggested computational scheme combines the good features of analytical and numerical methods and can be used with confidence for design of dynamically loaded journal bearing.non-dimensional axial co-ordinate time, s time increment, s non-dimensional time increment eccentricity ratio -1 journal radial velocity, s non-dimensional journal radial velocity, attitude angle between line of action of load and line of centers, rad non-dimensional journal tangential velocity, absolute viscosity of lubricant, Pa s non-dimensional co-ordinate in circumferential direction, tad angular position of the start and end of the pressure curve, rad angular velocities of bearing and journal respectively, rad/s MASS CONSERVING ALGORITHM average angular velocity of journal and sleeve relative to load vector, 0.5(COj @ (')B) &L non-dimensional average angular velocity, v/3 derivative with respect to time FIGUREAverage oil flow from, various predictive techniques(Paranjpe and Goenka, 1990).
FIGURE 2 Oil inflow comparison.

FIGURE 4
FIGURE 4 Power loss predicted by present study andCOMJOB.