III. Validation du logiciel OpenFOAM

Avant d'aborder les simulations numériques avec OpenFOAM sur les cas-tests proposés par Epsilon, il faut au préalable avoir au moins un point de validation du logiciel. Notre validation portera sur l'étude de la convection forcée entre une plaque plane chauffée et un écoulement d'air à température inférieure.

1. Analyse préliminaire

Cette validation concerne l'étude du transfert thermique entre une plaque plane chauffée et un courant d'air incident. Il s'agit par conséquent d'un cas de transfert thermique par convection forcée.

 

Présentation générale d'un écoulement sur plaque plane chauffée

Un écoulement incident sur une plaque plane chauffée sera principalement caractérisé par le développement d'une couche limite dynamique (en vitesse) et d'une couche limite thermique (en température).

  • L'origine du développement de la couche limite dynamique provient des effets de viscosité, c'est-à-dire de diffusion de quantité de mouvement, entre le fluide incident et la plaque. En effet, la vitesse du fluide est nulle sur la plaque, tandis qu'elle est non nulle à l'infini, ce qui implique qu'au voisinage de la plaque il y ait une zone de transition : la couche la plus lente freine la couche la plus rapide, qui, en retour, l'accélère. La couche limite dynamique est une zone où la variation de vitesse est la plus marquée, c'est-à-dire là où l'on trouve de forts gradients de vitesse. La vitesse à l'extrémité de la couche limite est égale à 99% de la vitesse à l'infini, on a donc :

$$U(\delta)=0,99  U_0$$

L'épaisseur de la couche limite varie selon la relation suivante[3] :

$$\delta(x)= 4,64 \sqrt \frac {\nu  x} {U_0}$$    

$\nu$ : viscosité cinématique du fluide incident en $m^2/s$

$x$ : abscisse de la plaque en $m$

$U_0$ : vitesse du fluide incident en $m/s$

 

  • De la même manière que la couche limite dynamique, la couche limite thermique est une zone proche de la plaque où la température varie de façon significative entre la température de la plaque et la température du fluide incident. La température à l'extrémité de la couche limite est telle que :

$$\frac{T(\delta_T)-T_P} {T_\infty - T_P} = 0,99$$

 

L'épaisseur de la couche limite thermique varie selon l'expression suivante[3] :

$$\delta_T= \frac { \delta} {1,026  Pr^{1/3}} =5,09 \sqrt \frac {\nu  x} {U_0}$$  

$\alpha$ : diffusivité thermique du fluide incident en $m^2/s$

Régimes d'écoulement sur plaque plane

 

Position du cas de validation

Notre cas de validation se présente tel que l'illustre la figure ci-dessous. De l'air incident à une vitesse de 50 cm/s et à une température de 293 K arrive sur une plaque plane chauffée à une température de 348 K.

Schéma du cas de validation du logiciel OpenFOAM

 

Avant de procéder à toute simulation, une analyse préliminaire est requise afin d'identifier les phénomènes et ordres de grandeur.

 

Hypothèses et mécanismes prédominants

  • écoulement stationnaire
  • écoulement établi
  • écoulement monophasique, pas de phénomène d'évaporation ou de condensation
  • problème à deux dimensions x et y
  • invariance par translation suivant l'axe z
  • convection prépondérante devant conduction et rayonnement
  • convection forcée
  • température considérée comme un scalaire passif (aucune influence sur l'hydrodynamique)

 

Échelles caractéristiques

  • longueur de la plaque : $L = 10  cm$
  • vitesse incidente : $U_0=50  cm/s$
  • iso-contour de température à l'épaisseur de la couche limite thermique : $T(\delta_T)=0,99  (T_{\infty}-T_P)+T_P=293,55  K$
  • iso-contour de vitesse à l'épaisseur de la couche limite dynamique :  $U(\delta)=0,495   m/s$

 

Détermination des nombres adimensionnels

  • nombre de Reynolds

Pour une plaque plane, le nombre de Reynolds de transition entre les régimes laminaire et turbulent est environ :

$${Re_{tr} \approx 5.10^5} $$

Dans notre cas de validation, le nombre de Reynolds à l'extrémité de la plaque s'écrit :

$${Re}=\frac{U_0 L} {\nu_{air}}=\frac {50.10^{-2} \times 10.10^{-2}} {1,495.10^{-5}}=3344,5$$

L'écoulement au-dessus de la plaque est donc laminaire en tout point de celle-ci.

  • nombre de Prandtl

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

 

Corrélation théorique

Dans notre cas de convection forcée sur plaque plane, d'après l'ouvrage Fundamentals of Heat and Mass Transfers, Sixth Edition de Frank P. Incropera et David P. DeWitt, la corrélation qui sera utilisée pour vérifier nos résultats est la suivante, qui relie le nombre de Nusselt local aux nombres de Prandtl et de Reynolds locaux :

 

$$Nu_x = \frac{h_x x}{\lambda} = f(Pr) {Re_x}^{0,5}$$

avec :                                            $f(Pr)=0,332  Pr^{1/3}$    pour    $0,6<Pr<50$

 

2. Création du maillage

Le maillage sous OpenFOAM a été réalisé avec l'outil blockMesh. Cet outil permet de mailler des géométries relativement simples puisqu'il s'agit d'un mailleur minimaliste. Dans le cas d'une géométrie beaucoup plus compliquée qu'une plaque plane, on lui préfèrera des mailleurs externes tels que Salomé, I-deas, Gmsh... Cependant, OpenFOAM dispose d'autres outils de maillage.

 

Caractéristiques du maillage

Avec blockMesh, le maillage doit être entièrement réalisé sans interface graphique. Ainsi, il faut définir point par point les différents éléments constitutifs de la géométrie dans un fichier texte, puis spécifier les paramètres de maillage.

Pour notre cas de validation, nous devons définir 2 zones : celle dont la surface inférieure sera la plaque, et celle en amont de la plaque. Nous raffinons ensuite au niveau de plaque pour une meilleure résolution des couches limites dynamique et thermique en utilisant un maillage géométrique avec un rapport 100 entre la maille la plus grande (le plus loin de la paroi) et la maille la plus petite (proche de la paroi).

Schéma explicatif de la géométrie, division du domaine en deux blocs

 

Maillage du domaine obtenu avec blockMesh après visualisation avec paraFoam

 

Les caractéristiques du maillage sont éditées dans le fichier blockMeshDict.

 

Capture du fichier blockMeshDict correspondant à notre maillage de la plaque plane

3. Simulations et résultats

Mise en place des simulations

Dans ce cas de validation, nous allons faire l'hypothèse que la température n'a aucune influence sur l'hydrodynamique, c'est à dire que la température est considérée comme scalaire passif. L'écoulement est incompressible et les vitesses considérées en entrée sont telles que la transition vers un état turbulent ne sera jamais observée le long de la plaque.

Dans ces conditions, nous avons choisi d'utiliser le solveur incompressible icoFoam et d'y intégrer l'équation de la chaleur (avec terme de convection), le nouveau solveur porte ainsi le nom de my_icoFoam.

$$ \frac {\partial T} {\partial t} + Div(\phi T) = \frac {\lambda} {\rho c_p} \Delta T $$

Cette modification du solveur a déjà été étudiée : un tutoriel en anglais est disponible à cette adresse. L'idée pour nous était de voir comment il est possible d'intégrer simplement de nouvelles variables et équations à un modèle déjà existant. Ceci nous permet d'évaluer le potentiel d'OpenFOAM en terme de flexibilité, un des atouts majeurs de tout logiciel "open source".

De plus, avant toute simulation, il faut veiller à avoir un nombre de Courant inférieur à 1 :

$$Co=\frac {U  \Delta t} {\Delta x} < 1$$

$U$ : vitesse dans la direction $x$

$\Delta t$ : intervalle de temps à choisir pour avoir un nombre de Courant au plus égal à 1

$\Delta x$ : distance entre deux mailles

Dans notre cas, le pas de temps à choisir est au plus :

$$\Delta t=\frac {Co  \Delta x} {U} =\frac{1 \times  1.10^{-3}} {50.10^{-2}}=2  ms$$

 

Visualisation des résultats

La procédure pour effectuer la simulation et visualiser les résultats se fait par le biais du terminal en trois parties principales :

0. Se placer dans le répertoire de travail contenant le cas de calcul

         

Contenu du cas de calcul "Validation" avant simulation

 

1. Mailler le domaine en appelant l'outil blockMesh

 

2. Résoudre les équations mises en jeu en appelant le solveur approprié, ici my_icoFoam

3. Afficher les résultats dans le logiciel Paraview en appelant l'outil de post-traitement paraFoam

Champ de température sur la plaque avec l'iso-contour $T(\delta_T)=293,55  K$ avec une vitesse d'entrée $U_0=50  cm/s$ en régime établi (t > 0,4 s)

 

  

Champ de vitesse sur la plaque avec l'iso-contour $U(\delta)=0,495  m/s$ avec une vitesse d'entrée $U_0=50  cm/s$ en régime établi (t > 0,4 s)

 

 

Visualisation des couches limites

Afin de comparer les épaisseurs des couches limites thermique et dynamique, on peut les visualiser sous Paraview en utilisant la fonction isosurface en spécifiant les valeurs correspondantes en vitesse et température.

En exportant les résultats sur Matlab, il est alors possible de superposer les profils numérique et théorique :

 

Tracé des couches limites thermique et dynamique sur la plaque plane d'après les iso-contours OpenFOAM

 

On remarque d'une part que la couche limite thermique se développe plus vite que la couche limite dynamique. Ceci s'explique par le fait que le nombre de Prandtl est inférieur à 1.

D'autre part, on constate un certain écart avec la théorie : la couche limite thermique est sur-estimée tandis que la couche limite dynamique est quant à elle sous-estimée. Deux raisons peuvent expliquer cet écart :

  • Vitesse et température loin de la plaque sont mal évaluées (par exemple, la vitesse loin de la plaque n'est pas exactement égale à la vitesse imposée en entrée car la plaque modifie l'écoulement)
  • Les paramètres $\nu$ et $\alpha$ ne sont pas estimés à la bonne température

 

 

Evolution du Nusselt local

Afin de calculer le Nusselt local le long de la plaque, il nous faut connaître l'évolution du gradient de température sur la plaque, perpendiculairement à celle-ci.

Ceci a été réalisé en ajoutant un nouveau champ dans le fichier icoFoam_gradT.C et calculé seulement le long de la plaque. Le code à rajouter doit s'effectuer sur le fichier précédent placé dans le dossier nommé icoFoam_gradT localisé dans les dossiers applications puis solvers dont le contenu est le suivant :

Les lignes de code à ajouter sont les suivantes :

Capture du fichier texte "icoFoam_gradT.C"

 

La nouvelle variable ainsi créée apparaît désormais dans Paraview, nous permettant de connaître l'évolution de la dérivée normale à la paroi de la température. En adimensionnalisant cette grandeur (on la multiplie par $ \frac{x}{\Delta T} $), on obtient le profil du Nusselt local.

La figure suivante montre les résultats obtenus pour différentes vitesses en entrée :

vitesses nombre de Reynolds
$U_0=5  mm/s$ $33,44$
$U_0=5  cm/s$ $334,4$
$U_0=50  cm/s$ $3344$
$U_0=5  m/s$ $33440$

 

Le tracé du Nusselt montre d'excellents résultats vis-à-vis de la théorie, même dans le cas d'une vitesse d'entrée de 5 m/s qui donne un nombre de Reynolds s'approchant de la valeur critique de transition laminaire/turbulente. On peut cependant observer quelques différences au niveau du bord d'attaque de la paroi : cela s'explique par la résolution du maillage qui n'est pas assez élevée dans cette zone où les gradients sont les plus forts.