Identification of dynamic coefficients of gas foil bearings by simulated excitation based on coupled fields

Abstract

The dynamic coefficients (dynamic stiffness coefficient and dynamic damping coefficient) of gas foil bearings (GFBs) are important for transient dynamics calculation and stability analysis of the GFBs-rotor system. Although the perturbation method can be used to calculate the dynamic coefficients, it can only consider a single whirl ratio. Moreover, the whirl ratio is assumed and not practical. Therefore, this paper proposes an identification method of dynamic coefficients of GFBs by simulated excitation. Firstly, based on single journal-single GFB model considering gas-structure interaction by coupled fields, harmonic excitations were applied to the rotor to obtain the shaft orbit. Next, frequency response function method and equivalence coefficient method were used to identify dynamic coefficients. The former method is suitable for the case where response and excitation are same frequency, namely linear problem. The latter method is fit for the situation in which the response has multiple frequency components, namely nonlinear problem. Finally, examples were carried out, and the results from different methods were compared and analyzed. This research work lays a theoretical foundation for the design of the GFBs-rotor system.

Keywords

Gas foil bearings harmonic excitation whirl ratio dynamic coefficient simulated excitation

Introduction

GFBs are widely used in high efficiency blowers, aircraft air circulation cooling systems, micro fuel cell booster, etc., duo to its significant advantages such as high speed, oil-free, high temperature resistance, and low maintenance cost.^1,2 If the journal is motivated externally, it will whirl around the static equilibrium position. To describe the above phenomenon, the dynamic coefficients are introduced, which can characterize the relationship between the gas film force and rotor displacement under small-disturbance. Dynamic coefficients are very important parameters in the analysis of rotor dynamics, which can be used to determine natural frequencies of the bearing-rotor system or to predict transient responses. Hence, how to calculate the dynamic coefficients seems very critical.³

In previous attempts to solve the dynamic coefficients of GFBs, Constantinescu⁴ and Sternlicht et al.⁵ utilized an approximate method where they neglected the partial derivative of the gas film pressure with respect to time in the Reynolds equation. Ausman,⁶ Castelli and Elrod,⁷ and Czolczynski⁸ used the perturbation method to calculate these coefficients. However, they ignored the fact that the coefficients depend on the disturbance frequency. Guo⁹ coupled the compressible Reynolds equation to the deformation equation of the foil, to calculate the dynamic coefficients by using the perturbation method. Literature¹⁰ introduced complex number and partial differential method to solve the dynamic Reynolds equation numerically. In literature,¹¹ test fig having high-speed heavy-load radial GFB with a cylinder loading device was designed to test the bearing capacity of the bearing. In Chen and Lee,¹² according to measured rotor position and Fourier transform, an equation including the dynamic coefficients was established, and these coefficients were estimated using the least squares method. Literature¹³ developed a program for assessing the dynamic parameters of large bearings by applying dynamic loads using a vibrator. However, this program is only applicable to low-speed bearings. In Chen and Lee¹⁴ the dynamic characteristics of rolling elements in flexible rotor-bearing systems were evaluated by the finite element method (FEM). In literature,¹⁵ the bearing is assumed to be a linear model. Both volterra method and harmonic response probing technique were applied to estimate the dynamic coefficients by measuring the response amplitude with disregarding non-linear factors. In literature,¹⁶ based on the numerical model of the bearing, a neural network model was trained to predict the stiffness and damping coefficients. Literature¹⁷ identified the dynamic characteristic coefficients of bearings using the impulse response method. This method allows for the calculation of quality coefficients of bearings based on previous work, resulting in more accurate calculations. Literature¹⁸ proposes a new method for calculating LFCM that produces robust Cannon Bell diagrams and associated modal stability diagrams. The study concludes that predicting modal properties in different cases is closely related to direct linearization. In literature,¹⁹ the low dynamic viscous coefficient of air means that airfoil bearings offer limited damping at high speeds and may be susceptible to hydrodynamic instability. Therefore, it is necessary to perform dynamic design and analysis to ensure safe rotor operation. Obtaining the dynamic characteristics of the bearings is a prerequisite for dynamic rotor design. Literature²⁰ compared two methods for predicting the onset of instability and found a significant difference between the rotor instability speed predicted based on the bearing dynamic coefficients and the onset of instability speed obtained through time-domain nonlinear dynamics. The difference increases as the foil stiffness decreases. The literature²¹ theoretically concluded that the greater the flexibility of the bump structure, the greater the difference between the two methods. Literature²² presents a new method for calculating the dynamic characteristic coefficients of bearings. This method resolves the inconsistency between the instability onset velocity predicted by the force coefficient model and that obtained by using nonlinear time domain analysis. Literature²³ presents a new dynamics model that includes gas foil bearings, labyrinth seals, and impellers. The study investigates the effect of each component on the critical speed of the system. The results show that the rotor dynamics coefficients generated by the gas foil bearings are much larger than those of the other components. Therefore, the effect of the other components can be neglected when calculating the critical speed.

Parameter identification methods have shown certain advancements in solving dynamic coefficients, but there are still some issues. Firstly, in order to simplify the solution for the impact method, the main stiffness and main damping are assumed be the same, and the cross stiffness and cross damping are supposed to be opposite number. Secondly, the perturbation method can only consider a single whirl ratio (the ratio of the whirl speed to the rotating speed). Thirdly, the linear harmonic excitation method requires constructing two sets of excitations, leading to complex computations.

To solve the above challenges, this paper proposes an identification method for the dynamic coefficients of GFBs by simulated excitation. This method can search for the actual static equilibrium position automatically, and enables to consider multiple whirl ratios, that is, the nonlinear gas film force can be reflected in these coefficients. Furthermore, the method is simple and demonstrates good identification performance.

Model of GFBs-rotor system based on time-domain trajectory method

Basic equation of GFBs

The dimensional state equation for isothermal gas lubrication is as follows²⁴:

\begin{matrix} \frac{\partial}{\partial x} (p h^{3} \frac{\partial p}{\partial x}) + \frac{\partial}{\partial z} (p h^{3} \frac{\partial p}{\partial z}) \\ = 6 μ U \frac{\partial}{\partial x} (ph) + 12 μ \frac{\partial}{\partial t} (ph) \end{matrix}

(1)

For the universality and the stability of the numerical calculation, the dimensionalization of the equation is necessary. The non-dimensional forms of each parameter are as follows: $φ = x / R$ , $λ = z / R$ , $\bar{h} = h / C_{0}$ , $\bar{p} = p / p_{a}$ , $U = ω R$ , $τ = ω t$ . The dimensionless form of equation (1) is

\begin{matrix} \frac{\partial}{\partial φ} (\bar{p} {\bar{h}}^{3} \frac{\partial \bar{p}}{\partial φ}) + \frac{\partial}{\partial λ} (\bar{p} {\bar{h}}^{3} \frac{\partial \bar{p}}{\partial λ}) \\ = Λ \frac{\partial}{\partial φ} (\bar{p} \bar{h}) + 2 Λ \frac{\partial}{\partial τ} (\bar{p} \bar{h}) \end{matrix}

(2)

where, $Λ = 6 μ ω {(R / C_{0})}^{2} / P_{a}$ the bearing number, $R$ the journal radius, $ω$ the rotor angular speed, $C_{0}$ the radius gap, $P_{a}$ the atmospheric pressure.

Figure 1 shows a cross-section of a GFB, where, $ϕ$ the bearing expansion angle, and rotor eccentricity and attitude angle are denoted by $e$ and $θ$ , respectively. The expression for the film thickness, neglecting the higher-order terms, is as follows:

h = C_{0} + e \cos (ϕ - θ)

(3)

Figure 1.

Shaft position under arbitrary perturbation.

Numerical solution to transient film pressure

Equation (1) is an expression of the unsteady Reynolds equation. Let $S = {\bar{p}}^{2}$ , then equation (1) can be rewritten as:

\begin{matrix} \frac{\partial}{\partial φ} ({\bar{h}}^{3} \frac{\partial S}{\partial φ}) + \frac{\partial}{\partial λ_{c}} ({\bar{h}}^{3} \frac{\partial S}{\partial λ_{c}}) \\ = 2 Λ \frac{\partial}{\partial φ} (\bar{p} \bar{h}) + 2 σ \frac{\partial}{\partial τ} (\bar{p} \bar{h}) \end{matrix}

(4)

where, $σ = 2 Λ$ . The above equation is solved by the alternating implicit iteration method, in which space and the time are discretized using central difference format and backward difference format, respectively. The time step is divided into two segments $τ^{n} ~ τ^{n + 0.5}$ and $τ^{n + 0.5} ~ τ^{n + 1}$ .

Once pressures of all node are determined, the load capacity of the bearing can be obtained by

{\begin{matrix} F_{gx} = \int_{- L / (2 R)}^{L / (2 R)} \int_{0}^{2 π} (\bar{p} - 1) \sin φ d φ d λ_{c} \\ F_{gy} = \int_{- L / (2 R)}^{L / (2 R)} \int_{0}^{2 π} (\bar{p} - 1) \cos φ d φ d λ_{c} \end{matrix}

(5)

Rotor dynamics equation and numerical solution process

For the symmetric rigid rotor system supported by GFB as shown in Figure 2, the dynamic model at the static equilibrium position is expressed as follows:

{\begin{matrix} m \overset{\cdot\cdot}{x} + F_{gx} = me ω^{2} \cos (ω t) \\ m \overset{\cdot\cdot}{y} + F_{gy} = me ω^{2} \sin (ω t) \end{matrix}

(6)

where, $m$ the rotor mass. $F_{g}$ the gas film pressure, and subscripts x and y represent the directions of the pressure.

Figure 2.

GFBs-rigid rotor system model.

In fact, the general form of equation (6) is written as

ma (x_{c}, y_{c}, t) = f (x_{c}, y_{c}, t)

(7)

where, a represents acceleration, x_c, y_c are displacements of the journal from bearing center, and t denotes time.

Verlet integration algorithm is employed to construct the relationship between acceleration, velocity and displacement, as shown in equation (9), in which $Δ t$ is time step.

{\begin{matrix} {x_{c}}_{(n + 1)} = x_{cn} + v_{n} Δ t + 0.5 a_{n} (Δ t)^{2} \\ v_{n + 1} = v_{n} + 0.5 (a_{n} + a_{n + 1}) Δ t \end{matrix}

(8)

Also, we can derive a formula for the displacement y that similar to equation (8). Then, the eccentricity $e_{n}$ and attitude angle $θ_{n}$ of the rotor at each time step n, are written as:

e_{n} = \sqrt{{x_{cn}}^{2} + {y_{cn}}^{2}}, θ_{n} = \arctan \frac{x_{cn}}{y_{cn}}

(9)

Trajectory method based on coupled fields

In order to understand transient behavior of the journal, the trajectory method based on coupled fields is applied. It couples gas film pressure and rotor motion. The iterative process is depicted in Figure 3, and the specific steps are listed by

(1) Given initial conditions (eccentricity and attitude angle), solving equation (4) to obtain the initial gap film pressure. Integrate the pressure to obtain the gap film resultant force, and then derive the acceleration through equation (8).

(2) The displacement and velocity of the rotor at this time are calculated by using Verlet algorithm and equation (9).

(3) The eccentricity and attitude angle of the rotor are computed by equation (10), and as the initial value for the next iteration. Repeating the above steps will get information about the shaft orbit of the rotor with time-varying.

Figure 3.

The calculation process of trajectory method.

Calculation of dynamic coefficient

In fact, GFBs-rotor system has strong nonlinear characteristics. Usually, the higher the speed, the more significant the nonlinear phenomenon. For the case where rotor speed is small, the frequencies of input and output signals are same, which means a linear problem. Thereupon, the frequency response function method can be used to solve the dynamic coefficients. However, for the large rotor speed, the output signal will contain many components that are different from the excitation frequency. Thus, equivalence coefficient method was used to consider the nonlinear phenomenon.

Identification method of dynamic coefficients based on frequency response function

Frequency response function based on dynamic equation

The motion equation of the rotor under the excitation force is written by

{\begin{matrix} m \overset{\cdot\cdot}{x} + C_{xx} \overset{\cdot}{x} + C_{xy} \overset{\cdot}{y} + K_{xx} x + K_{xy} y = f_{x} (t) \\ m \overset{\cdot\cdot}{y} + C_{yx} \overset{\cdot}{x} + C_{yy} \overset{\cdot}{y} + K_{yx} x + K_{yy} y = f_{y} (t) \end{matrix}

(10)

where, x, y represent the displacements from static equilibrium point, $f_{x} (t), f_{y} (t)$ the external excitation in two directions, separately. $K_{xx}$ , $K_{yy}$ are the main stiffness coefficients, $K_{xy}, K_{yx}$ the cross stiffness coefficients, $C_{xx}, C_{yy}$ are the main damping coefficients, and $C_{xy}, C_{yx}$ are the cross damping coefficients.

Dimensionless form the equation (10) is

{\begin{matrix} \bar{x} ″ + \frac{β_{xx}}{w} \bar{x}' + \frac{β_{xy}}{w} \bar{y}' + γ_{xx} \bar{x} + γ_{xy} \bar{y} = H_{0} {\tilde{F}}_{x} (τ) \\ \bar{y} ″ + \frac{β_{yx}}{w} \bar{x}' + \frac{β_{yy}}{w} \bar{y}' + γ_{yx} \bar{x} + γ_{yy} \bar{y} = H_{0} {\tilde{F}}_{y} (τ) \end{matrix}

(11)

where, the dimensionless quantization is perform as follows, ${ω_{0}}^{2} = g / C_{0}$ the angular frequency $τ = ω_{0} t$ , $\bar{x} = x / C_{0}$ , $\bar{y} = y / C_{0}$ , $w = ω / ω_{0},$ $F_{stat} = mg$ , $H_{0} = C_{0} / F_{stat}$ . The dimensionless coefficients are $γ_{xx} = K_{xx} C_{0} / F_{stat}$ , $γ_{xy} = K_{xy} C_{0} / F_{stat}$ , $γ_{yx} = K_{yx} C_{0} / F_{stat}$ , $γ_{yy} = K_{yy} C_{0} / F_{stat}$ , $β_{xx} = C_{xx} C_{0} ω / F_{stat}$ , $β_{xy} = C_{xy} C_{0} ω / F_{stat}$ , $β_{yx} = C_{yx} C_{0} ω / F_{stat}$ , $β_{yy} = C_{yy} C_{0} ω / F_{stat}$ .

In the simulated excitation, harmonic excitation as input has the following expression

{\begin{matrix} {\tilde{F}}_{x} (τ) = F_{x} \sin (ω_{n} t) = F_{x} \sin (\frac{ω_{n}}{ω_{0}} ω_{0} t) \\ = F_{x} \sin (η τ) = IM (F_{x} e^{i η τ}) \\ {\tilde{F}}_{y} (τ) = F_{y} \sin (ω_{n} t) = F_{y} \sin (\frac{ω_{n}}{ω_{0}} ω_{0} t) \\ = F_{y} \sin (η τ) = IM (F_{y} e^{i η τ}) \end{matrix}

(12)

where, $F_{x}, F_{y}$ the amplitude of the external excitation, $ω_{n}$ the external excitation angular speed, $η = ω_{n} / ω_{0}$ the non-dimensional frequency.

The corresponding response, namely, rotor displacement as output is expressed as

{\begin{matrix} x (τ) = IM (x^{*} e^{i η τ}) = IM (x e^{i ε} e^{i η τ}) \\ y (τ) = IM (y^{*} e^{i η τ}) = IM (y e^{i ε} e^{i η τ}) \end{matrix}

(13)

where, $x^{*} = x e^{i ε}$ , $y^{*} = y e^{i ε}$ , $ε$ is the phase difference between the external excitation and displacement, and $IM$ is an operator that choose the imaginary part of a complex number.

Thereupon, the frequency response function of GFBs-rotor system is described as:

{\begin{matrix} H_{xx} = IM (\frac{x^{*} e^{i η τ}}{F_{x} e^{i η τ}}) = \frac{x}{F_{x}} e^{i ε_{xx}} = V_{xx} e^{i ε_{xx}} \\ H_{yy} = IM (\frac{y^{*} e^{i η τ}}{F_{y} e^{i η τ}}) = \frac{y}{F_{y}} e^{i ε_{yy}} = V_{yy} e^{i ε_{yy}} \\ H_{xy} = IM (\frac{x^{*} e^{i η τ}}{F_{y} e^{i η τ}}) = \frac{x}{F_{y}} e^{i ε_{xx}} = V_{xy} e^{i ε_{xy}} \\ H_{yx} = IM (\frac{y^{*} e^{i η τ}}{F_{x} e^{i η τ}}) = \frac{y}{F_{x}} e^{i ε_{yx}} = V_{yx} e^{i ε_{yx}} \end{matrix}

(14)

where, $H_{xx}$ represents the response in the x direction caused by the force in the same direction, and $H_{xy}$ is the response in the x direction caused by the force in the y direction. $H_{yy}, H_{yx}$ have the same meaning. $V_{xy}$ represents the amplitude ratio.

Frequency response function based on trajectory method

In the research, we used numerical experimentation to calculate the frequency response function. According to the coupled field model proposed previously, the displacement of the rotor can be simulated by trajectory method. And some mathematical tools such as spectral analysis and Fourier transform will be used to construct frequency response function.

The external excitations with an angular speed $ω_{1}$ are assumed as

{\begin{matrix} F_{xx} = A \sin (ω_{1} t) = IM (A e^{i ω_{1} t}) \\ F_{yy} = A \sin (ω_{1} t) = IM (A e^{i ω_{1} t}) \end{matrix}

(15)

Corresponding displacements are defined by

{\begin{matrix} x_{1} = a_{1} \sin (ω_{1} t + {λ_{x}}_{1}) = IM (a_{1} e^{i (ω_{1} t + {λ_{x}}_{1})}) \\ y_{1} = b_{1} \sin (ω_{1} t + {λ_{y}}_{1}) = IM (b_{1} e^{i (ω_{1} t + {λ_{y}}_{1})}) \\ x_{2} = a_{2} \sin (ω_{1} t + {λ_{x}}_{2}) = IM (a_{2} e^{i (ω_{1} t + {λ_{x}}_{2})}) \\ y_{2} = b_{2} \sin (ω_{1} t + {λ_{y}}_{2}) = IM (b_{2} e^{i (ω_{1} t + {λ_{y}}_{2})}) \end{matrix}

(16)

In the simulated excitation, the $F_{xx}$ in equation (15) is first imposed, and the corresponding response are $x_{1}$ and $y_{1}$ which shown in equation (16). Then, the $F_{yy}$ is applied, and x₂ and y₂ are obtained.

The frequency response functions of numerical experimentation are defined as the followings

{\begin{matrix} h_{xx} = v_{xx} e^{i {λ_{x}}_{1}} \\ h_{yy} = v_{yy} e^{i λ_{y 2}} \\ h_{xy} = v_{xy} e^{i {λ_{y}}_{1}} \\ h_{yx} = v_{yy} e^{i λ_{y 2}} \end{matrix}

(17)

where, $v_{xx} = \frac{x_{1}}{F_{xx}} = IM (\frac{a_{1} e^{i (ω_{1} t + {λ_{x}}_{1})}}{A e^{i ω_{1} t}}) = \frac{a_{1}}{A} e^{i {λ_{x}}_{1}}$ ,

$v_{yy} = \frac{b_{2}}{A}$ , $v_{xy} = \frac{a_{2}}{A}$ , $v_{xy} = \frac{b_{1}}{A}$ have same definition, and indicate amplitude-frequency characteristics, and $λ_{x 1}$ , $λ_{y 2}$ , etc., are phase-frequency characteristics.

Fast Fourier Transform (FFT) is utilized to analyze the spectrum of the displacement. In order to get analytical formula of the displacement, we rewrite the discrete displacement signals by a Fourier series as follow

x^{*} (t) = A_{0} + \sum_{n = 1}^{\infty} [A_{n} \cos (ω_{xn} t) + B_{n} \sin (ω_{xn} t)]

(18)

where, $A_{0}, A_{n}, B_{n}$ are arbitrary constants, and can be expressed by

{\begin{matrix} A_{0} = \frac{1}{T} \int_{0}^{T} x (t) dt \\ A_{n} = \frac{2}{T} \int_{0}^{T} x (t) \cos (ω_{xn} t) dt \\ B_{n} = \frac{2}{T} \int_{0}^{T} x (t) \sin (ω_{xn} t) dt \end{matrix}

(19)

here, $T$ is whirling period of the rotor.

Consequently, ${x^{*}}_{1} (t)$ can be reconstructed as

\begin{matrix} {x^{*}}_{1} (t) = A_{1} \cos (ω_{x 1} t) + B_{1} \sin (ω_{x 1} t) \\ = \sqrt{{A_{1}}^{2} + {B_{1}}^{2}} \sin (ω_{x 1} t + λ_{x 1}) \end{matrix}

(20)

where, $ω_{x 1} = ω_{1}$ , $λ_{x 1} = \arctan A_{1} / B_{1}$ .

Identification of dynamic coefficients

The identification of dynamic coefficients is actually a curve fitting process of two frequency response functions. The goal is that finding appropriate coefficients to get best consistency. We define an error function, and use the particle swarm optimization (PSO) algorithm to search the optimal solution (minimum value). Applying the external excitation in two directions in equation (11), there are

{\begin{matrix} \bar{x} ″ + \frac{β_{xx}}{w} \bar{x}' + \frac{β_{xy}}{w} \bar{y}' + γ_{xx} \bar{x} + γ_{xy} \bar{y} = H_{0} {\tilde{F}}_{x} (τ) \\ \bar{y} ″ + \frac{β_{yx}}{w} \bar{x}' + \frac{β_{yy}}{w} \bar{y}' + γ_{yx} \bar{x} + γ_{yy} \bar{y} = 0 \\ \bar{x} ″ + \frac{β_{xx}}{w} \bar{x}' + \frac{β_{xy}}{w} \bar{y}' + γ_{xx} \bar{x} + γ_{xy} \bar{y} = 0 \\ \bar{y} ″ + \frac{β_{yx}}{w} \bar{x}' + \frac{β_{yy}}{w} \bar{y}' + γ_{yx} \bar{x} + γ_{yy} \bar{y} = H_{0} {\tilde{F}}_{y} (τ) \end{matrix}

(21)

By substituting equation (21) into equations (12) and (13), the first two equations of equation (22) are divided by ${\tilde{F}}_{x} (τ)$ , then the latter two equations are divided by ${\tilde{F}}_{y} (τ)$ . Finally, according to equation (14), the frequency response functions will be represented as non-dimensional functions

{\begin{matrix} {\bar{H}}_{xx} = \frac{H_{xx}}{H_{0}} = {\bar{V}}_{xx} e^{i ε_{xx}} = (γ_{yy} - η^{2} + i \frac{β_{yy}}{w} η) / N (η) \\ {\bar{H}}_{yy} = \frac{H_{yy}}{H_{0}} = {\bar{V}}_{yy} e^{i ε_{yy}} = (γ_{xx} - η^{2} + i \frac{β_{xx}}{w} η) / N (η) \\ {\bar{H}}_{xy} = \frac{H_{xy}}{H_{0}} = {\bar{V}}_{xy} e^{i ε_{xy}} = - (γ_{xy} + i \frac{β_{xy}}{w} η) / N (η) \\ {\bar{H}}_{yx} = \frac{H_{yx}}{H_{0}} = {\bar{V}}_{yx} e^{i ε_{yx}} = - (γ_{yx} + i \frac{β_{yx}}{w} η) / N (η) \end{matrix}

(22)

where, $\begin{matrix} N (η) = (γ_{yy} - η^{2} + i \frac{β_{yy}}{w} η) (γ_{xx} - η^{2} + i \frac{β_{xx}}{w} η) \\ - (γ_{xy} + i \frac{β_{xy}}{w} η) (γ_{yx} + i \frac{β_{yx}}{w} η) \end{matrix}$

The error function G is defined by using (17) and (22)

G = \sum_{i = x, y} \sum_{k = x, y} | {\bar{H}}_{ik} - h_{ik} |

(23)

In the PSO method, the solution to the error function is represented by a group of particles ${K_{1} . . . K_{n} . . . K_{N}}$ in which each element meets the constraints.

K_{n} = {[{K^{n}}_{xx}, {K^{n}}_{xy}, {K^{n}}_{yx}, {C^{n}}_{xx}, {C^{n}}_{xy}, {C^{n}}_{yx}, {C^{n}}_{yy}]}^{T}

Each particle is successively substituted into the equation (24) to obtain the error value vector.

G = [G_{1}, G_{2}, \dots, G_{N}]

Finding the minimum value of the error function of N particles and compare it with the next iteration. The specific process is listed as

(1) Initialize the position and velocity of the particle swarm, generate initial solution set as the current population;

(2) Calculate the fitness (error function value) of each particle in the current population;

(3) The fitness of particles is compared with each other to find the smallest value as the optimal value, record it to prepare for the next iteration;

(4) The optimal value of step 3 is compared with the optimal value of the previous iteration record, and the minimum value is selected as the global optimal value;

(5) Update the velocities and positions of the particle and return to the step 2 until the accuracy requirement is achieved (Figure 4).

Figure 4.

PSO calculation flowchart.

Equivalence coefficient method

Under external excitation, the system response may contain multiple frequencies. The gas film force and displacement obtained through trajectory method also contain nonlinear factors. If the harmonic component in the response is too large, the equivalent coefficient method can be used to solve the dynamic characteristic coefficients. The coefficients are determined by using the data from the trajectory method, including displacement, velocity, gas foil force, and acceleration. Then the actual air foil force can be calculated by using equation (25). Although the coefficients are linear, the gas film forces calculated by equation (25) contain nonlinear factors. From this point of view, the equivalent coefficient method can consider nonlinear effects. Further, substituting the data within all-time steps in the trajectory method into equation (25) yields N equations containing eight unknowns. These equations form an overdetermined system of equations, and the eight dynamic characteristic coefficients are further identified by the nonlinear least squares method with damping. This method considers the influences of all harmonic frequencies on the coefficients, which means a nonlinear system can be equivalent to a linear system. The implementation of the equivalent coefficients method is similar to section “Frequency response function based on trajectory method,”, but only one-time excitation is required.

The motion of the rotor without mass eccentricity are written as

{\begin{matrix} \overset{\cdot\cdot}{x} (τ_{i}) = \frac{1}{M} (f_{x} (τ_{i}) - f_{gas_x} (τ_{i})) \\ \overset{\cdot\cdot}{y} (τ_{i}) = \frac{1}{M} (f_{y} (τ_{i}) - f_{gas_y} (τ_{i})) \end{matrix}

(24)

Where $τ_{i} = T / N$ represents the time step, T represents the single cycle of the rotor, and N represents the number of steps. $f_{gas_x}$ , $f_{gas_y}$ are equivalent gas film force represented by linearized stiffness and damping coefficients, and can be expressed as

{\begin{matrix} C_{xx} \overset{\cdot}{x} (τ_{i}) + C_{xy} \overset{\cdot}{y} (τ_{i}) + K_{xx} x (τ_{i}) + K_{xy} y (τ_{i}) = f_{gas_x} (τ_{i}) \\ C_{yx} \overset{\cdot}{x} (τ_{i}) + C_{yy} \overset{\cdot}{y} (τ_{i}) + K_{yx} x (τ_{i}) + K_{yy} y (τ_{i}) = f_{gas_y} (τ_{i}) \end{matrix}

(25)

In equation (25), $f_{gas_x}$ , $f_{gas_y}$ are nonlinear force from the trajectory method based on coupled fields. But the coefficients, such as $C_{xx}$ , $k_{xx}$ , and so on, are linearized. This means that the linear coefficient contains non-linear effects. How to solve these coefficients is very critical. By substituting the data from the trajectory method at each moment equation (25) yields an overdetermined system of equations can be obtained. This system contains eight unknowns, which are the dynamic characteristic coefficients. To solve this system, a nonlinear least squares algorithm with damping can be used. The following is the specific procedure:

The error between the equivalent gas film force and the gas film force $f_{gx} (τ_{i})$ , $f_{gy} (τ_{i})$ by trajectory method is defined as:

e_{i}^{j} (K^{j}, C^{j}) = C^{j} [\begin{matrix} \overset{\cdot}{x} (τ_{i}) \\ \overset{\cdot}{y} (τ_{i}) \end{matrix}] + K^{j} [\begin{matrix} x (τ_{i}) \\ y (τ_{i}) \end{matrix}] - [\begin{matrix} f_{gx} (τ_{i}) \\ f_{gy} (τ_{i}) \end{matrix}]

(26)

where $C = [\begin{matrix} C_{xx} (τ_{i}) & C_{xy} (τ_{i}) \\ C_{yx} (τ_{i}) & C_{yy} (τ_{i}) \end{matrix}]$ $K = [\begin{matrix} K_{xx} (τ_{i}) & K_{xy} (τ_{i}) \\ K_{yx} (τ_{i}) & K_{yy} (τ_{i}) \end{matrix}]$ ,

$e_{i}^{j}$ denotes the error after iterating $f_{gas_x} (τ_{i})$ , $f_{gas_y} (τ_{i})$ , and $f_{gx} (τ_{i})$ , $f_{gy} (τ_{i})$ for j steps, when the rotor is at time i. i = 1⋯N, representing the rotor operating time, j = 1,2,3⋯ representing the number of iterations.

The overall error for the system of equations is obtained by summing the errors for each equation:

E^{j} = \sum_{i = 1}^{N} (| e_{i}^{j} |)

(27)

where $E^{j}$ denotes the overall error produced after iterating over j steps at all-time steps.

By constantly changing the assignments of $K$ and $D$ , $E^{j}$ satisfies the predefined convergence conditions. The iterative calculation process is as follows:

(1) To begin the algorithm, assign the values $K_{1}$ , $D_{1}$ , and the damping coefficient $ς_{0}$ and truncation error $ε$ . Next, substitute the given equation into equation (26) and calculate the initial error $E^{1}$ .

(2) After completing k iterations, the current Jacobi matrix $J_{k}$ and the error $e_{1 . . . N}^{k} (K^{k}, C^{k})$ can be determined:

J_{k} = \frac{\partial [e_{1}^{k} \cdot \cdot \cdot e_{N}^{k}]}{\partial [K_{xx}, K_{xy}, K_{yx}, K_{yy}, D_{xx}, D_{xy}, D_{yx}, D_{yy}]}

(3) Solve the incremental equation $Δ_{k}$ after k iterations:

Δ^{k} = - (J_{k}^{T} J_{k} + ς_{k} I)^{- 1} J_{k}^{T} [e_{1}^{k} \cdot \cdot \cdot e_{N}^{k}]

where $I$ represents the unit matrix, $ς^{k}$ denotes the damping factor of the algorithm.

(4) If $E^{k} < E^{k - 1}$ and $E^{k} < ε$ , then stop the operation and return the result of the required iteration. If $E^{k} < E^{k - 1}$ but $E^{k} > ε$ , update the damping factor $ς^{k + 1} = ς^{k} / 10$ and the dynamic coefficients $K^{k + 1} = K^{k} + Δ^{k}$ , $D^{k + 1} = D^{k} + Δ^{k}$ and return to step (2). If $E^{k} > E^{k - 1}$ , update the damping coefficients $ς^{k + 1} = 10 ς^{k}$ and recalculate the dynamic coefficients $K^{k + 1} = K^{k}$ , $D^{k + 1} = D^{k}$ at this point in time, then return to step 2.

The specific flowchart is given in Figure 5. As follows:

Figure 5.

Calculation procedure diagram of algorithm.

Example and analysis

The parameters of GFBs-rotor system are shown in Table 1. Two cases with different rotor speed are considered.

Table 1.

The main parameters of the GFB.

Parameters/unit	Case 1	Case 2
Bearing length $b / mm$	60.0	60.0
Bearing radius $R_{b} / mm$	25.0	25.0
Radius gap $C_{0} / mm$	0.05	0.05
External excitation/N	$\begin{matrix} F_{xx} = 5 \sin (500 t) \\ F_{yy} = 5 \sin (500 t) \end{matrix}$	$\begin{matrix} f_{x} = 3 \sin (300 t) \\ f_{y} = 3 \sin (350 t) \end{matrix}$
Atmospheric pressure $P_{a} / Pa$	$1.0 \times 10^{5}$	$1.0 \times 10^{5}$
Aerodynamic viscosity $μ (Pa \cdot s)$	$1.8 \times 10^{- 5}$	$1.8 \times 10^{- 5}$
Speed $ω / (r \cdot min^{- 1})$	$3.0 \times 10^{4}$	$6.0 \times 10^{4}$
Bearing number	1.7	1.7

where the external excitations applied in the frequency response function are represented by $F_{xx}$ , $F_{yy}$ , while $f_{x}$ , $f_{y}$ represent the external excitations applied in the equivalent coefficient method.

Case 1: The results of linear system

In case 1, rotor speed is 30,000 rpm, and shaft orbits of numerical simulation are shown in Figure 6. For the coast down of the rotor, when discarding and considering external excitation (ignoring unbalanced mass force), the shaft orbits are shown in Figure 6(a) and (b), respectively. It shows that the shaft orbit in Figure 6(b) are larger than that of Figure 6(a) because of the centrifugal effect.

Figure 6.

(a) Shaft orbit without unbalanced mass force and external excitation. (b) The shaft orbit without unbalanced mass force and with external excitation with $ω_{1} = 500 Hz$ .

Displacements of the rotor relative to the static equilibrium point under harmonic excitation are illustrated in Figure 7. It shows that the displacements seem a harmonic vibration. The spectrum of displacements is shown in Figure 8. Among these frequency components, the amplitude of 500 Hz (equal to the excitation frequency) is maximum. Meanwhile, there are other components, such as 79.5 Hz, 155.883 Hz, 182.167 Hz, and their amplitudes are very small. In summary, the case 1 seems a linear problem.

Figure 7.

(a) Displacement comparison in x direction with external excitation. (b) Displacement comparison in y direction with external excitation.

Figure 8.

(a) x-direction displacement spectrum analysis. (b) y-direction displacement spectrum analysis.

Using equation (21), the 500 Hz component $x^{*} (t)$ in the original displacement signal is extracted, and its waveforms are shown in Figures 9 and 10. Through Figure 9, the phase difference $λ_{x 1}$ in Equation (21) can be conveniently determined. Figure 10 shows that the reconstructed displacement exhibits good consistency with the original displacement, although there is a slight deviation in the vicinity of the amplitude. The reason is due to other frequency components with low-amplitude. The above results verify the correctness of the theory in section “Frequency response function based on trajectory method” and lay a foundation for determining the frequency response function.

Figure 9.

(a) The displacement in the x direction and the phase difference of the external excitation. (b) The displacement in the y direction and the phase difference of the external excitation.

Figure 10.

(a) Comparison of the x-direction reconstruction displacement and the original displacement. (b) Comparison of the x-direction reconstruction displacement and the original displacement.

The dynamic coefficients identified by the frequency response function are compared with the perturbation method, as shown in Table 2. With the results of the perturbation method as a reference benchmark, the conclusion can be summarized as follows.

Table 2.

Dynamic coefficients of GFB with 30,000 rpm.

Method	Dynamic stiffness coefficients (/ $\times 1 0^{6} N \cdot m^{- 1}$ )				Dynamic damping coefficients (/ $\times 1 0^{2} N \cdot s \cdot m^{- 1}$ )
	$K_{xx}$	$K_{xy}$	$K_{yx}$	$K_{yy}$	$C_{xx}$	$C_{xy}$	$C_{yx}$	$C_{yy}$
Method of perturbation	5.1419	2.7162	2.4425	18.468	10.742	3.7776	−6.1390	14.641
Frequency response function	3.8090	1.5041	1.5036	15.707	24.114	3.0386	−2.7804	16.921

(1) The average main stiffness coefficients $K_{xx}$ and $K_{yy}$ identified by the frequency response function are 3.8090 × 10⁶ N/m and 1.5707 × 10⁷ N/m and the relative errors of the perturbation method are respectively 25.92% and 14.95%. 2) The relative errors of the average cross stiffness $K_{xy}, K_{yx}$ are 44.62% and 38.44%, respectively. 3) The relative errors of the average cross damping coefficients $C_{xy}$ , average main damping coefficients $C_{yy}$ are 19.56% and 15.57%, respectively. The other damping coefficients are in the same order of magnitude. 4) The comparison of these data can be shown that method is reasonable to identify the dynamic coefficients in the linear system.

The reasons for the difference between the above two methods (perturbation method, frequency response function method) are as follows. The perturbation method is based on solving the dynamic Reynolds equation, ignoring the higher order terms of the nonlinear gas film force. But frequency response function method applies the coupled fields, and can consider more practical gas film force and rotor motion. In addition, the perturbation method can only consider a single whirl ratio, which is man-made assumptions and lacks practicality. However, the frequency response function method considers the actual whirl frequency. The frequency response function is fitted by PSO algorithm, and the dynamic coefficient is finally identified.

Case 2: The results of nonlinear system

In case 2, rotor speed is set as 60,000 rpm, and shaft orbits of numerical simulation are shown in Figure 11.

For the coast down of the rotor, when discarding and considering external excitation (ignoring unbalanced mass force), the shaft orbits are shown in Figure 11(a) and (b), respectively. Figure 11(a) illustrates the rotor has an elliptic shaft orbit for the absence of external excitation. However, under the effect of external excitation, the rotor whirls irregularly near the static equilibrium position, as shown in Figure 11(b). Compared with Figures 6 and 11, the form of rotor whirling will change with the rotation speed. The curves of the steady-state displacement over time are shown in Figure 12, and it can be observed that the effect of the external excitation on the displacement is very significant. To investigate the effect, the spectrum of rotor displacement under external excitation is plotted, as shown in Figure 13. For the case 1, the external excitation frequency (500 Hz) is main frequency as shown in Figure 8. Nevertheless, as to the case 2, the whirl frequency (89 Hz) is the main frequency as shown in Figure 13. It can be concluded from the results that the amplitude of the whirl frequency component increases with the rise of the speed, which makes the system tend to unstability. Figure 14 shows that the reconstructed displacement exhibits good consistency with the original displacement.

Figure 11.

(a) Shaft orbit without unbalanced mass force and external excitation. (b) The shaft orbit without unbalanced mass force and with external excitation with 60,000 rpm.

Figure 12.

(a) The time domain diagram of x-displacement. (b) The time domain diagram of y-displacement.

Figure 13.

(a) Spectrum of x-direction displacement. (b) Spectrum of y-direction displacement.

Figure 14.

(a) Comparison of the x-direction reconstruction displacement and the original displacement. (b) Comparison of the y-direction reconstruction displacement and the original displacement.

The calculation results of the perturbation method and the equivalence coefficient method are shown in Table 3. With the results of the perturbation method as a reference benchmark, the conclusion can be summarized as follows.

Table 3.

Dynamic coefficients of GFB with 60,000 rpm.

Method	Dynamic stiffness coefficients (/ $\times 1 0^{6} N \cdot m^{- 1}$ )				Dynamic damping coefficients (/ $\times 1 0^{2} N \cdot s \cdot m^{- 1}$ )
	$K_{xx}$	$K_{xy}$	$K_{yx}$	$K_{yy}$	$C_{xx}$	$C_{xy}$	$C_{yx}$	$C_{yy}$
Method of perturbation	6.415	0.8599	5.29	8.565	2.148	8.722	−9.59	7.543
Equivalence coefficient method	6.11	1.53	1.92	15.48	8.726	12.250	13.430	14.926

The stiffness coefficients obtained by both methods, namely, the equivalence coefficient method and the perturbation method, are relatively close. With the results of the perturbation method as a reference benchmark, the difference between K_xx and K_yy are 4.7% and 16%, respectively. The damping coefficients in the x and y directions are quite different, and the former is about seven times and two times of the latter, respectively.

The equivalent coefficient method includes all whirl effects, in other words, the nonlinear force at a certain rotational speed is forcibly expressed as an equivalent by linear coefficients. In a sense, although the obtained dynamic stiffness and damping coefficients are equivalent linearization coefficients, they also can reflect nonlinearity.

Modal and stability analysis of rotor

The stability of a gas foil bearing rotor system depends on the dynamic characteristic coefficients. Assuming that the stiffness is much greater than the damping and that the damping is negligible, the dynamic equations of the rotor system supported by the gas bearing can be expressed as follows:

M \overset{. .}{x} + Kx = 0

(28)

The rotor undergoes simple harmonic motion, as described by the characteristic equation:

(- M λ^{2} + K) x = 0

(29)

If the mass matrix $M$ and stiffness matrix $K$ are known, the characteristic root of the system, represented by $λ$ , can be calculated using the provided formula. This will allow for the natural frequency of the system to be determined, which can be used to assess the system’s stability.

Taking 30,000 rpm as an example, the dynamic characteristic coefficients can be substituted into (29) to obtain eight characteristic roots:

\begin{matrix} \begin{matrix} λ_{1} = - 857.31 i & λ_{2} = 857.31 i \end{matrix} \\ \begin{matrix} λ_{3} = - 1067.39 i & λ_{4} = 1067.39 i \end{matrix} \\ \begin{matrix} λ_{5} = - 1672.54 i & λ_{6} = 1672.54 i \end{matrix} \\ \begin{matrix} λ_{7} = - 2114.55 i & λ_{8} = 2114.55 i \end{matrix} \end{matrix}

By substituting the characteristic roots obtained into equation (28), four feature vectors can be derived:

\begin{matrix} u_{1} = [\begin{matrix} 1 \\ 0.76 \\ - 0.76 \\ - 0.51 \end{matrix}], u_{2} = [\begin{matrix} 1 \\ 0.89 \\ 185.61 \\ 123.06 \end{matrix}] \\ u_{3} = [\begin{matrix} 1 \\ - 1.32 \\ - 0.64 \\ 0.90 \end{matrix}], u_{4} = [\begin{matrix} 1 \\ - 1.32 \\ 180.29 \\ - 245.71 \end{matrix}] \end{matrix}

The natural frequency calculated in finite elements is (Table 4):

Table 4.

Modal comparison between analytical method and finite element method.

Result	Analyze/Hz	Finite element/Hz	Relative error
Rigid rotor 1 order	136.44( $λ_{1}$ )	136.14	0.226%
Rigid rotor 2 order	169.88( $λ_{2}$ )	162.03	4.847%
Rigid rotor 3 order	266.19( $λ_{3}$ )	274.75	−3.114%
Rigid rotor 4 order	336.54( $λ_{4}$ )	326.19	3.174%
First-order bending	–	952.46
Two-order bending	–	2613.41

Figures 15(a) to (d) show the rigid body modes of the rotor, while Figures 15(e) and (f) show the bending modes. The rotor experiences parallel vortex (resonance) when the external excitation frequency reaches 136.44 Hz and 169.88 Hz, and conical vortex (resonance) when the external excitation frequency reaches 266.19 Hz and 336.54 Hz.

Figure 15.

Rotor mode at 30,000 rpm.

Conclusions

Based on the trajectory method, the fluid-solid coupling dynamic model of GFBs-rotor system is established, and the transient trajectory of the rotor is solved. The dynamic coefficients identification method for simulated excitation is proposed. The following conclusions are drawn:

(1) The average main stiffness coefficients $K_{xx}$ and $K_{yy}$ identified by the frequency response function are $3.8090 \times 1 0^{6} N / m$ and $1.5707 \times 1 0^{7} N / m$ , respectively, and the relative errors with the perturbation method are 25.92% and 14.95%. The relative errors of the average cross stiffness $K_{xy}, K_{yx}$ are 44.62% and 38.44%, respectively. The relative errors of the average cross damping coefficients $C_{xy}$ , average main damping coefficients $C_{yy}$ are 19.56% and 15.57%, respectively. The other damping coefficients are in the same order of magnitude. The comparison of these data can be shown that method is reasonable to identify the dynamic coefficients in the linear system.

(2) The stiffness coefficients obtained by both methods, namely, the equivalence coefficient method and the perturbation method, are relatively close. With the results of the perturbation method as a reference benchmark, the difference between K_xx and K_yy are 4.7% and 16%, respectively. The damping coefficients in the x and y directions are quite different, and the former is about seven times and two times of the latter, respectively.

(3) The equivalent coefficient method includes all whirl effects, in other words, the nonlinear force at a certain rotational speed is forcibly expressed as an equivalent by linear coefficients. In a sense, although the obtained dynamic stiffness and damping coefficients are equivalent linearization coefficients, they also can reflect nonlinearity.

The two methods proposed in this paper are to use the numerical simulation results of the coupled field to identify the dynamic parameters, so it is closer to the actual situation than the perturbation method.

Footnotes

Handling Editor: Sharmili Pandian

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by the National Natural Science Foundation of China (No. 51705413, 52275271, 11972282, 6110119004) and the Shaanxi Natural Science Foundation Research Plan (No. 2022JM-194).

ORCID iDs

Wenjie Cheng

Yuming Zhu

References

Zhao

. Theoretical and experimental research on new high damping gas foil bearings and their rotor systems. PhD Thesis, Hunan University, Hunan, China, 2017.

Zhang

. Dynamic analysis and experimental investigation on a new type of active bump-type foil bearing-rotor system. PhD Thesis, Hunan University, Hunan, China, 2018.

Wenjie

Xiao

, et al. Identification of dynamic characteristic coefficient of gas foil bearing based on the trajectory method. J Bearing 2022; 515: 4047.

Constantinescu

. Stability of the motion of circular bearings lubricated with gas. J Studii Si Cercetari Mecania Aplicata 1959: 10: 117–140.

Sternlicht

Poritsky

Arwas

. Dynamic stability aspects of cylindrical journal bearings using compressible and incompressible fluids. In: First international symposium on gas-lubricated bearings, Washington, DC, 1959, pp.119–160.

Ausman

. Linearized ph stability theory for translatory half speed whirl of long, self-acting gas-lubricated journal bearings. J Basic Eng 1963; 85: 611–619.

Castelli

Elrod

. Solution of the stability problem for 360 deg self-acting, gas-lubricated bearings. J Basic Eng 1965; 87: 199–212.

Czolczynski

. How to obtain stiffness and damping coefficients of gas bearings. J Wear 1996; 201: 265–275.

Guo

. Theoretical analysis and experimental investigation on nonlinear dynamic behaviors of rotor supported on gas foil bearings. PhD Thesis, Hunan University, Hunan, China, 2018.

10.

Lihua

. Performance study of tilting pad and elastic foil hydrodynamic gas bearings. PhD Thesis, Xi’an Jiao Tong University, Xi’an, China, 2009.

11.

Kai

Huang

, et al. Test of high speed and heavy load gas foil bearing and set up of test rig. J Hunan Univ Nat Sci 2020; 47: 1–5.

12.

Chen

Lee

. Estimation of linearized dynamic characteristics of bearings using synchronous response. Int J Mech Sci 1995: 37: 197–219.

13.

Morton

. Measurement of the dynamic characteristics of a large sleeve bearing. J Lubr Technol 1971; 93: 143–150.

14.

Chen

Lee

. Identification of linearized dynamic characteristics of rolling element bearings. J Vib Acoust 1997; 119: 60–69.

15.

Chatterjee

Vyas

. Nonlinear parameter estimation in rotor-bearing system using volterra series and method of harmonic probing. J Vib Acoust 2003; 125: 299–306.

16.

Khamari

Kumar

Behera

. Numerical investigation of influence sensitivity of a gas foil bearing parameters on the dynamic coefficients. J Braz Soc Mech Sci Eng 2021; 43: 1–19.

17.

Feng

Liu

, et al. Calculation of dynamic coefficients of air foil journal bearings using time-domain identification. Lubricants 2023; 11: 294.

18.

Breńkacz

. Identification of stiffness, damping and mass coefficients of rotor-bearing system using impulse response method. J Vibroeng 2015; 17: 2272–2282.

19.

Bonello

Pourashraf

. A comparison of modal analyses of foil-air bearing rotor systems using two alternative linearisation methods. Mech Syst Signal Process 2022; 170: 108714.

20.

Zhou

Ren

, et al. Modeling and stability characteristics of bump-type gas foil bearing rotor systems considering stick–slip friction. Int J Mech Sci 2022; 219: 107091.

21.

Larsen

Santos

von Osmanski

. Stability of rigid rotors supported by air foil bearings: comparison of two fundamental approaches. J Sound Vib 2016; 381: 179–191.

22.

Ren

Zhou

. A novel method for calculating the dynamic force coefficients of gas foil bearings and its application in the rotor dynamic analysis. J Sound Vib 2021; 515: 116466.

23.

Ying

Liu

Zhang

, et al. Impact of gas foil bearings, labyrinth seals, and impellers on the critical speed of centrifugal compressors for fuel cell vehicles: a comprehensive investigation. Lubricants 2023; 11: 532.

24.

Dimofte

. Fast methods to numerically integrate the Reynolds equation for gas fluid films. In: STLE-ASME joint tribology conference, St. Louis, MO, 1992.