A New Geometrical Design to Overcome the Asymmetric Pressure Problem and the Resulting Response of Rotor-Bearing System to Unbalance Excitation

Abstract

Improving the bearing design helps in reducing the negative consequences related to errors in installation, manufacturing, deflections under severe loading conditions, progressive wear of machine elements, and many other aspects. One of the methods of such a design improvement effort is changing the bearing profile along the bearing width to compensate for the reduction in the geometrical gap between the shaft and the bearing inner surface due to the aforementioned causes. Since in all rotating machinery, unbalance usually exists at some level, this paper deals with the response of this modified bearing to unbalanced excitation to evaluate the effectiveness of such geometrical design on the dynamic characteristics of the rotor-bearing system. The numerical solution is performed using the finite difference method by assuming Reynolds boundary conditions to determine the cavitation limits, and the 4th-order Range-Kutta method is used to determine the time responses resulting from the unbalance excitation. The time responses to this type of excitation show that the rotor-bearing with the improved geometrical design is more stable, particularly at high speeds. In addition, this modification leads to an improvement in the lubricant layer thickness and the reduction in the levels of the generated pressure between the surfaces despite the presence of large deviations from the perfectly aligned bearing system. Furthermore, the suggested geometrical design overcomes the problem of asymmetricity in the pressure field resulting from the shaft deviation to a large extent. The results of this work (the enhancement in the level of the film thickness and the improvement in the dynamic response of the system as well as the reduction of the maximum pressure value) extend the range of misalignment in which the rotor bearing systems can operate safely which represents a significant step in designing the rotor-bearing system.

1. Introduction

Journal bearings are used in industrial applications to support the rotors in rotating machineries. This includes the use of this type of bearings in compressors and various types of turbines, such as gas, hydroelectric, and steam turbines, in addition to many other applications. Such a wide range of use is attributed to their high load, damping, operating speed, and relatively low manufacturing and maintenance cost. The industrial demand for high-power output leads to a more increase in the rotating machinery speed. However, operating at a high level of speed may result in large unstable oscillations related to the self-excited vibrations. This type of instability may lead to severe damage to the bearing. Therefore, predicting the onset of instability represents an essential step in designing this type of bearing for a safe operation [1]. Unbalance excitation is one of the rotor-bearing system’s most important causes of vibration. Unbalance, in fact, is a common machinery fault that results in due to the uneven mass distribution of any rotating machine element about its axis [2]. The vibration resulting from this type of excitation usually causes noise and may decrease the life of the bearings as well as result in unsafe working conditions [3]. These important aspects of the unbalance in the rotor bearing system have resulted in considerable literature being published on this topic over the last decades. One of the early works related to the effect of unbalance on the performance of a short bearing was presented by [4], but the results of the theory of short bearing involve significant errors when the bearing length to diameter ratio of the bearing is larger than 0.5 [5]. Brancati et al. [6] also analyzed an unbalanced rotor supported by journal bearings based on the use of the short bearing theory for different values of dimensionless unbalance of the rotor. As the unbalance is inherently present in any rotor and it is important to be considered in the modeling, Adiletta et al. [7,8,9] used numerical and experimental investigations to study the effect of unbalance on rotor response and examined the conditions for which chaotic motion may exist. However, the analyses were limited to specific bearing and rotor conditions. Chang-Jiang and Chen [10] used a simplified model to study the effect of rotating unbalance on the system’s response using long journal bearings. El-Saeidy and Sticher [11] presented a formulation for the dynamic analysis of rigid rotors subject to the mass imbalance in addition to base excitations. Sghir and Chouchane [1] investigated a wide domain of rotor bearing (short bearing) conditions where the effect of unbalance on the system response is investigated in each case. Their results illustrated that, in comparison with the perfectly balanced system, the unbalance might cause periodic oscillations at different speed ranges at multiple periods of rotation and even chaotic motion. Recently, Eling et al. [12] presented an interesting work on predicting the unbalance response of rotors supported on journal bearings considering unbalance force. They explained the interaction between the rotating shaft and the oil film that can result in unstable dynamic behavior characterized by sub-synchronous motion, the oil whirl, where the whirl frequency in journal bearing is close to half the rotation speed.

In addition to the problem of unbalance excitation, the system’s general performance is also greatly influenced by the presence of misalignment. In the typical use of journal bearing, the shaft is subjected to some degree of deviation with respect to the bearing axis due to many causes, such as the shaft deflection and the errors resulting from inappropriate installation or manufacturing, as well as the bearing wear. Researchers have illustrated the negative effects of the misalignment on the journal bearing performance [13,14,15,16,17,18] in terms of reducing the thickness of the oil layer that separates the journal and the bearing and the increase in the level of the pressure field in addition to the possibility of reducing bearing life as a result of the increase in the friction coefficient. These negative consequences can be reduced by modifying the bearing profile to compensate for the reduction in the separated clearance between the journal and the bearing. This geometrical design development is investigated by researchers considering different aspects. Optimizing the bearing geometry, for example, was experimentally investigated by Nacy [19] in order to control the side leakage using chamfered bearing edges. Another attempt in this direction was performed by [20], where they used defects in the bearing geometry in order to study its effects on reducing misalignment consequences on the bearing performance. Instead of performing localized modification of the bearing profile, Strzelecki [21] changed the bearing profile over the whole length using a hyperboloidal profile in order to increase the load-carrying capacity under misalignment. Chasalevris and Dohnal [22] showed that the stability of the journal-bearing system can be improved based on the use of variable geometry. Recently, Ren et al. [23] used a numerical investigation to study the effects of profile parameters on the bearing performance, where the results showed that the quadratic profile improves the bearing performance significantly. Jamali et al. [24,25] showed that the bearing performance can be improved by using a variable bearing profile. More recently, the misalignment problem in journal bearing has also been investigated by [26,27,28]. However, all these works did not consider the effect of changing the geometrical bearing profile parameters on the system response to unbalance excitations.

This work investigates the problem of unbalanced excitation of a rotor supported on finite-length journal bearings considering variable geometrical design parameters of the bearing profile. The effects of varying these parameters on the system’s characteristics in the case of 3D shaft deviation will be examined first in terms of the resulting level of film thickness and severity of the pressure spikes resulting from this 3D form of deviation. Then, the model will evaluate these geometrical modifications’ effectiveness on the rotor-bearing system’s unbalance response under a wide range of unbalanced excitations and rotational speeds. The finite difference method is considered in the discretization scheme of the governing equations using the Reynolds boundary conditions method. On the other hand, the 4th order Range Kutta method is used for solving the journal equations of motion that required identifying the shaft center trajectories under different operating and design conditions.

2. Governing Equations and the Mathematical Model

The required equations for the solution of finite length misaligned bearing with the consideration of bearing profile modification are presented in this section, and their numerical solution will be illustrated later. Figure 1 shows schematics drawings for the adopted model of the solution, whereas Figure 1a illustrates two views for the ideal case (without misalignment) of the bearing. Figure 1c shows the misaligned model, and Figure 1c explains the geometrical modification of the bearing profile. All the variables in this figure will be explained in more detail and the related equations.

Reynolds equation is the main equation that governs the solution of the hydrodynamic problem of journal bearing. This equation and the film thickness equation are given by [17,29],

$$\frac{\partial }{\partial x}\left(\frac{\rho h^3}{12\eta }\frac{\partial p}{\partial x}\right)+\frac{\partial }{\partial z}\left(\frac{\rho h^3}{12\eta }\frac{\partial p}{\partial z}\right)=U_m \frac{\partial \rho h}{\partial x}+\frac{\partial \rho h}{\partial t}$$

(1)

$$h=c(1+{\epsilon }_r \cos \left(\theta -\varnothing \right))$$

(2)

where,

$\eta , \rho$: viscosity and density of the lubricant, respectively.

$U_m$: mean entraining velocity, which is given by $U_m=\frac{U_j+U_b}{2}$, $U_j$ and $U_b$ are the journal and bearing speed where $U_{b=0}$ and $U_j=R\omega .$

p: pressure

h: film thickness

$\varnothing$: attitude angle

t: time

c: radial clearance

${\epsilon }_r$: eccentricity ratio, which is given by, ${\epsilon }_r=e/c$ where $e$ is eccentricity distance

Equation (1) is solved using Reynolds boundary conditions where the following conditions are used [30]:

$$P=0 \mathrm{at} \theta =0$$

$$\frac{\partial P}{\partial \theta }=P=0 \mathrm{at} \theta ={\theta }_{cav}$$

where ${\theta }_{cav}$ is the cavitation angle which is determined using an iterative solution [30,31].

Using the following variables, Equations (1) and (2) can be written in dimensionless forms:

$$P=\frac{p-p_o}{6\eta \omega }\left(\frac{c^2}{R^2}\right), H=\frac{h}{c}, x=R\theta \mathrm{and} Z=\frac{z}{L}$$

where,

The dimensionless forms of Equations (1) and (2) are:

$$\frac{\partial }{\partial \theta }\left(H^3\frac{\partial P}{\partial \theta }\right)+\alpha \frac{\partial }{\partial Z}\left(H^3\frac{\partial P}{\partial Z}\right)-\frac{\partial H}{\partial \theta }=0$$

(3)

$$H=1+{\epsilon }_r \cos \left(\theta -\varnothing \right)$$

(4)

where:

$$\alpha =\frac{R^2}{L^2}=\frac{1}{4{\left(L/D\right)}^2}$$

The dimensionless supported load is given by

$$\overline{W}=\sqrt{{\overline{W_r}}^2+{\overline{W_t}}^2}$$

(5)

where,

$$\overline{W_r}={{\displaystyle \int }}_0^1{{\displaystyle \int }}_0^{{\theta }_{cav}}P \cos \theta d\theta dz$$

(6)

$$\overline{W_t}={{\displaystyle \int }}_0^1{{\displaystyle \int }}_0^{{\theta }_{cav}}P \sin \theta d\theta dz$$

(7)

and $\overline{W}=\frac{w}{6\eta \omega RL}{\left(\frac{c}{R}\right)}^2$

And the attitude angle is [32].

$$\varnothing ={\tan }^{-1}\left(\frac{W_t}{W_r}\right)$$

(8)

The 3D misalignment model (Figure 1b) is adapted from the work of the first author [17], where the following equations govern this model:

$$\begin{gathered} \Delta v\left(z\right)=\Delta v_o \left(1-2Z\right) for Z\le 1/2 \\ \Delta v\left(z\right)=\Delta v_o \left(2Z-1\right) for Z>1/2 \\ \Delta h\left(z\right)=\Delta h_o \left(1-2Z\right) for Z\le 1/2 \\ \Delta h\left(z\right)=\Delta h_o \left(2Z-1\right) for Z>1/2 \end{gathered}$$

(9)

where, and $\mathrm{Z}=z/L$ and $\mathsf{\Delta }=\delta /c$ (dimensionless variables)

This 3D misalignment model represents a more realistic model for the vertical ($\Delta v\left(z\right))$ and horizontal $\left(\Delta h\left(z\right)\right)$ deviations of the shaft.

In contrast with the ideal case of the journal bearing, the attitude angle, as well as the eccentricity in the misaligned case, are not constant along the bearing length, and they are functions of the position Z [17],

$$\varnothing \left(z\right)={\tan }^{-1}\frac{e\sin \varnothing +\delta h\left(z\right)}{e\cos \varnothing -\delta v\left(z\right)} \mathrm{for} z\le L/2$$

$$e\left(z\right)=\sqrt{{\left(e \cos \varnothing -\delta v\left(z\right)\right)}^2+{\left(e \sin \varnothing +\delta h\left(z\right)\right)}^2}$$

(10)

$$\varnothing \left(z\right)={\tan }^{-1}\frac{e\sin \varnothing -\delta h\left(z\right)}{e\cos \varnothing +\delta v\left(z\right)} \mathrm{for} z>L/2$$

$$e\left(z\right)=\sqrt{{\left(e \cos \varnothing +\delta v\left(z\right)\right)}^2+{\left(e \sin \varnothing -\delta h\left(z\right)\right)}^2}$$

where,

$\varnothing$: attitude angle at $z=L/2$

$e$: eccentricity at $z=L/2$

Using these equations, the resulting gap due to the presence of 3D misalignment can be calculated.

The variation of the bearing profile shown previously in Figure 1c will change the gap between the shaft’s surface and the bearing’s inner surface, which is essentially used to compensate for the reduction in the gap due to misalignment. Detailed steps of this method were explained in [24].

The modification of the bearing profile consists of removing material from the inner surface of the bearing in order to overcome the negative consequences of the 3D misalignment (misalignment in the vertical and horizontal directions), which mainly causes a thinning in the gap between the shaft (journal) and the bearing. This modification is governed by two parameters which are the length along the longitudinal bearing direction where the modification is performed ($b$) and the distance in the radial direction ($a$). These two parameters are illustrated in Figure 1c and will be written in a dimensionless form for generality, as will be explained later.

The resulting gap due to bearing profile modification (function of Z) can be given by,

$$\begin{array}{ll}G\left(z\right)=A \left(1-Z\frac{1}{B}\right)& \mathrm{for} Z\le B\\ G\left(z\right)=A \left(1+\frac{1}{B}\left(Z-1\right)\right)& \mathrm{for} Z\ge 1-B\\ G\left(z\right)=0& \mathrm{for} B<Z<1-B\end{array}$$

(11)

where $A$ and $B$ are:

$A$=$a/\mathrm{C}$ and $B=b/L$(a and b are illustrated in Figure 1c)

The advantages of using a variable bearing profile can be evaluated based on the A and B parameters. These parameters (dimensionless) are given in terms of the $C$ (the radial clearance) and $L$ (the width of the bearing), which gives a more obvious understanding of the required geometrical change of the profile to minimize the 3D deviation negative effects.

The total gap considering the misalignment and the variable bearing profile can be determined by coupling the related equations (Equations (4), (9) and (11)).

3. The Dynamic Characteristics of the Finite Length Bearing

The stiffness and damping coefficients can be determined using linear analysis for the system’s stability. This analysis involves linearizing the nonlinear forces around the equilibrium position of the journal.

The Reynolds equation with time depending term is used to derive these coefficients as follows,

$$\frac{\partial }{\partial x}\left(\frac{h^3}{12\eta }\frac{\partial p}{\partial x}\right)+\frac{\partial }{\partial z}\left(\frac{h^3}{12\eta }\frac{\partial p}{\partial z}\right)=\frac{U_j}{2} \frac{\partial h}{\partial x}+\frac{\partial h}{\partial t}$$

(12)

The corresponding film thickness equation is [33]:

$$h=h_0+\Delta x\cos \theta +\Delta y\sin \theta$$

(13)

The term $\frac{\partial h}{\partial t}$ in Equation (12) can therefore be written as,

$$\frac{\partial h}{\partial t}=\Delta \dot{x}\cos \theta +\Delta \dot{y}\sin \theta$$

(14)

Using dimensionless forms for the variables and using Equations (12) and (14) gives,

$$\frac{\partial }{\partial \theta }\left(H^3\frac{\partial P}{\partial \theta }\right)+\alpha \frac{\partial }{\partial Z}\left(H^3\frac{\partial P}{\partial Z}\right)=\frac{\partial H}{\partial \theta }+2(\Delta \dot{Y}\sin \theta +\Delta \dot{X }\cos \theta )$$

(15)

where,

$$\dot{X}=\frac{R\dot{x}}{Uc} , \dot{Y}=\frac{R\dot{y}}{Uc}$$

(16)

The hydrodynamic forces are functions of the displacement and the velocity in the x and y directions [32,34], which are given by,

$$F_x=F_x (x, y, \dot{x}, \dot{y})$$

$$F_y=F_y (x, y, \dot{x}, \dot{y})$$

$$F_x=\underset{0}{\overset{1}{{\displaystyle \int }}}\underset{0}{\overset{{\theta }_{cav}}{{\displaystyle \int }}}P\cos \theta d\theta dZ$$

$$F_y=\underset{0}{\overset{1}{{\displaystyle \int }}}\underset{0}{\overset{{\theta }_{cav}}{{\displaystyle \int }}}P\sin \theta d\theta dZ$$

(17)

And the total force is

$$F=\sqrt{F_x{}^2+F_y{}^2}$$

The eight dynamic coefficients are [35]

$$\left[k\right]=\left[\begin{array}{cc}{k}_{xx}& k_{xy}\\ k_{yx}& k_{yy}\end{array}\right]=\left[\begin{array}{cc}\frac{\partial F_x}{\partial X}& \frac{\partial F_x}{\partial Y}\\ \frac{\partial F_y}{\partial x}& \frac{\partial F_y}{\partial Y}\end{array}\right]$$

(18)

$$\left[c\right]=\left[\begin{array}{cc}{c}_{xx}& c_{xy}\\ c_{yx}& c_{yy}\end{array}\right]=\left[\begin{array}{cc}\frac{\partial F_x}{\partial \dot{X}}& \frac{\partial F_x}{\partial \dot{Y}}\\ \frac{\partial F_y}{\partial \dot{X}}& \frac{\partial F_y}{\partial \dot{Y}}\end{array}\right]$$

(19)

The following form of the coefficients suggested by [34] is used,

$$K_{xx}=\frac{c k_{xx}}{F}, K_{xy}=\frac{c k_{xy}}{F} , K_{yx}=\frac{c k_{yx}}{F} , K_{yy}=\frac{c k_{yy}}{F}$$

(20)

$$C_{xx}=\frac{c \omega c_{xx}}{F}, C_{xy}=\frac{c \omega c_{xy}}{F} , C_{yx}=\frac{c \omega c_{yx}}{F} , C_{yy}=\frac{c \omega c_{yy}}{F}$$

(21)

Using Equations (3), (18) and (19) gives,

$$\begin{gathered} K_{xx}={{\displaystyle \int }}_0^1{{\displaystyle \int }}_0^{2\pi }{P}_x \cos \theta d\theta dz \\ {K}_{xy}={{\displaystyle \int }}_0^1{{\displaystyle \int }}_0^{2\pi }{P}_y \cos \theta d\theta dz \\ {K}_{yx}={{\displaystyle \int }}_0^1{{\displaystyle \int }}_0^{2\pi }{P}_x \sin \theta d\theta dz \\ {K}_{yy}={{\displaystyle \int }}_0^1{{\displaystyle \int }}_0^{2\pi }{P}_y \sin \theta d\theta dz \\ {C}_{xx}={{\displaystyle \int }}_0^1{{\displaystyle \int }}_0^{2\pi }{P}_{\dot{x}} \cos \theta d\theta dz \\ {C}_{xy}={{\displaystyle \int }}_0^1{{\displaystyle \int }}_0^{2\pi }{P}_{\dot{y}} \cos \theta d\theta dz \\ {C}_{yx}={{\displaystyle \int }}_0^1{{\displaystyle \int }}_0^{2\pi }{P}_{\dot{x}} \sin \theta d\theta dz \\ {C}_{yy}={{\displaystyle \int }}_0^1{{\displaystyle \int }}_0^{2\pi }{P}_{\dot{y}} \sin \theta d\theta dz \end{gathered}$$

where,

$$P_x=\frac{\partial P}{\partial X} , P_y=\frac{\partial P}{\partial Y}, P_{\dot{x}}=\frac{\partial P}{\partial \dot{X}} , P_{\dot{y}}=\frac{\partial P}{\partial \dot{Y}}$$

The following derivatives are also required in determining the eight coefficients,

$$\begin{array}{l}\frac{\partial H}{\partial t}=\Delta \dot{X}\cos \theta +\Delta \dot{Y}\sin \theta \\ \frac{\partial H}{\partial X}=\cos \theta \\ \frac{\partial H}{\partial Y}=\sin \theta \\ \frac{\partial H}{\partial \theta }=-\Delta X \sin \theta +\Delta Y\cos \theta \end{array}$$

The differentiation with respect to $X$, $Y$, $\dot{X}$ and $\dot{Y}$ yields,

$$\frac{\partial }{\partial \theta }\left(H^3\frac{\partial P_x}{\partial \theta }\right)+\alpha \frac{\partial }{\partial Z}\left(H^3\frac{\partial P_x}{\partial Z}\right)=-\frac{\partial }{\partial \theta }\left(3H^2\cos \theta \frac{\partial P}{\partial \theta }\right)-\alpha \frac{\partial }{\partial Z}\left(3H^2\cos \theta \frac{\partial P}{\partial Z}\right)-\sin \theta$$

(22)

$$\frac{\partial }{\partial \theta }\left(H^3\frac{\partial P_y}{\partial \theta }\right)+\alpha \frac{\partial }{\partial Z}\left(H^3\frac{\partial P_y}{\partial Z}\right)=-\frac{\partial }{\partial \theta }\left(3H^2\sin \theta \frac{\partial P}{\partial \theta }\right)-\alpha \frac{\partial }{\partial Z}\left(3H^2\sin \theta \frac{\partial P}{\partial Z}\right)-\cos \theta$$

(23)

$$\frac{\partial }{\partial \theta }\left(H^3\frac{\partial P_{\dot{x}}}{\partial \theta }\right)+\alpha \frac{\partial }{\partial Z}\left(H^3\frac{\partial P_{\dot{x}}}{\partial Z}\right)=\cos \theta$$

(24)

$$\frac{\partial }{\partial \theta }\left(H^3\frac{\partial P_{\dot{y}}}{\partial \theta }\right)+\alpha \frac{\partial }{\partial Z}\left(H^3\frac{\partial P_{\dot{y}}}{\partial Z}\right)=\sin \theta$$

(25)

These last four equations are solved using numerical solutions to determine the pressure derivatives used in the calculations of the eight stiffness and damping coefficients.

4. Unbalance Excitation and Linear Stability Analysis

The equations of motions of the system illustrated in Figure 2 related to the unbalance excitation can be solved after determining the dynamic coefficients. As the radial clearance between the journal and the bearing is in the order of microns, these equations for small perturbation around the static equilibrium position of the journal inside the bearing are given by [36].

$$m {\ddot{x}}^{\prime }=-F_x-f\sin \Omega t$$

(26)

$$m {\ddot{y}}^{\prime }=-F_y-f\cos \Omega t+W$$

(27)

where [36],

$\Omega$: journal speed

$F_x, F_y$: bearing forces

$f$: unbalance force

W: supported load and $x^{\prime }$ and $y^{\prime }$ are whirling axes

The dimensionless forms of Equations (26) and (27) are

$$\overline{M}{\ddot{X}}^{\prime }=-{\overline{\mathrm{F}}}_x-\overline{\mathrm{Ru}} \sin \mathrm{T}$$

(28)

$$\overline{M}{\ddot{Y}}^{\prime }=-{\overline{\mathrm{F}}}_y-\overline{\mathrm{Ru}} \cos \mathrm{T}+1$$

(29)

where,

$${\overline{F}}_x=\frac{F_x}{W} , {\overline{F}}_Y=\frac{F_Y}{W} \overline{M}=\frac{m c{ \mathsf{\Omega }}^2}{W} \mathrm{and} \overline{\mathrm{Ru}}=m_u \mathrm{r}\frac{{\mathsf{\Omega }}^2}{\mathrm{W}}$$

The response to the unbalance excitation is determined by solving these two equations.

The critical speed of the rotor is calculated as follows,

$$\overline{M}{\ddot{X}}^{\prime }+{\overline{\mathrm{F}}}_x=0$$

(30)

$$\overline{M}{\ddot{Y}}^{\prime }+{\overline{\mathrm{F}}}_y=0$$

(31)

Using linear stability analysis where the dynamic coefficients are determined after obtaining the equilibrium position of the journal center, the bearing forces are linearized as follows [33]:

$$\overline{F}x=K_{xx} X^{\prime }+K_{xy} Y^{\prime }+C_{xx} {\dot{X}}^{\prime }+C_{xy} {\dot{Y}}^{\prime }$$

(32)

$$\overline{F}y=K_{yx} X^{\prime }+K_{yy} Y^{\prime }+C_{yx} {\dot{X}}^{\prime }+C_{yy} {\dot{Y}}^{\prime }$$

(33)

Substitution of Equations (32) and (33) in Equations (30) and (31) and using the following solution [33] gives,

$$X^{\prime }=\mathrm{A}{e}^{i\lambda t} , Y^{\prime }=\mathrm{B}{e}^{i\lambda t}$$

(34)

$$\lambda =\sqrt{\frac{\left(keq-K_{xx}\right)\left(Keq-K_{yy}\right)-K_{xy} K_{yx}}{C_{xx} C_{yy}-C_{xy}{C}_{yx}}}$$

(35)

$$Keq=\frac{K_{xx} C_{yy}+K_{yy}{C}_{xx}-K_{yx}{C}_{xy}-K_{xy}{C}_{yx}}{C_{xx}+C_{yy}}$$

(36)

The critical speed (${\Omega }_{crit}$) is given by,

$${\Omega }_{crit}=\frac{\sqrt{keq}}{\lambda }$$

(37)

5. Numerical Solution

At first, the pressure field is obtained by discretizing the Reynolds equation, using the finite difference method, and the related equations for determining the gap between the two surfaces: the equations of film thickness, the 3D deviation, and the profile modification. The method of successive over-relaxation is considered in the solution under the use of the Gauss-Sedial method. These pressure fields are obtained by determining the stiffness and dynamic coefficients by numerical integration. Then, the equations of motion are solved numerically using a fourth-order Runge-Kutta method.

Discretizing the governing equations yields,

$$P_{\left(i,j\right)}=\frac{1}{\beta }\left[H_b{}^3 P_{\left(i+1,j\right)}+H_a{}^3 P_{\left(i-1,j\right)}+\alpha C_2{H}_c{}^3 P_{\left(i,j+1\right)}+\alpha C_2{H}_d{}^3 P_{\left(i,j-1\right)}-C_1{H}_{\left(i+1,j\right)}+C_1{H}_{\left(i-1,j\right)}\right]$$

(38)

$$H\left(i,j\right)=\left(1+{\epsilon }_r\left(Z\right)\cos ({\theta }_{\left(i,j\right)}-\varnothing \right))$$

(39)

where,

$C_1=\frac{\Delta \theta }{2}$, $C_2=\frac{\Delta {\theta }^2}{\Delta Z^2}$, $\beta =H_b{}^3+H_a{}^3+\alpha C_2{H}_c{}^3+\alpha C_2{H}_d{}^3$,$\Delta \theta$ and $\Delta Z$ are the mesh steps.

More details about the discretization and the solution steps can be found in [17]. After discretizing Equations (22)–(25), the pressure derivative required to calculate the dynamic coefficients is given by,

$${\overline{P}}_{\left(i.j\right)}=\frac{1}{\psi }\left[{\left(\Delta \theta \right)}^2 RHS-H_b{}^3 {\overline{P}}_{\left(i+1.j\right)}-H_a{}^3 {\overline{P}}_{\left(i-1.j\right)}-\alpha C_2{H}_c{}^3 {\overline{P}}_{\left(i.j+1\right)}-\alpha C_2{H}_d{}^3 {\overline{P}}_{\left(i.j-1\right)}+C_1{H}_{\left(i+1.j\right)}-C_1{H}_{\left(i-1.j\right)}\right]$$

(40)

where, $\psi =-H_b{}^3-H_a{}^3-\alpha C_2{H}_c{}^3-\alpha C_2{H}_d{}^3$ and the other constants are as defined previously.

The RHS of Equations (22)–(25) is calculated numerically as follows:

$$\begin{gathered} RHS\left(22\right)=\frac{\left(3\cos {\theta }_b{H}_b^2+3\cos {\theta }_a{H}_a^2\right)P_{\left(i.j\right)}}{{\left(\Delta \theta \right)}^2}-\frac{3\cos {\theta }_b{H}_b^2{P}_{\left(i+1.j\right)}}{{\left(\Delta \theta \right)}^2}-\frac{3\cos {\theta }_a{H}_a^2{P}_{\left(i-1.j\right)}}{{\left(\Delta \theta \right)}^2}+ \\ \alpha \frac{\left(3\cos {\theta }_c{H}_c^2+3\cos {\theta }_d{H}_d^2\right)P_{\left(i.j\right)}}{{\left(\Delta Z\right)}^2}-\alpha \frac{3\cos {\theta }_c{H}_c^2{P}_{\left(i.j+1\right)}}{{\left(\Delta Z\right)}^2}-\alpha \frac{3\cos {\theta }_d{H}_d^2{P}_{\left(i.j-1\right)}}{{\left(\Delta Z\right)}^2}-\sin \theta \end{gathered}$$

$$\begin{gathered} RHS\left(23\right)=\frac{\left(3\sin {\theta }_b{H}_b^2+3\sin {\theta }_a{H}_a^2\right)P_{\left(i.j\right)}}{{\left(\Delta \theta \right)}^2}-\frac{3\sin {\theta }_b{H}_b^2{P}_{\left(i+1.j\right)}}{{\left(\Delta \theta \right)}^2}-\frac{3\sin {\theta }_a{H}_a^2{P}_{\left(i-1.j\right)}}{{\left(\Delta \theta \right)}^2}+ \\ \alpha \frac{\left(3\sin {\theta }_c{H}_c^2+3\sin {\theta }_d{H}_d^2\right)P_{\left(i.j\right)}}{{\left(\Delta Z\right)}^2}-\alpha \frac{3\sin {\theta }_c{H}_c^2{P}_{\left(i.j+1\right)}}{{\left(\Delta Z\right)}^2}-\alpha \frac{3\sin {\theta }_d{H}_d^2{P}_{\left(i.j-1\right)}}{{\left(\Delta Z\right)}^2}-\cos \theta \end{gathered}$$

$$RHS\left(24\right)=\cos {\theta }_{\left(i.j\right)}$$

$$RHS\left(25\right)=\sin {\theta }_{\left(i.j\right)}$$

After discretizing all the related equations, the numerical solution is performed where the converges criterion for the pressure is, $\frac{\sum \left|P_{\left(i,j\right)new}-P_{\left(i,j\right)old}\right|}{\sum P_{\left(i,j\right)old}}<{10}^{-7}$.

The hydrodynamic load is calculated by numerical integration after determining the pressure field. This load is compared with the actual load, and the solution is accepted if an accuracy of $\mp {10}^{-5}$ is achieved ($W_{\mathrm{calculated} }$ is within the range of $W_{\mathrm{actual} }\mp {10}^{-5}{W}_{\mathrm{actual} }$). If this limit is <${10}^{-5}$, the solution requires a very long computational time without any significant change in the static and dynamic results. If this accuracy is not achieved ($W_{\mathrm{actual} }\mp {10}^{-5}{W}_{\mathrm{actual} })$, the eccentricity ratio ${\epsilon }_r$ is updated, and the whole process is repeated where the pressure distribution and the gap between the surfaces are recalculated; this step results in new values for the dynamic coefficients. After achieving both convergences (P and W), the equations of motion are solved using the 4th-order Range-Kutta method, as explained previously. A flow chart for the main steps of the solution procedure is illustrated in Figure 3.

6. Results and Discussions

The results in this work are presented considering finite length journal bearing where the length-to-diameter ratio $\left(L/D\right)$ is 1.25, and the load is corresponding to the case of a perfectly aligned bearing with a typical value of the eccentricity ratio $\left({\epsilon }_r\right)$ as 0.6. In order to better understand the geometrical design importance on the dynamic behavior of the rotor-bearing system, at first, the results will be focused on the resulting film thickness and pressure distributions due to the bearing variable geometry. Then, the investigation will be further advanced to explain in detail the effect of geometrical design parameters on the shaft trajectories under different unbalanced conditions. This sequence of presenting will give a clearer picture of the significance of the bearing profile design parameters on the system’s static and dynamic characteristics. The first step in performing the numerical solution, the independency of the results on the number of nodes, is examined as shown in Figure 4. In this figure, the values of the dimensionless minimum film thickness (H_min) and the dimensionless critical speed (CS) are plotted against the number of nodes. It can be seen that increasing the number of nodes $\left(N\times M\right)$ beyond 16,471 has a trivial effect on both results despite the narrow range of the y-axes of H_min and CS. Therefore, this number is adopted in this analysis. Furthermore, the current model is also verified using the results presented by reference [22], as shown in

Table 1. Comparison with the results of Ref. [22].

SF No.	Kyy		Kxx
	Ref. [22]	Results of the Presented Model	Ref. [22]	Results of the Presented Model
0.319	2.10	1.99	3.35	3.34
1.220	2.30	2.21	1.62	1.69

. A relatively low and high value of the Summerfield numbers is used in the comparison shown in this table. It can be seen that very good agreements are found for the values of $K_{XX}$ and $K_{YY}$ for the two Summer field numbers.

As the ideal case of a perfectly aligned bearing rarely exists in the typical applications of the journal bearings, as explained previously, the presence of misalignment is well known to have negative consequences on the general system performance. The misalignment significantly affects the film thickness, the pressure distributions, and their shape variation. Using a variable geometry for the bearing profile has the advantage of minimizing these negative effects on the system characteristics. Figure 5 shows the dimensionless film thickness and the dimensionless pressure distribution for three cases of the finite-length journal bearing. Figure 5a–c show the results for the ideal case of aligned bearing, the misaligned bearing under 3D misalignment, and the case of misaligned bearing with a modified geometry, respectively. Sever levels of misalignment are used in the case of 3D misalignment in order to illustrate how the modification can improve the bearing performance. The film thickness in the aligned case is prismatic, where the value of the film thickness at any circumferential position is constant along the bearing width. On the other hand, the pressure distribution is symmetric about the middle of the bearing width. The H_min and P_max values for this case are 0.4 and 0.645, respectively. The presence of the misalignment increases P_max by 63.8% to 1.057 and reduces H_min by 85.7% to 0.0572, representing a significant change in the designed values of the system. In addition, the film thickness distribution is no longer prismatic, and the pressure distribution is not symmetric about the middle axis of the bearing. Such a level of minim film thickness reduces the bearing life due to the expected high friction and wear ranges. In order to overcome this problem, the load-carrying capacity of the bearing should be reduced significantly to increase the value of the minimum film thickness. However, this is not an ideal solution as the bearing is usually designed for a specified supported load. Therefore, any other solution to this problem keeps the same designed load-carrying capacity without sacrificing acceptable levels of the minimum film thickness represents an important outcome. The values of the geometrical design parameters of the bearing profile can play a significant role in this direction. Improving the geometry of the bearing by removing material from the inner surface of the bearing to compensate for the reduction in the film thickness related to the misalignment has significantly positive effects on film thickness and pressure values. Figure 5c illustrates these effects where the dimensionless minimum film thickness becomes 0.244, and the maximum pressure reduces to 0.724. This represents several times improvement in the value of H_min and a decrease of 31.5% in P_max. The distribution shapes have also improved, particularly the pressure distribution shape, where the pressure spikes are reduced and shifted away from the bearing edges due to the improvement in the film thickness levels at the zones where the misalignment has the most negative effect.

The previous paragraphs emphasize the importance of changing the bearing profile parameters in enhancing the system performance in terms of reducing the value of P_max and elevating the level of the H_min for the same designed load-carrying capacity. Therefore, the investigation proceeded to examine the effect of the geometrical parameters on the main dynamic characteristics and the response of the rotor-bearing system to an unbalanced excitation. Figure 6 shows the variation of the equivalent stiffness coefficient for the modified and unmodified bearing geometry with a wide range of 3D misalignment parameters. Acceding this misalignment parameter range is insufficient for practical considerations as the minimum film thickness becomes extremely thin. It can be seen that for the ideal case (without misalignment) where $\Delta v_{max}=\Delta h_{max}=0$, the value of $K_{eq}$ is 1.347. As the misalignment parameters increase, the values of $K_{eq}$ also, increase with a maximum value of 3.161 when $\Delta v_{max}=\Delta h_{max}=0.59$. This increase means, in other words, an elevation in the level of the critical speed (which will be explained later) of the rotor bearing system, but large drawbacks accompany it in terms of H_min and P_max, as illustrated in the previous figure. Therefore, it cannot be considered an enhancement of the system’s performance. On the other hand, the geometrical modification of the profile reduces the $K_{eq}$ in comparison with the misalignment case. The reduction in $K_{eq}$ at the extreme case when $\Delta v_{max}=\Delta h_{max}=0.59$ is 32.64% ($K_{eq}=2.129)$. Despite this reduction in $K_{eq}$ as a result of introducing the geometrical change in the bearing profile, the equivalent stiffness values for the whole considered range of misalignment parameters remain higher than the stiffness of the ideal case where $\Delta v_{max}=\Delta h_{max}=0$. The improvement in $K_{eq}$ for the cases of $\Delta v_{max}=\Delta h_{max}=0.59$ is 58.1% in comparison with the ideal case. This outcome represents a significant improvement in the system performance due to the geometrical design changes, the previously mentioned reduction in P_max, and the significant increase in H_min. The modification parameters for the results presented in this figure are $A=B=0.25$. The effect of the design parameter $B$ on the equivalent stiffness is shown in Figure 7, and the optimum value of the other design parameter (parameter $A$) was found to be 0.25. It can be seen that increasing the length of modification along the axial direction of the bearing (parameter $B$) reduces the values of $K_{eq}$. However, the equivalent stiffness coefficient is not continued to change significantly when $B\ge 0.25$. This value of $B \left(0.25\right)$ also gives the best outcome in terms of Pmax and Hmin, which are therefore adopted for the main analysis in this work.

The effect of changing the geometry of the bearing on the dimensionless critical speed is shown in

Table 2. Dimensionless critical speed for the modified and unmodified system.

Case	As It Is	$\mathbf{Modified} \left(\mathit{A} = \mathit{B} = 0.25\right)$
Ideal	2.6789	2.8636
Mis. $\Delta v_{max}=\Delta h_{max}=0.5$	3.4366	3.2539
Mis. $\Delta v_{max}=\Delta h_{max}=0.59$	4.0347	3.4121

. This table compares the critical speed for two levels of misalignment which are $\Delta v_{max}=\Delta h_{max}=0.5$ and $\Delta v_{max}=\Delta h_{max}=0.59$ in addition to the perfectly aligned case. The critical speed of the ideal case is 2.6789, and the misalignment increases the critical speed by 28.28% and 50.61% (3.4366 and 4.0347) for the misalignment’s first and last levels, respectively. Introducing the bearing profile modification reduces these values, but in all modified cases, the dimensionless critical speed is stills greater than the corresponding value of the unmodified ideal bearing. The dimensionless critical speed for the misaligned and modified case when $\Delta v_{max}=\Delta h_{max}=0.59$ is 3.4121 compared with 2.6789 for the perfectly aligned case, representing an increase of 27.37%. This is another improvement in the dynamic performance of the rotor-bearing system.

After illustrating the importance of the bearing geometrical design on the key results related to the performance of the rotor-bearing system, the following results are related to the effect of unbalance excitations on the system’s dynamic behavior considering the bearing profile modification. Unbalance is a common machinery fault that results in due to the uneven distribution of a rotating component mass about its axis. Reference [36] considered a dimensionless (scaled to the load) unbalance force $\left({\overline{R}}_u=r{m}_u{\mathsf{\Omega }}^2/W\right)$ of 0.2, corresponding to an eccentricity ratio of 0.6 and a rotational speed of 3000 rpm ($\mathsf{\Omega }=314.159 \mathrm{rad}/\mathrm{s})$. This means changing the rotational speed resulting in a different unbalance force. Therefore, the following factor, $\mathsf{\Gamma }=r{m}_u/W$ is introduced to represent the unbalance force in terms of the rotational speed where ${\overline{R}}_u={\mathsf{\Gamma }\mathsf{\Omega }}^2$. Considering the value of 0.2 for ${\overline{R}}_u$ and $\mathsf{\Omega }=314.159 \mathrm{rad}/\mathrm{s}$ as used by reference [36], results in a value of $2.026423 \times {10}^{-6}$ for $\mathsf{\Gamma }$. In this work, three values for $\mathsf{\Gamma }$ are considered in the analysis, which are $0.5 \mathsf{\Gamma }$, $\mathsf{\Gamma },$ and $1.5 \mathsf{\Gamma }$ to evaluate the effect of unbalance excitation on the response of the system under different operating speeds. These values of $\mathsf{\Gamma }$ ($0.5 \mathsf{\Gamma }$, $\mathsf{\Gamma },$ and $1.5 \mathsf{\Gamma }$) mean that ${\overline{R}}_u=0.1, 0.2$ and 0.3 respectively when $\mathsf{\Omega }=3000 \mathrm{rpm} \left(314.159 \mathrm{rad}/\mathrm{s}\right)$ and give different values for ${\overline{R}}_u$ as $\mathsf{\Omega }$ changes.

Figure 8 shows the trajectories of the journal center under different rotational speeds and unbalanced excitations ($0.5 \mathsf{\Gamma }$, $\mathsf{\Gamma },$ and $1.5 \mathsf{\Gamma }$) for unmodified bearing. Figure 8a–c illustrate the trajectories when $\mathsf{\Omega }=3000 \mathrm{rpm},$ which is about half the critical speed, $\mathsf{\Omega }=0.75$ of the critical speed, and when $\mathsf{\Omega }$ is equal to the critical speed, respectively. The dimensionless critical speed of the system is 2.6789 (6540.93 rpm). It can be seen that when $\mathsf{\Omega }$ about half the critical speed of the system (Figure 8a), the journal centers rotate about the steady-state position in almost the same path of rotation for each of the three values of $\mathsf{\Gamma }$. The rotation amplitude increases as $\mathsf{\Gamma }$ increases, but in all cases, the trajectories do not cause the journal to be close to the bearing walls, and all the paths are uniform. As the rotational speed is increased, the trajectories follow different paths, as shown in Figure 8b where $\mathsf{\Omega }=0.75$ of the critical speed of the system. The higher value of $\mathsf{\Gamma }$ results in a non-safe trajectory of the journal center as the journal becomes very close to the bearing wall, and any error related to installation or the designed operating conditions may bring the surfaces of the shaft and the bearing much closer. This means in other words, that the presence of the unbalance forces in relatively high levels causes unsafe operation of the system despite that the operating speed is far from the value of the critical speed as shown in Figure 8b where $\mathsf{\Omega }=0.75$ of the critical speed. Figure 8c illustrates the trajectories when the operating speed is equal to the critical speed. It can be seen that even for the lower value of $\mathsf{\Gamma }\left(0.5\mathsf{\Gamma }\right)$, the trajectory is not uniform, and the amplitude of rotation about the steady state position varies with time. This variation is increased as $\mathsf{\Gamma }$ increased as can be seen for the case of $1.0 \mathsf{\Gamma }$ shown in blue, and the amplitude of the rotation increased rapidly for the case of 1.5$\mathsf{\Gamma }$ where the journal touches the bearing walls.

The trajectories of the journal center are further investigated by considering the bearing modification under an operating speed of 3000 rpm ($1.0 \mathsf{\Gamma }),$ as shown in Figure 9. Different values of the geometrical parameters are considered in this figure. Figure 9a shows the trajectory for the case of unmodified bearing, and the other figures illustrate the trajectories for the modified bearing where (Figure 9b): $A=B=0.1$, (Figure 9c): $A=B=0.2$, (Figure 9d): $A=B=0.25$ and (Figure 9e): $A=B=0.3$. It can be seen for this operating speed, the trajectories are uniform for all values of the modification parameters, and the journal is far away from the bearing wall. This means that modifying the bearing geometry under $1.0 \mathsf{\Gamma }$ keeps the operation of the system in the safe range.

The effect of the modification parameter in the axial direction, $B$ on the journal trajectory, is investigated for a constant value of the other modification parameter ($A=0.25$). Figure 10 shows the results of these cases when the rotational speed is 3000 rpm (1.0 $\mathsf{\Gamma }$). Figure 10a shows the case of the unmodified profile, which is repeated here from the previous figure for the purpose of comparison. Figure 10b–f illustrate the trajectories for the cases when $B=0.1$, $B=0.2$, $B=0.25, B=0.3$ and $B=0.4$. Again changing the $B$ parameter also maintains the safe operation of the system under this operating speed.

The previous figures show the trajectories relative to the bearing center using the $X,Y$ system of coordinates. To obtain a more clear comparison between the trajectories of the modified and unmodified system, the $X^{\prime }$, $Y^{\prime }$ system of coordinates is used in producing the results shown in Figure 11. This figure illustrates the journal trajectories for the unmodified and modified bearing profile relative to the steady-state position (in terms of $X^{\prime }$ and $Y^{\prime }$ coordinates) when the rotational speed is 3000 rpm (1.0 $\mathsf{\Gamma })$ and $A=B=0.25$. It can be seen that the two trajectories are very similar, with a slightly higher amplitude of the rotation in the case of modified bearing. However, as illustrated in the previous two figures, the two trajectories keep the journal surface away from the bearing wall.

Figure 12 illustrates another comparison between the journal trajectories for the unmodified and modified bearing profile relative to the steady-state position (in terms of $X^{\prime }$ and $Y^{\prime }$ coordinates) when the rotational speed is equal to the critical speed of the modified system, where the dimensionless critical speed is 2.8636 (6991.917 rpm) and the modification parameters are $A=B=0.25$. These results corresponding to the case when 1.0 $\mathsf{\Gamma }$. Despite that both trajectories are non-uniform due to the relatively high operating speed, the journal, in the case of modified bearing, tends to rotate around the steady state position with less amplitude for most of the time. This difference can be seen clearly in Figure 13, which illustrates the variation of $X^{\prime }$ and $Y^{\prime }$ with the time T (in rad) where T = $\Omega t$ and the units of $\Omega$ and t are rad/s and second, respectively. It can be seen that as T increased, the values of $X^{\prime }$ and $Y^{\prime }$ for the modified bearing are less than the corresponding values for the unmodified bearing for most of the duration of the response to the unbalance excitation.

It is clear from the above results that the modification helps in elevating the level of H_min and reducing the levels of P_max resulting from the 3D misalignment case and also resulting in higher critical speed in comparison with the ideally aligned bearing system. Furthermore, analyzing the system response under an unbalance excitation showed that the modification does not cause higher transient amplitude for the rotation around the steady state position at the high operating speed of rotation. This important outcome can be further advanced to consider the 3D misalignment in analyzing the response to unbalance excitation, which is an extremely complex case as the model requires a solution in space for the equations of motions. The authors intend to address this important issue in the near future work as both the misalignment and the unbalance excitation are common machinery problems.

7. Conclusions and Remarks

This research evaluates the effectiveness of novel bearing profile geometrical design parameters on the response of a rotor-bearing system to an unbalance excitation. A finite-length journal bearing operating at a range of speeds and a range of unbalanced excitation is considered in this analysis. The model considers a 3D misalignment case in evaluating the geometrical design parameters on the pressure and film thickness levels where the Reynolds boundary condition is used in the numerical solution to determine the cavitation zone limits. It has been found that the unbalance excitation causes the non-uniform transient path of the journal center trajectory around the steady state position, and the amplitude of the excitation increases with the increase of the operating speed as well as the increasing in the magnitude of the unbalance excitation. An improved bearing profile achieved through the optimization of the geometrical design parameters results in lower pressure values and higher minimum film thickness despite the high level of misalignment. Furthermore, the dynamic behavior is also enhanced in comparison with the ideal case of the bearing, where the equivalent stiffness is improved in addition to the improvement in the time response of the system to the unbalance excitation at high rotational speeds. The outcome of this work shows the possibility of improving the dynamic behavior of the rotor bearing system and enhancing the lubricant layer thickness, which is significantly affected by the presence of the misalignment, by optimizing the bearing profile parameters. Such improvements are expected to positively impact the bearing life and the design limits of the rotor bearing system, which will be further investigated in future work.