Alain Batailly

Mathias Legrand

Antoine Millecamps

François Garcin

Conjectural bifurcation analysis of the contact-induced vibratory response of an aircraft engine blade

Keywords: rotor/stator structural interaction, rubbing, nonlinear dynamics, unilateral contact dynamics, bifurcation analysis interaction rotor/stator, frottement aube/carter, dynamique non-linéaire, dynamique du contact, analyse de stabilité

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.

Introduction

In turbomachinery, fan, compressor and turbine blades are highly prone to structural failures [START_REF] Carter | Common failures in gas turbine blades[END_REF]. Such undesirable events shall be induced by either unavoidable high and low cycle fatigue, corrosion or creep mechanisms to name a few. Other damages stemming from structural unilateral contact occurrences between components, for instance, have always been considered as accidental in conventional engine designs. However, modern geometries with tighter operating clearances for improved energy efficiency are such that blade-tip/casing or blade-tip/abradable coating interaction events are becoming commonly accepted in standard operating conditions with possibly dangerous consequences on blades service life. Accordinlgy, rotor/stator structural interactions within turbomachines [START_REF] Morfey | Rotating pressure patterns in ducts: their generation and transmission[END_REF] were targetted in various numerical and experimental investigations related to the aerospace [START_REF] Sinha | Rotordynamic analysis of asymmetric turbofan rotor due to fan blade-loss event with contact-impact rub loads[END_REF][START_REF] Kou | Rub-induced non-linear vibrations of a rotating large deflection plate[END_REF] and electric power industries [START_REF] Torkhani | Light, medium and heavy partial rubs during speed transients of rotating machines: numerical simulation and experimental observation[END_REF]. Intricate and possibly diverging dynamics going from a single blade [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF][START_REF] Padova | Experimental results from controlled blade tip/shroud rubs at engine speed[END_REF] or a full bladed disk assembly [START_REF] Schmiechen | Travelling wave speed coincidence[END_REF] to several compressor stages with precessional shaft motions [START_REF] Sinha | Rotordynamic analysis of asymmetric turbofan rotor due to fan blade-loss event with contact-impact rub loads[END_REF] are found and several partially unelucidated industrial incidents are reported [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF][START_REF] Williams | Simulation of blade casing interaction phenomena in gas turbines resulting from heavy tip rubs using an implicit time marching method[END_REF]. Invariably, the engine angular velocity is identified as a central parameter and the prediction of critical speeds yielding large vibratory amplitudes is of primary importance for the design of blades robust to unilateral contact conditions.

This contribution brings further developments to previous numerical and experimental investigations limited to the framework of a single blade [START_REF] Batailly | Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF][START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF]. Possible contacts between neighbouring blades [START_REF] Krack | Reduced order modeling based on complex nonlinear modal analysis and its application to bladed disks with shroud contact[END_REF][START_REF] Chu | Impact vibration characteristics of a shrouded blade with asymmetric gaps under wake flow excitations[END_REF] are not accounted for. Even though it now seems possible to conduct simulations able to accurately predict critical angular velocities as well as the corresponding abradable coating wear patterns, experimental data sets show unexpected and sudden variations in the blade vibratory responses [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF] that remain unelucidated. A scenario which assumes a brutal and significant alteration of the abradable material properties is investigated in [START_REF] Batailly | Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF] where the blade response does feature sudden variations similar to experimental ones. Nonetheless, this scenario seems unsatisfactory because (1) it implies permanent contact separation after the interaction which is inconsistent with experimental facts and (2) the scale of the proposed alteration is unrealistic.

A first key ingredient of the present work lies in the post-processing strategy. As opposed to the previous solution methods [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF][START_REF] Batailly | Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF] essentially based on Fourier analyses of presumed periodic steady-states, the analysis is here expanded to transient time histories through a modified Empirical Mode Decomposition (EMD) procedure. The EMD is a recent signal processing method [START_REF] Huang | The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF][START_REF] Smith | The local mean decomposition and its application to eeg perception data[END_REF] which is well-suited to non-stationary and nonlinear1 time responses even though its mathematical foundation remains partially open [START_REF] Rilling | On empirical mode decomposition and its algorithms[END_REF]. Time-domain signals are projected onto a basis of Intrinsic Mode Functions (IMF) featuring specific mathematical properties such as completeness and orthogonality. In the field of turbomachinery, EMD-based methods such as the Local Mean Decomposition (LMD) [START_REF] Smith | The local mean decomposition and its application to eeg perception data[END_REF] have been extensively applied for monitoring and diagnosis purposes [START_REF] Yang | An ensemble local means decomposition method and its application to local rub-impact fault diagnosis of the rotor-systems[END_REF] focusing essentially on their applicability to nonlinear complex time responses. The EMD is here implemented with an adjusted procedure for the computation of the IMF based on B-spline interpolation [START_REF] Huang | The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] of time signal local extrema. This allows for a computationally efficient analysis of very long time responses.

A second key ingredient is the in-depth bifurcation analysis of the blade dynamics undergoing structural unilateral contacts. The numerical strategy previously introduced in [START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF] is used over a properly selected angular speed 2 . The numerical simulations carried out are twofold: on one hand, interactions between the blade and the surrounding casing are simulated under vacuum, consistently with the experimental setup detailed in [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF]. On the other hand, identical interactions are simulated by incorporating a simplified aerodynamic loading on the blade. Direct confrontation between the results obtained in both configurations sheds a new light on the possible connection between the vacuum chamber used experimentally and the sudden variations in the stress signals. The combination of unilateral contact forces and periodic external aerodynamic loading on the blade is systematically explored through bifurcation diagrams [START_REF] Virgin | Introduction to Experimental Nonlinear Dynamics[END_REF], phase diagrams and Poincaré maps (which are one snapshot of the phase space per forcing period).

In the first section, the selected interaction scenario [START_REF] Batailly | Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF] is briefly summarized along with the relevant experimental observations [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF]. Unpublished experimental results are also introduced to give a full picture of the phenomena. The EMD post-processing procedure is then detailed. In the third section, the blade model together with assumptions regarding the modeling of the aerodynamic loading are provided. The last two sections display various numerical results with and without aerodynamic loading. Extensive time and frequency domain analyses of the blade response are provided and a plausible scenario explaining the sudden drops of amplitude is suggested. When aerodynamic loading is accounted for, attention is paid to the nature of the simulated interactions with in-depth investigations of the bifurcations.

Note: All frequencies and angular speed are in Hertz. However, for confidentiality purposes, numerical data are normalized and units are omitted for numerical results.

Experimental data and interaction scenario

The experimental setup described in [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF] involves the last stage of a full-scale aircraft engine lowpressure compressor operating under vacuum. One blade is slightly longer to initiate contact with the abradable coating deposited on the casing. It is also instrumented with three strain gauges as depicted in Fig. 1. Contact between the blade-tip and the abradable coating is initiated through the centrifugal load acting on the blade. From t = 0 s to t = 113 s, the compressor rotates at very low speed. At t = 113 s, the compressor angular speed undergoes a sudden increase to reach the targeted velocity Ω exp . Between t = 133 s and t = 195 s the compressor rotates at Ω exp . The corresponding gauge time histories are displayed in Figs. 2a, 2b and 2c. Stress responses are all normalized with respect to the mean stress induced by the static centrifugal load σ ω .

Six phases are distinguished during the experiment: Phase I (from t = 113 s to t = 133 s): the compressor angular speed increases linearly with time. Phase II (from t = 133 s to t = 151 s): even though no significant vibration is detected, it is assumed contact is initiated during this phase [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF]. Phase III (from t = 151 s to t = 156 s): a complex diverging blade displacement is observed, significant increases of amplitudes are captured on each gauge. Phase IV (from t = 156 s to t = 162 s): no more significant dynamical activity is noticed and average amplitudes of phase II are retrieved. Phase V (from t = 162 s to t = 165 s): unstable behaviour similar to the one observed in phase III is witnessed: the amplitude of the stress signals suddendly increases to similar levels as those in phase III. Phase VI (from t = 165 s to t = 192 s3 ): stresses decrease again to levels as low as in phases II and IV before an unstable behaviour is established leading to the structural failure of the blade. Reference [START_REF] Batailly | Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF] describes an attempt to numerically reproduce the aforementioned experimental setup: the interaction angular speed as well as the abradable coating wear pattern and the areas of maximum stresses were accurately predicted with attention paid to phase VI. Recent studies conducted both experimentally [START_REF] Mandard | Interacting force estimation during blade/seal rubs[END_REF][START_REF] Cuny | Experimental measurement of dynamic forces generated during short-duration contacts: application to blade-casing interactions in aircraft engines[END_REF] and numerically [START_REF] Salvat | Modeling of abradable coating removal in aircraft engines through delay differential equations[END_REF] were motivated by the interaction phenomenon reported in [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF]. However, the unexpected jumps in amplitude were not addressed and the sequence from phase II to phase V is currently not understood.

Empirical mode decomposition

Various post-processing strategies are devoted to the analysis of nonlinear time histories stemming from nonlinear mechanical systems, such as wavelets [START_REF] Pratumnopharat | Wavelet transform-based stress-time history editing of horizontal axis wind turbine blades[END_REF][START_REF] Wang | Multi-scale enveloping order spectrogram for rotating machine health diagnosis[END_REF] or spectrograms. One disadvantage of the former lies in the fact that several mother wavelet functions may have to be tested for each nonlinear signal and a weakness of the latter is its sensitivity to the length of the time window. Promising EMD-based approaches [START_REF] Huang | The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] advantageously provide physically relevant results and have already been successfully used in the field of rotordynamics [START_REF] Yang | An ensemble local means decomposition method and its application to local rub-impact fault diagnosis of the rotor-systems[END_REF][START_REF] Peng | An improved hilbert-huang transform and its application in vibration signal analysis[END_REF]. That is the reason why they are considered in this work.

Basic principle

The EMD decomposes a time signal y(t) into a finite sum of N Intrinsic Mode Functions (IMF) denoted φ i (t):

y(t) = N i=1 φ i (t) + r(t) (1)
where r(t) is a residual function from which no IMF can be extracted, such as a monotonic function.

The IMF must satisfy two criteria in order to be valid. As stated in [START_REF] Huang | The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF]:

1. "in the whole data set, the number of extrema and the number of zero crossings must either equal or differ at most by one;" 2. "at any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero."

Mathematically, the IMF should also form a complete and orthogonal adaptative basis. While completeness is guaranteed by Eq. (1), orthogonality is not theoretically proven and must be carefully scrutinized based, for instance, on the criteria detailed in [START_REF] Huang | The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF].

Constructing an IMF φ i (t) typically relies on a double iterative process shown in Fig. 3 which involves (1) the sifting process where the computation of a mean function m k (t) based on the detection of the local extrema of the signal s k (t) (with s 0 (t) = y(t)) is undertaken [START_REF] Smith | The local mean decomposition and its application to eeg perception data[END_REF] and (2) collecting the envelope functions a k (t). The computation of the envelope functions depends on the

y(t) while y(t) = r s k (t) = y(t), while a k (t) = 1 m k (t) = mean of s k (t) h k (t) = s k (t) -m k (t) ⇒ a k (t) s k+1 (t) = h k (t)/a k (t) IMF i (t) = Π k j=1 a j (t)s k+1 (t) k = k + 1 y(t) = y(t) -IMF i (t) i = i + 1 Figure 3.
EMD principle as used in [START_REF] Batailly | Conjectural bifurcation analysis of an aircraft engine blade undergoing 3d unilateral contact constraints[END_REF] context in which the EMD is performed. As an example, their acquisition is associated to a moving average algorithm in the LMD [START_REF] Smith | The local mean decomposition and its application to eeg perception data[END_REF] while it relies on a B-spline interpolation of local extrema in the Hilbert-Huang transform [START_REF] Huang | The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF].

Similarly to the Hilbert-Huang transform [START_REF] Huang | The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF], it is here proposed to perform a B-spline interpolation of the local extrema of the signal while enforcing the symmetry of the two envelopes with respect to y = 0 (the one of the local maxima and the one of the local minima) as featured by the LMD; this is achieved by considering local maxima as well as local minima-if they are negative-for the B-spline interpolation of the envelopes.

One interesting feature of the EMD is its ability to discrimate locally the frequency content of a signal. By nature, the EMD extracts the highest frequency content in the first IMF and lower frequency contents are spread through the other IMF. Note that while the first IMF always refers to the highest frequency content, the frequency range to which it is actually associated is driven by the features (localized noise, sudden jumps, discontinuities. . . ) of the signal and may thus vary with time. In the end, and in contrast with the well-known Fourier transform, the EMD leads to a description of the signal with functions whose amplitude and frequency are time-dependent.

Validation

Nonstationary amplitude and frequency features in the signature of a vibratory response induced by complicated nonlinear dynamics can be efficiently tackled by the EMD. The proposed EMD implementation is succinctly validated with one simple illustrative example. The sampling frequency f is given and it is verified that there is a minimum sampling frequency f m < f for which the results of the EMD remain identical. The signal: y 2 (t) = (10 + t 2 e -0.2t ) sin(3t + t 2 e -0.1t )

(2) is sampled at f = 4 kHz and is plotted in Fig. 4a low-frequency content of y 2 (t) is stored in IMF φ 3 while its high-frequency content4 is collected in φ 1 and φ 2 .

The first two IF γ 1 and γ 2 are sensibly constant and are not shown for brevity. Figure 4e shows that IF γ 3 asymptotically tends to the theoretical limit f = 0.477 Hz when t → ∞. Here again, the sensitivity of the decomposition to edge effects is noticeable with oscillations in γ 3 for t > 90 s. The noisy portion of γ 3 for t < 20 s is induced by the numerical sensitivity of the Hilbert transform when the frequency content changes rapidly.

It is understood that this proposed test signal, yet relatively simple, validate the EMD implementation. Key features such as the accuracy of the sifting process for the extraction of the IMF and the proper detection of instantaneous frequencies is established. The procedure can further be employed for the analysis of more elaborate responses obtained from nonlinear contact simulations.

EMD analysis of experimental data

An EMD of the experimental time histories depicted in Fig. 2c is carried out in this section. The first five non negligible IMF are pictured in Fig. 5 where the six interaction phases mentioned in section 2 can be identified. Key interaction phases III, V and VI with highest amplitudes are mostly portrayed by IMF φ 1 and φ 2 . In particular, the last divergent phase VI is almost essentially described by IMF φ 1 which indicates that the highest frequency content during this phase is dominant.

The Hilbert spectrum of the third gauge is depicted in Fig. 6. It shows the amplitude of every major IMF in the frequency/time plane. The unavoidable experimental noise is first filtered out in order to obtain a comprehensible spectrum : high frequency content corresponding to negligible IMF are not pictured in this spectrum. Figure 6 indicates that the frequency content of the dominant IMF φ 1 (t) during phase VI lies in the vicinity of the blade first bending mode frequency. This is consistent with the Fourier analysis in [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF].

In the progression from phase II to phase VI, significant variations in γ 1 are noticeable during phases III and V. The transition between phases III and IV matches a peak of amplitude in IMF φ 1 also evidenced by the red dot visible along γ 1 at the transition between the two phases at the center of the blue circle in Fig. 6.

Overall, it seems worth mentioning that the instantaneous frequency of IMF φ 1 is located in the vicinity of the frequency of the first torsional mode during the first two interaction phases. Because the fluctuation of γ 1 is consistent with the sudden drop in the experimental signals, a critical element in the understanding of the sequence from phase II to phase VI lies in the modal participation of the first two modes (bending and torsion) in the blade dynamics.

Modeling of contact interactions

Contact simulations are highly nonlinear since contact areas are a priori unknown and the respective hybrid contact pressure/displacement boundary conditions are part of the solution. The equation to be solved may be written as follows:

Mẍ(t) + D ẋ(t) + K(Ω)x(t) + C λ(t) = F a (t) (3)
where x is the blade displacement vector, M, D and K(Ω) are respectively the mass, damping and angular speed (Ω) dependent stiffness matrices of the blade. F a stands for the aerodynamic loading applied on the blade. Contact related variables involve the contact constraints matrix C and the contact forces λ. The latter are computed using the plastic constitutive law introduced in [START_REF] Legrand | Numerical investigation of abradable coating removal in aircraft engines through plastic constitutive law[END_REF] when contact with the abradable coating is detected. Would the abradable coating be fully removed where the casing is impacted, the contact forces stem from the well known Kuhn-Tucker unilateral contact conditions:

λ(t) ≥ 0 g(x(t), Ω, θ(t)) = c(x(t), Ω) -α(t)d(θ(t)) + w(t, θ(t)) ≥ 0 λ(t) • g(x(t), Ω, θ(t)) = 0 (4)
where g is the gap function and θ(t) is the angular position of a contact node. Before contact occurs, θ(t) = Ωt but once contact has occured, the blade tangential vibration yield θ(t) = Ωt. As illustrated in Fig. 7, the gap function g involves three quantities: (1) c(x(t), Ω) the clearance between a contact node and the initial circular casing profile which accounts for the blade vibration as well as centrifugal stiffening, (2) α(t)d(θ(t)) the product of an exponential function (see section 4.5) and the angular casing distortion both depicted in Fig. 10 and(3) the wear function w(t, θ(t)) which continuously reflects the actual contact interface with abradable removal stemming from previous blade/abradable coating contacts at θ(t).

Computation of the contact forces -parameter λ(t)

In order to solve Eq. (3), the employed numerical strategy relies on a central finite difference time integration scheme combined with: (1) a Lagrange multiplier approach [START_REF] Carpenter | Lagrange constraints for transcient finite element surface contact[END_REF] for blade/casing contact management and (2) a plastic constitutive law [START_REF] Legrand | Numerical investigation of abradable coating removal in aircraft engines through plastic constitutive law[END_REF] for blade/abradable coating contact treatment. At each time step, the displacement field is first linearly predicted before it is corrected if penetrations are detected within the abradable coating or the casing. Modeling assumptions and the unilateral contact procedure are not recalled here for the sake of brevity, details may be found in [START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF][START_REF] Legrand | Numerical investigation of abradable coating removal in aircraft engines through plastic constitutive law[END_REF].

Blade/casing clearance function -parameter c(x(t), Ω)

The experimental observations described in section 2 are related to a low-pressure compressor blade. However, other unpublished data suggest that the alternance between low and high amplitudes of vibrations in the vicinity of a critical speed is not design specific. Accordingly, a more generic blade Hilbert spectrum of the time response obtained from strain gauge 3. The two first eigenfrequencies of the blade f1B and f1T are superimposed in order to ease the interpretation of the spectrum

θ(t) c(x(t), Ω) α(t -δt)d(θ(t)) α(t)d(θ(t)) θ = 2kπ, k ∈ N circular casing profile (t = 0) Ω contact node blade (partial view) distorted casing profile (t) g(x(t), Ω, θ(t)) w(t, θ(t)) z (air flow direction)

Casing distortion -parameter d(θ(t))

Empirical industrial data suggest that the assembly of the casing with other components of the stator may lead to a globally ovalized profile of its inner surface. Accordingly, the casing considered in this study features two symmetrical (θ 2θ 1 = π in Fig. 10) contact areas-referred to as lobes in the following-along its circumference as pictured in Fig. 10. No significant amplitude of vibration were found on the casing during experimenal measurements [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF]. Consequently, the casing is modeled as a perfectly rigid mathematical profile on which is deposited the abradable coating. The prescribed distortion of the casing was identified in [START_REF] Batailly | Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF] as one fundamental ingredient for the occurence of blade/abradable coating interactions.

Modeling of abradable removal -parameter w(t, θ(t))

The modeling of the abradable coating relies on massless plastic elements introduced in [START_REF] Legrand | Numerical investigation of abradable coating removal in aircraft engines through plastic constitutive law[END_REF][START_REF] Batailly | Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF].

The clearance between the blade and the abradable coating at rest is constant (0.15 % of the blade height) along its tip.

Contact scenario -parameter α(t)

The following contact scenario applies to all the simulations presented in this article. At t = 0, the casing is perfectly circular, the blade is at rest and the separating clearance is positive. During the first 150 revolutions (until t 99% in Fig. 10), the casing undergoes distortion, the clearance is progressively closed and contact is initiated. The function α(t) smoothly distorts the casing in order to ensure that no blade/casing penetration exist at t = 0. From revolution 150 to the last one, the blade is forced to interact with the abradable coating within the statically distorted casing with α(t > t 99% ) = α. Moreover, this is not an intrinsic limitation of the proposed approach, but friction is not accounted for in this study.

4.6 Aerodynamic loading -parameter F a (t) By design, the pressure of the air flow inside the compressor increases from a stage to the next one. This variation of pressure between the air upstream of the blade and downstream of the blade is here achieved by the introduction of an aerodynamic pressure field F a (t) on the blade pressure side as schematically depicted in Fig. 11. A precise modeling of the aerodynamic loading applied on the blade goes beyond the scope of this article. Instead, the focus is made on the use of a simplified-yet representative-pressure field in order to better apprehend the consequences of having a forcing term in Eq. ( 3). It is assumed that the pressure field is only dependent on the surrounding downstream guide vanes pictured in Fig. 11. The pressure field F a (t) is also assumed to be sinusoidal in time:

F a (t) = F a sin(2πf a t) (5)
of frequency f a = Ω × N gv where N gv = 60 is the number of guide vanes at the considered stage.

As shown in Fig. 11b, the pressure loading is geometrically concentrated around the blade tip and the leading edge in agreement with industrial CFD simulations. In order to fully benefit from the blade reduced order model, the pressure field F a (t) is projected on the associated basis. The subsequent reduction error 5 ε achieved with this reduced-order model is ε < 0.01 %. The linear forced response of the blade is depicted in Fig. 12. It features a succession of vibratory resonances at specific critical angular speeds where amplitudes of vibration reach a maximum. One aim of the paper is to explore how these critical speeds predicted within a linear framework, see Fig. 12, compare with their unilateral contact induced counterparts. Equation.

(3) describes a nonlinear oscillator. The consequences of adding the periodic forcing term F a (t) in such equation is a wide field of research in many areas such as mathematics [START_REF] Gross | On oscillation types in forced nonlinear oscillators close to harmonic resonance[END_REF], electrical systems [START_REF] Jabli | Qualitative analysis of periodically forced nonlinear oscillators responses and stability areas in the vicinity of bifurcation cascade[END_REF] and biology [START_REF] Keener | Global bifurcations of a periodically forced nonlinear oscillator[END_REF] to name a few. In particular, the forcing term may be responsible for additional bifurcations and thus requires particular attention.

Vibratory response under vacuum

In this section, aerodynamic forces are omitted, F a (t) = 0 and vibrations solely stem from the interaction of the blade with the abradable coating or the casing. Simulations are carried out during 750 blade revolutions-a duration that is assumed long enough to reach steady state-over the angular speed range of section 4.2: Ω ∈ [0.159 ; 0.397]. The resolution of the time series is 1,000 points per revolution. In total, 1,500 evenly distributed angular speeds are considered over the aforementioned angular speed range. The contact forcing period T is half the revolution period T Ω of the blade since there are two symmetric lobes on the casing.

Bifurcation diagrams

Thus the simulation duration is T d = 750T Ω = 1500T . In the following, u LE (t) refers to the radial displacement of the blade leading edge. For each angular speed, the value u LE (T 0 ) is computed:

u LE (T 0 ) = max t∈[1450T ;1451T ] u LE (t) (6)
This value is plotted along with u LE (T 0 +iT ), i = 1, . . . , 49 as red dots in the bifurcation diagram 13a: this diagram is thus plotted using the 50 last periods of each simulation, which is assumed sufficient to accurately detect each type of motion. The corresponding final circumferential profile w(t = 750T Ω , θ) of the abraded coating is depicted in Fig. 13b: the colour code indicates the depth of the wear lobes from white (no removal) to red (maximum removal). For almost all velocities, a single dot is plotted in Fig. 13a which means that a perfectly periodic steady state is reached. There are, -in Fig. 14a, amplitudes undergo brutal variations at Ω = 0.1815 for angular speeds slightly higher than the first linear resonance (Ω = 0.178 = f 1B /6, see Fig. 13c). As mentioned in previous works, angular speeds Ω = f 1B /n or Ω = f 1T /n, n ∈ N * , are expected to induce vibratory resonances [START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF] together with n privileged contact areas along the circumference of the casing. This relationship is referred to as the linear interaction condition in the following. Higher angular speeds feature a response whose period is twice longer until a period-halving flip bifurcation occurs at Ω = 0.1835 after which the solution period matches the forcing period again. In Fig. 14b, this bifurcation is hardly noticeable in the abradable pattern; -in Fig. 15a, the solution features a flip bifurcation at Ω = 0.3105 with a transition from periodic solutions to a period-doubling sequence before a drop in amplitude is realized with deeper and wider lobes in the abradable coating, see Fig. 15b. Complementary blade eigenanalysis information is graphically provided in Fig. 13c. Crossings between the line Ω = f and sub-harmonics of the first two natural frequencies are found for angular speeds close to each bifurcation. The position of these crossings is shifted at angular speeds slightly lower than the aforementioned interactions which is a consequence of unilateral-contact stiffening.

Attention is now paid to the interaction diagram in interval Ω ∈ [0.307 ; 0.344]. Corresponding results are shown in Fig. 16a and Fig. 16b. The maximum amplitude is found at Ω 0.326 and unexpected low-amplitude periodic orbits with no contact separation for angular speeds slightly above the interaction frequency are found. The comparative inspection of Figs. 16a and16b underlines that the depth and width of the wear lobes are strongly dependent on the solution: in particular, low-amplitude steady states feature deeper and wider lobes. Coexisting stable periodic solutions call

Interactions of interest

The two full time responses and the corresponding steady states are shown in Figs. 17a and17b for Ω c and in Figs. 18a and18b for Ω * c . Steady states are detected through an autocorrelation procedure which automatically finds when periodicity is established in the response. Figs. 17a and 18a collectively expose two distinct vibratory responses. At Ω c , the amplitude grows after contact is initiated until a steady state is reached. However, at Ω * c the amplitude grows after contact initiation, then stabilizes before undergoing a rapid drop about revolution 460. The contact areas maps in Figs. 17c and18c provide additional exploitable information.

It is first evident that contact occurs at the leading and trailing edges during the same revolution with two contact areas per revolution between revolutions 0 and 100. While identical at rest, the leading edge and trailing edge clearances may differ once centrifugal forces are applied. This explains why contact initiation is not simultaneous at these two locations in Figs. 17c and18c. As mentioned in section 5.1, the following relationship is satisfied: At Ω = Ω c , the interaction scenario can be described as follows: contact is prescribed twice per revolution within the 2-lobe casing distortion and the blade dynamics is essentially driven by its first bending mode since contacts occur evenly along the blade tip. Thus, and in agreement with Eq. ( 7), blade vibrations are synchronous and four oscillations are observed per revolution. During the first 100 revolutions, amplitudes are such that a full clearance consumption is only reached where the casing is distorted and only two contact areas are visible-around angular positions θ 1 and θ 2 consistently with the lobes angular location in Fig. 10-in Fig. 17c. As the amplitudes increase, see Fig. 17a, there is a time for which contact is activated for each oscillation of the blade and two secondary contact areas appear in Fig. 17c around π and 2π after 100 revolutions for the trailing edge. Shortly after, around revolution 150, secondary contact areas are also detected at the leading edge. With four contact locations per revolution and four oscillations along the first blade bending mode, the linear interaction criterion is satisfied: a periodic steady state is reached with high amplitudes of vibration.

Ω c ≈ Ω * c ≈ f 1B /4 (7)
During the first 100 revolutions, similar observations are made for Ω = Ω * c : two contact areas are observed in Fig. 18c both at the trailing and the leading edges and the amplitudes of vibration increase as depicted in Fig. 18a. It is noticeable however that the amplitudes of vibration are lower in Fig. 18a than in Fig. 17a. Also, the time for which two secondary contact areas are detected on the casing is delayed for Ω * c (revolution 150 instead of revolution 100 for Ω c ). These secondary contact areas never arise at the leading edge. Between revolutions 150 and 400, secondary contact areas sporadically emanate, at the trailing edge only, and it is assumed that the blade dynamics is twofold: (1) bending motions are privileged when contacts occur simultaneously at the leading and trailing edges (for the main contact areas) whereas (2) torsional modes become dominant when contacts only occur at the trailing edge (for the secondary contact areas). This combination of bending and torsional modes implies that the linear interaction condition (7) cannot be satisfied. However, it remains unelucidated at this stage why such low amplitude of vibrations are predicted numerically.

Blade response analysis

Fourier transform

A Fourier transform of the two periodic steady states is performed and the corresponding spectra are pictured in Fig. 19a and Fig. 19b, respectively. The fundamental excitation frequency due to unilateral contact is f c = N Ω where N = 2 stands for the number of lobes in the casing distortion. In Figs. 19a and19b, high (possibly super-) harmonics in the two spectra are consistent with the fundamental excitation frequency: peaks of amplitude are located at frequencies f = kf c , k ∈ N only. As an example, the eighth super-harmonic noticeable in Fig. 18b is underlined in the spectrum depicted in Fig. 19b. One should notice the order of magnitude difference between the two vertical axes though. In agreement with the interaction scenario detailed above for Ω c , the dominant peak in Fig. 19a is found for f = 2f c = 2 × 2Ω = 4Ω, which is precisely the first blade eigenfrequency as stated in Eq. [START_REF] Padova | Experimental results from controlled blade tip/shroud rubs at engine speed[END_REF]. For Ω * c , the participation of higher frequency torsional modes leads to significant peaks of amplitude for f > f 1B . However, these Fourier spectra do not explain the emergence of distinct vibratory responses and the EMD of the full time responses is carried out with emphasis

Empirical Mode Decomposition

Numerical time responses are sampled at 1,000 points per blade revolution and the computation time required for their EMD is about 25 minutes 6 . The sifting process to obtain the first IMF represents about 80 % of the computation time and requires 331 iterations. Higher IMF are computed in a few seconds with less than 100 iterations.

The EMD of the time responses in Figs. 17a and 18a yields the first five IMF displayed in Figs. 20 and 21. Higher IMF exhibit negligible amplitudes and are therefore discarded. Interaction at Ω c is almost entirely represented by the second IMF while at Ω * c , the dominant IMF is φ 1 at least for the first 300 blade revolutions. Be reminded that the first IMF are essentially representative of high frequency contents. In light of this, the two interactions significantly differ in their high frequency content. At Ω c , the IMF amplitudes are remarkably constant throughout the considered time interval: only minor variations due to contact initiation are visible during the first 100 revolutions and small amplitudes are detected once steady state is activated, from φ 3 to φ 5 . Similar observations can be drawn at Ω * c for the first 310 revolutions just before the amplitude of φ 1 drops down and significant components are subsequently detected in φ 2 to φ 5 . This event and the initiation of secondary contact areas at the trailing edge, see Fig. 18c, are simultaneous. The EMD provides here a very useful insight on the blade dynamics. The first secondary contact areas at the trailing edge (between revolutions 150 and 170 in Fig. 18c) are found to have no effect on the amplitude of the first IMF. This indicates that the participation of torsional modes in the blade dynamics only becomes significant starting at revolution 310. Once the participations of torsional modes (and possibly higher frequency bending modes) are activated, their respective influence is revealed through IMF φ 2 to φ 5 . The combination of abradable removal, contact separation at the trailing edge and structural damping leads to lower amplitudes of vibration between revolutions 310 and 475.

Hilbert spectrum

The EMD continuously monitors the blade dynamics in the time/frequency domain and is thus highly suited to the analysis of nonlinear signals. The participation of each IMF in the blade response is calculated via the Hilbert spectra of the time responses at Ω c and Ω * c . These spectra are depicted in Fig. 22a and Fig. 22b respectively. Results are normalized with respect to the maximum of the leading edge displacement u * LE = max t ( u LE (t) ). At Ω c , the Hilbert spectrum features two main frequency components; that is f = f c and f = 2f c associated with φ 1 as well as a third minor high frequency component f 15f c also apparent in the Fourier spectrum, Fig. 19a. The latter is noticeably absent in the Hilbert spectrum at Ω * c in Fig. 22b. This spectrum displays two main components f = f c and f = 2f c during the first 300 revolutions. The amplitude of φ 2 drops here between revolutions 300 and 500 where a cloudy content with no dominant frequency is noticeable. Above 500 revolutions, the spectrum shows five major frequency components: f c , 2f c , 3f c , 4f c , and 8f c in good agreement with the Fourier spectrum depicted in Fig. 19b.

Overall, the spectra computed for each angular speed emphasize the crucial role of the first bending mode in both interactions since the harmonic f = 2f c is dominant at Ω c and Ω * c . EMD output confirms that the key difference between the two interactions is the possibility of synchronous blade vibrations along the first bending mode. However, the likelihood of such occurrence does not solely depend on the angular speed since, as pictured in Fig. 16a over Ω ∈ [Ω c ; Ω * c ], high and very low amplitude steady states alternate even for neighbor angular speeds. This points out the fact that the simulated rotor/stator interaction scenarios probably also depend on parameters other than the angular speed. The determination of these parameters is addressed in the following sections.

Conjectural bifurcation analysis

Besides the angular speed Ω, one parameter seems to have a critical influence in the blade amplitudes of vibration: the blade tip clearance configuration at rest (CCR) [START_REF] Salvat | Modeling of abradable coating removal in aircraft engines through delay differential equations[END_REF], which implicitely drives the clearance c(x(t), Ω) in Eq. ( 4).

In previous simulations, the CCR is constant along the blade tip and the major discrepancies in the contact scenarios at Ω c and Ω * c lie in the number of contact areas at the leading edge (see Figs. 23b with the same colour code 7 . For a given angular speed, the sensitivity to the CCR is obvious and a clearance modification of 5 % on the leading edge strongly affects the steady state amplitude. However, no correlation between δc and the level of vibration may be found. Two distinct branches for post-critical angular speeds are visible: one featuring a combination of very low amplitudes of vibration (LAV) including trajectories with contact separation and one with high amplitudes of vibration (HAV) that may jeopardize the structural integrity of the blade.

Summary

The progressive removal of the abradable coating dynamically drives the clearance configuration and the jumps in amplitude mentioned in section 2 and reported in [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF]-for another type of blade-may be interpreted in light of the aforementioned observations. The jumps in amplitude could be a consequence of a minor clearance evolution that initiates a jump from a stable solution (featuring low or high amplitudes of vibration) to another. Complementary investigations on the role of the CCR in the interaction are needed but it seems likely that it plays an important role for angular speeds slightly higher than the identified critical angular speeds where coexisting stable branches are found. In particular, the influence of the CCR should be analysed when an aerodynamic loading is applied on the blade accounting for the coupling between the CCR and the aerodynamic loading along the blade tip.

Vibratory response with aerodynamic loading

In this section, aerodynamic loading is accounted for and F a (t) = F a sin(2πf a t) in Eq. ( 3). An identical contact scenario to the previous section is considered. -a sequence of period-doublings/period-halvings is displayed in Fig. 25a, which resembles an incomplete doubling cascade (complete doubling cascades are known as Feigenbaum trees [START_REF] Thompson | Nonlinear dynamics and chaos[END_REF]). Nevertheless, the final abradable profile seems to be identical to one another throughout the selected angular speed range; -a flip bifurcation is observed in Fig. 26a: a periodic solution undergoes a period-doubling sequence. This flip bifurcation appears at Ω = 0.3095 and quite likely stems from the linearly predicted potential interaction point in Fig. 13c located at the crossing of the ninth sub-harmonic of the first blade torsional mode with the line Ω = f . This assumption is confirmed by the fact that nine lobes about 1.6 mm deep are systematically predicted on the casing (at the trailing edge in Fig. 26b) "during" the period-doubling sequence while only two lobes about 2 mm deep are found for Ω > Ω B ; -in Fig. 27a, a sudden increase of the blade amplitude of vibration is captured for Ω C with practically no consequence on the abradable profile, see Fig. 27b. The nature of this bifurcation is unclear but it could be a shift between two stable branches or related to internal resonance phenomenon. At this point, the lack of experimental data prevents any meaningful numerical/ experimental comparison. It is instead proposed to explore the nature of the solution. Previous results [START_REF] Batailly | Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF][START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF] show that the probability of existing steady-state periodic orbits is high. Additionally, results exposed in section 5 reveal the possible coexistence of periodic solutions. The comparison between bifurcation diagrams 13a and 24a illustrates that aerodynamic loading leads to even more sophisticated dynamics and the remainder qualitatively explores the dynamics with and without aerodynamic forces.

Predicted trajectories and their nature

From the extensive post-processing analysis of the results with aerodynamic loading, a few angular speeds are found for which state space phase diagrams and Poincaré maps with unexpected features could be identified. Four of these angular speeds are considered in the following, Ω 1 = 0.205, Ω 2 = 0.306, Ω 3 = 0.310, and Ω 4 = 0.383. For each of them, four quantities are depicted in Fig. 28: (1) the displacement u LE (t), (2) the displacement u LE (t) over the last period8 , (3) the phase diagram in the plane (u LE (t), uLE (t)) over the last 50 periods and (4) the associated Poincaré map. Periodic motion after contact separation At Ω = Ω 1 , the steady state is perfectly periodic and the amplitude of the blade leading edge axial displacement matches the amplitude predicted from the blade forced response due to the sole aerodynamic loading: this is an indication that 28c, at Ω = Ω 3 , the final level of vibration is much higher than the one computed solely with aerodynamic loading. Indeed, the blade dynamics seems mostly guided by contact and abradable removal conditions. The phase diagram features a more complex trajectory but the motion is still perfectly periodic with two dots in the Poincaré map. Likely to be chaotic motion Finally, when Ω = Ω 4 , the blade displacement is constrained by the abradable coating and the final amplitude is lower than the linearly forced one. The last period of this time response exhibits a very fluctuating amplitude of the signal which leads to a cloud-like phase diagram. The Poincaré map underlines the likely to be chaotic nature of the motion as no clear orbit appears. Eventhough a very large number of revolutions was considered, one may not completely rule out the possibility that a particularly long transient occurs for this particular regime.

Conclusion

Unanticipated fluctuations of measured stress signals during the experimental simulation of a blade/casing structural interaction form the core motivation for the numerical analysis carried out in this work. The proposed numerical simulations rely on an in-house tool which handles unilateral contact conditions between the reduced-order order finite element model of a rotating turbomachinery blade and the surrounding casing assumed perfectly rigid. Both centrifugal stiffening and abradable coating removal are accounted for and explicit time integration combined to a Lagrange-based contact procedure is preferred. The numerical strategy is improved along two distinct avenues: (1) simplified aerodynamic loading is accounted for and (2) the post-processing procedure makes use of empirical mode decomposition of the time histories: this brings a new insight on the predicted transient phenomena.

With and without external aerodynamic loading, introduced results exhibit very sophisticated dynamics with flip bifurcations and evidence of period-doubling and period-halving sequences Under vacuum, fundamental disparities divide interactions stemming from torsional and interactions stemming from bending modes: the former is found to initiate flip bifurcations with minor effect on the amplitude while the latter is responsible for very high peaks of vibration in the vicinity of which several stable solutions may coexist. Numerical results suggest that under vacuum, several periodic limit cycles coexist in the vicinity of an linearly predicted potential interaction point. The essential role of the operating clearance at rest separating the blade and the casing is highlighted: from a scan through various clearance configurations at rest emerge two distinct stable branches in the bifurcation diagram. By means of an empirical mode decomposition of the time domain responses, the interaction scenario is drawn. The outbreak of significant contributions of torsional modes within the blade dynamics is responsible for a sudden decrease in the amplitudes of vibration.

The incorporation of aerodynamic loading is key in order to give a plausible explanation for the fluctuations found in experimentally measured stresses. It seems worth mentioning that the coexistence of stable solutions for a given angular speed are not found with aerodynamic loading. Aerodynamic loading leads to even richer blade dynamics: both quasi-periodic and likely to be chaotic motions are predicted while such motion had never been observed before with the employed methodology.

In the end, the proposed numerical analysis draws a plausible interaction scenario: as the blade rotates and impacts the abradable coating, the progressive removal of the coating may initiate a jump from a branch to another with distinct amplitudes of vibration.

Acknowledgement

Thanks go to Snecma for its technical and financial support.

Figure 1 .

 1 Figure 1. Instrumented blade with three strain gauges

Figure 2 .

 2 Figure 2. Experimental stress histories: (a) strain gauge 1, (b) strain gauge 2 and (c) strain gauge 3. The experimental results presented in Fig. 2c are taken from [6], all other experimental data have never been published.

Figure 4 .

 4 Figure 4. EMD of the second test signal y2(t): (a) test signal, (b) IMF φ1, (c) IMF φ2, (d) IMF φ3 and (e) IF γ3

5 Figure 5 .

 55 Figure 5. EMD of the signal experimentally obtained from strain gauge 3: (a) IMF φ1(t), (b) IMF φ2(t), (c) IMF φ3(t), (d) IMF φ4(t) and (e) IMF φ5(t)

Figure 6 .

 6 Figure 6. Hilbert spectrum of the time response obtained from strain gauge 3. The two first eigenfrequencies of the blade f1B and f1T are superimposed in order to ease the interpretation of the spectrum

Figure 7 .Figure 8 .Figure 9 .

 789 Figure 7. Planar representation of the blade/casing contact configuration for a given contact node

Figure 10 .

 10 Figure 10. Normalized casing distortion with 2 lobes (dark colour) superimposed with the exponential function used for its application over time (light colour)

Figure 11 .

 11 Figure 11. Aerodynamic considerations : (a) partial view of a high-pressure compressor and (b) zoom over a blade

Figure 12 .

 12 Figure 12. Blade forced response of the leading edge radial displacement under aerodynamic loading

Figure 13 .

 13 Figure 13. (a) Bifurcation diagram, (b) corresponding abraded lobes, and (c) linearly predicted potential interaction points: red gradients show the influence of contact stiffening on the actual localization of the interaction

Figure 14 .

 14 Figure 14. (a) Bifurcation diagram (area A), and (b) associated abraded lobes however, a few exceptions for which distinct dots are visible in the bifurcation diagram. Exemplarily, zooms in Figs. 14a and 15a depict a few speeds of interest:-in Fig.14a, amplitudes undergo brutal variations at Ω = 0.1815 for angular speeds slightly higher than the first linear resonance (Ω = 0.178 = f 1B /6, see Fig.13c). As mentioned in previous works, angular speeds Ω = f 1B /n or Ω = f 1T /n, n ∈ N * , are expected to induce vibratory resonances[START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF] together with n privileged contact areas along the circumference of the casing. This relationship is referred to as the linear interaction condition in the following.Higher angular speeds feature a response whose period is twice longer until a period-halving flip bifurcation occurs at Ω = 0.1835 after which the solution period matches the forcing period again. In Fig.14b, this bifurcation is hardly noticeable in the abradable pattern; -in Fig.15a, the solution features a flip bifurcation at Ω = 0.3105 with a transition from periodic solutions to a period-doubling sequence before a drop in amplitude is realized with deeper and wider lobes in the abradable coating, see Fig.15b. Complementary blade eigenanalysis information is graphically provided in Fig.13c. Crossings between the line Ω = f and sub-harmonics of the first two natural frequencies are found for angular speeds close to each bifurcation. The position of these crossings is shifted at angular speeds slightly lower than the aforementioned interactions which is a consequence of unilateral-contact stiffening.Attention is now paid to the interaction diagram in interval Ω ∈ [0.307 ; 0.344]. Corresponding results are shown in Fig.16aand Fig.16b. The maximum amplitude is found at Ω 0.326 and unexpected low-amplitude periodic orbits with no contact separation for angular speeds slightly above the interaction frequency are found. The comparative inspection of Figs.16a and 16bunderlines that the depth and width of the wear lobes are strongly dependent on the solution: in particular, low-amplitude steady states feature deeper and wider lobes. Coexisting stable periodic solutions call

Figure 15 .Figure 16 .

 1516 Figure 15. (a) Bifurcation diagram (zoom on area B), and (b) associated wear lobes on the casing in front of the trailing edge

Figure 17 .

 17 Figure 17. Radial displacement of the leading edge at Ωc: (a) full time response transient ( ), periodic steady state ( ), (b) periodic steady state, (c) contact areas, left: leading edge, right: trailing edge

Figure 18 .

 18 Figure 18. Radial displacement of the leading edge at Ω * c : (a) transient response ( ), periodic steady state ( ), (b) periodic steady state, (c) contact areas, left: leading edge, right: trailing edge

Figure 19 .

 19 Figure 19. Fourier transform computed over the last 50 revolutions (a) at Ωc, (b) at Ω * c

5 Figure 20 .

 520 Figure 20. Time response and first five IMF calculated at Ωc

  17c

5 Figure 21 .Figure 22 .Figure 23 .

 5212223 Figure 21. Time response and first five IMF calculated at Ω * c

Figure 24 .Figure 25 .

 2425 Figure 24. (a) Bifurcation diagram and (b) associated wear lobes at the trailing edge

Figure 26 .

 26 Figure 26. (a) Bifurcation diagram (zoom on area B) and (b) associated wear lobes at the trailing edge

Figure 27 .

 27 Figure 27. (a) Bifurcation diagram (zoom on area C) and (b) associated wear lobes at the trailing edge

Figure 28 .

 28 Figure 28. Superimposition of time response, a zoom over its last period, phase diagram and Poincaré map with aerodynamic loading. ( ) amplitude of the blade forced response due to the aerodynamic loading without structural contacts.

A discrete time signal is said to be nonlinear when it is associated to a nonlinear mechanical system governed by nonlinear ordinary differential equations.

In the article, the angular speed refers to the rotational frequency of the rotor.

The end of phase VI slightly varies from a figure to another depending on when the signal acquisition was lost on each gauge.

The discrimination between low and high frequencies applies locally in the context of an EMD as underlined in Fig.4d.

The reduction error is the relative comparison of the strain energy associated with the static deformation of the finite element model undergoing the pressure loading with the strain energy associated with the static deformation of the reduced order model undergoing the reduced pressure field.

Approximated computation times on a PC with an i7-2600 3.4 GHz processor.

All simulations lead to a periodic steady state and for every angular speed, a single dot is visible in Fig.23bfor each CCR.

The last period is half a revolution due to the two symmetric contact areas along the casing circumference, here the last half-revolution is shown: from revolution 749.5 to revolution 750