Nonlinear Dynamic Response of an Unbalanced Flexible Rotor Supported by Elastic Bearings Lubricated with Piezo-viscous Polar Fluids

On the basis of the V. K. Stokes micro-continuum theory, the effects of couple stresses on the nonlinear dynamic response of the unbalanced Jeffcott's flexible rotor supported by layered hydrodynamic journal bearings is presented in this paper. A nonlinear transient modified Reynolds' equation is derived and discretized by the finite element method to obtain the fluid-film pressure field as well as the film thickness by means of the implicit Euler method. The nonlinear orbits of the rotor center are determined by solving the nonlinear differential equations of motion with the explicit Euler's scheme taking into account the flexibility of rotor. According to the obtained results, the combined effects of couple stresses due to the presence of polymer additives in lubricant and the pressure dependent viscosity on the nonlinear dynamic response of the rotor-bearing system are significant and cannot be ignored or overlooked. As expected, these effects are more noticeable for polymers characterized by higher length molecular chains.


Introduction
Rotating machines such as turbo-machines are the most important class of machinery that are extensively used in diverse engineering applications such as power stations, aircraft propulsion systems, etc.These machines include at least one rotor generally supported by oil-lubricated journal bearings which must not be considered to be passive elements, but elements which sensitively affect the dynamic behavior and stability of the rotor.
The design trend of such mechanical systems in modern engineering is towards lower weight and operating at super critical speeds.
The role of rotor is to transmit or to transform the mechanical power.It is often of very much complex realization and includes bladed disks or impellers, gears, for example.
The correct prediction of dynamic behavior is extremely essential when the rotor rotating at high speeds is made very flexible.Further, the dynamic behavior of a rotor-bearing system largely depends on the nonlinear dynamic characteristics of oil-lubricated bearings which may be a source of self-induced vibrations popularly known as oil whirl and oil whip phenomena which occur especially when the rotor-bearing system is lightly loaded or operates at low values of eccentricity ratio.This vibratory motion can cause considerable mechanical problems, like rubbing between shaft and bearing, blades and stator in turbo-machines, or more generally vibrations of the whole rotating machinery.In 1924, Newkirk and Taylor [1] first demonstrated that the oil whirl is caused by the dynamic oil film forces in the bearings.Since then, a number of researchers have focused their theoretical and experimental studies on the stability, bifurcation, and chaos of rotor-bearing systems, and a significant progress in rotor-dynamics field has been made [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21].
Of the many published works, the most extensive portion of the literature on rotor dynamics is concerned with determining critical speeds, natural whirl frequencies, instability thresholds, and imbalance response.In these works generally based on the linearized theory, the authors used in their analysis several assumptions among them the rotor and bearings are assumed to be rigid.Moreover, the linear stability analysis predicts the stability limits only under small disturbances, it does not give any information about the transient phenomena of a rotor-bearing system under large disturbances.Those phenomena should be studied by means of a nonlinear stability theory.
On the other hand, the rheological behavior of mineral oils used as lubricants in industrial machinery is influenced by the presence of additives such as Viscosity Index (VI) improver polymers which are characterized by long-chains.Thus, Oils containing VI additives such as multi-grade engine oils must be considered as non-Newtonian fluids and can be modeled as polar fluids.Experimentally, it was found that the presence of dissolved polymer in the lubricant increases the load carrying capacity of the lubricating film and decreases the friction coefficient [22,23].
In order to better describe the rheological behavior of this kind of non-Newtonian lubricant so-called couple stress fluid, different micro-continuum theories have been developed [24,25].The Stokes micro-continuum theory [26] is the simplest theory of fluids proposed in the technical literature, which allows the polar effects such as the presence of couple-stresses mi and body couples noted ℓi in addition to the body forces and surface forces Figure 1.The iso-volume couple-stress fluids are characterized by two constants, namely  and  whereas only one parameter appears for a Newtonian iso-volume fluid which is the dynamic viscosity  .The new material constant  is responsible for couple-stress property.In the literature, the effects of couple-stresses on the behavior of journal bearings are studied by defining the dimensionless couple-stress parameter which has the dimension of length and can be thought of as a fluid property depending on the size of the high polymer molecule.
Owing to its relative mathematical simplicity, the couple-stress fluid model has been widely applied to analyze various hydrodynamic lubrication problems: hydrostatic thrust bearings, slider bearings, rolling bearings, squeeze film bearings, layered journal bearings, and dynamically loaded engine journal bearings [27][28][29][30][31][32][33].The theoretical results obtained showed that the presence of the couple-stress provides an enhancement in the load carrying capacity, improve the dynamic performance characteristics and the stability of journal bearings, and lengthens the response time of the squeeze film action of the system as compared to the Newtonian lubricant case.It was also shown that the effects of couple stress are more pronounced for high values of the couple-stress parameter ℓ.
According to the Stokes micro-continuum theory of couple stress fluids [26] together with the variation of viscosity with pressure, the combined effects of piezo-viscous dependency and non-Newtonian couple stresses on wide parallel rectangular-plate squeeze-film characteristics have recently been studied by Lin et al. [34].The viscosity-pressure relationship used by the authors is that of Barus [35] which is largely used in elastohydrodynamic lubrication analysis.They have found that the combined effects of piezo-viscosity and couple stresses cause an increase in the load carrying capacity and the squeeze film time.
As far as we know, very few research works have been devoted to analyze the nonlinear dynamic behavior of flexible rotors mounted in layered journal bearings lubricated with complex non-Newtonian fluids taking into account the combined effects of fluid couple-stresses, pressure-viscosity dependency, and compliance of the bearing structure.
The main objective of the present research is to theoretically investigate the nonlinear dynamic response of a flexible rotor supported by layered bearings using piezo-viscous couple-stress fluids as lubricants.The elastic deformation of the liner due to the fluid film pressure is calculated using a simplified elasticity model known as the "thin liner model".So, we assume that the radial deformation at the fluid-liner interface is proportional to the hydrodynamic pressure.On the other hand, the Barus law is used in order to take into account the viscosity-pressure dependency at constant temperature: (1) where is the atmospheric dynamic viscosity, and α is the pressure-viscosity coefficient which can be obtained by plotting the natural logarithm of dynamic viscosity versus pressure p.The slope of the graph corresponds to the value of α.The pressure-viscosity coefficient is a function of the molecular structure of the lubricant and its physical characteristics.If α = 0, the viscosity is then constant and the fluid is considered as iso-viscous.
There are various formulae available to calculate the pressure-viscosity coefficient.One of the early ones was derived by Wooster [37]: where α is the pressure-viscosity coefficient in (Pa −1 ), and is the atmospheric dynamic viscosity in (cP) or (mPa•s).There are also other equations for the calculation of the pressure-viscosity coefficient available in the literature.Some of these equations are accurate for certain fluids and inaccurate for others.One of the problems associated with available formulae is that they only allow the accurate calculation of pressure-viscosity coefficients at low shear rates.An accurate value of this coefficient can be determined experimentally.
The piezo-viscosity effect varies between oils; and it is more considerable for naphthenic oils than paraffinic oils.Water; by contrast shows only a small rise; almost negligible; in viscosity with pressure.
Note that the value of the pressure-viscosity coefficient is in general reduced at higher temperatures, with naphthenic oils being the most severely affected.In some cases even at moderate temperatures, there is a substantial reduction in the pressure-viscosity coefficient.
There are many other formulae for viscosity-pressure relationships.A short review of some of the empirical formulae for the viscosity-pressure relationships is given in [38].These formulae allow for the calculations of viscosity changes with pressure under various conditions and to various degrees of accuracy.
As a first approximation, we consider the simplest flexible rotor model popularly known as the Jeffcott model (1919).In spite of its age of over 90 years, the Jeffcott rotor is still widely used.
Figure 2 shows a typical Jeffcott rotor.It consists of a simply supported flexible massless shaft with a rigid disc mounted at the mid-span.The disc center of rotation, C, and its center of gravity, G, is offset by a distance, e which is called the unbalance eccentricity.The shaft spin speed is ω and the shaft whirls about the bearing axis with a whirl frequency ν.For present study, synchronous condition has been assumed, i.e., ν = ω.According to the beam theory, the stiffness of the shaft kr can be expressed as load deflection

48
(3 where Er is the Young's modulus of rotor, is the second moment of area of the rotor (shaft) cross-section, and Lr denotes the span of the rotor (shaft).The finite element method is used to approach the nonlinear transient pressure equation called in the present study the modified Reynolds' equation derived from momentum and mass conserving principles for an incompressible and piezo-viscous couple stress fluid using the V. K. Stokes micro-continuum theory [26].The first order differential equations system resulting from the discretization of the pressure equation is solved by the implicit Euler's scheme and the relaxed substitution iterative method in order to determine the instantaneous hydrodynamic forces acting on the rotor (shaft), i.e., at each time step.The nonlinear trajectories of the shaft center are obtained by solving the rotor-dynamics equations with a direct integration procedure, namely the explicit Euler's scheme.
In the parametric study, there are three key parameters dominating the combined effects of the pressure-viscosity dependency, the presence of couple-stresses, and the compliance of the bearing-liner which are the dimensionless viscosity-pressure coefficient , the dimensionless couple stress parameter ℓ , and the normalized elasticity parameter ∁, respectively.

Momentum Equations of the Polar or Couple Stress Fluid
On the basis of the V. K. Stokes micro-continuum theory (26), when the body forces i b and body couples i  are neglected Figure 1, the motion of the iso-volume polar fluids can be governed by the following equations: which is the material time derivative, and The constitutive equations for the polar or couple-stresses fluids are given as 1 , , 3 4 4 where 1 ω 2 is the ith component of the vorticity vector.
If we now consider the case where  and  are constants, the conservation of momentum is given by [26]: where  is the fluid density, i u is the velocity component, t is the time, p is the pressure,  is the classical shear viscosity and  is a new material parameter responsible for the couple-stress property.
Note that the dimension of the material constant  is that of viscosity (ML −1 T −1 ) whereas the dimensions of  is that of momentum (MLT −1 ).The ratio   has dimension of length square and, in the following, we denote this material constant by  , where . Some experiments for determining the material constants  and  for incompressible fluids are given in Equation (1).
Physically, the quantity  can be regarded as a characteristic length of the additives which are added to the base oil and which can be polymers or co-polymers.
In thin film theory, the dimension across the film thickness is small compared to the others.From this assumption, we can determine the order of magnitude for the different terms of Equation (9).
The field equations governing the motion of the lubricating oil in the Cartesian coordinates system are: where V  is the flow velocity field.According to the Barus law, the second term appearing on the left hand side of Equation ( 10) can be expressed as follows: according to the second equation of Equation (10).
The boundary conditions at the bearing surfaces (y = 0) and the journal surface (y = h) are: Conditions Equations ( 12) and ( 14) are the no-slip velocity conditions and conditions Equations ( 13), (15) result from the fact that the couple-stresses ( ) vanish at the solid boundary.
Using the above boundary conditions, the fluid velocity components u and w for an iso-volume piezo-viscous couple stress fluid are: h being the fluid film thickness depending on x, z and t variables.

Modified Transient Nonlinear Piezo-Viscous Reynolds' Equation and Boundary Conditions
Substituting the expressions of the velocity components u and w from Equation ( 16) in the continuity Equation (11) and integrating across the film thickness, i.e., with respect to y variable, using the boundary conditions for v given below: yields to the following transient nonlinear modified Reynolds' equation: where   , and being the volume flow rate components per unit length in x and z directions, respectively.According to Equations ( 17) and ( 19), Equation ( 18) becomes where Uj = ɷR is the linear velocity of the journal surface.
Equation (20a) is known as the modified Reynolds' equation which is a nonlinear partial differential equation.
For an iso-viscous and nonpolar (Newtonian) fluid for which 0     , Equation (20a) reduces to the standard linear Reynolds' equation since  Using the non-dimensional quantities (indicated with a superscript asterisk (*)) In Equation ( 20), * , * , / are the Cartesian instantaneous co-ordinates of the journal centre and the parameter ∁ represents the compliance factor of the elastic bearing-liner which is defined as where , The normalized pressure field must satisfy the modified Reynolds' Equation ( 21) in * 0,2 1 2 , 1 2 and the following boundary conditions on * : In addition, the Gümbel's condition also known as the half-Sommerfeld condition [39], suggested in the early 1920s, is used to take into account the film rupture occurring in the divergent region of the bearing.This condition neglects the sub-ambient pressure completely when calculating bearing performance and it is frequently used in the numerical simulations of hydrodynamic lubrication problems owing to its simplicity.
Although this condition does not respect the mass flow continuity between the over-ambient pressure (leading edge) region and the sub-ambient pressure (trailing edge) region, it does give a theoretical pressure curve comparable to that obtained experimentally or by using the mass conserving algorithms especially for highly loaded short journal bearings (L/D < 1) with low supply pressures.
Dowson and Taylor [40] presented a good review of cavitation phenomenon in hydrodynamic journal bearings, where different cavitation models were largely discussed.
For an aligned bearing, the pressure field is symmetric about the bearing mid-plane section.For computational efficiency, we calculate only half of the pressure field.

Rotor-Dynamic Equations
When the external load acting on the rotor or the journal varies both in direction and magnitude, the journal center describes a trajectory within the bearing.The determination of this trajectory requires the solution of the rotor-dynamic equations coupled with those governing the hydrodynamic behavior of lubricating oil films.
Generally, the external loads acting on the rotor are: the weight of the rotor 2W0 = 2 mg; -the dynamic load components due an unbalance mass characterized by its eccentricity e; -the hydrodynamic forces FX and FY due to the presence of the lubricating oil film.
The rotor-bearings system is modeled as a flexible Jeffcott rotor symmetrically supported by two identical layered hydrodynamics bearings ignoring the gyroscopic effects as depicted in Figure 2. We assume that the rotor mass is lumped at the midpoint with massless shaft, central load of 2m, rotor damping of 2br, and rotor stiffness of 2kr.
For each bearing are attributed a mass m of the rotor, a static applied load W0 = mg, and a synchronous dynamic excitation due to an unbalance mass 2 ) ( The application of the dynamic fundamental principle gives where FX and FY are the hydrodynamic forces in X and Y directions that are nonlinear functions of the displacement components (X,Y) and the velocities , of the journal center Oj.They are calculated by integrating the hydrodynamic pressure over the bearing surface.This latter is obtained by solving the transient modified Reynolds' equation by finite element method.
In dimensionless variables, Equation ( 25) is written as where Solving the normalized equations of motion ( 27) ahead of time numerically using the direct integration methods of second order differential equations such as the explicit Euler's method, the transient motion of the journal center Oj for balanced and unbalanced shafts is determined Figure 4 (See Appendix A for more details).

Weak Integral Formulation and Finite Element Discretization
In steady-state conditions, the normalized modified Reynolds' equation becomes The Galerkin's weighted residual method is employed to build up the integral formulation associated to Equation (28a), i.e.,  After writing Equation (30) over a subdomain * and discretisation by the finite element method, we obtain a system of nonlinear algebraic equations Over the parent element, the bilinear shape functions Ni (Lagrangian polynomials) are expressed as where  and  are the local coordinates of an element such that Accordingly, static pressure and static film thickness over a parent element are then approached by

The coefficients
i ij e e f k and are numerically evaluated over the parent element using the Gauss-Legendre quadrature with ( 2 2 ) integration points as depicted in Figure 5.After assembling of different elementary matrices, finite element approximation leads to the following nonlinear system: where     since the element is isoparametric, i.e., the element geometrical nodes are also the interpolation nodes.The lift force components Equation ( 35) are also calculated by means of the Gauss-Legendre quadrature with ( 2 2 ) integration points.

Method of Solution of Steady-State Nonlinear Modified REYNOLDS' Equation
The steady-state solution of the modified Reynolds' Equation (28a) which is a highly nonlinear elliptic partial differential equation is obtained by the relaxed substitution iterative method.This method consists of building up a series of solutions  k P by solving the linear system: We can write this in incremental form by introducing the residual vector . : Geometric nodes : Interpolation nodes where Ω0 is an under-relaxation factor which ensures and accelerates the convergence of the iterative process.
To obtain the steady-state fluid pressure, the required steps of the computational procedure are presented in detail in Appendix A.

Iterative Research of the Steady-State Equilibrium Position
The steady-state equilibrium position of shaft center (X0,Y0) due to the application of static load is determined when the integration of the resulting steady-state hydrodynamic pressure obtained from Equation (28a) balances with the applied load components   In this case, the bearing load is known as a constant and the position of the shaft center is to be calculated from numerical iterations.Therefore, an inverse solution of the steady-state modified Reynolds' equation is required (See Appendix A for details).

Weak Integral Formulation and Finite Element Discretization
Applying the Galerkin weighted residual method and the Green's theorem, the weak integral form associated to the transient modified Reynolds' Equation ( 21 According to Equation (40), Equation (39) becomes: (43)

Method of Solution of Transient Nonlinear Modified Reynolds' Equation
We found in the preceding section that solving the modified Reynolds' equation under dynamic conditions by the finite element reduces ultimately to the solution of simultaneous nonlinear first order ordinary differential equations of the form Equation (43).
After formation of matrix Equation (43) the matrices must be restructured or modified to account for any essential boundary conditions Equation (24).In addition to the boundary conditions, we must also know the initial conditions, i.e., 0 .There are many general methods and several special techniques for solving first-order matrix differential equations.We consider the direct numerical integration in time to solve the original set of coupled equations.The procedure relies on deriving recursion formulas that relate the values of at one instant of dimensionless time ̃ to the values of at a later time ̃ ∆ , where ∆ ̃ is the dimensionless time step.The recursion formulas make it possible for the solution to be marched out in time, starting from the initial conditions at time ̃ 0 and continuing step by step until reaching the desired duration.where the * ∆ * on the left-hand side of the equation are unknowns, and all of the terms on the right-hand side are known.Equation (46) represents a general family of recurrence relations.If β = 0, we get the explicit Euler algorithm, i.e., and if β = 1 we obtain the implicit Euler method If 0 0 , we get the nonlinear algorithm for the steady-state modified Reynolds' Equation.

Transient Solution vs. Steady-State Solution
The nonlinear dynamic analysis, described in the Section 5, has been incorporated into a computer program.
In order to check the correctness of the results obtained from the nonlinear transient analysis of hydrodynamic journal bearings, a separate 2-D finite difference computer program was developed to calculate the steady-state position by solving the inverse hydrodynamic lubrication by means of the relaxed iterative Newton-Raphson method as stated in Section 4. For (2-D) finite difference analysis, which in this paper will be referred to as the exact solution, a grid size of (61 × 16) nodes is used in θ and z directions, respectively.The number of nodes is chosen in order to ensure accurate results with minimum CPU time.The rotor and bearing parameter values used in the calculations as well as the finite element meshing characteristics and data of calculations are reported in Tables 1 and 2. In this section, the calculations were made for two different cases: Case 1: A balanced rigid journal bearing lubricated with Newtonian and iso-viscous lubricant: W0 = 340 kN (S = 0.18, * = 1.8), α * = 0, ℓ * = 0, Cd = 0, ε = 0.
For both cases, the geometrical characteristics such as the bearing length, the bearing radius and the radial clearance as well as the rotational velocity are the same and are given in Table 1.  Figure 6 shows the shaft center trajectory obtained using the transient analysis in the case 1, i.e., when the bearing is rigid and the lubricant is Newtonian.We can notice in this case that this trajectory tends to reach the shaft equilibrium position obtained from the steady state analysis which locates at ( * = 0.431287, * = 0.499599).These values are almost close to those predicted by the transient analysis which are (X * = 0.439241, Y * = 0.501445).In Figure 7, we represent a comparison of film thickness and pressure profile obtained by the transient and steady state methods in the case 1.A very good agreement is observed.Figure 8 depicts the shaft center trajectory calculated with the transient analysis in the case 2, i.e., when considering the compliance of the bearing-liner and the lubricant as a non-Newtonian couple stress fluid.Similar trends for case 2 are obtained as in case 1, i.e., the path of shaft center gradually approaches the static equilibrium position whose the co-ordinates are ( * = 0.3143, * = 0.3268).The values predicted by the transient analysis are (X * = 0.3272, Y * = 0.3136).As in the case 1, the discrepancy between the two results is very small.Figures 9 and 10 illustrate the differences between the results for a rigid bearing lubricated with an iso-viscous Newtonian fluid and for a compliant bearing using a piezo-viscous polar fluid as lubricant.

Flexible Rotor with Large Mass
Figure 9 represents the trajectories of the shaft center in given running conditions: W0 = 340 kN, th = 10 −2 m, E = 0.9 GPa, σ = 0.35, br = 0, kr = 5 MN/m, N = 3000 rpm or ω = 100 × π rad/s, and ε = 0.80.Note that for the rigid bearing-liner, the Young's modulus of liner E tends to infinity and the deformation coefficient is then equal to zero because it is inversely proportional to the elasticity modulus.
The relative unbalance eccentricity ε = 0.80 corresponding to e = 280 microns generates a dynamic load W = meω 2 = 940 kN at synchronous frequency (ν/ω = 1).For this operating condition, the value of W is greater than the static load W0 (W/W0 ≈ 3) and is representative of some emergency conditions in turbomachinery when a blade is lost for example.
It is noted that only the final form of shaft center orbits will be presented in the following, i.e., the results corresponding to the transient numerical effect due to initial conditions will be omitted.
The trajectories of the shaft center within the compliant bearing * 1.908, 0.35 * 0.04 are given in Figure 9a for two values of couple stress parameter l * = 0.0 (Newtonian case) and l * = 0.0 (non-Newtonian couple stress fluid), and different values of the piezo-viscosity coefficient ranging from 0.0 (iso-viscous case) to 0.50.In these operating conditions, the rotor has a very large amplitude of circular or pseudo-circular motion and the nonlinear dynamic behavior appears clearly.This is due to the fact that the dynamic loading due to a large unbalance mass is very important compared to the static one as afore-mentioned above.
On the other hand, we observe in the Newtonian case l * = 0 that increasing the pressure coefficient α * shorten the shaft trajectories.This is due to the pressure rise leading to a higher load carrying capacity which reacts and reduces the trajectories size.We can notice that the operating eccentricity of the journal bearing can be greater than the radial clearance in this case where the unbalance mass is large.Furthermore, the orbits described by the shaft center are widely modified by the bearing compliance especially for the non-Newtonian case as clearly illustrated in Figure 9a.They exactly follow the shape of the deformed bearing.
Figure 9b represents the case of rigid bearing-liner.It is seen that the influence of pressure viscosity α * on the orbits is weak in both Newtonian and non-Newtonian cases.That means in this case the distortion plays a major role in the bearing response.Moreover, the non-Newtonian orbits are smaller than those obtained in the Newtonian case, i.e., when the couple stress parameter l * of the lubricating fluid increases.

Flexible Rotor with Small Unbalance Mass
The relative unbalance eccentricity ε = 0.20 corresponding to e = 70 microns generates a dynamic load W = meω 2 = 235 kN at synchronous frequency (ν/ω= 1) which corresponds to 0.70 times the static load.This defect may be attributed to a large residual unbalance which could exist in the shaft.
The trajectories of the shaft center within the rigid and compliant bearings are given in Figure 10 for the same values of the couple stress parameter l * and the piezo-viscosity coefficient.As expected, the shaft center moves around the equilibrium position in both cases since the dynamic load is smaller than the static one.Moreover, the co-ordinates of the equilibrium position change when increasing the piezo-viscosity coefficient especially for higher values of this coefficient in both Newtonian and non-Newtonian cases.So, it results a shift of orbits towards the bearing center as depicted in the figure .For small unbalance value as the pressure magnitude is lower the impact of bush distortion is less pronounced than for a large unbalanced case.
For the rigid case, we observe in Figure 10b the same tendencies with smaller size trajectories.

Conclusions
In this paper, a transient nonlinear analysis was developed and presented with details to investigate the nonlinear dynamic response of an unbalanced Jeffcott flexible rotor mounted in layered journal bearings lubricated with a piezo-viscous polar fluid.The analysis is based on the V. K. Stokes micro-continuum theory for describing the flow of piezo-viscous polar lubricants blended with additives.Using the classical assumptions of hydrodynamic lubrication, a transient and nonlinear modified Reynolds equation was derived in order to take into consideration the combined effects of pressure dependent viscosity and couple stresses resulting from the presence of polymer additives in the base oil.The trajectories of the shaft center were obtained numerically by solving the rotordynamics equations with the explicit Euler's scheme.
The two dominant parameters in the present analysis are the viscosity-pressure coefficient and the couple stress parameter.Such parameters must be considered as key parameters to design adequately bearings of rotating machineries.
The obtained results have been compared with the iso-viscous and Newtonian case.
According to the theoretical results presented, the combined effects of couple stresses due to the presence of polymer additives in lubricant and the pressure dependent viscosity provide lift load enhancement and produce higher oil-film thickness and more contracted trajectories of the shaft center even at severe running conditions (e.g., higher loading, large unbalance mass, etc.).In these circumstances, the destructive metal-to-metal contact between the shaft and bearing surfaces may also be avoided because of the elasticity of the bearing.Qualitatively, these results agree very well with those obtained by the same authors [33] in the case of internal combustion engine connecting-rod bearings lubricated with an iso-viscous polar fluid.As shown in Figure 4, the steps of the computation procedure may be summarized as follows:

Nomenclature
Step 1: At time * 0 , for an initial position of the shaft center * * and initial velocities ′ * * 0 0 , we solve the set of nonlinear algebraic equations resulting from discretization of the normalized steady-state modified Reynolds Equation ( 21) (i.e., without the transient term) by the finite element method using the relaxed substitution iterative method in order to obtain the hydrodynamic pressure field * and the film thickness distribution * , and by integration of this latter over the bearing surface, we get the oil film force components: Step 2: The acceleration components of the shaft center ′′ * * are then determined from dynamic Equation (27).The choice of non-dimensional time increment, ∆ * , is important in this kind of computation procedure, i.e., a small time step would imply longer computation times whereas a larger value may involve calculation instabilities.The time increment value depends on the nature of the problem under consideration.For an unbalanced rotor response, it is generally ranged in the interval , , and it is difficult to guess it.

Method of Solution for Steady-State Nonlinear Modified Reynolds' Equation
Step 1: Select the input parameters of the problem , under-relaxation factor 0  whose the value ranges from 0 to 1, convergence criterion p  and maximum number of iterations, max k for the steady-state pressure solution.
Step 2: Initialize the iteration number k to 0, the norm n to 1, and the global vector containing nodal dimensionless steady-state pressures     0 0  k P .
Step 3: Calculate the dimensionless static film thickness profile using Equation (31b) for each node of the finite element grid.
and the relative least square norm of

End do while
Step 4: Apply the Gümbel's rupture film conditions by vanishing all the negative terms of calculated pressure.
Step 5: Calculate the steady-state lift force components using Equation (35).

Iterative Research of the Steady-State Equilibrium Position
The relaxed iterative Newton-Raphson method is used to solve Equations (38a) and (38b).Therefore, these equations can be rewritten in residual form.
In the Newton-Raphson method, the (k + 1)th trial solution is where Ω is the relaxation factor ranged in the interval (0, 1).
When the left-hand side of above equation equals zero and the higher terms are truncated, the shaft The analytical inversion of the Jacobean matrix gives: The steps of the Newton-Raphson procedure may be summarized as follows: Step 1. Select the input parameters of the procedure which are: The convergence criterion of the iterative procedure 1  , the maximum number of iterations max Step 3. Calculate the steady-state hydrodynamic lift components: where  * = 10 −6 . Step

Figure 1 .
Figure 1.Balance of forces and couples acting on elementary volume according to V. K. Stokes theory.

Figure 3 .
Figure 3. Photograph and schematic representation of a coated journal bearing.

Figure 4 .
Figure 4. Flow chart of the computational procedure.
Green's theorem (integration by parts), the weak form of Equation (29) is obtained further.As a result, the weak form contains lower order partial derivatives of * 0 p and includes two types of integration.Namely, domain integration over computational domain and contour integration along the boundary of * D , i.e., * D  , the latter can be neglected because * 0 p  is zero at both ends of the housing

)
Isoparametric four-noded rectangular elements with sides parallel to the global axes are used in the discretization of the fluid film domain * D such that

Figure 5 .
Figure 5. Representation of integration points over the iso-parametric parent element.
containing the term appearing on the RHS of the stationary modified Reynolds' equation.1  e A being the assembling operator.The static pressure field is integrated over the half bearing domain to generate the dimensionless lift force components, i.e., ) is * * , ℓ * , * ,

Ch 3 mz 1 *
bearing radial clearance, m E Young's modulus of the bearing-liner, Pa fluid-film thickness, m * h dimensionless fluid-film thickness C h  , 0 h static fluid-film thickness, m * 0 h dimensionless static fluid-film thickness C h dimensionless couple-stress parameter= C  , ∁ scalar compliance operator, N m mass of rotor per bearing, kg N rotation velocity of the rotor, rpm p fluid-film pressure, Pa * on the journal bearing X,Y displacement components of the journal centre, axial co-ordinate measured from middle section plane of the bearing, m * z non-dimensional axial co-ordinate L z  ,  pressure-viscosity coefficient, Pa −non-dimensionless pressure-viscosity coefficient material constant responsible for couple-stresses, kg•m•s −1  absolute viscosity of lubricant, Pa•s 0 absolute viscosity at atmospheric pressure, Pa•s  Poisson's ratio of the bearing-liner  bearing angle, rad  angular velocity of the rotor (shaft)= 60 2 N  , rad/s    square matrix  line vector, =  T Computation Procedure of the Journal Bearing Nonlinear Dynamic Response

Step 3 :Step 4 :Step 5 :Step 6 :
The new position together with the velocity components of the journal center are predicted at dimensionless time ̃ ∆ ̃ from the explicit Euler's scheme, that is: The new position and velocities allow us to calculate new pressure and film thickness distributions * ∆ * * and * ∆ * * by solving the first order differential equations system resulting from the discretization of the transient modified Reynolds' equation by means of the implicit Euler's scheme (β = 1) and the substitution iterative method.Thus, new components of the hydrodynamic (lift) force * * ∆ * * * ∆ * can be calculated through integration of * ∆ * * .The acceleration components at time * ∆ * are then calculated from dynamic Equations (27) considering the new values of the lift force components * * ∆ * and * * ∆ * .The next time interval is then considered while * * , the values obtained earlier play the role of * * and * * in Step 3.

3 . 3 .F 3 . 5 .F by introducing the essential boundary conditions 3 . 6 .
Initialize global matrices [K] and {F} to 0. 3.4.For each element: 3.4.1.Extract the elementary vector  the elementary global co-ordinates arrays of each node by means of the connectivity array 3.4.2.Compute the elementary matrices Form the restructured matrices   r K and   r Solve the restructured linear system     by means of the first-order two-variable Taylor series expansions of the functions given by

Initialize the iteration number k to 0. Step 2 .
facteur Ω, and the initial approximation of the steady-state position of shaft center Solve either the steady-state HD lubrication problem governed by Equation (31a) for * 0 p , or the steady-state EHD lubrication problem described by the coupled Equations (31a) and (31b) by iterations for

Step 5 .
Evaluate Jacobean matrix coefficients by numerical differentiation:

6 .Step 7 .
Calculate the corrections   Calculate the new approximations of the nonlinear system: is reached, i.e., the values of return to Step 2 for another iteration.

Table 1 .
Rotor and bearing parameter values used in the calculations.

Table 2 .
Finite element meshing characteristics and data of calculations.