# Transient Analysis of Speed-Varying Rotor with Uncertainty Based on Interval Approaches

Authors: Tao Wei; WenXuan Gou; Chao Fu; YongFeng Yang
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

The nonintrusive Chebyshev interval method is used to analyze the transient response of a rotor system. The main idea is to analyze the uncertainty in speed-varying rotor systems and reduce the error in the calculation. The hybrid analysis procedure, which combines the black-box modeling and overestimation controlling, is proposed via the Chebyshev approximation and searching the extreme values of the approximation function for the uncertain response. It can both increase accuracy and improve the computational efficiency of the proposed method. Effects of the uncertain parameters on response bounds and estimation errors of the system are investigated, e.g., supporting stiffness and damping. The accuracy and efficiency of the proposed method are verified via comparing with the scanning method. Overestimation effects caused by the interval arithmetic can be controlled, particularly the sensitive domains near the critical speed. It can provide guidance for general uncertain transient mechanical systems, particularly those that have severe fluctuations near resonances.


---

## Body

## 1. Introduction

Speed-varying of rotor systems is common due to startup, shutdown, and the maneuver operations. Useful information is contained in the speed-varying transient vibrations. In the past few decades, various studies on transient dynamics of rotor systems have been carried out [1–8]. In practice, complex mechanical systems are subject to variability and uncertainty. For rotor systems, the material properties, geometries, and exterior environments are not always deterministic because of deterioration and errors. Uncertainty propagations may lead to major deviations of the dynamic characteristics [9].Uncertainties in dynamical systems have already aroused the attention of researchers in recent years [10, 11]. In order to evaluate the effects of uncertainty on rotor dynamics, researchers worldwide have carried out many studies concerning various aspects including faults [11, 12], modeling techniques [13], and dynamical balancing [14]. Probabilistic methods are widely employed such as sampling methods [15, 16] and the Polynomial Chaos Expansion (PCE) [11, 17–19]. The traditional sampling method can provide accurate results with sufficient samples but great computational efforts are required. The PCE, which represents the uncertainties by truncated polynomials, is nonintrusive and frequently used. Fuzzy-set algorithms were also adopted by some researchers [20, 21]. Liao [22] proposed a global resonance optimization method for uncertainty quantification and applied it to rotor systems. Recently, the interval-based methods have aroused wide attention due to their minimum requirements of distributions of uncertain parameters compared with probabilistic methods and fuzzy-set methods. Qiu et al. [23, 24] studied both static and dynamic structural problems with uncertain parameters using Taylor interval methods and perturbation techniques. Qi et al. [25] proposed the collocation interval analysis method, which is suitable for finite element models and controls the overestimation effectively. Wu et al. [26] introduced the Chebyshev inclusion function into nonlinear dynamic problems, which achieved high efficiency. In this scheme, the interval uncertainties are approximated by orthogonal Chebyshev polynomials and a surrogate model is established for the original dynamic input-output system. Many applications can been found due to its efficiency and nonintrusive property [27, 28]. As suggested in literature [26], the Chebyshev interval method (CIM) leads to some overestimation introduced by the interval arithmetic and further algorithms should be incorporated to reduce it. Liu et al. [29] developed a hybrid interval method for uncertain dynamic systems using the derivative information, which exhibited better performance. Jiang et al. [30, 31] proposed time-varying uncertainty analysis approaches based on interval methods. Up until now, the interval methods have not been applied to uncertain transient rotor dynamics. Early effort was made by Dimarogonas [32]. Ma et al. [33, 34] studied the steady state response of rotor systems using interval modal superposition method and perturbation skills. Fu et al. [35] studied the dynamics of a hollow shaft overhung rotor system based on a nonintrusive interval method. The uncertain transient dynamics of rotor system with interval parameters remain misunderstood. Rigorous enclosure solutions of uncertain transient problems are hard to define by nonintrusive interval methods. The motive of this paper is to develop a hybrid interval procedure based on the Chebyshev approximation for transient calculation of rotors. In such a method, the overestimation is effectively controlled and the efficiency is satisfactory.This paper is organized as follows. The time-varying rotor model and deterministic transient response are described in Section2. Considering parametric uncertainties, the interval transient problem is explained in Section 3. In Section 4, the nonintrusive efficient Chebyshev interval method is explained and the corresponding numerical results are given. A hybrid procedure is proposed in Section 5 for overestimation control and the recalculated responses are presented. Conclusions are drawn in the last section.


## 2. Transient Analysis for Deterministic Rotor Model

### 2.1. Rotor System

The rotor system considered in this paper is comprised of a shaft supported by two bearings and a horizontal disc, as shown in Figure1(a). The lumped mass disc is located between the two supports with mass m, transverse mass moment of inertia I
d, and polar mass moment of inertia I
p. Detailed physical parameters and values are given in Table 1. The FEM is used to discretize the rotor into elements. Each element has two nodes and each node has four degrees of freedom in fixed rotating frame shown in Figure 1(b). Figure 1(b) is reproduced from Fu et al. (2017).Table 1
Parameters of the rotor system.
Symbol
Description
Value
E
Young’s modulus of elasticity
2.1
×
1
0
11 N/m2
L
1
Length of shaft 1
0.05 m
L
2
Length of shaft 2
0.96 m
L
3
Length of shaft 3
0.24 m
L
4
Length of shaft 4
0.03 m
D
s
Diameter of shafts
0.03 m
K
1
Stiffness of support 1
1
×
1
0
8 N/m
K
2
Stiffness of support 2
1
×
1
0
6 N/m
C
Viscous damping
200 N⋅s/m
m
Disc mass
20 kg
ρ
Density
7800 kg/m3
I
p
Polar mass moment of inertia
0.6250 kg⋅m2
I
d
Transverse mass moment of inertia
0.3125 kg⋅m2Figure 1
Schematic diagram of the rotor system: (a) rotor model and (b) shaft element coordinates in fixed frame.
(a)
(b)For the rotor system, the deterministic general equation of motion can be derived by assembling all the elements together. Generally, it can be expressed in second-order ODEs as(1)
M
U
¨
t
+
C
^
U
˙
t
+
K
U
t
=
F
t
,where M, K, and F represent the mass, stiffness, and external excitations matrices of the rotor system, respectively. Generalized damping matrix C
^
=
C
+
ω
G represents the damping and gyroscopic effects. Symbol ω represents the revolution speed of the rotor, which is time-variant. U, U
˙, and U
¨ are the displacement, velocity, and acceleration vectors, respectively. The excitations caused by mass imbalance can be expressed by(2)
f
t
=
m
e
ω
2
cos
⁡
ω
t
+
φ
0
+
a
sin
⁡
ω
t
+
φ
0
0
m
e
ω
2
sin
⁡
ω
t
+
φ
0
-
a
cos
⁡
ω
t
+
φ
0
0
,with φ
0 being the initial angle of the imbalance and a being the angular acceleration of the rotating system.


### 2.2. Deterministic Accelerating of Transient Response

Assume the initial state of the system is static and the initial angleφ
0
=
90
°. In order to clearly illustrate the uncertain effects in the forthcoming sections, the disk mass imbalance eccentricity is set to a relatively large value, i.e., e
=
0.001 m. Herein, we are interested in the vibration of the disk geometric center. Based on the deterministic parameters shown in Table 1, the accelerating of dynamic response of the rotor system with a constant angular acceleration a
=
40 rad/s2 is plotted in Figure 2. The transient vibration peak appears at ω
=
123.8 rad/s with an amplitude of 3.324
×
1
0
-
3 m.Figure 2
Deterministic accelerating of transient response of the rotor.


## 2.1. Rotor System

The rotor system considered in this paper is comprised of a shaft supported by two bearings and a horizontal disc, as shown in Figure1(a). The lumped mass disc is located between the two supports with mass m, transverse mass moment of inertia I
d, and polar mass moment of inertia I
p. Detailed physical parameters and values are given in Table 1. The FEM is used to discretize the rotor into elements. Each element has two nodes and each node has four degrees of freedom in fixed rotating frame shown in Figure 1(b). Figure 1(b) is reproduced from Fu et al. (2017).Table 1
Parameters of the rotor system.
Symbol
Description
Value
E
Young’s modulus of elasticity
2.1
×
1
0
11 N/m2
L
1
Length of shaft 1
0.05 m
L
2
Length of shaft 2
0.96 m
L
3
Length of shaft 3
0.24 m
L
4
Length of shaft 4
0.03 m
D
s
Diameter of shafts
0.03 m
K
1
Stiffness of support 1
1
×
1
0
8 N/m
K
2
Stiffness of support 2
1
×
1
0
6 N/m
C
Viscous damping
200 N⋅s/m
m
Disc mass
20 kg
ρ
Density
7800 kg/m3
I
p
Polar mass moment of inertia
0.6250 kg⋅m2
I
d
Transverse mass moment of inertia
0.3125 kg⋅m2Figure 1
Schematic diagram of the rotor system: (a) rotor model and (b) shaft element coordinates in fixed frame.
(a)
(b)For the rotor system, the deterministic general equation of motion can be derived by assembling all the elements together. Generally, it can be expressed in second-order ODEs as(1)
M
U
¨
t
+
C
^
U
˙
t
+
K
U
t
=
F
t
,where M, K, and F represent the mass, stiffness, and external excitations matrices of the rotor system, respectively. Generalized damping matrix C
^
=
C
+
ω
G represents the damping and gyroscopic effects. Symbol ω represents the revolution speed of the rotor, which is time-variant. U, U
˙, and U
¨ are the displacement, velocity, and acceleration vectors, respectively. The excitations caused by mass imbalance can be expressed by(2)
f
t
=
m
e
ω
2
cos
⁡
ω
t
+
φ
0
+
a
sin
⁡
ω
t
+
φ
0
0
m
e
ω
2
sin
⁡
ω
t
+
φ
0
-
a
cos
⁡
ω
t
+
φ
0
0
,with φ
0 being the initial angle of the imbalance and a being the angular acceleration of the rotating system.


## 2.2. Deterministic Accelerating of Transient Response

Assume the initial state of the system is static and the initial angleφ
0
=
90
°. In order to clearly illustrate the uncertain effects in the forthcoming sections, the disk mass imbalance eccentricity is set to a relatively large value, i.e., e
=
0.001 m. Herein, we are interested in the vibration of the disk geometric center. Based on the deterministic parameters shown in Table 1, the accelerating of dynamic response of the rotor system with a constant angular acceleration a
=
40 rad/s2 is plotted in Figure 2. The transient vibration peak appears at ω
=
123.8 rad/s with an amplitude of 3.324
×
1
0
-
3 m.Figure 2
Deterministic accelerating of transient response of the rotor.


## 3. Transient Dynamic Problem of Rotor System including Uncertainty

By taking the uncertainties in the rotor system into consideration, the general equation of motion expressed in (1) should be further modified. For the sake of clarity, we collect all the interval parameters and define an uncertain parameter vector q
=
[
q
1
,
q
2
,
…
,
q
r
]
T with r representing the uncertain dimension. If the uncertainties lie in the initial state of the rotor, they can also be included in q and the solution procedures are valid. Then, we can rewrite the motion equations of the system(3)
M
q
U
¨
q
,
t
+
C
^
q
U
˙
q
,
t
+
K
q
U
q
,
t
=
F
q
,
t
.It can be interpreted that the solution of (3) is that of (1) subject to parameter constraints of vector q. Due to the presence of uncertainty, the deterministic response will evolve into a response range defined by lower bound and upper bound. The uncertain response can be expressed as(4)
U
q
,
t
=
U
:
M
U
¨
+
C
U
˙
+
K
U
=
F
,
U
t
0
=
U
0
,
q
∈
q
I
,where the superscript I defines an interval character; U
0 is the initial state vector. The parameter interval vector is expressed as(5)
q
I
=
q
_
,
q
-
=
q
i
:
q
_
i
≤
q
i
≤
q
-
i
,
i
=
1,2
,
…
,
r
,in which underline denotes the lower bound of the character and overbar denotes the upper bound, respectively. The transient response range can be defined by(6)
U
_
=
min
q
∈
q
_
,
q
-
⁡
U
:
M
U
¨
+
C
U
˙
+
K
U
=
F
,
U
t
0
=
U
0
,and(7)
U
-
=
max
q
∈
q
_
,
q
-
⁡
U
:
M
U
¨
+
C
U
˙
+
K
U
=
F
,
U
t
0
=
U
0
.At this stage, we have expressed that uncertain dynamic problem and its solution goals. In general, the ODEs are solved by numerical methods, which discretize the integration process into small time steps. Given the initial state and the time step, we consider the equations as an initial value problem.U
1
I is the solution of the ODEs with U
0 and q
I. For every time step forward, the transient response interval U
t
k
+
1
I only depends on U
t
k
I and uncertain parameter vector q
I. When U
t
k
I is already obtained, U
t
k
+
1
I is the interval solution of initial integration problem within time t to t
k
+
1. To get the minimum and maximum responses, global optimization methods could be incorporated which shall introduce large amount of computational efforts. In the next section, we extend the Chebyshev interval method for prediction of the response range U
t
k
+
1
I using already obtained U
t
k
I under constraints q
I.


## 4. Transient Analysis Using the Chebyshev Interval Method

The Chebyshev orthogonal polynomials are used to approximate the response model with simple mathematical functions. It is similar to a surrogate model and the response range is obtained from the approximated function. Uncertain dynamic problems are transformed into a set of deterministic ODEs with a few interpolation points. In this section, we take a brief review of the CIM [26] and then give the results by using the interval method.


### 4.1. Nonintrusive Chebyshev Interval Method

For one uncertain dimension problem defined within interval[
x
_
,
x
-
], the single-dimensional n-order Chebyshev approximation of the response model can be expressed as(8)
P
n
x
=
a
0
2
+
∑
i
=
1
n
a
i
C
i
x
,
x
∈
x
_
,
x
-
,where a
i are the expansion coefficients. C
i represents the Chebyshev basis which can be defined by(9)
C
i
x
=
cos
⁡
i
θ
,
θ
=
arccos
⁡
2
x
-
x
_
+
x
-
x
-
-
x
_
∈
0
,
π
.The expansion coefficients can be calculated by Mehler integration. It is expressed as(10)
a
i
=
2
λ
∑
j
=
1
λ
U
~
t
k
cos
⁡
θ
j
cos
⁡
i
θ
j
,in which λ is the total number of interpolations. U
~
t
k
(
cos
⁡
θ
j
) denotes the deterministic response calculated from the original rotor model at the interpolation point θ
j which is defined as the zeros of Chebyshev series(11)
θ
j
=
2
j
-
1
2
λ
π
,
j
=
1,2
,
…
,
λ
.When multiple uncertain parameters are considered, the interval approximation should be expressed in tensor product form. TheL-dimensional n-order approximation of (4) is given by(12)
U
t
k
+
1
I
=
∑
k
1
=
0
n
⋯
∑
k
r
=
0
n
1
2
p
Γ
k
1
,
…
,
k
r
C
k
1
,
…
,
k
r
q
,where p denotes the number of zeros that appeared in the subscripts; n is the expansion order. C
k
1
,
…
,
k
r
(
q
I
) represents the multidimensional Chebyshev polynomials which can be defined as(13)
C
k
1
,
…
,
k
r
q
I
=
cos
⁡
k
1
θ
1
cos
⁡
k
2
θ
2
⋯
cos
⁡
k
r
θ
r
.Notation Γ
k
1
,
…
,
k
r is the expansion coefficient, which can be calculated by the Mehler integration. The following expression can be obtained:(14)
Γ
k
1
,
…
,
k
r
=
∑
j
1
=
0
λ
⋯
∑
j
r
=
0
λ
2
λ
r
U
~
t
k
cos
⁡
k
1
θ
j
1
⋯
cos
⁡
k
r
θ
j
r
.Here λ represents the number of the interpolation points used. It is required that λ
≥
n
+
1 in order to be accurate for the numerical integration. U
~
t
k denotes the deterministic transient response when the uncertain parameters take the values as the interpolation points.Considering the bounded characteristics of trigonometric functions, the final transient response range is given by(15)
U
t
k
+
1
I
=
1
2
r
Γ
0
,
…
,
0
+
-
1,1
∑
0
≤
k
1
,
…
,
k
r
≤
r
k
1
+
⋯
+
k
r
≥
1
1
2
p
Γ
k
1
,
…
,
k
r
.This solution procedure is nonintrusive and with high efficiency. But as demonstrated in [26], there is overestimation due to the interval arithmetic of the expansion coefficients in (15). It should be noted that, for some parameters, the overestimation in the transient uncertain dynamic problems is significant in high-sensitive areas, which will be presented next.


### 4.2. Uncertain Transient Response Simulation Using the CIM

In this subsection, we present the numerical results of several cases with different uncertain parameters by using the CIM. To provide a reference for the accuracy and efficiency, the results are compared with those obtained from the scanning method. The expansion order in the CIM is 3 and the interpolation number is 4. The sampling points are 1000 for the scanning method (SM). Three kinds of uncertainties are included here for different occasions in engineering, namely, uncertain damping, mass property, and support stiffness. A typical uncertain deviation coefficient of 10% is adopted for all cases in this paper. The uncertain intervals of parameters are obtained by taking symmetric percentage dispersion around the nominal values, which can be expressed as(16)
q
j
I
=
q
_
j
,
q
-
j
=
q
j
c
-
β
q
j
c
,
q
j
c
+
β
q
j
c
,
j
=
1,2
,
…
,
r
,where q
j is the jth member of the uncertain parameter vector q and β is the deviation coefficient. q
_
j, q
-
j, and q
j
c are the lower bound, upper bound, and nominal value, respectively. Figures 3–5 give the uncertain transient response ranges of the rotor with different interval parameters and the corresponding errors considering the scanning results for reference. From these figures, we can observe that the uncertain effects of interval parameters are quite significant. Moreover, the mechanisms are quite different as well. The uncertain mass unbalance may cause a linear deviation of the deterministic curve while the others produce a large-amplitude vibration band [35].Figure 3
(a) Transient response of the rotor with uncertain mass unbalance and (b) calculation error.
(a)
(b)Figure 4
(a) Transient response of the rotor with uncertain support stiffness 2 and (b) calculation error.
(a)
(b)Figure 5
(a) Transient response of the rotor with uncertain disk mass and (b) calculation error.
(a)
(b)We use Lenovo Thinkpad T430 to calculate the response of speed-varying rotor. The CPU is Intel Core i7-3520 M @ 2.9 GHz, RAM is 12.0 GB, and the OS is Windows 10. The average computing time for the CIM is 6.9 s for each uncertain case and 25 minutes for the scanning method, and the advantage will be more obvious in higher uncertain dimensions. We then can give credit to the efficiency of the CIM in solving uncertain dynamic problems. Figure3 shows that the CIM can provide satisfactory estimated bounds when the damping is uncertain comparing with the SM and the calculation error is small. In that case, we consider the CIM to be fairly accurate. That is true for parameters that do not influence the intrinsic dynamic characteristics such as critical speeds and modes of the rotor system and the disk imbalance eccentricity, for another example. However, for the sensitive parameters such as the bearing stiffness and disk mass property, the CIM causes relatively large overestimations as demonstrated in Figures 4 and 5. Multipeaks are observed which are not in accordance with facts. CIM even gives negative value estimations in Figure 5 due to the high sensitivity of the transient response to the disk mass.Though the uncertainty in viscous damping itself does not affect much of the calculation accuracy, in order to investigate the influence of damping magnitude on the overestimation effect, several cases are chosen where different damping is considered when the disk mass is uncertain given its high sensitivity. The response bounds calculated by both the CIM and the SM are listed in Table2 for different damping at a specific revolution speed 128 rad/s near the resonant peak. In Table 2, lower bound is abbreviated as LB and UB is short for upper bound.Table 2
Response bounds of the system at 128 rad/s via different methods (×10−3 m).
Viscous damping (N⋅s/m)
CIM
SM
Error
LB
UB
LB
UB
LB
UB
100
1.915
6.152
1.982
5.090
0.0678
-1.062
200
1.044
4.025
1.109
3.301
0.0647
-0.723
300
1.085
2.742
1.156
2.370
0.0706
-0.373
400
1.116
2.018
1.160
1.824
0.0436
-0.194
500
1.083
1.586
1.107
1.477
0.0234
-0.109It can be observed form Table2 that the magnitude of damping can influence the response range directly when the disk mass is uncertain. Increasing damping decreases the magnitude of the resonant peak, uncertain response range width, and the calculation error. The error is getting smaller as the damping increases. Another fact shown in Table 2 is that the errors in the upper bounds are greater than that of the lower ones and the lower bounds fell less than the upper bounds. The bounds may not be strictly accurate because of finite sampling times, so that in some occasions, the lower bounds in greater damping are higher than the ones in smaller damping, but the absolute value is minor.Note a general fact that the overestimations are mostly distributed in the resonant regions. The largest overestimation point leads to nearly 1 mm, which is unacceptable for engineering practice. It is stated before that the overestimation is introduced by the interval arithmetic in the CIM, which is the key of the high efficiency of the interval method. To reduce the overestimation of the transient response prediction in the sensitive regions while keeping much of the efficiency, further algorithms must be incorporated.


## 4.1. Nonintrusive Chebyshev Interval Method

For one uncertain dimension problem defined within interval[
x
_
,
x
-
], the single-dimensional n-order Chebyshev approximation of the response model can be expressed as(8)
P
n
x
=
a
0
2
+
∑
i
=
1
n
a
i
C
i
x
,
x
∈
x
_
,
x
-
,where a
i are the expansion coefficients. C
i represents the Chebyshev basis which can be defined by(9)
C
i
x
=
cos
⁡
i
θ
,
θ
=
arccos
⁡
2
x
-
x
_
+
x
-
x
-
-
x
_
∈
0
,
π
.The expansion coefficients can be calculated by Mehler integration. It is expressed as(10)
a
i
=
2
λ
∑
j
=
1
λ
U
~
t
k
cos
⁡
θ
j
cos
⁡
i
θ
j
,in which λ is the total number of interpolations. U
~
t
k
(
cos
⁡
θ
j
) denotes the deterministic response calculated from the original rotor model at the interpolation point θ
j which is defined as the zeros of Chebyshev series(11)
θ
j
=
2
j
-
1
2
λ
π
,
j
=
1,2
,
…
,
λ
.When multiple uncertain parameters are considered, the interval approximation should be expressed in tensor product form. TheL-dimensional n-order approximation of (4) is given by(12)
U
t
k
+
1
I
=
∑
k
1
=
0
n
⋯
∑
k
r
=
0
n
1
2
p
Γ
k
1
,
…
,
k
r
C
k
1
,
…
,
k
r
q
,where p denotes the number of zeros that appeared in the subscripts; n is the expansion order. C
k
1
,
…
,
k
r
(
q
I
) represents the multidimensional Chebyshev polynomials which can be defined as(13)
C
k
1
,
…
,
k
r
q
I
=
cos
⁡
k
1
θ
1
cos
⁡
k
2
θ
2
⋯
cos
⁡
k
r
θ
r
.Notation Γ
k
1
,
…
,
k
r is the expansion coefficient, which can be calculated by the Mehler integration. The following expression can be obtained:(14)
Γ
k
1
,
…
,
k
r
=
∑
j
1
=
0
λ
⋯
∑
j
r
=
0
λ
2
λ
r
U
~
t
k
cos
⁡
k
1
θ
j
1
⋯
cos
⁡
k
r
θ
j
r
.Here λ represents the number of the interpolation points used. It is required that λ
≥
n
+
1 in order to be accurate for the numerical integration. U
~
t
k denotes the deterministic transient response when the uncertain parameters take the values as the interpolation points.Considering the bounded characteristics of trigonometric functions, the final transient response range is given by(15)
U
t
k
+
1
I
=
1
2
r
Γ
0
,
…
,
0
+
-
1,1
∑
0
≤
k
1
,
…
,
k
r
≤
r
k
1
+
⋯
+
k
r
≥
1
1
2
p
Γ
k
1
,
…
,
k
r
.This solution procedure is nonintrusive and with high efficiency. But as demonstrated in [26], there is overestimation due to the interval arithmetic of the expansion coefficients in (15). It should be noted that, for some parameters, the overestimation in the transient uncertain dynamic problems is significant in high-sensitive areas, which will be presented next.


## 4.2. Uncertain Transient Response Simulation Using the CIM

In this subsection, we present the numerical results of several cases with different uncertain parameters by using the CIM. To provide a reference for the accuracy and efficiency, the results are compared with those obtained from the scanning method. The expansion order in the CIM is 3 and the interpolation number is 4. The sampling points are 1000 for the scanning method (SM). Three kinds of uncertainties are included here for different occasions in engineering, namely, uncertain damping, mass property, and support stiffness. A typical uncertain deviation coefficient of 10% is adopted for all cases in this paper. The uncertain intervals of parameters are obtained by taking symmetric percentage dispersion around the nominal values, which can be expressed as(16)
q
j
I
=
q
_
j
,
q
-
j
=
q
j
c
-
β
q
j
c
,
q
j
c
+
β
q
j
c
,
j
=
1,2
,
…
,
r
,where q
j is the jth member of the uncertain parameter vector q and β is the deviation coefficient. q
_
j, q
-
j, and q
j
c are the lower bound, upper bound, and nominal value, respectively. Figures 3–5 give the uncertain transient response ranges of the rotor with different interval parameters and the corresponding errors considering the scanning results for reference. From these figures, we can observe that the uncertain effects of interval parameters are quite significant. Moreover, the mechanisms are quite different as well. The uncertain mass unbalance may cause a linear deviation of the deterministic curve while the others produce a large-amplitude vibration band [35].Figure 3
(a) Transient response of the rotor with uncertain mass unbalance and (b) calculation error.
(a)
(b)Figure 4
(a) Transient response of the rotor with uncertain support stiffness 2 and (b) calculation error.
(a)
(b)Figure 5
(a) Transient response of the rotor with uncertain disk mass and (b) calculation error.
(a)
(b)We use Lenovo Thinkpad T430 to calculate the response of speed-varying rotor. The CPU is Intel Core i7-3520 M @ 2.9 GHz, RAM is 12.0 GB, and the OS is Windows 10. The average computing time for the CIM is 6.9 s for each uncertain case and 25 minutes for the scanning method, and the advantage will be more obvious in higher uncertain dimensions. We then can give credit to the efficiency of the CIM in solving uncertain dynamic problems. Figure3 shows that the CIM can provide satisfactory estimated bounds when the damping is uncertain comparing with the SM and the calculation error is small. In that case, we consider the CIM to be fairly accurate. That is true for parameters that do not influence the intrinsic dynamic characteristics such as critical speeds and modes of the rotor system and the disk imbalance eccentricity, for another example. However, for the sensitive parameters such as the bearing stiffness and disk mass property, the CIM causes relatively large overestimations as demonstrated in Figures 4 and 5. Multipeaks are observed which are not in accordance with facts. CIM even gives negative value estimations in Figure 5 due to the high sensitivity of the transient response to the disk mass.Though the uncertainty in viscous damping itself does not affect much of the calculation accuracy, in order to investigate the influence of damping magnitude on the overestimation effect, several cases are chosen where different damping is considered when the disk mass is uncertain given its high sensitivity. The response bounds calculated by both the CIM and the SM are listed in Table2 for different damping at a specific revolution speed 128 rad/s near the resonant peak. In Table 2, lower bound is abbreviated as LB and UB is short for upper bound.Table 2
Response bounds of the system at 128 rad/s via different methods (×10−3 m).
Viscous damping (N⋅s/m)
CIM
SM
Error
LB
UB
LB
UB
LB
UB
100
1.915
6.152
1.982
5.090
0.0678
-1.062
200
1.044
4.025
1.109
3.301
0.0647
-0.723
300
1.085
2.742
1.156
2.370
0.0706
-0.373
400
1.116
2.018
1.160
1.824
0.0436
-0.194
500
1.083
1.586
1.107
1.477
0.0234
-0.109It can be observed form Table2 that the magnitude of damping can influence the response range directly when the disk mass is uncertain. Increasing damping decreases the magnitude of the resonant peak, uncertain response range width, and the calculation error. The error is getting smaller as the damping increases. Another fact shown in Table 2 is that the errors in the upper bounds are greater than that of the lower ones and the lower bounds fell less than the upper bounds. The bounds may not be strictly accurate because of finite sampling times, so that in some occasions, the lower bounds in greater damping are higher than the ones in smaller damping, but the absolute value is minor.Note a general fact that the overestimations are mostly distributed in the resonant regions. The largest overestimation point leads to nearly 1 mm, which is unacceptable for engineering practice. It is stated before that the overestimation is introduced by the interval arithmetic in the CIM, which is the key of the high efficiency of the interval method. To reduce the overestimation of the transient response prediction in the sensitive regions while keeping much of the efficiency, further algorithms must be incorporated.


## 5. Hybrid Interval Procedure for Overestimation Control

### 5.1. Response Bounds Determination by Searching the Min/Max Values

The Chebyshev approximation is obtained at every time step in the numerical implementation. For some cases, the direct interval arithmetic in the CIM may lead to large unacceptable overestimation. Then, on the premise of the good accuracy of the Chebyshev approximation itself, we can search the min/max values of the approximation function to approximate the bounds of the transient response. The first-order derivative information is used to control the overestimation for ODEs in [29]. When the approximation function is monotonic, it gives tighter bounds in those time steps. However, its applicability and easiness to implementation for complex rotor system should be tested further. In [35], the bounds are obtained by the maximin theorem of continuous function on closed intervals for the finite element model. According to this scheme, we intend to control the overestimation to an acceptable degree by searching the min/max values when the complete approximation expression is constructed. To illustrate the method more clearly, the Chebyshev approximation of the transient response in this section is rewritten in power basis given that the Chebyshev series can be expressed in both the trigonometric and polynomials expressions. The three-term recurrence relation of the Chebyshev series in polynomial form is given by(17)
C
0
x
=
1
C
n
+
1
x
=
2
x
C
n
x
-
C
n
-
1
x
,
n
=
1,2
,
…
.By following the same pattern, the one-dimensional approximation can be expressed as(18)
P
n
x
=
1
λ
∑
j
=
1
λ
U
~
x
j
+
2
λ
∑
i
=
1
n
∑
j
=
1
λ
U
~
x
j
C
i
x
j
C
i
x
.The notations in (18) have similar meaning as in (8) and (10). x
j is the interpolation point and x
j
=
cos
⁡
(
θ
j
). According to the maximin theorem of continuous function on closed intervals, there exist maximum and minimum values of P
n
(
x
). The roots of the derivative equation P
n
′
(
x
)
=
0 can be obtained and then we can compare the extreme values P
n
(
x
r
o
o
t
s
) with the two end points P
n
(
1
) and P
n
(
-
1
). The bounds of the approximated transient response subject to one-dimensional uncertainty are defined as the maximum and minimum values of P
n
(
x
).When multiple uncertain parameters are considered, the multidimensional polynomial form Chebyshev approximation can also be derived:(19)
P
n
x
1
,
…
,
x
r
=
∑
k
1
=
0
n
⋯
∑
k
r
=
0
n
1
2
p
Γ
k
1
,
…
,
k
r
C
k
1
,
…
,
k
r
x
1
,
…
,
x
r
.The power form of C
k
1
,
…
,
k
r
(
x
1
,
…
,
x
r
) is the tensor product of single-dimensional Chebyshev series. For example, the two-dimensional Chebyshev components of C
(
x
1
,
x
2
) can be expressed as(20)
C
0,0
x
1
,
x
2
=
1
C
1,0
x
1
,
x
2
=
x
1
,
C
0,1
x
1
,
x
2
=
x
2
⋯
C
k
1
,
k
2
x
1
,
x
2
=
C
k
1
x
1
C
k
2
x
2
,
k
1
+
k
2
=
n
.Thus, when the multidimensional approximation expression is obtained, the response bounds can be found by searching the min/max of the multivariable function. Many approaches such as derivative methods, optimization methods, and the scanning methods can do the job since the target function is a simple polynomial expression with all variables defined at the standard interval. In this paper, we use the scanning approach to find the min/max values of the approximation function. The scanning method scatters equal stepped points within the standard interval for each dimension and then evaluates the polynomial at those points to find the min/max values. Unlike iterations for the original complex model, the scanning method here is not too computationally expensive due to simplicity of the function. Suppose N points are used for each uncertain dimension denoted as x
^
i
j
,
i
=
1,2
,
…
,
r
,
j
=
1,2
,
…
,
N, the estimated bounds are defined as(21)
U
_
t
k
+
1
=
P
:
min
⁡
P
n
x
^
1
,
…
,
x
^
r
U
¯
t
k
+
1
=
P
:
max
⁡
P
n
x
^
1
,
…
,
x
^
r
.


### 5.2. Hybrid Analysis for Uncertain Transient Rotor Dynamics

As stated previously, despite the high efficiency, the overestimation of the transient response range near the resonant peaks calculated by the CIM is large. The solution method expressed in Section4.1 gives tighter bounds but it is relatively more time-consuming than the former due to the derivative operation in the searching of main/max values, though more efficient than the traditional scanning method from an overall view. The overestimation is significant only near the resonant and attenuation areas while being weak in other speed intervals. Thus, we propose the hybrid procedure that in the nonsensitive ranges the transient response range is calculated by the CIM. In the high-sensitive regions, the derivative-based method is employed where the efficiency is sacrificed a little for accuracy.If the required calculation accuracy is that the overall error is less than0.0001m, then all the speed ranges, where the response bounds estimated by the CIM are further deviated from the accurate curves than ε, should be calculated by the procedure provided in Section 5.1. The speed range will be 115.7~150.4 rad/s if the uncertain support stiffness K
2 is considered with 10% deviation from the nominal value and 125.5~145.6 rad/s for the uncertain disk mass m. For engineering practice, if the response lower bounds in the attenuation area are not interested or the uncertain deviation coefficient is small, then the speed range can be shortened to about 10% intervals around the system resonant peak. It should be noted that the deterministic response U
~ at collocation interpolations is the same as long as the expansion order and the number of interpolation points are consistent in the two algorithms explained in Sections 4.1 and 5.1, which will save the computational efforts by reusing the response data. Figures 6 and 7 are the recalculated response ranges and the errors for uncertain stiffness K
2 and uncertain disk mass with the same uncertain coefficient of 10% using the hybrid procedure (HP). When the disk mass and stiffness K
2 both are uncertain, the estimated response curves are given in Figure 8 in which for each parameter 10 scanning points are used. To avoid huge computation cost, the sampling points are cut to 200 for each uncertain dimension in the SM.Figure 6
(a) Response ranges with uncertain support stiffness 2 calculated by the hybrid procedure and (b) errors.
(a)
(b)Figure 7
(a) Response ranges with uncertain disk mass calculated by the hybrid procedure and errors.
(a)
(b)Figure 8
(a) Response ranges with multiple interval parameters calculated by the hybrid procedure and errors.
(a)
(b)As we can see, by using the hybrid procedure, the overestimation of the estimated response bounds is well controlled in Figures6–8 for situations of both single and multiple sources uncertainties. Compared with the results of the CIM, the errors here are small for most cases though there are some fluctuations in the attenuation areas. Generally, the hybrid procedure gives the tight bounds close to the scanning method. The computation time for hybrid procedure in the single-dimensional problem is ten minutes. In the multidimensional case, the total computation amount is enormous in the scanning method which takes many hours while for the hybrid analysis procedure it is about half an hour. From a practical point of view, the computation burden is acceptable given the accuracy.


## 5.1. Response Bounds Determination by Searching the Min/Max Values

The Chebyshev approximation is obtained at every time step in the numerical implementation. For some cases, the direct interval arithmetic in the CIM may lead to large unacceptable overestimation. Then, on the premise of the good accuracy of the Chebyshev approximation itself, we can search the min/max values of the approximation function to approximate the bounds of the transient response. The first-order derivative information is used to control the overestimation for ODEs in [29]. When the approximation function is monotonic, it gives tighter bounds in those time steps. However, its applicability and easiness to implementation for complex rotor system should be tested further. In [35], the bounds are obtained by the maximin theorem of continuous function on closed intervals for the finite element model. According to this scheme, we intend to control the overestimation to an acceptable degree by searching the min/max values when the complete approximation expression is constructed. To illustrate the method more clearly, the Chebyshev approximation of the transient response in this section is rewritten in power basis given that the Chebyshev series can be expressed in both the trigonometric and polynomials expressions. The three-term recurrence relation of the Chebyshev series in polynomial form is given by(17)
C
0
x
=
1
C
n
+
1
x
=
2
x
C
n
x
-
C
n
-
1
x
,
n
=
1,2
,
…
.By following the same pattern, the one-dimensional approximation can be expressed as(18)
P
n
x
=
1
λ
∑
j
=
1
λ
U
~
x
j
+
2
λ
∑
i
=
1
n
∑
j
=
1
λ
U
~
x
j
C
i
x
j
C
i
x
.The notations in (18) have similar meaning as in (8) and (10). x
j is the interpolation point and x
j
=
cos
⁡
(
θ
j
). According to the maximin theorem of continuous function on closed intervals, there exist maximum and minimum values of P
n
(
x
). The roots of the derivative equation P
n
′
(
x
)
=
0 can be obtained and then we can compare the extreme values P
n
(
x
r
o
o
t
s
) with the two end points P
n
(
1
) and P
n
(
-
1
). The bounds of the approximated transient response subject to one-dimensional uncertainty are defined as the maximum and minimum values of P
n
(
x
).When multiple uncertain parameters are considered, the multidimensional polynomial form Chebyshev approximation can also be derived:(19)
P
n
x
1
,
…
,
x
r
=
∑
k
1
=
0
n
⋯
∑
k
r
=
0
n
1
2
p
Γ
k
1
,
…
,
k
r
C
k
1
,
…
,
k
r
x
1
,
…
,
x
r
.The power form of C
k
1
,
…
,
k
r
(
x
1
,
…
,
x
r
) is the tensor product of single-dimensional Chebyshev series. For example, the two-dimensional Chebyshev components of C
(
x
1
,
x
2
) can be expressed as(20)
C
0,0
x
1
,
x
2
=
1
C
1,0
x
1
,
x
2
=
x
1
,
C
0,1
x
1
,
x
2
=
x
2
⋯
C
k
1
,
k
2
x
1
,
x
2
=
C
k
1
x
1
C
k
2
x
2
,
k
1
+
k
2
=
n
.Thus, when the multidimensional approximation expression is obtained, the response bounds can be found by searching the min/max of the multivariable function. Many approaches such as derivative methods, optimization methods, and the scanning methods can do the job since the target function is a simple polynomial expression with all variables defined at the standard interval. In this paper, we use the scanning approach to find the min/max values of the approximation function. The scanning method scatters equal stepped points within the standard interval for each dimension and then evaluates the polynomial at those points to find the min/max values. Unlike iterations for the original complex model, the scanning method here is not too computationally expensive due to simplicity of the function. Suppose N points are used for each uncertain dimension denoted as x
^
i
j
,
i
=
1,2
,
…
,
r
,
j
=
1,2
,
…
,
N, the estimated bounds are defined as(21)
U
_
t
k
+
1
=
P
:
min
⁡
P
n
x
^
1
,
…
,
x
^
r
U
¯
t
k
+
1
=
P
:
max
⁡
P
n
x
^
1
,
…
,
x
^
r
.


## 5.2. Hybrid Analysis for Uncertain Transient Rotor Dynamics

As stated previously, despite the high efficiency, the overestimation of the transient response range near the resonant peaks calculated by the CIM is large. The solution method expressed in Section4.1 gives tighter bounds but it is relatively more time-consuming than the former due to the derivative operation in the searching of main/max values, though more efficient than the traditional scanning method from an overall view. The overestimation is significant only near the resonant and attenuation areas while being weak in other speed intervals. Thus, we propose the hybrid procedure that in the nonsensitive ranges the transient response range is calculated by the CIM. In the high-sensitive regions, the derivative-based method is employed where the efficiency is sacrificed a little for accuracy.If the required calculation accuracy is that the overall error is less than0.0001m, then all the speed ranges, where the response bounds estimated by the CIM are further deviated from the accurate curves than ε, should be calculated by the procedure provided in Section 5.1. The speed range will be 115.7~150.4 rad/s if the uncertain support stiffness K
2 is considered with 10% deviation from the nominal value and 125.5~145.6 rad/s for the uncertain disk mass m. For engineering practice, if the response lower bounds in the attenuation area are not interested or the uncertain deviation coefficient is small, then the speed range can be shortened to about 10% intervals around the system resonant peak. It should be noted that the deterministic response U
~ at collocation interpolations is the same as long as the expansion order and the number of interpolation points are consistent in the two algorithms explained in Sections 4.1 and 5.1, which will save the computational efforts by reusing the response data. Figures 6 and 7 are the recalculated response ranges and the errors for uncertain stiffness K
2 and uncertain disk mass with the same uncertain coefficient of 10% using the hybrid procedure (HP). When the disk mass and stiffness K
2 both are uncertain, the estimated response curves are given in Figure 8 in which for each parameter 10 scanning points are used. To avoid huge computation cost, the sampling points are cut to 200 for each uncertain dimension in the SM.Figure 6
(a) Response ranges with uncertain support stiffness 2 calculated by the hybrid procedure and (b) errors.
(a)
(b)Figure 7
(a) Response ranges with uncertain disk mass calculated by the hybrid procedure and errors.
(a)
(b)Figure 8
(a) Response ranges with multiple interval parameters calculated by the hybrid procedure and errors.
(a)
(b)As we can see, by using the hybrid procedure, the overestimation of the estimated response bounds is well controlled in Figures6–8 for situations of both single and multiple sources uncertainties. Compared with the results of the CIM, the errors here are small for most cases though there are some fluctuations in the attenuation areas. Generally, the hybrid procedure gives the tight bounds close to the scanning method. The computation time for hybrid procedure in the single-dimensional problem is ten minutes. In the multidimensional case, the total computation amount is enormous in the scanning method which takes many hours while for the hybrid analysis procedure it is about half an hour. From a practical point of view, the computation burden is acceptable given the accuracy.


## 6. Conclusion

Uncertain dynamic response of a speed-varying rotor system is analyzed in this study considering interval uncertainties. The efficient nonintrusive Chebyshev interval method is employed for global uncertain response calculation. A derivative-based interval approach is proposed in the high-sensitive speed ranges where the overestimation is significant. By comparing with the scanning method, the accuracy of the proposed computation scheme is verified. The hybrid analysis procedure can reduce the errors significantly and the computation burden is acceptable. It can be easily adapted to similar uncertain transient dynamic problems due to the simplicity of mathematical formulation and the transparency in calculation.



---

*Source: 5904724-2018-07-18.xml*