Vibration analysis of mistuned damped blades with nonlinear friction and contact

Abstract

In general, only in tangential direction, friction motion between blade dampers is concerned for vibration analysis of mistuned bladed disk. However, in practical, excitation acting on blades inevitably causes normal movement at friction interface due to the existence of angle between excitation force and contact surface. This fact leads to a variation of normal force or even a result of contact separation, which determine the maximum static friction for tangential frictional vibration. In order to assess the realistic nonlinear forced vibration of mistuned disk assemblies subjected to actual excitation with nonlinear friction and contact, an efficient method is developed by using large-scale finite element model and two-dimensional friction model. In the proposed paper, attentions are paid on the realistic coupled vibration and the impact of variable normal vibration on forced response to a traveling wave excitation. Three types of mistuned parameters, including tangential stiffness of friction surface, initial shroud gap (corresponding to preload normal force), and structural stiffness of blades are taken into account for analyzing the effect of mistuning on the coupling vibration, especially for cases of contact separation. Furthermore, by means of curve fitting of amplitude variation of normal forces, the dependence of friction motion of bladed disk studied at circumferential resonance on variation of axial displacement, even with contact separation involved, is addressed. The results show a complex nonlinear vibration with friction and contact, especially a complex multiple unstable vibration observed with repeated nonlinear snap at circumferential resonance for the cases of contact separation. With regard to a small initial normal force, a large amplitude value of force variation can greatly reduce the maximum response amplitude for mistuned blades, which is worth to pay attentions for the vibration analysis of bladed disk with nonlinear friction and contact.

Introduction

Blades are the key component of turbine machinery and responsible for the energy transfer from heat energy to mechanical kinetic energy in operating conditions that including centrifugal force, gas excitation force, and others. Friction damping designed on turbine blades effectively dissipates energy and reduces vibration for bladed disk assembles in operating. The dynamical characteristics of turbine bladed disk with dampers have been well studied through a large amount of numerical and experimental investigations though the behavior of dry friction is nonlinear.^1–7

Blade mistuning is one of the problems that hardly averted in industrial production processes of realistic turbine blades, such as blade manufacturing, installation tolerances, and wear, etc.^8–10 Indeed, the vibration amplitude and stress level of blades are highly sensitive to mistuning, especially in frequency steering regions. One evident effect of mistuning is scatters of natural frequencies, mode shapes, and properties in contact interfaces for dampers. The study of mistuned bladed disk has important practical significance for security and reliability. ¹¹ In order to accelerate solving process in practice, the vibration of damped blades is usually assumed to be a linear system in earlier studies.^12,13 The achievements including analytical methods and parameter impacts have been applied in the industrial field to enhance unit performance. However, the contact interaction force is complexly nonlinear for actual bladed disk assemblies in operation. The huge errors, induced between the actual vibration and results predicted simply by using linear models, are difficult to tolerate, especially for some precise requirements in blade design nowadays.^14,15

In order to reveal the real nonlinear dynamic of mistuned bladed disk system with dampers, researchers perform numerous analyses of mistuned forced vibration with nonlinear friction extensively. Sextro et al. ¹⁶ developed a method to calculate statistically envelopes of frequency response functions for a nonlinear mistuned damped system based on Weibull distribution of the vibration amplitudes. Later, Panning et al. ¹⁷ pointed out the advantages and disadvantages of both underplatform dampers and exhibited the influence of mistuning effects of mistuned bladed disk in terms of statistical natural frequencies. Kumar and Narayanan ¹⁸ presented extensions to a finite difference technique for solution of multidimensional Fokker–Planck equation of mistuned bladed disk assembly subjected to white noise excitation and investigated the effects of stiffness and damping mistuning on the forced response. Mitra et al. ⁹ analyzed random mistuning patterns of contact stiffness on nonlinear contact responses and observed the variation of amplification factors. Joannin et al. ¹⁹ explained the nonlinear complex modal properties that vary with the vibration amplitude in given cases. The effect of mistuning magnitude on the response of bladed disk assemblies was evaluated.

These studies assumed that the normal force on friction interface is invariable.^20,21 However, as bladed disk is subjected to complex variable external loads, especially excitation acting hardly along the tangential direction of contact interface, normal forced vibration is inevitably aroused. Actually, blades on service are discrepant across full circle and not always keeping locked, which even results in the problem of contact separation. The displacement in the normal direction is a real issue, even with a strong nonlinearity due to contact separation. Therefore, forced response of mistuned blades obtained by conventional nonlinear solution methods may be not accurate, because of failing to account for a complex nonlinear effect caused by the variable normal displacement, such as snap-through and snapback, etc. Petrov and Ewins ²² proposed a method for solving friction function and tangential Hessian matrix of two-dimensional friction model and investigated forced response of contact point subjected to tangential and normal coupling excitation. This research provides modeling contact interaction and fast calculation of interface characteristics for cases with partial separation of contact surfaces. Based on this friction model, vibration analysis of large-scale finite element model of practical mistuned blades with dampers with nonlinear friction and contact could be performed.

Differ from the works Petrov has done, another effective method for analyzing forced response of mistuned blades is developed in this paper based on modal substructure reduction method on the whole structure and not on a single sector as in Petrov and Ewins. ²³ The advantage of this method is that various mistuned factors in the system can be introduced comprehensively with accurate interactions obtained between sectors. Similarly, two-dimensional contact interface elements are built between shrouds to establish the high fidelity finite element model of the mistuned damped blades. But the coupling effect of multidirectional vibration on whole response of blades in all directions is further explored in this paper. Firstly, the forced vibration of the tuned bladed disk subjected to periodic excitation in axial and circumferential directions is analyzed and a complex unstable frequency regime along with normal contact separation is founded. Secondly, different mistuned parameters are introduced to observe the effect of mistuning on the vibration in each direction. The response curves of mistuned blades are scattered, and the amplitude amplification factors of bladed disk system mistuned by different parameters are calculated under variable mistuning level. In addition, the normal forces on interfaces are fitted by curves to simulate a case of coincident resonances for axial and circumferential vibrations and to explore the impact of variable normal force on the tangential resonance.

Method for nonlinear vibration

Multi-harmonic balance equation

For bladed disk assemblies including nonlinear friction interfaces, the equation of motion generally is expressed in the following differential form

$$\mathbf{M}\ddot{\mathbf{q}}\left(t\right)+\mathbf{C}\dot{\mathbf{q}}\left(t\right)+\mathbf{Kq}\left(t\right)={\mathbf{F}}_{\text{N}}\left(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}}\right)+{\mathbf{F}}_{\text{L}}(t)$$

(1)where q(t) is the displacement vector for all degrees of freedom (DOF) of analyzed structure; M, C, and K are the mass, viscous damping, and stiffness matrices of all DOFs; F_L is the vector of external excitation force acting on bladed structure; F_N is the vector of nonlinear interface force with respect to displacement, velocity, or acceleration. It should be noted that stiffness matrix K comprise elastic stiffness matrix and two other geometric stiffness matrices respectively reflecting effects of centrifugal forces and spin-softening. The excitation force F_L is assumed usually as a sinusoidal waveform function and satisfy ${\mathbf{F}}_{\text{L}}\left(t\right)={\mathbf{F}}_{\text{L}}\left(t+T\right)$ , in which T denotes vibration periodicity. As vibration response is periodic, the steady-state solution of equation (1) in time domain can be expressed as a restricted Fourier series. The exact solution can be approximated accurately by selecting a finite specific harmonic term, given by the following equation

$$\mathbf{q}\left(t\right)=Q_0+{\displaystyle \sum _{i=1}^{\mathit{Nh}}{Q}_i^c\cos \hspace{0.17em}\left(i\omega t\right)+Q_i^s\cos \hspace{0.17em}\left(i\omega t\right)}$$

(2)where i is the number of harmonic order, Nh is number of harmonic, Q is harmonic coefficient of displacement, and ω is principal angular frequency corresponding to the vibration periodicity T in the form: $\omega =2\pi /T$ . Meanwhile, the motion vector in frequency domain can be expressed by the following formula

$$\mathbf{q}(t)={\mathbf{H}}^{\text{T}}(t)\mathbf{Q}=\left\{1,\cos \hspace{0.17em}\omega t,\sin \omega t,\dots ,\cos \hspace{0.17em}\mathit{Nh}\omega t,\sin \hspace{0.17em}\mathit{Nh}\omega t\right\}\mathbf{Q}$$

(3)

$$\mathbf{Q}={\left\{Q_0,Q_1^c,Q_1^s,Q_2^c,\dots ,Q_{\mathit{Nh}}^c,Q_{\mathit{Nh}}^s\right\}}^{\text{T}}$$

(4)where H is harmonic function vector for translating terms from frequency domain to time domain. The formulation of motion in frequency domain can be determined with a substituted of the expansion in equation (3) into equation (1). The harmonic components can be written as

$$\mathbf{A}(\omega )\mathbf{Q}={\mathbf{F}}_{\text{N}}(\mathbf{Q})+{\mathbf{F}}_L$$

(4)

$$\mathbf{A}=\overline{\mathbf{M}}{\mathbf{T}}_1+\overline{\mathbf{C}}{\mathbf{T}}_2+\overline{\mathbf{K}}$$

(5)where A(ω), F_N(Q), and F_L are accordingly linear dynamic stiffness matrix and vectors of harmonic components of nonlinear force and external exciting force, respectively. Dynamic stiffness matrix A(ω) is dependent on mass, viscous damping, and stiffness matrices. The items in equation (5) are written as

$$\overline{\mathbf{M}}=\hspace{0.17em}\text{diag}(\mathbf{M},\dots ,\mathbf{M})$$

(6)

$$\overline{\mathbf{C}}=\hspace{0.17em}\text{diag}(\mathbf{C},\dots ,\mathbf{C})$$

(7)

$$\overline{\mathbf{K}}=\hspace{0.17em}\text{diag}(\mathbf{K},\dots ,\mathbf{K})$$

(8)

$${\mathbf{T}}_1=\text{diag}\left(0,-\omega \mathbf{I},-\omega \mathbf{I},-{\left(2\omega \right)}^2\mathbf{I},-{\left(2\omega \right)}^2\mathbf{I},\dots ,-{\left(\mathit{Nh}\omega \right)}^2\mathbf{I},-{\left(\mathit{Nh}\omega \right)}^2\mathbf{I}\right)$$

(9)

$${\mathbf{T}}_2=\text{diag}\left(0,\omega \left[\begin{array}{cc}0& \mathbf{I}\\ -\mathbf{I}& 0\end{array}\right],\dots ,\mathit{Nh}\omega \left[\begin{array}{cc}0& \mathbf{I}\\ -\mathbf{I}& 0\end{array}\right]\right)$$

(10)

Reduced order method

The 3D finite element model of bladed disk system possesses a great amount of DOFs, which significantly increases the demand for computational hardware. In order to improving the computational efficiency, fixed interface modal synthesis method is applied here to reduce the total DOFs. By neglecting damping matrix $\overline{\mathbf{C}}$ in equation (5) and nonlinear force vector F_N in equation (4), the time-invariable linear equation of disk structure is determined in the following form

$$\overline{\mathbf{M}}{\mathbf{T}}_1\mathbf{Q}+\overline{\mathbf{K}}\mathbf{Q}={\mathbf{F}}_{\text{L}}$$

(11)

The full-bladed disk can be customary dived into n + 1 segments comprising n segments of blades and 1 segment of disk, as in Lim et al. ²⁴ The DOFs of each segment are divided into slave displacement vector Q_s and main displacement vector Q_m, which mainly includes the DOFs of shroud contacts and root contacts, etc. During the system freely vibration with none external exciting effect, the nodal forces are 0 except for interface. The stiffness matrix in equation (8) and the vibration equation for ith segment can be written as respectively

$$\mathbf{K}=\text{diag}\left({\mathbf{K}}_1,\cdots ,{\mathbf{K}}_i,\cdots ,{\mathbf{K}}_n,{\mathbf{K}}_{\text{d}}\right)$$

(12)

$$\left(\left[\begin{array}{cc}{\mathbf{M}}_{i,\text{ss}}& {\mathbf{M}}_{i,\text{sm}}\\ {\mathbf{M}}_{i,\text{ms}}& {\mathbf{M}}_{i,\text{mm}}\end{array}\right]{\mathbf{T}}_{i,1}+\left[\begin{array}{cc}{\mathbf{K}}_{i,\text{ss}}& {\mathbf{K}}_{i,\text{sm}}\\ {\mathbf{K}}_{i,\text{ms}}& {\mathbf{K}}_{i,\text{mm}}\end{array}\right]\right)\left\{\begin{array}{c}{\mathbf{Q}}_{i\text{,s}}\\ {\mathbf{Q}}_{i,\text{m}}\end{array}\right\}=\left\{\begin{array}{c}0\\ {\mathbf{F}}_{i,\text{L}}\end{array}\right\}$$

(13)

Therefore, the reduced equation of system DOFs can be expressed as a combination of Ritz vector component Φ_rom and a reduced modal basis u _{i, k} in the following form

$$\left\{\begin{array}{c}{\mathbf{Q}}_{i,\text{s}}\\ {\mathbf{Q}}_{i,\text{m}}\end{array}\right\}=\left[\begin{array}{cc}{\mathbf{\Phi }}_{i,k}& {\mathbf{\Phi }}_{i,\text{sm}}\\ 0& \mathbf{I}\end{array}\right]\left\{\begin{array}{c}{\mathbf{u}}_{i,k}\\ {\mathbf{Q}}_{i,\text{m}}\end{array}\right\}={\mathbf{\Phi }}_{\text{rom}}{\mathbf{u}}_{i,k}$$

(14)where the main modal matrix Φ _{i , k} denotes the modes of first k orders retained by the main matrix and is obtained by solving the following eigenvalue equation

$$\left({\mathbf{M}}_{i,\text{ss}}{\mathbf{T}}_{i,\text{1}}+{\mathbf{K}}_{i,\text{ss}}\right){\mathbf{Q}}_{i\text{,s}}=0$$

(15)

It can be seen that equation (15) is determined by fixing the substructure master displacement Q _{i ,s}, and that is why the reduction method is called the fixed interface modal synthesis method. In terms of the technique provided by Bladh et al., ²⁵ the constrained modal matrix Φ _{i ,sm} can be written as

$${\mathbf{\Phi }}_{i\text{,sm}}=-{\left({\mathbf{K}}_{i\text{,ss}}\right)}^{-1}{\mathbf{K}}_{i,\text{sm}}$$

(16)

The reduced dynamic equation from equation (4) can be finally rewritten in the form

$$\overline{\mathbf{A}}\mathbf{Q}=\hspace{0.17em}\left(\overline{\mathbf{M}}{\mathbf{T}}_1+\overline{\mathbf{C}}{\mathbf{T}}_2+\overline{\mathbf{K}}\right)\mathbf{Q}={\overline{\mathbf{F}}}_{\text{N}}\left(\mathbf{Q}\right)+{\overline{\mathbf{F}}}_{\text{L}}$$

(17)

$$\tilde{\vee }={\mathbf{\Phi }}_{\text{rom}}^{\text{T}}\vee {\mathbf{\Phi }}_{\text{rom}},\vee =\overline{\mathbf{M}},\overline{\mathbf{C}},\overline{\mathbf{K}}$$

(18)

$$\tilde{\wedge }={\mathbf{\Phi }}_{\text{rom}}^{\text{T}}\wedge ,\wedge ={\mathbf{F}}_{\text{N}}\left(\mathbf{Q}\right),{\mathbf{F}}_{\text{L}}$$

(19)

Arclength method and Newton–Raphson method

Generally, due to some calculating limitations at turning points of response curve, it is difficult to deal with the convergence for the curvature of vibration response. To overcome these limitations, a continuation parameter, ds, along the arclength is chosen to trace the loci of the solutions. ²⁶ The amplitude solution can be obtained by adding a special equation into motion equation of system. The procedure can be regarded as searching solution on the hypersphere space centered on the final converged solution and radius ds. The displacement vector Q can be normalized by using arclength ds in the form

$${\left(\frac{\Vert d\mathbf{Q}\Vert }{\mathit{ds}}\right)}^2+{\left(\frac{d\omega }{\mathit{ds}}\right)}^2=1$$

(20)where variable term $d\mathbf{Q}=\mathbf{Q}-{\mathbf{Q}}_0$ , $d\omega =\omega -{\omega }_0$ . ${\mathbf{Q}}_0$ , and ${\omega }_0$ are the last converged solution and last input frequency respectively. Furthermore, a new variable $\mathbf{Z}={\{\omega ,{\mathbf{Q}}^{\text{T}}\}}^{\text{T}}$ is derived to contribute the solution of vibration with arclength. To help ensuring that the solution Z converges within a restricted step, appropriate initialization need to be applied for response parameters. The initial prediction ${\mathbf{Z}}_j^{\text{P}}$ at jth frequency is generally depended on the last two results and expressed in the following form

$${\mathbf{Z}}_j^{\text{P}}={\mathbf{Z}}_{j-1}-\mathit{ds}\frac{{\mathbf{Z}}_{j-1}-{\mathbf{Z}}_{j-2}}{\Vert {\mathbf{Z}}_{j-1}-{\mathbf{Z}}_{j-2}\Vert }$$

(21)

In addition, special processing is performed in calculating prediction Z in frequency step j = 1, which generally originates from the steady-state response of linear system undergoing excitation. It should be noted that, for a solution, convergence rate is very dependent on the choice of the initial value. The arc length in curvature is variable, with the controlling formula as follows

$$d{s}_j=d{s}_{j-1}\frac{\mathit{ite}{r}_{\text{p}}}{\min \hspace{0.17em}\left\{\mathit{ite}{r}_{j-1},10\right\}}$$

(22)where iter denotes iterative step and subscript p denotes expectancy of iterative step. Generally, iter_p is set as 4. However, this iterative calculation method performs not well when contact friction turns to be strongly nonlinear, such as slip. So, additional judgments are required to improve the convergence speed of calculating response curve. It indicates the trend of curve changes drastically when initial predictions deviate from convergence results, and arc length would be reduced to be reasonable, for instance, about half, to obtain a better refined curve.

Newton–Raphson method is one of most efficient methods for solution of nonlinear equation. The displacement harmonic vector Z and residual vector R are updated continuously by means of Fast Fourier Transform (FFT) and its inverse Fast Fourier Transform (FFT⁻¹) at each iterations of Newton–Raphson solver, according to the following steps

$${\mathbf{Z}}^{k-1},\hspace{1em}{\mathbf{R}}^{k-1}\stackrel{{\text{FFT}}^{\text{-1}}}{\to }{q}^{k-1},\hspace{1em}{f}_{\text{N}}^{k-1}\stackrel{\text{FFT}}{\to }{\mathbf{F}}_{\text{N}}^{k-1},\hspace{1em}{\mathbf{J}}^{k-1}\overrightarrow{\hspace{0.17em}}{\mathbf{Z}}^k,\hspace{1em}{\mathbf{R}}^k$$

where superscript k is the number of current iteration at jth frequency. The Jacobi matrix $\mathbf{J}$ is a matrix of derivatives of the Fourier coefficients for the residuals with respect to the Fourier coefficients for relative displacements. The iterative process of Newton-Raphson solver is expressed in the form

$${\mathbf{Z}}_j^k={\mathbf{Z}}_j^{k-1}-\frac{\mathbf{R}\left({\mathbf{Z}}_j^{k-1}\right)}{{\mathbf{J}}_j^{k-1}}$$

(23)

$$\mathbf{R}\left({\mathbf{Z}}_j^{k-1}\right)=\left[{\mathbf{R}}_{\omega }\left({\omega }_j^{k-1}\right);{\mathbf{R}}_Q\left({\mathbf{Z}}_j^{k-1}\right)\right]$$

(24)

$${\mathbf{R}}_{\omega }={\Vert {\mathbf{Z}}_j^{k-1}-{\mathbf{Z}}_{j-1}\Vert }^2-d{s}_{j}2$$

(25)

$${\mathbf{R}}_Q=\left(\tilde{\mathbf{M}}{\mathbf{T}}_1+\tilde{\mathbf{C}}{\mathbf{T}}_2+\tilde{\mathbf{K}}\right){\mathbf{Z}}_j^{k-1}-{\tilde{\mathbf{F}}}_{\text{N}}\left({\mathbf{Z}}_j^{k-1}\right)-{\tilde{\mathbf{F}}}_{\text{L}}$$

(26)

$$J_j^{k-1}=\frac{\partial \mathbf{R}\left({\mathbf{Z}}_j^{k-1}\right)}{\partial \mathbf{Q}}=\left[\begin{array}{cc}2\left({\omega }_j^{k-1}-{\text{ω}}_{j-1}\right)& 2\left({\mathbf{Q}}_j^{k-1}-{\mathbf{Q}}_{j-1}\right)\\ \frac{\partial \tilde{\mathbf{A}}}{\partial \omega }{\mathbf{Q}}_j^{k-1}-\frac{\partial {\mathbf{F}}_{\text{N}}\left({\mathbf{Z}}_j^{k-1}\right)}{\partial \omega }& \tilde{\mathbf{A}}-\frac{\partial {\mathbf{F}}_{\text{N}}\left({\mathbf{Z}}_j^{k-1}\right)}{\partial \mathbf{Q}}\end{array}\right]$$

(27)

Friction interface elements

Nonlinear friction interface elements are applied to build contact relations of relative motion between blade contacts, as shown in Figure 1. Springs with stiffness k_x and k_y are arranged in tangential, x, and normal, y, directions respectively between surface 1 and contact point A, and the motion between point A and surface 2 follows the Coulomb’s Friction Law. The tangential relative displacement of contact point A to surface 2 is w. The properties of contact interface are also described by a friction coefficient μ and an initial normal preload N₀ (a negative value is possible, corresponding to an initial gap g₀=－N₀/k_y). When surfaces remain contact, the relative frictional force between surfaces is f =－k_x(x-w) during motion cycle, otherwise f = 0 when surfaces separate.

Figure 1.

Nonlinear friction interface element.

The tangential and normal relative displacements of surface 1 to surface 2 are x and y. The displacement y in normal direction determines whether the contact interface is in contact. Two contact states, slip and stick, can be observed for surface 1 in contact, and the displacement w of contact point A is limited by the maximum static frictional force μN. The initial frictional force and initial displacement are $f_x^{\mathrm{\hspace{0.17em}}0}$ and x₀ when stick occurs. Expressions for nonlinear interaction forces are written in the following form

$$f_x=\{\begin{array}{cc}{f}_x^{\mathrm{\hspace{0.17em}0}}+k_x\left(x-x_0\right)& \text{stick}\\ \delta \mu f_y& \text{slip}\\ 0& \text{separation}\end{array},f_y=\{\begin{array}{cc}{N}_0+k_{y}y& \text{contact}\\ 0& \text{separation}\end{array}$$

(28)where δ= ±1 is a sign function representing the direction of tangential force interaction $f_x$ . The force $f_x$ reaches the critical value when transition of contact from stick to slip occurs. Moreover, slip to stick transition occurs only at conditions when derivation of tangential force ${\dot{f}}_x=\mu {\dot{f}}_y$ is satisfied.

For a friction interface element, stiffness matrix K _f can be derived analytically by differentiating vectors of Fourier expansion coefficients for tangential, F _x , and normal, F _y , forces with respect to vectors X and Y, namely, the term $\partial {\mathbf{F}}_{\text{N}}/\partial \mathbf{Q}$ in equation (27). As there is no relation between force F _y and displacement vectors X, the expression of stiffness matrix K _f is given in the form

$${\mathbf{K}}_f=\left[\begin{array}{cc}\frac{\partial {\mathbf{F}}_x}{\partial \mathbf{X}}& \frac{\partial {\mathbf{F}}_x}{\partial \mathbf{Y}}\\ 0& \frac{\partial {\mathbf{F}}_y}{\partial \mathbf{Y}}\end{array}\right]$$

(29)

Furthermore, in stick types, the terms in equation (29) can be written in the following form

$$\frac{\partial {\mathbf{F}}_x}{\partial \mathbf{X}}={\displaystyle \sum _{j=1}^{n_t}\left(k_x{\mathbf{W}}_j+w_j{\left(\frac{\partial c_j}{\partial \mathbf{X}}\right)}^{\text{T}}\right)},\frac{\partial {\mathbf{F}}_x}{\partial \mathbf{Y}}={\displaystyle \sum _{j=1}^{n_t}{\mathbf{w}}_j{\left(\frac{\partial c_j}{\partial \mathbf{Y}}\right)}^{\text{T}}},\frac{\partial {\mathbf{F}}_y}{\partial \mathbf{Y}}={\displaystyle \sum _{j=1}^{n_t}{k}_y{\mathbf{W}}_j}$$

(30)where matrix ${\mathbf{W}}_j$ is integral of product of vector $\mathbf{H}$ used for transformation from frequency domain to time domain in equation (3) and vector $\mathbf{H}\prime$ for transformation backward, vector ${\mathbf{w}}_j$ is integral of vector $\mathbf{H}\prime$ , $n_t$ is number of time steps of Fourier expansion, and $c_j$ is a constant calculated by $c_j=f_x^{\mathrm{\hspace{0.17em}0}}\left(t_j\right)-k_{x}x\left(t_j\right)$ . ²² Correspondingly, in slip types, equation (30) can be rewritten in the form

$$\frac{\partial {\mathbf{F}}_x}{\partial \mathbf{X}}=0,\frac{\partial {\mathbf{F}}_x}{\partial \mathbf{Y}}={\displaystyle \sum _{j=1}^{n_t}\delta \mu k_y{\mathbf{W}}_j},\frac{\partial {\mathbf{F}}_y}{\partial \mathbf{Y}}={\displaystyle \sum _{j=1}^{n_t}{k}_y{\mathbf{W}}_j}$$

(31)

In addition, all terms are equal to 0 when contacts separated. The derivatives of constant term $c_j$ with respect to $\mathbf{X}$ and $\mathbf{Y}$ are devised by differentiation of slip to stick transition conditions as in Petrov and Ewins. ²²

It should be noted that, during calculating the solution of interaction forces in time domain, the state, contact or separation, is necessary to be judged in the periodic motion in terms of normal displacement y. When surfaces remain contact state, the convergence values of constant terms in equation (28) are determined in time domain in terms of the slipstick transition, to calculate the variable frictional force in the vibration period interval. However, for the state of contact separation, the time intervals of contact and separation need to be ascertained accurately. The friction interface force and stiffness are determined in contact period by the initial displacement based on contact point at initial contact time.

Finite element model of mistuned blades

In order to analyze the impact of mistuning on vibration response of bladed disk with friction and contact, a finite element model of turbine straight-plate bladed disk comprising 30 blades with damped shrouds is analyzed, as demonstrated in Figure 2. Due to mistuning factor, the calculation of nonlinear forced response and complex modes is performed on the full blade structure. The DOFs of blade roots are considered to be fixed to that of disk due to locked coupling effect of friction on negligible motion of root under huge centrifugal force. The main dissipation of energy for bladed disk system comes from the friction interaction at adjacent shroud contacts of blades. Node–node friction interface elements are built at the interaction contacts of damped shrouds, and relative displacements of surfaces subjected to external exciting force can be determined. The corresponding nodes of contact elements are assumed to be applied with mutual equivalent frictional force. The initial normal force acting on a contact interface is averaged on all frictional elements during calculation for convenience. In present cases, the top 300 orders of mode shape are taken as main DOFs remained and the residual condition ǁR(Q)ǁ ≤ 10⁻⁵ is used to check the convergence of the iteration process.

Figure 2.

Finite element model of mistuned bladed disk.

The stiffness matrix K _f derived from the friction model need to be transformed into a 3D derivative matrix ${\mathbf{\text{K}}}_f\text{`}$ by means of coordinate transformation

${\mathbf{\text{K}}}_f\text{`}=\mathbf{R}{\mathbf{K}}_f^{}{\mathbf{R}}^{\mathbf{T}}$

(32)

$$\mathbf{R}=\left[\begin{array}{cc}\sin \hspace{0.17em}{\theta }_1\cos \hspace{0.17em}{\theta }_2& \cos \hspace{0.17em}{\theta }_1\cos \hspace{0.17em}{\theta }_2\\ \sin \hspace{0.17em}{\theta }_1\sin \hspace{0.17em}{\theta }_2& \cos \hspace{0.17em}{\theta }_1\sin \hspace{0.17em}{\theta }_2\\ \cos \hspace{0.17em}{\theta }_1& -\sin \hspace{0.17em}{\theta }_1\end{array}\right]$$

(33)where ${\theta }_1$ is the angle between friction interface and axis, and ${\theta }_2$ is the coordinate angle of interaction element in cylindrical coordinate system.

The full circle finite element model consists of 129,960 elements and 137,857 nodes. Material parameters of the blades are set as follows: the elastic modulus is 210 GPa, the density is 7850 kg/m³, and Poisson’s ratio is 0.3. Constraints on the structure are set as follows: axial and tangential displacement constraints are set on the DOFs of disk inlet, and tangential displacement and axial coupling constraints are set on the DOFs of disk outlet. Excitation by a first engine-order harmonic is taken in account and the linear damping loss factor was taken to be 0.003 due to material and aerodynamic damping, which is a typical value in reality. ²⁴

Results and discussion

By using the method proposed in this paper for solving forced vibration of the large-scale finite element model of damped blades, the forced response of bladed disk assembly in Figure 2 is analyzed. In general, the higher the nodal diameter, the greater phase difference in the vibration of adjacent blades, and the calculation results for comparing are more obvious. For the sake of generality, a traveling sine wave excitation with only one nodal diameter is applied on surface profile of blades uniformly. The exciting force can be decomposed into components that are tangential and normal to the friction surface, and the component in the normal direction also arouses relative motion between damped shrouds. In the following sections, the coupled vibration characteristics of the mistuned turbine blades subjected to excitation is focused on. The forced response of blades is accompanied by cases of frictional slipstick, normal separation–contact transition, and other nonlinear behaviors. Besides, in order to evaluate the influence of circumferential motion on axial vibration, a reasonable fitting equation of variable normal force is addressed based on the results of normal motion.

Number of harmonic

The first thing for expand the nonlinear response solution in equation (2) by multi-harmonic balance method is the selection of appropriate number Nh of retained harmonic determined by efficiency and computational accuracy. For a nonlinearity period such as friction, a sufficient number of harmonics and time step must be maintained to accurately describe the state transitions that occur at the interface. The calculation setting is as follows: a constant normal force on shroud interface is F_y = 100 N, the tangential friction stiffness of the contact K_x = 1×10⁷ N/m, the normal stiffness K_y = 1.4×10⁷ N/m, the axial excitation force F_z = 10 N. The relative motions between two frictional DOFs for different values of Nh are plotted in Figure 3. The increase of the number of Nh causes that the convergence of the relative displacement calculated by HBM tends to the real solution yielded by time integration. The accuracy of the relative motion approximated directly affects the quality of response amplitude, as shown in Figure 3(b). For comprehensive consideration, 15 harmonics are retained to obtain higher calculation accuracy in this study.

Figure 3.

Influence of harmonics number retained on the accuracy of (a) relative displacement and (b) amplitude.

The harmonic spectrum of the multi-harmonic vibration (for the case in Figure 3) is shown in Figure 4. Overall, the first-order harmonics coefficient takes the dominant position, and the other harmonics coefficients maintain significant amplitudes in the resonance region, which demonstrates the importance of keeping multi-order harmonics for the accuracy of nonlinear resonance solution. In addition, the zero harmonics reflects the relative displacement of contact points of adjacent shrouds during motion. The relative displacement, regarded as a variable gap in this paper, reflecting the variation of normal force at shroud friction interface simultaneously. The corresponding formula is as follows: F_y = −g₀ · K_y (g ≤ 0), but F_y = 0 when initial gap g₀ > 0. In the following, an initial gap g₀ is taken into account to prescribe the preload normal force.

Figure 4.

Amplitudes of harmonic components.

Nonlinear behavior

In the realistic damped bladed disk, an addition of normal motion of shroud contact interfaces can alter, to some extent, the dynamic response of friction process, thus giving rise to more complex friction nonlinearity. Here, the nonlinear vibration in ideal tuned state is analyzed first to be as a reference for cases of mistuning. The normal direction of shroud interface is the circumferential direction of system coordinate, and the tangential direction is the axial direction of system coordinate. The calculation conditions are set as follows: initial contact gap g₀ = −0.002 mm (corresponding to preload normal force F_y = 28 N), excitation force decomposed into axial force F₁ = 10 N and tangential force F₂ = 2 N. The result of forced response is compared with two linear cases: (1) for independent freestanding vibration without dampers and (2) for a full locked contact between shrouds, as shown in Figure 5. It can be seen from the results for g₀ = −0.002 mm, friction damping is demonstrated to be significantly helpful for reducing the amplitude of axial resonance, and a nonlinear jump phenomenon occurs in the response curve of circumferential vibration. The first vibration mode belongs to axial vibration and the second vibration mode belongs to circumferential vibration. The normal motion on contact surface directly affects the variation of the normal force, and further contributes to the axial vibration in the range of circumferential resonance.

Figure 5.

Nonlinear forced responses with coupled vibration.

Figure 6 illustrates the coupled axial and circumferential forced responses for initial gap g₀ = −0.002 mm and g₀ = 0 mm, respectively. For axial vibration, the response curve for the case of constant normal force has nearly no significant differences compared with that for the case of variable force, which mainly because the circumferential amplitude is too small. In Figure 6(b), instead of the separation of blade interfaces for the case of constant normal force, a coexistence state of periodic transformation between contact and separation of interfaces is remained for the case of g₀ = 0 mm. Besides, a complex vibration behavior is exhibited in a low frequency range (100–350 Hz), known for existence of super-harmonic resonances due to partial separation, and a high frequency range (950–1080 Hz) caused by the effect of coupling vibration.

Figure 6.

Comparison of forced responses for nonlinear variable and constant normal forces (a) for g₀ = −0.002 mm and (b) for g₀ = 0 mm.

In order to obtain the optimal value of preloaded normal force, Figure 7 shows frequency response curves calculated at axial vibration with different initial values of initial shroud gaps, from −0.05 mm to −0.0002 mm. Finally, a minimum resonance peak is obtained with an initial shroud gap g₀ = −0.005 mm, which is retained in further analysis in consideration of its particularity.

Figure 7.

Forced responses for different shroud gaps.

Mistuned parameter

Contact mistuning occurs when contact stiffness, normal load, and other nonlinear contact state variables change for some blades, and eventually leads to the redistribution of vibration energy of the system. For a mistuned contact case, each tuned contact stiffness or shroud gap and other variable is multiplied by a scaling factor to obtain the mistuned variable.

Random and scattered are the typical characteristics of realistic mistuned parameters of bladed disk, and are hard to eliminate. Due to contact of dampers, vibrations of damped blades are coupled with others, leading to a complex nonlinear performance for vibration of the whole bladed disk. The effect of coupling vibration with friction and contact is the main target of this research. In this section, some mistuning factors for parameters are added to the finite element model by a random deviation pattern on structure. For a given parameter in structure, the vector of mistuned parameter ${\mathbf{P}}_{\text{m}}$ is derived from the tuned value $P_{\text{t}}$ by adding a random error submitting Gaussian distribution, which can be written as

$${\mathbf{P}}_{\text{m}}=P_{\text{t}}\left(1+L\mathbf{\xi }\right)$$

(34)where ξ is a perturbation vector following a Gaussian distribution with an average of 0 and a standard deviation of 1, and L is a coefficient regarded as the mistuning level by multiplication. In the following analysis, a random pattern of deviation vector ξ is taken and the arrangement is shown in Figure 8. Based on these deviations, three kinds of structure parameters, including tangential friction stiffness, contact initial gap, and structural stiffness of blade, are assumed to be mistuned for a detailed research on the complex coupled vibration of mistuned bladed disk with friction dampers.

Figure 8.

A random pattern of deviation vector ξ obtained.

Mistuning by a scatter of tangential friction stiffness

Frequency response curves of mistuned blades due to scatter of tangential friction stiffness for initial contact gap g₀ = −0.005 mm and g₀ = 0 mm are demonstrated in Figure 9. The mistuning level L is setting to 0.3, and the variation range of mistuned values calculated for tangential friction stiffness is [5.53×10⁶–1.71×10⁷] N/m. It can be seen that hardly any obvious scatter of response amplitude can be observed in the range of whole frequency. To illustrate the detailed motion of mistuned blades, Figure 9 also shows the hysteresis curves of blades at axial resonance under conditions of different initial gaps. The mistuning by scatter of tangential friction stiffness brings about the change of relative slip displacement between adjacent shrouds, but has little influence on response amplitude of mistuned blades. In addition, for the case of initial gap g₀ = 0, it can be clearly seen that the frictional force becomes 0 in half cycle of the motion due to the separation of contact surfaces.

Figure 9.

Forced responses of all blades with mistuning by tangential friction stiffness (a) for g₀ = −0.002 mm and (b) for g₀ = 0 mm.

Mistuning by a scatter of initial shroud gap

The calculation conditions for the case of mistuning by a scatter of shroud gaps are set as follows: the mean value of initial gap g₀ = −0.005 mm and mistuning level L = 0.3. In accordance with the above conditions, the range of initial gap value is obtained as [−0.0085 mm to −0.0028 mm]. Figure 10 shows the forced response curves of the bladed disk and also demonstrates the hysteresis curve at axial resonance frequency A and variation range of normal force between contact surfaces of blade shrouds at circumferential resonance frequency B, respectively. It can be seen that the response curves of blades shows obvious difference at axial and circumferential resonances. In order to intuitively observe the change of response curves for shroud gap mistuning, two frequency response functions (FRFs), for 18th and tuned blades respectively, are also plotted at circumferential resonance region. A bifurcation phenomenon for the circumferential resonance occurs due to coupled nonlinear characteristic, or more precisely, the partial separation of shroud contact surfaces at resonance. The nonlinear bifurcation behavior from one blade arouse the chain reaction for remaining blades in whole system vibration, and intuitively, it seems that curve appears fold back for multiple times.

Figure 10.

Forced responses of blades mistuned by shroud gap (a) response curves, (b) hysteresis curves at frequency A, and (c) variation range of normal force at frequency B.

Mistuning by a scatter of structural stiffness

The deviation of natural frequencies of blades generated by structural stiffness mistuning of blades also gives rise to a complex coupled vibration. Figure 11 shows the forced FRFs of mistuned blades with mistuning level L = 0.05 and initial shroud gap g₀ = −0.005 mm. Obviously, a scatter occurs for response curves not only at tangential and normal resonances, but also in the range of 720–750 Hz and 820–850 Hz. A major influence of structural stiffness mistuning is the significant amplification of resonance amplitude of partial blades, which is illustrated through a comparison of maximum mistuned and tuned responses in Figure 11. Moreover, the contact state of shroud interface may turn to separate on account of a more large mistuning level or a small absolute value of initial gap.

Figure 11.

Forced responses for all blades mistuned by structural stiffness.

Amplification factors for mistuning

Amplification factor, also called normalized amplitude, is defined as the ratio of maximum amplitude of mistuned blades at resonance to the maximum amplitude of tuned blades, and enormously dependent on the mistuning level of mistuned parameters. In order to dissipate the vibrational energy caused by excitation, the component of excitation force along the friction interface is often designed to be maximized. The frictional bending vibration of blades, that is, axial vibration in this paper, greatly changes the amplification factor due to the disturbance of the mistuning factors for parameters. In this section, in order to assess the impact of the variation of mistuning level, an analysis of axial response of damped blades is performed under different degrees of mistuning level. Three above mistuned parameters are remained to obtain the corresponding amplification factors, and the results are shown in Figure 12.

Figure 12.

Amplification factor for different types of mistuned parameter.

For the mistuning by tangential friction stiffness, the increase of mistuning level has almost a negligible enlargement on the amplification of system vibration. However, for the mistuning by initial shroud gap, amplification factor increases continuously, more specifically, vibration localization continues to deteriorate. The current acceptance of the localization mechanism is generally based on the wave-transfer theory. Mistuning causes vibration reflections at the structural interface. The result of multiple reflections can limit vibration to a localized area of the structure, creating vibration localization phenomena. As the mistuning level increases, the reflected energy on the localized blade will also increase, resulting in a larger amplification. Moreover, for the mistuning by structure stiffness of blades, the amplitude amplification factor rises sharply and then decreases as mistuning level increases, and it can be explained by the fact that the separation of natural frequency causes an increase for resonance region, resulting in dispersing of vibration energy throughout the whole resonance region.

Effect of normal motion on frictional vibration

In Figures 9 to 11, the effect of circumferential motion of mistuned blades on axial vibration is less obvious, mainly due to the too low circumferential displacement caused by farness of circumferential resonance. Actually, because of the single problem of model established, the frequencies of first-order axial resonance and circumferential resonance are unchangeable, thus these two resonance frequencies cannot be close together. However, for the bladed disk working assemblies, one of the high-order circumferential resonance frequencies generally may be close to one of axial resonance. So, the possibility of an extreme circumferential movement at axial resonance should be taken into account. In the next context, a fitting equation is approved to simulate approximately the normal force (considering contact and separation). By coupling to axial vibration, the fitting equation can be used to analyze the impact of circumferential motion of mistuned blades on axial vibration. Meanwhile, a variate to represent normal displacement intensity is introduced for further analysis.

By using the method of curve fitting, the equation of variation of normal force can be achieved in time domain directly. And other influences of contact force, such as phase, on vibration are also taken into consideration to ensure consistent with actual results. Since the vibration response near circumferential resonance is mainly concerned, it can be assumed that the phases of contact forces between shrouds remain unchanged in this small resonance region. Figure 13 shows the amplitude variation of the normal force in the time domain at circumferential resonance (in Figure 2), which look like sine vibrations with phase differences. The expression of contact normal force N_m between mth and next blade shrouds can be given by equation (35)

$$N_m=N_0+A\sin \hspace{0.17em}\left(\omega t-\frac{2\left(m-k\right)\pi }{\text{30}}\right)m=1,\dots ,30$$

(35)where N₀ is the preload value of normal force, A is the variation amplitude, and k is the redundant phase difference. After fitting, k = 1.5 and N₀ = 70 is obtained. Besides, N_m = 0 when it is less than 0 on account of contact separation.

Figure 13.

Variation of contact force between blade shrouds.

Similarly, for other initial shroud gap, the variation of normal force in time domain at circumferential resonance can be fitted to simulate and calculate the response when resonances coincide. For g₀ = 0, the normal forces in Figure 9(b) is fitted. The preload normal force, N₀, determines the initial vibration with nonlinear contact, and the amplitude of load variation, A, affects the existence and translation of slip and viscous states of friction motion. To evaluate this impact of the variation of circumferential vibration, different amplitudes of load variation (A = 0, A = 0.5 N₀ and A = 1.5 N₀) with two different initial gap g₀ = −0.005 mm and g₀ = 0 mm are taken to calculate responses of mistuned blades at the axial resonance range, 550—750 Hz, respectively. Here, mistuning by the scatter of structural stiffness of blades is only taken. Envelopes of forced response extracted from results of blades are plotted, as shown in Figure 14. For comparison, FRFs of tuned systems are also plotted in both figures to evaluate the influence of mistuning on the vibration characteristics.

Figure 14.

Envelopes of maximum response amplitudes of mistuned bladed disk compared with tuned responses (a) for g₀ = −0.005 mm and (b) for g₀ = 0 mm.

For initial gap g₀ = −0.005 mm, regardless of mistuned or tuned blades, envelope of maximum response amplitudes is almost unchanged with the increase of A from 0 to 1.5N₀. This is because, in the presence of a large normal force with a finite range, the contact interface keeps in viscous state over the entire period of vibration. Even if the contact separation occurs, such as for the case of A = 1.5N₀, the impact on the axial vibration of blades is extremely limited due to the longstanding coupled vibration. However, for initial gap g₀ = 0 mm, the change of the amplitude of load variation brings about a large difference to response results. Regardless of whether the disk structure is tuned or mistuned, the resonance frequency gradually increases and the resonance amplitude decreases significantly with increasing amplitude A. The reason is that, as A increases, the nonlinear coupling effect between shroud surfaces enhances, and hence, effectively reduces the vibration. For this reason, the partial vibration induced by the mistuning of structure stiffness is also gradually weakened. In addition, it should be noted that, if the impact is analyzed for other cases, such as mistuning by initial shroud gap, the fitting curve needs to keep multi-harmonic to ensure the accuracy due to the chain reaction of nonlinear bifurcation.

In conclusion, the effect of circumferential vibration on axial vibration varies in different cases, including different initial normal forces and variation amplitudes. For the mistuned blades with nonlinear friction and contact, the turning time for dampers to contact under the centrifugal force means a low normal force generated on contact surfaces. The circumferential component of excitation force would have a great influence on the axial vibration in some cases, and which also can be reduced effectively by a large amplitude of normal force variation. While when the coupling effect between blades is enhanced further, that is, a higher initial normal force appears, the effect of circumferential vibration on axial vibration would be reduced to a very limited extent.

Conclusions

An effective method for analyzing forced response of a large finite element model of mistuned bladed disk with nonlinear friction and contact considered is proposed. The method is based on modal reduction method to reduce the analyzed model to be a matrix formulation including the main DOFs (for the nodes acted with friction, constraint, and excitation) and remaining slave DOFs. Meanwhile, through a process of coupling friction interface model to adjacent blade contacts, a nonlinear friction and contact connection suitable for the complex multidirectional excitation is added into the matrix. By using other efficient techniques for solution of nonlinear equations, such as harmonic balance method and Newton–Raphson method, the forced response of the mistuned bladed disk is calculated.

For compromising calculation accuracy and cost time, first 15 harmonics are selected to carry out multi-harmonic frequency domain analysis. The coupling of circumferential vibration (perpendicular to shroud interfaces) and axial vibration (parallel to the shroud interfaces) of bladed disk all shows typical nonlinear behavior, such as friction damping and contact separation.

Tangential friction stiffness, initial shroud gap, and the structural stiffness are chose as the mistuned parameters for turbine damped blades, and the effects of mistuning by these parameters on the axial and circumferential vibration characteristics of damped blades are analyzed. The results are as follows: under the normal excitation, the interfaces between shrouds may turn to separate at normal resonance, and a nonlinear snapback also appears for the frequency response curve. With the mistuning by friction tangential stiffness considered, response curves of blades for axial vibration turn to disperse. With the mistuning by initial shroud gap considered, contact separation of partial blade shrouds results in complex unstable vibrations at circumferential resonance. Meanwhile, the scattered distribution of initial normal force brings the scatter of axial response curves and the increase of maximum amplitude. The mistuning introduced by structure stiffness causes the splitting of resonance frequencies for blades, leading to a scatter of response curves of blades and the increase of maximum response amplitude. It is worth noting that an unstable nonlinear circumferential vibration may also be induced under a low normal force on shroud interface. Moreover, a further analysis is performed to reveal the influence of mistuning level of mistuned parameters on the amplification factor.

Finally, the influence of amplitude variation of normal motion on the axial vibration of mistuned blades is addressed. The results show that, regardless of mistuned or tuned cases, amplitude variation of normal force has little effect on resonance amplitude of axial vibration when the preload normal force is large enough. However, when the preload normal force is low or the contact between shrouds is separated, the maximum response amplitude of mistuned blades can be greatly reduced due to the enhanced friction damping.