Whole Engine Interaction in a Bladed Rotor-to-Stator Contact

Reducing the clearances between rotating and ﬁxed parts is an important factor in increasing the performances of turbomachines. The physical counterpart however is an evolution in possible rotor-stator contacts capable of causing unstable dynamic behavior. A proper prediction of the rotor-stator contact occurrences and associated induced phenomena, has therefore become of a great interest for aero-engine mechanical engineers. Most numerical simulations involving rotor-stator contact can be divided into two types of physical behavior. The ﬁrst focuses on contact induced blade/casing interactions, in only taking into account the blades and casing ﬂexibility. The second type of behavior takes into account the shaft dynamic while neglecting blade ﬂexibility. Future designs of aircraft engines will however raise the need to combine these two types of models. Since, the structural components are more ﬂexible, the dynamic coupling between engine modules is increased. This paper proposes a study based on a structure representative of the whole aircraft engine, including the contacts that may arise between the fan-blade tips and fan casing. We have introduced a fully-coupled phenomenological model with ﬂexible blades, shaft and casing. Furthermore, this model includes an elastic link between shaft and casing to simulate the fan frame behavior. We begin by explaining the linear results, which highlight the dynamic couplings between these various model components. During a second step, this paper presents the nonlinear results obtained by introducing a con-∗


INTRODUCTION
Turbomachine performance depends to a great extent on the clearance between rotating and fixed parts.Unfortunately, while a closer clearance improves efficiency, it also increases the probability of generating blade tip-casing contacts capable of leading to unstable dynamic behavior.Consequently, in order to reduce clearances, aero-engine engineers must ensure that the possible rotor-stator contacts do not trigger damaging phenomena.
The study of turbomachine dynamics is commonly divided into two distinct fields.On the one hand, bladed assembly vibration takes into account the blades and disk flexibility, while shaft vibration is neglected.On the other hand, rotor dynamics focus on lateral and torsional vibrations of the shaft in neglecting the blades and disk flexibility.Most of the literature reviewed on blade tip-casing contact is split between this two areas of study [1]: -In the bladed assembly vibration category, some numerical and experimental investigations focus on the travelling wave speed instability [2][3][4].This phenomenon (also called modal interaction) is related to the specific modal properties of axisymmetrical structures and may occur when the mode shape geometry of the bladed disk and the casing match and their travelling wave speeds coincide.

1
-In rotor dynamics, the investigations focus on the behavior of the flexible shaft.Most studies are found to use simplified models whereby the blade tip-casing contacts are introduced using nonlinear springs.Such models typically consider contacts between two cylindrical surfaces [5][6][7][8].
Besides tip clearance, another governing factor of jet engine design is weight.Reducing the mass generally leads to an increased flexibility and/or a greater ability to transmit vibrations.Future designs of aircraft engines will thus expose the need to extend the dynamics included in the prediction models so as to consider the couplings between engine components.In this context, a recent study [9] conducted an analysis on an industrial Whole Engine Model (WEM) including casing deformations and assuming rigid blades.
The present work focuses on the contacts that may occur between the fan blade tips and the fan casing of a turbofan jet engine.In the fan module, the blades are particularly slender and the shaft precession orbit may be significant.The blades and shaft dynamics must therefore be considered together.Sinha [10] and Lesaffre [11] developed flexible, bladed shaft models for this type of blade tip-casing contact investigation.The present work follows their approaches.The model developed in this paper has been based on the described model [11] into which we have introduced rigid-body displacements of the casing and a viscoelastic link between shaft and casing.
The first two sections of this paper will describe the phenomenological model and contact formulation.Linear results will then be presented to highlight the natural couplings and effect of permanent contact on the various modal properties.Lastly, transient analyses will investigate unstable scenarios involving two different modes, one of which enables an atypical contact configuration.These analyses will consider light intermittent contacts with rubbing.

MODEL DESCRIPTION
This section will describe the assumptions and the equations sequences used to develop the phenomenological model considered in this study.The global model architecture is diagrammed in Fig. 1; it consists of a flexible bladed rotor, a suspended flexible casing and a viscoelastic link between the shaft and the casing.To avoid time-dependent terms resulting from the periodicity of the rotating structure, the entire model has been developed in the rotating frame.

Flexible bladed rotor
The bladed rotor model used in the present paper has been derived from the model built by Sinha [10].It was first developed in the rotating frame through use of an energy-based method by Lesaffre [11] and then extended to a dual shaft by Gruin [12].The shaft is modeled by a Euler-Bernoulli beam suspended by isotropic viscoelastic bearings and connected to two rigid disks modeled by a point mass and rotational inertia.A set of flexible blades also modeled by Euler-Bernoulli beams, are clamped to the representative fan disk (see Fig. 1).
The shaft beam has a hollow circular cross-section and includes two displacements defined in the rotating frame: the two orthogonal translations x s and y s in the crosssection plane.The blade beams feature a constant rectangular cross-section and can deflect along their more flexible direction: x b .The displacements are formulated with Rayleigh-Ritz approximation functions that respects the geometric boundary conditions.
The shaft deformations are given by: where z is the position along the shaft axis and L s its length.X 0 (t) and Y 0 (t) are the rigid-body translations of the shaft.W m (z) is the shape function chosen to approximate its motion and m tot the number of modes under consideration.Since the shaft is supported by bearings, its shape function has no geometric boundary conditions to verify.The chosen shape function has been used in [12], it allows motion at both ends of the shaft which is not possible with the shape function used in [10] and [11].
The blade deformation is given by: where l b is the blade length and s the position along the blade (s = 0 indicates the blade foot and s = l b is the blade tip).n tot is the number of modes considered to represent the j th blade deformation.Y n (s) is the shape function of the blades; it respects the geometric boundary conditions: Y n (0) = 0 and Y ′ n (0) = 0 (where prime denotes differentiation versus space coordinates).
The bladed shaft energies are then developed as described in [13] and the system of equations is obtained by means of Lagrange equations.The bladed shaft model takes into account the gyroscopic effect, spin softening and stress stiffening of blades.

Suspended flexible casing
The system of equations governing the casing behavior will first be developed in the inertial frame and then transposed into the rotating frame.The casing model consists of an elastic ring with radial and tangential deformations, as in [11].It also contains four rigid-body degrees of freedom (dof) contributed by the casing suspension: two orthogonal translations and rotations in the cross-section plane of the shaft.The rigid-body dofs and deformation dofs of the casing are dissociated.This disconnection is related to both the deformation definition (Eqn.(4)) and kinetic energy formulation.Moreover, it stems from the orthogonality property of trigonometric functions as well as a linearization of the inertia effects.
We begin by developing the system of equations related to the four rigid-body dofs of the casing.In the inertial frame, this system comprises three 4x4 diagonal matrices: mass matrix M c , damping matrix D c and stiffness matrix K c .This system is then transposed in the rotating frame under the hypothesis of a constant rotation speed [14].The rotation-related matrices are defined as follows: P c is the transition matrix used to project the inertial frame dofs into the rotating frame.G c is called the gyroscopic matrix, N c the spin softening matrix, and C c the damping circulation matrix.
The equations of motion related to rigid-body displacements can now be established.Next, we focus on the casing deformations.The fan casing is modeled by an axisymmetric elastic ring with radial and tangential deformations.These are expressed in considering the nodal diameter modes of the axisymmetric structure.Each mode is introduced by a pair of general coordinates (A n , B n ).The tangential displacement expression is: Where θ is the angular position on the casing.n d is the number of nodal diameters considered: n d ≥ 2 because the n d = 0 mode shape is incompatible with the inextensible property and n d = 1 is provided by the rigid-body dofs.
The inextensible property of the elastic ring implies that the radial displacement can be expressed from the tangential displacement: The equations of kinetic energy, internal deformation energy and dissipation function are defined in [15] : where: E c the Young modulus of the casing, I z the second moment of area of the casing cross-section along the z axis, R c the radius, S c the cross-sectional area, ρ c the density, η c the viscous damping coefficient.
Substituting Eqn.(4) and Eqn.(5) into Eqn.( 6) and applying Lagrange equation leads to the system of equations related to the deformation dofs of the casing in the inertial frame.The mass, damping and stiffness matrices are all block diagonal matrices, where the blocks related to the n d nodal diameter are given by: The transition matrix P n between the inertial and rotating frames depends on the nodal diameters under consideration.The block relative to the n d nodal diameter is given by: Assuming a constant rotation speed, the transition from inertial frame to rotating frame creates three rotation-related matrices all of which are block diagonal with blocks related to the n d nodal diameter, as given by:

Elastic link between shaft and casing
Once the bladed shaft and casing have been modeled, we can now insert a link to represent the first bearings and fan-frame behavior.
To model this coupling, we introduce an isotropic viscoelastic link between the rigid-body dofs of the casing and the displacements of the appropriate shaft cross-section (see Fig. 1).This link offers the following translation and rotation stiffness and damping properties: (k, k φ ) and (c, c φ ).

CONTACT FORMULATION
A simple contact model is considered herein.A gap function has been developed to evaluate the blade-to-casing distance.Should a contact be detected, then opposing contact and friction forces will be applied to the rotating and fixed parts.The normal contact forces are proportional to the penetration: the radial stiffness k r is introduced.Tangential contact forces are obtained using Coulomb's law with friction coefficient µ.
In this section, we will begin by describing the gap function created to detect contact occurrences.Afterwards, the contact force vector entered into the system of equations will be presented.

Gap function
Contact occurrences are tracked by measuring the minimum distance between each blade tip point and the casing inner surface.This gap function takes into account: the disk and casing translations and rotations, blade flexure, and casing radial deformation.To simplify the formulation, we propose the following assumptions: -Disk and casing rotations are assumed to be sufficiently small so as to linearize their trigonometric functions.
-The influence of the casing tangential deformation is neglected.
-The influence of global blade tip displacement on the angular position is also neglected: the normal casing deformation introduced into the gap function is measured at the initial blade tip position.
To obtain gap function g j , we first determine the position of the j th blade tip in a frame attached to the casing: where: The transformation matrices P are given in the Appendix.
Once the blade tip position is known, the gap function can be obtained as sketched in Fig. 2: where g 0 is the initial tip clearance, and u(α j , t) the normal casing deformation.

Contact forces
If g j is negative, then contact is detected and the normal and tangential contact reactions are introduced: The relative speed between contact surfaces used to define the friction force direction is assumed to be constant and opposite to the sliding speed direction.These reactions are then projected onto the general coordinates in order to obtain the contact force vector included in the system of equations:

RESULTS AND DISCUSSION
The results obtained are presented in three sections.First, a modal analysis of the model without contacts shows

Natural couplings
Performing a modal analysis without considering contacts permits to detect the interactions that naturally exists within the model components.Because of its formulation, the model presents indeed blade/shaft and shaft/casing couplings that appear without necessitating contacts.
The Campbell diagram in Fig. 3 displays the evolution of the first eigenfrequencies vs. shaft rotation speed.In the rotating frame, downward sloping lines depict the forward modes (precession motion in the same direction as shaft's own rotation), while upward sloping lines indicate backward modes.
Because the model includes 10 flexible blades, the first bending deformation of the blade corresponds to 10 modes of the bladed disk arranged in Nodal Diameters (ND): from 0ND to 5ND.These bladed disk modes are indicated by the letter A in Fig. 3.A closer look at this group of modes is presented in Fig. 3 lower diagram and Fig. 4 displays the A(0ND), A(1ND) and A(2ND) mode shapes.A(1ND) mode shape displays a coupling between shaft and blades deformation: a portion of the A(1ND) deformation energy belongs indeed in the shaft and is transmitted to the casing through the elastic link.Meanwhile, the other bladed disk modes imply solely blades deformation.
A locus-veering (denoted L1 in Fig. 3) is also observed between the C and A(1ND) forward modes: the eigenvalues are nearly about to cross, but instead they veer away from one another and exchange mode shapes [16].The same qualitative observations were made by Gruin [12] and Lesaffre [11].Blade/shaft coupling is possible solely with the 1ND bladed-disk modes; this finding pertains to shaft cross-section kinematics.The shaft is modeled by a Euler-Bernoulli beam and can only deflect into its cross-section plane (i.e.torsional and longitudinal deformations are neglected).The projection of bladed disk mode kinematics at the disk center thus reveals that only a 1ND mode is compatible with the shaft cross-section kinematics.
A global observation of the modes also indicates couplings between blades, shaft and casing.The upper three mode pairs in Fig. 3, as indicated by the letters C, D and E, are involved in two locus-veering L2 and L3.Furthermore, even though the main deformation observed in their mode shapes is shaft bending, these modes also involve 1ND blade deformation and casing displacement.
The blades/shaft coupling described above is provided by gyroscopic forces while the shaft/casing coupling is due to the stiffness introduced to tie these two component displacements together.
The only modes with no coupling are B modes: the 2ND casing modes that deform the casing into an elliptical ring.This absence of coupling is consistent with the system of equation built for the casing.Moreover, as explained in regard to the blade/shaft coupling, the only modes capable of coupling with shaft bending are 1ND modes.This compatibility criterion also applies to the casing/shaft coupling: the only casing modes that can couple with shaft bending modes are the rigid-body casing modes, i.e. its 1ND modes.

A single blade in permanent contact
To study the global effect of blade/casing contact on modal properties, we have carried out an eigenvalue analysis in considering that one of the blades remains in permanent contact.To simplify the analysis, the results plotted in Fig. 5 exclude viscous damping and friction.The upper graph illustrates eigenfrequency evolution vs. rotation speed, while the lower diagram tracks the evolution of eigenvalues real parts yielding stability information: the eigensolution is unstable as long as it displays a positive real part [14].
The blade contact causes two kinds of instabilities.Divergence is a form of instability at zero frequency; its appearance is indicated by the letters D in Fig. 5. Secondly, flutter instability appears at eigenvalue crossover; these are indicated by the letter F in the same figure.These phenomena are due to the asymmetric elastic and inertial couplings related to the blade/casing contact and are described in [17][18][19].

Intermittent contact simulations
The preceding results are consistent with the literature and underscore our model's steady behavior.These results however have been obtained under the restrictive assumption of permanent contact, whereas other authors [2][3][4] showed that blade/casing contact related phenomena often involve intermittent contacts.In order to consider such intermittent contacts, transient simulations have therefore been run.The transient results presented herein were derived by applying the explicit central differences algorithm, as in [3,4].Unlike the previous results, friction is included whenever contact occurs.Viscous damping is also introduced in the casing, the suspensions, the shaft/casing link and between each blade.Contact is triggered by considering imbalance at the bladed disk location (see Fig. 1).The imbalance force is applied over a 1 sec ramp and then maintained throughout the simulation.In the rotating frame, imbalance creates a static force given by: Where m b is the imbalance mass, e b the eccentricity and γ b its angular position in the rotating frame.f b is then projected onto the general coordinates to be included in the system of equations.
In rotor-stator contacts, friction usually induces backward directed phenomena on the rotating parts.When flexible blades are introduced, friction forces trigger a backward rotating wave deformation on the bladed disk [3].When shaft dynamics are considered, contacts generally occur on the shaft deflection side and induce a backward precession [20].
The proposed model allows for atypical contact configuration and we have chosen to focus on two specific modes: -Mode C concentrates most of its deformation energy in the shaft and especially in the fan region.Its mode shape, drawn in Fig. 6 (a), presents a regular contact configuration: the minimum gap function values are found on the side of shaft deflection.
-Mode D is also a shaft-bending mode.The main deformations are localized in the turbine area.Its mode shape,   (see Fig. 3).
The transient simulation results for mode C are displayed in Fig. 7.The upper diagram plots the evolution of normal contact forces at the blade tips.The lower table describes the deformed shape during the transient simulation: This presentation starts by describing the various contact configurations obtained during the simulation.Next, the disk and casing center-of-gravity trajectories are characterized.These trajectories are described in the rotating frame for the purpose of detecting the effect of contact forces on displacements.During most of the simulations, these trajectories are in fact mainly driven by the imbalance forces.Hence, if they were to be observed in inertial frame, the forward imbalance-related precessions would counter the impact of contacts on the disk and casing displacements.Moreover, this table presents the nodal diameter analysis of displacements observed on the blades and casing.The colored stripes illustrate the evolution in the Discrete Fourier Transform (DFT) amplitude of the respective nodal diameters.Each one of these stripe results is then normalized, and the color scale is presented in the lower part of the figure.For the bladed disk, two nodal diameter analyses have been conducted.The first analysis applies a Fast Fourier Transform (FFT) algorithm to the radial displacement of all 10 blade tips.In adopting blade modeling assumptions, the radial displacement is exclusively due to the disk translations and rotations.The only possible deformation thus contains just one nodal diameter (as explained above, this is the only deformation compatible with the shaft cross-section kinematics).A second analysis proceeds by applying an FFT algorithm to the blade deflection measurements.In this case, the possible nodal deformation extend from 0ND (i.e. the 10 blades deflect in the same direction as displayed in Fig. 4) to 5ND (each consecutive blade deflects in a different direction).Only the 1ND and 2ND stripes are displayed because the other nodal diameter contributions are low.For the casing, the nodal diameter analysis is carried out on the radial displacements measured at various locations distributed over its circumference.These radial displacements are correlated with the casing's own deformation, as well as with its rigidbody translations and rotations.Only 1ND and 2ND deformations are observed: the 1ND deformation is solely due to the rigid-body displacements, and the 2ND deformation is related to the casing deformation (with the 2ND shape being the only deformation included in the model).D mode results are similarly presented in Fig. 10.
The two transient results presented in this paper differ, yet both follow the same scenario.First, the imbalance forces causes a quasi-static deformation (in the rotating frame) until the deformations are sufficient to absorb the initial clearance and contact occurs on two consecutive blades.Contact forces excite the model components, and the system vibrates with increasing amplitude, thus leading to intermittent contacts.This scenario accelerates amplitude increase until the simulation is stopped because the deformations have become incompatible with the modeling assumptions.
The permanent contact phase on the C mode simulation occurs from 0.98 to 1.07 sec.As expected by the mode shape observation, contact occurs on the side of shaft deflection, which mean that contact forces are creating a backward precession of the disk and a forward precession of the casing.The deformed shape at 1 sec is drawn in Fig. 8 (a).The bladed disk reveals a 1ND deformation fixed in the rotating frame.This synchronous 1ND shape results from: 1ND radial displacements due to disk motion, and 1ND bending deformations of the blades.On the casing, a synchronous 1ND co-rotating wave can be observed, whereas the 2ND deformation is nearly absent.
The permanent contact phase on the D mode simulation occurs from 0.94 to 1.16 sec.Unlike the C mode, these contacts occur on the opposite side of shaft deflection.Hence, contact forces cause a forward precession of the disk and a backward precession of the casing.The qualitative observations on the bladed disk and casing are identical to those found in C mode observations.The deformation shape at 1.16 sec is drawn in Fig. 9(a).
The deformed shapes obtained during the permanent contact phase are similar to the respective forced responses of the system subjected to an imbalance.These contacts are indeed soft and the nonlinear forces are insufficient to overcome the imbalance effect.
For the C mode simulation, the intermittent contacts phase begins at 1.07 sec.The backward circular precession of the disk center-of-gravity increases while the forward precession of the casing becomes irregular.On the bladed disk, the nodal deformations that had been synchronous during the permanent contact phase (i.e.fixed in the rotating frame) now behave like rotating waves.The radial displacements that were mainly driven by the shaft deformed-shape due to imbalance now become influenced by the rising backward precession due to contacts.The nodal deformation due to blade flexure behave as counter-rotating waves with an increasing amplitude.On the casing, the 1ND contribution stemming from its rigid-body displacement for the most part remain defined by the imbalance-related precession in the inertial frame.The 2ND co-rotating wave on the casing increases substantially and becomes the main contributor to radial casing displacements.This casing deformation caused by contact influences the blades, and the 2ND counter-rotating wave on the bladed disk increases to match the casing shape.The deformed shape at 1.15 sec is drawn in Fig. 8(b).
For the D mode simulation, the intermittent contact phase begins at 1.16 sec.The forward precession of the shaft becomes irregular and reverses to a backward motion.Meanwhile, the backward circular precession of the casing rises.On the bladed disk and casing, the evolution is qualitatively similar to the C mode results; however, the 2ND deformations do not overcome the 1ND deformation as it is the case in the C mode simulation.The D mode deformation shape at 1.16 sec is drawn in Fig. 9 (b).
In returning to the global observation, we can identify the consecutive phenomena that led to an unstable system for both analyses.The C and D modes are both whole engine 1ND modes, as reflected by: shaft bending creating a 1ND-like displacement of the disk; shaft/blade coupling due to the gyroscopic effect causing a 1ND deformation of the blades; and the elastic link between shaft and casing leading to a 1ND-like displacement of the casing.The imbalance forces thus naturally create a 1ND global interaction.Moreover, the deformations due to imbalance absorb the initial tip clearance, which leads to blade/casing contacts.Contact forces increase the existing 1ND global interaction while at the same time excite the 2ND casing mode, which adds a 2ND interaction between the casing and the bladed disk.
The disk and casing center-of-gravity trajectories illustrate the complexity of system motion.This complex behavior is correlated not only with the 1ND and 2ND combination, but also with the linked displacements of the shaft and casing.The elastic link introduced between the casing and shaft does actually prevent the respective trajectories from extending in opposite directions.

CONCLUSION
This study has presented a phenomenological model developed to investigate the set of phenomena induced by blade-casing contacts occurring in the fan stage of an aircraft engine.It has involved a flexible bladed shaft and flexible casing.In addition to the blade-casing coupling due to contacts, this model has included the natural blade-shaft and shaft-casing couplings.Consequently, the considered model, which contains fewer than 40 dofs, can still yield complex mode shapes encompassing all the components.
Linear results first displayed the natural couplings and the model's consistent behavior under permanent contact.Then, nonlinear transient simulations were run for the purpose of investigating unstable scenarios related to bladecasing contacts.We chose to focus on two specific modes, one of which presented an inverted contact configuration.Results showed the influence of each component dynamic on the unstable scenarios.The main conclusion drawn from both analyses was that a combination was generated with one nodal diameter (1ND) interaction involving the whole model and a 2ND interaction localized on the blades and casing.

FIGURE 2 .
FIGURE 2. SKETCH OF THE GAP DEFINITION r d : the fan-disk radius, z c : the casing center of gravity position along the z axis, x d , y d : disk translations, x c , y c : rigid-body casing translations, x b (l b , t): blade deflection at its tip, P β and P αj : blade positioning matrices, P φx d and P φy d : disk rotation matrices, P φx c and P φy c : rigid-body casing rotation matrices.

FIGURE 5 .
FIGURE 5. MODAL RESULTS CONSIDERING ONE BLADE IN PERMANENT CONTACT (FRICTION AND DAMPING EXCLUDED)