Friction-induced traveling wave coupling in tuned bladed-disks

Flutter is a major constraint on modern turbomachines; as the designs move toward more slender, thinner, and loaded blades, they become more prone to experience high cycle fatigue problems. Dry friction, present at the root attachment for cantilever configurations, is one of the main sources of energy dissipation. It saturates the flutter vibration amplitude growth, producing a limit cycle oscillation whose amplitude depends on the balance between the energy injected and dissipated by the system. Both phenomena, flutter and friction, typically produce a small correction of the purely elastic response of the structure. A large number of elastic cycles is required to notice their effect, which appears as a slow modulation of the oscillation amplitude. Furthermore, even longer time scales appear when multiple traveling waves are aerodynamically unstable and exhibit similar growth rates. All these slow scales make the system time integration very stiff and CPU expensive, bringing some doubts about whether the final solutions are properly converged. In order to avoid these uncertainties, a numerical continuation procedure is applied to analyze the solutions that set in, their traveling wave content, their bifurcations and their stability. The system is modeled using an asymptotic reduced order model and the continuation results are validated against direct time integrations. New final states with multiple traveling wave content are found and analyzed. These solutions have not been obtained before for the case of microslip friction at the blade attachment; only solutions consisting of a single traveling wave have been reported in previous works.


Introduction
Flutter onset is currently a very important limitation in modern turbomachinery, where blade designs are getting more slender and closer to their mechanical limit. It is an aeroelastic instability where the gas flowing around the blades tends to amplify the small elastic oscillations of the blades. As a consequence, the blade vibration amplitude grows exponentially until nonlinear effects, typically due to friction forces at the blade root in the case of Low Pressure Turbines (LPT), become relevant. Hopefully, the dissipation produced by nonlinear friction is strong enough to limit the growth of the aeroelastic instability and avoid a catastrophic blade failure. In this case, there is a balance between nonlinear friction damping and unstable flutter growth such that a final Limit Cycle Oscillation (LCO) state sets in with a bounded blade vibration amplitude.
The determination of this final vibration amplitude is crucial to estimate the blade fatigue level. This typ-ically requires a CFD description of the aerodynamic flow coupled with a FEM for the motion of the bladeddisk with nonlinear friction effects. These simulations can include a number of simplifications (linearized CFD, partial fluid-structure coupling, reduced number of modes for the motion of the structure,...), but they are always quite involved and CPU costly.
In a tuned configuration, if the final state is made of a single Traveling Wave (TW), then the calculations can be reduced to just one flow passage and one sector of the bladed-disk. This constitutes a huge computational reduction in the usual case of large blade counts.
For this reason, there has recently been some research activity looking into the expected TW composition of the final LCO state. Aeroelastic models of the bladed-disk with different levels of detail have been used to try to decide whether, in a tuned rotor, the final LCO state is always made of a single TW or if there can also be multi-TW final states.
In [4], a mass-spring model with microslip friction is used to analyze the problem. The system is integrated in time with fully coupled and semi-uncoupled linearized aerodynamic effects. Both formulations produce very similar results, indicating that, in the usual case of small vibration amplitude, aerodynamic forces can be safely assumed to be linear. Time evolution of the system first showed the slow growth of the unstable TWs and, afterward, the nonlinear interaction of unstable TWs through the effect of friction at the firtree, which takes place on an even longer time scale. In all cases reported, this interaction always ended up selecting a final state composed of just one single TW, the most aeroelastically unstable one.
A multiple scales method is applied in [9,10] to further simplify the mass-spring model. The resulting asymptotically reduced equations describe the evolution of the system directly in the slow time scale (associated with the small effects of flutter and friction), with the fast elastic oscillations filtered out and can be integrated with a much lower computational cost. Again, the simulations performed always produced single TW final states, although other unstable TWs close to the most unstable one were also found as possible final states of the system. The selection of the final single TW state is sensitive to the initial conditions. These initial conditions were picked from a random distribution [9] or both from a random distribution and starting from a pure TW [10].
Multi-TW solutions were reported as another possible final state in [8]. The model used consisted of 7 blades with linear aerodynamic effects and Coulomb friction forces to describe the effect of the contact between adjacent blades. In order to obtain the multi-TW states, the linear aeroelastic characteristics of the system were selected to have two unstable TWs with a quite high negative aerodynamic damping value. A realistic model of a LPT bladed-disk with linear aerodynamic forces and contact interaction at the tip shroud was considered in [7]. Here, it was shown that the inclusion of strong nonlinear coupling effects at the tip shroud can produce multi-TW final states. The analysis was continued in [2], to show that nonlinear tip contact effects change the vibration modes and frequencies of the system with respect to those obtained with linearized contact conditions. These changes in frequencies and mode shapes have to be considered in order to correctly estimate the amplitude of the resulting limit cycle oscillations.
In the present work, we consider the case of flutter vibrations with microslip contact forces at the blade fir-tree, where the effect of friction does not substantially change the natural vibration modes. As it was mentioned above, previous works on this configuration [4,9,10] have only found final vibration states consisting of just one single TW. The idea is to apply numerical continuation techniques [6] to analyze the stability of the solutions and search for multi-TW states. This procedure eliminates the need for direct time integration of the system, which requires very long integration times and always leaves some uncertainty about whether the final state is fully converged.
To this end, first, we briefly introduce the asymptotic model used, study the stability of the zero solution (no blade vibration) and compute the single TW solutions that appear when the magnitude of the flutter instability is increased. Then, we analyze the linear stability of single TW solutions and perform numerical continuation to follow new branches of multi-TW solutions that bifurcate from single TW solutions when they lose stability. New multi-TW states are obtained, with very rich TW content. They have two different frequencies and propagate around the rotor with nonuniform blade to blade amplitude.

Asymptotic model
A simple mass-spring model is used to represent the bladed-disk (Fig. 1). The system has 2 degrees of freedom per sector, one for the blade motion x j and one for the friction displacement y j . It includes: elastic coupling among the different blades to account for the effect of the disk K c , linear aerodynamic forces F a , and microslip friction forces F f (y) at the blade attachment. Similar models have been successfully used in previous works for the description of flutter vibration of a nearly flat modal family (see [9,10]).
The equations of motion of the mass-spring system are given by where m b is the mass of the blade and m f is the mass of the nonlinear friction Degree of Freedom (DOF) in each sector. The elastic coupling between both DOF is k b . It is convenient to rewrite this set of equations in dimensionless variables in order to clearly compare the relative importance of the different terms. Dividing Eq. (1) by the characteristic elastic force of the system wherex = x/(F c /k b ),ỹ = y/y c (y c is the microslip characteristic displacement), the time is made dimensionless T = ω b t with the blade alone natural fre- The nondimensional parameter θ is usually very small because the displacement of the fir-tree is typically much smaller than the modal displacement of the blade (or, in other words, the stiffness of the friction is much larger than that of the blade elastic mode).
There are two very different time scales in this problem. One corresponds to oscillations with the elastic frequency of the blade and, the other, which is much slower, is associated with the small effects of flutter and nonlinear friction. A multiple scales method (see, e.g., [1]) can be applied to obtain, in the limit θ 1, an asymptotic model with the fast elastic oscillations filtered out. The derivation of the asymptotic model is completely analogous to that explained in detail in [9,10] and is omitted here.
The asymptotic model expressed in terms of complex TW amplitudes A j and the slow time scale τ = θ T , can be written in the form where N = 24 is the number of blades of the rotor. The elastic and aerodynamic matrices are diagonal in the TW basis and the friction matrix is diagonal in the displacement basis E is the change matrix from TW amplitudes to blade displacements, and E H denotes the conjugate transpose of E. This matrix represents the discrete Fourier transform and is given by The complex blade displacements are represented by X j , and they are related to the dimensionless displacementsx j bỹ which corresponds to the first order in the asymptotic expansion [9,10]. Furthermore, the complex displacements X j are related to the TW amplitudes A j by X = EA.
The coefficients ω j in the elastic matrix represent the small deviations of the natural frequencies of the modal family from the blade alone normalized frequency (see Fig. 2). The material damping coefficient of the structure has the value ξ mat = c mat /2 √ k b m b = 2 × 10 −4 (much smaller than the aerodynamic damping).
The Olofsson model [11] is used to describe the microslip friction at the blade attachment. The dimensionless friction forceF f j is related to the DOFỹ j by a loading-unloading cycle, which is shown in Fig. 3. The magnitude of the friction force corresponding to the unloading part of the cycle is given bỹ and, the subsequent loading part is obtained as The complex friction coefficient Q(|X j |) in the friction matrix M friction used in the asymptotic model is related to the first term of the Fourier series expansion ofỹ j (see [9]). This coefficient accounts for the nonlinear dissipation and frequency change produced by friction at the fir-tree and depends on the amplitude of the blade elastic cycle |X j |. Real and imaginary parts of Q(|X j |) are shown in Fig. 4. The microslip regime is bounded by |X | ≤ 0.5.
The coefficients in M aero are the aerodynamic damping and frequency correction of the TW modes. A sinusoidal shape (blade aerodynamic interaction with itself and with its adjacent blades only) has been selected to  mimic the typical data in a realistic LPT (see, e.g., [3]) where j is the TW wavenumber, ξ a 0 = 0.25, ξ a 1 = 0.75, ξ a 2 = −0.05 , η a 0 = 1.4 and η a 1 = 1. The aerodynamic coefficients are shown in Fig. 5. The aerodynamically unstable TW modes (i.e., ξ k a < 0) are represented with a shaded area. Note that the mean aerodynamic damping is stable and there are 9 unstable TWs. Also, as it is usual in the turbomachinery field we also use the term number of Nodal Diameters (ND) to refer to the wavenumber of a TW.
The coefficient ξ a 1 will be used later on as the bifurcation parameter to explore the solutions of the system. By changing this parameter, we explore different levels of the flutter instability increasing or decreasing the number of unstable TW modes and their damping.

Traveling wave solutions
Equation (3) admits the trivial solution, which corresponds to no blade motion. Linearizing the problem around this solution A j = 0 + a j yields d dτ where a j represents the perturbed TW amplitude, with θ . This is a diagonal constant coefficient system, and the eigenvalues of the matrix are given by As the flutter intensity is increased through the parameter ξ a 1 [see Eq. (11)], the trivial solution becomes unstable for the first time atξ j = ξ j a +ξ mat θ = 0 for any j. As it can be seen in Fig. 5, this happens first for j = 6, which is the most aeroelastically unstable TW mode. For higher values of ξ a 1 , the adjacent TW modes become also unstable. The effect of friction is just a shift in frequency of magnitude R[Q(0)]/2 in the aeroelastic modes.
There is another type of simple solutions of the system which corresponds to a single TW with wavenumber r , with all blades vibrating with the same amplitude and frequency where √ N R r is the modulus of the TW amplitude, m r is the frequency correction of the single TW solution and α represents an arbitrary phase. Note that the magnitude of the blade displacements |X j | = R r for every j, so all the entries in the friction matrix in Eq. (3) are equal to Q(R r ). Therefore, after introducing the expression (15) in system (3), the following nonlinear equations are obtained which determines the value of R r for a given aerodynamic levelξ r , and which directly gives the frequency of this periodic solution once the value R r is known. Equation (16) indicates that the amplitude of the single TW periodic solution is the result of the balance between the growth rate of the aeroelastic instability and the nonlinear damping of the friction [10].
Since J[Q(|X |)] ≤ 0, the existence of single TW solutions requiresξ r < 0. Therefore, for each unstable TW mode with wavenumber j = r, there is a single TW solution with complex amplitude A r that bifurcates from the trivial solution precisely when the real part of the eigenvalue λ r from (14) becomes positive. The stability of these single TW solutions is analyzed in the next section

Stability of the traveling wave solutions
Single TW solutions [Eq. (15)] are periodic solutions of system (3). To study their stability, we introduce a change of variable to make these solutions steady in the new formulation. By factoring out the frequency of single TW solutions, m r , we have Introducing the change of variable (18) in Eq. (3), and substituting the value of m r from (17), gives d dτ where single TW solutions are now steady solutions of this autonomous system. The stability characteristics of single TW solutions can now be obtained by computing the eigenvalues of the Jacobian matrix of Eq. (19), evaluated at one of these solutions. To this end, we linearize the system writing the solution as a single TW with wavenumber r plus a small perturbation ⎛ Expression (20) is introduced in the system of eqs. (19) and only linear terms in the perturbation are retained. The resulting linear system for the time evolution of the TW perturbations can be expressed as d dτ where Q (R r ) is the derivative of the complex friction coefficient Q(|X |) evaluated at R r . System (21) has some particular characteristics that are worth noting.
The complex conjugate of the perturbation,ā j , appears through the linearization of the modulus of the blade displacement |X j |, present in the friction terms Q(|X j |). A perturbation with wavenumber r + k on a TW with wavenumber r produces also a complex conjugate component with wavenumber r − k. As a consequence, despite the fact that the system is tuned (that is, all sectors of the bladed-disk are identical), different TW perturbations are not uncoupled. Because of the nonlinearity of the friction terms, the TW perturbations a r +k and a r −k are coupled in pairs, and the perturbation a r is decoupled from the rest.
The second term in system (21) contains the damping differencesξ j −ξ r . Near the most aeroelastically unstable TW, the aerodynamic damping values are very similar (see Fig. 5). Their differences are thus very small and generate even longer time scales, associated with the selection of a particular TW as the final state.
This long time scale has been noticed in previous works [4,9,10], and it can be clearly seen in Fig. 6, where system (3) is integrated in time starting from a small initial condition (A j = 0.01 for every j) and with ξ a 1 = 0.75. Initially, there is an exponential growth of the aeroelastically unstable TW modes (straight lines due to logarithmic scale). The slope of this growth (blued dashed line) is practically given by the most unstable aerodynamic damping ξ 6 a (the material damping of the structure is much smaller). However, the slope of the TW decay in the nonlinear interaction region (red dashed line), which is related to the difference of aerodynamic dampings ξ j a − ξ r a , is much smaller and it takes a very long time to select the final TW with wavenumber 6. Direct time integration, even for this asymptotic model, requires long computation times, and there is always a certain level of doubt about whether the convergence to the final state has been reached. Now, expanding the small perturbation a = u + iv, Eq. (21) can be written in real and imaginary part as d dτ where J(R r ) is the Jacobian matrix evaluated at the amplitude R r of a single TW solution with wavenumber r . Therefore, a single TW solution branch is stable if the real part of all the eigenvalues of the Jacobian matrix are less or equal to zero or unstable if the real part of any of the eigenvalues is positive. As an example, for wavenumber 8, the single TW solution curve is computed for values of the bifurcation parameter ξ a 1 ∈ [0.2, 1.0]. Its stability is determined by the eigenvalues of the Jacobian matrix at each point, as shown in Fig. 7. For values close to ξ a 1 = 0.2, only the trivial solution exists, which becomes unstable at the pointξ 6 = 0, since the TW with wavenumber 6 is the most unstable. When the bifurcation parameter is increased andξ 8 = 0, the single TW8 branch emerges from the trivial solution, and it is initially unstable. Then, after a Hopf Bifurcation (HB) point, the LCO corresponding to TW8 becomes a stable solution.
The same process is repeated again for every TW mode that becomes aerodynamically unstable in the interval ξ a 1 ∈ [0.2, 1.0]. Gathering the information of each individual bifurcation diagram (as in Fig. 7), the stability of every single TW solution and the trivial solution is represented in Fig. 8. The trivial solution corresponds to the horizontal plane and is stable when there is no aerodynamically unstable TW mode present in the system. The straight dotted line corresponds tõ ξ 6 = 0, where the trivial solution becomes unstable. The single TW solutions emerge from the trivial solution when a TW mode becomes unstable (This condition is represented by the curved dotted line and is given byξ r = 0).
The single TW solution with wavenumber 6 (which corresponds to the most unstable aerodynamic damping) is always stable, but adjacent ones (with slightly less aerodynamic damping) only become stable as the bifurcation parameter ξ a 1 is increased. Every single TW that changes stability does it through a HB. Also, for a given value of the bifurcation parameter ξ a 1 , even though multiple single TW solutions exist, not every one of them has a stable LCO.
Previous works only reported the most aerodynamically unstable TWs as final states of single TW solutions [10] after performing numerical integration. Also, the most unstable mode seemed to have the largest basin of attraction. In this work, the stability results show that, for example, when ξ a 1 = 1, only 5 out of the 10 aerodynamically unstable TW modes have a stable LCO. These TW modes are precisely the most aerodynamically unstable ones, and their stability region is larger the more unstable a TW mode is.

Multi-traveling wave solutions
We now look for more complicated stable solutions that can bifurcate from the HB in single TW solutions. The approach to track these branches and compute their stability is to use numerical continuation. In particular, the numerical continuation software AUTO [6] is used. AUTO is a package that can continue steady solutions of ordinary differential equations and determine its bifurcation points using one (or more) free parameters. Furthermore, it continues periodic solutions, determining the Floquet multipliers at each step to compute the stability of these branches. The algorithms implemented in this software are explained in detail in [5].
Even though single TW solutions are periodic, recall that they became steady in Eq. (19) after performing a The change in stability of single TW solutions occurred through a HB, so we apply numerical continuation to track the branches that emerge from these points. The results show that a stable state only emerges from TW4. From the other single TW solutions, the new branch was always unstable.
The new solutions are represented in Fig. 9. Here, from the single TW with wavenumber 4 a multi-TW solution appears from the HB point. This multi-TW solution is stable for a very narrow range of values of the bifurcation parameter (ξ a 1 ∼ 0.5 − 0.515), but its TW components grow to a considerable amplitude (see the inset in Fig. 9 where the TW modes with largest amplitudes are represented). As one continues decreasing the bifurcation parameter, a Torus Bifurcation (TB) occurs and the multi-TW solution loses stability. In addition, the new frequency associated with this multi-TW branch is represented in Fig. 10. This frequency grows along the multi-TW branch and coexist with the single TW frequency m 4 .
The TW content of the stable multi-TW solution found is illustrated in Fig. 11 for different values of the bifurcation parameter ξ a 1 . As the bifurcation parameter moves closer to the end of their stability region, there is more modal content of different TWs, except for the  Numerical time integration of system (3) has also been performed in order to check the results of the numerical continuation. The time evolution of the different TW components of the solution is presented in Fig. 12 for ξ a 1 = 0.502, showing the convergence of the system to a multi-TW state. The initial condition for this simulation is the single TW4 solution plus a small perturbation of 0.01 to every TW mode. Nevertheless, for this value of the aerodynamic damping there are 4 stable single TW LCOs that coexist with the multi-TW solution, so for different initial conditions the final state can be a single TW. Also, notice that, in contrast to the case where the final state was a single TW solution (Fig.  6), the required time integration for multi-TW solutions is now roughly an order of magnitude larger.
Finally, the displacements of the blades |X j | of the new multi-TW solution for ξ a 1 = 0.502 are plotted   Fig. 13. The multi-TW solution takes the form of a traveling wave with non-uniform amplitude that rotates around the bladed-disk. The solution represented in this figure is actually the envelope of the fast elastic oscillation (which was filtered in the asymptotic model) of blade j. The period of this envelope (roughly ∼ 300) is the one associated with the new frequency present in the multi-TW solution (see Fig. 10).
Despite the complexity of the multi-TW solutions that are obtained in this work, the maximum blade displacement that they produce is always smaller than that associated with the single TW solution with highest amplitude (corresponding to the most aerodynamically unstable wavenumber 6). As presented in Table 1, the maximum blade displacement of the TW solution is 17% higher than that of the multi-TW solution. There-  In order to have a wider range of existence of stable multi-TW solutions, we have explored a number of modifications of the parameters in system (3) Modifications (i) and (ii) changed the stability regions of single TW solutions. Decreasing the frequency correction of both the structure and the aerodynamic matrix enlarged the single TW instability regions. However, by doubling the frequency terms, the stability region of TW solutions grew with respect to the case shown in Fig. 8. Nevertheless, after performing numerical continuation it was found that there were no stable multi-TW states in neither of these cases. The changes in aerodynamic damping characteristics are represented in Fig. 14. The baseline configuration is the one already used throughout this paper. Two new configurations are explored: (a) Wide, where the damp-  ing coefficients of the unstable region are more similar, and (b) Peak, with more variation in the unstable aerodynamic dampings.
It was found that the only modification that ended up producing more stable multi-TW solutions was the one that made the unstable aerodynamic damping values more similar (Wide in Fig. 14).
For the Wide case, the multi-TW solution bifurcates from the TW solution with wavenumber 2 through a HB and connects again with the same branch at a second HB point (Fig. 15). With respect to the baseline case, this new solution is now stable for a much larger interval of values of the bifurcation parameter ξ a 1 . Furthermore, as before, the maximum displacement of the blades produced by the multi-TW solution remains below that of the single TW solution with highest amplitude (see Table 2), even though their difference is now much smaller.

Conclusions
Multi-TW stable states are found in a simplified model of a bladed-disk with microslip friction at the bladeddisk attachment. The model describes the vibration of a nearly flat modal family with several aerodynamically unstable modes and is derived through the application of a multiple scales method that allows to filter out the fast elastic blade oscillation. Numerical continuation techniques are applied to follow the branches of new multi-TW solutions that bifurcate from single TW solutions at the stability change. Previous works reported only single TW final states, which were obtained through direct time integrations of the model. The continuation method directly computes the final states and their stability, avoiding the need of performing very long time integrations to reach a converged final state.
The following final remarks summarize the main results of this work: 1. The stability analysis of the zero solution (no blade vibration) indicates that nonlinear friction only produces a frequency shift in the aeroelastic modes. 2. Single TW solutions bifurcate from the trivial solution as the flutter intensity is increased. The stability analysis of these solutions shows that there is a long time scale associated with the difference between the aerodynamic damping values of unstable modes, which is responsible for the very slow selection of the final state. This long time scale had been detected in time integrations presented in previous works [4,9,10]. 3. There is always a stable single TW solution, which corresponds to the most unstable aerodynamic mode and has the largest amplitude of all single TW solutions. The adjacent TW solutions are first unstable, and they can become stable as the intensity of flutter is increased. There may be several stable single TW limit cycles for the same parameter values. 4. At the point of stability change of single TW solutions, new multi-TW solutions appear. Multi-TW solutions have two different frequencies, with non-uniform blade vibration amplitude, and rotate around the bladed-disk. The range of stability of these multi-TW solutions is seen to be related to the presence of unstable aeroelastic modes with similar aerodynamic damping values. The closer the unstable aerodynamic damping values, the larger the region of stability of the multi-TW solutions. 5. Multi-TW solutions give rise to blade vibration amplitudes that are always below that of the single TW solution. Therefore, from the point of view of blade fatigue calculation, it seems that it is conservative to consider only the highest amplitude single TW solution.