# An analysis of rotor-stator interaction

Celio Duran, Lionel Manin, Marie-Ange Andrianoely, Céline Bordegaray, Régis Dufour

► To cite this version:

Celio Duran, Lionel Manin, Marie-Ange Andrianoely, Céline Bordegaray, Régis Dufour. An analysis of rotor-stator interaction. 10th International Conference on Vibrations in Rotating Machinery (VIRM10), IMechE, London, Sept. 2012, Sep 2012, London, United Kingdom. pp.623-633, 10.1533/9780857094537.9.623. hal-05138638

HAL Id: hal-05138638

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

Submitted on 1 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
# An analysis of rotor-stator interaction

C Duran¹,², L Manin¹, M-A Andrianoely¹, C Bordegaray², R Dufour¹

¹ Université de Lyon, CNRS, INSA-Lyon, LaMCoS UMR5259, France

² SAFRAN Turboméca, Integration Department, France

## ABSTRACT

The aim of this paper is to investigate the influence of rotor-stator interaction on the dynamic behaviour of a rotor-bearing system. First the model of a SDOF system with two stop-ends is established. The analytical solution permits testing the stability and the robustness of the penalty contact law combined with two numerical schemes. Second, a three-disk rotor in bending equipped with two bearings is subject to a disk-stator contact. The two numerical algorithms analytically validated are again used and associated to four contact models: three types of smooth penalty methods and the Lagrange multiplier method. Spectrograms of time history responses are presented and permits analysing the spectral evolution of the rotor behaviour.

## 1. INTRODUCTION

In the development of the means of transport and of energy, rotating machines hold a central place. Their technologies have to be innovative in order to increase performance and to reduce costs and noise while continuing to maintain the component's reliability, durability and safety. To this end, a solution is to reduce the disk-stator gap. Therefore, particular attention should be paid to the influence of possible rotor-stator interactions on the rotordynamics and the integrity of the turbomachinery, [1], [2], [3], [4], [6] and [7]. Rotor-stator interaction can be initiated by a mass-unbalance increase, high acceleration or variation of speed of rotation, reduced clearance, external excitation or damping loss. Several behaviours have been experimented and simulated. Partial and full annular rub in forward whirl (FW) generate a periodic or quasi-periodic motions which are not too dangerous for machine safety. Dry Whirl and Dry Whip are backward whirl motions. Dry Whirl occurs in the configuration of rolling without sliding, natural frequency increase with rotation speed. The Dry Whip is the most destructive phenomenon, its amplitude highly grows up and rub-motion is self-excited. Rotating machine usually adopt different behaviours before starting Dry Whip. Most of the works in the literature focusing on the rotor/stator interaction depict the different behaviour induced by the nonlinear problem. Few of them are able to predict with accuracy the global behaviour of a rotor-bearing system.

This paper focuses on the influence of a disk-stator interaction on the rotordynamics. The approach is theoretical and numerical. First a forced SDOF model containing two stop-ends is established to study in detail penalty contact combined with two step-by-step time integration algorithms. The numerical solutions are compared with an analytical response. Second, a three-disk rotor bearing-system is subject to a sudden disk-stator contact. The two numerical algorithms (Linear Acceleration by Newmark and Method of Central Differences) are
retained and associated to four contact models: three types of smooth penalty methods and the Lagrange multiplier method. Spectrograms of responses in the time domain are presented in order to analyse the spectral evolution of the rotor.

## 2. SDOF WITH TWO STOP-ENDS

In order to test the association of a penalty contact model and time integration schemes, the equation of motion of the SDOF system are derived and solved. Analytical and numerical results are compared [8].

### 2.1. Presentation of the SDOF system

Let the SDOF system shown by Fig. 1 be defined by M, K and C, its mass, stiffness and damping parameters, respectively. It is able to impact the two stop-ends when subject to the harmonic force $F(t) = F_0 \cos(\omega t)$. Stiffness $k_c$ and gap $g$ are the two parameters which define the contact.

<figure><img src="image_2.png" /><figcaption>Fig. 1: SDOF with two stop-ends</figcaption></figure>

The parameter values are listed on Table 1. A contact stiffness $k_c=1 \times 10^{10}$ N.m$^{-1}$ is chosen to simulate interaction between the mass and the two stop-ends.

<table><caption>Tab. 1: Parameters of SDOF system</caption><thead><tr><th>K [N.m<sup>-1</sup>]</th><th>C [N.s.m<sup>-1</sup>]</th><th>M [kg]</th><th>ω<sub>dsc</sub> [rad.s<sup>-1</sup>]</th><th>ω [rad.s<sup>-1</sup>]</th><th>F<sub>0</sub> [N]</th><th>g [m]</th><th>simulation time [s]</th></tr></thead><tbody><tr><td>1E9</td><td>50</td><td>100</td><td>3162.3</td><td>1256.6</td><td>201240</td><td>4x10<sup>-4</sup></td><td>0.5</td></tr></tbody></table>

### 2.2. Analytical solution

The equations of motion of the SDOF system are the following:

$$M\ddot{x} + C\dot{x} + \tilde{K}x = F(t) + \tilde{\varepsilon}.k_c.g \quad (1)$$

with

$$\tilde{K} = K \text{ and } \tilde{\varepsilon} = 0 \text{ if } |x| \le g \quad (2)$$

$$\tilde{K} = K + k_c \text{ and } \tilde{\varepsilon} = \text{sign}(x) \text{ else,} \quad (3)$$
The analytical displacement is expressed as:

$$ 
\begin{gathered}
x(t) = e^{-\frac{1}{\tau}} (A_i \cos(\tilde{\omega}_d \cdot t) + B_i \sin(\tilde{\omega}_d \cdot t)) + \tilde{\alpha} \cos(\omega \cdot t) + \tilde{\beta} \sin(\omega \cdot t) + \varepsilon \frac{k_c}{\tilde{K}} \\
\forall t \in [tc_i; tc_{i+1}]
\end{gathered}
\quad (4)
$$

Time constant, natural frequency, damping ratio and damped frequency are expressed respectively in the following:

$$ \tau = \frac{2.M}{C} \quad (5) $$

$$ \tilde{\omega}_0 = \sqrt{\frac{\tilde{K}}{M}} \quad (6) $$

$$ \tilde{\xi} = \frac{C}{2\sqrt{\tilde{K}.M}} \quad (7) $$

$$ \tilde{\omega}_d = \tilde{\omega}_0 (1 - \tilde{\xi}^2) \quad (8) $$

The contact changes parameter $\tilde{K}$ and natural frequency $\tilde{\omega}_0$ of the system. Consequently, parameters $A_i$, $B_i$, $\tilde{\omega}_d$, $\tilde{\alpha}$ and $\tilde{\beta}$, have to be updated at the beginning of each phase. The particular solution, such as:

$$ 
\left\{
\begin{aligned}
x_p(t) &= \tilde{\alpha} \cos(\omega t) + \tilde{\beta} \sin(\omega t) \\
\dot{x}_p(t) &= -\omega_i [\tilde{\alpha} \sin(\omega t) - \tilde{\beta} \cos(\omega t)] \\
\ddot{x}_p(t) &= -\omega^2 [\tilde{\alpha} \cos(\omega t) + \tilde{\beta} \sin(\omega t)]
\end{aligned}
\right.
\quad (9)
$$

is introduced in Equation (10):

$$ M \ddot{x}_p(t) + C \dot{x}_p(t) + \tilde{K} x_p(t) = F_0 \cos(\omega t) \quad (10) $$

The identification of the trigonometric terms yields:

$$ \begin{bmatrix} \tilde{\alpha} \\ \tilde{\beta} \end{bmatrix} = \begin{bmatrix} \tilde{K} - \omega^2 \cdot M & C \cdot \omega \\ -C \cdot \omega & \tilde{K} - \omega^2 \cdot M \end{bmatrix}^{-1} \begin{bmatrix} F_0 \\ 0 \end{bmatrix} \quad (11) $$

There are only two possible values of $\tilde{\alpha}$ and $\tilde{\beta}$ which satisfy conditions: $|x| \le g$ and $|x| > g$. The time sequence $(tc_i)_{i \in \mathbb{N}}$ defines dates of system changes:

$$ tc_{i+1} = \min\{t \in \mathbb{R}_+ / (t > tc_i) \& (x(t) = \pm g)\} \quad (12) $$

The condition $\{t \in \mathbb{R}_+ / (t > tc_i) \& (x(t) = \pm g)\} = \{\}$ means that the SDOF system is not restrained by the stop-ends. The analytical solution is a piecewise function. It's assumed that displacement and velocity are continuous especially at each transition phase. The velocity law is expressed as:

$$ \dot{x}(t) = e^{-\frac{1}{\tau}} \left[ \left( \tilde{\omega}_d \cdot B_i - \frac{A_i}{\tau} \right) \cos(\tilde{\omega}_d \cdot t) - \left( \tilde{\omega}_d \cdot A_i + \frac{A_i}{\tau} \right) \sin(\tilde{\omega}_d \cdot t) \right] + \dot{x}_p(t) \quad (13) $$
$(A_i)_{i \in N}$ and $(B_i)_{i \in N}$ are real sequences and are computed using:

$$ \begin{bmatrix} A_i \\ B_i \end{bmatrix} = \left[ - \begin{array}{c} \cos(\tilde{\omega}_d \cdot t) \\ (\tilde{\omega}_d \cdot \sin(\tilde{\omega}_d \cdot t) + \frac{1}{\tau} \cos(\tilde{\omega}_d \cdot t)) (\tilde{\omega}_d \cdot \sin(\tilde{\omega}_d \cdot t) - \frac{1}{\tau} \cos(\tilde{\omega}_d \cdot t)) \end{array} \right]^{-1} \begin{bmatrix} V_i \\ W_i \end{bmatrix} \quad (14) $$

$$ \begin{bmatrix} V_i \\ W_i \end{bmatrix} = e^{i\omega/c} \left[ \begin{array}{c} x(t_c) - \tilde{\alpha} \cdot \cos(\omega \cdot t) - \tilde{\beta} \cdot \sin(\omega \cdot t) - \tilde{\epsilon} \cdot \frac{k_c}{K} \cdot g \\ \dot{x}(t_c) + \omega_c (\tilde{\alpha} \cdot \sin(\omega \cdot t) - \tilde{\beta} \cdot \cos(\omega \cdot t)) \end{array} \right] \quad (15) $$

The chosen set of initial condition is:

$$ \begin{cases} tc_0 = t_0 = 0 \\ x藝 = x(0) = 0 \\ \dot{x}(tc_0) = \dot{x}(0) = 0 \end{cases} \quad (16) $$

## 2.3. Numerical Solution

Numerical results are calculated with Linear Acceleration by Newmark (LAN) and Method of the Central Differences (MCD) [9], [10] and [11]. Contact is simulated with penalty method and appears during transient motion. A constant step-by-step time of $1 \times 10^{-6}$ s is used during simulations.

## 2.4. Analytical vs numerical results

The displacement of the SDOF is plotted on Fig. 2. The two stop-ends are represented by two horizontal dotted lines. In order to make the analysis easier, only last 0.05 s time history is shown: analytical and all numerical solutions perfectly match. A contact stiffness of $1 \times 10^{11}$ N.m$^{-1}$ and $5 \times 10^9$ N.m$^{-1}$ were also computed and results perfectly match too.

The Linear Acceleration by Newmark method (LAN) and the Method of the Central Difference (MCD) give very satisfactory results in comparison with the analytical result. The LAN algorithm is retained in the following investigation.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Time history of analytical & numerical displacements</figcaption></figure>
### 3. THREE-DISKS ROTOR WITH ONE CONTACT

In this section the dynamics of a three-rigid-disk rotor system subject to one contact is investigated. The objective is to compare the influence of four following types of contact models on the rotordynamics [12], [13] and [14]: a- Damped penalty, b- Smooth penalty, exponential law, c- Smooth penalty, power law, d- Lagrange Multiplier.

Particular attention is paid to two different behaviours generated by the contact: partial and full annular rubs in forward precession. Spectrograms of each simulation are then shown.

#### 3.1. Presentation of the system

The academic three-disk rotor in bending, presented in [15], is made of a steel shaft with a circular constant section and a 1.3m length. It is supported by two linear bearings, see Fig. 3, which their damping and stiffness parameters do not depend on the speed of rotation. Disks and bearing properties are described respectively in Tables 2 and 3. The operating rate is 6000 rpm (100 Hz). The stator is assumed to be rigid with no motion.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Three-Disk rotor scheme [15]</figcaption></figure>

<table><caption>Tab. 2: Disks properties</caption><thead><tr><td>Disk</td><td>Mass [kg]</td><td>Polar inertia [kg.m<sup>2</sup>]</td><td>Diametral inertia [kg.m<sup>2</sup>]</td></tr></thead><tbody><tr><td>D1</td><td>14.58</td><td>0.123</td><td>0.0646</td></tr><tr><td>D2</td><td>45.95</td><td>0.976</td><td>0.498</td></tr><tr><td>D3</td><td>55.13</td><td>1.176</td><td>0.602</td></tr></tbody></table>

<table><caption>Tab. 3: Bearing properties</caption><tr><td>k<sub>xx</sub>=k<sub>zz</sub> [N.m<sup>-1</sup>]</td><td>c<sub>xx</sub>=c<sub>zz</sub> [N.s.m<sup>-1</sup>]</td></tr><tr><td>7x10<sup>7</sup></td><td>1x10<sup>3</sup></td></tr></table>

The bending FE model is composed with 13 identical Euler-Bernoulli shaft elements and 56 DOF, (two lateral displacements and two slopes at each node). A 200mm.g unbalance mass is located on the middle disk. Free of contact, the two first natural frequencies of the rotor are associated to the first and second bending modes. Fig. 4 shows these two modes at 6000 rpm i.e. 63.4 Hz & 182.2 Hz.
<figure><img src="image_5.png" /><figcaption>Fig. 4: Campbell Diagram, first and second bending mode shapes</figcaption></figure>

## 3.2. Contact models

The contact is introduced by a sudden increase of mass unbalance. It is assumed that the contact occurs at the middle disk location. Let the rotor and stator sections be considered. Fig. 5 shows a rotor and stator penetration. Positions of their centres are defined by $(X_S; Z_S)$ and $(X_R; Z_R)$ respectively.

<figure><img src="image_6.png" /><figcaption>Fig. 5: Rotor section penetrating stator section</figcaption></figure>
Let $\vec{n}$ and $\vec{t}$ be respectively the normal and tangential contact vectors. They are
defined by:

$$
\left\{
\begin{aligned}
\vec{n} &= \frac{1}{\sqrt{(X_R - X_S)^2 + (Z_R - Z_S)^2}} \begin{bmatrix} X_R - X_S \\ 0 \\ Z_R - Z_S \end{bmatrix} \\
\vec{t} &= \frac{1}{\sqrt{(X_R - X_S)^2 + (Z_R - Z_S)^2}} \begin{bmatrix} Z_R - Z_S \\ 0 \\ -(X_R - X_S) \end{bmatrix}
\end{aligned}
\right.
\qquad (17)
$$

At any time, the dynamic gap is:

$$
g_n = \delta = \sqrt{(X_R - X_S)^2 + (Z_R - Z_S)^2} - g \quad (18)
$$

Let $\mu_{\text{dry}}$ be the dry friction coefficient and $\vec{V}_R$ the peripheral speed of the disk.
Knowing that the friction coefficient is defined such as:

$$
\tilde{\mu} = \mu_{\text{dry}} \cdot \text{sign}(\vec{V}_R \cdot \vec{t}), \quad (19)
$$

the contact models tested during simulations are the following:

(a) Damped penalty with $k_c = 1 \times 10^{10}$ N/m, $c_c = 1 \times 10^4$ N.s.m⁻¹

$$
\vec{F}_{\text{contact}} = -(k_c \cdot g_n + c_c \cdot \dot{g}_n) [\vec{n} + \tilde{\mu} \cdot \vec{t}] \quad (20)
$$

(b) Smooth penalty (exponential law) with $k_c = 1 \times 10^{10}$ N.m⁻¹

$$
\vec{F}_{\text{contact}} = k_c \cdot g_n \left( e^{-g_n / g_{\text{ref}}} - 1 \right) \left[ \vec{n} + \tilde{\mu} \cdot \vec{t} \right] \quad (21)
$$

(c) Smooth penalty (Hertz contact) with $k_c = 1 \times 10^{12} \text{ N}^{2/3} \cdot \text{m}^{-2/3}$

$$
\vec{F}_{\text{contact}} = -(k_c \cdot g_n)^{3/2} [\vec{n} + \tilde{\mu} \cdot \vec{t}] \quad (22)
$$

(d) Lagrange multiplier is able to respect condition of no-penetration. A constraint equation Eq. (30) is used and derivate to obtain expression of the Jacobian:

$$
h_n = g_n^2 = (X_R - X_S)^2 + (Z_R - Z_S)^2 - g = 0 \quad (23)
$$

Lagrange multiplier $\lambda_n$ is introduced with the Jacobian matrix $J$ in the
equation of motion such as:

$$
\left\{
\begin{array}{l}
M \cdot \ddot{q}_n + J^T \cdot \lambda_n = F_n^{\text{ext}} - C \cdot \dot{q}_n - K \cdot q_n \\
\\
\begin{bmatrix} M & J^T \\ J & 0 \end{bmatrix} \{ \ddot{q}_n \} = \begin{Bmatrix} F_n^{\text{ext}} - C \cdot \dot{q}_n - K \cdot q_n \\ D(q_n, \dot{q}_n) \end{Bmatrix}
\end{array}
\right.
\qquad (24)
$$

Contact force is valued with Jacobian and Lagrange multiplier [14]:

$$
\vec{F}_{\text{contact}} = -J^T \cdot \lambda_n \quad (25)
$$

The force vector $\vec{F}_{\text{contact}}$ is added in the second member of the equation of motion.
### 3.3. Numerical results

The motion of the 3-disk rotor-bearing system is described at time $t_n$ by the following set of ordinary differential equations:

$$ \underline{M} \dot{\underline{q}}_n + (\underline{G}(\Omega) + \underline{C}) \dot{\underline{q}}_n + \underline{K} \underline{q}_n = \underline{F}_n^{\text{ext}} \quad (26) $$

with $M$ the mass matrix, $G$ the gyroscopic matrix, $C$ the matrix due to the shaft and bearing damping, $K$ the matrix due to the shaft and bearing stiffness, $F^{\text{ext}}$, the vector containing the unbalance mass excitation and the contact force. Numerical results are calculated with the LAN method. A constant step-by-step time of $1 \times 10^{-5}$ s is selected. The two contact cases, generating partial and full annular rubs, use parameters presented in Table 4. Reduction of initial gap generate full annular rub in FW.

<table><caption>Tab. 4: Contacts properties</caption><thead><tr><th>Contact behaviours</th><th>k<sub>c</sub> [N.m<sup>-1</sup>]</th><th>g[m]</th></tr></thead><tbody><tr><td>Impact series</td><td>1x10<sup>10</sup></td><td>5.6x10<sup>-5</sup></td></tr><tr><td>Full annular rub</td><td>1x10<sup>10</sup></td><td>4x10<sup>-5</sup></td></tr></tbody></table>

#### 3.3.1. Impact Series

After 2 seconds, the rotor impacts few times the rigid stator. Force contact values and number of rebounds vary between different contact methods. The first impact is the most violent. Rebound duration is nearly 0.3 seconds before adopting a linear behaviour. The rotor does not impact the rigid stator anymore.

<figure><img src="image_7.png" /><figcaption>Fig. 6: Spectrograms. (a) Damping penalty, (b) Smooth penalty 1, (c) Smooth penalty 2, (d) Lagrange multiplier methods</figcaption></figure>
Spectrograms of Fig. 6 exhibit a large frequency spectrum during 0.3 s after the first impact. Then the rotor is mainly led by its first natural bending mode and by the mass unbalance excitation. Even if the frequency of the 2^nd^ mode appears in the first instants after the first impact, it is the first mode shape which remains the motion. The energy levels depend on the choice of contact models.

### 3.3.2. Full Annular Rub

The orbits of the rotor section centre at the middle disk location during the contact are calculated. In each case, the disk starts bouncing into the rigid stator. Then it follows full annular rub behaviour in FW. The penetration norm depends on the contact model chosen. The smooth penalty methods permit a deeper penetration than the damped penalty method. Lagrange multiplier does not authorize any penetration.

This full annular rub has a 2200 N constant force which does not depend on the type of contact law.

Spectrograms of the time-history responses are displayed in Fig. 7. Before the sudden mass unbalance generated at t=2s, the spectrograms show the presence of the first mode of the rotor and the mass unbalance frequency, while after they show mostly the order 2/3, 1 and 4/3 of the speed of rotation. At the end of the numerical simulation, the rotor reaches a steady state motion composed only with the order one of the speed of rotation. The transient state just after the rotor-stator impact (2s), exhibits a non linear behaviour: there is the presence of multiple third harmonics. All spectrograms are similar, and describe the same behaviour regardless the contact method chosen.

<figure><img src="image_8.png" /><figcaption>Fig. 7: Spectrograms. (a) Damping penalty, (b) Smooth penalty 1, (c) Smooth penalty 2, (d) Lagrange multiplier</figcaption></figure>
## 4. CONCLUSION

A review of integration time schemes and 4 contact methods have been presented. The SDOF system permits investigating robustness and stability of different integration time schemes regarding step-time. Numerical and analytical algorithms give the same responses. Linear Acceleration by Newmark has been chosen.

After a rotor/stator contact, two different behaviours have been observed with two gap values. An initial gap of $5.6 \times 10^{-5}$ m leads to a finite impacts series of the rotor within its casing. A $4 \times 10^{-5}$ gap generates a full annular rub in FW led by the mass unbalance excitation. All contact methods reviewed, simulate roughly the same behaviour with various transient motions.

## 5. ACKNOWLEDGEMENTS

The author would like to thank SAFRAN-Turbomeca for its support and permission to publish these results.

## REFERENCE LIST

[1] JIANG J., « The Analytical Solution and The Existence Condition of Dry Friction Backward Whirl in Rotor-to-Stator Contact Systems », *Journal of Vibrations and Acoustics*, Vol. 129, pp. 260-264, 2007.

[2] CHILDS D. W. and BHATTACHARYA A., « Prediction of Dry-Friction Whirl and Whip Between a rotor and a Stator », *Journal of Vibration and Acoustics*, Vol. 129, pp. 355-362, 2007.

[3] AHMAD S., « Rotor Casing Contact Phenomenon in Rotor Dynamics -- Literature Survey », *Journal of Vibration and Control* (2010), Vol. 16(9), pp. 1369-1377, 2010.

[4] GAO JIN JI and QI QI MIN, « Rotor-to-stator rub vibration in centrifugal compressor », *Dynamics and Control*, N86-30178, pp. 235-244, China.

[5] BARTHA A. R., *Dry friction. Backward whirl of rotors*. Thesis, Swiss Federal Institute of Technology Zurich, Zurich, Switzerland, 2000.

[6] DONALD E., BENTLY, JOHN J. YU, GOLDMAN P. and MUSZYNSKA A., « Full Annular Rub in Mechanical Seals Part I: Experimental Results », *International Journal of Rotating Machinery*, Vol. 8(5), pp. 319-328, 2002.

[7] DONALD E., BENTLY, JOHN J. YU and GOLDMAN P., « Full Annular Rub in Mechanical Seals Part II: Analytical Study », *International Journal of Rotating Machinery*, Vol. 8(5), pp. 319-328, 2002.

[8] BROGLIATO B. and ORHANT P., « Contact Stability Analysis of a One Degree-of-Freedom Robot », *Dynamics and Control*, Vol. 8, pp. 37-53, 1998.

[9] ZHANG W-M., MENG G., CHEN D., ZHOU J-B. and CHEN J-Y., « Nonlinear dynamics of a rub-impact micro-rotor system with scale-dependent friction model », *Journal of Sound and Vibration* (2007), Vol. 309(3-5), pp. 756-777, 2008.

[10] ROQUES S., LEGRAND M., CARTRAUD P., STOISSER C. and PIERRE C., « Modeling of a rotor speed transient, response with radial rubbing », *Journal of Sound and Vibration*, Vol. 329, pp. 527-546, 2007.

[11] XIE Y. M., « An Assessment of Time Integration Schemes for Non-Linear Dynamic Equations », *Journal of Sound and Vibration*, Vol. 192(1), pp. 321-331, 1996.

[12] ULOVIC S., ZIVKOVIC M. and GRUJOVIC N., « Contact Problems Based on the Penalty Method », *Scientific Technical Review*, Vol. LVIII, No. 3-4, pp. 33-37, 2008.

[13] YUAN Z., CHU F. and HAO R., « Simulation of rotor's axial rub-impact in full degrees of freedom », *Mechanism and Machine Theory*, Vol. 42, pp. 763-775, 2007.
[14] SAHINKAYA N., « Multiple Sliding and Rolling Contact Dynamics for a Flexible Rotor/Magnetic Bearing System », *ASME Transactions on Mechatronics*, Vol. 12(2), pp. 179-189, 2007.

[15] LALANNE M. and FERRARIS G., Rotordynamics Prediction in Engineering: 2nd Edition. Baffins Lane, Chichester, England: John Wiley & Sons Ltd., 1998.