# Coupled Lateral and Torsional Vibration of Rub-Impact Rotor during Hovering Flight

Authors: Nan Zheng; Moli Chen; Guihuo Luo; Zhifeng Ye
Journal: Shock and Vibration
 (2021)

---

## Abstract

When aircraft make a maneuvering during flight, additional loads acting on the engine rotor system are generated, which may induce rub-impact faults between the rotor and stator. To study the rub-impact response characteristics of the rotor system during hovering flight, the dynamic model of a rub-impact rotor system is established with lateral-torsional vibration coupling effect under arbitrary maneuvering flight conditions using the finite element method and Lagrange equation. An implicit numerical integral method combining the Newmark-β and Newton–Raphson methods is used to solve the vibration response. The results indicate that the dynamic characteristics of the rotor system will change during maneuvering flight, and the subharmonic vibrations are amplified in both lateral and torsional vibrations due to maneuvering overload. The form of the rub-impact is different during level and hovering flight conditions: the rub-impact may occur at an arbitrary phase of the whole cycle under the condition of level flight, while only local rub-impact occurs during hovering flight. Under the both flight conditions, the rub-impact has a large effect on the spectral characteristics, periodicity, and stability of the rotor system.


---

## Body

## 1. Introduction

Maneuverability is an important index for evaluating the performance of military aircraft. Especially for future military drone, the tolerant limit of pilots need not to be considered, and hence, the maneuvering overload will be increased substantially in the pursuit of higher maneuverability. However, additional loads acting on the rotor system of aeroengine will be generated during maneuvering flight, which change the vibration characteristics of the rotor [1, 2].In recent years, the vibration characteristics of rotor systems during maneuvering flight have been extensively investigated. Sakata et al. [3, 4] studied the steady-state response of a rotor system under the maneuvering flight actions of pitching and hovering, considering the gyroscopic moment generated by the rotation of foundation. Sherf et al. [5] investigated the vibration response of an F15 fighter engine system in the conditions of bench test, level flight, and maneuvering flight, respectively; they found that compared with the bench test, some parts of the engine undergo larger vibrations during level flight, and these vibrations are amplified again during maneuvering flight. Xu and Liao [6] studied the vibration response characteristics of a Jeffcott rotor with squeeze film dampers (SFDs) in maneuvering flight environment, and the numerical analysis indicated that the SFDs effectively suppress the transient vibrations caused by maneuvering flight. Based on this, Wang et al. [7] designed a rotor test rig that applies time-varying centrifugal force and gyroscopic moment on the rotor system according to maneuvering flight. The influence of the oil film clearance, the unbalanced force, and maneuvering flight parameters on the vibration reduction performance of SFDs was then analyzed. Lin and Meng [8, 9] studied the dynamic characteristics of a Jeffcott rotor with initial bending, stiffness asymmetry, and rub-impact faults in maneuvering flight environment. Ananthan and Leishman [10] investigated the vibration response of the rotor system of a turbo shaft engine under the maneuvering actions of diving, jumping, rolling, reverse rolling, and pulling up during descent. The results showed that an unstable additional load was generated and the vibration amplitude increased during these maneuvers. Besides, Abhishek et al. [11] studied the aerodynamics and structural dynamics of the engine rotor system of a UH-60A helicopter during maneuvering flight and investigated the mechanism of the unsteady rotating stall phenomenon observed during the pulling up maneuver. Zhu and Chen [12, 13] established a general dynamic model of a rotor system under arbitrary maneuvering flight conditions using the Lagrange equation and numerically analyzed the influence of the maneuvering load on the vibration characteristics of the rotor. Yang et al. [14] studied the nonlinear dynamical response of a cracked rotor system with rigid supports during maneuvering flight and observed three routes whereby the motion of the rotor would become chaotic. Using the finite element method and the Lagrange equation, Das et al. [15, 16] derived the differential equations of motion for a flexible rotor during maneuvering flight and proposed an active control strategy using electromagnetic actuators to suppress the vibration amplitude caused by maneuvering flight. Hou and Chen [17, 18] established the motion equations for a rotor system during maneuvering flight, and the rotor system with nonlinear elastic supports and the system with rolling were studied using numerical and analytic methods; the subharmonic vibration, irregular motion, and bifurcation of the rotor system were also investigated. Chen et al. [19] established the dynamic equations of a flexible rotor system with bearing under a time-varying base motion using the finite element method and the Lagrange equation. The dynamic characteristics of the rotor system were studied using numerical and experimental methods, and some nonlinear vibration behavior due to the rotation of the base was observed. Based on this, the vibration characteristics of an asymmetric rotor system with time-varying base motion were investigated [20]. Han and Ding [21] analyzed the transient response characteristics of a rotor system with SFDs under maneuvering flight using the finite element method. Gao et al. [22] established the motion equations of an asymmetric rotor system with intermediate bearing and SFDs during maneuvering flight and studied the dynamic response characteristics of the rotor system under the conditions of acceleration, rotation, and pitching of the aircraft by numerical method. Chen et al. [23] presented the differential equations of motion for a rotor system with base motion considering the nonlinear forces of the bearings and SFDs. The transient response, steady-state response, and periodic stability of the whirl of the rotor under base excitation were analyzed. Considering the initial static eccentricity of the SFDs, Chen et al. [24] studied the transient dynamic response of a rotor system under maneuvering flight conditions.The abovementioned studies have found that maneuvering has a strong influence on the dynamic characteristics of the rotor system. First, the nonlinear characteristics of the bearings and SFDs in the rotor system will change during maneuvering flight. Second, the additional centrifugal force and gyroscopic moment due to a particular maneuver will change the center and trajectory of the rotor whirl, which in turn greatly increases the probability of rub-impact faults between the rotor and stator. Hou et al. [25, 26] studied the vibration characteristics of a rub-impact rotor with nonlinear support during maneuvering flight using the harmonic balance method combined with time-domain method and examined the periodicity and stability of the rotor system through bifurcation diagrams and the Lyapunov exponent spectrum. Considering the effects of the swing of the rigid disk and stagger angles of blades, Ma et al. [27, 28] established a rotor-blade dynamic model using the finite element method and analyzed the rubbing vibration response of the rotor system. Li et al. [29] studied the rub-impact vibration response and stability of a rotor-blade system with nonlinear support at both ends by multiscale method, and the effects of normal rub-impact force, friction coefficient, damping, and bearing stiffness on the nonlinear vibration response were investigated.The lateral vibration and torsional vibration coexist in real rotor systems, and there is a coupling relationship between them due to the unbalanced mass of the rotor, resulting in complex nonlinear dynamical characteristics. Cohen and Porat [30] analyzed the stability of a rotor system considering the lateral and torsional coupling vibration due to the unbalanced force. Mohiuddin and Khulief [31] and Al-Bedoor [32] studied the lateral-torsional coupling vibration characteristics of a rotor system using numerical methods, and some nonlinear phenomena were observed, such as the softening effect and energetic interaction between the lateral and torsional vibrations. Darpe et al. [33] modified the Timoshenko beam element theory in establishing the differential equations of motion for a cracked rotor considering lateral, axial, and torsional vibrations; the coupled vibration response of the cracked rotor was analyzed. Yuan et al. [34] presented a dynamic model of a rotor system with six degrees of freedom at each node and proposed two coupling modes of internal coupling and external coupling. The stability of the system was analyzed by the harmonic balance method and Floquet theory, and the lateral-torsional coupling effect was investigated through numerical simulations. Roques et al. [35] modeled the dynamics model of a rotor system with rubbing fault due to blade loss, the possibility of lateral-torsional coupling vibration induced by rub-impact was discussed, and the load acting on the bearing was analyzed. Patel and Darpe [36] studied the dynamic characteristics of a rotor system with fatigue crack and rub-impact faults, and they investigated the influence of local friction excitation and nonlinear stiffness variation due to crack breathing on the lateral and torsional vibration responses of the rotor. Cao et al. [37] established a bearing-rotor model that considers the nonlinear forces of the bearings and SFDs and examined the nonlinear transient response of the combined lateral and torsional vibrations. Hong et al. [38] derived the differential equations for a rotor system considering the coupling effect of lateral and torsional vibrations due to rotor unbalance and analyzed the modal characteristics of the rotor using Floquet theory and Hill method.In the studies described above, it can be noted that there have been few studies on the dynamic characteristics of rub-impact rotors during maneuvering flight. Besides, only lateral vibration has been analyzed in the research on rotor dynamic characteristics during maneuvering flight, with the torsional vibration being completely neglected. Actually, during maneuvering flight of aircraft, the lateral-torsional coupling may be amplified, which in turn changes the vibration characteristics of the rotor system. Moreover, the lateral-torsional coupling vibration characteristics can provide useful information for fault diagnosis of the rotor system and improve the accuracy of diagnosis. In this study, the dynamic motion equations are derived using the finite element method and the Lagrange equation for a rub-impact rotor system considering lateral-torsional coupling under arbitrary maneuver flight conditions. Numerical integration is employed to solve the dynamic response of the rotor system during hovering flight, and the spectrum characteristics, periodicity, and stability of the rotor motion are analyzed through spectrum diagrams, bifurcation diagrams, and Poincaré maps.


## 2. Dynamic Modeling of Rotor System

### 2.1. Finite Element Formulation

The lateral-torsional coupled vibration equations are now derived using the Lagrange equation for the disk element and beam element under arbitrary maneuvering flight conditions.


#### 2.1.1. Disk Element

The coordinate system of an arbitrary disk of the rotor system is shown in Figure1, where OXYZ is the ground coordinate system for describing the spatial position, velocity, and acceleration of the aircraft, and oxyz is the rotor coordinate system for describing the translation and rotation of the rotor relative to the aircraft. Let us suppose that the engine rotor is installed on the aircraft such that one end coincides with the center of gravity of the aircraft, and thus, oxyz is also the relative coordinate system of the aircraft. The generalized displacement vector of the five degrees-of-freedom (5-DOFs) disk is given as(1)qd=xyθxθyθzT.Figure 1
Coordinate system of disk.Then, the kinetic energy of the disk can be expressed as [13, 15] follows:(2)T=12mr˙Tr˙+vBTvB−rTω⌣B2r+2r˙Tω⌣Br+2vBTr˙+2vBTω⌣Br+12Idθ˙x+ωB,X2+12Idθ˙y+ωB,Y2+12IpΩ2+θ˙z2+2θ˙zΩ+2θy+θB,Yθ˙x+ωB,Xθ˙z+Ω,(3)r=x+ecosΩt+θz+φy+esinΩt+θz+φz,vB=X˙BY˙BZ˙B,ω⌣B=0−ωB,ZωB,YωB,Z0−ωB,X−ωB,YωB,X0,where m denotes the quality of the disk, Id and Ip are the diameter and polar moment of inertia, respectively, Ω is the rotational speed of the rotor, r denotes the coordinate position of the disk in the oxyz coordinate system, e is the eccentricity distance, φ represents the initial phase angle, vB and ωB are the velocity and angular velocity of the aircraft in the OXYZ coordinate system, X˙B, Y˙B, Z˙B, and ωB,X, ωB,Y, and ωB,Z denote the velocity and angular velocity of the aircraft in the OXYZ coordinate system, respectively, and θB,Y is the angular displacement of the aircraft around the OY axis, and(4)θy+θB,Y≈θy,θy+θB,Y′≈θ˙y+ωB,Y.Denoting the damping matrix and stiffness matrix of the left and right beam elements adjacent to the disk asCl, Cr, Kl, and Kr, respectively, whose dimensions are all 10 × 10. Then, the damping matrix Cd and stiffness matrix Kd of the disk can be expressed as(5)Cd=C6,6l+C1,1rC6,7l+C1,2rC6,8l+C1,3rC6,9l+C1,4rC6,10l+C1,5rC7,6l+C2,1rC7,7l+C2,2rC7,8l+C2,3rC7,9l+C2,4rC7,10l+C2,5rC8,6l+C3,1rC8,7l+C3,2rC8,8l+C3,3rC8,9l+C3,4rC8,10l+C3,5rC9,6l+C4,1rC9,7l+C4,2rC9,8l+C4,3rC9,9l+C4,4rC9,10l+C4,5rC10,6l+C5,1rC10,7l+C5,2rC10,8l+C5,3rC10,9l+C5,4rC10,10l+C5,5r,(6)Kd=K6,6l+K1,1rK6,7l+K1,2rK6,8l+K1,3rK6,9l+K1,4rK6,10l+K1,5rK7,6l+K2,1rK7,7l+K2,2rK7,8l+K2,3rK7,9l+K2,4rK7,10l+K2,5rK8,6l+K3,1rK8,7l+K3,2rK8,8l+K3,3rK8,9l+K3,4rK8,10l+K3,5rK9,6l+K4,1rK9,7l+K4,2rK9,8l+K4,3rK9,9l+K4,4rK9,10l+K4,5rK10,6l+K5,1rK10,7l+K5,2rK10,8l+K5,3rK10,9l+K5,4rK10,10l+K5,5r.The dissipation energy and potential energy can be expressed as(7)D=12q˙dTCdq˙d,V=12qdTKdqd.Generally, the angular displacement amplitude of torsional vibration is very small. Under this assumption, the following relations are obtained:(8)cosΩt+θz+φ=cosΩt+φcosθz−sinΩt+φsinθz≈cosΩt+φ−θzsinΩt+φ,sinΩt+θz+φ=sinΩt+φcosθz+cosΩt+φsinθz≈sinΩt+φ+θzcosΩt+φ.Substituting equations (1)–(8) into the Lagrange equation, the motion differential equations of the 5-DOFs disk under arbitrary flight conditions can be obtained as follows:(9a)Md+Mduq¨d+Cd+ΩGd+Cdu+Cdbq˙d+Kd+Kdu+Kdbqd=Fdu+Fdb+Qd,(9b)Md=m00000m00000Id00000Id00000Ip,(9c)Gd=0000000000000Ip000−Ip0000000,(9d)Mdu=0000−mesinΩt+φ0000mecosΩt+φ0000000000−mesinΩt+φmecosΩt+φ00me2,(9e)Cdu=−2Ω0000mecosΩt+φ0000mesinΩt+φ000000000000000,(9f)Kdu=0000meΩ2sinΩt+φ0000meΩ2cosΩt+φ000000000000000,(9g)Fdu=meΩ2cosΩt+φmeΩ2sinΩt+φ000,(9h)Cdb=0−2mωB,Z0002mωB,Z00000000IpωB,Y0000−IpωB,X00IpωB,YIpωB,X0,(9i)Kdb=−mωB,Y2+ωB,Z2mωB,XωB,Y−ω˙B,Z000mωB,XωB,Y+ω˙B,Z−mωB,X2+ωB,Z20000000000000000Ipω˙B,Y0,(9j)Fdb=−mX¨B−ωB,ZY˙B+ωB,YZ˙B−mzωB,XωB,Z+ω˙B,Y−mY¨B−ωB,ZX˙B−ωB,XZ˙B−mzωB,YωB,Z−ω˙B,X−IpΩωB,Y−Idω˙B,XIpΩωB,X−Idω˙B,Y−IpωB,XωB,Y−IpθB,Yω˙B,X,where Md and Gd are the traditional mass and gyro matrices without lateral-torsional coupling or maneuvering flight conditions. Fdu is the unbalanced force, and Qd is the generalized external force acting on the disk. Mdu, Cdu, and Kdu denote the coupling items due to the lateral-torsional coupling effect; Cdb, Kdb, and Fdb represent the additional damping, stiffness, and excitation force resulting from maneuvering flight, respectively.In general, the additional stiffness induced by maneuvering flight is very small compared with the inherent stiffness of the rotor system, and the additional damping is related to the rotational angular velocity of the aircraft around each coordinate axis in theOXYZ coordinate system. If the rolling angular velocity of the aircraft is large enough when the aircraft performs the maneuver of rolling, the additional damping may reach the level of the inherent damping of the rotor system.Denoting the additional excitation force asFdb1 when the aircraft performs a hovering flight in the horizontal plane, Fdb2 when the aircraft performs a diving and pulling action in the vertical plane and then(10)Fdb1=−mωB,YZ˙B0−IpΩωB,Y00,(11)Fdb2=0mωB,XZ˙B0IpΩωB,X0.It can be seen that under the two special maneuvering flight conditions, the first two items of the additional excitation forces are the centrifugal force caused by the maneuvering flight action, and the 3rd and 4th items are the maneuvering gyro torque.


#### 2.1.2. Shaft

Figure2 shows a schematic diagram of a Timoshenko beam element with both lateral and torsional DOFs. A cylindrical segment with infinite small thickness is taken for analysis at an arbitrary axial position of the beam element. Suppose the thickness of the cylinder is dξ, and the distance between the cylinder and one end of the beam element is ξ, then the kinetic energy of the cylindrical segment can be written as(12)dTe=12ρeAer˙eTr˙e+vBTvB−reTω⌣B2re+2r˙eTω⌣Bre+2vBTr˙e+2vBTω⌣Bredξ+12dIpΩ2+θ˙ze2+2θ˙zeΩ+2θ˙xe+ωB,Xθye+θB,Yθ˙ze+Ω+12dIdθ˙xe+ωB,X2+12dIdθ˙ye+ωB,Y2,(13)re=xeyeze,(14)xeyeθxeθyeθze=N100N20N300N400N1−N2000N3−N4000−N⌣1N⌣2000−N⌣3N⌣400N⌣100N⌣20N⌣300N⌣400000N50000N⌣5qe,(15)qe=x1y1θx1θy1θz1x2y2θx2θy2θz2T,where qe is the nodal displacement vector, Nii=1−5 and N⌣ii=1−5 are the shape functions of the beam element [39], and re denotes the coordinate position of the cylindrical segment in the oxyz coordinate system. The kinetic energy of the beam element can be obtained by integrating dTe along the oz axis. Similar to the derivation process for the disk, the dynamic differential equations of the beam element under arbitrary maneuvering flight conditions can be obtained as(16a)Meq¨e+Ce+ΩGe+Cebq˙e+Ke+Kebqe=Feb+Qe,(16b)Ceb=2ρeAeωB,ZH1+2ρeIeωB,XH2+2ρeIeωB,YH3,(16c)Keb=ρeAeωB,Y2+ωB,Z2H4+ρeAeωB,X2+ωB,Z2H5+ρeAeωB,XωB,Y−ω˙B,ZH6+ρeAeωB,XωB,Y+ω˙B,ZH7+2ρeIeω˙B,XH8,(16d)Feb=ρeAeY˙BωB,Z−Z˙BωB,Y−ω˙B,Yze−X¨B−ωB,XωB,ZzeA1−X˙BωB,Z+Z˙BωB,X+ω˙B,Xze−Y¨B−ωB,YωB,ZzeA1X˙BωB,Z−Z˙BωB,X−ω˙B,Xze+Y¨B+ωB,YωB,ZzeA2Y˙BωB,Z−Z˙BωB,Y−ω˙B,Yze−X¨B−ωB,XωB,ZzeA20Y˙BωB,Z−Z˙BωB,Y−ω˙B,Yze−X¨B−ωB,XωB,ZzeA3−X˙BωB,Z+Z˙BωB,X+ω˙B,Xze−Y¨B−ωB,YωB,ZzeA3X˙BωB,Z−Z˙BωB,X−ω˙B,Xze+Y¨B+ωB,YωB,ZzeA4Y˙BωB,Z−Z˙BωB,Y−ω˙B,Yze−X¨B−ωB,XωB,ZzeA40+ρeIe2ΩωB,X−ω˙B,YB12ΩωB,Y+ω˙B,XB1−2ΩωB,Y+ω˙B,XB22ΩωB,X−ω˙B,YB2−2θB,Yω˙B,X+ωB,XωB,YA52ΩωB,X−ω˙B,YB32ΩωB,Y+ω˙B,XB3−2ΩωB,Y+ω˙B,XB42ΩωB,X−ω˙B,YB42θB,Yω˙B,X+ωB,XωB,YB5,(16e)Ai=∫0leNidξ,Bi=∫0leN⌣idξ,i=1∼5,where Me, Ge, and Ke are the traditional mass, gyro, and stiffness matrices, respectively (their expressions are omitted due to space limitations and refer Friswell et al. [39]). Ce denotes the external linear damping, and Qe represents the generalized external force. Ceb, Keb, and Feb are the additional damping, stiffness, and excitation force due to maneuvering flight, respectively. The matrices H1–H8 are related to the shape functions and geometric dimension of the beam (see the appendix).Figure 2
Diagram of Timoshenko beam element.Based on finite element theory, the dynamic differential equations of the rotor system under arbitrary maneuvering flight conditions can be written as follows by rearranging the matrices of the disk elements and beam elements [40].(17)M+Muq¨+C+ΩG+Cu+Cbq˙+K+Ku+Kbq=Fu+Fb+Q,where M, C, G, and K are the inertia, damping, gyro moment, and stiffness matrices of the rotor system without lateral-torsional coupling effect or maneuvering flight conditions. Mu, Cu, and Ku denote the additional inertia, damping, and stiffness due to lateral-torsional coupling effect, and Cb, Kb, and Fb represent the additional damping, stiffness, and excitation force vector caused by maneuvering flight. Fu is the unbalanced force, and Q represents the generalized external force acting on the rotor, including the rub-impact force and the nonlinear forces of the supports.


### 2.2. Nonlinear Forces of the Supports

Based on Hertz contact theory, the nonlinear forces of a rolling ball bearing can be written as [41](18a)fxbearing=−kn∑j=1NbuθjξHuθjcosθj,fybearing=−kn∑j=1NbuθjξHuθjsinθj,(18b)θj=2πj−1Nb+ωbmt,(18c)uθj=xir−xorcosθj+yir−yorsinθj−γ2,(18d)Huθj=0,uθj≤0,uθj,uθj>0,(18e)ωbm=ωinr1+ωoutr2r1+r2,where the superscripts “ir” and “or” denote the inner and outer rings of the bearing, respectively. kn is the contact stiffness between the roller and the ring, Nb is the number of rollers, and θj is the circumferential position of the jth roller at time t. uθj represents the elastic contact deformation of the jth roller in the radial direction. γ is the radial clearance. r1 and r2 denote the radius of the inner and outer rings, respectively, ωin and ωout are the rotational speeds of the inner and outer rings, and ωbm is the rotational speed of the revolution of the rollers. Generally, ξ=3/2 for rolling ball bearings.As for SFD, during the process of whirling of the rotor, the nonlinear oil film forces of the SFD can be written as [42, 43](19a)fxsfd=−μRL3c2x2+y21/2xε˙I2+εψ˙I1+yε˙I1+εψ˙I3,fysfd=−μRL3c2x2+y21/2yε˙I2+εψ˙I1−xε˙I1+εψ˙I3,(19b)ε=x2+y2c,(19c)ε˙=xx˙+yy˙cx2+y2,(19d)ψ˙=yx˙−xy˙x2+y2,(19e)tanψ=xy,where x and y are the eccentric positions of the journal, L, R, and c denote the length, radius, and radial clearance of the SFD, respectively, μ represents the dynamic viscosity, and the Sommerfeld integral Ijj=1−3 can be expressed as(20)I1=∫θ1θ2cos2θ1+εcosθ3dθ,I2=∫θ1θ2sinθcosθ1+εcosθ3dθ,I3=∫θ1θ2sin2θ1+εcosθ3dθ,where θ1=tan−1−ε˙/εψ˙, θ2=θ1+π.


### 2.3. Rub-Impact Forces

Considering the deformation and energy consumption in the process of rubbing, the normal rub-impact forceFN can be expressed as follows based on Hertz contact theory [44]:(21a)FN=kδq+cδ˙,(21b)c=1.5αkδq,(21c)δ=xi2+yi2+rin−rout,(21d)δ˙=x˙i2+y˙i2,(21e)k=1.5681−v12E1+1−v22E2−4/5rin2δ2/5m2δ˙21/5,where c is the viscous damping coefficient in the normal direction, α denotes the energy dissipation factor, δ is the invasion displacement between the rotor and stator (when δ>0, rub-impact occurs), δ˙ is the relative velocity of the rotor and stator in the vertical direction of the contact surface, xi and yi denote the displacement of node i of the rotor, and x˙i, y˙i the velocity. rin and rout represent the outer diameter of the rotor and the inner diameter of the casing, respectively, and k is the rubbing stiffness coefficient. v1, v2, E1, and E2 are the Poisson’s ratio and elastic moduli of the rotor and stator, respectively, and m2 is the equivalent quality of the rotor when rub-impact occurs.Based on Coulomb friction theory, the tangential friction due to rub-impact can be expressed as [45](22a)Fτ=ΘμfFN,(22b)Θ=−1,ψ<0,0,ψ=0,1,ψ>0,(22c)ψ=uv˙−vu˙v2+u2+rinω,where μf is the friction coefficient, Θ is employed to judge the direction of the tangential friction, and u and v denote the displacements of the rub-impact node in x and y directions, respectively.Figure3 shows the schematic diagram of rub-impact forces between the rotor and stator considering the deformation of the stator under maneuvering flight, where o is the initial centroid of the rotor and stator, and o′(x′, y′) the centroid of the stator after deformation due to maneuvering overload, and o1(u, v) the centroid of the rotor when rub-impact occurs. The rotor rotates counterclockwise. Then, the rub-impact forces in the rotor coordinate system can be written as(23)Fx=−FNcosθ+Fτsinθ,Fy=−FNsinθ−Fτcosθ,Mz=Fτrin,where Mz is the torque acting on the rotor due to rub-impact, and θ=arctanv−y′/u−x′. The translational displacement of the stator x′ and y′ can be determined according to the maneuvering overload and the deformation stiffness of the stator.Figure 3
Schematic diagram of rub-impact forces.


## 2.1. Finite Element Formulation

The lateral-torsional coupled vibration equations are now derived using the Lagrange equation for the disk element and beam element under arbitrary maneuvering flight conditions.


### 2.1.1. Disk Element

The coordinate system of an arbitrary disk of the rotor system is shown in Figure1, where OXYZ is the ground coordinate system for describing the spatial position, velocity, and acceleration of the aircraft, and oxyz is the rotor coordinate system for describing the translation and rotation of the rotor relative to the aircraft. Let us suppose that the engine rotor is installed on the aircraft such that one end coincides with the center of gravity of the aircraft, and thus, oxyz is also the relative coordinate system of the aircraft. The generalized displacement vector of the five degrees-of-freedom (5-DOFs) disk is given as(1)qd=xyθxθyθzT.Figure 1
Coordinate system of disk.Then, the kinetic energy of the disk can be expressed as [13, 15] follows:(2)T=12mr˙Tr˙+vBTvB−rTω⌣B2r+2r˙Tω⌣Br+2vBTr˙+2vBTω⌣Br+12Idθ˙x+ωB,X2+12Idθ˙y+ωB,Y2+12IpΩ2+θ˙z2+2θ˙zΩ+2θy+θB,Yθ˙x+ωB,Xθ˙z+Ω,(3)r=x+ecosΩt+θz+φy+esinΩt+θz+φz,vB=X˙BY˙BZ˙B,ω⌣B=0−ωB,ZωB,YωB,Z0−ωB,X−ωB,YωB,X0,where m denotes the quality of the disk, Id and Ip are the diameter and polar moment of inertia, respectively, Ω is the rotational speed of the rotor, r denotes the coordinate position of the disk in the oxyz coordinate system, e is the eccentricity distance, φ represents the initial phase angle, vB and ωB are the velocity and angular velocity of the aircraft in the OXYZ coordinate system, X˙B, Y˙B, Z˙B, and ωB,X, ωB,Y, and ωB,Z denote the velocity and angular velocity of the aircraft in the OXYZ coordinate system, respectively, and θB,Y is the angular displacement of the aircraft around the OY axis, and(4)θy+θB,Y≈θy,θy+θB,Y′≈θ˙y+ωB,Y.Denoting the damping matrix and stiffness matrix of the left and right beam elements adjacent to the disk asCl, Cr, Kl, and Kr, respectively, whose dimensions are all 10 × 10. Then, the damping matrix Cd and stiffness matrix Kd of the disk can be expressed as(5)Cd=C6,6l+C1,1rC6,7l+C1,2rC6,8l+C1,3rC6,9l+C1,4rC6,10l+C1,5rC7,6l+C2,1rC7,7l+C2,2rC7,8l+C2,3rC7,9l+C2,4rC7,10l+C2,5rC8,6l+C3,1rC8,7l+C3,2rC8,8l+C3,3rC8,9l+C3,4rC8,10l+C3,5rC9,6l+C4,1rC9,7l+C4,2rC9,8l+C4,3rC9,9l+C4,4rC9,10l+C4,5rC10,6l+C5,1rC10,7l+C5,2rC10,8l+C5,3rC10,9l+C5,4rC10,10l+C5,5r,(6)Kd=K6,6l+K1,1rK6,7l+K1,2rK6,8l+K1,3rK6,9l+K1,4rK6,10l+K1,5rK7,6l+K2,1rK7,7l+K2,2rK7,8l+K2,3rK7,9l+K2,4rK7,10l+K2,5rK8,6l+K3,1rK8,7l+K3,2rK8,8l+K3,3rK8,9l+K3,4rK8,10l+K3,5rK9,6l+K4,1rK9,7l+K4,2rK9,8l+K4,3rK9,9l+K4,4rK9,10l+K4,5rK10,6l+K5,1rK10,7l+K5,2rK10,8l+K5,3rK10,9l+K5,4rK10,10l+K5,5r.The dissipation energy and potential energy can be expressed as(7)D=12q˙dTCdq˙d,V=12qdTKdqd.Generally, the angular displacement amplitude of torsional vibration is very small. Under this assumption, the following relations are obtained:(8)cosΩt+θz+φ=cosΩt+φcosθz−sinΩt+φsinθz≈cosΩt+φ−θzsinΩt+φ,sinΩt+θz+φ=sinΩt+φcosθz+cosΩt+φsinθz≈sinΩt+φ+θzcosΩt+φ.Substituting equations (1)–(8) into the Lagrange equation, the motion differential equations of the 5-DOFs disk under arbitrary flight conditions can be obtained as follows:(9a)Md+Mduq¨d+Cd+ΩGd+Cdu+Cdbq˙d+Kd+Kdu+Kdbqd=Fdu+Fdb+Qd,(9b)Md=m00000m00000Id00000Id00000Ip,(9c)Gd=0000000000000Ip000−Ip0000000,(9d)Mdu=0000−mesinΩt+φ0000mecosΩt+φ0000000000−mesinΩt+φmecosΩt+φ00me2,(9e)Cdu=−2Ω0000mecosΩt+φ0000mesinΩt+φ000000000000000,(9f)Kdu=0000meΩ2sinΩt+φ0000meΩ2cosΩt+φ000000000000000,(9g)Fdu=meΩ2cosΩt+φmeΩ2sinΩt+φ000,(9h)Cdb=0−2mωB,Z0002mωB,Z00000000IpωB,Y0000−IpωB,X00IpωB,YIpωB,X0,(9i)Kdb=−mωB,Y2+ωB,Z2mωB,XωB,Y−ω˙B,Z000mωB,XωB,Y+ω˙B,Z−mωB,X2+ωB,Z20000000000000000Ipω˙B,Y0,(9j)Fdb=−mX¨B−ωB,ZY˙B+ωB,YZ˙B−mzωB,XωB,Z+ω˙B,Y−mY¨B−ωB,ZX˙B−ωB,XZ˙B−mzωB,YωB,Z−ω˙B,X−IpΩωB,Y−Idω˙B,XIpΩωB,X−Idω˙B,Y−IpωB,XωB,Y−IpθB,Yω˙B,X,where Md and Gd are the traditional mass and gyro matrices without lateral-torsional coupling or maneuvering flight conditions. Fdu is the unbalanced force, and Qd is the generalized external force acting on the disk. Mdu, Cdu, and Kdu denote the coupling items due to the lateral-torsional coupling effect; Cdb, Kdb, and Fdb represent the additional damping, stiffness, and excitation force resulting from maneuvering flight, respectively.In general, the additional stiffness induced by maneuvering flight is very small compared with the inherent stiffness of the rotor system, and the additional damping is related to the rotational angular velocity of the aircraft around each coordinate axis in theOXYZ coordinate system. If the rolling angular velocity of the aircraft is large enough when the aircraft performs the maneuver of rolling, the additional damping may reach the level of the inherent damping of the rotor system.Denoting the additional excitation force asFdb1 when the aircraft performs a hovering flight in the horizontal plane, Fdb2 when the aircraft performs a diving and pulling action in the vertical plane and then(10)Fdb1=−mωB,YZ˙B0−IpΩωB,Y00,(11)Fdb2=0mωB,XZ˙B0IpΩωB,X0.It can be seen that under the two special maneuvering flight conditions, the first two items of the additional excitation forces are the centrifugal force caused by the maneuvering flight action, and the 3rd and 4th items are the maneuvering gyro torque.


### 2.1.2. Shaft

Figure2 shows a schematic diagram of a Timoshenko beam element with both lateral and torsional DOFs. A cylindrical segment with infinite small thickness is taken for analysis at an arbitrary axial position of the beam element. Suppose the thickness of the cylinder is dξ, and the distance between the cylinder and one end of the beam element is ξ, then the kinetic energy of the cylindrical segment can be written as(12)dTe=12ρeAer˙eTr˙e+vBTvB−reTω⌣B2re+2r˙eTω⌣Bre+2vBTr˙e+2vBTω⌣Bredξ+12dIpΩ2+θ˙ze2+2θ˙zeΩ+2θ˙xe+ωB,Xθye+θB,Yθ˙ze+Ω+12dIdθ˙xe+ωB,X2+12dIdθ˙ye+ωB,Y2,(13)re=xeyeze,(14)xeyeθxeθyeθze=N100N20N300N400N1−N2000N3−N4000−N⌣1N⌣2000−N⌣3N⌣400N⌣100N⌣20N⌣300N⌣400000N50000N⌣5qe,(15)qe=x1y1θx1θy1θz1x2y2θx2θy2θz2T,where qe is the nodal displacement vector, Nii=1−5 and N⌣ii=1−5 are the shape functions of the beam element [39], and re denotes the coordinate position of the cylindrical segment in the oxyz coordinate system. The kinetic energy of the beam element can be obtained by integrating dTe along the oz axis. Similar to the derivation process for the disk, the dynamic differential equations of the beam element under arbitrary maneuvering flight conditions can be obtained as(16a)Meq¨e+Ce+ΩGe+Cebq˙e+Ke+Kebqe=Feb+Qe,(16b)Ceb=2ρeAeωB,ZH1+2ρeIeωB,XH2+2ρeIeωB,YH3,(16c)Keb=ρeAeωB,Y2+ωB,Z2H4+ρeAeωB,X2+ωB,Z2H5+ρeAeωB,XωB,Y−ω˙B,ZH6+ρeAeωB,XωB,Y+ω˙B,ZH7+2ρeIeω˙B,XH8,(16d)Feb=ρeAeY˙BωB,Z−Z˙BωB,Y−ω˙B,Yze−X¨B−ωB,XωB,ZzeA1−X˙BωB,Z+Z˙BωB,X+ω˙B,Xze−Y¨B−ωB,YωB,ZzeA1X˙BωB,Z−Z˙BωB,X−ω˙B,Xze+Y¨B+ωB,YωB,ZzeA2Y˙BωB,Z−Z˙BωB,Y−ω˙B,Yze−X¨B−ωB,XωB,ZzeA20Y˙BωB,Z−Z˙BωB,Y−ω˙B,Yze−X¨B−ωB,XωB,ZzeA3−X˙BωB,Z+Z˙BωB,X+ω˙B,Xze−Y¨B−ωB,YωB,ZzeA3X˙BωB,Z−Z˙BωB,X−ω˙B,Xze+Y¨B+ωB,YωB,ZzeA4Y˙BωB,Z−Z˙BωB,Y−ω˙B,Yze−X¨B−ωB,XωB,ZzeA40+ρeIe2ΩωB,X−ω˙B,YB12ΩωB,Y+ω˙B,XB1−2ΩωB,Y+ω˙B,XB22ΩωB,X−ω˙B,YB2−2θB,Yω˙B,X+ωB,XωB,YA52ΩωB,X−ω˙B,YB32ΩωB,Y+ω˙B,XB3−2ΩωB,Y+ω˙B,XB42ΩωB,X−ω˙B,YB42θB,Yω˙B,X+ωB,XωB,YB5,(16e)Ai=∫0leNidξ,Bi=∫0leN⌣idξ,i=1∼5,where Me, Ge, and Ke are the traditional mass, gyro, and stiffness matrices, respectively (their expressions are omitted due to space limitations and refer Friswell et al. [39]). Ce denotes the external linear damping, and Qe represents the generalized external force. Ceb, Keb, and Feb are the additional damping, stiffness, and excitation force due to maneuvering flight, respectively. The matrices H1–H8 are related to the shape functions and geometric dimension of the beam (see the appendix).Figure 2
Diagram of Timoshenko beam element.Based on finite element theory, the dynamic differential equations of the rotor system under arbitrary maneuvering flight conditions can be written as follows by rearranging the matrices of the disk elements and beam elements [40].(17)M+Muq¨+C+ΩG+Cu+Cbq˙+K+Ku+Kbq=Fu+Fb+Q,where M, C, G, and K are the inertia, damping, gyro moment, and stiffness matrices of the rotor system without lateral-torsional coupling effect or maneuvering flight conditions. Mu, Cu, and Ku denote the additional inertia, damping, and stiffness due to lateral-torsional coupling effect, and Cb, Kb, and Fb represent the additional damping, stiffness, and excitation force vector caused by maneuvering flight. Fu is the unbalanced force, and Q represents the generalized external force acting on the rotor, including the rub-impact force and the nonlinear forces of the supports.


## 2.1.1. Disk Element

The coordinate system of an arbitrary disk of the rotor system is shown in Figure1, where OXYZ is the ground coordinate system for describing the spatial position, velocity, and acceleration of the aircraft, and oxyz is the rotor coordinate system for describing the translation and rotation of the rotor relative to the aircraft. Let us suppose that the engine rotor is installed on the aircraft such that one end coincides with the center of gravity of the aircraft, and thus, oxyz is also the relative coordinate system of the aircraft. The generalized displacement vector of the five degrees-of-freedom (5-DOFs) disk is given as(1)qd=xyθxθyθzT.Figure 1
Coordinate system of disk.Then, the kinetic energy of the disk can be expressed as [13, 15] follows:(2)T=12mr˙Tr˙+vBTvB−rTω⌣B2r+2r˙Tω⌣Br+2vBTr˙+2vBTω⌣Br+12Idθ˙x+ωB,X2+12Idθ˙y+ωB,Y2+12IpΩ2+θ˙z2+2θ˙zΩ+2θy+θB,Yθ˙x+ωB,Xθ˙z+Ω,(3)r=x+ecosΩt+θz+φy+esinΩt+θz+φz,vB=X˙BY˙BZ˙B,ω⌣B=0−ωB,ZωB,YωB,Z0−ωB,X−ωB,YωB,X0,where m denotes the quality of the disk, Id and Ip are the diameter and polar moment of inertia, respectively, Ω is the rotational speed of the rotor, r denotes the coordinate position of the disk in the oxyz coordinate system, e is the eccentricity distance, φ represents the initial phase angle, vB and ωB are the velocity and angular velocity of the aircraft in the OXYZ coordinate system, X˙B, Y˙B, Z˙B, and ωB,X, ωB,Y, and ωB,Z denote the velocity and angular velocity of the aircraft in the OXYZ coordinate system, respectively, and θB,Y is the angular displacement of the aircraft around the OY axis, and(4)θy+θB,Y≈θy,θy+θB,Y′≈θ˙y+ωB,Y.Denoting the damping matrix and stiffness matrix of the left and right beam elements adjacent to the disk asCl, Cr, Kl, and Kr, respectively, whose dimensions are all 10 × 10. Then, the damping matrix Cd and stiffness matrix Kd of the disk can be expressed as(5)Cd=C6,6l+C1,1rC6,7l+C1,2rC6,8l+C1,3rC6,9l+C1,4rC6,10l+C1,5rC7,6l+C2,1rC7,7l+C2,2rC7,8l+C2,3rC7,9l+C2,4rC7,10l+C2,5rC8,6l+C3,1rC8,7l+C3,2rC8,8l+C3,3rC8,9l+C3,4rC8,10l+C3,5rC9,6l+C4,1rC9,7l+C4,2rC9,8l+C4,3rC9,9l+C4,4rC9,10l+C4,5rC10,6l+C5,1rC10,7l+C5,2rC10,8l+C5,3rC10,9l+C5,4rC10,10l+C5,5r,(6)Kd=K6,6l+K1,1rK6,7l+K1,2rK6,8l+K1,3rK6,9l+K1,4rK6,10l+K1,5rK7,6l+K2,1rK7,7l+K2,2rK7,8l+K2,3rK7,9l+K2,4rK7,10l+K2,5rK8,6l+K3,1rK8,7l+K3,2rK8,8l+K3,3rK8,9l+K3,4rK8,10l+K3,5rK9,6l+K4,1rK9,7l+K4,2rK9,8l+K4,3rK9,9l+K4,4rK9,10l+K4,5rK10,6l+K5,1rK10,7l+K5,2rK10,8l+K5,3rK10,9l+K5,4rK10,10l+K5,5r.The dissipation energy and potential energy can be expressed as(7)D=12q˙dTCdq˙d,V=12qdTKdqd.Generally, the angular displacement amplitude of torsional vibration is very small. Under this assumption, the following relations are obtained:(8)cosΩt+θz+φ=cosΩt+φcosθz−sinΩt+φsinθz≈cosΩt+φ−θzsinΩt+φ,sinΩt+θz+φ=sinΩt+φcosθz+cosΩt+φsinθz≈sinΩt+φ+θzcosΩt+φ.Substituting equations (1)–(8) into the Lagrange equation, the motion differential equations of the 5-DOFs disk under arbitrary flight conditions can be obtained as follows:(9a)Md+Mduq¨d+Cd+ΩGd+Cdu+Cdbq˙d+Kd+Kdu+Kdbqd=Fdu+Fdb+Qd,(9b)Md=m00000m00000Id00000Id00000Ip,(9c)Gd=0000000000000Ip000−Ip0000000,(9d)Mdu=0000−mesinΩt+φ0000mecosΩt+φ0000000000−mesinΩt+φmecosΩt+φ00me2,(9e)Cdu=−2Ω0000mecosΩt+φ0000mesinΩt+φ000000000000000,(9f)Kdu=0000meΩ2sinΩt+φ0000meΩ2cosΩt+φ000000000000000,(9g)Fdu=meΩ2cosΩt+φmeΩ2sinΩt+φ000,(9h)Cdb=0−2mωB,Z0002mωB,Z00000000IpωB,Y0000−IpωB,X00IpωB,YIpωB,X0,(9i)Kdb=−mωB,Y2+ωB,Z2mωB,XωB,Y−ω˙B,Z000mωB,XωB,Y+ω˙B,Z−mωB,X2+ωB,Z20000000000000000Ipω˙B,Y0,(9j)Fdb=−mX¨B−ωB,ZY˙B+ωB,YZ˙B−mzωB,XωB,Z+ω˙B,Y−mY¨B−ωB,ZX˙B−ωB,XZ˙B−mzωB,YωB,Z−ω˙B,X−IpΩωB,Y−Idω˙B,XIpΩωB,X−Idω˙B,Y−IpωB,XωB,Y−IpθB,Yω˙B,X,where Md and Gd are the traditional mass and gyro matrices without lateral-torsional coupling or maneuvering flight conditions. Fdu is the unbalanced force, and Qd is the generalized external force acting on the disk. Mdu, Cdu, and Kdu denote the coupling items due to the lateral-torsional coupling effect; Cdb, Kdb, and Fdb represent the additional damping, stiffness, and excitation force resulting from maneuvering flight, respectively.In general, the additional stiffness induced by maneuvering flight is very small compared with the inherent stiffness of the rotor system, and the additional damping is related to the rotational angular velocity of the aircraft around each coordinate axis in theOXYZ coordinate system. If the rolling angular velocity of the aircraft is large enough when the aircraft performs the maneuver of rolling, the additional damping may reach the level of the inherent damping of the rotor system.Denoting the additional excitation force asFdb1 when the aircraft performs a hovering flight in the horizontal plane, Fdb2 when the aircraft performs a diving and pulling action in the vertical plane and then(10)Fdb1=−mωB,YZ˙B0−IpΩωB,Y00,(11)Fdb2=0mωB,XZ˙B0IpΩωB,X0.It can be seen that under the two special maneuvering flight conditions, the first two items of the additional excitation forces are the centrifugal force caused by the maneuvering flight action, and the 3rd and 4th items are the maneuvering gyro torque.


## 2.1.2. Shaft

Figure2 shows a schematic diagram of a Timoshenko beam element with both lateral and torsional DOFs. A cylindrical segment with infinite small thickness is taken for analysis at an arbitrary axial position of the beam element. Suppose the thickness of the cylinder is dξ, and the distance between the cylinder and one end of the beam element is ξ, then the kinetic energy of the cylindrical segment can be written as(12)dTe=12ρeAer˙eTr˙e+vBTvB−reTω⌣B2re+2r˙eTω⌣Bre+2vBTr˙e+2vBTω⌣Bredξ+12dIpΩ2+θ˙ze2+2θ˙zeΩ+2θ˙xe+ωB,Xθye+θB,Yθ˙ze+Ω+12dIdθ˙xe+ωB,X2+12dIdθ˙ye+ωB,Y2,(13)re=xeyeze,(14)xeyeθxeθyeθze=N100N20N300N400N1−N2000N3−N4000−N⌣1N⌣2000−N⌣3N⌣400N⌣100N⌣20N⌣300N⌣400000N50000N⌣5qe,(15)qe=x1y1θx1θy1θz1x2y2θx2θy2θz2T,where qe is the nodal displacement vector, Nii=1−5 and N⌣ii=1−5 are the shape functions of the beam element [39], and re denotes the coordinate position of the cylindrical segment in the oxyz coordinate system. The kinetic energy of the beam element can be obtained by integrating dTe along the oz axis. Similar to the derivation process for the disk, the dynamic differential equations of the beam element under arbitrary maneuvering flight conditions can be obtained as(16a)Meq¨e+Ce+ΩGe+Cebq˙e+Ke+Kebqe=Feb+Qe,(16b)Ceb=2ρeAeωB,ZH1+2ρeIeωB,XH2+2ρeIeωB,YH3,(16c)Keb=ρeAeωB,Y2+ωB,Z2H4+ρeAeωB,X2+ωB,Z2H5+ρeAeωB,XωB,Y−ω˙B,ZH6+ρeAeωB,XωB,Y+ω˙B,ZH7+2ρeIeω˙B,XH8,(16d)Feb=ρeAeY˙BωB,Z−Z˙BωB,Y−ω˙B,Yze−X¨B−ωB,XωB,ZzeA1−X˙BωB,Z+Z˙BωB,X+ω˙B,Xze−Y¨B−ωB,YωB,ZzeA1X˙BωB,Z−Z˙BωB,X−ω˙B,Xze+Y¨B+ωB,YωB,ZzeA2Y˙BωB,Z−Z˙BωB,Y−ω˙B,Yze−X¨B−ωB,XωB,ZzeA20Y˙BωB,Z−Z˙BωB,Y−ω˙B,Yze−X¨B−ωB,XωB,ZzeA3−X˙BωB,Z+Z˙BωB,X+ω˙B,Xze−Y¨B−ωB,YωB,ZzeA3X˙BωB,Z−Z˙BωB,X−ω˙B,Xze+Y¨B+ωB,YωB,ZzeA4Y˙BωB,Z−Z˙BωB,Y−ω˙B,Yze−X¨B−ωB,XωB,ZzeA40+ρeIe2ΩωB,X−ω˙B,YB12ΩωB,Y+ω˙B,XB1−2ΩωB,Y+ω˙B,XB22ΩωB,X−ω˙B,YB2−2θB,Yω˙B,X+ωB,XωB,YA52ΩωB,X−ω˙B,YB32ΩωB,Y+ω˙B,XB3−2ΩωB,Y+ω˙B,XB42ΩωB,X−ω˙B,YB42θB,Yω˙B,X+ωB,XωB,YB5,(16e)Ai=∫0leNidξ,Bi=∫0leN⌣idξ,i=1∼5,where Me, Ge, and Ke are the traditional mass, gyro, and stiffness matrices, respectively (their expressions are omitted due to space limitations and refer Friswell et al. [39]). Ce denotes the external linear damping, and Qe represents the generalized external force. Ceb, Keb, and Feb are the additional damping, stiffness, and excitation force due to maneuvering flight, respectively. The matrices H1–H8 are related to the shape functions and geometric dimension of the beam (see the appendix).Figure 2
Diagram of Timoshenko beam element.Based on finite element theory, the dynamic differential equations of the rotor system under arbitrary maneuvering flight conditions can be written as follows by rearranging the matrices of the disk elements and beam elements [40].(17)M+Muq¨+C+ΩG+Cu+Cbq˙+K+Ku+Kbq=Fu+Fb+Q,where M, C, G, and K are the inertia, damping, gyro moment, and stiffness matrices of the rotor system without lateral-torsional coupling effect or maneuvering flight conditions. Mu, Cu, and Ku denote the additional inertia, damping, and stiffness due to lateral-torsional coupling effect, and Cb, Kb, and Fb represent the additional damping, stiffness, and excitation force vector caused by maneuvering flight. Fu is the unbalanced force, and Q represents the generalized external force acting on the rotor, including the rub-impact force and the nonlinear forces of the supports.


## 2.2. Nonlinear Forces of the Supports

Based on Hertz contact theory, the nonlinear forces of a rolling ball bearing can be written as [41](18a)fxbearing=−kn∑j=1NbuθjξHuθjcosθj,fybearing=−kn∑j=1NbuθjξHuθjsinθj,(18b)θj=2πj−1Nb+ωbmt,(18c)uθj=xir−xorcosθj+yir−yorsinθj−γ2,(18d)Huθj=0,uθj≤0,uθj,uθj>0,(18e)ωbm=ωinr1+ωoutr2r1+r2,where the superscripts “ir” and “or” denote the inner and outer rings of the bearing, respectively. kn is the contact stiffness between the roller and the ring, Nb is the number of rollers, and θj is the circumferential position of the jth roller at time t. uθj represents the elastic contact deformation of the jth roller in the radial direction. γ is the radial clearance. r1 and r2 denote the radius of the inner and outer rings, respectively, ωin and ωout are the rotational speeds of the inner and outer rings, and ωbm is the rotational speed of the revolution of the rollers. Generally, ξ=3/2 for rolling ball bearings.As for SFD, during the process of whirling of the rotor, the nonlinear oil film forces of the SFD can be written as [42, 43](19a)fxsfd=−μRL3c2x2+y21/2xε˙I2+εψ˙I1+yε˙I1+εψ˙I3,fysfd=−μRL3c2x2+y21/2yε˙I2+εψ˙I1−xε˙I1+εψ˙I3,(19b)ε=x2+y2c,(19c)ε˙=xx˙+yy˙cx2+y2,(19d)ψ˙=yx˙−xy˙x2+y2,(19e)tanψ=xy,where x and y are the eccentric positions of the journal, L, R, and c denote the length, radius, and radial clearance of the SFD, respectively, μ represents the dynamic viscosity, and the Sommerfeld integral Ijj=1−3 can be expressed as(20)I1=∫θ1θ2cos2θ1+εcosθ3dθ,I2=∫θ1θ2sinθcosθ1+εcosθ3dθ,I3=∫θ1θ2sin2θ1+εcosθ3dθ,where θ1=tan−1−ε˙/εψ˙, θ2=θ1+π.


## 2.3. Rub-Impact Forces

Considering the deformation and energy consumption in the process of rubbing, the normal rub-impact forceFN can be expressed as follows based on Hertz contact theory [44]:(21a)FN=kδq+cδ˙,(21b)c=1.5αkδq,(21c)δ=xi2+yi2+rin−rout,(21d)δ˙=x˙i2+y˙i2,(21e)k=1.5681−v12E1+1−v22E2−4/5rin2δ2/5m2δ˙21/5,where c is the viscous damping coefficient in the normal direction, α denotes the energy dissipation factor, δ is the invasion displacement between the rotor and stator (when δ>0, rub-impact occurs), δ˙ is the relative velocity of the rotor and stator in the vertical direction of the contact surface, xi and yi denote the displacement of node i of the rotor, and x˙i, y˙i the velocity. rin and rout represent the outer diameter of the rotor and the inner diameter of the casing, respectively, and k is the rubbing stiffness coefficient. v1, v2, E1, and E2 are the Poisson’s ratio and elastic moduli of the rotor and stator, respectively, and m2 is the equivalent quality of the rotor when rub-impact occurs.Based on Coulomb friction theory, the tangential friction due to rub-impact can be expressed as [45](22a)Fτ=ΘμfFN,(22b)Θ=−1,ψ<0,0,ψ=0,1,ψ>0,(22c)ψ=uv˙−vu˙v2+u2+rinω,where μf is the friction coefficient, Θ is employed to judge the direction of the tangential friction, and u and v denote the displacements of the rub-impact node in x and y directions, respectively.Figure3 shows the schematic diagram of rub-impact forces between the rotor and stator considering the deformation of the stator under maneuvering flight, where o is the initial centroid of the rotor and stator, and o′(x′, y′) the centroid of the stator after deformation due to maneuvering overload, and o1(u, v) the centroid of the rotor when rub-impact occurs. The rotor rotates counterclockwise. Then, the rub-impact forces in the rotor coordinate system can be written as(23)Fx=−FNcosθ+Fτsinθ,Fy=−FNsinθ−Fτcosθ,Mz=Fτrin,where Mz is the torque acting on the rotor due to rub-impact, and θ=arctanv−y′/u−x′. The translational displacement of the stator x′ and y′ can be determined according to the maneuvering overload and the deformation stiffness of the stator.Figure 3
Schematic diagram of rub-impact forces.


## 3. Calculation Results and Discussion

The nonlinear forces of the supports and the rub-impact forces are substituted into the vectorQ in equation (17) according to the position of the bearing and rub-impact nodes, and the dynamic finite element model of the rub-impact rotor system with lateral-torsional coupling vibration under arbitrary maneuvering flight conditions is then obtained, considering the nonlinear forces of the SFD and intermediate bearing. An implicit numerical integral method combining the Newmark-β and Newton–Raphson methods is employed to solve the resulting nonlinear equations.


### 3.1. Physical Model Introduction

A schematic diagram of a dual-rotor system is shown in Figure4. The rotor system is composed of two shafts with four disks and four supports. Disks 1 and 4 are installed on the inner rotor to simulate the low-pressure compressor and low-pressure turbine of a real aeroengine, whereas disks 2 and 3 are located on the outer rotor to simulate the high-pressure compressor and high-pressure turbine. A supporting scheme combining an elastic support, a rolling ball bearing, and SFDs is adopted for supports I, II, and III, and support IV is the intermediate support, with the inner ring connected to the inner rotor and the outer ring linked to the outer rotor. The radial stiffness of the elastic supports is listed in Table 1, the inertia properties and unbalance of the disks are given in Table 2, and the parameters of the intermediate bearing and SFDs are presented in Tables 3 and 4 , respectively. The density of the shafts is 7810 kg/m3, the elastic modulus is 196 GPa, and the shear modulus is 75.5 GPa. The disks are made of C45E. The rotational speed ratio of the dual-rotor system is λ=ω2/ω1=1.6, where ω1 and ω2 are the rotational speeds of the inner and outer rotors, respectively.Figure 4
Schematic diagram of dual-rotor system.Table 1
Radial stiffness of the elastic supports.
NameIIIIIIStiffness (N/m)2.9 × 1074.4 × 1072 × 107Table 2
Inertial attributes and unbalanced configuration of the disks.
Inner rotorOuter rotorDisk 1Disk 4Disk 2Disk 3Mass (kg)2.412.383.171.57Polar moment of inertia (kg·m2)0.00820.00820.01460.0059Unbalance (10−5 kg·m)221.53Table 3
Parameters of the intermediate bearing.
NameRadius of inner ring (mm)Radius of outer ring (mm)Number of rollersRadial clearance (μm)Value9.3514.1596Table 4
Parameters of SFDs.
NameInner rotorOuter rotorSFD ISFD IISFD IIIRadius (mm)251835Length (mm)131314Radial clearance (mm)0.250.180.35Dynamic viscosity (10−2 Pa·s)1.0752Table5 presents the node information of the finite element model after meshing. The model contains 54 elements with 46 beam elements, 4 disk elements, and 4 supporting elements.Table 5
Node information of the finite element model.
Node numberAxial location (m)Bearing/diskElement numberOuter diameter (m)Inner diameter (m)10.110.018020.1320.018030.1630.018040.1940.018050.2250.018060.2547960.018070.2887970.018080.3227980.018090.35879Disk 190.0180100.38369100.0180110.40859110.0180120.43349120.0180130.44479Bearing 1130.0180140.49140.0220150.56150.0220160.65025160.0220170.75298170.0220180.82180.0220190.95844190.0220201.03200.0220211.06517Bearing 4210.0220221.087220.0220231.1230.0220241.14Disk 4240.0220251.17250.0220261.19260.0220271.21088270.0170281.21488Bearing 2280.0140291.22290.0140301.23038End of the inner rotor310.642300.0350.03320.655310.0350.03330.6793Bearing 3320.0350.03340.6865330.0380.03350.7117340.0380.03360.7369350.0380.03370.7621Disk 2360.0380.03380.80784370.0380.03390.85358380.0380.03400.89932390.0380.03410.94506400.0380.03420.96Disk 3410.0380.03430.98420.0380.03441430.070.03451.02440.060.03461.04450.060.03471.0643Bearing 4460.060.03481.0738Based on the dynamic differential equations of the rotor system under arbitrary maneuvering flight conditions in Section2.1, the motion differential equations of the inner and outer rotors can be written as(24)ML+MuLq˙L+CL+ΩGL+CuL+CbLq˙L+KL+KuL+KbLqL=FuL+FbL+fL,bearing+fL,sfd+FL,rub,(25)MH+MuHq˙H+CH+ΩGH+CuH+CbHq˙H+KH+KuH+KbHqH=FuH+FbH+fH,bearing+fH,sfd+FH,rub,where the superscripts “L” and “H” denote the inner and outer rotors.Denoting the node number of the inner rotor and outer linked to the intermediate bearing asξ and η, respectively. qξbearing and qηbearing represent the translational displacement of node ξ and node η. Then, the displacement boundary conditions are obtained according to equation (18a)–(18e), and the coupling relationship between the inner and outer rotors at the position of the intermediate bearing is given as(26)fξL,bearing=−fηH,bearing=fbearingqξbearing,qηbearing.Equations (24)–(26) are the motion differential equations of the dual rotor system.Linearize the dynamics model of the dual rotor system and solve the eigenvalues for the linearized equations of motion using the state space method, and the Campbell diagram is obtained as shown in Figure5. The radial stiffness of the intermediate bearing is set to 2.5 × 106 N/m. The first three orders of critical speed with the rotational speed of the inner rotor ω1 as excitation frequency are 233.4 rad/s, 403.3 rad/s, and 986.1 rad/s. While the first three orders of critical speed with ω2 as excitation frequency are 143 rad/s, 248.1 rad/s, and 610.4 rad/s. The natural frequency of the 1st torsional vibration is 111.7 Hz. Besides, the mode shapes corresponding to the critical speeds are also shown in Figure 5.Figure 5
Campbell diagram of the dual-rotor system.To investigate the effects on the rotor dynamic characteristics of maneuvering flight, the actions of level flight and hovering flight are taken as examples to study. Besides, the magnitude of maneuvering overload is an important parameter affecting the vibration response of the rotor system. In general, the maneuvering overload does not exceed 9 g due to the tolerance limit of the human body, where g represents the acceleration due to gravity. However, the future high maneuverability military drone is selected to study in this article, and the maneuvering overload is set to 10 g. The maneuvering flight parameters are listed in Table6, and all flight parameters not explicitly mentioned are set to 0.Table 6
Parameters of the level and hovering flight.
ActionsZ˙B (m/s)ωB,Y(rad/s)Overload (g)Level flight9800Hovering flight98110


### 3.2. Implicit Newmark-β Method

The differential equations of motion of the dual rotor system can be solved using the numerical integration method. The DOFs of the rotor system are divided into two parts according to whether they are subjected to nonlinear forces or not,q=qkqjT, where qk and qj represent the linear and nonlinear DOFs vector, respectively.An implicit numerical integral method combining the Newmark-β and Newton–Raphson methods are adopted in this work. The flow chart of the solution is shown in Figure 6. First, the vibration response of the rotor system at time tn are qkn, qjn, q˙kn, q˙jn, q¨kn, and q¨jn, assuming that the displacement response of nonlinear DOFs at tn+1 is qjn+1; then, the q˙jn+1 can be obtained by Newmark hypothesis. Substituting qjn+1 and q˙jn+1 into equations (18a)–(18e) and (19a)–(19e), the nonlinear forces can be solved. Then, the equations transform into linear equations, which can be solved by Newmark-β method. Through Newton–Raphson iterative calculation, the convergence solution at tn+1 can be obtained. Then save the vibration response at tn+1 and start the calculation of the next time interval. For the Newmark-β method, we set α=0.3 and β=0.6.Figure 6
Flow chart of the implicit newmark-β method.


### 3.3. Response Analysis

Compared with the outer rotor, the inner rotor is slender in configuration and thus suffers from a higher vibration level in the torsional direction, and the spectrum characteristics of disks 3 and 4 in the lateral direction are basically consistent. Thus, the vibration response characteristics of disk 4 are the primary focus in this work.We assume that rub-impact occurs at the position of disk 4. The lateral vibration amplitude of disk 4 without considering rub-impact faults under the level flight and hovering flight conditions is calculated to determine the relevant parameters in the rubbing response analysis, such as the calculating rotational speeds of the rotor and the initial radial clearance between the rotor and stator. Figure7 shows that the lateral vibration amplitude of disk 4 varies with the rotational speed under the conditions of level flight and hovering flight. As the speed of the inner rotor increases from 10 rad/s to 550 rad/s with interval of 10 rad/s, three peaks in vibration amplitude are observed corresponding to the rotational speeds of 140 rad/s, 250 rad/s, and 410 rad/s, respectively. Compared with the Campbell diagram of the dual rotor system (Figure 5), the error of the value of the critical speeds is due to the selection of the rotational speed interval for the vibration response calculation and the linearization for the rotor system when calculating the eigenvalues. The rotational speeds corresponding to the three peaks of the vibration response from small to large are the first-order–second-order critical speeds with excitation frequency ω2, and second-order critical speed with excitation frequency ω1, respectively. The peak of the first-order critical speed with excitation frequency ω1 is not observed because the rotational speed interval for the vibration response calculation and the first-order critical speed with excitation frequency ω1 is close to the second-order critical speed with excitation frequency ω2.Figure 7
Vibration response of disk 4 with respect to rotational speed. (a) Horizontal direction. (b) Vertical direction.
(a)(b)Besides, the action of hovering flight causes a static offset of the rotor in both the horizontal and vertical directions. In the horizontal direction, the rotor system is subjected to a constant centrifugal force generated by maneuvering flight, and the static offset does not change with the rotational speed. While an additional gyroscopic moment related to the rotational speed is generated in the vertical direction, which causes the static offset to vary with the rotational speed.The waterfall spectrum plots for the lateral vibration response of disk 4 during level flight and hovering flight conditions are shown in Figure8. Under the level flight condition, the frequency components of the vibration response mainly include the rotational speeds of the inner and outer rotors, i.e., ω1, ω2, and the combinational frequency components due to the coupling of the inner and outer rotors through the intermediate bearing, such as 2ω2−ω1, 3ω2−2ω1, and 4ω2−3ω1, and the like. The two vibration peaks with frequency ω2 correspond to the first- and second-order critical speeds with excitation frequency ω2 (140 rad/s and 250 rad/s), respectively. While the two vibration peaks with frequency ω1 correspond to the first- and second-order critical speeds with excitation frequency ω1 (240 rad/s and 410 rad/s). Additionally, a subharmonic vibration frequency component ωs1 induced by the bearing nonlinearity can be observed. Compared with the condition of level flight, the combinational frequency components of the inner and outer rotors are weakened, the amplitude of the subharmonic vibration is clearly enhanced, and a new subharmonic component ωs2 appears during hovering flight. The subharmonic vibration causes significant damage to the bearing of the rotor system in engineering and therefore requires serious consideration.Figure 8
Waterfall spectrum plots for the lateral vibration of disk 4 in vertical direction. (a) Level flight. (b) Hovering flight.
(a)(b)Figure9 shows the waterfall spectrum plots for the torsional vibration response of disk 4 during level flight and hovering flight. In this work, the amplitude of the torsional vibration is defined as the relative torsional angle of the two disks on the inner or outer rotor, that is, the torsional vibration amplitude of disk 4 is the difference between the torsional angles of disks 1 and 4. Under the level flight condition, the main frequency components of torsional vibration are ω2−ω1, 2ω2−2ω1, 3ω2−3ω1, and so on. However, significant changes take place in the frequency components during hovering flight; the subharmonic vibration component in the lateral vibration amplified by maneuvering flight is transmitted to the torsional direction due to the coupling effect of lateral and torsional vibrations, and the main frequency components change to ω1, ω2, ω2−ω1, the subharmonic vibration frequency ωs1, ωs2, and the combinational frequency ω1+ωs1.Figure 9
Waterfall spectrum plots for the torsional vibration of disk 4. (a) Level flight. (b) Hovering flight.
(a)(b)Assuming that the rub-impact occurs at the second-order critical speed with excitation frequencyω2 (250 rad/s). Thus, to analyze the rub-impact response, the inner rotor speeds of 200 rad/s–300 rad/s are examined in detail, and the rotational speed interval is set to 2 rad/s. The relevant parameters of the casing are listed in Table 7. The friction coefficient is 0.1, and the energy dissipation coefficient is 0.16 s/m. The initial radial clearance between the rotor and casing at the position of disk 4 is 0.5 × 10−4 m under the condition of level flight. According to the static offset of the rotor and casing due to maneuvering flight, the initial clearance during hovering is set to 1 × 10−4 m.Table 7
Parameters of the casing.
NameMass (kg)Elastic modulus (GPa)Poisson ratioDeformation stiffness (N/m)Value6001960.258 × 10−7Figure10 compares the lateral vibration amplitude of disk 4 with or without rub-impact under the two flight conditions. In the level flight condition, the rotor whirls around the central axis of the casing. The whirling range is limited by the casing, and rub-impact may occur at an arbitrary phase of the full cycle. Thus, the vibration amplitudes of disk 4 with rub-impact are significantly different from those without rub-impact (see Figure 10(a)). In the case of hovering flight, the additional centrifugal force and gyroscopic moment are generated on the rotor, while the casing is only subjected to centrifugal force. Besides, the deformation stiffness of the casing is different from that of the rotor in general. These two factors cause different static offsets of the rotor and casing. Thus, full-cycle rub-impact does not occur during hovering flight, but rub-impact will take place at local positions where the radial clearance of the rotor and casing is relatively small. According to the simulation results in Figure 10(b), the casing does not completely limit the whirling range of the rotor, and the vibration amplitude of disk 4 increases due to rub-impact near the critical speed of the system.Figure 10
Lateral vibration response of disk 4 with or without rub-impact. (a) Level flight. (b) Hovering flight.
(a)(b)The comparison of the torsional vibration amplitude of disk 4 varying with the rotational speed is shown in Figure11, when considering rub-impact or not under the two flight conditions. It can be observed that the torsional vibration amplitude increases near the critical speed of the system due to the tangential force produced by rub-impact and the lateral-torsional coupling effect, and the vibration amplitude of rub-impact is much larger than that without the rub-impact in the torsional direction.Figure 11
Torsional vibration response of disk 4 with or without rub-impact. (a) Level flight. (b) Hovering flight.
(a)(b)Figure12 shows the waterfall spectrum plots for the lateral vibration response of disk 4 during level flight. Comparing Figure 12(a) with Figure 12(b), the frequency components of the vibration response change significantly when rub-impact occurs. The frequency component 2ω1, and some combinational frequencies are produced by the rub-impact, such as 2ω1−ω2, 3ω2−4ω1, 2ω2−ω1, and 3ω2−2ω1, etc.Figure 12
Waterfall spectrum plots for the lateral vibration of disk 4 during level flight. (a) Without rub-impact. (b) Rub-impact.
(a)(b)The waterfall spectrum plots are shown in Figure13 for the lateral vibration response of disk 4 during hovering flight. Compared with the level flight condition, rub-impact produces only small changes in the frequency components of the rotor system during hovering flight. It has been analyzed that the rub-impact can only occur at the local position where the radial clearance is relatively small during hovering flight; thus, the frequency component of the rubbing response is relatively simple, and only an additional frequency component of ωr is observed, which varies with the rotational speed of the rotor between 0.75ω1 and 0.9ω1.Figure 13
Waterfall spectrum plots for the lateral vibration of disk 4 during hovering flight. (a) Without rub-impact. (b) Rub-impact.
(a)(b)Figures14 and 15 show the waterfall spectrum plots for the torsional vibration response of disk 4 with or without rub-impact under the two flight conditions. The torsional vibration amplitude of the rotor increases significantly in the presence to rub-impact, and the vibration frequency component is different in the two flight conditions. During level flight, the frequency components mainly include the zero frequency and ω2−ω1 when rub-impact does not occur. While some combinational frequencies are observed due to rub-impact, such as 2ω2−2ω1, 3ω2−3ω1, 4ω2−4ω1, and so on. Under the condition of hovering flight, the frequency components include the zero frequency, ω1, ω2, and ω2−ω1 without rub-impact. When rub-impact occurs, the frequency components change to the zero frequency, ω2−ω1, ω2, and 2ω2.Figure 14
Waterfall spectrum plots for the torsional vibration of disk 4 during level flight. (a) Without rub-impact. (b) Rub-impact.
(a)(b)Figure 15
Waterfall spectrum plots for the torsional vibration of disk 4 during hovering flight. (a) Without rub-impact. (b) Rub-impact.
(a)(b)Figures16 and 17 show the orbit and Poincaré maps of disk 4 under different rotational speeds near the occurrence of rubbing during level flight, and the orbit and Poincaré maps during hovering flight are shown in Figures 18 and 19 . Before and after the occurrence of rub-impact, the shapes of the orbits are similar under the two flight conditions (see Figures 16(a), 16(b), 16(f), 18(a), 18(b), and 18(f)); however, the orbits have clearly different shapes when rub-impact occurs (see Figures 16(c)–16(e) and 18(c)–18(e)). Besides, under both flight conditions, the motion of the rotor system is in the state of quasiperiodic before and after rub-impact due to the nonlinearity of the intermediate bearing (see Figures 17(a), 17(b), 17(f), 19(a), 19(b), and 19(f)). When rub-impact occurs, the motion of the rotor system tends to chaos with the increase in the degree of rub-impact (see Figures 17(c)–17(e) and 19(c)–19(e)).Figure 16
Orbit of disk 4 with rub-impact during level flight. (a) 220 rad/s. (b) 226 rad/s. (c) 236 rad/s. (d) 252 rad/s. (e) 256 rad/s. (f) 270 rad/s.
(a)(b)(c)(d)(e)(f)Figure 17
Poincaré maps of disk 4 with rub-impact during level flight. (a) 220 rad/s. (b) 226 rad/s. (c) 236 rad/s. (d) 252 rad/s. (e) 256 rad/s. (f) 270 rad/s.
(a)(b)(c)(d)(e)(f)Figure 18
Orbit of disk 4 with rub-impact during hovering flight. (a) 220 rad/s. (b) 226 rad/s. (c) 236 rad/s. (d) 252 rad/s. (e) 256 rad/s. (f) 270 rad/s.
(a)(b)(c)(d)(e)(f)Figure 19
Poincaré maps of disk 4 with rub-impact during hovering flight. (a) 220 rad/s. (b) 226 rad/s. (c) 236 rad/s. (d) 252 rad/s. (e) 256 rad/s. (f) 270 rad/s.
(a)(b)(c)(d)(e)(f)Figure20 shows the bifurcation diagrams for the lateral vibration response of disk 4 during level flight and hovering flight. The full-cycle rubbing clearly has a greater effect on the periodicity of the motion of the rotor under the level flight condition. In the range of rotation speeds considered herein, the motion of the rotor in the lateral direction moves from quasiperiodic to chaotic and back to quasiperiodic. During hovering flight, the motion of the rotor is relatively regular, except for the region of the critical speed of the system.Figure 20
Bifurcation diagrams for the lateral vibration of disk 4 with rub-impact. (a) Level flight. (b) Hovering flight.
(a)(b)The bifurcation diagrams for the torsional vibration response of disk 4 during level flight and hovering flight are shown in Figure21. Similar to the lateral vibration, the motion of the rotor is quasiperiodic in the torsional direction before and after rub-impact. When rub-impact occurs, the motion becomes chaotic. Compared with the condition of level flight, the motion exhibits better periodicity during hovering flight.Figure 21
Bifurcation diagrams for the torsional vibration of disk 4 with rub-impact. (a) Level flight. (b) Hovering flight.
(a)(b)


## 3.1. Physical Model Introduction

A schematic diagram of a dual-rotor system is shown in Figure4. The rotor system is composed of two shafts with four disks and four supports. Disks 1 and 4 are installed on the inner rotor to simulate the low-pressure compressor and low-pressure turbine of a real aeroengine, whereas disks 2 and 3 are located on the outer rotor to simulate the high-pressure compressor and high-pressure turbine. A supporting scheme combining an elastic support, a rolling ball bearing, and SFDs is adopted for supports I, II, and III, and support IV is the intermediate support, with the inner ring connected to the inner rotor and the outer ring linked to the outer rotor. The radial stiffness of the elastic supports is listed in Table 1, the inertia properties and unbalance of the disks are given in Table 2, and the parameters of the intermediate bearing and SFDs are presented in Tables 3 and 4 , respectively. The density of the shafts is 7810 kg/m3, the elastic modulus is 196 GPa, and the shear modulus is 75.5 GPa. The disks are made of C45E. The rotational speed ratio of the dual-rotor system is λ=ω2/ω1=1.6, where ω1 and ω2 are the rotational speeds of the inner and outer rotors, respectively.Figure 4
Schematic diagram of dual-rotor system.Table 1
Radial stiffness of the elastic supports.
NameIIIIIIStiffness (N/m)2.9 × 1074.4 × 1072 × 107Table 2
Inertial attributes and unbalanced configuration of the disks.
Inner rotorOuter rotorDisk 1Disk 4Disk 2Disk 3Mass (kg)2.412.383.171.57Polar moment of inertia (kg·m2)0.00820.00820.01460.0059Unbalance (10−5 kg·m)221.53Table 3
Parameters of the intermediate bearing.
NameRadius of inner ring (mm)Radius of outer ring (mm)Number of rollersRadial clearance (μm)Value9.3514.1596Table 4
Parameters of SFDs.
NameInner rotorOuter rotorSFD ISFD IISFD IIIRadius (mm)251835Length (mm)131314Radial clearance (mm)0.250.180.35Dynamic viscosity (10−2 Pa·s)1.0752Table5 presents the node information of the finite element model after meshing. The model contains 54 elements with 46 beam elements, 4 disk elements, and 4 supporting elements.Table 5
Node information of the finite element model.
Node numberAxial location (m)Bearing/diskElement numberOuter diameter (m)Inner diameter (m)10.110.018020.1320.018030.1630.018040.1940.018050.2250.018060.2547960.018070.2887970.018080.3227980.018090.35879Disk 190.0180100.38369100.0180110.40859110.0180120.43349120.0180130.44479Bearing 1130.0180140.49140.0220150.56150.0220160.65025160.0220170.75298170.0220180.82180.0220190.95844190.0220201.03200.0220211.06517Bearing 4210.0220221.087220.0220231.1230.0220241.14Disk 4240.0220251.17250.0220261.19260.0220271.21088270.0170281.21488Bearing 2280.0140291.22290.0140301.23038End of the inner rotor310.642300.0350.03320.655310.0350.03330.6793Bearing 3320.0350.03340.6865330.0380.03350.7117340.0380.03360.7369350.0380.03370.7621Disk 2360.0380.03380.80784370.0380.03390.85358380.0380.03400.89932390.0380.03410.94506400.0380.03420.96Disk 3410.0380.03430.98420.0380.03441430.070.03451.02440.060.03461.04450.060.03471.0643Bearing 4460.060.03481.0738Based on the dynamic differential equations of the rotor system under arbitrary maneuvering flight conditions in Section2.1, the motion differential equations of the inner and outer rotors can be written as(24)ML+MuLq˙L+CL+ΩGL+CuL+CbLq˙L+KL+KuL+KbLqL=FuL+FbL+fL,bearing+fL,sfd+FL,rub,(25)MH+MuHq˙H+CH+ΩGH+CuH+CbHq˙H+KH+KuH+KbHqH=FuH+FbH+fH,bearing+fH,sfd+FH,rub,where the superscripts “L” and “H” denote the inner and outer rotors.Denoting the node number of the inner rotor and outer linked to the intermediate bearing asξ and η, respectively. qξbearing and qηbearing represent the translational displacement of node ξ and node η. Then, the displacement boundary conditions are obtained according to equation (18a)–(18e), and the coupling relationship between the inner and outer rotors at the position of the intermediate bearing is given as(26)fξL,bearing=−fηH,bearing=fbearingqξbearing,qηbearing.Equations (24)–(26) are the motion differential equations of the dual rotor system.Linearize the dynamics model of the dual rotor system and solve the eigenvalues for the linearized equations of motion using the state space method, and the Campbell diagram is obtained as shown in Figure5. The radial stiffness of the intermediate bearing is set to 2.5 × 106 N/m. The first three orders of critical speed with the rotational speed of the inner rotor ω1 as excitation frequency are 233.4 rad/s, 403.3 rad/s, and 986.1 rad/s. While the first three orders of critical speed with ω2 as excitation frequency are 143 rad/s, 248.1 rad/s, and 610.4 rad/s. The natural frequency of the 1st torsional vibration is 111.7 Hz. Besides, the mode shapes corresponding to the critical speeds are also shown in Figure 5.Figure 5
Campbell diagram of the dual-rotor system.To investigate the effects on the rotor dynamic characteristics of maneuvering flight, the actions of level flight and hovering flight are taken as examples to study. Besides, the magnitude of maneuvering overload is an important parameter affecting the vibration response of the rotor system. In general, the maneuvering overload does not exceed 9 g due to the tolerance limit of the human body, where g represents the acceleration due to gravity. However, the future high maneuverability military drone is selected to study in this article, and the maneuvering overload is set to 10 g. The maneuvering flight parameters are listed in Table6, and all flight parameters not explicitly mentioned are set to 0.Table 6
Parameters of the level and hovering flight.
ActionsZ˙B (m/s)ωB,Y(rad/s)Overload (g)Level flight9800Hovering flight98110


## 3.2. Implicit Newmark-β Method

The differential equations of motion of the dual rotor system can be solved using the numerical integration method. The DOFs of the rotor system are divided into two parts according to whether they are subjected to nonlinear forces or not,q=qkqjT, where qk and qj represent the linear and nonlinear DOFs vector, respectively.An implicit numerical integral method combining the Newmark-β and Newton–Raphson methods are adopted in this work. The flow chart of the solution is shown in Figure 6. First, the vibration response of the rotor system at time tn are qkn, qjn, q˙kn, q˙jn, q¨kn, and q¨jn, assuming that the displacement response of nonlinear DOFs at tn+1 is qjn+1; then, the q˙jn+1 can be obtained by Newmark hypothesis. Substituting qjn+1 and q˙jn+1 into equations (18a)–(18e) and (19a)–(19e), the nonlinear forces can be solved. Then, the equations transform into linear equations, which can be solved by Newmark-β method. Through Newton–Raphson iterative calculation, the convergence solution at tn+1 can be obtained. Then save the vibration response at tn+1 and start the calculation of the next time interval. For the Newmark-β method, we set α=0.3 and β=0.6.Figure 6
Flow chart of the implicit newmark-β method.


## 3.3. Response Analysis

Compared with the outer rotor, the inner rotor is slender in configuration and thus suffers from a higher vibration level in the torsional direction, and the spectrum characteristics of disks 3 and 4 in the lateral direction are basically consistent. Thus, the vibration response characteristics of disk 4 are the primary focus in this work.We assume that rub-impact occurs at the position of disk 4. The lateral vibration amplitude of disk 4 without considering rub-impact faults under the level flight and hovering flight conditions is calculated to determine the relevant parameters in the rubbing response analysis, such as the calculating rotational speeds of the rotor and the initial radial clearance between the rotor and stator. Figure7 shows that the lateral vibration amplitude of disk 4 varies with the rotational speed under the conditions of level flight and hovering flight. As the speed of the inner rotor increases from 10 rad/s to 550 rad/s with interval of 10 rad/s, three peaks in vibration amplitude are observed corresponding to the rotational speeds of 140 rad/s, 250 rad/s, and 410 rad/s, respectively. Compared with the Campbell diagram of the dual rotor system (Figure 5), the error of the value of the critical speeds is due to the selection of the rotational speed interval for the vibration response calculation and the linearization for the rotor system when calculating the eigenvalues. The rotational speeds corresponding to the three peaks of the vibration response from small to large are the first-order–second-order critical speeds with excitation frequency ω2, and second-order critical speed with excitation frequency ω1, respectively. The peak of the first-order critical speed with excitation frequency ω1 is not observed because the rotational speed interval for the vibration response calculation and the first-order critical speed with excitation frequency ω1 is close to the second-order critical speed with excitation frequency ω2.Figure 7
Vibration response of disk 4 with respect to rotational speed. (a) Horizontal direction. (b) Vertical direction.
(a)(b)Besides, the action of hovering flight causes a static offset of the rotor in both the horizontal and vertical directions. In the horizontal direction, the rotor system is subjected to a constant centrifugal force generated by maneuvering flight, and the static offset does not change with the rotational speed. While an additional gyroscopic moment related to the rotational speed is generated in the vertical direction, which causes the static offset to vary with the rotational speed.The waterfall spectrum plots for the lateral vibration response of disk 4 during level flight and hovering flight conditions are shown in Figure8. Under the level flight condition, the frequency components of the vibration response mainly include the rotational speeds of the inner and outer rotors, i.e., ω1, ω2, and the combinational frequency components due to the coupling of the inner and outer rotors through the intermediate bearing, such as 2ω2−ω1, 3ω2−2ω1, and 4ω2−3ω1, and the like. The two vibration peaks with frequency ω2 correspond to the first- and second-order critical speeds with excitation frequency ω2 (140 rad/s and 250 rad/s), respectively. While the two vibration peaks with frequency ω1 correspond to the first- and second-order critical speeds with excitation frequency ω1 (240 rad/s and 410 rad/s). Additionally, a subharmonic vibration frequency component ωs1 induced by the bearing nonlinearity can be observed. Compared with the condition of level flight, the combinational frequency components of the inner and outer rotors are weakened, the amplitude of the subharmonic vibration is clearly enhanced, and a new subharmonic component ωs2 appears during hovering flight. The subharmonic vibration causes significant damage to the bearing of the rotor system in engineering and therefore requires serious consideration.Figure 8
Waterfall spectrum plots for the lateral vibration of disk 4 in vertical direction. (a) Level flight. (b) Hovering flight.
(a)(b)Figure9 shows the waterfall spectrum plots for the torsional vibration response of disk 4 during level flight and hovering flight. In this work, the amplitude of the torsional vibration is defined as the relative torsional angle of the two disks on the inner or outer rotor, that is, the torsional vibration amplitude of disk 4 is the difference between the torsional angles of disks 1 and 4. Under the level flight condition, the main frequency components of torsional vibration are ω2−ω1, 2ω2−2ω1, 3ω2−3ω1, and so on. However, significant changes take place in the frequency components during hovering flight; the subharmonic vibration component in the lateral vibration amplified by maneuvering flight is transmitted to the torsional direction due to the coupling effect of lateral and torsional vibrations, and the main frequency components change to ω1, ω2, ω2−ω1, the subharmonic vibration frequency ωs1, ωs2, and the combinational frequency ω1+ωs1.Figure 9
Waterfall spectrum plots for the torsional vibration of disk 4. (a) Level flight. (b) Hovering flight.
(a)(b)Assuming that the rub-impact occurs at the second-order critical speed with excitation frequencyω2 (250 rad/s). Thus, to analyze the rub-impact response, the inner rotor speeds of 200 rad/s–300 rad/s are examined in detail, and the rotational speed interval is set to 2 rad/s. The relevant parameters of the casing are listed in Table 7. The friction coefficient is 0.1, and the energy dissipation coefficient is 0.16 s/m. The initial radial clearance between the rotor and casing at the position of disk 4 is 0.5 × 10−4 m under the condition of level flight. According to the static offset of the rotor and casing due to maneuvering flight, the initial clearance during hovering is set to 1 × 10−4 m.Table 7
Parameters of the casing.
NameMass (kg)Elastic modulus (GPa)Poisson ratioDeformation stiffness (N/m)Value6001960.258 × 10−7Figure10 compares the lateral vibration amplitude of disk 4 with or without rub-impact under the two flight conditions. In the level flight condition, the rotor whirls around the central axis of the casing. The whirling range is limited by the casing, and rub-impact may occur at an arbitrary phase of the full cycle. Thus, the vibration amplitudes of disk 4 with rub-impact are significantly different from those without rub-impact (see Figure 10(a)). In the case of hovering flight, the additional centrifugal force and gyroscopic moment are generated on the rotor, while the casing is only subjected to centrifugal force. Besides, the deformation stiffness of the casing is different from that of the rotor in general. These two factors cause different static offsets of the rotor and casing. Thus, full-cycle rub-impact does not occur during hovering flight, but rub-impact will take place at local positions where the radial clearance of the rotor and casing is relatively small. According to the simulation results in Figure 10(b), the casing does not completely limit the whirling range of the rotor, and the vibration amplitude of disk 4 increases due to rub-impact near the critical speed of the system.Figure 10
Lateral vibration response of disk 4 with or without rub-impact. (a) Level flight. (b) Hovering flight.
(a)(b)The comparison of the torsional vibration amplitude of disk 4 varying with the rotational speed is shown in Figure11, when considering rub-impact or not under the two flight conditions. It can be observed that the torsional vibration amplitude increases near the critical speed of the system due to the tangential force produced by rub-impact and the lateral-torsional coupling effect, and the vibration amplitude of rub-impact is much larger than that without the rub-impact in the torsional direction.Figure 11
Torsional vibration response of disk 4 with or without rub-impact. (a) Level flight. (b) Hovering flight.
(a)(b)Figure12 shows the waterfall spectrum plots for the lateral vibration response of disk 4 during level flight. Comparing Figure 12(a) with Figure 12(b), the frequency components of the vibration response change significantly when rub-impact occurs. The frequency component 2ω1, and some combinational frequencies are produced by the rub-impact, such as 2ω1−ω2, 3ω2−4ω1, 2ω2−ω1, and 3ω2−2ω1, etc.Figure 12
Waterfall spectrum plots for the lateral vibration of disk 4 during level flight. (a) Without rub-impact. (b) Rub-impact.
(a)(b)The waterfall spectrum plots are shown in Figure13 for the lateral vibration response of disk 4 during hovering flight. Compared with the level flight condition, rub-impact produces only small changes in the frequency components of the rotor system during hovering flight. It has been analyzed that the rub-impact can only occur at the local position where the radial clearance is relatively small during hovering flight; thus, the frequency component of the rubbing response is relatively simple, and only an additional frequency component of ωr is observed, which varies with the rotational speed of the rotor between 0.75ω1 and 0.9ω1.Figure 13
Waterfall spectrum plots for the lateral vibration of disk 4 during hovering flight. (a) Without rub-impact. (b) Rub-impact.
(a)(b)Figures14 and 15 show the waterfall spectrum plots for the torsional vibration response of disk 4 with or without rub-impact under the two flight conditions. The torsional vibration amplitude of the rotor increases significantly in the presence to rub-impact, and the vibration frequency component is different in the two flight conditions. During level flight, the frequency components mainly include the zero frequency and ω2−ω1 when rub-impact does not occur. While some combinational frequencies are observed due to rub-impact, such as 2ω2−2ω1, 3ω2−3ω1, 4ω2−4ω1, and so on. Under the condition of hovering flight, the frequency components include the zero frequency, ω1, ω2, and ω2−ω1 without rub-impact. When rub-impact occurs, the frequency components change to the zero frequency, ω2−ω1, ω2, and 2ω2.Figure 14
Waterfall spectrum plots for the torsional vibration of disk 4 during level flight. (a) Without rub-impact. (b) Rub-impact.
(a)(b)Figure 15
Waterfall spectrum plots for the torsional vibration of disk 4 during hovering flight. (a) Without rub-impact. (b) Rub-impact.
(a)(b)Figures16 and 17 show the orbit and Poincaré maps of disk 4 under different rotational speeds near the occurrence of rubbing during level flight, and the orbit and Poincaré maps during hovering flight are shown in Figures 18 and 19 . Before and after the occurrence of rub-impact, the shapes of the orbits are similar under the two flight conditions (see Figures 16(a), 16(b), 16(f), 18(a), 18(b), and 18(f)); however, the orbits have clearly different shapes when rub-impact occurs (see Figures 16(c)–16(e) and 18(c)–18(e)). Besides, under both flight conditions, the motion of the rotor system is in the state of quasiperiodic before and after rub-impact due to the nonlinearity of the intermediate bearing (see Figures 17(a), 17(b), 17(f), 19(a), 19(b), and 19(f)). When rub-impact occurs, the motion of the rotor system tends to chaos with the increase in the degree of rub-impact (see Figures 17(c)–17(e) and 19(c)–19(e)).Figure 16
Orbit of disk 4 with rub-impact during level flight. (a) 220 rad/s. (b) 226 rad/s. (c) 236 rad/s. (d) 252 rad/s. (e) 256 rad/s. (f) 270 rad/s.
(a)(b)(c)(d)(e)(f)Figure 17
Poincaré maps of disk 4 with rub-impact during level flight. (a) 220 rad/s. (b) 226 rad/s. (c) 236 rad/s. (d) 252 rad/s. (e) 256 rad/s. (f) 270 rad/s.
(a)(b)(c)(d)(e)(f)Figure 18
Orbit of disk 4 with rub-impact during hovering flight. (a) 220 rad/s. (b) 226 rad/s. (c) 236 rad/s. (d) 252 rad/s. (e) 256 rad/s. (f) 270 rad/s.
(a)(b)(c)(d)(e)(f)Figure 19
Poincaré maps of disk 4 with rub-impact during hovering flight. (a) 220 rad/s. (b) 226 rad/s. (c) 236 rad/s. (d) 252 rad/s. (e) 256 rad/s. (f) 270 rad/s.
(a)(b)(c)(d)(e)(f)Figure20 shows the bifurcation diagrams for the lateral vibration response of disk 4 during level flight and hovering flight. The full-cycle rubbing clearly has a greater effect on the periodicity of the motion of the rotor under the level flight condition. In the range of rotation speeds considered herein, the motion of the rotor in the lateral direction moves from quasiperiodic to chaotic and back to quasiperiodic. During hovering flight, the motion of the rotor is relatively regular, except for the region of the critical speed of the system.Figure 20
Bifurcation diagrams for the lateral vibration of disk 4 with rub-impact. (a) Level flight. (b) Hovering flight.
(a)(b)The bifurcation diagrams for the torsional vibration response of disk 4 during level flight and hovering flight are shown in Figure21. Similar to the lateral vibration, the motion of the rotor is quasiperiodic in the torsional direction before and after rub-impact. When rub-impact occurs, the motion becomes chaotic. Compared with the condition of level flight, the motion exhibits better periodicity during hovering flight.Figure 21
Bifurcation diagrams for the torsional vibration of disk 4 with rub-impact. (a) Level flight. (b) Hovering flight.
(a)(b)


## 4. Conclusions

Considering the lateral-torsional coupling effect, a dynamic finite element model of a rub-impact rotor system under arbitrary maneuvering flight conditions is presented to study the rub-impact response characteristics of the rotor system during hovering flight. An implicit numerical integration method combining the Newmark-β and Newton–Raphson methods is employed to solve the nonlinear equations. The spectral characteristics, periodicity, and stability of the motion of the rotor are analyzed using spectrum diagrams, bifurcation diagrams, and Poincaré maps. From the results obtained in this study, it may be concluded that(1)
Without considering rub-impact faults, the action of hovering flight causes a static offset of the rotor system in both the horizontal and vertical directions. The combinational frequency components of the inner and outer rotors are weakened, and the subharmonic vibrations are amplified in both lateral and torsional vibrations due to maneuvering overload.(2)
The form of the rub-impact during level and hovering flight conditions is different. Under the condition of level flight, the whirling range of the rotor is limited by the engine casing, and the rub-impact may occur at an arbitrary phase of the full cycle. The lateral vibration amplitude of the rotor in the case of rub-impact is significantly smaller than that without rub-impact. During hovering flight, only local rub-impact occurs, and the lateral vibration amplitude of the rotor increases due to rub-impact. Under both flight conditions, the torsional vibration amplitude increases significantly due to the rub-impact tangential force and the lateral-torsional coupling effect when the rub-impact occurs.(3)
Under level and hovering flight conditions, the frequency components of both lateral and torsional vibrations of the rotor will change when rub-impact occurs. Compared with the condition of level flight, the change of the frequency components during hovering flight is relatively small. The frequency characteristics when rub-impact occurs can provide useful information for fault diagnosis of the rotor system.(4)
Under both the conditions of level flight and hovering flight, the motion of the rotor system is in the state of quasiperiodic before and after the occurrence of rub-impact and becomes chaotic when rub-impact occurs. However, full-cycle rubbing under the level flight condition has a greater effect on the periodicity of the motion of the rotor, and the periodicity of the rubbing rotor remains largely regular during hovering flight, except for the region of the critical speed of the system.



---

*Source: 4077556-2021-10-13.xml*