# Amélioration de la prédiction de l'interaction rotor/stator dans un moteur d'avion

Mathias Legrand, Bernard Peseux, Christophe Pierre, Eric Seinturier

► To cite this version:

Mathias Legrand, Bernard Peseux, Christophe Pierre, Eric Seinturier. Amélioration de la prédiction de l'interaction rotor/stator dans un moteur d'avion. 6e Colloque national en calcul des structures, École polytechnique Palaiseau, 2003, Giens, France. hal-01349530

HAL Id: hal-01349530

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

Submitted on 27 Jul 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 CC BY 4.0 - Attribution - International License
# Amélioration de la prédiction de l'interaction rotor/stator dans un moteur d'avion

M. Legrand<sup>a</sup>, B. Peseux<sup>a</sup>, C. Pierre<sup>b</sup> et E. Seinturier<sup>c</sup>

<sup>a</sup> LMM, Ecole Centrale de Nantes, 1 Rue de la Noë, 44321 Nantes cedex 3
{mathias.legrand,bernard.peseux}@ec-nantes.fr

<sup>b</sup> University of Michigan, Ann Arbor, U.S.A.
pierre@umich.edu

<sup>c</sup> Snecma Moteurs, site de Villaroche, 77550 Moissy Cramayel
eric.seinturier@snecma.fr

## Résumé

Dans ce travail, un outil numérique est développé : il permet de simuler une interaction modale rotor-stator, coïncidence vibratoire entre un mode de la roue aubagée et un mode du carter qui lui fait face. Lors de ce phénomène, le contact entre les deux structures entraîne un transfert d'énergie qui peut se traduire par une instabilité et une destruction du moteur. Les premiers résultats numériques confirment la relation analytique existante qui prédit les vitesses de rotation pour lesquelles il y a risque d'interaction modale.

## 1 Introduction

Dans le cadre du dimensionnement des carters fan, le bureau d'études doit concevoir des structures capables de résister à différents types de sollicitations dynamiques, du contact modal rotor-stator à la rétention d'une aube de soufflante. La tenue des carters à la perte d'aube nécessite leur rigidification par des nervures externes. Il est montré que la raideur du carter joue un rôle non négligeable dans le phénomène d'interaction modale rotor-stator. Le contact entre les deux structures, dû à un balourd par exemple, entraîne un transfert d'énergie entre les deux structures qui peut se traduire par une instabilité (une augmentation exponentielle des amplitudes de vibration) et une ruine des composants ou des brides boulonnées qui les retiennent. Malgré la pose d'un matériau abradable à faible tenue mécanique sur le carter afin d'éviter le contact direct métal-métal, il est nécessaire de trouver le meilleur compromis possible sur la rigidité de la pièce pour que soient assurées à la fois la rétention et l'insensibilité à l'interaction. Un outil simplifié permettant numériquement de détecter une situation d'interaction rotor-stator a été développé. La première partie concerne la modélisation des deux structures en jeu, la seconde partie explique le phénomène d'interaction modale tel qu'il est connu actuellement et la troisième et dernière partie détaille l'algorithme de résolution développé et présente quelques résultats.
# 2 Modélisation des structures

## 2.1 Modèle de la roue aubagée

Pour cette étude, la roue aubagée a été modélisée le plus simplement possible. Elle est constituée d'un nombre défini de poutres rigides fixées au disque du rotor par l'intermédiaire d'une raideur et connectées entre elles par des ressorts de type spiral : sur la figure 1, les raideurs $k_a$ agissent directement sur les paramètres $\alpha_i$, positions angulaires des aubes respectives. Ici, la raideur centrifuge due à la vitesse de rotation de la roue aubagée n'est pas prise

<figure><img src="image_2.png" /><figcaption>Figure 1 – Schéma d'une partie de la roue aubagée</figcaption></figure>

en compte. L'équation de la mécanique associée à $\alpha_i$ est :

$$m_a \ddot{\alpha}_i + (k_b + 2k_a)\alpha_i - k_a(\alpha_{i+1} + \alpha_{i-1}) = 0 \quad \text{avec} \quad m_a = \frac{ml_a^2}{3} \quad (1)$$

La généralisation de cette équation à toute la roue aubagée permet d'écrire un système matriciel où la matrice raideur est circulante et la matrice masse diagonale :

$$[M]\{\ddot{\alpha}\} + [K]\{\alpha\} = \{0\} \quad (2)$$

Les modes propres de la roue aubagée sont calculés à partir de ce modèle.

## 2.2 Modèle du carter

Le carter est directement discrétisé sur ses modes propres et les masses et raideurs modales respectives sont calculées grâce à un modèle éléments finis très simple d'un carter.

Lors du calcul, des termes d'amortissement ont été ajoutés dans les équations pour prendre en compte l'amortissement structurel de la roue aubagée et du carter. En plus, le système matriciel respectif de chaque structure n'est pas normé par rapport à sa matrice masse pour conserver le vrai rapport massique entre le stator et le rotor.
# 3 L'interaction modale

## 3.1 Définitions

L'axisymétrie de la roue aubagée, qui se traduit généralement par des matrices masse et raideur circulantes, confère à ses modes propres une propriété intéressante : ils ont des fréquences de vibration égales deux à deux [1]. Quelques modes font exception :

- le premier si le nombre de degrés de liberté est impair.

- le premier et le dernier si le nombre de degrés de liberté est pair.

Les modes à fréquence égale peuvent être combinés de façon judicieuse pour former des ondes de déformation rétrograde et avance. Pour ces deux modes, les angles de chaque aube peuvent être décrits par les équations :

$$ \alpha(\theta) = \cos(n_d(\theta - \theta_0)) \text{ et } \alpha(\theta) = \sin(n_d(\theta - \theta_0)) $$

où $\theta_0$ est un paramètre indéterminé et $n_d$ le nombre de diamètres des deux modes.

<figure><img src="image_3.png" /><figcaption>Figure 2 – Un des modes à 3 diamètres de chaque structure</figcaption></figure>

L'onde qui se propage dans le sens direct est alors définie par :

$$ \alpha_r(\theta, t) = \cos(\omega t) \cos(n_d(\theta - \theta_0)) + \sin(\omega t) \sin(n_d(\theta - \theta_0)) $$

$$ \alpha_r(\theta, t) = \cos(\omega t - n_d(\theta - \theta_0)) $$

et l'onde qui se propage dans le sens indirect, par :

$$ \alpha_a(\theta, t) = \cos(\omega t) \cos(n_d(\theta - \theta_0)) - \sin(\omega t) \sin(n_d(\theta - \theta_0)) $$

$$ \alpha_a(\theta, t) = \cos(\omega t + n_d(\theta - \theta_0)) $$

La vitesse de propagation de ces deux ondes vaut $\pm\omega/n_d$ dans le repère lié à la structure et $\Omega \pm \omega/n_d$ dans un repère fixe si la roue aubagée est en rotation dans le sens direct à la vitesse angulaire $\Omega$. Ainsi, toute déformation de la roue aubagée est vue, non pas par une combinaison linéaire de modes propres mais pas une combinaison linéaire d'ondes avance et rétrograde.
Le carter constituant lui aussi une structure axisymétrique, toutes les notions utilisées pour la mise en équations de la roue aubagée sont valables pour son étude. De la même façon donc, la vitesse de propagation de ces deux ondes pour un diamètre donné vaut $\pm\omega/n_d$ dans le repère fixe. Ces nouvelles définitions permettent de mieux caractériser le phénomène d'interaction.

## 3.2 Interaction modale

Les amplitudes de vibration d'une structure peuvent devenir importantes si celle-ci est excitée sur un de ses modes. Pour qu'il y ait interaction modale, il faut donc que :

- les deux structures acquièrent des déformées propices à un échange d'énergie, c'est-à-dire qu'elles vibrent toutes les deux sur un mode à même diamètre.

- chaque structure vibre à sa fréquence propre.

- les vitesses de propagation des ondes de déformation dans le repère fixe coïncident.

Physiquement, quatre solutions remplissent les conditions précédentes mais en considérant que la vitesse de rotation de la roue aubagée est importante face aux vitesses de propagation des ondes de déformation, ce sont les forces de frottement entre les deux structures, dirigées dans le sens inverse de rotation du rotor sur les aubes et dans le sens de rotation du rotor sur le carter, qui nous permettent de dire que seule la vitesse de rotation du moteur $\Omega$ vérifiant l'équation suivante est dangereuse vis-à-vis de l'interaction :

$$ \omega_c = n_d \Omega - \omega_a \quad (3) $$

avec :

- $n_d$: nombre de diamètres du mode considéré

- $\omega_a$: pulsation propre du mode à $n_d$ diamètres de la roue aubagée

- $\omega_c$: pulsation propre du mode à $n_d$ diamètres du carter

- $\Omega$: vitesse de rotation de la roue aubagée.

Cette relation analytique représente une condition nécessaire pour l'apparition d'une interaction modale.

# 4 Modélisation du contact et intégration temporelle

Le but de cette partie est de mettre au point une simulation numérique qui permette de vérifier les résultats théoriques de la section 3.

## 4.1 Algorithme de résolution

### 4.1.1 cas général

Le système matriciel découlant de la modélisation des structures est mis sous la forme :

$$ [M]\{\ddot{u}\} + [C]\{\dot{u}\} + [K]\{u\} = \{F_{NL}\} \quad (4) $$
où {$u$} est le vecteur déplacement comprenant les degrés de liberté du carter et de la roue aubagée et {$F_{NL}$}, celui des forces de frottements entre ces deux structures. La matrice {$C$} satisfait la condition de Basile. Pour gérer convenablement le contact entre les aubes et le carter, la méthode des lagrangiens est couplée à une méthode explicite d'intégration en temps basée sur les différences finies centrées qui permettent de discrétiser en temps le vecteur déplacement :

$$ egin{cases} 
\ddot{u}_n = \frac{u_{n+1} - 2u_n + u_{n-1}}{h^2} \\ 
\dot{u}_n = \frac{u_{n+1} - u_{n-1}}{2h} 
\end{cases} 
\qquad (5) $$

Il est démontré que la seule formulation consistante pour gérer le contact dans un algorithme d'intégration temporelle explicite est de type prédiction - correction [4]. La prédiction est une résolution des équations de la mécanique qui ne tient pas compte des éventuels contacts et la correction intervient à chaque fois qu'une pénétration entre les deux structures a été détectée.

Pour interdire toute pénétration, il est nécessaire d'appliquer le principe d'Hamilton à la fonctionnelle habituelle augmentée de l'énergie de contact [2]. Ensuite, la méthode de résolution est la suivante :

1. prédiction du vecteur déplacement {$u_{n+1}$} sans prendre en compte les éventuels contacts.

2. test sur le contact pour chaque aube.

3. si une pénétration est détectée, correction de {$u_{n+1}$}.

4. passage au pas de temps suivant.

### 4.1.2 cas particulier avec deux aubes en contact

Dans cette partie, on étudie un système composé de six aubes dont deux sont en contact avec le carter à un pas de temps donné. Ce choix permet de mettre en lumière la structure des équations et des simplifications qui seront imposées par la suite.

En faisant l'approximation $\theta_i = \Omega t + \alpha_i + (i - 1)\frac{\pi}{3}$ (ce qui est pertinent puisque l'excéntricité $e$ est très faible) et en posant $\beta_i = n_d \theta_i$, les pénétrations deviennent d'après la figure 3 :

$$ egin{cases} 
d_1 = l_c + a \cos \beta_1 + b \sin \beta_1 - l_a \left( 1 + \frac{e^2}{2l_a^2} - \frac{e}{l_a} \cos \theta_1 \right) \\ 
d_2 = l_c + a \cos \beta_2 + b \sin \beta_2 - l_a \left( 1 + \frac{e^2}{2l_a^2} - \frac{e}{l_a} \cos \theta_2 \right) 
\end{cases} 
\qquad (6) $$

où $a$ et $b$ sont les contributions des deux modes du carter. L'énergie de contact s'écrit alors :

$$ \prod_{\text{contact}} (d_1, d_2, \lambda_1, \lambda_2) = \lambda_1 d_1 + \lambda_2 d_2 \qquad (7) $$

$\lambda_1$ et $\lambda_2$ sont les multiplicateurs de Lagrange associés aux pénétrations $d_1$ et $d_2$. L'application du principe d'Hamilton permet ensuite d'écrire le système d'équations pour les phases de prédiction et de correction. La phase de prédiction est basée sur le schéma des différences
<figure><img src="image_4.png" /><figcaption>Figure 3 – Schéma de le pénétration aube/carter pour un excentrement e donné</figcaption></figure>

finies centrées et ne pose pas de difficulté particulière. La phase de correction donne les équations suivantes :

$$
\left\{
\begin{array}{l}
\displaystyle \frac{m_c}{h^2} (a^{n+1} - a^{n+1*}) + \lambda_1 \cos \beta_1^{n+1} + \lambda_2 \cos \beta_2^{n+1} = 0 \\[2ex]
\displaystyle \frac{m_c}{h^2} (b^{n+1} - b^{n+1*}) + \lambda_1 \sin \beta_1^{n+1} + \lambda_2 \sin \beta_2^{n+1} = 0 \\[2ex]
\displaystyle \frac{m_a}{h^2} (\alpha_1^{n+1} - \alpha_1^{n+1*}) + \lambda_1 (-e \sin \theta_1^{n+1} + n_d (b_{n+1} \cos \beta_1^{n+1} - a_{n+1} \sin \beta_1^{n+1})) = 0 \\[2ex]
\displaystyle \frac{m_a}{h^2} (\alpha_2^{n+1} - \alpha_2^{n+1*}) + \lambda_2 (-e \sin \theta_2^{n+1} + n_d (b_{n+1} \cos \beta_2^{n+1} - a_{n+1} \sin \beta_2^{n+1})) = 0 \\[2ex]
\displaystyle l_c + a_{n+1} \cos \beta_1^{n+1} + b_{n+1} \sin \beta_1^{n+1} - l_a \left(1 + \frac{e^2}{2l_a^2} - \frac{e}{l_a} \cos \theta_1^{n+1}\right) = 0 \\[2ex]
\displaystyle l_c + a_{n+1} \cos \beta_2^{n+1} + b_{n+1} \sin \beta_2^{n+1} - l_a \left(1 + \frac{e^2}{2l_a^2} - \frac{e}{l_a} \cos \theta_2^{n+1}\right) = 0
\end{array}
\right.
\tag{8}
$$

L'exposant $^{n+1}$ indique le pas de temps en cours et l'exposant $^{n+1*}$ indique que les variables ont été calculés pendant la prédiction : elles sont donc connues durant cette étape. La phase de correction revient à résoudre un système de six équations non linéaires. Plus le nombre d'aubes simultanément en contact est élevé, plus la taille du système précédent augmente. Cette procédure est donc très pénalisante en temps de calcul. La solution adoptée pour diminuer ce temps de calcul est, dans un premier temps, de considérer les termes soulignés dans le système précédent au pas de temps $n$ au lieu du pas de temps $n + 1$, et dans un deuxième
temps, de calculer les multiplicateurs de Lagrange indépendamment les uns des autres : pour chaque multiplicateur, la contribution de ses homologues sur la déformation du carter est alors prise nulle. Ensuite, ces multiplicateurs sont sommés pour calculer $a_{n+1}$ et $b_{n+1}$. Une comparaison avec la méthode de la pénalité permet de défendre ces simplifications. Même si les équations de la mécanique ne sont pas vérifiées exactement, il n'y a pas besoin de choisir un coefficient de pénalité *a priori* et les phénomènes de rebonds caractéristiques de cette méthode sont filtrés.

Cet algorithme a donc été testé sur un cas d'interaction modale.

## 4.2 Prédiction d'une interaction modale

Un premier calcul confirme que la roue aubagée doit elle aussi être réduite sur ses modes propres pour pouvoir simuler une interaction modale. La méthode utilisée se rapproche donc de celle décrite dans la section précédente mais les déformations des deux structures sont décomposées sur deux modes à même diamètre, choisis avant le calcul. Ce sont les modes à deux diamètres dans notre cas. Le contact entre les deux structures est assuré par l'intermédiaire d'un excentrement progressif de l'axe de rotation jusqu'à une valeur maximale.

Ce processus a été utilisé en faisant varier la vitesse de rotation de la roue aubagée et en ne traçant que l'amplitude maximale des vibrations du premier mode du carter. La figure 4 montre les résultats obtenus.

<figure><img src="image_5.png" /><figcaption>Figure 4 – Amplitude maximale des vibrations du carter en fonction de la vitesse de rotation du moteur (rad/s)</figcaption></figure>

Les deux cercles mettent en évidence les vitesses de rotation pour lesquelles l'amplitude des vibrations est devenue très importante. Le premier juste avant $\Omega = 150\,\text{rad/s}$ n'a pour l'instant pas d'explication. Le deuxième encadre la vitesse de rotation du moteur pour laquelle il peut y avoir, théoriquement, une interaction. En utilisant l'équation (3), cette vitesse de
rotation vaut 274.95rad/s. La méthode donne donc des résultats plutôt encourageants même si elle est très sensible à la valeur maximale de l'excentricité *e*. Si cette valeur est trop faible, il n'y a jamais interaction, si elle est trop importante, il y a interaction quelle que soit la vitesse de rotation du moteur.

La taille du pas de temps a été un problème important lors des simulations numériques. Plusieurs tests, pour des pas de temps allant de $10^{-5}s$ à $10^{-8}s$ ont montrés que la convergence numérique était difficile à atteindre. Un compromis a du être trouvé entre la précision et le temps de calcul, qui passe de quelques secondes pour un pas de temps de $10^{-5}s$ à plusieurs heures pour un pas de temps de $10^{-8}s$. Finalement, un pas de temps de $10^{-6}s$ a été choisi pour le calcul d'interaction modale précédent

# 5 Conclusion

Les premiers résultats obtenus au niveau de l'interaction modale dans un moteur d'avion sont encourageants. Même s'il est difficile de tirer des conclusions actuellement, la théorie existante et le modèle numérique très simple développé s'accordent : il faut que les deux structures vibrent sur un mode à même diamètre pour amorcer une interaction. Cependant, les temps de calcul sont assez importants pour une précision acceptable, le modèle numérique ne tient pas compte de l'abridable présent entre les aubes et le carter, le modèle retenu pour la roue aubagée est beaucoup trop rigide et enfin, l'excentrement du rotor doit être généralisé dans les équations de la mécanique. Le nombre de contacts simultanés entre les deux structures lors d'une simulation numérique rend difficile la convergence de l'algorithme de résolution vers la solution vraie. Il est donc envisagé de résoudre ce problème de transfert d'énergie lors du contact par une méthode multi-harmonique qui a donné de très bons résultats pour des cas de friction pure [3] et qui permettrait de réduire la résolution du système d'équations, pour un modèle de type éléments finis des structures, aux seuls degrés de liberté non linéaires, à savoir ceux impliqués dans le contact rotor/stator [5].

# Références

[1] BLADH, R. *Efficient predictions of the vibratory response of mistuned bladed disk by reduced order modeling*. Thèse de Doctorat, University of Michigan (2000).

[2] GIRAUDON-GUILLOTEAU, I. *Modélisation du contact en dynamique: construction d'un élément simplifié de contact et application à l'interaction rotor/stator*. Thèse de Doctorat, Ecole Centrale de Nantes (1999).

[3] GUILLEN, J. *Studies of the dynamics of dry-friction-damped blade assemblies*. Thèse de Doctorat, University of Michigan (1999).

[4] N. CARPENTER, R. TAYLOR ET M. KATONA. Lagrange constraints for transient finite element surface contact. *8th international seminar on elastomers, 32, 281-290 (1991).*

[5] S. NACIVET, C. PIERRE, F. THOUVEREZ ET L. JEZEQUEL. Analysis of periodic frictional contact in finite elements problems. *ASME Design Engineering Technical Conferences. Proceedings of DETC'01 - CDROM - DETC2001/VIB-21734 (2001).*