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

Nicolas Salvat, Alain Batailly, Mathias Legrand

► To cite this version:

Nicolas Salvat, Alain Batailly, Mathias Legrand. Modeling of abreadable coating removal in aircraft engines through delay differential equations. Journal of Engineering for Gas Turbines and Power, 2013, 135 (10), 10.1115/1.4024959 . hal-00976915v2

HAL Id: hal-00976915

https://hal.science/hal-00976915v2

Submitted on 22 Apr 2014

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modélisation de l'enlèvement de matériau abradable par équations différentielles à retard

<figure><img src="image_2.png" /></figure>

## Nicolas Salvat

Laboratoire de Dynamique des Structures et Vibrations, Département de Génie Mécanique, Université McGill, 817 rue Sherbrooke Ouest, Montréal, Québec H3A 0C3

## Alain Batailly

Laboratoire de Dynamique des Structures et Vibrations, Département de Génie Mécanique, Université McGill, 817 rue Sherbrooke Ouest, Montréal, Québec H3A 0C3

## Mathias Legrand

Laboratoire de Dynamique des Structures et Vibrations, Département de Génie Mécanique, Université McGill, 817 rue Sherbrooke Ouest, Montréal, Québec H3A 0C3

## Résumé

L'augmentation du rendement des turbomachines aéronautiques conduit les concepteurs à privilégier deux stratégies distinctes : (1) l'utilisation de matériaux plus légers tels que les matériaux composites et (2) la réduction des jeux entre parties tournantes (telles que les roues aubagées) et le carter environnant. La réduction des jeux a notamment pour effet de favoriser l'apparition de contacts structuraux qui sont dorénavant considérés comme des événements normaux, bien qu'indésirables, pouvant se produire au cours du fonctionnement de la turbomachine. Afin de limiter les conséquences de ces contacts sur l'intégrité structurelle des composants de la turbomachine, il est fréquent d'utiliser un revêtement sacrificiel, appelé revêtement abradable, qui est disposé sur le carter et qui doit satisfaire le compromis suivant : être suffisamment souple pour ne pas endommager les aubes lors de contacts mais également suffisamment résistant à son environnement (haute température, corrosion...).

Longtemps vus comme une solution robuste pour empêcher l'apparition de hauts niveaux de vibrations lors de contacts structuraux, les matériaux abradables n'empêchent toutefois pas certains phénomènes d'interaction de se produire. Des observations expérimentales et numérique ont notamment mis en évidence des phénomènes d'interaction de type aube/abradable pouvant conduire à une perte d'aube. Du fait de l'absence de données sur l'origine de ces phénomènes, les premières études réalisées ont favorisé des stratégies avec intégration temporelle. Toutefois des limitations, en termes de temps de calculs notamment, conduisent naturellement à la recherche d'une stratégie d'étude qualitative afin d'obtenir une compréhension plus globale de ces phénomènes d'interaction.

Il est possible de faire une analogie entre le procédé de fraisage et les contacts structuraux entre aubes et revêtement abradable environnant, car ces deux phénomènes d'interaction—pièce/outil dans le premier cas et aube/abradable dans le deuxième—se traduisent par un enlèvement de matière comme illustré par la figure 1.

Dans le domaine de l'usinage à enlèvement de matière, tel que le tournage et le fraisage, il est primordial de limiter les vibrations entre outil et pièce, afin d'obtenir des états de surface convenables et de réduire l'usure de l'outil. Dans ces procédés, un phénomène connu sous le nom de *broutement* produit des vibrations auto-entretenues du fait des surfaces ondulées générées lors de la coupe.
<figure><img src="image_3.png" /><figcaption>**Figure 1:** Analogie entre le fraisage et l'enlèvement d'abradable dans une turbomachine</figcaption></figure>

Lors de la modélisation du broutement, considéré comme particulièrement dangereux et indésirable, la nature régénérative du phénomène est souvent prise en compte à l'aide d'équations différentielles à retard (EDR), i.e. l'état du système dépend de son état courant ainsi que d'états précédents. Dans le cas du fraisage (c.f. Fig. 1), l'épaisseur du copeau $h(t)$, déterminée à partir de l'avance par dent $f_t$ et de l'angle de coupe $\phi$, devient dépendante de la position de l'outil à l'état courant $u(t)$ et de sa position lorsque la dent précédente était en coupe $u(t - \tau)$. L'analyse des EDR issues de cette modélisation permet de faire un choix de paramètres de coupe maximisant les taux d'enlèvement de matière sans atteindre des zones d'instabilités

Ainsi, il est proposé de modéliser l'enlèvement de matériau abradable, du fait de contacts aubes/ abradable, par des équations différentielles à retard (c.f. Fig. 1). Ces équations se distinguent des équations différentielles ordinaires en ce sens que l'évolution du système à un instant $t$ dépend de certains états précédents du système à un instant $t - \tau_i$. Les termes de retard peuvent être : des constantes positives, dépendants du temps ou de l'état du système, ou encore des variables continues. Dans le cadre de cette étude, seul le cas avec un terme de retard en déplacement constant $\tau$ sera considéré :

$$ \dot{z}(t) = f(t, z(t), z(t - \tau)) \qquad (1) $$

où $z = [u \quad \dot{u}]^T$ est le vecteur d'état inconnu et la constante $\tau$ est le retard à prendre en compte.

Dans le cas du fraisage et du tournage, du fait de l'avance, il est possible de considérer que les deux parties sont en contact permanent, condition permettant de réaliser des études de stabilité. Seules quelques études permettent la gestion du décollement pièce/outil. Dans le cadre de ce modèle d'interaction aube/abradable, il est proposé de prendre en compte le décollement aube/abradable, en plus des analyses de stabilité classiques, et de réaliser une étude de bifurcations des EDR construites afin d'identifier des régimes de fonctionnement potentiellement dangereux.
L'étude se décompose en trois parties: (1) dans un premier temps le modèle avec équations différentielles à retard est détaillé, introduisant notamment la possibilité de prise en compte du décollement. (2) La stratégie proposée est ensuite validée par comparaison avec des résultats obtenus par intégration temporelle. Enfin, (3) une étude paramétrique est réalisée afin de mettre en évidence les applications potentielles de la méthode.

Les premiers résultats obtenus avec la méthode proposée sont présentés sous forme de cartes de stabilité mettant en avant des vitesse de rotation critiques. Il est notamment mis en évidence que les deux premiers modes propres—premier mode de flexion et premier mode de torsion—jouent un rôle fondamental dans la définition des zones d’instabilité. Des instabilités en torsion ou en flexion sont distinguées suivant la vitesse de rotation considérée. Par ailleurs, l’influence de la localisation du contact—au bord de fuite ou au bord d’attaque du sommet d’aube—a également été étudiée, celle-ci semble grandement affecter le comportement de la structure. La cohérence de ces cartes avec des résultats obtenus par un modèle d’intégration temporelle, renforce la validité de la stratégie avec équations différentielles à retard.
# Modeling of abradable coating removal in aircraft engines through delay differential equations

<figure><img src="image_4.png" /></figure>

**Nicolas Salvat**

Structural Dynamics and Vibration Laboratory, Department of Mechanical Engineering, McGill
University, 817 Sherbrooke St West, Montréal, Québec H3A 2K6, Canada
e-mail: nicolas.salvat@mail.mcgill.ca

**Alain Batailly**

Structural Dynamics and Vibration Laboratory, Department of Mechanical Engineering, McGill
University, 817 Sherbrooke St West, Montréal, Québec H3A 2K6, Canada

**Mathias Legrand**

Structural Dynamics and Vibration Laboratory, Department of Mechanical Engineering, McGill
University, 817 Sherbrooke St West, Montréal, Québec H3A 2K6, Canada

## 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.
Nomenclature

Γ external force matrix

Λ eigenvalue matrix

Φ monodromy matrix

Ψ<sub>c</sub> generalized constraint modes

Ψ<sub>s</sub>(0) static modes computed at rest

0<sub>n</sub> null matrix of dimension n×n

C damping matrix

D time discretization matrix

I<sub>n</sub> identity matrix of dimension n×n

K stiffness matrix

M mass matrix

V eigenvectors matrix

f external force vector

h abradable coating thickness

u displacement vector

x modal coordinate vector

z state-space vector

Ω rotational speed

τ time delay

ξ modal damping coefficient

f frequency

K<sub>e</sub> specific force coefficient (N·m<sup>-1</sup>)

N time discretization parameter

t time
# 1 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:

$$ \dot{z}(t) = g(t, z(t), z(t - \tau)) \qquad (2) $$

where $t$ is the time, $\tau$ is the so-called time delay and $\dot{z}$ denotes the derivative of $z$ with respect to time. Accordingly, initial conditions of a DDE are a function defined over the interval $[-\tau; 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:

$$ \ddot{M}(t) + C\dot{u}(t) + Ku(t) = f(t, u(t), u(t - \tau)) \qquad (3) $$

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

<figure><img src="image_5.png" /><figcaption>Figure 2: Analogy between abradable removal and milling</figcaption></figure>

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-incurring 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<fnref n="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 $\dot{z}(t)$ not only depends on its current state $z(t)$, but also on a past state $z(t - \tau)$ as specified in Eq. (2). In the particular case of mechanical systems, the forcing term $f$ in Eq. (3) is often written as:

$$f(t) = \Gamma(t)(u(t - \tau) - u(t)) \qquad (4)$$

where $\Gamma(t)$ is an external forcing matrix which may be time dependent. By introducing state-space coordinates $z^T(t) = (u^T(t) \quad \dot{u}^T(t))$ in order to obtain a first order DDE, Eqs. (3) and (4) yield:

$$\begin{aligned} \dot{z}(t) &= A(t)z(t) + B(t)z(t - \tau) \\ z(t) &= \varphi(t), \quad t \in [-\tau; 0] \end{aligned} \qquad (5)$$

with,

$$\begin{aligned} A(t) &= \begin{bmatrix} 0_n & I_n \\ -M^{-1}(K + \Gamma(t)) & -M^{-1}C \end{bmatrix} \\ B(t) &= \begin{bmatrix} 0_n & 0_n \\ M^{-1}\Gamma(t) & 0_n \end{bmatrix} \end{aligned} \qquad (6)$$

<footnote n="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].</footnote>
where, $z(t)$ is the $2n$-dimensional state-space vector, $A(t)$ and $B(t)^2$ are $2n \times 2n$ matrices dependent on the system parameters and external forcing, and $\varphi$ is the initial condition function defined over the interval $[-\tau; 0]$.

Two cases are of particular interest: (1) DDE with constant coefficients, where $A(t) = A$ and $B(t) = B$ are constant matrices, and (2) DDE with time-periodic coefficients of period $T$, where $A(t+T) = A(t)$ and $B(t+T) = B(t)$. In the present paper, both cases are considered, limiting the theoretical developments to $T = \tau$.

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 = 0$ solution, also referred to as equilibrium solution.

## 2.2 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 $\Phi$ which connects the state of the system at time $t_i$ to the state at time $t_i - \tau$, i.e. $z(t_i) = \Phi z(t_i - \tau)$. The monodromy matrix $\Phi$ is a finite-dimensional 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 $\lambda$ of $\Phi$. Different bifurcation mechanisms may be predicted such as Hopf, period-doubling or flip, and cyclic-fold bifurcations [2].

## 2.3 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 - \tau)$ into a piece-wise constant function over one time period $\tau = N\Delta t$, while all the other time-dependent terms remain unchanged. For $\Delta t$ sufficiently small (or $N$ sufficiently large), the time-step interval is defined as $\Delta t = [t_i; t_{i+1}]$, with $t_i \in [0; \tau]$ and $i = 0, 1, \dots (N-1)$.

This way, $\forall t \in [t_i; t_{i+1}]$, Eq. (5) may be approximated by:

$$ \dot{z}(t) = A(t)z(t) + B(t)z_{\tau,i} \qquad (7) $$

<footnote n="2">²Note that the last $n$ columns of $B(t)$ are zeros, accordingly, the resulting DDE is of retarded type.</footnote>
where $z_{\tau,i}$ is a weighted linear combination of the delayed vectors $z(t_{i-N})$ and $z(t_{i-N+1})$:

$$z_{\tau,i} = \frac{1}{2}z_{i-N+1} + \frac{1}{2}z_{i-N} \approx z(t - \tau) \qquad (8)$$

The delayed term is thus assumed to be constant over $[t_i; t_{i+1}]$ yielding a finite-dimensional ODE approximation of (5) for which a general analytical solution is available in the form:

$$z(t) = \exp(A_i(t - t_i))(z_i + A_i^{-1}B_i z_{\tau,i}) - A_i^{-1}B_i z_{\tau,i} \qquad (9)$$

where $z(t_i) = z_i$ is a chosen initial condition, and the discrete linear coefficients are approximated by their mean value over the chosen interval as:

$$A_i = \frac{1}{\Delta t} \int_{t_i}^{t_{i+1}} A(t)dt \quad \text{and} \quad B_i = \frac{1}{\Delta t} \int_{t_i}^{t_{i+1}} B(t)dt \qquad (10)$$

Merging Eqs. (8) and (9) at $t = t_{i+1}$ leads to:

$$z_{i+1} = P_i z_i + \frac{1}{2} R_i z_{i-N+1} + \frac{1}{2} R_i z_{i-N} \qquad (11)$$

with $P_i = \exp(A_i \Delta t)$ and $R_i = (\exp(A_i \Delta t) - I_{2n})A_i^{-1}B_i$. Accordingly, a discrete map can be built which defines a $n(N+1)$ column state vector $y_i$ as:

$$y_i = [u_i \quad \dot{u}_i \quad u_{i-1} \quad u_{i-2} \quad \dots \quad u_{i-N}]^T \quad \text{and} \quad y_{i+1} = D_i y_i \qquad (12)$$

with:

$$D_i = \begin{bmatrix} P_{i,1:n} & P_{i,n+1:2n} & 0_n & \cdots & 0_n & \frac{1}{2}R_{i,1:n} & \frac{1}{2}R_{i,1:n} \\ I_n & 0_n & 0_n & \cdots & 0_n & 0_n & 0_n \\ 0_n & 0_n & I_n & & 0_n & 0_n & 0_n \\ \vdots & & \ddots & & \ddots & & \vdots \\ 0_n & 0_n & & I_n & 0_n & 0_n & 0_n \\ 0_n & 0_n & \cdots & 0_n & I_n & 0_n & 0_n \end{bmatrix} \qquad (13)$$

where $P_{i,1:n}$ 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 $\Phi$ can then be built:

$$y_N = \Phi y_0 \quad \text{with} \quad \Phi = D_{N-1}D_{N-2}\cdots D_1D_0 \qquad (14)$$

The number of multiplications required to build $\Phi$ 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 $\Phi$ are less than 1 or not.

### 3 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.
## 3.1 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.

<figure><img src="image_6.png" /><figcaption>Figure 3: Aircraft engine with blade of interest: (1) fan, (2) low pressure compressor.</figcaption></figure>

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<fnref n="3" /> [15]. The equation of motion is projected onto the reduced-order space as follows:

$$ \mathbf{u} = \begin{pmatrix} \mathbf{u}_b \\ \mathbf{u}_i \end{pmatrix} = \begin{bmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{\Psi}_s(0) & \mathbf{\Psi}_c \end{bmatrix} \begin{pmatrix} \mathbf{u}_b \\ \mathbf{q} \end{pmatrix} \quad (15) $$

The boundary displacement vector $\mathbf{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 $\mathbf{u}_i$ is reduced to $\mathbf{q}$ modal participations of static and constraint modes computed at different rotational speeds (and stored in matrix $\mathbf{\Psi}_c$) in order to account for centrifugal stiffening<fnref n="4" />. More details may be found in [9, 10].

## 3.2 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 $\phi$ and the feed per tooth $f_t$. Since the cutting force $F$ is often considered as proportional to the chip thickness, it is

<footnote n="3">³The blade stiffness matrix is dependent on the rotational speed: $K = K(\Omega)$.</footnote>

<footnote n="4">⁴For conciseness, structural matrices — such as the mass and stiffness matrices $M$ and $K$ — are always the reduced version obtained with the aforementioned component mode synthesis method.</footnote>
also time-periodic:

$$F(t) = K_e h(t) \quad \text{with} \quad h(t) = f_t \sin(\phi(t)) \qquad (16)$$

where $K_e$ is the specific force coefficient and depends on the workpiece material and depth of cut [16].

<figure><img src="image_7.png" /><figcaption>Figure 4: Schematic of milling with chatter</figcaption></figure>

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 $\mathbf{u}(t)$ and $\mathbf{u}(t - \tau)$ — may be written as:

$$h(t) = f_t \sin(\phi(t)) + u_r(t - \tau) - u_r(t) \qquad (17)$$

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].

### 3.3 DDEs for abradable material removal

In this section, only radial displacements of the blade boundary nodes are considered and stored in the vector $\mathbf{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:

$$\mathbf{f}(\mathbf{h}(t)) = K_e \mathbf{h}(t) = K_e [\mathbf{u}_b(t - \tau) - \mathbf{u}_b(t)] \qquad (18)$$
<figure><img src="image_8.png" /><figcaption>Figure 5: Schematic of abradable material removal</figcaption></figure>

where the thickness of the removed abradable material is stored in a vector h, dependent on the displacement of the boundary nodes u_b from one revolution to the next and the contact configuration modeled.

Therefore, the governing equation of motion are of the form of Eqs. (3) and (4), where $\Gamma(t) = \Gamma$ becomes a constant matrix — due to the permanent contact assumption<fnref n="5" /> — storing the cutting force coefficient $K_e$ applied on the selected boundary nodes.

A convenient use of modal coordinates $\mathbf{u} = \mathbf{V}\mathbf{x}$ 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:

$$ \ddot{\mathbf{x}}(t) + 2\xi\mathbf{\Lambda}\dot{\mathbf{x}}(t) + \mathbf{\Lambda}^2\mathbf{x}(t) = \mathbf{V}^T\mathbf{\Gamma}\mathbf{V}[\mathbf{x}(t-\tau) - \mathbf{x}(t)] \quad (19) $$

where the underlying linear system is then recast in the state space $\mathbf{z}^T(t) = (\mathbf{x}^T(t) \quad \dot{\mathbf{x}}^T(t))$ and yields a system of linear coupled DDE as Eq. (5), where $\mathbf{A}$ and $\mathbf{B}$ become constant matrices:

$$ \mathbf{A} = \begin{bmatrix} \mathbf{0}_n & \mathbf{I}_n \\ -\mathbf{\Lambda}^2 - \mathbf{V}^T\mathbf{\Gamma}\mathbf{V} & -2\xi\mathbf{\Lambda} \end{bmatrix} \quad \text{and} \quad \mathbf{B} = \begin{bmatrix} \mathbf{0}_n & \mathbf{0}_n \\ \mathbf{V}^T\mathbf{\Gamma}\mathbf{V} & \mathbf{0}_n \end{bmatrix} \quad (20) $$

# 4 Preliminary results

Stability analysis with respect to $K_e$, introduced in Eq. (18), is conducted within the interval $K_e \in [0; 2.5] \times 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 $\Omega \in [0; 0.45]$.

## 4.1 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

<footnote n="5">⁵In addition to this assumption, unlike machining processes, the abradable coating removal in turbomachinery does not involve a feed.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 6: Test points: • stable domain ($K_e = 2 \cdot 10^3$ N/m); • unstable domain ($K_e = 2 \cdot 10^4$ N/m)</figcaption></figure>

position $z(t) = 0$ of the blade is stable or not. These results are compared to direct time integration of the DDE (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:

$$ \forall t \in [-\tau; 0], \quad z^T(t) = (u^T(t) \quad \dot{u}^T(t)) = (u_{b,LE}^0 \quad 0 \quad \dots \quad 0) \qquad (21) $$

with $u_{b,LE}^0 = 2 \cdot 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.

<figure><img src="image_10.png" /><figcaption>Figure 7: Tip radial displacement; $K_e = 2 \cdot 10^3$ N/m (—) and $K_e = 2 \cdot 10^4$ N/m (—)</figcaption></figure>
## 4.2 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 τ = *N*Δ*t*, since the time-step Δ*t* must be sufficiently small for the method to converge. Therefore, it is for lower rotational speeds Ω — as τ = 1/Ω — that a very large *N* may be required.

<figure><img src="image_11.png" /><figcaption>**Figure 8:** Time convergence: *N* = 40 (—); *N* = 60 (---), and *N* = 80 (---)</figcaption></figure>

Results pictured in Fig. 8 for increasing *N* highlight that for large rotational speeds (Ω > 0.25), *N* = 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* = 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.

## 4.3 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
<figure><img src="image_12.png" /><figcaption>**Figure 9:** Modal convergence: mode 1 (—); modes 1 and 2 (---); and modes 1 to 5 (—-)</figcaption></figure>

lobes while thinner ones are essentially unmodified. Centrifugal stiffening is accounted for in all the simulations presented in the following.

<figure><img src="image_13.png" /><figcaption>**Figure 10:** Instability lobes with (---) and without (—) centrifugal effects</figcaption></figure>

# 5 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.

**5.1 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

<figure><img src="image_14.png" /><figcaption>**Figure 11:** Instability lobes versus blade-tip contact location: trailing edge (—) and leading edge (---)</figcaption></figure>

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><img src="image_15.png" /><figcaption>**Figure 12:** Abradable wear map: contact at blade trailing edge</figcaption></figure>

amplitude, create a wear map where darker points indicate higher levels of removal. In order to initiate the blade/abradable coating contact the casing is ovalized, which results in two main grooves — for $\theta \in [1.4; 2]$ and $\theta \in [4.4; 5]$ — observed in the abradable profiles over the full rotational speed range in Figs. 12 and 13.

Wear maps depicted in Figs. 12 and 13 are obtained when contact occurs only at the trailing edge and leading edge respectively (both figures have the same color code for comparison). On the trailing edge, the abradable profile is uniform throughout the rotational speed range: it features two main lobes as shown in Fig. 12.

<figure><img src="image_16.png" /><figcaption>**Figure 13:** Abradable wear map: contact at blade leading edge</figcaption></figure>

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 $\Omega \simeq 0.37$ and seven lobes can be distinguished for $\Omega \simeq 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).

## 5.2 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.

<figure><img src="image_17.png" /><figcaption>**Figure 14:** Instability lobes: bending (---) and torsion (—) modes</figcaption></figure>

Figure. 15 shows the spectrum obtained from the time simulations results, associated with the radial displacement of the leading edge for each rotational speed. Two sets of peaks are distinguishable, around the natural frequencies of both the first bending and first torsion modes. Time simulations carried out in this section account for possible contacts both on the trailing edge and the leading edge of the blade.

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],
<figure><img src="image_18.png" /><figcaption>**Figure 15:** Spectrum of time-domain response: first torsion mode (—) and first bending mode (---); instability lobes from Fig. 14 are superimposed</figcaption></figure>

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$.

# 6 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.

Acknowledgments

Thanks go to Snecma for its technical and financial support. This work takes place in the framework of
the MAIA mechanical research and technology program sponsored by CNRS, ONERA and SAFRAN
Group.

References

[1] T. Erneux. *Applied delay differential equations*, volume 3 of *Surveys and Tutorials in the Applied Mathematical Sciences Series*. Springer-Verlag, 2009.

[2] B. Balachandran, T. Kalmar-Nagy, and D. Gilsinn. *Delay differential equations: Recent advances and new directions*. Springer Verlag, 2009.

[3] W. Wang and L. Chen. A predator-prey system with stage-structure for predator. *Computers Mathematics with Applications*, 33(8):83–91, 1997.

[4] M. Roussel. The use of delay differential equations in chemical kinetics. *Journal of Physical Chemistry*, 100(20):8323–8330, 1996.

[5] J. Cook and B. Powell. Modeling of an internal combustion engine for control analysis. *Control Systems Magazine*, 8(4):20–26, 1988.

[6] G. Stépán, R. Szalai, B. Mann, P. Bayly, T. Insperger, J. Gradisek, and E. Govekar. Nonlinear dynamics of high-speed milling—analyses, numerics, and experiments. *Journal of Vibration and Acoustics*, 127(2):197–203, 2005.

[7] T. Insperger, G. Stépán, and J. Turi. State-dependent delay in regenerative turning processes. *Nonlinear Dynamics*, 47(1):275–283, 2007.

[8] A. Millecamps, J.-F. Brunel, P. Dufrénoy, F. Garcin, and M. Nucci. Influence of thermal effects during blade-casing contact experiments. In *ASME Conference Proceedings: 22nd Biennial Conference on Mechanical Vibration and Noise*, number 48982, pages 855–862, San Diego, USA, 2009.

[9] M. Legrand, A. Batailly, and C. Pierre. Numerical investigation of abradable coating removal through plastic constitutive law in aircraft engine. *Journal of Computational and Nonlinear Dynamics*, 7, 2012.
[10] A. Batailly, M. Legrand, A. Millecamps, and F. Garcin. Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact. *Journal of Engineering for Gas Turbines and Power*, 134(8), 2012.

[11] R. Williams. Simulation of blade casing interaction phenomena in gas turbines resulting from heavy tip rubs using an implicit time marching method. In *ASME Conference Proceedings: Turbo Expo 2011*, number 45495, Vancouver, Canada, 2011.

[12] Z. Dombovari, D. Barton, R. Wilson, and G. Stépán. On the global dynamics of chatter in the orthogonal cutting model. *International Journal of Non-Linear Mechanics*, 46(1):330–338, 2011.

[13] E. Jarlebring. *The spectrum of delay-differential equations: numerical methods, stability and perturbation*. PhD thesis, Technische Universität Carolo-Wilhelmina zu Braunschweig, Vännäs, Schweden, 2008.

[14] T. Insperger and G. Stépán. *Semi-discretization for time-delay systems: stability and engineering applications*. Springer Science, 2011.

[15] A. Sternchüss and E. Balmès. On the reduction of quasi-cyclic disks with variable rotation speeds. *Proceedings of the International Conference on Advanced Acoustics and Vibration Engineering (ISMA)*, pages 3925–3939, 2006.

[16] J. Tlusty. *Manufacturing processes and equipment*. Prentice Hall, 2000.