Etudes de cas

Nous allons étudier 3 cas afin de valider au fur et à mesure les modèles utilisés dans le code NEPTUNE_CFD.

Tout d'abord un cas dit "Mono-Froid", le terme 'mono' signifiant monophasique, car dans le cas monophasique il n'y a pas de vapeur ni au début ni au cours de l'expérience. Dans ce cas l'eau reste à sa température d'entrée durant toute l'expérience car il n'y a aucun apport d'énergie.

Ensuite un cas dit "Mono-Chaud", de même que pour le premier cas on restera en monophasique tout au long de l'expérience. Mais cette fois-ci on imposera une densité de flux pariétal de manière à avoir de l'eau qui voit sa température augmenter en parcourant le tube mais tout en restant sous la température d'ébullition (ici 373.15K car dans les conditions atmosphériques).

Puis un cas dit "Ebullition", dans ce dernier cas on va faire monter la température de paroi au delà de la température de saturation du liquide et on verra l'apparition de bulle de vapeur.

Dans les 3 cas nous entrerons dans le tube avec de l'eau à 90°C. Ceci afin de pouvoir mettre l'eau en ébullition sans avoir à insérer une énergie trop importante.

Géométrie

Afin de se placer dans le cas réel d’une cuve enveloppant la cuve du réacteur tout en gardant le moins de complexité possible, le choix de la géométrie est celui d’un tube cylindrique vertical où l’écoulement progresse de bas en haut.

     

Propriétés du fluide

Le fluide considéré est de l’eau liquide à la température de 90°C et à la pression de 1 bar.

Cas "Mono-Froid"

Le cas mono froid signifie que la paroi ne voit pas apparaître de flux de chaleur et donc le liquide a un comportement de liquide dans une conduite. Nous voulons au traver de cet exemple valider le code pour des écoulement turbulent en conduite.

Le maillage

Ici le maillage choisi est intitulé « partiel2.des » avec 17050 cellules (cf partie maillage). Le cas du maillage « Partiel.des » aurait été suffisant mais le maillage « Partiel2.des » est choisi car il va être testé pour pouvoir être ré-utilisé ensuite en cas de validation. Le pas d’espace vertical  ΔY=6.45e-4 m, a été choisi de manière à permettre de capturer les phénomènes physiques (rapport de 3 avec la maille ΔX) sans pour autant avoir un rapport de 1 (maille cubique) avec le pas d’espace ΔX. L’intérêt d’avoir lâché du pas d’espace vertical est de permettre de diminuer le nombre de cellules et donc le temps de calcul. Ce temps de calcul va également diminuer en fixant le nombre de courant à  (cf partie paramètres « input » dans EDAMOX) puisque en prenant un pas d’espace plus grand, le pas de temps va augmenter pour un nombre de courant de 1 à une vitesse de 1,7 m/s.

L’intérêt de prendre un pas d’espace petit, c’est-à-dire réduire la maille, est que cela permet d’observer les gradients (on est en cas isotrope). Cependant étant donné que l’écoulement principal est suivant Y, on n’a besoin d’être plus précis dans cette direction et dans la direction Z. C’est donc pourquoi le maillage partiel2 a été choisi avec ce pas d’espace ΔY et dans le but de gagner en temps de calcul.

Les conditions initiales et limites

Dans ce cas simple, les conditions limites sont constantes et peuvent être directement définies dans l’interface de réglages des paramètres d’EDAMOX, sans devoir utiliser la routine FORTRAN usclim.F. Les valeurs ont été choisies :

  • Phase continue liquide en entrée :

 

  • Fraction volumique : 1
  • Vitesse axiale : 1,7  m/s
  • Turbulence :
  • Energie cinétique turbulente k= 0,0124m2/s2
  • Taux de dissipation de la turbulence ε = 0,9613m2/s3

 

  • Phase vapeur dispersée : il n’y en a pas.

A l’état initial, la conduite est complètement remplie d’eau, les vitesses sont nulles et la pression est constante et égale à la pression atmosphérique. La vitesse de 1,7 m/s a été choisie sur la base d’un cas test d’EDF et cette vitesse a été gardée dans le cas MONO_FROID. De plus le fluide en sortie est à pression atmosphérique.

Paramètres de configuration

La définition des paramètres de calcul est réalisée à l’aide du logiciel Edamox. Celui-ci est constitué de neuf interfaces dans lesquelles il va falloir paramétrer les valeurs afin de lancer la simulation.

Special modules

Ce bouton permet de choisir entre un cas eau/vapeur, eau/non condensable et une option aucuns. C’est la dernière qui sera cochée puisque pour le moment il n’y a pas formation de vapeur.

Fluid and flow properties

Dans cette rubrique nous définissons le nombre de phase à 1, l’eau liquide, puis remplissons ses caractéristiques comme vues précédemment, à savoir, sa masse volumique, sa viscosité dynamique (en Pa.s),… Par défaut un diamètre des particules est fixé à 0,003 mètres. Cette valeur n’est pas prise en compte dans l’algorithme puisque l’on ne travaille pas avec des particules mais avec de l’eau liquide dans laquelle on ne définit pas la taille des molécules d’eau. Cette valeur aura du sens lorsque l’on voudra définir des bulles de gaz en entrée par exemple. Le coefficient d’élasticité est également propre aux particules donc il n’intervient pas et sa valeur sera celle par défaut.  

Le modèle de turbulence choisi est un modèle k-ε à deux équations de transports. On ne choisit pas le modèle à 0 équations de la longueur de mélange donc celle-ci peut être laissée à 0.

Enfin la dernière option concerne les lois de parois avec l’option « WALL BC ». En cochant l’option friction, on considère qu’il y a frottement à la paroi, que la vitesse tangentielle est non nulle et que les gradients de vitesse radiale sont nuls.

                                     

Input and Output Control

Allocation des matrices

On peut définir la taille maximum des tableaux alloués au calcul (« Controlled memory allocation »). Il ne faut pas que la taille de la matrice d’allocation soit trop petite car le calcul ne se lancera pas. Si cette valeur est trop importante, le calcul prendra beaucoup de place. Nous fixons 100 000 valeurs d’entiers et 100 000 valeurs de nombres réels.

Temps de simulation

La gestion des entrées et sorties permet de contrôler la fréquence d’enregistrement des données dans des tables exploitables sous Paraview ou Ensight (option chrono qui va créer un fichier exploitable au format Ensight Gold), dans des fichiers « listing » (option listing qui permet de remplir un fichier texte au fur et à mesure du calcul).

Critères de simulation

Concernant le pas de temps, celui-ci peut être soit constant, soit variable en espace et en temps, ou bien variable en temps (« time dep »). Ce-dernier mode sera choisi et va dépendre des nombres de Courant et  Fourrier fixés.

On impose un nombre de Courant au maximum égal à 1 :

Pour ΔY=6.45e-4 m,  ΔtCourant=3,79e-4 s.

 

Pour le nombre de Fourrier,  soit  ΔtFourrier=12,8 s.

C’est dans l’onglet « CFL, FOURRIER LIMITS » qu’on fixe ces valeurs.

Sondes

Les sondes permettent de situer des points précis où seront mesurées les variables désirées et choisies dans l’onglet « Output Control ».

Moyenne temporelle

Cette option permet de calculer les moyennes temporelles ainsi que les fluctuations des moyennes de toutes les grandeurs calculées (activation dans la rubrique « average »).

Reprise d'un calcul

Il est possible de reprendre un calcul dans NEPTUNE grâce au fichier suiava. Celui-ci, présent dans le dossier RESU, doit alors être renommé suiamo puis être copié dans le dossier DATA.

                         

Generalities

On renseigne la pesanteur qui est dans la direction Y (-9.81 m/s2).  L’échelle L maximale des tourbillons se calcule via la formule suivante :

$L = \kappa \times \frac{D}{2}$

Où Kappa est la constante de Von Karman égale à 0.41 dans le cas d’une conduite circulaire de diamètre D=7 mm.

Il est possible de faire appel à des fichiers utilisateurs mais ce ne sera pas le cas ici (sinon il faut copier les routines qu’il y a dans /SRC/USERS dans /SRC).

Enfin dans notre cas la pression de référence a été fixée à 0, ce qui permet de travailler en pression relative et d’améliorer la précision des calculs. Lorsque l’on fera appel aux tables thermodynamiques Cathare dans la suite nous fixerons une pression absolue égale à 1 bar.

                  

Boundary conditions

On définit les 5 zones : entrée, sortie, mur, symétrie des murs, symétrie axiale ainsi que les références associées dans SIMAIL (entrée :10 , sortie : 90 , wall : 50, sym_axi : 41 , sym_wall = 40).

      

Une autre rubrique importante est le choix de la condition de pression dans chaque zone. Pour l’entrée et le mur, nous fixons une conditions « extr.P (i,w) » qui permet de faire une extrapolation du gradient de pression. Pour la sortie on fixe  « ddP=0 », c’est-à-dire qu’un profil de pression est imposée avec recalage sur la valeur de référence « Pref BC » (ici nulle). Enfin pour les deux symétries, on choisit une condition de flux nul (« dP=0 »).

La vitesse verticale est fixée à 1,7 m/s et les coefficients k et epsilon sont rentrés partout sauf au mur.

Pour les déterminer, on utilise calcule dans l’ordre suivant :

  • Le nombre de Reynolds : $Re = \frac{\rho U D}{\mu}$
  • L’échelle de longueur des grosses structures telle que: $L = \kappa \frac{D}{2}$
  • Le facteur de frottement à la paroi : $ f=0,316\times Re^{-0,25}$
  • Le coefficient de frottement pariétal : $C_{f} = \frac{f}{4}$
  • La contrainte de cisaillement à la paroi:  $\tau_{w}=0,5 \times C_{f}\times \rho \times U^{2}$
  • La vitesse Ustar définie par rapport à la couche limite, $U_{star}=\sqrt{\frac{\tau_{w}}{\rho}}$
  • L’énergie cinétique turbulente k approchée (U identique dans chaque direction et environ égal à ), soit: $k=\frac{3\times U_{star}^{2}}{2}$
  • Le taux de dissipation turbulente ε tel que : $\varepsilon = \frac{k^{1,5}}{L}$

Cela nous donne ainsi des valeurs de k=0,0124 m2/s2 et $\varepsilon$ = 0,9613 m2/s3.

Scalars

Les scalaires sont des grandeurs transportées par une phase porteuse (enthalpie, fraction massique, nombre de gouttes,…). Dans le cas froid, on n’injecte aucune température à la paroi, donc on ne signale aucunes grandeurs scalaires.

Variable Output Control

                                

Nous choisissons d’enregistrer la pression, les vitesses, l’énergie cinétique trubulente, le taux de dissipation et la viscosité turbulente. Pour choisir d’obtenir les données là où se trouve la sonde définie dans « INPUT » on rajoute son numéro dans la ligne de la grandeur à relever. Les résultats apparaîtront dans un fichier texte « hist ».

Run

Le code demande de fixer le nombre de cœur. Sur nos machines nous avons accès à 4 cœurs. Lorsqu’on ajoute un deuxième cœur pour faire un calcul en parallèle les processeurs échangent des données entre eux afin de se partager le travail. Il existe un nombre de mailles minimales par cœur qu’il faut atteindre afin de pouvoir ajouter un autre cœur au calcul. Ce nombre est de 10000 cellules par cœur. S’il y a trop de cœurs et pas assez de mailles, les processeurs passeront plus de temps à communiquer entre eux (échanges de résultats) qu’à faire le calcul. Ici nous sommes environ à 17 000 cellules, donc on choisit 2 cœurs.

Avant de lancer le calcul, il faut sauver notre fichier de paramètres et le nommer  « param ». Puis on peut lancer le calcul avec l’onglet RUN.

Exploitation des résultats

L’exploitation des résultats se fait avec Paraview et Matlab.

Les vitesses

Sous Paraview, nous regardons le champ de vitesse à mi-hauteur au temps t=1,55 s (les calculs convergeant au bout de 0,5 secondes). En passant à l’option Cell Data ton Point Datas, on obtient la figure d’en dessous, beaucoup plus homogénéisée.

Le choix de se placer à mi-hauteur est justifié par le fait qu’à cette hauteur la couche limite est pleinement développée et l’écoulement s’est stabilisé (cf profil d’énergie cinétique turbulente).

                            

                            

Afin de visualiser qualitativement les champs de vitesse, nous avons tracé les profils de vitesses radiales à différentes hauteurs y (y=0,2 m, y=0,5 m, y=1 m).

La vitesse de l’écoulement de 1,7 m/s imposée en entrée se retrouve dans le vert clair en frontière avec le jaune. La couche limite serait donc de l’ordre du tiers de la conduite.

Ces profils sont adimensionnés. La vitesse est normalisée par rapport à la vitesse moyenne (vitesse en entrée) et l’axe vertical x par rapport au rayon R=0,0035 m. Les profils apparaissent en général dans l’autre sens, à savoir que la paroi (ici en x/R=1) est normalement en x/R=0 (le maillage a été construit en prenant l’axe de symétrie comme référence (x= 0m)).

Nous pouvons faire plusieurs remarques :

  • Tout d’abord concernant la précision de la première maille. En x/R=1, la valeur de la vitesse normalisée est inférieure à 0,7.  Le raffinement du maillage à la paroi permettrait certes d’augmenter cette précision mais ce n’est pas l’aspect hydrodynamique qui est étudié ici.
  • Concernant la vitesse normalisée $\frac{\U}{\U_{moy}}=1$ (début de couche limite) on retrouve bien qu’à cette valeur on se trouve à $\frac{x}{R}=0,75$, c’est-à-dire au quart de la conduite. Cette couche limite vaut donc $\delta=2,6$ mm.
  • On observe que les résultats varient peu entre y=0,5 m et y=1 m, ce qui permet de dire que le calcul pourrait être mené sur une hauteur deux fois plus petite.

Nous avons cherché à comparer ces valeurs de vitesse avec les lois analytiques trouvés dans la littérature. L’équation de  Reichardt est :

$v(x)=U_{star} \left \{ \frac{1}{\kappa }ln\left ( 1+ \kappa\frac{xU_{star}}{\nu } \right )+c\left [ 1-exp\left ( -\frac{\frac{xU_{star}}{\nu}}{\chi } \right )-\frac{\frac{xU_{star}}{\nu}}{\chi }exp\left ( \frac{-0,33xU_{star}}{\nu} \right ) \right ] \right \}$

Avec $\kappa=0,41$ la constante de Von Karman, $\chi=11$ et $c=7,4$.

 

Nous avons tracé cette loi en inversant l’axe des abscisses de manière à faire la comparaison avec la figure ci-dessous. La différence étant que l’axe possède des valeurs négatives.

On observe des similitudes et notamment pour la première maille.

La première maille du maillage partiel 2 avait un pas d’espace Δx= 2,16e-4 m, soit un $x_{plus}$  d’environ 21. Ici en se plaçant à un $x_{plus}$  de 21 soit (-300+21=-279), on trouve une valeur de vitesse égale à 0,71 m/s. On retrouve exactement cette valeur sur le graphique avec les profils à différentes hauteurs. Cela permet de valider la précision de ce maillage pour les mailles proches de la paroi et de dire que l’espacement choisit est suffisant pour étudier le cas de chauffage de la paroi.

 

L'énergie cinétique turbulente

Comme pour la vitesse on observe que les résultats ne varient pas entre y=0,5 m et y=1 m, ce qui permet de dire que le calcul pourrait être mené sur une maille deux fois plus petite.

Comme on peut le voir l’énergie cinétique turbulente est plus importante en proche paroi qu’au centre du tube, ce qui est parfaitement normal car les fluctuations de vitesse en régime turbulent y sont plus élevées du fait des contraintes de cisaillement en proche paroi.

 

Remarque:

  • Les profils des figures 47 et 50 ont été tracés à partir de la première maille et non à la paroi.
  • Les profils des figures 48 et 49 contiennent en abscisses des signes négatifs dont il ne faut pas tenir compte.

 

Conclusion cas Mono Froid :

Au vu des résultats précédents, ce cas peut être validé, à savoir que les paramètres rentrés dans l’interface EDAMOX sont corrects et que le maillage partiel2.des peut être réutilisé dans les cas suivants.

Cas "Mono Chaud"

Dans le cas mono chaud on va désormais avoir un flux pariétal en paroi afin de voir la température du fluide s'élever de quelques degrés sans pour autant déclencher la crise d'ébullition. On pourra alors regarder la couche limite thermique turbulente et comparer aux exemple de la littérature.

Choix du flux pariétal

La grande différence entre le cas mono-froid et le cas mono-chaud est le fait que la paroi est chauffée. Cela implique de calculer au préalable le flux d'énergie que l'on décide d'injecter en paroi. Ce flux doit être judicieusement calculé afin de ne pas déclencher l'ébullition. Nous avons décider de se placer à un flux de 30 000 W/m² car le flux critique est environ de 37 000 W/m². Le calcul de ce flux critique est détaillé dans la partie ébullition. Mais il a été validé car l'observation de la fraction volumique de gaz qui reste égale à 0 dans tout le tube durant toute l'expérience avec un tel flux le permet.

Modification de la vitesse d'entrée

Des premiers tests ont montré qu'il est important de choisir une vitesse d'entrée du fluide de manière judicieuse. En effet si la vitesse du fluide est trop importante ce dernier n'a pas le temps de s'échauffer. Cela peut se comprendre si l'on compare deux temps caractéristiques. D'une part le temps de passage du fluide dans le tube qui est simplement la longueur du tube divisée par la vitesse d'entrée:

$t_{passage} = \frac{h}{U}$

D'autre part le temps de chauffage du fluide au repos qui peut être calculé par un bilan thermique. Ce temps est un temps de diffusion thermique étant en racine du coefficient de diffusion thermique divisé par le rayon du tube:

$t_{chauffage} = \sqrt{\frac{R}{D_{th}}}$

Avec h la hauteur du tube, U la vitesse en entrée, R le rayon du tube et $D_{th}$ le coefficient du diffusion thermique.

On peut donc dire que pour que le fluide ait le temps de se réchauffer il faut que le temps de passage soit supérieur au temps de chauffage. On trouve alors que plus la vitesse est petite plus il est possible de chauffer le fluide. Ce qui paraît tout à fait logique.

D'autre part nous allons par la suite calculer des nombre de Nusselt. Or la correlation utilisée pour calculer le Nusselt (Gnielinski) impose d'avoir un nombre de Reynolds supérieur à 10 000.

$Re= \frac{\rho D U}{\mu}$

Ce qui donne une borne inférieure pour la vitesse de 0.47 m/s afin d'être au delà de 10 000.

De plus plus la vitesse est faible plus les calculs seront rapide car le pas de temps est calculé afin de garder un nombre de Courant égale à 1 (Selon ce qu'on lui impose et pour garder un code stable). Le pas de temps est donc calculer de telle manière:

$\frac{C_{o}dZ}{U}$

Avec Co le nombre de Courant souvent égale à 1 ou moins, $\Delta Z$ le pas d'espace selon la direction de l'écoulement et U la vitesse du fluide.

Plus le pas de temps est grand moins il y a d'itérations pour un même temps de simulation. Il est donc important de chercher à toujours augmenter ce pas de temps.

Le choix le plus raisonnable est donc de prendre une vitesse de 0.5 m/s en entrée pour avoir un nombre de Reynolds de 10 723 tout en s'assurant que le fluide sera réchauffé au cours de l'expérience.

Paramètres de configuration

Special modules

Comme nous imposons un flux pariétal il faut une enthalpie en scalaire. De plus, nos voulons vérifier la fraction volumique de fluide et de gaz dans le tube pour s'assurer qu'il n'y a pas d'apparition de vapeur au cours de l'expérience. Pour ces deux raisons nous décidons de nous placer dans le module eau/vapeur qui permet d'avoir accès directement à tous ces paramètres.

Scalars

Comme nous imposons un flux pariétal il faut une enthalpie en scalaire. De plus, nos voulons vérifier la fraction volumique de fluide et de gaz dans le tube pour s'assurer qu'il n'y a pas d'apparition de vapeur au cours de l'expérience. Pour ces deux raisons nous décidons de nous placer dans le module eau/vapeur qui permet d'avoir accès directement à tous ces paramètres.

Boundary conditions

Enfin il suffit de mettre un flux pariétal différent de 0 en paroi (ici 30 000 afin de ne pas déclencher l'ébullition) pour le fluide 1 et rien pour le fluide 2.

Profils de température

Une fois les paramètres choisit dans edamox nous avons tracé l'évolution de la température pour un temps très grand afin d'être sur d'avoir passé la phase transitoire (typiquement 1.5s suffit nous avons pris 2s pour plus de sûreté).

Voici les profils de température le long du tube:

Nous observons que la température est plus chaude près de la paroi chauffante ce qui reste logique et que les profils de température sont de plus en plus haut lorsqu'on s'approche de la sortie de du tube. On peut alors tirer deux conclusions de cela:

-La formation de bulle se fera pour la suite préférentiellement proche paroi

-La formation de bulle se fera pour la suite préférentiellement proche de la sortie de tube

Pour valider cela voici deux profils axiaux de température, un proche paroi et l'autre au coeur du tube.

  

Nous observons qu'il y a très vite dans le tube un ΔT (différence de température entre la paroi et le fluide au coeur) constant (au delà de 0.2 m de haut) ce qui impose d'avoir un coefficient d'échange thermique constant et donc un nombre de Nusselt constant également. Tant que l'écoulement est développé en temps et en espace.

Validation avec le Nusselt

Selon cet écart de température on peut remonter au coefficient d'échange thermique en connaissant le flux pariétal imposé en paroi.

$h=\frac{\varphi }{\Delta T}$

Avec h le coefficient d'échange thermique, ΔT l'écart de température entre la paroi et le coeur et φ la densité de flux imposée en paroi.

Ensuite par définition nous pouvons calculer le nombre de Nusselt:

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

Où D est le diamètre du tube et λ la conductivité thermique du fluide.

Il est possible afin de valider le code NEPTUNE_CFD de comparer la correlation empirique de Gnielinski avec les résultats expériementaux.

La relation de Gnielinski valable pour: $0,5<Pr<2000$ et $2300<Re<5 000 000$ est:

$Nu=\frac{(f/8)(Re-1000)Pr}{1+12,7(f/8)^{1/2}(Pr^{2/3}-1)}$

Cette relation dépend du nombre de Reynolds. Nous avons donc lancé une campagne de tests paramétriques afin de comparer nos résultats avec cette relation.

                                   

Dans le tableau ci dessus nous avons recueilli les informations de cette étude paramétrique. Ensuite nous avons calculé les nombre de Nusselt avec la corrélation de Gnielinski et nous avons mis les résultats sous forme de graphe que voici.

    

Nous remarquons que nous obtenons bien le même ordre de grandeur. De plus pour les deux premier points et pour le dernier point il y a une bonne concordance. Par contre il y a une grosse différence entre la théorie et les résultats pour les deux points centraux.

Nous n'avons pas tout à fait compris cette erreur mais il se peut que l'erreur provienne du fait que nous n'avons pas tout à fait attendu l'état stationnaire pour ces deux tests (car pris par un manque de temps pour avoir des résultats).

En conclusion de cette partie nous pouvons dire que le code est validé par la corrélation de Gnielinski.

 

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.