Abstract
The bifurcation behaviors of a flexible rotor supported by gaslubricated bearing system with porous bushing are analyzed by a novel numerical method combining the finite difference method and differential transformation method. The results obtained by proposed method are verified with traditional finite difference method, and the analytical results by these two methods are consistent and also in good agreement. Furthermore, the dynamic orbits, power spectra, bifurcation diagram, Poincaré maps and maximum Lyapunov exponent are used to confirm the changes of rotor behavior as the rotor mass is increased. The results show that the rotor center reveals complex dynamic behaviors including periodic, subharmonic and chaotic motions. Especially, the rotor bearing system behaves chaos over the ranges of rotor mass 10.7 $\le {m}_{s}<$ 13.8 kg. The current results provide an effective means of the gas bearing systems and further understanding of the nonlinear dynamic behavior of bearing systems characterized by different rotor masses.
1. Introduction
Gaslubricated bearing system with porous bushing (GBSP) is different from other air bearing system due its particularly attractive features such as lownoise rotation, high load capacity, improved damping properties and zero friction. GBSP is simpler in construction and cheaper than the externally pressurized bearings. The comparison with oil bearings and GBSP reveals that GBSP generates less heat and provides higher accuracy. GBSP are also applied in high speed spindles, and machine tools.
In 1964, Sneck and Yen [1] considered onedimensional flow for radial direction in a porous medium and developed a perturbation solution for a finite journal bearing. Then, they verified theoretical solutions experimentally in 1965 [2].
In 2001, Sinha et al. [3] analyzed externally pressurized conical bearings with the porous constant gap when the slider rotated with uniform angular velocity. The system governing equations including coupled momentum and energy equations are analyzed by using finite difference method and various bearing characteristics are also determined. The results show that the inlet pressure decreases remarkably for highly porous surfaces, resulting in reduced load capacity of the bearing. But, the torque of system with respect to variation of permeability remains unaffected.
For the rotor dynamic study, the nonlinear dynamic behavior of rigid rotor supported by noncircular air journal bearing system including two, three and fourlobe bearings are solved by Rashidi et al. [4, 5] The Reynolds’ equation for air pressure distribution is analyzed by using finite element method and the rotor dynamic equations are calculated by RungeKutta method. The numerical results show that bifurcation phenomenon with different bearing numbers and rotor masses including periodic and nonperiodic motions at different operational situations. It’s found that bearing number and rotor mass are the major parameters for the bearing system.
In 2010 to 2013 [69], Wang et al. studied various kinds of gas bearing and analyzed the bifurcation and dynamics behaviors avoiding the occurrence of fatigue or unexpected situations. For example, as united gaslubricated bearings, the rotor center behaves a complex dynamic behavior under different values of the rotor mass and bearing number. The maximum Lyapunov exponent is also applied to distinguish different motions comprising periodic, subharmonic, quasiperiodic and chaotic responses. Furthermore, the dynamic behavior of noncircular air bearing system is studied by Wang et al. and a hybrid method including differential transformation method and finite difference method are proposed and used to solve the governing equations with various gas viscosity and rotational speed.
The remainder of this study is organized as follows. Section 2 develops a mathematical model of GBSP. Section 3 determines the Reynolds’ equation to obtain the gas pressure distribution by a hybrid method combining the finite difference method (FDM) and the differential transformation method (DTM). The numerical solutions are then compared with those obtained by traditional FDM. Section 4 presents the simulation results obtained using the proposed hybrid method for the vibrations of the rotor center for various rotor masses.
2. Governing equations
The GBSP system is used to form a closed cavity around the porous bushing with an appropriate casing and the porous material acts as a restrictor. The compressed gas is supplied and fed at a constant pressure, ${\widehat{P}}_{f}$, and the gas flows through the porous bushing into the bearing clearance. As the supply pressure drops to $\widehat{P}$, it exhausts into the atmosphere at pressure ${P}_{a}$. The GBSP model incorporates the ideal assumptions including the gas flow is assumed isothermal, the side flow of gas of the bearing is neglected, and the gas viscosity is constant. It is considered that a GBSP system comprising a perfectlybalanced flexible rotor supported symmetrically on two identical GBSPs mounting in turn on rigid pedestals. Due to the rotor is assumed as balancing perfectly and the GBSP is symmetric about its central axes, the system is confined to a single bearing supporting a rotor with two degrees of translatory oscillation in the transverse plane.
The gas pressure distribution between the rotor and the bushing is modeled by the nondimensionalized Reynolds’ equation as follows:
$=\frac{\partial}{\partial \theta}\left({\widehat{h}}^{3}\frac{\partial {\widehat{P}}^{2}}{\partial \theta}\right)+{\left(\frac{D}{L}\right)}^{2}{\widehat{h}}^{3}\frac{{\partial}^{2}{\widehat{P}}^{2}}{\partial {\widehat{z}}^{2}},$
where $\widehat{P}$ is the nondimensionalized pressure, ${\widehat{P}}_{f}$ is the nondimensionalized pressure in porous media, $\widehat{h}$ is the nondimensionalized thickness of gas film, $\mathrm{\Lambda}$ is the bearing number and $\theta $, $\widehat{z}$ are the nondimensionalized coordinates.
3. Simulation analysis by a novel numerical method
The simulation analysis is completed by two different numerical methods including the proposed method in this paper and traditional FDM. Accordingly, the FDM is applied for Eqs. (1) and (2) to discretize the $\theta $ and $\widehat{z}$ directions by the centraldifference scheme and then the time domain $\tau $ by the implicitbackdifference scheme. Also, a uniform mesh size is used. In order to compare the calculation results and increase the accuracy of system, a novel numerical method is proposed in this study. It is commenced by applying the DTM to discretize the Reynolds’ equation given in Eqs. (1) and (2) with respect to time and thus Eqs. (1) and (2) become as Eqs. (3)(6). Then, The FDM is used to discretize Eqs. (3)(6) with respect to the $\theta $ and $\widehat{z}$ directions. Note that Eqs. (3)(6) are discretized using the secondorderaccurate centraldifference scheme for both the first and the second derivatives. Substituting Eqs. (5)(6) into Eqs. (3)(4) yields Eqs. (7), (8):
$\begin{array}{c}4\varsigma \mathrm{\Lambda}\frac{\partial \varphi}{\partial \tau}\otimes \frac{\partial \widehat{h}}{\partial \theta}\otimes \widehat{P}+4\mathrm{\Lambda}\varsigma \frac{\partial \widehat{P}}{\partial \tau}\otimes \widehat{h}+4\varsigma \mathrm{\Lambda}\frac{\partial \widehat{h}}{\partial \tau}\otimes \widehat{P}\end{array}$
$=\frac{\partial J}{\partial \theta}\otimes \frac{\partial Q}{\partial \theta}+J\otimes \frac{{\partial}^{2}Q}{{\partial}^{2}\theta}+{\left(\frac{D}{L}\right)}^{2}J\otimes \frac{{\partial}^{2}Q}{\partial {\widehat{z}}^{2}},$
where:
The governing equations become:
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}2\varsigma \mathrm{\Lambda}\beta {\sum}_{l=0}^{k}\left(\frac{k+1}{\stackrel{~}{H}}\right)\left(\frac{{\left({\widehat{P}}_{f}\right)}_{i+1,j}\left(kl\right){\left({\widehat{P}}_{f}\right)}_{i1,j}\left(kl\right)}{2\mathrm{\Delta}\theta}\right){\varphi}_{i,j}\left(l+1\right)$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}=\left(\frac{{\left({\widehat{P}}_{f}\right)}_{i,j+1}\left(m\right)2{\left({\widehat{P}}_{f}\right)}_{i,j}\left(m\right)+{\left({\widehat{P}}_{f}\right)}_{i,j1}\left(m\right)}{{\left(\mathrm{\Delta}\widehat{z}\right)}^{2}}\right),$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+2\mathrm{\Lambda}{\sum}_{l=0}^{k}{\widehat{h}}_{i,j}\left(kl\right)\left(\frac{{\widehat{P}}_{i+1,j}\left(l\right){\widehat{P}}_{i1,j}\left(l\right)}{2\mathrm{\Delta}\theta}\right)+2\mathrm{\Lambda}{\sum}_{l=0}^{k}{\widehat{P}}_{i,j}\left(kl\right)\left(\frac{{\widehat{h}}_{i+1,j}\left(l\right){\widehat{h}}_{i1,j}\left(l\right)}{2\mathrm{\Delta}\theta}\right)$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}4\mathrm{\Lambda}\varsigma {\sum}_{l=0}^{k}{\widehat{h}}_{i,j}\left(kl\right){\sum}_{m=0}^{l}\left[\left(\frac{k+1}{\stackrel{~}{H}}\right)\left(\frac{{\widehat{P}}_{i+1,j}\left(m\right){\widehat{P}}_{i1,j}\left(m\right)}{2\mathrm{\Delta}\theta}\right){\varphi}_{i,j}\left(lm\right)\right]$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}4\mathrm{\Lambda}\varsigma {\sum}_{l=0}^{k}{\widehat{P}}_{i,j}\left(kl\right){\sum}_{m=0}^{l}\left[\left(\frac{k+1}{\stackrel{~}{H}}\right)\left(\frac{{\widehat{h}}_{i+1,j}\left(m\right){\widehat{h}}_{i1,j}\left(m\right)}{2\mathrm{\Delta}\theta}\right){\varphi}_{i,j}\left(lm\right)\right]$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}4\mathrm{\Lambda}\varsigma {\sum}_{l=0}^{k}\left[\left(\frac{k+1}{\stackrel{~}{H}}\right){\widehat{h}}_{i,j}\left(kl\right){\widehat{P}}_{i,j}\left(l+1\right)\right]+4\mathrm{\Lambda}\varsigma {\sum}_{l=0}^{k}\left[\left(\frac{k+1}{\stackrel{~}{H}}\right){\widehat{P}}_{i,j}\left(kl\right){\widehat{h}}_{i,j}\left(l+1\right)\right]$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}={\sum}_{l=0}^{k}\left(\frac{{J}_{i+1,j}\left(kl\right){J}_{i1,j}\left(kl\right)}{2\mathrm{\Delta}\theta}\right)\left(\frac{{Q}_{i+1,j}\left(l\right){Q}_{i1,j}\left(l\right)}{2\mathrm{\Delta}\theta}\right)$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+{\left(\frac{D}{L}\right)}^{2}{\sum}_{l=0}^{k}{J}_{i,j}\left(kl\right)\left(\frac{{Q}_{i,j+1}\left(l\right)2{Q}_{i,j}\left(l\right)+{Q}_{i,j1}\left(l\right)}{{\left(\mathrm{\Delta}\widehat{z}\right)}^{2}}\right).$
The calculation of iterative procedure can be integrated as follows. First, the new values of acceleration, velocity, and displacement of rotor are obtained by following a time increment $\mathrm{\Delta}\tau $. Second, from first step the rotor center displacements can then be determined and the corresponding change of the gas film gap ($\widehat{h}$) can be analyzed. Then the new value of $\widehat{h}$ is substituted into governing dequations to give the new pressure distribution of the gap between the shaft and the bushing. Third, the pressure distribution obtained from second step is integrated to calculate the internal force caused by the gaslubricated film. Finally, the displacement and velocity values calculated in the first step, the gas pressure computed in the second step, and the internal force analyzed in the third step are taken as the new initial conditions. Applying this new set of conditions, the calculation procedure returns to the first step to calculate the changes in the GBSP system during the time interval $\mathrm{\Delta}\tau \to 2\mathrm{\Delta}\tau $. These data including the orbital paths and velocity of the rotor center are used to generate bifurcation diagrams and Poincaré maps.
4. Numerical results
Table 1 presents the Poincaré maps obtained by the FDM and DTM&FDM for the orbits of the rotor center. It is observed that a good agreement by DTM&FDM method exists between the two sets of results at different rotor mass values. It also compares with different values of the time step. It can be seen that the orbits are in agreement to approximately 4 decimal places for the different time steps especially obtained by DTM&FDM method as shown in Table 1 and 2.
Table 1Comparison of rotor center orbits calculated by FDM and DTM&FDM methods, respectively
Displacement conditions  X2(nT)  
Time step: 0.001  Time step: 0.01  
FDM  ${m}_{r}=$7.2 kg $\omega =$1980 rad/s  0.008823820432  0.008961092741 
DTM&FDM  0.008825114313  0.008824802704  
FDM  ${m}_{r}=$10.1 kg $\omega =$1980 rad/s  0.059034858544  0.058062840042 
DTM&FDM  0.059136319514  0.059136084996 
Table 2Poincaré maps with different time increments and rotor mass by DTM&FDM
${m}_{r}=$ 7.2 kg  ${m}_{r}=$ 10.1 kg  
$\tau $  X2(nT)  Y2(nT)  $\tau $  X2(nT)  Y2(nT) 
$\pi $/300  0.008823943  –0.060185357  $\pi $/300  0.059136649  –0.075228092 
$\pi $/600  0.008823325  –0.060185087  $\pi $/600  0.059136640  –0.075228387 
Fig. 1Dynamic orbits of rotor center
a)${m}_{r}=$ 9.7 kg
b)${m}_{r}=\mathrm{}$9.96 kg
c)${m}_{r}=\mathrm{}$10.8 kg
d)${m}_{r}=\mathrm{}$13.5 kg
Figs. 1 and 2 show that the orbits and power spectra of the rotor center are regular at low values of the rotor mass (${m}_{r}=$9.7 and 9.96 kg), but become irregular at ${m}_{r}=$10.8 kg. This irregular behavior diverge as nonsymmetric and nonperiodic motion. When the rotor mass is increased at ${m}_{r}=$13.5 kg, the system converges to periodic motion.
Fig. 2Power spectra of rotor center
a)${m}_{r}=$ 9.7 kg
b)${m}_{r}=$9.96 kg
d)${m}_{r}=$10.8 kg
e)${m}_{r}=$13.5 kg
Fig. 3Bifurcation diagrams versus rotor mass
a) X2(nT)
b) Y2(nT)
Fig. 3 presents the bifurcation diagrams and shows the rotor center displacement of GBSP against the rotor mass. It is observed that qualitatively different behavior is varied at different values of ${m}_{r}$ within the range 6.0 to 14.0 kg. Fig. 4 present the Poincaré maps at ${m}_{r}=$ 9.7, 9.96, 10.8, and 13.5 kg, respectively. The dynamic motion of the rotor center is $T$periodic at lower values of the rotor mass, i.e. ${m}_{r}<$9.7 kg. However, the periodic motion loses its stability at ${m}_{r}=$9.7 kg and is replaced by subharmonic motion of 2$T$ as shown in Fig. 3 and Fig. 4(a). This subharmonic is persisted over the interval 9.7$\le {m}_{r}<$9.96, but when the mass is increased to ${m}_{r}=$9.96 kg, the subharmonic motion is changed to a $T$periodic motion in the $x$ and $y$directions shown in Fig. 3 and Fig. 4(b).
When the mass is increased to ${m}_{r}=$10.8 kg, system is changed to chaotic motion and maintained over the interval 10.8$\le {m}_{r}<$13.5 kg. The chaotic situation is shown in Fig. 3 and proved in Fig. 4(c). However, the chaotic motion transferred its motion at a rotor mass of 13.5 kg and is bifurcated into a 2$T$periodic motion shown in Fig. 3 and Fig. 4(d). Then, the subharmonic of 2$T$ motion is behaved over the interval 13.5$\le {m}_{r}\le $14.0 kg. It can be seen that two discrete points in the Poincaré maps at ${m}_{r}=$ 9.7 and 13.5 kg mean 2$T$ subharmonic motion and nonorder discrete points at ${m}_{r}=$10.8 kg reveal chaos.
In order to further verify the occurrence of chaotic motion, the maximum Lyapunov exponent is applied and it is shown that a positive value (as ${m}_{r}$ equals 10.8) presents chaos over the intervals 10.8$\le {m}_{r}<$13.5 kg shown in Fig. 5(a) and equals zero which means system behaves periodic motion. From the discussions above, it is evident that the behavior of the rotor center is dependent on the rotor mass. Table 3 summarizes the motions performed by the rotor center for rotor mass values in the range 6.0$\le {m}_{r}\le $14.0 kg.
Fig. 4Poincaré maps of rotor center trajectories
a)${m}_{r}=$ 9.7 kg
b)${m}_{r}=$9.96 kg
c)${m}_{r}=$10.8 kg
d)${m}_{r}=$13.5 kg
Fig. 5Maximum Lyapunov exponents of system at different values of rotor mass
a)${m}_{r}=$10.8 kg
b)${m}_{r}=$ 9.7 kg
Table 3Variation of rotor center response with rotor mass over interval 6.0 ≤mr≤ 14.0 kg
Rotor mass  [6.0, 9.7)  [9.7, 9.96)  [9.96, 10.8)  [10.8, 13.5)  [13.5, 14.0) 
Behavior  $T$  2$T$  $T$  Chaotic  2$T$ 
5. Conclusions
In this paper, the bifurcation behaviors of a flexible rotor supported by gaslubricated bearing system with porous bushing (GBSP) are analyzed and the chaotic behavior has also been found by utilizing a hybrid numerical scheme comprising the differential transformation method and the finite difference method in this paper. The system orbits, power spectra, bifurcation diagrams, Poincaré maps and maximum Lyapunov exponent have revealed the presence of a complex dynamic behavior comprising periodic, subharmonic and chaotic responses of the rotor center. From Table 1 the results obtained by the FDM and proposed hybrid method for the orbits of the rotor center prove that a good agreement exists between different sets of results. For the results of nonlinear dynamic motion, this paper provide a fully understanding of GBSP systems characterized by different rotor masses. Specifically, the results have shown that system exists chaotic motion over the ranges of rotor mass 10.8$\le {m}_{r}<$13.5 kg under rotational speed of 1980 rad/s.
References

Sneck H. J., Yen K. T. The externally pressurized, porous wall, gaslubricated journal bearing. Part 1. ASLE Transactions, Vol. 7, 1964, p. 288298.

Sneck H. J., Elwell R. C. The externally pressurized, porous wall, gaslubricated journal bearing. Part 2. ASLE Transactions, Vol. 8, 1965, p. 339345.

Sinha P., Chandra P., Bhartiya S. S. Thermal effects in externally pressurized porous conical bearings with variable viscosity. ACTA Mechanica, Vol. 149, 2001, p. 215227.

Rashidi R., Karami Mohammadi A., Bakhtiarinejad F. Preload effect on nonlinear dynamic behavior of a rigid rotor supported by noncircular gaslubricated journal bearing systems (three and four lobe). Nonlinear Dynamics, Vol. 60, 2009, p. 231253.

Rashidi R., Karami Mohammadi A., Bakhtiarinejad F. Bifurcation and nonlinear dynamic analysis of a rigid rotor supported by twolobe noncircular gaslubricated journal bearing system. Nonlinear Dynamics, Vol. 61, 2010, p. 783802.

Wang C. C., Kuo C. L. Nonlinear dynamic analysis of a relatively short spherical gas journal bearing system. Journal of Mechanical Science and Technology, Vol. 24, Issue 8, 2010, p. 15651571.

Wang C. C. Bifurcation and nonlinear dynamic analysis of united gaslubricated bearing system. Computers and Mathematics with Applications, Vol. 64, Issue 5, 2012, p. 729738.

Wang C. C., Wang C. C. Bifurcation and nonlinear dynamic analysis of noncircular aerodynamic journal bearing system. Nonlinear Dynamics, Vol. 72, Issue 1, 2013, p. 477489.

Wang C. C. Bifurcation analysis of bearing number in ultrashort gas bearing system. Smart Science, Vol. 1, Issue 1, 2013, p. 1824.
Cited by
About this article
The financial support of this research by Ministry of Science and Technology, under the Projects Nos. MOST 1032622E167015CC3 and 1032221E167011 are greatly appreciated.