Le TP doit être réalisé en binôme ou individuellement et doit être rendu sous la forme d'un Notebook jupyter en respectant la nomenclature suivante:
Tous les résultats, discussions, analyses, doivent donc être inclus dans le fichier.
import numpy as np
import matplotlib.pyplot as plt
import time
fs=20
plt.style.use('seaborn-dark')
plt.rc('xtick',labelsize=fs)
plt.rc('ytick',labelsize=fs)
#plt.rc('text', usetex=True)
%matplotlib inline
On se propose dans ce TP de simuler numériquement l'écoulement induit par viscosité dans une cavité de section carré. Pour pouvoir simuler les phénomènes visqueux, il est nécessssaire d'utiliser la quadrature d'ordre 3 et donc le réseau $D2Q9$.
Le modèle à 9 vitesses de la méthode de Boltzmann sur réseau assure l'égalitité des moments jusqu'à l'ordre 2 et permet donc la simulation des écoulements visqueux régit par les équations de Navier-Stokes avec une limitation en nombre de Mach. Les paramètres du modèle sont résumés dans le tableau suivant:
| $\widehat{c_\alpha}$ | $\omega_\alpha$ | $\widehat{c_0}$ | $\tau_g$ |
| $(0,0)\\ (1,0) (0,1) (-1,0) (0,-1)\\ (1,1) (-1,1) (-1,-1) (1,-1)$ | $\dfrac{4}{9},\\ \dfrac{1}{9},\dfrac{1}{9},\dfrac{1}{9},\dfrac{1}{9},\\ \dfrac{1}{36},\dfrac{1}{36},\dfrac{1}{36},\dfrac{1}{36}$ | $\dfrac{1}{\sqrt{3}}$ | $\dfrac{1}{2}+\dfrac{\tilde{\nu}}{\widehat{c_0}^2}$ |
La fonction d'équilibre du modèle s'écrit: $$ \displaystyle{g_{\alpha}^{eq}=\rho\omega_\alpha\left(1+\dfrac{c_{\alpha,i}u_i}{c_0^2}+\dfrac{(c_{\alpha,i}u_i)^2}{2c_0^4}-\dfrac{u_i^2}{2c_0^2}\right)} $$
et son algorithme général: $$ \displaystyle{g_{\alpha}^{coll} = g_{\alpha}-\dfrac{1}{\tau_g}[g_{\alpha}-g_{\alpha}^{eq}]\\ g_{\alpha}(x,y) = g_{\alpha}^{coll}(x-c_{\alpha,x},y-c_{\alpha,y})} $$
Les variables macroscopiques sont calculées à partir des moments des fonctions de distribution: $$ \rho=\sum_{\alpha=1}^9 g_{\alpha} $$
$$ \rho u_x=\sum_{\alpha=1}^9 c_{\alpha,x}g_{\alpha} $$$$ \rho u_y=\sum_{\alpha=1}^9 c_{\alpha,y}g_{\alpha} $$L'écoulement dans une cavité entrainée consite a imposer une vitesse $U_0$ constante sur le bord supérieur d'une cavité de largeur $L$ et de hauteur $L$. Le nombre de Reynolds de cet écoulement est basé sur la largeur $L$ et la vitesse d'entrainement $U_0$: $$ Re=\dfrac{U_0 L}{\nu}=\dfrac{\tilde{U_0} \tilde{L}}{\tilde{\nu}} $$
Figure 1: Géométrie du domaine de calcul
Dans ce TP nous allons nous focaliser sur les solutions stationnaires de cet écoulement pour des nombres de Reynolds compris entre 100 et 1000.
Pour les paramètres du calcul, on prendra un domaine unitaire $(L_x, L_y)=(1,1)$ de résolution $(n_x,n_y)=(50,50)$ mailles (on fera varier la résolution pour les dernières questions). On souhaite réaliser plusieurs simulations en faisant varier le nombre de Reynolds. Le nombre de Mach sera fixé à $M_0=0.2$.
1. Completer la cellule suivante à partir des données fournies et des définitions du cours.
ca=...
w=...
c0=...
# Domaine
nx=...
ny=...
# Taille de maille:
dx=...
# Vitesse du son réelle à 293K:
c0_real=...
# Pas de temps:
dt=...
# Mach:
M0=...
2. Compléter chacune des fonctions suivantes correspondant aux étapes de l'algorithme
Initialisation: Définir la fonction permettant d'initialiser le domaine à une densité uniforme $\rho=1$, à vitesse nulle dans la cavité et à la vitesse $U_0$ sur le bord supérieur. Calculer également le paramètre de relaxation $\tau_g$ en fonction du nombre de Reynolds souhaité:
def init(M0,Re):
# Initialisation du domaine
...
...
return geq,rho,ux,uy,taug,U0
Equilibre: Dans ce TP nous avons besoin de calculer explicitement les fonctions de distribution à l'équilibre. Définir la fonction permettant la mise à jour de ces distributions en fonction des moments $\rho$, $u_x$ et $u_y$.
def eq(rho,ux,uy,geq):
# Mise à jour de geq
...
...
...
Collision: Ecrire simplement l'étape de collision:
def collide(gcoll,g,taug):
# Etape de collision
...
...
...
Propagation: Pour cette étape on pourra utiliser du slicing plutôt qu'une boucle sur l'espace.
def propagate(g,gcoll):
# Etape de propagation
...
...
...
...
Calcul des moments: Calcul les variables macroscopiques $\rho$,$u_x$,$u_y$
def macro(g,rho,ux,uy):
# calcul des variables macro
...
...
...
Condition aux limites solide: Définir la fonction permettant de mettre à jour les distributions inconnues sur les parois solides. On utilisera ici la technique du Bounce-Back présentée dans le cours:
def wall(gcoll,g,mask):
# Paroi solide sur les bords latéraux et inférieur: Bounce back
...
...
...
Condition aux limites fluide: Définir la fonction permettant de mettre à jour les moments $u_x, \rho_{top}$ et les distributions inconnues sur le bord supérieur. On utilisera les approches vues dans le cours pour le bounce-back hors équilibre:
$$ g_?-g_?^{eq}=g_\overline{?}-g_\overline{?}^{eq} $$Ici $\overline{?}$ désigne la vitesse opposée à la vitesse $?$. En identifiant les distributions inconnues sur le bord supérieur, simplifier l'expression ci-dessus en fonction des moments connu sur le bord supérieur:
$g_?=...$
$g_?=...$
$g_?=...$
Exprimer la densité sur le bord supérieur $\rho_{top}$ en fonction des distributions connues et de la vitesse
$\rho_{top}=...$
def lid(g,rho,ux,uy,M0):
# Conditions de vitesse au bord superieur
# Macro: On impose ux et on en déduit rho_w:
...
...
# Distributions: On impose les distributions inconnues avec le Bounce-Back hors equilibre
...
...
...
3. Une fois les paramètres et les fonctions définis, fixer le nombre de Reynolds (on commencera par $Re=100$), initialiser les distributions à leur valeur d'équilibre puis définir les mailles concernées par une condition aux limites solide. Pour cela on pourra utiliser une variable de type mask de la même taille que le domaine et valant 0 pour une maille fluide et 1 pour une maille solide.
# Nombre de Reynolds:
Re=...
# initialisation:
geq,rho,ux,uy,taug,U0=init(M0,Re)
g,gcoll=geq.copy(),geq.copy()
# Marquage des conditions aux limites:
# Fluide: mask=0
# Solide: mask=1
mask=np.zeros((nx,ny));
mask[?,?]=1 # left
mask[?,?]=1 # right
mask[?,?]=1 # bottom
On recherche une solution stationnaire pour les paramètres choisis, c'est à dire une solution dont les quantités n'évoluent plus au cours du temps. Il faut donc utiliser une boucle while basée sur un critère d'arrêt.
4. Définir une boucle de calcul reprenant l'ensemble des fonctions de la question 2. et basée sur la convergence de l'énergie cinétique moyenne dans la cavité:
$$ \displaystyle{E_c=\dfrac{1}{V}\int_\mathcal{V}\rho\left(\dfrac{u_x^2}{U_0^2}+\dfrac{u_y^2}{U_0^2} \right) d\mathcal{V}} $$On prendra comme critère de convergence $\epsilon=\left|Ec^{n+1}-Ec^{n} \right|<10^{-5}$
Ec=...
eps=1
nt=0
start=time.time()
while (eps>1e-5):
nt+=1
...
...
...
...
...
Ec=...
eps=...
tcal=time.time()-start
print(str(nt)+" itérations en "+str(tcal)+"s: Performances: "+str(...)+" MLUPS")
5. Comparer les performances du code obtenu à celles du TP1. Commenter
6. Tracer les lignes de courant de l'écoulement et commenter la figure obtenue. (On pourra utiliser la fonction streamplot)
fig=plt.figure(figsize=(16,16))
plt.streamplot(?,?,?,?)
plt.xlabel(r'$X/L$',fontsize=fs);plt.ylabel(r'$Y/L$',fontsize=fs)
plt.title(r'$Re='+str(Re)+'$',fontsize=fs);plt.axis('equal');
Les fichiers Ghia82data_ux.csv et Ghia82data_uy.csv (à télécharger sur le Moodle) contiennent des résultats de référence enregistrés le long des lignes $x=L/2$ et $y=L/2$ respectivement. Les fichiers sont constitués de plusieurs colonnes. La première correspond à la coordonnée du point de mesure et les autres aux différentes valeurs de Reynolds:
| $x/y$ | $u_{x/y}$ $Re=100$ | $u_{x/y}$ $Re=400$ | $u_{x/y}$ $Re=1000$ | $u_{x/y}$ $Re=3200$ | $u_{x/y}$ $Re=5000$ | $u_{x/y}$ $Re=7500$ | $u_{x/y}$ $Re=10000$ |
|---|
Les grandeurs de vitesse sont normalisées par la vitesse d'entrainement $U_0$.
Ghia=np.loadtxt('Ghia82data_ux.csv')
7. Comparer les profils obtenus pour $u_x$ et $u_y$. On pourra faire varier le critère de convergence et commenter son influence sur les résultats.
fig=plt.figure(figsize=(16,16))
plt.subplot(121)
...
plt.subplot(122)
...
8. Etudier la précision des résultats en fonction du nombre de Mach pour $R_e=100$. Commenter
Commentaires
8. Compléter le tableau suivant en indiquand le nombre d'itération nécessaire à la convergence en faisant varier le Reynolds et la résolution. Que constatez-vous pour les faibles résolutions et les hauts Reynolds. Expliquez
| Nx\Re | 100 | 400 | 1000 | 3200 |
|---|---|---|---|---|
| 50 | ? | ? | x | x |
| 100 | ? | ? | x | x |
| 200 | ? | ? | ? | ? |
| 400 | ? | ? | ? | ? |
Commentaires
Présenter ici la synthèse de votre TP en décrivant les points importants et les principaux résultats.
from IPython.core.display import HTML
style=open('notebooks.css', "r").read()
HTML(style)