International Journal of Rotating Machinery

Volume 2014 (2014), Article ID 126843, 8 pages

http://dx.doi.org/10.1155/2014/126843

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

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

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.