DYNAMIC ANALYSIS OF A BLADED DISK WITH FRICTION AND FRETTING-WEAR IN BLADE

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 loading is proposed. The method is based on the Dynamical Lagrangian 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 performed to calculate the non-linear forces. Then, a Hybrid Powell’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. tact interfaces illustrate the performances of this method and show the coupling between dynamical and tribological phenomena.


INTRODUCTION
One of the most important sources of nonlinearities in assembled 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-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 nonlinear 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. Figure 1: Description of the problem with two solids

Continuous formulation
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 Ω l with smooth boundaries ∂Ω l as drawn in Fig. 1. Each boundary ∂Ω l is divided into three disjoint parts Γ l t , Γ l u and Γ c (common to both solids). Traction forces T l are imposed on Γ l t , displacements on Γ l u and frictional contact with fretting-wear conditions on Γ c .
The constitutive law for the interface is given by Signorini's unilateral contact conditions and Coulomb's law of friction, with µ 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 = 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)(2)(3)(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 u for each point x of both solids (with u N and u T its normal and tangential components) and the contact pressure field p (with p N and p T its normal and tangential components). W designates the dual of w i.e. the driving force for wear. It has the same physical meaning as p N . v, p and W designate the test-fields associated respectively with u, p and W . with is a closed convex set describing the friction and wear limit criterion. It is such that In (2), g is the initial gap between both solids. In (3), the components of each vector are expressed in an orthonormal basis n α , perpendicular to n c .ẇ is obtained through the following equation:ẇ 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: where M, C and K respectively designate the mass, damping and stiffness matrices. For the sake of simplicity, exponent l will be dropped afterwards. F F F ex is the vector of external forces it provides here a periodic excitation. F F F c c 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:

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 U U can be expressed as: where τ is the variable of the fast time scale and η the one of the slow time scale. In the frequency domain, after condensation on the interface nodes, the equation of motion becomes: or whereŨ U U r r r ,F F F r r r and Z r designate respectively the reduced multiharmonic relative displacement vector, reduced stiffness matrix and reduced external excitation.λ λ λ 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λ λ λ. Unfortunately it is not possible to calculate it directly in the Galerkin procedure, indeed it depends on the state of each contact nodestick, 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 theses 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λ λ λ is formulated as a pseudo-penalization of the equations of motion on the tangential and normal directions: ε T and ε N are pseudo-penalty coefficients,X X X r r r is a new vector of relative displacements, which is computed in the time domain. The pair λ λ λ,X X X r r r is determinated through an AFT proce-dure. Equation (11) becomes: The convergence ensures that the time domainX X X r r r and frequency domainŨ U U r r 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 2: Computation of the Lagrangian vector with wear 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: The period is split into N time steps.λ λ λ,λ λ λ u u u andλ λ λ x x x have λ λ λ n n n n=1..N , λ λ λ n n n u u u n=1..N and λ λ λ n n n x x x 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: λ λ λ n,T x x x = λ λ λ n−1,T x x x and λ λ λ n,N x x x = 0. The predicted contact forces are: The corrected contact forces will be: λ λ λ n n n = λ λ λ n n n u u u − λ λ λ n n n x x x , and λ λ λ n n n x x x will we be calculated to satisfy the contact and friction laws.
1. Separation: λ n,N pre ≥ 0 The contact is lost and the forces should be zero.
λ λ λ n n n x x x = λ λ λ n n n u u u .
2. Stick: λ n,N pre < 0 and λ λ λ n,T pre < µ λ n,N pre In this case, the prediction verifies the contact conditions: 3. Slip: λ n,N pre < 0 and λ λ λ n,T pre ≥ µ λ n,N pre 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: 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 δW M 1 , by wear rate integration: 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.

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 ∆N is: If l denotes the considered step of the external loop -still performed on the slow time scale -and η l the time of the beginning of this step, the vector of nodal wears at the beginning of the next step, ∆N(η) cycles later, is obtained by: If ∆N(η) 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 it's 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 re-spectivelyŨ U U r and W are both considered with the introduction of the corresponding Jacobians J and J w .
The Jacobian with respect toŨ U U r is defined as follows: f consists of two parts, a linear one and a nonlinear one: First,λ λ λ 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: where m is the number of the different contact states. τ = τ k and τ = τ k+1 are the transition moments. Using Leibniz's rule, the derivatives of these functions with respect toŨ i are, for the cosine terms (the expressions of derivatives ofλ Since λ x is a continuous and periodic function, proportional to displacement X r , a simplification occurs: (27) And finally the gradients of (25) are: In time domain, the gradient of λ λ λ x x x is calculated according to the contact states defined by Eqs. (17)-(19). The three components of vector λ λ λ n x are λ N,n x in the normal direction and λ λ λ T,n x in the tangential direction (with λ x,n x and λ y,n x its components in an orthonormal basis of the tangential contact plane, perpendicular to n c ).
where I designates the identity matrix, with the same size as Z r .
Finally, the introduction of this expression in (24), gives the Jacobian of function f : The gradient of f with respect to W , denoted J w is: where 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Ũ U U l+1 r r ,W l+1 and Ũ U U l r ,W l are supposed to be solutions of the nonlinear problem respectively at times η l+1 and η l . Then δW W W being known from the wear rate integration, the expression for δŨ U U r is: The prediction-correction algorithm is then used, with U l+1 r 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: 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: The error estimator based on the dissipated energy is then defined as:

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 µ = 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 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. 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  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].

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 multiharmonic 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. 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.