# 3D interaction in bladed rotor-to-stator contact

Marie-Océane Parent, Fabrice Thouverez, Fabrice Chevillot

► To cite this version:

Marie-Océane Parent, Fabrice Thouverez, Fabrice Chevillot. 3D interaction in bladed rotor-to-stator contact. EURODYN 2014, Jun 2014, Porto, Portugal. hal-01223447

HAL Id: hal-01223447

https://hal.science/hal-01223447

Submitted on 6 Dec 2016

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
3D interaction in bladed rotor-to-stator contact

Marie-Océane PARENT¹,², Fabrice THOUVEREZ¹, Fabrice CHEVILLOT²

¹LTDS, UMR CNRS 5513, École Centrale de Lyon, 36 avenue Guy de Collongue, 69134 Écully Cedex, France

²Snecma, SAFRAN Group, Site de Villaroche, Rond-Point René Ravaud, 77550 Moissy-Cramayel, France

email: marie-oceane.parent@ec-lyon.fr, fabrice.thouverez@ec-lyon.fr, fabrice.chevillot@snecma.fr

**ABSTRACT:** Reducing the clearance between blade tips and casing is an important factor in increasing the performance of turbomachines such as turbofan engines. However, small clearances lead to more chances of rotor-stator contact, which can cause unstable vibrations. To improve the prediction of such contact-related phenomena, this paper proposes a model that includes potential contacts between the fan blade tips and their surrounding casing. This phenomenological model consists of a fully flexible bladed rotor and a flexible casing. The bladed rotor model is built according to the method developed by Lesaffre [1] and based on the approximations used by Sinha [2]. The casing model is composed of an elastic ring with isotropic suspension. Two different contact formulations have been developed: a 3D contact detection law, which incorporates 3D model kinematics along with 3D local geometry of the contact area, is compared to a simpler formulation. For both formulations, once the clearance has been consumed, contact forces are introduced with normal contact stiffness, and friction is considered using Coulomb's Law. In an initial step, linear harmonic results serve to compare these contact formulations in their ability to detect contact. Next, the stability of balanced static solutions is studied and transient analyses are carried out to confirm the static results. Study results highlight the influence of 3D kinematics and 3D local geometry on both contact detection and the dynamic stability of the system subjected to blade-casing contact.

**KEY WORDS:** rotor dynamics, bladed assembly, nonlinear dynamics, blade-tip/casing contact

## INTRODUCTION

The clearance between rotating and fixed parts is a governing factor of turbomachine performance. By reducing the leakage flows, smaller clearances lead indeed to a more efficient turbomachine. The physical disadvantage however is an evolution of the blade tip-casing contacts capable of generating unstable vibrations. Consequently, in order to reduce clearances, aero-engine mechanical engineers must properly predict the occurrences of contact and ensure that the possible contacts do not trigger damaging phenomena.

Most of the literature reviewed on blade tip-casing contact is divided into two areas of study [3]: bladed assembly vibration and rotor dynamics.

- The study of bladed assembly vibration takes into account the blades and disk flexibility while neglecting shaft dynamics [4], [5]. As part of this area of study, we encountered extended finite element models of the bladed disk and casing [5]. Appropriate contact models were built to handle the discretization due to meshing [6], [7]. Other investigations included the mechanical properties of the abradable coating that covers the inner surface of the casing [8] or the thermal effects of friction [9].

- The field of rotor dynamics focuses on shaft vibrations while neglecting both the blades and disk flexibility. Most numerical investigations were found to use simplified models, whereby the blade tip-casing contact is introduced

using nonlinear springs. Such models typically consider contact between two cylindrical surfaces [10–14].

Beside tip clearance, weight is a recurring governing factor in aircraft engine design. Reducing the mass generally increases the global flexibility and/or eases the dynamic couplings between the various components of the engine. The contact events can thus jointly affect and be influenced by the global engine dynamics. Hence, lighter designs bring out the need to extend the dynamic considered in the prediction model. Moreover, the present work focuses on the contact that may occur between the fan blade tips and fan casing of a turbofan engine. This contact localization underscores the need to consider the fan blades and fan casing dynamics. Modern turbofan designs present indeed a high bypass ratio i.e. a wide fan casing and slender fan blades. These design characteristics imply including the blades, shaft and casing dynamics into the prediction model. The model presented in this paper thus merges both areas of study : bladed assembly vibration and rotor dynamics. It is based on the bladed-shaft models proposed by [1, 2, 15].

The phenomenological model proposed introduces these components while minimizing the number of degrees of freedom (dofs) considered. The reasonable size of the model limits the computation time and allows conducting influence analyses on various parameters. However, a light model necessitates focusing on a proper approximation of the dynamics properties while the geometric aspects of the
<figure><img src="image_2.png" /><figcaption>Figure 1. Model components</figcaption></figure>

components can be locally altered. Since contact is a local event that exerts a strong effect on global dynamics, we propose herein a 3D contact formulation so as to regain accuracy in contact detection without increasing the model size.

Once the model assumptions and the contact formulations described, the results provided in this paper will be exhibited in the final three sections. Linear results will be presented first so as to draw a comparison of the ability to detect contact with the developed contact formulations. Then, stability of the balanced static solutions will be analyzed, before concluding with transient analyses that confirm the static results.

## 1. MODEL DESCRIPTION

The model consists of a flexible bladed rotor and suspended flexible casing. To eliminate time-dependent terms resulting from the periodicity of the bladed disk, the entire model has been developed in the rotating frame.

**Bladed rotor** - The bladed rotor model discussed in the present paper has been based on the same approximations as the one described by Sinha [2]; it has been developed through application of an energy approach by Lesaffre [1] and Gruin [15]. The shaft is modeled by a Euler-Bernoulli beam suspended by isotropic viscoelastic bearings. The shaft beam has a hollow circular cross-section and includes two displacements defined in the rotating frame: the two orthogonal translations in the cross-section plane. Rigid disks connected to the shaft are modeled by a point mass and rotational inertia. Lastly, a set of ten flexible blades are clamped to the representative fan disk; they are all modeled by Euler-Bernoulli beams with a constant rectangular cross-section and can deflect along their more flexible direction. The model sketch with some rough dimensions is drawn in Figure 1.

The shaft and blade deformations are discretized through use of the appropriate Rayleigh-Ritz functions (the detailed equations can be found in [1] and [15]). The gyroscopic effect is taken into account along with the stress stiffening of the blades. Two degrees of freedom (dof) are assumed to model the deformation of each blade. The shaft contains

ten dofs: two rigid-body dofs related to the suspensions, and eight dofs introduced to model its deformation.

**Casing** - The casing model consists of an elastic ring with radial and tangential deformations [16], whose discretization considers the nodal diameter modes of the axisymmetric structure, as described in [1]. We have chosen to include the modes with two and three nodal diameters because for the considered settings, the natural frequencies of these modes lie within the studied range. In addition to its four deformation dofs, the casing contains four rigid-body dofs, contributed by its viscoelastic isotropic suspension: the two orthogonal translations and rotations in the cross-section plane of the shaft.

## 2. CONTACT FORMULATIONS

This section presents the 2D and 3D contact formulations developed. They both follow the same sequence of steps: first, a minimum distance function is developed in order to evaluate the blade-to-casing distance; afterwards, if contact is detected, then the opposing contact and friction forces are applied to the rotating and fixed parts. For both formulations, the normal contact forces are proportional to penetration: the radial stiffness $k_r$ is introduced. Tangential friction forces are obtained using Coulomb's Law: $-\operatorname{sgn}(V_{rel})\mu F_N$ where $\mu$ is the friction coefficient and $F_N$ the normal contact force. The sign of $V_{rel}$, the relative speed between contact surfaces, is assumed to be constant and opposite to the sliding speed direction due to shaft spin. We can thus define a constant tangential stiffness parameter: $k_t = -\operatorname{sgn}(V_{rel})\mu k_r$.

The contact formulations developed differ in their minimum distance function and in the distribution of contact forces on the blade tip. For the 2D formulation, the contact reactions are concentrated on a point-like blade tip, while the 3D formulation allows distributing contact reactions along the chord.

**2D contact formulation** - The 2D formulation considers a point-like blade tip, and its minimum distance function solely takes into account the disk and casing translations and the radial deformation of the casing:

$$g_j^{(2D)} = g_0 + l_b + r_d + u(\alpha_j, t) + (x_c(t) - x_d(t)) \cos \alpha_j + (y_c(t) - y_d(t)) \cos \alpha_j \quad (1)$$

where:

* $\alpha_j$: angular position of the $j^{th}$ blade,

* $g_0$: initial tip clearance,

* $l_b$: blade length,

* $r_d$: fan disk radius,

* $u(\alpha_j, t)$: radial deformation of the casing,

* $x_d, y_d$: disk translations,

* $x_c, y_c$: casing translations.

If $g_j^{(2D)} < 0$, then contact is detected and the normal and tangential contact reactions are given by:

$$\text{if } g_j^{(2D)} < 0 \Rightarrow \begin{cases} F_{jN}^{(2D)} = k_r g_j^{(2D)} \\ F_{jT}^{(2D)} = k_t g_j^{(2D)} \end{cases} \quad (2)$$
<figure><img src="image_3.png" /><figcaption>Figure 2. Contact region geometry</figcaption></figure>

The contact reactions relative to each blade tip j are then projected onto the general coordinates of the model in order to include the contact forces in the system of equations.

**3D formulation** - The 3D minimum distance function takes into account more model kinematics: in addition to the disk and casing translations and normal casing deformation, this function introduces the disk and casing rotations as well as blade flexure. The 3D contact formulation also includes the 3D geometry of the contact region, i.e.: chord length and inclination of the casing inner surface (see Figure. 2). To simplify the formulation, we propose the following kinematic assumptions:

*   Disk and casing rotations are assumed to be sufficiently small so as to linearise their trigonometric functions.

*   The influence of the casing tangential deformation on the minimum distance value is neglected.

*   The influence of casing thickness on deformations along the casing inner surface is neglected: the normal deformation is calculated at the mid-surface.

*   The influence of overall blade tip displacement on the angular position is neglected: the normal casing deformation input into the minimum distance function is measured at the initial blade tip point position.

A number of simplifications to the 3D local geometry have also been considered, namely:

*   The minimum distance function along the chord is obtained by linear interpolation: the gap is only calculated at the Leading Edge (LE) and Trailing Edge (TE) points.(see Figure 3).

*   The inner surface of the casing is considered to be perfectly conical.

To obtain the minimum distance function at both the LE and TE points of the blade tip chord, we begin by determining these positions in a frame attached to the casing. The selected points P(LE) and P(TE) are identified in Figure 2 and 3, with their final position being

<figure><img src="image_4.png" /><figcaption>Figure 3. 3D view of the clearance values interpolated along the chord line</figcaption></figure>

given by:

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

where:

*   P: initial coordinates of the selected point (P(LE) or P(TE)),

*   P': final coordinates of the selected point,

*   z_c: the casing center of gravity position along the z axis,

*   x_d, y_d: disk translations,

*   x_c, y_c: rigid-body casing translations,

*   x_b(l_b, t): blade deflection at its tip,

*   P_β and P_α_j: blade positioning matrices,

*   P_φ_{xd} and P_φ_{yd}: disk rotation matrices,

*   P_φ_{xc} and P_φ_{yc}: rigid-body casing rotation matrices.

Once the P'(LE) and P'(TE) values are known, the minimum distance to the LE and TE points of the blade tip chord is obtained by using the following expression:

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

where:

*   R_b, h: the base radius and height of the cone representing the casing inner surface,

*   θ: angular characteristic of the cone aperture (see Figure 2),

*   P' = {P'_x, P'_y, P'_z}^T: the selected point coordinates defined by (3),

*   α(P): initial angular position of the selected point,
• $u(\alpha(P), t)$: normal deformation of the casing at $\alpha(P)$.

The contact reactions are then obtained by integrating the contact pressure distributed along the chord line. If $g_j^{(3D)}(LE)$ and/or $g_j^{(3D)}(TE)$ is negative, then contact occurs and the normal and tangential pressures distributed along the chord line are given by:

$$ 
\begin{aligned}
p_{jN}(c) &= \frac{k_r}{l_c}(\zeta c + \eta) \quad \text{and} \quad p_{jT}(c) = \frac{k_t}{l_c}(\zeta c + \eta) \\
\text{where } \zeta &= \frac{g_j^{(3D)}(TE) - g_j^{(3D)}(LE)}{l_c} \qquad (5) \\
\text{and } \eta &= \frac{g_j^{(3D)}(TE) + g_j^{(3D)}(LE)}{2}
\end{aligned}
$$

with

• $c$: position along the chord,

• $l_c$: chord length.

Contact reactions are thus given by:

$$ F_{jN}^{(3D)} = \int_{c_1}^{c_2} p_{jN}(c)dc \quad \text{and} \quad F_{jT}^{(3D)} = \int_{c_1}^{c_2} p_{jT}(c)dc \quad (6) $$

Distribution along the chord also creates reaction torques, given by:

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

The limit values for integration ($c_1$ and $c_2$) depend on the minimum distances results:

$$ 
\left\{ 
\begin{array}{l}
g_j^{(3D)}(TE) \le 0 \text{ and } g_j^{(3D)}(LE) \le 0 \Rightarrow c_1 = -\frac{l_c}{2}; c_2 = \frac{l_c}{2} \\
g_j^{(3D)}(TE) \le 0 \text{ and } g_j^{(3D)}(LE) \ge 0 \Rightarrow c_1 = -\frac{l_c}{2}; c_2 = -\frac{\eta}{\zeta} \\
g_j^{(3D)}(TE) \ge 0 \text{ and } g_j^{(3D)}(LE) \le 0 \Rightarrow c_1 = \frac{\eta}{\zeta}; c_2 = \frac{l_c}{2}
\end{array}
\right. 
\qquad (8) $$

As for the 2D formulation, the contact-related reaction forces and torques are then projected onto the general coordinates to create the nonlinear force vector included in the system of equations.

### 3. RESOLUTION METHODS

This section describe the methods used to obtain the results presented in this paper. Three different sets of results will be presented in the following sections: contact detection, permanent contact stability analyses, and nonlinear transient analyses.

**Contact detection** - The contact detection results are obtained by including the linear imbalance response into the minimum distance functions defined above. In the rotating frame, the displacement field of the system subjected to imbalance is obtained by solving the following static equation:

$$ 
\begin{aligned}
\tilde{\mathbf{K}}(\Omega)\mathbf{q} &= \mathbf{F}(\Omega^2) \\
\tilde{\mathbf{K}} &= \mathbf{K} + \mathbf{K}_s(\Omega^2) + \mathbf{N}(\Omega) + \mathbf{C}(\Omega)
\end{aligned}
\qquad (9) $$

where:

• **K**: elastic stiffness matrix,

• $K_s(\Omega^2)$: stress stiffening matrix of the blades,

• $N(\Omega)$: spin softening matrix,

• $C(\Omega)$: viscous damping circulation matrix,

• **q**: general coordinates vector,

• **F**($\Omega^2$): vector of imbalance forces projected onto the general coordinates.

The spin speed $\Omega$ is assumed to be constant at each considered step. The resulting displacements are then introduced into the 2D and 3D minimum distance functions (1) and (4). This method neglects the nonlinear aspect of the contact: contact detection brings indeed contact forces affecting the displacement field and in turn, modifying the contact configuration. The results output must not be used in a refined analysis of the contact configurations. This method does however provide a quick way to obtain contact detection limits: it allows to compare the ability of each formulation to detect contact and identify the parameters that influence the contact detection range.

**Permanent contact stability analyses** - Unlike the previous method, the permanent contact analyses include nonlinear contact forces in the static system of equations. The static solution $\mathbf{q}_0$ is found by solving the following nonlinear system of equations:

$$ \tilde{\mathbf{K}}(\Omega)\mathbf{q}_0 = \mathbf{F}_{nl}(\mathbf{q}_0) + \mathbf{F}(\Omega^2) \qquad (10) $$

where:

• **M**: mass matrix,

• **D**: viscous damping matrix,

• **G**($\Omega$): gyroscopic effects matrix,

• **F**$_{nl}$($\mathbf{q}_0$): nonlinear contact forces vector projected onto the general coordinates of the system.

In order to study the local stability of the solution $\mathbf{q}_0$, we introduce a time-dependant perturbation $\tilde{\mathbf{q}}(t)$. Stability information is obtained by performing an eigen-analysis on the perturbed system: the static solution $\mathbf{q}_0$ being unstable if at least one of its eigen-values has a positive real part. The perturbation is assumed to be sufficiently small to allow the linearisation of contact forces. The linearised system of equations of the perturbed system is given by:

$$ M\ddot{\mathbf{q}}(t) + (\mathbf{D} + \mathbf{G}(\Omega))\dot{\mathbf{q}}(t) + (\tilde{\mathbf{K}}(\Omega) - \mathbf{J}_q\mathbf{F}_{nl}(\mathbf{q}_0))\bar{\mathbf{q}}(t) = \mathbf{0} \qquad (11) $$

where $\mathbf{J}_q\mathbf{F}_{nl}(\mathbf{q}_0)$ is the Jacobian matrix of the contact force vector at $\mathbf{q}_0$.

**Non-linear transient analyses** - The nonlinear transient analyses enable including the nonlinear contact forces in the dynamic system of equations:

$$ M\ddot{\mathbf{q}}(t) + (\mathbf{D} + \mathbf{G}(\Omega))\dot{\mathbf{q}}(t) + \tilde{\mathbf{K}}(\Omega)\mathbf{q}(t) = \mathbf{F}_{nl}(\mathbf{q}(t)) + \mathbf{F}(\Omega^2) \qquad (12) $$

The transient results reported herein were obtained by applying the explicit central differences algorithm, as presented in [5, 17]. The imbalance considered at the fan disk location for inducing contact is applied over a 1 sec ramp.
<figure><img src="image_5.png" /><figcaption>Figure 4. Imbalance response at the fan disk</figcaption></figure>

## 4. CONTACT DETECTION RESULTS

A reasonable imbalance is considered at the fan disk coordinate; the response recorded at the same location is plotted in Figure 5. The peaks observed in the diagram indicate the presence of two critical speeds. Imbalance indeed introduces a harmonic excitation of the same frequency as the shaft's rotation. Consequently, when the rotational frequency equals one of the system's natural frequencies, the corresponding mode will resonate and this spin speed will be referred to as critical speed. The first peak encountered at the normalized spin speed $\Omega \approx 3.5$ displays an amplitude more than five times greater than the peak observed at $\Omega \approx 8$. This finding means that the mode excited at $\Omega \approx 3.5$ contains more deformation energy in the fan module than the mode whose critical speed equals $\Omega \approx 8$. The most critical contact events can thus be expected near $\Omega \approx 3.5$.

By introducing the harmonic response into the minimum distance formulations, we determine whether or not contact would have occurred at the considered spin speed $\Omega$ and for the chosen initial tip clearance $g_0$ (see Section 3). Then, by varying these parameters, we map the contact areas: i.e. the parameter sets ($\Omega, g_0$) for which at least one blade is touching the casing. The maps obtained with the 2D and 3D formulations are plotted in Figure 6: the darker area corresponds to 3D contact detection and the clearer area plotted on top corresponds to 2D contact detection. For $g_0 = 0$ mm, contact occurs for all the considered spin speeds. Then, as the $g_0$ parameter increases, contact is detected for less spin speeds. The upper contours of the maps indicate thus the maximal initial tip clearance for which contact occurs at a given spin speed, i.e. the maximum clearance closure due to imbalance. As expected by the harmonic response at the fan disk, the maximum clearance closure is reached around $\Omega \approx 3.5$ for both formulations. The diagram also shows the efficiency of the 3D minimum distance function with respect to the 2D formulation: the area covered by the 3D formulation is wider, and the difference between clearance values reaches 1.3 mm at the first critical speed.

As described in Section 2, the 3D formulation enriches the minimum distance function on two aspects: the model kinematics and the local geometry of contact. To determine which aspect is more relevant, we considered them individually. Figure 7 presents the difference with

<figure><img src="image_6.png" /><figcaption>Figure 5. Maps of the contact detected with the 2D and 3D minimum distance function, for various spin speed and initial tip clearance values</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Difference between 2D and 3D clearance closures for various configurations of the 3D minimum distance function</figcaption></figure>

the 2D clearance closures for three configurations of 3D minimum distance functions:

* a full 3D configuration as described in Section 2,

* a configuration that includes 3D kinematics and neglects chord length and the conical property of the casing inner surface ($\theta = 0^\circ$ and $l_c = 1$ mm),

* a configuration that includes 3D geometry yet neglects the disk and casing rotations and blade flexure. The results plotted in figure 7 indicate that for the considered parameters, both degraded 3D configurations are of no utility in comparison with the full 3D version. The clearance difference between 2D and full 3D actually amounted to 1.3 mm, while the difference remained less than 0.1 mm for the other 3D configurations. These findings highlight the combined effect of 3D geometry and kinematics: to consider 3D geometry without enriching the kinematics would be useless and vice versa.

Lastly, the contact detection analyses allow comparing the influence of geometric parameters on the occurrences of contact. Figure 8 shows the contact areas according to $\theta$, the inclination of the casing inner surface, while figure 9 presents the contact areas obtained for various chord lengths. When unspecified, the default settings are: $\theta = 5^\circ$ and $l_c = 20$ cm. Both parameters enlarge the
<figure><img src="image_8.png" /><figcaption>Figure 7. Contact detection maps vs the inclination of the casing inner surface: from $\theta = 0^\circ$ to $30^\circ$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Contact detection maps vs chord length at the blade tip: from $l_c = 1$ cm to 30 cm</figcaption></figure>

contact detection area. However, for the parameter ranges considered herein, the conical property of the casing inner surface seems to exert more influence than the blade length on the contact detected using the 3D formulation.

## 5. PERMANENT CONTACT STABILITY ANALYSES

By introducing an asymmetric elastic and inertial coupling, the rotor-stator contacts can create dynamics instabilities [18–20]. Stability then depends on the contact stiffnesses brought to the system vs. its damping properties. The results presented in this section analyze the stability of balanced static solutions, as described in Section 3.

Figure 10 maps the contact configurations resulting from balanced static solutions calculated with the 2D and 3D contact formulations. The outer map contours match the contact areas obtained in the previous section. Within the contact area, the color scale indicates how many of the 10 blades are touching the casing for each balanced static solution. Figure 11 provides stability information on the balanced static solutions. The color scale represents the sum of all positive real parts of the system eigenvalues at each balanced solution. Summing the real parts of each unstable mode offers one way to indicate the degree of destabilization created by the contacts. A small positive real part of the eigenvalue indicates that the instability can be compensated by a reasonable damping supply.

Conversely, various unstable modes with a high real part suggest a strong contact-related asymmetry.

As could be expected, the stability maps globally resemble the contact configuration maps: the greater the number of blades touching the casing, the greater the number of unstable modes with high positive real parts. Some differences in appearance are observed however between the stability results of the two contact formulations. For the 2D configuration, the unstable region covers almost the entire contact detection area, while contact does not necessarily lead to an unstable solution for the 3D configuration. A strip of stable contact configuration surrounds the unstable area in the 3D formulation solutions and another stable strip exists within the peak of the first critical speed. The 2D formulation stability map clearly reveals delimited colored patches whose contours match those of the contact configuration observed in figure 10.

For the 2D contact formulation, each new touching blade contributes a new equivalent spring with radial and tangential stiffness ($k_r$ and $k_t$). For the 3D formulation on the other hand, stability results are smoother due to the linear distribution of contact pressure along the chord. Blade involvement is indeed more gradual for the 3D formulation: contact begins at the LE or TE point of the chord, and the length of chord in contact gradually increases as does the equivalent stiffness input into the system.

Since the contacts are more easily detected with the 3D formulation, the unstable region is, on the whole, wider with this formulation. The stable strip that appears within the peak of the first critical speed corresponds however to an unstable zone in the 2D formulation. This strip indicates a detected area with three contacts for both minimum distance functions. Since the 2D formulation will include three "full" contacts, the contacts detected by the 3D function do not involve an entire chord length and will thus introduce the equivalent "partial" stiffness.

## 6. NONLINEAR TRANSIENT SIMULATIONS

To confirm the accuracy of permanent contact stability results, we carried out transient simulations. The analyses are carried out at $\Omega = 3.3$ and for an initial tip clearance of 0.8 mm. This configuration lies within the strip described in the previous section: the balanced solution is unstable for the 2D formulation and stable for the 3D one. The evolution of normal contact forces at the blade tips ($F_{jN}^{(2D)}$ and $F_{jN}^{(3D)}$) is plotted in figure 12. The 2D blade tip reactions diverge after the full imbalance is applied, whereas the 3D simulation leads to a stable static solution. These results confirm those predicted by the stability analysis on static balanced solutions. For both simulations, three blades are in contact at 1 sec although for the 3D formulation, the third blade touches at 88% of its chord length. This difference is sufficient to maintain stability for the 3D results while the 2D simulation, which introduces a stiffness equivalent to all three blades, is unstable. Lastly, we can remark that at 1 sec, the contact reaction at the blade tip is greater for the 3D formulation than for the 2D.
<figure><img src="image_10.png" /><figcaption>Figure 9. Contact configurations resulting from the static nonlinear resolution</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10. Contact detection and stability results</figcaption></figure>

This finding illustrates the fact that dynamic stability is determined by the asymmetric stiffness link input into the system and not by the contact force amplitude.

## CONCLUSION

A phenomenological model has been developed to investigate phenomena related to blade-to-casing contacts in turbofan aircraft engines. The model developed contains 38 dofs and comprises a flexible bladed shaft and a flexible suspended casing. Two contact formulations have been tested: the first one includes 3D model kinematics and local geometry of the contact region while the second one

<figure><img src="image_12.png" /><figcaption>Figure 11. Contact detection and stability results</figcaption></figure>

is a simpler 2D formulation. Three distinct results were presented:

* The linear results allow comparing the ability of each minimum distance function to detect contacts. For the considered settings, these linear results prove that the 3D minimum distance function is more efficient than the 2D. Other investigations have indicated that the 3D function loses its benefit if one of its 3D attributes is removed, i.e. considering a 3D local geometry without 3D kinematics would be useless and *vice-versa*. Moreover, the analysis of geometric parameter influence indicates that the inclination of the casing inner surface heavily influences contact detection while chord length affects detection results to a lesser extent.

* The stability of a balanced static solution was also studied and the results obtained with each formulation compared. Thanks to overall contact detection efficiency, the range of unstable static solutions is generally broader with the 3D formulation, yet some configurations are stable for 3D solutions and unstable with the 2D solutions. The 3D formulation actually distributes contact pressure along the chord, hence implying a more gradual evolution in the asymmetric stiffness responsible for system destabilization.

* As a final item, nonlinear transient analyses were carried out at a spin speed and initial tip clearance value at which the balanced solution is unstable for the 2D formulation and stable for the 3D. The results generated confirm the stability analysis findings on balanced static solutions.

The 3D contact formulation presented in this paper constitutes an efficient means for improving the predictability of this phenomenological model. The enriched minimum distance function has proven to be more efficient at contact detection and has expanded the parameter set to be investigated. In addition to geometric contact detection, the 3D formulation affects stability results. By allowing the contact to be distributed along the chord, the 3D
formulation more effectively models the real involvement
of the blade tip in the casing. In closing, the proposed
formulation has the significant advantage of offering these
improvements without increasing the model size.

ACKNOWLEDGMENTS

The authors would like to thank the Snecma firm for
its technical and financial support. This work has been
conducted within the framework of the MAIA Mechanical
Research and Technology Program, sponsored by CNRS,
ONERA and the SAFRAN Group.

REFERENCES

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

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

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

[4] Padova, C., Dunn, M., Barton, J., Steen, T., and Turner, K., 2011. "Controlled fan blade tip/shroud rubs at engine conditions". In ASME Turbo Expo, Vancouver, Canada, June 6-10.

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

[6] Arnoult, E., Guilloteau, I., Peseux, B., and Bonini, J., 2000. "A new contact finite element coupled with an analytical search of contact". In European Congress on Computational Methods in Applied Science and Engineering, Barcelona.

[7] Batailly, A., Magnain, B., Legrand, M., and Pierre, C., 2010. "Validation of a 3d contact algorithm for the study of blade-tip/casing contacts in turbomachines". In 8th IFTOMM International Conference on Rotordynamics, Seoul : Korea, Republic Of.

[8] Batailly, A., Legrand, M., and Pierre, C., 2011. "Influence of abradable coating wear mechanical properties on rotor stator interaction". In ASME Turbo Expo, Vancouver, Canada, June 6-10.

[9] Millecamps, A., Brunel, J.-F., Dufrenoy, P., Garcin, F., and Nucci, M., 2009. "Influence of thermal effects during blade-casing contact experiments". *ASME Conference Proceedings, 2009*(48982), pp. 855-862.

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

[11] Childs, D., and Bhattacharya, A., 2007. "Prediction of dry-friction whirl and whip between a rotor and a stator". *Journal of Vibration and Acoustics, Transactions of the ASME*, **129**(3), pp. 355-362.

[12] Popprath, S., and Ecker, H., 2007. "Nonlinear dynamics of a rotor contacting an elastically suspended stator". *Journal of Sound and Vibration*, **308**(3-5), pp. 767–784.

[13] Cole, M., 2008. "On stability of rotordynamic systems with rotorstator contact interaction". *Proceedings of the Royal Society A*, **464**, pp. 3353-3375.

[14] Shang, Z., Jiang, J., and Hong, L., 2011. "The global responses characteristics of a rotor/stator rubbing system with dry friction effects". *Journal of Sound and Vibration*, **330**(10), pp. 2150–2160.

[15] Gruin, M., Thouverez, F., and Blanc, L., 2011. "Nonlinear dynamics of a bladed dual shaft". *European Journal of Computational Mechanics*, **20**, pp. 207-225.

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

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

[18] Crandall, S., 1995. *Nonlinear Dynamics and Stochastic Mechanics*. CRC Press, ch. Rotordynamics, pp. 1–44.

[19] Vangipuram Canchi, S. V., and Parker, R. G., 2005. "Parametric instability of a circular ring subjected to moving springs". *Journal of Sound and Vibration*, **293**, pp. 360-379.

[20] Lesaffre, N., Sinou, J.-J., and Thouverez, F., 2007. "Stability analysis of rotating beams rubbing on an elastic circular structure". *Journal of Sound and Vibration*, **299**(4-5), pp. 1005-1032.