Effect of wear on the dynamic characteristics of a rigid rotor supported by journal bearings

The purpose of this work is to examine the effects of wear severity on the static and dynamic characteristics of journal bearings, and on the vibration response of a rigid rotor supported by journal bearings. Numerical simulations using MATLAB was conducted for three different operating regimes, namely low loaded operating regime (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _0 = 0.15$$\end{document}ϵ0=0.15), moderately loaded operating regime (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _0 = 0.45$$\end{document}ϵ0=0.45) and highly loaded operating regime (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _0 = 0.75$$\end{document}ϵ0=0.75) with wear depth parameter ratio (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta$$\end{document}δ) varied from 0 to 0.5 at increments of 0.1. Numerical results showed that the vibration response of the rotor generally increases with the increase of the wear depth for all cases of low, moderately and highly loaded operating regimes of the bearings. For the values of parameters considered in this work, it was shown that the vibration response amplitude of the rotor in worn journal bearings may be six times larger compared to the response amplitude of the rotor in non-worn bearings.


Wear model
The wear model used in this paper is the model proposed by Dufrane et al. 13 and it is shown schematically in Fig. 1 5 .The region of positive pressure in the state of worn journal bearing can be divided into three sub regions 20 .The derivation of the equation presented herein follows the work of Papanikolaou et al. 5 and Jamil et al. 20 .
where, h 0 is the fluid film thickness in the non-worn region (m).C is the radial clearance (m).ǫ 0 is the eccentricity between the journal center and bearing center at equilibrium position (m).θ is the angular coordinate measured from the position of maximum fluid-film thickness in the direction of rotor angular velocity (counter-clockwise), (rad).
For small amplitude motions about the journal equilibrium position, the fluid-film thickness, h, can be rewritten as shown in the following equation: where, h is the fluid-film thickness for small amplitude motions about the journal equilibrium position (m).θ = γ − ϕ 0 ; γ is the angle measured from the vertical axis, X, in the direction of rotor angular velocity (counter-clockwise), (rad).
ϕ 0 is the attitude angle measured from the vertical axis, X, to the position of maximum fluid-film thickness in the direction of rotor angular velocity (counter-clockwise), (rad).
The time-dependent eccentricity ǫ(t) and attitude angle ϕ(t) are given by: (1) e and ϕ are small radial displacement and small variation of attitude angle respectively.Substituting Eq. (3) into Eq.( 2) yields: For small amplitude of motion, the following relationship is valid.Substituting Eq. (5) into Eq.(4) gives the following expression for the instantaneous fluid-film thickness in the non-worn region.
Differentiating Eq. ( 6) with respect to θ results in Eq. (7), which describes the change of fluid-film thickness in the circumferential direction of the bearing.Differentiating Eq. ( 6) with respect to time, t, results in Eq. (8), which describes the change, with time, of the fluid-film thickness at a fixed point in the bearing.
where, h 0w is the worn fluid film thickness at equilibrium position (m); C is the journal bearing radial clearance (m).d 0 is the wear depth (m).Substituting Eq. (3) into Eq.( 9) yields: www.nature.com/scientificreports/Substituting Eq. ( 5) into Eq.(10) gives the following expression for the instantaneous fluid-film thickness in worn region.
Differentiating Eq. ( 11) with respect to θ yields Eq. ( 12), which describes the change of fluid-film thickness in the circumferential direction of the bearing.
Differentiating Eq. ( 11) with respect to time, t, yields Eq. ( 13), which describes the change, with time, of the fluid-film thickness at a fixed point in the bearing.
In this region the thickness of lubricant oil fluid is as in the 1 st non-worn region.The starting and final angles of the region where the wear takes place (refer to Fig. 1) are given by: where, δ = d 0 c ; δ is the wear depth parameter.θ s and θ f are represented by the following equations 20 :

Solution to bearing forces and dynamic coefficients
The fluid film pressure is obtained by solving the Reynolds equation for infinitely short bearing (L/D < 0.5).The bearing forces are computed by integrating the fluid film pressure over the bearing area using Simpson's 1/3 rule.

Reynolds equation
For convenience purposes, polar coordinates is employed instead of Cartesian coordinates to solve the Reynolds equation.The polar coordinates used are the eccentricity ratio, ǫ 0 and the attitude angle ϕ 0 .
With the following assumptions: (i) short journal bearing, (ii) laminar fluid flow, (iii) iso-viscous fluid, (iv) incompressible fluid, and (v) π-film cavitation model, the Reynolds equation is written as follows 20 : where, is the rotational speed of the rotor (rad/s) µ is the fluid viscosity (Pa.s)By integrating the Reynolds equation with respect to Z, the pressure distribution for the non-worn and worn regions are obtained 20 : Non-worn region: The pressure distribution in the non-worn region, P, is expressed by Eq. (18).
where, L is the journal bearing length (m).Z is the axial distance in the direction of Z-axis Worn region: The pressure distribution in the worn region, P w , is expressed by Eq. (19).
where, P w is the pressure in the worn region of the journal bearing; h w is the fluid film thickness in the worn region of the journal bearing.
Vol.:(0123456789) Substituting the equation for pressure distribution in the non-worn regions and worn region of fluid film journal bearing into Eq.(20) gives: Substituting Eq. ( 21) with Eqs. ( 7), ( 8), ( 12) and ( 13) and integrating the result with respect to Z gives the following expression: where, The solution of Eq. ( 22) has the following form 22 : Simpsons 1/3 rule is used to solve Eq. ( 23) and the solution can be described as: www.nature.com/scientificreports/The dynamic coefficients K ij and C ij refer to the rotating coordinate system and not to the Cartesian coordi- nate system.The relationship between the dynamic coefficients in both coordinate systems can be determined by the following transformations 5 : For the worn journal bearing, the attitude angle ϕ 0 can be determined from Eq. (26) However, in the case of non-worn bearing, ϕ 0 can be related to ǫ 0 , by the following expression 22 .

Equation of motion of a rigid rotor in journal bearing
In the formulation of the equations of motion, the rotor is assumed to be rigid, symmetric and with a centrallylocated mass.Therefore, only the lowest mode (cylindrical) can be excited by the applied unbalance.Rotor motion in the axial direction is neglected, and only the motion in the vertical (X) and horizontal (Y) directions are considered in the formulation of the equations of motion.Let the rotor mass be 2M and the journal center amplitudes x and ȳ .Hence, the equations of motion can be written 22 : Here F x and F y are the journal bearing film forces.The equations are made dimensionless by setting 23 : where modified Sommerfeld number: µ is the lubricant viscosity, ω is the angular speed, D is the bearing diameter, L is the bearing length, W is the static load on the bearing, C is the radial bearing clearance.
The modified Sommerfeld number, σ , can be related to the eccentricity ratio, ǫ 0 , by the following expression.

Effect of wear on the static equilibrium position
The effect of wear on the static equilibrium position was examined for values of ǫ 0 between 0.1 and 0.9.The wear depth parameter ratio δ was varied from 0 to 0.5, in increments of 0.1.Figure 2a shows the relationship between ϕ 0 and ǫ 0 for bearing without wear.It is observed that ǫ 0 has an inverse relationship with ϕ 0 whereby that ǫ 0 increases as ϕ 0 decreases.To understand the relationship between ϕ 0 and ǫ 0 , the relationship between σ and ǫ 0 needs to be understood first as expressed in Eq. (31).Equation (31) provides a relationship to determine ǫ 0 required to generate the fluid film reaction force.Large σ represents low operating load, high ω or high fluid viscosity, µ , which results in small ǫ 0 or nearly centered operation, ǫ 0 → 0 , ϕ 0 → π/2 .That is, the journal eccentricity vec- tor is nearly orthogonal to the applied load.Small σ represents high operating load, low ω or small µ , which results in large ǫ 0 , ǫ 0 → 1.0 , ϕ 0 → 0 .That is, the journal eccentricity vector is nearly parallel to the applied load.
The relationship between ϕ 0 and ǫ 0 is expressed in Eq. ( 27) and this equation is only valid to calculate ϕ 0 for the case of bearing without wear.However, for the case of worn bearing, ϕ 0 is calculated iteratively using Eq. ( 26) 22 .
Figure 2b shows the relationship between ϕ 0 and ǫ 0 for worn bearing.The variation of ϕ 0 is demonstrated as a function of ǫ 0 and δ .It is observed that, for a constant value of ǫ 0 , ϕ 0 decreases as δ increases.For example, ǫ 0 is set constant at 0.01 and δ is varied from 0 to 0.5.When δ = 0, ϕ 0 = 1.558 and when δ = 0.5, ϕ 0 = 1.111, giving a percentage difference of 28.72% between a non-worn and worn bearing.Moreover, the percentage difference between a non-worn and worn bearing also increases as ǫ 0 increases 24 .For example,ǫ 0 is set constant at 0. δ is varied from 0 to 0.5.When δ = 0, ϕ 0 = 0.676 and when δ = 0.5, ϕ 0 = 0.007, calculating a percentage difference of 98.95% between a non-worn and worn bearing.The static equilibrium position from the present work in the condition when there is no wear ( δ = 0 ) matches closely to those that were presented by Jamil et al. 20 .The fluid film thickness in the journal bearing, given in Eqs.
(1) and ( 9), respectively for the non-worn and worn journal bearing, is a function of the attitude angle ϕ 0 .For a given value of ǫ 0 , the attitude angle varies with the depth of wear, as illustrated in Fig. 2b.Therefore, the angular position ( γ ) where the minimum film thickness occurs changes with the depth of wear.The fluid film thickness profile for the case of ǫ 0 =0.4 is shown in Fig. 3 for various wear depth values.It is observed that as the wear depth increases, the angular position ( γ ) where the minimum film thickness occurs decreases.

Effect of wear on the dynamic coefficients of journal bearings
The effect of wear on the dynamic coefficients of journal bearings was investigated for values of ǫ 0 between 0.1 and 0.75.The wear depth parameter ratio δ was varied from 0 to 0.5.Figures 4 and 5 show the dimensionless stiffness and damping coefficient as a function of σ for (a) bearing with δ = 0 (b) bearing with δ = 0.5.It is observed that wear alters the relation between the direct and the cross coupled stiffness coefficients in both the X-direction and the Y-direction, especially at low σ .As reported by Machado and Cavalca 19 , the wear increases the anisotropic characteristics of these coefficients.The dimensionless linearized dynamic coefficient from the present work in the condition when there is no wear ( δ = 0 ) closely matches those that were presented by Jamil et al. 20 .
Figure 6 shows the dimensionless dynamic stiffness coefficient against δ ranging from δ = 0 to δ = 0.5 when the rotor is operated under high operating regime ( ǫ 0 = 0.61).It is observed that K xx increases by 38.37%, K xy reduces by 51.25%, K yx reduces by 309.72%, and K yy reduces by 75.46%.
Dynamic direct stiffness coefficient relates the change in force in one direction to a displacement in the same direction as shown in Eqs. ( 34) and ( 35) 25 : The nature of the dynamic direct stiffness is to provide a restoring force that pushes the journal back toward its steady-state equilibrium position 25 .In Fig. 6, it is seen that the dynamic direct stiffness coefficients K yy decreases with increasing δ whereas the direct stiffness coefficient K xx increases with increasing δ .This is because when δ increases, ϕ 0 decreases and ǫ 0 increases which causes the rotor to move closer to the X-axis and further in the X-direction.As such, the displacement in the X-direction, x and the restoring force in the X-direction, F x , increases whereas the displacement in the Y-direction, y , and the restoring force in the Y-direction, F y , reduces.As a result, an increase in K xx and a drop in K yy is observed.
The cross-coupled stiffness relates to a displacement resulted in force component perpendicular to this displacement as given in Eqs. ( 36) and (37) 25 : (34) www.nature.com/scientificreports/where, x is the absolute magnitude of rotor vibration response in the X-direction.ȳ is the absolute magnitude of rotor vibration response in the Y-direction.φ x is the phase angle of vibration response in the X-direction measured counter-clockwise from the X-axis.φ y is the phase angle of vibration response in the Y-direction measured counter-clockwise from the Y-axis.
It is seen that the area can be a positive or negative number with positive being for counter-clockwise whirling and negative for clockwise whirling.The energy input, E, produced by the bearing coefficients is calculated by integrating the multiplication product of force and displacement around the closed curve of the ellipse.For the cross coupled stiffness coefficients, Zeidan et al. 27 has shown that: From Eq. ( 39), it can be noted that the cross-coupling's destabilizing effect is directly proportional to the area of the whirl orbit and with the net difference between K xy and K yx .It is also noted that E = 0, if K xy and K yx are equal and have identical sign.The worst combination is when K xy and K yx are equal and opposite in sign which produces a circular orbit 27 .Rotor systems are unique in that K xy = K yx and usually K xy > 0, K yx < 0 25 .It can be observed in Fig. 6 that K xy = K yx and K xy > 0 for values of δ =0 to δ=0.1.It is also noted that, when δ ≤ 0.05 , K xy and K yx have the same positive sign and when δ > 0.05 , K xy remain positive whereas K yx became negative which makes the net difference between K xy and K yx larger for δ > 0.05 as compared to when δ ≤ 0.05 .A larger net difference between K xy and K yx contributes to a larger energy input as per Eq. (39).Also, when the signs are the same, the cross-coupling has the effect of increasing the apparent asymmetry in K xx and K yy 26 .Thus, it can be deduced that when δ ≤ 0.05 the apparent asymmetry in K xx and K yy increases (rotor-bearing system becomes more stable) whereas when δ > 0.05 the apparent asymmetry in K xx and K yy decreases (rotor-bearing system becomes more unstable).
He et al. 25 mentioned that large asymmetry of K xx and K yy is the main cause of split critical speeds and non- circular (elliptical) orbit shapes.Incorporating asymmetry into the direct stiffnesses has the well-known effect of producing flat, elongated whirl orbits in the vertical and horizontal direction.This reduces the area of the orbit and thereby lessens the energy input by cross-coupled stiffness.Sufficient asymmetry flattens the orbit completely to a straight line, and the unstable system is rendered stable as the energy input becomes zero 26 .
Figure 7 shows the dimensionless damping dynamic coefficient against δ for values of δ = 0 to δ = 0.5 when the rotor is operated under high operating regime ( ǫ = 0.61 ).It is observed that C xx decreases by 17.38 %, C xy decreases by 69.10 %, C yx decreases by 69.10 %, and C yy decreases by 37.78 %.
Damping, which is generally assumed to be positive, dissipates energy from the dynamic motion of the rotorbearing system and thus promotes stability.The direct damping coefficient is related to the change in force due to a small change in velocity as shown in Eqs. ( 40) and (41) 25 .
Due to the rotor's whirling motion, the combination of the two direct damping coefficients, C xx and C yy produces a force that is tangential to the vibration orbit.This direct damping force acts against the whirling motion, helping to retard or slow it 25 .As for cross-coupled damping, it is often found that the effect of cross-coupled stiffness outweighs the effect of cross-coupled damping.This is because the resulting energy from cross-coupled damping is identically zero.Hence, in the absence of cross-coupled stiffness, cross-coupled damping cannot drive the rotor-bearing system unstable by itself 26 .It can be concluded that in general, the effect of wear on the dynamic coefficients of journal bearings is such that the dynamic coefficients except for K xx decrease with increasing δ .This result is consistent with the work of Zeidan & Herbage 27 .

Effect of wear on the vibration response of the rotor
The numerical simulation to study the effect of wear on the vibration response of the rotor was undertaken for 3 different loaded operating regime namely low loaded operating regime ( 0.1 ≤ ǫ 0 ≤ 0.2 ), moderately loaded oper- ating regime ( 0.3 ≤ ǫ 0 ≤ 0.5 ) and highly loaded operating regime ( 0.6 ≤ ǫ 0 ≤ 0.75) 28 .The unbalance parameter Low loaded operating regime ( 0.1 ≤ ǫ 0 ≤ 0.2) The effect of wear on the vibration response of rotor when the rotor is operated under low loaded operating regime ( ǫ 0 = 0.1 ) whereby m = 2 is shown from Fig. 8 to Fig. 10. Figure 8 shows the rotor whirl orbit for low loaded operating regime ( ǫ 0 = 0.1 ); (a) with a boundary view of unit circle (b) enlarged view.Figure 8a and b show that the rotor whirl orbit increases in size with increasing δ .Figure 9 shows the vibration response in the X-direction for low loaded operating regime ( ǫ 0 = 0.1 ) for (a) δ = 0 (b) δ = 0.5 .The magnitude of X response has increased from 0.018 in Fig. 9a to 0.046 in Fig. 9b which yields a net increase of 161% when δ increases from δ = 0 to δ = 0.5 .Figure 10 shows the vibration response in the Y-direction for low loaded operating regime ( ǫ 0 = 0.1 ) for (a) δ = 0 (b) δ = 0.5 .The magnitude of Y has increased from 0.006 in Fig. 10a to 0.018 in Fig. 10b which yields a net increase of 179% when δ increases from δ = 0 to δ = 0.5.
Moderately loaded operating regime ( 0.3 ≤ ǫ 0 ≤ 0.5) When the rotor is operated under moderately loaded operating regime ( ǫ 0 = 0.4 ) whereby m = 12.5 , the effect of wear on the vibration response of rotor can be seen in Figs.11, 12 and 13. Figure 11 shows the rotor whirl orbit for moderately loaded operating regime ( ǫ 0 = 0.4 ); (a) with a boundary view of unit circle (b) enlarged view.Figure 11a and b show that the rotor whirl orbit increases in size with increasing δ .Figure 12 shows the vibration response in the X-direction for moderately loaded operating regime ( ǫ 0 = 0.4 ) for (a) δ = 0 (b) δ = 0.5 .The magnitude of X response increased from 0.056 in Fig. 12a to 0.154 in Fig. 12b, which gives a net increase of 176% when δ increases from δ = 0 to δ = 0.5 .Figure 13 shows the vibration response in the Y-direction for moderately loaded operating regime ( ǫ 0 = 0.4 ) for (a) δ = 0(b) δ = 0.5 .The magnitude of Y response increased from 0.031 in Fig. 13a to 0.037 in Fig. 13b which yields a net increase of 20.4% when δ increases from δ = 0 to δ = 0.5.
Highly loaded operating regime ( 0.6 ≤ ǫ 0 ≤ 0.75) The effect of wear on the vibration response of rotor when the rotor is operated under highly loaded operating regime ( ǫ 0 = 0.75 ) whereby m = 73 is observed in Figs.14, 15 and 16. Figure 14 shows the rotor whirl orbit for highly loaded operating regime ( ǫ 0 = 0.75 ); (a) with a boundary view of unit circle (b) enlarged view.Figure 14a and b show that the rotor whirl orbit increases in size with increasing δ .Figure 15 shows the vibration response in the X-direction for highly loaded operating regime ( ǫ 0 = 0.75 ) for (a) δ = 0 (b) δ = 0.5 .The magnitude of X response increased from 0.059 in Fig. 15a to 0.443 in Fig. 15b which yields a net increase of 650.8% when δ increases from δ = 0 to δ = 0.5 .Figure 16 shows the vibration response in the Y-direction for highly loaded operating regime ( ǫ 0 = 0.75 ) for (a) δ = 0 (b) δ = 0.5 .The magnitude of Y response dropped from 0.059 in Fig. 16a to 0.035 in Fig. 16b which yields a net drop of 41.0% when δ increases from δ = 0 to δ = 0.5.
In general, the effect of wear on the vibration response of the rotor is such that when δ increases, the vibration response of the rotor in the X-and Y-direction under all 3 loaded operating regime increases with the exception of the response in the Y-direction for the highly loaded operating regime.For a constant value of ǫ 0 , ϕ 0 decreases as δ increases from δ = 0 to δ = 0.5 .This means that the rotor is moving closer to the vertical axis and when δ increases, the fluid-film thickness increases.As the fluid-film thickness increases, the dynamic stiffness reduces and when subjected to the same unbalance force, the amplitude of the rotor's displacement increases.However, the vibration response of the rotor in Y-direction when the rotor is operating under highly loaded operating regime decreases with increasing δ .This is because as δ increases, the clearance and fluid-film thickness in the Y-direction reduces.As such, the dynamic stiffness in that direction increases and when subjected to the same unbalance force, the amplitude of the rotor's displacement reduces as δ is increased from 0 to 0.5.

Conclusion
The effect of wear in journal bearings on the behavior of a rigid rotor, particularly on the static equilibrium position, the dynamic coefficients, and the vibration response of the rotor was investigated.The effect of wear on the static equilibrium position of the rotor was examined for values of ǫ 0 between 0.1 and 0.9 with δ varied from 0 to 0.5, in increments of 0.1.The results showed that ǫ 0 has an inverse relationship with ϕ 0 , where increase in the value of ǫ 0 cause the value of ϕ 0 to decrease.The effect of wear on the dynamic coefficients of journal bearings was investigated for values of ǫ 0 between 0.1 and 0.75 with δ varied from 0 to 0.5.The numerical results showed that  the dynamic coefficients, except for K xx , decrease with increasing δ .The effect of wear on the vibration response of a rigid rotor was undertaken for three different operating regimes of the journal bearings, namely low loaded operating regime, moderately loaded operating regime and highly loaded operating regime.The unbalance parameter u was set at 0.05.The dimensionless journal mass m was varied from 2 to 73. δ was varied from 0 to 0.5 in increments of 0.1.The numerical results showed that increasing the wear depth in the journal bearings generally caused an increase in the vibration response amplitude of the rigid rotor in the X-and Y-directions, for all three cases of operating regimes, with the exception of the response amplitude of the rotor in the Y-direction for the highly loaded operating regime.

Figure 3 .
Figure 3. Fluid film thickness profile for the case of ǫ 0 = 0.4 for various wear depth values.

Figure 8 .
Figure 8. Rotor whirl orbit for low loaded operating regime ( ǫ 0 = 0.1 ) (a) with a boundary view of unit circle (b) enlarged view.

Figure 11 .
Figure 11.Rotor whirl orbit for moderately loaded operating regime ( ǫ 0 = 0.4 ) (a) with a boundary view of unit circle (b) enlarged view.

Figure 14 .
Figure 14.Rotor whirl orbit for highly loaded operating regime ( ǫ 0 = 0.75 ) (a) with a boundary view of unit circle (b) enlarged view.