# Research on the Models of Coupling Dynamics and Damage Classification for Vehicle-Engine Vibration

Authors: Ao Lei; Chuan-xue Song; Yu-long Lei; Yao Fu; Shang Zheng
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Many researchers have designed the dynamic models to study the vehicle-engine vibration. However, the existing mechanical models are relatively simple, and the analysis of engine vibration damage is discussed rarely. In this paper, we proposed the models of coupling dynamics and damage classification of vehicle-engine vibration. The key advantages of these proposed models are (1) the finite elements method is adopted for the rotor and casing system, and the complex structure with multirotor and multicasing is modeled by defining support system and linking methods; (2) the hybrid numerical integral method is used to obtain the inherent frequency of the nonlinear dynamic system; and (3) the algorithms based on backpropagation (BP) neural network and radial basis function (RBF) neural network are chosen to construct the damage classification model of rotors. Experimental results based on the engine rotor tester prove that the proposed models are not only more robust than the existing works but also show that the classification algorithms can support engine damage analysis effectively.


---

## Body

## 1. Introduction

Currently, the main methods of calculating the critical speed and imbalance steady state response of the rotor-support system are the conventional transfer matrix, finite element (FEM), and substructure modal synthesis [1]. Sun and Dai used the substructure transfer matrix to calculate and analyze the vibration of the engine [2]. Zheng et al. applied the overall transfer matrix to discuss the critical speed and strain energy distribution of a certain engine rotor-support-casing installation system [3]. In [4], the model synthesis method and FEM were combined to study the vibration characteristics of the engine rotor-support-casing system. The influence of rotor stiffness on the dynamic characteristics was analyzed [5], and the static stiffness, dynamic stiffness, and FEM model were used to calculate the characteristics. A flexible double-rotor FEM model with the nonlinearity of rolling bearing and squeeze film damper was proposed [6], and then the transient response simulation of blade loss was implemented. Sun et al. proposed a FEM model of two-rotor gas turbine engine [7]. In the work, the numerical integration method was used to calculate the sudden unbalance response generated by the blade loss. Meanwhile, the touch nonlinearity of the rolling bearing, the nonlinearity of the squeezed oil film force, and the thermal growth effect of bearing components during blade loss were also discussed. Hai et al. built two-rotor dynamic model with nonlinear squeeze film damper bearing [8]. They utilized NASTRAN FEM software to obtain the model parameters of the linear undamped system, and then the nonlinear numerical simulation analysis was performed by MATLAB software. Finally, the calculation speed was greatly improved compared with the traditional method. Finley et al. provided an analytical approach for expeditiously understanding and solving the vibration problems [9]. They utilized the proper data collection and analysis techniques, and then the true source of the vibration can be discovered. In [10], the work considered the free vibration of a Jeffcott rotor whose shaft has a strong nonlinear elastic property. Wang et al. focused on the dynamic responses of whole aeroengine with blade-casing rubbing. They proposed several methods [11–15] to improve energy efficiency and compact structure. These results can provide important theoretical and engineering references for the safe operation of dual-rotor system and the exact identification of coupling faults.However, the methods of modal synthesis, the FEM, and the impedance coupling are limited in practical applications due to the model complexity and the large computational workload. For example, the transfer matrix method can only handle simple boundary conditions and the calculation accuracy of the substructure transfer method is related to the model simplification error. In addition, the existing methods focus on the analysis of linear systems. Actually, there are many nonlinear factors, such as rolling bearing clearance, nonlinear contact force, variable flexibility (VC) vibration, and squeeze film damper. While the rotor produces rubbing, looseness, misalignment, and bearing failure, the rotor-support-casing coupling system will show strong nonlinearity, which can make the existing methods noneffective. In [16], a new rotor-rolling bearing-mechanical coupling dynamics model was established. Considering the squeeze film damping effect, the rolling bearing nonlinearity, and the rubbing fault, the rotor was seen as the equal section free Euler beam. Then, the model was analyzed by the modal truncation method, and the numerical response method was used to obtain the system response. Herisanu et al. analyzed the nonlinear dynamic behavior of a rotating electrical machine rotor-bearing system using the optimal auxiliary functions method [17]. The approximate analytical solution is in very good agreement with the numerical simulation results, which prove the reliability of this procedure.In order to solve the complicated structure of the engine, this paper describes a general model of rotor-support-casing coupling dynamics for engine vibration, which uses the finite element method to model the rotor and casing system firstly. Secondly, this work takes the support system as a lumped parameter model, and then the various support and connection types were defined for the complex rotor-support-casing system modeling. Thirdly, the system nonlinear response is obtained by the numerical value method based on the nonlinear factors. In order to verify the method, the experiments have been implemented on the rotor testers.Besides the analysis of model coupling dynamics, this paper also studied the rotor damage detection. There are several intelligent methods that have been used, such as the BP neural network method [18–21], the radical basis function (RBF) neural network method [22–24], genetic algorithm [25], the decision tree method [26], the k-nearest neighbour method [27], and the support vector machine [28]. Among the above methods, the BP method is a kind of multilayered feed-forward neural network, which includes the signal forward transfer and the error backward propagation. With abilities of the distributed storage, nonlinear mapping, high accuracy, and parallel computing, it is adaptive for the dynamic damage detection of the rotor. Compared to BP, RBF neural network is an efficient feed-forward network with the best approximation performance and global optimal characteristics that other forward networks do not have, and the structure is simple and the training speed is fast. According to the advantages of the two algorithms, BP and RBF are selected to construct the detection models. Numerical analysis indicates that both the two classification models are feasible to build damage classification model, and RBF performs better than BP.The paper is organized as follows. Section2 presents the proposed coupled dynamics model. In Section 3, the rotor damage detection model is introduced. Section 4 shows the experimental results. In Section 5, we conclude this paper.


## 2. Complex Rotor-Support-Casing Dynamic Model

### 2.1. Rotor Model

Figure1 shows the proposed FEM rotor dynamic model. In the model, the rotor includes a number of supports and turntables, which are divided into common beam elements by the finite element method. Meanwhile, the shear deformation, gyro moment, and moment of inertia of the rotor system are also considered. Any rotor is coupled with other rotors, casings, and supports by nonlinear forces and moment of force. In addition, the rotor node will also be subjected to external excitation.Figure 1
Finite element rotor dynamic model.The rotor contains distributed mass, elastic shaft, and stiff plate, which haveN nodes and M plates. For each shaft unit, E, I, G, μ, L, ρ, and A represent the elastic modulus, area moment of inertia, shear modulus, Poisson ratio, shaft length, shaft density, and shaft cross-sectional area, respectively. For disc Pi, mrpi, Jddi, and Jpdi are disc mass, equatorial moment of inertia, and pole moment of inertia; Fxi and Fyi are the force of the ith node; Mxi and Myi are the moment of force of the ith node of the rotor.In the system,x, y, and z are the fixed coordinate system. In the deformed state, the position of any section can be determined as follows: the position of the elastic center line can be obtained by the displacements of x(s, t) and y(s, t), the orientation of the cross-sectional is decided by wrapping (s, t) angles of x and y, and the cross section also rotates its center line.


#### 2.1.1. Rigid Plate Element Equation of Motion

The plate mass, moment of inertia of the equator, the moment of inertia of the pole, and rotational angular velocity of plate are set asmp, Jdd, Jpd, and ω, respectively. According to the Lagrangian equation, the equation of motion of rigid plate with fixed coordinate system can be set as(1)Mtd+Mrdq¨d−ωGdq˙d=Qd,where Qd is generalized external force vector; Mtd and Mrd are mass matrix and mass inertia matrix; Gd is gyromatrix; and qd = [x, y, ϕ, ψ] is generalized displacement vector. They can be described as follows:(2)Mtd=mp0000mp0000000000,Mrd=0000000000Jdd0000Jdd,Gd=00000000000−Jpd00Jpd0.


#### 2.1.2. Ordinary Beam Element Equation of Motion

Let the element modulus of elasticity, the shear modulus, the Poisson ratio, the inner diameter, the outer diameter, and the length beE, G, μ, d, D, and L; then, I can be calculated:(3)I=π64D4−d4.Cross-sectional areaA:(4)A=π4D2−d2.Effective shear areaAs:(5)As=A10/91+1.6D×d/D2+d2As=A7+6μ/61+μ1+20+12μ/7+6μD×d/D2+d22.Shear deformation coefficientϕs:(6)ϕs=12EIGAsL2.Each beam element has 2 nodes and 8 degrees of freedom. Each node has 4 degrees of freedom, i.e., the coordinate system ofx and y, and the angles around them. The cross-sectional displacement of the unit as a function of time is also a function of the position along the axis of the unit. The generalized displacement of the unit endpoint with time change is(7)qet=q1eq2eq3eq4eq5eq6eq7eq8eT.According to Lagrange equation, the equation of motion of the beam element relative to the fixed coordinate is obtained:(8)Mte+Mreq¨e+−ωGeq˙e+Kbe−Kceqe=Qe,where Qe is the generalized external force vector, Mte and Mre are mass matrix and mass inertia matrix, Ge is gyromatrix, Kbe is shearing stiffness matrix, and Kae is the unit tensile stiffness matrix. So, they can be calculated as follows:(9)Mte=ρL1+ϕs2Mt100Mt4Mt300−Mt50Mt1−Mt400Mt3Mt500−Mt4Mt200−Mt5Mt60Mt400Mt2Mt500Mt6Mt300Mt5Mt500−Mt40Mt3−Mt500Mt1Mt400Mt5Mt600Mt4Mt20−Mt500Mt6−Mt400Mt2,(10)Mt1=1315+710ϕs+13ϕs2,Mt2=1105+160ϕs+1120ϕs2L2,Mt3=970+310ϕs+16ϕs2,Mt4=11210+11120ϕs+124ϕs2L,Mt5=13420+340ϕs+124ϕs2L,Mt6=−1140+160ϕs+1120ϕs2L2,(11)Mre=ρL1+ϕs2rρL2Mr100Mr4−Mr100Mr40Mr1−Mr400−Mr1−Mr400−Mr4Mr200Mr4Mr30Mr400Mr2−Mr400Mr3−Mr100−Mr2Mr100−Mr40−Mr1Mr400Mr1Mr400−Mr4Mr300Mr4Mr20Mr400Mr3−Mr400Mr2,(12)Mr1=65,Mr2=215+16ϕs+13ϕs2L2,Mr3=−130−16ϕs+16ϕs2L2,Mr4=110−12ϕsL,rρ=IA,(13)Ge=ρ15Lrρ1+ϕs20−G1G200G1G20G100G2−G100G2−G200−G4G200G30−G2G400G2−G300G1−G200−G1−G20−G100−G2G100−G4−G200G3G200−G40−G2−G300G2G40.In equation (13),(14)G1=36,G2=3L−15LϕS,G3=L2+5L2ϕs−5L2ϕs2,G4=4L2+5L2ϕs+10L2ϕs2,Kbe=EIL3Kb100Kb4−Kb100Kb40Kb1−Kb400−Kb1−Kb400−Kb4Kb200Kb4Kb30Kb400Kb2−Kb400Kb3−Kb100−Kb2Kb100−Kb30−Kb1Kb400Kb1Kb400−Kb4Kb300Kb4Kb20Kb400Kb3−Kb400Kb2,Kb1=121+ϕs;Kb2=4+ϕs1+ϕsL2;Kb3=2−ϕs1+ϕsL2;Kb4=61+ϕsL,Kae=H30L1+ϕs2Ka100Ka4−Ka100Ka40Ka1−Ka400−Ka1−Ka400−Ka4Ka200Ka4Ka30Ka400Ka2−Ka400Ka3−Ka100−Ka4Ka100−Ka40−Ka1Ka400Ka1Ka400−Ka4Ka300Ka4Ka20Ka400Ka3−Ka400Ka2,Ka1=36+60ϕs+30ϕs2,Ka2=4L2+5L2ϕs+2.5L2ϕs,Ka3=−L2+5L2ϕs+2.5L2ϕs,Ka4=3L.


#### 2.1.3. Motion Equation of Rotor System

By assembling the equation of motion of the element, the motion equation of the rotor system can be obtained:(15)Msq¨s+Cs−ωGsq˙s+Ksqs=Qs,where Qs is the system generalized external force vector, Ms is the system mass matrix, Gs is the system gyromatrix, Ks is the system stiffness matrix, and Cs is the system damping matrix.


### 2.2. Casing Model

There are three different treatment methods available for the casing: beam elements, cone shell elements, and curved shell elements. Although the casing is a kind of shell structure, its vibration mode is in various forms, including a mode with a circumferential wave number of 0, 1, 2, …. However, when it is coupled with the rotor, it can only be one circumferential wave number. A bending moment is generated at a coupling point with the rotor and is coupled to the rotor. As for the vibration modes of other circumferential wave numbers, the torque generated at the coupling node with the rotor is self-balancing and is not coupled with the rotor bending. Therefore, in the rotor dynamics analysis, the vibration mode of the casing with a circumferential wave number of 1 is usually considered. At this time, the cross section of the casing is not deformed, and its axial direction is a curved mode. Accordingly, the current general processing method is to treat the casing as a beam unit, which is equivalent to the nonrotating shaft (beam unit structure) and also considers the shear effect and moment of inertia [1].Therefore, this paper treats the casing as a nonrotating beam, which is the same as the rotor model. The FEM can be used to obtain the motion differential equation of the casing:(16)Mcq¨c+Ccq˙c+Kcqc=Qc,where Qc is the generalized external force vector of the casing system, Mc is the stiffness matrix of the casing system, Kc is the stiffness matrix of the casing system, and Cc is the damping matrix of the casing system and is assumed to be proportional damping.


### 2.3. Discrete Support Model

In order to model the coupling system for multirotor and multimachine of the actual engine, it is necessary to fully consider the connection and support between the rotor and the casing. Therefore, a variety of support and connection for rotor and casing are, namely, the rotor-casing support connection, the intermediate support connection between the rotors and the coupling between the rotors, the support connection between the casings, the elastic support between the casing and the foundation. The combination of these supports and connections can be used to model the multirotor and multicasing of engine.


#### 2.3.1. Supporting Connection between Rotor and Casing

ForRCi(i = 1, 2, …, N) between the rotor and the casing, rolling bearings, squeeze film dampers, bearing housings, and other components are included. This research assumes that mwi is the quality of the outer ring of the rolling bearing, mbi is the mass of the bearing seat, kti is the elastic bearing stiffness between the outer ring of the rolling bearing and the bearing seat, cti is the damping coefficient between the outer ring of the rolling bearing and the bearing housing and it is nonlinear damping if there is a squeeze film damper, kfi and cfi are the support stiffness and damping between the casing and the bearing chock, Fyri and Fxri are the force acting on the support for the rotor, and Fyci and Fxci are the force acting on the support for the casing.In this paper, the outer ring of the bearing is fixed on the bearing chock and the inner ring is fixed on the rotating shaft. The displacement of them node of the rotor is set to xrm and yrm, where x = xm − xwi, y = yrm − ywi. According to [29, 30], the bearing force of the rotor acting on the i support can be obtained:(17)Fxri=∑j=1NCbxcosθj+ysinθj−r03/2Hxcosθj+ysinθj−r0cosθj,Fyri=∑j=1NCbxcosθj+ysinθj−r03/2Hxcosθj+ysinθj−r0sinθj,where Cb is the Hertz contact stiffness obtained by Hertzian elastic analysis between the inner and outer of ring and the ball. H(∗) is the Heaviside function. When the function variable is greater than 0, the function value is 1; otherwise it is 0. θj is the angular position of the first ball and θj=ωcoge×t+2π/Nbj−1,j=1,2,…,Nb,Nb is the number of balls. ωcoge is the rotation speed of the frame, R is the radius of the outer raceway, r is the radius of the inner raceway, and ωcoge=ω×r/R+r, where ω is the angular velocity of the rotary shaft. Therefore, the differential equation of motion of the outer ring is(18)mbix¨bi+ktixbi−xwi+ctix˙bi−x˙wi=Fxroi,mbiy¨bi+ktiybi−ywi+ctiy˙bi−y˙wi=Fyroi−mbig,where Fdxi and Fdyi are the damping forces.


#### 2.3.2. Intermediary Bearing Support between Rotors

For the intermediate bearing support between the rotors, RRMi(i = 1, 2, …, N) includes rolling bearings, bearing housings, and other components. Here, the outer ring mass of the rolling bearing is set to mwi, mbi is the mass of the bearing seat, kti and cti are the elastic bearing stiffness and damping coefficient between the outer ring of the bearing and the bearing chock, and kfi and cfi are the bearing stiffness and damping between the outer rotor and the bearing chock. We suppose Fyrii and Fxr1i are the force acting on the inner rotor and the support and Fyroi and Fxroi are the force acting on the outer rotor on the support. RRMi(i = 1, 2, …, N) is connected to the inner rotor and the outer rotor.The outer ring of the intermediate bearing is fixed on the outer rotor and rotates with the outer rotor, and the inner ring is fixed on the inner rotor and rotates with the inner rotor. The rotation speeds are set toωout and ωin. The calculation formula of the bearing force is the same as formula (8). Since both the inner and outer rings of the bearing rotate, the rotation frequency of the cage is calculated:(19)ωcage=ωin×r+ωout×RR+r.Set the displacement of the rolling bearing toxrm and yrm, where x = xm − xwi, y = yrm − ywi. According to formula (17), the bearing force of the rotor acting on the support can be obtained as Fxηi and Fyrii(i = 1, 2, …, N). Therefore, the differential motion equation of the outer ring of the rolling bearing is(20)mwix¨wi+ktixwi−xbi+ctix˙wi−x˙bi=Fxrii,mwiy¨wi+ktiywi−ybi+ctijwi−y˙bi=Fyrii−mwig,i=1,2,…,N.Furthermore, the displacement of the outer rotor is set asxron and yron and the displacement of the bearing chock for RRMi is xbi and ybi. Then, the force acting on the outer rotor is(21)Fxroi=kfixron−xbi+cfix˙ron−x˙bi,Fyroi=kfiyron−ybi+cfiy˙ron−y˙bi,i=1,2,…,N.Therefore, the differential equation of motion of the bearing chock for is(22)mbix¨bi+ktixbi−xwi+ctix˙bi−x˙wi=Fxroi,mbiy¨bi+ktiybi−ywi+ctiy˙bi−y˙wi=Fyroi−mbig,i=1,2,…,N.


#### 2.3.3. Rotor-Rotor Coupling Connection

For the coupling connection RRCK(k = 1, 2, …, N) between the rotors, let the ith node of the left rotor connect the jth node of the right rotor, krr be the radial stiffness of the coupling, kra be the angular stiffness, crr be the radial damping, and cra be the angular damping. Meanwhile, xrli, yrli, ϕrli, ψrli, ψrli are set to the displacement of the ith node of the left rotor and x˙rli,y˙rli,ϕ˙rli,ψ˙rli are the velocity; xrrj, yrrj, ϕrrj, ψrrj are set to the displacement of the jth node of the ith rotor and x˙mj,y˙rvj,ϕ˙mj,ψmy are the velocity. Then, the forces and moments acting on the left rotor node i and the right rotor node j, Frxi, Fryi, Mrxi, Mryi and Frxj, Fryi, Mrxj, Mryj, can be obtained:(23)Frxi=krrxryj−xr1i+crrx˙ryj−x˙r1i,Fryi=krryryj−yr1i+crry˙ry−y˙r1i,Mrxi=krαϕryj−ϕrit+craϕ˙ryj−ϕ˙r1i,Mryi=krαψrrj−ψr1i+craψ˙ryj−ψ˙r1i,Fny=−Fmi,Mmy=−Mmi,Fny=−Fyi,Mny=−Myi.


#### 2.3.4. Elastic Connection between the Casing

For the elastic connection between the casingsCCk(k = 1, 2, …, N), let the ith node of casing 1 connect the jth node of casing 2 by bolts, kcr be the radial stiffness of the bolt, kcα be the angular stiffness, ccr be the radial damping, and ccα be the angular damping. Meanwhile, let xc1i, yc11, ϕc1i, ψc1i be the displacement of the ith node of casing 1 and x˙c1i,y˙c1i,ψ˙c1i,ψ˙cli be the velocity; let xc2j, yc2j, ϕc2j, ψc2j be the displacement of the jth node of casing 2 and x˙c2j,y˙c2j,ϕ˙c2j,ψ˙c2j be the velocity. Then, the forces and moments acting on the node i of casing 1 are Fxi, Fyi, Mxi, Myi, and the forces and moments acting on the node j of the casing 2 are Fxj, Fyj, Mxy, Myi. The corresponding formulas can be obtained:(24)Fc1xi=kcrxc2j−xc1i+ccrx˙c2j−x˙c1i,Fc1yi=kcryc2j−yc1i+ccry˙c2j−y˙c1i,Mc1xi=kcαϕc2j−ϕc1i+ccαϕ˙c2j−ϕ˙c1i,Mc1yi=kcαψc2j−ψc1i+ccαψ˙c2j−ψ˙c1i,Fc2xj=−Fc1xi,Fc2yj=−Fc1yi,Mc2xj=−Mc1xi,Mc2yj=−Mc1yi.


### 2.4. Numerical Simulation Solution Method

Since the dynamic model of the complex rotor-support-case coupling system has many free numbers and a lot of nonlinear factors, the effective method is nonlinear iteration. In this paper, the differential equations are solved by the combination of the Newmark method and an improved method. As shown in Figure2, the whole process is described. The Newmark method is used to solve the formed matrix of the rotor and the finite element model of casing. The support members that do not need to form matrix are solved by hybrid integral method. The method is characterized by the fact that only the dynamic matrix of a single rotor or casing component needs to be assembled, without forming a bulky matrix of the entire system.Figure 2
Complex rotor-support-casing coupled dynamics solution process.Let the system’s dynamic equation be the following form:(25)MX¨+KX=P,where [M], [K] are system inertia and stiffness matrix; {X} is the generalized displacement vector of the system; X¨ is the generalized acceleration vector of the system; and {P} is the generalized load vector of the system.The steps of the proposed method are as follows:(i)
Initial calculation:(1)
Generate stiffness matrix [K] and mass matrix [M].(2)
ObtainX0,V0,A0.(3)
Choose step Δt and parameters α, β and then calculate constant:(26)a0=1αΔt2,α1=βαΔt,a2=1αΔt,a3=12α−1,a4=βα−1,a5=Δt2βα−2,a6=Δt1−β,a7=βΔt.(4)
Generate effective stiffness matrix:(27)K˜=K+a0M+a1C.(5)
Obtain inverse matrixK˜−1.(ii)
Calculate each time step:(1)
The load vector atn + 1:(28)P˜n+1=Pn+Ma0Xn+a2Vn+a3An+Ca1Xn+a4Vn+a5An.(2)
Get the displacement, velocity, and acceleration atn + 1:(29)Xn+1=K˜−1P˜n+1,An+1=a0Xn+1−Xn−a2Vn−a3An,Vn+1=Vn+a6An+a7An+1.For equation (25), the integral of hybrid method is(30)Xn+1=Xn+VnΔt+12+ψAnΔt2−ψAn−1Δt2,Vn+1=Vn+1+ϕAnΔt−ϕAn−1Δt,where Δt is the time step, n represents t = nΔt instant, n + 1 describes t = (n + 1)Δt instant, n − 1 is the instant of t = (n − 1)Δt, ψ, ϕ are instant parameters, and ψ = ϕ = 1/3. Take ϕ, for example, the reason for choosing parameters as follows:(31)An+1/2=32An−12An−1,An+1=2An−An−1,Vn+1=Vn+ΔtAn+1+4An+1/2+An/6=Vn+43AnΔt−13An−1Δt.For the two methods, we can find that the Newmark-β method needs to form the dynamic matrix, but does not require the diagonalization of the mass matrix. The improved hybrid integral method (31) does not require the formation of the dynamic matrix, and it can be solved directly from the differential equation. But the mass matrix is required to be the angular array, so this paper combines the implicit integral Newmark-β and explicit hybrid integral method to avoid the formation of large dynamic matrix, which could improve the efficiency of system modeling and the speed of solution.


## 2.1. Rotor Model

Figure1 shows the proposed FEM rotor dynamic model. In the model, the rotor includes a number of supports and turntables, which are divided into common beam elements by the finite element method. Meanwhile, the shear deformation, gyro moment, and moment of inertia of the rotor system are also considered. Any rotor is coupled with other rotors, casings, and supports by nonlinear forces and moment of force. In addition, the rotor node will also be subjected to external excitation.Figure 1
Finite element rotor dynamic model.The rotor contains distributed mass, elastic shaft, and stiff plate, which haveN nodes and M plates. For each shaft unit, E, I, G, μ, L, ρ, and A represent the elastic modulus, area moment of inertia, shear modulus, Poisson ratio, shaft length, shaft density, and shaft cross-sectional area, respectively. For disc Pi, mrpi, Jddi, and Jpdi are disc mass, equatorial moment of inertia, and pole moment of inertia; Fxi and Fyi are the force of the ith node; Mxi and Myi are the moment of force of the ith node of the rotor.In the system,x, y, and z are the fixed coordinate system. In the deformed state, the position of any section can be determined as follows: the position of the elastic center line can be obtained by the displacements of x(s, t) and y(s, t), the orientation of the cross-sectional is decided by wrapping (s, t) angles of x and y, and the cross section also rotates its center line.


### 2.1.1. Rigid Plate Element Equation of Motion

The plate mass, moment of inertia of the equator, the moment of inertia of the pole, and rotational angular velocity of plate are set asmp, Jdd, Jpd, and ω, respectively. According to the Lagrangian equation, the equation of motion of rigid plate with fixed coordinate system can be set as(1)Mtd+Mrdq¨d−ωGdq˙d=Qd,where Qd is generalized external force vector; Mtd and Mrd are mass matrix and mass inertia matrix; Gd is gyromatrix; and qd = [x, y, ϕ, ψ] is generalized displacement vector. They can be described as follows:(2)Mtd=mp0000mp0000000000,Mrd=0000000000Jdd0000Jdd,Gd=00000000000−Jpd00Jpd0.


### 2.1.2. Ordinary Beam Element Equation of Motion

Let the element modulus of elasticity, the shear modulus, the Poisson ratio, the inner diameter, the outer diameter, and the length beE, G, μ, d, D, and L; then, I can be calculated:(3)I=π64D4−d4.Cross-sectional areaA:(4)A=π4D2−d2.Effective shear areaAs:(5)As=A10/91+1.6D×d/D2+d2As=A7+6μ/61+μ1+20+12μ/7+6μD×d/D2+d22.Shear deformation coefficientϕs:(6)ϕs=12EIGAsL2.Each beam element has 2 nodes and 8 degrees of freedom. Each node has 4 degrees of freedom, i.e., the coordinate system ofx and y, and the angles around them. The cross-sectional displacement of the unit as a function of time is also a function of the position along the axis of the unit. The generalized displacement of the unit endpoint with time change is(7)qet=q1eq2eq3eq4eq5eq6eq7eq8eT.According to Lagrange equation, the equation of motion of the beam element relative to the fixed coordinate is obtained:(8)Mte+Mreq¨e+−ωGeq˙e+Kbe−Kceqe=Qe,where Qe is the generalized external force vector, Mte and Mre are mass matrix and mass inertia matrix, Ge is gyromatrix, Kbe is shearing stiffness matrix, and Kae is the unit tensile stiffness matrix. So, they can be calculated as follows:(9)Mte=ρL1+ϕs2Mt100Mt4Mt300−Mt50Mt1−Mt400Mt3Mt500−Mt4Mt200−Mt5Mt60Mt400Mt2Mt500Mt6Mt300Mt5Mt500−Mt40Mt3−Mt500Mt1Mt400Mt5Mt600Mt4Mt20−Mt500Mt6−Mt400Mt2,(10)Mt1=1315+710ϕs+13ϕs2,Mt2=1105+160ϕs+1120ϕs2L2,Mt3=970+310ϕs+16ϕs2,Mt4=11210+11120ϕs+124ϕs2L,Mt5=13420+340ϕs+124ϕs2L,Mt6=−1140+160ϕs+1120ϕs2L2,(11)Mre=ρL1+ϕs2rρL2Mr100Mr4−Mr100Mr40Mr1−Mr400−Mr1−Mr400−Mr4Mr200Mr4Mr30Mr400Mr2−Mr400Mr3−Mr100−Mr2Mr100−Mr40−Mr1Mr400Mr1Mr400−Mr4Mr300Mr4Mr20Mr400Mr3−Mr400Mr2,(12)Mr1=65,Mr2=215+16ϕs+13ϕs2L2,Mr3=−130−16ϕs+16ϕs2L2,Mr4=110−12ϕsL,rρ=IA,(13)Ge=ρ15Lrρ1+ϕs20−G1G200G1G20G100G2−G100G2−G200−G4G200G30−G2G400G2−G300G1−G200−G1−G20−G100−G2G100−G4−G200G3G200−G40−G2−G300G2G40.In equation (13),(14)G1=36,G2=3L−15LϕS,G3=L2+5L2ϕs−5L2ϕs2,G4=4L2+5L2ϕs+10L2ϕs2,Kbe=EIL3Kb100Kb4−Kb100Kb40Kb1−Kb400−Kb1−Kb400−Kb4Kb200Kb4Kb30Kb400Kb2−Kb400Kb3−Kb100−Kb2Kb100−Kb30−Kb1Kb400Kb1Kb400−Kb4Kb300Kb4Kb20Kb400Kb3−Kb400Kb2,Kb1=121+ϕs;Kb2=4+ϕs1+ϕsL2;Kb3=2−ϕs1+ϕsL2;Kb4=61+ϕsL,Kae=H30L1+ϕs2Ka100Ka4−Ka100Ka40Ka1−Ka400−Ka1−Ka400−Ka4Ka200Ka4Ka30Ka400Ka2−Ka400Ka3−Ka100−Ka4Ka100−Ka40−Ka1Ka400Ka1Ka400−Ka4Ka300Ka4Ka20Ka400Ka3−Ka400Ka2,Ka1=36+60ϕs+30ϕs2,Ka2=4L2+5L2ϕs+2.5L2ϕs,Ka3=−L2+5L2ϕs+2.5L2ϕs,Ka4=3L.


### 2.1.3. Motion Equation of Rotor System

By assembling the equation of motion of the element, the motion equation of the rotor system can be obtained:(15)Msq¨s+Cs−ωGsq˙s+Ksqs=Qs,where Qs is the system generalized external force vector, Ms is the system mass matrix, Gs is the system gyromatrix, Ks is the system stiffness matrix, and Cs is the system damping matrix.


## 2.1.1. Rigid Plate Element Equation of Motion

The plate mass, moment of inertia of the equator, the moment of inertia of the pole, and rotational angular velocity of plate are set asmp, Jdd, Jpd, and ω, respectively. According to the Lagrangian equation, the equation of motion of rigid plate with fixed coordinate system can be set as(1)Mtd+Mrdq¨d−ωGdq˙d=Qd,where Qd is generalized external force vector; Mtd and Mrd are mass matrix and mass inertia matrix; Gd is gyromatrix; and qd = [x, y, ϕ, ψ] is generalized displacement vector. They can be described as follows:(2)Mtd=mp0000mp0000000000,Mrd=0000000000Jdd0000Jdd,Gd=00000000000−Jpd00Jpd0.


## 2.1.2. Ordinary Beam Element Equation of Motion

Let the element modulus of elasticity, the shear modulus, the Poisson ratio, the inner diameter, the outer diameter, and the length beE, G, μ, d, D, and L; then, I can be calculated:(3)I=π64D4−d4.Cross-sectional areaA:(4)A=π4D2−d2.Effective shear areaAs:(5)As=A10/91+1.6D×d/D2+d2As=A7+6μ/61+μ1+20+12μ/7+6μD×d/D2+d22.Shear deformation coefficientϕs:(6)ϕs=12EIGAsL2.Each beam element has 2 nodes and 8 degrees of freedom. Each node has 4 degrees of freedom, i.e., the coordinate system ofx and y, and the angles around them. The cross-sectional displacement of the unit as a function of time is also a function of the position along the axis of the unit. The generalized displacement of the unit endpoint with time change is(7)qet=q1eq2eq3eq4eq5eq6eq7eq8eT.According to Lagrange equation, the equation of motion of the beam element relative to the fixed coordinate is obtained:(8)Mte+Mreq¨e+−ωGeq˙e+Kbe−Kceqe=Qe,where Qe is the generalized external force vector, Mte and Mre are mass matrix and mass inertia matrix, Ge is gyromatrix, Kbe is shearing stiffness matrix, and Kae is the unit tensile stiffness matrix. So, they can be calculated as follows:(9)Mte=ρL1+ϕs2Mt100Mt4Mt300−Mt50Mt1−Mt400Mt3Mt500−Mt4Mt200−Mt5Mt60Mt400Mt2Mt500Mt6Mt300Mt5Mt500−Mt40Mt3−Mt500Mt1Mt400Mt5Mt600Mt4Mt20−Mt500Mt6−Mt400Mt2,(10)Mt1=1315+710ϕs+13ϕs2,Mt2=1105+160ϕs+1120ϕs2L2,Mt3=970+310ϕs+16ϕs2,Mt4=11210+11120ϕs+124ϕs2L,Mt5=13420+340ϕs+124ϕs2L,Mt6=−1140+160ϕs+1120ϕs2L2,(11)Mre=ρL1+ϕs2rρL2Mr100Mr4−Mr100Mr40Mr1−Mr400−Mr1−Mr400−Mr4Mr200Mr4Mr30Mr400Mr2−Mr400Mr3−Mr100−Mr2Mr100−Mr40−Mr1Mr400Mr1Mr400−Mr4Mr300Mr4Mr20Mr400Mr3−Mr400Mr2,(12)Mr1=65,Mr2=215+16ϕs+13ϕs2L2,Mr3=−130−16ϕs+16ϕs2L2,Mr4=110−12ϕsL,rρ=IA,(13)Ge=ρ15Lrρ1+ϕs20−G1G200G1G20G100G2−G100G2−G200−G4G200G30−G2G400G2−G300G1−G200−G1−G20−G100−G2G100−G4−G200G3G200−G40−G2−G300G2G40.In equation (13),(14)G1=36,G2=3L−15LϕS,G3=L2+5L2ϕs−5L2ϕs2,G4=4L2+5L2ϕs+10L2ϕs2,Kbe=EIL3Kb100Kb4−Kb100Kb40Kb1−Kb400−Kb1−Kb400−Kb4Kb200Kb4Kb30Kb400Kb2−Kb400Kb3−Kb100−Kb2Kb100−Kb30−Kb1Kb400Kb1Kb400−Kb4Kb300Kb4Kb20Kb400Kb3−Kb400Kb2,Kb1=121+ϕs;Kb2=4+ϕs1+ϕsL2;Kb3=2−ϕs1+ϕsL2;Kb4=61+ϕsL,Kae=H30L1+ϕs2Ka100Ka4−Ka100Ka40Ka1−Ka400−Ka1−Ka400−Ka4Ka200Ka4Ka30Ka400Ka2−Ka400Ka3−Ka100−Ka4Ka100−Ka40−Ka1Ka400Ka1Ka400−Ka4Ka300Ka4Ka20Ka400Ka3−Ka400Ka2,Ka1=36+60ϕs+30ϕs2,Ka2=4L2+5L2ϕs+2.5L2ϕs,Ka3=−L2+5L2ϕs+2.5L2ϕs,Ka4=3L.


## 2.1.3. Motion Equation of Rotor System

By assembling the equation of motion of the element, the motion equation of the rotor system can be obtained:(15)Msq¨s+Cs−ωGsq˙s+Ksqs=Qs,where Qs is the system generalized external force vector, Ms is the system mass matrix, Gs is the system gyromatrix, Ks is the system stiffness matrix, and Cs is the system damping matrix.


## 2.2. Casing Model

There are three different treatment methods available for the casing: beam elements, cone shell elements, and curved shell elements. Although the casing is a kind of shell structure, its vibration mode is in various forms, including a mode with a circumferential wave number of 0, 1, 2, …. However, when it is coupled with the rotor, it can only be one circumferential wave number. A bending moment is generated at a coupling point with the rotor and is coupled to the rotor. As for the vibration modes of other circumferential wave numbers, the torque generated at the coupling node with the rotor is self-balancing and is not coupled with the rotor bending. Therefore, in the rotor dynamics analysis, the vibration mode of the casing with a circumferential wave number of 1 is usually considered. At this time, the cross section of the casing is not deformed, and its axial direction is a curved mode. Accordingly, the current general processing method is to treat the casing as a beam unit, which is equivalent to the nonrotating shaft (beam unit structure) and also considers the shear effect and moment of inertia [1].Therefore, this paper treats the casing as a nonrotating beam, which is the same as the rotor model. The FEM can be used to obtain the motion differential equation of the casing:(16)Mcq¨c+Ccq˙c+Kcqc=Qc,where Qc is the generalized external force vector of the casing system, Mc is the stiffness matrix of the casing system, Kc is the stiffness matrix of the casing system, and Cc is the damping matrix of the casing system and is assumed to be proportional damping.


## 2.3. Discrete Support Model

In order to model the coupling system for multirotor and multimachine of the actual engine, it is necessary to fully consider the connection and support between the rotor and the casing. Therefore, a variety of support and connection for rotor and casing are, namely, the rotor-casing support connection, the intermediate support connection between the rotors and the coupling between the rotors, the support connection between the casings, the elastic support between the casing and the foundation. The combination of these supports and connections can be used to model the multirotor and multicasing of engine.


### 2.3.1. Supporting Connection between Rotor and Casing

ForRCi(i = 1, 2, …, N) between the rotor and the casing, rolling bearings, squeeze film dampers, bearing housings, and other components are included. This research assumes that mwi is the quality of the outer ring of the rolling bearing, mbi is the mass of the bearing seat, kti is the elastic bearing stiffness between the outer ring of the rolling bearing and the bearing seat, cti is the damping coefficient between the outer ring of the rolling bearing and the bearing housing and it is nonlinear damping if there is a squeeze film damper, kfi and cfi are the support stiffness and damping between the casing and the bearing chock, Fyri and Fxri are the force acting on the support for the rotor, and Fyci and Fxci are the force acting on the support for the casing.In this paper, the outer ring of the bearing is fixed on the bearing chock and the inner ring is fixed on the rotating shaft. The displacement of them node of the rotor is set to xrm and yrm, where x = xm − xwi, y = yrm − ywi. According to [29, 30], the bearing force of the rotor acting on the i support can be obtained:(17)Fxri=∑j=1NCbxcosθj+ysinθj−r03/2Hxcosθj+ysinθj−r0cosθj,Fyri=∑j=1NCbxcosθj+ysinθj−r03/2Hxcosθj+ysinθj−r0sinθj,where Cb is the Hertz contact stiffness obtained by Hertzian elastic analysis between the inner and outer of ring and the ball. H(∗) is the Heaviside function. When the function variable is greater than 0, the function value is 1; otherwise it is 0. θj is the angular position of the first ball and θj=ωcoge×t+2π/Nbj−1,j=1,2,…,Nb,Nb is the number of balls. ωcoge is the rotation speed of the frame, R is the radius of the outer raceway, r is the radius of the inner raceway, and ωcoge=ω×r/R+r, where ω is the angular velocity of the rotary shaft. Therefore, the differential equation of motion of the outer ring is(18)mbix¨bi+ktixbi−xwi+ctix˙bi−x˙wi=Fxroi,mbiy¨bi+ktiybi−ywi+ctiy˙bi−y˙wi=Fyroi−mbig,where Fdxi and Fdyi are the damping forces.


### 2.3.2. Intermediary Bearing Support between Rotors

For the intermediate bearing support between the rotors, RRMi(i = 1, 2, …, N) includes rolling bearings, bearing housings, and other components. Here, the outer ring mass of the rolling bearing is set to mwi, mbi is the mass of the bearing seat, kti and cti are the elastic bearing stiffness and damping coefficient between the outer ring of the bearing and the bearing chock, and kfi and cfi are the bearing stiffness and damping between the outer rotor and the bearing chock. We suppose Fyrii and Fxr1i are the force acting on the inner rotor and the support and Fyroi and Fxroi are the force acting on the outer rotor on the support. RRMi(i = 1, 2, …, N) is connected to the inner rotor and the outer rotor.The outer ring of the intermediate bearing is fixed on the outer rotor and rotates with the outer rotor, and the inner ring is fixed on the inner rotor and rotates with the inner rotor. The rotation speeds are set toωout and ωin. The calculation formula of the bearing force is the same as formula (8). Since both the inner and outer rings of the bearing rotate, the rotation frequency of the cage is calculated:(19)ωcage=ωin×r+ωout×RR+r.Set the displacement of the rolling bearing toxrm and yrm, where x = xm − xwi, y = yrm − ywi. According to formula (17), the bearing force of the rotor acting on the support can be obtained as Fxηi and Fyrii(i = 1, 2, …, N). Therefore, the differential motion equation of the outer ring of the rolling bearing is(20)mwix¨wi+ktixwi−xbi+ctix˙wi−x˙bi=Fxrii,mwiy¨wi+ktiywi−ybi+ctijwi−y˙bi=Fyrii−mwig,i=1,2,…,N.Furthermore, the displacement of the outer rotor is set asxron and yron and the displacement of the bearing chock for RRMi is xbi and ybi. Then, the force acting on the outer rotor is(21)Fxroi=kfixron−xbi+cfix˙ron−x˙bi,Fyroi=kfiyron−ybi+cfiy˙ron−y˙bi,i=1,2,…,N.Therefore, the differential equation of motion of the bearing chock for is(22)mbix¨bi+ktixbi−xwi+ctix˙bi−x˙wi=Fxroi,mbiy¨bi+ktiybi−ywi+ctiy˙bi−y˙wi=Fyroi−mbig,i=1,2,…,N.


### 2.3.3. Rotor-Rotor Coupling Connection

For the coupling connection RRCK(k = 1, 2, …, N) between the rotors, let the ith node of the left rotor connect the jth node of the right rotor, krr be the radial stiffness of the coupling, kra be the angular stiffness, crr be the radial damping, and cra be the angular damping. Meanwhile, xrli, yrli, ϕrli, ψrli, ψrli are set to the displacement of the ith node of the left rotor and x˙rli,y˙rli,ϕ˙rli,ψ˙rli are the velocity; xrrj, yrrj, ϕrrj, ψrrj are set to the displacement of the jth node of the ith rotor and x˙mj,y˙rvj,ϕ˙mj,ψmy are the velocity. Then, the forces and moments acting on the left rotor node i and the right rotor node j, Frxi, Fryi, Mrxi, Mryi and Frxj, Fryi, Mrxj, Mryj, can be obtained:(23)Frxi=krrxryj−xr1i+crrx˙ryj−x˙r1i,Fryi=krryryj−yr1i+crry˙ry−y˙r1i,Mrxi=krαϕryj−ϕrit+craϕ˙ryj−ϕ˙r1i,Mryi=krαψrrj−ψr1i+craψ˙ryj−ψ˙r1i,Fny=−Fmi,Mmy=−Mmi,Fny=−Fyi,Mny=−Myi.


### 2.3.4. Elastic Connection between the Casing

For the elastic connection between the casingsCCk(k = 1, 2, …, N), let the ith node of casing 1 connect the jth node of casing 2 by bolts, kcr be the radial stiffness of the bolt, kcα be the angular stiffness, ccr be the radial damping, and ccα be the angular damping. Meanwhile, let xc1i, yc11, ϕc1i, ψc1i be the displacement of the ith node of casing 1 and x˙c1i,y˙c1i,ψ˙c1i,ψ˙cli be the velocity; let xc2j, yc2j, ϕc2j, ψc2j be the displacement of the jth node of casing 2 and x˙c2j,y˙c2j,ϕ˙c2j,ψ˙c2j be the velocity. Then, the forces and moments acting on the node i of casing 1 are Fxi, Fyi, Mxi, Myi, and the forces and moments acting on the node j of the casing 2 are Fxj, Fyj, Mxy, Myi. The corresponding formulas can be obtained:(24)Fc1xi=kcrxc2j−xc1i+ccrx˙c2j−x˙c1i,Fc1yi=kcryc2j−yc1i+ccry˙c2j−y˙c1i,Mc1xi=kcαϕc2j−ϕc1i+ccαϕ˙c2j−ϕ˙c1i,Mc1yi=kcαψc2j−ψc1i+ccαψ˙c2j−ψ˙c1i,Fc2xj=−Fc1xi,Fc2yj=−Fc1yi,Mc2xj=−Mc1xi,Mc2yj=−Mc1yi.


## 2.3.1. Supporting Connection between Rotor and Casing

ForRCi(i = 1, 2, …, N) between the rotor and the casing, rolling bearings, squeeze film dampers, bearing housings, and other components are included. This research assumes that mwi is the quality of the outer ring of the rolling bearing, mbi is the mass of the bearing seat, kti is the elastic bearing stiffness between the outer ring of the rolling bearing and the bearing seat, cti is the damping coefficient between the outer ring of the rolling bearing and the bearing housing and it is nonlinear damping if there is a squeeze film damper, kfi and cfi are the support stiffness and damping between the casing and the bearing chock, Fyri and Fxri are the force acting on the support for the rotor, and Fyci and Fxci are the force acting on the support for the casing.In this paper, the outer ring of the bearing is fixed on the bearing chock and the inner ring is fixed on the rotating shaft. The displacement of them node of the rotor is set to xrm and yrm, where x = xm − xwi, y = yrm − ywi. According to [29, 30], the bearing force of the rotor acting on the i support can be obtained:(17)Fxri=∑j=1NCbxcosθj+ysinθj−r03/2Hxcosθj+ysinθj−r0cosθj,Fyri=∑j=1NCbxcosθj+ysinθj−r03/2Hxcosθj+ysinθj−r0sinθj,where Cb is the Hertz contact stiffness obtained by Hertzian elastic analysis between the inner and outer of ring and the ball. H(∗) is the Heaviside function. When the function variable is greater than 0, the function value is 1; otherwise it is 0. θj is the angular position of the first ball and θj=ωcoge×t+2π/Nbj−1,j=1,2,…,Nb,Nb is the number of balls. ωcoge is the rotation speed of the frame, R is the radius of the outer raceway, r is the radius of the inner raceway, and ωcoge=ω×r/R+r, where ω is the angular velocity of the rotary shaft. Therefore, the differential equation of motion of the outer ring is(18)mbix¨bi+ktixbi−xwi+ctix˙bi−x˙wi=Fxroi,mbiy¨bi+ktiybi−ywi+ctiy˙bi−y˙wi=Fyroi−mbig,where Fdxi and Fdyi are the damping forces.


## 2.3.2. Intermediary Bearing Support between Rotors

For the intermediate bearing support between the rotors, RRMi(i = 1, 2, …, N) includes rolling bearings, bearing housings, and other components. Here, the outer ring mass of the rolling bearing is set to mwi, mbi is the mass of the bearing seat, kti and cti are the elastic bearing stiffness and damping coefficient between the outer ring of the bearing and the bearing chock, and kfi and cfi are the bearing stiffness and damping between the outer rotor and the bearing chock. We suppose Fyrii and Fxr1i are the force acting on the inner rotor and the support and Fyroi and Fxroi are the force acting on the outer rotor on the support. RRMi(i = 1, 2, …, N) is connected to the inner rotor and the outer rotor.The outer ring of the intermediate bearing is fixed on the outer rotor and rotates with the outer rotor, and the inner ring is fixed on the inner rotor and rotates with the inner rotor. The rotation speeds are set toωout and ωin. The calculation formula of the bearing force is the same as formula (8). Since both the inner and outer rings of the bearing rotate, the rotation frequency of the cage is calculated:(19)ωcage=ωin×r+ωout×RR+r.Set the displacement of the rolling bearing toxrm and yrm, where x = xm − xwi, y = yrm − ywi. According to formula (17), the bearing force of the rotor acting on the support can be obtained as Fxηi and Fyrii(i = 1, 2, …, N). Therefore, the differential motion equation of the outer ring of the rolling bearing is(20)mwix¨wi+ktixwi−xbi+ctix˙wi−x˙bi=Fxrii,mwiy¨wi+ktiywi−ybi+ctijwi−y˙bi=Fyrii−mwig,i=1,2,…,N.Furthermore, the displacement of the outer rotor is set asxron and yron and the displacement of the bearing chock for RRMi is xbi and ybi. Then, the force acting on the outer rotor is(21)Fxroi=kfixron−xbi+cfix˙ron−x˙bi,Fyroi=kfiyron−ybi+cfiy˙ron−y˙bi,i=1,2,…,N.Therefore, the differential equation of motion of the bearing chock for is(22)mbix¨bi+ktixbi−xwi+ctix˙bi−x˙wi=Fxroi,mbiy¨bi+ktiybi−ywi+ctiy˙bi−y˙wi=Fyroi−mbig,i=1,2,…,N.


## 2.3.3. Rotor-Rotor Coupling Connection

For the coupling connection RRCK(k = 1, 2, …, N) between the rotors, let the ith node of the left rotor connect the jth node of the right rotor, krr be the radial stiffness of the coupling, kra be the angular stiffness, crr be the radial damping, and cra be the angular damping. Meanwhile, xrli, yrli, ϕrli, ψrli, ψrli are set to the displacement of the ith node of the left rotor and x˙rli,y˙rli,ϕ˙rli,ψ˙rli are the velocity; xrrj, yrrj, ϕrrj, ψrrj are set to the displacement of the jth node of the ith rotor and x˙mj,y˙rvj,ϕ˙mj,ψmy are the velocity. Then, the forces and moments acting on the left rotor node i and the right rotor node j, Frxi, Fryi, Mrxi, Mryi and Frxj, Fryi, Mrxj, Mryj, can be obtained:(23)Frxi=krrxryj−xr1i+crrx˙ryj−x˙r1i,Fryi=krryryj−yr1i+crry˙ry−y˙r1i,Mrxi=krαϕryj−ϕrit+craϕ˙ryj−ϕ˙r1i,Mryi=krαψrrj−ψr1i+craψ˙ryj−ψ˙r1i,Fny=−Fmi,Mmy=−Mmi,Fny=−Fyi,Mny=−Myi.


## 2.3.4. Elastic Connection between the Casing

For the elastic connection between the casingsCCk(k = 1, 2, …, N), let the ith node of casing 1 connect the jth node of casing 2 by bolts, kcr be the radial stiffness of the bolt, kcα be the angular stiffness, ccr be the radial damping, and ccα be the angular damping. Meanwhile, let xc1i, yc11, ϕc1i, ψc1i be the displacement of the ith node of casing 1 and x˙c1i,y˙c1i,ψ˙c1i,ψ˙cli be the velocity; let xc2j, yc2j, ϕc2j, ψc2j be the displacement of the jth node of casing 2 and x˙c2j,y˙c2j,ϕ˙c2j,ψ˙c2j be the velocity. Then, the forces and moments acting on the node i of casing 1 are Fxi, Fyi, Mxi, Myi, and the forces and moments acting on the node j of the casing 2 are Fxj, Fyj, Mxy, Myi. The corresponding formulas can be obtained:(24)Fc1xi=kcrxc2j−xc1i+ccrx˙c2j−x˙c1i,Fc1yi=kcryc2j−yc1i+ccry˙c2j−y˙c1i,Mc1xi=kcαϕc2j−ϕc1i+ccαϕ˙c2j−ϕ˙c1i,Mc1yi=kcαψc2j−ψc1i+ccαψ˙c2j−ψ˙c1i,Fc2xj=−Fc1xi,Fc2yj=−Fc1yi,Mc2xj=−Mc1xi,Mc2yj=−Mc1yi.


## 2.4. Numerical Simulation Solution Method

Since the dynamic model of the complex rotor-support-case coupling system has many free numbers and a lot of nonlinear factors, the effective method is nonlinear iteration. In this paper, the differential equations are solved by the combination of the Newmark method and an improved method. As shown in Figure2, the whole process is described. The Newmark method is used to solve the formed matrix of the rotor and the finite element model of casing. The support members that do not need to form matrix are solved by hybrid integral method. The method is characterized by the fact that only the dynamic matrix of a single rotor or casing component needs to be assembled, without forming a bulky matrix of the entire system.Figure 2
Complex rotor-support-casing coupled dynamics solution process.Let the system’s dynamic equation be the following form:(25)MX¨+KX=P,where [M], [K] are system inertia and stiffness matrix; {X} is the generalized displacement vector of the system; X¨ is the generalized acceleration vector of the system; and {P} is the generalized load vector of the system.The steps of the proposed method are as follows:(i)
Initial calculation:(1)
Generate stiffness matrix [K] and mass matrix [M].(2)
ObtainX0,V0,A0.(3)
Choose step Δt and parameters α, β and then calculate constant:(26)a0=1αΔt2,α1=βαΔt,a2=1αΔt,a3=12α−1,a4=βα−1,a5=Δt2βα−2,a6=Δt1−β,a7=βΔt.(4)
Generate effective stiffness matrix:(27)K˜=K+a0M+a1C.(5)
Obtain inverse matrixK˜−1.(ii)
Calculate each time step:(1)
The load vector atn + 1:(28)P˜n+1=Pn+Ma0Xn+a2Vn+a3An+Ca1Xn+a4Vn+a5An.(2)
Get the displacement, velocity, and acceleration atn + 1:(29)Xn+1=K˜−1P˜n+1,An+1=a0Xn+1−Xn−a2Vn−a3An,Vn+1=Vn+a6An+a7An+1.For equation (25), the integral of hybrid method is(30)Xn+1=Xn+VnΔt+12+ψAnΔt2−ψAn−1Δt2,Vn+1=Vn+1+ϕAnΔt−ϕAn−1Δt,where Δt is the time step, n represents t = nΔt instant, n + 1 describes t = (n + 1)Δt instant, n − 1 is the instant of t = (n − 1)Δt, ψ, ϕ are instant parameters, and ψ = ϕ = 1/3. Take ϕ, for example, the reason for choosing parameters as follows:(31)An+1/2=32An−12An−1,An+1=2An−An−1,Vn+1=Vn+ΔtAn+1+4An+1/2+An/6=Vn+43AnΔt−13An−1Δt.For the two methods, we can find that the Newmark-β method needs to form the dynamic matrix, but does not require the diagonalization of the mass matrix. The improved hybrid integral method (31) does not require the formation of the dynamic matrix, and it can be solved directly from the differential equation. But the mass matrix is required to be the angular array, so this paper combines the implicit integral Newmark-β and explicit hybrid integral method to avoid the formation of large dynamic matrix, which could improve the efficiency of system modeling and the speed of solution.


## 3. Rotor Damage Detection Model

This section discusses how to build a classification model for rotor damage detection. According to the previous work’s experience, this paper applies the neural network methods of BP and RBF to train the frequency data obtained by the built couple dynamics model in Section2, and then the classification models are generated. The details are described in the following sections.


### 3.1. The Model Based on BP Neural Network

The neural network has various forms, and the backpropagation (BP) model is the most widely used. It is a one-way propagation multilayer forward neural network. Its typical structure is shown in Figure3. BP has one or more hidden nodes. Because there is no coupling on the same node, the output of each node only affects the output of the next node. Meanwhile, BP network is a highly nonlinear mapping from input to output, which is based on the gradient descent method. The learning of BP network is decomposed into two processes in this paper: in Phase 1 (forward propagation process), it gives input information and calculates the actual output value of each unit layer by layer through the hidden layer and the input layer. When the output layer fails to produce the desired output, it will switch to the backpropagation. In Phase 2 (backpropagation process), the difference between the actual output and the expected output is recursively calculated layer by layer, and then the weight is adjusted based on the difference. Repeat the above two processes until the error is adjusted to the error tolerance and stop learning.Figure 3
The basic structure of BP network.The neural network needs to be trained during the processes of function approximation and pattern recognition. The data sample is first required before training, including the input vectorp and the corresponding expected output vector D. Then, the weight and threshold are constantly adjusted during the training in order to minimize the error function of the neural network. Normally, the BP neural network error function is defined as the mean square error of the expected output vector D and the network output A.The hidden layer output is(32)yj=f∑iwjixi−θj=fnetj.The actual output of output layer is(33)zk=f∑iwkjyj−θl=fnetk,where xi, yj, and zk are input nodes, hidden layer nodes, and output nodes, respectively. Wji is the network weight of the input node and the hidden layer node, and Wkj is the network weight of the hidden layer node and the output node.If the expected output istk, the error is(34)E=12∑ktk−zk2=12∑ktk−f∑jwkjf∑iwjixi−θj−θk2.Next, it is implemented by correcting the weights and thresholds and deriving the hidden layer and output layer weights, respectively.(35)∂E∂Wkj=∑l=1n∂E∂Zl∂Zl∂Wkj=∂E∂Zk∂Zk∂Wkj,∂E∂Wji=∑l∑j∂E∂l∂Zl∂yj∂yj∂wji.Since the correction of the weightsΔWkj and ΔWji are proportional to the error function falling along the gradient,(36)Δwjk=−η∂E∂wkj,Δwji=−η′∂E∂wji,where η is the learning rate, and then the weight correction formula is(37)wkjk+1=wkjk+Δwjk,wjik+1=wjik+Δwji.The above learning process reflects that the correction of the network weight and threshold should follow the direction in which the performance function declines the fastest—the negative gradient direction.


### 3.2. The Model Based on RBF Neural Network

There are three layers including the input layer, the hidden layer with nonlinear neurons, and the output layer with linear neurons in RBF neural network [22–24]. It uses the radial basis function (RBF) as the basis of the hidden layer unit to form the hidden layer space. The hidden layer projects the input vector of the low-dimensional data into the high-dimensional space, which makes it linearly separable in high-dimensional space.The RBF network and the BP network are basically the same in structure, but the difference is that the RBF network hidden layer node transfer function is a radial basis function. The radial basis function appears in many forms, and the common one is a Gaussian function:(38)φix=exp−pi−ci2σi2,where ci is the center of i, σi is the parameter that controls the size of the receiving domain, and ‖⋅‖ is European norm.Thekth output node of the network completes the linear combination of the output of the hidden layer node:(39)Yk=∑iωk,iφi−bk,where bk is the domain value of the kth output node and ωk,j is the output weight of φi to Yk.Next, this paper carries out the RBF model by two stages, including the learning of the hidden layer and the learning of the output layer. Firstly, the hidden layer learning is performed on the determination of the basis function, the center, and width. It is usually unsupervised learning. This paper selects the center vectorCi and the normalized parameter σi with reference to the input sequence. Secondly, the learning of the output layer is reflected in the study of weights. The connection weight of the neurons between the hidden layer and the output layer can be directly calculated by the least square method, namely, the partial derivative of w is solved for the loss function.


### 3.3. The Damage Detection Process

With the support of BP and RBF, the damage detection process is shown in Figure4. The whole process is generally divided into three parts. Firstly, the training data need to be preprocessed by considering the wear range, which will be introduced in the next section. Secondly, the learning algorithms of BP and RBF are applied to build the damage detection models. Finally, the testing data are used to evaluate the models, and the better model can also be chosen in comparison with the results.Figure 4
Damage detection process.


### 3.4. Wear Range of Rotor Setting

Actually, the rotor damage detection can be recognized as two classification problems. We need to find whether the rotor is damaged. In this paper, the damaged situation is seen as the positive samples and vice versa, and it falls into negative samples. However, rotor wear is normal when the it is running in the engine. Hence, the wear range of the rotor needs to be set in order to ensure the positive samples and the negative samples in this paper. Considering the practical situation, the research designs the wear range of the rotor artificially. Generally speaking, the wear range of the rotor is between 3% and 5%, and the rotor ends are easy to wear. Here, the paper sets the lower limits of the loss of the rotor ends and other parts are 95% and 98%, respectively.


## 3.1. The Model Based on BP Neural Network

The neural network has various forms, and the backpropagation (BP) model is the most widely used. It is a one-way propagation multilayer forward neural network. Its typical structure is shown in Figure3. BP has one or more hidden nodes. Because there is no coupling on the same node, the output of each node only affects the output of the next node. Meanwhile, BP network is a highly nonlinear mapping from input to output, which is based on the gradient descent method. The learning of BP network is decomposed into two processes in this paper: in Phase 1 (forward propagation process), it gives input information and calculates the actual output value of each unit layer by layer through the hidden layer and the input layer. When the output layer fails to produce the desired output, it will switch to the backpropagation. In Phase 2 (backpropagation process), the difference between the actual output and the expected output is recursively calculated layer by layer, and then the weight is adjusted based on the difference. Repeat the above two processes until the error is adjusted to the error tolerance and stop learning.Figure 3
The basic structure of BP network.The neural network needs to be trained during the processes of function approximation and pattern recognition. The data sample is first required before training, including the input vectorp and the corresponding expected output vector D. Then, the weight and threshold are constantly adjusted during the training in order to minimize the error function of the neural network. Normally, the BP neural network error function is defined as the mean square error of the expected output vector D and the network output A.The hidden layer output is(32)yj=f∑iwjixi−θj=fnetj.The actual output of output layer is(33)zk=f∑iwkjyj−θl=fnetk,where xi, yj, and zk are input nodes, hidden layer nodes, and output nodes, respectively. Wji is the network weight of the input node and the hidden layer node, and Wkj is the network weight of the hidden layer node and the output node.If the expected output istk, the error is(34)E=12∑ktk−zk2=12∑ktk−f∑jwkjf∑iwjixi−θj−θk2.Next, it is implemented by correcting the weights and thresholds and deriving the hidden layer and output layer weights, respectively.(35)∂E∂Wkj=∑l=1n∂E∂Zl∂Zl∂Wkj=∂E∂Zk∂Zk∂Wkj,∂E∂Wji=∑l∑j∂E∂l∂Zl∂yj∂yj∂wji.Since the correction of the weightsΔWkj and ΔWji are proportional to the error function falling along the gradient,(36)Δwjk=−η∂E∂wkj,Δwji=−η′∂E∂wji,where η is the learning rate, and then the weight correction formula is(37)wkjk+1=wkjk+Δwjk,wjik+1=wjik+Δwji.The above learning process reflects that the correction of the network weight and threshold should follow the direction in which the performance function declines the fastest—the negative gradient direction.


## 3.2. The Model Based on RBF Neural Network

There are three layers including the input layer, the hidden layer with nonlinear neurons, and the output layer with linear neurons in RBF neural network [22–24]. It uses the radial basis function (RBF) as the basis of the hidden layer unit to form the hidden layer space. The hidden layer projects the input vector of the low-dimensional data into the high-dimensional space, which makes it linearly separable in high-dimensional space.The RBF network and the BP network are basically the same in structure, but the difference is that the RBF network hidden layer node transfer function is a radial basis function. The radial basis function appears in many forms, and the common one is a Gaussian function:(38)φix=exp−pi−ci2σi2,where ci is the center of i, σi is the parameter that controls the size of the receiving domain, and ‖⋅‖ is European norm.Thekth output node of the network completes the linear combination of the output of the hidden layer node:(39)Yk=∑iωk,iφi−bk,where bk is the domain value of the kth output node and ωk,j is the output weight of φi to Yk.Next, this paper carries out the RBF model by two stages, including the learning of the hidden layer and the learning of the output layer. Firstly, the hidden layer learning is performed on the determination of the basis function, the center, and width. It is usually unsupervised learning. This paper selects the center vectorCi and the normalized parameter σi with reference to the input sequence. Secondly, the learning of the output layer is reflected in the study of weights. The connection weight of the neurons between the hidden layer and the output layer can be directly calculated by the least square method, namely, the partial derivative of w is solved for the loss function.


## 3.3. The Damage Detection Process

With the support of BP and RBF, the damage detection process is shown in Figure4. The whole process is generally divided into three parts. Firstly, the training data need to be preprocessed by considering the wear range, which will be introduced in the next section. Secondly, the learning algorithms of BP and RBF are applied to build the damage detection models. Finally, the testing data are used to evaluate the models, and the better model can also be chosen in comparison with the results.Figure 4
Damage detection process.


## 3.4. Wear Range of Rotor Setting

Actually, the rotor damage detection can be recognized as two classification problems. We need to find whether the rotor is damaged. In this paper, the damaged situation is seen as the positive samples and vice versa, and it falls into negative samples. However, rotor wear is normal when the it is running in the engine. Hence, the wear range of the rotor needs to be set in order to ensure the positive samples and the negative samples in this paper. Considering the practical situation, the research designs the wear range of the rotor artificially. Generally speaking, the wear range of the rotor is between 3% and 5%, and the rotor ends are easy to wear. Here, the paper sets the lower limits of the loss of the rotor ends and other parts are 95% and 98%, respectively.


## 4. Model Verification

### 4.1. Coupled Dynamics Model Verification

#### 4.1.1. Rotor-Support Bearing Chock Fault Tester

The rotor-rolling bearing fault tester designed and manufactured can effectively simulate the rotor imbalance and the common faults of rolling bearings. The rotor-rolling bearing fault tester physical map and dynamics model are shown in Figure5. The rotor fault tester mainly comprises a rotating shaft, a rotor disc, a flange connecting plate, a bearing seat, a speed regulating motor, a gear speed increasing device, and the sensors of vibration displacement, rotational speed, and acceleration which can be flexibly mounted on the experimental device in order to conduct the comprehensive vibration test. The experimental device is a single-span and double-disc rotor. The vibration displacement of the turntable is measured by eddy current displacement sensors in both horizontal and vertical directions, and the rotational speed is measured by an eddy current sensor installed at the connection between the speed increaser and the rotating shaft. The sensor senses the displacement pulse caused by the rotation and obtains the rotation speed by counting. In the rotor-support coupling system model, there are two turntables, which include rotor disk P1 and flange disk P2 connected by gear output shaft. S1, S2 are the supports and L1, L2, L3, L4 are the position of the supports on the rotating shaft.Figure 5
Rotor-rolling bearing tester.


#### 4.1.2. Analysis and Verification

(i)
Experimental model analysis: in this paper, the modal analysis of the tester under natural support state is carried out by the hammering method. The main instruments include NI9234 dynamic signal acquisition module of NI (National Instrument) company, 30927 Hammer of ENDEVCO Company, and 4508 type ICP accelerometer from B&K Company. A single measurement point and multipoint incentive program is applied. The tester shaft is evenly tapped from the support to the support by 20 points, recorded as 1, …, 20. The acceleration sensor is placed at the fourth point, and the acceleration response signal is measured. Finally, the force signal and the acceleration signal are put into the collector, and then the modal parameters of the tester are obtained through the modal analysis software.
(ii)
Computational model analysis: for the FEM rotor-support coupling dynamics model, the rotation axis is divided into 48 units from the left end to the right end, and then 49 nodes are obtained, where support is at node 7, support is at node 43, and the hammer point is fixed at node 13. Then, the measurement will include the change from node 1 to node 49. In the calculation, the tapping method of the simulation experiment is to apply a pulse force on node 13 in order to obtain the acceleration response of each node of the rotor through simulation, and then the response values are input into the model analysis software to obtain the intrinsic frequency and the intrinsic modal of the system.
(iii)
The analysis results of calculation modal and experimental modal are compared as shown in Tables1–3. Table 1 shows rotor calculation parameters, Table 2 shows rolling bearing parameters, and Table 3 shows the support parameters.Table 1
Main calculation parameters of the rotor.
Parameters
Value
Elastic modulusE (1011 Pa)
2.1
Shaft diameterD (m)
0.019
Densityρ/103kg/m3
7.8
Poisson’s ratioμ
0.3
Proportional damping coefficientα0
5
Proportional damping coefficientα/10−5
1.35
L
1 (mm)
100
L
2 (mm)
342
L
3 (mm)
370
L
4 (mm)
130
PlateP1 mass mp (kg)
2.4
PlateP1 polar inertia moment Jdp1kg⋅m2
0.0125
PlateP1 equatorial inertia moment Jdd1kg⋅m2
0.00625
PlateP2 mass mp (kg)
0.45
PlateP2 polar inertia moment Jdp1kg⋅m2
0.00025
PlateP2 equatorial inertia moment Jdd1kg⋅m2
0.000125Table 2
Rolling bearing parameters.
Parameters
Value
Pitch diameterDmmm
36
Ball diameterd (mm)
9.6
Outer radiusr (mm)
22.8
Inner radiusr (mm)
13.2
Numbers of ballsNb
7
StiffnessCb/109N/m3/2
11.67
Bearing gapr0μm
0
Outer ring massmwkg
0.08
Axle-bearing massmbkg
76Table 3
Support parameters.
k
t/108 (N/m)
c
t (N·s/m)
k
f
/
10
7
N
/
m
c
f (N·s/m)
1.0
500
3.0
2000In Figure6, the blue point presents the experimental results obtained by the slapping method, the red line presents the simulation results obtained by the proposed coupled dynamics model, and the number values present the mean results of twenty groups’ experiment. The x-axis presents the order of the frequency and the y-axis presents the value of the frequencies. The comparative analysis of parameters for experimental modal and calculated modal shows that the simulation values of the coupled system are close to the experimental values, and the calculated vibration modes of the rotor are closer to the experimental vibration modes. The simulated values are similar to the experimental values. It proves the correct validity of the modeling method in this paper.Figure 6
Comparison of frequency identification results.


### 4.2. Rotor Damage Detection Model Verification

In this section, the numerical results are given to compare accuracies of the BP neural network and the RBF neural network. The two models were implemented in Matlab 2014a environment, and experiments were conducted on an Intel(R) Core(TM) i76700HQ 8-core CPU (main frequency: 2.60 GHz for each core) and 16 GB RAM.By the above mentioned couple dynamics model, we firstly gather the frequencies that can be divided into qualified samples and unqualified samples. Then, 10000 sets of data are used in the experiment. Next, 100 times’ ten-fold cross validation (70% of data are randomly selected as the training data and 30% of data as the test data) is used to train the models. Finally, the model performance is mainly validated by accuracy that presents the ratio of the number of correctly predicted samples to the total number of predicted samples, which can be calculated as follows:(40)accuracy=TP+TNTP+TN+TF+FN,where TP presents the prediction result is damaged (positive) and the practical result is damaged (positive), FP indicates the prediction result is damaged (positive) and the practical result is undamaged (negative), TN describes the prediction result is undamaged (negative) and the practical result is undamaged (negative), and FN shows the prediction result is undamaged (negative) and the practical result is damaged (positive).In practical work, the actual metrical error is less than 1%, and this paper adds 5% white noise for evaluating the robustness of the problem. In order to implement these models, we apply functions newrb()and sim()in the Matlab toolbox to create an approximate BP and RBF neural networks and predict outputs of network, respectively.In Table4, the results show the accuracy comparisons of RBF neural network and BP neural network under five kinds of white noise interference, and the RBF neural network performs slightly better detection than BP neural network when the noise data increase.Table 4
Accuracy comparisons of RBF neural network and BP neural network under 1–5% white noise interference.
1%
2%
3%
4%
5%
RBF neural network
99.26%
96.50%
91.23%
85.16%
73.42%
BP neural network
98.35%
94.93%
84.83%
70.44%
58.29%


## 4.1. Coupled Dynamics Model Verification

### 4.1.1. Rotor-Support Bearing Chock Fault Tester

The rotor-rolling bearing fault tester designed and manufactured can effectively simulate the rotor imbalance and the common faults of rolling bearings. The rotor-rolling bearing fault tester physical map and dynamics model are shown in Figure5. The rotor fault tester mainly comprises a rotating shaft, a rotor disc, a flange connecting plate, a bearing seat, a speed regulating motor, a gear speed increasing device, and the sensors of vibration displacement, rotational speed, and acceleration which can be flexibly mounted on the experimental device in order to conduct the comprehensive vibration test. The experimental device is a single-span and double-disc rotor. The vibration displacement of the turntable is measured by eddy current displacement sensors in both horizontal and vertical directions, and the rotational speed is measured by an eddy current sensor installed at the connection between the speed increaser and the rotating shaft. The sensor senses the displacement pulse caused by the rotation and obtains the rotation speed by counting. In the rotor-support coupling system model, there are two turntables, which include rotor disk P1 and flange disk P2 connected by gear output shaft. S1, S2 are the supports and L1, L2, L3, L4 are the position of the supports on the rotating shaft.Figure 5
Rotor-rolling bearing tester.


### 4.1.2. Analysis and Verification

(i)
Experimental model analysis: in this paper, the modal analysis of the tester under natural support state is carried out by the hammering method. The main instruments include NI9234 dynamic signal acquisition module of NI (National Instrument) company, 30927 Hammer of ENDEVCO Company, and 4508 type ICP accelerometer from B&K Company. A single measurement point and multipoint incentive program is applied. The tester shaft is evenly tapped from the support to the support by 20 points, recorded as 1, …, 20. The acceleration sensor is placed at the fourth point, and the acceleration response signal is measured. Finally, the force signal and the acceleration signal are put into the collector, and then the modal parameters of the tester are obtained through the modal analysis software.
(ii)
Computational model analysis: for the FEM rotor-support coupling dynamics model, the rotation axis is divided into 48 units from the left end to the right end, and then 49 nodes are obtained, where support is at node 7, support is at node 43, and the hammer point is fixed at node 13. Then, the measurement will include the change from node 1 to node 49. In the calculation, the tapping method of the simulation experiment is to apply a pulse force on node 13 in order to obtain the acceleration response of each node of the rotor through simulation, and then the response values are input into the model analysis software to obtain the intrinsic frequency and the intrinsic modal of the system.
(iii)
The analysis results of calculation modal and experimental modal are compared as shown in Tables1–3. Table 1 shows rotor calculation parameters, Table 2 shows rolling bearing parameters, and Table 3 shows the support parameters.Table 1
Main calculation parameters of the rotor.
Parameters
Value
Elastic modulusE (1011 Pa)
2.1
Shaft diameterD (m)
0.019
Densityρ/103kg/m3
7.8
Poisson’s ratioμ
0.3
Proportional damping coefficientα0
5
Proportional damping coefficientα/10−5
1.35
L
1 (mm)
100
L
2 (mm)
342
L
3 (mm)
370
L
4 (mm)
130
PlateP1 mass mp (kg)
2.4
PlateP1 polar inertia moment Jdp1kg⋅m2
0.0125
PlateP1 equatorial inertia moment Jdd1kg⋅m2
0.00625
PlateP2 mass mp (kg)
0.45
PlateP2 polar inertia moment Jdp1kg⋅m2
0.00025
PlateP2 equatorial inertia moment Jdd1kg⋅m2
0.000125Table 2
Rolling bearing parameters.
Parameters
Value
Pitch diameterDmmm
36
Ball diameterd (mm)
9.6
Outer radiusr (mm)
22.8
Inner radiusr (mm)
13.2
Numbers of ballsNb
7
StiffnessCb/109N/m3/2
11.67
Bearing gapr0μm
0
Outer ring massmwkg
0.08
Axle-bearing massmbkg
76Table 3
Support parameters.
k
t/108 (N/m)
c
t (N·s/m)
k
f
/
10
7
N
/
m
c
f (N·s/m)
1.0
500
3.0
2000In Figure6, the blue point presents the experimental results obtained by the slapping method, the red line presents the simulation results obtained by the proposed coupled dynamics model, and the number values present the mean results of twenty groups’ experiment. The x-axis presents the order of the frequency and the y-axis presents the value of the frequencies. The comparative analysis of parameters for experimental modal and calculated modal shows that the simulation values of the coupled system are close to the experimental values, and the calculated vibration modes of the rotor are closer to the experimental vibration modes. The simulated values are similar to the experimental values. It proves the correct validity of the modeling method in this paper.Figure 6
Comparison of frequency identification results.


## 4.1.1. Rotor-Support Bearing Chock Fault Tester

The rotor-rolling bearing fault tester designed and manufactured can effectively simulate the rotor imbalance and the common faults of rolling bearings. The rotor-rolling bearing fault tester physical map and dynamics model are shown in Figure5. The rotor fault tester mainly comprises a rotating shaft, a rotor disc, a flange connecting plate, a bearing seat, a speed regulating motor, a gear speed increasing device, and the sensors of vibration displacement, rotational speed, and acceleration which can be flexibly mounted on the experimental device in order to conduct the comprehensive vibration test. The experimental device is a single-span and double-disc rotor. The vibration displacement of the turntable is measured by eddy current displacement sensors in both horizontal and vertical directions, and the rotational speed is measured by an eddy current sensor installed at the connection between the speed increaser and the rotating shaft. The sensor senses the displacement pulse caused by the rotation and obtains the rotation speed by counting. In the rotor-support coupling system model, there are two turntables, which include rotor disk P1 and flange disk P2 connected by gear output shaft. S1, S2 are the supports and L1, L2, L3, L4 are the position of the supports on the rotating shaft.Figure 5
Rotor-rolling bearing tester.


## 4.1.2. Analysis and Verification

(i)
Experimental model analysis: in this paper, the modal analysis of the tester under natural support state is carried out by the hammering method. The main instruments include NI9234 dynamic signal acquisition module of NI (National Instrument) company, 30927 Hammer of ENDEVCO Company, and 4508 type ICP accelerometer from B&K Company. A single measurement point and multipoint incentive program is applied. The tester shaft is evenly tapped from the support to the support by 20 points, recorded as 1, …, 20. The acceleration sensor is placed at the fourth point, and the acceleration response signal is measured. Finally, the force signal and the acceleration signal are put into the collector, and then the modal parameters of the tester are obtained through the modal analysis software.
(ii)
Computational model analysis: for the FEM rotor-support coupling dynamics model, the rotation axis is divided into 48 units from the left end to the right end, and then 49 nodes are obtained, where support is at node 7, support is at node 43, and the hammer point is fixed at node 13. Then, the measurement will include the change from node 1 to node 49. In the calculation, the tapping method of the simulation experiment is to apply a pulse force on node 13 in order to obtain the acceleration response of each node of the rotor through simulation, and then the response values are input into the model analysis software to obtain the intrinsic frequency and the intrinsic modal of the system.
(iii)
The analysis results of calculation modal and experimental modal are compared as shown in Tables1–3. Table 1 shows rotor calculation parameters, Table 2 shows rolling bearing parameters, and Table 3 shows the support parameters.Table 1
Main calculation parameters of the rotor.
Parameters
Value
Elastic modulusE (1011 Pa)
2.1
Shaft diameterD (m)
0.019
Densityρ/103kg/m3
7.8
Poisson’s ratioμ
0.3
Proportional damping coefficientα0
5
Proportional damping coefficientα/10−5
1.35
L
1 (mm)
100
L
2 (mm)
342
L
3 (mm)
370
L
4 (mm)
130
PlateP1 mass mp (kg)
2.4
PlateP1 polar inertia moment Jdp1kg⋅m2
0.0125
PlateP1 equatorial inertia moment Jdd1kg⋅m2
0.00625
PlateP2 mass mp (kg)
0.45
PlateP2 polar inertia moment Jdp1kg⋅m2
0.00025
PlateP2 equatorial inertia moment Jdd1kg⋅m2
0.000125Table 2
Rolling bearing parameters.
Parameters
Value
Pitch diameterDmmm
36
Ball diameterd (mm)
9.6
Outer radiusr (mm)
22.8
Inner radiusr (mm)
13.2
Numbers of ballsNb
7
StiffnessCb/109N/m3/2
11.67
Bearing gapr0μm
0
Outer ring massmwkg
0.08
Axle-bearing massmbkg
76Table 3
Support parameters.
k
t/108 (N/m)
c
t (N·s/m)
k
f
/
10
7
N
/
m
c
f (N·s/m)
1.0
500
3.0
2000In Figure6, the blue point presents the experimental results obtained by the slapping method, the red line presents the simulation results obtained by the proposed coupled dynamics model, and the number values present the mean results of twenty groups’ experiment. The x-axis presents the order of the frequency and the y-axis presents the value of the frequencies. The comparative analysis of parameters for experimental modal and calculated modal shows that the simulation values of the coupled system are close to the experimental values, and the calculated vibration modes of the rotor are closer to the experimental vibration modes. The simulated values are similar to the experimental values. It proves the correct validity of the modeling method in this paper.Figure 6
Comparison of frequency identification results.


## 4.2. Rotor Damage Detection Model Verification

In this section, the numerical results are given to compare accuracies of the BP neural network and the RBF neural network. The two models were implemented in Matlab 2014a environment, and experiments were conducted on an Intel(R) Core(TM) i76700HQ 8-core CPU (main frequency: 2.60 GHz for each core) and 16 GB RAM.By the above mentioned couple dynamics model, we firstly gather the frequencies that can be divided into qualified samples and unqualified samples. Then, 10000 sets of data are used in the experiment. Next, 100 times’ ten-fold cross validation (70% of data are randomly selected as the training data and 30% of data as the test data) is used to train the models. Finally, the model performance is mainly validated by accuracy that presents the ratio of the number of correctly predicted samples to the total number of predicted samples, which can be calculated as follows:(40)accuracy=TP+TNTP+TN+TF+FN,where TP presents the prediction result is damaged (positive) and the practical result is damaged (positive), FP indicates the prediction result is damaged (positive) and the practical result is undamaged (negative), TN describes the prediction result is undamaged (negative) and the practical result is undamaged (negative), and FN shows the prediction result is undamaged (negative) and the practical result is damaged (positive).In practical work, the actual metrical error is less than 1%, and this paper adds 5% white noise for evaluating the robustness of the problem. In order to implement these models, we apply functions newrb()and sim()in the Matlab toolbox to create an approximate BP and RBF neural networks and predict outputs of network, respectively.In Table4, the results show the accuracy comparisons of RBF neural network and BP neural network under five kinds of white noise interference, and the RBF neural network performs slightly better detection than BP neural network when the noise data increase.Table 4
Accuracy comparisons of RBF neural network and BP neural network under 1–5% white noise interference.
1%
2%
3%
4%
5%
RBF neural network
99.26%
96.50%
91.23%
85.16%
73.42%
BP neural network
98.35%
94.93%
84.83%
70.44%
58.29%


## 5. Conclusions

This paper studied the coupling dynamics model of the engine vibration and damage detection model of rotor. For the coupling dynamics model, the FEM is used to model the rotor and the casing system firstly, and then a complex structure with multirotor and multi-causing is proposed by defining support system and linking methods. Next, the numerical integral method is used to obtain the inherent frequency of the nonlinear dynamic system. Finally, the model experiment is carried out by the rotor tester, and the results prove the correctness and effectiveness of the couple dynamics model.After the coupled dynamics model is established, two damage detection models of the rotor by BP and RBF neural network are built. The experimental results show that the damage detection models have good robustness and RBF method is better while the data noise increases.In the future, we will optimize the initialized connection weights and threshold values and apply the deep learning techniques to solve the fault diagnosis of practical engine.



---

*Source: 5907613-2020-03-31.xml*