2. Simulations OpenFOAM

Dans cette partie, nous allons simuler l'écoulement de convection naturelle autour du cylindre à l'aide du logiciel OpenFOAM. Comme nous allons le voir plus en détails par la suite, le phénomène de convection naturelle étudié ici va imposer certaines contraintes au niveau de la résolution numérique. Un des objectifs de l'étude est de comparer les résultats obtenus pour différents modèles de turbulence : standard k-$\epsilon$, realizable  k-$\epsilon$, k-$\omega$, k-$\omega$  SST (Shear Stress Transport) et Spalart Allmaras.

 

Considérations numériques

  • Une des difficultés rencontrées lorsque l'on souhaite simuler une configuration de type convection naturelle réside dans la définition de conditions aux limites entrée/sortie. En fait, le fluide, initialement au repos, est mis en mouvement par l'apparition de gradients de densité, eux-mêmes dus aux variations de température. L'air chauffé par le cylindre monte, puis se refroidit , faisant apparaître un courant descendant car plus dense que l'air chaud. Imaginons que l'on fixe une condition d'entrée au niveau de la frontière inférieure et de sortie au niveau de la frontière supérieure, en imposant seulement une valeur de pression. Du fait de l'apparition de rouleaux de convection dans le domaine, les 2 frontières considérées voient à la fois du fluide sortir et rentrer simultanément. D'un point de vue numérique, ceci n'est pas acceptable pour une grande partie des codes de calcul CFD, y compris OpenFOAM ou du moins nous n'avons pas réussi à trouver les conditions aux limites adéquates. Nous sommes donc obligés d'isoler complètement le domaine, en imposant une condition de type wall sur toutes les frontières.

 

  • Le fait d'isoler complètement le cylindre va avoir une conséquence sur le choix du solveur sous OpenFOAM. En effet, nous avons le choix pour ce problème thermique entre 4 solveurs :
    • buoyantBoussinesqPimpleFoam
    • buoyantBoussinesqSimpleFoam
    • buoyantPimpleFoam
    • buoyantSimpleFoam

Les deux premiers solveurs (Boussinesq) fonctionnent en écoulement incompressible sous l'hypothèse de Boussinesq, tandis que les deux derniers fonctionnent en écoulement compressible. Étant données les hypothèses formulées précédemment, nous devons choisir parmi les deux solveurs Boussinesq.

Pimple se réfère à l'étude du régime transitoire, Simple le régime permanent. Comme nous considérons une géométrie complètement fermée, le cylindre va chauffer continuellement la pièce jusqu'à ce qu'à atteindre un état permanent où la température du domaine sera uniforme et égale à celle du cylindre. Ce n'est pas ce que nous souhaitons observer a priori, nous écartons donc également le solveur en régime permanent.

Nous allons donc effectuer nos simulations avec le solveur buoyantBoussinesqPimpleFoam, le but étant d'étendre la géométrie au maximum (attention au temps de calcul cependant) et de constater l'établissement de la couche limite autour du cylindre ainsi que le panache qui s'en dégage, avant les que re-circulations ne perturbent l'ensemble.

 

Maillage de la géométrie

La géométrie étant assez simple à réaliser, nous avons utilisé l'outil blockMesh pour ce problème également. Afin de simuler correctement la couche limite thermique proche du cylindre, nous raffinons le maillage dans cette zone. La gravité est selon l'axe horizontal et orientée vers la gauche (le fluide chaud "monte" vers la droite). Nous avons ici, ainsi que pour le second cas-test pour OpenFOAM, seulement créé la moitié de la géométrie afin de simplifier l'écriture des coordonnées dans le fichier blockMeshDict.

 

Maillage du domaine avec raffinement autour de la paroi du cylindre

 

Exploitation des résultats

 

Visualisation des résultats

Comme énoncé précédemment, le solveur que nous utilisons simule le régime transitoire. Nous avons donc dû estimer, pour une simulation donnée, le temps à partir duquel le régime permanent est atteint. Cela est réalisé à l'aide de Paraview : en traçant l'évolution du nombre de Nusselt le long du cylindre et en faisant varier le temps de simulation, on peut observer que la courbe se stabilise pendant un certain temps, avant de subir des changements sous l'effet des re-circulations provoquées par le fermeture complète du domaine. Ce laps de temps pendant lequel la courbe se stabilise nous donne alors une assez bonne approximation du régime permanent au voisinage du cylindre (seule l'évolution du nombre de Nusselt le long du cylindre nous intéresse).

Nous avons utilisé plusieurs valeurs de gravité ($g=10  m/s^2$, $g=100  m/s^2$, $g=1000  m/s^2$) afin de faire varier le nombre de Rayleigh. Parallèlement, plusieurs modèles de turbulence sont utilisés.

 

Résultats obtenus avec le modèle de turbulence $k-\omega$ :

 

Champs de température

 

 

 

 

 

 

 

 

Résultat de simulation pour le modèle k-$\omega$ et pour g=10 m/s2

 

Résultat de simulation pour le modèle k-$\omega$ et pour g=100 m/s2

 

Résultat de simulation pour le modèle k-$\omega$ et pour g=1000 m/s2

 

Champs de vitesse

Résultat de simulation pour g=10 m/s2

 

Résultat de simulation pour g=100 m/s2

 

Résultat de simulation pour g=1000 m/s2

 

L'augmentation de la gravité, et donc du nombre de Rayleigh, s'accompagne :

  • d'une réduction de l'épaisseur de la couche limite thermique autour du cylindre : les gradients thermiques sont plus élevés
  • d'une augmentation de la vitesse du fluide particulièrement visible au coeur du panache (les échelles de vitesses ne sont pas les mêmes sur chaque image). Pour la dernière valeur de gravité, le fluide atteint une vitesse de plus d'un mètre par seconde, alors que le cylindre et le mur de droite (correspondant à la frontière supérieure en réalité) ne sont séparés que de 15 cm. Il faut donc être très précis quant au moment à choisir pour l'estimation du régime permanent.

Des simulations pour des plus grandes valeurs de gravité ($g=10  000  m/s^2$ et $g=100  000  m/s^2$) ont été réalisées afin d'étendre encore plus l'intervalle de Rayleigh. Seulement, et cela rejoint la remarque faite précédemment, le domaine n'est plus assez étendu, la solution étant vite contaminée par les re-circulations (cf l'évolution du Nusselt sur la figure qui suit). De ce fait, il n'a pas été possible d'étendre l'intervalle de nombre de Rayleigh, les calculs étant trop volumineux.

 

Calcul du nombre de Nusselt

Afin de pouvoir visualiser l'évolution du Nusselt local le long du cylindre, nous avons choisi d'utiliser l'outil de post-traitement wallHeatFlux. Cet outil permet, une fois la simulation effectuée, de calculer les flux thermiques et les Nusselt associés sur toutes les frontières de type wall. Cependant, il est pré-programmé pour ne fonctionner qu'avec les solveurs simulant le phénomène de combustion en écoulement compressible (par exemple : reactingFoam). Il a fallu donc l'adapter pour le rendre compatible au solveur buoyantBoussinesqPimpleFoam qui, rappelons le, fonctionne en incompressible.

Ce problème se résout notamment en changeant le modèle thermophysique de l'outil : celui-ci se doit d'être le même que celui utilisé par le solveur buoyantBoussinesqPimpleFoam. Une fois correctement re-compilé et exécuté, l'outil rajoute dans chaque répertoire de temps un nouveau champ, NusseltNumber, calculé le long de toutes les parois.

 

Évolution du nombre de Nusselt local

La figure qui suit montre l'évolution du nombre de Nusselt local autour du cylindre, pour différentes valeurs du nombre de Rayleigh associées aux valeurs de gravité décrites précédemment. Le modèle de turbulence utilisé est le modèle k-$\epsilon$.

Évolution du nombre de Nusselt en fonction de l'angle autour du cylindre en régime turbulent avec le modèle k-$\epsilon$

 

Les courbes correpondantes aux valeurs $Ra_D=4.10^3$, $4.10^4$ et $4.10^5$ montrent d'excellents résultats qualitativement parlant, au vu de la courbe théorique présentée dans la partie préliminaire. Comme prévu, les 2 dernières courbes ($Ra_D=4.10^6$ et $4.10^7$) s'écartent de la courbe théorique, mais peuvent nous donner un ordre de grandeur du Nusselt local sur la partie inférieure du cylindre ($ \theta = - \pi /2$), qui est moins perturbée par les re-circulations.

 

Calcul du nombre de Nusselt moyen

Nous nous intéressons maintenant au nombre de nusselt moyen, $Nu_D$, qui s'obtient en intégrant le nombre de Nusselt local le long du cylindre, tracé en fonction du nombre de Rayleigh pour les différents modèles de turbulence considérés. Ces résultats sont comparés à la corrélation de Churchill et Chu :

 

Évolution du nombre de Nusselt moyen en fonction du nombre de Rayleigh pour divers modèles de turbulence

On peut remarquer que pour les "petites" valeurs de $Ra_D$, c'est-à-dire pour $Ra_D<10^5$, les valeurs de $Nu_D$ sont très proches de la courbe théorique. Cela est particulièrement visible pour les modèles k-$\omega$ et k-$\omega$ SST. Néanmoins, il semble que globalement, pour ces valeurs de $Ra_D$, le nombre de Nusselt moyen soit légèrement sur-estimé, ce qui se traduit par une sur-estimation du flux dégagé par le cylindre.

Lorsque le nombre de Rayleigh augmente encore, la tendance inverse est observée (sous-estimation du Nusselt moyen). Les modèles standard  k-$\epsilon$ et Spalart Allmaras, qui donnent l'impression de globalement sur-estimer le nombre de Nusselt moyen pour les faibles valeurs du nombre de Rayleigh comparativement aux autres modèles, sont ceux qui donnent les meilleurs résultats à haut nombre de Rayleigh.