Mathias Legrand Gem

S Ébastien Roques
email: sebastien.roques@ec-nantes.fr

Bernard Peseux
email: berbard.peseux@ec-nantes.fr

Christophe Pierre
email: pierre@umich.edu

AIRCRAFT ENGINE STRUCTURAL ROTOR-STATOR MODAL INTERACTION

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 k-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 of the numerical tool and theory show good agreement in the prediction of rotational speed to be avoided. To conclude, the mathematical statements of a multifrequency domain-method are proposed. This method is to be used to circumvent numerical issues inherent to time-marching procedures.

INTRODUCTION

One of the main objectives of turbo-machine engineering is to safely increase the efficiency of the machines. This efficiency, simply defined as the ratio of energy output to energy input, strongly depends on the clearance between the rotating and stationary parts: the wider the clearance, the less efficient the machine. Shape optimization and reduction of the clearance are critical because they increase the possibility of contact between the two 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 [START_REF] Arnoult | Modélisation numérique et approche expérimentale du contact en dynamique: application au contact aubes/carter de turboréacteur[END_REF].

The primary goal of the present paper is to give a better understanding of the last mentioned rotor-to-stator interaction phenomenon generally called 'modal interaction': to the best of our knowledge, one particularly interesting relevant dissertation [START_REF] Schmiechen | Travelling wave speed coincidence[END_REF] and a paper [START_REF] Legrand | Amélioration de la prédiction de l'interaction rotor/stator dans un moteur d'avion[END_REF] deal with this travelling wave speed coincidence. This non linear interaction may lead to devastating 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 numerical tool has been developed. Each structure is represented in terms of its two n dnodal diameter vibration modes which allow for travelling wave motion. The kinetic energy of the rotor is transformed into vibra-1 tory 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 timestepping method based on the explicit central differences scheme more relevant for transitional dynamics dealing with contact constraints [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF].

The paper is organized as follows. The first section explains the concept of modal interaction in simple mathematical statements. In a second section, the modelling of the structures used in the numerical tool is presented. The interaction in itself is investigated by an explicit time integration method in a third part and the respective results are compared to the theory.

Finally, a frequency-domain method is proposed to study extensively the particular behavior of such a mechanical system: this requires an extension of the usual harmonic balance method and a particular attention is given to the mathematical tools needed for the implementation.

MODAL INTERACTION DEFINITION Cyclic symmetry

The cyclic-symmetry of a turbo-machine bladed disk illustrated in Fig. 1 mathematically leads to circulant mass and stiffness matrices whose associated modal analysis statement in physical coordinates u is written as [START_REF] Bladh | Efficient predictions of the vibratory response of mistuned bladed disks by reduced order modeling[END_REF]:

Yu = 0 (1)
with:

Y = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Y 0 Y 1 0 ••• 0 Y 1 T Y 1 T Y 0 Y 1 0 ••• 0 . . . . . . . . . . . . . . . . . . 0 ••• 0 Y 1 T Y 0 Y 1 Y 1 0 ••• 0 Y 1 T Y 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (2)
In Eqn. [START_REF] Schmiechen | Travelling wave speed coincidence[END_REF], the number of blocks on the diagonal is equal to the number of sectors, N, and

Y i = K i -ω 2 M i for i = [0, 1].
K 0 and M 0 come from the mass and stiffness matrices of a fundamental sector. M 1 et K 1 are built by considering displacement compatibility between adjacent sectors. In the context of symmetrical components and based on block-circulant matrix theory, it can be shown that the initial eigenvalue problem (1) reduces to a series of smaller eigenvalue problems mostly depending on a parameter k in real-cyclic coordinates. Each block associated with this parameter will have the following partitions:

Y k = Y 0 + Y 2 cos kα Y 3 sin kα -Y 3 T sin kα Y 0 + Y 2 cos kα ( 3
)
where

Y 2 = Y 1 + Y 1 T and Y 3 = Y 1 -Y 1 T .
α stands for the fundamental interblade phase shift defined as 2π/N. By noting x, the counterpart of the physical displacement vector u in cyclic coordinates, the series of new eigenvalue problems becomes (k = 0 is not of interest):

Y k x k = 0, k = 1,...,K with K = ⎧ ⎪ ⎨ ⎪ ⎩ N -1 2 if N is odd N -2 2 if N is even (4)
The skew-symmetric structure of Y k gives rise to mode pairs with identical natural frequencies [START_REF] Meirovitch | Analytical Methods in Vibrations[END_REF]. By partitioning the eigenvector matrix such as x k = [x c,k , x s,k ] T , both respective modeshapes take the form:

x k ω,1 = x c,k x s,k ω,1 = x c,k ω,1 x s,k ω,1 x k ω,2 = x c,k x s,k ω,2 = ± -x s,k ω,1 x c,k ω,1 (5)
They are similar and rotated around the axis of symmetry by π/(2k). Once the modeshapes are known in real-cyclic coordinates, it is possible to come back to the physical coordinates for the n th sector thanks to the relation:

u n = 2 N K ∑ k=1
x c,k cos n -1)kα + x s,k sin(n -1)kα

+ 1 √ N x 0 + (-1) n-1 √ N x N 2 (6)
In Eqn. [START_REF] Meirovitch | Analytical Methods in Vibrations[END_REF], the last term only exists if there is an even number of sectors. In this equation, the number of nodal diameters, denoted by k for a given linear mode, that are characteristic of such structures and depicted in Fig. 1, is explained. Both shapes of sector n along the two consecutive k-nodal diameter modes are given in

ω c k ω bd k - ω c k - ω bd k Ω wave speed in a stationary frame condition of interaction Ω - ω bd k Ω + ω bd k Figure 2. THEORETICAL WAVE SPEED COINCIDENCE DIAGRAM
physical coordinates by combining Eqns. ( 5) and ( 6) :

u 1 n = x c,k cos(n -1)kα -x s,k sin(n -1)kα u 2 n = x s,k cos(n -1)kα + x c,k sin(n -1)kα (7)
It is then possible to combine these two modes into forward ( f ) and backward ( b ) travelling waves at pulsation ω:

u f n (t) = u 1 n cosωt + u 2 n sin ωt u b n (t) = u 1 n cos ωt -u 2 n sin ωt (8)

Axi-symmetry

Structures featuring perfect axi-symmetry, such as our casing, can be considered as limit cases of structures that exhibit cyclic-symmetry: the number of sectors tends to infinity. Consequently, they share the same properties and have double harmonic modes as well [START_REF] Leissa | Vibration of Shells[END_REF]. These mathematical characteristics are intrinsically linked to the structural invariance with respect to the revolution axis. Any modeshape can be subjected to a rotation and is necessarily a linear combination of two orthogonal modes with identical shape and frequency.

Interaction phenomenon: basic theory

For a motionless structure, the velocity of these two waves in a stationary frame is ±ω/k, depending on the direction of rotation. For a rotating structure like the bladed disk, the rotational velocity Ω (selecting the counter-clockwise direction as positive) has to be added. The propagation velocity of these two waves in a stationary frame becomes then Ω+ ω/k for the forward rotating wave and Ω -ω/k for the backward rotating wave.

The two natural frequencies are denoted ω c for the casing and ω bd for the bladed disk. We keep in mind that for a given number of nodal diameters, ω c is usually greater than ω bd for the first flexion mode of the blades. If both structures vibrate at their own natural frequency, two cases of travelling wave speed coincidence can be stated for Ω as illustrated in Fig. 2:

ω c = kΩ ± ω bd (9)
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 the bladed disk if the rotor rotates counter-clockwise -, only one of these equations can be considered as dangerous:

ω c = kΩ -ω bd ( 10
)
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.

STRUCTURAL MODELS

The design of the structures in the numerical tool we developed is based on the concluding remarks of a previous study [START_REF] Legrand | Amélioration de la prédiction de l'interaction rotor/stator dans un moteur d'avion[END_REF]: normal and tangential displacements of the blade tips have to be coupled together and the design must enable the modelling of friction forces between the structures.

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:

N 1 (s) = 1 - 3s 2 l e 2 + 2s 3 l e 3 N 2 (s) = s - 2s 2 l e + s 3 l e 2 N 3 (s) = 3s 2 l e 2 - 2s 3 l e 3 N 4 (s) = - s 2 l e + s 3 l e 2 ( 11
)
where l e is the length of a finite element. The discretized displacement field becomes for an element whose nodes are denoted 1 and 2:

u(s) = N 1 (s)u 1 + N 2 (s)u 1,s + N 3 (s)u 2 + N 4 (s)u 2,s v(s) = N 1 (s)v 1 + N 2 (s)v 1,s + N 3 (s)v 2 + N 4 (s)v 2,s (12)

Bladed disk

The blades are discretized by the usual Euler-Bernouilli 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 - x l e M 2 (x) = x l e ( 13
)
In a local frame fixed to a beam-bar finite element, the discretized displacement field becomes then:

u(x) = M 1 (x)u 1 + M 2 (x)u 2 v(x) = N 1 (x)v 1 + N 2 (x)θ 1 + N 3 (x)v 2 + N 4 (x)θ 2 (14)
The coupling between deflection and normal displacements, respectively v and u, at the blade tip location is achieved by curving the geometry of the blade. This represents a critical point of the structural modelling. Finally, a network of stiffnesses connect them together to represent the disk.

Modal reduction

In order to focus on a modal interaction between the two structures, the latter are reduced to their respective two k-nodal diameter vibration modes. As a consequence of the curved design of the blades, a k-nodal diameter modeshape of the casing fits its counterpart of the bladed disk as showed in Fig. 4.

In what follows, M, D and K represent respectively the mass, damping and stiffness matrices coming from the discretization in space. Their 4 × 4 (2 × 2 for each structure) modal counterparts are denoted by M m , D m and K m and related to each other by: M m = P T MP with P = P c 0 0 P bd (15) where P c (resp. P bd ) is constituted of the two retained n d -nodal diameter modes of the casing (resp. bladed disk).

CONTACT DYNAMICS

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 within the kinematically linear framework following, in essence, the procedure described in [START_REF] Laursen | Computational contact and impact mechanics[END_REF].

It is first convenient to arbitrary choose one surface subjected to contact as the master one so that the second one, commonly called the slave surface, can be parameterized. It is then possible to find for any material point x belonging to the master surface Γ (m) c , its closest counterpart ȳ on the slave surface Γ (s)

c : ȳ = arg min y∈Γ (s) c x -y (16)
According to these notations, the gap function between the two structures can be stated:

g(x) = g 0 (x) + u (m) (x) -u (s) (ȳ(x)) • n ( 17
)
where g 0 (x) represents the initial positive gap between the two structures and n, the outward normal to Γ

t N ≥ 0, g ≥ 0, t N g = 0 ( 18
)
where t N stands for the contact pressure, assumed positive, acting on Γ (s) c . To these unilateral contact conditions, we add the well-known Coulomb friction law:

t T ≤ µt N t T < µt N ⇒ v T = 0 t T = µt N ⇒ ∃α > 0 such as v T = α t T t T (19)
for which µ is the coefficient of friction, v T , the tangential slip and t T , the tangential stress vector.

In this work, the so-called impenetrability condition is enforced by the Lagrange multiplier method. In a very general way, it is necessary to construct an admissible contact pressure field C N and its counterpart relative to friction tractions, C T (t N ):

C N = t N : Γ (s) c → R \ t N ≥ 0 C T (t N ) = t T : Γ (s) c → R 3 \ t T • n = 0, t T ≤ µt N (20)
An alternative equivalent form to Eqs. ( 18) is:

∀s N ∈ C N , Γ (s) c g(s N -t N ) dS ≥ 0 ∀s T ∈ C T (t N ), Γ (s) c v T • (s T -t T ) dS ≥ 0 (21)
The weak form of the contact problem can be written in the following manner:

find the displacement field u such as for all admissible virtual displacement δ δ δu:

Ω ρ ü • δ δ δu dV + Ω σ σ σ : δ δ δ ε ε ε dV = Γ σ t d • δ δ δu dS + Γ (s) c (t N δg + t T • δ δ δu T ) dS + Ω f d • δ δ δu dV ( 22
)
where t N and t T are constrained by conditions (21). This formulation is first discretized in space using the shape functions described above for u and δ δ δu. Eqs. (21) do not need to be discretized because we are dealing with 'node-to-line' contact.

The obtained matrix form of ( 22) is solved using the explicit central differences scheme. Denoting the time step by h, vectors in velocity and acceleration become:

⎧ ⎪ ⎨ ⎪ ⎩ ün = u n+1 -2 u n + u n-1 h 2 un = u n+1 -u n-1 2h (23)

Solution algorithm for modal interaction

Many different time-marching procedures dealing with contact have been developed in the past few years depending on the type of correction used (displacement or velocity) [START_REF] Vola | Consistent time discretization for a dynamical contact problem and complementarity techniques[END_REF]. Because of its simplicity, the Forward Increment Lagrange Method [START_REF] Carpenter | Lagrange constraints for transcient finite element surface contact[END_REF] is used in our study. The algorithm is divided into three steps:

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

u n+1,p m = M m h 2 + D m 2h -1 2M m h 2 -K m u n m + D m 2h - M m h 2 u n-1 m + F ext,n m (24)
2. projection of the predicted modal displacements onto the physical coordinates:

u n+1,p = P u n+1,p m ( 25
)
and determination of the gap function vector d p between the two structures. Each gap function where a penetration has been detected is kept in d p , all other coordinates of the latter being zero.

3. correction of the displacements through the calculation of the Lagrange multipliers. This step implies that the gap functions (linearized when necessary) vanish:

d n+1 = C N T u n+1,c + d p = 0 (26)
where the superscript c means that the correction of the displacements is being calculated. C N 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 NT containing the normal and tangential constraints is built by considering that, because of high relative velocities between the two structures, only sliding occurs:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ λ λ λ = C N T M h 2 + D 2h -1 C NT -1 d p u n+1 = u n+1,p + M h 2 + D 2h -1 C NT λ λ λ (27)
At the end of this step, the physical displacements verifying the contact constraints u n+1 are projected onto the modal coordinates.

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

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 by a 100 µs forcing pulse on the first mode of the casing at the very beginning of the simulation.

Contact constraints and time-step size study

It is well-known [START_REF] Wriggers | Computational contact mechanics[END_REF] that non linearities like contact strongly constraint the time-step size of an explicit algorithm. Moreover, the two degrees-of-freedom modal discretization forces any wave's propagation speed of a local structural deformation to become infinite: the contact detection is then never precise enough and a mathematical singularity is arising. Decreasing the time-step gives different quantitative results as illustrated in Fig. 5.

Anyways, in this work, h = 10 -6 s has been kept 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.

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 for k = 4:

1. Ω Ω Ω < Ω Ω Ω c c c : the response is characterized by several impacts between the two structures at the beginning of the simulation due to the initial pulse on the casing. They are then followed by a decrease of the amplitudes of vibration to zero stemming from the structural damping as shown in Fig. 6; 2. Ω Ω Ω Ω Ω Ω c c c : 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. Both structures vibrate around their own natural frequency. The phase shift between the two curves in Fig. 7 depicting the time evolution of the two

×10 -5 ×10 -4
Ω = 350.2 RAD/S (Ω < Ω c ): MODE 1 ( ), MODE 2 ( )
modes of the casing is an illustration of a wave propagation whose mathematical statement is given by equation ( 8). 3. Ω Ω Ω > Ω Ω Ω c c c : the behavior of the system is unstable as pictured in Fig. 8 for the casing. The amplitudes of structural vibration become extremely large after about a hundred of rotor rounds.

Theoretical (Ω ct solution of Eq. ( 10)) and numerical (Ω cn ) values of the critical rotational velocity Ω c can be compared for k = 4:

Ω ct = 343.6 rad/s Ω cn = 351.4 rad/s (28)
Calculations for k = 3 and k = 5 have been carried out in order to ensure that our algorithm was predictive. Respective results are depicted in Fig. 9. This figure shows good agreement be- tween the theory and the numerical tool although theoretical and numerical values of Ω c are not identical. This can be explained by the structural damping that is considered in the numerical tool whereas not taken into account in the theory. Moreover, the theory predicts only one dangerous rotational velocity whereas the numerical tool predicts a range rotational velocities for which the interaction becomes unstable:

Ω c ≥ Ω cn Ω ct = ω c + ω bd k (29)
To conclude, a rigorous parametric study is rendered difficult because of time consumption and time-step sensitivity that represent the two main drawbacks of our method. In order to circumvent the latter and to better characterize and understand the temporal results described above, a frequency-domain approach is used. This method is still under development and only the theory will be presented hereinafter: the aperiodic aspect of the motion implies a generalization of the usual harmonic-balance method.

k = 3 Ω c = 301.56 k = 4 Ω c = 343.58 k = 5 Ω c =

N-DIMENSIONAL HARMONIC BALANCE METHOD Preliminary discussions

The application of the usual harmonic balance method is limited to predicting periodic behaviors only. Around a critical rotational velocity highlighted above (Ω Ω c ), both structures vibrate around their own natural frequencies even when any external excitation has disappeared and the system exhibits a complicated aperiodic behavior illustrated in Fig. 10 and made up of two incommensurable frequencies, ω c and ω bd in addition to Ω. In [START_REF] Lau | Incremental harmonic balance method with multiple time scales for aperiodic vibration of nonlinear systems[END_REF], the solution being sought is expanded into a multiple Fourier series containing a number of frequencies that are incommensurable with one another. Theoretically, there is no limitation on the number of time scales used in the n-HBM, and hence a very general aperiodic solution of any accuracy can be obtained if the proper frequencies and a sufficient number of corresponding harmonic terms are included. This motivates, for several reasons, the use of this method to capture the behavior of such a system:

1. the response and the non linear restoring forces are expanded into n-dimensional Fourier series: time integration has shown (almost) pure harmonic responses. 2. n-HBM is not sensitive to the initial conditions as well as the time-step size which represent a critical issue in the timemarching procedure: the system is studied intrinsically. 3. n-HBM transforms a non linear differential system of equations in a non linear algebraic system of equations that can be faster to solve in order to get the response. 4. n-HBM allows for fast parametric and stability studies.

Mathematical statements

In a general manner, the governing equations of any autonomous system take the following form:

Mü + Du + Ku + F nl (u, u) = 0 (30)
First, true time t is adimensionalized with respect to the M s frequencies ω m of the response [START_REF] Pusenjak | Incremental harmonic balance method with multiple time variables for dynamical systems with cubic non-linearities[END_REF]:

τ m = ω m t for m = 1, 2,...,M s (31)
so as to construct an hyper-time domain:

0 ≤ τ m ≤ 2π (32)
Assuming that the aperiodic solution for the generalized coordinate u j ( j th coordinate of u) exists, the questioned solution is sought in the form of truncated multiple Fourier series:

u j = N h ∑ j 1 =-N h N h ∑ j 2 =-N h ••• N h ∑ j Ms =-N h a j j 1 , j 2 ,..., j Ms cos M s ∑ m=1 j m τ m +b j j 1 , j 2 ,..., j Ms sin M s ∑ m=1 j m τ m ( 33
)
where N h denotes the highest order of the harmonic terms in the truncated Fourier series. By keeping only non-redundant frequency combinations thanks to the parity of the trigonometric functions, u j acquires a matrix form, where N is the number of degrees-of-freedom of the structure:

u j = Ta j , j = 1, 2,...,N (34)
with the notations:

T = [T cos , T sin ], a j = [a j cos , a j sin ] T (35)
and superscript T denotes matrix transposition. The total number of harmonic terms is N h = N cos + N sin . T is a one-row matrix made up of N h trigonometric functions. Submatrix T cos is filled up with N cos elements of cosine terms:

T cos = cos M s ∑ m=1 j 1 m τ m , cos M s ∑ m=1 j
= N sin + 1.
Vector a j is a column of harmonic coefficients belonging to the j th coordinate and is built in a consistent manner with T:

a j cos = a 1 cos , a 2 cos ,... ,a N cos cos , a j sin = a 1 sin , a 2 sin ,... ,a N sin sin (37)
By introducing the matrix Y as follows:

Y = ⎡ ⎢ ⎢ ⎢ ⎣ T 0 ••• 0 0 T ••• 0 . . . . . . . . . . . . 0 0 ••• T ⎤ ⎥ ⎥ ⎥ ⎦ (38)
and by assembling vectors a j as follows:

a = a T 1 , a T 2 ,... ,a T N T ( 39
)
the generalized displacement u, velocity u and acceleration ü can be expressed as:

u = Ya, u = M s ∑ n=1 ω n ∂Y ∂τ m a, ü = M s ∑ m=1 ω m M s ∑ n=1 ω n ∂ 2 Y ∂τ m ∂τ n a (
40) These new expressions (40) are subsequently plugged into equation (30) and a Galerkin procedure is applied over the entire hyper-time domain (32) resulting in a vector of nonlinear equations in a that are to be solved simultaneously:

R = Ha + F a nl (a) (41)
with the notations:

H = Y, M s ∑ m=1 ω m M s ∑ n=1 ω n M ∂ 2 Y ∂τ m ∂τ n + D ∂Y ∂τ m + KY F a nl (a) = Y, F nl (X, Ẋ) (42
) where the following scalar product U , V is used:

U , V = 2π 0 2π 0 ••• 2π 0 U T V dτ 1 dτ 2 ... dτ M s ( 43
)
Frequency/time domain alternance The vanishing of the residual vector R of Eqn. (41) whose unknowns are harmonic coefficients means that the equilibrium is reached:

R = 0 (44)
It has to be pointed out that the nonlinearities involved in our system, such as friction or unilateral contact, can only be accurately captured in the time-domain. Invoking a pioneering work [START_REF] Cameron | An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamics systems[END_REF], solving Eqn. (44) requires a frequency-domain/time-domain alternance at each iteration of a nonlinear solver. In this work, this is achieved by a n-dimensional Fast Fourier Transform in conjunction with a Newton-Raphson-like method.

Autonomous systems Because the system is autonomous, a = 0 is solution of Eqn. (44) but is not of interest. Physically, this solution is equivalent to no structural interaction. With little changes, it is possible to find a non-trivial solution by substituting the M s incommensurate frequencies unknown in advance in Eqn. (41) and consequently prescribing M s harmonic coefficients to be constant in the nonlinear solver: this approach has no influence on the almost periodic response except in the appearance of the phase shift. The new set of equations to be solved is:

H(ω ω ω)a + F a nl (ω ω ω, a) = 0 (45)
where ω ω ω is the vector of the Ms unknown frequencies and a, the vector of harmonic coefficients to be sought whose size is reduced Ms coordinates (respective to the Ms prescribed harmonic coefficients). This is, in essence, equivalent to the usual procedure used in vibration analysis for determining the linear normal modes of a mechanical system.

CONCLUSION

The emphasis of the study has been placed on the understanding of the modal interaction that can occur between a bladed disk and a casing in any aircraft engine. In order to capture this phenomenon and get interesting and realistic results, it is necessary (1) to couple normal and tangential displacements of both structures and (2) to include friction forces in the contact law.

In our tool, the first condition is achieved by curving the blade geometry and the second one, by an explicit time-stepping procedure in conjunction to the method of Lagrange multipliers where the usual Coulomb friction law is considered. Theory and results show good agreement for the prediction of a critical rotational velocity in terms of modal interaction. This is why this numerical tool is very promising for future investigations about the importance of: the initial gap distance between the structure, the curvature of the blades, an external engine order excitation due to aerodynamic forcing. . . Nevertheless, our algorithm, dealing with highly nonlinear terms, is very sensitive to the time-step size. In order to circumvent this sensitivity and due to the (almost) harmonic behavior of the structural set, a frequency-domain method is being developed. It is based on a general implementation of the harmonic balance method coupled to an alternating frequency/time domain procedure.

Figure 1 . 3

 13 Figure 1. 3 NODAL-DIAMETER MODES FOR DIFFERENT KINDS OF AXI-SYMMETRIC STRUCTURES

Figure 3 .

 3 Figure 3. SCHEMATIC OF THE BLADED DISK AND CASING MODELS

Figure 4 . 3 -

 43 Figure 4. 3-NODAL DIAMETER MODES OF THE CASING AND THE BLADED DISK

  contact conditions, referred to as the Kuhn-Tucker optimality conditions in the parlance commonly used in the literature, all x∈ Γ (m) c :

= 7 •

 7 10 7 E bd = 2.1 • 10 11 density (kg•m -3 ) ρ c = 2800 ρ bd = 7800 height (m) h c = 0.01 h bd = 0.01 width (m) l c = 0.1 l bd = 0.1 radius (m) r c = 0.4895 r bd = 0.4795 d.o.f. 160 660 modal damping ξ c = 0.03 ξ bd = 0.005 frequency (Hz) (n d = 3) ω c = 457.3 ω bd = 447.4

2 0d

 2 p (h = 10 -6 s)d c (h = 10 -6 s) d c (h = 10 -7 s) d p (h = 10 -7 s)

Figure 5 .Figure 6 .

 56 Figure 5. PREDICTED (d p ) AND CORRECTED (d c ) DISTANCES BE-TWEEN BLADE 1 AND THE CASING FOR TWO DIFFERENT TIME-STEPS

Figure 7 .Figure 8 .

 78 Figure 7. VIBRATIONS OF THE TWO MODES OF THE CASING FOR Ω = 351.4 RAD/S (Ω Ω c ): MODE 1 ( ), MODE 2 ( )

Figure 9 .Figure 10 .

 910 Figure 9. CRITICAL ROTATIONAL VELOCITIES VERSUS k WITH ON THE LEFT, THE THEORETICAL VALUE AND ON THE RIGHT, THE NU-MERICAL VALUE

Table 1 .

 1 CHARACTERISTICS OF THE MODEL

  τ m ≥ 0. Submatrix T sin is built analogously. No sine term at the zero frequency in T is needed and therefore N cos

	2 m τ m ,...,cos	M s m=1 ∑	j N cos m τ m
			(36)
	where j i m are integers. ∑ M s m=1 j i m τ m represent all combinations which satisfy ∑ M s m=1 j m

ACKNOWLEDGEMENT

This research was carried out for the Société Nationale d' Étude et de Construction de Moteurs d'Aviation. The authors gratefully acknowledge SNECMA for its support.


Mathias Legrand

Bernard Peseux

Christophe Pierre

Investigation of structural rotor-stator modal interaction in turbomachinery

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.

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 [START_REF] Arnoult | Modélisation numérique et approche expérimentale du contact en dynamique: application au contact aubes/carter de turboréacteur[END_REF].

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 jetengine generally called "modal interaction": to the best of our knowledge, one particularly interesting relevant dissertation [START_REF] Schmiechen | Travelling Wave Speed Coincidence[END_REF]) and a paper [START_REF] Legrand | Amélioration de la prédiction de l'interaction rotor/stator dans un moteur d'avion[END_REF] 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 [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF].

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.

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 [START_REF] Bladh | Efficient predictions of the vibratory response of mistuned bladed disk by reduced order modeling[END_REF]. Both respective modeshapes are similar and rotated around the axis of symmetry by =.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  be an angular position of a typical material point of the structure and ˛, an amplitude of displacement of this material point, these two rotating waves are simply expressed as:

˛f .t;  / D ˛0 cos.!t n d .  0 // ˛b.t;  / D ˛0 cos.!t C n d .  0 // (1)
where ! is the pulsation and  0 an undetermined parameter of angular position. For a motionless structure, the propagation velocity of these two waves in a stationary frame is ˙!=n d , depending on the direction of rotation. For a rotating structure like a bladed disk in a turbo-machine, the angular rotational velocity (selecting the counter-clockwise direction as positive) has to be added. The propagation velocity of these two waves in a stationary frame becomes then C !=n d for the forward rotating wave and !=n d for the backward rotating wave.

The two natural frequencies are denoted ! c for the casing and ! bd for the bladed disk. Keeping in mind that ! c is usually greater than ! 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 as illustrated in Fig. 2:

! c D n d ˙!bd (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 the bladed disk if the rotor rotates counter-clockwise -, only one of these equations can be considered as dangerous:

! c D n d ! bd (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.

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 [START_REF] Legrand | Amélioration de la prédiction de l'interaction rotor/stator dans un moteur d'avion[END_REF]: 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.

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/ D 1 3s 2 l e 2 C 2s 3 l e 3 N 2 .s/ D s 2s 2 l e C s 3 l e 2 N 3 .s/ D 3s 2 l e 2 2s 3 l e 3 N 4 .s/ D s 2 l e C s 3 l e 2
where l e is the length of a finite element. The discretized displacement field becomes for an element

u.s/ D N 1 u 1 C N 2 u 1;s C N 3 u 2 C N 4 u 2;s v.s/ D N 1 v 1 C N 2 v 1;s C N 3 v 2 C N 4 v 2;s (4) 3.2 Bladed disk
The blades are discretized by the usual Euler-Bernouilli 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/ D 1 x l e M 2 .x/ D x l e (5)
In a local frame fixed to a rod-bar finite element, the discretized displacement field becomes then :

u.x/ D M 1 u 1 C M 2 u 2 v.x/ D N 1 v 1 C N 2 Â 1 C N 3 v 2 C N 4 Â 2 (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 d -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 d -nodal diameter modeshape of the casing fits its counterpart of the bladed disk as showed in Fig. 4.

CONTACT DYNAMICS

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 decribed in [START_REF] Laursen | Computational contact and impact mechanics[END_REF]. 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 N :

C N D ˚tN W c ! R n t N 0 « (7)
and another one, respective to the admissible friction tractions C T .t N /:

C T .t N / D ˚tT W c ! R 3 n t T n D 0; kt T k Ä t N « ( 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:

Z R uıu dV C Z N N W ı N N " dV D Z t d ıu dS C Z c t N ıg C t T ıu T dS C Z f d ıu dV (9)
where t N and t T 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:

8 < : R u n D u nC1 2 u n C u n 1 h 2 P u n D u nC1 u n 1 2h (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 timestep, smaller than for linear cases, has to be chosen.

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 mentionned by [START_REF] Vola | Consistent time discretization for a dynamical contact problem and complementarity techniques[END_REF]. Because of its simplicity, the Forward Increment Lagrange Method developed by [START_REF] Carpenter | Lagrange constraints for transcient finite element surface contact[END_REF] 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 m , D m and K m . For instance the modal mass matrix of the casing, M c m , is calculated in the following manner:

M c m D P c T M c P c ( 11
)
where the two columns of P c 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 m of the structures without considering any contact. This predicted displacement, denoted with a superscript p , is analytically expressed as:

u nC1;p m D Ä M m h 2 C D m 2h 1  2M m h 2 K m à u n m C  D m 2h M m h 2 à u n 1 m C F ext;n m à (12)
2. projection of the predicted modal displacements onto the physical coordinates system:

u nC1;p D P u nC1;p m ( 13
)
and determination of the gap function vector d p between the two structures. Each gap function where a penetration has been detected is kept in d p , 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:

d nC1 D C N T u nC1;c C d p D 0 ( 14
)
where the superscript c means that the correction of the displacements is being calculated. C N 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 NT containing the normal and tangential constraints is built by considering that, because of high relative velocities between the two structures, only sliding occurs: 8 < :

D C N T Ä M h 2 C D 2h 1 C NT ! 1 d p u nC1 D u nC1;p C Ä M h 2 C D 2h 1 C NT (15)
At the end of this step, the physical displacements verifying the contact constraints u nC1 are projected onto the modal coordinates.

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

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 by a 100µs forcing pulse on the first mode of the casing at the very beginning of the simulation.

Contact algorithm and time step size

It is well-known [START_REF] Wriggers | Computational contact mechanics[END_REF]) 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. Moreoever, 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 D 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.

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 D 3 and c , solution of equation ( 3) such as: c D 301:56 rad/s (16) Fig. 6 shows that there is no dangerous interaction in terms of modal displacement if < 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 > 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 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:

! c C ! bd n d (17)
Instead, equation (3) predicts only one dangerous rotational velocity which is the lower limit of those given by the numerical tool.

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 [START_REF] Pusenjak | Incremental harmonic balance method with multiple time variables for dynamical systems with cubic non-linearities[END_REF]).

  Figure 1. 3 nodal-diameter modes of axi-symmetric structures

Figure 2 .

 2 Figure 2. Potential wave speed coincidence

Figure 3 .

 3 Figure 3. Schematic of the model of the bladed disk and casing

Figure 4

 4 Figure 4. 3-nodal diameter modes of the casing and the bladed disk

  ) c D 2800 bd D 7800 height (m) h c D 0:01 h bd D 0:01 width (m) l c D 0:1 l bd D 0:1 radius (m) r c D 0:4895 r bd D 0:4795 d.o.f. 160 660 modal damping c D 0:03 bd D 0:005 frequency (Hz) (n d D 3) ! c D 457:3 ! bd D 447:4

  h D 10 6 s) d c (h D 10 6 s) d c (h D 10 7 s) d p (h D 10 7 s) Predicted p ) and corrected (d c ) distances of penetration between blade 1 and the casing

Figure 9 .

 9 Figure 9. Critical rotational speeds versus n d with on the left, the theoretical value and on the right, the numerical value

Table 1

 1

	. Characteristics of the model
		casing	bladed disk
	Y. modulus (MPa)	E c D 7 10 7 E bd D 2:1 10 11
	density (kg m 3

ACKNOWLEDGEMENTS

This research was proposed and supported by SNECMA.


Mathias Legrand
email: mathias.legrand@ec-nantes.fr

Christophe Pierre
email: pierre@umich.edu

Stephen P Timoshenko

Bernard Peseux
email: bernard.peseux@ec-nantes.fr

INVESTIGATION OF ROTOR-STATOR INTERACTIONS IN TURBOMACHINERY

In modern turbo machines such as aircraft jet engines, contact between the casing and bladed disks 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 under what circumstances they occur. In this study, we focus on the phenomenon of interaction caused by modal coincidence. Simple two-dimensional models of the casing and bladed disk structures are introduced in order to predict the occurrence of contact and the rotation speeds at which response becomes dangerous. 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 motion. These two pairs of nodal diameter modes are made to interact through direct contact and the equations of motion are solved using an explicit time integration scheme in conjunction with the Lagrange multiplier method on the one hand. On the other hand, the Harmonic Balance Method is used for a better understanding of the solutions found by the time integration method. Both methods generally agree well and exhibit two distinct zones of completely different behaviors of the system.

INTRODUCTION

One of the main objectives of turbo-machines engineering is to increase the efficiency of the machines safely. 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. After this type of contact, several different dynamic phenomena could appear: forced excitation of the modeshapes of a blade resulting in the fragmentation of the tip of this latter or coupling between the vibration modes of the casing and those of the bladed disk, generally called "modal interaction": to the best of our knowledge, one particularly interesting relevant dissertation by Schmiechen (1997) and a paper by [START_REF] Legrand | Amélioration de la prédiction de l'interaction rotor/stator dans un moteur d'avion[END_REF] deal with this travelling wave speed coincidence. In none of the other literature on the topic is there any consideration of the interaction between two flexible structures in contact.

The primary goal of the present paper is to give a better understanding of the modal interaction between the stator and the rotor in a jet-engine. This non linear interaction may lead to devastating 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.

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 modelization of the structures used in the numerical tool is presented. The interaction is investigated by an explicit time integration method, in a third part, and by the Harmonic Balance Method in a fourth section. This method has shown very promising results for strong non linearities [START_REF] Poudou | Hybrid frequency-time domain methods for the analysis of complex structural systems with dry friction damping[END_REF]). Finally, the new results found are explained.

MODAL INTERACTION DEFINITION

The axi-symmetry of a structure, like a bladed disk or a casing, which mathematically leads to circulant mass and stiffness matrices, gives rise to mode pairs with identical natural frequencies, except for zero nodal diameter modes with a simple natural frequency as explained by [START_REF] Bladh | Efficient predictions of the vibratory response of mistuned bladed disk by reduced order modeling[END_REF]. Both respective modeshapes are similar and rotated around the axis of symmetry by =.2n d / where n d is the number of nodal diameters. In that special case, it is possible to combine these two modes into forward and backward travelling waves. Assuming that  represents an angular position and ˛, an amplitude of deflection, these two rotating waves are expressed as:

˛f .t; Â/ D ˛0 cos.!t n d .  0 // ˛b.t; Â/ D ˛0 cos.!t C n d .  0 // ( 1
)
where ! is the pulsation and  0 an undetermined parameter of angular position. For a static structure, the propagation velocity of these two waves in a stationary frame is ˙!=n d , considering the direction of rotation. For a rotating structure like a bladed disk in a turbo-machine, the angular rotational velocity (positive in a counter-clockwise direction) has to be added. The propagation velocity of these two waves in a stationary frame becomes then C !=n d for the forward rotating wave and !=n d for the backward rotating wave. If both structures vibrate at their own natural frequencies, denoted ! c for the casing and ! bd for the bladed disk, four cases of travelling wave speed coincidence can be stated for :

˙!c D n d ˙!bd (2)
Taking into account physical considerations on the direction of the contact force between the two structures -forward on the casing and backward on the bladed disk if the rotor rotates counterclockwise -, only one of these equations can be considered as dangerous:

! c D n d ! bd (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.

MODELIZATION OF THE STRUCTURES 2.1 BLADED DISK

For this study, the bladed disk is modelled as simply as possible.

A given number N b of rigid bodies linked to the disk and between each other by a stiffness constitutes the rotating structure. The structure is depicted in Fig.

m b R ˛i C .k b C 2k ib /˛i k ib .˛i C1 C ˛i 1 / D 0 (4)
where ˛i represents the vibration angle of the i th blade. This equation becomes, when generalized to the entire structure:

M bd R ˛C K bd ˛D 0 (5)
The mass matrix contains the moment of inertia of a rigid beam about one of its ends m b on its diagonal:

M bd D 2

CASING

The casing is directly discretized in the modal coordinates by keeping only its two n d -nodal diameter modes (Fig . 3). Their contri-Fig. 3: Discretization of the casing on its 3-nodal diameter modes butions are denoted u c and u s .During the calculations, structural damping is added in the equations for both structures.

CONTACT EQUATIONS AND TIME INTE-GRATION 3.1 THEORY

The forces of particular interest in this study are the forces of contact acting between the blade tips and the casing. No friction is considered. The model is thus conservative and Hamilton's Principle can be used for the derivation of the equations. The so-called impenetrability condition is enforced by the Lagrange multiplier method. Due to the chosen configuration space, this method has no physical interpretation and may be viewed as a simple mathematical correction of the displacements of the structures in order to avoid any penetration between them. Hamilton's Principle is written as follows:

Z t
T D 1 2 P u T M c P u C 1 2 P ˛T M bd P Ų D 1 2 u T K c u C 1 2 ˛T K bd Ęcont D T g.u; ˛/
where u and ˛are the displacement vectors of the casing and the bladed disk, respectively. and g contain the Lagrange multipliers and the gap functions between the blade tips and the casing. This formulation is discretized in time using the explicit central differences scheme. Denoting the displacement vector of the two structures by x, one has: 8 < : R

x n D x nC1 2 x n C x n 1 h 2 P x n D x n x n 1 h (9)
This scheme is stable for a very small time step h bounded by the Courant Criterion in linear cases. The velocity vector P x n is approximated this way to guarantee the robustness of the following algorithm.

ALGORITHM

Carpenter et al. (1991) explained that an explicit time integration scheme in conjunction with the Lagrange multiplier method dealing with contact problems is generally divided into three steps and is fairly simple:

1. prediction of the displacements without considering any possible contact.

2. determination of the penetrations between the structures.

3. correction of the displacements through the calculation of the Lagrange multipliers in order to eliminate all the penetrations.

The time step size allows for the linearization of the equations with respect to time between the prediction and the correction steps.

The algorithm becomes:

˛ N b 1 ˛ N b 9 > > > > > = > > > > > ; D 2
g i D l c l a C u c cos.ˇ i / C u s sin.ˇ i / (11)
where l a is the length of a blade, l c , the radius of the casing and ˆi the initial angle of the i th blade in a stationary frame.

correction of the displacements Lagrange multipliers have to be calculated such as nulling the gap between the guessed penetrating blades and the casing. Due to the linearization explained above, this correction consists in the resolution of a .2N C 2/ .2N C 2/ linear system where N is the number of blades in contact with the casing during the on-going time step. An example of the linear system is given for a case with two blades in contact, the first and fourth ones:

K C D 2 6 6 6 6 6 6 4
m c 0 0 0 cos ˇ 1 cos ˇ 4 0 m c 0 0 sin ˇ 1 sin ˇ 4 0 0 m bd 0 A 1 0 0 0 0 m bd 0 A 4 cos ˇ 1 sin ˇ 1 A 1 0 0 0 cos ˇ 4 sin ˇ 4 0 A 4 0 0 3 7 7 7 7 7 7 5 U C D 8 < : u nC1 c u c u nC1 s u s ˛nC1 1 ˛ 1 ˛nC1 4 ˛ 4 1 4 9 > > > > > > > > > = > > > > > > > > > ; ; F C D 8 < : 0 0 0 0 g 1 g 4 9 > > > > > > = > > > > > > ;
where

A i D n d . u c sin ˇ i C u s cos ˇ i /. U C is unknown
and is calculated by inverting the system K C U C D F C . The correction for the bladed disk is carried out in the physical coordinates and is projected onto the modal coordinates. Finally, time is incremented and a new prediction is calculated.

RESULTS

In this section, we present some numerical results illustrating the behavior of the mechanical system previously presented due to the contact between the rotor and the stator initiated by a very short forcing pulse of fifty time steps (50 s) on the first mode u c of the casing at the very beginning of the simulation. For sake of simplicity, the rotor is constituted of eight blades and both structures are discretized on their two three-nodal diameter modes. It has to be noted that the rotor of the chosen model is perfectly balanced and no eccentricity between the two structures is introduced.

VALIDATION

The validation of the algorithm has been carried out in two different ways. First, a convergence study of the time step has shown that h D 10 6 s is sufficient to get significant results while minimizing computational efforts. Next, a focus on the distances between the blade tips and the casing depicted in Fig. 5 confirms that the contact algorithm is successful for the chosen time step: the gap distance never becomes negative and the algorithm remains stable for all times. It is interesting to note the accuracy of the contact-algorithm based on a four-degrees of freedom model.

RESPONSE OF THE STRUCTURES

Browsing the space of initial conditions, two main types of responses appear. The first one is characterized by several impacts between the two structures, followed by a decrease of the amplitudes of vibration to zero due to structural damping. This response is not of interest. In contrast, an exchange of energy between the stator and the rotor occurs in the second type of responses. This interesting case features, after the transient response, a steady state where several blades remain in permanent contact with the casing without vibrating, while the other blades do not impact the casing anymore. Solid and dashed bold black lines in Fig. 5 show the two blades in permanent contact for a numerical simulation of an 8-bladed disk interacting with the casing. The propagation of the deformation wave of the casing has the same angular velocity as the bladed disk, and the latter reaches a permanent static-only deformation state as observed in Fig. 6. This new phenomenon is referred to as locked motion and is studied in details in what follows.

CONDITIONS OF INTERACTION

It appears that several parameters (rotational velocity of the rotor, initial impact amplitude, structural damping...) and the combi- ciently large to make the interaction phenomenon appear. Below ! c =n d whatever the initial conditions, the locked motion is impossible. A deeper study taking into account the structural damping of the casing and the bladed disk is rendered difficult because of the addition of two new parameters. Then, a clear conclusion on the time integration results is early and a new approach is used and explained in section four. Although the few previous findings have shown the importance of the natural frequency of the bladed disk ! bd (see Eq. 3), Fig. 6 demonstrates that the role of this structure in the phenomenon vanishes because of its permanent static-only deformation state behavior during the interaction.

STABILITY OF THE "LOCKED-MOTION"

After the locked motion is reached, a new forcing pulse of fifty time steps is applied on the first mode u c of the casing. For a very high amplitude of this pulse, the motion jumps from a locked motion 1 to a locked motion 2 instantaneously. Even if the blades in permanent contact are different, their number remains the same. The motion is then stable in amplitude but is able to take many states as illustrated in Fig. 9 where the red arrows indicate the blades in permanent contact before and after the forcing pulse.

before the pulse after the pulse In order to avoid the high sensitivity of the results to initial conditions, and particularly to the amplitude of the initial forcing pulse, the method of Harmonic Balance (HB) is used. This frequencydomain method is capable of predicting efficiently the steady state periodic behavior of nonlinear mechanical systems without considering the initial conditions [START_REF] Pierre | Multi harmonic analysis of dry friction damped systems using an incremental harmonic balance method[END_REF]). The HB method also allows for parametric studies that would be carry prohibitive cost using direct time integration, and results in a better physical understanding of the solutions.

GOVERNING EQUATIONS

The time integration model is kept in this section. The two blades in permanent contact with the casing are denoted by i 1 and i 2 . The resulting equations are: for the casing 8 <

:

m c R u c C k c u c C 2 X j D1 j cos ˇij D 0 m c R u s C k c u s C 2 X j D1 j sin ˇij D 0 (12)
for the two blades in contact

8j D 1; 2 .k a C 2k b /˛i j k b .˛i j 1 C ˛ij C1 / C j n d . u c sin ˇij C u s cos ˇij / D 0 (13)
for the remaining blades

8 < : .k a C 2k b /˛1 k b .˛N C ˛2/ D 0 .k a C 2k b /˛i k b .˛i 1 C ˛iC1 / D 0 8 i 2 OE2; N 1 fi 1 ; i 2 g .k a C 2k b /˛N k b .˛N 1 C ˛1/ D 0 (14) for blades i 1 and i 2 8j D 1; 2 u c cos ˇij C u s sin ˇij D l a l c (15)
In order to decrease the size of the problem, the bladed disk displacements are written in the modal coordinates. This formulation results in a non linear differential system of six equations and six unknowns u c , u s , 1 , 2 , ˛c, ˛s that cannot be solved analytically.

SOLUTION METHOD

The HB method is used to find a solution to this non linear differential system for several reasons: this approach is a classical technique for studying non linear system oscillations. in this method, the response being sought and the non linear restoring forces are expanded in two Fourier series: the time integration has shown almost pure harmonic or static responses HBM is not sensitive to the initial conditions which represent a critical issue in the time integration method. HBM transforms a non linear differential system of equations in a non linear algebraic system of equations that can be easier to solve.

Observing time integration results, a simple harmonic is kept for the casing:

( u c D u c c cos.n d t/ C u s c sin.n d t/ u s D u c s cos.n d t/ C u s s sin.n d t/ (16)
The two contact forces and the two bladed disk modal contributions are assumed to be constant. A projection of Eqs 12, 13, 14 and 15 onto OE1; cos.n d t/; sin.n d t/ basis followed by an integration over the domain t D OE1; 2 =.n d / (Galerkin method) leads to a non linear algebraic system of eight equations and eight unknowns. The vector solution is obtained iteratively by a non linear algorithm such as the quasi-Newton method of Matlab in our case.

RESULTS

The previous equations have been solved with as a parameter. All values are -dependent except the vibrations of the casing which are bounded by the initial clearance between the blade tips and the casing in itself. D ! c =n d seems to be a particular rotational velocity for this system: it is confirmed by modifying the mechanical characteristics of the casing and thus its natural frequency ! c . The sign of the Lagrange multipliers changes exactly for this value. The solution given by the HB method in Fig. 10 is physically unacceptable for < ! c =n d though numerically correct. Further simple calculations on the contact forces would show that the casing is pushing the bladed disk for < ! c =n d which is impossible. This main result is consistent with the pre-supposed interaction zone of the time integration results depicted in Fig. 8: the locked-motion is not reachable for Ä ! c =n d whatever the initial forcing amplitudes are. The introduction of an engine order

CONCLUSION

In this paper, a numerical tool has been developed in order to study the modal interaction between a rotor and a stator in a turbomachine. The model used is based on a exact mathematical formulation of the contact in conjunction with a modal description of the motion. A zone of possible interaction referred to as locked motion in this study, because of the characteristics of its behavior, has been found by both HB and time integration methods. Contrary to previous findings, the bladed disk natural frequency plays no role in the occurrence of the locked motion phenomenon. Moreover, the initial gap between the blade tips and the casing limits the amplitudes of vibration of the casing and forbids any amplitude divergence. In future work, it would be of interest to include a small eccentricity between the structures in order to apply a rotating excitation due to contact on the bladed disk as well as the casing.

1

 1

Fig. 1 :

 1 Fig. 1: Schematic of the bladed disk i th blade (or rigid body) is:

Noting k 1 D

 1 k b C 2k ib and k 2 D k ib the circulant stiffness matrix becomes: : 0 k 2 k 1 k 2 k 2 0 : : : 0 k 2 k 1 eigenshapes of the bladed disk are calculated from this model. Their two n d -nodal diameter modal contributions are denoted ˛c and ˛s.

Fig. 2 :

 2 Fig. 2: 3-nodal diameter modes of the bladed disk

predictionFig. 4 :

 4 Fig. 4: Schematic of the gap function g i

Fig. 5 :Fig. 6 :Fig. 7 :(Fig. 8 :

 5678 Fig. 5: Distances between the casing and the blade tips

Fig. 9 :

 9 Fig. 9: States of the two structures before and after the second forcing pulse

Fig. 10 :

 10 Fig. 10: Amplitudes of the contact forces versus

  / D ı˛.t 2 / D ıu.t 1 / D ıu.t 2 / D 0. The different energies take the forms:

	2
	ıT ıU C ıE cont dt D 0 8.ı˛; ıu; ı /	(8)
	t 1
	satisfying ı˛.t 1

P. Schmiechen. Travelling Wave Speed Coincidence. PhD thesis, Imperial College of Science, Technology and Medecine -University of London, 1997.