Non-probabilistic analysis of a double-disk rotor system with uncertain parameters

Vibration is a major issue in rotor systems. Due to the presence of material property dispersions, manufacture or assembling errors and time-varying working status, rotor systems are always subject to uncertainties. The uncertainties should be taken into consideration to understand the dynamic characteristics of rotors more thoroughly. In this study, interval analysis is carried out to investigate the non-probabilistic characteristics of a double-disk rotor with uncertain parameters. The uncertainties are modeled as uncertain-but-bounded variables due to insufficient essential information to define their precise probabilistic distributions. The deterministic analysis model is derived by the finite element method (FEM). The accuracy and effectiveness of the proposed method in solving uncertain rotor problems are validated by comparative study with the Monte Carlo simulation (MCS). Several cases where different physical parameters are regarded as uncertain are investigated and the dynamic response bounds are obtained. Simulations suggest that uncertainties have significant influence on the dynamic characteristics of the rotor system. Multi-source uncertainties propagation can cause heavy vibrations in mechanical systems.


Introduction
Rotor system is a crucial part of rotating machinery such as aero-engines and electric generators.One of the fundamental and significant issues encountered in rotor systems is the heavy vibration aroused by the mass unbalance.The dynamic characteristics are frequently investigated by designers and researchers to understand the performance of a rotor system.Many efforts have been put in during the past few decades both theoretically and experimentally [1][2][3][4][5][6][7][8].As can be seen, these contributions were mostly focused on rotor systems in which the parameters and loads are all deterministic.The fact is that uncertainty is unavoidable in practice, which is especially true in rotor systems.Material property dispersions, varying working status and assembly or measurement errors are all possible sources of uncertainty.In uncertain problems, the traditional analysis methods fail to give reasonable estimations of the response ranges induced by uncertain inputs.Thus, new solution methods should be incorporated.
Pioneering works have been done to study rotordynamics involved with various kinds of uncertainties [9].Liao [10] proposed a global resonance optimization method and applied it to the uncertainty quantification of rotor systems.Didder and Sinou et al. [11,12] had carried out many researches concerning uncertain rotor systems with nonlinearity and faults in which the uncertainties were represented by the Polynomial Chaos Expansion (PCE).The PCE was also applied to analyze the performance of a H-Darrieus rotor under uncertainty in transient turbulent flow simulations recently [13].Sinou et al. [14] further studied the effect of the PCE order on the dynamic response of an asymmetric uncertain rotor using a recursive stochastic harmonic balance method.Gan et al. [15] studied the sensitivity of the first order critical speed of an uncertain Jeffcott rotor to uncertain parameters via the MCS.It was found that the uncertain effect cannot be neglected in some circumstances.The proper orthogonal decomposition method can be applied to model reduction for deterministic and random rotor systems [16].Li et al. [17] investigated the nonlinear stochastic response of an angular misalignment rotor considering the fluid induced random forces.Nonparametric stochastic modelling method for uncertain rotor systems in unbalance and balancing aspects was proposed by Murthy et al. [18].It should be noted that all the works above were done in the probabilistic frame, which means that the precise distributions of the uncertainties must be known beforehand.That is a requirement hard to meet due to lack of sufficient statistical information.For example, the epistemic uncertainty introduced by the designer's inadequate essential knowledge is generally difficult to define.In that case, the probabilistic methods are sometimes inapplicable.However, the varying ranges of uncertain parameters are usually easy to obtain, which can be described by uncertain-but-bounded variables.So non-probabilistic based methods are advocated in general uncertain problems and the interval analysis method is a typical one among them.A universal review of the non-probabilistic methods for the FEM analysis in uncertain problems was contributed by Moens et al. [19].Qiu et al. [20] had proposed many interval approaches including perturbation method, Taylor interval expansion method and interval collocation method to calculate static and dynamic problems in structures.Qi et al. [21] further compared these interval methods with probabilistic ones in the aspects of efficiency, accuracy and overestimations.Wu et al. [22] proposed a non-intrusive interval method for dynamic systems with uncertain parameters by introducing the Chebyshev inclusion function, which has been successfully applied to uncertain vehicle systems [23].It can achieve high efficiency and accuracy with low expansion order and control the overestimation effectively.So far, few researches can be found using interval methods for uncertain rotor problems.Ma et al. [24][25][26][27] studied the variations of the critical speeds, dynamic responses and modes in uncertain rotor systems using combination of perturbation method, the interval modal superposition method and interval Taylor expansion method.Fang et al. [28] applied the interval method to a simple rotor model for modal characteristic investigations.These methods were validated by different rotor systems with considerable accuracy.However, an important feature of them is that modifications need to be made to the original procedure.In other words, they are intrusive.
In this study, the non-probabilistic dynamic response analysis of a double-disk rotor with uncertain-but-bounded parameters is carried out based on interval analysis in a non-intrusive way.The FEM is used to derive the deterministic analysis model.In order to verify the accuracy and efficiency of the interval method in solving uncertain rotordynamics, the results are compared with those obtained by the MCS.Several simulations are done considering different uncertain parameters.Some conclusions are drawn based on the simulation results.

Rotor model and equations of motion
In this paper, a double-disk rotor system supported by two journal bearings is considered.The two disks are located in the middle of the elastic shaft and the layout of the system is shown in Fig. 1(a).The FEM is used to discretize the system into elements with two nodes each.The nodal displacement vector, which consists of two rotating angles and two lateral displacements, is defined by [        ] in fixed coordinates as demonstrated in Fig. 1(b) [11].The deterministic values of physical parameters of the rotor are given in Table 1.
The equations of motion of the rotor system can be obtained by assembling all of the shaft elements: where  ,  ,  and  are the global mass, damping, gyroscopic and stiffness matrices of the system, respectively. ,  and  represent the displacement, velocity and acceleration vectors, respectively.() is the excitation vector and  corresponds to the rotating speed.

Interval analysis of the uncertain rotor response
The interval analysis methods are widely used nowadays in uncertain dynamic investigations.Approximation theory and derivative information are all introduced into the calculation process [21,22].These kind of interval procedures are non-intrusive and capable of controlling the large overestimation compared with Taylor series based methods.High accuracy can be achieved via low expansion orders and a few interpolations.Moreover, it is convenient to obtain higher accuracy by higher order and more interpolation points while the computation effort is not increased enormously comparing with the MCS.
Firstly, we define an interval variable  in mathematic form: in which superscript  denotes an interval character,  and ̅ are the lower and upper bounds of the uncertain variable .The first class Chebyshev orthogonal polynomials can be described as: where  is a nonnegative integer.The one-dimensional Chebyshev approximation is defined as: where  are the expansion coefficients.They can be expressed as: in which () is the function being approximated and  () is the deterministic value of ().√1 −  is the weight function.The integration can be further calculated by the Gauss-Chebyshev integration formulas: with  being the number of interpolation points.In order to use the boundedness of trigonometric functions, the variable  can be replaced as: where  ∈ [, ̅ ].
As stated previously, we consider that the rotor system under study subject to uncertainties and then the parameters such as the support stiffness and density are all uncertain [24].No sufficient data is available to define their probabilistic distributions.Therefore, the uncertainties are modeled as uncertain-but-bounded parameters.For clarity, we define an uncertain parameter vector, which is composed of all the uncertain parameters in the rotor system: in which  is the total number of uncertain parameters.The th member of  can be expressed in interval form as: where  and  are the lower and upper bounds, respectively. is the mean value of the parameter and  represents the deviation coefficient.The matrices in the general equations of motion are all related to one or more members of the uncertain parameter vector, thus Eq. ( 1) can be rewritten as: () (, ) + (, ) (, ) + (, )(, ) = (, ).
The dynamic response calculation of the system can be interpreted as initial value problems of Eq. ( 10) under the constraints of Eq. ( 8) with respect to time :  (, ) = :  +  +  = , ( ) =  ,  ∈  , (11) in which  is the initial value of the response displacement vector,  represents the possible interval parameter vector comprised of individuals expressed by Eq. ( 9).The multi-dimensional Chebyshev orthogonal approximation of Eq. ( 11) is: where  is the number of zeros in subscripts,  is the expansion order. ⋯ ( ) represents the -dimensional Chebyshev orthogonal series which is defined as [22]: Γ ⋯ denotes the -dimensional expansion coefficient, which can be obtained by Mehler integration at Chebyshev interpolation points.The expression is: with  representing the Chebyshev interpolation point,  being the total interpolation point number which should be no less than  + 1 for the benefit of accuracy. (, cos , ⋯ , cos ) denotes the deterministic rotor response values calculated at interpolation points.They can be numerically solved by regular numerical integration method as deterministic ODEs.Generally, the interpolation points are the zeros of the Chebyshev series: The response interval of the double-disk rotor system is finally given by: Now, we have transformed the uncertain dynamic problem into a set of deterministic problems at the interpolation points and the structures of the governing motion equations remain unchanged compared with Eq. ( 1).The dynamic response range can then be explicitly calculated by interval arithmetic as denoted by Eq. ( 16).

Numerical results and discussions
The non-probabilistic dynamic response of the uncertain rotor system detailed in Section 2 is obtained based on interval arithmetic and the Newmark- numerical integration method which has been proved to be unconditionally stable as long as the integral parameters are chosen properly [29].To verify the effectiveness and accuracy of the interval method, results are compared with those calculated by the MCS [30].Then the uncertain responses in single or multi-dimensional uncertainties are investigated, respectively.

Deterministic dynamic response
Based on the deterministic parameters of the rotor model, the dynamic characteristic analysis of the double-disk rotor system is carried out.Fig. 2 illustrates the steady state responses of the rotor for deterministic parameters in a traditional sense.The first two critical speeds are 1920 revolutions per minute (r/min) and 5435 r/min.The resonant peaks of the two disks are given by Table 2.It can be observed from Fig. 2 that the vibration of the second resonance is higher than the first one which is especially true for disk 1.In all of the uncertain dynamic response figures presented in the following text, the deterministic steady state response will be plotted firstly for direct reference.

Method verification
In this subsection, we use the MCS to verify the accuracy and efficiency of the interval method in solving non-probabilistic uncertain dynamic problems.For simplicity and clarity, the dynamic response ranges under one-dimensional uncertainty at several rotating speeds are investigated.The support stiffness 1 is assumed to be uncertain which can be expressed in interval form as: where  and  are the mean value and the deviation coefficient of  .The response bounds under deviation coefficient 10 % obtained from the MCS and the interval method are given in Table 3.Total sampling times in the MCS are 100 and 500.The results calculated with different expansion order and interpolation points are also given in the table for comparisons.We refer the interested readers to Ref. [21] for the selection principles of the calculation parameters based on the truncation error of the Chebyshev approximation.To illustrate the computational efforts taken in the two approaches, the computing costs (CPU time) are listed in Table 4.In Tables 3 and 4, RS is short for rotating speed and IA is the abbreviation of interval analysis.
From Table 3, it is observed that the response ranges calculated by the interval method are in good agreement with those generated by the MCS, which is considered to be accurate on the premise of sufficient samples.The response bounds calculated from the interval method with different expansion orders and interpolation numbers are almost identical when two decimal digits are kept.It can be concluded that the results are accurate adequately when the expansion order is 5 and the interpolation number is 10.There is still some discrepancy between the ranges obtained from the interval method and the MCS, which is small and acceptable.It can be explained that the sampling times is finite and minor errors could be introduced by the interval arithmetic.As demonstrated in Table 4, the computing time for the MCS to obtain accurate results is much longer than the interval method which indicates the efficiency of the latter.Moreover, increasing the expansion order introduces little computation cost.The results validate the accuracy and effectiveness of the interval method.Hereafter, we used the expansion order 5 and 10 interpolation points in the interval analysis.

Non-probabilistic analysis under different uncertainties
In this part, some parameters are regarded as uncertain-but-bounded variables and the corresponding dynamic responses are analyzed.The deviation coefficient is set to 10 % and the displacement of disk 1 with respect to rotating speed is presented in the single dimensional cases.For the multi-dimensional case, the response ranges of disk 1 and disk 2 are investigated.It should be noted that any location of interest in the rotor system with different deviation coefficients can be analyzed via the solution procedure if necessary.

Uncertain support stiffness
The effects of the uncertain stiffness of support 1 are investigated in this subsection which is expressed as Eq.(17).The dynamic response bound curves under 10 % deviation coefficient are plotted in Fig. 3.As we can see from Fig. 3, the response deviates from the original deterministic one in the presence of uncertainty and the deviation is especially obvious in the resonant areas.However, the deviation is not symmetric about the deterministic curve with a larger offset in the upper bound before and after the resonant peak while near the peak point the lower bound deviates further.It can be explained that the solution algorithm searches the maximum and minimum response in the uncertain parameter interval and the deterministic resonant peak is given in the near vicinity of the critical speed while the furthest and also lowest response is set as the lower bound.The resonant peak is shifted due to inherent stiffness matrix changes which influence the critical speeds of the system directly.The response of this rotor is sensitive to the stiffness of support 1 and a resonance band appears near the second critical speed where the displacement is large.It implies that the resonant region is expanded, and the displacement may stay high in a relatively wide rotating speed range around the critical speeds.Local fluctuations of the response upper curve in the resonant region are observed which could be caused by the errors introduced in the interval algorithm and it may be influenced by the expansion orders.

Uncertain density
To account for the uncertainty introduced by material dispersion, the density of the system is considered to be uncertain and treated as interval variable as:  = , ̅ =  −   ,  −   , (18) where  and  are the nominal value and the deviation coefficient of the density, respectively.The response varying range under 10 % uncertainty of the density is plotted in Fig. 4. It can be noted that the dynamic response is also sensitive to the density and a large displacement area is observed.The resonant peak is shifted as illustrated in Fig. 4. Also, we can see that there are local fluctuations in the top band of the second critical speed which assure again the high sensitivity of the dynamic response to the uncertain parameter.

Effect of multi-source uncertainties
In this subsection, the effects of multi-dimensional uncertainties on the uncertain dynamic characteristics are analyzed.Considering multi-source of uncertainties, the stiffness of support 1 and the density are assumed to be uncertain-but-bounded parameters.The non-probabilistic response ranges are obtained using the interval analysis procedure under a typical uncertain deviation coefficient 10 %.Figs. 5 and 6 are the response bounds of disk 1 and disk 2. Multi-source uncertainties can result in unexpected high displacement vibration in a large rotating speed range.

Conclusions
This work aims at analyzing the dynamic characteristics of a double-disk rotor system with non-probabilistic uncertain parameters included.The interval method is applicable when the precise probabilistic distribution is unknown and insufficient prior information to implement stochastic methods.Two parameters in this paper are modeled as uncertain-but-bounded variables and the corresponding dynamic response ranges are investigated.The MCS is utilized to illustrate the accuracy and effectiveness of the interval method.From the results, it is clear that the interval method is efficient in solving uncertain rotor dynamics and satisfactory accuracy can be achieved via low expansion orders and a few interpolation points.Based on the simulations in the previous sections, several conclusions can be summarized as follows: 1) The uncertainties in rotor systems have significant effects on the dynamic response which evolved from deterministic value to a response range.The influence of uncertainties is notable in the resonant regions while weak in other areas.The interval method can also be used to analyze the sensitivities of the dynamic response to system parameters.
2) A high displacement band appears in several uncertain cases.The displacement is high in a rotating speed range close to the critical speeds in comparison with the original one.
3) In some cases, there are obvious shifts of the resonant peaks due to the variations of the critical speeds.4) Multiple source uncertainties can cause large offset of the dynamic response and small range uncertainties propagation has major impacts on the dynamic characteristics of the rotor systems.To understand the uncertain dynamic characteristics more thoroughly, it is advised that the uncertainties be taken into consideration in the design or maintenance of rotors.

Table 2 .Fig. 2 .
Fig. 2. The steady state response of the deterministic rotor model

Fig. 4 .
Fig. 4. Dynamic response under the density uncertainty The local fluctuations are more obvious compared with former cases.The resonant peaks of the upper bounds are shifted to the right and the lower ones to the left due to the changes in critical speeds caused by multiple uncertainties.

Table 1 .
Parameter values of the rotor system

Table 3 .
Displacement bounds of disk 1 obtained from the MCS and the interval method (μm)

Table 4 .
The average cost of computing time (s)