A numerical approach for the analysis of deformable journal bearings

This paper presents a numerical approach for the analysis of hydrodynamic radial journal bearings. The effect of shaft and housing elastic deformation on pressure distribution within oil film is investigated. An iterative algorithm that couples Reynolds equation with a plane finite elements structural model is solved. Temperature and pressure effects on viscosity are also included with the Vogel-Barus model. The deformed lubrication gap and the overall stress state were calculated. Numerical results are presented with reference to a typical journal bearing configuration at two different inlet oil temperatures. Obtained results show the great influence of elastic deformation of bearing components on oil pressure distribution, compared with results for ideally rigid components obtained by Raimondi and Boyd solution.


INTRODUCTION
ournal bearings are machine elements in which the applied force is entirely supported by an oil film pressure.They are used in many different engineering applications, for example as supports of rotating shafts.They are considered superior to roll bearings because of their higher load-bearing capacity, higher operating angular speed, lower cost and easier manufacturing.Furthermore, a proper design can assure very large service lives.The fluid-dynamic motion of oil film within the lubrication gap is described by the well-known Reynolds equation [1].Explicit analytical solutions for the pressure distribution can be obtained only for asymptotic configurations (e.g.infinitely short and long journal bearings) [1,2], while for other journal bearing configurations numerical solution of Reynolds equation is required.The early studies on journal bearing performance under different operating conditions, based on the numerical solution of Reynolds equation, date back to the fundamental work of Raimondi and Boyd (R&B) [3,4].They summarized their results in useful dimensionless charts ready for design, which are nowadays accepted also in code standards [5].Raimondi and Boyd analysis is based on some simplified assumptions: constant viscosity of oil film, independency of viscosity on pressure and finally the postulation of perfectly rigid components (shaft and support).Such assumptions, however, can be somewhat oversimplified, considering for example that the deformation of journal bearing components under the imposed oil film pressure is expected to produce a change in the real lubrication gap and thus a modification in the resultant pressure distribution.Moreover, also the assumption of constant viscosity and its independence on pressure should be critically reviewed, as it is experimentally known that viscosity varies, other than with temperature, also with pressure, as summarized by many constitutive models [6].It would be then of interest to investigate in a more detail the correlation existing between all the above-mentioned aspects and journal bearing performance and design.

J
In light of the above considerations, the present paper aims to present a general numerical approach to study the behavior of hydrodynamic radial journal bearings, by including in the analysis the effect of the aforementioned aspects.Attention will focus on the computation of pressure distribution as a function of both temperature variation within lubrication gap and on viscosity-to-pressure sensitivity (according to the Vogel-Barus constitutive model [6]), as well as on components elastic flexibility.An iterative algorithm using a finite difference scheme will be developed to solve the Reynolds equation, based on the deformed lubrication gap calculated by a structural finite elements (FE) model coupled with the hydrodynamic equation.The numerical approach will be able to compute the pressure distribution and the local stress field by including shaft and support elastic deformation.Results will clearly emphasize the strong influence of component flexibility on journal bearing performance, with a significant reduction of the oil pressure distribution peak caused by components deformation, compared to the case of perfectly rigid elements.

JOURNAL BEARING: BASIC CONCEPTS
typical configuration of a radial journal bearing under a vertical load (see Fig. 1) consists of a shaft rotating inside a fixed support (choke), where it is usually fitted a bush.The nominal radial clearance between shaft (diameter d=2r) and choke (diameter D=2R) is c=R-r.The steady-state response of a journal bearing is governed by the fundamental equation of lubrication theory (Reynolds equation) [1]: where is the oil film thickness as a function of angular coordinate θ, symbol e is the eccentricity, U=ωr is the tangential velocity of shaft, ω is its angular velocity, p(θ) is the resultant oil pressure distribution over angle β (attitude angle), μ is the oil dynamic viscosity.The numerical solution of Reynolds equation gives the pressure distribution p(θ) within the lubrication gap together with other system operating parameters (eccentricity, minimum lubrication gap, force resultant components, etc.), as summarized for example in R&B diagrams [3,4].

A
During service, due to the relative velocity between shaft and support, the oil generates a pressure p(θ) over the attitude angle β, where pmax is the peak pressure that occurs at angle θpmax.The system moves in a new equilibrium configuration, where the eccentricity e characterizes the position of shaft axis with respect to the axis of fixed support, along the direction defined by the angle θ h0 (which also specify the direction of minimum oil thickness h 0 ).Several design charts are available in literature [3,4], which provide journal bearing operation parameters as a function of Sommerfeld number S=(r/c) 2 (μN/p m ), defined in terms of shaft radius r and rotational speed N, while p m =F/(LD) is the average (specific) pressure defined as the ratio of the applied radial force F and the nominal projected area (L is the length of journal bearing).Such diagrams were determined by R&B through numerical solution of Reynolds equation under the hypothesis of constant temperature (and thus viscosity) of lubrication film and also under the assumption of perfectly rigid components (shaft and support).An improvement of the analysis can be obtained by including in the Reynolds equation a more sophisticated constitutive model for viscosity.As it is well known, viscosity of lubricating oils is very sensitive to the operating temperature, showing a quite rapid decrease with increasing temperature.In literature several viscosity-temperature equations are available (for example, the most commonly used are those by Reynolds, Slotte, Walther, Vogel, see [6]).While the most widely used viscosity-temperature chart is ASTM chart based on Walther's equation, the most accurate model is that of Vogel, which can be written as μ=a•exp(b/(T-c)), where a, b, c are oil characteristic parameters.The lubricant viscosity also depends on pressure, especially for high pressures as in heavily loaded concentrated contacts (elastohydrodynamic lubrication).A number of attempts have been made to propose explicit formulae to synthesize lubricant pressure sensitivity.The best known equation for moderate pressures is the Barus equation μ=μ0•exp(αp), in which μ 0 is the viscosity at ambient atmospheric pressure and α [MPa -1 ] is a pressure-viscosity coefficient related to oil film pressure (typical values are α=0.01÷0.02MPa -1 ).In accordance to this constitutive model, an increase in dynamic viscosity occurs for high pressures, with a solid-like behavior for very high pressures.This effect, well-known in elastohydrodynamic studies (e.g.lubricated contacts), has not been investigated in the field of journal bearings.The simultaneous coupled effect of temperature and pressure can then be evaluated by combining two of the above mentioned relationships.For example, a very simple model is the Vogel-Barus equation μ=μ 0 •exp(αp), where now the pressure-independent viscosity term μ 0 is only function of temperature according to the Vogel equation: A further improvement in journal bearing study can be obtained by including in the numerical solution of Reynolds equation the deformed geometry of lubrication gap, caused by elastic deformation of shaft and support under the imposed oil pressure p(θ).This would allow a more realistic estimate of the overall stress distribution on journal bearing components (shaft and housing), compared to other approaches (see for example [7]) that are based on approximate analytical (asymptotic) solutions of Reynolds equation.This paper will develop a general numerical approach to solve Reynolds equation and to compute the resultant oil pressure distribution by including both temperature and pressure effects on viscosity, as well as the effect of components elastic flexibility.A typical journal bearing configuration (see Tab. 1), operating at two different inlet oil temperatures (T in =40 and 70 °C), will be investigated.A viscosity-temperature curve typical of oil ISO VG 680 will be used [6].

NUMERICAL SIMULATIONS
n numerical simulations a plane model for the journal bearing is adopted.In the first part of this paper the hypotheses used are perfectly rigid components and viscosity function of both temperature and pressure according to the Vogel-Barus equation.In the second part of the paper, the pressure effect will be neglected, while shaft and I support elastic deformation will be included into the analysis.

Temperature and pressure effect (with rigid components)
The Reynolds Eq. ( 1) is solved by using the finite difference method based on a central difference scheme.The numerical algorithm is sketched in Fig. 2. The input data are: journal bearing geometry, applied force F, bearing rotational speed N, inlet oil temperature Tin and pressure-viscosity coefficient α, other than the Vogel-Barus model μ(T, p).The unknown function in Eq. ( 1) is the pressure p(θ) that, upon integration, must equal the resultant applied load F. However, the problem at hand is actually non-linear for several reasons.Although the pressure p(θ) is the unknown function, Eq. ( 1) does not explicitly depends on the input load F (that is, the pressure resultant), but on the eccentricity e through the lubrication gap h(θ)=c−e•cos(θ).Therefore, at the beginning of the analysis a guess value of eccentricity e0 (not of the force F) must first be imposed to obtain a tentative lubrication gap and the resultant first-attempt pressure distribution.The Newton-Raphson rule is then applied to compute (based on the resultant of pressure Fi) an eccentricity increment, δei, that is next used to update the eccentricity value for next iteration, ei+1=ei+δei, and to compute an updated lubrication gap geometry h i+1 (θ) for solving again Eq. ( 1).At each iteration, the numerically calculated pressure distribution p(θ) must also be checked for negative values, which have no physical meaning and are set to zero.This adjustment inevitably modifies the overall pressure resultant that balances the applied force, giving rise to another source of nonlinearity in problem solution.Several iterations are required to converge at the correct eccentricity value (a tolerance is checked for both eccentricity and force), which gives the correct pressure distribution p(θ) that solves Eq. ( 1) and also balances (as a resultant of pressure) the applied input force F.

END
Fluid-dynamic analysis (Reynolds eq.)To evaluate the effect of temperature on viscosity (and then on pressure), the journal bearing in Tab. 1 was studied at two operating conditions (JB1, JB2) characterized by two different inlet temperatures (T in =40, 70 °C).Two hypotheses were then adopted to compute the pressure-independent viscosity term μ 0 as a function of oil temperature: in the first, using an average constant temperature Tm resulting by a thermal balance inside the oil film (as in R&B approach) [1], in the second using, as a first approximation, a linear variation from inlet Tin to outlet temperature Tout (that has been calculated by previous thermal balance); note that in both cases the same average oil film temperature T m is obtained.For both temperature distributions within lubrication gap (constant Tm, linear Tin-Tout), the Vogel-Barus equation has been implemented with two different pressure factors (α=0 and α=0.01).Tab. 2 shows an overall comparison of obtained results, while Fig. 3 compares the pressure distribution for different pressure sensitivity values for viscosity (assuming a linear temperature variation within oil film).The effect of temperature variation of oil film is now commented first.Referring to JB1 configuration in Tab. 2, a negligible difference is observed between the case of constant and linearly varying temperature, for both α=0 and α=0.01 values.Instead, larger differences (with a 10-12% increase of pmax value) are observed for JB2 configuration, considering both α=0 and α=0.01 values.This emphasizes how the variation of oil film temperature could have some effect on pressure distribution, at least for high temperature values.Considering the viscosity-temperature strong correlation, this seems to confirm that pressure distribution is more sensitive to a change of small (rather than high) viscosity values within lubrication gap.In any case, the constant temperature assumption used in R&B calculations seems too simplified.Numerical solutions for constant Tm and α=0 were also compared with results given by R&B charts, showing a good agreement only for JB1 configuration, while some difference characterizes JB2 configuration.The observed discrepancy can be attributed to the very low Sommerfeld number (S=0.00298) characterizing JB2 configuration, which makes difficult using R&B design charts and thus can be source of interpolation errors.In general, a non-zero viscosity-to-pressure sensitivity (α=0.01)determines a variation in the overall pressure distribution (change of attitude angle β) and in its maximum value pmax, depending on the general pressure levels attained.Limiting the attention to the case with T in -T out linear temperature variation, for peak pressures p max <100 MPa (case JB1), the pressure effect is actually negligible, as shown in Fig. 3(a)-(b), with only a small decrease of the maximum pressure of about 3.5%.For larger pressure levels (case JB2), an increment of pmax of about 12% is observed, see Fig. 3(c)-(d).The minimum oil thickness increment (h 0 =c-e) produced by the pressure effect is relevant in both cases, with a variation respectively of 28% and 32%.The obtained results can be summarized by saying that, if the influence of pressure on viscosity is taken into consideration, when α increases the peak pressure pmax increases, while the eccentricity e decreases.However, the pressureto-viscosity effect is smaller compared to temperature influence, at least for the maximum pressure peaks encountered in the examples studied.Accordingly, pressure dependence will be intentionally neglected in the rest of this paper.

Effect of component elastic deformation (with α=0, T linear)
In the second part of this work, the pressure distribution will be calculated by considering the real geometry of lubrication gap resulting from component elastic deformation.The pressure values calculated by solving the Reynolds equation ( 1) are applied as input mechanical loads in a FE model to compute the real meatus after deformation, which is next used to solve again equation ( 1) with an iterative analysis scheme.
A fluid-dynamic and structural coupled approach is developed in Matlab environment; the flow chart in Fig. 2 has been integrated with a structural FE analysis toolbox, see Fig. 4. For the structural analysis, the plane FE model and the global stiffness matrices for shaft and support, [Kshaft] and [Ksupp], are first calculated once at the beginning of the analysis.

START
Fluid-dynamic analysis (Reynolds eq.)At first iteration, a guess value of eccentricity e0 is entered into Reynolds equation (1) to compute the pressure distribution p(θ) and the eccentricity e for the case of not deformable components (with the algorithm of Fig. 2).The calculated pressure distribution is next applied as a boundary load in FE model, to calculate components elastic deformations; the relative difference of radial displacements between shaft and support is used to define a nominal gap as g(θ)=ushaft-usupp.A new oil film geometry h'(θ)=c-e•cos(θ)+g(θ) that incorporates mechanical deformation (thus it differs from the case of perfectly rigid components) can thus be calculated.At second iteration step, this updated gap geometry h'(θ) is entered again in Eq. ( 1) to get a new pressure distribution p'(θ).This iterative procedure is repeated until convergence is achieved with respect to an imposed threshold tolerance on the maximum pressure, calculated at each iteration.
The plane FE models of both shaft and housing used in the analysis are shown in Fig. 5.The shaft is modeled by a mapped mesh with 4-nodes isoparametric linear elements, while the housing is free meshed using 3-nodes CST triangular elements.Shaft and support are loaded by the same oil pressure distribution p(θ) applied on the outer and inner surfaces, respectively.The analysis assumes small displacements and a plane strain condition.Material has linear elastic behavior, with properties typical of a structural steel.It is worth noting that the use of a plane FE model for the structural analysis of a journal bearing requires a special attention in modeling mechanical constraints.In fact, in a real journal bearing the applied load F and the resulting pressure distribution are actually applied along different longitudinal locations along the shaft axis.Instead, in the plane FE model here adopted the external load F that balances the oil pressure is replaced by an appropriate constrain on shaft geometry.For this purpose, the shaft has been modeled with a central hole and all nodes on the inner circumference have imposed zero radial displacements; the support, instead, has all the external edges constrained.This modeling strategy, however, affects the shaft structural stiffness: a large inner radius determines an anomalous increment of shaft stiffness, while a very small inner hole gives rise to very large deformations and abnormally high reaction forces at constrained nodes.A proper sensitivity analysis has been preliminary carried out, in order to find the optimal radius of inner hole.As an example, the coupled fluid-structural numerical approach that includes journal bearing elastic deformation was applied for the analysis of JB1 configuration for the case of α=0 and linear temperature variation in the range Tin=40°C−Tout=80°C.Fig. 6(a) shows the calculated pressure distribution, which has to be compared to that of perfectly rigid components shown in Fig. 3(a).The comparison emphasizes that elastic deformation contributes to a reduction of about 48% (from 83 MPa to 43 MPa) of the maximum peak pressure and, accordingly, an increase in the attitude angle β (since the resultant of pressure distribution is always the same applied force F).The pressure profile, more uniform compared to the result for rigid components (R&B solution), seems to support the idea of using the average specific pressure pm as a structural design parameter, as suggested in some design codes [5].
The elastic deformation also induces an increase in eccentricity, from e=0.2392 mm (rigid components) to e=0.2953 (elastic components).Fig. 6(b) compares the geometry of lubrication meatus for the case of deformable and rigid components (angles are referred to the position of minimum oil gap θh0).In particular, it is observed that for deformable components the meatus is not symmetric and that eccentricity can assume values greater than the nominal clearance c, as the elastic deformation can increase the gap between shaft and support.Since the stress distribution of Fig. 7(c) is partly hydrostatic, the maximum von Mises stress (σ vm =21 MPa) is shown to be significantly lower than the maximum absolute radial stress σr=43 MPa.In particular, the elastic deformation of journal bearing elements gives a maximum von Mises stress on the hole surface that is smaller than the maximum oil pressure p max and that is comparable with the static strength of materials usually employed (for instance, white metal generally used as internal coating has a yield stress of about 50 MPa [8]).Instead, in the region underneath the hole surface the overall stress state becomes "less hydrostatic" and a larger von Mises stress (σ vm =35 MPa) then develops.

CONCLUSIONS
n this paper, a numerical procedure for the static analysis of hydrodynamic radial journal bearings has been developed.Influence of temperature and pressure on viscosity and thus on resultant pressure distribution were studied.A mechanical plane finite element model, coupled with solution of Reynolds equation, was also developed to study journal bearing structural behavior and its influence on pressure distribution.The main findings of the work can be summarized as follows:  a temperature increase was shown to give a decrease of attitude angle β and an increase in pressure peak;  an increase of viscosity-to-pressure sensitivity (α value) gives a general increase of peak pressure, especially for pressure peaks greater than about 100 MPa;  the temperature effect was shown to be generally larger than pressure effect on pressure distribution;  the elastic deformation of journal bearing elements gives a more uniform pressure distribution, with a reduced maximum pressure peak compared to the case of perfectly rigid components.In addition, the overall stress state on the surface hole in the support is partly hydrostatic, so that the maximum von Mises stress is lower than the maximum radial stress modulus (i.e. the maximum peak pressure).I

Figure 1 :
Figure 1: Sketch of a hydrodynamic journal bearing and parameters used in numerical simulations

Figure 2 :
Figure 2: Sketch of the numerical algorithm for fluid-dynamic analysis of journal bearing (rigid components)

Figure 4 :
Figure 4: Sketch of the numerical algorithm for fluid-dynamic and structural numerical analysis of journal bearing.

Figure 5 :
Figure 5: Finite element model of (a) shaft and (b) housing

Figure 7 :
Figure 7: (a) Radial and (b) von Mises stresses in the support (MPa units); (c) stress components on inner surface of support.For what concerns the calculated mechanical stresses, Fig.7(a)-(b) show the overall distribution of radial and von Mises stresses within the housing, while Fig.7(c) plots the values of stress components (radial σ r , hoop σ θ , axial σ z ) and von Mises stress on the hole surface of the housing, as a function of angle θ (the vertical symmetry axis is at angle θ=0).All stress values are compressive and show a similar trend with θ.The maximum absolute radial stress is 43 MPa, exactly equal to the maximum oil pressure (p max =43 MPa) applied on the inner hole.Hoop and axial stress components have approximately similar values, with the axial stress under plane strain condition calculated as σz=ν(σθ+σr), where ν is the Poisson ratio.Since the stress distribution of Fig.7(c) is partly hydrostatic, the maximum von Mises stress (σ vm =21 MPa) is shown to be significantly lower than the maximum absolute radial stress σr=43 MPa.In particular, the elastic deformation of journal bearing elements gives a maximum von Mises stress on the hole surface that is smaller than the maximum oil pressure p max and that is comparable with the static strength of materials usually employed (for instance, white metal generally used as internal coating has a yield stress of about 50 MPa[8]).Instead, in the region underneath the hole surface the overall stress state becomes "less hydrostatic" and a larger von Mises stress (σ vm =35 MPa) then develops.

Table 1 :
Parameters used in numerical simulations

Table 2 :
Overall comparison of results for numerical simulations with rigid components.