# Rotor-stator interaction: formulation, reduction and simulation of a nonsmooth thermo-elastic model

Anders Thorin, Nicolas Guérin, Mathias Legrand

► To cite this version:

Anders Thorin, Nicolas Guérin, Mathias Legrand. Rotor-stator interaction: formulation, reduction and simulation of a nonsmooth thermo-elastic model. 17th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery, Dec 2017, Maui, United States. hal-01871945

HAL Id: hal-01871945

https://hal.science/hal-01871945

Submitted on 11 Sep 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 4.0 International License
# Rotor-stator interaction: formulation, reduction and simulation of a nonsmooth thermo-elastic model

Anders Thorin<fnref n="1" />*, Nicolas Guérin<fnref n="2" />, Mathias Legrand<fnref n="1" />

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

ISROMAC 2017

International Symposium on Transport Phenomena and Dynamics of Rotating Machinery

Maui, Hawaii

December 16-21, 2017

**Abstract**

In turbomachines, heat fluxes are generated during unilateral contact and frictional rotor-stator interactions. Though they are expected to play a significant role in the overall dynamics of engines, their precise influence in the dynamics is not well understood. This contribution proposes a blade-casing reduced-order thermomechanical model involving unilateral contact and frictional terms formulated as a Measure Differential Inclusion together with various numerical methods able to properly capture the dynamics. Comparison with less advanced explicit time-marching schemes is provided. It is shown that the MDI discretized via an implicit scheme is more robust. A sensitivity analysis to the thermal coupling and model reduction is provided. It is found that the proposed solution method sets a sound basis for efficient numerical simulations of thermomechanical rub occurrences in turbomachinery.

**Keywords**

turbomachines, rotor-stator interaction, thermomechanics, unilateral contact, friction, nonsmooth dynamics

## INTRODUCTION

Due to the need for more efficient aircraft and rotorcraft engines, unilateral contact induced structural blade-to-casing interactions due to minimal operating clearances are known to be conceivable events. These rotor-stator contact interactions have been investigated exclusively from the purely mechanical standpoint where heat generation is neglected. However, heat fluxes generated during the contact occurrences are believed to affect the dynamics to an important degree.

Rotor-stator contact interactions are mathematically described by unilateral constraints preventing penetration between the rotating blades and stationary casings. Such complementarity conditions can be simply incorporated into the model via regularization [1,2]. However, this introduces difficulties such as numerical stiffness requiring artificially small time steps in the time-marching solution method, low local order of consistency, low global order of accuracy, as well as stability issues [3]. Another possibility is to use Carpenter's numerical scheme [4] consisting in a prediction step to detect interpenetration and a projection step to satisfy the unilateral constraint when required. Though acceptable in general, this method does not rely on a strong mathematical ground and faces stability issues. For example, Carpenter's scheme is explicit while such schemes should be avoided when solving the heat equation as they lead to stiff numerical problems [5] and therefore require very small time steps in order to achieve numerical stability. Additionally, Carpenter's scheme hides the impact law, required in a space

semi-discrete setting, offering no control on it. The authors have also experienced artificial energy increase due to the projection step during contact phases.

This contribution investigates the dynamics of a simplified blade impacting a rigid body. The dynamics is described by a Measure Differential Inclusion [3]. The scalability of this approach to industrial problems involving a large number of degrees-of-freedom is explored via a Craig-Bampton-like model reduction technique. Both the full and the reduced-order models are solved in the time domain using dedicated nonsmooth numerical schemes [3]. The effect of thermomechanical coupling, the number of degrees-of-freedom of the reduced model as well as numerical robustness are investigated.

## 1. THERMOELASTICITY FORMULATION

The spring-mass system shown in Fig. 1 is considered. The

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Simplified blade-casing thermoelastic model.</figcaption></figure>

presented methodology can be straightforwardly applied to more complex systems including Finite Element models and 3D geometry, provided the mass matrix is lumped on the contacting degree-of-freedom (dof). The nth mass, repre-

<footnote n="1">¹Departement of Mechanical Engineering, McGill University, Montreal, Canada</footnote>

<footnote n="2">²Laboratoire de Tribologie et Dynamique des Systèmes, École Centrale de Lyon, Ecully, France & Safran Helicopter Engines, Bordes, France</footnote>
senting the tip of the blade, undergoes a unilateral contact condition corresponding to an obstacle, reflecting the casing. Its time-dependent thermoelastic behaviour is described by

$$ \left\{ \begin{array}{l} u_{,tt} - c^2 u_{,xx} + \theta = f \\ \theta_{,t} - \theta_{,xx} + u_{,xt} = g \end{array} \right. \qquad (1a) $$

(1b)

complemented with the contact conditions, where $u$ is the relative displacement around the resting position $x$, and $\theta$ the temperature rise from an initial absolute temperature $T$. In all this work, the following finite element discretization of the aforementioned continuous problem is considered:

$$ \left\{ \begin{array}{l} M\ddot{u} + C^{uu}\dot{u} + K^{uu}u + K^{u\theta}\theta = f^u + f_{\text{contact}} \\ C^{\theta\theta}\dot{\theta} + K^{\theta\theta}\theta = f^\theta + f_{\text{tip}} \end{array} \right. \qquad (2a) $$

$$ C^{\theta\theta}\dot{\theta} + K^{\theta\theta}\theta = f^{\theta} + f_{\text{tip}} \qquad (2b) $$

where $u$ now denotes the vector of displacements, $\theta$ the vector of temperatures, $M$, $C^{uu}$ and $K^{uu}$ are the mass, stiffness and damping matrices, $C^{\theta\theta}$ and $K^{\theta\theta}$ are the heat capacity and heat conductivity matrices, $f_{\text{contact}}$, $f^u$ and $f^\theta$ are the unilateral and frictional contact forces, external structural loading and external heat loading, respectively. The coupling between mechanical and thermal quantities is twofold:

* through the coupling matrix $K^{u\theta}$, accounting for the temperature-induced mechanical expansions;

* through the heat production at the blade tip $f_{\text{tip}}$. Indeed, when the blade and the casing are in contact, friction generates heat. This is here incorporated via a heat flux on the last node which is proportional to the normal contact force: $f_{\text{tip}} = \alpha f_{\text{contact}}$.

Contact acts only on the last node which is unilaterally restricted: $f_{\text{contact}} = \lambda e^{(n)}$ where $\lambda$ is the normal contact force satisfying the Signorini conditions $\lambda \ge 0$, $d - u^{(n)} \ge 0$ and $\lambda(u^{(n)} - d) = 0$ which reflect no sticking, no penetration and no simultaneously vanishing normal contact force and separating clearance. The vector $e^{(n)}$ is the nth vector of the canonical basis (i.e. $[0 \dots 0 1]^T$), superscripts $\bullet^{(n)}$ denote the nth component obtained from the inner product with $e^{(n)}$, $d$ is the initial gap between the two contacting bodies. A Newton impact law is chosen: $v^{(n)+} = -e \cdot v^{(n)-}$ where $e$ is the coefficient of restitution, $v$ the vector of velocities ($v = \dot{u}$); $v^{(n)-}$ and $v^{(n)+}$ denote the pre- and post-impact velocities of the contacting dof. Eq. (2) is to be understood in a weak sense, as the velocity of mass $n$ is discontinuous at impact times, while $\lambda$ and the heat flux involve Dirac impulses. The complementarity conditions together with the impact law can be rigorously expressed by means of a Measure Differential Inclusion (MDI) [3]:

$$ -\lambda \in \partial \psi_{\mathcal{T}_{\mathbb{R}^+}(d-u^{(n)})} \left( - \frac{v^{(n)+} + e v^{(n)-}}{1 + e} \right) \qquad (3) $$

where $\partial \psi_{\mathcal{T}_{\mathbb{R}^+}(d-u^{(n)})}$ is the subderivative of the indicator function of the tangent cone to $\mathbb{R}^+$ evaluated at $d - u^{(n)}$. This concept and terminology can be confusing and a few explanations to better grasp the underlying meaning of this MDI

are now given; more details can be found in [3]. The tangent cone of $\mathbb{R}^+$ at $d - u^{(n)}$ is by definition [6]

$$ \mathcal{T}_{\mathbb{R}^+}(d - u^{(n)}) = \begin{cases} \mathbb{R} & \text{if } d - u^{(n)} > 0 \\ \mathbb{R}^+ & \text{if } d - u^{(n)} = 0 \end{cases} \qquad (4a) $$

$$ \text{ if } d - u^{(n)} = 0 \qquad (4b) $$

and it can be shown that $\partial \psi_{\mathcal{T}_{\mathbb{R}^+}(d-u^{(n)})}(y)$ is equal to the following normal cone:

$$ \mathcal{N}_{\mathcal{T}_{\mathbb{R}^+}(d-u^{(n)})}(y) = \begin{cases} \{0\} & \text{if } d - u^{(n)} > 0, \\ \{0\} & \text{if } d - u^{(n)} = 0 \text{ and } y > 0, \\ \mathbb{R}^{-} & \text{if } d - u^{(n)} = 0 \text{ and } y \le 0. \end{cases} \qquad (5a) $$

$$ \text{ if } d - u^{(n)} = 0 \text{ and } y > 0, \qquad (5b) $$

$$ \text{ if } d - u^{(n)} = 0 \text{ and } y \le 0. \qquad (5c) $$

If contact is not activated, $d - u^{(n)} > 0$ and from (3) and (5a), $\lambda \in \{0\}$: there is no contact force and the velocity is continuous since $\lambda = m_n(v^{(n)+} - v^{(n)-}) = 0$ by definition of the differential measure [3]. During a contact phase, $d - u^{(n)} = 0$ and $y = -v^{(n)} = 0$ so $\lambda \in \mathbb{R}^+$: the contact force is such that there is no sticking in the normal direction. At the end of a contact phase, $d - u^{(n)} = 0$ and $y = -v^{(n)} > 0$ so $\lambda \in \{0\}$: the contact force vanishes. The most interesting case is at the beginning of a contact phase (impact): then $d - u^{(n)} = 0$ and $v^{(n)-} > 0$. The impulse is proportional to the velocity discontinuity, the left-hand side (LHS) can be multiplied by $(m_n(1+e))^{-1}$ in the inclusion because the left-hand side (RHS) is a cone, and it follows from (3) and (4b) that

$$ \frac{v^{(n)-} - v^{(n)+}}{1 + e} = v^{(n)-} - \frac{v^{(n)+} + e v^{(n)-}}{1 + e} \in \partial \psi_{\mathbb{R}^{+}} \left( - \frac{v^{(n)+} + e v^{(n)-}}{1 + e} \right). \qquad (6) $$

It can be shown through a case-by-case discussion that this condition is satisfied if and only if

$$ \frac{v^{(n)+} + e v^{(n)-}}{1 + e} = \text{prox}(\mathbb{R}^{-}, v^{(n)-}) \qquad (7) $$

where by definition $\text{prox}(K, y) = \underset{z \in K}{\argmin} (\|z - y\|^2)$ with $K$ convex. At a gap closure time, $v^{(n)-} > 0$ so the LHS is zero and $v^{(n)+} = -e v^{(n)-}$ and the impact law is retrieved. To summarize, Eq. (3) combines the unilateral contact conditions with the impact law.

Altogether, the nonsmooth thermomechanical model is expressed in a first-order form as:

$$ A \begin{bmatrix} \dot{u} \\ \dot{v} \\ \dot{\theta} \end{bmatrix} = B \begin{bmatrix} u \\ v \\ \theta \end{bmatrix} + \begin{bmatrix} 0 \\ f^u + \lambda e^{(n)} \\ f^\theta + \alpha \lambda e^{(n)} \end{bmatrix} \qquad (8a) $$

$$ -\lambda \in \partial \psi_{\mathcal{T}_{\mathbb{R}^+}(d-u^{(n)})} \left( \frac{v^{(n)+} + e v^{(n)-}}{1 + e} \cdot e^{(n)} \right) \qquad (8b) $$

with

$$ A = \begin{bmatrix} I & 0 & 0 \\ 0 & M & 0 \\ 0 & 0 & C_{\theta\theta} \end{bmatrix}, \quad B = \begin{bmatrix} 0 & I & 0 \\ -K_{uu} & -C_{uu} & -K_{u\theta} \\ 0 & -C_{\theta u} & -K_{\theta\theta} \end{bmatrix}. \qquad (9) $$
## 2. MODEL REDUCTION

Industrial Finite Element Models of turbomachines might include thousands of degrees-of-freedom (dofs), leading to high numerical costs. The objective of model reduction is to approximate the dynamics with a “reasonable” accuracy using a reduced number of dofs. The present section describes a first-order reduction of the thermomechanical model (8).

Various reduction methods are available [7–10]. The Craig-Bampton incarnation [7] is recognized as a powerful technique when dealing with unilateral constraints, as they can be properly embedded in the reduced basis. In short, this reduction method defines the contact nodes as *boundary nodes* and projects the dynamics of the other *interior nodes* on a small family of low-frequency mode shapes.

A variation of the Craig-Bampton reduction technique devoted to thermoelasticity exists [11]. It merges the second-order dynamics and the first-order<fnref n="1" /> heat equation into one second-order equation, resulting in a singular generalized mass matrix. Instead, the present work suggests a similar reduction however directly applied to the first order system (8).

Because of the contact conditions, the description of the boundary is essential, therefore the boundary dofs of the complete model $u^{(n)}$, $v^{(n)}$, $\theta^{(n)}$ are kept as dofs in the reduced basis. This distinction between *boundary nodes* and *interior nodes* is reflected by rearranging the state vector $X$ through the permutation matrix $P$ chosen such that:

$$X = \begin{bmatrix} u \\ v \\ \theta \end{bmatrix} = P \begin{bmatrix} u_b \\ v_b \\ \theta_b \\ u_i \\ v_i \\ \theta_i \end{bmatrix} \quad (10)$$

where the subscripts $\bullet_b$ and $\bullet_i$ denote the boundary and the interior, respectively.

The behaviour of the internal nodes is assumed to be governed by a chosen number of shapes of two kinds:

* *static shapes*, describing the response of the internal nodes to a prescribed unit displacement, velocity or temperature applied to the boundary;

* *dynamic shapes*, describing the modal dynamics of the internal nodes while keeping the boundary fixed in its resting position.

The static shapes are obtained by successively computing the static response $\Psi$ to a unit prescribed displacement, velocity and temperature at the boundary, is successively computed. The equilibrium of the internal nodes yields:

$$[B_{ib} \quad B_{ii}] \begin{bmatrix} I \\ \Psi \end{bmatrix} = [0] \quad (11)$$

where $B_{ib}$ and $B_{ii}$ are the blocks (2, 1) and (2, 2) of $P^{-1}BP$. It follows that $\Psi = -B_{ii}^{-1}B_{ib}$, which is easily computed from $B$ in (9).

Then, the mechanical and thermal mode shapes, computed on the clamped-boundary problem such as they should not alter the behaviour of the boundary dofs, are first calculated assuming the mechanical and thermal problems are uncoupled. The reduction is operated by selecting only a small number of such modes: the $n^{uv} = 2n - 2$ mechanical modes with the eigenvalue of lower modulus are selected; similarly, the $n^{\theta} = n - 1$ thermal modes with eigenvalues of the lower modulus are selected. Though not true in general, one may notice that in a unidimensional problem such as the one presently considered, mechanical and thermal mode shapes are identical. In a second step, following [11], thermomechanical coupling is taken into account through a correction term $\Phi^{u\theta}$. The thermal-induced forces $-K_{ii}^{u\theta}\Phi^{\theta\theta}$ are balanced with internal forces satisfying $K_{ii}^{uu}\Phi^{u\theta} = -K_{ii}^{u\theta}\Phi^{\theta\theta}$ which yields the correction term $\Phi^{u\theta}$. Altogether, the full model is approximated by:

$$\begin{bmatrix} u_b \\ v_b \\ \theta_b \\ u_i \\ v_i \\ \theta_i \end{bmatrix} \approx R \begin{bmatrix} u_b \\ v_b \\ \theta_b \\ q^{uv} \\ q^{\theta} \\ n^{uv} \end{bmatrix} \quad \begin{array}{l} \} 1 \\ \} 1 \\ \} 1 \\ \} n^{uv} \\ \} n^{\theta} \\ \end{array} \quad := R\bar{X} \quad (12)$$

where the $3n \times (3 + n^{uv} + n^{\theta})$ reduction matrix $R$ is given by

$$R = \begin{bmatrix} I_3 & 0 \\ \Psi & \Phi \end{bmatrix} \quad \text{with} \quad \Phi = \begin{bmatrix} \Phi^{uu} & \Phi^{u\theta} \\ 0 & \Phi^{\theta\theta} \end{bmatrix}. \quad (13)$$

Instead of solving Eq. (8a) of the form $A\dot{X} = BX + F$, the goal is to solve $APR\dot{X} = BPR\bar{X} + F + r$ where $r$ is the residual term stemming from the approximation. This equation is then projected in a vector subspace orthogonal to the residual, by multiplying to the left with a matrix $\tilde{R}^T P^T$. While in the standard mechanical reduction process, the matrices $A$ and $B$ are both symmetric and the projection and basis transfer matrices $R$ and $\tilde{R}^T$ are equal, this no longer holds in the thermoelastic case. The matrix $\tilde{R}^T$ is chosen following the empirical formula proposed in [11, Eq. (15)]. The projection yields $\tilde{R}^T P^T APR\dot{X} \approx \tilde{R}^T P^T BPR\bar{X} + \tilde{R}^T P^T F$, that is:

$$\bar{A}\dot{X} = \bar{B}\bar{X} + \bar{F} \quad (14)$$

with $\bar{A} := \tilde{R}^T P^T APR$, $\bar{B} := \tilde{R}^T P^T BPR$, $\bar{F} := \tilde{R}^T P^T F$. Both $\bar{A}$ and $\bar{B}$ are square matrices with $3 + n^{uv} + n^{\theta}$ rows. Finally, the reduced-order problem reads

$$\begin{cases} \begin{bmatrix} \dot{u}_b \\ \dot{v}_b \\ \dot{\theta}_b \\ \dot{q}^{uv} \\ \dot{q}^{\theta} \end{bmatrix} = \bar{B} \begin{bmatrix} u_b \\ v_b \\ \theta_b \\ q^{uv} \\ q^{\theta} \end{bmatrix} + \begin{bmatrix} 0 \\ f^u + \lambda \\ f^\theta + \alpha\lambda \\ 0 \\ 0 \end{bmatrix} \\ \lambda \in \partial\psi_{T_{\mathbb{R}^+}(d-u_b)} \left( \frac{v_b^+ + e v_b^-}{1 + e} \right). \end{cases} \quad (15a) \quad (15b)$$

For implementation purposes, it may be convenient to split arbitrarily $q^{uv}$ and rearrange the equations in a form similar

<footnote n="1">¹The heat equation being a parabolic PDE, its space discretization using the Finite Element Method yields a first order time dependent ODE.</footnote>
to Eq. (8). The same algorithm can then be used to solve the full and the reduced models. However, this may lead to nonzero blocks $A_{21}$ and $A_{23}$.

### 3. TIME DISCRETIZATION OF THE MDI

We now proceed with the time discretization of Measure Differential Inclusions of the form (8). The $\theta$-method scheme is chosen for the time discretization, with an integration parameter $\gamma = 0.5$. This scheme is implicit, avoiding stability issues and numerical stiffness due to thermomechanical coupling, and less dissipative than a simple Euler implicit scheme. Introducing the time step $h > 0$, the time domain is discretized with $t_k = kh$ such that $X_k \approx X(t_k)$ with $X_k = [u_k^\top v_k^\top \theta_k^\top]^\top$ and $f_k \approx f(t_k)$ with $f_k = [0 f_k^{u\top} f_k^{\theta\top}]^\top$. We also introduce the vectors $C_{N\theta}$ such that $C_N^\top X = v^{(n)}$ and $C_{N\theta} = [0 e^{(n)\top} \alpha e^{(n)\top}]^\top$. The time discretization of Eq. (8) then reads:

$$ \frac{1}{h} A(X_{k+1} - X_k) = \gamma(BX_{k+1} + f_{k+1}) + (1 - \gamma)(BX_k + f_k) + C_{N\theta}\lambda_{k+1} \quad (16a) $$

and

$$ -\lambda_{k+1} \in \partial\psi_{\mathcal{T}_{\mathbb{R}^+}(d-u_{k+1}^{(n)})} \left( \frac{v_{k+1}^{(n)} + ev_k^{(n)}}{1 + e} \right). \quad (16b) $$

This can be reformulated under the form

$$ X_{k+1} - X_{k+1,p} = D^{-1}C_{N\theta}\lambda_{k+1} \quad (17) $$

where $X_{k+1,p} = X_k + D^{-1}h[BX_k + \gamma f_{k+1} + (1 - \gamma)f_k]$ is the predicted state vector neglecting contact constraints, and $D = A - \gamma hB$ is the state iteration matrix [3] supposedly invertible.

Left-multiplying equation (17) by $C_N^\top$ and solving for $\lambda_{k+1}$ yields

$$ \lambda_{k+1} = [C_N^\top D^{-1} C_{N\theta}]^{-1} (v_{k+1}^{(n)} - v_{k+1,p}^{(n)}). \quad (18) $$

Since (16b) is a cone inclusion, $\lambda$ can be multiplied by the scalar $C_N^\top D^{-1} C_{N\theta}/(1 + e)$ without changing the RHS. Using this property and the two following equalities:

$$ \frac{1}{1 + e} (v_{k+1}^{(n)} - v_k^{(n)}) = \frac{v_{k+1}^{(n)} + ev_k^{(n)}}{1 + e} - v_k^{(n)}, \quad (19) $$

$$ v_{k+1,p}^{(n)} = v_k^{(n)} + C_N^\top D^{-1} h [BX_k + \gamma f_{k+1} + (1 - \gamma)f_k] \quad (20) $$

together with Eq. (18), it finally comes that inclusion (16b) can be written as

$$ \frac{v_{k+1}^{(n)} + ev_k^{(n)}}{1 + e} - b_k \in -\partial\psi_{\mathcal{T}_{\mathbb{R}^+}(d-u_{k+1}^{(n)})} \left( \frac{v_{k+1}^{(n)} + ev_k^{(n)}}{1 + e} \right) \quad (21) $$

for some $b_k$. Using the property of the prox operator that $x-y \in -\partial\psi_K(x) \iff x = \text{prox}(K, y)$ [3, Eq. (1.36)], Eq. (21) becomes

$$ \frac{v_{k+1}^{(n)} + ev_k^{(n)}}{1 + e} = \text{prox}(\mathcal{T}_{\mathbb{R}^+}(d - u_{k+1}^{(n)}), -b_k). \quad (22) $$

The prox operator can be made explicit, using its definition and Eq. (4):

$$ \operatorname{prox}(\mathcal{T}_{\mathbb{R}^+}(u), v) = \begin{cases} v & \text{if } u > 0 \quad (23) \\ \min(0, v) & \text{if } u = 0. \quad (24) \end{cases} $$

Problem (8) can be solved iteratively from:

$$ \left\{ \begin{array}{l} \text{Eq. (16a) except row } 2n \hfill (25a) \\ v_{k+1}^{(n)} = -ev_k^{(n)} + (1+e)\text{prox}(\mathcal{T}_{\mathbb{R}^+}(d-u_{k+1}^{(n)}), -b_k) \quad (25b) \\ \lambda_{k+1} \text{ from Eq. (18)} \hfill (25c) \end{array} \right. $$

yielding $3n+1$ unknowns ($u_{k+1}$, $v_{k+1}$, $\theta_{k+1}$, $\lambda_{k+1}$) for $3n+1$ equations. The normal force $\lambda_{k+1}$ can be eliminated from (25a) using (25c), and the remaining can be solved in a semi-smooth Newton loop, at each time step.

### 4. RESULTS

Model parameters are tuned to match a typical blade of compressor<fnref n="2" /> using $n = 20$. The contacting node is subjected to a periodic sinusoidal loading of amplitude 7 N, rotation speed: 2000 rpm (rotations per minute) representing a cyclic loading. The clearance at rest is chosen as 15 µm. Results are presented with attention paid to the effect of thermal coupling, on the model reduction and on the robustness of the proposed model and simulation methods. Unless otherwise specified, the time step is $h = 10^{-4}$ s.

#### 4.1 Effect of thermal coupling

The blade tip position, the temperature and the normal contact force are depicted in Fig. 2 for two models: the one described in Sec. 1 with $K^{u\theta} = \beta K^{uu}$ and $\beta = -4.5 \times 10^{-7} \text{ m K}^{-1}$, and its uncoupled version obtained by setting $K^{\theta u} = 0$ in Eq. (2a). In the uncoupled model, the positions and the contact force are periodic in time; the temperature increases at each contact due to $f_{tip}$ and does not have time to diffuse at the considered time scale—the characteristic time of diffusion is $(l/n)^2 \rho c / \lambda_{th} \sim 3$ min. In the coupled model, the heat increase generates an expansion of the blade, which tends to close the gap and increase the magnitude of the normal contact force, increasing the heat flux, and so on: the coupling leads to an instability. In a more realistic model, the gap would open because of blade deflection for some level of contact force; this is out of the scope of the present work.

The heat flux $f_{tip} = \alpha\lambda$ at the blade tip reproduces the heat created in a turbomachine when a blade interacts with the surrounding casing. The coefficient $\alpha$ accounts for the friction coefficient as well as the tangential contact velocity; in the case of a Coulomb friction model, $\alpha$ is the product of Coulomb's friction coefficient and the tangential velocity. Three simulations are depicted in Fig. 3 with various

<footnote n="2">²Cross section $A = 10^{-6} \text{ m}^2$ – length $L = 30 \times 10^{-2} \text{ m}$ – Young’s modulus $E = 110 \times 10^9 \text{ Pa}$ – mass density: 4500 kg m$^{-3}$ – specific heat capacity $c = 585 \text{ J kg}^{-1} \text{ K}^{-1}$ – thermal conductivity $\lambda_{th} = 3.35 \text{ W m}^{-1} \text{ K}^{-1}$. Total mass $m = 0.1 \text{ kg}$.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Thermal coupling: coupled [—] and uncoupled [—] models.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Influence of the coefficient $\alpha$ in the heat flux $f_{\text{tip}} = \alpha\lambda$ created during contact phases. [—] $\alpha = 1 \, \text{m s}^{-1}$. [—] $\alpha = 3 \, \text{m s}^{-1}$. [—] $\alpha = 7 \, \text{m s}^{-1}$.</figcaption></figure>

$\alpha$. The above-described thermal instability is observed with $\alpha = 7 \, \text{m s}^{-1}$: after 0.3 s, as a result of the thermomechanical coupling, the contact remains closed while the normal

contact force and the temperature diverge, as expected. It is worth mentioning this physical instability induces no numerical issue or artefact.

Overall, the thermomechanical coupling is shown to have a significant effect on the mechanical response of the model.

## 4.2 Model reduction

We now proceed with the comparison between the complete coupled model and two reduced models, following the procedure described in Sec. 2. The first reduced model (RM1) has $n^{\text{uv}} = 3$ mechanical shapes and $n^{\theta} = 12$ thermal shapes. The second reduced model (RM2) is even smaller in size with $n^{\text{uv}} = 2$ and $n^{\theta} = 5$. It is recalled that the full model corresponds to $n^{\text{uv}} = 2n - 2 = 38$ and $n^{\theta} = n - 1 = 19$, so the size of the problem to solve is divided by approximately 3 for RM1 and 6 for RM2. It is worth recalling that during the reduction, some coupling terms appear in the matrix $A$, making it hard to foresee how the computation times will be reduced.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison between the complete model and the two reduced models. [—] Complete model. [—] RM1. [—] RM2.</figcaption></figure>

The simulations using the three models are depicted in Fig. 4. As long as the increase of temperature is limited, the three models give similar positions and contact forces, meaning that the dynamics is governed by a very limited number of degrees-of-freedom, typically three here. After a few contact phases, RM1 estimates the temperature pretty accurately while RM2 significantly underestimates the increase of temperature, leading to an underestimation of the contact force. This indicates that about a dozen of thermal modes is needed to accurately describe the solution of the problem. The number of thermal modes to ensure convergence exceeds the
number of mechanical modes: this is presumably due to the
rapid changes of temperature compared to the characteristic
time of diffusion.

All in all it has been shown that the reduction technique
proposed in Sec. 2 is compatible with the nonsmooth ther-
momechanical contact problem developed in Sec. 1.

4.3 Sensitivity to numerical parameters

In Carpenter’s scheme [4], contact is dealt with by an explicit
numerical procedure inconveniently stiff for the heat equa-
tion. Small time steps are required for stability purposes and
with the set of parameters of the present work, the stability
condition corresponds approximately to h < 2 × 10⁻⁵ s. For
such time steps, results similar to the ones generated with the
nonsmooth model are observed, see Fig. 5. A close-up view of

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison between Carpenter’s scheme and the proposed scheme. The three curves cannot be distinguished. [•••] Carpenter with $h = 10^{-5}$ s. [—] Proposed scheme with $h = 10^{-5}$ s. [—] Proposed scheme with $h = 10^{-6}$ s.</figcaption></figure>

the first contact phase is provided in Fig. 6. At these space and
time scales, the proposed scheme with $h = 10^{-6}$ s or $10^{-5}$ s
is still very similar to Carpenter’s scheme with $h = 10^{-5}$ s.
Significant differences become visible with $h = 10^{-4}$ s, how-
ever the predicted level of force and positions remain correct,
and these small-scale discrepancies does not significantly
alter the long-term behaviour of the forced system as seen
in Fig. 5. Another difference is visible in the position of the
contacting node. Carpenter’s scheme consists in projecting
the penetrating blade position $u^{(n)}$ on the contact interface:
this step might add superfluous potential energy. In contrast,
the nonsmooth model does not prevent penetration. This
is particularly visible in Fig 6 with curve corresponding to

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Close-up of the position of the blade tip and the contact force during a contact phase for two different time steps. [—] Carpenter's scheme with $h = 10^{-5}$ s. [•••] Proposed scheme with $h = 10^{-5}$ s. [••••] Proposed scheme with $h = 10^{-6}$ s. [•••••] Proposed scheme with $h = 10^{-4}$ s. The first three curves are nearly indistinguishable.</figcaption></figure>

the largest time step. When penetration occurs, contact is
maintained at the exact same level as long as the normal
force is positive. This residual penetration directly depends
on *h* and *ū*<sup>(n)</sup>, and can thus be easily controlled: it will never
exceed *h*ū<sup>(n)</sup>. As soon as unilateral contact is activated, the
normal contact force becomes non-zero, see the first time
step after 4 ms, magenta curve. The simulation with the
smallest time step shows two Delta-like spikes near 4 and
5 ms: the height of the spikes multiplied by the time step
corresponds to the momentum of the contacting mass just
before impact (m<sub>n</sub>Δ*ū*<sup>(n)</sup> = *λ*Δ*t*); the impulsive feature of *λ*
appears as Δ*t* decreases to 0, illustrating the non-smooth
character of contact.

The nonsmooth simulations rely on an implicit scheme,
unconditionally stable in the linear case. There is no regu-
larization parameter and the nonsmooth character can be
observed even with very coarse time steps, see Fig. 6, top. De-
spite the large time step (the first contact phase is covered by
about 65 steps with $h = 10^{-4}$ s while around 650 for Carpen-
ter), the position of the last node and the normal contact force
are estimated with reasonable accuracy. On a larger scale, as
shown in Fig. 7, using a large time step ($h = 5 \times 10^{-4}$ s leads
to a slight underestimation of the contact force at the contact
interface. Still, the position is accurately described and the
qualitative behaviour is recovered.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Influence of the time step in the proposed formulation. [—] $h = 10^{-5}$ s. [···] $h = 10^{-4}$ s. [---] $h = 5 \times 10^{-4}$ s. The first two curves are indistinguishable.</figcaption></figure>

## 5. CONCLUSION

Carpenter's scheme can be used to simulate rotor-stator interactions. However, it relies on an explicit numerical scheme which is not suitable for thermomechanical models, as it leads to constrained numerical stability. In this method, complementarity conditions are dealt with by projecting any predicted penetration on the contact interface. This strategy hides an impact law and can lead to an increase of the total energy. In this work, a sounder basis is proposed for the investigation of thermomechanical effects in rotor-stator contact interaction, through the use of nonsmooth methods. A simple 20-dof blade impacting a fixed carter is considered and modelled using a Measure Differential Inclusion accounting for thermomechanical coupling and friction-induced heat fluxes at the blade tip. The numerical methods used to simulate the model are detailed.

As expected, thermal effects can significantly modify the dynamics of the blade. While Carpenter's scheme requires time steps smaller than $2 \times 10^{-5}$ s to achieve numerical stability, the proposed nonsmooth solver is robust and allows to perform faster simulations thanks to larger time steps. No significant difference in the prediction of the contact force, position and temperature was observed with time steps of $5 \times 10^{-4}$ s. For large time steps (approx. $5 \times 10^{-3}$ s and greater), the temperature and the normal contact force is underestimated. This could be corrected by increasing artificially the friction-induced heat flux to avoid the error in temperature propagating to displacements on longer time

scales.

Industrial turbomachine models include thousands of nodes, requiring model reduction. In this context, a first-order Craig-Bampton-like reduction is adapted to linear thermomelastic systems. While the dynamics can be described with a very small number of mechanical shapes (two or three), due to the rapid changes of temperature induced by contact, the number of thermal shapes cannot be reduced to the same degree. In the end, in the presented results consisting in eight successive contact phases, the positions, temperatures and normal force were well estimated using about a third of the full model size.

Altogether, the proposed methodology offers a sound ground for nonsmooth thermomechanical simulations. Future works include the simulation of more challenging models with intricate contact interfaces and the incorporation of blade deflection and damage mechanisms to reflect the presence of abradable coating on the casing.

## REFERENCES

[1] Marie-Océane Parent and Fabrice Thouverez. Phenomenological model for stability analysis of bladed rotor-to-stator contacts. In *16th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery*, Honolulu, 2016. [hal-01537643].

[2] Zhiyong Shang, Jun Jiang, and Ling Hong. The global responses characteristics of a rotor/stator rubbing system with dry friction effects. *Journal of Sound and Vibration*, 330(10):2150-2160, 2011. [hal-01537669].

[3] Vincent Acary and Bernard Brogliato. *Numerical methods for nonsmooth dynamical systems: Applications in mechanics and electronics*. Springer Science & Business Media, 2008.

[4] Nicholas Carpenter, Robert Taylor, and Michael Katona. Lagrange constraints for transient finite element surface contact. *International journal for numerical methods in engineering*, 32(1):103-128, 1991. [hal-01389918].

[5] Lawrence Shampine and Charles Gear. A user's view of solving stiff ordinary differential equations. *SIAM Review*, 21(1):1-17, 1979.

[6] Tyrrell Rockafellar. Clarke's tangent cones and the boundaries of closed sets in $\mathbb{R}^n$. *Nonlinear Analysis: Theory, Methods & Applications*, 3(1):145-154, 1979.

[7] Roy Craig and Mervyn Bampton. Coupling of substructures for dynamic analyses. *AIAA Journal*, pages 1313-1319, 1968. [hal-01537654].

[8] Richard MacNeal. A hybrid method of component mode synthesis. *Computers & Structures*, 1(4):581-601, 1971. [hal-01537661].

[9] Matthew Castanier, Yung-Chang Tan, and Christophe Pierre. Characteristic constraint modes for component mode synthesis. *AIAA Journal*, 39(6):1182-1187, 2001. [hal-01537673].
[10] Duc-Minh Tran. Component mode synthesis methods using partial interface modes: Application to tuned and mistuned structures with cyclic symmetry. *Computers & Structures*, 87(17-18):1141-1153, 2009. [hal-01537665].

[11] Philippe Nachtergaele, Daniel Rixen, and Alexander Steenhoek. Efficient weakly coupled projection basis for the reduction of thermo-mechanical models. *Journal of computational and applied mathematics*, 234(7):2272-2278, 2010. [hal-01517953].