IV. Premier cas-test : convection naturelle autour d'un câble électrique

Ce premier cas-test traite de l'étude de la convection naturelle autour d'un câble électrique modélisé par un cylindre. Ce câble est soumis à un fort courant électrique et dissipe la chaleur qui est évacuée par convection naturelle entre le câble et l'air ambiant.

 

Cahier des charges

Les données utiles au traitement du problème, fournies par Epsilon, sont les suivantes :

 

 

1. Analyse préliminaire

 

Présentation générale d'un cas de convection naturelle autour d'un cylindre

En convection naturelle, lorsqu'un cylindre horizontal, de température $T_P=75 °C$ est immergé dans un fluide de température $T_{\infty}=20 °C$, une couche limite thermique se développe autour de celui-ci. D'après la figure ci-dessous[1], on peut constater que l'épaisseur de la couche limite thermique autour du cylindre est une fonction de l'angle $\theta$.

Panache de convection et couche limite thermique autour d'un cylindre chauffé

 

Analyse du cas-test

 

Hypothèses et mécanismes prédominants

  • écoulement stationnaire
  • écoulement établi
  • écoulement monophasique, pas de phénomène d'évaporation ou de condensation lors du transfert
  • problème à deux dimensions x et y
  • invariance par translation suivant l'axe z
  • symétrie selon le plan (y0z)
  • convection prépondérante devant conduction et rayonnement
  • convection naturelle
  • hypothèse de Boussinesq  $\rho=\rho_0 [1-\beta(T-T_0)]$   avec   $\beta= - \frac{1} {\rho} \left (\frac{\partial \rho} {\partial T} \right)_{P=P_0} $ le coefficient de dilatation thermique

 

Échelles caractéristiques

  • le diamètre du cylindre est $L = 1  cm$

 

Détermination des nombres adimensionnels

  • nombre de Prandtl

$$Pr=\frac {\nu_{air}} {\alpha_{air}}=0,698$$

 

  • nombre de Rayleigh basé sur le diamètre D du cylindre

$$Ra_D=\frac{g \beta (T_P-T_{\infty})D^3} { \nu \alpha}=5,26.10^3$$

avec, en considérant l'air comme un gaz parfait :

$\beta=\frac {1} { \overline T}=\frac {2} {T_P + T_{\infty}} = 3,12.10^{-3} K^{-1}$

 

  • nombre de Nusselt : déterminé par la corrélation théorique

 

Corrélation théorique

Dans le cas de la convection naturelle autour d'un cylindre, la corrélation à vérifier lors de nos simulations numériques est la suivante (Churchill et Chu - 1975) :

                              $\overline {Nu_D}=  \left \{0,60 + \Large{\frac {0,387 Ra_D^{1/6}} { \left [1+ \left (\frac{0,559} {Pr} \right)^{9/16} \right ]^{8/27}}} \right \}^2$       pour       $Ra_D<10^{12}$

Dans nos analyses, nous tracerons donc le nombre de Nusselt en fonction du nombre de Rayleigh afin de comparer les résultats numériques d'OpenFOAM et StarCCM+ avec cette corrélation.

 

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.

 

3. Simulations StarCCM+

Sous le logiciel commercial StarCCM+, les simulations sont beaucoup plus évidentes à mettre en place que sous OpenFOAM. Il faut par contre avoir les bonnes conditions aux limites et les bons modèles.

 

Géométrie du domaine et conditions aux limites

Afin de pouvoir observer le panache de convection, nous avons créé une géométrie rectangulaire beaucoup plus haute que large, en tenant compte des dimensions voulues par Epsilon pour le cylindre, celui-ci étant représenté par une condition limite de type wall à température fixée.

Schéma de la géométrie considérée et des conditions aux limites introduites dans StarCCM+

 

Maillage de la géométrie

Compte tenu du développement de la couche limite thermique sur le cylindre, il semble judicieux de raffiner le maillage d'autant plus que l'on se rapproche de la paroi du cylindre. Nous avons réalisé le maillage sous Ansys IcemCFD. Les détails du nombre de mailles et des lois de maillage sont présentées ci-dessous.

Maillage du domaine avec raffinement au niveau de la paroi du cylindre

 

Zoom sur la partie proche de la paroi du cylindre

 

Mise en place des simulations

Afin d'avoir des résultats les plus proches possible de ceux attendus, nous avons sélectionné plusieurs modèles dans StarCCM+ comme le montre la capture d'écran suivante. Nous avons entre autres sélectionné l'option segregated flow par opposition à coupled flow. La résolution de la vitesse avec segregated flow s'opère indépendamment de la pression, ce qui produit des simulations plus rapides, mais n'est pas adapté dans le cas d'écoulements supersoniques. Étant donné que nous avons des conditions de pressure outlet en entrée et en sortie du domaine, puisqu'il n'y a pas de vitesse en entrée, nous avons fixé la pression à 0 Pa. Nous avons imposé au cylindre une température de 348 K (75 °C).

Capture des modèles utilisés dans StarCCM+ pour la convection naturelle autour du cylindre

 

Résultats et vérification de la corrélation théorique

Nous avons ensuite réalisé les simulations pour différents modèles de turbulence qui seront présentés plus loin afin de voir l'influence qu'ils peuvent avoir sur la couche limite thermique et/ou sur la panache de convection.

 

Résultats en régime laminaire

Dans un premier temps, les simulations ont été effectuées en régime laminaire.

Résultat de simulation pour $T_P=348  K$  et  $T_\infty=293 K$ avec un modèle laminaire

 

L'objectif est de vérifier si les simulations apportent le résultat attendu. La figure ci-dessous montre l'évolution du nombre de Nusselt en fonction de l'angle $\theta$ autour du cylindre. Nous pouvons remarquer un maximum du nombre de Nusselt sur la partie inférieure du cylindre, c'est-à-dire pour un angle $\theta=-\frac {\pi} {2}$ car c'est sur cette partie que la couche limite est la plus fine, il y a donc des gradients thermiques beaucoup plus importants que sur la partie supérieure du cylindre, le nombre de Nusselt y est donc d'autant plus élevé.

Évolution du nombre de Nusselt en fonction de l'angle autour du cylindre en régime laminaire

 

Vérification de la corrélation

Maintenant que l'allure du nombre de Nusselt précédent est cohérent par rapport aux considérations physiques, nous allons dans un second temps vérifier la corrélation théorique de Churchill et Chu. Nous avons tracé les nombres de Nusselt théoriques puis ceux obtenus numériquement.

Rappelons la corrélation à vérifier ainsi que l'expression du nombre de Rayleigh :

                                            $\overline {Nu_D}=  \left \{0,60 + \Large{\frac {0,387 Ra_D^{1/6}} { \left [1+ \left (\frac{0,559} {Pr} \right)^{9/16} \right ]^{8/27}}} \right \}^2$       pour       $Ra_D<10^{12}$

 

Pour faire varier le nombre de Rayleigh, nous avons choisi de faire varier la gravité ainsi que la différence de température entre le cylindre et l'air environnant.

  • Pour la première figure ci-dessous, nous avons fait varier le nombre de Rayleigh pour les gravités suivantes : $g=9,81  m/s^2$, $g=98,1  m/s^2$, $g=981  m/s^2$ et $g=9810  m/s^2$.

Nombres de Nusselt numérique et théorique en fonction du nombre de Rayleigh calculé pour différentes gravités

 

  • Pour la seconde figure ci-dessous, nous avons fait varier le nombre de Rayleigh pour les différences de température : $\Delta T=55 °C$, $\Delta T=65 °C$, $\Delta T=75 °C$, $\Delta T=85 °C$, $\Delta T=95 °C$, $\Delta T=200 °C$  et  $\Delta T=1000 °C$.

Nombres de Nusselt numérique et théorique en fonction du nombre de Rayleigh calculé pour différents $\Delta T$

 

Les courbes numériques ci-dessous font apparaître des résultats relativement proches de la corrélation théorique avant de diverger de manière assez considérable à partir de nombres de Rayleigh relativement faibles. Dans le cas où nous avons fait varier le nombre de Rayleigh avec des gravités différentes (première figure), l'origine de cette divergence pourrait s'expliquer par l'aspect peu physique d'une gravité très importante, puisque nous somme allés jusqu'à une accélération de la pesanteur de $9810  m/s^2$.

 

Comparaison des modèles de turbulence

Nous avons également procédé à des simulations pour divers modèles de turbulence afin de voir leur influence sur la couche limite thermique ainsi que sur le panache de convection.

 

Résultats de simulation pour $T_P=348  K$  et  $T_\infty=293 K$ avec différents modèles de turbulence

Les simulations ci-dessous font apparaître des différences bien plus marquées au niveau du panache de convection qu'au niveau de la couche limite thermique. En effet, c'est dans ce panache que les effets de la turbulence peuvent se faire davantage ressentir du fait des instabilités qui peuvent se produire ainsi que les recirculations.

 

Le tracé du nombre de Nusselt ci-dessous permet aussi d'apprécier les différences qui apparaissent selon le modèle de turbulence utilisé. Ceux-ci sont quasi-confondus sur la partie supérieure du cylindre en $\theta=\frac {\pi} {2}$ tandis qu'ils diffèrent d'une unité pour les modèles Reynolds Stress (RSM) et realizable k-$\epsilon$ sur la partie inférieure de celui-ci. On peut donc en conclure que, quel que soit le modèle de turbulence utilisé, le panache qui s'élève du cylindre n'a aucune influence sur l'évolution du nombre de Nusselt.

Évolution du nombre de Nusselt numérique en fonction de l'angle autour du cylindre pour différents modèles de turbulence