Calcul numérique ERCOFTAC : Mise en place de la simulation numérique

Nos premières simulations on été effectuées pour représenter uniquement la géométrie introduite par le site ERCOFTAC (cas 09).

Ainsi, notre domaine de calcul se réduira à la zone où la couche limite turbulente est pleinement développée.

 

 

 

 

 

 

Analyse théorique et étude préliminaire

Dans cette partie, nous présentons les équations analytiques qui permettent de décrire ce problème de couche limite turbulente.

Pour simplifier le problème, nous ferons les hypothèses les hypothèses suivantes :

  • On considère que le problème est 2D plan $ \Rightarrow$ $  \frac{\partial} {\partial z} =0 $ et Uz=0.
  • On considère que le problème est stationnaire $ \Rightarrow$ $ \frac{\partial} {\partial t} =0 $.
  • On suppose que l'air est un gaz parfait.
  • On considère que l'écoulement suit l'hypothèse de Boussinesq, i.e les propriétés du fluide (l'air) sont constantes, sauf la masse volumique $ \rho $ dans le terme des forces de volume (la gravité) qui ne dépend que de la température T.

          $ \frac{\rho}{\rho_r} = 1 - \beta (T-T_r) $

Avec :   $T_r$  : La température de référence, pour notre cas $T_r = 15°C $

              $\rho_r$ : La masse volumique de référence.

              $\beta $ : Le coefficient de dilatation thermique (en K-1)
 
Par définition, $\beta = - \frac{1}{\rho} \frac{d \rho}{d T}(T=T_r)$ , pour un gaz parfait (notre cas) $\beta =\frac{1}{T_r} $ .
A.N     $\beta = 3.47 1e-3 $
  • On suppose que $ \frac{\partial} {\partial y} >>  \frac{\partial} {\partial x} $ .

Avec ces hypothèses, les équations de conservation en couche limite en convection naturelle sont :

L'équation de continuité :

               $ \frac{\partial U} {\partial x} +  \frac{\partial V} {\partial y} =0 $

L'équation de conservation de quantité de mouvement :

$ U \frac{\partial U} {\partial x} + V \frac{\partial U} {\partial y}  =  \frac{\partial }{\partial y} (\nu \frac{\partial U}{\partial y})+ g \beta (T-T_r) $

Quant à l'équation de l'énergie (ou de la température), puisqu'en convection naturelle, l'écoulement est à faible nombre de Mach (vitesses faibles associées à la convection naturelle), donc pas de terme de dissipation, par suite cette équation donne :

L'équation de l'énergie :

$ U \frac{\partial T} {\partial x} + V \frac{\partial T} {\partial y}  = \frac{\partial  }{\partial y} (\alpha ​\frac{\partial T}{\partial y}) $

On remarque bien que le terme de flottabilité couple le problème hydrodynamique et le problème thermique. Certes le modèle de Boussinesq complique le problème par rapport à celui du problème incompressible, mais c'est beaucoup plus simple que celui compressible.

Nous pouvons à ce stade, mener une étude dimensionnelle pour quantifier les paramètres adimensionnels qui contrôlent notre problème. En adimensionnalisant les équations précédentes, on voit bien l'apparition du fameux nombre de Grashof   : $ Gr_L=\frac{\beta g L^3 (T_w-T_a)}{\alpha^2} $. Ce nombre qui compare l'effet des forces de flottabilité à celui des forces visqueuses, joue le même rôle que celui du nombre de Reynolds en convection forcée. Et donc pour comparer l'effet de la convection naturelle à celui de la convection forcée, il suffit de comparer le nombre de Grashof et la racine du Reynolds.

 

Modélisation de la turbulence :

Pour modéliser la turbulence, nous allons nous baser dans un premier temps sur la résolution RANS. Cette approche consiste à moyenner les équations précédentes (moyenne d'ensemble) en utilisant la décomposition de Reynolds.

​Avec les mêmes hypothèses précédentes, on obtient :

L'équation de continuité :

$ \frac{\partial \mathbf{U}} {\partial x} +  \frac{\partial \mathbf{V}} {\partial y} =0 $

L'équation de conservation de quantité de mouvement :

$ \mathbf{U} \frac{\partial \mathbf{U}​ } {\partial x} + \mathbf{V} \frac{\partial \mathbf{U}} {\partial y}  =  \frac{\partial}{\partial y} (\nu \frac{\partial \mathbf{U}}{\partial y} - \overline{u' v'})+ g \beta (\mathbf{T}-T_r) $​

L'équation de l'énergie :

$ \mathbf{U} \frac{\partial \mathbf{T}} {\partial x} + \mathbf{V} \frac{\partial \mathbf{T}} {\partial y}  = \frac{\partial  }{\partial y} (\alpha ​\frac{\partial \mathbf{T}}{\partial y} - \overline{v' T'}) $

Les grandeurs en gras représentent des grandeurs moyennes, les termes en ' représentent les fluctuations.

Le long de ce projet, nous allons adopter l'approche RANS. Cependant, nous allons tester plusieurs modèles de turbulences : k-epsilon, k-epsilon RNG, RSM,...etc

Création du maillage

Maillage sous SALOME

La plate-forme SALOME est principalement développée par le CEA et EDF, celle ci permet de créer et mailler des géométries, ainsi que de modifier les propriétés numériques et physiques d'éléments géométriques et enfin de visualiser et post-traiter des champs de résultats. Cet outil est open-source, c'est pour cette raison que nous réalisons notre maillage avec, puisque l'integralité de la chaîne de calcul doit être libre.

On réalise ainsi un premier maillage régulier et uniforme (1cm la maille) de la géométrie présentée précédemment. 

Ensuite nous effectuerons un raffinement, avec un ratio 10, au niveau de la zone proche paroi et ce afin de bien représenter la couche limite qui se développe autour de la plaque chauffée.

Maillage sous OpenFOAM

Ce même maillage peut se faire directement sous OpenFoam puisque la géométrie est assez simple à construire, ceci en modifiant le fichier BlockMeshDict du cas test Cavity, mais cette fois-ci en raffinant aux parois avec un ratio 10.

   

 

Pour nos simulations, nous allons utiliser plusieurs maillages avec des ratios différents, afin d'analyser d'une part améliorer les résultats, et d'autre part, analyser son effet sur la qualité de ces résultats.

Paramètres de simulation

Conditions initiales

L'entrée

La difficulté principale de ce calcul est liée au choix de configuration, où l'entrée n'est plus libre, puisque elle est positionnée à 1.4 mètres en dessus du bord d'attaque. Il faut donc imposer en entrée, des profils non uniformes en vitesse ainsi qu'en température, issus des mesures expérimentales que nous représentons ci-dessous.

Il est à noter que ces données représentent un intervalle allant de 0 jusqu'à 12 cm à partir de la plaque verticale, sans pour autant donner des plus informations sur le reste du domaine d'entrée. Il va donc falloir choisir un profil qui représente plus ou moins ce qui se passe d'un point de vu dynamique et thermique.

Ainsi, nous imposons des profils linéaires qui relient les dernières valeurs données expérimentalement en y=0.12 m jusqu'à T=288.15K et U=1e-4m/s, en y=0.5 m.

NB: Le profil imposé en entrée doit coïncider avec les positions des points du maillage. Pour ce, nous effectuons des opérations d'interpolation via un programme Matlab qui permet de réaliser cette opération quelque soit la maillage mis en plaque.

 

Conditions limites

Le choix des conditions aux limites aux niveaux de la sortie et du bord de gauche, n'a pas été évident. Plusieurs essais ont été mis en place avant d'avoir choisit les bonnes BC qu'on résume dans le tableau ci-dessous.

  Entrée Plaque Sortie Bord
U (m/s) nonuniform List<vector>

fixedValue uniform (0 0 0)

zeroGradient zeroGradient
p (Pa) zeroGradient zeroGradient

fixedValue uniform

fixedValue uniform
T (K) nonuniform List<scalar>

fixedValue uniform 333.15K

zeroGradient zeroGradient

 

Mise en place des calculs

La mise en place de ce calcul, nécessite d'abord la création d'un répertoire contenant les fichiers où seront décrits les différents paramètres physiques et numériques, indispensables pour la simulation.

Pour ce, nous nous place dans le répertoire /heatTransfer/buoyantBoussinesqSimpleFoam où nous pouvons copier le tutoriel hotRoom sous un nouveau nom  > cp  -r  hotRoom  plaque, et comme pour tout calcul sous OpenFOAM, nous modifions les fichiers correspondants aux paramètres physiques, contenus dans les répertoires 0 > gedit 0/* & , ainsi que les paramètres numériques contenus dans /systeme, sans oublier de vérifier le maillage, les propriétés thermophysiques, le sens de la gravité et le modèle de turbulence, dans le répertoire /constant.

 

Une première tentative de mise en place du calcul, avec les profils de vitesse et de température non-uniforms en entrée, n'aboutit finalement pas et l'on obtient le message d'erreur ci-dessous:

--> FOAM FATAL ERROR : Continuity error cannot be removed by adjusting the outflow.

Please check the velocity boundary conditions and/or run potentialFoam to initialise the outflow.

Cette erreur est due à une discontinuité entre le profil des vitesses, non nul, imposé en entrée avec la condition dans le le domaine intérieur: internalField U (0 0 0).

Pour résoudre ce problème de discontinuité, nous passons par une étape intermédiaire, qui revient à lancer un nouveau calcul via le solveur potentialFoam.

 

Calcul sous PotentialFoam

PotentielFoam est un solveur d'écoulement potentiel qui peut être utilisé pour initialiser le domaine et résoudre tout problème de discontinuité.

Il suffit donc de se placer dans basic/potentialFoam/, de copier le cas "cylinder" sous un nouveau nom et ensuite récupérer le fichier des vitesses U crée précédemment.

Les résultats de ce calcul doivent être copiés dans le répertoire 0 du calcul sous buoyantBoussinesqSimpleFoam.

 

Calcul sous buoyantBoussinesqSimpleFoam

  • Calcul en parallèle

Il est relativement facile de faire fonctionner un calcul en parallèle, sur plusieurs processeurs. Il convient tout d'abord de copier le fichier decomposeParDict dans le dossier system du cas de calcul

> cd system

>  cp OpenFOAM/user/run/tutorials/multiphase/interFoam/ras/damBreak/system/decomposeParDict

Le fichier decomposeParDict peut être détaillé de la façon suivante :

Pour obtenir davantage de renseignements, le site suivant permettra de fournir davantage d'explications (notamment sur la manière de faire du calcul en parallèle en utilisant plusieurs machines d'un même réseau) :

http://www.openfoam.org/docs/user/running-applications-parallel.php#x12-820003.4

Dans notre cas, nous utiliserons directement le fichier copié tel quel, puisque nous voulons travailler en utilisant pour une unique simulation les 4 processeurs d'une même machine.  Ainsi, il convient de faire :

> cd nom_du_cas

> decomposePar

A ce stade, des dossiers processor*/ ont du apparaître. On vient de partager les domaines voulus sur les quatre processeurs. Enfin, on lace le calcul avec la commande générale :

mpirun --hostfile <machines> -np <nProcs> <foamExec> <otherArgs> -parallel > log &

Soit, dans notre cas simple :

> mpirun -np 4 buoyantBoussinesqSimpleFoam -parallel > log &

 

On peut surveiller l'évolution du calcul en surveillant l'évolution du fichier "log", en entrant :

> tail -f log

Une fois que le calcul est terminé, il faut rassembler les résultats :

> reconstructPar

 

On obtient maintenant les divers répertoires pour les différents temps physiques voulus, et on peut exécuter le post-traitement habituel avec :

> paraFoam

 

  • Convergence de la simulation:

Afin de vérifier la convergence de la simulation, un calcul des résidus peut se faire via la commande:

> pyFoamPlotRunner.py buoyantBoussinesqSimpleFoam

 

Présentation des résultats

Après convergence du calcul et atteinte du régime permanent, nous visualisons nos résultats sur paraview.

 

Vitesse (m/s)


Température (°C)


Lignes de courant (vitesse)

 

La visualisation sur paraview permet donc de valider ne serait ce que qualitativement les résultats de nos calculs.

Difficultées rencontrées lors de l'implantation des conditions limites

Il est tout d'abord extrêmement difficile de bien comprendre la signification des conditions limites proposées sur OpenFOAM. Dans notre cas, seuls de long essais nous ont permis de déterminer une fois pour toutes les conditions limites à utiliser pour nos simulations. L'implémentation des conditions limites fut la partie la plus longue du projet.

Cependant, les simulations dans leur ensemble étaient particulièrement difficiles à mettre en place, notamment pour les profils à l'entrée. Les valeurs expérimentales données permettaient d'établir un profil de vitesse ou de température, pour certains points donnés seulement. Notre maillage ne s'adaptant pas forcément exactement avec ces points, il a fallut effectuer une interpolation à partir des données expérimentales afin de trouver un profil à injecter qui convienne à notre maillage.

Malheureusement, ce problème en entraîna un autre : avec un profil de vitesse non-nulle en entrée et une vitesse nulle dans tout le reste du domaine, OpenFOAM s'avère incapable de s'adapter, et finit par s'arrêter sur ce qu'il nomme une : "discontinuity error". Cela impliquera dans notre cas de passer par un autre solveur, potentialFoam, pour obtenir des champs de vitesse dans tout le domaine, champs que OpenFOAM pourra accepter. Il nous faudra finalement injecter ces données dans le cas en cours pour parvenir à faire marcher la simulation.

 

Finalement, tout semble marcher pour le mieux lorsque l'on trace les profils de vitesse à l'entrée :

 

 

Il faut finalement accepter le fait que potentialFoam déforme le profil de vitesse à l'entrée (qui, après interpolation, se superpose exactement avec le valeurs expérimentales qui sont les conditions limites imposées).

Ce dernier problème n'a pu être résolu, en partie parce qu'il ne fut découvert que tardivement. Il faut donc garder à l'esprit que les simulations suivantes présenteront à la base ce genre de défaut ...