# Modeling of a rotor speed transient response with radial rubbing

Sébastien Roques, Mathias Legrand, Patrice Cartraud, Carlo Stoisser, Christophe Pierre

► To cite this version:

Sébastien Roques, Mathias Legrand, Patrice Cartraud, Carlo Stoisser, Christophe Pierre. Modeling of a rotor speed transient response with radial rubbing. Journal of Sound and Vibration, 2010, 329 (5), pp.527-546. 10.1016/j.jsv.2009.09.016. hal-00431223

HAL Id: hal-00431223

https://hal.science/hal-00431223v1

Submitted on 10 Nov 2009

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 CC BY 4.0 - Attribution - International License
# Modeling of a rotor speed transient response with radial rubbing

**Sébastien Roques**<fnref n="1" />

Institut de Recherche en Génie civil et Mécanique (GeM), UMR CNRS 6183, École Centrale de Nantes, BP 92101, 44321 Nantes cedex 3, France

EDF Research & Development, Acoustics and Mechanical Analysis Department, 1 avenue du Général de Gaulle, 92141 Clamart, France

**Mathias Legrand**

Structural Dynamics and Vibration Laboratory, Department of Mechanical Engineering, McGill University, 817 Sherbrooke St West, Montreal, Quebec H3A 2K6, Canada

**Patrice Cartraud**

Institut de Recherche en Génie civil et Mécanique (GeM), UMR CNRS 6183, École Centrale de Nantes, BP 92101, 44321 Nantes cedex 3, France

**Carlo Stoisser**

EDF Research & Development, Acoustics and Mechanical Analysis Department, 1 avenue du Général de Gaulle, 92141 Clamart, France

**Christophe Pierre**

Structural Dynamics and Vibration Laboratory, Department of Mechanical Engineering, McGill University, 817 Sherbrooke St West, Montreal, Quebec H3A 2K6, Canada

## Abstract

A rotor-stator model of a turbogenerator is introduced in order to investigate speed transients with rotor-to-stator rubbing caused by an accidental blade-off imbalance. In order to assess the angular deceleration of the rotor due to rubbing, the angular position of its cross-section is considered as an unknown of the problem. Displacement fields are discretized through a finite element formulation. The highly nonlinear equations due to contact conditions are solved through an explicit prediction-correction time-marching procedure combined with the Lagrange multiplier approach dealing with a node-to-line contact strategy. The developed numerical tool is suitable for analyzing rotor-stator interactions in turbomachines as the system passes through critical speeds during an accidental shutdown. The sensitivity of the system response to modeling, physical and numerical parameters is investigated. The results highlight the significant role of the friction coefficient together with the diaphragm modeling, from rigid to fully flexible, in the interaction phenomenon. Rigid models have the advantage of simplicity and provide reasonable estimations of the overall response of the turbine. A flexible model, however, may be more computationally intensive but is more appropriate in order to accurately capture quantities of interest such as shaft eccentricity and bearing loads.

**Keywords:** Speed transient, rotor, stator, contact, rubbing, finite element method, explicit time-marching technique

# 1 Introduction

In nuclear power plant turbosets similar to the one pictured in figure 1, the reference design-basis accident consists of a blade-off in the last stage of the low pressure turbine. During an accidental shutdown,

<footnote n="1">¹Corresponding author: sebastien.roques@edf.fr</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 1: A real power plant turboset. Courtesy of Électricité de France</figcaption></figure>

a severe rotor-stator interaction may occur at critical speeds due to large shaft line displacements orig-
inated by the lost blade high imbalance excitation. The potentially induced rubbing between the shaft
and the stator, as illustrated in figure 2, leads to an important angular deceleration rate and highly modi-
fies the dynamics of the turbogenerator. Risks of failure of the contacting components may occur due to
the heavy friction torque. It is then of primary importance to ensure that the adopted turbine design is
capable of going through critical speeds without catastrophic consequences for the shaft line.

Rubbing, known to be a serious malfunction in turbomachinery, has been the subject of a large
amount of research and a detailed overview is provided in [25, 23]. First mathematical models dedi-
cated to rubbing issues were as simple as Jeffcott rotors [28]. They were then extended to flexible rotors
through finite element approaches and/or modal synthesis techniques [15] allowing for more realistic
descriptions. However, these studies were limited to constant angular velocity steady states, which im-
plies an increase of the driving torque in case of rubbing, for instance. In [15], equations of motion were
rewritten considering a non-constant but known angular velocity law. Only Dai [11] has proposed a rigid
rotor model where an additional equation computes the instantaneous rotational velocity depending on
the contact forces.

The present work is a preliminary study of a design-basis accident in a nuclear power plant. A large-scale turbine operating at nominal conditions is suddenly disconnected from the electrical network after a blade-off. The main objective of this work is to suggest a predictive numerical tool in order to assess the ability of the designed turbine to go through critical speeds without catastrophic consequences. The present paper focuses on an accidental shutdown originated by a blade-off yielding radial rubbing between the shaft and the diaphragm as displayed in figure 2. Quite similar situations involving blade-tips to casing interaction have already been investigated in [18, 20]. Nevertheless, it is here assumed that the latter is negligible compared to the shaft-diaphragm interaction. Previous works [21, 24, 7, 8, 29, 2, 10, 22, 11] yield the following classification of rubbing phenomena:

* The duration of the contact involves either partial, i.e. intermittent, or full, i.e. continuous, annular rubbing. Stable vibrations with synchronous, sub-synchronous or super-synchronous responses and chaotic behaviors of the shaft-line with destructive damages of its components can be observed.

* The rotor features forward and reverse whirl motions, as well as oil and dry whip depending on running conditions.

* The contact load can cause heavy or light rubbing.
A reliable description of rubbing is only possible with a detailed finite element model together with a numerical approach capable of tackling transient dynamics and contact constraints. In that view, this paper contributes to a better description of shaft dynamics by incorporating rubbing with unknown angular velocity mainly affected by aerodynamical and contact forces. The related model includes gyroscopic effects and torsional displacements. The stator model is restricted in this study to the diaphragm part with an increasing level of complexity, from rigid to fully flexible.

<figure><img src="image_3.png" /><figcaption>FIGURE 2: Contact locations in a real turbine</figcaption></figure>

In the governing equations, the dynamics of the rotor and diaphragm are coupled through contact forces. Available approaches in contact mechanics [4] are usually the penalty method, the Lagrange multiplier method or the augmented Lagrangian method [31, 16, 33]. An accurate contact treatment strategy is strongly connected to the integration method. As our interest lies in the characterization of the transient response of the shaft during rubbing, approaches such as harmonic balance or shooting methods, more intended to forced and steady-state responses are discarded, and time marching procedures stand as a natural choice. Two main families are usually found in the literature: implicit versus explicit formulations [26, 27]. It is observed in [14] that the results highly depend on Newmark parameters for problems involving strong nonlinear terms such as direct contact constraints. Even though implicit methods with automatic control of the time step size [6] may be used, conditionally stable explicit methods appear to be more relevant to non-differentiable terms such as contact. Accordingly, the well-known central finite difference scheme is adopted here. A previous study [18] has shown that the penalty method is not always adapted to contact detection in crash analysis since residual penetrations between components is allowed: cumulative error and dependency of the results to the penalty parameter may be expected. Lagrange multipliers and augmented Lagrange multipliers are prone to Uzawa type algorithms that do not fit in explicit time marching techniques because of the required CPU. Nevertheless, the prediction-correction algorithm forward increment Lagrange method developed in [32] embeds the Lagrange multiplier approach in an explicit technique keeping the advantages of both. It has been proved reliable for contact-impact problems [1] by properly satisfying contact detection and ensuring displacement compatibility and is therefore preferred in this study.

## 2 Rotor speed transient modeling

The turbine model introduced in this work consists in Timoshenko beams for the rotating flexible shaft complemented by circular rigid disks for the bladed disks. Imbalances are originated by concentrated masses and the bearing behavior is linearized. It is assumed that the rotor is simply supported at its ends. This model is pictured in figure 3 and each of its components is thoroughly detailed in the sequel.

The geometrical and material parameters of the shaft line are chosen so that the first eigenfrequencies of the rotor match the corresponding ones of the real structures as explained later.
<figure><img src="image_4.png" /><figcaption>FIGURE 3: FE model schematic</figcaption></figure>

## 2.1 Shaft line components

An accidental blade-off in a nuclear power plant turboset leads to a controlled shutdown of the shaft line and subsequent high deceleration rates due to fluids friction and shaft-to-diaphragm resisting torques. Accordingly, attention must be paid to the modeling and formulation of such phenomena since the angular position of a cross-section becomes an unknown of the problem. In-plane as well as out-of-plane bending vibrations, torsional vibrations together with axial vibrations are considered. Gyroscopic terms are included and the following assumptions are also accounted for:

1. The shaft has a uniform circular cross-section along its length.

2. The shaft is initially balanced.

3. External torques with constant direction along the reference centroidal axis are applied at each end of the shaft.

### 2.1.1 Rigid disk

The kinetic energy $T_{rd}$ of a rigid circular disk rotating at angular velocity $\dot{\varphi}_r$ takes the form:

$$2T_{rd} = M_{rd} [\dot{u}_r^2 + \dot{v}_r^2 + \dot{w}_r^2] + I_{rd_x} [\dot{\theta}_{u_r}^2 \cos^2 \theta_{v_r} + \dot{\theta}_{v_r}^2] \\ + I_{rd_z} [\dot{\theta}_{u_r}^2 \sin^2 \theta_{v_r} + (\dot{\varphi}_r + \dot{\beta}_r)^2 + 2(\dot{\varphi}_r + \dot{\beta}_r) \dot{\theta}_{u_r} \sin \theta_{v_r}] \quad (1)$$

where $M_{rd}$, $I_{rd_x}$ and $I_{rd_z}$ respectively stand for the mass, the second moment of inertia and the polar moment of inertia of the considered rigid disk. Referring to figure 4, the notations used for the rotor kinematics are $(u_r, \theta_{v_r})$ for the bending displacements in the $(e_x, e_z)$-plane, $(v_r, \theta_{u_r})$ for the bending displacements in the $(e_y, e_z)$-plane, $w_r$ for the axial displacement and $\beta_r$ for the torsional twist, both along $e_z$. These quantities are considered in equation (1) at the disk center. In the small perturbation framework, the kinetic energy, including gyroscopic effects, is rewritten as follows:

$$2T_{rd} = M_{rd} [\dot{u}_r^2 + \dot{v}_r^2 + \dot{w}_r^2] + I_{rd_x} [\dot{\theta}_{u_r}^2 + \dot{\theta}_{v_r}^2] + I_{rd_z} [(\dot{\varphi}_r + \dot{\beta}_r)^2 + 2(\dot{\varphi}_r + \dot{\beta}_r) \dot{\theta}_{u_r} \theta_{v_r}] \quad (2)$$

### 2.1.2 Shaft

The kinetic energy $T_s$ of the shaft corresponds to the integration of the disk energy along its longitudinal direction thus yielding:

$$2T_s = \rho_s \int_0^{l_s} [S_s (\dot{u}_r^2 + \dot{v}_r^2 + \dot{w}_r^2) + I_{sx} (\dot{\theta}_{u_r}^2 + \dot{\theta}_{v_r}^2) + I_{sp} ((\dot{\varphi}_r + \dot{\beta}_r)^2 + 2(\dot{\varphi}_r + \dot{\beta}_r) \dot{\theta}_{u_r} \theta_{v_r})] dz \quad (3)$$
<figure><img src="image_5.png" /><figcaption>FIGURE 4: Notations for a disk (EDF convention for Euler angles)</figcaption></figure>

where $\rho_s$, $S_s$, $I_{sx}$ and $I_{sp}$ respectively refer to the mass density, the cross-section area, the moment of inertia and the polar moment of inertia of the shaft.

The potential energy $U_s$ of a spinning Timoshenko beam is equal to:

$$2U_s = \int_0^{I_s} \left[ E_s S_s w_{r,z}^2 + E_s I_{sx} (\theta_{u_r,z}^2 + \theta_{v_r,z}^2) + G_s I_{sp} \beta_{r,z}^2 + k_s G_s S_s ((u_{r,z} - \theta_{u_r})^2 + (v_{r,z} + \theta_{v_r})^2) \right] dz \quad (4)$$

where $E_s$, $G_s$ and $k_s$ are respectively the Young modulus, the shear modulus and the transverse shear form factor of the shaft.

### 2.1.3 Linearized bearing

In real turbines, oil film bearings support the shaft and involve nonlinear terms described by Reynolds' equations. The virtual work $\delta W$ of these external nonlinear forces acting on the shaft, $F_{u_r}$ and $F_{v_r}$ along $u_r$ and $v_r$ respectively, can be expressed in a general manner as:

$$\delta W = [F_{u_r} \ F_{v_r}] \begin{pmatrix} \delta u_r \\ \delta v_r \end{pmatrix} \qquad (5)$$

For small displacements with respect to the equilibrium position of the shaft in the bearings, stiffness and damping coefficients of the oil film, respectively $k_{o_{ij}}$ and $d_{o_{ij}}$, can be calculated by linearizing Reynolds' equations [17] leading to:

$$\begin{pmatrix} F_{u_r} \\ F_{v_r} \end{pmatrix} = - \begin{bmatrix} k_{o_{xx}} & k_{o_{xy}} \\ k_{o_{yx}} & k_{o_{yy}} \end{bmatrix} \begin{pmatrix} u_r \\ v_r \end{pmatrix} - \begin{bmatrix} d_{o_{xx}} & d_{o_{xy}} \\ d_{o_{yx}} & d_{o_{yy}} \end{bmatrix} \begin{pmatrix} \dot{u}_r \\ \dot{v}_r \end{pmatrix} \qquad (6)$$
### 2.1.4 Imbalance

In order to reflect a blade-off, a heavy imbalance is introduced through a concentrated mass $m_i$ located at a distance $r_i$ from the geometric center of the shaft. Its kinetic energy is:

$$2T_i = m_i [\dot{u}_r^2 + \dot{v}_r^2 + \dot{w}_r^2 + r_i^2 \phi_r^2 + 2r_i \phi_r (\dot{v}_r \cos \phi_r - \dot{u}_r \sin \phi_r)] \quad (7)$$

## 2.2 Diaphragm

<figure><img src="image_6.png" /><figcaption>FIGURE 5: Schematic representation of a real diaphragm</figcaption></figure>

As mentioned previously, rotor-to-stator contact corresponds to a shaft-diaphragm interaction in this study. Therefore, the stator modeling reduces to the diaphragm representation. The real diaphragm or flow straightener comprises inner and outer rings and twisted blades as depicted in figure 5 and three different types of diaphragm have been developed:

1. Diaphragm D1 is shown in figure 6(a): it is fully rigid and can be seen as a mathematical boundary. Its main purpose is the validation of the contact algorithm.

2. Diaphragm D2 is a rigid ring with flexibility and damping as depicted in figure 6(b). This model is an extension of the first model with a formulation based on a previous EDF study [9]: it is augmented with stiffnesses, dampers and the mass of the inner ring through a mass matrix $m_d$.

3. Diaphragm D3 is a flexible structure illustrated in figure 6(c). The inner ring is now flexible with the use of curved beams and the blades are considered as straight beams.

<figure><img src="image_7.png" /><figcaption>FIGURE 6: Investigated diaphragm models: (a) diaphragm D1, (b) diaphragm D2, (c) diaphragm D3</figcaption></figure>
### 2.2.1 Inner ring: curved beams

The inner ring, initially in the (e<sub>x</sub>, e<sub>y</sub>)-plane, is depicted in figure 7. Dimensions of the cross-section are small compared to the mean line length, so that Euler-Bernoulli theory holds. Strains are then written as follows [12]:

$$ 
\begin{aligned}
\epsilon_{ss} &= \frac{1}{r_{ir} + r} (v_{ir} + (r_{ir} + r)u_{ir,s} - r r_{ir} v_{ir,ss} + z (\theta_{ir,s} - r_{ir} w_{ir,ss})) \\
\gamma_{rs} &= \frac{z}{r_{ir} + r} (r_{ir} \theta_{ir,s} + w_{ir,s}) \\
\gamma_{zs} &= \frac{-r}{r_{ir} + r} (r_{ir} \theta_{ir,s} + w_{ir,s})
\end{aligned} 
\quad (8) $$

where $u_{ir}$, $v_{ir}$ and $w_{ir}$ (resp. $\theta_{ir,s}$, $\theta_{ir,r}$ and $\theta_{ir,z}$) respectively stand for radial, tangential and axial displacements (resp. rotations) of the centroidal line as displayed in figure 7. The path variable is $s$ and $r_{ir}$ is the

<figure><img src="image_8.png" /><figcaption>FIGURE 7: Notations for a curved beam</figcaption></figure>

average radius of the inner ring. Assuming $r \ll r_{ir}$ and a first-order expansion in $r$, $\epsilon_{ss}$ in equation (8) becomes:

$$ \epsilon_{ss} = \left( \frac{v_{ir}}{r_{ir}} + u_{ir,s} \right) - r \left( \frac{v_{ir}}{r_{ir}^2} + v_{ir,ss} \right) + z \left( \frac{\theta_{ir,s}}{r_{ir}} - w_{ir,ss} \right) + O(r^2) \quad (9) $$

Assuming a linear elastic behavior for the diaphragm, i.e. $\sigma_{ss} = E_d \epsilon_{ss}$, $\tau_{rs} = G_d \gamma_{rs}$ and $\tau_{zs} = G_d \gamma_{zs}$ (Hooke's law) where $E_d$ and $G_d$ are the Young's and shear moduli together with a uniform rectangular cross-section of the curved beam, the strain energy $U_{ir}$ of a finite element becomes:

$$ 2U_{ir} = E_d \int_{s_1}^{s_2} \left[ S_{ir} \left( \frac{v_{ir}}{r_{ir}} + u_{ir,s} \right)^2 + I_{irz} \left( \frac{v_{ir}}{r_{ir}^2} + v_{ir,ss} \right)^2 + I_{irr} \left( \frac{\theta_{ir,s}}{r_{ir}} - w_{ir,ss} \right)^2 + \frac{G_d J_{ir}}{E_d} \left( \theta_{ir,s,s} + \frac{w_{ir,s}}{r_{ir}} \right)^2 \right] ds \quad (10) $$

where $I_{irr}$, $I_{irz}$, $J_{ir}$ and $S_{ir}$ respectively stand for the two moments of inertia, the polar moment of inertia and the cross-section area of the inner ring. Similarly, the kinetic energy $T_{ir}$ is defined as:

$$ 2T_{ir} = \rho_d \int_{s_1}^{s_2} \left[ S_{ir} (\dot{u}_{ir}^2 + \dot{v}_{ir}^2 + \dot{w}_{ir}^2) + I_{irr} \dot{w}_{ir,s}^2 + I_{irz} \left( \frac{\dot{u}_{ir}}{r_{ir}} - \dot{v}_{ir,s} \right)^2 + I_{irp} \dot{\theta}_{ir,s}^2 \right] ds \quad (11) $$

with $I_{irp}$ denoting the moment of inertia.
<figure><img src="image_9.png" /><figcaption>FIGURE 8: Notations for a straight beam in its local reference frame</figcaption></figure>

### 2.2.2 Blades: straight beams

The Euler Bernoulli theory is also used for the blades as slender structures. Their energies collapse to equations (3) and (4) where the angular velocity is equal to zero so that the shear energy vanishes. Consequently, the notations introduced in figure 8 together with the two usual kinematic relationships $\theta_{w_b} = v_{b,x}$ and $\theta_{v_b} = -w_{b,x}$ yield:

$$ 
\begin{aligned}
2T_b &= \rho_d \int_0^{l_b} \left[ S_b (\dot{u}_b^2 + \dot{v}_b^2 + \dot{w}_b^2) + I_{b_z} \dot{\theta}_{w_b}^2 + I_{b_y} \dot{\theta}_{v_b}^2 + I_{b_x} \dot{\theta}_{u_b}^2 \right] dx \\
2U_b &= \int_0^{l_b} \left[ E S_b u_{b,x}^2 + E I_{b_z} \dot{\theta}_{w_b,x}^2 + E I_{b_y} \dot{\theta}_{v_b,x}^2 + G_d J_b \dot{\theta}_{u_b,x}^2 \right] dx
\end{aligned}
\quad (12) $$

where $I_{b_z}$, $I_{b_y}$, $I_{b_x}$, $J_b$ and $S_b$ respectively stand for the two second moments of inertia, the polar moment of inertia, the cross-sectional polar moment of inertia and the cross-section area of the rectangular straight beam.

### 2.2.3 Curved beam–straight beam compatibility

The assembling procedure between the straight beam finite elements of the blade and curved beam finite elements of the inner ring at the connecting nodes (red and green finite elements in figure 6(c)) in order to compute the global matrices of the system of interest requires a special attention. This procedure stems from the compatibility conditions along the generalized displacements and forces to be satisfied at each node. As shown in figures 7 and 8, in the present study these conditions reduce to:

$$ 
\begin{aligned}
u_b &= v_{ir} & \theta_{u_b} &= \theta_{ir_r} \\
v_b &= u_{ir} & \theta_{v_b} &= \theta_{ir_s} \\
w_b &= w_{ir} & \theta_{w_b} &= \theta_{ir_z}
\end{aligned}
\quad (13) $$

From equation (13), the rotation of a cross–section of a Euler–Bernoulli curved beam yields:

$$ \theta_{w_b} = \frac{u_{ir}}{r_{ir}} - v_{ir,s} \quad \text{and} \quad \theta_{u_b} = w_{ir,s} \quad (14) $$

condition that has to be accounted for in the finite element discretization.

## 3 Finite element discretization

The different energies detailed above are discretized according to the finite element technique formulated in displacement. The elements used are detailed hereafter.

### 3.1 Shaft line

Traction and torsion are discretized using the usual linear shape functions. Denoting $\xi = \frac{z}{l}$ and $\mathbf{q} = (w_r, \beta_r)$ yields:

$$ \mathbf{q}(z) = [1 - \xi \xi] \begin{pmatrix} \mathbf{q}_1 \\ \mathbf{q}_2 \end{pmatrix} \quad (15) $$
The modified Hermite functions $N_{ij}$ used in [5] for bending are selected. Degrees of freedom $u_r(z)$, $\theta_{v_r}(z)$, $v(z)$ and $\theta_{u_r}(z)$ are thus written as follows:

$$ \begin{pmatrix} u_r(z) \\ \theta_{v_r}(z) \end{pmatrix} = \begin{bmatrix} N_{11}(z) & N_{12}(z) & N_{13}(z) & N_{14}(z) \\ N_{21}(z) & N_{22}(z) & N_{23}(z) & N_{24}(z) \end{bmatrix} \begin{pmatrix} u_{r1} \\ \theta_{v_{r1}} \\ u_{r2} \\ \theta_{v_{r2}} \end{pmatrix} \qquad (16) $$

and:

$$ \begin{pmatrix} v_r(z) \\ \theta_{u_r}(z) \end{pmatrix} = \begin{bmatrix} N_{11}(z) & -N_{12}(z) & N_{13}(z) & -N_{14}(z) \\ N_{21}(z) & -N_{22}(z) & N_{23}(z) & -N_{24}(z) \end{bmatrix} \begin{pmatrix} v_{r1} \\ \theta_{u_{r1}} \\ v_{r2} \\ \theta_{u_{r2}} \end{pmatrix} \qquad (17) $$

Since the angular position $\phi_r$ is unknown, additional boundary conditions are necessary and detailed in section 4.3.

## 3.2 Diaphragm D3

Equations of motion are derived in a very general three-dimensional fashion. However, a numerical investigation conducted in a preliminary phase with a planar diaphragm model showed that conditions $w_{ir} = \theta_{ir_s} = 0$ for the inner ring and $w_b = \theta_{u_b} = \theta_{v_b} = 0$ for the blades hold in our study. The blades are discretized as usual straight Euler-Bernoulli beams FE with three degrees of freedom (dof) per node. The inner ring is discretized with curved beam FE with four dof per node ($u_{ir}$, $u_{ir,s}$, $v_{ir}$, $v_{ir,s}$) with cubic polynomials in $u_{ir}$ and $v_{ir}$. Defining $\zeta = \frac{s}{l_{ire}}$, where $l_{ire}$ is the curvilinear length of the finite element, the shape functions are:

$$ \begin{aligned} N_1(s) &= 1 - 3\zeta^2 + 2\zeta^3 & N_2(s) &= l_{ire}\zeta(1 - 2\zeta + \zeta^2) \\ N_3(s) &= \zeta(3\zeta^2 - 2\zeta^3) & N_4(s) &= l_{ire}\zeta(-\zeta + \zeta^2) \end{aligned} \qquad (18) $$

The discretized displacement field becomes, for an element whose nodes are 1 and 2:

$$ \begin{aligned} u_{ir}(s) &= N_1(s)u_{ir1} + N_2(s)u_{ir1,s} + N_3(s)u_{ir2} + N_4(s)u_{ir2,s} \\ v_{ir}(s) &= N_1(s)v_{ir1} + N_2(s)v_{ir1,s} + N_3(s)v_{ir2} + N_4(s)v_{ir2,s} \end{aligned} \qquad (19) $$

# 4 Contact forces

The forces of particular interest in this study are the unilateral contact and friction forces acting between the shaft and the diaphragm. In order to simplify the contact detection, it is assumed that only one point of the beam cross-section area of the rotor comes into contact with one curved element of the diaphragm, as depicted in figure 10.

## 4.1 General framework

To derive the contact dynamic equations, the master-slave approach [4] is used. Figure 9 introduces the notations: the configurations of the master and slave bodies are denoted $\Omega^{(m)}$ and $\Omega^{(s)}$ with corresponding boundaries $\Gamma^{(m)}$ and $\Gamma^{(s)}$ and contact interface $\Gamma_c$. The gap function $g(\mathbf{x})$ between any point $\mathbf{x}$ of the master component and its closest counterpart $\bar{\mathbf{x}}$ of the slave one can be computed such as:

$$ g(\mathbf{x}) = g_0(\mathbf{x}) + (\mathbf{u}^{(m)}(\mathbf{x}) - \mathbf{u}^{(s)}(\bar{\mathbf{x}})) \cdot \mathbf{n} \qquad (20) $$

where $g_0$ denotes the initial gap, $\mathbf{n}$ the unit outward vector normal to the master surface, $\mathbf{u}^{(m)}(\mathbf{x})$ and $\mathbf{u}^{(s)}(\bar{\mathbf{x}})$ displacements respective to the master and slaved structures. Every material point $\mathbf{x} \in \Gamma_c$ must satisfy the following Kuhn-Tucker optimality conditions:

$$ \lambda_N \ge 0, \quad g(\mathbf{x}) \ge 0 \quad \text{and} \quad \lambda_N \cdot g(\mathbf{x}) = 0 \qquad (21) $$
<figure><img src="image_10.png" /><figcaption>FIGURE 9: Notations for a contact problem formulation</figcaption></figure>

where $\lambda_N$ stands for the positive contact force acting on the slave surface in the normal direction. These unilateral contact conditions are augmented with the Coulomb friction law assuming only sliding occurs:

$$ \|\lambda_T\| = \mu|\lambda_N| \Rightarrow \exists \alpha \in \mathbb{R} \neq 0 \text{ such as } \mathbf{v}_T = \alpha \frac{\lambda_T}{\|\lambda_T\|} \qquad (22) $$

for which $\mu$ is the coefficient of friction, $\mathbf{v}_T$ the tangential slip velocity and $\lambda_T$ the contact force acting along the tangential direction. Since the finite element method is applied, the equations of motion of a mechanical model respective to contact dynamics can be written in a very general manner as follows:

$$ \mathbf{M}\ddot{\mathbf{u}} + \mathbf{D}\dot{\mathbf{u}} + \mathbf{K}\mathbf{u} + \mathbf{C}^T\lambda = 0 \qquad (23) $$

where $\mathbf{M}, \mathbf{D}, \mathbf{K}$ and $\mathbf{u}$ respectively stand for the mass, damping, stiffness matrices and displacement vector of the global system. Lagrange multipliers are stored in vector $\lambda$ and $\mathbf{C}$ is the contact constraint matrix in the normal and tangential directions. Embedded in an explicit time marching approach, a solution method for solving equation (23) together with conditions (21) and (22) is proposed in [32, 19].

## 4.2 Application to rotor-diaphragm system

By choice in this study, the master component is the rotor whereas the diaphragm is the slave structure. Within the small perturbation framework with planar diaphragm models, it is legitimately assumed that a cross-section of the shaft always remains in the $(\mathbf{e}_x, \mathbf{e}_y)$ plane. Accordingly, the determination of $\bar{\mathbf{x}}$ of equation (20) can then be obtained through the following considerations:

1. diaphragms D1 or D2: the contact occurs at hot spot C in figure 10, i.e. highest rotor eccentricity. Consequently the gap function in the radial direction is equal to the distance $\bar{\mathbf{C}}\bar{\mathbf{D}}$:

$$ \bar{\mathbf{C}}\bar{\mathbf{D}} = \mathbf{O}^7\bar{\mathbf{D}} - \mathbf{O}^7\bar{\mathbf{C}} \qquad (24) $$

2. diaphragm D3: in figure 10, contact on the diaphragm occurs at point D with same angular position as point C. This approximation allows for CPU time savings together with accurate predictions of rotational velocities and bearing loads [30].

## 4.3 Boundary conditions and initial conditions

The initial conditions are given as follows: operating at normal conditions, the turbine is suddenly disconnected after the blade-off. The latter gives rise to a heavy mass imbalance while no driving torque holds anymore. Only aerodynamical and Newtonian fluid frictions act on the shaft thus initiating a slow-down with possible contact interaction:
<figure><img src="image_11.png" /><figcaption>FIGURE 10: Contact detection</figcaption></figure>

* newtonian fluid friction torque: $C_{newt} = -A_{newt} \dot{\varphi}_r$;

* aerodynamical friction torque: $C_{aero} = -A_{aero} \dot{\varphi}_r^2$.

where $A_{newt}$ and $A_{aero}$ coefficients were identified in [13]. Fluids forces acting on the turbine are uniformly distributed along the shaft as a first approach. The additional boundary conditions relative to the angular position are:

* Constant driving torque $C_{max} \ge 0$ at one end of the shaft.

* $\dot{\varphi}$-dependent alternator resisting torque $-C_{max} \frac{\dot{\varphi}_r}{\Omega}$ at the other end of the shaft.

It is noteworthy saying that $C_{max} > 0$ refers to the run-up of the shaft line whereas $C_{max} = 0$ refers to the accidental shutdown.

# 5 Governing equations

## 5.1 Rotor

Energy-based Hamilton's principle is used for the derivation of the system's equations of motion. Gyroscopic terms ($\dot{\varphi}_r + \dot{\beta}_r \dot{\theta}_{u_r} \theta_{v_r}$ in equation (2) give rise to a strong nonlinear coupling between the flexural, torsional displacements and the angular position. Similarly, the term ($\dot{\varphi}_r + \dot{\beta}_r^2$) highly couples the angular position to the torsion angle whereas the contact constraints add even more complexity to the system dynamics. Neglecting torsional vibration, the governing equations are as follows:

$$ \mathbf{M}_1 \ddot{\mathbf{X}}_r + (\mathbf{D}_o + \dot{\varphi}_r \mathbf{D}_2) \dot{\mathbf{X}}_r + (\mathbf{K}_o + \mathbf{K}_1 + \ddot{\varphi}_r \mathbf{K}_2) \mathbf{X}_r + \mathbf{C}_r^T \mathbf{\lambda} = + \sum_{j=1}^{N_j} m_{i_j} r_{i_j} \left[ (\ddot{\varphi}_r \sin \varphi_r + \dot{\varphi}_r^2 \cos \varphi_r) \mathbf{F}_{|u_r} + (-\ddot{\varphi}_r \cos \varphi_r + \dot{\varphi}_r^2 \sin \varphi_r) \mathbf{F}_{|v_r} \right] \quad (25) $$

where $\mathbf{X}_r$ is the vector of the rotor generalized displacements (bending, traction and torsion when considered), and $\mathbf{M}_1, \mathbf{D}_2, \mathbf{K}_1, \mathbf{K}_2$ respectively the mass, gyroscopic and two stiffness matrices. Matrix $\mathbf{C}_r$ stands for the rotor contribution to the contact matrix $\mathbf{C}$. The linearized forces acting in the oil-film bearings appear through matrices $\mathbf{K}_o$ and $\mathbf{D}_o$. The right-hand side of equation (25) stand for the imbalance forces
originated by Nᵢ eccentric masses. Finally, the scalar equation governing the angular position φᵣ takes the form:

$$ \left( \rho_s I_{sp} l_s + \sum_{j=1}^{N_i} m_{ij} r_{ij}^2 + \sum_{j=1}^{N_{rd}} I_{rd} z_{j} \right) \ddot{\varphi}_r + \ddot{\mathbf{X}}_r^T \bar{\mathbf{M}} \mathbf{X}_r + \dot{\mathbf{X}}_r^T \bar{\mathbf{M}} \dot{\mathbf{X}}_r + \sum_{j=1}^{N_i} m_{ij} r_{ij} \left( \ddot{u}_r^i \cos \varphi_r - \ddot{u}_r^i \sin \varphi_r \right) \Big|_{\|j\} = F_{\|\varphi_r} = C(t) - C_f \quad (26) $$

where Nᵢ is the number of rotor imbalances and uᵣᵢ and vᵣᵢ stand for their nodal vertical and horizontal displacements. Nrd is the number of rigid disks, C(t) corresponds to the driving torque and Cf to the sum of the fluids friction and contact friction torques whose expression is Cfric = μrₛλN, where rₛ is the outer radius of the shaft. When torsion is considered, terms that cannot be written in a matrix form and not mentioned here for the sake of simplicity, are added to the governing equations.

## 5.2 Diaphragm

The equation of motion respective to the diaphragm dynamics (diaphragm D2 or D3) is:

$$ \mathbf{M}_d \ddot{\mathbf{X}}_d + \mathbf{D}_d \dot{\mathbf{X}}_d + \mathbf{K}_d \mathbf{X}_d - \mathbf{C}_d^T \boldsymbol{\lambda} = 0 \quad (27) $$

where X_d represents the vector of the generalized displacements of the diaphragm with consistent matrices M_d, D_d and K_d for mass, damping and stiffness of the system. Matrix C_d is the diaphragm contribution to the contact matrix C.

## 6 General algorithm

Coupled nonlinear equations of motion (25), (26) and (27) are now solved using the central finite difference approach. Using notation X = {X_r, X_d} storing the rotor and diaphragm generalized displacements (X reduces to X_r for diaphragm D1), and a ≡ X or φᵣ, the discretization in time yields:

$$ \dot{a} = \frac{a^{n+1} - a^{n-1}}{2\Delta t} \quad \text{and} \quad \ddot{a} = \frac{a^{n+1} - 2a^n + a^{n-1}}{\Delta t^2} \quad (28) $$

where a^n refers to quantity a at time t^n and Δt to the time step. The proposed time marching procedure simultaneously includes the calculation of the displacements and the Lagrange multipliers within three main steps [19]:

1. **Prediction:** the governing equations are solved assuming that there is no contact between the two structures.

2. **Contact detection:** the gap function g^p(x)² is computed with the predicted displacements. A correction of the displacements is required if a penetration is detected between the contacting structures.

3. **Correction:** if a penetration is detected (the gap value is negative), the corresponding Lagrange multiplier in the normal direction is calculated in order to satisfy the non-penetrability condition. The tangential contact force is then deduced from the Coulomb friction law eventually allowing for the calculation of the contact resisting torque within quasi-Newton loop.

At each time step, the matrices depending on the spinning speed and acceleration are updated. If a contact is detected, the rotational velocity and acceleration change according to the contact friction torque. Details of the time stepping procedure are provided in algorithm 1. The time step must guarantee the

<footnote>superscript p for prediction</footnote>
<caption>**Algorithm 1 Proposed solution method**</caption>

Initialization of X and $\varphi_r$ at $t_0$ and $t_1$:
$X^0, X^1$ and $\varphi_r^0, \varphi_r^1$
**for** $n = 2, \dots, n_{end}$ **do**
  prediction of $(X^{n+1}, \varphi_r^{n+1})$ assuming a linear angular acceleration and no contact
  search for $(X^{n+1}, \varphi_r^{n+1})$ satisfying equations (25), (26) and (27)
  **while** residual $> \varepsilon$ **do**
    **if** no penetration **then**
      $\lambda_N^{(k)} = 0$
    **else**
      Lagrange multipliers calculation
      $\lambda_N^{(k)} \leftarrow \lambda_N^p$ so that (21) is satisfied<fnref n="2" />
      correction of $(X^{n+1}, \varphi_r^{n+1})$ through $\lambda_N, \lambda_T$ and $C_{fric}$
    **end if**
    new solution guess $(X^{n+1}, \varphi_r^{n+1})$
  **end while**
  **if** $\varphi_r^{n+1} < \varphi_r^n$ **then**
    break<fnref n="3" />
  **end if**
  $(X^{n-1}, \varphi_r^{n-1}) \leftarrow (X^n, \varphi_r^n)$
  $(X^n, \varphi_r^n) \leftarrow (X^{n+1}, \varphi_r^{n+1})$
  $t \leftarrow t + \Delta t$
**end for**

numerical stability of the integration scheme and satisfy the Courant–Friedrichs–Lewy (CFL) criterium for an explicit technique: a linear analysis proves that the Lagrange multiplier do not modify the linear time step size. However, numerical experience suggests that contact conditions do slightly reduce the time step [3], consideration which is accounted for in the present study.

# 7 Time integration results

The developed algorithm is validated through the computation of a linear speed transient where $\varphi$ is known. Results were successfully compared with a commercial code. Consequently, calculations of transient run-down responses during accidental shutdowns are exclusively considered in what follows.

It is recalled from section 4.3 that the initial conditions in displacement and velocity of the shut-down correspond to the steady state at nominal operating conditions preceding the shutdown. From this initial state, the shutdown is initiated by setting $C_{max} = 0$ and by introducing the Newtonian fluid and aerodynamical friction resisting torques together with a heavy blade-off mass imbalance. To authors' knowledge, there is neither analytical nor numerical reference solutions. Accordingly, subsequent simulations are validated through a time step size convergence analysis. The mechanical parameters of the study are listed in table 1: they are adopted so that the first eigenfrequency of the real investigated turbine is retrieved.

## 7.1 Diaphragm D1

### 7.1.1 Without torsion

In a preliminary study, the behavior of diaphragm D1 is explored without torsional vibrations. The imbalance mass is set to 45 kg. Figure 11(a) depicts the computed rotational velocity. It illustrates the convergence with respect to the time step. Initially spinning at its nominal speed, the turbine slows down

<footnote n="2">²The Lagrange multipliers are calculated assuming a constant angular position in the correction step. This approach is valid since the solution guess takes place in the neighborhood of the converged solution.</footnote>

<footnote n="3">³This condition means that the shaft is now rotating in a counter-clockwise direction: this is physically possible but not of interest in the present study.</footnote>
<table><caption>TABLE 1: Numerical characteristics of the shaft line. The following equality for diaphragm D3 holds: radius of the inner ring = outer radius of the shaft + gap. Coupling terms of stiffness and damping matrices of the bearing are such as $k_{oxy} = k_{oyx} = d_{oxy} = d_{oyx} = 0$. Diaphragm D2 model has been constructed based on a previous EDF study and simplifying assumptions namely $k_{dxy} = k_{dyx} = d_{dxy} = d_{dyx} = 0$; conditions $m_{dxy} = m_{dyx} = 0$ are also considered. (a) rigid disk, (b) shaft, (c) diaphragm D3, (d) imbalance, (e) oil-film bearing, (f) diaphragm D2</caption><thead><tr><th></th><th>(a)</th><th>(b)</th><th>(c)</th></tr></thead><tbody><tr><td>Inner radius</td><td>0.5 m</td><td>E<sub>s</sub></td><td>E<sub>d</sub></td></tr><tr><td>Thickness</td><td>0.1 m</td><td>G<sub>s</sub></td><td>G<sub>d</sub></td></tr><tr><td>Outer radius</td><td>1 m</td><td>Mass</td><td>Mass</td></tr><tr><td>Mass density</td><td>7860 kg/m<sup>3</sup></td><td>density ρ<sub>s</sub></td><td>density ρ<sub>d</sub></td></tr><tr><td></td><td></td><td>Total</td><td>gap</td></tr><tr><td></td><td></td><td>length l<sub>s</sub></td><td>Blade length l<sub>b</sub></td></tr><tr><td></td><td></td><td>Radius</td><td>Inner ring cross section</td></tr><tr><td></td><td></td><td>k<sub>s</sub></td><td>Width</td></tr><tr><td></td><td></td><td></td><td>Height</td></tr><tr><td></td><td></td><td></td><td>Stiffener cross section</td></tr><tr><td></td><td></td><td></td><td>Width</td></tr><tr><td></td><td></td><td></td><td>Height</td></tr></tbody></table>

<table><thead><tr><th></th><th>(d)</th><th>(e)</th><th>(f)</th></tr></thead><tbody><tr><td>Mass m<sub>i</sub></td><td>45 kg</td><td>k<sub>oxx</sub> 2·10<sup>5</sup> N/m</td><td>k<sub>dxx</sub> 4·10<sup>9</sup> N/m</td></tr><tr><td>Distance to the rotor centerline r<sub>i</sub></td><td>1 m</td><td>k<sub>oyy</sub> 5·10<sup>5</sup> N/m</td><td>k<sub>dyy</sub> 4·10<sup>9</sup> N/m</td></tr><tr><td>Initial phase shift</td><td>0 rad</td><td>d<sub>oxx</sub> η × 2·10<sup>5</sup> N·s/m</td><td>d<sub>dxx</sub> 3.8·10<sup>5</sup> N·s/m</td></tr><tr><td></td><td></td><td>d<sub>oyy</sub> η × 5·10<sup>5</sup> N·s/m</td><td>d<sub>dyy</sub> 3.8·10<sup>5</sup> N·s/m</td></tr><tr><td></td><td></td><td>η 2·10<sup>-4</sup></td><td>m<sub>dxx</sub> 10<sup>4</sup> N·s<sup>2</sup>/m</td></tr><tr><td></td><td></td><td></td><td>m<sub>dyy</sub> 10<sup>4</sup> N·s<sup>2</sup>/m</td></tr></tbody></table>

due to fluids friction with no contact detected between rotating and stationary parts. When approaching the first critical speed (19.6 Hz) at $t = 3$ s, the initial 8 mm gap is completely consumed as shown in figure 12(a) and re-opens at $t = 6$ s. During interaction, the angular deceleration of the shaft is clearly affected by the contact forces.

The clearance between the rotor and the diaphragm, depicted in figure 12(a), shows that the contact is well treated: even for a large time step still ensuring convergence, the residual penetration<fnref n="3" /> is negligible (lower than 0.2 µm, value to be compared with the initial gap of 8 mm). Further results which are given in Section 7.4 show that the residual penetration decreases with the time-step: this confirms the stability and consistency properties of the proposed algorithm.

In order to analyze the frequency content of the rotor vibrations, a fast Fourier transform (FFT) is performed. As expected, the first bending modes of the rotor are excited by the imbalance: the FFT of the rotor response before contact, from $t = 0$ to $t = 3$ s, depicted in figure 15(a), shows synchronous components together with natural frequencies of the shaft line in agreement with Campbell diagram 14(a).

The sensitivity of the response to the friction coefficient is depicted in figure 16. For a small value, only the forward whirl motion of the rotor is excited during interaction. When the friction coefficient reaches a threshold ($\mu \approx 0.2$), a backward whirl phenomenon is observed and the shaft line is violently stopped. Due to this specific reason, the friction coefficient is now set to $\mu = 0.1$. Furthermore, only results for $\Delta t = 10^{-5}$ s are now presented since convergence of the quantities of interest with respect to the time step is achieved (see gap function in figure 13(a) and displacements in figure 11(b) for instance).

<footnote n="3">³The residual penetration depends on the numerical precision of the computer as well as the numerical error in inverting matrices and the precision $\epsilon$ in the nonlinear solver.</footnote>
<figure><img src="image_12.png" /><figcaption>FIGURE 11: Rotational velocity for $\dot{\phi}_r$ versus time in different configurations: (a) D1 without torsion: $\Delta t = 10^{-5}$s (- -), $\Delta t = 5 \cdot 10^{-6}$s (· - ·), $\Delta t = 10^{-6}$s (· · -), $\Delta t = 5 \cdot 10^{-7}$s (---), (b) D1 with torsion (---) and without torsion (- - -), (c) D1 (—), D2 (----), D3 (- - -) without torsion</figcaption></figure>

### 7.1.2 With torsion

When torsion $\beta_r$ is accounted for, new torsional modes arise in the modal content of the shaft as illustrated in the Campbell diagram of figure 14(b), such as $f = 147.1$ Hz for the first mode of torsion: they do not modify the flexural mode frequencies and do not depend on the rotational velocity. Torsional angle $\beta_r$ is shown figure 17. Its frequency spectrum is exhibited in figure 15(b) before interaction and in figure 15(c) during interaction, is complex to analyze. However the amplitude level are relatively low and the dominant frequency corresponds to the first torsional mode, see figure 14(b).

The angular speed, the clearance separating the shaft from the diaphragm as well as the vertical displacement are given in figure 11(b), figure 12(b) and figure 13(b) respectively. Comparison with the case-study involving no torsion shows that these results are almost identical. Similar conclusions hold for other quantities of interest such as the load on bearing, see Section 7.4. Since this paper aims at providing a numerical tool capable to predict the overall response of the turbine, in the sequel, torsion is not accounted for. This allows for simplifying the model and saving subsequent computation time while preserving a good approximation of quantities of interest.
<figure><img src="image_13.png" /><figcaption>FIGURE 12: Clearance at node 4 between the shaft and surrounding diaphragm in different configurations: (a) D1 without torsion: $\Delta t = 10^{-5}$s (—) and $\Delta t = 5 \cdot 10^{-6}$s (—), (b) D1 with torsion, (c) D2 without torsion, (d) D3 without torsion</figcaption></figure>

## 7.2 Diaphragm D2

Simulations are now performed with the more realistic diaphragm D2. The rotational velocity is depicted in figure 11(c) and compared with the one obtained with diaphragms D1 and D3⁴. Simulations indicate a softer contacting behavior of the structures due to the new mechanical flexibility. Although this diaphragm D2 is flexible, the model leads to results which are very close to those obtained with model D1 as displayed in figure 12(c) in comparison to figure 12(a) for the gap function, and figure 13(c) in comparison to figure 13(a) for the rotor vertical displacement.

Figure 15(d) depicts the spectrum of the vertical displacement of the diaphragm during interaction. The natural frequency of the diaphragm, equal to 100.7 Hz, is contained in the response. Concerning the other frequencies of the diaphragm response, these are common with the frequency content of the rotor during interaction (see figure 15(d) in comparison to figure 15(a)) and result from the coupling between rotor and diaphragm dynamics during contact.

<footnote n="4">⁴The model of diaphragm D3 is introduced in the next section</footnote>
<figure><img src="image_14.png" /></figure>

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

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

<figure><img src="image_17.png" /><figcaption>FIGURE 13: Vertical displacement of the rotor at node 2 for different configurations: (a) D1 without torsion: $\Delta t = 10^{-5}$s (---) and $\Delta t = 5 \cdot 10^{-6}$s (—), (b) D1 with torsion, (c) D2 without torsion, (d) D3 without torsion</figcaption></figure>

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

<figure><img src="image_19.png" /><figcaption>FIGURE 14: Campbell diagrams: (a) without torsion, (b) with torsion</figcaption></figure>
<figure><img src="image_20.png" /></figure>

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

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

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

<figure><img src="image_24.png" /><figcaption>FIGURE 15: Frequency spectra in different configurations: FFT performed from $t = 0$ to $t = 3.5$ s (a) rotor before contact, (b) torsion angle $\beta_r$ before interaction with D1, (c) torsion angle $\beta_r$ during interaction with D1, (d) D2 (---) and rotor (—) during interaction, (e) D3 (---) and rotor (—) during interaction</figcaption></figure>

## 7.3 Diaphragm D3

Once again, the rotational velocity evolution during shut-down is very similar to those obtained for diaphragm models D1 and D2, see figure 11(c). By contrast, the vibration levels of the shaft appear to be
<figure><img src="image_25.png" /><figcaption>FIGURE 16: Influence of friction coefficient on the computed rotational velocity: μ = 0.1 (---·), μ = 0.2 (· - ·), μ = 0.3 (---), μ = 0.5 (—), μ = 1 (—)</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>FIGURE 17: Torsion angle β<sub>r</sub> for the diaphragm D1</figcaption></figure>

significantly higher than in the previous models, as shown in figure 13(d). From the gap function evolu-
tion displayed in figure 12(d), it is noticeable that the flexible diaphragm vibrates more during a longer
contact interaction, up to 8 s. As illustrated in figures 15(d) and 15(e), the frequency content of the rotor
is similar to the one interacting with diaphragm D2. Also, the frequency content of the new diaphragm
D3 contains higher harmonics because of the chosen discretization but the magnitude of the main har-
monics is very similar. Moreover, a beating vibration phenomenon is observed in the rotor responses
during and after interaction. This phenomenon is more obvious in figure 13(d) and is peculiar to situa-
tions where comparable frequency contents of different structures are combined, through contact efforts
in the present study.

**7.4 Influence of the diaphragm model**

In a general manner, the rotor response during accidental situations depends on the geometrical and
material properties of the diaphragm. Since the load acting in the bearing is the main quantity of interest
when designing a turbine, its time evolution with respect to the diaphragm model is now investigated.
A linear calculation, thus discarding contact constraints in the formulation, of the rotor speed transient
provides the maximal possible load in bearings equal to 2 · 10<sup>6</sup> N. However, a better estimation of this
value is obtained when including the rotor-to-diaphragm contact constraints and obviously depends on the diaphragm models detailed in table 2.

Diaphragm D1 corresponds to the fastest calculation in term of CPU time but underestimates the bearing load. Also, the influence of torsion appears to be negligible and only increases the computational time. Since an enhanced version of this model, diaphragm D2 represents a good compromise between CPU time and approximation of the bearing load but leads to a bouncing behavior of the rotor because of a non physical contact treatment inherent to rigid body modelings. Finally, flexible diaphragm D3 provides the most accurate estimation of the loads exerting in the bearings. Due to the flexibility of this model, the shaft and diaphragm displacements reach larger amplitudes directly affecting the bearing loads. The numerical cost is inevitably a significant increase of the CPU time.

<table><caption>TABLE 2: Sensitivity of the results to the diaphragm models — column 1: $\Delta t$ (s), column 2: CPU time (s), column 3: maximal load on bearing (10⁵ N) and column 4: residual penetration (µm). (a) diaphragm D1 without torsion, (b) diaphragm D1 with torsion, (c) diaphragm D2 without torsion, (d) diaphragm D3 without torsion</caption><thead><tr><th></th><th colspan="4">(a)</th><th colspan="4">(b)</th></tr></thead><tbody><tr><td>10<sup>-5</sup></td><td>2,007</td><td>3.5424</td><td>0.76720</td><td>10<sup>-5</sup></td><td>2,738</td><td>3.4427</td><td>0.69885</td></tr><tr><td>5 &middot; 10<sup>-6</sup></td><td>4,389</td><td>3.5490</td><td>0.22468</td><td>5 &middot; 10<sup>-6</sup></td><td>5,456</td><td>3.4632</td><td>0.24114</td></tr><tr><td>10<sup>-6</sup></td><td>24,292</td><td>3.5309</td><td>0.01492</td><td>10<sup>-6</sup></td><td>26,684</td><td>3.4244</td><td>0.00807</td></tr><tr><td>5 &middot; 10<sup>-7</sup></td><td>47,980</td><td>3.5296</td><td>0.01636</td><td>5 &middot; 10<sup>-7</sup></td><td>53,111</td><td>3.5539</td><td>0.00205</td></tr></tbody></table><table><thead><tr><th></th><th colspan="4">(c)</th><th colspan="4">(d)</th></tr></thead><tbody><tr><td>10<sup>-5</sup></td><td>2,012</td><td>3.7041</td><td>0.74075</td><td>10<sup>-5</sup></td><td>7,742</td><td>5.4562</td><td>58.8643</td></tr><tr><td>5 &middot; 10<sup>-6</sup></td><td>4,431</td><td>3.7437</td><td>0.14734</td><td>5 &middot; 10<sup>-6</sup></td><td>15,918</td><td>5.4513</td><td>15.3205</td></tr><tr><td>10<sup>-6</sup></td><td>22,458</td><td>3.7343</td><td>0.00836</td><td>10<sup>-6</sup></td><td>80,629</td><td>5.4522</td><td>0.62922</td></tr><tr><td>5 &middot; 10<sup>-7</sup></td><td>44,186</td><td>3.7295</td><td>0.00451</td><td>5 &middot; 10<sup>-7</sup></td><td>161,265</td><td>5.4547</td><td>0.15778</td></tr></tbody></table>

# 8 Conclusions and prospects

The emphasis of the present study was twofold.

The purpose was first to establish a new theoretical formulation for a speed transient analysis of a rotor supported by linearized journal bearings together with gyroscopic effects, torsional and flexural couplings and contact/rubbing dynamics. The rotational velocity is not known and has to be simultaneously solved with the equations of motion. Discretizations in time and space are performed with the explicit and well-known central finite difference scheme and the usual finite element approach, respectively. Contact constraints are treated through the Lagrange multiplier technique in order to avoid unwanted residual penetrations between contacting structural components. Several diaphragm models were explored, starting from a fixed rigid one, up to a fully flexible one which incorporates the main features of the true structure. The related node-to-line contact algorithm has shown to be reliable since convergence with respect to the time step is achieved. Accordingly, a tool for analyzing rotor-diaphragm interactions is now available for the study of the turbine behavior, especially as it passes through critical speeds after a blade-off.

Second, the developed numerical tool is capable of computing global quantities, such as bearing loads which are of great interest in the design process of a turbine, as well as usual structural displacements. Main results are: (1) torsional twist only slightly influences the maximum load acting in the bearings but allows for a reliable estimation of the torsional stresses along the shaft. As a matter of fact, bearing loads appear to be more dependent on the diaphragm modeling. (2) Rigid diaphragms overestimate contact loads with the shaft and consequently underestimate loads in the bearings. (3) Displacements of flexible diaphragms due to rubbing are larger than in the case of rigid models. This implies an augmentation of the shaft vibrations and of the bearing loads. The use of flexible diaphragm is therefore mandatory for an accurate prediction of the bearing loads.
References

[1] L. Baillet, V. Linck, S. D'Errico, B. Laulagnet, and Y. Berthier. Finite element simulation of dynamic instabilities in frictional sliding contact. *Journal of Tribology*, 127(3):652–657, 2005.
doi:10.1115/1.1866160.

[2] A.R. Bartha. *Dry friction backward whirl of rotors*. PhD thesis, Swiss Federal Institute of Technology of Zurich, 2000.

[3] T. Belytschko and M. O. Neal. Contact-impact by the pinball algorithm with penalty and lagrangian methods. *International Journal for Numerical Methods in Engineering*, 31:547–572, 1991.
doi:10.1002/nme.1620310309.

[4] T. Belytschko, W. K. Liu, and B. Moran. *Nonlinear Finite Elements for Continua and Structures*. John Wiley and Sons, New York, 2000.

[5] L.-W. Chen and D.-M. Ku. Dynamic stability of a cantilever shaft-disk system. *Journal of Vibration and Acoustics*, 114:326–329, 1992.

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

[7] F. Chu and Z. Zhang. Periodic quasi-periodic and chaotic vibration of a rub impact rotor system supported on oil film bearing. *International Journal of Engineering Science*, 35(10–11):963–973, 1997.
doi:doi:10.1016/S0020-7225(97)89393-7.

[8] F. Chu and Z. Zhang. Bifurcation and chaos in a rub-impact jeffcott rotor system. *Journal of Sound and Vibration*, 210:1–18, 1998.
doi:doi:10.1006/jsvi.1997.1283.

[9] J.-C. Clément. Projet macadam — influence de la prise en compte du contact rotor-stator sur l'estimation des efforts aux paliers des groupes turbo-alternateurs de 900mw et 1300 mw en situation accidentelle (macadam project — influence of the rotor-to-stator contact forces on the bearing loads of 900mw and 1300mw turbogenerators in accidental situations). Technical report, 2004.

[10] X. Dai, X. Zhang, and X. Jin. The partial and full rubbing of a flywheel rotor-bearing-stop system. *International Journal of Mechanical Sciences*, 43:505–519, 2001.
doi:10.1016/S0020-7403(99)00121-6.

[11] X. Dai, Z. Jin, and X. Zhang. Dynamic behavior of the full rotor/stop rubbing : numerical simulation and experimental verification. *Journal of Sound and Vibration*, 251:807–822, 2002.
doi:doi:10.1006/jsvi.2001.3998.

[12] S. Dubigeon. *Mécanique des milieux continus (Continuum mechanics)*. Lavoisier, Tec & Doc, Nantes, 1998. ISBN 978-2-912-98500-2.

[13] T. Fortin. Modélisation du ralentissement d'un groupe turbo-alternateur (modelling of a turbogenerator slowdown). Technical report, 1993.

[14] I. Giraudon-Guilloteau. *Modélisation du contact en dynamique: construction d'un élément simplifié de contact et application à l'interaction rotor/stator (Contact modeling in structural dynamics: construction of a simplified contact finite element and application to rotor/stator interaction)*. PhD thesis, École Centrale de Nantes, 1999.

[15] M. Lalanne and G. Ferraris. *Rotordynamics Prediction in Engineering*. John Wiley, second edition, 1998. ISBN 978-0-471-97288-4.
[16] T. A. Laursen. *Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis*. Springer, Berlin, 2003. ISBN 978-3-540-42906-7.

[17] C. W. Lee. *Vibration Analysis of Rotors*. Springer, 1993. ISBN 978-0-792-32300-6.

[18] M. Legrand. *Modèles de prédiction de l'interaction rotor/stator dans un moteur d'avion (Models of rotor/stator interaction in aircraft engines)*. PhD thesis, École Centrale de Nantes, 2005. URL http://tel.archives-ouvertes.fr/tel-00011631/fr/.

[19] 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:doi:10.1016/j.jsv.2008.06.019. URL http://hal.archives-ouvertes.fr/hal-00328186/fr/.

[20] N. Lesaffre, J.-J. Sinou, and F. Thouverez. Model and stability analysis of a flexible bladed rotor. *International Journal of Rotating Machinery*, 2006, 2006. doi:10.1155/IJRM/2006/63756.

[21] A. Muszynska. Synchronous and self-excited rotor vibrations caused by a full annular rub. In *Machinery Dynamics—Eightth Seminar*, pages 22.1–22.21, Halifax, Nova Scotia, Canada, 1984. Canadian Machinery Association.

[22] A. Muszynska. Rotor-to-stationary part full annular contact modelling. In *Proceedings of The Ninth International Symposium on Transport Phenomena and Dynamics of Rotating Machinery*, Honolulu, Hawaii, February 2002.

[23] A. Muszynska. *Rotordynamics*. CRC Taylor & Francis Group, New-York, 2005. ISBN 0-8247-2399-6.

[24] A. Muszynska and P. Goldman. Chaotic responses of unbalanced rotor/bearing/stator systems with looseness or rubs. *Choas, Solitons & Fractals*, 5:1683–1704, 1995. doi:doi:10.1016/0960-0779(94)00171-L.

[25] A. Muszynska, D. E. Bently, W. D. Franklin, R. D. Hayashida, L. M. Kingsley, and A. E. Curry. Influence of rubbing on rotor dynamics – part i & ii. Technical report, NASA Contract N\'r NAS8-36179, 1989.

[26] L. Noels. *Contributions aux algorithmes d'intégration temporelle conservant l'énergie en dynamique non-linéaire des structures (Contributions to energy-conserving time integration algorithms for nonlinear dynamics)*. PhD thesis, Université de Liège, 2004.

[27] L. Noels, L. Stainier, and J.-P. Ponthot. Combined implicit/explicit algorithms for crashworthiness analysis. *International Journal of Impact Engineering*, 30:1161–1177, 2004. doi:doi:10.1016/j.ijimpeng.2004.03.004.

[28] E. E. Pavlovskaia, E. V. Karpenko, and M. Wiercigroch. Nonlinear dynamic interactions of a jeffcott rotor with preloaded snubber ring. *Journal of Sound and Vibration*, 276(1-2):361–379, 2004. ISSN 0022-460X. doi:10.1016/j.jsv.2003.07.033.

[29] H.C. Piccoli and H.I. Weber. Experimental observation of chaotic motion in a rotor with rubbing. *Nonlinear Dynamics*, 16:55–70, 1998. doi:10.1023/A:1008284317724.

[30] S. Roques, M. Legrand, C. Stoisser, P. Cartraud, and C. Pierre. Development of beam-to-beam contact detection algorithms for rotor-stator rubbing applications. In *Proceedings of IDETC/CIE 2007 — 21st ASME Biennial Conference on Mechanical Vibration and Noise*, Las Vegas, Nevada, September 2007.

[31] J. C. Simo and T. A. Laursen. An augmented lagrangian treatment of contact problems involving friction. *Computers & Structures*, 42(1):97–116, 1992. doi:10.1016/0045-7949(92)90540-G.
[32] R. Taylor, N. Carpenter, and M. Katona. Lagrange constraints for transient finite element surface contact. *International Journal for Numerical Methods in Engineering*, 32:103–128, 1991.
doi:10.1002/nme.1620320107.

[33] P. Wriggers. *Computational Contact Mechanics*. Wiley, 2002.