Rub-Impact Dynamics of Shrouded Blades under Bending-Torsion Coupling Vibration

Abstract

Shroud devices which are typical cyclic symmetric structures are widely used to reduce the vibration of turbine blades in aero engines. Asymmetric rub-impact of adjacent shrouds or aerodynamic excitation forces can excite the bending-torsion coupling vibration of shrouded blades, which will lead to complex contact motions. The aim of this paper is to study the rub-impact dynamic characteristics of bending-torsion coupling vibration of shrouded blades using a numerical method. The contact-separation transition mechanism under complex motions is studied, the corresponding boundary conditions are set up, and the influence of moments of contact forces and aerodynamic excitation forces on the motion of the blade is considered. A three-degree-of-freedom mass-spring model including two mass blocks with the same size and shape is established to simulate the bending-torsion coupling vibration, and the dynamic equations of shrouded blades under different contact conditions are derived. An algorithm based on the fourth-order Runge–Kutta method is presented for simulations. Variation laws of the forced response characteristics of shrouded blades under different parameters are studied, on the basis of which the method to evaluate the vibration reduction characteristics of the shrouded blade system when the motion of the blade is chaotic is discussed. Then, the vibration reduction law of shrouded blades under bending-torsion coupling vibration is obtained.

1. Introduction

High-cycle fatigue caused by high vibratory stress is one of the major causes of service failures of aero engine turbine blades. The friction constrained structure, which is simple and insensitive to temperature, is effective to reduce the vibration of turbine blades. The shroud device is a typical dry friction structure and has been widely used to reduce the vibration of turbine blades [1,2,3]. The main dynamic research progress on the system is shown as follows, including the dynamic model, the friction contact model, the response characteristics, and the vibration reduction characteristics.

1.1. About the Dynamic Model and Friction Contact Model

Based on one-dimensional motion, Iwan [4] proposed the bilinear hysteresis model, a classical dry friction force model in which contact elastic deformation was considered. Considering the normal motion and tangential complex curve motion of contact bodies, Yang et al. [5] improved the bilinear hysteresis model and examined the jump phenomenon. With the improved fiction model in [5], a 3D friction contact model was employed and the resonant response of the shrouded blade with a three-dimensional shroud constraint was studied in [6]. Petrov [7,8] developed the finite element (FE) model for blades with nonlinear dry friction contact, which is closer to engineering practice, and in [8] a new friction model in which time-varying friction contact parameters could be taken into account was provided. Cigeroglu [9] presented a two-dimensional micro slip friction model between two elastic structures with normal load variation, characterized the stick-slip-separation of the contact interface, and determined the resulting friction force. In [10], a friction contact stiffness model was proposed in which the contact stiffness was calculated using the Hertz contact theory and fractal geometry. Afzal [11] developed an analytical expression to compute the Jacobian matrix for the 3D friction contact modeling of a turbine blade with a shroud contact interface having an arbitrary 3D relative displacement, which drastically reduced the computation time of the Jacobian matrix with respect to the classical finite difference method. Considering the effects of the centrifugal stiffening, spin softening, and Coriolis force, Ma et al. [12] proposed a dynamic model of rotating shrouded blades with impacts and analyzed the response characteristics of shrouded blades. In [13], taking the twisted shrouded blade as a cantilever beam with a tip mass, the analytical model of the twisted shrouded blade with shroud contact was established by using Timoshenko beam theory and Coulomb’s friction law, with which the effects of the stagger angles and the twist angles on the vibration response of the damped blade system were discussed and validated by the finite element method. In [14], a dynamic model for the turbine blade with dry friction damper was set up, in which the convective inertial force along the tangential direction and Coriolis inertial force along the normal direction caused by the disc rotation were considered. Combining an experimental study with a numerical study, Umer [15] highlighted the relevance of an accurate representation of the constraints induced by friction contacts and discussed the adequacy of state-of the-art contact models. In [16], to study the nonlinear vibration caused by impact and friction between adjacent shrouds, a new dynamic model of shrouded blades with elastic support, which can consider the inertia effects of passive blades, the change of shroud position, and the variation of blade root support stiffness, was established and verified by the finite element model.

1.2. About the Response and Vibration Reduction Characteristics

In [17], the routes of periodic impacts to chaos in two kinds of strong resonance cases of a two-degree-of-freedom system were obtained by numerical simulations. In [18], the dynamic characteristics of a shrouded single blade of the shrouded cascade impact experiment rig were studied by finite element simulation and validated by experimental results. In [19], using an enhanced multiterm harmonic balance method, the forced vibration response of a single degree of freedom torsion system with Coulomb friction was studied when the normal load varied periodically. In [20], an approximate approach was proposed to compute the dynamic response of a structure constrained by friction interfaces due to tip-rub, and the nonlinear vibration characteristics of the blade were investigated using this approach in terms of the Poincare graph, the frequency spectrum of the responses, and the amplitude-frequency curves. Sayed et al. [21] explained the decrease of vibration amplitude by changes in boundary conditions caused by stick/slip behavior, while demonstrating the contribution of various energy dissipation and contact state changes to the peak level. In [22], an FEA method was introduced by which an actual modal analysis of the shrouded blade was made, and the dynamic characteristics of the system were discussed. In [23], the vibro-impact mechanism of a rotating shrouded blade with asymmetric gaps was investigated, on the basis of which the influence of the blade shroud gaps on the vibration amplitude of the system was studied. In [24], two types of the impact vibration between adjacent shrouded blades were analyzed, based on which the influence of the parameters on the dynamic behavior of the shrouded blade was studied. Santhosh [25] studied the multivalued frequency response and jump phenomenon of the shrouded blade with rub and impact. In [26], the mechanism of stick-slip transition of multicontacts was studied with Coulomb friction, and the effects of relevant parameters on the response characteristics of shrouded blades were obtained. In [27], a method based on the finite element method was proposed to analyze the damping characteristics of blades with shroud considering the change of equivalent contact stiffness of the shroud.

From the above literature reviews, the influence of contact forces and aerodynamic excitation forces on the translation of shrouded blades are considered, while the influence of moments of these forces on the motion of shrouded blades are not considered. The contact-separation transition boundary conditions under bending-torsion coupling vibration and the corresponding reduction characteristics of shrouded blades should be investigated further. In addition, it is clear that the stick-slip-separation transition may lead to chaos; thus, how to evaluate the vibration reduction characteristics of the shrouded blade system when the steady response of the blade is chaotic needs to be considered.

In this paper, the bending-torsion coupling vibration of shrouded blades caused by the asymmetric rub-impact and aerodynamic excitation forces is considered. The contact-separation transition boundary conditions under bending-torsion coupling vibration are set up, and moments of contact forces and aerodynamic excitation forces on the motion of the system are considered. On this foundation, a 3-DOF mass-spring model including two mass blocks with the same shape and size is established, and the dynamic equations of the system in different contact conditions are derived. A fourth-order Runge–Kutta method algorithm is presented for simulations. Characteristics of the forced response of the system in different conditions are investigated in detail, and the method to evaluate the vibration reduction characteristics of the shrouded blade system when the steady response of the blade is chaotic is studied. On the basis of the above, the vibration reduction characteristics of the bending-torsion coupling system are studied.

2. Dynamic Modeling of the Bending-Torsion Coupling Vibration System

2.1. Structure Model and the Dynamic Model of the System

In engineering, shrouded blades are set up in a circle around the disc, and the whole structure can be considered to have a cyclic symmetry. In this paper, the structure model including two shrouded blades is discussed, as shown in Figure 1, where the orthogonal coordinate system attached to the disc is defined by three orthogonal directions, namely the tangential direction ($x$), axial direction ($y$), and radial direction ($z$).

In Figure 1, asymmetric collisions and the moments of aerodynamic excitation forces can excite the torsional vibration around $z$ axis. The first-order bending vibration of the blades along the $x$ direction and $y$ direction and the first-order torsion vibration of the blades around the $z$ axis will be coupled by the contact forces and aerodynamic excitation forces when adjacent shrouds are in contact. As the displacement along the $z$ direction is very small compared with that along the $x$ direction and $y$ direction, it is ignored in this paper.

In order to analyze the complex collision motion between adjacent shrouds, adjacent shrouded blades are simplified as two homogeneous mass blocks which are two parallelograms with the same size and shape, and the dynamic model of the shrouded blade system, which is a three-degree-of-freedom mass-spring model, is established based on [26] and shown in Figure 2.

In Figure 2, when the two mass blocks are in contact, the normal load $N$ caused by impact is simulated by a linear spring, while the friction force $f$ is modelled by the bilinear hysteresis model with normal load varying. It should be pointed out that $N$ and $f$ are forces exerted by the right mass on the left mass, and their action points which are determined by the kinematics analysis in the following section keep changing. The moments of $N$ and $f$ about points $o_1$ and $o_2$ can be calculated when the magnitude, direction, and acting point of the contact forces are determined. Meanwhile, the geometric dimensions of the mass blocks are identified in Figure 2: one of the vertex angles of the parallelogram shroud is $\alpha$; the lengths of the two line segments $ap$ and $ob$ are denoted as $l_{ap}$ and $l_{ob}$, $l_{ed}$ and $l_{cg}$ are denoted similarly.

The dynamic parameters of the system are denoted as follows. The equivalent masses of the two adjacent blades are denoted as $m_1$ and $m_2$, while the moments of inertia about the corresponding $z$ axis are denoted as $I_1$ and $I_2$. The first-order bending stiffnesses along the $x$ and $y$ directions are denoted as $k_{1x}$, $k_{1y}$ and $k_{2x}$, $k_{2y}$ respectively, and the corresponding linear damping coefficients are denoted as $c_{1x}$, $c_{1y}$ and $c_{2x}$, $c_{2y}$, respectively. The torsional stiffnesses about the corresponding $z$ axis are denoted as $k_{T1}$ and $k_{T2}$, and the corresponding torsional damping coefficients are denoted as $c_{T1}$ and $c_{T2}$. The centers of the two mass blocks are $o_1$ and $o_2$, $Q_1$ and $Q_2$ are principle vectors of aerodynamic excitation forces acting on the two mass blocks which pass through the center of the corresponding mass blocks, and the angles between $Q_1$, $Q_2$ and the $x$ axis are ${\beta }_1$ and ${\beta }_2$, respectively. $M_1$ and $M_2$ are the moments of the aerodynamic excitation forces about points $o_1$ and $o_2$, and the impact stiffness is denoted as $k$, while the initial gap between two adjacent shrouds is $\Delta$.

According to Figure 2 and the theory of vibration, the dynamic equation of the shrouded blade system can be written as

$$\left\{\begin{array}{l}{m}_1{\ddot{x}}_1+c_{1x}{\dot{x}}_1+k_{1x}{x}_1=Q_{1x}-N_x-f_x\\ m_1{\ddot{y}}_1+c_{1y}{\dot{y}}_1+k_{1y}{y}_1=Q_{1y}+N_y-f_y\\ m_2{\ddot{x}}_2+c_{2x}{\dot{x}}_2+k_{2x}{x}_2=Q_{2x}+N_x+f_x\\ m_2{\ddot{y}}_2+c_{2y}{\dot{y}}_2+k_{2y}{y}_2=Q_{2y}-N_y+f_y\\ I_1{\ddot{\theta }}_1+c_{T1}{\dot{\theta }}_1+k_{T1}{\theta }_1=M_1+M_{N1}+M_{f1}\\ I_2{\ddot{\theta }}_2+c_{T2}{\dot{\theta }}_2+k_{T2}{\theta }_2=M_2+M_{N2}+M_{f2}\end{array}\right.$$

(1)

In Equation (1), the projections of displacements of point $o_1$ along the $x$ and $y$ directions are $x_1$ and $y_1$, and $x_2$, $y_2$ are denoted similarly. The torsional angles around the corresponding $z$ axis of the two mass blocks are ${\theta }_1$ and ${\theta }_2$. The projections of $N$ and $f$ along the $x$ and $y$ axes are $N_x$ and $N_y$ and $f_x$ and $f_y$, respectively. The moments of $N$ and $f$ about point $o_1$ are $M_{N1}$ and $M_{f1}$, while the moments of the reaction of $N$ and $f$ about point $o_2$ are $M_{N2}$ and $M_{f2}$. The determination of $N_x$, $N_y$, $f_x$, $f_y$, $M_{N1}$, $M_{f1}$, $M_{N2}$, and $M_{f2}$ will be analyzed in Section 2.3.

2.2. The Dry Friction Force Model and Stick-Slip Transition Analysis

As shown in Figure 3, in this paper, the relative contact motion between two blades is point to surface contact. Normal collision motion leads to the variation of the normal load, while tangential contact movement leads to stick-slip transition, and the torsion motion makes the direction of the contact force vary. The bilinear hysteresis model including the varying normal load [5,8], which is suitable for simulating this kind of dry friction contact, is adopted. Note that the typical hysteretic behavior characterizing dry-friction phenomena can be simulated by adopting more accurate hysteretic models available in the literature [28,29].

In Figure 3, the friction force is modelled by a spring along the tangential direction with its stiffness $k_d$, and the spring has no initial length and can yield, while the normal load is modelled by a spring along the normal direction with its stiffness $k$. The dynamic friction coefficient is denoted as $\mu$. The stick-slip transition mechanism at the contact face is illustrated as follows. When the two mass blocks are in contact, point 1 represents the tip of one of the mass blocks, while point 2 is always attached to the other mass block in contact with point 1. Point s is the sliding contact point which can slide relative to point 2 with a limiting friction force $\mu N$. When contact happens, points 1, 2, and s coincide initially, contact surface is in a viscous state, and point s is attached to point 2. The displacement of point 1 relative to the contact surface along normal directions is denoted as $v_1$, while the displacements of points 1 and s relative to point 2 are denoted as $u_1$ and ${\mu }_s$, respectively. In contact motion, point s remains static with point 2 when the distance between points 1 and s is less than $\mu N/k_d$; otherwise, point s slides relative to point 2 and keeps the distance between point 1 and s as a constant $\mu N/k_d$. The friction force can be determined by Equation (2), whether it is viscous or sliding, and it is important to trace the position of point s.

$$f=k_d(u_1-u_s)$$

(2)

2.3. The Separation-Contact Motion Analysis and Determination of Contact Forces and Their Moments

When the aero engine is working, the relative motions of two adjacent shrouds are very complex, as the bending-torsion coupling vibration of the blades is considered. The bending-torsion coupling vibration of the blades can be simplified as the planar motion of two mass blocks, which consists of the translation with the corresponding center of the mass block and the rotation around the corresponding central $z$ axis. Displacements caused by torsional motion of the mass blocks are shown in Figure 4.

The torsional angle is assumed to be positive along the clockwise direction. As shown in Figure 4, the two mass blocks rotate at an angle of ${\theta }_1$ and ${\theta }_2$ separately around the z-axis, which passes through points $o_1$ and $o_2$, respectively, while points $a$, $b$, $c$, and $d$ move to a new position, $a^{\prime }$, $b^{\prime }$, $c^{\prime }$, and $d^{\prime }$. Thus, the displacement vectors of the four points are obtained. The projections of the torsional displacements of points $a$, $b$, $c$, and $d$ along the $x$ and $y$ axes can be described as Equations (3) and (4):

$$\left\{\begin{array}{c}{x}_a=l_{ap}\sin \frac{{\theta }_1}{2}\sin \frac{2{\phi }_1-{\theta }_1}{2}\\ y_a=-l_{ap}\sin \frac{{\theta }_1}{2}\cos \frac{2{\phi }_1-{\theta }_1}{2}\\ x_b=-l_{ob}\sin \frac{{\theta }_1}{2}\sin \frac{2{\phi }_2+{\theta }_1}{2}\\ y_b=-l_{ob}\sin \frac{{\theta }_1}{2}\cos \frac{2{\phi }_2+{\theta }_1}{2}\end{array}\right.$$

(3)

In Equation (3), $x_a$ and $x_b$ are the projections of the torsion displacements of points $a$ and $b$ along the $x$ direction, while $y_a$ and $y_b$ are the projections of points $a$ and $b$ along the $y$ direction.

$$\left\{\begin{array}{c}{x}_c=l_{cg}\sin \frac{{\theta }_2}{2}\sin \frac{2{\phi }_2+{\theta }_2}{2}\\ y_c=l_{cg}\sin \frac{{\theta }_2}{2}\cos \frac{2{\phi }_2+{\theta }_2}{2}\\ x_d=-l_{ed}\sin \frac{{\theta }_2}{2}\sin \frac{2{\phi }_1-{\theta }_2}{2}\\ y_d=l_{ed}\sin \frac{{\theta }_2}{2}\cos \frac{2{\phi }_1-{\theta }_2}{2}\end{array}\right.$$

(4)

In Equation (4), $x_c$ and $x_d$ are the projections of the torsion displacements of points $c$ and $d$ along the $x$ direction, while $y_c$ and $y_d$ are the projections of points $c$ and $d$ along the $y$ direction.

Based on the analysis above, the total displacement vectors of points $a$, $b$, $c$, and $d$, denoted as ${\overrightarrow{r}}_a$, ${\overrightarrow{r}}_b$, ${\overrightarrow{r}}_c$, and ${\overrightarrow{r}}_d$, are composed of the translational displacements and the torsional displacements. The computing formulas can be expressed as

$$\left\{\begin{array}{c}{r}_{ax}=x_1+x_a\\ {\mathrm{r}}_{ay}=y_1+y_a\\ {\mathrm{r}}_{bx}=x_1+x_b\\ {\mathrm{r}}_{by}=y_1+y_b\end{array}\right.$$

(5)

$$\left\{\begin{array}{c}{r}_{cx}=x_2+x_c\\ {\mathrm{r}}_{cy}=y_2+y_c\\ {\mathrm{r}}_{dx}=x_2+x_d\\ {\mathrm{r}}_{dy}=y_2+y_d\end{array}\right.$$

(6)

In Equations (5) and (6), $r_{ax}$, $r_{bx}$, $r_{cx}$, and $r_{dx}$ are the projections of the total displacements of points $a$, $b$, $c$, and $d$ along the $x$ direction, while $r_{ay}$, $r_{by}$, $r_{cy}$, and $r_{dy}$ are the projections of points $a$, $b$, $c$, and $d$ along the $y$ direction.

Kinematics analysis shows that there are four kinds of contact-separation transition boundary conditions under bending-torsion coupling vibration. In each case, the friction force $f$ is calculated with reference to Equation (2).

• The First Case: ${\theta }_1>{\theta }_2$, $r_{ay}-r_{cy}<0$

As shown in Figure 5, when point $a$ happens to be in contact with the right mass block, the difference between projections of displacements of points $a$ and $c$ along the direction perpendicular to $cd$ is equal to the projection of the initial gap $\Delta$ along this direction. The first contact-separation transition boundary condition can be described as Equation (7):

$$(r_{ax}-r_{cx})\sin (\alpha -{\theta }_2)+(r_{ay}-r_{cy})\cos (\alpha -{\theta }_2)=\Delta \sin (\alpha -{\theta }_2)$$

(7)

Obviously, there are two states of separation and contact.

• Separation: when $(r_{ax}-r_{cx}-\Delta )\sin (\alpha -{\theta }_2)+(r_{ay}-r_{cy})\cos (\alpha -{\theta }_2)<0$, the two masses are separated, and the contact forces and corresponding moments about points $o_1$ and $o_2$ are all equal to zero.
•  
• Contact: when $(r_{ax}-r_{cx}-\Delta )\sin (\alpha -{\theta }_2)+(r_{ay}-r_{cy})\cos (\alpha -{\theta }_2)\ge 0$, the two masses are in contact.

In the process of the contact motion under bending-torsion coupling vibration, the magnitude, direction, and acting point of contact forces change constantly, which can be determined by kinematics and dynamics analysis, and the corresponding calculation formulas are derived in Equations (8)–(10):

$$\left\{\begin{array}{l}{v}_1=(r_{ax}-r_{cx}-\Delta )\sin (\alpha -{\theta }_2)+(r_{ay}-r_{cy})\cos (\alpha -{\theta }_2)\\ N=k{v}_1\end{array}\right.$$

(8)

$$\left\{\begin{array}{c}{f}_x=f\sin (\alpha -{\theta }_2)\\ f_y=f\cos (\alpha -{\theta }_2)\\ N_x=N\sin (\alpha -{\theta }_2)\\ N_y=N\cos (\alpha -{\theta }_2)\end{array}\right.$$

(9)

$$\left\{\begin{array}{c}{M}_{N1}=-N\times \frac{l_{ap}}{2}\cos (\alpha -{\phi }_1+{\theta }_1-{\theta }_2)\\ M_{f1}=f\times \frac{l_{ap}}{2}\sin (\alpha -{\phi }_1+{\theta }_1-{\theta }_2)\\ M_{N2}=N\times (\frac{l_{cd}}{2}\cos (\alpha +{\phi }_2)+\frac{r_{cy}-r_{ay}}{\sin (\alpha -{\theta }_2)})\\ M_{f2}=f\times \frac{l_{ed}\sin (\alpha -\phi 1)}{2}\end{array}\right.$$

(10)

• The Second Case: ${\theta }_1>{\theta }_2$, $r_{ay}-r_{cy}\ge 0$

As shown in Figure 6, when point $c$ happens to be in contact with the left mass block, the difference between projections of displacements of points $a$ and $c$ along the direction perpendicular to $ab$ is equal to the projection of the initial gap $\Delta$ along this direction. The second contact-separation transition boundary condition is described as Equation (11):

$$(r_{ax}-r_{cx})\sin (\alpha -{\theta }_1)+(r_{ay}-r_{cy})\cos (\alpha -{\theta }_1)=\Delta \sin (\alpha -{\theta }_1)$$

(11)

Obviously, there are two states of separation and contact.

• Separation: when $(r_{ax}-r_{cx}-\Delta )\sin (\alpha -{\theta }_1)+(r_{ay}-r_{cy})\cos (\alpha -{\theta }_1)<0$, the two masses are separated, and the contact forces and corresponding moments about points $o_1$ and $o_2$ are all equal to zero.
•  
• Contact: when $(r_{ax}-r_{cx}-\Delta )\sin (\alpha -{\theta }_1)+(r_{ay}-r_{cy})\cos (\alpha -{\theta }_1)\ge 0$, the two masses are in contact.

The magnitude, direction, and acting point of contact forces vary during the contact motion and can be determined by kinematics and dynamics analysis; the corresponding calculation formulas are derived in Equations (12)–(14):

$$\left\{\begin{array}{l}{v}_1=(r_{ax}-r_{cx}-\Delta )\sin (\alpha -{\theta }_1)+(r_{ay}-r_{cy})\cos (\alpha -{\theta }_1)\\ N=k{v}_1\end{array}\right.$$

(12)

$$\left\{\begin{array}{c}{f}_x=-f\sin (\alpha -{\theta }_1)\\ f_y=-f\cos (\alpha -{\theta }_1)\\ N_x=N\sin (\alpha -{\theta }_1)\\ N_y=N\cos (\alpha -{\theta }_1)\end{array}\right.$$

(13)

$$\left\{\begin{array}{l}{M}_{N1}=-N\times (\frac{l_{ap}}{2}\cos (\alpha -{\phi }_1)-\frac{r_{ay}-r_{cy}}{\sin (\alpha -{\phi }_1)})\\ M_{f1}=-f\times \frac{l_{ap}\sin (\alpha -{\phi }_1)}{2}\\ M_{N2}=-N\times \frac{l_{cg}}{2}\cos (\alpha +{\phi }_2-{\theta }_1+{\theta }_2)\\ M_{f2}=-f\times \frac{l_{cg}}{2}\sin (\alpha +{\phi }_2-{\theta }_1+{\theta }_2)\end{array}\right.$$

(14)

• The Third Case: ${\theta }_1\le {\theta }_2$, $r_{by}-r_{dy}>0$

As shown in Figure 7, when point $b$ happens to be in contact with the right mass block, the difference between projections of displacements of points $b$ and $d$ along the direction perpendicular to $cd$ is equal to the projection of the initial gap $\Delta$ along this direction. The third contact-separation transition boundary condition is described as Equation (15):

$$(r_{bx}-r_{dx})\sin (\alpha -{\theta }_2)+(r_{by}-r_{dy})\cos (\alpha -{\theta }_2)=\Delta \sin (\alpha -{\theta }_2)$$

(15)

Obviously, there are two states of separation and contact.

• Separation: when $(r_{bx}-r_{dx}-\Delta )\sin (\alpha -{\theta }_2)+(r_{by}-r_{dy})\cos (\alpha -{\theta }_2)<0$, the two masses are separated, and the contact forces and corresponding moments about points $o_1$ and $o_2$ are all equal to zero.
•  
• Contact: when $(r_{bx}-r_{dx}-\Delta )\sin (\alpha -{\theta }_2)+(r_{by}-r_{dy})\cos (\alpha -{\theta }_2)\ge 0$, the two masses are in contact.

During the contact motion, the magnitude, direction, and acting point of contact forces change and can be determined by kinematics and dynamics analysis; the corresponding calculation formulas are derived in Equations (16) and (17) ($N_x$, $N_y$ and $f_x$, $f_y$ are still obtained by Equation (9)):

$$\left\{\begin{array}{l}{v}_1=(r_{bx}-r_{dx}-\Delta )\sin (\alpha -{\theta }_2)+(r_{by}-r_{dy})\cos (\alpha -{\theta }_2)\\ N=k{v}_1\end{array}\right.$$

(16)

$$\left\{\begin{array}{l}{M}_{N1}=-N\times \frac{l_{ob}}{2}\cos (\alpha +{\phi }_2-{\theta }_2+{\theta }_1)\\ M_{f1}=f\times \frac{l_{ob}}{2}\sin (\alpha +{\phi }_2-{\theta }_2+{\theta }_1)\\ M_{N2}=-N\times (\frac{l_{ap}}{2}\cos (\alpha -{\phi }_1)-\frac{r_{by}-r_{dy}}{\sin (\alpha -{\theta }_2)})\\ M_{f2}=f\times \frac{l_{ed}\sin (\alpha -{\phi }_1)}{2}\end{array}\right.$$

(17)

• The Fourth Case: ${\theta }_1\le {\theta }_2$, $r_{by}-r_{dy}\le 0$

As shown in Figure 8, when point $d$ happens to be in contact with the left mass block, the difference between projections of displacements of points $b$ and $d$ along the direction perpendicular to $ab$ is equal to the projection of the initial gap $\Delta$ along this direction. The fourth contact-separation transition boundary condition can be described as Equation (18):

$$(r_{bx}-r_{dx})\sin (\alpha -{\theta }_1)+(r_{by}-r_{dy})\cos (\alpha -{\theta }_1)=\Delta \sin (\alpha -{\theta }_1)$$

(18)

Obviously, there are two states of separation and contact.

• Separation: when $(r_{bx}-r_{dx}-\Delta )\sin (\alpha -{\theta }_1)+(r_{by}-r_{dy})\cos (\alpha -{\theta }_1)<0$, the two masses are separated, and the contact forces and corresponding moments about points $o_1$ and $o_2$ are all equal to zero.
•  
• Contact: when $(r_{bx}-r_{dx}-\Delta )\sin (\alpha -{\theta }_1)+(r_{by}-r_{dy})\cos (\alpha -{\theta }_1)\ge 0$, the two masses are in contact.

Contact forces change during the contact motion, and their magnitude, direction, and acting point can be determined by kinematics and dynamics analysis; the corresponding calculation formulas are derived in Equations (19) and (20) ($N_x$, ${\phi }_1$ and $f_x$, $f_y$ are still obtained by Equation (13)):

$$\left\{\begin{array}{l}{v}_1=(r_{bx}-r_{dx}-\Delta )\sin (\alpha -{\theta }_1)+(r_{by}-r_{dy})\cos (\alpha -{\theta }_1)\\ N={kv}_1\end{array}\right.$$

(19)

$$\left\{\begin{array}{l}{M}_{N1}=N\times (\frac{l_{cg}}{2}\cos (\alpha +{\phi }_2)+\frac{r_{dy}-r_{by}}{\sin (\alpha -{\theta }_2)})\\ M_{f1}=-f\times \frac{l_{ap}\sin (\alpha -{\phi }_1)}{2}\\ M_{N2}=-N\times \frac{l_{ed}}{2}\cos (\alpha -{\phi }_1+{\theta }_2-{\theta }_1)\\ M_{f2}=-f\times \frac{l_{ed}}{2}\sin (\alpha -{\phi }_1+{\theta }_2-{\theta }_1)\end{array}\right.$$

(20)

3. Numerical Simulation of the Dynamics of the System

In this section, the rub-impact dynamics of shrouded blades is studied in detail. Firstly, the characteristics of the forced response of the system is analyzed, which is the base of the vibration reduction characteristics of the shrouded blade system. Secondly, the method to evaluate the vibration reduction characteristics of the shrouded blade system when the motion of the blade is chaotic or has fractional harmonics is studied. Finally, the vibration reduction characteristics of the shrouded blade system is discussed, which is important for the design of the dry friction damper.

3.1. The Characteristics of the Forced Response

Select $m_1=m_2=m$, $I_1=I_2=I$, $c_{1x}=c_{2x}=c_x$, $c_{1y}=c_{2y}=c_y$, $k_{1x}=k_{2x}=k_x$, $k_{1y}=k_{2y}=k_y$, $c_{T1}=c_{T2}=c_T$, $k_{T1}=k_{T2}=k_T$, $Q_i=F_0\sin (\omega t+{\varphi }_i),i=1,2$, $M_i=M_0\sin (\omega t+{\varphi }_i),i=1,2$. Denote ${\omega }_x^2=\frac{k_x}{m}$, ${\omega }_y^2=\frac{k_y}{m}$, ${\epsilon }_1=\frac{c_x}{2m{\omega }_x}$, ${\epsilon }_2=\frac{c_y}{2m{\omega }_y}$, ${\omega }_T^2=\frac{k_{\mathrm{T}}}{I}$, ${\epsilon }_T=\frac{c_T}{2I{\omega }_T}$, $X_i=\frac{x_i}{\Delta },i=1,2$, $Y_i=\frac{y_i}{\Delta },i=1,2$, ${\mathrm{\Theta }}_i=\frac{{\theta }_i}{\alpha },i=1,2$, $\tau =\omega t$, ${(•)}^{\prime }=\frac{d•}{d\tau }$, ${(•)}^{″}=\frac{d^2•}{d{\tau }^2}$, $\overline{f_x}=\frac{f_x}{m\Delta {\omega }_x^2}$, $\overline{f_y}=\frac{f_y}{m\Delta {\omega }_y^2}$, $\overline{Q_{ix}}=\frac{Q_i\cos \beta }{m\Delta {\omega }_x^2},i=1,2$, $\overline{Q_{iy}}=\frac{Q_i\sin \beta }{m\Delta {\omega }_y^2},i=1,2$, $\overline{M_{iT}}=\frac{M_i}{I\alpha {\omega }_T^2},i=1,2$, $\overline{N_x}=\frac{N_x}{m\Delta {\omega }_x^2}$, $\overline{N_y}=\frac{N_y}{m\Delta {\omega }_y^2}$, $\overline{M_{Ni}}=\frac{M_{Ni}}{I\alpha {\omega }_T^2},i=1,2$, $\overline{M_{fi}}=\frac{M_{fi}}{I\alpha {\omega }_T^2},i=1,2$. $F_0$ is the amplitude of $Q_1$ or $Q_2$, $M_0$ is the amplitude of $M_1$ or $M_2$, while ${\varphi }_1$ and ${\varphi }_2$ are the corresponding phase angles. Defining these will simplify Equation (1) as Equation (21):

$$\left\{\begin{array}{c}{X}_1^{″}+2{\epsilon }_1{X}_1^{\prime }+X_1=\overline{Q_{1x}}-\overline{N_x}-\overline{f_x}\hfill \\ Y_1^{″}+2{\epsilon }_2{Y}_1^{\prime }+Y_1=\overline{Q_{1y}}+\overline{N_y}-\overline{f_y}\hfill \\ X_2^{″}+2{\epsilon }_1{X}_2^{\prime }+X_2=\overline{Q_{2x}}+\overline{N_x}+\overline{f_x}\hfill \\ Y_2^{″}+2{\epsilon }_2{Y}_2^{\prime }+Y_2=\overline{Q_{2y}}-\overline{N_y}+\overline{f_y}\hfill \\ {\mathrm{\Theta }}_1^{″}+2{\epsilon }_T{\mathrm{\Theta }}_1^{\prime }+{\mathrm{\Theta }}_1=\overline{M_{1T}}-\overline{M_{N1}}+\overline{M_{f1}}\hfill \\ {\mathrm{\Theta }}_2^{″}+2{\epsilon }_T{\mathrm{\Theta }}_2^{\prime }+{\mathrm{\Theta }}_2=\overline{M_{2T}}+\overline{M_{N2}}+\overline{M_{f2}}\hfill \end{array}\right.$$

(21)

Based on the dynamic equation of the system and the transition analysis of the separation-contact motion and stick-slip motion, a numerical integration scheme implementing a fourth-order Runge–Kutta algorithm [30] is set up to analyze the characteristics of the forced response. In order to improve the simulation accuracy, the bisection method is used to capture the transition points of stick-slip motion and contact-separation motion. The influences of stiffness ratio, initial gap, and amplitude of the external excitation force on the response characteristics of the shrouded blade system are discussed. Denote $f_e=\omega /2\pi$.

As the characteristics of response of the blocks along all directions are identical, only the characteristics of $X_1$ are shown. Referring to [26], the common parameters are shown in

Table 1. Common parameters.

Parameters	Values	Parameters	Values
$m$	0.08 kg	$c_x$	3.57 N·s/m
$k_x$	1 × 105 N/m	$c_y$	19.6 N·s/m
$k_y$	3 × 106 N/m	$\beta$	$\pi /6$
$k_d$	1 × 106 N/m	$\mu$	0.5
${\varphi }_1$	0	${\varphi }_2$	$\pi$
$I$	2.46 × 10−3 kg·m2	$c_T$	0.506 N·m·s/rad
$k_T$	2 × 103 N·m/rad	$M_0$	100 N·m
$\alpha$	$\pi /3$	$l_{oa}$	0.05 m
$l_{ab}$	0.035 m	$\omega$	536 rad/s

, and specific parameters are given in the corresponding section.

3.1.1. Influence of Stiffness Ratio on Characteristics of the Response of the Shrouded Blade

The stiffness ratio is denoted as $\gamma =k/k_x$. The specific parameters are shown in

Table 2. The specific parameters.

Parameters	Values
$\Delta$	0.022 mm
$F_0$	60 N

and the simulation results are shown in Figure 9 and Figure 10.

Figure 9 shows the bifurcation diagram of $X_1$ vs. $\gamma$, which is illustrated by Figure 10. Analysis of the Poincare map, the spectrum, and the phase diagram shows the following conclusions. In Figure 10a, where $\gamma =0.3$, only higher harmonics exist, and the response is period-1 motion. In Figure 10b, where $\gamma =2.5$, there are fractional harmonic frequencies ($f_e/2,3f_e/2\dots$), and period-2 motion occurs. In Figure 10c, where $\gamma =4.2$, there are fractional harmonic frequencies ($f_e/3,4f_e/3\dots$), and period-3 motion occurs. In Figure 10d, where $\gamma =7.5$, chaos occurs. In summary, with the increase of $\gamma$, the responses follow a sequence of period-1, period-2, period-3, and chaotic motion.

3.1.2. Influence of the Initial Gap on Characteristics of the Response of the Shrouded Blade

The specific parameters and simulation results are shown in

Table 3. The specific parameters.

Parameters	Values
$\gamma$	2.6
$F_0$	60 N

and Figure 11 and Figure 12, respectively.

Figure 11 shows the bifurcation diagram of $X_1$ vs. $\Delta$, which is illustrated by Figure 12. Analysis of the Poincare map, the spectrum, and the phase diagram show the following conclusions. In Figure 12a, where $\Delta =0.007\mathrm{mm}$, the motion of the blade is chaotic. In Figure 12b, where $\Delta =0.02\mathrm{mm}$, there are fractional harmonic frequencies ($f_e/4,3f_e/2\dots$), and the response is period-4 motion. In Figure 12c, where $\Delta =0.06\mathrm{mm}$, there are fractional harmonic frequencies ($f_e/2,3f_e/2\dots$), and period-2 motion occurs. In Figure 12d, where $\Delta =0.15\mathrm{mm}$, the impact and dry friction forces become very small, and period-1 motion occurs. In summary, for the factor of $\Delta$, with the increase of its value, the responses follow a sequence of chaotic, period-4, period-2, and period-1 motion.

3.1.3. Influence of the Amplitude of Aerodynamic Excitation Forces on Characteristics of the Response of the Shrouded Blade

The specific parameters are shown in

Table 4. The specific parameters.

Parameters	Values
$\Delta$	0.022 mm
$\gamma$	2.6

, while the simulation results are shown in Figure 13 and Figure 14.

Figure 13 shows the bifurcation diagram of $X_1$ vs. $F_0$, which is illustrated by Figure 14. Analysis of the Poincare map, the spectrum, and the phase diagram show the following conclusions. In Figure 14a, where $F_0=10 \mathrm{N}$, only higher harmonics exist, and period-1 motion occurs. In Figure 14b, where $F_0=37 \mathrm{N}$, fractional harmonic frequencies ($f_e/4,3f_e/2\dots$) exist, and period-4 motion occurs. In Figure 14c, where $F_0=50 \mathrm{N}$, the response is chaotic. In Figure 14d, where $F_0=67 \mathrm{N}$, fractional harmonics frequencies ($f_e/2,3f_e/2\dots$) exist, and period-2 motion occurs. In summary, with the increase of the value of $F_0$, the responses follow a sequence of period-1, period-4, chaotic, and period-2 motion.

3.2. Evaluation of the Vibration Reduction Characteristics of the System

In Section 3.1, the response characteristics of the shrouded blade system were analyzed. It is clear that fractional harmonics or chaos may appear in the forced response of the blade under some parameters. In this section, how to evaluate the vibration reduction characteristics of the chaotic response is proposed.

The average power of the response is introduced to describe the vibration reduction characteristics of the system. In order to denote the average power of the response of a bending-torsion coupling system, a fixed orthogonal coordinate system $o_{1}XY$, which is parallel to the coordinate system $o_{1}xy$, is established in the centroid of the left mass block, and point $L$ is an arbitrary point on the mass block.

The total displacement of point $L$ is denoted as ${\overrightarrow{u}}_L\left(X,Y,t\right)$ whose projections along the $X$ direction and $Y$ direction are $r_{LX}$ and $r_{LY}$. According to the response of the bending-torsion coupling system, $r_{LX}$ and $r_{LY}$ can be described as Equation (22). Denote $T=2\pi /\omega$. The average power of the steady-state response of the bending-torsion coupling system is described as Equation (23), in which $P$ and $P_0$ correspond to the average power of impact occurrence state and no impact occurrence state during the whole movement of the shrouded blade system, respectively. In Equation (23), the area covered by the left mass block is denoted as $A$, while $n$ (a positive integer) is the number of the time periods of integration which should be selected when the response of the system is steady. $P_d$ is the percentage of the average steady-state power reduction of the system, which is defined as Equation (24). When no impact occurs between the two mass blocks during the whole motion, the total displacement of the arbitrary point on the mass block is ${\overrightarrow{u}}_{0L}\left(X,Y,t\right)$.

$$\left\{\begin{array}{l}{r}_{LX}(X,Y,t)=x_1(t)+2\sqrt{X^2+Y^2}\sin (\frac{{\theta }_1(t)}{2})\sin (\phi -\frac{{\theta }_1(t)}{2})\\ r_{LY}(X,Y,t)=y_1(t)+2\sqrt{X^2+Y^2}\sin (\frac{{\theta }_1(t)}{2})\cos (\phi -\frac{{\theta }_1(t)}{2})\end{array}\right.$$

(22)

$$\left\{\begin{array}{c}P=\frac{{\displaystyle {\mathsf{\int }}_{nT}{\displaystyle \underset{A}{\mathsf{\iint }}{\left|{\overrightarrow{u}}_L(X,Y,t)\right|}^2}dXdYdt}}{nT}\\ P_0=\frac{{\displaystyle {\mathsf{\int }}_{nT}{\displaystyle \underset{A}{\mathsf{\iint }}{\left|{\overrightarrow{u}}_{0L}(X,Y,t)\right|}^2}dXdYdt}}{nT}\end{array}\right.$$

(23)

$$P_d=\frac{P_0-P}{P_0}\times 100\%$$

(24)

When no impact occurs during the whole motion, according to the linear vibration theory, the steady-state response of the system must be periodic motion without fractional harmonics. Based on analysis of the system response characteristics in Section 3.1, when impact occurs, the steady state responses of the bending-torsion coupling system can be divided into three types, periodic motion with higher harmonics, periodic motion with fractional harmonics, and chaotic motion. According to the three different motions, how to determine the value of $n$, which is very important for the evaluation of the vibration reduction characteristics of the system, is discussed:

• When impact occurs and the steady state response of the blade is periodic motion with higher harmonics, $n$ can be taken as any integer greater than or equal to 1;
• When impact occurs and the steady-state response is periodic motion and has $f_e/i$ ($i$ is a positive integer) harmonics in the frequency spectrum, $n$ can be taken as any positive integral multiple of $i$;
• When impact occurs and the steady-state response of the blade is chaotic motion, we discuss whether or not the value of $P_d$ approach a stable value with the increase of the value of $n$.

As the simulation results under selected parameters are consistent, only the simulation results under two parameters are shown here. The two sets of parameters are taken from

Table 5. Two parameters of the system.

Parameters	Values	Values
$\Delta$	0.022 mm	0.022 mm
$\gamma$	7	2.6
$F_0$	60 N	50 N

, and the simulation results are shown in Figure 15.

In Figure 15, it can be seen that the value of $P_d$ eventually stabilizes with the increase of the value of $n$. Therefore, when the chaos occurs, we can select the stable value of $P_d$ for the vibration reduction analysis of the system.

3.3. The Vibration Reduction Characteristics of the Bending-Torsion Coupling System

In this section, the influences of some key parameters are discussed, including stiffness ratio, the initial gap, and external excitation amplitude on the vibration reduction characteristics of the bending-torsion coupling system. Some main parameters have been shown in Table 1, while other parameters are given in the following simulation.

3.3.1. Influence of the Stiffness Ratio on Vibration Reduction Characteristics of the System

The parameters needed are taken from

Table 6. Parameters needed.

Parameters	Values
$\Delta$	0.022 mm
$F_0$	60 N

. The simulation results of the influence of the stiffness ratio $\gamma$ on the vibration reduction characteristics of the system are shown in Figure 16.

In Figure 16, when $\gamma =0$ and $k=0$, the normal load and the friction force are equal to zero, therefore, $P_d$ is equal to zero. With the increase of the value of $\gamma$, the responses follow a sequence of period motion with higher harmonics, period motion with fractional harmonics, and chaotic motion, while the value of $P_d$ first increases and then decreases to an approximate stable value as the normal load and friction force increase. There exists a value of $\gamma$ with which the average power reduction of the system is the best and the system is in stick-slip motion. With the larger value of $\gamma$, $P_d$ stabilizes and chaotic motion occurs, the reason for which being that viscous motion plays a dominant role in stick-slip motion. The law of variation of $P_d$ with $\gamma$ is mainly related to the complex separation-stick-slip motion.

3.3.2. Influence of the Initial Gap on the Vibration Reduction Characteristics of the System

The parameters needed are taken from

Table 7. Parameters needed.

Parameters	Values
$\gamma$	4
$F_0$	60 N

. The simulation results of the influence of the initial gap on the vibration reduction characteristics of the system is shown in Figure 17.

In Figure 17, with the increase of the value of $\Delta$ from zero, the responses follow a sequence of chaotic motion, period motion with fractional harmonics, and period motion with higher harmonics, while the value of $P_d$ keeps decreasing as the amplitude of normal load and friction force decreases. When the value of $\Delta$ is larger and larger, the two masses of the system are separated during the whole motion, and the value of $P_d$ is equal to zero. In the range of the selected values of $\Delta$, the steady system has no pure viscous motion.

3.3.3. The Influence of the Amplitude of Aerodynamic Excitation Forces on the Vibration Reduction Characteristics of the System

The parameters needed are taken from

Table 8. Parameters needed.

Parameters	Values
$\gamma$	4
$\Delta$	0.022 mm

. The simulation results of the influence of the amplitude of excitation force on the vibration reduction characteristics of the system is shown in Figure 18.

In Figure 18, with the increase of the value of $F_0$, the responses follow a sequence of period motion with higher harmonics, period motion with fractional harmonics, and chaotic motion, while the value of $P_d$ approximately increases to an approximate stable value. The law of variation of $P_d$ with $F_0$ is mainly due to the increase of the amplitude of normal load and friction force as $F_0$ increases. The reason for the stabilization of $P_d$ is that viscous motion plays a dominant role when the system is steady, as the normal load and friction force increase to a certain extent.

4. Conclusions

In this paper, the bending-torsion coupling vibration of shrouded blades is considered, and the complex contact-separation transition boundary conditions are set up. On the basis of the analysis above, the response characteristics, how to evaluate the vibration reduction characteristics of the bending-torsion coupling system, and the vibration reduction characteristics of the system are discussed. Major conclusions can be drawn as follows:

• When impact occurs in the motion, with parameters varying, the steady-state response of the bending-torsion coupling system displays three types of motions which are periodic motion with higher harmonics, periodic motion with fractional harmonics, and chaotic motion;
• When chaotic motion occurs, the average power of the bending-torsion coupling system eventually stabilizes with the increase of the integral period number;
• When the normal load and friction force increase to a certain extent, viscous motion will play a dominant role in the steady system, and the effect of the bending-torsion vibration reduction is obvious and tends to be approximately stable. The best steady-state vibration reduction effect of the bending-torsion coupling system relates to complex separation-stick-slip motion.

The dynamic model in this paper is a spring-mass model. In future research, the complex contact-separation transition boundary conditions under bending-torsion vibration can be developed and introduced into the finite element model to obtain more meaningful conclusions for engineering design.