# Méthode Arlequin pour la modélisation multi-échelles en temps des machines tournantes

Assaf Ghanem, Mohamed Torkhani, Thouraya Nouri Baranger, Guillaume Altmeyer

## ► To cite this version:

Assaf Ghanem, Mohamed Torkhani, Thouraya Nouri Baranger, Guillaume Altmeyer. Méthode Arlequin pour la modélisation multi-échelles en temps des machines tournantes. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717113

HAL Id: hal-01717113

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

Submitted on 25 Feb 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
# Méthode Arlequin pour la modélisation multi-échelles en temps des machines tournantes

Assaf Ghanem<fnref n="1" />, Mohamed Torkhani<fnref n="1" />, Thouraya Nouri Baranger²³, Guillaume Altmeyer<fnref n="3" />

**Résumé** — L'extension du cadre Arlequin à la dynamique des machines tournantes nécessite une réflexion théorique particulière. En effet, pour un raccord 1D-3D, les modèles de poutre écrits par défaut dans un repère Galiléen fixe sont couplés à des modélisations tridimensionnelles pouvant être formulées dans le repère tournant. Dans ce papier nous étudions ce couplage multi-repères dans un contexte multi-schémas/multi-échelles en temps dans le cadre Arlequin. Ces travaux sont illustrés par des simulations numériques menées sur un exemple de contact entre un rotor et un stator d'une machine tournante.

**Mots clés** — Dynamique transitoire, analyse multi-échelles/multi-repères, couplage Arlequin, machines tournantes.

## 1 Introduction

La simulation de type machines tournantes en régime dynamique est généralement basée sur des modélisations de type poutre (structures filaires et élancées). La prise en compte d'un phénomène localisé (de type fissure, contact, impact) en se basant sur ces modèles globaux de structures minces, n'est pas suffisante, et une analyse sous l'angle du modèle tridimensionnel qui les sous-tend est alors nécessaire. Par ailleurs, l'étude de ces phénomènes locaux dans un régime dynamique peut nécessiter la prise en compte d'une échelle temporelle représentative. Ce papier vise à mettre en place les ingrédients d'un cadre rigoureux, permettant d'effectuer des études de machines tournantes qui combinent la rapidité des simulations avec la précision et la fiabilité des résultats. En se basant sur les travaux en multi-échelles en temps présentés dans [1, 2, 4], cette méthode permet de traiter simultanément un phénomène local formulé dans le repère tournant [5] sur un modèle 3D fin en temps et en espace, ainsi que le comportement global de la structure 1D contenant ce défaut local sur des échelles grossières en temps et en espace. La mise en place de cette approche nécessite des développements théoriques spécifiques. Les points-clés suivants ont été pris en compte lors de la mise en équation du problème à résoudre :

- Les équations du mouvement de chacun des modèles 1D et 3D doivent être écrites dans les repères fixe et tournant, respectivement.

- Les multiplicateurs de Lagrange $\Lambda$ sont toujours traités dans le repère fixe car ils sont définis sur le même espace que celui du modèle 1D.

- La condition de collage Arlequin est écrite dans le repère fixe.

## 2 Formulation Arlequin en dynamique de rotation

On considère une ligne d'arbre horizontale en rotation autour de l'axe $(\vec{X})$. Elle est modélisée avec un modèle de poutre de Timoshenko à six degrés de liberté par noeud, de section circulaire constante, et en appui simple sur ses deux extrémités. Un disque est placé aux 2/3 de sa longueur et un patch 3D constitué d'éléments volumiques à 20 noeuds représentant une partie volumique du rotor est introduit au 1/3 de la longueur. La partie 3D est constituée de deux zones de collage (64 éléments) et d'une zone de recouvrement (64 éléments). Chaque élément de poutre est couplé avec 16 éléments 3D. Le couplage 3D-1D est présenté dans la figure 1. La discrétisation en temps suit les différentes propriétés associés à chacun des modèles couplés :

<footnote n="1">¹LaMSID, UMR CNRS-EDF 2832, France, (Assaf.Ghanem, Mohamed.Torkhani)@edf.fr</footnote>

<footnote n="2">²Université de Lyon, CNRS, LaMCoS UMR5259, France, (Thouraya.Baranger,Guillaume.Altmeyer)@insa-lyon.fr</footnote>

<footnote n="3">³ Université Lyon1, France, (Thouraya.Baranger@univ-lyon1.fr)</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1 – Couplage Arlequin 3D-1D</figcaption></figure>

* Modèles 3D : échelle temporelle fine ($\Delta t$), repère tournant $\mathbf{R}_r(x_r, y_r, z_r)$, schéma de Newmark explicite ou implicite (selon l'application).

* Modèles 1D : échelle temporelle grossière ($\Delta T$), repère fixe (galiléen) $\mathbf{R}_f(x_f, y_f, z_f)$, schéma de Newmark implicite (conservatif ou bien dissipatif).

Sur le domaine $\Omega_{1d}$, on considère un pas de temps grossier $\Delta T$ qui varie de $t_0$ jusqu'à $t_m$ tel que $t_m = t_0 + \Delta T$ comme le montre la figure 2 et on lui associé les paramètres de Newmark $\gamma_{1d}$ et $\beta_{1d}$. De façon similaire, la discrétisation sur $\Omega_{3d}$ est faite aux instants $t_j \forall j \in 1, 2, \dots, m$ avec les paramètres de Newmark $\gamma_{3d}$ et $\beta_{3d}$.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Discrétisation temporelle sur les repères fixe $\mathbf{R}_f$ et rotationnel $\mathbf{R}_r$</figcaption></figure>

La formulation discrétisée en espace et en temps du problème bi-modèle Arlequin présenté sur la figure 1, s'écrit de la manière suivante :

$$ \mathbf{M}_{1d} \ddot{\mathbf{u}}_m^{1d} + \mathbf{C}_{1d} \dot{\mathbf{u}}_m^{1d} + \mathbf{K}_{1d} \mathbf{u}_m^{1d} + \mathbf{L}_{1d}^t \boldsymbol{\lambda}_m = \mathbf{f}_m^{1d}, \quad (1a) $$

$$ \mathbf{M}_{3d} \ddot{\mathbf{u}}_j^{3d} + \mathbf{C}_{3d} \dot{\mathbf{u}}_j^{3d} + \mathbf{K}_{3d} \mathbf{u}_j^{3d} + \mathbf{L}_{3d}^t \boldsymbol{\lambda}_j^r = \mathbf{f}_j^{3d}, \quad (1b) $$

$$ \mathbf{L}_{1d} \mathbf{u}_m^{1d} + \mathbf{\Phi}_m \mathbf{L}_{3d} \mathbf{u}_m^{3d} = 0 \quad (1c) $$

L'amortissement global $\mathbf{C}_i$ est constitué du cumul des matrices d'amortissement structural et de gyroscopie. Les matrices de raideur globales $\mathbf{K}_i$ cumulent, en plus de la raideur classique de la structure, les contributions des termes de raideur centrifuge et de précontrainte. Les matrices $\mathbf{L}_{1d}$ et $\mathbf{L}_{3d}$ représentent le couplage via les zones de collages.

Les multiplicateurs de Lagrange $\lambda_0^r$ et $\lambda_m^r$ liés au repère en rotation doivent respecter les conditions suivantes :

$$ \lambda_0^r = \mathbf{\Phi}_0^t \boldsymbol{\lambda}_0, \quad \lambda_m^r = \mathbf{\Phi}_m^t \boldsymbol{\lambda}_m. \quad (2) $$

La méthode multi-échelles/multi-schemas en temps développée dans [6] dans le but de prendre en compte le changement de repères tout en préservant la conservation de l'énergie globale [1, 2], permet d'exprimer $\lambda_j^r$ par une interpolation linéaire entre $\mathbf{\Phi}_0^t \boldsymbol{\lambda}_0$ et $\mathbf{\Phi}_m^t \boldsymbol{\lambda}_m$ en fonction d'un paramètre $\alpha_j$ :

$$ \lambda_j^r = (1 - \alpha_j) \mathbf{\Phi}_0^t \boldsymbol{\lambda}_0 + \alpha_j \mathbf{\Phi}_m^t \boldsymbol{\lambda}_m \quad \forall j \in [0, m] \quad \text{où } \alpha_j = \tilde{\gamma} (1 - a^j) \quad (3) $$

L'équation (3) représente l'interpolation, sur l'échelle fine en temps $t_j$, de la projection des multiplicateurs de Lagrange du repère fixe au repère tournant. La construction du paramètre de Newmark $\tilde{\gamma}$ dans la zone de collage $\Omega_0$ se fait comme suit :

$$ \tilde{\gamma} = \frac{1}{1 - a^m} \quad \text{avec } a^m \neq 1 \text{ et } a = \frac{\gamma_{3d} - 1}{\gamma_{3d}} \quad (4) $$

Cette pondération des multiplicateurs de Lagrange sur l'échelle fine en tenant compte de la rotation des repères formulée dans (3), constitue un cadre général qui couvre l'approche en dynamique sans
$$ \gamma = \gamma_{1d} \text{ sur } \Omega_{1d} \setminus \Omega_0 \qquad ; \qquad \beta = \beta_{1d} \text{ sur } \Omega_{1d} \setminus \Omega_0 $$

$$ \gamma = \gamma_{3d} \text{ sur } \Omega_{3d} \setminus \Omega_0 \qquad ; \qquad \beta = \beta_{3d} \text{ sur } \Omega_{3d} \setminus \Omega_0 $$

$$ \gamma = \tilde{\gamma} \text{ sur } \Omega_0 \qquad ; \qquad \beta = \tilde{\beta} \text{ sur } \Omega_0 $$

avec

$$ \tilde{\gamma} = \frac{\gamma_{3d}^n}{\gamma_{3d}^n - (1 - \gamma_{3d})^m} \text{ sur } \Omega_0 \text{ pour ; } \gamma_{3d} > 0.5 \text{ et } m \ge 1. \quad ; \quad \tilde{\beta} = 0.25(0.5 + \tilde{\gamma})^2 \text{ sur } \Omega_0. $$

<caption>Tableau 1 – Distribution des paramètres de Newmark - Couplage Arlequin multi-schémas/multi-échelles en temps</caption>

rotation d'ensemble. En effet, pour des angles de rotation nuls, les matrices $\Phi_0$ et $\Phi_m$ sont équivalentes à une matrice identité, ce qui nous ramène du cadre général (3) à un cadre plus particulier [6, 1]. La distribution des paramètres de Newmark sur les différentes zones du domaine est résumée dans le tableau 1. L'algorithme de résolution est décrit en détail dans les travaux présentés dans [6, 1].

## 3 Application au machines tournantes

Le modèle sur la figure 1 est soumis à un effet de balourd au niveau du disque. On calcule la réponse au balourd de la ligne d'arbre sur base physique et on surveille l'écart en déplacement entre le rotor et le stator. Les données caractérisant la poutre sont présentées dans le tableau 3.

<table><caption>Tableau 2 – Caractéristiques du rotor</caption><thead><tr><th>Rayon</th><th>Longueur</th><th>Module de Young (N.m<sup>2</sup>)</th><th>Masse volumique (kg.m<sup>3</sup>)</th><th>Coefficient de Poisson</th></tr></thead><tbody><tr><td>r = 0.01 m</td><td>L = 0.4 m</td><td>E = 2.10<sup>11</sup></td><td>ρ = 7800</td><td>ν = 0.3</td></tr></tbody></table>

L'application du contact dans le cadre multi-repères (figure 3) nécessite des modifications et des développements spécifiques supplémentaires dans l'algorithme de résolution. Il existe une variété de

<figure><img src="image_4.png" /><figcaption>Fig. 3 – Illustration du contact rotor/stator</figcaption></figure>

méthodes numériques permettant de gérer le contact entre deux solides dont la méthode de pénalisation, la méthode des multiplicateurs de Lagrange et les méthodes hybrides. Dans le cadre de cette application, nous avons fait le choix de mettre en place la méthode de pénalisation. C'est une méthode relativement simple à implémenter par rapport aux autres approches.
La validation numérique et la pertinence de la méthode est mise en valeur à travers deux applications en régime dynamique non-linéaire. Cette étude concerne le cas d'un contact entre le patch 3D représentant une partie de la ligne d'arbre du rotor et un stator. Cette étude est réalisée dans le cadre du couplage multi-schémas / multi-échelles en temps / multi-repères, et permet de valider ce formalisme global dans le cadre d'une application non linéaire pour une vitesse de rotation constante de 3000 tr/min. La deuxième étude concerne le comportement du rotor lors d'une montée en vitesse de type linéaire, suivie d'une stabilisation à 5000 tr/min. Dans les deux cas de figure nous nous concentrons, pour le traitement du contact, sur la couronne de noeuds présents au centre de la zone libre (figure 4). On note que le patch 3D n'est pas très fin en espace. Ce choix est induit par des contraintes de temps de calcul et de mémoire CPU disponible pour une maquette Matlab. Ceci justifie et motive le développement de cette approche dans un code industriel tel que code-ASTER.

<figure><img src="image_5.png" /><figcaption>Fig. 4 – Zoom sur le patch 3D</figcaption></figure>

## 3.1 Vitesse nominale

Cette première application s'effectue pour un ratio $m = 15$ entre les échelles de temps, et pour un couplage d'un schéma de Newmark dissipatif sur le modèle 3D ($\gamma_{3d} = 0.7$ ; $\beta_{3d} = 0.36$), avec un autre conservatif d'accélération moyenne sur la modélisation poutre ($\gamma_{1d} = 0.5$ ; $\beta_{1d} = 0.25$).

La distribution générale des paramètres du schéma de Newmark sur la poutre, le patch 3D, et la zone de collage est décrite dans le tableau 1. Les figures 5 et 6 montrent respectivement les déplacements dans le plan selon les axes $(\vec{Y})$ et $(\vec{Z})$ d'un noeud de la couronne de contact (figure 4) en fonction du temps. Dans ce cas de figure, nous avons fait le choix de tenir compte de la présence du contact à partir de l'instant $t_0 = 0.25$ s dans le but de se placer dans le domaine stabilisé. Dans ces deux graphiques nous représentons les déplacements nodaux pour la configuration sans contact, ainsi que ceux pour la configuration qui tient compte du contact. Nous remarquons qu'à partir de $t_0$ le contact est détecté correctement et que les déplacements respectifs selon les axes $(\vec{Y})$ et $(\vec{Z})$ tiennent compte de la présence d'un seuil de contact. La figure 7 montre l'orbite du même noeud de la couronne utilisé pour tracer les figures 5 et 6. Dans ce graphe on distingue trois courbes : l'orbite du nœud dans une configuration qui ne tient pas compte du contact, la même orbite corrigée en présence du contact, ainsi que l'illustration du stator. Nous pouvons distinguer la différence entre la trajectoire sans contact et celle qui tient compte de la correction des efforts de contact. Il est par ailleurs possible de remarquer l'interpénétration autorisée entre la trajectoire des nœuds de contact et la position du stator. Cette application à vitesse stabilisée constitue une première validation de l'approche globale multi-échelles / multi-schémas / multi-repères.
<figure><img src="image_6.png" /><figcaption>Fig. 5 – Déplacement nodal 3D selon l'axe Y avec/sans contact</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6 – Déplacement nodal 3D selon l'axe Z avec/sans contact</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7 – Comparaison de la trajectoire d'un noeud avec et sans contact</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8 – Couplage Arlequin 3D-1D - Modèle avec palier</figcaption></figure>

## 3.2 Transitoire de vitesse

Il s'agit de la même modélisation décrite dans la section 3.1 en ajoutant un palier avec amortissement visqueux aux environs du patch 3D du rotor (Figure 8). Les coefficients du palier sont définis comme suit :

$$ K_{yy} = 2.10^5 N.m^{-1} ; K_{zz} = 5.10^5 N.m^{-1} ; K_{zy} = K_{yz} = 0 \quad (5) $$

$$ C_{yy} = 4.10^1 N/m.s^{-1} ; C_{zz} = 1.10^2 N/m.s^{-1} ; C_{zy} = C_{yz} = 0 \quad (6) $$

Le couplage multi-échelles en temps s'effectue pour $\Delta T = 20\Delta t$. Le patch 3D est intégré avec un schéma de Newmark dissipatif ($\gamma_{3d} = 0.6$ ; $\beta_{3d} = 0.3036$), et la modélisation poutre avec un schéma de Newmark conservatif d'accélération moyenne ($\gamma_{1d} = 0.5$ ; $\beta_{1d} = 0.25$).

<figure><img src="image_10.png" /><figcaption>Fig. 9 – $w(t) = 5000t$</figcaption></figure>

Nous nous intéressons à présent à la détection du contact. Nous déterminons le seuil de contact à 25% de l'amplitude maximale du déplacement dans la configuration sans contact soit à $0,95 \cdot 10^{-5}$ m. On note que pour un instant $t_0$ où un contact est détecté, la prise en compte de l'effort de contact correspondant s'effectue à l'instant grossier suivant $t_m = t_0 + \Delta T$. Cette correction à posteriori permet de reproduire et d'analyser le comportement de la ligne d'arbre quand des phénomènes non-linéaires similaires prennent lieu. La figure 9 montre l'amplitude de la montée linaire en vitesse de la ligne d'arbre pendant 1 seconde de 0 tr/min à 5000 tr/min avant de stabiliser à 5000 tr/min :

$$ w(t) = 5000t \quad (7) $$

Le post-traitement s'effectue au niveau des noeuds 3D situés sur la couronne définie dans la figure 4. Nous pouvons constater qu'à l'instant $t = 0.78$ s un premier contact est établi. À ce niveau, plusieurs remarques peuvent être faites. En effet, en comparant la solution observée au niveau du patch 3D avec la solution de référence de poutre au même noeud il est possible de constater une diminution de l'amplitude globale des déplacements dans le plan (YOZ) où la rotation a lieu. Cette correction entraîne une perte de contact aux instants suivants. Ensuite, l'effet de rotation provoque un deuxième contact à l'instant $t = 0.814$ s. Ce dernier est aussi pris en compte et corrigé. Une fois les amplitudes des déplacements ne dépassent plus le seuil de contact, il est possible de remarquer que le comportement général du bi-modèle Arlequin suit l'allure de la réponse du modèle sans contact. Les figures 10, 11 montrent respectivement les déplacements selon les axes $(\vec{Y})$, $(\vec{Z})$. Le jeu de contact défini à 25% de l'amplitude maximale est à l'origine des décalages observés par rapport à la courbe de référence et montre la prise en compte du contact. Nous pouvons par ailleurs remarquer l'effet du premier contact à $t = 0.78$ s au niveau des
<figure><img src="image_11.png" /><figcaption>Fig. 10 – $\tilde{w}_1 = 5000t$ - Déplacement selon l'axe Y avec/sans contact</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11 – $\tilde{w}_1 = 5000t$ - Déplacement selon l'axe Z avec/sans contact</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12 – $\tilde{w}_1 = 5000t$ - Orbite d'un noeud de contact du patch 3D</figcaption></figure>
déplacements en (Y) dans figure 10, et celui du second contact en (Z) à $t = 0.814$ s dans la figure 11.
Ces mêmes constations sont visibles à travers la figure 12 qui trace la trajectoire au niveau d'un noeud
de contact 3D. Ce graphique permet de mieux visualiser la modification de l'orbite causée par le contact
entre le patch 3D et le stator.

## Conclusion

Nous avons présenté le formalisme général autorisant un couplage multi-modèles / multi-échelles / multi-schémas dans le cadre Arlequin pour les machines tournantes. À partir d'un cas test numérique représentatif d'une application industrielle nous avons montré l'efficacité de cette approche. Par l'intermédiaire de la non-linéarité de contact nous avons montré la faisabilité et l'utilité d'une telle approche dans le cadre d'une analyse avancée de machines tournantes. En se basant sur ces résultats encourageants, il est important de viser l'implémentation de cette approche dans un code industriel dans l'optique de réaliser des analyses avancées de machines tournantes telles que l'étude de la propagation de fissures en 3D et d'utiliser des algorithmes plus performants pour le traitement du contact.

# 4 Références bibliographiques

## Références

[1] A. Ghanem, M. Torkhani, N. Mahjoubi., T.N. Baranger, A. Combescure, *Arlequin framework for multi-model, multi-time scale and heterogeneous time integrators for structural transient dynamics*, Computer Methods in Applied Mechanics and Engineering, 254, 292-308, 2013.

[2] A. Combescure, A. Gravouil, *A numerical scheme to couple subdomains with different time steps for predominantly linear transient analysis*, Computer Methods in Applied Mechanics and Engineering, 191, 1129-1157, 2002.

[3] A. Ghanem, T.N. Baranger, N. Mahjoubi, M. Torkhani, *Arlequin framework for structural dynamics*, Proceedings of the 3rd International Conference on Computational Methods in Structural Dynamics & Earthquake Engineering (COMPDYN 2011), 26-28 May, Corfu, Greece.

[4] H.B. Dhia, G. Rateau, *The Arlequin method as a flexible engineering design tool*, International Journal for Numerical Methods in Engineering, 62, 1442-1462, 2005.

[5] T. GMUR, *Dynamique des structures : Analyse modale numérique*, Presses Polytechniques et Universitaires Romandes (PPUR), 2009.

[6] A. Ghanem, *Contributions à la modélisation avancée des machines tournantes en dynamique transitoire dans le cadre Arlequin*, Thèse, INSA-Lyon, 2013.