# Jumeau numérique d'un turbocompresseur à paliers à gaz pour la détection d'anomalies

Dimitri Goutaudier, Jürg Schiffmann, Fabio Nobile

► To cite this version:

Dimitri Goutaudier, Jürg Schiffmann, Fabio Nobile. Jumeau numérique d'un turbocompresseur à paliers à gaz pour la détection d'anomalies. 15e colloque national en calcul des structures, Université Polytechnique Hauts-de-France [UPHF], May 2022, Hyères-les-Palmiers, France. hal-03719273

HAL Id: hal-03719273

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

Submitted on 11 Jul 2022

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
# Jumeau numérique d'un turbocompresseur à paliers à gaz pour la détection d'anomalies

D. Goutaudier<fnref n="1" />,<fnref n="2" />, J. Schiffmann<fnref n="1" />, F. Nobile<fnref n="2" />

**Résumé** — Les turbocompresseurs à paliers à gaz présentent une efficacité énergétique prometteuse pour de nombreuses applications en chauffage domestique ou industriel. Le couplage non-linéaire des différents composants rend néanmoins leur surveillance opérationnelle difficile. L'objectif de ce travail est de développer un jumeau numérique pour la détection d'anomalies en utilisant des méthodes de réduction de modèle enrichies par les données. L'application visée est l'identification en temps-réel de précurseurs d'instabilités de rotation pouvant conduire à la rupture.

**Mots clefs** — jumeau numérique, modèle réduit, enrichissement par les données, détection d'anomalies

## 1. Introduction

La présence d'huile dans les circuits des pompes à chaleur classiques complique l'amélioration de leur efficacité énergétique (joints d'étanchéité, viscosité de l'huile, etc.). Pour répondre à cette difficulté, les turbocompresseurs à paliers à gaz ont été identifiés comme étant une technologie prometteuse, notamment dans les domaines du chauffage domestique et industriel [1]. Le principal lubrifiant utilisé dans les paliers à gaz est l'air qui est cent à mille fois moins visqueux que les fluides lubrifiants usuels. L'avantage est que la puissance dissipée par frottements est ainsi considérablement diminuée. En revanche, la couche d'air doit être très mince (< 30µm) et la vitesse de rotation très élevée (100-200krpm) pour supporter le rotor par effort aérodynamique [2]. Ainsi cette technologie convient particulièrement pour des turbocompresseurs de petite taille avec un diamètre de rotor de 10-30mm.

Cependant, les couplages entre les différents composants d'un turbocompresseur (voir Figure 1) agissent à différentes échelles de temps et plusieurs anomalies de fonctionnement résultant de ces couplages peuvent survenir. En particulier, le couplage fluide-structure entre le rotor et les couches fluides dans les paliers peut générer des orbites instables du rotor [3]. Une rupture catastrophique se produit au moindre contact du rotor avec l'un des paliers du fait de la vitesse de rotation très élevée. L'objectif du projet présenté dans cette contribution est de développer le jumeau numérique d'un turbocompresseur à paliers à gaz de petite taille pour détecter de telles anomalies en temps réel. Des lois de contrôle correctives pourraient alors être appliquées à temps pour les éviter.

Ce document est organisé comme suit. Dans la section 2 nous décrivons la modélisation physique du couplage fluide-structure entre le rotor et la couche fluide dans les paliers à gaz. Dans la section 3 nous présentons les techniques de réduction de modèle utilisées pour diminuer les temps d'évaluation du modèle physique discrétisé en espace. Dans la section 4 nous présentons des techniques d'enrichissement par les données que nous allons explorer dans la suite de ce projet pour améliorer la capacité prédictive du modèle en conditions opérationnelles. Enfin, dans la section 5 nous présentons brièvement un exemple d'utilisation de ce modèle hybride qui sera étudié dans ce projet.

<footnote n="1">¹LAMD, Ecole Polytechnique Fédérale de Lausanne, Suisse, {dimitri.goutaudier,jurg.schiffmann}@epfl.ch</footnote>

<footnote n="2">²CSQI, Institut de Mathématiques, Ecole Polytechnique Fédérale de Lausanne, Suisse, fabio.nobile@epfl.ch</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Composants principaux d'un turbocompresseur de petite taille à paliers à gaz (diamètre du rotor : 10-20mm, vitesse de rotation : 100-200krpm)</figcaption></figure>

## 2. Modélisation physique du système dynamique

Les instabilités de rotation résultent d'un couplage fluide-structure entre deux dynamiques, celle du rotor et celle de la couche fluide dans les paliers à gaz. Des vibrations indésirables du système couplé peuvent en effet s'activer dans certaines conditions, augmentant dangereusement l'amplitude des orbites du rotor [3]. Etant donnée la vitesse de rotation très élevée de ce dernier (100-200krpm), une rupture catastrophique survient au moindre contact entre le rotor et un des paliers.

### 2.1. Dynamique du rotor

Le turbocompresseur est dimensionné de sorte que la première fréquence propre du rotor ne soit pas excitée par les conditions nominales de fonctionnement [4]. Le rotor est alors en très bonne approximation modélisé par un arbre rigide en rotation (Figure 2) supporté par la force aérodynamique correspondant à la résultante de pression appliquée par la couche de fluide dans les paliers :

$$ \begin{bmatrix} F_x \\ F_y \end{bmatrix} = -2 \iint (P-1) \begin{bmatrix} \cos \theta \\ \sin \theta \end{bmatrix} d\theta dz \quad (1) $$

La rotation axiale étant imposée par un moteur et la translation axiale étant bloquée par des butées, quatre degrés de liberté suffisent à décrire ce solide rigide, par exemple en choisissant les excentricités $e_1 = [e_{x1}, e_{y1}]$ et $e_2 = [e_{x2}, e_{y2}]$ définies par les intersections de l'axe du rotor avec les plans médians de chaque palier. En discrétisant la double intégrale dans (1) avec une règle de quadrature, et en notant $y_1 = [e_1, e_2, \dot{e}_1, \dot{e}_2]$, on obtient une représentation en variables d'état de la dynamique du rotor de la forme suivante :

$$ \frac{dy_1}{dt} = A_1 y_1 + B_1 u_1 + C_1(y_1)y_2 \quad (2) $$

où les matrices $A_1$ et $B_1$ sont indépendantes du temps, $C_1$ est non-linéaire en $y_1$, et où $y_2$ représente des variables d'état de la couche fluide.
<figure><img src="image_3.png" /><figcaption>Figure 2: Configuration géométrique du système rotor-paliers [3]</figcaption></figure>

## 2.2. Dynamique de la couche fluide

La Figure 3 illustre la configuration géométrique d'un rotor en mouvement dans un palier à gaz. La vitesse de rotation du rotor est supposée constante dans tout ce qui suit. La dynamique de la couche fluide située entre chaque palier et le rotor est gouvernée par l'équation de Reynolds, obtenue à partir des équations de Navier-Stokes en faisant les trois hypothèses suivantes : (i) la couche fluide est mince, (ii) l'écoulement est laminaire, et (iii) le fluide est newtonien. En coordonnées cylindriques, cette équation est donnée par :

$$ \partial_{\theta}[h^3 P \partial_{\theta} P] + \partial_{z}[h^3 P \partial_{z} P] = \Lambda \partial_{\theta}[Ph] + \sigma \partial_{t}[Ph] \quad (3) $$

où $P(\theta, z, t)$ et $h(\theta, z, t)$ sont respectivement les distributions adimensionnées de pression du fluide et de jeu entre le rotor et le palier. Cette équation aux dérivées partielles non-linéaire est usuellement discrétisée en espace avec un schéma de différences finies sur une grille rectangulaire [5]. Pour cela, on réalise le changement de variable $\psi = Ph$ faisant disparaître la dépendance cubique en $h$, puis on désigne par une lettre en gras les vecteurs formés des variables $\psi$ et $h$ (et de leurs dérivées partielles) évaluées aux points de discrétisation $(\theta_i, z_j)$. En notant $\cdot$ le produit composante par composante entre deux vecteurs, il vient :

$$ \begin{aligned} & h \cdot [\partial_{\theta} \psi \cdot \partial_{\theta} \psi + \partial_{z} \psi \cdot \partial_{z} \psi] + h \cdot \psi \cdot [\partial_{\theta}^{2} \psi + \partial_{z}^{2} \psi] \\ & - \psi \cdot [\partial_{\theta} h \cdot \partial_{\theta} \psi + \partial_{z} h \cdot \partial_{z} \psi] - \psi \cdot \psi \cdot [\partial_{\theta}^{2} h + \partial_{z}^{2} h] \\ & = \Lambda \partial_{\theta} [\psi] + \sigma \partial_{t} [\psi] \end{aligned} \quad (4) $$

Lorsque des solutions transitoires sont recherchées, la distribution de jeu entre le rotor et les paliers dépend du mouvement du rotor par une relation de la forme suivante :

$$ h(\theta, z, t) = C[1 - l(\theta, z, t) \cdot y_1(t)] \quad (5) $$

où $y_1$ représente les variables d'état du rotor (voir section 2.1) et où $l$ désigne le produit scalaire usuel. En discrétisant les opérateurs de dérivation partielle dans (4) avec un schéma de différences finies, en tenant compte des conditions aux frontières, et en notant $y_2 = \psi$, on obtient en utilisant (5) une représentation en variables d'état du modèle dynamique de la couche fluide de la forme suivante :

$$ \frac{dy_2}{dt} = A_2 y_2 + f_2(y_2) + C_2(y_2)y_1 \quad (6) $$
où la matrice $A_2$ est indépendante du temps et où $f_2$ et $c_2$ sont non-linéaires en $y_2$. La dynamique de la couche fluide est ainsi couplée avec celle du rotor via le dernier terme dans le membre droit.

<figure><img src="image_4.png" /><figcaption>Figure 3: Configuration géométrique du rotor dans un palier [6]</figcaption></figure>

## 2.3. Système couplé

La dynamique du système rotor-paliers est gouvernée par les équations (2) et (6) en prescrivant des conditions initiales. En posant $\mathbf{y} = [\mathbf{y}_1, \mathbf{y}_2]$ on obtient une représentation en variables d'état du système couplé de la forme générale :

$$ \left\{ \begin{array}{l} \frac{dy}{dt} = g(\mathbf{y}) \\ \mathbf{y}(0) = \mathbf{y}_0 \end{array} \right. \qquad (7) $$

Cette équation différentielle ordinaire en temps est usuellement intégrée avec un solveur implicite itératif utilisant un pas de temps adaptatif [7]. La Figure 4 reproduit des trajectoires du centre d'un rotor obtenues dans [5] en résolvant un système couplé similaire.

<figure><img src="image_5.png" /><figcaption>Figure 4: Exemples de trajectoires du centre d'un rotor avec une discrétisation spatiale par différence finies (courbes en noir) pour deux vitesses de rotation du rotor [5]</figcaption></figure>
## 3. Réduction du modèle physique

La résolution numérique du modèle physique (7), c'est-à-dire l'évaluation de l'état $\mathbf{y}$ du système couplé sur un intervalle de temps et un jeu de paramètres donnés, nécessite un temps de calcul trop élevé avec les ressources de calcul que nous souhaitons utiliser pour l'application visée. L'expérience montre en effet que des instabilités de rotation apparaissent typiquement quelques secondes avant la rupture. Il est donc nécessaire de considérablement réduire les temps d'évaluation pour identifier à temps d'éventuels précurseurs d'instabilité. Pour cela, nous proposons en première étape du projet d'utiliser des techniques de réduction de modèle par projection dans une base réduite.

### 3.1. Proper Orthogonal Decomposition

La Proper Orthogonal Decomposition (POD) est une technique de construction d'une base réduite permettant d'approximer la solution d'un système dynamique [8]. Le principe consiste à collecter des snapshots de la solution à différents instants $t_1, ..., t_{n_s}$ où $n_s$ désigne le nombre de snapshots $\mathbf{y}(t_k)$. La matrice $\mathbf{Y}$ des snapshots est alors décomposée avec une Singular Value Decomposition (SVD) :

$$ \mathbf{Y} = \mathbf{V}\mathbf{\Sigma}\mathbf{W}^T \quad (8) $$

où $\mathbf{V}$ et $\mathbf{W}$ sont des matrices carrées unitaires et $\mathbf{\Sigma}$ est une matrice rectangulaire diagonale. Les coefficients diagonaux de $\mathbf{\Sigma}$, nommées valeurs singulières, sont positifs et rangés par ordre décroissant. En particulier, les colonnes de $\mathbf{V}$ sont les vecteurs propres de la matrice d'autocorrélation des snapshots $\mathbf{Y}\mathbf{Y}^T$. On réalise ensuite une troncature de cette décomposition en annulant les valeurs singulières à partir d'un certain rang $r$, choisi tel que la somme des carrés des valeurs singulières restantes représente au moins 99% de la somme totale. La matrice $\mathbf{V}_r$ de taille $N \times r$ obtenue en ne retenant que les $r$ premières colonnes de $\mathbf{V}$ est alors choisie comme base réduite. Cette troncature s'interprète comme une sélection des modes du système les plus énergétiques permettant de reconstruire les snapshots avec une très bonne précision. En fait, la base réduite $\mathbf{V}_r$ obtenue par ce procédé est celle qui minimise l'erreur de reconstruction des snapshots parmi toutes les bases réduites de taille $r$. En pratique, on observe pour de nombreux systèmes dynamiques que ce rang $r$ est généralement petit (e.g. 10-50) en comparaison de la taille initiale $N$ du système (e.g. $10^3 - 10^5$).

On souhaite à présent rechercher une approximation de $\mathbf{y}$ dans le sous-espace vectoriel engendré par l'image de $\mathbf{V}_r$, c'est-à-dire :

$$ \mathbf{y} \approx \mathbf{V}_r \tilde{\mathbf{y}} \quad (9) $$

avec $\tilde{\mathbf{y}}$ de taille $r$. On utilise pour cela une projection de Galerkin en projetant le système (7) dans l'espace engendré par la base POD :

$$ \frac{d\tilde{\mathbf{y}}}{dt} = \mathbf{V}_r^T A \mathbf{V}_r \tilde{\mathbf{y}} + \mathbf{V}_r^T g_{NL}(\mathbf{V}_r \tilde{\mathbf{y}}) \quad (10) $$

où l'on a pris soin de décomposer au préalable $g(\mathbf{y})$ en une partie linéaire $Ay$ et une partie non-linéaire $g_{NL}(\mathbf{y})$. Le produit $\mathbf{V}_r^T A \mathbf{V}_r$ de taille réduite $r \times r$ peut être calculé avant l'intégration en temps, si bien que la dimension initiale $N$ du système n'intervient plus, ce qui réduit considérablement le temps de calcul. En revanche, la dimension $N$ intervient toujours dans le calcul du terme non-linéaire projeté dans la base réduite (produits matrice-vecteur). Une autre technique de réduction doit être utilisée en complément de la POD.
## 3.2. Discrete Empirical Interpolation Method

L'objectif de la Discrete Empirical Interpolation Method (DEIM) est d'approximer une fonction vectorielle non-linéaire $f(t) = g(y(t))$ de taille $N$ à l'aide d'une base réduite [9]. Pour cela, on utilise des snapshots de $y(t)$ afin de former la matrice $F$ dont les colonnes sont les images par $g$ des snapshots $y(t_k)$. Une POD est ensuite appliquée à cette matrice afin d'extraire une base réduite $U$ de taille $N \times m$. L'objectif de la DEIM est alors de trouver une approximation de $f(t)$ dans le sous-espace engendré par cette base réduite. Pour cela, l'algorithme DEIM utilise la matrice $F$ en entrée et identifie $m$ indices distincts deux-à-deux $\sigma(1), ..., \sigma(m)$ tels que :

$$f(t) \approx U(P^T U)^{-1} P^T f(t) \qquad (11)$$

où $P$ est une matrice de taille $N \times m$ dont la i-ème colonne est le vecteur $e_{\sigma(i)}$ dont toutes les composantes sont nulles sauf la composante $\sigma(i)$ qui est égale à 1. Le produit matrice-vecteur $P^T f(t)$ consiste ainsi en pratique à n'évaluer que les $m$ composantes $\sigma(1), ..., \sigma(m)$ de $f(t)$ et ne fait donc pas intervenir la dimension $N$. Ainsi, en appliquant une DEIM à $g_{NL}(V_r \tilde{y})$ dans (10) on obtient :

$$V_r^T g_{NL}(V_r \tilde{y}) \approx V_r^T U (P^T U)^{-1} \tilde{g}_{NL}(V_r \tilde{y}) \qquad (12)$$

où $\tilde{g}_{NL}(V_r \tilde{y}) = P^T g_{NL}(V_r \tilde{y})$ est une évaluation de $m$ composantes de $g_{NL}$. De plus, la matrice $V_r^T U (P^T U)^{-1}$ de taille $r \times m$ est indépendante du temps et peut donc être calculée avant l'intégration en temps du système réduit. Enfin, l'approximation (11) satisfait la contrainte d'être exacte à tout instant pour les indices identifiés, d'où l'appellation d'interpolation discrète.

## 3.3. Système couplé réduit

Le système réduit POD-DEIM s'écrit finalement sous la forme suivante :

$$\begin{cases} \frac{d\tilde{y}}{dt} = V_r^T A V_r \tilde{y} + V_r^T U (P^T U)^{-1} \tilde{g}_{NL}(V_r \tilde{y}) \\ y(0) = y_0 \end{cases} \qquad (13)$$

Des comparaisons de performance (précision, temps de calcul, domaine de validité) seront présentées entre le modèle complet et le modèle réduit. En particulier, une étude sera menée sur l'influence de la sélection des snapshots sur les deux bases réduites $V_r$ et $U$ à manipuler.

# 4. Enrichissement par les données

Le modèle réduit développé a pour objectif de diminuer les temps d'évaluation du système sans pour autant dégrader la précision par rapport au modèle complet initial. Néanmoins, ce modèle est par construction insuffisant pour l'application visée. En effet, les hypothèses physiques retenues peuvent devenir fausses en condition opérationnelles et les paramètres du système peuvent varier en fonction du temps, compromettant la précision de l'approximation réalisée avec les deux bases réduites utilisées. Suivant l'avancement des travaux, plusieurs stratégies fondées sur l'assimilation de données simulées et/ou mesurées seront présentées pour enrichir ce modèle. On considérera notamment : l'interpolation de bases réduites locales et la réduction de modèle paramétrique [8], l'utilisation de filtres UKF ou PF [10] pour l'identification en temps réel des paramètres du système, ou encore le développement à partir de données d'un modèle de déviation prédictions-mesures en utilisant la DMD [11] ou des modèles autorégressifs [12].
# 5. Détection des instabilités de rotation

Un simple critère fondé sur le franchissement d'un seuil en amplitude permettrait en principe de détecter une anomalie de rotation. Cependant, une fois franchi ce seuil en amplitude, la progression de l'instabilité pourrait s'avérer trop rapide pour envisager d'appliquer une loi de contrôle sur le système pour la corriger. De plus, il serait difficile d'établir un tel seuil, car le système étudié pourrait présenter des cycles limites stables dont l'excentricité statique associée serait en dessous de ce seuil, mais dont la trajectoire permettant de l'atteindre franchirait une ou plusieurs fois ce seuil (voir Figure 4-a pour une intuition de ce phénomène). Suivant l'avancement des travaux, nous présenterons donc des techniques de classification fondées sur des algorithmes de machine learning pour détecter des précurseurs d'instabilité plus complexes qu'un seuil en amplitude [13].

# 6. Conclusion

Nous avons développé dans ce travail un modèle réduit de turbocompresseur à paliers à gaz que nous avons proposé d'enrichir avec des techniques d'assimilation de données simulées et mesurées. Ce modèle hybride tient compte du couplage non-linéaire entre la dynamique de la couche fluide dans les paliers et de celle du rotor. Il s'agit d'une première étape vers le développement d'un jumeau numérique pouvant identifier ses propres paramètres à partir des mesures sur le système réel qu'il représente, et de détecter des précurseurs d'instabilité de rotation du rotor. La suite de ce travail consistera à développer une hiérarchie de jumeaux numériques représentant différentes combinaisons de composants du turbocompresseur, et agissant à différentes échelles de temps, afin de décomposer l'influence des différents phénomènes modifiant les paramètres du système.

# Références

[1] L. Gu, E. Guenat, J. Schiffmann. A Review of Grooved Dynamic Gas Bearings, Applied Mechanics Reviews, ASME, Vol. 72, pp. 010802, 2020.

[2] M. Arghir, S. Le Lez, J. Frene. Butées et paliers aérodynamiques, Techniques de l'ingénieur, 2009.

[3] W. Liu, P. Bättig, P. H. Wagner, J. Schiffmann. Mechanical Systems and Signal Processing, Vol. 146, pp. 106983, 2021.

[4] J. Schiffmann. Integrated Design, Optimization and Experimental Investigation of a Direct Driven Turbocompressor for Domestic Heat Pumps, EPFL, Thèse n°4126, 2008.

[5] P. Bonello, H. M. Pham. The efficient computation of the nonlinear dynamic response of a foil-air bearing rotor system, Journal of Sound and Vibration, Vol. 333, pp. 3459-3478, 2014.

[6] E. Guenat, J. Schiffmann. Real-gas effects on aerodynamic bearings, Tribology International, Vol. 120, pp. 358-368, 2018.

[7] L. F. Shampine, M. Reichelt. The Matlab ODE suite, SIAM J. Sci. Comput., Vol. 18, pp. 1-22, 1997.

[8] P. Benner, S. Gugercin, K. Willcox. A Survey of Projection-Based Model Reduction Methods for Parametric Dynamical Systems, SIAM Review, Vol. 57, pp. 483-531, 2015.

[9] S. Chaturantabut, D. Sorensen, Nonlinear model reduction via Discrete Empirical Interpolation, SIAM J. Sci. Comput., Vol. 32, pp. 2737-2764, 2010.

[10] E. N. Chatzi, A. W. Smyth. The unscented Kalman filter and particle filter methods for nonlinear structural system identification with non-collocated heterogeneous sensing, Structural Control and Health Monitoring, Vol. 16, pp. 99-123, 2009.
[11] P. J. Schmid. Dynamic mode decomposition of numerical and experimental data, J. Fluid Mech., Vol. 656, pp. 5-28, 2010.

[12] M.D. Spiridonakos, E. Chatzi. Metamodeling of dynamic nonlinear structural systems through polynomial chaos NARX models, Computers & Structures, Vol. 157, pp. 99-113, 2015.

[13] T. Hastie, R. Tibshirani, J. Friedman. The elements of statistical learning, data mining, inference, and prediction, Second edition, Springer Series in Statistics, 2016.