# Dynamic analysis of a bladed disk with friction and fretting-wear in blade attachments

Loïc Salles, Laurent Blanc, Fabrice Thouverez, Aleksander M Gouskov, Pierrick Jean

► To cite this version:

Loïc Salles, Laurent Blanc, Fabrice Thouverez, Aleksander M Gouskov, Pierrick Jean. Dynamic analysis of a bladed disk with friction and fretting-wear in blade attachments. ASME Turbo Expo 2009: Power for Land, Sea, and Air, Jun 2009, Orlando, United States. pp.465-476, 10.1115/GT2009-60151. hal-02537384

HAL Id: hal-02537384

https://hal.science/hal-02537384

Submitted on 8 Apr 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
DYNAMIC ANALYSIS OF A BLADED DISK WITH FRICTION AND FRETTING-WEAR
IN BLADE ATTACHMENTS

Loïc Salles<sup>a,b,c*</sup>, Laurent Blanc<sup>a</sup>, Fabrice Thouverez<sup>a</sup>, Aleksander M. Gouskov<sup>b</sup> and Pierrick Jean<sup>c</sup>

(a) École Centrale de Lyon, Laboratoire de Tribologie et Dynamique des Systèmes
36, avenue Guy de Collongue, 69134 Ecully Cedex, France

(b) Moscow State Technical University named after Bauman, Moscow, Russia

(c) Snecma – Safran group, 77550 Moissy-Cramayel, France

E-mail: Loic.Salles@ec-lyon.fr

ABSTRACT

Assembled bladed disks have many contact interfaces (blade-disk joint, blade shrouds, friction dampers...). Because of relative displacements at these interfaces, fretting-wear occurs, which affects negatively the lifetime of the structure. Methods exist to predict fretting-wear in quasi-static analysis. However they don't predict all the phenomena observed in blade attachments on real industrial plants. This paper studies the assumption of a responsibility of dynamics for fretting-wear damage.

A numerical treatment of fretting-wear under vibratory load-
ing is proposed. The method is based on the Dynamical La-
grangian Frequency Time method. It models unilateral contact
through Coulomb's friction law. The basic idea is to separate
time in two scales, slow scale for tribological phenomena and fast
scale for dynamics. For a chosen number of periods of vibration,
a steady state is assumed and the variables are decomposed in
Fourier series. An Alternating Frequency Time procedure is per-
formed to calculate the non-linear forces. Then, a Hybrid Pow-
ell's algorithm is used as solver. A quasi-analytical expression of
the Jacobian matrix decreases the duration of calculations. This
expression is also used to predict new relative displacement at
the interfaces due to the increase of wear depth. This method is
similar to a prediction-correction method, with wear depth as the
term of continuation.

Numerical investigations on a bladed-disk with friction con-

tact interfaces illustrate the performances of this method and
show the coupling between dynamical and tribological pheno-
mena.

NOMENCLATURE

~ Multiharmonic vectors.

-T, -N Tangential and normal directions.

$p_N, p'_n$ Normal pressures.

$p_T, p'_T$ Tangential shears.

u, v Displacements.

g Initial gap.

w, **W** Wear depth, nodal wear depth vector.

-r Relative value.

**U** Time-domain nodal displacement vector.

**M, C, K** Mass, viscous damping and stiffness matrices.

**Z<sub>r</sub>** Reduced dynamic stiffness.

$F_c$ Nodal contact force.

λ, λ̃ Lagrangian multiplier in time and frequency domains.

INTRODUCTION

One of the most important sources of nonlinearities in as-
sembled bladed disks is the frictional contact at the interfaces.
These nonlinearities have often a positive aspect because they
introduce damping in the vibratory response. But relative dis-

<footnote>*Address all correspondence to this author.</footnote>
placements at the interface also have a damaging effect, since fretting-wear appears and can induce cracking.

Fretting-wear studies in blade roots are usually performed on the assumption of quasi-static loading. Thus, in the case of plane engines, only the take off or landing configurations are calculated [1] or tested [2]. Nevertheless, some industrial damage surveys lead one to think that wear also occurs at cruising rate, in spite of the fact that only micro-slidings due to vibrations are observed in such conditions. Wear could be explained by a coupling between vibrations and wear, and by the fact that fretting cycles are repeated during a much longer time than the durations of take off and landing stages.

Methods exist to study the vibratory response of bladed disks in the presence of frictional contact [3, 4]. But studies of fretting-wear under dynamical loading - including inertial effects - are seldom [5, 6]. Recently a new method [7] was proposed to couple the calculations of the vibratory response and of the wear kinetics based on the Dynamical Lagrangian Frequency Time method (DLFT) [4] and on a multiscale approach (these developments were also used [8] to calculate a kind of "modal wear" in synergy with the Nonlinear Normal Modes concept [9]). The present paper presents improvements to the method introduced in [7].

To model the coupling between fretting-wear and vibration, the continuous formulation of the problem is discretized through the finite elements method. The main assumption is that a periodic steady state is reached and that wear will modify this state only a little. Thus, it is possible to use methods based on Multi Harmonic Balance and Alternate Frequency Time procedures [10]. A multiscale formalism is introduced to distinguish between the "fast" phenomena caused by vibrations and the "slow" phenomena due to wear.

Iterative numerical techniques are necessary because a non-linear system must be solved, at each step of the slow scale. The first improvement to the DLFT-with-wear method [7] is the choice to proceed by "wear steps" instead of "long time steps", as previously. This way, the number of non-linear systems to solve is better controlled.

Here the strategy is based on a Hybrid Powell's solver. This algorithm requires to evaluate a Jacobian matrix. This step is very time consuming if the usual finite difference method is employed. An analytical estimation of the Jacobian is more efficient [10, 11]. The second improvement to the method is the implementation of a quasi-analytical Jacobian inspired by literature's analytical Jacobians. It is obtained by derivation of the non-linear forces expressed in the time domain according to the contact state. The dependancy of the Jacobian on wear is taken into account in the same manner. The two Jacobian matrices resulting from this operation allow one to predict an updated worn geometry in a very time-saving way.

The third improvement to the DLFT-with-wear method is

the proposal of a specific error estimator to control the number of required harmonics.

There is nothing new about wear quantification. Wear is a complex phenomenon because many causes can change wear debris creation: hardness, plasticity, grain structure, temperature... Here the Archard's model [12] is used on the local scale in the interface. It is assumed that the friction coefficient is unaffected by the evolution of wear. A model based on the dissipated energy [13] could have been used as well.

The improved DLFT-with-wear method is applied here to the industrial case of a bladed disk with dovetail joints. Its vibratory response and wear kinetics are calculated. The results provide a numerical evidence that worrying about the coupling between fretting-wear and dynamics is of particular relevance to designers.

# THE REFERENCE PROBLEM

## Continuous formulation

<figure><img src="image_2.png" /><figcaption>Figure 1: Description of the problem with two solids</figcaption></figure>

The formalism is the one described by Strömberg [14]. The constitutive equations are built for two elastic bodies with a frictional contact interface in situation of fretting-wear. Each solid $l$ occupies a domain $\Omega^l$ with smooth boundaries $\partial\Omega^l$ as drawn in Fig. 1. Each boundary $\partial\Omega^l$ is divided into three disjoint parts $\Gamma_l^f, \Gamma_l^u$ and $\Gamma_c$ (common to both solids). Traction forces $T^l$ are imposed on $\Gamma_l^f$, displacements on $\Gamma_l^u$ and frictional contact with fretting-wear conditions on $\Gamma_c$.

The constitutive law for the interface is given by Signorini's unilateral contact conditions and Coulomb's law of friction,
with $\mu$ the frictional coefficient. Wear is characterised by loss of matter. Archard's law is used. It is supposed to be true for an elementary volume. Then, a local relation is derived between normal pressure, tangential velocity and wear rate. An internal state variable $w$, the wear gap, is introduced so that the interface law is derived from a generalized standard material formulation. The friction coefficient $K_w$ is such that $K_w = \frac{K_a}{3H}$, where $K_a$ is the wear rate coefficient of Archard's law and $H$ is material's hardness in the contact zone. $K_w$ is deduced from experiments. For example, McColl [15] determines an averaged wear coefficient from the measured wear profile for a cylinder-flat contact (this averaged coefficient is then introduced in a discretized Archard's law to simulate the evolution of wear in a finite element "wear box" approach). Here $K_w$ stems from a flat punch contact configuration. Experiments have shown that the coefficient of friction changes with wear evolution. Indeed, it actually depends on the chimical composition of the contact area. Dick et al. [16] have proposed a model to take this into account. Here, material damage is neglected and $K_w$ is supposed not to evolve as wear increases.

Variationally, the problem is defined by three integral equations and one constitutive law (1-4). These equations stem from the principle of virtual work, the weak formulation of Signorini's unilateral contact conditions, the complementary law of Coulomb and the local formulation of Archard's law. More detailed explanations can be found in [14].

For each moment $t$ of the studied time interval $[0, T]$, are sought the displacement field $\mathbf{u}$ for each point $\mathbf{x}$ of both solids (with $\mathbf{u}_N$ and $\mathbf{u}_T$ its normal and tangential components) and the contact pressure field $\mathbf{p}$ (with $\mathbf{p}_N$ and $\mathbf{p}_T$ its normal and tangential components). $\mathcal{W}$ designates the dual of $w$ i.e. the driving force for wear. It has the same physical meaning as $\mathbf{p}_N$. $\mathbf{v}$, $\mathbf{p}'$ and $\mathcal{W}'$ designate the test-fields associated respectively with $\mathbf{u}$, $\mathbf{p}$ and $\mathcal{W}$.

$$ \int_{\Omega^l} \rho \ddot{u}_i v_i dV + \int_{\Omega^l} E_{ijkl} \frac{\partial u_k}{\partial x_l} \frac{\partial v_i}{\partial x_j} dV + \int_{\Gamma_c} p_i v_i dA - \int_{\Gamma_c^l} t_i v_i dA = 0 \quad \forall \mathbf{v} \in \mathcal{V}, \quad (1) $$

$$ \int_{\Gamma_c} (u_N - w - g)(p'_N - p_N) dA \le 0 \quad \forall p'_N \in \mathcal{K}_N, \quad (2) $$

$$ \int_{\Gamma_c} (\dot{u}_{T\alpha} (p'_{T\alpha} - p_{T\alpha}) + \dot{w} (\mathcal{W}' - \mathcal{W})) dA \le 0 \quad \forall (p'_{T}, \mathcal{W}') \in \mathcal{F}(p_N) \quad (3) $$

with

$$ \begin{aligned} \mathcal{V} &= \{ \mathbf{v} | \mathbf{v}(\mathbf{x}) = \mathbf{0}, \mathbf{x} \in \Gamma_u^l \}, \\ \mathcal{K}_N &= \{ p_N | p_N(\mathbf{x}) \ge 0, \mathbf{x} \in \Gamma_c \}, \\ \mathcal{F}(p_N) &= \{ (\mathbf{p}_T, \mathcal{W}) | (\mathbf{p}_T(\mathbf{x}), \mathcal{W}(\mathbf{x})) \in \mathcal{F}(p_N), \mathbf{x} \in \Gamma_c \}. \end{aligned} $$

$F(p_N)$ is a closed convex set describing the friction and wear limit criterion. It is such that

$$ F(p_N) = \{ \| \mathbf{p}_T \| - \mu | p_N | + K_w p_N \mathcal{W} - K_w p_N^2 \le 0 \}. $$

In (2), $g$ is the initial gap between both solids. In (3), the components of each vector are expressed in an orthonormal basis $\mathbf{n}_\alpha$, perpendicular to $\mathbf{n}_c$. $\dot{\mathbf{w}}$ is obtained through the following equation:

$$ \dot{\mathbf{w}} = K_w |p_N| \|\dot{\mathbf{u}}_T\| \quad (4) $$

## Finite element discretization

A FEM discretization of (1) is performed. In practice, this strategy can be coupled together with a component mode synthesis to reduce the size of the problem. Capital letters designate the FEM counterpart of the variables named by small letters before: they are vectors of nodal quantities. When both structures are in contact, the equations of motion for each structure are:

$$ \mathbf{M}\ddot{\mathbf{U}}^l + \mathbf{C}\dot{\mathbf{U}}^l + \mathbf{K}\mathbf{U}^l + \mathbf{F}_c^l(\mathbf{U}^l, \dot{\mathbf{U}}^l, \mathbf{W}^l) = \mathbf{F}_{ex}^l, \quad (5) $$

where $\mathbf{M}$, $\mathbf{C}$ and $\mathbf{K}$ respectively designate the mass, damping and stiffness matrices. For the sake of simplicity, exponent $l$ will be dropped afterwards.

$\mathbf{F}_{ex}$ is the vector of external forces it provides here a periodic excitation. $\mathbf{F}_c$ represents the non-linear contact forces due to friction and impact; they also depend on the wear and on the materials in contact.

Wear is calculated at each interface node by:

$$ \dot{\mathbf{W}}^M = \frac{K_w}{\alpha^M} |\mathbf{P}_N^M| \|\dot{\mathbf{U}}_T^M\|. \quad (6) $$

$\alpha^M$ is the weighting factor for node $M$. The weighting factors depend on the quadrature rule used to calculate the integrals on each elementary contact area. Exponent $M$ in $\mathbf{P}_N^M$ and $\dot{\mathbf{U}}_T^M$ designate the nodal quantities.

The constraints introduced by Eqs. (2) and (3) become:

$$ {}^t(\mathbf{U}_N - \mathbf{W} - \mathbf{G})(\mathbf{P}'_N - \mathbf{P}_N) \le 0 \quad \forall \mathbf{P}'_N \in \mathcal{K}_N^h, \quad (7) $$

$$ {}^t\dot{\mathbf{U}}_T (\mathbf{P}'_T - \mathbf{P}_T) \le 0 \quad \forall \mathbf{P}'_T \in \mathcal{F}^h(\mathbf{P}_N), \quad (8) $$
where $K_N^h = \{P_N^M | P_N^M \ge 0\}$ is the approximation of $K_N$ and $F^h(P_N) = \{P_T^M | \|P_T^M\| - \mu|P_N^M| \le 0\}$ is the approximation of $F(p_N)$ with $W^M = P_N^M$.

# BASIC ASPECTS OF THE DLFT-WITH-WEAR METHOD

## A Harmonic Balance Method with two time scales

The main idea is to split time into two different scales: a fast one associated with vibratory phenomena and a slow one associated with wear. Indeed, at the scale of a few cycles wear appears as an almost constant interface gap. It is then assumed that wear doesn't change the aspect of the periodic response during a short lapse of time. On this period it is then possible to describe displacements and contact forces with Fourier series. On a longer duration, Fourier coefficients will evolve as functions of a time variable linked with the fretting-wear. It is important to mention that, here, the depth of wear is very small compared with the characteristic dimensions of the structures in contact. That is why the modifications of the mass and stiffness matrices due to wear are neglected.

This procedure belongs to the family of multiscale methods, described by Meirovitch [17] for a single harmonic balance. The Fourier series of $U$ can be expressed as:

$$ U(\tau, \eta) = \tilde{U}_0(\eta) + \sum_{j=1}^{N_{ft}} (\tilde{U}_{n,c}(\eta)\cos(j\tau) + \tilde{U}_{n,s}(\eta)\sin(j\tau)), \quad (9) $$

where $\tau$ is the variable of the fast time scale and $\eta$ the one of the slow time scale. In the frequency domain, after condensation on the interface nodes, the equation of motion becomes:

$$ Z_r \tilde{U}_r(\eta) + \tilde{\lambda}(\eta) = \tilde{F}_r, \quad (10) $$

or

$$ f(\tilde{U}_r) = Z_r \tilde{U}_r(\eta) + \tilde{\lambda}(\eta) - \tilde{F}_r = 0, \quad (11) $$

where $\tilde{U}_r$, $\tilde{F}_r$ and $Z_r$ designate respectively the reduced multi-harmonic relative displacement vector, reduced stiffness matrix and reduced external excitation. $\tilde{\lambda}$ is the vector of the Lagrange multipliers which represent the contact forces in the frequency domain.

Algorithmically, the solution over the time interval $[0, T]$ is originally calculated through a step-by-step procedure involving an incrementation of the slow time variable: for each step, Eq. 11 is solved through a Newton-like solver; this process constitutes the external loop of the method.

## Modeling of the contact forces

Solving (10) requires to know the expression of $\tilde{\lambda}$. Unfortunately it is not possible to calculate it directly in the Galerkin procedure, indeed it depends on the state of each contact node - stick, slip or separation - which is a priori unknown. To undertake this difficulty it is common to use the Alternating Frequency Time method (AFT). Displacements and velocities are calculated in the frequency domain and transformed into the time domain using an inverse DFT procedure (iDFT). In the time domain, contact forces could be evaluated through different methods. The easiest one is to regularize the sign function - depending on the velocity in the evaluation of Coulomb's forces - by another function which is continuous. It allows a direct computation of the non-linear friction forces [18]. The use of a penalty method is another popular method [10, 19]. The additional stiffnesses may then represent a damper's stiffness or the contact asperities stiffness. In the case of fretting-wear these stiffnesses could change with wear process, to take the modification of material properties in the contact area into account.

Another method has been proposed by Nacivet et al. [4]: the Dynamic Lagrangian Frequency-Time method (DLFT). It uses augmented Lagrangians which allow to calculate without any softening of the non-smooth frictional contact law. A time-marching procedure in the time domain is also required. Compared to the conventional contact penalty method, the main advantage of the DLFT method is that, at convergence, results don't depend on any penalty coefficient (hence the term "pseudo-penalty" coefficient to designate the coefficient used below to enforce the matching between the time and frequency descriptions of displacements). This method has been successfully used to predict friction damping in blade attachments [3] and to quantify the efficiency of friction ring dampers [19]. In a precedent paper [7] this method has been coupled with a polynomial expansion of wear to calculate directly the wear evolution for a two-degrees-of-freedom model. Here the formalism for contact forces calculation is the same as in [7] (but wear depth is evaluated in a less complex manner, explained farther).

In the frequency domain, the Lagrange multiplier $\tilde{\lambda}$ is formulated as a pseudo-penalization of the equations of motion on the tangential and normal directions:

$$ \tilde{\lambda}^T = \tilde{F}_r^T - (Z_r \tilde{U}_r)^T + \epsilon_T (\tilde{U}_r^T - \tilde{X}_r^T), \quad (12a) $$

$$ \tilde{\lambda}^N = \tilde{F}_r^N - (Z_r \tilde{U}_r)^N + \epsilon_N (\tilde{U}_r^N - W^N - \tilde{X}_r^N). \quad (12b) $$

$\epsilon_T$ and $\epsilon_N$ are pseudo-penalty coefficients, $\tilde{X}_r$ is a new vector of relative displacements, which is computed in the time domain. The pair $(\tilde{\lambda}, \tilde{X}_r)$ is determined through an AFT procedure.
dure. Equation (11) becomes:

$$f(\tilde{U}_r) = \varepsilon (\tilde{U}_r - W - \tilde{X}_r). \quad (13)$$

The convergence ensures that the time domain $\tilde{X}_r$ and frequency domain $\tilde{U}_r$ match with respect for contact conditions.

The AFT is based on a prediction/correction procedure in the time domain (summarized in Fig. 2). Algorithmically, this is the internal loop of the method.

<figure><img src="image_3.png" /><figcaption>Figure 2: Computation of the Lagrangian vector with wear</figcaption></figure>

The contact forces are calculated in the time domain, where the transition criteria between the three possible states are easily formulated. Equation (12) is reformulated as:

$$\tilde{\lambda} = \underbrace{\tilde{F}_r - Z_r \tilde{U}_r + \varepsilon (\tilde{U}_r - W)}_{\tilde{\lambda}_u} - \underbrace{\varepsilon \tilde{X}_r}_{\tilde{\lambda}_x}, \quad (14)$$

The period is split into N time steps. $\tilde{\lambda}$, $\tilde{\lambda}_u$ and $\tilde{\lambda}_x$ have $\{\tilde{\lambda}^n\}_{n=1..N}$, $\{\tilde{\lambda}_u^n\}_{n=1..N}$ and $\{\tilde{\lambda}_x^n\}_{n=1..N}$ as respective time domain counterparts. These vectors are obtained from the frequency domain vectors through an iDFT procedure. A prediction/correction is then used to compute the contact forces. At each time increment it first assumes that the contact node is in stick situation, thus the node doesn't move: $\tilde{\lambda}_x^{n,T} = \tilde{\lambda}_x^{n-1,T}$ and

$\tilde{\lambda}_x^{n,N} = 0$. The predicted contact forces are:

$$\tilde{\lambda}_{pre}^{n,T} = \tilde{\lambda}_{u}^{n,T} - \tilde{\lambda}_{x}^{n-1,T}, \quad \tilde{\lambda}_{pre}^{n,N} = \tilde{\lambda}_{u}^{n,N}. \quad (15)$$

The corrected contact forces will be:

$$\tilde{\lambda}^n = \tilde{\lambda}_u^n - \tilde{\lambda}_x^n, \quad (16)$$

and $\tilde{\lambda}_x^n$ will we be calculated to satisfy the contact and friction laws.

1. Separation: $\tilde{\lambda}_{pre}^{n,N} \ge 0$
The contact is lost and the forces should be zero.

$$\tilde{\lambda}_x^n = \tilde{\lambda}_u^n. \quad (17)$$

2. Stick: $\tilde{\lambda}_{pre}^{n,N} < 0$ and $\|\tilde{\lambda}_{pre}^{n,T}\| < \mu |\tilde{\lambda}_{pre}^{n,N}|$
In this case, the prediction verifies the contact conditions:

$$\tilde{\lambda}_x^{n,N} = 0, \quad \tilde{\lambda}_x^{n,T} = \tilde{\lambda}_x^{n-1,T}. \quad (18)$$

3. Slip: $\tilde{\lambda}_{pre}^{n,N} < 0$ and $\|\tilde{\lambda}_{pre}^{n,T}\| \ge \mu |\tilde{\lambda}_{pre}^{n,N}|$
Again, there is no normal relative displacement. The correction is made assuming that the tangential contact force has the same direction as the tangential predicted force. The definition of relative velocity and the respect of the Coulomb's law give:

$$\tilde{\lambda}_x^{n,N} = 0, \quad \tilde{\lambda}_x^{n,T} = \tilde{\lambda}_x^{n-1,T} + \tilde{\lambda}_{pre}^{n,T} \quad 1 - \mu \frac{|\tilde{\lambda}_{pre}^{n,N}|}{\|\tilde{\lambda}_{pre}^{n,T}\|} \quad (19)$$

The final step consists of transforming back the time domain updated Lagrangian in the frequency domain using the DFT algorithm.

At the end of this internal loop, it becomes possible to calculate the nodal wear for a single fretting cycle, denoted $\delta W_1^M$, by wear rate integration:

$$\delta W_1^M(\eta) = \frac{K_w}{\alpha^M} \int_{\eta}^{\eta+T_{exc}} |P_N^M(\tau, \eta)| \|\dot{U}_T^M(\tau, \eta)\| d\tau, \quad (20)$$

where $T_{exc}$ is the time period of the harmonic excitation force (it is also the period at which the system is supposed to respond). The interface gap - i.e. the geometry - can then be updated. At the level of the contact, equilibrium conditions are then modified: this imposes a return back to the beginning of the prediction/correction process to re-equilibrate them. Then, a new loop - the wear updating loop or intermediate loop - appears in the algorithm. This is specific to the DLFT-with-wear method.
# IMPROVEMENTS TO THE DLFT-WITH-WEAR METHOD

## About the wear depth calculation

Previously, an arbitrary number of cycles was chosen before the geometry was updated. To decrease the number of the slow time steps, it is necessary to choose a more efficient cycle jump strategy.

Here the proposal consists of incrementing the wear depth instead of time. A maximum authorized total wear depth, denoted by $W_{tot}$, is heuristically set for the whole considered duration up to $T$. $W_{tot}$ is split into $k$ intermediate wear steps. A too important wear step would provoke such an unbalance in the contact zone that a high number of loops in the prediction/correction process would be required to compensate for it. $k$ is chosen heuristically as the probably best compromise between the total number of wear steps and the number of required prediction-correction loops.

The relation between $W_{tot}/k$ and the number of jumped cycles $\Delta N$ is:

$$ \Delta N(\eta) = \frac{W^{tot}}{k \max_{M}(\delta W_{1}^{M}(\eta))}. \qquad (21) $$

If $l$ denotes the considered step of the external loop - still performed on the slow time scale - and $\eta^l$ the time of the beginning of this step, the vector of nodal wears at the beginning of the next step, $\Delta N(\eta)$ cycles later, is obtained by:

$$ \mathbf{W}(\eta^{l+1}) = \mathbf{W}(\eta^l) + \Delta N(\eta^l) \delta \mathbf{W}_1(\eta^l) \qquad (22) $$

If $\Delta N(\eta)$ becomes "very" important, the algorithm stops. Physically, it means that the wear rate becomes too small and that a stationary state may have been reached.

## About the Jacobian matrix calculation

In the external loop, the use of a Newton-like solver - here a Hybrid Powell solver - requires the evaluation of a Jacobian matrix. It is possible to evaluate it numerically by finite differences but its time-consuming (this strategy was originally used in the DLFT-with-wear method). Petrov [10] has shown that it's possible to calculate analytically the Jacobian in the case of a penalty formulation of frictional contact. Here, it's proposed to evaluate quasi-analytically the Jacobian matrix in the case of the Dynamic Lagrangian formulation. The dependencies of $f$ on respectively $\tilde{U}_r$ and $W$ are both considered with the introduction of the corresponding Jacobians $\mathbf{J}$ and $\mathbf{J}_w$.

The Jacobian with respect to $\tilde{U}_r$ is defined as follows:

$$ \mathbf{J} = \frac{\partial f(\tilde{U}_r)}{\partial \tilde{U}_r}. \qquad (23) $$

$f$ consists of two parts, a linear one and a nonlinear one:

$$ \mathbf{J} = \mathbf{Z}_r + \frac{\partial \tilde{\lambda}}{\partial \tilde{U}_r} = \mathbf{Z}_r + \frac{\partial \tilde{\lambda}_u}{\partial \tilde{U}_r} - \frac{\partial \tilde{\lambda}_x}{\partial \tilde{U}_r}. \qquad (24) $$

First, $\tilde{\lambda}_x$ is formulated using an exact integration on the vibratory period. This period is split in several time intervals according to the state of the contact. The idea is the same as in [10], but the present approach differs from Petrov's one because it will finally use a time discrete algorithm. The strategy is as follows:

$$ \tilde{\lambda}_{x}^{(0)} = \sum_{k=1}^{m} \int_{\tau_k}^{\tau_{k+1}} \lambda_{x}(\tau) d\tau, \qquad (25a) $$

$$ \tilde{\lambda}_{x}^{(j,c)} = \sum_{k=1}^{m} \int_{\tau_k}^{\tau_{k+1}} \lambda_{x}(\tau) \cos(j\tau) d\tau, \qquad (25b) $$

$$ \tilde{\lambda}_{x}^{(j,s)} = \sum_{k=1}^{m} \int_{\tau_k}^{\tau_{k+1}} \lambda_{x}(\tau) \sin(j\tau) d\tau, \qquad (25c) $$

where $m$ is the number of the different contact states. $\tau = \tau_k$ and $\tau = \tau_{k+1}$ are the transition moments. Using Leibniz's rule, the derivatives of these functions with respect to $\tilde{U}_i$ are, for the cosine terms (the expressions of derivatives of $\tilde{\lambda}_x^{(0)}$ and $\tilde{\lambda}_x^{(j,s)}$ are similar):

$$ \frac{\partial \tilde{\lambda}_x^{(j,c)}}{\partial \tilde{U}_i} = \sum_{k=1}^{m} \left( \frac{\partial \tau_{k+1}}{\partial \tilde{U}_i} \lambda_x(\tau_{k+1}) \cos(j\tau_{k+1}) - \frac{\partial \tau_k}{\partial \tilde{U}_i} \lambda_x(\tau_k) \cos(j\tau_k) + \int_{\tau_k}^{\tau_{k+1}} \frac{\partial \lambda_x(\tau)}{\partial \tilde{U}_i} \cos(j\tau) d\tau \right) \qquad (26) $$

Since $\lambda_x$ is a continuous and periodic function, proportional to displacement $\mathbf{X}_r$, a simplification occurs:

$$ \sum_{k=1}^{m} \left( \frac{\partial \tau_{k+1}}{\partial \tilde{U}_i} \lambda_x(\tau_{k+1}) \cos(j\tau_{k+1}) - \frac{\partial \tau_k}{\partial \tilde{U}_i} \lambda_x(\tau_k) \cos(j\tau_k) \right) = 0. \qquad (27) $$

And finally the gradients of (25) are:

$$ \frac{\partial \tilde{\lambda}_x^{(j,0)}}{\partial \tilde{U}_i} = \sum_{k=1}^{m} \int_{\tau_k}^{\tau_{k+1}} \frac{\partial \lambda_x(\tau)}{\partial \tilde{U}_i} d\tau \qquad (28a) $$

$$ \frac{\partial \tilde{\lambda}_x^{(j,c)}}{\partial \tilde{U}_i} = \sum_{k=1}^{m} \int_{\tau_k}^{\tau_{k+1}} \frac{\partial \lambda_x(\tau)}{\partial \tilde{U}_i} \cos(j\tau) d\tau \qquad (28b) $$

$$ \frac{\partial \tilde{\lambda}_x^{(j,s)}}{\partial \tilde{U}_i} = \sum_{k=1}^{m} \int_{\tau_k}^{\tau_{k+1}} \frac{\partial \lambda_x(\tau)}{\partial \tilde{U}_i} \sin(j\tau) d\tau \qquad (28c) $$
In time domain, the gradient of $\lambda_x$ is calculated according to the contact states defined by Eqs. (17)-(19). The three components of vector $\lambda_x^n$ are $\lambda_x^{N,n}$ in the normal direction and $\lambda_x^{T,n}$ in the tangential direction (with $\lambda_x^{x,n}$ and $\lambda_x^{y,n}$ its components in an orthonormal basis of the tangential contact plane, perpendicular to $n_c$).

$$ \frac{\partial \lambda_x^{N,n}}{\partial \tilde{U}_r} = \begin{cases} \frac{\partial \lambda_u^{N,n}}{\partial \tilde{U}_r} & (separation) \\ 0 & (contact) \end{cases} \qquad (29a) $$

$$ \frac{\partial \lambda_x^{x,n}}{\partial \tilde{U}_r} = \begin{cases} \frac{\partial \lambda_u^{x,n}}{\partial \tilde{U}_r} & (separation) \\ \frac{\partial \lambda_x^{x,n-1}}{\partial \tilde{U}_r} & (stick) \\ \frac{\partial \lambda_x^{x,n-1}}{\partial \tilde{U}_r} + \frac{\partial \lambda_{pre}^{x,n}}{\partial \tilde{U}_r} & 1 - \mu \frac{|\lambda_u^{N,n}|}{\|\lambda_{pre}^T\|} \end{cases} + \frac{\mu}{\|\lambda_{pre}^T\|} \frac{\partial \lambda_u^{N,n}}{\partial \tilde{U}_r} \frac{\lambda_{pre}^{x,n}}{\|\lambda_{pre}^{T,n}\|} + \lambda_{pre}^{x,n} \frac{\lambda_x^{x,n} \cdot \frac{\partial \lambda_u^{x,n}}{\partial \tilde{U}_r} + \lambda_x^{y,n} \cdot \frac{\partial \lambda_u^{y,n}}{\partial \tilde{U}_r}}{\|\lambda_{pre}^{T,n}\|^3} |\lambda_u^{N,n}| \quad (29b) $$

$$ \frac{\partial \lambda_x^{y,n}}{\partial \tilde{U}_r} = \begin{cases} \frac{\partial \lambda_u^{y,n}}{\partial \tilde{U}_r} & (separation) \\ \frac{\partial \lambda_x^{y,n-1}}{\partial \tilde{U}_r} & (stick) \\ \frac{\partial \lambda_x^{y,n-1}}{\partial \tilde{U}_r} + \frac{\partial \lambda_{pre}^{y,n}}{\partial \tilde{U}_r} & 1 - \mu \frac{|\lambda_u^{N,n}|}{\|\lambda_{pre}^T\|} \end{cases} + \frac{\mu}{\|\lambda_{pre}^T\|} \frac{\partial \lambda_u^{N,n}}{\partial \tilde{U}_r} \frac{\lambda_{pre}^{y,n}}{\|\lambda_{pre}^{T,n}\|} + \lambda_{pre}^{x,n} \frac{\lambda_x^{x,n} \cdot \frac{\partial \lambda_u^{x,n}}{\partial \tilde{U}_r} + \lambda_x^{y,n} \cdot \frac{\partial \lambda_u^{y,n}}{\partial \tilde{U}_r}}{\|\lambda_{pre}^{T,n}\|^3} |\lambda_u^{N,n}| \quad (29c) $$

with

$$ \frac{\partial \lambda_{pre}^{T,n}}{\partial \tilde{U}_r} = \frac{\partial \lambda_u^{T,n}}{\partial \tilde{U}_r} - \frac{\partial \lambda_x^{T,n-1}}{\partial \tilde{U}_r}. \quad (30) $$

Then a Discrete Fourier Transform (DFT) of $\left\{ \frac{\partial \lambda_x^n}{\partial \tilde{U}_r} \right\}_{n=1..N}$ gives

$\frac{\partial \tilde{\lambda}_x}{\partial \tilde{U}_r} \cdot \frac{\partial \lambda_u^n}{\partial \tilde{U}_r}$ is obtained by:

$$ \frac{\partial \lambda_u^n}{\partial \tilde{U}_r} = i \text{DFT}(-\mathbf{Z}_r + \varepsilon \mathbf{I}). \quad (31) $$

where $\mathbf{I}$ designates the identity matrix, with the same size as $\mathbf{Z}_r$.

Finally, the introduction of this expression in (24), gives the Jacobian of function $f$:

$$ \mathbf{J} = \varepsilon \mathbf{I} - \frac{\partial \tilde{\lambda}_x}{\partial \tilde{U}_r}. \quad (32) $$

The gradient of $f$ with respect to $W$, denoted $\mathbf{J}_W$ is:

$$ \mathbf{J}_W = -\varepsilon \mathbf{I}_{M_{max}} + \frac{\partial \tilde{\lambda}_x}{\partial W} \quad (33) $$

where $\mathbf{I}_{M_{max}}$ is identity matrix, with the same size as the total number of interface nodes $M_{max}$.

These gradients are used here to predict the multiharmonic displacement vector for an updated worn geometry. First, $f$ is expanded as a Taylor series of the first order, with $\tilde{U}_r^{l+1} = \tilde{U}_r^l + \delta \tilde{U}_r$ and $W^{l+1} = W^l + \delta W$. $(\tilde{U}_r^{l+1}, W^{l+1})$ and $(\tilde{U}_r^l, W^l)$ are supposed to be solutions of the nonlinear problem respectively at times $\eta^{l+1}$ and $\eta^l$. Then

$$ f(\tilde{U}_r^{l+1}, W^{l+1}) = f(\tilde{U}_r^l, W^l) + \mathbf{J} \delta \tilde{U}_r + \mathbf{J}_W \delta W = 0, \quad (34) $$

$\delta W$ being known from the wear rate integration, the expression for $\delta \tilde{U}_r$ is:

$$ \delta \tilde{U}_r = -\mathbf{J}^{-1} \mathbf{J}_W \delta W \quad (35) $$

The prediction-correction algorithm is then used, with $U_r^{l+1}$ as initial vector. This method leads to a significant decrease of the number of iterations in the Hybrid Powell's algorithm.

The program used for this paper is written in Matlab language. The Hybrid Powell used here is the one implemented in the `fsolve` function of the `Optimization Toolbox`.

## About the error estimation

The accuracy of the results with HBM/AFT methods depends on the number of harmonics, $N_h$. The bigger $N_h$ is, the slower the calculation is. Moreover, an increase of $N_h$ doesn't necessarily improve the accuracy a significant way.

Before the implementation of wear, $N_h$ was controlled, as proposed by Laxalde [20], by an estimation of the error between the contact forces corrected in the time domain and the contact forces expanded in time domain from the multiharmonic force
vector. The definition of the error in force for the DLFT method is:

$$ \epsilon_{F_c}(N_h) = \frac{\int_0^T |\lambda(\tau) - \sum_{j=0}^{N_h} (\tilde{\lambda}^c \cos(j\tau) + \tilde{\lambda}^s \sin(j\tau))| d\tau}{\int_0^T |\lambda(\tau)| d\tau} \quad (36) $$

The present work proposes to expand this method to the dissipated energy. Indeed, in case of fretting-wear the accuracy of the dissipated energy is the most important: it conditions the evolution of wear profiles.

The dissipated power $P_d$ is obtained by:

$$ P_d(\tau, N_h) = \sum_{j=0}^{N_h} (\tilde{\lambda}^c \cos(j\tau) + \tilde{\lambda}^s \sin(j\tau)) \times \sum_{j=0}^{N_h} (\dot{\tilde{U}}_r^c \cos(j\tau) + \dot{\tilde{U}}_r^s \sin(j\tau)). \quad (37) $$

The error estimator based on the dissipated energy is then defined as:

$$ \epsilon_{E_d}(N_h) = \frac{\int_0^T |\lambda(\tau) \dot{U}_r(\tau) - P_d(\tau, N_h)| d\tau}{\int_0^T |\lambda(\tau) \dot{U}_r(\tau)| d\tau} \quad (38) $$

## NUMERICAL EXAMPLE

The algorithm developed for the coupled calculation of both the wear kinetics and the vibratory response has been applied to a compressor's bladed disk. The fretting-wear zones are the interfaces between disk and blades in the dovetail attachments. The disk (Fig. 3a) has 47 sectors. Its natural frequencies (normalized) are shown in Fig. 3b. All the parts are made of titanium.

It is also assumed that the wear process respects the cyclic symmetry. Therefore, only one sector is studied. This hypothesis is common for bladed-disk studies in the presence of contact nonlinearities. It hasn't been demonstrated but Petrov [21] shows numerical results for which contact nonlinearities don't cause localization phenomena. Vakakis [22] has encountered cases of localizations which break the cyclic symmetry. Indeed, he has revealed a case of mistuning due to nonlinearities between blades. However, here, nonlinearities occur inside each sector, hence the hypothesis.

The nonlinear interaction occurs at the nodes shown in Fig. 4. The friction coefficient is $\mu = 0.5$, corresponding to a contact between parts made of titanium without coating. The wear rate $K_w$ is chosen equal to $1.110^{-11} Pa^{-1}$ as in [14]. The first mode for 0 diameter (a flexural mode) is excited through an harmonic force applied perpendicularly to the top of the blade

<figure><img src="image_4.png" /><figcaption>Figure 3: Analyzed bladed disk: finite element geometry (a) and natural frequencies (b)</figcaption></figure>

(at the level of the isolated red point in Fig. 4). Its amplitude is $0.5N$.

From a purely numerical point of view, 3 harmonics have been used, with 128 time-step in the DFT procedures. Calculations have lasted 1.630s with the improved Jacobian evaluation, against 16.000s with the former approach.

A preliminary calculation is necessary to evaluate the static forces and displacements due to the centrifugal loading. This nonlinear problem can be solved either by a commercial finite element software or by a quasi-static formulation of the algorithm presented in this paper. The second option has been chosen here. The obtained contact pressure distribution is shown in Fig. 5. Without any consideration of wear, the frequency response of the blade displacement amplitude at the level of the excited point is drawn in Fig. 6. For comparison's sake, the vibratory level if the contact were glued (the problem is then linear) is shown besides the solution of the non-linear problem: this is a good example of the damping effect of friction. The average power dissipated during a single cycle in the contact zone is shown in Fig. 7.
<figure><img src="image_5.png" /><figcaption>Figure 4: Calculated bladed sector: finite element geometry (a) and contact nodes (b) on the intrados (top) and on the extrados (bottom). The left side of the grid is at the leading edge of the blade and the right side is at the trailing edge.</figcaption></figure>

For the first fretting cycle, the wear rate profile for each side of the dovetail is shown in Fig. 8. The values are normalized by the maximum value of wear rate in both contact areas. In Fig. 9 are drawn wear depth evolutions for different contact nodes. Wear kinetics is very complex and, of course, very depending the observed location. Nevertheless it is very interesting to remark, on this example, a tendency to reach a steady state, characterized by a null wear rate. This is more evident in Fig. 10. For some nodes, after an initial acceleration, wear rates reach a peak, decrease and then tends to zero. For others, wear rates only decrease.

The final worn geometry for both sides of the dovetail are shown in Fig. 11. These profiles are very asymmetrical.

To go further, the proposed strategy is used to simulate the interaction between wear and vibration in the vicinity of the resonance frequency. Indeed, a very important information is the evolution of the resonance amplitude due to wear (Fig. 12 and Fig. 13). When going on to excite the structure precisely at the initial resonance frequency, wear evolution provokes a decrease of the amplitude. It is a positive aspect of wear for this frequency. But actually the resonance peak has slid towards lower frequencies. And there, the amplitude is higher than at the unworn state. This can be very dangerous if the excitation frequency is not perfectly known and controlled and above all if the structure has not

<figure><img src="image_6.png" /><figcaption>Figure 5: Pressure distribution: on the extrados (a) and on the intrados (b)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Frequency response of the top of the blade</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Average power dissipated per cycle</figcaption></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: Wear rate profile: (a) in extrados (b) in intrados</figcaption></figure>

been designed to endure the corresponding loads. For this example, the results show clearly that wear evolution is coupled with dynamics phenomena. Results found in this study are similar to the results found for a two-degrees-of-freedom fretting-wear model proposed in [7].

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 9: Nodal wear depth evolution at the resonance on the extrados (a) and on the intrados (b)</figcaption></figure>

## CONCLUSIONS AND PROSPECTS

A method has been proposed to calculate simultaneously the wear and the vibrations of bladed disks. Steady-state periodic vibrations are analysed in the frequency domain using a multi-harmonic representation of the structure's displacements. Wear evolution is studied on a slow time scale. This method allows the calculation of forced responses as functions of the wear depth. Conversely it enables one to predict the worn geometry for a given number of cycles, in the presence of inertial forces.

During the calculations, a Jacobian matrix must be evaluated to solve the nonlinear equations of motion and to trace the resonance characteristics as functions of wear depth. Its former evaluation process was time consuming. Here, quasi-analytical expressions have been derived for this matrix to accelerate the calculation.

An efficient jump cycle strategy has been tested to decrease the number of iterations in the wear evolution loop. It is based on a control of the prediction-correction procedure through a constant wear depth target increment.

A new way to calculate the error committed by the Fourier series truncation has been proposed.
<figure><img src="image_13.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b)</figcaption></figure>

Figure 10: Nodal wear rate evolution at the resonance: on the extrados (a) and on the intrados (b)

The new approach has been applied to the finite element model of an assembled bladed disk with dovetail joints. This example has revealed a coupling between dynamics and wear. Wear depths are very small (a few microns) but they modify a lot the vibratory behaviour of the disk. That is why it is strongly recommended that designers take this parameter into account when they dimension dovetail attachments. They can use the method presented here, nevertheless, calibration tests must be performed to guarantee that the life expectancy calculations are faithful to reality.

In the future, the considerations about asymptotic states, that have emerged in this paper, may give rise to another method that would be especially valuable: it would enable one to know rapidly an asymptotic worn geometry, if it exists. Even if the calculation is so direct that it doesn't give the history of the wear process, such a method would take its place in the design loops of optimal dovetails. In the global design process, it would come before the method developed in this paper and work in synergy with it: indeed, the complete dynamical history on the vicinity of a target frequency would only be calculated on the *a priori* "best" candidates - from the point of view of their ultimate profile.

<figure><img src="image_15.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b)</figcaption></figure>

Figure 11: Worn geometry profile: (a) in extrados (b) in intrados

<figure><img src="image_17.png" /></figure>

Figure 12: Evolution of the frequency response of the top of the blade during the fretting cycles

## ACKNOWLEDGMENT

Thanks go to Snecma for its technical and financial support. This work takes place in the framework of the MAIA mechanical research and technology program sponsored by CNRS, ONERA and SAFRAN Group.
<figure><img src="image_18.png" /><figcaption>Figure 13: Excerpt of Fig. 12 for three given frequencies</figcaption></figure>

References

[1] Gallego, L., Nélias, D., and Jacq, C., 2006. “A Comprehensive Method to Predict Wear and to Define the Optimum Geometry of Fretting Surfaces”. *Journal of Tribology*, **128**, pp. 476–485.

[2] Paulin, C., Fouvry, S., and Deyber, S., 2005. “Wear kinetics of Ti-6Al-4V under constant and variable fretting sliding conditions”. *Wear*, **259**(1-6), pp. 292–299.

[3] Charleux, D., Gibert, C., Thouverez, F., and Dupeux, J., 2006. “Numerical and Experimental Study of Friction Damping in Blade Attachments of Rotating Bladed Disks”. *International Journal of Rotating Machinery*, **2006**, pp. 1–13.

[4] Nacivet, S., Pierre, C., Thouverez, F., and Jézéquel, L., 2003. “A dynamic Lagrangian frequency–time method for the vibration of dry-friction-damped systems”. *Journal of Sound and Vibration*, **265**(1), pp. 201–219.

[5] Levy, G., 1980. “Modeling of Coulomb Damping and Wear of Vibrating Systems”. *Wear*, **64**(1), pp. 57–82.

[6] Sextro, W., 2002. *Dynamical Contact Problems With Friction: Models, Methods, Experiments, and Applications*. Springer.

[7] Salles, L., Blanc, L., Thouverez, F., and Gouskov, A., 2008. “Dynamic Analysis of Fretting-Wear in Friction Contact Interfaces”. In *Proceedings of ASME Turbo Expo, 2008*.

[8] Laxalde, D., Salles, L., Blanc, L., and Thouverez, F., 2008. “Non-Linear Modal Analysis for Bladed Disks with Friction Contact Interfaces”. In *Proceedings of ASME Turbo Expo, 2008*.

[9] Laxalde, D., and Thouverez, F., 2009. “Complex Non-Linear Modal Analysis for Mechanical Systems: Application to Turbomachinery Bladings With Friction Interfaces”. *Journal of Sound and Vibration*, p. In press.

[10] Petrov, E. P., and Ewins, D. J., 2003. “Analytical Formula-
tion of Friction Interface Elements for Analysis of Nonlin-
ear Multi-Harmonic Vibrations of Bladed Disks”. Journal

of Turbomachinery, **125**, pp. 364–371.

[11] Borrajo, J. M., Zucca, S., and Gola, M. M., 2006. “Analytical formulation of the Jacobian matrix for non-linear calculation of the forced response of turbine blade assemblies with wedge friction dampers”. *International Journal of Non-Linear Mechanics*, **41**(10), pp. 1118–1127.

[12] Archard, J., 1953. "Contact and rubbing of flat surfaces". *Journal of Applied Physics*, **24**(8), pp. 981–988.

[13] Fouvry, S., Liskiewicz, T., Kapsa, P., Hannel, S., and Sauger, E., 2003. “An energy description of wear mechanisms and its applications to oscillating sliding contacts”. *Wear*, **255**(1-6), pp. 287–298.

[14] Stromberg, N., 1997. "An augmented Lagrangian method for fretting problems". *European Journal of Mechanics. A. Solids*, **16**(4), pp. 573–593.

[15] McColl, I. R., Ding, J., and Leen, S. B., 2004. “Finite element simulation and experimental validation of fretting wear”. *Wear*, **256**(11), pp. 1114–1127.

[16] Dick, T., and Cailletaud, G., 2006. “Analytic and FE based estimations of the coefficient of friction of composite surfaces”. *Wear*, **260**(11-12), pp. 1305–1316.

[17] Meirovitch, L., 2004. *Methods of Analytical Dynamics*. Dover Publications.

[18] Petrov, E. P., and Ewins, D. J., 2002. “Models of friction damping with variable normal load for time-domain analysis of vibrations”. *Proceedings of the International Conference on Noise and Vibration Engineering (ISMA), Leuven*.

[19] Laxalde, D., Thouverez, F., Sinou, J. J., and Lombard, J. P., 2007. “Qualitative analysis of forced response of blisks with friction ring dampers”. *European Journal of Mechanics/A Solids*, **26**(4), pp. 676–687.

[20] Laxalde, D., 2007. “Study of non-linear dampers applied to bladed disks and multi-stage systems”. PhD thesis, École Centrale de Lyon.

[21] Petrov, E. P., 2004. “A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks”. *Journal of Turbomachinery*, **126**, pp. 175–183.

[22] Vakakis, A. F., 1992. “Dynamics of a nonlinear periodic structure with cyclic symmetry”. *Acta Mechanica*, **95**(1), pp. 197–226.