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.