Nonsmooth Thermoelastic Simulations of Blade–Casing Contact Interactions

In turbomachinery, it is well known that tighter operating clearances improve the efﬁciency. However, this leads to unwanted potential unilateral and frictional contact occurrences between the rotating (blades) and stationary components (casings) together with attendant thermal excitations. Unilateral contact induces discontinuities in the velocity at impact times, hence the terminology nonsmooth dynamics. Current modeling strategies of rotor-stator interactions are either based on regularizing penalty methods or on explicit time-marching methods derived from Carpenter’s forward Lagrange multiplier method. Regularization introduces an artiﬁcial time scale in the formulation corresponding to numerical stiffness which is not desirable. Carpenter’s scheme has been successfully applied to turbomachinery industrial models in the sole mechanical framework, but faces serious stability issues when dealing with the additional heat equation. This work overcomes the above issues by using the Moreau–Jean nonsmooth integration scheme within an implicit (cid:18) -method. This numerical scheme is based on a mathematically sound description of the contact dynamics by means of measure differential inclusions and enjoys attractive features. The procedure is unconditionally stable opening doors to quick preliminary simulations with time-steps one hundred times larger than with previous algorithms. It can also deal with strongly coupled thermomechanical problems.


Introduction
Turbomachinery equipment manufacturers tend to improve the aerodynamic performance by reducing operating gaps between rotating and stationary components.This leads to frequent contact occurrences, referred to as rotor-stator interaction, in Address all correspondence to this author 1 compressor and turbine stages of modern aircraft and helicopter engines [1].To reduce the possibly harmful effects of such events on the global dynamics and subsequent damage, casings are commonly coated with abradable materials [1].This strategy has proven efficient, however the resulting dynamics, coupling contact and frictional heating, is still not well understood and modelled.
Various developments have been dedicated to the improvement of turbomachines design to account for blade-tip rubbing, modal interaction or whip-whirl motions [1].In these works, contact is dealt with using either a regularized approach [2] or the forward Lagrange multiplier formulation [3,4].The regularized approach consists in modeling the contact between the rotor and the stator through a stiff spring.Though very simple to implement, the additional and somehow artificial spring leads to stiff numerical problems which are prone to stability issues [5,6].The Lagrange multiplier formulation preserves the nonsmooth framework induced by contact conditions but uses an explicit procedure in time for which stability is guaranteed only for sufficiently small time-steps.This is particularly problematic since thermal fluxes generated by frictional heating seem to play a significant role in rotor-stator interactions [7] and are thus considered in the model.Indeed, explicit schemes are not well-suited for the heat equation as they lead to stiff numerical problems [8].More specifically, the stability is typically governed by a relationship of the form t Ä x 2 =2, implying that a refinement of the space discretization has drastic consequences on the necessary refinement in time.
A simple thermo-mechanical model of a bladed-disk sector section 1 is first considered in this work.The coupling is two-fold: from heat diffusion to mechanics via thermal expansion, and from mechanics to heat through a frictional heating law.The governing equations are solved via dedicated nonsmooth solvers which rely on an implicit numerical scheme [6] as described in section 2. Time-domain histories, systematic comparison with the forward Lagrange multiplier method as well as sensitivity analysis to time-step are discussed in section 3. Finally, the methodology is adapted to an industrial-size model in section 4.
1 Description of the simplified model Due to the high computational costs induced by non-linear phenomena in large-scale problems, model-order reduction methods are often implemented for realizable computations, as exposed in section 4.However, the simple one-sector model illustrated in fig. 1  and prescribed temperature of 25 ı C) on its internal bore and contact nodes on the blade tip.The structural dynamics is governed by an equation of the form where u is the vector of generalized displacements, Â is the vector of generalized temperatures, M uu , C uu , and K uu are respectively the mass, damping and stiffness matrices stemming from the finite element discretization.The temperatureinduced expansions are accounted for through the coupling matrix K uÂ .The external forces are split into the external load f u and the contact forces f u c .In a more elaborate model, f u could include centrifugal or aerodynamic effects while contact could be induced by the static deflection or the vibrations of the stator.In this work, the load f u , described further, activates contact.To every contact node i D 1; : : : ; 7 corresponds a gap function g i .u/measuring the algebraic distance to the stator and an impulsive contact force i in the inward-pointing normal direction.The following so-called Signorini conditions complements eq. ( 1): The first two inequalities enforce the non-penetration and the non-sticking condition, respectively.Equality (2c) translates the fact that the gap is open (g i > 0) if and only if the associated contact force is zero, except in the very specific case g i D 0 and i D 0. The seven Signorini conditions (one for each contact node) are usually gathered in the compact form 0 Ä g ?0 where g D OEg 1 ; : : : ; g 7 > , D OE 1 ; : : : ; 7 > and the operators Ä, ? and are defined component-wise.In addition to Signorini conditions, as in any space semi-discretized contact formulation, an impact law is required to ensure uniqueness of the solution [9].In order to mimic lasting contact phases, we choose an inelastic Newton impact law: y C i D 0 where y C i D P g i .u.t C // denotes the normal velocity right after the impact time.In the tangential direction, f u c includes Coulomb friction, therefore proportional to the normal forces.This can be written as f u c D C u NT where C NT is a rectangular matrix transferring the normal impulses in the local coordinates to the generalized coordinates in the normal direction, as well as in the tangential direction through the Coulomb friction coefficient.In the present case, the local nodal reference frames at the contact nodes are oriented according to the global cylindrical directions so that C u NT takes the form where R is the stator radius and the angular velocity-the contribution of the blade vibration in the tangential contact velocity is neglected.Note that because of the high tangential contact velocity, there is no stick-slip transitions and friction is therefore not source of nonsmoothness here.Also, eq. ( 1) does not include centrifugal stiffening or spin softening.They could however be incorporated in a straightforward manner.The thermal dynamics is governed by an equation of the form where C ÂÂ and K Â Â are the heat capacity and heat conductivity matrices.As for the structure, the external load is separated into external fluxes f Â and frictional heating due to contact f Â c .Each of the seven frictional heat fluxes Â i follows a simple law . This means that each contact occurrence generates a heat flux proportional to the normal contact force.This simple law [10] is sufficient to illustrate the capabilities of the proposed numerical method.However it is worth mentioning that it might induce temperature discontinuities.Indeed, if one i is impulsive, eq. ( 4) shows that P Â will also be corresponding to a discontinuous Â. Integrating more complex laws is straightforward.Also, in contrast to other available investigations [11,12], contact-induced wear in rotor-stator interactions are not considered in the present model.
Altogether, eq. ( 1) and eq.( 4) are compactly recast in the form and The matrices M and K were obtained using SOLID226 coupled-field elements in ANSYS R .The damping matrix verifies C uu D 10 5 K uu .Contact is activated via f u , chosen here as a sinusoidal function of frequency 33 Hz and amplitude 100 N, in the radial direction and on all contact nodes.The frequency corresponds to D 2000 rpm, assuming the blade hits the stator once per rotation.For comparison, the first flexural mode of the blade has a frequency of 1300 Hz.The contact frequency could also be due to the vibration of the blade, as in section 4.
A TA6V titanium alloy, widely used for the production of aircraft compressor components, is considered.The material properties, estimated from several alloys properties, read as follows: density 4430 kg m 3 , Young's modulus 110 GPa, Poisson's ratio 0:3, heat capacitance 520 J K 1 kg 1 , heat conductivity 6:7 J K 1 m 1 , dilatation coefficient 9 µm m 1 K 1 .The coefficient of frictional heating is ˛D 0:1 and the coefficient of friction in Coulomb's law is D 0:15.
2 Nonsmooth numerical methods 2.1 Formulation using Differential Measures Equation (1) may seem to be an Ordinary Differential Equation (ODE).More rigorously, it should be described by means of weaker mathematical objects such as distributions or measures.Indeed, f u c contains impulsive terms: every time a gap closes, the corresponding contact node instantly looses all its kinetic energy.This instant loss of kinetic energy is mirrored by an impulse that is proportional to a Dirac distribution (which is not a function).Equation ( 1) is therefore an abuse of notation in the context of nonsmooth dynamics.In this subsection, the model is reformulated in a sounder mathematical manner, using measure differential inclusions, making the use of efficient dedicated numerical methods possible.
Introducing differential measures dv and dÂ [6] and the Lebesgue measure dt, eq. ( 5) can be expressed as In short, the prefactors of dt are smooth terms, while dv and dr are measures which, when integrated between two arbitrarily close times, can yield non-zero values: they contain impulses, which are proportional to Dirac delta distributions.The measure dr stores all contact forces and frictional heating.Loosely speaking, the reader who is not familiar with measures can compare eq.(8a) to eq. ( 5) multiplied by dt.Using convex analysis, the Signorini conditions eq. ( 2) together with the inelastic impact law mentioned above can be formulated using the proj operator and the tangent cone T in R C [6] as: where the equality is componentwise, and the vector y stacks the normal contact velocities OE P g 1 ; : : : ; P g 7 > .The contact impulses are related to y through Let us inspect eq. ( 9) for the i th unilateral constraint.By definition, the tangent cone reads as and the following can be said.
During a free flight, g i .u/> 0 so y C D proj.R; y / D y and from eq. ( 11), i D 0 follows: there is no reaction impulse and velocity is continuous.At the beginning of a contact phase, g i .u/D 0 and y i Ä 0 so that y C i D 0 and the impact law is recovered.The corresponding impulse is D O W 1 y 0. In the middle of a contact phase, g i .u/D 0 and y D 0 so that y C D 0 and there is no reaction impulse.At the end of a contact phase, g i .u/D 0 and y i < 0 so that y C D proj.R C ; 0/ D 0 and i D 0. The contact force vanishes and the gap opens.
To summarize, the problem is governed by eq. ( 8), eq. ( 9), the frictional heating law f Â c D C Â NT as well as the geometric relationships f u c D C u NT and y D r u g P x.

Time discretization
We now proceed with the Moreau-Jean discretisation [13,14] of eq. ( 8), using an linearly implicit solver.For n 2 N C , let t n denote the nth time-step such that t n D nh where h is the time-step.Integrating eq. ( 8) between t n and t nC1 yields: Smooth quantities are discretized by means of the Â -method, where to avoid confusion with the temperature Â, the numerical parameter is denoted by : and similarly for x and f, so that eq. ( 12a) becomes We now introduce a predictor step, corresponding to the dynamics without the unilateral constraints.The purpose of this step is to estimate whether the contact status is about to change in the next time-step.Two new quantities are defined: the predicted gap chosen as Q g nC1 D g.u n C hv n =2/ and the predicted velocity The predicted normal velocities in the local coordinates are given by Q y nC1 D r u g Q v nC1 .Equation ( 9) is discretized as Then, the reactions in the global frame p nC1 are computed from the normal contact reactions in the local frame nC1 D O W 1 .ynC1 y n / through the geometric relationship p nC1 D C u NT nC1 .The velocity is finally updated by adding the effects of contact to the contact-less prediction Q v nC1 : 3 Validation using the simple model In this section, the reference solution is generated via the methodology described in the previous section, with a time-step h D 1=3 10 6 .The error of a data series is defined as the integral of the absolute value of the difference with the reference solution, divided by the integral of absolute value of the reference solution, over the time interval OE0 s; 1 s.For example, a time series of the temperature of the first node Â 1 shows the error Although imperfect, this definition is useful to quantify discrepancies in the remaining.

Time-step sensitivity analysis
Figure 2 displays the time evolution of the radial displacement, the temperature and the normal contact force for the middle of the blade tip, that is, contact node 4 in fig. 1, and different time-steps.For h D 10 4 s, the error is about 3 %, while temperature and contact force show 16 % and 25 % errors, respectively.After one second, the temperature is underestimated by about 5 %, resulting from an underestimation of the normal contact force.This can lead to significant errors on longer simulations, but might be sufficient if contact events are expected to be of short durations.With h D 10 5 s, the errors reduce to 0:21 %, 0:18 % and 1:72 %: the curves cannot be distinguished from the reference.Due to the strong thermomechanical coupling, the temperature oscillates at the loading frequency of 33 Hz.Heat accumulates and the temperature globally increases.As a result, the blade expands through the coupling term K uÂ and the gap tends to reduce while the normal force tend to increase.A close-up view on the last time range 0:9 s to 1 s is provided in fig. 3. Note that the simulated temperatures of thousands of degrees go far beyond the material capabilities: that is an obvious shortcoming of the simple chosen model, which does not include nonlinear effects (other than contact) nor damage, but is sufficient for the present scope: validating the methodology.A more realistic model is presented in section 4.
The normal contact force, always non-negative, is seen to be non-zero only when the gap is closed, in accordance with the Signorini conditions eq. ( 2).Even with relatively large time-steps, the position exhibits clear kinks when the gap opens and closes, corresponding to velocity discontinuities induced by the inelastic impact law, features which would not be accurately captured via a classical time-domain integration scheme dedicated to smooth dynamics.
The proposed algorithm does not strictly prevent penetration, as illustrated in fig. 4.However, the residual penetration, which is small compared to the geometrical tolerances of the machine components, tends to zero as the time-step reduces.Also, once contact is activated, penetration never increases.
To conclude on the convergence analysis, the error is depicted as a function of the time-step in fig. 5. Given that the Carpenter's algorithm is unstable with h > 5 10 7 s (see section 3.3), the proposed method enable the use of time-steps larger by several of orders of magnitude.

Effect of coupling on contact geometry
The previous figures were depicted for contact node 4 (middle of the blade tip), which is the first one to close the gap.In this section, it is shown that thermomechanical coupling can have an effect on the contact geometry.Figure 6 compares the simulation results between node 4 and its neighbour node 3, as defined in fig. 1.
For both nodes, the contact force is non-zero if and only if the gap is closed.The middle node touches the stator on a longer interval than node 3, resulting in a normal force of larger amplitude, and higher temperature due to frictional heating.At the beginning, all boundary nodes impact the stator.After a few contact occurrences, the middle of the blade has expanded so much due to heating that the nodes located at the end of the blade tip (nodes 1 and 7) separate from the stator.This is  opens, the trajectory of node 3 displays a small kink because of the structural coupling M uu and K uu .In contrast, when thermomechanical coupling is ignored, the discrepancies between normal forces and displacements tend to reduce, see fig. 7.  The effect of temperature on the normal force appears clearly when comparing the bottom plot with that of fig.6: thermal expansion further magnifies gap closure and attendant contact force.

Comparison with Carpenter's scheme
Results are now compared to the ones obtained with Carpenter's forward Lagrange multiplier method.Figure 8 shows that displacements and normal forces match well, with overall relative errors of 0:82 % on nodal displacement, 5:1 % on temperature and 7:3 % on contact force.Figure 9    The major benefit of the proposed method is its unconditional stability while Carpenter's scheme is unstable whenever h > 4:6 10 7 s.The presented nonsmooth approach enables the use of time-step hundreds of time larger with similar error.

Industrial application
We now show that the above simulation methods are not limited to simple models such as section 1 by using a more realistic model of an axial compressor sector with a twisted blade.

Model
The model illustrated in fig. 10 comprises 5686 nodes and a total of 18768 DOFs (including 4692 thermal DOFs), with 9 contact nodes located on the middle plane of the blade tip.In order to obtain reasonable computational times (of the order of The contact simulations were performed by periodically forcing the blade on its first flexural mode at frequency 507 Hz to reflect an aerodynamic excitation.The forcing amplitude was chosen in order to obtain approximately a steady-state radial displacement amplitude of 10 4 m on node 1 without contact.

Sensitivity to model parameters
To emphasize the effect of friction-induced heating during rotor-stator interactions in turbomachinery, a contact simulation is performed with and without frictional heating.The structure is initially at rest.The results are depicted for node 1 in fig.11.The response curves oscillate very fast because the simulation time is much larger than the excitation period.During the first contacts, little difference is observed between the two responses.Then, the gap in the simulation with heating tends to reduce, while the contact force increases.At the end of the simulation, discrepancies of the order of 25 % on the displacement and 30 % on the contact force are predicted.This qualitative analysis shows that thermal effects have noticeable influence during the contacts, even for short durations.
The nonlinear nature of the system, stemming from the contact conditions, can be emphasized by doubling the excitation force amplitude, see fig.12.
Magnification on the time interval [0 s,0:05 s] is provided in fig.13.As expected, the first contact occurs earlier when the external load is magnified, due to higher vibration amplitudes.The higher excitation level also leads to an increased contact force and thereby to higher temperatures, inducing slightly smaller displacements.The temperature seems to diverge with the doubled load, probably because frictional heating generates more heat than the blade can diffuse.A thermoelastic instability emerges: it could advance the wear of machine components or even lead to more extensive damages.

Sensitivity to time-step
For this model, the stability condition required by Carpenter's algorithm is h D 5 10 7 s.The Moreau-Jean and Carpenter algorithms were implemented with this time-step.The response for node 1 is displayed in fig.14.Again, they are in very good agreement.
Figure 15 shows results for various time-steps with Moreau-Jean's procedure.The results with h D 10 5 s and with h D 5 10 7 s are almost identical, meaning that convergence is achieved for h D 10 5 s.Therefore, the stability condition of Carpenter's algorithm is unnecessarily limiting.Moreover, for the same time-step h D 5 10 7 s, Carpenter's algorithm took about twice as long as Moreau-Jean's algorithm in our implementation-this might be due to implementation difference and may not be true in general.More computation times are reported in fig.16.All computations were performed using MATLAB R on a 4-core computer (2:3 GHz clock frequency).

Conclusion
Conventional time-domain integration schemes to simulate rotor-stator interaction thermomechanical formulations face numerical stability issues: regularized contact constraints lead to numerical stiffness while explicit solution methods such as Carpenter's forward Lagrange multiplier scheme are not suited for the heat equation.The methodology proposed in the present paper, based on the Moreau-Jean procedure for measure differential inclusions addresses these shortcomings.It is shown to be very robust on two models of different sizes.While Carpenter's stability requires time-steps of the order of

NomenclatureM
uu , C uu , K uu Structural mass, damping and stiffness matrices C ÂÂ , K ÂÂ Heat capacity and heat conductivity matrices K uÂ Thermoelastic coupling matrix I Identity matrix 0 Zero matrix u, Â Nodal displacement and temperature vectors x Generalized degrees of freedom vector v Generalized velocity vector M; C; K Generalized mass, damping, stiffness matrices f u , f Â Nodal force and heat flow vectors Normal contact force (Lagrange multiplier)

Fig. 2 :
Fig. 2: Sensitivity to time-step in terms of radial displacement u r , temperature Â and normal contact force at node 4. Reference [ ], h D 10 4 s [ ] and h D 10 5 s [ ].

Fig. 5 :
Fig. 5: Error as a function of the of time-steps at contact node 4 for one second of simulation: u r [ ], [ ] and Â [ ]
is a close-up of fig. 8.