Misalignment and Rub-Impact Coupling Dynamics of Power Turbine Rotor with Offset Disk

Abstract

When the dual rotor system of the aircraft engine is operating, the mass eccentricity of the power turbine rotor and the misalignment of the shaft coupling or the bearing will cause too large vibration of the rotor; this vibration leads to the rub-impact between the rotor and the casing. The power turbine rotor from the dual rotor system is taken as the research object in this paper. Considering the misalignment, the resulting rub-impact faults, the imbalance of rotor and the disk offset, the equation of motion for the system is developed according to the Lagrangian Equation, and then the Range-Kutta Method is adopted to solve the equation. The influence of the key parameters such as the rotating speed, the misalignment angle and the rub-impact clearance on the dynamics of the system is studied; the finite element analysis was carried out to validate the correctness of the theoretical modeling method. The results show that the rub-impact increases the stiffness of the system; the Hopf bifurcation occurs in the misalignment and rub-impact coupling system; the vibrational stability near the half of the switching speed slumps with the increase of the misalignment angle; with increasing of the stiffness, the number of the chaotic zone increases, and the range of the chaos is widening; enlarging the rub-impact clearance is beneficial to reduce the degree of the rub-impact system and enhance the stability of the system.

1. Introduction

The dual-rotor system of aircraft engines is supported by multiple bearings. The misalignment of the deviation angle for the bearing sleeve will cause periodic changes in the contact angle, clearance and speed of the rolling ball, and these changes enhance the vibration of the system [1]. The study of the dual-rotor system [2] shows that the misalignment fault characteristics of one rotor will be transmitted to another rotor through the intermediate bearing and that the system will resonate when the harmonic frequency caused by the misalignment is the same as the natural frequency of the dual-rotor system. When the system resonates, the sudden increase of the vibrational amplitude will generate the rub-impact between the rotor and the stator, which often results in wear and even the fracture of the blades. Therefore, it is of great significance to study the coupling fault of the misalignment and the rub-impact [3,4] for the safe operation of the engine and the development of aviation technology.

The misalignment of the rotor often leads to the rubbing of the rotor and the stator. Scholars have done a lot of research on the misalignment, the rubbing and their coupling dynamic characteristics. Lu et al. [5] studied the vibrational behaviors of a dual-rotor bearing system with rub-impact coupling faults. Yang et al. [6] adopted an approach to evaluating the reliability of a rotor system. Hou et al. [7] investigated the gap of bearing and Hertzian contact exponent under dynamic loads, and the results show that the dynamic loads have differences of vibrational modes to the system at different rotating speeds. Ma et al. [8,9,10,11] have obtained much progress in the rub-impact of rotor and analytical methods. Chen et al. [12,13] constructed a complicated blade-disk-rotor-bearing coupling model and studied the frictional response characteristic by means of the numeric simulation and the experiments. Prabith et al. [14,15] conducted a multi-disc rubbing-impact analysis of a dual-rotor model of the aero-engine using the approximate time variation method. Xie et al. [16] analyzed the vibrational responses of a dual-rotor system under combined misalignment faults. Liu et al. [17] developed a fault detecting method on the basis of the coupling fault characteristic of a rotor system supported by a sliding bearing. Lu et al. [3] established a mathematical model of a dual-shaft rotor-bearing-gear coupling system to investigate how misalignment impacts the rub-collision. Begg [18] found that when the ratio of rub-impact clearance to structural damping was less than a critical value, the system with full-cycle rubbing was in a stable state. Zhang et al. [19] calculated the stochastic vibrational responses of a rotor system with the rub-impact by using gPCE on the basis of the probabilistic model. Sun et al. [20] studied the effect of the rotor imbalance, the oil film temperature and the rotating speed on the stability of the rotor bearing system under consideration of factors such as the oil film forces, the temperature and the imbalanced forces. The misalignment fault will lead to an increase in the radial displacement of the rotor. When the radial displacement of the disk is greater than the design clearance between the impeller and the casing, the misalignment fault will cause a rub-impact fault. The misalignment-rub-impact coupling fault system often exhibits rub-impact characteristics. In the early stage of rub-impact, the frequency doubling occurs earlier and the amplitude increases faster. The frequency ratio is the index for judging the level of faults [21,22]. Han et al. [23] proposed a semi-supervised learning approach for intelligent fault diagnosis of wind turbines and carried out experiments based on a wind turbine dataset, revealing the efficacy of the method. The research on the vibration response of the coupled fault shows that the influence of misalignment on the dynamic characteristics is concentrated in the high speed region, and the high frequency rubbing is the cause of the unstable state of the rotor [3]. Saeed et al. [24] studied the stability, bifurcation characteristics, and vibration control of a discontinuous dynamic model. Jin et al. [25] analyzed the friction characteristics of the rotor-blade-shell system caused by the reduced asphyxiation of the Hess coupling. The results show that the closer the blade is to the shell, the more serious the friction is. Abouelanouar et al. [26] investigated the effect of misalignment on energy consumption using vibration analysis.

The researchers have carried out much work on the misalignment, the rubbing and the coupling of the two faults of rotors; however, the coupling dynamics of the misalignment, the rub-impact, the imbalance and the disk offset of rotor are rarely studied. In this paper, the power turbine rotor is taken as the research object, considering the misalignment of the bearing and shaft coupling, the resulting rub-impact fault, the imbalance of rotor and the disk offset, the influence mechanism of key parameters such as the rotating speed, the misalignment angle and the rub-impact clearance on the dynamic characteristics of the system is investigated.

2. Misalignment and Rub-Impact Coupling Dynamic Modeling

2.1. Misalignment Model

The misalignment of the two ends of the rotor supports leads to an angle between the driving torque and the shaft, which brings an additional bending moment to the rotor system. The schematic diagram of the misalignment model is shown in Figure 1. The airflow excitation force moment T is horizontally left, the misalignment angle is α; and the angle between the T_s moment and the y axis is β. T can be decomposed into T_a that coincides with the shaft and T_s that is perpendicular to the shaft, and T_s can be decomposed into T_x and T_y. According to the geometric relationship,

$$T_a=T\cos \alpha$$

(1)

$$T_s=T\sin \alpha$$

(2)

$$T_x=T_s\cos \beta =T\sin \alpha \sin \beta$$

(3)

$$T_y=T_s\sin \beta =T\sin \alpha \cos \beta$$

(4)

Supposing that the displacement difference of the bearing misalignment at both ends of the rotor is Δx and Δy, respectively, then α and β are

$$\alpha =arct\frac{\sqrt{\Delta x^2+\Delta y^2}}{l}$$

(5)

$$\beta =arct\frac{\Delta y}{\Delta x}$$

(6)

Since the radial displacement of the rotor at the bearing is much smaller than the shaft length, α and β can be regarded as fixed values. According to Euler’s Equation for the rotation around a fixed point, T_x, T_y, T_z are expressed as

$$\{\begin{array}{c}{T}_x=I_x{\dot{\omega }}_x+{\omega }_y{\omega }_z(T_y-T_z)\\ T_y=I_y{\dot{\omega }}_y+{\omega }_x{\omega }_z(T_z-T_x)\\ T_z=I_z{\dot{\omega }}_z+{\omega }_x{\omega }_y(T_x-T_y)\end{array}$$

(7)

where ω_i is the angular velocity of the rotor; I_i is the principal moments of the inertia of the rotor shaft, i = x, y, z; due to the rotation of the rotor around the z axis, ω_x, ω_y, ${\dot{\omega }}_x$, ${\dot{\omega }}_y$ are all zero. Therefore, Euler’s Equation of rotating around a fixed point can be simplified as

$$T_z=I_z{\epsilon }_z$$

(8)

where ε_z is the rotor angular acceleration. After deriving, the angular velocity and the misalignment angle satisfy the following relationship

$$\omega /{\omega }_d=M/(1+N\cos 2{\theta }_d)$$

(9)

where ω is the rotational angular velocity of the rotor; ω_d is the angular velocity of the air-excited turbine disk; and θ_d is the angular displacement of the turbine disk. Under the small angle misalignment, the turbine disk and the rotor speed are equal, i.e., ω_d ≈ ω. M and N are the coefficients related to the misalignment angle. M = 4cosα/(3 + cos2α), N = (1 − cos2α)/(3 + cos2α). After differentiating Equation (9), one can get

$${\epsilon }_z=\frac{2M\cdot N\sin (2\omega t)}{{\left[1+N\cos (2\omega t)\right]}^2}{\omega }^2$$

(10)

Substituting Equation (10) into Euler’s Equation, one can get

$$T=\frac{2I_z\cdot {\omega }^2\cdot M\cdot N\sin (2\omega t)}{\cos \alpha {\left[1+N\cos (2\omega t)\right]}^2}$$

(11)

By decomposing the misalignment torque, it can be obtained

$$T_y=\frac{2I_z\cdot {\omega }^2\cdot M\cdot N\sin (2\omega t)}{\cos \alpha \cdot {\left[1+N\cos (2\omega t)\right]}^2}\sin \alpha \cos \beta$$

(12)

$$T_x=\frac{2I_z\cdot {\omega }^2\cdot M\cdot N\sin (2\omega t)}{\cos \alpha \cdot {\left[1+N\cos (2\omega t)\right]}^2}\sin \alpha \sin \beta$$

(13)

2.2. Rub-Impact Model

The rub-impact model is shown in Figure 2, where point O, O₁ and O₂ are the geometric center of the casing, the geometric center of the rotor and the mass center of the rotor (material defect and manufacturing tolerances are essentially reflected in the position of the mass center), respectively. The rub-impact force includes the radial collision force P_r and the tangential friction force P_t. The rubbing force conforms to the Coulomb Friction Law, and the rub-impact force can be expressed as [27]

$$\{\begin{array}{l}{P}_r=k(r-\delta )\\ P_t=\mu p_r\end{array}$$

(14)

By decomposing P_r and P_t into x-axis and y-axis, one can get:

$$P_x=\frac{k(r-\delta )}{r}(-x+\mu y)$$

(15)

$$P_y=\frac{k(r-\delta )}{r}(-y-\mu x)$$

(16)

where μ is the friction coefficient; k is the stiffness of the rub-impact, δ is the rub-impact clearance, x and y are the displacements in the x direction and y direction, respectively; and r is the radial displacement of the rotor disk center, $r=\sqrt{x^2+y^2}$.

2.3. Dynamic Model of the Rotor System

The dynamic model of the rotor system is developed based on the following assumption: torsional and axial displacement are ignored and the transverse bending vibration of the rotor system is considered; since the operating temperature of the rotor system is not high enough, the influences of temperature on material properties and system behaviors are small and can be ignored; after estimation, it is found that the gyro effect has a weak influence on the vibration of the rotor. Accordingly, the gyro effect is not considered in this work; the power turbine rotor is the whole forging rotor, the impeller and the shaft are integral, and there is no interaction between the components in the rotor system; the hub and blade are simplified as a rigid disk considering the rotation effect; and the shaft is simplified as a massless elastic shaft. A schematic diagram of the offset rotor system is shown in Figure 3. The mass of the disk is m_d, the eccentricity of the imbalanced disk is e, and the viscous damping of the disk in all directions is c; the clearance between the rotor and the casing is δ; the distance between the left support and the disk is a, the distance between the right support and disk is b, and the length of the shaft is l. These offset parameters have an influence on the vibration characteristics of the system. The dynamic characteristics under the determination offset will be studied in this work; the supporting mass at the left and right ends is m_a and m_b, the supporting stiffness is k_a and k_b, and the supporting viscous damping is c_a and c_b, respectively. In addition, the elastic modulus of the rotating shaft is E, and the moment of inertia of the cross section of the rotating shaft is I; the polar moment of inertia of the disk is J_p, and the moment of inertia of the disk diameter is J_d. The parameter values of the rotor system studied in this work are from a real rotor, as shown in

Table 1. The related parameters of the rotor system.

Parameter	Value
Mass of disk, left and right support md, ma, mb (kg)	34.6, 2, 2
Damping of disk, left and right support c, ca, cb (N·s/m)	2100, 1050, 1050
Support stiffness at left and right end ka, kb (N/m)	1 × 107, 1 × 107
Shaft length and wheel position l, a, b (m)	0.5, 0.5/3, 1/3
Polar and diameter moment of inertia of the disk Jp, Jd (kg·m2)	0.7, 0.35
Wheel eccentricity e (m)	3 × 10−5
Elastic modulus of shaft E (kg·m2)	2.09 × 1011
Moment of inertia of cross section of shaft I (m4)	1.2566 × 10−7

.

The equation of motion for the rotor system is established by using the Lumped Mass Method and the Lagrange Energy Equation. Lagrange equation is

$$\frac{d}{d_t}\left(\frac{\partial L}{\dot{q}}\right)-\frac{\partial L}{\partial q}+\frac{\partial \mathsf{\Phi }}{\partial \dot{q}}=F\left(t\right)$$

(17)

where L = T − U; q represents the displacement of each degree of freedom; and F(t) is the external load of the system, the unbalanced force is considered in this study; compared with the unbalanced force, the influence of the aerodynamic force is negligible. The yaw angles of the wheel around the x axis and the y axis are θ_x and θ_y, respectively; the displacements of the bearing at the left and right ends in the x and y directions are x_a, y_a and x_b, y_b, respectively; the kinetic energy of the rotor system is

$$\begin{array}{l}T=\frac{1}{2}\left[J_d({\dot{\theta }}_x^2+{\dot{\theta }}_y^2)+J_d{\omega }^2\right]-J_p\omega {\dot{\theta }}_y{\theta }_x+\frac{1}{2}{m}_d({\dot{x}}^2+{\dot{y}}^2)\\ \begin{array}{ccc}\begin{array}{cc}& \end{array}& & +\frac{1}{2}{m}_a({\dot{x}}_a^2+{\dot{y}}_a^2)+\frac{1}{2}{m}_b({\dot{x}}_b^2+{\dot{y}}_b^2)\end{array}\end{array}$$

(18)

The generalized coordinates are

$$q=(q_x,q_y)=([\begin{array}{cccc}x& {\theta }_y& x_a& x_b\end{array}],[\begin{array}{cccc}y& {\theta }_x& y_a& y_b\end{array}])$$

(19)

Potential energy of the rotor system

$$U=\frac{1}{2}{q}_x{K}_x{q}_x^T+\frac{1}{2}{q}_y{K}_y{q}_y^T$$

(20)

where K_x and K_y are the stiffness matrices of the o’xz plane and the o’yz plane of the rotor system, respectively.

$$K_x=K_y=\left[\begin{array}{cc}\begin{array}{c}{K}_c\\ -{\mathit{\Phi }}_1^T{K}_c\end{array}& \begin{array}{c}-K_c{\mathit{\Phi }}_1\\ {\mathit{\Phi }}_1^T{K}_c{\mathit{\Phi }}_1\end{array}\end{array}\right]=\left[\begin{array}{cccc}{k}_{11}& k_{12}& k_{13}& k_{14}\\ k_{21}& k_{22}& k_{23}& k_{24}\\ k_{31}& k_{32}& k_{33}& k_{34}\\ k_{41}& k_{42}& k_{43}& k_{44}\end{array}\right]$$

(21)

where K_c is the stiffness matrix considering the axial bending deformation without elastic support, which can be obtained by the Stiffness Influence Coefficient Method and the Flexibility Influence Coefficient Method, and Φ₁ is the biased matrix.

$$K_c=3EIl\left[\begin{array}{c}\begin{array}{cc}\frac{a^2-ab+b^2}{a^3{b}^3}& \frac{a-b}{a^2{b}^2}\end{array}\\ \begin{array}{cc}\begin{array}{cc}& \frac{a-b}{a^2{b}^2}\end{array}& \frac{1}{ab}\end{array}\end{array}\right]$$

(22)

$${\mathit{\Phi }}_1=\left[\begin{array}{cc}\begin{array}{c}\frac{b}{l}\\ -\frac{1}{l}\end{array}& \begin{array}{c}\frac{a}{l}\\ \frac{1}{l}\end{array}\end{array}\right]$$

(23)

The dissipation energy of the disk is

$$\begin{array}{c}{\mathit{\Phi }}_d={\mathit{\Phi }}_{d1}+{\mathit{\Phi }}_{d2}\\ =\frac{\mathrm{l}}{2}c\left[(\dot{x}-\frac{b}{l}{\dot{x}}_a-\frac{a}{l}{\dot{x}}_b)+({\dot{\theta }}_y+\frac{1}{l}{\dot{x}}_a-\frac{1}{l}{\dot{x}}_b)\right]+\frac{\mathrm{l}}{2}c\left[(\dot{y}-\frac{b}{l}{\dot{y}}_a-\frac{a}{l}{\dot{y}}_b)+({\dot{\theta }}_x+\frac{1}{l}{\dot{y}}_a-\frac{1}{l}{\dot{y}}_b)\right]\end{array}$$

(24)

The dissipative function supported at both ends is

$${\mathit{\Phi }}_b=\frac{\mathrm{l}}{2}{c}_a({\dot{x}}_a^2+{\dot{y}}_a^2)+\frac{\mathrm{l}}{2}{c}_b({\dot{x}}_b^2+{\dot{y}}_b^2)$$

(25)

The dissipative energy of the rotor system

$$\mathit{\Phi }={\mathit{\Phi }}_d+{\mathit{\Phi }}_b$$

(26)

The kinetic energy, the potential energy and the dissipation energy of the above rotor system are brought into the Lagrange Equation to obtain the dynamic equation of the coupled rotor system.

$$m_d\ddot{x}+k_{11}x+k_{12}{\theta }_y+k_{13}{x}_a+k_{14}{x}_b+c\dot{x}-\frac{bc}{l}{\dot{x}}_{\mathrm{a}}-\frac{ac}{l}{\dot{x}}_b=p_x+m_{d}e{\omega }^2\cos (\omega t)$$

(27)

$$m_d\ddot{y}+k_{11}y+k_{12}{\theta }_x+k_{13}{y}_a+k_{14}{y}_b+c\dot{y}-\frac{bc}{l}{\dot{y}}_{\mathrm{a}}-\frac{ac}{l}{\dot{y}}_b=m_{d}g+p_y+m_{d}e{\omega }^2\sin (\omega t)$$

(28)

$$J_d{\ddot{\theta }}_y+k_{21}x+k_{22}{\theta }_y+k_{23}{x}_a+k_{24}{x}_b-J_p\omega {\dot{\theta }}_x+c{\dot{\theta }}_y+\frac{c}{l}{\dot{x}}_{\mathrm{a}}-\frac{c}{l}{\dot{x}}_b=T_x$$

(29)

$$J_d{\ddot{\theta }}_x+k_{21}y+k_{22}{\theta }_x+k_{23}{y}_a+k_{24}{y}_b+J_p\omega {\dot{\theta }}_y+c{\dot{\theta }}_x+\frac{c}{l}{\dot{y}}_{\mathrm{a}}-\frac{c}{l}{\dot{y}}_b=T_y$$

(30)

$$m_a{\ddot{x}}_a+k_{31}x+k_{32}{\theta }_y+\left(k_a+k_{33}\right)x_a+k_{34}{x}_b-\frac{b{c}_1}{l}{\dot{x}}_1+\frac{c_1}{l}{\dot{\theta }}_y+\left(c_a+\frac{b^2{c}_1+c_1}{l^2}\right){\dot{x}}_a+\frac{\left(abc-c\right)}{l^2}{\dot{x}}_b=0$$

(31)

$$m_a{\ddot{y}}_a+k_{31}y+k_{32}{\theta }_x+\left(k_a+k_{33}\right)y_a+k_{34}{y}_b-\frac{bc}{l}\dot{y}+\frac{c}{l}{\dot{\theta }}_x+\left(c_a+\frac{b^{2}c+c}{l^2}\right){\dot{y}}_a+\frac{\left(abc-c\right)}{l^2}{\dot{y}}_b=m_{a}g$$

(32)

$$m_b{\ddot{x}}_b+k_{41}x+k_{42}{\theta }_y+k_{43}{x}_a+\left(k_b+k_{44}\right)x_b-\frac{ac}{l}{\dot{x}}_1-\frac{c}{l}{\dot{\theta }}_y+\frac{\left(abc-c\right)}{l^2}{\dot{x}}_a+\left(c_b+\frac{b^{2}c+c}{l^2}\right){\dot{x}}_b=0$$

(33)

$$m_b{\ddot{y}}_b+k_{41}y+k_{42}{\theta }_x+k_{43}{y}_a+\left(k_b+k_{44}{y}_b\right)-\frac{ac}{l}\dot{y}-\frac{c}{l}{\dot{\theta }}_x+\frac{\left(abc-c\right)}{l^2}{\dot{y}}_a+\left(c_b+\frac{b^{2}c+c}{l^2}\right){\dot{y}}_b=m_{b}g$$

(34)

3. Results and Analyses

3.1. Rotating Speed

The rotating speed will change when the rotor system changes the operation mode and conditions or starts and stops. As an important parameter of rotating machinery, rotational speed has a significant impact on the dynamic characteristics of the system. When the rubbing clearance δ = 8 $\times$ 10⁻⁵ m, the rub-impact stiffness K = 2 $\times$ 10⁷ N/m, the friction coefficient μ = 0.1, and the misalignment angle α = 3°, the bifurcation diagram of the vibration amplitude for the rotor with the change of the rotating speed under the misalignment fault is shown in Figure 4. It can be seen from Figure 4 that the vibration of the first-order critical speed occurs in the system, and the critical speed is 499 rad/s. After the speed increases to the critical speed, the amplitude of the vibration jumps from the positive to the negative direction. In addition, due to the small misalignment angle, the vibration at half the critical speed is not excited. The bifurcation diagram under the misalignment-rubbing coupling fault is shown in Figure 5, and the motions of the different speed intervals are shown in

Table 2. The motion types for the different speed ranges.

Speed Range (rad/s)	Motion Type	Speed Range (rad/s)	Motion Type
1~584	Period-1	584~770	Quasi-periodic
770~1018	Period-1	1028~1225	Period-2
1225~1571	Period-1	1571~1726	Period-3
1726~1832	Quasi-periodic	1832~3000	Period-1

. It can be seen from Figure 5 that the system undergoes the period-1 motion, the quasi-period motion, the period-1 motion, the period-2 motion, the period-3 motion, the quasi-periodic motion and the period-1 motion in turn under the coupling fault. By comparing Figure 4 and Figure 5, it can be found that when the misalignment and rubbing faults exist at the same time, the critical speed of the system increases by 54% from 499 rad/s to 770 rad/s, and the Hopf bifurcation occurs; comparing the misalignment fault system, a new quasi-periodic motion of the coupling fault system appears in the range of 1726 rad/s~1832 rad/s.

When the rotating speed of the coupled fault system ω is 400 rad/s, 640 rad/s, 1100 rad/s, 1650 rad/s and 1780 rad/s respectively, the time domain, the spectrum, Poincare section and the axis trajectory diagrams are shown in Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10. It can be seen from Figure 6d that the system has a rub-impact at this rotating speed. Due to the weak rub-impact, there are no other frequency components except the rotational frequency, as shown in Figure 6b. The system is in period-1, which is consistent with Figure 6a,c. The frequency components of 0.25X, 0.5X, 0.75X, 1.5X, 2X and 1.27X in Figure 7b are not conventions. The beat vibration is obvious in Figure 7a,c presents a similar circular curve; the corresponding axis orbit is shown in Figure 7d. It can be seen from Figure 8b that there are two dominant frequency components, 0.5X and 1X, which are in accord with Figure 8c; the amplitude of 0.5X is larger than that of 1X. From Figure 8d, the axis orbit is evolved from the 8-type trajectory; there are two peaks in Figure 8a, which correspond to the two closed circle-like curves in Figure 8d. There is a X/3 frequency division in Figure 9a; the Poincare section is divided into three independent scatter points, as seen in Figure 9c. However, when the rotating speed is 1785 rad/s, the system enters the quasi-periodic motion again. Figure 10c presents a curve like shamrock, and the axis trajectory is not disorderly, as shown in Figure 10d.

3.2. Misalignment Angle

The misalignment of the power turbine rotor system will also affect the vibration response. In order to study the influence of the misalignment of the rotor on the vibration, the misalignment angle α = 5°, 8° and 10° are selected as the control parameter, and the bifurcation and frequency waterfall diagrams are shown in Figure 11, Figure 12 and Figure 13. In Figure 11a, when the rotating speed is 586 rad/s, the system enters the quasi-periodic motion from the period-1 motion, i.e., the state of motion changes, the speed is accordingly called the switching speed, and the switching speeds for α = 8° and 10° are similarly 588 rad/s and 593 rad/s respectively, as shown in Figure 12a and Figure 13a. It can be seen from Figure 11a, Figure 12a and Figure 13a that when the misalignment angle increases, the vibration stability near half of the switching speed slumps, the maximum amplitude of the vibration response decreases, and the range of the quasi-periodic motion widens. Figure 11b, Figure 12b and Figure 13b show that with the increase of the misalignment angle, the maximum amplitude of the resonance decreases slightly.

Chaos has intrinsic randomness, initial value sensitivity, and the characteristics of irregular order; therefore, a variety of tools is applied to judge chaos. Under the misalignment and rub-impact fault, when the speed ω is 600 rad/s, the bifurcation diagram of the amplitude with the misalignment angle is shown in Figure 14, which shows that the rotor system has experienced the quasi-periodic, the period-3, the period-6, the period-12, the period-4 and the chaotic motions. The time domain, the spectrum, Poincare and the axis orbit diagrams under the misalignment angles α = 7°, 9°, 11°, 14° and 16° respectively are shown in Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19. It can be shown from Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19 that with the increase of the misalignment angle, the amplitude of the 2X frequency gradually magnifies up to the first largest peak from the third largest peak in the spectrum diagrams. The other diagrams, such as the time domain, Poincare and the axis orbit, verify this point. The results show that the frequency of 2X is the characteristic frequency of the misalignment fault. When α = 11°, Figure 17c,d shows the system is in chaos, which is so unpredictable that this misalignment angle should be avoided in design.

3.3. Stiffness of Rub-Impact

In the misalignment-rubbing coupling system, the rubbing force is one of the nonlinear sources of the rotor system. As a factor affecting the rubbing force, rubbing stiffness has an important influence on the vibration response. Keeping the rub-impact clearance unchanged, the misalignment angle is 3°, and the bifurcation diagrams and spectrum diagrams under stiffness K = 1.3 $\times$ 10⁷ N/m, 1.5 $\times$ 10⁷ N/m and 1.8 $\times$ 10⁷ N/m are drawn, as shown in Figure 20, Figure 21 and Figure 22. It can be seen from Figure 20a, Figure 21a and Figure 22a that with the increase of the stiffness for the rub-impact, the number of the chaotic zone increases from zero to two; and the range of the chaos is widening, i.e., the larger the stiffness, the larger the unstable quasi-periodic motion interval. In addition, the corresponding critical rotational speeds are 561 rad/s, 567 rad/s and 579 rad/s respectively. This is because increasing the stiffness of the rub-impact will lead to an increase in the critical speed. It is found from Figure 20b, Figure 21b and Figure 22b that with the increase of the rub-impact stiffness, the maximum amplitude of the speed frequency decreases.

The bifurcation diagram of the rub-impact stiffness under the misalignment-rubbing coupling fault is drawn with the misalignment angle = 6°, as shown in Figure 23. It can be seen that with the increase of stiffness, the system goes through the period-1, the quasi-period, the period-8, the quasi-period, the period-5, the quasi-period, the period-2, the period-4, the multi-period, the period-3, the chaotic motions, intermittent inverted bifurcation, paroxysmal bifurcation, and so on.

When the stiffness is less than 2 × 10⁷ N/m, the system is in the period-1 motion because of the weak rub-impact characteristics under the small stiffness; the axis orbit is an ellipse, and there is a point in the Poincare section, as shown in Figure 24. With the increase of stiffness up to 3.5 $\times$ 10⁷ N/m, the rub-impact characteristics are more and more obvious, and the vibration response is more and more complex. The system enters the quasi-periodic motion, and the beat vibration phenomenon is obvious, as shown in Figure 25; when the rubbing stiffness is 5.5 $\times$ 10⁷ N/m, 6.3 $\times$ 10⁷ N/m, 6.8 $\times$ 10⁷ N/m and 7.3 $\times$ 10⁷ N/m, the system is in the period-8, the period-5, the period-2 and the period-4 motion respectively, the axis trajectory is shown in Figure 26; when the rub-impact stiffness increases up to 8.8 $\times$ 10⁷ N/m the system enters chaotic state whose characteristics are shown in Figure 27.

3.4. Clearance

The clearance between the rotor and the stator is an important design parameter of the rotor system. Smaller cleara nce of the rub-impact is beneficial to improve mechanical efficiency, but it also brings the risk of the rub-impact fault. The clearance in the real world can be changed in some cases. To investigate the effect of the clearance on the dynamic characteristics, speed bifurcation and frequency waterfall diagrams under the different clearances are drawn, as shown in Figure 28, Figure 29 and Figure 30. By comparing Figure 28a, Figure 29a and Figure 30a, it can be seen that with the increase of clearance, the intervals of the quasi-periodic motion in the supercritical and critical speed regions decrease, and the maximum amplitude of the vibration decreases. Therefore, increasing the rubbing clearance is beneficial for enhancing the stability of the system. It can be seen from Figure 28b, Figure 29b and Figure 30b that the larger the clearance, the smaller the number of the rub-impact characteristic frequencies, and the weaker the degree of the rub-impact.

The bifurcation diagram with the change in the rub-impact clearance is shown in Figure 31. It can be seen from Figure 31 that with the increase of the rub-impact clearance, the fluctuation range of the amplitude gradually narrows. The research shows that the system undergoes the quasi-period, the period-3, the quasi-period and the period-1 motions in turn. With increasing of the clearance up to 0.15 mm, the rub-impact disappears and there is a jumping phenomenon in the bifurcation diagram from the positive amplitude to negative amplitude.

In order to study the specific nonlinear dynamic characteristics of the system, the rub-impact clearance 3 × 10⁻⁵ m, 6 $\times$ 10⁻⁵ m and 12 $\times$ 10⁻⁵ m are chosen to draw the corresponding time domain, the spectrum, Poincare section and the axis trajectory diagrams, as shown in Figure 32, Figure 33 and Figure 34. It can be seen from Figure 32 that when the rub-impact clearance is 3 $\times$ 10⁻⁵ m, there is an obvious beat vibration phenomenon in Figure 32a, which shows a circular closed curve, indicating that the system is in the quasi-periodic motion state. In Figure 32d, there are annular curves; the amplitude of the frequency 1.35X generated by the rub-impact is significantly higher than that of the frequency 2X generated by the misalignment, as seen in Figure 32b. As the clearance increases, the degree of rub-impact is reduced, as shown in Figure 33.

When the rub-impact clearance is 12 $\times$ 10⁻⁵ m, the vibration response is shown in Figure 34. In Figure 34b, no other frequency components except for the rotating speed frequency and the 2X frequency, and the amplitude of the 2X frequency is far less than that of the rotating speed frequency, indicating that the system is approximatively in period-1, Figure 34c,d verify this point.

4. Verification of the Theoretical Modeling Method

In order to verify the correctness of the modeling approach developed in this paper, the finite element model of the rotor system was established according to the same parameters and constraints as the theoretical modeling based on the Ansys Workbench platform, as shown in Figure 35.

The specific steps of the finite element analysis process are as follows: establishing the geometric model of the rotor system in Creo; importing the model into Ansys Workbench for finite element analysis; defining material property to keep it consistent with the data used in the theoretical calculation; meshing the model of the rotor system. The number of elements is 15,070, the number of nodes is 27,179, the type of the element is tet10, the mean jacobian ratio is 1.039, the mean element quality is 0.729, the mean aspect ratio is 2.27, the mean skewness is 0.407, the mean maximum turning angle is 1.815, these parameters are chosen as the quality check standards; imposing the constraints. The rotor is limited at the two ends by bearings; the bearings are equivalent to springs; the spring stiffness is 1.0 $\times$ 10⁷ N/m which is consistent with the stiffness of theoretical calculation; the axial translation degrees of free (z direction) are limited; the vibrations of x direction and y direction are focused on; applying loads. The displacement response of the rotor vibration is calculated after applying the unbalanced load to the rotor.

After calculating, the time history curve and frequency domain diagram of the vibrational displacement at the center of mass for the disk under unbalanced load are drawn, and a comparison of the Finite Element Method (FEM) result and the Theoretical Modeling Method (TMM) result is carried out, as shown in Figure 36a,b. Because all the model parameters and boundary conditions are consistent with the theoretical calculation, the results of the finite element calculation and the theoretical calculation results are comparable. It can be shown from the figures that the theoretical and finite element calculations agree well, which validates the correctness of the theoretical modeling method.

5. Conclusions

The power turbine rotor from the dual rotor system of the aircraft engine is taken as the research object. Considering the misalignment, the rub-impact, the imbalance and the disk offset of rotor, the equation of motion is established to study the influence of the rotating speed, the misalignment angles and the rub-impact clearance on the dynamic characteristics of the system; the finite element analysis was carried out to validate the correctness of the theoretical modeling method. The main conclusions are as follows:

• The existence of the rub-impact increases the stiffness of the system, thereby increasing the critical speed. Compared with the system of the misalignment fault, the critical speed of the coupling fault system increases by 54%, and Hopf bifurcation occurs in the system.
• The vibration stability near half of the switching speed slumps with the increase of the misalignment angle, the maximum amplitude of the vibration response decreases, the range of the quasi-periodic motion is widening; the frequency 2X is the characteristic frequency of the misalignment. Besides, when α = 11°, the system is in the chaotic state, which is so unpredictable that this misalignment angle should be avoided in design.
• Increasing the rub-impact stiffness reduces the stability of the system. The rotor system with rub-impact is actually a nonlinear system with a non-smooth clearance. Rub-impact stiffness has important effects on the nonlinear characteristics of the system. With the increase of stiffness, the number of the chaotic zone increases; and the range of the chaos is widening; the maximum amplitude of the frequency at the critical speed decreases.
• Increasing the rubbing clearance is beneficial to reduce the degree of the rub-impact and enhance the stability of the system. With an increase in clearance, the interval of the quasi-periodic motion decreases, the maximum amplitude of the vibration decreases; the fluctuation range of the vibration gradually narrows. With an increase in clearance up to 0.15 mm, the rub-impact disappears and there is a reverse jumping phenomenon.