Nonlinear dynamic response and bifurcation of asymmetric rotor system supported in axial-grooved gas-lubricated bearings

Abstract

Dynamic characteristics of the asymmetric rotor system supported in axial-grooved gas-lubricated bearings are studied. In order to solve nonlinear dynamic response of rotor system effectively, a hybrid numerical model is established by coupling the motion equation of rotor with the rational function model of the gas film forces. The rational function model of the gas film forces of gas-lubricated bearing is established based on vector fitting theory. By using the hybrid numerical model, the repeated calculations of the unsteady Reynolds equation and gas film forces are avoided; the continuous rotor trajectory and the dynamic gas film forces can be calculated simultaneously; and for the rotor system supported in the same bearings, the computing cost can be saved effectively. The nonlinear dynamic responses of asymmetric rotor system supported in axial-grooved gas-lubricated bearings are investigated by trajectory diagrams, frequency spectrum, Poincaré maps, and time series. The bifurcations are analyzed by the bifurcation diagrams with different rotating speeds and mass eccentricities. The dynamic behaviors of the asymmetric rotor system appear complex nonlinear dynamic phenomenon and specific bifurcation characteristics.

Introduction

Gas-lubricated bearing has many outstanding advantages in the rotating machinery, such as high speed, high precision, low friction loss, low noise, and so on. The rotor system supported in gas-lubricated bearing has better performance. There are many studies about the stability and dynamic characteristics of the gas-lubricated bearing–rotor system.^1-7 Wang et al. ⁸ analyzed the bifurcation and nonlinear dynamic behaviors of a flexible rotor supported by relative short gas film bearing, and the dynamic behaviors of the rotor system in the horizontal and vertical directions were investigated by the state trajectory, Poincaré maps, power spectra, and bifurcation diagrams. Li et al. ⁹ investigated the influence of the surface waviness on the nonlinear dynamic performance of a gas bearing–rotor system, and the stability of the rotor system was analyzed with different directions, the amplitudes, and the numbers of waves. Li et al. ¹⁰ established a dynamic model of the rotor–bearing system with bolted-disk joint by the proposed joint element and lumped mass modeling method, and the dynamic vibration characteristics and stability of rotor under different tangential stiffness and transition point were investigated. The gas-lubricated bearings above mentioned are cylindrical bearings. The rotor system supported in the gas-lubricated cylindrical bearing usually generates instability at higher speed.

In order to ensure stable operation and obtain higher performance, the grooved gas-lubricated bearings are employed. We ¹¹ studied the unbalanced dynamic responses of the rotor system supported by three axial-grooved finitely long gas bearing. A time-dependent mathematical model was developed and solved by the differential transformation method. The influence of the nonlinear gas film forces on the stability of rotor system was investigated via analyzing the bifurcation and chaotic behaviors. We ¹² proposed a calculation method coupling with the precise integration method for analyzing nonlinear dynamic behaviors of the self-acting axial groove gas bearing–rotor system with double time delays. The stability of the system can be improved via choosing the proper feedback control gains and the time delays. The whirl stability of the herringbone-grooved gas journal bearing was investigated by Fujita. ¹³ The critical stability of the bearing–rotor system with the change of bearing numbers was analyzed, which was determined by the rotor mass. The effects of the length-to-diameter ratio and the parameters of the groove on the whirl of a herringbone-grooved gas bearing were studied. Du et al. ¹⁴ analyzed the complicated nonlinear dynamic behaviors of the rotor–bearing system with a spiral-grooved opposed-hemisphere gas bearing support. The nonsynchronous excitation responses were analyzed differently from the way that the synchronous excitation responses were done. Han et al. ¹⁵ proposed a new parallel elastic hydrodynamic lubrication numerical algorithm. The lubrication performance of the misaligned herringbone-grooved axial bearing was investigated, and it can be improved via designing appropriate bottom shapes of herringbone grooves. Zhang et al. ¹⁶ investigated the effect of air rarefaction on air bearing forces and the contact behaviors of the air lubricated spiral-groove thrust micro-bearing. The bearing forces decreased and the asperity contact forces increase significantly by taking the air rarefaction into consideration. We ¹⁷ established a lubrication model of a micro-grooved three-pad fixing pad aerodynamic journal bearing. The effects of the parameters of the grooves in the bearing pad surfaces on the load performance of the three-pad aerodynamic bearing were investigated. By two-dimensional narrow groove theory, a model of a rigid rotor with herringbone-grooved journal gas bearing support was established by Liu et al. ¹⁸ The onset speed of the sub-synchronous vibration of the rotor system was predicted by the theoretical model and experimental data. In order to reduce the gas film vortex and oscillation and improve the operation stability of rotor system supported in spherical spiral groove hybrid gas bearing, Jia et al. ¹⁹ investigated the dynamic change of the gas film stiffness and damping coefficients under different motion conditions. The above studies focus on the dynamic characteristics of gas bearing–rotor system.

In order to calculate the dynamic response of bearing–rotor system more accurately and quickly, many novel models and methods were developed by scholars in the studies of the rotor dynamics. González et al. ²⁰ presented a methodology to identify nonlinear instabilities of gas bearing–rotor system. The wavelets and phase diagram were obtained to estimate the nonlinear parameters and evaluate the dynamic behaviors. Hassini and Arghir ²¹ proposed a simplified method to investigate the nonlinear transient behaviors of a flexible rotor with the gas bearing support. The large nonlinear displacements, stability diagrams, and Campbell diagrams were obtained by their proposed method which approximates the impedances using rational functions. Bonello ²² developed a linearized method of the foil air bearing symmetric rigid rotor system based on the system Jacobi, and the effects of air film constraints and top foil detachment on the bearing–rotor system were investigated by using the proposed method. By considering the micro-deformation and interaction of various foils, Zhou et al. ²³ presented a semi-directly experimental method of the bump-type foil journal bearing. The dynamic behaviors of the rotor–bearing system were predicted by analyzing the journal orbits. Baum et al. ²⁴ established a novel physical and mathematical model to investigate the rotor dynamic characteristics of a rigid rotor with the foil air bearing support; the dimension of the model was reduced by Galerkin’s method, and the computational cost of the dynamic characteristics can be reduced by the proposed model.

Many different analysis methods of the nonlinear stability have been developed. Yang et al. ²⁵ investigated the unbalance response of micro-gas bearing–rotor system; the results showed that stability threshold speed of micro-rotor system increased by considering the gas rarefaction effect. He et al. ²⁶ adopted a new technique of the nonlinear expanded frequency to attain the stability criteria; the nonlinear stability of the Kelvin Helmholtz instability saturated in porous media with heat and mass transfer was investigated. Based on the homotopy perturbation method, He et al. ²⁷ proposed a reducing rank approach to solve the strongly damping nonlinear Klein–Gordon equation; the stability conditions were established for the first time; and the stability analysis was discussed. Sun et al. ²⁸ investigated the stability of a magnetic bearing supported rotor system with significant gyroscopic effects by using the double-frequency Bode plot and the Nyquist criterion. Guo et al. ²⁹ proposed a dynamic model of a hydrodynamic floating ring bearing to study the influence of lubricant temperature–viscosity on the performance, and the stability criterion was obtained by the Routh–Hurwitz method.

In this paper, a hybrid numerical model of asymmetric rotor system supported in two axial-grooved gas-lubricated bearings is established to calculate the nonlinear dynamic responses and bifurcation of gas-lubricated bearing–rotor system efficiently. The rational function model of gas-lubricated bearing is established by the rational function approximation theory and vector fitting method. The hybrid numerical model is coupling of the rational function model with motion equation of rotor. The nonlinear dynamic responses of gas-lubricated bearing–rotor system are obtained by solving the hybrid numerical model, and the bifurcations are calculated under the conditions of different rotating speeds and mass eccentricities.

Model and Method

Model of gas-lubricated bearing

The schematic diagram of the gas-lubricated bearing is shown in Figure 1. O _b and O _j are the bearing and journal centers, respectively, $e\hspace{0.17em}\left(e=\left|\overrightarrow{O_b{O}_j}\right|\right)$ is the eccentricity; θ is the attitude angle; R is the radius of the journal; W is the external force; ϕ and φ are the angles calculated from direction y and O _b O _j to the current position; h is the gas film thickness; ω is the rotating speed of the rotor; c is the radius clearance between the bearing and rotor; and γ and γ₀ are the arc bushing angle and groove width angle.

Figure 1.

Schematic diagram of two axial-grooved gas-lubricated bearings.

For gas-lubricated bearing, the nondimensional steady compressible Reynolds equation is

$$\frac{\partial }{\partial \phi }\left(P{H}^3\frac{\partial P}{\partial \phi }\right)+\frac{\partial }{\partial \lambda }\left(P{H}^3\frac{\partial P}{\partial \lambda }\right)=\mathrm{\Lambda }\frac{\partial \left(PH\right)}{\partial \phi }$$

(1)

where φ and λ are the circumferential and axial coordinates in the Cartesian coordinate system; P(φ, λ) is the nondimensional gas film pressure; H(φ) is the nondimensional gas film thickness; and Λ is the compressible bearing number.

The nondimensional boundary conditions of the Reynolds equation are

$$\{\begin{array}{l}P(\phi ,\lambda )=P(\phi +2\pi ,\lambda )\\ P(\phi ,\lambda )=P(\phi ,-\lambda )=0\\ {\frac{\partial P}{\partial \lambda }|}_{\lambda =0}=0\end{array}$$

(2)

The steady Reynolds equation is solved by FDM (finite difference method).

Dynamic coefficients of gas-lubricated bearing

The static gas film pressure P₀ of gas-lubricated bearing at steady equilibrium position is solved by computing iteratively the steady Reynolds equation. The dynamic gas film pressure P _d is caused by the journal deviating from steady equilibrium position due to the small perturbation. The position of the journal at any time can be expressed as

$$\{\begin{array}{c}\epsilon ={\epsilon }_0+\mathit{E}={\epsilon }_0+{\mathit{E}}_0{e}^{\text{j}\mathrm{\Omega }T}\\ \theta ={\theta }_0+\mathrm{\Theta }={\theta }_0+{\mathrm{\Theta }}_0{e}^{j\mathrm{\Omega }T}\end{array}$$

(3)

where E₀ is the perturbation amplitude of the eccentricity; Θ₀ is the perturbation amplitude of the attitude angle; ε ₀ and θ₀ are the eccentricity and attitude angle at steady equilibrium position; j is an imaginary number; Ω is the nondimensional perturbation frequency, Ω=ν/ω; and ν is the perturbation frequency.

With the small perturbation of the journal, the dynamic gas film thickness and pressure between gas-lubricated bearing and journal can be written as

$$\{\begin{array}{c}H=H_0+H_{d0}{e}^{\mathrm{j}\mathrm{\Omega }T}=\left(1+{\epsilon }_0\cos \phi \right)+\left(E_0\cos \phi +{\mathrm{\Theta }}_0{\epsilon }_0\sin \phi \right)e^{\mathrm{j}\mathrm{\Omega }T}\\ P=P_0+P_{d0}{e}^{\mathrm{j}\mathrm{\Omega }T}\end{array}$$

(4)

where H₀ and P₀ are the static gas film thickness and gas film pressure at steady equilibrium position; H_d0 is the perturbation amplitude of gas film thickness; and P_d0 is the perturbation amplitude of gas film pressure.

The differential of dynamic gas film forces with respect to perturbation amplitude of eccentricity and attitude angle can be expressed as

$$\{\begin{array}{l}{P}_E=\frac{\partial P_{d0}}{\partial E_0}\hfill \\ P_{\mathrm{\Theta }}=\frac{1}{{\epsilon }_0}\frac{\partial P_{d0}}{\partial {\mathrm{\Theta }}_0}\hfill \\ H_E=\frac{\partial H_{d0}}{\partial E_0}\hfill \\ H_{\mathrm{\Theta }}=\frac{1}{{\epsilon }_0}\frac{\partial H_{d0}}{\partial {\mathrm{\Theta }}_0}\hfill \end{array}$$

(5)

By substituting the equation (4) into the Reynolds equation, and taking the derivative with respect to E₀ and Θ₀, the partial differential equations (6) and (7) can be obtained. P _E and P_Θ can be calculated by solving the partial differential equations

$$\begin{array}{l}\frac{\partial }{\partial \phi }\left(P_0{H}_0^3\frac{\partial P_E}{\partial \phi }\right)+\frac{\partial }{\partial \lambda }\left(P_0{H}_0^3\frac{\partial P_E}{\partial \lambda }\right)+\frac{\partial }{\partial \phi }\left(P_E{H}_0^3\frac{\partial P_0}{\partial \phi }\right)+\frac{\partial }{\partial \lambda }\left(P_E{H}_0^3\frac{\partial P_0}{\partial \lambda }\right)+\frac{\partial }{\partial \phi }\left(3H_0^2{P}_0\frac{\partial P_0}{\partial \phi }{H}_E\right)+\frac{\partial }{\partial \lambda }\left(3H_0^2{P}_0\frac{\partial P_0}{\partial \lambda }{H}_E\right)=\mathrm{\Lambda }\frac{\partial }{\partial \phi }\left(H_0{P}_E+P_0{H}_E\right)+2\mathrm{j}\mathrm{\Lambda }\mathrm{\Omega }\left(H_0{P}_E+P_0{H}_E\right)\end{array}$$

(6)

$$\begin{array}{l}\frac{\partial }{\partial \phi }\left(P_0{H}_0^3\frac{\partial P_{\mathrm{\Theta }}}{\partial \phi }\right)+\frac{\partial }{\partial \lambda }\left(P_0{H}_0^3\frac{\partial P_{\mathrm{\Theta }}}{\partial \lambda }\right)+\frac{\partial }{\partial \phi }\left(P_{\mathrm{\Theta }}{H}_0^3\frac{\partial P_0}{\partial \phi }\right)+\frac{\partial }{\partial \lambda }\left(P_{\mathrm{\Theta }}{H}_0^3\frac{\partial P_0}{\partial \lambda }\right)\\ +\frac{\partial }{\partial \phi }\left(3H_0^2{P}_0\frac{\partial P_0}{\partial \phi }{H}_{\mathrm{\Theta }}\right)+\frac{\partial }{\partial \lambda }\left(3H_0^2{P}_0\frac{\partial P_0}{\partial \lambda }{H}_{\mathrm{\Theta }}\right)=\mathrm{\Lambda }\frac{\partial }{\partial \phi }\left(H_0{P}_{\mathrm{\Theta }}+P_0{H}_{\mathrm{\Theta }}\right)+2\mathrm{j}\mathrm{\Lambda }\mathrm{\Omega }\left(H_0{P}_{\mathrm{\Theta }}+P_0{H}_{\mathrm{\Theta }}\right)\end{array}$$

(7)

where $H_E=\cos \phi$ and $H_{\mathrm{\Theta }}=\sin \phi$ .

The dynamic stiffness and damping coefficients of gas-lubricated bearing can be expressed as equation (8) by the differential of dynamic gas film forces P _E and P_Θ. The dynamic stiffness and damping coefficients in the x and y directions can be obtained by the transformation matrix

$$\{\begin{array}{l}-{\displaystyle {\int }_{-\frac{B}{2R}}^{\frac{B}{2R}}{\displaystyle {\int }_0^{2\pi }{P}_E}\cos \phi d\phi d\lambda =K_{y\epsilon }+\mathrm{j}\mathrm{\Omega }{D}_{y\epsilon }}\\ -{\displaystyle {\int }_{-\frac{B}{2R}}^{\frac{B}{2R}}{\displaystyle {\int }_0^{2\pi }{P}_E}\sin \phi d\phi d\lambda =K_{x\epsilon }+\mathrm{j}\mathrm{\Omega }{D}_{x\epsilon }}\\ -{\displaystyle {\int }_{-\frac{B}{2R}}^{\frac{B}{2R}}{\displaystyle {\int }_0^{2\pi }{P}_{\mathrm{\Theta }}}\cos \phi d\phi d\lambda =K_{y\theta }+\mathrm{j}\mathrm{\Omega }{D}_{y\theta }}\\ -{\displaystyle {\int }_{-\frac{B}{2R}}^{\frac{B}{2R}}{\displaystyle {\int }_0^{2\pi }{P}_{\mathrm{\Theta }}}\sin \phi d\phi d\lambda =K_{x\theta }+\mathrm{j}\mathrm{\Omega }{D}_{x\theta }}\end{array}$$

(8)

$$\{\begin{array}{l}\left[\begin{array}{c}{K}_{xx}\\ K_{xy}\\ K_{yx}\\ K_{yy}\end{array}\right]=\left[\begin{array}{cccc}\sin {\theta }_0& \cos {\theta }_0& 0& 0\\ \cos {\theta }_0& -\sin {\theta }_0& 0& 0\\ 0& 0& \sin {\theta }_0& \cos {\theta }_0\\ 0& 0& \cos {\theta }_0& -\sin {\theta }_0\end{array}\right]\left[\begin{array}{c}{K}_{x\epsilon }\\ K_{x\theta }\\ K_{y\epsilon }\\ K_{y\theta }\end{array}\right]\\ \left[\begin{array}{c}{D}_{xx}\\ D_{xy}\\ D_{yx}\\ D_{yy}\end{array}\right]=\left[\begin{array}{cccc}\sin {\theta }_0& \cos {\theta }_0& 0& 0\\ \cos {\theta }_0& -\sin {\theta }_0& 0& 0\\ 0& 0& \sin {\theta }_0& \cos {\theta }_0\\ 0& 0& \cos {\theta }_0& -\sin {\theta }_0\end{array}\right]\left[\begin{array}{c}{D}_{x\epsilon }\\ D_{x\theta }\\ D_{y\epsilon }\\ D_{y\theta }\end{array}\right]\end{array}$$

(9)

The impedance of gas film forces can be obtained by a complex number form of dynamic stiffness and damping coefficients

$$Z\left(s\right)=Z\left(\mathrm{j}\mathrm{\Omega }\right)=K+\mathrm{j}\mathrm{\Omega }D$$

(10)

Rational function approximation of impedance

Impedance of the gas-lubricated bearing

According to the definition of impedance, the impedance Z(s) can be expressed as the differential of the perturbation gas film force versus small perturbation amplitude

$$Z\left(s\right)=\frac{\mathrm{\Delta }F\left(s\right)}{\mathrm{\Delta }\mathrm{\Xi }\left(s\right)}$$

(11)

where ΔF(s) is the increment of perturbation gas film force; ΔΞ(s) is the increment of perturbation displacement; and s is the variable of the Laplace transform, s=jΩ.

Fitting of the rational function coefficients

According to the response relationship between ΔF(s) and ΔΞ(s), the transfer function can be expressed as rational fractions

$$G\left(s\right)=\frac{\mathrm{\Delta }\mathrm{\Xi }\left(s\right)}{\mathrm{\Delta }F\left(s\right)}\approx \frac{{\displaystyle \sum _{i=0}^n{b}_i{s}^i}}{{\displaystyle \sum _{i=0}^n{a}_i{s}^i}}$$

(12)

where a _i and b _i are the unknown rational function coefficients.

By combining equation (11) with equation (12), the approximate expression of the impedance Z(s) using rational function is as follows

$$Z\left(s\right)=\frac{\mathrm{\Delta }F\left(s\right)}{\mathrm{\Delta }\mathrm{\Xi }\left(s\right)}=\frac{1}{G\left(s\right)}\approx \frac{{\displaystyle \sum _{i=0}^n{a}_i{s}^i}}{{\displaystyle \sum _{i=0}^n{b}_i{s}^i}}$$

(13)

According to partial fraction expansion, the rational function can be expressed as the sum of several rational fractional functions, and the impedance Z(s) is in the following form

$$Z\left(s\right)\approx {\displaystyle \sum _{i=1}^n\frac{c_i}{s-p_i}}+d$$

(14)

where the unknown coefficients c _i and d are rational function coefficients and remainder, respectively, and p _i are the unknown poles of the rational fractional functions. The coefficients a _i and b _i of equation (13) can be obtained by calculating unknown parameters c _i , d, and p _i of equation (14), which are obtained by the vector fitting method.

For the rational fractional functions, taking n = 2, the impedance Z(s) is approximated as a second-order rational function in the following form

$$Z\left(s\right)\approx {\displaystyle \sum _{i=1}^2\frac{c_i}{s-p_i}}+d=\frac{c_1}{s-p_1}+\frac{c_2}{s-p_2}+d$$

(15)

The calculation process of rational function unknown coefficients is as follows:

1. A set of initial poles ${\tilde{p}}_i$ is assumed to solve the real poles p_i of impedance Z(s).

The initial poles ${\tilde{p}}_i$ should be taken as conjugate complex numbers with real part α and imaginary part β, ³⁰ which distribute linearly over the given frequency range s as follows

$$\{\begin{array}{c}{\tilde{p}}_i=-\alpha +\mathrm{j}\beta \\ {\tilde{p}}_i=-\alpha -\mathrm{j}\beta \end{array}$$

(16)

where α=β/100.

• 2. Calculation of the unknown coefficients of impedance Z(s).

The rational function δ(s) with the same poles ${\tilde{\mathit{p}}}_i$ can be expressed as

$$\delta \left(s\right)\approx {\displaystyle \sum _{i=1}^2\frac{{\tilde{c}}_i}{s-{\tilde{p}}_i}}+1$$

(17)

where the unknown coefficients ${\tilde{\mathit{c}}}_i$ are the rational function coefficients, which are different from c _i .

By multiplying the rational functions δ(s) and Z(s), the following equation can be obtained

$$\delta \left(s\right)Z\left(s\right)\approx \left({\displaystyle \sum _{i=1}^2\frac{{\tilde{c}}_i}{s-{\tilde{p}}_i}}+1\right)Z\left(s\right)$$

(18)

where Z(s) is a series of impedances at the given frequency range s .

δ(s)Z(s) can be approximated to a new rational function which is similar to equation (15)

$$\delta \left(s\right)Z\left(s\right)\approx {\displaystyle \sum _{i=1}^2\frac{c_i}{s-{\tilde{p}}_i}}+d$$

(19)

Equation (19) can be written as

$$\tilde{Z}\left(s\right)={\displaystyle \sum _{i=1}^2\frac{c_i}{s-{\tilde{p}}_i}}+d$$

(20)

where equation (20) can be solved because the poles ${\tilde{\mathit{p}}}_i$ in the denominator are known.

Equation (20) can be expressed as the overdetermined equation

$$\left[\begin{array}{ccccc}\frac{1}{\mathit{s}-{\tilde{p}}_1}& \frac{1}{\mathit{s}-{\tilde{p}}_2}& \frac{-\mathit{Z}\left(\mathit{s}\right)}{\mathit{s}-{\tilde{p}}_1}& \frac{-\mathit{Z}\left(\mathit{s}\right)}{\mathit{s}-{\tilde{p}}_2}& 1\end{array}\right]{\left\{\begin{array}{ccccc}{c}_1& c_2& {\tilde{c}}_1& {\tilde{c}}_2& d\end{array}\right\}}^T=\mathit{Z}\left(\mathit{s}\right)$$

(21)

Z(s) is known, and the unknown coefficients c _i , ${\tilde{\mathit{c}}}_i$ , and d can be solved by the least square method.

By setting the initial poles ${\tilde{\mathit{p}}}_{\mathit{i}}{=\tilde{z}}_{\mathit{i}}$ , the real poles p _i can be calculated iteratively. The coefficients of impedance Z(s) are all known. (The poles p _i of equation (15) are equal to the zeros ${\tilde{z}}_i$ of δ(s), and the corresponding derivation is shown in Appendix 1).

3) To solve the rational function coefficients of impedance Z(s).

According to equation (13), the second-order rational function fraction of impedance Z(s) is

$$Z\left(s\right)=\frac{a_2{s}^2+a_{1}s+a_0}{b_2{s}^2+b_{1}s+b_0}$$

(22)

By solving equations (15) and (22) simultaneously, the following equation can be obtained

$${\displaystyle \sum _{i=1}^2\frac{c_i}{s-p_i}}+d=\frac{a_2{s}^2+a_{1}s+a_0}{b_2{s}^2+b_{1}s+b_0}$$

the rational function coefficients a _i and b _i (i=0, 1, 2) can be calculated.

Z(s)=Z(jΩ) is approximated by second-order rational fractions, according to equations (10) and (22), and the stiffness and damping coefficients can be expressed as

$$K+\mathrm{j}\mathrm{\Omega }D\approx \frac{{\displaystyle \sum _{i=0}^2{a}_i{\left(\mathrm{j}\mathrm{\Omega }\right)}^i}}{{\displaystyle \sum _{i=0}^2{b}_i{\left(\mathrm{j}\mathrm{\Omega }\right)}^i}}=\frac{a_0+a_1\mathrm{j}\mathrm{\Omega }-a_2{\mathrm{\Omega }}^2}{b_0+b_1\mathrm{j}\mathrm{\Omega }-b_2{\mathrm{\Omega }}^2}$$

(23)

where the normalized rational function coefficient b₂ = 1.

The stiffness coefficient K and damping coefficient D ³¹ can be expressed as the following

$$K\left(\mathrm{\Omega }\right)=\frac{-a_2{b}_{\text{0}}{\mathrm{\Omega }}^2+a_2{\mathrm{\Omega }}^4+a_1{b}_1{\omega }^2+a_0{b}_0-a_0{\mathrm{\Omega }}^2}{{\left(b_0-{\mathrm{\Omega }}^2\right)}^2+{\left(b_1\mathrm{\Omega }\right)}^2}$$

(24)

$$D\left(\mathrm{\Omega }\right)=\frac{a_2{b}_1{\mathrm{\Omega }}^2+a_1{b}_{\text{0}}-a_1{\mathrm{\Omega }}^2-a_0{b}_1}{{\left(b_0-{\mathrm{\Omega }}^2\right)}^2+{\left(b_1\mathrm{\Omega }\right)}^2}$$

(25)

Hybrid model of the bearing–rotor system

Motion equation of the rotor

Figure 2 shows the schematic diagram of the asymmetric rotor system supported in axial-grooved gas-lubricated bearings. The asymmetric rotor is composed of five segments of different lengths and diameters; the gas-lubricated bearings 1 and 2 are mounted on the segments 1 and 5 separately. O₁ and O₂ are the journal centers of gas-lubricated bearings 1 and 2 separately; O is the mass center of asymmetric rotor; m is the mass of rotor; l₁ and l₂ are the distances between the mass center and the journal centers at bearing 1 and 2 stations separately, l = l₁+l₂.

Figure 2.

Schematic diagram of gas-lubricated bearing–asymmetric rotor system.

The motion equation of asymmetric rotor system is

$$\{\begin{array}{l}\frac{m{l}_2}{l}{\ddot{x}}_1+\frac{m{l}_1}{l}{\ddot{x}}_2=m{e}_x{\omega }^2\cos \left(\omega t\right)+m{e}_y{\omega }^2\sin \left(\omega t\right)-f_{x1}-f_{x2}\\ \frac{m{l}_2}{l}{\ddot{y}}_1+\frac{m{l}_1}{l}{\ddot{y}}_2=m{e}_y{\omega }^2\cos \left(\omega t\right)-m{e}_x{\omega }^2\sin \left(\omega t\right)+mg-f_{y1}-f_{y2}\\ -\frac{I_y}{l}{\ddot{x}}_1+\frac{I_y}{l}{\ddot{x}}_2+\frac{I_G\omega }{l}{\dot{y}}_1-\frac{I_G\omega }{l}{\dot{y}}_2=f_{x1}{l}_1-f_{x2}{l}_2\\ -\frac{I_x}{l}{\ddot{y}}_1+\frac{I_x}{l}{\ddot{y}}_2-\frac{I_G\omega }{l}{\dot{x}}_1+\frac{I_G\omega }{l}{\dot{x}}_2=f_{y1}{l}_1-f_{y2}{l}_2\end{array}$$

(26)

where m is the mass of rotor; x₁, x₂, y₁, and y₂ are the displacements of the journal centers in the x and y directions at bearing 1 and 2 stations separately; ${\ddot{\mathit{x}}}_1$ , ${\ddot{\mathit{x}}}_2$ , ÿ₁, and ÿ₂ are the accelerated velocities of the journal centers in the x and y directions at bearing 1 and 2 stations separately; ${\dot{\mathit{x}}}_1$ , ${\dot{\mathit{x}}}_2$ , ${\dot{\mathit{y}}}_1$ , and ${\dot{\mathit{y}}}_2$ are the velocities of the journal centers in the x and y directions at bearing 1 and 2 stations separately; e _x and e _y are the mass eccentricities of rotor in the x and y directions separately; f_x1, f_x2, f_y1, and f_y2 are the gas film forces of the bearings 1 and 2 in the x and y directions separately; I _x and I _y are the equatorial moments of the inertia of the rotor in the x and y directions separately; and I _G is the polar moment of the inertia of the rotor. g is the gravity acceleration and τ is the time.

By introducing the following nondimensional variables, P_a is the ambient pressure and R is the radius of the journal

$$\overline{x}=\frac{x}{c},\hspace{0.17em}\overline{y}=\frac{y}{c},\hspace{0.17em}{\overline{x}}^{\prime }=\frac{\dot{\mathit{x}}}{\omega c},\hspace{0.17em}{\overline{y}}^{\prime }=\frac{\dot{\mathit{y}}}{\omega c},\hspace{0.17em}{\overline{x}}^{″}=\frac{\ddot{\mathit{x}}}{{\omega }^{2}c},\hspace{0.17em}{\overline{y}}^{″}=\frac{\ddot{y}}{{\omega }^{2}c},\hspace{0.17em}{\overline{\epsilon }}_y=\frac{e_y}{c},\hspace{0.17em}{\overline{\epsilon }}_x=\frac{e_x}{c},\hspace{0.17em}{\overline{I}}_x=\frac{I_x{\omega }^{2}c}{P_a{R}^2},\hspace{0.17em}{\overline{I}}_y=\frac{I_y{\omega }^{2}c}{P_a{R}^2},\hspace{0.17em}{\overline{I}}_G=\frac{I_G{\omega }^{2}c}{P_a{R}^2},\hspace{0.17em}{\overline{f}}_x=\frac{f_x}{P_a{R}^2},\hspace{0.17em}{\overline{f}}_y=\frac{f_y}{P_a{R}^2},\hspace{0.17em}\overline{m}=\frac{m{\omega }^{2}c}{P_a{R}^2},\overline{g}=\frac{g}{{\omega }^{2}c},\tau =\omega t$$

The nondimensional equation of asymmetric rotor system can be obtained as follows

$$\{\begin{array}{l}\frac{\overline{m}{l}_2}{l}{{\overline{x}}^{″}}_1+\frac{\overline{m}{l}_1}{l}{{\overline{x}}^{″}}_2=\overline{m} {\overline{e}}_x\cos \tau +\overline{m} {\overline{e}}_y\sin \tau -{\overline{f}}_{x1}-{\overline{f}}_{x2}\\ \frac{\overline{m}{l}_2}{l}{{\overline{y}}^{″}}_1+\frac{\overline{m}{l}_1}{l}{{\overline{y}}^{″}}_2=\overline{m} {\overline{e}}_y\cos \tau -\overline{m} {\overline{e}}_x\sin \tau +\overline{m} \overline{g}-{\overline{f}}_{y1}-{\overline{f}}_{y2}\\ -\frac{{\overline{I}}_y}{l}{{\overline{x}}^{″}}_1+\frac{{\overline{I}}_y}{l}{{\overline{x}}^{″}}_2+\frac{{\overline{I}}_G}{l}{{\overline{y}}^{\prime }}_1-\frac{{\overline{I}}_G}{l}{{\overline{y}}^{\prime }}_2={\overline{f}}_{x1}{l}_1-{\overline{f}}_{x2}{l}_2\\ -\frac{{\overline{I}}_x}{l}{{\overline{y}}^{″}}_1+\frac{{\overline{I}}_x}{l}{{\overline{y}}^{″}}_2-\frac{{\overline{I}}_G}{l}{{\overline{x}}^{\prime }}_1+\frac{{\overline{I}}_G}{l}{{\overline{x}}^{\prime }}_2={\overline{f}}_{y1}{l}_1-{\overline{f}}_{y2}{l}_2\end{array}$$

(27)

In order to calculate the displacements of the journal centers of motion equation (27), the unsteady Reynolds equation and nonlinear gas film forces need to be solved at each step. The multiple iteration results in higher computing cost.

Rational function model of gas film forces

Equation (13) is the rational function equation in frequency domain, and the relation between disturbance increment and response increment can be transformed into time domain by inverse Fourier transform; the equation in the time domain can be obtained as follows

$${\displaystyle \sum _{i=0}^n{b}_i^{\alpha \beta }}\mathrm{\Delta }{f}_{\alpha \beta }{}^{\left(i\right)}-{\displaystyle \sum _{i=0}^n{a}_i^{\alpha \beta }}\mathrm{\Delta }{\xi }^{\left(i\right)}=0$$

(28)

where a _i ^αβ and b _i ^αβ represent the coefficients of the direct and the cross coupled impedances, α, β∈{x, y}; Δξ is the displacement increment, which represents Δx or Δy; Δf _αβ is the gas film force increment; and (i) represents the ith order differential.

Because the direct and cross coupled impedances share same poles, b _i ^xx =b _i ^xy =b _i ^yx =b _i ^yy =b _i . By expanding equation (28) in the x and y directions, equation (29) can be obtained

$$\{\begin{array}{l}{\displaystyle \sum _{i=0}^n\left[b_i\left(\mathrm{\Delta }{f}_{xx}{}^{\left(i\right)}+\mathrm{\Delta }{f}_{xy}{}^{\left(i\right)}\right)-a_i^{xx}\mathrm{\Delta }{x}^{\left(i\right)}-a_i^{xy}\mathrm{\Delta }{y}^{\left(i\right)}\right]}=0\\ {\displaystyle \sum _{i=0}^n\left[b_i\left(\mathrm{\Delta }{f}_{\text{y}x}{}^{\left(i\right)}+\mathrm{\Delta }{f}_{yy}{}^{\left(i\right)}\right)-a_i^{yx}\mathrm{\Delta }{x}^{\left(i\right)}-a_i^{yy}\mathrm{\Delta }{y}^{\left(i\right)}\right]}=0\end{array}$$

(29)

According to the relationship of gas film force components in the x and y directions, equation (30) can be expressed as

$$\{\begin{array}{l}\mathrm{\Delta }{f}_x=\mathrm{\Delta }{f}_{xx}+\mathrm{\Delta }{f}_{xy}\\ \mathrm{\Delta }{f}_y=\mathrm{\Delta }{f}_{yx}+\mathrm{\Delta }{f}_{yy}\end{array}\to \{\begin{array}{l}{\displaystyle \sum _{i=0}^n\left(b_i\mathrm{\Delta }{f}_x{}^{\left(i\right)}-a_i^{xx}\mathrm{\Delta }{x}^{\left(i\right)}-a_i^{xy}\mathrm{\Delta }{y}^{\left(i\right)}\right)}=0\\ {\displaystyle \sum _{i=0}^n\left(b_i\mathrm{\Delta }{f}_y{}^{\left(i\right)}-a_i^{yx}\mathrm{\Delta }{x}^{\left(i\right)}-a_i^{yy}\mathrm{\Delta }{y}^{\left(i\right)}\right)}=0\end{array}$$

(30)

The rational function model of gas-lubricated bearing is established by equation (30), which is coupled with journal displacements. The second-order expansion equations at bearing 1 and 2 stations are

$$\{\begin{array}{l}{f^{¨}}_{x1}+b_{1\text{,}1}{f^{\dot{}}}_{x1}+b_{\text{0,}1}\left(f_{x1}-f_{sx1}\right)-a_{2\text{,}1}^{xx}{x^{¨}}_1-a_{1\text{,}1}^{xx}{x^{\dot{}}}_1-a_{0\text{,}1}^{xx}\left(x_1-x_{s1}\right)-a_{2\text{,}1}^{xy}{y^{¨}}_1-a_{1\text{,}1}^{xy}{y^{\dot{}}}_1-a_{0\text{,}1}^{xy}\left(y_1-y_{s1}\right)=0\\ {f^{¨}}_{x2}+b_{1\text{,}2}{f^{\dot{}}}_{x2}+b_{\text{0,}2}\left(f_{x2}-f_{sx2}\right)-a_{2\text{,}2}^{xx}{x^{¨}}_2-a_{1\text{,}2}^{xx}{x^{\dot{}}}_2-a_{0\text{,}2}^{xx}\left(x_2-x_{s2}\right)-a_{2\text{,}2}^{xy}{y^{¨}}_2-a_{1\text{,}2}^{xy}{y^{\dot{}}}_2-a_{0\text{,}2}^{xy}\left(y_2-y_{s2}\right)=0\\ {f^{¨}}_{y1}+b_{1\text{,}1}{f^{\dot{}}}_{y1}+b_{\text{0,}1}\left(f_{y1}-f_{sy1}\right)-a_{2\text{,}1}^{yx}{x^{¨}}_1-a_{1\text{,}1}^{yx}{x^{\dot{}}}_1-a_{0\text{,}1}^{yx}\left(x_1-x_{s1}\right)-a_{2\text{,}1}^{yy}{y^{¨}}_1-a_{1\text{,}1}^{yy}{y^{\dot{}}}_1-a_{0\text{,}1}^{yy}\left(y_1-y_{s1}\right)=0\\ {f^{¨}}_{y2}+b_{1\text{,}2}{f^{\dot{}}}_{y2}+b_{\text{0,}2}\left(f_{y2}-f_{sy2}\right)-a_{2\text{,}2}^{yx}{x^{¨}}_2-a_{1\text{,}2}^{yx}{x^{\dot{}}}_2-a_{0\text{,}2}^{yx}\left(x_2-x_{s2}\right)-a_{2\text{,}2}^{yy}{y^{¨}}_2-a_{1\text{,}2}^{yy}{y}^{\dot{}}{}_2-a_{0\text{,}2}^{yy}\left(y_2-y_{s2}\right)=0\end{array}$$

(31)

where x _s , y _s , f _sx , and f _sy are the steady displacements and gas film forces in the x and y directions; Δf is replaced by the difference between transient and steady gas film forces; and Δξ is replaced by the difference between transient and steady displacements.

Hybrid model of gas-lubricated bearing–rotor system

Equation (31) is converted into nondimensional form and coupled with the nondimensional motion equation of rotor. A hybrid numerical model of the rotor supported in gas-lubricated bearing is established. The hybrid dynamic equation is

$$\mathit{M}\ddot{\mathit{u}}+\mathit{G}\dot{\mathit{u}}+\mathit{K}\mathit{u}=\mathit{P}=\mathit{Q}+\mathit{F}+\mathit{W}$$

(32)

where the variable vector is

$$\mathit{u}={\left\{{\overline{x}}_1,{\overline{y}}_1,{\overline{x}}_2,{\overline{y}}_2,{\overline{f}}_{x1},{\overline{f}}_{y1},{\overline{f}}_{x2},{\overline{f}}_{y2}\right\}}^T$$

M is the nondimensional mass matrix; G is the nondimensional gyroscopic matrix; K is the nondimensional stiffness matrix; and P is the nondimensional resultant force vector, which is composed of the external excitation force vector Q , gas film forces vector F , and gravity vector W . The entities of matrixes are shown as follows

$$\mathit{M}=\left[\begin{array}{cccccccc}\frac{\overline{m}{l}_2}{l}& 0& \frac{\overline{m}{l}_1}{l}& 0& 0& 0& 0& 0\\ 0& \frac{\overline{m}{l}_2}{l}& 0& \frac{\overline{m}{l}_1}{l}& 0& 0& 0& 0\\ -\frac{{\overline{I}}_y}{l}& 0& \frac{{\overline{I}}_y}{l}& 0& 0& 0& 0& 0\\ 0& -\frac{{\overline{I}}_x}{l}& 0& \frac{{\overline{I}}_x}{l}& 0& 0& 0& 0\\ -a_{\mathrm{2,1}}^{xx}{\omega }^{2}c& -a_{\mathrm{2,1}}^{xy}{\omega }^{2}c& 0& 0& P_a{R}^2& 0& 0& 0\\ 0& 0& -a_{\mathrm{2,2}}^{xx}{\omega }^{2}c& -a_{\mathrm{2,2}}^{xy}{\omega }^{2}c& 0& P_a{R}^2& 0& 0\\ -a_{\mathrm{2,1}}^{yx}{\omega }^{2}c& -a_{\mathrm{2,1}}^{yy}{\omega }^{2}c& 0& 0& 0& 0& P_a{R}^2& 0\\ 0& 0& -a_{\mathrm{2,2}}^{yx}{\omega }^{2}c& -a_{\mathrm{2,2}}^{yy}{\omega }^{2}c& 0& 0& 0& P_a{R}^2\end{array}\right]$$

(33)

$$\mathit{G}=\left[\begin{array}{cccccccc}0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& \frac{{\overline{I}}_G}{l}& 0& -\frac{{\overline{I}}_G}{l}& 0& 0& 0& 0\\ -\frac{{\overline{I}}_G}{l}& 0& \frac{{\overline{I}}_G}{l}& 0& 0& 0& 0& 0\\ -a_{\mathrm{1,1}}^{xx}\omega c& -a_{\mathrm{1,1}}^{xy}\omega c& 0& 0& b_{\mathrm{1,1}}{P}_a{R}^2& 0& 0& 0\\ 0& 0& -a_{\mathrm{1,2}}^{xx}\omega c& -a_{\mathrm{1,2}}^{xy}\omega c& 0& 0& b_{\mathrm{1,2}}{P}_a{R}^2& 0\\ -a_{\mathrm{1,1}}^{yx}\omega c& -a_{\mathrm{1,1}}^{yy}{\omega }^{2}c& 0& 0& 0& b_{\mathrm{1,1}}{P}_a{R}^2& 0& 0\\ 0& 0& -a_{\mathrm{1,2}}^{yx}\omega c& -a_{\mathrm{1,2}}^{yx}\omega c& 0& 0& 0& b_{\mathrm{1,2}}{P}_a{R}^2\end{array}\right]$$

(34)

$$\mathit{K}=\left[\begin{array}{cccccccc}0& 0& 0& 0& 1& 0& 1& 0\\ 0& 0& 0& 0& 0& 1& 0& 1\\ 0& 0& 0& 0& -l_1& 0& l_2& 0\\ 0& 0& 0& 0& 0& -l_1& 0& l_2\\ -a_{\mathrm{0,1}}^{xx}c& -a_{\mathrm{0,1}}^{xy}c& 0& 0& b_{\mathrm{0,1}}{P}_a{R}^2& 0& 0& 0\\ 0& 0& -a_{\mathrm{0,2}}^{xx}c& -a_{\mathrm{0,2}}^{xy}c& 0& 0& b_{\mathrm{0,2}}{P}_a{R}^2& 0\\ -a_{\mathrm{0,1}}^{yx}c& -a_{\mathrm{0,1}}^{yy}c& 0& 0& 0& b_{\mathrm{0,1}}{P}_a{R}^2& 0& 0\\ 0& 0& -a_{\mathrm{0,2}}^{yx}c& -a_{\mathrm{0,2}}^{yy}c& 0& 0& 0& b_{\mathrm{0,2}}{P}_a{R}^2\end{array}\right]$$

(35)

$$\mathit{P}=\left[\begin{array}{c}\overline{m} \overline{g} {\overline{e}}_x\cos \tau +\overline{m} \overline{g} {\overline{e}}_y\sin \tau \\ \overline{m} \overline{g} {\overline{e}}_y\cos \tau -\overline{m} \overline{g} {\overline{e}}_x\sin \tau \\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ b_{\mathrm{0,1}}{P}_a{R}^2{\overline{f}}_{sx1}+a_{\mathrm{0,1}}^{xx}c{\overline{x}}_{s1}+a_{\mathrm{0,1}}^{xy}c{\overline{y}}_{s1}\\ b_{\mathrm{0,2}}{P}_a{R}^2{\overline{f}}_{sx2}+a_{\mathrm{0,2}}^{xx}c{\overline{x}}_{s2}+a_{\mathrm{0,2}}^{xy}c{\overline{y}}_{s2}\\ b_{\mathrm{0,1}}{P}_a{R}^2{\overline{f}}_{sy1}+a_{\mathrm{0,1}}^{yx}c{\overline{x}}_{s1}+a_{\mathrm{0,1}}^{yy}c{\overline{y}}_{s1}\\ b_{\mathrm{0,2}}{P}_a{R}^2{\overline{f}}_{sy2}+a_{\mathrm{0,2}}^{yx}c{\overline{x}}_{s2}+a_{\mathrm{0,2}}^{yy}c{\overline{y}}_{s2}\end{array}\right]+\left[\begin{array}{c}0\\ \overline{m} \overline{g}\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right]=\left[\begin{array}{c}\overline{m} \overline{g} {\overline{e}}_x\cos \tau +\overline{m} \overline{g} {\overline{e}}_y\sin \tau \\ \overline{m} \overline{g} {\overline{e}}_y\cos \tau -\overline{m} \overline{g} {\overline{e}}_x\sin \tau +\overline{m} \overline{g}\\ 0\\ 0\\ b_{\mathrm{0,1}}{P}_a{R}^2{\overline{f}}_{sx1}+a_{\mathrm{0,1}}^{xx}c{\overline{x}}_{s1}+a_{\mathrm{0,1}}^{xy}c{\overline{y}}_{s1}\\ b_{\mathrm{0,2}}{P}_a{R}^2{\overline{f}}_{sx2}+a_{\mathrm{0,2}}^{xx}c{\overline{x}}_{s2}+a_{\mathrm{0,2}}^{xy}c{\overline{y}}_{s2}\\ b_{\mathrm{0,1}}{P}_a{R}^2{\overline{f}}_{sy1}+a_{\mathrm{0,1}}^{yx}c{\overline{x}}_{s1}+a_{\mathrm{0,1}}^{yy}c{\overline{y}}_{s1}\\ b_{\mathrm{0,2}}{P}_a{R}^2{\overline{f}}_{sy2}+a_{\mathrm{0,2}}^{yx}c{\overline{x}}_{s2}+a_{\mathrm{0,2}}^{yy}c{\overline{y}}_{s2}\end{array}\right]$$

(36)

The trajectories of journal centers and the gas film forces in the x and y directions at bearing 1 and 2 stations can be calculated by solving the hybrid numerical model simultaneously when the rational function coefficients a _i and b _i are known. In the computational process, the multiple iterations are avoided and the computing cost is reduced greatly. Because the rational function coefficients are same as long as the structural parameters of the bearing are unchanged, for the rotor system of different structures or parameters supported in same bearing, the rational function coefficients do not need to be recalculated.

Stability analysis

The motion equation of the rotor system can be written as

$$\{\begin{array}{l}\frac{m{l}_2}{l}{\ddot{x}}_1+\frac{m{l}_1}{l}{\ddot{x}}_2+k_{xx}{x}_1+k_{xy}{y}_1+k_{xx}{x}_2+k_{xy}{y}_2+d_{xx}{x^{\dot{}}}_1+d_{xy}{y^{\dot{}}}_1+d_{xx}{x^{\dot{}}}_2+d_{xy}{y^{\dot{}}}_2=0\\ \frac{m{l}_2}{l}{\ddot{y}}_1+\frac{m{l}_1}{l}{\ddot{y}}_2+k_{yx}{x}_1+k_{yy}{y}_1+k_{yx}{x}_2+k_{yy}{y}_2+d_{yx}{x^{\dot{}}}_1+d_{yy}{y^{\dot{}}}_1+d_{yx}{x^{\dot{}}}_2+d_{yy}{y^{\dot{}}}_2=0\\ -\frac{I_y}{l}{\ddot{x}}_1+\frac{I_y}{l}{\ddot{x}}_2+\frac{I_G\omega }{l}{y^{\dot{}}}_1-\frac{I_G\omega }{l}{y^{\dot{}}}_2-\left(k_{xx}{x}_1+k_{xy}{y}_1+d_{xx}{x^{·}}_1+d_{xy}{y^{·}}_1\right)l_1+\left(k_{xx}{x}_2+k_{xy}{y}_2+d_{xx}{x^{·}}_2+d_{xy}{y^{·}}_2\right)l_2=0\\ -\frac{I_x}{l}{\ddot{y}}_1+\frac{I_x}{l}{\ddot{y}}_2-\frac{I_G\omega }{l}{x^{\dot{}}}_1+\frac{I_G\omega }{l}{x^{\dot{}}}_2-{\left(k_{yx}{x}_1+k_{yy}{y}_1+d_{yx}{x^{·}}_1+d_{yy}{y}^{·}\right)}_1{l}_1+{\left(k_{yx}{x}_2+k_{yy}{y}_2+d_{yx}{x^{·}}_2+d_{yy}{y}^{·}\right)}_2{l}_2=0\end{array}$$

(37)

where k _xx , k _yy , k _xy , k _yx , d _xx , d _yy , d _xy , and d _yx are the stiffness and damping coefficients of bearing 1, respectively, which are same as bearing 2. The solution of motion equation (37) can be expressed as

$$\{\begin{array}{c}{x}_1=x_{10}{e}^{\kappa t}\\ y_1=y_{10}{e}^{\kappa t}\\ x_2=x_{20}{e}^{\kappa t}\\ y_2=y_{20}{e}^{\kappa t}\end{array}$$

(38)

where $\kappa =\gamma +\mathrm{j}{\omega }_n$ , κ is the complex frequency, γ is the real part, and ω _n is the imaginary part.

By substituting equation (38) into equation (37), equation (39) can be obtained

$$\left[\begin{array}{cccc}\frac{m{l}_2}{l}{\kappa }^2+k_{xx}+d_{xx}\kappa & k_{xy}+d_{xy}\kappa & \frac{m{l}_1}{l}{\kappa }^2+k_{xx}+d_{xx}\kappa & k_{xy}+d_{xy}\kappa \\ k_{yx}+d_{yx}\kappa & \frac{m{l}_2}{l}{\kappa }^2+k_{yy}+d_{yy}\kappa & k_{yx}+d_{yx}\kappa & \frac{m{l}_1}{l}{\kappa }^2+k_{yy}+d_{yy}\kappa \\ -\frac{I_y}{l}{\kappa }^2-k_{xx}{l}_1-d_{xx}\kappa l_1& \frac{I_G\omega }{l}\kappa -k_{xy}{l}_1-d_{xy}\kappa l_1& \frac{I_y}{l}{\kappa }^2+k_{xx}{l}_2+d_{xx}\kappa l_2& -\frac{I_G\omega }{l}\kappa +k_{xy}{l}_2+d_{xy}\kappa l_2\\ -\frac{I_G\omega }{l}\kappa -k_{yx}{l}_1-d_{yx}\kappa l_1& -\frac{I_x}{l}{\kappa }^2-k_{yy}{l}_1-d_{yy}\kappa l_1& \frac{I_G\omega }{l}\kappa +k_{yx}{l}_2+d_{yx}\kappa l_2& \frac{I_x}{l}{\kappa }^2+k_{yy}{l}_2+d_{yy}\kappa l_2\end{array}\right]\left[\begin{array}{c}{x}_{10}\\ y_{10}\\ x_{20}\\ y_{20}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right]$$

(39)

Equation (39) has untrivial solutions when the determinant equation of equation (39) equals to 0, as follows

$$A_8{\kappa }^8+A_7{\kappa }^7+A_6{\kappa }^6+A_5{\kappa }^5+A_4{\kappa }^4+A_3{\kappa }^3+A_2{\kappa }^2+A_1\kappa +A_0=0$$

(40)

where A₁∼ A₈ are the coefficients of determinant equation. According to the Routh–Hurwitz method, the stability of rotor system can be determined by the subdeterminant of the coefficient matrix of determinant equation (40). The coefficients of each order and the stability criterions are shown in Appendix 2.

Numerical results

In order to evaluate the performance of the gas-lubricated bearing–rotor system under different rotating speeds and mass eccentricities of rotor, the nonlinear dynamic responses of the rotor system are obtained and analyzed. The asymmetric rotor system supported in two axial-grooved gas-lubricated bearing is shown in Figure 2. The parameters of bearing 1 are identical with 2. The parameters of the two axial-grooved gas-lubricated bearing and asymmetric rotor are shown in Table 1.

Table 1.

The parameters of gas-lubricated bearing–rotor system.

·Parameters	Unit	Value
Journal radius, R	m	1×10−2
Width–diameter ratio, L/D		1
Radius clearence, c	m	7.5×10−6
Arc bushing angle, γ	deg	150
Groove width angle, γ0	deg	30
Mass eccentricities, e x =e y	m	5×10−7
Eccentricity ratio, e		0.28
Radiuses of each axis segment, R1, R2, R3, R4, R5	m	0.01, 0.012, 0.02, 0.13, 0.01
Lengths of each axis segment, l1, l2, l3, l4, l5	m	0.08, 0.1, 0.06, 0.09, 0.08
Density, ρ	Kg/m3	7.8×103
Ambient pressure, P a	MPa	0.1
Coefficient of kinetic viscosity, μ	Pa·s	1.8×10−5

Verification

Figure 3 is the time series of the journal centers at bearing 1 and 2 stations solved by the proposed hybrid numerical model and the ones obtained by solving the transient Reynolds equation and the rotor motion equation (27) by the Newmark-β method at rotating speed ω=1000 rad/s. By comparison of the trajectories of journal centers, the result obtained by the hybrid numerical model is in good agreement with the Newmark-β method. The computing time required by the hybrid numerical model and Newmark-β method is 521 s and 1972 s, respectively. It can be seen that the computing cost is reduced by the proposed hybrid numerical model.

Figure 3.

Comparison of time series of Y of journal centers at bearing 1 and 2 stations: (a) At bearing 1 station and (b) at bearing 2 station.

Dynamic response versus rotating speed

The nonlinear dynamic responses of asymmetric rotor system supported by two axial-grooved gas-lubricated bearings are complex and diverse with the change of rotating speed. The bifurcation characteristics of gas-lubricated bearing–rotor system for ω=500–3700 rad/s are shown in Figure 4. It can be seen that the bearing–rotor system has periodic, period-doubling, and quasi-periodic responses. When ω=500–1050 rad/s, the motions of the journal center are quasi-periodic at bearing 1 station. The amplitudes of the trajectories decrease with the increase of rotating speed. When ω=1100 rad/s, the quasi-periodic motion of journal center becomes periodic motion. When the rotating speed increases to 1850 rad/s, the periodic motion bifurcates to period-doubling motion. The period-doubling motion bifurcates inversely to periodic motion for ω=2250 rad/s. When ω=3160 rad/s, the periodic motion bifurcates to quasi-periodic motion, but when ω=3150, 3340, and 3350 rad/s, the motions of the journal center are period-14 motions. When the rotating speed increases to 3475 rad/s, the quasi-periodic motion of gas-lubricated bearing–rotor system turns into period-doubling motion.

Figure 4.

When ω = 500–3700 rad/s, the bifurcation diagram of journal center at bearing 1 station.

Figure 5(a) and (b) shows the motion bifurcations of gas-lubricated bearing–rotor system, which are from quasi-periodic motion to periodic motion and from quasi-periodic motion to period-doubling motion, respectively. It is worth noting the quasi-periodic motion bifurcates to multi-periodic motion and then turns into quasi-periodic motion; the quasi-periodic motion bifurcates to period-doubling motion with the increase of rotating speed.

Figure 5.

The bifurcation diagram of journal center at bearing 1 station: (a) ω = 900–1500 rad/s and (b) ω = 3400–3600 rad/s.

When ω=1000 rad/s, the quasi-periodic motions of the journal centers at bearing 1 and 2 stations are shown in Figure 6. Figure 6(a) and (b) shows the nondimensional trajectories of journal centers; as can be seen, the motions of journal centers are different at bearing 1 and 2 stations due to the gyroscopic effect. Figure 6(c) shows that the projection of the Poincaré map on the X-Y plane is torus attractor, and the spectrum is discrete in Figure 6(d). The time series of journal centers at bearing 1 and 2 stations are shown in Figure 6(e). For ω=1500 rad/s, the dynamic response of gas-lubricated bearing–rotor system is shown in Figure 7(a) at bearing 1 and 2 stations. As shown in Figure 7(b) and (c), the trajectory of journal center is periodic motion, and the spectrum is discrete.

Figure 6.

When ω = 1000 rad/s, the quasi-periodic motions of journal centers at bearing 1 and 2 stations: (a) Trajectories of journal centers; (b) comparison of trajectories of journal centers; (c) Poincaré map of journal center; (d) spectrum diagram of journal center; and (e) time series of journal centers.

Figure 7.

When ω = 1500 rad/s, the periodic motions of journal centers at bearing 1 and 2 stations: (a) Trajectories of journal centers; (b) trajectory of journal center at bearing 1 station; (c) spectrum diagram of journal center; and (d) time series of journal centers.

When ω=2000 rad/s, the periodic motions of journal centers bifurcate to period-doubling motions at bearing 1 and 2 stations, which are shown in Figure 8(a) and (b). The Poincaré map has two projection points on the X-Y plane, as shown in Figure 8(c). The spectrum diagram contains two spectral components, as shown in Figure 8(d). For ω=2500 rad/s, it can be seen that the period-doubling motions of the journal centers bifurcate inversely to periodic motions at bearing 1 and 2 stations in Figure 9. As shown in Figure 9(a)–(c), the trajectory of journal center is periodic motion, and the spectrum diagram contains one spectral component.

Figure 8.

When ω = 2000 rad/s, the period-doubling motions of journal centers at bearing 1 and 2 stations: (a) Trajectories of journal centers; (b) trajectory of journal center; (c) Poincaré map of journal center; (d) spectrum diagram of journal center; and (e) time series of journal centers.

Figure 9.

When ω = 2500 rad/s, the periodic motions of journal centers at bearing 1 and 2 stations: (a) Trajectories of journal centers; (b) trajectory of journal center; (c) spectrum diagram of journal center; and (d) time series of journal centers.

For ω=3400 rad/s, the dynamic response of gas-lubricated bearing–rotor system is shown in Figure 10. The periodic motions of journal centers bifurcate to quasi-periodic motions at bearing 1 and 2 stations, which are shown in Figure 10(a). The trajectory of journal center at bearing 1 station is shown in Figure 10(b), and the trajectory of journal center at bearing 2 station is similar with the trajectory at bearing 1 station. The Poincaré maps of journal centers at bearing 1 and 2 stations are shown in Figure 10(c)–(d), and the projections of the Poincaré maps on the X-Y plane are torus attractors. The spectrum diagram shows as discrete spectrum in Figure 10(e). When the rotating speed increases to 3420 and 3440 rad/s, the quasi-periodic motions become to period-10 motions shown in Figure 11. Figure 11(a)–(d) show the trajectories and Poincaré maps of journal centers at bearing 1 station; the Poincaré maps have 10 projection points on the X-Y plane in Figures 11(c) and (d), and the spectrum diagram is shown as discrete spectrum in Figure 11(e).

Figure 10.

When ω = 3400 rad/s, the quasi-periodic motions of journal centers at bearing 1 and 2 stations: (a) Trajectories of journal centers; (b) trajectory of journal center; (c) Poincaré map of journal center at bearing 1 station; (d) Poincaré map of journal center at bearing 2 station; (e) spectrum diagram of journal center; and (f) time series of journal centers.

Figure 11.

When ω=3420 and 3440 rad/s, the period-10 motions of journal center at bearing 1 station: (a) Trajectory of journal center for ω=3420 rad/s; (b) trajectory of journal center for ω = 3440 rad/s; (c) Poincaré map of journal center for ω=3420 rad/s; (d) Poincaré map of journal center for ω=3440 rad/s; and (e) spectrum diagram of journal center for ω = 3420 rad/s.

When ω=3500 rad/s, the dynamic response bifurcates inversely to period-doubling motion from quasi-periodic motion, as shown in Figure 12(a). The trajectories of journal centers at bearing 1 and 2 stations are shown in Figures 12(b) and (c). Figures 12(d) and (e) depict the Poincaré map and spectrum diagram of journal center at bearing 1 station; the projections of the Poincaré map on the X-Y plane are two fixed points; and the spectrum diagram contains two spectral components.

Figure 12.

When ω = 3500 rad/s, the period-doubling motions of journal centers at bearing 1 and 2 stations: (a) Trajectories of journal centers; (b) trajectory of journal center at bearing 1 station; (c) trajectory of journal center at bearing 2 station; (d) Poincaré map of journal center; (e) spectrum diagram of journal center; and (f) time series of journal centers.

Dynamic response versus mass eccentricity

As can be seen in Figure 13, the bifurcation diagram shows the nondimensional vertical motion of journal center at bearing 1 station, which is plotted versus unbalance mass eccentricities. The dynamic responses of gas-lubricated bearing–rotor system are periodic motions when the mass eccentricities e _x and e _y are less than 4.5×10⁻⁷ m, and the periodic motions of journal center at bearing 1 station bifurcate to quasi-periodic motions when e _x =e _y =4.5×10⁻⁷ m, with the increase of mass eccentricities.

Figure 13.

When ω = 1000 rad/s, e _x =e _y =0–9×10⁻⁷ m, the bifurcation diagram of journal center at bearing 1 station.

In Figures 14–18, when ω=1000 rad/s, the mass eccentricities e _x and e _y are 3×10⁻⁷, 5×10⁻⁷, 7×10⁻⁷, 9×10⁻⁷, and 1×10⁻⁶ m respectively; the dynamic responses of gas-lubricated bearing–rotor system are shown as periodic and quasi-periodic motions. When e _x =e _y =3×10⁻⁷ m, the trajectories of periodic motion of journal center are shown in Figure 14(a) and (b), and the Poincaré map has one projection point on the X-Y plane in Figure 14(c). In Figures 15–17, when e _x =e _y =5×10⁻⁷, 7×10⁻⁷, and 9×10⁻⁷ m, respectively, the dynamic responses are quasi-periodic motions. Figures 15(a), 16(a), and 17(a) show the trajectories of journal centers at bearing 1 and 2 stations with different eccentricity masses. The projection points of the Poincaré maps on the X-Y plane are torus attractors which are shown in Figures 15(c), 16(d)–(e) and 17(c). Figures 15(d), 16(f) and 17(d) show that the spectrum diagrams are discrete spectrums. Figure 18(a) and (b) show that the quasi-periodic motion bifurcates inversely to periodic motion when e _x =e _y =1×10⁻⁶ m. The spectrum diagram contains one spectral component shown in Figure 18(c).

Figure 14.

When e _x =e _y =3×10⁻⁷ m, the periodic motions of journal centers at bearing 1 and 2 stations: (a) Trajectories of journal centers; (b) trajectory of journal center; and (c) spectrum diagram of journal center.

Figure 15.

When e _x =e _y =5×10⁻⁷ m, the quasi-periodic motions of journal centers at bearing 1 and 2 stations: (a) Trajectories of journal centers; (b) trajectory of journal center; (c) Poincaré map of journal center; and (d) spectrum diagram of journal center.

Figure 16.

When e _x =e _y =7×10⁻⁷ m, the quasi-periodic motions of journal centers at bearing 1 and 2 stations: (a) Trajectories of journal centers; (b) trajectory of journal center at bearing 1 station; (c) trajectory of journal center at bearing 2 station; (d) Poincaré map of journal center at bearing 1 station; (e) Poincaré map of journal center at bearing 2 station; and (f) spectrum diagram of journal center.

Figure 17.

When e _x =e _y =9×10⁻⁷ m, the quasi-periodic motions of journal centers at bearing 1 and 2 stations: (a) Trajectories of journal centers; (b) trajectory of journal center; (c) Poincaré map of journal center; and (d) spectrum diagram of journal center.

Figure 18.

When e _x =e _y =1×10⁻⁶ m, the periodic motions of journal centers at bearing 1 and 2 stations: (a) Trajectories of journal centers; (b) trajectory of journal center; and (c) spectrum diagram of journal center.

Stability of rotor system

By calculating the subdeterminant of the coefficient matrix of determinant equation, the stability criterions can be obtained. The stability criterions Δ₁∼Δ₄ are shown in Figure 19(a) and the stability criterions Δ₅∼Δ₈ are shown in Figure 19(b). The stability criterions Δ₁∼Δ₄ are greater than 0 for ω=500–5000 rad/s, and the stability criterions Δ₅∼Δ₈ are greater than 0 for ω=500–3700 rad/s. When ω=3800–5000 rad/s, it can be seen in Figure 19(b), the stability criterions Δ₅∼Δ₈ are smaller than 0, and the rotor system shows instability.

Figure 19.

When ω = 500–5000 rad/s, the stability criterions: (a) stability criterions Δ₁∼Δ₄ and (b) stability criterions Δ₅∼Δ₈.

Conclusion

The nonlinear dynamic responses and bifurcation of asymmetric rotor system supported in two axial-grooved gas-lubricated bearings are investigated. A hybrid numerical model is established based on the rational function approximation and motion equation of rotor. The dynamic responses are obtained by solving the proposed model versus different rotating speed, and the computing cost is reduced more than the Newmark-β method. The bifurcation characteristics of gas-lubricated bearing–rotor system are analyzed with the change of the rotating speed and mass eccentricity.

1. The complex number of dynamic stiffness and damping coefficients of the gas-lubricated bearing is approximated as a rational function, which contains a set of rational function coefficients. The rational function coefficients are calculated by the vector fitting method. The rational function model of gas-lubricated bearing is established by using rational function coefficients, which shows the response relationship between nonlinear gas film forces of bearing and motion trajectory of rotor.

2. A hybrid numerical model is established by coupling the rational function model of gas-lubricated bearing with the motion equation of rotor, and the nonlinear dynamic responses of gas-lubricated bearing–rotor system are investigated by solving the proposed hybrid model. The hybrid numerical model avoids the multiple calculations of unsteady Reynolds equation in the dynamics analysis process, and the nonlinear gas film forces can be obtained simultaneously when the responses are solved. The computing cost can be saved.

3. The nonlinear dynamic responses of rotor system supported in gas-lubricated bearings are investigated by trajectory diagrams, Poincare maps, spectrum diagrams, and time series diagrams of the journal centers at bearing 1 and 2 stations. The bifurcation of the nonlinear dynamic response is studied by taking the rotating speed and mass eccentricity as bifurcation parameters. The results show that the unbalance responses of the rotor system have complex dynamic characteristics under different rotating speeds and mass eccentricities.

4. By the Routh–Hurwitz method, the stability of rotor system can be determined by the subdeterminant of the coefficient matrix of determinant equation, and the rotor system shows instability for 3800–5000 rad/s.