Full Three-Dimensional Rotor/Stator Interaction Simulations in Aircraft Engines With Time-Dependent Angular Speed

Modern aircraft engine designs feature reduced clearances that may initiate structural contacts between rotating and static components. A numerical strategy dedicated to the simulation of such interactions is here enriched in order to account for time-dependent angular speeds. This contribution first details the evolution of the numerical strategy before validating the developments by comparing numerical results with experimental observations made on an industrial test bench. Further numerical investigations allow to assess the sensitivity of numerical results to acceleration and deceleration rates. Results, obtained with and without abradable coating, underline the fundamental nonlinear nature of the analysed system. It is found that lower acceleration rates favour the arisal of interaction phenomena and that amplitudes of vibration at a given angular speed are generally lower when the blade decelerates.


Introduction
Increasingly tight operating clearances between rotating and static components of aircraft engines advantageously increase their overall efficiency but also lead to more frequent structural contacts. Consequently the latter cannot be only considered in accidental configurations anymore. These contact events must be accounted for as early as the design stage of bladed components [1] since they might initiate potentially damaging interaction phenomena such as rubbing [2], modal interaction [3] or whirl motions [4]. Accordingly, engine manufacturers are focusing on the development of robust numerical strategies [4,5,6,7] in order to get a better understanding of the physical phenomena at play when such contacts occur. A review of existing numerical strategies dedicated to the simulation of rotor/stator interactions-with a focus on the blade-tip/casing interface-may be found in [8]. Recent numerical investigations [9,10] led to promising results with a good agreement between predicted results and experimental observations both on low-and high-pressure compressor blades. In an attempt to provide a more general numerical strategy for the prediction of rotor/stator interactions initiated by structural contacts, this paper focuses on the numerical simulation of such interactions featuring a time-dependent angular speed. Both the stiffness of the rotor and the blade tip/casing clearances are updated live during the simulation. Accounting for a time-dependent angular speed allows for an extension of previous confrontations between numerical results and experimental observations, as more realistic interaction scenarios can be numerically simulated. The impact of this new feature is assessed with and without considering the deposition of an abradable coating [11] on the casing. Noticeable vibration behaviours of the blade are underlined in each case depending on whether it undergoes an acceleration or a deceleration. All numerical simulations are carried out with the industrial finite element model of a low-pressure compressor blade of a civil aircraft engine-see Fig. 1-that was previously investigated from both an experimental [2] and a numerical standpoint [9]. The interaction configuration of interest involves a single blade clamped on its foot and a pefectly rigid casing, thus it may be assimilated to a rubbing event [12,13,14]. Even though the proposed numerical developments are applied in the context of an aircraft engine, the proposed methodology is general and could be applied with no effort to an industrial steam turbine for which the simulation of rubbing events during speed transients is of particular interest [15,16]. The first section of this article introduces the numerical developments required to account for a time-dependent angular speed. Numerical results are then confronted to experimental observations before extensive exploratory investigations are carried out in the third section.

Existing numerical strategy
This article focuses on an evolution of the numerical strategy previously introduced in [17]. Thus, the modeling assumptions inherent to this numerical strategy apply: in particular, the simulations are carried out within the framework of linear small perturbation theory. Also, as a first approximation, angular acceleration stiffening terms are neglected. The reduced-order models computed for the rotor (be it a blade or a full bladed disk) account for centrifugal effects over a given angular speed range [0 ; Ω max ]. The corresponding stiffness matrix is dependent on the angular speed Ω: where matrices K 0 , K 1 , and K 2 are a combination of stiffness matrices computed at rest and for Ω = Ω max /2 and Ω = Ω max [18]. Potential structural contacts may only occur along the blade tip, where the boundary nodes used for contact management are located, as depicted in Fig. 1. The evolution of blade tip/casing clearances due to the angular speed is also accounted for. The clearances of boundary nodes along the blade tip are gathered in the vector c.

Proposed evolution
It is assumed that the angular speed Ω(t) is time dependent and is known at all time. At each time step of the explicit time integration procedure, two quantities must then be updated: (1) the stiffness matrix K(Ω(t)) and (2) the blade/tip clearance vector c(t).

Time integration procedure
The explicit time intergration procedure relies on the central finite difference method. At each time step n + 1, displacements u n+1 are linearly predicted as follows: If contact is detected at this time step, the following correction is applied: where λ stands for the contact forces computed as a function of the penetrations numerically predicted between the blade and the casing. Under the assumption that the damping matrix D is not dependent on the angular speed, only matrix A of the time integration procedure depends on Ω. The stiffness matrix appears as a multiplying term in Eq. (2) that may be written as: The explicit definition of K(Ω) in Eq. (1) yields: Because Ω(t) is known at each time step, the matrix A(Ω) may be updated with minimum impact on computation times.

Computation of clearances
In order to obtain c(t) at all time, the following procedure is considered: 1. before initiating the time integration procedure, boundary nodes displacements are computed for Ω = Ω max /3, Ω = 2Ω max /3, and Ω = Ω max (in addition, displacements of the boundary nodes for Ω = 0 are obviously null), 2. for each node and each angular speed, the blade tip/casing clearance is computed based on the aforementioned displacements, 3. values obtained for Ω = 0, Ω = Ω max /3, Ω = 2Ω max /3, and Ω = Ω max allow for a cubic interpolation of the blade tip/clearance for each boundary node depending on Ω, 4. inside the time integration procedure, at every time step, the value of Ω(t) is used in the cubic interpolation to obtain corresponding blade tip/casing clearances. The added computation time related to this procedure is relatively limited. It is estimated that numerical simulations with a time-dependent angular speed are about 15 % longer than simulations featuring a constant angular speed.

Numerical convergence
A detailed convergence analysis of the employed numerical strategy goes beyond the scope of this article. Indeed, the convergence of numerical results must be checked with respect to a large number of parameters such as: (1)  Convergence analyses for each of the aforementioned parameters may be found in previous publications [11,17,9], it is only recalled here with respect to the time step of the time integration scheme since the proposed developments impact the time integration procedure.
The time history of the angular speed considered for this convergence analysis is plotted in Fig. 2. It is centered around an angular speed Ω c related to the experimental scenario investigated in section 3. Contact between the blade-tip and the casing is initiated through a deformation (ovalization) of the casing over the first 20 % of the simulation. Results presented in the following are obtained considering four time steps: h = 5 · 10 −7 s ( ), h = 2 · 10 −7 s ( ), h = 1 · 10 −7 s ( ) and h = 5 · 10 −8 s ( ). First of all, the radial displacement of the blade leading edge over time is pictured in Fig. 3. The full time history is depicted in Fig. 3(a) and a zoom over the green area is made in Fig. 3

(b).
A typical all or nothing type of convergence is observed. For a too large time step, the numerical simulation diverges soon after the first contacts, see the red line ( ) in Fig. 3(a), while for small enough time steps, time histories are almost perfectly superimposed. The zoom in Fig. 3(b) underlines that a non negligible shift exists between the results obtained from time simulations with h = 2 · 10 −7 s and the results obtained with smaller time steps. This shift results is the consequence of a variation in the detection of the first contact due to the modification Batailly et al.  Similar conclusions may be drawn from the contact force time histories plotted in Fig. 4. The numerical convergence of contact forces is fundamental in order to provide meaningful results, it is here achieved for h = 1 · 10 −7 s ( ). In the remainder, all numerical simulations are carried out with the time step h = 1 · 10 −7 s.

Numerical/experimental validation
Results obtained experimentally in [2] and numerically in [9] are used as a reference point for the validation carried out in this section. The interaction scenario is a rubbling event between a single blade-depicted in Fig. 1-and a casing on which is deposited an abradable coating. Because no significant level of vibration was measured on the casing during the interaction, it is numerically modeled as a slightly ovalized perfectly rigid casing. The evolution of the angular speed is pictured in Fig. 5 where Ω c is the targeted critical angular speed mentioned in [2].
For the sake of brevity, experimental results are only briefly recalled here. During the blade acceleration, structural contacts are initiated due to centrifugal stiffening. Soon after the blade angular speed reaches Ω c , increased amplitudes of vibration are measured on the blade: distinct phases featuring higher and lower amplitudes of vibration alternate until the amplitudes reach a critical point leading to blade failure soon before the blade decelerates. A close inspection of the casing after the interaction reveals a very specific wear pattern of the abradable   [2] coating along the circumferential direction: two large lobes are measured in front of the leading edge while six deeper lobes are measured in front of the trailing edge.
In a first attempt to predict this interaction, numerical simulations carried out with a constant angular speed Ω Ω c [9] provided consistent results with experimental observations: (1) predicted amplitudes of vibration were much higher at the trailing edge than at the leading edge, (2) two wear lobes were predicted in front of the leading edge while six deeper lobes were predicted in front of the trailing edge, and (3) maximum stress areas within the blade profile matched cracked areas experimentally observed.
With the same blade model as the one used in [9], and accounting for a time-dependent angular speed, the time history Ω(t) depicted in Fig. 5 is used for the validation of the proposed numerical strategy.  Fig. 6 for the boundary node at the leading edge ( ) and at the trailing edge ( ). When the normalised distance is lower than 5, the boundary node penetrates the abradable coating. When the normalised distance is 0, the boundary nodes impacts the casing surface which means that the abradable coating was locally fully removed. Soon after the blade reaches Ω c , a noticeable increase of the amplitudes of vibration is predicted at the trailing edge while amplitudes of vibration at the leading edge remain constant throughout the interaction. In order to push forward the analysis of the numerical results, the wear profiles predicted on the casing are pictured in Fig. 7.
Wear profiles predicted at different instants of the interaction are depicted with a colour code from white to dark blue. The black wear profile is the wear profile at the end of the simulation. In agreement with experimental observations, two very distinct types of wear patterns are predicted at the leading edge-where two small lobes are predicted-and at the trailing edge-where six deeper lobes are predicted. Wear areas between the main lobes that are visible on the final wear profile at the trailing edge ( ) result from the fast deceleration of the blade leading to a non synchronous vibration featuring sporadic abradable removal along the casing circumference.
Overall, the presented numerical results are in good agreement with experimental observations and validate the proposed numerical developments. The possibility to account for a time-dependent angular speed opens avenue for finer numerical investigations that will complement previous studies regarding the co-existence of stable solutions [19]. In particular, in the remainder of this article, the focus will be made on distinct vibration behaviors predicted for an acceleration or a deceleration of the blade.

Exploratory investigations. . .
Numerical investigations are all carried out with the same low-pressure compressor blade as in the previous sections. Presented results are separated in two categories: (1) those obtained without abradable coating deposited on the casing and (2) those obtained with an abradable coating. From a numerical standpoint, the distinction between these two configurations must be underlined: when an abradable coating is modeled, contacts that occur at a given Batailly et al.  time and a given angular speed will modify the contact surface,-since some particles of the abradable coating are removed-and thus affect the contact configuration for the subsequent revolutions. To the contrary, without a coating, the casing surface is perfectly rigid and the contact configuration is independent of time. When it is accounted for, the removal of abradable coating is modeled using the plastic constitutive law [11] previously used in [9]. Simulations are carried out over a large time interval t ∈ [0 ; T ] over which the blade can do several hundred revolutions. For the sake of confidentiality, T = 1 in the remainder. Three angular speeds are considered in the following: Ω L < Ω c , Ω H,1 > Ω c and Ω H,2 > Ω H,1 . Based on these three angular speeds, four interaction scenarios are considered: two accelerations and two decelerations of the blade as depicted in Fig. 8. Finally, for each simulation, a spectrogram of the leading edge radial displacement is plotted. The colour code is identical for all the spectrograms shown in the article: from white (no amplitude of vibration) to black (maximum amplitude of vibration).

. without abradable coating
Considering angular speeds Ω L and Ω H,1 , the spectrograms of the blade's radial displacement over time are pictured in Fig. 9(a) and in Fig. 9(b) respectively for an acceleration and a deceleration of the blade. Superposed onto these spectrograms, four transparent colored bands indicate the location of the blade first four eigenfrequencies (1B, 1T, 2T and 2B) between the angular speeds considered in the interaction scenario. Highest amplitudes of vibration, pictured as black areas in the spectrograms are located in the vicinity of these colored bands. In particular, most elevated amplitudes are found around colored bands related to modes 1B and 2B for both interaction scenarios. Overall, three interactions are found in Fig. 9(a): 1. for t ∈ [0.05 ; 0.25], high amplitudes of vibration are predicted around the 2B colored band (area A), 2. for t ∈ [0.65 ; 0.8], here again, high amplitudes of vibration are predicted around the 2B colored band. These significant amplitudes of vibration are found soon after the blade angular speed has crossed the critical value Ω c which happens for t 0.6, see Fig. 8 (area B), 3. for t ∈ [0.9 ; 0.95], a narrow black spot is visible in the vicinity of the 1T colored band which indicate an elevated amplitude of vibration related the first torsion mode (area C). When looking at the spectrogram plotted in Fig. 9(b), it is noticeable that the interactions areas are significantly different (areas D, E and G). The torsion induced interaction around the 1T colored band is still visible (area D). Though it is now located around t 0.05 since the blade decelerates. Also, the elevated amplitudes of vibration around the 2B colored band for lower angular speeds (area G) are still predicted, but the interactions areas are larger and located over t ∈ [0.78 ; 1]. However, no significant amplitudes of vibration are predicted after the blade angular speeds crosses the critical value Ω c which, in this deceleration scenario, happens for t 0.4, see Fig. 8 This way, this first comparison between acceleration and deceleration of the blade underlines the sensitivity of the predicted interactions to the initial conditions, which is characteristic of a highly nonlinear system. In order to assess the sensitivity of the presented results to the acceleration or deceleration of the blade, additional spectrograms are plotted in Fig. 10(a) and 10(b) for angular speeds Ω L and Ω H,2 . The duration of this contact scenario is identical to the previous one but the gap between the two considered angular speeds is higher thus blade acceleration and deceleration are greater. Because of this larger gap between the angular speeds, the frequency width of the colored band is larger in Figs. 10(a) and 10(b) than it was in Figs. 9(a) and 9(b). For both acceleration and deceleration interaction scenarios, there is no noticeable amplitudes of vibration for lower angular speeds: for t ∈ [0 ; 0.6] in Fig. 10(a) and for t ∈ [0.4 ; 1] in Fig. 10(b). For higher angular speeds, a few elevated amplitudes of vibration are visible around the 1T and 2B colored bands (areas H, I, K and L). Again, the crossing of the blade angular speed with the critical speed Ω c , leads to small yet significant amplitudes of vibration along the 1T colored band (area J) when the blade accelerates-for t 0.44 in Fig. 10(a)-while no peak can be found when the blade decelerates-for t 0.6 in Fig. 10(b).
The numerical results presented in this section suggest that the acceleration and deceleration rates of the blade have a very strong impact on the arisal of rotor/stator interactions: the amplitudes of vibration predicted for low acceleration and deceleration rates are indeed much higher than in the case of high acceleration and deceleration rates.

. . . with abradable coating
An abradable coating is now considered along the casing contact surface. The interaction scenarios are identical to the ones considered in the previous section. The two first spectrograms are related to acceleration and deceleration between Ω L and Ω H,1 . Same as in the previous section, the casing is slightly ovalized, meaning that there are two privileged contact areas along its circumference. The difference between spectrograms plotted in Fig. 11(a) and Fig. 11(b) is patent: when the blade accelerates, very significant amplitudes of vibration are predicted after its angular speed crosses the critical speed Ω c around t 0.6 (areas M and N) while only non negligible amplitudes of vibration are found (area O), around t 0.4, when the blade decelerates. Corresponding amplitudes are more than one order of magnitude lower than when the blade accelerates.
Another difference between the two interaction scenarios lies in the fact that at the end of the simulation, contact has been lost when the blade accelerates while contact occur throughout the simulation when the blade decelerates. Indeed, the very high amplitudes of vibration predicted over a large time interval when the blade accelerates leads to an excessive abradable removal and thus opens blade tip/casing clearances which may lead to contact loss. To the contrary when the blade decelerates, the fact that only small amplitudes of vibration are predicted hint that abradable coating removal is very small and contact is maintained throughout the simulation.
Instead of focusing on the acceleration or deceleration rates, it is here proposed to analyse the impact of a different casing deformation over the predicted results. To this end, the casing is now distorted along a 3-nodal diameter free vibration mode and thus features three symmetric privileged contact areas along its circumference.
Spectrograms obtained for this configuration are plotted in Figs. 12(a) and 12(b). Similarly to the previous interaction scenario with two privileged contact areas, contact is lost at the end of the simulation when the blade accelerates while it is maintained throughout the simulation when the blade decelerates. Beside of this observation, results pictured in Figs. 12(a) and 12(b) significantly differ from previous observations. In fact, significant amplitudes of vibration are predicted both when the blade accelerates and decelerates (areas P, Q, R and S). In particular,  when the blade decelerates, high amplitudes of vibration predicted around the 1B and 2B colored bands are likely to be related with the crossing of the blade angular speed with the critical value Ω c . Overall, the results obtained with abradable coating confirm results obtained without the coating and underline the distinct vibration behaviour of the blade depending on the fact that it accelerates or decelerates. The last spectrogram plotted in Fig. 13 is obtained while the blades undergoes a very sudden acceleration from Ω L to Ω H,1 for t ∈ [0.45 ; 0.55].
No significant amplitudes of vibration are predicted during the acceleration, which confirms that a higher acceleration rate is less propitious to the arisal of rotor/stator interactions. The distinct frequency content of the blade's dynamics response may then be observed between Ω L , for which low frequencies are dominant, and Ω H,1 , for which the blade response is located around colored bands 1B and 1T.  through the confrontation of numerical results with experimental observations. Further numerical investigations are then carried out with the 3D finite element model of a low-pressure compressor blade. Predicted results underline the great sensitivity of the results to acceleration and deceleration rates. The highly non-linear nature of the simulated interactions is highlighted since interaction areas witnessed for an acceleration or a deceleration are different. In the end, the proposed developments open avenues for the numerical simulation of more sophisticated interaction scenarios in order to get a better understanding of the physical phenomena at play.