Investigation of Rotor-Casing Interactions in the Centrifugal Compressor of a Helicopter Engine

The phenomenology of rotor-casing setups experiencing contact interactions is still poorly understood, particularly when complex geometries such as centrifugal compressors are involved. Although interaction phenomena have been witnessed and recorded during industrial experiments, the physical understanding of what occurs during these interactions is limited. The usual design approach is to consider possible modal interaction points in a linear framework and move these outside of normal operating conditions by means of minor geometric changes. Based on this linear approach, no information on the severity of these interactions is available to the designer. Besides, a possible interaction point appearing in the linear framework may not produce any harmful interactions, thereby increasing design restrictions. Based on an in-house numerical strategy previously presented, contact interactions for a flexible centrifugal compressor from a helicopter engine and rigid casing setup are investigated. By imposing a small deformation on the casing geometry, blade/casing contact is initiated and subsequent interactions feature complex phenomena that are analyzed. In comparison to previous interaction simulations involving axial compressors, a higher degree of complexity of the numerical simulations stems from a strong curvature of the blade and very significant bladedisk coupling. This coupling presents itself towards the trailing edge where compressor mode shapes indicate a significant component normal to the casing surface. Accordingly, these modes may lead to large amplitude contact forces. Time simulation results are confronted with experimental observations, and the consistency of the behavior of the numerical model with respect to industrial observations is underlined. A frequency domain post-processing of the results reveals specific engine order interactions and frequency spectra are plotted in order to interpret the phenomenon of interest. Such methodology will enable designers to more efficiently discriminate potential critical interaction speeds as compared to the classical linear frequency approach. INTRODUCTION Rotor/stator interactions [1] involve a large variety of complex phenomena. In modern turbo-machinery clearances between the rotor and casing are usually designed at normal operating conditions, allowing for optimal performance without significant interactions between the two parts. Nevertheless, in certain cases an engine may experience operating conditions for which clearances can be overcome and interactions can become non-negligible. As a result, contacts, such as blade/casing contacts, may occur. Subsequent interaction phenomena — which are a specific type of rotor/stator interactions — that were also observed for axial compressors both experimentally [2, 3] and numerically [4] are still under investigation. Several references may be found regarding the study of modal dynamics of centrifugal compressor geometries [5–8]. More recently, a detailed numerical model of the geometry of interest in this paper was presented in [9]. Along with other interaction phenomena such as impeller-vaned diffuser interaction [10], the understanding of interaction phenomena arising from structural contacts is a growing field of investigation. Consequences of these phenomena may result in strong energy vibration between rotor and stator which may cause damage to the structure, as seen in [11]. The proposed study focuses on blade/casing contacts within the centrifugal compressor of a helicopter engine. One specificity of centrifugal compressors lies in a complicated blade-tip geometry. Due to the high curvature of almost 90 between inlet and outlet as well as very high stiffness blade areas such as trailing edge, as compared to axial compressors, the computational challenge of contact simulations is significantly increased. Based on an explicit time integration procedure combined with Lagrange multiplier contact treatment [12], a specific contact configuration between the centrifugal compressor and its surrounding casing is investigated throughout a wide rotational speed range. It is assumed that the casing is rigid and subject to a time-invariant deformation along a two-nodal diameter shape — which means it is ovalized, see Fig. 1. The interest lies in the vibratory level of the impeller depending on its rotational speed. The modeling of the system is detailed in the first section of the paper: both the computation of the reduced order model of the impeller and the contact algorithm are outlined. The second section presents the modal analysis of the impeller and the typical linear criterion for the definition of critical rotational speeds is recalled. In the third section, convergence in time and space is addressed. Results are post-processed in the frequency domain with


INTRODUCTION
Rotor/stator interactions [1] involve a large variety of complex phenomena. In modern turbo-machinery clearances be-tween the rotor and casing are usually designed at normal operating conditions, allowing for optimal performance without significant interactions between the two parts. Nevertheless, in certain cases an engine may experience operating conditions for which clearances can be overcome and interactions can become non-negligible. As a result, contacts, such as blade/casing contacts, may occur. Subsequent interaction phenomena -which are a specific type of rotor/stator interactions -that were also observed for axial compressors both experimentally [2,3] and numerically [4] are still under investigation.
More recently, a detailed numerical model of the geometry of interest in this paper was presented in [9]. Along with other interaction phenomena such as impeller-vaned diffuser interaction [10], the understanding of interaction phenomena arising from structural contacts is a growing field of investigation. Consequences of these phenomena may result in strong energy vibration between rotor and stator which may cause damage to the structure, as seen in [11].
The proposed study focuses on blade/casing contacts within the centrifugal compressor of a helicopter engine. One specificity of centrifugal compressors lies in a complicated blade-tip geometry. Due to the high curvature of almost 90 • between inlet and outlet as well as very high stiffness blade areas such as trailing edge, as compared to axial compressors, the computational challenge of contact simulations is significantly increased. Based on an explicit time integration procedure combined with Lagrange multiplier contact treatment [12], a specific contact configuration between the centrifugal compressor and its surrounding casing is investigated throughout a wide rotational speed range. It is assumed that the casing is rigid and subject to a time-invariant deformation along a two-nodal diameter shape -which means it is ovalized, see Fig. 1. The interest lies in the vibratory level of the impeller depending on its rotational speed.
The modeling of the system is detailed in the first section of the paper: both the computation of the reduced order model of the impeller and the contact algorithm are outlined. The second section presents the modal analysis of the impeller and the typical linear criterion for the definition of critical rotational speeds is recalled. In the third section, convergence in time and space is addressed. Results are post-processed in the frequency domain with a 2D-Fourier transform. The two dimensions involved -time and space -are referred to with a specific vocabulary: nodal diameters refer to harmonics in space while harmonics refer to harmonics in time.

MODELING
Note 1: symbols and notations used in the following are all referenced in the nomenclature at the end of the paper.
Note 2: for the sake of confidentiality, all values in the space, time and frequency domains presented have been normalized.
In this section, the methodology used to treat unilateral contact between the flexible rotor and rigid casing are outlined. Contact between the two structures is initiated by a static ovalization of the casing, restricting the tip clearance equally along the entire blade tip for certain angular positions. The reference clearance is about 3% at the LE to 2.2% at the TE normalized to the TE blade height. The ovalization of the casing is implemented by using a Gaussian bell distortion on top of the reference geometry. The restriction amplitude is computed uniformly along the chord at 4.25% the TE blade height. The penetration depth will hence vary from 1.25% at the LE to 2.1% at the TE. The abradable coating that is frequently deposited on the casing is not accounted for in our study. The focus is on direct structural contacts in order to get a better understanding of possible interaction phenomena.

Bladed Disk Model
The high-pressure radial compressor stage of a modern highperformance helicopter engine is investigated. It consists of ten sectors containing each a main and a splitter blades. Contact is only managed on the main blade of each sector. Contact is managed on ten nodes evenly spaced along the chord of each main blade as pictured in Fig. 2. Consequently, a total of a hundred contact nodes are accounted for on the impeller.
The challenge in facing the impeller geometry in comparison with axial compressors lies in the varying blade tip contact normal along the chord. That, in combimation with impeller design specificities such as a low height to chord ratio and significant influence of the disk for very low frequency free-vibration modes, leads to more complex numerical simulations.

Model Reduction
The finite element model of an elementary sector of the impeller depicted in Fig. 2 contains about 120,000 degrees of freedom (DOF) and the full impeller contains more than a million DOF. Consequently, time integration on such large models would lead to cumbersome computation times. Accordingly, component mode synthesis (CMS) is used in order to reduce the dimension of the models. It is critical to avoid costly mappings between the finite element and the reduced spaces for contact management. For that reason, the CMS method used in this studywhich is based on the Craig-Bampton method [13] -features a mixed reduced space with both physical and modal coordinates.
The combination of a modified Craig-Bampton procedure [14] with cyclic symmetry [15] allows for a fast computation of a reduced order model that accounts for centrifugal stiffening. The key steps of this computation are detailed in the following. The procedure requires three stiffness matrices computed at distinct rotational frequencies: K(0), K( Ω 1 2 ) and K(Ω 1 ) as well as the mass matrix M of an elementary sector of the finite element model.
2. computation of nodal diameter matrices: each nodal diameter stiffness and mass matrix is computed from the finite element matrices, see [15].
3. centrifugal stiffening: the three stiffness matrices account-ing for centrifugal stiffening are: 4. modal reduction: for each of the three rotational speeds (Ω = 0, Ω = Ω 1 /2 and Ω = Ω 1 ), the modal reduction basiŝ Φ (i) (Ω) is: For each harmonic, the global reduction basisΥ (i) may then be computed as: where Ψ (i), * contains an orthonormal basis generated from the static and constraint modes given in Eq. (4). It yields the computation of the reduced harmonic mass and stiffness matrices: 5. recomposition of the reduced model at a given Ω: the final reduced mass and stiffness matrices are: For a given impeller, matrices K r are computed only once. Those matrices are relatively small -around a hundred times smaller than the finite element matrices of a sector of the impeller (K and M) depending on the level of modal reduction -thus making the computation of the reduced order model given in Eq. (7) particularly fast and easy. The stiffness matrix K r (Ω) of the reduced model is simply denoted K r in the following.

Contact Treatment
Time simulations are carried out based on the strategy presented in [12] which is briefly recalled here for the sake of clarity. Though compatible with the following contact treatment procedure, friction is not accounted for in this study. The considered time-marching procedure involves the explicit central differences scheme combined with a Lagrange multiplier based contact algorithm [16], assumed on 10 contacting nodes, spread equidistantly along the chord line. At each time step q, the procedure is divided into four steps: 1. prediction at time step q + 1 of the displacements u r : 2. determination of the gap function and detection of the contacting nodes of the blade. 3. correction of the predicted displacements leading to a vanishing of the gap function: Lagrange multipliers (or contact forces) and updated displacements are obtained as follows: 4. time increment.

MODAL ANALYSIS
The modal analysis of the impeller is carried out considering that the shaft is perfectly rigid and hence the DOF of the impeller that lie on the impeller-shaft interface are clamped, see  Only these first two families of free vibration modes are represented in this paper. Since the first mode families generally have the highest vibration energy participation it is assumed that the first two modal families are dominant in the impeller response even when contact occurs.
In order to identify potential critical velocities, including harmonics, it is convenient to plot the Campbell diagram in the absolute frame, as pictured in Fig. 5, where each eigenfrequency evolution with respect to the rotational speed is given by (11), a and f (i) are the absolute and relative frame rotor modal frequencies of the i-th nodal diameter at a rotational speed of Ω respectively: Since the casing is perfectly rigid in this study, interactions are considered at the intersection between f cur, the rotational speed -in the absolute frame -of the modal shape of the impeller must be null since the casing is perfectly rigid. Under these conditions, in a linear framework, a modal coincidence can be achieved for a rigid and ovalized casing.
Interaction points associated with other modal families are not considered here for two reasons: 1. because the other modal families involve higher frequency modes, the interaction points occur outside of the rotational speed range of interest, 2. as mentioned above, even when contacts occur, recent studies suggest that the first free vibration modes or families of free vibration modes are predominant in the dynamics of the structure.

NUMERICAL SIMULATIONS
The rotational speed range of interest covers the typical rotational speed range of the impeller with Ω ∈ [1.969 ; 5.951]. As mentioned above, numerical simulations on the full finite element model are not conceivable because of its very large dimension. As a consequence, a reference finite element solution cannot be obtained. Asymptotic space and time convergence are respectively assessed versus the dimension of the modal reduction basis and the time step.

Convergence
Computation times grow nonlinearly with the dimension of the reduced order model and it is thus critical to find the smallest reduction basis allowing to obtain space convergence of the simulations. Figure 6 pictures the radial displacement on the leading edge of the first sector of the impeller over the first five milliseconds of the simulation for different reduced order models. It appears quite clearly that too small reduction bases (330 and 540 DOF) lead to inaccurate results early in the interaction. Larger models (900, 1140 and 1560 DOF) lead to almost perfectly superimposed displacements. Identical observations are made for the tangential displacement of the contact nodes on the trailing edge in Fig. 7.
Similar results are obtained over a longer time interval. A reduced order model containing 900 DOF is well suited and allows for space converged results throughout the rotational speed range of interest. With 900 DOF, the gap between the ten first eigenfrequencies for each nodal diameter of the reduced order model and those of the full finite element model is less than 1%. Simultaneously, time convergence is assessed by refining the time step size and h = 10 −4 allows for convergence.

Case Study
First, a contact simulation for Ω = 5.4 is considered. Results feature displacements, speeds, contact forces and observed contact areas for which both time and frequency domain analyses may be carried out. In addition, displacement and stress fields within a sector at any time step can easily be retrieved. Figure 8 pictures the radial displacement of the trailing edge contact node during the simulation. Over the simulation, the impeller makes 50 rotations and it is visible that a steady state -a response with stable amplitude -is reached for t > 20. This is also visible in the frequency map Fig. 10, as peaks of the FFT signal lie on EO lines.
Once a steady state is reached, a fast Fourier transform (FFT) of the signal is performed, see Fig. 9. It should first be noticed that only the significant peaks are detected for f < 50. Among these peaks, highest amplitudes are observed for f ≃ 10.9, 16.3, 21.8 and 32.7. However, because of the cyclic symmetry of the structure, the first modal families exhibit a clustering in eigenfrequency, particularly around f ≃ 10.9 where the gap in frequencies for the first mode family is negligible, see Fig. 4. Consequently, the precision of the FFT may not be sufficient for a clear identification of the dominant free vibration mode at a given frequency.
It is thus of interest to analyse the results in terms of nodal diameter content. The displacement vector u(t) may be projected in the cyclic symmetry space using the Fourier matrix F: The displacement vector in the cyclic symmetry spaceû(t) contains terms associated with each nodal diameter i, i = 0, 1 . . . 5: The frequency content of nodal diameter i may then be analysed with a FFT of the componentû (i) (t). Results provided in the following feature both FFT of finite element displacements u(t) and nodal diameter displacementsû (i) (t). The combination of the analysis of both types of displacements allows for an indepth investigation of the interaction phenomenon. Figure 10 shows the radial displacement response of the trailing edge node in the frequency domain. This interaction map is drawn putting side by side the spectrum obtained for each rotational speed once a steady state is reached. Highest peaks of amplitude are depicted in red and engine order lines are pictured as black dashed lines defined by f = kΩ, k = 1, 2 . . . 18.

Interaction Map
Consistently with the casing ovalization -that implies an excitation frequency of the impeller double the rotational speed Ω -the peaks highest in magnitude are detected around the even engine order lines k = 2, 4 . . . 18. The impeller response varies significantly with Ω and some critical velocities are identified in Tab. 1. For these rotational speeds, significantly higher displacements are found.  Rotational Speed Ω  Some of the critical rotational speeds identified with a modal analysis of the impeller given in Fig. 5 do match with slightly higher amplitudes of vibration in the frequency map. For instance, around Ω ≃ Ω 5,1 = 2.12, significant peaks are visible around the 12 th and 14 th engine order lines (see highlighted areas A and B in Fig. 10). Nevertheless, none of these predicted critical rotational speeds match the highest peak of amplitudes.
Overall, the discrepancy between the linear prediction of the critical rotational speeds -mentioned in Fig. 5 -and those identified with numerical simulations is patent. For highly nonlinear phenomena such as blade-tip/casing contacts, linear interaction criteria often fail [2,4]. In particular, contact simulations on axial compressors have revealed the key role of the structure super harmonics in potential interaction areas: rotational speeds Ω = f 1 /n, n = 1, 2 . . . are privileged critical speeds. Accordingly, it seems of great interest to consider super harmonics of the first free vibration modes of the impeller in order to identify more critical rotational speeds. Super harmonics n = 2 of the two first stages of modes and the associated predicted critical rotational speeds are depicted in the Campbell diagram in Fig. 11. For the radial compressors, the number of critical speeds increases significantly: these speeds are listed in Tab. 2.
In addition, critical speeds associated with super harmonics of high free vibration modes are prone to fall within the rotational speed range of interest. Only the critical rotational speeds associated with super harmonics of the two first families of modes are

Frequency Split Study
Looking closely at the results in Fig. 10 around area (E), the odd engine order lines feature a response that was previously unknown to the authors in the speed range ∆Ω = [5.33 ; 5.58]. Firstly, a response on the odd engine order lines does not comply with the linear interaction considerations, as well as secondly, a symmetric frequency split around the 1st, 3rd and 5th engine order is observed for Ω > 5.50 (see area E in Fig. 10). A refinement of the aforementioned speed range -with an increased resolution since 150 time simulations 1 are carried out over ∆Ω -is  depicted in Fig. 12(a). This section aims at identifying the nature of this frequency split. Figure 12(a) represents the frequency domain representation of the radial displacement of trailing edge contact node of the first sector and Figs. 12(b)-(g) show the frequency content for each harmonic component. It is visible that each harmonic dominates on distinct engine order lines. More precisely, the response of harmonic i = 0 is dominant at the zero frequency line, the response of harmonic i = 2 dominates along the 2nd engine order line, the response of harmonic i = 4 is dominant along the 4th and 6th engine order lines due to the aliasing effect, while odd engine orders show an almost negligible participation.
In order to better understand the nature of the response on the odd engine order lines, two zooms over the frequency spectrum obtained for each harmonic for two different rotational speeds Ω a = 5.33 and Ω b = 5.53 are respectively pictured in Figs. 13 and 14. The spectrum depicted in Fig. 13 is obtained for a rotational speed below the one for which the frequency split occurs: the three peaks of interest on the 1st, 3rd and 5th engine order lines solely involve even harmonics.  Similarly, when looking at the spectrum given in Fig. 14 computed for the rotational speed Ω b , depicted in Fig. 10, for which the frequency split is observed -it is visible that the residual peaks only involve even harmonics.
Because of their position on odd engine order lines, these peaks cannot be linked with fundamental harmonics of free vi- bration modes of the impeller. Accordingly, their very existence is an illustration of the incorrectness of linear considerations for the prediction of critical rotational speeds. Taking into account previous results obtained with axial compressors [12], it seems reasonable to assume that these peaks stem from super harmonics influence. Supporting evidence of this assumption lies in the agreement between a few critical rotational speeds predicted for super harmonics n = 2 and n = 3 and the presented numerical results. The nature of the frequency split remains to be investigated.

CONCLUSION
This study focuses on the numerical simulation of blade/casing structural contacts for the centrifugal compressor of a helicopter engine. The use of cyclic symmetry and component mode synthesis method allows for efficient time simulations over a wide rotational speed range. The contact configuration involves a rigid ovalized casing and unilateral contact is accounted for on ten contact nodes on each main blade of the impeller.
Similarly to what was already witnessed with axial compressors [4] asymptotic space and time convergence are shown. A modal analysis of the impeller leads to the computation of its Campbell diagram and the usual linear criterion -defined by the intersection between frequency lines and engine order lines -is used to predict potential critical rotational speeds.
Results are presented both in time and frequency domains. It is first highlighted that linearly predicted critical rotational speeds do not match with the numerically predicted critical rotational speeds. Based on an in-depth frequency domain analysis as well as previous observations with axial compressors, the assumption is made that super harmonics of the free vibration modes of the impeller play a key role in the occurrence of critical regimes.
Work is in progress for the precise identification of the split frequency phenomenon observed around odd engine order lines for high rotational speeds. Also, further numerical studies will incorporate abradable coating on the casing. Also, accounting for a flexible casing within the presented strategy will allow for the detection of travelling modal coincidence [17] in the absolute frame.

ACKNOWLEDGEMENT
Thanks go to Turbomeca and Snecma for their technical and financial support. This work takes place in the framework of the MAIA mechanical research and technology program sponsored by CNRS, ONERA and SAFRAN Group.