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.