# Modelling of rotor speed transient with rotor to stator contact

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

► To cite this version:

Sébastien Roques, Carlo Stoisser, Patrice Cartraud, Mathias Legrand, Bernard Peseux, et al. Modelling of rotor speed transient with rotor to stator contact. 7th IFToMM International Conference on Rotor Dynamics, Sep 2006, Vienne, Austria. hal-01008503

HAL Id: hal-01008503

https://hal.science/hal-01008503

Submitted on 17 Aug 2016

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
Modelling of rotor speed transient with rotor-to stator contact

Sébastien Roques

GeM, Pôle Structure et Couplage
École Centrale de Nantes
1, rue de la Noë
44321 Nantes cedex 3, France
sebastien.roques@ec-nantes.fr

Carlo M. Stoisser

EDF Research & Development
1, avenue du Général de Gaulle
92140 Clamart, France
carlo.stoisser@edf.fr

Patrice Cartraud

GeM, Pôle Structure et Couplage
École Centrale de Nantes
1, rue de la Noë
44321 Nantes cedex 3, France
patrice.cartraud@ec-nantes.fr

Mathias Legrand

GeM, Pôle Structure et Couplage
École Centrale de Nantes
1, rue de la Noë
44321 Nantes cedex 3, France
mathias.legrand@ec-nantes.fr

Bernard Peseux

GeM, Pôle Structure et Couplage
École Centrale de Nantes
1, rue de la Noë
44321 Nantes cedex 3, France
bernard.peseux@ec-nantes.fr

Christophe Pierre

McGill University
Pavillon McDonald
817, rue Sherbrooke ouest
Montréal H3A 2K6, Canada
christophe.pierre@mcgill.ca

ABSTRACT

A rotor-stator model of a turbogenerator has been developed to solve rotational velocity transients: to this end, the angular position of a beam cross section is considered as a variable, allowing for the angular deceleration due to rotor-to-stator friction to be calculated in a more realistic representation of the interaction. The equations of motion derived from an energetic formulation are highly nonlinear and require a specific time integration procedure. Calculations were run considering different types of diaphragms.

KEY WORDS

Speed transient, heavy rotor-to-stator contact, rubbing, Lagrange constraints.

NOMENCLATURE

<table><tr><td>w</td><td>traction (Z axis)</td><td>β</td><td>torsion (Z direction)</td></tr><tr><td>(u, θ<sub>v</sub>)</td><td>bending in (X,Z) plane</td><td>(v, θ<sub>u</sub>)</td><td>bending in (Y,Z) plane</td></tr><tr><td>φ</td><td>angular position</td><td>s</td><td>curvilinear abscissa</td></tr><tr><td>(E<sub>c</sub>, E<sub>d</sub>)</td><td>respectively kinetic and strain energies</td><td>(k<sub>ij</sub>, d<sub>ij</sub>)</td><td>stiffness and damping coefficients</td></tr><tr><td>φ̇</td><td>angular spinning velocity</td><td>Ω<sub>nom</sub></td><td>nominal rotating speed</td></tr><tr><td colspan="4"><em>Disk nomenclature</em></td></tr><tr><td>M<sub>D</sub></td><td>disk mass</td><td>I<sub>DX</sub></td><td>moment of inertia</td></tr><tr><td>I<sub>DZ</sub></td><td>polar moment of inertia</td><td></td><td></td></tr><tr><td colspan="4"><em>Beam nomenclature</em></td></tr><tr><td>E</td><td>Young's modulus</td><td>ν</td><td>Poisson's ratio</td></tr><tr><td>G</td><td>shear modulus</td><td>S</td><td>cross section area</td></tr><tr><td>k<sub>x</sub></td><td>transverse shear form factor</td><td>λ = <sup>12EI<sub>x</sub></sup>&frasl;<sub>k<sub>x</sub>GSl<sup>2</sup></sub></td><td>transverse shear effect</td></tr><tr><td>I<sub>x</sub></td><td>moment of inertia</td><td>I<sub>p</sub></td><td>polar moment of inertia</td></tr><tr><td>ρ</td><td>mass density</td><td>l</td><td>element length</td></tr><tr><td>z</td><td>axial distance</td><td></td><td></td></tr><tr><td colspan="4"><em>Subscripts</em></td></tr><tr><td>D</td><td>disk</td><td>b</td><td>imbalance</td></tr><tr><td>c</td><td>casing</td><td></td><td></td></tr></table>
<figure><img src="image_2.png" /><figcaption>Figure 1: Nomenclature</figcaption></figure>

# 1 INTRODUCTION

In nuclear power plant turbosets, the design basis accident consists of the loss of a blade on the low pressure turbine last stage. During the accidental shut down, a severe rotor-casing interaction may occur at critical speeds due to large shaft line displacements originated by a high unbalance excitation. The contact between the shaft and the diaphragm induces an important deceleration rate and modifies the turbogenerator dynamics including the amplitude of the loads in the bearings. Risks of failure of the contacting components may also occur due to the heavy friction torque. Therefore the main objective is to verify that the designed turbine is capable of going through critical speeds without catastrophic consequences for the shaft line.

In a general manner, contact dynamic is highly nonlinear and many studies have been carried out to understand the rubbing phenomenon in order to prevent serious damages. Most of papers usually make the assumption that the angular velocity is known, even during contact [1]. As far as we know, a few models have considered the angular position as a variable and the rotor was assumed to be rigid ([2] and [3]).

In this paper, only the shaft-diaphragm interaction is considered: the blades-to-diaphragm contact [4] is assumed to be negligible in comparison with the shaft-stator interaction to compute the angular deceleration due to friction. Moreover the nonlinear behavior of the bearings is linearized. The Lagrange multipliers method in conjunction with an explicit time-stepping procedure is used to solve equations of motion [5].

# 2 ROTOR SPEED TRANSIENT MODELLING

## 2.1 Shaft line components

A 1-D model has been developed considering the angular position of a cross section as a variable. Flexural vibrations in two orthogonal planes, torsional vibrations of a straight shaft and axial vibrations are considered. The following assumptions are used throughout the paper:

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

2. Timoshenko theory is used.

3. The shaft is balanced (no mass eccentricity).

4. Gyroscopic terms are included.

5. External torques with constant direction along the undeformed centroidal axis are applied at each extremity of the shaft.

### 2.1.1 Rigid disk

Considering a rigid circular disk with an unknown angular speed, and due to its geometric symmetry, the kinetic energy $E_c^D$ takes the form:

$$2E_c^D = M_D[\dot{u}^2 + \dot{v}^2 + \dot{w}^2] + I_{DX}[\dot{\theta}_u^2 \cos^2 \theta_v + \dot{\theta}_v^2] + I_{DZ}[\dot{\theta}_u^2 \sin^2 \theta_v + (\dot{\varphi} + \dot{\beta})^2 + 2(\dot{\varphi} + \dot{\beta})\dot{\theta}_u \sin \theta_v] \quad (1)$$
where $M_D, I_{DX}, I_{DZ}$ respectively stand for the mass, the second moment of inertia and the polar moment of inertia of the disk. Assuming small rotations and small displacements, the kinetic energy is rewritten as follows:

$$2E_c^D = M_D [\dot{u}^2 + \dot{v}^2 + \dot{w}^2] + I_{DX} [\dot{\theta}_u^2 + \dot{\theta}_v^2] + I_{DZ} [\dot{\phi} + \dot{\beta}]^2 + 2 (\dot{\phi} + \dot{\beta}) \dot{\theta}_u \dot{\theta}_v \quad (2)$$

This formula is exact up to the third order and equations of motion will be written up to the second order.

### 2.1.2 Shaft

The kinetic energy $E_c$ of a rotating beam is written by considering that the energy of the shaft corresponds to the integration of the disk energy along the longitudinal direction:

$$2E_c = \rho \int_0^l \left[ S (\dot{u}^2 + \dot{v}^2 + \dot{w}^2) + I_x (\dot{\theta}_u^2 + \dot{\theta}_v^2) + I_p \left( (\dot{\phi} + \dot{\beta})^2 + 2 (\dot{\phi} + \dot{\beta}) \dot{\theta}_u \dot{\theta}_v \right) \right] dz \quad (3)$$

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

$$2E_d = \int_0^l \left\{ E S w_{,z}^2 + E I_x (\theta_{u,z}^2 + \theta_{v,z}^2) + k_x G S \left[ (u,z - \theta_v)^2 + (v,z + \theta_u)^2 \right] \right\} dz \quad (4)$$

Bending and torsion are coupled together because of gyroscopic terms and to the unknown angular position. The equations of motion derived from Hamilton's principle are nonlinear even without contact.

### 2.1.3 Linearized bearing

The dynamical study of a shaft line requires to describe the nonlinear behavior of the oil film. Assuming small displacements of the rotor, stiffness and damping coefficients of the oil film can be calculated by linearizing Reynolds' equations with respect to the equilibrium position. Calculating the virtual work $\delta W$ of external forces acting on the shaft yields:

$$\delta W = [F_u \ F_v] \cdot \begin{Bmatrix} \delta u \\ \delta v \end{Bmatrix}$$

where $F_u$ and $F_v$ are the components of the generalized forces acting in the bearings. After linearizing, they can be written in a matrix form:

$$\begin{Bmatrix} F_u \\ F_v \end{Bmatrix} = - \begin{bmatrix} k_{xx} & k_{xy} \\ k_{yx} & k_{yy} \end{bmatrix} \cdot \begin{Bmatrix} u \\ v \end{Bmatrix} - \begin{bmatrix} d_{xx} & d_{xy} \\ d_{yx} & d_{yy} \end{bmatrix} \cdot \begin{Bmatrix} \dot{u} \\ \dot{v} \end{Bmatrix}$$

### 2.1.4 Imbalance

The imbalance is a concentrated mass. Then writing the kinetic energy, it comes:

$$2E_c^B = m_b [\dot{u}^2 + \dot{v}^2 + \dot{w}^2 + r_b^2 \dot{\phi}^2 + 2r_b \dot{\phi} (\dot{v} \cos \varphi - \dot{u} \sin \varphi)] \quad (5)$$

### 2.1.5 Diaphragm

Different types of casing have been considered.

**Casing C1: rigid and fixed ring** Casing C1 (see figure 2(a)) is completely rigid and can be seen as a mathematical boundary: it was first developed to validate the contact algorithm.
<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic of diaphragm models</figcaption></figure>

**Casing C2: rigid ring with flexibility and damping** This model is an extension of the first model and based on a previous EDF study [6]: stiffnesses and dampers have been added (see figure 2(b)).

**Casing C3: flexible diaphragm** A planar flexible diaphragm (see figure 2(c)) has been developed to enrich the modelling: the inner ring is discretized in curved beam FE with four dof per node ($u_c, u_{c,s}, v_c$ and $v_{c,s}$ where s denotes path variable). The blades are discretized in straight Euler-Bernoulli beams FE with three dof per node.

## 2.2 FE modelling

**Shaft line** Traction and torsion are discretized using the usual shape functions. Writing $\xi = \frac{z}{l}$ and $q = w, \beta$ yields:

$$q(z) = [1 - \xi \xi] \begin{Bmatrix} q_1 \\ q_2 \end{Bmatrix}$$

The modified Hermite shape functions for bending used in [7] are chosen. The degrees of freedom $u(z), \theta_v(z), v(z)$ and $\theta_u(z)$ are written as follow:

$$\left\{ \begin{aligned} u(z) &= [N_{11}(z) \ N_{12}(z) \ N_{13}(z) \ N_{14}(z)] \begin{Bmatrix} u_1 \\ \theta_{v1} \\ u_2 \\ \theta_{v2} \end{Bmatrix} & \text{and} & v(z) &= [N_{11}(z) - N_{12}(z) \ N_{13}(z) - N_{14}(z)] \begin{Bmatrix} v_1 \\ \theta_{u1} \\ v_2 \\ \theta_{u2} \end{Bmatrix} \\ \theta_v(z) &= [N_{21}(z) \ N_{22}(z) \ N_{23}(z) \ N_{24}(z)] \begin{Bmatrix} u_1 \\ \theta_{v1} \\ u_2 \\ \theta_{v2} \end{Bmatrix} & \theta_u(z) &= [N_{21}(z) - N_{22}(z) \ N_{23}(z) - N_{24}(z)] \begin{Bmatrix} v_1 \\ \theta_{u1} \\ v_2 \\ \theta_{u2} \end{Bmatrix} \end{aligned} \right.$$

**Diaphragm C3** The inner ring of casing C3 is discretized with cubic polynomials in $u$ and $v$. Defining $\zeta = \frac{s}{l_c}$ the shape functions are:

$$N_1(s) = 1 - 3\zeta^2 + 2\zeta^3 \qquad N_2(s) = l_c\zeta(1 - 2\zeta + \zeta^2)$$

$$N_3(s) = \zeta(\lambda + 3\zeta - 2\zeta^2) \qquad N_4(s) = l_c\zeta(-\zeta + \zeta^2)$$

where $l_c$ denotes the length of a finite element. The discretized displacement field becomes for an element whose nodes are denoted 1 and 2:

$$\left\{ \begin{aligned} u_c(s) &= N_1(s)u_1 + N_2(s)u_{1,s} + N_3(s)u_2 + N_4(s)u_{2,s} \\ v_c(s) &= N_1(s)v_1 + N_2(s)v_{1,s} + N_3(s)v_2 + N_4(s)v_{2,s} \end{aligned} \right.$$

Equations of motion are then derived using Hamilton's principle and the stiffness, gyroscopic and imbalance matrices depend on the angular speed and acceleration.
## 2.3 Contact forces

The forces of particular interest in this study are the contact forces acting between the shaft and the diaphragm where friction is considered. As a first approach, contact is assumed to be punctual at the hot spot, i.e. highest value of rotor eccentricity (point C, see figure 3). With the assumption of small rotations and displacements, the radial gap is equal to the sum of the eccentricity of the cross section's center (distance OR) and the radius of the shaft.

<figure><img src="image_4.png" /><figcaption>**Figure 3:** Schematic of contact detection</figcaption></figure>

The radial contact force component $F_N$ is obtained by using Lagrange multipliers method, which guarantees the impenetrability condition. The tangential contact force $F_T$ is deduced using the Coulomb friction law $F_T = \mu F_N$. Finally the friction torque takes the form $C_{fric} = \mu RF_N$.

## 2.4 Boundary conditions and external loads

As the angular speed is a variable, boundary conditions are necessary for the angular speed law:

- Driven torque: assumed to be constant and equal to $C_{max}$;

- Resisting torque (alternator): equal to $-C_{max} \frac{\dot{\varphi}}{\Omega_{nom}}$ (model use to obtain an exponential shape).

During the accidental shut down, forces acting on the turbogenerator correspond to aerodynamical and fluids forces:

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

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

$A_{newt}$, $A_{aero}$ are positive constants. These friction coefficients have been identified by a EDF's study [8]. The accidental shut down is only of interest for our studies: as the shaft line is disconnected from the rest of the network there is no driven torque ($C_{max} = 0$).

## 2.5 General algorithm

A previous study [4] shows that the penalty method is not adapted for the contact detection. A sensitivity study [9] on implicit time scheme with contact constraints shows a dependency on Newmark parameters. It turns out that the central finite differences based on Carpenter's work [5] (prediction-correction algorithm) exactly satisfies the contact detection and ensures the compatibility of the speed and the acceleration. Consequently, equations of motion are discretized in time according to central finite differences. Because of the nonlinear system of equations, the following time procedure is developed as follows:
**Initialization** of variables X and $\varphi$ for $t_0$ and $t_1$: $X_0, X_1$ and $\varphi_0, \varphi_1$ are estimated

For $n$ from $n_2$ to $n_{end}$ do

**Prediction** of $(X^{n+1}, \varphi^{n+1})$ assuming that the acceleration is locally linear.

**Search of the solution** $(X^{n+1}, \varphi^{n+1})$ **satisfying the governing equations** (nonlinear solver)

While ($\|\text{Residual vector of the discretized equations of motion}\|^2 \le \varepsilon$) do

If (No penetration) then

$$| \lambda_N = 0$$

else

**Calculation of the Lagrange multipliers** $\lambda_N = \lambda_N(\text{predicted solution})^{**}$

**Correction** by radial and tangential (according to Coulomb law) contact forces and then the friction torque $C_{fric} = \mu R \lambda_N$

end If

New solution guess $(X^{n+1}, \varphi^{n+1})$

done

If $(\varphi_{n+1} \mid \varphi_n)$ then

| Computation stopped*

end If

**Results are recorded**

Next time step increment $a^{n-1} \leftarrow a^n$ and $a^n \leftarrow a^{n+1}$

end For

*: when the angular speed become negative, which means the shaft is counter rotating, the calculation is stopped because this case is not of interest.

**: the Lagrange multipliers is calculated assuming that the angular position, which is estimated in the neighborhood of the converged solution, is correct.

At each time step, the rotor displacements are predicted. The matrices depending on the spinning speed and acceleration are updated. If a contact is detected, the rotating speed and acceleration change according to the friction torque.

## 3 SIMPLIFIED MODEL OF SHAFT LINE

In order to check the developed algorithm, a simplified turbine model has been considered (see figure 4(a)). FE models for the turbogenerator and the diaphragm are chosen so that the first eigenfrequencies match the corresponding ones of the real structures. Campbell diagram of the rotor is shown in figure 4(b).

<figure><img src="image_5.png" /><figcaption>**Figure 4:** Model of shaft line</figcaption></figure>
<table><caption>**Table 1:** Numerical characteristics of the shaft line</caption><tbody><tr><td colspan="2">Disk</td><td></td><td></td></tr><tr><td>Inner radius</td><td>0.5m</td><td>Thickness</td><td>0.1m</td></tr><tr><td>Outer radius</td><td>1m</td><td>Mass density</td><td>7860kg/m³</td></tr><tr><td colspan="2">Shaft</td><td></td><td></td></tr><tr><td>E</td><td>200GPa</td><td>ν</td><td>0.3</td></tr><tr><td>Total length</td><td>10 m</td><td>Mass density</td><td>7860kg/m³</td></tr><tr><td>Radius</td><td>0.5m</td><td>k<sub>x</sub></td><td>6/7</td></tr><tr><td colspan="2">Imbalance</td><td></td><td></td></tr><tr><td>Imbalance mass</td><td>5kg for C1 and 45kg for C2 and C3</td><td>Phase shift</td><td>0rad</td></tr><tr><td>Distance between imbalance and center of the disk</td><td></td><td></td><td>1m</td></tr><tr><td colspan="2">Bearing</td><td></td><td></td></tr><tr><td>β</td><td>0.0002</td><td>(0.02% of the stiffness coefficients )</td><td></td></tr><tr><td>k<sub>xx</sub></td><td>2 · 10<sup>5</sup>N/m</td><td>k<sub>xy</sub></td><td>0N/m</td></tr><tr><td>k<sub>yx</sub></td><td>0N/m</td><td>k<sub>yy</sub></td><td>5 · 10<sup>5</sup>N/m</td></tr><tr><td>c<sub>xx</sub></td><td>β × 2 · 10<sup>5</sup>N.s/m</td><td>c<sub>xy</sub></td><td>0N.s/m</td></tr><tr><td>c<sub>yx</sub></td><td>0N.s/m</td><td>c<sub>yy</sub></td><td>β × 5 · 10<sup>5</sup>N.s/m</td></tr><tr><td colspan="2">Casing C3</td><td></td><td></td></tr><tr><td>E</td><td>200GPa</td><td>ν</td><td>0.3</td></tr><tr><td>Mass density</td><td>7860kg/m³</td><td></td><td></td></tr><tr><td colspan="4">Radius of the inner ring = outer radius of the shaft + gap (gap=8mm)</td></tr></tbody></table>

# 4 TIME INTEGRATION RESULTS

The developed algorithm was first checked and validated simulating a speed transient. The speed law was given and results have been compared to the ones obtained from a commercial code.

## 4.1 Simulations with casing C1

**Case without torsion** As a preliminary study case, the mathematical diaphragm C1 is considered. The imbalance mass is set to 5 kg in order to excite the rotor. Figure 5(a) depicts the computed spinning velocity. A zoom of this figure during the contact phase is depicted in Figure 5(b)) and the angular acceleration is plotted in Figure 5(c). According to simulations, convergence of the displacements is achieved.

<figure><img src="image_6.png" /><figcaption>**Figure 5:** Rotating speed and angular deceleration with casing C1</figcaption></figure>

The eccentricity of the node in front of the diaphragm (see figure 6(a)) shows that the contact is well considered: even for a high value of time step that ensures convergence, the residual penetration is negligible (lower than 0.2μm, value to be compared with the gap 8mm). The rotor response spectrum (see figure 6(b)) shows components at the rotating frequency and at the natural frequencies of the model.

The influence of the friction coefficient is shown in figures 7(a) and 7(b). For a small value, just the forward whirl motion is excited. With μ greater than 0.5, the reverse whirl phenomenon is observed. Backward whirl phenomena are pointed out with a lower friction coefficient with a higher imbalance mass too. These two figures show the trivial dependency on the friction coefficient which is directly linked with the contact friction torque.
<figure><img src="image_7.png" /><figcaption>Figure 6: Eccentricity of the node 4 and frequency content of the shaft line</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Influence of the friction coefficient $\mu$</figcaption></figure>

**Case with torsion** When torsion $\beta$ is added to the model, results become more complex. Figure 8(a) shows that there is convergence except during contact: during interaction, results locally depend on the time step, which may be due to high frequency phenomena. The spatial discretization of the shaft might be too poor and could explain it. Even with a high value of time step that ensures convergence, the residual penetration is negligible (lower than $0.2\mu m$, value to be compared with the gap 8mm). Further analysis must be performed to understand the results. Figure 8(c) depicts the angle of torsion $\beta$. The frequency spectrum of the latter is complex to analyze due to the nonlinear coupling between bending and torsion and the contact nonlinearity.

<figure><img src="image_9.png" /><figcaption>Figure 8: Rotating speed, deceleration and the torsion d.o.f. with casing C1</figcaption></figure>

As a intermediate conclusion, the developed procedure has been validated simulating rotor-stator interactions. It has been shown also that additional analysis concerning the d.o.f. of torsion have to be done. However flexibility and damping of the diaphragm are not included in such a casing. The unbalance mass is now set to 45kg and only simulations without considering the torsion d.o.f. are considered.
## 4.2 Simulations with casing C2

Simulations are performed with the diaphragm C2, which is more realistic but simplified. Analogously figure 9(a) depicts the calculated rotating speed and figure 9(b) the corresponding angular deceleration. The results can be compared to computations with casing C1 (see — on figures 9(a) and 9(c)): without flexibility (casing C1), the contact is harder. The contact detection is still well considered as figure 9(c) shows. Figure 10(a) depicts the spectrum of the

<figure><img src="image_10.png" /><figcaption>**Figure 9:** Rotating speed, deceleration and eccentricity with casing C2</figcaption></figure>

diaphragm displacements during interaction. Some frequencies of the stator are common with the frequency content of the rotor: it seems that these frequencies are harmonics of the rotating speed and related to the number of impacts per revolution.

<figure><img src="image_11.png" /><figcaption>**Figure 10:** Frequency spectrum when contact occurs</figcaption></figure>

## 4.3 Simulations with casing C3

Time integration results are shown figures 11(a) and 11(b). Conclusions remain the same: no dependency regarding the time step, negligible residual penetration.

Loads on bearings directly depend on the casing stiffness. If coupled with a nonlinear solver Reynolds equations solver, this algorithm could compute the real effects of nonlinear bearings and verify that the turbine passes critical speeds without catastrophic consequences.

Authors start to study a 3-D or pseudo 2-D casing FE model in order to add axial deformations of the diaphragm.

# 5 CONCLUSION AND PROSPECTS

A numerical tool has been developed to analyze speed transients of a turbogenerator with rotor-stator interactions. The developed contact algorithm have been checked according to the first time integration results. With the use of Lagrange multipliers, the computations don't depend on time step and the contact detection is well considered so that there is just numerical residual penetration.
<figure><img src="image_12.png" /><figcaption>**Figure 11:** Rotating speed and eccentricity of node 4 with casing C3</figcaption></figure>

## References

[1] Muszynska, A. (2005): *Rotordynamics*. Taylor & Francis edition.

[2] Dai, X. Jin, Z. and Zhang, X. (2002): Dynamic behavior of the full rotor/stop rubbing: numerical simulation and experimental verification. *Journal of Sound and Vibration*, 251: 807-822.

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

[4] Legrand, M. (2005): *Modèles de prédiction de l'interaction rotor/stator dans un moteur d'avion*. PhD thesis, École Centrale de Nantes.

[5] Carpenter, N., Taylor, R. and Katona, M. (1991): Lagrange constraints for transient finite element surface contact. *Int. Journal for Numerical Methods in Engineering*, 32: 103-128.

[6] Clement, J.-C. (2004): Projet Macadam – influence de la prise en compte du contact rotor-stator sur l’estimation des efforts aux paliers des GTA CP2 et 1300 MW en situation accidentelle. *Internal report EDF R&D – HT-65/03/038/A*.

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

[8] Fortin, T. (1993): Modélisation de ralentissement de GTA. *Internal report EDF R&D – HP-65/93/189*.

[9] Giraudon-Guilloteau, I. (1999) *Modélisation du contact en dynamique: Construction d'un élément simplifié de contact et application à l'interaction rotor/stator*. PhD thesis, École Centrale de Nantes.