Projets des groupes         Autres travaux               Imft                     EDF




 

Présentation des différentes méthodes numériques utilisées

Bien que notre problème fasse intervenir un écoulement incompressible, nous le résolvons à l'aide du code NSMB, conçu pour modéliser des écoulements compressibles. Nous allons donc dans la suite détailler différentes méthodes numériques utilisées dans le code de calcul pour résoudre les équations de Navier Stokes dans le cas d'écoulements compressibles.

Equations de Navier-Stokes pour les écoulements compressibles


Equation de conservation de la masse :



Equation de conservation de la quantité de mouvement :



Equation de conservation de l'énergie totale :



Nous pouvons écrire ces équations sous la forme d'un système vectoriel de la forme
suivante :



U est le vecteur d'état des variables conservatives, Fi est le flux advectif et Fivis est le flux diffusif. Le nombre de Reynolds utilisé dans ce système est défini comme tel : Re0=(ρ0
L0V0)/μ0
Ces équations sont ensuite discrétisées en temps et en espace. La discrétisation spatiale de ces équations peut se faire de plusieurs manières différentes. Le choix du schéma numérique dépend de la nature de l'écoulement à modéliser. Dans ce projet, les termes diffusifs des équations de Navier-Stokes seront toujours discrétisés avec un schéma centré. Nous décrirons dans les paragraphes suivants les schémas de discrétisation centré et amont.

Schéma amont


L'idée générale de cette approche par rapport au schéma centré est la prise en compte des propriétés physiques de l'écoulement, en considérant notamment les directions privilégiées de la propagation de l'information dans l'écoulement. Dans la pratique, ce type de schéma présente l'avantage d'ajouter une dissipation qui permet d'éviter les oscillations au voisinage des discontinuités, comme cela peut être le cas pour le schéma centré.
Si on considère un point de coordonnées (a,b,c), les termes convectifs ∂Fi/
∂xi s'écrivent sous la forme discrétisée suivante :

avec


Où ∆xi représente le pas de discrétisation dans la direction i. Aabc représentent les matrices de Roe, dont le détail est donné par Bouhadji. Les expressions des flux numériques aux interfaces (a-1/2,b,c), (a,b-1/2,c ), (a,b,c-1/2 ) se déduisent des expressions précédantes.

Schéma centré


Comme pour le schéma amont présenté dans la section précédante, si on considère un point de coordonnées (a,b,c), les termes convectifs
∂Fi/∂xi s'écrivent sous la forme discrétisée suivante :


avec


Par contre, contrairement au schéma amont, la discrétisation des termes suivants de
fait de manière centrée :


 

Présentation des différentes approches de la turbulence

Les approches LES


Dans les simulation aux grandes échelles (Large Eddy Simulation), seuls les processus de grande taille sont simulés alors que les effets des petites structures sont modélisés par une loi de "sous-maille". On néglige donc la partie du spectre qui correspond aux petites échelles. Celles-ci obéissent au hypothèses d'équilibre de la turbulence homogène et isotrope, et peuvent être modélisées par un terme assurant la dissipation de l'énergie provenant des structures résolues de plus grande taille.
La mise en place de la simulation aux grandes échelles implique la mise en place d'un filtre spatial permettant de distinguer les structures résolues de celles dont les effets devront être modélisés. Soit v une variable dépendant de x et de t, sa valeur filtrée s'exprime de la manière suivante :


Où G(x, x', ∆) représente l'opérateur de filtrage spatial au point x et ∆ est un paramètre
caractéristique de la coupure spectrale de ce filtre.
L'ensemble des variables de l'écoulement peuvent être décomposées de la manière suivante:


où v' désigne la fluctuation de v par rapport à la moyenne filtrée.
Les approches LES nécessitent un maillage très fin dans les zones proche paroi pour lesécoulements à haut nombre de Reynolds. Dans ce projet, nous n'utiliserons donc pas cette approche. Cette succinte présentation permet simplement de situer les différentes approches les unes par rapport aux autres.

Les approches de type URANS


Ces deux approches utilisent un traitement statistique de la turbulence à priori menant à une résolution des équations de Navier-Stokes moyennées :

     
où Ũi représente une grandeur moyennée et ui une grandeur fluctuante.
Les différentes approches diffèrent alors selon la séparation des variables effectuée. La séparation des variables en 3 composantes a été proposée par Hussain et Reynolds :

                                                                                            
où 
Ūi représente la composante moyenne de l'écoulement, Űi représente la fluctuation périodique et ũi la fluctuation aléatoire. Les approches de modélisation se distingue par l'utilisation d'une moyenne d'ensemble (URANS) ou d'une moyenne de phase (OES).


  • Les approches URANS classiques :
La présence d'instationnarités et de structures organisées dans les écoulements turbulents conduit à ne plus considérer l'ensemble du mouvement fluctuant comme aléatoire, et ainsi à adopter des approches instationnaires. L'approche URANS, Unsteady Reynolds Average Navier-Stokes, consiste en une séparation de la moyenne statistiques de l'ensemble des fluctuations, chaotiques comme organisées telle que: Ui = Ūi + ui' où ui' regroupe les fluctuations périodiques et aléatoires.
Les équations du mouvement moyen sont les mêmes que les équations RANS, au terme temporel près. Bien qu'elle soit robuste dans les zones proche parois pour les écoulements à grand nombre de Reynolds, cette approche implique souvent une forte dissipation et conduit à de faibles prédictions en ce qui concerne les aspects de non-équilibre statistique induit par les instationnarités. Ceci s'explique par le manque d'adaptation des échelles de vitesse et de longueur du mouvement turbulent utilisées dans ces modèles.

  • Les approches OES :

La méthode de macrosimulation OES (Organised Eddy Simulation) est basée sur la résolution des équations de Navier Stokes en moyenne de phase, c'est à dire la partie cohérente des fluctuations n'est plus modélisée mais simulée menant à une décomposition en moyenne de phase :
Ui = <Ui> + ui'  où <Ui> regroupe la moyenne d'ensemble Ũi et les fluctuations cohérentes Űi .
Cette approche a été proposée afin de modéliser différemment le terme des contraintes turbulentes. Ceci est nécessaire car du point de vue spectral l'interaction non-linéaire entre les structures cohérentes et la turbulence aléatoire induit une modification de la pente de la partie continue du spectre des fluctuations turbulentes par rapport à la pente (-5/3) observée dans les écoulements turbulents pleinement développés en équilibre, en accord avec théorie statistique de Kolmogorov dans la zone inertielle.

Il faut donc reconsidérer les échelles de temps et de longueur de la turbulence pour cette
approche et reconsidérer par conséquent les modèles de turbulence à deux équations à l'aide notamment de la modélisation au second ordre afin de bien prédire l'instationnarité du décollement et des structures cohérentes dans le sillage. Ainsi, un nouveau coefficient de diffusivité a été déterminé : il a été montré que les valeurs de Cμ étaient plus faibles (de l'ordre de 0,02) que les valeurs utilisées pas les modèles à deux équations (Cμ =0,09).
La présence de structures cohérentes implique la présence d'un ou plusieurs pics sur les spectres de quantités fluctuantes. L'approche OES isole ces pics du spectre afin de calculer le mouvement cohérent leur correspondant, et ainsi de modéliser la partie continue du spectre qui correspond à tout le mouvement aléatoire, des basses jusqu'aux hautes fréquences.

Du fait de l'étendue du spectre à modéliser dans toute la gamme des nombres d'onde,
elle permet d'utiliser des concepts de la modélisation statistique, moyennant une forme de spectre modifié par rapport au spectre en équilibre statistique. On parle donc de l'élaboration de modèles statistiques avancés.
Ainsi cette modélisation permet une bonne prédiction du développement des structures organisées ainsi que de l'amplification des instabilités et des fréquences des oscillations apparaissant sous l'effet de la viscosité.


La figure précédente montre que la distinction entre les
structures résolues (1) et celles dont le comportement doit être modélisé (2). Cette distinction est basée sur le comportement organisé ou aléatoire des différentes structures.

Les approches hybrides


Les approches hybrides permettent de concillier les avantages des deux approches présentées ci-dessus. En effet, la modélisation statistique de la turbulence est utilisée en région proche paroi et la modélisation de type LES est utilisée dans les zones de d'écoulement détaché. Ceci est mis en oeuvre par le même système d'équations que la méthode URANS, en choisissant l'échelle de longueur de la turbulence dans chaque volume élémentaire de fluide selon la relation suivante : l=min(lURANS , lLES ). Cette échelle de longueur apparaît dans le terme de dissipation de l'équation de transport de l'énergie cinétique turbulente.



Présentation des différents modèles de turbulence utilisés

Les modèles de fermeture du premier ordre sont basés sur l'hypothèse de Boussinesq qui relie linéairement le tenseur de Reynolds au champ moyen des vitesses, formulée par Prandtl sous la forme suivante :

avec νt la viscosité turbulente, νt = Cμ κ² / ε .

Les modèles à deux équations reposent sur le
transport de l'énergie cinétique turbulente κ et celui de la dissipation de la turbulence ou d'une grandeur contenant ce terme.

Lorsque l'écoulement est de type pariétal, il convient d'atténuer la turbulence (damp
ing) dans la zone de proche paroi. Les approches d'amortissement sont basées sur des développements asymptotiques des grandeurs turbulentes dans cette zone.

Nous présentons dans la suite les quatres modèles utilisés au cours de notre étude :

- le modèle κ − ε Chien.
    - les modèles κ − ε Chien/OES et κ − ω /OES.
    - le modèle κ − ω Baseline.
    - le modèle κ − ω SST.


Modèle k-ε Chien


Le modèle κ−
ε est un modèle de fermeture du premier ordre de turbulence et fait intervenir deux équations de transport : une équation pour la dissipation de l'énergie cinétique turbulente ε, et une équation pour l'énergie cinétique turbulente κ, dont la structure est la suivante :


La viscosité turbulente est donnée par :

Pour fermer complètement ce système il reste à déterminer les constantes CμC1C2, σ. Ces constantes sont déterminées à partir de résultats expérimentaux, et pour nos calculs nous prendrons Cμ = 0.09, C1 =1.35, C2 =1.8, σ =1.3.

Le modèle de Chien utilise la méthode des développements limités en série de Taylor pour étudier le comportement des contraintes turbulentes, de l'énergie cinétique et du taux de dissipation près de la paroi.
On obtient une nouvelle forme pour κ et ε respectivement énergie cinétique turbulente et dissipation isotrope (pseudo-dissipation) :


ll est à noter que ce n'est pas le taux de dissipation réel :


En raison de la condition de non-glissement à la paroi et de l'équation de continuité pour les fluides incompressibles, on obtient :


Il faut donc ajouter un terme additionnel représentant le taux de dissipation réel à la paroi afin de compenser le terme de diffusion moléculaire de κ.
On obtient à la paroi :


avec ν la viscosité cinématique.
Nous remarquons alors que D est égal à la diffusion
moléculaire pour y=0.
D correspond alors à la dissipation pariétale. Le terme de
dissipation dans l'équation de transport de l'énergie cinématique sera alors noté ε+( 2νκ / y² )
On écrit alors la forme bas Reynolds de base de κ :


L'équation de transport de la dissipation du modèle bas Reynolds de Jones et Launder revient à prendre N = 0 dans l'approximation d'ε. Le terme de dissipation dans l'équation de
transport de ε à grand nombre de Reynolds est alors fini pour y = 0 et il équilibre
le terme de la diffusion moléculaire pariétale. De plus la forme limite de νt donnée par
les équations précedentes est automatiquement satisfaite. Il semble toutefois nécessaire
d'ajouter aux équations de transport des termes empiriques afinn de notamment prendre
en compte l'effet d'atténuation dû à la présence d'une paroi solide.
On obtient alors le modèle suivant :

Modèles k-ε / OES et k-ω /OES


Lors d'une modélisation OES, l'écoulement moyen et les fluctuations organisées sont simulées alors que les fluctuations chaotiques sont modélisées. Ceci a pour effet de rompre l'égalité entre la production et la dissipation d'énergie cinétique. Il faut donc réévaluer certains paramètres.
Après calculs, il apparaît qu'il faut utiliser un Cμ équiva
lent ayant pour valeur Cμ=0.02 au lieu de 0.09 comme nous l'avons vu dans les modèles U-RANS classiques.

Modèle k-ω baseline


Il s'agit d'un modèle du premier ordre à deux équations : une équation de transport pour l'énergie cinétique turbulente k, et une équation de transport pour une fréquence caractéristique de la turbulence ω .
Ce modèle dérive du modèle k-ω de Wilcox. Cependant ce dernier a un tres grande sensibilité aux conditions "free stream". L'idée a donc été de coupler ce modèle au modèle k-
ε pour obtenir le modèle k-ω baseline. La transformation consiste essentiellement en un ajout d'un terme de diffusion et une modification des constantes des modèles.

Équation de k :


Équation de ω :


où la fonction F1 est définie telle que F 1 = tanh(arg14 )

Les différentes constantes de ce modèle ont été obtenues de la même façon que pour le modèle k-ε.

Modèle k-ω SST (Shear Stress Transport)


Ce modèle combine les avantages de deux modèles : d'une part, il permet d'avoir de bons résultats près des parois solides, et il peut être utilisé sans fonction de "damping", et d'autre part il limite le surplus de viscosité. En effet, pour remédier au surplus de viscosité généré par les approches URANS classiques, Menter a proposé d'introduire des limiteurs de viscosité turbulente en présence d'un gradient de pression adverse.

Modèle Spalart-Allmaras


Le modèle Spalart Allmaras est un modèle du premier ordre à une équation pour la viscosité turbulente. Dans ce type de modèle, la viscosité de turbulence s'obtient par une fonction auxiliaire, solution d'une équation de transport spécifique. La viscosité turbulente est donnée par νtν~ .fν
avec le terme d'amortissement visqueux fν = χ3 / (Cν33) où χ = ν~/ν
L'équation de transport de ν~ est donnée par :

Avec les fonctions auxiliaires suivantes :



Le modèle de Spalart-Allmaras est considéré comme un bon compromis entre les modèles
algébriques et les modèles à deux équations de transport.

Interaction fluide-structure

Nous avions initialement prévu d'effectuer des simulations numériques pour décrire le phénomène d'instabilité fluide élastique, cependant les résultats obtenus dans les études préalables n'étaient pas suffisament satisfaisants pour mener une étude dynamique. La bonne prédiction, dans le cas statique, des instationnarités induites par la turbulence (amplitudes des paramètres globaux) est essentielle à une bonne étude dynamique.
En effet ces instationnarités interagiront non-linéairement avec les vibrations dues aux mouvements de la structure dans l'étude complète, et doivent être préalablement bien définies. Comme nous le verrons dans la section "Résultats", les résultats obtenus avec les calculs statiques convergent mal, et une étude dynamique ne ferait qu'amplifier cette mauvaise convergence. Il faudrait améliorer les conditions limites et le maillage avant de se lancer dans l'étude dynamique de ce système.

Nous présenterons donc simplement le phénomène physique théorique des instabilités fluide-élatiques.
Selon les observations d'EDF, les tubes sont sensibles aux vibrations quand ils sont soumis à des écoulements transverses. Pour décrire le phénomène d'interaction fluide-structure il faut tout d'abord introduire la vitesse réduite Ur =
f0 D / U  où U désigne la vitesse de l'écoulement, f0 la fréquence d'oscillation du tube et D son diamètre. Quand la vitesse réduite est nulle, le fluide est au repos. Pour les faibles valeurs de la vitesse réduite, les tremblements turbulents sont dominants, et quand Ur augmente, les instabilités dynamiques apparaissent.



La figure précédente présente l'évolution de l'amplitude des vibrations en fonction de la vitesse réduite, et détaille les différentes phases de cette évolution.
- Quand le solide est entouré par un fluide au repos, l'interaction fluide-structure est caractérisée par un ajout de masse et un effet d'amortissement visqueux résultant des forces d'inertie et des forces dissipatives.
- Pour une large gamme de vitesses le tube est excité par la turbulence ("Broad-band excitation"). Dans cette zone, du bruit aléatoire induit des variations de pression aléatoires sur la surface des cylindres et produit les vibrations du tube de faible amplitude. Ces vibrations sont principalement indépendentes du mouvement du tube.
- Dans la zone de "lock-in", les vibrations sont soudainement amplifiées puis diminuent rapidement, ce qui caractérise une synchronisation entre les détachements tourbillon-
aires et les vibrations du tube.
- Enfin, pour les écoulements à grande vitesse, passée une certaine vitesse critique, l'amplitude des vibrations augmente considérablement, pouvant fragiliser et casser le tube. Cette zone correspond à l'instabilité fluide-élastique où les forces dominantes du fluide dépendent du mouvement du tube. Il s'agit d'interactions complexes entre le mouvement de
la structure et la dynamique des tourbillons.

L'instabilité fluide élastique peut être statique ou dynamique et dans ce dernier cas elle
peut être provoquée par deux phénomènes distincts :
- Pour les faibles nombres de Scruton c'est un mécanisme d'amortissement négatif qui génère cette instabilité ("damping dominated instability"). Il s'agit de forces qui dépendent de la vitesse, et l'effet de la transition entre le mouvement de la structure et les changements du fluide joue un rôle important.
- Pour les grands nombres de Scruton, ce sont les déplacements qui provoquent l'instabilité. Ils dépendent de la rigidité du fluide, ou de la constante du ressort qui modélise le forçage.



L'instabilité fluide élastique est liée aux mouvements de la structure au sein du fluide au cours desquels le cylindre extrait de l'énergie au fluide. De nombreux modèles ont été développés pour tenter d'expliquer ce phénomène. Ceux-ci font intervenir un temps de retard entre la réponse du fluide et la vibration de la structure. Tout le problème de l'explication de l'instabulité se résume à quantifier ce temps de retard. En effet, lorsque la vitesse de l'écoulement transverse augmente, le déphasage entre la force appliquée par le fluide sur le tube et le mouvement du tube change de signe, ce qui entraîne un changement de signe de l'amortissement apparent et donc l'instabilité.