Multibody dynamic simulation and vibration transmission characteristics of dualrotor system for aeroengine with rubbing coupling faults
Pingping Ma^{1} , Jingyu Zhai^{2} , Hao Zhang^{3} , Qingkai Han^{4}
^{1, 2, 3, 4}Dalian University of Technology, Dalian, China
^{4}Corresponding author
Journal of Vibroengineering, Vol. 21, Issue 7, 2019, p. 18751887.
https://doi.org/10.21595/jve.2019.20206
Received 10 September 2018; received in revised form 21 February 2019; accepted 21 March 2019; published 15 November 2019
JVE Conferences
In this paper, a dualrotor system multibody dynamic model with rubbing coupling faults is established for practical aeroengine. In the model, the rubbing fault simulation method is introduced, the coupling effect between the internal and external rotor is considered. The numerical simulations of rotor vibration are accomplished by the utilization of multibody dynamic platform, where the simulation model consists of discs unbalances and local rubimpact between discs and casing shells. The timedomain responses, the frequency spectra and the shaftcenter trajectories of dualrotor with different unbalance and different rubbing positions are obtained. The vibration and its transmission characteristics of the inner and outer rotors are calculated. Finally, the simulation results are compared with the measured vibration of a dual rotor tester with rubbing fault. The simulation results are consistent with the measured results, which confirms the feasibility of the established model and the multibody dynamics simulation method in this paper. The application of multibody dynamics simulation method in aeroengine can deepen the understanding of the internal operation nature and laws of aeroengine, reduce the repetition of physical tests, greatly improve the efficiency and quality of development, accelerate the research and manufacture process.
 A multibody dynamic model of rotorbearingrotor coupling system is established for aeroengine, and a rubbing coupling faults simulation method is proposed.
 The influence of rubbing position and unbalance on the vibrations of dualrotor with rubbing fault is analyzed, and the vibration and its transmission characteristics of dualrotor are calculated.
 The simulation results are compared with the measured vibration of a dual rotor tester with rubbing fault.
Keywords: aeroengine dualrotor system, rubimpact, multibody dynamics model, vibration coupling and transfer.
1. Introduction
When it comes to aeroengine rotor system, rotor rubbing fault becomes one of the most common faults [1]. It not only brings about increase of the rotor stator gap, wearing, cracking or even breaking of a blade bearing, but may result in motion of rotor system instability, deformation of casing shells caused by vibration, which will seriously affect the efficiency of engine once exceeding the standard [2, 3]. Hence, the research on the influence of dualrotor system rotorstator rubbing fault to system response characteristics occupies an important role in the diagnosis and identification of aeroengine rubbing fault.
In recent years, the research on the dynamic characteristics of the dualrotor system with rubimpact fault has been conducted by relevant scholars. Qingkai Han [4] built a rigidflexible multibody model of the dualrotor system with local rubimpacts with MSC. Adams, based on rigidflexible multibody simulations, investigated the nonlinear dynamical characteristic. Guihuo Luo [5] established a nonlinear dynamic model of dualrotor system with highdimensional rubbing fault, and analyzed the dynamic response characteristics of counterrotating dualrotor system. The results show that the friction fault will make the axis orbit of the system become messy and except the two unbalance excitation frequency, the combination frequency also appeared in frequency Spectrum. Huiqun Yuan [6] established a multidisc double rotorstator coupling system dynamic model, taking the coupling of lowpressure rotor and highpressure rotor through the intermediate bearing into consideration together with highpressure rotor and casing shells, in which the fourthorder partial differential equations of motion of the continuum double rotor multidiscstator system were formed. Besides, he also concentrated on analyzing the effect of speed and intermediate bearing’s stiffness to the rub response, the results show that appropriate settings of intermediate bearing stiffness can reduce the amplitude of rubimpact fault. Junhong Zhang et al. [7] had implemented tests on misalignment and singlepoint rubimpact fault based on the experimental platform of aeroengine rotor with casing shells, the results show that the rubbing fault causes a large number of doubling and combination frequency. Yang Yang etc. [8] established a dualrotor system capable of describing the mechanical vibration resulting from two imbalances and fixed point rubbing, studied the dynamic characteristics of dualrotor system with fixed point rubbing by the means of theoretical analysis and experimental. The results indicate that for the case of fixed point rubbing, several combination frequencies in the responses of the dualrotor system are identified as the particular fault frequencies. Chen G. [9] established a new rotorball bearingssupportstator coupling system dynamic model and studied the rubbing coupling faults used the method of combining theory with experiment. Hui Ma, etc. [10] summarized the research results of rubimpact faults home and abroad and analyzed the typical fault characteristics for three kinds of rubimpact forms, such as singlepoint rubbing, local rubbing and circular rubbing.
Aero  engine structure is complex, working state is changeable, and it is difficult to directly reproduce vibration fault on the prototype machine. Multibody dynamic analysis technology based on virtual prototyping has been widely used in the research of rotor dynamics. ADAMS is based on the theory of multibody system dynamics and contains many professional modules. It can be used to establish complex kinematics and dynamics models of mechanical system and the simulation of system motion and dynamic characteristics under actual working conditions can be realized. Based on the theories of dynamics of dualrotor system, together with the rubimpact mechanism involving its stiffness and damping and friction, a multibody dynamic model of the aeroengine dualrotor system is built on the software platform of ADAMS. Based on the multibody dynamic model, the unbalance and rubimpact faults of dualrotor are simulated, and the vibration transmission law of unbalance and rubimpact caused by vibration coupling is revealed through the analysis of the results. It provides a verification method for the theoretical research of aeroengine rotor with rubbing faults.
2. Establishment of the dynamic model for dualrotor system
2.1. A simplified model of dualrotor system
Fig. 1 is a typical aeroengine dualrotor system structure, the lowpressure rotor consists of four fans and one turbine, supported by bearings of 1, 2, 3 and 6. The high pressure rotor is composed of nine compressor discs and a highpressure turbine, which is supported by bearings of 3 and 5.
Fig. 1. The supporting structure of dualrotor system
According to the principle of dynamic similarity, the dualrotor system of figure 1 is simplified. The rotating shaft is simplified as beam element. The discs and blades of the compressor and the turbine are considered as rigid bodies with rotating effect, which are simplified as discs in the form of concentrated mass and moment of inertia in their centroid position. In order to facilitate calculation, the coupling in the prototype machine is simplified as a rotating shaft and the supports are simplified as isotropic linear elastic supports. The specific simplification methods are as follows.
The fourstage fans in the lowpressure rotor is simplified as the disc 1 (referred to as LPC), the lowpressure turbine disc is equivalent to disc 2 (referred to as LPT). Nine compressor discs of highpressure rotor are simplified as disc 3 (referred to as HPC), highpressure turbine disc equivalent to disc 4 (referred to as HPT). The spring stiffness ${K}_{1}$ is the equivalent stiffness of bearing 1, squirrel cage, bearing cone and auxiliary plate. The stiffness ${K}_{2}$ is the combination stiffness of bearings 2 and 3. The stiffness of bearings 4, 5 and 6 are equivalent to ${K}_{3}$, ${K}_{4}$ and ${K}_{5}$ respectively. The dynamic model of Fig. 2 is obtained through simplification,
Fig. 2. The structure schematic diagram of dualrotor system
2.2. Equations of motions
The system analytic model is established by using the Lagrange equation as follows:
where, ${q}_{j}$ and ${\dot{q}}_{j}$ are generalized coordinates and generalized velocities respectively. $T$ and $U$ are the kinetic energy and potential energy of the system respectively. $D$ is the energy dissipation function of the system. ${Q}_{j}\left(t\right)$ is the generalized external excitation force. If the work done by some exciting force has been expressed as the kinetic energy and potential energy form of the vibration system, or the form of energy dissipation function, then these exciting forces are no longer considered on the right side of the equation. In this rotor system, the influence of energy dissipation function $D$ is neglected. Function Eq. (1) can be written as:
Considering the motion equation of the rotor, and ignoring the external excitation. The differential equations of the system motion are derived and arranged as follows:
${m}_{1}{\ddot{z}}_{1}+{k}_{cz}\left({\theta}_{y1}{c}_{1}{\theta}_{y2}{c}_{2}+{z}_{1}{z}_{2}\right)+{k}_{b1z}\left({z}_{1}+{a}_{1}{\theta}_{y1}\right)=0,$
${J}_{d1}{\ddot{\theta}}_{y1}+{J}_{p1}{\mathrm{\Omega}}_{1}{\dot{\theta}}_{z1}+{k}_{cz}{\theta}_{y1}{{c}_{1}}^{2}{k}_{cz}\left({\theta}_{y2}{c}_{2}+{z}_{1}{z}_{2}\right){c}_{1}$
$+\left({k}_{b1z}{{a}_{1}}^{2}+{k}_{c\theta y}\right){\theta}_{y1}+{k}_{b1z}{a}_{1}{z}_{1}{k}_{c\theta y}{\theta}_{y2}=0,$
${J}_{d1}{\ddot{\theta}}_{z1}{J}_{p1}{\mathrm{\Omega}}_{1}\dot{\theta}{y}_{1}+{k}_{cy}{{c}_{1}}^{2}{\theta}_{z1}+{k}_{cy}\left({\theta}_{z2}{c}_{2}+{y}_{1}{y}_{2}\right){c}_{1}$
$+\left({k}_{b1y}{{a}_{1}}^{2}+{k}_{c\theta z}\right){\theta}_{z1}{k}_{b1y}{a}_{1}{k}_{c\theta z}{\theta}_{z2}=0,$
$+{k}_{cy}\left({\theta}_{z1}{c}_{1}{y}_{1}+{y}_{2}\right)+\left({k}_{b2y}+{k}_{b5y}\right){y}_{2}=0,$
${m}_{2}{\ddot{z}}_{2}+\left({k}_{cz}{c}_{2}+{k}_{b2z}{a}_{2}+{k}_{b4z}{a}_{41}{k}_{b5z}{a}_{5}\right){\theta}_{y2}+{k}_{b4z}\left({\theta}_{y3}{a}_{4}+{z}_{2}{z}_{3}\right)$
$+{k}_{cz}\left({\theta}_{y1}{c}_{1}{z}_{1}+{z}_{2}\right)+\left({k}_{b2z}+{k}_{b5z}\right){z}_{2}=0,$
${J}_{d2}{\ddot{\theta}}_{y2}+{J}_{p2}{\Omega}_{1}{\dot{\theta}}_{z2}+\left({k}_{cz}{{c}_{2}}^{2}+{k}_{b2z}{{a}_{2}}^{2}+{k}_{b4z}{{a}_{41}}^{2}+{k}_{b5z}{{a}_{5}}^{2}+{k}_{c\theta y}+{q}_{b4\theta y}\right){\theta}_{y2}$
$+{k}_{b4z}\left({\theta}_{y3}{a}_{4}+{z}_{2}{z}_{3}\right){a}_{41}{k}_{cz}\left({\theta}_{y1}{c}_{1}+{z}_{1}{z}_{2}\right){c}_{2}$
$+\left({k}_{b2z}{a}_{2}{k}_{b5z}{a}_{5}\right){z}_{2}{k}_{c\theta y}{\theta}_{y1}{k}_{b4\theta y}{\theta}_{y3}=0,$
${J}_{d2}{\ddot{\theta}}_{z2}{J}_{p2}{\mathrm{\Omega}}_{1}{\ddot{\theta}}_{y2}+\left({k}_{cy}{{c}_{2}}^{2}+{k}_{b2y}{{a}_{2}}^{2}+{k}_{b4y}{{a}_{41}}^{2}+{k}_{b5y}{{a}_{5}}^{2}+{k}_{c\theta z}+{k}_{b4\theta z}\right){\theta}_{z2}$
$+{k}_{b4y}\left({\theta}_{z3}{a}_{4}{y}_{2}+{y}_{3}\right){a}_{41}+{k}_{cy}\left({\theta}_{z1}{c}_{1}+{y}_{1}{y}_{2}\right){c}_{2}$
$+\left({k}_{b2y}{a}_{2}+{k}_{b5y}{a}_{5}\right){y}_{2}{k}_{c\theta z}{\theta}_{z1}{k}_{b4\theta z}{\theta}_{z3}=0,$
$\left({m}_{3}+{m}_{4}\right){\ddot{y}}_{3}+{m}_{4}{a}_{34}{\ddot{\theta}}_{z3}+{k}_{b4y}\left({\theta}_{z2}{a}_{41}+{\theta}_{z3}{a}_{4}{y}_{2}+{y}_{3}\right)$
${k}_{b3y}{\theta}_{z3}{a}_{3}+{k}_{b3y}{y}_{3}=0,$
$\left({m}_{4}{{a}_{34}}^{2}+{J}_{d3}+{J}_{d4}\right){\ddot{\theta}}_{y3}{m}_{4}{a}_{34}{\ddot{z}}_{3}+\left({J}_{p3}+{J}_{p4}\right){\Omega}_{3}{\dot{\theta}}_{z3}+{{a}_{4}}^{2}{k}_{b4z}{\theta}_{y3}$
$+{k}_{b4z}\left({\theta}_{y2}{a}_{41}+{z}_{2}{z}_{3}\right){a}_{4}+\left({k}_{b3z}{{a}_{3}}^{2}+{k}_{b4\theta y}\right){\theta}_{y3}+{k}_{b3z}{a}_{3}{z}_{3}{k}_{b4\theta y}{\theta}_{y2}=0,$
$\left({m}_{4}{{a}_{34}}^{2}+{J}_{d3}+{J}_{d4}\right){\ddot{\theta}}_{z3}+{m}_{4}{a}_{34}{\ddot{y}}_{3}{\mathrm{\Omega}}_{3}\left({J}_{p3}+{J}_{p4}\right){\dot{\theta}}_{y3}+{k}_{b4y}{{a}_{4}}^{2}{\theta}_{z3}$
$+{k}_{b4y}\left({\theta}_{z2}{a}_{41}{y}_{2}+{y}_{3}\right){a}_{4}+\left({k}_{b3y}{{a}_{3}}^{2}+{k}_{b4\theta z}\right){\theta}_{z3}{k}_{b3y}{a}_{3}{y}_{3}{k}_{b4\theta z}{\theta}_{z2}=0.$
3. Establishment and verification of multibody dynamic model for dualrotor system with rubbing faults
3.1. Model verification
The basic parameters of the dualrotor structure are shown in Table 1. ${L}_{d}$ is the diameter of the lowpressure shaft, while ${H}_{1d}$ and ${H}_{2d}$ are respective illustration of the inner and outer diameter of highpressure shaft. In order to facilitate the calculation, the four discs use the same size, $d$ and $h$ are the diameter and thickness of the discs respectively.
Table 1. The basic parameters of dual rotor system
Parameters

Values

Parameters

Values

${L}_{AB}$

185 mm

${L}_{AI}$

1500 mm

${L}_{AC}$

330 mm

$d$

200 mm

${L}_{AD}$

500 mm

$h$

30 mm

${L}_{AE}$

685 mm

${L}_{d}$

65 mm

${L}_{AF}$

1015 mm

${H}_{1d}$

80 mm

${L}_{AG}$

1200 mm

${H}_{2d}$

90 mm

${L}_{AH}$

1315 mm

–

–

Table 2. Mode comparison
Based on finite element

Based on ADAMS

1st 118.2 Hz

1st 124.3 Hz

2st 171.7 Hz

2st 167.8 Hz

In order to verify the correctness of the model, the first two order natural frequencies of dualrotor based on ADAMS are compared with the results based on the finite element method.
The results show that the first two order modes are the same, and the frequency difference is less, which can verify the rationality of the model.
3.2. The multibody dynamic model of rubimpact
The ‘Impact’ is used to calculate the rubimpact force in ADAMS, which is based on the impact function to calculate the contact force between the two objects. The contact force is composed of two parts, the first part is the elastic force due to the mutual penetration between the two objects, the other part is the damping force produced by relative speed.
Rubimpact mechanics expression:
Among them, $k$ is the stiffness coefficient of rubimpact, ${\delta}_{0}$ and $x$ are the initial clearance and actual distance between the two contact objects in the rubbing process. ${C}_{max}$ is the maximum damping coefficient which used to characterize the collision energy loss. $N$ is the rotating shaft speed (inner and outer pressure rotor speed are ${N}_{2}$ and ${N}_{1}$ respectively). $e$ is rubimpact index which used to represent the nonlinear degree of material. $d$ is the penetration depth.
The multibody dynamic model of dualrotor rubimpact based on ADAMS is illustrated in Fig. 3. In the model, the elastic supportings are use the “spring” in the horizontal and vertical direction to simulate and the parameter setting of stiffness like the following Table 3. The dualrotor unbalance is simulated by adding a mass block on the lowpressure compressor disc and the highpressure compressor disc. A rubbing block is arranged at a certain distance from the lowpressure disc and the different rubbing conditions are simulated the material parameters are set according to steel. The dynamic friction between the turbine disc and the rubbing block is 0.1.
Fig. 3. The dualrotor system model based on Adams
Table 3. Parameter setting of the simulation model base on ADAMS
Parameter

Data

Parameter

Data

Rubimpact stiffness

1e6 N/m

Lowrotor motion ${N}_{1}$

8800 r/min

Rubimpact damping

100 Ns/m

Highrotor motion ${N}_{2}$

12000 r/min

Force exponent

1.5

${K}_{1}$

1e7 N/m

Penetration depth

0.1 mm

${K}_{2}$

5e7 N/m

Rubbing clearance

0.1 mm

${K}_{3}$

1e7 N/m

LPC unbalance quality ${m}_{e1}$

24.5 g

${K}_{4}$

1e8 N/m

HPC unbalance quality ${m}_{e3}$

20 g

${K}_{5}$

1e7 N/m

The software ADAMS provides a variety of integrators, this paper selects the GSTIFF integrator with high speed and high solve precision. There are three kinds of GSTIFF integrator integral format: I3, SI2 and SI1, the speed of solving I3 format is faster, but error will be generated when it comes to speed, especially the acceleration, the SI1 format is not suitable for processing contact problems with friction [11], so this paper chooses SI2 integral method, the integral error is set to 0.001.
4. Simulation results of dualrotor system
As what is shown in Table 4, the vibration characteristics of the dualrotor system under different operating conditions are analyzed. In this paper, timedomain responses, frequency spectra and shaftcenter trajectories of the high and low pressure discs are obtained via simulation analysis under four kinds of working conditions.
Table 4. Four kinds of working conditions
Parameter

Condition Ⅰ

Condition Ⅱ

Condition Ⅲ

Condition Ⅳ

${m}_{e1}$

√

√

√

√

${m}_{e3}$

√

–

–

√

LPTRubbing

–

√

–

–

HPTRubbing

–

–

√

√

4.1. Unbalance response of dualrotor system
The vibration of the dualrotor unbalance is simulated by applying me1 on low pressure compressor disc and ${m}_{e3}$ on high pressure compressor disc under working condition Ⅰ. The timedomain responses, frequency spectra and shaftcenter trajectories of the high and low pressure discs are obtained as shown in Fig. 4.
From the Fig. 4, we can find that the vibration of the dualrotor unbalance is in a certain vibration state. The four discs are in a certain range of periodic reciprocating vibration. The main frequencies of unbalanced response are ${N}_{1}$ and ${N}_{2}$. The time domain diagram and axis trajectory of the low and high pressure turbine discs are under zero bits by gravity.
The time domain diagram of low and high pressure compressor discs are regular, ${N}_{1}$ and ${N}_{2}$ occupies the absolute priority in the frequency domain respectively without other frequency components, and the amplitude of ${N}_{2}$ for the highpressure compressor disc reaches approximately half of the ${N}_{1}$ for lowpressure compressor disc. The trajectory line of the discs are in round shape. There are multiple peaks in time domain of the low and high turbine disc. ${N}_{1}$ and ${N}_{2}$ existing in the frequency domain and ${N}_{2}$ enjoys the absolute advantage. The shape of the discs trajectory appears some confusion within a certain range and the low pressure turbine disc more obvious than high pressure turbine disc.
From the Table 5, we can find that the amplitude of frequencies ${N}_{1}$ and ${N}_{2}$ for four discs decreases step by step along the transmission paths LPCLPTHPTHPC and HPCHPTLPTLPC respectively.
Table 5. Vibration amplitude of four discs
Parameter

LPC

HPC

HPT

LPT


Value (μm)

${N}_{1}$

${N}_{2}$

${N}_{1}$

${N}_{2}$

${N}_{1}$

${N}_{2}$

${N}_{1}$

${N}_{2}$

148.5

7.344

2.169

64.28

3.164

33.47

6.096

30.5

With the help of analytical method, the vibration of the four discs is measured during the unbalance of the LPC from 200 g.mm to 2000 g.mm. From Fig. 5(a) we can find that the amplitude of ${N}_{1}$ is in the order of LPC > LPT > HPT > HPC. The amplitude variation diagram of frequency ${N}_{2}$ of the four discs during the change of HPC unbalance from 200 g.mm to 2000 g.mm is obtained, as shown in Fig. 5(b). As can be seen from Fig. 5(b), the amplitude of ${N}_{2}$ is in the order of HPCHPTLPTLPC. Transmission law of unbalanced vibration of dualrotor based on analytical method is consistent with the simulation results based on multibody dynamics.
Fig. 4. Case Ⅰ, the vibration characteristic of the dualrotor system
a) LPC timedomain
b) LPC frequency spectra
c) LPC shaftcenter trajectories
d) HPC timedomain
e) HPC frequency spectra
f) HPC shaftcenter trajectories
g) HPT timedomain
h) HPT frequency spectra
i) HPT shaftcenter trajectories
j) LPT timedomain
k) LPT frequency spectra
l) LPT shaftcenter trajectories
Fig. 5. Variation curves of fundamental frequency of four discs vibration with unbalance
a)${N}_{1}$
b)${N}_{2}$
4.2. LPTrubbing
The vibration response characteristics of dualrotor with lowpressure turbine disc rubbing faults are simulated under the lowpressure compressor disc unbalance by using the operating condition Ⅱ. The timedomain responses, frequency spectra and shaftcenter trajectories of the high and low pressure discs are obtained as shown in Fig. 6.
Fig. 6. Case Ⅱ, the vibration characteristic of the dualrotor system
a) LPC timedomain
b) LPC frequency spectra
c) LPC shaftcenter trajectories
d) HPC timedomain
e) HPC frequency spectra
f) HPC shaftcenter trajectories
g) HPT timedomain
h) HPT frequency spectra
i) HPT shaftcenter trajectories
j) LPT timedomain
k) LPT frequency spectra
l) LPT shaftcenter trajectories
It can be seen from the figure that the frequency spectrum of each disc mainly falls to (2${N}_{2}{N}_{1}$)/2 and there are combination of frequency components such as 2${N}_{2}{N}_{1}$, ${N}_{2}{N}_{1}$, $\text{(4}{N}_{1}+{N}_{2}\text{)/2.}$ The amplitude of the principal frequency of the four discs is LPC > LPT > HPC > HPT.
The center orbit of the highpressure disc and lowpressure turbine discs rebounds along the rubbing direction. The trajectory lines of the lowpressure compressor disc and the turbine disc are circular, and disorder appears at the top of the trajectory of the lowpressure turbine disc. The trajectory lines of the highpressure compressor disc and turbine disc are circular and appear disorder at the bottom of the trajectory.
4.3. HPTrubbing
The vibration response of dualrotor and the transmission of vibration between high and low pressure rotor are analyzed under the situation of rotor unbalance and highpressure turbine rubbing. Simulating with the case Ⅲ and case Ⅳ in the Table 4. Their vibration features are shown in Fig. 7 and Fig. 8 respectively.
Fig. 7. Case Ⅲ, the vibration characteristic of the dualrotor system
a) LPC timedomain
b) LPC frequency spectra
c) LPC shaftcenter trajectories
d) HPC timedomain
e) HPC frequency spectra
f) HPC shaftcenter trajectories
g) HPT timedomain
h) HPT frequency spectra
i) HPT shaftcenter trajectories
j) LPT timedomain
k) LPT frequency spectra
l) LPT shaftcenter trajectories
Fig. 8. Case Ⅳ, the vibration characteristic of the dualrotor system
a) LPC timedomain
b) LPC frequency spectra
c) LPC shaftcenter trajectories
d) HPC timedomain
e) HPC frequency spectra
f) HPC shaftcenter trajectories
g) HPT timedomain
h) HPT frequency spectra
i) HPT shaftcenter trajectories
j) LPT timedomain
k) LPT frequency spectra
l) LPT shaftcenter trajectories
From the Fig. 7, we can find that the time domain diagram of lowpressure compressor disc is regular, the frequency spectrum is mainly dominated by the ${N}_{1}$ and the combination frequencies ${N}_{1}+{N}_{2}$, 1/4${N}_{2}$ and 2${N}_{2}{N}_{1}$ appeared. The trajectory line of LPC is in round shape. The time domain diagram of the highpressure compressor disc is under zero bits, the frequency spectrum are dominated by the combination frequencies ${N}_{2}$, ${N}_{1}$, ${N}_{1}+{N}_{2}$, 2${N}_{2}$ and the frequency ${N}_{2}$ occupies the absolute priority. The trajectory line of HPC appears some confusion within a certain range. The frequency spectrum of the highpressure turbine disc is dominated by the combination frequencies ${N}_{2}$ and 2${N}_{2}$, the trajectory line appears some confusion within a certain range. The frequency spectrum of the lowpressure turbine disc is mainly dominated by the frequencies ${N}_{2}$ and ${N}_{1}$, the combination frequencies ${N}_{1}+{N}_{2}$, 1/4${N}_{2}$ and 2${N}_{2}{N}_{1}$ appeared. The trajectory line of LPT appears complex like ‘petal’.
Compared the vibration response with LPTrubbing, the frequency spectrum of the four discs became more complex with HPT rubbing faults, the combination frequencies ${N}_{1}$, 2${N}_{2}{N}_{1}$, ${N}_{1}+{N}_{2}$, 2${N}_{2}$ and 1/4${N}_{2}$ occupied the priority. The changes of frequency spectrum amplitude of the highpressure turbine disc is most obvious, has significantly improved. The timedomain changes of lowpressure turbine disc is maximum, appeared many different peaks. The trajectory line of lowpressure turbine disc becomes the most complex like ‘petal’.
It can be seen from Fig. 8 that the time domain diagram of lowpressure compressor disc is regular, the frequency spectrum is dominated by ${N}_{1}$, together with the combination frequencies ${N}_{1}+{N}_{2}$, 1/4${N}_{2}$ and 2${N}_{2}{N}_{1}$, the trajectory line is in round shape. The frequency spectrum of highpressure compressor disc dominated by ${N}_{2}$, 2${N}_{2}$ and ${N}_{1}$, the trajectory line appears confusion within a small certain range. The frequency spectrum of the highpressure turbine disc and lowpressure turbine disc are dominated by ${N}_{2}$, 2${N}_{2}$ and ${N}_{1}$, in which ${N}_{2}$ take the absolute priority. Their trajectory lines appear some confusion within a small certain range and the trajectory lines of lowpressure turbine disc appeared more complex than the highpressure turbine disc.
Comparison Fig. 8 with Fig. 7, we can draw the conclusion: after adding the unbalance value of highpressure compressor disc, the frequency spectra of each disc are dominated by ${N}_{1}$ and ${N}_{2}$, the combination frequencies 2${N}_{2}{N}_{1}$, 2${N}_{2}$ and 1/4${N}_{2}$ appeared. The vibration characteristics of LPC not change. ${N}_{1}$ in the frequency spectrum of the highpressure compressor disc disappears, which is mainly composed of ${N}_{2}$ and 2${N}_{2}$. The amplitude of ${N}_{1}$, ${N}_{2}$ and 2${N}_{2}$ in the frequency spectrum of highpressure turbine disc are increased. The Lowpressure turbine disc’s frequency spectrum changes is the largest, the main frequency components from ${N}_{1}$ and ${N}_{2}$ in Fig. 6 into ${N}_{2}$ in Fig. 7, and the combination of frequency components reduced on the whole.
5. Verification of simulation results
Wen Bangchun, Jin Yezhuang and others have developed an aeroengine testing machine with similar structure, as shown in Fig. 9. Based on this testing machine, the rubbing fault of rotor is studied.
Fig. 9. Structure of aeroengine rubimpact rotor tester: 1 – motor, 2 – coupling, 3 – engine model, 4 – front horizontal support, 5 – eccentric transmission and locking mechanism, 6 – hanging, 710 – displacement sensor, 1112 – acceleration sensor
The time domain and frequency spectrum of vibration displacement in the vertical and horizontal directions of the dualrotor are obtained based on the dualrotor experiment rig. Through the analysis of the test results of dualrotor, it can be seen that when rubbing occurs in the rotor system, the frequency spectrum of the vibration signal contains not only the fundamental frequency components of the high and low pressure rotors, but also one or several combination frequencies components. Under the above different rubbing conditions, the test results show that the combination of the components of the signal in the vertical direction is shown in Table 6. It provides a verification method for simulation results.
Table 6. The frequency components of experimental and simulation results
Frequency components

Experimental results

Simulation results

${N}_{2}{N}_{1}$

√

√

2${N}_{1}{N}_{2}$

√

–

${N}_{1}$

√

√

${N}_{2}$

√

√

2${N}_{2}{N}_{1}$

√

√

2${N}_{2}$

√

√

The simulation result is also consistent with the measured vibration of an aeroengine happening rubimpact practically, which confirms the feasibility of the established model and simulation method.
6. Conclusions
In the dualrotor system of the aeroengine, the vibration of the inner and the outer rotor will be transmitted through intermediate bearing. The vibration of the inner rotor and the outer rotor will have coupling effect. The vibration caused by this coupling effect under rotor imbalance and rubbing fault has different transmission laws. The following conclusions are obtained through the simulation in this paper:
1) The unbalance response of the rotor is different due to the location and size of the unbalance. In the unbalance response of the dual rotor system, the unbalance response frequencies ${N}_{1}$ and ${N}_{2}$ of inner and outer rotors are absolutely dominant, and the axis trajectory is basically regular. The vibration of the dual rotor caused by unbalance of inner rotor and outer rotor decreases step by step on the transmission path LPCLPTHPTHPC and HPCHPTLPTLPC respectively.
2) In the rubimpact response of the dualrotor system, besides the inner and outer rotor rotation frequencies, the combined frequency of the inner and outer rotor rotation frequencies also appears. Due to the coupling effect of the intermediate bearings, the response of the system becomes complicated and is no longer the sum of the responses caused by simple unbalanced excitation.
3) In the vibration transmission of inner and outer rotors under the rubbing fault, the highpressure and lowpressure turbine disc are the most sensitive and the lowpressure compressor disc is the most stable. Vibration response of dualrotor system caused by highpressure rotor rubbing is the most complex, with many vibration frequency components and more resonance points, which is very unfavorable to the safe operation of the rotor. Therefore, in the dual rotor system, we should try our best to avoid the rubbing of the highpressure rotor or reduce the rubbing degree to reduce the harm caused by the highpressure rotor rubimpact.
In this paper, the research on unbalance fault and rubbing fault of $f$ dualrotor provides a verification means for further research on rubbing theory of aeroengine rotors.
References
 Xue S., Long X. U., Zhao Y., et al. Comparative analysis of Hertz contact simulation based on Adams and RecurDyn. Journal of Changchun University of Science and Technology, Vol. 39, Issue 4, 2016, p. 7377. [Search CrossRef]
 Chen Yuxu, Zhang Huabiao Research progress and Prospect of aircraft engine dynamics. Acta Aeronauticaet Et Astronautica Sinica, Vol. 32, Issue 8, 2011, p. 13711391. [Search CrossRef]
 Lahriri S., Weher H. I., Santos I. F., et al. Rotorstator contact dynamics using a nonideal driveTheoretical and experimental aspects. Journal of Sound and Vibration, Vol. 331, 2012, p. 45184536. [Publisher]
 Han Qingkai, Luo Haitao, Wen Bangchun Simulations of a dualrotor system with local rubimpacts based on rigidflexible multibody model. Key Engineering Materials, Vol. 413, Issue 414, 2009, p. 677682. [Publisher]
 Luo Guihuo, Yang Xiguan, Fei Wang Study on the rub impact response of high dimension double rotor system. Journal of Vibration Engineering, Vol. 28, Issue 1, 2015, p. 100107. [Search CrossRef]
 Yuan Huiqun, Wei He, Han Qinkai Analysis of rub impact fault of engine double rotor casing coupling system. Journal of Aerospace Power, Vol. 26, Issue 11, 2011, p. 24012408. [Search CrossRef]
 Zhang Junhong, Liang Ma, Ma Wenming, et al. Dynamic analysis and experimental study of a flexible rotor bearing system under the action of unbalance, rub and misalignment faults. Journal of Tianjin University (Science and Technology), Vol. 45, Issue 10, 2012, p. 885864. [Search CrossRef]
 Yang Yang, Cao Dengqing, Yu Tianhu, Wang Deyou, Li Chenggang Prediction of dynamic characteristics of a dualrotor system with fixed point rubbing – theoretical analysis and experimental study. International Journal of Mechanical Sciences, Vol. 115, Issue 116, 2016, p. 253261. [Publisher]
 Ma Hui, Yang Jian, Rongze Song, et al. Research progress and Prospect of rub impact fault of rotor system. Journal of Vibration and Shock, Vol. 33, Issue 6, 2014, p. 212. [Search CrossRef]
 Han Qinkai, Tao Yu, Wang Deyou, Tao Qu Nonlinear vibration analysis and diagnosis of Fault Rotor System. Science Press, Bei Jing, 2010. [Search CrossRef]
 Lee K. C., Hong D. K., Jeong Y. H., Kim C. Y. Dynamic simulation of radial active magnetic bearing system for high speed rotor using ADAMS and MATLAB cosimulation. Automation Science and Engineering, 2012, p. 880885. [Publisher]
 Jin Yezhuang, Wang Deyou, Wen Bangchun Rubbing characteristic test based on an aeroengine rotorstator testing rig. Nosie and Vibration Control, Vol. 3, Issue 35, 2015, p. 204207. [Search CrossRef]