# Méthode de l'équilibrage harmonique généralisé

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

► To cite this version:

Sébastien Roques, Mathias Legrand, Christophe Pierre, Bernard Peseux, Patrice Cartraud. Méthode de l'équilibrage harmonique généralisé. 7e Colloque national en calcul des structures, May 2005, Giens, France. hal-01399609

HAL Id: hal-01399609

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

Submitted on 19 Nov 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
# Méthode de l'équilibrage harmonique généralisé

## Application à l'interaction modale rotor/stator

Sébastien Roques* — Mathias Legrand* — Christophe Pierre** — Bernard Peseux* — Patrice Cartraud*

* GeM, Pôle Structure et Couplage - École Centrale de Nantes
1 Rue de la Noë, 44321 Nantes cedex 3
{sebastien.roques,mathias.legrand,bernard.peseux,patrice.cartraud}@ec-nantes.fr

** Vibrations and Acoustics Laboratory - University of Michigan
3112 G.G. Brown Building, Ann Arbor MI 48109-2125, USA
pierre@umich.edu

RÉSUMÉ. La méthode de l'équilibrage harmonique permet l'analyse de problèmes fortement non-linéaires. Lorsque ces derniers sont oscillants mais apériodiques, il est possible de généraliser la formulation usuelle en utilisant le concept d'hyper-temps. Dans cet article, cette méthode est appliquée à un système impliquant une non-linéarité qui ne s'exprime pas de façon explicite. Les transformées directe et inverse de Fourier sont alors nécessaires pour alterner domaines temporel et fréquentiel. Les résultats montrent l'efficacité et la précision de la méthode.

ABSTRACT. The harmonic balance method is one of the methods for the analysis of strongly non-linear problems. When the latter are oscillating but aperiodic, it is possible to generalize the usual formulation by using the concept of hyper-time domain. In this article, the method is applied to a system with a non-linearity which cannot be expressed in the frequency domain. The direct and inverse fast Fourier transforms are then necessary to alternate between time and frequency domains. The results show the efficiency and the precision of the method.

MOTS-CLÉS : interaction modale, méthode de l'équilibrage harmonique généralisé, contact direct

KEYWORDS: modal interaction, generalized harmonic balance method, direct contact
## 1. Introduction

Dans le domaine des turbomachines comme les moteurs d'avion, les contacts structurels entre le carter et la roue aubagée peuvent avoir plusieurs origines : coïncidence vibratoire entre des modes, apparition d'un balourd au niveau du rotor, déformation thermique du carter, ... Ces interactions non linéaires peuvent fortement endommager les structures et il est important de comprendre l'origine de ces mécanismes. Cet article trouve son origine dans les travaux de M. Legrand concernant l'analyse de l'interaction modale rotor/stator [LEG 03]. Il s'agit ici d'étudier ce phénomène dans le domaine fréquentiel pour s'affranchir des déficits inhérents aux méthodes temporelles.

La nature apériodique du comportement de l'ensemble suggère une généralisation de la méthode usuelle de l'équilibrage harmonique traitant tout type de non-linéarité : c'est ce qui est présenté dans un premier temps. Afin de valider l'implantation numérique de la méthode, le cas simple d'une poutre en contact unilatéral est ensuite considéré. Enfin, on s'interesse aux modifications nécessaires pour traiter le cas des systèmes autonomes, notamment dans le cas de l'interaction modale rotor/stator.

## 2. Méthode de l'équilibrage harmonique à n dimensions

De manière générale, l'équation à résoudre pour un système non-linéaire s'écrit :

$$ M\ddot{X} + D\dot{X} + KX + F_{nl}(X, \dot{X}) = F_{ext} \quad [1] $$

où M, D, K représentent respectivement les matrices masse, amortissement et rai-deur; $F_{ext}$ et $F_{nl}$ les vecteurs force d'excitation et force non-linéaire. L'application de la méthode classique de l'équilibrage harmonique se limite à l'analyse de vibrations harmoniques mais se généralise à des systèmes apériodiques contenant $M_s$ fréquences sans dénominateur commun. Le lecteur pourra se référer à [LAU 83], à [EWI 01] ou [PUS 04] pour plus de détails.

Introduisons $M_s$ variables temporelles : $\tau_m = \omega_m t$ avec $m = 1, \dots, M_s$ appartenant au domaine hyper-temps $0 \le \tau_m \le 2\pi$ où $\omega_m$ sont les fréquences sans dénominateur commun. Le déplacement est discrétisé dans le domaine multi-fréquentiel comme suit :

$$ X = \sum_{j_1=-N_h}^{N_h} \sum_{j_2=-N_h}^{N_h} \cdots \sum_{j_m=-N_h}^{N_h} \left[ \mathbf{a}_{j_1, j_2, \dots, j_{M_s}} \cos \left( \sum_{m=1}^{M_s} j_m \tau_m \right) + \mathbf{b}_{j_1, j_2, \dots, j_{M_s}} \sin \left( \sum_{m=1}^{M_s} j_m \tau_m \right) \right] \quad [2] $$

En veillant à ne conserver que les combinaisons non redondantes grâce à la parité des fonctions trigonométriques, la composante $X_j$ du mouvement peut alors se mettre sous forme matricielle :

$$ X_j = T\mathbf{a}_j, \quad j = 1, 2, \dots, N \quad [3] $$
où $\mathbf{T} = [\mathbf{T}^{\cos}, \mathbf{T}^{\sin}]$, $\mathbf{a}_j = [\mathbf{a}_{\cos}^j, \mathbf{a}_{\sin}^j]^T$ avec $N_h = N_{\cos} + N_{\sin}$, le nombre total d'harmoniques conservés. Le vecteur ligne $\mathbf{T}$ est une combinaison de deux sous-vecteurs lignes. $\mathbf{T}^{\cos}$ se compose des termes en cosinus avec $j_m^i$ entiers :

$$ \mathbf{T}^{\cos} = \left[ \cos \left( \sum_{m=1}^{M_s} j_m^1 \tau_m \right), \cos \left( \sum_{m=1}^{M_s} j_m^2 \tau_m \right), \dots, \cos \left( \sum_{m=1}^{M_s} j_m^{N_{\cos}} \tau_m \right) \right] \quad [4] $$

et $\mathbf{T}^{\sin}$ a une expression similaire avec des termes en sinus. $\mathbf{a}_j$ est un vecteur colonne de coefficients harmoniques associés à la $j^e$ composante du vecteur déplacement et assemblé de façon cohérente avec $\mathbf{T}$ :

$$ \mathbf{a}_{\cos}^j = [a_{\cos}^1, a_{\cos}^2, \dots, a_{\cos}^{N_{\cos}}] \quad \text{et} \quad \mathbf{a}_{\sin}^j = [a_{\sin}^1, a_{\sin}^2, \dots, a_{\sin}^{N_{\sin}}] \quad [5] $$

En considérant les conventions suivantes avec le produit de Kronecker $\otimes$,

$$ \mathbf{Y} = \mathbf{T} \otimes \mathbf{I}_N \quad \text{et} \quad \mathbf{a} = [\mathbf{a}_1^T, \mathbf{a}_2^T, \dots, \mathbf{a}_N^T]^T \quad [6] $$

on construit les vecteurs déplacement $\mathbf{X}$ et vitesse $\dot{\mathbf{X}}$ par assemblage :

$$ \mathbf{X} = \mathbf{Y}\mathbf{a} \quad \text{et} \quad \dot{\mathbf{X}} = \left( \sum_{n=1}^{M_s} \omega_n \frac{\partial \mathbf{Y}}{\partial \tau_m} \right) \mathbf{a} \quad [7] $$

L'introduction des expressions [7] dans [1] suivie de la procédure de Galerkin sur le domaine de l'hyper-temps permet de réécrire les équations du mouvement dans le domaine fréquentiel :

$$ \mathbf{R} = \mathbf{H} \cdot \mathbf{a} + \mathbf{F}_{nl}^a(\mathbf{a}) - \mathbf{F}_{ext}^a \quad [8] $$

avec les notations suivantes :

$$ \mathbf{H} = \langle \mathbf{Y}, \left[ \sum_{m=1}^{M_s} \omega_m \left( \sum_{n=1}^{M_s} \omega_n \mathbf{M} \frac{\partial^2 \mathbf{Y}}{\partial \tau_m \partial \tau_n} + \mathbf{D} \frac{\partial \mathbf{Y}}{\partial \tau_m} \right) + \mathbf{K} \mathbf{Y} \right] \rangle \\ \mathbf{F}_{nl}^a(\mathbf{a}) = \langle \mathbf{Y}, \mathbf{F}_{nl}(\mathbf{X}, \dot{\mathbf{X}}) \rangle \\ \mathbf{F}_{ext}^a = \langle \mathbf{Y}, \mathbf{F}_{ext} \rangle $$

où l'opérateur $\langle \cdot, \cdot \rangle$ est défini comme suit :

$$ \langle \mathbf{M}_1, \mathbf{M}_2 \rangle = \int_0^{2\pi} \int_0^{2\pi} \dots \int_0^{2\pi} \mathbf{M}_1^T \mathbf{M}_2 \, d\tau_1 d\tau_2 \dots d\tau_M $$

L'annulation du résidu $\mathbf{R}$ avec un solveur non-linéaire de type quasi-Newton par exemple signifie que l'état d'équilibre est atteint.

### 3. Exemple d'application

Afin de valider cette nouvelle démarche pour une forte non-linéarité non explique dans le domaine fréquentiel, une poutre, discrétisée par quatre éléments finis de type Euler-Bernoulli et illustrée sur la figure 1, est soumise à une excitation externe apériodique composée de deux fréquences sans dénominateur commun d'expression : $f_{ext} = 100 \cdot (\cos(\pi t) + \cos(e^1 t))$. Le déplacement central est contraint par un mur.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schéma de la poutre contactant un mur</figcaption></figure>

En tenant compte des conditions aux limites en appui simple, le vecteur déplacement devient $\mathbf{X} = [\theta_1, v_2, \theta_2, v_3, \theta_3, v_4, \theta_4, \theta_5]^T$.

L'équation du mouvement se met sous la forme :

$$ \begin{cases} \mathbf{M}\ddot{\mathbf{X}} + \mathbf{D}\dot{\mathbf{X}} + \mathbf{K}\mathbf{X} = \mathbf{F}_{\text{ext}} \\ v_3 \le v_{\text{max}}, \quad \forall t \end{cases} \quad [9] $$

Pour ce cas précis, la prise en compte de la condition de contact dans la procédure fréquentielle est la suivante : si une violation est détectée dans le domaine temporel, on impose la valeur $v_{\text{max}}$ au déplacement du nœud 3 avant une projection dans le domaine fréquentiel.

Le résultat de la simulation est comparé à celui issu d'une intégration temporelle explicite fondée sur le schéma des différences finies centrées.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Comparaison entre la méthode fréquentielle pour $N_h = 12$ et la méthode temporelle</figcaption></figure>

Pour bien évaluer ce type de non-linéarité dans le domaine fréquentiel, il est nécessaire de considérer un grand nombre d'harmoniques. Par conséquent, douze harmoniques dans chaque direction sont conservés dans la simulation. Le temps de calcul de la méthode fréquentielle reste rédhibitoire à cause des transformées de Fourier directe et inverse en deux dimensions exécutées à chaque itération lors de la résolution non-linéaire. Néanmoins, ce calcul montre que la méthode fréquentielle est capable
de résoudre ce type de problème. Les deux solutions, temporelle et fréquentielle, sont quasiment indiscernables comme illustré sur les figures 2(a) et 2(b).

## 4. Phénomène d'interaction modale rotor/stator

L'intégration temporelle fait apparaître [LEG 04], sous certaines conditions portant notamment sur la vitesse de rotation du moteur, qu'un échange d'énergie, par l'intermédiaire de contacts structuraux, associé à un comportement oscillant apériodique s'instaure entre le carter et la roue aubagée même si toute sollicitation extérieure a disparu. Ce phénomène apparaît lorsque la condition de coïncidence vibratoire est satisfaite, condition qui sera retenue par la suite :

$$ \omega_{ra} = n_d \Omega - \omega_c \quad [10] $$

où $\omega_{ra}$, $\omega_c$, $\Omega$ et $n_d$ représentent respectivement les fréquences de réponse de la roue aubagée et du carter, la fréquence de rotation de la roue aubagée et le nombre de modes à diamètres considérés.

La figure 3 illustre le comportement de notre système pour $n_d = 3$. Le petit carré noir ■ suit artificiellement la propagation du mode tournant avance sur le carter alors que l'aube noire indique la progression du mode tournant rétrograde sur la roue aubagée. On y remarque que quand le mode avance du carter a effectué un tour complet (une durée de $n_d T_c$ s'est alors écoulée), l'aube noire n'est pas revenue à sa position initiale, ce qui résume le caractère apériodique de l'ensemble.

<figure><img src="image_4.png" /><figcaption>Figure 3. Évolution de l'ensemble roue aubagée - carter sur trois périodes $T_c$</figcaption></figure>

Dans l'optique de s'affranchir, à la fois de l'influence des conditions initiales et de la taille du pas de temps, lesquelles sont inhérentes à une résolution temporelle et particulièrement problématiques dans ce cas présent, une résolution multi-fréquentielle s'impose, ce qui nécessite des modifications étant donné que le système est autonome.

## 5. Cas des systèmes autonomes

Dans le cas de systèmes autonomes, la méthode de l'équilibrage harmonique généralisée évoquée précédemment doit être adaptée. En effet, la force d'excitation extérieure étant nulle, la solution nulle est triviale mais d'intérêt limité. Afin de trouver une
solution non nulle, il suffit d'inclure les fréquences dans les inconnues du problème en imposant les phases correspondantes à ces fréquences afin que la formulation soit bien posée.

Ceci revient à modifier l'équation [8] de la façon suivante :

$$ \mathbf{H} \cdot \Delta\mathbf{a} = \mathbf{H} \cdot \mathbf{a} + \mathbf{F}_{nl}^{a}(\mathbf{a}) - \sum_{m=1}^{M_s} \eta \Delta\omega_m \mathbf{Q}_m \quad [11] $$

où les notations sont celles utilisées pour l'équation [8] avec, en plus :

$$ \mathbf{Q}_m = \langle \mathbf{Y}, \left[ \sum_{n=1}^{M_s} 2\omega_n \mathbf{M} \frac{\partial^2 \mathbf{Y}}{\partial \tau_m \partial \tau_n} + \mathbf{D} \frac{\partial \mathbf{Y}}{\partial \tau_m} \right] \rangle \cdot \mathbf{a} \quad [12] $$

Le vecteur **a** est construit de la même façon qu'auparavant à ceci près que les composantes en première harmonique en sinus, en nombre égal aux fréquences sans dénominateur commun, sont imposées. Toutes les fréquences du système doivent être déterminées lors de la recherche d'une solution non nulle. La convergence est obtenue dès lors que le vecteur résidu composé des $\Delta\mathbf{a}$ et des $\Delta\omega_m$ s'annule.

## 6. Conclusions et perspectives

La méthode de l'équilibrage harmonique généralisé a été adaptée pour la résolution de systèmes fortement non-linéaires. En considérant le cas test de la poutre venant en contact avec un mur, les comparaisons de résultats entre une intégration temporelle et l'équilibrage harmonique généralisé sont satisfaisantes et montrent la précision de cette nouvelle méthode. Des travaux sont en cours de réalisation pour pouvoir mettre en œuvre cette méthode dans le cas de l'interaction modale rotor/stator.

## 7. Bibliographie

[EWI 01] EWINS D. J., VON GROLL G., « The harmonic balance method with arc-length continuation in rotor/stator contact problems », *Journal of Sound and Vibration*, vol. 241, n° 2, 2001, p. 223-233.

[LAU 83] LAU S. L., CHEUNG Y. K., WU S. Y., « Incremental harmonic balance method with multiple time scales for aperiodic vibration of nonlinear systems », *Journal of Applied Mechanics*, vol. 50, 1983, p. 871-876.

[LEG 03] LEGRAND M., PESEUX B., PIERRE C., « Amélioration de la prédiction de l'interaction rotor/stator dans un moteur d'avion », *Actes du 6e Colloque National en Calcul des Structures*, Giens, juin 2003, p. 53-68.

[LEG 04] LEGRAND M., PESEUX B., PIERRE C., « Étude de l'interaction modale rotor/stator dans un moteur d'avion », 14e *Colloque Vibrations, Chocs et Bruit*, Lyon, juin 2004.

[LEG 05] LEGRAND M., « Modèles de prédiction de l'interaction rotor/stator dans un moteur d'avion », Thèse de doctorat, École Centrale de Nantes, 2005.

[PUS 04] PUSENJAK R. R., OBLAK M. M., « Incremental harmonic balance method with multiple time variables for dynamical systems with cubic non-linearities », *International Journal for Numerical Methods in Engineering*, vol. 59, 2004, p. 255-292.