Modélisation des transferts par ébullition dans les réservoirs cryogéniques de lanceurs Ariane V-A5ME

Modélisation des transferts par ébullition dans les réservoirs cryognéniques de lanceurs Ariane V-A5ME


(source de l'image : Air Liquide)

Abstract :  

This industrial study main objective is to present a model capable of predicting the vapor creation by boiling of a cryogenic fluid on a wall. This model has to be valid for different gravitational acceleration values and easily set up in the numerical code NEPTUNE_CFD.

Le principal objectif de ce bureau d'étude industrielle est de fournir un modèle capable de prédire l'ébullition nucléée d'un fluide cryogénique sur une paroi. Ce modèle doit être valable pour différentes valeur d'accélération gravitationelle et facilement implémentable dans le code numérique NEPTUNE_CFD.

 

Etudiants en charge du projet : 

Marylou GARNIER - 3ème année à l'ENSEEIHT - option Energétique et Procédés

Zacharie VATIMBELLA - 3ème année à l'ENSEEIHT - option Mécanique des fluides numérique

Encadrants du projet : 

Catherine COLIN - Professeur INP/ENSEEIHT - Institut de Mécanique des Fluides de Toulouse

Fabrice MATHEY - CFD Team manager at Air Liquide Advanced Technologies

               

Préface

Contexte industriel

Notre projet est apparu suite à la demande de l'industriel Air Liquide pour comprendre le phénomène d'ébullition en microgravité dans les réservoirs supérieurs de la fusée Arianne V qui contiennent de l'hydrogène et de l'oxygène.

 

D'où vient ce problème d'ébullition ?
Les phases de lancement d'une fusée se composent d'une phase propulsée et d'une phase balistique. Lors de cette dernière phase, la gravité est réduite et les seuls flux de chaleur qui sont transmis aux réservoirs supérieurs proviennent du rayonnement solaire et de la dissipation thermique des moteurs. Ces flux se trouvent être assez importants pour pouvoir mettre en ébullition les liquides cryogéniques. Il faut ainsi être capable de contrôler la pression et la température dans le réservoirs dans ce régime. Il devient dans ce cas nécessaire d'établir des modèles de prédiction de l'ébullition dans ce régime de microgravité. Il s'agit ainsi du cadre de notre projet dans lequel nous nous intéresserons plus particulièrement à l'oxygène liquide.

Objectif

L'objectif de ce BEI est d'établir un modèle d'ébullition nucléée en microgravité afin de l'implanter dans le code NEPTUNE qui sera exploité dans une autre partie de ce projet (se référer au projet  Ebullition dans les reservoirs cryogéniques de lanceurs Ariane V-A5ME, RENAUDIERE DE VAUX  Sébastien et DAVY Guillaume). Le travail va d'abord consister en l'extraction de données expérimentales fournies par Air Liquide pour en tirer des diamètres et de fréquence de détachement de bulles pour différentes gravités et différents flux. Par ailleurs une étude bibliographique sur l'ébullition en microgravité va permettre d'aboutir à l'étude de plusieurs modèles afin de comparer les résultats obtenus aux données expérimentales. Le modèle se rapprochant le plus des données expérimentales sera gardé. 

Introduction

Avant de rentrer dans le contexte de notre problème, il est nécessaire au lecteur d'avoir un aperçu général sur les transferts de chaleur avec changement de phase. Nous pouvons allons donc ici effectuer un bref rappel sur les différents régimes d'ébullition qui existent avant de nous attarder à notre problème.

 

  1. Régimes d'ébullition

La figure suivante illustre l'existence de différents régimes d'ébullition. Cette courbe d'ébullition issue de l'expérience de Drew & Müller montre l'évolution du flux de chaleur sortant d'une plaque plane horizontale lorsque l'on fait varier le paramètre $Tp-T_{sat}$ avec $T_p$ la température de la paroi et $T_{sat}$ la température de saturation du liquide qui est fixe pour un problème donné. C'est ainsi la température de la paroi que l'on fait varier.  

 

Source : cours "Transfert  changements de phase" , Catherine COLIN

Cette courbe met en évidence l'existence de différents régimes d'ébullition que l'on va expliquer ici : 

Dans le cadre de ce projet, nous nous intéressons exclusivement au régime d'ébullition nucléée. En effet chaque régime possède des modélisations qui lui sont spécifiques avec différentes problématiques industrielles associées. 

  1. Régime d'ébullition nucléée​.

Dans ce régime on a apparition de bulles qui jouent un rôle important dans le transfert de chaleur car en plus du régime de convection naturelle, le mouvement des bulles fait apparaître d'autres flux de chaleur. Les principaux paramètres qui vont nous intéresser dans ce projet concernant le régime d'ébullition nucléée sont les suivants :

Ces paramètres doivent être déterminés soit à partir de formules analytiques soit de manière expérimentale. Il est dans tous les cas nécessaire de les déterminer car ils joueront un rôle direct dans la modélisation des flux de chaleur du régime d'ébullition nucléée. Une modélisation de ce régime consiste à diviser la contribution du flux de chaleur en trois parties : 

Il va ainsi falloir modéliser ces flux et regarder leur évolution pour les différentes gravités. Ces problématiques feront l'objet des parties qui suivent.

Modèle d'ébullition nuclée

Dans cette partie est présentée un modèle d'ébullition nucléée (Kurul et Podowski) et une corrélation permettant de transposer le flux total à la paroi selon différentes valeurs de gravité (Raj et Kim).

 

Introduction

Le modèle de Kurul et Podowski est utilisé comme modèle de fermeture dans le logiciel Neptune. Un rappel de la modélisation des écoulements diphasique faite par Neptune est rapidement faite ici.

Le code NEPTUNE_CFD est un code de calcul polyphasique basé sur une approche Euler/Euler : l'évolution des grandeurs locales est décrite par trois équations de conservation moyennées dans le temps, chaque phase étant décrite à l'aide de variables continues.

Ces trois équations sont:

  • Équations de conservation de la masse,
  • Équation de conservation de la quantité de mouvement,
  • Équation de conservation d'enthalpie locale.

Ce système d'équations n'est pas fermé et nécessite des lois de fermeture pour des termes tels que les transferts interfaciaux de masse, de quantité de mouvement, d'enthalpie, le tenseur de Reynolds, etc...
Nous nous intéresserons dans ce travail à l'équation de conservation de la masse et à la fermeture du transfert interfacial correspondant.

L'équation de conservation de masse moyennée pour la phase $k$ est la suivante:
\begin{equation}
\frac{\partial}{\partial t}(\alpha_k \rho_k) + \frac{\partial}{\partial x_j} (\alpha_k \rho_k U_{k,j})=\Gamma_k
\end{equation}

Avec:
- $U_{k,j}$ : les composantes de la vitesse moyenne de la phase $k$
- $\alpha_k$ : le taux de présence de la phase $k$
- $\rho_k$ : la masse volumique de la phase $k$
- $\Gamma_k$ : le transfert interfacial de masse à la phase $k$

Le transfert interfacial  $\Gamma_k$ peut être décomposé en deux termes:

\begin{equation}
\Gamma_k = \Sigma_{p \ne k}\left( \Gamma^c_{(p \rightarrow k)}+\Gamma^e_{(\omega \rightarrow k)} \right)
\end{equation}

Avec :
          - $\Gamma^c_{(p \rightarrow k)}$ : transfert interfacial de masse au coeur de l'écoulement, de la phase $p$ à la phase $k$
          - $\Gamma^e_{(\omega \rightarrow k)}$ : création de la masse de vapeur à la paroi

Dans le cas étudié, le liquide est sous-refroidi (et le sous-refroidissement est important, de $5K$) donc la vapeur n'est créée que sur la paroi. Une fois la bulle formée et détachée, elle se recondense dans le liquide. Donc le terme $\Gamma^c_{(p \rightarrow k)}$ ne concerne que le transfert de la phase gazeuse à la phase liquide. La résolution de ce transfert est réalisée à partir du terme de transfert d'enthalpie de la phase gazeuse à la phase liquide, ce qui n'est pas détaillé ici.

Le transfert correspondant à $\Gamma^e_{(\omega \rightarrow k)}$ est résolu par un modèle d'ébullition nucléée. Le modèle utilisé par Neptune_CFD est celui de Kurul et Podowski.

Modèle à trois flux

1. Modèle à trois flux

Le modèle de Kurul et Podowski est un modèle d'ébullition nucléée présentée en 1990 très largement cité dans la littérature. C'est le modèle utilisé par Neptune_CFD comme fermeture du transfert de chaleur en paroi.

Ce modèle décompose la paroi chauffée en deux zones distinctes qui sont traitées séparément comme le décrit la figure suivante :

Dans la zone grise de surface $A_{nb}$ aucune bulle ne se forme. Ainsi, le seul flux de chaleur dans cette zone est un flux de convection : toute la chaleur sert à chauffer le liquide sous-refroidi.

Dans les zones blanches de surface totale $A_{eb}$ les bulles se forment puis se détachent à une fréquence $f$. Lors de la formation de la bulle, la totalité de la chaleur sert à évaporer le liquide à $T_{sat}$. C'est un flux d'évaporation net. Lorsque la bulle se détache, le liquide froid remouille la paroi à l'endroit où était la bulle. Un flux de conduction instationnaire modélise l'énergie nécessaire pour réchauffer le fluide sous-refroidi jusqu'à la température $T_{sat}$ à partir de laquelle une nouvelle bulle se forme.

L'aire $A_{eb}$ correspond à l'aire de la paroi occupée par la surface adimensionelle projetée des bulles.
\begin{equation*}
A_{eb}=K \pi r_d^2 N_a
\end{equation*}

Où :
       - $K$ est un facteur de projection
       - $N_a$ est la densité de site de nucléation, en site par mètre carré.

La somme des aires adimensionelles $A_{eb}$ et $A_{nb}$ doit donc être égale à l'unité :
\begin{equation*}
A_{eb}+A_{nb}=1
\end{equation*}

 

Finalement le flux total de chaleur transmis au fluide par la paroi, $\Phi_{\omega}$, est une somme de trois flux.
\begin{equation}
\Phi_{\omega}=\Phi_e + \Phi_c + \Phi_{ci}
\end{equation}

Où :
       - $\Phi_e$ est le flux d'évaporation net
       - $\Phi_c$ le flux de convection
       - $\Phi_{ci}$ le flux de conduction instationnaire

 

      2. Zone d'aire $A_{nb}$.

Le flux de convection entré par défaut dans Neptune_CFD est un flux de convection forcée. Or, dans la cas étudié, l'ébullition se fait dans un fluide au repos donc le phénomène de convection observé est un phénomène de convection naturelle.

La loi de Newton donne l'expression du flux de chaleur surfacique en convection (l'expression est ici pondérée par $A_{nb}$ car ce flux n'a lieu que dans les zones où il n'y a pas de bulle) :

\begin{equation}
\Phi_c=h (T_{\omega}-T_{l}) A_{nb}
\end{equation}

Afin de modéliser le coefficient de transfert $h$, la corrélation de Dropkin-Sommerscales à été validée par le deuxième groupe travaillant sur ce projet (Renaudière de Vaux Sébastien et Davy Guillaume) pour l'oxygène liquide.

\begin{equation*}
Nu=0.069 \ Ra^{1/3} \ Pr^{0.074}
\end{equation*}

Où :

  • $Pr$ est le nombre de Prandtl du liquide.

\begin{equation}
Pr=\frac{\nu_l}{D_l}
\end{equation}

  • $Ra$ est le nombre de Rayleigh du liquide
    \begin{equation}
    Ra=\frac{g \beta (T_{\omega}-T_l) L_c^3}{\nu_l D_l}
    \end{equation}
     
  • $Nu$ est le nombre de Nusselt
    \begin{equation}
    h=Nu \frac{L_c}{k_l}
    \end{equation}
     

Avec :
          - $\mathbf{k_l }$ : conductivité thermique du liquide
          - $\mathbf{g }$ :  gravité
          - $\mathbf{\beta}$ : coefficient d'expension volumétrique
          - $\mathbf{\nu_l}$  : viscosité cinématique et dynamique du liquide
          - $\mathbf{T_{\omega}}$ : Température de la paroi
          - $\mathbf{T_l}$ : Température du liquide
          - $\mathbf{D_l }$ : diffusivité thermique du liquide

La grandeur caractéristique $L_c$ choisie ici est la taille de l'élément chauffant.

   

             2. Zone d'aire $A_{eb}$ :

Dans cette zone, deux phénomènes ont lieu tour à tour : la formation de la bulle qui correspond à un flux d'évaporation net et le remouillage de la zone par le liquide froid ce qui correspond à un flux de conduction instationnaire.

Évaporation :

 

Le flux d'évaporation net est modélisé dans Neptune_CFD par:
\begin{equation}
\Phi_e=f \ \rho_g \ h_{lg} \ V  N_a
\end{equation}

Avec :
          - $f$ : fréquence de décollement des bulles
          - $\mathbf{\rho_g}$ : masse volumique du gaz
          - $\mathbf{h_{lg}}$ : chaleur latente de vaporisation
          - $V$ : volume de la bulle
                 $V=\frac{4}{3}\pi r^3$
         - $\mathbf{N_a}$ : densité de site de nucléation [$m^{-2}$]
 

Conduction instationnaire :

Le flux de conduction instationnnaire est modélisé dans Neptune par :

\begin{equation}
\Phi_{ci}= \sqrt{f} \ \frac{2 k_l (T_{\omega}-T_l)}{\sqrt{\pi D_l}} A_{eb}
\end{equation}

 

 

Fermeture

Ainsi le modèle de Kurul et Podowski permet d'apporter une fermeture quant au transfert de masse interfacial à la paroi. C'est un modèle d'ébullition nucléée, il n'est donc pas applicable dans les autres régimes d'ébullition (ébullition en film par exemple).

Malheureusement, le modèle de Kurul et Podowski à besoin lui aussi de trois lois de fermetures pour les termes introduits :

  • $N_a$ la densité de site de nucléation par unité de surface
  • $f$ la fréquence à la quelle les bulles se détachent de la paroi
  • $r_d$ le rayon des bulles lorsqu'elles se détachent de la paroi

Différents modélisations de ces paramètres ont été faites afin de fermer le modèle de Kurul et Podowski et sont présentées ci-dessous :

1. Densité de sites de nucléation.

L'hypothèse faite ici est que la variation de la gravité n'a aucune influence sur la densité de sites de nucléation, qui est déterminée par l'état de surface de la paroi et du mouillage du liquide sur cette paroi. La corrélation déjà présente dans Neptune_CFD est donc gardée :
\begin{equation}
N_a=\left[ 210 (T_{\omega}-T_{sat}) \right]^{1.8}
\end{equation}

2. Fréquence de détachement des bulles.

La corrélation utilisée par Neptune prend en compte la gravité et le diamètre de détachement des bulles.
\begin{equation}
f=\sqrt{\frac{4}{3}g\frac{\rho_l - \rho_g}{2 \rho_l \ r_d}}
\end{equation}

3. Diamètre de détachement des bulles.

La corrélation d'Ünal utilisée dans Neptune ne prend pas en compte la gravité : elle n'est valable que pour une gravité terrestre de $g=9.81$. Une autre corrélation doit être trouvée pour modéliser $r_d$, valable pour d'autres valeur de $g$.

La corrélation de Fritz donne une valeur du rayon de détachement en fonction de la gravité :
\begin{equation}
r_d= \sqrt{\frac{3}{8}} \ \sqrt{\frac{\sigma}{(\rho_l-\rho_g)g}} \ \beta
\end{equation}

Il est à noter que cette corrélation diverge lorsque la valeur de $g$ tend vers $0$. Ce qui est rassurant puisqu'en totale apesanteur les bulles ne se décollent pas de la paroi et donc la notion de diamètre de détachement n'a plus de sens.

 

Ces modèles seront comparés aux résultats expérimentaux tirés des vidéos fournies.

Corrélation de Raj et Kim

Dans un article publié en 2011, Raj et Kim établissent que deux modes d'ébullition différents existent, selon la taille de l'élément chauffant et de la longueur capillaire $l_c$:
\begin{equation*}
l_c=\sqrt{\frac{\sigma}{g(\rho_l-\rho_g)}}
\end{equation*}
Soit $L_h$ la taille de l'élément chauffant. Le critère permettant de déterminer le  régime d'ébullition concerné est :
\begin{equation*}
\frac{L_h}{l_c}=2.1
\end{equation*}
Si le rapport $L_h/l_c$ est inférieur à cette valeur, l'ébullition est gouvernée par les forces capillaires. Sinon elle est  gouvernée par les forces de flotabilité.
Dans notre cas,
\begin{equation*}
2.48  \le \frac{L_h}{l_c} \le 7.84
\end{equation*}
Ainsi nous nous trouvons toujours dans le cas où l'ébullition est gouvernée par la flotabilité. Dans ce cas, Raj et Kim proposent une corrélation permettant de déduire le flux total à une certaine gravité connaissant le flux à une gravité différente:
\begin{equation*}
\frac{\Phi_{\omega}}{\Phi_{\omega}^{ref}}=\left( \frac{g}{g^{ref}} \right)^{m}
\end{equation*}

Avec:
\begin{equation*}
m=\frac{0.9 \ T^*}{1+2.6\ T^*}
\end{equation*}
\begin{equation*}
T^*=\frac{T_{\omega}-T_{ONB}}{T_{CHF}-T_{ONB}}
\end{equation*}
Où $T_{CHF}$  est la température de la paroi à laquelle correspond le flux critique et $T_{ONB}$  la température de la paroi à laquelle se forme la première bulle.

Données expérimentales

Nous avons précédemment introduit les modèles que nous allons utiliser afin de prévoir les fréquence et rayons de détachement des bulles ainsi que la densité de site de nucléations. Il faut à présent s'attarder sur l'extraction de données expérimentales ce qui fait l'objet de cette partie. Il va d'abord s'agir de montrer comment les données vont être extraites à partir de vidéos fournies par Air Liquide puis de les traiter afin d'apporter une conclusion quant à l'évolution des paramètres qui nous intéressent. 

Protocole expérimental

Le dispositif expérimental utilisé par Air Liquide afin de mener leur expérience d'ébullition en microgravité est présenté sur la figure suivante : 

Source de l'image : Air Liquide

Pour des raisons de confidentialité, les dimensions du dispositif ne seront pas précisées. 

Description du dispositif : 
Le réservoir contient de l'oxygène liquide sous-refroidi. Le contrôle de la gravité est fait via un champs magnétique. Les parois du réservoir sont adiabatiques, une sortie contrôlée par une valve permet de maintenir une pression constante dans l'enceinte. Un flux de chaleur est imposé via un élément chauffant dans la paroi inférieure. Pour chaque vidéo et chaque mesure, un flux est imposé puis la température de la paroi inférieur est mesurée une fois le régime établi.

Vidéos d'ébullition en microgravité

  1. Vidéos d'ébullition

     Pour des raisons de confidentialité, les valeurs numériques des gravités correspondantes aux vidéos ne sont pas précisées ici.  Voici donc deux extraits de vidéos à flux de chaleur imposés identiques mais pour deux gravités différentes.

(source des vidéos : Air Liquide)

  1. Extraction des données

    De ces vidéos il faut extraire des rayons de bulles et des fréquences de détachement. La méthode utilisée est décrite ci-dessous.

  • Rayon de détachement
    Pour extraire des rayons de détachement, on isole une image à un instant donné et le plus grand nombre de rayons de bulles possible est mesuré sur celle-ci. Cette opération est répétée sur un grand nombre d'images afin de pouvoir construire un histogramme donnant la densité de population de bulle. Il faudra donc construire des classes de taille de bulles, ce qui sera l'objet de la prochaine partie.
    Montrons un exemple de mesure d'un diamètre de bulle. La figure suivante présente un image extraite d'une des vidéos. Un gros plan de la zone de détachement (cadre bleu) est utilisé :


    (source de l'image : Air Liquide)
    La figure suivante montre le zoom sur ce gros plan : 

    (source de l'image : Air Liquide)

    Sur cette dernière figure est tracée en rouge la dimension de la paroi sur laquelle les bulles se forment dont la taille est connue, ce qui permet d'avoir une référence de taille en terme de pixels. En bleu sont tracés deux exemples de mesure de diamètre de bulle. En mesurant un diamètre horizontal et un diamètre vertical puis en en prenant la moyenne, un diamètre en pixels est obtenu, que l'on peut convertir en millimètre grâce à la référence de la taille de la paroi connue.
    Cette opération est effectuée sur le plus grand nombre de bulles possibles et sur le plus grand nombres d'images possible afin d'obtenir un histogramme réaliste. 
     

  • Fréquence de détachement
    Cette donnée est plus délicate à extraire car il faut se placer sur un site de nucléation précis et déterminer le temps passé entre la formation de deux bulles successives. Cette opération doit être faite pour un maximum de cycles et de sites de nucléation. Un exemple de mesure est présenté ci-dessous. La série d'images suivantes montre un cycle complet de formation et détachement d'une bulle sur un site de nucléation :

    Cette série d'images est focalisée sur un site de nucléation pour lequel on voit apparaître une bulle. On  compte ainsi 9 images pour voir un cycle de naissance de détachement de la bulle jusqu'à apparition d'une nouvelle. Connaissant le nombre d'images par seconde enregistrées par la caméra, on déduit aisément la fréquence de détachement. Cette opération doit être répétée le plus de fois possibles et pour un nombre de site de nucléation le plus grand possible ainsi que pour toutes les gravités. 

Traitement des données

1. Gravité constante, variation du flux

Un premier set de vidéo fournit par Air Liquide présente les régimes d'ébullition pour différents flux de chaleur imposés à la paroi, la gravité restant constante.

De ce vidéos ont étés extraits des rayons de bulle au détachement, ce qui a permi d'obtenir les histogrammes suivants :

    

2. Flux constant, variation de la gravité

Un deuxième set de vidéo fournit par Air Liquide présente les régimes d'ébullition pour différentes valeurs de gravité, le flux de chaleur imposé à la paroi restant constant.

De la même façon des histogrammes ont été obtenus pour chacune des valeurs de gravité.

 

 

3. Affinement des histogrammes

Les histogrammes obtenus sont parfois très étalés. Ceci est dû au fait que certaines bulles coalescent après s'être détachées de la paroi et le diamètre mesuré correspond donc à plusieurs diamètres de détachement.

Un raffinement des histogrammes est fait, afin de ne pas prendre en compte ces diamètres trop importants. Un exemple d'affinage est présenté ci-dessous, pour la gravité G2.

Finalement le rayon moyen correspondant à l'histogramme affiné est utilisé comme valeur de rayon de détachement pour être comparé aux modèles choisis.

 

Résultats et Analyse

Cette partie est consacrée à la présentation des résultats ainsi qu'à leur interprétation.

Flux critiques en fonction de la gravité

Les premiers résultats que l'on a obtenu au cours de ce projet concernent d'une part l'évolution du flux critique en fonction de la gravité ainsi que l'évolution du flux en fonction de la température de paroi pour les trois gravités. Pour rappel le flux critique correspond au flux maximal atteint lors de la phase d'ébullition nucléée avant de passer en ébullition de transition. 

A noter que les résultats sont de nouveau présentés sans dimension, ainsi les flux, gravité et température sont sans unité.

 

  1. Évolution du flux critique de paroi

    La figure suivante montre l'évolution du flux critique en fonction de la gravité. La première observation que l'on peut faire est que le flux critique augmente avec la gravité. Cela signifie que diminuant la gravité, le domaine d'ébullition nucléée diminue également et que pour une température donnée en ébullition nucléée, le système est plus proche de l'ébullition de transition. 

 

  1. Évolution des flux de paroi en fonction de la température de paroi

    Les trois figures qui suivent montrent l'évolution des flux de paroi en fonction de la température adimensionnalisée $\frac {T_{paroi}-T_{saturation}}{T_{critique}}$ pour des gravités respectivement croissantes $G1$, $G2$ et $G3$. A noter que le paramètre que l'on fait varier ici est $T_{paroi}$.
    Ces trois courbes on été tracé pour un domaine de température allant au delà de la phase d'ébullition nucléée. En d'autres termes les phases d'ébullition nucléée, de transition ainsi de d'ébullition en film sont montrées ici. La phase d'ébullition nucléée apparaît clairement et correspond à la phase d'évolution quasi linéaire du flux en fonction de la température, c'est à dire sur l'intervalle [0.96 1]. C'est cette phase qui nous intéresse tout au long du projet, les autres phases faisant l'objet de modélisations différentes. Les flux ont été adimensionnalisés par le flux critique.
     

 

 

 

Traitement des vidéos

1. Corrélation donnant $\theta$ en fonction de la température

L'angle de mouillage dépend de la paroi (rugosité), du liquide mouillant, de la température et de la gravité. Comme ni le liquide mouillant ni la paroi ne changent, $\theta$ dépend de la température de paroi et de la gravité. Ce paragraphe présente une modélisation de la variation de $\theta$ avec la température.

Un premier set de vidéo fournit par Air Liquide donne les régimes d'ébullition pour différentes valeur de flux imposé, donc de température de paroi. Les trois vidéos sont faites pour une même valeur de gravité. Les vidéos sont traitées afin d'obtenir un rayon moyen de détachement des bulles pour chaque température de paroi.

La relation de Fritz est inversée afin d'obtenir l'angle de mouillage $\theta$ connaissant le rayon de détachement des bulles :

\begin{equation*}
\theta=\frac{180}{\pi} \ \sqrt{\frac{8}{3}} \ \sqrt{\frac{g(\rho_l-\rho_g)}{\sigma}} \ r_d
\end{equation*}
NB : l'angle ainsi calculé est exprimé en degré.

La figure suivante montre les résultats obtenus (en noir) et l'interpolation linéaire de ces résultats (en magenta).

Une interpolation linéaire est faite donnant une corrélation entre $\theta$ et la température de la paroi :

\begin{equation*}
\theta = 0.27276 \ T_{\omega} -17.748
\end{equation*}

 

2. Fréquence de détachement

Un deuxième set de vidéo, également fourni par Air Liquide, donne les régimes d'ébullition pour différentes valeurs de gravité, avec un flux imposé identique pour chaque expérience.

La figure suivante montre la comparaison entre les fréquences mesurées expérimentalement et celles donné par la relation utilisée par Neptune :

\begin{equation*}
f=\sqrt{\frac{4}{3}g\frac{\rho_l - \rho_g}{2 \rho_l \ r_d}}
\end{equation*}

La fréquence calculée analytiquement est du même ordre de grandeur que celle mesurée expérimentalement. L'erreur y est cependant assez importante ($\varepsilon=91\%$) ce qui est dû à la difficulté de mesurer correctement la fréquence de détachement étant donné la qualité des vidéos traitées.

 

3. Rayons de détachement

Finalement les rayons de détachement des bulles mesurés expérimentalement sont comparés avec les rayons donnés par la corrélation de Fritz :

\begin{equation*}
r_d= \sqrt{\frac{3}{8}} \ \sqrt{\frac{\sigma}{(\rho_l-\rho_g)g}} \ \theta
\end{equation*}

L'angle de mouillage utilisé pour G1 est celui donné par la corrélation trouvée précédemment. Comme cette corrélation n'est valable que pour G1, l'angle de mouillage a été déterminé pour G3 (ce qui explique que les deux rayons soient les mêmes pour cette gravité) et utilisé pour G2.

 

Les résultats sont très satisfaisants, l'erreur calculée est $\varepsilon=7\%$. La corrélation de Fritz est donc validée pour le fluide cryogénique et va être utilisée dans le modèle de Kurul et Podowski.

Prédiction des flux

1. Modèle de Kurul et Podowski

Dans un premier temps le modèle de Kurul et Podowski a été utilisé pour calculé les différents flux en gardant les expressions de la densité de site de nucléation $N_a$ et de la fréquence de détachement $f$ telles qu'utilisées dans Neptune_CFD.

Les courbes suivantes illustrent les résultats de cette simulation pour la plus faible gravité G1 et la plus importante G3.

On peut remarquer ici que le modèle de Kurul et Podowski ne parvient pas à reproduire le flux total tel que mesuré expérimentalement. Ce modèle doit donc être modifié.

2. Modèle de Kurul et Podowski optimisé

La modification a porté sur les paramètres $N_a$ et $f$ puisque les expressions utilisées dans Neptune_CFD ne sont pas utilisables pour des températures cryogéniques. Ainsi un facteur multiplicatif a été ajouté aux expressions de la densité de site de nucléation et de la fréquence de détachement et ce sont ces facteurs multiplicatifs qui ont été soumis à une optimisation.

L'optimisation a été faite pour la gravité médium, G2. On cherche les valeurs des facteurs multiplicatifs telles que le flux de chaleur à la paroi calculé par le modèle de Kurul et Podowski soit le plus proche possible du flux de chaleur à la paroi expérimental.

La figure suivante illustre cette optimisation :

Cette optimisation donne les expressions de $N_a$ et $f$ suivantes :

\begin{equation*}
 N_a= 2.7 \left[ 210 (T_{\omega}-T_{sat}) \right]^{1.8}
\end{equation*}
\begin{equation*}
f=2.5 \ \sqrt{\frac{4}{3}g\frac{\rho_l - \rho_g}{2 \rho_l \ r_d}}
\end{equation*}

Ces expressions ont ensuite été utilisées pour calculer les flux de chaleur donnés par le modèle de Kurul et Podowski aux deux autres gravités, G1 et G3.

Les deux figures précédentes montrent qu'après optimisation le modèle d'ébullition nucléée de Kurul et Podowski donne des résultats cohérents pour l'ébullition d'un fluide cryogénique en micro-gravité.

3. Corrélation de Raj et Kim

La corrélation de Raj et Kim permet de déduire le flux total de chaleur à la paroi pour une valeur de gravité connaissant le flux total à la paroi pour une gravité de référence.

Le flux utilisé comme référence est le flux expérimental correspondant à G3. Une interpolation quadratique est réalisée pour modéliser ce flux comme une fonction de la température. La Figure suivante illustre cette interpolation. L'interpolation donne:
\begin{equation*}
\Phi_{\omega}=409.5 \ T_{\omega}^2-8.00171 \ 10^{5} \ T_{\omega} +3.9044 \ 10^7
\end{equation*}
 

 

Les figures suivantes comparent les flux donnés par la corrélation de Raj et Kim et les résultats expérimentaux.

Cette corrélation de Raj et Kim donne des résultats satisfaisant pour modéliser le flux total passant à la paroi en fonction de la température. Cependant cette corrélation ne donne pas la répartition des différents flux; elle n'est pas capable de prévoir la quantité de chaleur utilisée pour évaporer le liquide, par exemple.

Comparaison des modèles

Comparaison des modèles

La figure suivante compare les trois modèles testés (Kurul et Podowski non optimisé, optimisé et Raj et Kim) pour la gravité minimale G1.

Les erreurs moyennées sur les trois gravités donnent :

- pour le modèle de Kurul et Podowski : $\varepsilon=80\%$

- pour le modèle de Kurul et Podowski optimisé : $\varepsilon=21.8\%$

-pour la corrélation de Raj et Kim : $\varepsilon=21.7\%$

La corrélation de Raj et Kim donne les meilleurs résultats mais ne permet pas de différentier les différents flux, notamment le flux d'évaporation qui permet de déduire la quantité de vapeur créée sur la paroi.

Le modèle de Kurul et Podowki optimisé est bien meilleur que celui utilisé pour l'instant dans Neptune_CFD et donne des résultats intéressants. Cependant il est notable que lorsque la gravité diminue, un palier apparaît pour le flux de chaleur calculé avec ce modèle, lorsque les flux deviennent importants. Ceci est dû au fait qu'à partir d'une certaine valeur de flux de chaleur à la paroi, la densité de site de nucléation atteint un maximum (toute la paroi est recouverte de bulles qui se forment puis se détachent) et ne peut plus augmenter. Ainsi ce modèle ne pourra pas reproduire la situation réelle pour des gravités faibles et des flux de chaleurs imposés à la paroi importants.

 

Conclusion et perspectives

Conclusion sur les résultats 

Les résultats obtenus nous ont donc permis de montrer que seul le modèle de fermeture du rayon de détachement $r_d$ donnait de bons résultats comparés aux résultats expérimentaux. En revanche les deux autres modèles de fermeture $Na$ et $f$ ont du être modifiés et optimisés afin d'obtenir des flux de paroi qui approchaient au mieux le flux total expérimental.

On peut donc conclure que le modèle de Kurul & Podowsky associé au modèle brut sur $r_d$ et aux modèles de fermeture optimisés sur  $Na$ et $f$ donnent des résultats qui sont cohérents et qui permettent de donner un premier modèle de prédiction sur l'ébullition nucléée en microgravité. 

 

Conclusion sur le projet

Étant donné le temps imparti pour un projet si ambitieux, on peut être satisfait des résultats obtenus. Néanmoins il s'agit d'un sujet compliqué et idéalement il faudrait effectuer un travail bibliographique plus approfondi sur l'ébullition nucléée afin de mieux comprendre le phénomène et tirer des modélisations plus évoluées. D'autre part, une amélioration de l'optimisation des paramètres $Na$ et $f$ est à envisager. De plus, la mesure des rayons et des fréquences de détachement expérimentaux présente beaucoup d'incertitude et pourrait être à l'origine des écarts observés avec les modèles. Enfin, d'autres séries d'expériences, ou des vidéos de plus grande qualité permettraient d'affiner les résultats en permettant d'obtenir une loi de fermeture pour l'angle de mouillage $\theta$ en fonction de la température de paroi et de la gravité.  

Bibliographie

Olivier Kannengieser, "Etude de l'ébullition sur plaque plane en microgravité, application aux reservoirs cryogéniques des fusées Arianne V", Thèse de doctorat de l'université de Toulouse, sous la direction de Catherine Colin et Wladimir Bergez, Toulouse,IMFT, 2009, 228 p.

Michaël MONTOUT, "Contribution au développement d'une Approche Prédictive Locale de la crise d'ébullition", Thèse de doctorat de l'université de Toulouse sous la direction de Catherine COLIN et Pierre-Antoine Haynes, Toulouse, IMFT, 2009, 257 p.

Rishi RAJ, Jungho KIM, John McQUILLEN, "On the scaling of pool boiling heat flux with gravity and heater size", January 2012

Kurul, N., Podowski, M. "Multidimensional effects in forced convection subcooled boiling". In Proceedings of the 9th Heat Transfer Conference (1990), pp. 21–26.

Fritz. Berechnung des maximalvolumens von dampfblasen. Physik Zeitschr., 36 :379–384, 1935.

Catherine COLIN, "Transferts avec changements de phase", cours master DET, Toulouse