Simulation du mouvement d’un train avec le logiciel OpenFOAM - Comparaison de modèles de turbulence

Simulation du mouvement d’un train avec le logiciel OpenFOAM

Comparaison de modèles de turbulence

Adrien Emery & Nicolas Verdin

Encadrants INP :  Dominique Legendre & Annaig Pedrono

Encadrants industriels : Sydney Tekam (Tech-Am Ingénierie) & André-Marie Dogbo (Setec)

 


      

Introduction et objectifs

Ce sujet consiste à modéliser sous OpenFoam le passage d'un train dans un tunnel et d'en étudier les différents phénomènes physiques. Le projet est commun aux groupes 23 et à notre groupe (24). Ici, nous nous concentrerons sur l'aspect turbulent de l'écoulement.

 

Présentation de l'entreprise

Ce projet est réalisé en collaboration avec la société Tech-Am Ingénierie dont le seul membre est Sydney Tekam, le créateur de cette entreprise.

C'est une société de conseil en ingénierie qui réalise des missions pour de plus grands groupes. Ainsi, la branche BTP de l'entreprise Setec a mandaté Tech-Am Ingénierie pour réaliser une étude sur le mouvement d'un train dans un tunnel.

 

Présentation du sujet

Le passage d'un train dans un tunnel engendre des écoulements d'air complexes (structures tourbillonnaires, fluctuations de pression importantes dans le tunnel...). Or les fortes variations de pression au passage d'un tunnel sont douloureuses pour les passagers, notamment au niveau des tympans. Ces fluctuations de pression engendrent également des efforts sur les parois du tunnel et les cabines du train et en augmentent la fatigue mécanique.

La bonne prédiction des écoulements d'air dans le tunnel au passage d'un train est donc un enjeu essentiel tout autant pour dimensionner les trains (acoustique, aération...) que pour optimiser la construction des tunnels. L'objectif final étant d'atteindre une sécurité maximale et un confort optimal des passagers.

 

Présentation du logiciel OpenFoam

OpenFoam (Open Field Operation and Manipulation) est un code de calcul CFD, sans interface graphique, développé en C++. C'est un logiciel OpenSource depuis 2004.

L'aspect objet du C++ permet une très grande modularité d'utilisation. Contrairement aux autres solveurs commerciaux, l'utilisateur a la possibilité de contrôler tous les différents aspects de la simulation : solveurs, opérateurs, grandeurs et de les assembler selon ses besoins.

OpenFoam est donc un outil puissant et adaptable à de nombreuses applications.

Mais cette liberté de contrôle conduit à une certaine opacité du code de par sa prise en main complexe.

 

Objectifs du projet

Nous allons donc réaliser l'étude CFD de ce problème sous OpenFoam et notre objectif principal sera la modélisation de la turbulence. On pourra de plus évaluer les performances d'OpenFoam par rapport à des logiciels payants comme Fluent.

Nous avons organisé notre étude de la manière suivante :

  1. Choix des modèles de turbulence
  2. Mise en place d'un calcul test
  3. Application au mouvement du train dans un tunnel
  4. Analyse des résultats
  5. Étude paramétrique
  6. Simulation avec maillage mobile

 

Choix des modèles de turbulence

La première étape de notre travail a été de choisir les modèles de turbulence avec lesquels nous allons travailler.

Tout d'abord, il faut déterminer si nous nous trouvons en compressible ou en incompressible.

La limite de compressibilité se trouve pour un nombre de Mach proche de $0,3$. En considérant la vitesse du son à $340 m.s^{-1}$, on trouve une vitesse limite de $102 m.s^{-1}$ à partir de laquelle on doit considérer la compressibilité de l'air. Mais après contact avec les industriels, il s'avère que la vitesse du train que nous aurons à simuler sera bien inférieure à cette limite.

Notre cadre d'étude sera donc incompressible.

Ensuite, il nous a fallu choisir l'approche nécessaire à la résolution de la turbulence, à savoir RANS, LES, DES ou DNS. Le choix a été très rapide, compte tenu des ressources de calcul disponibles à l'ENSEEIHT : c'est donc l'approche RANS qui a été choisie et où la turbulence est modélisée.

Une première recherche sur le site officiel d'OpenFoam nous a permis de connaître l'ensemble des modèles de turbulence disponibles sous OpenFoam. Pour l'approche RANS en incompressible, notre cas d'étude, il y en a quinze.

Pour faire un premier tri parmi cette multitude de modèles, nous avons réalisé une recherche bibliographique sur ces différents modèles.

L'article Evaluation of Turbulence Models for Prediction of Flow Separation at a Smooth Surface de Eric Furbo, Janne Harju et Henric Nilsson de l'université d'Uppsala [B1] nous a été d'une grande utilité.

En effet, ils comparent tous les modèles de turbulence d'OpenFoam sur un cas proche du nôtre : l'écoulement confiné derrière une bosse.

Configuration de l'étude d'Uppsala [1]

 

Ils constatent que 8 des modèles de turbulence présentent des problèmes de convergence voire divergent totalement et sont donc inadaptés à ce type de problème. Nous avons donc décidé de ne pas les tester dans notre étude.

Ils nous restent alors 7 modèles potentiellement utilisables dans ce projet. Dans le tableau suivant, nous allons résumer quelques caractéristiques de chacun :

 

Modèle Spécificité Pertinence
Non-linear Shih $k-\epsilon$ Modèle non linéaire complexe Modèle complexe, très peu d'informations en ligne
$k-\epsilon$ Standard

Facile d'utilisation, convergence assurée, "couteau suisse" de la turbulence.

Modèle à 2 équations de transport.

Obtenir un premier résultat facilement, pas forcément très précis. Permet la comparaison avec une étude analogue sous Fluent
$k-\epsilon$ RNG

Prend en compte de plus petites échelles de turbulence que $k-\epsilon$ standard.

Certaines constantes du modèle $k-\epsilon$ ne le sont plus.

Donne de meilleurs résultats pour les écoulements tournants et pour modéliser les cavités entraînées.

Pas d'intérêt particulier pour les autres configurations dont notre cas
Realizable $k-\epsilon$ Modèle différent pour l'équation de dissipation par rapport à $k-\epsilon$ standard Modèle $k-\epsilon$ amélioré mais on préférera utiliser le modèle standard pour la simplicité de ses résultats.
$k-\omega$ standard Corrige les problèmes de séparation des couches limites du modèle $k-\epsilon$ standard Intéressant à tester
$k-\omega$ SST

Combinaison du modèle $k-\omega$ près des murs et de $k-\epsilon$ au coeur de l'écoulement.

Modèle à 2 équations de transport.

Plus intéressant à tester que le $k-\omega$ standard.

Modèle complet et adapté à notre cas d'écoulement confiné.

Spalart-Allmaras

Conçu au départ pour l'aérodynamique externe.

Modèle à une équation de transport.

Facile d'utilisation (modèle à une équation) et peut donner de bons résultats. Conserver un regard critique car modèle non conçu pour étudier les écoulements confinés au départ.

 

On retiendra les 3 modèles suivants :

  •  $k - \epsilon$ standard ($1972$)
  •  $k - \omega$ SST​ ($1993$)
  • Spalart-Allmaras ($1992$)

 

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.

Application au mouvement du train dans un tunnel

Le but de notre étude va être d'étudier les phénomènes de turbulence à l'intérieur de cet écoulement.

1. Paramètres de l'étude et maillage

Dimension du problème

Nous allons simuler le mouvement d'un train de 80 mètres dans un tunnel d'une hauteur de 7 mètres et long de 300 mètres :

Le train se déplace de la droite vers la gauche.

Le maillage sera décomposé en trois zones principales : la zone en amont du train, la zone du train et le zone en aval du train.

Chacun de ces blocs comprendra des conditions aux limites murs pour les parois du toit et du sol du tunnel ainsi qu'une zone fluide interne.

La structure du train sera considérée comme une condition aux limites mur.

L'étude du mouvement d'un train dans un tunnel est différente d'une étude en soufflerie car le tunnel entraîne un effet "piston". L'air est confiné et ne peut se déplacer que le long des murs. Le train chasse l'air devant lui créant une succion : de l'air rentre dans le tunnel et s'écoule vers l'avant du train.

Cet effet piston est très prononcé puisque le train occupe quasiment entièrement la section du tunnel. Cela oblige à mettre une condition de vitesse en entrée, contrairement à ce que l'on pouvait premièrement penser par analogie avec les études en soufflerie. Mais ce phénomène ne peut pas être reproduit par un maillage fixe. On y reviendra dans la dernière partie de ce projet où l'on travaillera avec le maillage mobile.

Dans notre étude on sera donc dans le cas classique de la soufflerie. On imposera donc la vitesse en inlet et la pression en outlet.

​De plus, la personne travaillant chez Setec nous a demandé de considérer que l'air à l'intérieur du tunnel n'est pas au repos et qu'il y a un écoulement d'air de $15 m.s^{-1} à 23 m.s^{-1}$ de la sortie du tunnel vers l'entrée. Dans un premier temps, on fixera cet écoulement d'air à $20 m.s^{-1}$.

Setec souhaite une vitesse de $144 km.h^{-1}$ ( $40 m.s^{-1}$) pour le train. Pour simuler correctement ce problème, il faudra donc que l'on injecte en entrée du tunnel de l'air à la vitesse $40 - 20 = 20 m.s^{-1}$. 

On remarque qu'on est bien loin de la vitesse limite de compressibilité de $102 m.s^{-1}$ calculée en introduction.

Calculons le nombre de Reynolds de l'écoulement : 

$\boxed{Re = \frac{U_{\infty} . L}{\nu} = \frac{20 . 7}{1,56.10^{-5}} = 8,9 . 10^{6} }$

 

Qualité du maillage

La société Tech-Am Ingénierie nous a fournis un maillage que l'on a importé sous OpenFoam grâce à la commande fluentMeshToFoam.

La commande checkMesh d'OpenFoam nous a fournis les informations suivantes : 

 

 

Le maillage comporte 87 460 cellules dont 87 402 hexahèdres et 58 triangles. Les mailles triangulaires sont généralement de faibles qualités, c'est donc une bonne chose qu'elles soient très rares. Elles sont présentes près du train où les formes sont complexes.

Le checkMesh confirme également que l'importation du maillage a été bien faite, les dix conditions aux limites ayant bien été reconnues. Et la vérification du fichier constant/PolyMesh/boundary confirme qu'elles ont été définies par le bon type (Patch pour entree_tunnel et sortie_tunnel...)

D'autres informations sur la qualité du maillage sont également présentes, notamment sur la dissymétrie des mailles (skewness) ou la non orthogonalité.

Avec Paraview, plusieurs de ces aspects de la qualité du maillage ont pu être visualisés : 

- La distorsion des mailles grâce au Scaled Jacobian

Le Jacobien informe sur la distorsion des mailles, il est compris entre $0$ et $1$. Ainsi il vaut $1$ pour les carrés et rectangles et sa valeur diminue quand la distorsion de la maille augmente.

Un maillage est considéré comme de bonne qualité si toutes ses cellules ont un Jacobien supérieur à $0,6$.

On remarque que seuls les mailles triangulaires ont un Jacobien inférieur à $0,6$. Leur faible nombre permet de valider la bonne qualité du maillage.

 

- Le volume des mailles

 

Il y a une grande discontinuité du volume des mailles entre la zone en amont du train et la zone du train. Cette discontinuité nuit à l'interpolation des résultats sur cette interface. Un raccordement plus progressif améliorerait la précision des résultats.

Ce problème apparaît moins à la jonction avec la zone en aval du train où le raccordement entre les deux blocs est moins abrupt. C'est une bonne chose car c'est dans cette région à l'arrière du train que les phénomènes les plus complexes (traînée, recirculation...) pourront être observés.

Malgré ce problème de différence de taille des mailles et la présence de quelques mailles triangulaires autour du train, il est bien raffiné autour de la géométrie complexe du train, c'est un maillage de bonne qualité.

2. Mise en place sous OpenFoam

Pour le modèle de turbulence $k-\epsilon$

La structure des dossiers et des fichiers sera quasiment identique à celle du cas test de l'écoulement derrière une marche.

Le maillage importé est bien évidemment différent et donc le dossier polyMesh aussi comme détaillé précédemment.

La valeur de la vitesse en entrée est différente, ce qui va également modifier les paramètres de turbulence. Ainsi pour une vitesse en entrée du tunnel de $U_{\infty} = 20m.s^{-1}$, les grandeurs turbulentes k et $\epsilon$ seront telles que :

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

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

avec $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=1\%$, on obtient le couple de valeurs suivant :

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

 

Ce sont les seules modifications sous OpenFoam par rapport au cas test de l'écoulement derrière une marche.

 

Pour le modèle de turbulence $k-\omega  SST$

Il n'y a pas de tutoriel sur ce modèle avec le solveur PimpleDyMFoam. Mais dans le tutoriel tutorials/incompressible/SimpleFoam/motorbike il est utilisé avec le solveur SimpleFoam. On constate que le modèle $k-\omega  SST$ utilise trois fichiers k, nut et omega dans le dossier 0.

On copie alors le fichier omega dans le cas $k-\epsilon$ précèdent pour mettre en place le modèle $k-\omega  SST$. Comme pour le cas $k-\epsilon$, OpenFoam va créer automatiquement le fichier nut à partir des données des fichiers k et omega.

On modifie les conditions aux limites de omega en les calquant sur celles de k. Leurs lois de parois diffèrent, il faut utiliser omegaWallFunction  pour omega.

Il faut également modifié le fichier RASProperties en mettant KOmega SST pour RASModel.

On reprend la valeur de $k$ de l'étude avec le modèle $k-\epsilon$ (il se calcule de la même manière).

$\omega$ est calculé de la manière suivante : 

$\omega = \frac{\epsilon}{k.\beta^{*}}$ avec $\beta^{*} = C_{\mu} = 0,09$

$\omega= \frac{0,35}{0,06.0,09}$

D'où : $\boxed{\omega = 64,8 s^{-1}}$

 

Une fois ces modifications faites, le modèle de turbulence $k-\omega  SST$ peut être simulé. 

 

 

Pour le modèle de turbulence $Spalart-Allmaras$

Ce modèle n'est pas non plus utilisé dans un tutoriel avec le solveur PimpleDyMFoam mais il est présent dans le tutoriel tutorials/incompressible/SimpleFoam/airfoil2D

On constate que le modèle de $Spalart-Allmaras$ utilise deux fichiers dans le dossier 0nut et nuTilda. On copie donc ces fichiers dans le dossier  de notre cas de travail Spalart Allmaras.

Le modèle $Spalart Allmaras$ consiste en une équation de transport d'une viscosité turbulente modifiée $\tilde{\nu}$ aussi appelée variable de Spalart-Allmaras. 

​On estime la valeur de $\tilde{\nu}$ grâce à la formule suivante : 

$\tilde{\nu} = \sqrt{\frac{3}{2}} . U.I.l$

avec $I$ l'intensité turbulente prise égale à $1$%

et l'échelle de longueur turbulente $l=0,007 m$ déjà calculée lors de la partie précédente

D'où : $\tilde{\nu} = 0,0017 m².s^{-1}$

 

Puis la définition du modèle $Spalart-Allmaras$ nous donne :

$\nu_t = f_{v1} . \tilde{\nu}$

$f_{v1} = \frac{\chi^{3}}{\chi^{3} + C_{v1}^{3}}$

avec $\chi = \frac{\tilde{\nu}}{\nu}$ et $C_{v1} = 7,1$ une constante du modèle

Pour l'air $\nu = 1,56.10^{-5} m².s^{-1}$, d'où  $f_{v1} = 0,999$

On obtient donc :  $\boxed{\nu_t = \tilde{\nu} = 0,0017 m².s^{-1}}$

 

 

Pour finaliser la mise en place du modèle $Spalart-Allmaras$, il faut modifier les conditions aux limites des fichiers nut et nuTilda​.

Aux murs, on utilise la loi de paroi nutUWallFunction​ avec comme valeur celle calculée en entrée.

Il ne faut également pas oublier d'adapter le fichier constant/RASPropeties  en précisant SpalartAllmaras  en face de RASModel.

 

Convergence des modèles

Afin d'augmenter la stabilité des solutions, nous avons ajouté, par rapport au cas de la marche la sous-relaxation des grandeurs physiques avec les facteurs suivant :

U ($0,7$), k , epsilon, omega, nuTilda ($0,8$).

(La pression étant déjà sous-relaxée depuis le cas de la marche, avec le facteur $0,3$).

Au niveau du choix des schémas de discrétisations spatiales, nous avions précédemment utilisés des schémas centrés. Pour le cas du train nous avons utilisé en grande partie des schémas amonts et ce pour plusieurs raisons :

- Même si ces schémas sont moins précis que les schémas centrés, ils introduisent naturellement de la diffusion numérique qui vient stabiliser le calcul,

- Les schémas centrés sont plus dispersifs et introduisent des oscillations parasites hautes fréquences ("wriggles"),

- Légèrement plus rapide (car stencil plus petit), même si cela ne se ressent pas beaucoup sur le temps total de la simulation,

- Dans un soucis de comparaison, les schémas amonts ont été utilisés car ce sont ceux qui sont en commun avec le logiciel Fluent. C'est la raison principale de notre choix, car le souci de stabilité évoqué au premier point est enlevé ici du fait qu'on utilise un schéma implicite pour la discrétisation temporelle.

 

Pour vérifier cela, nous avons lancé une simulation avec le modèle $k-\epsilon$ et des schémas centrés, dont voici les résidus :

On peut voir que le schéma n'est pas très stable (résidus pour des schémas amonts sous les mêmes conditions, avec le modèle $k-\omega$).

Même si en moyenne la convergence est obtenue pour les schémas centrés, les schémas amonts sont choisis.

 

Alors que sur le cas de la marche, nous avions fixé la valeur du nombre de Courant à $1,5$, ce qui nous permettait de lancer nos simulations avec un pas de temps variable, ici nous avons choisi de fixer le pas de temps afin d'avoir un contrôle maximale sur la simulation.

Après plusieurs tests, nous avons fixé la valeur du pas de temps à $0,001 s$, valeur permettant la convergence des trois modèles pour un temps de calcul raisonnable.

Pour le modèle $k-\omega  SST$, on a les fichiers suivants : fvSchemes et fvSolution.

 

Résumé des paramètres des simulations : 

 

Reynolds en entrée : $8000$

 

Modèles de turbulence :

- $k- \epsilon$ avec $k = 0,06 m².s^{-2}$ et $\epsilon = 0,35 m².s^{-3}$

- $k - \omega  SST$ avec $k = 0,06 m².s^{-2}$ et $\omega = 64,8 s^{-1}$

- $Spalart-Allmaras$ avec $\nu_{t} = \tilde{\nu} = 0,0017 m².s{-1}$

 

Pas de temps : $1.10^{-3} s$ pour tous les modèles06Pas de temps : kNombre : 

Temps physique simulé : $20s$

 

3. Résultats généraux et vérification physique

Résultats généraux

 

Dans cette partie on ne va pas chercher à comparer les résultats obtenus suivant les différents modèles de turbulence mais on va seulement s'attacher à décrire le comportement physique général.

Ci dessous, on peut observer l'établissement du champ de vitesse à l'arrière du train :

On observe le développement de tourbillons et d'une traînée à l'arrière du train. Aux alentours de $11s$, l'écoulement est établi et le sillage derrière le train est fixe.

Par la suite on se placera toujours en régime établi, au temps t = $20s$,  pour observer les différentes grandeurs de l'écoulement

 

Ainsi, à partir du régime établi, on obtient le champ de pression suivant au niveau du nez du train :

Une forte dépression apparaît au niveau de la courbure qui délimite le nez du train au début du toit, c'est donc ici que les vitesses sont les plus élevées. A contrario, au niveau du nez du train (zone la plus à gauche), une surpression est présente, laissant deviner un point d'arrêt (comme on peut le visualiser sur le tracé des vecteurs vitesses ci-dessous).

Un gradient de pression se met alors en place autour du nez du train (environ $1600 Pa$ sur $4m$, soit $400Pa.m^{-1}$) ce qui reste assez faible dans nos considérations (écoulement incompressible et vitesse de $20m.s^{-1}$) mais qui peut devenir important dans un cas plus réel (où on est proche de la limite de compressibilité de l'air et où ce gradient de pression va créer des chocs, en particulier lors de l'entrée dans le tunnel). C'est pour cela que sur les trains à grande vitesse, les ingénieurs travaillent sur la forme du nez du train afin d'obtenir une meilleure répartition de ce gradient pression (qui en plus des avantages aérodynamiques, permet aussi de diminuer le bruit).

Au niveau des encoches (séparation entre wagons), des zones de recirculation se mettent aussi en place. Des zones de vitesses quasiment nulles cohabitent avec des zones de survitesses, d'où un intense gradient de vitesse local.

A l'arrière du train, on observe l'établissement du champ d'énergie cinétique suivant :

Pendant les 5 premières secondes, on observe une évolution aléatoire de ce champ avec la création de tourbillons. Ensuite, le champ d'énergie cinétique turbulente se stabilise jusqu'à la formation d'un zone plus large. On peut en voir les effets sur le tracé des lignes de courant ci-dessous : 

On observe l'établissement de trois zones de recirculation disposées de telle manière que deux points de recollements apparaissent sur la paroi inférieur.

Les zones de recirculation après un obstacle apparaissent le long des parois à cause des gradients de pression. Ce sera une zone intéressante à étudier, notamment pour confronter les modèles de turbulence et les résultats de Fluent. En effet, on va pouvoir par la suite comparer les abscisses des deux points de recollement pour les différents modèles/logiciels. 

 

Pour donner une idée des profils de vitesses nous avons utilisé la fonction Plot over line disponible sous Paraview en plusieurs abscisses le long du tunnel, une fois le champ de vitesse établi. Nous avons exporté ces données au format .csv et tracé ces profils sous Malab.

 

Les abscisses et ordonnées ont été translaté de manière à ce que le premier point de mesure soit situé en $(0 ;0)$. Les profils de vitesse ont été adimensionnalisés par la valeur $U_{\infty} = 40m.s^{-1}$ (par souci de visualisation). On alors obtient les profils suivants à l'amont du train :

On constate une augmentation de la vitesse au dessus et en dessous du train qui s'explique par la conservation du débit

Le nez du train pousse l'écoulement d'air vers l'avant. La séparation entre wagons entraîne également une recirculation de l'air.

En ce qui concerne le raccordement du maillage entre la zone amont et la zone du train (avec la discontinuité du volume des mailles évoquée précédemment), on peut évaluer son impact sur un profil de vitesse axial sur une ligne tracée entre ces deux raccordements :

On observe une discontinuité de $0,04 m.s^{-1}$ lors du passage de la zone à grandes mailles vers la zone à mailles plus petites. Cela est respectable dans la mesure où cela ne représente que $0,2\%$ de la valeur de la vitesse autour de cette discontinuité (à cette hauteur là soit $y=4 m$), mais cela reste une perte de précision que l'on aurait pu éviter.

 

Vérification physique

 

Pour vérifier la physique de nos résultats, nous avons uniquement trouvé le comportement de la pression à la paroi supérieure du tunnel lors du passage d'un train.

Sous Paraview on définit donc la ligne qui correspond à cette paroi du tunnel (outil Plot over line) pour obtenir le profil suivant :

 

On observe un saut de pression à l'avant du train puis une baisse de pression tout le long du train. C'est un profil de pression typique d'un train dans un tunnel. On retrouve ce comportement et son explication physique dans la littérature ([B5], [B6]).

On retrouve ainsi le profil de pression obtenu, pour un train allant en sens inverse du nôtre : 

Profil de pression le long du tunnel issu de l'article [B5]

 

Néanmoins, on peut seulement valider un comportement qualitatif. Il n'y a pas de mesures expérimentales correspondant aux paramètres de notre étude, on ne peut donc pas valider l'ordre de grandeur de nos résultats mais seulement l'aspect qualitatif.

L'article [B6] apporte des explications physiques sur la forme de ce profil de vitesse en reliant pression et traînée. Ainsi on distingue des régions où se produit essentiellement de la traînée de forme (ou de pression) et des régions où se produit essentiellement de la traînée de friction :

Profil de pression le long du tunnel issu de l'article [B6]

 

La traînée de forme provient de la distribution de pression générée par l'écoulement autour du train. L'action de ces pressions se fait de façon normale à la surface du train et se ressent par deux grands sauts de pression au niveau de l'avant et de l'arrière du train.

Le passage de l'air le long des parois du train crée directement de la friction. C'est cette friction qui est à l'origine de la baisse de pression le long du train. Ce gradient de pression a pour conséquence une différence de pression entre l'amont et l'aval du train dont l'intensité peut atteindre plusieurs kilopascals (600 Pascals dans notre cas). On observe bien ce gradient de pression entre l'avant et l'arrière du train.

À noter qu'on retrouve dans le profil que nous avons obtenu l'influence des trois encoches matérialisant la séparation des wagons et qui entraînent des fluctuations de pression.

Les résultats obtenus rendent donc bien compte de la physique du problème et on retrouve un comportement décrit à plusieurs reprises dans la littérature.

Analyse des résultats

Les résultats obtenus pour les trois modèles de turbulence choisis vont être comparés à ceux obtenus avec le logiciel commercial Fluent.

​Nous analyserons aussi les performances des trois modèles de turbulence étudiés.

1. Comparaison avec le logiciel commercial Fluent

Paramètres du calcul sous Fluent

 

Le solveur 2D Transient  a été utilisé pour cette simulation.

Nous nous sommes placés dans les mêmes conditions qu'une simulation OpenFoam où le modèle $k-\epsilon$ était utilisé avec une intensité de turbulence de $1$%. Nous avons essayé d'utiliser des paramètres numériques les plus proches possible de ceux utilisé sous OpenFoam (discrétisation temporelle implicite à l'ordre 1, schémas amonts...). Le pas de temps à été fixé à $0,001 s$ sous les 2 solveurs.

Un résumé des paramètres utilisés sous Fluent est à consulter ici.

Pour simuler $20$ secondes de temps physique, Fluent a mis un peu plus de $10$ heures contre $9$ heures sous OpenFoam. Cette différence de temps de calcul peut s'expliquer par le fait que les méthodes numériques utilisées sous les deux solveurs ne sont pas rigoureusement les mêmes (résolution de la pression sous Fluent via l'algorithme PISO et une méthode multi-grille, algorithme PIMPLE, inexistant sous Fluent et gradient conjugué pour la résolution sous OpenFoam) ainsi que les critères de convergence. Mais il faut aussi noter que les temps de calcul dépendent de l'occupation de la machine au moment de la simulation.

 

Résultats

 

Les résultats obtenus sont comparables à ceux d'OpenFoam :

 

On observe de nouveau le décollement des couches limites à l'aval du train et une zone de recirculation. La gamme de valeur est cohérente, la valeur maximale sous Fluent est de $56,13 m.s^{-1}$ contre $55,98 m.s^{-1}$ sous OpenFoam.

 

Pour la pression, on retrouve également le même comportement que sous OpenFoam. L'ordre de grandeur est également identique, même si les valeurs extrêmes sont plus importantes sous Fluent, $-1522,91 Pa$ contre $-1247 Pa$ et $568,83 Pa$ contre $493Pa$.

 

Le comportement et les ordres de grandeurs trouvés sont donc à peu près identiques pour OpenFoam et Fluent.

Mais on va maintenant comparer plusieurs grandeurs plus quantitativement pour mettre en évidence les différences obtenues selon les logiciels.

Pour ce faire, on a comparé les profils de différentes grandeurs (vitesse, pression, contrainte de cisaillement et énergie cinétique turbulente) à différents endroits du tunnel : 

 

 

 

Comparaison des profils de vitesse à l'arrière du train

 

 

On observe bien la recirculation identique pour les deux logiciels mais OpenFoam surestime légèrement la vitesse par rapport à Fluent (de moins de $9$%). Ce sont donc des écarts faibles, acceptables et qui ne permettent pas de départager OpenFoam et Fluent sur ce point. 

 

Comparaison des profils de pression

 

 

Sur le gradient de pression le long du train observe un écart assez important entre OpenFoam et Fluent de l'ordre de $25$%.

Mais dans le reste du tunnel, les deux profils de pression sont quasiment identiques.

Une explication à cette différence assez importante, mais présente seulement le long du train, peut être une gestion différente des lois de parois pour une valeur de $y^{+}$ importante. En effet, le rétrécissement de la zone de passage du fluide entraîne des survitesses au dessus du train. Le passage brusque de vitesse (de $20 m.s^{-1}$ à plus de $55 m.s^{-1}$) entraîne également une variation brutale du $y^{+}$ qui dépasse $500$ dans la zone au dessus du train. Ceci nous amène à penser que les lois de parois des deux codes ne se comportent pas de la même manière pour des valeurs extrêmes de $y^{+}$. Une solution pour vérifier cette hypothèse aurait été de raffiner le maillage dans cette zone critique mais nous ne disposions pas des fichiers nécessaires à sa modification.

Mais il est aussi possible que ce problème proviennent des choix des lois de parois que nous avons faits sous OpenFoam...

 

Comparaison de la contrainte de cisaillement

 

 

Le tracé de la contrainte de cisaillement le long du sol permet de déterminer les points de recollement de l'écoulement derrière le train (qui correspondent à une contrainte de cisaillement nulle).

Avec les deux logiciels on obtient bien deux points de recollement mais la position du second varie suivant le logiciel utilisé.

Le premier point de recollement situé à $10$ mètres derrière le train est capté de la même manière par Fluent et OpenFoam. Mais pour le second point on constate une différence de $2,3$ mètres entre les deux logiciels.

Nous ne disposons pas de données expérimentales correspondant à notre cas d'étude, il est donc délicat de dire qui d'OpenFoam ou de Fluent est le plus proche du second point de recollement réel.

On peut quand même noter que sous Fluent la contrainte de cisaillement n'est pas scrupuleusement nulle (elle descend jusqu'à $0,01$ Pascals).

 

Comparaison de l'énergie cinétique turbulente

 

Les résultats obtenus pour l'énergie cinétique turbulente sont différents avec Fluent. L'allure général du champ est la même comme on peut le voir ci-dessous. 

Mais lorsque l'on regarde plus en détail les profils d'énergie cinétique turbulente $k$ derrière le train et au niveau de la première encoche après le wagon de tête, on constate des différences entre les résultats de Fluent et d'OpenFoam : 

 

Derrière le train, le maximum d'énergie cinétique obtenue sous Fluent est inférieure de près de $10$% à OpenFoam.

Cette différence est encore plus marquée au niveau de la séparation des deux premiers wagons :

 

Cette fois-ci c'est avec Fluent que les valeurs d'énergie cinétique les plus importantes sont obtenues. Pour le pic principal, la valeur de Fluent est près de $30$% plus importante que celle d'OpenFoam et l'on approche les $50$% de différence sur le second pic.

 

Certes l'allure générale est la même mais il y a des différences significatives au niveau des valeurs obtenues.

De plus, OpenFoam majore l'énergie cinétique turbulente à l'arrière du train alors qu'au contraire il sous-estime fortement la turbulence entre les wagons par rapport à Fluent. C'est assez surprenant mais cela va dans le même sens que notre hypothèse lors de la comparaison des profils de pression le long du tunnel : les encoches matérialisant la séparation entre les wagons sont des zones où une recirculation de l'air a lieu et donc où des vitesses plus faibles subsistent. Encore une fois ce fort gradient de vitesses (vitesses quasiment nulles dans les encoches et des vitesses de plus de $50 m.s^{-1}$ au dessus du train) entraîne un fort gradient de $y^{+}$ dont la valeur maximale est suffisamment élevée pour que la limite de validité de la loi-log de la couche limite soit atteinte.

 

Conclusion

Les résultats obtenus sous OpenFoam sont donc globalement en accord avec ceux obtenus sous Fluent mais une comparaison avec un troisième code de calcul aurait permis de départager Fluent d'OpenFoam et de savoir lequel se rapproche le plus de la réalité au niveau des valeurs extrêmes.

2. Comparaison des différents modèles de turbulence

En terme de temps de calcul et de convergence

 

Le modèle $k-\omega  SST$ a été le plus lent en terme de temps de calcul, $13$ heures contre $8$ heures pour le modèle $k-\epsilon$ et un peu plus de $9$ heures pour le modèle $Spalart-Allmaras$.

​Malgré tout, le temps de calcul pour le modèle $Spalart-Allmaras$, qui ne résout qu'une équation turbulente est plus important que celui du modèle $k-\epsilon$ qui doit résoudre 2 équations de turbulence. De plus il y a 5 heures de différence entre les deux modèles $k-\epsilon$ et $k-\omega  SST$ qui ont le même nombre d'équations à résoudre (et exactement les mêmes conditions de simulations).

Mais ces temps de calcul sont à relativiser car ils dépendent de l'occupation de la machine au moment de la simulation.

 

La convergence des résidus a été observée sur les 3 modèles de turbulence, comme on peut le voir ci-dessous, pour le modèle $k-\omega  SST$ :

 

En terme de résultats

 

Pour comparer les trois modèles de turbulence choisis, nous avons tracé les profils de pression, de vitesse et déterminé les abscisses des points de recollement.

 

Ainsi, pour les modèles $k-\epsilon$, $k-\omega  SST$ et $Spalart-Allmaras$, le profil de pression le long du tunnel a le comportement suivant :

Les trois modèles permettent d'obtenir globalement le même comportement en pression avec un saut de pression à l'avant du train et un gradient de pression le long de celui-ci.

Néanmoins des différences sont présentes entre les modèles, tout d'abord en terme de valeurs obtenues où l'on peut avoir jusqu'à $7$% de différence entre les modèles $k-\omega  SST$ et $Spalart-Allmaras$.

On constate également que ces différences apparaissent à partir de l'avant du train, $Spalart- Allmaras$ majorant le saut de pression.

En l'absence de données expérimentales auxquelles comparer nos résultats, il n'est pas possible à partir de ces profils de pression de déterminer lequel des modèles est le plus proche de la réalité.

 

Intéressons nous maintenant aux profils de vitesse à l'arrière du train :

Le comportement général est encore une fois le même pour les trois modèles. Mais des différences importantes apparaissent principalement près de la paroi inférieure du tunnel. En effet le modèle $Spalart-Allmaras$ capture différemment les zones de recirculation contrairement aux modèles $k\omega  SST$ et $k-\epsilon$.

En zoomant près de la paroi, ce constat est flagrant (profils en $x=15 m$): 

Le modèle $k-\epsilon$ est peu performant, il capte moins bien la zone de recirculation. C'est un des défauts connus et attendus du modèle $k-\epsilon$, qui a tendance à lisser les courbes.

Pour le modèle $Spalart-Allmaras$, les vitesses en proche paroi ont un sens opposé à celui des deux autres modèles. En effet, comme on peut le voir avec le tracé ci-dessous des contraintes de cisaillement pariétal à l'aval du train, on voit que le modèle $Spalart-Allmaras$ est celui qui prédit le plus tôt le premier point de recollement. Ainsi, en $x=15 m $, on se situe avec le modèle $Spalart-Allmaras$  dans la zone entre les 2 recollements, d'où les vitesses de sens opposés, alors que pour les deux autres modèles, le premier point de recollement n'a pas encore été capté. Le fait que $Spalart-Allmaras$ soit conçu pour des écoulements externes nous laisse à penser que ce sont les résultats de ce modèle qui s'éloignent le plus de la réalité ici.

Au contraire, $k-\omega  SST$ capte très bien les recirculations en proche paroi. Une des raisons du choix de ce modèle et qu'il combine le modèle $k-\epsilon$ au coeur de l'écoulement et le modèle $k-\omega$ en proche paroi où il est très performant. Cela se retrouve dans nos résultats de profil de vitesse et on peut donc affirmer que parmi les trois modèles choisis $k-\omega  SST$ est le plus adapté pour traiter l'écoulement autour d'un train dans un tunnel.

 

Ci-dessous, nous avons tracé la contrainte de cisaillement pariétal derrière le train pour comparer l'abscisse des points de recollement suivant les modèles :

Comme évoqué ci-dessus, le modèle $Spalart-Allmaras$ est celui qui a le comportement le plus anecdotique : il prédit un peu plus tôt le premier point de recollement que les deux autres modèles et bien plus tard le second (d'environ $8 m$). De plus, au niveau du second point de recollement, l'évolution de la courbe est différente, plus lissée que pour les deux autres modèles.

Pour le modèle $k-\epsilon$, on observe également une évolution différente du profil par rapport aux deux autres modèles sur la partie $[5,6 ; 12] m $ où un pic apparaît, certainement produit par le comportement particulier de ce modèle en proche paroi. Ensuite, le modèle se rapproche plus du modèle $k-\omega  SST$ que de $Spalart-Allmaras$.

 

Conclusion

 
On retrouve le même comportement général pour les trois modèles de turbulence étudiés avec néanmoins quelques différences. Les valeurs de pression peuvent ainsi atteindre jusqu'à $7$% d'écart suivant le modèle utilisé. Mais la différence la plus importante se trouve au niveau des zones de recirculation à l'arrière du train que seul le modèle $k-\omega$SST parvient à capter correctement.
Cette étude tend donc à montrer que le modèle $k-\omega  SST$ est le plus performant pour modéliser la turbulence de l'écoulement dans un tunnel lors du passage d'un train.
Mais la validation de nos résultats n'est pas parfaitement rigoureuse, surtout en termes d'ordres de grandeurs des champs obtenus. Il n'y a en effet aucun article, aucunes données expérimentales qui correspondent au cas de notre simulation.
 

Étude paramétrique

Influence du vent dans le tunnel

 

L'air dans le tunnel n'est pas au repos et se déplace de la droite vers la gauche avec une vitesse comprise entre $15 m.s^{-1}$ et $23 m.s^{-1}$ d'après l'ingénieur de l'entreprise Setec.

Jusqu'ici nous avons considéré que cet écoulement d'air se déplaçait à $20 m.s^{-1}$. Pour étudier l'influence de ce déplacement d'air dans le tunnel, on va réaliser une étude paramètrique en comparant les résultats obtenus pour une vitesse de $15 m.s^{-1}$, $20 m.s^{-1}$ et $23 m.s^{-1}$.

La vitesse du train reste fixée à $40 m.s^{-1}$, la condition de vitesse en entrée du tunnel sera donc respectivement de $25 m.s^{-1}$, $20 m.s^{-1}$ et $17 m.s^{-1}$.

En utilisant le modèle le plus performant $k-\omega  SST$ et en conservant une intensité de turbulence de $1$%, nous avons déjà réalisé les simulations pour une vitesse de $20 m.s^{-1}$ qui est notre condition d'entrée de base.

Pour lancer les simulations pour des vitesses entrantes de $U_\infty = 25 m.s^{-1}$ et $U_\infty = 17 m.s^{-1}$, il faut recalculer les paramètres de turbulence correspondants : 

 

pour $U_\infty = 25 m.s^{-1}$  $\boxed{k = 0,094 m^{2}.s^{-2}  et   \omega = 80 s^{-1}}$

 

pour $U_\infty = 17 m.s^{-1}$ : $\boxed{k = 0,043 m^{2}.s^{-2}  et  \omega = 54,1 s^{-1}}$

 

Les deux simulations ont bien été effectuées, mais malheureusement, le temps qui nous était imparti a été insuffisant pour exploiter les résultats.

Simulation avec maillage mobile

L'étape finale de ce projet et de mettre en commun notre étude sur les modèles de turbulence avec le travail sur le maillage mobile réalisé par le groupe 23.

 

Caractéristiques du maillage mobile

 

L'équipe maillage mobile n'a pas utilisé le même maillage que nous, mais à mailler son propre domaine de calcul comme suivant :

- une zone à l'arrière du train où les mailles sont étirées,

- la zone du train où les mailles ne se déforme pas, cette zone est juste translatée,

- une zone à l'avant où les mailles sont compressées.

Ils ont ensuite défini un critère de compression des mailles qui impose de remailler les domaines aval et amont avant de reprendre les calculs. Toutes les tâches ont été automatisées à l'aide d'un script en bash.

On ne peut plus imposer les mêmes conditions aux limites qu'en maillage fixe. En maillage mobile, on n'est plus dans le cas d'un écoulement en soufflerie, la condition en vitesse imposée est celle du train et non celle d'un écoulement.

La vitesse du train est donc fixé à $40m.s^{-1}$, le train se déplaçant de la gauche vers la droite, maintenant. De plus, l'air n'étant pas au repos dans le tunnel, un écoulement de $20m.s^{-1}$ dans le même sens que le déplacement du train est initialisé à l'aval de celui-ci, afin de se rapprocher de la réalité du phénomène.

Plus de renseignements sont disponible ici.

 

Mise en place de la turbulence

 

Ce qui ressort de notre étude est que le modèle $k - \omega  SST$ donne les meilleurs résultats. De plus c'est un modèle complet et qui est capable de résoudre la sous-couche visqueuse pour des maillages raffinés comme le maillage crée par le groupe $23$, plus raffiné ($y_{min}^{+} = 2,5$) que celui utilisé pour notre étude. 

Pour la partie nous concernant, c'est à dire les fichiers $k$ et $\omega$, nous avons dû modifier les conditions aux limites afin de les adapter aux noms des patchs choisis par le groupe Maillage mobile.

De plus, lors du cas avec le maillage fixe, nous avions parlé de la création automatique d'un fichier nut lors de la simulation. Ici, il est nécessaire de paramètrer ce fichier. En effet, le $y^{+}$ étant suffisamment faible aux parois du train, il convient de placer la première maille dans la sous-couche visqueuse afin de résoudre correctement la couche limite. Pour cela on applique aux parois du train la condition limite suivante : nutLowReWallFunction, afin de prendre en compte l'approche Low Reynolds aux parois. La valeur de $\nu_t$  aux parois sera fixée à 0.

Pour les valeurs du modèle $k-\omega  SST$, nous avons choisi une intensité turbulente de $5$% soit $k = 1,5 m^{2}.s^{-2}$ et $\omega = 319,4 s^{-1}$. La valeur de nut en entrée et en sortie sera calculé automatiquement grâce à la fontion calculated.

Voici, par exemple, le fichier nut utilisé lors de la simulation avec le maillage mobile.

 

Résultats

Voici l'établissement du champ de vitesse à l'arrière du train avec une caméra fixe au niveau du train mobile :

Ou à un instant donné :

On remarque cependant qu'à l'arrière du train le maillage n'est pas assez raffiné en proche paroi, le comportement de la couche limite est donc peu fiable.

Pour le champ de pression, on retrouve quelque chose de comparable au cas du maillage fixe. Cependant, la gamme de valeur n'est pas la même que précédemment,  le régime stationnaire n'étant pas atteint ici. Nous n'avons en effet pas pu simuler suffisamment de temps physique pour l'atteindre ( limitation par la longueur finie du domaine et temps de calcul long).

Du fait de certains problèmes survenus lors des calculs (arrêt inexpliqué en plein milieu de la simulation...) et du temps important pour obtenir un résultat sur les machines à notre disposition, nous n'avons pas pu exploiter entièrement le cas du maillage mobile.
 

Conclusion

Conclusion

 

Le comportement en pression obtenu coïncide avec la littérature. Nos simulations permettent donc d'obtenir des résultats vérifiant la physique du problème. Malgré tout, l'absence de données expérimentales se rapprochant des conditions de notre étude ne nous a pas permis de vérifier l'exactitude des ordres de grandeurs de nos résultats.

Ce projet a permis de déterminer que $k-\omega  SST$ est le modèle le plus adapté pour simuler les effets de la turbulence lors du passage d'un train dans un tunnel.

Nous avons également mené une étude comparative des résultats obtenus avec les logiciels OpenFoam et Fluent. Cela a permis de montrer qu'OpenFoam obtient des résultats très proches de ceux de Fluent et qu'il est donc tout à fait possible d'utiliser ce logiciel OpenSource pour ce type de simulation. La license Fluent étant onéreuse, cette conclusion est très intéressante d'un point de vue économique.

Encore une fois, l'absence de données expérimentales ne nous a pas permis de trancher sur lequel des deux logiciels est le plus proche de la réalité. Il pourrait donc être intéressant de comparer avec un troisième code de calcul (StarCCM+ par exemple).

En ce qui concerne le carctère post-dictif de la turbulence avec l'approche RANS, un moyen d'éviter cela aurait été de mettre en place ce calcul avec les méthodes LES (Large Eddy Simulation). Cette approche est disponible sous OpenFoam mais est certainement plus compliqué d'utilisation et nécessite des moyens de calculs plus importants que ceux disponibles à l'ENSEEIHT.

Ce sujet d'étude est très intéressant et reste ouvert sur de nombreux points. Pour aller plus loin, il serait ainsi intéressant de simuler le problème en trois dimensions pour se rapprocher encore plus du cas  réel. Dans la même optique, la prise en compte de la compressibilité de l'air (et donc simuler un train à grande vitesse) pourrait faire l'objet d'une prochaine étude.

 

Bibliographie / Webographie

Bibliographie

 

[B1]  Eric Furbo, Janne Harju and Henric Nilsson -  Evaluation of Turbulence Models for Prediction of Flow Separation at a Smooth Surface - Uppsala University - June 2009

[B2] Luke Jongebloed - Numerical Study using Fluent of the Separation and the Reattachment Points for Backwards-Facing Step - Rensselaer Polytechnic Institute, Hartford, Connecticut - December 2008

[B3] Rodler J., Hagenah B. - Determination of Aerodynamic Burden in Rail Tunnels Using Measurements and Simulation -  Gruner GmH Ingenieure und Planer, Austria - $6^{th}$ International Conference "Tunnel Safety and Ventilation" 2012, Graz

[B4] Jakub Novák - Single Train Passing through a Tunnel - Skoda Research, Dep. of Fluid Mechanics - ECOMACS CFD 2006

[B5] Rudol Bopp, Bernd Hagenah - Aerodynamics, Ventialtion and Tunnel Safety for High-Speed Rail Tunnels - Gruner

[B6] Patrick Haas - Aérodynamique des trains rapides, Exposés des phénomènes relevant de la circulation en tunnel - CMEFE - Mars 2000

 

 

Webographie

 

[W1] http://www.openfoam.com/

[W2] http://www.cfd-online.com/

[W3] http://hmf.enseeiht.fr/travaux/bei/beiep/content/g21/5-les-logiciels-cfd-utilises-openfoam-et-starccm