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

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.

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 10⁶ 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 10⁶ 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 𝑂_𝑏(𝑋₁,𝑌₁),𝑂_𝑗(𝑋₂,𝑌₂),and𝑂_𝑟(𝑋₃,𝑌₃) 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: 𝜕𝜕𝜃₁𝐻³_𝑖𝜕𝑃²_𝑖𝜕𝜃₁+𝜕𝜕𝜍₁𝐻³_𝑖𝜕𝑃²_𝑖𝜕𝜍₁=2Λ_𝑖𝜕𝑃_𝑖𝐻_𝑖𝜕𝜃₁+2Λ_𝑖Ω𝜕𝑃_𝑖𝐻_𝑖𝜕𝜃₁+4Λ_𝑖𝛾𝜕𝑃_𝑖𝐻_𝑖,𝜕𝜏 (1)

Outer film: 𝜕𝜕𝜃₂𝐻³_𝑖𝑖𝜕𝑃²_𝑖𝑖𝜕𝜃₂+𝜕𝜕𝜍₂𝐻³_𝑖𝑖𝜕𝑃²_𝑖𝑖𝜕𝜍₂=2Λ_𝑖𝑖𝜕𝑃_𝑖𝑖𝐻_𝑖𝑖𝜕𝜃₂+4Λ_𝑖𝑖𝛾𝜕𝑃_𝑖𝑖𝐻_𝑖𝑖,𝜕𝜏 (2) 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.(b)Gas pressure 𝑃 is an even function for 𝜍, 𝑃(𝜃,𝜍)=𝑃(𝜃,−𝜍).(c)Gas pressure 𝑃 is continuous at 𝜍 = 0, 𝜕𝑃/𝜕𝜍|_𝜍=0=0.(d)Gas pressure 𝑃 is a periodic function for 𝜃, 𝑃(𝜃,𝜍)=𝑃(𝜃+2𝜋,𝜍), 𝜕𝑃/𝜕𝜃|_𝜃=𝜕𝑃/𝜕𝜃|_𝜃+2𝜋.

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 𝐹^∗_elx−𝐹^∗_gfx=𝑀_𝑟𝑑²𝑋^∗₂𝑑𝑡²,𝐹^∗_ely−𝐹^∗_gfy=𝑀_𝑟𝑑²𝑌^∗₂𝑑𝑡², (3) 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: 𝑋₂=𝑋^∗₂𝐶_𝑟𝑌,𝑌=^∗2𝐶_𝑟𝑑,𝜏=𝜔𝑡,𝑑𝑑𝑡=𝜔.𝑑𝜏 (4) Furthermore, let the following nondimensional groups be defined: 𝐹_elx=𝐹^∗_elx𝑃_𝑎𝑅𝐿,𝐹_ely=𝐹^∗_ely𝑃_𝑎,𝐹𝑅𝐿_gfx=𝐹^∗_gfx𝑃_𝑎𝑅𝐿,𝐹_gfy=𝐹^∗_gfy𝑃_𝑎𝑅𝐿,𝑚_𝑟=𝑀_𝑟𝐶_𝑟𝜔²𝑃_𝑎.𝑅𝐿 (5)

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: 𝑚_𝑟𝑑²𝑋₂𝑑𝜏²=𝐹_elx−𝐹_gfx,𝑚_𝑟𝑑²𝑌₂𝑑𝜏²=𝐹_ely−𝐹_gfy. (6)

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: 3𝐻_𝑖^𝑛+1_𝑖,𝑗²𝐻_𝑖^𝑛+1_𝑖+1,𝑗−𝐻_𝑖^𝑛+1_{𝑖−1,𝑗}2Δ𝜃₁×𝑃²_𝑖^𝑛+1_𝑖+1,𝑗−𝑃²_𝑖^𝑛+1_{𝑖−1,𝑗}2Δ𝜃₁+𝐻_𝑖^𝑛+1_𝑖,𝑗³𝑃²_𝑖^𝑛+1_𝑖+1,𝑗−2𝑃²_𝑖^𝑛+1_𝑖,𝑗+𝑃²_𝑖^𝑛+1_{𝑖−1,𝑗}Δ𝜃₁²+3𝐻_𝑖^𝑛+1_𝑖,𝑗²𝐻_𝑖^𝑛+1_𝑖,𝑗+1−𝐻_𝑖^𝑛+1_{𝑖,𝑗−1}2Δ𝜍₁×𝑃²_𝑖^𝑛+1_𝑖,𝑗+1−𝑃²_𝑖^𝑛+1_{𝑖,𝑗−1}2Δ𝜍₁+𝐻_𝑖^𝑛+1_𝑖,𝑗³𝑃²_𝑖^𝑛+1_𝑖,𝑗+1−2𝑃²_𝑖^𝑛+1_𝑖,𝑗+𝑃²_𝑖^𝑛+1_{𝑖,𝑗−1}Δ𝜍₁²=2Λ_𝑖𝐻_𝑖^𝑛+1_𝑖+1,𝑗−𝐻_𝑖^𝑛+1_{𝑖−1,𝑗}2Δ𝜃₁𝑃_𝑖^𝑛+1_𝑖,𝑗+2Λ_𝑖𝑃_𝑖^𝑛+1_𝑖+1,𝑗−𝑃_𝑖^𝑛+1_{𝑖−1,𝑗}2Δ𝜃₁𝐻_𝑖^𝑛+1_𝑖,𝑗+2Λ_𝑖Ω𝐻_𝑖^𝑛+1_𝑖,𝑗+1−𝐻_𝑖^𝑛+1_{𝑖,𝑗−1}2Δ𝜍₁𝑃_𝑖^𝑛+1_𝑖,𝑗+2Λ_𝑖Ω𝑃_𝑖^𝑛+1_𝑖,𝑗+1−𝑃_𝑖^𝑛+1_{𝑖,𝑗−1}2Δ𝜍₁𝐻_𝑖^𝑛+1_𝑖,𝑗+4Λ_𝑖𝛾𝐻_𝑖^𝑛+1_𝑖,𝑗𝑃_𝑖^𝑛+1_𝑖,𝑗−𝑃_𝑖^𝑛_𝑖,𝑗Δ𝜏+4Λ_𝑖𝛾𝑃_𝑖^𝑛+1_𝑖,𝑗𝐻_𝑖^𝑛+1_𝑖,𝑗−𝐻_𝑖^𝑛_𝑖,𝑗,Δ𝜏 (7)

Outer film: 3𝐻_𝑖𝑖^𝑛+1_𝑖,𝑗²𝐻_𝑖𝑖^𝑛+1_𝑖+1,𝑗−𝐻_𝑖𝑖^𝑛+1_{𝑖−1,𝑗}2Δ𝜃₂×𝑃²_𝑖𝑖^𝑛+1_𝑖+1,𝑗−𝑃²_𝑖𝑖^𝑛+1_{𝑖−1,𝑗}2Δ𝜃₂+𝐻_𝑖𝑖^𝑛+1_𝑖,𝑗³𝑃²_𝑖𝑖^𝑛+1_𝑖+1,𝑗−2𝑃²_𝑖𝑖^𝑛+1_𝑖,𝑗+𝑃²_𝑖𝑖^𝑛+1_{𝑖−1,𝑗}Δ𝜃₂²+3𝐻_𝑖𝑖^𝑛+1_𝑖,𝑗²𝐻_𝑖𝑖^𝑛+1_𝑖,𝑗+1−𝐻_𝑖𝑖^𝑛+1_{𝑖,𝑗−1}2Δ𝜍₂×𝑃²_𝑖𝑖^𝑛+1_𝑖,𝑗+1−𝑃²_𝑖𝑖^𝑛+1_{𝑖,𝑗−1}2Δ𝜍₂+𝐻_𝑖𝑖^𝑛+1_𝑖,𝑗³𝑃²_𝑖𝑖^𝑛+1_𝑖,𝑗+1−2𝑃²_𝑖𝑖^𝑛+1_𝑖,𝑗+𝑃²_𝑖𝑖^𝑛+1_{𝑖,𝑗−1}Δ𝜍₂²=2Λ_𝑖𝑖𝐻_𝑖𝑖^𝑛+1_𝑖+1,𝑗−𝐻_𝑖𝑖^𝑛+1_{𝑖−1,𝑗}2Δ𝜃₂(𝑃_𝑖𝑖)^𝑛+1_𝑖,𝑗+2Λ_𝑖𝑖𝑃_𝑖𝑖^𝑛+1_𝑖+1,𝑗−𝑃_𝑖𝑖^𝑛+1_{𝑖−1,𝑗}2Δ𝜃₂𝐻_𝑖𝑖^𝑛+1_𝑖,𝑗+4Λ_𝑖𝑖𝛾𝐻_𝑖𝑖^𝑛+1_𝑖,𝑗𝑃_𝑖𝑖^𝑛+1_𝑖,𝑗−𝑃_𝑖𝑖^𝑛_𝑖,𝑗Δ𝜏+4Λ_𝑖𝑖𝛾𝑃_𝑖𝑖^𝑛+1_𝑖,𝑗𝐻_𝑖𝑖^𝑛+1_𝑖,𝑗−𝐻_𝑖𝑖^𝑛_𝑖,𝑗.Δ𝜏 (8)

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: 3𝐼⊗𝜕𝐻_𝑖𝜕𝜃₁⊗𝜕𝑄𝜕𝜃₁𝜕+𝐽⊗²𝑄𝜕𝜃²₁+3𝐼⊗𝜕𝐻_𝑖𝜕𝜍₁⊗𝜕𝑄𝜕𝜍₁𝜕+𝐽⊗²𝑄𝜕𝜍²₁=2Λ_𝑖𝑃_𝑖⊗𝜕𝐻_𝑖𝜕𝜃₁+2Λ_𝑖𝐻_𝑖⊗𝜕𝑃_𝑖𝜕𝜃₁+2Λ_𝑖Ω𝑃_𝑖⊗𝜕𝐻_𝑖𝜕𝜃₁+2Λ_𝑖Ω𝐻_𝑖⊗𝜕𝑃_𝑖𝜕𝜃₁+4Λ_𝑖𝛾𝑃_𝑖⊗𝜕𝐻_𝑖𝜕𝜏+4Λ_𝑖𝛾𝐻_𝑖⊗𝜕𝑃_𝑖,𝜕𝜏 (9)

Outer film: 3𝐼⊗𝜕𝐻_𝑖𝑖𝜕𝜃₂⊗𝜕𝑄𝜕𝜃₂𝜕+𝐽⊗²𝑄𝜕𝜃²₂+3𝐼⊗𝜕𝐻_𝑖𝑖𝜕𝜍₂⊗𝜕𝑄𝜕𝜍₂𝜕+𝐽⊗²𝑄𝜕𝜍²₂=2Λ_𝑖𝑖𝑃_𝑖𝑖⊗𝜕𝐻_𝑖𝑖𝜕𝜃₂+2Λ_𝑖𝑖𝐻_𝑖𝑖⊗𝜕𝑃_𝑖𝑖𝜕𝜃₂+4Λ_𝑖𝑖𝛾𝑃_𝑖𝑖⊗𝜕𝐻_𝑖𝑖𝜕𝜏+4Λ_𝑖𝑖𝛾𝐻_𝑖𝑖⊗𝜕𝑃_𝑖𝑖,𝜕𝜏 (10) where 𝑄(𝑘)=𝑃²=𝑃⊗𝑃=^𝑘_𝑙=0𝑃_𝑖,𝑗(𝑘−𝑙)𝑃_𝑖,𝑗(𝑙),𝐼(𝑘)=𝐻²=𝐻⊗𝐻=^𝑘_𝑙=0𝐻_𝑖,𝑗(𝑘−𝑙)𝐻_𝑖,𝑗(𝑙),𝐽(𝑘)=𝐻³==𝐻⊗𝐻⊗𝐻_𝑘𝑙=0𝐻_𝑖,𝑗(𝑘−𝑙)^𝑙_𝑚=0𝐻_𝑖,𝑗(𝑙−𝑚)𝐻_𝑖,𝑗(𝑚). (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: 3^𝑘_𝑙=0𝐼_𝑖,𝑗(𝑘−𝑙)^𝑙_𝑚=0𝐻_{𝑖𝑖,𝑗}()⋅𝑄𝑙−𝑚_𝑖+1,𝑗(𝑚)−2𝑄_𝑖,𝑗(𝑚)+𝑄_{𝑖−1,𝑗}(𝑚)Δ𝜃₁²+^𝑘_𝑙=0𝐽_𝑖,𝑗()𝑄𝑘−𝑙^𝑖+1,𝑗(𝑙)−2𝑄_𝑖,𝑗(𝑙)+𝑄_{𝑖−1,𝑗}(𝑙)Δ𝜃₁²+3^𝑘_𝑙=0𝐼_𝑖,𝑗(𝑘−𝑙)^𝑙_𝑚=0𝐻_{𝑖𝑖,𝑗}()⋅𝑄𝑙−𝑚_𝑖,𝑗+1(𝑚)−2𝑄_𝑖,𝑗(𝑚)+𝑄_{𝑖,𝑗−1}(𝑚)Δ𝜍₁²+^𝑘_𝑙=0𝐽_𝑖,𝑗()𝑄𝑘−𝑙^𝑖,𝑗+1(𝑙)−2𝑄_𝑖,𝑗(𝑙)+𝑄_{𝑖,𝑗−1}(𝑙)Δ𝜍₁²=2Λ_𝑖^𝑘_𝑙=0𝑃_{𝑖𝑖,𝑗}()𝐻𝑘−𝑙^{𝑖𝑖+1,𝑗}(𝑙)−𝐻_{𝑖𝑖−1,𝑗}(𝑙)2Δ𝜃₁+2Λ_𝑖^𝑘_𝑙=0𝐻_{𝑖𝑖,𝑗}()𝑃𝑘−𝑙^{𝑖𝑖+1,𝑗}(𝑙)−𝑃_{𝑖𝑖−1,𝑗}(𝑙)2Δ𝜃₁+2Λ_𝑖Ω^𝑘_𝑙=0𝑃_{𝑖𝑖,𝑗}()𝐻𝑘−𝑙^{𝑖𝑖+1,𝑗}(𝑙)−𝐻_{𝑖𝑖−1,𝑗}(𝑙)2Δ𝜃₁+2Λ_𝑖Ω^𝑘_𝑙=0𝐻_{𝑖𝑖,𝑗}()𝑃𝑘−𝑙^{𝑖𝑖+1,𝑗}(𝑙)−𝑃_{𝑖𝑖−1,𝑗}(𝑙)2Δ𝜃₁+4Λ_𝑖𝛾^𝑘_𝑙=0𝑘+1𝐻𝑃_{𝑖𝑖,𝑗}(𝑘−𝑙)𝐻_{𝑖𝑖,𝑗}(𝑙+1)+4Λ_𝑖𝛾^𝑘_𝑙=0𝑘+1𝐻𝐻_{𝑖𝑖,𝑗}(𝑘−𝑙)𝑃_{𝑖𝑖,𝑗}(𝑙+1), (12) Outer film: 3^𝑘_𝑙=0𝐼_𝑖,𝑗(𝑘−𝑙)^𝑙_𝑚=0𝐻_{𝑖𝑖𝑖,𝑗}()⋅𝑄𝑙−𝑚_𝑖+1,𝑗(𝑚)−2𝑄_𝑖,𝑗(𝑚)+𝑄_{𝑖−1,𝑗}(𝑚)Δ𝜃₂²+^𝑘_𝑙=0𝐽_𝑖,𝑗()𝑄𝑘−𝑙^𝑖+1,𝑗(𝑙)−2𝑄_𝑖,𝑗(𝑙)+𝑄_{𝑖−1,𝑗}(𝑙)Δ𝜃₂²+3^𝑘_𝑙=0𝐼_𝑖,𝑗(𝑘−𝑙)^𝑙_𝑚=0𝐻_{𝑖𝑖𝑖,𝑗}()⋅𝑄𝑙−𝑚_𝑖,𝑗+1(𝑚)−2𝑄_𝑖,𝑗(𝑚)+𝑄_{𝑖,𝑗−1}(𝑚)Δ𝜍₂²+^𝑘_𝑙=0𝐽_𝑖,𝑗()𝑄𝑘−𝑙^𝑖,𝑗+1(𝑙)−2𝑄_𝑖,𝑗(𝑙)+𝑄_{𝑖,𝑗−1}(𝑙)Δ𝜍₂²=2Λ^𝑘_{𝑖𝑖𝑙=0}𝑃_{𝑖𝑖𝑖,𝑗}()𝐻𝑘−𝑙^{𝑖𝑖𝑖+1,𝑗}(𝑙)−𝐻_{𝑖𝑖𝑖−1,𝑗}(𝑙)2Δ𝜃₂+2Λ^𝑘_{𝑖𝑖𝑙=0}𝐻_{𝑖𝑖𝑖,𝑗}()𝑃𝑘−𝑙^{𝑖𝑖𝑖+1,𝑗}(𝑙)−𝑃_{𝑖𝑖𝑖−1,𝑗}(𝑙)2Δ𝜃₂+4Λ_𝑖𝑖𝛾^𝑘_𝑙=0𝑘+1𝐻𝑃_{𝑖𝑖𝑖,𝑗}(𝑘−𝑙)𝐻_{𝑖𝑖𝑖,𝑗}(𝑙+1)+4Λ_𝑖𝑖𝛾^𝑘_𝑙=0𝑘+1𝐻𝐻_{𝑖𝑖𝑖,𝑗}(𝑘−𝑙)𝑃_{𝑖𝑖𝑖,𝑗}(𝑙+1), (13) where 𝐻 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 (𝑋₂₀,𝑌₂₀) 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 acceleration, velocity, and displacement can be estimated as 𝑑²𝑋𝑑𝜏²=𝐹_elx−𝐹_gfx𝑚_𝑟=𝐴𝑐_𝑋,𝑑²𝑌𝑑𝜏²=𝐹_ely−𝐹_gfy𝑚_𝑟=𝐴𝑐_𝑌,𝑉_𝑋=𝑉_𝑋𝑜+𝐴𝑐_𝑋𝑉⋅Δ𝜏,_𝑌=𝑉_𝑌𝑜+𝐴𝑐_𝑌⋅Δ𝜏,𝑋=𝑋_𝑜+𝑉_𝑋1⋅Δ𝜏+2𝐴𝑐_𝑋⋅Δ𝜏²,𝑌=𝑌_𝑜+𝑉_𝑌1⋅Δ𝜏+2𝐴𝑐_𝑌⋅Δ𝜏². (14)(15)(16)

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 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

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.

Table 1 

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

Tables 2 and 3 compare the Poincaré map data calculated by the hybrid method using different time step values, 𝐻, 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, 𝐻.

Table 2 

Comparison of Poincaré maps of rotor center with different values of 𝑚_𝑟 and 𝐻.

Table 3 

Comparison of Poincaré maps of rotor center with different values of Λ and 𝐻.

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

3.2.1. 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.27 kg. 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.

(a)

(b)

(c)

(d)

Figure 2 

Trajectories of rotor center at 𝑚𝑟=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).

3.2.2. Power Spectra

Figures 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.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 𝑚_𝑟 = 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

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.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≤𝑚_𝑟<13.179 kg. However, at 𝑚_𝑟=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≤𝑚_𝑟<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 has a 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) and 5(b).

(a)

(b)

Figure 3 

Bifurcation diagrams: (a) 𝑋2(𝑛𝑇) and (b) 𝑌2(𝑛𝑇) versus rotor mass 𝑚𝑟 at 𝜔=3.5×106 rpm.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 4 

Poincaré maps of rotor center trajectory at: (a) 𝑚𝑟=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)

Figure 5 

Maximum Lyapunov exponents of system at different values of rotor mass at: (a) 𝑚𝑟=14.74 (b) 15.12 kg.

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.

Table 4 

Behavior of rotor center at different rotor masses in interval 0.1≤𝑚_𝑟≤16.5kg.

3.3. 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.

3.3.1. 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.

(a)

(b)

(c)

(d)

Figure 6 

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 𝑚𝑟=3.0 kg).

3.3.2. 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.

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

Figure 7 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 𝑇-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𝑇→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.

(a)

(b)

Figure 7 

Bifurcation diagrams: (a) 𝑋2(𝑛𝑇) and (b) 𝑌2(𝑛𝑇) versus bearing number Λ over interval 1.0≤Λ≤7.8 at 𝑚𝑟=3.0 kg.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 8 

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)

Figure 9 

Maximum Lyapunov exponents of system at different values of bearing number at: (a) Λ=5.35 (b) 7.28.

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.

Table 5 

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

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

Acknowledgment

The financial support provided to this research by the National Science Council of China under Grant no. NSC 100-2628-E-269-016-MY2 is greatly appreciated.