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
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.
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
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.
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.
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
The computational time for each method is listed in Table 5.
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.