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

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 systemover representative ranges of the rotormass 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.


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 10 6 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 gaslubricated 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 dampermounted 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 arclength 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 timedependent 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.

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 10 6 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   ( 1 ,  1 ),   ( 2 ,  2 ), and   ( 3 ,  3 ) are the center of bearing, rotor, and floating ring, respectively.
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 timedependent motions of the rotor center are modeled with the following form.
Inner film: Outer film: where  , is the dimensionless pressure corresponding to the atmospheric pressure   ;  , is the dimensionless film thickness between the rotating shaft and the bushing, corresponding to the radial clearance   ;  is the diameter of bearing,  is the length of bearing, Λ , 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 pressure   , (, ±) = 1.
The analysis presented in this study considers a FRGB system comprising a perfectly balanced rigid rotor of mass   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   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 in which  * elx and  * ely are the components of the external loading in the horizontal and vertical directions, respectively.
And the transformations are introduced as follows: Furthermore, let the following nondimensional groups be defined: 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:

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-backdifference 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: × ( + ((  ) Outer film: × (

Δ𝜏
) . ( 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: Outer film: where (11) 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: Outer film: where H is time step value.From (9),  , () is obtained for each time interval, where  and  are the coordinates of the node position, and  indicates the th 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 ( 20 , 20 ) is specified as the static equilibrium position of the shaft and defines the gap  , () 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  elo to  el .After a time increment Δ, the new values of the rotor Step 2. The rotor displacement causes a change in the gap  between the rotor and the bushing.Substituting the new value of  into (1) and ( 2) gives the new pressure distribution in the gap .
Step 3. The internal force is estimated by integrating the pressure distribution obtained from Step 2.
Step 4. The displacement and velocity values from Step 1, 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 andLyapunov 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.

Numerical Analysis. Table 1 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.
Tables 2 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.

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 × 10 6 rpm, and the rotor mass   was varied in the range 0.1 ≤   ≤ 16.5 kg.

Dynamic Orbits and Phase Trajectories. Figure 2(a)
shows that the dynamic orbits of the rotor center are regular at a low value of the rotor mass (  = 3.0 kg) but become irregular at a rotor mass of   = 11.27kg.For rotor mass values of   = 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   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.2(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   = 3.0 kg.However, when the rotor mass is increased to   = 11.27kg,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   = 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.

Bifurcation Diagrams, Poincaré Maps, and Maximum
Lyapunov Exponents.Figures 3(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   .Meanwhile, Figures 4(a)-4(j) present the Poincaré maps of the rotor center trajectories at   = 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,   < 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   = 11.27kg, 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 ≤   < 13.179 kg.However, at   = 13.17kg, 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 ≤   < 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 ≤   ≤ 16.5 kg (see Figure 4(j)).Figure 5, corresponding to rotor mass, show that the maximum Lyapunov exponent hasa value of positive, which indicates that the system has a chaotic response over the intervals 14.74 ≤   < 15.03 and 15.12 ≤   ≤ 16.5, shown in Figures 5(a From the discussions above, it is evident that the rotor center behavior is significantly dependent on the rotor mass.Table 4 summarizes the motions performed by the rotor center for rotor mass values in the interval 0.1 ≤   ≤ 16.5 kg.

Dynamic Analysis: Situation II.
In the second series of analyses, the rotor mass was specified as   = 3.0 kg, and the bearing number Λ was increased over the range 1.0 ≤ Λ ≤ 7.8.

Dynamic Orbits and Phase Trajectories. Figure 6(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.

Power Spectra. Figures 6(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.

Bifurcation Diagrams, Poincaré Maps, and Maximum
Lyapunov Exponents..17,(d) 13.32, (e) 14.13, (f) 14.24, (g) 14.61, (h) 14.74, (i) 15.03, (j) 15.12 kg.  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 -periodic motion over the bearing number range 1.0 ≤ Λ < 5.13 and proved by Figure 8(a).The -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-.For bearing numbers in the range 5.29 ≤ Λ ≤ 7.8, the rotor center transits through 7 → Chaotic →  → 2 → Chaotic motion (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.In other words, these bearing number intervals, which correspond to the typical operating conditions of a real-world FRGB system, are characterized by -, 2-, 7-, multi- 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.

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, subharmonic, 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 masses   and bearing numbers Λ.Specifically, the results have shown that at a rotor mass of   = 11.27kg, 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 -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 and 2motion.

Figure 1 :
Figure 1: Cross-section of floating ring gas bearing system.

Table 1 :
Comparison of rotor center orbits calculated by FDM and hybrid method, respectively.

Table 2 :
Comparison of Poincaré maps of rotor center with different values of   and H.

Table 3 :
Comparison of Poincaré maps of rotor center with different values of Λ and H.

Table 4 :
Behavior of rotor center at different rotor masses in interval 0.1 ≤   ≤ 16.5 kg.

Table 5 :
Behavior of rotor center at different bearing numbers in interval 1.0 ≤ Λ ≤ 7.8.