Interval analysis of rotor dynamic response based on Chebyshev polynomials

Ynhong MA, Yongfeng WANG, Cun WANG, Jie HONG,*

a School of Energy and Power Engineering, Beihang University, Beijing 100083, China

b Collaborative Innovation Center of Advanced Aero-Engine, Beijing 100083, China

KEYWORDS Aero-engine;Chebyshev expansion;Interval analysis methods;Rotor dynamics;Uncertainty

Abstract Uncertainty is extensively involved in the rotor systems of rotating machinery, which may cause an unstable vibrational response. To take the uncertainty into consideration for the uncertain rotor-bearing system,an improved unified interval analysis method based on the Chebyshev expansion is established in this paper. Firstly, the Chebyshev Interval Method (CIM) to calculate not only the critical speeds but also the dynamic response of rotor with uncertain parameters is introduced. Then, the numerical investigation is carried out based on the developed double disk rotor model and computation procedure,and the results demonstrate the validity.But when the uncertainty is sufficiently large to influence critical speeds,the upper and lower bounds are far from the actual bounds.In order to overcome the defects,a Bound Correction Interval analysis Method (BCIM) is proposed based on the Chebyshev expansion and the modal superposition. In use of the improved method, the bounds of the interval responses, especially the upper bound,are corrected, and the comparison with other methods demonstrates that the higher accuracy and a wider application range.

1. Introduction

The analysis of rotor dynamic response is very important in the design of aero-engines. Various methods are employed to obtain the rotor response accurately; the three-dimensional finite element method is widely used currently.However,most traditional methods are only suitable for rotors with specified parameters, thereby ignoring the influence of uncertainties.The certainty hypothesis is suitable for most rotors, but the dynamic responses of rotors in aero-engines are becoming increasingly sensitive to changes in structural parameters,1for increasingly heavy loads and lightweight structures. It is necessary to take different kinds of the variability and/or uncertainty into account for dynamic analysis of aero-engine rotors.

Probabilistic approaches have been intensively studied and widely adopted to study the dynamic of uncertain structures,and different methods have been developed,such as the Monte Carlo simulation,2Stochastic finite element methods,3the perturbation method, and the stochastic collocation method.However, sufficient statistical information is usually need to be identified beforehand to get the parameter probability distributions.4Meanwhile,the computational costs of probabilistic approaches like the Monte Carlo simulation are still too high for the design of aero-engine rotor systems. In this case,non-probabilistic methods can be alternatively employed.The interval analysis method5is one of the typical nonprobabilistic methods, which only requires bounds on the parameters,and it seems to be much more effective to analyze the dynamic characteristic of uncertain structures by the lower and upper bounds. Thus, the interval analysis method has drawn more and more attention in recent years.

The interval concept and arithmetic first appeared in the 1960s, but it is 1990s when Dimarogonas6and Rao7et al.made early effort to consider the effects of uncertain parameters use of the interval calculus in the mechanical vibrating system. Then, Moens8,9and Rao10et al. applied the interval method to frequency response envelope analysis of uncertain structures. Commonly, the interval arithmetic method will experience overestimation in the computation because of its intrinsic wrapping effect.Massa et al.11has proposed an extension of Moens’works and develops a general uncertainty propagation method for interval or fuzzy eigenvalue and superposition analysis12to avoid classical problems with interval arithmetic.Many other schemes are required to be incorporated into the interval analysis methods to control the overestimation, such as the Interval Taylor series method5,13and Taylor model method.14Qiu et al.15,16conducted a number of studies in structural analysis with uncertainty. The parameter perturbation method17and Taylor’s expansion method18were developed in their studies in the case of the uncertainty being slight, while second-order Taylor’s expansion19and an updated second-order Taylor series expansion20are also used to achieve higher accuracy.Hu et al.21introduced a new Taylor expansion-based Interval Pattern Analysis(TIPA) to analyze the pattern tolerance of the linear array antenna with excitation amplitudes uncertainty, which is simpler and can produce more accurate bounds of the antenna pattern and the pattern features.

Under the weak uncertainty condition,the interval analysis approach based on Taylor expansion or perturbation method is sufficiently accurate. However, the range of the dynamic response of the structure may be overestimated when the uncertainty increases, or a large number of uncertain parameters are combined. Meanwhile, the Taylor expansion method must be applied on the partial differential form of the original equation,which makes it almost impossible to be incorporated into commercial finite element software.

To overcome the overestimation problem in Taylor-seriesbased algorithms,Qiu and Qi22further proposed a collocation interval analysis method that is based on the Chebyshev polynomials, which was developed by Wei et al.,23,24and used for the torsional vibrations of a wind-turbine geared transmission system.Wu et al.25also proposed a Chebyshev inclusion function method for multibody systems with interval uncertainties to control overestimations much better. They introduced the Chebyshev inclusion function into nonlinear systems with large uncertainty to obtain the dynamic response ranges,26,27which can offer a higher numerical accuracy in the approximation of solutions. A mixed probability and interval reliability analysis method was proposed by Gao et al.28to calculate the mean and variance of random interval variables for engineering problems, and sparse regression was also used to improve the efficiency and accuracy.29However, the control of the overestimation remains unsatisfactory, especially on the lower bound of the response for the system with multiple uncertain parameters.

The interval analysis methods are widely used to solve the uncertainty problem in practical engineering. In structural analysis, Qiu et al.30,31have used interval set models in the study of the static response and eigenvalue problems of structures with bounded uncertain parameters. Chen and Yang32have presented one interval finite element method based on interval mathematics to solve uncertain problems of the beam structures. The methods were also used to study the response statistics of engineering structures with interval parameters under stationary random excitation.22Moreover, the interval analysis technique can be proposed for structural damage identification33or structural optimization34as well.Dimarogonas6introduced interval modal analysis,and interval solution of the eigenvalue problem for interval rotor-bearing systems and solutions were developed using interval calculus. Ma et al.35-37extended the interval analysis method to rotor dynamics, and they studied the modal characteristics and steady and transient responses via a perturbation method, the interval modal superposition method and Taylor expansion, respectively.While, Han et al.38,39used the Chebyshev interval method in the analysis of a parametrically excited system, as well as in the study of the start-up transient process of an overhung rotor.40Most of analyses based on Chebyshev polynomials are applied on a simple example to get one specific dynamic characteristic,the unified interval analysis method for studying both the critical speeds and response of rotors still need to be developed.

In this paper, a Bound correction Chebyshev Interval Method (BCIM) based on the interval analysis of eigenvalues and modal superposition is proposed,which not only increases the accuracy but also unifies the critical speed analysis and the dynamic response analysis of the rotor.The nonintrusive Chebyshev algorithm is presented, which can be used on complicated models of engineering structures. The improved method is applied to a finite-element model of a double disk rotor system and the critical speeds and dynamic response under parameter uncertainty are analyzed, which are compared with the results of classical methods. Based on the numerical results,the validity of the improved interval analysis of the rotor dynamic response based on Chebyshev polynomials is demonstrated, and it can be used in the design of aeroengine rotors with multiple large uncertainties.

2. The Chebyshev interval method (CIM) for dynamic analysis of uncertain rotors

2.1. Uncertain parameters and the interval equation of rotorbearing systems

Variability and/or uncertainty exist inevitably in aero-engine rotor-bearing systems, due to manufacturing and assembly errors,material dispersion and changes in working conditions,lubricant properties, and other factors. The varieties of joints in rotor and support structures, as shown in Fig.1(a), are one of the main source of uncertainty. The dynamic properties of joints are functions of the geometry, material, surface condition,load,etc.,41and the sparse experimental data on the joint geometries that have been studied show huge variability.These data were collected from a few specimens; testing on larger data sets has demonstrated even higher variability.Meanwhile,the wear in joint interfaces and rubbing between the rotor and sealing coatings cause the imbalance to be an uncertain parameter.

In addition, the stiffness and damping of the supports typically vary in different operating conditions, especially when a squeeze film damper is used.42The support stiffness is a significant parameters affecting the dynamic characteristics of rotors, which generally consists of three parts38: the support structure stiffness(Kstructure),the stiffness of squeeze film damper (KSFD)and the stiffness of bearing (Kbrg). Considering the different assembly state and loads,as well as the wear of interfaces,the Kstructureand Kbrgwill change.The lubricant properties such as density and viscosity vary with the oil temperature and loads, which cause the changes of KSFD. The parameters aforementioned are typical ‘‘uncertain but bounded” ones,the interval method is suitable and efficient for the dynamic analysis of the rotor-bearing systems.

The typical high-pressure rotor-bearing systems of aeroengine could be simplified as a rotor system with two offset disks and supported at both ends,as shown in Fig.1(b),where E and ρ are the elastic modulus and density of the material,the m1,2and g1,2are the mass and gyroscopic of different Disk,k1,2and c1,2are the stiffness and damping of different supporting structures. The differential equation of the free vibration of a rotor-bearing system with n degrees of freedom can be represented as

where M=(mij) is the mass matrix; C=(cij) is the damping matrix; G=(gij) is the gyroscopic matrix, which is a function of the rotation speed Ω;K=(kij)is the stiffness matrix; and z is the displacement vector. The dynamic analysis of the rotor system with uncertain parameters can be described by Eq. (1) with constraint:

where, mij, cij, gij, and kijare mass, damping, gyroscopic and stiffness for different parts, while - and _ are the upper-value and lower-value of the parameters.

2.2. Chebyshev method for obtaining the interval natural frequencies

Introduce a 2n-dimensional vector y:

Fig. 1 Schematic diagram of rotor-bearing system and simplified dynamic model.

2.3. Algorithm for interval analysis of the rotor dynamic response

The differential equation that describes the forced vibration of the rotor can be represented as

The construction of Chebyshev inclusion function is similar to that described in Section 2.2. The interval solutions are obtained by combining deterministic solutions Z at interpolation points. The definitions of the interval parameters and interpretation points are same as in Eqs. (14)-(16).

3. Numerical analysis of the rotor dynamic response based on CIM

As shown in Fig.2,the rotorg system is modeled by beam element, mass element and spring element with interval parameters. Detailed deterministic geometric and physical parameters of the rotor studied are given by Table 1.

3.1. Interval analysis of rotor critical speeds

The rotor critical speeds that were calculated via the Chebyshev method will be presented and compared to those that were obtained via combination methods and the perturbation(adopted in the Ref. [38]), to demonstrate the advantages and disadvantages of the method.The support stiffness k1,the density ρ, and the elastic modulus E of rotor are set as interval variables, while the other parameters are deterministic. The interval support stiffness k1=[k1c-βk1c, k1c+βk1c], the interval density ρ=[ρc-βρc, ρc+βρc], and the interval elastic modulus E=[Ec-βEc, Ec+βEc], where the nominal values are k1c=1.0×107N.m-1, ρc=7.8×103kg.m-3and Ec=200 GPa. The uncertain coefficient β is set as 0.05 and 0.15.

The first two mode shapes of the rotor system are shown in Fig.3,while the Campbell diagram(as shown in Fig.4)shows the variation of the eigenfrequencies that correspond to the mode shapes with respect to the rotation speed,which are calculated via the Chebyshev method with the elastic modulus E as an interval variable. The critical speeds correspond to the intersection points between frequency curves and the line f=ω/60,where ω is the circular frequency.The critical speeds are no longer fixed values due to the variation of elastic modulus E, while the 2nd-order critical speed with the bending vibration of rotor are more sensitive to the change of E.

The critical speeds calculated with different uncertain parameters(E,ρ,k1)and different uncertain coefficient values(β=0.05,0.15)are listed in Table 2,and compared in Figs.5 and 6. As is known, the combination method could get the accurate interval results in the linear rotor system by combining the lower or upper value of different parameters, in case the influence of parameters is monotonous,7which is chosen as the benchmark.

Fig.5 shows the interval of critical speeds affected by different uncertain parameters. The results show that, the elastic modulus E has a substantial influence on the second critical speed for this rotor, while the support stiffness k1almost has no influence on it,and the variation of the density ρ has significant influence on both critical speeds. As the degree of uncertainties increases, the errors among the results of the methods increase,while the Chebyshev method obtains higher accuracy in almost all situations compared to perturbation method, as shown in Fig. 6, when ρ is an uncertain parameter. The order and the integration points of the Chebyshev method have almost no effect on the results, in the interval eigenfrequency analysis,lower-order and fewer integration points yield results that satisfy the accuracy requirements.

3.2. The dynamic response of the uncertain rotor

In addition to the analysis of critical speeds,the corresponding interval analysis method for the response with imbalance excitation was also established, and the computation procedure was developed, as shown in Fig. 7, based on the algorithm in Section 2.3. In this part, direct interval analysis on the model that is illustrated in Fig. 2 will be performed with a single uncertain parameter. The constant structural damping coefficient γ for the system is 0.01; therefore, the damping matrix C in Eq. (24) is written as C=2γK/Ω.

The value of uncertain coefficient β is set as 0.05.According to the interval results for the critical speeds in Table 2,the harmonic response is calculated between 50 Hz and 350 Hz, and the step is set as 0.2 Hz. The step has a substantial influence on the accuracy of the results because it determines whether the results that are calculated at the integration points(the system has various parameter values) yield the real response amplitudes or not.

Shown in Fig. 8 are the interval responses at the Disk 2,when k1is an uncertain parameter with uncertain coefficient β=0.05, and the parameters m in Eq. (15) and n in Eq. (20)of the Chebyshev method are set as 3 and 2, respectively.The Monte Carlo method is employed to evaluate the effectiveness of the Chebyshev method.As we can see,the response of rotor is not fixed but an interval value due to the variation of parameters,especially at the frequency closed to the critical speeds when the response is huge.For the conciseness,only the dynamic response around the critical speeds are given out and analyzed below.

Fig. 2 Finite-element model of double disk rotor system.

Table 1 Values of physical parameters of double disk rotor.

By the Chebyshev Interval Method (CIM), the uncertainty characteristics due to the uncertain coefficient can be acquired,and the intervals of the critical speeds are similar with the results get by Monte Carlo. But, the accuracy of the dynamic response still needs to be improved, especially for the lower bounds. As we can see, the dynamic responses are negative due to the error of the algorithm, which won’t appear to happen in practice.

The dynamic response calculated with different uncertain parameters (k1, E, ρ) are gave out in Figs. 8-10. The influence of the uncertainty varies with different parameters,and it is different at different critical speed of the rotor as well. Similar to the critical speed, the support stiffness k1,as well as E and ρ, has a substantial influence on the response around the second critical speed for the rotor. Meanwhile,the accurate of results by CIM related to the sensitivity of critical speed to Parameters. If the parameter has a small influence on the critical speeds, the response estimate is accurate, but if the critical speeds have a wide range, the upper and lower bounds are far from the actual bounds due to the error of the algorithm.

Take the influence of uncertain support stiffness k1as an example, to illustrate the defects of the algorithm. Fig. 11(a)shows the response of the system when the support stiffness k1is evaluated at three integration points, while Fig. 11(b)shows the bounds, the constant term in Eq. (29) and the response when the elastic modulus is set to the nominal value,which corresponds to the result when j1=2 in Fig. 11(a).

According to Eq.(16),when m=3,the support stiffness k1is set as 0.95×107, 1.0×107and 1.05×107N.m-1.

As shown in Eq.(28),the interval result at various rotation speeds is the linear combination of the results with various values of k1. The result at 272 Hz in Fig. 11(b) can be used to illustrate the defects of the algorithm: The responses at 272 Hz that are calculated with various values of k1are all small(Fig.11(a)).Therefore,their linear combination is small.Increasing the number of integration points may improve the accuracy. The parameters m in Eq. (15) and n in Eq. (20) of the Chebyshev method are set as 10 and 8, respectively. The results are recalculated and shown in Fig. 12.

The order of the Chebyshev method and the number of integration points have small influences on the interval response near the 1st critical speed, while the accuracy of the upper bound near the 2nd critical speed is improved.The computational cost increases dramatically, especially when there are multiple uncertain parameters, while the accuracy remains unsatisfactory. In the next part, a correction method for the upper bound will be proposed.

Fig. 3 Mode shapes of rotor.

Fig. 4 Critical speeds with interval elastic modulus (β=0.15).

Table 2 Critical speeds of rotor with different uncertain parameters calculated via combination,the perturbation,and the Chebyshev method.

Fig. 5 Critical speeds of rotor with different uncertain parameters.

4. The bound correction method and application in multiple uncertainty systems

4.1. The bound correction Chebyshev interval method (BCIM)

Substantial attention should be paid to the bound of the response, especially the upper bound, because the vibration amplitude may affect the safe operation of rotating machinery. Therefore, a method for estimating the upper bound based on the modal superposition method will be proposed to overcome the disadvantages of the classical Chebyshev algorithm.

In Appendix B, the modal superposition method for calculating the response of system is presented. The dynamic response of the damped rotor system is expressed as follows with the first k modes:

Fig. 6 Relative error of Critical speeds calculated via different methods when ρ is an uncertain parameter.

Fig. 7 Computation procedure of CIM for rotor dynamic response analysis.

The notations are defined in the Appendix. In Eq. (31),ai/bi=λi,aiis purely imaginary number and biis real number.When λi=iΩ, the system is in the resonance of the ith mode.If the frequencies of other modes are far from the frequency of the ith mode, the response when λi=iΩ can be simplified as follows:

Consider the uncertainty of the parameters and assume that the interval parameters have little influence on the modes of the system. Hence, the eigenvectors φiare the same as those when the parameters are set to the nominal values.

If M,K,and G are interval matrices,ai,biand λiwill be interval numbers(Eq.(B.3)).The upper and lower bounds of eigenfrequency λiwere obtained in Section 3. When iΩ ∈λiI, the upper bound of the response can be solved via Eq.(32),where riis the modal damping.If the modal damping of the ith mode or the complex load is uncertain,the bound of the response can be obtained easily due to the monotony of Eq.(32).

Meanwhile, the lower bound is corrected by the interval critical speeds. When the uncertain parameters cause the critical speeds to reach the bounds, the response that is calculated with these parameter values can form the lower bound. For a single uncertain parameter that has a monotonous influence on the critical speeds, combining minimum value of the response with the bounds of the uncertain parameter will yield an accurate lower bound.If there are multiple uncertain parameters in the system, the parameters that yield the most accurate approximation to the interval critical bounds among all integration points will be used to calculate the lower bound.

The bound correction method is applied to the results in Figs.13-15 when iΩ ∈λiI,while the interval response far away from the critical speed interval is solved via the Chebyshev method with m=3 and n=2. The Monte Carlo results may differ from those in Section 3.2 due to the randomness.

Figs. 13-15 show the bounds after correction, and the results shows that the upper and lower bounds are more accurate compared to the results in Figs.8-10 via classical CIM.It seems that the interval upper bound exceeds the real upper bound (Fig. 13(a)) in some frequencies. This is caused by the frequency step, which is not small enough to capture the real amplitude, and in Figs. 14 and 15, we calculated the results with more steps. Meantime, the interval bound is exceeded by the Monte Carlo results at various speeds, which is caused by the assumption that the mode of the system is unchanged and the use of only the dominant mode to calculate response.Therefore,Eq.(32)can only be used when there is a wide margin between the critical speeds.

4.2. Application in a rotor system with multiple uncertain parameters

4.2.1. The critical speeds of rotor

The critical speeds of the rotor with multiple interval parameters are listed in Table 3, with parameters m in Eq. (15) and n in Eq.(20)of the Chebyshev method are set as 3 and 2,respectively. As the number of uncertain parameters increases, the difference in efficiency increases. The results show that, the Chebyshev method estimates the upper bound with higher accuracy in all situations, and is better than the perturbation method for the lower bound in most cases.

Fig. 8 Amplitudes of disk at various rotation speeds (k1 with β=0.05).

Fig. 9 Amplitudes of disk at various rotation speeds (E with β=0.05).

Fig. 10 Amplitudes of disk at various rotation speeds (ρ with β=0.05).

Fig. 11 Amplitude of disk (k1 with β=0.05).

Fig. 12 Amplitudes of disk at various rotation speeds (E with β=0.05).

Fig. 13 Bounds of amplitude after correction (k1 with β=0.05).

Fig. 14 Bounds of amplitude after correction (E with β=0.05).

Fig. 15 Bounds of amplitude after correction (ρ with β=0.05).

It also should be noted that the eigenvectors are assumed to be constant in perturbation method,which means it could not be applied to the system when the uncertain parameter has a great influence on vibration mode, but the Chebyshev method is adoptable in this situation.Meanwhile,the intrusive calculation of perturbation method is based on MATLAB or other computational languages.Once the amount of degrees of freedoms is huge,the efficiency will decrease sharply.However,the Chebyshev method is based on the results which could be acquired from commercial finite element software directly, it is more appropriate for the complicated model.

4.2.2. Dynamic responses of rotor with multiple uncertainties

The proposed method can also be applied to analyze the dynamic response of a rotor system with multiple uncertain parameters, while the interval critical speeds that are calculated in Table 3 can be used to correct the bound. Figs. 16 and 17 show the interval response when the elastic modulus E and density ρ are interval parameters with β=0.05. The parameters m in Eq. (15) and n in Eq. (20) of the Chebyshev method are again set as 3 and 2.

Table 3 Critical speeds of rotor with multiple interval parameters

Fig. 16 Bounds of amplitude near the 1st-order critical speed.

Fig. 17 Bounds of amplitude near the 2nd-order critical speed.

The results demonstrate that the accuracy of the Chebyshev method decreases as the number of uncertain parameters increases,while the proposed correction method(BCIM)could provide more precise bounds and be suitable for scenarios with multiple uncertain parameters. Combined with the classical correction method (CIM), a low-order Chebyshev method with few integration points could be used between critical speeds to improve the efficiency of the interval response analysis.

5. Conclusions and discussion

This paper established an interval analysis method for the rotor dynamics based on the Chebyshev expansion, including the interval analysis of critical speeds and an interval algorithm for determining the dynamic response. Different from previously proposed analysis methods, the critical speeds and the response are unified into a single framework for analysis,and the bounds of the response were corrected. The following conclusions could be drawn:

(1) The Chebyshev method for uncertain rotor-bearing systems obtains higher accuracy on critical speeds, especially on the upper bound, compared to the perturbation method. Meanwhile, the proposed method is much more appropriate as structure models are complicated, because it is based on results that can be acquired directly from commercial finite-element software.

(2) The double disk rotor model and computation procedure based on Chebyshev polynomials were developed.The classical Chebyshev method(CIM)is suitable under the weak uncertainty condition, but when the uncertainty is strong and the parameters have a great influence on the dynamics characteristics (such as the critical speeds), the method failed to obtain a satisfactory result for the dynamic response, and the upper and lower bounds are far from the actual bounds, even with a high expansion order and additional integration points.

(3) A corresponding bound correction method (BCIM)based on modal superposition is proposed in this paper,which could provide more precise bounds.The presented algorithm and the Monte Carlo approach are compared and the numerical results are in good agreement, which demonstrates the validity of the correction method.

The proposed interval analysis method, which is based on the Chebyshev expansion, is a valid and practical method for the dynamic study of rotors with uncertain parameters. And the proposed bound correction method is a positive development in rotor dynamics and can be used in scenarios with multiple uncertain parameters and substantial uncertainty.However,several aspects still needed to be improved in future,especially the applicability when the margin between critical speeds is small. Meanwhile, the application in complex rotors with large DOFs, as well as the influences of the local uncertain parameters, local stiffness and/or damping of the bolted joint for example, are still need deeply investigated.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This study was supported by the National Natural Science Foundation of China (Nos. 51575022, 11772022 and 11672017).

Appendix A

A.1. Chebyshev polynomial expansion theory