Mise en place d'un calcul test

 

Avant de nous lancer dans la simulation de l'écoulement du train dans un tunnel, nous avons décidé de réaliser une simulation préliminaire sur un cas simple dont on connaît le comportement physique.

Cela va nous permettre de valider les solveurs et modèles choisis sous OpenFoam et surtout de nous assurer que nous avons pris en main  correctement ce logiciel.

 

1. Description physique

Physique du problème

Pour ce calcul test, nous avons décidé de simuler l'écoulement d'air bidimensionnel derrière une marche :

La conduite est de section carrée et on doit observer à l'arrière de la marche une zone de recirculation puis un recollement de la couche limite à la paroi.

 

Modélisation de la turbulence

Cet écoulement est turbulent et on utilise le modèle $k-\epsilon$.

Afin de se rapprocher de la publication, on a choisi un $Re=8000$ soit une vitesse entrante de :

$U=\frac{Re \nu}{D}=\frac{ 8.10^{3}.1,56.10^{-5}}{0,1}=1,256 m.s^{-1}$.

On calcule les paramètres de turbulence :

$k=\frac{3 }{2} (U_{\infty} I)^2$     

$\epsilon=\frac{C_{\mu}^{0,75} k^{1,5}}{l}$

avec l'échelle de longeur turbulente $l=0,07 . D_h = 0,07 . \frac{4A}{P} = 0,007 m$ ici

et $C_{\mu}$ une constante du modèle de turbulence égale à $0,09$.

En prenant une intensité turbulente $I=5\%$, on obtient le couple de valeurs suivant :

$ \boxed{k=0,0059 m^{2}.s^{-2}   ,   \epsilon=0,0106 m^{2}.s^{-3}} $

 

Maillage du domaine

On a utilisé IcemCFD pour réaliser le maillage de l'écoulement derrière une marche. Nous avons choisi un $y^{+}$ (nombre adimensionnel de la distance à la paroi) de $35$. Ce choix est fait pour se mettre dans les mêmes conditions que l'étude du train. En effet, le maillage que nous ont fourni les industriels partenaires pour le cas du train a un $y^{+}$ supérieur à $30$ et nécessitera donc l'utilisation de lois de paroi (zone log de la couche limite). Grâce à l'application pointwisenous avons déterminé que pour un $y^{+}$ de $35$ il nous faut une première maille à $0,007m$ de la paroi.

Le maillage final est le suivant :

 

 

Le maillage est cartésien et comporte $5000$ cellules. Il est légèrement raffiné dans les zones d'intérêt (proche parois et zone de recirculation).

 

On a défini les conditions aux limites suivantes :

-  Velocity Inlet en entrée avec une vitesse de $1,256 m.s^{-1}$

- Pressure Outlet en sortie avec une pression fixée à zéro.

- Wall pour les murs haut et bas

 

La qualité du maillage a été contrôlée via IcemCFD. L'aspect ratio du maillage varie entre $1$ et $8,8$. Cette valeur maximale est très convenable.

2. Simulation avec OpenFoam

On va maintenant voir les différentes étapes à réaliser pour mettre en place notre cas avec OpenFoam.

 

Structure générale

Tout d'abord, il faut respecter une certaine structure des dossiers. Ainsi l'utilisateur doit créer un fichier de calcul où tous les résultats seront générés.

L'organisation du code est la suivante :

source : BEI FEP 2012-2013 Gr 21

 

Tout d'abord, il faut choisir le solveur que l'on va utiliser. C'est un choix que l'on doit faire en concertation avec le groupe 23 qui travaille sur la partie maillage mobile de la simulation du mouvement du train dans un tunnel.

Le solveur choisi est pimpleDyMFoam qui permet de gérer les maillages mobiles (DyM pour Dynamic Mesh).

C'est un solveur pour les écoulements transitoires incompressibles qui utilise l'algorithme PIMPLE pour la résolution des équations de Navier Stokes basées sur la pression (combinaison des algorithmes PISO et SIMPLE).

OpenFoam propose des tutoriels, calculs déjà prêts à être lancés. Il est donc utile de s'appuyer sur un tutoriel correspondant au solveur choisi. Pour le solveur pimpleDyMFoam, nous avons utilisé le tutoriel movingCone.

Par contre, aucun modèle de turbulence n'est utilisé dans ce tutoriel, on va donc devoir le mettre en place. Il va également falloir importer le maillage de l'écoulement derrière une marche cité précédemment.

 

Importation du maillage

Sous OpenFoam on peut importer un maillage fait sous IcemCFD via la commande fluentMeshToFoam. Cela génère alors tous les fichiers de description du maillage dans le dossier constant/polyMesh.

Pour vérifier que l'importation s'est bien déroulée et que la qualité du maillage importée est satisfaisante, on peut utiliser la commande checkMesh qui fournit notamment les informations suivantes sur notre maillage :

Les différentes conditions aux limites ont bien été reconnues. Il faut cependant vérifier dans le fichier constant/PolyMesh/boundary qu'elles ont été définies avec le bon type (INLET et OUTLET doivent être de type Patch...)

On notera que sous OpenFoam, les maillages 2D sont extrudés dans la troisième direction et une zone frontAndBackPlanes correspondante est alors créee.

On retrouve ensuite des informations similaires quant à la qualité du maillage que celles fournies par IcemCFD.

Nous allons maintenant détailler la modification des différents fichiers pour la mise en place de notre calcul.

 

Mise en place de la turbulence

Dans le tutoriel movingCone, la turbulence n'est pas utilisée. Pour la simuler dans notre calcul, il va falloir copier les fichiers turbulenceProperties et RASProperties dans le dossier constant  de notre cas.

On a copié ces fichiers à partir d'un autre tutoriel prenant en compte la turbulence. On les modifie ensuite de manière à ce qu'ils décrivent bien le modèle $k-\epsilon$ que l'on a décidé d'utiliser. Ainsi le module RASModel doit être choisi dans le fichier turbulenceProperties et le module kEpsilon dans RASProperties.

Enfin, il faut créer des fichiers k et epsilon dans le dossier temporel 0 pour que ces grandeurs soient calculées.

 

Mise en place des conditions initiales

On se place maintenant dans le dossier temporel 0 et on modifie les fichiers correspondant à chaque champ en imposant les modules suivants aux limites :

 

CL / fichiers U p k epsilon
Velocity Inlet "fixedValue" "zeroGradient" "fixedValue" "fixedValue"
Pressure Outlet "zeroGradient" "fixedValue" "zeroGradient" "zeroGradient"
Wall "fixedValue" "zeroGradient" "kqrWallFunction" "epsilonWallFunction"

- "fixedValue" permet à l'utilisateur de fixer une valeur (on rentre par exemple pour k la valeur calculée précédemment, une valeur nulle pour la vitesse aux murs...)

- "zeroGradient" permet d'imposer un gradient perpendiculaire nul sur la condition limite

- "kqrWallFunction" et "epsilonWallFunction" sont des lois de paroi. Notre $y^{+}$ étant supérieur à $30$, l'utilisation de lois de parois est nécessaire. Leurs valeurs doivent normalement être fixées à zéro. Mais d'après les informations de la communauté OpenFoam, donner pour ces lois de parois la même valeur qu'en entrée permet d'augmenter la convergence. C'est ce que nous avons fait par la suite.

- On remarque que OpenFoam va automatiquement générer un fichier nut lors du calcul à partir des grandeurs k et $\epsilon$ fournies.

- La zone frontAndbackPlanes dûe à l'extrusion 3D du maillage doit être imposée comme "empty".

 

Paramètres numériques

Ceux-ci sont à régler dans le dossier system, qui comprend quatre fichiers :

  1. fvSchemes : fichier dans lequel l'utilisateur choisit les méthodes de résolution des différents opérateurs mathématiques présents. Ces choix sont déterminants pour la stabilité et la précision des résultats :

- Le premier choix à faire est celui de la discrétisation temporelle. Ici on choisit le schéma implicite d'Euler.

- Ensuite, il faut choisir les discrétisations spatiales associées aux différents termes résolus par les différentes équations.

Dans un premier temps, nous avons choisi d'utiliser des schémas centrés d'ordre 2 pour la plupart des opérateurs, plus précis que les schémas amonts. (à titre de comparaison, sous Fluent la plupart des schémas sont par défaut des schémas amonts). Ici, on sélectionne des schémas centrés via le dictionnaire Gauss linear.

Pour ce calcul, notre fichier fvSchemes  est le suivant.

 

  1. fvSolution : fichier dans lequel l'utilisateur choisit les critères de convergence des différentes méthodes numériques utilisées. Ces paramètres sont importants pour augmenter la convergence du calcul et sa rapidité.

Ici,  pour la résolution de pression on a choisit d'utiliser la méthode des gradients conjugués pré-conditionnés,  (PCG sous OpenFoam). Comme pré-conditionneur on a choisi FDIC (pour Faster diagonal incomplete-Cholesky)  qui est supposé plus efficace lors de calculs parallèles que le pré-conditionneur de type DIC. Cette méthode numérique étant valable uniquement pour les matrices symétriques, la résolution en vitesse (matrices asymétriques) sera effectué via une résolution par gradients bi-conjugués (PBiCG sous OpenFoam).

Par rapport aux différents tutoriels OpenFoam dont on a pu s'inspirer, les tolérances des différentes méthodes numériques ont été abaissées afin d'assurer la convergence des calculs. En effet, les paramètres de tolérance absolue tolerance et de tolérance relative relTol ont été fixés respectivement à $10^{-10}$ et $0,01$.

De plus, afin de stabiliser le calcul, nous avons décidé de sous-relaxer la pression lors de la résolution de cette dernière via l'algorithme PIMPLE, et fixé le facteur de sous-relaxation à $0,3$ (au détriment du temps de calcul).

Le contenu du fichier fvSolution est visible ici.

 

  1. controlDict : fichier dans lequel l'utilisateur donne le solveur utilisé, les grandeurs temporelles du calcul telles que le temps initial, le temps final, le pas de temps et l'intervalle d'écriture des solutions.

 

Pour ce calcul, notre fichier controlDict est de la forme suivante :

La fonction adjustTimeStep permet de fixer un critère maximum de stabilité et ajuste le pas de temps en conséquence.

Après plusieurs simulations, le critère de stabilité à été fixé ici à $1,5$ afin d'obtenir un bon compromis entre temps de calcul et précision des résultats.

 

  1. decomposeParDict : fichier qui permet de gérer le découpage du maillage et du champ initial en plusieurs sous-domaines de calcul pour lancer la simulation en parallèle.
    Ici, deux sous domaines de calcul ont été définis. Cela permet dans notre cas de diminuer le temps de calcul par environ $1,75$, ce qui sur notre maillage à faible nombre de cellules ne représente qu'un peu plus d'une minute mais qui sera considérable par la suite sur le maillage autour du train.

 

On lance alors le calcul via les commandes suivantes (dans le dossier courant):

  • decomposePar, qui répartit les sous-domaines aux différents processeurs (2 dossiers sont alors crées)
  • mpirun -np 2 pimpleDyMFoam -parallel >> log & pour lancer le calcul en parallèle et stocker les résultats dans un fichier log
     

 

3. Résultats et validation

Résumé des paramètres de la simulation : 

Reynolds en entrée : $8000$

Modèle de turbulence $k - \epsilon$ avec :

$k=0,0059 m².s^{-2}$ et $\epsilon = 0,0106 m².s^{-3}$ 

CFL : $1,5$

Pas de temps : $0,0065 s$

Temps physique simulé : $20s$

Temps CPU : $100s$

 

Résultats 

On peut visualiser les résidus via les outils des librairies pyFoam grâce à la commande pyFoamPlotWatcher.py log :

(Pour comparaison, voici les résidus sous les mêmes conditions mais avec des schémas amonts)

On reconstruit ensuite les champs complets des différentes grandeurs avec reconstructPar pour pouvoir les visualiser sous paraFoam.

On s'assure tout d'abord que l'écoulement en sortie du canal d'entrée ($x=0.9 m$) est pleinement turbulent.

 

 

 

 

Le sommet du profil est aplati contrairement au profil parabolique de Poiseuille pour un écoulement laminaire. En effet, en sortie du canal d'entré, l'écoulement turbulent est pleinement établi, la couche limite est plus épaisse qu'en régime laminaire. Ceci est dû à la formation de petits tourbillons à l'intérieur de la couche limite qui entraîne un brassage de l'écoulement et a tendance à uniformiser les vitesses.
On observe de plus que la vitesse maximale est supérieure à la vitesse imposée de $1.256 m.s^{-1}$ à cause des fluctuations de vitesses induites par le régime turbulent.

 

 

 

 

 

Ensuite, on peut observer, sur les résultats ci-dessous, l'établissement du champ de vitesse vers l'état stationnaire :

On retrouve le comportement physique attendu avec une zone de recirculation puis une zone de recollement. En effet, le gradient de pression adverse dû au changement de section induit un écoulement en sens opposé qui va provoquer la séparation de la couche limite. C'est la zone de recirculation (voir courbe ci-dessous, en $x=0,1 m$ après le changement de section).

 

 

À l'aval de cette zone, la couche limite se rattache à la paroi en un point appelé point de recollement (noté $x_1$) comme on peut le voir avec le tracé des lignes de courant ci-dessous.

 

Validation

Pour valider de manière plus précise les résultats obtenus avec OpenFoam, nous allons comparer la longueur de recollement obtenue avec la littérature.

D'après Numerical Study using Fluent of the Separation and the Reattachment Points for Backwards-Facing Step Flow de Luke Jongebloed [B2], pour les mêmes conditions que notre cas ($Re = 8000$, hauteur de la marche $S=0,094 m$), la grandeur caractéristique adimensionnelle de recollement est $\frac{x_1}{S}=6,8$.

 

Pour déterminer numériquement le point de recollement, il faut trouver l'endroit où le coefficient de frottement s'annule à la paroi.

Tous calculs faits, on utilise la commande wallShearStress dans le dossier courant, pour obtenir le champ de cette grandeur.

Ainsi, on obtient en régime stationnaire :

 

 

 

 

On trouve un point de recollement situé à $x_1=0,63 m$ du changement de section.

Les résultats de notre simulation donnent une valeur de $\frac{x_1}{S}$ de $6,7$. 

Cette valeur est très proche de la valeur théorique ce qui valide notre étude d'autant plus que le modèle $k-\epsilon$ n'est pas le meilleur pour la prédiction de la  séparation/du recollement de couches limites. Cet écart à la théorie peut aussi provenir de l'aspect numérique du problème. En effet, l'interpolation des différentes grandeurs a pu être faussée de par la taille des cellules, relativement élevée.

 

 

 

 

 

 

Les résultats de ce cas test sont donc validés. La simulation de la turbulence avec le modèle de $k-\epsilon$ pour le solveur pimpleDyMFoam a bien fonctionné.

De plus, le choix d'utiliser des schémas centrés était judicieux puisqu'à titre de comparaison, avec les schémas amont d'ordre deux,  on obtient une valeur de $7,02$ pour la longueur adimensionnelle de recollement. Cette valeur atteste de la diffusion numérique prononcé de ces schémas..

Ce premier cas test nous a aussi permis de prendre en main en main le logiciel et on peut maintenant se lancer dans la simulation du mouvement d'un train dans un tunnel.