# Nonsmooth thermoelastic simulations of blade-casing contact interactions

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

► To cite this version:

Anders Thorin, Nicolas Guérin, Mathias Legrand, Fabrice Thouverez, Patricio Almeida. Nonsmooth thermoelastic simulations of blade-casing contact interactions. ASME Turbo Expo 2018: Turbomachinery Technical Conference and Exposition, Jun 2018, Oslo, Norway. 10.1115/GT2018-75959. hal-01713919

HAL Id: hal-01713919

https://hal.science/hal-01713919

Submitted on 24 Feb 2018

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
NONSMOOTH THERMOELASTIC SIMULATIONS OF
BLADE–CASING CONTACT INTERACTIONS

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

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

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

(c) Safran Helicopter Engines, Bordes, France

Email: anders.thorin@mcgill.ca

ABSTRACT

In turbomachinery, it is well known that tighter operating clearances improve the efficiency. However, this leads to unwanted potential unilateral and frictional contact occurrences between the rotating (blades) and stationary components (casings) together with attendant thermal excitations. Unilateral contact induces discontinuities in the velocity at impact times, hence the terminology nonsmooth dynamics. Current modeling strategies of rotor-stator interactions are either based on regularizing penalty methods or on explicit time-marching methods derived from Carpenter's forward Lagrange multiplier method. Regularization introduces an artificial time scale in the formulation corresponding to numerical stiffness which is not desirable. Carpenter's scheme has been successfully applied to turbomachinery industrial models in the sole mechanical framework, but faces serious stability issues when dealing with the additional heat equation.

This work overcomes the above issues by using the Moreau-Jean nonsmooth integration scheme within an implicit θ-method. This numerical scheme is based on a mathematically sound description of the contact dynamics by means of measure differential inclusions and enjoys attractive features. The procedure is unconditionally stable opening doors to quick preliminary simulations with time-steps one hundred times larger than with previous algorithms. It can also deal with strongly coupled thermomechanical problems.

**Keywords:** nonsmooth dynamics, unilateral contact, friction, rotor-stator interaction

INTRODUCTION

Turbomachinery equipment manufacturers tend to improve the aerodynamic performance by reducing operating gaps between rotating and stationary components. This leads to frequent contact occurrences, referred to as rotor-stator interaction, in compressor and turbine stages of modern aircraft and helicopter engines [1]. To reduce the possibly harmful effects of such events on the global dynamics and subsequent damage, casings are commonly coated with abradable materials [1]. This strategy has proven efficient, however the resulting dynamics, coupling contact and frictional heating, is still not well understood and modelled.

Various developments have been dedicated to the improvement of turbomachines design to account for blade-tip rubbing, modal interaction or whip-whirl motions [1]. In these works, contact is dealt with using either a regularized approach [2] or the forward Lagrange multiplier formulation [3, 4]. The regularized approach consists in modeling the contact between the rotor and the stator through a stiff spring. Though very simple to implement, the additional and somehow artificial spring leads to stiff

numerical problems which are prone to stability issues [5, 6]. The Lagrange multiplier formulation preserves the nonsmooth framework induced by contact conditions but uses an explicit procedure in time for which stability is guaranteed only for sufficiently small time-steps. This is particularly problematic since thermal fluxes generated by frictional heating seem to play a significant role in rotor-stator interactions [7] and are thus considered in the model. Indeed, explicit schemes are not well-suited for the heat equation as they lead to stiff numerical problems [8]. More specifically, the stability is typically governed by a relationship of the form $\Delta t \le \Delta x^2/2$, implying that a refinement of the space discretization has drastic consequences on the necessary refinement in time.

A simple thermo-mechanical model of a bladed-disk sector section 1 is first considered in this work. The coupling is two-fold: from heat diffusion to mechanics via thermal expansion and from mechanics to heat through a frictional heating law. The governing equations are solved via dedicated nonsmooth solvers which rely on an implicit numerical scheme [6] as described in section 2. Time-domain histories, systematic comparison with the forward Lagrange multiplier method as well as sensitivity analysis to time-step are discussed in section 3. Finally, the methodology is adapted to an industrial-size model in section 4.

1 DESCRIPTION OF THE SIMPLIFIED MODEL

Due to the high computational costs induced by non-linear phenomena in large-scale problems, model-order reduction methods are often implemented for realizable computations, as exposed in section 4. However, the simple one-sector model illustrated in fig. 1 is first used to validate the methodology. The stator is rigid

<figure><img src="image_2.png" /><figcaption>FIGURE 1: Simplified finite element model. Contact nodes [●], constrained nodes [■] and corresponding contact forces [→].</figcaption></figure>

and the rotor includes 2028 degrees of freedom (507 free nodes in each space direction and temperature) and seven contact points. It possesses constrained nodes (zero displacement and prescribed
temperature of 25 °C) on its internal bore and contact nodes on the blade tip. The structural dynamics is governed by an equation of the form

$$M^{uu}\ddot{u} + C^{uu}\dot{u} + K^{uu}u + K^{u\theta}\theta = f^u + f_c^u \quad (1)$$

where **u** is the vector of generalized displacements, **θ** is the vector of generalized temperatures, **M**<sup>uu</sup>, **C**<sup>uu</sup> and **K**<sup>uu</sup> are respectively the mass, damping and stiffness matrices stemming from the finite element discretization. The temperature-induced expansions are accounted for through the coupling matrix **K**<sup>uθ</sup>. The external forces are split into the external load **f**<sup>u</sup> and the contact forces **f**<sub>c</sub><sup>u</sup>. In a more elaborate model, **f**<sup>u</sup> could include centrifugal or aerodynamic effects while contact could be induced by the static deflection or the vibrations of the stator. In this work, the load **f**<sup>u</sup>, described further, activates contact. To every contact node i = 1, ..., 7 corresponds a gap function g<sub>i</sub>(**u**) measuring the algebraic distance to the stator and an impulsive contact force λ<sub>i</sub> in the inward-pointing normal direction. The following so-called Signorini conditions complements eq. (1):

$$g_i \ge 0 \quad (2a)$$

$$\lambda_i \ge 0 \quad (2b)$$

$$g_i \lambda_i = 0. \quad (2c)$$

The first two inequalities enforce the non-penetration and the non-sticking condition, respectively. Equality (2c) translates the fact that the gap is open ($g_i > 0$) if and only if the associated contact force is zero, except in the very specific case $g_i = 0$ and $\lambda_i = 0$. The seven Signorini conditions (one for each contact node) are usually gathered in the compact form $0 \le g \perp \lambda \ge 0$ where $g = [g_1, \dots, g_7]^T$, $\lambda = [\lambda_1, \dots, \lambda_7]^T$ and the operators $\perp$, $\perp$ and $\ge$ are defined component-wise. In addition to Signorini conditions, as in any space semi-discretized contact formulation, an impact law is required to ensure uniqueness of the solution [9]. In order to mimic lasting contact phases, we choose an inelastic Newton impact law: $y_i^+ = 0$ where $y_i^+ = \dot{g}_i(u(t^+))$ denotes the normal velocity right after the impact time. In the tangential direction, $f_c^u$ includes Coulomb friction, therefore proportional to the normal forces. This can be written as $f_c^u = C_{NT}^u \lambda$ where $C_{NT}$ is a rectangular matrix transferring the normal impulses $\lambda$ in the local coordinates to the generalized coordinates in the normal direction, as well as in the tangential direction through the Coulomb friction coefficient. In the present case, the local nodal reference frames at the contact nodes are oriented according to the global cylindrical directions so that $C_{NT}^u$ takes the form

$$C_{NT}^u = I_7 \otimes \begin{bmatrix} -1 \\ -\mu R \Omega \\ 0 \end{bmatrix} \quad (3)$$

where R is the stator radius and Ω the angular velocity—the contribution of the blade vibration in the tangential contact velocity is neglected. Note that because of the high tangential contact velocity, there is no stick-slip transitions and friction is therefore not source of nonsmoothness here. Also, eq. (1) does not include centrifugal stiffening or spin softening. They could however be incorporated in a straightforward manner.

The thermal dynamics is governed by an equation of the form

$$C^{\theta\theta}\dot{\theta} + K^{\theta\theta}\theta = f^{\theta} + f_c^{\theta} \quad (4)$$

where $C^{\theta\theta}$ and $K^{\theta\theta}$ are the heat capacity and heat conductivity matrices. As for the structure, the external load is separated into external fluxes $f_c^\theta$ and frictional heating due to contact $f_c^u$. Each of the seven frictional heat fluxes $\lambda_i^\theta$ follows a simple law $\lambda_i^\theta = \alpha\lambda_i$ which are included in $f_c^\theta$ through a relationship of the form $f_c^\theta = C_{NT}^\theta \lambda$ where $C_{NT}^\theta = \alpha I_7$. This means that each contact occurrence generates a heat flux proportional to the normal contact force. This simple law [10] is sufficient to illustrate the capabilities of the proposed numerical method. However it is worth mentioning that it might induce temperature discontinuities. Indeed, if one $\lambda_i$ is impulsive, eq. (4) shows that $\dot{\theta}$ will also be corresponding to a discontinuous $\theta$. Integrating more complex laws is straightforward. Also, in contrast to other available investigations [11, 12], contact-induced wear in rotor-stator interactions are not considered in the present model.

Altogether, eq. (1) and eq. (4) are compactly recast in the form

$$M\ddot{x} + C\dot{x} + Kx = f + f_c \quad (5)$$

with

$$M = \begin{bmatrix} M^{uu} & 0 \\ 0 & 0 \end{bmatrix}, \quad C = \begin{bmatrix} C^{uu} & 0 \\ 0 & C^{\theta\theta} \end{bmatrix}, \quad K = \begin{bmatrix} K^{uu} & K^{u\theta} \\ 0 & K^{\theta\theta} \end{bmatrix} \quad (6)$$

and

$$x = \begin{bmatrix} u \\ \theta \end{bmatrix}, \quad f = \begin{bmatrix} f^u \\ f^\theta \end{bmatrix}, \quad f_c = \begin{bmatrix} f_c^u \\ f_c^\theta \end{bmatrix}. \quad (7)$$

The matrices **M** and **K** were obtained using SOLID226 coupled-field elements in ANSYS®. The damping matrix verifies $C^{uu} = 10^{-5}K^{uu}$. Contact is activated via $f^u$, chosen here as a sinusoidal function of frequency 33 Hz and amplitude 100 N, in the radial direction and on all contact nodes. The frequency corresponds to $\Omega = 2000$ rpm, assuming the blade hits the stator once per rotation. For comparison, the first flexural mode of the blade has a frequency of 1300 Hz. The contact frequency could also be due to the vibration of the blade, as in section 4.

A TA6V titanium alloy, widely used for the production of aircraft compressor components, is considered. The material properties, estimated from several alloys properties, read as follows: density 4430 kg m⁻³, Young's modulus 110 GPa, Poisson's ratio 0.3, heat capacitance 520 J K⁻¹ kg⁻¹, heat conductivity 6.7 J K⁻¹ m⁻¹, dilatation coefficient 9 µm m⁻¹ K⁻¹. The coefficient of frictional heating is $\alpha = 0.1$ and the coefficient of friction in Coulomb's law is $\mu = 0.15$.

# 2 NONSMOOTH NUMERICAL METHODS

## 2.1 Formulation using Differential Measures

Equation (1) may seem to be an Ordinary Differential Equation (ODE). More rigorously, it should be described by means of weaker mathematical objects such as distributions or measures. Indeed, $f_c^u$ contains impulsive terms: every time a gap closes, the corresponding contact node instantly loses all its kinetic energy. This instantaneous loss of kinetic energy is mirrored by an impulse that is proportional to a Dirac distribution. Equation (1) is therefore an abuse of notation in the context of nonsmooth dynamics. In this subsection, the model is reformulated in a sounder mathematical manner, using measure differential inclusions, making the use of efficient dedicated numerical methods possible.
Introducing differential measures dv and dθ [6] and the Lebesgue measure dt, eq. (5) can be expressed as

$$ \{ \mathbf{M} \mathrm{d}v + (\mathbf{C}v + \mathbf{K}x - \mathbf{f}) \mathrm{d}t = \mathrm{d}\mathbf{r} \quad (8a) $$

$$ \{ \mathbf{v} = \dot{\mathbf{x}}. \quad (8b) $$

In short, the prefactors of dt are smooth terms, while dv and dr are measures which, when integrated between two arbitrarily close times, can yield non-zero values: they contain impulses, which are proportional to Dirac delta distributions. The measure dr stores all contact forces and frictional heating. Loosely speaking, the reader who is not familiar with measures can compare eq. (8a) to eq. (5) multiplied by dt.

Using convex analysis, the Signorini conditions (2) together with the inelastic impact law mentioned above can be formulated using the proj operator and the tangent cone $\mathcal{T}$ in $\mathbb{R}^+$ [6] as:

$$ \mathbf{y}^{+} = \mathrm{proj}(\mathcal{T}_{\mathbb{R}^{+}}(\mathbf{g}(\mathbf{u})), \mathbf{y}^{-}) \quad (9) $$

where the equality is componentwise and the vector y stacks the normal contact velocities $[\dot{g}_1, \dots, \dot{g}_7]^T$. The contact impulses are related to y through

$$ \lambda = \hat{\mathbf{W}}^{-1}(\mathbf{y}^{+} - \mathbf{y}^{-}) \quad \text{and} \quad \hat{\mathbf{W}} = (\nabla_{\mathbf{u}}\mathbf{g})^{\top}\mathbf{M}^{-1}(\nabla_{\mathbf{u}}\mathbf{g}). \quad (10) $$

Let us inspect eq. (9) for the ith unilateral constraint. By definition, the tangent cone reads as

$$ \mathcal{T}_{\mathbb{R}^+}(g_i(\mathbf{u})) = \begin{cases} \mathbb{R} & \text{if } g_i(\mathbf{u}) > 0 \\ \mathbb{R}^+ & \text{if } g_i(\mathbf{u}) = 0 \end{cases} \quad (11) $$

and the following can be said.

* During a free flight, $g_i(\mathbf{u}) > 0$ so $\mathbf{y}^+ = \mathrm{proj}(\mathbb{R}, \mathbf{y}^-) = \mathbf{y}^-$ and from eq. (11), $\lambda_i = 0$ follows: there is no reaction impulse and velocity is continuous.

* At the beginning of a contact phase, $g_i(\mathbf{u}) = 0$ and $y_i^- \le 0$ so that $y_i^+ = 0$ and the impact law is recovered. The corresponding impulse is $\lambda = -\hat{\mathbf{W}}^{-1}\mathbf{y}^- \ge 0$.

* In the middle of a contact phase, $g_i(\mathbf{u}) = 0$ and $\mathbf{y}^- = \mathbf{0}$ so that $\mathbf{y}^+ = \mathbf{0}$ and there is no reaction impulse.

* At the end of a contact phase, $g_i(\mathbf{u}) = 0$ and $y_i^- < 0$ so that $\mathbf{y}^+ = \mathrm{proj}(\mathbb{R}^+, \mathbf{0}) = \mathbf{0}$ and $\lambda_i = 0$. The contact force vanishes and the gap opens.

To summarize, the problem is governed by eq. (8), eq. (9), the frictional heating law $\mathbf{f}_c^\theta = \mathbf{C}_{\mathrm{NT}}^\theta \lambda$ as well as the geometric relationships $\mathbf{f}_c^u = \mathbf{C}_{\mathrm{NT}}^u \lambda$ and $\mathbf{y} = \nabla_{\mathbf{u}}\mathbf{g}\dot{\mathbf{x}}$.

## 2.2 Time discretization

We now proceed with the Moreau-Jean discretisation [13, 14] of eq. (8), using a linearly implicit solver. For $n \in \mathbb{N}^*$, let $t_n$ denote the nth time-step such that $t_n = nh$ where h is the time-step. Integrating eq. (8) between $t_n$ and $t_{n+1}$ yields:

$$ \int_{t_n}^{t_{n+1}} \mathbf{M} \mathrm{d}v + \int_{t_n}^{t_{n+1}} (\mathbf{C}v + \mathbf{K}x - \mathbf{f}) \mathrm{d}t = \int_{t_n}^{t_{n+1}} \mathrm{d}\mathbf{r} \quad (12a) $$

$$ \mathbf{v} = \dot{\mathbf{u}}. \quad (12b) $$

Smooth quantities are discretized by means of the $\theta$-method, where to avoid confusion with the temperature $\theta$, the numerical parameter

is denoted by $\gamma$:

$$ \int_{t_n}^{t_{n+1}} \mathbf{v} \mathrm{d}t = h(\gamma \mathbf{v}_{n+1} + (1-\gamma)\mathbf{v}_n) := h \mathbf{v}_{n+\gamma} \quad (13) $$

and similarly for $\mathbf{x}$ and $\mathbf{f}$, so that eq. (12a) becomes

$$ \{ \mathbf{M}(\mathbf{v}_{n+1} - \mathbf{v}_n) + h(\mathbf{C}\mathbf{v}_{n+\gamma} + \mathbf{K}\mathbf{x}_{n+\gamma} - \mathbf{f}_{n+\gamma}) = \mathbf{p}_{n+1} \quad (14a) $$

$$ \mathbf{v}_n = \dot{\mathbf{u}}_n \quad (14b) $$

with $\mathbf{p}_{n+1} = \int_{t_n}^{t_{n+1}} \mathrm{d}\mathbf{r}$. Introducing $\hat{\mathbf{M}} = (\mathbf{M} + h\gamma\mathbf{C} + h^2\gamma^2\mathbf{K})$, eq. (14a) is rewritten as

$$ \hat{\mathbf{M}}(\mathbf{v}_{n+1} - \mathbf{v}_n) = -h(\mathbf{C}\mathbf{v}_n + \mathbf{K}\mathbf{u}_n - \mathbf{f}_{n+\gamma}) - h^2\gamma\mathbf{K}\mathbf{v}_n + \mathbf{p}_{n+1}. \quad (15) $$

We now introduce a predictor step, corresponding to the dynamics without the unilateral constraints. The purpose of this step is to estimate whether the contact status is about to change in the next time-step. Two new quantities are defined: the predicted gap chosen as $\tilde{\mathbf{g}}_{n+1} = \mathbf{g}(\mathbf{u}_n + h\mathbf{v}_n/2)$ and the predicted velocity $\tilde{\mathbf{v}}_{n+1}$ such that

$$ \hat{\mathbf{M}}(\mathbf{v}_{n+1} - \tilde{\mathbf{v}}_{n+1}) = \mathbf{p}_{n+1}. \quad (16) $$

The predicted normal velocities in the local coordinates are given by $\tilde{\mathbf{y}}_{n+1} = \nabla_{\mathbf{u}}\mathbf{g}\tilde{\mathbf{v}}_{n+1}$. Equation (9) is discretized as

$$ \mathbf{y}_{n+1} = \mathrm{proj}(\mathcal{T}_{\mathbb{R}^+}(\tilde{\mathbf{g}}_{n+1}), \tilde{\mathbf{y}}_{n+1}). \quad (17) $$

Then, the reactions in the global frame $\mathbf{p}_{n+1}$ are computed from the normal contact reactions in the local frame $\lambda_{n+1} = \hat{\mathbf{W}}^{-1}(\mathbf{y}_{n+1} - \mathbf{y}_n)$ through the geometric relationship $\mathbf{p}_{n+1} = \mathbf{C}_{\mathrm{NT}}^u \lambda_{n+1}$. The velocity is finally updated by adding the effects of contact to the contact-less prediction $\tilde{\mathbf{v}}_{n+1}$:

$$ \mathbf{v}_{n+1} = \tilde{\mathbf{v}}_{n+1} + \hat{\mathbf{M}}^{-1}\mathbf{p}_{n+1}. \quad (18) $$

## 3 VALIDATION USING THE SIMPLE MODEL

In this section, the reference solution is generated via the methodology described in the previous section, with a time-step $h = 1/3 \times 10^{-6}$. The error of a data series is defined as the integral of the absolute value of the difference with the reference solution, divided by the integral of absolute value of the reference solution, over the time interval $[0 \text{ s}, 1 \text{ s}]$. For example, a time series of the temperature of the first node $\theta_1$ shows the error

$$ \eta = \frac{\int_{[0,1]} |\theta_1 - \theta_1^{\mathrm{ref}}|}{\int_{[0,1]} |\theta_1^{\mathrm{ref}}|}. \quad (19) $$

Although imperfect, this definition is useful to quantify discrepancies in the remaining.

### 3.1 Time-step sensitivity analysis

Figure 2 displays the time evolution of the radial displacement, the temperature and the normal contact force for the middle of the blade tip, that is, contact node 4 in fig. 1, and different time-steps. For $h = 10^{-4}$ s, the error is about 3%, while temperature and contact force show 16% and 25% errors, respectively. After one second, the temperature is underestimated by about 5%, resulting from an underestimation of the normal contact force. This
<figure><img src="image_3.png" /><figcaption>**FIGURE 2:** Sensitivity to time-step at node 4 in terms of radial displacement $u_r$, temperature $\theta$ and normal contact force $\lambda$. Reference [—], $h = 10^{-4}$ s [—] and $h = 10^{-5}$ s [—].</figcaption></figure>

can lead to significant errors on longer simulations, but might be sufficient if contact events are expected to be of short durations. With $h = 10^{-5}$ s, the errors reduce to 0.21 %, 0.18 % and 1.72 %: the curves cannot be distinguished from the reference. Due to the strong thermomechanical coupling, the temperature oscillates at the loading frequency of 33 Hz. Heat accumulates and the temperature globally increases. As a result, the blade expands through the coupling term $K u^\theta$ and the gap tends to reduce while the normal force tend to increase. A close-up view on the last time range 0.9 s to 1 s is provided in fig. 3. Note that the simulated temperatures of

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

thousands of degrees go far beyond the material capabilities: that is an obvious shortcoming of the simple chosen model, which does not include nonlinear effects (other than contact) nor damage, but is sufficient for the present scope: validating the methodology. A more realistic model is presented in section 4.

The normal contact force, always non-negative, is seen to be non-zero only when the gap is closed, in accordance with the Signorini conditions eq. (2). Even with relatively large time-steps, the position exhibits clear kinks when the gap opens and closes, corresponding to velocity discontinuities induced by the inelastic impact law, features which would not be accurately captured via a classical time-domain integration scheme dedicated to smooth dynamics.

The proposed algorithm does not strictly prevent penetration, as illustrated in fig. 4. However, the residual penetration, which

<figure><img src="image_5.png" /><figcaption>**FIGURE 4:** Sensitivity of the residual penetration and contact force to time-step. Reference [—], $h = 10^{-4}$ s [—] and $h = 10^{-5}$ s [—].</figcaption></figure>

is small compared to the geometrical tolerances of the machine components, tends to zero as the time-step reduces. Also, once contact is activated, penetration never increases.

To conclude on the convergence analysis, the error is depicted as a function of the time-step in fig. 5. Given that the Carpenter's algorithm is unstable with $h > 5 \times 10^{-7}$ s (see section 3.3), the proposed method enable the use of time-steps larger by several of orders of magnitude.

<figure><img src="image_6.png" /><figcaption>**FIGURE 5:** Error at node 4: $u_r$ [—], $\lambda$ [· · ·] and $\theta$ [· · ·]</figcaption></figure>

## 3.2 Effect of coupling on contact geometry

The previous figures were depicted for contact node 4 (middle of the blade tip), which is the first one to close the gap. In this section, it is shown that thermomechanical coupling can have an effect on the contact geometry. Figure 6 compares the simulation results between node 4 and its neighbour node 3, as defined in fig. 1. For
<figure><img src="image_7.png" /><figcaption>**FIGURE 6:** Time histories for contact nodes 3 [---] and 4 [---].</figcaption></figure>

both nodes, the contact force $\lambda$ is non-zero if and only if the gap is closed. The middle node touches the stator on a longer interval than node 3, resulting in a normal force of larger amplitude and higher temperature due to frictional heating. At the beginning, all boundary nodes impact the stator. After a few contact occurrences, the middle of the blade has expanded so much due to heating that the nodes located at the end of the blade tip (nodes 1 and 7) separate from the stator. This is a logical consequence of the simple model and tend to validate the simulation methods. Every time the gap $g_4$ closes and opens, the trajectory of node 3 displays a small kink because of the structural coupling $M^{uu}$ and $K^{uu}$. In contrast, when thermomechanical coupling is ignored, the discrepancies between normal forces and displacements tend to reduce, see fig. 7. The

<figure><img src="image_8.png" /><figcaption>**FIGURE 7:** Time histories for contact nodes 3 [---] and 4 [---] without thermomechanical coupling.</figcaption></figure>

effect of temperature on the normal force appears clearly when comparing the bottom plot with that of fig. 6: thermal expansion further magnifies gap closure and attendant contact force.

## 3.3 Comparison with Carpenter's scheme

Results are now compared to the ones obtained with Carpenter's forward Lagrange multiplier method. Figure 8 shows that displacements and normal forces match well, with overall relative errors of 0.82% on nodal displacement, 5.1% on temperature and 7.3% on

contact force. Figure 9 is a close-up of fig. 8. The major benefit

<figure><img src="image_9.png" /><figcaption>**FIGURE 8:** Comparison of time histories for node 4. Moreau-Jean [---] and Carpenter [---].</figcaption></figure>

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

of the proposed method is its unconditional stability while Carpenter's scheme is unstable whenever $h > 4.6 \times 10^{-7}$ s. It enables the use of time-steps larger by several orders of magnitude.

# 4 INDUSTRIAL APPLICATION

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

## 4.1 Model

The model illustrated in fig. 10 comprises 5686 nodes and a total of 18768 DOFs (including 4692 thermal DOFs), with 9 contact
<figure><img src="image_11.png" /><figcaption>**FIGURE 10:** Finite element model of simplified industrial compressor sector. Contact nodes [•] and constrained nodes [■]. Colouring refers to the radial displacement.</figcaption></figure>

nodes located on the middle plane of the blade tip. In order to obtain reasonable computational times (of the order of the minute), the model is first reduced using a thermomechanical model-order reduction method developed by the authors [15]. The contact simulations were performed by periodically forcing the blade on its first flexural mode at frequency 507 Hz to reflect an aerodynamic excitation. The forcing amplitude was chosen in order to obtain approximately a steady-state radial displacement amplitude of $10^{-4}$ m on node 1 without contact.

## 4.2 Sensitivity to model parameters

To emphasize the effect of friction-induced heating during rotor-stator interactions in turbomachinery, a contact simulation is performed with and without frictional heating. The structure is initially

<figure><img src="image_12.png" /><figcaption>**FIGURE 11:** Time histories for node 1 with [—] and without [—] frictional heating.</figcaption></figure>

at rest. The results are depicted for node 1 in fig. 11. The response curves oscillate very fast because the simulation time is much larger than the excitation period.

During the first contacts, little difference is observed between the two responses. Then, the gap in the simulation with heating tends to reduce, while the contact force increases. At the end of the simulation, discrepancies of the order of 25% on the displacement and 30% on the contact force are predicted. This qualitative analysis shows that thermal effects have noticeable influence during the contacts, even for short durations.

The nonlinear nature of the system, stemming from the contact conditions, can be emphasized by doubling the excitation force amplitude, see fig. 12. Magnification on the time interval [0 s, 0.05 s]

<figure><img src="image_13.png" /><figcaption>**FIGURE 12:** Time histories for node 1. Effect of excitation amplitude: forcing of 1 [—] and 2 [—] normalized amplitude.</figcaption></figure>

is provided in fig. 13. As expected, the first contact occurs earlier

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

when the external load is magnified, due to higher vibration amplitudes. The higher excitation level also leads to an increased contact force and thereby to higher temperatures, inducing slightly smaller displacements. The temperature seems to diverge with the doubled load, probably because frictional heating generates more heat than the blade can diffuse. A thermoelastic instability emerges: it could advance the wear of machine components or even lead to more extensive damages.
## 4.3 Sensitivity to time-step

For this model, the stability condition required by Carpenter's algorithm is $h = 5 \times 10^{-7}$ s with which Moreau-Jean and Carpenter algorithms were implemented. The responses for node 1 are displayed in fig. 14. Again, they are in very good agreement.

<figure><img src="image_15.png" /><figcaption>**FIGURE 14:** Contact simulation results for node 1. Carpenter algorithm [---] and Moreau algorithm [---] for $h = 5 \times 10^{-7}$ s.</figcaption></figure>

Figure 15 shows results for various time-steps with Moreau-Jean's procedure only. The responses with $h = 10^{-5}$ s and

<figure><img src="image_16.png" /><figcaption>**FIGURE 15:** Moreau-Jean's scheme sensitivity to time-step for node 1: $h = 5 \times 10^{-7}$ s [---], $10^{-6}$ s [---], $5 \times 10^{-6}$ s [---] and $10^{-5}$ s [---].</figcaption></figure>

$h = 5 \times 10^{-7}$ s are almost identical, meaning that convergence is achieved for $h = 10^{-5}$ s. Therefore, the stability condition of Carpenter's algorithm is unnecessarily limiting. Moreover, for the same time-step $h = 5 \times 10^{-7}$ s, Carpenter's algorithm took about

twice as long as Moreau-Jean's algorithm in our implementation— this might be due to implementation difference and may not be true in general. More computation times are reported in fig. 16. All computations were performed using MATLAB® on a 4-core

<figure><img src="image_17.png" /><figcaption>**FIGURE 16:** Computation times for a 1 s simulation. Carpenter algorithm [---] versus Moreau-Jean algorithm [---].</figcaption></figure>

computer (2.3 GHz clock frequency).

## CONCLUSION

Conventional time-domain integration schemes to simulate rotor-stator interaction thermomechanical formulations face numerical stability issues: regularized contact constraints lead to numerical stiffness while explicit solution methods such as Carpenter's forward Lagrange multiplier scheme are not suited for the heat equation. The methodology proposed in the present paper, based on the Moreau-Jean procedure for measure differential inclusions addresses these shortcomings. It is shown to be very robust on two models of different sizes. While Carpenter's stability requires time-steps of the order of $5 \times 10^{-7}$ s, the thermomechanical Moreau-Jean scheme generates reasonable responses with time-steps of the order of $10^{-4}$ s for the simplified model and $10^{-5}$ s for the industrial model, offering significant gains in the computational effort. With $h = 10^{-5}$ s, the error in displacement, temperature and normal contact force is approximately 2%. The methodology is also compatible with model reduction [15] offering additional computation gains of several orders of magnitude and application to industrial-size models, including those with complex geometries.

## 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>Structural 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</td><td>Identity matrix</td></tr><tr><td>0</td><td>Zero matrix</td></tr><tr><td>u, θ</td><td>Nodal displacement and temperature vectors</td></tr><tr><td>x</td><td>Generalized degrees of freedom vector</td></tr><tr><td>v</td><td>Generalized velocity vector</td></tr><tr><td>M, C, K</td><td>Generalized mass, damping, stiffness matrices</td></tr><tr><td>f<sup>u</sup>, f<sup>θ</sup></td><td>Nodal force and heat flow vectors</td></tr><tr><td>λ</td><td>Normal contact force (Lagrange multiplier)</td></tr><tr><td>g</td><td>Gap vector</td></tr></table>
<table><tr><td>y</td><td>Relative velocity vector</td></tr><tr><td>μ</td><td>Coulomb friction coefficient</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></td><td>Normal nodal contact force</td></tr><tr><td>h</td><td>Time-step size</td></tr><tr><td>γ</td><td>θ-method numerical parameter</td></tr><tr><td>proj(K, y)</td><td>Projection operator of y on K</td></tr><tr><td>T<sub>K</sub>(x)</td><td>Tangent cone to K evaluated at x</td></tr><tr><td>⊗</td><td>Kronecker tensor product</td></tr></table>

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

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

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

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

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

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

[8] Shampine, L., and Gear, C., 1979. "A user's view of solving stiff ordinary differential equations". *SIAM Review*, **21**(1), pp. 1–17. [hal-01711390].

[9] Glocker, C., 2006. "An introduction to impacts". In *Nonsmooth mechanics of solids*, Vol. 485 of CISM International Centre for Mechanical Sciences book series. Springer, pp. 45–101. [hal-01716872].

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

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

[12] Berthoul, B., Batailly, A., Legrand, M., Stainier, L., and Cartraud, P., 2015. "Abradable coating removal in turbomachines: a macroscopic approach accounting for various wear mechanisms". In ASME Turbo Expo: Turbomachinery Technical Conference and Exposition, Vol. 7B: Structures and Dynamics. [hal-01188846].

[13] Moreau, J.-J., 1988. "Unilateral contact and dry friction in finite freedom dynamics". In *Nonsmooth Mechanics and Applications*, Vol. 302 of International Centre for Mechanical Sciences book series. Springer, pp. 1–82. [hal-01713847].

[14] Jean, M., 1999. "The non-smooth contact dynamics method". *Computer Methods in Applied Mechanics and Engineering*, **177**(3), pp. 235–257. [hal-01390459].

[15] Guérin, N., Thorin, A., Thouverez, F., Legrand, M., and Almeida, P., 2018. "Thermomechanical model reduction for efficient simulations of rotor-stator contact interaction". In ASME Turbo Expo: Turbomachinery Technical Conference and Exposition. [hal-01713849].