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

HAL
open science

# Phenomenological model for stability analysis of bladed rotor-to-stator contacts

Marie-Océane Parent, Fabrice Thouverez

► To cite this version:

Marie-Océane Parent, Fabrice Thouverez. Phenomenological model for stability analysis of bladed rotor-to-stator contacts. International Symposium on Transport Phenomena and Dynamics of Rotating Machinery, Apr 2016, Honolulu, United States. hal-01537643

HAL Id: hal-01537643

https://hal.science/hal-01537643

Submitted on 12 Jun 2017

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Phenomenological model for stability analysis of bladed rotor-to-stator contacts

Marie-Océane Parent<fnref n="1" />*, Fabrice Thouverez<fnref n="2" />

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

ISROMAC 2016

International Symposium on Transport

Phenomena and Dynamics of

Rotating Machinery

Hawaii, Honolulu

April 10-15, 2016

**Abstract**

The present study aims to predict the dynamic stability of a turbofan engine submitted to light contacts between fan blades tips and their surrounding casing. The study is based on a hybrid model which introduces a simplified bladed wheel and a flexible suspended casing to a rotor-shaft model. Inertial coupling is considered between blades and shaft dofs and an elastic link is introduced between the shaft and the suspension dofs of the casing. A 3D contact formulation has been implemented to recover accuracy in the contacts detection and the application of contact forces; it considers the 3D model kinematic and introduces the local geometry of the contact area. The model behavior with blade-to-casing contacts is analyzed through two approaches: the first one assumes permanent contacts while the other one allows intermittent contacts. Good agreement is obtained between the two types of results. The results highlight the importance of couplings in the initiation of unstable phenomena and the relevance of the 3D contact formulation in predicting the stability of the system.

**Keywords**

Rotordynamics — Bladed assembly — Blade-tip/casing contacts

## INTRODUCTION

By limitating the leakage flows, reducing the clearance between blade tips and surrounding casing improves the performance and operability of turbomachines. The aero-engine designers would thus reduce theses distances at their minimum, but the counterpart lays in an increase of the possible contacts between rotating and fixed parts, which can cause high amplitude vibrations or even lead to unstable dynamic behaviours. The phenomena related to these contacts represent an important research field which aims to provide the numerical tools to design high performance turbomachine while ensuring to avoid dangerous dynamic behaviour.

Contacts between rotating and fixed parts can be induced by many phenomena: vibration levels associated with unbalance (residual or accidental), aircraft maneuver, or casing thermal deformation. Furthermore, they may appear in various locations and involve very different organs of the engine. This variety of contact configurations leads to a wide spectrum of phenomena. Depending on the contact configuration and the studied phenomena, most of the analyses splits into two areas of study:

Rotor-to-stator contact (or rubbing) focuses on the contacts which do not involve the bladed assembly dynamics. The analyses are based on rotordynamics, they includes a rich phenomenology described and investigated by numerous authors. The reviews presented in [1, 2, 3] summarize the contributions made both numerically and experimentally.

The dynamic behavior resulting from blade-to-casing contacts is subject to a more recent literature. The traveling wave

speed coincidence was first studied by Schmiechen [4]. This phenomenon consists of a modal interaction between two rotating modes: one on the bladed wheel and the other on the casing [5, 6, 7]. The models used on these investigations take into account the bladed disks and casing flexibility while neglecting the shaft dynamics.

In this paper, we focus on the engine dynamic stability when the response to a residual unbalance lead to small touches of the fan blade tips on their surrounding casing (see Figure 1 (a)). The study proposed in this paper involve an hybrid model including some aspects of the two approaches presented above: blades and casing flexibility are considered as well as the Low Pressure (LP) shaft dynamics. Sinha [8] and Lesaffre [9] developed flexible, bladed shaft models for this type of blade tip-casing contact investigation. The study presented in this paper follows their approaches.

The model developed in [9] has been completed to introduce rigid-body displacements of the casing and the coupling with the fan frame struts and the first bearings (see Figure 1). A 3D contact formulation has also been implemented to improve the accuracy in the touches detection and the application of contact forces.

Results obtained with the same model architecture, but with a simpler contact formulation, were presented in [10]. They highlighted the effect of shaft-casing and blades-shaft coupling on the system behaviour during contacts. Another proceeding [11] presented the advantage of the 3D contact formulation regarding touches detection and stability results. In the present paper, the model formulation is fully described with all parameters given. The compliance between the per-

<footnote n="1">¹Snecma, SAFRAN Group, Site de Villaroche, France</footnote>

<footnote n="2">²LTDS, École Centrale de Lyon, France</footnote>

<footnote>*Corresponding author: marieoceane.parent@snecma.fr</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 1. (a) schematic section view of turbofan engine, (b) section view of the model</figcaption></figure>

manent contact results and the transient results is clarified and unusual behaviour (precession orientation switches) is described.

# 1. NUMERICAL MODEL

## 1.1 Model formulation

The bladed rotor model used in the present paper derives from the model built by Sinha [8]. It consists of a flexible shaft connected to a rigid disk to which a set of flexible blades are clamped. Lesaffre et al. [12] formulated the bladed shaft model in the rotating frame by the energetic method. In [9], a flexible casing modeled by an non-extensible elastic ring completed the model to carry out stability analysis considering blade-tip to casing contacts.

The models built in [8, 9] consider the LP shaft dynamics but limits its kinematics outside from the fan area: the shaft has indeed a clamped boundary condition at its rear cross-section. In the present model, the rear part of the rotor is supported by an isotropic suspension representative of the turbine rear frame struts and connected to a rigid disk representative of the LP turbine. In addition, the rigid body displacement of the casing are considered and an elastic link between these suspension dofs and the shaft cross-section is introduced to describe the existing mechanical link going through the bearings and the fan frame (See Figure 1).

To eliminate time-dependent terms resulting from the periodicity of the bladed disk, the entire model has been developed in the rotating frame under the assumption of a constant spin speed (denoted $\Omega$). All parameters definitions and values are presented in appendix.

### 1.1.1 Flexible bladed rotor

The shaft allows two displacements defined in the rotating frame: the two orthogonal translations in the cross-section plane, denoted $x_s(z, t)$ and $y_s(z, t)$. Two rigid disks con-

<figure><img src="image_5.png" /><figcaption>Figure 2. 3D visualisation of the model</figcaption></figure>

nected to the shaft are modeled by their mass and rotational inertia. Lastly, a set of ten flexible blades are clamped to the representative fan disk; they are modeled by Euler-Bernoulli beams with a rectangular cross-section. The blades can only deflect along their more flexible direction; the bending at the coordinate $s$ along the $j^{th}$ blade is denoted: $x_{b_j}(s, t)$. The shaft and blade deformations are discretized through use of Rayleigh-Ritz functions. As in [12] and [17], the system of equation is build through use of the energetic method. The Rayleigh-Ritz functions and all energies expressions are given in appendix.

### 1.1.2 Casing

The casing being fixed, we choose to develop its formulation in the inertial frame and then transpose in the rotating frame.

The casing suspension provides four rigid body dofs: the two translations and two rotations in the shaft cross-section plane. They are denoted $X_c, Y_c, \Phi_{X_c}, \Phi_{Y_c}$ in the inertial frame. The casing model consists of an elastic ring with radial and tangential displacements $(u, w)$. The rigid-body rotations and the deformation orientations are sketched in figure 3. From the inertial frame ($X, Y, Z$) whose origin coincide with the casing gravity center, the rotations $\Phi_{X_c}$ and $\Phi_{Y_c}$ create successively the frames $(X_1, Y_1, Z_1)$ and $(X_2, Y_2, Z_2)$. In this last frame, the displacements at the angular position $\theta$ are defined by $u$ et $w$.

The tangential displacement at $\theta$ is given by [13, 14]:

$$w(\theta, t) = \sum_{n_d} A_{n_d}(t) \cos(n_d \theta) + B_{n_d}(t) \sin(n_d \theta) \quad (1)$$

The inextensible property of the elastic ring implies a relation between the radial and tangential displacement:

$$u(\theta, t) = \frac{\partial w(\theta, t)}{\partial \theta} \quad (2)$$

The two and three nodal diameter deformations are introduced ($n_d = 2, 3$). The zero nodal diameter is indeed not compatible with the inextensibility property and the one nodal diameter mode is brought by the rigid body displacements $X_c(t)$ et $Y_c(t)$.

The casing kinetic energy is given by:

$$\mathcal{T}_c = \frac{1}{2} \int_0^{2\pi} \rho_c S_c \mathbf{V}_c^T \mathbf{V}_c R_c d\theta + \frac{1}{2} \int_0^{2\pi} \rho_c \mathbf{\Omega}_c^T \mathbf{I}_c \mathbf{\Omega}_c R_c d\theta$$
<figure><img src="image_6.png" /><figcaption>Figure 3. Casing rotations and frames definition</figcaption></figure>

(3)

$V_c$ and $\Omega_c$ are the instantaneous speeds at the center of gravity of the $\theta$ cross-section:

$$ \begin{cases} \mathbf{V}_c = \frac{d}{dt}(X_c(t)\mathbf{X} + Y_c(t)\mathbf{Y} + R_c\mathbf{u} + u(\theta, t)\mathbf{u} + w(\theta, t)\mathbf{w}) \\ \mathbf{\Omega}_c = \dot{\Phi}_{xc}\mathbf{X}_1 + \dot{\Phi}_{yc}\mathbf{Y}_2 \end{cases} \quad (4) $$

The potential energy and dissipation function related to the elastic ring dofs are given by [14]:

$$ \begin{aligned} \mathcal{U}_c &= \frac{1}{2} \int_0^{2\pi} \frac{E_c I_z}{R_c^4} \left( \frac{\partial^2 u(\theta, t)}{\partial \theta^2} + \frac{\partial w(\theta, t)}{\partial \theta} \right)^2 R_c d\theta \\ \mathcal{F}_c &= \frac{1}{2} \int_0^{2\pi} \frac{\eta_c E_c I_z}{R_c^4} \left( \frac{\partial^2 \dot{u}(\theta, t)}{\partial \theta^2} + \frac{\partial \dot{w}(\theta, t)}{\partial \theta} \right)^2 R_c d\theta \end{aligned} \quad (5) $$

Due to the orthogonal property of the trigonometric functions and the linearisation of the kinetic energy terms, the matrices resulting from the application of Lagrange equations, dissociate the dofs of the casing. The mass matrix related to the rigid-body dofs is denoted $M_c$ and the one related to the $n_d$ deformation mode is denoted $M_{n_d}$. Similarly, the damping and stiffness matrices related to the $n_d$ deformation mode are denoted $D_{n_d}$ and $K_{n_d}$; and the visco-elastic suspension provides potential and dissipation energy respectively leading to the matrices $D_c$ and $K_c$ related to the rigid-body dofs. The casing matrices are given in appendix (Tab 3).

The casing and suspension system of equation must then be transposed into the rotating frame in which the bladed-shaft model is written. We can take advantage of the decoupling of the rigid-body dofs and nodal diameter deformation dofs to transpose the system of equations successively.

Let's begin by the rigid-body dofs related to the mass matrix $M_c$, damping matrix $D_c$ and stiffness matrix $K_c$. The transposition of the sytem of equation into the rotating frame creates a gyroscopic matrix $G_c$, a spin softening matrix $N_c$ and a circulation matrix $C_c$, given by:

$$ G_c = 2P^T M_c \dot{P}; \quad N_c = P^T M_c \ddot{P}; \quad C_c = P^T D_c \dot{P} \quad (6) $$

<figure><img src="image_7.png" /><figcaption>Figure 4. Angular position of a given point on the casing seen from the fixed and rotating frame</figcaption></figure>

The transposition matrix $P$ is detailed in appendix (Tab 3).

To apply the same transformation to the dofs related to the deformation of the $n_d$ nodal diameter, we must establish the relation between the dofs ($A_{n_d}$, $B_{n_d}$) defined in the inertial frame and ($A_{n_d}^*$, $B_{n_d}^*$) defined in the rotating frame. As sketched in figure 4, by expressing the physical displacements $u$ and $w$ of the same cross section depicted by the angular position $\theta$ in the fixed frame and $\psi$ in the rotating frame, we obtain the transposition matrix $P_{n_d}$ given in appendix. As for the rigid body dofs, the transposition in the rotating frame introduces the matrices $G_{n_d}$, $N_{n_d}$ and $C_{n_d}$ by the same relations as equ. (6).

To complete the model formulation, it remains to describe the link considered between the shaft and the casing (see figure 1). The link is modeled by a visco-elastic element with isotropic properties, it connects the rigid body dofs of the casing with the displacement of the facing cross-section of the shaft (positionned by $z_l$). Its potential energy and dissipation function are given in appendix.

The total number of dofs of the model depends on the Ritz function considered to describe the blade, shaft and casing deformations:

$$ \text{ndof} = 2 + 2 \times m_{\text{tot}} + N_b \times n_{\text{tot}} + 4 + 2 \times n_{\text{d,tot}} \quad (7) $$

The model contains 10 blades with 2 dofs each. The bladed wheel brings thus two families of ten modes from 0 to 5 nodal diameters, based on the two first deflexion modes of a clamped-free beam. The casing deformations are introduced by 4 dofs describing the 2 and 3 nodal diameters deformations, and plus 4 rigid-body dofs. Lastly, 4th order Ritz functions are considered to describe the shaft which contribute to 10 dofs. The full model described in this paper contains thus 38 dofs.

## 1.2 Contact formulation

The rotor-to-stator contacts studied in this paper, involve the fan blade tips and the inner surface of the casing coated with abradable material. The fan blades of modern jet engines present a wide chord and a complex 3D shape. In addition, the inner surface of the casing is not cylindrical but slightly conical. These specificities of the local geometry of the contact area sketched in figure 5, do not appear in our
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Geometry of the contact region and 3D view of the clearance values interpolated along the chord line</figcaption></figure>

simplified model, where the casing is a cylindrical ring and the blades, rectangular cross-section beams. This section presents a contact formulation that enables to improve accuracy in the contact detection and application of contact forces without increasing the model size.

### 1.2.1 Minimum distance function

In order to improve accuracy, the contact formulation introduces elements of the local geometry and considers all the kinematics brought by the model. The clearance is evaluated at both ends of the blade tip chord. These points are depicted in figure 5: P(LE) for the Leading Edge and P(TE) for the Trailing Edge 3D coordinates of these two points as well as the inclination of the casing inner surface $\hat{\theta}$ are known parameters of design.

The clearance value at the selected points are obtained by first determining their positions based on the system displacements. The final position $P'$ of the point with initial position P, is given in the frame attached to the casing, by:

$$P' = \mathbf{P}_{\phi_{yc}} \mathbf{P}_{\phi_{xc}} \left( \begin{array}{c} \mathbf{P}_{\phi_{xd}}^{T} \mathbf{P}_{\phi_{yd}}^{T} \mathbf{P}_{\alpha_{j}}^{T} (P \\ + \mathbf{P}_{\beta}^{T} \left\{ \begin{array}{c} x_{b}(L_{b}, t) \\ 0 \\ 0 \end{array} \right\} ) + \left\{ \begin{array}{c} x_{d}(t) - x_{c}(t) \\ y_{d}(t) - y_{c}(t) \\ -z_{c} \end{array} \right\} \end{array} \right) \quad (8)$$

Once the position $P' = \{P'_x, P'_y, P'_z\}^T$ known, the clearance is evaluated by calculating the minimum distance between a point and a cone surface and introducing the radial

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Scheme of the gap line interpolated on the chord and lineic contact pressure applied</figcaption></figure>

deformation of the casing:

$$g_j(P) = \left( -\frac{R_b}{h} P'_z + R_b - \sqrt{P_x'^2 + P_y'^2} + u(\alpha(P), t) \right) \cos \hat{\theta} \quad (9)$$

Note that the minimum distance function considers the radial deformation evaluated at the initial angular position of the blade tip point. The displacement are assumed small enough to allow this approximation.

### 1.2.2 Contact forces

To study the model behavior under contacts, we must write the nonlinear forces acting on the system when the blade tips touch the casing. The penalty method has been chosen to express the contact forces: its implementation is simpler than a Lagrangian multiplier approach and some industrial values of radial stiffness parameters corresponding the abradable coatings were available to our analysis. The normal contact force $F_N$ is thus proportional to penetration: the radial stiffness parameter is denoted $k_r$. The tangential forces due to rubbing are obtained using Coulomb's Law: $-\operatorname{sgn}(V_{rel})\mu F_N$ where $\mu$ is the friction coefficient and $V_{rel}$ the relative speed between contact surfaces. $\operatorname{sgn}(V_{rel})$ is assumed to be constant and opposite to the sliding speed direction due to shaft spin.

The clearances evaluated at LE and TE will thus allow to interpolate the clearance along the blade tip chord line. A corresponding lineic pressure is then applied on the portion of chord touching the casing. This strategy is sketched in figure 6. The normal and tangential contact pressures at the position c of the chord line are given by:

$$p_{jN}(c) = \frac{k_r}{l_c}(\zeta c + \eta); \quad p_{jT}(c) = -\operatorname{sgn}(V)\mu p_{jN}(c) \\ \text{with } \zeta = \frac{g_j(TE) - g_j(LE)}{l_c}; \quad \eta = \frac{g_j(TE) + g_j(LE)}{2} \quad (10)$$

The linear interpolation chosen here assume the blade tip chord straight. This restrictive assumption is compensated by the additional assumption of the casing portion facing the blade tip being also straight.

The contact reactions at the $j^{th}$ blade tip are obtained by integrating the pressures defined in eq. (10):

$$F_{jN} = \int_{c_1}^{c_2} p_{jN}(c)dc \quad \text{and} \quad F_{jT} = \int_{c_1}^{c_2} p_{jT}(c)dc \quad (11)$$
An uneven distribution along the chord can also creates reaction torques given by:

$$M_{jN} = \int_{c_1}^{c_2} c \, p_{jN}(c) \, dc \quad \text{and} \quad M_{jT} = \int_{c_1}^{c_2} c \, p_{jT}(c) \, dc \quad (12)$$

$c_1$ and $c_2$ indicates the limit coordinates of the portion of chord in contact.

The forces acting on the blade tip and the casing can then be projected onto the general coordinates of the system.

## 1.3 Analysis

Resolution is performed in the rotating frame. The system of equations of the model subjected to unbalance and blade to casing contact can thus be written as follow:

$$M\ddot{q}(t) + \ddot{D}(\Omega)\dot{q}(t) + \tilde{K}(\Omega)q(t) = F_b(\Omega) + F_{nl}(q(t)) \quad (13)$$

An initial gap of 0.5 mm is considered over the entire chord length and a 0.225 m.kg unbalance is applied at the disk representative of the LP turbine. This unbalance amplitude and location has been chosen to excite enough the two global mode having critical speeds in the studied speed range. The unbalance amplitude would be too high to be considered as a residual unbalance at the LP turbine, however we focus here on the displacements seen at the fan. The chosen unbalance permits, with only one set of parameters, to excite two different modes at sufficient amplitude to reach blade-to-casing contacts.

### 1.3.1 Permanent contacts: static resolution

The system is build in the rotating frame, which means the synchronous unbalance forces constitute a static loading. Considering only permanent contact allows thus to assume a static solution, the non-linear equation can be written:

$$\tilde{K}(\Omega)q_0 = F_b(\Omega) + F_{nl}(q_0) \quad (14)$$

This equation is solved by using a Newton-Raphson algorithm. Once the static solution $q_0$ found, the analysis of the local stability of the solution comes down to the following eigenvalues problem:

$$\det(\lambda^2 M + \lambda \tilde{D}(\Omega) + K(\Omega) - J_q F_{nl}(q_0)) = 0 \quad (15)$$

where $J_q F_{nl}(q_0)$ is the jacobian matrix of the contact forces vector at $q_0$. The eigenvalues real parts give stability information of the static solution: the solution being unstable if it has a positive real part [15].

### 1.3.2 Intermittent contacts: transient analysis

The transient results presented in this paper are obtained with an explicit central differences scheme as in [5, 16]. The chosen time step has been $4 \times 10^{-5}$ sec.

The higher eigen-frequency within the studied speed range is 7500 Hz: for a linear system, the time step ensuring the numerical scheme stability would be $2/\omega \approx 4.2 \times 10^{-5}$. The studied problem being non-linear, this stability criteria do not necessarily apply: the time step size convergence has thus

been verified by reducing the step up to $1 \times 10^{-5}$ sec for some transient simulations (divergent or not) with similar loading and model configuration.

The simulation are performed considering a constant spin speed and applying unbalance over a 1 sec ramp. Their maximal duration is fixed to 5 sec and the simulation is stopped if the normal contact load at one blade tip reaches 8000 N. In order to compare the static and transient results, FFT analyses are performed on the obtained transient signals. The frame chosen for the FFT depends of the signal result: if no contact occurs during the simulation, the frame begins at 1.5 sec (i.e. 0.5 sec after the full unbalance is applied); if contact occurs, the frame begins 0.5 sec after the first contact.

## 2. RESULTS AND DISCUSSION

### 2.1 Model natural properties

The Campbell diagram on figure 7 presents the evolution of the eigen-frequencies with the shaft spin speed. Four families of modes are highlighted in the diagram:

* The green lines denoted C depict the eigen-frequencies of the casing deformation modes: the two and three nodal diameters. The variation of their eigen-frequencies is only due to the description in the rotating frame, the corresponding lines of the co and contra rotative modes would be perfectly horizontal if described in the inertial frame.

* The red lines denoted B depict the bladed-disk modes of the first flexion family. The lines of the 0 and 2 to 5 nodal diameters modes form a group with very similar frequencies which increase with spin speed due to stress stiffening and spin softening. The B($1n_d$) lines are not among the group of the other nodal diameters because of the coupling of this particular modes with the shaft dynamics. This inertial coupling has been described in [17, 12]; the explanation lies in the compatibility of the $1n_d$ bladed disk mode resulting effort on the disk, with the kinematic of the shaft cross-section. The blades-shaft coupling results in various locus-veering on the Campbell diagram: the eigenvalues are nearly about to cross, but instead they veer away from one another and exchange modeshapes. It is also visible on the modeshapes sketched besides the graph: where we can see the contribution of the blades deformation to the global modeshapes.

* The dashed lines denoted W depict the eigen-frequencies obtained with a degraded configuration of the model: the blades are assumed rigid and considered only through their inertia properties. The W modes thus form the Whole Engine Model modes usually observed on rotordynamics centerline models. These WEM modes encompass the shaft deflexion and casing rigid-body displacements which are connected by the elastic link.

* The blue lines denoted W' depict the eigen-frequencies of these WEM modes with the contribution of the blades flexibility. The W'(1,2) and W'(7,8) modeshapes at 0 rpm are sketched besides the Campbell. Their forward modes W'(1) and W'(7) have their critical speed among the studied range; and they are respectively called Fan and LP Turbine modes because they concentrate most of their deformation energy in
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Campbell diagram plotted in the rotating frame and modeshapes of the Fan mode:W'(1,2) and the LP turbine mode:W'(7,8) sketched at 0 rpm</figcaption></figure>

these areas. The influence of the blades-shaft coupling on the WEM modes can be observed on the Campbell diagram, where the eigen-frequencies evolution of the W' differs from the corresponding W modes and various locus-veering appear. The modeshapes sketched illustrates this coupling and illustrate the effect on the W'(1) mode stronger than on the W'(7) (visible also on the Campbell). The W'(1) concentrates indeed most of its shaft deformation at the fan and thus tends more to couples with the blade deflexion (B(1n<sub>d</sub>) modes).

The phenomenological model, containing only 38 dofs, presents an elaborated behaviour due to the various couplings considered. Its 1 nodal diameter modeshape contains deformation energy in the entire model: involving the blades, the shaft and the casing displacement. In the following sections, the role played by these couplings and the additional links brought by the contact will be analysed.

## 2.2 Permanent contact results

The eigenvalue calculated at each static equilibria found from 3000 to 8000 rpm, are plotted in figure 8. Each rpm step correspond to a new resolution of the static equation, the contact configuration can thus change drastically at each new rpm step, leading to important variations of the corresponding eigenvalues.

The speed ranges related to three results status are highlighted on the graph :

* When the response to unbalance have not led to sufficient displacement to consume the initial clearance, there is no contact. We encounter then the modal properties of the linear model. On these speed ranges, the background diagram is white.

* The shaded areas depicted *contact and stable*, illustrate the speed ranges for which contacts occur, but the system damping is sufficient to compensate the effect of

the assymetries brought by the contact and the real parts of all eigenvalues remain negative.

* The colored areas depicted *contact and unstable*, indicate the speeds for which the static solutions are unstable; at least one of the eigenvalue calculated on the equilibrium have a positive real part.

To better perceive the deviations from the model behavior without contact, these results are superimposed in dotted lines in the diagrams. Note also that the eigenvalues of unstable modes are marked with a larger dot whose color reflects the value of the corresponding real part: a high real part corresponding to a darker point. It allows to simply illustrate the degree of instability. A positive real part but of low amplitude may indeed be easily compensated by a slight increase of damping. Conversely, a high positive real part will be more difficult to compensate.

The general observation of the diagram indicate that the distribution of the speed ranges (*Contact-less*, *Contact and stable* and *Contact and unstable*) is consistent with the system response to unbalance. The two critical speeds contained in the studied speed range (4450 rpm for the fan mode W'(1) and 7150 rpm for the LP turbine mode W'(7) ) are indeed associated to the highest response levels and thus to the greatest clearance reduction. The residual penetrations reach -0.8 mm at the first critical speed and -0.15 mm at the second (the clearance reductions calculated without considering contact forces lead respectively to the penetration depth of -1.3 mm and -0.3 mm). Near these speeds, the contact configurations are the most critical, with many blades in contact over large portion of their chord; unstable ranges are therefore naturally present around these critical speeds. The unstable regions are then surrounded by speed ranges, for which the clearance consumption is still sufficient to lead to contacts, but their effect are outweighed by the viscous damping of the system.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Eigenvalues calculated around the static equilibria considering permanent contact due to unbalance: stability analysis</figcaption></figure>

The instability obtained on the range located around the first critical speed is related to the eigenvalues real part of two backward modes: the B'(1$n_d$) mode and also the W'(4) mode. Their modeshapes are influenced by the introduction of contact leading to fully coupled modes including the participation of bladed disk modes different from the 1 nodal diameter and also of the nodal deformation of the casing. For the speed range around the second critical speed, these same modes are destabilised but the W'(4) presents here higher real parts.

Even if they stay stable, the nodal diameters modes of the bladed disk B and the casing C have their natural frequencies changing the most due to contact. These modes concentrate their deformation energy in the blades deflection or the casing deformation: at the direct contact interface, they are thus strongly affected when contact occurs. For the bladed disk modes, the changes of their natural frequencies alternate between stiffening and softening with the evolution of contact configuration. Whereas the casing modes are mode evenly softened on the speed ranges with contacts. This can be explained by their respective modeshapes. In this study, contact are introduced by the response to unbalance, which means the touches will establish on one side of the bladed wheel; on one or two blades and up to about half of the bladed wheel as the speed steps calculated approach the critical speeds. From the casing point of view, a new or a lost contact leads to a change in the distribution of the touches on its inner surface. From the bladed wheel point of view, the change is more radical for

its nodal diameter modes which see a new blade connected to the casing or conversely be released. Likewise, if to the 2 and $3n_d$ casing modes, contacts provide supports leading to a lower eigenfrequency; for the bladed disk modes, the stiffening or softening effect will depend on the disposition of the touches regarding the concerned modeshape.

This static results allow us to predict the system stability when submitted to contact and indicates the important role of couplings to the dynamic stability and global behaviour under contacts. The participation of the couplings naturally present in the model is highlighted as well as the effects of the couplings brought by the blade-to-casing contacts.

The assumption of a static solution is however a strong restriction of the system behaviour. Due to its non-linearity, the nature or number of solutions can not be known a priori. The static equilibria encountered are therefore not necessarily the unique solutions; and by considering the dynamic effects, one unstable equilibria may correspond to anharmonic solution (stable or not) or to a chaotic solution. Another strong disadvantage of the assumption of a static solution is the incompatibility with intermittent contacts. The transient analysis presented in the next section will thus provide a validation for the static results.

## 2.3 Transient results

Transient simulation are carried out at regular intervals of 200 rpm, from 3200 to 7600 rpm. The three types of results mentioned above are encountered: solutions for which no contact is detected; the simulations that despite contact, result in a stable equilibrium and finally the unstable solutions for which the displacement field diverges.

Figure 9 presents the frequency spectra obtained for each transient simulation performed from 3200 to 7600 rpm. Some informations from the static analysis have been plotted on the graph. The speed ranges corresponding to the different types of results (*Contactless*, *Contact and stable* and *Contact and unstable*) are indicated with coloured areas on the graph basis. Note that the type of result obtained through static and transient analysis is identical. The equilibria corresponding levels are also plotted in the zero frequency plane of the graph. Lastly, the frequencies of the most unstable modes (related to the eigenvalue having the highest real parts) are indicated by dots whose color reflects the value of the corresponding real part.

Figure 9 allows thus to compare the results obtained through each strategy, in terms of spectral content and zero frequency contribution. Note that for an unstable system (colored speed ranges), the peak of the signal spectrum obtained by transient analysis, coincides with the frequency of the unstable mode identified by the static approach. Conversely, when the system is stable (with or without contact), the transient signal is constant and the related frequency spectrum exhibits only a component at 0 Hz. These 0 Hz components of the transient signal spectrum follow the amplitudes of static equilibria plotted in the same plane.

The good agreement encountered between permanent con-
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Frequency spectrum of the transient simulations carried out from 3200 to 7600 rpm: comparison with the static results (see figure 8)</figcaption></figure>

tact and transient analyses, verifies the stability conclusions obtain with the first method. The static resolution has indeed a computational cost substantially lower : less than 1 min calculation for the 500 rpm steps presented in figure8, while each transient simulation necessitates 5 min calculation. The permanent contact strategy allows to perform extensive influence analyses by varying: unbalance amplitude, contact parameters or clearance levels ... and then give the suitable design recommendations.

From another point of view, the agreement between the unstable frequencies encountered with the permanent contact analyses and the spectral content of the transient signal, verifies also the consistency of transient results (chosen time-step).

The transient simulations contains also more information about system behavior during contact and stability loss. The following results study the casing and disk trajectories during the unstable scenario obtained at 7000 rpm.

### 7000rpm scenario

The results plotted on figure 9, have allowed to compare the two approaches. In this section, the density of informations available in the transient simulation is exploited to monitor an unstable scenario, from the initiation of contacts to the divergence of the system. The chosen scenario is the one obtained at 7000 rpm.

Figure 10 shows the evolution of normal contact reactions and the trajectories of the centers of gravity of the fan disk and the casing described in the rotating frame.

The trajectories from 0 to 1.5 sec [T1] illustrates the system static behaviour during application of the unbalance and the permanent contact state. This static contribution corresponds

in the inertial frame, to a synchronous forward circular precession widening as the unbalance progresses. The blades gradually come into contact, but the trajectories stay circular and of constant amplitude, corresponding to the static position observed from the rotating frame.

From 1.5 to 2.5 sec [T2], a backward whirl in the rotating frame widen for the disk and casing. The disk trajectory is nearly circular while the casing trajectory is less regular: its flattens and the precession center moves. The trajectories observations around 3.25 sec [T3], indicate that they both diverge in a forward orientation: the whirling direction changed during the simulation. The switches are not simultaneous: for the casing it occurs at circ. 2.9 sec; the trajectory becomes more irregular, flattened and form an 8 to eventually widen again in the forward direction. Disk and casing then maintain opposite whirling until circ. 3.2 sec, when the disk whirling orientation switches to forward.

## CONCLUSION

This study has presented a phenomenological model developed to investigate phenomena related to blade-to-casing contacts in turbofan aircraft engines.

The model contains only 38 dofs and comprises a flexible bladed shaft and a flexible suspended casing. A shaft-casing coupling representative of the fan frame is introduced and the blades-shaft inertial coupling considered. The model can thus yields complex mode shapes encompassing all the components. A 3D contact formulation has been implemented to introduce the blade-to-casing contacts. Its permits to recover accuracy in the touches detection and application of contact reactions without increasing the computational cost.
<figure><img src="image_13.png" /><figcaption>T1: Trajectories from 0 to 1.5 sec</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>T2: Trajectories from 1.5 to 2.5 sec</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>T3: Trajectories from 2.5 to 3.25 sec</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 10.** Transient analysis at 7000 rpm, trajectories of the center of gravity of the fan disk and casing described in the rotating frame</figcaption></figure>
Static results obtained by assuming permanent contact has been presented. They highlight the relevance of the global modes and the additional coupling brought by the contacts. The stability status obtained through static analysis is then confirmed by the transient simulations carried out under similar conditions. For unstable equilibria, the spectral content of the diverging displacement field exhibits the unstable eigen-frequencies identified through the static approach.

Transient analysis also permits to follow the disk and casing trajectories during an unstable scenario. The trajectories take complex shape, unusual orientation and even switch direction during the simulation. This rich evolution of trajectories highlight the combined influence of 3D contact formulation, the blades flexibility and shaft-casing link introduced to the model which complicates the kinematics at play in the unstable phenomenon.

Further investigations could focus on the evolution of the nodal diameters deformation observed on the bladed wheel and the casing (some results can be found in [18]). In order to fully study the unstable system, the model could also be completed to take into account the removal of the abradable material.

## ACKNOWLEDGMENTS

The authors would like to thank the Snecma firm and SAFRAN Group for its technical and financial support.

## REFERENCES

[1] A. Muszynska. Rotor-to-stationary element rub-related vibration phenomena in rotating machinery: Literature survey. *The Shock and Vibration digest*, 21:3–11, 1989.

[2] S. Ahmad. Rotor casing contact phenomenon in rotor dynamics — literature survey. *Journal of Vibration and Control*, 16:1369–1377, 2010.

[3] G. Jacquet-Richardet, M. Torkhani, P. Cartraud, F. Thouverez, T. Nouri, M. Baranger, M. Herran, C. Gibert, S. Baguet, P. Almeida, and L. Peletan. Rotor to stator contacts in turbomachines. review and application. *Mechanical Systems and Signal Processing*, 40:401–420, 2013.

[4] P. Schmiechen. *Travelling Wave Speed Coincidence*. PhD thesis, Imperial College London, 1997.

[5] M. Legrand, C. Pierre, P. Cartraud, and J.-P. Lombard. Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction. *Journal of Sound and Vibration*, 319(1-2):366–391, 2009.

[6] M. Meingast, A. Batailly, M. Legrand, and J.-P. Osty. Investigation of rotor-casing interactions in the centrifugal compressor of a helicopter engine. In *ASME Turbo Expo, June 3-7, San Antonio (USA)*. American Society of Mechanical Engineers, 2013.

[7] P. Almeida, C. Gibert, F. Thouverez, X. Leblanc, and J.-P. Osty. Experimental analysis of dynamic interaction between a centrifugal compressor and its casing. In *ASME Turbo Expo, Düsseldorf, June 16-20, 2014*.

[8] S. K. Sinha. Dynamic characteristics of a flexible bladed-rotor with coulomb damping due to tip-rub. *Journal of Sound and Vibration*, 273:875–919, 2004.

[9] N. Lesaffre, J.-J. Sinou, and F. Thouverez. Contact analysis of a flexible bladed-rotor. *European Journal of Mechanics - A/Solids*, 26(3):541–557, 2007.

[10] M-O. Parent, F. Thouverez, and F. Chevillot. Whole engine interaction in bladed rotor-to-stator contact. In *ASME Turbo Expo, Düsseldorf, Germany, June 15-20, 2014*.

[11] M-O. Parent, F. Thouverez, and F. Chevillot. 3d interaction in bladed rotor-to-stator contact. In *EURODYN, June 30- July 2, Porto, Portugal*, 2014.

[12] N. Lesaffre, J.-J. Sinou, and F. Thouverez. Model and stability analysis of a flexible bladed rotor. *International Journal of Rotating Machinery*, 2006:1–16, 2006.

[13] K. F. Graff. *Wave motion in elastic solids*. Dover Publications, 1975.

[14] Z. Celep. In-plane vibrations of circular rings on a tensionless foudation. *Journal of Sound and Vibration*, 143:461–471, 1990.

[15] T. Gmur. *Dynamique des Structures, Analyse modale numerique*. Presses Polytechniques et Universitaires Roumandes, 2007.

[16] M. Legrand, A. Batailly, B. Magnain, P. Cartraud, and C. Pierre. Full three-dimensional investigation of structural contact interactions in turbomachines. *Journal of Sound and Vibration*, 331(11):2578–2601, 2012.

[17] M. Gruin, F. Thouverez, and L. Blanc. Nonlinear dynamics of a bladed dual shaft. *European Journal of Computational Mechanics*, 20:207–225, 2011.

[18] M-O. Parent. *Prédiction de la stabilité en contact rotor-stator dans les turboréacteurs d'avion*. PhD thesis, Ecole Centrale de Lyon, 2015.

## APPENDIX

### Energies expressions for the bladed-shaft

The shaft deformations are given by:

$$ 
\begin{aligned}
x_s(z,t) &= X_0(t) + \sum_{m=1}^{m_{tot}} X_m(t)W_m(z) \\
y_s(z,t) &= Y_0(t) + \sum_{m=1}^{m_{tot}} Y_m(t)W_m(z) \\
\text{with } W_m(z) &= \left(\frac{z}{L_s}\right)^m
\end{aligned} 
\quad (16) 
$$

The shape function $W_m(z)$ has been used in [17]; it allows motion at both ends of the shaft.
The blade deformation is given by:

$$
\begin{align}
x_{b_j}(s,t) &= \sum_{n=1}^{n_{tot}} X_{nj}(t)Y_n(s) \nonumber \\
Y_n(s) &= \sin(a_n s) + b_n s + c_n s^2 \nonumber \\
a_n &= \frac{(2n-1)\pi}{2L_b}; \quad b_n = -a_n; \quad c_n = \frac{b_n^2}{2} \sin(b_n L_b) \tag{17}
\end{align}
$$

The shape function of the blades, $Y_n(s)$, respects the geometric boundary conditions: $Y_n(0) = 0$ and $Y'_n(0) = 0$ and has been used in [17] and [8].

The shaft, disks and blades general expression of the ki-
netic energy is given by:

$$
\begin{equation}
\begin{aligned}
\mathcal{T}_{bs} ={}& \frac{1}{2} \rho_s S_s \int_0^{L_s} \mathbf{V}_s^T(z,t) \mathbf{V}_s(z,t) dz \\
& + \frac{1}{2} \rho_s \int_0^{L_s} \mathbf{\Omega}_s^T(z,t) \mathbf{I}_s \mathbf{\Omega}_s(z,t) dz \\
& + \sum_{d=1,2} \left( \frac{m_d}{2} \mathbf{V}_s^T(z_d,t) \mathbf{V}_s(z_d,t) \right. \\
& \qquad \left. + \frac{1}{2} \mathbf{\Omega}_s^T(z_d,t) \mathbf{I}_d \mathbf{\Omega}_s(z_d,t) \right) \\
& + \sum_{j=1}^{10} \frac{1}{2} \rho_b S_b \int_0^{L_b} \mathbf{V}_b^T(s,t) \mathbf{V}_b(s,t) ds \\
& + \frac{1}{2} \rho_b \int_0^{L_b} \mathbf{\Omega}_b^T(s,t) \mathbf{I}_b \mathbf{\Omega}_b(s,t) ds
\end{aligned}
\tag{18}
\end{equation}
$$

The formulation of the translation and rotation instantaneous
speed : $\mathbf{V}_s$ and $\mathbf{\Omega}_s$ for the shaft cross-section and $\mathbf{V}_b$ and $\mathbf{\Omega}_b$
for the blade cross-section are given in [12]. The correct for-
mulation of the blade instantaneous speeds ensure the inertial
link between shaft’s and blade’s dofs.

The potential energy brought by the isotropic bearing and
the elastic deformations of the bladed-shaft is given by:

$$
\begin{equation}
\begin{aligned}
\mathcal{U}_{bs} ={}& \frac{k_p}{2} (x_s^2(z_p, t) + y_s^2(z_p, t)) \\
& + \frac{E_s I_s}{2} \int_0^{L_s} (x_s''^2(z, t) + y_s''^2(z, t)) dz \\
& + \sum_{j=1}^{10} \frac{E_b I_b / z_b}{2} \int_0^{L_b} x_b''^2(s, t) ds \\
& + \frac{1}{2} k_b (x_{bj+1}(L_b, t) - x_{bj}(L_b, t))^2
\end{aligned}
\tag{19}
\end{equation}
$$

To organise the blade modes according to nodal diameters, an
elastic link ($k_b$) is introduced between each consecutive blade
to represent the effect of disk flexibility [12, 17].

Centrifugal stiffening effect is considered, the correspond-
ing potential on the j^th blade is defined by:

$$
\mathcal{U}_{b,stif} = \frac{\Omega^2 \rho_b S_b}{2} \int_0^{L_b} \frac{(r_d + L_b)^2 - (r_d + s)^2}{2} (x_b''(s,t))^2 ds
$$

(20)

Lastly, the dissipation function brought by the bearings,
the shaft and the blades is given by:

$$
\begin{equation}
\begin{split}
\mathcal{F}_{bs} ={}& \frac{c_p}{2} \left( (\dot{x}_s(z_p, t) - \Omega y_s(z_p, t))^2 \right. \\
& \qquad \left. + (\dot{y}_s(z_p, t) + \Omega x_s(z_p, t))^2 \right) \\
& + \frac{\eta_s E_s I_s}{2} \int_0^{L_s} \left( \dot{x}_s''^2(z, t) + \dot{y}_s''^2(z, t) \right) dz \\
& + \sum_{j=1}^{10} \frac{\eta_b E_b I_b / z_b}{2} \int_0^{L_b} \dot{x}_b''^2(s, t) ds
\end{split}
\tag{21}
\end{equation}
$$

Energies expressions for the shaft-casing link

The deformation energy of the isotropic elastic link consid-
ered between the casing and the shaft is given by:

$$
\begin{equation}
\begin{aligned}
\mathcal{U}_l ={}& \frac{1}{2} k_l \left( (x_c(z_l, t) - x_s(z_l, t))^2 \right. \\
& \qquad \left. + (y_c(z_l, t) - y_s(z_l, t))^2 \right) \\
& + \frac{1}{2} k_{l,\phi} \left( (\phi_{x_c}(t) - \phi_{x_s}(z_l, t))^2 \right. \\
& \qquad \left. + (\phi_{y_c}(t) - \phi_{y_s}(z_l, t))^2 \right)
\end{aligned}
\tag{22}
\end{equation}
$$

The dissipation function related to its viscous damping is
given by:

$$
\begin{equation}
\begin{split}
\mathcal{F}_l ={}& \frac{1}{2} c_l \left( (\dot{x}_c(z_l, t) - \Omega y_c(z_l, t) - \dot{x}_s(z_l, t) \right. \\
& \qquad \left. + \Omega y_s(z_l, t) \right)^2 + (\dot{y}_c(z_l, t) + \Omega x_c(z_l, t) \right. \\
& \qquad \left. - \dot{y}_s(z_l, t) - \Omega x_s(z_l, t) \right)^2 \right) \\
& + \frac{1}{2} c_{l,\phi} \left( (\dot{\phi}_{x_c}(t) - \Omega \phi_{y_c}(t) - \dot{\phi}_{x_s}(z_l, t) \right. \\
& \qquad \left. + \Omega \phi_{y_s}(z_l, t) \right)^2 + (\dot{\phi}_{y_c}(t) + \Omega \dot{\phi}_{y_c}(t) \right. \\
& \qquad \left. - \dot{\phi}_{y_s}(z_l, t) - \Omega \dot{\phi}_{x_s}(z_l, t) \right)^2 \right)
\end{split}
\tag{23}
\end{equation}
$$

The rigid body translation of the casing at the link position
are given by:

$$
\begin{align}
x_c(z_l, t) &= x_c(t) + z_l \tan \phi_{y_c}(t); \nonumber \\
y_c(z_l, t) &= y_c(t) - z_l \tan \phi_{x_c}(t) \tag{24}
\end{align}
$$

Note that the deformation dofs of the casing are not consid-
ered in the energies expressions (22) et (23), they are indeed
filtered by the shaft kinematic.

**Terminology and data sets**
Blades

<table><caption>**Table 1.** Flexible blades properties</caption><tr><td>η<sub>b</sub></td><td>Viscous damping factor</td><td>4 × 10<sup>-5</sup></td></tr><tr><td>ρ<sub>b</sub></td><td>Density</td><td>4450 kg/m<sup>3</sup></td></tr><tr><td>E<sub>b</sub></td><td>Young modulus</td><td>114 × 10<sup>9</sup> Pa</td></tr><tr><td>β</td><td>Orientation of rectangular cross-section</td><td>-10°</td></tr><tr><td>L<sub>b</sub></td><td>Length</td><td>0.52 m</td></tr><tr><td>r<sub>df</sub></td><td>Fan disk radius</td><td>0.23 m</td></tr><tr><td>b</td><td>Small edge of the cross-section</td><td>0.0165 m</td></tr><tr><td>h</td><td>Large edge of the cross-section</td><td>0.1165 m</td></tr><tr><td>I<sub>b</sub></td><td>Matrix of area moment of inertia</td><td><sup>1</sup>&frasl;<sub>12</sub> [ <sup>bh<sup>3</sup></sup>&frasl;<sub>hb<sup>3</sup> + bh<sup>3</sup></sub> <sup>hb<sup>3</sup></sup>&frasl;<sub>hb<sup>3</sup></sub> ]</td></tr><tr><td>k<sub>b</sub></td><td>Inter-blade stiffness</td><td>100 N/m</td></tr><tr><td>n<sub>tot</sub></td><td>Order of the Ritz function</td><td>2</td></tr><tr><td>N<sub>b</sub></td><td>Number of blades</td><td>10</td></tr></table>

Shaft

<table><caption>**Table 2.** Shaft properties</caption><tr><td>η<sub>s</sub></td><td>Viscous damping factor</td><td>0</td></tr><tr><td>ρ<sub>s</sub></td><td>Density</td><td>7800 kg/m<sup>3</sup></td></tr><tr><td>E<sub>s</sub></td><td>Young modulus</td><td>119 × 10<sup>9</sup> Pa</td></tr><tr><td>L<sub>s</sub></td><td>Length</td><td>2 m</td></tr><tr><td>m<sub>tot</sub></td><td>Order of the Ritz function</td><td>4</td></tr><tr><td>R<sub>e</sub></td><td>Cross-section external radius</td><td>0.065 m</td></tr><tr><td>R<sub>i</sub></td><td>Cross-section internal radius</td><td>0.055 m</td></tr><tr><td></td><td>Spin orientation</td><td>-z</td></tr><tr><td>I<sub>s</sub></td><td>Matrix of area moment of inertia</td><td>[ I<sub>d</sub> I<sub>i</sub> ]</td></tr><tr><td>I<sub>d</sub></td><td>Polar area moment of inertia</td><td>π(R<sub>e</sub><sup>4</sup> - R<sub>i</sub><sup>4</sup>)</td></tr><tr><td>I<sub>p</sub></td><td>Diametral area moment of inertia</td><td>π(R<sub>e</sub><sup>4</sup> - R<sub>i</sub><sup>4</sup>)</td></tr><tr><td>z<sub>p</sub></td><td>Bearing coordinate</td><td>2 m</td></tr><tr><td>k<sub>p</sub></td><td>Bearing stiffness</td><td>1 × 10<sup>8</sup> N/m</td></tr><tr><td>c<sub>p</sub></td><td>Bearing viscous damping property</td><td>8 × 10<sup>3</sup> N.s/m</td></tr><tr><td>z<sub>l</sub></td><td>Shaft-casing link coordinate</td><td>0.5 m</td></tr><tr><td>k<sub>l</sub></td><td>Stiffness</td><td>1.75 × 10<sup>8</sup> N/m</td></tr><tr><td>k<sub>l,φ</sub></td><td>Angular stiffness</td><td>2.5 × 10<sup>8</sup> N/rad</td></tr><tr><td>c<sub>l</sub></td><td>Viscous damping</td><td>14 × 10<sup>3</sup> N.s/m</td></tr><tr><td>c<sub>l,φ</sub></td><td>Angular viscous damping</td><td>5 × 10<sup>3</sup> N.s/m</td></tr></table>

Casing

<table><caption>**Table 3.** Casing properties</caption><tr><td>n<sub>d</sub></td><td>Nodal diameters considered</td><td>{2, 3}</td></tr><tr><td>n<sub>d,tot</sub></td><td>Total number of nodal diameters</td><td>2</td></tr><tr><td>η<sub>c</sub></td><td>Viscous damping factor</td><td>4 × 10<sup>-4</sup></td></tr><tr><td>ρ<sub>c</sub></td><td>Density</td><td>2710 kg/m<sup>3</sup></td></tr><tr><td>E<sub>c</sub></td><td>Young modulus</td><td>69 × 10<sup>9</sup> Pa</td></tr><tr><td>R<sub>c</sub></td><td>Radius</td><td>0.76 m</td></tr><tr><td>I<sub>z</sub></td><td>Area moment of inertia for flexure</td><td>5.9 × 10<sup>-6</sup> m<sup>4</sup></td></tr><tr><td>z<sub>c</sub></td><td>Coordinate</td><td>0.2 m</td></tr><tr><td>m<sub>c</sub></td><td>Mass</td><td>165 kg</td></tr><tr><td>I<sub>d,c</sub></td><td>Diametral moment of inertia</td><td>54.4 kg.m<sup>2</sup></td></tr><tr><td>k<sub>c</sub></td><td>Suspension stiffness</td><td>1.3 × 10<sup>8</sup> N/m</td></tr><tr><td>k<sub>c,φ</sub></td><td>Suspension angular stiffness</td><td>3 × 10<sup>8</sup> N/rad</td></tr><tr><td>c<sub>c</sub></td><td>Viscous damping suspension</td><td>10 × 10<sup>3</sup> N.s/m</td></tr><tr><td>c<sub>c,φ</sub></td><td>Viscous damping angular property</td><td>6 × 10<sup>3</sup> N.rad/m</td></tr><tr><td>M<sub>c</sub></td><td>Mass matrix related to rigid-body dofs</td><td>[<sup>m<sub>c</sub></sup>&frasl;<sub>I<sub>d,c</sub></sub> <sup>m<sub>c</sub></sup>&frasl;<sub>I<sub>d,c</sub></sub>]</td></tr><tr><td>M<sub>n<sub>d</sub></sub></td><td>Mass matrix related to n<sub>d</sub> nodal deformation mode</td><td><sup>m<sub>c</sub></sup>&frasl;<sub>2</sub>(1 + n<sub>d</sub><sup>2</sup>) [ <sup>1</sup>&frasl;<sub>0</sub> <sup>1</sup>&frasl;<sub>1</sub> ]</td></tr><tr><td>D<sub>c</sub></td><td>Viscous damping matrix related to rigid-body dofs (suspension)</td><td>[<sup>c<sub>c</sub></sup>&frasl;<sub>c<sub>c,φ</sub></sub> <sup>c<sub>c</sub></sup>&frasl;<sub>c<sub>c,φ</sub></sub>]</td></tr><tr><td>D<sub>n<sub>d</sub></sub></td><td>Viscous damping matrix related to n<sub>d</sub> nodal deformation mode</td><td><sup>η<sub>c</sub>E<sub>c</sub>I<sub>z</sub>π</sup>&frasl;<sub>R<sub>c</sub><sup>2</sup></sub>n<sub>d</sub><sup>2</sup>(n<sub>d</sub><sup>2</sup> - 1)<sup>2</sup> [ <sup>1</sup>&frasl;<sub>0</sub> <sup>1</sup>&frasl;<sub>1</sub> ]</td></tr><tr><td>K<sub>c</sub></td><td>Stiffness matrix related to rigid-body dofs (suspension)</td><td>[<sup>k<sub>c</sub></sup>&frasl;<sub>k<sub>c,φ</sub></sub> <sup>k<sub>c,φ</sub></sup>&frasl;<sub>k<sub>c,φ</sub></sub>]</td></tr><tr><td>K<sub>n<sub>d</sub></sub></td><td>Stiffness matrix related to n<sub>d</sub> nodal deformation mode</td><td><sup>E<sub>c</sub>I<sub>z</sub>π</sup>&frasl;<sub>R<sub>c</sub><sup>2</sup></sub>n<sub>d</sub><sup>2</sup>(n<sub>d</sub><sup>2</sup> - 1)<sup>2</sup> [ <sup>1</sup>&frasl;<sub>0</sub> <sup>1</sup>&frasl;<sub>1</sub> ]</td></tr><tr><td>P</td><td>Transposition matrix for the rigid-body dofs</td><td>[<sup>cos Ωt - sin Ωt</sup>&frasl;<sub>sin Ωt</sub> <sup>cos Ωt</sup>&frasl;<sub>sin Ωt</sub>]</td></tr><tr><td>P<sub>n<sub>d</sub></sub></td><td>Transposition matrix for the n<sub>d</sub> nodal deformation dofs</td><td>[<sup>cos(n<sub>d</sub>Ωt) - sin(n<sub>d</sub>Ωt)</sup>&frasl;<sub>sin(n<sub>d</sub>Ωt) cos(n<sub>d</sub>Ωt)</sub>]</td></tr></table>
Fan disk

<table><tr><td>z<sub>df</sub></td><td>Coordinate</td><td>0 m</td></tr><tr><td>M<sub>df</sub></td><td>Mass</td><td>100 kg</td></tr><tr><td>I<sub>diam,df</sub></td><td>Diametral moment of inertia</td><td>1.86 kg.m<sup>2</sup></td></tr><tr><td>I<sub>pol,df</sub></td><td>Polar moment of inertia</td><td>3.14 kg.m<sup>2</sup></td></tr></table>

LP turbine representative disk

<table><caption>**Table 4.** Disks properties</caption><tr><td>z<sub>d<sub>t</sub></sub></td><td>Coordinate</td><td>1.95 m</td></tr><tr><td>M<sub>d<sub>t</sub></sub></td><td>Mass</td><td>170 kg</td></tr><tr><td>I<sub>diam,d<sub>t</sub></sub></td><td>Diametral moment of inertia</td><td>13.40 kg.m<sup>2</sup></td></tr><tr><td>I<sub>pol,d<sub>t</sub></sub></td><td>Polar moment of inertia</td><td>31.66 kg.m<sup>2</sup></td></tr><tr><td>m<sub>b</sub> × e</td><td>Unbalance applied at the LP turbine disk</td><td>0.225 m.kg</td></tr></table>

<table><caption>**Table 5.** Contact formulation properties</caption><tr><td><b>P</b><sub>φ<sub>xc</sub></sub></td><td>Rotation matrix related to φ<sub>xc</sub></td><td><b>[</b> 1 0 0 0 cos φ<sub>xc</sub> sin φ<sub>xc</sub> 0 −sin φ<sub>xc</sub> cos φ<sub>xc</sub> cos φ<sub>yc</sub> 0 −sin φ<sub>yc</sub> 0 1 0 sin φ<sub>yc</sub> 0 cos φ<sub>yc</sub> <b>]</b></td></tr><tr><td><b>P</b><sub>φ<sub>yc</sub></sub></td><td>Rotation matrix related to φ<sub>yc</sub></td><td><b>[</b> cos β 0 −sin β 0 1 0 sin β 0 cos β <b>]</b></td></tr><tr><td><b>P</b><sub>β</sub></td><td>Rotation matrix related to the orientation of the blades cross-section</td><td><b>[</b> cos α<sub>j</sub> sin α<sub>j</sub> 0 −sin α<sub>j</sub> cos α<sub>j</sub> 0 0 0 1 <b>]</b></td></tr><tr><td><b>P</b><sub>α<sub>j</sub></sub></td><td>Rotation matrix related to the j<sup>th</sup> blade position on the disk</td><td><b>[</b> cos α<sub>j</sub> sin α<sub>j</sub> 0 −sin α<sub>j</sub> cos α<sub>j</sub> 0 0 0 1 <b>]</b></td></tr><tr><td><b>P</b><sub>φ<sub>xd</sub></sub></td><td>Rotation matrix related to φ<sub>xd</sub></td><td><b>[</b> 1 0 0 0 cos φ<sub>xd</sub> sin φ<sub>xd</sub> 0 −sin φ<sub>xd</sub> cos φ<sub>xd</sub> cos φ<sub>yd</sub> 0 −sin φ<sub>yd</sub> 0 1 0 sin φ<sub>yd</sub> 0 cos φ<sub>yd</sub> <b>]</b></td></tr><tr><td><b>P</b><sub>φ<sub>yd</sub></sub></td><td>Rotation matrix related to φ<sub>yd</sub></td><td><b>[</b> 1 0 0 0 cos φ<sub>xd</sub> sin φ<sub>xd</sub> 0 −sin φ<sub>xd</sub> cos φ<sub>xd</sub> cos φ<sub>yd</sub> 0 −sin φ<sub>yd</sub> 0 1 0 sin φ<sub>yd</sub> 0 cos φ<sub>yd</sub> <b>]</b></td></tr><tr><td>θ̂</td><td>angular characteristic of the cone aperture</td><td>Industrial value</td></tr><tr><td>R<sub>b</sub></td><td>base radius of the cone representing the casing inner surface</td><td>Industrial value</td></tr><tr><td>h</td><td>height on the cone</td><td>Industrial value</td></tr><tr><td>l<sub>c</sub></td><td>chord length at the blade tip: distance between P(LE) and P(TE)</td><td>Industrial value</td></tr><tr><td>μ</td><td>Friction coefficient</td><td>0.2</td></tr><tr><td>k<sub>r</sub></td><td>Normal contact stiffness</td><td>5 × 10<sup>6</sup> N/m</td></tr><tr><td>g<sub>0</sub></td><td>Initial gap</td><td>0.5 mm</td></tr></table>