Dynamical Behaviors Analysis of the Rotor Model with Coupling Faults and Applications of the TPOD Method

: The transient proper orthogonal decomposition (TPOD) method is applied for order reduction in the rotor-bearing system with the coupling faults in this paper. A 24 degrees of freedom (DOFs) rotor model supported by a pair of sliding bearings with both crack and rub-impact faults is established by the discrete modeling method. The complexity of dynamic behaviors of the rotor system with the coupling faults is discussed via the comparison of the rotor system with the single fault (crack or rub-impact). The proper orthogonal mode (POM) energy method is proposed to conﬁrm the DOF number of the reduced model. The TPOD method is used in the coupling faults system to obtain the optimal order reduction model based on the POM energy. The e ﬃ ciency of the order reduction method is veriﬁed by comparing the bifurcation behaviors between the original and the reduced system. The TPOD method provides the optimal order reduction model to study the non-linear dynamic characteristics of the complex rotor system with the coupling faults.


Introduction
The influence of the faults in the rotating machines on the associated dynamic behavior has been a focus of attention for many researchers. The presence of the faults may lead to a dangerous and catastrophic effect on the dynamic behavior of rotating structures and cause serious damage to the rotating machineries. In general, the faults in the rotor system contain the crack fault [1][2][3][4][5], the rub-impact of rotor-to-stator fault [6][7][8][9], the pedestal looseness fault [10][11][12][13], the ball bearing fault [14,15], the misalignment fault, etc. [16][17][18][19]. On the basis of the common faults in the rotor systems, the coupling faults (the looseness and rub-impact, the crack and rub-impact, etc.) usually occur in the actual rotor systems. The amplitude varies when the cracked rotor is operating, which will lead to the rub-impact with other units of the rotating machine. However, if the rub-impact fault occurs, the crack of shaft or impeller may occur on account of the harsh operating conditions and the strong influence. Few researchers studied the rotor systems with the coupling faults, since they are more complex than the systems with only single fault [20]. The dynamical characteristics of the rotor systems with the complex coupling faults should be further explored.
The construction processes are listed as follows: (1) Provide the initial conditions, proceed numerical simulation and obtain displacement information from the transient process of various DOFs, which is denoted by z 1 (t), z 2 (t), . . . z M (t). M is the number of DOFs in the system, recording transient time interval displacement sequence of all DOFs z i = (z i (t 1 ), z i (t 2 ), . . . z i (t N )) T , i = 1, . . . , M. The number of time interval is N and the interval is equal to each other, the time series form the matrix X = [z 1 , z 2 , . . . , z M ], the order of X is N × M. When calculating the self-correlation matrix T = X T X, the order of matrix is M × M.
The eigenvector of the matrix can therefore be denoted by ϕ 1 , ϕ 2 , . . . , ϕ M , and the corresponding eigenvalues are λ 1 > λ 2 > . . . > λ M . (2) U is the matrix formed by the first n orders of T = X T X and it can also indicate that the matrix U contains the first n largest eigenvalues of T. The order of U is M × n, so we know that the order of U T U is n × n. Taking the coordinates transformation on system coordinates Z, acquiring the new coordinates P, Z = UP, substituting Z into Equation (1), gaining the Equation (2): Due to the column vector of U is orthogonal and nonzero, and the matrix U T U is n order diagonal and full rank, hence there exits inverse matrix of U T U. For the Equation (2), taking (U T U) −1 U T left multiplication on both sides, we can obtain (3): ..
Setting C R = −(U T U) −1 U T CU, K R = (U T U) −1 U T KU, F R = (U T U) −1 U T F, then we receive the Formula (4): ..
In this way, the original system is transformed into the n-DOFs reduced system by the TPOD method.

The Rotor Model with Coupling Faults
The high-pressure rotor system supported by a pair of sliding bearings with both crack and rub-impact faults is modeled as a 24-DOFs system by the Newton's second law, as shown in Figure 1. Here we assume that the axial and torsional vibrations of the system and the gyroscopic moment are neglectful. O i (i = 1 . . . 12) are the geometric centers of the discs. m i (i = 1 . . . 12) are the equivalent lumped masses. c i (i = 1 . . . 12) are the equivalent damping coefficients at the position of the lumped masses. k i (i = 1 . . . 11) are the equivalent stiffness of the corresponding discs. Assume that the crack occurs on the shaft between the first and second disc. The rotor model is supported by a pair of liquid-film bearings on both ends.
The schematic diagram of the cracked section is shown in Figure 2a. α is the depth of crack, β is the intersection angle between the unbalance and the crack normal vector. Law of open and close describes the effects of crack intersection angle to the shaft stiffness. There are three classic models that characterizing the time-varying property of a crack when the rotor is operating: the cosine model [53], square wave model and hybrid one [54]. In Figure 2a, the initial intersection angle of the cracked normal vector and x-axis is 0. We consider the case of gravity dominance, the open and close function of crack f (ϕ) can be expressed as the function of rotating angle ωt. k is the cracked stiffness of rotating shaft, ∆k ξ and ∆k η are the stiffness variations of cracked normal and tangential vectors. The stiffness of rotating shaft with crack is expressed as: ∆k ξ cos 2 (ωt) + ∆k η sin 2 (ωt) ∆k ξ − ∆k η sin(ωt) cos(ωt) ∆k ξ − ∆k η sin(ωt) cos(ωt) ∆k ξ sin 2 (ωt) + ∆k η cos 2 (ωt) The expression of f (ϕ) is shown in Formula (6), this model can select the square wave and cosine model based on the depth of the crack respectively.
The rub-impact model is shown in Figure 2b, it is combined with linear contact force and Coulomb friction [55], the formula is expressed as: where r is radial displacement, r 0 is clearance between rotor and stator, k b is rub-impact stiffness, µ 0 is frictional coefficient, P x and P y represent the rub-impact forces along x and y directions. The dynamical equation is shown in Formula (8) and the dimensionless process is expressed as follows: . .. .. . ..
, M 12 = m 12 cω 2 sP X i (i = 1, . . . 12) and Y i (i = 1, . . . 12) represent the vibration displacement of x and y directions. The dimensionless equation is expressed as Formula (9). f x and f y are the model of dimensionless are the S(x, y, α) = x cos α+y sin α 1−(x cos α+y sin α) 2 , α directional components of bearing non-linear oil-film force. c is the bearing clearance, is the Sommerfeld number. µ is the lubricating oil viscosity, L is the bearing length, R is the radius of the bearing, ω is the external excitation, P is the loading, and τ is the dimensionless time. The detailed values of the corresponding parameters of the rotor system are shown as in the Appendix A. The non-linear oil-film force [56] along x and y directions can be found in Appendix B. We also provide the formulas of the corresponding parameters in Appendix B. .. ..
Equation (9) is the dimensionless form of Equation (8). The dimensionless Equation (9) can be written as Equation (10), in which z = (x 1 , y 1 , · · · , x 12 , y 12 ) T , the expression of c, k, f are shown as above: To facilitate the theory analysis, we implement the Taylor series expansion of the oil-film force, thus α can be rewritten as: For the convenience of calculation, Equation (10) is written as Equation (12) briefly: ..
In Equation (12), C is the damping matrix, K is the stiffness matrix, F is the force vector, which includes oil-film and external excitation etc. Z = [z 1 z 2 . . . z 24 ] T corresponds to [x 1 y 1 . . . x 12 y 12 ] T in Equation (10).

Discussion of the Dynamical Characteristics and Effects of Systematic Parameters
In this section, we will study the non-linear dynamical characteristics of the rotor with the coupling faults by comparing the coupling faults with the single fault (crack or rub-impact). Then we will discuss the effects of the systematic parameters to the coupling fault. Here we define the initial values as z i = 0(i = 1, 2, . . . , 24), z i = 0.001(i = 1, 2, . . . , 24, i 5, 6), and the integral step is π/256. Figure 3 shows the four amplitude-frequency curves of the rotor systems with different faults: the original system without faults; the cracked fault; the rub-impact fault and the coupling fault. 1/2 sub-harmonic vibration occurs in the cracked and coupling system at the frequency of 2ω c , ω c is the natural frequency. The amplitude of the coupling system is apparently higher than the system with only crack fault, which indicates that the rub-impact will increase the response amplitude of the crack. The amplitude of the system with rub-impact bumps up at the frequency of 1000 rad/s(nearby 3ω c ), this result illustrates rub-impact will aggravate the complex motion of the system.  It is clear that the bifurcation behaviors vary with different rotor faults. The obvious period-doubling bifurcation occurs at the second order natural frequency of the cracked fault system. The period-doubling bifurcation also occurs in the coupling fault system, but the amplitude is much larger. When the rotating speed is 1000 rad/s (nearby 3ω c ), the amplitude of the rub-impact fault system suddenly increases, which indicates that the rub-impact fault can give rise to the complex motion of the system. When ω = 350 rad/s (nearby the primary resonance), the time history curves of the four systems are drawn clearly in Figure 5. The curves contain both the transient and steady-state processes. The transient process of the oil-filmed system (Figure 5a) vanishes in about two periods, the other three systems (Figure 5b-d) need three periods. This indicates that the stability of the system with faults can't be kept in comparison to the normal system. Meanwhile, in Figure 5d, the steady-state signal of the coupling fault is not standard sinusoidal, which has a visible difference in the other three systems. In Figure 6, it shows the trajectories of the orbit of shaft center of the four systems. The motion curves in (a) and (c) are all elliptical. The major axis and minor axis of ellipse of the oil-filmed system is longer than the rub-impact system. The trajectories of the orbit of shaft center of cracked (b) and coupling (d) faults are the closed curves similar to ellipse, but the bottom is a little narrow. The crack fault is the main reason to change the ellipse of trajectories of the orbit of shaft center as closed curves similar to the ellipse near the primary resonance.

Discussion of the Dynamical Characteristics
We use the comparison of the phase portraits ( Figure 7) to verify the results in Figure 6 based on the same initial values and rotating speed. In Figure 7a,c, it is easy to discover that the phase portraits are both the closed curves similar to circle, and the diameter of the oil-filmed system is larger. In the other two figures (Figure 7a,c), the curves are both irregularly closed.     As is clearly seen in Figure 10a,c, the phase portraits of the oil-filmed and rub-impact systems are the curves similar to ellipse. The curves in Figure 10b,d are both double-loop, the two loops are closer in the cracked system relative to the coupling system. This result also verifies that the crack fault can lead to double-loop curves, the rub-impact fault affects the topological structures of phase portraits.  Figures 8-10, the crack fault affects more to the rotor system at the frequency of 2ω c . The signal of cracked fault system is more complex, the double-loop curves occur in both trajectories of the orbit of shaft center and phase portraits, the period-doubling bifurcation occurs in Figure 4b. The rub-impact fault will not affect the natural structures of the cracked characteristics. Hence, the dynamical characteristics near the frequency of 2ω c can be regarded as the important judgments to verify the crack fault of the rotor-bearing system.

Remark 2. In
In Figure 11, the four figures (Figure 11a-d) show the time histories of the rotor-bearing models at the frequency of 3ω c . There exists a certain difference between the four faults of the rotor system. The difference of the time history can provide theory guidance to the fault detection of the rotor system. As is shown in Figure 12, it shows the orbit of shaft center diagrams of the oil-film, crack, rub-impact, coupling fault when ω = 1020 rad/s (the third order natural frequency). Based on the constraint of the displacement of rub-impact fault, the orbit of shaft center is easier than the other three curves. However, the displacement is the largest, which indicates that the rub-impact fault arouses the larger motion near the third order frequency. The orbit of shaft center of the crack fault is more complex than that of the oil-film fault. The complexity of the orbit of shaft center of the coupling fault is situated between the oil-film fault and crack fault.  Figure 13 shows phase diagrams of oil-film, crack, rub-impact and the coupling systems at the frequency of ω = 1020 rad/s. Oil-film rotor causes complex motion at this time. The phase diagram of the rub-impact fault is annular narrow band, it is a quasi-periodic attractor. The sort of the four fault cases can be given based on the complexity of the attractor: crack, coupling, oil-film and rub-impact. Figures 11-13, we can obtain the conclusions as follows: the crack fault can arouse the 1/2 sub-harmonic vibration, the violent vibration at the frequency of 3ω c can be aroused by the rub-impact fault. At the frequency of primary resonance, the rub-impact and coupling fault may give rise to the reducing of the first order critical speed, while the crack and coupling fault may cause the structural variation of trajectories of the orbit of shaft center and phase portraits. The crack can restrain the vibration amplitude at the frequency of 3ω c in the coupling system, the crack also arouses the trajectories of the orbit of shaft center and attractor more complex.

The Effects of Systematic Parameters
We will discuss the effects of systematic parameters of the rotor system to the dynamical characteristics of the coupling fault. The systematic parameters include cracked depth, clearance between rotor and stator, stator stiffness, eccentricity. The bifurcation and amplitude-frequency behaviors will be provided as the values of the systematic parameters vary.

The Cracked Depth
The effect of the cracked depth variation on dynamical characteristics of the coupling fault is investigated in this section. As is shown in Figure 14a,b, they reflect the bifurcation diagrams of the coupling fault systems with the cracked depth a = 0.6R and a = R. Figure 15 shows the comparison of frequency-amplitude curves of two cracked depths. As the cracked depth decreases, the stiffness variation decreases correspondingly, so the first order critical speed of the coupling fault rotor slightly increases. The period-doubling bifurcation vanishes at the frequency of 2ω c , the 1/2 sub-harmonic vibration amplitude of the relevant amplitude-frequency characteristic decreases at the same time ( Figure 15). Since the crack fault characteristic does not play the main role and even vanishes, the rub-impact fault emerges, the vibration amplitude suddenly increases at the frequency of 3ω c . This change demonstrates the crack is the main factor to generate 1/2 sub-harmonic vibration in the coupling system again, there will be more effect with the increase of the cracked depth. The rub-impact fault will result in the amplitude of the coupling fault peaking at the frequency of 3ω c , there is reversible inhibition of the crack to the amplitude aroused by the rub-impact. The restrain will be stronger as the cracked depth increases. When ω ∈ (1100, 1200), the cracked depth will have little impact to the amplitude based on the comparison of the amplitude-frequency curves ( Figure 15). In comparison to the bifurcation diagrams (Figure 14), the complex motion state of the system changes.

Clearance between Rotor and Stator
This section discusses the effects of clearance between rotor and stator to the coupling fault. Similarly, Figures 16 and 17 present the bifurcation and amplitude-frequency when the clearance is 0.12 and 0.18, respectively. When the clearance varies, the clearance has little impact on the vibration amplitude and critical speed of the coupling system. Near the frequency of 2ω c , the amplitude of 1/2 sub-harmonic resonance decreases as the clearance decreases, the amplitudes of the systems with different clearance are similar when ω ∈ [2ω c , 3ω c ]. The amplitude will have a sharp increase as the clearance decreases at the frequency of 3ω c , the rub-impact fault characteristic enhanced obviously, which once verifies that the sudden increase of the 3ω c amplitude of the coupling fault system is primarily caused by the rub-impact fault. At the high frequency (over 1100 rad/s), the amplitude increases a little when the clearance decreases. The topological structures of the bifurcation diagram will be more complex in comparison to Figure 16b.

Stator Stiffness
The effects of the stator stiffness variation to the bifurcation behaviors and amplitude-frequency are studied in this section. We select two different parameter values to discuss the comparison of bifurcation ( Figure 18) and amplitude-frequency (Figure 19) of the coupling fault system: the first is 500 N/m, the second is 1000 N/m. In Figure 19, we can see clearly that the variation of the stator stiffness has little impact on the primary frequency, vibration amplitude and critical speed. 1/2 sub-harmonic resonance amplitude increases as the stiffness decreases. At the frequency of 3ω c , the amplitude changes obviously when the stiffness decreases. The complex motion region is simpler when the stiffness is 500 N/m.

Eccentricity
Finally, we study the effects of the eccentricity variation to the coupling fault system. Figures 20 and 21 present the comparison of bifurcation and amplitude-frequency of two different eccentricities (0.00015 mm and 0.00012 mm). As is shown in Figure 21, the amplitude increases when the eccentricity magnifies at the frequency of primary resonance. 1/2 sub-harmonic resonance is more obvious near the frequency of 2ω c , it is clear that the annulus of the period-doubling bifurcation is larger. When the frequency is over 2ω c , the amplitude decreases as the eccentricity increases, this result indicates that increasing the eccentricity appropriately can improve the stability of the system under a certain speed. However, in the high frequency region, the bifurcation behaviors will be more complex when the eccentricity increases, in other words, the non-linear dynamic characteristics will be more complex.
As is discussed above, the eccentricity of the system has more impact to dynamical characteristic of the amplitude-frequency and bifurcation at the frequencies of primary resonance, 2ω c , 3ω c , high frequency. The eccentricity may affect the visualizations of different faults at each frequency range.  Remark 5. The systematic parameters play a key role in the rotor models, the variations of the parameters will affect the dynamical characteristics of the rotor-bearing system. We make a preliminary study on the effects of cracked depth, clearance between rotor and stator, stator stiffness, eccentricity to the bifurcation diagrams and amplitude-frequency curves. This research is only a beginning to discuss the variations of the dynamical characteristics with the systematic parameters. We can also use the frequency-spectrum analysis to describe the effects of the variations of the parameters.

The Optimal Reduced Model Based on the POM Energy
The first few order POMs pertaining to the rotating speed of the rotor will be discussed to confirm the dimension of the reduced model in this section.
In Figure 22, it shows the first three order POM energy curve of the 24-DOFs rotor model with the coupling fault. It is clear that the three order POMs occupy 98.7% energy when the rotating speed is about 1120 rad/s, so we can obtain the optimal order reduction model at this speed based on the POM energy method. The POM energy method has been applied in Ref. [50], then we use the TPOD method in Section 2.1 to reduce the 24-DOFs rotor model to a 3-DOFs one. The POM energy method is applied to confirm the dimension of the order reduction model. Then we use the TPOD method to get the optimal reduced model at the corresponding rotating speed.

The Efficiency of the Order Reduction Method
In this section, we will discuss the efficiency of the order reduction method via the bifurcation analysis of the original and the reduced system.
On the basis of the POM energy curve (Figure 22), we can obtain the 3-DOFs reduced model when the rotating speed is 1120 rad/s. As is shown in Figure 23, we can see clearly that the reduced system (Figure 23b) preserves the period-doubling bifurcation, the rotating speed of complex motion and the topological structure of the original one (Figure 23a). The comparison of the bifurcation diagrams verifies the efficiency of the TPOD method, we can also use other dynamical behaviors (phase portraits, amplitude-frequency curve, etc.) to discuss the efficiency [49,50], we don't state one by one here.

Remar 7.
The reduced order reduction model obtained by the TPOD method preserves the dynamical characteristics of the original one well. The TPOD method is proposed by the present authors, it is applied in the rotor models supported by the sliding bearings or ball bearings with looseness. The efficiency of the order reduction method applied in the rotor system with coupling fault further generalizes this method. The Runge-Kutta method is used to calculate the dynamic response, and we can obtain the results in Figures 3-23. It is promising to adapt this method to the uncertainty quantification of faulted rotor systems with parametric variabilities [57][58][59][60]. Fortran Language is used to solve the dynamic equation.

Conclusions and Outlooks
The non-linear dynamic characteristics of the rotor system with the coupling faults have been discussed in this paper. A 24-DOFs rotor system supported by a pair of sliding bearings with both crack and rub-impact faults has been established by the Newton's second law. The coupling faults have been studied via comparing the single fault (crack or rub-impact). The effects of the systematic parameters (cracked depth, clearance between rotor and stator, stator stiffness, eccentricity) to the coupling faults have also been presented by comparing with the bifurcation and amplitude-frequency characteristics. The POM energy method has been used to confirm the DOF number of the reduced system. The TPOD has been applied to reduce the original system to a 3-DOFs one, the efficiency of the optimal order reduction model has been verified via the comparison of the bifurcation diagrams. The POM energy is a direct method to obtain the optimal reduced system of high-dimensional non-linear dynamic system. The dynamic characteristics analysis of the rotor system with coupling faults and the applications of the TPOD method can provide qualitative guidance to the common faults in rotor system and prior information of fault diagnosis.
Further studies on this subject are being carried out by the present authors in the three aspects: the first is to study the dynamical characteristics of the ball bearing rotor model with the coupling faults; the second is to study the coupling faults of the rotor system with uncertainties based on the polynomial dimensional decomposition and other methods.

Conflicts of Interest:
The authors declare that they have no conflict of interest. geometric centers of the discs m i (i = 1 . . . 12) equivalent lumped masses c i (i = 1 . . . 12) equivalent damping coefficients k i (i = 1 . . . 11) equivalent stiffness ω rotating speed α depth of crack β intersection angle between unbalance and crack normal vector k stiffness of cracked rotating shaft ∆k ξ ,∆k η stiffness variations of cracked normal and tangential vectors r radial displacement r 0 clearance between rotor and stator k b rub-impact stiffness µ 0 frictional coefficient P x , P y rub-impact forces along x and y directions X i , Y i (i = 1, . .