# Méthode VOF pour la simulation numérique de l'écoulement de l'air comprimé et du transfert thermique associé dans un piston liquide

El Mehdi Gouda, Mustapha Benaouicha, Thibault Neu, Lingai Luo, Yilin Fan

► To cite this version:

El Mehdi Gouda, Mustapha Benaouicha, Thibault Neu, Lingai Luo, Yilin Fan. Méthode VOF pour la simulation numérique de l'écoulement de l'air comprimé et du transfert thermique associé dans un piston liquide. 24e Congrès Français de Mécanique (CFM2019), Aug 2019, Brest, France. hal-02478668

HAL Id: hal-02478668

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

Submitted on 28 Sep 2020

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 VOF pour la simulation numérique de l'écoulement de l'air comprimé et du transfert thermique associé dans un piston liquide

E. GOUDA<sup>a,b</sup>, M. BENAOUICHA<sup>a,c,*</sup>, T. NEU<sup>a</sup>, L. LUO<sup>b</sup>, Y. FAN<sup>b</sup>

a. Segula Technologies. Cellule de Recherche et Innovation en ingénierie Navale et Energie. 19 rue d'Arras, 92000 Nanterre, France.

b. Université de Nantes, CNRS, Laboratoire de thermique et énergie de Nantes, LTeN, UMR 6607, F-44000 Nantes, France

c. Université de Caen Normandie, Laboratoire Universitaire des Sciences Appliquées de Cherbourg, LUSAC, EA 4253, 50130 Cherbourg-En-Cotentin, France.

* auteur correspondant : mustapha.benaouicha@segula.fr

## Résumé :

Cette étude porte sur la modélisation 2D d'un processus de compression d'air dans un piston liquide. Afin d'étudier numériquement cet écoulement, le choix a porté sur l'utilisation de la méthode VOF (Volume of Fluid) qui est adaptée à ce type d'écoulement. Le régime de l'écoulement varie au cours du processus de compression passant de laminaire à transitoire puis turbulent, ainsi la méthode LES (Large Eddy Simulation) a été choisie pour la modélisation de la turbulence. Des améliorations ont été apportées afin de s'approcher des résultats expérimentaux obtenus précédemment et de mieux illustrer les différents phénomènes physiques présents lors de la compression.

## Abstract :

This study focuses on 2D modeling of a process of air compression in a liquid piston. In order to study this flow numerically, the choice was made to use the VOF (Volume of Fluid) method which is adapted to this type of flow. The flow regime switches from laminar to transient, so the LES (Large Eddy Simulation) method was chosen for turbulence modeling. We tried to get closer to the experimental results and to better illustrate the different physical phenomena present during the compression using different approaches

**Mots clefs :** Piston liquide, Écoulement diphasique, Méthode VOF, air comprimé

## 1 Introduction

Le stockage d'énergie est un enjeu majeur pour permettre l'utilisation étendue des énergies renouvelables intermittentes. Dans ce contexte d'exploitation à grande échelle, le stockage d'énergie sous forme
d'air comprimé s'avère une solution prometteuse. Les technologies exploitées actuellement, basées sur des compressions adiabatiques, souffrent d'un rendement assez faible. La voie de la compression d'air quasi-isotherme pourrait être une alternative permettant d'améliorer le rendement et la rentabilité des systèmes de stockage d'énergie par air comprimé. L'utilisation d'un piston liquide pour la compression de l'air s'avère la solution la plus pertinente pour atteindre une compression quasi-isotherme. L'importante surface de contact parois/air en raison du petit diamètre, le contact eau/air, eau/parois, la faible vitesse d'avancement de la colonne d'eau et le travail réduit permettent au piston liquide de favoriser les échanges de chaleur et ainsi de limiter la variation de la température au cours de la compression/détente et donc d'assurer un processus quasi-isotherme [4].

L'intérêt de la compression de gaz par piston liquide a été mis en évidence par Van de Ven et Li [1] dans une étude numérique comparative entre un piston classique et un piston liquide en effectuant notamment plusieurs hypothèses sur les échanges thermique, cette étude a démontré que le rendement du piston pouvait atteindre 83% bien supérieur au rendement du piston classique qui est 70%. Plusieurs hypothèses sur l'écoulement et le transfert de chaleur ont été faites. D'autres études numériques et expérimentales ont suivi après, et ont essayé de s'approcher des conditions expérimentales de l'écoulement.

Le présent travail de recherche consiste à modéliser numériquement le processus de compression d'air dans un piston liquide dans le cadre du projet REMORA®(Système de stockage d'énergie électrique sous forme d'air comprimé en mer) [2] développé par Segula Technologies. Une étude sur l'optimisation énergétique du système a été présentée par Maisonnave et al. [3].

L'écoulement dans la chambre de compression est un écoulement diphasique (avec une interface gaz-liquide), non homogène et compressible auquel sont associés des transferts thermiques dus à la compression de l'air. Il existe toujours une méconnaissance sur l'écoulement induit dans le piston liquide et les phénomènes associés. Une compréhension de la nature de l'écoulement et des différents paramètres associés permettra une meilleure optimisation du processus de compression. L'objectif de cette étude est de mieux comprendre la physique de l'écoulement dans la chambre de compression et les transferts thermiques associés et de pouvoir reproduire l'écoulement observé lors de les études expérimentales menées précédemment par Neu dans [4], [5], [6] et [7]. Afin d'étudier numériquement cet écoulement, le choix a porté sur l'utilisation de la méthode VOF (Volume of Fluid) qui est adaptée au suivi de l'interface dans le cas de ce type d'écoulement. Le régime de l'écoulement varie au cours du processus de compression passant de laminaire à transitoire puis turbulent, ainsi le modèle de turbulence LES (Large Eddy Simulation) a été utilisé pour prendre en compte ces trois régimes de l'écoulement. Plusieurs paramètres et configurations sont ainsi testés et comparés aux résultats expérimentaux disponibles.

# 2 Méthodologie

## 2.1 Géométrie

Il a été choisi de simuler le problème en utilisant une géométrie simplifiée en 2D représentée sur la Figure 1. La hauteur de la chambre de compression est celle du cas expérimental [6], $L_{ini} = 2m$, et la distance entre les parois est égale au diamètre de la chambre, $D = 30mm$.
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Géométrie 2D du domaine choisi et conditions aux limites de l’écoulement</figcaption></figure>

Il est ainsi supposé que l'écoulement et l'aspect qualitatif des échanges thermiques associés ne dépendent pas de la direction normale au plan de simulation et que le problème est équivalent à une compression effectuée entre deux plans infinis. En effet, la quantité de chaleur dans l'air échangée avec l'extérieur dépend de l'aire des surfaces des parois. Or, pour le même volume d'air compressé, la surface de parois disponibles pour les échanges thermiques dans la modélisation cartésienne est différente du cas expérimental. Si on nomme $C_{2D}$ et $C_{3D}$ les rapports de compacité des cas 2D et 3D, on a :

$$C_{2D} = \frac{V_{2D}}{S_{2D}} = \frac{DL}{2L} = \frac{D}{2} \qquad (1)$$

$$C_{3D} = \frac{V_{3D}}{S_{3D}} = \frac{\pi \left(\frac{D}{2}\right)^2 L}{2\pi \frac{D}{2} L} = \frac{D}{4} \qquad (2)$$

La paroi supérieure de la chambre de compression des deux géométries (2D et 3D) est négligée puisqu'elle correspond respectivement à 0.75% et 0.375% de la surface pariétale totale. La compacité de cette modélisation est donc deux fois moins grande que celle du cas expérimental. Par conséquent, une quantité de chaleur moins importante sera échangée avec l'extérieur [14].

Le temps de calcul induit par le nombre important de simulations préalables à l'obtention de résultats numériques corrects est une ressource limitée du projet. Ainsi, il a été choisi de mener cette étude sur une géométrie 2D, ce choix peut engendrer des différences sur l'écoulement et sur les transferts thermiques, ces différences sont dues à la réduction des surfaces de contact entre l'air et les parois. L'influence de la compacité du domaine sur la dynamique de l'écoulement est pour le moment difficile à estimer. Une réduction des échanges de chaleur pourrait entraîner une variation de température plus importante par rapport à une compression 3D.

## 2.2 Conditions aux limites et conditions initiales

Les conditions aux limites imposées sont les mêmes que celles du banc d'essai de l'étude expérimentale de Neu [6]. Une vitesse d'entrée du fluide (eau) est imposée en bas du piston. Le mouvement du fluide compresse alors l'air situé au dessus. La distance parcourue par le fluide est 1.5m. La figure 1 et le
tableau 1 récapitulent les différentes conditions aux limites et initiales de l'écoulement.

Où $U_{pist}$ est la vitesse de déplacement du piston liquide, D est le diamètre de la chambre et $T_{ini}$ et $T_{ref}$

<table><caption>Tableau 1 – Conditions aux limites et conditions initiales</caption><thead><tr><td>Conditions aux limites</td><td>Conditions initiales</td></tr></thead><tbody><tr><td>$T_{ref} = 293K$</td><td>$T_{ini} = 293K$</td></tr><tr><td>$p_{ref} = 100600Pa$</td><td>$p_{ini} = 101325Pa$</td></tr><tr><td>Vitesse d'entrée $U_{pist} = 0.5m/s$</td><td>$U_{y,ini} = 10^{-4}m/s$</td></tr><tr><td>Adhérence aux parois $U(x, y, t) = 0$</td><td>$U_{x,ini} = 10^{-5}m/s$</td></tr></tbody></table>

sont respectivement les températures initiales et les températures imposées aux parois, $p_{ini}$ et $p_{ref}$ sont respectivement les pressions initiales de l'air et la pression atmosphérique de référence.

## 2.3 Critères de convergence

Pour calculer précisément les échanges thermiques entre les parois et l'air et modéliser la couche limite visqueuse, les bords du domaine doivent être raffinés [12]. La grandeur adimensionnalisée $y^+$ est utilisée pour calculer la dimension des mailles adjacentes aux parois selon la direction normale :

$$y^+ = \frac{yu_{\tau}}{\nu} \qquad u_{\tau} = \sqrt{\frac{\tau_w}{\rho}} \qquad (3)$$

Où $\nu$ est la viscosité cinématique du fluide, $\rho$ sa densité volumique, $\tau_w$ est le taux de cisaillement aux parois et $u_{\tau}$ est la vitesse de friction.

Pour une modélisation utilisant le modèle de turbulence LES, il est recommandé pour $y^+$ d'être inférieure à 1 [11]. On obtient pour cette modélisation une largeur de maille adjacente $y = 5.8e^{-5}m$, la valeur de $y^+$ a été fixée pour les différents maillages, ainsi, on assure un traitement égal de la couche limite dans tous les calculs réalisés.

Le schéma temporel utilisé étant implicite, la stabilité des résultats ne dépend pas, en théorie, de la valeur de la condition CFL [16]. Dans notre étude réalisée en 2D, la condition CFL qui est le nombre de Courant $Co$ est égale à :

$$Co = \frac{U_x \Delta t}{\Delta x} + \frac{U_y \Delta t}{\Delta y} \qquad (4)$$

Où $U = (U_x, U_y)$ est la vitesse de l'air, $(\Delta x, \Delta y)$ sont les pas d'espaces et $\Delta t$ est le pas de temps. L'écoulement de l'air a été observé comme globalement unidirectionnel, même si des recirculations ont été observées dans l'écoulement [4]. Par conséquent la composante $U_x$ de la vitesse est supposée nulle dans le calcul de la condition CFL. Dans les essais expérimentaux, la composante $U_y$ de la vitesse a été estimée supérieure à celle du piston [4]. Il a été choisi de calculer le CFL sur les valeurs moyennes de $\Delta y$ et de $\Delta t$. Ainsi, les études de convergence de maillages ont permis d'estimer que pour un calcul stable et rapide, la condition CFL s'établit à $Co = 0.1$. En première approximation, la composante $U_y$ sera ici supposée égale à celle du piston liquide, $U_{int} = 0.5m.s^{-1}$ dans le calcul de la condition CFL.

$$\frac{\Delta t}{\Delta y} = 0,2$$
## 2.4 Propriétés physiques de l'écoulement

L'écoulement dans la chambre de compression est un écoulement diphasique gaz-liquide. Le gaz comprimé est l'air, et la colonne du piston est composé d'eau.

$$ \rho_{eau} = 998\text{kg.m}^{-3} \qquad (5) $$

$$ \mu_{eau} = 1 \times 10^{-3}\text{kg.m}^{-1}\text{s}^{-1} \qquad (6) $$

$$ c_{p_{eau}} = 4.185 \times 10^{3}\text{J.kg}^{-1}\text{K}^{-1} \qquad (7) $$

Ainsi, le nombre de Reynolds $Re$ de l'écoulement dans la colonne liquide du piston est :

$$ Re = \frac{\rho_{eau}U_{pist}D}{\mu_{eau}} \approx 15 \times 10^{3} \qquad (8) $$

L'écoulement de l'eau dans le piston est donc transitoire, d'où la nécessité d'utiliser un modèle turbulence adapté. La température de l'eau au cours de la compression est supposée constante.

La compression de l'air est quasi-isotherme, toutefois, la variation de la température modifie ses propriétés physiques. Les propriétés initiales de l'air sont :

$$ \rho_{air} = 1.222\text{kg.m}^{-3} \qquad (9) $$

$$ \mu_{air} = 1.8 \times 10^{-5}\text{kg.m}^{-1}\text{s}^{-1} \qquad (10) $$

$$ c_{p_{air}} = 1 \times 10^{3}\text{J.kg}^{-1}\text{K}^{-1} \qquad (11) $$

Le champs de vitesse de l'air lors de la compression varie localement, la potentielle existence de zones de recirculation rend sa caractérisation difficile. On se basera sur la vitesse moyenne de l'air et la vitesse maximale obtenue pour estimer le régime de l'écoulement. Lors de la compression, la température de l'air évolue, l'approximation de Boussinesq est utilisée pour modéliser la variation de la densité de l'air en fonction de la température suivant l'équation 12.

$$ \rho - \rho_{ref} = -\rho_{ref}\beta(T - T_{ref}) \qquad (12) $$

$$ \text{Où : } \beta = -\frac{1}{\rho} \frac{\partial \rho}{\partial T} \bigg|_{p} $$

Il est supposé qu'il n'y a aucun transfert de masse entre l'air et l'eau, les transferts thermiques entre l'eau et l'air sont négligés.

## 2.5 Méthode VOF

Pour modéliser l'écoulement multiphasique de la compression d'air par un piston liquide, il a été choisi de réaliser la reconstruction de l'interface eau-air par la méthode VOF (Volume of fluid) sous le code commercial Ansys Fluent. Le couplage pression/vitesse est assuré par un algorithme PISO plus adapté aux écoulement fortement compressibles. L'algorithme PRESTO! est choisi pour la discrétisation spatiale de la pression. L'équation permettant de suivre l'interface dans son mouvement (équation d'advection sur le taux de présence) a été traitée numériquement soit à l'aide d'un schéma temporel implicite, ainsi l'algorithme "Modified HRIC (High Resolution Interface Capturing)" spécifique aux problèmes impliquant plusieurs corps et des forces de volume importantes a été retenu pour la reconstruction de l'interface eau-air.
Le paramètre "Volume Fraction cutoff" définit le seuil à partir duquel les fractions volumiques seront considérées comme entières, égales à 0 ou 1. Dans ce cas, toutes les cellules du domaine dont la fraction volumique est inférieure à $1e^{-6}$ ou supérieure à $1 - 1e^{-6}$ seront arrondies à 0 ou 1 respectivement. Il est supposé que la tension de surface a un impact faible sur l'écoulement et la gravité est négligeable dans l'air par rapport aux forces de pression, ainsi un traitement pour les forces de gravité ou la tension de surface n'est pas nécessaire.

## 2.6 Maillage

Deux types de maillages sont utilisés afin de choisir le plus adapté à notre cas. Ainsi un maillage rectangulaire structuré et puis un triangulaire non-structuré sont représentés dans les figures 2a et 2b.

Afin de mener une étude comparative, les deux domaines sont maillés avec le même nombre de mailles.

Les mailles au voisinage des parois sont aussi raffinées du même ordre de grandeur que pour le maillage structuré et les cellules ont une forme proche d'un triangle équilatéral.

Pour le maillage rectangulaire, au centre du domaine, le rapport d'aspect des mailles est supérieur à 10 et atteint 20 sur les bords du domaine. Ainsi, il a été choisi de s'intéresser et raffiner plus le maillage triangulaire, en respectant une qualité des éléments géométriques (skewness) proche de 1.

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Présentation des deux types de maillages utilisés</figcaption></figure>

Une étude de convergence de maillage a été menée et est présentée dans le tableau 2, celle-ci a été couplée à une adaptation du pas de temps pour satisfaire les critères fixées dans le paragraphe 2.3.

<table><caption>Tableau 2 – Étude de convergence des différents maillages testés</caption><thead><tr><td>Type de maillage</td><td>Nombre de mailles</td><td>Δt en s</td><td>CFL</td><td>Température moyenne en fin de compression (K)</td><td>Temps de simulation (h)</td></tr></thead><tbody><tr><td>Rectangulaire</td><td>75000</td><td>5.0 * 10<sup>-4</sup></td><td>0.085</td><td>427</td><td>13.5</td></tr><tr><td>Triangulaire</td><td>125000</td><td>5.0 * 10<sup>-4</sup></td><td>0.125</td><td>408</td><td>22</td></tr><tr><td>Rectangulaire</td><td>127000</td><td>5.0 * 10<sup>-4</sup></td><td>0.129</td><td>423</td><td>19</td></tr><tr><td>Triangulaire</td><td>233000</td><td>2.0 * 10<sup>-4</sup></td><td>0.133</td><td>396</td><td>61</td></tr><tr><td>Triangulaire</td><td>375000</td><td>1.6 * 10<sup>-4</sup></td><td>0.131</td><td>392</td><td>83</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Comparaison de l'évolution de la température au cours de la compression pour les différents maillages</figcaption></figure>

Pour un même nombre de mailles, les résultats obtenus avec un maillage triangulaire se rapprochent plus des résultats expérimentaux [6] et [5].

Le choix d'un maillage de 233k éléments au lieu de 125k éléments a permis de se rapprocher des résultats expérimentaux et donc de réduire l'écart de 6%.

En fin de compression, la température finale de l'air pour le maillage comptant 375k éléments est de $T = 392K$, contre $T = 396K$ pour le maillage à 250k éléments. Par rapport aux résultats expérimentaux, cela représente une progression globale d'environ 7%.

L'écart relatif entre les courbes numériques des maillages à 250k et 375k de la figure 3 sur l'ensemble de la compression est inférieur à 1%, pour une augmentation du temps de calcul de 20h. Par rapport aux ressources disponibles, il a été estimé qu'il n'était pas nécessaire de prolonger l'étude de sensibilité à la discrétisation en temps et en espace. En effet, cela engendrerait des simulations d'une durée supérieure à 80h pour un gain en précision marginal.

Il a été conclu que l'étude s'arrêtera à un nombre de mailles de $375 \times 10^3$ éléments et un pas de temps de $\Delta t = 1.6 \times 10^{-4}$.
<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Courbes de température moyenne de l’air obtenues pour les maillages cartésien et triangulaire</figcaption></figure>

Le maillage triangulaire montre une convergence plus rapide (figure 3), ainsi il a été choisi de raffiner plus ce maillage.

La courbe de température de la figure 4 obtenue avec un maillage triangulaire présente une caractéristique qui n'est pas présente sur la courbe du maillage rectangulaire. En effet, à $t \simeq 2.3s$ apparaît une variation importante de la température sous la forme d'une convexité, elle pourrait correspondre à la position relative de transition du régime laminaire au régime transitoire, une étude détaillé sur ce phénomène est présenté dans le paragraphe 3.3.

## 2.7 Choix du modèle de turbulence

La difficulté de caractériser le régime d'écoulement nous a amené à tester plusieurs modèles, ainsi plusieurs simulations ont été réalisées en laminaire, avec le modèle RANS choisi qui est $k - \epsilon$ RNG et LES (Large Eddy Simulation).

L'utilisation du modèle LES parait judicieuse du fait qu'il a été observé expérimentalement plusieurs phases d'évolution de l'écoulement. En fait le modèle LES permet de simuler les différents régimes d'écoulement. La turbulence est naturellement 3D [12], toutefois, pour certaines applications, une simulation 2D permet de prendre en compte l'effet de la turbulence sur l'écoulement[13].

La figure 5 compare l'évolution de la température moyenne de l'air dans le piston liquide pour les trois modèles de turbulence. Pour un maillage identique de 125k éléments triangulaires.
<figure><img src="image_6.png" /><figcaption>FIGURE 5 – Évolution temporelle de la température moyenne de l'air au cours de la compression pour les 3 modèles turbulences</figcaption></figure>

La moyennisation issue des modèles de type RANS (Reynolds Averaged Navier-Stokes) entraîne un lissage trop important de l'écoulement pour faire apparaître les structures tourbillonnaires que le cas expérimental a mis en évidence en fin de compression [4]. Les modèles LES et laminaire présentent des résultats similaires beaucoup plus fidèles à ce qui a été observé expérimentalement. Il est d'ailleurs difficile de déterminer à partir de ces figures et des observations empiriques le modèle le plus adapté au problème puisque la différence entre les résultats issus des modèles laminaire et LES est plus faible que ce qui a pu être distingué dans les essais expérimentaux.

# 3 Résultats

## 3.1 Champs de vitesse et de température

Les premières simulations montrent l'évolution de la vitesse et de la température en fin de compression suivant le modèle de turbulence.

La capacité du modèle laminaire à reproduire numériquement les tourbillons aperçus dans le cas empirique au moment où l'écoulement semblait turbulent est d'ailleurs une question intéressante quant à la nature même de l'écoulement. Les conclusions de l'étude expérimentale conjecturent la présence de deux régimes d'écoulement, d'abord laminaire puis turbulent [6] et [7]. Dans l'objectif de se rapprocher des résultats expérimentaux, le modèle LES est finalement choisi dans cette étude pour simuler le problème.
<figure><img src="image_7.png" /><figcaption>FIGURE 6 – Contour de température en fin de compression pour les différents modèles de turbulence</figcaption></figure>

Les résultats de l'évolution de la température présentés dans la figure 6 sont en adéquation avec ceux de la figure 5, la température est uniforme quand le modèle *k* — *ε* est choisi. Les résultats de l'étude expérimentale concluent à l'existence de zones localement plus chaudes. Le modèle laminaire et LES présentent des résultats quasi-identiques, ainsi, la température locale maximale atteint 478K. Les résultats expérimentaux disponibles dans [6] ont utilisé une méthode inverse pour le calcul de la température. La même méthode de calcul sera utilisée numériquement pour pouvoir mener une étude comparative.
<figure><img src="image_8.png" /><figcaption>FIGURE 7 – Contour de la composante $U_y$ en fin de compression pour les différents modèles de turbulence</figcaption></figure>

## 3.2 Reproduction du cas expérimental

Afin de reproduire numériquement les phénomènes physiques constatés expérimentalement, et donc de se rapprocher des courbes de températures empiriques, plusieurs paramètres sont introduits dans ce paragraphe.

La pertinence de la comparaison entre les résultats expérimentaux [6] et [5] et ceux issus de notre étude numérique repose sur la capacité à reproduire des conditions d’écoulement identiques. En effet, l’importance de paramètres tels que la géométrie de la chambre de compression et la vitesse de déplacement du piston liquide sur l’évolution de la température moyenne a été prouvée expérimentalement [6].

### 3.2.1 Vitesse de compression du piston

Le profil de vitesse expérimental présente de nombreuses oscillations qui proviennent principalement de l'utilisation d'un réseau d'eau d'un laboratoire soumis à des perturbations non-contrôlées. Il fût donc difficile de réaliser les essais expérimentaux avec un profil de vitesse constant. L'amplitude importante de ces oscillations présentes dans la figure 8b et la décroissance de la vitesse en fin de compression en font cependant une estimation relativement approximative. De plus, la vitesse de compression est un paramètre de premier plan pour les transferts thermiques et donc pour l'évolution de la température de l'air [9].

L'approximation d'une vitesse constante a une influence non-négligeable sur l'évolution de la température moyenne. Premièrement, on observe, dans la figure 8b en début de compression des oscillations, probablement dues aux variations importantes de la vitesse au départ de la compression, qui n'étaient pas présentes sur la courbe de température moyenne obtenue précédemment (figure 5). Puis, jusqu'à $t = 1.5s$, la température obtenue avec le profil de vitesse expérimentale est plus importante que celle pour la vitesse de compression constante. Il apparaît dans le graphe que l'approximation de la vitesse de
<figure><img src="image_9.png" /><figcaption>(a) Comparaison des profils de vitesse expérimentale et numérique</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Courbes de température obtenues pour une vitesse de compression constante et expérimentale</figcaption></figure>

FIGURE 8 – Profil de vitesse et son influence sur l'évolution de la température

compression par une constante est imprécis, ainsi, pour reproduire fidèlement les essais expérimentaux, il est nécessaire d'utiliser le même profil de vitesse de compression.

### 3.2.2 Mesures numériques des grandeurs

La mesure directe de la température moyenne d'une phase n'est pas proposée par Fluent pour les modélisations multiphasiques. On cherche donc à l'obtenir par une méthode indirecte.

Dans l'étude expérimentale, la méthode inverse calcule la température d'un champ homogène équivalent au champ de température réel. Ainsi, si on appelle $T_{mel}(t)$ la température du mélange eau/air que Fluent est capable de mesurer, $T_{moy-air}(t)$ et $T_{moy-eau}(t)$ la température moyenne des deux phases distinctes air et eau, $V_{air}(t)$ et $V_{eau}(t)$ les volumes occupés par les deux fluides et $L(t)$ la position du piston liquide à l'instant $t$, on a la relation barycentrique suivante :

$$T_{mel}(t) = \frac{T_{moy-air}(t) * V_{air}(t) + T_{moy-eau}(t) * V_{eau}(t)}{V_{air}(t) + V_{eau}(t)} \quad (13)$$

La température de l'eau est supposée constante au cours du temps et égale à 293K. De plus, $V_{air}(t) + V_{eau}(t) = V_{total}$. En divisant les volumes par la section du piston et en réarrangeant l'équation 13, on obtient :

$$T_{moy-air}(t) = \frac{T_{mel}(t) * L_{ini} - 293 * L(t)}{L_{ini} - L(t)} \quad (14)$$

L'écart maximal mesuré entre la position numérique et expérimentale du piston est de $2.10^{-3}m$, ce qui est satisfaisant. Cette mesure de la position du piston est donc adoptée et la méthode barycentrique pour calculer la température moyenne de l'air est appliquée pour les résultats présentés ci-dessous. La courbe obtenue est présentée sur la figure 9.
<figure><img src="image_11.png" /><figcaption>FIGURE 9 – Comparaison des températures moyennes de l'air mesurées sur la diagonale et par la méthode barycentrique</figcaption></figure>

## 3.3 Changement de régime

Comme évoqué dans le paragraphe sur le choix de maillage (2.6), la courbe de température obtenue avec le maillage triangulaire de la figure 4 présente une caractéristique qui n'est pas présente sur la courbe du maillage rectangulaire, celle-ci pourrait correspondre à une phase de transition du régime laminaire au régime transitoire. Afin de vérifier ceci, une comparaison sera faite avec les données expérimentales disponibles. Par analogie avec ce qui a été réalisé dans l'étude expérimentale, le quotient $q$ défini comme le quotient entre l'énergie interne et le travail de compression est calculé :

$$q = \frac{U_{int}}{W_{comp}} = \frac{mC_v(T(t + \Delta t) - T(t))}{\int_{t}^{t+\Delta t} P(V)dV} \quad (15)$$

Le calcul de l'énergie interne est effectué en différenciant la quantité d'énergie emmagasinée sous forme de chaleur par l'air entre les instants $t$ et $t+\Delta t$, soit $U_{int} = mC_v(T(t+\Delta t) - T(t))$ et le travail en calculant l'énergie nécessaire à l'avancée du piston liquide entre $t$ et $t+\Delta t$, soit $W_{comp}(t) = \int_{t}^{t+\Delta t} P(V)dV$. L'évolution de $q$ au cours du temps est décrite par la figure 10.

Une attention particulière est portée sur les paramètres présents dans la corrélation empirique obtenue dans l'étude expérimentale [7] et déterminant la position du piston liquide lors de transition entre le régime laminaire et l'apparition des tourbillons dans l'écoulement :

$$L_{tr}^* = \left( -0.0344L_{ini} + 109U_{pist}D^2 + \frac{0.0227}{D} \right) \left( \frac{p_{ini}}{p_{ref}} \right)^{-0.645*\sqrt{U_{pist}}} \quad (16)$$

Où $L^*$ est la position adimensionnalisée et normalisée du piston liquide où la transition de régime apparaît, $L_{ini}$ est la hauteur de la chambre de compression.

On s'aperçoit que la proportion d'énergie convertie en énergie interne augmente considérablement entre $L = 1m$ et $L = 1.1m$ puis diminue jusqu'à $L = 1.2m$. Cette brusque variation de l'énergie thermique emmagasinée par l'air explique la convexité observée" dans la courbe de température relevée précédemment de la figure 4.
La formule 16 permet de calculer la position théorique relative de la transition, dans notre cas ($L_{ini} = 2m, U_{pist} = 0.5m.s^{-1}, D = 30mm, p_{ini} = 101325Pa$), la position de transition est atteinte pour $L_{théorique}^* = 0.734$ ce qui correspond à $L_{théorique} = 1.469m$.

L'observation de ce phénomène physique remarqué dans les résultats expérimentaux [6] confirme la pertinence de l'utilisation d'un maillage triangulaire de bonne qualité par rapport au maillage cartésien composé de cellules étirées.

La figure 10 montre que, numériquement, cette transition (pic du quotient $q$) survient autour de $L = 1.1m$, alors que d'après la formule (16), celle-ci devrait être à $L = 1.469m$. Cette différence pourrait être expliquée par l'évolution plus importante de la température due à la géométrie 2D, en fait, l'augmentation de la température modifie les propriétés de l'air et donc précipite cette transition. Il faut aussi prendre en compte la marge d'erreur de la corrélation (16) qui est de 10% [4]. On peut alors en déduire que le pic présent sur la figure 10 correspond bien à une transition du régime convectif, en fait celle-ci correspond à un passage d'une régime laminaire à un régime transitoire/turbulent. Après ce pic, des fluctuations apparaissent sur la courbe du quotient $q$ ce qui confirme Toutefois, cette transition ne coïncide pas avec la position théorique qui prédit le surgissement du changement de régime convectif.

<figure><img src="image_12.png" /><figcaption>FIGURE 10 – Évolution de l'énergie interne relative au travail mécanique de compression</figcaption></figure>

## 3.4 Comparaison avec les résultats expérimentaux

Pour une modélisation de compression utilisant le profil de vitesse expérimental pour le déplacement du piston liquide, la courbe de température moyenne de l'air obtenue s'est rapprochée de celle du cas expérimental, ce qui était attendu. En fin de compression, l'écart entre la courbe expérimentale et numérique est d'environ $17K$, soit un progrès relatif de 42%. Les différences les plus importantes apparaissent au début et à la fin de la compression. En fin de compression, les oscillations expérimentales ont des amplitudes relativement plus faibles que celles de la courbe numérique. Cela peut être dû aux incertitudes de mesure, dont l'influence est plus significative en fin de compression, mais également à la difficulté de reproduire précisément le profil expérimental de vitesse au moment où son impact sur les échanges thermiques est le plus important.

La transition de régime thermique est difficile à distinguer sur la courbe de température obtenue avec
la vitesse de compression expérimentale compte tenu des variations engendrées par le profil de vitesse oscillant. Elle semble se produire un peu avant celle du cas à vitesse constant, à $t = 2s$.

<figure><img src="image_13.png" /><figcaption>FIGURE 11 – Courbes de température moyenne de l'air obtenues pour la vitesse constante et expérimentale de compression</figcaption></figure>

La courbe de température moyenne (spatiale) de l'air, calculée avec la méthode barycentrique présente des variations similaires à celles de la courbe expérimentale. Ainsi, on peut conclure que cette méthode permet de bien évaluer l'évolution de la température de l'air durant la compression. A partir de la valeur critique de transition, de petites variations apparaissent sur les valeurs de température, ceci permet de confirmer l'hypothèse émis plus tôt sur la transition d'un régime laminaire à un régime turbulent. Les trois courbes présentent des écarts assez conséquents en fin de compression. En fait, la présence de forts gradients de vitesses et de températures et le choix d'une géométrie 2D simplifiée semblent réduire les échanges de chaleur avec les parois (plus froides), ainsi la température de l'air est plus élevée que celle observée expérimentalement.

L'évolution de la température de l'air lors du processus de compression est de $70K$, ce qui correspond bien à une compression quasi-isotherme, comparée à une compression d'un piston solide où l'évolution de la température peut atteindre $350K$ [1].

## 4 Conclusion

La méthode VOF pour le suivi de l'interface couplée avec le modèle de turbulence LES sont un bon choix pour simuler l'écoulement et les transferts thermiques de la compression d'air par piston liquide. Ainsi, l'étude a permis, numériquement, de mettre en évidence le phénomène de changement de régime convectif apparu dans l'étude expérimentale [7], celui-ci correspond au passage d'un régime laminaire à un régime transitoire/turbulent.

La comparaison avec les résultats expérimentaux en suivant la même procédure de traitement des résultats a permis de valider la cohérence entre l'étude numérique et l'étude expérimentale et ainsi reproduire numériquement l'évolution de la température lors du processus de compression quasi-isotherme de l'air. Toutefois, l'écart en fin de compression est assez important, ceci est probablement dû au choix de sim-
plification de géométrie (2D).

Les perspectives de cette étude consistent donc à améliorer la modélisation 2D, en choisissant une géométrie 2D axisymétrique, se servir des corrélations d'échange thermique établies dans la thèse expérimentale [4] et construire une modélisation 3D pour caractériser précisément l'écoulement et les échanges thermiques associés. Pour cela une étude est en cours de réalisation et utilise cette fois la méthode de Lattice-Boltzmann (LBM) qui permet de mieux coupler les transferts thermiques aux écoulements fortement compressibles.

Références

[1] James D. Van de Ven, Perry Y. Li, Liquid piston gas compression, Applied Energy, 2009, doi:10.1016/j.apenergy.2008.12.001

[2] Thibault, NEU, Device and method for converting and storing electrical energy in the form of compressed air, U.S. Patent Application No. 15/578,061.

[3] O. Maisonnave, L. Moreau, R. Aubrée, M-F. Benkhoris, T. Neu, D. Guyomarc'h. Optimal energy management of an underwater compressed air energy storage station using pumping systems, Energy conversion and management, 165, 1:771-782, 2018. doi:10.1016/j.enconman.2018.04.007

[4] T. Neu, Étude expérimentale et modélisation de la compression quasi isotherme d'air pour le stockage d'énergie en mer, Thèse, École nationale supérieure Mines-Télécom Atlantique, 2017.

[5] T. Neu, D. Guyomarc'h, A. Bignon, C. Solliec. Modélisation de la compression d'air par piston liquide - application au stockage d'énergie renouvelable intermittente, Colloque Interuniversitaire Franco-Québécois sur la Thermique des Systèmes, 2015.

[6] T. Neu, D. Guyomarc'h, A. Bignon, C. Solliec. Elaboration d'une loi d'échange convectif dans le cas d'une compression d'air par piston liquide, AFM, Association Française de Mécanique. 2015.

[7] T. Neu, D. Guyomarc'h, A. Bignon, C. Solliec. Analyse de la transition du régime d'échange convectif lors de la compression isotherme d'air par piston liquide : Application au stockage d'énergie, SFT, Société Française de la Thermique, 2015.

[8] C. Zhang, B. Yan, J. Wieberdink, P.Y. Li, J. D. Van de Ven, E. Loth, & T.W. Simon, Thermal analysis of a compressor for application to compressed air energy storage, Applied Thermal Engineering, 73(2), 1402-1411, 2014.

[9] M.Mutlu, M.Kiliç, Effects of piston speed, compression ratio, and cylinder geometry on system performance of a liquid piston, Thermal Science, 20(6). 1953-1961, 2016.

[10] J.Sommeria, Sur la turbulence bidimensionnelle : une approche par la magnétohydrodynamique. Thèse, Université scientifique et Médicale de Grenoble; Institut National Polytechnique de Grenoble, 1985.

[11] Davidson, Large Eddy Simulations : how to evaluate resolution. International Journal of Heat and Fluid Flow, 30(5), 1016-1025, 2009.

[12] Yang Zhiyin. Large-eddy simulation : Past, present and the future. Chinese Journal of Aero-nautics, 28(1), 11-24, 2015.

[13] M. De Marchis, B. Milici, and E. Napoli. "Numerical observations of turbulence structure modification in channel flow over 2D and 3D rough walls." International Journal of Heat and Fluid Flow 56 (2015): 108-123.
[14] R. K. Shah, D.P.Sekulic, Fundamentals of heat exchanger design, John Wiley & Sons, p. 8-10, 2003.

[15] J.C. Neu, The dynamics of stretched vortices. Journal of Fluid Mechanics, 143, 253-276. 1984.

[16] R.Courant, K.Friedrichs, & H.Lewy. On the Partial Difference Equations of Mathematical Physics. IBM Journal of Research and Development, 11(2), 215-234, 1967.