Validation du solveur open-source OpenFOAM pour la résolution de problèmes thermiques

Validation du solveur OpenFOAM pour la résolution de problèmes thermiques

 

 

 

 

Equipe :
L'équipe est composée de Antoine Piedfert et de Yiming Gu, élèves en 3ème année à l'ENSEEIHT, filière Hydraulique et Mécanique des Fluides, et spécialisés en Energétique et Procédés.

Encadrants :
Le projet est encadré par Mme Annaig Pedrono et M. Dominique Legendre, enseignants à l'INP-ENSEEIHT et chercheurs à l'IMFT.

Contacts Industriels :
Ce projet est réalisé grâce aux informations fournies par MM. Doucet et Charalabidis, ingénieurs au sein du groupe Altran.

Introduction

Validation du solveur OpenFOAM pour la résolution de problèmes thermiques

 

 

 

 

Equipe :
L'équipe est composée de Antoine Piedfert et de Yiming Gu, élèves en 3ème année à l'ENSEEIHT, filière Hydraulique et Mécanique des Fluides, et spécialisés en Energétique et Procédés.

Encadrants :
Le projet est encadré par Mme Annaig Pedrono et M. Dominique Legendre, enseignants à l'INP-ENSEEIHT et chercheurs à l'IMFT.

Contacts Industriels :
Ce projet est réalisé grâce aux informations fournies par MM. Doucet et Charalabidis, ingénieurs au sein du groupe Altran.

Présentation du projet

Contexte industriel

Contexte Industriel

 

La simulation d'écoulements trouve énormément d'applications industrielles. Cependant, de nombreux marchés restent inexploités en raison du coût des logiciels de CFD. Pour cette raison, le groupe Altran a commencé un projet de validation du solveur libre OpenFOAM, ce qui permettrait à des entreprises plus petites et avec moins de moyens d'avoir accès à un logiciel gratuit.

La validation thermique d'OpenFOAM passe par la validation de nombreux cas tests, effectués lors d'expériences ou de simulations en laboratoires validés par l'Université de Manchester, et disponibles sur le site Ercoftac. Il conviendra de valider des solveurs incompressibles, laminaires et turbulents et de modéliser des écoulements confinés, en 2D et en 3D.

Nous avons donc choisi un cas, le numéro 45 de la base de donnée Ercoftac, simulant un écoulement turbulent entre deux plaques parallèles chauffées.

Présentation d'OpenFOAM

Présentation d'OpenFOAM

OpenFOAM est un logiciel de CFD libre développé par OpenCFD Ltd au Groupe ESI et distribué par OpenFOAM Foundation. Il utilise la méthode des volumes finis pour résoudre les équations aux dérivées partielles connues dans le domaine de la mécanique des fluides.

Un des principaux avantages d'OpenFOAM est qu'il propose  de très nombreux solveurs codés en C, modifiable par l'utilisateur. Celui ci peut alors calculer de nouvelles variables en fonction de nouvelles équations différentielles. La difficulté de ce logiciel est qu'il ne possède pas d'interface graphique : toutes les informations sont contenues dans des fichiers texte, ce qui peut le rendre difficile à prendre en main.

La génération du maillage peut soit se faire à partir du logiciel libre Salome, soit directement dans les fichiers texte du cas en utilisant une simple ligne de code "blockMesh". Un autre outil permet de faire le post-processing : paraFoam, qui est une variante de Paraview adapté à la lecture de résultats obtenus avec OpenFOAM.

Après une prise en main et une recherche des fonctions à utiliser pour les conditions aux limites, voici le premier résultat obtenu à partir d'icofoam, solveur d'écoulements incompressibles et laminaires :

Bilan de la phase préliminaire

Bilan de la phase préliminaire et de la première semaine

Premier contact
Lors de la phase préliminaire, nous avons été en contact avec les ingénieurs d'Altran, qui nous ont expliqué le but de ce projet et son contexte. Nous avons pu échanger avec eux pour savoir ce qu'ils attendraient de nous pendant ces 6 semaines.

Prise en main d'OpenFOAM
Nous nous sommes directement lancés dans les tutoriels d'OpenFOAM afin d'en comprendre le fonctionnement, et ce jusqu'à la fin de la première semaine de la phase intensive. La première étude consistait en une cavité 2D, entourée de murs, dont un était un mur mobile. Nous avons pu d'abord observer le champ de vitesse crée. Le solveur utilisé est icoFOAM, solveur d'écoulements laminaires incompressibles.

Ensuite, nous avons pu apprendre à modifier le solveur afin qu'il calcule aussi le champ de température crée si le mur mobile est en plus chauffé. à flux constant.

Adapter ces connaissances aux futurs solveurs devrait être faisable sans trop de difficultés maintenant.

Etude du cas n°45

Le cas de l'écoulement entre deux plaques parallèles chauffées a déjà fait l'objet de plusieurs études et simulations dont certains résultats sont en libre accès sur internet. Nous avons pu commencer à étudier la théorie à l'aide de ces articles ainsi que des corrélations sur les échanges thermiques connues.

Détail du cas test sélectionné

Cas test sélectionné

 

Le cas que nous avons selectionné est un écoulement entre deux plaques parallèles chauffées.

 

Hypothèses

La température est imposée aux parois. On negligera tout effet de compressibilité dans l'étude, ce qui fait de la température un scalaire passif.

On est en 2D : $\frac{d}{dz}=0$ et $U_z=0$
Les effets de compressibilité sont négligés : $\rho=constante$
On considère l'écoulement stationnaire : $\frac{\partial}{\partial t}=0$
On néglige les effets de la gravité.

Equations gouvernantes

Conservation de la masse 
$ \vec \nabla .\vec U = 0 $

Conservation de la quantité de mouvement
$ (\vec U. \vec \nabla) \vec U = -\frac{1}{\rho} \vec{grad} P +\nu \Delta^2 \vec U$

Conservation de l'énergie
$ (\vec U. \vec \nabla) T=\alpha \Delta^2 T$

 

Les résultats expérimentaux seront les résultats obtenus par une DNS du cas par les professeurs Kasagi et Kuroda, du département d'ingénierie mécanique de l'Université de Tokyo. Ils fournissent dans leurs données les valeurs des fluctuations de vitesse et de température, le nombre de Prandtl turbulent... etc

 

Résultats

Détails des simulations

1. Les solveurs

Nous avions beaucoup de solveurs capables de simuler la turbulence. Voici les solveurs que nous avons choisis d'étudier :

pisoFoam : Méthodes LES et RAS (k-epsilon, Spalart Allmaras...)
dnsFoam : Méthode de Direct Numerical Simulation
channelFoam : Un solver propre au cas de l'écoulement entre plaques parallèles.

Tous ces solveurs résolvent les équations de Navier-Stokes incompressibles, plus d'autres équations modélisant la turbulence (sauf dnsFoam, qui calcule tous les effets sans modèle de turbulence).

Les trois solveurs ont dû être modifiés afin de calculer la température comme un scalaire passif, résolvant l'équation de conduction de la chaleur décrite plus loin dans ce site.

2. Les constantes

L'écoulement se caractérise par les nombre adimensionels suivants :
$Re=\frac {2 \delta U_m}{\nu}=4560$
$Re_\tau =\frac {\delta U_\tau}{\nu}=150$

$Pr=\frac{\alpha}{\nu}=0.71$

Avec :
$\alpha$ le coefficient de diffusion de température de l'air en $m^2 s^{-1}$
$\nu$ la viscosité cinématique, $\nu=1.0*10^{-5} m^2 s^{-1}$ pour de l'air.
$\delta$ la demi-largeur de l'écoulement
 

Tous les résultats présentés seront adimensionalisés par les constantes suivantes :

$ u_\tau = \sqrt{\nu \frac{dU}{dy}} $ la vitesse de friction,
$ T_\tau = \frac{q_w}{C_p u_\tau \rho} $ la température de friction, où $q_w=\lambda \frac{dT}{dx} $
$ \nu $ la viscosité cinématique de l'air.

3. Le maillage

Le maillage représente un écoulement de dimension $60\delta.2\delta.2\delta$. Nous avons choisi de modéliser un écoulement de taille $60\delta$ suivant l'axe des x car c'est environ la distance nécessaire à l'établissement d'un écoulement turbulent instationnaire pour les vitesses ET pour la température.
Le maillage possède 400*60*1=24000 points. Il est raffiné au niveau de la paroi tel que le dy des cellules en paroi soit 10 fois plus petit que le dy au centre.

PisoFoam

PisoFoam

Le solveur PisoFoam permet de simuler plusieurs modèles de turbulence, RAS (k-epsilon, k-oméga, Spalart Allmaras...) ou LES. Nous avons choisi de comparer les modèles suivants :

  • Modèle k-epsilon
  • Modèle Spalart Allmaras
  • Modèle Large Eddy Simulation

Ces modèles figurent parmi les plus courants en CFD.

Les constantes en entrée des modèles sont obtenues par les calculs suivants :

$I = 0.16 Re^{-1/8}=5.6$% l'intensité turbulente.
$k=1.5 (U_m I)^2 = 2.5 .10^{-6} m^2 s^{-2} $
$ \epsilon = 0.09^{3/4} k^{3/2} I^{-1} = 8.9 .10^{-9}m^2 s^{-3}$

Après adimensionalisation, voici les résultats obtenus :

Pour la vitesse, c'est le modèle k-epsilon qui est le plus précis. Les deux autres modèles sont assez peu satisfaisants.

Encore une fois, le modèle k-epsilon est le plus satisfaisant. Cependant, on note que les écarts de température ne sont pas aussi importants que les écarts de vitesse, ce qui est un point positif pour la validation des problèmes thermiques. Cette information doit toutefois être prise avec des pincettes, car il est probable que les erreurs en température varient avec le nombre de Prandtl.

Influence du maillage

Influence du maillage

Le modèle k-$\epsilon$ étant parmi les plus rapides, nous l'avons utilisé pour observer l'influence de différents paramètres, à commencer par le maillage.

Le premier maillage 300x60x1 est un maillage 2D pour lequel les mailles proche parois sont 10 fois plus fines que celles au centre de l'écoulement.

Le second est identique, mais est moins raffiné dans la direction de l'écoulement.

Le troisième est un maillage identique au premier, mais non raffiné au parois. Toutes les mailles sont de dimensions égales.

Le maillage raffiné n'est pas adapté à notre cas. Il était peu probable que le résultat soit cohérent pour le modèle k-$\epsilon$, mais cela expliquera certains résultats à venir.

C'est bien sûr le premier qui semble le plus précis. Le second maillage est très précis aussi, bien que plus grossier. La visualisation de la température moyenne confirme la tendance :

L'erreur sur le maillage régulier n'est pas négligeable. Pour des soucis de précision, nous garderons le maillage 300x60x1, qui reste le plus précis même pour la température.

Influence du paramètre k

Influence du paramètre k

Nous avons testé le calcul en surestimant la valeur de k, c'est-à-dire en créant trop de turbulence par rapport à la réalité, pour voir s'il y a un effet sur la température. Nous avons calculé la valeur de $ \overline{u'^2} $ pour voir l'évolution de la turbulence :

Nous voyons que $ \overline{u'^2} $ est plus grand et passe par un maximum plus proche du mur lorsque k est surestimé, ce qui est logique. Voyons maintenant les effets sur la vitesse et sur la température :

La différence est relativement faible malgré le grand écart entre les valeurs de k, mais on observe quand même une meilleure approche lorsque l'on conserve la bonne valeur de k.

 

DnsFoam

Le solveur permet de calculer tous les effets sans avoir besoin de modèle pour la turbulence en utilisant la méthode Direct Numerical Simulation.

Ce solveur a été la source de beaucoup de problèmes.

Premièrement, dnsFoam est un solveur qui utilise des transformées de Fourier. Pour cette raison, il impose plusieurs conditions sur le maillage. Celui-ci doit être complètement régulier, la méthode de fft n'étant pas valable sur des maillages irréguliers. C'est-à-dire qu'il nous a fallu adapter notre maillage, que nous utilisions avec pisoFoam, pour qu'il soit régulier, mais aussi suffisamment fin.  Il a fallu du temps pour déterminer que c'était l'irrégularité du maillage qui empêchait le code de tourner. Le temps de calcul a de plus été considérablement augmenté par ce problème.

Deuxièmement, dnsFoam nous demande d'entrer certaines valeurs pour modéliser la turbulence. Bien que cela semble étonnant au premier abord, car le principe même de la DNS est qu'elle n'a besoin d'aucune valeur de modèle pour modéliser correctement la turbulence. Pourtant, lorsque l'on garde les valeurs proposées dans le tutorial, voilà ce que l'on obtient :

La turbulence n'est pas vraiment réaliste, très peu diffuse et de grande amplitude. Nous avons effectué beaucoup de recherche pour comprendre l'utilité de ces 4 valeurs (nommées UOkappa, UOsigma, UOKupper et UOKlower). Les forums dédiés à OpenFoam ne nous ont pas été d'une grande aide, le peu d'utilisateurs semblant les connaître restant très flou sur la signification de ces valeurs. Nous avons trouvé que ce sont les 4 constantes nécessaires à la création d'une loi de probabilité uniforme, servant à forcer la turbulence (valeurs minimale et maximale, écart-type...) dans le cadre d'un procédé dit de Ornstein–Uhlenbeck. Nous ne somme malheureusement pas en mesure d'expliquer à quoi correspondent ces valeurs aléatoires. De plus amples informations peuvent être trouvées sur la page du site Wikipedia. Pour les simulations suivantes, nous avons mis toutes ces valeurs à 0, en espérant que la turbulence serait correctement modélisée.

Malgré une erreur dans la vitesse, on observe que la température correspond bien au profil, en particulier loin de la paroi. Il faut aussi prendre en compte que la DNS (comme la LES) n'a pas vraiment de signification physique lorsqu'elle est en 2D, et que l'erreur peut venir de là. Cependant, les contraintes de maillages ne nous ont pas permis de tenter un calcul avec davantage de mailles en 3D.

ChannelFoam

ChannelFoam est un solveur sensé être conçu spécialement pour les cas d'écoulements entre plaques parallèles. Cependant, son fonctionnement n'est pas du tout intuitif. Notamment à cause  d'un double dossier "0", c'est-à-dire qu'il existe deux dossiers dans lesquels il faut rentrer les conditions initiales. Nous n'avons pas été en mesure de déterminer l'interêt de chacun des dossiers.

 

Voici les résultats obtenus avec ChannelFoam, qui a fonctionné en 3D :

Les résultats sont assez satisfaisants, bien qu'on observe comme sur les autres solveurs quelques erreurs dans une zone peu éloignée du mur.