Modeling of abradable coating removal in aircraft engines through delay differential equations

In modern turbomachinery, abradable materials are implemented on casings to reduce operating tip clearances and mitigate direct unilateral contact occurrences between rotating and stationary components. However, both experimental and numerical investigations revealed that blade/abradable interactions may lead to blade failures. In order to comprehend the underlying mechanism, an accurate modeling of the abradable removal process is required. Time-marching strategies where the abradable removal is modeled through plasticity are available but another angle of attack is proposed in this work. It is assumed that the removal of abradable liners shares similarities with machine tool chatter encountered in manufacturing. Chatter is a self-excited vibration caused by the interaction between the machine and the workpiece through the cutting forces, and the corresponding dynamics are efﬁciently captured by delay differential equations. These equations differ from ordinary differential equations in that previous states of the system are involved in the formulation. This mathematical framework is here employed for the exploration of the blade stability during abradable removal. The proposed tool advantageously features a reduced computational cost and consistency with existing time-marching solution methods. Potentially dangerous interaction regimes are accurately predicted and instability lobes match both ﬂexural and torsional modal responses. Essentially, the regenerative nature of chatter in machining processes can also be attributed to abradable coating removal in turbomachinery.


Abstract
In modern turbomachinery, abradable materials are implemented on casings to reduce operating tip clearances and mitigate direct unilateral contact occurrences between rotating and stationary components. However, both experimental and numerical investigations revealed that blade/abradable interactions may lead to blade failures.
In order to comprehend the underlying mechanism, an accurate modeling of the abradable removal process is required. Time-marching strategies where the abradable removal is modeled through plasticity are available but another angle of attack is proposed in this work. It is assumed that the removal of abradable liners shares similarities with machine tool chatter encountered in manufacturing. Chatter is a self-excited vibration caused by the interaction between the machine and the workpiece through the cutting forces, and the corresponding dynamics are efficiently captured by delay differential equations. These equations differ from ordinary differential equations in that previous states of the system are involved in the formulation. This mathematical framework is here employed for the exploration of the blade stability during abradable removal.
The proposed tool advantageously features a reduced computational cost and consistency with existing time-marching solution methods. Potentially dangerous interaction regimes are accurately predicted and instability lobes match both flexural and torsional modal responses. Essentially, the regenerative nature of chatter in machining processes can also be attributed to abradable coating removal in turbomachinery.

Introduction
Delay differential equations (DDE) differ from ordinary differential equations (ODE) in that the system of interest involves lags [1]. Retarded DDE, that only involve the delayed state of the system, are usually distinguished from advanced and neutral DDE, for which the equation involves both delayed states and delayed rates of change. Furthermore, the delays may either be positive constants, time dependent, state dependent or distributed variables [2]. Only retarded DDE with a single constant delay will be considered in this paper: where t is the time, is the so-called time delay and P z denotes the derivative of z with respect to time.
Accordingly, initial conditions of a DDE are a function defined over the interval OE I 0. These types of equations are used in many fields such as species dynamics [3], chemistry [4] and torque control within combustion engines [5] to name a few. In manufacturing, DDE have proved to be very suitable for the modeling of the cutter-workpiece interaction phenomena occurring in various manufacturing processes such as milling [6] and turning [7]. In these research areas, it is highly desirable to avoid self-excited vibrations between the cutter and the workpiece, also known as chatter effect, in order to obtain smooth surfaces without major wear of the tool. The equations of motion are written as: As depicted in Fig. 2, milling shares similarities with aircraft engine structural interaction events. Blade-tips may slightly cut the surrounding abradable coating due to vibrations of the blades or the casing. Abradable coatings are implemented as a sacrificial material deposited on the casing in order abradable casing bladed-disk workpiece material to be removed cutter Figure 2: Analogy between abradable removal and milling to minimize operating clearances. They should be sufficiently resilient to endure severe thermal conditions and hostile constraints but also adequately soft not to alter the blades' structural integrity. Both experimental [8] and numerical [9] investigations have shown that blade/abradable contacts may lead to blade structural failure thus jeopardizing safety. An accurate modeling of the abradable coating removal is of prime importance to better understand the initiation of the divergence. To the best of our knowledge, two models are proposed in the literature [10,11] and a good agreement with experimental data was found. However, the proposed solution methods involve costly time-stepping techniques that are not well adapted to fast sensitivity analyses.
An improved understanding of the interaction calls for a qualitative model to quickly identify critical rotational speeds or critical material parameters for instance. On that account, the present study aims at taking advantage of the DDE framework to explore the stability of the equilibrium position of the abradable-inccuring blades. It is assumed that the axis of rotation of the shaft is perfectly rigid and contact initiation between the blade-tip and surrounding abradable coating originates from casing distortion or blade vibration only. Even if actual blade/abradable contacts are not permanent, it is here assumed that the blade never separates from the abradable layer during the revolution considered for the stability analysis, in agreement with what is assumed in the modeling of the milling and turning processes 1 .
First, the respective mathematical background is briefly recalled. The second section is devoted to the industrial model of the targeted compressor blade, also, main assumptions, limitations, and properties of the developed model of abradable coating removal through DDE are introduced. A convergence study is carried out in the third section in order to ensure the accuracy of the results. Finally, the results are compared with existing time-domain simulations.
2 Delay differential equations 2.1 Theoretical background As previously mentioned, DDE are a mathematical tool commonly used to model the dynamical behavior of systems involving lags, i.e. the evolution of the system P z.t/ not only depends on its current state z.t/, but also on a past state z.t / as specified in Eq. (2). In the particular case of mechanical systems, the forcing term f in Eq. (3) is often written as: where .t/ is an external forcing matrix which may be time dependant. By introducing state-space in order to obtain a first order DDE, Eqs. (3) and (4) yield: with, 1 The relative displacement between the tool and the workpiece usually makes it possible to assume that, at any time t , the cutter is in contact with the workpiece. Only a few studies consider cutter/workpiece detachment [12]. When the aforementioned conditions are satisfied, the first order system given in Eqs. (5) and (6) is propitious to a stability analysis of the z D 0 solution, also referred to as equilibrium solution.

Stability analysis
A first family aims at deriving a nonlinear eigenvalue problem assuming an ansatz exponential solution to obtain the characteristic equation. The delay introduces an exponential term within the characteristic equation and the corresponding delay eigenvalue problem becomes infinite-dimensional [13]. Dedicated solution methods often consider a rational or polynomial approximation of the characteristic equation in order to solve it such as the Solution Operator Discretization (SOD) and the Infinitesimal Generator Discretization (IGD) methods.
A second class attempts to build the monodromy matrixˆwhich connects the state of the system at time t i to the state at time t i , i.e. z.t i / Dˆz.t i /. The monodromy matrixˆis a finitedimensional approximation of the monodromy operator, which acts on the space of continuous functions and is infinite-dimensional. Different approaches have been developed such as, the Semi-Discretization Method (SDM), the Full Discretization Method (FDM), collocation techniques or temporal finite element analysis. From Floquet theory, the stability of the zero solution is addressed through the computation of the eigenvalues ofˆ. Different bifurcation mechanisms may be predicted such as Hopf, period-doubling or flip, and cyclic-fold bifurcations [2].

Semi-Discretization Method
The Semi-Discretization Method (SDM) is a recent development in the stability analysis of linear time-periodic DDE. It builds an approximation of the previously introduced monodromy operator and computes its eigenvalues. Floquet theory is then invoked for the stability analysis. SDM has been used for the analysis of a delayed Mathieu equation, as well as in turning and milling applications [14].
For a linear time-periodic 2n-dimensional DDE such as Eq. (5), the SDM is based on the discretization of the delayed term z.t / into a piece-wise constant function over one time period D Nt, while all the other time-dependent terms remain unchanged. For t sufficiently small (or N sufficiently large), the time-step interval is defined as t D OEt i I t i C1 , with t i 2 OE0 I and i D 0; 1 : : : .N 1/.
This way, 8t 2 OEt i I t iC1 , Eq. (5) may be approximated by: where z ;i is a weighted linear combination of the delayed vectors z.t i N / and z.t i N C1 /: The delayed term is thus assumed to be constant over OEt i I t i C1 yielding a finite-dimensional ODE approximation of (5) for which a general analytical solution is available in the form: where z.t i / D z i is a chosen initial condition, and the discrete linear coefficients are approximated by their mean value over the chosen interval as: Merging Eqs. (8) and (9) at t D t i C1 leads to: Accordingly, a discrete map can be built which defines a n.N C 1/ column state vector y i as: with: where P i;1Wn denotes the 1 to n columns of the matrix P i . Because the delay applies on displacements, only the first n columns of matrix R i are accounted for and the delayed terms in speed are not included in the discrete map. An approximation of the monodromy operatorˆcan then be built: The number of multiplications required to buildˆsolely depends on the time discretization N whereas the size of the system is dependent on both time and spatial discretizations. Finally, the stability analysis tests whether the eigenvalues modulii ofˆare less than 1 or not.

Modeling
In agreement with experimental and numerical investigations reported in [8] and [10] it is assumed that the interaction phenomenon of interest involves a single blade and the surrounding abradable coating.

Blade model
The blade of interest belongs to the last stage of the low-pressure compressor of an aircraft engine depicted in Fig. 3. The finite element mesh of the blade comprises 11,339 quadratic tetrahedron elements and 22,898 connecting nodes. Contact between blade and casing is treated on carefully selected boundary nodes on the blade tip. The number and location of these nodes (between the leading edge and the trailing edge) vary in this work depending on the contact configuration of interest. The size of the finite element model does not allow for acceptable computation times. Accordingly, it is reduced using a component mode synthesis method derived from the Craig-Bampton method that accounts for the linear contribution of centrifugal stiffening 3 [15]. The equation of motion is projected onto the reduced-order space as follows: The boundary displacement vector u b is defined in such a way that the targeted contact forces can directly be handled in the reduced space. The internal displacement vector u i is reduced to q modal participations of static and constraint modes computed at different rotational speeds (and stored in matrix ‰ c ) in order to account for centrifugal stiffening 4 . More details may be found in [9,10].

Cutting law in milling
During the cutting process in milling, when vibrations between cutter and workpiece are negligible, the evolution of the nominal chip thickness h is periodic and depends on the cutting angle and the feed per tooth f t . Since the cutting force F is often considered as proportional to the chip thickness, it is also time-periodic: where K e is the specific force coefficient and depends on the workpiece material and depth of cut [16]. When chatter occurs, either arising from regenerative waviness or mode-coupling, the uneven surface produced during the cut results in a variation of the instantaneous chip thickness, and the corresponding variations of the cutting force give rise to self-excited vibrations [16]. The dependence of the chip thickness h on the cutter's radial position u r -quantity directly extracted from the displacements vectors u.t/ and u.t / -may be written as: as illustrated in Fig. 4. Accounting for the vibrations around the teeth nominal paths leads to a system of coupled DDEs, and the typical state-space model used in milling is similar to the one defined in Eq. (5). It usually involves one or two degree-of-freedom (dof) models of the cutter and/or the workpiece. The stability analysis of these equations is used to maximize material removal rates without reaching instability zones [6].

DDEs for abradable material removal
In this section, only radial displacements of the blade boundary nodes are considered and stored in the vector u b . Depending on the number of nodes and their location, different contact configurations may be reproduced.
Similarly to the milling derivations, a cutting force proportional to the penetration now approximates the contact force arising from the incursion of the blade-tip in the abradable layer. This penetration is determined by the difference between the blade-tip displacement from one revolution to the succeeding one as shown in Fig. 5, while assuming that there is no separation between the blade and the abradable layer; the radial contact force f becomes:  3) and (4), where .t/ D becomes a constant matrix -due to the permanent contact assumption 5 -storing the cutting force coefficient K e applied on the selected boundary nodes.
A convenient use of modal coordinates u D Vx allows to have direct access to the contribution of a specific number of modes n and avoid the inversion of the mass matrix. This results in a projection of Eqs. (3) and (4) as: where the underlying linear system is then recast in the state space z T .t/ D x T .t/ P x T .t/ and yields a system of linear coupled DDE as Eq. (5), where A and B become constant matrices:

Preliminary results
Stability analysis with respect to K e , introduced in Eq. (18), is conducted within the interval K e 2 OE0 I 2:5 10 4 N/m so that the cutting forces are consistent with results reported in [10]. All frequencies in this section are normalized with respect to the first eigenfrequency of the blade at rest. The rotation speed range is 2 OE0 I 0:45.

Validation
Cutting conditions are only handled at the leading edge. The respective stability diagram is shown in Fig. 6 where the grey lobes display unstable configurations. This diagram determines whether the rest  (5) carried out with the MATLAB dde23 function, into which a constant radial displacement of the leading edge over the first revolution is prescribed as an initial condition: with u 0 b;LE D 2 10 3 m. Time integration is performed over 30 revolutions. As depicted in Fig. 7 both stable and unstable solutions picked from diagram 6 are accurately reproduced with time-integration of the DDE. In Fig. 7a, the solution diverges faster than in Fig. 7b since the maximum eigenvalue modulus of the respective monodromy matrix is higher for the explored configuration.

Convergence
In order to ensure accuracy, two types of convergence are assessed: time and modal convergence. The first type entails the convergence of the uniform mesh into which the time-period is discretized while the second focuses on the number of modes retained for the simulation. In these convergence studies, centrifugal stiffening effects are not accounted for and only contact at the leading edge is considered. The number N of time steps per period depends directly on the value of the delay D Nt, since the time-step t must be sufficiently small for the method to converge. Therefore, it is for lower rotational speeds -as D 1= -that a very large N may be required. Results pictured in Fig. 8 for increasing N highlight that for large rotational speeds . > 0:25/, N D 40 is sufficient to obtain accurate results. However, for lower rotational speeds, an accurate superimposition of the stability lobes is only achieved for N 60. Accordingly, N D 60 is kept for all subsequent simulations.
In addition, stability lobes are also dependent on the n modes retained in the projection of Eq. (19). Figure 9 underlines that the first eigenmode (first flexural mode) is not sufficient to suitably capture the blade dynamics within the entire targeted frequency range. Additional instability lobes appear when the second eigenmode (first torsion mode) is included. Also, considering more eigenmodes have little influence since the lobes obtained with five modes superimpose well with those obtained with two modes in Fig. 9. Therefore, only the first bending and torsional modes are retained for the remainder.

Centrifugal stiffening
Eigenfrequencies of the blade are affected by centrifugal stiffening and so are the corresponding instability lobes as depicted in Fig. 10. Interestingly, the stable domain gets larger at high rotational speeds. Also, centrifugal stiffening is substantially visible on the geometry and position of larger

Results
Time-domain simulation results are provided below for comparison purposes. They rely on previous developments [9] which are not reminded here for the sake of brevity. The abradable coating removal is modeled through a constitutive plastic law and contact conditions are handled on the trailing or/and leading edges to remain consistent with the DDE simulations. When blade/abradable coating contacts occur, the abradable elements -which are mono dimensional plastic bars -may undergo plastic deformations which are stored for the abradable profile update. At the end of the time simulations, the obtained abradable profile exhibits predominant contact areas. This comparison focuses on two critical aspects of usual interaction studies: 1. the identification of the most critical type of mode (usually bending or torsion) for a given blade design and a contact configuration, 2. the influence of the blade-tip/casing clearance distribution from leading to trailing edge to minimize the amplitude of vibration of the blade.

Sensitivity to blade-tip contact location
When a worst-case interaction scenario is identified on a turbomachine stage, it is usually dealt with by opening the clearance at this specific stage which alters the overall engine efficiency. Possible finer strategies would only locally adjust the clearance distribution since contact initiation on the trailing or leading edges may have distinct consequences. As an example, two ideal fictitious configurations are considered here: (1) contact only occurs on the leading edge boundary node and the trailing edge boundary node may vibrate with no restriction, and the symmetrical configuration (2) where contact only occurs on the trailing edge and the leading edge may vibrate with no restriction. The corresponding instability lobes are depicted in Fig. 11. This figure indicates that the unstable domain is significantly larger when contact occurs at the leading edge. Considering that -for the design of interest -the blade stiffness on the trailing edge is much higher than on the leading edge, these results hint that the stiffer the blade, the larger the stability domain. In practice, the blade stiffness is limited by aerodynamic and weight considerations thus leading to necessary design compromises.
In order to corroborate the results obtained with the DDE, time simulations from tool derived in [9], are carried out over the full rotational speed range and during 50 revolutions of the blade for each idealized configuration. Attention is paid to predicted profiles of removed abradable coating at the end of each simulation. Unfolded side by side, these profiles, directly related to the blade vibratory   Figure. 13 highlight a much more complex blade behavior. In this case, the abradable profile is highly dependent on the rotational speed of the blade. As an example, four main lobes are visible when ' 0:37 and seven lobes can be distinguished for ' 0:14. Dark areas in Fig. 13 underline that the worn lobes are deeper than in Fig. 12 which means that the blade amplitude of vibration is much higher in the present configuration. Results given by the proposed DDE strategy are confirmed by time simulations since the wear maps exhibit much higher levels of penetration within the abradable coating when contact occurs on the leading edge (Fig. 13).

Bending/torsion instability
As mentioned above, the first two eigenmodes of the blade (first bending and first torsion modes) allow for accurately capturing the instability of the trivial solution. In this section, the stability analysis is carried out retaining a single modal contribution. The two sets of instability lobes are superimposed in Fig. 14, where the bending mode instability lobes are wider and structurally more critical than the torsional lobes since they arise for smaller K e . Also, contrary to torsion instability, bending instability arises also at very low rotational speeds. Along the dashed blue line in Fig. 15, the position of the peaks matches very well the one of the instability lobes associated with the first bending mode. The main lobe of instability associated with the first bending mode of the blade predicted by the DDE approach clearly lines up with the resonance of the same mode in the time integration.
Overall, results provided with DDE show that the blade is more sensitive to bending instability than to torsion instability. This is confirmed by the higher peaks detected around the first bending mode in the frequency domain analysis of the time simulations, and is consistent with previous work [10], where the dominance of the first bending mode during a blade/abradable coating interaction was both experimentally and numerically observed.
The resonances predicted by the simulations in the time-domain are very consistent with the instability lobes predicted by the DDE despite of the fact that the latter lies on a significantly simplified assumption of the interaction phenomenon. The fact that only two modes (first bending and first torsion modes) are sufficient for an accurate description of the stability is also confirmed with the frequency domain analysis of time simulations results since no significant peaks of vibration amplitude are detected for vibratory frequencies f > 5.

Conclusion
This study focuses on a qualitative analysis of abradable coating removed by an incurring blade within turbomachines based on delay differential equations. Preliminary results highlight the convergence of the method and the accuracy of the instability lobes detected even when a relatively small number of blade modes is considered. The provided comparison shows consistency with simulations previously developed in the time domain. The proposed strategy accurately predicts peaks of resonance of the blade for both the first bending and torsion modes. Also, the sensitivity of stability areas to the contact locations on the tip of the blade is addressed and is in agreement with time simulations.
The proposed approach advantageously involves reduced computation times. It could be included in the blade design stage to account for possible interaction criterion that would involve a minimal specific force coefficient K e under which instability is unacceptable. Designers may take advantage of the presented strategy to discriminate blade designs poorly robust to the interaction.
Among future roads of investigation, work is in progress to assess (1) the influence of friction in the interaction phenomenon, (2) the effects of blade/abradable separation events and (3) the coupling with tangential tip-displacements yielding state-dependent DDE. Also, the implementation of a flexible casing through coupled DDE is under investigation.