Nonlinear Dynamics of a Blade Rotor with Coupled Rubbing of Labyrinth Seal and Tip Seal

*e dynamic response and its stability of a blade rotor with coupled rubbing in the labyrinth seal and tip seal are investigated.*e dynamic equations are established based on the Hertz contact rubbing force at the labyrinth seal and the tip rubbing force considering both the contact deformation of the tip seal and the bending deformation of the blade. Numerical simulations show that the coupled rubbing response includes periodic motions, almost periodic motions, and chaotic motions. Compared with the single rubbing fault, coupled rubbing increases the range of rotation velocity of contact. A new continuation shooting method is used in the solution and stability analysis of the periodic response to give the bifurcation diagrams. *e paths of the system for entering and exiting chaos are analyzed.


Introduction
To improve the operational efficiency of the turbo machinery, the clearance between the rotor and stator is designed to be smaller and smaller, which also increases the possibility of rubbing. Severe rubbing may cause breakage of the blade and instability of the rotor, which eventually lead to serious operational accidents [1].
ere are usually two kinds of seals for turbo machines: labyrinth seals and tip seals. e rubbing of labyrinth seals is usually described as a rubbing problem between two cylindrical surfaces. e contact forces may be expressed as a linear or nonlinear function of the deformation of interaction. e friction forces are given by the Coulomb law of friction. Many researchers focused on related problems, such as forward periodic motion [2][3][4], reverse precession [5,6], chaos [7,8], and transient response [9]. Other interesting rubbing problems have also attracted people's attention. Patel [10,11] investigated a rotor with coupled faults of fatigue crack and rotor-stator rub. Chang-Jian et al. [12] focused on the response of a flexible rub-impact rotor supported by oil film bearings. Hou et al. [13,14] presented the effects of maneuver flight on the rotor/seal rubbing response. Cao et al. [15] investigated a rub-impact rotor system with fractional order damping. Yang et al. [16] studied the dynamic behaviors of a rotor-casing system subjected to axial load and radial rub. Sun et al. [17] considered a dual-rotor system induced by impact rub. Tai et al. [18] studied the periodic response of the rotor system with single-point impact rub. Roques et al. [19] investigated speed transients with rotor-to-stator rubbing caused by an accidental blade-off imbalance. Jin et al. [20] studied a dualrotor-bearing-coupling misalignment system with bladecasing rubbing numerically and experimentally to reveal the nonlinear vibration characteristics of an actual aero-engine rotor system. Wang et al. [21] considered the stick-slip transitions, which are related to the dynamical scenarios of sliding bifurcations, of the self-excited dry friction backward whirls of a general rotor/stator rubbing system. Hildebrandt et al. [22] discussed the heat flux distribution problem between bristle package and rotor during brush seal rubbing event based on a 3D CFD model.
When rubbing happens at the tip seal, the number of blades in contact and the deformation of each blade is different at each moment, which is more likely to lead to a complex dynamic response. For the tip rubbing, by assuming the blade as a cantilever beam, Choy [9] and Padovan [23] derived the relationship between the normal contact force and the radial deformation of the blade and investigated the nonlinear dynamics of rotor/casing rub interactions. Jiang et al. [24] proposed a rubbing contact force model in the form of periodic pulse force and derived the expression of the maximum contact force through theoretical calculation. Lesaffre et al. [25] studied the blademachine impact friction stability, and the rubbing was achieved by using an additional stiffness matrix and damping matrix. Padova [26] conducted an experimental study to analyze the interaction between the blade and casing with different incursion depths. Kim [27] studied the rubbing of the asymmetric blade rotor and the anisotropic casing. Sinha [28] investigated the dynamics of an asymmetric turbofan rotor with fan blade-loss event.
e fan blades are modeled as pretwisted thin shallow shells. Chen [29] carried out numerical simulation and experimental research on the rotor dynamics of the blade tip rubbing, in which the frictional force of the blade and the casing adopts the linear rubbing model. Ma et al. [30] studied a rotor system with fixed-point rubbing fault. ey established a point-point contact model to simulate the rotor-stator rubbing. e augmented Lagrangian method is applied to deal with contact constraint conditions. Ma et al. [31] investigated blade disk dovetail structure under blade tip rubbing based on ANSYS software and discussed the effects of rubbing under different rotating speeds and penetration depths. e local rubbing phenomenon was modeled by a pulse force. e conditions of no rubbing, single-blade rubbing, and four-blade rubbing are considered in the numerical simulations. Yang et al. [32] analyzed the vibration of a dual-rotor-bearing-double casing system with pedestal looseness and multistage turbine blade-casing rub considering the distribution of soft and hard coatings.
Since the clearance of the labyrinth seal and the tip seal is very small, when the blade rotor is excessively vibrating, rubbing may occur at both two kinds of seals. Coupled rubbing may result in a more complex dynamic response.
us, it is necessary to investigate the dynamic response of the blade rotor with coupled rubbing. Furthermore, the existing research mainly adopted numerical simulation and experimental methods. Few studies focused on the stability and bifurcation of periodic motion.
In this paper, we propose a study on the coupled rubbing response of the blade rotor and the bifurcation behavior of the rubbing periodic solutions. e structure of this paper is as follows. In Section 2, dynamic equations of coupling rubbing of the blade rotor are established. e complex dynamic behaviors appearing at the rotor's rising up and slowing down are analyzed by numerical simulation. In Section 3, a new continuation shooting method is introduced. And it is employed to solve the periodic response of the coupled rubbing and analyze the stability of the solutions.
e influence of rotation velocity on the dynamic response of the coupled rubbing is analyzed. en, the paths of entering and exiting chaos of the system are discussed. Finally, the conclusion of this paper is summarized.

Numerical Simulation of the Blade Rotor
System with Coupled Rubbing Figure 1 is the model of the blade rotor with coupled rubbing. e rotor is rigidly supported, and the mass of the rotor is concentrated at the disks. e shaft of the rotor is a mass less Euler-Bernoulli beam, and the bending stiffness is EI. m 1 and m 2 are the lumped masses of the labyrinth disk and the blade disk, respectively. Since the diameter of the labyrinth disk is small, the moment of inertia can be negligible; thus, J p1 � J t1 � 0. J t2 and J p2 are the equivalent moment of inertia and the polar moment of inertia of the blade disk, respectively. e equation of motion of the system can be given as

Dynamical Equations of the Blade Rotor.
where x 1 and y 1 are the lateral displacements of the labyrinth disk, and x 2 , y 2 , θ x2 , and θ y2 are the lateral displacements and the deflection angles of the blade disk, respectively. We assume that the eccentricity of the rotor only exists at the blade disk. e 2 is the eccentricity, ω � 2πn/60 is the angular velocity of the rotor, and n is the rotation velocity. k i (i � 1, 2, . . . , 6) are the stiffness coefficients, and their expressions can be given as e damping is supposed to be proportional to the stiffness; there are c i � βk i (i � 1, 2, . . . , 6), where β is the proportional damping coefficient. F x1 and F y1 are horizontal and vertical components of rubbing forces of the labyrinth seal, respectively. And their expressions can be given as follows: where the contact force is given by F Nl � (r 1 − r 10 /α 1 ) 3/2 Hertz contact theory (see Figure 2), where , and μ 1s are Young's modulus and the Poisson ratio of the rotor and the seal, respectively, r 1c is the radius of the labyrinth seal, and r 10 is the clearance between the seal and disk, and the friction force is F tl � F Nl μ 1 .

e Model of Tip Rubbing Forces.
For the tip rubbing, we consider two aspects of deformation: first is the contact deformation of the tip seal Δ H , which is caused by the contact of blade tip and casing and can be described by the Hertz contact model, and the second is the radial deformation of the blade Δ B , which is caused by the bending of the blade (shown in Figure 3(a)) under the action of rubbing forces. e blade is regarded as an Euler beam, and F N is the contact force of a single blade. e friction force is set to obey the Coulomb friction law, and there is F f � F N μ. l b and R d are the span length of the blade and the radius of the disk, respectively. s is the span coordinate of the blade. e bending moment equation of the blade can be given as where w is the deflection of the blade. ρ, A, and E b I b are the density, cross-sectional area, and bending stiffness of the blade, respectively. Considering the boundary conditions w(0) � 0 and w ′ (0) � 0, the deflection function is assumed to be a high-order polynomial as where a i (i � 2 . . . N) are undetermined coefficients. Substituting (5) into (3), the expression of w can be obtained. Figure 4 shows the deflection curves of the beam with N � 3 ∼ 6 . . e calculation results of N � 4 ∼ 6 are essentially the same. Hence, we set N � 4. e radial deformation of the blade tip can be obtained as And based on Hertz contact theory, the deformation of the seal can be given as where α 2 � 0.8 9 16 where E b , E 2s , μ b , and μ 2s are Young's modulus and the Poisson ratio of the blade and seal at the tip seal, respectively, r 2c is the radius of the seal, r 20 is the tip clearance, and there is r 2c � R d + l b + r 20 . We can find the relationship between the total deformation Δ and the contact force F N as Shock and Vibration 3 Hence, the contact force can be written as Figure 3(b) shows the diagrammatic sketch of tip rubbing of the blade rotor. is paper ignores the influence of the initial position angle of the blade. Assuming that the position angle of the first blade is 0, then the coordinates of the first blade tip are where ϕ i � ωt + 2(i − 1)π/N b is the position angle of i th blade and N b is the number of blades. en, we can find the radius of the tip of i th blade as If r ti > r c , the blade will contact with the seal. Considering (9), the rubbing forces of ith blade can be written as and if r ti ≤ r c , Projecting the rubbing forces of ith blade in the x and y directions, then there are   Shock and Vibration By summing up the rubbing forces of all blades, the tip rubbing forces of the rotor can be given as

Numerical Simulation of the Rubbing Response of Blade
Rotor System. By setting (1) can be written in a nondimensional form. en, the Runge-Kutta method is used to obtain the dynamic response. e parameters of the system are given as follows (Table 1). Figure 5 shows the bifurcation diagrams and Lyapunov exponents of the rubbing response, where Figures 5(a) and 5(b) are the bifurcation diagrams considering only the labyrinth seal rubbing and the tip rubbing, respectively. e rubbing response of the labyrinth seal is dominated by the period 1 motion. e tip rubbing response is dominated by a chaotic response. And the coupled rubbing response includes both chaotic response and complex periodic response. At the rotor's rising up, the range of rotation velocity of coupled rubbing, in which contact occurs, is greater than single rubbing fault.
Furthermore, for the coupled rubbing, the contact of the labyrinth seal can maintain at the whole range of rotation velocity of rubbing, and the contact of tip seal ends at n � 5150rpm. Figure 6 shows the partial enlargement of the bifurcation diagram and Lyapunov exponents of coupled rubbing. Figure 7 shows the typical periodic coupled rubbing responses, and Figure 8 gives the aperiodic coupled rubbing responses. It can be seen that the response of the system includes period 1-8, period 12, and period 16 motions, almost periodic motions, and chaotic motions. e transformation between various types of response is also very complicated. us, it is necessary to study the change of response with rotation velocity, especially the path for the periodic response to chaos.

Stability Analysis of the Periodic Solutions of the Coupled Rubbing
In this section, the stability of the coupled rubbing response of the blade rotor is studied by the continuation shooting method. e principle of the continuation shooting method is shown in Figure 9. From a regular point of the solution curve on the Poincare map, we predict the next point in the tangent direction of the solution curve and then make the correction by the Newton-Raphson method on the hyperplane, which is orthogonal to the tangent direction, to get the next regular point of the solution curve. e calculation process is as follows. Consider the following nonlinear differential equations: It can be rewritten in the form of nonlinear algebraic equations as follows: where y � (λ, u 0 ), Σ is the Poincare map, u 0 is the intersection of the periodic solution and Poincare map corresponding to parameter λ, and T is the period of the system. Deriving y i � (λ (i) , u 0(i) ) as a regular point of the solution curve, the next prediction point is where δ i > 0 is the step. v i ∈ R n+1 is the direction of prediction, and it can be determined by where H ′ (y i ) � [zH((λ (i) , u 0(i) )/zλ), zH((λ (i) , u 0(i) )/zu 0 )] and 0 is n dimensional zero vector. However, it is difficult to solve the value of H ′ (y i ). e Secant line can be used instead of the tangent as the prediction direction in practical Parameters Values e masses of the labyrinth disk m 1 2 kg e masses of the blade disk m 1 8 kg e moment of inertia of the blade disk J t2 0.07 kg·m 2 e polar moment of inertia of the blade disk J p2 0.14 kg·m 2 e clearance of the labyrinth seal    Shock and Vibration  calculations. Assuming (λ (i) , u 0(i) ) and (λ (i−1) , u 0(i−1) ) are two adjacent regular points, then e correction process is carried out by the following formula: and H ′ (z k ) � [zH(λ p(k) , u p0(k) )/zλ, zH(λ p(k) , u p0(k) )/zu 0 ]. zH(λ p(k) , u p0(k) )/zu 0 can be solved by the following initial value problem: When t � T, zH(λ p(k) , u p0(k) )/zu 0 � z(T). zH(λ p(k) , u p0(k) )/zλ can be given as where Δ is a small perturbation. Iterate from k � 1until it is satisfied ‖z k+1 − z k ‖ ≤ ε, and set y i+1 � z k+1 . For this method, zH(z k )/zu 0 obtained in the last iteration is the Floquet multiplier matrix of the periodic solution, and the stability of the periodic solution can be determined by its eigenvalues.
Combining the numerical simulation with the results of periodic response calculation and the stability analysis, Figure 10 shows the bifurcation diagram of the periodic response of blade rotor with coupled rubbing. e response is period 1 motion at the beginning of rubbing. As the rotation velocity increases, period-doubling bifurcation occurs at point a and the response becomes period 2 motion. en, the response enters almost periodic motion through Hopf bifurcation at point b. e response breaks into chaos through almost periodic torus. e response returns to the period 2 motion at the c point through Hopf bifurcation. Further, the period 2 motion enters chaos through the period-doubling bifurcation near the point d, and the chaotic motion changes into the periodic motion at about n � 5033rpm. en, the response goes back to period 2 motion (shown in Figures 10(a) and 10(b)). e period 2 motion enters chaos again through the Hopf bifurcation at point e and then leaves chaos by the saddle node bifurcation at point f. e response of the system becomes period 6 motion, and the period 6 motion jumps into period 2 motion through a saddle node bifurcation at point g. en, period 2 motion enters chaos through the saddle node bifurcation at point h (shown in Figure 10(c)). e chaotic response of the system changes into period 6 motion by an inverse period-doubling bifurcation near point i and then changes to period 2 motion by saddle node bifurcation at point j. Period 2 motion becomes period 1 motion again through saddle node bifurcation at point k (shown in Figure 10(d)). Point l and m are period-doubling bifurcation points, and the response between the two points is period 2 motion (shown in Figure 10(e)). With further increase in rotating velocity, the motion of period 1 enters chaos by saddle node bifurcation at point n. e chaotic response turns into period 6 motion through saddle node bifurcation at point o, and then period 6 motion becomes unstable through period-doubling bifurcation at p point and becomes stable again through period-doubling bifurcation at point q. After period-doubling inverse bifurcation at point r, the response of the system is period 3 motion (shown in Figures 10(f ) and 10(g)). Period 3 motion enters chaos through intermittent at point s, then exits chaos at point t through Hopf bifurcation, and enters chaos again through intermittent at point u (shown in Figure 10(h)). e response of the system moves into period 5 motion at the point v through the saddle node bifurcation. e period 5 motion enters chaos through period-doubling bifurcation at w point, exits chaos through period-doubling bifurcation at point x, and then enters chaos again through Hopf bifurcation at y point. It is worth noting that the path to chaos here is not almost periodic torus. Figure 11 is the Poincare map of the chaotic response at n � 5670.8rpm. As the fixed points of period 5 motion become unstable, almost periodic torus disappears, and finally the response is  With further increase in rotation velocity, the chaotic response returns to period 5 motion through the perioddoubling bifurcation (near z point) and then enters chaos through the saddle node bifurcation at point a 1 .
e response of the system becomes period 1 motion through the Hopf bifurcation point b 1 and then jumps to period 3 motion through the saddle node bifurcation at the point c 1 . e motion of period 3 enters chaos at the point d 1 through the period-doubling bifurcation, then exits chaos through the period-doubling bifurcation at point e 1 , and enters chaos again through the saddle node bifurcation at point f 1 . At the point g 1 , the response of the system returns to period 1 motion through the Hopf bifurcation and then jumps through saddle node bifurcation at point h 1 , and the rotor and the seal are out of contact.
When the rotor slows down, the response of the system jumps up near n � 5300rpm and rubbing occurs. e following changes of the response are mostly the same as the process of rotor's rising up, and the difference is in the range of rotational velocity of 5100-5140 rpm. When the rotation velocity is reduced, the response changes into period 2 motion through a period-doubling bifurcation at the point i 1 . en, saddle node bifurcation occurs at the point j 1 , the period 2 response jumps to the upper and lower solution branches, and then the periodic response becomes unstable through the Hopf bifurcation at the point k 1 .
In addition, it is shown that there are coexistence of multiple solutions corresponding to some rotation velocities, such as the coexistence of periodic solutions (see Figures 10(d) and 10(k)), coexistence of periodic solution and almost periodic solution (see Figure 10(j)), and coexistence of periodic solution and chaos (see Figure 10(i)). If the system is disturbed in the process of work, the response may change into other motions.

Conclusions
In this paper, the dynamic response and its stability of a blade rotor with coupled rubbing of the labyrinth seal and the tip seal are investigated. Dynamic equations of the blade rotor with coupled rubbing are established. e rubbing force of the labyrinth seal adopts the Hertz contact model. e tip rubbing model considers the Hertz contact deformation of the tip seal and the bending deformation of the blade. e bending deformation of the blade is described more accurately by using a high-order polynomial. e dynamic equations of the blade rotor with coupled rubbing are numerically simulated. It is found that the labyrinth seal rubbing response is dominated by periodic motion, the tip rubbing response is dominated by chaos, and the coupled rubbing response includes both chaotic motion and complex periodic motions. Moreover, the range of rotational speed at which the coupled rubbing rotor contacts is greater than a single rubbing fault.
A new continuation shooting method is employed to calculate the periodic solutions and analyze their stability. Considering the results of numerical simulation, the bifurcation diagrams of the periodic response of coupled rubbing are given, and the paths of the system entering and exiting chaos are analyzed. It shows that the paths of the system entering chaos include period-doubling bifurcation and intermittent, and the paths exiting chaos include Hopf bifurcation, intermittent, and period-doubling bifurcation. In addition, the investigation finds the coexistence of multiple solutions corresponding to some rotational velocities, such as coexistence of periodic solutions, coexistence of periodic solution and almost periodic solution, and coexistence of periodic solution and chaos.
Data Availability e data that support the findings of this study are available from the corresponding author (Huabiao Zhang) upon reasonable request.

Conflicts of Interest
e author(s) declare that there are no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.  Shock and Vibration 13