Thermomechanical transient dynamics of radial
rotor-stator contact

Coline Jacobs, Fabrice Thouverez, Mathias Legrand, Patricio Almeida

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

To cite this version:

Coline Jacobs, Fabrice Thouverez, Mathias Legrand, Patricio Almeida. Thermomechanical transient dynamics of radial rotor-stator contact. 2024. hal-04558059

HAL Id: hal-04558059

https://hal.science/hal-04558059

Preprint submitted on 24 Apr 2024

**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.
# Thermomechanical transient dynamics of radial rotor-stator contact

Coline Jacobs, Fabrice Thouverez, Mathias Legrand, Patricio Almeida

► To cite this version:

Coline Jacobs, Fabrice Thouverez, Mathias Legrand, Patricio Almeida. Thermomechanical transient dynamics of radial rotor-stator contact. Journal of Sound and Vibration, 2025, 10.1016/j.jsv.2025.119207. hal-04558059v2

HAL Id: hal-04558059

https://hal.science/hal-04558059v2

Submitted on 8 Jul 2025

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 transient dynamics of radial rotor-stator contact

Coline Jacobs<sup>a,b,c</sup>, Fabrice Thouverez<sup>a</sup>, Mathias Legrand<sup>b</sup>, Patricio Almeida<sup>c</sup>

<sup>a</sup>Laboratoire de Tribologie et Dynamique des Systèmes, Ecole Centrale de Lyon, 36 Avenue Guy de Collongue, Ecully, France

<sup>b</sup>Structural Dynamics and Vibration Laboratory, McGill University, 817 Sherbrooke West, Montreal, Quebec, Canada

<sup>c</sup>Safran Helicopter Engines, Avenue Joseph Szydlowski, Bordes, France

## Abstract

The present paper discusses the numerical investigation of structural interactions initiated radially between a turbine shaft and an annular stator. The rotorshaft is modeled under the rigid bodies assumption while the stator is discretized in space with curved Bernoulli beam elements. The turbomachine is assumed to be driven by two torques: the first one is prescribed from aerodynamics while the second one is induced by sliding friction from Coulomb's model and must be assessed. The study is motivated by the need to predict the dynamic behavior of the rotorshaft in the presence of radial contact between the labyrinth and the annular stator. If rotordynamics accounting for unilateral and frictional contact is widely developed in literature, the combination with an unknown rotational velocity and thermomechanical coupling must still be investigated. Indeed, heat generation and subsequent thermal expansion is expected to modify the contact interface and the dynamic response. The system is solved in the time domain and contact treatment is achieved though a modified version of the Carpenter algorithm. The identification of the conditions affecting the rotational speed is achieved through a sensitivity analysis on the stator properties, friction coefficient and coefficient of thermal expansion. Results confirm that the higher the friction coefficient and stator stiffness are, the lower the maximum rotational velocity value is. It is also shown that high temperatures are located in specific regions on the stator and that the temperature peaks are amplified with thermal expansion.

**Keywords:** time-marching techniques, rotor-stator interaction, torque, thermomechanical coupling, labyrinth, Coulomb's friction, transient dynamics, thermal expansion

## 1. Introduction and background material

Helicopter engines are powered by gas turbines that perform classical thermodynamic cycles of compression, combustion and expansion. The latter is achieved through the blading of possibly multiple turbine stages. Expansion is characterized by a conversion of kinetic energy of the burnt gas into mechanical energy at the rotorshaft level. In nominal operating conditions, this mechanical energy is used to drive the compressor stages and maintain the cycle. The present paper considers that the turbine is in standalone operation and its dynamics is governed by two torques. The first one drives the turbine and comes from aerodynamics and is assumed to be constant. The second torque is resistive and stems from frictional contact occurrences: it is an unknown of the problem. The rotor-stator interaction of interest occurs during a radial contact between the labyrinth and its sealing envelope. The labyrinth is made of a succession of several annular teeth as illustrated in Figure 1. Their purpose is to make the air difficult to flow. Because of turbine vibration, the gap of several hundreds of microns is consumed so that unilateral contact and frictional occurrences are initiated radially between the teeth and the seal envelope. For the sake of simplicity, a single tooth is considered, called *rotor* in the remainder, while the envelope is called *stator*. Rubbing within the seal will generate contact efforts, a braking torque and thus heat. The objective of the paper is to investigate how friction affects the rotational velocity of the rotor and the temperature distribution in the stator.

Experiments and models provided in literature show that a rotor confined within a stator might generate rubbing mechanisms and might exhibit various and complex dynamical responses [42, 34, 15, 47]. Commonly, the motion of the rotor is described by (1) the lateral vibration of the structure (characterized through precessional speed $v_{prec}$) and (2) the rotation with angular velocity $\Omega$ [36]. When the structure is orbiting in the same direction as $\Omega$, the response is called *forward precessional* motion. In other words, the circumferential, or precessional, velocity is realized in the direction of $\Omega$. In contrast, when the precessional velocity acts against $\Omega$, *backward precessional* takes place [30]. When contact is fully annular, backward motion is divided into two subcategories as follows.
# Thermomechanical transient dynamics of radial rotor-stator contact

Coline Jacobs$^{a,b,c}$, Fabrice Thouverez$^{a}$, Mathias Legrand$^{b}$, Patricio Almeida$^{c}$

$^{a}$Laboratoire de Tribologie et Dynamique des Systèmes, Ecole Centrale de Lyon, Ecully, France

$^{b}$Structural Dynamics and Vibration Laboratory, McGill University, Montreal, Quebec, Canada

$^{c}$Safran Helicopter Engines, Bordes, France

## Abstract

The present paper discusses the numerical investigation of structural interactions initiated radially between a turbine shaft and an annular stator. The rotorshaft is modeled under the rigid bodies assumption while the stator is discretized in space with curved Bernoulli beam elements. The turbomachine is assumed to be driven by two torques: the first one is prescribed from aerodynamics while the second one is induced by sliding friction from Coulomb's model and must be assessed. The study is motivated by the need to predict the dynamic behavior of the rotorshaft in the presence of radial contact between the labyrinth and the annular stator. If rotordynamics accounting for unilateral and frictional contact is widely developed in literature, the combination with an unknown rotational velocity and thermomechanical coupling must still be investigated. Indeed, heat generation and subsequent thermal expansion is expected to modify the contact interface and the dynamic response. The system is solved in the time domain and contact treatment is achieved though a modified version of the Carpenter algorithm. The identification of the conditions affecting the rotational speed is achieved through a sensitivity analysis on the stator properties, friction coefficient and coefficient of thermal expansion. Results confirm that the higher the friction coefficient and stator stiffness are, the lower the maximum rotational velocity value is. It is also shown that high temperatures are located in specific regions on the stator and that the temperature peaks are amplified with thermal expansion.

**Keywords:** time-marching techniques, rotor-stator interaction, torque, thermomechanical coupling, labyrinth, Coulomb's friction, transient dynamics, thermal expansion

## 1. Introduction and background material

Helicopter engines are powered by gas turbines that perform classical thermodynamic cycles of compression, combustion and expansion. The latter is achieved through the blading of possibly multiple turbine stages. Expansion is characterized by a conversion of kinetic energy of the burnt gas into mechanical energy at the rotorshaft level. In nominal operating conditions, this mechanical energy is used to drive the compressor stages and maintain the cycle. The present paper investigates the accidental situation when the shaft fails such that the turbine is no longer coupled to its tare load, i.e., the compressors. This scenario is also known as *loss of load* in the literature [50]. Since burnt gas is still flowing through the blading, the turbine undergoes an overspeed transient and, if not managed properly, may blow up and release high-energy debris [5]. Recent works have demonstrated that turbine speed limitation is possible through rotor-stator contacts occurring axially [43, 22, 50, 28], or radially [26]. However, little information is available regarding radial rub during a shaft failure. As noted above, contact occurrences in aircraft engines, which are the most explored, generally have an axial nature and occur at the blading level. In helicopter engines, the rotor-stator interaction of interest occurs during a radial contact between the labyrinth and its sealing envelope. In practice, the labyrinth is made of a succession of several annular teeth as illustrated in Figure 1. Their purpose is to make the air difficult to flow. Because of turbine vibration, the gap of several hundreds of microns is consumed so that unilateral contact and frictional occurrences are initiated radially between the teeth and the seal envelope. Figure 1 aims to provide a brief overview of a labyrinth seal. Research on rotor-stator contact within this component is limited, as it is often overlooked in aircraft engines during shaft failure scenarios. However, in helicopter engines, contact in this area plays a critical role in preventing overspeed. Achieving initial insights requires certain simplifications to the seal model. A single tooth is considered, called *rotor* in the remainder. This paper does not seek to implement a detailed model of the labyrinth seal; rather, its goal is to precisely define the contact region on the fractured shaft and explore the phenomenology of radial rubbing. While future research should address multiple
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Example of a labyrinth seal located at the shaft level [49]</figcaption></figure>

• The term *dry whip* is used when the rotor rolls and slips, i.e. the sliding speed $v_T$ is non-zero. This response is commonly to be avoided in rotordynamics because of its instability, coming with high contact efforts that may damage the structures in contact [27, 37].

• The term *dry whirl* is used when rolling without sliding motion occurs. During this response, the kinematic constraint $||v_T|| = 0$, must be satisfied where $v_T$ stands for the relative tangential velocity at contact level.

In practice, fully annular rub is not always taking place. A third response type is *partial contact*, characterized by a non-circular trajectory of the rotor. The response with partial contact can take place in both forward and backward precessional motions [3]. Actually, based on a former model developed by the authors [26], it is expected that this type of transition with contact separation takes place.

Many parameters are known to govern the type and switches of precession. For instance, previous works showed that increased friction coefficients promote backward precessional motion [11, 29, 13]. Friction versus rotational speed diagrams suggested in [29] could show the existence and co-existence domains for precession types. A similar diagram where the rotor-to-stator stiffness ratio is considered was also suggested. It turned out that, for a sufficiently high stiffness ratio (e.g. a softer stator), dry whip could not be triggered by rotor unbalance. Various authors also investigated the transition between dry whip and dry whirl. Black [7] showed analytically that there exist specific rotational speed ranges for dry whirl to occur. These results were confirmed numerically and experimentally on a Jeffcott rotor in [50], then extended to multiple degrees-of-freedom in [12]. The authors showed that the precessional frequency is proportional to the rotational speed with a factor equal to the radius to initial gap ratio. Dry whip is characterized by constant precessional frequency and transition between dry whip and dry whirl comes with a frequency jump. The rotational velocity at which the jump happens depends on whether the rotational velocity is increasing or decreasing. This hysteresis phenomenon was also highlighted in [13]. Bartha [5] identified a critical radial velocity function of the impact incidence angle. Above this critical velocity, backward precessional motion occurs. Literature is extensive on radial frictional contact occurrences caused by unbalance. Ma & al. [34] showed that out-of-phase unbalance on two-disc rotors leads to a quasi-periodic response with partial contacts. Instead, in-phase unbalance yields forward full annular rub. Zhang [53] studied periodicity and stability of a turbine shaft when coupled frictional contacts at labyrinth seal and blade tip levels take place initiated by unbalance on the blisk.

Based on the conclusions of a previous work by the authors [26], large contact efforts are expected. When combined with high rotational velocity, unilateral contact with rubbing leads to high frictional power and temperature elevation. Friction creates an energy loss in the system mainly in the form of heat which induces thermal dilation, wear mechanisms and possibly plastic deformation in the components [39, 6, 44]. They all modify the clearance between the contacting bodies. Experimental investigation of wear is performed for a labyrinth-abradable contact with a single tooth in [16]. It is thus clear that the coupling between the structural dynamics and temperature field has a significant influence on the contact interface geometry. Solving the complete multiphysics problem on the contact interface is not a trivial task and should be simplified for a preliminary investigation such as the present work, which is grounded on a thermoelastic constitutive law and does not account for wear. The heat distribution between the bodies is
teeth, this study serves as an important initial milestone. The present paper proposes a simplified scenario: the turbine works in standalone operation and its dynamics is governed by two torques. The first one drives the turbine and comes from aerodynamics and is assumed to be constant. The second torque is resistive and stems from frictional contact occurrences: it is an unknown of the problem. Rubbing within the seal will generate contact efforts, a braking torque and thus heat. The objective of the paper is to investigate how friction affects the rotational velocity of the rotor and the temperature distribution in the stator, i.e. the seal envelope.

<figure><img src="image_2.png" /><figcaption>**Figure 1:** Example of a labyrinth seal located at the shaft level [54]</figcaption></figure>

Experiments and models provided in literature show that a rotor confined within a stator might generate rubbing mechanisms and might exhibit various and complex dynamical responses [46, 38, 16, 52]. Commonly, the motion of the rotor is described by (1) the lateral vibration of the structure (characterized through precessional speed $v_{prec}$) and (2) the rotation with angular velocity $\Omega$ [40]. When the structure is orbiting in the same direction as $\Omega$, the response is called *forward precessional* motion. In other words, the circumferential, or precessional, velocity is realized in the direction of $\Omega$. In contrast, when the precessional velocity acts against $\Omega$, *backward precessional* takes place [34]. When contact is fully annular, backward motion is divided into two subcategories as follows.

* The term *dry whip* is used when the rotor rolls and slips, i.e. the sliding speed $v_T$ is non-zero. This response is commonly to be avoided in rotordynamics because of its instability, coming with high contact efforts that may damage the structures in contact [31, 41].

* The term *dry whirl* is used when rolling without sliding motion occurs. During this response, the kinematic constraint $||v_T|| = 0$, must be satisfied where $v_T$ stands for the relative tangential velocity at contact level.

In practice, fully annular rub is not always taking place. A third response type is *partial contact*, characterized by a non-circular trajectory of the rotor. The response with partial contact can take place in both forward and backward precessional motions [3]. Actually, based on a former model developed by the authors [30], it is expected that this type of transition with contact separation takes place.

Many parameters are known to govern the type and switches of precession. For instance, previous works showed that increased friction coefficients promote backward precessional motion [12, 33, 14]. Friction versus rotational speed diagrams suggested in [33] could show the existence and co-existence domains for precession types. A similar diagram where the rotor-to-stator stiffness ratio is considered was also suggested. It turned out that, for a sufficiently high stiffness ratio (e.g. a softer stator), dry whip could not be triggered by rotor unbalance. Various authors also investigated the transition between dry whip and dry whirl. Black [8] showed analytically that there exist specific rotational speed ranges for dry whirl to occur. These results were confirmed numerically and experimentally on a Jeffcott rotor in [55], then extended to multiple degrees-of-freedom in [13]. The authors showed that the precessional frequency is proportional to the rotational speed with a factor equal to the radius to initial gap ratio. Dry whip is characterized by constant precessional frequency and transition between dry whip and dry whirl comes with a frequency jump. The rotational velocity at which the jump happens depends on whether the rotational velocity is increasing or decreasing. This hysteresis phenomenon was also highlighted in [14]. Bartha [6] identified a critical radial velocity function of the impact incidence angle. Above
also a complex multiphysics mechanism. The thermal resistance depends on friction of course, but also on the surface state since conduction takes place through surface roughness. Convection in the gas pockets trapped between the bodies interfaces and radiation also plays a key role in heat transfer [51]. The vast majority of works simplify the heat source computation by distributing the frictional power based on some arbitrary heat transfer ratio [17, 47] or based on the conductance of the bodies [25, 4]. In the present paper, it is assumed that heat is equally distributed between the contacting bodies [25, 17]. The temperature transient is of primary interest for the stator only and thermodynamics mechanisms are ignored for the rotor.

The system of interest and its model are presented in Section 2. The next section introduces the modified version of the Carpenter algorithm selected to solve the nonlinear transient thermomechanics in time domain. Eventually, Section 4 analyzes results for a reference case. Sensitivity studies are performed to identify how friction coefficient, stator properties and thermal expansion affect the response for both physics. Conclusions and prospects are proposed in Section 5.

## 2. System of interest and modeling

### 2.1. Geometry and components

As described in Figure 2, the seal is located upstream the turbine disc. The rotor-stator contact occurs

<figure><img src="image_3.png" /><figcaption>**Figure 2:** Rotor-stator system between static assembly and the turbine (□) supported by a rear bearing (■). Simplification of the assembly ring (○), component A (■) and component B (■) into a ring supported by equivalent springs located at screw points.</figcaption></figure>

radially between a tooth, indicated on the schematic, and the envelope of the surrounding labyrinth seal. In practice, the stator is an annular ring screwed to other components of the engine, hereafter named A and B. The number of screws is defined by $N_{sc}$. In order to account properly for the structures supporting the ring, the static assembly is simplified to a point mass with equivalent springs connected at screw points. Their stiffness is provided by load case simulations performed on the real structures in ANSYS.

Component B, see Figure 2, is modeled by as many stiffness pairs ($k_{Bu}, k_{Bv}$) (blue pairs in Figure 2) as the number of screws ($N_{sc} = 8$ for illustration purposes). Component A contributes to the dynamics by adding a point mass $m_A$ connected at the screws with a radial stiffness $k_A$ (in red). In Figure 2, the point mass is centred in the stator while in reality, the component is located outside the stator. This modification is made for easier springs visualization. The suggested simplifications translate a 3D problem into a 2D model.

### 2.2. Turbine

As illustrated in Figure 3, the geometry of the turbine is also simplified to a rotorshaft including the turbine and the rotor discs. The latter is located on the shaft line at point R. In the present document, the turbine is considered as a perfectly rigid body. The bearing is reduced to a linear and symmetric spring-damper system connected to point P on the rotorshaft. At rest, the shaftline coincides with the Z-axis of the
this critical velocity, backward precessional motion occurs. Literature is extensive on radial frictional
contact occurrences caused by unbalance. Ma & al. [38] showed that out-of-phase unbalance on two-disc
rotors leads to a quasi-periodic response with partial contacts. Instead, in-phase unbalance yields forward
full annular rub. Zhang [58] studied periodicity and stability of a turbine shaft when coupled frictional
contacts at labyrinth seal and blade tip levels take place initiated by unbalance on the blisk.

Based on the conclusions of a previous work by the authors [30], large contact efforts are expected. When combined with high rotational velocity, unilateral contact with rubbing leads to high frictional power and temperature elevation. Friction creates an energy loss in the system mainly in the form of heat which induces thermal dilation, wear mechanisms and possibly plastic deformation in the components [43, 7, 48]. They all modify the clearance between the contacting bodies. Experimental investigation of wear is performed for a labyrinth-abradable contact with a single tooth in [17]. It is thus clear that the coupling between the structural dynamics and temperature field has a significant influence on the contact interface geometry. Solving the complete multiphysics problem on the contact interface is not a trivial task and should be simplified for a preliminary investigation such as the present work, which is grounded on a thermoelastic constitutive law and does not account for wear. The heat distribution between the bodies is also a complex multiphysics mechanism. The thermal resistance depends on friction of course, but also on the surface state since conduction takes place through surface roughness. Convection in the gas pockets trapped between the bodies interfaces and radiation also plays a key role in heat transfer [56]. The vast majority of works simplify the heat source computation by distributing the frictional power based on some arbitrary heat transfer ratio [18, 52] or based on the conductance of the bodies [29, 4]. In the present paper, it is assumed that heat is equally distributed between the contacting bodies [29, 18]. The temperature transient is of primary interest for the stator only and thermodynamics mechanisms are ignored for the rotor.

The system of interest and its model are presented in Section 2. The next section introduces the
modified version of the Carpenter algorithm selected to solve the nonlinear transient thermomechanics in
time domain. Eventually, Section 4 analyzes results for a reference case. Sensitivity studies are performed
to identify how friction coefficient, stator properties and thermal expansion affect the response for both
physics. Conclusions and prospects are proposed in Section 5.

**2. System of interest and modeling**

**2.1. Static assembly and assumptions**

In the previous section, it was noted that, due to shaft break, the compressors are disconnected from
the turbine stage. The turbine rotordynamics thus becomes unbalanced, which has two consequences.
Initially, the turbine’s rotational speed tends to increase, due to the airflow from the combustion chamber.
Additionally, the shaft loses the upstream bearing located at the compressors level. Therefore, the support
conditions of the isolated turbine are modified. The turbine remains supported only by its rear bearing and
by the labyrinth seal where radial rub emerges. As described in Figure 2, the seal is located upstream the
turbine disc. The rotor-stator contact occurs radially between a tooth, indicated on the schematic, and
the envelope of the surrounding labyrinth seal. In practice, the stator is an annular ring screwed to other
components of the engine, hereafter named A and B. The real assembly of components is presented on the
left hand side of Figure 2, which indicates that the screws are aligned with the engine rotation axis. The
present paper suggests an academic version of this system. Since component A is only connected to the
ring by screws, it is simplified as a point mass mA (in red) with equivalent springs of stiffness kA. Note
that in Figure 2 the point mass mA is centered in the stator while the real component is located outside.
This modification is made for easier visualization. The number of screws is defined by Nsc (Nsc = 8 for
illustration purposes). Component B is fixed to the left flange, and is considered to provide stiffness only
to the assembly. It is modeled by Nsc spring pairs (kBu, kBv) (blue pairs in Figure 2). The equivalent
stiffness of components A and B are provided by load case simulations performed on the corresponding
3D FEA models in ANSYS. Radial and tangential efforts are applied on a given screw hole while clamping
neighbor holes. To summarize, the simplified mass and springs models reflect mechanical properties of
components A and B while screw points just refer to locations. The suggested simplifications translate a
3D problem into a 2D model, with a smaller size and thus a more reasonable computation time.

**2.2. Turbine**

As illustrated in Figure 3, the rotating components after shaft failure only account for a standalone turbine
supported by the remaining rear bearing. The geometry is simplified into a rotorshaft including the turbine
<figure><img src="image_4.png" /><figcaption>**Figure 3:** Simplified turbine and corresponding degrees-of-freedom: two pairs of translations at point R (x_r, y_r) and point P (x_p, y_p) and the rotation along shaft line θ</figcaption></figure>

coordinate system frame. The mechanical properties of the turbine reduce to its mass *m*, along with its polar and diametrical inertia *J*<sub>p</sub> and *J*<sub>d</sub>. Its center of gravity is denoted by *G* located at a distance *d* from the rear bearing *P*. As previously said, contact is restricted to the labyrinth tooth (rotor) located along the Z axis with distance *L* from the rear bearing. A common assumption in the literature is to consider the rotation speed as known and constant or at least as a known function of time [30]. Instead, when frictional contact occurrences are involved, the torque is prescribed and the rotational velocity Ω(*t*) becomes an unknown of the problem, which is the case in the present paper and in [41, 15] as well. Because of high rotational speeds in helicopter engines (up to 40 000 rpm), gyroscopic terms cannot be ignored. They are known to couple even further the nonlinear mechanisms in the response [23] and the governing equation involves terms in the stiffness matrix that depend on angular velocity Ω(*t*) and on angular acceleration Ω̇(*t*) [23, 45, 33]. In the present document, since the turbine is perfectly rigid, stiffening effects are ignored and only gyroscopic and inertial effects are included.

Also, ignoring the translational motion along the Z axis by assumption, five degrees-of-freedom should be defined, following [26]. first, the titling angles θ<sub>x</sub>, θ<sub>y</sub> about the X and Y axes, respectively, are assumed to be small and their expression thus reduces to translations at points P and R in the form

$$ \theta_x \approx \frac{y_p - y_r}{L} \quad \text{and} \quad \theta_y \approx \frac{x_r - x_p}{L}. \qquad (1) $$

These translations are defined in the fixed frame (X, Y, Z). The turbine motion is thus fully determined by the vector of generalized coordinates

$$ \mathbf{q}_r = (x_p, y_p, x_r, y_r, \theta)^\top, \qquad (2) $$

where θ is the spinning angle satisfying θ̇(*t*) = Ω(*t*). Displacements of the center of gravity *G* are expressed by linear interpolation of the rotor and bearing displacements, that is *q*<sub>G</sub> = η*q*<sub>r</sub> + (1 - η)*q*<sub>p</sub> for *q* = *x*, *y*, where the length ratio η = *d*/*(L - δ)* is used. A similar definition of the generalized coordinates is found in [15] for the particular case where δ = 0 and η = 1/2. The corresponding governing equation for the turbine dynamics reads

$$ \mathbf{M}_r(\mathbf{q}_r)\ddot{\mathbf{q}}_r + (\dot{\theta}\mathbf{G}_r + \mathbf{D}_b)\dot{\mathbf{q}}_r + \mathbf{K}_b\mathbf{q}_r = \mathbf{f}^{\text{ext}}. \qquad (3) $$

The terms **K**<sub>b</sub> and **D**<sub>b</sub> are the bearing stiffness and damping matrices. The mass and gyroscopic matrices are, respectively

$$ \mathbf{M}_r(\mathbf{q}_r) = \frac{1}{2} \begin{pmatrix} 2M_{11} & 0 & 2M_{12} & 0 & -H(y_r - y_p) \\ 2M_{11} & 0 & 2M_{12} & H(x_r - x_p) & 0 \\ 2M_{22} & 0 & 0 & H(y_r - y_p) & 0 \\ \text{sym.} & 0 & 2M_{22} & -H(x_r - x_p) & 2J_p \end{pmatrix} \quad \text{and} \quad \mathbf{G}_r = \begin{pmatrix} 0 & H & 0 & -H & 0 \\ -H & 0 & H & 0 & 0 \\ 0 & -H & 0 & H & 0 \\ H & 0 & -H & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix} \qquad (4) $$
<figure><img src="image_3.png" /><figcaption>**Figure 2:** Rotor-stator system between static assembly and the turbine (□) supported by a rear bearing (■). Simplification of the assembly ring (□), component A (■) and component B (■) into a ring supported by equivalent springs located at screw points.</figcaption></figure>

and the rotor discs. The latter is located on the shaft line at point R. The rear bearing is reduced to a linear and symmetric spring-damper system connected to point P on the rotorshaft. At rest, the shaft line coincides with the Z-axis of the coordinate system frame. In the present document, the turbine is

<figure><img src="image_4.png" /><figcaption>**Figure 3:** Simplified turbine and corresponding degrees-of-freedom: two pairs of translations at point R (x<sub>r</sub>, y<sub>r</sub>) and point P (x<sub>p</sub>, y<sub>p</sub>) and the rotation along shaft line θ</figcaption></figure>

considered as a perfectly rigid body. The mechanical properties of the turbine reduce to its mass *m*, along with its polar and diametrical inertia *J*<sub>p</sub> and *J*<sub>d*. Its center of gravity is denoted by *G* located at a distance *d* from the rear bearing *P*. As previously said, contact is restricted to the rotor (i.e. labyrinth tooth) located along the *Z* axis with distance *L* from the rear bearing. A common assumption in the literature is to consider the rotation speed as known and constant or at least as a known function of time [34]. Instead, when frictional contact occurrences are involved, the torque is prescribed and the rotational velocity *Ω*(*t*) becomes an unknown of the problem, which is the case in the present paper and in [45, 16] as well. Because of high rotational speeds in helicopter engines (up to 40 000 rpm), gyroscopic terms cannot be ignored. They are known to couple even further the nonlinear mechanisms in the response [25] and the governing equation involves terms in the stiffness matrix that depend on angular velocity *Ω*(*t*) and on angular acceleration *Ω*̇(*t*) [25, 49, 37]. In the present document, since the turbine is perfectly rigid, stiffening effects are ignored and only gyroscopic and inertial effects are included.

Also, ignoring the translational motion along the Z axis by assumption, five degrees-of-freedom should be defined, following [30]. first, the titling angles θ<sub>x</sub>, θ<sub>y</sub> about the X and Y axes, respectively, are
in which $M_{11} = (1 - \eta)^2 m + J_d/(L - \delta)^2$, $M_{12} = \eta(1 - \eta)m - J_d/(L - \delta)^2$, $M_{22} = \eta^2 m + J_d/(L - \delta)^2$, and $H = J_p/(L - \delta)^2$. The external force vector $\mathbf{f}^{\text{ext}} = (0, -(1 - \eta)mg, 0, -\eta mg, \tau_a)^\top$ accounts for gravity and the driving torque from burnt gas blowing through the turbine blading $\tau_a$, which is set constant.

Bearing $P$ is initially at rest but the information on the initial rotational velocity of the rotor are not provided for proprietary reasons. Still, it should be reminded that the first impact takes place in forward precessional direction with a non-vanishing relative tangential velocity.

*Remark.* It is possible to rewrite Equation (3) in order to keep the mass matrix constant and explicitly show the angular acceleration (or Coriolis) stiffness matrix $\mathbf{P}_r = \mathbf{G}_r/2$, which gives

$$ \mathbf{M}_r \ddot{\mathbf{q}}_r + (\mathbf{D}_b + \dot{\theta} \mathbf{G}_r) \dot{\mathbf{q}}_r + (\mathbf{K} + \ddot{\theta} \mathbf{P}_r) \mathbf{q}_r = \mathbf{f}^{\text{ext}} \quad (5) $$

Nevertheless, Equation (3) is preferred since it fits the general formulation of a second order Lagrangian system[2, 1].<fnref n="1" />

## 2.3. Stator

### 2.3.1. Continuum Mechanics and simplifying assumptions

The stator is modeled as a circular ring with constant rectangular cross-sectional area of width $b$ and half-thickness $a$, see Figure 4. The neutral fiber of the ring has a radius $R_s$. The ring is assumed to be thin and satisfies the condition

$$ 2a \ll R_s. \quad (6) $$

The Euler-Bernoulli beam theory is considered. Any point one the ring bears the six degrees-of-freedom ($u, v, w, \theta_s, \theta_r, \theta_w$), see Figure 4. From the small deformation assumption combined with condition (6), it

<figure><img src="image_5.png" /><figcaption>Figure 4: Stator geometry, reference frame and degrees-of-freedom for a curved Bernoulli beam</figcaption></figure>

is possible to simplify the expressions of the deformation in 3D by involving only the three translations ($u, v, w$), the torsion angle ($\theta_s$) and their derivatives with respect to the curvilinear abscissa $s$, denoted with subscripts ,s and ,ss for compactness [18]:

$$ \begin{aligned} \epsilon_{ss} &= \left(\frac{v}{R_s} + u, s\right) - \left(\frac{v}{R_s^2} + v, ss\right) r + \left(\frac{\theta_s}{R_s} - w, ss\right) z \\ \gamma_{rs} &= \frac{z}{R_s} (R_s \theta_{s,s} + w, s) \\ \gamma_{zs} &= \frac{-r}{R_s} (R_s \theta_{s,s} + w, s). \end{aligned} \quad (7) $$

However, in the present paper, contact is only investigated in the $(X, Y)$ plane in Figure 4. As a result, the torsion $\theta_s$ and axial displacement $w$ are neglected. The shear deformations $\gamma_{rs}$ and $\gamma_{zs}$ thus vanish and do not contribute to the deformation energy so that Hooke's law $\sigma_{ss} = E\epsilon_{ss}$ reads [19, 32]

$$ 2E_d = \int_{s_1}^{s_2} ES \left(\frac{v}{R_s} + u, s\right)^2 + EI_z \left(\frac{v}{R_s^2} + v, ss\right)^2 ds \quad (8) $$

<footnote n="1">¹For more information, see also the Siconos user guide.</footnote>
assumed to be small and their expression thus reduces to translations at points P and R in the form

$$ \theta_x \approx \frac{y_p - y_r}{L} \quad \text{and} \quad \theta_y \approx \frac{x_r - x_p}{L}. \qquad (1) $$

These translations are defined in the fixed frame (X, Y, Z). The turbine motion is thus fully determined by the vector of generalized coordinates

$$ \mathbf{q}_r = (x_p, y_p, x_r, y_r, \theta)^{\mathrm{T}}, \qquad (2) $$

where $\theta$ is the spinning angle satisfying $\dot{\theta}(t) = \Omega(t)$. Displacements of the center of gravity G are expressed by linear interpolation of the rotor and bearing displacements, that is $q_G = \eta q_r + (1 - \eta)q_p$ for $q = x, y$, where the length ratio $\eta = d/L$ is used. A similar definition of the generalized coordinates is found in [16] for the particular case where $\eta = 1/2$. The corresponding governing equation for the turbine dynamics reads

$$ \mathbf{M}_r(\mathbf{q}_r)\ddot{\mathbf{q}}_r + (\dot{\theta}\mathbf{G}_r + \mathbf{D}_b)\dot{\mathbf{q}}_r + \mathbf{K}_b\mathbf{q}_r = \mathbf{f}^{\mathrm{ext}}. \qquad (3) $$

The terms $\mathbf{K}_b$ and $\mathbf{D}_b$ are the bearing stiffness and damping matrices. The mass and gyroscopic matrices are, respectively

$$ \mathbf{M}_r(\mathbf{q}_r) = \frac{1}{2} \begin{pmatrix} 2M_{11} & 0 & 2M_{12} & 0 & -H(y_r - y_p) \\ 2M_{11} & 0 & 2M_{12} & H(x_r - x_p) & 0 \\ 2M_{22} & 0 & H(y_r - y_p) & 0 & 0 \\ \text{sym.} & 0 & 2M_{22} & -H(x_r - x_p) & 0 \\ 0 & 0 & 0 & 0 & 2J_p \end{pmatrix} \quad \text{and} \quad \mathbf{G}_r = \begin{pmatrix} 0 & H & 0 & -H & 0 \\ -H & 0 & H & 0 & 0 \\ 0 & -H & 0 & H & 0 \\ H & 0 & -H & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix} \qquad (4) $$

in which $M_{11} = (1 - \eta)^2 m + J_d/L^2$, $M_{12} = \eta(1 - \eta)m - J_d/L^2$, $M_{22} = \eta^2 m + J_d/L^2$, and $H = J_p/L^2$. Note the dependency of mass matrix with components of displacement vector $\mathbf{q}_r$ (cf Equation (2)) reflecting the coupling between transversal displacement and angular acceleration.

$$ \mathbf{f}^{\mathrm{ext}} = (0, -(1 - \eta)mg, 0, -\eta mg, \tau_a)^{\mathrm{T}} \qquad (5) $$

accounts for gravity and the driving torque from burnt gas blowing through the turbine blading $\tau_a$, which is set constant.

Bearing P is initially at rest but the information on the initial rotational velocity of the rotor are not provided for proprietary reasons. Still, it should be reminded that the first impact takes place in forward precessional direction with a non-vanishing relative tangential velocity.

*Remark.* It is possible to rewrite Equation (3) in order to keep the mass matrix constant and explicitly show the angular acceleration (or Coriolis) stiffness matrix $\mathbf{P}_r = \mathbf{G}_r/2$, which gives

$$ \mathbf{M}_r \ddot{\mathbf{q}}_r + (\mathbf{D}_b + \dot{\theta} \mathbf{G}_r) \dot{\mathbf{q}}_r + (\mathbf{K} + \ddot{\theta} \mathbf{P}_r) \mathbf{q}_r = \mathbf{f}^{\mathrm{ext}} \qquad (6) $$

Nevertheless, Equation (3) is preferred since it fits the general formulation of a second order Lagrangian system[1, 2].<fnref n="1" />

## 2.3. Stator

### 2.3.1. Continuum Mechanics and simplifying assumptions

The stator is modeled as a circular ring with constant rectangular cross-sectional area of width b and half-thickness a, see Figure 4. The neutral fiber of the ring has a radius $R_s$. The ring is assumed to be thin and satisfies the condition

$$ 2a \ll R_s. \qquad (7) $$

The Euler-Bernoulli beam theory is considered. Any point on the ring bears six degrees-of-freedom ($u, v, w, \theta_s, \theta_r, \theta_w$), see Figure 4. From the small deformation assumption combined with condition (7), it

<footnote n="1">¹for more information, see also: https://nonsmooth.gricad-pages.univ-grenoble-alpes.fr/siconos/users guide/dynamical systems.html</footnote>
where $E, S = 2ab$ and $I_z = b(2a)^3/12$ represent Young's modulus, the cross-section area and the quadratic moment of inertia, respectively. The kinetic energy of the curved beam is

$$2E_c = \rho \int_{s_1}^{s_2} S(\dot{u}^2 + \dot{v}^2) + I_z \left(\frac{\dot{u}}{R_s} - \dot{v}, s\right)^2 ds. \quad (9)$$

Finally, the second term in Equation (9), stemming from the rotation of a straight cross-section, is neglected as in [19, 31, 40].

### 2.3.2. Space discretization

The flexible ring is discretized using curved beam finite elements [42]. It is divided into $N_e$ elements (equal to the number of nodes) of length $l_e = 2\pi R_s / N_e$. In order to guarantee continuity of the second derivatives in Equation (8), cubic Hermite polynomial shape functions are employed: $N_1(\chi) = 1 - 3\chi^2 + 2\chi^3$, $N_2(\chi) = \ell_e(\chi - 2\chi^2 + \chi^3)$, $N_3(\chi) = 3\chi^2 - 2\chi^3$ and $N_4(\chi) = \ell_e(-\chi^2 + \chi^3)$ where the dimensionless variable $\chi = s/l_e$ stands for the curvilinear coordinate to element length ratio. From the eight generalized coordinates $\mathbf{q}_s^e(t) = (u_1(t), u_{1,s}(t), v_1(t), v_{1,s}(t), u_2(t), u_{2,s}(t), v_2(t), v_{2,s}(t))^\top$ for each element$^2$, the displacement vector $\bar{\mathbf{q}}(\chi, t) = (u(\chi, t), v(\chi, t))^\top$ reads

$$\bar{\mathbf{q}}(\chi, t) \approx \mathbf{N}(\chi)\mathbf{q}_s^e(t) \quad (10)$$

with

$$\mathbf{N}(\chi) = \begin{bmatrix} N_1(\chi) & N_2(\chi) & 0 & 0 & N_3(\chi) & N_4(\chi) & 0 & 0 \\ 0 & 0 & N_1(\chi) & N_2(\chi) & 0 & 0 & N_3(\chi) & N_4(\chi) \end{bmatrix}. \quad (11)$$

The elemental stiffness matrix is computed from the discretized form of the potential energy in Equation (8)

$$2E_d \approx \mathbf{q}_s^{e\top} \left[ \int_0^1 (\mathbf{D}\mathbf{N}(\chi))^\top \mathbf{H}(\mathbf{D}\mathbf{N}(\chi)) \mathrm{d}\chi \right] \mathbf{q}_s^e = \mathbf{q}_s^{e\top} \mathbf{K}_e \mathbf{q}_s^e \quad (12)$$

where $\mathbf{D}$ stores differentiation operators and $\mathbf{H}$ is the elasticity matrix:

$$\mathbf{D} = \begin{bmatrix} \frac{1}{l_e} \frac{d}{d\chi} & \frac{1}{R_s} \\ 0 & \frac{1}{l_e^2} \frac{d^2}{d\chi^2} + \frac{1}{R_s^2} \end{bmatrix}, \quad \mathbf{H} = \begin{bmatrix} ES & 0 \\ 0 & EI \end{bmatrix}. \quad (13)$$

The same procedure is applied with the kinetic energy, that is

$$2E_c \approx \dot{\mathbf{q}}_s^{e\top} \left[ \int_0^1 \rho S \mathbf{N}(\chi)^{\top} \mathbf{N}(\chi) \mathrm{d}\chi \right] \dot{\mathbf{q}}_s^e = \dot{\mathbf{q}}_s^{e\top} \mathbf{M}_e \dot{\mathbf{q}}_s^e. \quad (14)$$

The elemental mass and stiffness matrices $\mathbf{M}_e$ and $\mathbf{K}_e$ are expanded in Appendix A. The global assembly procedure is verified by comparing the free ring eigenvalues for two and three nodal diameters. For a nodal diameter $n$, the analytical expression is [8]

$$f_n = \frac{n(n^2 - 1)}{2\pi R_s^2} \sqrt{\frac{EI_z}{2a(n^2 + 1)\rho b}} \quad (15)$$

which yields, for $n = 2$ and $n = 3$, $f_2 = 797$ Hz and $f_3 = 2255$ Hz, respectively. A convergence analysis of these two frequencies shows that $N_e \ge 24$ is required. The number of elements is always a multiple of the screws number so that there are always nodes coinciding with them. Note that these eigenfrequencies are already high reflecting a stiff structure because of the small annular radius $R_s$ of several centimetres$^3$. In the remainder, the reference case is $N_e = 36$. All the mechanical degrees of freedom for the stator (ring and point mass) are collected in the vector $\mathbf{q}_s$

### 2.3.3. Boundary conditions

As already said, component $A$ is a point mass in the simplified model and thus adds two degrees-of-freedom ($x_A, y_A$) to the system, which means that $\mathbf{q}_s$ is of size $4N_e + 2 = 146$ for reference case. In Figure 5, the

<footnote>$^2$For instance, $u_{1,s}$ corresponds to curvilinear derivative of displacement $u$ at node 1.</footnote>

<footnote>$^3$The accurate data is not provided because of the confidentiality policy.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4: Stator geometry, reference frame and degrees-of-freedom for a curved Bernoulli beam</figcaption></figure>

is possible to simplify the expressions of the deformation in 3D by involving only the three translations $(u, v, w)$, the torsion angle $(\theta_s)$ and their derivatives with respect to the curvilinear abscissa $s$, denoted with subscripts $s$ and $ss$ for compactness [19]:

$$ 
\begin{aligned}
\epsilon_{ss} &= \left(\frac{v}{R_s} + u_s\right) - \left(\frac{v}{R_s^2} + v_{,ss}\right)r + \left(\frac{\theta_s}{R_s} - w_{,ss}\right)z \\
\gamma_{rs} &= \frac{z}{R_s}(R_s \theta_{s,s} + w_{,s}) \\
\gamma_{zs} &= -\frac{r}{R_s}(R_s \theta_{s,s} + w_{,s}).
\end{aligned} 
\quad (8) $$

However, in the present paper, contact is only investigated in the $(X, Y)$ plane (see Figure 4). As a result, the torsion $\theta_s$ and axial displacement $w$ are neglected. The shear deformations $\gamma_{rs}$ and $\gamma_{zs}$ thus vanish and do not contribute to the deformation energy so that Hooke's law $\sigma_{ss} = E\epsilon_{ss}$ reads [20, 36]

$$ 2E_d = \int_{s_1}^{s_2} E S \left(\frac{v}{R_s} + u_s\right)^2 + EI_z \left(\frac{v}{R_s^2} + v_{,ss}\right)^2 ds \quad (9) $$

where $E, S = 2ab$ and $I_z = b(2a)^3/12$ represent Young's modulus, the cross-section area and the quadratic moment of inertia, respectively. The kinetic energy of the curved beam is

$$ 2E_c = \rho \int_{s_1}^{s_2} S(\dot{u}^2 + \dot{v}^2) + I_z \left(\frac{\dot{u}}{R_s} - \dot{v}, s\right)^2 ds. \quad (10) $$

Finally, the second term in Equation (10), stemming from the rotation of a straight cross-section, is neglected as in [20, 35, 44].

### 2.3.2. Space discretization

The flexible ring is discretized using curved beam finite elements [46]. It is divided into $N_e$ elements (equal to the number of nodes) of length $l_e = 2\pi R_s/N_e$. In order to guarantee continuity of the second derivatives in Equation (9), cubic Hermite polynomial shape functions are employed: $N_1(\chi) = 1 - 3\chi^2 + 2\chi^3$, $N_2(\chi) = l_e(\chi - 2\chi^2 + \chi^3)$, $N_3(\chi) = 3\chi^2 - 2\chi^3$ and $N_4(\chi) = l_e(-\chi^2 + \chi^3)$ where the dimensionless variable $\chi = s/l_e$ stands for the curvilinear coordinate to element length ratio. From the eight generalized coordinates $\mathbf{q}_s^e(t) = (u_1(t), u_{1,s}(t), v_1(t), v_{1,s}(t), u_2(t), u_{2,s}(t), v_2(t), v_{2,s}(t))^\top$ for each element$^2$, the displacement vector $\bar{\mathbf{q}}(\chi, t) = (u(\chi, t), v(\chi, t))^\top$ reads

$$ \bar{\mathbf{q}}(\chi, t) \approx \mathbf{N}(\chi)\mathbf{q}_s^e(t) \quad (11) $$

with

$$ \mathbf{N}(\chi) = \begin{bmatrix} N_1(\chi) & N_2(\chi) & 0 & 0 & N_3(\chi) & N_4(\chi) & 0 & 0 \\ 0 & 0 & N_1(\chi) & N_2(\chi) & 0 & 0 & N_3(\chi) & N_4(\chi) \end{bmatrix}. \quad (12) $$

The elemental stiffness matrix is computed from the discretized form of the potential energy in Equation (9)

$$ 2E_d \approx \mathbf{q}_s^{e\top} \left[ \int_0^1 (\mathbf{D}\mathbf{N}(\chi))^\top \mathbf{H}(\mathbf{D}\mathbf{N}(\chi)) \mathrm{d}\chi \right] \mathbf{q}_s^e = \mathbf{q}_s^{e\top} \mathbf{K}_e \mathbf{q}_s^e \quad (13) $$

<footnote>$^2$For instance, $u_{1,s}$ corresponds to curvilinear derivative of displacement $u$ at node 1.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5:** Reduction of components A and B into springs and mass connected on curved beams assembly (●) at screw nodes (●). Node positioning under deformed configuration with local displacements ($u_j, v_j$)</figcaption></figure>

screw nodes are emphasized as cyan nodes while the other nodes are black. The first node $j = 0$ lies on the X-axis and nodal identifiers $j = 0, \dots, N_e - 1$ are assigned counter-clockwise.

Figure 5 also shows that the position of a node $j$ in the deformed configuration, in the fixed frame (O, X, Y), is evaluated by the sum of the position at rest ($x_{0j}, y_{0j}$) = $R_s(\cos \alpha_j, \sin \alpha_j)$ (where $\alpha_j = 2\pi j / N_e$) and the projection of nodal displacement vector

$$ \begin{pmatrix} x_{sj}(\mathbf{q}_s) \\ y_{sj}(\mathbf{q}_s) \end{pmatrix} = \begin{pmatrix} x_{0j} \\ y_{0j} \end{pmatrix} + \begin{pmatrix} -\sin \alpha_j & \cos \alpha_j \\ \cos \alpha_j & \sin \alpha_j \end{pmatrix} \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} u_j \\ u_{j,s} \\ v_j \\ v_{j,s} \end{pmatrix}. \quad (16) $$

Note that $\mathbf{q}_{sj} \equiv (u_j, u_{j,s}, v_j, v_{j,s})^\top \in \mathbf{q}_s$ and that each of the stator coordinates reads in the condensed form

$$ \begin{cases} x_{sj}(\mathbf{q}_s) = x_{0j} + \mathbf{P}_x \mathbf{U} \mathbf{q}_{sj} \\ y_{sj}(\mathbf{q}_s) = y_{0j} + \mathbf{P}_y \mathbf{U} \mathbf{q}_{sj} \end{cases} \quad \text{where} \quad \mathbf{U} = \begin{pmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix} \quad (17) $$

along with $\mathbf{P}_x = (-\sin \alpha_j \cos \alpha_j)$ and $\mathbf{P}_y = (\cos \alpha_j \sin \alpha_j)$.

A screw node is characterized by $j = \frac{nN_e}{N_{sc}}$ where $n = 0, N_{sc} - 1$. In order to avoid confusion with classical nodes, screw nodes are identified by the additional variable $n$. The potential energy at a screw node $n$ has two contributions provided by the connecting components:

* Component B: $E_{p,B} = \frac{1}{2}(k_{Bu}u_n^2 + k_{Bv}v_n^2)$ and the stiffness matrix $\mathbf{K}_B$ is thus diagonal.

* Component A: $E_{p,A} = \frac{1}{2}k_A(v_n - x_A \cos \alpha_n - y_A \sin \alpha_n)^2$, where $\alpha_n = 2\pi n / N_{sc}$ is the angular position of screw $n$ on the ring. Therefore, the stiffness matrix $\mathbf{K}_A$ is non-zero only on the elements involving the displacements ($v_n, x_A, y_A$). In other words, the condensed stiffness matrix reflecting component A at screw node $n$ is

$$ \mathbf{K}_{A,n} = k_A \begin{bmatrix} 1 & -\cos \alpha_n & -\sin \alpha_n \\ -\cos \alpha_n & \cos^2 \alpha_n & \cos \alpha_n \sin \alpha_n \\ -\sin \alpha_n & \cos \alpha_n \sin \alpha_n & \sin^2 \alpha_n \end{bmatrix}. \quad (18) $$

Note that coupling terms $\cos \alpha_n \sin \alpha_n$ exist between the displacements $x_A$ and $y_A$. During the matrix assembly, matrices from Equation (18) are extended and concatenated<fnref n="4" />

$$ \mathbf{K}_A = \bigcup_{n=0}^{N_{sc}-1} \mathbf{K}_{A,n}. \quad (19) $$

<footnote n="4">⁴Concatenation is defined by the union symbol in the present paper.</footnote>
where **D** stores differentiation operators and **H** is the elasticity matrix:

$$ \mathbf{D} = \begin{bmatrix} \frac{1}{\ell_e} \frac{d}{d\chi} & \frac{1}{R_s} \\ 0 & \frac{1}{\ell_e^2} \frac{d^2}{d\chi^2} + \frac{1}{R_s^2} \end{bmatrix}, \quad \mathbf{H} = \begin{bmatrix} ES & 0 \\ 0 & EI \end{bmatrix}. \qquad (14) $$

The same procedure is applied with the kinetic energy, that is

$$ 2E_c \approx \dot{\mathbf{q}}_s^{e\top} \left[ \int_0^1 \rho S \mathbf{N}(\chi)^{\top} \mathbf{N}(\chi) d\chi \right] \dot{\mathbf{q}}_s^e = \dot{\mathbf{q}}_s^{e\top} \mathbf{M}_e \dot{\mathbf{q}}_s^e. \qquad (15) $$

The elemental mass and stiffness matrices $\mathbf{M}_e$ and $\mathbf{K}_e$ are expanded in Appendix A. The global assembly procedure is verified by comparing the free ring eigenvalues for two and three nodal diameters. For a nodal diameter $n$, the analytical expression is [9]

$$ f_n = \frac{n(n^2 - 1)}{2\pi R_s^2} \sqrt{\frac{EI_z}{2a(n^2 + 1)\rho b}} \qquad (16) $$

which yields, for $n = 2$ and $n = 3$, $f_2 = 797$ Hz and $f_3 = 2255$ Hz, respectively. A convergence analysis of these two frequencies shows that $N_e \ge 24$ is required. The number of elements is always a multiple of the screws number so that there are always nodes coinciding with them. Note that these eigenfrequencies are already high reflecting a stiff structure because of the small annular radius $R_s$ of several centimeters<fnref n="3" />. In the remainder, the reference case is $N_e = 36$. All the mechanical degrees of freedom for the stator (ring and point mass) are collected in the vector $\mathbf{q}_s$

### 2.3.3. Boundary conditions

As already said, component A is a point mass in the simplified model and thus adds two degrees-of-freedom $(x_A, y_A)$ to the system, which means that $\mathbf{q}_s$ is of size $4N_e + 2 = 146$ for reference case. In Figure 5, the

<figure><img src="image_6.png" /><figcaption>**Figure 5:** Reduction of components A and B into springs and mass connected on curved beams assembly (●) at screw nodes (●). Node positioning under deformed configuration with local displacements $(u_j, v_j)$</figcaption></figure>

screw nodes are emphasized as cyan nodes while the other nodes are black. The first node $j = 0$ lies on the X-axis and nodal identifiers $j = 0, \dots, N_e - 1$ are assigned counter-clockwise.

The Figure 5 also shows that the position of a node $j$ in deformed configuration, in the fixed frame $(O, X, Y)$, is evaluated by the sum of the position at rest $(x_{0j}, y_{0j}) = R_s(\cos \alpha_j, \sin \alpha_j)$ (where $\alpha_j = 2\pi j / N_e$) and the projection of nodal displacement vector $\mathbf{q}_{sj} = (u_j, u_{j,s}, v_j, v_{j,s})^\top$

$$ \begin{pmatrix} x_{sj}(\mathbf{q}_s) \\ y_{sj}(\mathbf{q}_s) \end{pmatrix} = \begin{pmatrix} x_{0j} \\ y_{0j} \end{pmatrix} + \begin{pmatrix} -\sin \alpha_j & \cos \alpha_j \\ \cos \alpha_j & \sin \alpha_j \end{pmatrix} \begin{pmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix} (u_j, u_{j,s}, v_j, v_{j,s})^\top. \qquad (17) $$

Note that $\mathbf{q}_{sj} \in \mathbf{q}_s$ and that each of the stator coordinates writes under the condensed form

$$ \left\{ \begin{array}{l} x_{sj}(\mathbf{q}_s) = x_{0j} + \mathbf{P}_x \mathbf{U} \mathbf{q}_{sj} \\ y_{sj}(\mathbf{q}_s) = y_{0j} + \mathbf{P}_y \mathbf{U} \mathbf{q}_{sj} \end{array} \right. . \qquad (18) $$

<footnote n="3">³The accurate data is not provided because of the confidentiality policy.</footnote>
Yet, for $N_{sc}$ evenly spaced screws on the ring, these elements linking $x_A$ and $y_A$ are summed and the result is zero since

$$ \sum_{n=0}^{N_{sc}-1} \cos \frac{2\pi n}{N_{sc}} \sin \frac{2\pi n}{N_{sc}} = 0 \qquad (20) $$

These extra-diagonal elements in Equation (18) vanish in the assembled matrix $K_A$.

The global mass matrix $M_s$ and stiffness matrix $K_s$ are built by assembling elemental mass and stiffness matrices for the ring (resp. $M_e$ and $K_e$) along with the stiffness matrices of the $N_{sc}$ connecting springs $K_B$, $K_A$. With a 36-node mesh on the ring and a point mass located by a translation pair $x_A$, $y_A$, the whole stator model includes a total of 146 degrees-of-freedom. The damping matrix is set proportional to the stiffness matrix, that is $D_s = r_k K_s$ where $r_k = 10^{-6}$ is the Rayleigh coefficient. The computation of eigenvectors shows that the damping ratio is about 4% for the first then modes which is in agreement with the observed measurements.

## 2.4. Thermomechanical coupling

In contrast to the turbine, the numerical investigation of the thermal transient is performed for the stator. Piezothermal effects are neglected and the heat source only stems from the flux at the contact surface $S_c$. The finite element method applied to the heat equation for the stator of volume $V_s$ leads to the integral form

$$ \left[ \int_{V_s} \rho c N_{\theta}^{\top} N_{\theta} \, dV \right] \dot{T}_e + \left[ \int_{V_s} \lambda \nabla N_{\theta}^{\top} \cdot \nabla N_{\theta} \, dV \right] T_e = \int_{S_c} N_{\theta}^{\top} \dot{Q}_s \, dS \qquad (21) $$

where $c$ is the heat capacity, $\lambda$, the heat conductivity, $\dot{Q}_s$, the thermal flux density on $S_c$ (in W/m$^2$) and $N_{\theta}$, the shape functions for the temperature field. It is now required to properly quantify $\dot{Q}_s$. The present work aims to only provide a preliminary insight of the temperature transient, in line with [4, 47] where the authors considered the heat transfer formulation proposed in [25]. They accounted for conduction, Coulomb's friction model and Archard's law for wear. The essential mathematical developments are described in the present paragraph. Let us consider two bodies, named rotor and stator, of respective volumes $V_r$ and $V_s$. These bodies are in contact on a surface $S_c$. The general energy balance states that

$$ \left( \mu p_N(\mathbf{x}, t) + \frac{k_a}{3p_s} p_N^2(\mathbf{x}, t) \right) \|\mathbf{v}_T(\mathbf{x}, t)\| - \dot{Q}_r(\mathbf{x}, t) - \dot{Q}_s(\mathbf{x}, t) = 0 \quad \forall \mathbf{x} \text{ on } S_c \qquad (22) $$

where $k_a$ stands for the Archard coefficient, $p_s$ is the material hardness and $\dot{Q}_r$ and $\dot{Q}_s$ are the heat flux (unit W/m$^2$) entering in the rotor and stator through the contact interface $S_c$. $p_N$ is the normal contact pressure applied on $S_c$ and $\mathbf{v}_T$ stands for the relative tangential velocity (or sliding velocity) between the bodies. The present paper does not account for wear mechanisms and it is possible to reduce the expression of frictional powers for rotor and stator, respectively,

$$ \left\{ \begin{array}{l} \dot{Q}_r(\mathbf{x}, t) = \frac{\vartheta_r}{\vartheta_r + \vartheta_s} (\vartheta_s(T_s(\mathbf{x}, t) - T_r(\mathbf{x}, t)) + \mu \|\mathbf{v}_T(\mathbf{x}, t)\|) p_N(\mathbf{x}, t) \\ \\ \dot{Q}_s(\mathbf{x}, t) = \frac{\vartheta_s}{\vartheta_r + \vartheta_s} (\vartheta_r(T_r(\mathbf{x}, t) - T_s(\mathbf{x}, t)) + \mu \|\mathbf{v}_T(\mathbf{x}, t)\|) p_N(\mathbf{x}, t) \end{array} \right. \quad \forall \mathbf{x} \text{ on } S_c \qquad (23) $$

where $\vartheta$ is the thermal conductance of contact, $T_r$ and $T_s$ represent the rotor and stator temperatures. In reality, the conductance depends on material properties and the surface deformation but is commonly linearized with respect to $p_N$ [22], leading to Equation (23). The first term in this equation is the contribution of conduction to the heat source and the second term represents frictional power. Since temperature distribution in the rotor is not investigated, the first term is ignored and Equation (23) simplifies to

$$ \dot{Q}_s = \frac{\vartheta_s}{\vartheta_r + \vartheta_s} \mu \|\mathbf{v}_T\| p_N. \qquad (24) $$

Eventually, it is assumed that the rotor and the stator have identical contact conductances [4, 25] which implies that frictional power is equally split between the rotor and stator and that $\frac{\vartheta_s}{\vartheta_r + \vartheta_s} = 0.5$. One thus retrieves the heat partitioning coefficient used in [17].

In the present work, the stator is discretized by beam elements with constant section $S$. This means that the heat equation can only be solved on the neutral fiber along the curvilinear coordinate $s$, see Figure 4.
where $\mathbf{P}_x = (-\sin \alpha_j \cos \alpha_j)$, $\mathbf{P}_y = (\cos \alpha_j \sin \alpha_j)$ and $\mathbf{U} = \begin{pmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix}$.

A screw node is characterized by $j = \frac{nN_e}{N_{sc}}$ where $n = 0, N_{sc} - 1$. In order to avoid confusion with classical nodes, screw nodes are identified by the additional variable $n$. The potential energy at a screw node $n$ has two contributions provided by the connecting components:

* Component B: $E_{p,B} = \frac{1}{2}(k_{Bu}u_n^2 + k_{Bv}v_n^2)$ and the stiffness matrix $\mathbf{K}_B$ is thus diagonal.

* Component A: $E_{p,A} = \frac{1}{2}k_A(v_n - x_A \cos \alpha_n - y_A \sin \alpha_n)^2$, where $\alpha_n = 2\pi n/N_{sc}$ is the angular position of screw $n$ on the ring. Therefore, the stiffness matrix $\mathbf{K}_A$ is non-zero only on the elements involving the displacements $(v_n, x_A, y_A)$. In other words, the condensed stiffness matrix reflecting component A at screw node $n$ is

$$ \mathbf{K}_{A,n} = k_A \begin{bmatrix} 1 & -\cos \alpha_n & -\sin \alpha_n \\ -\cos \alpha_n & \cos^2 \alpha_n & \cos \alpha_n \sin \alpha_n \\ -\sin \alpha_n & \cos \alpha_n \sin \alpha_n & \sin^2 \alpha_n \end{bmatrix}. \qquad (19) $$

Note that coupling terms $\cos \alpha_n \sin \alpha_n$ exist between the displacements $x_A$ and $y_A$. During the matrix assembly, matrices from Equation (19) are extended and concatenated<fnref n="4" />

$$ \mathbf{K}_A = \bigcup_{n=0}^{N_{sc}-1} \mathbf{K}_{A,n}. \qquad (20) $$

Yet, for $N_{sc}$ evenly spaced screws on the ring, these elements linking $x_A$ and $y_A$ are summed and the result is zero since

$$ \sum_{n=0}^{N_{sc}-1} \cos \frac{2\pi n}{N_{sc}} \sin \frac{2\pi n}{N_{sc}} = 0 \qquad (21) $$

These extra-diagonal elements in Equation (19) vanish in the assembled matrix $\mathbf{K}_A$.

The global mass matrix $\mathbf{M}_s$ and stiffness matrix $\mathbf{K}_s$ are built by assembling elemental mass and stiffness matrices for the ring (resp. $\mathbf{M}_e$ and $\mathbf{K}_e$) along with the stiffness matrices of the $N_{sc}$ connecting springs $\mathbf{K}_B$, $\mathbf{K}_A$. With a 36-node mesh on the ring and a point mass located by a translation pair $x_A, y_A$, the whole stator model includes a total of 146 degrees-of-freedom. The damping matrix is set proportional to the stiffness matrix, that is $\mathbf{D}_s = r_k \mathbf{K}_s$ where $r_k = 10^{-6}$ is the Rayleigh coefficient. The computation of eigenvectors shows that the damping ratio is about 4% for the first then modes which is in agreement with the observed measurements.

## 2.4. Thermomechanical coupling

In contrast to the turbine, the numerical investigation of the thermal transient is performed for the stator. Piezothermal effects are neglected and the heat source only stems from the flux at the contact surface $S_c$. The finite element method applied to the heat equation for the stator of volume $V_s$ leads to the integral form

$$ \left[ \int_{V_s} \rho c \mathbf{N}_{\theta}^{\top} \mathbf{N}_{\theta} \, dV \right] \dot{\mathbf{T}}_e + \left[ \int_{V_s} \lambda \nabla \mathbf{N}_{\theta}^{\top} \cdot \nabla \mathbf{N}_{\theta} \, dV \right] \mathbf{T}_e = \int_{S_c} \mathbf{N}_{\theta}^{\top} \dot{\mathbf{Q}}_s \, dS \qquad (22) $$

where $c$ is the heat capacity, $\lambda$, the heat conductivity, $\dot{\mathbf{Q}}_s$, the thermal flux density on $S_c$ (in W/m²) and $\mathbf{N}_{\theta}$, the shape functions for the temperature field. It is now required to properly quantify $\dot{\mathbf{Q}}_s$. The present work aims to only provide a preliminary insight of the temperature transient, in line with [4, 52] where the authors considered the heat transfer formulation proposed in [29]. They accounted for conduction, Coulomb's friction model and Archard's law for wear. The essential mathematical developments are described in the present paragraph. Let us consider two bodies, named rotor and stator, of respective volumes $V_r$ and $V_s$. These bodies are in contact on a surface $S_c$. The general energy balance states that

$$ \left( \mu p_N(\mathbf{x}, t) + \frac{k_a}{3p_s} p_N^2(\mathbf{x}, t) \right) \|\mathbf{v}_T(\mathbf{x}, t)\| - \dot{\mathbf{Q}}_r(\mathbf{x}, t) - \dot{\mathbf{Q}}_s(\mathbf{x}, t) = 0 \quad \forall \mathbf{x} \text{ on } S_c \qquad (23) $$

where $k_a$ stands for the Archard coefficient, $p_s$ is the material hardness and $\dot{\mathbf{Q}}_r$ and $\dot{\mathbf{Q}}_s$ are the heat flux (unit W/m²) entering in the rotor and stator through the contact interface $S_c$. $p_N$ is the normal contact

<footnote n="4">⁴Concatenation is defined by the union symbol in the present paper.</footnote>
As a result, temperature is discretized in space by a basic linear interpolation, involving the shape function $N_{\theta}(\chi) = (1 - \chi, \chi)$, $\chi \in [0, 1]$. The discretized integral form of heat equation then reads

$$ \left[ \int_{0}^{1} \rho c S \ell_{e} N_{\theta}^{\top} N_{\theta} d\chi \right] \dot{T}_{e} + \left[ \int_{0}^{1} \frac{\lambda S}{\ell_{e}} N_{\theta, \chi}^{\top} N_{\theta, \chi} d\chi \right] T_{e} = \int_{S_{c}} N_{\theta}^{\top} \dot{Q}_{s} dS. \quad (25) $$

Simplifying Equation (25) leads to the ODE

$$ C_e^{\theta\theta} \dot{T} + K_e^{\theta\theta} T = \dot{Q}_{s,e} \quad \text{with} \quad C_e^{\theta\theta} = \frac{\rho c S \ell_e}{6} \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix} \quad \text{and} \quad K_e^{\theta\theta} = \frac{\lambda S}{\ell_e} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}. \quad (26) $$

The expression for assembled nodal heat vector $\dot{Q}_s$ depends on the contact treatment methodology discussed in the next section.

As already mentioned, friction is not the only phenomenon coupling mechanics and thermodynamics. Thermal expansion is also of interest and the thermal dilation matrix must be computed. Its integral form for a single element is [46]

$$ K_e^{q\theta} = - \int_{V_s} \alpha (DN)^{\top} H I_d N_{\theta} dV \quad (27) $$

where $I_d$ is the vector reflecting that temperature only affects normal deformation.<fnref n="5" /> In the case of a Bernoulli curved beam, only the temperature contributes to the elongation $\epsilon_{ss}$ of the neutral fiber [43]. Consequently $I_d = (1 \ 0)^{\top}$. With $D$, $N$, and $H$ already known, the elementary expansion matrix is easily computed and its expression is found in Appendix A. Thermal expansion adds a term in the stiffness matrix of the form $K^{q\theta}(T_s - T_{ref})$ with $T_{ref}$, the reference temperature.

Eventually, the governing equations for the turbine dynamics and the stator thermomechanics including contact conditions is

$$ \begin{bmatrix} \mathbf{M}_r(\mathbf{q}_r) & 0 & 0 \\ 0 & \mathbf{M}_s & 0 \\ 0 & 0 & 0 \end{bmatrix} \begin{pmatrix} \ddot{\mathbf{q}}_r \\ \ddot{\mathbf{q}}_s \\ \dot{\mathbf{T}}_s \end{pmatrix} + \begin{bmatrix} \dot{\theta}\mathbf{G}_r + \mathbf{D}_b & 0 & 0 \\ 0 & \mathbf{D}_s & 0 \\ 0 & 0 & \mathbf{C}^{\theta\theta} \end{bmatrix} \begin{pmatrix} \dot{\mathbf{q}}_r \\ \dot{\mathbf{q}}_s \\ \dot{\mathbf{T}}_s \end{pmatrix} + \begin{bmatrix} \mathbf{K}_b & 0 & 0 \\ 0 & \mathbf{K}_s & \mathbf{K}^{q\theta} \\ 0 & 0 & \mathbf{K}^{\theta\theta} \end{bmatrix} \begin{pmatrix} \mathbf{q}_r \\ \mathbf{q}_s \\ \mathbf{T}_s \end{pmatrix} = \begin{pmatrix} \mathbf{f}_r^{\text{ext}} \\ \mathbf{K}^{q\theta} \mathbf{T}_{\text{ref}} \\ 0 \end{pmatrix} + \begin{pmatrix} \mathbf{f}_r^c \\ \mathbf{f}_s^c \\ \dot{\mathbf{Q}}_s \end{pmatrix}. \quad (28) $$

In the remainder, contact efforts and torque applied to the turbine and the stator are collected in the vector $\mathbf{f}_c^q = (\mathbf{f}_r^{c\top} \ \mathbf{f}_s^{c\top})^\top$ and the generalized coordinates vector is

$$ \mathbf{q} = (\mathbf{q}_r, \mathbf{q}_s, \mathbf{T}_s)^{\top} \quad (29) $$

with turbine related unknowns ($\mathbf{q}_r$) and stator-related unknowns ($\mathbf{q}_s, \mathbf{T}_s$). For the reference case (36 nodes), the matrices are of size $187 \times 187$.

The initial conditions for the overall system state that the stator is at rest and no thermal pre-stress is applied, implying that $\mathbf{q}_s(0) = \dot{\mathbf{q}}_s(0) = \mathbf{0}$ and $\mathbf{T}_s(0) = \mathbf{T}_{\text{ref}}$. Since the heat equation is linear, imposing a homogeneous initial condition for temperature and suppressing term $\mathbf{K}^{q\theta}\mathbf{T}_{\text{ref}}$ in the second member of Equation (28) provides the same solution just offset by $\mathbf{T}_{\text{ref}}$. Consequently, the equivalent initial condition writes $\mathbf{T}_s(0) = \mathbf{T}_{\text{ref}} = \mathbf{0}$ for the sake of simplicity.

# 3. Solution method

## 3.1. Unilateral contact and friction treatment

Unilateral contact conditions prevent mechanical bodies to penetrate each other and always create a reaction force of magnitude $F_N$ in direction $\mathbf{n}$ normal to the contact surface. Defining $d_N$ as the distance between the bodies, a classical formulation of unilateral contact constraints is given by the Signorini conditions expressed at the displacement level [52]: the non-penetration constraint $d_N \ge 0$, the non-traction constraint $F_N \ge 0$ and $d_N F_N = 0$, enforcing the fact that both constraints cannot be active simultaneously.

A large majority of works and industrial tools simplify the Signorini conditions by regularization techniques [11, 34, 20]. However, the present work aims to strictly enforce the non-smooth law and to

<footnote n="5">⁵With tensor formulation, the thermoelastic law for an isotropic material is $\sigma^{\otimes 2} = E^{\otimes 4} : [\epsilon^{\otimes 2} - \alpha(T - T_{\text{ref}})I_d^{\otimes 2}]$, where $I_d$ is the identity second order tensor. If the deformation and constraints tensors are condensed through the Voigt notation as column arrays $\epsilon$ and $\sigma \in \mathbb{R}^6$, then $I_d = (1, 1, 1, 0, 0, 0)$.</footnote>
pressure applied on $S_c$ and $v_T$ stands for the relative tangential velocity (or sliding velocity) between the bodies. The present paper does not account for wear mechanisms and it is possible to reduce the expression of frictional powers for rotor and stator, respectively,

$$ egin{cases} \dot{Q}_r(\mathbf{x}, t) = \frac{\vartheta_r}{\vartheta_r + \vartheta_s} (\vartheta_s(T_s(\mathbf{x}, t) - T_r(\mathbf{x}, t)) + \mu \|v_T(\mathbf{x}, t)\|) p_N(\mathbf{x}, t) \\ \dot{Q}_s(\mathbf{x}, t) = \frac{\vartheta_s}{\vartheta_r + \vartheta_s} (\vartheta_r(T_r(\mathbf{x}, t) - T_s(\mathbf{x}, t)) + \mu \|v_T(\mathbf{x}, t)\|) p_N(\mathbf{x}, t) \end{cases} \quad \forall \mathbf{x} \text{ on } S_c \quad (24) $$

where $\vartheta$ is the thermal conductance of contact, $T_r$ and $T_s$ represent the rotor and stator temperatures. In reality, the conductance depends on material properties and the surface deformation but is commonly linearized with respect to $p_N$ [24], leading to Equation (24). The first term in this equation is the contribution of conduction to the heat source and the second term represents frictional power. Since temperature distribution in the rotor is not investigated, the first term is ignored and Equation (24) simplifies to

$$ \dot{Q}_s = \frac{\vartheta_s}{\vartheta_r + \vartheta_s} \mu \|v_T\| p_N. \quad (25) $$

Eventually, it is assumed that the rotor and the stator have identical contact conductances [4, 29] which implies that frictional power is equally split between the rotor and stator and that $\frac{\vartheta_s}{\vartheta_r + \vartheta_s} = 0.5$. One thus retrieves the heat partitioning coefficient used in [18].

In the present work, the stator is discretized by beam elements with constant section $S$. This means that the heat equation can only be solved on the neutral fiber along the curvilinear coordinate $s$, see Figure 4. As a result, temperature is discretized in space by a basic linear interpolation, involving the shape function $N_{\vartheta}(\chi) = (1 - \chi, \chi)$, $\chi \in [0, 1]$. The discretized integral form of heat equation then reads

$$ \left[ \int_{0}^{1} \rho c S \ell_{e} N_{\vartheta}^{\top} N_{\vartheta} \, d\chi \right] \dot{T}_{e} + \left[ \int_{0}^{1} \frac{\lambda S}{\ell_{e}} N_{\vartheta, \chi}^{\top} N_{\vartheta, \chi} \, d\chi \right] T_{e} = \int_{S_{c}} N_{\vartheta}^{\top} \dot{Q}_{s} \, dS. \quad (26) $$

Simplifying Equation (26) leads to the ODE

$$ C_e^{\vartheta\vartheta} \dot{T} + K_e^{\vartheta\vartheta} T = \dot{Q}_{s,e} \quad \text{with} \quad C_e^{\vartheta\vartheta} = \frac{\rho c S \ell_e}{6} \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix} \quad \text{and} \quad K_e^{\vartheta\vartheta} = \frac{\lambda S}{\ell_e} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}. \quad (27) $$

The expression for assembled nodal heat vector $\dot{Q}_s$ depends on the contact treatment methodology discussed in the next section.

As already mentioned, friction is not the only phenomenon coupling mechanics and thermodynamics. Thermal expansion is also of interest and the thermal dilation matrix must be computed. Its integral form for a single element is [51]

$$ K_e^{q\vartheta} = - \int_{V_s} \alpha (DN)^{\top} H I_d N_{\vartheta} \, dV \quad (28) $$

where $I_d$ is the vector reflecting that temperature only affects normal deformation.<fnref n="5" /> In the case of a Bernoulli curved beam, only the temperature contributes to the elongation $\epsilon_{ss}$ of the neutral fiber [47]. Consequently $I_d = (1 \ 0)^{\top}$. With $D$, $N$, and $H$ already known, the elementary expansion matrix is easily computed and its expression is found in Appendix A. Thermal expansion adds a term in the stiffness matrix of the form $K^{q\vartheta}(T_s - T_{ref})$ with $T_{ref}$, the reference temperature.

Eventually, the governing equations for the turbine dynamics and the stator thermomechanics including contact conditions is

$$ \begin{bmatrix} \mathbf{M}_r(\mathbf{q}_r) & 0 & 0 \\ 0 & \mathbf{M}_s & 0 \\ 0 & 0 & 0 \end{bmatrix} \begin{pmatrix} \ddot{\mathbf{q}}_r \\ \ddot{\mathbf{q}}_s \\ \dot{T}_s \end{pmatrix} + \begin{bmatrix} \dot{\theta}\mathbf{G}_r + \mathbf{D}_b & 0 & 0 \\ 0 & \mathbf{D}_s & 0 \\ 0 & 0 & \mathbf{C}^{\vartheta\vartheta} \end{bmatrix} \begin{pmatrix} \dot{\mathbf{q}}_r \\ \dot{\mathbf{q}}_s \\ \dot{T}_s \end{pmatrix} + \begin{bmatrix} \mathbf{K}_b & 0 & 0 \\ 0 & \mathbf{K}_s & \mathbf{K}^{q\vartheta} \\ 0 & 0 & \mathbf{K}^{\vartheta\vartheta} \end{bmatrix} \begin{pmatrix} \mathbf{q}_r \\ \mathbf{q}_s \\ \mathbf{T}_s \end{pmatrix} = \begin{pmatrix} \mathbf{f}_r^{\text{ext}} \\ \mathbf{K}^{q\vartheta} \mathbf{T}_{\text{ref}} \\ 0 \end{pmatrix} + \begin{pmatrix} \mathbf{f}_r^c \\ \dot{\mathbf{Q}}_s \\ \dot{\mathbf{Q}}_s \end{pmatrix}. \quad (29) $$

<footnote n="5">⁵With tensor formulation, the thermoelastic law for an isotropic material is $\sigma^{\otimes 2} = E^{\otimes 4} : [\epsilon^{\otimes 2} - \alpha(T - T_{ref})I_d^{\otimes 2}]$, where $I_d$ is the identity second order tensor. If the deformation and constraints tensors are condensed through the Voigt notation as column arrays $\epsilon$ and $\sigma \in \mathbb{R}^6$, then $I_d = (1, 1, 1, 0, 0, 0)$.</footnote>
prevent non-physical and residual penetration. A result, preference is given to solution techniques involving constraints management and Lagrange multipliers. A widely known procedure allowing to solve dynamics with contact occurrence in time domain is the Carpenter algorithm [9]. Various published works [42, 4, 32] selected the Carpenter algorithm for its simplicity of implementation while giving results in agreement with experiments [48]. Note that other solutions methods with Lagrange multipliers are also available such as the Moreau-Jean [35, 28] or the Paoli-Schatzman [38] algorithms. They include the non-smooth Newton impact law [1] in the Signorini conditions.

Friction is accounted for by introducing a force $\mathbf{F}_f$ oriented in the direction opposed to the relative speed $\mathbf{v}_T$ at the contact interface between the interacting bodies. Its magnitude is commonly dictated by the Coulomb dry friction model stating that the friction force must satisfy

$$ \|\mathbf{F}_f\| = \mu F_N, \quad \text{if} \quad \mathbf{v}_T \neq \mathbf{0} \qquad (30a) $$

$$ \|\mathbf{F}_f\| \le \mu F_N, \quad \text{if} \quad \mathbf{v}_T = \mathbf{0} \qquad (30b) $$

Besides, the implementation of Coulomb's friction law requires careful attention because of the non-smoothness of the sticking friction phase governed by Equation (30b). If one desires to strictly enforce the Coulomb friction law as written in Equations (30a) and (30b), the exact approach for frictional contact treatment is based on mathematical concepts from convex analysis and defines contact efforts in the Coulomb cone [1]. Otherwise, it is possible to regularize the Coulomb friction law [11, 17] or to consider a sliding-only friction law characterized by

$$ F_f = \mathbf{F}_f \cdot \mathbf{t} = -\mu F_N \mathrm{sign}(\mathbf{v}_T) \quad \text{with} \quad \mathbf{v}_T = \mathbf{v}_T \cdot \mathbf{t}. \qquad (31) $$

where $\mathbf{t}$ is the unit tangential vector in local contact frame. The above relation implies that the relative velocity is known and its assessment is further developed in the present section. This formulation is commonly considered when the rotor operates at high rotational speeds [32, 4, 42].

In [26], Carpenter and Moreau-Jean algorithms were compared on rotor-stator rubbing between rigid bodies with various friction models. Good agreement is reported and pure rolling motion (i.e. sticking friction) can take place during the simulation. More interestingly, rotational velocity transient was mainly affected during sliding friction phases. This explains why sliding friction only, see Equation (31), is considered in the present model. Accordingly, it is decided to manage contact occurrences with the Carpenter algorithm. All unknowns are grouped into the vector $\mathbf{q} = (\mathbf{q}_r, \mathbf{q}_s, \mathbf{T}_s)$ and global mass, damping, gyroscopic and stiffness matrices are named $\mathbf{M}(\mathbf{q})$, $\mathbf{D}$, $\mathbf{G}$, and $\mathbf{K}$, respectively, such that Equation (28) becomes

$$ \mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} + (\mathbf{D} + \dot{\theta}\mathbf{G})\dot{\mathbf{q}} + \mathbf{K}\mathbf{q} = \mathbf{f}^{\text{ext}} + \mathbf{f}_c. \qquad (32) $$

Dependency of the mass matrix on the turbine dofs is shown to stress that it is not constant.

## 3.2. Modified Carpenter algorithm

In Carpenter algorithm, Equation (32) is discretized in time using an integrator from the Newmark family [24]. In order to avoid confusion with other symbols, numerical damping and stability parameters of the Newmark schemes are respectively denoted $\beta_1$ and $\beta_2$. Simulation time is divided into $N$ intervals of duration $h$. For studies involving mechanics only, the centred finite-difference scheme with $\beta_1 = 0.5$ and $\beta_2 = 0$ is appropriate but Thorin & al. [47] showed that $\beta_1$ must be slightly offset because of stability issues to solve the heat equation with thermomechanical coupling. Consequently, the integrator is chosen in this work with $\beta_1 = 0.501$ and $\beta_2 = 0$. In the present paper, to reduce the amount of notation, it is decided to reflect dependencies regarding unknown quantities only. In discrete time, the unknown vector is $\mathbf{q}_{i+1}$ at time step $i + 1$ and the contact-free equations of motion become

$$ \bar{\mathbf{M}}(\mathbf{q}_{i+1})\mathbf{q}_{i+1} + \bar{\mathbf{D}}(\mathbf{q}_{i+1})\mathbf{q}_i + \bar{\mathbf{K}}(\mathbf{q}_{i+1})\mathbf{q}_{i-1} = \bar{\mathbf{f}}_i, \quad i = 1, \dots, N \qquad (33) $$

where

$$ \begin{aligned} \bar{\mathbf{M}}(\mathbf{q}_{i+1}) &:= \frac{\mathbf{M}(\mathbf{q}_i)}{h^2} + \beta_1 \frac{\mathbf{D} + \dot{\theta}_i(\mathbf{q}_{i+1})\mathbf{G}}{h} \\ \bar{\mathbf{D}}(\mathbf{q}_{i+1}) &:= \frac{2\mathbf{M}(\mathbf{q}_i)}{h^2} + \frac{2\beta_1 - 1}{h} (\mathbf{D} + \dot{\theta}_i(\mathbf{q}_{i+1})\mathbf{G}) - (0.5 + \beta_1)\mathbf{K} \\ \bar{\mathbf{K}}(\mathbf{q}_{i+1}) &:= \frac{\mathbf{M}(\mathbf{q}_i)}{h^2} + \frac{\beta_1 - 1}{h} (\mathbf{D} + \dot{\theta}_i(\mathbf{q}_{i+1})\mathbf{G}) + (0.5 - \beta_1)\mathbf{K} \end{aligned} \qquad (34) $$

$$ \bar{\mathbf{f}}_i := (1 - \beta_1)\mathbf{f}_{i-1}^{\text{ext}} + (1 + \beta_1)\mathbf{f}_i^{\text{ext}} $$
In the remainder, contact efforts and torque applied to the turbine and the stator are collected in the vector $\mathbf{f}_c^q = (\mathbf{f}_r^{c\top} \mathbf{f}_s^{c\top})^\top$ and the generalized coordinates vector is

$$ \mathbf{q} = (\mathbf{q}_r, \mathbf{q}_s, \mathbf{T}_s)^\top \qquad (30) $$

with turbine related unknowns ($\mathbf{q}_r$) and stator-related unknowns ($\mathbf{q}_s, \mathbf{T}_s$). For the reference case (36 nodes), the matrices are of size 187 × 187.

The initial conditions for the overall system state that the stator is at rest and no thermal pre-stress is applied, implying that $\mathbf{q}_s(0) = \dot{\mathbf{q}}_s(0) = \mathbf{0}$ and $\mathbf{T}_s(0) = \mathbf{T}_{\text{ref}}$. Since the heat equation is linear, imposing a homogeneous initial condition for temperature and suppressing term $\mathbf{K}^{q\theta}\mathbf{T}_{\text{ref}}$ in the second member of Equation (29) provides the same solution just offset by $\mathbf{T}_{\text{ref}}$. Consequently, the equivalent initial condition writes $\mathbf{T}_s(0) = \mathbf{T}_{\text{ref}} = \mathbf{0}$ for the sake of simplicity.

### 3. Solution method

#### 3.1. Unilateral contact and friction treatment

Unilateral contact conditions prevent mechanical bodies to penetrate each other and always create a reaction force of magnitude $F_N$ in direction $\mathbf{n}$ normal to the contact surface. Defining $d_N$ as the distance between the bodies, a classical formulation of unilateral contact constraints is given by the Signorini conditions expressed at the displacement level [57]: the non-penetration constraint $d_N \ge 0$, the non-traction constraint $F_N \ge 0$ and $d_N F_N = 0$, enforcing the fact that both constraints cannot be active simultaneously.

A large majority of works and industrial tools simplify the Signorini conditions by regularization techniques [12, 38, 21]. However, the present work aims to strictly enforce the non-smooth law and to prevent non-physical and residual penetration. A result, preference is given to solution techniques involving constraints management and Lagrange multipliers. A widely known procedure allowing to solve dynamics with contact occurrence in time domain is the Carpenter algorithm [10]. Various published works [46, 4, 36] selected the Carpenter algorithm for its simplicity of implementation while giving results in agreement with experiments [53]. Note that other solutions methods with Lagrange multipliers are also available such as the Moreau-Jean [39, 32] or the Paoli-Schatzman [42] algorithms. They include the non-smooth Newton impact law [2] in the Signorini conditions.

Friction is accounted for by introducing a force $\mathbf{F}_f$ oriented in the direction opposed to the relative speed $\mathbf{v}_T$ at the contact interface between the interacting bodies. Its magnitude is commonly dictated by the Coulomb dry friction model stating that the friction force must satisfy

$$ \|\mathbf{F}_f\| = \mu F_N, \quad \text{if} \quad \mathbf{v}_T \neq \mathbf{0} \qquad (31a) $$

$$ \|\mathbf{F}_f\| \le \mu F_N, \quad \text{if} \quad \mathbf{v}_T = \mathbf{0} \qquad (31b) $$

Besides, the implementation of Coulomb's friction law requires careful attention because of the non-smoothness of the sticking friction phase governed by Equation (31b). If one desires to strictly enforce the Coulomb friction law as written in Equations (31a) and (31b), the exact approach for frictional contact treatment is based on mathematical concepts from convex analysis and defines contact efforts in the Coulomb cone [2]. Otherwise, it is possible to regularize the Coulomb friction law [12, 18] or to consider a sliding-only friction law characterized by

$$ F_f = \mathbf{F}_f \cdot \mathbf{t} = -\mu F_N \text{sign}(\mathbf{v}_T) \quad \text{with} \quad \mathbf{v}_T = \mathbf{v}_T \cdot \mathbf{t}. \qquad (32) $$

where $\mathbf{t}$ is the unit tangential vector in local contact frame. The above relation implies that the relative velocity is known and its assessment is further developed in the present section. This formulation is commonly considered when the rotor operates at high rotational speeds [36, 4, 46].

In [30], Carpenter and Moreau-Jean algorithms were compared on rotor-stator rubbing between rigid bodies with various friction models. Good agreement is reported and pure rolling motion (i.e. sticking friction) can take place during the simulation. More interestingly, rotational velocity transient was mainly affected during sliding friction phases. This explains why sliding-only friction is considered in the present model, see Equation (32). The dry-whirl response (i.e. pure rolling motion) can still be captured but corresponds to an almost vanishing sliding velocity with rapid sign changes [31]. This friction model is simple to implement but keeps physical limitations for the dry-whirl modeling. Chattering of sliding speed takes place during dry whirl because the magnitude of friction force is ignores Equation (31b)
<figure><img src="image_7.png" /><figcaption>Figure 6: Computation of the nodal distance between de rotor center and the 4th node on the stator (j = 4)</figcaption></figure>

The rotational velocity at time step i is a function of $q_{i+1}$ because of the time discretization

$$ \dot{\theta}_i(\mathbf{q}_{i+1}) \approx \frac{\theta_{i+1} - \theta_{i-1}}{2h} \qquad (35) $$

along with the fact that $\theta_{i+1}$ is the fifth element of $q_{i+1}$, see Equations (2) and (29). Strictly speaking, matrices $\bar{M}$, $\bar{D}$, $\bar{K}$ and scalar $\dot{\theta}_i$ also depend on $q_i$ and $q_{i-1}$. However, these vectors are already known and merely imply an update of the system at each time step.

The modified version of Carpenter's algorithm [26] performs the following steps.

**Step 1: prediction** The objective is to find a constraint-free candidate $q_{i+1}^p$, where exponent $p$ stands for predicted. Solution to Equation (33) is performed by calling the Python built-in nonlinear solver `fsolve` that requires an initial guess, chosen as $q_i$. This function is a FORTRAN wrapper of a modified version of the Powell hybrid method [10] in order to find a root close to the initial guess. The computational time is reduced by providing the Jacobian matrix of Equation (33)

$$ \mathbf{J}(\mathbf{q}_{i+1}) = \left( \frac{\mathbf{M}(\mathbf{q}_i)}{h^2} + \beta_1 \frac{\mathbf{D}}{h} \right) + \frac{\beta_1}{2h^2} (\theta_{i+1} - \theta_{i-1})\mathbf{G} + \frac{\mathbf{G}}{2h^2} (\beta_1 \mathbf{q}_{i+1} - (2\beta_1 - 1)\mathbf{q}_i + (\beta_1 - 1)\mathbf{q}_{i-1}), \quad (36) $$

which is also function of $q_{i+1}$. Please note that the last term in Equation (36) is a column vector to be added on the 5th column of $\mathbf{J}$ (i.e. the column corresponding to dof $\theta_{i+1}$ of the turbine).

**Step 2: penetration check** The approach selected is to evaluate the nodal rotor-stator distances

$$ d_N^j(\mathbf{q}) = \delta_N^j(\mathbf{q}) - R_r - a \quad \text{for } j = 0, \dots, N_{e-1} \qquad (37) $$

with

$$ \delta_N^j(\mathbf{q}) = \sqrt{(x_{sj}(\mathbf{q}) - x_r)^2 + (y_{sj}(\mathbf{q}) - y_r)^2}. \qquad (38) $$

The couple $(x_{sj}, y_{sj})$ embeds the $j^{th}$ stator node coordinates in the fixed frame $(O, X, Y)$ as shown in Figures 5 and 6. Their expression is given in Equation (17) where $q_s$ can be easily replaced by $q$ as function variable based on Equation (29), and $(x_r, y_r)$ are the rotor displacements already defined in Figure 3. Nodes that violate the Signorini conditions are grouped in the set of active constraints denoted $\mathcal{A} = \{j \in [0, N_e - 1] | d_N^j < 0\}$. In the rest of the procedure, the index of node j in the set $\mathcal{A}$ is noted k and goes from 1 to $a = |\mathcal{A}|$, cardinal of $\mathcal{A}$. The associated negative penetrations are collected in the vector $g_{i+1}^p$.

**Step 3: correction** If $\mathcal{A}$ is empty, there is no correction to apply and discrete time step is incremented. In contrast, when $\mathcal{A}$ contains $a$ nodes, $a$ Lagrange multipliers $\lambda^{(k)}$ ($k = 1, \dots, a$) should be computed. Note that these multipliers correspond to the normal unilateral contact forces $F_N^{(k)}$ when solution is converged. The computation of the resultant contact effort implies the definition of the following vectors for each active constraint $k$, involving node j:
from Coulomb law. Therefore, $F_f$ from Equation (32) can only take two values $-\mu F_N$ or $\mu F_N$, while Coulomb law indicates that $F_f \in [-\mu F_N, \mu F_N]$. In the present paper, contact occurrences are handled with Carpenter's algorithm.

All unknowns are grouped into the vector $\mathbf{q} = (\mathbf{q}_r, \mathbf{q}_s, \mathbf{T}_s)$ and global mass, damping, gyroscopic and stiffness matrices are named $\mathbf{M}(\mathbf{q}), \mathbf{D}, \mathbf{G}$, and $\mathbf{K}$, respectively, such that Equation (29) becomes

$$ \mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} + (\mathbf{D} + \dot{\theta}\mathbf{G})\dot{\mathbf{q}} + \mathbf{K}\mathbf{q} = \mathbf{f}^{\text{ext}} + \mathbf{f}_c. \quad (33) $$

Dependency of the mass matrix on the turbine dofs is shown to stress that it is not constant.

## 3.2. Modified Carpenter algorithm

In Carpenter algorithm, Equation (33) is discretized in time using an integrator from the Newmark family [27]. In order to avoid confusion with other symbols, numerical damping and stability parameters of the Newmark schemes are respectively denoted $\beta_1$ and $\beta_2$. Simulation time is divided into $N$ intervals of duration $h$. For studies involving mechanics only, the centered finite-difference scheme with $\beta_1 = 0.5$ and $\beta_2 = 0$ is appropriate but Thorin & al. [52] showed that $\beta_1$ must be slightly offset because of stability issues to solve the heat equation with thermomechanical coupling. Consequently, the integrator is chosen in this work with $\beta_1 = 0.501$ and $\beta_2 = 0$. In the present paper, to reduce the amount of notation, it is decided to reflect dependencies regarding unknown quantities only. In discrete time, the unknown vector is $\mathbf{q}_{i+1}$ at time step $i + 1$ and the contact-free equations of motion become

$$ \bar{\mathbf{M}}(\mathbf{q}_{i+1})\mathbf{q}_{i+1} + \bar{\mathbf{D}}(\mathbf{q}_{i+1})\mathbf{q}_i + \bar{\mathbf{K}}(\mathbf{q}_{i+1})\mathbf{q}_{i-1} = \bar{\mathbf{f}}_i, \quad i = 1, \dots, N \quad (34) $$

where

$$ \begin{align} \bar{\mathbf{M}}(\mathbf{q}_{i+1}) &:= \frac{\mathbf{M}(\mathbf{q}_i)}{h^2} + \beta_1 \frac{\mathbf{D} + \dot{\theta}_i(\mathbf{q}_{i+1})\mathbf{G}}{h} \nonumber \\ \bar{\mathbf{D}}(\mathbf{q}_{i+1}) &:= \frac{2\mathbf{M}(\mathbf{q}_i)}{h^2} + \frac{2\beta_1 - 1}{h} (\mathbf{D} + \dot{\theta}_i(\mathbf{q}_{i+1})\mathbf{G}) - (0.5 + \beta_1)\mathbf{K} \tag{35} \\ \bar{\mathbf{K}}(\mathbf{q}_{i+1}) &:= \frac{\mathbf{M}(\mathbf{q}_i)}{h^2} + \frac{\beta_1 - 1}{h} (\mathbf{D} + \dot{\theta}_i(\mathbf{q}_{i+1})\mathbf{G}) + (0.5 - \beta_1)\mathbf{K} \nonumber \\ \bar{\mathbf{f}}_i &:= (1 - \beta_1)\mathbf{f}_{i-1}^{\text{ext}} + (1 + \beta_1)\mathbf{f}_i^{\text{ext}}. \nonumber \end{align} $$

The rotational velocity at time step $i$ is a function of $\mathbf{q}_{i+1}$ because of the time discretization

$$ \dot{\theta}_i(\mathbf{q}_{i+1}) \approx \frac{\theta_{i+1} - \theta_{i-1}}{2h} \quad (36) $$

along with the fact that $\theta_{i+1}$ is the fifth element of $\mathbf{q}_{i+1}$ (see Equations (2) and (30)). Strictly speaking, matrices $\bar{\mathbf{M}}, \bar{\mathbf{D}}, \bar{\mathbf{K}}$ and scalar $\dot{\theta}_i$ also depend on $\mathbf{q}_i$ and $\mathbf{q}_{i-1}$. However, these vectors are already known and merely imply an update of the system at each time step.

The modified version of Carpenter's algorithm [30] performs the following steps.

**Step 1: prediction** The objective is to find a constraint-free candidate $\mathbf{q}_{i+1}^p$, where exponent $p$ stands for predicted. Solution to Equation (34) is performed by calling the Python built-in nonlinear solver `fsolve` that requires an initial guess, chosen as $\mathbf{q}_i$. This function is a FORTRAN wrapper of a modified version of the Powell hybrid method [11] in order to find a root close to the initial guess. Computation time is reduced by providing the Jacobian matrix of Equation (34)

$$ \mathbf{J}(\mathbf{q}_{i+1}) = \left( \frac{\mathbf{M}(\mathbf{q}_i)}{h^2} + \beta_1 \frac{\mathbf{D}}{h} \right) + \frac{\beta_1}{2h^2} (\theta_{i+1} - \theta_{i-1})\mathbf{G} + \frac{\mathbf{G}}{2h^2} (\beta_1 \mathbf{q}_{i+1} - (2\beta_1 - 1)\mathbf{q}_i + (\beta_1 - 1)\mathbf{q}_{i-1}), \quad (37) $$

which is also function of $\mathbf{q}_{i+1}$. Please note that the last term in Equation (37) is a column vector to be added on the 5th column of $\mathbf{J}$ (i.e. the column corresponding to dof $\theta_{i+1}$ of the turbine).
* Normal direction components $C_N^{(k)}$ which corresponds to the gradient of normal distance $\nabla_q d_N^j(q)$ evaluated for predicted displacement $\mathbf{q}_{i+1}^p$. Applying the gradient operator to Equation (37) yields

$$C_N^{(k)} = \nabla_q \delta_N^j(\mathbf{q}_{i+1}^p) \qquad (39)$$

Therefore $C_N^{(k)}$ must be computed at each time step $i + 1$ but index is omitted for simplicity.

* The components $C_T^{(k)}$ in the tangential direction are obtained in a similar way to $C_N^{(k)}$, however by deriving a tangential distance $d_T$, which writes for node $j$

$$d_T^j(\mathbf{q}) = \delta_N^j \frac{y_{sj}(\mathbf{q}) - y_r}{x_{sj}(\mathbf{q}) - x_r} - R_r \theta. \qquad (40)$$

In the expression above, $\delta_N^j$ is considered constant in the gradient computation. The analytic expressions for $C_N^{(k)}$ and $C_T^{(k)}$ are given in Appendix B.

For sliding-only friction, each multiplier contributes to contact efforts according to the Coulomb friction law

$$\mathbf{f}_{c,i+1}^q = \sum_{k=1}^{A} (C_N^{(k)} - \mu \operatorname{sign}(v_{T,i}^{(k)}) C_T^{(k)}) \lambda_{i+1}^{(k)} \qquad (41)$$

where $v_{T,i}^{(k)}$, in the $k^{th}$ active constraint, defines the tangential component of the relative speed in the frame of node $j \in \mathbb{A}$. Based on the predicted speed vector $\mathbf{v}_i^p = (\mathbf{q}_{i+1}^p - \mathbf{q}_i)/h$, the sliding velocity is updated at each time step $i$ using

$$v_{T,i}^{(k)} = C_T^{(k)\top} \mathbf{v}_i^p \qquad (42)$$

but index $i$ is dropped to simplify notations. One should notice that each column vector $C_N^{(k)}$ or $C_T^{(k)}$ is almost empty. Only degrees-of-freedom of the rotor and of node $j$ belonging to the stator are affected by the $k^{th}$ constraint. In discrete time, the $a$ multipliers are collected in the vector $\Lambda_{i+1}$. A matrix $C_{NT}$ may also be defined with $a$ columns, where each column equals $C_N^{(k)} - \mu \operatorname{sign}(v_T^{(k)}) C_T^{(k)}$ such that Equation (41) becomes

$$\mathbf{f}_{c,i+1}^q = C_{NT} \Lambda_{i+1}. \qquad (43)$$

The same approach is achieved for the heat power vector $\dot{\mathbf{Q}}_s$, see Equation (28). Based on Equation (24), on the right member of Equation (25) and knowing that there exist heat sources for the stator only at the $k$ activated nodes $j \in \mathbb{A}$, we have

$$\dot{\mathbf{Q}}_{s,i+1}^{(k)} = \frac{1}{2} \mu |v_T^{(k)}| \lambda_{i+1}^{(k)}. \qquad (44)$$

The above expression can be included in Carpenter's algorithm and for each activated constraint $k$ is assigned a heat contribution vector $C_{\theta}^{(k)}$ that is zero except on its $j^{th}$ component

$$C_{\theta_j}^{(k)} = \frac{1}{2} \mu |v_T^{(k)}| \qquad (45)$$

resulting in $\dot{\mathbf{Q}}_{s,i+1} = C_{\theta} \Lambda_{i+1}$. Eventually, the contributions of Lagrange multipliers to both physics are grouped in the matrix $C_{NT\theta} = (C_{NT}, C_{\theta})^\top$. Therefore, the contribution $\mathbf{f}_c$ to the dynamics and heat equations, see Equation (32), is

$$\mathbf{f}_{c,i+1} = C_{NT\theta} \Lambda_{i+1}. \qquad (46)$$

They are computed at step $i + 1$ based on the relation

$$\Lambda_{i+1} = -(C_N^T \bar{\mathbf{M}}^{-1}(\mathbf{q}_{i+1}^p) C_{NT\theta})^{-1} \mathbf{g}_{i+1}^p \qquad (47)$$

It should be noted that the Lagrange multipliers are computed based on the simplifying assumption that the rotational velocity used to compute gyroscopic effects is the predicted one, that is $\dot{\theta}_i^p$ which is an entry of $\mathbf{v}_i^p$. Displacement and temperature encompassed in vector $\mathbf{q}_{i+1}$, see Equation (29), are corrected via

$$\mathbf{q}_{i+1} = \mathbf{q}_{i+1}^p + \bar{\mathbf{M}}(\mathbf{q}_{i+1}^p)^{-1} C_{NT\theta} \Lambda_{i+1} \qquad (48)$$

and eventually discrete time is incremented.

The above procedure is summarized in Algorithm 1 in Appendix C.
<figure><img src="image_7.png" /><figcaption>**Figure 6:** Computation of the nodal distance between de rotor center and the 4th node on the stator (j = 4)</figcaption></figure>

**Step 2: penetration check** The approach selected is to evaluate the nodal rotor-stator distances

$$d_{N}^{j}(\mathbf{q}) = \delta_{N}^{j}(\mathbf{q}) - R_{r} - a \quad \text{for } j = 0, \dots, N_{e-1} \qquad (38)$$

with

$$\delta_{N}^{j}(\mathbf{q}) = \sqrt{(x_{sj}(\mathbf{q}) - x_{r})^{2} + (y_{sj}(\mathbf{q}) - y_{r})^{2}}. \qquad (39)$$

The couple $(x_{sj}, y_{sj})$ embeds the jth stator node coordinates in the fixed frame $(O, X, Y)$ as shown in Figures 5 and 6. Their expression is given in Equation (18) where $q_s$ can be easily replaced by $q$ as function variable based on Equation (30), and $(x_r, y_r)$ are the rotor displacements already defined in Figure 3. Nodes that violate the Signorini conditions are grouped in the set of active constraints denoted $\mathcal{A} = \{j \in [0, N_e - 1] | d_N^j < 0\}$. In the rest of the procedure, the index of node j in the set $\mathcal{A}$ is noted k and goes from 1 to $a = |\mathcal{A}|$, cardinal of $\mathcal{A}$. The associated negative penetrations are collected in the vector $g_{i+1}^p$.

**Step 3: correction** If $\mathcal{A}$ is empty, there is no correction to apply and discrete time step is incremented. In contrast, when $\mathcal{A}$ contains $a$ nodes, $a$ Lagrange multipliers $\lambda^{(k)}$ ($k = 1, \dots, a$) should be computed. Note that these multipliers correspond to the normal unilateral contact forces $F_N^{(k)}$ when solution is converged. The computation of the resultant contact effort implies the definition of the following vectors for each active constraint $k$, involving node $j$:

* Normal direction components $C_N^{(k)}$ which corresponds to the gradient of normal distance $\nabla_q d_N^j(\mathbf{q})$ evaluated for predicted displacement $q_{i+1}^p$. Applying the gradient operator to Equation (38) yields

$$C_N^{(k)} = \nabla_q \delta_N^j (q_{i+1}^p) \qquad (40)$$

Therefore $C_N^{(k)}$ must be computed at each time step $i + 1$ but index is omitted for simplicity.

* The components $C_T^{(k)}$ in the tangential direction are obtained in a similar way to $C_N^{(k)}$, however by deriving a tangential distance $d_T$, which writes for node $j$

$$d_T^j(\mathbf{q}) = \delta_N^j \frac{y_{sj}(\mathbf{q}) - y_r}{x_{sj}(\mathbf{q}) - x_r} - R_r \theta. \qquad (41)$$

In the expression above, $\delta_N^j$ is considered constant in the gradient computation. The analytic expressions for $C_N^{(k)}$ and $C_N^{(k)}$ are given in Appendix B.

For sliding-only friction, each multiplier contributes to contact efforts according to the Coulomb friction law

$$f_{c,i+1}^q = \sum_{k=1}^{A} (C_N^{(k)} - \mu \operatorname{sign}(v_{T,i}^{(k)}) C_T^{(k)}) \lambda_{i+1}^{(k)} \qquad (42)$$
### 3.3. Choice of time step

One last remark is about the time step $h$ which should satisfy the (linear) stability criterion

$$h < \frac{\sqrt{2}}{\max(2\pi f_n)\sqrt{\beta_1 - \beta_2}} \qquad (49)$$

based on the eigenvalues of the system without damping and the parameters of the Newmark integrator. It is clear that the time scale is governed by the eigenfrequencies of the stator and Equation (49) with 36 nodes is respected for $h \le 4.5 \times 10^{-7}$ s. However, Equation (49) is restricted to linear domain and stability proofs exist for contact algorithms only in particular cases. Simulations were tested with $h = 4 \times 10^{-7}$ s and showed converged results, as detailed in the following section.

## 4. Results

### 4.1. Reference case

The reference case involves a ring connected to the casing through 12 screws. The friction coefficient is set to $\mu = 0.1$ and simulations are run until $t_{ref} = 300$ ms. The largest rotational velocity $\Omega_{max} = \max_{t \in [0, t_{ref}]}(\Omega(t))$ on the considered time interval is of primary interest in this research. It is desired that $\Omega_{max} < \Omega_{lim}$, where $\Omega_{lim}$ is a speed limit not to be exceeded. A convergence analysis is performed and speed transients for 36 and 72 nodes are displayed in Figure 7(a). It turns out that for both cases, the rotational velocity evolution is similar and $\Omega_{max}$ remains bounded by $\Omega_{lim}$.

The frictional torque presented in Figure 7(b) sometimes vanishes at the beginning of the simulation, meaning that contact is open. Indeed, the frictional torque is computed using

$$\tau_f(t) = R_r \sum_{j=1}^{N_e} \mu \lambda_j(t) \qquad (50)$$

where a vanishing $\lambda_j(t)$ implies that node $j$ does not violate the non-penetration constraint (i.e. $j \notin A$).

<figure><img src="image_8.png" /><figcaption>**Figure 7:** Transient response for the turbine rotation: 36 nodes (—) and 72 nodes (- -). Rotational velocity remains bounded since frictional torque overcomes aerodynamic torque.</figcaption></figure>

Moreover, it becomes clear that the frictional torque progressively overcomes the driving aerodynamic torque, which explains the rotational speed decrease. One also notices that the gap transient at node 9, i.e. the 4th screw in Figure 8(a), has the same profile for both meshes. The same remark holds for the rotor vibration, which is illustrated along the x-direction in Figure 8(b). The corrected gap never becomes negative which confirms the absence of residual penetration and that the Signorini conditions are correctly satisfied. After about 20 ms, vibration magnitudes and gap start to slowly shift between the meshes but this lag seems to be limited. Consequently, it is considered that the nonlinear dynamics is correctly captured with 36 nodes only. The question is now to identify the origin of this torque increase. Precessional speed is thus investigated and is defined here as the tangential speed of the rotor. It implies that the vector $v_{prec}$ is aligned with unit tangential vector $t$, and that its component $v_{prec}$ along $t$ is, see Figure 9(b)

$$v_{prec} \approx \frac{x_r \dot{y}_r - y_r \dot{x}_r}{\sqrt{x_r^2 + y_r^2}} \qquad (51)$$
where $v_{T,i}^{(k)}$, in the $k^{th}$ active constraint, defines the tangential component of the relative speed in the frame of node $j \in A$. Based on the predicted speed vector $v_i^p = (q_{i+1}^p - q_i)/h$, the sliding velocity is updated at each time step $i$ using

$$v_{T,i}^{(k)} = C_{T}^{(k)\top} v_{i}^{p} \qquad (43)$$

but index $i$ is dropped to simplify notations. One should notice that each column vector $C_N^{(k)}$ or $C_T^{(k)}$ is almost empty. Only degrees-of-freedom of the rotor and of node $j$ belonging to the stator are affected by the $k^{th}$ constraint. In discrete time, the $a$ multipliers are collected in the vector $\Lambda_{i+1}$. A matrix $C_{NT}$ may also be defined with $a$ columns, where each column equals $C_N^{(k)} - \mu \text{sign}(v_T^{(k)})C_T^{(k)}$ such that Equation (42) becomes

$$\mathbf{f}_{c,i+1}^{q} = \mathbf{C}_{\text{NT}} \mathbf{\Lambda}_{i+1}. \qquad (44)$$

The same approach is achieved for the heat power vector $\dot{\mathbf{Q}}_s$, see Equation (29). Based on Equation (25), on the right member of Equation (26) and knowing that there exist heat sources for the stator only at the $k$ activated nodes $j \in A$, we have

$$\dot{\mathbf{Q}}_{s,i+1}^{(k)} = \frac{1}{2} \mu |v_T^{(k)}| \lambda_{i+1}^{(k)}. \qquad (45)$$

The above expression can be included in Carpenter's algorithm and for each activated constraint $k$ is assigned a heat contribution vector $C_{\theta}^{(k)}$ that is zero except on its $j^{th}$ component

$$C_{\theta_j}^{(k)} = \frac{1}{2} \mu |v_T^{(k)}| \qquad (46)$$

resulting in $\dot{\mathbf{Q}}_{s,i+1} = \mathbf{C}_{\theta} \mathbf{\Lambda}_{i+1}$. Eventually, the contributions of Lagrange multipliers to both physics are grouped in the matrix $\mathbf{C}_{NT\theta} = (\mathbf{C}_{NT}, \mathbf{C}_{\theta})^{\top}$. Therefore, the contribution $\mathbf{f}_c$ to the dynamics and heat equations, see Equation (33), is

$$\mathbf{f}_{c,i+1} = \mathbf{C}_{\text{NT}\theta} \mathbf{\Lambda}_{i+1}. \qquad (47)$$

They are computed at step $i + 1$ based on the relation

$$\mathbf{\Lambda}_{i+1} = -(\mathbf{C}_{\text{N}}^{\top} \bar{\mathbf{M}}^{-1} (\mathbf{q}_{i+1}^{p}) \mathbf{C}_{\text{NT}\theta})^{-1} \mathbf{g}_{i+1}^{p} \qquad (48)$$

It should be noted that the Lagrange multipliers are computed based on the simplifying assumption that the rotational velocity used to compute gyroscopic effects is the predicted one, that is $\dot{\theta}_i^p$ which is an entry of $v_i^p$. Displacement and temperature encompassed in vector $\mathbf{q}_{i+1}$ (cf. Equation (30)) are corrected via

$$\mathbf{q}_{i+1} = \mathbf{q}_{i+1}^{p} + \bar{\mathbf{M}}(\mathbf{q}_{i+1}^{p})^{-1} \mathbf{C}_{\text{NT}\theta} \mathbf{\Lambda}_{i+1} \qquad (49)$$

and eventually discrete time is incremented.

The above procedure is summarized in Algorithm 1 in Appendix C.

### 3.3. Choice of time step

One last remark is about the time step $h$ which should satisfy the (linear) stability criterion

$$h < \frac{\sqrt{2}}{\max(2\pi f_n)\sqrt{\beta_1 - \beta_2}} \qquad (50)$$

based on the eigenvalues of the system without damping and the parameters of the Newmark integrator. It is clear that the time scale is governed by the eigenfrequencies of the stator and Equation (50) with 36 nodes is respected for $h \le 4.5 \times 10^{-7}$ s. However, Equation (50) is restricted to linear domain and stability proofs exist for contact algorithms only in particular cases. Simulations were tested with $h = 4 \times 10^{-7}$ s and showed converged results, as detailed in the following section.
<figure><img src="image_9.png" /><figcaption>(a) Gap function at node 9 (4th screw)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Rotor horizontal displacement</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(a) Transient response</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Definition by projection of rotor speed on unit tangential vector</figcaption></figure>

<caption>**Figure 9:** Rotor precessional speed: forward precessional motion at the very beginning of the simulation and then negative speed on average, reflecting partial contact in backward precession.</caption>

This precessional speed tends to decrease with a negative mean value in Figure 9(a) reflecting that backward precessional motion becomes dominant. Hence, the response for the reference case is characterized by a partial rub with backward precessional motion. This behavior is known to come with higher amplitude of the rebounds [21] which is confirmed by the slight increase of the vibration amplitude of the rotor with time, leading to higher contact efforts. However, contact does not become fully annular and dry whip is not initiated in the simulation time for the reference case.

The main difference between the two meshes lies in the thermodynamics response. Both cases give the same time evolution for the frictional energy, as shown in Figure 10, computed according to

$$Q_f(t) = \sum_{j=0}^{N_e-1} \int_{0}^{t} \mu |v_{T,j}(\tau)| \lambda_j(\tau) d\tau. \qquad (52)$$

However, the temperature profiles at the end of the simulation differ. The configuration with 36 elements shows negative temperature for nodes that are not screwed, which is non-physical since the frictional power at these nodes is positive, so that there are always heat sources and not sinks. This is a purely numerical phenomenon related to mesh refinement. With $N_e = 72$, the temperature at the end of the simulation and far from screw nodes lies between 689 °C and 1055 °C with 904 °C on average. For a stator in stainless steel, the melting point is at $T_m = 1540$ °C. For zero initial temperature, it implies that the stator does not melt excepting at the screw points. Assuming that the initial temperature was around 300 °C, which is a realistic value for a gas turbine in nominal operation, the above range would represent a temperature elevation. However, the stator temperature away from the screws would remain under the melting point. Most of the contact efforts (and so the attendant frictional energy) are concentrated at the screws as confirmed by the Lagrange multipliers for nodes 10, 11 (no screw) and 9 (screw) in Figure 11. As the precessional speed lies between -2 m/s and 2 m/s, it leads to a precessional frequency which could exceed 400 Hz. As a consequence, the thermal influence zone around the screws is very narrow and if the
<figure><img src="image_8.png" /><figcaption>Figure 7: Transient response for the turbine rotation: 36 nodes (—) and 72 nodes (---). Rotational velocity remains bounded since frictional torque overcomes aerodynamic torque.</figcaption></figure>

# 4. Results

## 4.1. Reference case

The reference case involves a ring connected to the casing through 12 screws. The friction coefficient is set to $\mu = 0.1$ and simulations are run until $t_{ref} = 300$ ms. The largest rotational velocity $\Omega_{max} = \max_{t \in [0, t_{ref}]} (\Omega(t))$ on the considered time interval is of primary interest in this research. It is desired that $\Omega_{max} < \Omega_{lim}$, where $\Omega_{lim}$ is a speed limit not to be exceeded. A convergence analysis is performed and speed transients for 36 and 72 nodes are displayed in Figure 7(a). It turns out that for both cases, the rotational velocity evolution is similar and $\Omega_{max}$ remains bounded by $\Omega_{lim}$.

The frictional torque presented in Figure 7(b) sometimes vanishes at the beginning of the simulation, meaning that contact is open.

Indeed, the frictional torque is computed using

$$ \tau_f(t) = R_r \sum_{j=1}^{N_e} \mu \lambda_j(t) \qquad (51) $$

where a vanishing $\lambda_j(t)$ implies that node j does not violate the non-penetration constraint (i.e. $j \notin A$).

Moreover, it becomes clear that the frictional torque progressively overcomes the driving aerodynamic torque, which explains the rotational speed decrease. One also notices that the gap transient at node 9, i.e. the 4th screw in Figure 8(a), has the same profile for both meshes. The same remark holds for the rotor vibration, which is illustrated along the x-direction in Figure 8(b). The corrected gap never becomes

<figure><img src="image_9.png" /><figcaption>Figure 8: Convergence study of structural vibration: 36 nodes (—) and 72 nodes (---). Solution is converged for $t < 30$ ms. Then, a constant phase shift takes place in the rest of the simulation.</figcaption></figure>

negative which confirms the absence of residual penetration and that the Signorini conditions are correctly satisfied. After about 20 ms, vibration magnitudes and gap start to slowly shift between the meshes but this lag seems to be limited. Consequently, it is considered that the nonlinear dynamics is correctly captured with 36 nodes only. The question is now to identify the origin of this torque increase. Precessional speed
<figure><img src="image_13.png" /><figcaption>**Figure 10:** Thermal response: 36 nodes (—), 72 nodes (---). Total frictional energy is identical with mesh refinement but temperature distribution shows extremely local peaks at screw nodes.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11:** Evidence of contact efforts and frictional energy concentration at screw nodes with $N_e = 36$: $j = 9$ (4th screw) (—), $j = 10$ (---) and $j = 11$ (----).</figcaption></figure>

mesh is not sufficiently refined within this layer, results may be non-physical. The issue can be solved by increasing the number of mesh elements and the configuration with 72 nodes case indeed provides more reliable information.

## 4.2. Sensitivity to the friction coefficient

As mentioned in the literature review, multiple mechanical parameters are known to influence the emergence of backward precessional motions, one of them being the friction coefficient. One seeks to investigate the variation of the rotational speed max values when the friction coefficient is increased. Speed transients in Figure 12 confirm that $\Omega_{max}$ decreases with $\mu$. However, $\mu$ must be chosen carefully since an abrupt

<figure><img src="image_15.png" /><figcaption>**Figure 12:** Sensitivity of the rotational speed transient to the friction coefficient: $\mu = 0.1$ (ref) (—), $\mu = 0.12$ (---), $\mu = 0.15$ (----), and $\mu = 0.2$ (---).</figcaption></figure>

decrease is obtained for $\mu \ge 0.15$. This phenomenon reflects that *dry whip* is taking place, i.e. fully annular backward precessional motion. The dynamical response is further investigated for $\mu = 0.15$ with
is thus investigated and is defined here as the tangential speed of the rotor. It implies that the vector $v_{prec}$ is aligned with unit tangential vector $t$, and that its component $v_{prec}$ along $t$ is (see Figure 9(b))

$$v_{prec} \approx \frac{x_r \dot{y}_r - y_r \dot{x}_r}{\sqrt{x_r^2 + y_r^2}}. \qquad (52)$$

This precessional speed tends to decrease with a negative mean value in Figure 9(a) reflecting that backward precessional motion becomes dominant. Hence, the response for the reference case is characterized by a

<figure><img src="image_10.png" /><figcaption>**Figure 9:** Rotor precessional speed: forward precessional motion at the very beginning of the simulation and then negative speed on average, reflecting partial contact in backward precession.</figcaption></figure>

partial rub with backward precessional motion. This behavior is known to come with higher amplitude of the rebounds [23] which is confirmed by the slight increase of the vibration amplitude of the rotor with time, leading to higher contact efforts. However, contact does not become fully annular and dry whip is not initiated in the simulation time for the reference case.

The main difference between the two meshes lies in the thermodynamics response. Both cases give the same time evolution for the frictional energy, as shown in Figure 10, computed according to

$$Q_f(t) = \sum_{j=0}^{N_e-1} \int_{0}^{t} \mu |v_{T,j}(\tau)| \lambda_j(\tau) d\tau. \qquad (53)$$

However, the temperature profiles at the end of the simulation differ. The configuration with 36 elements

<figure><img src="image_11.png" /><figcaption>**Figure 10:** Thermal response: 36 nodes (—), 72 nodes (- -). Total frictional energy is identical with mesh refinement but temperature distribution shows extremely local peaks at screw nodes.</figcaption></figure>

shows negative temperature for nodes that are not screwed, which is non-physical since the frictional power at these nodes is positive, so that there are always heat sources and not sinks. This is a purely numerical phenomenon related to mesh refinement. Indeed, with $N_e = 72$, the nodal temperatures are always positive but their accuracy must still be discussed. Towards the end of the simulation, the temperature at screw nodes exceeds 10000° C which is unrealistic. High values are expected since contact efforts (and so
mesh refinement. It turns out that $\Omega_{max}$ is identical for all mesh sizes, in Figure 13(a), and velocity profiles are similar until $t = 135$ ms. Again, perfect match is predicted for $t < 30$ ms for the rotor displacement $x_r(t)$ in Figure 13(b). Then, a phase lag takes place and dry whip is initiated in a time frame between 135 and 160 ms, which is more accurate than for a model with rigid stator and rotor. Fully annular contact is

<figure><img src="image_16.png" /><figcaption>**Figure 13:** Mesh convergence analysis for $\mu = 0.15$: $N_e = 36$ (—), $N_e = 48$ (—), $N_e = 60$ (—).</figcaption></figure>

confirmed by frictional torque that never vanishes after 150 ms.

<figure><img src="image_17.png" /><figcaption>**Figure 14:** Dynamical response for $\mu = 0.15$ and $N_e = 36$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 15:** Precessional velocity for dry whip: $N_e = 36$ and $\mu = 0.15$.</figcaption></figure>

Violent dry whip also emerges with an abrupt increase of the precessional velocity in absolute value as shown in Figure 15. In other words, kinetic energy is transferred from rotation to vibration through frictional mechanisms. Let us consider the front view of a rotorshaft with circular section depicted in Figure 16. In case of backward precessional motion, in Figure 16(b), the friction force $F_f$ acts against $\Omega$ but in favour of $v_{prec}$. The consequence is that the rotor displacement, Figure 13(b), is amplified by a factor 10. Figure 17 also shows that the stator suffers significant deformation that could indicate a response out of the elastic domain for material properties. However, once $t > 160$ ms, frictional torque, presented in Figure 14(a) for $N_e = 36$, tends exponentially to zero and precessional speed stabilizes at -20 m/s. This
the attendant frictional energy) are mostly concentrated at the screws, as confirmed by the Lagrange multipliers for nodes 10, 11 (no screw) and 9 (screw) in Figure 11. However, if discretizing parts A and B by springs has little effect on the structural dynamics, it does introduce singularities in the heat equation. These physical phenomena operate on different time and length scales. Firstly, heat sources are localized to individual nodes, and the thermal boundary conditions set by the screws are not fully realistic. This represents a limitation of the academic model: the temperature at the screw level cannot be precisely quantified. Despite its limitations, the model offers valuable insights expected from a phenomenological approach, highlighting areas most sensitive to high temperatures and most susceptible to melt wear. Notably, with precessional speeds ranging from -2 m/s and 2 m/s, the heat source frequency exceeds 400 Hz. Such a high frequency suggests that thermal information evolves too rapidly to penetrate deeply into the beam. This is indeed the case: the thermal influence zone around the screws remains very narrow, and without sufficient mesh refinement in this layer, results can become unrealistic. This issue can be addressed by increasing the number of mesh elements, as seen in the configuration with 72 nodes, which yields more reliable information far from the screw nodes. The temperature there lies between 689 °C and 1055 °C with 904 °C on average. For a stator in stainless steel, the melting point is at $T_m = 1540$ °C. For zero initial temperature, it implies that the stator does not melt excepting at the screw points. To go further, it is realistic to assume initial temperature around 300 °C for a gas turbine in nominal operation. The above range would then represent a temperature elevation but the stator temperature away from the screws would still remain under the melting point.

<figure><img src="image_12.png" /><figcaption>**Figure 11:** Evidence of contact efforts and frictional energy concentration at screw nodes with $N_e = 36$: $j = 9$ (4th screw) (—), $j = 10$ (—) and $j = 11$ (---).</figcaption></figure>

## 4.2. Sensitivity to the friction coefficient

As mentioned in the literature review, multiple mechanical parameters are known to influence the emergence of backward precessional motions, one of them being the friction coefficient. One seeks to investigate the variation of the rotational speed max values when the friction coefficient is increased. Speed transients in Figure 12 confirm that $\Omega_{max}$ decreases with $\mu$. However, $\mu$ must be chosen carefully since an abrupt

<figure><img src="image_13.png" /><figcaption>**Figure 12:** Sensitivity of the rotational speed transient to the friction coefficient: $\mu = 0.1$ (ref) (—), $\mu = 0.12$ (—), $\mu = 0.15$ (---), and $\mu = 0.2$ (----).</figcaption></figure>

decrease is obtained for $\mu \ge 0.15$. This phenomenon reflects that *dry whip* is taking place, i.e. fully
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 16:** Description of velocities and efforts during a dry whip instability because friction force acts in the direction of vibration speed $v_{prec}$ and opposed to rotational velocity $\Omega$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 17:** Stator deformation (●) and contact interface during dry whip phase at $t = 160$ ms: stator initial shape (—), contact points (x) and screws nodes (●)</figcaption></figure>

reflects a transition from dry whip to dry whirl, i.e. rolling motion without sliding. Indeed, the friction force is still non-zero but the tangential relative speed at contact points is nearly zero and friction no longer produces work. As a result, global frictional power (in Watts)

$$ \dot{Q}_f(t) = \sum_{j=1}^{N_e} \mu |v_{T,j}(t)| \lambda_j(t) \qquad (53) $$

reduces to zero with time. Note that during dry whirl, frictional torque shows high frequency oscillations between negative and positive values, see Figure 14(a). Literature [27] explains that these sign switches are caused by the Coulomb friction law chosen restricted to sliding, see Equation (31). Although the relative velocities at nodes $v_{T,j}$ for $j = 1, \dots, N_e$ tend to zero, they do not strictly equal zero and slightly change of sign, while the normal forces ($\lambda_j$) remain always positive.

To conclude, a higher friction coefficient implies a decrease of $\Omega_{max}$. For $\mu > 0.12$, dry whip is triggered, involving an exponential increase of efforts and braking torque. A transition from this unstable regime to dry whirl takes place once the relative velocity between the rotor and the stator goes to zero. A friction force still exists but it no longer produces work and structural vibrations are mitigated.

### 4.3. Sensitivity to the number of screws and magnitude of springs stiffness

The structural stiffness of the stator is known to affect the maximum rotational velocity. Depending on the engine type, the fixing conditions and stator stiffness may vary, and one must ensure that the stator design is in agreement with the speed limitation objective. Besides, it is important to identify if there exist stator configurations that lead to melting far from the screws. In the present work, the stiffness is modified in the following ways:

* the number of screws is reduced but stiffness of the connected springs is unchanged, or
annular backward precessional motion. The dynamical response is further investigated for $\mu = 0.15$ with mesh refinement. It turns out that $\Omega_{max}$ is identical for all mesh sizes, in Figure 13(a), and velocity profiles are similar until $t = 135$ ms. Again, perfect match is predicted for $t < 30$ ms for the rotor displacement $x_r(t)$ in Figure 13(b). Then, a phase lag takes place and dry whip is initiated in a time frame between 135 and 160 ms, which is more accurate than for a model with rigid stator and rotor. Fully annular contact is

<figure><img src="image_14.png" /><figcaption>**Figure 13:** Mesh convergence analysis for $\mu = 0.15$: $N_e = 36$ (—), $N_e = 48$ (—), $N_e = 60$ (—).</figcaption></figure>

confirmed by frictional torque that never vanishes after 150 ms. Violent dry whip also emerges with an

<figure><img src="image_15.png" /><figcaption>**Figure 14:** Dynamical response for $\mu = 0.15$ and $N_e = 36$ with emphasis of transition to dry whip.</figcaption></figure>

abrupt increase of the precessional velocity in absolute value as shown in Figure 14(c). In other words, kinetic energy is transferred from rotation to vibration through frictional mechanisms. Let us consider the front view of a rotorshaft with circular section depicted in Figure 15. In case of backward precessional motion Figure 15(b), friction force $F_f$ acts against $\Omega$ but in favor of $v_{prec}$. The consequence is that the rotor displacement (see Figure 13(b)) is amplified by a factor 10. Figure 16 also shows that the stator suffers significant deformation that could indicate a response out of the elastic domain for material properties. However, once $t > 160$ ms, frictional torque, presented in Figure 14(a) for $N_e = 36$, tends exponentially to zero and precessional speed stabilizes at -20 m/s. This reflects a transition from dry whip to dry whirl,
* the number of screws is kept at 12 but the springs stiffness $k_{Bu}$, $k_{Bv}$ and $k_A$ are all increased/decreased by 10%.

First, the number of screws, i.e. the number of springs, is set to 5 so that the whole stator becomes softer. It is important to keep in mind that these springs are a condensed model of plain parts that aim to simplify the boundary conditions applied to the stator. In order to work with the same number of elements as in the reference case, $N_e$ is set to 35. The friction coefficient is $\mu = 0.1$. Figure 18(a) shows that with fewer screws, the braking torque becomes less efficient since $\Omega_{max}$ is higher by 12%. Besides, the temperature

<figure><img src="image_22.png" /><figcaption>**Figure 18:** Sensitivity of thermomechanical response to screws number: reference $N_{sc} = 12$ (—) and $N_{sc} = 5$ (---).</figcaption></figure>

distributions in Figure 18(b) demonstrates that temperature peaks always correspond to the locations of the screws. Temperature in the vicinity of the screws presents again negative values related to the

<figure><img src="image_23.png" /><figcaption>**Figure 19:** Influence of screws number on contact response at simple node located on $\phi = \pi/9$: $N_{sc} = 12$ (ref) (—) and $N_{sc} = 5$ (---).</figcaption></figure>

lack of nodes. However, with same mesh refinement ($N_e = 70$ close to the 72 nodes for reference case), the temperature far from the screw nodes is lower when $N_{sc}$ decreases. Indeed, temperature in these regions ranges between 582 °C and 1014 °C, as shown in Figure 10(b). The average temperature value drops to 815 °C, which represents a reduction of 9.84 %. As already mentioned, this solution represents a temperature elevation if a non-zero initial condition is considered. For instance, if $T_0 \approx 300$ °C, fusion would not take place. This result may seem paradoxical since the rotational velocity, and so the relative sliding velocity ($v_T$), are higher by 12% when $N_{sc}$ drops from 12 to 5. That increase of sliding velocity is compensated by lower normal contact efforts. This is an expected outcome when the structure is softer and is confirmed by transient of normal forces. For instance, $F_N$ at the node with polar coordinate $\phi = \pi/9$ on the stator, see Figure 19(a), is smaller when $N_{sc}$ decreases during the whole simulation. One can also notice that the number of contact occurrences drops from 23 to 16 when the screws number decreases but they last longer. Contact is maintained during roughly 5 ms while this duration never exceeds 3 ms when $N_{sc} = 12$. As a result, slightly less heat is generated at nodes, as confirmed in Figure 19(b), but the lower contact efforts predicted with fewer screws lead to a less efficient braking mechanism.

In the second case study, transients of rotational velocity in Figure 20(a) indicate that when the stiffness at the screws is increased, then $\Omega_{max}$ decreases. This statement is in agreement with literature [14] who also found that contact efforts increased with stiffer stator supports. Indeed, Figure 20(b) confirms that
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15:** Description of velocities and efforts during a dry whip instability because friction force acts in the direction of vibration speed $v_{prec}$ and opposed to rotational velocity $\Omega$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 16:** Stator deformation (●-●) and contact interface during dry whip phase at $t = 160$ ms: stator initial shape (---), contact points (x) and screws nodes (●)</figcaption></figure>

i.e. rolling motion without sliding. Although friction force is non-zero (and proportional to normal force in sliding-only Coulomb friction model, cf Equation (32)), the tangential relative speed at contact points is nearly zero and thus friction produces almost zero work. As a result, global frictional power (in Watts)

$$ \dot{Q}_f(t) = \sum_{j=1}^{N_e} \mu |v_{T,j}(t)| \lambda_j(t) \qquad (54) $$

tends to zero with time. Note that during dry whirl, frictional torque shows high frequency oscillations between negative and positive values, see Figure 14(a). Literature [31] explains that these sign switches are caused by the Coulomb friction law chosen restricted to sliding. Even if the relative velocities at nodes $v_{T,j}$ for $j = 1, \dots, N_e$ tend to zero, they do not strictly equal zero and slightly change of sign, see Figure 14(d), while the normal forces ($\lambda_j$) remain always positive.

To conclude, a higher friction coefficient implies a decrease of $\Omega_{max}$. For $\mu > 0.12$, dry whip is triggered, involving an exponential increase of efforts and braking torque. A transition from this unstable regime to dry whirl takes place once the relative velocity between the rotor and the stator goes to zero. A friction force still exists but it no longer produces work and structural vibrations are mitigated.

### 4.3. Sensitivity to the number of screws and magnitude of springs stiffness

The structural stiffness of the stator is known to affect the maximum rotational velocity. Depending on the engine type, the fixing conditions and stator stiffness may vary, and one must ensure that the stator design is in agreement with the speed limitation objective. Besides, it is important to identify if there exist stator configurations that lead to melting far from the screws. In the present work, the stiffness is modified in the following ways:

* the number of screws is reduced but stiffness of the connected springs is unchanged, or
<figure><img src="image_24.png" /><figcaption>**Figure 20:** Sensitivity of thermomechanical response to springs stiffness with $N_e = 36$: reference (—), -10% (—) and +10% (---).</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 21:** Influence of supports stiffness on the temperature profile: reference (—), -10% (—), +10% (---)</figcaption></figure>

contact effort $F_N$ for a screw node (here node $j = 9$ corresponding to the 4th screw) depends proportionally on stiffnesses $k_{Bu}$, $k_{Bv}$ and $k_A$. This dependency is also reflected in the algorithmic of the solution method. For a given predicted penetration, a softer structure, or softer supports in this case, implies that the required effort to correct the gap is smaller. Stator stiffness also affects temperature profile. It seems that the relative sliding velocity decrease involved in the frictional energy is not dominant with regard to the increase of the contact effort. As a result, peaks at screws and average temperature increase with the stator support stiffness, see Figure 21. Indeed, with stiffness reduced (resp. increased) by 10%, the average temperature equals 828 °C (resp. 1032 °C). For the stiffer stator, temperature reaches 1226 °C locally. One should thus consider that melting wear takes place elsewhere that at the screws.

Therefore attention must be paid in the modeling of the stator supports. They should be calibrated with experiments since the rotordynamics response may not fulfill the objective $\Omega_{max} < \Omega_{lim}$ in a stiffness margin of 10%.

### 4.4. Thermal effect analysis

In the above investigation, the only coupling parameter between structural dynamics and thermodynamics is friction. The objective is now to include thermal expansion and investigate how it affects the speed transient. The thermal expansion coefficient for steel is set to $\alpha = 1.05 \times 10^{-5} \text{ K}^{-1}$. According to Figure 22(a), it turns out that $\Omega_{max}$ is slightly decreased by 2% and that the deceleration phase becomes steeper with thermal dilatation. Again, the torque generated by friction, shown in Figure 22(b), exceeds its counterpart from aerodynamics. However, it is not obvious that the frictional torque is higher when thermal dilation is considered although one can see in Figure 23 that the corresponding frictional energy is effectively greater, mechanism which should be understood. Figure 24 shows portions of the frictional energy provided individually to a screw node and a neighbor classical node. For any screw node, see for instance Figure 24(a) for node 19, the frictional energy is greater when thermal dilation is accounted for. For node 18, in Figure 24(b), the amount of thermal energy received decreases. It also seems that the energy reaches a plateau, which means there is no more heat source at the considered node. This comment holds for any classical node. Actually, the gap function at the classical node 18, indicated in Figure 25(b), shows that there is permanent contact separation after 220 ms. The thermal expansion is
* the number of screws is kept at 12 but the springs stiffness $k_{Bu}$, $k_{Bv}$ and $k_A$ are all increased/decreased by 10%.

First, the number of screws, i.e. the number of springs, is set to 5 so that the whole stator becomes softer. It is important to keep in mind that these springs are a condensed model of plain parts that aim to simplify the boundary conditions applied to the stator. In order to work with the same number of elements as in the reference case, $N_e$ is set to 35. The friction coefficient is $\mu = 0.1$. Figure 17(a) shows that with fewer screws, the braking torque becomes less efficient since $\Omega_{max}$ is higher by 12%. Besides,

<figure><img src="image_19.png" /><figcaption>**Figure 17:** Sensitivity of thermomechanical response to screws number: reference $N_{sc} = 12$ (—) and $N_{sc} = 5$ (---).</figcaption></figure>

the temperature distributions in Figure 17(b) demonstrate that temperature peaks always correspond to the screws. Temperature in the vicinity of the screws presents again negative values related to the

<figure><img src="image_20.png" /><figcaption>**Figure 18:** Influence of screws number on contact response at simple node located on $\phi = \pi/9$: $N_{sc} = 12$ (ref) (—) and $N_{sc} = 5$ (---).</figcaption></figure>

lack of nodes. However, with same mesh refinement ($N_e = 70$ close to the 72 nodes for reference case), the temperature far from the screw nodes is lower when $N_{sc}$ decreases. Indeed, temperature in these regions ranges between 582 °C and 1014 °C, as shown in Figure 17(b). The average temperature value drops to 815 °C, which represents a reduction of 9.84 %. As already mentioned, this solution represents a temperature elevation if a non-zero initial condition is considered. For instance, if $T_0 \approx 300$ °C, fusion would not take place. This result may seem paradoxical since the rotational velocity, and so the relative sliding velocity ($v_T$), are higher by 12% when $N_{sc}$ drops from 12 to 5. That increase of sliding velocity is compensated by lower normal contact efforts. This is an expected outcome when the structure is softer and is confirmed by transient of normal forces. For instance, $F_N$ at the node with polar coordinate $\phi = \pi/9$ on the stator, see Figure 18(a), is smaller when $N_{sc}$ decreases during the whole simulation. One can also notice that the number of contact occurrences drops from 23 to 16 when the screws number decreases but they last longer. Contact is maintained during roughly 5 ms while this duration never exceeds 3 ms when $N_{sc} = 12$. As a result, slightly less heat is generated at nodes, as confirmed in Figure 18(b), but the lower contact efforts predicted with fewer screws lead to a less efficient braking mechanism.

In the second case study, transients of rotational velocity in Figure 19(a) indicate that when the stiffness at the screws is increased, then $\Omega_{max}$ decreases. This statement is in agreement with literature [15] who also found that contact efforts increased with stiffer stator supports. Indeed, Figure 19(b) confirms that
<figure><img src="image_26.png" /><figcaption>**Figure 22:** Sensitivity of mechanical response: without thermal dilatation (—) and with thermal dilatation (---).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 23:** Frictional energy generated on the entire contact interface without dilatation (—) and with dilatation (---).</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Figure 24:** Energy source provided locally to stator nodes by friction, increases for any screw node with thermal dilatation, while it decreases for other nodes: $\alpha = 0$ (—) and $\alpha = 1.05 \times 10^{-5} \text{ K}^{-1}$ (---).</figcaption></figure>

already observable on the gap function. While the maximum rotor-stator distance equals about twice the initial gap in the reference case, this distance increases to more than 3 mm, which is close to $4g_0$.

The main conclusion of this analysis is that the stator expands due to the increased temperature but not homogeneously within the structure. This is illustrated in Figure 26 where the first quarter of the deformed ring at the end of the simulation is plotted. Because of the petal-shape of the stator, contact is concentrated at screws that were already carrying most of the contact efforts. These petals appear because the structure is stiffer at the screws. One could imagine that they are almost clamped. Heat affects more the rest of the ring that bends to the outer side. Higher efforts at the screws and smaller efforts on the other ones does not have much effect on the rotational velocity: the friction-induced braking is not altered by thermal expansion. However, the temperature magnitude predicted at the screw nodes is not realistic. More complex phenomena could take place such as wear and melting and thermal expansion could amplify metal fusion in the vicinity of the screws.
<figure><img src="image_21.png" /><figcaption>**Figure 19:** Sensitivity of thermomechanical response to springs stiffness with $N_e$ = 36: reference (—), -10% (—) and +10% (---).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 20:** Influence of supports stiffness on the temperature profile: reference (—), -10% (—), +10% (---)</figcaption></figure>

contact effort $F_N$ for a screw node (here node $j = 9$ corresponding to the 4th screw) depends proportionally on stiffness $k_{Bu}$, $k_{Bv}$ and $k_A$. This dependency is also reflected in the algorithmic of the solution method. For a given predicted penetration, a softer structure, or softer supports in this case, implies that the required effort to correct the gap is smaller. Stator stiffness also affects temperature profile. It seems that the relative sliding velocity decrease involved in the frictional energy is not dominant with regard to the increase of the contact effort. As a result, peaks at screws and average temperature increase with the stator support stiffness (see Figure 20). Indeed, with stiffness reduced (resp. increased) by 10%, the average temperature equals 828 °C (resp. 1032 °C). For the stiffer stator, temperature reaches 1226 °C locally. One should thus consider that melting wear takes place elsewhere that at the screws.

Therefore attention must be paid in the modeling of the stator supports. They should be calibrated with experiments since the rotordynamics response may not fulfill the objective $\Omega_{max} < \Omega_{lim}$ in a stiffness margin of 10%.

### 4.4. Thermal effect analysis

In the above investigation, the only coupling parameter between structural dynamics and thermodynamics is friction. The objective is now to include thermal expansion and investigate how it affects the speed transient. The thermal expansion coefficient for steel is set to $\alpha = 1.05 \times 10^{-5} \text{ K}^{-1}$. According to Figure 21(a), it turns out that $\Omega_{max}$ is slightly decreased by 2% and that the deceleration phase becomes steeper with thermal dilatation. Again, the torque generated by friction, shown in Figure 21(b), exceeds its counterpart from aerodynamics. However, it is not obvious that the frictional torque is higher when thermal dilation is considered although one can see in Figure 22 that the corresponding frictional energy is effectively greater, mechanism which should be understood. Figure 23 shows portions of the frictional energy provided individually to a screw node and a neighbor classical node. For any screw node, see for instance Figure 23(a) for node 19, the frictional energy is greater when thermal dilation is accounted for. For node 18, in Figure 23(b), the amount of thermal energy received decreases. It also seems that the energy reaches a plateau, which means there is no more heat source at the considered node. This comment holds for any classical node. Actually, the gap function at the classical node 18, indicated in Figure 24(b), shows that there is permanent contact separation after 220 ms. The thermal expansion is
<figure><img src="image_29.png" /><figcaption>**Figure 25:** Gap function at node 17: $\alpha = 0$ (—) and $\alpha = 1.05 \times 10^{-5} \text{ K}^{-1}$ (—)</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 26:** Rotor-Stator contact interface with thermal dilatation at the end of simulation $t = 300$ ms for a 90°-sector: Rotor rigid disc (—), stator deformed shape (●), contact points (x) and screws (●). Amplitudes are magnified by a scale factor 5 for illustration purposes. Blue points represent screw nodes and black ones, conventional nodes. A node marked with a cross indicates closed contact.</figcaption></figure>

**5. Conclusions**

The present work proposed to investigate how unilateral and frictional contact occurrences combined with thermomechanical coupling affect the unknown rotational velocity transient of a free turbine. A constant driving torque is initially applied to the turbine which as a consequence tends to accelerate. However, it is expected that frictional mechanisms overcome that torque in such a way that the turbine eventually slows down. Such a response is expected to involve strongly non-linear dynamics because the unknown rotational velocity that governs gyroscopic and inertial effects, and must be addressed properly. The rotordynamics is nonlinear even without structural interaction and required an adaptation of Carpenter's algorithm. Moreover, unilateral contact with friction couples the structural dynamics to the thermodynamics though the generation of heat sources at the contact points. Linear thermodynamics and dilation were developed in the governing equations. Both physics are simultaneously solved with strong coupling within the same time-stepping solution method.

The investigation of the transient response is conducted on a rotor-stator system where the rotor was assumed to be perfectly rigid. The stator assembly, made of three parts in practice, was simplified into a flexible circular ring, connected to a point mass, and supported by springs at the screws level. The ring was meshed with curved Bernoulli beams elements. Results showed that for a given reference configuration, the rotational velocity remains bounded, reflecting the ability of the frictional torque to overcome the prescribed aerodynamic torque and brake the turbine. The origin of such high contact efforts is that the dynamic response is governed by backward precessional motion. Boundary conditions applied on the stator with stiff springs placed on a finite number of nodes, reflecting the presence of screws, caused a major disparity of contact effort mainly concentrated on these nodes. For the reference case, contact remained partial but the vibration of the structures is amplified which are the characteristics of a response in backward precessional motions. The rotorshaft entered into unstable dry whip operating conditions when the friction coefficient was increased. The thermal solution showed that a huge frictional power is generated, reaching hundreds of kiloWatts at screw nodes, because of the high rotational velocity and high contact efforts. As a result, temperature peaks beyond the melting point of the considered material
<figure><img src="image_23.png" /><figcaption>**Figure 21:** Sensitivity of mechanical response: without thermal dilatation (—) and with thermal dilatation (---).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 22:** Frictional energy generated on the entire contact interface without dilatation (—) and with dilatation (---).</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 23:** Energy source provided locally to stator nodes by friction, increases for any screw node with thermal dilatation, while it decreases for other nodes: $\alpha = 0$ (—) and $\alpha = 1.05 \times 10^{-5} \text{ K}^{-1}$ (---).</figcaption></figure>

already observable on the gap function. While the maximum rotor-stator distance equals about twice the initial gap in the reference case, this distance increases to more than 3 mm, which is close to $4g_0$.

The main conclusion of this analysis is that the stator expands due to the increased temperature but not homogeneously within the structure. This is illustrated in Figure 25 where the first quarter of the deformed ring at the end of the simulation is plotted. Because of the petal-shape of the stator, contact is concentrated at screws that were already carrying most of the contact efforts. These petals appear because the structure is stiffer at the screws. One could imagine that they are almost clamped. Heat affects more the rest of the ring that bends to the outer side. Higher efforts at the screws and smaller efforts on the other ones does not have much effect on the rotational velocity: the friction-induced braking is not altered by thermal expansion. However, the temperature magnitude predicted at the screw nodes is not realistic. More complex phenomena could take place such as wear and melting and thermal expansion could amplify metal fusion in the vicinity of the screws.
are found on the screw nodes. This extremely localized heat distribution reflects that the discretization mesh must be sufficiently refined to avoid non-physical results. Sensitivity studies provided responses in agreement with the existing literature and demonstrated that the stator spring supports must be designed carefully. Many screw points with stiffer springs certainly make the braking property of friction more efficient, but this comes with more heat generated at nodes that are more likely to wear by melting. Eventually, thermal expansion, when accounted for, does not significantly affect the speed transient, but the stator presented a non-uniform outward deformation. The obtained petal-shaped topology forces contact events to occur at the screw nodes only. The predicted screws temperature is rather qualitative since the implemented thermodynamics model is limited in its scope.

## Acknowledgements

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

## Appendix A. Elementary matrices for the stator

The matrices below represent elementary mass, stiffness and thermal dilatation matrices for a curvilinear Bernoulli beam element:

$$ M_e = \frac{\rho S l_e}{420} \begin{pmatrix} 156 & 22l_e & 0 & 0 & 54 & -13 & 0 & 0 \\ 4l_e^2 & 0 & 0 & 13l_e & -3l_e^2 & 0 & 0 & 0 \\ & 156 & 22l_e & 0 & 0 & 54 & -13l_e & 0 \\ & 4l_e^2 & 0 & 0 & 13l_e & -3l_e^2 & 0 & 0 \\ & & 156 & -22l_e & 0 & 0 & 0 & 0 \\ & & & 4l_e^2 & 0 & 0 & 0 & 0 \\ & & & & 156 & -22l_e & 0 & 0 \\ & & & & & 4l_e^2 & 0 & 0 \end{pmatrix} \quad (A.1) $$

$$ K_e = \rho S l_e \begin{pmatrix} \frac{6S}{5l_e} & \frac{S}{10} & -\frac{S}{2R_s} & -\frac{S l_e}{10 R_s} & \frac{6S}{5l_e} & \frac{S}{10} & -\frac{S}{2R_s} & \frac{S l_e}{10 R_s} \\ \frac{2S l_e}{15} & \frac{S l_e}{10 R_s} & 0 & -\frac{S l_e}{10} & -\frac{S l_e}{30} & -\frac{S l_e}{10 R_s} & C & F \\ & A & B & \frac{S}{2R_s} & -\frac{S l_e}{60 R_s} & -\frac{S l_e}{60 R_s} & -F & G \\ & & D & \frac{S l_e}{10 R_s} & -\frac{S l_e}{10} & \frac{S}{2R_s} & \frac{S l_e}{10 R_s} & 0 \\ & & & \frac{6S}{5l_e} & \frac{S}{10} & \frac{2S l_e}{15} & -\frac{S l_e}{10 R_s} & -B \\ & & & & & & D \end{pmatrix}, \quad K^{uT,e} = \frac{\alpha E S}{60} \begin{pmatrix} 30 & 30 \\ -5l_e & 5l_e \\ -\frac{21l_e}{R_s} & -\frac{9l_e}{R_s} \\ -\frac{3l_e^2}{R_s} & -\frac{2l_e^2}{R_s} \\ -30 & -30 \\ 5l_e & -5l_e \\ -\frac{9l_e}{R_s} & -\frac{21l_e}{R_s} \\ \frac{2l_e^2}{R_s} & \frac{3l_e^2}{R_s} \end{pmatrix} $$

with

$$ A = \frac{13l_e^4 SR_s^2 + 13l_e^4 I_z - 84I_z l_e^2 R_s^2 + 420I_z R_s^4}{35l_e^3 R_s^4} $$

$$ C = 3 \frac{3l_e^4 SR_s^2 + 3l_e^4 I_z + 56I_z l_e^2 R_s^2 - 280I_z R_s^4}{70l_e^3 R_s^4} $$

$$ F = \frac{2520I_z R_s^4 - 13l_e^4 SR_s^2 - 13l_e^4 I_z - 84I_z l_e^2 R_s^2}{420l_e^2 R_s^4} $$

$$ B = \frac{11l_e^4 SR_s^2 + 11l_e^4 I_z - 252I_z l_e^2 R_s^2 + 1260I_z R_s^4}{210l_e^2 R_s^4} $$

$$ D = \frac{l_e^4 SR_s^2 + l_e^4 I_z - 28I_z l_e^2 R_s^2 + 420I_z R_s^4}{105l_e R_s^4} $$

$$ G = \frac{840I_z R_s^4 - 3l_e^4 SR_s^2 - 3l_e^4 I_z + 28I_z l_e^2 R_s^2}{420l_e^2 R_s^4} $$

## Appendix B. Formulation of normal and tangential contribution matrices of Lagrange multipliers

The matrices $C_N$ and $C_T$ define respectively the Jacobians of the normal, and tangential distances associated to the active nodes $j \in A$ of the stator. The index of node $j$ in the set $A$ is defined by variable $k$. In other words, each column vector $C_N^{(k)}$ (resp. $C_T^{(k)}$) is computed as $\nabla_q d_N^j(q)$ (resp. $\nabla_q d_T^j(q)$) and is empty except for degrees of freedom of the rotor and of node $j$ belonging to the stator. Portion of contributions vectors for the rotor are stressed by $r$ index and come directly from Equations (37) and (40):

$$ \left\{ \begin{aligned} \mathbf{C}_{\mathrm{N},r}^{(k)}(\mathbf{q}) &= \frac{1}{\delta_{\mathrm{N}}^{j}(\mathbf{q})} (0, 0, -(x_{sj}(\mathbf{q}) - x_r), -(y_{sj}(\mathbf{q}) - y_r), 0)^{\mathrm{T}} \\ \mathbf{C}_{\mathrm{T},r}^{(k)}(\mathbf{q}) &= \frac{1}{\delta_{\mathrm{N}}^{j}(\mathbf{q})} (0, 0, (y_{sj}(\mathbf{q}) - y_r), -(x_{sj}(\mathbf{q}) - x_r), -R_r)^{\mathrm{T}} \end{aligned} \right. \quad (B.1) $$
<figure><img src="image_26.png" /><figcaption>**Figure 24:** Gap function at node 18: $\alpha = 0$ (—) and $\alpha = 1.05 \times 10^{-5} \text{ K}^{-1}$ (—)</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 25:** Rotor-Stator contact interface with thermal dilatation at the end of simulation $t = 300$ ms for a 90°-sector: Rotor rigid disc (—), stator deformed shape (●), contact points (x) and screws (●). Amplitudes are magnified by a scale factor 5 for illustration purposes. Blue points represent screw nodes and black ones, conventional nodes. A node marked with a cross indicates closed contact.</figcaption></figure>

**5. Conclusions**

The present work proposed to study how unilateral and frictional contact occurrences combined with thermomechanical coupling affect the unknown rotational velocity transient of a free turbine. This standalone configuration results from shaft failure, which is a topic still to explore and involving complex phenomena from aerodynamics and rotordynamics. The present research focused on a simplified scenario. The investigation of the transient response was conducted on a rotor-stator system where the rotor was assumed to be perfectly rigid. The stator assembly, made of three parts in practice, was simplified into a flexible circular ring, connected to a point mass, and supported by springs at the screws level. The ring was meshed with curved Bernoulli beams elements. A constant driving torque was initially applied to the turbine which, consequently, tends to accelerate. However, it was expected that frictional mechanisms overcome that torque in such a way that the turbine eventually slows down. Such a response was expected to involve strongly non-linear dynamics because the unknown rotational velocity that governs gyroscopic and inertial effects, and must be addressed properly. The rotordynamics is nonlinear even without structural interaction and required an adaptation of Carpenter's algorithm. Moreover, unilateral contact with friction couples the structural dynamics to the thermodynamics though the generation of heat sources at the contact points. Linear thermodynamics and dilatation were developed in the governing equations. Both physics were simultaneously solved with strong coupling within the same time-stepping solution method.

Results showed that for a given reference configuration, the rotational velocity remained bounded, reflecting the ability of the frictional torque to overcome the prescribed aerodynamic torque and brake the turbine. The origin of such high contact efforts is that the dynamic response is governed by backward precessional motion. Boundary conditions applied on the stator with stiff springs placed on a finite number of nodes, reflecting the presence of screws, caused a major disparity of contact effort mainly concentrated on these nodes. For the reference case, contact remained partial but the vibration of the structures is amplified which are the characteristics of a response in backward precessional motions. The rotorshaft entered into unstable dry whip operating conditions when the friction coefficient was increased. The thermal solution showed that a huge frictional power is generated, reaching hundreds of kiloWatts at screw
The analytic expressions of $C_N^{(k)}$ and $C_T^{(k)}$ portions for stator (stressed by extra index s) require to apply a composed derivation on the jth node cartesian coordinates ($x_{sj}$, $y_{sj}$). Injecting Equation (17) into Equations (37) and (40) and proceeding to gradient operation thus leads to the following matricial form

$$ \left\{ \begin{aligned} C_{N,s}^{(k)}(\mathbf{q}) &= \frac{1}{\delta_N^j(\mathbf{q})} \left( (x_{sj}(\mathbf{q}) - x_r)\mathbf{P}_x + (y_{sj}(\mathbf{q}) - y_r)\mathbf{P}_y \right) \mathbf{U} \\ C_{T,s}^{(k)}(\mathbf{q}) &= \frac{1}{\delta_T^j(\mathbf{q})} \left( -(y_{sj}(\mathbf{q}) - y_r)\mathbf{P}_x + (x_{sj}(\mathbf{q}) - x_r)\mathbf{P}_y \right) \mathbf{U} \end{aligned} \right. \quad (B.2) $$

## Appendix C. Adaptation of Carpenter’s algorithm to thermomechanical applications

Algorithm 1 is the modified Carpenter algorithm accounting for nonlinear rotordynamics because of inertial and gyroscopic effects, multi-constraint management and thermomechanical coupling.

### Algorithm 1 Carp-oneVC

Require: $\mathbf{q}_0, \mathbf{q}_1$

for $i = 1, \dots, N-1$ do

{Initial guess for contact-free solution}

$\mathbf{q}_{i+1}^{p,0} \leftarrow \mathbf{q}_i$

{Nonlinear solver: finds solution $\mathbf{q}_{i+1}^p$ of Equation (33) with initial guess $\mathbf{q}_{i+1}^0$}

Solve $\bar{\mathbf{M}}(\mathbf{q}_{i+1}^p)\mathbf{q}_{i+1}^p + \bar{\mathbf{D}}(\mathbf{q}_{i+1}^p)\mathbf{q}_i + \bar{\mathbf{K}}(\mathbf{q}_{i+1}^p)\mathbf{q}_{i-1} - \bar{\mathbf{f}}_i = \mathbf{0}$

{Gap computation through Equation (37)}

for $j = 1, \dots, N_e$ do

$$ d_N^j \leftarrow \delta_N^j(\mathbf{q}_{i+1}^p) - R_r - a $$

end for

{Search of active constraints}

$$ \mathcal{A} \leftarrow \{j | d_N^j < 0\} $$

if $\mathcal{A} = \emptyset$ then

{No correction}

$$ \mathbf{\Lambda}_{i+1} \leftarrow \mathbf{0} $$

$$ \mathbf{q}_{i+1} \leftarrow \mathbf{q}_{i+1}^p $$

else

{Correction}

$$ \mathbf{v}_i^p \leftarrow (\mathbf{q}_{i+1}^p - \mathbf{q}_i)/h $$

$$ \mathbf{g}_{i+1}^p \leftarrow [d_N^j | j \in \mathcal{A}] $$

for $k = 1, \dots, |\mathcal{A}|$ do

$$ C_N^{(k)} \leftarrow \nabla_q d_N^j $$

$$ C_T^{(k)} \leftarrow \nabla_q d_T^j $$

$$ \mathbf{v}_T^{(k)} \leftarrow C_T^{(k)} \mathbf{v}_i^p $$

$$ C_{\theta}^{(k)} \leftarrow \frac{\theta_s}{\theta_s + \theta_r} \mu |\mathbf{v}_T^{(k)}| $$

end for

$$ C_{NT} \leftarrow \cup_{k=1}^{|\mathcal{A}|} C_N^{(k)} - \mu \mathrm{sign}(\mathbf{v}_T^{(k)}) C_T^{(k)} $$

{Block merge of mechanical and heat contributions of multipliers}

$$ C_{NT\theta} \leftarrow (C_{NT}C_{\theta})^{\top} $$

{Correction}

$$ \mathbf{\Lambda}_{i+1} \leftarrow -(C_N^{\mathrm{T}} \bar{\mathbf{M}}(\mathbf{q}_{i+1}^p)^{-1} C_{NT\theta})^{-1} \mathbf{g}_{N,i+1}^p $$

$$ \mathbf{q}_{i+1} \leftarrow \mathbf{q}_{i+1}^p + \bar{\mathbf{M}}(\mathbf{q}_{i+1}^p)^{-1} C_{NT\theta} \mathbf{\Lambda}_{i+1} $$

end if

end for

## References

[1] V. Acary and B. Brogliato. *Numerical methods for nonsmooth dynamical systems*, volume 35 of *Lecture Notes in Applied and Computational Mechanics*. Springer, 2008. ISBN 978-3-540-75391-9. [doi]. [oai].

[2] V. Acary, O. Bonnefon, M. Brémond, O. Huber, F. Pérignon, and al. An introduction to Siconos. Technical Report RT-0340, INRIA, 2019. [oai].

[3] S. Ahmad. Rotor casing contact phenomenon in rotor dynamics — literature survey. *Journal of Vibration and Control*, 16(9):1369–1377, 2010. ISSN 1077-5463. [doi].

[4] P. Almeida, C. Gibert, F. Thouverez, and J.-P. Ousty. On some physical phenomena involved in blade-casing contact. In 9th International Conference on *Structural Dynamics*, pages 1–9, 2014. [oai].
nodes, because of the high rotational velocity and high contact efforts. As a result, temperature peaks beyond the melting point of the considered material were found on the screw nodes. This extremely localized heat distribution reflected that the discretization mesh must be sufficiently refined to avoid non-physical results. Sensitivity studies provided responses in agreement with the existing literature and demonstrated that the stator spring supports must be designed carefully. Many screw points with stiffer springs certainly make the braking property of friction more efficient, but this comes with more heat generated at nodes that are more likely to wear by melting. Eventually, thermal expansion, when accounted for, did not significantly affect the speed transient, but the stator presented a non-uniform outward deformation. The obtained petal-shaped topology forced contact events to occur at the screw nodes only. The predicted screws temperature were rather qualitative since the implemented thermodynamics model was limited in its scope.

## Acknowledgments

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

## Appendix A. Elementary matrices for the stator

The matrices below represent elementary mass, stiffness and thermal dilatation matrices for a curvilinear Bernoulli beam element:

$$ M_e = \frac{\rho S l_e}{420} \begin{pmatrix} 156 & 22l_e & 0 & 0 & 54 & -13 & 0 & 0 \\ 4l_e^2 & 0 & 0 & 13l_e & -3l_e^2 & 0 & 0 & 0 \\ & 156 & 22l_e & 0 & 0 & 54 & -13l_e & 0 \\ & & 4l_e^2 & 0 & 0 & 13l_e & -3l_e^2 & 0 \\ & & & 156 & -22l_e & 0 & 0 & 0 \\ & & & & 4l_e^2 & 0 & 0 & 0 \\ & & & & & 156 & -22l_e & 0 \\ & & & & & & 4l_e^2 & 0 \end{pmatrix} \quad (A.1) $$

$$ K_e = \rho S l_e \begin{pmatrix} \frac{6S}{5l_e} & \frac{S}{10} & -\frac{S}{2R_s} & -\frac{S l_e}{10 R_s} & \frac{6S}{5l_e} & \frac{S}{10} & -\frac{S}{2R_s} & \frac{S l_e}{10 R_s} \\ \frac{2S l_e}{15} & \frac{S l_e}{10 R_s} & 0 & -\frac{S}{10} & -\frac{S l_e}{30} & -\frac{S l_e}{10 R_s} & C & F \\ A & B & \frac{S}{2R_s} & -\frac{S l_e}{10 R_s} & \frac{S l_e}{60 R_s} & -F & G & D \\ D & \frac{S l_e}{10 R_s} & \frac{6S}{5l_e} & -\frac{S}{10} & \frac{S l_e}{2S l_e}{15} & \frac{S}{2R_s} & \frac{S l_e}{10 R_s} & 0 \\ & & & & & & A & -B \\ & & & & & & & D \end{pmatrix}, \quad K^{uT,e} = \frac{\alpha E S}{60} \begin{pmatrix} 30 & 30 \\ -5l_e & 5l_e \\ -\frac{21l_e}{R_s} & -\frac{9l_e}{R_s} \\ -\frac{3l_e^2}{R_s} & -\frac{2l_e^2}{R_s} \\ -30 & -30 \\ 5l_e & -5l_e \\ -\frac{9l_e}{R_s} & -\frac{21l_e}{R_s} \\ \frac{2l_e^2}{R_s} & \frac{3l_e^2}{R_s} \end{pmatrix} $$

with

$$ A = \frac{13l_e^4 SR_s^2 + 13l_e^4 I_z - 84I_z l_e^2 R_s^2 + 420I_z R_s^4}{35l_e^3 R_s^4} $$

$$ C = 3 \frac{3l_e^4 SR_s^2 + 3l_e^4 I_z + 56I_z l_e^2 R_s^2 - 280I_z R_s^4}{70l_e^3 R_s^4} $$

$$ F = \frac{2520I_z R_s^4 - 13l_e^4 SR_s^2 - 13l_e^4 I_z - 84I_z l_e^2 R_s^2}{420l_e^2 R_s^4} $$

$$ B = \frac{11l_e^4 SR_s^2 + 11l_e^4 I_z - 252I_z l_e^2 R_s^2 + 1260I_z R_s^4}{210l_e^2 R_s^4} $$

$$ D = \frac{l_e^4 SR_s^2 + l_e^4 I_z - 28I_z l_e^2 R_s^2 + 420I_z R_s^4}{105l_e R_s^4} $$

$$ G = \frac{840I_z R_s^4 - 3l_e^4 SR_s^2 - 3l_e^4 I_z + 28I_z l_e^2 R_s^2}{420l_e^2 R_s^4} $$

## Appendix B. Formulation of normal and tangential contribution matrices of Lagrange multipliers

The matrices $C_N$ and $C_T$ define respectively the Jacobians of the normal, and tangential distances associated to the active nodes $j \in A$ of the stator. The index of node $j$ in the set $A$ is defined by variable $k$. In other words, each column vector $C_N^{(k)}$ (resp. $C_N^{(k)}$) is computed as $\nabla_q d_N^j(q)$ (resp. $\nabla_q d_T^j(q)$) and is empty except for degrees of freedom of the rotor and of node $j$ belonging to the stator. Portion of contributions vectors for the rotor are stressed by $r$ index and come directly from Equations (38) and (41):

$$ \left\{ \begin{array}{l} C_{N,r}^{(k)}(\mathbf{q}) = \frac{1}{\delta_N'(\mathbf{q})} (0, 0, -(x_{sj}(\mathbf{q}) - x_r), -(y_{sj}(\mathbf{q}) - y_r), 0)^\top \\ C_{T,r}^{(k)}(\mathbf{q}) = \frac{1}{\delta_N'(\mathbf{q})} (0, 0, (y_{sj}(\mathbf{q}) - y_r), -(x_{sj}(\mathbf{q}) - x_r), -R_r)^\top \end{array} \right. \quad (B.1) $$
[5] A. Bartha. *Dry friction backward whirl of rotors*. PhD thesis, ETH Zurich, 2000. [doi].

[6] B. Berthoul, A. Batailly, L. Stainier, M. Legrand, and P. Cartraud. Phenomenological modeling of abradable wear in turbomachines. *Mechanical Systems and Signal Processing*, 98:770–785, 2018. [doi], [oai].

[7] H. Black. Interaction of a whirling rotor with a vibrating stator across a clearance annulus. *Journal of Mechanical Engineering Science*, 10(1):1–12, 1968. ISSN 0022-2542, 2058-3389. [doi].

[8] R. Blevins. Formulas for natural frequency and mode shape. Technical report, Van Nostrand Reynhold Company, 1979.

[9] N. Carpenter, R. Taylor, and M. Katona. Lagrange constraints for transient finite element surface contact. *International Journal for Numerical Methods in Engineering*, 32(1):103–128, 1991. [doi], [oai].

[10] H.-S. Chen and M. Stadtherr. A modification of Powell's dogleg method for solving systems of nonlinear equations. *Computers & Chemical Engineering*, 5(3):143–150, 1981. [doi].

[11] S. Chen and M. Géradin. Finite element simulation of non-linear transient response due to rotor-stator contact. *Engineering Computations*, 14(6):591–603, 1997. [doi].

[12] D. Childs and A. Bhattacharya. Prediction of dry-friction whirl and whip between a rotor and a stator. *Journal of Vibration and Acoustics*, 129(3):355–362, 2007. ISSN 1048-9002. [doi].

[13] Y. Choi. Investigation on the whirling motion of full annular rotor rub. *Journal of Sound and Vibration*, 258 (1):191–198, 2002. ISSN 0022-460X. [doi].

[14] F. Choy and J. Padovan. Non-linear transient analysis of rotor-casing rub events. *Journal of Sound and Vibration*, 113(3):529–545, 1987. ISSN 0022-460X. [doi].

[15] X. Dai, Z. Jin, and X. Zhang. Dynamic Behavior of the full rotorstop rubbing: numerical simulation and experimental verification. *Journal of Sound and Vibration*, 251(5):807–822, 2002. [doi].

[16] C. Delebarre, V. Wagner, J. Paris, G. Dessein, J. Denape, and J. Gurt-Santanach. The wear mechanisms occurring in a labyrinth seal/abradable contact depending on the incursion depth parameter. *Mechanics & Industry*, 17(6):601, 2016. [doi], [oai].

[17] A. Draganis. Finite element modeling of transient thermomechanical rolling contact featuring mixed control of the rigid body motion. *Journal of Tribology*, 139(1):011503, 2017. ISSN 0742-4787, 1528-8897. [doi].

[18] S. Dubigeon. *Mécanique des Milieux Continus*. Lavoisier, Tec & Doc, Nantes, 1998.

[19] C. Dumartineix. *Modélisation et étude de la dynamique complexe d'un système bi-rotor aubagé couplé*. PhD thesis, Ecole Centrale de Lyon, 2019. [oai].

[20] C. Duran, L. Manin, M.-A. Andrianoely, C. Bordegaray, F. Battle, and R. Dufour. Effect of rotor-stator contact on the mass unbalance response. In *9th IFToMM International Conference on Rotor Dynamics*, pages 1965–1975, 2015. [doi], [oai].

[21] F. Fatarella. *On the Dynamics of Reverse Whirl Due to Rotor/Stator Interaction*. PhD thesis, Imperial College of Science, Medicine & Technology, 1999.

[22] E. Fried. *Thermal Conduction Contribution to Heat Transfer at Contacts*. Thermal Conductivity. Academic Press, 2 edition, 1969.

[23] G. Genta. *Dynamics of Rotating Systems*. Springer, 2005. ISBN 978-0-387-20936-4. [doi].

[24] M. Géradin and D. Rixen. *Mechanical vibrations: theory and application to structural dynamics*. Wiley, 3rd edition, 2014. ISBN 978-1-118-90020-8.

[25] P. Ireman, A. Klarbring, and N. Strömberg. Finite element algorithms for thermoelastic wear problems. *European Journal of Mechanics - A/Solids*, 21(3):423–440, 2002. [doi], [oai].

[26] C. Jacobs, M. Legrand, F. Thouverez, and P. Almeida. Turbomachinery transient dynamics of radial rotor-stator contact occurrences with friction. In *ASME Turbo Expo*, pages V11BT26A007-1-12, 2023. [doi], [oai].

[27] G. Jacquet-Richardet, M. Torkhani, P. Cartraud, F. Thouverez, T. Nouri Baranger, M. Herran, C. Gibert, S. Baguet, P. Almeida, and L. Peletan. Rotor to stator contacts in turbomachines. Review and application. *Mechanical Systems and Signal Processing*, 40(2):401–420, 2013. ISSN 0888-3270. [doi], [oai].

[28] M. Jean. The non-smooth contact dynamics method. *Computer Methods in Applied Mechanics and Engineering*, 177(3-4):235–257, 1999. [doi], [oai].

[29] J. Jiang. Determination of the global responses characteristics of a piecewise smooth dynamical system with contact. *Nonlinear Dynamics*, 57(3):351–361, 2009. ISSN 0924-090X, 1573-269X. [doi].

[30] M. Lalanne and G. Ferraris. *Rotordynamics Prediction in Engineering*. Wiley, 2 edition, 1998. ISBN 0-585-32863-3 978-0-585-32863-8.

[31] M. Legrand. *Modèles de prédiction de l'interaction rotor/stator dans un moteur d'avion*. PhD thesis, Ecole Centrale de Nantes, 2006. [oai].

[32] M. Legrand, C. Pierre, P. Cartraud, and J.-P. Lombard. Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction. *Journal of Sound and Vibration*, 319(1-2):366–391, 2009. [doi], [oai].

[33] N. Lesaffre, J.-J. Sinou, and F. Thouverez. Contact analysis of a flexible bladed-rotor. *European Journal of Mechanics - A/Solids*, 26:541–557, 2007. [doi], [oai].

[34] H. Ma, Q. Zhao, X. Zhao, Q. Han, and B. Wen. Dynamic characteristics analysis of a rotor-stator system under different rubbing forms. *Applied Mathematical Modelling*, 39(8):2392–2408, 2015. [doi].
The analytic expressions of $C_N^{(k)}$ and $C_T^{(k)}$ portions for stator (stressed by extra index s) require to apply a composed derivation on the jth node cartesian coordinates ($x_{sj}, y_{sj}$). Injecting Equation (18) into Equations (38) and (41) and proceeding to gradient operation thus leads to the following matricial form

$$ egin{cases} C_{N,s}^{(k)}(\mathbf{q}) = \frac{1}{\delta_N^j(\mathbf{q})} ((x_{sj}(\mathbf{q}) - x_r)\mathbf{P}_x + (y_{sj}(\mathbf{q}) - y_r)\mathbf{P}_y) \mathbf{U} \\ C_{T,s}^{(k)}(\mathbf{q}) = \frac{1}{\delta_T^j(\mathbf{q})} ((-y_{sj}(\mathbf{q}) - y_r)\mathbf{P}_x + (x_{sj}(\mathbf{q}) - x_r)\mathbf{P}_y) \mathbf{U} \end{cases} \quad (B.2) $$

## Appendix C. Adaptation of Carpenter’s algorithm to thermomechanical applications

Algorithm 1 is the modified Carpenter algorithm accounting for nonlinear rotordynamics because of inertial and gyroscopic effects, multi-constraint management and thermomechanical coupling.

### Algorithm 1 Carp-oneVC

**Require:** $\mathbf{q}_0, \mathbf{q}_1$

for $i = 1, \dots, N-1$ do

{Initial guess for contact-free solution}

$$ \mathbf{q}_{i+1}^{p,0} \leftarrow \mathbf{q}_i $$

{Nonlinear solver: finds solution $\mathbf{q}_{i+1}^p$ of Equation (34) with initial guess $\mathbf{q}_{i+1}^0$}

Solve $\bar{\mathbf{M}}(\mathbf{q}_{i+1}^p)\mathbf{q}_{i+1}^p + \bar{\mathbf{D}}(\mathbf{q}_{i+1}^p)\mathbf{q}_i + \bar{\mathbf{K}}(\mathbf{q}_{i+1}^p)\mathbf{q}_{i-1} - \bar{\mathbf{f}}_i = \mathbf{0}$

{Gap computation through Equation (38)}

for $j = 1, \dots, N_e$ do

$$ d_N^j \leftarrow \delta_N^j(\mathbf{q}_{i+1}^p) - R_r - a $$

end for

{Search of active constraints}

$$ \mathcal{A} \leftarrow \{j | d_N^j < 0\} $$

if $\mathcal{A} = \emptyset$ then

{No correction}

$$ \mathbf{\Lambda}_{i+1} \leftarrow \mathbf{0} $$

$$ \mathbf{q}_{i+1}^p \leftarrow \mathbf{q}_{i+1}^p $$

else

$$ \mathbf{v}_i^p \leftarrow (\mathbf{q}_{i+1}^p - \mathbf{q}_i)/h $$

$$ \mathbf{g}_{i+1}^p \leftarrow [d_N^j | j \in \mathcal{A}] $$

for $k = 1, \dots, |\mathcal{A}|$ do

$$ C_N^{(k)} \leftarrow \nabla_q d_N^j $$

$$ C_T^{(k)} \leftarrow \nabla_q d_T^j $$

$$ \mathbf{v}_T^{(k)} \leftarrow C_T^{(k)} \mathbf{v}_i^p $$

$$ C_{\theta}^{(k)} \leftarrow \frac{\partial_s}{\partial_s + \partial_r} \mu |\mathbf{v}_T^{(k)}| $$

end for

$$ C_{NT} \leftarrow \cup_{k=1}^{|\mathcal{A}|} C_N^{(k)} - \mu \mathrm{sign}(\mathbf{v}_T^{(k)}) C_T^{(k)} $$

{Block merge of mechanical and heat contributions of multipliers}

$$ C_{NT\theta} \leftarrow (C_{NT}C_{\theta})^{\top} $$

{Correction}

$$ \mathbf{\Lambda}_{i+1} \leftarrow -(C_N^{\top} \bar{\mathbf{M}}(\mathbf{q}_{i+1}^p)^{-1} C_{NT\theta})^{-1} \mathbf{g}_{N,i+1}^p $$

$$ \mathbf{q}_{i+1} \leftarrow \mathbf{q}_{i+1}^p + \bar{\mathbf{M}}(\mathbf{q}_{i+1}^p)^{-1} C_{NT\theta} \mathbf{\Lambda}_{i+1} $$

end if

end for

## References

[1] Acary, V., Bonnefon, O., Brémond, M., Huber, O., Pérignon, F., Sinclair, S., 2019. An introduction to Siconos. Technical Report RT-0340. INRIA. oai:inria-00162911.

[2] Acary, V., Brogliato, B., 2008. Numerical methods for nonsmooth dynamical systems. Number 35 in Lecture Notes in Applied and Computational Mechanics, Springer. doi:10.1007/978-3-540-75392-6. oai:hal-00423530.

[3] Ahmad, S., 2010. Rotor casing contact phenomenon in rotor dynamics — literature survey. Journal of Vibration and Control 16, 1369–1377. doi:10.1177/1077546309341605.

[4] Almeida, P., Gibert, C., Thouverez, F., Ousty, J.P., 2014. On some physical phenomena involved in blade-casing contact, in: 9th International Conference on Structural Dynamics, pp. 1–9. oai:hal-02121521.
[35] J.-J. Moreau. Unilateral contact and dry friction in finite freedom dynamics. In J.-J. Moreau and P. Panagiotopoulos, editors, *Nonsmooth Mechanics and Applications*, volume 302 of *International Centre for Mechanical Sciences*, pages 1-82. Springer, 1988. [doi], [oai].

[36] A. Muszyńska. Rotor-to-stationary element rub-related vibration phenomena in rotating machinery - literature survey. *The Shock and Vibration Digest*, 21:3-11, 1989. [doi].

[37] A. Muszyńska. *Rotordynamics*, volume 188 of *Mechanical engineering*. Taylor & Francis, 2005. ISBN 978-0-8247-2399-6.

[38] L. Paoli and M. Schatzman. A numerical scheme for impact problems I: The one-dimensional case. *SIAM Journal on Numerical Analysis*, 40(2):702-733, 2002. [doi], [oai].

[39] A. Psarra. *Gas Turbine shaft failure Modelling: friction and wear modelling of turbines in contact*. PhD thesis, Cranfield University, 2010. [hdl].

[40] P. Raveendranath, G. Singh, and B. Pradhan. Free vibration of arches using a curved beam element based on a coupled polynomial displacement field. *Computers & Structures*, 78(4):583-590, 2000. ISSN 0045-7949. [doi].

[41] S. Roques. *Modélisation du comportement dynamique couplé rotor-stator d'une turbine en situation accidentelle*. PhD thesis, Ecole Centrale de Nantes, 2007. [oai].

[42] S. Roques, M. Legrand, P. Cartraud, C. Stoisser, and C. Pierre. Modeling of a rotor speed transient response with radial rubbing. *Journal of Sound and Vibration*, 329(5):527-546, 2010. [doi], [oai].

[43] J. Salençon. *Mécanique des milieux continus : milieux curvilignes*, volume 3. Editions de l'Ecole Polytechnique, 2016. ISBN 978-2-7302-1644-9.

[44] L. Salles. *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, Université Technique d'Etat Bauman de Moscou, 2010. [oai].

[45] N. Salvat. *Modeling of shaft precessional motions induced by unilateral and frictional blade/casing contacts in aircraft engines*. PhD thesis, McGill University, 2015. [oai].

[46] N. Strömberg. Finite element treatment of two-dimensional thermoelastic wear problems. *Computer Methods in Applied Mechanics and Engineering*, 177(3-4):441-455, 1999. [doi], [oai].

[47] A. Thorin, N. Guérin, M. Legrand, F. Thouverez, and P. Almeida. Nonsmooth thermoelastic simulations of blade-casing contact interactions. *Journal of Engineering for Gas Turbines and Power*, 141(2):022502-1-7, 2018. [doi], [oai].

[48] M. Torkhani, L. May, and P. Voinis. Light, medium and heavy partial rubs during speed transients of rotating machines: Numerical simulation and experimental observation. *Mechanical Systems and Signal Processing*, 29:45-66, 2012. ISSN 0888-3270. [doi].

[49] Turbomeca Group. Sealing: Gas turbine engines. Technical report, SNECMA, 2004.

[50] J. Wilkes, D. Childs, B. Dyck, and S. Phillips. The numerical and experimental characteristics of multimode dry-friction whip and whirl. *Journal of Engineering for Gas Turbines and Power*, 132(5):052503-1-9, 2010. ISSN 0742-4795. [doi].

[51] P. Wriggers. Contact constraints within coupled thermomechanical analysis—A finite element model. *Computer Methods in Applied Mechanics and Engineering*, 113(3-4):301-319, 1994. ISSN 0045-7825. [doi].

[52] P. Wriggers and T. Laursen, editors. *Computational Contact Mechanics*, volume 498 of CISM International Centre for Mechanical Sciences. Springer, 2007. ISBN 978-3-211-77297-3 978-3-211-77298-0. [doi].

[53] H. Zhang, X. Li, D. Wang, and T. Liu. Nonlinear dynamics of a blade rotor with coupled rubbing of labyrinth seal and tip seal. *Shock and Vibration*, 2021:1-14, 2021. [doi].
[5] Australian Transport Safety Bureau, 2013. In-flight uncontained engine failure - overhead Batam Island, Indonesia - 4 November 2010 - VH-OQA, Airbus A380-842. Technical Report.

[6] Bartha, A., 2000. Dry friction backward whirl of rotors. Ph.D. thesis. ETH Zurich.

[7] Berthoul, B., Batailly, A., Stainier, L., Legrand, M., Cartraud, P., 2018. Phenomenological modeling of abradable wear in turbomachines. Mechanical Systems and Signal Processing 98, 770-785. doi:10.1016/j.ymssp.2017.05.020. oai:hal-01569316.

[8] Black, H., 1968. Interaction of a whirling rotor with a vibrating stator across a clearance annulus. Journal of Mechanical Engineering Science 10, 1-12. doi:10.1243/JMES_JOUR_1968_010_003_02.

[9] Blevins, R., 1979. Formulas for natural frequency and mode shape. Technical Report. Van Nostrand Reynold Company.

[10] Carpenter, N., Taylor, R., Katona, M., 1991. Lagrange constraints for transient finite element surface contact. International Journal for Numerical Methods in Engineering 32, 103-128. doi:10.1002/nme.1620320107. oai:hal-01389918.

[11] Chen, H.S., Stadtherr, M., 1981. A modification of Powell's dogleg method for solving systems of nonlinear equations. Computers & Chemical Engineering 5, 143-150. doi:10.1016/0098-1354(81)85003-X.

[12] Chen, S., Géradin, M., 1997. Finite element simulation of non-linear transient response due to rotor-stator contact. Engineering Computations 14, 591-603. doi:10.1108/02644409710180356.

[13] Childs, D., Bhattacharya, A., 2007. Prediction of dry-friction whirl and whip between a rotor and a stator. Journal of Vibration and Acoustics 129, 355-362. doi:10.1115/1.2731412.

[14] Choi, Y., 2002. Investigation on the whirling motion of full annular rotor rub. Journal of Sound and Vibration 258, 191-198. doi:10.1006/jsvi.2002.5091.

[15] Choy, F., Padovan, J., 1987. Non-linear transient analysis of rotor-casing rub events. Journal of Sound and Vibration 113, 529-545. doi:10.1016/S0022-460X(87)80135-9.

[16] Dai, X., Jin, Z., Zhang, X., 2002. Dynamic Behavior of the full rotorstop rubbing: numerical simulation and experimental verification. Journal of Sound and Vibration 251, 807-822. doi:10.1006/jsvi.2001.3998.

[17] Delebarre, C., Wagner, V., Paris, J., Dessein, G., Denape, J., Gurt-Santanach, J., 2016. The wear mechanisms occurring in a labyrinth seal/abradable contact depending on the incursion depth parameter. Mechanics & Industry 17, 601. doi:10.1051/meca/2015118. oai:hal-02135086.

[18] Draganis, A., 2017. Finite element modeling of transient thermomechanical rolling contact featuring mixed control of the rigid body motion. Journal of Tribology 139, 011503. doi:10.1115/1.4033048.

[19] Dubigeon, S., 1998. Mécanique des Milieux Continus. Lavoisier, Tec & Doc, Nantes.

[20] Dumartineix, C., 2019. Modélisation et étude de la dynamique complexe d'un système bi-rotor aubagé couplé. Ph.D. thesis. Ecole Centrale de Lyon. oai:tel-03052668.

[21] Duran, C., Manin, L., Andrianoely, M.A., Bordegaray, C., Battle, F., Dufour, R., 2015. Effect of rotor-stator contact on the mass unbalance response, in: Pennacchi, P. (Ed.), 9th IFTOMM International Conference on Rotor Dynamics, pp. 1965-1975. doi:10.1007/978-3-319-06590-8_163. oai:hal-01075974.

[22] Eryilmaz, I., Pachidis, V., 2019. Turbine thermomechanical modelling during excessive axial movement and overspeed. The Aeronautical Journal 123, 248-264. doi:10.1017/aer.2018.162.

[23] Fatarella, F., 1999. On the Dynamics of Reverse Whirl Due to Rotor/Stator Interaction. Ph.D. thesis. Imperial College of Science, Medicine & Technology.

[24] Fried, E., 1969. Thermal conduction contribution to heat transfer at contacts. volume 2 of *Thermal Conductivity*. Tye R.P. Academic Press, London.

[25] Genta, G., 2005. Dynamics of rotating systems. Mechanical engineering series, Springer.

[26] Gonzalez, A., Pachidis, V., 2014. On the numerical simulation of turbine blade tangling after a shaft failure, in: Volume 7B: Structures and Dynamics, American Society of Mechanical Engineers, Düsseldorf, Germany. p. V07BT33A026. doi:10.1115/GT2014-27061.

[27] Géradin, M., Rixen, D., 2014. Mechanical vibrations: theory and application to structural dynamics. 3rd ed., Wiley.

[28] Haake, M., Fiola, R., Staudacher, S., 2011. Multistage compressor and turbine modeling for the prediction of the maximum turbine speed resulting from shaft breakage. Journal of Turbomachinery doi:10.1115/1.4001188.

[29] Ireman, P., Klarbring, A., Strömberg, N., 2002. Finite element algorithms for thermoelastic wear problems. European Journal of Mechanics - A/Solids 21, 423-440. doi:10.1016/S0997-7538(02)01208-1.

[30] Jacobs, C., Legrand, M., Thouverez, F., Almeida, P., 2023. Turbomachinery transient dynamics of radial rotor-stator contact occurrences with friction, in: ASME Turbo Expo, p. V11BT26A007. doi:10.1115/ GT2023-102198. oai:hal-04197339.

[31] Jacquet-Richardet, G., Torkhani, M., Cartraud, P., Thouverez, F., Nouri Baranger, T., Herran, M., Gibert, C., Baguet, S., Almeida, P., Peletan, L., 2013. Rotor to stator contacts in turbomachines. Review and application. Mechanical Systems and Signal Processing 40, 401-420. doi:10.1016/j.ymssp.2013.05.010. oai:hal-00934050.

[32] Jean, M., 1999. The non-smooth contact dynamics method. Computer Methods in Applied Mechanics and Engineering 177, 235-257. doi:10.1016/S0045-7825(98)00383-1. oai:hal-01390459.

[33] Jiang, J., 2009. Determination of the global responses characteristics of a piecewise smooth dynamical system with contact. Nonlinear Dynamics 57, 351-361. doi:10.1007/s11071-008-9446-z.
[34] Lalanne, M., Ferraris, G., 1998. Rotordynamics prediction in engineering. 2 ed., Wiley.

[35] Legrand, M., 2006. Modèles de prédiction de l'interaction rotor/stator dans un moteur d'avion. Ph.D. thesis. Ecole Centrale de Nantes. oai:tel-00011631.

[36] Legrand, M., Pierre, C., Cartraud, P., Lombard, J.P., 2009. Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction. Journal of Sound and Vibration 319, 366-391. doi:10.1016/j.jsv.2008.06.019. oai:hal-00328186.

[37] Lesaffre, N., Sinou, J.J., Thouverez, F., 2007. Contact analysis of a flexible bladed-rotor. European Journal of Mechanics - A/Solids 26, 541-557. doi:10.1016/j.euromechsol.2006.11.002. oai:hal-00322887.

[38] Ma, H., Zhao, Q., Zhao, X., Han, Q., Wen, B., 2015. Dynamic characteristics analysis of a rotor-stator system under different rubbing forms. Applied Mathematical Modelling 39, 2392-2408. doi:10.1016/j.apm.2014.11.009.

[39] Moreau, J.J., 1988. Unilateral contact and dry friction in finite freedom dynamics, in: Moreau, J.J., Panagiotopoulos, P. (Eds.), Nonsmooth Mechanics and Applications. Springer. volume 302 of *International Centre for Mechanical Sciences*, pp. 1-82. doi:10.1007/978-3-7091-2624-0_1. oai:hal-01713847.

[40] Muszyńska, A., 1989. Rotor-to-stationary element rub-related vibration phenomena in rotating machinery – Literature survey. The Shock and Vibration Digest 21, 3-11.

[41] Muszyńska, A., 2005. Rotordynamics. 1 ed., Taylor & Francis.

[42] Paoli, L., Schatzman, M., 2002. A numerical scheme for impact problems I: The one-dimensional case. SIAM Journal on Numerical Analysis 40, 702-733. doi:10.1137/S0036142900378728. oai:hal-01885729.

[43] Psarra, A., 2010. Gas Turbine shaft failure Modelling: friction and wear modelling of turbines in contact. Ph.D. thesis. Cranfield University. URL: http://dspace.lib.cranfield.ac.uk/handle/1826/8348.

[44] Raveendranath, P., Singh, G., Pradhan, B., 2000. Free vibration of arches using a curved beam element based on a coupled polynomial displacement field. Computers & Structures 78, 583-590. doi:10.1016/S0045-7949(00)00038-9.

[45] Roques, S., 2007. Modélisation du comportement dynamique couplé rotor-stator d'une turbine en situation accidentelle. Ph.D. thesis. Ecole Centrale de Nantes. oai:tel-00353938.

[46] Roques, S., Legrand, M., Cartraud, P., Stoisser, C., Pierre, C., 2010. Modeling of a rotor speed transient response with radial rubbing. Journal of Sound and Vibration 329, 527-546. doi:10.1016/j.jsv.2009.09.016. oai:hal-00431223.

[47] Salençon, J., 2016. Mécanique des milieux continus : milieux curvilignes. volume 3. Editions de l'Ecole Polytechnique.

[48] Salles, L., 2010. Etude de l'usure par fretting sous chargements dynamiques dans les interfaces frottantes: application aux pieds d'aubes de turbomachines. Ph.D. thesis. Ecole Centrale de Lyon, Université Technique d'Etat Bauman de Moscou. oai:tel-00600613.

[49] Salvat, N., 2015. Modeling of shaft precessional motions induced by unilateral and frictional blade/casing contacts in aircraft engines. Ph.D. thesis. McGill University. oai:tel-01259136.

[50] Soria, C., 2014. Gas Turbine Shaft Over-speed/Failure Modelling. Aero/Thermodynamics Modelling and Overall Engine System Response. Ph.D. thesis. Cranfield University. URL: http://dspace.lib.cranfield.ac.uk/handle/1826/12131.

[51] Strömberg, N., 1999. Finite element treatment of two-dimensional thermoelastic wear problems. Computer Methods in Applied Mechanics and Engineering 177, 441-455. doi:10.1016/S0045-7825(98)00392-2. oai:hal-01574166.

[52] Thorin, A., Guérin, N., Legrand, M., Thouverez, F., Almeida, P., 2018. Nonsmooth thermoelastic simulations of blade-casing contact interactions. Journal of Engineering for Gas Turbines and Power 141. doi:10.1115/1.4040857. oai:hal-01989188.

[53] Torkhani, M., May, L., Voinis, P., 2012. Light, medium and heavy partial rubs during speed transients of rotating machines: Numerical simulation and experimental observation. Mechanical Systems and Signal Processing 29, 45-66. doi:10.1016/j.ymssp.2012.01.019.

[54] Turbomeca Snecma Group, 2004. Sealing.

[55] Wilkes, J., Childs, D., Dyck, B., Phillips, S., 2010. The numerical and experimental characteristics of multimode dry-friction whip and whirl. Journal of Engineering for Gas Turbines and Power 132. doi:10.1115/1.3204658.

[56] Wriggers, P., 1994. Contact constraints within coupled thermomechanical analysis—A finite element model. Computer Methods in Applied Mechanics and Engineering 113, 301-319. doi:10.1016/0045-7825(94) 0051-5.

[57] Wriggers, P., Laursen, T. (Eds.), 2007. Computational Contact Mechanics. volume 498 of CISM International Centre for Mechanical Sciences. Springer. doi:10.1007/978-3-211-77298-0.

[58] Zhang, H., Li, X., Wang, D., Liu, T., 2021. Nonlinear dynamics of a blade rotor with coupled rubbing of labyrinth seal and tip seal. Shock and Vibration 2021, 1-14. doi:10.1155/2021/3826429.