THERMOMECHANICAL COMPONENT MODE SYNTHESIS FOR BLADE CASING INTERACTION PREDICTION

Increasing the efﬁciency of turbomachines is a major concern as it directly translates into lower environmental impact and improved operational costs. One solution is to reduce the blade-casing operating clearance in order to mitigate aerodynamic losses at the unavoidable cost of increased structural unilateral contact and friction occurrences. In centrifugal compressors, the dynamic behaviour of the structures interacting through unilateral contact and friction is not yet fully understood. In fact, the heat generated during such events may affect the dynamics through thermal stresses. This paper presents a complete thermomechanical modelling strategy of impeller rotor and casing, and of blade-tip/casing contact events. A fully coupled thermomechanical modal synthesis technique is introduced and applied to turbomachinery-related models. The blisk is reduced via a hybrid modal synthesis technique combining the Craig-Bampton method and the characteristic constraint mode method. The casing model is reduced using an axisymmetric harmonic modal synthesis. Both strategies involve thermomechanical modes embedding thermal dilatation effects. The contact modelling algorithm is then introduced. It handles unilateral contact and friction occurrences together with heating effects. This algorithm uses the above mentioned reduced-order models as input data to avoid CPU-intensive simulations. The results show that the thermomechanical behaviour of the structures is well preserved by the reduction strategy proposed. Contact simulations on simple cases show qualitative results in accordance with expectations. Further work is needed in order to validate the strategy based on experimental results. However, this methodology opens the way to extended multiphysics simulations of contact events in turbomachinery.


INTRODUCTION
With a growing need for performance and power-to-weight ratios of aircraft and helicopter engines, designers make their best to improve the efficiency of every component.A solution is to reduce blade tip-casing operating clearances in order to optimize pressure ratios in the compressor and turbine stages [1].Due to such clearance reductions, the occurrence of unilateral and friction contact phenomena between blades and casing becomes a common situation during the lifetime of the engines.It is now recognized that such phenomena have to be accounted for in the design in order to avoid premature maintenance operation and therefore excessive operational costs.Several investigations, both numerical [2,3] and experimental [4,5] are devoted to the modelling and understanding of such events.A detailed summary [6] lists recent studies regarding material, thermal and mechanical aspects of rotor-stator rubbing.
Recent experimental campaigns [4,5] have shown the importance of thermomechanical coupling effects during the interaction.However, very little documentation is available [7,8] and thorough numerical investigations have yet to be undertaken.Still, thermomechanical contact models exist [9,10] and related issues are explored in other industrial fields, for instance braking systems [11] or nuclear reactors [12].This paper deals with the development of a reduction strategy taking into account the thermomechanical response of a bladed-disk and its matching casing during rubbing events, while maintaining low computational costs.The proposed strategy is based on [13], and takes advantage of cyclic and axisymmetric structures properties.A time-marching algorithm is also presented.It has been developed in order to model both contact forces and heat flows during rotor-stator interactions.First simulation results are then presented.

THERMOMECHANICAL COMPONENT MODE SYN-THESIS 1.Linear thermoelastic solids
This paper focuses on the modelling of linear thermoelastic solids.Other physical phenomena are neglected, as well as the variation of the material properties with temperature.Based on the first and second laws of thermodynamics for thermoelastic solids, the Clausius-Duhem inequality writes where ψ = e − T s is the Helmholtz energy, σ , the stress tensor, ε, the strain tensor, ρ, the material density, s, the entropy, and T , the absolute temperature.Vector q is an external thermal flux.Under the assumption of small perturbations of the system, and considering ψ only dependent on the strain tensor and the temperature, the state laws of thermoelasticity are formulated using Clausius-Duhem's inequality.A first order Taylor series expansion of the laws of thermoelasticity and Fourier's conduction law gives the expressions of the linear thermoelastic behaviour law and of the heat equation where Ē stands for the tensor of elastic moduli and K for the thermal expansion tensor.T 0 is the reference temperature, c is the material heat capacity, λ its thermal conductivity, and r an external heat source.From the variational formulation of the above Partial differential equations, the following spatially discretized formulation is defined: It is commonly used in commercial finite element software.The coupling of thermal and elasticity terms is dictated by the K uθ and C θ u matrix blocks, respectively standing for dilatation and thermoelastic damping effects.The thermoelastic damping is often neglected in mesoscopic scaled systems.In the case of microscopic systems, such as microelectromechanical systems (MEMS), it plays a greater role and has to be accounted for [14].

Craig-Bampton thermomechanical reduction
In order to perform a Craig-Bampton coupled reduction of the system in Eq. ( 5), matrices need to be partitioned into internal (• i ) and boundary (• b ) degrees-of-freedom (DOF).This applies to both thermal (θ θ θ ) and mechanical (u) DOF, and leads to the following matrix basis change: As in the standard Craig-Bampton technique [15], the method proposed in [13] uses both static and dynamic modeshapes in the reduction basis.Static modeshapes are calculated using unitary nodal displacements and temperatures successively imposed on each boundary node while blocking the other boundary degreesof-freedom: where Ψ Ψ Ψ ib is defined as This reduction basis contains thermomechanical coupling through the K uθ blocks.Calculating the boundary constrained modes of the system is less obvious.Two issues can be raised: • Which modeshapes should be incorporated to capture the transient thermomechanical behaviour of the system?• How could these modes be evaluated at a low computational cost?
In a standard mechanical reduction method, some eigenvectors would be selected based on a specific criterion, for example the eigenvalue.The reduction is then operated by only retaining the modeshapes associated with the lowest eigenvalues.In the present work, such a strategy would be too expensive.In fact, the eigenvector calculation of thermomechanical systems is ill-conditioned.
To overcome this situation, a formulation of a coupled reduction basis from the uncoupled sub-problems is proposed in [13].The dynamic reduction matrix is built by means of purely mechanical vibration modes Φ Φ Φ uu m , as well as purely thermal heat conduction modes Φ Φ Φ θ θ m .The coupling is operated using so-called "thermaldriven modes" that may be interpreted as the displacements induced by thermal modeshapes due to dilatation effects.The dynamic reduction basis is then The full reduction matrix is built by concatenating the static and dynamic reduction matrices in Eq. ( 7) and (10).The above described thermomechanical reduction is based on the approach developed by Craig and Bampton, involving fixed interface and displacement static modes [15].However, other component mode synthesis methods could be considered for reduction purposes.

TURBOENGINE IMPELLER
Turboengine impellers feature peculiar geometries leading to large finite element models.In order to reduce such models, a double reduction methods is implemented.

Proposed strategy
The reduction strategy developed is illustrated in Fig. 1.The first step consists in the application of the above thermomechanical Craig-Bampton reduction on the finite element model of a sector of the impeller, see Fig. 1(a).For this reduction, the retained degrees-of-freedom are the thermal and mechanical DOF on the cyclic symmetry boundaries of the sector as well as DOF of nodes on the blade tip.This reduction leads to the construction of a sector superelemen (SE) as shown in Fig. 1(b).This superelement is then replicated to reproduce all sectors of the blisk.Later on, superelements are assembled using displacement compatibility conditions between adjacent sectors.This operation already cuts half of the degrees of freedom of the attach faces, further reducing the model as illustrated in Fig. 1(c).However, since the amount of degrees-of-freedom on these faces is in general very high (especially when considering impeller rotors), a second reduction pass has to be performed.
The second reduction step lies on the interface or characteristic constraint modes [16,17] reduction method.It is used to reduce interface DOF resulting from the assembly of superelements and is based on the eigenvectors of the interface subsystem.In this particular case, interfaces are the attached faces remaining once sector superelements have been assembled.This method was applied to thermomechanical problems using the thermal driven modes Eq. ( 9).Compared to the reduction of a full blisk finite element model, the proposed strategy is expected to be more efficient.In fact, the matrices handled during calculations are much tinier, leading to both lower RAM and CPU needs.

Application model
The finite element model displayed in Fig. 2 and used to create the blisk superelement was built using ANSYS software and Solid226 coupled-field elements.In order to obtain a CPU efficient model, a simplified geometry as well as a coarse mesh were studied.Each blade of this model possesses a boundary node (i.e.retained during reduction)on its tip with three nodal displacements and a nodal temperature.The system is constrained both mechanically and thermally on its internal bore.All other degrees-of-freedom are considered free of any constraint.Thus, external boundaries are perfectly insulated.Such conditions do not represent actual turbomachinery applications.However, test rigs dedicated to the study of bladecasing contact events on actual engine components work under low pressure conditions in order to avoid aerodynamic perturbations [8].Such an environment leads to low convection exchange between the remaining fluid and the contacting structures.The constant temperature on the internal bore models the thermal exchange with a shaft.

Modal analysis
As stated before, it may be difficult to calculate fully coupled thermomechanical modes of large finite element models.Moreover, a thorough time-stepping comparison between a reduced and a full model is a tedious task.In order to check the validity of the proposed reduction strategy, the eigenvalues of the full and reduced models thermal and mechanical subproblems are first compared.The comparison is performed with respect to two parameters: the number of modes retained during the first and second reductions, respectively.Both parameters vary from 5 to 30 by steps of 5.As observed in Fig. 3, the eigenfrequencies show a maximum deviation under 1 %.The convergence is sensitive to both reduction steps and depends on the considered mode.An identical modal comparison is performed on the thermal block of the complete and reduced models, displayed in Fig. 4, and leads to the same conclusions.Figures 3 and 4 show that the reduc-

Transient analysis
Time domain transient analysis is now performed under various external load cases so as to ensure the quality of the thermomechanical coupling following the second reduction step.This transient analysis is performed on a blisk superelement with 10 structural interface modes and 10 thermal interface modes.A  Nodal displacements and temperatures are reported in Fig. 6.The displacements of boundary nodes seem to be well described under such mechanical loads.Discrepancies arise in nodal temperatures though.They are probably related to the coupling approximation.However, the effect of a mechanical load on boundary temperatures is practically zero in both cases.A second validation is com-

CASING THERMOMECHANICAL MODAL SYNTHESIS
In order to perform relevant contact simulations, it may be interesting to model the casing as an elastic component.As this structure exhibits an axisymmetric geometry, a modal harmonic synthesis of both nodal displacements and temperatures is proposed.

Axis ans cyclically symmetric structures
Axisymmetric structures have geometrical properties such that problem dimensions in a Finite Element framework can be reduced to a minimum.In particular, such structures may be represented with a continuous spatial Fourier transform.It is then possible to travel between physical (x) and Fourier (x) coordinates via the transformation where β is the angular position and k is the harmonic index.This approach can be exploited with thermal and thermomechanical models.Therefore, a thermomechanical equivalent of Eq. ( 11) is An analogous formulation can be employed for cyclically symmetric structures by replacing the continuous spatial Fourier transform with a discrete one In this formulation, the angular position is replaced by the sector number n multiplied with α, the sector angle.

Modal analysis
According to [18], the eigenproblem of the whole structure is equivalent to the set of independent eigenproblems associated to each harmonic index k.For each of them, the uncoupled thermal and mechanical modes are calculated.The thermal-driven modes in Eq. ( 9) are inserted to relate both physics.Therefore, the modal shapes for a given k are The method provides access to the thermomechanical modal response of an axisymmetric structure based on a single crosssection.Using the previous modal basis onto which the dynamics of the casing is projected, Fourier coordinates are expressed as Finally, via Eq.( 12), physical coordinates are expressed as a linear combination of time-dependent contributions of modeshapes: Equation ( 16) is displayed graphically in Fig. 8.

FIGURE 8:
Casing reduced-order model strategy

Application: impeller casing model
The casing on which the method was applied is displayed in Fig. 9

ROTOR-STATOR INTERACTION SIMULATION 4.1 Continuous formulations
A modelling of the thermomechanical behaviour of a contact interface has been proposed by several authors [9,10,19].This formulation completes the equations governing the contacting solids by introducing an interface potential ψ I and another Clausius-Duhem inequality governing the thermodynamics of the interface: In this case, the interface law only accounts for unilateral contact, Coulomb friction and frictional heating.The phenomena are described by the normal (f N , u N ) and tangential (f T , u T ) contact forces and relative displacements, the interface absolute temperature Θ Θ Θ I and entropy s I ; Θ Θ Θ 1 and Θ Θ Θ 2 are the absolute temperatures of the contacting bodies, q c,1 and q c,2 are the contact heat flux.Other phenomena may be considered in the contact law, such as adhesive or abrasive wear [20].

Resolution scheme
Contact problems are mainly solved using two strategies: time-marching methods or harmonic balance methods.Direct integration methods allow to study transient phenomena.Such methods are already in use to study rotor-stator contact problems, e.g.central finite differences [21,3].Harmonic methods are used to model perioddic steady-state phenomena and have been applied to study blade-root / disk contacts [22] and friction phenomena in general [23] .
Regarding thermomechanically coupled problems, two approaches can be highlighted, either tightly or loosely coupled.In the first case, the fully coupled thermoelastic system is solved in one step.In the other case, successive resolutions of the mechanical and thermal problems are performed.A recent evaluation of both methods [12] shows that simulation performances greatly depend on the size and type of problem explored.In particular, for thermomechanical contact problems, the tight coupling seems to be more efficient.
In this study, a strongly coupled resolution is carried out, in accordance with the reduction strategy presented.In order to solve the dynamic system, a scheme based on the weighed residuals method is used.It matches the central finite difference method used in [2], when the integration parameter α is of 1/2.However, such a parameter value leads to an unconditionally unstable scheme due to the thermal part of the model.For values above 1/2 the scheme is conditionally stable.Practically, α is chosen close to 1/2 to avoid excessive numerical damping.The time discretisation leads to where ∆t is the time increment and n is the increment number.

Contact event handling
In structural dynamics, these problems are often tackled using penalty or Lagrange multiplier methods [24,25,2], the latter having the advantage of preventing residual penetration between the contacting bodies.The algorithm used is based on a forward Lagrange multiplier prediction-correction [24,2].The gap g p separating the two structures is first calculated without contact conditions.If a negative gap is detected, Lagrange multipliers (i.e.normal contact forces) are calculated according to The C N,p matrix contains the predicted direction of the normal contact force.The C NT,p matrix stores the predicted direction of the contact force taking into account both normal and tangential forces.The displacement corrections δ u are then calculated according to In the present work, the friction law is regularized using where µ is the friction coefficient, λ λ λ N , the normal contact force, and uT , the relative velocity between a pair of points located on the rotor and on the stator.The regularization is controlled by ε.
Thermal considerations in a contact simulation lead to the formulation of heat flows generated at the contact interface.They are induced by heat creation at the interface (depending on the contact law) and by conduction between the two structures.Therefore, they are very much influenced by the surface roughness of the two structures (impacting the thermal contact conductance), the friction coefficient and local temperatures.The contact heat flows take the following form: where ν 1 and ν 2 are the thermal contact conductances of the bodies coming into contact.Based on the displacements correction Eq. ( 20) in the purely mechanical case, the following formulation is proposed for the thermomechanical case, where C NT θ ,p gathers the participations of the normal and tangential forces as well as the contact heat flux.Thus, this strategy consists in an isothermal prediction step followed by a thermomechanical correction step.The

Excitation-driven contact
This first simulation is performed with a undistorted rigid stator.The rotor spins at 2000 RPM and is excited with a 1000 Hz sinusoidal force.This excitation force is imposed on one blade, in order to close the gap between the two structures and generate a radial contact between the two structures, as illustrated in figure 11.The excitation level increases along an exponential law.

Displacement-driven contact
This second simulation is performed with a flexible rotor and a rigid stator displaying a one diameter shape.Thus, the initial contact area is localized along the circumference of the stator.At t = 4 s, the gap is opened to observe the free behaviour of the rotor.As before, the rotor spins at 2000 RPM.Simulation results are gathered in Figs 14 and 15.It is observed that the three blades are rubbing on the casing due to the rotation.Therefore the tip temperatures increase, but with lower magnitudes than in the previous simulation.The normal displacement exhibits a white-noise like frequency spectrum due to successive contact occurrences.The tangential displacement spectrum shows a major participation of the first eigenfrequency of the rotor.

CONCLUSION
Thermomechanical model reduction methods have been proposed and applied on turbomachinery related structures in order to provide efficient models for thermomechanical rub events simulations.These methods show accurate modal characteristics, compared with reference models, while displaying a minimal amount of degrees-of-freedom.Results of time integration simulations  show that the thermal behaviour and thermomechanical coupling are well described by the reduction methods, thus allowing to link both physics during simulations.The currently developed contact algorithm, based on a forward increment Lagrange multiplier method is able to capture mechanical and thermal phenomena emerging during contact events.Further studies are however needed to extensively compare the results with other numerical tools and experimental measurements.Applying such methods on large finite element models may lead to error propagation issues and instability of the reducedorder models.Therefore, such applications require a careful implementation in order to obtain light but trustworthy models.However, low number of remaining degrees-of-freedom, the CPU cost remains limited when compared to initial full models.
Current developments attempt to extend the proposed approach to industrial models, opening the way to multiphysics simulations of rotor-stator rub events in turbomachinery.

ACKNOWLEDGMENT
The authors are grateful to Safran Helicopter Engines for providing the financial support for this project, and for giving permission to publish this work.

4 M
od e nu m be r (f re qu en cy or de re d) modes in first reduction

FIGURE 3 :
FIGURE 3: Convergence analysis for the mechanical subproblem

4 MFIGURE 4 :
FIGURE 4: Convergence analysis for the thermal subproblem tion strategy proposed allows to reduce models while preserving modal behaviours.

FIGURE 6 :FIGURE 7 :
FIGURE 6: Blade tip nodal quantities under tip load: CB ( ) and CB+IM ( ) pleted with heat flows prescribed on the blade tip boundary node.The results are presented in Fig.7.Nodal temperatures show a good correlation between the two compared models.Moreover, nodal displacements induced by the thermal loads are also well (a).Red dots indicate fixed degrees-of-freedom.As an example, Fig. 9(b) displays the first three nodal-diameters modeshapes.

FIGURE 10 :
FIGURE 10: Algorithm organisation corresponding algorithm is presented in Fig. 10.

FIGURE 11 :
FIGURE 11: Contact Configuration: Stator ( ) and Rotor ( ) At t = 4 s, the excitation force is stopped.Responses presented in Figs 12 and 13 show the gap closing due to the increasing excitation level.At t ≈ 0.3 s, the gap closes.When the external excitation stops, nodal displacements converge towards a non-zero value due to the thermomechanical coupling.The spectrogram of the excited node axial displacement in Fig 13(a) shows that the frequency content of the displacement is directly related to the excitation

FIGURE 13 :
FIGURE 13: Time-frequency analysis of blade 1 tip displacement frequency.As soon as contact occurs, harmonics of the excitation frequency are generated due to the non-linear behaviour.When the excitation stops, the gap reopens and the frequency content of the nodal displacements follows the natural frequencies of the system.In this very case, the first bending mode frequency is 1306 Hz.The spectrogram of the tangential displacement of the contacting node supports the above statements.As contact is maintained during the simulation, the contacting node temperature increases rapidly as shown in Fig12(c).As soon as the gap reopens, heat generation stops and the blade tip begins to cool down due to conductive effects within the blade.The temperature levels displayed highly depend on the contact parameters: friction coefficient, relative speed, and contact conductance.