The Influences of Parameters on the Dynamic Characteristics of a Multi-Foil Aerodynamic Journal Bearing with Bump-Backing Foils: Model Predictions

Abstract

In this work, the development and implementation of a dynamic characteristics model for a specific multi-foil aerodynamic journal bearing with bump-backing foils (MFJB) is considered. Based on the previously established static characteristics model, the elastohydrodynamic influence is carefully considered, and the perturbation method is adopted, as this model is more effective and computationally efficient. The effects of the operational, structural, and geometric parameters on stiffness and damping coefficients are emphasized. The results show that the eccentricity ratio effects are more intensive when the bearing speed is at a moderately high level, which is no more than approximately 30,000 rpm. The foil thickness has obvious effects on dynamic characteristics, whereas the influence of the elastic modulus is limited. Within the research scope, the eight-foils bearing exhibits a better performance than the four-foils. This paper is designed to provide effective methods and supply theoretical guidance for improving the engineering design and operational stability of bearings.

1. Introduction

As an essential type of mechanical supporting element, it is common for the foil aerodynamic journal bearing to be investigated and applied. It offers inherent advantages such as being self-acting, reliable, adaptable, and non-polluting [1,2]. Thus, it plays a crucial role in supporting rotating machinery, for instance, turbo-expanders [3,4], turbochargers [5,6], small micro-turbines [7,8], and turboshaft propulsion engines [9]. Consequently, the investigation of the foil aerodynamic journal bearing has attracted many scholars’ attention [10,11,12,13].

Generally speaking, the dynamic characteristics are one of the most dominant research domains of foil aerodynamic journal bearings [6,14,15]. The Finite Difference Method (FDM) and the small perturbation method were adopted by Li et al. [16] to simulate the influence of installation position on aerodynamic journal bearings’ dynamic performances. They found that a sensitive area may exist, positioned approximately complementary to the bearing attitude angle numerically, which could obviously affect the dynamic performances. This provides guidelines for the installation slot of the bearing. Then, Zhao et al. [17] used the Finite Element Method (FEM) and beam elements to stimulate the interactions between top and bump foils. The friction contacts were taken into consideration, and they suggested that this could be the cause of hysteretic behavior and energy dissipation of the bearings. Meanwhile, the deformation model of the bump foil, using the elastic theory, was established by Xu et al. [18] based on the compressible Reynolds equation. The influences of structural parameters of the bump foil on the dynamic characteristics were investigated, and they found that the increase in bump pitch could result in the growth of stiffness and damping coefficients. The aforementioned results supply valuable guidelines for the engineering design and practical applications of aerodynamic foil journal bearings. However, most of them are concerned about the comparatively simple Hydresil single-ring bump foil type of aerodynamic journal bearings (depicted in Figure 1a). Nevertheless, the corresponding explorations of the relatively more complicated multi-foil aerodynamic journal bearing with backing foils (MFJB) (exhibited in Figure 1b) require improvement.

However, with the introduction of multi-foil and bump-backing foil structures, elastic deformations of the foils in the bearing will be generated while operating at high speeds. In addition, the structure of the flow path in the very narrow and variable cross-section clearance will be quite complex [19,20,21]. Thus, more challenges will be faced in the prediction of the dynamic characteristics. According to the previous literature, various scholars attempted to solve these problems through different methods [22,23,24,25,26]. Some of them evaluated the advantages and disadvantages of various numerical models, as well as their calculation algorithms. Andres et al. [23,24] successively put forward two FEM models of the top foils of the bearings. The first one is the single two-dimensional shell model, and the second one is the one-dimensional beam model that contains the foil elastic deformation, deflection, and lubrication fluid dynamics. Comparatively speaking, the second one exhibited better conformity with the experimental data. It proved the availability and superiority of the beam model in the simulation of the bearing. Later, Leister et al. [25] proposed the runtime-efficient beam-based approach to simulate the foil structure and found that the generally applied simple elastic foundation model could lead to the overestimation of the bearing force. Because the compliance of the top foil could result in an adverse impact on the load capacity. Arghir et al. [26] introduced a structural model based on the contact mechanics with the gaps and friction being included. In that model, the top-foil elasticity and three different kinds of gaps are considered, which are those between the rotor-top foils, the top-bump foils, and the bump-foils sleeves. That model was thought to not only conform to reality but also control the machining errors of foils. Wu et al. [27] put forward a model for the dynamic characteristics of gas micro-bearings while considering the elastic deformation, as well as gaseous rarefaction effects. They claimed that the rarefaction effects might happen within ultra-thin lubrication gas film and supplied some information about the dynamic characteristics under effective viscosity. Larsen [28] emphasized the significance of considering both easy programming and simple use, which offer valuable suggestions for promoting the efficiency and effectiveness of bearing simulation models in the subsequent investigation. Later, Osmanski and Larsen et al. [29] put forward a fully-coupled model that contained the truss-based foil model, as well as incorporated the foil mass model and dynamic friction model. On account of directly measurable quantities and no “engineering assumptions”, this model can predict the natural frequencies and mode shapes. Nevertheless, it needed to be further developed in terms of capturing the unbalanced response. Gu et al. [30] introduced a novel, fully-coupled elastohydrodynamic model for aerodynamics based on the two-dimensional beam elements, which considered the separation and friction behaviors of the structures of the foils. It exhibited good accordance with engineering facts; however, the dynamic characteristics still need to be further explained. Recently, Zywica et al. [31] adopted static and kinetic friction coefficients into the calculation of damping coefficients and compared the numerical results with the experimental data, which exhibited a similar tendency. However, the dynamic properties under different bearing loads need to be further explored. Later, Zywica et al. [32] established a FEM numerical model for the dynamic characteristics of a nonlinear Hydresil-type foil bearing. In this model, the effect of the excitation force amplitude is emphasized, as well as the original introduction of the assembly preload. This contributed to the development of an in-depth investigation of the entire foil bearing.

Nevertheless, those studies mostly concentrated on the comparatively simple Hydresil-type ring foil aerodynamic journal bearings, and the relevant insights of the comparatively more complicated multi-foil aerodynamic journal bearings call for further development. Further, when dealing with the MFJB, to the best knowledge of the authors, most literature was reported on the investigation of their static characteristics. Previous investigation was mainly focused on the lubrication performance of gas under various parameters. There is still a lack of clear understanding of these relatively more complicated types of bearings. In addition, some previous dynamic models required a large amount of computation and possessed finite precision. Thus, the accuracy and computational efficiency of dynamic characteristics numerical models for this type of bearing need to be reinforced, as well as the parameter effects and influencing rules on the dynamic characteristics requiring further explanation.

In the current study, the dynamic characteristics of the MFJB are investigated. The numerical model is built based on the static characteristic model from the preceding literature [20] by the authors. The perturbation method is adopted, and the FDM is applied in the solving procedure. This model has the advantages of programming simplicity, easy use of the program, and satisfactory accuracy. The effects of the operating parameters, foil structure parameters, geometric parameters, and disturbance frequency on stiffness coefficients and damping coefficients are systematically and specifically studied. This paper aims to lay an effective and efficient theoretical foundation and research tools for identifying the dynamic characteristics of this specific type of bearings and improving their engineering design and performance optimization.

2. Numerical Models

The schematic diagram of an MFJB is illustrated in Figure 2. In the current investigation, the compressible Reynolds equation incorporated with the gas film thickness equation was applied to predict the distributions of lubricating gas film, and the foil deformation was considered via a curved beam model. Simultaneously, both the perturbation method and the FDM were adopted to solve the dynamic characteristics, with the stiffness coefficients and damping coefficients being included.

2.1. Reynolds Equation

In the current research, the lubricating gas is presumed to be the isothermal ideal gas. The dimensionless Reynolds equation could be adopted to monitor the pressure distributions, which is as follows [33]:

$${\displaystyle \frac{\partial }{\partial \theta }}\left(P{H}^3{\displaystyle \frac{\partial P}{\partial \theta }}\right)+{\left({\displaystyle \frac{2R}{L}}\right)}^2{\displaystyle \frac{\partial }{\partial Z}}\left(P{H}^3{\displaystyle \frac{\partial P}{\partial Z}}\right)=\Lambda {\displaystyle \frac{\partial \left(PH\right)}{\partial \theta }}+2\Lambda \gamma {\displaystyle \frac{\partial }{\partial {\tau }_s}}\left(PH\right)$$

(1)

where

$$\theta ={\displaystyle \frac{x}{R}}; Z={\displaystyle \frac{z}{0.5L}}; P={\displaystyle \frac{p}{pa}}; H={\displaystyle \frac{h}{C}};\Lambda ={\displaystyle \frac{6\mu \omega }{pa}}{\left({\displaystyle \frac{R}{C}}\right)}^2; \gamma ={\displaystyle \frac{{\omega }_s}{\omega }}; \overline{t}={\omega }_{s}t$$

2.2. Gas Film Thickness Equation

The dimensionless gas film thickness could be expressed as follows:

$$H={\displaystyle \raisebox{1ex}{$h$}\!\left/ \!\raisebox{-1ex}{$C$}\right.}=\left(h_0+u\right)/C=H_0+U$$

(2)

Utilizing the stiffness superposition [20], the relation between the overall stiffness of the top foil K_top and the deformation u under the action of force F can be obtained by the following equation:

$${\mathit{F}=\mathit{K}}_{top}\mathit{u}$$

(3)

Additionally, the overall stiffness matrix K_bump of the bump foils can be built via integration. By matrixing and combining to each relevant node of the K_top, the top-bump foil global stiffness matrix K_all can be built.

$${\mathit{K}}_{all}={\mathit{K}}_{top}+{\mathit{K}}_{bump}$$

(4)

On account of the restrained foil displacement at the fixed end, it can be expressed that [34]:

$${\mathit{K}}_{all}{\mathit{U}}_{all}={\mathit{F}}_{all}$$

(5)

where ${\mathit{F}}_{all}=\left\{0, \right.\left.F_{p1},0,F_{p2}, \dots \right\}$, F_pi (i = 1, 2, …), and it corresponds to the pressure of the gas film at each node. Additionally, ${\mathit{U}}_{all}=\left\{u_{fx1}, \right.\left.{ u}_{fx2},{ \theta }_{z1},{\theta }_{z2}, \dots \right\}$, u_fxi, u_fyi(i = 1, 2, …), which is related to the movement of each node. By extracting the lateral movement u_fxi of each node, the radial deformation of the foil u_f can be derived [20].

2.3. Perturbation of Pressure Governing Equation

When the shaft is rotating smoothly at the bearing speed of ω, the balance position is set as (x₀, y₀). Once the small perturbation has occurred near the balance position, the displacement, velocity, and acceleration are given as (Δx, Δy), ($∆\dot{x}$, $∆\dot{y}$), and ($∆\ddot{x}$, $∆\ddot{y}$), respectively. The dimensionless forms of the aforementioned parameters are expressed by the following forms:

$$\left.\begin{array}{c}∆X={\displaystyle \frac{∆x}{C}}=\left|∆X\right|e^{i\overline{t}}\\ ∆Y={\displaystyle \frac{∆y}{C}}=\left|∆Y\right|e^{i\overline{t}}\end{array}\right\}$$

(6)

By derivation of Equation (6), the velocity and acceleration can be transformed into the following forms:

$$\left\{\begin{array}{c}∆\dot{X}=\left|∆X\right|e^{i\overline{t}}i=∆X_i\\ ∆\dot{Y}=\left|∆Y\right|e^{i\overline{t}}i=∆Y_i\\ ∆\ddot{X}=\left|∆X\right|e^{i\overline{t}}i·i=-∆X\\ ∆\ddot{Y}=\left|∆Y\right|e^{i\overline{t}}i·i=-∆Y\end{array}\right.$$

(7)

If the dimensionless gas film thickness H, gas film pressure P, and foil deformation U are Taylor expanded at the equilibrium positions H₀, P_0, and U₀, then:

$$\left.\begin{array}{c}P=P_0+P_x∆X+P_x^{\prime }∆\dot{X}+P_y∆y+P_y^{\prime }∆\dot{Y}\\ H=H_0+H_x∆X+H_x^{\prime }∆\dot{X}+H_y∆y+H_y^{\prime }∆\dot{Y}\\ U=U_0+U_x∆X+U_x^{\prime }∆\dot{X}+U_y∆y+U_y^{\prime }∆\dot{Y}\end{array}\right\}$$

(8)

Substituting Equation (8) into the dimensionless Reynolds Equation (1) by omitting the higher derivative term, merging the similar term, and sorting out, the following equations can be derived:

$${\displaystyle \frac{\partial }{\partial \theta }}\left(P{H}^3{\displaystyle \frac{\partial P_0}{\partial \theta }}\right)+{\left({\displaystyle \frac{2R}{L}}\right)}^2{\displaystyle \frac{\partial }{\partial Z}}\left(P{H}^3{\displaystyle \frac{\partial P_0}{\partial Z}}\right)=\Lambda {\displaystyle \frac{\partial \left(P_0{H}_0\right)}{\partial \theta }}$$

(9)

$$\begin{gathered} {\displaystyle \frac{\partial }{\partial \theta }}\left(P_0{H}_0^3{\displaystyle \frac{\partial P_x}{\partial \theta }}+P_x{H}_0^3{\displaystyle \frac{\partial P_0}{\partial \theta }}+3P_0{H}_0^2{H}_x{\displaystyle \frac{\partial P_0}{\partial \theta }}\right)+{\left({\displaystyle \frac{2R}{L}}\right)}^2{\displaystyle \frac{\partial }{\partial Z}}\left(P_0{H}_0^3{\displaystyle \frac{\partial P_x}{\partial Z}}+P_x{H}_0^3{\displaystyle \frac{\partial P_0}{\partial Z}}+3P_0{H}_0^2{H}_x{\displaystyle \frac{\partial P_0}{\partial Z}}\right) \\ =\Lambda {\displaystyle \frac{\partial }{\partial \theta }}\left(P_0{H}_x+P_x{H}_0\right)-2\Lambda \gamma (P_0{H}_x^{\prime }+P_x^{\prime }{H}_0) \end{gathered}$$

(10)

$$\begin{gathered} {\displaystyle \frac{\partial }{\partial \theta }}\left(P_0{H}_0^3{\displaystyle \frac{\partial P_x^{\prime }}{\partial \theta }}+P_x^{\prime }{H}_0^3{\displaystyle \frac{\partial P_0}{\partial \theta }}+3P_0{H}_0^2{H}_x^{\prime }{\displaystyle \frac{\partial P_0}{\partial \theta }}\right)+{\left({\displaystyle \frac{2R}{L}}\right)}^2{\displaystyle \frac{\partial }{\partial Z}}\left(P_0{H}_0^3{\displaystyle \frac{\partial P_x^{\prime }}{\partial Z}}+P_x^{\prime }{H}_0^3{\displaystyle \frac{\partial P_0}{\partial Z}}+3P_0{H}_0^2{H}_x^{\prime }{\displaystyle \frac{\partial P_0}{\partial Z}}\right) \\ =\Lambda {\displaystyle \frac{\partial }{\partial \theta }}\left(P_0{H}_x^{\prime }+P_x^{\prime }{H}_0\right)+2\Lambda \gamma (P_0{H}_x+P_x{H}_0) \end{gathered}$$

(11)

$$\begin{gathered} {\displaystyle \frac{\partial }{\partial \theta }}\left(P_0{H}_0^3{\displaystyle \frac{\partial P_y}{\partial \theta }}+P_y{H}_0^3{\displaystyle \frac{\partial P_0}{\partial \theta }}+3P_0{H}_0^2{H}_y{\displaystyle \frac{\partial P_0}{\partial \theta }}\right)+{\left({\displaystyle \frac{2R}{L}}\right)}^2{\displaystyle \frac{\partial }{\partial Z}}\left(P_0{H}_0^3{\displaystyle \frac{\partial P_y}{\partial Z}}+P_y{H}_0^3{\displaystyle \frac{\partial P_0}{\partial Z}}+3P_0{H}_0^2{H}_y{\displaystyle \frac{\partial P_0}{\partial Z}}\right) \\ =\Lambda {\displaystyle \frac{\partial }{\partial \theta }}\left(P_0{H}_y+P_y{H}_0\right)-2\Lambda \gamma (P_0{H}_y^{\prime }+P_y^{\prime }{H}_0) \end{gathered}$$

(12)

$$\begin{gathered} {\displaystyle \frac{\partial }{\partial \theta }}\left(P_0{H}_0^3{\displaystyle \frac{\partial P_y^{\prime }}{\partial \theta }}+P_y^{\prime }{H}_0^3{\displaystyle \frac{\partial P_0}{\partial \theta }}+3P_0{H}_0^2{H}_y^{\prime }{\displaystyle \frac{\partial P_0}{\partial \theta }}\right)+{\left({\displaystyle \frac{2R}{L}}\right)}^2{\displaystyle \frac{\partial }{\partial Z}}\left(P_0{H}_0^3{\displaystyle \frac{\partial P_y^{\prime }}{\partial Z}}+P_y^{\prime }{H}_0^3{\displaystyle \frac{\partial P_0}{\partial Z}}+3P_0{H}_0^2{H}_y^{\prime }{\displaystyle \frac{\partial P_0}{\partial Z}}\right) \\ =\Lambda {\displaystyle \frac{\partial }{\partial \theta }}\left(P_0{H}_y^{\prime }+P_y^{\prime }{H}_0\right)+2\Lambda \gamma (P_0{H}_y+P_y{H}_0) \end{gathered}$$

(13)

where Equation (9) is the steady state pressure governing equation at the equilibrium position, Equations (10) and (11) are the perturbation pressure governing equations for the x direction, whereas Equations (12) and (13) are the those for the y direction.

2.4. Relationship Between Gas Film Thickness and Pressure

During the rotation process of the bearing system, if the shaft is perturbated, which can be expressed as (Δx, Δy), at the equilibrium position (x₀, y₀), the dimensionless expression of the film thickness Equation (2) can be transformed into [20]:

$$H=H_0+∆X\mathrm{s}\mathrm{i}\mathrm{n}\theta -∆Yco\mathrm{s}\theta +∆U$$

(14)

By Substituting the Taylor expansion of Equation (8), the above formula, Equation (14), can be simplified as follows:

$$\left\{\begin{array}{c}{H}_x=U_x+sin\theta \\ H_x^{\prime }=U_x^{\prime }\\ H_y=U_y-cos\theta \\ H_y^{\prime }=U_y^{\prime }\end{array}\right.$$

(15)

Currently, Equation (5) can be transferred into the following form:

$$\left[F\right]={\left[K\right]}_{top}\left[u\right]+{\left[D\right]}_{top}{\displaystyle \frac{\partial \left[u\right]}{\partial t}}$$

(16)

where ${\left[D\right]}_{top}$ is the foil-damping matrix [20]. Further, the dimensionless form of the foil deformation equation can be expressed in the following form:

$$\left[\overline{F}\right]={\left[\overline{K}\right]}_{top}\left[U\right]+{\left[\overline{D}\right]}_{top}{\displaystyle \frac{\partial \left[U\right]}{\partial \overline{t}}}$$

(17)

where ${\left[\overline{K}\right]}_{top}$ and ${\left[\overline{D}\right]}_{top}$ are the dimensionless stiffness and damping matrixes of the foils, respectively. $\left[\overline{F}\right]$ is the dimensionless gas film pressure matrix. It exists as the following:

$$\left\{\begin{array}{c}{\left[\overline{K}\right]}_{top}={\displaystyle \frac{{\left[K\right]}_{top}C}{p_a}}\\ {\left[\overline{D}\right]}_{top}={\displaystyle \frac{{\left[D\right]}_{top}C\omega }{p_a}}\\ \left[\overline{F}\right]={\displaystyle \frac{{\left[F\right]}_G}{p_a}}\end{array}\right.$$

(18)

By substituting the Taylor expansion Formula (8) into the foil force deformation Equation (16) and considering the relationship between the partial differential terms shown as Formula (15), the following formula can be obtained after combining and sorting out:

$$\left[P_x\right]={\left[\overline{K}\right]}_{top}\left[H_x\right]-{\left[\overline{K}\right]}_{top}\mathrm{s}\mathrm{i}\mathrm{n}\theta -{\left[\overline{D}\right]}_{top}\gamma \left[H_x^{\prime }\right]$$

(19)

$$\left[P_x^{\prime }\right]={\left[\overline{K}\right]}_{top}\left[H_x^{\prime }\right]-{\left[\overline{D}\right]}_{top}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\theta +{\left[\overline{D}\right]}_{top}\gamma \left[H_x\right]$$

(20)

$$\left[P_y\right]={\left[\overline{K}\right]}_{top}\left[H_y\right]-{\left[\overline{K}\right]}_{top}\mathrm{s}\mathrm{i}\mathrm{n}\theta -{\left[\overline{D}\right]}_{top}\gamma \left[H_y^{\prime }\right]$$

(21)

$$\left[P_y^{\prime }\right]={\left[\overline{K}\right]}_{top}\left[H_y^{\prime }\right]-{\left[\overline{D}\right]}_{top}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta +{\left[\overline{D}\right]}_{top}\gamma \left[H_y\right]$$

(22)

2.5. Solution Method and Boundary Conditions

On account of the solution process of the disturbance pressure in the y direction and x direction being similar, the y direction is taken as an example to clarify the solving procedure.

Defining the following:

$$\left\{\begin{array}{c}\overline{KK}={\left({\left[\overline{K}\right]}_{top}·{\left[\overline{K}\right]}_{top}+{\gamma }^2{\left[\overline{D}\right]}_{top}·{\left[\overline{D}\right]}_{top}\right)}^{-1}{\left[\overline{D}\right]}_{top}\\ \overline{CC}={\left({\left[\overline{K}\right]}_{top}·{\left[\overline{K}\right]}_{top}+{\gamma }^2{\left[\overline{D}\right]}_{top}·{\left[\overline{D}\right]}_{top}\right)}^{-1}{\left[\overline{D}\right]}_{top}\gamma \end{array}\right.$$

(23)

Then, the disturbance terms of the gas film thickness in Equations (21) and (22) can be expressed as:

$$\left\{\begin{array}{c}{H}_y=\overline{KK}{P}_y+\overline{CC}{P}_y^{\prime }+sin\theta \\ H_y^{\prime }=\overline{KK}{P}_y^{\prime }-\overline{CC}{P}_y\end{array}\right.$$

(24)

For the convenience of programming and usability, the central difference form of the FDM is introduced, and the terms of $P_y$ and $P_y^{\prime }$ can be expressed as:

$$\left\{\begin{array}{c}{\left({\displaystyle \frac{\partial P}{\partial \theta }}\right)}_{i, j}={\displaystyle \frac{P_{i+1, j}-P_{i-1, j}}{2∆\theta }}\\ {\left({\displaystyle \frac{\partial P}{\partial Z}}\right)}_{i, j}={\displaystyle \frac{P_{i, j+1}-P_{i, j-1}}{2∆Z}}\\ {\left({\displaystyle \frac{{\partial }^{2}P}{\partial {\theta }^2}}\right)}_{i, j}={\displaystyle \frac{P_{i+1, j}-2P_{i, j}+P_{i-1, j}}{{\left(∆\theta \right)}^2}}\\ {\left({\displaystyle \frac{{\partial }^{2}P}{\partial Z^2}}\right)}_{i, j}={\displaystyle \frac{P_{i, j+1}-2P_{i, j}+P_{i, j-1}}{{\left(∆Z\right)}^2}}\end{array}\right.$$

(25)

The boundary conditions for the dynamic solution can be shown as:

$$\left\{\begin{array}{c}{P}_x\left(1:n\theta +\mathrm{1,0}\right)=P_x\left(\mathrm{0,1}:nz+1\right)=P_x\left(n\theta ,0:nz+1\right)=0\\ P_x^{\prime }\left(1:n\theta +\mathrm{1,0}\right)=P_x^{\prime }\left(\mathrm{0,1}:nz+1\right)=P_x^{\prime }\left(n\theta ,0:nz+1\right)=0\\ P_y\left(1:n\theta +\mathrm{1,0}\right)=P_y\left(\mathrm{0,1}:nz+1\right)=P_y\left(n\theta ,0:nz+1\right)=0\\ P_y^{\prime }\left(1:n\theta +\mathrm{1,0}\right)=P_y^{\prime }\left(\mathrm{0,1}:nz+1\right)=P_y^{\prime }\left(n\theta ,0:nz+1\right)=0\end{array}\right.$$

(26)

Because the gas film pressure of the bearing clearance is distributed symmetrically along the axial midsection, for better calculation efficiency and economy, the half-length of bearing pressure distribution is considered. Then, the symmetric boundary condition exists as the following:

$$\left\{\begin{array}{c}{P}_x\left(1:n\theta ,nz-1\right)=P_x\left(0:n\theta ,nz+1\right)\\ P_x^{\prime }\left(1:n\theta ,nz-1\right)=P_x^{\prime }\left(0:n\theta ,nz+1\right)\\ P_y\left(1:n\theta ,nz-1\right)=P_y\left(0:n\theta ,nz+1\right)\\ P_y^{\prime }\left(1:n\theta ,nz-1\right)=P_y^{\prime }\left(0:n\theta ,nz+1\right)\\ {\displaystyle \frac{\partial P_{x\left(i,j\right)}}{\partial Z}}={\displaystyle \frac{\partial P_{x\left(i,j\right)}^{\prime }}{\partial Z}}={\displaystyle \frac{\partial P_{y\left(i,j\right)}}{\partial Z}}={\displaystyle \frac{\partial P_{y\left(i,j\right)}^{\prime }}{\partial Z}}=0 (i=0:n\theta ,j=nz)\end{array}\right.$$

(27)

Through processing, both the perturbation pressures in the x direction and the y direction can be obtained. Then, the stiffness coefficients and the damping coefficients of this type of bearing can be derived and expressed in the following forms:

$$\left[\begin{array}{c}{K}_{xx}\\ K_{yx}\end{array}\right]=\left[\begin{array}{c}{\displaystyle \frac{\partial F_x}{\partial x}}\\ {\displaystyle \frac{\partial F_y}{\partial x}}\end{array}\right]={\displaystyle \frac{p_{a}RL}{C}}\left[\begin{array}{c}\overline{K_{xx}}\\ \overline{K_{yx}}\end{array}\right]=-{\displaystyle \frac{p_{a}RL}{C}}{\int }_0^1{\int }_0^{2\pi }{P}_x\left[\begin{array}{c}sin\theta \\ -cos\theta \end{array}\right]d\theta dZ$$

(28)

$$\left[\begin{array}{c}{K}_{xy}\\ K_{yy}\end{array}\right]=\left[\begin{array}{c}{\displaystyle \frac{\partial F_x}{\partial y}}\\ {\displaystyle \frac{\partial F_y}{\partial y}}\end{array}\right]={\displaystyle \frac{p_{a}RL}{C}}\left[\begin{array}{c}\overline{K_{xy}}\\ \overline{K_{yy}}\end{array}\right]=-{\displaystyle \frac{p_{a}RL}{C}}{\int }_0^1{\int }_0^{2\pi }{P}_y\left[\begin{array}{c}sin\theta \\ -cos\theta \end{array}\right]d\theta dZ$$

(29)

$$\left[\begin{array}{c}{D}_{xx}\\ D_{yx}\end{array}\right]=\left[\begin{array}{c}{\displaystyle \frac{\partial F_x}{\partial \dot{x}}}\\ {\displaystyle \frac{\partial F_y}{\partial \dot{x}}}\end{array}\right]={\displaystyle \frac{p_{a}RL}{\tau C}}\left[\begin{array}{c}\overline{D_{xx}}\\ \overline{D_{yx}}\end{array}\right]=-{\displaystyle \frac{p_{a}RL}{\tau C}}{\int }_0^1{\int }_0^{2\pi }{P}_x^{\prime }\left[\begin{array}{c}sin\theta \\ -cos\theta \end{array}\right]d\theta dZ$$

(30)

$$\left[\begin{array}{c}{D}_{xy}\\ D_{yy}\end{array}\right]=\left[\begin{array}{c}{\displaystyle \frac{\partial F_x}{\partial \dot{y}}}\\ {\displaystyle \frac{\partial F_y}{\partial \dot{y}}}\end{array}\right]={\displaystyle \frac{p_{a}RL}{\tau C}}\left[\begin{array}{c}\overline{D_{xx}}\\ \overline{D_{yx}}\end{array}\right]=-{\displaystyle \frac{p_{a}RL}{\tau C}}{\int }_0^1{\int }_0^{2\pi }{P}_y^{\prime }\left[\begin{array}{c}sin\theta \\ -cos\theta \end{array}\right]d\theta dZ$$

(31)

where K_xx and K_yy represent the direct stiffness coefficients in the x direction and y direction, respectively, as well as K_xy and K_yx representing cross-coupled stiffness coefficients. Additionally, D_xx and D_yy represent direct damping coefficients in the x and y directions, respectively. Further, D_xy and D_yx represent the cross-coupled damping coefficients.

2.6. Solution Flow Chart of Dynamic Characteristic

The solving process of the stiffness and damping coefficients of the MFJB is illustrated in Figure 3. The mesh structure and computation domains are exhibited in Figure 4 [20]. Based on calculated static pressure and film thickness distributions, the governing Equations (10) and (11), as well as (19) and (20), are solved simultaneously via FDM, respectively. The disturbance pressure P_x and $P_x^{\prime }$ in the x direction can be obtained. Similarly, the governing Equations (12) and (13) can be solved, meanwhile, with the film thickness Equations (21) and (22). The disturbance pressure P_y and $P_y^{\prime }$ in the y direction can also be derived. Consequently, the stiffness coefficients and the damping coefficients of the MFJB are derived by integration of the corresponding formulas from Equations (28) to (31).

The convergence criteria can be shown as follows:

$$\left\{\begin{array}{c}\sqrt{\sum _{i=0}^{n\theta }\sum _{j=0}^{nZ}{\left({\displaystyle \frac{{\delta }_{i,j}}{P_{i,j}}}\right)}^2}\le {10}^{-6}\\ \sqrt{\sum _{i=0}^{n\theta }\sum _{j=0}^{nZ}{\left({\displaystyle \frac{{∆P}_{i,j}}{P_{i,j}}}\right)}^2}\le {10}^{-6}\end{array}\right.$$

(32)

2.7. Validation of Model

To ensure the accuracy and efficiency of the established models, they were selected to calculate the relationships between the bearing stiffness and damping coefficients, and bearing speed. In Figure 5, the simulated results are compared with the experimental results and simulated data in the literature [35]. The parameters and operating conditions are kept identical, which can be easily looked up in the literature [35]. Specifically, as depicted in Figure 5, the curves with small solid balls represent the calculated results by the current model. The curves with hollow squares and hollow circles each represent the simulated and experimental results in the literature [35], respectively. Notably, the predicted results by current models exhibit high consistency with both the experimental data and the simulated results in the literature [35]. Further, the maximum deviation between the experimental and simulated data is no more than 8.7%. This indicates that the current model has good reliability and accuracy.

To further prove the validity of the numerical model, the simulated results from the model in this current research are compared with those in the literature [36]. The relevant results are exhibited in Figure 6. In addition, the thin metal shims in the literature [36] are inserted under the bump foils, which are in contact with the bearing housing, thus creating the multiple-lobe clearance cross-section. The relative bearing parameters are kept the same as those listed in the literature [36]. As shown in Figure 6, the tendency of the numerical results from this model accords well with those in the literature [36]. Additionally, the maximum deviation between them is no more than 6.4%. This indicates that the current model possesses good reliability and accuracy. On account of the actual operational state, the model in the current research can be regarded as meeting the requirements of investigation under most operating conditions. In addition, the corresponding bearing parameters for the current research in the following sections are listed in

Table 1. Relevant bearing parameters for current research.

Parameter Name	Value
Shaft radius (Rj)	23.49 × 10−3 m
Bearing length (L)	30.0 × 10−3~70.0 × 10−3 m
Bearing radius (Rb)	26.25 × 10−3 m
Incircle radius (Ri)	23.5 × 10−3 m
Foil radius (Rf)	25.25 × 10−3 m
Foil thickness (t)	1.0 × 10−4~2.5 × 10−4 m
Foil elastic modulus (Eb)	2.00 × 1011~2.20 × 1011 Pa
Foil number (N)	4~8
Ambient pressure (pa)	1.01325 × 105 Pa
Foil Poisson’s ratio (υb)	0.3
Bearing speed (ω)	3.0 × 104~1.2 × 105 rpm
Bump foil thickness (tb)	1.016 × 10−4 m
Span of bump foil (s)	4.2 × 10−3 m
Half-length of bump foil (l)	1.75 × 10−3 m
Dynamic viscosity of gas (μ)	1.932 × 10−5 Pa·s
Disturbance frequency (fd)	50–1500 Hz
Eccentricity ratio (ε)	0.1~0.7

.

3. Results and Discussion

In this section, the impacts of the operational parameters (including eccentricity ratio, disturbance frequency, and bearing number), foil structural parameters (such as foil thickness, foil elastic modulus, and foil number), and bearing geometric parameters (for instance, the length-to-diameter ratio) on the dynamic characteristics, such as the stiffness coefficients and damping coefficients, are investigated. This section is designed to provide more comprehensive insights into the dynamic characteristics of MFJBs under varying operating and design parameters.

3.1. Eccentricity Ratio

The dynamic characteristics play a significant role in the operational stability of MFJBs. Therefore, as a significant operating parameter, the influences of the eccentricity ratios on the stiffness coefficients and damping coefficients were first simulated. As depicted in Figure 7, the eccentricity ratio increased from ε = 0.1 to ε = 0.7, while the other conditions were kept at a bearing speed of ω = 2.0 × 10⁴ rpm, a disturbance frequency of f_d = 1000 Hz, and a bearing number of Λ = 13.2. Specifically, as seen in Figure 7a, both the direct stiffness coefficients K_xx and K_yy increase, whereas the amplification of K_yy is larger and is accelerated when the eccentricity ratio surpasses 0.6. In contrast, the cross-coupled stiffness coefficients K_xy and K_yx are orders of magnitude, less and more stable, respectively. As depicted in Figure 7b, with the rise of the eccentricity ratio, the direct damping coefficient D_xx declines gradually, and the other coefficients remain steady overall.

Subsequently, keeping the same disturbance frequency, the bearing speed and bearing number were increased to ω = 3.0 × 10⁴ rpm and Λ = 19.8. The relevant results are illustrated in Figure 8. Compared with Figure 7a, the variation trends in the stiffness coefficients are similar; however, with the acceleration of bearing speed and bearing number, the distinction between the direct stiffness coefficients K_xx and K_yy diminishes. Meanwhile, the cross-coupled stiffness coefficients K_xy and K_yx approach zero, which stays the same as the traditional type of foil aerodynamic journal bearing. Although it is known that the bearing load capacity can be improved with the increase in bearing speed, the tendency of stiffness coefficient is not obvious. This was confirmed in the investigation results in the literature [37]. It can be explained that due to the elasticity of the supporting foils and the compressibility of the lubrication gas film, the further enhancement potential of the stiffness coefficient is limited. Hence, too large a bearing speed is not beneficial to the improvement of the bearing stiffness coefficient.

As depicted in Figure 9, the disturbance frequency is enlarged to 1500 Hz, with other parameters being kept the same as those of Figure 8. The overall trend of Figure 9 is similar to that of Figure 8. However, the direct stiffness coefficient of different directions in Figure 9 stays closer. As for the damping coefficient, the values in Figure 9 are lower than those in Figure 8. It could be suggested that under relatively higher bearing speed and eccentricity ratio, the hydrodynamic effect of the lubrication gas will be strengthened. Then, the contact area of adjacent foils will be larger. As a result, the stiffness will be increased, whereas the damping coefficient will be decreased.

3.2. Disturbance Frequency

In this section, the effects of the disturbance frequency on the dynamic characteristics are explored. While ensuring the same bearing capacity, the typical different combinations of bearing speed, eccentricity ratio, and bearing number were chosen. For instance, one condition was ω = 4.0 × 10⁴ rpm, ε = 0.7, and Λ = 26.4, and the other was ω = 6.0 × 10⁴ rpm, ε = 0.5, and Λ = 39.7. The relevant variation laws are depicted in Figure 10 and Figure 11, respectively. As illustrated in Figure 10a, with the disturbance frequency lifting from f_d = 50 Hz to f_d = 1500 Hz, the direct stiffness coefficient K_xx is increased steadily. Nevertheless, the K_yy of the y direction is firstly declined and then elevated slightly. In the meantime, the cross-coupled stiffness coefficients K_xy and K_yx originally increase, then become stable and tend to be close to zero. This differs from the variation in the stiffness coefficient with disturbance frequency in Hydresil type bearings, which is illustrated in Figure 1b [38]. More specifically, though, the original stiffness coefficient of the multi-foil type bearing, which is illustrated in Figure 1a, could be less than that of the Hydresil-type bearing. Nevertheless, the multi-foil type with backing foils possesses a wider, relatively constant stiffness property range, as well as an expected damping property within the high-frequency scope [37]. As seen from Figure 10b, when the disturbance frequency becomes larger, the direct damping coefficients D_xx and D_yy are decreased, whereas the former exhibited a bigger reducing range. One interesting thing is that, firstly, D_yy is larger than D_xx; however, as the disturbance frequency exceeds about f_d = 750 Hz, D_xx will surpass D_yy at an enlarging disparity. In contrast, the cross-coupled damping coefficients D_xy and D_yx are, at first, slightly increased and then stable.

Comparing Figure 11a with Figure 10a, some differences can be observed. First of all, the stiffness coefficient K_yy of the y direction rapidly decreases after the disturbance frequency surpasses f_d = 1200 Hz. Further, as seen from Figure 11b, the direct damping coefficient of the x direction D_xx is more stable, and the D_yy of the y direction shows a reduced decline amplitude than in Figure 11b. Interestingly, the cross-coupled damping coefficients D_xy and D_yx show different variation trends.

3.3. Foil Thickness

The variations in bearing dynamic characteristics under different foil thicknesses (from t = 1.0 × 10⁻⁴ m to t = 2.5 × 10⁻⁴ m) are illustrated in Figure 12. The other parameters are set as typical operation conditions, which are a bearing speed of ω = 8.0 × 10⁴ m, an eccentricity ratio of ε = 0.6, a bearing number of Λ = 52.9, and a disturbance frequency of f_d = 800 Hz. With the increase in foil thickness, the direct stiffness coefficients K_xx and K_yy, as well as the direct damping coefficients D_xx and D_yy, are gradually enlarged and then flattened out. As explained within the scope of this research, when the foil thickness increases, the minimum gas film thickness will decline. Then, the compression influence of the gas film is strengthened, as well as the maximum gas film pressure is increased. Thus, the direct stiffness coefficient appears to be uptrend. However, if the foil thickness is further increased, the strengthening effect will be limited, and the variation in stiffness coefficient will level off. Meanwhile, it also confirms that the smaller cross-stiffness coefficient contributes to the operational stability of the bearing [39,40].

Figure 13 depicts the variation in dynamic characteristics under different foil thicknesses at a relatively small eccentricity ratio of ε = 0.3, and the other parameters are the same as in Figure 12. As illustrated in Figure 13a, the direct stiffness coefficient at x direction K_xx will be enlarged at first with the increase in foil thickness and become stable at t = 2.0 × 10⁻⁴ m and above. As for the cross-coupled stiffness coefficients K_xy and K_yx, the maximum peak values appear at t = 2.0 × 10⁻⁴ m. Additionally, as seen from Figure 13b, the direct damping coefficient at y direction D_yy is prominently increased at t = 2.0 × 10⁻⁴ m.

Furthermore, to explore the variation law of dynamic characteristics with foil thickness at different bearing speeds, the typical conditions of ω = 4.0 × 10⁴ rpm, ε = 0.6, Λ = 52.9, and f_d = 800 Hz are chosen, and the relevant results are illustrated in Figure 14. As shown in Figure 14a, with the rise in foil thickness, the direct stiffness coefficients K_xx and K_yy will originally rapidly increase and subsequently become stable. This indicates that the bearing dynamic characteristics are obviously affected by the bearing speed. When under a relatively lower bearing speed, the variation in dynamic characteristics of foil thickness becomes smoother. As illustrated in Figure 14b, the direct damping coefficients D_xx and D_yy will originally increase and then gradually decrease. In contrast, the cross-coupled damping coefficients D_xy and D_yx are overall decreased when the foil thickness increases.

3.4. Foil Elastic Modulus

In this section, the effects of foil elastic modulus on the dynamic characteristics will be considered. The variations in stiffness coefficients and damping coefficients are depicted in Figure 15. While the elastic modulus is increased from E_b = 1.90 × 10¹¹ Pa to E_b = 2.25 × 10¹¹ Pa, and the other parameters are set as ω = 8.0 × 10⁴ rpm, ε = 0.6, Λ = 52.9, and f_d = 800 Hz. As seen from Figure 15a, the damping coefficients K_xx and K_yy are slightly increased. Meanwhile, the K_xy and K_yx are almost kept the same. Additionally, merely the direct damping coefficient at x direction D_xx is marginally enlarged after the elastic modulus surpasses E_b = 2.15 × 10¹¹ Pa, whereas the other damping coefficients exhibit no obvious variation. It can be explained that although the bearing capacity will be added to the growth of the elastic modulus, nevertheless, the effect is too limited when compared with other factors.

To explore the influence of foil elastic modulus under a larger disturbance frequency, the typical large disturbance frequency of f_d = 1500 Hz is selected, and other parameters are kept the same. The relevant results are shown in Figure 16. As shown in Figure 16a,b, the varying trends in the stiffness coefficients and damping coefficients under different foil elastic modulus are stable. Thus, it reflects the limited effects of foil elastic modulus on the dynamic characteristics.

3.5. Foil Number

This section will emphasize the influences of foil number on the dynamic characteristics of the bearing. As depicted in Figure 17, when under four-foils and eight-foils conditions, the eccentricity ratio is increased from ε = 0.15 to ε = 0.6, and the other parameters are kept as ω = 8.0 × 10⁴ rpm, Λ = 52.9, and f_d = 800 Hz. As illustrated in Figure 17a, the overall stiffness coefficients under eight foils are larger than that of four foils. The direct stiffness coefficient, K_xx, is higher than K_yy, with less variation in amplitude. This indicates that within the research scope, the bearing operational stability of eight foils is better than that of four foils. As seen from Figure 17b, a relatively low peak value of direct damping coefficient at y direction D_yy is exhibited at ε = 0.6. Additionally, apart from D_xy, the damping coefficients under eight foils are larger than those of four foils. This also proves that within the research scope, the bearing of eight foils obtains better operational stability.

The variations in dynamic characteristics under different bearing numbers at four foils and eight foils are exhibited in Figure 18, and the other parameters are kept as ω = 8.0 × 10⁴ rpm, ε = 0.6, and f_d = 800 Hz. As illustrated in Figure 18a, with the increase in bearing number, the direct stiffness coefficient at y direction K_yy under four foils is decreased, whereas the direct stiffness coefficient at x direction K_xx is steadily increased. On the contrary, the direct stiffness coefficients K_xx and K_yy under eight foils remain increasing. When the bearing number surpasses Λ = 32, the K_yy under eight foils is larger than that of four foils. It suggests that a relatively larger bearing number is beneficial to the bearing operational stability, and the eight-foils type is more suitable for larger bearing number conditions. As illustrated in Figure 18b, apart from the cross-coupled damping coefficient D_yx under four foils, which is larger than that of eight foils, the damping coefficients under eight foils are higher in total. This indicates that within the research scope, the bearing of eight foils obtains better operational stability.

Additionally, the variation in dynamic characteristics at different foil numbers at various disturbance frequencies is investigated. The disturbance frequency is increased from f_d = 450 Hz to f_d = 1500 Hz, and the other parameters remain as ω = 8.0 × 10⁴ rpm, ε = 0.6, and Λ = 52.9. As depicted in Figure 19a, the direct stiffness coefficients K_xx and K_yy of eight foils are approximately linearly increased with the lift of disturbance frequency. In contrast, merely the direct stiffness coefficient at the y direction K_yy of four foils exhibits a similar variation tendency. When under eight or four foils, both the cross-coupled damping coefficients D_xy and D_yx are slightly decreased. As seen in Figure 19b, the damping coefficients of eight foils are larger than those of four foils and are not affected by the disturbance frequency. It can be inferred that within the research scope, the bearing of eight foils obtains better dynamic characteristics and operational stability.

3.6. Length-to-Diameter Ratio

In the current section, the effects of the length-to-diameter ratio (L/D) will be focused on. Firstly, the L/D is increased from 0.43 to 1.69, and the other parameters remain as ω = 8.0 × 10⁴ rpm, ε = 0.6, f_d = 800 Hz, and Λ = 52.9. The relevant changing rules for stiffness and damping coefficients are shown in Figure 20. As seen from Figure 20a, the direct stiffness coefficients K_xx and K_yy steadily increase with the rise of L/D, and the locally small peak values appear at L/D = 1.48. Additionally, the cross-coupled stiffness coefficients K_xy and K_yx are slightly increased and decreased, respectively. As depicted in Figure 20b, the direct damping coefficients D_xx and D_yy are overall increased, whereas the locally small peak values are also distinct at L/D = 1.48. With the rise of L/D, the cross-coupled damping coefficients D_xy and D_yx are decreased and increased, respectively.

Figure 20 depicts the influence of L/D on dynamic characteristics of bearing under a relatively large disturbance frequency of f_d = 1500 Hz. As depicted in Figure 21a, the K_xx and K_yy are at first increased and then decreased when they surpass L/D = 1.48. Additionally, the distinction between cross-coupled stiffness coefficients K_xy and K_yx is exhibited when it surpasses L/D = 1.48, which is increased and decreased, respectively. As shown in Figure 21b, the overall variation trends are similar to Figure 20b, whereas the locally small peak values appeared in advance at L/D = 1.27.

To explore the effect of L/D on dynamic characteristics under different bearing speeds, the typical parameters of ω = 4.0 × 10⁴ rpm, ε = 0.6, f_d = 800 Hz, and Λ = 52.9 are selected. The L/D is increased from 0.43 to 1.69, and the relevant stiffness and damping coefficients are simulated and depicted in Figure 22. As illustrated in Figure 22a, a locally small value of direct stiffness coefficient at y direction K_yy has appeared at L/D = 1.28, which is less than that of direct stiffness coefficient at x direction K_xx. Additionally, the locally small peak values of cross-coupled stiffness coefficients K_xy and K_yx appeared at L/D = 1.48. In contrast, the variation law of the damping coefficient in Figure 22b is similar to that of Figure 20b, which is under a relatively larger bearing speed of ω = 8.0 × 10⁴ rpm.

Now, this section will emphasize the influence of L/D on dynamic characteristics under different eccentricity ratios. The L/D is increased from 0.43 to 1.69, while the parameters are kept as ω = 4.0 × 10⁴ rpm, ε = 0.3, f_d = 800 Hz, and Λ = 52.9, and the relevant results are exhibited in Figure 23. As illustrated in Figure 23a, the variation in the stiffness coefficient under a lower eccentricity ratio is more obvious than that of Figure 22a, which is under a higher eccentricity ratio. For instance, the increasing tendency in direct stiffness for coefficient a y direction K_yy slows down. As depicted in Figure 23b, compared with Figure 22b, the variation in the damping coefficient is smoother, and the locally small value is not depicted. It shows that the effects of eccentricity ratio are more obvious than that of bearing speed.

4. Conclusions

In the current investigation, a resultful numerical model with high computational efficiency for the dynamic characteristics of the MFJB has been built. The one-dimensional curved beam model is introduced to describe the aeroelasticity driving foil deformation. Additionally, the compressible Reynolds equation coupled with the modified gas film thickness equation is solved by the adoption of the perturbation method. The significant influencing laws of parameters of operation, structure, and geometry on dynamic characteristics are comprehensively investigated and discussed. The main conclusions that can be drawn are as follows:

(1)

The established numerical model enables dynamic characteristics analysis of the MFJB. The effects of the parameters, especially, can be directly and accurately obtained by numerical results with high computational efficiency.

(2)

When the bearing speed is less than 30,000 rpm, the stiffness is strengthened with the increase of the eccentricity ratio at a declining magnitude. However, when the bearing speed reaches as high as 80,000 rpm, that enhancement effect created by the eccentricity ratio is confined. Similarly, too large of a foil thickness can restrict the strengthening effect of stiffness.

(3)

Compared with bearing speed, the bearing eccentricity ratio more effectively affects the dynamic characteristics when under different L/D. The foil thickness will have an obvious influence on the dynamic characteristics, whereas the influence of the elastic modulus is very limited.

(4)

Within the research scope, compared with the four-foil type, the eight-foil bearing exhibits overall better dynamic characteristics, which could result in its better operational stability in engineering applications.

In summary, these findings emphasize the parameter effects on the dynamic characteristics. Furthermore, researching the influence of lubrication temperature performance could help acquire insight into its impact to achieve better stiffness and damping performance.