# Simulation d'écoulements de gaz denses en turbines

Pietro Marco Congedo, Christophe Corre

► To cite this version:

Pietro Marco Congedo, Christophe Corre. Simulation d'écoulements de gaz denses en turbines. CFM 2009 - 19e Congrès Français de Mécanique, Aug 2009, Marseille, France. hal-03378579

HAL Id: hal-03378579

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

Submitted on 14 Oct 2021

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
# Simulation d'écoulements de gaz denses en turbines

P.M. Congedo, Christophe Corre

Laboratoire des Écoulements Géophysiques et Industriels (LEGI)
1023 rue de la Piscine, Domaine Universitaire, 38400 SAINT MARTIN D'HERES

**Résumé :**

Des gaz denses particuliers, dits gaz BZT, présentent en régime transsonique des propriétés atypiques qui peuvent être exploitées pour obtenir des rendements élevés dans des turbines les utilisant comme fluides moteurs. La difficulté d'étudier expérimentalement ces écoulements motive le développement d'outils de simulation efficaces. Un code turbulent (RANS) non-structuré sera appliqué à la simulation et à l'analyse de l'écoulement d'un fluide BZT dans une turbine en considérant différents points de fonctionnement.

**Abstract :**

Specific dense gases, so-called BZT gases, display non-classical properties when flowing in the transonic regime. These properties can be exploited to achieve high levels of efficiency when using such gases as working fluids in Organic Rankine Cycle-based turbines. The complexity of setting up experimental studies for such dense gas flows motivates the development of efficient numerical tools. A turbulent (RANS) code is applied to the simulation and analysis of the flow of a BZT fluid in a turbine for a range of thermodynamic conditions.

**Mots-clefs :**

simulation d'écoulement compressible, gaz denses, régime turbulent

## 1 Introduction

Les gaz denses sont des vapeurs monophasées dont les propriétés dévient de façon significative de celles prédites par la loi des gaz parfaits lorsqu'ils opèrent à des températures et des pressions proches des valeurs critiques. Un intérêt particulier s'est développé depuis quelques années pour une classe spécifique de gaz denses, connus sous le nom de gaz de Bethe-Zel'dovich-Thompson ou, en abrégé, gaz BZT. Pour de tels gaz BZT, des comportements non-classiques sont théoriquement prévus pour des conditions thermodynamiques situées au-dessus de la courbe de coexistence liquide / vapeur, telles que la dérivée fondamentale de la dynamique des gaz, définie comme $\Gamma = 1 + \frac{\rho}{c} \left( \frac{\partial c}{\partial \rho} \right)_s$, avec $\rho$ la masse volumique du fluide, $c$ la vitesse du son et $s$ l'entropie, devient négative. Lorsque $\Gamma < 0$, pour une gamme spécifique de température et de pression dans la phase vapeur, les classiques chocs de compression de la théorie des gaz parfaits violent l'inégalité d'entropie et ne sont donc plus admissibles. Les fluides BZT peuvent être des fluides utilisés commercialement dans des machines frigorifiques notamment. Les phénomènes non-classiques observés avec de tels fluides lorsqu'ils s'écoulent à grande vitesse (régime transsonique) ont plusieurs conséquences pratiques parmi lesquelles figure la possibilité d'utiliser ces fluides pour réduire les pertes causées par la traînée d'onde et les interactions choc / couche limite dans les turbomachines et les tuyères, plus particulièrement dans le contexte de cycles de Rankine Organiques (ou Organic Rankine Cycles, ORC en abrégé). La difficulté de mise en place d'études expérimentales pour de tels gaz (Colonna et al (2007)) a motivé le développement d'outils numériques pour analyser les performances associées aux écoulements de fluides
BZT, évaluer l'intérêt de leur utilisation comme fluides moteurs et définir leurs conditions op-
timales de mise en oeuvre. Transformer un code de simulation d'écoulements compressibles
initialement développé pour calculer des écoulements transsoniques de gaz parfaits en un sol-
veur d'écoulements de gaz denses exige de faire évoluer les traitements numériques de façon
à pouvoir traiter une équation d'état gaz complexe, telle que l'équation d'état de Martin-Hou
par exemple, indispensable pour décrire de façon correcte le comportement du gaz en écoule-
ment. Le choix de l'équation d'état impacte à la fois la formule de flux numérique utilisée pour
discrétiser le système des équations d'Euler, de Navier-Stokes ou de Navier-Stokes moyennées
(RANS) selon le degré de complexité retenu pour la modélisation, mais aussi les conditions
aux limites loin des parois solides où sont utilisés en général des invariants de Riemann qui
dépendent directement de l'équation d'état retenue. Des travaux précédemment réalisés par les
auteurs ont permis de développer un solveur non-structuré de type volumes finis pour des écou-
lements de gaz denses Congedo et al (2007) ; cet outil numérique a été récemment étendu à la
simulation de configurations tridimensionnelles et appliqué par les auteurs pour étudier la per-
formance aérodynamique d'écoulements de gaz denses en régime transsonique et en l'absence
d'effets de la viscosité. Un autre axe d'évolution de l'outil numérique a porté sur la prise en
compte des effets visqueux en régime laminaire puis en régime turbulent. L'objet du présent
travail est précisément de décrire les développements effectués pour modéliser les effets de la
turbulence dans des écoulements compressibles de gaz denses puis d'analyser des simulations
d'écoulements turbulents de gaz denses pour différents points de fonctionnement afin d'évaluer
de façon plus réaliste les gains d'efficacité potentiellement réalisables grâce à l'utilisation de
tels gaz.

## 2 Description physique des écoulements

Les écoulements de gaz denses étudiés sont régis par les équations de Navier-Stokes, exprimées
pour des écoulements à l'équilibre et non-réactifs. Dans le présent travail, on considère ces
équations sous leur forme intégrale, écrite pour un volume de contrôle Ω de frontière ∂Ω :

$$ \frac{\partial}{\partial t} \int_{\Omega} \mathbf{w} \, d\Omega + \int_{\partial\Omega} (\mathbf{F}_c(\mathbf{w}) - \mathbf{F}_d(\mathbf{w}, \nabla \mathbf{w})) \cdot \mathbf{n} \, d\Omega = 0, \quad (1) $$

avec $\mathbf{w} = (\rho, \rho\mathbf{V}, \rho E)$ le vecteur des variables conservatives ($\rho$ la masse volumique, $\mathbf{V}$ le vecteur vitesse, $E$ l'énergie totale spécifique), $\mathbf{F}_c = (\rho\mathbf{V}, \rho(\mathbf{V} \otimes \mathbf{V}) + p\mathbf{I}, (\rho E + p)\mathbf{V})$ les flux convectifs ou flux Euler ($p$ désigne la pression) et $\mathbf{F}_d = (0, \tau + \tau_R, (\tau + \tau_R) \cdot \mathbf{V} - (\mathbf{q} + \mathbf{q}_t))$ les flux diffusifs ou visqueux ($\tau$ désigne le tenseur des contraintes visqueuses, $\mathbf{q}$ le vecteur flux de chaleur, $\tau_R$ et $\mathbf{q}_t$ sont respectivement le tenseur des contraintes de Reynolds et le vecteur flux de chaleur turbulent). Le système (1) est complété par une équation d'état thermique et une équation d'état calorifique de la forme :

$$ p = p(\rho(\mathbf{w}), T(\mathbf{w})) \quad (2) $$

$$ e = e(\rho(\mathbf{w}), T(\mathbf{w})) \quad (3) $$

où $e$ désigne l'énergie interne spécifique, $T$ est la température.

### 2.1 Modélisation thermodynamique

L'équation d'état considérée dans ce travail est celle dite de Martin-Hou Cinnella et Congedo (2005) qui utilise cinq termes du développement du viriel et satisfait dix contraintes thermodynamiques. Cette équation est l'un des meilleurs modèles disponibles pour prendre en compte des
effets de gaz denses tout en maintenant un niveau de complexité de mise en oeuvre numérique raisonnable Cinnella et Congedo (2005) ; elle fournit une description précise du comportement du gaz au voisinage de la courbe de saturation en utilisant un nombre réduit de données d'en-trée thermodynamiques. Il est cependant établi que cette équation d'état tend à sous-estimer l'étendue de la zone dite d'inversion dans laquelle la dérivée fondamentale de la dynamique des gaz devient négative ; ce défaut n'est pas réellement problématique puisque les résultats de simulation obtenus en utilisant cette équation d'état auront simplement tendance à sous-évaluer (légèrement) les bénéfices attendus de l'utilisation de fluides BZT.

Dans le cas considéré ici de la simulation d'écoulements de gaz denses avec prise en compte des effets de la viscosité, des modèles thermodynamiques sont également nécessaires pour lier la viscosité dynamique $\mu$ et la conductivité thermique $k$ à la température et à la pression du gaz. Ces propriétés physiques des fluides BZT sont calculées en utilisant les modèles proposés par Chung et al (2005) :

$$ \mu = 40.785 \frac{F_c M^{1/2} T^{1/2}}{V_c^{2/3} \Omega_v} \qquad (4) $$

$$ k = 3.75 \times 10^{-3} R \mu \Psi \qquad (5) $$

Dans ces formules, $M$ est la masse molaire du fluide (en $g/mol$), $V_c$ est le volume critique en $cm^3/mol$, $\Omega_v$ est l'intégrale de collision, $R$ est la constante du gaz, $F_c$ et $\Psi$ sont des coefficients propres au fluide considéré qui dépendent, notamment, du facteur acentrique du fluide et de son moment dipolaire. Les modèles (4) et (5) ont été utilisés dans de précédents travaux (Cramer et al (1996)), (Cramer et Park (1999)).

Les effets de la turbulence sont pris en compte à l'aide du modèle de Spalart-Allmaras, dont les détails d'implémentation sont fournis dans la section qui suit. Il est important de souligner que l'hypothèse sous-jacente qui justifie ce choix est de supposer que, en première approximation, la structure de la turbulence n'est pas affectée par les effets de gaz denses, pour des conditions thermodynamiques suffisamment éloignées du point critique liquide / vapeur en tout cas. Dans le cadre d'une étude menée parallèlement à ce travail, nous nous intéressons à l'étude de la validité de cette hypothèse en menant des simulations des grandes échelles pour une configuration simple d'écoulement de gaz dense en canal et en faisant évoluer sur la base du comportement obtenu le modèle statistique que nous souhaitons mettre en oeuvre pour la simulation d'écoulements en géométries complexes (turbines 2D puis 3D).

## 2.2 Modélisation de la turbulence

Le modèle de Spalart-Allmaras est un modèle à une équation de transport pour la viscosité turbulente Spalart et Allmaras (1994), qui présente l'intérêt de ne pas exiger une résolution de maillage supérieure à celle associée à la mise en oeuvre de modèles de turbulence algébriques ; la version utilisée dans la présente étude est une extension du modèle initial à la simulation d'écoulements compressibles Deck et al (2002). L'équation de transport pour la variable turbulente de travail, $\rho\tilde{\nu}$, est donnée par :

$$ \frac{\partial}{\partial t} \int_{Omega} (\rho \tilde{\nu}) d\Omega + \int_{\partial \Omega} \left( (\rho \tilde{\nu}) \mathbf{V} + \left( \frac{1}{\sigma} (\mu + \rho \tilde{\nu}) \nabla \tilde{\nu} \right) \right) \mathbf{n} \, d\Sigma = \int_{\Omega} T(\mathbf{w}, \nabla \mathbf{w}) d\Omega \quad (6) $$

où le terme-source de (6) s'exprime comme :

$$ T = \left( c_{b_1} (1 - f_{t_2}) \tilde{S} \rho \tilde{\nu} + \frac{c_{b_2}}{\sigma} \nabla (\rho \tilde{\nu}) \cdot \nabla \tilde{\nu} - \left( \rho c_{\omega_1} f_{\omega} - \frac{c_{b_1}}{\kappa^2} f_{t_2} \right) \left( \frac{\tilde{\nu}}{d} \right)^2 \right) + \rho \tilde{\nu} f_{t_1} \Delta U^2 \quad (7) $$
La viscosité turbulente $\mu_t = \rho\nu_t$ qui apparaît dans l'expression du tenseur des contraintes de Reynolds est liée à la variable de travail $\rho\tilde{\nu}$ par la relation :

$$\mu_t = \rho\tilde{\nu}f_{\nu_1}$$

où la fonction $f_{\nu_1}$ est donnée par : $f_{\nu_1} = \frac{\lambda^3}{\lambda^3 + c_{\nu_1}^3}$ avec $\lambda = \frac{\rho\tilde{\nu}}{\mu}$.

Les différentes fonctions qui apparaissent dans (7) sont elles-mêmes définies par :

$$\tilde{S} = f_{\nu_3}(\lambda)S + \frac{\tilde{\nu}}{\kappa^2 d^2}f_{\nu_2}(\lambda) \quad (8)$$

avec

$$f_{\nu_2}(\lambda) = \left(1 + \frac{\lambda}{c_{\nu_2}}\right)^{-3} , \quad f_{\nu_3}(\lambda) = \frac{(1 + \lambda f_{\nu_1})(1 - f_{\nu_2})}{\lambda}$$

$$f_{\omega}(g) = g \left( \frac{1 + c_{\omega_3}^6}{g^6 + c_{\omega_3}^6} \right)^{1/6} , \quad g = r + c_{\omega_2}(r^6 - r) , \quad r = \frac{\tilde{\nu}}{\tilde{S} \kappa^2 d^2} \quad (9)$$

$$f_{t_1} = c_{t_1} g_t \exp \left( -c_{t_2} \frac{\omega_t^2}{\Delta U^2} [d^2 + g_t^2 d_t^2] \right) , \quad f_{t_2} = c_{t_3} \exp(-c_{t_4} \lambda^2) \quad (10)$$

Dans les formules ci-dessus, $S$ désigne la magnitude de la vorticité, $d$ est la distance entre le point courant et la paroi solide la plus proche, $d_t$ est la distance entre le point courant et la localisation de l'abcisse de transition (trip) à la paroi - supposée spécifiée -, $\omega_t$ est la vorticité pariétale au niveau de cette abcisse de transition et $\Delta U$ est la différence entre la vitesse au point courant et à la vitesse au niveau de l'abcisse de transition. La quantité $g_t$ est définie par $g_t = \min(0.1, \Delta U/(\omega_t \Delta x_t))$ où $\Delta x_t$ est le pas d'espace à la paroi au niveau de l'abcisse de transition. Enfin, les différentes constantes du modèle sont fixées aux valeurs suivantes :

$$\sigma = 2/3, \quad c_{b_1} = 0.1355, \quad c_{b_2} = 0.622, \quad \kappa = 0.41$$

$$c_{\omega_1} = \frac{c_{b_1}}{\kappa^2} + \frac{1 + c_{b_2}}{\sigma}, \quad c_{\nu_2} = 5, \quad c_{\omega_2} = 0.3, \quad c_{\omega_3} = 2$$

$$c_{t_1} = 1, \quad c_{t_2} = 2, \quad c_{t_3} = 1.2, \quad c_{t_4} = 0.5$$

Le flux de chaleur turbulent est classiquement décrit en introduisant un nombre de Prandlt turbulent égal à 0.9.

Comme on peut l'observer à l'examen des équations (6) et (7) l'équation du modèle de turbulence est résolue de façon totalement découplée des équations de conservation (masse, quantité de mouvement, énergie) du champ moyen.

### 3 Discrétisation en temps et en espace

Le solveur non-structuré s'appuie sur une discrétisation de type volumes finis du système de lois de conservation pour le champ moyen (1), formulée dans des maillages bidimensionnels généraux de type triangulaires, quadrangulaires ou mixtes. Les flux Euler sont calculés au niveau de chaque face d'une cellule du maillage en utilisant un solveur de Riemann approché de type HLL Luo et al (2005); une précision d'ordre 2 en espace est assurée grâce à une technique de reconstruction polynomiale (approche dite MUSCL) des variables conservatives. La reconstruction linéaire utilisée fait appel au gradient des variables conservatives dans chaque cellule du maillage; ce gradient est estimé à l'aide d'une technique de moindres carrés. Dans la
mesure où le solveur est appliqué à des écoulements compressibles en régime transsonique qui sont susceptibles de contenir des discontinuités, la formule de reconstruction inclut une étape de limitation du gradient, qui s'appuie sur l'approche désormais standard proposée dans Barth et Jespersen (1989) et revisitée dans Venkatakrishnan (1995), qui assure une capture de choc dénuée d'oscillation. Une convergence rapide vers l'état stationnaire est obtenue en implicant le schéma, suivant une procédure inspirée de Luo et al (1998), dans laquelle une simple phase implicite d'ordre 1 basée sur le schéma de Rusanov est associée à la phase explicite d'ordre 2 basée sur le schéma HLL ; ce choix est suffisant pour assurer la stabilité de la méthode lorsque de grands CFL sont utilisés. Comme, de surcroît, la phase implicite peut être résolue à faible coût à l'aide d'une méthode de relaxation par point, le schéma ainsi construit permet d'obtenir efficacement des solutions stationnaires des écoulements étudiés. Le flux numérique au niveau des faces frontières est calculé en utilisant des conditions basées sur la théorie des caractéristiques et l'introduction de cellules fantômes ; il est important de souligner que ces conditions doivent être adaptées en fonction des équations d'état décrivant le fluide étudié.

L'équation du modèle de turbulence est résolue en mettant en oeuvre les mêmes ingrédients : dans la mesure où une simple équation scalaire est à résoudre, le schéma de discrétisation spatiale retenu est un schéma décentré en fonction de la vitesse du champ moyen, étendu à l'ordre 2 par technique de reconstruction MUSCL. Il est nécessaire, pour le succès de la simulation RANS à l'aide du modèle de SA, de garantir des valeurs non-négatives de la viscosité turbulente en tout point de maillage et à tout instant. Pour satisfaire cette contrainte, nous avons suivi l'approche proposée dans (Spalart et Allmaras (1994)) qui introduit en particulier une linéarisation adaptée des termes d'advection, de diffusion et des termes sources dans l'équation du modèle. Au niveau d'une paroi solide et en champ lointain la variable $\tilde{v}$ est mise à zéro.

## 4 Simulations effectuées

Le code mis au point sur la base des modèles et des méthodes décrits ci-dessus est actuellement appliqué à la simulation d'écoulements turbulents de gaz denses en 2 dimensions d'espace, pour des configurations externes (profils d'aile) et internes (canaux inter-aubes). Ces simulations sont menées en faisant systématiquement varier les conditions thermodynamiques de l'état générateur afin de couvrir une gamme complète de fonctionnement ; les résultats obtenus sont comparés à ceux associés à un classique écoulement de gaz parfait dans les mêmes conditions de Mach, de Reynolds et d'incidence. Une analyse complète des résultats obtenus sera présentée lors de la Conférence, qui incluera également une comparaison entre certains des résultats obtenus et ceux récemment publiés dans Cinnella et Congedo (2008).

## Références

* P. Colonna, A. Guardone, N.R. Nannan 2007 Siloxanes : A new class of candidate BZT fluids. *Phys Fluids* **19** 086102

* P.M. Congedo, C. Corre, P. Cinnella 2007 Airfoil Shape Optimization for Transonic Flows of BZT Fluids. *AIAA J.* **45** 1303-1316

* P. Cinnella, P.M. Congedo 2005 Aerodynamic performance of transonic dense gas flows past an airfoil. *AIAA J.* **43**(2) 370-378

* T.H. Chung, M. Ajlan, L.L. Lee, K.E. Starling 1988 Generalized multiparameter correction of nonpolar and polar fluid transport properties. *Ind Eng Chem Res* **27** 671-679.
M.S. Cramer, S.T. Whitlock, G.M. Tarkenton 1996 Transonic and boundary layer similarity laws in dense gases. *J Fluid Eng* **118** 481-485

M.S. Cramer, S. Park 1999 On the suppression of shock-induced separation in Bethe-Zel'dovich-Thompson fluids. *J Fluid Mech* **393** 1-21.

P.R. Spalart, S.R. Almarras 1994 A one equation turbulence model for aerodynamic flows *La Recherche Aerospatiale* **1** 5-21

S. Deck, P. Duveau, P. d'Espiney, P. Guillen 2002 Development and application of Spalart-Allmaras one equation turbulence model to three dimensional supersonic complex configurations. *Aerospace Science and Technology* **6** 171-183

H. Luo, J. Baum, R. Löhner 2005 Extension of Harten-Lax-van Leer Scheme for Flows at All Speeds. *AIAA J.* **43**(6) 1160-1166

T.J. Barth, D.C. Jespersen 1989 The Design and Application of Upwind Schemes on Unstructured Meshes *AIAA Paper 89-0366*

V. Venkatakrishnan 1995 Convergence to Steady-State Solutions of the Euler Equations on Unstructured Grids with Limiters *Journal of Computational Physics* **118** (1) 120-130

H. Luo, J. Baum, R. Löhner 1998 A Fast, Matrix-Free Implicit Method for Compressible Flows on Unstructured Grids *Journal of Computational Physics* **146**(2) 664-690

P. Cinnella, P.M. Congedo 2008 Optimal airfoil shapes for viscous transonic flows of BZT fluids. *Comp Fluids* **37**(3) 250-264