PYL 12/10 MF204
http://www.lmm.jussieu.fr/~lagree

Enoncé et corrigé de la PC ENSTA MF204

PC n°5 Résolution de Poiseuille Turbulent

Le problème de départ est l'écoulement turbulent en 2D entre deux plans infinis distants de 2h. Un gradient de pression pousse le fluide. Les hypothèses de moyenne, d'invariance par translation aboutissent à la résolution de l'équation de Poiseuille turbulent en 2D pour le champ moyen U (et V est toujours nulle), U donnée nulle en y=0, condition de symétrie en y=h. Consulter la Petite Classe pour plus de détails sur les hypothèses, consulter aussi le cours sur les transferts turbulents.

Dans le cas turbulent, on a le tenseur de Reynolds qui donne la contribution -<u'v'> que l'on modélise par la longueur de mélange de Prandtl ( κ=0.41 constante de Kármán)
-<u'v'> = ρ κy²(∂U/∂y)²

On peut alors montrer que l'équation de Navier Stokes moyenne turbulente se réduit à
0 = (1/ρ) ∂( μ ∂U/∂y + κy²(∂U/∂y)²)/∂y -(1/ρ) dP₀/dx
U(0)=0 et ∂U/∂y(h)=0

On a une boîte de taille h (0<y<h) et (-h/2<x<h/2)
Ici on a pour source le gradient de pression constant -(1/ρ)dP₀/dx que l'on pose égal à u_*²/h
On utilise u_* et h pour adimensionner. Sans dimension on a:
0 = ∂( 1/R_* ∂u/∂y + κy²(∂u/∂y)²)/∂y + 1
u(0)=0 et ∂u/∂y(1)=0
avec R_*=ρu_*h/μ

Pour les valeurs de y d'ordre 1, le terme 1/R_* est négligeable, on a (car - <u'v'> = κy²(∂u/∂y)²) :
- <u'v'> = 1 - y

Pour les petites valeurs de y on pose y⁺ =R_*y et u⁺ =u, le problème sans dimension (dans les nouvelles échelles, dites de paroi) est après avoir négligé 1/R_* et intégré une première fois
∂u⁺/∂y⁺ + κy⁺ ²(∂u⁺/∂y⁺ )² = 1.
reste alors u⁺ =y⁺ pour les petits y⁺. Pour les valeurs quelconques, on extrait le gradient, le problème sans dimension est :
∂u⁺/∂y⁺ = (-1 + √(1 + 4 κy⁺ ²))/(2 κy⁺ ² )
qui s'intègre en
u⁺ = (1 - √[1 + 4 κy⁺ ²])/ ( 2 κy⁺ ²) + log[2 κy⁺ + √[1 + 4 κy⁺ ²]]/κ
qui est une "solution exacte" à grand R_*. Cette expression permet de retrouver u⁺ =y⁺ pour les petits y⁺, et pour les y⁺ assez grands, on a la loi log approchée
u⁺ = log[y⁺]/κ + (-1 + log[4.] + log[κ])/κ soit u⁺ = 2.4390 log[y⁺] -1.23245

Nous allons vérifier ces résultats avec gerris

gerris2D: résout

∂U/∂t + u ∂U/∂x + v ∂U/∂y = alpha { -∇ p + ∇. (   SourceViscosity (∇ U + ∇U^T) ) } + Source U

où U= (u,v) et alpha est 1/ρ ici ρ=1
SourceViscosity est "normalement" μ la viscosité dynamique, mais nous allons la redéfinir car ici la viscosité dépend du gradient de vitesse.

Dans le fichier gerris2D turb.gfs, on définit Rstar 500 et on code en dur la solution analytique
Define ulog(y) log(0.82*y + pow(1 + 4*pow(0.41,2)*pow(y,2),0.5))*pow(0.41,-1) + (pow(0.41,-2)*pow(y,-1)*(1 - pow(1 + 4*pow(0.41,2)*pow(y,2),0.5)))/2.
On rafine pour avoir au moins un point en dessous de y=1
      Refine {
      if (fabs(y) <= 0.05 )
    return 9;
      if (fabs(y) > 0.05 && fabs(y) <= 0.1)
    return 8;
    if (fabs(y) > 0.10 && fabs(y) <= 0.2)
    return 7;
      if (fabs(y) > 0.2 && fabs(y) < 0.4)
      return 6;
      return 5;
}
Source U est mis à 1
On ne met pas de terme de dérivée totale pour accélérer le calcul
# use backward Euler to avoid Crank-Nicholson oscillations in time
# we need this so that acceleration can be balanced by viscous stress
    AdvectionParams { gc = 1 }
    AdvectionParams { scheme = none }

on adapte la viscosité en mettant dans SourceViscosity:
Mu = 1./Rstar + 0.41*0.41*sqrt(2.)*(y*y)*D2;

Conditions aux limites,
GfsBox {
    bottom = Boundary {
        BcDirichlet U 0
        BcDirichlet V 0
    }
    top = Boundary {
        BcNeumann U 0
        BcDirichlet V 0
    }
dont périodicité annoncée par
1 GfsSimulation GfsBox GfsGEdge {
est conclue par
1 1 right

Mise en oeuvre de gerris2D:
Le script run.sh lance le calcul.
Le fichier gerris est turb.gfs .
Le fichier gfsview est v.gfv .
Les courbes suivantes sont obtenues par gnuplot grâce au script à gauche en semi log le profil calculé comparé à la solution exacte et à la solution linéaire. A droite -<u'v'> lèche bien la droite.

Pour aller plus loin on codera une longueur de mélange de Van Driest et de Nikuradse.
On ajoutera aussi l'équation de la chaleur: fichier gerris est turbT.gfs qui donne l'image suivante, profil de vitesse en vert, de température en rouge (T=0 à la paroi flux unité en y=h).

Pour comparer, une résolution numérique par différences finies de l'équation, codée en Java.
Mise en oeuvre de Java une A-pyl-quette qui calcule le profil turbulent a aussi été codée.
on retrouve à peu près la même chose que gerris, l'interface ressemble à ceci:

OK 12/2010, page ouaibe de cette page http://www.lmm.jussieu.fr/~lagree/COURS/ENSTA/HTMLPC/PC5/index.html