An interval precise integration method for transient unbalance response analysis of rotor system with uncertainty

Mechanical Systems and Signal Processing 107 (2018) 137–148

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

An interval precise integration method for transient unbalance response analysis of rotor system with uncertainty Chao Fu a, Xingmin Ren a, Yongfeng Yang a,⇑, Yebao Xia a, Wangqun Deng b a b

Institute of Vibration Engineering, Northwestern Polytechnical University, Xi’an 710072, China AECC Hunan Aviation Powerplant Research Institute, Zhuzhou 412002, China

a r t i c l e

i n f o

Article history: Received 3 May 2017 Received in revised form 15 January 2018 Accepted 18 January 2018

Keywords: Transient response Rotors Precise integration Uncertainty Interval analysis

a b s t r a c t A non-intrusive interval precise integration method (IPIM) is proposed in this paper to analyze the transient unbalance response of uncertain rotor systems. The transfer matrix method (TMM) is used to derive the deterministic equations of motion of a hollow-shaft overhung rotor. The uncertain transient dynamic problem is solved by combing the Chebyshev approximation theory with the modified precise integration method (PIM). Transient response bounds are calculated by interval arithmetic of the expansion coefficients. Theoretical error analysis of the proposed method is provided briefly, and its accuracy is further validated by comparing with the scanning method in simulations. Numerical results show that the IPIM can keep good accuracy in vibration prediction of the start-up transient process. Furthermore, the proposed method can also provide theoretical guidance to other transient dynamic mechanical systems with uncertainties. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction High speed and heavy power are the development directions of modern rotating machineries. Transient unbalance response analysis is a crucial and indispensable part in rotor dynamic researches because rich information is contained in the start-up or other speed varying process. Many works can be found on the evaluations of rotor transient dynamic behaviors during the past few decades. Bouaziz et al. [1] studied the transient behaviors of an active magnetic bearings-rotor system with angular misalignment using the Newmark numerical method. To reduce the transient imbalance vibrations in rotating machines, Kim and Na [2] designed a new ball-type automatic balancer to improve the performance of traditional ones. The nonlinear response of rub effects in rotor system was investigated by lots of researchers [3–6]. Zapomeˇl et al. [7] analyzed the transient behaviors of a flexibly supported rigid rotor considering unbalance and damping of short magnetorheological squeeze film dampers. Crack and misalignment effects in rotor dynamic characteristics are frequently investigated [8,9]. Other aspects including clearance [10], base motions [11], damping ratio identification [12] and viscoelastic properties [13] were also studied. The common feature of the above works is that only deterministic models were considered. However, due to manufacturing and assembly errors, material dispersion and changes in working conditions, uncertainty exists inevitably in rotor systems. In order to obtain reasonable evaluations of transient characteristics, it is highly recommended that uncertainty should be taken into consideration. Presently, the probabilistic methods have been intensively studied and widely employed in

⇑ Corresponding author. E-mail address: [email protected] (Y. Yang). https://doi.org/10.1016/j.ymssp.2018.01.031 0888-3270/Ó 2018 Elsevier Ltd. All rights reserved.

138

C. Fu et al. / Mechanical Systems and Signal Processing 107 (2018) 137–148

uncertain rotordynamics. Didier et al. have done many profound researches [14–16] in which rotor dynamic response is approximated by Fourier series and the Polynomial Chaos Expansion (PCE). A comparative study was carried out by Dourado et al. [17] for uncertainty quantification in rotating systems. Focused on the local nonlinearity, Sinou et al. [18] studied the stochastic nonlinear response of an uncertain flexible rotor system. The vibration signatures in diagnosis of chordal crack were verified to be still dominant in the presence of uncertainties [19]. Furthermore, an insightful study of the influence of the expansion order for the polynomial chaos on an asymmetric rotor was conducted [20]. Based on experimental data, Faverjon et al. [21] used the PCE for model validation of a stochastic damped structure with the measurement uncertainties. Jacquelin et al. [22] extended the PCE for representing the response of dynamical systems with mixed random and fuzzy variables. Murthy et al. [23] studied the influence of nonparametric stochastic uncertainties on dynamic characteristics at low rotating speeds. Gan et al. [24] used Monte Carlo method to analyze the sensitivity of the first order critical speed to uncertain parameters. It should be pointed out that probabilistic methods are applicable only when the precise parameter probability distributions are given. However, as frequently encountered, the complete statistic information of uncertain parameters is often not available. In practical, on the other hand, the varying ranges of parameters are easier to define which can be modeled by unknown-but-bounded interval variables. Qiu et al. [25,26] have carried out plenty of studies on interval uncertain structural dynamics using Taylor expansion and perturbation techniques. Qi and Qiu [27] further proposed a collocation interval analysis method based on the Chebyshev polynomials to cope with the overestimation problem in Taylor series based algorithms. Wu et al. [28,29] introduced the Chebyshev inclusion function into multibody mechanical problems and high accuracy and efficiency were achieved. Moreover, it is not restricted to small uncertainties compared with Taylor series and has excellent controlling in the wrapping effect. It has been used to solve uncertain dynamic problems such as gear systems [30]. Up until now, few researches apply interval methods to rotor transient dynamic analysis. Ma et al. [31,32] studied the dynamic response of rotors using perturbation method, Taylor expansion and the interval modal superposition method. Another dilemma encountered in transient dynamic research is the accuracy and stability problems of numerical integration methods such as Newmark-b and Runge-Kutta methods. Very small integration step is required to obtain reasonable results occasionally and these methods are not capable of dealing with singular or ill conditioned matrices. Zhong [33,34] proposed the PIM for dynamic response calculations, which can get exact solution at each time step. It has been proved that it is unconditionally stable and its precision is independent of step length [35]. Very wide applications of PIM can be found for problems governed by constant coefficient differential equations [36,37] and it has been successfully applied to rotordynamics [38]. Further, some applications in time-varying problems emerged based on modified PIMs [39,40]. Yue [41] extended PIM to time varying structures with rheonomic transition matrices and implemented Magnus expansion [42] to compute the exponential matrix. It is shown that it possesses high accuracy and efficiency. Yet, it has not been applied to uncertain rotordynamics. In this work, we are devoted to introducing interval analysis into the uncertain transient problem of an overhung rotor system and the IPIM method is proposed by coupling with the modified precise integration method. The hollow shaft rotor system and the derivation of the deterministic equations of motion are illustrated in Section 2. The proposed IPIM method is explained in detail in Section 3. Based on the procedure above, numerical results are given for the rotor system and some discussions are made in Section 4. At last, conclusions are summarized in Section 5.

Case

Shaft 2

Shaft 1

L

Bearing 1

I

Disk

Bearing 2

II

III,IV

Fig. 1. Schematic diagram of the casing rotor. (I: left end of the rotor, II: left cross-section of the disc, III: right cross-section of the disc, IV: right end of the rotor.)

C. Fu et al. / Mechanical Systems and Signal Processing 107 (2018) 137–148

139

2. Rotor model and deterministic equations of motion The rotor system under study is comprised of a hollow shaft supported by two elastic bearings and an overhung horizontal disc, as shown in Fig. 1. The lumped-mass disc is located at the right end of the shaft and the shaft is of step diameters. The TMM [43–45] is vastly utilized in rotating machines for critical speed calculation and building the general equations of motion due to its simplicity in structure and convenience to programming. Different sub-structures are correspondent to certain pattern of transfer matrices. General equations of motion and displacements of each node can be derived by considering boundary conditions. For the rotor shown in Fig. 1, a small length of shaft (DL
L) is added at the ends of the rotor to use the free end boundary conditions. Assume the shaft is isotropic and there is no torsion in cross-section plane, the transfer mode for shaft elements in the X direction is illustrated in Fig. 2. It is similar in the Y direction. In Fig. 2, notations x; b; M; Q are the displacement, slope angle, bending moment and shear force of the respective crosssection within the XOZ plane, respectively. The hollow property of the shaft is reflected in the polar moment of inertia Is ¼ pðD4  d Þ=64 where D and d are the outer diameter and the inner diameter, respectively. The disc element is considered as equivalent external excitations and the increment of the cross-section displacement vector can be expressed as 4

DZ ¼ ½0; 0; DMx ; DQ x ; 0; 0; DM y ; DQ y T where

ð1Þ

8 € þ Ip xa_ > DM x ¼ I d b > > > > > DQ x ¼ m€x  cx_ þ F x > > > < DM ¼ I a _ y d €  I p xb €  cy_ þ F y > DQ y ¼ my > > > > 2 > _ € sinð/ðtÞt _ > F ¼ me½ x cosð/ðtÞt þ /Þ þ / þ /Þ x > > : 2 _ € _ F y ¼ me½x sinð/ðtÞt þ /Þ  / cosð/ðtÞt þ /Þ þ mg

where b; a are the angular displacement in plane XOZ and YOZ, x is the rotating speed, m is the disc mass, Ip is the polar mass _ moment of inertia, Id is the transverse mass moments of inertia, c is the damping coefficient, e is the eccentricity, /, /ðtÞ and _ € /ðtÞ are the initial unbalance angle, angular velocity and acceleration, respectively. Generally, we assume x ¼ /ðtÞ . A short shaft of length 0.001 m is added at both ends of the rotor. The transfer matrix from the left end of the rotor to the left cross-section of the disc is expressed by T1 and T2 represents the transfer matrix for the right cross-section of the disc to the right end of the rotor, let T ¼ T2 T1 . When the physical parameters are given, these transfer matrices are determined. According to the transfer flow, we have

Zr ¼ TZl þ T2 DZ

ð2Þ

Zd ¼ T1 Zl þ DZ

ð3Þ

where subscripts r, l and d denote the cross-sections of the right, left end of the rotor and the disc, respectively. For clarity in the forth-coming description, we define matrices D, S and P which stand for the transfer matrices T, T1 and T2 , respectively. Divide them into blocks of 2  2 matrices as

2

D11 6 . 6 D ffi T ¼ 4 .. D41

3 2    D14 S11 6 .. .. 7 7; S ffi T1 ¼ 6 .. . . 5 4 .    D44 S41

3 2    S14 P11 6 .. .. 7 7; P ffi T2 ¼ 6 .. . . 5 4 .    S44 P41

Fig. 2. Sketch map of the shaft element.

3    P14 .. .. 7 7 . . 5    P44

ð4Þ

140

C. Fu et al. / Mechanical Systems and Signal Processing 107 (2018) 137–148

Considering the free boundary conditions at the two ends (M ¼ 0; Q ¼ 0), the deflection and angular displacements of the disc can be expressed as

½x b y aT ¼ Te ½DMx DQ x DM y DQ y T where

ð5Þ

8   S11 S13 > 1 > T ¼ T T T ; T ¼ > e a c a b < S31 S33     > D P22 P24 D 21 23 > > ; Tc ¼ : Tb ¼ D41 D43 P42 P44

At this stage, substituting the transfer matrices into Eq. (5) and combining with Eq. (1), the final differential equations of motion of the disc is as follows

€ þ CðtÞUðtÞ _ MUðtÞ þ KðtÞUðtÞ ¼ FðtÞ

ð6Þ

where M, CðtÞ and KðtÞ represent the mass, damping and stiffness matrices of the rotor system, respectively. FðtÞ denotes the _ € excitation forces. UðtÞ, UðtÞ, UðtÞ are the displacement, velocity and acceleration vectors, respectively. They can be expressed as

2

c 60 6 M ¼ diagðm; Id ; m; Id Þ; CðtÞ ¼ 6 40

0 0 0

3 0 0 0 Ip xðtÞ 7 7 7 c 0 5

0 Ip xðtÞ 0 2

0

1

6 1 0 6 KðtÞ ¼ K0 T1 e ; K0 ¼ 6 4 0 0 0

0 0 0

0

0

3

07 7 T T 7; UðtÞ ¼ ½x; b; y; a ; FðtÞ ¼ ½F x ; 0; F y ; 0 15

0 1 0

where x; b; y; a are the deflection and angular displacements of the disc geometric center in x and y direction, respectively. F x and F y represent the external excitations and are detailed in Eq. (1). Taking the mass property of the shafts into consideration, Eq. (6) is time varying when the rotor system is accelerating up as the generalized damping and stiffness matrices are changing with time. The gyroscopic effect is integrated into the damping matrix. 3. Interval precise integration method for transient analysis of uncertain rotor system In this section, the IPIM based on the extension of precise integration and Chebyshev orthogonal polynomial approximation is explained for transient unbalance response bounds determination of rotor systems subject to interval uncertainties. The major difference in the precise integration compared with previous PIMs is that the coefficients are time variant [39,41]. By Chebyshev basis representation, the non-probabilistic uncertainties will be modeled without probabilistic distributions. Moreover, it is not constrained to small uncertainties compared with perturbation technique or Taylor expansion and the overestimation caused by the wrapping effect is well controlled [25,27]. The proposed mythology will preserve these advantages and solve the uncertain transient dynamic problem effectively. 3.1. General procedure of the IPIM Suppose the parameters of the rotor system such as support stiffness, unbalance eccentricity and damping are uncertain and there is no sufficient statistical data to define their specific probabilistic distributions, thus the stochastic methods in probabilistic frame will be inapplicable. In this manner, interval uncertainties are advocated. We do not consider the uncertain initial conditions for the rotor system, which is solvable with the proposed method but has little practical application background. By collecting all the possible uncertain parameters, an uncertain parameter vector can be defined

q ¼ ½q1 ; q2 ; . . . ; qr T

ð7Þ

where r denotes the number of uncertain parameters. The ith member of q can be expressed in interval form as

i  ¼ ½qci  bi qci ; qci þ bi qci ; i ¼ 1; 2; . . . ; r qIi ¼ ½qi ; q

ð8Þ

i are the lower and upper bounds, respectively. qci is the mean value of the parameter and bi represents the where qi and q uncertain degree. The general equations of motion then can be rewritten as

C. Fu et al. / Mechanical Systems and Signal Processing 107 (2018) 137–148

€ _ MðqÞUðq; tÞ þ Cðq; tÞUðq; tÞ þ Kðq; tÞUðq; tÞ ¼ Fðq; tÞ

141

ð9Þ

The transient response of the system can be taken as the solution of Eq. (9) subject to additional constraints

n o € þ CU _ þ KU ¼ F; Uðt 0 Þ ¼ U0 ; q 2 qI UI ðq; tÞ ¼ U : MU

ð10Þ

in which U0 is the initial state of the transient response vector, qI represents the uncertain parameter vector set expressed by Eq. (8). The goal here is to search the bounds to determine the response range. Chebyshev series is used to deal with the interval uncertainties. Original mechanical input-output system is approximated as a surrogate with a few interpolation points. Firstly, we present the basics of Chebyshev polynomial approximation. The first class Chebyshev polynomial basis is defined by

C n ðxÞ ¼ cosðn arccos xÞÞ; x 2 ½1; 1

ð11Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi where n denotes the nonnegative integer. fC n ðxÞg are orthogonal on interval [1,1] with weight function qðxÞ ¼ 1  x2 . In order to use the bounded characteristics of trigonometric functions, the following replacement is made

h ¼ arccos x 2 ½0; p

ð12Þ

Eqs. (11) and (12) suggest that the Chebyshev polynomials are all bounded in the interval [1,1] within their definition range in the light of j cosðnhÞj 6 1, which will be used in the interval arithmetic in the last calculation step. The one-dimensional Chebyshev polynomial approximation can be expressed as

Pn ðhÞ ¼

n a0 X ai C i ðnhÞ; h 2 ½0; p þ 2 i¼1

ð13Þ

where ai represents the coefficient of Chebyshev polynomial expansion which can be calculated by Mehler integration. They can be expressed as

ai ¼

2

p

Z p

ðUðt; cos hÞ cos ihÞdh; i ¼ 1; 2; . . . ; n

ð14Þ

0

in which Uðt; cos hÞ is the deterministic value at interpolation points. The step will associate the approximation function with the deterministic dynamic model. The integration can be further calculated by Gauss-Chebyshev integration formulas

Z p

ðUðt; cos hÞ cos ihÞdh ¼

0

k X p k¼1

k

Uðt; cos hi ÞC i ðcos hi Þ; i ¼ 1; 2; . . . ; n

ð15Þ

with k being the number of interpolation points. We adopt the zeros of the Chebyshev polynomials as the interpolation points

hj ¼

2j  1 p; j ¼ 1; 2; . . . ; k 2k

ð16Þ

where h is the interpolation point, k represents the number of the interpolation points used. It is required that k P n þ 1 in order to be accurate for the numerical integration. The transient response range can be calculated by searching the bounds of the simple approximation expression. In multi-dimensional form, the approximation of Eq. (10) can be expressed as

UI ðqI ; tÞ ¼

n X



k1 ¼0

n X 1 I p Ck1 ;;kr C k1 ;;kr ðq Þ 2 kr ¼0

ð17Þ

where p denotes the number of zeros appeared in the subscripts, n is the expansion order. C k1 ;...;kr ðqI Þ represents the multidimensional Chebyshev polynomials which can be defined as

C k1 ;...;kr ðqI Þ ¼ cosðk1 h1 Þ cosðk2 h2 Þ    cosðkr hr Þ

ð18Þ

Ck1 ;...;kr is the expansion coefficient which can be calculated by Gauss-Chebyshev quadrature. The following expression can be obtained

Ck1 ;;kr ¼

k X j1 ¼0



k  r X 2 ~ Uðt; cos hj1 ; . . . ; cos hjr Þ cosðk1 hj1 Þ    cosðkr hjr Þ k j ¼0

ð19Þ

r

~ cos hj ; . . . ; cos hj Þ is the deterministic transient response when the uncertain parameters take the value as the where Uðt; 1 r interpolation points. For real cases, the uncertain parameter within arbitrary interval can be transferred to standard interval via

hI ¼ arccos

2xI  ðx þ xÞ I ; h 2 ½0; p; xI ¼ ½x; x x  x

ð20Þ

142

C. Fu et al. / Mechanical Systems and Signal Processing 107 (2018) 137–148

where superscript I defines an interval character, x and  x are the lower and upper bounds of xI . In this way, the convex set of the independent interval variables in multi-dimensional problem are transformed from super rectangular to super cube model which has the same number of interpolation points on each side. Now, we have transformed the uncertain transient problem into a set of deterministic motion equations at the interpolation points and the structure of them remain unchanged compared with Eq. (6). In order to obtain accurate response results and improve the solution stability, a modified PIM [41] can be used to deal with the time-varying coefficients for the deterministic differential equations at the specified interpolation points. In state space, the deterministic ODEs at interpolation points is expressed as

m_ ¼ AðtÞm þ fðtÞ in which



AðtÞ ¼

ð21Þ

0

I

M1 KðtÞ M1 CðtÞ

  ; fðtÞ ¼

   U ; m ¼ M1 FðtÞ U_ 0

where I is unit matrix. AðtÞ is changing with time which is different from normal cases. For two consecutive time steps, the general recursive solution of Eq. (21) is given as [35]

mkþ1 ¼ gkþ1 mk þ

Z

tþdt t

gkþ1 fðsÞds

ð22Þ

where gðtÞ being the time-varying state transition matrix, dt is the integration time step, mkþ1 and mk are the transient response at times t kþ1 and tk , respectively. A matrix exponential term is used to solve Eq. (22) numerically, which can be expressed as

gk ¼ expðdðtÞÞ

ð23Þ

In order to avoid the non-homogeneous terms, the increment form is introduced

3 U AðtÞ fðtÞ 6 7 ; m~ ¼ 4 U_ 5 HðtÞ ¼ 0 0 1 

2



ð24Þ

Then, Eq. (21) is combined with Eq. (24) as

m~_ ¼ HðtÞm~

ð25Þ

To cope with the time-varying features, the fourth order Magnus expansion [42] is introduced

8 pffiffiffi > H1 ¼ Hðt þ ð1=2  3=6ÞdtÞ > > pffiffiffi > > > > < H2 ¼ Hðt þ ð1=2 þ 3=6ÞdtÞ pffiffiffi ~dðtÞ ¼ dt=2ðH1 þ H2 Þ þ 3ðdtÞ2 =12½H2 ; H1  > > > > þðdtÞ3 =80½H2  H1 ; ½H2 ; H1  > > > : gk ¼ expð~dðtÞÞ

ð26Þ

in which operation ½A; B ¼ AB  BA: Thus, the time-varying transition matrix is represented by the function of the two specific times. The most important part of the PIM is the calculation of the exponential term and it is the key to achieve high accuracy. In actual calculation, the time step is discretized into even intervals. The 2N algorithm is used to compute it precisely, which applies the addition theorem 1=l expð~dðtÞÞ ¼ ½expð~dðtÞlÞ

ð27Þ

where l ¼ 1=2N , N is the cycle parameter which is suggested to be 20 typically [39]. In this way, the calculation error can be neglected reasonably. To improve the stability and avoid unnecessary matrix inversion, the state transfer matrix is expressed in block matrices during the cycle computing

g~ kþ1 ¼ expð~dðt þ dtÞÞ ¼ I þ Ta ; Ta ¼



TaA

TaF fðtÞ

0

0



ð28Þ

Then Ta is expanded into Taylor series up to fourth order

(

2

3

4

TaA ¼ AðtÞl þ ðAðtÞlÞ =2! þ ðAðtÞlÞ =3! þ ðAðtÞlÞ =4! 2

TaF ¼ Il þ AðtÞl =2! þ AðtÞ2 L3 =3! þ AðtÞ3 L4 =4! In this pattern, we perform the iteration for N times. The iteration formula is as follows

ð29Þ

143

C. Fu et al. / Mechanical Systems and Signal Processing 107 (2018) 137–148



TaFðiþ1Þ ¼ 2TaFðiÞ þ TaAðiÞ  TaFðiÞ TaAðiþ1Þ ¼ 2TaAðiÞ þ TaAðiÞ  TaAðiÞ

; i ¼ 1; 2; . . . ; N

ð30Þ

After the cycle calculation completed, left out the last line, the final response vector can be given as

mkþ1 ¼ ðI þ TaA Þmk þ TaF

ð31Þ

By transferring back in the expression of Eq. (24), the exact transient response values are obtained for deterministic interpolation points

~ cos hj ; . . . ; cos hj Þ ¼ m ð1 : 4Þj Uðt; cos hj 1 r

1

;...;cos hjr

ð32Þ

Combining Eq. (32) and Eq. (17), we can obtain the interval solution at every time step based on the Chebyshev inclusion function [29]

UI ¼

X 1 1 jCk ;...;k j r C0;...;0 þ ½1; 1 2 2p 1 r 06k1 ;...;kr 6r

ð33Þ

k1 þ...þkr P1

where p is the number of zeros in the subscripts. In this paragraph, we take a general view of the whole analysis procedure. The uncertain governing equations of motion are transformed into a set of deterministic equations based on Chebyshev series with original response values at interpolation points. At each interpolation points, the modified PIM is utilized to solve the time-varying system numerically and then construct the Chebyshev coefficients. From here, we can see that the size and structure of all the differential equations solved in the analysis procedure are constant and it is convenient for programming calculation. Another advantage is that it can achieve higher accuracy via only adding the number of interpolation points and higher expansion orders without reconstructing the differential equations. It can deal with large uncertainties without probabilistic distributions. Moreover, it is of general interest to other time-varying uncertain problems. The algorithm of the IPIM for uncertain transient rotor systems subject to uncertainties is demonstrated in Fig. 3.

3.2. Error analysis In this subsection, we give a theoretical analysis of the calculation error of the IPIM. Its accuracy verification in numerical simulations will be demonstrated in the next section. The multi-dimensional problems can be regarded as the tensor form of one-dimensional cases. Therefore, single dimensional problem will be considered here without loss of generality. Three sources of errors may occur besides the error caused by the computer data storage, which can be expressed as

n ¼ n1 þ n2 þ n3

ð34Þ

where n1 , n2 and n3 are the inherent errors of the Gauss-Chebyshev truncation, the PIM and the interval arithmetic denoted as Eq. (33). The Gauss-Chebyshev quadrature truncation error n1 can be given by [27]

n1 ¼

2p 22k ð2kÞ!

where f ðxÞ and f f

ðnþ1Þ

f

ð2kÞ

ð2kÞ

ðxÞ; x 2 ½1; 1

ð35Þ

ðxÞ are the derivable distribution function and the 2k-order derivative function. It can be found that

ðxÞ will be zero because f ðxÞ is of order n. Then n1 will be zero as long as k P ðn þ 1Þ=2 for an odd n or k P ðn þ 2Þ=2

for an even n [27]. Thus, accuracy improves as the order n increases. In the PIM, the accuracy is high due to the 2N algorithm 4

and the possible error n2 comes from the truncation of the Taylor expansion. It is minor and can be estimated as ðTa lÞ =5! based on the next term of the expansion [35]. The maximum error will be introduced by n3 in the interval arithmetic at the last calculation step. It is difficult to describe from the mathematical aspect and the truncation errors are submerged in n3 .

4. Numerical results In this section, the transient response based on the deterministic parameters will be given firstly. Then several typical cases will be investigated where different parameters are regarded as interval uncertain variables. In the first uncertain case, the accuracy and efficiency of the proposed method are validated by comparative study with the scanning method, which scatters equal distanced sampling points in the uncertain interval. Some discussions are made according to the results.

144

C. Fu et al. / Mechanical Systems and Signal Processing 107 (2018) 137–148

...

Fig. 3. Computation procedure of the IPIM for uncertain transient analysis.

4.1. Deterministic transient response of the rotor Detailed deterministic geometric and physical parameters of the overhung rotor under study are given by Table 1. Herein, we are interested in the start-up transient process of the rotor with constant angular acceleration in an unbalanced state. In € the start-up process, the angular acceleration /ðtÞ is 55 rad/s2 and the initial rotating speed x0 ¼ 0. Set integration time step dt ¼ 1  103 s, the deterministic transient unbalance response of the disc is shown in Fig. 4. The transient vibration peak is 1:03  103 m at 3.82 s. 4.2. Transient response variability in presence of uncertainties In this part, the transient response ranges will be presented when different parameters are assumed to be uncertain with a typical uncertain degree 10%. To illustrate the variability caused by uncertainties, the deterministic response curve will be plotted in all cases for direct reference. The expansion order and the interpolation number are n ¼ 3 and k ¼ 4.

C. Fu et al. / Mechanical Systems and Signal Processing 107 (2018) 137–148

145

Table 1 Values of physical parameters of the rotor. Symbol

Description

Value

E

Young’s modulus of elasticity

L1 L2 D1 D2 D0 K1

Length of shaft 1 Length of shaft 2 Outer diameter of shaft 1 Outer diameter of shaft 2 Inner diameter of shafts Stiffness of support 1

2:1  1011 N=m2 0.3 m 0:12m 0:02m 0:03m 0.01 m

K2

Stiffness of support 2

C m e

Viscous damping Disk mass Eccentricity

q

Density Polar mass moment of inertia Transverse mass moment of inertia

Ip Id

1  108 N/m 1  106 N/m 120 N s=m 8.4 kg 8  105 m 7800 kg=m3 0:0695 kg m2 0:0357 kg m2

Fig. 4. Deterministic accelerating up transient unbalance response of the rotor.

Fig. 5. Transient response range considering uncertain eccentricity and viscous damping.

4.2.1. Uncertainties in the disc mass imbalance and viscous damping In this subsection, the eccentricity and the viscous damping at the disc are taken as uncertain parameters to account for uncertainties in external excitation and damping effects. After calculation, the uncertain transient response range is shown in Fig. 5. To validate the proposed method, the transient vibration bounds are compared with those calculated by the scanning method, which are assumed to be of satisfactory accuracy. Detailed data of the lower bounds (LB) and upper bounds (UB) at several time steps is given in Table 2. The scanning point number for each uncertain parameters is 10 and 20.

146

C. Fu et al. / Mechanical Systems and Signal Processing 107 (2018) 137–148

Table 2 Deflection bounds of the disc obtained from the scanning method and the IPIM (lm). Time (s)

2 3 3.5 3.8 4 5

Scanning (10 points)

IPIM (n ¼ 3; k ¼ 4)

Scanning (20 points)

LB

UB

LB

UB

LB

UB

34.05 160.08 501.36 860.60 523.37 162.13

41.67 197.82 653.87 1225.63 722.53 198.67

34.05 160.08 501.35 860.60 523.36 162.12

41.67 197.83 653.87 1225.63 722.54 198.67

34.04 159.88 497.28 839.96 510.92 162.08

41.67 197.85 653.87 1225.64 722.56 198.68

It can be observed from Fig. 5 that when the rotor is subject to uncertainty, the transient vibration response is significantly affected which is especially evident near the resonant area. The deterministic response curve is evolved to a response range defined by the upper and lower bounds due to the presence of uncertainties. The vibration peak deviates with a large offset. Other than the main resonant area, the attenuation region is also influenced obviously. The results from the scanning method with 20 points are reference for the scanning method itself with 10 scanning points and the IPIM. When 10 scanning points are used for each uncertain parameter, the rotor analysis model should be solved 102 ¼ 100 times in this two-dimensional problem. The computing time elapsed was 15.22 min. While in the IPIM, the iterative times will be 42 ¼ 16 and 2.78 min was cost for computing. Thus, we can see the advantage of the IPIM. For higher uncertain dimensions and complex mechanical systems engaging large analysis models, it will be more evident as the computational efforts grow geometrically in the scanning method. As the response bounds given in Table 2 illustrated, the proposed method can achieve good accuracy with low expansion order and a few interpolation points. The bounds estimated by the IPIM agree well with those of the scanning method in non-resonant areas while near the resonant peak there is overestimation. This is caused by the interval arithmetic in the last step of the IPIM in response bounds determination. However, as indicated in Table 2, it is less than 2.4% which is acceptable given the efficiency. 4.2.2. Uncertainties in the bearing stiffness and Young’s modulus The support stiffness in rotor system is always difficult to be precisely defined and the material property is also in dispersion. In that case, we will analyze the effects of the uncertainties in the stiffness of support 1 and Young’s modulus. The variability of the transient response is demonstrated in Fig. 6. From Fig. 6, the same influence mechanism is observed that the uncertainties have greater effects on the resonant region and the attenuation areas than other parts. A phenomenon, which is different form the case in Section 4.2.3, is that there are peak shifts in the resonant region with the upper resonant peak shifted to the right and the lower one to the left. This indicates the variation of the system critical speed due to the change of sensitive parameters. Specifically, there are two resonant peaks in the upper response curve. Similar effects were found in the vibration response dynamic system subject to interval uncertainties [29,30]. It is mostly because of the high sensitivity of the natural frequencies to uncertain parameters considered. 4.2.3. Uncertainties in the inner shaft diameter and shaft length Affected by errors in the manufacture process, some geometric parameters in the rotor system may subject to uncertainty. Thus, in this subsection, the inner diameter of the elastic shaft and the length of shaft 1 are modeled as unknown but bounded uncertain parameters.

Fig. 6. Transient response range considering uncertain bearing stiffness and Young’s modulus.

C. Fu et al. / Mechanical Systems and Signal Processing 107 (2018) 137–148

147

Fig. 7. Transient response range considering uncertain inner shaft diameter and shaft length.

Fig. 7 shows the uncertain transient response range of the overhung rotor when the inner shaft diameter and the length of shaft 1 are uncertain. It is indicated that the deviation of the response curve is weak compared to former cases, which suggests less impacts of the chosen uncertain parameters on the transient start-up vibration. However, there are still slight shifts of the resonant peak, which indicates variation of the natural frequency. 5. Conclusions In this paper, the non-intrusive IPIM for transient vibration prediction of rotor system in the start-up process is proposed by coupling the interval analysis method with the modified precise integration method. Uncertainties in this study are treated as interval variables, which overcomes the constraints of probabilistic methods when precise probabilistic distributions are unavailable. The uncertain time-varying equations are transformed into a set of deterministic ODEs with the same size of the original ones using Gauss-Chebyshev quadrature points. The modified PIM is used to obtain the deterministic transient response values. The accuracy and efficiency are validated by comparing with the scanning method. It can obtain satisfactory results with low expansion order. High uncertain dimension analysis can be performed conveniently and with high efficiency at the same time. Sampling methods like the scanning method will be awkward in high dimensional problems because of large iterations of the deterministic analysis model, which requires a great deal of computational efforts. Some cases have been studied where different parameters are considered uncertain. Some parameters have significant effects on the uncertain response ranges. The proposed method can also be used to analyze the sensitivity of the system response to parameters. As can be observed from the results, the IPIM can also lead to some overestimations. Future work may concentrate on developing improved algorithms, which can obtain more tight wrapping bounds without losing the efficiency. Acknowledgements We thank Prof. Dayi Zhang from Beihang University and Dr. Kuan Lu from Northwestern Polytechnical University for their valuable suggestions during the preparation of this paper. This work acknowledges the financial support from National Natural Science Foundation of China (No. 11272257), Aeronautical Science Foundation of China (No. 2013ZB08001) and the Fundamental Research Funds for the Central Universities (No. 3102016ZY016). References [1] S. Bouaziz, N.B. Messaoud, J.Y. Choley, M. Maatar, M. Haddar, Transient response of a rotor-AMBs system connected by a flexible mechanical coupling, Mechatronics 23 (2013) 573–580. [2] T. Kim, S. Na, New automatic ball balancer design to reduce transient-response in rotor system, Mech. Syst. Signal Process. 37 (2013) 265–275. [3] F. Chu, W. Lu, Experimental observation of nonlinear vibrations in a rub-impact rotor system, J. Sound Vib. 283 (2005) 621–643. [4] S. Roques, M. Legrand, P. Cartraud, C. Stoisser, C. Pierre, Modeling of a rotor speed transient response with radial rubbing, J. Sound Vib. 329 (2010) 527– 546. [5] H. Ma, C. Shi, Q. Han, B. Wen, Fixed-point rubbing fault characteristic analysis of a rotor system based on contact theory, Mech. Syst. Signal Process. 38 (2013) 137–153. [6] C. Wang, D. Zhang, Y. Ma, Z. Liang, J. Hong, Theoretical and experimental investigation on the sudden unbalance and rub-impact in rotor system caused by blade off, Mech. Syst. Signal Process. 76 (2016) 111–135. [7] J. Zapomeˇl, P. Ferfecki, P. Forte, A computational investigation of the transient response of an unbalanced rigid rotor flexibly supported and damped by short magnetorheological squeeze film dampers, Smart Mater. Struct. 21 (2012) 105011. [8] Y. Yang, H. Chen, T. Jiang, Nonlinear response prediction of cracked rotor based on EMD, J. Franklin I (352) (2015) 3378–3393. [9] P. Pennacchi, A. Vania, S. Chatterton, Nonlinear effects caused by coupling misalignment in rotors equipped with journal bearings, Mech. Syst. Signal Process. 30 (2012) 306–322. [10] H. Li, G. Meng, Z. Meng, B. Wen, Effects of boundary conditions on a self-excited vibration system with clearance, Int. J. Nonlin. Sci. Num. 8 (2007) 571– 580.

148

C. Fu et al. / Mechanical Systems and Signal Processing 107 (2018) 137–148

[11] Q. Han, F. Chu, Dynamic behaviors of a geared rotor system under time-periodic base angular motions, Mech. Mach. Theory. 78 (2014) 1–14. [12] W. Wang, Q. Li, J. Gao, J. Yao, P. Allaie, An identification method for damping ratio in rotor systems, Mech. Syst. Signal Process. 68–69 (2016) 536–554. [13] M.I. Friswell, J.K. Dutt, S. Adhikari, A.W. Lees, Time domain analysis of a viscoelastic rotor using internal variable models, Int. J. Mech. Sci. 52 (2010) 1319–1324. [14] J. Didier, J.-J. Sinou, B. Faverjon, Nonlinear vibrations of a mechanical system with non-regular nonlinearities and uncertainties, Commun. Nonlinear Sci. Numer. Simul. 18 (2013) 3250–3270. [15] J. Didier, J.-J. Sinou, B. Faverjon, Multi-dimensional harmonic balance with uncertainties applied to rotor dynamics, J. Vib. Acoust. 134 (2012) 061003. [16] J. Didier, J.-J. Sinou, B. Faverjon, Study of the non-linear dynamic response of a rotor system with faults and uncertainties, J. Sound Vib. 331 (2012) 671– 703. [17] A. Dourado, A. Cavalini Jr, V. Steffen Jr, Uncertainty quantification techniques applied to rotating systems: a comparative study, J. Vib. Control. 1077546317698556 (2017). [18] J.-J. Sinou, J. Didier, B. Faverjon, Stochastic non-linear response of a flexible rotor with local non-linearities, Int. J. Nonlin. Mech. 74 (2015) 92–99. [19] J.-J. Sinou, B. Faverjon, The vibration signature of chordal cracks in a rotor system including uncertainties, J. Sound Vib. 331 (2012) 138–154. [20] J.-J. Sinou, E. Jacquelin, Influence of Polynomial Chaos expansion order on an uncertain asymmetric rotor system response, Mech. Syst. Signal Process. 50–51 (2015) 718–731. [21] B. Faverjon, P. Ladevèze, F. Louf, An updating method for structural dynamics models with uncertainties, Shock Vib. 15 (2013) 245–256. [22] E. Jacquelin, M.I. Friswell, S. Adhikari, O. Dessombz, J.-J. Sinou, Polynomial chaos expansion with random and fuzzy variables, Mech. Syst. Signal Process. 75 (2016) 41–56. [23] R. Murthy, J.C. Tomei, X.Q. Wang, M.P. Mignolet, A. El-Shafei, Nonparametric stochastic modeling of structural uncertainty in rotordynamics: unbalance and balancing aspects, J. Eng. Gas Turb. Power. 136 (2014) 062506. [24] C. Gan, Y. Wang, S. Yang, Y. Cao, Nonparametric modeling and vibration analysis of uncertain Jeffcott rotor with disc offset, Int. J. Mech. Sci. 78 (2014) 126–134. [25] Z. Qiu, L. Ma, X. Wang, Non-probabilistic interval analysis method for dynamic response analysis of nonlinear systems with uncertainty, J. Sound Vib. 319 (2009) 531–540. [26] Z. Qiu, X. Wang, Parameter perturbation method for dynamic responses of structures with uncertain-but-bounded parameters based on interval analysis, Int. J. Solids Struct. 42 (2005) 4958–4970. [27] W. Qi, Z. Qiu, A collocation interval analysis method for interval structural parameters and stochastic excitation, Sci. China Phys. Mech. 55 (2012) 66– 77. [28] J. Wu, Z. Luo, Y. Zhang, N. Zhang, L. Chen, Interval uncertain method for multibody mechanical systems using Chebyshev inclusion functions, Int. J. Numer. Methods Eng. 95 (2013) 608–630. [29] J. Wu, Y. Zhang, L. Chen, Z. Luo, A Chebyshev interval method for nonlinear dynamic systems under uncertainty, Appl. Math. Model. 37 (2013) 4578– 4591. [30] S. Wei, J. Zhao, Q. Han, F. Chu, Dynamic response analysis on torsional vibrations of wind turbine geared transmission system with uncertainty, Renew. Energy 78 (2015) 60–67. [31] Y. Ma, Z. Liang, M. Chen, J. Hong, Interval analysis of rotor dynamic response with uncertain parameters, J. Sound Vib. 332 (2013) 3869–3880. [32] C. Wang, Y. Ma, D. Zhang, J. Hong, Interval analysis on aero-engine rotor system with misalignment, ASME Turbo Expo 2015: Turbine Technical Conference and Exposition, Montréal, Canada, V07AT30A002. [33] W. Zhong, J. Zhu, X. Zhong, On a new time integration method for solving time dependent partial differential equations, Comput. Method Appl. M. 130 (1996) 163–178. [34] W. Zhong, J. Lin, Q. Gao, The precise computation for wave propagation in stratified materials, Int. J. Numer. Methods Eng. 60 (2004) 11–25. [35] W. Zhong, On precise time-integration method for structural dynamics, J. Dalian Univ. Technol. 34 (1994) 131–136 (in Chinese). [36] C.C. Caprani, A modal precise integration method for the calculation of footbridge vibration response, Comput. Struct. 128 (2013) 116–127. [37] J. Zhang, Q. Gao, S. Tan, W. Zhong, A precise integration method for solving coupled vehicle–track dynamics with nonlinear wheel–rail contact, J. Sound Vib. 331 (2012) 4763–4773. [38] W. Li, H. Lv, C. Pei, S. Xia, Precise numerical integration of nonlinear rotor system with multi-degree of freedom, J. Vib. Eng. 4 (2004) 55–60 (in Chinese). [39] R.E. Bartels, A time integration algorithm based on the state transition matrix for structures with time varying and nonlinear properties, Comput. Struct. 81 (2003) 349–357. [40] R. Zhao, K. Yu, An efficient transient analysis method for linear time-varying structures based on multi-level substructuring method, Comput. Struct. 146 (2015) 76–90. [41] C. Yue, X. Ren, Y. Yang, W. Deng, A modified precise integration method based on Magnus expansion for transient response analysis of time varying dynamical structure, Chaos Solitons Fract. 89 (2016) 40–46. [42] S. Blanes, F. Casas, J.A. Oteo, J. Ros, The Magnus expansion and some of its applications, Phys. Rep. 470 (2009) 151–238. [43] S.T. Choi, S.Y. Mau, Dynamic analysis of geared rotor-bearing systems by the transfer matrix method, J. Mech. Des. 123 (2001) 562–568. [44] H. Cao, L. Niu, S. Xi, X. Chen, Mechanical model development of rolling bearing-rotor systems: a review, Mech. Syst. Signal Process. 102 (2018) 37–58. [45] G. Genta, Dynamics of Rotating Systems, Springer, New York, 2005.