Nonlinear Dynamics and Vibration Localization of Shrouded Blisk with Contact and Friction Effects

Abstract

Shrouds have been widely used to reduce the level of excessive vibration of blisks. However, complicated nonlinear motion can be induced by the contact and friction between shrouds. Even worse, harmful localization of vibration can be encountered due to nonsmooth behaviour at the shroud contact interfaces in the tuned disk. In this paper, the nonlinear dynamics and localization of vibration of a shrouded blisk are studied considering the spin softening. The continuous parametric model of a shrouded blisk is established, and the transition boundaries between different status (i.e., stick, slip and separation) of shrouds are determined based on the Coulomb friction model. The steady-state responses of the blisk are analysed using variable rotation speed, and the primary resonance, beat and quasi-periodic vibration are presented in connection with the non-smooth behaviour of contact. A particular type of vibration localization initiated by the contact and friction effect of the shrouds in the form of asymmetric vibration in the tuned blisk is discovered. It is found that the blades-disk coupling has strong influences on the level of the localization. The effects of the contact stiffness and localization on the blade vibration are demonstrated through the change in the powers of the blade motions.

1. Introduction

High-speed rotating disks are commonly used in gas turbine and steam turbines that operate in aeroengines and power generators. Under various circumstances, excessive vibration takes place due to the action of large aerodynamic and mechanical forces, which may lead to the eventual fatigue of the blades [1]. Shrouds have been widely used by steam turbines developers to dissipate energy and reduce stress of the blades, so as to improve the life of the blade [2]. However, additional effect of structural coupling can be introduced to adjacent shrouds through contact and friction on the mating interfaces, which leads to nonlinear characteristics such as non-smoothness and hence makes accurate prediction of vibration response of the disk difficult to achieve.

Extensive research has been devoted to the modeling of shrouded blisks including structurally modeling of blisks and the corresponding contacts of adjacent shrouds. A tuned, cyclically symmetrical blisk can be simplified as a single bladed sector [3], integrally shrouded group blade [4,5] or blade packet [6]. To admit the dynamical characteristics of the disk, Petrov [7,8] used the surface elements to simulate the contact between adjacent shrouds of the full structure. Beam-like blades and a rigid disk were coupled to a continuous parameter model of a blisk by Shadmani [9]. As for the contact model, various methods have been proposed to characterize the impact and friction of shrouds. In early publications, Sgn-contact, macro-slip and micro-slip models were proposed to simulate the lumped mass model of shrouds by Iwan [10], Griffin [11] and Menq [12], respectively. In subsequent studies, Cigeroglu et al. [13] put forward a two-dimensional distributed parameter microslip friction model considering the change in normal pressure caused by motion in normal direction. Yang and Menq [14] proposed a three-dimensional contact model to predicted the resonant response of structures having 3D frictional constraint.

It was known previously that collision and friction are in their nature nonlinear or, at least, piecewise linear, due to the change in relative velocity of vibrating bodies on their mating interfaces [15,16], and are attributable to complicated dynamic behaviors such as bifurcation, chaos and switch of stability [17]. For a slender structure undergoing impact, Liu et al. [18] proposed a model of fractal geometry and analyzed the effect of the surface roughness and normal loads on dynamic responses of blades, Li et al. [19] established a macro-slip friction model of the contact interfaces at the root of the beam with a dovetail tenon to characterize the friction on the beam. As for the friction between shrouds, He et al. [20] established an integrally shrouded blades considering impact and friction between adjacent shrouds, and reported multi-periodic, quasi-periodic and chaotic vibration of blade due to nonsmooth behaviour at the shroud contact surfaces. As far as the nonlinear vibration of a full blisk is concerned, progresses have been achieved in veering and merging analysis [21,22], localization of vibration [23,24], robust analysis [25] and parametric sensitivity [26]. In addition, some scholars study the contact and friction of blades based on the micro-slip contact model, Pesaresi et al. [27] presented a new modelling approach for underplatform dampers and evaluated against the experimental data of a recently-developed test rig, Chen et al. [28] proposed a new contact slip modeling method which can preserve the pressure distribution of the joint’s contact surface to capture the micro-slip phenomenon.

As for the localization of vibration, most of the existing researches were carried out under the assumption that the vibration localization is caused by the mistuning of blisks [29,30,31,32]. A few publications were dedicated to the vibration localization of cyclically symmetric structures with nonlinearities. Vakakis et al. [33,34] found that nonlinear mode localization can occur in perfectly symmetric nonlinear periodic systems, and the only prerequisite for its existence is weak coupling between subsystems. Based on the researches of Vakakis, Grolet and Thouverez [35] computed the free and forced response of a system with cyclic symmetry under geometric nonlinearity, and found that some of the bifurcated solutions correspond to localized nonlinear modes in the branching point bifurcation of the system. In addition, Fontanela et al. [36] developed a fully numerical approach to compute quasi-periodic vibrations bifurcating and localized oscillations from nonlinear periodic states in cyclic and symmetric structures.

Despite many existent publications regarding vibration of rotating shrouded blisks, there still lacks research that focused on interactions of adjacent shrouds through contact and friction, and their possible contribution to the localization of vibration of the disk. Such localization is usually accompanied with harmful concentration of vibration energy in one or a few blades, and is considered one of the major inducements to the onset of fatigue failure.

In this paper, the nonlinear dynamics and the localization of vibration of a tuned shrouded blisk that has short blades of small installation angles are investigated considering the effect of contact and friction between adjacent shrouds. The paper is organized as follows: Firstly, the continuous parametric model of the whole shrouded blisk is presented where the interaction between shrouds is modeled using linear spring and velocity-dependent force. The transition of contact states (i.e., stick, slip and separation) of adjacent shrouds are determined. Secondly, the vibration characteristic of the blades such as primary resonance, beat phenomenon and quasi-periodic vibration are analysed with the increase rotation speed. Finally, a particular type of vibration localization due to nonsmooth behaviour at the shroud contact interfaces is discovered in the tuned blisk, the effects of the blades-disk coupling on the level of localization are illustrated, and the effects of the normal contact stiffness and the localization of vibration on the blade vibration are demonstrated.

2. Dynamics Modelling of Shrouded Blisk

As shown in Figure 1, a rotating shrouded blisk with short blades of small installation angles is concerned in this study, where $\mathsf{\Omega }$ represents rotation speed. Structurally, the blisk is comprised of a disk and multiple identical blades that crowned with shrouds.

One may notice that for such a blisk, the dynamical modeling can be simplified as an in-plane dynamic model for the following reasons:

• The torsion of the blade body in the spanwise direction is generally small in view of turbine design for short blades. When short blades are concerned, the change in the circumferential velocity in the blade’s spanwise direction can be assumed small. Consequently, the directions of the relative velocity of the fluid to the blade vary limitedly in the same direction, and thus the torsion angle is idealized to be small for short blades.
• The bending rigidity of the blade is generally large compared with a long blade, and the modal couplings between in-plane and out-of-plane motions due to geometrically nonlinearity is weak enough to be ignored.
• When the installation angle is small, the aerodynamic force normal to the blade surface stands mainly in the plane of blisk rotation [20], meaning the in-plane component of the force is much larger than its out-of-plane counterpart. Under the same assumption of small installation angle, the principal axes of the blade cross-section are arranged in such a way that the bending rigidity in the out-of-plane direction overwhelms that in the in-plane directions. Hence, it is reasonable not to consider the out-of-plane motions in the present study owing to the small load and large bending rigidity in the out-of-plane.

2.1. Dynamic Modeling

Figure 2 depicts the model of the shrouded blisk, where two cartesian frames are established to express the blade motions. Frame XYZ is inertial cartesian coordinate, where X-axis is parallel to the ground and Z-axis is parallel to the normal of the disk. The other frame, $x_i{y}_i{z}_i$, is a body coordinate fixed at the barycentre of the i-th blade root, and rotates synchronously with the blade. Let (I, J, K) and ($i_i,j_i,k_i$) be the unit base vectors of the two frames, respectively bounded by the following transformation

$$\left[\begin{array}{l}{i}_i\hfill \\ j_i\hfill \\ k_i\hfill \end{array}\right]=\left[\begin{array}{ccc}\cos {\theta }_i& \sin {\theta }_i& 0\\ -\sin {\theta }_i& \cos {\theta }_i& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{l}I\hfill \\ J\hfill \\ K\hfill \end{array}\right]$$

(1)

Two degrees of freedom, $x_{\mathrm{d}}$ and $y_{\mathrm{d}}$ are adopted to locate the instantaneous center of disk in the X- and Y- directions, respectively. The contribution of the shaft to the disk is described by stiffness coefficients $k_x$ and $k_y$, and viscous damping coefficients $c_x$ and $c_y$, respectively. $d_{\mathrm{r}}$ is the radius of the disk. The angular orientation of the i-th blade is

$${\theta }_i=\mathsf{\Omega }t+\frac{2\pi }{n}\left(i-1\right)$$

(2)

where n is the total number of blades. In Figure 2b, $\alpha$ is the tilt angle of shrouds and L is the length of the blade. ${\delta }_{\mathrm{s}}$ and $k_{\mathrm{s}}$ are the gap and the normal contact stiffness between adjacent shrouds. $N_{\mathrm{s},i}$, $F_{\mathrm{s},i}$ and $Q_i$ are forces of normal contact, friction and aerodynamic [20]. $N_{\mathrm{s},i}$ and $F_{\mathrm{s},i}$ are perpendicular and parallel to the contact interfaces, respectively. $Q_i$ is parallel to the negative direction of the $y_i$-axis of the i-th blade.

$$Q_i=Q_0\sin \left(\lambda \mathsf{\Omega }t+{\phi }_i\right)$$

(3)

where $Q_0$ is the forcing amplitude, $\lambda$ the number of inlet vanes in the upstream of the blades [20] and ${\phi }_i=i\pi$ is the phase of $Q_i$.

The position of the i-th blade $R_{\mathrm{b},i}$, shroud $R_{\mathrm{s},i}$ and barycenter of the disk $R_{\mathrm{d}}$ at the current moment are

$$\begin{array}{l}{R}_{\mathrm{b},i}=x_{\mathrm{d}}I+y_{\mathrm{d}}J+d_{\mathrm{r}}{i}_i+\left(x_i+u_i\right)i_i+v_i{j}_i,\\ R_{\mathrm{s},i}=R_{\mathrm{b},i}\delta \left(x_i-L\right),\\ R_{\mathrm{d}}=x_{\mathrm{d}}I+y_{\mathrm{d}}J\end{array}$$

(4)

where $u_i$ and $v_i$ are spanwise and flapwise displacements of the blade. $\delta \left(x_i-L\right)$ is the Dirac-delta function.

As for the contact and friction forces between adjacent shrouds, it is noticed that the actual friction between the adjacent shrouds lie in a 2D surface defined in the three-dimensional space. Nonetheless, it is appropriate to model the shroud as a concentrated mass since it is considerably smaller than the blade. With this modeling, the friction force on the contact surface can be vectorially decomposed into two components in- and out-of- the blisk’s rotation plane. However, the out-of-plane component of the friction force can be effectively ignored since the out-of-plane movement has been assumed negligible in the first place, and hence ideally there will be neither relative movement nor trends of relative movement between two adjacent shrouds. Further, the state of contact between two adjacent shrouds during operation one of the three: stick, slip and separation, depending on the relative motion of the shrouds. In this particular scenario, the normal distance between the shrouds is expressed as

$$d_i=\left(u_{i-1}-u_i\right)\cos \alpha +\left(v_{i-1}-v_i-{\delta }_{\mathrm{s}}\right)\sin \alpha$$

(5)

Hence, the normal force imposed on the right side of the i-th shroud is

$$N_{\mathrm{s},i}=\{\begin{array}{cc}{k}_{\mathrm{s}}{d}_i& d_i>0\\ 0& d_i\le 0\end{array}$$

(6)

For the friction on the solid, dry contact surface, a velocity-dependent model [37] is adopted where friction asymptotically decreases as the relative velocity increases, as

$$\mu \left(s_{\mathrm{r},i}\right)={\mu }_{\mathrm{m}}+\left({\mu }_{\mathrm{s}}-{\mu }_{\mathrm{m}}\right)e^{-\beta \left|s_{\mathrm{r},i}\right|}$$

(7)

where ${\mu }_{\mathrm{m}}$ and ${\mu }_{\mathrm{s}}$ are coefficients of minimum kinetic friction and maximum static friction. $\beta$ is the tuning parameter. $s_{\mathrm{r},i}$ is the relative velocity aligning with the contact surfaces of the i-th and the (i-1)-th shrouds, which can be expressed as

$$s_{r,i}=-\left({\dot{u}}_{i-1}-{\dot{u}}_i\right)\sin \alpha +\left({\dot{v}}_{i-1}-{\dot{v}}_i\right)\cos \alpha$$

(8)

Thus, $F_{\mathrm{s},i}$ governed by the stick-slip model can be written as

$$F_{\mathrm{s},i}=\{\begin{array}{cc}\mu \left(s_{\mathrm{r},i}\right)N_{\mathrm{s},i}\mathrm{sgn}\left(s_{\mathrm{r},i}\right)& s_{\mathrm{r},i}\ne 0\\ -{\mu }_{\mathrm{s}}{N}_{\mathrm{s},i}\le F_{\mathrm{s},i}\le {\mu }_{\mathrm{s}}{N}_{\mathrm{s},i}& s_{\mathrm{r},i}=0\end{array}$$

(9)

2.2. Governing Equations

The governing equations of the shrouded blisk can be derived through the Hamilton’s principle. For the spanwise direction, one has

$$\begin{array}{c}{\rho }_{\mathrm{l}}\left[{\ddot{u}}_i+{\ddot{x}}_{\mathrm{d}}\cos {\theta }_i+{\ddot{y}}_{\mathrm{d}}\sin {\theta }_i-2\mathsf{\Omega }{\dot{v}}_i-{\mathsf{\Omega }}^2{u}_i-{\mathsf{\Omega }}^2\left(x_i+d_r\right)\right]-{\left(EA{u^{\prime }}_i\right)}^{\prime }+c_u{\dot{u}}_i\\ +m_{\mathrm{s}}\left[{\ddot{u}}_i+{\ddot{x}}_{\mathrm{d}}\cos {\theta }_i+{\ddot{y}}_{\mathrm{d}}\sin {\theta }_i-2\mathsf{\Omega }{\dot{v}}_i-{\mathsf{\Omega }}^2{u}_i-{\mathsf{\Omega }}^2\left(x_i+d_{\mathrm{r}}\right)\right]\delta \left(x_i-L\right)=\\ \left(-N_{\mathrm{s},i+1}\cos \alpha +F_{\mathrm{s},i+1}\sin \alpha +N_{\mathrm{s},i}\cos \alpha -F_{\mathrm{s},i}\sin \alpha \right)\delta \left(x_i-L\right)\end{array}$$

(10)

For the flapwise direction, the equation becomes

$$\begin{array}{c}{\rho }_{\mathrm{l}}\left({\ddot{v}}_i-{\ddot{x}}_{\mathrm{d}}\sin {\theta }_i+{\ddot{y}}_{\mathrm{d}}\cos {\theta }_i+2\mathsf{\Omega }{\dot{u}}_i-{\mathsf{\Omega }}^2{v}_i\right)+{\left(EI{v}^{″}\right)}^{″}+c_v{\dot{v}}_i\\ +m_{\mathrm{s}}\left({\ddot{v}}_i-{\ddot{x}}_{\mathrm{d}}\sin {\theta }_i+{\ddot{y}}_{\mathrm{d}}\cos {\theta }_i+2\mathsf{\Omega }{\dot{u}}_i-{\mathsf{\Omega }}^2{v}_i\right)\delta \left(x_i-L\right)=\\ \left(-Q_i-N_{\mathrm{s},i+1}\sin \alpha -F_{\mathrm{s},i+1}\cos \alpha +N_{\mathrm{s},i}\sin \alpha +F_{\mathrm{s},i}\cos \alpha \right)\delta \left(x_i-L\right)\end{array}$$

(11)

One the other hand, the equations of motion for the disk are

$$\begin{array}{c}{\displaystyle \sum _{i=1}^n{\rho }_{\mathrm{l}}\left[\begin{array}{l}{\ddot{x}}_{\mathrm{d}}+{\ddot{u}}_i\cos {\theta }_i-{\ddot{v}}_i\sin {\theta }_i-2\mathsf{\Omega }{\dot{u}}_i\sin {\theta }_i-2\mathsf{\Omega }{\dot{v}}_i\cos {\theta }_i\\ +{\mathsf{\Omega }}^2{v}_i\sin {\theta }_i-{\mathsf{\Omega }}^2{u}_i\cos {\theta }_i-{\mathsf{\Omega }}^2\left(x_i+d_{\mathrm{r}}\right)\cos {\theta }_i\end{array}\right]}\\ +{\displaystyle \sum _{i=1}^n{m}_{\mathrm{s}}\left[\begin{array}{l}{\ddot{x}}_{\mathrm{d}}+{\ddot{u}}_i\cos {\theta }_i-{\ddot{v}}_i\sin {\theta }_i-2\mathsf{\Omega }{\dot{u}}_i\sin {\theta }_i-2\mathsf{\Omega }{\dot{v}}_i\cos {\theta }_i\\ +{\mathsf{\Omega }}^2{v}_i\sin {\theta }_i-{\mathsf{\Omega }}^2{u}_i\cos {\theta }_i-{\mathsf{\Omega }}^2\left(x_i+d_{\mathrm{r}}\right)\cos {\theta }_i\end{array}\right]\delta \left(x_i-L\right)}\\ +\left(m_{\mathrm{d}}{\ddot{x}}_{\mathrm{d}}+c_x{\dot{x}}_{\mathrm{d}}+k_x{x}_{\mathrm{d}}\right)\delta \left(x_i-0\right)=0\end{array}$$

(12)

and

$$\begin{array}{c}{\displaystyle \sum _{i=1}^n{\rho }_{\mathrm{l}}\left[\begin{array}{l}{\ddot{y}}_{\mathrm{d}}+{\ddot{u}}_i\sin {\theta }_i-2\mathsf{\Omega }{\dot{v}}_i\sin {\theta }_i+{\ddot{v}}_i\cos {\theta }_i+2\mathsf{\Omega }{\dot{u}}_i\cos {\theta }_i\\ -{\mathsf{\Omega }}^2{v}_i\cos {\theta }_i-{\mathsf{\Omega }}^2{u}_i\sin {\theta }_i-{\mathsf{\Omega }}^2\left(x_i+d_{\mathrm{r}}\right)\sin {\theta }_i\end{array}\right]}\\ +{\displaystyle \sum _{i=1}^n{m}_{\mathrm{s}}\left[\begin{array}{l}{\ddot{y}}_{\mathrm{d}}+{\ddot{u}}_i\sin {\theta }_i-2\mathsf{\Omega }{\dot{v}}_i\sin {\theta }_i+{\ddot{v}}_i\cos {\theta }_i+2\mathsf{\Omega }{\dot{u}}_i\cos {\theta }_i\\ -{\mathsf{\Omega }}^2{v}_i\cos {\theta }_i-{\mathsf{\Omega }}^2{u}_i\sin {\theta }_i-{\mathsf{\Omega }}^2\left(x_i+d_{\mathrm{r}}\right)\sin {\theta }_i\end{array}\right]\delta \left(x_i-L\right)}\\ +\left(m_{\mathrm{d}}{\ddot{y}}_{\mathrm{d}}+c_y{\dot{y}}_{\mathrm{d}}+k_y{y}_{\mathrm{d}}\right)\delta \left(x_i-0\right)=0\end{array}$$

(13)

where ${\rho }_{\mathrm{l}}$ is the linear mass density of the blades and $m_{\mathrm{s}}$ the mass of the shrouds. EA and EI stand for tensile and bending stiffness of the blade, respectively. $c_u$ and $c_v$ are viscous damping coefficients in the spanwise and flapwise directions, respectively.

In the followings, the Galerkin’s method is adopted to decompose the spanwise and flapwise displacements in the modal space, as

$$\begin{array}{cc}{u}_i\left(x,t\right)={\displaystyle \sum _{j=1}{\mathsf{\Psi }}_j\left(x\right)p_{ji}\left(t\right)},& v_i\left(x,t\right)={\displaystyle \sum _{j=1}{\mathsf{\Phi }}_j\left(x\right)q_{ji}\left(t\right)}\end{array}$$

(14)

where ${\mathsf{\Psi }}_j\left(x\right)$ and ${\mathsf{\Phi }}_j\left(x\right)$ are the j-th mode shapes of cantilever beam in the spanwise and flapwise directions, respectively [38]. $p_{ji}\left(t\right)$ and $q_{ji}\left(t\right)$ are the j-th modal coordinates of the i-th blade, respectively. Considering only the first mode in both two directions, Equations (10)–(13) are reduced

$$\begin{array}{c}{m}_{i,i}{\ddot{p}}_i+m_{i,2n+1}{\ddot{x}}_{\mathrm{d}}+m_{i,2n+2}{\ddot{y}}_{\mathrm{d}}+c_{i,i}{\dot{p}}_i+g_{i,n+i}{\dot{q}}_i+k_{i,i}{p}_i\\ =a_i-N_{\mathrm{s},i+1}\cos \alpha +F_{\mathrm{s},i+1}\sin \alpha +N_{\mathrm{s},i}\cos \alpha -F_{\mathrm{s},i}\sin \alpha ,i=1,2,\dots ,n,\end{array}$$

(15)

$$\begin{array}{c}{m}_{n+i,n+i}{\ddot{q}}_i+m_{n+i,2n+1}{\ddot{x}}_{\mathrm{d}}+m_{n+i,2n+2}{\ddot{y}}_{\mathrm{d}}+g_{n+i,i}{\dot{p}}_i+c_{n+i,n+i}{\dot{q}}_i+k_{n+i,n+i}{q}_i\\ =-Q_i-N_{\mathrm{s},i+1}\sin \alpha -F_{\mathrm{s},i+1}\cos \alpha +N_{\mathrm{s},i}\sin \alpha +F_{\mathrm{s},i}\cos \alpha ,i=1,2,\dots ,n,\end{array}$$

(16)

$$\begin{array}{c}{\displaystyle \sum _{i=1}^n\left(m_{2n+1,i}{\ddot{p}}_i+m_{2n+1,n+i}{\ddot{q}}_i\right)}+m_{2n+1,2n+1}{\ddot{x}}_{\mathrm{d}}+{\displaystyle \sum _{i=1}^n\left(g_{2n+1,i}{\dot{p}}_i+g_{2n+1,n+i}{\dot{q}}_i\right)}\\ +c_{2n+1,2n+1}{\dot{x}}_{\mathrm{d}}+{\displaystyle \sum _{i=1}^n\left(k_{2n+1,i}{p}_i+k_{2n+1,n+i}{q}_i\right)}+k_{2n+1,2n+1}{x}_{\mathrm{d}}=a_{2n+1},\end{array}$$

(17)

$$\begin{array}{c}{\displaystyle \sum _{i=1}^n\left(m_{2n+2,i}{\ddot{p}}_i+m_{2n+2,n+i}{\ddot{q}}_i\right)}+m_{2n+2,2n+2}{\ddot{y}}_{\mathrm{d}}+{\displaystyle \sum _{i=1}^n\left(g_{2n+2,i}{\dot{p}}_i+g_{2n+2,n+i}{\dot{q}}_i\right)}\\ +c_{2n+2,2n+2}{\dot{y}}_{\mathrm{d}}+{\displaystyle \sum _{i=1}^n\left(k_{2n+2,i}{p}_i+k_{2n+2,n+i}{q}_i\right)}+k_{2n+2,2n+2}{y}_{\mathrm{d}}=a_{2n+2}\end{array}$$

(18)

where notations m, c, g, k and a are defined in Equations (A1)–(A5), Appendix B.

The non-dimensional forms of Equations (15)–(18) can be rewritten by introducing the following variables:

$$\begin{array}{l}{\overline{p}}_i=\frac{p_i}{{\delta }_{\mathrm{s}}},{\overline{q}}_i=\frac{q_i}{{\delta }_{\mathrm{s}}},{\overline{x}}_{\mathrm{d}}=\frac{x_{\mathrm{d}}}{{\delta }_{\mathrm{s}}},{\overline{y}}_{\mathrm{d}}=\frac{y_{\mathrm{d}}}{{\delta }_{\mathrm{s}}},\\ {\omega }_i^2=\frac{k_{i,i}}{m_{i,i}},{\omega }_{n+i}^2=\frac{k_{n+i,n+i}}{m_{n+i,n+i}},{\omega }_{2n+1}^2=\frac{k_{2n+1}}{m_{2n+1}},{\omega }_{2n+2}^2=\frac{k_{2n+2}}{m_{2n+2}},\\ {\tau }_q={\omega }_{n+i}t,\stackrel{\ast }{\left(\right)}=\frac{d\left(\right)}{d{\tau }_q},\\ {\overline{a}}_i=\frac{a_i}{{\omega }_{n+i}^2{\delta }_s{m}_{i,i}},{\overline{k}}_{\mathrm{s}}=\frac{k_{\mathrm{s}}}{k_{i,i}},{\overline{Q}}_0=\frac{Q_0}{{\omega }_{n+i}^2{\delta }_{\mathrm{s}}{m}_{n+i,n+i}}\end{array}$$

(19)

Moreover, $d_i$, $N_{\mathrm{s},i}$ and $F_{\mathrm{s},i}$ can be rewritten as

$${\overline{d}}_i=\mathsf{\Psi }\left(L\right)\left({\overline{p}}_{i-1}-{\overline{p}}_i\right)\cos \alpha +\left[\mathsf{\Phi }\left(L\right)\left({\overline{q}}_{i-1}-{\overline{q}}_i\right)-1\right]\sin \alpha$$

(20)

$${\overline{N}}_{\mathrm{s},i}=\{\begin{array}{cc}{\overline{k}}_{\mathrm{s}}{\overline{d}}_i& {\overline{d}}_i>0\\ 0& {\overline{d}}_i\le 0\end{array}$$

(21)

$${\overline{F}}_{\mathrm{s},i}=\{\begin{array}{cc}\mu \left(s_{\mathrm{r},i}\right){\overline{N}}_{\mathrm{s},i}\mathrm{sgn}\left(s_{\mathrm{r},i}\right)& s_{\mathrm{r},i}\ne 0\\ -{\mu }_{\mathrm{s}}{\overline{N}}_{\mathrm{s},i}\le {\overline{F}}_{\mathrm{s},i}\le {\mu }_{\mathrm{s}}{\overline{N}}_{\mathrm{s},i}& s_{\mathrm{r},i}=0\end{array}$$

(22)

Using Equation (19), the non-dimensional form of the governing equations are

$$\begin{array}{c}{\stackrel{**}{\overline{p}}}_i+{\overline{m}}_{i,2n+1}{\stackrel{\ast \ast }{\overline{x}}}_{\mathrm{d}}+{\overline{m}}_{i,2n+2}{\stackrel{\ast \ast }{\overline{y}}}_{\mathrm{d}}+{\overline{c}}_{i,i}{\stackrel{\ast }{\overline{p}}}_i+{\overline{g}}_{i,n+i}{\stackrel{\ast }{\overline{q}}}_i+{\overline{\omega }}_{i,i}^2{\overline{p}}_i={\overline{a}}_i\\ +{\overline{m}}_{qp}\left(-{\overline{N}}_{s,i+1}\cos \alpha +{\overline{F}}_{\mathrm{s},i+1}\sin \alpha +{\overline{N}}_{\mathrm{s},i}\cos \alpha -{\overline{F}}_{\mathrm{s},i}\sin \alpha \right),i=1,2,\dots ,n\end{array}$$

(23)

$$\begin{array}{c}{\stackrel{\ast \ast }{\overline{q}}}_i+{\overline{m}}_{n+i,2n+1}{\stackrel{\ast \ast }{\overline{x}}}_{\mathrm{d}}+{\overline{m}}_{n+i,2n+2}{\stackrel{\ast \ast }{\overline{y}}}_{\mathrm{d}}+{\overline{g}}_{n+i,i}{\stackrel{\ast }{\overline{p}}}_i+{\overline{c}}_{n+i,n+i}{\stackrel{\ast }{\overline{q}}}_i+{\overline{q}}_i=-{\overline{Q}}_i\\ -{\overline{N}}_{\mathrm{s},i+1}\sin \alpha -{\overline{F}}_{\mathrm{s},i+1}\cos \alpha +{\overline{N}}_{\mathrm{s},i}\sin \alpha +{\overline{F}}_{\mathrm{s},i}\cos \alpha ,i=1,2,\dots ,n\end{array}$$

(24)

$$\begin{array}{c}{\stackrel{\ast \ast }{\overline{x}}}_{\mathrm{d}}+{\displaystyle \sum _{i=1}^n\left({\overline{m}}_{2n+1,i}{\stackrel{\ast \ast }{\overline{p}}}_i+{\overline{m}}_{2n+1,n+i}{\stackrel{\ast \ast }{\overline{q}}}_i\right)}+{\displaystyle \sum _{i=1}^n\left({\overline{g}}_{2n+1,i}{\stackrel{\ast }{\overline{p}}}_i+{\overline{g}}_{2n+1,n+i}{\stackrel{\ast }{\overline{q}}}_i\right)}\\ +{\overline{c}}_{2n+1,2n+1}{\stackrel{\ast }{\overline{x}}}_{\mathrm{d}}+{\displaystyle \sum _{i=1}^n\left({\overline{k}}_{2n+1,i}{\overline{p}}_i+{\overline{k}}_{2n+1,n+i}{\overline{q}}_i\right)}+{\overline{\omega }}_{2n+1}^2{\overline{x}}_{\mathrm{d}}={\overline{a}}_{2n+1}\end{array}$$

(25)

$$\begin{array}{c}{\stackrel{\ast \ast }{\overline{y}}}_{\mathrm{d}}+{\displaystyle \sum _{i=1}^n\left({\overline{m}}_{2n+2,i}{\stackrel{\ast \ast }{\overline{p}}}_i+{\overline{m}}_{2n+2,n+i}{\stackrel{\ast \ast }{\overline{q}}}_i\right)}+{\displaystyle \sum _{i=1}^n\left({\overline{g}}_{2n+2,i}{\stackrel{\ast }{\overline{p}}}_i+{\overline{g}}_{2n+2,n+i}{\stackrel{\ast }{\overline{q}}}_i\right)}\\ +{\overline{c}}_{2n+2,2n+2}{\stackrel{\ast }{\overline{y}}}_{\mathrm{d}}+{\displaystyle \sum _{i=1}^n\left({\overline{k}}_{2n+2,i}{p}_i+{\overline{k}}_{2n+2,n+i}{q}_i\right)}+{\overline{\omega }}_{2n+2}^2{\overline{y}}_{\mathrm{d}}={\overline{a}}_{2n+2}\end{array}$$

(26)

The coefficients of Equations (23)–(26) are provided in Equations (A6)–(A9), Appendix B. A numerical simulation method implementing the Runge-Kutta method with variable steps is developed using MATLAB to solve the equations of motion of the shrouded blades in this study considering the nonsmooth contact and friction behavior. Given the significance of the contact behavior between adjacent shrouds on the blade dynamics, the bisection method is used to capture the precise times where the contact switches from one state to another in Appendix D. The transition of the contact states (i.e., stick, slip and separation) of adjacent shrouds and corresponding friction forces are provided in Appendix C.

3. Nonlinear Dynamics and Vibration Localization of Shrouded Blisk

This section includes the results of nonlinear dynamics study and vibration localization analysis of a shrouded, tuned blisk. The shrouded blisk with twelve blades is used for numerical investigation. The main parameters of the shrouded blisk [20,39] are shown in

Table 1. Gross parameters of the shrouded blisk.

Notation	Value	Notation	Value
$c_u,c_v$	1·N·s·m−1	$c_x,c_y$	5 N·s·m−1
$d_{\mathrm{r}}$	0.5 m	EI	1.45 × 103 N·m2
EA	1.74 × 108 N	$k_x,k_y$	1.0 × 107 N·m−1
L	0.5 m	$m_{\mathrm{d}}$	14.04 kg
$m_{\mathrm{s}}$	0.135 kg	$Q_0$	50 N
$\alpha$	$\frac{\pi }{3}$ rad	$\beta$	5 s·m−1
${\delta }_{\mathrm{s}}$	2.0 × 10−5 m	$\lambda$	1
${\mu }_{\mathrm{s}}$	0.5	${\mu }_{\mathrm{m}}$	0.3
${\rho }_{\mathrm{l}}$	6.75 kg·m−1	Ω	6000 r·min−1

.

A finite element analysis containing two shrouded blades was carried out using ANSYS software to determine the normal contact stiffness ($k_{\mathrm{s}}$, caused by the elastic deformation of the shroud) [5] at various rotation speed in this paper. As shown in Figure 3a, contact pair was used to simulate the connecting force between adjacent shrouds, the solid185-type elements were used to model the structural parts of the blades and shrouds. The parameters of the model are presented in Appendix E. Then, the contact stiffness was calculated based on the elastic deformations of the shrouds with respect to the reaction forces of the blades in the normal direction of the shroud. Afterwards, the contact stiffness was fitted to generate $k_{\mathrm{s}}$ as a function of $\mathsf{\Omega }$, which is depicted in Figure 3b.

The vibration response can be solved once the normal contact stiffness determined. The steady-state responses with different rotation speeds are illustrated. Then, the localization of vibration due to the contact and friction of shrouds is studied.

3.1. Steady-State Response, Primary Resonance and Beat

The rotation of the disk creates centrifugal force and subsequently the effect of spin softening through which the blade’s axial stiffness is reduced, as well as changes the state of contact between adjacent shrouds as shown in Figure 3b. In the following, Equations (23)–(26) are solved with various rotation speeds to demonstrate how steady-state response is affected by these effects. Since the disk is structurally tuned, the response will be presented hereafter with only one sector including the first shrouded blade (i.e., i = 1 in Figure 2b). To facilitate the discussion in next subsections, two forcing parameters are introduced to represent the “net” contact and friction:

$$N_1={\overline{N}}_{\mathrm{s},1}-{\overline{N}}_{\mathrm{s},2},F_1={\overline{F}}_{\mathrm{s},1}-{\overline{F}}_{\mathrm{s},2}.$$

The frequency spectrum of the flapwise motion ${\overline{q}}_1$ is presented in Figure 4 where peaks A, B and C mark the frequency of the aerodynamic excitation ($f_{\mathrm{Q}}$) and the frequencies of the first flapwise mode ($f_{\mathrm{q}}$) and the first spanwise mode ($f_{\mathrm{p}}$). When $\mathsf{\Omega }$ escalates to 4300 r/min, primary resonance (point D) in the blade motion occurs as $f_{\mathrm{Q}}$ approaches $f_{\mathrm{q}}$. This demonstrates the significant influence of the aerodynamic force through the flapwise movement of the blade. The magnitude of the resonance decays as the two frequencies departs, which in turn reduces the effect of contact and friction on the blade shroud.

Further examinations are carried out focusing on the vibration response under four selected rotation speeds, i.e., 4300, 4750 and 9000 r/min, and the results are presented in Figure 5. As shown in Figure 5a–d, the spanwise motion ${\overline{p}}_1$ appears close to periodic, whereas the flapwise ${\overline{q}}_1$ is not. Further, the primary resonance in the flapwise direction occurs when $\mathsf{\Omega }=4300$ r/min. Based on the fast Fourier transform analysis, ${\overline{p}}_1$ is dominated by $f_{\mathrm{p}}$ in term of amplitude, and is slightly modulated by weak disturbances that carry low frequencies such as $f_{\mathrm{q}}$ and $f_{\mathrm{Q}}$. For the flapwise motion ${\overline{q}}_1$, all of the frequencies are incommensurable whose amplitudes are in the same scale, Figure 5f. It should be noted that free vibration will be produced by the initial disturbance during each contact along with the forced response resulted from the aerodynamic force on the blade. To be specific $f_{\mathrm{p}}$ and $f_{\mathrm{q}}$ are frequencies of the free vibration, whereas $f_{\mathrm{Q}}$ and its odd multiples are frequencies related to the forced vibration.

As far as the contact and friction on the shroud are concerned, both $N_1$ and $F_1$ are non-smooth periodic and spectrally governed by $f_{\mathrm{Q}}$ and its odd multiples (e.g., $3f_{\mathrm{Q}}$ and $5f_{\mathrm{Q}}$) that can be identified in Figure 5g. Further, the direction of the friction force reverses every time the normal force reaches its maximum value, Figure 5h. Both $N_1$ and $F_1$ change periodically and vanish when the shroud is completely separated from the other shroud, or the forces by the two adjacent shrouds cancel with each other during contact.

As $\mathsf{\Omega }$ arrives at 4750 r/min, a slight difference of 0.02 Hz between $f_{\mathrm{Q}}$ and $f_{\mathrm{q}}$ is found in Figure 6a, and the phenomenon of beat is observed from Figure 6b. The peak-to-peak amplitude of ${\overline{q}}_1$ enlarges owing to the presence of beat, and so do the contact and friction forces between adjacent shrouds. Further, the net contact and friction forces appear periodic as shown in Figure 6c,d. Moreover, localization of vibration is likely to take places in this case which will be presented in the next subsection.

Finally, when the blade is accelerated to $\mathsf{\Omega }=9000$ r/min, both ${\overline{p}}_1$ and ${\overline{q}}_1$ become quasi-periodic, Figure 7a–d. The frequency $f_{\mathrm{Q}}$ can hardly be identified from Figure 7e,f. As a matter of fact, the normal contact stiffness becomes smaller at high rotation speed (c.f. Figure 3) and the actions of contact and friction are both weak, which explains the amplitudes associated with $f_{\mathrm{Q}}$ are minor in the spectra of ${\overline{p}}_1$ and ${\overline{q}}_1$. Further, as shown in Figure 7g, the spectra of normal and friction forces both contain the identified $f_{\mathrm{N}}$, the aerodynamic frequency $f_{\mathrm{Q}}$ and its odd multiples as well as combined frequencies of $3f_{\mathrm{N}}+f_{\mathrm{Q}}$, $5f_{\mathrm{N}}+f_{\mathrm{Q}}$, $3f_{\mathrm{N}}+3f_{\mathrm{Q}}$ and $5f_{\mathrm{N}}+3f_{\mathrm{Q}}$, etc. This clearly shows the nonlinear yet complicated nature of the vibration induced by the action of the normal and friction forces.

In this study it is found that the primary resonance, beat and quasi-periodic vibration appear successively due to the contact and friction between shrouds with the increase of rotation speed. On the whole, the blisk structures can be designed more reasonably based on the in-depth understandings of these phenomena and are more likely to perform better under the circumstance of vibration localization.

3.2. Localization of Vibration of Blades

Factor of Localization of Vibration

In this section, different form existing publications, we discover a particular type of vibration localization initiated by the contact and friction effect of the shrouds in the form of asymmetric vibration in the tuned bladed, cyclically symmetric disk, then reveal and discuss its influence on the vibration of the blisk. To demonstrate the unevenness in the blade deformations, the maximal displacements of each blade are selected over the period of the vibration and compared individually with their arithmetic means. One may define the following factors of localizations of vibration in the flapwise and spanwise directions, i.e.,

$$L_p={\displaystyle \sum _{i=1}^n\frac{\max \left({\overline{p}}_i\left(t\right)\right)-{\overline{p}}_{\mathrm{ave}}}{{\overline{p}}_{\mathrm{ave}}}},L_q={\displaystyle \sum _{i=1}^n\frac{\max \left({\overline{q}}_i\left(t\right)\right)-{\overline{q}}_{\mathrm{ave}}}{{\overline{q}}_{\mathrm{ave}}}},\forall t>0$$

(27)

where

$${\overline{p}}_{\mathrm{ave}}=\frac{1}{n}{\displaystyle \sum _{i=1}^n\max \left({\overline{p}}_i\left(t\right)\right),}{\overline{q}}_{\mathrm{ave}}=\frac{1}{n}{\displaystyle \sum _{i=1}^n\max \left({\overline{q}}_i\left(t\right)\right)}$$

(28)

Note that these factors will be zero if the displacements are constant or circumferentially periodic among the blades.

Next, the maximum displacements are presented in Figure 8 for three cases: the general vibration when $\mathsf{\Omega }=6000$ r/min, the case of beat when $\mathsf{\Omega }=4750$ r/min and the case of blades free of contact and friction when $\mathsf{\Omega }=4750$ r/min. The factors of the localization of vibration are tabulated in

Table 2. Factors of localization of vibration in different cases.

$\mathit{\xi }$	${\mathit{L}}_{\mathit{p}}$	${\mathit{L}}_{\mathit{q}}$
General case, $\mathsf{\Omega }=6000$ r/min	0.0	0.0
The beat case, $\mathsf{\Omega }=4750$ r/min	−3.902 × 10−5	−1.197 × 10−4
Contact- and friction free, $\mathsf{\Omega }=4750$ r/min	0.0	0.0

. Factors of localization of vibration in different cases for these three cases. Through examining Figure 8a–c one can see the uneven distribution of maximum displacements as the result of the localization of vibration resulted from the beat, as well as the contact and friction between the perfectly tuned blades. The level of localization is found much stronger in the flapwise direction than in the spanwise direction due to the high tensile stiffness of the blades.

3.3. Effect of Coupling Stiffness on the Localization of Vibration

The interaction among the blades enables the transfer of energy across the disk parts, and influence the contact and friction on the shroud. In the present investigation, such actions can be considered through adopting the cross-coupling stiffness ${\overline{k}}_{2n+1,i}$, ${\overline{k}}_{2n+1,n+i}$, ${\overline{k}}_{2n+2,i}$, ${\overline{k}}_{2n+2,n+i}$ in Equations (25) and (26). To demonstrate how the localization is affected by the cross coupling, Equations (23) through (26) are solved with various cross-coupling stiffness. For convenience of comparison, these coefficients are amplified by the identical factor denoted by $\xi$.

The maximum displacements of the blades are shown in Figure 9 with varying $\xi$, and the factors of localization are tabulated in

Table 3. Factors of the localization of vibration. $\mathsf{\Omega }=4750$ r/min.

$\mathit{\xi }$	${\mathit{L}}_{\mathit{p}}$	${\mathit{L}}_{\mathit{q}}$
0.01	−5.024 × 10−5	1.190 × 10−5
0.1	2.003 × 10−5	−3.704 × 10−5
1	−3.902 × 10−5	−1.197 × 10−5
10	2.429 × 10−5	0.0017
100	6.471 × 10−5	0.0063

. Based on these results, the factors of localization increase significantly relative to $\xi =1$ when the cross-coupling is strong, with $L_q$ being more sensitive than $L_p$. This can be explained by the fact that $L_p$ is closely related to the spanwise displacements which are governed mainly by the centrifugal force, and hence is insensitive to the cross-coupling effect. Conversely, $L_q$ is affected primarily by the contact forces, and hence more prone to the change of the cross-coupling effect.

Meanwhile, as shown in Figure 10, the increasing $\xi$ enlarges the normal and tangential disk-blade interactions at the root of the i-th blade (defined as $F_{\mathrm{N},i}\left(t\right)$ and $F_{\mathrm{T},i}\left(t\right)$ in Equations (A10) and (A11)), which boosts the energy transferal from the disk to the blades. Consequently, the normal and friction forces on the shroud interfaces will be enhanced, as well as the localization of vibration caused by these forces.

3.4. Power of Blade Vibration under Localization of Vibration

Next, the effect of the localization of vibration on the blade vibration is studied considering varying normal contact stiffness and gap between the adjacent shrouds. To measure the level of the blade motions, a generalized, average power over a period of steady-state vibration is defined

$$\begin{array}{cc}{E}_i=\frac{1}{T}{\displaystyle {\int }_0^T\left[p_i^2\left(t\right)+q_i^2\left(t\right)\right]}dt,& i=1,2,\dots ,n\end{array}$$

(29)

Using various normal contact stiffness between adjacent shrouds ($k_{\mathrm{s}}$), the generalized powers of all twelve blade motions are depicted in Figure 11a,b under the localization of vibration and compared to the ones from regular, non-localized vibration. As shown in Figure 11a,b, $E_i$ appear confined to two belts of narrowly spaced curves. These belts are separated under the localization of vibration throughout the entire range of $k_{\mathrm{s}}$, which is attributable to the distorted pattern of maximum displacements of the blades (c.f., Figure 8b). When it happens, the localization of vibration tends to produce larger difference in displacements among the blades, thus creating a unique pattern of two groups of $E_i$, each carrying nearly identical powers of vibration. Moreover, it is observed in Figure 11a that they both decrease with $k_{\mathrm{s}}$. As shown previously in Table 3, the localization of vibration can be enhanced by the increase of the normal contact stiffness, which further enlarges the difference in displacements as well as those $E_i$ among the blades.

In the non-localization case, the $E_i$ drop at first and then escalate when $k_{\mathrm{s}}$ increases, Figure 11b. On one hand, more energies are dissipated as friction on the shroud interface grows with the increasing contact forces, which reduces $E_i$. On the other hand, the distance between the adjacent shrouds decreases when $k_{\mathrm{s}}$ is large enough (see Figure 11c), which reduces the friction force (see Figure 11d) and consequently, increases the level of $E_i$. Thus, the powers $E_i$ decreases unanimously with the increasing normal contact stiffness, Figure 11a. In contrast, these distances are amplified when the localization of vibration is encountered.

4. Conclusions

The nonlinear dynamics and the vibration localization of a whole rotating blisk with the shrouded friction contacts are presented in this paper. A continuous parametric model of a shrouded blisk with short blades of small installation angles is established, and the transition boundaries between different contact statuses (i.e., stick, slip and separation) of the whole shrouded blisk are determined based on the relative motion status of the blades. Using this model, the nonlinear vibration responses with variable rotation speed are analyzed, the particular type of vibration localization initiated by the contact and friction effect of the shrouds is discovered, and its influence on the vibration of the blisk is revealed and discussed. The followings summarize the present study:

(1)

With the increase of rotation speed, the primary resonance, beat and quasi-periodic vibration appear successively due to the contact and friction between shrouds.

(2)

The localization of vibration occurs in tuned shrouded blisks only when the beat exists, i.e., frequency of free vibration slightly differs from the frequency of the forced vibration caused by the contact and friction of shrouds. Further, the level of such localization increases as the blade-disk coupling becomes strong.

(3)

The power of the blade motion drops consistently with the contact stiffness of the shrouds in the localization of vibration, which is different from the non-localization case.