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

This article deals with the numerical investigation of the unilateral contact-induced dynamics of a turbomachine blade rotating within a perfectly rigid yet distorted casing. This investigation is motivated by unelucidated vibratory behaviours observed experimentally. The simulations are based on an in-house time-marching strategy incorporating Lagrange multipliers for the unilateral contact treatment, as well as centrifugal stiffening and abradable coating removal. Significant extensions are proposed through the implementation of (1) aerodynamic loading on the blade and, (2) post-processing techniques involving the empirical mode decomposition which provides fruitful insights on important transient phenomena. A thorough bifurcation analysis with and without aerodynamic loading highlights the existence of flip bifurcations with period-doubling and period-halving sequences over a broad angular speed range. Numerical simulations with external aerodynamic loading yields quasi-periodic and likely to be chaotic motions that could not be observed under vacuum. The proposed numerical investigations underline the key role of the aerodynamic loading in the blade dynamics and suggest that unexplained experimental vibratory behaviours are related to the vacuum conditions of the experiment.


Introduction
In turbomachinery, fan, compressor and turbine blades are highly prone to structural failures [1]. 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 [2] were targetted in various numerical and experimental investigations related to the aerospace [3,4] and electric power industries [5]. Intricate and possibly diverging dynamics going from a single blade [6,7] or a full bladed disk assembly [8] to several compressor stages with precessional shaft motions [3] are found and several partially unelucidated industrial incidents are reported [6,9]. 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 [10,6]. Possible contacts between neighbouring blades [11,12] 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 [6] that remain unelucidated. A scenario which assumes a brutal and significant alteration of the abradable material properties is investigated in [10] 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 [6,10] 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 [13,14] which is well-suited to non-stationary and nonlinear 1 time responses even though its mathematical foundation remains partially open [15]. 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) [14] have been extensively applied for monitoring and diagnosis purposes [16] 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 [13] 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 [17] 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 [6]. 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 [18], phase diagrams and Poincaré maps (which are one snapshot of the phase space per forcing period).
In the first section, the selected interaction scenario [10] is briefly summarized along with the relevant experimental observations [6]. 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 1 A discrete time signal is said to be nonlinear when it is associated to a nonlinear mechanical system governed by nonlinear ordinary differential equations. 2 In the article, the angular speed refers to the rotational frequency of the rotor.
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 [6] 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 [6]. 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 s 3 ): 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 [10] 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 [19,20] and numerically [21] were motivated by the interaction phenomenon reported in [6]. 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 [22,23] 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 [13] advantageously provide physically relevant results and have already been successfully used in the field of rotordynamics [16,24]. 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): 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 [13]: 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 [13].
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 [14] and (2) collecting the envelope functions a k (t). The computation of the envelope functions depends on the EMD principle as used in [25] context in which the EMD is performed. As an example, their acquisition is associated to a moving average algorithm in the LMD [14] while it relies on a B-spline interpolation of local extrema in the Hilbert-Huang transform [13].
Similarly to the Hilbert-Huang transform [13], 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: is sampled at f = 4 kHz and is plotted in Fig. 4a for t ∈ [0 ; 30π]. For large t its frequency tends to f = 3/2π 0.477 Hz. Only the first three IMF plotted in Figs. 4b, 4c and 4d participate in the tested function. From the EMD procedure, a single non-negligible intrinsic mode function φ 3 is obtained as illustrated in Figs. 4b, 4c and 4d. The peculiarity of the IMF lies in the fact that they are well suited for a Hilbert transform [13]. From every IMF φ i (t) the associated instantaneous frequency (IF) γ i (t) is retrieved through a Hilbert transform: γ 3 (t) is pictured in Fig. 4e. It is worthy to notice that the  low-frequency content of y 2 (t) is stored in IMF φ 3 while its high-frequency content 4 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 [6].
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: 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 [26] 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: 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 [27] for blade/casing contact management and (2) a plastic constitutive law [26] 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 [17,26].

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 I  II  III  IV V  VI   110  120  130  140  150  160  170  180    engine high-pressure compressor blade. It is clamped on its root as pictured in Fig. 8. Its finite element model contains 13,803 nodes and 7,388 tetrahedron quadratic finite elements. It is reduced through modal synthesis accounting for centrifugal stiffening [28,10] based on the Craig-Bampton method [29]. Eight evenly spaced nodes are selected along the tip for contact purposes. The full-scale model of 42,000 degrees-of-freedom (DoF) reduces to a 95 DoF model including 24 physical DoF for contact management and 71 modal coefficients. Both space (modal) and time domain convergences of the proposed model were carefully checked but are not detailed.
The Campbell diagram pictured in Fig. 9 is calculated with the reduced-order model. Within the range of interest Ω ∈ [0.159 ; 0.397], centrifugal stiffening yields a very significant increase of the first two eigenfrequencies of the blade, respectively the first bending (1B) and first torsional (1T) modes.

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 [6]. 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 [10] 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 [26,10].
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.

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.  Fig. 11. The pressure field F a (t) is also assumed to be sinusoidal in time: 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 [30], electrical systems [31] and biology [32] 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. 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:

Bifurcation diagrams
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 [17] 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 and 16b 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 and 17b for Ω c and in Figs. 18a and 18b 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 and 18c 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 and 18c. 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.
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   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.

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 and 19b, 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. (7). 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) [21], 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 at Ω * c and 18c). Accordingly, departure from this reference configuration depicted as the horizontal straight line in Fig. 23a is realized by disturbing the CCR of a magnitude δc at the leading edge, δc/2 at the midchord (MC) and 0 at the trailing edge. A quadratic interpolation describes the clearance at the remaining contact nodes. As pictured in Fig. 23a, seven perturbed operating clearances from δc = −25 % to δc = +25 % are considered.  Fig. 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 [6]-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.

Time domain simulations
As in section 5.1, a bifurcation diagram is plotted in Fig. 24a along with the wear lobes in Fig. 24b. The comparison between Figs. 13a and 24a emphasizes the emergence of new resonance peaks when aerodynamic loading is accounted for. These peaks may be predicted both from linear considerations with the forced response plotted in Fig. 12 and nonlinear considerations with the Campbell diagram in Fig. 13c. For many angular speeds, more than a single dot are visible illustrating the potential existence of non periodic orbits. To support this observation, the blade dynamics is zoomed in on the diagrams provided in Figs. 25a, 26a, and 27a: -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 [33]). 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 [10,17] 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 period 8 , (3) the phase diagram in the plane (u LE (t),u LE (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 Quasi-periodic motion At Ω = Ω 2 , the blade response in Fig. 28b exhibits a quasi-periodic motion characterized by a closed orbit in the Poincaré map. The final amplitude of vibration is lower than the one under aerodynamic loading: the abradable coating constrains the blade displacement and intermittent contact occurs. Once a quasi-periodic steady state is established, the contact configuration periodically repeats itself and the abradable elements undergo elastic deformations only. Periodic motion As illustrated in Fig. 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)  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.