Bifurcation and Nonlinear Dynamic Analysis of Externally Pressurized Double Air Films Bearing System

This paper studies the chaotic and nonlinear dynamic behaviors of a rigid rotor supported by externally pressurized double air films (EPDAF) bearing system. A hybrid numerical method combining the differential transformation method and the finite difference method is used to calculate pressure distribution of EPDAF bearing system and bifurcation phenomenon of rotor center orbits.The results obtained for the orbits of the rotor center are in good agreement with those obtained using the traditional finite difference approach.The results presented summarize the changes which take place in the dynamic behavior of the EPDAF bearing system as the rotor mass and bearing number are increased and therefore provide a useful guideline for the bearing system.


Introduction
Air bearing systems are characterized by low noise under rotation and by their low frictional losses.As a result, they are frequently employed within precision instruments, where they yield zero friction when the instruments are used as null devices, and within high-speed electrical motors.Compared with traditional oil bearings, air bearings have the advantages of lower heat generation, less contamination, and a higher precision.However, their major disadvantage is that they tend to be rather unstable, and this frequently restricts their permissible range of application [1,2].
For externally pressurized double air films (EPDAF) bearing system, it provides double films of lubrication, higher supporting force caused by gas pressure, and better performance than other air bearings.Air film pressure distribution is governed by Reynolds equation and the nonlinearity of analysis is concerned by researchers.In 1985, Gero and Ettles [3] evaluated the relative precision of the finite difference method (FDM) and finite element method (FEM) when applied to a steady, isoviscous, and incompressible lubrication problem.In their study, it was assumed that the solution of a complicated coupled problem could be derived by solving a sequential series of simple, uncoupled, and steady problems.The results for two-dimensional bearings demonstrated that the relative errors of the FDM solutions were smaller than those associated with the FEM approach.Furthermore, it was shown that the FDM approach was more rapid than the FEM technique, with an average CPU time of 0.15 sec as compared to 0.17 sec for the FEM method.
Because the rotational speed of the rotor could reach 10 6 rpm, the stability of rotor dynamics is one of the most important issues for considering the design of EPDAF bearing systems.How to increase the stability of rotor systems and control the appearance of nonperiodic motion became the major execution of this paper.Zhao et al. [4] investigated the subharmonic and quasiperiodic 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 subharmonic and quasiperiodic motions generated at speeds of more than twice the system critical speed bifurcated from an unstable harmonic solution.In 1997, Adiletta et al. [5,6] studied that a rigid rotor supported by short bearings behaved with subharmonic, quasiperiodic, or chaotic motion at particular values of the system parameters and performed these results by theoretical and experimental evidence.
According to the recent research, nonlinear dynamic responses of rotor-bearing systems are analyzed and published.In 2009, Rashidi et al. [7,8] studied the bifurcation and dynamic behavior of rigid rotor supported by noncircular aero-lubricated journal bearing system including two-, three-, and four-lobe bearings.Finite element method is applied to solve Reynolds' equation and the dynamic equations are analyzed by Runge-Kutta method.The results show that the bearing number and rotor mass are the major parameters for bifurcation phenomenon and the rotor behaviors including periodic and nonperiodic motions at different operational situations.For MENS, Zhou et al. [9] analyzed the bifurcation behavior of ultrashort self-acting gas journal bearings.The bearing system is modeled as a rigid rotor supported by bearing forces as a result of gas viscosity and rotational speed.In 2010, Wang [10] analyzed the bifurcation behavior and nonlinear dynamics of rigid rotor supported by microgas journal bearings and showed that the rotors exhibited a complex dynamic behavior comprising periodic, subharmonic, and quasiperiodic responses at different values of the rotor mass and bearing number, respectively.Wang [11] analyzed the nonlinear behavior of a high speed spindle air 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.
More and more high precision rotational mechanisms are applied in the process of manufacture and for high speed rotor are used and controlled in gas lubricated bearing system.So, the present study analyzes the nonlinear dynamic response of a rigid rotor supported by two externally pressurized double air films bearings.To analyze the EPDAF bearing system, the motions of the rotor center are governed by Reynolds' equations and solved by using a hybrid method combining the differential transformation method (DTM) and the finite difference method.The validity of the hybrid method is confirmed by comparing the results obtained for the rotor center orbits with those using the FDM scheme.The proposed method is then applied to analyze the bifurcation phenomenon and dynamic response of the rotor for rotor mass and bearing number in the ranges 1.0∼12.0kg and 1.0∼12.0,respectively.

Theoretical Modeling
2.1.Reynolds' Equations.An externally pressurized double air films bearing provides two films lubrication and is different from other aero-dynamic bearing systems.The application for high rotational speed more than 10 6 rpm can be used by EPDAF bearing system and the stability of this system is focused to be analyzed.In general, the pressure distribution in the air film between the rotor, floating-ring, and the bushing in an externally pressurized double air films bearing system is modeled using Reynolds' equation for an ideal gas shown in Figure 1.In the present study, the timedependent motions of the rotor center are modeled with the following.
The first film between rotor and floating-ring: Mathematical Problems in Engineering 3 The second film between floating-ring and bushing: where  is nondimensional pressure corresponding to the atmospheric pressure   ,   is nondimensional film thickness between rotor and bushing, corresponding to the radial clearance   ,  is diameter of bearing, and  is length of bearing.Λ is bearing number and  and  are coordinates in the circumferential and axial directions, respectively.The EPDAF bearing system model incorporates the design assumptions including the gas flow in and out of the sides of the bearing (side flow) which is neglected, the flow is isothermal, and the gas viscosity is constant.The air film pressure distribution must fulfill the following boundary conditions: The analysis presented in this study considers an EPDAF bearing system comprising a perfectly balanced rigid rotor of mass supported symmetrically on two identical EPDAF bearings, mounted in turn on rigid pedestals.Because the floating-ring is steady, externally pressurized provides gas for the first film when rotor is rotating, and these two films also can be separated for providing different gas pressures to control the stiffness, damping, and other performance coefficients.At the same time, the second film can be a flexible supporting pad and can obtain the optimum supporting effect.In real case, bearing exists eccentricity (  ⇀   ) and the forces caused by the two films acting on the rotor are ⃗  1 and ⃗  2 , respectively: where  1 and  2 are the stiffness of the first and second films, respectively;  1 and  2 are the damping coefficients of two films, respectively; ⃗  1 and ⃗  2 are the eccentricities of rotor and floating-ring, respectively; ⃗ ė 1 and ⃗ ė 2 are the velocity. 1 and  2 are the mass of rotor and floating-ring, respectively.
In the transient state, the equation of motion of the rotor can be written as (4)

Numerical Simulations.
In this study, the hybrid method is proposed and integrated FDM and differential transformation method (DTM) [12,13] to solve Reynolds' equation given in ( 1) and ( 2).DTM is one of the most widely used schemes for solving both linear and nonlinear differential equations due to its rapid convergence rate and minimal calculation error.
For solving Reynolds' equation in the current EPDAF bearing system, DTM is used to transform Reynolds' equation with respect to the time domain , and thus ( 1) and ( 2) become as follows.
The first film between rotor and floating-ring: The second film between floating-ring and bushing: where The finite difference method (FDM) is then applied to discretize ( 5) and ( 6) with respect to the  and  directions.Note that ( 5) and ( 6) are discretized using the second-orderaccurate central-difference scheme for both the first and the second derivatives.Substituting ( 7) into ( 5) and ( 6) yields the following.
The first film between rotor and floating-ring: The second film between floating-ring and bushing: where H is time step value.From ( 8) and ( 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 bearing, and the velocity of the rotor is assumed to be zero.
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 EPDAF bearing system over representative ranges of the rotor mass and bearing number, respectively.

Results and Discussions
3.1.Numerical Results.The orbits of the rotor center are compared to obtain the results by the FDM and hybrid method (FDM and DTM) shown in Table 1.It is observed that a good agreement exists between the two sets of results.However, it can be seen that the FDM suffers numerical instability at specific values of the rotor mass and time step and the hybrid method converges under all the considered conditions and therefore represents a more appropriate method for analyzing the dynamic response of the EPDAF bearing system.In Tables 2 and 3, the Poincaré map data calculated by the hybrid method using different time step values, H, are compared for different rotor mass and bearing number values, respectively.For a given rotor mass and bearing number, the rotor center orbits are in agreement with approximately 4 decimal places for the different time steps, H.

Dynamic Analysis with Increasing Rotor Mass.
A constant rotational speed with a gradually increasing rotor mass of the EPDAF bearing system is considered.The gas bearing was loaded with a constant rotational speed of 6 × 10 6 rpm and the rotor mass  1 was varied in the range 1.0 ≤  1 ≤ 12.0 kg.
The results show that the dynamic orbits and phase trajectories of the rotor center are regular at a low value of the rotor mass ( 1 = 2.2 kg) but become irregular at a rotor mass of  1 = 8.17 and 8.22 kg as shown in Figures 2(a) and 2(b) (2.1-2.12),respectively.For rotor mass values of  1 = 8.29 and 9.1 kg, the orbits transfer to periodic characteristic.Then increasing rotor mass causes irregular and nonperiodic motion especially when  1 is greater than 10.08 kg.
Figures 2(c) and 2(d) (2.1-2.12)show the power spectra 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  1 = 2.2 kg.When the rotor mass is increased to 8.17 and 8.22 kg, the spectral response graphs (Figures 2(c) and 2(d) (2.2-2.3))displayed that the rotor center behaves quasiperiodic motion.As the rotor mass is increased to 8.29 and 9.1 kg, the system turns into a T-periodic motion.As  1 = 10.08 kg, system orbit is mutated to chaotic motion and then changes to subharmonic motion when rotor mass increased to 10.7 kg.This type of movement converted back into chaos at  1 = 10.75 kg.
Thereafter, when rotor mass is increased to 11.23, 11.27, 11.66, and 11.68 kg, the system shows multiperiodic motion and chaotic motions alternately.But the system still dominated by chaotic motion ultimately as  1 is greater than 11.68 kg.
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  1 .Meanwhile, Figures 4(a)-4(j) present the Poincaré maps of the rotor center trajectories at  1 = 2.2, 8.17, 8.22, 8.29, 9.1, 10.08, 10.7, 10.75, 11.23, 11.27, 11.66, and 11.68 kg, 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,  1 < 8.17 kg.This is confirmed by the Poincaré map shown in Figure 4(a) for a rotor mass of 2.2 kg.However, at  1 = 8.17 kg, the T-periodic motion is replaced by quasiperiodic motion (see Figure 4(b)).From Figure 3, it is seen that this quasiperiodic motion persists over the rotor mass range 8.17 ≤  1 < 8.18 kg .However, at  1 = 8.18 kg,          the maximum Lyapunov exponent has a positive value, which indicates that the system has a chaotic response, shown in Figures 5(a) and 5(b).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 1.0 ≤  1 ≤ 12.0 kg.

Dynamic Analysis with Increasing Bearing Number.
In the second situation, the rotor mass is specified as   = 2.3 kg and the bearing number Λ was increased over the range 1.0 ≤ Λ ≤ 12.0.Figures 6(a) and 6(b) (6.1-6.6)show that the dynamic orbits and phase trajectories of the rotor center are regular at a low value of the bearing number (Λ = 3.2, 5.6) but become irregular at Λ = 9.05 and 9.18.When the bearing number is changed to 9.31 and 11.2, the rotor center performs periodic motion.It is also observed that the results between dynamic orbits and phase trajectories are consistent at different bearing numbers.
From Figures 6(c) and 6(d) (6.1-6.6), the rotor center performs a periodic motion at a bearing number of Λ = 3.2 and 5.6.However, when the bearing number is changed to Λ = 9.05 and 9.18, the power spectra show that the rotor center performs quasiperiodic motion in both the horizontal and the vertical directions.For values of Λ equal to 9.31 and 11.2, the rotor center performs periodic motion.
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 12.0 are shown in Figure 7.And in order to further confirm the rotor behavior, the Poincaré maps of the rotor center trajectories at Λ = 3.2, 5.6, 9.05, 9.18, 9.31, and 11.2 are presented in Figures 8(a)-8(f).The results show that the rotor center performs T-periodic motion over the bearing number range 1.0 ≤ Λ < 9.05 and proved by Figures 8(a) and 8(b).The T-periodic motion becomes unstable at bearing numbers of Λ = 9.05 and is replaced by quasiperiodic motion (see Figure 8(c)).This quasiperiodic behavior persists over 9.05 ≤ Λ < 9.31 (see Figures 8(c) and 8(d)).This quasiperiodic motion is transferred to periodic motion at Λ = 9.31 (see Figure 8(e)) and behaves in the range 9.31 ≤ Λ ≤ 12.0 proved by Figure 8(f).Figure 9 shows that the maximum Lyapunov exponent has a zero value when Λ equals 9.05 and 9.18 and also indicates that the system has a periodic response.Meanwhile, over the intervals 1.0 ≤ Λ < 9.05 and 9.31 ≤ Λ ≤ 12.0, shown in  Figures 7 and 9, the system behaves in stable behavior and also should be operated under these parameters.These bearing number intervals, which correspond to the typical operating conditions of an EPDAF bearing system, are characterized by T-and quasiperiodic motions of the rotor center.However, at specific and higher bearing numbers (9.31 ≤ Λ), the rotor performs predominantly periodic 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 ≤ Λ ≤ 12.0.

Conclusions
The bifurcation and nonlinear dynamic behavior of an externally pressurized double air films bearing system has been analyzed by utilizing a hybrid numerical scheme comprising the differential transformation method and the finite difference method in this study.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, quasiperiodic, and chaotic responses of the rotor center.Although this hybrid method might be not the best candidate, from Table 1 the results obtained by the FDM and hybrid method for the orbits of the rotor center prove that a good agreement exists between different sets of results.Moreover, the hybrid method converges under all the considered conditions and therefore represents a higher accuracy and appropriate method for analyzing the dynamic response of the EPDAF bearing system.

Figure 5 ,Figure 5 :
Figure 5: Maximum Lyapunov exponents of system at different values of rotor mass at (a) 1 = 10.08 and (b) 11.27 kg.

Table 1 :
Comparison of rotor center orbits for steady floating-ring calculated by FDM and hybrid method, respectively.

Table 2 :
Comparison of Poincaré maps of rotor center for steady floating-ring with different values of  1 and H by hybrid method.

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

Table 4 :
Behavior of rotor center at different rotor masses in interval 1.0 ≤  1 ≤ 12.0 kg.

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