Dynamic Behavior of Twin-Spool Rotor-Bearing System with Pedestal Looseness and Rub Impact

Abstract

The twin-spool rotor-bearing system plays a crucial role in the aero-engine. The potential manufacturing defect, assembly error, and abnormal working loads in the rotor-bearing system can induce multiple rotor failures, such as bolt looseness and rub impact. However, the prediction of the fault rotor dynamic behavior for the aero-engine remains a difficult frontier in numerical modeling. We present a dynamic model of the twin-spool rotor-bearing system, the failure model of bearing seat loosening, and the failure model of rub impact by using second-type Lagrangian equations, finite element theory, and the Timoshenko beam theory. In particular, to improve the accuracy of the numerical model, the rotating speed control equation and the actual aero-engine parameter are taken into account. An analysis is conducted on the impact of critical failure parameters, such as looseness stiffness and rub impact initial gap, on the vibration behaviors of the essential components of the twin-spool rotor system and on the entire engine. Additionally, this paper examines the twin-spool rotor-bearing system affected by looseness–rub coupled failures. The obtained conclusions can serve as a theoretical foundation for optimizing the structure and diagnosing faults in the aero-engine rotor system.

1. Introduction

1.1. Dynamic Modeling for Two-Spool Rotor System

Since the 1950s, axial and multiaxial rotor configurations such as twin-spool and triple-spool rotors have been increasingly used in turbojet and turbofan engines. However, there is limited literature available on the dynamic characteristics of multiaxial rotor systems. With the advancement of computer processing capacity in the early 21st century, the complicated dynamic properties of the twin-rotor construction of aero-engines have garnered significant interest among scholars [1].

On the basis of Jeffcott and Green’s work on rotor dynamics [2,3], Childs [4] constructed the dynamic equation of a twin-spool rotor model by utilizing Momentum Theory. Childs’ model incorporates “1-3-0” and “1-1-1” configurations for the high- and low-pressure rotor supports. The bearings are simplified as springs with certain stiffness and damping properties. Additionally, the model accounts for the gyroscopic effect and the impact of the engine case. While the model can capture the reciprocal impact of vibration between high-pressure (HP) and low-pressure (LP) rotors, it has undergone significant simplification at the support. Hibner [5] utilized the transfer matrix method to create a twin-spool rotor-case model. The equivalent damping force at the bearings was determined by applying the Reynolds equation. The impact of this damping force on the system response was then analyzed. However, the above model derivation process is too complicated for rapid engineering analysis.

In the 21st century, as the finite element method has become more advanced and computer processing power has improved, scholars have developed more sophisticated models of the twin-spool rotor-bearing system. This refined model has allowed for a more comprehensive analysis of the dynamic characteristics.

In modeling, Chiang and Hsu [6] employed the finite element method and the Timoshenko beam theory to develop the dynamic equations for the twin-spool rotor system. They also investigated the impact of support stiffness on the critical speed. This approach facilitated the derivation and analysis process, ensuring model accuracy. Fei [7] examined the transient dynamics properties of a twin-spool rotor system during the acceleration phase. Fu [8] conducted a study to examine the impact of speed ratio, radial clearance of inter-shaft bearing, and damping ratio on the hysteresis phenomena of a twin-spool rotor system. Hong [9] examined the impact of unevenness, rotor bending deformation, and inertia load of the wheel disk on the dynamic load of the pivot point at various rotational velocities. Wang [10] and Wang [11] developed a robust finite element model for the twin-spool rotor system. They calculated the critical rotational speeds and Campbell diagrams, and the numerical values are in good agreement with the experimental results. Based on previous studies, more aero-engine components were considered in the numerical model. Jin [12,13] established a twin-spool rotor model with squeeze-film dampers and analyzed the effect of blade loss on the vibration response. The active magnetic bearings was used to control the vibration of a co-rotating twin-spool rotor system by Singh [14].

The dynamic model of a twin-spool rotor-bearing system now takes many elements into account, and it has become more sophisticated. Future modeling research should concentrate on achieving a better balance between the accuracy of the model and degrees of freedom, creating a more effective method for handling nonlinear elements, and progressively converting qualitative analysis to quantitative analysis [15,16]. In this article, based on the previous modeling methods, the rotor-bearing system of aero-engine PW6000 is modeled with considering the transient effect and the rotational speed control strategy to improve the simulation accuracy and convergence.

1.2. Pedestal Looseness Failure

Pedestal looseness failure occurs when the pre-tensioning force on pedestal bolts is lost, resulting in the relative displacement of the connecting parts and leading to failure [17].

In the 1990s, Goldman [18] and Muszynska [19] introduced a conventional model for failure in rotor bearing pedestal foundations. This model utilizes segmental linear stiffness and damping to explain the abrupt alteration in support stiffness and damping when failure occurs. Furthermore, Goldman [18] devised pertinent validation tests for this model, and the numerical and experimental findings exhibit strong concurrence. As a result, the majority of previous investigations on failures caused by the loosening of the rotor system foundation have utilized this model. Wang [20] developed a model to study the loosening of a bearing assembly. The model is based on the traditional foundation loosening model and analyzes the time–frequency domain characteristics of the system’s steady-state response. Wang [21] utilizes the traditional loosening model to examine a twin-spool rotor system and investigates the impact of rotational speed on the dynamic behavior of the rotor system with the loosening bolt.

Currently, the classical pedestal looseness failure model effectively simulates the actual bolt loosening in bearing seats. Comprehensive investigations of the twin-spool rotor-bearing system with looseness failure should be the primary focus of future research. Analyzing failure reaction patterns in various geometries, materials, and operating settings is part of this process. Conducting thorough experimental study and designing safe and feasible experiments for loosening failure are also crucial. In this article, we focus on comprehensively analyzing the transmission of vibration and the dynamic behavior variation of each component with different looseness parameters for a twin-spool rotor-bearing system.

1.3. Coupling Failure

Imbalances, looseness, and misalignments in the twin-spool rotor system of an aero-engine can increase vibration amplitude, which in turn can induce secondary defects like rotor–stator rub impact. Yang [22] developed a model of a twin-spool rotor system with rolling bearings that takes into account the coupling faults caused by pedestal looseness and rub impact. The study examined how the stiffness of the looseness affects the time–frequency characteristics of the HP and LP rotors. It also analyzed the time–frequency response characteristics at various rotational speeds and the impact of reducing stiffness on the time–frequency domain characteristics of rotors. Additionally, the progression of friction at various rotational speeds is examined. Chen Guo [23] developed a model for a twin-spool rotor-bearing system that takes into account the effects of coupling misalignment and rub impact. The model was used to analyze the fault frequency properties of the system when coupling faults occur. The complexity of the twin-spool rotor dynamics and the interplay of rotor failure pose significant challenges in analyzing the coupling fault in the turbofan engine. Furthermore, there is scarce literature on this topic.

The purpose of this research is to study the dynamic behavior of a twin-spool rotor-bearing system with coupled faults. It specifically aims to investigate how rotor geometry and fault parameters affect the time–frequency response of HP and LP rotors. We will examine the vibration properties of a twin-spool rotor-bearing system in an aero-engine with pedestal looseness and rub impact coupled faults. On the basis of the existing twin-spool rotor-bearing model, we further take into account transient effects and rotational speed control in modeling; the effect of vibration in the acceleration phase on the dynamic behavior in the stable rotation phase is considered in calculating the fault response in order to improve the accuracy and convergence of results.

The structure of the paper is organized as follows: Firstly, the dynamic modeling of the twin-spool rotor bearing system including multiple rotating shaft, bearing, inter-shaft bearing, pedestal looseness fault and rub impact fault is set up in the Section 2. Next, the effects of pedestal looseness and rub impact on the dynamic behavior of the twin-spool rotor system are analyzed in Section 3. Finally, some conclusions are summarized in Section 4.

2. Modeling

To create a simplified dynamic model of the twin-spool rotor system that fulfills technical requirements, it is essential to use the rotor structure of the real aero-engine as a foundation. The primary components should be retained while making appropriate simplifications to the actual structure.

In this paper, a simplified dynamics model will be established with PW6000 as an example. The aero-engine twin-spool rotor system comprises HP and LP rotors with the primary mechanical connection between the two pieces being the inte-shaft bearing. Twin-spool rotor systems employ hollow rotor shafts to minimize weight. The LP rotor shaft is positioned inside the HP rotor shaft, creating a connection between the LP compressor and the LP turbine at both ends. The HP rotor shaft is linked to both the HP compressor and HP turbine at both of its extremities. The transmission of vibrations between this system mostly occurs through inter-shaft bearings, where the outer ring is in contact with the HP shaft and the inner ring is in contact with the LP shaft. The HP and LP rotor shafts are equipped with ball bearings that are directly linked to the engine case, together with the inter-shaft bearing, to guarantee the stability of the support structure. The twin-spool rotor system consists of three primary components:

• Disks: radial size is much larger than the axial size of the rotating body, such as fans, compressors, and turbines.
• Rotating shafts: axial size is much larger than the radial size of the rotating body, such as LP shaft and HP drum.
• Bearings: including inter-shaft bearings and ball bearings.

The coordinate system for the three main components will be established based on the definition of Euler angles. Using this coordinate system, the second type of Lagrange equations and the finite element method will be employed to create the disk, rotating shaft, and bearing elements. By grouping these elements, the complete differential equations of the twin-spool rotor dynamics will be obtained.

2.1. Coordinate System

According to the definition of Euler angle [24], establish the generalized coordinate system shown in Figure 1

• $OXYZ$: Fixed coordinate system with origin O.
• $O{x}^{\prime }{y}^{\prime }{z}^{\prime }$: It is obtained by rotating $OXYZ$ counterclockwise about the $O{x}^{\prime }{y}^{\prime }{z}^{\prime }$ by an angle $\Gamma$.
• $O{x}^{″}{y}^{″}{z}^{″}$: It is obtained by rotating $O{x}^{\prime }{y}^{\prime }{z}^{\prime }$ around the $y^{\prime }$ axis by an angle B counterclockwise.
• $Oxyz$: It is obtained by rotating $O{x}^{″}{y}^{″}{z}^{″}$ around the $x^{″}$ axis by an angle $\Phi$ counterclockwise, in which case the positive direction in the x axis is the rotor rotation direction.

The change of position from the fixed coordinate system $OXYZ$ to the coordinate system $Oxyz$ is expressed in the generalized form:

$$OXYZ\stackrel{OZ,\Gamma }{\to }O{x}^{\prime }{y}^{\prime }{z}^{\prime }\stackrel{O{y}^{\prime },B}{\to }O{x}^{″}{y}^{″}{z}^{″}\stackrel{O{x}^{″},\Phi }{\to }Oxyz$$

Define the displacement of a rigid body in a fixed coordinate system as $(V,W,S)$. At this point, the general motion of the rigid body can be expressed in terms of three Euler angles $(\Gamma ,B,\Phi )$ and three translational coordinates $(V,W,S)$ in the fixed coordinate system.

2.2. Disk Element

According to the principle of equivalence of dynamic properties, the components, such as the fan, compressor, and turbine, which have substantially larger radial dimensions compared to their axial dimensions, can be represented as a disk element. This subsection derives the dynamic differential equations of the disk element by employing the Lagrangian equations of the second type. It also obtains the mass matrix, gyroscopic matrix, and transient matrix (the matrix associated with rotational acceleration) of the disk element.

To streamline the analysis procedure, this paper assumes the following:

• The influence of the rotor disk deformation on the vibration characteristics of the rotor system is small and negligible.
• Due to processing or assembly errors, the mass center of the rotor does not coincide with geometry center. There is an eccentricity.
• The effect of the axial displacement on the vibration characteristics of the rotor system is small and negligible when service environments such as maneuvering flights, which change the spatial attitude of the rotor system, are not taken into consideration [25].

After introducing the assumption that the disk is considered as a rigid body and the axial displacement is neglected, the dynamic differential equations of the disk are derived using the Lagrange equations of the second type. The main steps of the derivation are as follows:

• Project the rotation velocities $\dot{\Phi },\dot{B},\dot{\Gamma }$ described by the Euler angles to the rotation coordinate system $Oxyz$.
• Write the rotational kinetic energy in a rotational coordinate system.
• Write the translational kinetic and potential energies in a fixed coordinate system.
• Write the generalized forces in a fixed coordinate system. (The main generalized forces on the disk element are gravity and eccentric forces. The geometry relation between mass center and geometry center is shown in Figure 2).
• Substituting kinetic energy, potential energy and generalized force into the Lagrange equations yields the matrix form of dynamic equations.

$${\mathbf{M}}_d{\ddot{\mathit{q}}}_d+\dot{\Phi }{\mathbf{G}}_d{\dot{\mathit{q}}}_d+\ddot{\Phi }{\mathbf{N}}_d{\mathit{q}}_d={\mathit{Q}}_d$$

(1)

where $\mathit{q}$ is the displacement array: ${\mathit{q}}_d={\{V,W,B,\Gamma \}}^{\mathrm{T}}$, ${\mathbf{M}}_d$ is the inertia matrix of the disk, ${\mathbf{G}}_d$ is the gyroscopic matrix, ${\mathbf{N}}_d$ is the transient matrix, and ${\mathit{Q}}_d$ is the generalized force on the disk:

$$\begin{array}{cc}\hfill {\mathbf{M}}_d& =\left[\begin{array}{cccc}{m}_d& 0& 0& 0\\ 0& m_d& 0& 0\\ 0& 0& I_d& 0\\ 0& 0& 0& I_d\end{array}\right],{\mathbf{G}}_d=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& I_p\\ 0& 0& -I_p& 0\end{array}\right],{\mathbf{N}}_d=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& -I_p& 0\end{array}\right]\hfill \\ \hfill {\mathit{Q}}_d& ={\left\{\begin{array}{cccc}{F}_{dy}+G_{dy}& F_{dz}& 0& 0\end{array}\right\}}^{\mathrm{T}}\hfill \end{array}$$

2.3. Shaft Element

In order to reduce the weight, the shaft in the aero-engine usually is a hollow shaft with a drum structure [26], and the whole rotor shaft is mainly subjected to transverse loads (loads in the $OYZ$). Therefore, beam models are commonly used in engineering to simulate rotor shafts in aero-engine rotor systems. Considering the bending, shearing, gyroscopic effect, and transient vibration of the shaft, the dynamic differential equations of the shaft element are obtained by the theory of the Timoshenko beams, Lagrange’s equations, and the finite element method. Any cross-section of the rotary axis in the model produces only translational displacements $(V,W)$ in the $Y,Z$ directions and rotational displacements $(\Phi ,B,\Gamma )$ around $z^{\prime },y^{″},x$, as shown in Figure 3.

The matrix form of the differential equations for the shaft element is:

$$\left(\left[{\mathbf{M}}_{Ts}\right]+\left[{\mathbf{M}}_{Rs}\right]\right){\ddot{\mathit{q}}}_s-\dot{\Phi }\left[{\mathbf{G}}_s\right]{\dot{\mathit{q}}}_s+\left[{\mathbf{K}}_s\right]{\mathit{q}}_s-\ddot{\Phi }\left[{\mathbf{N}}_s\right]{\mathit{q}}_s={\mathit{Q}}_s$$

(2)

where the mass matrix $\left[{\mathbf{M}}_{Ts}\right]$ is:

$$\left[{\mathbf{M}}_{Ts}\right]=k_{MTs}\left[\begin{array}{cccccccc}{M}_{T1}& 0& 0& M_{T4}& M_{T3}& 0& 0& -M_{T5}\\ 0& M_{T1}& -M_{T4}& 0& 0& M_{T3}& M_{T5}& 0\\ 0& -M_{T4}& M_{T2}& 0& 0& -M_{T5}& M_{T6}& 0\\ M_{T4}& 0& 0& M_{T2}& M_{T5}& 0& 0& M_{T6}\\ M_{T3}& 0& 0& M_{T5}& M_{T1}& 0& 0& -M_{T4}\\ 0& M_{T3}& -M_{T5}& 0& 0& M_{T1}& M_{T4}& 0\\ 0& M_{T5}& M_{T6}& 0& 0& M_{T4}& M_{T2}& 0\\ -M_{T5}& 0& 0& M_{T6}& -M_{T4}& 0& 0& M_{T2}\end{array}\right]$$

(3)

$k_{MTs}$ is constant:

$$k_{MTs}=\frac{{\rho }_{l}l}{{(1+{\phi }_s)}^2}=\frac{\rho Al}{{(1+{\phi }_s)}^2}$$

where ${\rho }_l$ denotes linear density, $\rho$ is density, A denotes the cross-sectional area of the rotor shaft, l denotes the length of the shaft element; ${\phi }_s=\frac{12EI}{G{A}_s{l}^2}$ is constant, $E,I,G$ represent the elastic modulus, second axial moment of area, and shear modulus, and $A_s$ is the effective shear area. And the expression for each element in the matrix is shown below:

$$\begin{array}{cc}{M}_{T1}=\frac{13}{35}+\frac{7}{10}{\phi }_s+\frac{1}{3}{\phi }_s^2\hfill & M_{T2}=\left(\frac{1}{105}+\frac{1}{60}{\phi }_s+\frac{1}{120}{\phi }_s^2\right)l^2\hfill \\ M_{T3}=\frac{9}{70}+\frac{3}{10}{\phi }_s+\frac{1}{6}{\phi }_s^2\hfill & M_{T4}=\left(\frac{11}{210}+\frac{1}{120}{\phi }_s+\frac{1}{24}{\phi }_s^2\right)l\hfill \\ M_{T5}=\left(\frac{13}{420}+\frac{3}{40}{\phi }_s+\frac{1}{24}{\phi }_s^2\right)l\hfill & M_{T6}=-\left(\frac{1}{140}+\frac{1}{60}{\phi }_s+\frac{1}{120}{\phi }_s^2\right)l^2\hfill \end{array}$$

Mass matrix $\left[{\mathbf{M}}_{Rs}\right]$:

$$\left[{\mathbf{M}}_{Rs}\right]=k_{MRs}\left[\begin{array}{cccccccc}{M}_{R1}& 0& 0& M_{R4}& -M_{R1}& 0& 0& M_{R4}\\ 0& M_{R1}& M_{R4}& 0& 0& -M_{R1}& -M_{R4}& 0\\ 0& -M_{R4}& M_{R2}& 0& 0& M_{R4}& M_{R3}& 0\\ M_{R4}& 0& 0& M_{R2}& -M_{R4}& 0& 0& M_{R3}\\ -M_{R1}& 0& 0& -M_{R4}& M_{R1}& 0& 0& -M_{R4}\\ 0& -M_{R1}& M_{R4}& 0& 0& M_{R1}& M_{R4}& 0\\ 0& -M_{R4}& M_{R3}& 0& 0& M_{R4}& M_{R2}& 0\\ M_{R4}& 0& 0& M_{R3}& -M_{R4}& 0& 0& M_{R2}\end{array}\right]$$

(4)

$k_{MRs}$ is constant:

$$k_{MRs}=\frac{I_p}{l{(1+{\phi }_s)}^2}=\frac{{\rho }_{l}I}{lA{(1+{\phi }_s)}^2}=\frac{\rho I}{l{(1+{\phi }_s)}^2}$$

where I denotes the cross-section moment of inertia, and $I_p$ denotes the polar moment of inertia. And the expression for each element in the matrix is shown below:

$$\begin{array}{cc}{M}_{R1}=\frac{6}{5}\hfill & M_{R2}=\left(\frac{2}{15}+\frac{1}{6}{\phi }_s+\frac{1}{3}{\phi }_s^2\right)l^2\hfill \\ M_{R3}=\left(-\frac{1}{30}-\frac{1}{6}{\phi }_s+\frac{1}{6}{\phi }_s^2\right)l^2\hfill & M_{R4}=\left(\frac{1}{10}-\frac{1}{2}{\phi }_s\right)l\hfill \end{array}$$

Transient matrix $\left[{\mathbf{N}}_s\right]$:

$$\left[{\mathbf{N}}_s\right]=k_{Ns}\left[\begin{array}{cccccccc}0& -N_{s1}& N_{s2}& 0& 0& N_{s1}& N_{s2}& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& -N_{s2}& N_{s4}& 0& 0& N_{s2}& -N_{s3}& 0\\ 0& N_{s1}& -N_{s2}& 0& 0& -N_{s1}& -N_{s2}& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& -N_{s2}& -N_{s3}& 0& 0& N_{s2}& N_{s4}& 0\end{array}\right]$$

(5)

$k_{Ns}$ is constant:

$$k_{Ns}=\frac{I_p}{30l{(1+{\phi }_s)}^2}=\frac{{\rho }_{l}I}{15lA{(1+{\phi }_s)}^2}=\frac{\rho I}{15l{(1+{\phi }_s)}^2}$$

The expression for each element in the matrix is shown below:

$$\begin{array}{cc}{N}_{s1}=36\hfill & N_{s2}=3l-15l{\phi }_s\hfill \\ N_{s3}=l^2+5l^2{\phi }_s-5l^2{\phi }_s\hfill & N_{s4}=4l^2+5l^2{\phi }_s+10l^2{\phi }_s\hfill \end{array}$$

Gyroscopic matrix $\left[{\mathbf{G}}_s\right]=\left(\left[{\mathbf{N}}_s\right]-{\left[{\mathbf{N}}_s\right]}^{\mathrm{T}}\right)$:

$$\left[{\mathbf{G}}_s\right]=k_{Gs}\left[\begin{array}{cccccccc}0& -G_{s1}& G_{s2}& 0& 0& G_{s1}& G_{s2}& 0\\ G_{s1}& 0& 0& G_{s2}& -G_{s1}& 0& 0& G_{s2}\\ -G_{s2}& 0& 0& -G_{s4}& G_{s2}& 0& 0& G_{s3}\\ 0& -G_{s2}& G_{s4}& 0& 0& G_{s2}& -G_{s3}& 0\\ 0& G_{s1}& -G_{s2}& 0& 0& -G_{s1}& -G_{s2}& 0\\ -G_{s1}& 0& 0& -G_{s2}& G_{s1}& 0& 0& -G_{s2}\\ -G_{s2}& 0& 0& G_{s3}& G_{s2}& 0& 0& -G_{s4}\\ 0& -G_{s2}& -G_{s3}& 0& 0& G_{s2}& G_{s4}& 0\end{array}\right]$$

(6)

$k_{Gs}=k_{Ns}$ is constant. And the expression for each element in the matrix is shown below:

$$\begin{array}{cc}{G}_{s1}=36\hfill & G_{s2}=3l-15l{\phi }_s\hfill \\ G_{s3}=l^2+5l^2{\phi }_s-5l^2{\phi }_s\hfill & G_{s4}=4l^2+5l^2{\phi }_s+10l^2{\phi }_s\hfill \end{array}$$

Stiffness matrix $\left[{\mathbf{K}}_s\right]$:

$$\left[{\mathbf{K}}_s\right]=k_{Ks}\left[\begin{array}{cccccccc}{K}_{s1}& 0& 0& K_{s4}& -K_{s1}& 0& 0& K_{s4}\\ 0& K_{s1}& -K_{s4}& 0& 0& -K_{s1}& -K_{s4}& 0\\ 0& -K_{s4}& K_{s2}& 0& 0& K_{s4}& K_{s3}& 0\\ K_{s4}& 0& 0& K_{s2}& -K_{s4}& 0& 0& K_{s3}\\ -K_{s1}& 0& 0& -K_{s4}& K_{s1}& 0& 0& -K_{s4}\\ 0& -K_{s1}& K_{s4}& 0& 0& K_{s1}& K_{s4}& 0\\ 0& -K_{s4}& K_{s3}& 0& 0& K_{s4}& K_{s2}& 0\\ K_{s4}& 0& 0& K_{s3}& -K_{s4}& 0& 0& K_{s2}\end{array}\right]$$

(7)

$k_{Ks}$ is constant:

$$k_{Ks}=\frac{EI}{l^3(1+{\phi }_s)}$$

The expression for each element in the matrix is shown below:

$$\begin{array}{cc}{K}_{s1}=12\hfill & K_{s2}=l^2({\phi }_s+4)\hfill \\ K_{s3}=l^2(2-{\phi }_s)\hfill & K_{s4}=6l\hfill \end{array}$$

In this subsection, we have developed the differential equations that describe the dynamics of the rotating shaft element using the Timoshenko beam theory. We have also found specific expressions for the mass matrix, gyroscopic matrix, stiffness matrix, and N matrix (transient matrix).

2.4. Bearing Element

In an aero-engine, the bearing is the core component that connects the rotor to the stator as well as the HP rotor to LP rotor, and it plays a supporting role. According to the different connection objects, they can be divided into ordinary bearings connected with the case and inter-shaft bearings connected with the HP and LP rotor.

The inner ring of a ordinary bearing is connected to the rotor shaft and the outer ring is fixed to the bearing housing. The bearing housing is usually not directly connected to the case. Instead, the bearing housing is attached to elastic support (such as squirrel cage elastic support), which is then connected to the case.

The inter-shaft bearing’s inner ring is linked to the LP rotor shaft, while the outer ring is attached to the HP rotor shaft. In the aero-engine rotor system, the inclusion of inter-shaft bearings can enhance the load-bearing structure and enhance the engine’s overall performance. However, this also results in the transmission of vibrations between the HP and LP rotors through the inter-shaft bearings. This leads to a more complex vibration response, which can potentially cause excessive vibrations in certain load conditions. Hence, the twin-spool rotor-bearing model in this paper takes into account the impact of both typical bearings and inter-shaft bearings on the vibration characteristics of the rotor system.

During the operation of an aero-engine, the stiffness provided by the bearing fluctuates, but overall, it remains stable and is much greater than the support stiffness provided by the elastic support [27,28]. According to this characteristic, this paper ignores the motion of the internal structure of the bearings and focuses on the influence of the overall stiffness of the bearings, and it establishes the differential equations of the ordinary bearing element and the inter-shaft bearing element through force analysis and Newton’s second law, respectively.

Ordinary bearing element. An ordinary bearing element as shown in Figure 4 is established, and the effects of bearing stiffness/damping, elastic support stiffness/damping, and pedestal mass are considered. In Figure 4a, the rotating shaft is affected by the force $Q_i$ from the bearing, and the pedestal is affected by the forces $Q_{b1},Q_{b2}$ from the bearing and the elastic support at the same time.

If it is assumed that the bearing links node i on the rotating shaft to the case, we can establish the differential equation for the dynamics of the bearing housing using Newton’s second law:

$$\begin{array}{cc}\hfill m_b{\ddot{V}}_b& =Q_{bV}=\left[k_b(V_i-V_b)+c_b({\dot{V}}_i-{\dot{V}}_b)\right]+\left[k_s(-V_b)+c_s(-{\dot{V}}_b)\right]\hfill \\ \hfill m_b{\ddot{W}}_b& =Q_{bW}=\left[k_b(W_i-W_b)+c_b({\dot{W}}_i-{\dot{W}}_b)\right]+\left[k_s(-W_b)+c_s(-{\dot{W}}_b)\right]\hfill \end{array}$$

(8)

where $Q_{bV}$ and $Q_{bW}$ represent the external forces acting on the bearing housing in the Y and Z directions, respectively. $V_i$ and $W_i$ represent the displacements of the nodes i on the rotating shaft in the Y and Z directions, respectively. Similarly, $V_b$ and $W_b$ represent the displacements of the bearing housing in the Y and Z directions, respectively. The variables $m_b, k_b, c_b$ represent the mass, stiffness, and damping of the pedestal, respectively. Similarly, $k_s$ and $c_s$ denote the stiffness and damping of the elastic support. The bearing element has an impact on both the pedestal and the i node on the shaft when the shaft experiences external forces $Q_{iV}, Q_{iW}$ from the bearing in the $Y,Z$ direction:

$$\begin{array}{cc}\hfill Q_{iV}& =-k_b(V_i-V_b)-c_b({\dot{V}}_i-{\dot{V}}_b)\hfill \\ \hfill Q_{iW}& =-k_b(W_i-W_b)-c_b({\dot{W}}_i-{\dot{W}}_b)\hfill \end{array}$$

(9)

Equations (8) and (9) are written in matrix form to obtain:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\mathit{Q}}_i\\ \mathbf{0}\end{array}\right\}=\left[\begin{array}{cc}-{\mathbf{K}}_{in}& {\mathbf{K}}_{in}\\ {\mathbf{K}}_{in}& -{\mathbf{K}}_{in2}\end{array}\right]\left\{\begin{array}{c}{V}_i\\ W_i\\ V_b\\ W_b\end{array}\right\}+\left[\begin{array}{cc}-{\mathbf{C}}_{in}& {\mathbf{C}}_{in}\\ {\mathbf{C}}_{in}& -{\mathbf{C}}_{in2}\end{array}\right]\left\{\begin{array}{c}{\dot{V}}_i\\ {\dot{W}}_i\\ {\dot{V}}_b\\ {\dot{W}}_b\end{array}\right\}-\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{M}}_{in}\end{array}\right]\left\{\begin{array}{c}{\ddot{V}}_i\\ {\ddot{W}}_i\\ {\ddot{V}}_b\\ {\ddot{W}}_b\end{array}\right\}\end{array}$$

(10)

where:

$${\mathbf{K}}_{in}=\left[\begin{array}{cc}{k}_b& 0\\ 0& k_b\end{array}\right],{\mathbf{C}}_{in}=\left[\begin{array}{cc}{c}_b& 0\\ 0& c_b\end{array}\right],{\mathbf{M}}_{in}=\left[\begin{array}{cc}{m}_b& 0\\ 0& m_b\end{array}\right],{\mathit{Q}}_i=\left\{\begin{array}{c}{Q}_{iV}\\ Q_{iW}\end{array}\right\}$$

$${\mathbf{K}}_{in2}=\left[\begin{array}{cc}{k}_b+k_s& 0\\ 0& k_b+k_s\end{array}\right],{\mathbf{C}}_{in2}=\left[\begin{array}{cc}{c}_b+c_s& 0\\ 0& c_b+c_s\end{array}\right]$$

The bearing element is incorporated into the dynamic differential equation of the complete model as an external force term within the assembly. To simplify, the equation mentioned above is enlarged into a $six$-dimensional equation. The dynamic equation of the ordinary bearing in the full model is expressed as follows:

$$\left[{\mathbf{M}}_{b1}\right]{\ddot{\mathit{q}}}_{b1}+\left[{\mathbf{C}}_{b1}\right]{\dot{\mathit{q}}}_{b1}+\left[{\mathbf{K}}_{b1}\right]{\mathit{q}}_{b1}=0$$

(11)

$$\left[{\mathbf{M}}_{b1}\right]=\left[\begin{array}{cc}{\mathbf{M}}_{b1}^{11}& {\mathbf{M}}_{b1}^{12}\\ {\mathbf{M}}_{b1}^{21}& {\mathbf{M}}_{b1}^{22}\end{array}\right]=\left[\begin{array}{ccc}\mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{M}}_{in}\end{array}\right],\left[{\mathbf{C}}_{b1}\right]=\left[\begin{array}{cc}{\mathbf{C}}_{b1}^{11}& {\mathbf{C}}_{b1}^{12}\\ {\mathbf{C}}_{b1}^{21}& {\mathbf{C}}_{b1}^{22}\end{array}\right]=\left[\begin{array}{ccc}{\mathbf{C}}_{in}& \mathbf{0}& -{\mathbf{C}}_{in}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}\\ -{\mathbf{C}}_{in}& \mathbf{0}& {\mathbf{C}}_{in2}\end{array}\right]$$

$$\left[{\mathbf{K}}_{b1}\right]=\left[\begin{array}{cc}{\mathbf{K}}_{b1}^{11}& {\mathbf{K}}_{b1}^{12}\\ {\mathbf{K}}_{b1}^{21}& {\mathbf{K}}_{b1}^{22}\end{array}\right]=\left[\begin{array}{ccc}{\mathbf{K}}_{in}& \mathbf{0}& -{\mathbf{K}}_{in}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}\\ -{\mathbf{K}}_{in}& \mathbf{0}& {\mathbf{K}}_{in2}\end{array}\right],{\mathit{q}}_{b1}={\left\{\begin{array}{cccccc}{V}_i& W_i& B_i& {\Gamma }_i& V_b& W_b\end{array}\right\}}^{\mathrm{T}}$$

where the dimensions of the submatrices $\mathbf{0},{\mathbf{M}}_{in},{\mathbf{C}}_{in},{\mathbf{K}}_{in},{\mathbf{C}}_{in2},{\mathbf{K}}_{in2}$ are $2\times 2$ and the subscript i indicates that the bearing element is connected to the node i on the rotating shaft.

Inter-Shaft Bearing. An inter-shaft bearing element is built as shown in Figure 5, connecting the i node on the HP rotor shaft with the j node on the LP rotor shaft, taking into account the effect of bearing stiffness/damping. From Figure 5a, it can be seen that the force generated by the inter-shaft bearing is applied to both the HP rotor shaft and the LP rotor shaft, and the magnitude of the force is $Q_i,Q_j$, and $Q_i=-Q_j$.

The force on the rotor shaft in the $Y,Z$ direction is:

$$\begin{array}{cc}\left\{\begin{array}{c}{Q}_{iV}=k_t(V_j-V_i)+c_t({\dot{V}}_j-{\dot{V}}_i)\hfill \\ Q_{iW}=k_t(W_j-W_i)+c_t({\dot{W}}_j-{\dot{W}}_i)\hfill \end{array}\right.& \left\{\begin{array}{c}{Q}_{jV}=-Q_{iV}\hfill \\ Q_{jW}=-Q_{iW}\hfill \end{array}\right.\end{array}$$

(12)

where $Q_{iV}, Q_{iW}, Q_{jV}, Q_{jW}$ denote the force acting on node $i,j$ in the $Y,Z$ direction, $V_i, W_i, V_j, W_j$ denote the displacement of node $i, j$ in the $Y,Z$ direction, and $k_t, c_t$ denote the bearing stiffness and damping.

Organizing the above equation into matrix form gives:

$$\left\{\begin{array}{c}{Q}_{iV}\\ Q_{iW}\\ Q_{jV}\\ Q_{jW}\end{array}\right\}=\left[\begin{array}{cc}-{\mathbf{K}}_{im}& {\mathbf{K}}_{im}\\ {\mathbf{K}}_{im}& -{\mathbf{K}}_{im}\end{array}\right]\left\{\begin{array}{c}{V}_i\\ W_i\\ V_j\\ W_j\end{array}\right\}+\left[\begin{array}{cc}-{\mathbf{C}}_{im}& {\mathbf{C}}_{im}\\ {\mathbf{C}}_{im}& -{\mathbf{C}}_{im}\end{array}\right]\left\{\begin{array}{c}\dot{V_i}\\ \dot{W_i}\\ \dot{V_j}\\ \dot{W_j}\end{array}\right\}$$

(13)

where:

$${\mathbf{K}}_{im}=\left[\begin{array}{cc}{k}_t& 0\\ 0& k_t\end{array}\right],{\mathbf{C}}_{im}=\left[\begin{array}{cc}{c}_t& 0\\ 0& c_t\end{array}\right]$$

Similar to the ordinary bearing element, the above equation is enlarged into an $eight$-dimensional equation. The dynamic equation of the inter-shaft bearing in the full model is expressed as follows:

$$\left[{\mathbf{C}}_{b2}\right]{\dot{\mathit{q}}}_{b2}+\left[{\mathbf{K}}_{b2}\right]{\mathit{q}}_{b2}=0$$

(14)

$$\left[{\mathbf{K}}_{b2}\right]=\left[\begin{array}{cc}{\mathbf{K}}_{b2}^{11}& {\mathbf{K}}_{b2}^{12}\\ {\mathbf{K}}_{b2}^{21}& {\mathbf{K}}_{b2}^{22}\end{array}\right]=\left[\begin{array}{cccc}{\mathbf{K}}_{im}& \mathbf{0}& -{\mathbf{K}}_{im}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ -{\mathbf{K}}_{im}& \mathbf{0}& {\mathbf{K}}_{im}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\end{array}\right],\left[{\mathbf{C}}_{b2}\right]=\left[\begin{array}{cc}{\mathbf{C}}_{b2}^{11}& {\mathbf{C}}_{b2}^{12}\\ {\mathbf{C}}_{b2}^{21}& {\mathbf{C}}_{b2}^{22}\end{array}\right]=\left[\begin{array}{cccc}{\mathbf{C}}_{im}& \mathbf{0}& -{\mathbf{C}}_{im}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ -{\mathbf{C}}_{im}& \mathbf{0}& {\mathbf{C}}_{im}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\end{array}\right]$$

$${\mathit{q}}_{b2}={\left\{\begin{array}{cccccccc}{V}_i& W_i& B_i& {\Gamma }_i& V_j& W_j& B_j& {\Gamma }_j\end{array}\right\}}^{\mathrm{T}}$$

where the dimensions of the submatrices $\mathbf{0}, {\mathbf{M}}_{im}, {\mathbf{C}}_{im}, {\mathbf{K}}_{im}$ are $2\times 2$; subscripts i and j indicate that the inter-shaft bearing element is connected to node $i, j$ on the rotating shaft.

2.5. Assembly

Based on the dynamic differential equations of each element and the geometric structure of the dual-rotor system, this subsection describes the method of assembling the elements by taking the twin-spool rotor system shown in Figure 6 as an example. It should be noted that the structure of the twin-spool rotor system in Figure 6 is different from that of a real aero-engine, and it is only a case study to illustrate the methodology of assembly. According to the geometry of the rotor system, nodes are set up at critical locations, such as the rotor, the bearing location, and the center of the rotor shaft. Nodes 1∼5 are located at the LP rotating shaft, nodes 6∼10 are located at the HP rotating shaft, and nodes $11, 12$ are located at the bearing housing connected to nodes $2, 10$, respectively. In addition, the inter-shaft bearing connects nodes 4 and 9 on the HP and LP shafts, and the disk element is located at node 7.

In this case, nodes 1∼10 are located in the rotating shaft, each node corresponds to 4 degrees of freedom (DOF) $V_i, W_i, B_i,$ and ${\Gamma }_i$; nodes $11 ,12$ are located in the bearing housing, each node corresponds to 2-DOF $V_i, W_i$, where the subscript i denotes the node number. The overall matrix with $44\times 44$ dimensions is shown in Figure 7. The small squares in the diagram represent the elements in the matrix, and a blank square means that the position is 0. The blocks with different colors represent different elements: blue for the shaft element, yellow for bearing 1, brown for bearing 2, hollow pink for the disk element, and red for the inter-shaft bearing. The overlapping color blocks mean that the overlapping elements are superimposed. The following describes the method of grouping each element separately.

First, the rotating shaft element is assembled. Nodes 1∼5 divide the LP shaft into four segments, and nodes 6∼10 divide the HP shaft into four segments. Each segment is a rotating shaft element. The mass, stiffness, gyroscopic, and transient matrices of each shaft element are partitioned according to the DOF of node. Take the stiffness matrix as an example:

$$\left[{\mathbf{K}}_s\right]=\left[\begin{array}{cc}{\mathbf{K}}_s^{11}& {\mathbf{K}}_s^{12}\\ {\mathbf{K}}_s^{21}& {\mathbf{K}}_s^{22}\end{array}\right]=k_{Ks}\left[\begin{array}{cccccccc}{K}_{s1}& 0& 0& K_{s4}& -K_{s1}& 0& 0& K_{s4}\\ 0& K_{s1}& -K_{s4}& 0& 0& -K_{s1}& -K_{s4}& 0\\ 0& -K_{s4}& K_{s2}& 0& 0& K_{s4}& K_{s3}& 0\\ K_{s4}& 0& 0& K_{s2}& -K_{s4}& 0& 0& K_{s3}\\ -K_{s1}& 0& 0& -K_{s4}& K_{s1}& 0& 0& -K_{s4}\\ 0& -K_{s1}& K_{s4}& 0& 0& K_{s1}& K_{s4}& 0\\ 0& -K_{s4}& K_{s3}& 0& 0& K_{s4}& K_{s2}& 0\\ K_{s4}& 0& 0& K_{s3}& -K_{s4}& 0& 0& K_{s2}\end{array}\right]$$

The similar matrices are stacked diagonally from low to high according to the node numbers, as shown in Figure 7, and so on, to complete the assembling of HP shafts.

The overall differential equations of the twin-spool rotor-bearing system are obtained:

$$\left[\mathbf{M}\right]\ddot{\mathit{q}}+\left(\left[\mathbf{C}\right]-\left[\mathbf{G}\right]\right)\dot{\mathit{q}}+\left[\mathbf{K}\right]\mathit{q}-\left[\mathbf{N}\right]\mathit{q}=\mathit{Q}$$

(15)

$$\mathit{q}={\left\{\begin{array}{ccccccccccccc}{V}_1& W_1& B_1& {\Gamma }_1& V_2& W_2& B_2& {\Gamma }_2& \cdots & V_{11}& W_{11}& V_{12}& W_{12}\end{array}\right\}}^{\mathrm{T}}$$

where $\left[\mathbf{M}\right],\left[\mathbf{C}\right],\left[\mathbf{G}\right],\left[\mathbf{K}\right],$ and $\left[\mathbf{N}\right]$ are the overall mass, damping, gyroscopic, stiffness, and transient matrices,$\mathit{Q}$ is the external force vector, and $\mathit{q}$ is the system displacement, respectively. It should be noted that the rotational speeds of the HP and LP rotors are usually not the same, and when the rotational speed ratio is $r=\frac{{\dot{\Phi }}_{\mathrm{HP}}}{{\dot{\Phi }}_{\mathrm{LP}}}$ (${\dot{\Phi }}_{\mathrm{LP}},{\dot{\Phi }}_{\mathrm{HP}}$ are the HP and LP rotor speeds, respectively), the LP part of the matrix $\left[\mathbf{G}\right]$ needs to be multiplied by the LP shaft speed $\dot{\Phi }$ and the HP part by the HP shaft speed $r\dot{\Phi }$. The LP part of the matrix $\left[\mathbf{N}\right]$ needs to be multiplied by the LP shaft acceleration $\ddot{\Phi }$ and the HP part needs to be multiplied by the HP shaft acceleration $r\ddot{\Phi }$.

2.6. Rotational Speed Control

The rotational speed variation of the twin-spool rotor system when the aero-engine PW6000 is in service is divided into three phases: ramp-up, stable operation, and ramp-down. The LP and HP rotors rotate in opposite directions, and the speed ratio is relatively stable when the PW6000 is in service.Therefore, if the LP rotor shaft rotation speed function ${\dot{\Phi }}_{\mathrm{LP}}$ and the speed ratio r are determined, the HP rotor shaft rotation speed could be calculated ${\dot{\Phi }}_{\mathrm{HP}}=r{\dot{\Phi }}_{\mathrm{LP}}$. The rotation speed function of the LP rotor shaft ${\dot{\Phi }}_{\mathrm{LP}}$ is:

$${\dot{\Phi }}_{\mathrm{LP}}=\left\{\begin{array}{cc}at\hfill & t<t_1\hfill \\ V_{max}\hfill & t_1\le t\le t_2\hfill \\ V_{max}-a(t-t_2)\hfill & t_2<t\hfill \end{array}\right.$$

(16)

where $V_{max},a$ are control parameters, denoting the maximum rotation speed and rotation acceleration, respectively. When $t<t_1$, the rotor is in the acceleration phase; when $t_1<t<t_2$, the rotor is in the stable rotation phase with speed $V_{max}$; when $t_2<t$, the rotor is in the deceleration phase.

2.7. Looseness Fault Model

If there is a manufacturing or assembly error in the bearing pedestal, the preload of the bolt may be be lost after harsh working loads, resulting in a gap between the bolt and its attached interface, which is called a looseness fault, as depicted in Figure 8. When a bolt becomes loose, the pedestal may shift within the gap caused by the loosening. Consequently, the relative position between the bearing housing and the elastic support will also be altered.

When the looseness fault occurs, the relative motion between the pedestal and the elastic support is shown in Figure 9, and the motion can be divided into three stages:

• $W_b<\alpha$: the pedestal and the elastic support are close together.
• $0\le W_b\le \alpha$: the pedestal moves within the loosening gap.
• $W_b>\alpha$: the pedestal moves away from the elastic support and is beyond the loosening gap.

According to some scholars’ studies on the looseness faults of rotor systems [22,29,30], the effects of the above looseness faults on the dynamic performance of rotor systems can be described by the piecewise linear stiffness and damping:

$$k_{\mathrm{support}}=\left\{\begin{array}{cc}{k}_s\hfill & W_b<0\hfill \\ k_l\hfill & 0<W_b<\alpha \hfill \\ k_s\hfill & W_b>\alpha \hfill \end{array}\right.$$

(17)

$$c_{\mathrm{support}}=\left\{\begin{array}{cc}{c}_s\hfill & W_b<0\hfill \\ c_l\hfill & 0<W_b<\alpha \hfill \\ c_s\hfill & W_b>\alpha \hfill \end{array}\right.$$

(18)

where $k_s, c_s$ are the stiffness and damping provided by the elastic support when no looseness failure occurs, $k_l, c_l$ are the stiffness and damping provided by the connection structure when looseness failure occurs; $W_b$ is the longitudinal displacement of the pedestal, and $\alpha$ is the length of the loosening gap.

2.8. Rub Impact Fault

The vibration amplitude in certain disks can surpass the initial gap between the rotor and stator, leading to rotor/stator rub impact failure. This is typically caused by an excessive imbalanced mass or failure of the rotor system. This section establishes a model for rotor/stator rub impact failure by analyzing the force state of the rotor.

It is assumed that the rubbing occurs between the disk and the stator, so the rubbing force acts on the disk. When the radial rubbing occurs, the rubbing force is shown in Figure 10. Coordinate $OYZ$ is the fixed coordinate system of the rotor system. $\delta$ is the initial gap between the disk and the stator; $F_n$ and $F_{\tau }$ are the radial and tangential component of the rubbing force, respectively.

The stiffness of the case in this model is much smaller than the rubbing stiffness, so it is assumed that an elastic collision occurs. The rubbing force consists of a linear contact force and a Coulomb friction force. The expression of the friction force is [31]:

$$\left\{\begin{array}{c}{F}_n\hfill \\ F_{\tau }\hfill \end{array}\right\}=\left\{\begin{array}{cc}0& e<\delta \\ k_c(e-\delta )\left\{\begin{array}{c}1\\ \mu \end{array}\right\}& e\ge \delta \end{array}\right.$$

(19)

where $e=\sqrt{V^2+W^2}$ is the radial displacement of the disk, $(V,W)$ are the concentric coordinates of the disk, $k_c$ is the radial stiffness of the case stator, and $\mu$ is the Coulomb friction coefficient between the rotor/stator. Introducing the Heaviside function $H(·)$, the friction force is obtained by decomposing the force along the $Y,Z$ coordinate axes:

$${\mathbf{F}}_r=\left\{\begin{array}{c}{F}_x\left(V,V\right)\hfill \\ F_y\left(W,W\right)\hfill \end{array}\right\}=H\left(e-\delta \right)·k_c\left(\frac{\delta }{e}-1\right)\left[\begin{array}{cc}1& -\mu \\ \mu & 1\end{array}\right]\left\{\begin{array}{c}V\hfill \\ W\hfill \end{array}\right\}$$

(20)

The physical meaning of the above equation is that when $e<\delta$, no rubbing occurs and the force is 0; when $e\ge \delta$, rubbing occurs to produce friction force.

2.9. Twin-Spool Rotor-Bearing Model of a Typical Turbofan

Based on the actual structure of the PW6000 engine (Figure 11) and applying the rotor system modeling methodology described in the previous section, the corresponding simplified dynamics model is developed in this section.

The rotor structure of the PW6000 engine consists of a fan, a four-stage LP compressor, a six-stage HP compressor, a one-stage HP turbine, and a three-stage LP turbine; three bearings are located at the LP compressor, the HP compressor, and the LP turbine; and an inter-shaft bearing, which connects the LP rotor to the HP rotor, is located at the HP turbine. According to the data published by Pratt and Whitney [32], the PW6000 has a total length of about $2.75$ m and a total mass of about 2289 kg. In this paper, it assumes that the PW6000 has a total rotor system length of about $2.4$ m and a rotor system mass of about 400 kg. Using the location of the disk and bearings as nodes, the HP and LP shafts are divided into a number of elements, as shown in Figure 12. The location information of each node is shown in

Table 1. Node information for PW6000 twin-spool rotor system.

Node	Location	Distance from the Left End (m)	Components	Node Association
1	LP Shaft	0	Fan	/
2∼5	LP Shaft	0.4670, 0.5557, 0.6443, 0.7330	LP Compressor	2 to 17
6	LP Shaft	1.9600	/	6 to 16
7∼9	LP Shaft	2.1000, 2.2500,2.4000	LP Turbine	9 to 18
10∼15	HP Shaft	0, 0.0934, 0.1868, 0.2802, 03736	HP Compressor	10 to 19
16	HP Shaft	0.9600	HP Turbine	16 to 6
17∼18	Bearing	0.4670, 2.4000	/	17 to 2, 18 to 9
19	Bearing	0	/	19 to 10

.

After node division, the material and geometry of each element need to be determined to calculate the element matrix. The fan, compressor and turbine are simplified as homogeneous discs with a certain thickness, whose geometry and materials [33] are shown in

Table 2. Material and geometry for PW6000 twin-spool rotor system.

Component	Material	Length (m)	Inner and Outer Radius (m)	Density (Kg/m3)	Elastic Modulus (GPa)	Poisson’s Ratio	Eccentricity (m)
LP Shaft	GH4169 [34]	2.4	0.015, 0.035	8240	150 (700 °C)	0.3	/
HP Shaft	GH4169	0.96	0.04, 0.056	8240	150 (700 °C)	0.3	/
Fan	TC4	0.0005	0, 0.7	4510	/	/	$8\times {10}^{-5}$
LP Compressor	VT8	0.01	0, 0.3	4430	/	/	$5\times {10}^{-5}$
HP Compressor	VT8	0.015	0, 0.2	4430	/	/	$4\times {10}^{-5}$
LP Turbine	N07718	0.01	0, 0.35	8190	/	/	$5\times {10}^{-5}$
HP Turbine	N07718	0.02	0, 0.25	8190	/	/	$5\times {10}^{-5}$

. In addition, let the bearings 1∼3 have a mass of 5 kg, a support stiffness of ${10}^8\mathrm{N}/\mathrm{m}$, and damping of $300\mathrm{N}·\mathrm{s}/\mathrm{m}$; let the inter-shaft bearing have a stiffness of ${10}^9\mathrm{N}/\mathrm{m}$ and damping of $3000\mathrm{N}·\mathrm{s}/\mathrm{m}$.

The element matrix is calculated based on the material and geometrical dimensions of each component and assembled to obtain the differential equations for the twin-spool rotor-bearing system with 70 DOFs:

$$\left[\mathbf{M}\right]\ddot{\mathit{q}}+\left(\left[\mathbf{C}\right]-\dot{\Phi }\left[\mathbf{G}\right]\right)\dot{\mathit{q}}+\left[\mathbf{K}\right]\mathit{q}-\ddot{\Phi }\left[\mathbf{N}\right]\mathit{q}=\mathit{Q}$$

(21)

where

$$\mathit{q}={\left\{\begin{array}{ccccccccccccccc}{V}_1& W_1& B_1& {\Gamma }_1& V_2& W_2& B_2& {\Gamma }_2& \cdots & V_{17}& W_{17}& V_{18}& W_{18}& V_{19}& W_{19}\end{array}\right\}}^{\mathrm{T}}$$

The subscripts in the above equation indicate the node numbers, each node located on the rotating shaft has four DOFs, and the nodes located on the bearings have two DOFs.

3. Analysis

3.1. Dynamic Analysis of the System with Pedestal Looseness

The study aims to analyze the vibration response characteristics of the system at certain rotational speeds to further investigate the influence of the looseness failure on the vibration behaviors.

Figure 13a,b show the shaft centerline orbit of the system before and after the looseness failure at a rotor speed of $1428\mathrm{rad}/\mathrm{s}$, where l denotes the distance of the node from the left end of the rotor shaft, and $V,W$ denote the displacements of the node in the $Y,Z$ directions, respectively. The shaft centerline orbit of the nodes located in the LP rotor shaft without loosening fault are regular patterns, which indicate that the nodes are in cyclic motion at this time. The shaft centerline orbit of LP compressor No. 3 is a circle, which indicates that the node is in simple harmonic motion and has only one dominant frequency, while the shaft centerline orbit of the rest of the nodes are regular and complex, which indicates that the vibration signals of the nodes contain multiple frequencies. After the looseness failure, the shaft centerline orbit of the fan and LP compressor become complex and irregular, and the shaft centerline orbit of the LP turbine is an approximate circle, which indicates that the loosening of the bearing housing located in the LP compressor mainly affects the vibration of the fan and the LP compressor, and the effect on the LP turbine is relatively small.

Figure 14a,b show the shaft centerline orbit of the HP rotor prior to and subsequent to the act of loosening. Before loosening, the HP rotor is encountering simple harmonic motion at this speed. Once the pedestal in the LP rotor shaft is loosened, the intricate vibration response is transmitted to the HP rotor via the inter-shaft bearing. This causes the vibration response of the HP rotor to become intricate as well with an overall increase in amplitude. The proximity of the node to the inter-shaft bearing directly correlates with the magnitude of the impact caused by looseness failure.

The spectrum of the steady-state responses at the fan and LP compressor before and after the looseness failure can be derived by applying a rapid Fourier transform to the vibration data, as shown in Figure 15. The fan oscillates at $227\mathrm{Hz}$, which is equal to the LP rotor frequency $f_L$. After the failure, the frequency spectrum shows $\frac{1}{2}$ times the HP rotor frequency $\frac{1}{2}{f}_H$, and the components of the frequency are close to $f_L-\frac{1}{2}{f}_H$, i.e., $147.5\mathrm{Hz}$ and $80\mathrm{Hz}$. This characteristic frequency is caused by the failure, and thus it is more obvious in the node where the failure occurs. The spectrum of the LP compressor before the breakup contains $f_L$ and $f_H$; the spectrum after the breakup shows a significant $\frac{1}{2}{f}_H$ and an approximate $f_L-\frac{1}{2}{f}_H$ frequency component.

3.2. Dynamic Analysis of the System with Rub Impact Failure

The vibration behavior of the rub impact nodes will also impact the overall vibration behavior of the entire machine. Figure 16a,b show the shaft centerline orbit of the LP rotor before and after the occurrence of a rub impact. When the gap between the rotor and stator is wide, the system rotates at a speed of $1034\mathrm{rad}/\mathrm{s}$ without experiencing any rub impact faults. The orbit of the shaft’s centerline in the entire LP rotor forms a perfect circle. During this time, the system exhibits simple harmonic motion, and the vibration of the fan, intermediary bearing, and the No.1 LP turbine remains low. Occasionally, the rotor–stator clearance is decreased during the design process of an aero-engine to enhance its performance. However, this reduction in clearance might result in rub impact faults. Figure 16b displays the centerline orbit of the entire LP rotor following the occurrence of rub impact failure in the No.3 LP turbine. The No.1 LP compressor is enclosed by a bearing, which reduces its susceptibility to failure. However, the vibration amplitude at the other nodes experiences a substantial rise and no longer exhibits simple harmonic motion. The rub impact failure primarily impacts the inter-shaft bearing and the No.1 LP turbine, causing the vibration amplitude at these two locations to increase by almost twofold and enter a condition of intricate motion. The HP rotor experiences rub impact failure like the LP rotor. This results in an almost doubling of the overall amplitude with the inter-shaft bearing being more significantly affected as it is located closer.

In order to further analyze the influence of the rub impact on the system vibration response, this paper compares the time history and spectrum at the LP turbine before and after the rub impact, as shown in Figure 17. The LP turbine is in simple harmonic motion at the LP rotor rotation frequency $f_L=164.6\mathrm{Hz}$ before the rub impact fault occurs. After the rub impact failure, there is a significant $\frac{1}{4}{f}_L$ frequency component in the spectrum. This frequency component is caused by the local rub impact; the rotor intermittently collides with the case to change its trajectory, resulting in regular fluctuations in the vibration signal.

To investigate the impact of changes in the rub impact initial gap on the vibration characteristics of the system, this study examines the variation in the shaft centerline orbit at different rub impact gap values ($\delta =2.17,2.15,2.00,1.73\mathrm{mm}$), as depicted in Figure 18. When the clearance between the components is almost equal to the amplitude of vibration, Figure 18a displays the path of the shaft’s centerline during the rub impact. At this point, the amplitude of vibration of the disk is slightly greater than the gap between the components, indicating that the disk and the case remain in contact, resulting in a complete rub impact throughout the entire rotation of the disk. As the rub impact gap decreases, the shaft centerline orbit gradually transitions to Figure 18d. The amplitude fluctuations of the system’s vibration response resemble those shown in Figure 17c, indicating that the system is experiencing local rub impact.

3.3. Dynamic Analysis of the System with Coupled Failure

This subsection examines the impact of the looseness–rub coupled failure on the vibration response of the twin-spool rotor-bearing system by comparing it with the vibration response under the looseness failure. We analyze how the vibration response of the system changes under different rotational speeds and degrees of loosening, and we derive the vibration characteristics of the looseness–rub coupled fault.

The loosening gap is $0.1\mathrm{mm}$, the loosening stiffness is $2.5\times {10}^7\mathrm{N}/\mathrm{m}$, and the loosening damping is $400\mathrm{N}·\mathrm{s}/\mathrm{m}$ to simulate a slight looseness failure; the rub impact gap is $0.8\mathrm{mm}$ and the rub impact stiffness is $2\times {10}^7\mathrm{N}/\mathrm{m}$ to simulate a rub impact fault. The effect of the coupled failure on the vibration response of the system was comparatively analyzed under these fault parameters, as shown in Figure 19. Figure 19 shows the shaft centerline orbit and spectrum at the LP compressor at a LP rotor speed of 560 rad/s. The system was operated at a speed of fL = 0.5 rpm. At this speed, the system is in periodic motion with $f_L=89.13\mathrm{Hz}$ and $f_H=115.86\mathrm{Hz}$ as the main frequencies, and the change of the shaft centerline orbit of the system before and after the occurrence of the looseness failure is not obvious. But multiplier frequencies of the HP rotor rotating frequencies, $2f_H$ and $3f_H$, appeared in the frequency spectrum after the looseness failure, as shown in Figure 19b. Considering the influence of rub impact faults on the basis of looseness failures, the shaft centerline orbit and spectrum of the system under the influence of coupling faults are shown in Figure 19c,d. The red circle in the shaft centerline orbit diagram indicates the rub impact initial gap, and it can be seen that the shaft centerline orbit of the LP compressor partially exceeds the rub impact gap and a slight local rub impact fault occurs. At the same time by the impact of local rub impact, the area of the circle surrounded by the shaft centerline orbit also increased. In addition, after the coupled failure, there is an obvious and chaotic low-frequency component in the spectrum.

The loosening gap is $0.2\mathrm{mm}$, the loosening stiffness is $1\times {10}^7\mathrm{N}/\mathrm{m}$, and the loosening damping is $400\mathrm{N}·\mathrm{s}/\mathrm{m}$ to simulate a serious looseness failure; the rub impact gap is $1.56\mathrm{mm}$, and the rub impact stiffness is $1.8\times {10}^7\mathrm{N}/\mathrm{m}$ to simulate a rub impact fault. Figure 20 shows the shaft centerline orbit and spectrum at the LP compressor at $1034\mathrm{rad}/\mathrm{s}$ rotor speed. At this speed, the system is in simple harmonic motion with $f_L=164.57\mathrm{Hz}$, and the shaft centerline orbit at the LP compressor is no longer in a regular circular shape after the looseness failure occurs; there is an obvious $3f_L$ frequency component in the spectrum, and a relatively small amount of spurious wave $f_b$ appears in the low-frequency band (Figure 20b). After considering the influence of the looseness–rubbing coupling fault, the shaft centerline orbit and spectrum are shown in Figure 20a,d. The system has already produced the local rub impact phenomenon, the multiplier component in the spectrum increases slightly, and the $f_b$ frequency component and the low-frequency component caused by loosening increase significantly.

To examine the dynamic behavior of the twin-spool rotor system affected by the looseness–rub coupled failure, the evolution law of coupled fault and their impact on the vibration behaviors of the entire machine is investigated under various levels of loosening stiffness. The variable of loosening stiffness is chosen due to its heightened sensitivity within the system’s looseness failure parameters.

A looseness failure occurs at bearing 1 (node 2) with a loosening gap of $0.5\mathrm{mm}$, a loosening damping of $400\mathrm{N}·\mathrm{s}/\mathrm{m}$, and a loosening stiffness of $10\times {10}^7$∼$0.1\times {10}^7\mathrm{N}/\mathrm{m}$ to simulate the evolution from no fault to a severe looseness failure, a rub impact gap of $1.6\mathrm{mm}$ and a rub impact stiffness of $6\times {10}^6\mathrm{N}/\mathrm{m}$ to simulate the rub impact fault of the No. 4 LP compressor. The spectrum of each node with different loosening stiffness is shown in Figure 21 for the system affected by looseness failures and looseness–rub coupled failure and operated at a rotational speed of 1034 rad/s with a rotational speed ratio of $-1.3$.

Figure 21a,b show the frequency spectrum of the fan during the effects of looseness failure and looseness–rub coupled faults, respectively. It can be seen from the figure that the main frequency of the fan vibration is $f_L$ in the case of slight looseness failure, and as the loosening stiffness decreases, the frequency spectrum around $80\mathrm{Hz}$ has a high percentage of spurious waves, and there is a frequency component close to $\frac{3}{2}{f}_L$; when the loosening is coupled with the rub impact, the rub impact leads to an increase in the low-frequency component around $80\mathrm{Hz}$, and enhances the frequency component of $\frac{3}{2}{f}_L$. The subharmonic vibration at the fan following a significant looseness failure is particularly noteworthy. In this case, the vibration is primarily characterized by the $f_b$ component.

The frequency spectrum of response at the LP compressor is similar to that at the fan, except that there is no subharmonic vibration at the LP compressor, and the $f_b$ frequency component is relatively small. Figure 21c,d show the spectrum of the rotational displacement response of the LP compressor, respectively. Compared with the translational displacement, the vibration response of the rotational displacement is more susceptible to coupled faults, and as the degree of faults becomes higher, the $f_b$ and $\frac{3}{2}{f}_L$ components in the frequency spectrum increase rapidly and eventually exceed the $f_L$ frequency component. The vibration response at the LP turbine is similar to that of the LP compressor.

The frequency spectrum of the response at the HP compressor is shown in Figure 21e,f. The spectrum of the HP compressor contains both LP and HP rotor frequencies, $f_L$ and $f_H$, when there is no fault, and there are obvious $2f_L$ and $3f_L$ multiplier components in the spectrum just after the system has been loosened. As the fault level increases, the $f_b$ frequency component begins to appear in the spectrum, and the rub impact enhances this component and the multiplier components. The HP turbine is located in the vicinity of the inter-shaft bearing and is more affected by the vibration of the LP rotor, so the changes in the frequency spectrum are similar to those of the LP rotor.

4. Conclusions

In this paper, the dynamic behavior of the twin spool rotor-bearing system with looseness–rub coupled failure is investigated based on the proposed simulation model, and the conclusions are as follows:

The impact of looseness failure on the system varies depending on the speed. At low speeds, the system is susceptible to multiplier components, but the centerline orbit remains mostly unchanged. At high speeds, the system tends to generate a complex motion spectrum with noticeable low-frequency components. As the degree of looseness increases, the system may transition into chaotic motion once the loosening stiffness diminishes to a specific threshold value. If the case stiffness is low or the amplitude is near the touching gap, the system may experience complete rub. As the rub impact stiffness increases, the system will transition from complete rub to local rub, resulting in the appearance of frequency components below the main frequency in the spectrum.

Based on the results of the single failure vibration characteristic analysis, the influence of the looseness–rub coupled failure on the twin-spool rotor-bearing system vibration characteristic is further considered. The frequency components that remain constant regardless of rotational speed are observed in the frequency spectrum around 50∼80 Hz during the occurrence of the looseness–rub coupled failure. The rub impact noticeably amplifies the presence of 50∼80 Hz in the spectrum.

Further research in this field should extend to (1) the comprehensive experiments for the failure twin-spool rotor-bearing system and (2) the more general dynamic simulation model with a case, squeezed-film damper, suspension component, and bolted disk-drum structure.