Dynamic Response of Dual-Disk Rotor System with Uncertainties Based on Chebyshev Convex Method

Abstract

In this paper, a non-probabilistic Chebyshev convex method (CCM) for the transient dynamics of a dual-disk rotor system with uncertain parameters is proposed. The dynamic equation of the dual-disk rotor system is derived by the finite element method (FEM), and the deterministic response is obtained. Then the CCM is used to obtain the uncertain transient responses of the dual-disk rotor system. The amplitude ranges and response characteristics are quite close to the calculation results of the convex Monte Carlo simulation (CMCS), which verifies the accuracy and validity of the CCM. The experimental results demonstrate that the proposed method has good performance in the uncertainty analysis of the rotor system. More importantly, the results will be helpful to understand the dynamic behavior of dual-disk rotor systems with uncertainties and provide guidance for robust design and analysis.

1. Introduction

The rotor system is the core component of rotating machinery. To optimize and improve the dynamic performance, the dual-disk rotor system is widely employed in modern aero-engines [1,2]. It is more complicated than single-rotor systems in both the structural and dynamical regimes [3,4]. The analysis of traditional rotor dynamics is mostly based on deterministic models, which means that the parameters of the analysis object remain consistent. However, the physical parameters of the practical engineering mechanical always behave inherently in an uncertain way [5,6]. Component machining errors, installation gaps, uneven material distribution and fatigue during service are all typical sources of uncertainty [7,8,9]. Dynamic responses of rotor systems must be non-deterministic due to the uncertainty. Therefore, it is necessary to consider the influence of uncertain factors in the design, development, production, and use of rotor systems.

In fact, researchers have paid attention to the uncertainty quantifications in the field of rotor dynamics. Zhu [10] analyzes the stress-strain response induced by probabilistic plasticity with uncertainty. Giannella [11] proposes a stochastic approach, which takes multiple sources of uncertainty into account, to predict the probability distribution of the residual fatigue life of a defective railway axle. Sinou [12] applies the Kriging element model to predict the uncertain dynamic behaviors of flexible rotor system. Garoli and Castro [13] consider the influence of temperature uncertainty on the bearing dynamic coefficient, and use polynomial chaos expansion to solve the oil film whirl instability problem of rotor-bearing system. The polynomial chaos expansion and the harmonic balance method is used to quantify the effects of different random parametric uncertainties on the linear and non-linear dynamical characteristics [14]. The stochastic finite element method and the random variable model are adopted by Didier [15] to obtain the statistical characteristics of the unbalanced response of the flexible rotor. The non-parametric modeling technology is discussed by Murthy [16] on the unbalanced response and dynamic balance of the flexible rotor.

The probability density distribution of uncertainty is important to the dynamic research of rotor system with uncertainty. Whereas, for aero-engine rotor systems it is often difficult to obtain an accurate probability distribution of system uncertainty. In order to avoid the inability, scholars have proposed the concept of cognitive uncertainty and establish a variable probability model [17]. A series of methods such as evidence theory, interval theory, and convex model are established by Jiang [18,19,20] and Qiu [21]. To quantify the uncertainty of dynamic systems with non-probabilistic interval numbers, Wu [22] proposes Chebyshev inclusion function method by interval theory. The improved nonlinear Chebyshev method is applied by Ma [23,24] to carry out theoretical and experimental analyses on the bolted rotor system. Yang [25] proposes a stable non-embedded interval method based on orthogonal polynomial approximation theory and precise integration method. Transient dynamics on the accelerating start process of the rotor system with interval uncertainty is studied, and vibration characteristics with cognitive uncertainties is obtained. Fu [26] investigates the non-probabilistic steady-state dynamics of a dual-rotor system with parametric uncertainties under two-frequency excitations. However, the problem of uncertain rotor dynamics remains to be further studied, especially in the case of small samples. In addition, the calculation burden also needs to be reduced.

This paper introduces a convex model to describe the correlated interval variables of a dual-disk rotor system. And a non-probabilistic Chebyshev convex method is applied to obtain dynamic responses of dual-disk rotor system. The rest of the content is organized as follows. First, the finite element method for dual-disk rotor model and deterministic equations of motion will be described in Section 2. Then, in Section 3, the Chebyshev convex method (CCM) for dynamic analysis of uncertain rotors is explained. In Section 4, the dynamic responses with uncertainty are achieved by CCM as well as experimental verification is carried out. Finally, the concluding remarks are drawn in Section 5.

2. Dynamic Equation of Motion of the Rotor System

Finite element method (FEM) [27] is an analytical method that simulates real system by mathematical approximation. At present, FEM has been widely used in the field of rotor dynamics. It plays an important role in solving the critical speed, unbalanced response, and harmonic response analysis of the rotor system. And it is also widely used in the design and evaluation of its structural integrity [28,29]. The basic idea of FEM is discretizing continuous elastic body into finite numbers of interactive elements. When the finite element model converges, the approximate solution will converge to the exact solution. In the finite element model, shape function is used to approximate the unknown field function in each element, adjacent elements are connected to each other by nodes, and force is also transferred from one element to another by the nodes interconnecting them, and deformation relationship between elements is determined by deformation coordination relationship. By establishing the mechanical relationship between nodal force and nodal displacement of elements, a group of equations with nodal displacement as an unknown variable is obtained, thereby the node displacement component can be obtained by solving the equations.

A typical rotor system is usually composed of some discrete discs, shaft sections, and bearing supports. When building the finite element model, nodes are generated at the center of the disc, the position of the bearing, the center of the journal, and some positions on the elastic shaft. Then the rotor system is divided into several elements along the axis.

2.1. Modeling the Elastic Shaft Element

The main response mode of rotor system is bending vibration. Hence, the influence of the axial force and torsional vibration of the elastic shaft is not considered. The elastic shaft can be divided into multiple Euler-Bernoulli beam elements that do not consider the influence of the shear deformation. Figure 1 shows an elastic shaft element. The nodal displacement vector can be defined as $\left\{U_s\right\}={\left\{x_A,-{\theta }_{xA},y_A,{\theta }_{yA},x_B,-{\theta }_{xB},y_B,{\theta }_{yB}\right\}}^{\mathrm{T}}$ in fixed coordinates frame.

Based on Lagrange equation, the rotor motion equation can be established by FEM. The differential equation motion of elastic shaft element is

$$\left[M_s\right]\left\{{\ddot{U}}_s\right\}+\mathsf{\Omega }\left[J_s\right]\left\{{\dot{U}}_s\right\}+\left[K_s\right]\left\{U_s\right\}=\left[F_s\right]$$

(1)

where, [M_s] is the uniform mass matrix, Ω is the rotate speed of the shaft element, [J_s] is the rotary matrix and [K_s] is the rigidity matrix, [F_s] is the generalized force matrix.

It should be noted that the different coordinate order will cause different [M_s], [J_s] and [K_s], but the same result of U can be obtained after transformation.

2.2. Modeling the Rigid Disk

When the disc on the shaft only affects inertia characteristics but does not affect stiffness characteristics of the element, it can be treated as a rigid disk or a concentrated mass element.

Assuming that the axis of the disk coincides with the center of gravity, the nodal displacement vector of the disk element in generalized coordinate is

$$\left\{U_d\right\}=\left\{x,{\theta }_x,y,{\theta }_y\right\}$$

(2)

Based on Lagrange equation, the differential equation motion of rigid disk is

$$\left[M_d\right]\left\{{\ddot{U}}_d\right\}+\mathsf{\Omega }\left[J_p\right]\left\{{\dot{U}}_d\right\}=\left\{F_d\right\}$$

(3)

where, [M_d] is mass matrices for rigid disks, [J_p] is the rotary matrix, and {F_d} is the corresponding generalized concentrated force on the generalized concentrated mass element.

2.3. Mathematical Model of Flexible Rotor System

A rotor system structure diagram is shown in Figure 2. It consists of a shaft, two elastic bearings, and several disks. According to the FEM, the high-speed flexible rotor system is divided along the axial direction into n−1 shaft elements, k disks, and two bearings.

The rotor system with n nodes connected by n − 1 finite elements. The displacement vector of the system is

$$\left.\begin{array}{c}\left\{U_1\right\}={\left[x_1,-{\theta }_{x1},x_2,-{\theta }_{x2},\cdots ,x_N,-{\theta }_{xN}\right]}^{\mathrm{T}}\\ \left\{U_2\right\}={\left[y_1,{\theta }_{y1},y_2,{\theta }_{y2},\cdots ,y_N,{\theta }_{yN}\right]}^{\mathrm{T}}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\end{array}\right\}$$

(4)

Establish the differential equation motion of rotor system by Equations (1) and (3), and shown as follows:

$$\left.\begin{array}{c}\left[M_1\right]\left\{{\ddot{U}}_1\right\}+\left(\left[C_1\right]+\mathsf{\Omega }\left[J_1\right]\right)\left\{{\dot{U}}_2\right\}+\left[K_1\right]\left\{U_1\right\}=\left\{F_1\right\}\\ \left[M_2\right]\left\{{\ddot{U}}_2\right\}-\left(\left[C_2\right]+\mathsf{\Omega }\left[J_2\right]\right)\left\{{\dot{U}}_1\right\}+\left[K_2\right]\left\{U_2\right\}=\left\{F_2\right\}\end{array}\right\}$$

(5)

When,

$$\left\{U\right\}=\left\{\begin{array}{c}{U}_1\\ U_2\end{array}\right\}$$

(6)

The differential equation of motion of the rotor system changes to

$$\left[M\right]\left\{\ddot{U}\right\}+\left(\left[C\right]+\mathsf{\Omega }\left[J\right]\right)\left\{\dot{U}\right\}+\left[K\right]\left\{U\right\}=\left\{F\right\}$$

(7)

where, [M] is the global mass matrix, [C] is the damping matrix, Ω[J] is the rotary matrix, and [K] is the rigidity matrix. They are symmetric banded sparse matrix with order of 4N × N. {F} is the generalized external force matrix, and is the column matrix with order of 4N × 1.

3. The Chebyshev Convex Method for Dynamic Analysis

In this section, the ellipsoidal convex model is used to quantify the uncertainty. And the Chebyshev expansion function is introduced to analyze the dynamic response of the rotor system. Convex model is different from traditional interval model. It does not require precise probability distribution models or fuzzy membership functions of parameters, but only the upper and lower boundary of uncertain parameters. Meanwhile, the correlation between uncertain parameters is considered by convex model. Compared with perturbation technique or Taylor expansion, it is not limited by small uncertainty. The overestimation caused by the wrapping effect is also controlled effectively [30]. Chebyshev expansion function is an approximate proxy model. It uses Chebyshev series expansion method to obtain the solutions of differential equations at interpolation points. It can be used to characterize the response function of the original system, and has advantages of short calculation time, good calculation efficiency, and high calculation accuracy. Combining the ellipsoidal convex set model and Chebyshev expansion function, the Chebyshev convex method (CCM) is proposed to analyze the multi-convex uncertain system response.

3.1. The Ellipsoid Convex Method for Uncertain Parameters of Rotor Systems

The ellipsoid model is the most extensive among the convex models with the advantages of rounded boundaries and compact uncertainty regions. The multi-ellipsoid model is applied to describe the finite-dimensional distribution of the convex model process.

The n-dimensional interval variable is usually denoted as $X^{\mathrm{I}}=\left[X_1^{\mathrm{I}},X_2^{\mathrm{I}},\cdots ,X_n^{\mathrm{I}}\right]$. X represents an uncertain parameter. In FEM model, it is the value in the [M], [C], [K] and [J] of Equation (7). The uncertain region of $X^{\mathrm{I}}$ is expressed as a region contained by an n-dimensional hyper-ellipsoid as follows:

$$E_{X^{\mathrm{I}},W_e}=\left\{\left.X\right|{\left(X-X^m\right)}^{\mathrm{T}}{W}_e\left(X-X^m\right)\le 1\right\}$$

(8)

The region contained by $E_{X^{\mathrm{I}},W_e}$ is an n-dimensional ellipsoid model, assuming that all possible values of uncertain parameters are contained in ellipsoid. X is the ellipsoid convex variable, X^m is the median value of a convex variable, and W_e is the characteristic matrix of ellipsoid model. It determines the size and shape of the ellipsoid model and can represent the correlation between variables.

$$W_e=\left[\begin{array}{cccc}{W}_{11}& W_{12}& \cdots & W_{1n}\\ W_{21}& W_{22}& \cdots & W_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ W_{n1}& W_{n2}& \cdots & W_{nn}\end{array}\right]$$

(9)

When the axis of the ellipsoid is along the direction of the coordinate axis, W_e becomes a diagonal matrix.

As an example, the geometric characteristics of the ellipsoid model are studied, taking the two-dimensional ellipse in Figure 3. The geometric quantity is required to determine an ellipse: elliptic central coordinates $\left[X_1^m,X_2^m\right]$; length of semi-major axis $r_1$; length of short half axis $r_2$; and the angle between the major axis of the ellipse and the X₁ axis is donated as θ(−45° < θ < 45°). To facilitate the solution, the mean value, variance, and covariance of the ellipse are usually used to describe the geometric characteristics of the ellipse [31,32].

The ellipsoid convex model describes the uncertain parameters in the rotor system. The differential equation of the convex uncertain rotor system is

$$M\ddot{U}\left(\mathrm{x},\mathrm{t}\right)+\left(C+G\right)\dot{U}\left(\mathrm{x},\mathrm{t}\right)+KU\left(\mathrm{x},\mathrm{t}\right)=p\left(t\right),\hspace{1em}\hspace{1em}\mathrm{x}\in E_{X^{\mathrm{I}},W_e}$$

(10)

From the convex variable X, the response interval form of the uncertain rotor system is

$$U^I=\left[\underset{\_}{U},\overline{U}\right]$$

(11)

where

$$\begin{array}{l}\underset{\_}{U}=\min \left\{U|U=U\left(\mathrm{x}\right),\mathrm{x}\in E_{X^I,W_e}\right\}\\ \overline{U}=\max \left\{U|U=U\left(\mathrm{x}\right),\mathrm{x}\in E_{X^I,W_e}\right\}\end{array}$$

(12)

3.2. Chebyshev Method for Convex Analysis of the Rotor Dynamic Response

In this section, a Chebyshev interval expansion function based on the Chebyshev series expansion is proposed.

Considering functions in the ellipsoid as follows:

$$F=f\left(\mathrm{x}\right),\hspace{1em}\hspace{1em}{\mathrm{x}=[x_1,x_2,\dots ,x_n]}^{\mathrm{T}}\in E_{{\mathrm{X}}^I,W_e}=\left\{\mathrm{x}\left|{\left(\mathrm{x}-{\mathrm{x}}^c\right)}^{\mathrm{T}}{W}_x\left(\mathrm{x}-{\mathrm{x}}^c\right)\le 1\right.\right\}$$

(13)

Function f (x) can be approximated by Chebyshev polynomials.

$$f\left(\mathrm{x}\right)\approx p_p\left(\mathrm{x}\right)=\frac{1}{2}{f}_0+{\displaystyle \sum _{i=1}^k{f}_i{C}_i\left(\mathrm{x}\right)},\hspace{1em}\hspace{1em}\mathrm{x}\in E_{X^I,W_e}=\left\{\mathrm{x}\left|{\left(\mathrm{x}-{\mathrm{x}}^c\right)}^{\mathrm{T}}{W}_e\left(\mathrm{x}-{\mathrm{x}}^c\right)\right.\le 1\right\}$$

(14)

where, P is the order of Chebyshev polynomials, f_i is the constant coefficient of Chebyshev polynomials C_i.

For an n-dimensional problem, the Chebyshev polynomials are defined as

$$C_{k_1,k_2,\cdots ,k_n}\left(x_1,x_2,\cdots ,x_n\right)=\cos \left(k_1{\theta }_1\right)\cos \left(k_2{\theta }_2\right)\cdots \cos \left(k_n{\theta }_n\right)$$

(15)

where, k_i (I = 1, 2, …, n) is the nonnegative integer, θ_i (I = 1, 2, …, n) is a function of the convex variable.

$${\theta }_i=\arccos \left(\frac{2x_i-\left(x_i^R+x_i^L\right)}{x_i^R-x_i^L}\right)\in \left[0,\pi \right],\hspace{1em}\hspace{1em}i=1,2,\dots ,n$$

(16)

The truncated Chebyshev polynomial of function f(x) can be rewritten as

$$\begin{array}{l}f\left(x\right)\approx p_p\left(x\right)=\frac{1}{2}{f}_0+{\displaystyle \sum _{k_i=0}^k\cdots }{\displaystyle \sum _{k_n=0}^k{f}_{k_1,\dots ,k_n}{C}_{k_1,\dots ,k_n}\left(\mathrm{x}\right)}\\ \hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}=\frac{1}{2}{f}_0+{\displaystyle \sum _{k_i=0}^k\cdots }{\displaystyle \sum _{k_n=0}^k{f}_{k_1,\dots ,k_n}\cos i\left[{\theta }_i\right]},\hspace{1em}\theta \in \left[0,\pi \right]\end{array}$$

(17)

According to the Mehler integral method, f_i can be calculated by Equation (18).

$$f_i=\frac{2}{\pi }{\displaystyle {\int }_{-1}^1\frac{1}{\sqrt{1-x^2}}f\left(x\right)C_i\left(x\right)dx}\approx \frac{2}{m}{\displaystyle \sum _{j=1}^{m}f\left(x_j\right)C_i\left(x_j\right)=}\frac{2}{m}{\displaystyle \sum _{j=1}^{m}f\left(\cos {\theta }_j\right)\cos i{\theta }_j}$$

(18)

Coefficient $f_{k_1,\dots ,k_n}$ is expressed as

$$\begin{array}{l}{f}_{k_1,\dots ,k_n}={\left(\frac{2}{\pi }\right)}^n{\displaystyle {\int }_{-1}^1\cdots }{\displaystyle {\int }_{-1}^1\frac{f\left(\mathrm{x}\right)C_{k_1,\dots ,k_n}\left(\mathrm{x}\right)}{\sqrt{1-x_1^2}\cdot \dots \cdot \sqrt{1-x_1^2}}d{x}_1\dots d{x}_n}\\ \hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\approx {\left(\frac{2}{m}\right)}^n{\displaystyle \sum _{j_1=1}^m\cdots }{\displaystyle \sum _{j_n=1}^{m}f\left(\cos {\theta }_{j_1},\dots ,\cos {\theta }_{j_n}\right)\cos \left(k_1{\theta }_{j_1}\right)\dots }\cos \left(k_n{\theta }_{j_n}\right)\end{array}$$

(19)

$$\begin{array}{l}{\theta }_{j_i}=\frac{2j_i-1}{m}\cdot \frac{\pi }{2},\hspace{1em}{j}_i=1,2,\dots ,m,\hspace{1em}i=1,2,\dots ,n\\ x_{i,j_i}=\frac{x_i^R+x_i^L}{2}+\frac{x_i^R-x_i^L}{2}\cos {\theta }_{j_i},\hspace{1em}{j}_i=1,2,\dots ,m,\hspace{1em}i=1,2,\dots ,n\end{array}$$

(20)

where, m is the number of sampling points that need interpolation in each dimension, l denotes the total number of zeros in subscript $k_1,\dots ,k_n$.

Chebyshev polynomial surrogate model of n-dimensional function f(x) are defined as

$$\begin{array}{c}f\left(x\right)=f\left(x_1,x_2,\cdots ,x_n\right)\approx p_p\left(x\right)={\displaystyle \sum _{k_1=0}^p\cdots }{\displaystyle \sum _{k_n=0}^p{\left(\frac{1}{2}\right)}^l}{f}_{k_1,k_2,\cdots ,k_n}{C}_{k_1,k_2,\cdots ,k_n}\left(x\right),\\ x\in E_{X^{\mathrm{I}},W_e}=\left\{\left.x\right|{\left(x-x^m\right)}^{\mathrm{T}}{W}_e\left(x-x^m\right)\right\}\end{array}$$

(21)

The differential Equation (10) of the convex uncertain rotor system can be solved based on Equation (21).

$$U=U\left(x\right)\approx p_p\left(x\right)={\displaystyle \sum _{k_1=0}^p\cdots }{\displaystyle \sum _{k_n=0}^p{\left(\frac{1}{2}\right)}^l{U}_{k_1,\dots ,k_n}{C}_{k_1,\dots ,k_n}\left(x\right)},\hspace{1em}x\in E_{X^I,W_e}$$

(22)

The extreme values of the response U can be found as follows:

$$\begin{array}{l}{U}^L=\underset{\mathrm{x}\in E_{X^{\mathrm{I}},W_e}}{\min }\left\{U\left(x\right)\right\}\approx \underset{\mathrm{x}\in E_{X^{\mathrm{I}},W_e}}{\min }\left\{p_p\left(x\right)\right\}\\ U^R=\underset{\mathrm{x}\in E_{X^{\mathrm{I}},W_e}}{\max }\left\{U\left(x\right)\right\}\approx \underset{\mathrm{x}\in E_{X^{\mathrm{I}},W_e}}{\max }\left\{p_p\left(x\right)\right\}\end{array}$$

(23)

The computational flowchart of the framework is demonstrated in Figure 4.

4. Numerical Analysis of the Rotor Dynamic Response Based on CCM

In this section, some simulations are carried out taking the high-speed flexible rotor as an example. The characteristic of the model is introduced and the deterministic response of the rotor system will be given first. To verify the reliability and stability of the CCM, the transient response of rotor system with correlated uncertain parameters is calculated and compared with the results calculated by CMCS. After that, experiment is conducted to verify the correctness of the simulation results. According to the results of numerical simulation and experiment, some preliminary discussions are made.

4.1. Numerical Simulation Model and Deterministic Response

As shown in Figure 5, the rotor system is modeled by beam element, mass element, and spring element with convex variable parameters. Detailed deterministic geometric and physical parameters of the rotor system studied are given in

Table 1. Values of physical parameters of the dual-disk rotor.

Symbol	Description	Value
E (GPa)	Young’s modulus	210
R (m)	Radius of the shaft	0.005
r (m)	Radius of disk 1, 2	0.0375
L (m)	Length of shaft	0.55
L1, L2, L3, L4, L5 (m)	Length of shaft 1, 2, 3, 4, 5	0.10, 0.10, 0.15, 0.10, 0.10
K1, K2 (N/m)	Stiffnesses of supports 1, 2	1.20 × 105, 1.0 × 105
C (N·s/m)	Damping coefficient	32
m1, m2 (kg)	Mass of disk 1, 2	0.483, 0.481
ρ (kg/m3)	Density	7800
Ip1, Ip2 (kg·m2)	Polar moment of inertia	3.242 × 10−4, 3.228 × 10−4
e1,e2 (m)	Eccentricity of disk 1, 2	4 × 10−5, 4 × 10−5

.

The deterministic response is illustrated to characterize the dual-disk rotor system without uncertainty and to give an initial impression of its dynamic behaviors. Based on the deterministic model and numerical simulations, the first three critical speeds of the rotor are shown in the

Table 2. Calculated value of the first three critical speeds of high−speed flexible rotors.

Critical Speeds /(rpm)
first-order	second-order	third-order
2838.67	6406.81	9985.07

and the transient-state vibration amplitude time when the two disks rotate at the same time, history of the disk geometric center is shown in Figure 6.

By observing the value of the first three critical speeds and deterministic response of the dual-disk rotor, it can be found that the rotor system has two critical speeds in the working speed range, and the first two peaks appear at 2838.67 rpm and 6406.81 rpm for both of the disks, respectively. Meanwhile, the second-order peak of disk 1 is significantly greater than that of disk 2. This phenomenon is analyzed by researchers and could be attributed to the different stiffness of bearings [33,34,35]. Moreover, the response law of the two disks is basically similar. The vibration state at the critical speed satisfies the actual operating law of the rotor system. In order to simplify the analysis process, the subsequent uncertainty analysis only analyzes disk 2.

4.2. Single Parameter Uncertainty

When there is a single uncertain parameter, the convex model evolves into an interval model. By the Chebyshev Convex Method (CCM), the uncertainty caused by the uncertain coefficient can be acquired. Several uncertain parameters are considered below, and the response boundary of the critical speeds calculated by CCM are compared with those by CMCS. In order to compare the influence of different single parameters of the rotor system, parameters are defined in same uncertain degree. According to the actual needs of the project, the degree of uncertainty is set to 10%.

Firstly, the damping is viewed as an uncertain parameter, and the degree of uncertainty is 10%. Using CCM, the uncertain dynamic response affected by uncertain parameters can be obtained, and the results are shown in Figure 7.

Compared with CMCS, the response boundary in the critical speed range calculated by CCM is similar, and the deterministic transient response curve of the system is symmetrically shifted up and down. The resonance peak position does not change and falls within the critical speed range.

Generally, it is difficult to define the accurate value of the stiffness of a support. In this case, the stiffness of bearing 2 is taken as a convex variable to cover its variability and the uncertain degree of support stiffness K₂ is 10%. The response range of the center of disk 2 is shown in Figure 8.

As can been seen from Figure 8, we can find that there is a peak shift in second-order resonance peaks under the influence of the support stiffness uncertainty at no such phenomenon in the first-order. This means the different sensitivities of the system in different speed ranges to the first two critical speeds. Meanwhile, the upper and lower boundaries of the deflection amplitude increase at the peak of the second-order critical speed, which suggests that the inherent properties of the high-speed flexible rotor system are sensitive to the support stiffness. In addition, due to the simplification of the model and the calculation method in the theoretical analysis, the theoretical modeling cannot fully simulate the actual structure of the original high-speed flexible rotor system, resulting in small errors and overestimation of response.

In engineering practice, the rotor system may be affected by manufacturing error, material degradation, and wear. The actual mass often deviates from the theoretical value after long-term work. Therefore, the uncertainty of mass should be considered in the analysis. To consider the uncertainty of mass, the uncertain degree of mass of disk 2 is taken to be 10%. The ellipsoidal convex model is used to solve the analysis interval response, and the results are shown in Figure 9.

The uncertainty of rotor disk mass leads to the deviation of system response in the whole speed range as show in Figure 9. More specifically, compared with the deterministic peak, the first peak of the uncertainty response shifts to the left and the second peak shifts to the right. That means the double rotors are very sensitive to the disk mass, which can be used as a key factor in the design and maintenance of such engineering systems. The results of the CCM also show that the upper and lower boundary gauges completely encapsulate the deterministic response values.

Figure 7, Figure 8 and Figure 9 show that the accuracy of CCM results is related to the sensitivity of critical speed to parameters. In the whole response range, the interval response results calculated by CCM are stable. Compared with the deterministic response results, the actual operating conditions of the rotor system can be well reflected. The resonance peaks under different uncertain parameters are calculated and recorded in

Table 3. Calculation of peaks under different uncertainty parameters.

Uncertainty and Nominal Peaks		Peaks(m)
		First-Order [lower, upper]	Second-Order [lower, upper]
Nominal value uncertainty 10%		1.740 × 10−4	5.500 × 10−5
	C	[1.570 × 10−4,1.936 × 10−4]	[4.654 × 10−5, 6.495 × 10−5]
	K2	[1.718 × 10−4, 1.768 × 10−4]	[4.360 × 10−5, 6.601 × 10−5]
	m2	[1.678 × 10−4, 1.805 × 10−4]	[5.133 × 10−5, 6.063 × 10−5]
		Maximum error in interval
	C	11.3%	18.1%
	K2	1.7%	20.7%
	m2	3.7%	10.2%

. Comparing the error caused by the parameter uncertainty at the peak of the first-order critical speed and the second-order critical speed under the same uncertainty, it can be seen that the error at the second order is significantly greater than that at the first order. This demonstrates that the peaks in the second-order critical speed range are more sensitive to parameters than the first-order peaks. And compared with the results of CMCS, CCM can have a good interval wrapping effect. Meanwhile, in order to verify the computational efficiency of CCM, the calculation of the mass uncertainty response results is taken as an example, and the computer time consumption of the two methods at the same calculation environment is given in

Table 4. Calculation time of transient response of disk 2 with different methods.

	CMCS (50 Samples)	CMCS (100 Samples)	CMCS (500 Samples)	CCM (n = 5, p = 8)
calculation time (min)	3.26	6.01	32.84	0.57

. Obviously, it can be seen that the calculation time of Chebyshev method based on CCM is much shorter than that of CMCS and the computing efficiency is significantly improved.

4.3. Multi-Parameters Uncertainty

This subsection pays attention to the influences of multi uncertain parameters [33] on the dynamic behaviors of the dual-disk rotor. Consider the uncertainties in the supporting stiffness K₂, damping C and mass m₂ of disk 2 that studied in the previous subsections, in which the supporting stiffness K₂ and damping C are correlated uncertain parameters, and the disk mass m₂ is an independent parameter. The uncertain degrees are 10% for the uncertain parameters. The deflection curves of rotor system calculated by CCM and CMCS are shown in Figure 10.

Comparing the curves of Figure 7 and Figure 10 to Figure 9, at the same uncertainty degree, the overall envelope range caused by multi-parameters uncertainty is significantly wider than that caused by single parameter uncertainty. The peak shifts are observed at the critical speeds. It shows that the dynamic response is significantly affected by the multiple uncertain parameters.

4.4. Critical Speeds Uncertainty

This subsection pays attention to the critical speeds influenced by uncertainty of the rotor system. The elastic modulus E, material density ρ, and support stiffness K₂ are regarded as uncertain parameters. Among them, the elastic modulus E and material density ρ are related. The first three critical speeds influenced by different uncertain parameters are calculated when the uncertain degree is 10%. And the calculation results are recorded in

Table 5. Calculation of critical speeds under different uncertainty parameters.

Uncertainty and Nominal Critical Speeds		Critical Speeds (rpm)
		First-Order [Lower, Upper]	Second-Order [Lower, Upper]	Third-Order [Lower, Upper]
Nominal value uncertainty 10%		2838.67	6406.81	9985.07
	K2	[2796.70, 2873.55]	[6255.72, 6544.52]	[9894.02, 10,079.23]
	E	[2762.54, 2905.78]	[6369.19, 6437.51]	[9717.22, 10,246.43]
	ρ	[2800.36, 2864.39]	[6397.32, 6444.52]	[9901.79, 10,035.64]
	Multiple	[2752.25, 2918.21]	[6188.02, 6624.86]	[9444.11, 10,593.81]
		Maximum error in interval
	K2	1.48%	2.36%	0.94%
	E	2.68%	0.59%	2.68%
	ρ	1.30%	0.58%	0.89%
	Multiple	3.04%	3.41%	6.10%

.

It can be concluded from Table 5 that the critical speeds of the rotor system are no longer fixed values, but change within an interval. Comparing the errors of critical speeds influenced by different uncertain parameters, it can be seen that errors effected by multi- parameters are greater than those caused by single parameters.

The first three critical speeds influenced by uncertain multi-parameters are calculated when the uncertain degrees are 10% and 15% respectively. And the calculation results are recorded in

Table 6. Calculation of critical speeds under different uncertainty.

Uncertainty and Nominal Critical Speeds		Critical Speeds
Multi-Parameters Uncertainty (E, ρ, K2)		First-Order (rpm) [Lower, Upper]	Second-Order (rpm) [Lower, Upper]	Third-Order (rpm) [Lower, Upper]
		2838.67	6406.81	9985.07
Uncertainty	10%	[2752.25, 2918.21]	[6188.02, 6624.86]	[9444.11, 10,593.81]
	15%	[2704.75, 2959.06]	[6071.80, 6739.12]	[9196.52, 10,909.02]
		Maximum error in interval
	10%	3.04%	3.41%	6.10%
	15%	4.72%	5.23%	9.25%

. Comparing the critical speed interval errors under different uncertainty, the error on the critical speeds increases accordingly with the increase of uncertainty.

4.5. Experimental for Uncertainty Response

The dual-disk rotor system shown in Figure 11 is used for the uncertainty response experimental, and the physical parameters of the rotor are tabulated in Table 1. The rotor system consists of a straight shaft and two discs overhung on the right. The shaft is fixed on the bracket by two bearing supports, and the left end of the shaft is connected with the drive motor by the coupling. The maximum speed of the motor is 10,000 r/min. The two displacement sensors (NO. 3 and 4 in Figure 11) are arranged horizontally on the right side of the bearing supports of the discs to measure the dynamic displacement of the rotor. And the rotor speed is measured by two displacement sensors (NO. 1 and 2 in Figure 11) which receive the displacement information on the coupling. The deflections of the discs are obtained by Hilbert transform under the condition of neglecting the anisotropy of rotor.

In order to verify the effectiveness of CCM, the uncertainty experiment of mass parameter is carried out. The simulation of mass uncertainty is realized by adding screws to the disk. The screws shown in Figure 12 with the weight of 4.73 g are inserted into symmetrical holes on the disk 2 to simulate the change of mass (4 screws correspond to the mass uncertainty of 4% and 8 screws correspond to the mass uncertainty of 8%). Due to the existence of uncertainty, the dynamic response of the system can be any possible value in the response interval. Compare the vibration response of disk 2 measured by the experiment with the numerical simulation results, and the results are shown in Figure 12.

Comparing the experimental results of dual-disk rotor system with the simulation results shown in Figure 13.

The red solid lines in the Figure 13 are the upper and lower bounds of the rotor system simulation results when the uncertainty is 15%. It can be seen that when the uncertainty is 4%, 8%, and no uncertainty, the peaks of the experimental deflection curves are surrounded by the simulation results. Meanwhile, the experimental results show that with the increase of mass uncertainty, the resonance peaks of the deflection curves shift to the left in the first critical speed range and shift to the right in the second critical speed range. The trend of simulation results is consistent with experimental result when the mass is uncertain. In addition, the slight frequency shift and interval overestimation can be observed at the second-order critical speed. This is because theoretical modeling cannot fully simulate the actual structure of the high-speed flexible rotor system, resulting in errors in response calculations. In general, the experimental rotor response is basically contained in the simulation uncertainty interval, which proves the reliability and stability of CCM.

5. Conclusions

The uncertain dynamics responses of a dual-disk rotor system with uncertainties are analyzed by Chebyshev convex method. The rotor motion differential equation is established by finite element method. An ellipsoid convex model is proposed to describe the uncertain parameters with correlation. And a non-probabilistic CCM is developed to obtain uncertain dynamic response of dual-disk rotor system. The critical speeds and dynamic responses are calculated by the CCM. Results show that the responses of the rotor system are no longer deterministic values, but interval values. And there are phenomena such as the left-right shift and up-down fluctuation of the resonance peak. It is caused by multi-parameters uncertainty. The calculation accuracy and efficiency are verified by CMCS. Compared with the CMCS, CCM computing time is much less. At the same time, CCM provides accurate upper and lower bounds. Moreover, experimental results of dual-disk rotor uncertainty verify the accuracy and efficiency of CCM. Future work will be carried out by the present authors on optimizing the CCM, which could provide effective help to the rotor dynamics designer.