# Comparison of Common Methods in Dynamic Response Predictions of Rotor Systems with Malfunctions

Authors: Hongliang Yao; Qian Zhao; Qi Xu; Bangchun Wen
Journal: International Journal of Rotating Machinery
 (2014)

---

## Abstract

The efficiency and accuracy of common time and frequency domain methods that are used to simulate the response of a rotor system with malfunctions are compared and analyzed. The Newmark method and the incremental harmonic balance method are selected as typical representatives of time and frequency domain methods, respectively. To improve the simulation efficiency, the fixed interface component mode synthesis approach is combined with the Newmark method and the receptance approach is combined with the incremental harmonic balance method. Numerical simulations are performed for rotor systems with single and double frequency excitations. The inherent characteristic that determines the efficiency of the two methods is analyzed. The results of the analysis indicated that frequency domain methods are suitable single and double frequency excitation rotor systems, whereas time domain methods are more suitable for multifrequency excitation rotor systems.


---

## Body

## 1. Introduction

Numerical simulation plays an important role in the study of rotor systems with malfunctions. With complex structures of the modern rotor systems (e.g., multidisc or parallel shafts) and the inherent strong nonlinear characteristics of malfunctions [1] (e.g., the self-excitation property of oil whirl and oil whip, the parametric vibration characteristics of the crack rotor, and the nonlinear stiffness of rotor to stator rub), few analytical solutions work well in analyzing these strongly nonlinear systems with large degree of freedom (DOF). Therefore, numerical simulation may be the only way to predict the response of these systems accurately.The numerical simulation methods commonly used in rotor dynamics can be divided into two types: time domain methods and frequency domain methods. The time domain methods mainly include the Runge-Kutta method [2], Newmark-β method [3], and the shooting method [4]. The frequency domain methods mainly include the harmonic balance (HB) method [5], describing function (DF) method [6], incremental harmonic balance (IHB) method [7], and so forth.When analyzing large-scale rotor systems with complex structures, both the time domain and the frequency domain methods will encounter the challenge of solving a large system of equations with numerous interdependent variables. Therefore, a dimension reduction approach must be adopted. There are two kinds of dimension reduction approaches that are widely used: the dynamical substructure approach and the receptance-based approach. The former is often combined with time domain methods and the latter is often combined with frequency domain methods. The fixed interface component mode synthesis approach [8] (CMS) is the most widely used dynamic substructure approach, in which the structure is divided into the linear part and the nonlinear part. Modal truncation is applied to the linear portion and only the lower-order modes are retained, which reduces the dimension of the former structure. The frequency domain methods often achieve dimension reduction using linear receptance data [9, 10], in which the vibrations of the linear parts are substituted by those of the nonlinear parts. The dimension of the former system then reduces to equal the number of nonlinear DOFs.Both the types of methods have advantages and disadvantages. Frequency domain methods are highly efficient because they can skip the transient response and obtain the steady state response directly. However, these methods can only seek periodic and quasiperiodic responses. They also generally need prior information on the behavior of the system, such as which harmonic terms in the responses are dominant [11]. Time integration methods are capable of seeking periodic, aperiodic, and quasiperiodic solutions as well as transient responses. However, they have low efficiency when solving for steady state responses since the transient responses must be solved first.Modern rotor systems have various forms and excitation conditions. Therefore, selecting an appropriate numerical method with good precision and high efficiency is very important when analyzing a rotor system with fault. Few studies have compared the efficiency and accuracy of the time and frequency domain methods or discussed the application scope of these methods when applied to malfunctioning rotor systems. In this study, the two types of methods are compared and analyzed. Typical time domain and frequency domain methods are applied and compared under different conditions of excitation. The inherent characteristic that determines the efficiency of the two methods is analyzed. The scopes of application of the two kinds of methods are discussed.


## 2. Typical Time and Frequency Domain Methods

### 2.1. The Newmark-β Method and the Fixed Interface Component Mode Synthesis Approach

The Newmark-β method is a commonly used time integration method with good convergence and computational stability. Therefore it is chosen here to represent time domain methods.The equation of motion of a malfunctioning rotor system can be written as(1)
M
x
¨
+
D
x
˙
+
K
x
+
f
n
x
,
x
˙
=
F
t
,
where M, D, and K are the mass, damping (including linear viscous damping and gyroscopic moment), and stiffness matrices, respectively; x is the displacement vector; F
t is the excitation vector caused by the imbalance; f
n
(
x
,
x
˙
) is the nonlinear force vector, which is a function of x and x
˙; the dots above x in (1) denote derivatives with respect to time t.LettingN
x
,
x
˙
=
D
x
˙
+
K
x
+
f
n
x
,
x
˙, (1) can be changed to
(2)
M
x
¨
+
N
x
,
x
˙
=
F
t
.Based on known values ofx
n, x
˙
n, x
¨
n at time t
n, the approximate values x
¨
n
+
1 and x
˙
n
+
1 are [12]
(3)
x
˙
n
+
1
=
a
1
x
n
+
1
-
x
n
-
a
4
x
˙
n
-
a
5
x
¨
n
,
x
¨
n
+
1
=
a
0
x
n
+
1
-
x
n
-
a
2
x
˙
n
-
a
3
x
¨
n
,
where a
0
=
1
/
β
Δ
t
2, a
1
=
γ
/
β
Δ
t, a
2
=
1
/
β
Δ
t, a
3
=
1
/
2
β
-
1, a
4
=
γ
/
β
-
1, and a
5
=
Δ
t
/
2
γ
/
β
-
2.Let(4)
x
˙
~
n
+
1
=
-
a
1
x
n
-
a
4
x
˙
n
-
a
5
x
¨
n
,
x
¨
~
n
+
1
=
-
a
0
x
n
-
a
2
x
˙
n
-
a
3
x
¨
n
.Equation (3) can then be rewritten as
(5)
x
˙
n
+
1
=
x
˙
~
n
+
1
+
a
1
x
n
+
1
,
x
¨
n
+
1
=
x
¨
~
n
+
1
+
a
0
x
n
+
1
.Assuming thatx
n
+
1
i
+
1
=
x
n
+
1
i
+
Δ
x
n
+
1
i, where Δ
x
n
+
1
i is the increment, the above equation is given by
(6)
x
˙
n
+
1
i
+
1
=
x
˙
n
+
1
i
+
a
1
Δ
x
n
+
1
i
,
x
¨
n
+
1
i
+
1
=
x
¨
n
+
1
i
+
a
0
Δ
x
n
+
1
i
.Expanding (2) via a Taylor series around x
n
+
1
i yields
(7)
Δ
x
n
+
1
i
a
0
M
+
a
1
C
(
x
˙
n
+
1
i
)
+
K
x
n
+
1
i
=
F
n
+
1
t
-
M
x
¨
n
+
1
i
-
N
x
˙
n
+
1
i
,
x
n
+
1
i
.WithJ
n
+
1
=
a
0
M
+
a
1
C
(
x
˙
n
+
1
i
)
+
K
(
x
n
+
1
i
), (7) can be changed to
(8)
J
n
+
1
Δ
x
n
+
1
i
=
F
n
+
1
t
-
M
x
¨
n
+
1
i
-
N
x
˙
n
+
1
i
,
x
n
+
1
i
,
where J
n
+
1 is the instantaneous stiffness matrix. Δ
x
n
+
1
i can be obtained using (8) and x
n
+
1
i
+
1
=
x
n
+
1
i
+
Δ
x
n
+
1
i, x
n
+
1 can be obtained when convergence criteria are satisfied.There are four DOFs in one node and if the rotor system is made up ofN nodes, there are a total of 4
N DOFs in the system. The dimension of the instantaneous stiffness matrix J
n
+
1 is 4
N
×
4
N. As N increases, the calculation speed decreases rapidly. The fixed interface component mode synthesis approach can be used to reduce the number of dimensions of the system and thus the simulation efficiency can be increased.The DOFs of the rotor system are separated into two distinct groups. The linear groupx
1 contains the DOFs that are not related to the nonlinear forces. The nonlinear group x
2 contains the DOFs that are related to the nonlinear forces. Equation (2) is rearranged into the following form:
(9)
M
i
i
M
i
j
M
j
i
M
j
j
u
¨
i
u
¨
j
+
C
i
i
C
i
j
C
j
i
C
j
j
u
˙
i
u
˙
j
+
K
i
i
K
i
j
K
j
i
K
j
j
u
i
u
j
=
F
i
F
j
+
0
F
b
.The nonlinear DOFs are fixed and the modes of the system comprised of the remaining DOFs are calculated. The modes that contribute most to the responses of the former rotor system are chosen to construct the fixed interface normal mode set:(10)
ψ
m
=
ψ
l
0
T
.Assuming unit displacement in each nonlinear DOF and setting the movement of the linear DOFs to zero, the constrained mode of the nonlinear DOFs can be obtained:(11)
X
i
=
-
K
i
i
-
1
K
i
j
X
j
.The constrained mode set can then be constructed from all the constrained modes:(12)
ψ
c
=
-
K
i
i
-
1
K
i
j
I
.The modal matrix of the dimension-reduced system can then be obtained by combining the normal mode set and the constrained mode set:(13)
ψ
=
ψ
m
ψ
c
.Using the modal matrix, (9) can be rewritten as:
(14)
M
~
u
¨
+
D
~
u
˙
+
K
~
u
=
F
~
t
+
0
F
b
,
where M
~
=
ψ
T
M
i
i
M
i
j
M
j
i
M
j
j
ψ, D
~
=
ψ
T
D
i
i
D
i
j
D
j
i
D
j
j
ψ, K
~
=
ψ
T
K
i
i
K
i
j
K
j
i
K
j
j
ψ, and F
~
=
ψ
T
F
i
F
j.The responseu can be obtained using Newmark method, and the response x can then be obtained by x
=
ψ
u.In the derivation process above, the dimension reduction mainly depends on the reduction of the normal modes. However, this may result in a loss of accuracy.


### 2.2. Incremental Harmonic Balance Method and Dimension Reduction Using Receptance Data

In the three commonly used frequency domain simulation methods (HB, DF, and IHB), the IHB method is equivalent to the HB method plus the Newton-Raphson method [13]. The IHB method can handle strong nonlinearities and is therefore chosen to represent the frequency domain methods. Here, the multidimensional IHB is introduced because there are many multifrequency excited systems in engineering, such as the dual-rotor system in aeroengines. The detailed information of the single frequency IHB method can be seen in [14, 15].The equation of motion of a multifrequency excited rotor system is expressed as(15)
M
x
¨
+
C
x
˙
+
K
x
+
f
n
x
,
x
˙
=
∑
i
=
1
r
F
c
i
cos
⁡
ω
i
t
+
F
s
i
sin
ω
i
t
,
where ω
i are the exciting frequencies and F
c
i and F
s
i are the amplitudes of the exciting forces corresponding to ω
i.When the rotor system is in steady state, the responsex of the rotor system can be expanded to the form of a multiple Fourier series:
(16)
x
k
=
a
0
k
+
∑
l
1
=
-
R
1
R
1
⋯
∑
l
r
=
-
R
r
R
r
a
l
1
⋯
l
r
k
cos
⁡
∑
i
=
1
r
l
i
ω
i
t
+
b
l
1
⋯
l
r
k
sin
∑
i
=
1
r
l
i
ω
i
t
.Assumingω
=
ω
1
,
…
,
ω
r
T, l
i
=
l
1
,
…
,
l
r, and ν
i
=
l
i
·
ω, (10) can be simplified to
(17)
x
k
=
a
0
k
+
∑
i
=
1
p
c
i
k
cos
⁡
ν
i
t
+
d
i
k
sin
ν
i
t
=
C
s
A
k
,
where C
s
=
1
cos
⁡
ν
1
t
sin
ν
1
t
⋯
cos
⁡
ν
p
t
sin
ν
p
t and A
k
=
a
0
k
c
1
k
d
1
k
⋯
c
p
k
d
p
k
T.LetS
=
diag
⁡
C
s
C
s
⋯
C
s and A
=
A
1
T
A
2
T
⋯
A
N
T
T. Therefore, x
=
S
A. Assume that
(18)
x
=
x
0
+
ε
x
.Substitute (18) into (15) and neglecting the high-order harmonic components, one obtains
(19)
M
S
¨
ε
A
+
C
S
˙
ε
A
+
K
S
ε
A
+
C
n
S
˙
ε
A
+
K
n
S
ε
A
=
-
M
S
¨
A
0
+
C
S
˙
A
0
+
K
S
A
0
+
∑
i
=
1
r
F
c
i
cos
⁡
ω
i
t
+
F
s
i
sin
ω
i
t
-
f
n
x
0
,
where K
n
=
∂
f
n
/
∂
x and C
n
=
∂
f
n
/
∂
x
˙.Letτ
=
ω
1
t
,
…
,
ω
r
t
T
=
τ
1
,
…
,
τ
r
T and applying the Galerkin’s procedure, (19) yields
(20)
∫
0
2
π
⋯
∫
0
2
π
S
T
M
S
¨
+
C
S
˙
+
K
S
+
C
n
S
˙
+
K
n
S
d
τ
1
⋯
d
τ
r
·
ε
A
=
-
∫
0
2
π
⋯
∫
0
2
π
S
T
M
S
¨
+
C
S
˙
+
K
S
d
τ
1
⋯
d
τ
r
·
A
0
+
∫
0
2
π
⋯
∫
0
2
π
S
T
∑
i
=
1
r
F
c
i
cos
⁡
ω
i
t
+
F
s
i
sin
ω
i
t
-
f
n
x
0
∑
i
=
1
r
F
c
i
cos
⁡
ω
i
t
+
F
s
i
sin
ω
i
t
r
d
τ
1
⋯
d
τ
.Equation (20) can be rewritten as
(21)
K
m
Δ
A
=
R
m
1
A
0
+
R
m
2
.A
0 can be obtained from (21) iteratively, following which x can be obtained.K
m is the instantaneous stiffness matrix and has dimensions N
×
2
r
+
1
×
N
×
2
r
+
1. An increase in the DOFs and the number of harmonic components will greatly reduce the computational efficiency of the IHB method.The receptance data is used in the dimension reduction strategy in the IHB method. The DOFs of the rotor system are separated into two distinct groups. The nonlinear groupx
1 contains the DOFs that are related to the nonlinear forces and the linear group x
2 contains the DOFs that are not related to the nonlinear forces. Rearrange (15) to
(22)
M
11
M
12
M
21
M
22
x
¨
1
x
¨
2
+
C
11
C
12
C
21
C
22
x
˙
1
x
˙
2
+
K
11
K
12
K
21
K
22
x
1
x
2
+
f
n
1
x
1
,
x
˙
1
0
=
F
11
F
21
cos
⁡
ω
t
+
F
12
F
22
sin
ω
t
.When the nonlinear system is in steady state, the response vectorx is periodic and the nonlinear force vector f
n
x
,
x
˙ is also periodic, which can be expressed as
(23)
f
n
x
,
x
˙
=
S
P
,
where P
=
P
1
T
P
2
T
⋯
P
h
T
T and P
i
=
p
0
i
p
1
i
q
1
i
⋯
p
r
i
q
r
i
T.Assume thatC
s
e
=
1
e
j
ω
t
j
e
-
j
ω
t
⋯
e
j
r
ω
t
j
e
-
j
r
ω
t, S
e
=
diag
⁡
C
s
e
C
s
e
⋯
C
s
e, and x
-
=
S
e
A. Therefore x
=
Re
x
- and x
- satisfies the following equation:
(24)
M
x
-
¨
+
C
x
-
˙
+
K
x
-
+
f
-
n
=
F
1
e
j
ω
t
+
j
F
2
e
-
j
ω
t
,
where f
-
n
=
S
e
P and f
n
=
Re
f
-
n.Rewrite (24) as
(25)
M
11
M
12
M
21
M
22
x
-
¨
1
x
-
¨
2
+
C
11
C
12
C
21
C
22
x
-
˙
1
x
-
˙
2
+
K
11
K
12
K
21
K
22
x
-
1
x
-
2
+
f
-
n
1
0
=
F
11
F
21
e
j
ω
t
+
j
F
12
F
22
e
j
ω
t
.Assume thatx
-
=
∑
k
=
0
r
x
-
H
k, where x
-
k
=
x
-
H
1
k
x
-
H
2
k
⋯
x
-
H
N
k is the kth order harmonic term of x
-, and x
-
H
i
k
=
p
i
k
e
j
k
ω
t
+
j
q
i
k
e
-
j
k
ω
t. From (25), the kth harmonic term of x
-
H
2
k can be written as
(26)
x
-
H
2
k
=
-
K
22
-
M
22
k
2
ω
2
+
j
C
22
k
ω
-
1
×
K
21
-
M
21
k
2
ω
2
+
j
C
21
k
ω
x
-
H
1
k
.Substituting (26) into the first equation of (25) yields
(27)
M
e
k
x
-
¨
H
1
k
+
C
e
k
x
-
˙
H
1
k
+
K
e
k
x
-
H
1
k
+
f
-
H
n
1
k
=
0
,
where
(28)
M
e
k
=
M
11
-
M
12
K
22
-
M
22
k
2
ω
2
+
j
C
22
k
ω
-
1
×
K
21
-
M
21
k
2
ω
2
+
j
C
21
k
ω
,
C
e
k
=
C
11
-
C
12
K
22
-
M
22
k
2
ω
2
+
j
C
22
k
ω
-
1
×
K
21
-
M
21
k
2
ω
2
+
j
C
21
k
ω
,
K
e
k
=
K
11
-
K
12
K
22
-
M
22
k
2
ω
2
+
j
C
22
k
ω
-
1
×
K
21
-
M
21
k
2
ω
2
+
j
C
21
k
ω
and f
-
H
n
1
k is the kth harmonic term of f
-
n
1.For the fundamental harmonic termx
-
H
1, assume that x
-
H
1
=
x
-
L
+
Δ
x
-
H
1, where x
-
L
=
K
-
ω
2
M
+
j
ω
C
-
1
F
1
e
j
ω
t
+
K
-
ω
2
M
-
j
ω
C
-
1

j
F
2
e
-
j
ω
t is the response of the no-rub rotor system. Δ
x
-
H
1 then satisfies the following equation:
(29)
M
Δ
x
-
¨
H
1
+
C
Δ
x
-
˙
H
1
+
K
Δ
x
-
H
1
+
f
-
H
n
1
=
0
.From (29), we get
(30)
M
e
1
x
-
¨
H
11
+
C
e
1
x
-
˙
H
11
+
K
e
k
x
-
H
11
+
f
-
H
n
11
=
M
e
1
x
-
¨
L
1
+
C
e
1
x
-
˙
L
1
+
K
e
k
x
-
L
1
.From (27) and (30), we get
(31)
∑
k
=
0
r
M
e
k
x
-
¨
H
1
k
+
∑
k
=
0
r
C
e
k
x
-
˙
H
1
k
+
∑
k
=
0
r
K
e
k
x
-
H
1
k
+
f
-
n
=
M
e
1
x
-
¨
1
L
+
C
e
1
x
-
˙
1
L
+
K
e
1
x
-
1
L
.The dimensions of (31) are much less than that of the original system and only related to the nonlinear forces. This is because the responses of the linear part are represented by the responses of nonlinear parts using the receptance data. The real part of (31) can be written as
(32)
Re
∑
k
=
0
r
M
e
k
x
-
¨
H
1
k
+
∑
k
=
0
r
C
e
k
x
-
˙
H
1
k
+
∑
k
=
0
r
K
e
k
x
-
H
1
k
+
f
n
x
,
x
˙
=
Re
M
e
1
x
-
¨
1
L
+
C
e
1
x
-
˙
1
L
+
K
e
1
x
-
1
L
.Equation (32) can be solved using the IHB method. After x
-
1 is obtained, x
-
2 can be recovered from x
-
1 using (26).Since the characteristics of local nonlinearity are utilized, the harmonic terms of all DOFs of the system can be represented by the responses of the DOFs with nonlinearities. Therefore, only the nonlinear DOFs are retained and the computation speed increases greatly. In addition, the accuracy of the results is maintained.


## 2.1. The Newmark-β Method and the Fixed Interface Component Mode Synthesis Approach

The Newmark-β method is a commonly used time integration method with good convergence and computational stability. Therefore it is chosen here to represent time domain methods.The equation of motion of a malfunctioning rotor system can be written as(1)
M
x
¨
+
D
x
˙
+
K
x
+
f
n
x
,
x
˙
=
F
t
,
where M, D, and K are the mass, damping (including linear viscous damping and gyroscopic moment), and stiffness matrices, respectively; x is the displacement vector; F
t is the excitation vector caused by the imbalance; f
n
(
x
,
x
˙
) is the nonlinear force vector, which is a function of x and x
˙; the dots above x in (1) denote derivatives with respect to time t.LettingN
x
,
x
˙
=
D
x
˙
+
K
x
+
f
n
x
,
x
˙, (1) can be changed to
(2)
M
x
¨
+
N
x
,
x
˙
=
F
t
.Based on known values ofx
n, x
˙
n, x
¨
n at time t
n, the approximate values x
¨
n
+
1 and x
˙
n
+
1 are [12]
(3)
x
˙
n
+
1
=
a
1
x
n
+
1
-
x
n
-
a
4
x
˙
n
-
a
5
x
¨
n
,
x
¨
n
+
1
=
a
0
x
n
+
1
-
x
n
-
a
2
x
˙
n
-
a
3
x
¨
n
,
where a
0
=
1
/
β
Δ
t
2, a
1
=
γ
/
β
Δ
t, a
2
=
1
/
β
Δ
t, a
3
=
1
/
2
β
-
1, a
4
=
γ
/
β
-
1, and a
5
=
Δ
t
/
2
γ
/
β
-
2.Let(4)
x
˙
~
n
+
1
=
-
a
1
x
n
-
a
4
x
˙
n
-
a
5
x
¨
n
,
x
¨
~
n
+
1
=
-
a
0
x
n
-
a
2
x
˙
n
-
a
3
x
¨
n
.Equation (3) can then be rewritten as
(5)
x
˙
n
+
1
=
x
˙
~
n
+
1
+
a
1
x
n
+
1
,
x
¨
n
+
1
=
x
¨
~
n
+
1
+
a
0
x
n
+
1
.Assuming thatx
n
+
1
i
+
1
=
x
n
+
1
i
+
Δ
x
n
+
1
i, where Δ
x
n
+
1
i is the increment, the above equation is given by
(6)
x
˙
n
+
1
i
+
1
=
x
˙
n
+
1
i
+
a
1
Δ
x
n
+
1
i
,
x
¨
n
+
1
i
+
1
=
x
¨
n
+
1
i
+
a
0
Δ
x
n
+
1
i
.Expanding (2) via a Taylor series around x
n
+
1
i yields
(7)
Δ
x
n
+
1
i
a
0
M
+
a
1
C
(
x
˙
n
+
1
i
)
+
K
x
n
+
1
i
=
F
n
+
1
t
-
M
x
¨
n
+
1
i
-
N
x
˙
n
+
1
i
,
x
n
+
1
i
.WithJ
n
+
1
=
a
0
M
+
a
1
C
(
x
˙
n
+
1
i
)
+
K
(
x
n
+
1
i
), (7) can be changed to
(8)
J
n
+
1
Δ
x
n
+
1
i
=
F
n
+
1
t
-
M
x
¨
n
+
1
i
-
N
x
˙
n
+
1
i
,
x
n
+
1
i
,
where J
n
+
1 is the instantaneous stiffness matrix. Δ
x
n
+
1
i can be obtained using (8) and x
n
+
1
i
+
1
=
x
n
+
1
i
+
Δ
x
n
+
1
i, x
n
+
1 can be obtained when convergence criteria are satisfied.There are four DOFs in one node and if the rotor system is made up ofN nodes, there are a total of 4
N DOFs in the system. The dimension of the instantaneous stiffness matrix J
n
+
1 is 4
N
×
4
N. As N increases, the calculation speed decreases rapidly. The fixed interface component mode synthesis approach can be used to reduce the number of dimensions of the system and thus the simulation efficiency can be increased.The DOFs of the rotor system are separated into two distinct groups. The linear groupx
1 contains the DOFs that are not related to the nonlinear forces. The nonlinear group x
2 contains the DOFs that are related to the nonlinear forces. Equation (2) is rearranged into the following form:
(9)
M
i
i
M
i
j
M
j
i
M
j
j
u
¨
i
u
¨
j
+
C
i
i
C
i
j
C
j
i
C
j
j
u
˙
i
u
˙
j
+
K
i
i
K
i
j
K
j
i
K
j
j
u
i
u
j
=
F
i
F
j
+
0
F
b
.The nonlinear DOFs are fixed and the modes of the system comprised of the remaining DOFs are calculated. The modes that contribute most to the responses of the former rotor system are chosen to construct the fixed interface normal mode set:(10)
ψ
m
=
ψ
l
0
T
.Assuming unit displacement in each nonlinear DOF and setting the movement of the linear DOFs to zero, the constrained mode of the nonlinear DOFs can be obtained:(11)
X
i
=
-
K
i
i
-
1
K
i
j
X
j
.The constrained mode set can then be constructed from all the constrained modes:(12)
ψ
c
=
-
K
i
i
-
1
K
i
j
I
.The modal matrix of the dimension-reduced system can then be obtained by combining the normal mode set and the constrained mode set:(13)
ψ
=
ψ
m
ψ
c
.Using the modal matrix, (9) can be rewritten as:
(14)
M
~
u
¨
+
D
~
u
˙
+
K
~
u
=
F
~
t
+
0
F
b
,
where M
~
=
ψ
T
M
i
i
M
i
j
M
j
i
M
j
j
ψ, D
~
=
ψ
T
D
i
i
D
i
j
D
j
i
D
j
j
ψ, K
~
=
ψ
T
K
i
i
K
i
j
K
j
i
K
j
j
ψ, and F
~
=
ψ
T
F
i
F
j.The responseu can be obtained using Newmark method, and the response x can then be obtained by x
=
ψ
u.In the derivation process above, the dimension reduction mainly depends on the reduction of the normal modes. However, this may result in a loss of accuracy.


## 2.2. Incremental Harmonic Balance Method and Dimension Reduction Using Receptance Data

In the three commonly used frequency domain simulation methods (HB, DF, and IHB), the IHB method is equivalent to the HB method plus the Newton-Raphson method [13]. The IHB method can handle strong nonlinearities and is therefore chosen to represent the frequency domain methods. Here, the multidimensional IHB is introduced because there are many multifrequency excited systems in engineering, such as the dual-rotor system in aeroengines. The detailed information of the single frequency IHB method can be seen in [14, 15].The equation of motion of a multifrequency excited rotor system is expressed as(15)
M
x
¨
+
C
x
˙
+
K
x
+
f
n
x
,
x
˙
=
∑
i
=
1
r
F
c
i
cos
⁡
ω
i
t
+
F
s
i
sin
ω
i
t
,
where ω
i are the exciting frequencies and F
c
i and F
s
i are the amplitudes of the exciting forces corresponding to ω
i.When the rotor system is in steady state, the responsex of the rotor system can be expanded to the form of a multiple Fourier series:
(16)
x
k
=
a
0
k
+
∑
l
1
=
-
R
1
R
1
⋯
∑
l
r
=
-
R
r
R
r
a
l
1
⋯
l
r
k
cos
⁡
∑
i
=
1
r
l
i
ω
i
t
+
b
l
1
⋯
l
r
k
sin
∑
i
=
1
r
l
i
ω
i
t
.Assumingω
=
ω
1
,
…
,
ω
r
T, l
i
=
l
1
,
…
,
l
r, and ν
i
=
l
i
·
ω, (10) can be simplified to
(17)
x
k
=
a
0
k
+
∑
i
=
1
p
c
i
k
cos
⁡
ν
i
t
+
d
i
k
sin
ν
i
t
=
C
s
A
k
,
where C
s
=
1
cos
⁡
ν
1
t
sin
ν
1
t
⋯
cos
⁡
ν
p
t
sin
ν
p
t and A
k
=
a
0
k
c
1
k
d
1
k
⋯
c
p
k
d
p
k
T.LetS
=
diag
⁡
C
s
C
s
⋯
C
s and A
=
A
1
T
A
2
T
⋯
A
N
T
T. Therefore, x
=
S
A. Assume that
(18)
x
=
x
0
+
ε
x
.Substitute (18) into (15) and neglecting the high-order harmonic components, one obtains
(19)
M
S
¨
ε
A
+
C
S
˙
ε
A
+
K
S
ε
A
+
C
n
S
˙
ε
A
+
K
n
S
ε
A
=
-
M
S
¨
A
0
+
C
S
˙
A
0
+
K
S
A
0
+
∑
i
=
1
r
F
c
i
cos
⁡
ω
i
t
+
F
s
i
sin
ω
i
t
-
f
n
x
0
,
where K
n
=
∂
f
n
/
∂
x and C
n
=
∂
f
n
/
∂
x
˙.Letτ
=
ω
1
t
,
…
,
ω
r
t
T
=
τ
1
,
…
,
τ
r
T and applying the Galerkin’s procedure, (19) yields
(20)
∫
0
2
π
⋯
∫
0
2
π
S
T
M
S
¨
+
C
S
˙
+
K
S
+
C
n
S
˙
+
K
n
S
d
τ
1
⋯
d
τ
r
·
ε
A
=
-
∫
0
2
π
⋯
∫
0
2
π
S
T
M
S
¨
+
C
S
˙
+
K
S
d
τ
1
⋯
d
τ
r
·
A
0
+
∫
0
2
π
⋯
∫
0
2
π
S
T
∑
i
=
1
r
F
c
i
cos
⁡
ω
i
t
+
F
s
i
sin
ω
i
t
-
f
n
x
0
∑
i
=
1
r
F
c
i
cos
⁡
ω
i
t
+
F
s
i
sin
ω
i
t
r
d
τ
1
⋯
d
τ
.Equation (20) can be rewritten as
(21)
K
m
Δ
A
=
R
m
1
A
0
+
R
m
2
.A
0 can be obtained from (21) iteratively, following which x can be obtained.K
m is the instantaneous stiffness matrix and has dimensions N
×
2
r
+
1
×
N
×
2
r
+
1. An increase in the DOFs and the number of harmonic components will greatly reduce the computational efficiency of the IHB method.The receptance data is used in the dimension reduction strategy in the IHB method. The DOFs of the rotor system are separated into two distinct groups. The nonlinear groupx
1 contains the DOFs that are related to the nonlinear forces and the linear group x
2 contains the DOFs that are not related to the nonlinear forces. Rearrange (15) to
(22)
M
11
M
12
M
21
M
22
x
¨
1
x
¨
2
+
C
11
C
12
C
21
C
22
x
˙
1
x
˙
2
+
K
11
K
12
K
21
K
22
x
1
x
2
+
f
n
1
x
1
,
x
˙
1
0
=
F
11
F
21
cos
⁡
ω
t
+
F
12
F
22
sin
ω
t
.When the nonlinear system is in steady state, the response vectorx is periodic and the nonlinear force vector f
n
x
,
x
˙ is also periodic, which can be expressed as
(23)
f
n
x
,
x
˙
=
S
P
,
where P
=
P
1
T
P
2
T
⋯
P
h
T
T and P
i
=
p
0
i
p
1
i
q
1
i
⋯
p
r
i
q
r
i
T.Assume thatC
s
e
=
1
e
j
ω
t
j
e
-
j
ω
t
⋯
e
j
r
ω
t
j
e
-
j
r
ω
t, S
e
=
diag
⁡
C
s
e
C
s
e
⋯
C
s
e, and x
-
=
S
e
A. Therefore x
=
Re
x
- and x
- satisfies the following equation:
(24)
M
x
-
¨
+
C
x
-
˙
+
K
x
-
+
f
-
n
=
F
1
e
j
ω
t
+
j
F
2
e
-
j
ω
t
,
where f
-
n
=
S
e
P and f
n
=
Re
f
-
n.Rewrite (24) as
(25)
M
11
M
12
M
21
M
22
x
-
¨
1
x
-
¨
2
+
C
11
C
12
C
21
C
22
x
-
˙
1
x
-
˙
2
+
K
11
K
12
K
21
K
22
x
-
1
x
-
2
+
f
-
n
1
0
=
F
11
F
21
e
j
ω
t
+
j
F
12
F
22
e
j
ω
t
.Assume thatx
-
=
∑
k
=
0
r
x
-
H
k, where x
-
k
=
x
-
H
1
k
x
-
H
2
k
⋯
x
-
H
N
k is the kth order harmonic term of x
-, and x
-
H
i
k
=
p
i
k
e
j
k
ω
t
+
j
q
i
k
e
-
j
k
ω
t. From (25), the kth harmonic term of x
-
H
2
k can be written as
(26)
x
-
H
2
k
=
-
K
22
-
M
22
k
2
ω
2
+
j
C
22
k
ω
-
1
×
K
21
-
M
21
k
2
ω
2
+
j
C
21
k
ω
x
-
H
1
k
.Substituting (26) into the first equation of (25) yields
(27)
M
e
k
x
-
¨
H
1
k
+
C
e
k
x
-
˙
H
1
k
+
K
e
k
x
-
H
1
k
+
f
-
H
n
1
k
=
0
,
where
(28)
M
e
k
=
M
11
-
M
12
K
22
-
M
22
k
2
ω
2
+
j
C
22
k
ω
-
1
×
K
21
-
M
21
k
2
ω
2
+
j
C
21
k
ω
,
C
e
k
=
C
11
-
C
12
K
22
-
M
22
k
2
ω
2
+
j
C
22
k
ω
-
1
×
K
21
-
M
21
k
2
ω
2
+
j
C
21
k
ω
,
K
e
k
=
K
11
-
K
12
K
22
-
M
22
k
2
ω
2
+
j
C
22
k
ω
-
1
×
K
21
-
M
21
k
2
ω
2
+
j
C
21
k
ω
and f
-
H
n
1
k is the kth harmonic term of f
-
n
1.For the fundamental harmonic termx
-
H
1, assume that x
-
H
1
=
x
-
L
+
Δ
x
-
H
1, where x
-
L
=
K
-
ω
2
M
+
j
ω
C
-
1
F
1
e
j
ω
t
+
K
-
ω
2
M
-
j
ω
C
-
1

j
F
2
e
-
j
ω
t is the response of the no-rub rotor system. Δ
x
-
H
1 then satisfies the following equation:
(29)
M
Δ
x
-
¨
H
1
+
C
Δ
x
-
˙
H
1
+
K
Δ
x
-
H
1
+
f
-
H
n
1
=
0
.From (29), we get
(30)
M
e
1
x
-
¨
H
11
+
C
e
1
x
-
˙
H
11
+
K
e
k
x
-
H
11
+
f
-
H
n
11
=
M
e
1
x
-
¨
L
1
+
C
e
1
x
-
˙
L
1
+
K
e
k
x
-
L
1
.From (27) and (30), we get
(31)
∑
k
=
0
r
M
e
k
x
-
¨
H
1
k
+
∑
k
=
0
r
C
e
k
x
-
˙
H
1
k
+
∑
k
=
0
r
K
e
k
x
-
H
1
k
+
f
-
n
=
M
e
1
x
-
¨
1
L
+
C
e
1
x
-
˙
1
L
+
K
e
1
x
-
1
L
.The dimensions of (31) are much less than that of the original system and only related to the nonlinear forces. This is because the responses of the linear part are represented by the responses of nonlinear parts using the receptance data. The real part of (31) can be written as
(32)
Re
∑
k
=
0
r
M
e
k
x
-
¨
H
1
k
+
∑
k
=
0
r
C
e
k
x
-
˙
H
1
k
+
∑
k
=
0
r
K
e
k
x
-
H
1
k
+
f
n
x
,
x
˙
=
Re
M
e
1
x
-
¨
1
L
+
C
e
1
x
-
˙
1
L
+
K
e
1
x
-
1
L
.Equation (32) can be solved using the IHB method. After x
-
1 is obtained, x
-
2 can be recovered from x
-
1 using (26).Since the characteristics of local nonlinearity are utilized, the harmonic terms of all DOFs of the system can be represented by the responses of the DOFs with nonlinearities. Therefore, only the nonlinear DOFs are retained and the computation speed increases greatly. In addition, the accuracy of the results is maintained.


## 3. Numerical Example and Discussion

### 3.1. Single Frequency Excitation

A double-span rotor is comprised of two rotors that are connected by a membrane coupling as shown in Figure1. For the finite element model, the two rotors are divided into beam segments with dimensions that are listed in Table 1. The coupling is represented by a hollow shaft with an outer diameter of 160 mm, an inner diameter of 140 mm, and a length of 500 mm. The hollow shaft is located between nodes 16 and 17.Table 1
Shaft segment dimension of the two-span rotor system.
Segment
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
L(mm)
80
80
290
100
200
240
150
180
160
180
130
100
80
75
75
D (mm)
140
180
210
220
235
250
250
300
240
235
230
190
160
140
140
Segment
16
17
18
19
20
21
22
23
24
25
26
27
28
29
L (mm)
75
75
62
360
90
390
300
500
300
700
80
360
62
100
D (mm)
225
225
270
300
360
420
450
460
450
440
380
310
270
225Figure 1
The model of the two-span rotor system.Assuming that the elastic modulus of the steel material is2
×
1
0
11
Pa, the supporting stiffness of the left and right rotor is 1.5
×
1
0
8
N
/
m and 1
×
1
0
9
N
/
m, respectively. The eccentricities of the left rotor are at nodes 8 and 9, values of m
r
1
=
5
×
1
0
-
3 kg·m, respectively. The eccentricities of the right rotor are at nodes 24 and 25, values of m
r
2
=
4
×
1
0
-
2 kg·m, respectively. The rotor system is assumed to have proportional damping with a mass damping factor of α
0
=
0 and a stiffness damping factor of β
0
=
0.001.Assuming that the full angular rub occurs on node 24 and that the clearance is8
×
1
0
-
5
m, the friction coefficient is 0.15.Figure2(a) shows a comparison of the time response for a rotational speed of 450 rad/s and a contact stiffness of 1
×
1
0
9
N
/
m. The results of the Newmark method without dimension reduction are used as a benchmark. In Figure 2(a), “○” represents the results of the Newmark method with the CMS approach with 5 normal modes, “△” represents the results of the Newmark method with the CMS approach with 10 normal modes, and “×” represents the results of the IHB method with receptance dimension reduction. The time domain methods need a large number of normal modes to attain accuracy targets, whereas frequency domain methods attain very good accuracy more effectively.The rotor system response under single frequency excitation.
(a)
Comparison of the time response
(b)
Amplitude-frequency response curve
(c)
Stable and unstable time responsesThe frequency domain methods combined with the arc-length algorithm [16, 17] can obtain the unstable steady state responses, as shown in Figure 2(b). The curves are the results for contact stiffness values of 5
×
1
0
8
N
/
m, 1
×
1
0
9
N
/
m, and 2
×
1
0
9
N
/
m, respectively. a
b
¯ and c
b
¯ are the unstable steady state responses. Figure 2(c) shows the stable and unstable time responses for a rotational speed of 550 rad/s.A comparison of the efficiency of the two methods is listed in Table2. The computations were performed on a microcomputer operating at 2.83 GHz with 3.25 GB of RAM. The calculation time of the frequency domain method is significantly less than that of the time domain method.Table 2
CPU time for the two methods under single frequency excitation.
Method
Time (s)
Newmark, CMS order 5
14.08
Newmark, CMS order 10
16.11
Newmark, CMS order 20
19.76
IHB,r
=
1
0.29


### 3.2. Double Frequency Excitation

A model of the dual rotor system of an aeroengine is shown in Figure3 with the dimensions of the two rotors listed in Tables 3 and 4.Table 3
The dimensions of the inner rotor.
Segment
1
2
3
4
5
6
7
Node
1
2
2
3
3
4
4
5
5
6
6
7
7
8
Length (mm)
150
100
250
1000
250
100
150
Diameter (mm)
100
500
100
100
100
500
100Table 4
The dimensions of the outer rotor.
Segment
1
2
3
4
5
Node
9
10
10
11
11
12
12
13
13
14
Length (mm)
200
100
400
100
200
External diameter (mm)
200
400
200
400
200
Internal diameter (mm)
150
150
150
150
150Figure 3
Model of the dual-rotor system.Bearings I, II, and III are fixed with supports. Bearing IV is an intershaft bearing that connects the outer and inner rotors. The stiffness of the bearings is7
×
1
0
6
N
/
m each. Since there are imbalances in both rotors, the system experiences double frequency excitation. The eccentricity of the inner rotor is 0.02 kg·m and mounted on node 4; the eccentricity of the outer rotor is 0.05 kg·m and mounted on node 11. The rotor system is assumed to have proportional damping with a mass damping factor of α
0
=
0 and a stiffness damping factor of β
0
=
0.001. Assuming that the rubbing occurs at node 12 and that the clearance is e
=
2
×
1
0
-
4
m, the rubbing is partial rubbing with the rubbing angle of 0–90 degrees. The rotating speed of the outer rotor is 1.4 times that of the inner rotor. Assuming the rotating speed of the inner rotor is ω
1 and that of the outer rotor is ω
2, so ω
2
/
ω
1
=
1.4.The IHB method is compared with the Newmark method to predict the steady state responses of the dual rotor system. When the IHB method is applied, the harmonic terms are chosen based on the method by Ushida and Chua [18] as the following ranges are chosen for this study: -
3
≤
l
1
≤
3 and -
3
≤
l
2
≤
3, l
1
ω
1
+
l
2
ω
2
≥
0. After eliminating the frequencies with negative values, there are all together 25 harmonic terms in the responses.In Figure4, “–” represents the results obtained by the IHB method for a rotating speed of 100
rad
/
s, “△” represents the results obtained by the Newmark method for 5 normal modes, and “×” are the results obtained by the Newmark method for 10 normal modes. The results obtained by the IHB method are almost equal to those obtained by the Newmark with a large number of normal modes. However, the Newmark method can only produce time domain results, whereas the IHB method can produce both time domain and frequency results simultaneously. For example, when the IHB method is used to obtain the responses in a given speed range, the amplitude-frequency responses curve and the three-dimensional spectrum diagram can be readily obtained. When the rotating speed ω
1 varies from 50 rad/s to 350 rad/s and ω
2
/
ω
1
=
1.4, the amplitude-frequency response curve of node 4 is shown in Figure 4(b). The corresponding three-dimensional spectrum diagram is shown in Figure 4(c) and contains information on each harmonic term.The responses under double frequency excitation.
(a)
Comparison of time response
(b)
Amplitude-frequency response curve
(c)
Three-dimensional spectrum diagramThe computational time for each method is listed in Table5.Table 5
CPU time for the two methods under double frequency excitation.
Method
Time (s)
Newmark, CMS order 5
6.49
Newmark, CMS order 10
8.85
Newmark, CMS order 20
10.48
IHB,r
=
25
34.7


### 3.3. Discussion

The simulation results above show that the frequency domain methods are much more efficient than the time domain methods in single frequency excitation condition. This is because of the following.(1)
The frequency domain methods can skip the transient responses. They can obtain the steady state responses directly, while the time domain methods must devote computational time and resources to calculate the transient responses, especially when the damping is small.
(2)
The size of the instantaneous stiffness matrix is very small. When the number of the nonlinear DOFs and harmonic terms is small, the size of the instantaneous stiffness matrix in frequency domain methods is also small. In contrast, the size of the instantaneous stiffness matrix in time domain methods is mainly determined by the number of normal modes.The results also show that the efficiency of the IHB method is very high for single frequency excitation but descends quickly for double frequency excitation. This is because of the following.(1)
Calculating the multiple integral for double frequency excitation is time-consuming. The Galerkin integration procedure in the IHB method for the double frequency excitation involves a multiple integral process, which is very time-consuming. For example, for single frequency excitation, one period is divided into 100 subintervals in both the Newmark and IHB method. Therefore, the two methods complete the integration over the same time frame. However, for double frequency excitation, 10,000 subintervals must be created for the multiple integral of the IHB method, which results in much higher time costs than in the Newmark method.
(2)
The dimensions of instantaneous stiffness matrix are very high in double frequency excitation when using the IHB method. The dimensions of the instantaneous stiffness matrix increase with the increase in the number of harmonic terms, which is much higher for double frequency excitation than for single frequency excitation. For example, when-
5
≤
l
1
≤
5 and -
5
≤
l
2
≤
5, ω
2
/
ω
1
=
1.4, the number of harmonic terms is r
=
49 and the size of the instantaneous stiffness matrix is 99
N
×
99
N. This is much larger than the dimensions of 11
N
×
11
N for single frequency excitation. For the Newmark method, the dimensions of the instantaneous stiffness matrix remain the same.If the number of DOFs of the rotor system is larger with fewer nonlinear DOFs, and the damping is small, the frequency domain methods can still have a high efficiency. Otherwise, time domain methods are more efficient.There are some rotor systems in practice, such as the geared supercharger, that experience three or more than three frequencies excitations. If frequency domain methods are applied, the multiple integral and the large instantaneous stiffness matrix will require too much computational time to determine the response. Therefore, time domain methods are most appropriate in those scenarios.


## 3.1. Single Frequency Excitation

A double-span rotor is comprised of two rotors that are connected by a membrane coupling as shown in Figure1. For the finite element model, the two rotors are divided into beam segments with dimensions that are listed in Table 1. The coupling is represented by a hollow shaft with an outer diameter of 160 mm, an inner diameter of 140 mm, and a length of 500 mm. The hollow shaft is located between nodes 16 and 17.Table 1
Shaft segment dimension of the two-span rotor system.
Segment
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
L(mm)
80
80
290
100
200
240
150
180
160
180
130
100
80
75
75
D (mm)
140
180
210
220
235
250
250
300
240
235
230
190
160
140
140
Segment
16
17
18
19
20
21
22
23
24
25
26
27
28
29
L (mm)
75
75
62
360
90
390
300
500
300
700
80
360
62
100
D (mm)
225
225
270
300
360
420
450
460
450
440
380
310
270
225Figure 1
The model of the two-span rotor system.Assuming that the elastic modulus of the steel material is2
×
1
0
11
Pa, the supporting stiffness of the left and right rotor is 1.5
×
1
0
8
N
/
m and 1
×
1
0
9
N
/
m, respectively. The eccentricities of the left rotor are at nodes 8 and 9, values of m
r
1
=
5
×
1
0
-
3 kg·m, respectively. The eccentricities of the right rotor are at nodes 24 and 25, values of m
r
2
=
4
×
1
0
-
2 kg·m, respectively. The rotor system is assumed to have proportional damping with a mass damping factor of α
0
=
0 and a stiffness damping factor of β
0
=
0.001.Assuming that the full angular rub occurs on node 24 and that the clearance is8
×
1
0
-
5
m, the friction coefficient is 0.15.Figure2(a) shows a comparison of the time response for a rotational speed of 450 rad/s and a contact stiffness of 1
×
1
0
9
N
/
m. The results of the Newmark method without dimension reduction are used as a benchmark. In Figure 2(a), “○” represents the results of the Newmark method with the CMS approach with 5 normal modes, “△” represents the results of the Newmark method with the CMS approach with 10 normal modes, and “×” represents the results of the IHB method with receptance dimension reduction. The time domain methods need a large number of normal modes to attain accuracy targets, whereas frequency domain methods attain very good accuracy more effectively.The rotor system response under single frequency excitation.
(a)
Comparison of the time response
(b)
Amplitude-frequency response curve
(c)
Stable and unstable time responsesThe frequency domain methods combined with the arc-length algorithm [16, 17] can obtain the unstable steady state responses, as shown in Figure 2(b). The curves are the results for contact stiffness values of 5
×
1
0
8
N
/
m, 1
×
1
0
9
N
/
m, and 2
×
1
0
9
N
/
m, respectively. a
b
¯ and c
b
¯ are the unstable steady state responses. Figure 2(c) shows the stable and unstable time responses for a rotational speed of 550 rad/s.A comparison of the efficiency of the two methods is listed in Table2. The computations were performed on a microcomputer operating at 2.83 GHz with 3.25 GB of RAM. The calculation time of the frequency domain method is significantly less than that of the time domain method.Table 2
CPU time for the two methods under single frequency excitation.
Method
Time (s)
Newmark, CMS order 5
14.08
Newmark, CMS order 10
16.11
Newmark, CMS order 20
19.76
IHB,r
=
1
0.29


## 3.2. Double Frequency Excitation

A model of the dual rotor system of an aeroengine is shown in Figure3 with the dimensions of the two rotors listed in Tables 3 and 4.Table 3
The dimensions of the inner rotor.
Segment
1
2
3
4
5
6
7
Node
1
2
2
3
3
4
4
5
5
6
6
7
7
8
Length (mm)
150
100
250
1000
250
100
150
Diameter (mm)
100
500
100
100
100
500
100Table 4
The dimensions of the outer rotor.
Segment
1
2
3
4
5
Node
9
10
10
11
11
12
12
13
13
14
Length (mm)
200
100
400
100
200
External diameter (mm)
200
400
200
400
200
Internal diameter (mm)
150
150
150
150
150Figure 3
Model of the dual-rotor system.Bearings I, II, and III are fixed with supports. Bearing IV is an intershaft bearing that connects the outer and inner rotors. The stiffness of the bearings is7
×
1
0
6
N
/
m each. Since there are imbalances in both rotors, the system experiences double frequency excitation. The eccentricity of the inner rotor is 0.02 kg·m and mounted on node 4; the eccentricity of the outer rotor is 0.05 kg·m and mounted on node 11. The rotor system is assumed to have proportional damping with a mass damping factor of α
0
=
0 and a stiffness damping factor of β
0
=
0.001. Assuming that the rubbing occurs at node 12 and that the clearance is e
=
2
×
1
0
-
4
m, the rubbing is partial rubbing with the rubbing angle of 0–90 degrees. The rotating speed of the outer rotor is 1.4 times that of the inner rotor. Assuming the rotating speed of the inner rotor is ω
1 and that of the outer rotor is ω
2, so ω
2
/
ω
1
=
1.4.The IHB method is compared with the Newmark method to predict the steady state responses of the dual rotor system. When the IHB method is applied, the harmonic terms are chosen based on the method by Ushida and Chua [18] as the following ranges are chosen for this study: -
3
≤
l
1
≤
3 and -
3
≤
l
2
≤
3, l
1
ω
1
+
l
2
ω
2
≥
0. After eliminating the frequencies with negative values, there are all together 25 harmonic terms in the responses.In Figure4, “–” represents the results obtained by the IHB method for a rotating speed of 100
rad
/
s, “△” represents the results obtained by the Newmark method for 5 normal modes, and “×” are the results obtained by the Newmark method for 10 normal modes. The results obtained by the IHB method are almost equal to those obtained by the Newmark with a large number of normal modes. However, the Newmark method can only produce time domain results, whereas the IHB method can produce both time domain and frequency results simultaneously. For example, when the IHB method is used to obtain the responses in a given speed range, the amplitude-frequency responses curve and the three-dimensional spectrum diagram can be readily obtained. When the rotating speed ω
1 varies from 50 rad/s to 350 rad/s and ω
2
/
ω
1
=
1.4, the amplitude-frequency response curve of node 4 is shown in Figure 4(b). The corresponding three-dimensional spectrum diagram is shown in Figure 4(c) and contains information on each harmonic term.The responses under double frequency excitation.
(a)
Comparison of time response
(b)
Amplitude-frequency response curve
(c)
Three-dimensional spectrum diagramThe computational time for each method is listed in Table5.Table 5
CPU time for the two methods under double frequency excitation.
Method
Time (s)
Newmark, CMS order 5
6.49
Newmark, CMS order 10
8.85
Newmark, CMS order 20
10.48
IHB,r
=
25
34.7


## 3.3. Discussion

The simulation results above show that the frequency domain methods are much more efficient than the time domain methods in single frequency excitation condition. This is because of the following.(1)
The frequency domain methods can skip the transient responses. They can obtain the steady state responses directly, while the time domain methods must devote computational time and resources to calculate the transient responses, especially when the damping is small.
(2)
The size of the instantaneous stiffness matrix is very small. When the number of the nonlinear DOFs and harmonic terms is small, the size of the instantaneous stiffness matrix in frequency domain methods is also small. In contrast, the size of the instantaneous stiffness matrix in time domain methods is mainly determined by the number of normal modes.The results also show that the efficiency of the IHB method is very high for single frequency excitation but descends quickly for double frequency excitation. This is because of the following.(1)
Calculating the multiple integral for double frequency excitation is time-consuming. The Galerkin integration procedure in the IHB method for the double frequency excitation involves a multiple integral process, which is very time-consuming. For example, for single frequency excitation, one period is divided into 100 subintervals in both the Newmark and IHB method. Therefore, the two methods complete the integration over the same time frame. However, for double frequency excitation, 10,000 subintervals must be created for the multiple integral of the IHB method, which results in much higher time costs than in the Newmark method.
(2)
The dimensions of instantaneous stiffness matrix are very high in double frequency excitation when using the IHB method. The dimensions of the instantaneous stiffness matrix increase with the increase in the number of harmonic terms, which is much higher for double frequency excitation than for single frequency excitation. For example, when-
5
≤
l
1
≤
5 and -
5
≤
l
2
≤
5, ω
2
/
ω
1
=
1.4, the number of harmonic terms is r
=
49 and the size of the instantaneous stiffness matrix is 99
N
×
99
N. This is much larger than the dimensions of 11
N
×
11
N for single frequency excitation. For the Newmark method, the dimensions of the instantaneous stiffness matrix remain the same.If the number of DOFs of the rotor system is larger with fewer nonlinear DOFs, and the damping is small, the frequency domain methods can still have a high efficiency. Otherwise, time domain methods are more efficient.There are some rotor systems in practice, such as the geared supercharger, that experience three or more than three frequencies excitations. If frequency domain methods are applied, the multiple integral and the large instantaneous stiffness matrix will require too much computational time to determine the response. Therefore, time domain methods are most appropriate in those scenarios.


## 4. Conclusions

Time domain and frequency domain methods are used to predict the responses of rotor system with malfunctions. The efficiency and accuracy of the two types of methods are compared based on the excitation type. The main conclusions are as follows.(1)
The dimension reduction approach in time domain methods is a substructure approach. The high order modes are truncated, which may affect the accuracy of the methods. Frequency domain methods that use the receptance data for dimension reduction can retain accuracy.
(2)
The efficiency of the frequency domain methods is determined by the number of excitation frequencies and the number of harmonic terms. These two factors determine the complexity of the numerical integration process and the size of the instantaneous stiffness matrix. Frequency domain methods are highly efficient for single frequency excitations but suffer a decline in efficiency for multiple frequency excitations. The size of the instantaneous stiffness matrix remains the same in time domain methods. Therefore, the efficiency remains the same for both single and multifrequency excitation.
(3)
Frequency domain methods are recommended for single frequency excitation. Time domain methods are recommended in cases where there are three or more frequency excitations. For double frequency excitations, the choice of simulation method depends on the specific conditions.



---

*Source: 126843-2014-12-09.xml*