Thermomechanical model reduction for efficient simulations of rotor-stator contact interaction

Turbomachinery rotor-stator unilateral contact induced inter-actions play a growing role in lifecycle analysis and thus motivate the use of accurate numerical prediction tools. Recent literature conﬁrmed by ongoing in-house experiments have shown the impor-tance of thermomechanical coupling effects in such interactions. However, most available (possibly reduced-order) models are re-stricted to the sole mechanical aspects. This work describes a reduction technique of thermomechanical models involving unilateral contact and frictional contact occurrences between rotor and stator components. The proposed methodology is grounded on Guyan and Craig–Bampton methods for the reduction of the structural dynamics in conjunction with Krylov subspace techniques, and speciﬁcally the Craig–Hale approach, for the reduction of the thermal equations. The method has the capability to drastically reduce the size of the model while preserving accuracy. It stands as a reliable strategy to perform simulations of thermomechanical models with localized mechanical and thermal loads.


INTRODUCTION
As the requirements for engine performance become more and more stringent, multiphysics and nonlinear phenomena are paid a continued attention in the turbomachinery industry. Of interest in this work is the mechanism of rotor-stator rubbing [1]. Such interactions were investigated numerically using reduced-order finite element models [2,3,4] or simplified rotor models [5]. Except for reference [5] in which rotor-stator rub induced thermal expansion effects are studied, thermomechanical coupling is ignored in this context. However, thermomechanical coupling was observed to play a significant role in experimental measurements [6,3]. Indeed, heat is generated by friction and wear phenomena during rub events, leading to material expansion, diffusion and convection which induce strains and affect the contact configuration. Accordingly, the present work targets a reduction method capable of coping with thermomechanical models. Only a few investigations have been dedicated to reduction methods to for coupled problems combining Newton's second law of motion and the heat equation. They were mostly developed for spacecraft applications [7], electronic packages [8], or microelectromechanical systems [9].
Concerning turbomachinery applications, a reduced-order model was developed for the computation of thermal stresses based on a modal synthesis using thermal eigenmodes [10], taking into account convective exchange but ignoring structural dynamics. In the linear framework, a reduction technique for thermomechanical MEMS models has been developed, relying on the Craig-Bampton (CB) method [11] and describing both structural and thermal dynamics. This method was implemented to simulate rotor-stator thermomechanical interactions, using an explicit time-stepping algorithm [12].While it performed well on the small model described herein, some efficiency issues were encountered for models of larger orders: the convergence rate of the method is too slow to use it on industrial models.
Taking advantage of the Krylov subspace [13,14] and Craig-Bampton/Craig-Hale methods [15,16], this paper proposes a new reduction method for coupled problems yielding significant improvement in terms of computational efficiency. As for the reduction process, the simulation of thermomechanical contact problems poses modeling and solving issues. The robust treatment of unilateral contact conditions inspired by numerical methods developed in nonsmooth dynamics is detailed in a companion paper [17]. The model of interest is first introduced (section 1), then the reduction methodology is described (section 2). The results are illustrated and compared to other available methods by means of frequency response functions and time-stepping simulations (section 3).

ROTOR-STATOR MODEL
This section introduces the simplified rotor-stator model used throughout this study, with a flexible rotor and a rigid stator. It is supplemented with a more realistic model in section 4. experimental setups involving partial bladed-disks [18] or over-1 sized blades with contact occurrences on a single or few blades only [3]. The model embeds distinctive features of a turbomachinery compressor bladed-disk sector: a constrained internal bore and a contact interface along the blade tip. The finite element mesh is created using SOLID226 coupled-field elements in ANSYS R and possesses 507 free nodes (i.e. 2028 DOFs for the three directions in space and the temperature), seven of which are boundary nodes (blue dots in fig. 1). A standard TA6V titanium alloy, widely used in aircraft compressor components, is considered with the following estimated properties: mass density 4430 kg m 3 , Young's modulus 110 GPa, Poisson's ratio 0:3, heat capacity 520 J K 1 kg 1 , heat conductivity 6:7 J K 1 m 1 and dilatation coefficient 9 µm m 1 K 1 .

Linear thermoelastic model
A linear thermoelastic framework is adopted. In particular, centrifugal stiffening and temperature-dependence of material properties are excluded, as well as convection. The space semi-discretized governing equations are where M uu , C uu and K uu , are the mass, mechanical damping and stiffness matrices respectively. Matrices C Â Â and K Â Â reflect heat capacity and heat conductivity; K uÂ accounts for thermomechanical coupling effects due to thermal expansion. u stands for the relative nodal displacements, and Â the for the nodal temperatures. In order to activate contact, the external force vector f includes a sinusoidal term of amplitude 100 N in the radial direction, with a frequency of 33 Hz corresponding to 2000 rotations per minute, on each of the seven boundary nodes. It also stores normal forces and tangential friction. The external heat fluxes are imposed through q.

Thermomechanical contact
Numerically, rotor-stator rubbing events have essentially been dealt with using the penalty [19] and the forward Lagrange multiplier [20,1] methods. The penalty approach is known to induce numerical stiffness [21,22,1]. The forward Lagrange multiplier method relies on an explicit scheme which suffers from stability limitations, especially coping with a thermomechanical model [17]. However, a robust numerical method dedicated to nonsmooth problems, namely, the Moreau-Jean scheme [22], is employed in the present work. Its description is detailed in a companion paper [17]. Contact conditions are incorporated via a signed gap function g.u/. Contact occurs whenever the gap is closed, that is, g.u/ D 0. The associated reaction impulsive force in the inward normal direction is such that the so-called Signorini conditions hold: For each of the seven contacting nodes, friction is imposed in the tangential direction through a friction coefficient D 0:15 which multiplies . Frictional heating proportional to the normal reaction is included in q, with a coefficient of 0:1 m s 1 [23]. There are hence two thermal couplings: one linear through K uÂ and one nonlinear via q.

MODEL ORDER REDUCTION
The idea of reduction methods is to approximate the solution of the governing equations (1) with the solutions of another system of Ordinary Differential Equations with fewer DOFs. With the appropriate notations, eq. (1) simplifies to MR x C CP x C Kx D f. The reduction consists in assuming x RQ x, then left-multiplying the differential equation with a matrix P > such that: The choice of the reduction method is governed by the definition of R and P [24,25]. The proposed strategy consists in reducing mechanics and thermics independently, that is R and P are blockdiagonal matrices, so that the reduced equations keep the same structure as eq. (1) [26]. This is an inappropriate assumption when coupling terms are considered in R [11]. Additionally, we choose P D R, as in the Craig-Bampton (CB) reduction.
In the following, we detail two families of reduction methods, which will be combined: component mode synthesis (specifically CB) and block-Krylov (specifically CH). CB accurately captures the mechanical behaviour while CH describes better transient temperature evolutions. For the sake of conciseness, only the thermal part of eq. (1) is considered in this section, that is The inefficiency of CB to approach the solution to this heat equation will now be illustrated.

Craig-Bampton order-reduction method for conductive heat transfer: derivation and deficiencies
Component mode synthesis methods, such as MacNeal [27] or Craig-Bampton, consist in approximating the response within a subspace spanned by a selection of linear eigenmodes. CB is known for its ability to preserve stability properties and for its numerical robustness [28]. It is based on the definition of interior and boundary DOFs. The interior DOFs are approximated by chosen mode shapes R while boundary DOFs are preserved in the substitution x RQ x. This distinction is relevant when dealing with contact problems as the contact forces directly apply to DOFs of the reduced model [15,4]. The splitting between internal (subscripts i) and boundary (subscripts b) DOFs is mirrored in the reorganization of eq. (4) in the form where the superscript Â Â has been dropped in C and K to simplify the notations. The fixed-boundary eigenmodes are computed by solving the eigenproblem . C ii C K ii / D 0, and the m eigenmodes corresponding to the lowest eigenvalues are stored inˆi m . The CB reduction matrix also includes Guyan static constraint modes [29], which capture the static response of the complete system for prescribed boundary loads. More precisely, the j th Guyan mode corresponds to the response of the internal DOFs to a prescribed unit temperature on the j th boundary node and zero temperature on the other boundary nodes. Together with the abovementioned internal eigenmodes, CB method leads to a reduction matrix of the form Since boundary DOFs are not reduced, a larger number of boundary DOFs implies more CPU-intensive computations.
Despite its successful applications in thermal and or thermomechanical simulations involving slow excitations with respect to the characteristic time constants of the system [10,11], CB method provides poor results with the present application. The model presented in section 1 is a piecewise-linear thermoelastic model subject to localized and high-frequency mechanical and thermal loads. Indeed, frictional heating is proportional to the contact pressure, which exhibits fast variations (typically, a few milliseconds), while heat diffusion is a much slower mechanism (of the order of a second). The main challenge is therefore to build a compact Reduced-Order Model (ROM) that covers a broad range of excitation frequencies.
The gain of the thermal transfer function 1 of boundary node 1 is depicted in fig. 2, for the full model illustrated in fig. 1 and successive CB ROMs. Even with m Â D 400, that is with 80 % of the DOFs of the full model, the error for 1000 Hz reaches about 75 %. While not problematic for low-frequency applications, this inaccuracy leads to significant errors in time-integration with contact, see section 3. To summarize, CB method fails to describe response to high-frequency excitations efficiently.

Krylov subspace iteration techniques
For a given non-singular matrix A 2 R n n and a vector b 2 R n , a Krylov subspace of a n-dimensional vector space is a`-dimensional subspace spanned by the`vectors .b; Ab; A 2 b; : : : ; A` 1 b/, with Ä n. Krylov model reduction methods are based on the conservation of the transfer function of the system throughout the reduction process [14]. Assuming K Â Â is non-singular, eq. (4) implies with A D K ÂÂ 1 C Â Â . In the Laplace domain and denoting by s the Laplace variable, eq. (7) becomes .I As/Â D b. The transfer function H.s/ is defined by H.s/ D .I As/ 1 . Its Taylor expansion around s D 0 is The thermal transfer function was computed by means of a Fourier transform of the governing ODE (4) with a unit heat flux. so that The vectors m k are called the moments of the transfer function and the family .m 1 ; : : : ; m`/ forms the Krylov subspace associated to matrix A and vector b. It can be efficiently computed using Arnoldi [30] or Lanczos [31] algorithms or their recent improvements [14]. Krylov reduction methods use the matrix R D OEm 1 m 2 : : : m`. They can be extended to second-order systems [32]. The main limitation with these methods for the present application is that they do not preserve boundary nodes in the reduced basis. Dealing with contact conditions hence necessitates a systematic mapping between the full and reduced bases and thus requires additional computations.
The Taylor series of H.s/ can also be expanded around arbitrary points s and for various loads b. Krylov methods involving several expansion points are referred to as Rational Krylov methods [33]. They improve the precision of the reduced models over broader frequency ranges. Methods using multiple vectors b (stored in a matrix B) for the computation of the reduction basis, are called Block-Krylov Methods [33]. For each order of expansion, the moments are computed for B instead of b, resulting in matrix blocks instead of vectors: R D OEM 1 M 2 : : : M`. The block-Krylov subspace is then the span of the columns in .B; AB; A 2 B; : : : ; A` 1 B/.

Proposed method: Rational Craig-Hale
The reduction methodology proposed in this work combines the efficiency of boundary/interior splitting for contact with the broadfrequency capabilities of Krylov-subspace methods. The sought reduction matrix is thus of the form (6) and the objective is to find an appropriateˆi m using Krylov subspaces.
Moment-matching techniques were combined with boundary/interior splitting in the Craig-Hale (CH) method [16] and affiliates [25,34]. They differ from CB in that they use a block-Krylov family of moments inˆi m instead of component normal modes for the description of the internal dynamics of the system. In CH,ˆi m is built through a single expansion around s D 0, resulting in a limited frequency range of accuracy. Instead, the extension of CH to multiple expansion points s is implemented in the present work. Additional vectors will be included inˆi m . By analogy to Rational Krylov Methods, we refer to the following method as the Rational Craig-Hale method.
The Laplace transform of (5) is arranged as Unit temperatures are then successively enforced on each boundary DOF (i.e. Â b D OE1 0 0 : : : 0 > , Â b D OE0 1 0 : : : 0 > , etc.) and the internal equilibrium q i D 0 is guaranteed by an appropriate choice of the interior temperature Â i . Stacking the successive OEÂ b ; Â i > column-wise generates a matrix of the form OEI bbˆib > such that for an appropriate Q bb . The second block row of eq. (11) yieldŝ ib .s/ D .C ii s C K ii / 1 .C ib s C K ib /: As such, OEI bbˆib .s/ > is written as a static constraint mode matrix (see first column of right-hand side in eq. (6)). In order to compute the fixed interface modes fromˆi b .s/, the static constraint modes participation should be deducted, from the linear transformation leading to a structure of the usual CB reduction matrix (right hand side in eq. (6)).
Because it provides the temperature of the interior nodes in response to a loading on the boundary,ˆi b .s/ can be seen as a transfer function. It is now expanded using multiple points. Let s 0 be one of them. The Taylor expansion ofˆi b .s/ in eq. (12) takes the general form where In this formulation, H 0 ib , H 1 ib and H 2 ii correspond to zero-, first-and second-order terms, respectively. As in eq. (9), the block mode M k is given by the matrix coefficient of .s 0 s/ k , which depends on In the following, all expansion points are chosen as real numbers as proposed in [13, §6.2.2], corresponding to frequencies, providing a real reduction matrix and thus faster computations, as well as lower memory requirements.

Coupled thermomechanical reduction
We now construct the reduction of the full model (1) by combining the previous developments with the standard CB reduction for the mechanical system. Rearranging OEu Â > into OEu b Â b u i Â i > leads to the reduction basis R D ib are the static constraint modes. The variables Q u m ; Q Â m of the reduced-order model are defined as follows: This substitution is introduced in the governing equations which are then left-multiplied by R > , as in eq. (3), yielding the reduced system of order 4b C m u C m Â with m Â D bp.`C 1/.

SIMPLIFIED MODEL RESULTS
The CB method is efficient for reducing mechanics with unilateral constraints, but fails in thermal analyses due to the fast thermal excitations stemming from frictional heating, as explained in section 2.1. To overcome this limitation, the Rational Craig-Hale (RCH) method was proposed in section 2.4 to accurately reduce the dimension of the thermal problem. A validation analysis is now conducted through comparison with other existing methods. Unless stated otherwise, results are presented for a time step h D 10 5 s as in the companion paper [17] and m u D 10.

Validation of the sole thermal model reduction
We first consider a thermal model alone, with no mechanics and no contact, to validate the ability of the reduced model to capture rapid variations of heat fluxes. Since all computed transfer functions were found to follow the same trends, results are only shown for the boundary node 1 (see fig. 1). The input is chosen as a unitary thermal load on the boundary nodes.
The influence of the number of expansion points is depicted in terms of frequency response function in fig. 3. The expansion points are logarithmically spaced in the expected frequency range OE10 2 Hz; 10 5 Hz of the thermomechanical model.Each expansion is performed at order 0, that is`D 0 in eq. (14). As the test model possesses b D 7 boundary nodes, the size of the RCH models is 7.p C 1/. Figure 3 shows that increasing the amount of expansion points leads to a global and fast reduction of the transfer function error between RCH and full models. With p D 7, the RCH model comprises 56 thermal DOFs and its frequency response functions displays a peak error of 1:4 %. Similarly, fig. 4 depicts the convergence for one expansion point (p D 1) and multiple expansion orders`. As`increases, the error decreases quickly. In both cases, the convergence of the reduction method is fast, with maximum errors of less than 1 % with less than 100 DOFs. It is also possible to select multiple expansion points with multiple orders. For example, the peak error with p D`D 3 is 0:4 %, for a RCH model of 70 DOFs-to be compared to the 2028 DOFs of the full model. Though it might not hold true for other models, the convergence is faster when increasing the amount of expansion points rather than the truncation order. Accordingly, we select`D 0 in the following. Optimizing the frequencies of the expansion points in order to minimize the error is not discussed here but could increase the efficiency of RCH reduction even further. Figure 5 offers a comparison of the frequency response functions between CB, CH and three RCH models, all with m D 70. CB and CH models have been computed using the conventional Craig-Bampton and Craig-Hale methods. The three RCH models have been created using respectively a single Gram-Schmidt orthogonalization (sGS), the left-singular vectors obtained from a SVD performed on the completeˆi m matrix of eq. (15) (SVD), or multiple Gram-Schmidt orthogonalizations (mGS)-one every time a new vector is appended toˆi m . As observed before, the CB model is very inaccurate for moderate to high frequencies. Although better, the CH model displays smaller but still large errors at high frequencies, despite a ninthorder truncation. The RCH-mGS model is much better but still shows significant error ( 5 %) at high frequencies. In contrast, both RCH-sGS and RCH-SVD feature negligible errors, with a peak at 0:06 %. This shows that the RCH method works best with a single orthogonalization onˆi m or a SVD. The RCH-sGS is used in the remaining.

Thermomechanical contact simulations
Attention is now paid to the performance of the thermomechanical reduced models in practical use, that is simulations of time evolutions. Again, the numerical method employed to enforce the contact conditions are detailed in [17]. Results are provided for boundary node 3. Figure 6 depicts the impact of thermoelastic coupling on the radial displacement u r , temperature Â and contact force N . The only difference between the two simulations is that K uÂ D 0 when the formulation is uncoupled.
The model is clearly able to capture thermomechanical coupling and contact occurrences. It is worth noting that for the chosen node, coupling tends to increase the normal force because of local heat expansion effects. In practice, abradable coatings dissipate energy during contacts and limit the growth of the contact forces through material removal. During the time interval OE0:9 sI 1 s shown in fig. 7, the contact force begins to decrease, leading to a temperature drop. This load transfer from the considered node to another one may be a specificity of the studied configuration. With multiple nodes in contact, some may show more or less expansion under thermal stresses depending on their position on the blade tip and the behaviour could even be unstable.
The purpose of the present work is to demonstrate the relevance of the reduction and simulation methodologies, which should   be able to address larger-scale models. Here, 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 discussion. A more realistic model is presented in section 4. Figure 8 offers a comparison between the full coupled model and two RCH reduced-order models (p D 3 and p D 5). For p D 5, the curves cannot be distinguished. The RCH model calculated with 3 expansion points shows large errors 2 on the nodal displacement (11 %), temperature (28 %) and contact force (30 %). Since the temperature field is wrongly approximated, the temperature induced displacements are badly predicted, resulting in inaccurate contact forces. Due to the two-way coupling, the displacement error tends to grow during the simulation.
This can also be seen in the close-up view fig. 9: with p D 3, the contact phases are shorter, the temperature and the normal force 2 The error for x is computed as  The previous results have shown that the RCH model computed using reduction parameters (`D 0, p D 5) displays a relative error of 0:84 %, 1:86 % and 1:46 % in terms of nodal displacements, temperature and contact force. Similar errors are observed with the other nodes. Figure 11 compares the full model, the RCH ROM with`D 0 and p D 5 and a CB ROM comprising the same amount of DOFs (m Â D p b D 35 and m u D 10). The radial displacements and contact forces are all similar. The CB ROM however displays significant errors in terms of nodal temperatures. The temperature plot shows that it fails to describe the entire temperature variations. This is consistent with the frequency response functions discussed in section 2.1. Interestingly, the error in temperature does not generate large errors in u r and N during the first second, presumably because the average of the temperature is reasonably estimated. However the errors slowly increases with time, which could become problematic for longer simulations. The relative errors are of 2:1 %, 19 %, and 5:6 % for the nodal displacement, temperature and contact force.

INDUSTRIAL APPLICATION
Now that the convergence of the reduced models has been shown on a low order model, we focus on the reduction of a more realistic model, shown in fig. 12. This model comprises 5686 nodes and a total of 18768 DOF (4692 thermal), including 9 boundary nodes located on the middle plane of the blade tip.
To assess the convergence of RCH ROMs on this model, transfer functions were again computed. CB ROMs were also computed to serve as references. Figures 13 and 14  sponds to b D 9 boundary nodes and the following number p of expansion points of order`, in the form .p;`/: .3; 1/, .5; 1/, .3; 2/, .7; 1/, .9; 1/, .5; 2/, .7; 2/. As with the simplified blade model, errors with CB ROMs are large even with high order models. In contrast, the worst RCH model shows a maximum error of 18 % for a very small model order (m Â D 27). For greater orders, the accuracy of RCH-reduced models is such that their transfer functions cannot be distinguished from the RCH-full models ( fig. 14).
To verify all previous observations in the case of contact events, an actual contact simulation was performed on the full model and on one RCH and one CB ROMs of equal order (m u D 10, m Â D 126). During this simulation, the blade is excited on its first bending mode frequency (507 Hz) and modeshape. Simulation  results for node 1 are gathered in fig. 15. The RCH ROM is shown to be very accurate in terms of displacement, temperature and contact force, during the studied time range. The CB ROM shows a significant error in the temperature forecast: too low at the beginning and too high at the end of the simulation. Even though the contact forces do not show large differences with the full order computation, there is a growing error in the displacement amplitude. Such behaviour could lead to large response differences for longer simulations.
Timings for the computation of the thermal dynamic reduction matrixˆÂ Â im are indicated in fig. 16. All computations were  performed on identical hardware using MATLAB R software. The computational effort required to form the reduction basis is of the same order for both CB and RCH reduction method, and negligible compared to the simulation time of an industrial model.
Computation times required to perform the contact simulation shown in fig. 15 are given in table 1. The full order model computation was performed on 18 processors with clock frequency of 2:8 GHz, while the reduced model simulations were performed on 4 processors of 3:3 GHz clock frequency. The advantages of the proposed RCH ROM is clear: it is both fast and accurate.

CONCLUSION
A reduction method for thermomechanical finite element models with contact conditions was proposed. It advantageously combines  the Craig-Bampton technique, which preserves contacting nodes, and Krylov modes, which can be tuned to cover the appropriate range of frequencies. The mechanical problem is solved using Craig-Bampton reduction, while the thermal part is reduced with a modification of Craig-Hale method relying on the expansion of the frequency response functions of the system around multiple frequencies. The choice of these frequencies offers a large flexibility and can be adjusted to increase the accuracy in some targeted frequency ranges. In this work, the expansion points were equally spaced in a log scale. It was shown that the proposed method is much more accurate than existing reduction techniques for an equivalent reduced order, both on a simple model and on an industrial model with nearly 20000 DOFs. Further developments include optimizing the choice of the expansion frequencies for a more accurate approximation. Full scale thermomechanical rotor-stator interaction simulations, with a full rotor and stator model, are planned in a near future. Considering more realistic thermal boundary conditions and nonlinear thermal behavior would also be relevant. Fortunately, Krylov methods are still well suited since any convection load or imposed temperature constraint of known shape can be described easily with very few vectors. Moreover, nonlinear models can also be reduced using Krylov methods [35].