Accepted Manuscript Steady-state response analysis of cracked rotors with uncertain-but-bounded parameters using a polynomial surrogate method Chao Fu , Xingmin Ren , Yongfeng Yang , Kuan Lu , Weiyang Qin PII: DOI: Reference: S1007-5704(18)30265-X https://doi.org/10.1016/j.cnsns.2018.08.004 CNSNS 4618 To appear in: Communications in Nonlinear Science and Numerical Simulation Received date: Revised date: Accepted date: 14 May 2018 4 July 2018 14 August 2018 Please cite this article as: Chao Fu , Xingmin Ren , Yongfeng Yang , Kuan Lu , Weiyang Qin , Steady-state response analysis of cracked rotors with uncertain-but-bounded parameters using a polynomial surrogate method, Communications in Nonlinear Science and Numerical Simulation (2018), doi: https://doi.org/10.1016/j.cnsns.2018.08.004 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. ACCEPTED MANUSCRIPT Highlights A polynomial surrogate for the cracked rotor system with UBB parameters is constructed. The FE rotor analysis model is treated as a black box and only runs at sample points. Calculation accuracy of the surrogate is validated by comparisons with other methods. The surrogate method can deal with a large number of uncertain quantities with efficiency. AC CE PT ED M AN US CR IP T ACCEPTED MANUSCRIPT Steady-state response analysis of cracked rotors with uncertain-but-bounded parameters using a polynomial surrogate method Chao Fu, Xingmin Ren, Yongfeng Yang , Kuan Lu, Weiyang Qin Institute of Vibration Engineering, Northwestern Polytechnical University, Xi’an 710072, China Abstract CR IP T Uncertain-but-bounded (UBB) parameters are used to describe the non-probabilistic uncertainties in rotor systems. A general non-intrusive polynomial surrogate is constructed to quantify the uncertain effects of the UBB quantities on the dynamic responses. The surrogate is convenient to establish and able to deal with a large number of uncertain variables. The zeros of the Chebyshev series are used as sample points and the least square method (LSM) is employed to evaluate the coefficients. At the sample points of the polynomial surrogate, the harmonic balance method is AN US applied to obtain the sample responses (PS-HBM). During the surrogate modeling, the deterministic rotor system with a breathing crack is treated as a black box and no modifications should be made to the deterministic finite element (FE) analysis process. It needs no distribution laws and is especially helpful in small sample sized problems. Numerical simulations of the rotor system with various UBB parameters are carried out to demonstrate the effectiveness of the M surrogate. Moreover, its accuracy and efficiency are verified by comparisons with other classic sampling methods. 1. Introduction ED Keywords: rotor system; uncertain-but-bounded parameters; dynamic response; polynomial surrogate PT Uncertainty propagation has been one of the primary research interests in engineering mechanical systems [1-4]. Factors such as inherent dispersions in materials, variable working environment and fluctuating external excitations are all possible causes of parametric uncertainties [5-8]. Variations CE in physical parameters can influence the dynamic characteristics of the rotor system significantly. It is a common view that it will be beneficial for robust designing and dynamic assessment to take AC the uncertainties into consideration [9-11]. For this purpose, the challenges are finding the suitable way to model the uncertainties and developing advanced uncertainty quantification (UQ) methodologies. Generally, the uncertainties are modeled as random quantities and the UQ algorithms are often developed in stochastic ways, which use the probability distribution function (PDF) to characterize the distributions of uncertain quantities [12-16]. The most straightforward method may be the Monte Carlo Simulation and it is used as reference for other UQ algorithms in the probabilistic domain. Its major drawback is the poor convergence rate, which means the computation cost is huge. The Polynomial Chaos proposed and developed in Refs. [1, 17] are generalized into the Corresponding author. E-mail addresses: [email protected] (C. Fu), [email protected] (Y.F. Yang). ACCEPTED MANUSCRIPT uncertainty analysis of nonlinear rotor systems by Sinou et al. [18-21]. Other methods such as the perturbation skills [22], Neumann expansion [23], the stochastic response surface method [24] and the maximum entropy principle [25, 26] have been applied to UQ in rotordynamics as well. These methods all have their advantages and disadvantages. One point we should not overlook is that, to use a probabilistic method, the PDFs of uncertainties are assumed to be known beforehand or there are sufficient samples to define such distribution models. The task is not easy often due to insufficient prior data or it is too expensive to gather enough samples. In some circumstances, assumptions are made to apply the stochastic methods. However, it may be unreliable and CR IP T subjective. The interval methods or imprecise probability methods prevail in such conditions [27-31]. In interval methodologies, only the two bounds that define the extreme values of uncertainties are necessary and there is no need to clarify the evolutions in the varying ranges [32, 33]. The interval Taylor method (ITM) [7], the Chebyshev inclusion function (CIF) [34] and the interval scanning AN US method (ISM) can be used to perform interval analysis. The ITM is an intrusive and derivative-based method, which can hardly be extended to high orders because of the difficulty in obtaining the derivatives. The CIF proposed by Wu et al. [34, 35] is non-intrusive and has high efficiency. It is not suitable for high-dimensional problems, however. The ISM can be regarded as a counterpart of the Monte Carlo simulation and serving as reference roles in the non-probabilistic frame. Recently, researchers have made successful applications of the interval methods to the M uncertain rotordynamics [7, 29, 36-38]. It may be concluded that intrusive interval methods are more likely to obtain tight results than the non-intrusive ones, though they are all exposed to ED wrapping effects. However, the non-intrusive type methods can bring us much convenience in implementation and adaptation, especially when the rotor system is complicated or the solution process is complex. Non-intrusive metamodels or surrogate modeling such as the Kriging PT surrogate has attracted the attention of many researches in the past few years [39-42]. It constructs surrogate models for actual uncertain problems, which avoids direct handling of the complex CE uncertain problem. The procedure is often characterized as surrogate formulation, coefficient regression and optimal solution exploration. Only responses at some sample points are required to establish such surrogates. In some of the surrogate methods, the required number of samples may AC be high to maintain reliable and accurate in high-dimensional or high-order uncertain problems. That will restrict the application of such methods in UQ for complex rotor systems, which may be sensitive and has a large number of DOFs. In this paper, we are motivated to develop a non-intrusive polynomial-based surrogate modeling method for cracked rotor systems with UBB parameters. Instead of complicated theory or learning mechanism, a simple polynomial function is assumed to be the surrogate of the uncertain dynamic responses. The goal is to provide a direct and easy-to-follow UQ procedure for general nonlinear rotor system considering uncertainties, which are small sample sized. Another purpose of this study is to reduce the sample numbers in previous methods, which indicate the underlying ACCEPTED MANUSCRIPT computational cost. Thus, the established surrogate can deal with a large number of uncertain quantities. The reminder of this paper is as follows. The modeling of rotor-bearing-disk system with a crack based on the finite element method (FEM) is briefly described in Sect. 2. In Sect. 3, the theory and construction steps of the polynomial surrogate are presented. The numerical results and related discussions are given in Sect. 4. In Sect. 5, conclusions are summarized. 2. Modeling of rotor system based on the FEM CR IP T The FEM is widely employed in rotordynamics to establish the governing differential equations of motion. Generally, disks are considered as rigid bodies and their effects can be added by equivalent point mass and moments at proper nodes on the elastic shaft [43]. For a general rotating shaft, one can discretize it into N Timoshenko elements with N+1 nodes [44]. Four degrees of freedom [v, w, , ] in fixed coordinates, i.e. two lateral displacements and two rotation angles, AN US are present at each node with the torsional and axial vibration being ignored [40, 45]. Figure 1 gives the displacement relationships in a beam element. By assembling the elemental matrices to their respective locations, the FE equations of motion of the rotating system can be given as [43, 46] Mq(t ) (C G)q(t ) Kq(t ) F1 (t ) F2 (1) where q is the displacement vector of the system and the dot over it denotes the derivative with M respect to time. The matrices M, C, G and K are the global mass, damping, gyroscopic and stiffness matrices, respectively. The excitations on the system are expressed as the time-variant ED unbalance force F1 and the constant gravitational force F2 . Their expressions can be found in [43]. Faults are common in rotor systems, such as cubic bearing stiffness [47], misalignment [12, 48], PT rub-impact phenomenon [49, 50] and crack fault [51]. Comprehensive uncertainty analysis of rotor systems with multi faults was carried out in Ref. [48]. The crack fault has constantly attracted wide CE attentions in the field of rotordynamics [51-53]. Assuming that a crack occurs in the elastic shaft of the rotor under study and the cross section of a cracked element is shown in Fig. 2. To describe AC the size of the crack more generally, a non-dimensional crack depth is defined h R (2) where R and h are the shaft radius and the actual crack depth, respectively. For a cracked rotor system, the modeling process is ordinary except for the element bearing the crack. A typical breathing model proposed in Ref. [54] to simulate the time-varying behavior of the crack has been widely accepted and adopted by researchers [55-58]. It can be described as 1 cos( t ) g (t ) 2 where is the rotating speed of the shaft. The evolution states of the breathing crack in a full rotation of the rotor are plotted in Fig. 3, which is a cosine-type curve. The crack will be fully open (3) ACCEPTED MANUSCRIPT when the breathing function g (t ) equals to 1 while completely closed when it is 0. The neutral axis method [44] will be adopted to formulate the local stiffness matrix of the cracked element. Due to the presence of the crack, the moments of inertia of the cross-section is time-varying and then the local stiffness of the shaft is changing over time. It can be expressed as 1 k cj k j g (t )k loss k j (1 cos( t ))k loss (4) j j 2 where k j is the regular local stiffness of the beam element, k loss represents the stiffness loss of j CR IP T the element when the crack is fully open. The subscript j in Eq. (4) denotes that the crack appears at element j. Shear effect is included herein. The detailed deduction procedure will not be presented in order to be concise and the explicit expressions of them can be found in [44, 53]. Therefore, the global stiffness matrix of the rotor system considering the transverse crack will be noted by 1 K K (1 cos( t ))K c (5) 2 where K is the same as in Eq. (1), K c is the global stiffness loss matrix which is contributed AN US by k loss in the cracked element j and zeros in other elements. Replacing the original stiffness j matrix in Eq. (1), we obtain the general equations of motion of a cracked rotor-bearing-disk system Mq(t ) (C G)q(t ) Kq(t ) F1 (t ) F2 (6) The closed-form analytic solution of Eq. (6) is difficult or even impossible to get. It can be solved by the numerical integration methods or the HBM [59], which is an efficient and fast approach to M obtain the approximate analytic solutions. A lot of successful applications can be found using the HBM to predict the steady-state responses of cracked rotor systems [60-62]. ED 3. General polynomial surrogate method for uncertainty propagation analysis The modeling in the last section is performed in a deterministic sense. In this section, a general polynomial surrogate model for the dynamic response of the uncertain rotor system will be built PT taking the UBB parameters into consideration. It solves the problem in combination with the HBM (PS-HBM). The deterministic rotor system is treated as a black box and only the responses at the CE sample points, are used to construct the surrogate [42]. In view of this feature, the method has advantages of transparency, convenience in implementation and can be adapted to other uncertain engineering systems easily. Considering that the rotor system contains k uncertain parameters, the AC governing equations Eq. (6) should be rewritten as M()q(, t ) [C() G()]q(, t ) K()q(, t) F1(, t) F2 ( ) (7) where is an uncertain parameter vector, which belongs to the interval physical parameter set I . The superscript I defines an interval character. It can be further expressed as I c d x c d , (8) 2 2 xk ], x1 , , xk [ 1, 1] x [ x1 x2 in which c , d , and are the middle value, radius, lower bound and upper bound vectors ACCEPTED MANUSCRIPT of the interval vector I , respectively. The notation x denotes a k-dimensional standard interval vector. The above relationship links the actual UBB physical parameter set I to a standard mathematical variable vector x , which brings us convenience in the surrogate modeling. Based on the HBM, the uncertain steady-state dynamic response can be expressed by the finite-term Fourier series in interval form q I ( , t , x) A 0I (x) A kI (x)cos(k t ) BkI (x)sin(k t ) m (9) k 1 algebraic equations (IAEs) can be obtained LI xY I x I x where Y I is the interval Fourier coefficient vector Y I [A0I A1I B1I AmI BmI ]T CR IP T where m is the number of harmonics kept. Substituting Eq. (9) into Eq. (7), a set of interval (10) (11) The deterministic expression of the matrix and vectors in Eq. (10) can be found in Refs. [44, 53]. For the deterministic AEs, the Fourier coefficients can be determined by the Newton-Raphson AN US iteration method or direct division since Eq. (10) is linear. Considering the uncertainties, the IAEs cannot be directly solved. In the following context, we will establish the polynomial surrogate for the uncertain Fourier coefficients or the dynamic response. Firstly, we assume the uncertain response can be approximated by a continuous polynomial-form function and define a n-order k-dimensional surrogate model f x 1 i1 i2 ,ik 1 2 x x xkik , i1 , , ik 0, 1, ,n (12) M 0i1 ik n i , In Eq. (12), f x can be any of the interval Fourier coefficients appeared in Eq. (9) or all of (n k )! terms in the expression of f x . The function n !k ! has the highest order n and can be further expressed as f X T X (13) ED them in a row. There are in total N1 PT where X is the basis function vector and is the coefficient vector. They are defined as X [1 x1 xk x12 x1 x2 xk2 x1n xkn ]T CE 0, ,0 1, ,0 0, ,1 2, ,0 1,1, ,0 0, ,2 n, 0, ,n T ,0 (14) (15) The main task is to determine the coefficient vector to explicitly formulate the surrogate model. AC Moreover, the coefficients associate the surrogate function with the actual engineering problem and influence the accuracy of the former. It works similar to the response surface method [24] on the point that generating coefficients based on deterministic values calculated at some samples. In this study, we use the zeros of the Chebyshev polynomials as the candidate sample pool [3, 64]. The first-class Chebyshev polynomial can be expressed by T0 ( x ) 1, T1 ( x ) x , x [1, 1] Tn 1 ( x ) 2 xTn ( x ) Tn 1 ( x ) (16) Since the desired order of the surrogate function is n, the correspondent sample points can be chosen as the zeros of the (n+1)-order Chebyshev polynomial ACCEPTED MANUSCRIPT j , j 1, 2, , n 1 (2 j 1) j cos j , j 2( n 1) (17) In k-dimensional uncertain problem, the sample points will build a sample pool in tensor product 1 2 k (18) The number of sample points in the sample pool is N 2 (n 1)k . A sample means a call of the deterministic FE model. If all the sample points are used in the surrogate modeling, the sampling CR IP T method is similar to the CIF and works in tensor product way [3]. The underlying restriction is that the number of total samples N 2 increases exponentially in multi-dimensional interval problems, which is computationally expensive. The minimum number of points required to determine the surrogate function equals to the dimension of the unknown coefficient vector , which is N1 . Generally, N 2 N1 holds true in uncertain problems with several uncertain quantities. Then, the number of chosen samples can be varying in range [ N1 , N 2 ] . More samples can lead to higher AN US accuracy, but lower efficiency. There should be a trade-off between the accuracy and efficiency. It is suggested that when the sample number is N 3 2 N1 , the estimated results will be robust [63]. In this manner, the sampling method can be termed as the Chebyshev collocation method (CCM) [64] and the samples are randomly extracted from the sample pool. Since the sample number used is not equal to the number of the unknowns, the coefficient vector can be obtained by using the LSM. That is expressed as YN3 T M where Y Y1 Y2 (T 1 T Y (19) is the deterministic response vector obtained at the used sample set N k , which is also called collocation point set. At each collocation point, the uncertain parameters are all specified to determined values. Then, the deterministic HBM can be applied to obtain the sampled responses point by point i [ 1,i 2,i k ,i ]T (20) Yi i L i , i 1, 2, , N 3 PT ED 3 CE The matrix is of dimension N 3 N1 and is composed of the values of the basis function AC vector X at collocation points [ X1T , X T2 , , X TN 3 ] 1 1,1 1 1,2 1 1, N 3 k ,1 k ,2 1,12 1,12 1,1 2,1 1,2 2,2 k ,N 1,2N 1, N 2, N 3 3 3 3 k2,1 k2,2 1,1n n 1,2 k2, N 1,nN 3 3 kn,1 kn,2 n k , N (21) 3 Combining Eq. (12) and Eq. (19), the polynomial surrogate function is completely determined. The obtained surrogate is a simple polynomial expression and its maximum values define the desired interval Fourier coefficients ranges. For such a polynomial function in power basis, the optimization methods or the ISM can be applied to search its extreme values. Then, the interval harmonic dynamic responses of the rotor system can be expressed as ACCEPTED MANUSCRIPT Y I [Y, Y] = [min( f x max( f x] (22) As can be seen, the surrogate model function is constructed from polynomial basis and requires no complicated mathematical deductions. Its concept is simple and easy to follow. The computational procedure of the PS-HBM for dynamic response analysis is illustrated in Fig. 4. The accuracy and efficiency of the surrogate model will be tested in the numerical simulations. 4. Numerical simulations and discussions A 10-element rotor-bearing-disk system shown in Fig. 5 is studied in this paper. The crack is CR IP T located at element 3 and the system is supported at the two ends by isotropic bearings. Only the damping at the bearings is considered. The length of the rotor is 0.62 m and two rigid discs are assembled at nodes 3 and 8. A slight mass unbalance is attached to the right disc with angle / 3 . Its magnitude is given as a product of the disk mass and eccentricity. Other physical parameters are given in Table 1. In the simulations, the harmonic order in the HBM is kept as m=4 and the Ai2 Bi2 , i 1, 2, 3, 4 . Generally, AN US responses are drawn at node 2 in the form of deflection, i.e. four harmonics are sufficient for dynamics analysis of faulted rotor systems, which can be verified by Refs. [21, 48, 65]. Our extra calculations also indicate that the amplitudes of higher order harmonics are trivial compared with the first four components. To rigorously verify the harmonic solutions, a comparison with the numerical solutions obtained from the Newmark-β method is conducted. The HBM-4 solution is in good agreement with the numerical solution ( t 104 s is M required), which verifies the accuracy of the harmonic results. Suppose the rotor contains a crack =0.5 , the first backward and forward critical speeds of the system are 741.5 rad/s and 765 rad/s ED with such a crack size. In comparison, the first forward critical speed drops a little to 764.6 rad/s and the first backward critical speed changes to 737.5 rad/s when the crack stays fully open. The open crack case can be analyzed by specifying the breathing function in Eq. (3) to g (t ) 1 or PT using the procedure proposed in [44]. For reader comprehension, the simulations in the following text will be identified by case CE numbers. The variabilities of different UBB parameter sets as well as the corresponding crack depth and the surrogate order used for each case are provided in Table 2. Firstly, we consider two UBB parameters individually using the 3-order polynomial surrogate (Case 1 and 2). In single AC uncertain dimensional problems, the Chebyshev points in the candidate pool are all used to construct the polynomial surrogate. Therefore, the run times of the FE rotor model will be N3=4. Figures 6 and 7 plot the response ranges of the cracked rotor containing a crack =0.5 with 10% uncertainty in the non-dimensional crack depth and the stiffness of support #1, respectively. It can be seen that the response ranges in Fig. 6 are wider than in Fig. 7 which indicates that the uncertainty in the crack depth has greater impacts on the harmonic responses than in the stiffness of support #1. This can be viewed as sensitivity analysis of physical parameters on the response of the cracked system. In engineering context, the number of design parameters may be large and this clarification can help to prioritize those with high sensitivities. For multi-dimensional uncertain problems, we study Case 3 firstly, i.e. the uncertainties in the non-dimensional crack depth , the ACCEPTED MANUSCRIPT magnitude of the mass unbalance me and the bearing damping C. This is a three-dimensional uncertain problem with 5% variations of their nominal values and the interval parameters can be expressed as I [0.95, 1.05] c , meI [0.95, 1.05]mec , C I [0.95, 1.05]C c (23) where a superscript c denotes the nominal value. Figure 8 plots the interval harmonic responses of the rotor system using a 3-order polynomial surrogate when =0.3 . It can be observed that the responses ranges are smoothly estimated and no spurious peaks are noticed. The latter are usually observed due to overestimation or fake resonances of the IAEs. CR IP T Increase the variation abilities of these uncertain parameters to 10% and they become I [0.9, 1.1] c , meI [0.9, 1.1]mec , C I [0.9, 1.1]C c (24) The uncertain response ranges are given in Fig. 9 using the 3-order surrogate (Case 4). To provide comparison, the bounds calculated by the CIF are presented as well. The same order 3 is used in the CIF but the bounds are obtained by scanning the CIF rather than the interval arithmetic in Ref. AN US [3]. It is because that the interval arithmetic will introduce extra overestimations. For the first impression of Fig. 9, the bounds enclosed by the CIF agree well with the response ranges produced by the polynomial surrogate. As we can see, the interval response bounds are deviated significantly from the deterministic response lines with these large uncertainties. The first harmonic component is the least influenced and the other harmonic components are enveloped in wider ranges. To further illustrate the accuracy of the results in Fig. 9, the relative amplitude difference rate of (25) ED M the polynomial surrogate regarding to the ISM is presented in Fig. 10. The rate is defined as d dˆ 100% d where d̂ is the response amplitude calculated by the polynomial surrogate and the d is the response amplitude calculated by the ISM. In the ISM, 20 evenly distributed points are used for PT each uncertain dimension and the results have been verified being convergent by comparing with those of itself with higher number samples. From Fig. 10, we can find that the differences between CE the two methods are very trivial (peak values less than 1%), which certificates the excellent accuracy of the proposed surrogate method. The difference rate gets larger in the last two harmonic components and it is due to the small magnitudes of the responses where a very little error leads to AC a large rate. It should be noted that the deterministic rotor model will run N3=40 times in the proposed surrogate and N2=64 times in the CIF. In the ISM, the total sample number is 8000. From the above comparison, we can validate the good efficiency of the polynomial surrogate in solving uncertain dynamic system. Moreover, the CIF is more difficult to construct in mathematical formulation. Figure 11 shows the simulation results (Case 5) using 4-order surrogate with 1% uncertainty in the shaft density, 5% uncertainties in the non-dimensional crack depth and unbalance initial angle. We can see that the envelope pattern is different from the previous cases that resonance bands appear in the subcritical and critical speeds. In this case, the order 4 is needed. Weak unsteady ACCEPTED MANUSCRIPT fluctuations or spurious peaks will occur in the areas of resonance bands if order 3 is applied. It is because of the high sensitivity of the rotor system to the shaft density. The total sample number will be N3=70. With the same uncertainties as expressed in Eq. (23), we have calculated the interval dynamic responses, shown in Fig. 12, of the rotor system using order 5 when =0.5 (Case 6). The results verify the effectiveness of the proposed method in deeper crack problems of the rotor system with multiple UBB parameters. We can observe that the amplitudes are increased significantly at the 1/4, 1/3 and 1/2 subcritical speeds as well as at the first backward critical speed. Similarly, the bounds estimated by the CIF using the same order 5 are given for comparison. The CR IP T results generated by the two methods are almost identical. However, the efficiency is much higher in the polynomial surrogate than in the CIF, with N3=112 and N2=216. For a more comprehensive comparison, the evolution trends of the number of samples versus order in the two interval methods to solve the three-dimensional uncertain problems are illustrated in Fig. 13. We can see from Fig. 13 that the advantage of the surrogate in efficiency is dominant and will be more AN US significant with higher orders. To demonstrate the effectiveness of the proposed method in solving uncertain problems with a large number of uncertainties, we consider a case with five UBB parameters (Case 7), i.e. 5% uncertainties in the mass unbalance magnitude and the non-dimensional depth, 3% uncertainties in the stiffness of the two bearings and 2% uncertainty in Young’s modulus. The order is still 5 and the simulation results are plotted in Fig. 14. With these multiple uncertainties, the dynamic M responses of the rotor system are significantly influenced relative to the deterministic response curves. This phenomenon reflects the importance of taking uncertainties into consideration in ED designing or doing researches of rotor systems. The computation cost is reasonable since N3=504. It will be almost impossible for the ISM to perform such simulations with five UBB parameters due to the exponentially growing sampling number. The CIF is also very costly and the running PT time could be very long as N2=7776. Figure 15 gives the comparison of the computational efficiency of the proposed surrogate and the CIF in this 5-dimensional uncertain problem. It can be CE observed that the CIF is not suitable for uncertain problems with high-dimensions and high-orders, especially when the deterministic FE model is large-sized. The polynomial surrogate will demonstrate better performance in such cases. AC The above uncertain harmonic responses of the cracked rotor system can provide guidance for crack detections where non-probabilistic uncertainties are present. Compared with a linear rotor system without the crack, the typical crack signatures of the vibrations are dominant, i.e. the super harmonic response components [60, 61]. This is verified in a previous study of rotor system with a chordal crack and stochastic uncertainties [21]. Therefore, it is still valid to diagnosis crack faults in rotor systems by using the n amplitudes at subcritical speeds. On the other hand, the indicator of critical speed shifts seems unpractical as the uncertainties can cause frequency shifts as well. ACCEPTED MANUSCRIPT 5. Conclusion In this paper, the general polynomial surrogate model is formulated for the interval uncertain rotordynamics. The FE rotor model is viewed as a black box and used only at sample points, which are the zeros of the Chebyshev series. The HBM is used to obtain the deterministic solutions at the sample points. It is simple in mathematical formulation and the original solution procedure for the deterministic problem is not modified. The surrogate order required is low for the uncertain response estimations with satisfactory accuracy. Compared with the Chebyshev inclusion function CR IP T and the interval scanning method, it has higher efficiency as the number of sample points is significantly decreased. The methodology of this study will provide guidance for uncertainty analysis of small sample sized uncertain rotordynamics, especially with a large number of bounded uncertainties. The simulation results presented will be helpful for early crack detection of rotor systems. Moreover, the surrogate method developed is a generalized uncertainty propagation AN US analysis tool, which can be applied to rotor systems with other kinds of faults or multi faults. Acknowledgements We thank the anonymous reviewers sincerely for their insightful and valuable comments. This work is financially supported by the National Natural Science Foundation of China (No. 11272257) and the Fundamental Research Funds for the Central Universities (No. 3102018ZY016). M References [1] Ghanem RG, Spanos PD. Stochastic Finite Elements: A Spectral Approach. Springer; 1991. [2] Chipato E, Shaw AD, Friswell MI. Effect of gravity-induced asymmetry on the nonlinear vibration of ED an overhung rotor. Commun Nonlinear Sci Numer Simul 2018;62: 78-89. [3] Wu JL, Zhang YQ, Chen LP, Luo Z. A Chebyshev interval method for nonlinear dynamic systems under uncertainty. Appl Math Model 2013;37: 4578-91. PT [4] Elishakoff I, Ren YJ. Finite Element Methods for Structures with Large Stochastic Variations. Oxford University Press, Oxford; 2003. [5] Fu C, Ren XM, Yang YF, Qin WY. Dynamic response analysis of an overhung rotor with interval CE uncertainties. Nonlinear Dynam 2017;89: 2115-24. [6] Liao HT. Global resonance optimization analysis of nonlinear mechanical systems: Application to the uncertainty quantification problems in rotor dynamics. Commun Nonlinear Sci Numer Simul 2014; 19: AC 3323-45. [7] Ma YH, Liang ZC, Chen M, Hong J. Interval analysis of rotor dynamic response with uncertain parameters. J Sound Vib 2013;332: 3869-80. [8] Bai CQ, Zhang HY. Almost sure asymptotic stability of rotor systems subjected to stochastical axial loads. Mech Mach Theory 2012;58: 192-201. [9] Sampaio R, Soize C. On measures of nonlinearity effects for uncertain dynamical systems—Application to a vibro-impact system. J Sound Vib 2007;303: 659-74. [10] Zhang XN, Han QK, Peng ZK, Chu FL. A new nonlinear dynamic model of the rotor-bearing system considering preload and varying contact angle of the bearing. Commun Nonlinear Sci Numer Simul 2015;22: 821-41. [11] Liu BG. Eigenvalue problems of rotor system with uncertain parameters. J Mech Sci Technol 2012;26: ACCEPTED MANUSCRIPT 1-10. [12] Li ZG, Jiang J, Tian Z. Non-linear vibration of an angular-misaligned rotor system with uncertain parameters. J Vib Control 2016;22: 129-44. [13] Zhang YM, Wen BC, Liu QL. Uncertain responses of rotor-stator systems with rubbing. JSME Int J 2004;46: 150-4. [14] Dourado A, Cavalini AA, Steffen V. Uncertainty quantification techniques applied to rotating systems: A comparative study. J Vib Control 2017;107754631769855. [15] Zhou ST, Wu X, Li HG. Critical speed analysis of flexible rotor system with stochastic uncertain parameters. J Vib Eng Technol 2017;5: 319-28. Meas Diagn 2013;33: 450-5. CR IP T [16] Xiang L, Wang ZR, Tang GJ. Empirical parameter study of uncertain rotor coupling system. J Vib [17] Xiu DB, Karniadakis GE. Modeling uncertainty in flow simulations via generalized polynomial chaos. J Comput Phys 2003;187: 137-67. [18] Didier J, Sinou JJ, Faverjon B. Nonlinear vibrations of a mechanical system with non-regular nonlinearities and uncertainties. Commun Nonlinear Sci Numer Simul 2013;18: 3250-70. [19] Jacquelin E, Friswell MI, Adhikari S, Dessombz O, Sinou JJ. Polynomial chaos expansion with random AN US and fuzzy variables. Mech Syst Signal Process 2016;75: 41-56. [20] Sinou JJ, Didier J, Faverjon B. Stochastic non-linear response of a flexible rotor with local non-linearities. Int J Nonlin Mech 2015;74: 92-9. [21] Sinou JJ, Faverjon B. The vibration signature of chordal cracks in a rotor system including uncertainties. J Sound Vib 2012;331: 138-54. [22] Loève M. Probability theory. Springer-Verlag; 1978. Eng Mech 1988;114: 1335-54. M [23] Yamazaki F, Shinozuka M, Dasgupta G. Neumann expansion for stochastic finite element analysis. J [24] Isukapalli SS, Roy A, Georgopoulos PG. Stochastic response surface methods (SRSMs) for uncertainty ED propagation: Application to environmental and biological systems. Risk Anal 1998;18: 351-63. [25] Murthy R, Tomei JC, Wang XQ, Mignolet MP, El-Shafei A. Nonparametric stochastic modeling of structural uncertainty in rotordynamics: unbalance and balancing aspects. J Eng Gas Turb Power PT 2014;136: 062506. [26] Gan CB, Wang YH, Yang SX, Cao YL. Nonparametric modeling and vibration analysis of uncertain CE Jeffcott rotor with disc offset. Int J Mech Sci 2014;78: 126-34. [27] Rao SS, Berke L. Analysis of uncertain structural systems using interval analysis. AIAA J 1997;35: 727-35. AC [28] Moore RE. Methods and applications of interval analysis. SIAM Philadelphia, PA; 1979. [29] Dimarogonas AD. Interval analysis of vibrating systems. J Sound Vib 1995;183: 739-49. [30] Wang Z, Tian Q, Hu HY. Dynamics of spatial rigid–flexible multibody systems with uncertain interval parameters. Nonlinear Dynam 2016;84: 527-48. [31] Moens D, Hanss M. Non-probabilistic finite element analysis for parametric uncertainty treatment in applied mechanics: Recent advances Finite Elem Anal Des 2011;47: 4-16. [32] Khodaparast HH, Mottershead JE, Badcock KJ. Interval model updating with irreducible uncertainty using the Kriging predictor. Mech Syst Signal Process 2011;25: 1204-26. [33] Mao WG, Han X, Liu GP, Liu J. Bearing dynamic parameters identification of a flexible rotor-bearing system based on transfer matrix method. Inverse Probl Sci Eng 2016;24: 372-92. [34] Wu JL, Luo Z, Zhang YQ, Zhang N, Chen LP. Interval uncertain method for multibody mechanical systems using Chebyshev inclusion functions. Int J Numer Meth Eng 2013;95: 608-30. ACCEPTED MANUSCRIPT [35] Wu JL, Luo Z, Zhang N, Zhang YQ. A new uncertain analysis method and its application in vehicle dynamics. Mech Syst Signal Process 2015;50–51: 659-75. [36] Fu C, Ren XM, Yang YF, Xia YB, Deng WQ. An interval precise integration method for transient unbalance response analysis of rotor system with uncertainty. Mech Syst Signal Process 2018;107: 137-48. [37] Wang C, Ma YH, Zhang DY, Hong J. Interval analysis on aero-engine rotor system with misalignment, in: ASME Turbo Expo 2015: Turbine Technical Conference and Exposition, American Society of Mechanical Engineers 2015; V07AT30A002. [38] Shiau TN, Kang CH, Liu DS. Interval optimization of rotor-bearing systems with dynamic behavior CR IP T constraints using an interval genetic algorithm. Struct Multidiscip O 2008;36: 623-31. [39] Avendaño-Valencia LD, Chatzi EN, Spiridonakos MD. Surrogate modeling of nonstationary systems with uncertain properties, in: European Safety and Reliability Conference ESREL; 2015. [40] Sinou JJ, Nechak L, Besset S. Kriging metamodeling in rotordynamics: application for predicting critical speeds and vibrations of a flexible rotor. Complexity 2018: 1264619. [41] Xie LX, Liu J, Huang GL, Zhu WQ, Xia BZ. A polynomial-based method for topology optimization of phononic crystals with unknown-but-bounded parameters. Int J Numer Meth Eng 2018;114: 777-800. AN US [42] Wu JL, Luo Z, Zhang N, Zhang YQ. A new interval uncertain optimization method for structures using Chebyshev surrogate models. Comput Struct 2015;146: 185-96. [43] Friswell MI, Penny JET, Garvey SD, Lees AW. Dynamics of Rotating Machines. Cambridge University Press; 2010. [44] AL-Shudeifat MA. On the finite element modeling of the asymmetric cracked rotor. J Sound Vib 2013;332: 2795-807. M [45] Qin ZY, Han QK, Chu FL. Analytical model of bolted disk-drum joints and its application to dynamic analysis of jointed rotor. Proc Inst Mech Eng Part C J Mech Eng Sci 2014;228: 646-63. [46] Li CF, She HX, Tang QS, Wen BC. The effect of blade vibration on the nonlinear characteristics of ED rotor–bearing system supported by nonlinear suspension. Nonlinear Dynam 2017;89: 1-24. [47] Lu K. Statistical moment analysis of multi-degree of freedom dynamic system based on polynomial dimensional decomposition method. Nonlinear Dynam 2018; DOI: 10.1007/s11071-018-4303-1. PT [48] Didier J, Sinou JJ, Faverjon B. Study of the non-linear dynamic response of a rotor system with faults and uncertainties. J Sound Vib 2012;331: 671-703. CE [49] Ma H, Yu T, Han QK, Zhang YM, Wen BC, Chen XL. Time–frequency features of two types of coupled rub-impact faults in rotor systems. J Sound Vib 2009;321: 1109-28. [50] Chu FL, Lu WX. Experimental observation of nonlinear vibrations in a rub-impact rotor system. J AC Sound Vib 2005;283: 621-43. [51] Zeng J, Ma H, Zhang WS, Wen BC. Dynamic characteristic analysis of cracked cantilever beams under different crack types. Eng Fail Anal 2017;74: 80-94. [52] Al-Shudeifat MA, Butcher EA. New breathing functions for the transverse breathing crack of the cracked rotor system: approach for critical and subcritical harmonic analysis. J Sound Vib 2011;330: 526-44. [53] Lu ZY, Hou L, Chen YS, Sun CZ. Nonlinear response analysis for a dual-rotor system with a breathing transverse crack in the hollow shaft. Nonlinear Dynam 2016;83: 169-85. [54] Mayes IW, Davies WGR. Analysis of the response of a multi-rotor-bearing system containing a transverse crack in a rotor. J Vib Acoust 1984;106: 139–45. [55] Sinou JJ, Lees AW. The influence of cracks in rotating shafts. J Sound Vib 2005;285: 1015-37. [56] Friswell MI, Penny JET. Crack modeling for structural health monitoring. Struct Health Monit 2002;1: ACCEPTED MANUSCRIPT 139-48. [57] Sinou JJ, Lees AW. A non-linear study of a cracked rotor. Eur J Mech A Solids 2007;26: 152-70. [58] Gasch R. A survey of the dynamic behaviour of a simple rotating shaft with a transverse crack. J Sound Vib 1993;160: 313-32. [59] Nayfeh AH, Mook DT. Nonlinear Oscillations. John Wiley & Sons; 2008. [60] Sinou JJ. Detection of cracks in rotor based on the 2× and 3× super-harmonic frequency components and the crack–unbalance interactions. Commun Nonlinear Sci Numer Simul 2008; 13: 2024-40. [61] Han QK, Zhao JS, Lu WX, Peng ZK, Chu FL. Steady-state response of a geared rotor system with slant 1156-74. CR IP T cracked shaft and time-varying mesh stiffness. Commun Nonlinear Sci Numer Simul 2014;19: [62] Jun OS, Eun HJ, Earmme YY, Lee CW. Modelling and vibration analysis of a simple rotor with a breathing crack. J Sound Vib 1992;155: 273-90. [63] Isukapalli SS. Uncertainty Analysis of Transport-Transformation Models. The State University of New Jersey, New Brunswick; 1999. [64] Wu JL, Luo Z, Zhang N, Zhang YQ. A new sampling scheme for developing metamodels with the zeros of Chebyshev polynomials. Eng Optimiz 2015;47: 1264-88. AN US [65] Tai XY, Ma H, Liu FH, Liu Y, Wen BC. Stability and steady-state response analysis of a single rub-impact rotor system. Arch Appl Mech 2015;85: 133-48. AC CE PT ED M List of figures Fig. 1. Elemental displacements in fixed coordinates. ACCEPTED MANUSCRIPT Y Uncracked section ¦ ţ R O' h Crack section CR IP T X O AC CE PT ED M AN US Fig. 2. The cross-section of a cracked beam element. AN US CR IP T ACCEPTED MANUSCRIPT AC CE PT ED M Fig. 3. Cosine model of the breathing mechanism. ACCEPTED MANUSCRIPT Build the matrices M, C, K and F based on the FEM Determine the desired surrogate order n and the dimension of uncertainties k CR IP T Generate the sampling pool using Eqs. (17) and (18); Truncate the candidate samples to Calculate the sample responses at sample point i based on the HBM Yi i L i i=i+1 AN US Calculate the basis function vector at sample point X iT using Eq. (14) i < N3 Y N Generate the coefficient vector ( T 1 T Y CE PT ED M Formulate the surrogate f x and scan for the maximum values Ω< Ωend Y N Output interval harmonic response AC Fig. 4. The computational procedure of the surrogate method. ACCEPTED MANUSCRIPT Disk 1 Disk 2 Unbalance Crack 2 3 ... 8 9 Support 1 Support 2 CE PT ED M AN US Fig. 5. Computational model of the rotor system. AC 10 CR IP T 1 (b) (c) (d) ED M AN US (a) CR IP T ACCEPTED MANUSCRIPT AC CE PT Fig. 6. Response variations by scanning the 3-order polynomial surrogate considering 10% uncertainty in the non-dimensional crack depth: (a) 1 harmonic component, (b) 2 harmonic component, (c) 3 harmonic component, (d) 4 harmonic component. (b) (c) (d) PT ED M AN US (a) CR IP T ACCEPTED MANUSCRIPT AC CE Fig. 7. Response variations by scanning the 3-order polynomial surrogate considering 10% uncertainty in the stiffness of support #1: (red area-scanning the 3-order polynomial surrogate; black dotted linesdeterministic lines) (a) 1 harmonic component, (b) 2 harmonic component, (c) 3 harmonic component, (d) 4 harmonic component. (b) (c) (d) PT ED M AN US (a) CR IP T ACCEPTED MANUSCRIPT AC CE Fig. 8. Response variations by scanning the 3-order polynomial surrogate considering 5% uncertainties in the crack depth, the unbalance magnitude and the bearing damping: (a) 1 harmonic component, (b) 2 harmonic component, (c) 3 harmonic component, (d) 4 harmonic component. (b) (c) (d) ED M AN US (a) CR IP T ACCEPTED MANUSCRIPT AC CE PT Fig. 9. Response variations considering 10% uncertainties in the crack depth, the unbalance magnitude and the bearing damping (red area-scanning the 3-order polynomial surrogate; blue solid lines-bounds by scanning the 3-order CIF; black dotted lines- deterministic lines): (a) 1 harmonic component, (b) 2 harmonic component, (c) 3 harmonic component, (d) 4 harmonic component. (b) M AN US (a) CR IP T ACCEPTED MANUSCRIPT ED (c) (d) AC CE PT Fig. 10. Amplitude bounds difference rates between the surrogate method and the ISM considering 10% uncertainties in the crack depth, the unbalance magnitude and the bearing damping (blue solid line-lower bound; orange solid line-upper bound): (a) 1 harmonic component, (b) 2 harmonic component, (c) 3 harmonic component, (d) 4 harmonic component. (b) (c) (d) ED M AN US (a) CR IP T ACCEPTED MANUSCRIPT CE PT Fig. 11. Response variations by scanning the 4-order polynomial surrogate considering 1% uncertainty in the shaft density, 5% uncertainties in the non-dimensional crack depth and unbalance initial angle: (a) 1 harmonic component, (b) 2 harmonic component, (c) 3 harmonic component, (d) 4 AC harmonic component. (b) (c) (d) ED M AN US (a) CR IP T ACCEPTED MANUSCRIPT PT Fig. 12. Response variations considering 5% uncertainties in the crack depth, the unbalance magnitude and the bearing damping when =0.5 (red area-scanning the 5-order polynomial surrogate; blue solid AC CE lines-bounds by scanning the 5-order CIF): (a) 1 harmonic component, (b) 2 harmonic component, (c) 3 harmonic component, (d) 4 harmonic component. AN US CR IP T ACCEPTED MANUSCRIPT AC CE PT ED M Fig. 13. The evolutions of the sample numbers versus order in the surrogate method and the CIF for 3-dimensional uncertain problems. ACCEPTED MANUSCRIPT (d) ED M AN US (c) CR IP T (b) (a) AC CE PT Fig. 14. Response variations by scanning the 5-order polynomial surrogate considering 5% uncertainties in the crack depth and the unbalance magnitude, 3% uncertainties in the two bearing stiffness and 2% uncertainty in Young’s modulus: (a) 1 harmonic component, (b) 2 harmonic component, (c) 3 harmonic component, (d) 4 harmonic component. AN US CR IP T ACCEPTED MANUSCRIPT AC CE PT ED M Fig. 15. The evolutions of the sample numbers versus order in the surrogate method and the CIF for 5-dimensional uncertain problems. ACCEPTED MANUSCRIPT List of tables Table 1 Parameters of the rotor system. Description Value Description Value 0.62 m Young’s modulus 2.11011 N / m2 Radius of shaft 0.014 m Bearing stiffness 8 106 N / m Poisson’s ratio 0.3 Bearing damping Disk mass 0.6 kg Mass unbalance Radius of disk 0.078 m Shaft density me C 1 10% - - - 2 - - - 3 5% 5% 5% 4 10% 10% 10% 5 5% - - 6 5% 5% 5% 7 5% 5% - 7800 kg / m3 K1 K2 E Depth Order - - - - 0.5 3 - - 10% - - 0.5 3 - - - - - 0.3 3 - - - - - 0.3 4 1% 5% - - - 0.3 4 - - - - - 0.5 5 - - 3% 3% 2% 0.3 5 ED PT CE AC 1 105 kg m M Case 500 N s / m AN US Table 2 Variabilities of uncertain parameter sets. CR IP T Length of shaft

1/--страниц