Aircraft Engine Structural Rotor-Stator Modal Interaction

,


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 [1].
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 [2] and a paper [3] 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-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 [4].
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 [5]: with: In Eqn.(2), the number of blocks on the diagonal is equal to the number of sectors, N, and 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: where α 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): The skew-symmetric structure of Y k gives rise to mode pairs with identical natural frequencies [6].By partitioning the eigenvector matrix such as x k = [x c,k , x s,k ] T , both respective modeshapes take the form: 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: In Eqn.(6), 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 physical coordinates by combining Eqns.( 5) and ( 6) : It is then possible to combine these two modes into forward ( f ) and backward ( b ) travelling waves at pulsation ω:

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 [7].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: 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: 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 [3]: 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: 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:

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: In a local frame fixed to a beam-bar finite element, the discretized displacement field becomes then: 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).

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 [8].
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 Γ According to these notations, the gap function between the two structures can be stated: where g 0 (x) represents the initial positive gap between the two structures and n, the outward normal to Γ where t N stands for the contact pressure, assumed positive, acting on Γ (s) c .To these unilateral contact conditions, we add the well-4 known Coulomb friction law: 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 ): An alternative equivalent form to Eqs. ( 18) is: 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: 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:

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) [9].Because of its simplicity, the Forward Increment Lagrange Method [10] 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: 2. projection of the predicted modal displacements onto the physical coordinates: 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: 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: 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 [11] 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 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 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: 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: 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.

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 [12], 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: First, true time t is adimensionalized with respect to the M s frequencies ω m of the response [13]: so as to construct an hyper-time domain: 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: 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: with the notations: 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: Vector a j is a column of harmonic coefficients belonging to the j th coordinate and is built in a consistent manner with T: By introducing the matrix Y as follows: and by assembling vectors a j as follows: the generalized displacement u, velocity u and acceleration ü can be expressed as: 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: with the notations: ) where the following scalar product U , V is used: 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: 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 [14], 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: 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 3 .
Figure 3. SCHEMATIC OF THE BLADED DISK AND CASING MODELS

Figure 4 . 3 -
Figure 4. 3-NODAL DIAMETER MODES OF THE CASING AND THE BLADED DISK (m) c , its closest counterpart ȳ on the slave surface Γ (s) c : ȳ = arg min y∈Γ (s) c x − y (16) contact conditions, referred to as the Kuhn-Tucker optimality conditions in the parlance commonly used in the literature, all x∈ Γ (m) c :

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

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

Table 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