Analysis of Nonlinear Vibration Characteristics and Whirl Behavior of Dual-Rotor Systems with Inter-Shaft Rub Impact

Abstract

Previous studies on rub-impact faults have mainly focused on the rub-impact between rotors and stators, with less research on inter-rotor rub impact. The impact of inter-rotor rub impact on rotor nonlinear vibration is particularly significant. This study investigates the effects of inter-shaft rub impact on the vibration characteristics and whirl behavior of dual-rotor systems. Initially, a dual-rotor model with inter-shaft bearings is established using the finite element method, and inter-shaft rub-impact forces are derived based on contact mechanics. Next, the system response is solved using the Newmark method. Vibration characteristics are analyzed through Campbell diagrams, 3D waterfall plots, time-frequency domain plots, and steady-state rub-impact force plots. Finally, the influence of inter-shaft rub impact on the whirl behavior of the dual-rotor system is studied based on the theory of full-spectrum analysis. The study concludes that inter-shaft rub-impact faults shift the system’s resonance points backward, increase harmonic and combination frequency components, and significantly affect the system response under dual-rotor co-rotation. Excessive friction can lead to self-excited vibrations and sudden amplitude increases, particularly in the LP rotor frequency. Additionally, inter-shaft rub impact primarily affects the whirl behavior of the LP-compressor disk1, showing multiple cycles of forward and backward whirl alternation during acceleration due to combined unbalanced and rub-impact excitations.

1. Introduction

In the structural design of aircraft engines, reducing the clearance between components to pursue higher efficiency and better aerodynamic performance can easily lead to mutual impact and friction between two components due to nonlinear motion such as vibration [1]. This secondary fault is commonly referred to as a rub-impact fault. The occurrence of this fault causes the rotor to exhibit high-intensity nonlinear dynamic behaviors, including quasi-periodic and chaotic responses, sub-harmonic and super-harmonic responses, and reverse whirl, thereby affecting the security and reliability of the engine [2,3]. Therefore, studying the vibration characteristics of the rotor system under rub-impact faults and their effects on whirl behavior is of great significance.

To calculate the nonlinear dynamic response caused by rub-impact faults, it is crucial to establish an accurate rotor dynamic model. Scholars widely apply effective modeling methods, such as the transfer matrix method, lumped parameter method, and finite element method. Shan et al. [4] conducted a simulation study on the nonlinear characteristics of rub-impact forces in rotors. They established the rotor rub-impact dynamic equations using the global transfer coefficient method and obtained the fault characteristics of the rotor under local rub-impact conditions. Zhu et al. [5] established a dynamic model of a dual-rotor-bearing system incorporating rub-impact faults using the lumped parameter method. They investigated the system’s vibration characteristics across three stages: normal operation, minor rubbing, and severe rubbing. Ma et al. [6,7,8] have conducted extensive research on the modeling of complex rotor coupling structures using the finite element method. They have corrected the rub-impact force model between blade structures and casings. The results indicate that blade-tip rubbing can induce high-frequency components, and as the structure becomes more complex, the nonlinear behavior of the system also increases. Recently, some scholars have explored the mechanism of rub-impact faults from a more microscopic level. Yang et al. [9,10,11] have considered the effect of coatings on nonlinear vibration in rotor–stator contact. Through experiments, they have derived a function describing the effect of coatings on the contact stiffness of rotor–stator contact. Compared to the previous assumption of linear contact stiffness, this model more accurately describes the mechanism of nonlinear motion. Building upon this work, Kang et al. [12] considered the influence of coating wear during rubbing. Through finite element simulations and experimental comparisons, they validated the feasibility of the new model. In terms of the impact of rub-impact faults on rotor-whirl behavior, Wang et al. [13] conducted simulation and experimental analyses on a dual-rotor system under single-point rub-impact faults. The study indicated that the concentrated composite frequency components caused by fixed-point rubbing have different weight ratios in forward whirl and reverse whirl. Prabith et al. [14] analyzed the reverse whirl motion of a dual-rotor aircraft engine under multi-disc rub-impact excitations. They observed significant sub-harmonic components in the spectrum. Kang et al. [15,16] first studied the reverse whirl behavior of rotors under unbalanced excitation. They then analyzed the whirl characteristics of a dual-rotor system with single-disk and multi-disk rub impacts. The study found that rubbing can cause significant reverse whirl frequencies, which are particularly influenced by the rub-impact stiffness.

The above are all studies on rub-impact faults between rotors and stators. In dual-rotor aircraft engines, gas leakage is one of the key factors contributing to efficiency loss [17]. To avoid such problems when both rotors are operating in tandem, the gap between the two rotors is typically sealed [18]. This design feature reduces the clearance between the HP and LP rotors. However, if the relative displacement due to the bending vibration of the two rotors exceeds the gap range, rub-impact faults between the shafts may occur. As early as the 1980s, NASA studies pointed out that, under adverse conditions such as blade loss, attention should be paid not only to the rub impact between rotating and stationary components but also to the rub impact between rotors and rotor shafts [19]. He et al. [20] applied spectral analysis to investigate the measured vibration signals of rub-impact faults between HP and LP rotors in aircraft engines. They found that when rotor-to-rotor rub-impact occurs, the signal components are complex, showing significant 0.5× and 2× frequencies in the spectrum. However, the study used a relatively simple rotor model and did not consider the effects of HP rotor vibration. Yu et al. [21,22] found three forms of rub impact during inter-shaft rubbing: continuous rubbing, intermittent rubbing, and self-excited vibration. They analyzed each form and discovered that the dynamic coupling between HP and LP rotors due to inter-shaft rubbing can induce multiple sub-harmonic frequency components. At a certain speed, the self-excited vibration of the rotor is significantly influenced by the initial conditions. Furthermore, they studied the nonlinear vibration behavior of the system under the conditions of HP and LP rotors rotating in the same direction and the opposite direction. They found that the rotation direction significantly affects the inter-shaft rub impact in the rotor system. Under conditions of both rotors rotating in the same direction, the modal components of the system are more complex, and the response amplitude during the self-excited vibration process of the rotor is higher during acceleration. Building on the research by Yu et al., Ling et al. [23] proposed a computational framework suitable for analyzing the response of high-dimensional rotor systems to inter-shaft rub impact. They found that reducing the friction coefficient or rub-impact stiffness can effectively reduce the speed range where self-excited vibrations occur. The aforementioned studies represent the research on the nonlinear vibration characteristics of inter-shaft rub-impact faults. For the diagnosis of inter-shaft rub-impact faults, Bian et al. [24] proposed a parameterized adaptive variational mode decomposition (APVMD) method that effectively extracted characteristic frequencies from the signals of HP rotors when inter-shaft rub-impact occurred.

In previous studies on inter-shaft rub-impact faults, scholars established dual-rotor models where the HP rotor and the LP rotor were treated as independent entities, with dynamic coupling only occurring during rub-impact faults. However, in actual aircraft engines with dual-rotor configurations, the use of inter-shaft bearings between the HP and LP rotors is common. Therefore, building on previous work, this paper first establishes a finite element model of a dual-rotor structure with inter-shaft bearings suitable for studying inter-shaft rub-impact issues. The paper then analyzes the effects of rub-impact stiffness and friction coefficient on the nonlinear dynamic response. Furthermore, there is a gap in the research regarding the impact of inter-shaft rub impact on the whirl behavior in dual-rotor systems. Therefore, this paper, based on the full-spectrum theory, studies the effect of inter-shaft rub-impact faults on the whirl behavior of rotors.

The second section of the paper presents the establishment of a dual-rotor dynamic model using the finite element method, along with the derivation of the dynamic mechanism underlying inter-shaft rub-impact faults. The third section utilizes 3D waterfall plots, time-domain plots, frequency-domain plots, and steady-state rub-impact force plots to study the effects of inter-shaft rub-impact faults on the vibration characteristics of the dual-rotor system. The fourth section, based on the full-spectrum theory, analyzes the impact of inter-shaft rub-impact faults on the whirl behavior of the rotors. The fifth section presents the main conclusions of the entire paper.

2. Modeling of Inter-Shaft Rub-Impact Dynamics in Dual-Rotor Systems

In this section, the establishment of the finite element model for the dual-rotor system and the model for inter-shaft rub-impact forces, as well as the subsequent solution of the fault dynamics equations in the following sections, are all accomplished through the development of MATLAB programs. The simplified model of the dual-rotor structure of an aircraft engine, based on the structure shown in Figure 1, is depicted in Figure 2. This rotor structure consists of a low-pressure rotor, a high-pressure rotor, and bearings.

The LP rotor comprises two disks and two supports, with disk1 and disk2 representing the simplified fan and turbine components, respectively. Similarly, the HP rotor includes two disks and two supports, with disk3 and disk4 representing the simplified compressor and turbine components, respectively. The HP and LP rotors are supported by inter-shaft bearings. Based on reference [25], the specific parameters of the dual-rotor system are given in

Table 1. Structural parameters of components in the dual-rotor system.

Name	Parameter	Value
Disk	Mass md1, md2, md3, md4 (kg) Diameter moment of inertia Jd1, Jd2, Jd3, Jd4 (kg∙m2) Polar moment of inertia Jp1, Jp2, Jp3, Jp4 (kg∙m2)	9.683, 9.139, 9.683, 9.139 0.242, 0.228, 0.242, 0.228 0.484, 0.456, 0.484, 0.456
Shaft length	LP shaft (m) HP shaft (m)	0.6 0.3
Shaft diameter	DLP, (DHPi, DHPo) (m)	0.022, (0.026,0.038), (0.040,0.052)
Shaft material	Density(kg/m3) Young’s modulus (Gpa) Poisson’s ratio	7850 210 0.3
Support	Support stiffness kb1, kb2, kb3, kb4, kb5 (N/m) Support damping coefficient cb1, cb2, cb3, cb4, cb5(N∙s/m)	2.6 × 107,1.75 × 107,1.75 × 107,0.5 × 107,8.75 × 106 1.1 × 103, 1.1 × 103, 1.1 × 103, 1.1 × 103, 1.1 × 103

.

To prevent gas leakage and more closely align with the actual structure of aircraft engines, the front end of the HP rotor has a smaller inner and outer diameter. It is sealed with the LP rotor, and the intermediate section uses a variable-section shaft to transition to the rear end of the HP rotor. The LP rotor is slender and long, making it more flexible and prone to bending deformation. When there is significant relative vibration displacement between the two rotors, areas with small inter-rotor clearances are more likely to experience rub-impact faults, leading to highly nonlinear vibration phenomena.

Based on the finite element method (FEM), the dual-rotor system shown in Figure 2 is divided into multiple simple and uniform finite sub-elements through virtual segmentation. Nodes are defined in each element, and the displacements of these nodes are used to represent the displacements of any micro-element in the element. Through this operation, the continuous system with infinite degrees of freedom is simplified into a system with finite degrees of freedom. The low-pressure rotor is divided into 12 beam elements, each with a length of 50 mm, and 13 nodes in total due to the shared nodes between elements. The high-pressure rotor is discretized into five beam elements with six nodes in total. The lengths of the first, second, fourth, and fifth beam elements are 50 mm, while the third beam element’s length is 100 mm. All disks are simulated using 4-DOF rigid disk elements with mass eccentricity, and all bearing support components are simplified using linear spring-damping elements.

2.1. Finite Element Dynamic Equations

To solve the dynamic response of the entire system, the corresponding dynamic equations need to be derived. The dual-rotor system was discretized using finite element modeling in the preceding text, where the disks, shafts, and bearings were discretized into separate sub-elements. To obtain the dynamic equations for the entire system, it is necessary to derive the dynamic equations for each sub-element and obtain the parameter matrices for each element. These parameter matrices are then assembled in a defined order to form the overall parameter matrix, which is then substituted into the dynamic equations of the entire system.

The first step is to derive the dynamic equations and parameter matrices for the rigid disk element. The kinetic energy equation for the disk element [26] is shown in Equation (1).

$$T_{\mathrm{d}}=\frac{1}{2}{\mathrm{m}}_{\mathrm{d}}({\dot{y}}_d{ }^2+{\dot{z}}_d{ }^2)+\frac{1}{2}{J}_{\mathrm{p}\mathrm{d}}\{{\omega }^2+\omega ({\dot{\varphi }}_{yd}{\varphi }_{zd}-{\varphi }_{yd}{\dot{\varphi }}_{zd})\}+\frac{1}{2}{J}_{\mathrm{d}\mathrm{d}}({\dot{\varphi }}_{yd}{ }^2+{\dot{\varphi }}_{zd}{ }^2)$$

(1)

Here, ${\mathrm{m}}_{\mathrm{d}}$ is the mass of the disk, $J_{\mathrm{d}\mathrm{d}}$ is the moment of the inertia, and $J_{\mathrm{pd}}$ is the polar moment of the inertia. The displacement vector of the disk in the coordinate system is demoted by $u_d={[y_d,z_d,-{\varphi }_{yd},{\varphi }_{zd}]}^T$, where $y_d$ and $z_d$ are the translational degrees of freedom along the Y-axis and Z-axis, respectively, and $-{\varphi }_{yd}$ and ${\varphi }_{zd}$ are the rotational degrees of freedom around the Y-axis and Z-axis, respectively. $\omega$ represents the angular velocity component. Substituting the above equation into the Lagrange Equation (2) [26], we obtain dynamic Equation (3) for the rigid disk as follows:

$$\frac{\mathrm{d}}{\mathrm{dt}}\left(\frac{\partial Td}{\partial \dot{u}d}\right)-\frac{\partial Td}{\partial ud}=Qs$$

(2)

$$[m^{\mathrm{d}}]\left\{{\ddot{u}}_{\mathrm{d}}\right\}+[g^{\mathrm{d}}]\left\{u_{\mathrm{d}}\right\}=\left\{f_{\mathrm{d}}\right\}$$

(3)

The mass matrix $[m^{\mathrm{d}}]$ and gyroscopic matrix $[g^{\mathrm{d}}]$ of the disk are represented by Equation (4). $\left\{f_{\mathrm{d}}\right\}$ is the unbalanced force caused by the eccentricity of the disk element, as shown in Equation (5). Here, $e_d$ represents the eccentricity of the disk, $\mathrm{\Omega }$ represents the current rotational speed of the disk, and ${\phi }_d$ represents the initial phase of the disk.

$$[m^{\mathrm{d}}]=\left[\begin{array}{cccc}{\mathrm{m}}_{\mathrm{d}}& 0& 0& 0\\ 0& {\mathrm{m}}_{\mathrm{d}}& 0& 0\\ 0& 0& J_{\mathrm{d}\mathrm{d}}& 0\\ 0& 0& 0& J_{\mathrm{d}\mathrm{d}}\end{array}\right][g^{\mathrm{d}}]=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& J_{\mathrm{p}\mathrm{d}}\\ 0& 0& -J_{\mathrm{p}\mathrm{d}}& 0\end{array}\right]$$

(4)

$$\left\{f_{\mathrm{d}}\right\}={\left[\begin{array}{llll}{m}_d{e}_d{\mathrm{\Omega }}^2\cos (\mathrm{\Omega }t+{\phi }_d)\hfill & m_d{e}_d{\mathrm{\Omega }}^2\sin (\mathrm{\Omega }t+{\phi }_d)\hfill & 0\hfill & 0\hfill \end{array}\right]}^{\mathrm{T}}$$

(5)

The rotor shaft, as the main component of the dual-rotor system, has a total number of degrees of freedom that mainly depends on the number of finite element discretizations of the shaft. In this paper, three types of Timoshenko beam elements, which consider shear deformation, are used to model the shaft, as shown in Figure 3. Nodes 1 to 13 in the LP shaft are simulated using 12 traditional Timoshenko beam elements, as shown in Figure 3a. In the HP shaft, nodes 14 to 16 and nodes 17 to 19, totaling four beam elements, are constructed using equal cross-section hollow beam elements, as shown in Figure 3b. The variable cross-section beam element shown in Figure 3c is applied between nodes 16 and 17, serving as the transitional section between the front and rear ends of the HP shaft. These three beam elements have different geometric structures but share common assumptions about degrees of freedom. Each element has two nodes, and each node has four degrees of freedom: two radial displacements in the Z and Y directions and two angular displacements around the Z and Y axes. Axial displacement and torsional deformation are ignored. Therefore, each beam element has a total of eight degrees of freedom. The kinetic energy and potential energy of the beam element shown in Figure 3a with an equal section can be determined by integrating the kinetic energy and potential energy of the microelement $\mathrm{d}s$ along the length $l$ of the beam element. Similar to the disk element, the parameter matrices for the beam element can be derived using the Lagrange equation. The detailed derivation process and specific parameter matrices have been clearly outlined in reference [27,28,29].

Before introducing the spring-damping support elements, Rayleigh damping $\mathrm{C}i=\alpha \mathrm{M}i+\beta \mathrm{K}i(i=LP,HP)$ is applied to the rotor part (without damping). This process forms the corresponding damping matrix. Subsequently, support damping is incorporated into the damping matrix of the rotor part based on the position of the support nodes, resulting in the final damping matrix.

To obtain the finite element model of the dual-rotor system, the elements are condensed. The specific condensation method is shown in Figure 4. Taking the stiffness matrix as an example, the LP rotor has 13 nodes and 52 degrees of freedom, while the HP rotor has 6 nodes and 24 degrees of freedom. After condensation, the system has a total of 76 degrees of freedom. The LP-compressor disk is located at node 3, the LP-turbine disk is at node 11, the HP-compressor disk is at node 16, and the HP-turbine disk is at node 17. Nodes 1, 13, 14, and 19 represent regular spring-damping units, and their stiffness is added to the model based on their node positions. The blue part represents the inter-shaft bearing unit located between nodes 9 and 10. This unit couples the independent equations of the LP and HP rotors, introducing coupling terms in the stiffness matrix.

Inserting the condensed overall matrix into the dynamic equations yields the control equation for the entire system, as shown in Equation (6):

$$\mathrm{M}\ddot{\mathrm{q}}+(\mathrm{C}+\left[\begin{array}{cc}\mathrm{\Omega }LP\mathrm{G}LP& 0\\ 0& \mathrm{\Omega }HP\mathrm{G}HP\end{array}\right])\dot{\mathrm{q}}+\mathrm{Kq}=0$$

(6)

In the equation, $\mathrm{M}$, $\mathrm{K}$, and $\mathrm{C}$ represent the mass, stiffness, and damping matrices of the entire system, respectively. $\mathrm{G}LP$ and $\mathrm{G}HP$ represent the gyroscopic matrices of the LP and HP rotors, respectively. $\mathrm{\Omega }LP$ and $\mathrm{\Omega }HP$ are the rotational speeds of the LP and HP rotors, respectively. $\mathrm{q}$, $\dot{\mathrm{q}}$, and $\ddot{\mathrm{q}}$ represent the displacement column vector, velocity column vector, and acceleration column vector of the entire rotor system, respectively. If solving for the dynamic response of the system, $\mathrm{q}$, $\dot{\mathrm{q}}$, and $\ddot{\mathrm{q}}$ will encompass the displacements, velocities, and accelerations of all nodes in all degrees of freedom. It is important to emphasize that, because rigid disk elements are used in this study, considering only the mass and gyroscopic moment of the disk without considering deformation. Therefore, the displacement of the disk is equivalent to the displacement of the nodes on the shaft.

2.2. Inter-Shaft Rub-Impact Force Model

Due to the effect of disk imbalance forces, the system exhibits whirl phenomena. When the relative whirl amplitude between the HP and LP rotors exceeds the gap between them, inter-shaft rub-impact faults occur. In the aforementioned study, the finite element method discretizes the rotor structure into multiple elements, with each beam element having two nodes. When rubbing occurs, there is an interacting excitation force, which can be considered as acting on node 5 of the LP shaft and an equal and opposite reaction force acting on node 15 of the HP shaft.

Figure 5 illustrates the stable operation of the high and LP rotors and the occurrence of inter-shaft rub impact, where the two shafts come into contact with each other. In the stable operating state, the initial clearance between the two rotors is $\delta 0$, where $RLP$ and $RHP$ are the radii of the LP rotor and the inner radius of the HP rotor, respectively. $\mathrm{\Omega }LP$ and $\mathrm{\Omega }HP$ are the rotational speeds of the LP and HP rotors, respectively. When rub-impact occurs, $r_{LP}$ and $r_{HP}$ are the precession motion radii of the LP rotor and the HP rotor, respectively. In this paper, these two quantities are vectors, as shown in Equation (7). $z_{LP}$ and $y_{LP}$ are the displacements of the rub point of the LP rotor in the z and y directions, respectively, while $z_{HP}$ and $y_{HP}$ are the displacements of the rub point of the LP rotor in the z and y directions, respectively. $\mathsf{\omega }LP$ and $\mathsf{\omega }HP$ are their respective precession angular velocities. To accurately depict the kinematic mechanism under this condition, the linear model and Coulomb friction model are used to simulate the variation of forces during the rubbing process, leading to the derivation of Equation (8).

$$\begin{array}{l}\left|r_{LP}\right|=\sqrt{{z_{LP}}^2+{y_{LP}}^2}\\ \left|r_{HP}\right|=\sqrt{{z_{HP}}^2+{y_{HP}}^2}\end{array}$$

(7)

$$\begin{array}{l}\{\begin{array}{l}FnLP=H(\left|rLP-rHP\right|-{\delta }_0)\times krub(\left|rLP-rHP\right|-{\delta }_0)\\ FtLP=sign(vrel)\times \mu FnLP\end{array}\\ \{\begin{array}{l}FnHP=-FnLP\\ FtHP=-FtLP\end{array}\\ sign(vrel)=\{\begin{array}{c}-1\\ 0\\ 1\end{array}\begin{array}{c}vrel<0\\ vrel=0\\ vrel>0\end{array}\end{array}$$

(8)

$$vrel=(\mathrm{\Omega }LPRLP+\omega LPrLP)-(\mathrm{\Omega }HPRHP+\omega HPrHP)$$

(9)

The formula describes the normal contact force experienced by the LP rotor during inter-shaft rub impact, where $FnLP$ is the normal rub-impact force, $FtLP$ is the tangential rub-impact force, $FnHP$, and $FtHP$ are the forces acting on the HP and LP rotors, respectively, with equal magnitudes but opposite directions. $krub$ is the rub-impact stiffness. $\mu$ is the Coulomb friction coefficient. $H(\cdot )$ is the Heaviside function used to indicate whether rub impact occurs. When $\left|rLP-rHP\right|-\delta 0>0$, rub impact occurs. And $H(\left|rLP-rHP\right|-\delta 0)=1$, otherwise $H(\left|rLP-rHP\right|-\delta 0)=0$. $sign(\cdot )$ is a sign function, representing the direction of the friction force. Here, $vreL$ is the relative velocity at the contact point, considering both the rotor precession speed and rotor rotation speed. This is detailed in Equation (9).

The change in the direction of the friction force directly impacts the positive or negative work done on the whirl motion. If the LP rotor undergoes reverse whirl motion, and the amplitude is large, when the direction of the friction force is the same as the direction of the whirl motion, it increases the rotor’s energy and reduces the damping dissipation capacity, which may lead to system instability. Decomposing the normal force and tangential force in the radial degree of freedom direction from Figure 5 yields the following Equation (10):

$$\begin{array}{l}\{\begin{array}{l}frLPz=H(\left|rLP-rHP\right|-\delta 0)\cdot krub\left(1-\frac{\delta 0}{\left|rLP-rHP\right|}\right)\times (-(zLP-zHP)+sign(vrel)\cdot \mu (yLP-yHP))\\ frLPy=H(\left|rLP-rHP\right|-\delta 0)\cdot krub\left(1-\frac{\delta 0}{\left|rLP-rHP\right|}\right)\times (-sign(vrel)\cdot \mu (zLP-zHP)-(yLP-yHP)\end{array}\\ \{\begin{array}{l}frHPz=-frLPz\\ frHPy=-frLPy\end{array}\end{array}$$

(10)

The terms $frLPz$ and $frLPy$ represent the forces acting on node 5 of the LP rotor in the Y-axis and Z-axis directions, while $frHPz$ and $frHPy$ correspond to the forces acting on node 15 of the HP rotor in the Y-axis and Z-axis directions.

The overall dynamic equation of the dual-rotor system under inter-shaft rub impact is derived by incorporating the unbalanced forces and rub-impact forces into the right-hand side of the dynamic equation based on their respective force node positions. This yields Equation (11):

$$\mathrm{M}\ddot{\mathrm{q}}+(\mathrm{C}+\left[\begin{array}{cc}\mathrm{\Omega }LP\mathrm{G}LP& 0\\ 0& \mathrm{\Omega }HP\mathrm{G}HP\end{array}\right])\dot{\mathrm{q}}+\mathrm{Kq}=\mathrm{F}d+\mathrm{F}rub$$

(11)

where the $\mathrm{M}$, $\mathrm{K}$, $\mathrm{C}$, $\mathrm{G}LP$, and $\mathrm{G}HP$ matrices correspond to those in Equation (5), ${\mathrm{F}}_d$ represents the column vector of unbalanced forces of the disk, and $\mathrm{F}rub$ represents the column vector of inter-shaft rub-impact forces.

3. Analysis of Vibration Characteristics of Inter-Shaft Rub-Impact Fault in Dual-Rotor System

3.1. Natural Characteristics of the Dual-Rotor System

Figure 6 shows the Campbell diagram of the dual-rotor system. In the diagram, the red line represents the LP synchronous excitation line, the purple line represents the HP synchronous excitation line, the blue lines represent the various orders of reverse whirl frequency lines, and the yellow lines represent the various orders of forward whirl frequency lines. The intersection of the excitation lines and the dynamic frequency lines corresponds to the system’s critical speeds. Define the speed ratio between the HP and LP rotors as $\lambda =\mathrm{\Omega }HP/\mathrm{\Omega }LP$. When the HP and LP rotors rotate in the same direction, $\lambda >0$. When they rotate in opposite directions, $\lambda <0$. In this paper, it is assumed that $\lambda =1.5$ under the condition of same-direction rotation and $\lambda =-1.5$ under the condition of opposite-direction rotation. The critical speed values vary under different rotation conditions, indicating an influence of the rotor system’s mode when there is an inter-shaft bearing between the two rotors. Figure 6a depicts the Campbell diagram of the system when the two rotors rotate in the same direction. The first critical speed point, labeled as A, is caused by the unbalance of the HP disk, while the first critical speed point B is caused by the unbalance of the LP disk. Similarly, the second critical speed points, labeled as C and D, are caused by the unbalance of the HP and LP disks, respectively. The third critical speed points, labeled E and F, are also caused by the unbalance of the HP and LP disks, respectively, as shown in Figure 6b and similar diagrams. The specific critical speed data are presented in

Table 2. Critical speeds of the dual-rotor system under different rotation conditions.

Order	Point	Critical Speeds of Dual-Rotor System (rpm)		Excitation Source
		Co-Rotation	Counter-Rotation
1	A	$\mathrm{\Omega }LP=3489$	$\mathrm{\Omega }LP=3445$	HP rotor
	B	$\mathrm{\Omega }LP=5410$	$\mathrm{\Omega }LP=5299$	LP rotor
2	C	$\mathrm{\Omega }LP=4527$	$\mathrm{\Omega }LP=4411$	HP rotor
	D	$\mathrm{\Omega }LP=6872$	$\mathrm{\Omega }LP=6610$	LP rotor
3	E	$\mathrm{\Omega }LP=7425$	$\mathrm{\Omega }LP=7454$	HP rotor
	F	$\mathrm{\Omega }LP=\mathrm{11,194}$	$\mathrm{\Omega }LP=\mathrm{11,278}$	LP rotor

.

3.2. Analysis of the Influence of Inter-Shaft Rub-Impact Parameters on System Vibration Characteristics

This section investigates the nonlinear vibration response of a dual-rotor system with inter-shaft rub impact. The dynamic equations of the system are derived, and numerical simulations of the rub-impact process are conducted using the Newmark-β explicit numerical integration method. In the calculation, the integration step is chosen as $\mathrm{dt}={10}^{-4}$ s and the total integration set as 3 s, which ensures the system response is in a steady state. During the analysis, steady-state data were selected to analyze the response characteristics. These data were complemented by the rotor structural parameters shown in Table 1. The inter-rotor gap size was set to $\delta 0=0.01$ mm, the disk eccentricity distance was set to $e=0.0025$ mm, and both unbalanced phases were set to $0^{\circ }$.

In this model, the HP-compressor disk (disk3) is located close to the rubbing position, leading to a more pronounced response that can significantly reflect the system’s dynamic behavior under fault conditions. The magnitude of the rubbing excitation directly determines the amplitude of the system’s response. According to the derived rubbing force Formula (6) in Section 2.2, the magnitude of the normal collision force is mainly related to the rubbing stiffness $krub$, while the magnitude of the tangential friction force is mainly influenced by the friction coefficient $\mu$. This section will analyze the effects of these two parameters on the system’s vibration characteristics.

3.2.1. Influence of Rubbing Stiffness on System

To investigate the influence of rubbing stiffness on the vibration characteristics, three different rubbing stiffness values, $krub=5e6$, $krub=7e6$, and $krub=9e6$ are selected to calculate the system response during the acceleration process, with a friction coefficient of $\mu =0.3$. Figure 7 and Figure 8, respectively, depict the waterfall plots of the response at the disk3 position under the influence of rubbing stiffness for the $\lambda =1.5$ and $\lambda =-1.5$ conditions. When the speed ratio is $\lambda =1.5$, it can be observed that there are two prominent frequency amplitudes in all three operating conditions, $f_L$ and $f_H$, which are the rotational frequencies of the LP and HP rotors. In the absence of rubbing, the system exhibits only these two frequency components. The amplitude at the $f_H$ reaches its peak near the third critical speed of $\mathrm{\Omega }LP=7425$ rpm caused by the HP rotor. When the rub-impact fault occurs, with the increase of the rub-impact stiffness, the peak amplitude shifts noticeably to higher speeds and decreases. Considering the fault as part of the system structure, increasing the rub-impact stiffness is equivalent to increasing the system stiffness, resulting in the shift of the critical speed. It also has a certain inhibitory effect on high-frequency vibration displacement. In addition to the dual-rotor rotation excitation frequency, there are also a small number of other harmonic and combination frequency components. These frequencies are mainly multiples of the LP rotor excitation, as indicated in Figure 7a. As the rubbing stiffness increases, the amplitudes of the other harmonic components increase, except for the excitation frequency of the HP rotor, with the $f_H-f_L$ component showing the most significant increase. Figure 7d shows the waterfall plot of the displacement of disk3 as the rub-impact stiffness changes when $\mathrm{\Omega }LP=7680$. From this plot, it is more evident that $f_H$ decreases as the rub-impact stiffness increases. When the speed ratio is $\lambda =-1.5$, $f_L$ and $f_H$ still dominate. This is the same as when $\lambda =1.5$. Increasing the rubbing stiffness leads to a decrease in the amplitude of $f_H$, and the same applies to other frequency components. The difference lies in the fact that, under counter-rotation conditions, the overall amplitude of $f_L$ is lower than that under co-rotation conditions, while $f_H$ is generally higher, indicating more severe vibration of the HP rotor. The amplitudes of other harmonic and combination frequency components also increase with an increasing rubbing stiffness, but the peak amplitude of $f_H-f_L$ shifts noticeably and increases significantly compared to the co-rotation condition.

3.2.2. Impact of Friction Coefficient on the System

Compared to the rubbing stiffness, the friction coefficient has a more significant impact on the system. Figure 9 and Figure 10 show the waterfall plots of the system response at the disk3 position under conditions of $\lambda =1.5$ and $\lambda =-1.5$, respectively, with different rubbing friction coefficients of $\mu =0.5$ and $\mu =0.7$ and a stiffness of $krub=5e6$ in comparison with Figure 7a and Figure 8a. Under the condition of co−rotation, when $\mu =0.5$, the waterfall plot shows a more complex combination of frequency components, such as $2f_H+f_L$ and $f_H+3f_L$, indicating that, with the increase of the friction coefficient, the coupled vibration of the dual rotors becomes more diverse. Throughout the entire speed-increase process, there has been a vibration amplitude at a frequency of 120 Hz, which is caused by self-excited vibration between the shafts. When the rotor speed reaches the self-excited vibration frequency, resonance occurs, resulting in a significant vibration-displacement response. As shown in Figure 8a, when the LP rotor’s rotational frequency $f_L$ reaches 120 Hz, which corresponds to a speed of 7200 rpm, the vibration amplitude is very high, exceeding the amplitude of $f_H$ throughout the entire acceleration process. This frequency is close to the second critical speed caused by the LP rotor, as analyzed earlier. It indicates that the excitation energy of the LP rotor can be transmitted to the HP rotor through inter-shaft rubbing and inter-shaft bearing. When $\mu =0.7$, in the range of 7300–7500 rpm, the region of self-excited vibration expands. The positive work done by the friction force on the rotor may surpass the energy dissipation produced by the rotor’s original damping, resulting in a sharp increase in vibration amplitude and a high risk of instability. In other speed ranges, the response is similar to that when $\mu =0.5$, as shown in the yellow and green waterfall plots in Figure 8b. With an increase in the friction coefficient, the amplitude of other frequency components also increases accordingly. Under counter-rotation conditions, as the friction coefficient increases, the amplitude of each frequency changes, but not as significantly as under co-rotation conditions. There is no occurrence of self-excited vibration during the increase in speed, thus avoiding excessive vibration amplitudes caused by resonance. This indicates that, when the HP rotor rotates much faster than the LP rotor during counter-rotation, the relative velocity at the contact point is high, and the friction force always does negative work on both rotors, making the system less prone to instability. Therefore, overall, with the increase in the friction coefficient, the stability of the dual-rotor system under counter-rotation conditions is significantly better than under co-rotation conditions.

To further investigate the influence of the friction coefficient on the vibration response characteristics of a co-rotating dual-rotor system, a comparative analysis was conducted based on the time-domain waveforms, frequency spectra, steady-state rub-impact force plots, and relative speed at the rub-impact point. Figure 11 shows the rub-impact response at the disk2 and disk3 positions under co-rotation conditions at a fixed speed of 7200 rpm, with a rub-impact stiffness of $krub=5e6$ and friction coefficients of $\mu =0.3$, $\mu =0.5$, and $\mu =0.7$. It also depicts the rub-impact force acting on the LP rotor and the relative speed at the rub-impact point. From the spectrum, it can be observed that the frequency components of the HP and LP rotors are generally the same, with differences in the amplitude. With an increase in the friction coefficient, the vibration amplitude of the LP rotor increases significantly, whether observed at the LP or HP rotor. There are accompanying frequencies around each order, mainly caused by resonance with self-excited vibration. The time-domain waveform shows clear impact signals, indicating the system’s instability. The plot of the rubbing force vividly illustrates this instability. At 7200 rpm, the friction force acting on the LP rotor is negative. According to the derivation of the rubbing-force model in Section 2, when the friction force is positive, it does negative work on the rotor motion. Therefore, in this operating condition, the friction force does positive work on the LP rotor, which is consistent with the above analysis. When $\mu =0.3$, both the rubbing normal force and rubbing tangential force vary periodically over time, with periods of zero values. This indicates that rubbing is intermittent, with the rotor repeatedly contacting and separating. When $\mu =0.5$ and $\mu =0.7$, the rubbing occurs intermittently within the cycle, accompanied by sustained rubbing.

4. Analysis of Inter-Shaft Rub-Impact Fault Whirl Behavior Based on Full-Spectrum Method

Rotor whirl is the combination of the rotor’s self-rotation and precession. When the precession direction is the same as the rotation direction, it is called forward whirl. Otherwise, it is called backward whirl. When the rotor undergoes backward whirl, it can lead to more severe system vibration and decreased stability. Therefore, investigating the direction of the whirl during inter-shaft rub-impact faults is particularly important. This paper applies the full-spectrum method to study the whirl behavior of a dual-rotor system with inter-shaft rub-impact faults. The full-spectrum method was first proposed by Bently. For a detailed introduction to the theory and derivation process of the full-spectrum method, please refer to [17].

In the full-spectrum plot, for the LP rotor, the frequency component on the positive half-frequency axis corresponds to the forward precession component, while the frequency component on the negative half-frequency axis corresponds to the backward precession component. For the HP rotor, when both rotors rotate in the same direction, the determination method is the same as for the LP rotor; when both rotors rotate in opposite directions, the determination method is opposite. By comparing the amplitudes of the forward and backward precession components under the main frequency component, the rotor’s whirl direction can be determined.

4.1. Analysis of Whirl Behavior in Unbalanced Dual-Rotor Systems

To verify the feasibility of this method, a whirl direction analysis is conducted on a system with unbalanced excitation only. Figure 12 shows the full-spectrum waterfall plot of the LP-turbine disk2 for the high-low rotor-speed ratio of $\lambda =1.5$ and $\lambda =-1.5$. When there is no rub impact, the frequency components in Figure 11 are relatively simple, with N_L and N_H representing the rotational frequencies of the LP and HP rotors, respectively. When $\lambda =1.5$, NL and NH exist only on the positive half-frequency axis, and there are no frequency components on the negative half-frequency axis. This indicates that, throughout the rotor acceleration process, disk2 maintains a state of prograde precession. When $\lambda =-1.5$, NL is located on the positive half-frequency axis, and NH is located on the negative half-frequency axis. When the LP rotor speed is $\mathrm{\Omega }LP=9600$ rpm, the amplitude of NL exceeds that of NH. This indicates that, before 9600 rpm, the precession direction of disk2 is mainly influenced by the unbalanced excitation of the HP rotor, showing retrograde precession. After the speed reaches 9600 rpm, the precession direction of disk2 changes to prograde precession.

To further verify the accuracy of the above rules, significant rotation-speed conditions were selected to analyze the LP-turbine disk2 at constant rotation speeds. Figure 13 shows the fixed rotation-speed spectrum of disk2 at $\mathrm{\Omega }LP=7200$ rpm and $\mathrm{\Omega }LP=9600$ rpm under the co-rotation condition for the dual-rotor system. Figure 14 shows the orbit diagram of disk2 at $\mathrm{\Omega }LP=7200$ rpm and $\mathrm{\Omega }LP=9600$ rpm under the co-rotation condition for the dual-rotor system. The conclusion is consistent with the previous analysis. When $\lambda =1.5$, the rotor undergoes a forward whirl at a constant speed, but with different magnitudes of the main frequency component. This indicates a change in the dominant excitation, causing the rotor to undergo a forward whirl at $\mathrm{\Omega }LP=9600$ rpm, mainly due to NL. In the orbit diagram, red and green dots are marked on the trajectory, corresponding to the axial motion positions of the center at time t and t + 2dt, respectively, where dt is the calculation step size. The position of the marked points indicates that the center of disk2 moves counterclockwise, which is in the same direction as the LP rotor’s rotation, indicating a forward whirl.

Figure 15 shows the fixed rotation-speed spectrum of disk2 at $\mathrm{\Omega }LP=7200$ rpm and $\mathrm{\Omega }LP=9600$ rpm under the counter-rotation condition for the dual-rotor system. Figure 16 shows the orbit diagram of disk2 at $\mathrm{\Omega }LP=7200$ rpm and $\mathrm{\Omega }LP=9600$ rpm under the counter-rotation condition for the dual-rotor system. By comparing the full spectrum and orbit diagram at $\mathrm{\Omega }LP=7200$ rpm, it is found that the center of disk2 moves clockwise, which is opposite to the direction of rotation, indicating a reverse whirl. At $\mathrm{\Omega }LP=9600$ rpm, disk2 returns to a forward whirl. The conclusion of this study is consistent with the previous research by Kang et al. [16], confirming the effectiveness of this method.

4.2. Analysis of Whirl Behavior in Dual-Rotor Systems with Inter-Shaft Rub-Impact Faults

In aircraft engines, counter-rotating dual-rotor configurations significantly improve rotor gyroscopic torque effects and enhance maneuvering flight performance. They also allow for a reduction in the number of guide vanes in the LP turbine, theoretically enabling the elimination of the LP-turbine guide apparatus, thus reducing the overall weight of the engine. Therefore, counter-rotating dual-rotor engine technology is an important trend in the development of aircraft engines. This section focuses on studying the effect of fault excitation on whirl behavior in counter-rotating dual-rotor systems when inter-shaft rub-impact faults occur between the rotors. As is well known, inter-shaft rub-impact faults between rotors can lead to a reverse whirl of the rotor. Kang et al. [17] found in their study that reverse whirl caused by rotor–stator rub impact includes not only frequency components due to unbalanced excitation but also significantly large negative reverse whirl frequencies. However, in this study, the inter-shaft rub-impact excitation does not directly act on the disk position. The whirl of the disk is likely to be predominantly influenced by the unbalanced excitation rather than the inter-shaft rub-impact excitation. This paper verifies this using simulation methods.

Figure 17 shows the full-spectrum waterfall plots of the LP-compressor disk (disk1), LP-turbine disk (disk2), HP-compressor disk (disk3), and HP-turbine disk (disk4) when inter-shaft rub impact occurs. The speed ratio is $\lambda =1.5$, the rub-impact stiffness is $krub=5e6$, and the friction coefficient is $\mu =0.3$. The occurrence of rub-impact faults will generate other sub-harmonics and harmonics, as verified in Section 3. Upon conducting a full-spectrum analysis of the signals, it can be observed that frequency components caused by rubbing are present on both the positive and negative semi-frequency axes. As shown in Figure 16a, $f_{\mathrm{bw}}$ and $f_{\mathrm{fw}}$, respectively, represent the frequencies of backward whirl and forward whirl caused by rubbing. In Figure 16b,c, the amplitudes of $f_{\mathrm{bw}}$ and $f_{\mathrm{fw}}$ are relatively small, indicating that the excitation of the disk’s whirl motion at this location is still mainly dominated by the unbalanced excitation. The other frequency components only have a minor influence on the whirl direction. However, in Figure 16a, the LP-compressor disk (disk1) exhibits multiple changes in whirl direction throughout the acceleration process. Before $\mathrm{\Omega }LP=5280$ rpm, the amplitude of NH on the negative half-frequency axis is higher than the amplitude of NL on the positive half-frequency axis, indicating that disk1 is undergoing reverse whirl. During the acceleration from 5760 rpm to 7680 rpm, the unbalanced excitation of the LP rotor dominates, causing disk1 to whirl in the positive direction. Near $\mathrm{\Omega }LP=7680$ rpm, the amplitude of $f_{\mathrm{bw}}$, which has the same frequency ratio as NL but is located on the negative half-frequency axis, rises sharply and exceeds the NL. In this speed range, disk1 undergoes reverse whirl. In Figure 7a of Section 3, it can be seen that, when $\mathrm{\Omega }LP=7680$ rpm, the amplitude of NH reaches its peak. This peak frequency is attributed to the dynamic coupling of the dual-rotor system caused by inter-shaft rub impact. Subsequently, the amplitude of NL rises again, exceeding $f_{\mathrm{bw}}$, and disk1 returns to forward whirl. To provide a more intuitive representation of the whirl direction changes of disk1, fixed-speed full-spectrum and axial trajectory analyses are conducted for these selected speed conditions, as shown in Figure 18 and Figure 19. The changes from reverse whirl motion to forward motion, then back to reverse, and finally to the forward motion of disk1 can be observed from the variations in the peak values in the full-spectrum graph.

The analysis of the above data reveals that inter-shaft rub impact primarily affects the whirl behavior of the LP-compressor disk1, with its impact on the whirl behavior of other disks being smaller compared to the effect of unbalanced excitation.

5. Inter-Shaft Rub-Impact Test

Using the existing laboratory facilities, experimental research was conducted on inter-shaft rub-impact faults under a constant rotational speed. Figure 20 shows the dual-rotor test rig with inter-shaft rub-impact faults. The highlighted section in the figure is the HP disk, corresponding to disk3 in the simulation model. An eddy-current displacement sensor is installed above the disk to capture the vertical vibration-displacement signal of the disk. The signal-acquisition system is located at the bottom left of the figure. Signals are stored in the computer via a data-acquisition device. The high-pressure rotor is driven by the motor through a belt. This test rig can only accommodate the operating condition of co-rotation of the HP rotor and LP rotor.

When the LP rotor speed is 7200 rpm, with a speed ratio $\lambda =1.5$, the HP rotor speed is 10,800 rpm. The vibration-displacement response of the HP disk when a rub-impact fault occurs is shown in Figure 21.

In Figure 21a, due to the presence of environmental noise and equipment interference during the experimental process, the fault characteristics in the time-domain plot are not very distinct. Performing an FFT transformation on the signal yields the frequency-domain plot shown in Figure 21b. In the frequency-domain plot, the main frequency components have been identified, including the rotational frequencies $f_L$ and $f_H$ of the LP rotor and HP rotor, as well as the combination frequencies and harmonics, such as $f_H-f_L$, $2f_L$, $f_H+f_L$, $2f_H$, $f_H+2f_L$, $2f_H+f_L$, and $f_H+3f_L$. The experimental results are in good agreement with the simulation results.

6. Conclusions

This study focuses on the nonlinear vibration characteristics of a dual-rotor coupled system with inter-shaft rub-impact faults and their effects on rotor-whirl behavior. First, a dynamic model of the dual-rotor system is established using the finite element method, and the mechanism of inter-shaft rub-impact faults is analyzed to develop a mathematical model. Next, the impact of rub-impact stiffness and friction coefficient on the system’s vibration characteristics is analyzed using the established model, with a focus on Campbell diagrams, 3D waterfall plots, time-frequency plots, and steady-state rub-impact force plots. Furthermore, the whirl behavior of the dual-rotor system under different inter-shaft rub-impact conditions is analyzed based on full-spectrum theory. Finally, an experimental test was conducted under constant speed conditions on the dual-rotor system with inter-shaft rub impact, validating the key conclusions drawn from the vibration characteristics analysis in Section 3. The main conclusions are as follows:

(1)

The occurrence of inter-shaft rub-impact faults is equivalent to increasing the structural stiffness of the rotor system, leading to a shift in the system’s resonance speed points. This viewpoint is also explicitly reflected in reference [21]. In the displacement response spectrum of the HP-compressor disk, multiple frequency and combination frequency components, such as $f_H-f_L$, $2f_L$, $f_H+f_L$, $2f_H$, $f_H+2f_L$, $2f_H+f_L$, and $f_H+3f_L$ appear. In Section 5 of this paper, this finding was verified through the test;

(2)

As the inter-shaft rub-impact stiffness increases, the amplitude of the high-pressure rotor’s whirl frequency $f_H$ shows a decreasing trend. This indicates that, under stable rotor operation, inter-shaft rub-impact faults constrain the vibration displacement of the high-pressure rotor. Compared to the rub-impact stiffness, the friction coefficient has a more pronounced effect on the response variation. When the dual rotors are in a state of co-rotation, with an increase in the friction coefficient, inter-shaft rub impact is no longer limited to intermittent rubbing; instead, continuous rubbing will occur within each rub cycle. This can easily lead to a self-excited vibration of the rotor. When the self-excited vibration frequency is the same as the rotational frequency, the amplitude will increase sharply, directly causing the system to become unstable, which poses a significant risk to the system;

(3)

Under unbalanced excitation, the full-spectrum graph of the dual-rotor system shows only NL and NH frequency components. In a co-rotating dual-rotor system, no backward whirl behavior occurs during acceleration, while in a counter-rotating dual-rotor system, backward whirl behavior does occur. This finding is validated by comparing it with the shaft center trajectory graph. The results obtained in this study exhibit a high degree of consistency with those in reference [16];

(4)

The inter-shaft rub-impact fault has the most significant impact on the whirl behavior of the LP-compressor disk1, while the full-spectrum response of other disks is still dominated by unbalanced excitation. During the acceleration of the LP rotor from 1000 rpm to 12,000 rpm, disk1 undergoes multiple changes in the whirl direction. The inter-shaft rub-impact fault makes the dynamic coupling of the dual-rotor system more complex, which is a significant factor in generating high-amplitude backward whirl frequencies.