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.