Simulation du refroidissement extérieur d’un REP à l'aide du code Neptune

         

ZMELTY Aurélien & LIGEOIS Thomas

Encadré par Catherine COLIN  & Hervé NEAU
 


L'objectif de ce BEI est de simuler l'écoulement d'eau le long des parois externes d'un REP en cas de réchauffement brutal de la structure dû à un accident nucléaire.

Le but est de regarder l'évolution du flux évacué par la paroi lorsqu'une crise d'ébullition se produit dans les conditions atmosphériques.


Introduction

Le contexte

Le Bureau d'Etude Industrielle (BEI) est un projet s’inscrivant dans l'enseignement de l'ENSEEIHT. Avant la période de stage, ce projet de deux mois permet de mettre l’étudiant en suite d’ingénieur afin de résoudre une problématique réelle proposée par des industriels. L’étudiant (associé à un binôme) s'investit pendant six semaines à plein temps, avec au préalable une période d'imprégnation du sujet qui se résume à sept séances de deux heures (entre novembre et janvier). Dans ce sens, des réunions avec l’encadrant, et éventuellement le contact industriel ont lieu pour préciser le cadre de l'étude et la définition des objectifs.

Ce sujet a été proposé par EDF. Dès lors, ce projet sera réalisé avec la collaboration de Monsieur MINOUNI, ingénieur recherche chez EDF sur le site de CHATOU, et sera encadré par Madame Catherine COLIN et Monsieur Hervé NEAU, tous les deux ingénieurs chercheurs à l'IMFT et professeurs à l'INPT. En ce qui concerne, les étudiants chargés de l'étude, il s'agit de Monsieur Thomas LIGEOIS (élève en mécanique des fluides numériques) et Monsieur Aurélien ZMELTY (élève énergétique et procédés).

Le groupe EDF par son statut de premier électricien nucléaire mondial, se doit d'assurer la pérennité de ses ouvrages. Le parc nucléaire français est le plus fourni d'Europe et il réside des risques inhérents au bon fonctionnement de la centrale, parmi lesquels la crise d'ébullition dans le circuit primaire. En effet, ce phénomène déplorable détériore la surface des barreaux de combustible. Le but de cette étude serait donc d'améliorer les modèles locaux actuels de prédiction de ce phénomène.

 

La mission

Les réacteurs à eau pressurisé constituent l'essentiel du parc nucléaire français. Ainsi, l'énergie issue de la réaction sert à vaporiser de l'eau qui entraînera les turbines servant à produire l'électricité. La vaporisation de l'eau se fait au sein des générateurs de vapeur. Dans le circuit primaire, l'eau pressurisée n'est pas vaporisée dans des conditions normales de fonctionnement. Cependant, en cas de dépression brutale, ou d'augmentation ponctuelle de température une ébullition peut se déclencher sur la paroi. En effet, une ébullition nucléée peut évoluer en une ébullition en film. Dans ces conditions, les coefficients d'échange chutent brutalement. L'énergie contenue dans les barreaux n'est plus dissipée par le fluide; la température des parois augmente donc de façon drastique. Cet échauffement ponctuel risque d'endommager lesdites parois et de compromettre la longévité du matériel. Il est donc nécessaire de bien modéliser les transferts thermiques en proche paroi.  Ainsi, cette étude s'inscrit dans une problématique sécuritaire au niveau du fonctionnement du cycle primaire. 

Le projet et son organisation

Quelques rappels préalables

Le cas des REP

Les réacteurs à eau pressurisée (REP) ou bouillante (REB) constituent l'essentiel du parc nucléaire français. Deux circuits d'eau fermés (primaires et secondaires) et un circuit d'eau ouvert (circuit de refroidissement) sont nécessaires au fonctionnement d'un REP.

Ainsi, l'énergie issue de la réaction exothermique produite par la fission du combustible nucléaire situé dans le cœur du réacteur sert à vaporiser de l'eau à l'aide d'une source de chaleur (barreaux contenant les matériaux radioactifs) qui entraînera les turbines. Ces dernières feront tourner les pâles d'un alternateur, lequel produira de l'électricité. Enfin, la vaporisation de l'eau se fait au sein des générateurs de vapeur. 

Dans le cas où les barreaux se mettraient à chauffer anormalement plusieurs dispositifs sont prévus notamment l’ouverture de la vanne du condenseur. En effet cela  permet d’évacuer l’air chaud sous haute pression vers les cuves remplies d’eau du condenseur et qui à leurs tours transmettront l’eau froide vers la cuve du réacteur. Une deuxième solution en complément de la première, est de baigner la cuve du réacteur dans de l’eau et de la mettre en mouvement afin d’extraire et d’évacuer la chaleur. C’est sur ce sujet que repose notre étude.

 

Comment cela fonctionne-t-il ?

L’idée est de regarder dans le cas d’accidents nucléaires type Fukushima où se trouve une surchauffe importante des réservoirs des réacteurs.

Il est prévu de faire un arrosage extérieur pour refroidir par convection le réservoir par l’extérieur (ce qui n’a pas été fait à temps à Fukushima car les personnes en charge voulaient préserver la centrale un maximum de temps ; d’où une certaine hésitation à mettre de l’eau de mer au départ car cela compromettait l’existence de la centrale). A cause de ce retard, la température a commencé à augmenter et son état a atteint un point d’irréversibilité avec des dégagements d’hydrogène et des conséquences plus importantes par la suite.

Il faut donc réfléchir à la manière la plus efficace pour asperger de l’eau par un circuit extérieur autour de l’enveloppe du réacteur nucléaire pour refroidir au maximum à l’intérieur du réacteur (on regarde la fusion du cœur par exemple (accident de Tchernobyl)).

L'étude porte sur la partie du circuit primaire, et comme cela vient d’être dit sur le refroidissement extérieur du cœur du réacteur.

Il existe deux méthodes de refroidissement :

- Refroidissement par jet impactant : cette méthode est probablement moins efficace que la seconde car le cœur étant tellement chaud les gouttelettes ne touchent pas la paroi car elles sont vaporisées avant, on est au-delà du flux critique. L’augmentation de la puissance du jet pourrait être vue comme une bonne solution mais cela pourrait endommager les parois du réacteur.

- Un re-noyage par l’extérieur (écoulement d’eau qui va circuler dans une enveloppe à l’extérieur du cœur) semble la méthode la mieux adaptée. Nous sommes alors dans une configuration d’un écoulement dans un canal avec une ébullition nucléée sur la paroi.

Dans le cas des régimes d'ébullition en paroi, les lois qui sont mises en place correspondent à des régimes hautes pressions fortement sous refroidi à une température de 320°C. On est donc 40 à 50°C en-dessous de la température de saturation justement pour éviter qu'il y ait une forte ébullition dans ces circuits primaires de centrale (l’eau ne doit pas bouillir en fonctionnement normal). En effet les flux de chaleurs sont extrêmement forts, si l'ébullition commence à se déclencher, le risque est le développement d'une vaporisation très forte sur la paroi. Cette vaporisation créée alors une isolation entre l'écoulement liquide et les parois qui ne sont plus sous-refroidis, d'où une surchauffe importante.

Études menées et perspectives

De nombreuses études ont été menées, notamment la thèse de Michael Montou avec l’étude des modèles d’ébullition en paroi, puis les rentrer dans le code de simulation Neptune afin de tester ces modèles et améliorer ceux qu’il y a dans le code pour finalement s’approcher du flux critique et regarder les conditions.

Concernant notre étude, nous serons placés plus en amont, à savoir la recherche du bon modèle d’ébullition qui traduise ce qu’il se passe à pression atmosphérique avec de l’eau.

Donc la première partie de l’étude va consister à chercher qu’est-ce qu’on a comme résultats expérimentaux à notre disposition sur des écoulements bouillants avec de l’eau à pression atmosphérique (étape de recherche bibliographique et surtout ces expériences permettront (à pression atmosphérique) de regarder ce que donnera le code Neptune et de faire la validation).

 

Nous sommes un groupe de quatre et la séparation en deux binômes permettra de traiter séparément :

-       La Modélisation des flux de chaleur en paroi (Recherche des données, validation et comparaison expérimentale).

-       La Simulation numérique avec Neptune (groupe de Thomas LIGEOIS et d’Aurélien ZMELTY).

 

Cependant les deux approches sont assez complémentaires donc les deux groupes seront amenés à travailler ensemble (échange de données).

Malgré cela, notre groupe devra s’attacher à la compréhension du cours de diphasique, à savoir les modèles d’ébullition nucléées en paroi (notamment résultats issus des modèles de Yu et Wang). Attention tout de même car ici il s’agit d’une ébullition convective et les lois de fermeture sur les fréquences et les diamètres sont différentes du cas classique. Il va y avoir  de plus sur les aspects de modélisation, il y a des choses intéressantes à faire comme les diamètres de détachement des bulles, les fréquences de déformation, (ce n’est pas assez précis et adapté).

Une étudiante a regardé l’an passé la modification des modèles en paroi qu’il y a dans Neptune et  elle a fait un codage sous Matlab. Mais aucune comparaison avec des bases de données expérimentales. Il peut donc être intéressant de voir ce de voir ce qu’il y a dans Neptune. Dans ce cadre nous allons devoir nous auto former sur ce logiciel ainsi que le logiciel de maillage qu’il lui est associé : Simail.

De plus la collaboration avec Stéphane MINOUNI va s’avérer nécessaire en raison de son expertise dans ce domaine et pour l’échange de de maillage (maillage EDF). De même l’étude approfondie de  la thèse de Michael MONTOU va être nécessaire pour comprendre tous les aspects importants de ce sujet.

Nous lui demanderons également les résultats des simulations des écoulements en tube ainsi que l’envoi d’un cas test bouillant qui fonctionne correctement, donc il suffit de le lancer pour voir les résultats.

Démarche adoptée

Afin d’aborder ce problème la démarche a été d’étudier trois grands cas à pression atmosphérique :

  • « MONO_FROID » : cas monophasique où l’écoulement sur la paroi est à température de de 20°C et la paroi n’est pas chauffée.
  • « MONO_CHAUD » : cas monophasique avec une température en paroi inférieure à la température d’ébullition.
  • « EBULLITION » : cas diphasique avec une température en paroi supérieure à la température d’ébullition.

Les deux premiers cas vont permettre également une adaptation progressive avec le logiciel NEPTUNE_CFD et permettront de voir l’influence du maillage sur l’étude. En effet celui-ci à un rôle crucial en fonction du cas à étudier et va être amené à évoluer tout au long de l’étude.

Dans un second temps nous nous attacherons à la modélisation des transferts thermiques en proche paroi (interaction avec l’autre groupe et recherches propre au nôtre). Les deux étapes seront :

  • La modélisation des flux de chaleur qu’on va avoir dans un extérieur d’un réacteur à eau pressurisée.
  • Code Neptune : pour le moment ce-dernier n’est pas utilisé en simulation pour l’extérieur des cuves mais pour l’intérieur du réacteur ; il est utilisé et il validé pour des écoulements d’eau à haute pression (155 bars)  et haute température (320°C). La simulation sera donc réalisée à basse pression en tenant compte des modèles de flux de chaleur.

Les étapes clés :

  • Angle d’attaque du problème : données du problème
  • Initiation à Simail
  • Initiation à Neptune
  • Analyse des modèles les plus utilisés dans les différents codes de calculs (Neptune, CFX, Fluent, StarCCM+)
  • Proposition d’amélioration de ces modèles
  • Implémentation dans le code Neptune
  • Simulation numérique avec Neptune bouillant en tube et comparaison avec les résultats expérimentaux issus de la littérature et de la thèse de Michael MONTOUT.
  • Traitement de texte d’un rapport écrit.
  • Préparation de diapositive sous PowerPoint.

Organisation des tâches

Production de document et conditions de distribution

Environnement

Pour la transmission électronique des documents, le système adopté doit accepter les noms de fichiers long (ex:Windows98 ou Linux)

Logiciels informatiques

Afin de ne pas engendrer de conflits de versions informatiques, il est impératif de respecter les spécifications suivantes:

  • Pour le traitement de texte: Microsoft Word97
  • Pour la Planification de Projet: Microsoft Projet 98
  • Pour les informations tabulaires et graphiques: Microsoft Excel97
  • Pour les présentations: Microsoft PowerPoint97
  • Pour les mailler l'environnement : Simail
  • Pour analyser le problème diphasique : Neptune CFD
  • Pour traiter les données graphiques : Paraview
  • Pour traiter les données numériques : Matlab R2013a

Distribution Electronique

Les documents peuvent être échangés et distribués au moyen de systèmes d'email employés pour des transferts de document. (norme de la plupart des logiciels: (MS Exchange, MS Outlook, Netscape/Mozilla).

Diagramme de Gant et management des risques

Diagramme de Gant

 

Management des risques

Malgré la simplicité du projet, il entraîne quelques risques. Le tableau suivant résume les risques prévus et les actions à entreprendre :

Les maillages

Tout au long du BEI il a fallut adapter le maillage au problème. Dans ce chapitre nous allons voir comment le maillage a évolué au fur et à mesure. Nous proposons un tutoriel de tous les maillages conçu pour ce BEI à l'aide du logiciel Simail.

Le maillage cylindrique

Pour simuler un écoulement dans un tube la première solution qui vient à l'esprit est de faire le maillage complet du tube. On a donc le maillage que l'on appelle tube_3D.des avec le logiciel simail. Pour faire ce maillage il faut commencer par faire un O'Grid en 2D. Le maillage O'Grid est un maillage 2D définissant un cercle maillé de manière à ce que l'on puisse bien raffiner les maillages de paroi, c'est donc celui qui sera retenu. Un maillage comportant des maillages hexaédriques dans une forme de cercle n'est pas chose simple à réaliser. Pour ce faire les étapes d'acheminement du maillage complet sont explicitées ci-dessous:

1) Placer les points qui serviront de base pour notre cercle. Pour cela il faut placer 7 points. Nous avons placé un points au centre qui n'a pas grand intérêt sauf pour se repérer et être sûr de soi. Le point 1 correspond donc au centre du cercle. Les points 2 et 3 sont les points se trouvant sur le cercle dans 2 directions orthogonales (Ici selon Ox et Oy). Ensuite les points 4 et 5 se situent à 3/4 du rayon pour des raisons qui seront évoquées plus tard. Enfin les points 6 et 7 servent pour finir le carré central avec les points 4 et 5. Pour placer les points avec le logiciel simail il faut aller dans l'onglet GEOMETRIE puis POINTS. Ensuite il faut rentrer les points un par un en mettant le numéro correspondant mais aucun label et la référence doit être 0.

Voici le tableau récapitulatif de simail pour le positionnement des points ainsi que la représentation graphique :

2) Une fois les points placés il faut tracer les lignes. Il y a 7 lignes à tracer. La ligne 1 rejoint les points 2 et 3 pour former le premier quart de cercle. Cette ligne étant une paroi il faut la définir comme telle, nous retiendrons pour la suite la référence 50 pour les parois. Puis les lignes 3 et 4 qui relient le carré central avec le cercle. Puis il y a les lignes 2, 5, 6 et 7 qui forment le carré central. Une raison de 10 a été prise pour les lignes 3 et 4 afin de raffiner près de la paroi. Toutes les lignes ont la référence 0 qui correspond au fluide.

Voici le tableau récapitulatif de Simail pour le positionnement des lignes ainsi que la représentation graphique :

3) Après les lignes il est temps de passer au maillage. Il faudra faire les 2 maillages séparément puis les recoller ultérieurement. Le maillage 1 est celui du carré central, quand au 2 il s'agit du maillage correspondant au quart de cercle. Pour créer un maillage il faut aller dans l'onglet MAILLEURS puis SURFACIQUES. On choisit un maillage en QUADRULATION COURBE et BIDIMENSIONNEL. Pour le sous-domaine il faut sélectionner le nombre de sous-domaines que comportera notre maillage, ici 1. Pour le nombre de lignes du côté 1 il faut entrer le nombre maximum de lignes sur le même côté, ici 1. Enfin pour les lignes du contour il faut entrer les lignes qui entourent le domaine et de manière à ce qu'elles se suivent.

Voici la fenêtre Simail pour le maillage 1 ainsi que la représentation graphique des maillages 1 et 2 :

4) Maintenant que l'on a le maillage de base du quart de cercle on va le faire tourner 3 fois. De sorte que les maillages 3, 4 et 5 représentent les 3 autres quarts de cercle qui sont tournés respectivement de 90°, 180° et 270°. Pour cela il faut aller dans l'onglet TRANSFORMATION et GEOMETRIQUES et choisir la ROTATION. Dans cette partie il faut choisir le maillage d'entrée, ici 1. Les maillages de sortie seront 3, 4 ou 5. Le point invariant sera toujours 0,0,0.

Voici ce qui est affiché dans simail pour une rotation dans le cas du maillage 3 :

5) L'étape suivante est le recollage des maillages. En effet nous avons désormais 5 maillages mais il pour simplifier la suite nous allons travailler avec un seul maillage. Le maillage 6 est donc un recollage des maillages 1, 2, 3, 4 et 5. Pour cela il faut aller dans l'onglet TRANSFORMATION puis AUTRES puis dans RECOLLER et entrer les maillages d'entrée et celui de sortie. Le paramètre de précision est purement indicatif dans le cas de maillages qui se collent parfaitement mais celui joue un rôle essentiel dans le cas où les mailles des maillages que l'on veut recoller ne coïncident pas. Pour ne pas mettre 0 on mettra 0.001 mais cela n'aura pas d'influence à priori.

Voici ce qu'affiche Simail pour le recollage ainsi que le maillage 6 dans la fenêtre graphique :

6) Finalement on remarque qu'entre le carré central et les quart de cercle le maillage n'est pas très régulier ce qui peut poser des problèmes de calcul. Pour avoir un maillage plus régulier on procède à l'étape de régularisation qui va courber les lignes pour homogénéiser le maillage. Pour cela il faut aller dans l'onglet TRANSFORMATION puis AUTRES puis REGULARISER. On va donc régulariser le maillage 6 en 7 en faisant entre 15 et 20 itérations car en deçà le maillage n'est pas assez régulariser et au-delà on va trop modifier le maillage de base.

Voici ce qui est affiché dans simail pour cette étape ainsi que la fenêtre graphique du maillage 7 régularisé :

Une fois que le maillage O'Grid 2D est terminé il suffit de le passer en 3D. Cette dernière étape s'effectue en allant dans l'onglet MAILLEURS puis VOLUMIQUES puis ELEVATION. Après avoir choisi les maillages d'entrée et de sortie il faut sélectionner la montée verticale. Dans cette onglet on doit préciser, la côte de la première section, ici 0 car la géométrie impose que le bas du cylindre soit à la côte 0. La côte de la dernière section impose le haut du cylindre dans c'est indirectement la hauteur du cylindre, ici 1. La subdivision des éléments 2D doit toujours rester à 1 car on ne veut pas modifier notre maillage de base. On choisit de rester avec des quadrangles sur la hauteur et 500 sections équidistantes. Les références des faces correspondent aux références que l'on veut pour le bas et le haut du cylindre. Ici comme on veut que le fluide entre par le bas et sorte par le haut on impose les références 10 pour l'entrée et 90 pour la sortie. Enfin la transmission des références va définir les références des nouvelles faces créent. Donc comme ici on va créer uniquement des parois latérales alors on va choisir 50 50 pour signifier que tout ce qui est en référence 50 doit le rester tout le long du tube.

Voici ce qui s'affiche avec simail pour l'élévation verticale ainsi que la fenêtre graphique qui représente le maillage 3D complet :

Quelques chiffres pour ce maillage:

Il y a 20 mailles sur le quart de cercle (il y en a aussi 20 par côté du carré central). Or comme il y a 10 mailles sur le petit côté de longueur 1/4 rayon situé entre le carré et le cercle cela signifie qu'il y a 200 cellules par quart de cercle. De plus comme le carré comporte aussi 20 mailles de côté il y a donc 400 cellules dans le carré central. Cela nous donne donc 4 fois 200 plus 400 soit 1200 cellules dans la base 2D. Étant donné que l'élévation comporte 500 mailles alors on a un maillage 3D de 600 000 mailles.

Les pas d'espace $\Delta X$ et  $\Delta Y$ sont du même ordre de grandeur. Avec  $\Delta X= \frac{1}{20}$ et $l=\frac{2\pi R}{4}$,  la longueur du quart d'arc de cercle. Donc $\Delta X$ = 8.10-4 m et $\frac{\Delta Y}{\Delta X}$ = 1.

Le maillage partiel

Pour simplifier le cas du cylindre on va découper ce dernier en quartier. C'est un cas de symétrie axiale en 2D sur une maille d'épaisseur (typiquement un secteur angulaire compris entre 2 et 8° d'angle). On va donc avoir deux stratégies possibles pour faire un tel maillage.

1) On créer la forme de la base et on l'élève comme pour le cylindre.

2) On créer un rectangle que l'on va faire tourner autour d'un axe.

On gardera la seconde stratégie car elle est plus simple donc plus rapide. Voici le détail des opérations pour réaliser ce maillage:

  • La première étape de la conception du maillage consiste à faire un rectangle. Pour cela on va faire 4 points et les relier par des lignes. Nous nous plaçons de telle sorte que le point en bas à gauche du rectangle soit placé à l'origine et soit le point 1. Ensuite les points tournent dans le sens polaire. De même pour les lignes, on commence par la ligne 1 entre les points 1 et 2 puis on tourne dans les sens polaire pour les lignes 2, 3 et 4. Pour les références, on choisit que le fluide entre par le bas, sorte par le haut. Sur la droite on a un mur et sur la gauche une condition de symétrie. Donc on a 10 pour la ligne du bas (1), 50 pour la ligne de droite (2), 90 pour la ligne du haut (3) et 41 pour la ligne de gauche (4). Nous avons imposé une raison de 0.5 sur la ligne 1 et de 2 pour la ligne 3 (inverser car l'une va de gauche à droite et l'autre l'inverse) afin de raffiner proche paroi.

Voici ce que l'on observe avec Simail pour les points et les lignes :

Voici ce que l'on observe dans la fenêtre graphique :

Le tube étant filiforme on a du mal à observer mieux avec la fenêtre graphique. Pour la suite les images seront zoomées sur la base du tube.

  • Une fois que les lignes sont faites il faut faire le maillage exactement comme pour le O'Grid 2D du cylindre. Les caractéristiques sont données ci-dessous :

Voici le maillage vu avec la fenêtre graphique de Simail :

  • Avant de faire la rotation il faut s'assurer que le maillage soit à une certaine distance de l'axe de symétrie. Cette remarque vient du fait que si on fait tourner le maillage comme tel avec comme axe de symétrie le côté gauche du maillage alors on va se retrouver avec un maille centrale triangulaire. Ce qui n'est pas bon numériquement. EDF nous propose de décaler le maillage d'une distance égale à L/5 où L est la taille de la maille centrale. Il faut donc prendre cette valeur et la rajouter aux points 1 et 4 selon l'axe Ox.

Remarque 1 : Il ne faut pas décaler tout le maillage sinon le rayon de notre tube se verrai augmenter ce qui poserais des problèmes quant à la conservation du débit pour comparer par la suite avec les expériences.

Remarque 2 : Pour pouvoir gagner en efficacité cette étape peut être réalisé en fin de compilation du maillage dans le fichier .dat. En modifiant directement les coordonnées des points plutôt que de devoir revalider toutes les étapes du maillage une à une.

  • Le maillage étant décalé on peut désormais effectuer la rotation. Il faut au préalable calculer l'angle de rotation en gardant en tête que l'on désire garder une maille la plus cubique possible. Pour cela l'angle est donné par la relation suivante:

$\theta = \frac{360 \Delta X}{2 \pi R}$

avec $\theta$  en degrés d'angle, ΔX la taille de la plus petite maille et R le rayon du cylindre. Dans notre cas  et  ce qui nous donne un angle de 4.4° On choisira donc un angle de 5° pour avoir un rapport de 1.1 entre ΔX et ΔZ.

Voici ce qui est entré dans SIMAIL pour effectuer la rotation :

Le changement d'origine n'ayant pas lieu il faut tout de même mettre 0 partout. Pour l'angle Ox/axe de rotation il s'agit de l'angle entre le plan qui contient l'axe de rotation et l'axe Ox. Ici comme c'est l'axe Oy qui fait office d'axe de rotation il y a 90° entre les deux. Puis on entre l'angle calculé précédemment et enfin un point invariant se trouvant donc sur l'axe de rotation hormis l'origine. Pour les transmissions de références c'est comme pour le cas précédent. Il faut garder les références de toutes les lignes qui ont été construites aux préalables. Pour le nombre de sections on n'en choisira que 2 car il n'y a qu'une seule maille d'épaisseur.

Voici ce que l'on peut observer dans la fenêtre graphique une fois la rotation effectuée :

-Le problème est que contrairement au cas de l'élévation, ici on ne peut pas choisir la référence des faces que l'on a fait tourner. Pour cela il y a une astuce. Dans l'onglet TRANSFORMATION puis ATTRIBUTS on peut sélectionner l'onglet REFERENCE dans NUMERO. Après avoir indiqué que l'on passe du maillage 2 en 3 on entre les références que l'on désire changer. Comme l'intérieur du maillage 2D était du fluide toute les faces de bord ont la référence 1. Il faut donc passer les références de bord de 1 à 40 et toutes les cellules des deux faces créées seront des symétries.

Voici ce qu'affiche Simail pour cette étape :

On peut s'assurer que l'opération à bien été effectué en regardant le changement sur le maillage dans la fenêtre graphique :

Les faces sont vert clair ce qui représente une condition de symétrie.

Quelques chiffres pour ce maillage :

Il y a 10 mailles sur le rayon soit un ΔX égal à 2,7.10-4 m. Sur la hauteur on a choisit de mettre 3000 mailles ce qui donne un  ΔY =  3,2.10-4 m soit ΔY /ΔX = 1,23 . Pour la troisième dimension on a déjà calculé qu'avec un angle de 5°, on trouve un rapport ΔZ /ΔX = 1,1 .

Le maillage comporte donc en tout 30 000 cellules ce qui confère un gain de temps de calcul conséquent par rapport au cylindre complet.

Le maillage partiel raffiné

Les résultats obtenus avec le maillage partiel sont beaucoup plus rapide qu'avec le maillage cylindrique. De plus la précision des résultats est assez bonne sachant que l'on a eu l'occasion de faire de très petits pas de temps du au maille fine sur la hauteur du maillage. Malgré cela, nous pouvons tout de même remarquer que le maillage reste grossier près de la paroi. Le solution apportée est de raffiner le maillage dans la couche limite. Nous avons calculé que la couche limite coïncidait avec la troisième maille en partant de la paroi. Le choix a été de créer un sous domaine dans ses trois mailles et de raffiner précisément dans ce domaine alors que l'on peut se permettre de relâcher dans l'autre domaine loin de la paroi. Les explications détaillées de la construction du nouveau maillage sont ci-dessous:

  • la première étape est la construction de 2 nouveaux points. Un sur la ligne d'entrée et un sur la ligne de sortie tous deux situés au niveau de la troisième maille en partant de la paroi. Voici le récapitulatif de Simail pour les points:

                                 

  • Ensuite il faut absolument refaire les lignes car il doit y avoir 6 lignes au lieu de 4. il y en a désormais 2 en bas et 2 en haut. On introduit donc une raison de 0,25 en bas et 4 en haut sur la petite ligne afin d'avoir un fort raffinage proche paroi. Attention il faut prendre en compte le calcul de la nouvel taille de la petite maille et donc recalculer le nombre de mailles sur la hauteur. Voici le récapitulatif des lignes donné par SIMAIL:

                                 

  • Puis on procède de la même façon pour les autres maillages. C'est à dire un maillage bidimensionnel, la différence est que celui-ci doit contenir 2 sous-domaines et que la ligne qui comporte le plus de découpage en a 2 et non plus 1. Voici ce qu'affiche SIMAIL pour la création du maillage:

                                 

  • Une fois le maillage crée, on va pourvoir faire la rotation comme pour le maillage partiel et ajouter en attribut le changement de référence pour les faces de symétrie afin de les passer de 1 à 40. Voici ce qu'on observe en fenêtre graphique pour le maillage final:

                                 

Quelques chiffres pour ce maillage:

Il y a 17 mailles sur le rayon soit un un pas d'espace selon X, ΔX = 5,52.10-5 m . Sur la hauteur on a choisit de mettre 7000 mailles ce qui donne un pas d'espace suivant Y, ΔY =  1,43.10-4 m, $\frac{\Delta Y}{\Delta X}$ =2,59. Pour la troisième dimension on choisit un angle de 2.7° ce qui donne un  ΔZ = 1,65.10-4 m, on trouve alors un rapport $\frac{\Delta Z}{\Delta X}$ = 2,98. De cette manière on a un maillage qui satisfait les conditions de bon calcul des gradients normaux à l'écoulement avec des rapports inférieurs à 3.

Le maillage comporte donc en tout 119 000 cellules ce qui le place en meilleur maillage pour les temps de calcul que le maillage cylindrique. Et il devient beaucoup plus précis que le maillage partiel de base et donc pourra à priori capter de manière plus précise les phénomènes physiques qui se produisent proche paroi pour la suite.

 

 

Le maillage partiel 2

Le problème du maillage précédent est qu'il était trop raffiné. Après une étude et des tests sur des simulations, nous avons remarqué qu'il y avait les premières mailles se trouvant dans la sous couche visqueuse. Cela a eu pour effet de sur estimer la valeur de la vitesse dans tout le domaine. Pour avoir un compromis entre raffinement du maillage et cohérence des résultats on décide de placer la maille pariétal à y+ = 20. Sachant que les limites de validation de la loi Log commencent à  et s'étendent jusqu'à y+ = 200.

La maille pariétale est donc maintenant plus grossière que le maillage raffiné mais reste tout de même plus raffiné que le maillage partiel du début. C'est donc le maillage le plus réfléchi et adapté pour notre étude. Voici le récapitulatif des changements apportés au maillage partiel2.des.

Pour les points voici le récapitulatif des coordonnées affichées par SIMAIL:

                              

On remarque que les points 1 et 4 sont déjà décalés de ce qu'il faut (L/5) par rapport à 0 pour pouvoir effectuer la rotation. Ensuite on trace les mêmes lignes que pour le premier maillage partiel avec les conditions suivantes:

                              

Puis on crée le maillage 2D comme pour le maillage partiel, on fait tourner le maillage d'un angle de 5° et on attribue aux faces de symétrie la références 40. Voici ce que l'on peut observer avec la fenêtre graphique pour le maillage final:

                                

Quelques chiffres pour ce maillage:

Il y a 11 mailles sur le rayon soit un ΔX = 2,16.10-4 m. Sur la hauteur on a choisit de mettre 7000 mailles ce qui donne un  soit ΔY = 6,45.10-4 m, soit $\frac{\Delta Y}{\Delta X}$= 2,99. Pour la troisième dimension on choisit un angle de 5° ce qui donne ΔZ =3.10-4 m,   on trouve alors un rapport $\frac{\Delta Z}{\Delta X}$ = 1,39 .

De cette manière on a un maillage qui satisfait les conditions de bon calcul des gradients normaux à l'écoulement avec des rapports inférieurs à 3.

Le maillage comporte donc en tout 17 050 cellules ce qui le place en meilleur maillage que ce soit pour les temps de calcul ou pour la précision des résultats.

Pour la suite on effectuera tous les calculs avec ce maillage.

Remarque : le maillage partiel2.des a été également travaillé, conçu, de manière à avoir une la plus petite maille de l’ordre du diamètre des plus grosses bulles, à savoir un diamètre de 3,5.10-4 m (voir maillage dans la partie EBULLITION).

Voici un récapitulatif de l'évolution des maillages:

           

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.

Conclusion

—Conclusions sur l'étude:

Conclusions sur le plan personnel :

  • Expérience d’un nouveau code de calcul, NEPTUNE et d'un logiciel de maillage associé: SIMAIL.
  • Etude de la physique diphasique. Plus généralement, le sujet était assez large et a permis au binôme de revoir "ses gammes" de la mécanique des fluides enseignée à l'ENSEEIHT: transferts thermiques, écoulements laminaires et turbulents, mécanique des fluides numériques,  et découverte des échanges en diphasique.
  • Expérience d’un cas industriel concret et en relation directe avec des situations réelles.
  • Véritable expérience de travail en groupe, entre deux personnes d'un binôme n'ayant jamais travaillé dans le passé ensemble et aboutissant à l'élaboration d'un projet ainsi qu'à l'obtention de résultats.
  • Gain majeur en termes d’organisation des tâches : par exemple dans les résultats à noter directement après l’expérience ou après un calcul.
  • Méthodologie sur le projet : dé-zoomer fréquemment afin de se ressituer dans le contexte et afficher les résultats importants.
  • Gain également sur le plan des détections d’erreur : prendre le temps de lire les rapports d’erreur à mise en place de structures de vérifications et de contrôle dans les codes.
  • Bilan enfin positif sur les pics de travaux  à effectuer dans un laps très court.

Remerciements

Nous tenions à remercier les trois professeurs qui nous ont laissé de leur temps et leurs patience, à savoir, Hervé NEAU, Catherine COLIN et Renaud ANSART.
Merci pour vos coups de pouce et vos visites à l'ENSEEIHT sans lesquelles nous n'aurions pu mener ce projet correctement et efficacement. Vos conseils nous seront précieux pour la suite.

Thomas & Aurélien.

Références

  • Thèse de Michaël MONTOUT: Contribution au développement d'une Approche Prédictive Locale de la crise d'ébullition.
  • Thèse de Romain DENEFLE: Modélisation locale diphasique eau-vapeur des
    écoulements dans les générateurs de vapeur.
  • Thèse de Saliha SENHAJI: Simulation numérique de l’évaporation des alcools sous forme de
    film turbulent en présence d’un écoulement d’air dans un tube
    vertical.
  • Thèse de Catalin Viorel POPA: Étude théorique et expérimentale du comportement transitoire
    d’un écoulement laminaire en convection mixte dans un tube
    vertical.
  • Cours de Mécanique des fluides et transferts I; partie sur les écoulements turbulents: Gr´egoire WINCKELMANS et VincenT LEGAT.
  • Cours de thermique appliquée: Frédéric DOUMENC.
  • Cours d'introduction à la couche limite turbulente: ENSTA.
  • Cours d'écoulements diphasiques avec changement de phase, troisième année ENSEEIHT, parcours énergétique, Catherine COLIN.
  • Mise en oeuvre d'un cas de calcul NEPTUNE_CFD v1.08@Tlse à l'ENSEEIHT - RELEASE 2012, Hervé NEAU, Olivier SIMONIN.
  • Modèles de fermetures pour les simulations diphasiques, Renaud ANSART, Hervé NEAU.

Téléchargements

Ci-dessous, le lien vers:

  • les routines Matlab des 3 cas.
  • les fichiers paramètres de ces cas.
  • les fichiers utilisateurs à copier dans les répertoire USER des cas mono-chauds et ébullition.
  • le maillage partiel2.des et son fichier de création partiel2.dat pour modifier à sa guise.

/travaux/bei/beiep/sites/default/files/users/tligeois/fichier_projet_bei_ligeois_zmelty.zip