A Prediction Method with Altering Equivalent Stiffness for Damping Evaluation of Shrouded Bladed Disk Dynamic Systems

Abstract

In turbomachinery bladed disks are typical cyclic symmetric structures where high-cycle fatigue of the blades can easily occur. Increasing the shroud structure is an effective blade vibration reduction method using the dry friction of the shroud contact surface to dissipate vibration energy. During one vibration cycle of the blade, the contact surface of the shroud may have one or several states of stick, slip, and separation due to the different amplitude of blade, which affects not only the damping, but also the stiffness of the system. This paper proposes a method to analyze the damping characteristics of blades with shroud considering the change of equivalent contact stiffness of the shroud by using the finite element method. With a 2D contact model, this method calculates the equivalent damping and obtains the damping characteristic curve to evaluate the damping of shroud. During the design phase, the objective is that when resonance occurs, the shroud can provide the optimal damping ratio under the allowable vibration stress of the blade. The parameter sensitivity of the damping characteristic for one bladed disk with parallel shroud is investigated. It is shown that vibration phase angle, contact stiffness, contact force, friction coefficient significantly influence the damping characteristics.

1. Introduction

Turbine rotor blades in aero-engines are subjected to high centrifugal loads, static stresses caused by thermal loads, and vibration stresses caused by varying loads, which are likely to cause high-cycle fatigue (HCF) of the blades. This makes it necessary to reduce the static stress and control the vibration stress [1]. During the blade operation, the excitation frequency range is wide, and there are many engine orders (EOs), and it is difficult to avoid all resonances. The dry friction damping technology, which dissipates vibration energy through con-tact surface friction, is an effective means to reduce the vibration stress of rotor blades at critical positions [2]. Low-pressure turbine blades often adopt a shroud structure, which can not only reduce gas leakage to improve aerodynamic efficiency, but increase blade stiffness and damping to reduce blade vibration.

Two typical shroud configurations are zigzag and slant. The desired contact pressure distribution can be obtained through pre-twisting or interference during assembly. At the same time, centrifugal force will also change contact pressure during rotation. This paper adopts a slant shroud and each blade has the half-tenon. Two blades are assembled in pairs in the same mortise, as Figure 1, where the initial normal contact force between shrouds can be obtained through interference assembly and centrifugal force. The assessment of shroud depends on the grasp of the dynamics and contact of the bladed disk structure [3].

The contact model describes the interaction of an elastic object (usually metal) in dry friction contact, and its relationship with the relative motion of contact surfaces [4,5]. Currently, there are multiple contact models: According to the dimension of the relative tangential displacement of the contact surface, it can be divided into 1D model [6,7], 2Dmodel [8,9]; Model with constant [6,8] or variable normal force [7,9]. According to whether to consider local sliding and local stick, they are divided into macro-slip and micro-slip models [10]. Using the developed contact model, the dynamic characteristics of bladed disk can be calculated. The direct time integration (DTI) method is very time-consuming. In 2008, Phadke and Berger [11] used the finite element tool to calculate the steady-state response of the blade with underplatform dampers by DTI. The response calculation at one excitation frequency required 40 h at that time. Therefore, researchers have developed faster numerical calculation methods. Considering that the vibration response of bladed disk is periodic, the harmonic balance method (HBM) [12] expands the variables (displacement and force) into Fourier series, combined with the Galerkin method, and the coefficient of the corresponding harmonic term of the equation residual is zero. In this way, the equation of motion is converted from a differential equation in the time domain to an algebraic equation of the Fourier expansion coefficients of the displacement, excitation, and friction, and excitation frequency in the frequency domain, and then this nonlinear equation is solved in the frequency domain. In order to describe contact more effectively, Cameron and Griffin [13] proposed alternating frequency-time (AFT), where after converting the frequency domain displacement to the time domain, the contact forces are solved in the time domain, and then they are converted back to the frequency domain. In order to further reduce the calculation time, the researchers proposed reduced order models (ROM) suitable for different situations to calculate the nonlinear forced response [14,15,16], such as component mode synthesis (CMS) [17] and frequency response function (FRF) [18]. Some researchers used non-linear modal analysis to deal with the dynamics with dry friction from the perspective of non-conservative systems [19,20]. In addition, some other interesting vibration phenomena, such as random vibration by an equivalent linearization method [21] and rub-dynamic characteristics by macroslip friction model [22] of bladed disks were studied.

The amplitude-frequency response curve is often used to evaluate the damping effect of shroud [23]. However, the energy dissipated by shroud contact has been seldom studied. This paper proposes a method for evaluating the damping of turbine blade shroud by the energy method (EM). The goal in design stage is to provide greater damping ratio under the allowable vibration stress of the blade. In this method, the equivalent damping ratio is calculated by the proportional relationship between friction energy consumption and the max vibration kinetic energy of the bladed disk. The influence of the shroud at the blade tip on the blade stiffness is considered by the inter-shroud spring. In some vibration modes, the separation phenomenon is easy to occur in the contact state between the slant shrouds. In this paper, the effect of shroud separation on blade stiffness is considered.

In this paper, a method for analyzing the damping characteristics of shrouded turbine blades considering the variation of the contact equivalent stiffness of the shroud is presented. Firstly, the contact model, especially the calculation of the equivalent stiffness of the contact surface, is introduced. Then, the tangential relative displacement of the contact surface is determined by considering the modal analysis results with the equivalent stiffness in three directions. The hysteresis loops of the friction force and the relative displacement are obtained by the contact model, and the equivalent damping of different contact models is calculated using the energy method. Then, the curve of the equivalent stiffness with the amplitude is established, and the influence of the separation on the equivalent stiffness is considered. Finally, the suggestions in the design stage are given, and the parameter sensitivity analysis of damping characteristics is carried out from four aspects: normal and tangential vibration phase angle, contact stiffness, initial contact force and friction coefficient.

2. Dynamic Model of Shrouded Bladed Disks with Contact

2.1. Equation of Motion

In this paper, a disk with 94 slant shrouded blades is considered. Two blades are assembled in pairs in the same mortise of the disk through two half-mortise structures. By cyclic symmetric boundary conditions, the dynamics of the tuned bladed disk can be analyzed by a single sector [24,25], which greatly reduces the amount of calculation. The bladed disk structure is divided into 47 sectors, and each sector contains a pair of blades assembled in the same mortise. The equation of motion (EoM) of one sector in a finite element model with contact is:

$$M\ddot{u}+C\dot{u}+Ku=f_e-f_c$$

(1)

where M, C and K are the mass, damping and stiffness matrices, respectively. f_e is the excitation vector. f_c is the contact force vector. f_c is not zero only in the degree of freedom of the contact surface. This paper assumes that the damping provided by the root is small [26], and only the damping provided by the shroud is analyzed.

2.2. Contact Model

The contact law describes the relationship between the contact force of the contact surface, that is, the tangential friction force and the normal force, and the relative displacement of the contact surface, which is directly related to the friction energy consumption of the system. In this paper, a pair of contact nodes is used to model the contact between shrouds. In the contact local coordinate system, the 3D relative motion of the contact node pair can be decomposed into the normal relative motion out of the plane and the two mutually perpendicular relative motions in the plane. However, this paper focuses on the evaluation method for damping, therefore the simplified contact model is adopted to speed up the calculation.

The local coordinate system of the shroud contact is shown as Figure 2. Because of the 2D contact model, it is necessary to specify the normal and tangential directions of the contact surface. The normal direction of the surface z_n is the normal direction of the contact model, while x_t and y_t are the two tangential directions.

The 1D contact model assumes that the normal force is constant, and only one tangential relative motion at the contact surface is considered, as Figure 3. Based on 1D, the 2D contact model considers the change of normal force, as Figure 4. See [9] and Appendix A for 1D and Appendix B for 2D for the detailed contact law and the transition criteria of different contact states (stick, slip and separation).

2.3. Equivalent Stiffness of the Shroud Contact

The contact has an effect on the damping and stiffness of the system. The equivalent stiffness between the contact surfaces of the blade shroud is needed to calculate the dynamic characteristics of the bladed disk system. There is a large gap between the adjacent pairs so their shrouds will not contact during vibration, while two shrouds in the same pair can contact and the contact pressure can be adjusted by interference fitting. The equivalent stiffness is calculated by nonlinear static analysis with contact. The cyclic symmetry interface and the contact surfaces including the shrouds, the blade tenon and disk mortise are as Figure 5.

The nonlinear static analysis with the circumferential (y) and axial (z) displacement constraints in the global cylindrical coordinate system at some nodes of the disk, the slant constraints at the side of the blade tenon and the rotating speed to consider the compression effect of centrifugal force on the shroud. The equivalent stiffness of contact depends on the normal pressure and the shroud deformation in the analysis results.

The contact normal force N₀ first. Because the contact stress results of contact nonlinear static analysis depend on the mesh density and contact algorithm, the displacement results are relatively more robust. In order to obtain more accurate normal force, multiple linear static analyses are iteratively carried out by applying different pressure loads on the contact surface of the shroud until the displacement result is consistent with that of the contact nonlinear static analysis. The normal force on the whole surface at this time which is the product of pressure load and the area is the equivalent normal force N₀.

The deformation of the contact surface δ follows. Based on two static analyses with contact and without contact of the shrouds, the deformation of the shrouds under the two states is obtained. Taking the left shroud as an example, the deformed shape of the shroud in two states is drawn as Figure 6. The difference in normal deformation of the left shroud in contact with the free state is δ_l = |h_fr − h_co|, similarly δ_r for the right shroud-.

The nonlinear static analysis with the circumferential (y) and axial (z) displacement constraints in the global cylindrical coordinate system at some nodes of the disk, the slant constraints at the side of the blade tenon and the rotating speed to consider the compression effect of centrifugal force on the shroud. The equivalent stiffness of contact depends on the normal pressure and the shroud deformation in the analysis results.

The equivalent normal contact stiffness is calculated by:

$$k_{\mathrm{n}}=\frac{N_0}{{\delta }_{\mathrm{l}}+{\delta }_{\mathrm{r}}}$$

(2)

The tangential contact stiffness is obtained as follows: for the sphere-plate contact [27], the tangential and normal contact stiffness satisfies:

$$k_{\mathrm{t}}=\lambda k_{\mathrm{n}}=\frac{2(1-\nu )}{2-\nu }{k}_{\mathrm{n}}$$

(3)

where ν is Poisson’s ratio. For the plate-plate contact in this paper, λ is slightly smaller, and take λ = 1/3. The results of the normal force and the contact stiffness are shown as

Table 1. Equivalent normal load and stiffness.

Normal Force N0 (N)	Normal Contact Stiffness kn (×106 N/m)	Tangential Contact Stiffness kt (×106 N/m)
7.145	2.06	0.69

.

3. Equivalent Damping Analysis

3.1. Calculation Method of Equivalent Damping Ratio

The contact has an effect on the damping and stiffness of the system. The equivalent stiffness between the contact surfaces of the blade shroud is needed to calculate the dynamic characteristics of the bladed disk system. The mechanism of dry friction damping is that the friction generated by the relative motion between the contact surfaces dissipates the vibration energy. The ratio of dissipated energy to vibration energy can be used to evaluate the damping effect. The energy dissipated by friction in one vibration cycle can be obtained from the area enclosed by the hysteresis loop as:

$$E_{\mathrm{f}}={\displaystyle \oint f_{c,t}(x)\mathrm{d}x}$$

(4)

Loss factor is often used to compare the damping effect of engineering materials [28], which can be expressed as the proportion of the dissipated energy of each vibration cycle to the vibration energy of the system, i.e.:

$$\eta =\frac{E_f}{2\pi E}\simeq 2\zeta$$

(5)

The vibration energy of the system can be approximated by the max vibration kinetic energy in a period, and the equivalent damping ratio provided by the shroud is approximately calculated by:

$$\zeta =\frac{E_f}{4\pi E_k}$$

(6)

3.2. Parameters for Damping Evaluation

The matrix 27 elements are established between several pairs of nodes on the contact surface of the shroud, and their stiffness value are the contact stiffness, and the prestress modal analysis is carried out, which can be directly calculated with the help of commercial software, such as ANSYS^®. In the static analysis, the contact surface load is applied on the shroud contact surfaces, and the element stiffness of matrix 27 is set to 0; in the modal analysis, the surface load is deleted, and the stiffness of matrix 27 is set to the equivalent contact stiffness calculated above.

The modes of interest are determined by Campbell. In this section, the equivalent damping ratios are calculated for three nodal diameters (ND), first bend (1B) mode and five ND first edgewise (1E) mode.

Designers hope that the shroud can play a significant damping role in resonance. The steady-state forced response of bladed disk remains periodic under EO traveling excitation. Assuming that the response is simple harmonic and traveling, the motion of shrouds can be expressed as:

$$\{\begin{array}{l}{u}_{\mathrm{t}\mathrm{l}}(t)=X_{\mathrm{t}\mathrm{l}}\sin (\omega t)\\ u_{\mathrm{n}\mathrm{l}}(t)=X_{\mathrm{n}\mathrm{l}}\sin (\omega t)\end{array}$$

(7)

$$\{\begin{array}{l}{u}_{\mathrm{t}\mathrm{r}}(t)=X_{\mathrm{t}\mathrm{r}}\sin (\omega t+{\phi }_{\mathrm{k}})\\ u_{\mathrm{n}\mathrm{r}}(t)=X_{\mathrm{n}\mathrm{r}}\sin (\omega t+{\phi }_{\mathrm{k}})\end{array}$$

(8)

where the subscript t means tangential, n means normal, l represents left shroud, and r represents right shroud. φ_k is Inter-Blade Phase Angle (IBPA), and ω is the natural angular frequency.

Using the full expansion of cyclic symmetry modal analysis, the sector with the largest amplitude is selected to extract the parameters needed for damping analysis. Parameters for damping evaluation are given in

Table 2. Modal parameters for damping evaluation.

Parameters	3ND 1B	5ND 1E
Normal amplitude of the left shroud surface Xtl/mm	45.65	4.88
Tangential amplitude of the left shroud surface Xnl/mm	−1.00	45.41
Normal amplitude of the right shroud surface Xtr/mm	45.56	4.90
Tangential amplitude of the right shroud surface Xnr/mm	3.20	45.57
Max vibration stress of blade σref/MPa	5268	13,056
Max vibration kinetic energy of the selected sector Ek,ref/J	132	616
IBPA φk/°	12	19

.

where max vibration kinetic energy of the selected sector E_k,ref is:

$$E_{\mathrm{k},\mathrm{ref}}=\frac{\max \left\{A_i^2,i=1,2,\dots ,47\right\}}{{\displaystyle \sum _{i=1}^{47}{A}_i^2}}{E}_{\mathrm{k},\mathrm{wh}}$$

(9)

in which E_k,wh is the result of the full expansion. IBPA φ_k = 2π/N_b, where k is ND number, and N_b is the number of blades.

In the contact model, the relative displacement of the shrouds is required as:

$$\begin{array}{l}\Delta u_{\mathrm{t}}(t)=u_{\mathrm{tr}}-u_{\mathrm{tl}}=\sqrt{X_{\mathrm{t}\mathrm{l}}^2+X_{\mathrm{t}\mathrm{r}}^2-2X_{\mathrm{t}\mathrm{l}}{X}_{\mathrm{t}\mathrm{r}}\cos ({\phi }_{\mathrm{k}})}\sin (\omega t+{\theta }_{\mathrm{t}})=A\sin (\omega t+{\theta }_{\mathrm{t}})\\ \Delta u_{\mathrm{n}}(t)=u_{\mathrm{nr}}-u_{\mathrm{nl}}=\sqrt{X_{\mathrm{n}\mathrm{l}}^2+X_{\mathrm{n}\mathrm{r}}^2-2X_{\mathrm{n}\mathrm{l}}{X}_{\mathrm{n}\mathrm{r}}\cos ({\phi }_{\mathrm{k}})}\sin (\omega t+{\theta }_{\mathrm{t}})=B\sin (\omega t+{\theta }_{\mathrm{n}})\end{array}$$

(10)

1D contact model only needs tangential while in 2D model, not only the tangential and normal displacements but also the phase angle are needed. The relative displacements can be:

$$\begin{array}{l}\Delta u_{\mathrm{t}}(t)=A\sin (\omega t)\\ \Delta u_{\mathrm{n}}(t)=B\sin (\omega t+\theta )\end{array}$$

(11)

where A is the tangential amplitude, B is the normal amplitude and θ = θ_n − θ_t is the phase angle between normal and tangential motion, as

Table 3. Modal parameters of relative motion for damping evaluation.

Parameters	3ND 1B	5ND 1E
Normal amplitude Aref/mm	9.13	1.63
Tangential amplitude Bref/mm	4.19	15.13
phase angle θ/°	−92	0

.

3.3. Damping Characteristic Curve

On a given mode, different amplitudes will correspond to different vibration stresses. Under small vibration amplitude conditions, the following relationships are satisfied between the modal and real vibration amplitude of relative motion, vibration stress, and max kinetic energy:

$$\frac{X}{X_{\mathrm{ref}}}=\frac{\sigma }{{\sigma }_{\mathrm{ref}}}=\sqrt{\frac{E_{\mathrm{k}}}{E_{\mathrm{k},\mathrm{ref}}}}$$

(12)

Firstly, the nonlinear static analysis is conducted to calculate the contact stiffness, which prepares for the subsequent modal analysis. Then, the modal parameters such as displacements, stress and energy are obtained by prestress modal analysis. Given the vibration stress σ, using equation and the parameters of (11), the tangential displacement A and the max vibration kinetic energy E_k can be obtained; Then selecting the corresponding contact model, the corresponding friction energy dissipation E_f can be achieved, and then using (6) to calculate the equivalent damping ratio under the current vibration stress ζ; By changing the vibration stress, the curve of damping ratio versus vibration stress can be acquired. The process of damping evaluation of shrouds under constant equivalent stiffness without considering the influence of shroud separation is shown as Figure 7.

Taking the friction coefficient μ = 0.3, the 1D and 2D contact models are used to calculate the damping curve for 3ND 1B and 5ND 1E modes, as seen in Figure 8.

Generally, under the same vibration stress, the damping ratio for 5ND 1E is larger than that for 3ND 1B, because the tangential displacement of the former is larger than that of the latter. It is worth noting the difference between 1D contact model and 2D contact model. The results for 3ND 1B mode show that the damping ratio of 2D contact model is larger than that of 1D, and the decreasing trend of damping ratio is not as obvious as that of 1D with the increase of vibration stress. The normal amplitude is several times larger than the tangential amplitude, and the separation will occur when the vibration stress is small. Although the separation will reduce the friction energy consumption, the max static friction also increases with the increase of vibration stress, resulting in more friction energy dissipation than 1D contact model. Figure 9 shows the contact hysteresis loop when the vibration stress is 10 MPa. The area enclosed by hysteresis loop (friction energy dissipation) of 2D contact model increases with the vibration stress approximately to the second power, while that of 1D contact model increases with the first power, which leads to that the change of equivalent damping ratio of 2D contact model with the vibration stress is not obvious.

However, for 5ND 1E mode, the normal displacement and the normal pressure of the contact surface change little, and the separation occurs only when the vibration stress is large enough. The friction energy dissipation between the two contact models is almost the same, so there is little difference in the damping characteristic curve, as Figure 10.

4. Damping-Evaluation with Altering Equivalent Stiffness

It is easy to separate for the slant shrouds of blade pair, which causes the change of the equivalent stiffness of the contact surface. At this time, the condition of constant equivalent stiffness in the Section 3 is too strong. In this section, taking the 3ND 1B mode as an example, considering the change of equivalent stiffness caused by the separation state, the damping characteristics of the blade shroud are analyzed.

4.1. Equivalent Stiffness on the Displacement Amplitude with the Separation

The shroud surfaces of the blade pair may contact or separate. The equivalent stiffness is assumed to be the contact stiffness when contacting, and be zero when separating during vibration, as shown in Figure 11.

4.1.1. Normal Force and Relative Motion

If the shrouds do not separate during vibration, the normal relative displacement is Δu_n(θ) = Bsin(θ), and the normal force is f_c,n(θ) = N₀ + k_nBsin(θ). When the amplitude increases, leading to the separation of the shrouds, the normal force is zero, as follows:

$$f_{\mathrm{c},\mathrm{n}}(\theta )=\{\begin{array}{ll}{N}_0+k_{\mathrm{n}}B\sin \theta ,& 0<\theta \le {\theta }_1\mathrm{or}{\theta }_2<\theta \le 2\pi \\ 0,& {\theta }_1<\theta \le {\theta }_2\end{array}$$

(13)

where θ₁ = π − arcsin(−N₀/(k_nB)) and θ₂ = 2π + arcsin(−N₀/(k_nB)) that are phase when separation, which can be calculated by f_c,n(θ) = 0. The Fourier expansion of the normal force only retains the constant and the first order term, which can be expressed as:

$${\widehat{f}}_{\mathrm{c},\mathrm{n}}(\theta )=F_{\mathrm{c},\mathrm{n}}^{\left(0\right)}+F_{\mathrm{c},\mathrm{n}}^{\left(1,\mathrm{c}\right)}\cos \theta +F_{\mathrm{c},\mathrm{n}}^{\left(1,\mathrm{s}\right)}\sin \theta$$

(14)

in which:

$$\begin{gathered} F_{\mathrm{c},\mathrm{n}}^{\left(0\right)}=\frac{1}{2\pi }\left[N_0\left(2\pi -{\theta }_2+{\theta }_1\right)+k_{\mathrm{n}}B\left(\cos {\theta }_2-\cos {\theta }_1\right)\right], F_{\mathrm{c},\mathrm{n}}^{\left(1,\mathrm{c}\right)}=0, \\ {F}_{\mathrm{c},\mathrm{n}}^{\left(1,\mathrm{s}\right)}=\frac{1}{\pi }\left[N_0\left(\cos {\theta }_2-\cos {\theta }_1\right)+k_{\mathrm{n}}B\left(\pi -\frac{{\theta }_2-{\theta }_1}{2}-\frac{\sin 2{\theta }_1-\sin 2{\theta }_2}{4}\right)\right] \end{gathered}$$

(15)

The normal force and its Fourier series are shown as Figure 12.

Assuming that the normal amplitude with contact is γ times of that without contact (γ > 1, defined as amplitude factor), the normal displacement during one cycle (Figure 13) is:

$$\Delta u_{\mathrm{n}}(\theta )=\{\begin{array}{ll}B\sin \theta ,& 0<\theta \le {\theta }_1 \mathrm{or} {\theta }_2<\theta \le 2\pi \\ \gamma B\sin \left(\alpha \theta +\beta \right),& {\theta }_1<\theta \le {\theta }_2\end{array}$$

(16)

where α = [π/2 + arcsin(sinθ₁/γ)]/(1.5π − θ₁) and β = 1.5π(1 − α).

The Fourier expansion of the normal relative displacement as Figure 14 can be:

$$\Delta {\widehat{u}}_{\mathrm{n}}(\theta )=\Delta U_{\mathrm{n}}^{\left(0\right)}+\Delta U_{\mathrm{n}}^{\left(1,\mathrm{c}\right)}\cos \theta +\Delta U_{\mathrm{n}}^{\left(1,\mathrm{s}\right)}\sin \theta$$

(17)

in which:

$$\begin{array}{ll}\Delta U_{\mathrm{n}}^{\left(0\right)}=& \frac{1}{2\pi }\left[\frac{\gamma B}{\alpha }\left(\cos \left(\alpha {\theta }_1+\beta \right)-\cos \left(\alpha {\theta }_2+\beta \right)\right)+B\left(\cos {\theta }_2-\cos {\theta }_1\right)\right], \Delta U_{\mathrm{c},\mathrm{n}}^{\left(1,\mathrm{c}\right)}=0\\ \Delta U_{\mathrm{c},\mathrm{n}}^{\left(1,\mathrm{s}\right)}=& \frac{1}{\pi }[B\left(\pi -\frac{{\theta }_2-{\theta }_1}{2}-\frac{\sin 2{\theta }_1-\sin 2{\theta }_2}{4}\right)\\ & +\frac{\gamma B}{2}(\frac{\sin \left(\alpha {\theta }_2-{\theta }_2+\beta \right)}{\alpha -1}-\frac{\sin \left(\alpha {\theta }_2+{\theta }_2+\beta \right)}{\alpha +1}-\frac{\sin \left(\alpha {\theta }_1-{\theta }_1+\beta \right)}{\alpha -1}+\frac{\sin \left(\alpha {\theta }_1+{\theta }_1+\beta \right)}{\alpha +1})]\end{array}$$

(18)

in which, as an unknown quantity, γ needs to be determined. It is noted that γ depends on to the amplitude B of normal relative displacement.

4.1.2. Amplitude Factor γ

The amplitude factor of the separation and contact can be calculated by the lumped parameter EoM of the system:

$$\{\begin{array}{l}M\ddot{u}+C_{\mathrm{c}}\dot{u}+K_{\mathrm{c}}u=f_{\mathrm{e}}, \mathrm{contact}\\ M\ddot{u}+C_{\mathrm{s}}\dot{u}+K_{\mathrm{s}}u=f_{\mathrm{e}}, \mathrm{separate}\end{array}$$

(19)

in which a circumferential harmonic excitation force is applied to the node at the position of 2/3 blade height and 1/2 chord length. The subscript c means shrouds in the contact state while s means in the separation state.

By the mass-normalized modes in two states for 3ND 1B φ_c and φ_s, and further normalize them to φ_c/φ_cb and φ_s/φ_sb, where φ_cb and φ_sb represent the amplitude of relative displacement from φ_c and φ_s. The EoM is transformed into the modal space by u = (φ_c/φ_cb)q and u = (φ_s/φ_sb)q, so as to obtain the EoM expressed by the modal coordinates:

$$\{\begin{array}{ll}{m}_{\mathrm{c}}\ddot{q}+c_{\mathrm{c}}\dot{q}+k_{\mathrm{c}}q=f_{\mathrm{e},\mathrm{c}}\cos \theta , & q<N_0/k_{\mathrm{n}}\\ m_{\mathrm{s}}\ddot{q}+k_{\mathrm{s}}q=f_{\mathrm{e},\mathrm{s}}\cos \theta ,& q\ge N_0/k_{\mathrm{n}}\end{array}$$

(20)

where m, c and k are modal mass, damping coefficient and stiffness, f_e means the modal excitation. m_c/m_s = (φ_sb/φ_cb)², k_c/k_s = m_c/m_s(ω_c/ω_s)² in which ω_c and ω_s are the natural angular frequencies for 3ND 1B modes in two states, respectively. f_e,c/f_e,s =φ_cfφ_sb/(φ_sfφ_cb) where φ_cf and φ_cb are the modal displacements at excited position. The EoM is dimensionless as:

$$\{\begin{array}{ll}\ddot{\overline{q}}(\overline{\omega }t)+2\zeta (\overline{B})\dot{\overline{q}}(\overline{\omega }t)+\overline{q}(\overline{\omega }t)={\overline{f}}_{\mathrm{e},\mathrm{c}}\cos (\overline{\omega }t), & \overline{q}<1\\ \frac{m_{\mathrm{s}}}{m_{\mathrm{c}}}\ddot{\overline{q}}(\overline{\omega }t)+\frac{k_{\mathrm{s}}}{k_{\mathrm{c}}}\overline{q}(\overline{\omega }t)={\overline{f}}_{\mathrm{e},\mathrm{s}}\cos (\overline{\omega }t), & \overline{q}\ge 1\end{array}$$

(21)

where $\overline{q}=q/B_0$, $\overline{\omega }=\omega /{\omega }_{\mathrm{c}}$, ${\overline{f}}_{\mathrm{e},\mathrm{c}}=f_{\mathrm{e},\mathrm{c}}/\left(m_{\mathrm{c}}{\omega }_{\mathrm{c}}^2{B}_0\right)$, ${\overline{f}}_{\mathrm{e},\mathrm{s}}=f_{\mathrm{e},\mathrm{s}}/\left(m_{\mathrm{c}}{\omega }_{\mathrm{c}}^2{B}_0\right)$, $B_0=N_0/k_{\mathrm{n}}$ is defined as the critical amplitude of separation, and $\overline{B}=B/B_0$. $\zeta \left(\overline{B}\right)$ is damping ratio calculated above with the constant equivalent stiffness (red dashed line in Figure 8). The Newmark method can be used to solve the dimensionless lumped-parameter EoM:

(1)

For each amplitude and frequency of excitation, the corresponding displacement-time curve can be calculated. It should be noted that since ζ is related to the unknown $\overline{B}$, an iterative solution between is needed;

(2)

By changing the excitation frequency, the max response amplitude and its corresponding frequency can be obtained, which are respectively resonance amplitude and resonance frequency; the amplitude of separation and contact in a period of resonance can be extracted, and the amplitude factor γ under the current resonance amplitude can be calculated;

(3)

Change the excitation amplitude to obtain amplitude factor. Figure 15 shows the curve of amplitude factor γ with the dimensionless normal amplitude $\overline{B}$.

4.1.3. Equivalent Stiffness on Amplitude

Given B, the corresponding γ extracted from Figure 15, Fourier coefficients of normal displacement from (14) and (17), the equivalent normal and tangential stiffness can be calculated by:

$$\{\begin{array}{l}{k}_{\mathrm{n}}^{\mathrm{eq}}=F_{\mathrm{c},\mathrm{n}}^{\left(1,\mathrm{s}\right)}/\Delta U_{\mathrm{n}}^{\left(1,\mathrm{s}\right)}\\ k_{\mathrm{t}}^{\mathrm{eq}}=k_{\mathrm{t}}{k}_{\mathrm{n}}^{\mathrm{eq}}/k_{\mathrm{n}}\end{array}$$

(22)

Define the dimensionless normal stiffness as ${\overline{k}}_{\mathrm{n}}^{\mathrm{eq}}=k_{\mathrm{n}}^{\mathrm{eq}}/k_{\mathrm{n}}$, and the curve of ${\overline{k}}_{\mathrm{n}}^{\mathrm{eq}}$ on $\overline{B}$ is shown as Figure 16. With the increase of $\overline{B}$, ${\overline{k}}_{\mathrm{n}}^{\mathrm{eq}}$ will decrease gradually until it is half of the original. This is because the ratio of the separation time to the whole vibration period is gradually approaching 0.5.

4.2. Parameters for Damping Evaluation

According to Figure 16, given the normal relative amplitude, the corresponding equivalent normal and tangential stiffness can be calculated. Then, the parameters required for damping can be achieved following the method in Section 3.2. The difference is that the equivalent stiffness in this section changes with the amplitude. For every normal amplitude, a modal analysis is needed, which takes a long time. This paper selects several modes of a single blade to fit the modes of bladed disk with different equivalent stiffness, so as to simplify the computation.

Nine modes of the single left blade ψ_j, j = 1,…, N_r (1B, 1E, first torsion (1T), 2B, 3B, 4B, rigid body rotating around z-axis, translating along x-axis and translating along y-axis) normalized by amplitude to form the reduction basis Ψ and the right blade is similar. The columns in Ψ are the selected modes which are linearly independent. The 3ND 1B mode φ($k_{\mathrm{n}}^{\mathrm{eq}}$) of the bladed disk is still selected as the target mode to be fitted. Noting that φ($k_{\mathrm{n}}^{\mathrm{eq}}$) varies with equivalent stiffness, several typical equivalent stiffness values are selected for modal analysis, and the corresponding 3ND 1B modes φ($k_{\mathrm{n}}^{\mathrm{eq}}$) are obtained. The 3ND 1B under three cases of 1, 0.5, 0 times stiffness calculated in Section 2.3.

Table 4. Fitting coefficient of reduced modes of left blade.

Single Blade Reduction Basis ψj	χj of Different Stiffness
	φ(0)	φ(0.5)	φ(1)
1B	−0.56	−0.54	−0.53
1E	0.07	0.07	0.07
1T	0.03	0.04	0.04
2B	−0.02	−0.03	−0.04
3B	0.00	0.00	0.01
4B	0.00	−0.01	−0.01
RB Rotating z	−1.15	−1.25	−1.31
RB Translating x	0.69	0.75	0.78
RB Translating y	0.01	0.01	0.02
error ε	0.00537%	0.00585%	0.00656%

and

Table 5. Fitting coefficient of reduced modes of right blade.

Single Blade Reduction Basis ψj	χj of Different Stiffness
	φ(0)	φ(0.5)	φ(1)
1B	0.55	0.54	0.53
1E	0.07	0.08	0.08
1T	0.03	0.04	0.05
2B	0.02	0.03	0.04
3B	0.00	0.00	−0.01
4B	0.00	0.01	0.01
RB Rotating z	−1.18	−1.26	−1.32
RB Translating x	0.70	0.75	0.79
RB Translating y	0.05	0.05	0.05
error ε	0.00507%	0.00557%	0.00616%

respectively list the coefficients χ_j and errors ε of the selected reduced basis vectors of the left and right blades when fitting the 3ND 1B mode of the bladed disk under the above typical equivalent stiffness values.

where χ_j and ε [29] can be expressed as:

$${\chi }_j=\frac{{\phi }^T(k_{\mathrm{n}}^{\mathrm{eq}}){\psi }_j}{\Vert \phi (k_{\mathrm{n}}^{\mathrm{eq}})\Vert \Vert {\psi }_j\Vert }$$

(23)

The approximate mode after fitting is:

$$\tilde{\phi }={\displaystyle \sum _j^{N_r}{\chi }_j{\psi }_j}$$

(24)

and $\epsilon =\Vert \phi -\tilde{\phi }\Vert /\Vert \phi \Vert$.

For other equivalent stiffness values $k_{\mathrm{n}}^{\mathrm{eq}}$, the fitting coefficients ${\chi }_j^{k_{\mathrm{n}}^{\mathrm{eq}}}$ of each reduced mode are obtained by linear interpolation of data in Table 4 and Table 5 as:

$$\frac{{\chi }_j^{k_{\mathrm{n}}^{\mathrm{eq}}}-{\chi }_j^{0.5}}{k_{\mathrm{n}}^{\mathrm{eq}}/k_{\mathrm{n}}^{}-0.5}=\frac{{\chi }_j^{1.0}-{\chi }_j^{0.5}}{1.0-0.5}$$

(25)

Thus, the calculation steps of modal parameters for equivalent damping can be summarized as:

(1)

Extract the following modal data from each mode ψ_j in reduction basis Ψ: the displacements in three directions ${\left(X_{\mathrm{c},\mathrm{x},j}^{\mathrm{ro}},X_{\mathrm{c},\mathrm{y},j}^{\mathrm{ro}},X_{\mathrm{c},\mathrm{z},j}^{\mathrm{ro}}\right)}^{\mathrm{T}}$, the vibration stress at critical location ${\sigma }_j^{\mathrm{ro}}$, the max kinetic energy $E_{\mathrm{k},j}^{\mathrm{ro}}$;

(2)

Given the amplitude B, the new equivalent stiffness $k_{\mathrm{n}}^{\mathrm{eq}}=F_{\mathrm{c},\mathrm{n}}^{\left(1,\mathrm{s}\right)}/\Delta U_{\mathrm{n}}^{\left(1,\mathrm{s}\right)}$ and the fitting coefficients of the reduced modes under the current stiffness ${\chi }_j^{k_{\mathrm{n}}^{\mathrm{eq}}}$ are obtained;

(3)

Combined with Section 3.2, the displacements in three directions under this equivalent stiffness can be obtained by $X_{\mathrm{c},\mathrm{x}}^{k_{\mathrm{n}}^{\mathrm{eq}}}={\displaystyle \sum _{j=1}^{N_{\mathrm{r}}}{\chi }_j^{k_{\mathrm{n}}^{\mathrm{eq}}}{X}_{\mathrm{c},\mathrm{x},j}^{\mathrm{ro}}}$, for $X_{\mathrm{c},\mathrm{y}}^{k_{\mathrm{n}}^{\mathrm{eq}}}$ and $X_{\mathrm{c},\mathrm{z}}^{k_{\mathrm{n}}^{\mathrm{eq}}}$ in the same way; The relative amplitudes A_ref and B_ref in the local coordinate system are calculated; Vibration stress vector under this stiffness is ${\sigma }_{}^{k_{\mathrm{n}}^{\mathrm{eq}}}={\displaystyle \sum _{j=1}^{N_{\mathrm{r}}}{\chi }_j^{k_{\mathrm{n}}^{\mathrm{eq}}}{\sigma }_j^{\mathrm{ro}}}$, in which the max component is taken as σ_ref; The max vibration kinetic energy is ${\sigma }_{\mathrm{k},\mathrm{ref}}^{}={\displaystyle \sum _{j=1}^{N_{\mathrm{r}}}{\chi }_j^{k_{\mathrm{n}}^{\mathrm{eq}}}{E}_{\mathrm{k},j}^{\mathrm{ro}}}$.

4.3. Damping Characteristic Curve with the Altering Stiffness

Given normal amplitude B, the real and modal parameters: normal and tangential displacement, stress and kinetic energy satisfy:

$$\frac{\sigma }{{\sigma }_{\mathrm{ref}}}=\frac{B}{B_{\mathrm{ref}}}=\frac{A}{A_{\mathrm{ref}}}=\frac{\sqrt{E_k}}{\sqrt{E_{k,\mathrm{ref}}}}$$

(26)

Combined with the modal parameters calculated in Section 4.2, the corresponding friction energy consumption can be obtained by selecting the corresponding contact model, and then the equivalent damping ratio under different vibration stress can be calculated by Equation (6), as shown in Figure 17, where the upper right part is to calculate the change of equivalent stiffness with the amplitude of normal motion, and the lower right part is the step to fit the vibration mode. These two parts provide the necessary parameters for the calculation of equivalent damping with altering stiffness, which need to be conducted in advance. Given the vibration stress B, using equation and the parameters of (26), the tangential displacement A and the max vibration kinetic energy E_k can be obtained; Then using the 2D contact model, the corresponding friction energy dissipation E_f can be achieved, and then using (6) to calculate the equivalent damping ratio under the current vibration stress ζ.

For the shrouded bladed disk in this paper, Figure 18 shows the damping characteristic curves obtained by three methods. The damping ratio of 2D contact model with altering stiffness method is lower than that with constant stiffness method, but it is still higher than that of 1D contact model. It is worth noting that the calculation time of a damping characteristic curve is very short, which is very suitable for engineering applications. The calculation time of three methods that are 1D and 2D contact model with constant stiffness and 2D with altering stiffness for 3ND 1B under the conditions of a normal PC configuration (Intel^® Core™ i7−6700 CPU @ 3.40 GHz, RAM 16.0 GB) is 0.21, 0.08 and 1.15, respectively.

5. Discussion

5.1. Vibration-Reduction Design

The mechanism of dry friction damping is that the friction generated by the relative motion between the contact surfaces dissipates the vibration kinetic energy. From the perspective of energy, the ratio of friction dissipation to the max vibration kinetic energy of the system can be used to evaluate the damping effect. In this paper, the damping characteristic curve is used to study the vibration reduction, which is the curve of the damping ratio provided by the system with the vibration stress, as shown in Figure 19. The important parameters are: the critical vibration stress σ_cr which is the resonance vibration stress at critical location, the peak damping ratio ζ_max and the corresponding vibration stress σ_m. During the design phase, the objective is to increase the peak damping ratio as much as possible under the allowable vibration stress; at the same time, the critical vibration stress should be relatively small; in addition, the typical vibration stress of the blade should preferably fall in the left half of the damping characteristic curve, so that when the vibration stress increases, the damping ratio provided by the shroud will also increase.

The flow chart of vibration reduction design based on the damping characteristics of the blade shroud is shown in Figure 20, where σ_a is the allowable vibration stress of the system, which is often located in the blade body.

5.2. Sensitivity Analysis

The parameter sensitivity analysis of the damping curve for the 3ND 1B of the bladed disk is carried out that are phase angle between normal and tangential motion, the contact stiffness, friction coefficient and initial normal force.

5.2.1. Phase Angle between Normal and Tangential Motion

When the 2D contact model is used, the phase angle will greatly affect the hysteresis loop, and then the energy dissipation, thus influencing the damping characteristic curve. Figure 21 shows the damping characteristic curves when the phase angle is from 0° to 90°. The phase angle of black solid line is 90° and that of red solid line is 0° in the two figures.

5.2.2. Contact Stiffness

In the contact model, the contact stiffness is an important parameter, but difficult to determine. Its value directly affects the friction energy dissipation and vibration response. It is related to the contact surface pressure, geometry, friction coefficient, temperature and surface roughness. At present, there is no reliable and unified numerical calculation method [30], and it is usually obtained by experimental measurement as input parameters for simulating. Therefore, its value often come with variable degrees of uncertainty. Figure 22 shows the sensitivity of the damping curve to variations in the contact stiffness. 2D contact model, phase angle and normal force in Table 1 are used herein. The stiffness values in this figure are given in the form of (k_t, k_n) in which k₀ = 2.06 × 10⁶ N/m. The black color indicates k_t = k_n: when the contact stiffness decreases from 4k₀ to k₀, ζ_max decreases from 3.4% to 0.5%, and σ_cr increases from 0.93 MPa to 7.13 MPa; With k_t = k_n = 4k₀, when σ > σ_m, the equivalent damping ratio decreases obviously with the increase of vibration stress. The multicolor indicates k_t ≠ k_n: when the contact stiffness varies from k_t = k_n/2 to k_t = k_n/4, ζ_max decreases from 1.45% to 0.83%, and σ_cr increases from 2.61 MPa to 3.07 MPa.

5.2.3. Initial Contact Normal Force

The initial normal force is a design parameter, which is related to the interference between the shrouds and the rotating speed. The amplitude first decreases and then increases with the increase of normal force, so there is an optimal value to minimize the amplitude [30,31]. Figure 23 shows the effect of the normal force on the damping curves. 2D contact model, phase angle and contact stiffness in Table 1 are used herein. When the initial normal force increases, the sliding difficulty of sliding increases, and the critical vibration stress increases. However, the peak damping ratio remains basically unchanged, but its corresponding vibration stress increases.

5.2.4. Friction Coefficient

With the increase of blade working time, the contact surface of shroud will become smooth gradually, resulting in the change of friction coefficient. Take the friction coefficient μ = 0.1, 0.2, 0.3, 0.4, 0.5 to study its influence on the damping curves as in Figure 24. Initial normal force and contact stiffness in Table 1 are used herein. The critical vibration stress increases with the friction coefficient due to the increase of the Coulomb limit, which makes the contact surface more difficult to slip. The peak damping ratio also grows with the friction coefficient, which is due to the increase of friction force and energy dissipation. Both the critical vibration stress and peak damping ratio first increase and then stabilize with the friction coefficient as in Figure 25.

6. Conclusions

In this paper, a method for analyzing the damping characteristics of shrouded turbine blades is proposed, which takes into account the variation of the contact equivalent stiffness of the shroud. A 2D contact model with variable normal force is adopted, and the influence of separation state on equivalent stiffness is considered. The curve of equivalent stiffness with amplitude is calculated. Then, the damping characteristic curve of shrouds, that is, the relationship of the equivalent damping ratio with vibration stress is calculated by the energy method. The design objective of this method is to provide a larger damping ratio under the allowable vibration stress. Based on ordinary modal analysis, this method is simple and time-saving, which is effectively useful in the design stage.

The 3ND 1B mode of bladed disk is taken as an example. When considering the influence of shroud separation on equivalent stiffness, with the increase of blade amplitude, the equivalent stiffness will gradually decrease to half of the original. Under the same vibration stress, the damping ratio of 2D contact model with variable equivalent stiffness method is lower than that of 2D contact model with constant stiffness method because there is no friction energy consumption in the separation stage, but it is still higher than that of 1D contact model.

In the 2D contact model, the phase angle between normal and tangential motion will significantly affect the damping characteristic curve, and there is no obvious law. With the rising of contact stiffness, the critical vibration stress descends and the peak damping ratio ascends. Under the allowable vibration stress of the blade, the initial contact normal force has an optimal value, which makes the damping ratio of the system maximum. The critical vibration stress and peak damping ratio first go up and then tend to be stable with the growing of friction coefficient.

Despite the fact the variation of the equivalent stiffness of the blade shroud is considered in this method, the linear mode of the bladed disk system is used directly, which leads to the expression of vibration displacement being only a single harmonic. In addition, in the actual vibration of the blade, if the amplitude is large, there will be slip, stick, and separation in the contact state during one cycle, which appear alternately in turn. In the actual vibration process, it is difficult to determine the mode of the structure. As mentioned earlier, some researchers [19,20] put forward the concept of “nonlinear mode” to study the dynamic system with local nonlinearities. This paper uses the mode with altering stiffness to make some studies, which still make some assumptions. And the corresponding modes can be considered as one combination of the mode of the system with normal stiffness, that with three direction stiffness and free mode. This method does not consider the modes separately, but always uses the modes with stiffness in three directions, and through Fourier expansion, the influence of the three states is reflected as the variation of equivalent stiffness. However, the vibration form of the shrouded bladed disk still needs further study.