Réalisation du cas d'ébullition

Après avoir validé le code en convection naturelle, il est maintenant nécessaire de passer à la géométrie réelle et au cas d'ébullition. Dans un premier temps, nous réalisons de l'ébullition d'eau (grace au module eau/vapeur de NEPTUNE).

Géométrie et maillage

Encore une fois nous sommes obligés de passer par des étapes de validation et de prise en main. Nous adoptons la géométrie du banc OLGA, avec une sortie modifiée. Le maillage 2D est relativement fin. Au centre, un élément chauffant sert à chauffer le fluide et à provoquer l’ébullition.

 

Maillage

Mise en place du cas sous NEPTUNE_CFD

Après avoir réalisé le maillage, il faut mettre en place le cas d'ébullition sous l'interface Edamox. Avant de commencer, l'utilisateur doit se placer en mode expert. L'interface de NEPTUNE_CFD est composée de plusieurs boutons et doit être parcourue comme suit :

interface

Modules spéciaux

NEPTUNE_CFD possède par défaut des modules spéciaux qui permettent de traiter des cas particuliers prédéfinis. Dans notre cas il faut sélectionner le module eau/vapeur qui est conçu pour traiter les écoulements bouillants en convection forcée. Néanmoins, nous verrons que le modèle reste pertinent pour de la convection naturelle.

 

modules

Propriétés du fluide et de l'écoulement

Ensuite, les propriétés propres aux deux phases eau et vapeur doivent être définies. Par défaut, avec le choix du module eau/vapeur, les valeurs sont entrées. Si ce n'est pas le cas, il faut entrer les valeurs à la main.

fluid

Attention, même si l'écoulement bouillant comporte des bulles dont le volume peut varier à cause de variation de pression et température, il ne faut en aucun cas cocher compressible, car l'écoulement n'est pas compressible.

Il faut ensuite ajouter des modèles de turbulence et de transferts entre phases :

fluid

Contrôle des entrées et sorties

Cet onglet permet le choix du nombre de pas de temps, de la fréquence des sorties pour le post-traitement etc. La modification de ces données ne change pas le calcul.

in_out

Généralités

C'est la partie la plus importante pour mettre en place un cas d'ébullition. Nous utilisons les tables CATHARES standard, pour de l'eau à 1bar.

generalities

Dans le cas d'un écoulement bouillant il faut toujours travailler en pression absolue et ne pas oublier de reporter la pression aux conditions limites de sortie.

 generalities

Pour les cas d'ébullition, il faut modéliser le flux de chaleur pariétal $\Phi_W$. Le modèle 4 flux utilisé est un modèle de Kurul et Podowski étendu (modèle de Seiler) :

\begin{equation*} \Phi_W = \Phi_C + \Phi_Q + \Phi_E + \Phi_{C2} \end{equation*}

  • $\Phi_C$ le flux de chaleur convectif entre la paroi et la phase liquide

  • $\Phi_{C2}$  le flux de chaleur convectif entre la paroi et la phase vapeur

  • $\Phi_Q$ le flux de Quenching : flux de condction instationaire

  • $\Phi_E$ le flux d'évaporation

Schémas numériques

La partie schémas numériques est laissée par défaut.

Scalaires

Lorsqu'on choisit le module eau/vapeur, deux scalaires sont automatiquement ajoutés. Il servent à transporter l'enthalpie des deux phases. Pour l'eau le scalaire 1 est défini comme nul à ($T_0$, $P_0$) et pour la vapeur le scalaire 2 est défini à enthalpie de saturation nulle à $P_0$, c'est-à-dire $H_{sat}(P_0)=0$.

scalars

Lorsque l'on impose $H_{sat}(P_0)=0$, cette condition impose la température de la phase vapeur à la température de saturation.

Conditions limites

C'est ici que les conditions limites doivent être définies, comme pour un cas classique.

boundary

Notre objectif est de chauffer une cellule en bas de notre cavité. Le liquide étant sous refroidi, nous voulons imposer des parois à une certaine température en dessous de celle d'ébullition. Imposer un flux de chaleur et des parois adiabatiques ne posent pas de problème, ni même qu'imposer une température de paroi et des parois adiabatiques. Il nous a cependant été impossible d'imposer des conditions mixtes, à savoir une température en paroi et un flux de chaleur sur une autre paroi adjacente. Pour contourner le problème, il faut définir la paroi à laquelle on veut imposer une température (sous forme de Dirichlet) comme étant une entrée, et imposer une vitesse nulle.

Contrôle des variables de sortie

C'est ici que l'utilisateur va définir les variables qu'il souhaite avoir pour son post-traitement.

Lancement du calcul

Nous pouvons à présent lancer nos calculs.

Résultats de simulation

Sans refroidissement

Dans un premier temps, nous nous sommes intéressés au changement de gravité. En prenant la géométrie précédente, à flux fixé, nous avons fait varier la gravité de 0 à 1g. Toutes les autres parois (hormis la cellule chauffante) sont adiabatiques. Nous chauffons à flux constant de 450 000 W/m² sans refroidir. L'eau est initialisé à $5K$ en-dessous de la température de saturation, soit à $367.14K$ à 1bar.

Il s’agit ici d’une phase de test. En effet il faut refroidir la paroi autour de la cellule chauffante. On voit cependant bien l’impact de la gravité sur les forces de flottabilité favorisant ainsi le détachement de petites bulles. Celles-ci quand elles existent se recondensent très rapidement dans l’eau qui est sous-refroidie. Nous pouvons aussi observer le champ de température :

De même, les effets de convection naturelle sont très présents pour les gravités fortes, alors qu'ils sont inexistants à gravité zéro.

Avec refroidissement

En faible gravité, le paquet de bulles a tendance à s'étaler. Pour éviter cela et rester dans le cadre de l'expérience OLGA, nous imposons une température de $367.14K$, soit un sous refroidissement de $5K$ sur les parois de la partie basse autour de la cellule chauffante.

Globalement le taux de vide reste similaire par rapport au cas non refroidi, mais la position des bulles est légèrement différente. Refroidir les parois atténue le fait que les bulles soient emportées par les gros tourbillons de convection. Seul en gravité 0 les bulles ont tendance à s'étaler plutôt qu'à être emportées par l'écoulement.

Comparaison des modèles de transfert de chaleur en paroi

Deux modèles sont disponibles pour modéliser le flux en paroi :

  • modèle à 3 flux de Kurul & Podowski,
  • modèle de Kurul & Podowski étendu à 4 flux.

Modèle à 3 flux

Dans ce modèle le flux de chaleur est composé de trois contributions :

\begin{equation} \Phi_W = \Phi_C + \Phi_Q + \Phi_E  \end{equation}

où :

  • $\Phi_C$ est le flux convectif pour le liquide,
  • $\Phi_Q$ est le flux de quenching (flux conductif instationnaire),
  • $\Phi_E$ est le flux d'évaporation.

Dans ce modèle, toute la chaleur sert à chauffer le liquide. Il n'y a pas de prise en compte des bulles en paroi.

Modèle à 4 flux

Ce modèle considère un quatrième flux $\Phi_{C2}$, qui est le flux convectif pour la vapeur. Ce modèle prend donc en compte les bulles à la paroi.

\begin{equation} \Phi_W = \Phi_C + \Phi_Q + \Phi_E  +\Phi_{C2} \end{equation}

Courbe d'ébullition

Nous simulons pour différents flux et pour chaque modèle, afin d'obtenir une courbe d'ébullition. La température de paroi est prise au centre de la cellule chauffante, afin de s'épargner des effets de bord. La figure de droite est un agrandissement de la zone avant flux critique.

On peut voir que le flux critique n'est pas atteint avec le modèle à 3 flux. Si on impose un flux plus élevé, il n'y a pas de convergence. Le modèle à 4 flux semble mieux se comporter. Nous pouvons aussi regarder l'évolution temporelle de $T_W$ pour un flux donné :

En régime établi, le modèle à 4 flux ne varie plus et $T_W$ est constante. Au contraire, le modèle à 3 flux a des variations périodiques. Celles-ci pourraient être dues à des détachements de bulles périodiques. Nous n'avons cependant pas d'explication pour les trois pics qui apparaissent à chaque fois en un point. Globalement le modèle à 4 flux impose une température en paroi plus chaude et la vapeur est plus surchauffée qu'avec le modèle à 3 flux.

Ebullition avec de l'oxygène liquide

Cette étape consiste en l'implémentation dans les sources du code de tables ou de corrélation pour les grandeurs physiques de l'oxygène liquide. Les tables CATHARE ne seront désormais plus utilisées (désactivées dans l'interface Edamox).

Les tranferts inter-phases d'énergie, de masse et de quantité de mouvement ne sont pas modifiés. Il s'agit de coder des lois d'évolution des grandeurs physiques dans la routine usphyv.F. C'est le tableau PROHY(NCELET,NPROP) qui doit être modifié.

Description du tableau PROHY(NCELET,NPROP)

Le tableau PROPHY contient les données physiques des deux phases. Il est structuré comme cela :

Description du tableau PROPHY
Entr​ée Signification
PROPHY(IEL,ITEMPK(IPHAS)) Température de la phase $i$
PROPHY(IEL,ICP(IPHAS)) Capacité calorifique de la phase $i$
PROPHY(IEL,IHSAT(IPHAS)) Enthalpie de saturation de la phase $i$
PROPHY(IEL,ITSAT) Température de saturation
PROPHY(IEL,IHLAT) Chaleur latente
PROPHY(IEL,ITHSIG)) Tension de surface à saturation
N/A viscosité, diffusivité de la chaleur

Ajout des propriétés dans le code

Nous avons implémenté un modèle de Boussinesq pour la phase liquide. Les $Cp$ sont supposés constants, et nous avons fixé l'enthalpie et la température à saturation.

Partie du code source du fichier Fortran 77 usphyv.F pour l'Oxygène liquide où PROPHY est une table de propriétés physiques, IEL est le numéro de la cellule :


C        Enthalpie de saturation
    hlsat = 178438.8
    hvsat = 102577.419
    Tsat = hlsat / Cp1
    PROPHY(IEL,IHSAT(1)) = hlsat
    PROPHY(IEL,IHSAT(2)) = hvsat
    PROPHY(IEL,ITSAT) = Tsat
C        Chaleur latente
    PROPHY(IEL,IHLAT) = abs(hlsat - hvsat)
C        Tension de surface
    PROPHY(IEL,ITHSIG) = 11.0d-3

Premiers résutats

Nous avons pour l’instant choisi de ne pas implémenter de loi d’évolution des grandeurs, mais de les laisser constantes pour les premières simulations. Seule la densité est variable avec un modèle de Boussinesq. Les conditions initiales sont les mêmes qu’au cas précédent, à savoir du LOx dans la cavité et de la vapeur à la sortie. Nous conservons le modèle à 4 flux car le modèle Kurul & Podowski de base échoue. En revanche, nous avons diminué le diamètre de détachement des bulles, en le passant de 1mm à 0,4mm. Ce diamètre correspond aux observations effectuées par le groupe chargé de la modélisation. Nous chauffons la cellule sans refroidir les parois adjacentes. Nous observons donc les évolutions des taux de vide pour différentes gravités. Au delà d'un certain temps néanmoins, les calculs divergent. 

 

Dans le dernier cas, en gravité zéro, les bulles ne montent plus mais coalescent pour former une grande bulle accrochée à la paroi :

De manière générale, les bulles sont peu refroidies après détachement et ne se recondensent pas.

Analyse numérique

Ces simulations nous ont donné accès à un certain nombre de paramètres : les différents flux du modèle 4-flux de NEPTUNE_CFD, la température de paroi Tw, les températures, les densités, les fractions volumiques ...

En ce qui concerne les flux, le flux total se stabilise rapidement avec le temps vers le flux imposé dans NEPTUNE. Ce qui est étonnant en revanche, c'est le faible flux de convection liquide qui représente 10% du flux total, et le faible flux d'évaporation (inférieur à 5%). La grande majorité (+80%) du flux total est contenu dans le flux de Quenching.

Évolution de la température de paroi en fonction de la gravité
Gravité (m/s²) Température de paroi Tw (K)
9,81 122.6
4,905 123.5
0,981 124.5
0 186

Comme dit précédemment, en 0g, toutes les bulles coalescent la cellule chauffe uniquement la vapeur, et la température de paroi croit fortement.

Ci-dessous, les profils verticaux de taux de vide sont représentés à différentes gravités. Plus la gravité diminue, plus la fraction de gaz en paroi est forte. En gravité 0,5g, des paquets de bulles se détachent et remontent périodiquement.

interface