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

Authors: Huabiao Zhang; Xinye Li; Dongai Wang; Tingting Liu
Journal: Shock and Vibration
 (2021)

---

## Abstract

The dynamic response and its stability of a blade rotor with coupled rubbing in the labyrinth seal and tip seal are investigated. The 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. The paths of the system for entering and exiting chaos are analyzed.


---

## Body

## 1. 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].There are usually two kinds of seals for turbo machines: labyrinth seals and tip seals. The rubbing of labyrinth seals is usually described as a rubbing problem between two cylindrical surfaces. The contact forces may be expressed as a linear or nonlinear function of the deformation of interaction. The friction forces are given by the Coulomb law of friction. Many researchers focused on related problems, such as forward periodic motion [2–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 dual-rotor-bearing-coupling misalignment system with blade-casing 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 blade-machine 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. The 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. They established a point-point contact model to simulate the rotor-stator rubbing. The 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. The local rubbing phenomenon was modeled by a pulse force. The 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. Thus, 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. The structure of this paper is as follows. In Section2, dynamic equations of coupling rubbing of the blade rotor are established. The 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. The influence of rotation velocity on the dynamic response of the coupled rubbing is analyzed. Then, the paths of entering and exiting chaos of the system are discussed. Finally, the conclusion of this paper is summarized.


## 2. Numerical Simulation of the Blade Rotor System with Coupled Rubbing

### 2.1. Dynamical Equations of the Blade Rotor

Figure1 is the model of the blade rotor with coupled rubbing. The rotor is rigidly supported, and the mass of the rotor is concentrated at the disks. The shaft of the rotor is a mass less Euler–Bernoulli beam, and the bending stiffness is EI. m1 and m2 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, Jp1=Jt1=0. Jt2 and Jp2 are the equivalent moment of inertia and the polar moment of inertia of the blade disk, respectively. The equation of motion of the system can be given as(1)m1x¨1+k1x1+k4x2+k5θy2+c1x˙1+c4x˙2+c5θ˙y2+Fx1=0,m1y¨1+k1y1+k4y2+k5θx2+c1y˙1+c4y˙2+c5θ˙x2+Fy1=−m1g,m2x¨2+k4x1+k2x2+k6θy2+c4x˙1+c2x˙2+c6θ˙y2+Fx2=m2e2ω2cosωt,m2y¨2+k4y1+k2y2+k6θx2+c4y˙1+c2y˙2+c6θ˙x2+Fy2=m2e2ω2sinωt−m1g,Jt2θ¨x2+Jp2ωθ˙y2+k5y2+k6yx2+k3θx2+c5y˙1+c6y˙2+c3θ˙x2=0,Jt2θ¨y2+Jp2ωθ˙x2+k5x2+k6xy2+k3θy2+c5x˙1+c6x˙2+c3θ˙y2=0,where x1andy1 are the lateral displacements of the labyrinth disk, and x2,y2,θ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. e2 is the eccentricity, ω=2πn/60 is the angular velocity of the rotor, and n is the rotation velocity. kii=1,2,…,6 are the stiffness coefficients, and their expressions can be given as(2)k1=12EIn1+133+4n1n12l3,k2=3EI4n1n23+12n23+4n1+3n23l33+4n1,k3=3EI4n1n2+4n1+4n2+3n2l3+4n1,k4=−6EI2n12+6n1+3n13+4n1l3,k5=6EIn1+12n1+1n13+4n1l3,k6=−3EI4n1n22+6n22−4n1−3n22l33+4n1.Figure 1
Model of the blade rotor.The damping is supposed to be proportional to the stiffness; there areci=βkii=1,2,…,6, where β is the proportional damping coefficient. Fx1andFy1 are horizontal and vertical components of rubbing forces of the labyrinth seal, respectively. And their expressions can be given as follows:(3)Fx1=FNlx1−Ftly1r1,Fy1=Ftlx1+FNly1r1,r1>r10,Fx1=Fy1=0,r1⩽r10,where the contact force is given by FNl=r1−r10/α13/2 Hertz contact theory (see Figure 2), where α1=0.89/161−μ1s2/E1s+1−μ12/E12r10/r1c21/3, r1=x12+y12, E1,E1s,μ1,andμ1s are Young's modulus and the Poisson ratio of the rotor and the seal, respectively, r1c is the radius of the labyrinth seal, and r10 is the clearance between the seal and disk, and the friction force is Ftl=FNlμ1.Figure 2
Rubbing forces of the labyrinth seal rubbing.


### 2.2. The 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. The blade is regarded as an Euler beam, and FN is the contact force of a single blade. The friction force is set to obey the Coulomb friction law, and there is Ff=FNμ. lbandRd are the span length of the blade and the radius of the disk, respectively. s is the span coordinate of the blade. The bending moment equation of the blade can be given as(4)d2wds2=FNwlb−w+FNμ2lb−s−∫slbρARd+δω2wδ−wdδEbIb,where w is the deflection of the blade. ρ,A,andEbIb are the density, cross-sectional area, and bending stiffness of the blade, respectively. Considering the boundary conditions w0=0andw′0=0, the deflection function is assumed to be a high-order polynomial as(5)w=∑i=2Naisi,where aii=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.. The calculation results of N=4∼6 are essentially the same. Hence, we set N=4. The radial deformation of the blade tip can be obtained as(6)ΔB=∫0lbdwds2ds.Figure 3
The model of tip rubbing of the blade rotor: (a) rubbing forces of a single blade and (b) diagrammatic sketch of tip rubbing of the blade rotor.
(a)(b)Figure 4
The deflection curve of the bladeFN=5000Nandn=6000rpm.And based on Hertz contact theory, the deformation of the seal can be given as(7)ΔH=α2FN2/3,where(8)α2=0.89161−μ2s2E2s+1−μb2Eb2r20r2c21/3,where Eb,E2s,μb,andμ2s are Young's modulus and the Poisson ratio of the blade and seal at the tip seal, respectively, r2c is the radius of the seal, r20 is the tip clearance, and there is r2c=Rd+lb+r20. We can find the relationship between the total deformation Δ and the contact force FN as(9)Δ=ΔH+ΔB=fFN.Hence, the contact force can be written as(10)FN=f−1Δ=f¯Δ.Figure3(b) shows the diagrammatic sketch of tip rubbing of the blade rotor. This 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(11)xti=x2+Rd+lbcosϕi,yti=y2+Rd+lbsinϕi,where ϕi=ωt+2i−1π/Nb is the position angle of i th blade and Nb is the number of blades. Then, we can find the radius of the tip of i th blade as(12)rti=x2+Rd+lbcosϕi2+y2+Rd+lbsinϕi2.Ifrti>rc, the blade will contact with the seal. Considering (9), the rubbing forces of ith blade can be written as(13)FNi=f¯rti−rc,Ffi=μFNi,and if rti≤rc,(14)FNi=0,Ffi=0.Projecting the rubbing forces ofith blade in the x and y directions, then there are(15)Fxi=FNicosϕi−Ffisinϕi,Fyi=FNisinϕi+Fficosϕi.By summing up the rubbing forces of all blades, the tip rubbing forces of the rotor can be given as(16)Fx2=∑i=1NbFxi,Fy2=∑i=1NbFyi.


### 2.3. Numerical Simulation of the Rubbing Response of Blade Rotor System

By settingx1=X1r20,y1=Y1r20,x2=X2r20,y2=Y2r20,θx2=ψx2r20/l,θy2=ψy2r20/l, and τ=ωt, (1) can be written in a nondimensional form. Then, the Runge–Kutta method is used to obtain the dynamic response. The parameters of the system are given as follows (Table 1).Table 1
Parameters of the rotor system.
ParametersValuesThe masses of the labyrinth diskm12 kgThe masses of the blade diskm18 kgThe moment of inertia of the blade diskJt20.07 kg·m2The polar moment of inertia of the blade diskJp20.14 kg·m2The clearance of the labyrinth sealr100.1 mmThe clearance of the tip sealr200.1 mmThe proportional damping coefficient isβ0.0001 N/mThe length of the shaft between the labyrinth disk and the blade diskl60 mmYoung’s modulus of the labyrinth diskE12 × 1011 PaYoung’s modulus of the bladeEb2 × 1011 PaYoung’s modulus of the labyrinth sealE1s2 × 1011 PaYoung’s modulus of the tip sealE2s2 × 1011 PaThe diameter of the shaftd40 mmThe ratio of the length of the shaftn16The ratio of the length of the shaftn25.5The span length of the bladelb45 mmThe radius of the diskRd73 mmThe density of the bladeρ8000 kg/m3The cross-sectional areaA64 mm2Poisson ratio of the labyrinth diskμ10.1Poisson ratio of the blade diskμ20.1The bending stiffness of the bladeEbIb21 N/mThe number of the bladesNb20The eccentricity of the blade diske20.008 mmThe radius of the labyrinth sealr1c30 mmAcceleration of gravityg9.8 m/sFigure5 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. The rubbing response of the labyrinth seal is dominated by the period 1 motion. The 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.Figure 5
Bifurcation diagrams and Lyapunov exponents of the rotor with rubbing: (a) rubbing of labyrinth seal; (b) tip rubbing; (c) coupled rubbing at the rotor’s running up; (d) coupled rubbing at the rotor’s slowing up. RC is the range of rotation velocity of labyrinth rubbing; RB is the range of rotation velocity of tip rubbing; LE is the Lyapunov exponents.
(a)(b)(c)(d)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 atn=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. The transformation between various types of response is also very complicated. Thus, it is necessary to study the change of response with rotation velocity, especially the path for the periodic response to chaos.Figure 6
Enlargements of the bifurcation diagram and Lyapunov exponents of blade rotor with coupled rubbing: (a–c) at the rotor rising up; (d) at the rotor slowing down. LE is the Lyapunov exponents.
(a)(b)(c)(d)Figure 7
Orbits and poincare maps of periodic responses of the blade rotor with coupled rubbing: (a) Period 1 motion without rubbingn=4900rpm; (b) Period 1 motion with rubbing n=5013rpm; (c) Period 2 motion n=5017rpm; (d) Period 24 motion n=5035rpm; (e) Period 4 motion n=5090rpm; (f) Period 12 motion n=5104rpm; (g) Period 6 motion n=5108rpm; (h) Period 8 motion n=5510rpm; (i) Period 16 motion n=5528rpm; (j) Period 7 motion n=5534rpm; (k) Period 3 motion n=5608rpm; (l) Period 5 motion n=5640rpm.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 8
Orbits and Poincare maps of aperiodic responses of blade rotor with coupled rubbing: (a)n=5068rpm, (b) n=5603rpm, (c) n=5930rpm, (d) n=5016rpm, (e) n=5032rpm, (f) n=5078rpm, (g) n=5500rpm, (h) n=5612rpm, where (a)∼(c) are almost periodic motions and (d)∼(h) are chaos.
(a)(b)(c)(d)(e)(f)(g)(h)


## 2.1. Dynamical Equations of the Blade Rotor

Figure1 is the model of the blade rotor with coupled rubbing. The rotor is rigidly supported, and the mass of the rotor is concentrated at the disks. The shaft of the rotor is a mass less Euler–Bernoulli beam, and the bending stiffness is EI. m1 and m2 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, Jp1=Jt1=0. Jt2 and Jp2 are the equivalent moment of inertia and the polar moment of inertia of the blade disk, respectively. The equation of motion of the system can be given as(1)m1x¨1+k1x1+k4x2+k5θy2+c1x˙1+c4x˙2+c5θ˙y2+Fx1=0,m1y¨1+k1y1+k4y2+k5θx2+c1y˙1+c4y˙2+c5θ˙x2+Fy1=−m1g,m2x¨2+k4x1+k2x2+k6θy2+c4x˙1+c2x˙2+c6θ˙y2+Fx2=m2e2ω2cosωt,m2y¨2+k4y1+k2y2+k6θx2+c4y˙1+c2y˙2+c6θ˙x2+Fy2=m2e2ω2sinωt−m1g,Jt2θ¨x2+Jp2ωθ˙y2+k5y2+k6yx2+k3θx2+c5y˙1+c6y˙2+c3θ˙x2=0,Jt2θ¨y2+Jp2ωθ˙x2+k5x2+k6xy2+k3θy2+c5x˙1+c6x˙2+c3θ˙y2=0,where x1andy1 are the lateral displacements of the labyrinth disk, and x2,y2,θ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. e2 is the eccentricity, ω=2πn/60 is the angular velocity of the rotor, and n is the rotation velocity. kii=1,2,…,6 are the stiffness coefficients, and their expressions can be given as(2)k1=12EIn1+133+4n1n12l3,k2=3EI4n1n23+12n23+4n1+3n23l33+4n1,k3=3EI4n1n2+4n1+4n2+3n2l3+4n1,k4=−6EI2n12+6n1+3n13+4n1l3,k5=6EIn1+12n1+1n13+4n1l3,k6=−3EI4n1n22+6n22−4n1−3n22l33+4n1.Figure 1
Model of the blade rotor.The damping is supposed to be proportional to the stiffness; there areci=βkii=1,2,…,6, where β is the proportional damping coefficient. Fx1andFy1 are horizontal and vertical components of rubbing forces of the labyrinth seal, respectively. And their expressions can be given as follows:(3)Fx1=FNlx1−Ftly1r1,Fy1=Ftlx1+FNly1r1,r1>r10,Fx1=Fy1=0,r1⩽r10,where the contact force is given by FNl=r1−r10/α13/2 Hertz contact theory (see Figure 2), where α1=0.89/161−μ1s2/E1s+1−μ12/E12r10/r1c21/3, r1=x12+y12, E1,E1s,μ1,andμ1s are Young's modulus and the Poisson ratio of the rotor and the seal, respectively, r1c is the radius of the labyrinth seal, and r10 is the clearance between the seal and disk, and the friction force is Ftl=FNlμ1.Figure 2
Rubbing forces of the labyrinth seal rubbing.


## 2.2. The 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. The blade is regarded as an Euler beam, and FN is the contact force of a single blade. The friction force is set to obey the Coulomb friction law, and there is Ff=FNμ. lbandRd are the span length of the blade and the radius of the disk, respectively. s is the span coordinate of the blade. The bending moment equation of the blade can be given as(4)d2wds2=FNwlb−w+FNμ2lb−s−∫slbρARd+δω2wδ−wdδEbIb,where w is the deflection of the blade. ρ,A,andEbIb are the density, cross-sectional area, and bending stiffness of the blade, respectively. Considering the boundary conditions w0=0andw′0=0, the deflection function is assumed to be a high-order polynomial as(5)w=∑i=2Naisi,where aii=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.. The calculation results of N=4∼6 are essentially the same. Hence, we set N=4. The radial deformation of the blade tip can be obtained as(6)ΔB=∫0lbdwds2ds.Figure 3
The model of tip rubbing of the blade rotor: (a) rubbing forces of a single blade and (b) diagrammatic sketch of tip rubbing of the blade rotor.
(a)(b)Figure 4
The deflection curve of the bladeFN=5000Nandn=6000rpm.And based on Hertz contact theory, the deformation of the seal can be given as(7)ΔH=α2FN2/3,where(8)α2=0.89161−μ2s2E2s+1−μb2Eb2r20r2c21/3,where Eb,E2s,μb,andμ2s are Young's modulus and the Poisson ratio of the blade and seal at the tip seal, respectively, r2c is the radius of the seal, r20 is the tip clearance, and there is r2c=Rd+lb+r20. We can find the relationship between the total deformation Δ and the contact force FN as(9)Δ=ΔH+ΔB=fFN.Hence, the contact force can be written as(10)FN=f−1Δ=f¯Δ.Figure3(b) shows the diagrammatic sketch of tip rubbing of the blade rotor. This 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(11)xti=x2+Rd+lbcosϕi,yti=y2+Rd+lbsinϕi,where ϕi=ωt+2i−1π/Nb is the position angle of i th blade and Nb is the number of blades. Then, we can find the radius of the tip of i th blade as(12)rti=x2+Rd+lbcosϕi2+y2+Rd+lbsinϕi2.Ifrti>rc, the blade will contact with the seal. Considering (9), the rubbing forces of ith blade can be written as(13)FNi=f¯rti−rc,Ffi=μFNi,and if rti≤rc,(14)FNi=0,Ffi=0.Projecting the rubbing forces ofith blade in the x and y directions, then there are(15)Fxi=FNicosϕi−Ffisinϕi,Fyi=FNisinϕi+Fficosϕi.By summing up the rubbing forces of all blades, the tip rubbing forces of the rotor can be given as(16)Fx2=∑i=1NbFxi,Fy2=∑i=1NbFyi.


## 2.3. Numerical Simulation of the Rubbing Response of Blade Rotor System

By settingx1=X1r20,y1=Y1r20,x2=X2r20,y2=Y2r20,θx2=ψx2r20/l,θy2=ψy2r20/l, and τ=ωt, (1) can be written in a nondimensional form. Then, the Runge–Kutta method is used to obtain the dynamic response. The parameters of the system are given as follows (Table 1).Table 1
Parameters of the rotor system.
ParametersValuesThe masses of the labyrinth diskm12 kgThe masses of the blade diskm18 kgThe moment of inertia of the blade diskJt20.07 kg·m2The polar moment of inertia of the blade diskJp20.14 kg·m2The clearance of the labyrinth sealr100.1 mmThe clearance of the tip sealr200.1 mmThe proportional damping coefficient isβ0.0001 N/mThe length of the shaft between the labyrinth disk and the blade diskl60 mmYoung’s modulus of the labyrinth diskE12 × 1011 PaYoung’s modulus of the bladeEb2 × 1011 PaYoung’s modulus of the labyrinth sealE1s2 × 1011 PaYoung’s modulus of the tip sealE2s2 × 1011 PaThe diameter of the shaftd40 mmThe ratio of the length of the shaftn16The ratio of the length of the shaftn25.5The span length of the bladelb45 mmThe radius of the diskRd73 mmThe density of the bladeρ8000 kg/m3The cross-sectional areaA64 mm2Poisson ratio of the labyrinth diskμ10.1Poisson ratio of the blade diskμ20.1The bending stiffness of the bladeEbIb21 N/mThe number of the bladesNb20The eccentricity of the blade diske20.008 mmThe radius of the labyrinth sealr1c30 mmAcceleration of gravityg9.8 m/sFigure5 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. The rubbing response of the labyrinth seal is dominated by the period 1 motion. The 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.Figure 5
Bifurcation diagrams and Lyapunov exponents of the rotor with rubbing: (a) rubbing of labyrinth seal; (b) tip rubbing; (c) coupled rubbing at the rotor’s running up; (d) coupled rubbing at the rotor’s slowing up. RC is the range of rotation velocity of labyrinth rubbing; RB is the range of rotation velocity of tip rubbing; LE is the Lyapunov exponents.
(a)(b)(c)(d)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 atn=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. The transformation between various types of response is also very complicated. Thus, it is necessary to study the change of response with rotation velocity, especially the path for the periodic response to chaos.Figure 6
Enlargements of the bifurcation diagram and Lyapunov exponents of blade rotor with coupled rubbing: (a–c) at the rotor rising up; (d) at the rotor slowing down. LE is the Lyapunov exponents.
(a)(b)(c)(d)Figure 7
Orbits and poincare maps of periodic responses of the blade rotor with coupled rubbing: (a) Period 1 motion without rubbingn=4900rpm; (b) Period 1 motion with rubbing n=5013rpm; (c) Period 2 motion n=5017rpm; (d) Period 24 motion n=5035rpm; (e) Period 4 motion n=5090rpm; (f) Period 12 motion n=5104rpm; (g) Period 6 motion n=5108rpm; (h) Period 8 motion n=5510rpm; (i) Period 16 motion n=5528rpm; (j) Period 7 motion n=5534rpm; (k) Period 3 motion n=5608rpm; (l) Period 5 motion n=5640rpm.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 8
Orbits and Poincare maps of aperiodic responses of blade rotor with coupled rubbing: (a)n=5068rpm, (b) n=5603rpm, (c) n=5930rpm, (d) n=5016rpm, (e) n=5032rpm, (f) n=5078rpm, (g) n=5500rpm, (h) n=5612rpm, where (a)∼(c) are almost periodic motions and (d)∼(h) are chaos.
(a)(b)(c)(d)(e)(f)(g)(h)


## 3. 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. The principle of the continuation shooting method is shown in Figure9. 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. The calculation process is as follows.Figure 9
Principle of the continuation shooting method.Consider the following nonlinear differential equations:(17)u˙=fλ,u,t,u∈Rn,t>0.It can be rewritten in the form of nonlinear algebraic equations as follows:(18)Hy=Hλ,u0=uλ,u0,T−u0=0,u0∈Σ,where y=λ,u0, Σ is the Poincare map, u0 is the intersection of the periodic solution and Poincare map corresponding to parameter λ, and T is the period of the system. Deriving yi=λi,u0i as a regular point of the solution curve, the next prediction point is(19)z0=yi+δivi,where δi>0 is the step. vi∈Rn+1 is the direction of prediction, and it can be determined by(20)H′yiviTvi=0¯1vi−1Tvi>0,where H′yi=∂Hλi,u0i/∂λ,∂Hλi,u0i/∂u0 and 0¯ is n dimensional zero vector. However, it is difficult to solve the value of H′yi. The Secant line can be used instead of the tangent as the prediction direction in practical calculations. Assuming λi,u0i and λi−1,u0i−1 are two adjacent regular points, then(21)vi=λi−λi−1,u0i−u0i−1λi−λi−1,u0i−u0i−1.The correction process is carried out by the following formula:(22)wk=H′zkviT−1−Hzk0,zk+1=zk+wk,where zk=λpk,up0k and H′zk=∂Hλpk,up0k/∂λ,∂Hλpk,up0k/∂u0. ∂Hλpk,up0k/∂u0 can be solved by the following initial value problem:(23)z˙=Dfλ,u,tz,z∈Rn×n,z0=I.Whent=T, ∂Hλpk,up0k/∂u0=zT. ∂Hλpk,up0k/∂λ can be given as(24)∂Hλpk,up0k∂λ≈Hλpk+Δ¯,up0k−Hλpk,up0kΔ¯=uλpk+Δ¯,up0k,T−uλpk,up0k,TΔ¯,where Δ¯ is a small perturbation. Iterate from k=1until it is satisfied zk+1−zk≤ε, and set yi+1=zk+1. For this method, ∂Hzk/∂u0 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, Figure10 shows the bifurcation diagram of the periodic response of blade rotor with coupled rubbing. The 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. Then, the response enters almost periodic motion through Hopf bifurcation at point b. The response breaks into chaos through almost periodic torus. The 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. Then, the response goes back to period 2 motion (shown in Figures 10(a) and 10(b)). The 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. The 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. Then, period 2 motion enters chaos through the saddle node bifurcation at point h (shown in Figure 10(c)).Figure 10
Bifurcation diagrams of the periodic response of coupled rubbing of the blade rotor: (a–k) the bifurcation diagrams with the rotor rising up; (l) the bifurcation diagram with the rotor slowing down. US means the unstable solution, SP means the static bifurcation point, HP means the Hopf bifurcation point, and DP means the periodic doubling bifurcation point.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)The chaotic response of the system changes into period 6 motion by an inverse period-doubling bifurcation near pointi 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. The 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)).The response of the system moves into period 5 motion at the pointv through the saddle node bifurcation. The 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 separated from the position near the fixed points and moves towards the chaotic attractor.Figure 11
Poincare map of the chaotic response atn=5670.8rpm.With further increase in rotation velocity, the chaotic response returns to period 5 motion through the period-doubling bifurcation (nearz point) and then enters chaos through the saddle node bifurcation at point a1. The response of the system becomes period 1 motion through the Hopf bifurcation point b1 and then jumps to period 3 motion through the saddle node bifurcation at the point c1. The motion of period 3 enters chaos at the point d1 through the period-doubling bifurcation, then exits chaos through the period-doubling bifurcation at point e1, and enters chaos again through the saddle node bifurcation at point f1. At the point g1, the response of the system returns to period 1 motion through the Hopf bifurcation and then jumps through saddle node bifurcation at point h1, and the rotor and the seal are out of contact.When the rotor slows down, the response of the system jumps up nearn=5300rpm and rubbing occurs. The 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 i1. Then, saddle node bifurcation occurs at the point j1, 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 k1.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 Figures10(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.


## 4. 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. The rubbing force of the labyrinth seal adopts the Hertz contact model. The tip rubbing model considers the Hertz contact deformation of the tip seal and the bending deformation of the blade. The bending deformation of the blade is described more accurately by using a high-order polynomial. The 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.



---

*Source: 3826429-2021-10-23.xml*