# Application d'un modèle numérique de roulement non linéaire aux calculs d'un système birotor

Mikhail Gouskov, Jean-Jacques Sinou, Fabrice Thouverez, Oleg Naraïkine

► To cite this version:

Mikhail Gouskov, Jean-Jacques Sinou, Fabrice Thouverez, Oleg Naraïkine. Application d'un modèle numérique de roulement non linéaire aux calculs d'un système birotor. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01812937

HAL Id: hal-01812937

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

Submitted on 12 Jun 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC0 1.0 - Universal - International License
Application d'un modèle numérique de
roulement non linéaire aux calculs d'un
système birotor

Mikhail Gouskov a — J-Jacques Sinou b — Fabrice Thouverez c
— Oleg S. Naraïkine d

a, b et c Laboratoire de Tribologie et Dynamique des Systèmes CNRS-UMR 5513
Equipe D2S, Ecole Centrale de Lyon
36 avenue Guy de Collongue, 69134 Ecully Cedex, France.

a et d Département de Mécanique appliquée,
Université Technique d'Etat de Moscou Bauman
Russie, 107005, Moscou, rue 2-ya Baoumanskaya, 5

a mikhail.gouskov@ec-lyon.fr

b jean-jacques.sinou@ec-lyon.fr

c fabrice.thouverez@ec-lyon.fr

d oleg@hexa.ru

RÉSUMÉ. Une méthode permettant d'estimer les efforts non-linéaires dans les roulements lors de calculs en dynamique des rotors est proposée. Elle consiste à précalculer la caractéristique effort-déplacement du roulement afin de pouvoir ensuite l'utiliser dans des calculs dynamiques à travers une interpolation polylinéaire. Un exemple de calcul sur un système birotor comportant un roulement à rouleau avec jeu radial est traité.

ABSTRACT. A method for computing the non-linear bearing force terms in the rotordynamics calculations is proposed. It consists in computing the multidimensional non-linear bearing stiffness chart beforehand. While the dynamical calculations are performed, the polylinear interpolation is used. An example of a dual rotor system with a simplified roller bearing model is considered.

MOTS-CLÉS : dynamique des rotors, roulement non-linéaire, interpolation polylinéaire, birotor

KEYWORDS: rotordynamics, non-linear bearing, polylinear interpolation, dual rotor
# 1. Introduction

Le développement actuels des machines tournantes tend vers une augmentation de la puissance et des vitesses tout en réduisant leurs masse et gabarits. Cela entraîne une compréhension fine de la dynamique de tels systèmes qui n’était pris en compte que sommairement auparavant. En particulier, la modélisation des non-linéarités dues à la présence de roulement fait actuellement l’objet d’études de plus en plus détaillées. Par exemple, un modèle détaillé de roulement à billes, tenant compte de la raideur non-linéaire et des effets centrifuges et gyroscopiques est présenté par Jones en 1960. Il a été complété et généralisé au cas des roulements à rouleaux par de Mul et al, 1989 en appliquant une technique de découpage (slicing) présentée par Andréason, 1973. Ces développements ont ensuite été appliqués aux études de la dynamique d’un rotor supporté par des roulements à billes par Liew et al, 2000 par le biais de l’intégration temporelle. De plus, dans le cas de systèmes multirotors un roulement interarbre est souvent utilisé, représenté par un roulement à rouleaux. Couramment les études des systèmes multirotors ne considèrent alors que des modélisations linéaires des roulements interarbres (Ferraris et al, 1986; Gupta et al, 1996). La prise en compte de la non-linéarité de roulements (inter-arbre ou autre) augmente alors le temps de calcul afin d’obtenir le comportement vibratoire du système.

Aussi, cette étude se propose de calculer la dynamique non-linéaire d'un système bi-rotor en s'appuyant sur une interpolation polylinéaire des efforts non-linéaires au niveau du roulement afin de gagner en temps de calcul et d'éviter un calcul supplémentaire lié à la prise en compte des réactions non-linéaires des roulements.

# 2. Système bi-rotor et méthode de calcul

Le système étudié représente un modèle simplifié birotor, composé de deux rotors « ponctuels », reliés par un roulement interarbre. Les équations du système sont:

$$ 
\begin{align} 
m_i \ddot{x}_i + c_i \dot{x}_i + k_i x_i + f_x(x_1, y_1, \theta_1, \psi_1, x_2, y_2, \theta_2, \psi_2) &= (me)_i \omega_i^2 \cos \omega_i t \nonumber \\ 
m_i \ddot{y}_i + c_i \dot{y}_i + k_i y_i + f_y(x_1, y_1, \theta_1, \psi_1, x_2, y_2, \theta_2, \psi_2) &= (me)_i \omega_i^2 \sin \omega_i t - m_i g \tag{1} \\ 
J_i \ddot{\theta}_i + c_{i0} \dot{\theta}_i - \omega_i J_{pi} \psi_i + k_{0i} \theta_i + f_0(x_1, y_1, \theta_1, \psi_1, x_2, y_2, \theta_2, \psi_2) &= 0 \nonumber \\ 
J_i \dot{\psi}_i + c_{i0} \dot{\psi}_i + \omega_i J_{pi} \dot{\theta}_i + k_{0i} \psi_i + f_\psi(x_1, y_1, \theta_1, \psi_1, x_2, y_2, \theta_2, \psi_2) &= 0 \nonumber 
\end{align} 
$$

avec $i = 1,2$ numéros du rotor, $m_i$ masses des rotors, $J_i$ moments d'inertie diamétrale, $J_{pi}$ moments d'inertie polaire, $c_i$ et $c_{0i}$ coefficients de dissipation linéique et angulaire, $k_i$ et $k_{0i}$ raideurs linéique et angulaire, $(me)_i$ balourds, $x_i$ et $y_i$ déplacements transversaux selon les axes $x$ et $y$, $\theta_i$ et $\psi_i$ rotations autour de ces axes, $\omega_i$ vitesses de rotation des rotors, $f_x, f_y, f_0, f_\psi$ les composantes respectives du vecteur de l'effort non-linéaire, du au roulement non-linéaire. Ainsi, pour
chacune des quatre directions, le système peut être conçu comme un double masse-ressort où les masses sont reliées par un ressort non-linéaire (une telle représentation pour une direction transversale x ou y est donnée sur la Figure 1(a)).

En adimensionnant (1), on obtient

$$
\begin{align}
\xi_i'' + 2\zeta_i \xi_i' + \omega_i^2 \xi_i + \phi_{\xi} (\xi_1, \eta_1, \theta_1, \psi_1, \xi_2, \eta_2, \theta_2, \psi_2) &= \varepsilon_i \cos \Omega_i \tau \nonumber \\
\eta_i'' + 2\zeta_i \eta_i' + \omega_i^2 \eta_i + \phi_{\eta} (\xi_1, \eta_1, \theta_1, \psi_1, \xi_2, \eta_2, \theta_2, \psi_2) &= \varepsilon_i \sin \Omega_i \tau - \Gamma_i \tag{2} \\
\theta_i'' + 2\zeta_{i\theta} \theta_i' - \Omega_i \gamma_i \psi_i' + \omega_{i\theta}^2 \theta_i + \phi_0 (\xi_1, \eta_1, \theta_1, \psi_1, \xi_2, \eta_2, \theta_2, \psi_2) &= 0 \nonumber \\
\psi_i'' + 2\zeta_{i\phi} \psi_i' + \Omega_i \gamma_i \theta_i' + \omega_{i\phi}^2 \psi_i + \phi_{\psi} (\xi_1, \eta_1, \theta_1, \psi_1, \xi_2, \eta_2, \theta_2, \psi_2) &= 0 \nonumber
\end{align}
$$

Les efforts non-linéaires sont calculés par interpolation d'un tableau de valeurs, calculé préalablement dans les axes polaires ρ, φ, θρ, θφ (Figure 1(b)). Les composantes des réactions non-linéaires φρ, φφ, φθρ, φθφ sont dont stockées dans des tableaux définis par les $\bar{\phi}_{\rho}^{ijkl}$, $\bar{\phi}_{\phi}^{ijkl}$, $\bar{\phi}_{\theta\rho}^{ijkl}$ et $\bar{\phi}_{\theta\phi}^{ijkl}$ qui comportent quatre entrées ($\rho^i, \phi^j, \theta_{\rho}^k, \theta_{\phi}^l$). Le passage des variables $\xi_1, \eta_1, \theta_1, \psi_1, \xi_2, \eta_2, \theta_2, \psi_2$ aux variables $\rho, \phi, \theta_{\rho}, \theta_{\phi}$ se fait à partir des relations:

$$
\begin{align}
\rho &= \sqrt{(\xi_2 - \xi_1)^2 + (\eta_2 - \eta_1)^2} = (\xi_2 - \xi_1) \cos \theta_{\xi\eta} + (\eta_2 - \eta_1) \sin \theta_{\eta\xi} \nonumber \\
\varphi &= (\theta_{\xi\eta} + [\Omega_1 + \varpi(\Omega_2 - \Omega_1)]\tau) \mod (2\pi/Z) \tag{3} \\
\vartheta_{\rho} &= (\theta_2 - \theta_1) \cos \theta_{\xi\eta} + (\psi_2 - \psi_1) \sin \theta_{\eta\xi} \nonumber \\
\vartheta_{\phi} &= -(\theta_2 - \theta_1) \sin \theta_{\xi\eta} + (\psi_2 - \psi_1) \cos \theta_{\eta\xi} \nonumber \\
\phi_{\xi} &= \phi_{\rho} \cos \theta_{\xi\eta} - \phi_{\varphi} \sin \theta_{\xi\eta}; \quad \phi_{\eta} = \phi_{\rho} \sin \theta_{\eta\xi} + \phi_{\varphi} \cos \theta_{\eta\xi} \nonumber \\
\phi_{\varphi} &= \phi_{\theta\rho} \cos \theta_{\xi\eta} - \phi_{\theta\phi} \sin \theta_{\xi\eta}; \quad \phi_{\psi} = \phi_{\theta\rho} \sin \theta_{\eta\xi} + \phi_{\theta\phi} \cos \theta_{\eta\xi} \tag{4}
\end{align}
$$

L'interpolation s'effectue de façon polylinéaire. Soit un point $\rho_\circ, \phi_\circ, \vartheta_{\rho_\circ}, \vartheta_{\phi_\circ}$, pour lequel il faut calculer la caractéristique $\phi_\rho$ du roulement, avec $i, j, k$ et $l$ les numéros d'entrées du tableau $\bar{\phi}_\rho^{ijkl}$ tels que $\rho^i \le \rho_\circ < \rho^{i+1}$, $\varphi^j \le \varphi_\circ < \varphi^{j+1}$, $\vartheta_\rho^k \le \vartheta_{\rho_\circ}^k < \vartheta_\rho^{k+1}$, $\vartheta_\phi^k \le \vartheta_{\phi_\circ}^k < \vartheta_\phi^{k+1}$, la valeur interpolée $\tilde{\phi}_\rho^{ijkl}$ est alors recherchée sous la forme :

$$
\begin{align}
\tilde{\phi}_{\rho}^{ijkl} ={}& a_1 \bar{\phi}_{\rho}^{ijkl} + a_2 \bar{\phi}_{\rho}^{i+1 jkl} + \dots + a_5 \bar{\phi}_{\rho}^{ijkl+1} + a_6 \bar{\phi}_{\rho}^{i+1 j+1 kl} + \dots + a_{12} \bar{\phi}_{\rho}^{ijk+1 l+1} \nonumber \\
& + a_{13} \bar{\phi}_{\rho}^{i+1 j+1 k+1 l} + \dots + a_{15} \bar{\phi}_{\rho}^{ij+1 k+1 l+1} + a_{16} \bar{\phi}_{\rho}^{i+1 j+1 k+1 l+1} \tag{5}
\end{align}
$$

On définit alors une fonction quadrilinéaire $\Phi$ dans le domaine quadridimensionnel $[0, 1]^4$ de manière à satisfaire les conditions aux bornes du domaine :

$$
\Phi(\delta_1, \delta_2, \delta_3, \delta_4) = \bar{\phi}_{\rho}^{i+\delta_1, j+\delta_2, k+\delta_3, l+\delta_4}, \delta_{1,2,3,4} \in \{0,1\} \quad (6)
$$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Système birotor (a) modélisation (b) conventions de signe</figcaption></figure>

En substituant l’argument de $\Phi$ dans (6) $\rho_c/(\rho^{i+1}-\rho^i)$, $\phi_c/(\phi^{i+1}-\phi^i)$, $\vartheta_{\rho_c}/(\vartheta_{\rho}^{k+1}-\vartheta_{\rho}^k)$, $\vartheta_{\phi_c}/(\vartheta_{\phi}^{i+1}-\vartheta_{\phi}^i)$, on obtient 16 équations linéaires, permettant de déterminer les coefficients $a_{1...16}$ de l’équation (5). Cette procédure revient alors à inverser une matrice triangulaire.

Dans le modèle présenté dans cette étude, une caractéristique approchée et analytique du roulement est proposée, de la forme :

$$ 
\begin{aligned}
\phi_{\rho} &= (\phi_{\rho 0} + \kappa \phi_{\vartheta 0}) (1 + \varepsilon_{VC} \cos \varphi); & \phi_{\vartheta} &= (\phi_{\vartheta 0} - \kappa \phi_{\rho 0}) (1 + \varepsilon_{VC} \sin \varphi) \\
\phi_{\vartheta \rho} &= (\phi_{\vartheta 0} - \kappa \phi_{\varphi 0}) (1 + \varepsilon_{VC} \cos \varphi); & \phi_{\vartheta \varphi} &= (\phi_{\vartheta 0} + \kappa \phi_{\varphi 0}) (1 + \varepsilon_{VC} \sin \varphi)
\end{aligned}
\quad (7) $$

avec $\kappa$ le coefficient de couplage entre les dimensions, $\varepsilon_{VC}$ l'amplitude relative de la variation de la souplesse,

$$ 
\begin{aligned}
\phi_{\rho 0} &= k_{\rho} (\rho - \delta)^{\alpha} h(\rho - \delta); & \phi_{\vartheta 0} &= 0 \\
\phi_{\vartheta 0} &= k_{\vartheta} \left(\sqrt{\vartheta_{\rho}^{2} + \vartheta_{\varphi}^{2}} - \delta_{\vartheta}\right)^{\alpha} h\left(\sqrt{\vartheta_{\rho}^{2} + \vartheta_{\varphi}^{2}} - \delta_{\vartheta}\right) \vartheta_{\rho} / \sqrt{\vartheta_{\rho}^{2} + \vartheta_{\varphi}^{2}} & (8) \\
\phi_{\vartheta 0} &= k_{\vartheta} \left(\sqrt{\vartheta_{\rho}^{2} + \vartheta_{\varphi}^{2}} - \delta_{\vartheta}\right)^{\alpha} h\left(\sqrt{\vartheta_{\rho}^{2} + \vartheta_{\varphi}^{2}} - \delta_{\vartheta}\right) \vartheta_{\varphi} / \sqrt{\vartheta_{\rho}^{2} + \vartheta_{\varphi}^{2}}
\end{aligned}
$$

où $\delta, \delta_\vartheta$ sont les jeux radial et angulaire, $h$ est la fonction de Heavyside, $k_\rho, k_\vartheta$ sont les coefficients de raideur non-linéaire, $\alpha$ est l'exposant de la non-linéarité. Cette fonction a été calculée sur une plage de valeurs $0 \le \rho \le \rho_0$, $0 \le \varphi \le 2\pi/Z$, $0 \le \vartheta_\rho \le \vartheta_{\rho 0}$, $0 \le \vartheta_\varphi \le \vartheta_{\varphi 0}$. Elle est ensuite utilisée dans les calculs dynamiques lors de la recherche du comportement vibratoire par intégration de type Runge-Kutta.

### 3. Application et résultats

On s'intéresse maintenant à l'application de cette procédure au cas d'un système bi-rotor (comme illustré en Figure 1(a)) correspondant au jeu de données suivant : $\zeta_1 = \zeta_2 = \zeta_{10} = \zeta_{20} = 0,05$ ; $\omega_1 = \omega_2 = \omega_{10} = \omega_{20} = 1$ ; $\Omega_2 = 3\Omega_1$ ; $0 < \varepsilon_1 = \varepsilon_2 < 1$ ; $\Gamma_2 = 3\Gamma_1 = 3$ ; $0 < \delta < 1$ ; $\delta_\vartheta = 0$ ; $\kappa = 0,1$ ; $\varepsilon_{VC} = 0,1$.
Après intégration temporelle avec utilisation de l'interpolation polylinéaire, la réponse fréquentielle sur le degré de liberté horizontal pour la masse $m_1$ du système bi-rotor est donnée en Figure 2. Suivant les quatre directions les résonances excitées par le balourd du premier rotor tombent sur 2 pics situés approximativement à $\Omega_1=1$ et $\Omega_1=1.7$ (dans le cas linéaire $\Omega_1=1$ et $\Omega_1 = \sqrt{3}$).

<figure><img src="image_3.png" /><figcaption>Figure 2. Réponses à balourd

$a — \delta=0,5; \varepsilon_{1,2}=0,25 ; b — \delta=0,75; \varepsilon_{1,2}=0,75.$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. Trajectoires des rotors (cas de Fig. 2a) : a — rotor 1, $\Omega_1=2,5$, b — rotor 1, $\Omega_1=3$, c — rotor 2, $\Omega_1=2,5$, d — rotor 2 (dans le repère lié au rotor 1), $\Omega_1=2,6$</figcaption></figure>
Du fait de la présence d'une raideur non-linéaire, la position fréquentielle de la résonance dépend du balourd (niveau d'excitation) et décroît avec l'augmentation du jeu radial dans le roulement interarbre. Une résonance à la moitié de la fréquence de rotation du premier rotor est présente autour de Ω₁=3. La Figure 3 présente les trajectoires des rotors dans le repère absolu (a-c) et dans le repère lié au rotor 1 (d).

**4. Conclusion**

Une étude portant sur la dynamique non-linéaire des systèmes bi-rotor est présentée en considérant une estimation des efforts non-linéaires dans les roulements. Un abaque calculé préalablement permet alors d'éviter des itérations lors des calculs des réactions des corps roulants lors des études dynamiques. Une application sur un exemple birotor avec un roulement interarbre comportant un jeu et une raideur non-linéaire est présentée. Des phénomènes dus à la présence des non-linéarités inhérentes aux roulements sont observés.

**Remerciements**

Les auteurs tiennent à remercier l'ambassade de France en Russie pour le soutien apporté à ce travail.

**Références**

Andréason S. « Load distribution in a taper roller bearing arrangement considering misalignment », *Tribology*, 1973, 84-92.

Gupta K., Gupta K.D., Athre K. Unbalance Response of a Dual Rotor System: Theory and Experiment, *Journal of Vibration and Acoustics*, 1993, Vol. 115, 427-435.

Jones A.B., A General Theory for Elastically Constrained Ball and Radial Roller Bearings Under Arbitrary Load and Speed Conditions, *ASME Journal of Basic Engineering*, 1960, 309-320.

Harris T., *Rolling Bearing Analysis*, 4th edition — Wiley, 2001.

Ferraris G., Maisonnvuve V., Lalanne M. Prediction of the Dynamic Behavior of Non-Symmetric Coaxial Co- or Counter-rotating Rotors, *Journal of Sound and Vibration*, 1996, Vol. 195(4), 649-666.

Liew A., Feng N., Hahn E.J. Transient rotordynamic modelling of rolling element bearing systems, *Journal of Engineering for Gas Turbines and Power*, Vol. 124, 2002, 984-991.

De Mul J.M., Vree J.M., Maas D.A., Equilibrium and Associated Load Distribution in Ball and Roller Bearings Loaded in Five Degrees of Freedom While Neglecting Friction - Part I: General Theory and Application to Ball Bearings, *ASME Journal of Tribology*, 1989, Vol. 111, 142-148.

De Mul J.M., Vree J.M., Maas D.A., Equilibrium and Associated Load Distribution in Ball and Roller Bearings Loaded in Five Degrees of Freedom While Neglecting Friction - Part II: Application to Roller Bearings and Experimental Verification, *ASME Journal of Tribology*, 1989, Vol. 111, 149-155.