Analysis of the Performance of Chamfered Finite-Length Journal Bearings under Dynamic Loads

Abstract

Misalignment is one of the most common challenges that the normal operation of journal bearings faces. This type of problem may be the result of a wide range of reasons, such as bearing wear, shaft deformation, and errors related to the manufacturing and installation process. The main undesirable consequences of the misalignment, such as pressure rise and lubricant film reduction, are concentrated on the bearing edges. Therefore, chamfering the bearing edges reduces such misalignment-related drawbacks. This work presents a novel numerical solution to the problem of finite-length journal bearing considering edge chamfering. This solution involves the determination of the levels of lubricant layer thickness and pressure distribution in addition to the journal trajectory under impact load with the related stability limits. The finite difference method is used in this solution, and the equations of motion are also solved numerically using the Runge–Kutta method. The Results of this novel analysis show that chamfering the bearing edges increases the film thickness and reduces pressure spikes associated with the system operation under the case of 3D misalignment. Furthermore, the chamfered bearing shows a wide stability range under impact loads, where the normal bearing is unstable as the critical speed increases by 26.98%, which has positive consequences on the journal’s trajectory.

1. Introduction

Journal bearing is the most common type of bearing as a result of its characteristics in terms of load capacity, cost of manufacturing as well as maintenance, dynamic characteristics, etc. It is used in a machine that rotates at high speeds, such as the different types of turbines (hydroelectric, steam, and gas turbines), compressors, and others. The performance of these rotating types of machinery depends to a large extent on the performance of the bearing system. The demand for high power output in today’s applications increases the rotating machinery speed, increasing the instability problems in rotors supported on such bearings. Under high speeds, the journal bearing is subjected to dynamic instability, such as oil whirl and whip [1]. Shaft whipping, or whip, means that the shaft vibrates with a constant frequency when the rotation speed is double the critical speed [2]. This vibration displaces the center of the journal from its equilibrium position. Due to unbalanced forces, the center of the journal begins to whirl from the equilibrium position in the clearance of the bearing space and may cause instability in the journal bearing. In general, considerations of stability are significant in high-speed applications.

In the aligned case, the bearing and shaft axes remain parallel, despite the high rotational speed as well as the large load. This ideal situation is perfect, but in fact, it is difficult to obtain in the case of the actual use of journal bearings. During the operation of the bearing system, the rotating shaft will be under misalignment as a result of several causes, such as its large deformation under load, manufacturing errors in the designed tolerances, asymmetricity in the bearing load, etc.

It is well known that journal bearing performance is generally affected by the degree of misalignment. Under high levels of misalignment, the lubricant layer decreases significantly at the edges, which increases the wear rate of the bearing. The presence of misalignment increases critical speed to a certain limit [3], but on the other hand, it results in a significant increase in the maximum pressure and a sharp reduction in the thickness of the lubricant layer [4,5]. In other words, misalignment reduces the load-carrying capacity of the bearing with respect to a design limit for the film thickness. Therefore, the designer needs a critical balance between the stability limits and the designed load capacity of the bearing.

One of the experimental attempts of bearing chamfering was performed by Nacy [6] in order to reduce the side flow. Bouyer and Fillon [7] solved numerically a thermo-hydrodynamic misaligned journal bearing problem. They proposed the use of designed defects to reduce the misalignment drawbacks where a 60% improvement in the film thickness level was obtained. Strzelecki [8] presented a solution considering bearings with a variable axial profile. The modification was only performed over the whole bearing length. Chasalevris and Dohnal [9] reduced the vibration amplitude by using variable bearing geometry. Later, Chasalevris and Dohnal [10] analyzed the problem using two partial arcs bearing adjustable dynamic coefficients and varied the geometry to improve the turbine stability limits. The effect of bearing profile was investigated by Liu et al. [11] where the study was extended to consider the elasto-hydrodynamic regime in a journal bearing. Several studies proposed chamfering the bearing edges to reduce the misalignment effects without considering its impact on the system stability [5,12].

The numerical solution of journal bearing problems to identify stability limits was largely covered over the last few decades. Majumdar and Brewe [13], for example, analyzed the stability characteristics of the rigid rotor supported by journal bearings using the finite difference method. Pai and Majumdar [14] studied the stability of journal bearings under dynamic loads. Qiu and Tieu [15] solved Reynolds equation using the finite difference method under static loads where the static and dynamic characteristics as well as the stability limits of the bearings were calculated considering the effect of misalignment. Tieu and Qiu [16] compared the stability contour found based on linear theory and from a nonlinear simulation. Their results showed that the critical speeds were the same for the linear and nonlinear approaches. However, the trajectories of the journal under dynamic load were slightly different in the two solutions. In their paper, neither the effect of misalignment nor the bearing modification were considered in the solution. Papadopoulos et al. [17] investigated the stability limits of the bearings, considering the effect of misalignment. Boukhelef et al. [18] used the finite element method to calculate the dynamic and static characteristics as well as the stability of journal bearings. Nicoletti [19] used an optimization method to find bearing profiles that ensure better stability margins of rotor-bearing systems. Results showed the possibility of finding the required shapes for the bearing that considerably improved the stability margins for the case of perfectly aligned bearings. Mehrjardi et al. [20] analyzed the journal bearing stability using the 4th-order Runge–Kutta method to solve the equations of motion. Results showed that stability was enhanced by increasing the amount of bearing non-circularity. Jain and Sharma [1] analyzed journal bearings’ linear and nonlinear motions based on a numerical simulation for the journal center path. Song et al. [21] studied the influence of shaft and bearing shape errors on the stability of journal bearings. Yadav et al. [22] investigated the stability of the journal bearing using the 4th-order Runge–Kutta method in order to analyze the dynamic stability of the system. Dyk et al. [23] used an analytical method to calculate the bearing dynamic coefficients but with a limited π–film condition. Sun et al. [24] studied the effect of shaft shape errors on the dynamic characteristics of journal bearings and analyzed the stability of cylindrical journal bearings. It was found that when the shaft shape errors existed, there was a new threshold rotating speed of the system. Visnadi and Castro [25] investigated the effect of bearing clearance on the stability of short journal bearings, and the results showed that the variance in bearing clearance influenced the stability threshold. Zhang et al. [26] studied the effect of misalignment on the dynamic characteristics of journal bearings. Li et al. [27] investigated the effect of cylindricity error on the dynamic coefficients and stability of the system.

It is clear from the aforementioned work that bearing stability is quite an important subject for industry. This importance comes from today’s applications that require high levels of speed which raises problems of stability. The finite-length journal bearing is the most common type of this bearing in the industry [28]. This type generally requires numerical solutions to predict the system’s characteristics accurately. However, despite the extensive previous work on this topic, a comprehensive study for the effects of a chamfer on the stability of journal bearings is still required. Such a study will provide a clear picture of the influence of chamfer parameters on the stability of the journal bearing.

Therefore, in the present work, a numerical solution is used to solve the hydrodynamic problem of finite-length journal bearings to predict the stability limits of the rotor under classical hypotheses of hydrodynamic lubrication. Furthermore, a 3D misalignment model is incorporated into the solution scheme in order to evaluate the misalignment effects on the stability of the chamfered bearings. In addition, the presence of the dynamic load will also be examined in order to determine the journal trajectories and their reflections on the stability of the system. The stability of the journal bearing is analyzed by calculating the dynamic characteristics of the bearing, which are calculated from the pressure distribution by solving the Reynolds equation. In this work, the Reynolds equation is solved numerically by using a finite difference method. All the parameters are used in the dimensionless form in order to generalize the results. The stability analysis is performed for three cases: aligned, misaligned, and chamfered bearings. The dynamic behavior is solved by determining the journal trajectory (position of journal center) as the trajectory gives a better understanding of the dynamic behavior of the journal bearing system. The equations of motion are solved numerically by using the Runge–Kutta (4th-order) method to find the trajectories of the journal in the case of aligned bearing and chamfered bearing profiles.

2. The Mathematical Problem

The complete solution of the hydrodynamic lubrication problem of journal bearings in this work involves the following steps:

• Setting the main governing equations, i.e., the Reynolds and the film thickness equations.
• Deriving the equations required to calculate the dynamic coefficients.
• Incorporating the 3D misalignment and the chamfer equations into the film thickness equation.
• Setting the equations of motion in order to determine the critical speed and the journal trajectory equations.
• Performing a numerical solution using a finite difference method to determine the film thickness, pressure distribution, and dynamic coefficients.
• Using the Runge–Kutta method to solve the equations of motion required to determine the journal trajectory under impact loads.

3. Governing Equations

The problem of hydrodynamic journal bearings requires the solution of the Reynolds equation and the description of the film thickness function, which are represented by Equations (1) and (2), under the classical hypotheses of the isothermal hydrodynamic lubrication regime [29,30]:

$$\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)

$$\begin{array}{c}h=c(1+{\epsilon }_r \cos \left(\theta -\varnothing \right))\\ U_m=\frac{U_j+U_b}{2}\end{array}$$

(2)

where $\rho$: density, h: film thickness, $\eta$: viscosity, p: pressure, $U_m$: mean velocity, $U_j$: journal velocity, $U_b$: bearing velocity, t: time, c: radial clearance, ${\epsilon }_r$: eccentricity ratio, and $\varnothing$: attitude angle.

The bearing is stationary when $U_{b=0}$, and $U_j=R\omega$. The squeeze term is equal to zero ($\frac{\partial h}{\partial t}=0$) for the steady-state case, and $\rho$ is constant for incompressible flow.

The pressure distribution is obtained by solving Equation (2) using appropriate boundary conditions, such as the Reynolds boundary method, which requires [31]:

$$P=0$$

at

$$\theta =0$$

$$P=\frac{\partial P}{\partial \theta }=0$$

at

$$\theta ={\theta }_{cav}$$

The cavitation angle ${\theta }_{cav}$ can be determined by using an iteration method [31,32]. This involves finding the negative pressure resulting from the solutions of the Reynolds equation, starting from the last node before the boundary at the exit in the circumferential direction, setting its value to zero, and resolving the Reynolds equation based on this new condition. This process continues until all the pressure values become positive (or zero). It is worth mentioning that the identification of the locations of the negative pressure is performed in both directions $\left(z, \theta \right)$.

The Reynolds equation will be written in a dimensionless form for the purpose of the generality of the results.

The introduced dimensionless variables are:

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

where L: bearing length, R: radius, P: dimensionless pressure, p: pressure, $p_o$: the atmospheric pressure, and H: dimensionless film thickness.

Therefore, Equation (1) in a dimensionless form becomes:

$$\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)

where

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

Equation (2) describes the lubricant film thickness, and it can be written in a dimensionless form as

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

(4)

The dimensionless load components can be calculated by integrating the pressure field over the bearing surface domain [29]:

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

(5)

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

(6)

The integration can be evaluated up to the cavitation zone or continuing to include all the nodes. The pressure is zero over this zone, so the numerical integration will not be affected. However, the calculation time can be reduced by limiting the integration to consider only the area where the pressure is generated.

Therefore, the total load is:

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

(7)

where

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

The attitude angle, defined as the angle between the line of centers and the load (W) direction, can be written as [33]:

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

(8)

In the dynamic analysis, the nonlinear hydrodynamic fluid film forces should be linearized around a shaft-steady position allowing the determination of four stiffness and four damping coefficients. These coefficients are necessary for investigating the stability of the shaft around the equilibrium position. These coefficients are symbolized, respectively by $K_{xx }$, $K_{xy}$, $K_{yx}$, $K_{yy}$, $C_{xx}$, $C_{xy}$, $C_{yx}$, and $C_{yy}$.

4. Calculation of Dynamic Characteristics

The stiffness and damping coefficients equations are derived in this section, based on the following time-dependent on the Reynolds equation.

$$\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}$$

(9)

Under the dynamic condition, the equation of film thickness is [34]

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

(10)

Differentiation of Equation (9) with respect to time yields

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

(11)

Substituting Equation (11) into Equation (9) gives

$$\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}{2} \frac{\partial h}{\partial x}+\Delta \dot{x}\cos \theta +\Delta \dot{y}\sin \theta$$

(12)

Equation (12) can be written in a dimensionless form as

$$\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 )$$

(13)

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

It is worth mentioning that the x-axis is downward in the work of Lund and Thomson [34] and the same direction for this axis is considered in this section. However, the results of the dynamic coefficients are reversed in order to be consistent with the adopted system of coordinates in this work.

The resultant forces are a function of $x$,$y$,$\dot{x}$, and $\dot{y}$ as given by [33,34]:

$$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$$

(14)

where the resultant force is

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

The dynamic coefficient and matrices (stiffness and damping) can be written in the following form [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]$$

(15)

$$\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]$$

(16)

The stiffness coefficients are written using the form suggested by [34]:

$$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}$$

(17)

Similarly, the damping coefficients are

$$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}$$

(18)

Therefore, the differentiation of Equation (3) and using Equations (15) and (16) yields

$$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$$

where

$$P_x=\frac{\partial P}{\partial X} ,P_y=\frac{\partial P}{\partial Y}$$

In the same way, the damping coefficients are determined from the integration of the pressure derivatives:

$$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$$

where

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

The derivatives related to the calculation of the coefficients are evaluated as follows:

$$\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$$

The differentiation with respect to $X$ and Y yields

$$\begin{array}{l}\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 \end{array}$$

(19)

$$\begin{array}{l}\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 \end{array}$$

(20)

Similarly, the differentiation with respect to $\dot{X}$ and $\dot{Y}$ gives

$$\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$$

(21)

$$\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$$

(22)

Equations (19)–(22) will be solved numerically in order to find the pressure derivatives required to calculate the dynamic coefficients $K_{xx }$, $K_{xy}$, $K_{yx}$, $K_{yy}$, $C_{xx}$, $C_{xy}$, and $C_{yx}$. The resulting general equation after discretizing these equations is given in the section of the numerical solution.

5. Stability of Journal Bearing

The shaft whirl happened to the journal center around its equilibrium position at a certain angular speed [16]. The system is stable if the critical speed (function of dynamic coefficients) is higher than the journal rotating speed. The critical speed can be determined by solving the rotor equations of the motion. The equations of motion for a rigid unbalanced rotor supported by two bearings that are shown in Figure 1 are [36].

Figure 1. Bearing system: (a) side view of the journal bearing; (b) supported rotor [36] edited.

Figure 1. Bearing system: (a) side view of the journal bearing; (b) supported rotor [36] edited.

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

(23)

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

(24)

where [16]$f_{ex}$ and $f_{ey}$ represent the external loads, $f$ is the unbalanced force, while the bearing forces are represented by $F_x\mathrm{and}{F}_y$.

The axes of whirling of the journal center at an angular speed around its equilibrium position are ($x^{\prime }$,$y^{\prime }$) as shown in Figure 1a.

Equations (23) and (24) can be written as

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

(25)

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

(26)

where

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

The following function is used in this paper to describe the external force ${\overline{\mathrm{F}}}_{ey}$ as an impact load [16]:

$${\overline{\mathrm{F}}}_{ey}={\overline{\mathrm{F}}}_{eyo}{e}^{-a^2{\left(T/\Omega -t_o\right)}^2}$$

(27)

and

$${\overline{\mathrm{F}}}_{ex}=0$$

Neglecting the external and unbalanced forces, in order to determine the critical speed, the standard linear equations are

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

(28)

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

(29)

The bearing forces can be expressed for small displacement around the equilibrium position by [16]

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

(30)

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

(31)

Substitution of Equations (30) and (31) in Equations (28) and (29) gives

$$\overline{M} {\ddot{X}}^{\prime }+K_{xx} X^{\prime }+K_{xy} Y^{\prime }+C_{xx} {\dot{X}}^{\prime }+C_{xy} {\dot{Y}}^{\prime }=0$$

(32)

$$\overline{M} {\ddot{Y}}^{\prime }+K_{yx} X^{\prime }+K_{yy} Y^{\prime }+C_{yx} {\dot{X}}^{\prime }+C_{yy} {\dot{Y}}^{\prime }=0$$

(33)

The solution for the above equations is [37]

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

Substitution of these solutions in Equations (32) and (33) yields

$$\left(keq-K_{xx}\right)\left(keq-K_{yy}\right)-{\lambda }^2{C}_{xx}{C}_{yy}-K_{xy} K_{yx}+{\lambda }^2 C_{xy}{C}_{yx}=0$$

(34)

where

$$\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

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

(37)

6. The Misalignment Effects

The model of misalignment used in the current investigation is essentially adopted from a previous work of one of the authors of this paper [30], where a detailed illustration of the misalignment model can be found. A schematic drawing for the 3D misalignment model is illustrated in Figure 2. The use of the 3D model gives a more realistic representation of the shaft axis deviations. The deviations (dimensionless) along the bearing length in both directions (horizontal and vertical) are

$$\begin{array}{c}\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{array}$$

(38)

where $\mathsf{\Delta }=\delta /c$, and $\mathrm{Z}=z/L$.

The attitude angle and the eccentricity in the case of the misaligned bearing are functions of z position, which can be given by [30]

$$\begin{array}{c}\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}\end{array}$$

(39)

$$\begin{array}{c}\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}\end{array}$$

where $e$ and $\varnothing$ are the eccentricity distance and the attitude angle at the middle of the bearing width $\left(z=L/2\right),$ respectively. The eccentricity ratio ${\epsilon }_{r }$ is determined by dividing the eccentricity, $e\left(z\right),$ by the radial clearance, $c$. The effect of misalignment is incorporated into the film thickness, Equation (4), by calculating the new radial gap between the journal and the bearing at any $\theta ,Z$ position. This calculation can be achieved by determining the attitude angle and the eccentricity ratio corresponding to the misalignment parameters given by Equation (39).

7. Profile Modification of the Bearing

The chamfered bearing is schematically illustrated in Figure 3. More detail about this modification can be found in the previous reference [12]. The gap in the circumferential direction resulted from the bearing chamfer, which is a function of z along the bearing length as [12]

$$\begin{array}{c}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\end{array}$$

(40)

$$G\left(z\right)=0\mathrm{for}B\le Z\le 1-B$$

The dimensionless parameters $A$ and $B$ in these equations are

$$A =a/\mathrm{C}\mathrm{and}B=b/L$$

The values of these parameters are used to evaluate the adopted method of bearing chamfering. Scaling the design parameters to the clearance, $c$, and the bearing length, $L$, gives a clearer picture of the amount of change in terms of the bearing design parameters.

Figure 3. Bearing modification: (a) 3D unmodified bearing, (b) 3D modified bearing, and (c) schematic drawing of a chamfered bearing.

Figure 3. Bearing modification: (a) 3D unmodified bearing, (b) 3D modified bearing, and (c) schematic drawing of a chamfered bearing.

The coupling of Equations (4), (38), and (40) gives the total gap between the misaligned shaft and the chamfered bearing, and the resulting gap is used in the analyses of static and dynamic characteristics described previously.

8. Numerical Solution

The dynamic coefficients are calculated numerically by solving the governing equations (film thickness equation, Reynolds equation, and the equations related to the misalignment and modification of the bearing) mentioned previously. After determining the dynamic coefficients, the equations of motion used to determine the stability limits and the journal trajectory are solved numerically using a fourth-order Runge–Kutta method. The solution plane is assumed to be a rectangle that contains nodes distributed in axial and circumferential directions. The Gauss–Sedial method is used to determine the pressure distribution using the successive over-relaxation method to reduce the number of iterations.

Discretizing the governing equations gives the following equations:

$$\begin{array}{l}{P}_{\left(i,j\right)}=\frac{1}{\beta }[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)}]\end{array}$$

(41)

where $\alpha =\frac{R^2}{L^2}$, $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 Z$ and $\Delta \theta$ are the mesh steps.

The oil film thickness in terms of $i, j$ position is

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

(42)

More details about the discretizing scheme can be found in [5,30].

The convergence of the numerical solution is achieved by satisfying the convergence of the pressure values, which is given by $\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}$ and also by satisfying the load convergence criterion.

After obtaining the convergence of the solution based on the pressure, the bearing forces in both directions are calculated numerically, and the resultant hydrodynamic load can be calculated. The resulting load is then compared with the actual applied load to achieve an accuracy limit of $\mp {10}^{-5}$. If the difference is not within this limit, the eccentricity ratio is updated, and the pressure distribution is calculated again. This process continues until the load (as well as the pressure) convergence is obtained. This load corresponds to an eccentricity ratio of 0.6 in the perfectly aligned case. A flow chart for the solution steps is shown in Figure 4. After that, the dynamic coefficients can be calculated, and the equations of motion are solved using the Runge–Kutta method.

The numerical solution of Equations (19)–(22) is required to calculate the dynamic coefficients. Discretizing these equations gives the following general equation, which is also solved numerically as explained above, to find the corresponding pressure derivative (${\overline{P}}_{\left(i.j\right)}$):

$$\begin{array}{l}{\overline{P}}_{\left(i.j\right)}=\frac{1}{\psi }[{\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)}]\end{array}$$

(43)

where

$$\alpha =\frac{R^2}{L^2},C_1=\frac{\Delta \theta }{2},C_2=\frac{{(\Delta \theta )}^2}{{(\Delta Z)}^2}\mathrm{and}\psi =-H_b{}^3-H_a{}^3-\alpha C_2{H}_c{}^3-\alpha C_2{H}_d{}^3$$

The right side of Equations (19)–(22) can be determined numerically as follows:

$$\begin{array}{c}RHS\left(19\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 \\ RHS\left(20\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 \\ RHS\left(21\right)=\cos {\theta }_{\left(i.j\right)}\\ RHS\left(22\right)=\sin {\theta }_{\left(i.j\right)}\end{array}$$

The solution of the hydrodynamic lubrication problem with the consideration of the 3D misalignment and the bearing chamfer has faced considerable challenges. Examples of that are the determination of the cavitation zone boundaries, which requires an iterative solution, and derivation of the required equations to calculate the dynamic coefficients with the consideration of the misalignment and the chamfer models. Furthermore, the relaxation factor that achieved a fast convergence for the pressure distribution should be determined in addition to the finding of the faster algorithm for the load conveyance.

9. Results and Discussions

In the beginning, a series of tests are performed to find the optimal number of mesh points for this analysis as well as the number of time steps. Furthermore, the stiffness coefficients are compared with the results of reference [10] for the purpose of verification.

Table 1. Verification of the current work.

Summerfield No.	Kxx		Kyy
	Current Work	Reference [10]	Current Work	Reference [10]
0.10701	3.965	3.75	1.849	1.80
1.05856	1.317	1.45	2.149	2.00

shows this comparison using two values of the Summerfield number, where $L/D=0.75$. The agreement is excellent at the relatively low and high values of the Summerfield number.

The effect of the bearing chamfer on the pressure distribution is shown in Figure 5 for a finite-length bearing where $L/D=1.25$, using different misalignment and modification cases. The pressure distribution for the aligned case ($\Delta h_o=\Delta v_o=0$) is shown in Figure 5a where uniform and symmetric distribution about the line $Z=0.5$ (half of the bearing length) is the distinguishing feature of this case. Figure 5b,c show the results for the cases $\Delta h_o=\Delta v_o=0.55$ and $\Delta h_o=\Delta v_o=0$.58, respectively, where the pressure increases significantly due to these levels of misalignment and also loses its symmetricity about the middle of the bearing along the z direction. Introducing the chamfer at the bearing edges for the last extreme case ($\Delta h_o=\Delta v_o=0$.58) reduces the pressure spikes significantly as shown in Figure 5d,e. The chamfer parameters for these two figures are $A=B=0.1$ and $A=B=0.25,$ respectively. It can be seen that when $A=B=0.25,$ the pressure becomes closer to the aligned case distribution where the asymmetric distribution due to misalignment is overcome in addition to the reduction in its maximum value [12]. The maximum pressure is reduced by 26.34% due to the bearing chamfer compared to the misalignment case.

The corresponding film thickness of these cases is shown in Figure 6, where the chamfer improves the minimum film thickness by almost three times.

Figure 5. Effect of bearing chamfer on the dimensionless pressure distribution (L/D = 1.25). (a) $\Delta h_o=\Delta v_o=0$, (b) $\Delta h_o=\Delta v_o=0.55$, (c) $\Delta h_o=\Delta v_o=0.58$, (d) case (c) with chamfer ($A=B=0.1$) and (e) case (c) with chamfer ($A=B=0.25$).

Figure 5. Effect of bearing chamfer on the dimensionless pressure distribution (L/D = 1.25). (a) $\Delta h_o=\Delta v_o=0$, (b) $\Delta h_o=\Delta v_o=0.55$, (c) $\Delta h_o=\Delta v_o=0.58$, (d) case (c) with chamfer ($A=B=0.1$) and (e) case (c) with chamfer ($A=B=0.25$).

Figure 6. Three-dimensional film thickness distribution (L/D = 1.25). (a) $\Delta h_o=\Delta v_o=0$, (b) $\Delta h_o=\Delta v_o=0.58$, and (c) Case (b) with chamfer ($A=B=0.25$).

Figure 6. Three-dimensional film thickness distribution (L/D = 1.25). (a) $\Delta h_o=\Delta v_o=0$, (b) $\Delta h_o=\Delta v_o=0.58$, and (c) Case (b) with chamfer ($A=B=0.25$).

10. Effect of the Chamfer on the System Stability

As explained in the previous section, the misalignment has significant negative sequences on the system performance. The chamfer has shown a very clear positive effect on the levels of film thickness, the maximum pressure value, and the shape of the distribution. The results shown in this section focused on examining the chamfer effects on the system stability, which involve both the critical speed limits and the journal trajectory under impact loads. The chamfer effect on the critical speed is shown in

Table 2. Dimensionless critical speed and equivalent stiffness for different cases.

Case	Keq	Critical Speed	% Chang CS
Aligned	1.3468	2.6789	-
Misaligned	3.0040	3.9189	46.28
Chamfered	2.0927	3.4020	26.98

for the extreme case of misalignment where $\Delta h_o=\Delta v_o=0.58$ and chamfer parameters of $A=B=0.25$.

The critical speed of the aligned case is 2.678, and the corresponding critical speed for the misaligned case is 3.918. This increase in the critical speed is actually not beneficial as the misalignment increases the maximum pressure and reduces the minimum film thickness to an extremely thin level, as explained previously. Therefore, introducing the bearing edge chamfer should find a balance between the improvement in the pressure value (and the shape of distribution) and the film thickness levels and the improvement in the critical speed. Figure 7 gives a clearer picture of this outcome, whereas Figure 7a shows the effect of the chamfer parameter $A$ on the critical speed and the maximum pressure values. It can be seen that when $A\le 0.25$ the maximum pressure reduced significantly in comparison with the case of misalignment without chamfer ($A=0.0$). The corresponding critical speed is also decreased over this range of $A,$ but it is still larger than the critical speed of the aligned case (see previous table). When $A>0.25,$ the critical speed starts to increase, but the corresponding maximum pressure also increases, which represents a drawback for the chamfer. The same behavior can be seen in Figure 7b, where the corresponding values for H_min are illustrated (the critical speed values are repeated again for the purpose of clarity). This behavior can be explained as follows: for the misaligned case (the perfectly aligned bearing is a rarely existing case in practical applications of this type of bearing), the shaft is inclined due to misalignment, which reduces the gap between the shaft and the chamfered bearing close to the bearing edge. This situation helps generate the pressure between the two surfaces, i.e., the modified part of the bearing will support the load. On the other hand, as the chamfer height increases, the modified part will not significantly participate in supporting the load, which means that the bearing becomes shorter, which requires the pressure to be larger to support the same applied load.

It can be seen that for the chamfered bearing under the misalignment condition, when $A=B=0.25,$ the critical speed is increased by 26.98% compared to the aligned case. The chamfered bearing will therefore operate in a stable condition over a wider range of speed compared with the unchamfered bearing. Keeping in mind the previously mentioned improvement in performance, such a gain in critical speed represents an additional benefit and an important step in the right direction. This means, in other words, that decreasing the maximum pressure and increasing the minimum film thickness can maintain an acceptable improvement in the critical speed when the misaligned bearing is chamfered.

The variation in the critical speed occurred due to the misalignment, and the bearing chamfer is a result of the change in the dynamic coefficients of the bearing system. The dynamic coefficients for the cases illustrated above are shown in Figure 8. Figure 8a,b show the variations in the stiffness and damping coefficients, respectively. The results in Figure 8 illustrate that the chamfer returns, to some extent, the values of the coefficients back towards their values in the perfectly aligned case.

Figure 7. Effect of chamfer parameter A on the dimensionless (a) critical speed and Pmax and (b) critical speed and Hmin. (B = 0.25, L/D = 1.25, $\Delta h_o=\Delta v_o=0.58, \mathrm{and} {\epsilon }_r=0.6$).

Figure 7. Effect of chamfer parameter A on the dimensionless (a) critical speed and Pmax and (b) critical speed and Hmin. (B = 0.25, L/D = 1.25, $\Delta h_o=\Delta v_o=0.58, \mathrm{and} {\epsilon }_r=0.6$).

Figure 8. Variation of the dynamic coefficients due to misalignment and the bearing chamfer: (a) stiffness coefficients and (b) damping coefficients. (${\epsilon }_r=0.6,L/D=1.25,$ and for the modified and chamfered cases: $A=B=0.25$, $\Delta h_o=\Delta v_o=0.58$ ).

Figure 8. Variation of the dynamic coefficients due to misalignment and the bearing chamfer: (a) stiffness coefficients and (b) damping coefficients. (${\epsilon }_r=0.6,L/D=1.25,$ and for the modified and chamfered cases: $A=B=0.25$, $\Delta h_o=\Delta v_o=0.58$ ).

When the journal is hit with an impact hummer, impulse force is applied to the rotor. The trajectory of a journal relative to this force can reflect the stability of the system. The journal trajectories under impact load for the unchamfered and chamfered bearings are illustrated in Figure 9, where $L/D=1.25$, ${\epsilon }_r$= 0.6, and $A=B=0.25$ in the case of the chamfered bearing. In Figure 9a, the value of speed is 2500 rpm which is less than the critical speed of the unchamfered bearing. It can be seen that the journal returns to its steady-state position after its deviation due to the impact load. This case represents a stable condition. Similar behavior can be seen on the right side of this figure for the chamfered bearing, as the rotation speed is also less than its critical speed. Figure 9b shows the trajectories when the journal speed is 2.67 (5664.61 rpm), which is equal to the critical speed $({\Omega }_{crit})$ of the unchamfered bearing. It can be seen that at this speed, the shaft of the unchamfered bearing continues to rotate around its steady-state position on the same orbit. This trajectory represents the critical case. On the other hand, the chamfered bearing, where its journal trajectory is illustrated on the right side of Figure 9b, is stable, and the shaft returns to its original location.

In Figure 9c, ${\Omega }_{crit}$ is 3.40 (6055.17 rpm), which is equal to the critical speed of the chamfered bearing system. The shaft of the chamfered bearing rotates around the equilibrium position in a fixed orbit, but the chamfered journal bearing is unstable as the amplitude increases with time and finally reaches the bearing wall. It is clear that chamfering the bearing increases the speed range in which the journal can rotate in stable conditions under impact load.

The trajectory of the 3D misalignment journal under impact load and the bearing chamfer is an extremely complex case as the position of the journal center varies along the bearing length, which requires a different method of solution as the equation of motion should consider this variation. However, the authors intend to address this case in the next work on this path of research.

Figure 9. Effect of chamfer on the journal trajectory subjected to an impact load. Left: unchamfered, right: chamfered. (a) 2500 rpm, (b) 5664.61 rpm, (c) 6055.17 rpm. (L/D = 1.25, ${\epsilon }_r$= 0.6).

Figure 9. Effect of chamfer on the journal trajectory subjected to an impact load. Left: unchamfered, right: chamfered. (a) 2500 rpm, (b) 5664.61 rpm, (c) 6055.17 rpm. (L/D = 1.25, ${\epsilon }_r$= 0.6).

11. Conclusions

In this paper, the stability of chamfered journal bearings and the journal trajectory under an impact load are investigated in detail using a novel numerical approach. Chamfering the bearing profile and its consequences on the general performance of the journal bearing system is also presented. Furthermore, a 3D misalignment model is incorporated into the analysis. It has been found that the misalignment increases the levels of $P_{max}$ and decreases $H_{min}$ significantly, and the chamfer helps to reduce $P_{max }$ and increase $H_{min}$. In addition, the chamfer has obvious influences on ${\Omega }_{crit}$, although the ${\Omega }_{crit }$(in the case of misalignment) is greater than those of the corresponding aligned case (46.28% greater in the extreme misaligned case), and there are clear drawbacks of the misalignment on $P_{max}$ and $H_{min}$. Therefore, the chamfer is used to ensure acceptable levels of $P_{max}$ and $H_{min}$ and also to maintain the stability of the journal bearing. The chamfered bearing is stable within a range of speed in which the non-chamfered bearing is unstable. This represents an important outcome as the stability of the system is also improved in addition to the previous benefits of the chamfer, which are the reduction in $P_{max}$ and the increase in $H_{min}$. The trajectories of the journal center under impact load for the chamfered and unchamfered bearings provide a clearer understanding of the response of the system under dynamic load. It can be concluded that the bearing edge modification has significant advantages in designing the journal bearing, such as increasing the minimum film thickness, reducing the maximum pressure resulting from the misalignment, and improving in the stability of the system under impact loads.