Nonlinear analysis and simulation of active hybrid aerodynamic and aerostatic bearing system

Abstract

Active hybrid aerodynamic and aerostatic bearing system has gradually become more valued in recent years for its application in the precision machine field, especially for precision instruments and mechanisms that require high rotational speed, high precision, and high rigidity support. In this system, air lubricated bearing is mainly used for support. Although the carrying capacity is not as high as oil film bearings, air lubricated bearings can provide a work environment where the rotor (spindle) will not experience axial deformation at high rotational speed and low heat generation. In practice, spindle system dynamic problems include critical speed, spindle imbalance, and improper bearing design, which can cause the spindle system to produce aperiodic motion, instability, and even chaotic motion under certain parameters and conditions during its operation. If severe, these irregular motions may cause machine damage or delay production. In order to understand what type of work situation can produce aperiodic phenomenon, and to prevent irregular vibration and reduce instantaneous air hammer effect, the finite difference method and mixed method are used to explore the related characteristics for spindle system. Meanwhile, relevant theories including bifurcation diagram, Poincare map, spectrum response phenomenon, and maximum Lyapunov exponents are applied to analyze rotor non-linear dynamic behavior. In addition, we verified the non-linear motion parameter conditions to prevent spindle system from falling into the instable status during design. The objective is to lower the probability of chaotic phenomenon in the system and reduce damage caused by irregular vibrations in the system. The result of this study can serve as a reference when designing precision spindle systems or mechanisms.

Introduction

Active hybrid aerodynamic and aerostatic bearing (AHAAB) system is different than general bearing system. The rotation speed of general rotors (spindle) is only between 14,000 and 24,000 r/min, and has a DmN value (the spindle high-speed indicator, which is related to bearing material, structure design, carrying capacity, precision, and critical speed) upper limit of 2 × 10⁶. However, the spindle rotational speed of AHAAB can reach more than 3 × 10⁵ r/min and a DmN value greater than 2 × 10⁶. This is mainly because gas bearing has no contact; therefore, has very little friction. The bearing is almost operating in a no friction status. Thus, there is little wear between the rotor and the bearing, which results in longer spindle system life.^1,2 Gas is also compressible, has very low value of viscosity, which means it has very little fluid shear and does not require much power to achieve very high rotation speed. This means that the bearing does not heat up easily, and can operate in a high temperature environment without affecting its dynamic performance. Because the lubricant used by this type of spindle system is air, it can also operate at extremely high rotation speed, high temperature, or low temperature. The homogenizing effect of gas can effectively compensate for processing deviation on the working surface and also significantly improve the rotational and translation accuracy of spindle system. Thus, gas bearing can serve as a circularity reference, straight line reference, and as a main machine component for ultra-precision processing technology and micro-processing technology. As such, it has received the attention of domestic and international academic fields and industries. The production and manufacturing standard of gas bearings represent the research level of humans in the ultra-precision machining technology field. Currently, air lubricated bearings have replaced large number of traditional rolling component or liquid film lubricated bearings in instruments and equipment that require high precision axis positioning, high rotational speed, low friction, and have other special environmental requirements. ³

Studies on the performance of bearing and rotor (spindle) related dynamic characteristics of air lubricated spindle system show that the static and dynamic characteristics of high-speed rotor (spindle) directly affect the precision quality of the work piece. Therefore, when designing an spindle system, not only does the torsion or bending strength of the spindle must be calculated, but the static and dynamic stiffness should also be carefully designed. In air lubricated bearing studies, how to accurately calculate air film pressure distribution has always been an important topic in analyzing air lubricated bearings. In 1991, Nguyen^4,5 proposed using finite element method (FEM) to solve complex border gas bearing pressure distribution problems. The AHAAB used for this study have both externally supplied pressure and special groove or floating ring designs that provided superior support and stability. However, there are few studies conducted on this subject. For groove basic theory, scholars have used FEM to solve the pressure distribution of herringbone groove bearing and studied the effects produced by groove design. This includes stability analysis, pressure distribution, and the impact of critical parameters on the system. Study results showed that spindle systems with grooves can achieve superior stability.^6–8 For air lubricated bearings with externally supplied pressure, holes setup above the bearing can provide appropriate gas pressure and also produce greater rigidity. Consequently, the rotor (spindle) can operate at a higher rotation speed.9,10 The aforementioned reference shows that the air film dynamic characteristics has a significant impact on the overall bearing support system, especially for the air film pressure, stiffness, and damping effect. In 1974, Bently ¹¹ discovered that rotor-bearing system has 2nd order and 3rd order subharmonic vibration by experiments. Choi and Noah ¹² proved that rotor-bearing system has subharmonic vibrations by numerical simulation. Although these references conducted studies on oil film bearing, very few studies were conducted on the rotor behavior of air lubricated bearing systems.

For the analysis methods of nonlinear behaviors, a new method called the Exp-function method was first proposed by He and Wu ¹³ and was successfully applied to a KdV equation with variable coefficients ¹⁴ and to high-dimensional nonlinear evolution equations. This method was also applied to obtain generalized solitary solutions of the generalized Drinfel’d–Sokolov–Wilson (DSW) system. Then, some of the solitary solutions are converted to periodic solutions or hyperbolic function solutions by a simple transformation. The research results show that the Exp-function method is a powerful and convenient mathematical tool for solving nonlinear evolution equations with higher order nonlinearity. ¹⁵ In 2003, a new perturbation method was proposed and can analyze the rotor dynamic equation efficiently. In contrast to the traditional perturbation methods, this technique does not require a small parameter for the rotor equation. According to the homotopy technique, a homotopy with an imbedding parameter is considered and replaced the small parameter, so this method can be called the homotopy perturbation method. ¹⁶ In 2008, He ¹⁷ proposed He’s frequency formulation for analyzing the nonlinear oscillators and also presented the application of variational principle.^18–20 In 2017, He studied the well-known Hamilton’s principle and found that it is valid for the conditions prescribed at the beginning and at the end of the motion. ²¹ For AHAAB system, the initial position of rotor is known and can be applied the Hamilton’s principle to obtain the orbits caused by external forces and different operation conditions.

In recent years, Wang et al.^22–25 used mixed method to solve the air film pressure distribution of high-speed rotor air lubricated bearing with grooves. Meanwhile, the system state trajectories, Poincaré maps, power spectra, bifurcation diagrams and Lyapunov exponents were used to analyze the dynamic behavior of the rotor and journal centers. In this study, we discovered that under different operation conditions, the rotor and journal centers have periodic, quasi-periodic, subharmonic or nonperiodic motions. It is also found that within a specific operating interval, the system may even produce chaotic motion.

System modeling and research method

Dynamic analysis of active hybrid aerodynamic and aerostatic bearing system

In AHAAB system, the interior of air lubricated bearings often have excess gas capacity (GC), and this GC generally has lower natural frequency. When certain interference source appears inside the system, and its frequency comes close to the natural frequency of the bearing systems, excitation will occur, causing resonance in the excess GC. This phenomenon is similar to air hammer striking an object, and is called air hammer vibration phenomenon. This phenomenon often causes bearing damage and rotor instability. The main reasons of air hammer are as follows:

Internal GC of a certain size: If the bearing’s total internal GC is ${\displaystyle {\sum }_1^N{V}_i}$ , of which $V_i$ is each GC’s volume and N is the total GC content. The total air film volume of bearing is $V_T$ , and $V_r=\left({\displaystyle {\sum }_1^N{V}_i}\right)/V_T$ shows the corresponding GC. In terms of flat thrust bearing, air hammer effect will appear when $V_r\ge \left(2\sim 10\%\right)$ . In terms of radial bearing, air hammer will occur when $V_r\ge \left(5\sim 15\%\right)$ . Thus, the size of the internal GC must be considered and controlled when bearing system is designed.

Existence of an excitation source: Air hammer vibration is a self-excited vibration meaning an oscillation that occurs from the system’s internal excitation by internal gas pressure or volume changes. Thus, regardless of whether the system is static or rotating, or whether the pressure is normal or high pressure, air hammer still occurs. When the external interference frequency and internal excitation frequency is closed, a resonance phenomenon can be triggered.

The distribution of air film pressure: The bearing’s internal pressure is determined by externally supplied pressure and bearing rotation speed. The higher the internal pressure, the more severe the air hammer vibration. When the rotation speed of rotor is closed or equal to the system GC’s vibration frequency, the air hammer can be exacerbated. The air hammer vibration always exists at critical rotation speed.

To overcome and consider the aforementioned three factors, we analyzed the pure air guided type of AHAAB systems and described below.

Pure air guided type bearing mainly utilizes air guided method shown in Figure 1, as well as integrated bearing surface groove design. This bearing type is divided into single-way valve forward and reverse air guided type, all groove air guided type, and partial groove air guided type. Because this type of system has the advantage of both static pressure and dynamic pressure bearings, it can effectively inhibit air hammer effect and improve system stability.

Figure 1.

Pure air guided type bearing (a) with a single way valve, (b) with continued grooves and (c) with partial grooves.

This study mainly analyzes groove air guided bearing system. This spindle system has both dynamic and static pressure effects, and its dynamic pressure effects come from the groove design on the bearing surface. The static pressure effects come from the configuration of the air (Ps) from air guide holes. When analyzing the air film pressure distribution (P), designers must consider the Reynolds equation (1), as follows:

$$\frac{\partial }{\partial z}\left(h^{3}P\frac{\partial P}{\partial z}\right)+\frac{\partial }{\partial \theta }\left(h^{3}P\frac{\partial P}{\partial \theta }\right)={\Lambda }_h\frac{\partial \left(Ph\right)}{\partial \theta }+\frac{\partial \left(Ph\right)}{\partial t}$$

(1)

where non-dimensional air film thickness is $h=1+\epsilon \hspace{0.17em}\cos \hspace{0.17em}\theta$ ; $\epsilon$ is the eccentricity; ${\Lambda }_h=\frac{12\mu \omega }{\left(P_s+P_a\right)}{\left(\frac{R_o}{h_o}\right)}^2$ ; $R_o$ is the rotor radius; P_a is the atmospheric pressure. The boundary conditions of AHAAB must satisfy equation (2):

$$\{\begin{array}{ll}z=L/D,\hfill & P=1\hfill \\ z=0,\hfill & \frac{\partial P}{\partial z}=-f\left(\xi \right)\frac{{\Lambda }_s{P}_{s^2}}{z{h}^3{P}_i}\sqrt{\frac{1+{\delta }^2}{1+{\left(\delta /h\right)}^2}}\hfill \end{array}$$

(2)

And also be satisfied periodic condition equation (3):

$$\{\begin{array}{l}P\left(\theta ,z\right)=P\left(\theta +2\pi ,z\right)\hfill \\ {\frac{\partial P}{\partial \theta }|}_{\theta ,z}={\frac{\partial P}{\partial \theta }|}_{\theta +2\pi ,z}\hfill \end{array}$$

(3)

where $\{\begin{array}{cc}\frac{P_i}{P_s}\le {\left(\frac{2}{\xi +1}\right)}^{\frac{1}{\xi -1}},& f\left(\xi \right)=C_f\sqrt{\frac{2\xi }{\xi +1}}{\left(\frac{2}{\xi +1}\right)}^{\frac{1}{\xi -1}}\\ 1\ge \frac{P_i}{P_s}>{\left(\frac{2}{\xi +1}\right)}^{\frac{1}{\xi -1}},& f\left(\xi \right)=C_f\sqrt{\frac{2\xi }{\xi +1}}{\left(\frac{P_i}{P_s}\right)}^{\frac{1}{\xi }}\sqrt{1-{\left(\frac{P_i}{P_s}\right)}^{\frac{\xi -1}{\xi }}}\end{array}$ ; $\delta =\frac{r_a^2}{h_o{d}_a}$ ; $d_a$ is the diameter of air hole; ${\Lambda }_s=\frac{6\mu N{r}_{a^2}}{P_s{h}_o^3}\sqrt{\frac{RT}{1+{\delta }^2}}$  ; $r_a$ is the hole radius of inside section flow, μ is the gas viscosity coefficient, h_o is the average air film thickness, R is the gas constant, T is the temperature, and $C_f$ is the gas flow coefficient. Equation (4) is the neck shrinkage effect comprehensive corrective coefficient for gas flow friction and wear at the exit hole. This is mainly determined by the hole shape, pressure ratio η (proportion of throttle end exit pressure and provided air pressure, η=P_d/P_s), and the Reynolds number. The aforementioned parameter design can change the size of the GC inside the bearing and reduce the air hammer effect.

$$C_f=\{\begin{array}{ll}0.8+0.0171\eta ,\hfill & \eta <0.7\hfill \\ {\left[0.26-{\left(\eta -0.4\right)}^2\right]}^{0.5}+0.4,\hfill & 0.7\le \eta <0.9\hfill \\ -5\left(\eta -1\right)\hfill & 0.9\le \eta \le 1\hfill \end{array}$$

(4)

Numerical method was used to calculate the pressure distribution of the air film. Integration was conducted on lubrication area to obtain the relative carrying capacity:

$$\begin{array}{c}{F}_x^{=}\frac{1}{2_{\left(L/D\right)}}{\displaystyle \underset{\iint }{\hspace{0.17em}}{P}_{\cos }{}_{\theta }\text{d}}\theta \text{d}\varsigma ,\\ F_y^{=}\frac{1}{2\left(L/D\right)}{\displaystyle \underset{\iint }{\hspace{0.17em}}{P}_{\sin }\theta \text{d}}{\hspace{0.17em}}_{\theta }\text{d}\varsigma \end{array}$$

(5)

The air film total carrying capacity W is

$$W_{=}\sqrt{{\left({\hspace{0.17em}}_F^x\right)}^2+{\left({\hspace{0.17em}}_F^y\right)}^2}$$

(6)

To analyze the dynamic parameter of the air lubricated bearing system, the rotor mass is considered as m_r and supported in the gas bearing. At the same time, the bearing is fixed on a rigid base. In the transient state, the equations of motion of O₂ (i.e. the rotor center) and O₃ (i.e. the journal center) can be written in Cartesian coordinates as:

$$m_r\frac{d^2{x}_2}{d{t}^2}+K_p(x_2-x_3)=m_r\rho {\omega }^2\cos \omega t$$

(7)

$$m_r\frac{d^2{y}_2}{d{t}^2}+K_p(y_2-y_3)=m_r\rho {\omega }^2\sin \omega t$$

(8)

where ρ is the mass eccentricity of the rotor, K_p is the shaft stiffness, and ω is the rotational speed of the rotor.

The main advantage of the air guided type mixed dynamic/static pressure air lubricated bearing system is its high-speed stability, which is superior to general air film bearing. Thus relevant coefficients in equations (7) and (8) can be applied to confirm the stability of the AHAAB system. The motion behavior of rotor bearing system was also analyzed by the system state trajectories, Poincaré maps, power spectra, bifurcation diagrams and maximum Lyapunov exponents.

Proposed numerical method

In order to solve the Reynolds equation efficiently, we used two different numerical methods to analyze the bearing system, including traditional finite difference method (FDM) and proposed mixed method. These two sets of numerical results by different methods were confirmed and compared. The proposed mixed method mainly combines differential transformation method (DTM) and FDM. First, DTM was used to discretize time in the non-linear Reynolds equation. FDM central differential processing was used on position coordinates. The following describes the DTM in detail:

The differential transformation of the function of f(t) is defined as:

$$\overline{f}\left(k\right)=\frac{T^k}{k!}{\left[\frac{{\partial }^{k}f\left(t\right)}{\partial t^k}\right]}_{t=0}$$

(9)

where $\overline{f}\left(k\right)$ is the transformation function in the transformation domain. $f\left(t\right)$ is the original function in the time domain, and T is the time gap. k is the transformation parameter. The inverse transformation equation for $\overline{f}\left(k\right)$ is:

$$f\left(t\right)={\displaystyle \sum _{k=0}^{\infty }\overline{f}\left(k\right){\left(\frac{t}{T}\right)}^k}$$

(10)

This method also requires confirmation of accuracy during calculation. That is, when time increases to the second interval, then the following equations (11) and (12) shall be used for accuracy calculation.

$$f_{i-1}\left(0\right)=f_i\left(0\right)+{{\displaystyle \sum _{k=1}^{\infty }{\overline{f}}_i\left(k\right)\left(\frac{-T_{i-1}}{T_i}\right)}}^k$$

(11)

$$\text{If}{T}_o=T_1=T_2=\dots ,\text{then}\text{the}\text{above}\text{equation}\text{is}{f}_{i-1}\left(0\right)=f_i\left(0\right)+{\displaystyle \sum _{k=1}^{\infty }{\overline{f}}_i\left(k\right){\left(-1\right)}^k}$$

(12)

If the calculated value on the two sides of equation in equation (11) and (12) is too different, that means the accuracy was insufficient, and the time interval T must be reduced and the calculation redone.

The differential transformation basic algorithm is as shown in Table 1.

Table 1.

Differential transformation basic calculation formula.

Original function	Transformation function
$af\left(t\right)\pm bg\left(t\right)$	$a\overline{f}\left(k\right)\pm b\overline{g}\left(k\right)$
$f\left(t\right)g\left(t\right)$	${\displaystyle \sum _{l=0}^k\overline{f}\left(k-l\right)}\overline{g}\left(l\right)={\displaystyle \sum _{l=0}^k\overline{f}\left(l\right)}\overline{g}\left(k-l\right)$
$f\left(t\right)/g\left(t\right)$	$\frac{\overline{f}\left(k\right)-{\displaystyle \sum _{l=0}^{k-1}\left(\overline{f}\left(l\right)/\overline{g}\left(l\right)\right)\overline{g}\left(k-l\right)}}{\overline{g}\left(0\right)}=\frac{\overline{f}\left(k\right)-{\displaystyle \sum _{l=0}^{k-1}\left(\overline{f}\left(k-l\right)/\overline{g}\left(k-l\right)\right)\overline{g}\left(l\right)}}{\overline{g}\left(0\right)}$
$\cos \left(at+b\right)$ , $a\ne 0$	$\frac{{\left(aT\right)}^k}{k!}\cos \left(\frac{\pi k}{2}+b\right)$
$\frac{d^{n}f\left(t\right)}{d{t}^n}$	$\frac{\left(k+n\right)!}{k!T^n}\overline{f}\left(k+n\right)$
$f^m\left(t\right)$	${\displaystyle \sum _{l=0}^k\overline{f}\left(k-l\right)}{\displaystyle \sum _{n=0}^l{\hspace{0.17em}}_{\overline{f}}{\hspace{0.17em}}^{m-2}\left(n\right)}\overline{f}\left(l-n\right)$

For AHAAB system, the pressure distribution function can then be calculated by differential transformation with boundary conditions. The pressure function is divided into n number of subintervals and each subinterval is divided into T₀, T₁, T₂, T₃, …, then equation (10) shall show the function in the first interval as:

$$P_o\left(t\right)={\displaystyle \sum _{k=0}^{\infty }{\overline{P}}_o\left(k\right){\left(\frac{t}{T_o}\right)}^k}{,}_{\hspace{1em}}0\le t\le T_o$$

(13)

If initial value $P_o\left(0\right)$ is known, the above equation is used to obtain $P_o\left(0\right)={\overline{P}}_o\left(0\right)$ . The ${\overline{P}}_o\left(0\right)$ can be used to obtain other discrete value ${\overline{P}}_o\left(k\right){,}_{\hspace{0.17em}\hspace{0.17em}}k=1,2,\mathrm{\dots }$ in the interval. In the first subinterval, the end value of function $P_o\left(t\right)$ is $P_o\left(t=T_o\right)$ , which is the initial value $P_1\left(t=0\right)$ of the second subinterval function $P_1\left(t\right)$ . The second subinterval function $P_1\left(t\right)$ is:

$$P_1\left(t\right)={\displaystyle \sum _{k=0}^{\infty }{\overline{P}}_1\left(k\right){\left(\frac{t}{T_1}\right)}^k},0\le t\le T_1$$

(14)

Using ${\overline{P}}_1\left(0\right)$ can similarly obtain other ${\overline{P}}_1\left(k\right){,}_{\hspace{0.17em}\hspace{0.17em}}k=1,2,\mathrm{\dots }$ , and so on to obtain the ${\overline{P}}_i\left(k\right)$ value. Then, inverse differential transformation can be used to obtain its various time item range solution. Its series solution form is as follows:

$$P_i\left(t\right)={\displaystyle \sum _{k=0}^m{\overline{P}}_i\left(k\right){\left(\frac{t}{T_i}\right)}^k},0\le t\le T_i$$

(15)

The aforementioned numerical method can be used with traditional perturbation method to conduct accuracy confirmation. This, combined with comparison of numerical results from follow-up experiments, can further verify the critical performance coefficient.

In solving the Reynolds equation for the current gas bearing system, the DTM is used for taking transformation with respect to the time domain τ, and hence equation (1) becomes

$$\begin{array}{c}3I\otimes \frac{\partial H}{\partial \theta }\otimes \frac{\partial S}{\partial \theta }+J\otimes \frac{{\partial }^{2}S}{\partial {\theta }^2}+3I\otimes \frac{\partial H}{\partial \phi }\otimes \frac{\partial S}{\partial \phi }+J\otimes \frac{{\partial }^{2}S}{\partial {\phi }^2}=2\Lambda H\otimes \frac{\partial P}{\partial \theta }+2\Lambda P\otimes \frac{\partial H}{\partial \theta }\\ +4\Lambda H\otimes \frac{\partial P}{\partial \tau }+4\Lambda P\otimes \frac{\partial H}{\partial \tau }\end{array}$$

(16)

where

$$S(k)=P^2=P\otimes P={\displaystyle \sum _{l=0}^k{P}_{i,j}(k-l)P_{i,j}(l)}$$

(17)

$$I(k)=H^2=H\otimes H={\displaystyle \sum _{l=0}^k{H}_{i,j}(k-l)H_{i,j}(l)}$$

(18)

$$J(k)=H^3=H\otimes H\otimes H={\displaystyle \sum _{l=0}^k{H}_{i,j}(k-l){\displaystyle \sum _{m=0}^l{H}_{i,j}(l-m)}{H}_{i,j}(m)}$$

(19)

The FDM is then used to discretize (12) with respect to the θ and $\phi$ directions. Note that (16) is discretized using the second-order-accurate central-difference scheme for both the first and the second derivatives.

Substituting (17) to (19) into (16) yields

$$\begin{array}{c}3{\displaystyle \sum _{l=0}^k{I}_{i,j}(k-l){\displaystyle \sum _{m=0}^l\left(\frac{H_{i+1,j}(l-m)-H_{i-1,j}(l-m)}{2\Delta \theta }\right)}\cdot \left(\frac{S_{i+1,j}(m)-S_{i-1,j}(m)}{2\Delta \theta }\right)}\\ +{\displaystyle \sum _{l=0}^k{J}_{i,j}(k-l)\left(\frac{Q_{i+1,j}(l)-2Q_{i,j}(l)+Q_{i-1,j}(l)}{{\left(\Delta \phi \right)}^2}\right)}\\ \times 3{\displaystyle \sum _{l=0}^k{I}_{i,j}(k-l){\displaystyle \sum _{m=0}^l\left(\frac{H_{i,j+1}(l-m)-H_{i,j-1}(l-m)}{2\Delta \phi }\right)\cdot \left(\frac{S_{i,j+1}(m)-P_{i,j-1}(m)}{2\Delta \phi }\right)}}\\ +{\displaystyle \sum _{l=0}^k{J}_{i,j}(k-l)\left(\frac{Q_{i,j+1}(l)-2Q_{i,j}(l)+Q_{i,j-1}(l)}{{\left(\Delta \phi \right)}^2}\right)}=2\Lambda \cdot {\displaystyle \sum _{l=0}^k\left[H_{i,j}(k-l)\left(\frac{P_{i+1,j}(l)-P_{i-1,j}(l)}{2\Delta \theta }\right)\cdot \right]}\\ +2\Lambda {\displaystyle \sum _{l=0}^k\left[P_{i,j}(k-l)\left(\frac{H_{i+1,j}(l)-H_{i-1,j}(l)}{2\Delta \theta }\right)\cdot \right]}+4\Lambda {\displaystyle \sum _{l=0}^k\left[\left(\frac{l+1}{\tilde{H}}\right)P_{i,j}(k-l)\cdot H_{i,j}(l+1)\right]}\\ +4\Lambda {\displaystyle \sum _{l=0}^k\left[\left(\frac{l+1}{\tilde{H}}\right)H_{i,j}(k-l)\cdot P_{i,j}(l+1)\right]}\end{array}$$

(20)

From (20), P_i,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 (X_o,Y_o) is specified as the static equilibrium position of the shaft and defines the gap H_i,j(k) between the shaft and the journal bearing, and the velocity of the rotor is assumed to be zero.

Calculation steps are as follows:

Step 1: When time $\tau =0$ , use initial condition to define initial rotor displacement. After $\Delta \tau$ time, calculate the rotor’s new acceleration, speed, and displacement.

Step 2: Substitute the rotor center displacement obtained in step 1 into the rotor equilibrium equation to obtain the rotor center’s new displacement. As rotor changes the air film thickness H, the new displacement value forms the new H. Then substitute back the differential produced by the Reynolds equation to calculate the new pressure value.

Step 3: Use integration on pressure distribution obtained in step 2 which can produce new air film force, rigidity coefficient, and damping coefficient.

Step 4: Use the displacement and speed obtained in step 1, pressure distribution obtained in step 2, and air film force obtained in step 3 to form a new set of initial conditions. When time increases by $\Delta \tau$ , use the aforementioned new initial condition to once again calculate the change of system.

The aforementioned four steps can obtain the pressure distribution, air film thickness to position change correlation, change in air film support force, rigidity coefficient, damping coefficient, and rotor dynamic orbits with time increased. The maximum Lyapunov exponents are used as the main judgment criteria to predict periodic, quasi-periodic, aperiodic, and chaotic for different affecting factors such as bearing number, and rotor mass. High-speed air lubricated bearing’s rotation speed and rotor mass are used as the bifurcation parameter to construct the bifurcation diagram and analyze the dynamic stability of AHAAB system. We also analyzed the relationship between the maximum Lyapunov exponents and the key factors, as well as established the stable and unstable areas for different operating parameter of AHAAB system.

Result and discussion

Numerical analysis of bearing system

We used two different numerical methods to solve this bearing system for pure air guided type AHAAB system. Comparison shows that there is a difference in accuracy between the mixed method (DTM combined with FDM) and traditional FDM. Numerical test with different rotor mass at the same time step showed that the mixed method has superior accuracy of numerical simulation, and has good consistency with traditional FDM analysis results. Thus, the mixed method proposed in this study had superior convergence for this system. The comparison of rotor center orbits is shown in Tables 2 and 3.

Table 2.

Comparison of rotor center’s Poincare sectional value at different time interval τ and different rotor mass (results obtained with mixed method calculation).

mr =3.1 kg			mr =4.6 kg
τ	X2 (n T)	Y2 (n T)	τ	X2 (n T)	Y2 (n T)
π/300	−0.05826172659	0.03106521457	π/300	−0.08461595662	0.05876383874
π/600	−0.05821042909	0.03105903518	π/600	−0.08461366813	0.05878327213

Table 3.

Comparison of rotor center’s Poincare sectional value at different time interval τ and different bearing number (results obtained with mixed method calculation).

Λ = 3.5			Λ = 7.1
τ	X2 (n T)	Y2 (n T)	τ	X2 (n T)	Y2 (n T)
π/300	−0.0592348521	0.1728183665	π/300	0.2521848954	−0.0299801717
π/600	−0.0592552394	0.1728325412	π/600	0.2521451824	−0.0299723613

For mixed method stability numerical analysis, we tested the effects of different time intervals on values, as well as used Poincare sectional value to compare displacement changes under different operating conditions. Results are shown in Tables 2 and 3.

The numerical results show that the mixed method used in this study has superior convergence for AHAAB system and is more suitable for this system. The mixed method can effectively calculate the track of system regardless of rotor mass change or increase in bearing number. Tables 2 and 3 results show that time interval does not need to be overly meticulously to obtain sufficiently precise numerical results. Thus, for follow-up dynamic analysis, π/300 is used for the time interval calculation.

Rotor dynamic behavior analysis of AHAAB system

Effect of rotor mass on the system

The bearing number (rotation speed) is set to be Λ = 3.45, and rotor mass $m_r$ is used as the bifurcation parameter:

Phase diagram and dynamic trajectory diagram

Fig. 2(2.1(a) and (b), 2.2(a) and (b), … 2.12(a) and (b)) shows that the rotor center and journal center with less mass ( $m_r$ =3.7kg) appear to be regular. When the mass is increased to $m_r$ =9.14 kg, aperiodic phenomenon replaces regular motion. When the rotor mass reaches $m_r$ = 9.15 kg, dynamic track turns again to periodic motion. As the mass continue to be increased, the system exhibited periodic than aperiodic phenomenon. When the mass is 11.29 kg, aperiodic effects gradually expanded. When at $m_r$ = 12.67 kg, the bearing system then maintains an unstable and aperiodic behavior.

Figure 2.

Phase diagram (2.1a to 2.12a), journal center dynamic track (2.1b to 2.12b), spectrum response of rotor center in the horizontal direction (2.1c to 2.12c), and in the vertical direction (2.1d to 2.12d) when the rotor mass at the rotor center is $m_r$ = 3.7, 9.14, 9.15, 9.18, 9.2, 9.25, 9.28, 11.29, 11.99, 12.04, 12.56, 12.67 (system bearing number Λ = 3.45).

• b. Spectrum analysis

Fig. 2(2.1(c) and (d), 2.2(c) and (d), …, 2.12(c) and (d)) shows the dynamic effects of the rotor center at the horizontal and vertical directions. Study shows that when the rotor mass is at $m_r$  = 3.7 kg, the rotor center exhibits a T-periodic motion. At $m_r$ =9.14 kg, spectrum response figures (Fig. 2(2.2(c) and (d))) show that quasi-periodic motion occurs at the horizontal and vertical directions of rotor center. As rotor mass changes to $m_r$  = 11.29 kg, the motion alternates between T-periodic motion and quasi-periodic motion. As mass increases to $m_r$  = 12.04 and 12.67 kg, the dynamic track then changes to chaotic motion.

c. Bifurcation diagram

As shown in Figures 3 and 4, we used rotor mass $m_r$ as the analysis parameter to explore the bifurcation phenomenon in horizontal and vertical directions, respectively and analyzed the effects of different rotor mass on the system. We also set the range of rotor mass from 0.01 to 13.0 kg. Figures 3 and 4 show that when the rotor mass is smaller, the system is more stable. At the same time, Figure 3 shows that when $m_r$ <9.14 kg, the system rotor center exhibited periodic motion in both the horizontal and vertical directions. This motion behavior can be further proved with the Poincare sectional diagram (Figure 5(a)). However, this periodic motion of rotor center is replaced by quasi-periodic motion when the mass increases to $m_r$ =9.14 kg, as shown in Figure 5(b). This motion status is maintained within the interval of 9.14 $\le m_r$  < 9.15 kg. However, as mass increases to $m_r$ =9.15 kg, the quasi-periodic motion rapidly changed to T-periodic motion. This status can be further proved by Figure 5(c).

Figure 3.

Bifurcation diagram when the rotor center in the horizontal direction X2(nT) with different rotor mass; system bearing number Λ = 3.45.

Figure 4.

Bifurcation diagram of the rotor center in the vertical direction Y2(nT) with different rotor mass; system bearing number Λ = 3.45.

Figure 5.

Poincare sectional diagram when the rotor is at different rotor mass (a) $m_r$ = 3.7, (b) 9.14, (c) 9.15, (d) 9.18, (e) 9.2, (f) 9.25, (g) 9.28, (h) 11.29, (i) 11.99, (j) 12.04, (k) 12.56, and (l) 12.67 kg (system bearing number Λ = 3.45).

When the mass increased to $m_r$ = 9.18 kg, the system once again changes to quasi-periodic motion, as shown in Figure 5(d). This motion is only maintained within the interval of 9.18 $\le m_r$  < 9.2 kg. As mass continued to increase to $m_r$ =9.2 kg, the system motion already exhibits intermittent vibration. This vibration occurs between the T and quasi-periodic motion (at 9.2 $\le m_r$ <11.99 kg), as shown in Figure 5(e) to (h). When $m_r$ =11.99 kg, the original T-periodic motion changes to multiple T-periodic (sub-harmonic) motion. This motion is only maintained within the range of 11.99 $\le m_r$ <12.04 kg. Subsequently, the motion exhibited by the system shows chaos and 8 T vibration. This motion is only maintained within the interval of 12.04 $\le m_r\le$ 13.0 kg. When $m_r$ =12.04 and 12.67 kg, the system changes to chaotic motion, as proved by Figure 5(j) and (l). As $m_r$ =12.56 kg, system shows 8 T-motion, as proven by Figure 5(k). The results show the rotor motion type also exhibits a certain level of cycle sequence. Table 4 shows the motion form exhibited by the system at different rotor mass and also summarizes the findings.

Table 4.

Rotor center dynamic behavior for different rotor mass (system bearing number Λ = 3.45).

Rotor mass	[0.01,9.14]	[9.14,9.15]	[9.15,9.18]	[9.18,9.2]	[9.2,9.23]
Dynamic behavior	T	Quasi	T	Quasi	T
Rotor mass	[9.23,9.24]	[9.24,9.25]	[9.25,9.28]	[9.28,11.29]	[11.29,11.99]
Dynamic behavior	Quasi	T	Quasi	T	Quasi
Rotor mass	[11.99,12.04]	[12.04,12.56]	[12.56,12.59]	[12.59,12.6]	[12.6,12.65]
Dynamic behavior	Multi-T	Chaos	8T	Chaos	8T
Rotor mass	[12.65,12.66]	[12.66,12.67]	[12.67,13.0]
Dynamic behavior	Chaos	8T	Chaos

Figure 5 shows that when $m_r$ =3.7, 9.15, and 9.2 kg, the motion exhibited by the bearing is a single point on the Poincare sectional diagram. This means that the system has a single periodic motion. Quasi-periodic motion occurs when the motion status is a single closed curve composed of multiple discrete points. Figure 5(j) and (l) is irregular discrete point, which means the system is in chaotic motion.

To confirm chaotic behavior, we used maximum Lyapunov exponents for interpretation and further verification. Figure 6 shows that when the maximum Lyapunov exponent is close to zero, then the system is in non-chaotic behavior. When $m_r$ =12.04 and 12.67, as shown in Figure 6(j) and (l), the index is greater than zero. Thus, the system motion is in chaotic behavior.

Figure 6.

Maximum Lyapunov exponents of the rotor center at different rotor mass (a) $m_r$ = 3.7, (b) 9.14, (c) 9.15, (d) 9.18, (e) 9.2, (f) 9.25, (g) 9.28, (h) 11.29, (i) 11.99, (j) 12.04, (k) 12.56, and (l) 12.67 kg (system bearing number Λ = 4.5).

Effect of bearing number on the system

Because the rotational speed of rotor related to the bearing number, it affects the performance of the AHAAB system, the bearing number (Λ) is used as the bifurcation parameter and the rotor mass is set at m_r=2.7 kg to analyze the system dynamic behavior:

Phase diagram and dynamic trajectory diagram

Figure 7(7.1(a) and (b), 7.2(a) and (b), …, 7.6(a) and (b)) shows that with smaller bearing number (Λ = 2.3, 4.5), the rotor exhibited regular behavior. When Λ was increased to 5.65, the system exhibited irregular motion. When bearing number reached 5.82, the dynamic orbit changed to regular motion. When Λ = 5.97 and 6.4, the rotor center’s behavior became more stable and regular.

Figure 7.

Phase diagram (7.1a to 7.6a), journal center dynamic trajectory diagram (7.1b to 7.6b), spectrum response of rotor center in the horizontal direction (7.1c to 7.6c), and in the vertical direction (7.1d to 7.6d) when the rotor center bearing number is Λ = 2.3, 4.5, 5.65, 5.82, 5.97, 6.4 (rotor mass set at m_r=2.7 kg).

• b. Spectrum analysis

Figure 7(7.1(c) and (d), 7.2(c) and (d), …, 7.6(c) and (d)) shows the dynamic effect of different bearing number on the rotor center in the horizontal and vertical directions. Results show that when the bearing number is at Λ = 2.3 and 4.5, the rotor center exhibits a T-periodic motion. When Λ = 5.65, spectrum response figure (Figure 7(7.3(c) and (d))) shows the rotor center of horizontal and vertical directions both has aperiodic and chaotic motion. As bearing number increased to Λ = 5.82, the motion changes to 3 T-periodic motion. At Λ = 5.97 and 6.4, the system changes to T-periodic motion. This is the same in both horizontal and vertical directions.

c. Bifurcation diagram

As shown in Figures 8 and 9, we used different bearing number Λ as the primary analysis parameter to explore the effect of different bearing number Λ on the system. At the same time, we set the bearing number between 1.0 and 7.0 for the actual operating condition. At low bearing number Λ = 2.3 and 4.5, the system rotor center shows T-periodic motion in both vertical and horizontal directions. This phenomenon is verified by Figure 10(a) and (b). When bearing number is increased to Λ = 5.65, system changes to chaotic motion. Figure 10(c) clearly shows that at this time, the system produced multiple uneven discrete points. This chaotic motion is maintained only when 5.65 $\le$ Λ< 5.82. When bearing number is increased to Λ = 5.82, the system changes to 3 T-periodic motion. This status can be further verified by Figure 10(d) (three discrete points). This motion status is maintained at the interval of 5.82≤Λ < 5.97. As Λ increased to 5.97, the bearing system gradually shows stable T-periodic motion, as shown in Figure 10(e) and (f).

Figure 8.

Bifurcation diagram when the rotor center in the horizontal direction X2(nT) with different bearing number; rotor mass set at m_r=2.7 kg.

Figure 9.

Bifurcation diagram when the rotor center in the vertical direction Y2(nT) with different bearing number; rotor mass set at m_r=2.7 kg.

Figure 10.

Poincare sectional diagram when the rotor center is at different bearing number (a) Λ = 2.3, (b) 4.5, (c) 5.65, (d) 5.82, (e) 5.97, and (f) 6.4 (rotor mass set at m_r=2.7 kg).

The relationship between the rotor center’s motion and bearing number is very complex. Thus, we organized 1.0≤Λ ≤ 7.0 interval in Table 5, as shown.

Table 5.

Dynamic behavior of rotor center with different bearing number (1.0≤Λ ≤ 7.0, rotor mass set at m_r=2.7 kg).

Bearing number Λ	[1.0,5.65]	[5.65,5.82]	[5.82,5.97]	[5.97,7.0]
Dynamic behavior	T	Chaos	3T	T

For chaotic behavior confirmation, we used maximum Lyapunov exponents for interpretation and verification. Figure 11 shows that when maximum Lyapunov exponents are less than zero, the system exhibits non-chaotic behavior. When Λ = 5.65 (as shown in Figure 11(c)), the index is greater than zero, and the system motion is in chaotic behavior.

Figure 11.

Maximum Lyapunov exponents of the rotor center at different bearing number (a) Λ = 2.3, (b) 4.5, (c) 5.65, (d) 5.82, (e) 5.97, and (f) 6.4 (rotor mass set at m_r=2.7 kg).

Conclusions

This study has analyzed the dynamic behavior of an AHAAB system by utilizing a mixed numerical scheme comprising the DTM and the FDM. The system state trajectories, Poincaré maps, power spectra, bifurcation diagrams and Lyapunov exponents have revealed the presence of a complex dynamic behavior comprising periodic, sub-harmonic, quasi-periodic and chaotic responses of the rotor center. Meanwhile, the mixed method is the best candidate, and from Table 6 the results obtained by the FDM and DTM&FDM (mixed method) methods for the orbits of the rotor center prove that a good agreement exists between different sets of results. Moreover, the mixed method converges under all the considered conditions and therefore represents a more appropriate method for analyzing the nonlinear dynamic response of the AHAAB system. The results of this study provide an understanding of the nonlinear dynamic behavior of AHAAB systems characterized by different rotor masses and bearing numbers. Specifically, the results have shown that chaos appears in four rotor mass intervals, namely 12.04–12.56, 12.59–12.6, 12.65–12.66 and 12.67–13.0 kg. Regarding the influence of the bearing number on the dynamic response of the bearing system, the results have shown that the chaotic motion only occurs in the bearing number range of 5.65 $\le$ Λ<5.82. It is important to avoid these unstable conditions and guarantee the AHAAB system to maintain under stable situation. The results can also be the guideline for designing the AHAAB system in the future.

Table 6.

Comparison of various numerical calculation results for pure air guided type bearing system rotor center displacement.

DisplacementNumerical method and operating conditions			X2 (nT)		Y2 (nT)
			$\tilde{H}=0.001$	$\tilde{H}=0.01$	$\tilde{H}=0.001$	$\tilde{H}=0.01$
1	FDM	mr =3.7 kgΛ=2.3	−0.06087137902	−0.06045043725	0.03886728132	0.03849926128
	DTM and FDM		−0.06054473319	−0.06054351081	0.03865198187	0.03862333781
2	FDM	mr =6.9 kgΛ=2.3	−0.11980291451	−0.11902020847	0.11035462781	0.11005245447
	DTM and FDM		−0.11959102184	−0.11959114787	0.11027545894	0.11027703665
3	FDM	mr =2.4 kgΛ=3.7	0.03702541797	0.03710735332	0.1353711961	0.1351064247
	DTM and FDM		0.03733738069	0.03733607179	0.1352456436	0.1352297949
4	FDM	mr =3.1 kgΛ=3.7	0.37453014551	0.37403026441	0.2172084618	0.2178334053
	DTM and FDM		0.37426403842	0.37426014958	0.2171939163	0.2171950693

FDM: finite difference method; DTM: differential transformation method.