Cas "Ebullition"

Tous les résultats qui seront présentés avec des calculs de simulation sous NEPTUNE_CFD ont été pris pour des temps suffisamment long pour pouvoir considérer l'écoulement comme stationnaire. Mais dans la réalité certains calculs ont eu des soucis de convergence dans les dernières itérations et nous avons décider de tout de même afficher les dernières valeurs correctes afin de faire nos comparaisons. Dans ce cadre les temps de calculs ne sont pas affichés mais il faut considérer le système comme à l'équilibre.

Une partie sur les erreurs rencontrées lors des calculs explicitera plus précisément les choix qui ont été pris tout au long du BEI.

Flux pariétal minimum

Dans le cas précédent le flux imposé en paroi ne permettait pas de déclencher l'ébullition. En effet, pour savoir si on peut avoir apparition de bulles il faut connaître la température de paroi avec le flux imposé. Nous sommes dans les conditions atmosphériques donc sous 1 bar de pression. A cette pression là, la température d'ébullition est de 100°C soit 373.15K. Sachant que l'on entre dans le tube une eau déjà à 363.15K pour être proche de la crise d'ébullition il suffit de calculer le flux nécessaire pour chauffer le liquide d'un ΔT de 10.

D'après la théorie, nous nous trouvons dans un cas de convection forcée. Le flux imposé en paroi peut donc se retrouvé comme étant:

$\varphi =h \Delta T$

Avec h le coefficient de transfert thermique. Ce dernier peut être exprimé à l'aide du Nusselt comme vu précédemment:

$Nu =  \frac{hD}{\lambda}$

Où D est le diamètre du tube et λ la conductivité thermique du fluide (ici fixée à 0.669 W/m.K pour une eau à 90°C).

Nous avons vu dans la partie précédente que le Nusselt pouvait s'exprimer à partir du Reynolds et du Prantl grâce à la corrélation de Gnielinski car nous sommes toujours dans les même conditions de travail (Re=10 723 et Pr=2).

On en déduit alors à partir du Reynolds et du Prantl et de toutes caractéristiques physiques de notre fluide en entrée la valeur du flux minimum à imposée en paroi pour déclencher la crise d'ébullition.

Avec  cette méthode dans une routine Matlab (calcul_flux_min.m) on trouve que le flux minimum à imposer en paroi est de 37 792 W/m²

Tant que l'on reste sous ce flux l'ébullition ne peut pas se produire par contre dès que l'on dépasse ce flux il y a risque de formation des premières bulles de vapeur d'eau.

Paramètres de configuration

Les paramètres sont les mêmes que ceux fixés sous Edamox dans le cas "chaud" hormis la température de paroi qui sera différentes car nous imposons des flux de parois croissants à la paroi conduisant à de l'ébullition.

Influence du flux pariétal

Pour valider ce qui vient d'être expliqué, nous avons procédé à des tests pour différentes valeur de flux mais toujours dans les même conditions de Reynolds, de Prandtl et toutes autres caractéristiques physiques mis à part. L'étude paramétrique a été lancé avec 5 flux allant de 10000 W/m² à 50000 W/m² avec un pas de 10000 W/m². L'ébullition devant se déclencher aux alentours de 35000 cela nous donne deux valeurs en dessous de ce flux critique, deux valeurs au dessus et une valeur un peu limite.

A priori si la crise d'ébullition à lieu cela se situera plutôt en fin de tube. Effectivement, le fluide ayant plus de temps de contact avec la paroi chaude lorsqu'il arrive en bout de tube il sera plus propice à voir sa température dépasser la température de saturation. Cette hypothèse sera validée ultérieurement.

      

Ce graphe montre bien que pour les valeur de flux imposée en paroi supérieur à 30 000 W/m² il y a apparition de bulles de vapeur.

Les bulles sont considérés ici au sens de vapeur d'eau. Les bulles étant petites il est difficile de pouvoir les observer. La fraction volumique de gaz est calculé pour être le complémentaire de la fraction volumique de liquide. Soit les deux équations suivantes:

$\alpha _{V}=\frac{V_{V}}{V_{T}}$

$\alpha _{V}+\alpha _{E}$

Où l'indice V représente la vapeur, le T pour total et le E pour eau liquide.

 

Concernant le chargement des données ous Edamox, nous vous redirigeons vers le lien suivant qui contient le dossier EBULLITION dans lequel se trouve le fichier de définition "param":

http://hmf.enseeiht.fr/travaux/bei/beiep/content/g10/telechargements

Dans ce fichier, ce sont les valeurs des flux imposés à la paroi qui ont été modifiées ainsi que la vitesse en entrée, le pas de temps, les valeurs de k et ε.

Evolution de la fraction volumique

Pour justifier l'hypothèse émise précédemment sur le fait que la vapeur aurait plutôt tendance à ce former en sortie du tube on propose d'étudier deux profils radiaux, le premier à y=1m (donc en sortie) et le second à y=0,5m donc à mi-hauteur dans le tube.

   

On est ici dans le cas du flux pariétal de 50 000 W/m² on a donc la crise d'ébullition qui arrive très tôt car la paroi est à environ 105°C. L'hypothèse première qui suggère que la formation des bulles est plus accentuée proche de la sortie du tube est ici validée car on atteint difficilement les 2% de fraction volumique de gaz à mi-hauteur contre près de 23% en sortie de tube.

Par la suite on ne s'intéressera donc qu'aux phénomènes qui se déroulent proche de la sortie du tube.

Validation avec la loi de Chen

On va donc maintenant valider notre code NEPTUNE_CFD en comparant les résultats d'écart de température entre la paroi et le cœur proche de la sortie du tube avec différents flux. On devrait observer deux choses; d'une part la corrélation de Chen comporte deux composantes, la première correspondant à l'énergie du flux de convection forcée et la seconde la partie de flux qui concerne l'ébullition nucléée, d'autre part il faut remarquer que cette partie de flux d'ébullition nucléée n'apparaît que pour les flux imposée supérieur au flux minimum.

      

La figure précédente permet de valider toutes les hypothèses faite au préalable. En effet, la courbe en rouge qui représente les valeurs expérimentale de flux imposée en paroi en fonction de la différence de température entre la température de la paroi et la température au cœur du fluide. On remarque que cette même courbe coïncide très bien avec la courbe théorique de la corrélation de Chen, qui à l'instar de la courbe rouge détermine le flux en connaissant la différence de température entre la paroi et le cœur.

On pourra également remarquer que la courbe noire se détache à partir des flux d'environ 30 000 se qui est en accord avec ce qui était attendu.

On conclura que pour les calculer par NEPTUNE_CFD sont bien de deux natures diverses. D'une part un flux de convection forcé servant à réchauffer le fluide proche paroi, d'autre part un flux d'ébullition nucléée qui sert à fournir l'énergie nécessaire au fluide pour le faire passer en phase vapeur.

Résultats

      

Comme on peut le voir sur cette image il se peut qu'au cours du temps il y ai des fraction volumique de gaz qui prédomine dans notre canalisation. Ce phénomène va vite générer des problèmes de stabilité du code. Toutefois il faut constater que le code a su tourner sans erreur jusqu'à 1.44 s et d'ailleurs il est aller un peu plus loin avant de s'arrêter.

Problèmes de convergence

Dans cette partie nous allons évoquer les erreurs qui peuvent valoir la peine d'être signalées dans le but d'améliorer les simulations futures. Effectivement lorsque les calculs durent plus longtemps il commence à y avoir une fraction volumique importante. Cette apparition de vapeur à pour effet de générer un plantage dans le code de calcul.

Pour réussir à prolonger le temps de calcul nous avons ajouter une cinquantaine de cycle de vérification APH. Ce sont les cycles de pré-calcul qui vérifient que toutes les valeurs de pression et vitesse ont bien convergés. En augmentant le nombre de cycle cela permet de converger plus à chaque cycle (surtout ceux qui ne convergeait pas auparavant) et donc repousse la possibilité de faire des calculs sur des temps plus important.

Malgré les 60 cycles qui font que les calculs convergent plus longtemps lorsqu'il y a trop de vapeur les calculs s'arrêtent tout de même. Nous avons donc regarder les origines possible de ce problème. Nous nous sommes donc aperçu que pour conserver le débit massique lors de la formation de bulle, cette dernière phase qui possède une masse volumique 2000 fois plus faible que celle du liquide fait que la vitesse augmente violemment:

$\rho U S =constante$

Il existe alors un taux de cisaillement fort entre le liquide au cœur de l'écoulement qui est à 0.5 m/s environ et les bulles de gaz qui sont maintenant à environ 265 m/s comme le montre cet imprime écran fait sur le listing de l'étape de calcul qui plante.

Les critères permettant de dire que cette itération est un problème sont les suivants:

-le critère de convergence imposé dans Edamox est de 1.10-6. Or ici dans la dernière itération la valeur du cycle est de 1.10-1.

  • La vitesse du gaz est de 265 m/s
  • La somme des fractions volumique des deux phases ne vaut plus 1.
  • La fraction volumique maximale de la phase est supérieur à 1.
  • Le calcul s'arrête dès l'itération suivante

                                     

Une hypothèse à également été faite sur la méthode de mise en température de fluide. Car nous avons directement injecté un flux d'enthalpie au niveau de la paroi dans Edamox. Cette méthode à le défaut de créer une marche d'escalier se qui rend le modèle brusque à traiter pour le code.

Il est possible de contourner ce problème à l'aide de la routine fortran usclim.F qui permet d'insérer une rampe d'énergie. Par exemple pendant 2 s l'enthalpie injectée passe de 0 à la valeur voulue de manière linéaire. Cela rend le traitement plus simple pour NEPTUNE_CFD car la méthode est plus douce.