A Study of Lubricant Inertia Effects for Squeeze Film Dampers Incorporated into High-Speed Turbomachinery

Abstract

This work proposes a numerical model that incorporates the effect of lubricant inertia on the hydrodynamic pressure distribution, fluid film reaction forces, and the fluid velocity component profiles for finite-length open-ended squeeze film dampers (SFDs). Firstly, the thin film flow equations for the SFD in presence of fluid inertia effects are introduced. Furthermore, a small first-order perturbation by means of the expressions for the fluid film velocity components and the lubricant pressure distribution that are expanded in power series of the squeeze film Reynolds number is applied to the flow equations. Subsequently the developed lubricant flow equations are solved to develop expressions for the velocity component profiles and the hydrodynamic pressure distribution in SFDs. The pressure expression is numerically solved by using Gauss–Seidel method with finite difference discretization. Moreover, the fluid film reaction forces are determined by numerically integrating the hydrodynamic pressure expressions over the journal surface. Additionally, the proposed pressure distribution expression and the numerical SFD forces are incorporated into a simulation model and the simulation results are compared with the existing models in the literature under different operating conditions, including eccentricity ratios and inertia effects (i.e., Reynolds numbers). The simulation results demonstrate the significant influence of both convective and temporal (i.e., unsteady) lubricant inertia terms on the SFD hydrodynamic pressure distribution and the fluid film reaction forces. Furthermore, the proposed SFD model is incorporated into a multi-mass flexible rotordynamic model to evaluate the effect of SFD fluid inertia on the mass unbalance induced steady-state vibrations of the rotor and the nodal transient orbits by implementing finite element method and transient modal integration with predictor–corrector solver. The results of the analysis demonstrate the significant effect of fluid inertia on the resonance frequencies of the rotor and the steady-state vibration amplitudes and the transient orbits at the resonance zone.

1. Introduction

Unbalance induced vibrations are the main source of structural vibrations in high-speed turbomachinery. This mass unbalance is associated with the limitations and imperfections in manufacturing rotor systems and leads to a synchronous load cycle in the rotor. Squeeze film dampers are essential components in high-speed turbomachinery, including aircraft jet engines, high performance compressors, gas turbines, and automotive turbochargers, that are incorporated to attenuate or completely suppress the steady-state unbalance induced vibration amplitudes at the resonance frequencies, reduce the forces transmitted to the supports, and to ensure the stable operation of the system. Figure 1 demonstrates the geometry of a conventional SFD. A typical SFD consists of a stationary outer bearing (i.e., the bush) and an inner journal with approximately identical diameters. The journal is assembled on the outer surface of a rolling element and is prevented from rotation by using an anti-rotation mechanism. The annular region between the journal and the housing is filled with a lubricant. The precession motion of the journal is induced by residual unbalance of the rotor and generates a hydrodynamic squeeze film pressure distribution that applies reaction forces over the journal, providing the damping force to attenuate the transmitted forces and in turn reducing the rotor vibration. The dynamic force response of SFDs is determined by the damper geometry, operating speed, and lubricant properties. Furthermore, SFDs do not produce direct stiffness; a suitable degree of stiffness is introduced by using retaining springs parallel to the squeeze film or by using a pair of elastomer O-rings in radial disposition.

The pioneering work by Cooper [1] demonstrated the benefits associated with the rotor operation from coupling of a damping element, represented by an oil squeeze film, to an elastic element placed within the rotor support. Ever since, research works have focused on providing greater insight into different features of SFDs, including the effect of lubricant inertia. According to classic lubrication theory, the pressure distribution in the thin film region of a hydrodynamic journal bearing is determined by using the Reynolds equation, where it is assumed that the inertial forces are negligible relative to the viscous forces (i.e., $Re\approx 0$) [2]. In modern turbomachinery, increasing velocity and size of rotors as well as application of low-viscosity lubricants requires the fluid inertia effect to be included in the design and analysis of journal bearings. The effect of fluid inertia for thin films and hydrodynamic journal bearings has been considerably studied in the literature. Smith [3] used simplified journal bearing geometries, including the short and long bearing models, to determine the effect of fluid inertia on the dynamic characteristics of journal bearings. He concluded that the effect of fluid inertia introduces an added mass to the rotor system, which significantly influences the dynamic characteristics of short stiff rotors. Typically, the proposed fluid inertia models assume circular-centered orbits (CCOs) of the journal center. Circular-centered orbits of the journal center is a very common type of journal motion in industrial applications of squeeze film dampers, including vertical rotors mounted on SFDs, horizontal rotors mounted on SFDs with centralizing springs, rotors operating close to critical velocities, and for rotor response to large unbalance forces [4]. Additionally, several investigations [5] have assumed small amplitude motions of the journal center in their corresponding investigations. According to order of magnitude calculations, for small amplitude motions of the journal center, convective inertia terms in the flow equations are negligible relative to unsteady (temporal) inertia terms [6]. However, for large amplitude motions of the journal center, including the displacement at critical speeds, the effect of convective inertia is no longer negligible and should be necessarily incorporated into the calculations.

The force coefficients are the most common representation for the SFD forces in the presence of fluid inertia. In this technique, firstly, the thin film equations are integrated into an expression for the lubricant pressure distribution by adopting either one of the momentum approximation method (aka method of averaged inertia) [7,8,9], the perturbation method [10], or the energy approximation method. Han and Rogers [11,12] have provided a detailed comparison of three approximation methods, namely momentum approximation method, iterative method, and energy method, on the force coefficients for short, long, and finite-length cylindrical SFDs. Subsequently, geometric approximations (i.e., short bearing approximation and long bearing approximation) are applied to the pressure expression to provide an approximate closed-form representation of the lubricant pressure distribution. The closed-form pressure expression is integrated over the journal surface to obtain the fluid film reaction forces (i.e., the tangential and radial components). These forces are nonlinear functions of the velocity and the acceleration of the journal center. The Jacobian matrices of the journal forces with respect to the velocities and the accelerations of the journal center are computed to develop equivalent inertia and damping coefficients respectively. According to the momentum approximation method [13,14], firstly, the flow equations are integrated along the film thickness to represent the flow dynamics in terms of the mean flows, averaged inertia, and wall shear stress differences. Subsequently, the terms in the equations are approximated by assuming that the shape of the velocity profiles is not strongly influenced by the inertia forces [15]. San Andres and Vance [16] implemented the momentum approximation technique to determine the force coefficients for finite-length SFDs executing CCOs. They included convective fluid inertia as well as end seal effects in the calculations. Furthermore, they suggested that the proposed model is strictly valid for small Reynolds numbers (i.e., $Re\le 1$). San Andres and Vance [17] analytically investigated the effect of fluid inertia and turbulence on the force coefficients of short and long SFDs. They included both temporal and convective inertia in their analysis. They suggested that at small eccentricity ratios, where the effect of temporal inertia is dominant, an added mass is produced, which corresponds to the radial direct inertia coefficient, however, at large eccentricity ratios, where the effect of convective inertia is superior, this effect is completely reversed. Dousti et al. [18] developed an extended short bearing Reynolds equation that is applicable to both laminar and turbulent flow regimes. Alternatively, according to the perturbation technique [10,19,20], a small first-order perturbation by means of the expressions for the fluid film velocity components and the lubricant pressure distribution that are typically expanded in power series of the squeeze film Reynolds number is applied to the flow equations. This perturbation technique separates the flow equations into a set of inertialess equations that are characterized by using the classic Reynolds equation, and a set of first-order inertial correction equations. According to [16,17,21] it is suggested that the first-order perturbation is applicable for SFD flow regimes with $Re\le 25$. Finally, the energy approximation [14,22] firstly develops the expressions for the fluid kinetic energy. Subsequently, the Lagrange’s equations along with the Reynolds transport theorem are applied to the energy expressions to obtain the inertia forces. Similar to the momentum approximation technique, in order to develop the fluid expressions, it is assumed that the shape of the velocity profiles is not strongly influenced by the fluid inertia effects.

The force coefficient approach is particularly valuable for rotordynamic analysis, since it directly provides an accelerated estimation of the fluid film reaction forces. However, in order to model the effect of supplementary dynamic lubricant phenomenon on the damping characteristics of SFDs, including lubricant cavitation and lubricant temperature variation, it is additionally required to develop expressions for the lubricant hydrodynamic pressure distribution and the fluid film velocity components. Furthermore, the force coefficients are typically developed for limiting bearing geometries, which makes their predictions inaccurate for arbitrary bearing geometries.

Alternatively, numerical models [23,24,25,26] provide very accurate predictions for the SFD behavior. The bulk flow model [27,28,29] is an additional numerical technique that is used to study the SFD dynamics. In this technique, bulk flow variables are introduced by calculating the average lubricant velocities across the film thickness and are substituted into the flow equations. Subsequently, the bulk-flow model system of equations, including the continuity equation and the momentum transport equation, is solved for the hydrodynamic pressure and the velocity profiles by using finite volume method. The bulk flow modeling along with the finite volume numerical provide superior accuracy for the prediction of the SFD parameters; however, the process is generally computationally very expensive, especially for integration of the SFD model into rotordynamic systems, where the SFD parameters are calculated over a considerable number of iterations.

Additional analytical techniques have been proposed for calculating the closed-form pressure distribution in SFDs [4,21,26,30,31,32,33], however, either the results are strictly valid for specific SFD geometries and configurations, which makes them inapplicable to arbitrary bearing geometries, or they are extremely complicated for integration into rotordynamic models.

Furthermore, the design and application of SFDs in turbomachinery is expansively represented in [34,35]. Typically, the rotordynamic models incorporating SFDs are categorized as follows: (1) Rigid Rotor Model; (2) Simple Flexible Rotor Model; and (3) Complex Flexible Rotor Model. The rigid rotor model assumes that the flexibility of the rotor shaft is neglected and represents the equation of motion of the rotor based on the relative position, velocity, and acceleration of the center of the SFD bearing, the center of the SFD journal, and the mass center of the rotor. The dynamics of a rigid rotor incorporating SFDs and fluid-film bearings was investigated in several studies [36,37,38,39]. The simple flexible rotor model assumes that the total mass of the rotor is either concentrated in the center or in the middle plane of the rotor, or alternatively it is distributed between the rotor center of mass and the center of the bearings. The steady-state orbit whirls of a simple flexible rotor supported by SFDs was investigated in [40,41,42]. In many practical examples, the number of degrees of freedom (DOFs) represented by the simple flexible rotor model is insufficient to provide an accurate approximation of the system dynamics. In this case, the continuous rotor model is discretized and a finite number of DOFs are introduced. The dynamics of multi-mass and multi-degrees of freedom (MDOF) rotors supported with SFDs is theoretically studied in [43,44]. The precedent studies assumed that the effect of lubricant inertia on the fluid film reaction forces in SFDs is negligible and either used the complete or the approximate (i.e., long bearing and short bearing approximations) Reynolds equation to represent the SFD dynamics. For large propulsion turbines and aero-engines the operating SFD squeeze Reynolds number is moderately large, typically on the order of one to twenty [45] and the effect of fluid inertia can no longer be neglected for those applications.

This work proposes a numerical model that represents the effect of lubricant inertia on the hydrodynamic pressure distribution, fluid film reaction forces, and the fluid velocity component profiles for finite-length open-ended squeeze film dampers (SFDs). The proposed models in this work are powerful tools that provide precise and accelerated evaluation of the SFD hydrodynamic pressure distribution, velocity profiles, and fluid film reaction forces, for application in rotordynamic models as well as the study of the effect of lubricant cavitation and lubricant temperature variation on the damping effectiveness of SFDs. The following sections describe the derivation of the analytical pressure and velocity expressions. Subsequently, the proposed expressions are incorporated into a simulation model and the performance of the SFD models are evaluated under different operating parameters. Finally, the proposed SFD model is incorporated into a flexible multi-mass rotordynamics model, to determine the effect of SFD lubricant inertia on the transient orbits and the steady-state unbalance induced vibrations of rotor systems.

2. Governing Equations

The SFD configuration in this work is a symmetric damper about its mid plane with open ends (i.e., no seal). The geometry of the system is represented in Figure 2. An orthogonal Cartesian coordinate system {x,y,z} is fixed in the plane of the lubricant, where the z-axis is perpendicular to the plane of motion. Furthermore, an orthogonal Cartesian system {x′,y′,z′} translating with angular velocity Ω is introduced, where the x′-axis is perpendicular to the line connecting the centers of the inner and outer cylinders, and the y′-axis is in the direction of the minimum thickness. The angle θ′ starts from the origin of the fixed Cartesian system and the angle θ is measured at the maximum film thickness in the direction of the whirling motion. Finally, a fixed inertial coordinate system {X,Y} is defined at the center of the bearing.

The dynamic behavior of a viscous Newtonian fluid within boundaries is generally characterized by using the 3-dimensional continuity and Navier–Stokes equations as follows [2]:

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho V\right)=0,$$

(1)

$$\rho \left(\frac{\partial V}{\partial t}+V\cdot \nabla V\right)=-\nabla P+\nabla \cdot \left(\mu \nabla V\right)+\nabla \left(-\frac{2\mu }{3}\nabla \cdot V\right)+\rho g,$$

(2)

where Equation (1) is the continuity equation corresponding to the conservation of mass within the fluid boundaries; and Equations (2) correspond to the conservation of momentum within the fluid boundaries. The terms in Equations (1) and (2) are expanded as follows [46]:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u\right)}{\partial x}+\frac{\partial \left(\rho v\right)}{\partial y}+\frac{\partial \left(\rho w\right)}{\partial z}=0,$$

(3)

$$\begin{array}{l}\rho \left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}\right)=-\frac{\partial P}{\partial x}+\frac{\partial }{\partial y}\left[\mu \left(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\right)\right]+\frac{\partial }{\partial y}\left[\mu \left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)\right]\\ +\frac{2}{3}\frac{\partial }{\partial x}\left[\mu \left(\frac{\partial u}{\partial x}-\frac{\partial v}{\partial y}\right)\right]+\frac{2}{3}\frac{\partial }{\partial x}\left[\mu \left(\frac{\partial u}{\partial x}-\frac{\partial w}{\partial z}\right)\right]+\rho g,\end{array}$$

(4)

$$\begin{array}{l}\rho \left(\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}\right)=-\frac{\partial P}{\partial y}+\frac{\partial }{\partial x}\left[\mu \left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right]+\frac{\partial }{\partial z}\left[\mu \left(\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z}\right)\right]\\ +\frac{2}{3}\frac{\partial }{\partial y}\left[\mu \left(\frac{\partial v}{\partial y}-\frac{\partial u}{\partial x}\right)\right]+\frac{2}{3}\frac{\partial }{\partial y}\left[\mu \left(\frac{\partial v}{\partial y}-\frac{\partial w}{\partial z}\right)\right]+\rho g,\end{array}$$

(5)

$$\begin{array}{l}\rho \left(\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}\right)=-\frac{\partial P}{\partial z}+\frac{\partial }{\partial y}\left[\mu \left(\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z}\right)\right]+\frac{\partial }{\partial x}\left[\mu \left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)\right]\\ +\frac{2}{3}\frac{\partial }{\partial z}\left[\mu \left(\frac{\partial w}{\partial z}-\frac{\partial u}{\partial x}\right)\right]+\frac{2}{3}\frac{\partial }{\partial z}\left[\mu \left(\frac{\partial w}{\partial z}-\frac{\partial v}{\partial y}\right)\right]+\rho g.\end{array}$$

(6)

Furthermore, it is assumed that [46]:

• The body force terms are small compared to the viscous, inertia, and pressure terms.
• According to an order of magnitude analysis, the velocity gradients $\partial u/\partial y$ and $\partial w/\partial y$ are large compared to all other velocity gradients.
• The lubricant is Newtonian, incompressible (i.e., density gradient is zero), and iso-viscous (i.e., the viscosity gradient is zero).

Applying the above assumptions to Equations (3)–(6) gives:

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0,$$

(7)

$$\rho \left\{\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}\right\}=-\frac{\partial P}{\partial x}+\mu \frac{{\partial }^{2}u}{\partial y^2},$$

(8)

$$\rho \left\{\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}\right\}=-\frac{\partial P}{\partial y}+\mu \frac{{\partial }^{2}v}{\partial y^2},$$

(9)

$$\rho \left\{\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}\right\}=-\frac{\partial P}{\partial z}+\mu \frac{{\partial }^{2}w}{\partial y^2}.$$

(10)

Additionally, according to the thin film assumption in hydrodynamic lubrication theory, which is characterized by the small ratio of film thickness to the bearing’s other physical dimensions, i.e., $c<<R$, it is concluded that:

• The effect of the curvature of the film is negligible; hence a linear coordinate system is used to describe the lubricant dynamics.
• The variation of the pressure across the film is negligible (i.e., $\partial P/\partial y=0$).

Based on the preceding description, for an incompressible and iso-viscous lubricant, the flow equations in the SFD are reduced to:

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0,$$

(11)

$$\rho \left\{\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}\right\}=-\frac{\partial P}{\partial x}+\mu \frac{{\partial }^{2}u}{\partial y^2},$$

(12)

$$\frac{\partial P}{\partial y}=0,$$

(13)

$$\rho \left\{\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}\right\}=-\frac{\partial P}{\partial z}+\mu \frac{{\partial }^{2}w}{\partial y^2}.$$

(14)

Subsequently, in order to demonstrate the dependence of the fluid inertia terms on the Reynolds number, the flow equations are normalized by introducing dimensionless parameters as follows:

$$\begin{array}{c}\begin{array}{cccc}\theta ={\theta }^{\prime }-\varphi =\frac{x}{R}-\mathsf{\Omega }t,& \eta =\frac{y}{c},& \eta =\frac{y}{c},& \xi =\frac{z}{R},\end{array}\\ \begin{array}{cccc}\tau =\mathsf{\Omega }t,& \overline{u}=\frac{u}{R\mathsf{\Omega }},& \overline{v}=\frac{v}{c\mathsf{\Omega }},& \overline{w}=\frac{w}{R\mathsf{\Omega }},\end{array}\\ \begin{array}{cccc}\overline{F}=\frac{F{c}^2}{\mu \mathsf{\Omega }{R}^4},& \overline{P}=\frac{P{c}^2}{R^2\mathsf{\Omega }\mu },& H=\frac{h}{c}=1+\epsilon \cos \theta ,& \mathrm{Re}=\frac{\rho \mathsf{\Omega }{c}^2}{\mu }.\end{array}\end{array}$$

(15)

The dimensionless parameters in Equation (15) are substituted into Equations (11), (12) and (14) as follows:

$$\frac{\partial \overline{u}}{\partial \theta }+\frac{\partial \overline{v}}{\partial \eta }+\frac{\partial \overline{w}}{\partial \xi }=0,$$

(16)

$$\mathrm{Re}\left\{\frac{\partial \overline{u}}{\partial \tau }+\overline{u}\frac{\partial \overline{u}}{\partial \theta }+\overline{v}\frac{\partial \overline{u}}{\partial \eta }+\overline{w}\frac{\partial \overline{u}}{\partial \xi }\right\}=-\frac{\partial \overline{P}}{\partial \theta }+\frac{{\partial }^2\overline{u}}{\partial {\eta }^2},$$

(17)

$$\mathrm{Re}\left\{\frac{\partial \overline{w}}{\partial \tau }+\overline{u}\frac{\partial \overline{w}}{\partial \theta }+\overline{v}\frac{\partial \overline{w}}{\partial \eta }+\overline{w}\frac{\partial \overline{w}}{\partial \xi }\right\}=-\frac{\partial \overline{P}}{\partial \xi }+\frac{{\partial }^2\overline{w}}{\partial {\eta }^2}.$$

(18)

Additionally, the dimensionless velocity boundary conditions for the lubricant are defined as follows:

$$\{\begin{array}{cc}\overline{u}=0,\overline{v}=0,\overline{w}=0& \eta =0\\ \overline{u}=0,\overline{v}=\frac{\partial H}{\partial \tau },\overline{w}=0& \eta =H\end{array}.$$

(19)

Moreover, the boundary conditions for the hydrodynamic pressure in an open ended SFD are given by:

• The pressure is periodic and continuous in the circumferential direction (θ), i.e., $\overline{P}(\theta ,\xi )=\overline{P}(\theta +2\pi ,\xi )$
• The pressure equals atmospheric pressure at the axial ends of the bearing, i.e., $\overline{P}(\theta ,L/D)=\overline{P}(\theta ,-L/D)=0$
• The hydrostatic pressure must be above the liquid cavitation pressure, i.e., $\begin{array}{cc}\overline{P}\ge {\overline{P}}_{cav}& 0\le \theta \le 2\pi ,-L/D\le \xi \le L/D\end{array}$

where P_cav is the saturation pressure of the lubricant or the saturation pressure for release of entrapped gases, typically ambient pressure.

Subsequently, a small first-order perturbation by means of the expressions for the fluid film velocity components and the lubricant pressure distribution that are expanded in power series of the squeeze film Reynolds number is introduced as follows:

$$\overline{u}=\overline{u_0}+\mathrm{Re}\overline{u_1},$$

(20)

$$\overline{v}=\overline{v_0}+\mathrm{Re}\overline{v_1},$$

(21)

$$\overline{w}=\overline{w_0}+\mathrm{Re}\overline{w_1},$$

(22)

$$\overline{P}=\overline{P_0}+\mathrm{Re}\overline{P_1}.$$

(23)

The above approximation separates the pressure and the velocity components into a zeroth-order inertialess term and a first-order inertial correction component. Substituting Equations (20)–(23) into Equations (16)–(18) and assuming that the inertial expressions in the equations are approximated by using zeroth-order velocities [10] gives:

$$\frac{\partial \overline{u_0}}{\partial \theta }+\frac{\partial \overline{v_0}}{\partial \eta }+\frac{\partial \overline{w_0}}{\partial \xi }=-\mathrm{Re}\left[\frac{\partial \overline{u_1}}{\partial \theta }+\frac{\partial \overline{v_1}}{\partial \eta }+\frac{\partial \overline{w_1}}{\partial \xi }\right],$$

(24)

$$-\frac{\partial \overline{P_0}}{\partial \theta }+\frac{{\partial }^2\overline{u_0}}{\partial {\eta }^2}=\mathrm{Re}\left[\frac{\partial \overline{u_0}}{\partial \tau }+\overline{u}\frac{\partial \overline{u}}{\partial \theta }+\overline{v}\frac{\partial \overline{u}}{\partial \eta }+\overline{w}\frac{\partial \overline{u}}{\partial \xi }+\frac{\partial \overline{P_1}}{\partial \theta }-\frac{{\partial }^2\overline{u_1}}{\partial {\eta }^2}\right],$$

(25)

$$-\frac{\partial \overline{P_0}}{\partial \xi }+\frac{{\partial }^2\overline{w_0}}{\partial {\eta }^2}=\mathrm{Re}\left[\frac{\partial \overline{w_0}}{\partial \tau }+\overline{u}\frac{\partial \overline{w}}{\partial \theta }+\overline{v}\frac{\partial \overline{w}}{\partial \eta }+\overline{w}\frac{\partial \overline{w}}{\partial \xi }+\frac{\partial \overline{P_1}}{\partial \xi }-\frac{{\partial }^2\overline{w_1}}{\partial {\eta }^2}\right].$$

(26)

According to the expressions for the inertialess continuity and momentum transport equations, the left-hand side of Equations (24)–(26) is equal to zero. Consequently, the following set of zeroth-order inertialess equations and first-order inertial equations are introduced:

2.1. Zeroth-Order Equations

$$\frac{\partial \overline{u_0}}{\partial \theta }+\frac{\partial \overline{v_0}}{\partial \eta }+\frac{\partial \overline{w_0}}{\partial \xi }=0,$$

(27)

$$-\frac{\partial \overline{P_0}}{\partial \theta }+\frac{{\partial }^2\overline{u_0}}{\partial {\eta }^2}=0,$$

(28)

$$-\frac{\partial \overline{P_0}}{\partial \xi }+\frac{{\partial }^2\overline{w_0}}{\partial {\eta }^2}=0,$$

(29)

with the following velocity boundary conditions:

$$\{\begin{array}{cc}\overline{u_0}=0,\overline{v_0}=0,\overline{w_0}=0& \eta =0\\ \overline{u_0}=0,\overline{v_0}=\frac{\partial H}{\partial \tau },\overline{w_0}=0& \eta =H\end{array}.$$

(30)

2.2. First-Order Equations

$$\frac{\partial \overline{u_1}}{\partial \theta }+\frac{\partial \overline{v_1}}{\partial \eta }+\frac{\partial \overline{w_1}}{\partial \xi }=0,$$

(31)

$$\frac{\partial \overline{u_0}}{\partial \tau }+\overline{u}\frac{\partial \overline{u}}{\partial \theta }+\overline{v}\frac{\partial \overline{u}}{\partial \eta }+\overline{w}\frac{\partial \overline{u}}{\partial \xi }=-\frac{\partial \overline{P_1}}{\partial \theta }+\frac{{\partial }^2\overline{u_1}}{\partial {\eta }^2},$$

(32)

$$\frac{\partial \overline{w_0}}{\partial \tau }+\overline{u}\frac{\partial \overline{w}}{\partial \theta }+\overline{v}\frac{\partial \overline{w}}{\partial \eta }+\overline{w}\frac{\partial \overline{w}}{\partial \xi }=-\frac{\partial \overline{P_1}}{\partial \xi }+\frac{{\partial }^2\overline{w_1}}{\partial {\eta }^2},$$

(33)

with the following velocity boundary conditions:

$$\{\begin{array}{cc}\overline{u_1}=0,\overline{v_1}=0,\overline{w_1}=0& \eta =0\\ \overline{u_1}=0,\overline{v_1}=0,\overline{w_1}=0& \eta =H\end{array}.$$

(34)

The zeroth-order velocity components and pressure distribution are determined by integrating Equations (27)–(29) in the radial direction and applying the velocity and pressure boundary conditions as follows [24]:

$$\overline{u}\left(\theta ,\eta ,\xi \right)=\frac{1}{2}\frac{\partial \overline{P_0}}{\partial \theta }\left({\eta }^2-\eta H\right),$$

(35)

$$\overline{w}\left(\theta ,\eta ,\xi \right)=\frac{1}{2}\frac{\partial \overline{P_0}}{\partial \xi }\left({\eta }^2-\eta H\right),$$

(36)

$$\overline{v}\left(\theta ,\eta ,\xi \right)=-\frac{1}{2}\left[\left(\frac{{\partial }^2\overline{P_0}}{\partial {\theta }^2}+\frac{{\partial }^2\overline{P_0}}{\partial {\xi }^2}\right)\left(\frac{{\eta }^3}{3}-H\frac{{\eta }^2}{2}\right)-\frac{\partial \overline{P_0}}{\partial \theta }\frac{\partial H}{\partial \theta }\frac{{\eta }^2}{2}\right],$$

(37)

$$\frac{\partial }{\partial \theta }\left(H^3\frac{\partial \overline{P_0}}{\partial \theta }\right)+\frac{\partial }{\partial \xi }\left(H^3\frac{\partial \overline{P_0}}{\partial \xi }\right)=12\frac{\partial H}{\partial \tau },$$

(38)

where Equation (38) is referred to as Reynolds equation [2], which is adopted to determine the inertialess lubricant pressure field. Furthermore, according to Equation (38):

$$\frac{{\partial }^2\overline{P_0}}{\partial {\theta }^2}+\frac{{\partial }^2\overline{P_0}}{\partial {\xi }^2}=\frac{12}{H^3}\frac{\partial H}{\partial \tau }-\frac{3}{H}\frac{\partial \overline{P_0}}{\partial \theta }\frac{\partial H}{\partial \theta }.$$

(39)

Substituting Equation (39) into Equation (37) gives:

$$\overline{v}\left(\theta ,\eta ,\xi \right)=\frac{1}{2}\frac{\partial \overline{P_0}}{\partial \theta }\frac{\partial H}{\partial \theta }\left(\frac{{\eta }^3}{H}-{\eta }^2\right)-\left(\frac{2{\eta }^3-3{\eta }^{2}H}{H^3}\right)\frac{\partial H}{\partial \tau }.$$

(40)

Subsequently, the inertialess velocity components in Equations (35), (36) and (40) are substituted into Equations (32) and (33) to determine the first-order velocity profiles. Firstly, according to Equation (32):

$$\begin{array}{l}\frac{\partial }{\partial \tau }\left[\frac{1}{2}\frac{\partial \overline{P_0}}{\partial \theta }\left({\eta }^2-\eta H\right)\right]+\frac{1}{2}\frac{\partial \overline{P_0}}{\partial \theta }\left({\eta }^2-\eta H\right)\frac{\partial }{\partial \theta }\left[\frac{1}{2}\frac{\partial \overline{P_0}}{\partial \theta }\left({\eta }^2-\eta H\right)\right]\\ -\frac{1}{2}\left[\left(\frac{{\partial }^2\overline{P_0}}{\partial {\theta }^2}+\frac{{\partial }^2\overline{P_0}}{\partial {\xi }^2}\right)\left(\frac{{\eta }^3}{3}-H\frac{{\eta }^2}{2}\right)-\frac{\partial \overline{P_0}}{\partial \theta }\frac{\partial H}{\partial \theta }\frac{{\eta }^2}{2}\right]\frac{\partial }{\partial \eta }\left[\frac{1}{2}\frac{\partial \overline{P_0}}{\partial \theta }\left({\eta }^2-\eta H\right)\right]\\ +\frac{1}{2}\frac{\partial \overline{P_0}}{\partial \xi }\left({\eta }^2-\eta H\right)\frac{\partial }{\partial \xi }\left[\frac{1}{2}\frac{\partial \overline{P_0}}{\partial \theta }\left({\eta }^2-\eta H\right)\right]=-\frac{\partial \overline{P_1}}{\partial \theta }+\frac{{\partial }^2\overline{u_1}}{\partial {\eta }^2}.\end{array}$$

(41)

Equation (41) is integrated twice in the radial direction and is solved for the first-order circumferential velocity component. Furthermore, the integration constants are determined by applying the velocity boundary conditions in Equation (34), hence:

$$\begin{array}{l}\overline{u_1}=\left(\frac{{\eta }^2-H\eta }{2}\right)\frac{\partial \overline{P_1}}{\partial \theta }-\frac{\partial \overline{P_0}}{\partial \theta }\frac{\partial H}{\partial \tau }\left(\frac{{\eta }^6}{15H^3}-\frac{{\eta }^5}{5H^2}+\frac{{\eta }^4}{8H}+\frac{{\eta }^3}{12}-\frac{3H^2\eta }{40}\right)\\ +{\left(\frac{\partial \overline{P_0}}{\partial \theta }\right)}^2\frac{\partial H}{\partial \theta }\left(\frac{{\eta }^6}{60H}-\frac{{\eta }^5}{20}+\frac{{\eta }^{4}H}{24}-\frac{H^4\eta }{120}\right)+\frac{{\partial }^2\overline{P_0}}{\partial \tau \partial \theta }\left(\frac{{\eta }^4}{24}-\frac{{\eta }^{3}H}{12}+\frac{H^3\eta }{24}\right)\\ +\left(\frac{\partial \overline{P_0}}{\partial \theta }\frac{{\partial }^2\overline{P_0}}{\partial {\theta }^2}+\frac{{\partial }^2\overline{P_0}}{\partial \tau \partial \xi }\frac{\partial \overline{P_0}}{\partial \xi }\right)\left(\frac{{\eta }^6}{120}-\frac{{\eta }^{5}H}{40}+\frac{{\eta }^4{H}^2}{48}-\frac{H^5\eta }{240}\right).\end{array}$$

(42)

Similarly, the axial first-order velocity component is calculated as follows:

$$\begin{array}{l}\overline{w_1}=\left(\frac{{\eta }^2-H\eta }{2}\right)\frac{\partial \overline{P_1}}{\partial \xi }-\frac{\partial \overline{P_0}}{\partial \xi }\frac{\partial H}{\partial \tau }\left(\frac{{\eta }^6}{15H^3}-\frac{{\eta }^5}{5H^2}+\frac{{\eta }^4}{8H}+\frac{{\eta }^3}{12}-\frac{3H^2\eta }{40}\right)\\ +\frac{\partial \overline{P_0}}{\partial \theta }\frac{\partial \overline{P_0}}{\partial \xi }\frac{\partial H}{\partial \theta }\left(\frac{{\eta }^6}{60H}-\frac{{\eta }^5}{20}+\frac{{\eta }^{4}H}{24}-\frac{H^4\eta }{120}\right)+\frac{{\partial }^2\overline{P_0}}{\partial \tau \partial \xi }\left(\frac{{\eta }^4}{24}-\frac{{\eta }^{3}H}{12}+\frac{H^3\eta }{24}\right)\\ +\left(\frac{\partial \overline{P_0}}{\partial \theta }\frac{{\partial }^2\overline{P_0}}{\partial {\theta }^2}+\frac{{\partial }^2\overline{P_0}}{\partial \tau \partial \xi }\frac{\partial \overline{P_0}}{\partial \xi }\right)\left(\frac{{\eta }^6}{120}-\frac{{\eta }^{5}H}{40}+\frac{{\eta }^4{H}^2}{48}-\frac{H^5\eta }{240}\right).\end{array}$$

(43)

Subsequently, in order to eliminate the radial velocity component, Equation (31) is integrated along the film thickness and the velocity boundary conditions are applied:

$$\frac{\partial }{\partial \theta }\left({\displaystyle \underset{0}{\overset{H}{\int }}\overline{u_1}d\eta }\right)+\frac{\partial }{\partial \xi }\left({\displaystyle \underset{0}{\overset{H}{\int }}\overline{w_1}d\eta }\right)=0,$$

(44)

where, the first-order lubricant flows are defined as [2]:

$$\overline{q_{i_1}}={\displaystyle \underset{0}{\overset{H}{\int }}\overline{u_{i_1}}}d\eta .$$

(45)

Substituting Equation (44) into Equation (45) gives:

$$\frac{\partial \overline{q_{{\theta }_1}}}{\partial \theta }+\frac{\partial \overline{q_{{\xi }_1}}}{\partial \xi }=0.$$

(46)

Integrating the first-order velocity expressions in Equations (42) and (43) along the film thickness and substituting into Equation (46) gives:

$$\begin{array}{l}\frac{\partial }{\partial \theta }\left(H^3\frac{\partial \overline{P_1}}{\partial \theta }\right)+\frac{\partial }{\partial \xi }\left(H^3\frac{\partial \overline{P_1}}{\partial \xi }\right)=\\ \frac{\partial }{\partial \theta }\left\{-\frac{3H^7}{560}\frac{\partial }{\partial \theta }\left[{\left(\frac{\partial \overline{P_0}}{\partial \theta }\right)}^2+{\left(\frac{\partial \overline{P_0}}{\partial \xi }\right)}^2\right]-\frac{3H^6}{140}\frac{\partial H}{\partial \theta }{\left(\frac{\partial \overline{P_0}}{\partial \theta }\right)}^2+\frac{13H^4}{70}\frac{\partial H}{\partial \tau }\frac{\partial \overline{P_0}}{\partial \theta }+\frac{{\partial }^2\overline{P_0}}{\partial \tau \partial \theta }\left(\frac{H^5}{10}\right)\right\}\\ \frac{\partial }{\partial \xi }\left\{-\frac{3H^7}{560}\frac{\partial }{\partial \xi }\left[{\left(\frac{\partial \overline{P_0}}{\partial \theta }\right)}^2+{\left(\frac{\partial \overline{P_0}}{\partial \xi }\right)}^2\right]-\frac{3H^6}{140}\frac{\partial H}{\partial \theta }\frac{\partial \overline{P_0}}{\partial \theta }\frac{\partial \overline{P_0}}{\partial \xi }+\frac{13H^4}{70}\frac{\partial H}{\partial \tau }\frac{\partial \overline{P_0}}{\partial \xi }+\frac{{\partial }^2\overline{P_0}}{\partial \tau \partial \xi }\left(\frac{H^5}{10}\right)\right\}.\end{array}$$

(47)

Equation (47) provides a Reynolds-like expression that characterizes the first-order pressure distribution. The total fluid film pressure distribution is calculated by using Equation (23). Additionally, the total fluid film pressure distribution is numerically integrated over the journal surface to calculate the fluid film reaction forces as follows:

$$\left[\begin{array}{c}{\overline{F}}_r\\ \overline{F_t}\end{array}\right]={\displaystyle \underset{-L/D}{\overset{L/D}{\int }}{\displaystyle \underset{{\theta }_1}{\overset{{\theta }_2}{\int }}\overline{P}(\theta ,\xi )}}\left[\begin{array}{c}\cos (\theta )\\ \sin (\theta )\end{array}\right]d\theta d\xi .$$

(48)

where θ₁ to θ₂ is the positive range of the pressure distribution. The following section describes a numerical scheme that determines the zeroth-order and first-order pressure distributions.

3. Numerical Solution

This section describes the procedure to determine the pressure distribution and the fluid film forces in SFDs by numerically solving Equations (38) and (47). In order to determine the numerical solution, firstly, a solution domain is defined for the problem. Subsequently, the partial differential pressure expressions are discretized over the solution domain by using finite difference approximation. Finally, an iterative Gauss–Seidel numerical algorithm with successive over-relaxation (SOR) is developed to calculate the point-wise pressure distribution in the thin film domain.

3.1. Zeroth-Order Pressure Solution

In order to calculate the total fluid film pressure distribution, the zeroth-order pressure distribution must be first determined. The zeroth-order pressure distribution is characterized by Equation (38). This equation is expanded and discretized by applying finite difference approximation. The detailed description of the numerical procedure is provided in Appendix A. The discretized zeroth-order pressure expression is produced as follows:

$${\overline{P}}_0{}_{{}_{i,j}}=\frac{A_4+\left[A_1+A_2\right]{\overline{P_0}}_{{}_{i-1,j}}+\left[A_2\right]{\overline{P_0}}_{{}_{i+1,j}}+\left[A_3\right]\left({\overline{P_0}}_{{}_{i,j+1}}+{\overline{P_0}}_{{}_{i,j-1}}\right)}{\left[A_1+2A_2+2A_3\right]},$$

(49)

where

$$\begin{array}{l}{A}_1=-\frac{3H_i^2\left(\epsilon \sin {\theta }_i\right)}{\Delta \theta },\\ A_2=-\frac{H_i^3}{\Delta {\theta }^2},\\ A_3=-\frac{H_i^3}{\Delta {\xi }^2},\\ A_4=12\epsilon \sin {\theta }_i.\end{array}$$

(50)

In general, Reynolds equation is classified as an elliptical PDE. Assuming that the lubricant is incompressible and iso-viscous, and the journal center executes CCO whirls, the following Gauss–Seidel numerical procedure is used to determine the fluid film pressure distribution for a specified SFD eccentricity ratio:

• The boundary points are initialized to their prescribed values, and the interior points are adjusted to zero.
• Equation (49) is iteratively solved for the interior points.
• The iteration is only interrupted when the solution error reaches a convergence criterion.
• Finally, a SOR technique is used to accelerate the convergence of the solution:

  $$P_{0_{i,j}}^{(k)}=P_{0_{i,j}}^{(k-1)}+\lambda \left(P_{0_{i,j}}^{(k)}-P_{0_{i,j}}^{(k-1)}\right),$$

  (51)

  where k denotes the iteration number.

3.2. First-Order Pressure Solution

The numerical procedure that is employed to determine the first-order pressure distribution is very similar to the one that was described for the zeroth-order pressure. The first-order pressure distribution in the fluid film is characterized by Equation (47). Similarly, this equation is expanded and discretized by applying finite difference approximation. The detailed representation of the numerical procedure is provided in Appendix B. The discretized first-order pressure expression is produced as follows:

$$P_1{}_{{}_{i,j}}=\frac{\left[A_4+G_{1_{i,j}}+G_{2_{i,j}}\right]+\left[A_1+A_2\right]{\overline{P}}_{1_{i-1,j}}+\left[A_2\right]{\overline{P}}_{1_{i+1,j}}+\left[A_3\right]\left({\overline{P}}_{1_{i,j+1}}+{\overline{P}}_{1_{i,j-1}}\right)}{\left[A_1+2A_2+2A_3\right]}.$$

(52)

where

$$\begin{array}{cc}\hfill G_1\left(\theta ,\xi \right)=& -\frac{3H^6}{40}\frac{\partial H}{\partial \theta }\left[\frac{\partial \overline{P_0}}{\partial \theta }\frac{{\partial }^2\overline{P_0}}{\partial {\theta }^2}+\frac{\partial \overline{P_0}}{\partial \xi }\frac{{\partial }^2\overline{P_0}}{\partial \theta \partial \xi }\right]\hfill \\ & -\frac{3H^7}{280}\left[{\left(\frac{{\partial }^2\overline{P_0}}{\partial {\theta }^2}\right)}^2+\frac{\partial \overline{P_0}}{\partial \theta }\frac{{\partial }^3\overline{P_0}}{\partial {\theta }^3}+{\left(\frac{{\partial }^2\overline{P_0}}{\partial \theta \partial \xi }\right)}^2+\frac{\partial \overline{P_0}}{\partial \xi }\frac{{\partial }^3\overline{P_0}}{\partial {\theta }^2\partial \xi }\right]\hfill \\ & -\frac{9H^5}{70}{\left(\frac{\partial H}{\partial \theta }\right)}^2{\left(\frac{\partial \overline{P_0}}{\partial \theta }\right)}^2-\frac{3H^6}{140}\frac{{\partial }^{2}H}{\partial {\theta }^2}{\left(\frac{\partial \overline{P_0}}{\partial \theta }\right)}^2-\frac{3H^6}{70}\frac{\partial H}{\partial \theta }\frac{\partial \overline{P_0}}{\partial \theta }\frac{{\partial }^2\overline{P_0}}{\partial {\theta }^2}\hfill \\ & +\frac{26H^3}{35}\frac{\partial H}{\partial \theta }\frac{\partial H}{\partial \tau }\frac{\partial \overline{P_0}}{\partial \theta }+\frac{13H^4}{70}\frac{{\partial }^{2}H}{\partial \tau \partial \theta }\frac{\partial \overline{P_0}}{\partial \theta }+\frac{13H^4}{70}\frac{\partial H}{\partial \tau }\frac{{\partial }^2\overline{P_0}}{\partial {\theta }^2}\hfill \\ & +\frac{{\partial }^3\overline{P_0}}{\partial \tau \partial {\theta }^2}\left(\frac{H^5}{10}\right)+\frac{{\partial }^2\overline{P_0}}{\partial \tau \partial \theta }\frac{H^4}{2}\frac{\partial H}{\partial \theta },\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill G_2\left(\theta ,\xi \right)=& -\frac{3H^7}{280}\left[{\left(\frac{{\partial }^2\overline{P_0}}{\partial \theta \partial \xi }\right)}^2+\frac{\partial \overline{P_0}}{\partial \theta }\frac{{\partial }^3\overline{P_0}}{\partial \theta \partial {\xi }^2}+{\left(\frac{{\partial }^2\overline{P_0}}{\partial {\xi }^2}\right)}^2+\frac{\partial \overline{P_0}}{\partial \xi }\frac{{\partial }^3\overline{P_0}}{\partial {\xi }^3}\right]\hfill \\ & -\frac{3H^6}{140}\frac{\partial H}{\partial \theta }\frac{{\partial }^2\overline{P_0}}{\partial \theta \partial \xi }\frac{\partial \overline{P_0}}{\partial \xi }-\frac{3H^6}{140}\frac{\partial H}{\partial \theta }\frac{\partial \overline{P_0}}{\partial \theta }\frac{{\partial }^2\overline{P_0}}{\partial {\xi }^2}\hfill \\ & +\frac{13H^4}{70}\frac{\partial H}{\partial \tau }\frac{{\partial }^2\overline{P_0}}{\partial {\xi }^2}+\frac{{\partial }^3\overline{P_0}}{\partial \tau \partial {\xi }^2}\left(\frac{H^5}{10}\right),\hfill \end{array}$$

(54)

Subsequently, Equation (52) is iteratively solved by using the numerical procedure that was described in the previous section. Finally, the total pressure distribution is determined based on Equation (23) and the fluid film reaction forces are computed by numerically integrating the total pressure field over the journal surface:

$$\left[\begin{array}{c}{\overline{F}}_r\\ \overline{F_t}\end{array}\right]={\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=1}^M{\overline{P}}_{i,j}}}\left[\begin{array}{c}\cos ({\theta }_i)\\ \sin ({\theta }_i)\end{array}\right]\Delta \theta \Delta \xi .$$

(55)

4. Results and Discussion

This section represents the simulation results for the pressure distribution and the fluid film reaction forces for the proposed hydrodynamic SFD model. The simulation results are provided at different operating conditions, including eccentricity ratios and Reynolds numbers (i.e., inertia effects). The numerical algorithm that was developed in the previous section is incorporated into Matlab and Simulink^® to evaluate the effect of SFD operating parameters on the lubricant pressure distribution and the fluid film reaction forces in SFDs. In order to investigate the contribution of convective inertia on the SFD dynamics, the results of the simulations are compared against [24], where the effect fluid inertia is described for SFDs executing small-amplitude CCOs. According to an order of magnitude analysis, for small amplitude motions of the journal center, convective inertia terms in the flow equations are negligible relative to unsteady (temporal) inertia terms [6]. However, for large amplitude motions of the journal center, including the displacement at critical speeds, the effect of convective inertia is no longer negligible and should be necessarily incorporated into the calculations. Additionally, this work incorporates a π-film cavitation model to determine the fluid film reaction forces in the SFDs, however, in order to compare the pressure profiles at different SFD operating conditions, the negative intervals of the pressure calculations are included in the figures.

Furthermore, the fluid film reaction force components for the complete inertia model (i.e., both convective inertia and temporal inertia) and the temporal inertia model in [24] are compared with the force coefficient model developed by Vance [47]. Vance has developed the force coefficients for short-length open-ended SFDs by using the π-film assumption (i.e., film cavitation region is developed in half the damper circumference). The force coefficients are represented as follows:

$$C_{tt}=\frac{\pi \mu D}{4{\left(1-{\epsilon }^2\right)}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{\left(\frac{L}{c}\right)}^3,$$

(56)

$$C_{rt}=\frac{\mu \epsilon D}{{\left(1-{\epsilon }^2\right)}^2}{\left(\frac{L}{c}\right)}^3,$$

(57)

$$M_{rr}=\frac{\pi \rho D}{24}\left(\frac{L^3}{c}\right)\left[1-2{\left(1-{\epsilon }^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right]\left\{\frac{{\left(1-{\epsilon }^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-1}{{\epsilon }^2{\left(1-{\epsilon }^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}\right\},$$

(58)

$$M_{tr}=-\frac{27\rho D}{140\epsilon }\left(\frac{L^3}{c}\right)\left[2+\frac{1}{\epsilon }\ln \left(\frac{1-\epsilon }{1+\epsilon }\right)\right].$$

(59)

The fluid film reaction force components are calculated based on the force coefficients in Equations (56)–(59) as follows:

$$\begin{array}{l}{F}_r=-\left(C_{rt}{V}_t+M_{rr}{A}_r\right),\\ F_r=-\left(C_{tt}{V}_t-M_{tr}{A}_r\right),\end{array}$$

(60)

where

$$\begin{array}{l}{V}_t=\epsilon c\mathsf{\Omega },\\ A_r=-\epsilon c{\mathsf{\Omega }}^2,\end{array}$$

(61)

and

$$\begin{array}{l}\overline{F_r}=F_r/C_f,\\ \overline{F_t}=F_t/C_f,\end{array}$$

(62)

where

$$C_f=\mu \mathsf{\Omega }{R}^4/c^2.$$

(63)

Figure 3 represents the effect of fluid inertia on the dimensionless lubricant pressure distribution at the bearing midplane ($\xi =0$) at small eccentricity ratios for Reynolds equation (i.e., no inertia effect), the temporal inertia model in [24], and the proposed complete inertia model. In general, the fluid inertia effect causes a significant elevation in the pressure magnitude, a change in the shape of the pressure profile, and a phase shift of the pressure peak in the direction of the journal precision. At small Reynolds numbers, the influence of the viscous forces makes the pressure profile closer to a sinusoid, however, at moderate and large Reynolds numbers the pressure is in phase with the gap acceleration and transforms into a cosine wave shape. Furthermore, the temporal inertia model and the complete inertia model shows demonstrate a close agreement at small eccentricity ratios for a small to moderately large range of inertia effects.

Figure 4 illustrates the effect of fluid inertia on the dimensionless lubricant pressure distribution at the bearing midplane ($\xi =0$) at moderate eccentricity ratios. At small Reynolds numbers, the temporal inertia model and the complete inertia model demonstrate a very close agreement. However, by increasing the effect of fluid inertia (i.e., Reynolds number), the discrepancy between the results of the two models grows noticeably. At moderately small and moderate Reynolds numbers, the divergence between the results of the two models is more significant in terms of the pressure magnitude, while the shape and the phase of the profile is unaffected. However, at moderately large Reynolds numbers, the effect of convective inertia terms considerably influences the shape of the pressure profile.

Figure 5 demonstrates the effect of fluid inertia on the dimensionless lubricant pressure distribution at the bearing midplane ($\xi =0$) at large eccentricity ratios. For large eccentricity ratios of the journal center, the effect of convective inertia effects is very significant on the pressure magnitude even at moderately small Reynolds numbers. Further increasing the effect of fluid inertia results in the noticeable divergence of the two pressure profiles. The results imply that the effect of convective fluid inertia is very significant at large eccentricity ratios, including at the resonance frequencies of the rotor systems incorporating SFDs.

Figure 6 represents the effect of fluid inertia on the SFD fluid film reaction forces at small eccentricity ratios for four SFD models, namely no inertia effects, temporal inertia effects, complete inertia effects, and the force coefficient model [47]. The fluid forces are calculated by numerically integrating the lubricant hydrodynamic pressure distribution over the journal surface. For a π-film SFD model at small to moderate inertia effects, the contribution of the inertia forces to the radial force component is a positive value. This positive contribution is added to the negative viscous radial forces and reduces the magnitude of the radial force component relative to the inertialess model, which diminishes the likelihood of bi-stable rotor operation [16]. Furthermore, the magnitude of the inertialess radial forces is negative, meaning that the force is directed towards the center of the bearing (i.e., inwards). However, introducing the effect of fluid inertia initially changes the value of the forces to positive at small eccentricities, demonstrating the outward direction of the forces (i.e., outwards). Subsequently, at larger eccentricity ratios, the value of the radial forces switches back to negative, meaning that the force direction is once again towards the center of the bearing. Moreover, the contribution of the inertial forces to the tangential fluid film reaction forces is a negative value, and is added to the already negative purely viscous tangential forces, thus increasing the total magnitude of the tangential forces. Furthermore, the comparison between the temporal inertia model and the complete inertia model shows that at small Reynolds numbers the results of the two models are in a close agreement. At moderately small inertia effects, the discrepancy between the results of the two models is noticeable. At small eccentricity ratios, the radial forces suggested by the two models are in close agreement, however, increasing the eccentricity ratio leads to the divergence of the force calculations. The discrepancy between the force calculations is much more significant for the tangential force components, demonstrating the substantial influence of convective inertia effects on the SFD tangential forces. Further increasing the Reynolds numbers results in a more significant discrepancy between the radial force calculations at large eccentricity ratios and the tangential forces.

Additionally, the results of the complete inertia model and the force coefficient model are in close agreement for small to moderate Reynolds numbers. Increasing the Reynolds number leads to the divergence of the force calculations, since it is suggested that the force coefficient model is mostly accurate at small to moderate inertia effects. Furthermore, the force coefficient model is derived by incorporating the short bearing approximation, which reduces the accuracy of the calculations at moderate and large eccentricity ratios.

Subsequently, the proposed SFD model is incorporated into a multi-mass flexible rotor model to study the effect of SFD inertia forces on the steady-state mass unbalance induced vibrations and the transient orbits of the rotor system. The details of the proposed rotordynamics model is given in [48]. The rotor model is discretized by using finite element analysis (FEA) into Timoshenko beam elements [49] and the element mass matrix, stiffness matrix, and gyroscopic effects matrix are defined and assembled to provide a global rotor system model. Furthermore, the proposed nonlinear cavitated SFD model with centering springs, including lubricant inertia effects, is incorporated into the rotordynamic model. Finally, a predictor–corrector transient modal integration technique is implemented to determine the transient and the steady-state unbalance response of the rotor system at different operating conditions.

The main components of a conventional rotor system include the shaft, the disk, the mass unbalance, the rolling elements, the SFDs, and the retaining springs. This work applies finite element analysis [50] to determine the vibration amplitudes in the rotor system. In order to apply FEA to the rotor system, firstly, the structure is discretized into subdivisions of simple geometry that are referred to by elements. Subsequently, the elastic, inertia, damping, and external forces and moments on each element (local forces) are expressed in terms of the local coordinates (i.e., translations and rotations). Furthermore, the forces and the moments from every element are assembled together to produce the global forces in terms of the generalized coordinates. Finally, the discretized equation of motion is solved by using numerical techniques. In a rotor system, the shaft is discretized into several beam elements with nodes at both ends of each element. The disks and the bearings are assumed to be attached to the shaft at these nodes. In this work, the rotor vibrations in the lateral and transverse directions are considered. Consequently, each node has four generalized coordinates, including transverse displacements in the x and y directions as well as rotations about x-axis and y-axis. The properties of a rigid disk are determined by using Lagrange’s energy method. It is assumed that the disk strain energy is negligible with respect to the disk kinetic energy. Subsequently, the kinetic energy of the disk due to translation and rotation are calculated, and Lagrange’s energy method is applied to the energy equation to determine the element mass matrix and element gyroscopic matrix for the disk. The shaft contributes mass, stiffness, and gyroscopic effect to the rotor model. Similarly, the stiffness, inertia, and gyroscopic element matrices for the shaft are determined by using Lagrange’s energy method. The element matrices are derived for either the Bernoulli beam model or the Timoshenko beam model (includes both rotary inertia and shear effects). In order to determine the element matrices, firstly, the displacement along the shaft element is approximated by using shape functions. Subsequently, expressions for the kinetic energy and the strain energy are used to determine the element matrices. The strain energy of the shaft is calculated to approximate the element stiffness matrix. The kinetic energy due to the translation of the shaft gives the inertia element matrix and the kinetic energy due to the rotational motion of the shaft gives the gyroscopic effect matrix. The bearing element includes the stiffness and the damping contribution from the rolling elements as well as the stiffness contribution from the retaining springs in the SFDs. The bearings are flexible and absorb energy. It is assumed that the rotor is symmetric, and a disk is located at the axial center of the rotor. Furthermore, the rotor is supported with SFDs at either axial end. The SFDs are assembled on the outer race of two identical rolling elements that the rotor is assembled on. Moreover, the unbalance forces are present only on the centered disk. The equation of motion (EOM) of the rotor system is determined either by using the Lagrange’s energy method. For a conventional rotor system, the EOM is constructed as follows:

$$\left[M\right]\left\{\ddot{q}\right\}+\left(\left[C\right]+\mathsf{\Omega }\left[G\right]\right)\left\{\dot{q}\right\}+\left[K\right]\left\{q\right\}=\left\{F\left(t\right)\right\},$$

(64)

where q is a vector that contains the translations and the rotations at the nodes. The element matrices are 8-by-8 for stiffness, gyroscopic effect, and inertia. The effect of bearings, SFDs, and disks are directly incorporated into the corresponding nodes. Furthermore, the global mass, stiffness, and gyroscopic effect matrices are constructed by transforming the local element matrices into global coordinates and assembling the matrices together.

The transient response of the rotor system is determined to evaluate the effect of time variant nonlinear components (i.e., SFDs) on the performance of the rotor system. The most common and inclusive method for computing transient forced vibrations in rotor systems is the transient numerical integration approach. In the transient integration technique, the rotor system degrees of freedom are marched forward in time through a force balance at every time step. This technique allows any nonlinear elements (i.e., SFD forces) to be included directly as long as the forces can be represented as a function of position and velocity. Furthermore, transient modal integration is a powerful technique to determine the transient solution for dynamic systems with large number of degrees of freedom [51,52]. This technique is appealing since it facilitates the elimination of higher frequency modes with little impact on the system dynamics, which yields to significantly lower number of degrees of freedom compared to the direct integration method. Transformation of the rotor equation of motion into modal coordinates gives:

$$\left\{\ddot{p}\left(t\right)\right\}+\left[{\omega }^2\right]\left\{p\left(t\right)\right\}={\left[M_r\right]}^{-1}{\left[\overline{\varphi }\right]}^T\left\{F\left(t\right)\right\}-{\left[M_r\right]}^{-1}{\left[\overline{\varphi }\right]}^T\left(\left[C\right]+\mathsf{\Omega }\left[G\right]\right)\left[\overline{\varphi }\right]\left\{p\left(t\right)\right\},$$

(65)

where the model transformation is given by:

$$\left\{q\left(x,t\right)\right\}=\left[\overline{\varphi }\left(x\right)\right]\left\{p\left(t\right)\right\},$$

(66)

and the modal mass and the modal stiffness matrices are given as follows:

$$\left[M_r\right]={\left[\overline{\varphi }\right]}^T\left[M\right]\left[\overline{\varphi }\right],$$

(67)

$$\left[K_r\right]={\left[\overline{\varphi }\right]}^T\left[K\right]\left[\overline{\varphi }\right],$$

(68)

and F(t) is the sum of the external forces, including the unbalance forces and the SFD fluid film reaction forces. Assuming that the unbalance mass is attached to the disk, the unbalance force in the rotor is determined as:

$$F_{unb}\left(t\right)=\Re \left({\mathsf{\Omega }}^2{F}_0{e}^{j\mathsf{\Omega }t}\right),$$

(69)

where Ω²F₀ is the vector of forces and moments acting at the i^th node due to a disk with an unbalance, where F₀ is given as follows:

$$F_0=\left\{\begin{array}{c}{m}_{unb,i}\sigma e^{j\delta }\\ -j{m}_{unb,i}\sigma e^{j\delta }\\ j\left(I_{d,i}-I_{p,i}\right)\psi e^{j\gamma }\\ \left(I_{d,i}-I_{p,i}\right)\psi e^{j\gamma }\end{array}\right\}.$$

(70)

Assuming that no other unbalance forces are applied to the rotor system, the rest of the force vector entities are zero. Furthermore, the fluid film reaction forces are calculated by Equation (55). Subsequently, a transformation is applied to calculate the SFD forces in the fixed inertial coordinate (i.e., {X,Y}) as follows:

$$\begin{array}{l}{K}_{SFD,eq}=F_r/\left(\epsilon c\right),\\ C_{SFD,eq}=F_t/\left(\mathsf{\Omega }\epsilon c\right),\end{array}$$

(71)

and

$$\begin{array}{l}{F}_X=K_{SFD,eq}\cdot X+C_{SFD,eq}\cdot \dot{X},\\ F_Y=K_{SFD,eq}\cdot Y+C_{SFD,eq}\cdot \dot{Y}.\end{array}$$

(72)

Additionally, a numerical procedure integrating a predictor–corrector technique is incorporated into the simulation model to determine the rotor vibration amplitudes at different frequencies. In this study, the predictor uses backward finite difference method and the corrector uses Newmark–Beta method. The proposed transient modal integration algorithm is represented in Figure 7.

In the simulations, it is assumed that the unbalance forces are present only on the centered disk. Additionally, the rotor shaft is discretized into 20 elements.

Table 1. The rotor system parameters.

Parameter	Value	Unit
c11, c22	175.1268	N·s/m
c12, c21	0	N·s/m
$c/R$	0.00242
Dout, Din	0.0603 and 0.0413	m
Dd,out, Dd,in	0.2032 and 0.0603	m
D	0.1047	m
Ee	2.11 × 1011	pa
Gshear,e	8.1783 × 1010	pa
k11, k22	8.7563 × 106	N/m
k12, k21	0	N/m
Lshaft	0.5994	m
$L/D$	0.325
munbσ	3.6004 × 10−4, 1.4402 × 10−4, 7.2008 × 10−5	kg·m
Wd	0.0508	m
α, β	0.5 and 0.25
μ	0.0051 at 100 °C	Pa·s
ρe	7810	kg/m3
ρd	7810	kg/m3
ρ	1003.5	kg/m3
${\upsilon }_e$	0.29
Ω	100 to 15,000	rpm
Δt	10−5	s

summarizes the rotor and SFD parameters. A schematic diagram of the rotor in the simulations is demonstrated in Figure 8.

Figure 9 represents the effect of SFD viscous forces (i.e., neglecting inertia effect) on the steady-state vibration amplitudes of the rotor system at the disk node for the SFD radial clearance of 1.27 × 10⁻⁴, at the lubricant temperature of 100 °C, and for the three different unbalance magnitudes in Table 1. Firstly, the steady-state vibration amplitudes increase with the unbalance magnitude. Furthermore, even for the case where the effect of fluid inertia is neglected, the steady-state vibrations of the disk center is significantly attenuated at the first resonance frequency (i.e., the first bending mode) by incorporating the SFDs into the rotor system.

Figure 10 represents the effect of fluid inertia on the steady-state vibration amplitudes of the rotor system at the disk node for two different cases, namely: (1) Viscous SFD model (i.e., no inertia effects) and (2) complete fluid inertia effects. The results are represented at lubricant temperature of 100 °C and for the 3 unbalance magnitudes in Table 1. In general, the fluid inertia effects considerably attenuate the steady-state vibrations at the resonance zone and at large rotor velocities. At small to moderate rotor speeds, the squeeze Reynolds number is relatively small, and the effect of fluid inertia is negligible, consequently, the steady-state vibrations for the two cases is very similar. At the first resonance zone (i.e., first bending mode), the squeeze Reynolds number is moderately large, and the effect of fluid inertia is significant. Additionally, the vibration amplitudes at the resonance frequency zone are noticeably larger, resulting in larger SFD eccentricity ratios. Consequently, at the first resonance zone, the inertial correction component is large and the discrepancy between the results of the two models is noticeable, demonstrating the influence of fluid inertias effects on the suppression of the unbalance induced vibration amplitudes. At the post resonance zone, the rotor speeds are large, which results in the elevation of the Reynolds number and the inertia effects, hence, the inertia correction term remains significant and leads to the reduction of the steady-state vibration amplitudes relative to the no inertia model. Finally, the effect of fluid inertia leads to a shift in the first resonance frequency from 8250 rpm to approximately 8000 rpm. In general, the effect of SFD fluid inertia is represented by and added mass, which is supplemented into the overall rotordynamic system mass, which in turn reduces the first resonance frequency of the system.

Figure 11 represents the transient orbits of the SFD journal center at the 0.00242 radial clearance ratio and at 0.00036004 kg-m unbalance force at different rotor speeds for three different configurations: (1) No SFD; (2) Viscous SFD model (i.e., no inertia effects); (3) Complete inertia model. At small rotor speeds (1000 rpm and 5000 rpm), the squeeze Reynolds number is relatively small and the results for the three SFD configurations are in close agreement. Furthermore, including the SFD effect significantly accelerates the convergence of the orbits to the steady-state value. At the resonance zone (7310 rpm and 8250 rpm) the squeeze Reynolds number and the orbit radius are relatively larger, and the inertia models demonstrate both significantly faster convergence and smaller orbit radius. At larger rotor speeds (10,000 rpm and 12,500 rpm) the squeeze Reynolds number is moderately large while the orbit radius rapidly declines relative to the resonance zone.

5. Conclusions

This work represented expressions for the hydrodynamic pressure distribution and the lubricant velocity profiles, including fluid inertia effects, for finite-length open-ended squeeze film dampers incorporated into high-speed turbomachinery. The proposed pressure model provides very close agreement compared to the existing pressure models, including limiting geometry models. Moreover, a numerical procedure is proposed to solve the model for point-wise pressure distribution, which is significantly more computationally efficient relative to the bulk-flow technique. Finally, the numerical pressure field is numerically integrated over the journal surface to determine the fluid film reaction forces.

The analysis of the pressure model confirms the significant contribution of the fluid inertia to the dynamics of the lubricant. According to the results:

• Including the inertia effects increases the pressure magnitude, shifts the position of the pressure peak in the direction of the journal precision, and changes the shape of the pressure profile. The shape of the pressure is sinusoidal at small Reynolds numbers, while at moderate and large Reynolds numbers the pressure is in phase with the gap acceleration and transforms into a cosine wave shape.
• At small to moderate Reynolds numbers, the fluid inertia effect reduces the magnitude of the radial force component relative to the inertialess model, which diminishes the likelihood of bi-stable rotor operation. Furthermore, the direction of the radial forces initially changes to outwards at small eccentricity ratios, however, the direction switches back to inwards at moderate to large eccentricity ratios.
• The fluid inertia effect increases the magnitude of the tangential reaction forces; however, the direction of the forces remains unaffected.

Subsequently, the proposed SFD model was incorporated into a multi-mass flexible rotor model and the effect of SFD characteristics on the steady-state unbalance induced vibrations and the transient orbits of the rotor was investigated. According to the results of the simulations, the SFD lubricant inertia effects significantly attenuates the steady-state vibration amplitudes of the rotor system relative to the viscous SFD model (i.e., no inertia) at the resonance zone and at large rotor speeds. However, at small rotor speeds where the squeeze Reynolds number is small, the two models are in close agreement. Furthermore, the effect of lubricant inertia reduces the first resonance frequency of the rotor system. Increasing the fluid inertia effects further shifts the first resonance to lower frequencies.

The results of this work are especially valuable to the high-speed turbomachinery industry (i.e., jet engines and gas turbines) since it provides a method to quickly and accurately review damper designs and provide results as inputs to the system engineering team during the conceptual design phase of an engine/turbine. This could lead to significant development cost reduction through reduced system design iteration. Furthermore, the damper model is effectively integrated into the rotordynamic model of the complete system, providing a very powerful simulation tool to accurately predict the system vibrations during the development phase.