International Journal of Rotating Machinery
Volume 2014 (2014), Article ID 126843, 8 pages
http://dx.doi.org/10.1155/2014/126843

Research Article

Comparison of Common Methods in Dynamic Response Predictions of Rotor Systems with Malfunctions

Hongliang Yao, Qian Zhao, Qi Xu, and Bangchun Wen

School of Mechanical Engineering and Automation, Northeastern University, Shenyang 110819, China

Received 14 June 2014; Accepted 17 November 2014; Published 9 December 2014

Academic Editor: Tariq Iqbal

Copyright © 2014 Hongliang Yao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The efficiency and accuracy of common time and frequency domain methods that are used to simulate the response of a rotor system with malfunctions are compared and analyzed. The Newmark method and the incremental harmonic balance method are selected as typical representatives of time and frequency domain methods, respectively. To improve the simulation efficiency, the fixed interface component mode synthesis approach is combined with the Newmark method and the receptance approach is combined with the incremental harmonic balance method. Numerical simulations are performed for rotor systems with single and double frequency excitations. The inherent characteristic that determines the efficiency of the two methods is analyzed. The results of the analysis indicated that frequency domain methods are suitable single and double frequency excitation rotor systems, whereas time domain methods are more suitable for multifrequency excitation rotor systems.

1. Introduction

Numerical simulation plays an important role in the study of rotor systems with malfunctions. With complex structures of the modern rotor systems (e.g., multidisc or parallel shafts) and the inherent strong nonlinear characteristics of malfunctions [1] (e.g., the self-excitation property of oil whirl and oil whip, the parametric vibration characteristics of the crack rotor, and the nonlinear stiffness of rotor to stator rub), few analytical solutions work well in analyzing these strongly nonlinear systems with large degree of freedom (DOF). Therefore, numerical simulation may be the only way to predict the response of these systems accurately.

The numerical simulation methods commonly used in rotor dynamics can be divided into two types: time domain methods and frequency domain methods. The time domain methods mainly include the Runge-Kutta method [2], Newmark- method [3], and the shooting method [4]. The frequency domain methods mainly include the harmonic balance (HB) method [5], describing function (DF) method [6], incremental harmonic balance (IHB) method [7], and so forth.

When analyzing large-scale rotor systems with complex structures, both the time domain and the frequency domain methods will encounter the challenge of solving a large system of equations with numerous interdependent variables. Therefore, a dimension reduction approach must be adopted. There are two kinds of dimension reduction approaches that are widely used: the dynamical substructure approach and the receptance-based approach. The former is often combined with time domain methods and the latter is often combined with frequency domain methods. The fixed interface component mode synthesis approach [8] (CMS) is the most widely used dynamic substructure approach, in which the structure is divided into the linear part and the nonlinear part. Modal truncation is applied to the linear portion and only the lower-order modes are retained, which reduces the dimension of the former structure. The frequency domain methods often achieve dimension reduction using linear receptance data [9, 10], in which the vibrations of the linear parts are substituted by those of the nonlinear parts. The dimension of the former system then reduces to equal the number of nonlinear DOFs.

Both the types of methods have advantages and disadvantages. Frequency domain methods are highly efficient because they can skip the transient response and obtain the steady state response directly. However, these methods can only seek periodic and quasiperiodic responses. They also generally need prior information on the behavior of the system, such as which harmonic terms in the responses are dominant [11]. Time integration methods are capable of seeking periodic, aperiodic, and quasiperiodic solutions as well as transient responses. However, they have low efficiency when solving for steady state responses since the transient responses must be solved first.

Modern rotor systems have various forms and excitation conditions. Therefore, selecting an appropriate numerical method with good precision and high efficiency is very important when analyzing a rotor system with fault. Few studies have compared the efficiency and accuracy of the time and frequency domain methods or discussed the application scope of these methods when applied to malfunctioning rotor systems. In this study, the two types of methods are compared and analyzed. Typical time domain and frequency domain methods are applied and compared under different conditions of excitation. The inherent characteristic that determines the efficiency of the two methods is analyzed. The scopes of application of the two kinds of methods are discussed.

2. Typical Time and Frequency Domain Methods

2.1. The Newmark-β Method and the Fixed Interface Component Mode Synthesis Approach

The Newmark- method is a commonly used time integration method with good convergence and computational stability. Therefore it is chosen here to represent time domain methods.

The equation of motion of a malfunctioning rotor system can be written as

where , , and are the mass, damping (including linear viscous damping and gyroscopic moment), and stiffness matrices, respectively; is the displacement vector; is the excitation vector caused by the imbalance; is the nonlinear force vector, which is a function of and ; the dots above in (1) denote derivatives with respect to time .

Letting , (1) can be changed to

Based on known values of , , at time , the approximate values and are [12]

where , , , , , and .

Let

Equation (3) can then be rewritten as

Assuming that , where is the increment, the above equation is given by

Expanding (2) via a Taylor series around yields

With , (7) can be changed to

where is the instantaneous stiffness matrix. can be obtained using (8) and , can be obtained when convergence criteria are satisfied.

There are four DOFs in one node and if the rotor system is made up of nodes, there are a total of DOFs in the system. The dimension of the instantaneous stiffness matrix is . As increases, the calculation speed decreases rapidly. The fixed interface component mode synthesis approach can be used to reduce the number of dimensions of the system and thus the simulation efficiency can be increased.

The DOFs of the rotor system are separated into two distinct groups. The linear group contains the DOFs that are not related to the nonlinear forces. The nonlinear group contains the DOFs that are related to the nonlinear forces. Equation (2) is rearranged into the following form:

The nonlinear DOFs are fixed and the modes of the system comprised of the remaining DOFs are calculated. The modes that contribute most to the responses of the former rotor system are chosen to construct the fixed interface normal mode set:

Assuming unit displacement in each nonlinear DOF and setting the movement of the linear DOFs to zero, the constrained mode of the nonlinear DOFs can be obtained:

The constrained mode set can then be constructed from all the constrained modes:

The modal matrix of the dimension-reduced system can then be obtained by combining the normal mode set and the constrained mode set:

Using the modal matrix, (9) can be rewritten as:

where , , , and .

The response can be obtained using Newmark method, and the response can then be obtained by .

In the derivation process above, the dimension reduction mainly depends on the reduction of the normal modes. However, this may result in a loss of accuracy.

2.2. Incremental Harmonic Balance Method and Dimension Reduction Using Receptance Data

In the three commonly used frequency domain simulation methods (HB, DF, and IHB), the IHB method is equivalent to the HB method plus the Newton-Raphson method [13]. The IHB method can handle strong nonlinearities and is therefore chosen to represent the frequency domain methods. Here, the multidimensional IHB is introduced because there are many multifrequency excited systems in engineering, such as the dual-rotor system in aeroengines. The detailed information of the single frequency IHB method can be seen in [14, 15].

The equation of motion of a multifrequency excited rotor system is expressed as

where are the exciting frequencies and and are the amplitudes of the exciting forces corresponding to .

When the rotor system is in steady state, the response of the rotor system can be expanded to the form of a multiple Fourier series:

Assuming , , and , (10) can be simplified to

where and .

Let and . Therefore, . Assume that

Substitute (18) into (15) and neglecting the high-order harmonic components, one obtains

where and .

Let and applying the Galerkin’s procedure, (19) yields

Equation (20) can be rewritten as

can be obtained from (21) iteratively, following which can be obtained.

is the instantaneous stiffness matrix and has dimensions . An increase in the DOFs and the number of harmonic components will greatly reduce the computational efficiency of the IHB method.

The receptance data is used in the dimension reduction strategy in the IHB method. The DOFs of the rotor system are separated into two distinct groups. The nonlinear group contains the DOFs that are related to the nonlinear forces and the linear group contains the DOFs that are not related to the nonlinear forces. Rearrange (15) to

When the nonlinear system is in steady state, the response vector is periodic and the nonlinear force vector is also periodic, which can be expressed as

where and .

Assume that , , and . Therefore and satisfies the following equation:

where and .

Rewrite (24) as

Assume that , where is the th order harmonic term of , and . From (25), the th harmonic term of can be written as

Substituting (26) into the first equation of (25) yields

where

and is the th harmonic term of .

For the fundamental harmonic term , assume that , where is the response of the no-rub rotor system. then satisfies the following equation:

From (29), we get

From (27) and (30), we get

The dimensions of (31) are much less than that of the original system and only related to the nonlinear forces. This is because the responses of the linear part are represented by the responses of nonlinear parts using the receptance data. The real part of (31) can be written as

Equation (32) can be solved using the IHB method. After is obtained, can be recovered from using (26).

Since the characteristics of local nonlinearity are utilized, the harmonic terms of all DOFs of the system can be represented by the responses of the DOFs with nonlinearities. Therefore, only the nonlinear DOFs are retained and the computation speed increases greatly. In addition, the accuracy of the results is maintained.

3. Numerical Example and Discussion

3.1. Single Frequency Excitation

A double-span rotor is comprised of two rotors that are connected by a membrane coupling as shown in Figure 1. For the finite element model, the two rotors are divided into beam segments with dimensions that are listed in Table 1. The coupling is represented by a hollow shaft with an outer diameter of 160 mm, an inner diameter of 140 mm, and a length of 500 mm. The hollow shaft is located between nodes 16 and 17.

Table 1: Shaft segment dimension of the two-span rotor system.

Segment
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

(mm)
80
80
290
100
200
240
150
180
160
180
130
100
80
75
75

(mm)
140
180
210
220
235
250
250
300
240
235
230
190
160
140
140

Segment
16
17
18
19
20
21
22
23
24
25
26
27
28
29

(mm)
75
75
62
360
90
390
300
500
300
700
80
360
62
100

(mm)
225
225
270
300
360
420
450
460
450
440
380
310
270
225

Figure 1: The model of the two-span rotor system.

Assuming that the elastic modulus of the steel material is , the supporting stiffness of the left and right rotor is and , respectively. The eccentricities of the left rotor are at nodes 8 and 9, values of  kg·m, respectively. The eccentricities of the right rotor are at nodes 24 and 25, values of  kg·m, respectively. The rotor system is assumed to have proportional damping with a mass damping factor of and a stiffness damping factor of .

Assuming that the full angular rub occurs on node 24 and that the clearance is , the friction coefficient is 0.15.

Figure 2(a) shows a comparison of the time response for a rotational speed of 450 rad/s and a contact stiffness of . The results of the Newmark method without dimension reduction are used as a benchmark. In Figure 2(a), “○” represents the results of the Newmark method with the CMS approach with 5 normal modes, “△” represents the results of the Newmark method with the CMS approach with 10 normal modes, and “×” represents the results of the IHB method with receptance dimension reduction. The time domain methods need a large number of normal modes to attain accuracy targets, whereas frequency domain methods attain very good accuracy more effectively.

(a) Comparison of the time response

(b) Amplitude-frequency response curve

(c) Stable and unstable time responses

Figure 2: The rotor system response under single frequency excitation.

The frequency domain methods combined with the arc-length algorithm [16, 17] can obtain the unstable steady state responses, as shown in Figure 2(b). The curves are the results for contact stiffness values of , , and , respectively. and are the unstable steady state responses. Figure 2(c) shows the stable and unstable time responses for a rotational speed of 550 rad/s.

A comparison of the efficiency of the two methods is listed in Table 2. The computations were performed on a microcomputer operating at 2.83 GHz with 3.25 GB of RAM. The calculation time of the frequency domain method is significantly less than that of the time domain method.

Table 2: CPU time for the two methods under single frequency excitation.

Method
Time (s)

Newmark, CMS order 5
14.08

Newmark, CMS order 10
16.11

Newmark, CMS order 20
19.76

IHB,
0.29

3.2. Double Frequency Excitation

A model of the dual rotor system of an aeroengine is shown in Figure 3 with the dimensions of the two rotors listed in Tables 3 and 4.

Table 3: The dimensions of the inner rotor.

Segment
1
2
3
4
5
6
7

Node
1
2
2
3
3
4
4
5
5
6
6
7
7
8

Length (mm)
150
100
250
1000
250
100
150

Diameter (mm)
100
500
100
100
100
500
100

Table 4: The dimensions of the outer rotor.

Segment
1
2
3
4
5

Node
9
10
10
11
11
12
12
13
13
14

Length (mm)
200
100
400
100
200

External diameter (mm)
200
400
200
400
200

Internal diameter (mm)
150
150
150
150
150

Figure 3: Model of the dual-rotor system.

Bearings I, II, and III are fixed with supports. Bearing IV is an intershaft bearing that connects the outer and inner rotors. The stiffness of the bearings is each. Since there are imbalances in both rotors, the system experiences double frequency excitation. The eccentricity of the inner rotor is  kg·m and mounted on node 4; the eccentricity of the outer rotor is  kg·m and mounted on node 11. The rotor system is assumed to have proportional damping with a mass damping factor of and a stiffness damping factor of . Assuming that the rubbing occurs at node 12 and that the clearance is , the rubbing is partial rubbing with the rubbing angle of 0–90 degrees. The rotating speed of the outer rotor is 1.4 times that of the inner rotor. Assuming the rotating speed of the inner rotor is and that of the outer rotor is , so .

The IHB method is compared with the Newmark method to predict the steady state responses of the dual rotor system. When the IHB method is applied, the harmonic terms are chosen based on the method by Ushida and Chua [18] as the following ranges are chosen for this study: and , . After eliminating the frequencies with negative values, there are all together 25 harmonic terms in the responses.

In Figure 4, “–” represents the results obtained by the IHB method for a rotating speed of , “△” represents the results obtained by the Newmark method for 5 normal modes, and “×” are the results obtained by the Newmark method for 10 normal modes. The results obtained by the IHB method are almost equal to those obtained by the Newmark with a large number of normal modes. However, the Newmark method can only produce time domain results, whereas the IHB method can produce both time domain and frequency results simultaneously. For example, when the IHB method is used to obtain the responses in a given speed range, the amplitude-frequency responses curve and the three-dimensional spectrum diagram can be readily obtained. When the rotating speed varies from 50 rad/s to 350 rad/s and , the amplitude-frequency response curve of node 4 is shown in Figure 4(b). The corresponding three-dimensional spectrum diagram is shown in Figure 4(c) and contains information on each harmonic term.

(a) Comparison of time response

(b) Amplitude-frequency response curve

(c) Three-dimensional spectrum diagram

Figure 4: The responses under double frequency excitation.

The computational time for each method is listed in Table 5.

Table 5: CPU time for the two methods under double frequency excitation.

Method
Time (s)

Newmark, CMS order 5
6.49

Newmark, CMS order 10
8.85

Newmark, CMS order 20
10.48

IHB,
34.7

3.3. Discussion

The simulation results above show that the frequency domain methods are much more efficient than the time domain methods in single frequency excitation condition. This is because of the following.

(1)The frequency domain methods can skip the transient responses. They can obtain the steady state responses directly, while the time domain methods must devote computational time and resources to calculate the transient responses, especially when the damping is small.(2)The size of the instantaneous stiffness matrix is very small. When the number of the nonlinear DOFs and harmonic terms is small, the size of the instantaneous stiffness matrix in frequency domain methods is also small. In contrast, the size of the instantaneous stiffness matrix in time domain methods is mainly determined by the number of normal modes.

The results also show that the efficiency of the IHB method is very high for single frequency excitation but descends quickly for double frequency excitation. This is because of the following.

(1)Calculating the multiple integral for double frequency excitation is time-consuming. The Galerkin integration procedure in the IHB method for the double frequency excitation involves a multiple integral process, which is very time-consuming. For example, for single frequency excitation, one period is divided into 100 subintervals in both the Newmark and IHB method. Therefore, the two methods complete the integration over the same time frame. However, for double frequency excitation, 10,000 subintervals must be created for the multiple integral of the IHB method, which results in much higher time costs than in the Newmark method.(2)The dimensions of instantaneous stiffness matrix are very high in double frequency excitation when using the IHB method. The dimensions of the instantaneous stiffness matrix increase with the increase in the number of harmonic terms, which is much higher for double frequency excitation than for single frequency excitation. For example, when and , , the number of harmonic terms is and the size of the instantaneous stiffness matrix is . This is much larger than the dimensions of for single frequency excitation. For the Newmark method, the dimensions of the instantaneous stiffness matrix remain the same.

If the number of DOFs of the rotor system is larger with fewer nonlinear DOFs, and the damping is small, the frequency domain methods can still have a high efficiency. Otherwise, time domain methods are more efficient.

There are some rotor systems in practice, such as the geared supercharger, that experience three or more than three frequencies excitations. If frequency domain methods are applied, the multiple integral and the large instantaneous stiffness matrix will require too much computational time to determine the response. Therefore, time domain methods are most appropriate in those scenarios.

4. Conclusions

Time domain and frequency domain methods are used to predict the responses of rotor system with malfunctions. The efficiency and accuracy of the two types of methods are compared based on the excitation type. The main conclusions are as follows.

(1)The dimension reduction approach in time domain methods is a substructure approach. The high order modes are truncated, which may affect the accuracy of the methods. Frequency domain methods that use the receptance data for dimension reduction can retain accuracy.(2)The efficiency of the frequency domain methods is determined by the number of excitation frequencies and the number of harmonic terms. These two factors determine the complexity of the numerical integration process and the size of the instantaneous stiffness matrix. Frequency domain methods are highly efficient for single frequency excitations but suffer a decline in efficiency for multiple frequency excitations. The size of the instantaneous stiffness matrix remains the same in time domain methods. Therefore, the efficiency remains the same for both single and multifrequency excitation.(3)Frequency domain methods are recommended for single frequency excitation. Time domain methods are recommended in cases where there are three or more frequency excitations. For double frequency excitations, the choice of simulation method depends on the specific conditions.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors would like to gratefully acknowledge the National Basic Research Program of China (2011CB706504), the Natural Science Foundation of China (51475085), and the Fundamental Research Funds for the Central Universities (N120403007) for the financial support for this study.

References

M. Agnieszka, Rotor Dynamics, CRC Press, New York, NY, USA, 2005.

R. R. Meybodi, A. K. Mohammadi, and F. Bakhtiari-Nejad, “Numerical analysis of a rigid rotor supported by aerodynamic four-lobe journal bearing system with mass unbalance,” Communications in Nonlinear Science and Numerical Simulation, vol. 17, no. 1, pp. 454–471, 2012.

S. Bouaziz, N. B. Messaoud, J.-Y. Choley, M. Maatar, and M. Haddar, “Transient response of a rotor-AMBs system connected by a flexible mechanical coupling,” Mechatronics, vol. 23, no. 6, pp. 573–580, 2013.

C.-C. Wang, “Application of a hybrid method to the nonlinear dynamic analysis of a flexible rotor supported by a spherical gas-lubricated bearing system,” Nonlinear Analysis: Theory, Methods and Applications, vol. 70, no. 5, pp. 2035–2053, 2009.

M. Guskov, J.-J. Sinou, and F. Thouverez, “Multi-dimensional harmonic balance applied to rotor dynamics,” Mechanics Research Communications, vol. 35, no. 8, pp. 537–545, 2008.

F. Wei and G. T. Zheng, “Multiharmonic response analysis of systems with local nonlinearities based on describing functions and linear receptance,” Journal of Vibration and Acoustics, vol. 132, no. 3, Article ID 31004, 2010.

Y. Shen, S. Yang, and X. Liu, “Nonlinear dynamics of a spur gear pair with time-varying stiffness and backlash based on incremental harmonic balance method,” International Journal of Mechanical Sciences, vol. 48, no. 11, pp. 1256–1263, 2006.

S. Kawamura, T. Naito, H. M. Zahid, and H. Minamoto, “Analysis of nonlinear steady state vibration of a multi-degree-of-freedom system using component mode synthesis method,” Applied Acoustics, vol. 69, no. 7, pp. 624–633, 2008.

Y. Ren and C. F. Beards, “A new receptance-based perturbative multi-harmonic balance method for the calculation of the steady state response of non-linear systems,” Journal of Sound and Vibration, vol. 172, no. 5, pp. 593–604, 1994.

P. Bonello and P. Minh Hai, “A receptance harmonic balance technique for the computation of the vibration of a whole aero-engine model with nonlinear bearings,” Journal of Sound and Vibration, vol. 324, no. 1-2, pp. 221–242, 2009.

G. Jacquet-Richardet, M. Torkhani, P. Cartraud et al., “Rotor to stator contacts in turbomachines. Review and application,” Mechanical Systems and Signal Processing, vol. 40, no. 2, pp. 401–420, 2013.

W. Kim, J. Y. Lee, and J. Chung, “Dynamic analysis for a planetary gear with time-varying pressure angles and contact ratios,” Journal of Sound and Vibration, vol. 331, no. 4, pp. 883–901, 2012.

S. H. Chen, Y. K. Cheung, and H. X. Xing, “Nonlinear vibration of plane structures by finite element and incremental harmonic balance method,” Nonlinear Dynamics, vol. 26, no. 1, pp. 87–104, 2001.

S. L. Lau, Y. K. Cheung, and S. Y. Wu, “Incremental harmonic balance method with multiple time scales for nonlinear aperiodic vibrations,” ASME Journal of Applied Mechanics, vol. 50, pp. 871–876, 1983.

J. X. Zhou and L. Zhang, “Incremental harmonic balance method for predicting amplitudes of a multi-d.o.f. non-linear wheel shimmy system with combined Coulomb and quadratic damping,” Journal of Sound and Vibration, vol. 279, no. 1-2, pp. 403–416, 2005.

P. Sundararajan and S. T. Noah, “Dynamics of forced nonlinear systems using shooting/Arc-length continuation method—application to rotor systems,” Journal of Vibration and Acoustics, Transactions of the ASME, vol. 119, no. 1, pp. 9–20, 1997.

J. V. Ferreira and A. L. Serpa, “Application of the arc-length method in nonlinear frequency response,” Journal of Sound and Vibration, vol. 284, no. 1-2, pp. 133–149, 2005.

A. Ushida and L. O. Chua, “Frequency-domain analysis of nonlinear circuits driven by multitone signals,” IEEE Transactions on Circuits and Systems, vol. 31, no. 9, pp. 766–779, 1984.