# Certain Type Turbofan Engine Whole Vibration Model with Support Looseness Fault and Casing Response Characteristics

Authors: H. F. Wang; G. Chen
Journal: Shock and Vibration
 (2014)

---

## Abstract

Support looseness fault is a type of common fault in aeroengine. Serious looseness fault would emerge under larger unbalanced force, which would cause excessive vibration and even lead to rubbing fault, so it is important to analyze and recognize looseness fault effectively. In this paper, based on certain type turbofan engine structural features, a rotor-support-casing whole model for certain type turbofan aeroengine is established. The rotor and casing systems are modeled by means of the finite element beam method; the support systems are modeled by lumped-mass model; the support looseness fault model is also introduced. The coupled system response is obtained by numerical integral method. In this paper, based on the casing acceleration signals, the impact characteristics of symmetrical stiffness and asymmetric stiffness models are analyzed, finding that the looseness fault would lead to the longitudinal asymmetrical characteristics of acceleration time domain wave and the multiple frequency characteristics, which is consistent with the real trial running vibration signals. Asymmetric stiffness looseness model is verified to be fit for aeroengine looseness fault model.


---

## Body

## 1. Introduction

Looseness fault is a common fault in rotating machinery, which is caused by low installation quality or long-term vibration. When looseness fault exists in rotor-support-casing system, rotor will be lifted up periodically; if the unbalanced force in rotor system is greater than the gravitational force, it will cause severe vibration. Serious looseness fault may lead to rub-impact fault between the rotor and the stator. Therefore, effective identification of the looseness fault is of great significance.At present, many scholars have carried out a few works on the looseness fault, mostly based on a lumped-mass model. Muszynska and Goldman [1] studied one-lateral-mode unbalanced and radical side-loaded rotor with either a loose pedestal or occasional rotor-to-stator rubbing, which exhibited regular periodic vibrations of synchronous (1×) and subsynchronous (1/2×, 1/3×,…) orders, as well as chaotic vibration patterns of the rotor, all accompanied by higher harmonics. Goldman and Muszynska [2] developed an analytical algorithm for investigating local nonlinear effects in rotor systems. They used a specially developed variable transformation that smoothes discontinuities, and then they applied an averaging technique. Their results showed good agreement with experimentally observed typical behaviors and orbits of rubbing rotors. In the last two papers, the effects of pedestal looseness on the system response were also studied. Muszynska [3] studied a mathematical model of this phenomenon for loose rotating parts in rotating machines. Chu and Tang [4] established a nonlinear mathematical model, containing stiffness and considering damping forces with trilinear forms. The shooting method was used to obtain the periodic solutions of the system. Stability of these periodic solutions was analyzed by using the Floquet theory. Period-doubling bifurcation and Neimark-Sacker bifurcation were found. Finally, the governing equations were integrated by the fourth-order Runge-Kutta method. Three kinds of routes to or out of chaos, that is, period-to-chaos, quasiperiodic route, and intermittence, were found. Chen [5] established a new rotor-ball bearing-stator coupling dynamics model and the nonlinear dynamic response characteristics of rotor-ball bearing-stator system under unbalance and looseness coupling faults were obtained. Lu and Chu [6] studied the looseness fault of rotor system by experiments. The characteristics of multiple frequency and frequency division were found. Ji and Zu [7] analyzed the free and forced vibration of a nonlinear bearing system to illustrate the nonlinear effect on the free and forced vibrations of the system by the method of multiple scales. He et al. [8] and Lee and Choi [9] studied fault diagnosis of rotor systems with pedestal looseness with different analytical methods, such as genetic algorithm and Hilbert-Huang transform. Reddy and Sekhar [10] used artificial neural networks to identify unbalance and looseness. Ma et al. [11] set up a mechanical model of looseness of fastening bolt on the bearing pedestal and investigated the dynamic characteristics of rotor by using the nonlinear oil-film model.In recent years, the finite element method has been used to establish looseness fault model by domestic and foreign scholars. Behzad and Asayeshthe [12] proposed a finite element method for studying the effects of loose rotating disks on the rotor-bearing systems’ response. Lee et al. [13] studied a generalized finite element modeling method of a rotor-bearing system by using the state-space Newmark method based on the average velocity concept. Ma et al. [14] established a finite element model of a rotor system with pedestal looseness stemming from a loosened bolt and the effects of the looseness parameters on its dynamic characteristics were analyzed.However, at present, there are a few works on support looseness fault analysis in aeroengine. Due to the low bearing stiffness in aeroengine, the wide use of thin-walled structure in rotor and casing, its great flexibility, and the casing acceleration signal are important to distinguish aeroengine faults, so aimed at certain type turbofan engine, it is of great significance to establish a whole model of a rotor-bearing-casing coupling system and make research on the acceleration response characteristics of the looseness fault.In this paper, a finite element model of certain type turbofan engine with looseness fault is established. The looseness fault between the bearing and the casing in the model is studied. Numerical integration method is used to solve the response of the system. The response characteristics of the casing acceleration from the time domain and frequency domain are analyzed. Finally, the simulation is compared with the real trial running vibration signals of the turbofan engine, and connector looseness characteristics are analyzed. It is proved that the aeroengine support looseness fault model is valid.


## 2. A Dynamic Model for Certain Type Turbofan Engine

### 2.1. The Structure Sketch Map for Certain Type Turbofan Engine

Rotor-bearing-casing structure sketch map for certain type turbofan engine is shown in Figure1. The symbols in Figure 1 are described as follows: P
1 denotes fan disk; P
2 denotes motor disk; P
3 denotes compressor disk; P
4 denotes the first turbine disk; P
5 denotes the second turbine disk; C
1 denotes the intermediate casing; G
1 denotes the gear coupling between the fan shaft and the transmission shaft; G
2 denotes the gear coupling between the transmission shaft and the compressor shaft; G
3 denotes the gear coupling between the compressor shaft and the turbine shaft; S
1 denotes the fan support point bearings; S
2 denotes the compressor front support bearings; S
3 denotes the compressor rear support bearings; S
4 denotes the turbine support bearings; I
1 denotes the front installation node; I
2 denotes the back installation node; k
g is the mesh stiffness of a gear pump; k
f
1, k
f
2, k
f
3, and k
f
4 are the support stiffness between the rotor and the casing; and k
c is the connection stiffness between the casing and the base.Figure 1
Rotor-bearing-casing model sketch map of a type of aeroengine (unit: mm).The rotor-support-casing coupling dynamic model for a turbofan engine is established. Its concrete methods for modeling are described as follows.


### 2.2. Dynamics Model

#### 2.2.1. Rotor Model

The rotor is composed of a number of rigid rotating disks and elastic shaft. The elastic shaft is discreted by beam element model, and the shear deformations and gyroscopic moments and inertias of shaft are considered; the rigid rotating disk is modeled as discrete lumped mass, and the mass, inertias, and gyroscopic moments are considered. All parts are coupled by nonlinear forces and moments. The finite element rotor dynamic model is shown in Figure2.Figure 2
The finite element rotor dynamic model.The symbols in Figure2 are described as follows: for each rotor shaft beam element, E, I, G, μ, L, ρ, and A, respectively, are the elastic modulus, moment of inertia, shear modulus, Poisson’s ratio, shaft length, shaft density, and shaft cross-sectional area. P
i is the rigid rotating disks; F
x
i, F
y
i, M
x
i, and M
y
i are forces and moments which act on the ith node of rotor in the X and Y directions. X
Y
Z in Figure 3 is a fixed coordinate system; under deformation state, any cross-section’s position relative to the fixed coordinate system is determined as follows: elastic shaft centerline position is determined by the displacement x in the X direction and the displacement y in the Y direction; cross-section orientation ϕ is determined by the angle around the x-axis and the angle Ψ around the y-axis; in addition, the cross-section also rotates around the z-axis.Figure 3
Rotor-casing support.(1) Motion Equations of Rigid Disk Element. Suppose the circular disk’s mass is m
d, diameteral moment of inertia is J
d
d, and polar moment of inertia is J
P
d; polar moment of inertia J
P
d
=
2
J
d
d, andthe second term refers to “disk,” shown as P
1
,
P
2
,
…
,
P
M; ω is the disk’s rotating speed. Through the Lagrange equation we can obtain the motion equations of rigid disk relative to the fixed coordinate system:
(1)
M
T
d
+
M
R
d
q
¨
d
-
ω
G
d
q
˙
d
=
Q
d
,
where Q
d is the generalized external force vector; M
T
d and M
R
d are, respectively, the mass matrix and the mass inertia matrix; G
d is gyroscopic matrix; and q is generalized displacement vector, q
=
[
x
,
y
,
ϕ
,
ψ
]. Then,
(2)
M
T
d
=
m
d
0
0
0
0
m
d
0
0
0
0
0
0
0
0
0
0
,
M
R
d
=
0
0
0
0
0
0
0
0
0
0
J
d
d
0
0
0
0
J
d
d
,
G
d
=
0
0
0
0
0
0
0
0
0
0
0
-
J
p
d
0
0
J
p
d
0
.(2) Motion Equations of Beam Element. In this paper, the beam element is adopted, and each beam element has 2 nodes and 8 DOFs, and each node has 4 DOFs, which are, respectively, the displacements x and y in the X and Y directions, and rotating angles ϕ and Ψ around the x-axis and y-axis. The cross-section’s displacement of the element is the function of the time and the position along the element axis. The generalized displacement vector of the element is q
e
(
t
)
=
[
q
1
e
q
2
e
q
3
e
q
4
e
q
5
e
q
6
e
q
7
e
q
8
e
]
T. Through the Lagrange equation, the motion equations of beam element relative to the fixed coordinate system can be obtained:
(3)
M
T
e
+
M
R
e
q
¨
e
+
-
ω
G
e
q
˙
e
+
K
B
e
-
K
A
e
q
e
=
Q
e
,
where Q
e is the generalized external force vector; M
T
e and M
R
e are, respectively, the mass matrix and the mass inertia matrix; G
e is gyroscopic matrix; K
B
e are the element’s bending and shear stifness matrix; and K
A
e is the element’s tensile stifness matrix. These matrixes can be found in the literature [15].The rotor system’s motion equation through the element’s motion equations can be obtained, which is(4)
M
s
q
¨
s
+
C
s
-
ω
G
s
q
˙
s
+
K
s
q
s
=
Q
s
,
where Q
s is the generalized external force vector; M
s is the mass matrix; G
s is gyroscopic matrix; Κ
s is the stiffness matrix of the system; and C
s is the damping matrix of the system.In this paper,C
s is assumed to be proportional damping matrix, that is, C
s
=
α
0
M
s
+
α
1
Κ
s, of which, α
0 and α
1 are constants. Because the ith order damping ratio is
(5)
ξ
i
=
1
2
α
0
ω
i
+
α
1
ω
i
,
obviously, after any two natural frequencies and damping ratios of the rotor are acquired by modal experiment, α
0 and α
1 can be solved through (5), and the system’s damping ratio matrix C
s can be obtained.


#### 2.2.2. Casing Model

In general, the casing finite element model can be constructed by the beam element, conical shell element, and surface shell element. However the casing is a shell structure, and its vibration modal has many forms whose circumferential wave number is0,1
,
2
,
…, and so on. However, when the casing couples with the rotor, only the modal of wave number 1 appears, and the cross-section of rotor does not deform. Therefore, the casing can be constructed by using the common beam element which does not rotate [15].In this paper, the casing is modeled as a beam which does not rotate, like the rotor modeling, by means of finite element modeling method the differential equation of motion of casing is(6)
M
c
q
¨
c
+
C
c
q
˙
c
+
K
c
q
c
=
Q
c
,
where Q
c is the generalized excitation force vector; M
c is the mass matrix of casing; Κ
c is the stiffness matrix of casing; C
c is the damper matrix of casing; likewise, C
s is assumed to be proportional damping matrix.


#### 2.2.3. Discrete Support Model

(1) The Support Connection between Rotor and Casing. Each support connection between rotor and casing R
C
i (i
=
1,2
,
…
,
N) includes ball bearing and bearing housing. Suppose that m
w
i is the outer mass of bearing; m
b
i is the mass of bearing housing; k
t
i is the support stiffness between outer race and bearing house; c
t
i is the damping coefficient between outer race and the bearing housing; and k
f
i, c
f
i are, respectively, the support stiffness and damping coefficient between the casing and the bearing housing. As shown in Figure 3, F
y
R
i and F
x
R
i are the force of rotor acting on the support R
C
i and the F
y
C
i and F
x
C
i are the force of casing acting on support R
C
i. It is assumed that support R
C
i is connected with the mth node of rotor and the nth node of casing. In this paper, the looseness fault between the bearing house and casing is considered.In this paper, the bearing outer ring is supposed to be fixed on the bearing housing, and the inner ring is fixed on the rotating shaft. Assume that the displacements of themth node of rotor are x
R
m and y
R
m; let x
=
x
R
m
-
x
w
i, y
=
y
R
m
-
y
w
i; according to Chen et al. [16], the ball force can be expressed as
(7)
F
x
R
i
=
∑
j
=
1
N
C
b
x
cos
⁡
θ
j
+
y
sin
θ
j
-
r
0
3
/
2
×
H
x
cos
⁡
θ
j
+
y
sin
θ
j
-
r
0
cos
⁡
θ
j
,
F
y
R
i
=
∑
j
=
1
N
C
b
x
cos
⁡
θ
j
+
y
sin
θ
j
-
r
0
3
/
2
×
H
x
cos
⁡
θ
j
+
y
sin
θ
j
-
r
0
sin
θ
j
.In the formula,C
b is the Hertzian contact stiffness and it can be obtained from the nonlinear Hertzian contact elastic analysis of the inner race, the outer race and the balls H
(
·
) are the Heaviside function; when the function independent variable is greater than 0, function value is 1; otherwise it is 0. θ
j is the jth ball’s angle position, that is, θ
j
=
ω
Cage
×
t
+
(
2
π
/
N
b
)
(
j
-
1
), j
=
1,2
,
…
,
N
b, where N
b is the number of balls. ω
Cage is the cage’s rotating speed. Suppose that the outer race radius is R, the inner ring radius is r, ω
Cage
=
(
ω
×
r
)
/
(
R
+
r
), and ω is shaft rotating angular velocity.Therefore, the bearing outer race’s differential equation of motion is(8)
m
w
i
x
¨
w
i
+
k
t
i
x
w
i
-
x
b
i
+
F
d
x
i
=
F
x
R
i
,
m
w
i
y
¨
w
i
+
k
t
i
y
w
i
-
y
b
i
+
F
d
y
i
=
F
y
R
i
-
m
w
i
g
,
i
=
1,2
,
…
,
N
,
where F
d
x
i and F
d
y
i are damping forces; if the viscous damping is considered, then
(9)
F
d
x
i
=
c
t
i
x
˙
w
i
-
x
˙
b
i
,
F
d
y
i
=
c
t
i
y
˙
w
i
-
y
˙
b
i
.(2) The Coupling Connection between Rotor and Rotor. Assume that the coupling connection R
R
C
k (k
=
1,2
,
…
,
N) connects the ith node of left rotor and the jth node of right rotor, the radial stiffness of coupling is k
R
r, the angular stiffness is k
R
α, the radial damping is c
R
r, and the angular damping is c
R
α. Let the generalized displacements of the ith node of the left rotor be x
R
L
i
,
y
R
L
i
,
ϕ
R
L
i
,
ψ
R
L
i; the velocities are x
˙
R
L
i
,
y
˙
R
L
i
,
ϕ
˙
R
L
i
,
ψ
˙
R
L
i; let the jth node displacements of the right rotor be x
R
R
j
,
y
R
R
j
,
ϕ
R
R
j
,
ψ
R
R
j; the velocities are x
˙
R
R
j
,
y
˙
R
R
j
,
ϕ
˙
R
R
j
,
ψ
˙
R
R
j. Then, the forces and the moments which act on the ith node of left rotor F
R
x
i, F
R
y
i, M
R
x
i, M
R
y
i and the forces and the moments F
R
x
j, F
R
y
j, M
R
x
j, M
R
y
j for the jth node of right rotor are
(10)
F
R
x
i
=
k
R
r
x
R
R
j
-
x
R
L
i
+
c
R
r
x
˙
R
R
j
-
x
˙
R
L
i
,
F
R
y
i
=
k
R
r
y
R
R
j
-
y
R
L
i
+
c
R
r
y
˙
R
R
j
-
y
˙
R
L
i
,
M
R
x
i
=
k
R
α
ϕ
R
R
j
-
ϕ
R
L
i
+
c
R
α
ϕ
˙
R
R
j
-
ϕ
˙
R
L
i
,
M
R
y
i
=
k
R
α
ψ
R
R
j
-
ψ
R
L
i
+
c
R
α
ψ
˙
R
R
j
-
ψ
˙
R
L
i
,
F
R
x
j
=
-
F
R
x
i
,
F
R
y
j
=
-
F
R
y
i
,
M
R
x
j
=
-
M
R
x
i
,
M
R
y
j
=
-
M
R
y
i
.(3) The Bearing Housing Looseness Model. Generally, symmetric stiffness and asymmetrical stiffness models [4, 5] are considered in the looseness fault model. But presently the looseness stiffness model has no universal definition. In this paper, the simulation analysis is made in two conditions. k
f
0 is assumed to be the equivalent stiffness between the bearing housing and the casing. Under the condition of relative displacement, piecewise linear stiffness between the bearing housing and the casing is considered in this paper; δ is assumed to be the looseness amounts of the bearing housing. The piecewise nonlinear stiffness k
f can be expressed as follows.Asymmetric stiffness:(11)
k
f
=
k
f
0
(
x
b
-
x
c
>
δ
)
k
f
0
5
(
0
≤
x
b
-
x
c
≤
δ
)
5
k
f
0
x
b
-
x
c
<
0
;
symmetric stiffness:
(12)
k
f
=
k
f
0
(
x
b
-
x
c
>
δ
)
k
f
0
5
(
0
≤
x
b
-
x
c
≤
δ
)
k
f
0
x
b
-
x
c
<
0
,
where x
b is the displacement of the bearing housing and x
c is the displacement of the casing; the stiffness directions in x direction and y direction are the same. The formulas are appied in horizontal and vertical directions.


#### 2.2.4. Elastic Support (Installation Node) between Casing and Base

The elastic connection between casing and baseC
B
k (k
=
1,2
,
…
,
N) connects the ith node of casing and the base, the support stiffness is k
C, and the damping is c
C. Let the displacements of the ith node of casing be x
C
i
,
y
C
i and the velocities x
˙
C
i
,
y
˙
C
i; then the forces acting on the ith node of casing F
C
x
i and F
C
y
i are
(13)
F
C
x
i
=
-
k
C
x
C
i
-
c
C
x
˙
C
i
,
F
C
y
i
=
-
k
C
y
C
i
-
c
C
y
˙
C
i
.


### 2.3. Solution of Finite Element Rotor-Support-Casing Coupling Dynamic Model

Because the number of DOFs in the finite element rotor-support-casing coupling dynamic model is very large, in addition, there are a lot of strong nonlinear factors, in this paper a combination time integration approach is used to solve system’s responses, which combines the implicit Newmark-β method and the improved explicit Newmark-β method (Zhai [17]). This finite element rotor-support coupling dynamic solution procedure is shown in Figure 4.Figure 4
Solving flow for rotor-support-casing coupling dynamics.The merits of combination methods are described as follows: (1) the Newmark-β method is used to solve the rotor and casing finite element model, and the Zhai method is used to solve the lumped parameter support models; therefore, the combination method is very fit for solving coupling system; (2) the enormous matrix need not be formed, and the matrixes of the rotors and the casings need not be combined into a large matrix; therefore, the computation efficiency is improved greatly.


## 2.1. The Structure Sketch Map for Certain Type Turbofan Engine

Rotor-bearing-casing structure sketch map for certain type turbofan engine is shown in Figure1. The symbols in Figure 1 are described as follows: P
1 denotes fan disk; P
2 denotes motor disk; P
3 denotes compressor disk; P
4 denotes the first turbine disk; P
5 denotes the second turbine disk; C
1 denotes the intermediate casing; G
1 denotes the gear coupling between the fan shaft and the transmission shaft; G
2 denotes the gear coupling between the transmission shaft and the compressor shaft; G
3 denotes the gear coupling between the compressor shaft and the turbine shaft; S
1 denotes the fan support point bearings; S
2 denotes the compressor front support bearings; S
3 denotes the compressor rear support bearings; S
4 denotes the turbine support bearings; I
1 denotes the front installation node; I
2 denotes the back installation node; k
g is the mesh stiffness of a gear pump; k
f
1, k
f
2, k
f
3, and k
f
4 are the support stiffness between the rotor and the casing; and k
c is the connection stiffness between the casing and the base.Figure 1
Rotor-bearing-casing model sketch map of a type of aeroengine (unit: mm).The rotor-support-casing coupling dynamic model for a turbofan engine is established. Its concrete methods for modeling are described as follows.


## 2.2. Dynamics Model

### 2.2.1. Rotor Model

The rotor is composed of a number of rigid rotating disks and elastic shaft. The elastic shaft is discreted by beam element model, and the shear deformations and gyroscopic moments and inertias of shaft are considered; the rigid rotating disk is modeled as discrete lumped mass, and the mass, inertias, and gyroscopic moments are considered. All parts are coupled by nonlinear forces and moments. The finite element rotor dynamic model is shown in Figure2.Figure 2
The finite element rotor dynamic model.The symbols in Figure2 are described as follows: for each rotor shaft beam element, E, I, G, μ, L, ρ, and A, respectively, are the elastic modulus, moment of inertia, shear modulus, Poisson’s ratio, shaft length, shaft density, and shaft cross-sectional area. P
i is the rigid rotating disks; F
x
i, F
y
i, M
x
i, and M
y
i are forces and moments which act on the ith node of rotor in the X and Y directions. X
Y
Z in Figure 3 is a fixed coordinate system; under deformation state, any cross-section’s position relative to the fixed coordinate system is determined as follows: elastic shaft centerline position is determined by the displacement x in the X direction and the displacement y in the Y direction; cross-section orientation ϕ is determined by the angle around the x-axis and the angle Ψ around the y-axis; in addition, the cross-section also rotates around the z-axis.Figure 3
Rotor-casing support.(1) Motion Equations of Rigid Disk Element. Suppose the circular disk’s mass is m
d, diameteral moment of inertia is J
d
d, and polar moment of inertia is J
P
d; polar moment of inertia J
P
d
=
2
J
d
d, andthe second term refers to “disk,” shown as P
1
,
P
2
,
…
,
P
M; ω is the disk’s rotating speed. Through the Lagrange equation we can obtain the motion equations of rigid disk relative to the fixed coordinate system:
(1)
M
T
d
+
M
R
d
q
¨
d
-
ω
G
d
q
˙
d
=
Q
d
,
where Q
d is the generalized external force vector; M
T
d and M
R
d are, respectively, the mass matrix and the mass inertia matrix; G
d is gyroscopic matrix; and q is generalized displacement vector, q
=
[
x
,
y
,
ϕ
,
ψ
]. Then,
(2)
M
T
d
=
m
d
0
0
0
0
m
d
0
0
0
0
0
0
0
0
0
0
,
M
R
d
=
0
0
0
0
0
0
0
0
0
0
J
d
d
0
0
0
0
J
d
d
,
G
d
=
0
0
0
0
0
0
0
0
0
0
0
-
J
p
d
0
0
J
p
d
0
.(2) Motion Equations of Beam Element. In this paper, the beam element is adopted, and each beam element has 2 nodes and 8 DOFs, and each node has 4 DOFs, which are, respectively, the displacements x and y in the X and Y directions, and rotating angles ϕ and Ψ around the x-axis and y-axis. The cross-section’s displacement of the element is the function of the time and the position along the element axis. The generalized displacement vector of the element is q
e
(
t
)
=
[
q
1
e
q
2
e
q
3
e
q
4
e
q
5
e
q
6
e
q
7
e
q
8
e
]
T. Through the Lagrange equation, the motion equations of beam element relative to the fixed coordinate system can be obtained:
(3)
M
T
e
+
M
R
e
q
¨
e
+
-
ω
G
e
q
˙
e
+
K
B
e
-
K
A
e
q
e
=
Q
e
,
where Q
e is the generalized external force vector; M
T
e and M
R
e are, respectively, the mass matrix and the mass inertia matrix; G
e is gyroscopic matrix; K
B
e are the element’s bending and shear stifness matrix; and K
A
e is the element’s tensile stifness matrix. These matrixes can be found in the literature [15].The rotor system’s motion equation through the element’s motion equations can be obtained, which is(4)
M
s
q
¨
s
+
C
s
-
ω
G
s
q
˙
s
+
K
s
q
s
=
Q
s
,
where Q
s is the generalized external force vector; M
s is the mass matrix; G
s is gyroscopic matrix; Κ
s is the stiffness matrix of the system; and C
s is the damping matrix of the system.In this paper,C
s is assumed to be proportional damping matrix, that is, C
s
=
α
0
M
s
+
α
1
Κ
s, of which, α
0 and α
1 are constants. Because the ith order damping ratio is
(5)
ξ
i
=
1
2
α
0
ω
i
+
α
1
ω
i
,
obviously, after any two natural frequencies and damping ratios of the rotor are acquired by modal experiment, α
0 and α
1 can be solved through (5), and the system’s damping ratio matrix C
s can be obtained.


### 2.2.2. Casing Model

In general, the casing finite element model can be constructed by the beam element, conical shell element, and surface shell element. However the casing is a shell structure, and its vibration modal has many forms whose circumferential wave number is0,1
,
2
,
…, and so on. However, when the casing couples with the rotor, only the modal of wave number 1 appears, and the cross-section of rotor does not deform. Therefore, the casing can be constructed by using the common beam element which does not rotate [15].In this paper, the casing is modeled as a beam which does not rotate, like the rotor modeling, by means of finite element modeling method the differential equation of motion of casing is(6)
M
c
q
¨
c
+
C
c
q
˙
c
+
K
c
q
c
=
Q
c
,
where Q
c is the generalized excitation force vector; M
c is the mass matrix of casing; Κ
c is the stiffness matrix of casing; C
c is the damper matrix of casing; likewise, C
s is assumed to be proportional damping matrix.


### 2.2.3. Discrete Support Model

(1) The Support Connection between Rotor and Casing. Each support connection between rotor and casing R
C
i (i
=
1,2
,
…
,
N) includes ball bearing and bearing housing. Suppose that m
w
i is the outer mass of bearing; m
b
i is the mass of bearing housing; k
t
i is the support stiffness between outer race and bearing house; c
t
i is the damping coefficient between outer race and the bearing housing; and k
f
i, c
f
i are, respectively, the support stiffness and damping coefficient between the casing and the bearing housing. As shown in Figure 3, F
y
R
i and F
x
R
i are the force of rotor acting on the support R
C
i and the F
y
C
i and F
x
C
i are the force of casing acting on support R
C
i. It is assumed that support R
C
i is connected with the mth node of rotor and the nth node of casing. In this paper, the looseness fault between the bearing house and casing is considered.In this paper, the bearing outer ring is supposed to be fixed on the bearing housing, and the inner ring is fixed on the rotating shaft. Assume that the displacements of themth node of rotor are x
R
m and y
R
m; let x
=
x
R
m
-
x
w
i, y
=
y
R
m
-
y
w
i; according to Chen et al. [16], the ball force can be expressed as
(7)
F
x
R
i
=
∑
j
=
1
N
C
b
x
cos
⁡
θ
j
+
y
sin
θ
j
-
r
0
3
/
2
×
H
x
cos
⁡
θ
j
+
y
sin
θ
j
-
r
0
cos
⁡
θ
j
,
F
y
R
i
=
∑
j
=
1
N
C
b
x
cos
⁡
θ
j
+
y
sin
θ
j
-
r
0
3
/
2
×
H
x
cos
⁡
θ
j
+
y
sin
θ
j
-
r
0
sin
θ
j
.In the formula,C
b is the Hertzian contact stiffness and it can be obtained from the nonlinear Hertzian contact elastic analysis of the inner race, the outer race and the balls H
(
·
) are the Heaviside function; when the function independent variable is greater than 0, function value is 1; otherwise it is 0. θ
j is the jth ball’s angle position, that is, θ
j
=
ω
Cage
×
t
+
(
2
π
/
N
b
)
(
j
-
1
), j
=
1,2
,
…
,
N
b, where N
b is the number of balls. ω
Cage is the cage’s rotating speed. Suppose that the outer race radius is R, the inner ring radius is r, ω
Cage
=
(
ω
×
r
)
/
(
R
+
r
), and ω is shaft rotating angular velocity.Therefore, the bearing outer race’s differential equation of motion is(8)
m
w
i
x
¨
w
i
+
k
t
i
x
w
i
-
x
b
i
+
F
d
x
i
=
F
x
R
i
,
m
w
i
y
¨
w
i
+
k
t
i
y
w
i
-
y
b
i
+
F
d
y
i
=
F
y
R
i
-
m
w
i
g
,
i
=
1,2
,
…
,
N
,
where F
d
x
i and F
d
y
i are damping forces; if the viscous damping is considered, then
(9)
F
d
x
i
=
c
t
i
x
˙
w
i
-
x
˙
b
i
,
F
d
y
i
=
c
t
i
y
˙
w
i
-
y
˙
b
i
.(2) The Coupling Connection between Rotor and Rotor. Assume that the coupling connection R
R
C
k (k
=
1,2
,
…
,
N) connects the ith node of left rotor and the jth node of right rotor, the radial stiffness of coupling is k
R
r, the angular stiffness is k
R
α, the radial damping is c
R
r, and the angular damping is c
R
α. Let the generalized displacements of the ith node of the left rotor be x
R
L
i
,
y
R
L
i
,
ϕ
R
L
i
,
ψ
R
L
i; the velocities are x
˙
R
L
i
,
y
˙
R
L
i
,
ϕ
˙
R
L
i
,
ψ
˙
R
L
i; let the jth node displacements of the right rotor be x
R
R
j
,
y
R
R
j
,
ϕ
R
R
j
,
ψ
R
R
j; the velocities are x
˙
R
R
j
,
y
˙
R
R
j
,
ϕ
˙
R
R
j
,
ψ
˙
R
R
j. Then, the forces and the moments which act on the ith node of left rotor F
R
x
i, F
R
y
i, M
R
x
i, M
R
y
i and the forces and the moments F
R
x
j, F
R
y
j, M
R
x
j, M
R
y
j for the jth node of right rotor are
(10)
F
R
x
i
=
k
R
r
x
R
R
j
-
x
R
L
i
+
c
R
r
x
˙
R
R
j
-
x
˙
R
L
i
,
F
R
y
i
=
k
R
r
y
R
R
j
-
y
R
L
i
+
c
R
r
y
˙
R
R
j
-
y
˙
R
L
i
,
M
R
x
i
=
k
R
α
ϕ
R
R
j
-
ϕ
R
L
i
+
c
R
α
ϕ
˙
R
R
j
-
ϕ
˙
R
L
i
,
M
R
y
i
=
k
R
α
ψ
R
R
j
-
ψ
R
L
i
+
c
R
α
ψ
˙
R
R
j
-
ψ
˙
R
L
i
,
F
R
x
j
=
-
F
R
x
i
,
F
R
y
j
=
-
F
R
y
i
,
M
R
x
j
=
-
M
R
x
i
,
M
R
y
j
=
-
M
R
y
i
.(3) The Bearing Housing Looseness Model. Generally, symmetric stiffness and asymmetrical stiffness models [4, 5] are considered in the looseness fault model. But presently the looseness stiffness model has no universal definition. In this paper, the simulation analysis is made in two conditions. k
f
0 is assumed to be the equivalent stiffness between the bearing housing and the casing. Under the condition of relative displacement, piecewise linear stiffness between the bearing housing and the casing is considered in this paper; δ is assumed to be the looseness amounts of the bearing housing. The piecewise nonlinear stiffness k
f can be expressed as follows.Asymmetric stiffness:(11)
k
f
=
k
f
0
(
x
b
-
x
c
>
δ
)
k
f
0
5
(
0
≤
x
b
-
x
c
≤
δ
)
5
k
f
0
x
b
-
x
c
<
0
;
symmetric stiffness:
(12)
k
f
=
k
f
0
(
x
b
-
x
c
>
δ
)
k
f
0
5
(
0
≤
x
b
-
x
c
≤
δ
)
k
f
0
x
b
-
x
c
<
0
,
where x
b is the displacement of the bearing housing and x
c is the displacement of the casing; the stiffness directions in x direction and y direction are the same. The formulas are appied in horizontal and vertical directions.


### 2.2.4. Elastic Support (Installation Node) between Casing and Base

The elastic connection between casing and baseC
B
k (k
=
1,2
,
…
,
N) connects the ith node of casing and the base, the support stiffness is k
C, and the damping is c
C. Let the displacements of the ith node of casing be x
C
i
,
y
C
i and the velocities x
˙
C
i
,
y
˙
C
i; then the forces acting on the ith node of casing F
C
x
i and F
C
y
i are
(13)
F
C
x
i
=
-
k
C
x
C
i
-
c
C
x
˙
C
i
,
F
C
y
i
=
-
k
C
y
C
i
-
c
C
y
˙
C
i
.


## 2.2.1. Rotor Model

The rotor is composed of a number of rigid rotating disks and elastic shaft. The elastic shaft is discreted by beam element model, and the shear deformations and gyroscopic moments and inertias of shaft are considered; the rigid rotating disk is modeled as discrete lumped mass, and the mass, inertias, and gyroscopic moments are considered. All parts are coupled by nonlinear forces and moments. The finite element rotor dynamic model is shown in Figure2.Figure 2
The finite element rotor dynamic model.The symbols in Figure2 are described as follows: for each rotor shaft beam element, E, I, G, μ, L, ρ, and A, respectively, are the elastic modulus, moment of inertia, shear modulus, Poisson’s ratio, shaft length, shaft density, and shaft cross-sectional area. P
i is the rigid rotating disks; F
x
i, F
y
i, M
x
i, and M
y
i are forces and moments which act on the ith node of rotor in the X and Y directions. X
Y
Z in Figure 3 is a fixed coordinate system; under deformation state, any cross-section’s position relative to the fixed coordinate system is determined as follows: elastic shaft centerline position is determined by the displacement x in the X direction and the displacement y in the Y direction; cross-section orientation ϕ is determined by the angle around the x-axis and the angle Ψ around the y-axis; in addition, the cross-section also rotates around the z-axis.Figure 3
Rotor-casing support.(1) Motion Equations of Rigid Disk Element. Suppose the circular disk’s mass is m
d, diameteral moment of inertia is J
d
d, and polar moment of inertia is J
P
d; polar moment of inertia J
P
d
=
2
J
d
d, andthe second term refers to “disk,” shown as P
1
,
P
2
,
…
,
P
M; ω is the disk’s rotating speed. Through the Lagrange equation we can obtain the motion equations of rigid disk relative to the fixed coordinate system:
(1)
M
T
d
+
M
R
d
q
¨
d
-
ω
G
d
q
˙
d
=
Q
d
,
where Q
d is the generalized external force vector; M
T
d and M
R
d are, respectively, the mass matrix and the mass inertia matrix; G
d is gyroscopic matrix; and q is generalized displacement vector, q
=
[
x
,
y
,
ϕ
,
ψ
]. Then,
(2)
M
T
d
=
m
d
0
0
0
0
m
d
0
0
0
0
0
0
0
0
0
0
,
M
R
d
=
0
0
0
0
0
0
0
0
0
0
J
d
d
0
0
0
0
J
d
d
,
G
d
=
0
0
0
0
0
0
0
0
0
0
0
-
J
p
d
0
0
J
p
d
0
.(2) Motion Equations of Beam Element. In this paper, the beam element is adopted, and each beam element has 2 nodes and 8 DOFs, and each node has 4 DOFs, which are, respectively, the displacements x and y in the X and Y directions, and rotating angles ϕ and Ψ around the x-axis and y-axis. The cross-section’s displacement of the element is the function of the time and the position along the element axis. The generalized displacement vector of the element is q
e
(
t
)
=
[
q
1
e
q
2
e
q
3
e
q
4
e
q
5
e
q
6
e
q
7
e
q
8
e
]
T. Through the Lagrange equation, the motion equations of beam element relative to the fixed coordinate system can be obtained:
(3)
M
T
e
+
M
R
e
q
¨
e
+
-
ω
G
e
q
˙
e
+
K
B
e
-
K
A
e
q
e
=
Q
e
,
where Q
e is the generalized external force vector; M
T
e and M
R
e are, respectively, the mass matrix and the mass inertia matrix; G
e is gyroscopic matrix; K
B
e are the element’s bending and shear stifness matrix; and K
A
e is the element’s tensile stifness matrix. These matrixes can be found in the literature [15].The rotor system’s motion equation through the element’s motion equations can be obtained, which is(4)
M
s
q
¨
s
+
C
s
-
ω
G
s
q
˙
s
+
K
s
q
s
=
Q
s
,
where Q
s is the generalized external force vector; M
s is the mass matrix; G
s is gyroscopic matrix; Κ
s is the stiffness matrix of the system; and C
s is the damping matrix of the system.In this paper,C
s is assumed to be proportional damping matrix, that is, C
s
=
α
0
M
s
+
α
1
Κ
s, of which, α
0 and α
1 are constants. Because the ith order damping ratio is
(5)
ξ
i
=
1
2
α
0
ω
i
+
α
1
ω
i
,
obviously, after any two natural frequencies and damping ratios of the rotor are acquired by modal experiment, α
0 and α
1 can be solved through (5), and the system’s damping ratio matrix C
s can be obtained.


## 2.2.2. Casing Model

In general, the casing finite element model can be constructed by the beam element, conical shell element, and surface shell element. However the casing is a shell structure, and its vibration modal has many forms whose circumferential wave number is0,1
,
2
,
…, and so on. However, when the casing couples with the rotor, only the modal of wave number 1 appears, and the cross-section of rotor does not deform. Therefore, the casing can be constructed by using the common beam element which does not rotate [15].In this paper, the casing is modeled as a beam which does not rotate, like the rotor modeling, by means of finite element modeling method the differential equation of motion of casing is(6)
M
c
q
¨
c
+
C
c
q
˙
c
+
K
c
q
c
=
Q
c
,
where Q
c is the generalized excitation force vector; M
c is the mass matrix of casing; Κ
c is the stiffness matrix of casing; C
c is the damper matrix of casing; likewise, C
s is assumed to be proportional damping matrix.


## 2.2.3. Discrete Support Model

(1) The Support Connection between Rotor and Casing. Each support connection between rotor and casing R
C
i (i
=
1,2
,
…
,
N) includes ball bearing and bearing housing. Suppose that m
w
i is the outer mass of bearing; m
b
i is the mass of bearing housing; k
t
i is the support stiffness between outer race and bearing house; c
t
i is the damping coefficient between outer race and the bearing housing; and k
f
i, c
f
i are, respectively, the support stiffness and damping coefficient between the casing and the bearing housing. As shown in Figure 3, F
y
R
i and F
x
R
i are the force of rotor acting on the support R
C
i and the F
y
C
i and F
x
C
i are the force of casing acting on support R
C
i. It is assumed that support R
C
i is connected with the mth node of rotor and the nth node of casing. In this paper, the looseness fault between the bearing house and casing is considered.In this paper, the bearing outer ring is supposed to be fixed on the bearing housing, and the inner ring is fixed on the rotating shaft. Assume that the displacements of themth node of rotor are x
R
m and y
R
m; let x
=
x
R
m
-
x
w
i, y
=
y
R
m
-
y
w
i; according to Chen et al. [16], the ball force can be expressed as
(7)
F
x
R
i
=
∑
j
=
1
N
C
b
x
cos
⁡
θ
j
+
y
sin
θ
j
-
r
0
3
/
2
×
H
x
cos
⁡
θ
j
+
y
sin
θ
j
-
r
0
cos
⁡
θ
j
,
F
y
R
i
=
∑
j
=
1
N
C
b
x
cos
⁡
θ
j
+
y
sin
θ
j
-
r
0
3
/
2
×
H
x
cos
⁡
θ
j
+
y
sin
θ
j
-
r
0
sin
θ
j
.In the formula,C
b is the Hertzian contact stiffness and it can be obtained from the nonlinear Hertzian contact elastic analysis of the inner race, the outer race and the balls H
(
·
) are the Heaviside function; when the function independent variable is greater than 0, function value is 1; otherwise it is 0. θ
j is the jth ball’s angle position, that is, θ
j
=
ω
Cage
×
t
+
(
2
π
/
N
b
)
(
j
-
1
), j
=
1,2
,
…
,
N
b, where N
b is the number of balls. ω
Cage is the cage’s rotating speed. Suppose that the outer race radius is R, the inner ring radius is r, ω
Cage
=
(
ω
×
r
)
/
(
R
+
r
), and ω is shaft rotating angular velocity.Therefore, the bearing outer race’s differential equation of motion is(8)
m
w
i
x
¨
w
i
+
k
t
i
x
w
i
-
x
b
i
+
F
d
x
i
=
F
x
R
i
,
m
w
i
y
¨
w
i
+
k
t
i
y
w
i
-
y
b
i
+
F
d
y
i
=
F
y
R
i
-
m
w
i
g
,
i
=
1,2
,
…
,
N
,
where F
d
x
i and F
d
y
i are damping forces; if the viscous damping is considered, then
(9)
F
d
x
i
=
c
t
i
x
˙
w
i
-
x
˙
b
i
,
F
d
y
i
=
c
t
i
y
˙
w
i
-
y
˙
b
i
.(2) The Coupling Connection between Rotor and Rotor. Assume that the coupling connection R
R
C
k (k
=
1,2
,
…
,
N) connects the ith node of left rotor and the jth node of right rotor, the radial stiffness of coupling is k
R
r, the angular stiffness is k
R
α, the radial damping is c
R
r, and the angular damping is c
R
α. Let the generalized displacements of the ith node of the left rotor be x
R
L
i
,
y
R
L
i
,
ϕ
R
L
i
,
ψ
R
L
i; the velocities are x
˙
R
L
i
,
y
˙
R
L
i
,
ϕ
˙
R
L
i
,
ψ
˙
R
L
i; let the jth node displacements of the right rotor be x
R
R
j
,
y
R
R
j
,
ϕ
R
R
j
,
ψ
R
R
j; the velocities are x
˙
R
R
j
,
y
˙
R
R
j
,
ϕ
˙
R
R
j
,
ψ
˙
R
R
j. Then, the forces and the moments which act on the ith node of left rotor F
R
x
i, F
R
y
i, M
R
x
i, M
R
y
i and the forces and the moments F
R
x
j, F
R
y
j, M
R
x
j, M
R
y
j for the jth node of right rotor are
(10)
F
R
x
i
=
k
R
r
x
R
R
j
-
x
R
L
i
+
c
R
r
x
˙
R
R
j
-
x
˙
R
L
i
,
F
R
y
i
=
k
R
r
y
R
R
j
-
y
R
L
i
+
c
R
r
y
˙
R
R
j
-
y
˙
R
L
i
,
M
R
x
i
=
k
R
α
ϕ
R
R
j
-
ϕ
R
L
i
+
c
R
α
ϕ
˙
R
R
j
-
ϕ
˙
R
L
i
,
M
R
y
i
=
k
R
α
ψ
R
R
j
-
ψ
R
L
i
+
c
R
α
ψ
˙
R
R
j
-
ψ
˙
R
L
i
,
F
R
x
j
=
-
F
R
x
i
,
F
R
y
j
=
-
F
R
y
i
,
M
R
x
j
=
-
M
R
x
i
,
M
R
y
j
=
-
M
R
y
i
.(3) The Bearing Housing Looseness Model. Generally, symmetric stiffness and asymmetrical stiffness models [4, 5] are considered in the looseness fault model. But presently the looseness stiffness model has no universal definition. In this paper, the simulation analysis is made in two conditions. k
f
0 is assumed to be the equivalent stiffness between the bearing housing and the casing. Under the condition of relative displacement, piecewise linear stiffness between the bearing housing and the casing is considered in this paper; δ is assumed to be the looseness amounts of the bearing housing. The piecewise nonlinear stiffness k
f can be expressed as follows.Asymmetric stiffness:(11)
k
f
=
k
f
0
(
x
b
-
x
c
>
δ
)
k
f
0
5
(
0
≤
x
b
-
x
c
≤
δ
)
5
k
f
0
x
b
-
x
c
<
0
;
symmetric stiffness:
(12)
k
f
=
k
f
0
(
x
b
-
x
c
>
δ
)
k
f
0
5
(
0
≤
x
b
-
x
c
≤
δ
)
k
f
0
x
b
-
x
c
<
0
,
where x
b is the displacement of the bearing housing and x
c is the displacement of the casing; the stiffness directions in x direction and y direction are the same. The formulas are appied in horizontal and vertical directions.


## 2.2.4. Elastic Support (Installation Node) between Casing and Base

The elastic connection between casing and baseC
B
k (k
=
1,2
,
…
,
N) connects the ith node of casing and the base, the support stiffness is k
C, and the damping is c
C. Let the displacements of the ith node of casing be x
C
i
,
y
C
i and the velocities x
˙
C
i
,
y
˙
C
i; then the forces acting on the ith node of casing F
C
x
i and F
C
y
i are
(13)
F
C
x
i
=
-
k
C
x
C
i
-
c
C
x
˙
C
i
,
F
C
y
i
=
-
k
C
y
C
i
-
c
C
y
˙
C
i
.


## 2.3. Solution of Finite Element Rotor-Support-Casing Coupling Dynamic Model

Because the number of DOFs in the finite element rotor-support-casing coupling dynamic model is very large, in addition, there are a lot of strong nonlinear factors, in this paper a combination time integration approach is used to solve system’s responses, which combines the implicit Newmark-β method and the improved explicit Newmark-β method (Zhai [17]). This finite element rotor-support coupling dynamic solution procedure is shown in Figure 4.Figure 4
Solving flow for rotor-support-casing coupling dynamics.The merits of combination methods are described as follows: (1) the Newmark-β method is used to solve the rotor and casing finite element model, and the Zhai method is used to solve the lumped parameter support models; therefore, the combination method is very fit for solving coupling system; (2) the enormous matrix need not be formed, and the matrixes of the rotors and the casings need not be combined into a large matrix; therefore, the computation efficiency is improved greatly.


## 3. Looseness Fault Analysis

### 3.1. Dynamic Model Parameters

The finite element parameters of the rotor, the casing, and the connection parameters of rotor-casing-bearing system are shown in Tables1, 2, 3, 4, 5, 6, and 7.Table 1
Unit number of rotors and casings.
Fan rotor
Compressor rotor
Turbine rotor
Casing
11
10
11
24Table 2
Parameters of rotors.
Parameters
DiskP
1
DiskP
2
DiskP
3
DiskP
4
DiskP
5
Massm
p (kg)
3.88
1.41
5.17
10.28
10.28
Cross-polar inertiaJ
p
d (kg·m2)
0.03
0.003
0.03
0.05
0.05
Cross-equator inertiaJ
d
d (kg·m2)
0.015
0.0015
0.015
0.025
0.025
Elastic modulusE (Pa)
2.07
2.07
2.07
2.07
2.07
Poisson’s ratioμ
0.3
0.3
0.3
0.3
0.3
Densityρ (kg/m3)
7.8 × 103
7.8 × 103
7.8 × 103
7.8 × 103
7.8 × 103
Proportion damping ratioα
0
5
5
5
5
5
Proportion damping ratioα
1
1.35
1.35
1.35
1.35
1.35Table 3
Parameters of casings.
Wall thickness/mm
Elastic modulusE (Pa)
Densityρ (kg/m3)
Poisson’s ratioμ
Proportiondamping ratio α
0
Proportiondamping ratio α
1
15
2.07 × 1011
7.8 × 103
0.3
5
1.35 × 10−5Table 4
Parameters of ball bearing.
Rolling bearing
Outer raceway radiusR/mm
Inner raceway radiusr/mm
Ball numberN
b
Contact stiffnessC
b/(N/m3/2)
Bearing clearancer
0/um
Outer massm
w (kg)
Bearing carrier massm
b (kg)
B
1
39.5
29
13
12.4 × 109
0
2
10
B
2
39.5
29
13
12.4 × 109
0
2
10
B
3
32
17
14
11.9 × 109
0
2
10
B
4
32
17
14
11.9 × 109
0
2
10Table 5
Support parameters of rotor-casing.
Supports
Node of rotor
Casing (node)
k
t (N/m)
c
t (N·s/m)
k
f (N/m)
c
f (N·s/m)
R
C
1
3
2
1 × 108
2000
1 × 108
1000
R
C
2
1
9
1 × 108
2000
1 × 108
1000
R
C
3
11
16
1 × 108
2000
1 × 108
1000
R
C
4
8
22
1 × 108
2000
1 × 108
1000Table 6
Spring collection parameters of rotor-casing.
Collection
Node of rotor
Casing (node)
k
g
x (N/m)
c
g
x (N·s/m)
k
g
y (N/m)
c
g
y (N·s/m)
R
C
1
6
4
1 × 108
0
1 × 108
0Table 7
Collection parameters of casing-base.
Supports
Node of rotor
Casing (node)
k
t (N/m)
c
t (N·s/m)
k
f (N/m)
C
B
1
8
1 × 109
1 × 105
2000
0
C
B
2
23
1 × 109
1 × 105
2000
0


### 3.2. Calculation Condition

(1)
The looseness fault between the bearing housing and the casing in a horizontal and vertical direction is considered, and their looseness amounts are all 1.0μm.
(2)
The looseness fault in the fan fulcrum and the front compressor fulcrum is considered.
(3)
The output is the casing horizontal and the vertical vibration acceleration response in front of compressor fulcrum.
(4)
The speed range is 10000–30000 rpm.


### 3.3. Critical Speed Analysis

Figure5 shows the amplitude-speed curves of the bearing housing S
1 lateral acceleration, the casing lateral acceleration at the ninth node, and the compressor rotor lateral displacement at the first node under imbalance fault. As can be seen from the figures, the first-order and the second-order critical speeds are 17300 rpm and 18500 rpm, respectively.Amplitude-speed curves of bearing, casing lateral acceleration, and rotor displacement (without looseness).
(a)
The bearing housing lateral acceleration
(b)
The casing lateral acceleration
(c)
The rotor lateral displacement


### 3.4. Autocorrelation Method

Autocorrelation technology can be used in signal filter when the signal is weak and sound energy is large, which is effective to detect weak periodic signal ingredients in the rotor system and extract the interested frequency components and multiple frequency components.Because the aeroengine rotor vibration signal is periodic, autocorrelation time domain monitoring method can improve the ability of monitoring weak periodic signal effectively under the condition of the sinusoidal signal frequency position. The periodic signal which is buried in noise can be extracted and the detected signal amplitude has higher precision. The steps are described as follows.(1)
Intercepting2
N points from stationary random signal x
(
n
), x
2
N
(
n
) can be obtained and then using FFT to it, thus X
2
N
(
k
) can be obtained.
(2)
Dividing the square of the amplitude ofX
2
N
(
k
) by N, thus X
2
N
k
2
/
N can be obtained.
(3)
Using IFFT inX
2
N
k
2
/
N, its correlation function can be obtained.


### 3.5. Characteristics of the Bearing HousingS
1 Acceleration Response, the Casing Acceleration Response, and the Rotor Displacement under the Asymmetric Stiffness Model

Figure6 shows cascade plot of the bearing housing S
1 lateral acceleration response at different speeds; Figure 7 shows cascade plot of the casing lateral acceleration response at the ninth node at different speeds; Figure 8 shows cascade plot of the compressor rotor lateral displacement at the first node at different speeds, according to the asymmetric stiffness equation (11). Figure 6 shows that subharmonic resonance appears under the high speed and superharmonic resonance appears when the speed is higher than the second-order critical speed. Comparing Figures 6–8, we find that the multiple frequency components of the bearing housing S
1 acceleration response and the casing acceleration is obvious than the multiple frequency components of the rotor displacement response.Cascade plot of the bearing housing acceleration response at different speeds.
(a)
Speeds from 10000 rpm to 20000 rpm
(b)
Speeds from 20000 rpm to 30000 rpmCascade plot of the casing acceleration response at different speeds.
(a)
Speeds from 10000 rpm to 20000 rpm
(b)
Speeds from 20000 rpm to 30000 rpmCascade plot of the rotor displacement response at different speeds.
(a)
Speeds from 10000 rpm to 20000 rpm
(b)
Speeds from 20000 rpm to 30000 rpmFigures9, 10, and 11 show time domain waveform and spectrum of casing acceleration response at the speeds of 19500 rpm, 26000 rpm, and 30000 rpm, which are obtained before and after the noise reduction of autocorrelation. As can be seen from the figures, there are many multiple frequency components and dividing frequency components, which are caused by the looseness fault; the waveforms after the noise reduction of autocorrelation have a typical shock characteristic with the waveform of truncated shape, which is longitudinal asymmetrical one. The casing acceleration time domain waveforms could determine the characteristics of the looseness fault.Time domain waveform and frequency spectrum before and after the noise reduction when the rotating speed is 19500 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionThe time domain waveform and frequency spectrum before and after the noise reduction when the rotating speed is 26000 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionThe time domain waveform and frequency spectrum before and after the noise reduction when the rotating speed is 30000 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionFigure12 shows cascade plot of the bearing housing S
1 lateral acceleration response at different speeds, according to the symmetric stiffness equation (12). As can be seen from the figures, there are not many multiple frequency components and the looseness fault characteristics would not be obvious. Figure 13 shows the time domain waveform and the spectrum of casing acceleration response at the speed of 25000 rpm, which is obtained before and after the noise reduction of autocorrelation. As can be seen from the figures, the spectrum of acceleration response appears double frequency, caused by the looseness fault. As can be seen from the figures, the shock characteristics are not obvious, without the waveform of truncated shape, and many multiple frequency components could not be seen from the spectrum diagram.Cascade plot of the casing acceleration response at different speeds.
(a)
Speeds from 10000 rpm to 20000 rpm
(b)
Speeds from 20000 rpm to 30000 rpmThe time domain waveform and frequency spectrum before and after the noise reduction when the rotating speed is 25000 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reduction


### 3.6. The Influence of the Looseness Stiffness between the Bearing HousingS
1 and the Casing

In order to study the influence of looseness stiffness on the whole aeroengine vibration, the compressor looseness stiffness between the bearing housingS
1 and the casing is reduced by two times and increased by two times, respectively, and the characteristics of the casing acceleration signal are analyzed.Figure14 shows cascade plot of the casing lateral acceleration response at different speeds according to the asymmetric stiffness equation (11), when the looseness stiffness k
f
1 between the compressor bearing housing S
1 and the casing is reduced by two times. As can be seen from the figures, there are many multiple frequency components and dividing frequency components, which are caused by the looseness fault. Figure 15 shows time domain waveform and spectrum of casing acceleration response at the speed of 17000 rpm, which is obtained before and after the noise reduction of autocorrelation. As can be seen from the figures, the waveforms after the noise reduction of autocorrelation have a typical shock characteristic, with the waveform of truncated shape, which is longitudinal asymmetrical one. The looseness fault would be judged by the casing acceleration time domain waveforms.Cascade plot of the casing acceleration response at different speeds.
(a)
Speeds from 10000 rpm to 20000 rpm
(b)
Speeds from 20000 rpm to 30000 rpmThe time domain waveform and frequency spectrum before and after the noise reduction when rotating speed is 17000 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionFigure16 shows cascade plot of the casing lateral acceleration response at different speeds according to the asymmetric stiffness equation (11), when looseness stiffness k
f
1 between the compressor bearing housing S
1 and the casing is increased by two times. As can be seen from the figures, there are many multiple frequency components and dividing frequency components, which are caused by the looseness fault. Figure 17 shows the time domain waveform and spectrum of the casing acceleration response at the speed of 21000 rpm, which is obtained before and after the noise reduction of autocorrelation. As can be seen from the figures, the waveforms have the characteristic that is a typical shock characteristic, with the waveform of truncated shape, which is longitudinal asymmetrical one. As can be seen from the influence of support stiffness on the looseness fault characteristic, different support stiffness has little effect on the time domain characteristic of the looseness fault, but there is some difference in spectrum; namely, the higher support stiffness will produce much higher multiple frequency components.Cascade plot of the casing acceleration response at different speeds.
(a)
Speeds from 10000 rpm to 20000 rpm
(b)
Speeds from 20000 rpm to 30000 rpmThe time domain waveform and frequency spectrum before and after the noise reduction when rotating speed is 21000 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reduction


## 3.1. Dynamic Model Parameters

The finite element parameters of the rotor, the casing, and the connection parameters of rotor-casing-bearing system are shown in Tables1, 2, 3, 4, 5, 6, and 7.Table 1
Unit number of rotors and casings.
Fan rotor
Compressor rotor
Turbine rotor
Casing
11
10
11
24Table 2
Parameters of rotors.
Parameters
DiskP
1
DiskP
2
DiskP
3
DiskP
4
DiskP
5
Massm
p (kg)
3.88
1.41
5.17
10.28
10.28
Cross-polar inertiaJ
p
d (kg·m2)
0.03
0.003
0.03
0.05
0.05
Cross-equator inertiaJ
d
d (kg·m2)
0.015
0.0015
0.015
0.025
0.025
Elastic modulusE (Pa)
2.07
2.07
2.07
2.07
2.07
Poisson’s ratioμ
0.3
0.3
0.3
0.3
0.3
Densityρ (kg/m3)
7.8 × 103
7.8 × 103
7.8 × 103
7.8 × 103
7.8 × 103
Proportion damping ratioα
0
5
5
5
5
5
Proportion damping ratioα
1
1.35
1.35
1.35
1.35
1.35Table 3
Parameters of casings.
Wall thickness/mm
Elastic modulusE (Pa)
Densityρ (kg/m3)
Poisson’s ratioμ
Proportiondamping ratio α
0
Proportiondamping ratio α
1
15
2.07 × 1011
7.8 × 103
0.3
5
1.35 × 10−5Table 4
Parameters of ball bearing.
Rolling bearing
Outer raceway radiusR/mm
Inner raceway radiusr/mm
Ball numberN
b
Contact stiffnessC
b/(N/m3/2)
Bearing clearancer
0/um
Outer massm
w (kg)
Bearing carrier massm
b (kg)
B
1
39.5
29
13
12.4 × 109
0
2
10
B
2
39.5
29
13
12.4 × 109
0
2
10
B
3
32
17
14
11.9 × 109
0
2
10
B
4
32
17
14
11.9 × 109
0
2
10Table 5
Support parameters of rotor-casing.
Supports
Node of rotor
Casing (node)
k
t (N/m)
c
t (N·s/m)
k
f (N/m)
c
f (N·s/m)
R
C
1
3
2
1 × 108
2000
1 × 108
1000
R
C
2
1
9
1 × 108
2000
1 × 108
1000
R
C
3
11
16
1 × 108
2000
1 × 108
1000
R
C
4
8
22
1 × 108
2000
1 × 108
1000Table 6
Spring collection parameters of rotor-casing.
Collection
Node of rotor
Casing (node)
k
g
x (N/m)
c
g
x (N·s/m)
k
g
y (N/m)
c
g
y (N·s/m)
R
C
1
6
4
1 × 108
0
1 × 108
0Table 7
Collection parameters of casing-base.
Supports
Node of rotor
Casing (node)
k
t (N/m)
c
t (N·s/m)
k
f (N/m)
C
B
1
8
1 × 109
1 × 105
2000
0
C
B
2
23
1 × 109
1 × 105
2000
0


## 3.2. Calculation Condition

(1)
The looseness fault between the bearing housing and the casing in a horizontal and vertical direction is considered, and their looseness amounts are all 1.0μm.
(2)
The looseness fault in the fan fulcrum and the front compressor fulcrum is considered.
(3)
The output is the casing horizontal and the vertical vibration acceleration response in front of compressor fulcrum.
(4)
The speed range is 10000–30000 rpm.


## 3.3. Critical Speed Analysis

Figure5 shows the amplitude-speed curves of the bearing housing S
1 lateral acceleration, the casing lateral acceleration at the ninth node, and the compressor rotor lateral displacement at the first node under imbalance fault. As can be seen from the figures, the first-order and the second-order critical speeds are 17300 rpm and 18500 rpm, respectively.Amplitude-speed curves of bearing, casing lateral acceleration, and rotor displacement (without looseness).
(a)
The bearing housing lateral acceleration
(b)
The casing lateral acceleration
(c)
The rotor lateral displacement


## 3.4. Autocorrelation Method

Autocorrelation technology can be used in signal filter when the signal is weak and sound energy is large, which is effective to detect weak periodic signal ingredients in the rotor system and extract the interested frequency components and multiple frequency components.Because the aeroengine rotor vibration signal is periodic, autocorrelation time domain monitoring method can improve the ability of monitoring weak periodic signal effectively under the condition of the sinusoidal signal frequency position. The periodic signal which is buried in noise can be extracted and the detected signal amplitude has higher precision. The steps are described as follows.(1)
Intercepting2
N points from stationary random signal x
(
n
), x
2
N
(
n
) can be obtained and then using FFT to it, thus X
2
N
(
k
) can be obtained.
(2)
Dividing the square of the amplitude ofX
2
N
(
k
) by N, thus X
2
N
k
2
/
N can be obtained.
(3)
Using IFFT inX
2
N
k
2
/
N, its correlation function can be obtained.


## 3.5. Characteristics of the Bearing HousingS
1 Acceleration Response, the Casing Acceleration Response, and the Rotor Displacement under the Asymmetric Stiffness Model

Figure6 shows cascade plot of the bearing housing S
1 lateral acceleration response at different speeds; Figure 7 shows cascade plot of the casing lateral acceleration response at the ninth node at different speeds; Figure 8 shows cascade plot of the compressor rotor lateral displacement at the first node at different speeds, according to the asymmetric stiffness equation (11). Figure 6 shows that subharmonic resonance appears under the high speed and superharmonic resonance appears when the speed is higher than the second-order critical speed. Comparing Figures 6–8, we find that the multiple frequency components of the bearing housing S
1 acceleration response and the casing acceleration is obvious than the multiple frequency components of the rotor displacement response.Cascade plot of the bearing housing acceleration response at different speeds.
(a)
Speeds from 10000 rpm to 20000 rpm
(b)
Speeds from 20000 rpm to 30000 rpmCascade plot of the casing acceleration response at different speeds.
(a)
Speeds from 10000 rpm to 20000 rpm
(b)
Speeds from 20000 rpm to 30000 rpmCascade plot of the rotor displacement response at different speeds.
(a)
Speeds from 10000 rpm to 20000 rpm
(b)
Speeds from 20000 rpm to 30000 rpmFigures9, 10, and 11 show time domain waveform and spectrum of casing acceleration response at the speeds of 19500 rpm, 26000 rpm, and 30000 rpm, which are obtained before and after the noise reduction of autocorrelation. As can be seen from the figures, there are many multiple frequency components and dividing frequency components, which are caused by the looseness fault; the waveforms after the noise reduction of autocorrelation have a typical shock characteristic with the waveform of truncated shape, which is longitudinal asymmetrical one. The casing acceleration time domain waveforms could determine the characteristics of the looseness fault.Time domain waveform and frequency spectrum before and after the noise reduction when the rotating speed is 19500 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionThe time domain waveform and frequency spectrum before and after the noise reduction when the rotating speed is 26000 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionThe time domain waveform and frequency spectrum before and after the noise reduction when the rotating speed is 30000 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionFigure12 shows cascade plot of the bearing housing S
1 lateral acceleration response at different speeds, according to the symmetric stiffness equation (12). As can be seen from the figures, there are not many multiple frequency components and the looseness fault characteristics would not be obvious. Figure 13 shows the time domain waveform and the spectrum of casing acceleration response at the speed of 25000 rpm, which is obtained before and after the noise reduction of autocorrelation. As can be seen from the figures, the spectrum of acceleration response appears double frequency, caused by the looseness fault. As can be seen from the figures, the shock characteristics are not obvious, without the waveform of truncated shape, and many multiple frequency components could not be seen from the spectrum diagram.Cascade plot of the casing acceleration response at different speeds.
(a)
Speeds from 10000 rpm to 20000 rpm
(b)
Speeds from 20000 rpm to 30000 rpmThe time domain waveform and frequency spectrum before and after the noise reduction when the rotating speed is 25000 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reduction


## 3.6. The Influence of the Looseness Stiffness between the Bearing HousingS
1 and the Casing

In order to study the influence of looseness stiffness on the whole aeroengine vibration, the compressor looseness stiffness between the bearing housingS
1 and the casing is reduced by two times and increased by two times, respectively, and the characteristics of the casing acceleration signal are analyzed.Figure14 shows cascade plot of the casing lateral acceleration response at different speeds according to the asymmetric stiffness equation (11), when the looseness stiffness k
f
1 between the compressor bearing housing S
1 and the casing is reduced by two times. As can be seen from the figures, there are many multiple frequency components and dividing frequency components, which are caused by the looseness fault. Figure 15 shows time domain waveform and spectrum of casing acceleration response at the speed of 17000 rpm, which is obtained before and after the noise reduction of autocorrelation. As can be seen from the figures, the waveforms after the noise reduction of autocorrelation have a typical shock characteristic, with the waveform of truncated shape, which is longitudinal asymmetrical one. The looseness fault would be judged by the casing acceleration time domain waveforms.Cascade plot of the casing acceleration response at different speeds.
(a)
Speeds from 10000 rpm to 20000 rpm
(b)
Speeds from 20000 rpm to 30000 rpmThe time domain waveform and frequency spectrum before and after the noise reduction when rotating speed is 17000 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionFigure16 shows cascade plot of the casing lateral acceleration response at different speeds according to the asymmetric stiffness equation (11), when looseness stiffness k
f
1 between the compressor bearing housing S
1 and the casing is increased by two times. As can be seen from the figures, there are many multiple frequency components and dividing frequency components, which are caused by the looseness fault. Figure 17 shows the time domain waveform and spectrum of the casing acceleration response at the speed of 21000 rpm, which is obtained before and after the noise reduction of autocorrelation. As can be seen from the figures, the waveforms have the characteristic that is a typical shock characteristic, with the waveform of truncated shape, which is longitudinal asymmetrical one. As can be seen from the influence of support stiffness on the looseness fault characteristic, different support stiffness has little effect on the time domain characteristic of the looseness fault, but there is some difference in spectrum; namely, the higher support stiffness will produce much higher multiple frequency components.Cascade plot of the casing acceleration response at different speeds.
(a)
Speeds from 10000 rpm to 20000 rpm
(b)
Speeds from 20000 rpm to 30000 rpmThe time domain waveform and frequency spectrum before and after the noise reduction when rotating speed is 21000 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reduction


## 4. Verification and Analysis of the Real Trial Running Data

Figure18 shows cascade plot of the real trial running data. The measuring points are the compressor casing horizontal and vertical direction of the compressor fulcrum. As can be seen from the figures, there are many multiple frequency components and dividing frequency components. Figures 19, 20, and 21 show time domain waveform and spectrum of the casing acceleration response at the speeds of 24000 rpm, 27000 rpm, and 29000 rpm, respectively, which is obtained before and after the noise reduction of autocorrelation. As can be seen from the figures, the time domain waveform characteristics before the noise reduction of autocorrelation are not obvious; however, the waveforms after the noise reduction of autocorrelation have a typical shock characteristic with the waveform of truncated shape, which is longitudinal asymmetrical one. The real trial running data and the decomposition of the engine show that the turbofan engine presents the bearing outer ring obvious signs of wear and uneven wear phenomena, which proves the existence of support looseness. Based on more fault simulation analysis, only the asymmetric stiffness of the looseness fault excites the casing acceleration signal, which is very similar to the real trial running data. Not only are the multiple frequency characteristics similar in the frequency domain, but also the waveforms after the noise reduction of autocorrelation have a typical shock characteristic, with the waveform of truncated shape, which is longitudinal asymmetrical one. To some extent, the source of the engine fault is determined by the support looseness, which indicates that the looseness fault would be a main fault.Cascade plot of casing acceleration response at different measuring points.
(a)
Horizontal measuring point
(b)
Vertical measuring pointThe time domain waveform and frequency spectrum before and after the noise reduction when rotating speed is 24000 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionThe time domain waveform and frequency spectrum before and after the noise reduction when rotating speed is 27000 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionThe time domain waveform and frequency spectrum before and after the noise reduction when the rotating speed is 29000 rpm.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reduction


## 5. Connector Looseness Fault Experimental Verification

### 5.1. Connector Looseness Experiment Principle

In order to verify the nonsynchronous response characteristics of connector looseness fault, the experiment rig with looseness clearance is established, and the connector looseness experiment is carried out. Figure22 is the experimental site map of looseness fault, Figure 23 is the experimental schematic diagram of looseness fault, and Figure 24 is the experimental three-dimensional diagram of looseness fault.Figure 22
The experimental site map of looseness fault.Figure 23
Experimental schematic diagram of looseness fault.Figure 24
Experimental three-dimensional diagram of looseness fault.Figure24 shows that the second disk is connected to the third disk by spring k
1, the third disk is connected to vibrating table by bolts, and the first disk is connected to the third disk by using three polish rods. Then the first disk, the third disk, and vibrating table will vibrate together. The first disk and the second disk are connected by spring k
2, and the connector looseness fault is simulated by controlling the clearance.During the experiment, the first disk, the third disk, and vibration table vibrate together, and the basement excitation to the second disk is carried out. When the relative displacement between the second disk and the third disk is smaller, the springk
2 does not connect with the second disk, so the spring connection stiffness is only k
1. When the relative displacement between the second disk and the third disk is larger, the spring k
2 connects with the second disk, so the spring connection stiffness is k
1
+
k
2.The energy of vibrating table is input by power amplifier and the feedback of vibration signals is realized by controlling the software and the vibration acceleration sensor on the third disk. The vibration of the vibrating table is controlled at the designated frequency and amplitude by vibration controller. The vibration acceleration of the second disk is measured by the vibration acceleration sensor on the second disk. The collected digital signals are input into the computer for preservation by NI USB-9234 capture card.In order to obtain the vibration performance of the vibration system, the excitation on vibrating table is applied by using the linear sine sweep frequency method from 5 Hz to 100 Hz, and the looseness fault characteristics at different frequencies are obtained.


### 5.2. Looseness Fault Experiment Verification

#### 5.2.1. Analysis of the Second Disk Response at Different Contact Forms

In the experiment, the acceleration of the second disk is measured by three methods; namely, (1) when the second disk does not contact with the springk
2, the second disk only contacts with the third disk through spring k
1, the vibration system is a linear system, and the experiment result is shown in Figure 25(a); (2) when the second disk always contacts with the spring k
2, the second disk contacts with the third disk and the first disk through spring k
1 and k
2, the vibration system is also a linear system, its spring connection stiffness is k
1
+
k
2, and the experiment result is shown in Figure 25(b); (3) when there is clearance between the second disk and spring k
1, at different relative vibration displacement between the second disk and the first disk, the second disk may contact with spring k
1 ornot, the vibration system is a nonlinear system, and the experiment result is shown in Figure 25(c).The response of the second disk under different contact conditions.
(a)
No contact between the second disk and springk
2
(b)
Complete contact between the second disk and springk
2
(c)
Discontinuous contact between the second disk and springk
2Figure25(a) shows that when the second disk does not contact with the spring k
2, the vibration system is a linear system, and its natural frequency is 14.6 Hz. Figure 25(b) shows that when the second disk always contacts with the spring k
2, the vibration system is also a linear system, and its spring connection stiffness is k
1
+
k
2, so its natural frequency is 51.22 Hz. When there is looseness clearance, the second disk may contact with spring k
1 or not. When the vibration displacement is large, the second disk will contact spring k
2 only half period in a movement period. At this time, the natural frequency of the system is [18]
(14)
1
f
n
=
1
2
1
f
1
+
1
2
1
f
2
.According to (14), the calculated result f
n is 22.7 Hz, and the experimental result is 22.2 Hz, shown in Figure 25(c). Obviously, the analysis results are verified fully.


#### 5.2.2. Analysis of Looseness Characteristics

Figure26 is three-dimensional cascade plot of the second disk vibration acceleration changing with excitation frequency and response frequency. Figure 26 shows that the frequency multiplication appears when excitation frequency is 1
/
2
f
n, 1/2 frequency division appears when excitation frequency is 2
f
n, and 1/3 frequency division appears when excitation frequency is 3f
n.Figure 26
Cascade plot of the second disk response.Figure27 shows time domain waveform and spectrum of the second disk when the excitation frequency is 12.5 Hz, which is obtained before and after the noise reduction of autocorrelation. As shown in Figure 27, when the excitation frequency is equal to 1/2f
n, the larger double frequency components appeared, that is, the natural frequency of the system; the waveforms after the noise reduction of autocorrelation have a typical shock characteristic with the waveform of truncated shape, which is longitudinal asymmetrical one.Waveform and frequency spectrum before and after the noise reduction when frequency is1
/
2
f
n.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionFigure28 shows time domain waveform and spectrum of the second disk when the excitation frequency is 45 Hz, which is obtained before and after the noise reduction of autocorrelation. As shown in Figure 28, when the frequency of the vibration table is equal to 2f
n, the larger 1/2 dividing frequency components appeared, that is, the natural frequency of the system; the waveforms after the noise reduction of autocorrelation have a typical shock characteristic with the waveform of truncated shape, which is longitudinal asymmetrical one.Waveform and frequency spectrum before and after the noise reduction when frequency is2
f
n.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionFigure29 shows time domain waveform and spectrum of the second disk when the excitation frequency is 80 Hz, which is obtained before and after the noise reduction of autocorrelation. As shown in Figure 29, when the frequency of the vibrating table is about 3f
n, the larger 1/3 dividing frequency components appeared, that is, the natural frequency of the system; the waveforms after the noise reduction of autocorrelation have a typical shock characteristic, with the waveform of truncated shape, which is longitudinal asymmetrical one.Waveform and frequency spectrum before and after the noise reduction when frequency is3
f
n.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionBy comparison, in different frequencies of the vibrating table, for connector looseness fault, the acceleration waveforms after the noise reduction of autocorrelation have the characteristic that is a typical shock characteristic, with the waveform of truncated shape, which is longitudinal asymmetrical one, which could be identified as the waveforms after the noise reduction of autocorrelation have a typical shock characteristic, of the looseness fault.


## 5.1. Connector Looseness Experiment Principle

In order to verify the nonsynchronous response characteristics of connector looseness fault, the experiment rig with looseness clearance is established, and the connector looseness experiment is carried out. Figure22 is the experimental site map of looseness fault, Figure 23 is the experimental schematic diagram of looseness fault, and Figure 24 is the experimental three-dimensional diagram of looseness fault.Figure 22
The experimental site map of looseness fault.Figure 23
Experimental schematic diagram of looseness fault.Figure 24
Experimental three-dimensional diagram of looseness fault.Figure24 shows that the second disk is connected to the third disk by spring k
1, the third disk is connected to vibrating table by bolts, and the first disk is connected to the third disk by using three polish rods. Then the first disk, the third disk, and vibrating table will vibrate together. The first disk and the second disk are connected by spring k
2, and the connector looseness fault is simulated by controlling the clearance.During the experiment, the first disk, the third disk, and vibration table vibrate together, and the basement excitation to the second disk is carried out. When the relative displacement between the second disk and the third disk is smaller, the springk
2 does not connect with the second disk, so the spring connection stiffness is only k
1. When the relative displacement between the second disk and the third disk is larger, the spring k
2 connects with the second disk, so the spring connection stiffness is k
1
+
k
2.The energy of vibrating table is input by power amplifier and the feedback of vibration signals is realized by controlling the software and the vibration acceleration sensor on the third disk. The vibration of the vibrating table is controlled at the designated frequency and amplitude by vibration controller. The vibration acceleration of the second disk is measured by the vibration acceleration sensor on the second disk. The collected digital signals are input into the computer for preservation by NI USB-9234 capture card.In order to obtain the vibration performance of the vibration system, the excitation on vibrating table is applied by using the linear sine sweep frequency method from 5 Hz to 100 Hz, and the looseness fault characteristics at different frequencies are obtained.


## 5.2. Looseness Fault Experiment Verification

### 5.2.1. Analysis of the Second Disk Response at Different Contact Forms

In the experiment, the acceleration of the second disk is measured by three methods; namely, (1) when the second disk does not contact with the springk
2, the second disk only contacts with the third disk through spring k
1, the vibration system is a linear system, and the experiment result is shown in Figure 25(a); (2) when the second disk always contacts with the spring k
2, the second disk contacts with the third disk and the first disk through spring k
1 and k
2, the vibration system is also a linear system, its spring connection stiffness is k
1
+
k
2, and the experiment result is shown in Figure 25(b); (3) when there is clearance between the second disk and spring k
1, at different relative vibration displacement between the second disk and the first disk, the second disk may contact with spring k
1 ornot, the vibration system is a nonlinear system, and the experiment result is shown in Figure 25(c).The response of the second disk under different contact conditions.
(a)
No contact between the second disk and springk
2
(b)
Complete contact between the second disk and springk
2
(c)
Discontinuous contact between the second disk and springk
2Figure25(a) shows that when the second disk does not contact with the spring k
2, the vibration system is a linear system, and its natural frequency is 14.6 Hz. Figure 25(b) shows that when the second disk always contacts with the spring k
2, the vibration system is also a linear system, and its spring connection stiffness is k
1
+
k
2, so its natural frequency is 51.22 Hz. When there is looseness clearance, the second disk may contact with spring k
1 or not. When the vibration displacement is large, the second disk will contact spring k
2 only half period in a movement period. At this time, the natural frequency of the system is [18]
(14)
1
f
n
=
1
2
1
f
1
+
1
2
1
f
2
.According to (14), the calculated result f
n is 22.7 Hz, and the experimental result is 22.2 Hz, shown in Figure 25(c). Obviously, the analysis results are verified fully.


### 5.2.2. Analysis of Looseness Characteristics

Figure26 is three-dimensional cascade plot of the second disk vibration acceleration changing with excitation frequency and response frequency. Figure 26 shows that the frequency multiplication appears when excitation frequency is 1
/
2
f
n, 1/2 frequency division appears when excitation frequency is 2
f
n, and 1/3 frequency division appears when excitation frequency is 3f
n.Figure 26
Cascade plot of the second disk response.Figure27 shows time domain waveform and spectrum of the second disk when the excitation frequency is 12.5 Hz, which is obtained before and after the noise reduction of autocorrelation. As shown in Figure 27, when the excitation frequency is equal to 1/2f
n, the larger double frequency components appeared, that is, the natural frequency of the system; the waveforms after the noise reduction of autocorrelation have a typical shock characteristic with the waveform of truncated shape, which is longitudinal asymmetrical one.Waveform and frequency spectrum before and after the noise reduction when frequency is1
/
2
f
n.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionFigure28 shows time domain waveform and spectrum of the second disk when the excitation frequency is 45 Hz, which is obtained before and after the noise reduction of autocorrelation. As shown in Figure 28, when the frequency of the vibration table is equal to 2f
n, the larger 1/2 dividing frequency components appeared, that is, the natural frequency of the system; the waveforms after the noise reduction of autocorrelation have a typical shock characteristic with the waveform of truncated shape, which is longitudinal asymmetrical one.Waveform and frequency spectrum before and after the noise reduction when frequency is2
f
n.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionFigure29 shows time domain waveform and spectrum of the second disk when the excitation frequency is 80 Hz, which is obtained before and after the noise reduction of autocorrelation. As shown in Figure 29, when the frequency of the vibrating table is about 3f
n, the larger 1/3 dividing frequency components appeared, that is, the natural frequency of the system; the waveforms after the noise reduction of autocorrelation have a typical shock characteristic, with the waveform of truncated shape, which is longitudinal asymmetrical one.Waveform and frequency spectrum before and after the noise reduction when frequency is3
f
n.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionBy comparison, in different frequencies of the vibrating table, for connector looseness fault, the acceleration waveforms after the noise reduction of autocorrelation have the characteristic that is a typical shock characteristic, with the waveform of truncated shape, which is longitudinal asymmetrical one, which could be identified as the waveforms after the noise reduction of autocorrelation have a typical shock characteristic, of the looseness fault.


## 5.2.1. Analysis of the Second Disk Response at Different Contact Forms

In the experiment, the acceleration of the second disk is measured by three methods; namely, (1) when the second disk does not contact with the springk
2, the second disk only contacts with the third disk through spring k
1, the vibration system is a linear system, and the experiment result is shown in Figure 25(a); (2) when the second disk always contacts with the spring k
2, the second disk contacts with the third disk and the first disk through spring k
1 and k
2, the vibration system is also a linear system, its spring connection stiffness is k
1
+
k
2, and the experiment result is shown in Figure 25(b); (3) when there is clearance between the second disk and spring k
1, at different relative vibration displacement between the second disk and the first disk, the second disk may contact with spring k
1 ornot, the vibration system is a nonlinear system, and the experiment result is shown in Figure 25(c).The response of the second disk under different contact conditions.
(a)
No contact between the second disk and springk
2
(b)
Complete contact between the second disk and springk
2
(c)
Discontinuous contact between the second disk and springk
2Figure25(a) shows that when the second disk does not contact with the spring k
2, the vibration system is a linear system, and its natural frequency is 14.6 Hz. Figure 25(b) shows that when the second disk always contacts with the spring k
2, the vibration system is also a linear system, and its spring connection stiffness is k
1
+
k
2, so its natural frequency is 51.22 Hz. When there is looseness clearance, the second disk may contact with spring k
1 or not. When the vibration displacement is large, the second disk will contact spring k
2 only half period in a movement period. At this time, the natural frequency of the system is [18]
(14)
1
f
n
=
1
2
1
f
1
+
1
2
1
f
2
.According to (14), the calculated result f
n is 22.7 Hz, and the experimental result is 22.2 Hz, shown in Figure 25(c). Obviously, the analysis results are verified fully.


## 5.2.2. Analysis of Looseness Characteristics

Figure26 is three-dimensional cascade plot of the second disk vibration acceleration changing with excitation frequency and response frequency. Figure 26 shows that the frequency multiplication appears when excitation frequency is 1
/
2
f
n, 1/2 frequency division appears when excitation frequency is 2
f
n, and 1/3 frequency division appears when excitation frequency is 3f
n.Figure 26
Cascade plot of the second disk response.Figure27 shows time domain waveform and spectrum of the second disk when the excitation frequency is 12.5 Hz, which is obtained before and after the noise reduction of autocorrelation. As shown in Figure 27, when the excitation frequency is equal to 1/2f
n, the larger double frequency components appeared, that is, the natural frequency of the system; the waveforms after the noise reduction of autocorrelation have a typical shock characteristic with the waveform of truncated shape, which is longitudinal asymmetrical one.Waveform and frequency spectrum before and after the noise reduction when frequency is1
/
2
f
n.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionFigure28 shows time domain waveform and spectrum of the second disk when the excitation frequency is 45 Hz, which is obtained before and after the noise reduction of autocorrelation. As shown in Figure 28, when the frequency of the vibration table is equal to 2f
n, the larger 1/2 dividing frequency components appeared, that is, the natural frequency of the system; the waveforms after the noise reduction of autocorrelation have a typical shock characteristic with the waveform of truncated shape, which is longitudinal asymmetrical one.Waveform and frequency spectrum before and after the noise reduction when frequency is2
f
n.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionFigure29 shows time domain waveform and spectrum of the second disk when the excitation frequency is 80 Hz, which is obtained before and after the noise reduction of autocorrelation. As shown in Figure 29, when the frequency of the vibrating table is about 3f
n, the larger 1/3 dividing frequency components appeared, that is, the natural frequency of the system; the waveforms after the noise reduction of autocorrelation have a typical shock characteristic, with the waveform of truncated shape, which is longitudinal asymmetrical one.Waveform and frequency spectrum before and after the noise reduction when frequency is3
f
n.
(a)
The waveform before the noise reduction
(b)
The waveform after the noise reduction
(c)
The spectrum before the noise reduction
(d)
The spectrum after the noise reductionBy comparison, in different frequencies of the vibrating table, for connector looseness fault, the acceleration waveforms after the noise reduction of autocorrelation have the characteristic that is a typical shock characteristic, with the waveform of truncated shape, which is longitudinal asymmetrical one, which could be identified as the waveforms after the noise reduction of autocorrelation have a typical shock characteristic, of the looseness fault.


## 6. Conclusion

In this paper, certain type turbofan engine whole vibration model with support looseness fault is established and casing acceleration response is obtained. Some results are obtained as follows.(1)
Aiming at certain type turbofan engine, the dynamic model with the looseness fault is established. The nonlinear rolling bearing is considered in the model. A variety of support connections are used, that is, the support connection between the casing and the rotor; the coupling connection between two rotors; and the elastic connection between the rotor and the casing.
(2)
The looseness models of asymmetric stiffness and symmetric stiffness are established. The looseness fault exists between the bearing and the casing. Numerical simulation method is used to obtain the engine casing acceleration response under the looseness fault.
(3)
The simulation results are compared with the real trial running vibration signals and connector looseness fault experiment, and the simulated fault characteristics are very similar with the real trial running vibration characteristics and the connectors looseness characteristics. The waveforms have the characteristic that is a typical shock characteristic, with the waveform of truncated shape, which is longitudinal asymmetrical one, and there are many multiple frequency components and dividing frequency components, which could provide a qualitative judgment of the support looseness fault for the engine fault diagnosis. In addition, the comparison results of the simulation and the real trial running vibration signals also show the validity of the asymmetric stiffness looseness model.



---

*Source: 683469-2014-12-17.xml*