# Chaotic Behaviour Investigation of a Front Opposed-Hemispherical Spiral-Grooved Air Bearing System

Authors: Cheng-Chi Wang
Journal: Shock and Vibration
 (2017)

---

## Abstract

In recent years, spiral-grooved air bearing systems have attracted much attention and are especially useful in precision instruments and machines with spindles that rotate at high speed. Load support can be multidirectional and this type of bearing can also be very rigid. Studies show that some of the design problems encountered are dynamic and include critical speed, nonlinearity, gas film pressure, unbalanced rotors, and even poor design, all of which can result in the generation of chaotic aperiodic motion and instability under certain conditions. Such irregular motion on a large scale can cause severe damage to a machine or instrument. Therefore, understanding the conditions under which aperiodic behaviour and vibration arise is crucial for prevention. In this study, numerical analysis, including the Finite Difference and Differential Transformation Methods, is used to study these effects in detail in a front opposed-hemispherical spiral-grooved air bearing system. It was found that different rotor masses and bearing number could cause undesirable behaviour including periodic, subperiodic, quasi-periodic, and chaotic motion. The results obtained in this study can be used as a basis for future bearing system design and the prevention of instability.


---

## Body

## 1. Literature Review

The mathematical theory of gas lubrication was first introduced by Reynolds [1] in 1886 and he derived a partial differential equation related to pressure, density, relative motion, and speed, the well-known Reynolds equation. This established the foundation of fluid lubrication theory. In 1958, Whipple [2] was the first to study the properties of an air bearing with herringbone grooves. In 1959, Whitley and Williams [3] studied the model proposed by Whipple and also published a great number of methods for the design of the herringbone groove. In 1963, Vohr and Pan [4] studied air bearings with spiral grooves and derived a differential equation for the gas film pressure based on an assumption of the number of spiral grooves close to infinity and this became the “Narrow Groove Theory.” Our survey of the literature revealed that the dynamic characteristics, particularly the gas film pressure, gas film rigidity, and the damping effect of the gas film can have a significant impact on the entire bearing supporting system. In 1994, Bonneau and Absi [5] used the finite element method to study characteristics related to the herringbone groove, including stability analysis, pressure analysis, and the impact of several important parameters on the system. Recently, Hirayama et al. [6] discussed the optimization of groove dimensions in herringbone grooved journal bearings for the design of precision spindles with improved run-out characteristics. Schiffmann and Favrat [7] studied the comparison of the influences on the real gas and the perfect gas for herringbone grooved gas journal bearing system and suggested bearing designers to further consider the real gas effects. In 2011, Chen et al. [8] analyzed the groove-ridge discontinuity effect for herringbone grooved gas journal bearings and obtained the pressure distribution of the fluid film on the herringbone grooves. In 2013, Schiffmann [9] optimized the groove geometry of the gas bearing system and set up the design guidelines for enhanced herringbone grooved gas journal bearings. For the transient analysis for the gas bearing system, Hassini and Arghir [10, 11] proposed a novel simplified method for simulating large nonlinear displacements and stability analysis in gas-lubricated bearings. However, none of the references provided an analysis on the dynamic characteristics of a bearing. Furthermore, no series of complete and effective analytic methods for the motion behaviour of spiral groove air bearings and rotors could be found.In 2013, Grigor’ev and Smirnov [12] used finite element method for calculating the characteristics of gas-lubricated spiral-grooved bearings based on integrating the Reynolds equation and the thrust bearing is also considered as a specific example. The results on the dependence of the carrying capacity on the compressibility number and the number of grooves are presented and compared with the asymptotic Narrow Groove Theory. With regard to studies of the dynamic behaviour of rotors, Bently [13] discovered in 1974 that a rotary oil film bearing system could exhibit two- or even three-stage subharmonic vibration and Child et al. [14, 15] used an analytic method to prove the existence of subharmonic vibration in a rotor-bearing system. However, these investigations were all directed towards oil film bearings, and there have been very few studies of rotor behaviour in an air bearing system. From 2007 to 2011, Wang [16–18] used the Finite Difference Method to investigate gas film pressure in a herringbone air bearing system, and he also used orbit graphs, spectrograms, and bifurcation diagrams to analyze the dynamic behaviour of the centres of the rotor and the journal. He learned that, under different operational conditions, the centres of the rotor and the journal demonstrated periodic, quasi-periodic, and subharmonic motion, and he also discovered that the system could generate nonlinear chaotic motion. These results proved useful for this present study. In 2015, Guangwei et al. [19, 20] studied the spiral-grooved opposed-hemisphere gas bearing system with a rigid rotor and focused particular attention on its whirl motion. They used finite element method combined with the Finite Difference Method to solve the time-dependent Reynolds equation and analyze the complicated dynamic behaviour of the rotor-bearing system by phase portraits, power spectra, Poincare maps, and bifurcation diagrams obtained from the numerical procedure. The results revealed that the conical whirl instability appears earlier than the cylindrical whirl instability with increasing rotational speed for the rotor-bearing system with no unbalance mass.In addition, a spiral-grooved air bearing system is typically nonlinear. Conventionally, the steady-state condition of a deterministic nonlinear dynamic system can be regarded as being in one of three main states: equilibrium, periodic motion, and quasi-periodic motion. These may collectively be described as the “attractor” because, after a transient-state response, a system in a steady state can be attracted towards another state. However, recent research shows that a nonlinear dynamic system can exhibit not only regular behaviour but also random or chaotic motion. In a situation where a tiny uncertainty exists in the initial state of motion, the associated dynamic behaviour cannot easily be predicted. The geometric shape formed in the phase space may also differ significantly from that of a normal qualitative system. When chaotic behaviour arises it is very difficult to predict the outcome which may be totally unexpected and even include damage to the system. There is a definite need for a means to prevent the occurrence of such phenomena and the objective of this study is the analysis of the characteristics of a spiral-grooved air bearing and an evaluation of its dynamic behaviour in response to different operational criteria. The bifurcation characteristics associated with nonlinear behaviour generated by the rotor in the system are also discussed. Based on such analysis and study, determinations on whether the chaotic phenomena occur in the system and accurate prediction on the dynamic orbits of the system can be achieved.


## 2. Theoretical Analysis and Model Construction

### 2.1. System Governing Equation

Front opposed-hemispherical spiral-grooved air bearing (FOSAB) systems are used extensively in precision instruments, such as gyroscopes, and in other aerospace applications and are suitable for multidirectional load bearing. They can handle heavy loads making them extremely suitable for many practical engineering applications. In this study two different numerical analytic methods have been used to study a FOSAB system and the performance and dynamic characteristics of the bearing were investigated in some detail. The diagrams in Figures1 and 2 show the structure of a FOSAB in which the centre O
x
y
z and the centre of the sphere coincide. The z-axis is vertical and downwards and the corresponding spherical coordinate is (φ, θ, ξ).Figure 1
The structure and coordinate system of hemispherical spiral-grooved air bearing.
(a)
(b)
(c)Figure 2
Types of opposed-hemispherical spiral-grooved air bearing: (a) front opposed type and (b) rear opposed type.As shown in the Figures1 and 2, the spherical radius of the bearing is R
o, the average clearance is C
r, and the rotational speed is ω. Parameters of the spiral groove include angle of groove α, depth of groove δ, width of groove b
1, and separation between grooves b
2. The compressible Reynolds equation of the bearing of such type uses the spherical coordinate system in principle, which is modified as follows:(1)
1
sin
⁡
θ
∂
∂
θ
F
θ
+
∂
∂
φ
F
φ
=
2
Λ
∂
∂
τ
P
^
h
^wherein(2)
F
θ
=
P
^
η
1
sin
⁡
θ
∂
P
^
∂
θ
+
η
2
∂
P
^
∂
φ
+
2
λ
c
Λ
η
4
Ω
2
-
Ω
1
sin
2
⁡
θ
,
(3)
F
φ
=
P
^
η
2
sin
⁡
θ
∂
P
^
∂
θ
+
η
3
∂
P
^
∂
φ
+
2
λ
s
Λ
η
4
Ω
2
-
Ω
1
sin
2
⁡
θ
-
Λ
h
^
Ω
1
+
Ω
2
sin
2
⁡
θ
,
(4)
P
^
=
P
P
a
,
δ
^
1
=
δ
h
2
,
B
^
1
=
b
1
b
,
b
^
1
=
b
2
b
1
,
Λ
=
6
μ
ω
P
a
R
o
C
r
2
,
(5)
Ω
1
=
ω
1
ω
=
-
1
g
r
o
o
v
e
d
s
u
r
f
a
c
e
r
o
t
a
t
i
o
n
,
Ω
2
=
ω
2
ω
=
1
s
m
o
o
t
h
s
u
r
f
a
c
e
r
o
t
a
t
i
o
n
,
(6)
h
^
o
=
1
-
ε
z
cos
⁡
θ
,
h
^
=
h
^
o
-
ε
r
cos
⁡
φ
sin
⁡
θ
,
(7)
h
^
1
=
h
1
C
r
=
h
^
+
B
^
1
δ
^
1
,
h
^
2
=
h
2
C
r
=
h
^
+
δ
^
1
,
(8)
β
1
=
h
^
2
3
-
h
^
3
B
^
1
1
-
B
^
1
,
β
2
=
1
-
B
^
1
h
^
2
3
+
B
^
1
h
^
3
,
(9)
η
1
=
β
1
sin
2
⁡
α
+
h
^
2
3
h
^
3
β
2
,
η
2
=
-
β
1
β
2
,
(10)
η
3
=
β
1
cos
2
⁡
α
+
h
^
2
3
h
^
3
β
2
,
η
4
=
h
^
2
3
-
h
^
3
β
2
,
(11)
λ
c
=
1
-
B
^
1
B
^
1
δ
^
1
sin
⁡
α
cos
⁡
α
2
,
λ
s
=
1
-
B
^
1
B
^
1
δ
^
1
sin
2
⁡
α
2
.P
^ is the dimensionless gas film pressure, h
^ is the dimensionless gas film thickness, ρ is the gas film density, and μ is the gas film viscosity. To solve the bearing gas pressure distribution the traditional Finite Difference Method and a hybrid method [21, 22] were used. The hybrid method is a combination of the Differential Transformation Method (DTM) and the Finite Difference Method (FDM). For the nonlinear Reynolds equation time domain, the Differential Transformation Method was first used to discretize the time into intervals and the central difference scheme of the Finite Difference Method was then applied to the coordinates of the locations. The Differential Transformation Method was used to solve the nonlinear Reynolds equation, and the basic principle is as follows.The differential transformation of the function off
(
t
) is defined as(12)
f
-
k
=
T
k
k
!
∂
k
f
t
∂
t
k
t
=
0
.Whereinf
-
k refers to the transformation function in the transformation domain, f
t refers to the original function in the time domain, T refers to the time separation, and k refers to the conversion parameter. The reverse equation of f
-
k is(13)
f
t
=
∑
k
=
0
∞
f
-
k
t
T
k
.Calibration is needed to ensure precision during the calculation; therefore, when the second time interval is used, (14) and (15) are used to perform calibration.(14)
f
i
-
1
0
=
f
i
0
+
∑
k
=
1
∞
f
-
i
k
-
T
i
-
1
T
i
k
.IfT
o
=
T
1
=
T
2
=
⋯ then the above equation becomes(15)
f
i
-
1
0
=
f
i
0
+
∑
k
=
1
∞
f
-
i
k
-
1
k
.If the calculated values on each side of the equals’ signs in (14) and (15) vary greatly, precision is insufficient. Therefore, the time interval T needs to be reduced and recalculation is necessary.The basic arithmetic operation for the differential transformation is as shown in Table1.Table 1
Basic arithmetic formula for the Differential Transformation.
Before transformation
After transformation
a
f
t
±
b
g
t
a
f
-
k
±
b
g
-
k
f
t
g
t
∑
l
=
0
k
f
-
k
-
l
g
-
l
=
∑
l
=
0
k
f
-
l
g
-
k
-
l
f
t
g
t
f
-
k
-
∑
l
=
0
k
-
1
f
-
l
/
g
-
l
g
-
k
-
l
g
-
0
=
f
-
k
-
∑
l
=
0
k
-
1
f
-
k
-
l
/
g
-
k
-
l
g
-
l
g
-
0
cos
⁡
a
t
+
b, a
≠
0
a
T
k
k
!
cos
⁡
π
k
2
+
b
d
n
f
t
d
t
n
k
+
n
!
k
!
T
n
f
-
k
+
n
f
m
t
∑
l
=
0
k
f
-
k
-
l
∑
n
=
0
l
f
-
m
-
2
n
f
-
l
-
nThe differential transformation of the governing equation and the boundary condition are obtained and the pressure function theorem domain is divided inton number of subintervals; the pressure solution p
f for each subinternal is expressed as an inverse transformation equation. If the theorem domains of the subintervals are T
0
,
T
1
,
T
2
,
T
3
,
…, then the expression for the function in the first interval, based on (13), will be(16)
P
o
t
=
∑
k
=
0
∞
P
-
o
k
t
T
o
k
,
0
≤
t
≤
T
o
.If the initial value ofP
o
0 is known, then, from the above equation, P
o
0
=
P
-
o
0 can be obtained. From P
-
o
0, other discrete values of P
-
o
k, k
=
1,2
,
… in this interval can be obtained. In the first subinterval, the destination value P
o
t
=
T
o of the function of P
o
t refers to the initial value P
1
(
t
=
0
) of the function P
1
t in the second subinterval. The function P
1
t in the second subinterval is(17)
P
1
t
=
∑
k
=
0
∞
P
-
1
k
t
T
1
k
,
0
≤
t
≤
T
1
.FromP
-
1
0, other P
-
1
k, k
=
1,2
,
… can also be obtained. In such order, the value of P
-
i
k can be obtained. Then, the inverse differential transformation can be further used to obtain the solution at each time domain, and the equation for the series solution is as follows: (18)
P
i
t
=
∑
k
=
0
m
P
-
i
k
t
T
i
k
,
0
≤
t
≤
T
i
.


### 2.2. Numerical Analysis

The traditional Finite Difference Method is used first on the Reynolds equation and to apply the central difference scheme to the spatial coordinates. For the time domain, the implicit backwards difference method is used. Equation (1) is transformed into a differential equation, and the Finite Difference Method of Successive Over Relation (SOR) is used to reduce the five unknown numbers to three, following which iteration is used to obtain the correct pressure value at each time interval. Iteration at each time interval is performed in the θ direction only, and the three points of I, I
-
1, I
+
1 are unknown numbers. The two points of J
+
1 and J
-
1 that are originally unknown are replaced by the previous iteration values. Consequently, all the equations can be formed into a triple-diagonal matrix, which can be then solved by Gaussian Elimination.After obtaining the pressure distribution between the bearing and the rotor, the bearing force of the gas film can be obtained by integration. The numerical calculations of the hybrid method can be done using (12)–(18), as shown in the following:(19)
η
1
∂
P
^
∂
θ
⊗
∂
P
^
∂
θ
+
η
2
csc
⁡
θ
⊗
∂
P
^
∂
θ
⊗
∂
P
^
∂
φ
+
2
Λ
λ
c
η
4
Ω
2
-
Ω
1
sin
⁡
θ
⊗
∂
P
^
∂
θ
+
η
1
cot
⁡
θ
⊗
∂
P
^
∂
θ
⊗
P
^
+
η
1
P
^
⊗
∂
2
P
^
∂
θ
2
+
η
2
csc
⁡
θ
⊗
P
^
⊗
∂
2
P
^
∂
φ
∂
θ
+
4
Λ
λ
c
η
4
Ω
2
-
Ω
1
cos
⁡
θ
⊗
P
^
+
η
2
sin
⁡
θ
⊗
∂
P
^
∂
θ
⊗
∂
P
^
∂
φ
+
η
2
∂
P
^
∂
φ
⊗
∂
P
^
∂
φ
+
2
Λ
λ
s
η
4
Ω
2
-
Ω
1
sin
⁡
θ
⊗
sin
⁡
θ
⊗
∂
P
^
∂
φ
+
η
2
sin
⁡
θ
⊗
P
^
⊗
∂
2
P
^
∂
φ
∂
θ
+
η
3
P
^
⊗
∂
2
P
^
∂
φ
2
-
2
Λ
Ω
2
+
Ω
1
h
^
⊗
sin
⁡
θ
⊗
sin
⁡
θ
⊗
∂
P
^
∂
φ
+
Λ
Ω
2
+
Ω
1
P
^
⊗
sin
⁡
θ
⊗
sin
⁡
θ
⊗
∂
h
^
∂
φ
=
2
Λ
∂
P
^
∂
τ
⊗
h
^
+
2
Λ
∂
h
^
∂
τ
⊗
P
^
.Next, the Differential Transformation Method is applied to the time domain in (19) to produce discrete time intervals, and the central difference scheme of the Finite Difference Method is then applied to the coordinates of the locations. The following can then be derived:(20)
η
1
∑
l
=
0
k
P
^
i
+
1
,
j
k
-
l
-
P
^
i
-
1
,
j
k
-
l
2
Δ
θ
·
P
^
i
+
1
,
j
l
-
P
^
i
-
1
,
j
l
2
Δ
θ
+
η
2
csc
⁡
i
Δ
θ
·
∑
l
=
0
k
P
^
i
+
1
,
j
k
-
l
-
P
^
i
-
1
,
j
k
-
l
2
Δ
θ
·
P
^
i
,
j
+
1
l
-
P
^
i
,
j
-
1
l
2
Δ
φ
+
2
Λ
λ
c
η
4
Ω
2
-
Ω
1
∑
l
=
0
k
·
P
^
i
+
1
,
j
k
-
l
-
P
^
i
-
1
,
j
k
-
l
2
Δ
θ
sin
⁡
l
Δ
θ
+
η
1
cot
⁡
i
Δ
θ
·
∑
l
=
0
k
P
^
i
+
1
,
j
k
-
l
-
P
^
i
-
1
,
j
k
-
l
2
Δ
θ
·
P
^
i
,
j
l
+
η
1
∑
l
=
0
k
P
^
i
+
1
,
j
k
-
l
-
2
·
P
^
i
,
j
k
-
l
+
P
^
i
-
1
,
j
k
-
l
Δ
θ
2
·
P
^
i
,
j
l
+
η
2
csc
⁡
i
Δ
θ
·
∑
l
=
0
k
P
^
i
+
1
,
j
+
1
k
-
l
-
P
^
i
+
1
,
j
-
1
k
-
l
-
P
^
i
-
1
,
j
+
1
k
-
l
+
P
^
i
-
1
,
j
-
1
k
-
l
Δ
θ
Δ
φ
·
P
^
i
,
j
l
+
4
Λ
λ
c
η
4
Ω
2
-
Ω
1
∑
l
=
0
k
P
^
i
,
j
k
-
l
·
sin
⁡
l
Δ
θ
+
η
2
sin
⁡
i
Δ
θ
·
∑
l
=
0
k
P
^
i
+
1
,
j
k
-
l
-
P
^
i
-
1
,
j
k
-
l
2
Δ
θ
·
P
^
i
,
j
+
1
l
-
P
^
i
,
j
-
1
l
2
Δ
φ
+
η
2
∑
l
=
0
k
P
^
i
,
j
+
1
k
-
l
-
P
^
i
,
j
-
1
k
-
l
2
Δ
φ
·
P
^
i
,
j
+
1
l
-
P
^
i
,
j
-
1
l
2
Δ
φ
+
2
Λ
λ
s
η
4
Ω
2
-
Ω
1
sin
⁡
i
Δ
θ
·
∑
l
=
0
k
·
P
^
i
,
j
+
1
k
-
l
-
P
^
i
,
j
-
1
k
-
l
2
Δ
φ
sin
⁡
l
Δ
θ
+
η
2
sin
⁡
i
Δ
θ
·
∑
l
=
0
k
P
^
i
+
1
,
j
+
1
k
-
l
-
P
^
i
+
1
,
j
-
1
k
-
l
-
P
^
i
-
1
,
j
+
1
k
-
l
+
P
^
i
-
1
,
j
-
1
k
-
l
Δ
θ
Δ
φ
·
P
^
i
,
j
l
-
2
Λ
Ω
2
+
Ω
1
sin
⁡
i
Δ
θ
·
∑
l
=
0
k
h
^
i
,
j
k
-
l
∑
m
=
0
l
·
P
^
i
,
j
+
1
l
-
m
-
P
^
i
,
j
-
1
l
-
m
2
Δ
φ
sin
⁡
m
Δ
θ
+
Λ
Ω
2
+
Ω
1
sin
⁡
i
Δ
θ
·
∑
l
=
0
k
P
^
i
,
j
k
-
l
∑
m
=
0
l
·
h
^
i
,
j
+
1
l
-
m
-
h
^
i
,
j
-
1
l
-
m
2
Δ
φ
sin
⁡
m
Δ
θ
+
η
3
∑
l
=
0
k
P
^
i
,
j
+
1
k
-
l
-
2
·
P
^
i
,
j
k
-
l
+
P
^
i
,
j
-
1
k
-
l
Δ
φ
2
·
P
^
i
,
j
l
=
2
Λ
·
∑
l
=
0
k
l
+
1
H
~
P
^
i
,
j
k
-
l
·
h
^
i
,
j
l
+
1
+
2
Λ
·
∑
l
=
0
k
l
+
1
H
~
h
^
i
,
j
k
-
l
·
P
^
i
,
j
l
+
1
.


### 2.3. Dynamic Equations of the Rotor

In this system, the ends of a flexible rotor of massm
r are supported on two identical and symmetrical air bearings secured on a rigid base. Under such ideal conditions the motion of the rotor can be treated as two having two directions of movement. O
2 refers to the geometric centre of the rotor, O
3 to the geometric centre of the journal, ρ
¯ to the eccentricity of the rotor, K
p to the rigidity of the rotating axle, and ω to the rotational speed of the rotor.In a steady state, the equation of motion of the rotor centreO
2 is as follows:(21)
m
r
d
2
x
2
d
t
2
+
K
p
x
2
-
x
3
=
m
r
ρ
¯
ω
2
cos
⁡
ω
t
,
(22)
m
r
d
2
y
2
d
t
2
+
K
p
y
2
-
y
3
=
m
r
ρ
¯
ω
2
sin
⁡
ω
t
.The force equilibrium equation at the journal centreO
3 is as follows:(23)
f
g
x
=
K
p
x
2
-
x
3
2
,
(24)
f
g
y
=
K
p
y
2
-
y
3
2
.Upon dimensionless transformation:(25)
X
2
=
x
2
C
r
;
X
3
=
x
3
C
r
;
Y
2
=
y
2
C
r
;
Y
3
=
y
3
C
r
;
τ
=
ω
t
;
d
d
t
=
ω
d
d
τ
.In addition, the dimensionless parameters are introduced as follows:(26)
F
g
x
=
f
g
x
P
a
R
L
;
F
g
y
=
f
g
y
P
a
R
L
;
M
r
=
m
r
ω
2
K
p
;
ξ
=
m
r
ρ
ω
2
K
p
C
r
.By substituting (25) into (21)~(24) and introducing the dimensionless parameters (26), the following equations can be obtained: (27)
F
g
x
=
K
p
C
r
2
P
a
R
L
X
2
-
X
3
,
(28)
F
g
y
=
K
p
C
r
2
P
a
R
L
Y
2
-
Y
3
,
(29)
M
r
d
2
X
2
d
τ
2
+
X
2
-
X
3
=
ξ
cos
⁡
τ
,
(30)
M
r
d
2
Y
2
d
τ
2
+
Y
2
-
Y
3
=
ξ
sin
⁡
τ
.In this study iteration was used to calculate the motion behaviour of the rotor. First, the acceleration, speed, and displacement were calculated, step by step, with respect to time. The initial state of the system was taken as static, and the initial displacement (X
2
o, Y
2
o) of the rotor refers to the location of the rotor centre when the system was in static equilibrium; such equilibrium position also defines the initial gas film thickness value h
^. In addition, the initial speed was set to zero.The steps of calculation are as follows.Step 1.
At timeτ
=
0, the initial condition is used; after a period of time Δ
τ, (31)~(36) can be used to calculate the new acceleration, speed, and displacement of the rotor:(31)
d
2
X
2
d
τ
2
=
ξ
cos
⁡
τ
-
X
2
-
X
3
M
r
=
A
X
,
(32)
d
2
Y
2
d
τ
2
=
ξ
sin
⁡
τ
-
Y
2
-
Y
3
M
r
=
A
Y
,
(33)
V
X
=
V
X
o
+
A
X
·
Δ
τ
,
(34)
V
Y
=
V
Y
o
+
A
Y
·
Δ
τ
,
(35)
X
2
=
X
2
o
+
V
X
·
Δ
τ
+
1
2
A
X
·
Δ
τ
2
,
(36)
Y
2
=
Y
2
o
+
V
Y
·
Δ
τ
+
1
2
A
Y
·
Δ
τ
2
.Step 2.
The displacement of the rotor centre obtained from Step1 can be used to calculate the new displacement of the journal centre. The rotation of the journal causes changes in the gas film thickness h
^; consequently, after the new displacement value yields a new h
^, this is substituted into (2) to calculate a new pressure value, and the gas film bearing force can be derived by integration.Step 3.
From the displacement and velocity values obtained in Step1 and the pressure distribution from Step 2, a new set of initial conditions can be derived from the gas film force obtained. When the time is increased by Δ
τ, (i.e., Δ
τ
→
2
Δ
τ), the aforementioned new initial conditions are used again to further calculate all the variations in the system. These three steps allow determination of the pressure distribution, the relationship between the gas film thickness and change of position, variation of the gas film bearing force, the dynamic orbits of the rotor, and so on at each time point. In addition, the maximum Lyapunov exponents can be used as a primary determination factor to predict periodic, quasi-periodic, aperiodic, and chaotic motion caused by such factors as the bearing number and mass of the rotor. By using the bearing number and the rotor mass of the air bearing as bifurcation parameters, a full view bifurcation diagram can be constructed. An analysis of stability can be conducted and the relationship between the maximum Lyapunov exponents and each of the previously mentioned factors can be interpreted. The parameters of the stable and the nonstable region of the system can be determined.To further verify accuracy of the numerical analysis, two different numerical methods were used (as previously mentioned) to ensure that the hybrid method proposed in this paper was applicable and to increase the calculation precision.


## 2.1. System Governing Equation

Front opposed-hemispherical spiral-grooved air bearing (FOSAB) systems are used extensively in precision instruments, such as gyroscopes, and in other aerospace applications and are suitable for multidirectional load bearing. They can handle heavy loads making them extremely suitable for many practical engineering applications. In this study two different numerical analytic methods have been used to study a FOSAB system and the performance and dynamic characteristics of the bearing were investigated in some detail. The diagrams in Figures1 and 2 show the structure of a FOSAB in which the centre O
x
y
z and the centre of the sphere coincide. The z-axis is vertical and downwards and the corresponding spherical coordinate is (φ, θ, ξ).Figure 1
The structure and coordinate system of hemispherical spiral-grooved air bearing.
(a)
(b)
(c)Figure 2
Types of opposed-hemispherical spiral-grooved air bearing: (a) front opposed type and (b) rear opposed type.As shown in the Figures1 and 2, the spherical radius of the bearing is R
o, the average clearance is C
r, and the rotational speed is ω. Parameters of the spiral groove include angle of groove α, depth of groove δ, width of groove b
1, and separation between grooves b
2. The compressible Reynolds equation of the bearing of such type uses the spherical coordinate system in principle, which is modified as follows:(1)
1
sin
⁡
θ
∂
∂
θ
F
θ
+
∂
∂
φ
F
φ
=
2
Λ
∂
∂
τ
P
^
h
^wherein(2)
F
θ
=
P
^
η
1
sin
⁡
θ
∂
P
^
∂
θ
+
η
2
∂
P
^
∂
φ
+
2
λ
c
Λ
η
4
Ω
2
-
Ω
1
sin
2
⁡
θ
,
(3)
F
φ
=
P
^
η
2
sin
⁡
θ
∂
P
^
∂
θ
+
η
3
∂
P
^
∂
φ
+
2
λ
s
Λ
η
4
Ω
2
-
Ω
1
sin
2
⁡
θ
-
Λ
h
^
Ω
1
+
Ω
2
sin
2
⁡
θ
,
(4)
P
^
=
P
P
a
,
δ
^
1
=
δ
h
2
,
B
^
1
=
b
1
b
,
b
^
1
=
b
2
b
1
,
Λ
=
6
μ
ω
P
a
R
o
C
r
2
,
(5)
Ω
1
=
ω
1
ω
=
-
1
g
r
o
o
v
e
d
s
u
r
f
a
c
e
r
o
t
a
t
i
o
n
,
Ω
2
=
ω
2
ω
=
1
s
m
o
o
t
h
s
u
r
f
a
c
e
r
o
t
a
t
i
o
n
,
(6)
h
^
o
=
1
-
ε
z
cos
⁡
θ
,
h
^
=
h
^
o
-
ε
r
cos
⁡
φ
sin
⁡
θ
,
(7)
h
^
1
=
h
1
C
r
=
h
^
+
B
^
1
δ
^
1
,
h
^
2
=
h
2
C
r
=
h
^
+
δ
^
1
,
(8)
β
1
=
h
^
2
3
-
h
^
3
B
^
1
1
-
B
^
1
,
β
2
=
1
-
B
^
1
h
^
2
3
+
B
^
1
h
^
3
,
(9)
η
1
=
β
1
sin
2
⁡
α
+
h
^
2
3
h
^
3
β
2
,
η
2
=
-
β
1
β
2
,
(10)
η
3
=
β
1
cos
2
⁡
α
+
h
^
2
3
h
^
3
β
2
,
η
4
=
h
^
2
3
-
h
^
3
β
2
,
(11)
λ
c
=
1
-
B
^
1
B
^
1
δ
^
1
sin
⁡
α
cos
⁡
α
2
,
λ
s
=
1
-
B
^
1
B
^
1
δ
^
1
sin
2
⁡
α
2
.P
^ is the dimensionless gas film pressure, h
^ is the dimensionless gas film thickness, ρ is the gas film density, and μ is the gas film viscosity. To solve the bearing gas pressure distribution the traditional Finite Difference Method and a hybrid method [21, 22] were used. The hybrid method is a combination of the Differential Transformation Method (DTM) and the Finite Difference Method (FDM). For the nonlinear Reynolds equation time domain, the Differential Transformation Method was first used to discretize the time into intervals and the central difference scheme of the Finite Difference Method was then applied to the coordinates of the locations. The Differential Transformation Method was used to solve the nonlinear Reynolds equation, and the basic principle is as follows.The differential transformation of the function off
(
t
) is defined as(12)
f
-
k
=
T
k
k
!
∂
k
f
t
∂
t
k
t
=
0
.Whereinf
-
k refers to the transformation function in the transformation domain, f
t refers to the original function in the time domain, T refers to the time separation, and k refers to the conversion parameter. The reverse equation of f
-
k is(13)
f
t
=
∑
k
=
0
∞
f
-
k
t
T
k
.Calibration is needed to ensure precision during the calculation; therefore, when the second time interval is used, (14) and (15) are used to perform calibration.(14)
f
i
-
1
0
=
f
i
0
+
∑
k
=
1
∞
f
-
i
k
-
T
i
-
1
T
i
k
.IfT
o
=
T
1
=
T
2
=
⋯ then the above equation becomes(15)
f
i
-
1
0
=
f
i
0
+
∑
k
=
1
∞
f
-
i
k
-
1
k
.If the calculated values on each side of the equals’ signs in (14) and (15) vary greatly, precision is insufficient. Therefore, the time interval T needs to be reduced and recalculation is necessary.The basic arithmetic operation for the differential transformation is as shown in Table1.Table 1
Basic arithmetic formula for the Differential Transformation.
Before transformation
After transformation
a
f
t
±
b
g
t
a
f
-
k
±
b
g
-
k
f
t
g
t
∑
l
=
0
k
f
-
k
-
l
g
-
l
=
∑
l
=
0
k
f
-
l
g
-
k
-
l
f
t
g
t
f
-
k
-
∑
l
=
0
k
-
1
f
-
l
/
g
-
l
g
-
k
-
l
g
-
0
=
f
-
k
-
∑
l
=
0
k
-
1
f
-
k
-
l
/
g
-
k
-
l
g
-
l
g
-
0
cos
⁡
a
t
+
b, a
≠
0
a
T
k
k
!
cos
⁡
π
k
2
+
b
d
n
f
t
d
t
n
k
+
n
!
k
!
T
n
f
-
k
+
n
f
m
t
∑
l
=
0
k
f
-
k
-
l
∑
n
=
0
l
f
-
m
-
2
n
f
-
l
-
nThe differential transformation of the governing equation and the boundary condition are obtained and the pressure function theorem domain is divided inton number of subintervals; the pressure solution p
f for each subinternal is expressed as an inverse transformation equation. If the theorem domains of the subintervals are T
0
,
T
1
,
T
2
,
T
3
,
…, then the expression for the function in the first interval, based on (13), will be(16)
P
o
t
=
∑
k
=
0
∞
P
-
o
k
t
T
o
k
,
0
≤
t
≤
T
o
.If the initial value ofP
o
0 is known, then, from the above equation, P
o
0
=
P
-
o
0 can be obtained. From P
-
o
0, other discrete values of P
-
o
k, k
=
1,2
,
… in this interval can be obtained. In the first subinterval, the destination value P
o
t
=
T
o of the function of P
o
t refers to the initial value P
1
(
t
=
0
) of the function P
1
t in the second subinterval. The function P
1
t in the second subinterval is(17)
P
1
t
=
∑
k
=
0
∞
P
-
1
k
t
T
1
k
,
0
≤
t
≤
T
1
.FromP
-
1
0, other P
-
1
k, k
=
1,2
,
… can also be obtained. In such order, the value of P
-
i
k can be obtained. Then, the inverse differential transformation can be further used to obtain the solution at each time domain, and the equation for the series solution is as follows: (18)
P
i
t
=
∑
k
=
0
m
P
-
i
k
t
T
i
k
,
0
≤
t
≤
T
i
.


## 2.2. Numerical Analysis

The traditional Finite Difference Method is used first on the Reynolds equation and to apply the central difference scheme to the spatial coordinates. For the time domain, the implicit backwards difference method is used. Equation (1) is transformed into a differential equation, and the Finite Difference Method of Successive Over Relation (SOR) is used to reduce the five unknown numbers to three, following which iteration is used to obtain the correct pressure value at each time interval. Iteration at each time interval is performed in the θ direction only, and the three points of I, I
-
1, I
+
1 are unknown numbers. The two points of J
+
1 and J
-
1 that are originally unknown are replaced by the previous iteration values. Consequently, all the equations can be formed into a triple-diagonal matrix, which can be then solved by Gaussian Elimination.After obtaining the pressure distribution between the bearing and the rotor, the bearing force of the gas film can be obtained by integration. The numerical calculations of the hybrid method can be done using (12)–(18), as shown in the following:(19)
η
1
∂
P
^
∂
θ
⊗
∂
P
^
∂
θ
+
η
2
csc
⁡
θ
⊗
∂
P
^
∂
θ
⊗
∂
P
^
∂
φ
+
2
Λ
λ
c
η
4
Ω
2
-
Ω
1
sin
⁡
θ
⊗
∂
P
^
∂
θ
+
η
1
cot
⁡
θ
⊗
∂
P
^
∂
θ
⊗
P
^
+
η
1
P
^
⊗
∂
2
P
^
∂
θ
2
+
η
2
csc
⁡
θ
⊗
P
^
⊗
∂
2
P
^
∂
φ
∂
θ
+
4
Λ
λ
c
η
4
Ω
2
-
Ω
1
cos
⁡
θ
⊗
P
^
+
η
2
sin
⁡
θ
⊗
∂
P
^
∂
θ
⊗
∂
P
^
∂
φ
+
η
2
∂
P
^
∂
φ
⊗
∂
P
^
∂
φ
+
2
Λ
λ
s
η
4
Ω
2
-
Ω
1
sin
⁡
θ
⊗
sin
⁡
θ
⊗
∂
P
^
∂
φ
+
η
2
sin
⁡
θ
⊗
P
^
⊗
∂
2
P
^
∂
φ
∂
θ
+
η
3
P
^
⊗
∂
2
P
^
∂
φ
2
-
2
Λ
Ω
2
+
Ω
1
h
^
⊗
sin
⁡
θ
⊗
sin
⁡
θ
⊗
∂
P
^
∂
φ
+
Λ
Ω
2
+
Ω
1
P
^
⊗
sin
⁡
θ
⊗
sin
⁡
θ
⊗
∂
h
^
∂
φ
=
2
Λ
∂
P
^
∂
τ
⊗
h
^
+
2
Λ
∂
h
^
∂
τ
⊗
P
^
.Next, the Differential Transformation Method is applied to the time domain in (19) to produce discrete time intervals, and the central difference scheme of the Finite Difference Method is then applied to the coordinates of the locations. The following can then be derived:(20)
η
1
∑
l
=
0
k
P
^
i
+
1
,
j
k
-
l
-
P
^
i
-
1
,
j
k
-
l
2
Δ
θ
·
P
^
i
+
1
,
j
l
-
P
^
i
-
1
,
j
l
2
Δ
θ
+
η
2
csc
⁡
i
Δ
θ
·
∑
l
=
0
k
P
^
i
+
1
,
j
k
-
l
-
P
^
i
-
1
,
j
k
-
l
2
Δ
θ
·
P
^
i
,
j
+
1
l
-
P
^
i
,
j
-
1
l
2
Δ
φ
+
2
Λ
λ
c
η
4
Ω
2
-
Ω
1
∑
l
=
0
k
·
P
^
i
+
1
,
j
k
-
l
-
P
^
i
-
1
,
j
k
-
l
2
Δ
θ
sin
⁡
l
Δ
θ
+
η
1
cot
⁡
i
Δ
θ
·
∑
l
=
0
k
P
^
i
+
1
,
j
k
-
l
-
P
^
i
-
1
,
j
k
-
l
2
Δ
θ
·
P
^
i
,
j
l
+
η
1
∑
l
=
0
k
P
^
i
+
1
,
j
k
-
l
-
2
·
P
^
i
,
j
k
-
l
+
P
^
i
-
1
,
j
k
-
l
Δ
θ
2
·
P
^
i
,
j
l
+
η
2
csc
⁡
i
Δ
θ
·
∑
l
=
0
k
P
^
i
+
1
,
j
+
1
k
-
l
-
P
^
i
+
1
,
j
-
1
k
-
l
-
P
^
i
-
1
,
j
+
1
k
-
l
+
P
^
i
-
1
,
j
-
1
k
-
l
Δ
θ
Δ
φ
·
P
^
i
,
j
l
+
4
Λ
λ
c
η
4
Ω
2
-
Ω
1
∑
l
=
0
k
P
^
i
,
j
k
-
l
·
sin
⁡
l
Δ
θ
+
η
2
sin
⁡
i
Δ
θ
·
∑
l
=
0
k
P
^
i
+
1
,
j
k
-
l
-
P
^
i
-
1
,
j
k
-
l
2
Δ
θ
·
P
^
i
,
j
+
1
l
-
P
^
i
,
j
-
1
l
2
Δ
φ
+
η
2
∑
l
=
0
k
P
^
i
,
j
+
1
k
-
l
-
P
^
i
,
j
-
1
k
-
l
2
Δ
φ
·
P
^
i
,
j
+
1
l
-
P
^
i
,
j
-
1
l
2
Δ
φ
+
2
Λ
λ
s
η
4
Ω
2
-
Ω
1
sin
⁡
i
Δ
θ
·
∑
l
=
0
k
·
P
^
i
,
j
+
1
k
-
l
-
P
^
i
,
j
-
1
k
-
l
2
Δ
φ
sin
⁡
l
Δ
θ
+
η
2
sin
⁡
i
Δ
θ
·
∑
l
=
0
k
P
^
i
+
1
,
j
+
1
k
-
l
-
P
^
i
+
1
,
j
-
1
k
-
l
-
P
^
i
-
1
,
j
+
1
k
-
l
+
P
^
i
-
1
,
j
-
1
k
-
l
Δ
θ
Δ
φ
·
P
^
i
,
j
l
-
2
Λ
Ω
2
+
Ω
1
sin
⁡
i
Δ
θ
·
∑
l
=
0
k
h
^
i
,
j
k
-
l
∑
m
=
0
l
·
P
^
i
,
j
+
1
l
-
m
-
P
^
i
,
j
-
1
l
-
m
2
Δ
φ
sin
⁡
m
Δ
θ
+
Λ
Ω
2
+
Ω
1
sin
⁡
i
Δ
θ
·
∑
l
=
0
k
P
^
i
,
j
k
-
l
∑
m
=
0
l
·
h
^
i
,
j
+
1
l
-
m
-
h
^
i
,
j
-
1
l
-
m
2
Δ
φ
sin
⁡
m
Δ
θ
+
η
3
∑
l
=
0
k
P
^
i
,
j
+
1
k
-
l
-
2
·
P
^
i
,
j
k
-
l
+
P
^
i
,
j
-
1
k
-
l
Δ
φ
2
·
P
^
i
,
j
l
=
2
Λ
·
∑
l
=
0
k
l
+
1
H
~
P
^
i
,
j
k
-
l
·
h
^
i
,
j
l
+
1
+
2
Λ
·
∑
l
=
0
k
l
+
1
H
~
h
^
i
,
j
k
-
l
·
P
^
i
,
j
l
+
1
.


## 2.3. Dynamic Equations of the Rotor

In this system, the ends of a flexible rotor of massm
r are supported on two identical and symmetrical air bearings secured on a rigid base. Under such ideal conditions the motion of the rotor can be treated as two having two directions of movement. O
2 refers to the geometric centre of the rotor, O
3 to the geometric centre of the journal, ρ
¯ to the eccentricity of the rotor, K
p to the rigidity of the rotating axle, and ω to the rotational speed of the rotor.In a steady state, the equation of motion of the rotor centreO
2 is as follows:(21)
m
r
d
2
x
2
d
t
2
+
K
p
x
2
-
x
3
=
m
r
ρ
¯
ω
2
cos
⁡
ω
t
,
(22)
m
r
d
2
y
2
d
t
2
+
K
p
y
2
-
y
3
=
m
r
ρ
¯
ω
2
sin
⁡
ω
t
.The force equilibrium equation at the journal centreO
3 is as follows:(23)
f
g
x
=
K
p
x
2
-
x
3
2
,
(24)
f
g
y
=
K
p
y
2
-
y
3
2
.Upon dimensionless transformation:(25)
X
2
=
x
2
C
r
;
X
3
=
x
3
C
r
;
Y
2
=
y
2
C
r
;
Y
3
=
y
3
C
r
;
τ
=
ω
t
;
d
d
t
=
ω
d
d
τ
.In addition, the dimensionless parameters are introduced as follows:(26)
F
g
x
=
f
g
x
P
a
R
L
;
F
g
y
=
f
g
y
P
a
R
L
;
M
r
=
m
r
ω
2
K
p
;
ξ
=
m
r
ρ
ω
2
K
p
C
r
.By substituting (25) into (21)~(24) and introducing the dimensionless parameters (26), the following equations can be obtained: (27)
F
g
x
=
K
p
C
r
2
P
a
R
L
X
2
-
X
3
,
(28)
F
g
y
=
K
p
C
r
2
P
a
R
L
Y
2
-
Y
3
,
(29)
M
r
d
2
X
2
d
τ
2
+
X
2
-
X
3
=
ξ
cos
⁡
τ
,
(30)
M
r
d
2
Y
2
d
τ
2
+
Y
2
-
Y
3
=
ξ
sin
⁡
τ
.In this study iteration was used to calculate the motion behaviour of the rotor. First, the acceleration, speed, and displacement were calculated, step by step, with respect to time. The initial state of the system was taken as static, and the initial displacement (X
2
o, Y
2
o) of the rotor refers to the location of the rotor centre when the system was in static equilibrium; such equilibrium position also defines the initial gas film thickness value h
^. In addition, the initial speed was set to zero.The steps of calculation are as follows.Step 1.
At timeτ
=
0, the initial condition is used; after a period of time Δ
τ, (31)~(36) can be used to calculate the new acceleration, speed, and displacement of the rotor:(31)
d
2
X
2
d
τ
2
=
ξ
cos
⁡
τ
-
X
2
-
X
3
M
r
=
A
X
,
(32)
d
2
Y
2
d
τ
2
=
ξ
sin
⁡
τ
-
Y
2
-
Y
3
M
r
=
A
Y
,
(33)
V
X
=
V
X
o
+
A
X
·
Δ
τ
,
(34)
V
Y
=
V
Y
o
+
A
Y
·
Δ
τ
,
(35)
X
2
=
X
2
o
+
V
X
·
Δ
τ
+
1
2
A
X
·
Δ
τ
2
,
(36)
Y
2
=
Y
2
o
+
V
Y
·
Δ
τ
+
1
2
A
Y
·
Δ
τ
2
.Step 2.
The displacement of the rotor centre obtained from Step1 can be used to calculate the new displacement of the journal centre. The rotation of the journal causes changes in the gas film thickness h
^; consequently, after the new displacement value yields a new h
^, this is substituted into (2) to calculate a new pressure value, and the gas film bearing force can be derived by integration.Step 3.
From the displacement and velocity values obtained in Step1 and the pressure distribution from Step 2, a new set of initial conditions can be derived from the gas film force obtained. When the time is increased by Δ
τ, (i.e., Δ
τ
→
2
Δ
τ), the aforementioned new initial conditions are used again to further calculate all the variations in the system. These three steps allow determination of the pressure distribution, the relationship between the gas film thickness and change of position, variation of the gas film bearing force, the dynamic orbits of the rotor, and so on at each time point. In addition, the maximum Lyapunov exponents can be used as a primary determination factor to predict periodic, quasi-periodic, aperiodic, and chaotic motion caused by such factors as the bearing number and mass of the rotor. By using the bearing number and the rotor mass of the air bearing as bifurcation parameters, a full view bifurcation diagram can be constructed. An analysis of stability can be conducted and the relationship between the maximum Lyapunov exponents and each of the previously mentioned factors can be interpreted. The parameters of the stable and the nonstable region of the system can be determined.To further verify accuracy of the numerical analysis, two different numerical methods were used (as previously mentioned) to ensure that the hybrid method proposed in this paper was applicable and to increase the calculation precision.


## 3. Results and Discussion

### 3.1. Comparison between Results Generated by Different Numerical Methods

The results obtained by the hybrid method that combines both the Differential Transformation Method (DTM) and the Finite Difference Method (FDM) agree with the results of the traditional Finite Difference Method and Perturbation Method. This was confirmed by use of the two different methods as previously mentioned. The values of the Finite Difference Method and Perturbation Method obtained under certain parameters would cause unstable phenomena; in addition, the result also shows that the hybrid method (DTM and FDM) yields greater accuracy. Table2 shows a comparison of the values of the trace of the rotor centre (X
2, Y
2), wherein H
~ refers to the time step for calculation. With regard to the analysis on the stability of the values, the impact of different time step on the values is also completed, showing that the stability of the values obtained from the hybrid method is indeed greater than those obtained from the traditional Finite Difference Method. Furthermore, precision of the displacement calculation can be accurate to five decimal places by hybrid method.Table 2
Comparison of numerical calculation results by Poincaré map of numerical displacement of rotor centre.
Conditions
Poincaré map
X
2(n
T)
Y
2(n
T)
H
~
=
0.001
H
~
=
0.01
H
~
=
0.001
H
~
=
0.01
FDM
m
r = 5.1 kg
ω = 6 × 106 rpm
0.1037130279
0.1038378524
−0.2539313115
−0.2538991587
DTM and FDM
0.1038771861
0.1038756587
−0.2535875862
−0.2535826268
Perturbation Method
0.1032651536
0.1070747058
−0.2558305761
−0.2523962686
FDM
m
r =6.2 kg
ω = 6 × 106 rpm
0.1333187451
0.1330440082
−0.2452767681
−0.2454073653
DTM and FDM
0.1334811874
0.1334689894
−0.2457199581
−0.2457613815
Perturbation Method
0.1332782833
0.1331421019
−0.2458672612
−0.2456736825The results show that the hybrid method presented here demonstrates excellent convergence and applicability. The method can be used to effectively calculate the orbits of the system with different rotor mass and increased bearing number (rotational speed). Examination of the test results in Table3 shows that the time step selected does not need to be highly precise to obtain sufficiently accurate numerical results and this further reduces the bifurcation characteristics of subsequent calculations. Therefore, π
/
300 is used as the time step for dynamic analysis calculation.Table 3
Comparison of rotor centre displacement under different time steps by hybrid method.
m
r = 6.2 kg
m
r = 11.1 kg
τ
X
2
(
n
T
)
Y
2
(
n
T
)
τ
X
2
(
n
T
)
Y
2
(
n
T
)
π/300
0.1334893205
−0.2457499192
π/300
−0.0195491566
−0.4737438386
π/600
0.1334152991
−0.2457851961
π/600
−0.0195681361
−0.4737158324
Λ
=
2.0
Λ
=
5.0
τ
X
2
(
n
T
)
Y
2
(
n
T
)
τ
X
2
(
n
T
)
Y
2
(
n
T
)
π/300
−0.0923588407
−0.0567265833
π/300
0.5145289218
−0.0817031729
π/600
−0.0923990874
−0.0567583091
π/600
0.5145152184
−0.0817172361


### 3.2. Analysis of the Dynamic Behaviour of a FOSAB System: Using Rotor Mass as a Bifurcation Parameter

In this section, analysis of the impact of rotor mass on the air bearing system is conducted; here the bearing number of the FOSAB system is set toΛ
=
2.2, and the rotor mass m
r is used as the bifurcation parameter.


#### 3.2.1. Dynamic Orbits Analysis

An examination of Figures3(a1, b1), 3(a2, b2), …, and 3(a7, b7) show that the orbits of the centres of the rotor (X
2, Y
2) and the journal (X
3, Y
3) will follow regular paths when the mass is small (m
r = 5.1 kg). However, when the mass is increased to m
r = 7.06 kg, aperiodic motion occurs. This shows that, with a rotor mass of less than 7.06 kg, the system will remain relatively stable. When the rotor mass reaches m
r = 8.06 kg, aperiodic motion is reduced but comes back when m
r = 8.21 kg. However, the nonperiodicity stabilizes when m
r = 8.52 kg. This stable behaviour does not last long, and the rotor centre becomes unstable again when its mass reaches 9.30 kg. Finally, at m
r = 10.76 kg, the system becomes relatively stable.Figure 3
Dynamic orbits ((a1)–(a7)) of the rotor centre of the front opposed-hemispherical spiral-grooved air bearing system when the rotor mass ism
r = 5.1, 7.06, 8.06, 8.21, 8.52, 9.30, and 10.76 kg, respectively; dynamic orbits of the journal centre ((b1)–(b7)); spectrum response diagrams for the rotor centre in the horizontal direction ((c1)–(c7)) and the vertical direction ((d1)–(d7)), where the system bearing number is Λ
=
2.2.


#### 3.2.2. Spectrum Analysis

Figures3(c1, d1), 3(c2, d2), …, and 3(c7, d7) show the dynamic response of the rotor centre in the horizontal and vertical directions of rotation. The results show that when the rotor mass is m
r = 5.1 and 10.76 kg, motion of the rotor centre is periodic; when m
r = 7.06 and 9.30 kg, the spectrum response graphs (Figures 3(c2, d2), 3(c6, d6)) show that the motion of the rotor centre becomes chaotic. When the rotor mass is 8.06 and 8.52 kg, the system generates 2T subharmonic motion. However, motion of the system is quasi-periodic when the rotor mass is 8.21 kg.


#### 3.2.3. Bifurcation Analysis

The effect of different rotor mass on the system was studied using rotor massm
r as an analytic parameter, as shown in Figure 4. The mass was set within a range of 1.0 to 12.0 kg according to actual operating conditions. Figure 5 shows that when m
r < 7.06 kg, the rotor centre undergoes T periodic motion in both the horizontal and vertical directions, and this is illustrated by the Poincaré map shown in Figure 5(a). Nevertheless, when the mass was increased to m
r = 7.06 kg, the T periodic motion of the rotor centre was replaced by chaotic motion, as shown in Figure 5(b). When the mass reached m
r = 8.06 kg, the system returned to 2T subperiodic motion (see Figure 5(c)). Although this motion continued over the range 8.06 ≤ m
r < 8.21 kg it did not last to the end. When m
r increased to 8.21 kg, quasi-periodic motion appeared (see Figure 5(d)) and continued over the range 8.21 ≤ m
r < 8.52 kg. This is illustrated by the closed curve formed by multiple discrete points shown on the Poincaré map in Figure 5(d). As the mass increased to m
r = 8.52 kg, the quasi-periodic motion of the original system bifurcated into 2T subharmonic motion, as shown in Figure 5(e) where the map shows two discrete points and this type of motion is present throughout the range of 8.52 ≤ m
r < 9.30 kg. Next, as the mass increased further to m
r = 9.30 kg, the subharmonic motion of the original system became chaotic motion, as shown in Figure 5(f). Here the map generates two attraction subregions, and this is present throughout the range of 9.30 ≤ m
r < 10.76 kg. Finally, when the rotor mass is increased to m
r = 10.76 kg, the chaotic motion of the system has become a periodic motion that is relatively stable, as shown in Figure 5(g).Figure 4
Bifurcation diagrams of (a)X
2(n
T) and (b) Y
2(n
t) for the rotor centre of the front opposed-hemispherical spiral-grooved air bearing system with respect to different rotor masses, where the system bearing number is Λ
=
2.2.
(a)
(b)Figure 5
Poincaré maps of the rotor centre with respect to different rotor masses: (a)m
r = 5.1, (b) 7.06, (c) 8.06, (d) 8.21, (e) 8.52, (f) 9.30, and (g) 10.76 kg.
(a)
(b)
(c)
(d)
(e)
(f)
(g)The maximum Lyapunov exponents have been very useful for the analysis and verification of the chaotic behaviour of the system. From Figures6(a), 6(c), 6(d), 6(e), and 6(g), it can be seen that when m
r = 5.1, 8.06, 8.21, 8.52 and 10.76 kg, the maximum Lyapunov exponents are close to zero, and the system is nonchaotic. However, when m
r = 7.06 and 9.30, as shown in Figures 6(b) and 6(f), the exponents are greater than zero and the system state of motion is indeed chaotic. A full view of the maximum Lyapunov exponents of the system is shown in Figure 7, and the result indicates that the areas of chaotic behaviour (exponent greater than zero) occur at two intervals, 7.06 ≤ m
r < 8.06 kg and 9.30 ≤ m
r < 10.76 kg. This is in clear agreement with the results. However, these results show that the motion of the rotor can be a relatively complicated model. A summary of the different kinds of system motion, with rotors of different mass, is given in Table 4.Table 4
Dynamic behaviour of rotor centre with different rotor mass.
Rotor mass
[1.0, 7.06)
[7.06, 8.06)
[8.06, 8.21)
[8.21, 8.52)
[8.52, 9.3)
[9.3, 10.76)
[10.76, 12.0]
Behavior
T
Chaotic
2T
Quasi
2T
Chaotic
TFigure 6
Maximum Lyapunov exponents of the rotor centre with respect to different rotor masses: (a)m
r = 5.1, (b) 7.06, (c) 8.06, (d) 8.21, (e) 8.52, (f) 9.30, and (g) 10.76 kg.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Figure 7
Full maximum Lyapunov exponents of the rotor centre with different rotor mass.


### 3.3. Analysis on Dynamic Behaviour of FOSAB System: Using the Bearing Number as the Bifurcation Parameter

For an air bearing system, the rotational speed referred to that bearing number of the rotor-bearing has a direct effect on pressure distribution inside the bearing, which in turn has a direct effect on the relevant performance and stability of the entire system. Accordingly, in this section, the bearing numberΛ is used as the bifurcation parameter, and the rotor mass is set to m
r = 2.3 kg to analyze and study the dynamic behaviour of the system with respect to bearing number.


#### 3.3.1. Dynamic Orbits Analysis

An examination of Figures8(a1, b1), 8(a2, b2), …, and 8(a6, b6) shows that when the bearing number is small, Λ = 2.3, 3.12, the behaviour of the journal and the rotor will be relatively regular. However, when Λ is increased to 3.84, behaviour of the journal and the rotor centre will become irregular and asymmetrical. When the bearing number is further increased to 3.86, the dynamic orbits return to relatively regular periodic motion. However, irregular motion reappears when Λ
=
4.55 and only becomes regular again when the bearing number is increased to 4.89.Figure 8
Dynamic orbits ((a1)–(a6)) of the rotor centre when the bearing numberΛ = 2.3, 3.12, 3.84, 3.86, 4.55, and 4.89, respectively. Dynamic orbits of the journal centre ((b1)–(b6)). Spectrum response diagrams for the rotor centre in the horizontal direction ((c1)–(c6)) and the vertical direction of rotation ((d1)–(d6)), with a rotor mass m
r = 2.3 kg.


#### 3.3.2. Spectrum Analysis

Figures8(c1, d1), 8(c2, d2), …, and 8(c6, d6) show the dynamic responses of the rotor centre in the horizontal and vertical directions of rotation with different bearing numbers. The study shows that at bearing numbers of Λ = 2.3, 3.86, and 4.89 the rotor centre demonstrates T periodic motions. When Λ = 3.84, the spectrum response diagram shows aperiodic motion of the rotor centre, but when the bearing number is Λ = 3.12, the motion becomes 2T subharmonic motion. However, when Λ = 4.55, the rotor centre motion becomes aperiodic.


#### 3.3.3. Bifurcation Analysis

As shown in Figure9, the bearing number Λ can be used as a primary analytic parameter for studying the behaviour of an air bearing system and has been set within a range of 1.0 to 5.5 according to actual operational conditions. When the bearing number is small, such as Λ
=
2.3, the rotor centre exhibits T periodic motion in both the horizontal and rotational directions and this can be verified by the appearance of a single dot on the Poincaré map as shown in Figure 10(a). When the bearing number is increased to Λ
=
3.12, the system will bifurcate and generate 2T periodic motion, clearly shown by the presence of two discrete dots on the map as shown in Figure 10(b). When the bearing number increases to 3.84
≤
Λ
<
3.86, the system motion becomes chaotic, clearly shown by the irregular discrete points generated on the map in Figure 10(c). When Λ
=
3.86, the chaotic motion becomes T periodic motion, as shown in Figure 10(d), and the same occurs in the interval of 3.86
≤
Λ
<
4.55; see Figure 9. In addition, for the interval of 4.55
≤
Λ
<
4.89, the system generates quasi-periodic motion, as shown in Figure 10(e). When Λ
=
4.89, the system assumes T-periodic motion that is relatively stable and this continues up to Λ
=
5.5; see Figure 10(f).Figure 9
Bifurcation diagrams: (a)X
2(n
T) and (b) Y
2(n
t) for the rotor centre with respect to different bearing number Λ, with a rotor mass m
r = 2.3 kg.
(a)
(b)Figure 10
Poincaré maps of the rotor centre with respect to different bearing numbersΛ = (a) 2.3, (b) 3.12, (c) 3.84, (d) 3.86, (e) 4.55, and (f) 4.89.
(a)
(b)
(c)
(d)
(e)
(f)With regard to the question about whether a change of bearing number will cause the generation of chaotic behaviour, in this paper, similarly, the maximum Lyapunov exponents are also used for determination and further verification. From Figures11(a), 11(b), 11(d), 11(e), and 11(f), it can be seen that when Λ = 2.3, 3.12, 3.86, 4.55, and 4.89, the exponents obtained can all be smaller than or equal to zero and the system motion is nonchaotic. In view of the complexity of the relationship between the motion of the rotor centre and bearing number, all motions and behaviours generated with respect to the interval of 1.0
≤
Λ
≤
5.5 are summarized in Table 5 and a full view of the maximum Lyapunov exponents for the system is shown in Figure 12. The result indicates that chaotic regions (exponents greater than zero) only occur when 3.84
≤
Λ
<
3.86, which agrees with the results shown in Table 5.Table 5
Dynamic behaviour of the rotor centre with different bearing number (1.0
≤
Λ
≤
5.5).
Λ
[1.0, 3.12)
[3.12, 3.84)
[3.84, 3.86)
[3.86, 4.55)
[4.55, 4.89)
[4.89, 5.5]
Behavior
T
2T
Chaos
2T
Quasi
TFigure 11
Maximum Lyapunov exponents of the rotor centre with respect to different bearing numberΛ = (a) 2.3, (b) 3.12, (c) 3.84, (d) 3.86, (e) 4.55, and (f) 4.89.
(a)
(b)
(c)
(d)
(e)
(f)Figure 12
Full view of maximum Lyapunov exponents of the rotor centre with respect to different rotor-bearing numbers.


## 3.1. Comparison between Results Generated by Different Numerical Methods

The results obtained by the hybrid method that combines both the Differential Transformation Method (DTM) and the Finite Difference Method (FDM) agree with the results of the traditional Finite Difference Method and Perturbation Method. This was confirmed by use of the two different methods as previously mentioned. The values of the Finite Difference Method and Perturbation Method obtained under certain parameters would cause unstable phenomena; in addition, the result also shows that the hybrid method (DTM and FDM) yields greater accuracy. Table2 shows a comparison of the values of the trace of the rotor centre (X
2, Y
2), wherein H
~ refers to the time step for calculation. With regard to the analysis on the stability of the values, the impact of different time step on the values is also completed, showing that the stability of the values obtained from the hybrid method is indeed greater than those obtained from the traditional Finite Difference Method. Furthermore, precision of the displacement calculation can be accurate to five decimal places by hybrid method.Table 2
Comparison of numerical calculation results by Poincaré map of numerical displacement of rotor centre.
Conditions
Poincaré map
X
2(n
T)
Y
2(n
T)
H
~
=
0.001
H
~
=
0.01
H
~
=
0.001
H
~
=
0.01
FDM
m
r = 5.1 kg
ω = 6 × 106 rpm
0.1037130279
0.1038378524
−0.2539313115
−0.2538991587
DTM and FDM
0.1038771861
0.1038756587
−0.2535875862
−0.2535826268
Perturbation Method
0.1032651536
0.1070747058
−0.2558305761
−0.2523962686
FDM
m
r =6.2 kg
ω = 6 × 106 rpm
0.1333187451
0.1330440082
−0.2452767681
−0.2454073653
DTM and FDM
0.1334811874
0.1334689894
−0.2457199581
−0.2457613815
Perturbation Method
0.1332782833
0.1331421019
−0.2458672612
−0.2456736825The results show that the hybrid method presented here demonstrates excellent convergence and applicability. The method can be used to effectively calculate the orbits of the system with different rotor mass and increased bearing number (rotational speed). Examination of the test results in Table3 shows that the time step selected does not need to be highly precise to obtain sufficiently accurate numerical results and this further reduces the bifurcation characteristics of subsequent calculations. Therefore, π
/
300 is used as the time step for dynamic analysis calculation.Table 3
Comparison of rotor centre displacement under different time steps by hybrid method.
m
r = 6.2 kg
m
r = 11.1 kg
τ
X
2
(
n
T
)
Y
2
(
n
T
)
τ
X
2
(
n
T
)
Y
2
(
n
T
)
π/300
0.1334893205
−0.2457499192
π/300
−0.0195491566
−0.4737438386
π/600
0.1334152991
−0.2457851961
π/600
−0.0195681361
−0.4737158324
Λ
=
2.0
Λ
=
5.0
τ
X
2
(
n
T
)
Y
2
(
n
T
)
τ
X
2
(
n
T
)
Y
2
(
n
T
)
π/300
−0.0923588407
−0.0567265833
π/300
0.5145289218
−0.0817031729
π/600
−0.0923990874
−0.0567583091
π/600
0.5145152184
−0.0817172361


## 3.2. Analysis of the Dynamic Behaviour of a FOSAB System: Using Rotor Mass as a Bifurcation Parameter

In this section, analysis of the impact of rotor mass on the air bearing system is conducted; here the bearing number of the FOSAB system is set toΛ
=
2.2, and the rotor mass m
r is used as the bifurcation parameter.


### 3.2.1. Dynamic Orbits Analysis

An examination of Figures3(a1, b1), 3(a2, b2), …, and 3(a7, b7) show that the orbits of the centres of the rotor (X
2, Y
2) and the journal (X
3, Y
3) will follow regular paths when the mass is small (m
r = 5.1 kg). However, when the mass is increased to m
r = 7.06 kg, aperiodic motion occurs. This shows that, with a rotor mass of less than 7.06 kg, the system will remain relatively stable. When the rotor mass reaches m
r = 8.06 kg, aperiodic motion is reduced but comes back when m
r = 8.21 kg. However, the nonperiodicity stabilizes when m
r = 8.52 kg. This stable behaviour does not last long, and the rotor centre becomes unstable again when its mass reaches 9.30 kg. Finally, at m
r = 10.76 kg, the system becomes relatively stable.Figure 3
Dynamic orbits ((a1)–(a7)) of the rotor centre of the front opposed-hemispherical spiral-grooved air bearing system when the rotor mass ism
r = 5.1, 7.06, 8.06, 8.21, 8.52, 9.30, and 10.76 kg, respectively; dynamic orbits of the journal centre ((b1)–(b7)); spectrum response diagrams for the rotor centre in the horizontal direction ((c1)–(c7)) and the vertical direction ((d1)–(d7)), where the system bearing number is Λ
=
2.2.


### 3.2.2. Spectrum Analysis

Figures3(c1, d1), 3(c2, d2), …, and 3(c7, d7) show the dynamic response of the rotor centre in the horizontal and vertical directions of rotation. The results show that when the rotor mass is m
r = 5.1 and 10.76 kg, motion of the rotor centre is periodic; when m
r = 7.06 and 9.30 kg, the spectrum response graphs (Figures 3(c2, d2), 3(c6, d6)) show that the motion of the rotor centre becomes chaotic. When the rotor mass is 8.06 and 8.52 kg, the system generates 2T subharmonic motion. However, motion of the system is quasi-periodic when the rotor mass is 8.21 kg.


### 3.2.3. Bifurcation Analysis

The effect of different rotor mass on the system was studied using rotor massm
r as an analytic parameter, as shown in Figure 4. The mass was set within a range of 1.0 to 12.0 kg according to actual operating conditions. Figure 5 shows that when m
r < 7.06 kg, the rotor centre undergoes T periodic motion in both the horizontal and vertical directions, and this is illustrated by the Poincaré map shown in Figure 5(a). Nevertheless, when the mass was increased to m
r = 7.06 kg, the T periodic motion of the rotor centre was replaced by chaotic motion, as shown in Figure 5(b). When the mass reached m
r = 8.06 kg, the system returned to 2T subperiodic motion (see Figure 5(c)). Although this motion continued over the range 8.06 ≤ m
r < 8.21 kg it did not last to the end. When m
r increased to 8.21 kg, quasi-periodic motion appeared (see Figure 5(d)) and continued over the range 8.21 ≤ m
r < 8.52 kg. This is illustrated by the closed curve formed by multiple discrete points shown on the Poincaré map in Figure 5(d). As the mass increased to m
r = 8.52 kg, the quasi-periodic motion of the original system bifurcated into 2T subharmonic motion, as shown in Figure 5(e) where the map shows two discrete points and this type of motion is present throughout the range of 8.52 ≤ m
r < 9.30 kg. Next, as the mass increased further to m
r = 9.30 kg, the subharmonic motion of the original system became chaotic motion, as shown in Figure 5(f). Here the map generates two attraction subregions, and this is present throughout the range of 9.30 ≤ m
r < 10.76 kg. Finally, when the rotor mass is increased to m
r = 10.76 kg, the chaotic motion of the system has become a periodic motion that is relatively stable, as shown in Figure 5(g).Figure 4
Bifurcation diagrams of (a)X
2(n
T) and (b) Y
2(n
t) for the rotor centre of the front opposed-hemispherical spiral-grooved air bearing system with respect to different rotor masses, where the system bearing number is Λ
=
2.2.
(a)
(b)Figure 5
Poincaré maps of the rotor centre with respect to different rotor masses: (a)m
r = 5.1, (b) 7.06, (c) 8.06, (d) 8.21, (e) 8.52, (f) 9.30, and (g) 10.76 kg.
(a)
(b)
(c)
(d)
(e)
(f)
(g)The maximum Lyapunov exponents have been very useful for the analysis and verification of the chaotic behaviour of the system. From Figures6(a), 6(c), 6(d), 6(e), and 6(g), it can be seen that when m
r = 5.1, 8.06, 8.21, 8.52 and 10.76 kg, the maximum Lyapunov exponents are close to zero, and the system is nonchaotic. However, when m
r = 7.06 and 9.30, as shown in Figures 6(b) and 6(f), the exponents are greater than zero and the system state of motion is indeed chaotic. A full view of the maximum Lyapunov exponents of the system is shown in Figure 7, and the result indicates that the areas of chaotic behaviour (exponent greater than zero) occur at two intervals, 7.06 ≤ m
r < 8.06 kg and 9.30 ≤ m
r < 10.76 kg. This is in clear agreement with the results. However, these results show that the motion of the rotor can be a relatively complicated model. A summary of the different kinds of system motion, with rotors of different mass, is given in Table 4.Table 4
Dynamic behaviour of rotor centre with different rotor mass.
Rotor mass
[1.0, 7.06)
[7.06, 8.06)
[8.06, 8.21)
[8.21, 8.52)
[8.52, 9.3)
[9.3, 10.76)
[10.76, 12.0]
Behavior
T
Chaotic
2T
Quasi
2T
Chaotic
TFigure 6
Maximum Lyapunov exponents of the rotor centre with respect to different rotor masses: (a)m
r = 5.1, (b) 7.06, (c) 8.06, (d) 8.21, (e) 8.52, (f) 9.30, and (g) 10.76 kg.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Figure 7
Full maximum Lyapunov exponents of the rotor centre with different rotor mass.


## 3.2.1. Dynamic Orbits Analysis

An examination of Figures3(a1, b1), 3(a2, b2), …, and 3(a7, b7) show that the orbits of the centres of the rotor (X
2, Y
2) and the journal (X
3, Y
3) will follow regular paths when the mass is small (m
r = 5.1 kg). However, when the mass is increased to m
r = 7.06 kg, aperiodic motion occurs. This shows that, with a rotor mass of less than 7.06 kg, the system will remain relatively stable. When the rotor mass reaches m
r = 8.06 kg, aperiodic motion is reduced but comes back when m
r = 8.21 kg. However, the nonperiodicity stabilizes when m
r = 8.52 kg. This stable behaviour does not last long, and the rotor centre becomes unstable again when its mass reaches 9.30 kg. Finally, at m
r = 10.76 kg, the system becomes relatively stable.Figure 3
Dynamic orbits ((a1)–(a7)) of the rotor centre of the front opposed-hemispherical spiral-grooved air bearing system when the rotor mass ism
r = 5.1, 7.06, 8.06, 8.21, 8.52, 9.30, and 10.76 kg, respectively; dynamic orbits of the journal centre ((b1)–(b7)); spectrum response diagrams for the rotor centre in the horizontal direction ((c1)–(c7)) and the vertical direction ((d1)–(d7)), where the system bearing number is Λ
=
2.2.


## 3.2.2. Spectrum Analysis

Figures3(c1, d1), 3(c2, d2), …, and 3(c7, d7) show the dynamic response of the rotor centre in the horizontal and vertical directions of rotation. The results show that when the rotor mass is m
r = 5.1 and 10.76 kg, motion of the rotor centre is periodic; when m
r = 7.06 and 9.30 kg, the spectrum response graphs (Figures 3(c2, d2), 3(c6, d6)) show that the motion of the rotor centre becomes chaotic. When the rotor mass is 8.06 and 8.52 kg, the system generates 2T subharmonic motion. However, motion of the system is quasi-periodic when the rotor mass is 8.21 kg.


## 3.2.3. Bifurcation Analysis

The effect of different rotor mass on the system was studied using rotor massm
r as an analytic parameter, as shown in Figure 4. The mass was set within a range of 1.0 to 12.0 kg according to actual operating conditions. Figure 5 shows that when m
r < 7.06 kg, the rotor centre undergoes T periodic motion in both the horizontal and vertical directions, and this is illustrated by the Poincaré map shown in Figure 5(a). Nevertheless, when the mass was increased to m
r = 7.06 kg, the T periodic motion of the rotor centre was replaced by chaotic motion, as shown in Figure 5(b). When the mass reached m
r = 8.06 kg, the system returned to 2T subperiodic motion (see Figure 5(c)). Although this motion continued over the range 8.06 ≤ m
r < 8.21 kg it did not last to the end. When m
r increased to 8.21 kg, quasi-periodic motion appeared (see Figure 5(d)) and continued over the range 8.21 ≤ m
r < 8.52 kg. This is illustrated by the closed curve formed by multiple discrete points shown on the Poincaré map in Figure 5(d). As the mass increased to m
r = 8.52 kg, the quasi-periodic motion of the original system bifurcated into 2T subharmonic motion, as shown in Figure 5(e) where the map shows two discrete points and this type of motion is present throughout the range of 8.52 ≤ m
r < 9.30 kg. Next, as the mass increased further to m
r = 9.30 kg, the subharmonic motion of the original system became chaotic motion, as shown in Figure 5(f). Here the map generates two attraction subregions, and this is present throughout the range of 9.30 ≤ m
r < 10.76 kg. Finally, when the rotor mass is increased to m
r = 10.76 kg, the chaotic motion of the system has become a periodic motion that is relatively stable, as shown in Figure 5(g).Figure 4
Bifurcation diagrams of (a)X
2(n
T) and (b) Y
2(n
t) for the rotor centre of the front opposed-hemispherical spiral-grooved air bearing system with respect to different rotor masses, where the system bearing number is Λ
=
2.2.
(a)
(b)Figure 5
Poincaré maps of the rotor centre with respect to different rotor masses: (a)m
r = 5.1, (b) 7.06, (c) 8.06, (d) 8.21, (e) 8.52, (f) 9.30, and (g) 10.76 kg.
(a)
(b)
(c)
(d)
(e)
(f)
(g)The maximum Lyapunov exponents have been very useful for the analysis and verification of the chaotic behaviour of the system. From Figures6(a), 6(c), 6(d), 6(e), and 6(g), it can be seen that when m
r = 5.1, 8.06, 8.21, 8.52 and 10.76 kg, the maximum Lyapunov exponents are close to zero, and the system is nonchaotic. However, when m
r = 7.06 and 9.30, as shown in Figures 6(b) and 6(f), the exponents are greater than zero and the system state of motion is indeed chaotic. A full view of the maximum Lyapunov exponents of the system is shown in Figure 7, and the result indicates that the areas of chaotic behaviour (exponent greater than zero) occur at two intervals, 7.06 ≤ m
r < 8.06 kg and 9.30 ≤ m
r < 10.76 kg. This is in clear agreement with the results. However, these results show that the motion of the rotor can be a relatively complicated model. A summary of the different kinds of system motion, with rotors of different mass, is given in Table 4.Table 4
Dynamic behaviour of rotor centre with different rotor mass.
Rotor mass
[1.0, 7.06)
[7.06, 8.06)
[8.06, 8.21)
[8.21, 8.52)
[8.52, 9.3)
[9.3, 10.76)
[10.76, 12.0]
Behavior
T
Chaotic
2T
Quasi
2T
Chaotic
TFigure 6
Maximum Lyapunov exponents of the rotor centre with respect to different rotor masses: (a)m
r = 5.1, (b) 7.06, (c) 8.06, (d) 8.21, (e) 8.52, (f) 9.30, and (g) 10.76 kg.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Figure 7
Full maximum Lyapunov exponents of the rotor centre with different rotor mass.


## 3.3. Analysis on Dynamic Behaviour of FOSAB System: Using the Bearing Number as the Bifurcation Parameter

For an air bearing system, the rotational speed referred to that bearing number of the rotor-bearing has a direct effect on pressure distribution inside the bearing, which in turn has a direct effect on the relevant performance and stability of the entire system. Accordingly, in this section, the bearing numberΛ is used as the bifurcation parameter, and the rotor mass is set to m
r = 2.3 kg to analyze and study the dynamic behaviour of the system with respect to bearing number.


### 3.3.1. Dynamic Orbits Analysis

An examination of Figures8(a1, b1), 8(a2, b2), …, and 8(a6, b6) shows that when the bearing number is small, Λ = 2.3, 3.12, the behaviour of the journal and the rotor will be relatively regular. However, when Λ is increased to 3.84, behaviour of the journal and the rotor centre will become irregular and asymmetrical. When the bearing number is further increased to 3.86, the dynamic orbits return to relatively regular periodic motion. However, irregular motion reappears when Λ
=
4.55 and only becomes regular again when the bearing number is increased to 4.89.Figure 8
Dynamic orbits ((a1)–(a6)) of the rotor centre when the bearing numberΛ = 2.3, 3.12, 3.84, 3.86, 4.55, and 4.89, respectively. Dynamic orbits of the journal centre ((b1)–(b6)). Spectrum response diagrams for the rotor centre in the horizontal direction ((c1)–(c6)) and the vertical direction of rotation ((d1)–(d6)), with a rotor mass m
r = 2.3 kg.


### 3.3.2. Spectrum Analysis

Figures8(c1, d1), 8(c2, d2), …, and 8(c6, d6) show the dynamic responses of the rotor centre in the horizontal and vertical directions of rotation with different bearing numbers. The study shows that at bearing numbers of Λ = 2.3, 3.86, and 4.89 the rotor centre demonstrates T periodic motions. When Λ = 3.84, the spectrum response diagram shows aperiodic motion of the rotor centre, but when the bearing number is Λ = 3.12, the motion becomes 2T subharmonic motion. However, when Λ = 4.55, the rotor centre motion becomes aperiodic.


### 3.3.3. Bifurcation Analysis

As shown in Figure9, the bearing number Λ can be used as a primary analytic parameter for studying the behaviour of an air bearing system and has been set within a range of 1.0 to 5.5 according to actual operational conditions. When the bearing number is small, such as Λ
=
2.3, the rotor centre exhibits T periodic motion in both the horizontal and rotational directions and this can be verified by the appearance of a single dot on the Poincaré map as shown in Figure 10(a). When the bearing number is increased to Λ
=
3.12, the system will bifurcate and generate 2T periodic motion, clearly shown by the presence of two discrete dots on the map as shown in Figure 10(b). When the bearing number increases to 3.84
≤
Λ
<
3.86, the system motion becomes chaotic, clearly shown by the irregular discrete points generated on the map in Figure 10(c). When Λ
=
3.86, the chaotic motion becomes T periodic motion, as shown in Figure 10(d), and the same occurs in the interval of 3.86
≤
Λ
<
4.55; see Figure 9. In addition, for the interval of 4.55
≤
Λ
<
4.89, the system generates quasi-periodic motion, as shown in Figure 10(e). When Λ
=
4.89, the system assumes T-periodic motion that is relatively stable and this continues up to Λ
=
5.5; see Figure 10(f).Figure 9
Bifurcation diagrams: (a)X
2(n
T) and (b) Y
2(n
t) for the rotor centre with respect to different bearing number Λ, with a rotor mass m
r = 2.3 kg.
(a)
(b)Figure 10
Poincaré maps of the rotor centre with respect to different bearing numbersΛ = (a) 2.3, (b) 3.12, (c) 3.84, (d) 3.86, (e) 4.55, and (f) 4.89.
(a)
(b)
(c)
(d)
(e)
(f)With regard to the question about whether a change of bearing number will cause the generation of chaotic behaviour, in this paper, similarly, the maximum Lyapunov exponents are also used for determination and further verification. From Figures11(a), 11(b), 11(d), 11(e), and 11(f), it can be seen that when Λ = 2.3, 3.12, 3.86, 4.55, and 4.89, the exponents obtained can all be smaller than or equal to zero and the system motion is nonchaotic. In view of the complexity of the relationship between the motion of the rotor centre and bearing number, all motions and behaviours generated with respect to the interval of 1.0
≤
Λ
≤
5.5 are summarized in Table 5 and a full view of the maximum Lyapunov exponents for the system is shown in Figure 12. The result indicates that chaotic regions (exponents greater than zero) only occur when 3.84
≤
Λ
<
3.86, which agrees with the results shown in Table 5.Table 5
Dynamic behaviour of the rotor centre with different bearing number (1.0
≤
Λ
≤
5.5).
Λ
[1.0, 3.12)
[3.12, 3.84)
[3.84, 3.86)
[3.86, 4.55)
[4.55, 4.89)
[4.89, 5.5]
Behavior
T
2T
Chaos
2T
Quasi
TFigure 11
Maximum Lyapunov exponents of the rotor centre with respect to different bearing numberΛ = (a) 2.3, (b) 3.12, (c) 3.84, (d) 3.86, (e) 4.55, and (f) 4.89.
(a)
(b)
(c)
(d)
(e)
(f)Figure 12
Full view of maximum Lyapunov exponents of the rotor centre with respect to different rotor-bearing numbers.


## 3.3.1. Dynamic Orbits Analysis

An examination of Figures8(a1, b1), 8(a2, b2), …, and 8(a6, b6) shows that when the bearing number is small, Λ = 2.3, 3.12, the behaviour of the journal and the rotor will be relatively regular. However, when Λ is increased to 3.84, behaviour of the journal and the rotor centre will become irregular and asymmetrical. When the bearing number is further increased to 3.86, the dynamic orbits return to relatively regular periodic motion. However, irregular motion reappears when Λ
=
4.55 and only becomes regular again when the bearing number is increased to 4.89.Figure 8
Dynamic orbits ((a1)–(a6)) of the rotor centre when the bearing numberΛ = 2.3, 3.12, 3.84, 3.86, 4.55, and 4.89, respectively. Dynamic orbits of the journal centre ((b1)–(b6)). Spectrum response diagrams for the rotor centre in the horizontal direction ((c1)–(c6)) and the vertical direction of rotation ((d1)–(d6)), with a rotor mass m
r = 2.3 kg.


## 3.3.2. Spectrum Analysis

Figures8(c1, d1), 8(c2, d2), …, and 8(c6, d6) show the dynamic responses of the rotor centre in the horizontal and vertical directions of rotation with different bearing numbers. The study shows that at bearing numbers of Λ = 2.3, 3.86, and 4.89 the rotor centre demonstrates T periodic motions. When Λ = 3.84, the spectrum response diagram shows aperiodic motion of the rotor centre, but when the bearing number is Λ = 3.12, the motion becomes 2T subharmonic motion. However, when Λ = 4.55, the rotor centre motion becomes aperiodic.


## 3.3.3. Bifurcation Analysis

As shown in Figure9, the bearing number Λ can be used as a primary analytic parameter for studying the behaviour of an air bearing system and has been set within a range of 1.0 to 5.5 according to actual operational conditions. When the bearing number is small, such as Λ
=
2.3, the rotor centre exhibits T periodic motion in both the horizontal and rotational directions and this can be verified by the appearance of a single dot on the Poincaré map as shown in Figure 10(a). When the bearing number is increased to Λ
=
3.12, the system will bifurcate and generate 2T periodic motion, clearly shown by the presence of two discrete dots on the map as shown in Figure 10(b). When the bearing number increases to 3.84
≤
Λ
<
3.86, the system motion becomes chaotic, clearly shown by the irregular discrete points generated on the map in Figure 10(c). When Λ
=
3.86, the chaotic motion becomes T periodic motion, as shown in Figure 10(d), and the same occurs in the interval of 3.86
≤
Λ
<
4.55; see Figure 9. In addition, for the interval of 4.55
≤
Λ
<
4.89, the system generates quasi-periodic motion, as shown in Figure 10(e). When Λ
=
4.89, the system assumes T-periodic motion that is relatively stable and this continues up to Λ
=
5.5; see Figure 10(f).Figure 9
Bifurcation diagrams: (a)X
2(n
T) and (b) Y
2(n
t) for the rotor centre with respect to different bearing number Λ, with a rotor mass m
r = 2.3 kg.
(a)
(b)Figure 10
Poincaré maps of the rotor centre with respect to different bearing numbersΛ = (a) 2.3, (b) 3.12, (c) 3.84, (d) 3.86, (e) 4.55, and (f) 4.89.
(a)
(b)
(c)
(d)
(e)
(f)With regard to the question about whether a change of bearing number will cause the generation of chaotic behaviour, in this paper, similarly, the maximum Lyapunov exponents are also used for determination and further verification. From Figures11(a), 11(b), 11(d), 11(e), and 11(f), it can be seen that when Λ = 2.3, 3.12, 3.86, 4.55, and 4.89, the exponents obtained can all be smaller than or equal to zero and the system motion is nonchaotic. In view of the complexity of the relationship between the motion of the rotor centre and bearing number, all motions and behaviours generated with respect to the interval of 1.0
≤
Λ
≤
5.5 are summarized in Table 5 and a full view of the maximum Lyapunov exponents for the system is shown in Figure 12. The result indicates that chaotic regions (exponents greater than zero) only occur when 3.84
≤
Λ
<
3.86, which agrees with the results shown in Table 5.Table 5
Dynamic behaviour of the rotor centre with different bearing number (1.0
≤
Λ
≤
5.5).
Λ
[1.0, 3.12)
[3.12, 3.84)
[3.84, 3.86)
[3.86, 4.55)
[4.55, 4.89)
[4.89, 5.5]
Behavior
T
2T
Chaos
2T
Quasi
TFigure 11
Maximum Lyapunov exponents of the rotor centre with respect to different bearing numberΛ = (a) 2.3, (b) 3.12, (c) 3.84, (d) 3.86, (e) 4.55, and (f) 4.89.
(a)
(b)
(c)
(d)
(e)
(f)Figure 12
Full view of maximum Lyapunov exponents of the rotor centre with respect to different rotor-bearing numbers.


## 4. Conclusion

The objective of this study was an analysis of the dynamic behaviour of a flexible rotor supported by a front opposed-hemispherical spiral-grooved air bearing system. The Finite Difference Method, Perturbation Method, and a hybrid method were used to solve the pressure distribution at the highest nonlinearity in the system, after which dynamic equations of the flexible rotor centre were used to obtain the orbits of the rotor centre. Analysis was then conducted on the orbit data to generate spectrum diagrams, Poincaré maps, bifurcation diagrams, and maximum Lyapunov exponents. This showed that the orbits of the rotor centre change along with the mass and the bearing number and synchronously generate complicated motion in the horizontal and vertical directions and include normal and subharmonic vibration, as well as periodic, quasi-periodic, and chaotic motion. The hybrid method applied in this paper yields better accuracy and precision for the verification of values than the Finite Difference Method and Perturbation Method. Accurate solutions can be obtained without the need for very precise position and time step determination. Findings with respect to the motion and behaviour of the rotor revealed that the rotor orbits demonstrate unstable nonlinear chaotic behaviour when the mass of the rotor is high but the system is relatively stable at the other intervals and ranges. The bearing number affects the orbits of the rotor which show chaotic behaviour in the interval3.84
≤
Λ
<
3.86 while the orbits in the other intervals and ranges show relatively organized and stable motion. The results obtained in this study can be used as a basis for future bearing system design and the prevention of instability.



---

*Source: 8985040-2017-02-14.xml*