# Chaotic and Subharmonic Motion Analysis of Floating Ring Gas Bearing System by Hybrid Numerical Method

Authors: Cheng-Chi Wang; Her-Terng Yau; Chi-Chang Wang
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This paper studies the nonlinear dynamic behaviors including chaotic, subharmonic, and quasi-periodic motions of a rigid rotor supported by floating ring gas bearing (FRGB) system. A hybrid numerical method combining the differential transformation method and the finite difference method used to calculate pressure distribution of FRGB system and rotor orbits. The results obtained for the orbits of the rotor center are in good agreement with those obtained using the traditional finite difference approach. Moreover, the hybrid method avoids the numerical instability problem suffered by the finite difference scheme at low values of the rotor mass and computational time-step. Moreover, power spectra, Poincaré maps, bifurcation diagrams and Lyapunov exponents are applied to examine the nonlinear dynamic response of the FRGB system over representative ranges of the rotor mass and bearing number, respectively. The results presented summarize the changes which take place in the dynamic behavior of the FRGB system as the rotor mass and bearing number are increased and therefore provide a useful guideline for the bearing system.


---

## Body

## 1.  Introduction

Gas bearing systems have been extensively used for a variety of electromechanical system applications, and it is particularly valuable when used with precision instruments. This is due, in part, to low noise when there is rotation, and to zero friction when the instruments are used as null devices. This bearing system has a number of advantages compared to their rolling-element or oil-lubricated counterparts, including low friction losses and zero risk of contamination through lubricant leakage. As a result, they are widely applied in a diverse range of rotational systems.Floating ring gas bearing (FRGB) system is different from general gas bearings due to the double thins lubrication. Because the rotational speed of the rotor applied in precision control could reach 106 rpm, the stability of rotor dynamics is one of the most important factors for considering the design of FRGB systems. So, how to increase the stability of rotor systems and control the appearance of nonperiodic motion becomes the major execution of this paper.According to the recent research, nonlinear dynamic responses of rotor-bearing systems are analyzed and published. In 1994, Malik and Bert [1] studied the differential quadrature method (DQM) and applied it for the first time to the solution of steady-state oil and gas lubrication problems of self-acting hydrodynamic bearings. In that work, the quadrature solutions of the Reynolds equation for incompressible lubrication were compared with the exact solutions of finite-length bearings. The quadrature solutions of the compressible Reynolds equation for finite-length plain journal bearings were compared with the finite difference and finite element solutions. The work also included comparison of the CPU times of the quadrature solutions with those of the trigonometric series and finite element solutions of oil-lubricated plain slider and journal bearings. For the gas-lubricated journal bearing, the CPU times of the quadrature solutions were compared with those of the finite solutions. In all cases, it was found that the quadrature method was capable of yielding accurate solutions to the lubrication problems and that it was computationally more efficient than other methods of solution. Zhao et al. [2] investigated the subharmonic and quasi-periodic motions of an eccentric squeeze film damper-mounted rigid rotor system. The authors showed that for large values of the rotor unbalance and static misalignment, the sub-harmonic and quasi-periodic motions generated at speeds of more than twice the system critical speed bifurcated from an unstable harmonic solution. Sundararajan and Noah [3] utilized a simple shooting scheme integrated with an arc-length continuation algorithm to analyze the dynamics of periodically forced rotor systems. The results revealed the occurrence of periodic, quasi-periodic, or chaotic motion at different values of the rotor speed. In 2002, Wang et al. [4] investigated static and dynamic characteristics of a flexible rotor supported by externally pressurized porous gas journal bearings. In this work, a modified Reynolds equation was solved by the finite difference method, and a comparative stability of rotor center and journal center was done.In 2006, Rahmatabadi et al. [5, 6] studied the static and dynamic characteristics of noncircular gas journal bearings by considering the effect of mount angles and preload. They proved that noncircular bearings have better dynamic characteristic than circular bearings. Also, by using suitable value of mount angles, stability margin can be increased. Although previous works provide insight into the behavior of the system, the bifurcation and nonlinear dynamic behavior of the gas film in a FRGB has not been examined. Therefore, this paper presents study of nonlinear dynamic behavior of a rigid rotor supported by floating ring gas bearings.Wang [7, 8] analyzed the bifurcation behavior and nonlinear dynamics of flexible and rigid rotors supported by gas journal bearings (relative short gas journal bearings and relative short spherical gas bearing systems) and showed that the rotors exhibited a complex dynamic behavior comprising periodic, sub-harmonic, and quasi-periodic responses at different values of the rotor mass and bearing number, respectively. In 2009, Wang [9] analyzed the bifurcation behavior of a spherical gas journal bearing system by a hybrid method including differential transformation method and finite difference method. The bearing system is modeled as a rigid rotor supported by bearing forces as a result of gas viscosity and rotational speed.The present study analyzes the nonlinear dynamic response of a rigid rotor supported by two floating ring gas bearings. In analyzing the bearing system, the time-dependent motions of the rotor center are described using the Reynolds equation. The modified Reynolds equation is solved by using a hybrid method combining the differential transformation method (DTM) and the finite difference method (FDM). The validity of the hybrid method is confirmed by comparing the results obtained for the orbits of the rotor center with those obtained using the conventional FDM scheme. The proposed method is then applied to analyze the nonlinear dynamic response of the rotor for rotor masses and bearing number in the ranges 0.1 ~ 16.5 kg and 1.0 ~ 7.8, respectively.


## 2. Mathematical Modeling

### 2.1. Governing Equations

A floating ring gas bearing is different from general gas bearing system due to the double thins lubrication. The application for high rotational speed more than 106 rpm can be used by FRGB system, and the stability of this system is focused on to be analyzed. The FRGB system is shown in Figure 1, and Ob(X1,Y1),Oj(X2,Y2),andOr(X3,Y3) are the center of bearing, rotor, and floating ring, respectively.Figure 1
Cross-section of floating ring gas bearing system.In general, the pressure distribution in the gas film between the shaft, floating ring, and the bushing in a floating ring gas bearing (FRGB) system is modeled using the Reynolds equation for an ideal gas. In the present study, the time-dependent motions of the rotor center are modeled with the following form.Inner film:(1)∂∂θ1[Hi3∂Pi2∂θ1]+∂∂ς1[Hi3∂Pi2∂ς1]=2Λi∂(PiHi)∂θ1+2ΛiΩ∂(PiHi)∂θ1+4Λiγ∂(PiHi)∂τ,Outer film:(2)∂∂θ2[Hii3∂Pii2∂θ2]+∂∂ς2[Hii3∂Pii2∂ς2]=2Λii∂(PiiHii)∂θ2+4Λiiγ∂(PiiHii)∂τ,
where Pi,ii is the dimensionless pressure corresponding to the atmospheric pressure Pa; Hi,ii is the dimensionless film thickness between the rotating shaft and the bushing, corresponding to the radial clearance Cr; D is the diameter of bearing, L is the length of bearing, Λi,ii is the bearing number; θ and ς are the coordinates in the circumferential and axial directions, respectively.The FRGB system model incorporates the following design assumptions.(a)
The gas flow in and out of the sides of the bearing (side flow) is neglected.(b)
Assume that the flow is isothermal because the ability of the bearing materials to conduct away heat is greater than the heat generating capacity of the gas film.(c)
As gas viscosity is somewhat insensitive to changes in pressure and the temperature is virtually constant, we assume the gas viscosity to be constant.The gas film pressure distribution must fulfill the following boundary conditions.(a)
Gas pressure on both ends of the housing is equal to the atmospheric pressurePa, P(θ,±ς)=1.(b)
Gas pressureP is an even function for ς,  P(θ,ς)=P(θ,-ς).(c)
Gas pressureP is continuous at ς = 0, ∂P/∂ς|ς=0=0.(d)
Gas pressureP is a periodic function for θ, P(θ,ς)=P(θ+2π,ς), ∂P/∂θ|θ=∂P/∂θ|θ+2π.The analysis presented in this study considers a FRGB system comprising a perfectly balanced rigid rotor of massmr supported symmetrically on two identical floating ring gas bearings, mounted in turn on rigid pedestals. Since the rotor is perfectly balanced and the FRGB is symmetric about its central axes, the current analyses are confined to a single bearing supporting a rotating rotor of mass mr with two degrees of translatory oscillation in the transverse plane.In the transient state, the equation of motion of the rotor can be written as(3)Felx*-Fgfx*=Mrd2X2*dt2,Fely*-Fgfy*=Mrd2Y2*dt2,
in which Felx* and Fely* are the components of the external loading in the horizontal and vertical directions, respectively.And the transformations are introduced as follows:(4)X2=X2*Cr,Y=Y2*Cr,τ=ωt,ddt=ωddτ.
Furthermore, let the following nondimensional groups be defined:
(5)Felx=Felx*PaRL,Fely=Fely*PaRL,Fgfx=Fgfx*PaRL,Fgfy=Fgfy*PaRL,mr=MrCrω2PaRL.Substituting the transformations given in (4) into (3) and introducing the non-dimensional groups defined above, the transient equations of motion of the rotor center can be reformulated as follows:
(6)mrd2X2dτ2=Felx-Fgfx,mrd2Y2dτ2=Fely-Fgfy.


### 2.2. Mathematical Formulation of Numerical Simulations

In solving the modified Reynolds equation using the finite difference method, (1) and (2) are discretized initially using the central-difference scheme in the θ and ς directions and are then discretized once again using the implicit-back-difference scheme in the time domain τ. Note that for simplicity, a uniform mesh size is used. (1) and (2) can be transformed into the following form.Inner film:(7)3((Hi)i,jn+1)2((Hi)i+1,jn+1-(Hi)i-1,jn+12Δθ1)×((Pi2)i+1,jn+1-(Pi2)i-1,jn+12Δθ1)+((Hi)i,jn+1)3((Pi2)i+1,jn+1-2(Pi2)i,jn+1+(Pi2)i-1,jn+1(Δθ1)2)+3((Hi)i,jn+1)2((Hi)i,j+1n+1-(Hi)i,j-1n+12Δς1)×((Pi2)i,j+1n+1-(Pi2)i,j-1n+12Δς1)+((Hi)i,jn+1)3((Pi2)i,j+1n+1-2(Pi2)i,jn+1+(Pi2)i,j-1n+1(Δς1)2)=2Λi((Hi)i+1,jn+1-(Hi)i-1,jn+12Δθ1)(Pi)i,jn+1+2Λi((Pi)i+1,jn+1-(Pi)i-1,jn+12Δθ1)(Hi)i,jn+1+2ΛiΩ((Hi)i,j+1n+1-(Hi)i,j-1n+12Δς1)(Pi)i,jn+1+2ΛiΩ((Pi)i,j+1n+1-(Pi)i,j-1n+12Δς1)(Hi)i,jn+1+4Λiγ(Hi)i,jn+1((Pi)i,jn+1-(Pi)i,jnΔτ)+4Λiγ(Pi)i,jn+1((Hi)i,jn+1-(Hi)i,jnΔτ),Outer film:(8)3((Hii)i,jn+1)2((Hii)i+1,jn+1-(Hii)i-1,jn+12Δθ2)×((Pii2)i+1,jn+1-(Pii2)i-1,jn+12Δθ2)+((Hii)i,jn+1)3((Pii2)i+1,jn+1-2(Pii2)i,jn+1+(Pii2)i-1,jn+1(Δθ2)2)+3((Hii)i,jn+1)2((Hii)i,j+1n+1-(Hii)i,j-1n+12Δς2)×((Pii2)i,j+1n+1-(Pii2)i,j-1n+12Δς2)+((Hii)i,jn+1)3((Pii2)i,j+1n+1-2(Pii2)i,jn+1+(Pii2)i,j-1n+1(Δς2)2)=2Λii((Hii)i+1,jn+1-(Hii)i-1,jn+12Δθ2)(Pii)i,jn+1+2Λii((Pii)i+1,jn+1-(Pii)i-1,jn+12Δθ2)(Hii)i,jn+1+4Λiiγ(Hii)i,jn+1((Pii)i,jn+1-(Pii)i,jnΔτ)+4Λiiγ(Pii)i,jn+1((Hii)i,jn+1-(Hii)i,jnΔτ).The pressure distribution at each time step can then be obtained using an iterative calculation process.The hybrid method proposed in this study is commenced by using the differential transformation method (DTM) [9, 10] to discretize the modified Reynolds equation given in (1) with respect to time. DTM is one of the most widely used techniques for solving nonlinear differential equations due to its rapid convergence rate and minimal calculation error.In solving the Reynolds equation for the current FRGB system, DTM is used to transform the Reynolds equation with respect to the time domainτ, and thus (1) and (2) become as follows.Inner film:(9)3I⊗∂Hi∂θ1⊗∂Q∂θ1+J⊗∂2Q∂θ12+3I⊗∂Hi∂ς1⊗∂Q∂ς1+J⊗∂2Q∂ς12=2ΛiPi⊗∂Hi∂θ1+2ΛiHi⊗∂Pi∂θ1+2ΛiΩPi⊗∂Hi∂θ1+2ΛiΩHi⊗∂Pi∂θ1+4ΛiγPi⊗∂Hi∂τ+4ΛiγHi⊗∂Pi∂τ,Outer film:(10)3I⊗∂Hii∂θ2⊗∂Q∂θ2+J⊗∂2Q∂θ22+3I⊗∂Hii∂ς2⊗∂Q∂ς2+J⊗∂2Q∂ς22=2ΛiiPii⊗∂Hii∂θ2+2ΛiiHii⊗∂Pii∂θ2+4ΛiiγPii⊗∂Hii∂τ+4ΛiiγHii⊗∂Pii∂τ,
where
(11)Q(k)=P2=P⊗P=∑l=0kPi,j(k-l)Pi,j(l),I(k)=H2=H⊗H=∑l=0kHi,j(k-l)Hi,j(l),J(k)=H3=H⊗H⊗H=∑l=0kHi,j(k-l)∑m=0lHi,j(l-m)Hi,j(m).The finite difference method (FDM) is then used to discretize (9) and (10) with respect to the θ and ς directions. Note that (9) and (10) are discretized using the second-order accurate central-difference scheme for both the first and the second derivatives.Substituting (11) into (9)-(10) yields the following form.Inner film:(12)3∑l=0kIi,j(k-l)∑m=0l(Hi)i,j(l-m)·(Qi+1,j(m)-2Qi,j(m)+Qi-1,j(m)(Δθ1)2)+∑l=0kJi,j(k-l)(Qi+1,j(l)-2Qi,j(l)+Qi-1,j(l)(Δθ1)2)+3∑l=0kIi,j(k-l)∑m=0l(Hi)i,j(l-m)·(Qi,j+1(m)-2Qi,j(m)+Qi,j-1(m)(Δς1)2)+∑l=0kJi,j(k-l)(Qi,j+1(l)-2Qi,j(l)+Qi,j-1(l)(Δς1)2)=2Λi∑l=0k(Pi)i,j(k-l)((Hi)i+1,j(l)-(Hi)i-1,j(l)2Δθ1)+2Λi∑l=0k(Hi)i,j(k-l)((Pi)i+1,j(l)-(Pi)i-1,j(l)2Δθ1)+2ΛiΩ∑l=0k(Pi)i,j(k-l)((Hi)i+1,j(l)-(Hi)i-1,j(l)2Δθ1)+2ΛiΩ∑l=0k(Hi)i,j(k-l)((Pi)i+1,j(l)-(Pi)i-1,j(l)2Δθ1)+4Λiγ∑l=0k[(k+1H~)(Pi)i,j(k-l)(Hi)i,j(l+1)]+4Λiγ∑l=0k[(k+1H~)(Hi)i,j(k-l)(Pi)i,j(l+1)],
Outer film:
(13)3∑l=0kIi,j(k-l)∑m=0l(Hii)i,j(l-m)·(Qi+1,j(m)-2Qi,j(m)+Qi-1,j(m)(Δθ2)2)+∑l=0kJi,j(k-l)(Qi+1,j(l)-2Qi,j(l)+Qi-1,j(l)(Δθ2)2)+3∑l=0kIi,j(k-l)∑m=0l(Hii)i,j(l-m)·(Qi,j+1(m)-2Qi,j(m)+Qi,j-1(m)(Δς2)2)+∑l=0kJi,j(k-l)(Qi,j+1(l)-2Qi,j(l)+Qi,j-1(l)(Δς2)2)=2Λii∑l=0k(Pii)i,j(k-l)((Hii)i+1,j(l)-(Hii)i-1,j(l)2Δθ2)+2Λii∑l=0k(Hii)i,j(k-l)((Pii)i+1,j(l)-(Pii)i-1,j(l)2Δθ2)+4Λiiγ∑l=0k[(k+1H~)(Pii)i,j(k-l)(Hii)i,j(l+1)]+4Λiiγ∑l=0k[(k+1H~)(Hii)i,j(k-l)(Pii)i,j(l+1)],
where H~ is time step value.From (9), Pi,j(k) is obtained for each time interval, where i and j are the coordinates of the node position, and k indicates the kth term.The motions of the rotor center are computed using an iterative procedure which commences by determining the acceleration and then computes the velocity and the displacement on a step-by-step basis over time. In defining the initial conditions, the initial displacement (X20,Y20) is specified as the static equilibrium position of the shaft and defines the gap Hi,j(k) between the shaft and the journal bearing, and the velocity of the rotor is assumed to be zero.The iterative computation procedure can be summarized as follows.Step 1.
At timeτ=0, the external loading increases from Felo to Fel. After a time increment Δτ, the new values of the rotor acceleration, velocity, and displacement can be estimated as
(14)d2Xdτ2=Felx-Fgfxmr=AcX,d2Ydτ2=Fely-Fgfymr=AcY,(15)VX=VXo+AcX·Δτ,VY=VYo+AcY·Δτ,(16)X=Xo+VX·Δτ+12AcX·Δτ2,Y=Yo+VY·Δτ+12AcY·Δτ2.Step 2.
The rotor displacement causes a change in the gapH between the rotor and the bushing. Substituting the new value of H into (1) and (2) gives the new pressure distribution in the gap H.Step 3.
The internal force is estimated by integrating the pressure distribution obtained from Step2.Step 4.
The displacement and velocity values from Step1, the pressure distribution from Step 2, and the internal force from Step 3 are taken as the new initial conditions. Using these new conditions, the calculation procedure returns to Step 1 to compute the changes which take place in the time interval Δτ→2Δτ.In this study, the data generated via the iterative computation procedure described above are used to construct power spectra, Poincaré maps, bifurcation diagrams and Lyapunov exponents with which to examine the nonlinear dynamic response of the FRGB system over representative ranges of the rotor mass and bearing number, respectively. For analyzing the behavior of the FRGB system, the time-series data corresponding to the first 1000 revolutions are deliberately excluded in order to ensure that the analyzed data correspond to steady-state conditions.


## 2.1. Governing Equations

A floating ring gas bearing is different from general gas bearing system due to the double thins lubrication. The application for high rotational speed more than 106 rpm can be used by FRGB system, and the stability of this system is focused on to be analyzed. The FRGB system is shown in Figure 1, and Ob(X1,Y1),Oj(X2,Y2),andOr(X3,Y3) are the center of bearing, rotor, and floating ring, respectively.Figure 1
Cross-section of floating ring gas bearing system.In general, the pressure distribution in the gas film between the shaft, floating ring, and the bushing in a floating ring gas bearing (FRGB) system is modeled using the Reynolds equation for an ideal gas. In the present study, the time-dependent motions of the rotor center are modeled with the following form.Inner film:(1)∂∂θ1[Hi3∂Pi2∂θ1]+∂∂ς1[Hi3∂Pi2∂ς1]=2Λi∂(PiHi)∂θ1+2ΛiΩ∂(PiHi)∂θ1+4Λiγ∂(PiHi)∂τ,Outer film:(2)∂∂θ2[Hii3∂Pii2∂θ2]+∂∂ς2[Hii3∂Pii2∂ς2]=2Λii∂(PiiHii)∂θ2+4Λiiγ∂(PiiHii)∂τ,
where Pi,ii is the dimensionless pressure corresponding to the atmospheric pressure Pa; Hi,ii is the dimensionless film thickness between the rotating shaft and the bushing, corresponding to the radial clearance Cr; D is the diameter of bearing, L is the length of bearing, Λi,ii is the bearing number; θ and ς are the coordinates in the circumferential and axial directions, respectively.The FRGB system model incorporates the following design assumptions.(a)
The gas flow in and out of the sides of the bearing (side flow) is neglected.(b)
Assume that the flow is isothermal because the ability of the bearing materials to conduct away heat is greater than the heat generating capacity of the gas film.(c)
As gas viscosity is somewhat insensitive to changes in pressure and the temperature is virtually constant, we assume the gas viscosity to be constant.The gas film pressure distribution must fulfill the following boundary conditions.(a)
Gas pressure on both ends of the housing is equal to the atmospheric pressurePa, P(θ,±ς)=1.(b)
Gas pressureP is an even function for ς,  P(θ,ς)=P(θ,-ς).(c)
Gas pressureP is continuous at ς = 0, ∂P/∂ς|ς=0=0.(d)
Gas pressureP is a periodic function for θ, P(θ,ς)=P(θ+2π,ς), ∂P/∂θ|θ=∂P/∂θ|θ+2π.The analysis presented in this study considers a FRGB system comprising a perfectly balanced rigid rotor of massmr supported symmetrically on two identical floating ring gas bearings, mounted in turn on rigid pedestals. Since the rotor is perfectly balanced and the FRGB is symmetric about its central axes, the current analyses are confined to a single bearing supporting a rotating rotor of mass mr with two degrees of translatory oscillation in the transverse plane.In the transient state, the equation of motion of the rotor can be written as(3)Felx*-Fgfx*=Mrd2X2*dt2,Fely*-Fgfy*=Mrd2Y2*dt2,
in which Felx* and Fely* are the components of the external loading in the horizontal and vertical directions, respectively.And the transformations are introduced as follows:(4)X2=X2*Cr,Y=Y2*Cr,τ=ωt,ddt=ωddτ.
Furthermore, let the following nondimensional groups be defined:
(5)Felx=Felx*PaRL,Fely=Fely*PaRL,Fgfx=Fgfx*PaRL,Fgfy=Fgfy*PaRL,mr=MrCrω2PaRL.Substituting the transformations given in (4) into (3) and introducing the non-dimensional groups defined above, the transient equations of motion of the rotor center can be reformulated as follows:
(6)mrd2X2dτ2=Felx-Fgfx,mrd2Y2dτ2=Fely-Fgfy.


## 2.2. Mathematical Formulation of Numerical Simulations

In solving the modified Reynolds equation using the finite difference method, (1) and (2) are discretized initially using the central-difference scheme in the θ and ς directions and are then discretized once again using the implicit-back-difference scheme in the time domain τ. Note that for simplicity, a uniform mesh size is used. (1) and (2) can be transformed into the following form.Inner film:(7)3((Hi)i,jn+1)2((Hi)i+1,jn+1-(Hi)i-1,jn+12Δθ1)×((Pi2)i+1,jn+1-(Pi2)i-1,jn+12Δθ1)+((Hi)i,jn+1)3((Pi2)i+1,jn+1-2(Pi2)i,jn+1+(Pi2)i-1,jn+1(Δθ1)2)+3((Hi)i,jn+1)2((Hi)i,j+1n+1-(Hi)i,j-1n+12Δς1)×((Pi2)i,j+1n+1-(Pi2)i,j-1n+12Δς1)+((Hi)i,jn+1)3((Pi2)i,j+1n+1-2(Pi2)i,jn+1+(Pi2)i,j-1n+1(Δς1)2)=2Λi((Hi)i+1,jn+1-(Hi)i-1,jn+12Δθ1)(Pi)i,jn+1+2Λi((Pi)i+1,jn+1-(Pi)i-1,jn+12Δθ1)(Hi)i,jn+1+2ΛiΩ((Hi)i,j+1n+1-(Hi)i,j-1n+12Δς1)(Pi)i,jn+1+2ΛiΩ((Pi)i,j+1n+1-(Pi)i,j-1n+12Δς1)(Hi)i,jn+1+4Λiγ(Hi)i,jn+1((Pi)i,jn+1-(Pi)i,jnΔτ)+4Λiγ(Pi)i,jn+1((Hi)i,jn+1-(Hi)i,jnΔτ),Outer film:(8)3((Hii)i,jn+1)2((Hii)i+1,jn+1-(Hii)i-1,jn+12Δθ2)×((Pii2)i+1,jn+1-(Pii2)i-1,jn+12Δθ2)+((Hii)i,jn+1)3((Pii2)i+1,jn+1-2(Pii2)i,jn+1+(Pii2)i-1,jn+1(Δθ2)2)+3((Hii)i,jn+1)2((Hii)i,j+1n+1-(Hii)i,j-1n+12Δς2)×((Pii2)i,j+1n+1-(Pii2)i,j-1n+12Δς2)+((Hii)i,jn+1)3((Pii2)i,j+1n+1-2(Pii2)i,jn+1+(Pii2)i,j-1n+1(Δς2)2)=2Λii((Hii)i+1,jn+1-(Hii)i-1,jn+12Δθ2)(Pii)i,jn+1+2Λii((Pii)i+1,jn+1-(Pii)i-1,jn+12Δθ2)(Hii)i,jn+1+4Λiiγ(Hii)i,jn+1((Pii)i,jn+1-(Pii)i,jnΔτ)+4Λiiγ(Pii)i,jn+1((Hii)i,jn+1-(Hii)i,jnΔτ).The pressure distribution at each time step can then be obtained using an iterative calculation process.The hybrid method proposed in this study is commenced by using the differential transformation method (DTM) [9, 10] to discretize the modified Reynolds equation given in (1) with respect to time. DTM is one of the most widely used techniques for solving nonlinear differential equations due to its rapid convergence rate and minimal calculation error.In solving the Reynolds equation for the current FRGB system, DTM is used to transform the Reynolds equation with respect to the time domainτ, and thus (1) and (2) become as follows.Inner film:(9)3I⊗∂Hi∂θ1⊗∂Q∂θ1+J⊗∂2Q∂θ12+3I⊗∂Hi∂ς1⊗∂Q∂ς1+J⊗∂2Q∂ς12=2ΛiPi⊗∂Hi∂θ1+2ΛiHi⊗∂Pi∂θ1+2ΛiΩPi⊗∂Hi∂θ1+2ΛiΩHi⊗∂Pi∂θ1+4ΛiγPi⊗∂Hi∂τ+4ΛiγHi⊗∂Pi∂τ,Outer film:(10)3I⊗∂Hii∂θ2⊗∂Q∂θ2+J⊗∂2Q∂θ22+3I⊗∂Hii∂ς2⊗∂Q∂ς2+J⊗∂2Q∂ς22=2ΛiiPii⊗∂Hii∂θ2+2ΛiiHii⊗∂Pii∂θ2+4ΛiiγPii⊗∂Hii∂τ+4ΛiiγHii⊗∂Pii∂τ,
where
(11)Q(k)=P2=P⊗P=∑l=0kPi,j(k-l)Pi,j(l),I(k)=H2=H⊗H=∑l=0kHi,j(k-l)Hi,j(l),J(k)=H3=H⊗H⊗H=∑l=0kHi,j(k-l)∑m=0lHi,j(l-m)Hi,j(m).The finite difference method (FDM) is then used to discretize (9) and (10) with respect to the θ and ς directions. Note that (9) and (10) are discretized using the second-order accurate central-difference scheme for both the first and the second derivatives.Substituting (11) into (9)-(10) yields the following form.Inner film:(12)3∑l=0kIi,j(k-l)∑m=0l(Hi)i,j(l-m)·(Qi+1,j(m)-2Qi,j(m)+Qi-1,j(m)(Δθ1)2)+∑l=0kJi,j(k-l)(Qi+1,j(l)-2Qi,j(l)+Qi-1,j(l)(Δθ1)2)+3∑l=0kIi,j(k-l)∑m=0l(Hi)i,j(l-m)·(Qi,j+1(m)-2Qi,j(m)+Qi,j-1(m)(Δς1)2)+∑l=0kJi,j(k-l)(Qi,j+1(l)-2Qi,j(l)+Qi,j-1(l)(Δς1)2)=2Λi∑l=0k(Pi)i,j(k-l)((Hi)i+1,j(l)-(Hi)i-1,j(l)2Δθ1)+2Λi∑l=0k(Hi)i,j(k-l)((Pi)i+1,j(l)-(Pi)i-1,j(l)2Δθ1)+2ΛiΩ∑l=0k(Pi)i,j(k-l)((Hi)i+1,j(l)-(Hi)i-1,j(l)2Δθ1)+2ΛiΩ∑l=0k(Hi)i,j(k-l)((Pi)i+1,j(l)-(Pi)i-1,j(l)2Δθ1)+4Λiγ∑l=0k[(k+1H~)(Pi)i,j(k-l)(Hi)i,j(l+1)]+4Λiγ∑l=0k[(k+1H~)(Hi)i,j(k-l)(Pi)i,j(l+1)],
Outer film:
(13)3∑l=0kIi,j(k-l)∑m=0l(Hii)i,j(l-m)·(Qi+1,j(m)-2Qi,j(m)+Qi-1,j(m)(Δθ2)2)+∑l=0kJi,j(k-l)(Qi+1,j(l)-2Qi,j(l)+Qi-1,j(l)(Δθ2)2)+3∑l=0kIi,j(k-l)∑m=0l(Hii)i,j(l-m)·(Qi,j+1(m)-2Qi,j(m)+Qi,j-1(m)(Δς2)2)+∑l=0kJi,j(k-l)(Qi,j+1(l)-2Qi,j(l)+Qi,j-1(l)(Δς2)2)=2Λii∑l=0k(Pii)i,j(k-l)((Hii)i+1,j(l)-(Hii)i-1,j(l)2Δθ2)+2Λii∑l=0k(Hii)i,j(k-l)((Pii)i+1,j(l)-(Pii)i-1,j(l)2Δθ2)+4Λiiγ∑l=0k[(k+1H~)(Pii)i,j(k-l)(Hii)i,j(l+1)]+4Λiiγ∑l=0k[(k+1H~)(Hii)i,j(k-l)(Pii)i,j(l+1)],
where H~ is time step value.From (9), Pi,j(k) is obtained for each time interval, where i and j are the coordinates of the node position, and k indicates the kth term.The motions of the rotor center are computed using an iterative procedure which commences by determining the acceleration and then computes the velocity and the displacement on a step-by-step basis over time. In defining the initial conditions, the initial displacement (X20,Y20) is specified as the static equilibrium position of the shaft and defines the gap Hi,j(k) between the shaft and the journal bearing, and the velocity of the rotor is assumed to be zero.The iterative computation procedure can be summarized as follows.Step 1.
At timeτ=0, the external loading increases from Felo to Fel. After a time increment Δτ, the new values of the rotor acceleration, velocity, and displacement can be estimated as
(14)d2Xdτ2=Felx-Fgfxmr=AcX,d2Ydτ2=Fely-Fgfymr=AcY,(15)VX=VXo+AcX·Δτ,VY=VYo+AcY·Δτ,(16)X=Xo+VX·Δτ+12AcX·Δτ2,Y=Yo+VY·Δτ+12AcY·Δτ2.Step 2.
The rotor displacement causes a change in the gapH between the rotor and the bushing. Substituting the new value of H into (1) and (2) gives the new pressure distribution in the gap H.Step 3.
The internal force is estimated by integrating the pressure distribution obtained from Step2.Step 4.
The displacement and velocity values from Step1, the pressure distribution from Step 2, and the internal force from Step 3 are taken as the new initial conditions. Using these new conditions, the calculation procedure returns to Step 1 to compute the changes which take place in the time interval Δτ→2Δτ.In this study, the data generated via the iterative computation procedure described above are used to construct power spectra, Poincaré maps, bifurcation diagrams and Lyapunov exponents with which to examine the nonlinear dynamic response of the FRGB system over representative ranges of the rotor mass and bearing number, respectively. For analyzing the behavior of the FRGB system, the time-series data corresponding to the first 1000 revolutions are deliberately excluded in order to ensure that the analyzed data correspond to steady-state conditions.


## 3. Results and Discussions

### 3.1. Numerical Analysis

Table1 compares the results obtained from the FDM and hybrid method (FDM&DTM) for the orbits of the rotor center. It is observed that a good agreement exists between the two sets of results. Moreover, it can be seen that while the FDM suffers numerical instability at low values of the rotor mass and time step, the hybrid method converges under all the considered conditions and therefore represents a more appropriate method for analyzing the nonlinear dynamic response of the FRGB system.Table 1
Comparison of rotor center orbits calculated by FDM and hybrid method, respectively.
Displacement
Conditions
X
2
(
n
T
)
Y
2
(
n
T
)
H
~
=
0.001
H
~
=
0.01
H
~
=
0.001
H
~
=
0.01
SOR and FDM
m
r
=
2.1
kg  ω=3.5×106rpm
0.0171111035
0.0172835103
-
0.0247812154
-
0.0258745091
DTM and FDM
0.0172889076
0.0172858013
-
0.0251754812
-
0.0251713214
SOR and FDM
m
r
=
6.0
kg  ω=3.5×106rpm
0.0182283150
0.0181341276
-
0.1258257815
-
0.1285321157
DTM and FDM
0.0183243783
0.0183297430
-
0.1258701321
-
0.1258781121Tables2 and 3 compare the Poincaré map data calculated by the hybrid method using different time step values, H~, for rotor mass and bearing number values, respectively. For a given rotor mass and bearing number, the rotor center orbits are in agreement to approximately 4 decimal places for the different time steps, H~.Table 2
Comparison of Poincaré maps of rotor center with different values ofmr and H~.
τ
X
2
(
n
T
)
Y
2
(
n
T
)
m
r
=
3.0
kg
π
/
300
-
0.0123810007
-
0.0411895333
π/600
-
0.0123892510
-
0.0411955171
m
r
=
8.0
kg
π
/
300
-
0.0870651990
-
0.2553740891
π/600
-
0.0870179312
-
0.2553712540Table 3
Comparison of Poincaré maps of rotor center with different values of Λ andH~.
τ
X
2
(
n
T
)
Y
2
(
n
T
)
Λ
=
2.2
π
/
300
0.0250423417
-
0.1103104169
π
/
600
0.0250431948
-
0.1104462652
Λ
=
6.1
π
/
300
-
0.2414364192
0.01215789323
π
/
600
-
0.2414519731
0.01213502841


### 3.2. Dynamic Analysis: Situation I

The first dynamic analysis of the FRGB system considered a constant rotational speed with a gradually increasing rotor mass. The gas bearing was loaded with a constant rotational speed of 3.5 × 106 rpm, and the rotor mass mr was varied in the range 0.1≤mr≤16.5 kg.


#### 3.2.1. Dynamic Orbits and Phase Trajectories

Figure2(a) shows that the dynamic orbits of the rotor center are regular at a low value of the rotor mass (mr=3.0 kg) but become irregular at a rotor mass of mr=11.27 kg. For rotor mass values of mr = 13.17, 13.32, 14.13, 14.24, 14.61, 14.74, 15.03, and 15.12 kg, the orbits still exhibit non-periodic characteristic. So, increasing rotor mass causes irregular and non-periodic motion especially when mr is greater than 11.27 kg. Figure 2(b) shows the phase trajectories of the rotor center at different values of the rotor mass. It is observed that the tendencies of the phase trajectories with changes in the rotor mass are identical to those of the rotor center orbits.Trajectories of rotor center atmr=3.0, 11.27, 13.17, 13.32, 14.13, 14.24, 14.61, 14.74, 15.03, 15.12 kg (a); phase trajectories of rotor center (b); power spectra of rotor center in horizontal direction  (c) and vertical direction (d) (at ω=3.5×106 rpm).
(a)
(b)
(c)
(d)


#### 3.2.2. Power Spectra

Figures2(c) and 2(d) show the dynamic responses of the rotor center in the vertical and horizontal directions, respectively. It is seen that the rotor center exhibits T-periodic motion at a rotor mass of mr = 3.0 kg. However, when the rotor mass is increased to mr=11.27 kg, Figures 2(c) and 2(d), second panels, show that the rotor center performs non-periodic motion in the horizontal and vertical directions. For higher rotor mass values of mr = 13.17, 13.32, 14.13, 14.24, 14.61, 14.74, 15.03, and 15.12 kg, the rotor center exhibits sub-harmonic, quasi-periodic, or chaotic motion.


#### 3.2.3. Bifurcation Diagrams, Poincaré Maps, and Maximum Lyapunov Exponents

Figures3(a) and 3(b) present the bifurcation diagrams of the rotor center displacement in the horizontal and vertical directions as a function of the rotor mass mr. Meanwhile, Figures 4(a)–4(j) present the Poincaré maps of the rotor center trajectories at mr = 3.0, 11.27, 13.17, 13.32, 14.13, 14.24, 14.61, 14.74, 15.03, and 15.12, respectively. In Figures 3(a) and 3(b), it can be seen that the rotor center performs T-periodic motion in both the horizontal and the vertical directions at low values of the rotor mass, that is, mr<11.27 kg. This is confirmed by the Poincaré map shown in Figure 4(a) for a rotor mass of 3.0 kg. However, at mr=11.27 kg, the T-periodic motion is replaced by quasi-periodic motion (see Figure 4(b)). From Figure 3, it is seen that this quasi-periodic motion persists over the rotor mass range 11.27≤mr<13.179 kg. However, at  mr=13.17 kg, the quasi-periodic motion is replaced by multiperiodic motion in the horizontal and vertical directions. This motion is maintained over the rotor mass interval 13.17≤mr<13.32 kg. Specifically, the rotor center response transits through a cycle of Chaotic - Multi-T periodic behavior as the rotor mass is increased from 13.32 to 15.52 kg (see Figures 4(c)–4(i)). Then, the behavior of rotor center reverts from Multi-T periodic to chaotic motion for rotor mass values in the range 15.12≤mr≤16.5 kg (see Figure 4(j)). Figure 5, corresponding to rotor mass, show that the maximum Lyapunov exponent has a value of positive, which indicates that the system has a chaotic response over the intervals 14.74≤mr<15.03 and 15.12≤mr≤16.5, shown in Figures 5(a) and 5(b).Bifurcation diagrams: (a)X2(nT) and (b) Y2(nT) versus rotor mass  mr at ω=3.5×106 rpm.
(a)
(b)Poincaré maps of rotor center trajectory at: (a)mr=3.0, (b) 11.27, (c) 13.17, (d) 13.32, (e) 14.13, (f)  14.24, (g) 14.61, (h) 14.74, (i) 15.03, (j) 15.12 kg.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)Maximum Lyapunov exponents of system at different values of rotor mass at: (a)mr=14.74 (b) 15.12 kg.
(a)
(b)From the discussions above, it is evident that the rotor center behavior is significantly dependent on the rotor mass. Table4 summarizes the motions performed by the rotor center for rotor mass values in the interval 0.1≤mr≤16.5 kg.Table 4
Behavior of rotor center at different rotor masses in interval0.1≤mr≤16.5kg.
Rotor massmr  (kg)
Dynamic behavior
0
.
1
≤
m
r
<
11.27
T
11.27
≤
m
r
<
13.17
Quasi
13.17
≤
m
r
<
13.32
Multi-T (8T)
13.32
≤
m
r
<
14.13
Chaos
14.13
≤
m
r
<
14.24
Multi-T (12T)
14.24
≤
m
r
<
14.61
Multi-T
14.61
≤
m
r
<
14.74
Multi-T (14T)
14.74
≤
m
r
<
15.03
Chaos
15.03
≤
m
r
<
15.12
Multi-T (14T)
15.12
≤
m
r
≤
16.5
Chaos


### 3.3. Dynamic Analysis: Situation II

In the second series of analyses, the rotor mass was specified asmr=3.0 kg, and the bearing number Λ was increased over the range 1.0≤Λ≤7.8.


#### 3.3.1. Dynamic Orbits and Phase Trajectories

Figure6(a) shows that the dynamic orbits of the rotor center are regular at a low value of the bearing number (Λ = 3.1) but become irregular at Λ = 5.13. At a bearing number of Λ = 5.29, 5.99, and 6.24, the rotor center resorts to a regular periodic motion. When the bearing number is changed to 5.35 and 7.28, the rotor center performs an irregular motion. Figure 6(b) shows the phase trajectories of the rotor center at different values of the bearing number. It is observed that the results are consistent with those shown in Figure 6(a) for the rotor center orbits, namely a regular motion at Λ = 3.1, 5.29, 5.99, and 6.24, but a non-periodic motion at Λ = 5.13, 5.35, and 7.28.Trajectories of rotor center atΛ=3.1, 5.13, 5.29, 5.35, 5.99, 6.24, 7.28 (a); phase trajectories of rotor center (b), power spectra of rotor center in horizontal direction (c) and vertical direction (d) (at mr=3.0 kg).
(a)
(b)
(c)
(d)


#### 3.3.2. Power Spectra

Figures6(c) and 6(d) show that the rotor center performs a periodic motion at a bearing number of Λ = 3.1 and 5.99. However, when the bearing number is changed to Λ = 5.13, 5.29, and 6.24, the power spectra show that the rotor center performs sub-harmonic motion in both the horizontal and the vertical directions. For values of Λ equal to 5.35, and 7.28, the rotor center performs non-periodic motion.


#### 3.3.3. Bifurcation Diagrams, Poincaré Maps, and Maximum Lyapunov Exponents

Figure7 plots the bifurcation diagrams for the rotor center displacement in the horizontal and vertical directions as a function of the bearing number Λ in the range from 1.0 to 7.8. Figures 8(a)–8(g) present the Poincaré maps of the rotor center trajectories at Λ = 3.1, 5.13, 5.29, 5.35, 5.99, 6.24, and 7.28, respectively. In Figure 7, it can be seen that the rotor center performs T-periodic motion over the bearing number range 1.0≤Λ<5.13 and proved by Figure 8(a). The T-periodic motion becomes unstable at bearing numbers of Λ=5.13 and is replaced by multi-periodic motion (see Figure 8(b)). This multiperiodic behavior persists over 5.13≤Λ<5.29 and also changes to sub-harmonic with a period of multi-T. For bearing numbers in the range 5.29≤Λ≤7.8, the rotor center transits through 7T→Chaotic→T→2T→Chaoticmotion (see Figures 8(c)–8(g)). Figure 9 shows that the maximum Lyapunov exponent has a positive value when Λ equals 5.35 and 7.28 and also indicates that the system has a chaotic response. Meanwhile, over the intervals 5.35≤Λ<5.99, and 7.28≤Λ≤7.8, shown in Figures 7 and 9, the system behaves non-stable behavior and also should be avoided to operate under these parameters.Bifurcation diagrams: (a)X2(nT) and (b) Y2(nT) versus bearing number Λ over interval 1.0≤Λ≤7.8 at mr=3.0 kg.
(a)
(b)Poincaré maps of rotor center trajectory atΛ = (a) 3.1, (b) 5.13, (c) 5.29, (d) 5.35, (e) 5.99, (f) 6.24, (g) 7.28.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Maximum Lyapunov exponents of system at different values of bearing number at: (a)Λ=5.35 (b) 7.28.
(a)
(b)In other words, these bearing number intervals, which correspond to the typical operating conditions of a real-world FRGB system, are characterized byT-, 2T-, 7T-, multi-T and chaotic motions of the rotor center. However, at specific and higher bearing numbers, the rotor performs predominantly chaotic motion in the horizontal and vertical directions. Table 5 summarizes the motions performed by the rotor center for bearing number values in the interval 1.0≤Λ≤7.8.Table 5
Behavior of rotor center at different bearing numbers in interval1.0≤Λ≤7.8.
Bearing numberΛ
Dynamic behavior
1.0
≤
Λ
<
5.13
T
5.13
≤
Λ
<
5.29
Multi-T
5.29
≤
Λ
<
5.35
7T
5.35
≤
Λ
<
5.99
Chaos
5.99
≤
Λ
<
6.24
T
6.24
≤
Λ
<
7.28
2T
7.28
≤
Λ
≤
7.8
Chaos


## 3.1. Numerical Analysis

Table1 compares the results obtained from the FDM and hybrid method (FDM&DTM) for the orbits of the rotor center. It is observed that a good agreement exists between the two sets of results. Moreover, it can be seen that while the FDM suffers numerical instability at low values of the rotor mass and time step, the hybrid method converges under all the considered conditions and therefore represents a more appropriate method for analyzing the nonlinear dynamic response of the FRGB system.Table 1
Comparison of rotor center orbits calculated by FDM and hybrid method, respectively.
Displacement
Conditions
X
2
(
n
T
)
Y
2
(
n
T
)
H
~
=
0.001
H
~
=
0.01
H
~
=
0.001
H
~
=
0.01
SOR and FDM
m
r
=
2.1
kg  ω=3.5×106rpm
0.0171111035
0.0172835103
-
0.0247812154
-
0.0258745091
DTM and FDM
0.0172889076
0.0172858013
-
0.0251754812
-
0.0251713214
SOR and FDM
m
r
=
6.0
kg  ω=3.5×106rpm
0.0182283150
0.0181341276
-
0.1258257815
-
0.1285321157
DTM and FDM
0.0183243783
0.0183297430
-
0.1258701321
-
0.1258781121Tables2 and 3 compare the Poincaré map data calculated by the hybrid method using different time step values, H~, for rotor mass and bearing number values, respectively. For a given rotor mass and bearing number, the rotor center orbits are in agreement to approximately 4 decimal places for the different time steps, H~.Table 2
Comparison of Poincaré maps of rotor center with different values ofmr and H~.
τ
X
2
(
n
T
)
Y
2
(
n
T
)
m
r
=
3.0
kg
π
/
300
-
0.0123810007
-
0.0411895333
π/600
-
0.0123892510
-
0.0411955171
m
r
=
8.0
kg
π
/
300
-
0.0870651990
-
0.2553740891
π/600
-
0.0870179312
-
0.2553712540Table 3
Comparison of Poincaré maps of rotor center with different values of Λ andH~.
τ
X
2
(
n
T
)
Y
2
(
n
T
)
Λ
=
2.2
π
/
300
0.0250423417
-
0.1103104169
π
/
600
0.0250431948
-
0.1104462652
Λ
=
6.1
π
/
300
-
0.2414364192
0.01215789323
π
/
600
-
0.2414519731
0.01213502841


## 3.2. Dynamic Analysis: Situation I

The first dynamic analysis of the FRGB system considered a constant rotational speed with a gradually increasing rotor mass. The gas bearing was loaded with a constant rotational speed of 3.5 × 106 rpm, and the rotor mass mr was varied in the range 0.1≤mr≤16.5 kg.


### 3.2.1. Dynamic Orbits and Phase Trajectories

Figure2(a) shows that the dynamic orbits of the rotor center are regular at a low value of the rotor mass (mr=3.0 kg) but become irregular at a rotor mass of mr=11.27 kg. For rotor mass values of mr = 13.17, 13.32, 14.13, 14.24, 14.61, 14.74, 15.03, and 15.12 kg, the orbits still exhibit non-periodic characteristic. So, increasing rotor mass causes irregular and non-periodic motion especially when mr is greater than 11.27 kg. Figure 2(b) shows the phase trajectories of the rotor center at different values of the rotor mass. It is observed that the tendencies of the phase trajectories with changes in the rotor mass are identical to those of the rotor center orbits.Trajectories of rotor center atmr=3.0, 11.27, 13.17, 13.32, 14.13, 14.24, 14.61, 14.74, 15.03, 15.12 kg (a); phase trajectories of rotor center (b); power spectra of rotor center in horizontal direction  (c) and vertical direction (d) (at ω=3.5×106 rpm).
(a)
(b)
(c)
(d)


### 3.2.2. Power Spectra

Figures2(c) and 2(d) show the dynamic responses of the rotor center in the vertical and horizontal directions, respectively. It is seen that the rotor center exhibits T-periodic motion at a rotor mass of mr = 3.0 kg. However, when the rotor mass is increased to mr=11.27 kg, Figures 2(c) and 2(d), second panels, show that the rotor center performs non-periodic motion in the horizontal and vertical directions. For higher rotor mass values of mr = 13.17, 13.32, 14.13, 14.24, 14.61, 14.74, 15.03, and 15.12 kg, the rotor center exhibits sub-harmonic, quasi-periodic, or chaotic motion.


### 3.2.3. Bifurcation Diagrams, Poincaré Maps, and Maximum Lyapunov Exponents

Figures3(a) and 3(b) present the bifurcation diagrams of the rotor center displacement in the horizontal and vertical directions as a function of the rotor mass mr. Meanwhile, Figures 4(a)–4(j) present the Poincaré maps of the rotor center trajectories at mr = 3.0, 11.27, 13.17, 13.32, 14.13, 14.24, 14.61, 14.74, 15.03, and 15.12, respectively. In Figures 3(a) and 3(b), it can be seen that the rotor center performs T-periodic motion in both the horizontal and the vertical directions at low values of the rotor mass, that is, mr<11.27 kg. This is confirmed by the Poincaré map shown in Figure 4(a) for a rotor mass of 3.0 kg. However, at mr=11.27 kg, the T-periodic motion is replaced by quasi-periodic motion (see Figure 4(b)). From Figure 3, it is seen that this quasi-periodic motion persists over the rotor mass range 11.27≤mr<13.179 kg. However, at  mr=13.17 kg, the quasi-periodic motion is replaced by multiperiodic motion in the horizontal and vertical directions. This motion is maintained over the rotor mass interval 13.17≤mr<13.32 kg. Specifically, the rotor center response transits through a cycle of Chaotic - Multi-T periodic behavior as the rotor mass is increased from 13.32 to 15.52 kg (see Figures 4(c)–4(i)). Then, the behavior of rotor center reverts from Multi-T periodic to chaotic motion for rotor mass values in the range 15.12≤mr≤16.5 kg (see Figure 4(j)). Figure 5, corresponding to rotor mass, show that the maximum Lyapunov exponent has a value of positive, which indicates that the system has a chaotic response over the intervals 14.74≤mr<15.03 and 15.12≤mr≤16.5, shown in Figures 5(a) and 5(b).Bifurcation diagrams: (a)X2(nT) and (b) Y2(nT) versus rotor mass  mr at ω=3.5×106 rpm.
(a)
(b)Poincaré maps of rotor center trajectory at: (a)mr=3.0, (b) 11.27, (c) 13.17, (d) 13.32, (e) 14.13, (f)  14.24, (g) 14.61, (h) 14.74, (i) 15.03, (j) 15.12 kg.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)Maximum Lyapunov exponents of system at different values of rotor mass at: (a)mr=14.74 (b) 15.12 kg.
(a)
(b)From the discussions above, it is evident that the rotor center behavior is significantly dependent on the rotor mass. Table4 summarizes the motions performed by the rotor center for rotor mass values in the interval 0.1≤mr≤16.5 kg.Table 4
Behavior of rotor center at different rotor masses in interval0.1≤mr≤16.5kg.
Rotor massmr  (kg)
Dynamic behavior
0
.
1
≤
m
r
<
11.27
T
11.27
≤
m
r
<
13.17
Quasi
13.17
≤
m
r
<
13.32
Multi-T (8T)
13.32
≤
m
r
<
14.13
Chaos
14.13
≤
m
r
<
14.24
Multi-T (12T)
14.24
≤
m
r
<
14.61
Multi-T
14.61
≤
m
r
<
14.74
Multi-T (14T)
14.74
≤
m
r
<
15.03
Chaos
15.03
≤
m
r
<
15.12
Multi-T (14T)
15.12
≤
m
r
≤
16.5
Chaos


## 3.2.1. Dynamic Orbits and Phase Trajectories

Figure2(a) shows that the dynamic orbits of the rotor center are regular at a low value of the rotor mass (mr=3.0 kg) but become irregular at a rotor mass of mr=11.27 kg. For rotor mass values of mr = 13.17, 13.32, 14.13, 14.24, 14.61, 14.74, 15.03, and 15.12 kg, the orbits still exhibit non-periodic characteristic. So, increasing rotor mass causes irregular and non-periodic motion especially when mr is greater than 11.27 kg. Figure 2(b) shows the phase trajectories of the rotor center at different values of the rotor mass. It is observed that the tendencies of the phase trajectories with changes in the rotor mass are identical to those of the rotor center orbits.Trajectories of rotor center atmr=3.0, 11.27, 13.17, 13.32, 14.13, 14.24, 14.61, 14.74, 15.03, 15.12 kg (a); phase trajectories of rotor center (b); power spectra of rotor center in horizontal direction  (c) and vertical direction (d) (at ω=3.5×106 rpm).
(a)
(b)
(c)
(d)


## 3.2.2. Power Spectra

Figures2(c) and 2(d) show the dynamic responses of the rotor center in the vertical and horizontal directions, respectively. It is seen that the rotor center exhibits T-periodic motion at a rotor mass of mr = 3.0 kg. However, when the rotor mass is increased to mr=11.27 kg, Figures 2(c) and 2(d), second panels, show that the rotor center performs non-periodic motion in the horizontal and vertical directions. For higher rotor mass values of mr = 13.17, 13.32, 14.13, 14.24, 14.61, 14.74, 15.03, and 15.12 kg, the rotor center exhibits sub-harmonic, quasi-periodic, or chaotic motion.


## 3.2.3. Bifurcation Diagrams, Poincaré Maps, and Maximum Lyapunov Exponents

Figures3(a) and 3(b) present the bifurcation diagrams of the rotor center displacement in the horizontal and vertical directions as a function of the rotor mass mr. Meanwhile, Figures 4(a)–4(j) present the Poincaré maps of the rotor center trajectories at mr = 3.0, 11.27, 13.17, 13.32, 14.13, 14.24, 14.61, 14.74, 15.03, and 15.12, respectively. In Figures 3(a) and 3(b), it can be seen that the rotor center performs T-periodic motion in both the horizontal and the vertical directions at low values of the rotor mass, that is, mr<11.27 kg. This is confirmed by the Poincaré map shown in Figure 4(a) for a rotor mass of 3.0 kg. However, at mr=11.27 kg, the T-periodic motion is replaced by quasi-periodic motion (see Figure 4(b)). From Figure 3, it is seen that this quasi-periodic motion persists over the rotor mass range 11.27≤mr<13.179 kg. However, at  mr=13.17 kg, the quasi-periodic motion is replaced by multiperiodic motion in the horizontal and vertical directions. This motion is maintained over the rotor mass interval 13.17≤mr<13.32 kg. Specifically, the rotor center response transits through a cycle of Chaotic - Multi-T periodic behavior as the rotor mass is increased from 13.32 to 15.52 kg (see Figures 4(c)–4(i)). Then, the behavior of rotor center reverts from Multi-T periodic to chaotic motion for rotor mass values in the range 15.12≤mr≤16.5 kg (see Figure 4(j)). Figure 5, corresponding to rotor mass, show that the maximum Lyapunov exponent has a value of positive, which indicates that the system has a chaotic response over the intervals 14.74≤mr<15.03 and 15.12≤mr≤16.5, shown in Figures 5(a) and 5(b).Bifurcation diagrams: (a)X2(nT) and (b) Y2(nT) versus rotor mass  mr at ω=3.5×106 rpm.
(a)
(b)Poincaré maps of rotor center trajectory at: (a)mr=3.0, (b) 11.27, (c) 13.17, (d) 13.32, (e) 14.13, (f)  14.24, (g) 14.61, (h) 14.74, (i) 15.03, (j) 15.12 kg.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)Maximum Lyapunov exponents of system at different values of rotor mass at: (a)mr=14.74 (b) 15.12 kg.
(a)
(b)From the discussions above, it is evident that the rotor center behavior is significantly dependent on the rotor mass. Table4 summarizes the motions performed by the rotor center for rotor mass values in the interval 0.1≤mr≤16.5 kg.Table 4
Behavior of rotor center at different rotor masses in interval0.1≤mr≤16.5kg.
Rotor massmr  (kg)
Dynamic behavior
0
.
1
≤
m
r
<
11.27
T
11.27
≤
m
r
<
13.17
Quasi
13.17
≤
m
r
<
13.32
Multi-T (8T)
13.32
≤
m
r
<
14.13
Chaos
14.13
≤
m
r
<
14.24
Multi-T (12T)
14.24
≤
m
r
<
14.61
Multi-T
14.61
≤
m
r
<
14.74
Multi-T (14T)
14.74
≤
m
r
<
15.03
Chaos
15.03
≤
m
r
<
15.12
Multi-T (14T)
15.12
≤
m
r
≤
16.5
Chaos


## 3.3. Dynamic Analysis: Situation II

In the second series of analyses, the rotor mass was specified asmr=3.0 kg, and the bearing number Λ was increased over the range 1.0≤Λ≤7.8.


### 3.3.1. Dynamic Orbits and Phase Trajectories

Figure6(a) shows that the dynamic orbits of the rotor center are regular at a low value of the bearing number (Λ = 3.1) but become irregular at Λ = 5.13. At a bearing number of Λ = 5.29, 5.99, and 6.24, the rotor center resorts to a regular periodic motion. When the bearing number is changed to 5.35 and 7.28, the rotor center performs an irregular motion. Figure 6(b) shows the phase trajectories of the rotor center at different values of the bearing number. It is observed that the results are consistent with those shown in Figure 6(a) for the rotor center orbits, namely a regular motion at Λ = 3.1, 5.29, 5.99, and 6.24, but a non-periodic motion at Λ = 5.13, 5.35, and 7.28.Trajectories of rotor center atΛ=3.1, 5.13, 5.29, 5.35, 5.99, 6.24, 7.28 (a); phase trajectories of rotor center (b), power spectra of rotor center in horizontal direction (c) and vertical direction (d) (at mr=3.0 kg).
(a)
(b)
(c)
(d)


### 3.3.2. Power Spectra

Figures6(c) and 6(d) show that the rotor center performs a periodic motion at a bearing number of Λ = 3.1 and 5.99. However, when the bearing number is changed to Λ = 5.13, 5.29, and 6.24, the power spectra show that the rotor center performs sub-harmonic motion in both the horizontal and the vertical directions. For values of Λ equal to 5.35, and 7.28, the rotor center performs non-periodic motion.


### 3.3.3. Bifurcation Diagrams, Poincaré Maps, and Maximum Lyapunov Exponents

Figure7 plots the bifurcation diagrams for the rotor center displacement in the horizontal and vertical directions as a function of the bearing number Λ in the range from 1.0 to 7.8. Figures 8(a)–8(g) present the Poincaré maps of the rotor center trajectories at Λ = 3.1, 5.13, 5.29, 5.35, 5.99, 6.24, and 7.28, respectively. In Figure 7, it can be seen that the rotor center performs T-periodic motion over the bearing number range 1.0≤Λ<5.13 and proved by Figure 8(a). The T-periodic motion becomes unstable at bearing numbers of Λ=5.13 and is replaced by multi-periodic motion (see Figure 8(b)). This multiperiodic behavior persists over 5.13≤Λ<5.29 and also changes to sub-harmonic with a period of multi-T. For bearing numbers in the range 5.29≤Λ≤7.8, the rotor center transits through 7T→Chaotic→T→2T→Chaoticmotion (see Figures 8(c)–8(g)). Figure 9 shows that the maximum Lyapunov exponent has a positive value when Λ equals 5.35 and 7.28 and also indicates that the system has a chaotic response. Meanwhile, over the intervals 5.35≤Λ<5.99, and 7.28≤Λ≤7.8, shown in Figures 7 and 9, the system behaves non-stable behavior and also should be avoided to operate under these parameters.Bifurcation diagrams: (a)X2(nT) and (b) Y2(nT) versus bearing number Λ over interval 1.0≤Λ≤7.8 at mr=3.0 kg.
(a)
(b)Poincaré maps of rotor center trajectory atΛ = (a) 3.1, (b) 5.13, (c) 5.29, (d) 5.35, (e) 5.99, (f) 6.24, (g) 7.28.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Maximum Lyapunov exponents of system at different values of bearing number at: (a)Λ=5.35 (b) 7.28.
(a)
(b)In other words, these bearing number intervals, which correspond to the typical operating conditions of a real-world FRGB system, are characterized byT-, 2T-, 7T-, multi-T and chaotic motions of the rotor center. However, at specific and higher bearing numbers, the rotor performs predominantly chaotic motion in the horizontal and vertical directions. Table 5 summarizes the motions performed by the rotor center for bearing number values in the interval 1.0≤Λ≤7.8.Table 5
Behavior of rotor center at different bearing numbers in interval1.0≤Λ≤7.8.
Bearing numberΛ
Dynamic behavior
1.0
≤
Λ
<
5.13
T
5.13
≤
Λ
<
5.29
Multi-T
5.29
≤
Λ
<
5.35
7T
5.35
≤
Λ
<
5.99
Chaos
5.99
≤
Λ
<
6.24
T
6.24
≤
Λ
<
7.28
2T
7.28
≤
Λ
≤
7.8
Chaos


## 3.3.1. Dynamic Orbits and Phase Trajectories

Figure6(a) shows that the dynamic orbits of the rotor center are regular at a low value of the bearing number (Λ = 3.1) but become irregular at Λ = 5.13. At a bearing number of Λ = 5.29, 5.99, and 6.24, the rotor center resorts to a regular periodic motion. When the bearing number is changed to 5.35 and 7.28, the rotor center performs an irregular motion. Figure 6(b) shows the phase trajectories of the rotor center at different values of the bearing number. It is observed that the results are consistent with those shown in Figure 6(a) for the rotor center orbits, namely a regular motion at Λ = 3.1, 5.29, 5.99, and 6.24, but a non-periodic motion at Λ = 5.13, 5.35, and 7.28.Trajectories of rotor center atΛ=3.1, 5.13, 5.29, 5.35, 5.99, 6.24, 7.28 (a); phase trajectories of rotor center (b), power spectra of rotor center in horizontal direction (c) and vertical direction (d) (at mr=3.0 kg).
(a)
(b)
(c)
(d)


## 3.3.2. Power Spectra

Figures6(c) and 6(d) show that the rotor center performs a periodic motion at a bearing number of Λ = 3.1 and 5.99. However, when the bearing number is changed to Λ = 5.13, 5.29, and 6.24, the power spectra show that the rotor center performs sub-harmonic motion in both the horizontal and the vertical directions. For values of Λ equal to 5.35, and 7.28, the rotor center performs non-periodic motion.


## 3.3.3. Bifurcation Diagrams, Poincaré Maps, and Maximum Lyapunov Exponents

Figure7 plots the bifurcation diagrams for the rotor center displacement in the horizontal and vertical directions as a function of the bearing number Λ in the range from 1.0 to 7.8. Figures 8(a)–8(g) present the Poincaré maps of the rotor center trajectories at Λ = 3.1, 5.13, 5.29, 5.35, 5.99, 6.24, and 7.28, respectively. In Figure 7, it can be seen that the rotor center performs T-periodic motion over the bearing number range 1.0≤Λ<5.13 and proved by Figure 8(a). The T-periodic motion becomes unstable at bearing numbers of Λ=5.13 and is replaced by multi-periodic motion (see Figure 8(b)). This multiperiodic behavior persists over 5.13≤Λ<5.29 and also changes to sub-harmonic with a period of multi-T. For bearing numbers in the range 5.29≤Λ≤7.8, the rotor center transits through 7T→Chaotic→T→2T→Chaoticmotion (see Figures 8(c)–8(g)). Figure 9 shows that the maximum Lyapunov exponent has a positive value when Λ equals 5.35 and 7.28 and also indicates that the system has a chaotic response. Meanwhile, over the intervals 5.35≤Λ<5.99, and 7.28≤Λ≤7.8, shown in Figures 7 and 9, the system behaves non-stable behavior and also should be avoided to operate under these parameters.Bifurcation diagrams: (a)X2(nT) and (b) Y2(nT) versus bearing number Λ over interval 1.0≤Λ≤7.8 at mr=3.0 kg.
(a)
(b)Poincaré maps of rotor center trajectory atΛ = (a) 3.1, (b) 5.13, (c) 5.29, (d) 5.35, (e) 5.99, (f) 6.24, (g) 7.28.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Maximum Lyapunov exponents of system at different values of bearing number at: (a)Λ=5.35 (b) 7.28.
(a)
(b)In other words, these bearing number intervals, which correspond to the typical operating conditions of a real-world FRGB system, are characterized byT-, 2T-, 7T-, multi-T and chaotic motions of the rotor center. However, at specific and higher bearing numbers, the rotor performs predominantly chaotic motion in the horizontal and vertical directions. Table 5 summarizes the motions performed by the rotor center for bearing number values in the interval 1.0≤Λ≤7.8.Table 5
Behavior of rotor center at different bearing numbers in interval1.0≤Λ≤7.8.
Bearing numberΛ
Dynamic behavior
1.0
≤
Λ
<
5.13
T
5.13
≤
Λ
<
5.29
Multi-T
5.29
≤
Λ
<
5.35
7T
5.35
≤
Λ
<
5.99
Chaos
5.99
≤
Λ
<
6.24
T
6.24
≤
Λ
<
7.28
2T
7.28
≤
Λ
≤
7.8
Chaos


## 4. Conclusions

This study has analyzed the nonlinear dynamic behavior of a rigid rotor supported by floating ring gas bearing (FRGB) system and utilized a hybrid numerical scheme comprising the differential transformation method (DTM) and the finite difference method (FDM) to obtain the pressure distribution of FRGB system. The system dynamic orbits, phase trajectories, power spectra, bifurcation diagrams, Poincaré maps, and maximum Lyapunov exponents have revealed the presence of a complex dynamic behavior comprising periodic, sub-harmonic, quasi-periodic and chaotic responses of the rotor center. The proposed hybrid method avoids the numerical instability problem suffered by the finite difference scheme at low values of the rotor mass and computational time-step. Moreover, the DTM&FDM method converges under all the considered conditions and therefore represents a more appropriate method for analyzing the nonlinear dynamic response of the FRGB system. The results of this study provide an understanding of the nonlinear dynamic behavior of floating ring gas bearing systems characterized by different rotor massesmr and bearing numbers Λ. Specifically, the results have shown that at a rotor mass of mr=11.27 kg, the Poincaré map has the form of a closed curve, indicating that the rotor center performs quasi-periodic motion. And chaotic motion appears in three rotor mass intervals to be avoided. Regarding the influence of the bearing number on the dynamic response of the bearing system, the results have shown that at Λ=3.1, the Poincaré map contains one discrete points, indicating the presence of T-periodic motion. However, when the bearing number is increased to Λ=5.35, the system is transferred to chaotic motion. Then, chaotic motion appears in two intervals in the bearing number range 5.35≤Λ≤7.8. As shown in Table 5, the intermediate regions in the bearing number range are characterized by T- and 2T-motion.



---

*Source: 145716-2013-05-27.xml*