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

Authors: Cheng-Chi Wang; Chin-Chia Liu; Chi-Chang Wang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

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.


---

## Body

## 1. 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 106 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.0 kg and 1.0~12.0, respectively.


## 2. 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 106 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 time-dependent motions of the rotor center are modeled with the following.Cross-section of externally pressurized floating-ring air bearing system.
(a)
(b)The first film between rotor and floating-ring:(1)
∂
2
∂
θ
1
[
H
i
3
∂
P
i
2
∂
θ
1
]
+
∂
2
∂
ς
1
[
H
i
3
∂
P
i
2
∂
ς
1
]
-
Λ
i
∂
(
P
i
H
i
)
∂
θ
1
-
Λ
i
Ω
∂
(
P
i
H
i
)
∂
θ
1
-
2
Λ
i
γ
∂
(
P
i
H
i
)
∂
τ
=
0
.The second film between floating-ring and bushing:(2)
∂
2
∂
θ
2
[
H
i
i
3
∂
P
i
i
2
∂
θ
2
]
+
∂
2
∂
ς
2
[
H
i
i
3
∂
P
i
i
2
∂
ς
2
]
-
Λ
i
i
∂
(
P
i
i
H
i
i
)
∂
θ
2
-
2
Λ
i
i
γ
∂
(
P
i
i
H
i
i
)
∂
τ
=
0
,
where P is nondimensional pressure corresponding to the atmospheric pressure P
a, H
i is nondimensional film thickness between rotor and bushing, corresponding to the radial clearance C
r, D is diameter of bearing, and L 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:(a)
floating-ring is steady asω
2
=
0, and air pressure on both ends of the housing is equal to the atmospheric pressure P
a, P
(
θ
,
±
ς
)
=
1,
(b)
air pressure is an even function forς, P
(
θ
,
ς
)
=
P
(
θ
,
-
ς
), and is continuous at ς = 0, ∂
P
∂
ς
|
ς
=
0
=
0,
(c)
air pressure is a periodic function forθ, P
(
θ
,
ς
)
=
P
(
θ
+
2
π
,
ς
), ∂
P
∂
θ
|
θ
=
∂
P
∂
θ
|
θ
+
2
π.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 (e
i
⃑) and the forces caused by the two films acting on the rotor are F
→
r
1 and F
→
r
2, respectively:
(3)
F
→
r
1
=
C
s
1
(
1
-
ω
1
2
Ω
)
(
e
˙
→
1
-
e
˙
→
2
)
+
K
s
1
(
e
→
1
-
e
→
2
)
,
F
→
r
2
=
C
s
2
e
˙
→
1
+
K
s
2
e
→
2
,
where K
s
1 and  K
s
2 are the stiffness of the first and second films, respectively; C
s
1 and  C
s
2 are the damping coefficients of two films, respectively; e
→
1and  e
→
2 are the eccentricities of rotor and floating-ring, respectively; e
˙
→
1 and e
˙
→
2 are the velocity. m
1 and m
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)
m
1
e
¨
→
1
+
F
→
r
1
=
0
,
m
2
e
¨
→
2
+
F
→
r
2
-
F
→
r
1
=
0
,
v
v
v
v
v
v
v
v
v
v
v
v
v
m
1
e
¨
→
1
+
C
s
1
(
1
-
ω
1
2
ω
s
)
(
e
˙
→
1
-
e
˙
→
2
)
+
K
s
1
(
e
→
1
-
e
→
2
)
=
0
,
m
2
e
¨
→
2
+
C
s
2
e
˙
→
2
+
K
s
2
e
→
2
+
C
s
1
(
1
-
ω
1
2
ω
s
)
×
(
e
˙
→
2
-
e
˙
→
1
)
+
K
s
1
(
e
→
2
-
e
→
1
)
=
0
.


### 2.2. 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:(5)
3
2
T
⊗
∂
H
i
∂
θ
1
⊗
∂
S
∂
θ
1
+
1
2
U
⊗
∂
2
S
∂
θ
1
2
+
3
2
T
⊗
∂
H
i
∂
ς
1
⊗
∂
S
∂
ς
1
+
1
2
U
⊗
∂
2
S
∂
ς
1
2
=
Λ
i
P
i
⊗
∂
H
i
∂
θ
1
+
Λ
i
H
i
⊗
∂
P
i
∂
θ
1
+
Λ
i
Ω
P
i
⊗
∂
H
i
∂
θ
1
+
Λ
i
Ω
H
i
⊗
∂
P
i
∂
θ
1
+
2
Λ
i
γ
P
i
⊗
∂
H
i
∂
τ
+
2
Λ
i
γ
H
i
⊗
∂
P
i
∂
τ
.The second film between floating-ring and bushing:(6)
3
2
T
⊗
∂
H
i
i
∂
θ
2
⊗
∂
S
∂
θ
2
+
1
2
U
⊗
∂
2
S
∂
θ
2
2
+
3
2
T
⊗
∂
H
i
i
∂
ς
2
⊗
∂
S
∂
ς
2
+
1
2
U
⊗
∂
2
S
∂
ς
2
2
=
Λ
i
i
P
i
i
⊗
∂
H
i
i
∂
θ
2
+
Λ
i
i
H
i
i
⊗
∂
P
i
i
∂
θ
2
+
2
Λ
i
i
γ
P
i
i
⊗
∂
H
i
i
∂
τ
+
2
Λ
i
i
γ
H
i
i
⊗
∂
P
i
i
∂
τ
,
where
(7)
S
(
k
)
=
P
2
=
P
⊗
P
=
∑
l
=
0
k
P
i
,
j
(
k
-
l
)
P
i
,
j
(
l
)
,
T
(
k
)
=
H
2
=
H
⊗
H
=
∑
l
=
0
k
H
i
,
j
(
k
-
l
)
H
i
,
j
(
l
)
,
U
(
k
)
=
H
3
=
H
⊗
H
⊗
H
=
∑
l
=
0
k
H
i
,
j
(
k
-
l
)
∑
m
=
0
l
H
i
,
j
(
l
-
m
)
H
i
,
j
(
m
)
.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-order-accurate 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:(8)
3
2
∑
l
=
0
k
T
i
,
j
(
k
-
l
)
∑
m
=
0
l
(
H
i
)
i
,
j
(
l
-
m
)
·
(
S
i
+
1
,
j
(
m
)
-
2
S
i
,
j
(
m
)
+
S
i
-
1
,
j
(
m
)
(
Δ
θ
1
)
2
)
+
1
2
∑
l
=
0
k
U
i
,
j
(
k
-
l
)
(
S
i
+
1
,
j
(
l
)
-
2
S
i
,
j
(
l
)
+
S
i
-
1
,
j
(
l
)
(
Δ
θ
1
)
2
)
+
3
2
∑
l
=
0
k
T
i
,
j
(
k
-
l
)
∑
m
=
0
l
(
H
i
)
i
,
j
(
l
-
m
)
v
v
v
v
v
·
(
S
i
,
j
+
1
(
m
)
-
2
S
i
,
j
(
m
)
+
S
i
,
j
-
1
(
m
)
(
Δ
ς
1
)
2
)
+
1
2
∑
l
=
0
k
U
i
,
j
(
k
-
l
)
(
S
i
,
j
+
1
(
l
)
-
2
S
i
,
j
(
l
)
+
S
i
,
j
-
1
(
l
)
(
Δ
ς
1
)
2
)
=
Λ
i
∑
l
=
0
k
(
P
i
)
i
,
j
(
k
-
l
)
(
(
H
i
)
i
+
1
,
j
(
l
)
-
(
H
i
)
i
-
1
,
j
(
l
)
2
Δ
θ
1
)
+
Λ
i
∑
l
=
0
k
(
H
i
)
i
,
j
(
k
-
l
)
(
(
P
i
)
i
+
1
,
j
(
l
)
-
(
P
i
)
i
-
1
,
j
(
l
)
2
Δ
θ
1
)
+
Λ
i
Ω
∑
l
=
0
k
(
P
i
)
i
,
j
(
k
-
l
)
×
(
(
H
i
)
i
+
1
,
j
(
l
)
-
(
H
i
)
i
-
1
,
j
(
l
)
2
Δ
θ
1
)
+
Λ
i
Ω
∑
l
=
0
k
(
H
i
)
i
,
j
(
k
-
l
)
(
(
P
i
)
i
+
1
,
j
(
l
)
-
(
P
i
)
i
-
1
,
j
(
l
)
2
Δ
θ
1
)
+
2
Λ
i
γ
∑
l
=
0
k
[
(
k
+
1
H
~
)
(
P
i
)
i
,
j
(
k
-
l
)
(
H
i
)
i
,
j
(
l
+
1
)
]
+
2
Λ
i
γ
∑
l
=
0
k
[
(
k
+
1
H
~
)
(
H
i
)
i
,
j
(
k
-
l
)
(
P
i
)
i
,
j
(
l
+
1
)
]
.The second film between floating-ring and bushing:(9)
3
2
∑
l
=
0
k
T
i
,
j
(
k
-
l
)
∑
m
=
0
l
(
H
i
i
)
i
,
j
(
l
-
m
)
·
(
S
i
+
1
,
j
(
m
)
-
2
S
i
,
j
(
m
)
+
S
i
-
1
,
j
(
m
)
(
Δ
θ
2
)
2
)
+
3
2
∑
l
=
0
k
U
i
,
j
(
k
-
l
)
(
S
i
+
1
,
j
(
l
)
-
2
S
i
,
j
(
l
)
+
S
i
-
1
,
j
(
l
)
(
Δ
θ
2
)
2
)
+
3
2
∑
l
=
0
k
T
i
,
j
(
k
-
l
)
∑
m
=
0
l
(
H
i
i
)
i
,
j
(
l
-
m
)
·
(
S
i
,
j
+
1
(
m
)
-
2
S
i
,
j
(
m
)
+
S
i
,
j
-
1
(
m
)
(
Δ
ς
2
)
2
)
+
1
2
∑
l
=
0
k
U
i
,
j
(
k
-
l
)
(
S
i
,
j
+
1
(
l
)
-
2
S
i
,
j
(
l
)
+
S
i
,
j
-
1
(
l
)
(
Δ
ς
2
)
2
)
=
Λ
i
i
∑
l
=
0
k
(
P
i
i
)
i
,
j
(
k
-
l
)
(
(
H
i
i
)
i
+
1
,
j
(
l
)
-
(
H
i
i
)
i
-
1
,
j
(
l
)
2
Δ
θ
2
)
+
Λ
i
i
∑
l
=
0
k
(
H
i
i
)
i
,
j
(
k
-
l
)
(
(
P
i
i
)
i
+
1
,
j
(
l
)
-
(
P
i
i
)
i
-
1
,
j
(
l
)
2
Δ
θ
2
)
+
2
Λ
i
i
γ
∑
l
=
0
k
[
(
k
+
1
H
~
)
(
P
i
i
)
i
,
j
(
k
-
l
)
(
H
i
i
)
i
,
j
(
l
+
1
)
]
+
2
Λ
i
i
γ
∑
l
=
0
k
[
(
k
+
1
H
~
)
(
H
i
i
)
i
,
j
(
k
-
l
)
(
P
i
i
)
i
,
j
(
l
+
1
)
]
,whereH
~ is time step value.From (8) and (9), 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 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.


## 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 106 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 time-dependent motions of the rotor center are modeled with the following.Cross-section of externally pressurized floating-ring air bearing system.
(a)
(b)The first film between rotor and floating-ring:(1)
∂
2
∂
θ
1
[
H
i
3
∂
P
i
2
∂
θ
1
]
+
∂
2
∂
ς
1
[
H
i
3
∂
P
i
2
∂
ς
1
]
-
Λ
i
∂
(
P
i
H
i
)
∂
θ
1
-
Λ
i
Ω
∂
(
P
i
H
i
)
∂
θ
1
-
2
Λ
i
γ
∂
(
P
i
H
i
)
∂
τ
=
0
.The second film between floating-ring and bushing:(2)
∂
2
∂
θ
2
[
H
i
i
3
∂
P
i
i
2
∂
θ
2
]
+
∂
2
∂
ς
2
[
H
i
i
3
∂
P
i
i
2
∂
ς
2
]
-
Λ
i
i
∂
(
P
i
i
H
i
i
)
∂
θ
2
-
2
Λ
i
i
γ
∂
(
P
i
i
H
i
i
)
∂
τ
=
0
,
where P is nondimensional pressure corresponding to the atmospheric pressure P
a, H
i is nondimensional film thickness between rotor and bushing, corresponding to the radial clearance C
r, D is diameter of bearing, and L 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:(a)
floating-ring is steady asω
2
=
0, and air pressure on both ends of the housing is equal to the atmospheric pressure P
a, P
(
θ
,
±
ς
)
=
1,
(b)
air pressure is an even function forς, P
(
θ
,
ς
)
=
P
(
θ
,
-
ς
), and is continuous at ς = 0, ∂
P
∂
ς
|
ς
=
0
=
0,
(c)
air pressure is a periodic function forθ, P
(
θ
,
ς
)
=
P
(
θ
+
2
π
,
ς
), ∂
P
∂
θ
|
θ
=
∂
P
∂
θ
|
θ
+
2
π.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 (e
i
⃑) and the forces caused by the two films acting on the rotor are F
→
r
1 and F
→
r
2, respectively:
(3)
F
→
r
1
=
C
s
1
(
1
-
ω
1
2
Ω
)
(
e
˙
→
1
-
e
˙
→
2
)
+
K
s
1
(
e
→
1
-
e
→
2
)
,
F
→
r
2
=
C
s
2
e
˙
→
1
+
K
s
2
e
→
2
,
where K
s
1 and  K
s
2 are the stiffness of the first and second films, respectively; C
s
1 and  C
s
2 are the damping coefficients of two films, respectively; e
→
1and  e
→
2 are the eccentricities of rotor and floating-ring, respectively; e
˙
→
1 and e
˙
→
2 are the velocity. m
1 and m
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)
m
1
e
¨
→
1
+
F
→
r
1
=
0
,
m
2
e
¨
→
2
+
F
→
r
2
-
F
→
r
1
=
0
,
v
v
v
v
v
v
v
v
v
v
v
v
v
m
1
e
¨
→
1
+
C
s
1
(
1
-
ω
1
2
ω
s
)
(
e
˙
→
1
-
e
˙
→
2
)
+
K
s
1
(
e
→
1
-
e
→
2
)
=
0
,
m
2
e
¨
→
2
+
C
s
2
e
˙
→
2
+
K
s
2
e
→
2
+
C
s
1
(
1
-
ω
1
2
ω
s
)
×
(
e
˙
→
2
-
e
˙
→
1
)
+
K
s
1
(
e
→
2
-
e
→
1
)
=
0
.


## 2.2. 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:(5)
3
2
T
⊗
∂
H
i
∂
θ
1
⊗
∂
S
∂
θ
1
+
1
2
U
⊗
∂
2
S
∂
θ
1
2
+
3
2
T
⊗
∂
H
i
∂
ς
1
⊗
∂
S
∂
ς
1
+
1
2
U
⊗
∂
2
S
∂
ς
1
2
=
Λ
i
P
i
⊗
∂
H
i
∂
θ
1
+
Λ
i
H
i
⊗
∂
P
i
∂
θ
1
+
Λ
i
Ω
P
i
⊗
∂
H
i
∂
θ
1
+
Λ
i
Ω
H
i
⊗
∂
P
i
∂
θ
1
+
2
Λ
i
γ
P
i
⊗
∂
H
i
∂
τ
+
2
Λ
i
γ
H
i
⊗
∂
P
i
∂
τ
.The second film between floating-ring and bushing:(6)
3
2
T
⊗
∂
H
i
i
∂
θ
2
⊗
∂
S
∂
θ
2
+
1
2
U
⊗
∂
2
S
∂
θ
2
2
+
3
2
T
⊗
∂
H
i
i
∂
ς
2
⊗
∂
S
∂
ς
2
+
1
2
U
⊗
∂
2
S
∂
ς
2
2
=
Λ
i
i
P
i
i
⊗
∂
H
i
i
∂
θ
2
+
Λ
i
i
H
i
i
⊗
∂
P
i
i
∂
θ
2
+
2
Λ
i
i
γ
P
i
i
⊗
∂
H
i
i
∂
τ
+
2
Λ
i
i
γ
H
i
i
⊗
∂
P
i
i
∂
τ
,
where
(7)
S
(
k
)
=
P
2
=
P
⊗
P
=
∑
l
=
0
k
P
i
,
j
(
k
-
l
)
P
i
,
j
(
l
)
,
T
(
k
)
=
H
2
=
H
⊗
H
=
∑
l
=
0
k
H
i
,
j
(
k
-
l
)
H
i
,
j
(
l
)
,
U
(
k
)
=
H
3
=
H
⊗
H
⊗
H
=
∑
l
=
0
k
H
i
,
j
(
k
-
l
)
∑
m
=
0
l
H
i
,
j
(
l
-
m
)
H
i
,
j
(
m
)
.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-order-accurate 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:(8)
3
2
∑
l
=
0
k
T
i
,
j
(
k
-
l
)
∑
m
=
0
l
(
H
i
)
i
,
j
(
l
-
m
)
·
(
S
i
+
1
,
j
(
m
)
-
2
S
i
,
j
(
m
)
+
S
i
-
1
,
j
(
m
)
(
Δ
θ
1
)
2
)
+
1
2
∑
l
=
0
k
U
i
,
j
(
k
-
l
)
(
S
i
+
1
,
j
(
l
)
-
2
S
i
,
j
(
l
)
+
S
i
-
1
,
j
(
l
)
(
Δ
θ
1
)
2
)
+
3
2
∑
l
=
0
k
T
i
,
j
(
k
-
l
)
∑
m
=
0
l
(
H
i
)
i
,
j
(
l
-
m
)
v
v
v
v
v
·
(
S
i
,
j
+
1
(
m
)
-
2
S
i
,
j
(
m
)
+
S
i
,
j
-
1
(
m
)
(
Δ
ς
1
)
2
)
+
1
2
∑
l
=
0
k
U
i
,
j
(
k
-
l
)
(
S
i
,
j
+
1
(
l
)
-
2
S
i
,
j
(
l
)
+
S
i
,
j
-
1
(
l
)
(
Δ
ς
1
)
2
)
=
Λ
i
∑
l
=
0
k
(
P
i
)
i
,
j
(
k
-
l
)
(
(
H
i
)
i
+
1
,
j
(
l
)
-
(
H
i
)
i
-
1
,
j
(
l
)
2
Δ
θ
1
)
+
Λ
i
∑
l
=
0
k
(
H
i
)
i
,
j
(
k
-
l
)
(
(
P
i
)
i
+
1
,
j
(
l
)
-
(
P
i
)
i
-
1
,
j
(
l
)
2
Δ
θ
1
)
+
Λ
i
Ω
∑
l
=
0
k
(
P
i
)
i
,
j
(
k
-
l
)
×
(
(
H
i
)
i
+
1
,
j
(
l
)
-
(
H
i
)
i
-
1
,
j
(
l
)
2
Δ
θ
1
)
+
Λ
i
Ω
∑
l
=
0
k
(
H
i
)
i
,
j
(
k
-
l
)
(
(
P
i
)
i
+
1
,
j
(
l
)
-
(
P
i
)
i
-
1
,
j
(
l
)
2
Δ
θ
1
)
+
2
Λ
i
γ
∑
l
=
0
k
[
(
k
+
1
H
~
)
(
P
i
)
i
,
j
(
k
-
l
)
(
H
i
)
i
,
j
(
l
+
1
)
]
+
2
Λ
i
γ
∑
l
=
0
k
[
(
k
+
1
H
~
)
(
H
i
)
i
,
j
(
k
-
l
)
(
P
i
)
i
,
j
(
l
+
1
)
]
.The second film between floating-ring and bushing:(9)
3
2
∑
l
=
0
k
T
i
,
j
(
k
-
l
)
∑
m
=
0
l
(
H
i
i
)
i
,
j
(
l
-
m
)
·
(
S
i
+
1
,
j
(
m
)
-
2
S
i
,
j
(
m
)
+
S
i
-
1
,
j
(
m
)
(
Δ
θ
2
)
2
)
+
3
2
∑
l
=
0
k
U
i
,
j
(
k
-
l
)
(
S
i
+
1
,
j
(
l
)
-
2
S
i
,
j
(
l
)
+
S
i
-
1
,
j
(
l
)
(
Δ
θ
2
)
2
)
+
3
2
∑
l
=
0
k
T
i
,
j
(
k
-
l
)
∑
m
=
0
l
(
H
i
i
)
i
,
j
(
l
-
m
)
·
(
S
i
,
j
+
1
(
m
)
-
2
S
i
,
j
(
m
)
+
S
i
,
j
-
1
(
m
)
(
Δ
ς
2
)
2
)
+
1
2
∑
l
=
0
k
U
i
,
j
(
k
-
l
)
(
S
i
,
j
+
1
(
l
)
-
2
S
i
,
j
(
l
)
+
S
i
,
j
-
1
(
l
)
(
Δ
ς
2
)
2
)
=
Λ
i
i
∑
l
=
0
k
(
P
i
i
)
i
,
j
(
k
-
l
)
(
(
H
i
i
)
i
+
1
,
j
(
l
)
-
(
H
i
i
)
i
-
1
,
j
(
l
)
2
Δ
θ
2
)
+
Λ
i
i
∑
l
=
0
k
(
H
i
i
)
i
,
j
(
k
-
l
)
(
(
P
i
i
)
i
+
1
,
j
(
l
)
-
(
P
i
i
)
i
-
1
,
j
(
l
)
2
Δ
θ
2
)
+
2
Λ
i
i
γ
∑
l
=
0
k
[
(
k
+
1
H
~
)
(
P
i
i
)
i
,
j
(
k
-
l
)
(
H
i
i
)
i
,
j
(
l
+
1
)
]
+
2
Λ
i
i
γ
∑
l
=
0
k
[
(
k
+
1
H
~
)
(
H
i
i
)
i
,
j
(
k
-
l
)
(
P
i
i
)
i
,
j
(
l
+
1
)
]
,whereH
~ is time step value.From (8) and (9), 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 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.


## 3. 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 Table1. 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.Table 1
Comparison of rotor center orbits for steady floating-ring calculated by FDM and hybrid method, respectively.
Displacement
Conditions
X (nT)
Y (nT)
H
~
=
0.001
H
~
=
0.01
H
~
=
0.001
H
~
=
0.01
FDM
m
1
=
4.3 kg
ω
=
6
×
10
6 rpm
−0.1081349035
−0.1181808003
0.0517288907
0.0514134509
Hybrid method
−0.1082889076
−0.1082858013
0.0517174812
0.0517172143
FDM
m
1
=
7.0 kg
ω
=
6
×
10
6 rpm
−0.2582283150
−0.2501341276
0.0658257815
0.0685321157
Hybrid method
−025823243783
−0.2582657430
0.0658701321
0.0658701121In Tables2 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
~.Table 2
Comparison of Poincaré maps of rotor center for steady floating-ring with different values ofm
1 and H
~ by hybrid method.
m
1
=
2.2 kg
m
1
=
8.29 kg
τ
X (nT)
Y (nT)
τ
X (nT)
Y (nT)
π
/300
−0.0457107001
0.0221118339
π/300
−0.0563709651
0.0976387409
π/600
−0.0457891250
0.0221197511
π/600
−0.0563017129
0.0976312750Table 3
Comparison of Poincaré maps of rotor center with different values ofΛ and H
~ by hybrid method.
Λ
=
5.6
Λ
=
9.31
τ
X (nT)
Y (nT)
τ
X (nT)
Y (nT)
π/300
0.0948996672
−0.0373542936
π/300
−0.1656444141
0.0127485637
π/600
0.0948479922
−0.0373550067
π/600
−0.1656440338
0.0127006040


### 3.2. 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 × 106 rpm and the rotor mass m
1 was varied in the range 1.0
≤
m
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 (m
1
=
2.2 kg) but become irregular at a rotor mass of m
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 m
1
=
8.29 and 9.1 kg, the orbits transfer to periodic characteristic. Then increasing rotor mass causes irregular and nonperiodic motion especially when m
1 is greater than 10.08 kg.Trajectories of rotor center atm
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 (Figure 2(a) (2.1–2.12)); phase trajectories of rotor center (Figure 2(b) (2.1–2.12)); power spectra of rotor center in horizontal direction (Figure 2(c) (2.1–2.12)) and vertical direction (Figure 2(d) (2.1–2.12)) (at ω
=
6
×
10
6 rpm).Figures2(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 m
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 m
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 m
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 asm
1 is greater than 11.68 kg.Figures3(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 m
1. Meanwhile, Figures 4(a)–4(j) present the Poincaré maps of the rotor center trajectories at m
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, m
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 m
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
≤
m
1
<
8.18 kg . However, at m
1
=
8.18 kg, the quasiperiodic motion is replaced by T-periodic motion in the horizontal and vertical directions. This motion is maintained over the rotor mass interval 8.18
≤
m
1
<
8.19 kg. Specifically, the rotor center response transits through a cycle of T-Quasi-T periodic behavior as the rotor mass is increased in the range 1.0
≤
m
1
<
10.08 kg (see Figures 4(c)–4(e)). Then, the behavior of rotor center reverts from T periodic to chaotic motion for rotor mass values in the range 10.08
≤
m
1
<
10.7 kg (see Figure 4(f)). This chaotic motion becomes subharmonic motion of multi-T over the rotor mass range 10.7
≤
m
1
<
10.75 kg. The rotor center response transits through a cycle of chaos-multi-T periodic-chaos behavior as the rotor mass is increased from 10.08 to 12.0 kg (see Figures 4(f)–4(l)). Figure 5, corresponding to rotor mass, show that 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).Bifurcation diagrams: (a)X
(
nT
) and (b) Y
(
nT
) versus rotor mass m
1 at ω
=
6
×
10
6 rpm.
(a)
(b)Poincaré maps of rotor center trajectory at (a)m
1
=
2.2, (b) 8.17, (c) 8.22, (d) 8.29, (e) 9.1, (f) 10.08, (g) 10.7, (h) 10.75, (i) 11.23, (j) 11.27, (k) 11.66, and (l) 11.68 kg.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)Maximum Lyapunov exponents of system at different values of rotor mass at (a)m
1
=
10.08 and (b) 11.27 kg.
(a)
(b)From the discussions above, it is evident that the rotor center behavior is significantly dependent on the rotor mass. Table4 summarizes the motions performed by the rotor center for rotor mass values in the interval 1.0
≤
m
1
≤
12.0 kg.Table 4
Behavior of rotor center at different rotor masses in interval1.0
≤
m
1
≤
12.0 kg.
Rotor massm
r (kg)
[1.0, 8.17)
[8.17, 8.18)
[8.18, 8.19)
[8.19, 8.21)
[8.21, 8.22)
Rotor behavior
T
Quasi
T
Quasi
T
Rotor massm
1 (kg)
[8.22, 8.23)
[8.23, 8.26)
[8.26, 8.27)
[8.27, 8.28)
[8.28, 8.29)
Rotor behavior
Quasi
T
Quasi
T
Quasi
Rotor massm
1 (kg)
[8.29, 10.08)
[10.08, 10.7)
[10.7, 10.75)
[10.75, 11.23)
[11.23, 11.27)
Rotor behavior
T
Chaos
Multi-T (10T)
Chaos
Multi-T (8T)
Rotor massm
1 (kg)
[11.27, 11.66)
[11.66, 11.68)
[11.68, 12.0]
Rotor behavior
Chaos
Multi-T (14T)
Chaos


### 3.3. Dynamic Analysis with Increasing Bearing Number

In the second situation, the rotor mass is specified asm
r
=
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.Figure 6
Trajectories of rotor center atΛ = 3.2, 5.6, 9.05, 9.18, 9.31, and 11.2 (Figure 6(a) (6.1–6.8)), phase trajectories of rotor center (Figure 6(b) (6.1–6.8)), and power spectra of rotor center in horizontal direction (Figure 6(c) (6.1–6.8)) and vertical direction (Figure 6(d) (6.1–6.8)) (at m
1
=
2.3 kg).From Figures6(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.Bifurcation diagrams: (a)X
(
nT
) and (b) Y
(
nT
) versus bearing number Λ over interval 1.0
≤
Λ
≤
12.0 at m
1
=
2.3 kg.
(a)
(b)Poincaré maps of rotor center trajectory at (a) 3.2, (b) 5.6, (c) 9.05, (d) 9.18, (e) 9.31, and (f) 11.2.
(a)
(b)
(c)
(d)
(e)
(f)Maximum Lyapunov exponents of system at different values of bearing number at (a)Λ
=
9.05 and (b) 9.18.
(a)
(b)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.Table 5
Behavior of rotor center at different bearing numbers in interval1.0
≤
Λ
≤
12.0.
Bearing numberΛ
[1.0, 9.05)
[9.05, 9.31)
[9.31, 12.0]
Rotor behavior
T
Quasi
T


## 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 Table1. 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.Table 1
Comparison of rotor center orbits for steady floating-ring calculated by FDM and hybrid method, respectively.
Displacement
Conditions
X (nT)
Y (nT)
H
~
=
0.001
H
~
=
0.01
H
~
=
0.001
H
~
=
0.01
FDM
m
1
=
4.3 kg
ω
=
6
×
10
6 rpm
−0.1081349035
−0.1181808003
0.0517288907
0.0514134509
Hybrid method
−0.1082889076
−0.1082858013
0.0517174812
0.0517172143
FDM
m
1
=
7.0 kg
ω
=
6
×
10
6 rpm
−0.2582283150
−0.2501341276
0.0658257815
0.0685321157
Hybrid method
−025823243783
−0.2582657430
0.0658701321
0.0658701121In Tables2 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
~.Table 2
Comparison of Poincaré maps of rotor center for steady floating-ring with different values ofm
1 and H
~ by hybrid method.
m
1
=
2.2 kg
m
1
=
8.29 kg
τ
X (nT)
Y (nT)
τ
X (nT)
Y (nT)
π
/300
−0.0457107001
0.0221118339
π/300
−0.0563709651
0.0976387409
π/600
−0.0457891250
0.0221197511
π/600
−0.0563017129
0.0976312750Table 3
Comparison of Poincaré maps of rotor center with different values ofΛ and H
~ by hybrid method.
Λ
=
5.6
Λ
=
9.31
τ
X (nT)
Y (nT)
τ
X (nT)
Y (nT)
π/300
0.0948996672
−0.0373542936
π/300
−0.1656444141
0.0127485637
π/600
0.0948479922
−0.0373550067
π/600
−0.1656440338
0.0127006040


## 3.2. 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 × 106 rpm and the rotor mass m
1 was varied in the range 1.0
≤
m
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 (m
1
=
2.2 kg) but become irregular at a rotor mass of m
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 m
1
=
8.29 and 9.1 kg, the orbits transfer to periodic characteristic. Then increasing rotor mass causes irregular and nonperiodic motion especially when m
1 is greater than 10.08 kg.Trajectories of rotor center atm
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 (Figure 2(a) (2.1–2.12)); phase trajectories of rotor center (Figure 2(b) (2.1–2.12)); power spectra of rotor center in horizontal direction (Figure 2(c) (2.1–2.12)) and vertical direction (Figure 2(d) (2.1–2.12)) (at ω
=
6
×
10
6 rpm).Figures2(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 m
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 m
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 m
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 asm
1 is greater than 11.68 kg.Figures3(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 m
1. Meanwhile, Figures 4(a)–4(j) present the Poincaré maps of the rotor center trajectories at m
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, m
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 m
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
≤
m
1
<
8.18 kg . However, at m
1
=
8.18 kg, the quasiperiodic motion is replaced by T-periodic motion in the horizontal and vertical directions. This motion is maintained over the rotor mass interval 8.18
≤
m
1
<
8.19 kg. Specifically, the rotor center response transits through a cycle of T-Quasi-T periodic behavior as the rotor mass is increased in the range 1.0
≤
m
1
<
10.08 kg (see Figures 4(c)–4(e)). Then, the behavior of rotor center reverts from T periodic to chaotic motion for rotor mass values in the range 10.08
≤
m
1
<
10.7 kg (see Figure 4(f)). This chaotic motion becomes subharmonic motion of multi-T over the rotor mass range 10.7
≤
m
1
<
10.75 kg. The rotor center response transits through a cycle of chaos-multi-T periodic-chaos behavior as the rotor mass is increased from 10.08 to 12.0 kg (see Figures 4(f)–4(l)). Figure 5, corresponding to rotor mass, show that 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).Bifurcation diagrams: (a)X
(
nT
) and (b) Y
(
nT
) versus rotor mass m
1 at ω
=
6
×
10
6 rpm.
(a)
(b)Poincaré maps of rotor center trajectory at (a)m
1
=
2.2, (b) 8.17, (c) 8.22, (d) 8.29, (e) 9.1, (f) 10.08, (g) 10.7, (h) 10.75, (i) 11.23, (j) 11.27, (k) 11.66, and (l) 11.68 kg.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)Maximum Lyapunov exponents of system at different values of rotor mass at (a)m
1
=
10.08 and (b) 11.27 kg.
(a)
(b)From the discussions above, it is evident that the rotor center behavior is significantly dependent on the rotor mass. Table4 summarizes the motions performed by the rotor center for rotor mass values in the interval 1.0
≤
m
1
≤
12.0 kg.Table 4
Behavior of rotor center at different rotor masses in interval1.0
≤
m
1
≤
12.0 kg.
Rotor massm
r (kg)
[1.0, 8.17)
[8.17, 8.18)
[8.18, 8.19)
[8.19, 8.21)
[8.21, 8.22)
Rotor behavior
T
Quasi
T
Quasi
T
Rotor massm
1 (kg)
[8.22, 8.23)
[8.23, 8.26)
[8.26, 8.27)
[8.27, 8.28)
[8.28, 8.29)
Rotor behavior
Quasi
T
Quasi
T
Quasi
Rotor massm
1 (kg)
[8.29, 10.08)
[10.08, 10.7)
[10.7, 10.75)
[10.75, 11.23)
[11.23, 11.27)
Rotor behavior
T
Chaos
Multi-T (10T)
Chaos
Multi-T (8T)
Rotor massm
1 (kg)
[11.27, 11.66)
[11.66, 11.68)
[11.68, 12.0]
Rotor behavior
Chaos
Multi-T (14T)
Chaos


## 3.3. Dynamic Analysis with Increasing Bearing Number

In the second situation, the rotor mass is specified asm
r
=
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.Figure 6
Trajectories of rotor center atΛ = 3.2, 5.6, 9.05, 9.18, 9.31, and 11.2 (Figure 6(a) (6.1–6.8)), phase trajectories of rotor center (Figure 6(b) (6.1–6.8)), and power spectra of rotor center in horizontal direction (Figure 6(c) (6.1–6.8)) and vertical direction (Figure 6(d) (6.1–6.8)) (at m
1
=
2.3 kg).From Figures6(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.Bifurcation diagrams: (a)X
(
nT
) and (b) Y
(
nT
) versus bearing number Λ over interval 1.0
≤
Λ
≤
12.0 at m
1
=
2.3 kg.
(a)
(b)Poincaré maps of rotor center trajectory at (a) 3.2, (b) 5.6, (c) 9.05, (d) 9.18, (e) 9.31, and (f) 11.2.
(a)
(b)
(c)
(d)
(e)
(f)Maximum Lyapunov exponents of system at different values of bearing number at (a)Λ
=
9.05 and (b) 9.18.
(a)
(b)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.Table 5
Behavior of rotor center at different bearing numbers in interval1.0
≤
Λ
≤
12.0.
Bearing numberΛ
[1.0, 9.05)
[9.05, 9.31)
[9.31, 12.0]
Rotor behavior
T
Quasi
T


## 4. 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 Table1 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.The results of this study provide an understanding of the nonlinear dynamic behavior of EPDAF bearing systems characterized by different rotor massesm
1 and bearing numbers Λ. Specifically, the results have shown that system exists as two kinds of cycle motion including T-Quasi-T periodic and Chaos-Multi-periodic-Chaos over the ranges of rotor mass 1.0
≤
m
1
≤
10.08 and 10.08
≤
m
1
<
12.0 kg, respectively. Regarding the influence of the bearing number on the dynamic response of the bearing system, the results have shown that, at Λ
=
3.2, the Poincaré map contains one discrete point, indicating the presence of T-periodic motion. However, when the bearing number is increased to Λ
=
9.05, the system is transferred to quasiperiodic motion. Then, quasiperiodic motion appears in the bearing number range 9.05
≤
Λ
<
9.31.



---

*Source: 798374-2014-04-28.xml*