Simulation numérique

Présentation du code NSMB

Le code Navier-Stokes Multi-Block est un solveur Navier Stokes parallèle développé par l'EPFL (Ecole Polytechnique Fédérale de Louisiane), le KTH (Kungl. Tekniska Hoegskolan), le CERFACS et Airbus France. L'Institut de Mécanique des Fluides de Toulouse participe à son développement depuis 2001.

L'avantage principal de ce solveur est la parallélisation du calcul par la décomposition automatique du domaine en différents blocs. Il est compatible avec des maillages structurés très variés pour des configurations allant des écoulements subsoniques incompressibles jusqu'aux écoulements supersoniques.

Codé en Fortran 77, il utilise la méthode des volumes finis avec des schémas spatiaux centrés ou amont et une intégration temporelle explicite ou implicite. 

Il existe une grande diversité de modèles implémentés permettant des calculs RANS (Reynolds Average Navier-Stokes) ainsi que des calculs LES (Large Eddy Simulation).

 

Modèle Splalart Almaras

Les simulations mises en place sont des calculs RANS utilisant le modèle de Spalart-Almaras avec la correction d'Edward et Chandra. Ce dernier constitue un modèle à une équation de transport de la viscosité turbulente permettant la fermeture du système pour la résolution des équations de Navier-Stokes moyennées définies comme suit:

  • Equation de conservation de la masse:

\begin{eqnarray}\frac{\partial \bar{\rho} }{\partial t}+\frac{\partial \bar{\rho} \tilde{u_{j}}}{\partial x_{j}}=0\end{eqnarray}

  • Equation de quantité de mouvement

\begin{eqnarray}\frac{\partial \bar{\rho} \tilde{u_{i}}}{\partial t}+\frac{\partial \bar{\rho}  \tilde{u_{i}} \tilde{u_{j}}}{\partial x_{j}}+\frac{\partial \bar{p}}{\partial x_{j}}-\frac{\partial \bar{\tau} _{ij}-\bar{\rho}\widetilde{u''_{i}u''_{j}}}{\partial x_{j}}= 0\end{eqnarray} Avec le tenseur de Reynolds \begin{eqnarray}\tau_{r} = \bar{\rho} \widetilde{u''_{i}u''_{j}}\end{eqnarray}

  • Equation de conservation de l'énergie

\begin{eqnarray}\frac{\partial \bar{\rho} \hat{E}}{\partial t}+\frac{\partial \bar{\rho} \tilde{u_{j}}  \hat{H}}{\partial x_{j}}+\frac{\partial \bar{q_{j}}+\bar{\rho}\widetilde{e''u''_{j}}}{\partial x_{j}}-\frac{\partial (\bar{\tau} _{ij}-\bar{\rho}\widetilde{u''_{i}u''_{j}})\tilde{u}_{i}}{\partial x_{j}}=-T_{s}^{k}\end{eqnarray}

Avec

L'énergie totale: \begin{eqnarray}\hat{E}=\tilde{e}+\frac{1}{2}\tilde{u}_{k}\tilde{u}_{k}\end{eqnarray}

L'enthalpie totale: \begin{eqnarray}\hat{H}=\tilde{e}+\frac{\bar{p}}{\rho}\end{eqnarray}

Tsk un terme source d'énergie cinétique turbulente

L'hypothèse de Boussinesq donne \begin{eqnarray}-\widetilde{u''_{i}u''_{j}}=2 \nu_{t}S_{ij}\end{eqnarray}

Le modèle implémenté donne donc l'évolution de la viscosité turbulente.

Théorie du modèle

L'équation de transport de la viscosité turbulente est implémentée comme suit:

\begin{eqnarray}\frac{D \tilde{\nu}}{Dt}=c_{b1}[1-f_{t2}]\tilde{S}\tilde{\nu} +\frac{1}{\sigma}[\nabla ((\nu +\tilde{\nu})\nabla\tilde{\nu})+c_{b2}(\nabla\tilde{\nu})^{2}]-(c_{w1}f_{w}-\frac{c_{b1}}{\kappa}f_{t2})(\frac{\tilde{\nu}}{d})+f_{t1}\Delta U^{2}\end{eqnarray}

 

 

 

 

 

 

Avantages

 

 

 

 

 

Limites:

Mise en place des simulations

Les simulations ont été réalisées à partir d'un profil 2D d'aile super critique spécifique fournit par l'industriel Dassault. La configuration à partir de laquelle a été menée les simulations est un profil de corde 0.25m et le domaine s'étend sur huit fois la longueur de la corde. L'origine du système est prise au niveau du bord d'attaque de l'aile.

Profil V2C

Le maillage préalablement réalisé est de type structuré en forme C et raffiné en proche paroi afin de capter les échelles de la couche limite. Au total, le domaine a été discrétisé sur 160 000 mailles.

Maillage du profil V2C

Différentes configurations ont été définies:

Configurations

Transition (% de corde)

Mach

Angle d'attaque (°)

 
1 0 0.7 0 à 7
2 0 0.75 0 à 7
3 30 0.7 0 à 7
4 30 0.75 0 à 7

Différentes configurations étudiées

Le calcul est ensuite réalisé à partir d'un schéma spatial amont avec une intégration temporelle de type implicite. La résolution se fait par une méthode RANS en stationnaire et URANS en instationnaire à partir du modèle de Spalart-Almaras présenté dans la paragraphe précédent.

 

Démarche de calcul

 La démarche utilisée pour chaque configuration suit la logique suivante:

La vérification de la convergence du calcul stationnaire permet de repérer l'angle d'attaque critique à partir duquel le phénomène de tremblement apparaît. Elle se fait par le tracé du résidu du calcul au cours des itérations.

Cela permettra ensuite d'analyser l'influence des paramètres tels que le Mach, et le type de couche limite sur le déclenchement du buffet.

Si le buffet apparaît, le calcul instationnaire est réalisé. Pour ce type de simulations, il est important de garantir la qualité des résultats par une analyse de la mise en place du régime permanent. En effet, le régime transitoire n'a aucune valeur d'un point de vue physique et il est nécessaire de vérifier  la convergence vers un mode cyclique d'amplitude et de fréquence constante. Ceci est caractéristique de l'instationnarité du tremblement.

 

 

 

 

 

 

 

 

 

L'utilitaire forces2tec est utilisé pour calculer à partir des fichiers de sortie de NSMB, la courbe d'évolution temporelle des coefficients aérodynamiques adimensionnalisés (Cx,Cz,Cm). Cela permet d'observer l'atteinte du régime cyclique et de mesurer la fréquence du buffet. 

Ex. configuration 1: Couche limite turbulente, Mach 0.7, Alpha 5°

Le post-traitement complet se fait ensuite via Tecplot. Il s'agit d'un outil de visualisation de calculs CFD à partir duquel les champs et profils de multiples grandeurs peuvent être obtenues. L'utilitaire n2t permet de traduire les fichiers de sortie de NSMB en format lisible par Tecplot. C'est de cette façon que les champs de vitesse, de pression et le profil du coefficient de portance autour de l'aile ont été tracés.