# Thermomechanical Component Mode Synthesis for blade casing interaction prediction

Nicolas Guérin, Fabrice Thouverez, Claude Gibert, Mathias Legrand, Patricio Almeida

## ► To cite this version:

Nicolas Guérin, Fabrice Thouverez, Claude Gibert, Mathias Legrand, Patricio Almeida. Thermo-mechanical Component Mode Synthesis for blade casing interaction prediction. ASME Turbo Expo 2017: Turbomachinery Technical Conference and Exposition, Jun 2017, Charlotte, United States. 10.1115/ GT2017-64342. hal-01569918

HAL Id: hal-01569918

https://hal.science/hal-01569918

Submitted on 27 Jul 2017

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 4.0 International License
# THERMOMECHANICAL COMPONENT MODE SYNTHESIS FOR BLADE CASING INTERACTION PREDICTION

Nicolas Guerin<sup>(a,c)</sup>, Fabrice Thouverez<sup>(a)</sup>, Claude Gibert<sup>(a)</sup>, Mathias Legrand<sup>(b)</sup>, Patricio Almeida<sup>(c)</sup>

(a) Ecole Centrale de Lyon, Laboratoire de Tribologie et Dynamique des Systemes
UMR CNRS 5513, 36 avenue Guy de Collongue, 69134 Ecully CEDEX France

(b) Structural Dynamics and Vibration Laboratory, Department of Mechanical Engineering
McGill University, 817 Sherbrooke St West, Montreal, Quebec, Canada H3A 0C3

(c) Safran Helicopter Engines, 64511 Bordes CEDEX France

Email: nicolas.guerin@ec-lyon.fr

## ABSTRACT

Increasing the efficiency of turbomachines is a major concern as it directly translates into lower environmental impact and improved operational costs. One solution is to reduce the blade-casing operating clearance in order to mitigate aerodynamic losses at the unavoidable cost of increased structural unilateral contact and friction occurrences. In centrifugal compressors, the dynamic behaviour of the structures interacting through unilateral contact and friction is not yet fully understood. In fact, the heat generated during such events may affect the dynamics through thermal stresses.

This paper presents a complete thermomechanical modelling strategy of impeller rotor and casing, and of blade-tip/casing contact events. A fully coupled thermomechanical modal synthesis technique is introduced and applied to turbomachinery-related models. The blisk is reduced via a hybrid modal synthesis technique combining the Craig-Bampton method and the characteristic constraint mode method. The casing model is reduced using an axisymmetric harmonic modal synthesis. Both strategies involve thermomechanical modes embedding thermal dilatation effects. The contact modelling algorithm is then introduced. It handles unilateral contact and friction occurrences together with heating effects. This algorithm uses the above mentioned reduced-order models as input data to avoid CPU-intensive simulations.

The results show that the thermomechanical behaviour of the structures is well preserved by the reduction strategy proposed. Contact simulations on simple cases show qualitative results in accordance with expectations. Further work is needed in order to validate the strategy based on experimental results. However, this methodology opens the way to extended multiphysics simulations of contact events in turbomachinery.

**Keywords** component mode synthesis, thermomechanical coupling, impeller, unilateral contact, friction, turbomachines, explicit time-marching technique

## INTRODUCTION

With a growing need for performance and power-to-weight ratios of aircraft and helicopter engines, designers make their best to improve the efficiency of every component. A solution is to reduce blade tip-casing operating clearances in order to optimize pressure ratios in the compressor and turbine stages [1]. Due to such clearance reductions, the occurrence of unilateral and friction contact phenomena between blades and casing becomes

a common situation during the lifetime of the engines. It is now recognized that such phenomena have to be accounted for in the design in order to avoid premature maintenance operation and therefore excessive operational costs. Several investigations, both numerical [2, 3] and experimental [4, 5] are devoted to the modelling and understanding of such events. A detailed summary [6] lists recent studies regarding material, thermal and mechanical aspects of rotor-stator rubbing.

Recent experimental campaigns [4, 5] have shown the importance of thermomechanical coupling effects during the interaction. However, very little documentation is available [7, 8] and thorough numerical investigations have yet to be undertaken. Still, thermomechanical contact models exist [9, 10] and related issues are explored in other industrial fields, for instance braking systems [11] or nuclear reactors [12]. This paper deals with the development of a reduction strategy taking into account the thermomechanical response of a bladed-disk and its matching casing during rubbing events, while maintaining low computational costs. The proposed strategy is based on [13], and takes advantage of cyclic and axisymmetric structures properties. A time-marching algorithm is also presented. It has been developed in order to model both contact forces and heat flows during rotor-stator interactions. First simulation results are then presented.

# 1 THERMOMECHANICAL COMPONENT MODE SYNTHESIS

## 1.1 Linear thermoelastic solids

This paper focuses on the modelling of linear thermoelastic solids. Other physical phenomena are neglected, as well as the variation of the material properties with temperature. Based on the first and second laws of thermodynamics for thermoelastic solids, the Clausius-Duhem inequality writes

$$ \bar{\bar{\sigma}} : \dot{\bar{\bar{\epsilon}}} - \rho(\dot{\psi} + s\dot{T}) - \frac{\bar{q}}{T} \overrightarrow{\text{grad}(T)} \geq 0 \quad (1) $$

where $\psi = e - Ts$ is the Helmholtz energy, $\bar{\bar{\sigma}}$, the stress tensor, $\bar{\bar{\epsilon}}$, the strain tensor, $\rho$, the material density, $s$, the entropy, and $T$, the absolute temperature. Vector $\bar{q}$ is an external thermal flux. Under the assumption of small perturbations of the system, and considering $\psi$ only dependent on the strain tensor and the temperature, the state laws of thermoelasticity

$$ \bar{\bar{\sigma}} = \frac{\partial \psi}{\partial \bar{\bar{\epsilon}}} \qquad s = - \frac{\partial \psi}{\partial T} \quad (2) $$
are formulated using Clausius-Duhem's inequality. A first order Taylor series expansion of the laws of thermoelasticity and Fourier's conduction law gives the expressions of the linear thermomeastic behaviour law

$$ \bar{\bar{\sigma}} = \bar{\bar{E}} : \bar{\bar{\varepsilon}} - \bar{\bar{K}}(T - T_0) \quad (3) $$

and of the heat equation

$$ \rho c \dot{T} = r + \lambda \Delta T - \bar{\bar{K}} : \dot{\bar{\bar{\varepsilon}}} T_0 \quad (4) $$

where $\bar{\bar{E}}$ stands for the tensor of elastic moduli and $\bar{\bar{K}}$ for the thermal expansion tensor. $T_0$ is the reference temperature, c is the material heat capacity, $\lambda$ its thermal conductivity, and r an external heat source. From the variational formulation of the above Partial differential equations, the following spatially discretized formulation is defined:

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

It is commonly used in commercial finite element software. The coupling of thermal and elasticity terms is dictated by the $\mathbf{K}^{u\theta}$ and $\mathbf{C}^{\theta u}$ matrix blocks, respectively standing for dilatation and thermoelastic damping effects. The thermoelastic damping is often neglected in mesoscopic scaled systems. In the case of microscopic systems, such as microelectromechanical systems (MEMS), it plays a greater role and has to be accounted for [14].

## 1.2 Craig-Bampton thermomechanical reduction

In order to perform a Craig-Bampton coupled reduction of the system in Eq. (5), matrices need to be partitioned into internal ($\bullet_i$) and boundary ($\bullet_b$) degrees-of-freedom (DOF). This applies to both thermal ($\mathbf{\theta}$) and mechanical ($\mathbf{u}$) DOF, and leads to the following matrix basis change:

$$ \begin{pmatrix} \mathbf{u} \\ \mathbf{\theta} \end{pmatrix} \rightarrow \begin{pmatrix} \mathbf{u}_b \\ \mathbf{\theta}_b \\ \mathbf{u}_i \\ \mathbf{\theta}_i \end{pmatrix}. \quad (6) $$

As in the standard Craig-Bampton technique [15], the method proposed in [13] uses both static and dynamic modeshapes in the reduction basis. Static modeshapes are calculated using unitary nodal displacements and temperatures successively imposed on each boundary node while blocking the other boundary degrees-of-freedom:

$$ \begin{bmatrix} \mathbf{I}_b \\ \mathbf{\Psi}_{ib} \end{bmatrix} \quad (7) $$

where $\mathbf{\Psi}_{ib}$ is defined as

$$ - \begin{bmatrix} \mathbf{K}_{ii}^{uu} & \mathbf{K}_{ii}^{u\theta} \\ \mathbf{0} & \mathbf{K}_{ii}^{\theta\theta} \end{bmatrix}^{-1} \begin{bmatrix} \mathbf{K}_{ib}^{uu} & \mathbf{K}_{ib}^{u\theta} \\ \mathbf{0} & \mathbf{K}_{ib}^{\theta\theta} \end{bmatrix} \quad (8) $$

This reduction basis contains thermomechanical coupling through the $\mathbf{K}^{u\theta}$ blocks. Calculating the boundary constrained modes of the system is less obvious. Two issues can be raised:

* Which modeshapes should be incorporated to capture the transient thermomechanical behaviour of the system?

* How could these modes be evaluated at a low computational cost?

In a standard mechanical reduction method, some eigenvectors would be selected based on a specific criterion, for example the eigenvalue. The reduction is then operated by only retaining the modeshapes associated with the lowest eigenvalues. In the present work, such a strategy would be too expensive. In fact, the eigenvector calculation of thermomechanical systems is ill-conditioned. To overcome this situation, a formulation of a coupled reduction basis from the uncoupled sub-problems is proposed in [13]. The dynamic reduction matrix is built by means of purely mechanical vibration modes $\mathbf{\Phi}_m^{uu}$, as well as purely thermal heat conduction modes $\mathbf{\Phi}_m^{\theta\theta}$. The coupling is operated using so-called “thermal-driven modes”

$$ \mathbf{\Phi}_m^{u\theta} = -(\mathbf{K}_{ii}^{uu})^{-1} \mathbf{K}_{ii}^{u\theta} \mathbf{\Phi}_m^{\theta\theta} \quad (9) $$

that may be interpreted as the displacements induced by thermal modeshapes due to dilatation effects. The dynamic reduction basis is then

$$ \begin{bmatrix} \mathbf{0} \\ \mathbf{\Phi}_m \end{bmatrix} = \begin{bmatrix} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \\ \mathbf{\Phi}_m^{uu} & \mathbf{\Phi}_m^{u\theta} \\ \mathbf{0} & \mathbf{\Phi}_m^{\theta\theta} \end{bmatrix}. \quad (10) $$

The full reduction matrix is built by concatenating the static and dynamic reduction matrices in Eq. (7) and (10). The above described thermomechanical reduction is based on the approach developed by Craig and Bampton, involving fixed interface and displacement static modes [15]. However, other component mode synthesis methods could be considered for reduction purposes.

# 2 TURBOENGINE IMPELLER

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

## 2.1 Proposed strategy

The reduction strategy developed is illustrated in Fig. 1. The first step consists in the application of the above thermomechanical Craig-Bampton reduction on the finite element model of a sector of the impeller, see Fig. 1(a). For this reduction, the retained degrees-of-freedom are the thermal and mechanical DOF on the cyclic symmetry boundaries of the sector as well as DOF of nodes on the blade tip. This reduction leads to the construction of a sector superelemen (SE) as shown in Fig. 1(b). This superelement is then replicated to reproduce all sectors of the blisk. Later on, superelements are assembled using displacement compatibility conditions between adjacent sectors. This operation already cuts half of the degrees of freedom of the attach faces, further reducing the model as illustrated in Fig. 1(c). However, since the amount of degrees-of-freedom on these faces is in general very high (especially when considering impeller rotors), a second reduction pass has to be performed.

The second reduction step lies on the interface or characteristic constraint modes [16, 17] reduction method. It is used to reduce interface DOF resulting from the assembly of superelements and is based on the eigenvectors of the interface subsystem. In this particular case, interfaces are the attached faces remaining once sector superelements have been assembled. This method was applied to thermomechanical problems using the thermal driven modes Eq. (9). Compared to the reduction of a full blisk
finite element model, the proposed strategy is expected to be more efficient. In fact, the matrices handled during calculations are much tinier, leading to both lower RAM and CPU needs.

<figure><img src="image_2.png" /><figcaption>FIGURE 1: Impeller reduction strategy</figcaption></figure>

## 2.2 Application model

The finite element model displayed in Fig. 2 and used to create the blisk superelement was built using ANSYS software and Solid226 coupled-field elements. In order to obtain a CPU efficient model, a simplified geometry as well as a coarse mesh were studied. Each blade of this model possesses a boundary node (i.e. retained during reduction) on its tip with three nodal displacements and a nodal temperature.

<figure><img src="image_3.png" /><figcaption>FIGURE 2: Dummy blisk model: Constrained nodes (●), Interface nodes (■), and Contact nodes (○)</figcaption></figure>

The system is constrained both mechanically and thermally on its internal bore. All other degrees-of-freedom are considered free of any constraint. Thus, external boundaries are perfectly insulated. Such conditions do not represent actual turbomachinery applications. However, test rigs dedicated to the study of blade-casing contact events on actual engine components work under low pressure conditions in order to avoid aerodynamic perturbations [8]. Such an environment leads to low convection exchange between the remaining fluid and the contacting structures. The constant temperature on the internal bore models the thermal exchange with a shaft.

## 2.3 Modal analysis

As stated before, it may be difficult to calculate fully coupled thermomechanical modes of large finite element models. Moreover, a thorough time-stepping comparison between a reduced

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

<figure><img src="image_4.png" /><figcaption>FIGURE 3: Convergence analysis for the mechanical subproblem</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4: Convergence analysis for the thermal subproblem</figcaption></figure>

tion strategy proposed allows to reduce models while preserving modal behaviours.
## 2.4 Transient analysis

Time domain transient analysis is now performed under various external load cases so as to ensure the quality of the thermomechanical coupling following the second reduction step. This transient analysis is performed on a blisk superelement with 10 structural interface modes and 10 thermal interface modes. A

<figure><img src="image_6.png" /><figcaption>FIGURE 5: Comparison procedure</figcaption></figure>

first validation step is done loading blade tips with unitary forces. Nodal displacements and temperatures are reported in Fig. 6. The displacements of boundary nodes seem to be well described under such mechanical loads. Discrepancies arise in nodal temperatures though. They are probably related to the coupling approximation. However, the effect of a mechanical load on boundary temperatures is practically zero in both cases. A second validation is com-

<figure><img src="image_7.png" /><figcaption>FIGURE 6: Blade tip nodal quantities under tip load: CB (—) and CB+IM (···)</figcaption></figure>

pleted with heat flows prescribed on the blade tip boundary node. The results are presented in Fig. 7. Nodal temperatures show a good correlation between the two compared models. Moreover, nodal displacements induced by the thermal loads are also well

described by the fully reduced model, ensuring that the second reduction preserves thermomechanical coupling effects.

<figure><img src="image_8.png" /><figcaption>FIGURE 7: Blade tip nodal quantities under tip heat flux: CB (—) and CB+IM (···)</figcaption></figure>

# 3 CASING THERMOMECHANICAL MODAL SYNTHESIS

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

## 3.1 Axis ans cyclically symmetric structures

Axisymmetric structures have geometrical properties such that problem dimensions in a Finite Element framework can be reduced to a minimum. In particular, such structures may be represented with a continuous spatial Fourier transform. It is then possible to travel between physical (x) and Fourier (x̄) coordinates via the transformation

$$x(\beta) = \sum_{k=0}^{\infty} e^{ik\beta} \tilde{x}^k \quad (11)$$

where β is the angular position and k is the harmonic index. This approach can be exploited with thermal and thermomechanical models. Therefore, a thermomechanical equivalent of Eq. (11) is

$$\left(\frac{\mathbf{u}}{\boldsymbol{\theta}}\right)(\beta) = \sum_{k=0}^{\infty} e^{ik\beta} \left(\tilde{\mathbf{u}}^k, \tilde{\boldsymbol{\theta}}^k\right). \quad (12)$$

An analogous formulation can be employed for cyclically symmetric structures by replacing the continuous spatial Fourier transform with a discrete one

$$\left(\frac{\mathbf{u}}{\boldsymbol{\theta}}\right)(n) = \sum_{k=0}^{N-1} e^{ikn\alpha} \left(\tilde{\mathbf{u}}^k, \tilde{\boldsymbol{\theta}}^k\right). \quad (13)$$

In this formulation, the angular position is replaced by the sector number n multiplied with α, the sector angle.
## 3.2 Modal analysis

According to [18], the eigenproblem of the whole structure is equivalent to the set of independent eigenproblems associated to each harmonic index k. For each of them, the uncoupled thermal and mechanical modes are calculated. The thermal-driven modes in Eq. (9) are inserted to relate both physics. Therefore, the modal shapes for a given k are

$$ \tilde{\mathbf{\Phi}}^k = \begin{bmatrix} \tilde{\mathbf{\Phi}}^{uu,k} & \tilde{\mathbf{\Phi}}^{u\theta,k} \\ \mathbf{0} & \tilde{\mathbf{\Phi}}^{\theta\theta,k} \end{bmatrix}. \quad (14) $$

The method provides access to the thermomechanical modal response of an axisymmetric structure based on a single cross-section. Using the previous modal basis onto which the dynamics of the casing is projected, Fourier coordinates are expressed as

$$ \left( \begin{matrix} \tilde{\mathbf{u}}^k \\ \tilde{\mathbf{\theta}}^k \end{matrix} \right) (t) = \tilde{\mathbf{\Phi}}^k \tilde{\mathbf{q}}_k (t). \quad (15) $$

Finally, via Eq. (12), physical coordinates are expressed as a linear combination of time-dependent contributions of modeshapes:

$$ \left( \begin{matrix} \mathbf{u} \\ \mathbf{\theta} \end{matrix} \right) (t, \beta) = \sum_{k=0}^{K} \tilde{\mathbf{\Phi}}^k \tilde{\mathbf{q}}_k(t) e^{ik\beta} \quad (16) $$

Equation (16) is displayed graphically in Fig. 8.

<figure><img src="image_9.png" /><figcaption>FIGURE 8: Casing reduced-order model strategy</figcaption></figure>

## 3.3 Application: impeller casing model

The casing on which the method was applied is displayed in Fig. 9(a). Red dots indicate fixed degrees-of-freedom. As an example, Fig. 9(b) displays the first three nodal-diameters modeshapes.

<figure><img src="image_10.png" /><figcaption>FIGURE 9: Casing application model</figcaption></figure>

# 4 ROTOR-STATOR INTERACTION SIMULATION

## 4.1 Continuous formulations

A modelling of the thermomechanical behaviour of a contact interface has been proposed by several authors [9, 10, 19]. This formulation completes the equations governing the contacting solids by introducing an interface potential $\psi_I$ and another

Clausius-Duhem inequality governing the thermodynamics of the interface:

$$ \psi_I \leq \mathbf{f}_N \dot{\mathbf{u}}_N + \mathbf{f}_T \dot{\mathbf{u}}_T - s_I \dot{\mathbf{\Theta}}_I + \frac{\mathbf{q}_{c,1}}{\mathbf{\Theta}_1} (\mathbf{\Theta}_1 - \mathbf{\Theta}_I) + \frac{\mathbf{q}_{c,2}}{\mathbf{\Theta}_2} (\mathbf{\Theta}_2 - \mathbf{\Theta}_I) \quad (17) $$

In this case, the interface law only accounts for unilateral contact, Coulomb friction and frictional heating. The phenomena are described by the normal ($\mathbf{f}_N, \mathbf{u}_N$) and tangential ($\mathbf{f}_T, \mathbf{u}_T$) contact forces and relative displacements, the interface absolute temperature $\mathbf{\Theta}_I$ and entropy $s_I$; $\mathbf{\Theta}_1$ and $\mathbf{\Theta}_2$ are the absolute temperatures of the contacting bodies, $\mathbf{q}_{c,1}$ and $\mathbf{q}_{c,2}$ are the contact heat flux. Other phenomena may be considered in the contact law, such as adhesive or abrasive wear [20].

## 4.2 Resolution scheme

Contact problems are mainly solved using two strategies: time-marching methods or harmonic balance methods. Direct integration methods allow to study transient phenomena. Such methods are already in use to study rotor-stator contact problems, e.g. central finite differences [21, 3]. Harmonic methods are used to model periodic steady-state phenomena and have been applied to study blade-root / disk contacts [22] and friction phenomena in general [23].

Regarding thermomechanically coupled problems, two approaches can be highlighted, either tightly or loosely coupled. In the first case, the fully coupled thermoelastic system is solved in one step. In the other case, successive resolutions of the mechanical and thermal problems are performed. A recent evaluation of both methods [12] shows that simulation performances greatly depend on the size and type of problem explored. In particular, for thermomechanical contact problems, the tight coupling seems to be more efficient.

In this study, a strongly coupled resolution is carried out, in accordance with the reduction strategy presented. In order to solve the dynamic system, a scheme based on the weighed residuals method is used. It matches the central finite difference method used in [2], when the integration parameter $\alpha$ is of 1/2. However, such a parameter value leads to an unconditionally unstable scheme due to the thermal part of the model. For values above 1/2 the scheme is conditionally stable. Practically, $\alpha$ is chosen close to 1/2 to avoid excessive numerical damping. The time discretisation leads to

$$ \begin{align} \left( \frac{\mathbf{M}}{\Delta t^2} + \frac{\alpha \mathbf{C}}{\Delta t} \right) \left( \begin{matrix} \mathbf{u} \\ \mathbf{\theta} \end{matrix} \right)_{n+1} &= \left( \frac{1}{2} + \alpha \right) \left( \begin{matrix} \mathbf{f} \\ \mathbf{q} \end{matrix} \right)_n + \left( \frac{1}{2} - \alpha \right) \left( \begin{matrix} \mathbf{f} \\ \mathbf{q} \end{matrix} \right)_{n-1} \nonumber \\ &\quad + \left( \frac{2\mathbf{M}}{\Delta t^2} + \frac{2\alpha - 1}{\Delta t} \mathbf{C} - \left( \frac{1}{2} + \alpha \right) \mathbf{K} \right) \left( \begin{matrix} \mathbf{u} \\ \mathbf{\theta} \end{matrix} \right)_n \tag{18} \\ &\quad - \left( (\alpha - 1) \frac{\mathbf{C}}{\Delta t} + \frac{\mathbf{M}}{\Delta t^2} - \left( \frac{1}{2} - \alpha \right) \mathbf{K} \right) \left( \begin{matrix} \mathbf{u} \\ \mathbf{\theta} \end{matrix} \right)_{n-1} \nonumber \end{align} $$

where $\Delta t$ is the time increment and $n$ is the increment number.

## 4.3 Contact event handling

In structural dynamics, these problems are often tackled using penalty or Lagrange multiplier methods [24, 25, 2], the latter having the advantage of preventing residual penetration between the contacting bodies. The algorithm used is based on a forward Lagrange multiplier prediction-correction [24, 2]. The gap $g_p$ separating the two structures is first calculated without contact conditions. If a negative gap is detected, Lagrange multipliers (i.e. normal contact forces) are calculated according to

$$ \lambda_N = \left[ \mathbf{C}_{N,p} \left( \frac{\mathbf{M}}{\Delta t^2} + \frac{\alpha \mathbf{C}}{\Delta t} \right)^{-1} \mathbf{C}_{NT,p} \right]^{-1} \mathbf{g}_p. \quad (19) $$
The $C_{N,p}$ matrix contains the predicted direction of the normal contact force. The $C_{NT,p}$ matrix stores the predicted direction of the contact force taking into account both normal and tangential forces. The displacement corrections $\delta u$ are then calculated according to

$$ \delta \mathbf{u} = - \left( \frac{\mathbf{M}}{\Delta t^2} + \frac{\alpha \mathbf{C}}{\Delta t} \right)^{-1} \mathbf{C}_{NT,p} \boldsymbol{\lambda}_N. \quad (20) $$

In the present work, the friction law is regularized using

$$ \mathbf{f}_T = \mu \boldsymbol{\lambda}_N \tanh \left( \frac{\dot{\mathbf{u}}_T}{\varepsilon} \right) \quad (21) $$

where $\mu$ is the friction coefficient, $\boldsymbol{\lambda}_N$, the normal contact force, and $\dot{\mathbf{u}}_T$, the relative velocity between a pair of points located on the rotor and on the stator. The regularization is controlled by $\varepsilon$.

Thermal considerations in a contact simulation lead to the formulation of heat flows generated at the contact interface. They are induced by heat creation at the interface (depending on the contact law) and by conduction between the two structures. Therefore, they are very much influenced by the surface roughness of the two structures (impacting the thermal contact conductance), the friction coefficient and local temperatures. The contact heat flows take the following form:

$$ \begin{aligned} \mathbf{q}_{c,1} &= \frac{v_1}{v_1 + v_2} (v_2(\boldsymbol{\Theta}_2 - \boldsymbol{\Theta}_1) - \mu \|\dot{\mathbf{u}}_T\|) \boldsymbol{\lambda}_N \\ \mathbf{q}_{c,2} &= \frac{v_2}{v_1 + v_2} (v_1(\boldsymbol{\Theta}_2 - \boldsymbol{\Theta}_1) - \mu \|\dot{\mathbf{u}}_T\|) \boldsymbol{\lambda}_N \end{aligned} \quad (22) $$

where $v_1$ and $v_2$ are the thermal contact conductances of the bodies coming into contact. Based on the displacements correction Eq. (20) in the purely mechanical case, the following formulation is proposed

$$ \left( \frac{\delta \mathbf{u}}{\delta \boldsymbol{\theta}} \right) = - \left( \frac{\mathbf{M}}{\Delta t^2} + \frac{\alpha \mathbf{C}}{\Delta t} \right)^{-1} \mathbf{C}_{NT\theta,p} \boldsymbol{\lambda}_N \quad (23) $$

for the thermomechanical case, where $C_{NT\theta,p}$ gathers the participations of the normal and tangential forces as well as the contact heat flux. Thus, this strategy consists in an isothermal prediction step followed by a thermomechanical correction step. The

<figure><img src="image_11.png" /><figcaption>FIGURE 10: Algorithm organisation</figcaption></figure>

corresponding algorithm is presented in Fig. 10.

## 4.4 Excitation-driven contact

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

<figure><img src="image_12.png" /><figcaption>FIGURE 11: Contact Configuration: Stator (●) and Rotor (■)</figcaption></figure>

At $t = 4$ s, the excitation force is stopped. Responses presented in Figs 12 and 13 show the gap closing due to the increasing excitation level.

At $t \approx 0.3$ s, the gap closes. When the external excitation stops, nodal displacements converge towards a non-zero value due to the thermomechanical coupling. The spectrogram of the excited node axial displacement in Fig 13(a) shows that the frequency content of the displacement is directly related to the excitation

<figure><img src="image_13.png" /><figcaption>FIGURE 12: Excitation-driven contact
Blade 1 (●), 2 (■), 3 (■)</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIGURE 13: Time-frequency analysis of blade 1 tip displacement</figcaption></figure>

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

## 4.5 Displacement-driven contact

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

## 5 CONCLUSION

Thermomechanical model reduction methods have been proposed and applied on turbomachinery related structures in order to provide efficient models for thermomechanical rub events simulations. These methods show accurate modal characteristics, compared with reference models, while displaying a minimal amount of degrees-of-freedom. Results of time integration simulations

<figure><img src="image_15.png" /><figcaption>FIGURE 14: Excitation-driven contact
Blade 1 (—), 2 (—), 3 (—)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 15: Time-frequency analysis of blade 1 tip displacement</figcaption></figure>
show that the thermal behaviour and thermomechanical coupling are well described by the reduction methods, thus allowing to link both physics during simulations. The currently developed contact algorithm, based on a forward increment Lagrange multiplier method is able to capture mechanical and thermal phenomena emerging during contact events. Further studies are however needed to extensively compare the results with other numerical tools and experimental measurements.

Applying such methods on large finite element models may lead to error propagation issues and instability of the reduced-order models. Therefore, such applications require a careful implementation in order to obtain light but trustworthy models. However, low number of remaining degrees-of-freedom, the CPU cost remains limited when compared to initial full models.

Current developments attempt to extend the proposed approach to industrial models, opening the way to multiphysics simulations of rotor-stator rub events in turbomachinery.

## ACKNOWLEDGMENT

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

## NOMENCLATURE

σ̄, σ̄̄ Stress and strain tensors

ρ Material density

ψ Helmholtz energy

s Entropy

T Absolute temperature

q̄ Heat flow vector

K̄ Thermal coefficients tensor

T₀ Initial temperature of the material

c, λ Heat capacity and Heat conductivity

r Radiative source

M<sup>uu</sup>, C<sup>uu</sup>, K<sup>uu</sup> Mass, damping and stiffness matrices

C<sup>θθ</sup>, K<sup>θθ</sup> Heat capacity and heat conductivity matrices

C<sup>θu</sup> Thermoelastic damping matrix

K<sup>uθ</sup> Thermoelastic matrix

I Identity matrix

Ψ<sub>ib</sub> Thermomechanical static modeshapes matrix

Φ<sup>uu</sup> Mechanical modeshapes matrix

Φ<sup>uθ</sup> "Thermal-driven" modeshapes matrix

Φ<sup>θθ</sup> Thermal modeshapes matrix

u Nodal relative displacement vector

θ Nodal relative temperature vector

f Nodal force

q Nodal heat flow

•b Boundary DOF

•i Internal DOF

β Angular position on axisymmetric structure

k Harmonic index

• Fourier DOF

ψ<sub>I</sub> Interface Helmholtz energy

λ<sub>N</sub> Normal contact Lagrange multiplier

f<sub>N</sub> Normal contact force

f<sub>T</sub> Tangential contact force

u<sub>N</sub> Relative normal displacement

u<sub>T</sub> Relative tangential displacement

s<sub>I</sub> Interface entropy

T<sub>I</sub> Interface temperature

q<sub>c</sub> Contact heat flow

Δt Time step

n Time iteration

g<sub>p</sub> Vector of predicted gaps

δu Displacement correction due to contact constraints

δθ Temperature correction due to contact heat flux

ν Thermal contact conductance

μ Coulomb friction coefficient

ε Regularization parameter

## REFERENCES

[1] Danish, S. N., Qureshi, S. R., Imran, M. M., Khan, S. U.-D., Sarfraz, M. M., El-Leathy, A., Al-Ansary, H., and Wei, M., 2016. "Effect of tip clearance and rotor-stator axial gap on the efficiency of a multistage compressor". *Applied Thermal Engineering*, **99**, pp. 988–995. [hal-01574153].

[2] Legrand, M., 2005. "Modèles de prédiction de l'interaction rotor/stator dans un moteur d'avion". PhD thesis, Ecole Centrale de Nantes, Nantes, France. [tel-00011631].

[3] Batailly, A., Meingast, M., and Legrand, M., 2015. "Unilateral contact induced blade/casing vibratory interactions in impellers: Analysis for rigid casings". *Journal of Sound and Vibration*, **337**, pp. 244–262. [hal-01120157].

[4] Millecamps, A., 2010. "Interaction aube-carter: contribution de l'usure de l'abradable et de la thermomécanique sur la dynamique d'aube". PhD thesis, Université de Lille 1, Lille, France. [tel-01200536].

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

[6] Millecamps, A., Batailly, A., Legrand, M., and Garcin, F., 2015. "Snecma's viewpoint on the numerical and experimental simulation of blade-tip/casing unilateral contacts". In Proceedings of the ASME Turbo Expo 2015: Turbine Technical Conference and Exposition, no. GT2015-42682 in Volume 7B: Structures and Dynamics. [hal-01223582].

[7] Millecamps, A., Brunel, J.-F., Dufrenoy, P., Garcin, F., and Nucci, M., 2009. "Influence of thermal effects during blade-casing contact experiments". In Proceedings of the ASME 2009 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, no. DETC2009-86842 in Volume 1: 22nd Biennial Conference on Mechanical Vibration and Noise, Parts A and B, pp. 855–862. [hal-01223060].

[8] Almeida, P., 2014. "Couplage thermomécanique et vibra-toire d'un compresseur centrifuge lors d'un contact aube-carter". PhD thesis, Ecole Centrale de Lyon, Lyon, France. [tel-01289766].

[9] Wriggers, P., and Miehe, C., 1994. "Contact constraints within coupled thermomechanical analysis—A finite element model". *Computer Methods in Applied Mechanics and Engineering*, **113**(3), pp. 301–319. [hal-01574175].

[10] Strömberg, N., 1999. "Finite element treatment of two-dimensional thermoelastic wear problems". *Computer Methods in Applied Mechanics and Engineering*, **177**(3–4), pp. 441–455. [hal-01574166].
[11] Majcherczak, D., Dufrenoy, P., and Berthier, Y., 2007. "Tri-bological, thermal and mechanical coupling aspects of the dry sliding contact". *Tribology International*, **40**(5), pp. 834–843. [hal-00378119].

[12] Novascone, S. R., Spencer, B. W., Hales, J. D., and Williamson, R. L., 2015. "Evaluation of coupling approaches for thermomechanical simulations". *Nuclear Engineering and Design*, **295**, pp. 910–921. [hal-01574160].

[13] Nachtergaele, P., Rixen, D. J., and Steenhoek, A. M., 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].

[14] Duwel, A., Gorman, J., Weinstein, M., Borenstein, J., and Ward, P., 2003. "Experimental study of thermoelastic damping in MEMS gyros". *Sensors and Actuators A: Physical*, **103**(1–2), pp. 70–75. [hal-01574211].

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

[16] Castanier, M. P., Tan, Y.-C., and Pierre, C., 2001. "Characteristic Constraint Modes for Component Mode Synthesis". *AIAA Journal*, **39**(6), pp. 1182–1187. [hal-01537673].

[17] Tran, D.-M., 2009. "Component mode synthesis methods using partial interface modes: Application to tuned and mistuned structures with cyclic symmetry". *Computers & Structures*, **87**(17–18), pp. 1141–1153. [hal-01537665].

[18] Thomas, D. L., 1979. "Dynamics of rotationally periodic structures". *International Journal for Numerical Methods in Engineering*, **14**(1), pp. 81–102. [hal-01574169].

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

[20] Berthoul, B., Batailly, A., Legrand, M., Stainier, L., and Cartraud, P., 2015. "Abradable coating removal in turbomachines: A macroscopic approach accounting for several wear mechanisms". In *Proceedings of the ASME Turbo Expo 2015: Turbine Technical Conference and Exposition*, no. GT2015-42500 in Volume 7B: Structures and Dynamics. [hal-01188846].

[21] Legrand, M., Batailly, A., Magnain, B., Cartraud, P., and Pierre, C., 2012. "Full three-dimensional investigation of structural contact interactions in turbomachines". *Journal of Sound and Vibration*, **331**(11), May, pp. 2578–2601. [hal-00660863].

[22] Salles, L., 2010. "Etude de l'usure par fretting sous chargements dynamiques dans les interfaces frottantes : application aux pieds d'aubes de turbomachines". PhD Thesis, Ecole Centrale de Lyon, Lyon, France. [tel-00600613].

[23] Zucca, S., and Firrone, C. M., 2014. "Nonlinear dynamics of mechanical systems with friction contacts: Coupled static and dynamic Multi-Harmonic Balance Method and multiple solutions". *Journal of Sound and Vibration*, **333**(3), pp. 916–926. [hal-01574178].

[24] Carpenter, N. J., Taylor, R. L., and Katona, M. G., 1991. "Lagrange constraints for transient finite element surface contact". *International Journal for Numerical Methods in Engineering*, **32**(1), pp. 103–128. [hal-01389918].

[25] Wriggers, P., and Laursen, T., 2007. *Computational contact mechanics*. No. 498 in CISM courses and lectures. Springer.