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

Nicolas Guérin, Anders Thorin, Fabrice Thouverez, Mathias Legrand, Patricio Almeida

## ► To cite this version:

Nicolas Guérin, Anders Thorin, Fabrice Thouverez, Mathias Legrand, Patricio Almeida. Thermomechanical model reduction for efficient simulations of rotor-stator contact interaction. ASME Turbo Expo 2018: Turbomachinery Technical Conference and Exposition, Jun 2018, Oslo, Norway. 10.1115/ GT2018-75880. hal-01713849

HAL Id: hal-01713849

https://hal.science/hal-01713849

Submitted on 1 Mar 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# THERMOMECHANICAL MODEL REDUCTION FOR EFFICIENT SIMULATIONS OF ROTOR-STATOR CONTACT INTERACTION

Nicolas Guérin(a,c), Anders Thorin(b), Fabrice Thouverez(a), Mathias Legrand(b), Patricio Almeida(c)

(a) École Centrale de Lyon, Laboratoire de Tribologie et Dynamique des Systèmes, Écully, France

(b) Structural Dynamics and Vibration Laboratory, McGill University, Montreal, Canada

(c) Safran Helicopter Engines, Bordes, France

Email: nicolas.guerin@ec-lyon.fr

## ABSTRACT

Turbomachinery rotor-stator unilateral contact induced interactions play a growing role in lifecycle analysis and thus motivate the use of accurate numerical prediction tools. Recent literature confirmed by ongoing in-house experiments have shown the importance of thermomechanical coupling effects in such interactions. However, most available (possibly reduced-order) models are restricted 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 specifically 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.

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

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

## 1 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.

### 1.1 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

<figure><img src="image_2.png" /><figcaption>**FIGURE 1:** Finite element model of bladed-disk sector: boundary nodes [●], constrained nodes [■] and contact forces [---].</figcaption></figure>

experimental setups involving partial bladed-disks [18] or over-
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® 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⁻³, Young's modulus 110 GPa, Poisson's ratio 0.3, heat capacity 520 J K⁻¹ kg⁻¹, heat conductivity 6.7 J K⁻¹ m⁻¹ and dilatation coefficient 9 µm m⁻¹ K⁻¹.

## 1.2 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

$$ \begin{bmatrix} \mathbf{M}^{uu} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} \begin{pmatrix} \ddot{\mathbf{u}} \\ \ddot{\theta} \end{pmatrix} + \begin{bmatrix} \mathbf{C}^{uu} & \mathbf{0} \\ \mathbf{0} & \mathbf{C}^{\theta\theta} \end{bmatrix} \begin{pmatrix} \dot{\mathbf{u}} \\ \dot{\theta} \end{pmatrix} + \begin{bmatrix} \mathbf{K}^{uu} & \mathbf{K}^{u\theta} \\ \mathbf{0} & \mathbf{K}^{\theta\theta} \end{bmatrix} \begin{pmatrix} \mathbf{u} \\ \theta \end{pmatrix} = \begin{pmatrix} \mathbf{f} \\ \mathbf{q} \end{pmatrix} \quad (1) $$

where $\mathbf{M}^{uu}$, $\mathbf{C}^{uu}$ and $\mathbf{K}^{uu}$, are the mass, mechanical damping and stiffness matrices respectively. Matrices $\mathbf{C}^{\theta\theta}$ and $\mathbf{K}^{\theta\theta}$ reflect heat capacity and heat conductivity; $\mathbf{K}^{u\theta}$ accounts for thermomechanical coupling effects due to thermal expansion. $\mathbf{u}$ stands for the relative nodal displacements, and $\theta$ the for the nodal temperatures. In order to activate contact, the external force vector $\mathbf{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 $\mathbf{q}$.

## 1.3 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) = 0$. The associated reaction impulsive force $\lambda$ in the inward normal direction is such that the so-called Signorini conditions hold:

$$ \forall t \in \mathbb{R}_+, \quad g(u(t)) \ge 0, \quad \lambda(t) \ge 0, \quad g(u(t))\lambda(t) = 0. \quad (2) $$

For each of the seven contacting nodes, friction is imposed in the tangential direction through a friction coefficient $\mu = 0.15$ which multiplies $\lambda$. Frictional heating proportional to the normal reaction $\lambda$ is included in $\mathbf{q}$, with a coefficient of 0.1 m s⁻¹ [23]. There are hence two thermal couplings: one linear through $\mathbf{K}^{u\theta}$ and one nonlinear via $\mathbf{q}$.

# 2 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 $\mathbf{M}\ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} + \mathbf{K}\mathbf{x} = \mathbf{f}$. The reduction consists in assuming $\mathbf{x} \approx \mathbf{R}\bar{\mathbf{x}}$, then left-multiplying the differential equation with a matrix $\mathbf{P}^\top$ such that:

$$ (\mathbf{P}^\top \mathbf{M} \mathbf{R}) \ddot{\mathbf{x}} + (\mathbf{P}^\top \mathbf{C} \mathbf{R}) \dot{\mathbf{x}} + (\mathbf{P}^\top \mathbf{K} \mathbf{R}) \bar{\mathbf{x}} = \mathbf{P}^\top \mathbf{f}. \quad (3) $$

The choice of the reduction method is governed by the definition of $\mathbf{R}$ and $\mathbf{P}$ [24, 25]. The proposed strategy consists in reducing mechanics and thermics independently, that is $\mathbf{R}$ and $\mathbf{P}$ are block-diagonal 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 $\mathbf{R}$ [11]. Additionally, we choose $\mathbf{P} = \mathbf{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

$$ \mathbf{C}^{\theta\theta} \dot{\theta} + \mathbf{K}^{\theta\theta} \theta = \mathbf{q}. \quad (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 [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 $\mathbf{R}$ while boundary DOFs are preserved in the substitution $\mathbf{x} \approx \mathbf{R}\bar{\mathbf{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

$$ \begin{bmatrix} \mathbf{C}_{bb} & \mathbf{C}_{bi} \\ \mathbf{C}_{ib} & \mathbf{C}_{ii} \end{bmatrix} \begin{pmatrix} \dot{\theta}_b \\ \dot{\theta}_i \end{pmatrix} + \begin{bmatrix} \mathbf{K}_{bb} & \mathbf{K}_{bi} \\ \mathbf{K}_{ib} & \mathbf{K}_{ii} \end{bmatrix} \begin{pmatrix} \theta_b \\ \theta_i \end{pmatrix} = \begin{pmatrix} \mathbf{q}_b \\ \mathbf{q}_i \end{pmatrix} \quad (5) $$

where the superscript $\theta\theta$ has been dropped in $\mathbf{C}$ and $\mathbf{K}$ to simplify the notations. The fixed-boundary eigenmodes are computed by solving the eigenproblem $(\lambda\mathbf{C}_{ii} + \mathbf{K}_{ii})\phi = \mathbf{0}$, and the $m$ eigenmodes corresponding to the lowest eigenvalues are stored in $\mathbf{\Phi}_{im}$.

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 above-mentioned internal eigenmodes, CB method leads to a reduction matrix of the form

$$ \mathbf{R} = \begin{bmatrix} \mathbf{\Phi}_{\text{static}} & \mathbf{\Phi}_{\text{dynamic}} \end{bmatrix} = \begin{bmatrix} \mathbf{I}_{bb} & \mathbf{0}_{bm} \\ -\mathbf{K}_{ii}^{-1} \mathbf{K}_{ib} & \mathbf{\Phi}_{im} \end{bmatrix}. \quad (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 thermo-mechanical 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<fnref n="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_θ = 400$, that is with 80% of

<figure><img src="image_3.png" /><figcaption>FIGURE 2: Thermal CB ROM transfer function convergence. Full model [—], $m^θ = 1$ [—], 10 [—], 100 [—] and 400 [—]. For $m^θ = 500$, the error is of the magnitude of machine precision.</figcaption></figure>

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.

## 2.2 Krylov subspace iteration techniques

For a given non-singular matrix $A \in \mathbb{R}^{n \times n}$ and a vector $b \in \mathbb{R}^n$, a Krylov subspace of a $n$-dimensional vector space is a $l$-dimensional subspace spanned by the $l$ vectors ($b, Ab, A^2b, \dots, A^{l-1}b$), with $l \le n$. Krylov model reduction methods are based on the conservation of the transfer function of the system throughout the reduction process [14]. Assuming $K^{\theta\theta}$ is non-singular, eq. (4) implies

$$ \theta - A\dot{\theta} = b \qquad (7) $$

with $A = -K^{\theta\theta^{-1}}C^{\theta\theta}$. In the Laplace domain and denoting by $s$ the Laplace variable, eq. (7) becomes $(I - As)\theta = b$. The transfer function $H(s)$ is defined by $H(s) = (I - As)^{-1}$. Its Taylor expansion around $s = 0$ is

$$ H(s) = \sum_{k=0}^{l-1} A^k s^k + o(s^{l-1}) \qquad (8) $$

so that

$$ \theta(s) = H(s)b = \sum_{k=0}^{l-1} m_{k+1} s^k + o(s^{l-1}) \qquad (9) $$

with $m_{k+1} = A^k b$. The vectors $m_k$ are called the moments of the transfer function and the family ($m_1, \dots, m_l$) 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 = [m_1 m_2 \dots m_l]$. 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 = [M_1 M_2 \dots M_l]$. The block-Krylov subspace is then the span of the columns in $(B, AB, A^2B, \dots, A^{l-1}B)$.

## 2.3 Proposed method: Rational Craig-Hale

The reduction methodology proposed in this work combines the efficiency of boundary/interior splitting for contact with the broad-frequency capabilities of Krylov-subspace methods. The sought reduction matrix is thus of the form (6) and the objective is to find an appropriate $\Phi_{im}$ using Krylov subspaces.

Moment-matching techniques were combined with boundary/interior splitting in the Craig-Hale (CH) method [16] and affixiates [25, 34]. They differ from CB in that they use a block-Krylov family of moments in $\Phi_{im}$ instead of component normal modes for the description of the internal dynamics of the system. In CH, $\Phi_{im}$ is built through a single expansion around $s = 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 $\Phi_{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

$$ \left( \begin{bmatrix} C_{bb} & C_{bi} \\ C_{ib} & C_{ii} \end{bmatrix} s + \begin{bmatrix} K_{bb} & K_{bi} \\ K_{ib} & K_{ii} \end{bmatrix} \right) \begin{pmatrix} \theta_b \\ \theta_i \end{pmatrix} = \begin{pmatrix} q_b \\ q_i \end{pmatrix}. \qquad (10) $$

Unit temperatures are then successively enforced on each boundary DOF (i.e. $\theta_b = [1 \ 0 \ 0 \ \dots \ 0]^{\top}$, $\theta_i = [0 \ 1 \ 0 \ \dots \ 0]^{\top}$, etc.) and the internal equilibrium $q_i = 0$ is guaranteed by an appropriate choice of the interior temperature $\theta_i$. Stacking the successive $[\theta_b, \theta_i]^{\top}$ column-wise generates a matrix of the form $[I_{bb} \ \Phi_{ib}]^{\top}$ such that

$$ \left( \begin{bmatrix} C_{bb} & C_{bi} \\ C_{ib} & C_{ii} \end{bmatrix} s + \begin{bmatrix} K_{bb} & K_{bi} \\ K_{ib} & K_{ii} \end{bmatrix} \right) \begin{bmatrix} I_{bb} \\ \Phi_{ib} \end{bmatrix} = \begin{bmatrix} Q_{bb} \\ 0_{ib} \end{bmatrix} \qquad (11) $$

<footnote n="1">¹The thermal transfer function was computed by means of a Fourier transform of the governing ODE (4) with a unit heat flux.</footnote>
for an appropriate $Q_{bb}$. The second block row of eq. (11) yields

$$ \Phi_{ib}(s) = -(C_{ii}s + K_{ii})^{-1}(C_{ib}s + K_{ib}). \quad (12) $$

As such, $[I_{bb} \Phi_{ib}(s)]^\top$ 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 $\Phi_{ib}(s)$, the static constraint modes participation should be deducted, from the linear transformation

$$ \left[ \begin{matrix} I_{bb} & I_{bb} \\ -K_{ii}^{-1} K_{ib} & \Phi_{ib}(s) \end{matrix} \right] \rightarrow \left[ \begin{matrix} I_{bb} & \mathbf{0}_{bb} \\ -K_{ii}^{-1} K_{ib} & \Phi_{ib}(s) + K_{ii}^{-1} K_{ib} \end{matrix} \right] \quad (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, $\Phi_{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 $\Phi_{ib}(s)$ in eq. (12) takes the general form

$$ \Phi_{ib}(s) = -H_{ib}^0 - \sum_{k=1}^{\ell} (H_{ii}^2)^{k-1} H_{ib}^{\ell} (s_0 - s)^k + o((s_0 - s)^{\ell}) \quad (14) $$

where $H_{ib}^0 = (A_{ii}^0)^{-1}(C_{ib}s_0 + K_{ib})$, $H_{ib}^1 = (A_{ii}^0)^{-1}(C_{ii}H_{ib}^0 - C_{ib})$, $H_{ii}^2 = (A_{ii}^0)^{-1}C_{ii}$ and $A_{ii}^0 = (C_{ii}s_0 + K_{ii})$. In this formulation, $H_{ib}^0$, $H_{ib}^1$ and $H_{ii}^2$ 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_{ib}^0$, $H_{ib}^1$ and $H_{ii}^2$.

Such a formulation is valid for any order $\ell$ and any point of expansion $s_0$. For $p$ expansion points $(s_0^{(1)}, \dots, s_0^{(p)})$ and a common order expansion $\ell$ for simplicity, $p \times \ell$ block modes are obtained: $(M_1^{(1)}, \dots, M_\ell^{(1)}), \dots, (M_1^{(p)}, \dots, M_\ell^{(p)})$. Applying transformation (13) to each of the $M_1^{(j)}$ and defining $\tilde{M}_1^{(1)} = M_1^{(1)} + K_{ii}^{-1} K_{ib}$ yields

$$ \Phi_{im} = \left[ \left[ \tilde{M}_1^{(1)} M_2^{(1)} \dots M_\ell^{(1)} \right] \dots \left[ \tilde{M}_1^{(p)} M_2^{(p)} \dots M_\ell^{(p)} \right] \right] \quad (15) $$

which, inserted in the CB form (6), provides the reduction matrix $R$ of the Rational Craig-Hale reduction method. The assembly of the reduction basis $\Phi_{im}$ is summarized in Algorithm 1. The dimension of $\Phi_{im}$ is $i \times m$ where $m = bp(\ell + 1)$ is the product of the number of boundary degrees of freedom $b$, the order of expansion plus one $\ell + 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 $\Phi_{im}$ or the use of the left-singular vectors stemming from a SVD of $\Phi_{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.

## 2.4 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 $[u \theta]^\top$ into $[u_b \theta_b u_i \theta_i]^\top$ leads to the reduction basis

$$ R = \begin{bmatrix} I_{bb}^{uu} & 0 & 0 & 0 \\ 0 & I_{bb}^{\theta\theta} & 0 & 0 \\ \Psi_{ib}^{uu} & 0 & \Phi_{im}^{uu} & 0 \\ 0 & \Psi_{ib}^{\theta\theta} & 0 & \Phi_{im}^{\theta\theta} \end{bmatrix} \quad (16) $$

**Data:** Vector of expansion points $s_0$; Expansion order $\ell$ for $n \leftarrow 1$ to $p$ do

$$ \begin{aligned} & A_{ii}^0 \leftarrow (C_{ii}s_0(n) + K_{ii}) \\ & H_{ib}^0 \leftarrow (A_{ii}^0)^{-1}(C_{ib}s_0(n) + K_{ib}) \end{aligned} $$

$$ \begin{aligned} & M_1^{(n)} \leftarrow -H_{ib}^0 + K_{ii}^{-1} K_{ib} \\ & \text{if } \ell \ge 1 \text{ then} \end{aligned} $$

$$ \begin{aligned} & H_{ib}^1 \leftarrow (A_{ii}^0)^{-1}(C_{ii}H_{ib}^0 - C_{ib}) \\ & M_2^{(n)} \leftarrow -H_{ib}^\ell \end{aligned} $$

$$ \begin{aligned} & \text{if } \ell \ge 2 \text{ then} \\ & \quad H_{ii}^2 \leftarrow (A_{ii}^0)^{-1}C_{ii} \\ & \quad \text{for } k \leftarrow 2 \text{ to } \ell \text{ do} \\ & \quad \quad M_k^{(n)} \leftarrow (H_{ii}^2)^{k-1}H_{ib}^1 \end{aligned} $$

end

end

$$ \Phi_{im} \leftarrow [M_1^{(1)} M_2^{(1)} \dots M_\ell^{(1)} \dots \tilde{M}_1^{(p)} M_2^{(p)} \dots M_\ell^{(p)}] $$

**Result:** Dynamic reduction basis $\Phi_{im}$

### ALGORITHM 1: Computation of the RCH basis $\Phi_{im}$.

where $\Phi_{im}^{uu}$ stacks the chosen internal modes of vibration (eigenvectors of $M_{ii}^{uu-1} K_{ii}^{uu}$), $\Phi_{im}^{\theta\theta}$ is obtained from eq. (15), $\Psi_{ib}^{uu} = -K_{ii}^{uu-1} K_{ib}^{uu}$ and $\Psi_{ib}^{\theta\theta} = -K_{ii}^{\theta\theta-1} K_{ib}^{\theta\theta}$ are the static constraint modes. The variables $\tilde{u}_m, \tilde{\theta}_m$ of the reduced-order model are defined as follows:

$$ \begin{bmatrix} \mathbf{u}_b \\ \theta_b \\ \mathbf{u}_i \\ \theta_i \end{bmatrix} \approx R \begin{bmatrix} \mathbf{u}_b \\ \theta_b \\ \tilde{\mathbf{u}}_m \\ \tilde{\theta}_m \end{bmatrix}. \quad (17) $$

This substitution is introduced in the governing equations which are then left-multiplied by $R^\top$, as in eq. (3), yielding the reduced system of order $4b + m^u + m^\theta$ with $m^\theta = bp(\ell + 1)$.

## 3 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 = 10^{-5}$ s as in the companion paper [17] and $m^u = 10$.

### 3.1 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 $[10^{-2}$ Hz, $10^5$ Hz] of the thermomechanical model. Each expansion is performed at order 0, that is $\ell = 0$ in eq. (14). As the test model
<figure><img src="image_4.png" /><figcaption>**FIGURE 3:** RCH ROM transfer function convergence with order $l = 0$ and $p$ expansion points. Full model [---], $p = 1$ [---], $3$ [---] and $7$ [---].</figcaption></figure>

possesses $b = 7$ boundary nodes, the size of the RCH models is $7(p + 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 = 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 = 1$) and multiple expansion orders $l$. As $l$ increases, the error decreases quickly. In both cases,

<figure><img src="image_5.png" /><figcaption>**FIGURE 4:** RCH ROM transfer function convergence with $p = 3$ expansion points and a $l$-th order expansion. Full model [---], $l = 0$ [---], $1$ [---], $2$ [---] and $3$ [---].</figcaption></figure>

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 = l = 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 $l = 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 = 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 $\Phi_{im}$ matrix of eq. (15) (SVD), or multiple Gram-Schmidt orthogonalizations (mGS)—one every time a new vector is appended to $\Phi_{im}$.

<figure><img src="image_6.png" /><figcaption>**FIGURE 5:** Comparison of CB [---], CH [---], RCH-mGS [---], RCH-SVD [---] and RCH-sGS [---] ROMs transfer functions with 70 modes; full model [---].</figcaption></figure>

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 ninth-order 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 $\Phi_{im}$ or a SVD. The RCH-sGS is used in the remaining.

## 3.2 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 $\theta$ and contact force $\lambda_N$. The only difference between the two simulations is that $K^{u\theta} = 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 [0.9 s; 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
<figure><img src="image_7.png" /><figcaption>FIGURE 6: Thermomechanical model responses at node 3. Uncoupled [—] and coupled [—] models.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7: Close-up of fig. 6.</figcaption></figure>

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 = 3$ and $p = 5$). For $p = 5$, the curves cannot be distinguished. The RCH model calculated with 3 expansion points shows large errors<fnref n="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 = 3$, the contact phases are shorter, the temperature and the normal force

<figure><img src="image_9.png" /><figcaption>FIGURE 8: RCH ROM convergence in contact simulation. Full model [—], RCH ROM with $ℓ = 0$, $p = 3$ [---] and 5 [- -].</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9: Close-up view of fig. 8.</figcaption></figure>

are underestimated.

The previous results have shown that the RCH model computed using reduction parameters ($ℓ = 0, p = 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 $ℓ = 0$ and $p = 5$ and a CB ROM comprising the same amount of DOFs ($m^θ = p \times b = 35$ and $m^u = 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 $\lambda_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.

<footnote n="2">²The error for $x$ is computed as $\frac{\int_0^1 |x_{ROM} - x_{FULL}|dt}{\int_0^1 |x_{FULL}|dt}$.</footnote>
<figure><img src="image_11.png" /><figcaption>**FIGURE 10:** CB and RCH ROM comparison in contact simulation. Full model [---], CB ($m^θ = p \times b = 35$) [---] and RCH ($ℓ = 0, p = 5$) [---].</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**FIGURE 11:** Close-up view of fig. 11.</figcaption></figure>

## 4 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

<figure><img src="image_13.png" /><figcaption>**FIGURE 12:** Finite element model of simplified industrial compressor sector. Boundary nodes [•] and constrained nodes [---].</figcaption></figure>

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

<figure><img src="image_14.png" /><figcaption>**FIGURE 13:** CB ROM transfer function convergence: full model [---], $m^θ = 27$ [---], 45 [---], 54 [---], 63 [---], 81 [---], 90 [---], 126 [---].</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**FIGURE 14:** RCH ROM transfer function convergence. Full model [---], $m^θ = 27$ [---], 45 [---], 54 [---], 63 [---], 81 [---], 90 [---] and 126 [---].</figcaption></figure>

sponds to $b = 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^θ = 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 = 10, m^θ = 126$). During this simulation, the blade is excited on its first bending mode frequency (507 Hz) and modeshape. Simulation
<figure><img src="image_16.png" /><figcaption>FIGURE 15: Time histories for node 1. Full model [—], RCH [—], CB ROM [—] with $m^u = 10$ and $m^θ = 126$.</figcaption></figure>

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 $\Phi_{im}^{\theta\theta}$ are indicated in fig. 16. All computations were

<figure><img src="image_17.png" /><figcaption>FIGURE 16: Computation times for various thermal reduction orders: CB [●] and RCH [+].</figcaption></figure>

performed on identical hardware using MATLAB® 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

<table><caption>TABLE 1: Computation times for 1 s of simulation. Full and reduced models with $m^u = 10$ and $m^θ = 126$ modes. Error term defined in footnote 2, page 6.</caption><thead><tr><th>Model</th><th>Full</th><th>CB ROM</th><th>RCH ROM</th></tr></thead><tbody><tr><td>Computation Time</td><td>37 236 s</td><td>14.431 s</td><td>13.969 s</td></tr><tr><td>Computation Speed</td><td>×1</td><td>×2580</td><td>×2665</td></tr><tr><td>Error on u<sub>r</sub></td><td>0%</td><td>6.16%</td><td>1.43%</td></tr><tr><td>Error on θ</td><td>0%</td><td>8.40%</td><td>0.45%</td></tr></tbody></table>

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

## 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

<table><tr><td>M<sup>uu</sup>, C<sup>uu</sup>, K<sup>uu</sup></td><td>Mass, damping and stiffness matrices</td></tr><tr><td>C<sup>θθ</sup>, K<sup>θθ</sup></td><td>Heat capacity and heat conductivity matrices</td></tr><tr><td>K<sup>uθ</sup></td><td>Thermoelastic coupling matrix</td></tr><tr><td>I, 0</td><td>Identity matrix, zero matrix</td></tr><tr><td>Ψ<sup>uu</sup><sub>ib</sub></td><td>Matrix of mechanical static constraint modes</td></tr><tr><td>Ψ<sup>θθ</sup><sub>ib</sub></td><td>Matrix of thermal static constraint modes</td></tr><tr><td>Φ<sup>uu</sup><sub>im</sub></td><td>Mechanical modeshapes matrix</td></tr><tr><td>Φ<sup>θθ</sup><sub>im</sub></td><td>Thermal modeshapes matrix</td></tr><tr><td>u</td><td>Vector of nodal displacements</td></tr><tr><td>θ</td><td>Vector of nodal temperatures</td></tr><tr><td>f, q</td><td>Vectors of nodal forces and nodal heat flows</td></tr><tr><td>h</td><td>Transfer function</td></tr><tr><td>•<sub>b</sub>, •<sub>i</sub></td><td>Boundary and Internal DOF</td></tr><tr><td>•<sub>m</sub></td><td>Modal DOF</td></tr><tr><td>u<sub>r</sub></td><td>Radial displacement</td></tr><tr><td>θ</td><td>Nodal temperature</td></tr><tr><td>λ<sub>N</sub>, λ<sub>T</sub></td><td>Normal and tangential contact forces</td></tr><tr><td>g<sub>N</sub>, g<sub>T</sub></td><td>Relative normal and tangential displacements</td></tr></table>
$\mu$ Coulomb friction coefficient

$s$ Laplace variable

$s_0$ Expansion point

$\ell$ Taylor series order of truncation

$p$ Number of expansion points

$m^u$ Number of mechanical dynamic modeshapes

$m^\theta$ 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

REFERENCES

[1] Jacquet-Richardet, G., Torkhani, M., Cartraud, P., Thouverez, F., Nouri-Baranger, T., Herran, M., Gibert, C., Baguet, S., Almeida, P., and Peletan, L., 2013. "Rotor to stator contacts in turbomachines. Review and application". *Mechanical Systems and Signal Processing*, **40**(2), pp. 401–420. [hal-00934050].

[2] Legrand, M., Batailly, A., and Pierre, C., 2011. "Numerical investigation of abradable coating removal in aircraft engines through plastic constitutive law". *Journal of Computational and Nonlinear Dynamics*, **7**(1). [hal-00627526].

[3] Batailly, A., Agrapart, Q., Millecamps, A., and Brunel, J.-F., 2016. "Experimental and numerical simulation of a rotor/stator interaction event localized on a single blade within an industrial high-pressure compressor". *Journal of Sound and Vibration*, **375**, pp. 308–331. [hal-01342401].

[4] Batailly, A., Legrand, M., Cartraud, P., and Pierre, C., 2010. "Assessment of reduced models for the detection of modal interaction through rotor stator contacts". *Journal of Sound and Vibration*, **329**(26), pp. 5546–5562. [hal-00524762].

[5] Goldman, P., and Muszynska, A., 1995. "Rotor-to-stator, rub-related, thermal/mechanical effects in rotating machinery". *Chaos, Solitons & Fractals*, **5**(9), pp. 1579–1601. [hal-01711340].

[6] Almeida, P., Gibert, C., Thouverez, F., Leblanc, X., and Ousty, J.-P., 2014. "Experimental analysis of dynamic interaction between a centrifugal compressor and its casing". *Journal of Turbomachinery*, **137**(3). [hal-01574149].

[7] Murthy, R., Wang, X. Q., Matney, A., and Mignolet, M., 2017. "A construction of thermal basis functions for coupled structural-thermal reduced order models". In 58th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference. [hal-01712954].

[8] Zukowski, E., Wilde, J., Rudnyi, E.-B., and Korvink, J.-G., 2005. "Model reduction for thermo-mechanical simulation of packages". In International Workshop on Thermal Investigation of ICs and Systems, pp. 134–138. [hal-00189465].

[9] Bechtold, T., Salimbahrami, B., Rudnyi, E., Lohmann, B., and Korvink, J. G. G., 2003. "Krylov-subspace-based order reduction methods applied to generate compact-electro-thermal models for mems". In NSTI Nanotechnology Conference & Trade Show Nanotech. [hal-01615922].

[10] Botto, D., Zucca, S., and Gola, M., 2007. "Reduced-order models for the calculation of thermal transients of heat conduction/convection fe models". *Journal of Thermal Stresses*, **30**(8), pp. 819–839. [hal-01712937].

[11] Nachtergaele, P., Rixen, D., and Steenhoek, A., 2010. "Efficient weakly coupled projection basis for the reduction of thermo-mechanical models". *Journal of Computational and Applied Mathematics*, **234**(7), pp. 2272–2278. [hal-01517953].

[12] Guérin, N., Thouverez, F., Gibert, C., Legrand, M., and Almeida, P., 2017. "Thermomechanical component mode synthesis for blade cas-

ing interaction prediction". In ASME Turbo Expo: Turbomachinery Technical Conference and Exposition. [hal-01569918].

[13] Grimme, J. E., 1997. "Krylov projection methods for model reduction". PhD thesis, University of Illinois at Urbana-Champaign. [tel-01711328].

[14] Freund, R. W., 2000. "Krylov-subspace methods for reduced-order modeling in circuit simulation". *Journal of Computational and Applied Mathematics*, **123**(1-2), pp. 395-421. [hal-01615923].

[15] Craig, R. R., and Bampton, M., 1968. "Coupling of substructures for dynamic analyses". *AIAA Journal*, **6**(7), pp. 1313-1319. [hal-01537654].

[16] Craig, R. R., and Hale, A. L., 1988. "Block-Krylov component synthesis method for structural model reduction". *Journal of Guidance, Control, and Dynamics*, **11**(6), pp. 562-570. [hal-01615926].

[17] Thorin, A., Guérin, N., Legrand, M., Thouverez, F., and Almeida, P., 2018. "Nonsmooth thermoelastic simulations of blade-casing contact interactions". In ASME Turbo Expo: Turbomachinery Technical Conference and Exposition. [hal-01713919].

[18] Padova, C., Barton, J., Dunn, M., and Manwaring, S., 2006. "Experimental results from controlled blade tip/shroud rubs at engine speed". *Journal of Turbomachinery*, **129**(4), pp. 713-723. [hal-01333702].

[19] Parent, M.-O., and Thouverez, F., 2016. "Phenomenological model for stability analysis of bladed rotor-to-stator contacts". In 16th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery. [hal-01537643].

[20] Legrand, M., Pierre, C., Cartraud, P., and Lombard, J.-P., 2009. "Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction". *Journal of Sound and Vibration*, **319**(1-2), pp. 366-391. [hal-00328186].

[21] Stewart, D., 2000. "Rigid-body dynamics with friction and impact". *SIAM review*, **42**(1), pp. 3-39. [hal-01570533].

[22] Acary, V., and Brogliato, B., 2008. *Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics*. Springer. [inria-00423530].

[23] Ireman, P., Klarbring, A., and Strömberg, N., 2002. "Finite element algorithms for thermoelastic wear problems". *European Journal of Mechanics-A/Solids*, **21**(3), pp. 423-440. [hal-01574157].

[24] Besselink, B., Tabak, U., Lutowska, A., van de Wouw, N., Nijmeijer, H., Rixen, D., Hochstenbach, M., and Schilders, W., 2013. "A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control". *Journal of Sound and Vibration*, **332**(19), pp. 4403-4422. [hal-01711355].

[25] Qu, Z.-Q., 2004. *Model Order Reduction Techniques*. Springer. [hal-01713836].

[26] Schilders, W., and Lutowska, A., 2014. "A novel approach to model order reduction for coupled multiphysics problems". In *Reduced Order Methods for Modeling and Computational Reduction*, Vol. 9 of *Modeling, Simulation and Applications book series*. Springer, pp. 1-49. [hal-01713843].

[27] MacNeal, R., 1971. "A hybrid method of component mode synthesis". *Computers & Structures*, **1**(4), pp. 581-601. [hal-01537661].

[28] Bladh, R., Pierre, C., and Castanier, M., 2003. "Numerical instability of classical free-interface component mode synthesis techniques". *AIAA Journal*, **41**(8), pp. 1621-1624. [hal-01711566].

[29] Guyan, R. J., 1965. "Reduction of stiffness and mass matrices". *AIAA Journal*, **3**(2), pp. 380-380. [hal-01711552].

[30] Arnoldi, W. E., 1951. "The principle of minimized iterations in the solution of the matrix eigenvalue problem". *Quarterly of Applied Mathematics*, **9**(1), pp. 17-29. [hal-01712943].

[31] Lanczos, C., 1950. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators". *Journal of Research of the National Bureau of Standards*, **45**(4), pp. 255-282. [hal-01712947].

[32] Bai, Z., and Su, Y., 2005. "Dimension reduction of large-scale second-order dynamical systems via a second-order Arnoldi method".
*SIAM Journal on Scientific Computing*, **26**(5), pp. 1692–1709. [hal-01711570].

[33] Elfadel, I. M., and Ling, D. D., 1997. “A block rational Arnoldi algorithm for multipoint passive model-order reduction of multi-port RLC networks”. In *IEEE/ACM International Conference on Computer-aided Design*, pp. 66–71. [hal-01711349].

[34] Pohle, M., Panning-von Scheidt, L., and Wallaschek, J., 2012. “Reduced order model of mistuned bladed disks using the Krylov sub-

space combined with the Craig-Bampton reduction”. In 10th World Congress on Computational Mechanics, Vol. 1, pp. 899–908. [hal-01712851].

[35] Cardona, A., and Idelsohn, S., 1986. “Solution of non-linear thermal transient problems by a reduction method”. *International Journal for Numerical Methods in Engineering*, **23**(6), pp. 1023–1042. [hal-01711528].