Ajout de la turbulence

Pour obtenir une simulation réaliste, il est évidemment nécessaire de prendre en compte la turbulence dans le calcul, puisque le $Re$ est de $20 \cdot 10⁶$, ce qui correspond à un régime pleinement turbulent.

Nous nous sommes entièrement basés sur les conclusions du travail du groupe n°24 (Voir ici). Selon eux, deux modèles étaient à privilégier, le modèle $k-\epsilon$ et le modèle $k-\omega$. Nous nous attacherons ici à détailler brièvement ces deux modèles et les équations qui les composent. Les valeurs utilisées ainsi que les conditions limites appliquées pour les nouveaux paramètres de la turbulence seront explicités.

Rappel sur les modèles utilisés


Les modèles de turbulence sont basés sur des équations de transport. Il existe des modèles à une équation, à deux équations ($k-\epsilon$ et $k-\omega$ par exemple) et plus encore (Reynolds Stress à 7 équations etc...). Les modèles à deux équations sont basés sur l'hypothèse de Boussinesq qui suppose que le tenseur de Reynolds qui apparaît lorsque l'on réécrit les équations de Navier-Stokes en utilisant une décomposition de Reynolds des grandeurs et en les moyennant, est proportionnel au tenseur des vitesses de déformations.

Le modèle $k-\epsilon$


Le modèle $k-\epsilon$ est le modèle le plus commun et le plus utilisé en CFD, il est un peu le "couteau-suisse" des modèles de turbulence. C'est un modèle à deux équations de transport, une pour l'énergie cinétique turbulente $k$ et une pour le taux de dissipation $\epsilon$.

Dans ce modèle, la viscosité turbulente s'écrit ainsi :

$\nu_t=\rho C_{\mu} \frac{k²}{\epsilon}$

Ces deux équations sont les suivantes :

$\frac{\partial k}{\partial t}+ U_k \frac {\partial k}{\partial x_k}=\frac {\partial}{\partial x_k}\left[ \left( \nu+ \frac {C_{\mu} k²}{\sigma_k \epsilon}\right)\frac {\partial k}{\partial x_k}\right]+\frac {C_{\mu} k²}{\sigma_k \epsilon}\left(\frac {\partial U_i}{\partial x_k}  +\frac {\partial U_k}{\partial x_i}\right)\frac {\partial U_i}{\partial x_k}  -\epsilon$

$\frac{\partial \epsilon}{\partial t}+ U_k \frac {\partial \epsilon}{\partial x_k}=\frac {\partial}{\partial x_k}\left[ \left( \nu+ \frac {C_{\mu} k²}{\sigma_{\epsilon} \epsilon}\right)\frac {\partial \epsilon}{\partial x_k}\right]+C_{\epsilon_1}C_{\mu}k\left(\frac {\partial U_i}{\partial x_k}  +\frac {\partial U_k}{\partial x_i}\right)\frac {\partial U_i}{\partial x_k}  -C_{\epsilon_2} \frac{\epsilon²}{k}$

Les paramètres prennent les valeurs suivantes :

$C_{\mu}$ $\sigma_k$ $\sigma_{\epsilon}$ $C_{\epsilon_1}$ $C_{\epsilon_2}$
0.09 1.0 1.22 1.44 1.92

Ce modèle est connu pour avoir des comportements pathologiques dans certains types d'écoulements, notamment ceux où de forts gradients de pression adverses sont présents. Comme tous les modèles de type RANS, la qualité de la modélisation dépend aussi fortement des valeurs de l'énergie cinétique turbulente et du taux de dissipation injectées dans l'écoulement. En effet, il est incapable de générer de la turbulence, et par conséquent de mauvaises valeurs de $k$ ou d'$\epsilon$ vont soit dissiper trop de turbulence, soit en générer beaucoup trop. Il faut donc être très prudent avec les résultats et les paramètres utilisés.

Le modèle  $k-\omega$ SST


Le modèle $k-\omega$ SST est également un modèle à deux équations. Il s'agit toujours de deux équations de transport, une sur $k$, l'énergie cinétique turbulente, et une autre sur $\omega$ qui peut être vu comme une fréquence caractéristique de la turbulence. Il s'écrit aussi $\omega=\frac{\epsilon}{k}$.

La viscosité cinématique turbulente s'écrit :

$\nu _T = \frac{a_1 k}{max(a_1 \omega, S F_2)}  $

Et les deux équations de transport sont les suivantes :

$\frac{\partial k}{\partial t} + U_j \frac{\partial k}{\partial x_j } = P_k - \beta ^* k\omega + \frac{\partial}{\partial x_j } \left[ {\left( {\nu + \sigma_k \nu _T } \right){{\partial k} \over {\partial x_j }}} \right]$

${{\partial \omega } \over {\partial t}} + U_j {{\partial \omega } \over {\partial x_j }} = \alpha S^2 - \beta \omega ^2 + {\partial \over {\partial x_j }}\left[ {\left( {\nu + \sigma_{\omega} \nu _T } \right){{\partial \omega } \over {\partial x_j }}} \right] + 2( 1 - F_1 ) \sigma_{\omega 2} {1 \over \omega} {{\partial k } \over {\partial x_i}} {{\partial \omega } \over {\partial x_i}}$

Les différents termes introduits dans ces deux équations sont les suivants :

$F_2=\mbox{tanh} \left[ \left[ \mbox{max} \left( { 2 \sqrt{k} \over \beta^* \omega y } , { 500 \nu \over y^2 \omega } \right) \right]^2 \right]$

$P_k=\mbox{min} \left(\tau _{ij} {{\partial U_i } \over {\partial x_j }} , 10\beta^* k \omega \right)$

$F_1=\mbox{tanh} \left\{ \left\{ \mbox{min} \left[ \mbox{max} \left( {\sqrt{k} \over \beta ^* \omega y}, {500 \nu \over y^2 \omega} \right) , {4 \sigma_{\omega 2} k \over CD_{k\omega} y^2} \right] \right\} ^4 \right\}$

$CD_{k\omega}=\mbox{max} \left( 2\rho\sigma_{\omega 2} {1 \over \omega} {{\partial k} \over {\partial x_i}} {{\partial \omega} \over {\partial x_i}}, 10 ^{-10} \right )$

$\phi = \phi_1 F_1 + \phi_2 (1 - F_1)$

Les constantes du modèle sont définies dans le tableau suivant :

$\alpha_1$ $\alpha_2$ $\beta_1$ $\beta_2$ $\beta^*$ $\sigma_{k1}$ $\sigma_{k2}$ $\sigma_{\omega_1}$ $\sigma_{\omega_2}$
$0.55$ $0.44$ $0.075$ $0.0828$ $0.09$ $0.85$ $1$ $0.5$ 0.856

Ce modèle a été développé par Menter et cherche à bénéficier de l'avantage du modèle $k-\omega$ en proche-paroi et du modèle $k-\epsilon$ dans la partie principale de l'écoulement, évitant ainsi les comportements indésirables des deux modèles dans l'écoulement.


Paramètres de la simulation

La simulation Fluent réalisée par SETEC a été lancée avec une intensité turbulente de $1\%$. A partir de là, on peut en déduire les valeurs pour $k$ et $\epsilon$ pour obtenir une intensité équivalente :

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

$\epsilon=\frac{C_{\mu}^{\frac{3}{4}}k^{\frac{3}{2}}}{l}$

 

On obtient alors les valeurs suivantes pour $k$ et $\epsilon$ : $k=0.06~m^2.s^{-2}$ et $\epsilon=0.35~m^2.s^{-3}$.

Ces valeurs permettent d'initialiser les champs de $k$ et d'$\epsilon$ au début de la simulation. La vitesse initiale dans le domaine est initialisée à $20~m/s$ afin de simuler un courant d'air dans le tunnel.

Pour la simulation $k-\omega$ SST, on réalise les mêmes calculs, et on trouve $k=0.06~m^2.s^{-2}$ et $\omega=64.8s^{-1}$.

Conditions limites

On rappelle que l'intensité turbulente vaut $1\%$.

La figure suivante répertorie les différentes conditions limites utilisées dans cette simulation.

Pour la simulation avec $k-\epsilon$, les conditions limites sont les suivantes :

Champs Right Left farFieldLeft farFieldRight farFieldMoving Hull
k fixedValue 0.06 fixedValue 0.06 kqRWallFunction kqRWallFunction kqRWallFunction kqRWallFunction
$\epsilon$ fixedValue 0.35 fixedValue 0.35 epsilonWallFunction epsilonWallFunction epsilonWallFunction epsilonWallFunction
$\nu_t$ calculated calculated nutkWallFunction nutkWallFunction nutkWallFunction nutkWallFunction

Pour la simulation avec $k-\omega$ SST :

Champs Right Left farFieldLeft farFieldRight farFieldMoving Hull
k fixedValue 0.06 fixedValue 0.06 kqRWallFunction kqRWallFunction kqRWallFunction kqRWallFunction
$\omega$ fixedValue 64.8 fixedValue 64.8 omegaWallFunction omegaWallFunction omegaWallFunction omegaWallFunction
$\nu_t$

calculated

calculated nutLowReWallF* nutLowReWallF* nutLowReWallF*

nutLowReWallF*

nutLowReWallF* : nutLowReWallFunction

Notre $y+$ proche des parois étant toujours  élevé, des modèles de loi de paroi doivent être utilisés, d'où les conditions limites pour $k$,$\omega$,$\epsilon$ et $\nu_t$ aux parois.


Finalement, une fois la turbulence implantée, et le script prêt, il ne reste plus qu'à lancer les simulations finales.