# Investigation of Structural Rotor-Stator Modal Interaction in Turbomachinery

Mathias Legrand, Bernard Peseux, Christophe Pierre

► To cite this version:

Mathias Legrand, Bernard Peseux, Christophe Pierre. Investigation of Structural Rotor-Stator Modal Interaction in Turbomachinery. Sixth European Conference on Structural Dynamics (Eurodyn 2005), Sep 2005, Paris, France. hal-01008366

HAL Id: hal-01008366

https://hal.science/hal-01008366

Submitted on 16 Aug 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
# Investigation of structural rotor-stator modal interaction in turbomachinery

Mathias Legrand & Bernard Peseux

*École Centrale de Nantes, GeM, Équipe Calculs et Structures, Nantes, France*

Christophe Pierre

*Department of Mechanical Engineering, University of Michigan, Ann Arbor, USA*

**ABSTRACT:** In modern turbo machines such as aircraft jet engines, contact between the casing and bladed disk may occur through a variety of mechanisms: coincidence of vibration modes, thermal deformation of the casing, rotor imbalance, etc. These nonlinear interactions may result in severe damage to both structures and it is important to understand the physical mechanisms that cause them and the circumstances under which they occur. In this study, we focus on the phenomenon of interaction caused by modal coincidence. A simple two-dimensional model of the casing and bladed disk structures is introduced in order to predict the occurrence of the interaction phenomenon versus the rotation speed of the rotor. Each structure is represented in terms of its two $n_d$-nodal diameter vibration modes, which are characteristic of axi-symmetric structures and allow for travelling wave motions that may interact through direct contact. The equations of motion are solved using an explicit time integration scheme in conjunction with the Lagrange multiplier method where friction is considered. Results show good agreement between the theory and the numerical tool in the prediction of rotational speed to be avoided.

## 1 INTRODUCTION

One of the main objectives of turbo-machines engineering is to safely increase the efficiency of the machines. This efficiency, simply defined as the ratio of energy output to energy input, depends strongly on the clearance between the rotating and stationary parts: the wider the clearance, the less efficient the machine. The shapes optimization and the reduction of the clearance are critical because they increase the possibility of contact between both structures. This contact has many origins: gyroscopic effect under certain operating conditions, apparition of a rotor imbalance due to the uncertainties of the design or following bird ingestions. Theses interactions can give rise to long lived phenomena characterized by intermittent soft contacts between the structures that induce very high stress levels due to the excitation of the modeshapes of both structures (Arnoult 2000).

The primary goal of the present paper is to give a better understanding of the last mentioned interaction phenomenon between the stator and the rotor in a jet-engine generally called "modal interaction": to the best of our knowledge, one particularly interesting relevant dissertation (Schmiechen 1997) and a paper (Legrand et al. 2003) deal with this travelling wave speed coincidence. This non linear interaction may lead to dev-

astating effects like the destruction of the engine and jeopardizes passengers' safety. It is therefore important to get a full insight of the phenomenon and to give accurate predictions. To this end, a simple numerical tool, based on two elastic structures discretized on their own modeshapes has been developed. Each structure is represented in terms of its two $n_d$-nodal diameter vibration modes which allow for travelling wave motion. The kinetic energy of the rotor is transformed into vibratory energy thanks to direct contact between the two structures and may result in a case of interaction. In this numerical tool, in conjunction with the Lagrange multiplier method where friction is considered, the equations of motions are solved using a time-stepping method based on the explicit central differences scheme more relevant for transitional dynamics dealing with contact constraints (Belytschko et al. 2000).

The paper is organized as follows. In the first section, very simple mathematical statements explain the concept of modal interaction. Then, in the second section, the modelling of the structures used in the numerical tool is presented. The interaction is investigated by an explicit time integration method, in a third part. Finally, the results found are explained and compared to the theory.
## 2 MODAL INTERACTION DEFINITION

The axi-symmetry of a structure, like a bladed disk or a casing of a turbo-machinery illustrated in Fig. 1, which mathematically leads to circulant mass and stiffness matrices, gives rise to mode pairs with identical natural frequencies as explained by Blad¹, ², ³, ⁴

<figure><img src="image_2.png" /><figcaption>Figure 1. 3 nodal-diameter modes of axi-symmetric structures</figcaption></figure>

Both respective modeshapes are similar and rotated around the axis of symmetry by $\pi/(2n_d)$ where $n_d$ is the number of nodal diameters of the modes. In that special case, it is possible to combine these two modes into forward and backward travelling waves. Let $\theta$ be an angular position of a typical material point of the structure and $\alpha$, an amplitude of displacement of this material point, these two rotating waves are simply expressed as:

$$ 
\begin{aligned}
\alpha_f(t, \theta) &= \alpha_0 \cos(\omega t - n_d (\theta - \theta_0)) \\
\alpha_b(t, \theta) &= \alpha_0 \cos(\omega t + n_d (\theta - \theta_0))
\end{aligned}
\qquad (1) $$

where $\omega$ is the pulsation and $\theta_0$ an undetermined parameter of angular position. For a motionless structure, the propagation velocity of these two waves in a stationary frame is $\pm \omega/n_d$, depending on the direction of rotation. For a rotating structure like a bladed disk in a turbo-machine, the angular rotational velocity $\Omega$ (selecting the counter-clockwise direction as positive) has to be added. The propagation velocity of these two waves in a stationary frame becomes then $\Omega + \omega/n_d$ for the forward rotating wave and $\Omega - \omega/n_d$ for the backward rotating wave.

The two natural frequencies are denoted $\omega_c$ for the casing and $\omega_{bd}$ for the bladed disk. Keeping in mind that $\omega_c$ is usually greater than $\omega_{bd}$ for a given number of nodal diameters, if both structures vibrate at their own natural frequency, two cases of travelling wave speed coincidence can be stated for $\Omega$ as illustrated in Fig. 2:

$$ \omega_c = n_d \Omega \pm \omega_{bd} \qquad (2) $$

Taking into account physical considerations on the direction of the contact and friction forces between the two structures - forward on the casing and backward on

<figure><img src="image_3.png" /><figcaption>Figure 2. Potential wave speed coincidence</figcaption></figure>

the bladed disk if the rotor rotates counter-clockwise, only one of these equations can be considered as dangerous:

$$ \omega_c = n_d \Omega - \omega_{bd} \qquad (3) $$

At this speed, both structures are driven in resonance by the contact forces. This rotational speed must be avoided because resonance can cause undesired large amplitude vibrations: it has to be said that these statements are really simple compared to the nonlinear contact phenomena that can take place during a modal interaction and that is why they need to be confirmed by a predicting numerical tool.

## 3 MODELIZATION OF THE STRUCTURES

The design of the structures in the numerical tool we developed is based on the concluding remarks of a previous study (Legrand et al. 2003): in order to sum up, the degrees-of-freedom of the normal and tangential displacement of the bladed disk have to be coupled together. Moreover, the design must enable the modelling of friction forces between the structures.

### 3.1 Casing

The casing is discretized using curved beam finite elements with cubic polynomials in $u$ and $v$ and four degrees-of-freedom per node: $u, u, s, v$ and $v, s$ (where $s$ is the path variable) as depicted in Fig. 3. The shape functions are as follows:

$$ N_1(s) = 1 - \frac{3s^2}{l_e^2} + \frac{2s^3}{l_e^3} \qquad N_2(s) = s - \frac{2s^2}{l_e} + \frac{s^3}{l_e^2} $$

$$ N_3(s) = \frac{3s^2}{l_e^2} - \frac{2s^3}{l_e^3} \qquad N_4(s) = -\frac{s^2}{l_e} + \frac{s^3}{l_e^2} $$

where $l_e$ is the length of a finite element. The discretized displacement field becomes for an element
<figure><img src="image_4.png" /><figcaption>Figure 3. Schematic of the model of the bladed disk and casing</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. 3-nodal diameter modes of the casing and the bladed disk</figcaption></figure>

whose nodes are denoted 1 and 2:

$$ 
\begin{aligned}
u(s) &= N_1 u_1 + N_2 u_{1,s} + N_3 u_2 + N_4 u_{2,s} \\
v(s) &= N_1 v_1 + N_2 v_{1,s} + N_3 v_2 + N_4 v_{2,s}
\end{aligned} 
\quad (4) $$

## 3.2 Bladed disk

The blades are discretized by the usual Euler-Bernoulli beam elements with three degrees of freedom (u, v and θ) per node. In the v direction, the shape functions are the same as those of the casing but are linear along u:

$$ M_1(x) = 1 - \frac{x}{l_e} \quad M_2(x) = \frac{x}{l_e} \quad (5) $$

In a local frame fixed to a rod-bar finite element, the discretized displacement field becomes then:

$$ 
\begin{aligned}
u(x) &= M_1 u_1 + M_2 u_2 \\
v(x) &= N_1 v_1 + N_2 \theta_1 + N_3 v_2 + N_4 \theta_2
\end{aligned} 
\quad (6) $$

Moreover, the blades are initially geometrically curved in order to couple the deflection and normal displacement, respectively v and u. This represents a critical point of the modelization. Finally, a network of stiffnesses connect them together to represent the disk.

Each structure is then reduced to its two n<sub>d</sub>-nodal diameter vibration modes, which are characteristic of axi-symmetric structures and allow for travelling wave motion. As a consequence of the curved design of the blades, the n<sub>d</sub>-nodal diameter modeshape of the casing fits its counterpart of the bladed disk as showed in Fig. 4.

# 4 CONTACT DYNAMICS

## 4.1 General theory

The forces of particular interest in this study are the forces of contact acting between the blade tips and the

casing where friction is considered. Equations of motion are derived using the Principle of Virtual Work following the procedure described in Laursen (2002). The so-called impenetrability condition is enforced by the Lagrange multiplier method. In a very general way, it is necessary to construct a field of admissible contact pressures C<sub>N</sub>:

$$ C_N = \{t_N : \Gamma_c \to \mathbb{R} \setminus t_N \ge 0\} \quad (7) $$

and another one, respective to the admissible friction tractions C<sub>T</sub>(t<sub>N</sub>):

$$ C_T(t_N) = \{t_T : \Gamma_c \to \mathbb{R}^3 \setminus t_T n = 0, \|t_T\| \le \mu t_N\} \quad (8) $$

where n stands for the normal direction to the contact surface. The formulation of the contact problem can be written in the following manner: find the displacement field u such as for all admissible virtual displacement δu:

$$ 
\begin{aligned}
\int_{\Omega} \rho \ddot{u} \delta u \, dV + \int_{\Omega} \bar{\bar{\sigma}} : \delta \bar{\bar{\epsilon}} \, dV = \int_{\Gamma_{\sigma}} t_d \delta u \, dS \\
+ \int_{\Gamma_c} t_N \delta g + t_T \delta u_T \, dS + \int_{\Omega} f_d \delta u \, dV
\end{aligned} 
\quad (9) $$

where t<sub>N</sub> and t<sub>T</sub> are constrained by conditions (7) and (8). This formulation is discretized in time using the explicit central differences scheme. Keeping in mind notations of equation (9) and denoting the time step by h, vectors in velocity and acceleration become:

$$ 
\left\{ 
\begin{aligned}
\ddot{u}_n &= \frac{u_{n+1} - 2u_n + u_{n-1}}{h^2} \\
\dot{u}_n &= \frac{u_{n+1} - u_{n-1}}{2h}
\end{aligned}
\right. 
\quad (10) $$

This scheme is stable for a very small time step bounded by the Courant Criterion in linear cases. In non linear cases such as those involving contact, an empiric time-step, smaller than for linear cases, has to be chosen.
## 4.2 Algorithm for modal interaction

Many different algorithms dealing with contact problems in explicit dynamics have been developed in the past few years depending on the correction used (displacement or velocity) as mentioned by Vola et al. (1998). Because of its simplicity, the Forward Increment Lagrange Method developed by Carpenter et al. (1991) is used in our study. In what follows, **M**, **D** and **K** represent respectively the mass, damping and stiffness matrices coming from the discretization of the displacement field of the two structures described in section 3. Their 4 × 4 (2 × 2 for each structure) modal counterparts calculated from the two $n_d$-nodal diameter modes are denoted by **M**<sub>m</sub>, **D**<sub>m</sub> and **K**<sub>m</sub>. For instance the modal mass matrix of the casing, **M**<sub>m</sub><sup>c</sup>, is calculated in the following manner:

$$ \mathbf{M}_{\mathrm{m}}^{c}=\mathbf{P}_{\mathrm{c}}^{T} \mathbf{M}^{c} \mathbf{P}_{\mathrm{c}} \qquad (11) $$

where the two columns of **P**<sub>c</sub> are the two $n_d$-nodal diameter modes of the casing. All these matrices are necessary, the first ones in order to take into account the contact constraints and the second ones for the prediction of a *modal* interaction. Our algorithm is divided into three steps and is fairly simple:

1. prediction of the modal displacement **u**<sub>m</sub> of the structures without considering any contact. This predicted displacement, denoted with a superscript <sup>p</sup>, is analytically expressed as:

$$ \mathbf{u}_{\mathrm{m}}^{n+1, p}=\left[\frac{\mathbf{M}_{\mathrm{m}}}{h^{2}}+\frac{\mathbf{D}_{\mathrm{m}}}{2 h}\right]^{-1}\left(\left(\frac{2 \mathbf{M}_{\mathrm{m}}}{h^{2}}-\mathbf{K}_{\mathrm{m}}\right) \mathbf{u}_{\mathrm{m}}^{n}+\left(\frac{\mathbf{D}_{\mathrm{m}}}{2 h}-\frac{\mathbf{M}_{\mathrm{m}}}{h^{2}}\right) \mathbf{u}_{\mathrm{m}}^{n-1}+\mathbf{F}_{\mathrm{m}}^{\mathrm{ext}, n}\right) \quad (12) $$

2. projection of the predicted modal displacements onto the physical coordinates system:

$$ \mathbf{u}^{n+1,p}=\mathbf{P} \mathbf{u}_{\mathrm{m}}^{n+1,p} \qquad (13) $$

and determination of the gap function vector **d**<sup>p</sup> between the two structures. Each gap function where a penetration has been detected is kept in **d**<sup>p</sup>, all other coordinates of the latter being zero.

3. correction of the displacements through the calculation of the Lagrange multipliers in order to avoid any penetration. This step implies that the gap functions (linearized when necessary) vanish:

$$ \mathbf{d}^{n+1}=\mathbf{C}_{\mathrm{N}}^{T} \mathbf{u}^{n+1, c}+\mathbf{d}^{p}=\mathbf{0} \qquad (14) $$

where the superscript <sup>c</sup> means that the correction of the displacements is being calculated. **C**<sub>N</sub> is

the contact constraint matrix in the normal direction. The new equations of motion taking the contact forces into account and the contact constraints have to be solved simultaneously. To this end, a new contact matrix **C**<sub>NT</sub> containing the normal and tangential constraints is built by considering that, because of high relative velocities between the two structures, only sliding occurs:

$$ \left\{ \begin{array}{l} \boldsymbol{\lambda}=\left(\mathbf{C}_{\mathrm{N}}^{T}\left[\frac{\mathbf{M}}{h^{2}}+\frac{\mathbf{D}}{2 h}\right]^{-1} \mathbf{C}_{\mathrm{NT}}\right)^{-1} \mathbf{d}^{p} \\ \mathbf{u}^{n+1}=\mathbf{u}^{n+1, p}+\left[\frac{\mathbf{M}}{h^{2}}+\frac{\mathbf{D}}{2 h}\right]^{-1} \mathbf{C}_{\mathrm{NT}} \boldsymbol{\lambda} \end{array} \right. \qquad (15) $$

At the end of this step, the physical displacements verifying the contact constraints **u**<sup>n+1</sup> are projected onto the modal coordinates.

Finally, time is incremented before going back to the first step until final time is reached.

# 5 RESULTS

Table 1 sums up the characteristics of our model. The mechanical parameters of the structures have been chosen so that the eigenfrequencies of the casing are greater than those of the bladed disk in agreement with Fig. 2. The contact between the rotor and the stator is initiated

<table><caption>Table 1. Characteristics of the model</caption><thead><tr><th></th><th>casing</th><th>bladed disk</th></tr></thead><tbody><tr><td>Y. modulus (MPa)</td><td>E<sub>c</sub> = 7 · 10<sup>7</sup></td><td>E<sub>bd</sub> = 2.1 · 10<sup>11</sup></td></tr><tr><td>density (kg·m<sup>-3</sup>)</td><td>ρ<sub>c</sub> = 2800</td><td>ρ<sub>bd</sub> = 7800</td></tr><tr><td>height (m)</td><td>h<sub>c</sub> = 0.01</td><td>h<sub>bd</sub> = 0.01</td></tr><tr><td>width (m)</td><td>l<sub>c</sub> = 0.1</td><td>l<sub>bd</sub> = 0.1</td></tr><tr><td>radius (m)</td><td>r<sub>c</sub> = 0.4895</td><td>r<sub>bd</sub> = 0.4795</td></tr><tr><td>d.o.f.</td><td>160</td><td>660</td></tr><tr><td>modal damping</td><td>ξ<sub>c</sub> = 0.03</td><td>ξ<sub>bd</sub> = 0.005</td></tr><tr><td>frequency (Hz) (n<sub>d</sub> = 3)</td><td>ω<sub>c</sub> = 457.3</td><td>ω<sub>bd</sub> = 447.4</td></tr></tbody></table>

by a 100μs forcing pulse on the first mode of the casing at the very beginning of the simulation.

## 5.1 Contact algorithm and time step size

It is well-known (Wriggers 2002) that strong non linearities like contact play a important role on the size of the time-step of an explicit algorithm. The accuracy of the results in our study depends strongly on this size due to the linearization of the gap functions. Moreover, the two degrees-of-freedom structural discretization in modal coordinates forces the casing (and bladed disk)
wave's propagation speed of a local deformation to become infinite: the contact detection is then never precise enough. Decreasing the time step would give different quantitative results. In this work, $h = 10^{-6}$s has been kept as a time step in order to get significant results while minimizing CPU consumption. A focus on the distances between the blade tips and the casing depicted in Fig. 5 confirms that the correction is successful for the chosen time-step.

<figure><img src="image_6.png" /><figcaption>Figure 5. Predicted ($d^c$) and corrected ($d^p$) distances of penetration between blade 1 and the casing</figcaption></figure>

## 5.2 Modal interaction phenomenon

A study of the modal interaction has been carried out and the response of the two structures is analyzed with respect to two parameters: the rotational speed of the engine and the number of nodal diameters of the modes. For a given number of nodal diameters, the coupling by means of structural contacts between the two main components of our aircraft engine model exhibits three different types of response. They are depicted below for $n_d = 3$ and $\Omega_c$, solution of equation (3) such as:

$$ \Omega_c = 301.56 \text{ rad/s} \qquad (16) $$

Fig. 6 shows that there is no dangerous interaction in terms of modal displacement if $\Omega < \Omega_c$. This response is characterized by several impacts between the two structures at the beginning, followed by a decrease of the amplitudes of vibration to zero due to structural damping. In contrast, in Fig. 8, for which $\Omega > \Omega_c$, amplitudes of vibration of the casing become extremely large very quickly. Just in-between, an almost-periodic response is found (see Fig. 7). In that case, the modal shapes of the both structures fit together and structural contacts with friction allow for a forward rotating mode on the casing and a backward rotating mode on the bladed disk to occur. The phase shift between the two

<figure><img src="image_7.png" /><figcaption>Figure 6. Vibrations of the two modes of the casing for $\Omega = 250 \text{ rad/s}$ ($\Omega < \Omega_c$)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Vibrations of the two modes of the casing for $\Omega = 310 \text{ rad/s}$ ($\Omega \approx \Omega_c$)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Vibrations of the two modes of the casing for $\Omega = 400 \text{ rad/s}$ ($\Omega > \Omega_c$)</figcaption></figure>

curves depicting both modes time evolution is an illustration of a wave propagation whose mathematical state-
ment is given by equation (1).

Calculations for $n_d = 4$ and $n_d = 5$ have been carried out in order to ensure that our algorithm was predictive. Respective results are depicted in Fig. 9. This

<figure><img src="image_10.png" /><figcaption>Figure 9. Critical rotational speeds versus $n_d$ with on the left, the theoretical value and on the right, the numerical value</figcaption></figure>

figure shows good agreement between the theory and the numerical tool, even if the latter predicts a range of rotational speeds for which the interaction becomes unstable:

$$ \Omega \geq \frac{\omega_c + \omega_{bd}}{n_d} \quad (17) $$

Instead, equation (3) predicts only one dangerous rotational velocity which is the lower limit of those given by the numerical tool.

## 6 CONCLUSION

In this paper, we presented an evolution of a previous modal interaction modelling. In order to get interesting and more realistic results, in addition to friction, it seems necessary to couple normal and tangential displacements of both structures. To this end, blades have been geometrically curved. The theory and the numerical tool are in good agreement for the prediction of a modal interaction rotational speed. Nevertheless, our time-marching procedure, dealing with high nonlinear terms like contact constraints, is still very sensitive to the time-step size. In order to restrict this sensitivity and due to the almost periodic behavior of the structural set, an alternating time-multi frequency domain method is being developed. It is based on a general implementation of the harmonic balance method (Pusenjak and Oblak 2004).

## ACKNOWLEDGEMENTS

This research was proposed and supported by SNECMA.

## REFERENCES

Arnoult, E. 2000. Modélisation numérique et approche expérimentale du contact en dynamique: application au contact aubes/carter de turboréacteur. Ph. D. thesis: Université de Nantes.

Belytschko, T., Liu, W.K., & Moran, B. 2000. *Nonlinear Finite Elements for Continua and Structures*. John Wiley & Sons.

Bladh, R. 2000. *Efficient predictions of the vibratory response of mistuned bladed disk by reduced order modeling*. Ph. D. thesis: University of Michigan.

Carpenter, N., Taylor, R., & Katona, M. 1991. Lagrange constraints for transient finite element surface contact. *International Journal for Numerical Methods in Engineering 32* 103–128.

Laursen, T.A. 2002. *Computational contact and impact mechanics*. Springer.

Legrand, M., Peseux, B., & Pierre, C. 2003. Amélioration de la prédiction de l'interaction rotor/stator dans un moteur d'avion. *Sixième Colloque National en Calcul des Structures*: Giens, France.

Pusenjak, R.R. & Oblak, M.M. 2004. Incremental harmonic balance method with multiple time variables for dynamical systems with cubic non-linearities. *International Journal for Numerical Methods in Engineering 59* 255–292.

Schmiechen, P. 1997. *Travelling Wave Speed Coincidence*. Ph. D. thesis: Imperial College of Science, Technology and Medicine - University of London.

Vola, D., Pratt, E., Raous, M., & Jean, M. 1998. Consistent time discretization for a dynamical contact problem and complementarity techniques. *Revue Européenne des Éléments Finis* 7 149–162.

Wriggers, P. 2002. *Computational contact mechanics*. John Wiley & Sons.