#### Abstract

The asymmetric rotor and the rub-impact behavior are important sources of instability and may cause severe vibrations. However, the dynamics of the rotor-bearing system simultaneously considering the two factors has not gained sufficient attention in available investigations. In this paper, the steady-state response and stability of an asymmetric rotor with rub-impact were evaluated. The asymmetric rotor was modeled by beam elements with asymmetric cross section, and the nonlinear equations of motion were established in the rotating frame. The multiharmonic balance (MHB) method was employed to obtain the linearized form of the nonlinear equations of motion. Either the asymmetry of rotor or rub-impact can result in instability and make the problem difficult to solve. Thus, the arc-length method was utilized to trace the branch of the solutions. In order to improve the calculation speed and accurately predict the solution, the alternating frequency/time domain (AFT) was adopted to calculate the iteration of the arc-length method. Based on the proposed method, the effects of stator stiffness, gap size, unbalance, and asymmetric in shaft on the steady-state response and stability were obtained.

#### 1. Introduction

The rotors with two unequal principal moments of inertia are called asymmetric rotors, which widely exist in rotating machinery, such as two-pole generator rotors and cracked rotors [1, 2]. Besides, the clearance between rotor and stator in rotating machinery has been designed to be smaller and smaller to gain a higher efficiency. However, a smaller clearance makes the rub-impact easier to occur. Both the asymmetry of shaft and rub-impact can result in instability and alter the dynamic behavior. Thus, the investigation on steady-state response and stability of asymmetric rotors with rub-impact are very important to the design and operate the rotating machinery.

The asymmetry of shaft makes the asymmetric rotor-bearing system to become a typical parametric vibration system. The steady-state response and stability of the asymmetric rotor have received much attention [3–6]. If the rotor loses its axial symmetry, there will be unstable region during which the rotor will encounter instability problem [7]. Even if the rotor is well balanced, the stiffness of the rotor changes twice during a rotation period, which results in second-order vibration and resonant peak near half of the critical rotational speed [8]. The system composed of asymmetric rotor and anisotropic bearings is called general rotor [9], whose dynamic behavior is more complicated. The anisotropy of bearings can result in more unstable regions [10] and subcritical resonant peaks [11]. Moreover, the vibration amplitude of the asymmetric rotor-bearing system is usually higher than that of the symmetric one, which arouses higher alternating stress and accelerates the rotor’s damage.

The rub-impact between the rotor and stator is a common fault in rotating machinery and an important issuer of rotor dynamics [12]. The effect of rub-impact on the dynamic behavior have been investigated experimentally and numerically, and some results have been summarized by Muszynska [13], Ahmad [14], and Jacquet-Richardet et al. [15]. The rub-impact in the rotor-bearing system excited by periodic excitation can induce the quasiperiodic except periodic vibration [16, 17]. In order to predict the steady-state response of a rub-impact rotor, Groll and Ewins [18] solved the nonlinear equations of motion in frequency domain by utilizing a harmonic balance method. Besides, the arc-length method was used for searching the equilibrium results. The method was validated by a simple Jeffcott rotor, and the unstable region was determined. Then, Tai et al. [19] employed a similar technique presented by Groll and Ewins [18] to investigate the nonlinear dynamics of a rotor with two disks. Sun et al. [20] used AFT technique to accelerate the convergence of the MHB method. In their study, the nonlinear steady-state response of a dual-rotor system was investigated. This technique was also mentioned in the works of Kim et al. [21, 22], which was used to study the dynamics with other nonlinearities.

Various methods have been developed to study the dynamics of rotor-bearing system with different types of nonlinearity. Specifically, the rub-impact can be described as piecewise linear force with a clearance type nonlinearity [21]. In order to obtain the periodic response of such a nonlinear rotor-bearing system, the equations of motion can be solved in a linearized form by adopting the MHB method [22] or numerically through the time-integration method from the initial conditions [16, 17]. However, when applying the time-integration method, a short integration step is needed if the system has strong nonlinearity, which results in a very large computational effort. Besides, the application of the MHB method increases the dimension of the equation of motion of the nonlinear system, which also increases the computational effort. But it is cost-effective when compared with the time-integration method. The arc-length continuation method is usually adopted after linearization. However, while employing this method, many researchers calculated the derivatives through finite difference approximation [18, 23], which results in a very large computational effort and a loss of accuracy of the numerical evaluation of the derivatives and robustness. To overcome disadvantages of finite difference approximation, Kim et al. [21, 22], Petrov and Ewins [24], Petrov [25, 26], and Sun et al. [20] developed the AFT technique. This technique can guarantee the accuracy of numerical results.

The above literature shows that many researchers have studied the dynamics of asymmetric rotor and rub-impact rotor separately. However, there is little research on the dynamics of the asymmetric rotor with rub-impact. On the one hand, due to time-dependent coefficient of the equations of motion in the fixed frame, the dynamic characteristics of the asymmetric rotor are easier to be analyzed in the rotating frame while the rub-impact forces is usually established in the fixed frame. On the other hand, the rub-impact is of high nonlinearity. When it occurs in an asymmetric rotor, the equations of motion of the asymmetric rotor with rub-impact will become complicated, which is more difficult to solve.

In our research, the steady-state response and stability of an asymmetric rotor with rub-impact were evaluated numerically. The MHB method was used to obtain the linearized form of the nonlinear equations of motion of the system. To guarantee the convergence of the solution, the arc-length method was adopted. Besides, analytical formulations for calculating the vector of nonlinear rub-impact forces and tangent stiffness matrix were presented to accelerate the computation. Then, the effects of stator stiffness, gap size, unbalance, and asymmetry in shaft on the steady-state response and stability were evaluated.

#### 2. Governing Equations of Motion

The equations of motion of an asymmetric rotor with rub-impact can be established both in the fixed frame and the rotating frame. Though most of the equations of motion of rub-impact rotor are in the fixed frame, the asymmetry in shaft introduces time-variant coefficients to all degrees of freedom (DOFs) of the system when the fixed frame is adopted. In order to reduce the difficulty of modeling, the equations of motion were established in the rotating frame. Thus, the coefficients in the equations of motion became time-invariant, and the rub-impact only affects few DOFs of the system, which was easy to deal with.

A schematic diagram of an off-center rotor whose asymmetry was introduced by a shaft with an elliptic section was proposed to obtain the equations of motion (Figure 1). In the rotating frame, the equations of motion of an asymmetric rotor with rub-impact can be written aswhere , , and are the mass, damping, and stiffness matrices of the asymmetric rotor; is the Coriolis matrix; is the spin softening matrix; is the displacement vector in the rotating frame; is the unbalance force vector, which is linear and independent on displacements; and is the nonlinear force vector caused by rotor-stator rub-impact, which is dependent on displacements. By defining and , equation (1) can be written in a simpler form as

In the rotating frame, the direction of the unbalance force is constant. The unbalance vector can be expressed aswhere is the unbalance force vector under unit rotational speed.

The nonlinear rub-impact forces (Figure 2) were regarded as single-sided contact function. The Coulomb-type friction relationship was considered for describing the rub-impact forces. The nonlinear rub-impact forces in the rotating frame can be transformed from that in the fixed frame. The nonlinear rub-impact forces at the disk in the fixed frame can be expressed aswhere the superscript “” represents variables in the fixed frame; and are the normal and the tangential force, respectively; is the horizontal displacement of the disk in the fixed frame; is the stator stiffness; is the gap size between the disk and stator; and is the friction coefficient. In order to obtain the rub-impact forces in the rotating frame, the transformation of displacements between the rotating and fixed frame is introduced aswhere , , , and are displacements in the fixed and rotating frame, respectively. Similarly, the rub-impact forces in the rotating frame can be obtained by the transformation in equation (5):where and are the normal and tangential force in the rotating frame, respectively.

#### 3. Methodology Formulation

##### 3.1. Multiharmonic Balance Method

The asymmetric rotor with rub-impact is excited by periodic excitation. Therefore, the steady-state response is periodic as well, which satisfies . For the rotating machinery, the frequency of excitation force is equal to or multiple times as the rotational speed. The solution of the periodic and steady-state response can be expressed aswhere is the *j*th harmonic components of displacement and is the number of retained harmonic components in the MHB method. In general, the result will become more accurate if the number of retained harmonic components increases. However, the computational costs will increase at the same time. Similarly, the linear forces and nonlinear forces can be written aswhere and are the *j*th harmonic components of the linear forces and nonlinear forces, respectively, and can be calculated by Fourier transform as

Substituting equations (7) and (8) into equation (2) and maintaining the balance of each harmonics, the linearized form of equation (2) can be written aswhere is the residual vector in frequency domain; is the vector that consists of the harmonic components of displacement; and are the vectors that consist of the harmonic components of unbalance forces and nonlinear rub-impact forces of which the form is similar to that of , respectively; and is the dynamic stiffness matrix of the linear part of the system. The dimension of the system in the time domain is . After linearization by the MHB method, the dimension of the system in the frequency domain is and the dimension of the equation increases. The dynamic stiffness matrix can be written aswhere

##### 3.2. Arc-Length Method

The Newton–Raphson method is a common method for solving equation (10), which has quadratic convergence rate. The iteration form of the Newton–Raphson method can be written aswhere is the number of iterations. Convergence aside, however, the Newton–Raphson method fails to accurately predict the steady-state response of system with strong nonlinearity. In order to improve the robustness of the Newton–Raphson method, Riks [27] and Crisfield [28] presented the arc-length method by introducing an additional constrained equation as follows:where is the arc-length and is a scale parameter. When , the corresponding arc-length method is called the spherical arc-length method. and are the increments of displacements and rotational speed, respectively, i.e.,where and are the known solutions of the last iteration.

By applying Taylor series, the linearization form of equations (10) and (14) can be written as

Transform equations (16) and (17) into matrix form, the iteration form of arc-length method can be written aswhere and are the derivatives of the residual vector with respect to displacements and rotational speed.

##### 3.3. Alternating Frequency/Time-Domain Procedure

The calculation of the derivatives of the residual vector with respect to displacement and rotational speed in equation (18) is the most time-consuming process and should be conducted every iteration. However, in the application of arc-length method, the finite difference method is often adopted to obtain the two derivatives, which results in a very large computational effort [18, 23]. To overcome disadvantages of the finite difference method, the two derivatives were solved analytically in this paper to improve the speed and robustness of the arc-length method. First of all, the residual should be written in detail aswhere is the harmonic components of the unbalance load under unit rotational speed, which is constant and independent of displacements or rotational speed. The first derivative can be expressed aswhere is the derivative of the harmonic coefficients of the rub-impact forces with respect to the harmonic coefficients of displacements and can be written as

The derivative is composed of blocks. The block that represents the derivative of the th harmonic coefficient of rub-impact forces with respect to the th harmonic coefficient of displacements can be expressed aswhere is the derivation of the rub-impact forces with respect to the displacements in the rotating frame and . Comparing with equation (9), the calculation of in equation (22) can be treated as Fourier transform. In detail, the expression of can be obtained by transforming that in the fixed frame to the rotating frame, i.e.,where is the derivation of the rub-impact forces with respect to the displacements in the fixed frame.

The nonlinear rub-impact forces only act on disk node. Thus, the nonzero elements of derivatives only exist at the DOFs of the disk node while the rest elements equal to zero. If the th node represents the disk node, the DOFs which the rub-impact forces act on will be numbered from to . The corresponding block of the derivative where the nonzero elements appear can be given by

Then, the derivative which represents the derivative of the residual vector with respect to rotational speed can be expressed as

Both damping and stiffness matrix in the rotating frame are dependent on the rotational speed, so that the dynamic stiffness matrix is also related to the rotational speed. As a result, the derivatives of block of dynamic stiffness matrix with respect to the rotational speed can be expressed as

So far, the derivatives involved in the arc-length method are deduced analytically. Compared with the finite difference method, the speed and accuracy of the computational process are improved. In addition, there are some techniques to improve the speed. Firstly, the fast-Fourier transform (FFT) can be used to obtain the harmonic components of nonlinear rub-impact force and the two derivatives, i.e., and . In practice, the Fourier transformations defined in equations (9) and (22) can be carried out in MATLAB using the FFT package. Secondly, the derivatives in equation (18) can be expressed as the algebraic expression of the rotational speed, and then the coefficients of which can be determined in advance. Finally, for the models with a large amount of DOFs, the nonlinear rub-impact forces only act on limited DOFs. Consequently, when the nonlinear rub-impact forces need to be handled with, only the derivatives at those DOFs should be focused, which can also reduce the computational effort.

##### 3.4. Stability

Once the steady-state response of the asymmetric rotor with rub-impact is obtained, the stability equations can be established by introducing a perturbation motion. Firstly, a small disturbance was assumed near the known steady-state solution , then the vibration after the disturbance was . Substitute the disturbance into equation (2), the stability equations in the time domain can be expressed aswhere the superscript “” represents the derivatives that are calculated for the known solution. Equation (27) describes the free vibration of the nonlinear system. It is obvious that the vibration is periodic with the same period as that of the known solution . The disturbance can be expressed in an exponential function form aswhere is the eigenvector in the time domain. The derivatives of disturbance can be expressed as

Substituting equations (29) and (30) into equation (27), the stability equations can be given by

The above stability equations are in the time domain. The derivatives of nonlinear rub-impact forces with respect to displacement also vary with time, which is difficult to solve. Thus, the stability equations need to be transformed into the frequency domain. First of all, is known as periodic vector so that it can be expressed aswhere is th harmonic components of the eigenvector and is the number of retained harmonic components. Substituting equation (32) into equation (31) and maintaining the balance of each harmonics, the stability equations in the frequency domain can be given bywhere

Eigenvalues and corresponding eigenvectors can be obtained by solving equation (33) in state space. The eigenvalues are complex of which the imaginary parts are the whirl frequencies, and the real parts decide stability. The system is unstable if there is at least one of the eigenvalues with positive real part.

#### 4. Numerical Results and Discussion

##### 4.1. Model Validation

In order to validate the proposed method, a comparison between the results of the proposed method and previous method was presented. The validation model (Figure 1) was discretized by the finite element method. The corresponding finite element model (Figure 3) was composed of twenty beam elements [29] with equal length. The Young’s modulus , density , and Poisson’s ratio of the material of the rotor were 200 GPa, 7800 kg·m^{−3}, and 0.3. The rotor was rigidly supported. The mean area moment of the shaft cross section was 3.068 × 10^{7} m^{4}. The area moments of the shaft cross section along the axes *X* and *Y* were given as and , respectively, where was nondimensional coefficient that was used to parameterize the asymmetry of the rotor. In this research, the effects of stator stiffness, gap size, unbalance, and asymmetry of shaft on the nonlinear dynamics of the system were mainly focused. Therefore, the tangential part of the rub-impact forces was ignored.

In order to validate the propose method, the rotor-bearing system with rub-impact and the asymmetric rotor-bearing system were both investigated separately by conventional methods mentioned in the previous studies. Specifically, the response of the symmetric rotor with rub-impact rotor (Figure 4(a)) was obtained by the proposed method in the rotating frame and the conventional method in the fixed frame, respectively. Meanwhile, the response of asymmetric rotor (Figure 4(b)) was obtained by the proposed method and the method provided by Kang et al. [30], respectively. A good agreement can be observed (Figure 4) between the results of the proposed method and that of the conventional methods, which proves the accuracy of the proposed method. The gray shading speed range in Figure 4 represents the unstable region. Conclusion can be made that either the asymmetry of rotor or rub-impact results in instability. Within the unstable region caused by rub-impact (Figure 4(a)), three different steady-state solutions can be found at the same rotational speed and it seems that the resonant peak shifts to right, which is the so-called “jump phenomenon.” The existence of the asymmetry in shaft can result in two isolated resonant peaks, and there is an unstable region between them (Figure 4(b)). For simplicity, the upper and lower bounds of the unstable region are referred to as inlet and outlet rotational speeds. The points at which the rub-impact starts and ends are defined as starting point and ending point, respectively.

**(a)**

**(b)**

The horizontal amplitudes of disk node versus rotational speed of the asymmetric rotor with rub-impact (Figure 5) were obtained by varying the number of retained harmonics from 1 to 5. When the number of retained harmonics in the MHB method is greater than three, the response amplitudes are hardly changed. However, the computational costs increase as the number of retained harmonics increases. Therefore, in order to maintain the balance between the accuracy and the computational costs, the number of retained harmonics in the MHB method was chosen as five to gain the steady-state response and stability of the asymmetric rotor with rub-impact. The combination of the asymmetry in the shaft and rub-impact makes the two resonant peaks (Figure 4(b)) shift to right, and the unstable region become more complicated (Figure 5). Furthermore, the speed range of unstable region is wider than that caused by either the asymmetric rotor or the rub-impact rotor.

##### 4.2. Effect of Various Parameters

In this section, the stator stiffness, the gap size, the unbalance, and the asymmetry in shaft were chosen as the main factors that influence the nonlinear dynamics of asymmetric rotor with rub-impact. Thus, the effects of these factors on the steady-state response and stability were systematically investigated.

###### 4.2.1. Effect of the Stator Stiffness

The horizontal amplitudes of disk node versus the rotational speed with different stator stiffnesses (Figure 6), i.e., *k*_{n} = 1 × 10^{6} N/m, *k*_{n} = 5 × 10^{6} N/m, and *k*_{n} = 7 × 10^{6} N/m, were discussed in this section. No matter the rotor is symmetric or not, the larger the stiffness of the stator is, the more the resonant peak shifts to right which means the higher is the resonant rotational speed. When the rotor is symmetric with rub-impact (Figure 6(a)), there is no unstable region at *k*_{n} = 1 × 10^{6} N/m while there is one stable region at either *k*_{n} = 5 × 10^{6} N/m or *k*_{n} = 7 × 10^{6} N/m. For the symmetric rotor, if the stator stiffness can beconsidered large enough, the impact of the friction will cause instability, whereas if the stator stiffness is small enough, it will only lead to a slight increase in the resonant rotational speed. The amplitude of resonant response increases with the stator stiffness because the unbalance force is proportional to the square of the rotational speed. As the stator stiffness increases, the outlet rotational speed of the unstable region increases while the inlet rotational speed remains unchanged, which means the range of unstable region is widened by the stator stiffness. For the steady-state response of asymmetric rotor (Figure 6(b)), the “jump phenomenon” happens at both the two resonant peaks near the first critical rotational speed at *k*_{n} = 5 × 10^{6} N/m and *k*_{n} = 7 × 10^{6} N/m. But the “jump phenomenon” cannot be observed at *k*_{n} = 1 × 10^{6} N/m. Different from the symmetric rotor, the larger the stator stiffness, the larger the outlet rotational speed of the unstable region while there is no intuitive law for the inlet rotational speed of the unstable region.

**(a)**

**(b)**

###### 4.2.2. Effect of the Gap Size

The horizontal amplitudes of disk node versus the rotational speed with different gap sizes (Figure 7), i.e., = 0.01 mm, = 0.015 mm, and = 0.02 mm, were discussed in this section. For the symmetric rotor, as the gap size increases, the resonant rotational speed decreases, so does the amplitude of resonant peak because the unbalance force is proportional to the square of the rotational speed. The inlet and outlet rotational speeds of unstable region increase with the decrease of gap size. The reason is that the rub-impact is more likely to appear when the gap size is smaller. In addition, the width of unstable region is not obviously changed by the increase of gap size. For the asymmetric rotor, both the inlet and outlet rotational speeds of the unstable region are influenced by the gap size, which leads to the change of the unstable region. The fact is that the outlet rotational speed of the unstable region is not dominated by the ending point of rub-impact. Moreover, since the stator stiffness is constant, it can be found that the outlet rotational speeds with different gap sizes tend to be the same. The inlet rotational speed of the unstable region decreases as the gap size increases. Therefore, the width of the unstable region is widened by the increase of the gap size.

**(a)**

**(b)**

###### 4.2.3. Effect of the Unbalance

The horizontal amplitudes of disk node versus the rotational speed with different unbalances (Figure 8), i.e., me = 2.5 × 10^{−4} kg·m, me = 1 × 10^{−3} kg·m, and me = 2.5 × 10^{−3} kg·m, were discussed in this section. For the symmetric rotor, the response will be too small to induce rub-impact at me = 2.5 × 10^{−4} kg·m. As the unbalance increases, the starting rotational speed of rub-impact decreases while the ending rotational speed increases. Though the amplitudes of responses are enlarged by the increasing of unbalance, the stability does not get worse with the increase of unbalance. The fact is that there is no unstable region in the symmetric rotor with me = 2.5 × 10^{−3} kg·m. For the asymmetric rotor, the speed range where the rub-impact occurs is widened by the unbalance. The “jump phenomenon” happens near the first resonant peaks at the three investigated unbalances and near the second resonant peaks at me = 2.5 × 10^{−4} kg·m and me = 1 × 10^{−3} kg·m. However, the “jump phenomenon” disappears near the second resonant peak at me = 2.5 × 10^{−3} kg·m since the amplitude of response is enlarged by the unbalance. There is always one unstable region in the asymmetric rotor with rub-impact no matter the magnitude of unbalance. As the unbalance increases, the inlet rotational speed of the unstable region increases while the outlet rotational speed remains unchanged since the stator stiffness is constant, which means the range of unstable region is narrowed by the increase of the unbalance.

**(a)**

**(b)**

###### 4.2.4. Effect of the Asymmetry in Shaft

The horizontal amplitudes of disk node versus the rotational speed with different asymmetries in shaft (Figure 9), i.e., *μ* = 0.4, *μ* = 0.5, and *μ* = 0.6, were discussed in this section. For the asymmetrical rotor without rub-impact, the inlet rotational speed of the unstable region decreases when the asymmetry in shaft increases. Meanwhile, the outlet rotational speed increases. The range of the unstable region is widened by the increase of the asymmetry in the shaft. The “jump phenomenon” happens at the three investigated asymmetries in the shaft. The inlet rotational speed is dominated by the ending point near the first resonant peak and decreases with the asymmetry in the shaft. The outlet rotational speed increases as the asymmetry increases so that the range of the unstable region is widened by the asymmetry in the shaft which is similar to that of asymmetric rotor without rub-impact.

**(a)**

**(b)**

#### 5. Conclusions

In this paper, the effects of the stator stiffness, the gap size, the unbalance, and the asymmetry in shaft on the steady-state response and stability range were investigated by the MHB method combined with the arc-length method. Some conclusions are summarized below:(1)By comparing with the results obtained by conventional methods, the accuracy of the proposed method for obtaining the steady-state response and stability of the asymmetric rotor with rub-impact is proved.(2)The combination of the asymmetry in shaft and rub-impact makes the single resonant peak split into two isolated resonant peaks and both of them shift to right.(3)The increase of stator stiffness will make the resonant peak of the system shift to right. However, only if the stator stiffness is large enough, the rub-impact will lead to “jump phenomenon.” As the stator stiffness increases, the outlet rotational speed of the symmetric and asymmetric rotor increases, while the inlet rotational speed of the symmetric rotor remains unchanged.(4)If the gap size becomes smaller, the rub-impact is more likely to appear. As the gap size increases, the resonant rotational speed of symmetric rotor decreases while that of asymmetric rotor remains unchanged. The increase of gap size widens the range of the unstable region of the asymmetric rotor but has no significant influence on that of symmetric rotor.(5)The unbalance influences the amplitudes of responses directly. The larger the unbalance, the larger resonant rotational speed of the symmetric rotor, but the asymmetric rotor remains unchanged. In addition, the larger unbalance tends to make the system become more stable.(6)No matter there is rub-impact or not, as the asymmetry in shaft increases, the inlet rotational speed of the unstable region of the system decreases and the outlet rotational speed increases. The range of the unstable region expands with the increase of the asymmetry in shaft.

#### Data Availability

The data used to support the findings of this study are included within the article. The processed data are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This work was supported by the research program funded by the 111 Project (B16038), China.