Nicolas Guérin
email: nicolas.guerin@ec-lyon.fr

Anders Thorin

Fabrice Thouverez

Mathias Legrand

Patricio Almeida

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

Keywords: component mode synthesis, model order reduction, thermomechanical coupling, impeller, rotor-stator interaction, Krylov subspace method

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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 [START_REF] Jacquet-Richardet | Rotor to stator contacts in turbomachines. Review and application[END_REF]. Such interactions were investigated numerically using reduced-order finite element models [START_REF] Legrand | Numerical investigation of abradable coating removal in aircraft engines through plastic constitutive law[END_REF][START_REF] Batailly | Experimental and numerical simulation of a rotor/stator interaction event localized on a single blade within an industrial high-pressure compressor[END_REF][START_REF] Batailly | Assessment of reduced models for the detection of modal interaction through rotor stator contacts[END_REF] or simplified rotor models [START_REF] Goldman | Rotor-to-stator, rub-related, thermal/mechanical effects in rotating machinery[END_REF]. Except for reference [START_REF] Goldman | Rotor-to-stator, rub-related, thermal/mechanical effects in rotating machinery[END_REF] 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 [START_REF] Almeida | Experimental analysis of dynamic interaction between a centrifugal compressor and its casing[END_REF][START_REF] Batailly | Experimental and numerical simulation of a rotor/stator interaction event localized on a single blade within an industrial high-pressure compressor[END_REF]. 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 [START_REF] Murthy | A construction of thermal basis functions for coupled structural-thermal reduced order models[END_REF], electronic packages [START_REF] Zukowski | Model reduction for thermo-mechanical simulation of packages[END_REF], or microelectromechanical systems [START_REF] Bechtold | Krylov-subspace-based order reduction methods applied to generate compact-electro-thermal models for mems[END_REF].

Concerning turbomachinery applications, a reduced-order model was developed for the computation of thermal stresses based on a modal synthesis using thermal eigenmodes [START_REF] Botto | Reduced-order models for the calculation of thermal transients of heat conduction/convection fe models[END_REF], 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 [START_REF] Nachtergaele | Efficient weakly coupled projection basis for the reduction of thermomechanical models[END_REF] and describing both structural and thermal dynamics. This method was implemented to simulate rotor-stator thermomechanical interactions, using an explicit time-stepping algorithm [START_REF] Guérin | Thermomechanical component mode synthesis for blade cas-ing interaction prediction[END_REF].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 [START_REF] Grimme | Krylov projection methods for model reduction[END_REF][START_REF] Freund | Krylov-subspace methods for reduced-order modeling in circuit simulation[END_REF] and Craig-Bampton/Craig-Hale methods [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF][START_REF] Craig | Block-Krylov component synthesis method for structural model reduction[END_REF], 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 [START_REF] Thorin | Nonsmooth thermoelastic simulations of blade-casing contact interactions[END_REF]. 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.

System of interest

A single finite element sector of a simplified rotor is considered, as depicted in fig. 1. The stator is assumed rigid in line with rotor stator experimental setups involving partial bladed-disks [START_REF] Padova | Experimental results from controlled blade tip/shroud rubs at engine speed[END_REF] or over-sized blades with contact occurrences on a single or few blades only [START_REF] Batailly | Experimental and numerical simulation of a rotor/stator interaction event localized on a single blade within an industrial high-pressure compressor[END_REF]. 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

Ä M uu 0 0 0  R u R  à C Ä C uu 0 0 C    P u P  à C Ä K uu K u 0 K    u  à D  f q à (1)
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 [START_REF] Parent | Phenomenological model for stability analysis of bladed rotor-to-stator contacts[END_REF] and the forward Lagrange multiplier [START_REF] Legrand | Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction[END_REF][START_REF] Jacquet-Richardet | Rotor to stator contacts in turbomachines. Review and application[END_REF] methods. The penalty approach is known to induce numerical stiffness [START_REF] Stewart | Rigid-body dynamics with friction and impact[END_REF][START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF][START_REF] Jacquet-Richardet | Rotor to stator contacts in turbomachines. Review and application[END_REF]. The forward Lagrange multiplier method relies on an explicit scheme which suffers from stability limitations, especially coping with a thermomechanical model [START_REF] Thorin | Nonsmooth thermoelastic simulations of blade-casing contact interactions[END_REF]. However, a robust numerical method dedicated to nonsmooth problems, namely, the Moreau-Jean scheme [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF], is employed in the present work. Its description is detailed in a companion paper [START_REF] Thorin | Nonsmooth thermoelastic simulations of blade-casing contact interactions[END_REF]. 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:

8t 2 R C ; g.u.t
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 [START_REF] Ireman | Finite element algorithms for thermoelastic wear problems[END_REF].

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:

.P > MR/ R Q x C .P > CR/ P Q x C .P > KR/Q x D P > f: (3)
The choice of the reduction method is governed by the definition of R and P [START_REF] Besselink | A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control[END_REF][START_REF] Qu | Model Order Reduction Techniques[END_REF]. 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) [START_REF] Schilders | A novel approach to model order reduction for coupled multiphysics problems[END_REF]. This is an inappropriate assumption when coupling terms are considered in R [START_REF] Nachtergaele | Efficient weakly coupled projection basis for the reduction of thermomechanical models[END_REF]. 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

C Â Â P Â C K Â Â Â D q: (4)
The inefficiency of CB to approach the solution to this heat equation will now be illustrated.

2.1 Craig-Bampton order-reduction method for conductive heat transfer: derivation and deficiencies Component mode synthesis methods, such as MacNeal [START_REF] Macneal | A hybrid method of component mode synthesis[END_REF] 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 [START_REF] Bladh | Numerical instability of classical free-interface component mode synthesis techniques[END_REF]. 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 [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF][START_REF] Batailly | Assessment of reduced models for the detection of modal interaction through rotor stator contacts[END_REF]. The splitting between internal (subscripts i) and boundary (subscripts b) DOFs is mirrored in the reorganization of eq. ( 4) in the form

Ä C bb C bi C ib C ii  P  b P  i à C Ä K bb K bi K ib K ii   b  i à D  q b q i à (5)
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 ˆim .

The CB reduction matrix also includes Guyan static constraint modes [START_REF] Guyan | Reduction of stiffness and mass matrices[END_REF], 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

R D Ä ˆstatic ˆdynamic D Ä I bb 0 bm K 1 ii K ib ˆim : (6)
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 [START_REF] Botto | Reduced-order models for the calculation of thermal transients of heat conduction/convection fe models[END_REF][START_REF] Nachtergaele | Efficient weakly coupled projection basis for the reduction of thermomechanical models[END_REF], 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 [START_REF] Freund | Krylov-subspace methods for reduced-order modeling in circuit simulation[END_REF]. Assuming K Â Â is non-singular, eq. ( 4) implies

 A P  D b (7)
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

H.s/ D ` 1 X kD0 A k s k C o.s ` 1 / (8)
1 The thermal transfer function was computed by means of a Fourier transform of the governing ODE (4) with a unit heat flux. so that

Â.s/ D H.s/b D ` 1 X kD0 m kC1 s k C o.s ` 1 / (9)
with m kC1 D A k b. 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 [START_REF] Arnoldi | The principle of minimized iterations in the solution of the matrix eigenvalue problem[END_REF] or Lanczos [START_REF] Lanczos | An iteration method for the solution of the eigenvalue problem of linear differential and integral operators[END_REF] algorithms or their recent improvements [START_REF] Freund | Krylov-subspace methods for reduced-order modeling in circuit simulation[END_REF]. Krylov reduction methods use the matrix R D OEm 1 m 2 : : : m `. They can be extended to second-order systems [START_REF] Bai | Dimension reduction of large-scale second-order dynamical systems via a second-order Arnoldi method[END_REF]. 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 [START_REF] Elfadel | A block rational Arnoldi algorithm for multipoint passive model-order reduction of multiport RLC networks[END_REF]. 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 [START_REF] Elfadel | A block rational Arnoldi algorithm for multipoint passive model-order reduction of multiport RLC networks[END_REF]. 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 ˆim using Krylov subspaces.

Moment-matching techniques were combined with boundary/interior splitting in the Craig-Hale (CH) method [START_REF] Craig | Block-Krylov component synthesis method for structural model reduction[END_REF] and affiliates [START_REF] Qu | Model Order Reduction Techniques[END_REF][START_REF] Pohle | Reduced order model of mistuned bladed disks using the Krylov sub-space combined with the Craig-Bampton reduction[END_REF]. They differ from CB in that they use a block-Krylov family of moments in ˆim instead of component normal modes for the description of the internal dynamics of the system. In CH, ˆim 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 ˆim . 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

Ä C bb C bi C ib C ii s C Ä K bb K bi K ib K ii !   b  i à D  q b q i à : (10)
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

Ä C bb C bi C ib C ii s C Ä K bb K bi K ib K ii ! Ä I bb ˆib D Ä Q bb 0 ib ( 11
)
for an appropriate Q bb . The second block row of eq. ( 11) yields

ˆib .s/ D .C ii s C K ii / 1 .C ib s C K ib /: (12)
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 ˆib .s/, the static constraint modes participation should be deducted, from the linear transformation

Ä I bb I bb K 1 ii K ib ˆib .s/ ! Ä I bb 0 bb K 1 ii K ib ˆib .s/ C K 1 ii K ib ( 13
)
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, ˆib .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 ˆib .s/ in eq. ( 12) takes the general form

ˆib .s/ D H 0 ib X kD1 .H 2 ii / k 1 H ìb .s 0 s/ k C o..s 0 s/ `/ (14)
where

H 0 ib D .A 0 ii / 1 .C ib s 0 C K ib /, H 1 ib D .A 0 ii / 1 .C ii H 0 ib C ib /, H 2 ii D .A 0 ii / 1 C ii and A 0 ii D .C ii s 0 C K ii /.
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 H 0 ib , H 1 ib and H 2 ii . Such a formulation is valid for any order `and any point of expansion s 0 . For p expansion points .s .1/ 0 ; : : : ; s .p/ 0 / and a common order expansion `for simplicity, p `block modes are obtained: .M .1/ 1 ; : : : ; M .1/ `/; : : : ; .M .p/ 1 ; : : : ; M .p/ `/. Applying transformation [START_REF] Grimme | Krylov projection methods for model reduction[END_REF] to each of the M .j / 1 and defining Q

M .1/ 1 D M .1/ 1 C K 1 ii K ib yields ˆim D Ä h Q M .
1/ 1 M .1/ 2 : : : M .1/ `i : : :

h Q M .p/ 1 M .p/ 2 : : : M .p/ `i (15)
which, inserted in the CB form [START_REF] Almeida | Experimental analysis of dynamic interaction between a centrifugal compressor and its casing[END_REF], provides the reduction matrix R of the Rational Craig-Hale reduction method. The assembly of the reduction basis ˆim is summarized in Algorithm 1. The dimension of ˆim is i m where m D bp.`C 1/ is the product of the number of boundary degrees of freedom b, the order of expansion plus one `C 1 and the number of expansion points p. For numerical reasons, it is recommended to perform either a Gram-Schmidt orthogonalization on the columns of ˆim or the use of the leftsingular vectors stemming from a SVD of ˆim (see next section).

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 where ˆuu im stacks the chosen internal modes of vibration (eigenvectors of M uu ii 1 K uu ii ), ˆÂ im is obtained from eq. ( 15),

0 ii .C ii s 0 .n/ C K ii / H 0 ib .A 0 ii / 1 .C ib s 0 .n/ C K ib / M .n/ 1 H 0 ib C K 1 ii K ib if ` 1 then H 1 ib .A 0 ii / 1 .C ii H 0 ib C ib / M .n/ 2 H ìb if ` 2 then H 2 ii .A 0 ii / 1 C ii for k 2 to `do M .n/ k .H 2 ii / k
‰ uu ib D K uu ii 1 K uu ib and ‰ Â Â ib D K Â Â ii 1 K Â Â
ib are the static constraint modes. The variables Q u m ; Q Â m of the reduced-order model are defined as follows:

2 6 6 4 u b  b u i  i 3 7 7 5 R 2 6 6 4 u b  b Q u m Q  m 3 7 7 5 : (17)
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 [START_REF] Thorin | Nonsmooth thermoelastic simulations of blade-casing contact interactions[END_REF] 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 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 ˆim matrix of eq. ( 15) (SVD), or multiple Gram-Schmidt orthogonalizations (mGS)-one every time a new vector is appended to ˆim . 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 ˆim 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 [START_REF] Thorin | Nonsmooth thermoelastic simulations of blade-casing contact interactions[END_REF]. 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 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. total of 18768 DOF (4692 thermal), including 9 boundary nodes located on the middle plane of the blade tip.

R 1 0 jx ROM x FULL jdt R 1 0 jx FULL jdt .
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 and14 display the gain of node 1 and the maximal gain error measured over all boundary nodes, for the CB and RCH ROMs, respectively. The number of modes corre- 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.

3 1 FIGURE 1 :

 311 FIGURE 1: Finite element model of bladed-disk sector: boundary nodes [ ], constrained nodes [ ] and contact forces [ ].

1 ]FIGURE 2 :

 12 FIGURE 2: Thermal CB ROM transfer function convergence. Full model [ ], m  D 1 [ ], 10 [ ], 100 [ ] and 400 [ ]. For m  D 500, the error is of the magnitude of machine precision.

  Vector of expansion points s 0 ; Expansion order for n 1 to p do A

1 ]FIGURE 3 :

 13 FIGURE 3: RCH ROM transfer function convergence with order `D 0 and p expansion points. Full model [ ], p D 1 [ ], 3 [ ] and 7 [ ]. possesses b D 7 boundary nodes, the size of the RCH models is 7.p C 1/.Figure3shows 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.4depicts the convergence for one expansion point (p D 1) and multiple expansion orders `. As `increases, the error decreases quickly. In both cases,

1 ]FIGURE 4 :

 14 FIGURE 4: RCH ROM transfer function convergence with p D 3 expansion points and a `-th order expansion. Full model [ ], `D 0 [ ], 1 [ ], 2 [ ] and 3 [ ].

1 ]FIGURE 5 :

 15 FIGURE 5: Comparison of CB [ ], CH [ ], RCH-mGS [ ], RCH-SVD [ ] and RCH-sGS [ ] ROMs transfer functions with 70 modes; full model [ ].

FIGURE 6 :

 6 FIGURE 6: Thermomechanical model responses at node 3. Uncoupled [ ] and coupled [ ] models.

FIGURE 7 :

 7 FIGURE 7: Close-up of fig. 6.

FIGURE 8 :

 8 FIGURE 8: RCH ROM convergence in contact simulation. Full model [ ], RCH ROM with `D 0, p D 3 [ ] and 5 [ ].

FIGURE 9 :

 9 FIGURE 9: Close-up view of fig. 8. are underestimated.

  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).

FIGURE 10 :

 10 FIGURE 10: CB and RCH ROM comparison in contact simulation. Full model [ ], CB (m  D p b D 35) [ ] and RCH (`D 0, p D 5) [ ].

FIGURE 11 : 1 FIGURE 12 :

 11112 FIGURE 11: Close-up view of fig. 11.

FIGURE 13 :

 13 FIGURE 13: CB ROM transfer function convergence: full model [ ], m  D 27 [ ], 45 [ ], 54 [ ], 63 [ ], 81 [ ], 90 [ ], 126 [ ].

FIGURE 14 :

 14 FIGURE 14: RCH ROM transfer function convergence. Full model [ ], m  D 27 [ ], 45 [ ], 54 [ ], 63 [ ], 81 [ ], 90 [ ] and 126 [ ].

FIGURE 15 :

 15 FIGURE 15: Time histories for node 1. Full model [ ], RCH [ ], CB ROM [ ] with m u D 10 and m  D 126.

  [START_REF] Craig | Block-Krylov component synthesis method for structural model reduction[END_REF]. All computations were

	Coulomb friction coefficient Laplace variable Expansion point `Taylor series order of truncation s s 0 p Number of expansion points m u Number of mechanical dynamic modeshapes m  Number of thermal dynamic modeshapes CB Craig-Bampton CH Craig-Hale RCH Rational Craig-Hale ROM Reduced-Order Model DOF Degree Of Freedom SVD Singular Value Decomposition
	Computation times [s]	27 0 1 2	45	54	63	81	90	126
			Number of thermal modeshapes m Â

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 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 [START_REF] Cardona | Solution of non-linear thermal transient problems by a reduction method[END_REF].

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.

NOMENCLATURE