Dynamics Investigation on Axial-Groove Gas Bearing-Rotor System with Rod-Fastened Structure

Abstract

This research report discusses the dynamic behaviors of an axial-groove gas bearings-rotor system with rod-fastened structure. The time-based dependency-compressible Reynolds equation in the gas bearing nonlinear system is solved by the differential transformation method, and the continuous gas film forces of a three-axial-groove gas bearing are obtained. A dynamic mathematical model of the rotor system with rod-fastened structure supported in two- and three-axial-groove gas bearings with eight degrees of freedom is established. The dynamic motion equation of the rod-fastened rotor system is solved by the modified Newmark-β method based on disturbance compensation, which can reduce the computing error and improve computing stability. The dynamic characteristics of the rod-fastened rotor-gas bearing system are analyzed efficiently by the diversiform unbalance responses. The influence of the position angle of the pad on the nonlinear characteristics of the rod-fastened rotor system is also studied.

1. Introduction

Oil-free lubrication provides many unique characteristics for gas bearings, such as low friction, high positioning precision and low pollutant emission. As a result, the rotor system supported in a gas bearing is frequently employed within the high-speed rotating machinery, whether in civilian, agriculture, or military aerospace fields. The characteristic of rotor dynamics is one key factor for the stable operation of high-speed rotating machinery. It is necessary to carry out more research about the stability analysis and dynamic performance of the gas bearing-rotor system [1,2,3,4].

Li et al. [5] studied the effects of the circumferential and axial surface waviness on the nonlinear dynamic characteristics of a gas-lubricated bearing-rotor system; the dynamic stability of the bearing-rotor system can be improved with the increase of amplitude of circumferential waviness. Based on the finite difference method, perturbation method, and mixing method, Wang et al. [6] investigated the nonlinear dynamic performance of the flexible rotor system with opposed high-speed gas bearing support; by changing the parameters of the bearing, chaos phenomena and the system loss produced by the irregular vibration can be reduced. Yang et al. [7] analyzed numerically the unbalance response and dynamic stability of the rotor system supported on a micro gas bearing with a gas rarefaction effect; the influence of mass eccentricities of rotor on the stability of rotor system supported on the micro gas bearing was discussed. Wang et al. [8] analyzed the nonlinear dynamic behavior of a hybrid active aerostatic and aerodynamic bearing-rotor system by using a mixed numerical method, and the complex dynamic phenomena under different parameters and conditions were obtained. These gas bearings are cylindrical bearings, which have a simple structure and high-load-carrying capacity, but the rotor system supported by the cylindrical gas bearing is prone to instability under higher speed. The application of grooved gas bearings ensures the stable operation and improves the work performance of high-speed machinery; many research results were obtained on the lubrication performance of the grooved bearings. Feng et al. [9] focused on gas bearings with a micro spherical spiral groove, and the thermohydrodynamic performance of the bearings was analyzed by considering the surface roughness and air rarefaction, the effects of surface roughness and slip flow were obtained by using the Weierstrass–Mandelbrot function and Wu’s slip model. Zhang et al. [10] studied the contact characteristics of a gas-lubricated thrust micro-bearing with a spiral-groove by considering the gas rarefaction, and the effects of the standard deviation of asperity height and the groove depth on the surface contact forces of the gas-lubricated bearing were discussed. Jia et al. [11,12] built a mathematical model of spherical spiral-groove hybrid-gas bearings rotor system coupled with the analysis of nonlinear dynamic lubrication. The critical speed of the rotor-bearing system was calculated and the stability of the hybrid gas bearings was predicted; they then forced the dynamic variation rule of the gas film damping and stiffness under different states of motion, and the vortex motion and oscillation phenomenon of the gas film was analyzed. Liu et al. [13] studied the nonlinear dynamic characteristics of a rigid rotor system with herringbone-grooved-journal gas bearings support, the onset speed and the whirl frequency ratio of subsynchronous vibration in the gas bearings-rotor system were estimated and analyzed. Du et al. [14] investigated the whirl motion of a rigid rotor system with spiral-grooved opposed-hemisphere gas bearing support; the complicated rotor responses under the initial disturbance were analyzed, and the synchronous and nonsynchronous excitations were obtained. Based on time delays and feedback control gains, the dynamic behavior of the rotor supported by the self-acting three-axial-grooved gas-lubricated bearings were investigated by us [15], and a mathematical method used to obtain the dynamic characteristics of the nonlinear gas bearing-rotor system is presented. The unbalanced responses of a rigid rotor system supported by two-axial-groove gas bearings were analyzed by us [16]; a numerical model was established and utilized to acquire the rich nonlinear phenomena.

Continuous rotor models have been investigated thoroughly in rotor-gas bearing systems [5,6,7,8,11,12,13,14,15,16]. Li et al. [17] analyzed the rotor vibration characteristics of a gas bearing-rotor system, jointed by bolted-disk, by considering dynamic parameters of bolted-disk joints. Using the changes of the bending stiffness and tangential stiffness of the bolted-disk joints, the dynamic vibration characteristics of the rotor system were analyzed. An improved integrated modal technique with a free interface for reducing the number of degrees of freedom of a continuous flexible rotor system was presented by us [18]. The nonlinear dynamic behavior of a continuous rotor with multi-span bearing support was investigated by our proposed method. Chasalevris and Papadopoulos [19] proposed a new semi-analytical method for simulating the dynamic response of multi-segment continuous rotor bearing systems. Yang et al. [20] numerically suppressed the unbalance responses of the continuous flexible rotor systems supported on the tilting-pad gas bearings by using the proportional-derivative control model. Compared with the continuous rotors, the rod-fastened rotors are widely used to adapt the actual working conditions in the aircraft engine and gas turbines. Because of the effects of contacts among the components, the dynamic characteristics of the rod-fastened rotor have been studied increasingly. Liu et al. [21] proposed a method for investigating the dynamic stability and bifurcation behaviors of the rod fastening rotor and complete bearing system. The comparison results of the periodic motions and stability conditions in the rod fastening and continuous rotor-bearings system indicated that their bifurcation characteristics have a general resemblance. Wu et al. [22] built a mechanical model of a rod fastening rotor-bearing system by considering the contact effect between the disks. To characterize the contact interface, a nonlinear stiffness matrix containing stiffness coefficients was presented. Li et al. [23] studied the dynamic characteristics of an air bearing-rod fastening rotor system, which considered the heterogeneous bending stiffness of the contact interface; the linear and nonlinear dynamics behaviors were obtained. Wang et al. [24] analyzed the effects of the internal damping of a rod-fastened rotor-bearing nonlinear system on the dynamic responses and rotor stability. With the effect of internal damping, the amplitudes of dynamic response were decreased to some extent at lower speed and amplified significantly at higher speed. The dynamic responses and bifurcation phenomenon of a rod-fastened rotor supported by a fixed-tilting pad gas bearing were investigated by us [25]. The periodic orbits, pendulum angles, and Poincaré map with the different pad pivot ratios and preloads for the symmetrical and asymmetrical rotor-bearing system were analyzed. Our results showed that the rod-fastened rigid rotor system has better stability than the continuous integral rotor system. Hu et al. [26,27] investigated dynamic characteristics of the rod-fastened rotor bearing system, considering the rub-impact effect and influence of rotating velocities and radial stiffnesses of stator. Based on the D’Alembert principle, a nonlinear dynamic hybrid model of the rod-fastened rotor was built by considering different rub-impact and initial permanent deflection; the effects of rotating speed and the length of initial permanent deflection on the dynamic behaviors of the journal bearing rod-fastened rotor system were analyzed.

A dynamic mathematical model of a rotor system with a rod-fastened structure supported on three-axial-groove gas bearings is formulated in this report. The time-based dependency-compressible Reynolds equation in the gas bearing lubrication nonlinear system is solved by using the differential transformation method, and the nonlinear gas film forces of a three-axial-groove gas bearing are obtained. The eight degrees of freedom dynamic motion equation of the rod-fastened rotor-gas bearing system is solved by the modified Newmark-β method, which considers the disturbance compensation. The dynamic behaviors of the rod-fastened rotor-bearing nonlinear system are analyzed efficiently, and the dynamic unbalance responses of the rod-fastened rotor and continuous integral rotor, three-axial-grooved and two-axial-grooved gas bearings are compared. Furthermore, the effects of the position angle of the pad on the nonlinear dynamic behavior are also investigated.

2. Kinematic Model and Mathematical Method of Three-Axial-Groove Gas Bearing

2.1. Kinematic Model

The axial-groove gas bearing and the coordinate system are shown in Figure 1. Axial grooves are distributed uniformly in the circumferential direction. xO_by denotes the bearing coordinate system; O_b denotes the bearing center of the axial-groove gas bearing; O_j denotes the journal center; ϕ denotes the angle which is from the negative y direction to the gas film position along the rotating direction; e denotes the eccentricity of the gas bearing; ω denotes the rotation angular speed; R_b denotes the gas bearing radius; c denotes the radial clearance; α denotes the angle of groove width; α₀ denotes the angle from leading edge of the pad to the negative y direction; β denotes the arc bushing angle of the single pad; and F_x and F_y denote the gas film forces of the coordinate system xO_by in the x and y directions.

In order to obtain the gas film forces F_x, F_y of the axial-groove gas bearing, the gas pressure p needs to be calculated by solving the compressible Reynolds equation. Reynolds equation is the theoretical basis of fluid lubricated bearing research, which depicts the relationship between the gas pressure and gas thickness. Considering a gas journal bearing, the compressible lubrication Reynolds equation is:

$${\partial }_x\left(p{h}^3{\partial }_{x}p\right)+{\partial }_z\left(p{h}^3{\partial }_{z}p\right)=6\mu U{\partial }_x(ph)+12\mu {\partial }_t(ph)$$

(1)

where p denotes the gas film pressure in the coordinate system; h denotes the gas film thickness, for the ith pad, h_i = 1 + e_icos(ϕ − θ_i); z denotes the coordinate in the axial direction; x denotes the angle from the line of centers O_j and O_b to the gas film current position in the circumferential direction; μ denotes the average velocity of the gas; and U denotes the component of velocity in the x direction.

The dimensionless procedure is applied to Equation (1) by using the following transformational relations: p = P × Pa, h = H × c, z = λ × R_b, x = φ × R_b, τ = t × ω, U = ω × R. A unified equation can be then obtained as follows:

$${\partial }_{\phi }\left(P{H}^3{\partial }_{\phi }P\right)+{\partial }_{\lambda }\left(P{H}^3{\partial }_{\lambda }P\right)=\mathsf{\Lambda }{\partial }_{\phi }(PH)+2\mathsf{\Lambda }{\partial }_{\tau }(PH)$$

(2)

where Λ = 6μωR²/(Pac²) and denotes the bearing number of the gas bearing.

For Equation (2), the boundary conditions are:

• The gas film pressures at the both ends of bearing: P (φ, B_b/2R_b) = P (φ, −B_b/2R_b) = 1;
• The gas film pressures at the leading edge (φ_i1) and trailing edge (φ_i2) of single pad in the axial direction: P (φ₁, λ) = P (φ₂, λ) = 1.

By using the differential transformation method, the dimensionless Reynolds equation is solved efficiently.

2.2. Mathematical Method of Nonlinear Gas Film Force

The differential transformation method [28,29] has significant merits; rapid convergence can be achieved and the computational error can be reduced. The time is discretized in the nonlinear Reynolds equation by the differential transformation method. The gas film pressure P is a time-dependent function in the time domain, and the differential transform of the kth derivative of the gas film pressure P(t) can be expressed as:

$$P_t(k)=\frac{T^k}{k!}{\left[\frac{{\partial }^{k}P(t)}{\partial t^k}\right]}_{t=0}$$

(3)

where P_t(k) denotes the transform function in the transform domain. The differential inverse transform P_t(k) can be written as:

$$P(t)={\displaystyle \sum _{k=0}^m{P}_t(k){\left(\frac{t}{T}\right)}^k},0\le t\le T$$

(4)

where T denotes the time gap. While the value P(0) (k = 0) is known, P_t (0) can be obtained by Equation (4), and then the other discrete value P(k) can be obtained in the time interval T.

By using the differential transformation method, the compressible Reynolds equation of the axial-groove gas bearing can be solved, and the transform with respect to the time domain τ is made.

Introducing $S(k)=P\otimes P$, $I(k)=H\otimes H$, and $J(k)=H\otimes H\otimes H$, S(k) denotes the differential transform of P². I(k) and J(k) denote the differential transforms of H² and H³. The differential transformation of Equation (2) can be written as:

$$\begin{array}{l}3I\otimes {\partial }_{\phi }H\otimes {\partial }_{\phi }S+J\otimes {\partial }_{\phi }^{2}S+3I\otimes {\partial }_{\lambda }H\otimes {\partial }_{\lambda }S+J\otimes {\partial }_{\lambda }^{2}S=\\ 2\mathsf{\Lambda }{\partial }_{\phi }H\otimes P+2\mathsf{\Lambda }{\partial }_{\phi }P\otimes H+4\mathsf{\Lambda }{\partial }_{\tau }H\otimes P+4\mathsf{\Lambda }{\partial }_{\tau }P\otimes H\end{array}$$

(5)

where $\otimes$ denotes the convolution operator.

$$\mathrm{If}\{\begin{array}{c}y(k)=a{f}_1(t)\\ y(k)=f_1(t)f_2(t)\\ y(k)=f_1(t)f_2(t)f_3(t)\end{array},\mathrm{then}\{\begin{array}{c}Y(k)=a{\displaystyle \sum _{m=0}^k{F}_1(m)}\\ Y(k)={\displaystyle \sum _{m=0}^k{F}_1(m)F_2(k-m)}\\ Y(k)={\displaystyle \sum _{m=0}^k{F}_1(k-m){\displaystyle \sum _{n=0}^m{F}_2(n)F_3(m-n)}}\end{array}$$

Equation (5) is discretized in the circumferential and axial directions by the central difference method. The discretized equation can be obtained:

$$\begin{array}{l}3{\displaystyle \sum _{m=0}^k{I}_{i,j}(k-m)}\times {\displaystyle \sum _{n=0}^m\left[\left(\frac{S_{i+1,j}(n)-S_{i-1,j}(n)}{2\Delta \phi }\right)\times (\frac{H_{i+1,j}(m-n)-H_{i-1,j}(m-n)}{2\Delta \phi })\right]}+\\ {\displaystyle \sum _{m=0}^k\left[\left(\frac{S_{i+1,j}(m)-2S_{i,j}(m)+S_{i-1,j}(m)}{{(\Delta \phi )}^2}\right)J_{i,j}(k-m)\right]}+\\ 3{\displaystyle \sum _{m=0}^k{I}_{i,j}(k-m)}{\displaystyle \sum _{n=0}^m[\left(\frac{S_{i,j+1}(n)-S_{i,j-1}(n)}{2\Delta \lambda }\right)\times \left(\frac{H_{i,j+1}(m-n)-H_{i,j-1}(m-n)}{2\Delta \lambda }\right)]}\\ +{\displaystyle \sum _{m=0}^k\left[\left(\frac{S_{i,j+1}(m)-2S_{i,j}(m)+S_{i,j-1}(m)}{{(\Delta \lambda )}^2}\right)J_{i,j}(k-m)\right]}=2\mathsf{\Lambda }{\displaystyle \sum _{m=0}^k\left[P_{i,j}(m)\left(\frac{H_{i+1,j}(k-m)-H_{i-1,j}(k-m)}{2\Delta \phi }\right)\right]}+\\ 2\mathsf{\Lambda }{\displaystyle \sum _{m=0}^k\left[H_{i,j}(m)\left(\frac{P_{i+1,j}(k-m)-P_{i-1,j}(k-m)}{2\Delta \phi }\right)\right]}+4\mathsf{\Lambda }{\displaystyle \sum _{m=0}^k\left[\left(\frac{m+1}{\Delta \tau }\right)H_{i,j}(k-m)P_{i,j}(m+1)\right]}+\\ 4\mathsf{\Lambda }{\displaystyle \sum _{m=0}^k\left[\left(\frac{m+1}{\Delta \tau }\right)P_{i,j}(k-m)H_{i,j}(m+1)\right]}\end{array}$$

(6)

where Δφ and Δλ denote the steps in the circumferential and axial directions; Δτ denotes the time step; k, m and n denote the differential transformation orders (k ≥ m ≥ n); and i and j denote the coordinates of the nodal position in the φ and λ directions.

Because the centerline of the journal is parallel to the centerline of the gas bearing, the gas film thickness is constant in the λ direction, i.e., ∂_λH_i,j = 0. When k = m = n = 0, Equation (6) can be transformed as:

$$\begin{array}{l}3I_{i,j}(0)\left[\frac{S_{i+1,j}(0)-S_{i-1,j}(0)}{2\Delta \phi }\frac{H_{i+1}(0)-H_{i-1}(0)}{2\Delta \phi }\right]+\left(\frac{S_{i+1,j}(0)-2S_{i,j}(0)+S_{i-1,j}(0)}{{(\Delta \phi )}^2}\right)J_{i,j}(0)+\\ \left(\frac{S_{i,j+1}(0)-2S_{i,j}(0)+S_{i,j-1}(0)}{{(\Delta \lambda )}^2}\right)J_{i,j}(0)=2\mathsf{\Lambda }{P}_{i,j}(0)\left(\frac{H_{i+1}(0)-H_{i-1}(0)}{2\Delta \phi }\right)+\\ 2\mathsf{\Lambda }{H}_i(0)\left(\frac{P_{i+1,j}(0)-P_{i-1,j}(0)}{2\Delta \phi }\right)+4\mathsf{\Lambda }\left(\frac{1}{\Delta \tau }\right)H_i(0)P_{i,j}(1)+4\mathsf{\Lambda }\left(\frac{1}{\Delta \tau }\right)P_{i,j}(0)H_i(1)\end{array}$$

(7)

S_i,j(0), I_i,j (0), and J_i,j (0) are expressed as:

$$S_{i,j}(0)=P_{i,j}(0)P_{i,j}(0),I_{i,j}(0)=H_i(0)H_i(0),J_{i,j}(0)=H_i(0)H_i(0)H_i(0)$$

(8)

where H_i(0) denotes the gas film thickness in the last moment and P_i,j(0) denotes the gas film pressure in the last moment.

By substituting Equation (8) into Equation (7), Equation (9) follows:

$$\begin{array}{l}\frac{3H_i^2(0)\left(H_{i+1}(0)-H_{i-1}(0)\right)+4H_i^3(0)}{4{(\Delta \phi )}^2}{S}_{i+1,j}(0)+\frac{-3H_i^2(0)\left(H_{i+1}(0)-H_{i-1}(0)\right)+4H_i^3(0)}{4{(\Delta \phi )}^2}{S}_{i-1,j}(0)+\\ \left(\frac{H_i^3(0)}{{(\Delta \lambda )}^2}\right)S_{i,j+1}(0)+\left(\frac{H_i^3(0)}{{(\Delta \lambda )}^2}\right)S_{i,j-1}(0)+\left(-\frac{2H_i^3(0)}{{(\Delta \phi )}^2}-\frac{2H_i^3(0)}{{(\Delta \lambda )}^2}\right)S_{i,j}(0)=\\ \mathsf{\Lambda }\left(\frac{H_i(0)}{\Delta \phi }\right)P_{i+1,j}(0)+\mathsf{\Lambda }\left(-\frac{H_i(0)}{\Delta \phi }\right)P_{i-1,j}(0)+\mathsf{\Lambda }\left(\frac{H_{i+1}(0)-H_{i-1}(0)}{\Delta \phi }\right)P_{i,j}(0)\\ +4\mathsf{\Lambda }\left(\frac{1}{\Delta \tau }\right)H_i(0)P_{i,j}(1)+4\mathsf{\Lambda }\left(\frac{1}{\Delta \tau }\right)P_{i,j}(0)H_i(1)\end{array}$$

(9)

where H_i(1) and P_i,j(1) denote the relevant variables of the thickness and pressure of the gas film in the time interval τ.

When k = 1, then m = n = 0 and Equation (6) can be written as:

$$\begin{array}{l}3I_{i,j}(1)\left[\frac{S_{i+1,j}(0)-S_{i-1,j}(0)}{2\Delta \phi }\frac{H_{i+1}(0)-H_{i-1}(0)}{2\Delta \phi }\right]+3I_{i,j}(0)\left[\frac{S_{i+1,j}(0)-S_{i-1,j}(0)}{2\Delta \phi }\frac{H_{i+1}(1)-H_{i-1}(1)}{2\Delta \phi }\right]+\\ 3I_{i,j}(0)\left[\frac{S_{i+1,j}(1)-S_{i-1,j}(1)}{2\Delta \phi }\frac{H_{i+1}(0)-H_{i-1}(0)}{2\Delta \phi }\right]+J_{i,j}(1)\left(\frac{S_{i+1,j}(0)-2S_{i,j}(0)+S_{i-1,j}(0)}{{(\Delta \phi )}^2}\right)+\\ J_{i,j}(0)\left(\frac{S_{i+1,j}(1)-2S_{i,j}(1)+S_{i-1,j}(1)}{{(\Delta \phi )}^2}\right)+J_{i,j}(1)\left(\frac{S_{i,j+1}(0)-2S_{i,j}(0)+S_{i,j-1}(0)}{{(\Delta \lambda )}^2}\right)+\\ J_{i,j}(0)\left(\frac{S_{i,j+1}(1)-2S_{i,j}(1)+S_{i,j-1}(1)}{{(\Delta \lambda )}^2}\right)=2\mathsf{\Lambda }\left(\frac{H_{i+1}(0)-H_{i-1}(0)}{2\Delta \phi }\right)P_{i,j}(1)+\\ 2\mathsf{\Lambda }\left(\frac{H_{i+1}(1)-H_{i-1}(1)}{2\Delta \phi }\right)P_{i,j}(0)+2\mathsf{\Lambda }\left(\frac{P_{i+1,j}(0)-P_{i-1,j}(0)}{2\Delta \phi }\right)H_i(1)+\\ 2\mathsf{\Lambda }\left(\frac{P_{i+1,j}(1)-P_{i-1,j}(1)}{2\Delta \phi }\right)H_i(0)+4\mathsf{\Lambda }\left(\frac{1}{\Delta \tau }\right)P_{i,j}(1)H_i(1)+4\mathsf{\Lambda }\left(\frac{2}{\Delta \tau }\right)P_{i,j}(0)H_i(2)+\\ 4\mathsf{\Lambda }\left(\frac{1}{\Delta \tau }\right)H_i(1)P_{i,j}(1)+4\mathsf{\Lambda }\left(\frac{2}{\Delta \tau }\right)H_i(0)P_{i,j}(2)\end{array}$$

(10)

S_i,j(1), I_i,j (1), and J_i,j (1) can be expressed as:

$$S_{i,j}(1)=2P_{i,j}(0)P_{i,j}(1),I_{i,j}(1)=2H_i(0)H_i(1),J_{i,j}(1)=3H_i(0)H_i(0)H_i(1)$$

(11)

By substituting Equations (8) and (11) into Equation (10), one can obtain:

$$\begin{array}{l}\frac{S_{i+1,j}(0)}{4{(\Delta \phi )}^2}[6H_i(0)H_i(1)\left(H_{i+1}(0)-H_{i-1}(0)\right)+3H_i^2(0)\left(H_{i+1}(1)-H_{i-1}(1)\right)+12H_i^2(0)H_i(1)]+\\ \frac{S_{i-1,j}(0)}{4{(\Delta \phi )}^2}[-6H_i(0)H_i(1)\left(H_{i+1}(0)-H_{i-1}(0)\right)-3H_i^2(0)\left(H_{i+1}(1)-H_{i-1}(1)\right)+12H_i^2(0)H_i(1)]+\\ \frac{3H_i^2(0)H_i(1)}{{(\Delta \lambda )}^2}{S}_{i,j+1}(0)+\frac{3H_i^2(0)H_i(1)}{{(\Delta \lambda )}^2}{S}_{i,j-1}(0)+\left(-\frac{6H_i^2(0)H_i(1)}{{(\Delta \phi )}^2}-\frac{6H_i^2(0)H_i(1)}{{(\Delta \lambda )}^2}\right)S_{i,j}(0)+\\ \frac{3H_i^2(0)\left(H_{i+1}(0)-H_{i-1}(0)\right)+4H_i^3(0)}{4{(\Delta \phi )}^2}{S}_{i+1,j}(1)+\frac{-3H_i^2(0)\left(H_{i+1}(0)-H_{i-1}(0)\right)+4H_i^3(0)}{4{(\Delta \phi )}^2}{S}_{i-1,j}(1)\\ +\frac{H_i^3(0)}{{(\Delta \lambda )}^2}{S}_{i,j+1}(1)+\frac{H_i^3(0)}{{(\Delta \lambda )}^2}{S}_{i,j-1}(1)+\left(-\frac{2H_i^3(0)}{{(\Delta \phi )}^2}-\frac{2H_i^3(0)}{{(\Delta \lambda )}^2}\right)S_{i,j}(1)=\\ [\frac{\mathsf{\Lambda }{H}_i(1)}{\Delta \phi }{P}_{i+1,j}(0)+\frac{-\mathsf{\Lambda }{H}_i(1)}{\Delta \phi }{P}_{i-1,j}(0)+\left(\frac{\mathsf{\Lambda }\left(H_{i+1}(1)-H_{i-1}(1)\right)}{\Delta \phi }+\frac{8\mathsf{\Lambda }\left(H_i(2)\right)}{\Delta \tau }\right)P_{i,j}(0)]+\\ [\frac{\mathsf{\Lambda }{H}_i(0)}{\Delta \phi }{P}_{i+1,j}(1)+\frac{-\mathsf{\Lambda }{H}_i(0)}{\Delta \phi }{P}_{i-1,j}(1)+\left(\frac{\mathsf{\Lambda }\left(H_{i+1}(0)-H_{i-1}(0)\right)}{\Delta \phi }+\frac{8\mathsf{\Lambda }\left(H_i(1)\right)}{\Delta \tau }\right)P_{i,j}(1)]+\\ 4\mathsf{\Lambda }(\frac{2}{\Delta \tau })H_i(0)P_{i,j}(2)\end{array}$$

(12)

By substituting H_i(0), P_i,j(0), and Hi(1) into Equations (9) and (12), P_i,j (1) can be calculated from Equation (9) and P_i,j (2) can be calculated from Equation (12). P_i,j is calculated in the following form:

$$P_{i,j}={\displaystyle \sum _{k=0}^2{P}_{i,j}(k)}$$

(13)

The calculation flow chart of the differential transform method is shown in Figure 2.

3. Dynamic Model and Method of the Rod-Fastened Rotor-Gas Bearing System

3.1. Dynamic Motion of the Rod-Fastened Rotor-Gas Bearing System

The rod-fastened rotor-gas bearing system is depicted in Figure 3. The rod-fastened rotor consists of shafts, disks, and gas bearings. Two rigid disks are bolted by rods, which are circumferentially distributed on the disk. The mass and damping of the rod are ignored when it is equivalent to a nonlinear stiffness spring structure. The rod-fastened rotor is supported by two- and three-axial-groove gas bearings, and the mass of the axial-groove gas bearing is not considered. In Figure 3a, O₁ and O₂ denote the journal centers at the left- and right-bearing stations, respectively; m₁ and m₂ denote the lumped masses of the left and right segments, respectively; l₁ and l₂ denote the lengths of the left and right segments, respectively; O_d1 and O_d2 denote the centers of the rigid disks, respectively; and m_d1 and m_d2 denote the masses of two rigid disks, respectively. Figure 3b shows the circumferential distribution of the rods; R_d denotes the radius of rigid disks.

The dynamic motion model of the rod-fastened rotor system with two- and three-axial-groove gas bearings support follows:

$$M\ddot{u}+Ku=P$$

(14)

where u = [x₁, y₁, x₂, y₂, x_d1, y_d1, x_d2, y_d2]^T; M denotes the mass matrix of the dynamic motion; K denotes the stiffness matrix; and P denotes the dynamic resultant force matrix.

$$M=\left[\begin{array}{cccccccc}{m}_1& 0& 0& 0& 0& 0& 0& 0\\ 0& m_1& 0& 0& 0& 0& 0& 0\\ 0& 0& m_2& 0& 0& 0& 0& 0\\ 0& 0& 0& m_2& 0& 0& 0& 0\\ 0& 0& 0& 0& m_{d1}& 0& 0& 0\\ 0& 0& 0& 0& 0& m_{d1}& 0& 0\\ 0& 0& 0& 0& 0& 0& m_{d2}& 0\\ 0& 0& 0& 0& 0& 0& 0& m_{d2}\end{array}\right]$$

(15)

$$K=\left[\begin{array}{cccccccc}{k}_1& 0& 0& 0& -k_1& 0& 0& 0\\ 0& k_1& 0& 0& 0& -k_1& 0& 0\\ 0& 0& k_2& 0& 0& 0& -k_2& 0\\ 0& 0& 0& k_2& 0& 0& 0& -k_2\\ -k_1& 0& 0& 0& k_1-k_b& 0& k_b& 0\\ 0& -k_1& 0& 0& 0& k_1-k_b& 0& k_b\\ 0& 0& -k_2& 0& -k_b& 0& k_2+k_b& 0\\ 0& 0& 0& -k_2& 0& -k_b& 0& k_2+k_b\end{array}\right]$$

(16)

$$P=\left[\begin{array}{c}-f_{x1}\\ -f_{y1}+m_{1}g\\ -f_{x2}\\ -f_{y2}+m_{2}g\\ f_{rx1}+m_{d1}{e}_{d1}{\omega }^2\cos \left(\omega t\right)+m_{d1}{e}_{d1}{\omega }^2\sin \left(\omega t\right)\\ f_{ry1}+m_{d1}{e}_{d1}{\omega }^2\cos \left(\omega t\right)-m_{d1}{e}_{d1}{\omega }^2\sin \left(\omega t\right)+m_{d1}g\\ f_{rx2}+m_{d2}{e}_{d2}{\omega }^2\cos \left(\omega t\right)+m_{d2}{e}_{d2}{\omega }^2\sin \left(\omega t\right)\\ f_{ry2}+m_{d2}{e}_{d2}{\omega }^2\cos \left(\omega t\right)-m_{d2}{e}_{d2}{\omega }^2\sin \left(\omega t\right)+m_{d2}g\end{array}\right]$$

(17)

where f_x and f_y denote the gas film forces of the axial-groove gas bearing, respectively; e_d1 and e_d1 are the mass eccentricities of the rigid disks, respectively; x_d1 and y_d1 denote the displacements of disk 1 in the x and y directions, respectively; x_d2 and y_d2 denote the displacements of disk 2 in the x and y directions, respectively; k₁ and k₂ denote the bending stiffnesses of the left and right segments, respectively; k_b denotes the bending stiffness of the bolted rod; f_rx1 and f_ry1 denote the restoring forces of the disk 1, respectively; f_rx2 and f_ry2 denote the restoring forces of the disk 2, respectively; and k_r denotes the restoring stiffness.

$$\{\begin{array}{c}{f}_{rx1}=k_r{\left(x_{d2}-x_{d1}\right)}^3\\ f_{ry1}=k_r{\left(y_{d2}-y_{d1}\right)}^3\end{array},\{\begin{array}{c}{f}_{rx2}=-k_r{\left(x_{d2}-x_{d1}\right)}^3\\ f_{ry2}=-k_r{\left(y_{d2}-y_{d1}\right)}^3\end{array}$$

(18)

where x_r1, y_r1, x_r1, and y_r1 denote the motion displacements of disk centers in the x and y directions in the last moment.

The dynamic motion model of the rod-fastened rotor system can be solved by using the Newmark-β method; however, calculation errors may occur while a disturbance is produced by an increment of acceleration.

3.2. Newmark-β Method Based on Disturbance Compensation

For solving the dynamic motion model of the rod-fastened rotor-gas bearing system efficiently, the modified Newmark-β method is employed by considering the disturbance compensation. The Newmark-β method is based on the average constant acceleration theory; that is, the acceleration is assumed to be a constant between ${\ddot{u}}_{(i-1)\tau }$ and ${\ddot{u}}_{i\tau }$. The integration constants of Newmark-β method are given in

Table 1. Integration constant parameters of Newmark-β method.

Parameters	Value or Formula
α	0.25
δ, δ = 2α	0.5
a0	1/ατ2
a1	δ/ατ2
a2	1/ατ
a3	1/2α − 1
a4	δ/α − 1
a5	(1/2α − 1)
a6	(1 − δ)τ
a7	δτ

.

For the equation of motion, the dynamic response equation at time step t = iτ is expressed as follows:

$$M{\ddot{u}}_{i\tau }+K{u}_{i\tau }=P_{i\tau }$$

(19)

where ${\ddot{u}}_{i\tau }$ and $u_{i\tau }$ denote the acceleration and displacement at time step iτ. By the average constant acceleration theory, the average constant acceleration $\ddot{u}$ can be obtained as $\ddot{u}={\ddot{u}}_{(i-1)\tau }+\delta ({\ddot{u}}_{i\tau }-{\ddot{u}}_{(i-1)\tau })$, which is between previous time step (i − 1)τ and time step iτ, ${\ddot{u}}_{(i-1)\tau }$ denotes the acceleration at previous time step (i − 1)τ. ${\dot{u}}_{i\tau }$ denotes the acceleration at time iτ; ${\dot{u}}_{i\tau }$ and $u_{i\tau }$ can be written as:

$$u_{i\tau }=u_{(i-1)\tau }+\tau {\dot{u}}_{(i-1)\tau }+\frac{1}{2}{\tau }^2\ddot{u}=u_{(i-1)\tau }+\tau {\dot{u}}_{(i-1)\tau }+\frac{{\tau }^2}{2}\left[(1-\delta ){\ddot{u}}_{(i-1)\tau }+\delta {\ddot{u}}_{i\tau }\right]$$

(20)

$${\ddot{u}}_{i\tau }={\dot{u}}_{(i-1)\tau }+\tau \ddot{u}={\dot{u}}_{(i-1)\tau }+\tau (1-\delta ){\ddot{u}}_{(i-1)\tau }+\tau \delta {\ddot{u}}_{i\tau }={\dot{u}}_{(i-1)\tau }+a_6{\ddot{u}}_{(i-1)\tau }+a_7{\ddot{u}}_{i\tau }$$

(21)

where $u_{(i-1)\tau }$ and ${\dot{u}}_{(i-1)\tau }$ are the velocity and displacement at the previous time step (i − 1)τ. According to Equation (20), ${\ddot{u}}_{i\tau }$ can be written as:

$${\ddot{u}}_{i\tau }=\frac{1}{\alpha {\tau }^2}(u_{i\tau }-u_{(i-1)\tau })-\frac{1}{\alpha \tau }{\dot{u}}_{(i-1)\tau }-(\frac{1}{2\alpha }-1){\ddot{u}}_{(i-1)\tau }=a_1(u_{i\tau }-u_{(i-1)\tau })-a_4{\dot{u}}_{(i-1)\tau }-a_5{\ddot{u}}_{(i-1)\tau }$$

(22)

By substituting Equations (21) and (22) into Equation (19), the equivalent motion Equation (23) can be obtained:

$$\widehat{K}{u}_{i\tau }={\widehat{P}}_{i\tau }$$

(23)

where $\widehat{K}$ denotes the equivalent stiffness matrix, which does not change with time; ${\widehat{P}}_{i\tau }$ denotes the equivalent dynamic resultant force matrix which is relevant to time.

$$\widehat{K}=K+a_{0}M$$

(24)

$${\widehat{P}}_{i\tau }=P_{i\tau }+M\left[a_0{u}_{(i-1)\tau }+a_2{\dot{u}}_{(i-1)\tau }+a_3{\ddot{u}}_{(i-1)\tau }\right]$$

(25)

The displacement $u_{i\tau }$ at time iτ can be obtained by solving Equation (23).

$$u_{i\tau }=(K+a_{0}M{)}^{-1}\{P_{i\tau }+M\left[a_0{u}_{(i-1)\tau }+a_2{\dot{u}}_{(i-1)\tau }+a_3{\ddot{u}}_{(i-1)\tau }\right]\}$$

(26)

According to the constrained condition of motion, the ${\ddot{u}}_{i\tau }$, ${\dot{u}}_{i\tau }$, and $u_{i\tau }$ are obtained by the Newmark-β method, but according to the dynamic motion equation, the actual acceleration ${\tilde{\ddot{u}}}_{i\tau }$ is obtained by substituting ${\dot{u}}_{i\tau }$ and $u_{i\tau }$ into Equation (19).

$${\tilde{\ddot{u}}}_{i\tau }=M^{-1}(P_{i\tau }-K{u}_{i\tau })$$

(27)

It can be seen that an increment of acceleration $\Delta {\tilde{\ddot{u}}}_{i\tau }$ at time step iτ is produced as Equation (28), which produces a disturbance error at each time step. By substituting ${\tilde{\ddot{u}}}_{i\tau }={\ddot{u}}_{i\tau }-\Delta {\tilde{\ddot{u}}}_{i\tau }$ into Equation (27), Equation (29) can be obtained; the increment of acceleration gives rise to an increment of dynamic resultant force $\Delta P_{i\tau }=M\Delta {\tilde{\ddot{u}}}_{i\tau }$, which is the disturbance increment produced by the increment of acceleration $\Delta {\tilde{\ddot{u}}}_{i\tau }$.

$$\Delta {\tilde{\ddot{u}}}_{i\tau }={\ddot{u}}_{i\tau }-{\tilde{\ddot{u}}}_{i\tau }$$

(28)

$$M{\ddot{u}}_{i\tau }+K{u}_{i\tau }=P_{i\tau }+\Delta P_{i\tau }$$

(29)

For acquiring the displacements, velocities, and accelerations at time step t = iτ, which satisfy the constraint condition of motion and dynamic motion equation simultaneously, the increment of dynamic resultant force needs to be eliminated by the following process. According to the motion equation, an increment motion equation is established; $\Delta {\ddot{u}}_{i\tau }$, $\Delta {\dot{u}}_{i\tau }$, and $\Delta u_{i\tau }$ are the displacement increment, velocity increment and acceleration increment caused by $-\Delta P_{i\tau }$, respectively, which can be obtained by Equation (30):

$$M\Delta {\ddot{u}}_{i\tau }+K\Delta u_{i\tau }=-\Delta P_{i\tau }$$

(30)

$$\{\begin{array}{c}\Delta u_{i\tau }=-{\widehat{K}}^{-1}\Delta P_{i\tau }\\ \Delta {\dot{u}}_{i\tau }=a_1\Delta u_{i\tau }\\ \Delta {\ddot{u}}_{i\tau }=a_0\Delta u_{i\tau }\end{array}$$

(31)

The displacement, velocity, and acceleration can be obtained by adding Equations (29) and (30), which are considering the incremental compensation.

$$\{\begin{array}{c}{u}_{i\tau }{}^{*}=u_{i\tau }+\Delta u_{i\tau }\\ {\dot{u}}_{i\tau }{}^{*}={\dot{u}}_{i\tau }+\Delta {\dot{u}}_{i\tau }\\ {\ddot{u}}_{i\tau }{}^{*}={\ddot{u}}_{i\tau }+\Delta {\ddot{u}}_{i\tau }\end{array}$$

(32)

The displacement $u_{i\tau }{}^{*}$, velocity ${\dot{u}}_{i\tau }{}^{*}$ and acceleration ${\ddot{u}}_{i\tau }{}^{*}$ are the accurate solution of the motion equation, which take the incremental compensation into consideration in Newmark-β method. The solutions conform to the dynamic motion equation, meanwhile, satisfy the motion constraint condition of the Newmark-β method.

4. Results and Discussions

The dynamics model of the rod-fastened rotor-gas bearing system is established in Section 3.1. The geometric construction of the rod-fastened rotor system is shown in Figure 3a. The dynamic responses of the rod-fastened rotor-gas bearing system are obtained using the modified Newmark-β method, illustrated in Section 3.2. In the following, the unbalance responses of the rod-fastened rotor system with two- and three-axial-groove gas bearings support are calculated and analyzed by orbit diagrams, Poincaré maps, and spectrograms. The parameter values selected for the rod-fastened rotor and the three-axial-groove gas bearing are listed in

Table 2. Structure parameter values of the rod-fastened rotor-gas bearing system.

Parameter	Value	Unit
Structure parameters parameter values of the three-axial-groove gas bearing
Radius of bearing, Rb	0.005	m
Bearing width, Bb	0.01	m
Width-to-diameter ratio, Bb/2Rb	1
Radial clearance, c	5 × 10−6	m
Position angle of pad, α0	10	deg
Groove width angle, α	5	deg
Arc bushing angle, β	115	deg
Structure parameters parameter values of the rod-fastened rotor
Radius of shaft, R	0.005	m
Length of the left and right segments, l1 = l2	0.2	m
Radius of disk, Rd	0.01	m
Width of disk, Bd	0.02	m
Mass eccentricities of the disk, ed1 = ed2	1 × 10−6	m

.

4.1. Comparison of the Rod-Fastened Rotor and Continuous Integral Rotor

The orbits of the rod-fastened rotor and continuous integral rotor are compared by using the same parameter values. When ω = 1500 rad/s, Figure 4a shows the orbits comparison of the rod-fastened rotor and continuous integral rotor at the left bearing station, and Figure 4b shows the spectra of the continuous integral rotor. The dynamic response is period-doubling motion in the continuous integral rotor system, whereas, the orbit of the rod-fastened rotor is periodic motion. For the model of the gas bearing-rotor system in this report, when ω = 1500 rad/s, the unbalance responses of the rod-fastened rotor are more stable than the continuous integral rotor.

4.2. Bifurcation Characteristics

The bifurcation diagram depicts the projections of orbit of the journal center in the y direction, i.e., the y₁ of the vector u = [x₁, y₁, x₂, y₂, x_d1, y_d1, x_d2, y_d2]^T in Equation (14). When the rotating speed is at ω = 900~2500 rad/s, the bifurcation behavior of the rod-fastened rotor at the left bearing station are shown in Figure 5. The rod-fastened rotor system experiences quasiperiodic motion at ω = 900~1050 rad/s, and inversely bifurcates to periodic motion at ω = 1100 rad/s, and then maintains the stable motion state until ω = 1750 rad/s. When the rotating speed increases to 1775 rad/s, the response of the rod-fastened rotor-gas bearing system bifurcates to quasiperiodic motion again, and the amplitude of quasiperiodic motion increases with the amplification of the rotating speed. From the bifurcation diagram, the rod-fastened rotor-gas bearing system maintains a relatively stable operation at ω = 1100~1750 rad/s.

4.3. Unbalance Responses Versus Rotating Speed ω

At different rotating speeds, the journal center of the rod-fastened rotor experiences different motion states. The parameters of the rod-fastened rotor and gas bearing are the same as those given in Table 2. In Figure 6, when ω = 1000 rad/s, the dynamic motion of the rod-fastened rotor comes out of quasiperiodic motion. Figure 6a depicts the orbit and Poincaré map of the journal center at the left bearing station, and Figure 6b depicts the orbit and Poincaré map of journal center at the left disk station. The spectrograms of the rod-fastened at the left bearing and left disk stations are described in Figure 6c. Figure 6d shows the comparison of the time series diagram at the left bearing and left disk stations. It can be seen that the amplitude at the left disk station is greater than the left bearing station, as shown in Figure 6c,d. When ω = 1500 rad/s, the periodic motion of the rod-fastened rotor-gas bearing system is described in Figure 7. Figure 7a,b shows the orbit and Poincaré map at the left bearing and left disk stations; the orbit is a single circle, and the synchronous behavior is inferred from the single projection point in the Poincaré map. Figure 7c,d shows the comparison of the spectrogram and time series diagram at the left bearing and left disk stations; the amplitude at the left disk station is greater but similar to the quasiperiodic motion.

When ω = 2000 rad/s, the quasiperiodic motion of the rod-fastened rotor-gas bearing system is shown in Figure 8a–d. It can be seen in Figure 8a,b, the orbits of the journal center at the left bearing and left disk stations are periodic multirings, but not the standard circle; the Poincaré map presents a closed loop on the projection plane, the quasiperiodic motion is more unstable than at ω = 1000 rad/s. The spectrogram and time series diagram at the left bearing and left disk stations are described in Figure 8c,d; the amplitude of the left bearing is slightly greater than the one at the left disk station. In Figure 9, the rotating speed is 2500 rad/s, and the unbalance response of the rod-fastened rotor-gas bearing system is still quasiperiodic motion. Figure 9a,b shows the orbit diagrams and Poincaré map at the left bearing and left disk stations, Figure 9c,d shows the spectrogram and time series diagram at the left bearing and left disk stations. The amplitudes of the journal center of the rod-fastened rotor gas bearing system are much greater, and the dynamic stability of the rod-fastened rotor-gas bearing system is obviously reduced.

4.4. Effects of the Position Angle of Pad α₀

The change of the position angle of pad α₀ influences the load-carry duty, and then effects the unbalance responses of the rod-fastened rotor-gas bearing system. When the position angle of pad α₀ = 10°, the load carrying is on pad 1 (indicated in Figure 1a), when the position angle of pad α₀ = 62.5°, load carrying is between the pad 1 and pad 2. For ω = 1500 rad/s, the comparisons of the orbits of the rod-fastened rotor at the left bearing and disk stations at α₀ = 10° and 62.5° are shown in Figure 10a,b; the dynamic motions are both periodic motions. The amplitudes of the orbits at the left bearing and disk stations are both smaller at α₀ = 62.5° than α₀ = 10°; the nonlinear dynamic motion of the rod-fastened rotor-gas bearing system is slightly stable when the load carrying is between the pads. When ω = 2000 rad/s, the unbalance responses of the rod-fastened rotor-gas bearing system are quasiperiodic motion and periodic motion at α₀ = 10° and 62.5°, respectively. In Figure 11, the comparisons of the orbits of the rod-fastened rotor at the left bearing and disk stations at α₀ = 10° and 62.5° show that the motion is more stable at α₀ = 62.5°. When ω is 2500 rad/s, the dynamic motion of the rod-fastened rotor-gas bearing system becomes quasiperiodic at α₀ = 62.5°, as indicated in Figure 12, but the vibration amplitudes are slightly smaller than the one at α₀ = 10°.

4.5. Comparison of Three-Axial-Grooved and Two-Axial-Grooved Gas Bearings

Figure 13 and Figure 14 show the orbit comparisons of the rod-fastened rotor-gas bearing system with three-axial-groove and two-axial-groove gas bearings at ω = 1500 and 2000 rad/s. When ω is 1500 rad/s, the comparison of the orbits of the rod-fastened rotor at the left bearing and disk stations is shown in Figure 13a,b; the orbits of journal center under the two conditions are both periodic motions, but the orbit amplitudes of the rod-fastened rotor system with two-axial-groove gas bearings are much smaller than the orbit amplitudes of the rod-fastened rotor system with three-axial-groove gas bearings. Figure 14a,b shows the comparisons of the orbits of the rod-fastened rotor at the left bearing and left disk stations at ω = 2000 rad/s, and the dynamic motion of the rod-fastened rotor system with three-axial-groove gas bearings is quasiperiodic motion, and the dynamic motion of the rod-fastened rotor system with two-axial-groove gas bearings is periodic motion.

5. Conclusions

Different from the cylindrical bearing and continuous integral rotor structure, in this report, a rod-fastened rotor and two- and three-axial-groove gas bearings are applied in the bearing-rotor system. The continuous gas film forces of the three-axial-groove gas bearing are obtained from the time-based dependency compressible Reynolds equation in the gas bearing lubrication nonlinear system, which is solved by the differential transformation method. The Newmark-β method is used and modified to compensate the calculation error by considering the disturbance compensation. The modified Newmark-β method conforms to the dynamic motion equation of the rod-fastened rotor-gas bearing system, while satisfying the motion constraint condition of the Newmark-β method.

A dynamic mathematical method of a rod-fastened rotor system with two- and three-axial-groove gas bearings support is established. The dynamic characteristics of the rod-fastened rotor-gas bearing system under different rotating speeds are analyzed by the orbits, Poincaré maps, and spectrograms. The orbits and Poincaré maps of the rod-fastened rotor and continuous integral rotor, three-axial-grooved and two-axial-grooved gas bearings are obtained and compared. In this report, the comparison shows that the rod-fastened rotor has better stability than the continuous integral rotor at ω = 1500 rad/s, and that the rotor system with the two-axial-grooved gas bearing has better stability. The bifurcation characteristics of the rod-fastened rotor-gas bearing system are investigated; the bifurcation behavior shows that there are different periodic movements, including periodic and quasiperiodic motion. The effects of the position angle of the pad on the orbits and Poincaré map of the journal center of rod-fastened system are studied. The vibration amplitude can be decreased when the position angle of the pad is α₀ = 62.5° versus α₀ = 10°.