Groupe :
NOM Prénom des personnes du groupe:
Le TP doit être rendu sous la forme d'un Notebook jupyter en respectant la nomenclature suivante:
En remplaçant N par le numéro de votre groupe. Tous les résultats, discussions, analyses, doivent donc être inclus dans le fichier.
Vous pouvez vous familiariser avec les notions de Python utiles pour ce TP dans ce Notebook.
On se propose dans se TP de simuler numériquement l'écoulement d'un fluide autour d'un cylindre carré à l'aide d'un code Lattice Boltzmann utilisant le modèle $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:
$c_\alpha$ | $\omega_\alpha$ | $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}}{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 autour d'un cylindre carré est un cas test très utilisé en CFD et en particulier en Lattice Boltzmann car il présente l'avantage d'avoir une géométrie cartésienne simple. Le nombre de Reynolds de cet écoulement est basé sur le coté $D$ du carré: $$ Re=\dfrac{U_0 D}{\nu} $$
Au delà d'une valeur critique $R_e=48$, l'écoulement devient instable et on observe une allée de von Karmann se former, caractérisée par des tourbillons en opposition de phase en aval du cylindre. Pour les faibles nombres de Reynolds ($Re<500$), on peut estimer la distance $\delta$ moyenne entre les tourbillons: $\delta\sim 2.5 D$
Pour les paramètres du calcul, on basera les grandeurs sur le diamètre du cylindre D. On prendra un domaine de $(20D,10D)$, $M_0=0.3$ et les paramètres du modèle D2Q9 décrits plus haut. On placera le centre du cylindre au point $(x_0,y_0)=(8D,5D)$. Le diamètre du cylindre pourra prendre différentes valeurs pendant l'étude.
Le nombre d'itération sera choisi pour atteindre un temps normalisé $t^*=n_t \frac{U_0}{D}=150$.
On souhaite réaliser plusieurs simulations en faisant varier le nombre de Reynolds entre 20 et 150. En prenant 3 valeurs de Reynolds entre les bornes spécifiées: [20, 80, 150]
Calculer les temps de relaxation associés $\tau_g$ et completer les cellules suivantes pour pouvoir effectuer les calculs.
import numpy as np
import matplotlib.pyplot as plt
import time
#Pour pouvoir afficher les figures dans le Notebook
%matplotlib inline
# Parametres du modele D2Q9
Re=... # Nombre de Reynolds (à faire varier)
ca=...# vitesses discrètes
op=... # symetrique de chaque vitesse discrète
w=...# coefficients omega_alpha
c0=...# coefficient de vitesse du son
D=.... # Resolution du cylindre (nombre de point dans le diamètre)
nx,ny=
M0=...
taug=... # Temps de relaxation du modèle (défini en fonction du Reynolds)
nt=... # Nombre d'itération a définir
mask=np.zeros((nx,ny));
# Definition du cylindre => Choisir les bons indices pour placer
#le cylindre au centre du domaine en y et dans le premier tiers en x
cx1,cx2=?,?
cy1,cy2=?,?
mask[cx1:cx2,cy1:cy2]=1
Compléter les fonctions suivantes (remplacer les "...") correspondants aux étapes principales de l'algorithme LBM:
def init(M0):
# fonction d'initialisation des variables macroscopiques
# et de la fonction d'équilibre
rho=...
ux=...
uy=...
geq=np.zeros((nx,ny,9))
eq(rho,ux,uy,geq)
return geq,rho,ux,uy
def collide(gcoll,g,geq,taug):
# Etape de collision
gcoll[:,:,:]=...
def propagate(g,gcoll):
# Etape de propagation
g[:,:,0]=gcoll[...,0]
...
def macro(g,rho,ux,uy):
# calcul des variables macro
rho[:,:]=...
ux[:,:]=...
uy[:,:]=...
def eq(geq,rho,ux,uy):
# Calcul de la fonction d'équilibre
geq[:,:,0]=...
...
def boundary(g,gcoll):
# Conditions aux limites periodiques
# Entree => quelles sont les g inconnues en entree
g[0,:,?]=gcoll[...]
g[0,:,?]=gcoll[...]
g[0,:,?]=gcoll[...]
# Sortie => quelles sont les g inconnues en sortie
g[-1,:,?]=gcoll[...]
g[-1,:,?]=gcoll[...]
g[-1,:,?]=gcoll[...]
# Bas => quelles sont les g inconnues en bas
g[:,0,?]=gcoll[...]
g[:,0,?]=gcoll[...]
g[:,0,?]=gcoll[...]
# Haut => quelles sont les g inconnues en haut
g[:,-1,?]=gcoll[...]
g[:,-1,?]=gcoll[...]
g[:,-1,?]=gcoll[...]
def wall(g,mask):
# Conditions de parois
g[mask==1,1]=...
...
Une fois les paramètres et les fonctions définis, on procède au calcul:
geq,rho,ux,uy=init()
g,gcoll=geq.copy(),geq.copy()
prof=np.zeros(nt) #profil de pression en fonction du temps
start=time.time()
for t in range(nt):
collide(gcoll,g,geq,taug) # Collision
propagate(g,gcoll) # Propagation
boundary(g,gcoll) # Conditions aux limites domaine
wall(g,mask) # Conditions aux limites solides
macro(g,rho,ux,uy) # Calcul des moments
eq(geq,rho,ux,uy) # Calcul de l'equilibre
prof[t]=... # Ecriture des résultats
tcal=time.time()-start
print("Temps de calcul: "+str(tcal))
Réaliser 5 simulations avec des résolutions différentes pour le nombre de Reynolds le plus grand. On pourra prendre $D\in[4,6,8,10,15]$ Tracer dans chaque cas l'évolution du signal temporel de pression en fonction du temp normalisé $t^*$.
fig=plt.figure()
time=
plt.plot(...,label='D=?')
plt.legend()
plt.title();
Représenter le champ de pression et le champ de vitesse longitudinale obtenus à la fin de la simulation pour chacun des 3 Reynolds. (On prendra la résolution la plus grande)
fig=plt.figure()
plt.pcolormesh(...)
plt.title('Reynolds='+str(Re));
fig=plt.figure()
plt.pcolormesh(...)
plt.title('Reynolds='+str(Re));
A partir des résultats précédents, calculer la distance moyenne entre 2 tourbillons et comparer à la valeur théorique donnée plus haut.
Tracer l'évolution temporelle des signaux de pression dans un point du sillage et décrire le résultat obtenu.
Peut-on observer une fréquence particulière dans ces signaux. Commenter.
#Affichage du signal temporelle de pression dans un point du sillage:
fig = plt.figure();
Proposer une méthode pour calculer le coefficient de Trainée du cylindre et donner la valeur obtenue pour les 3 Reynolds choisis.
Résultats:
Re | $C_d$ |
---|---|
--- | --- |
--- | --- |
from IPython.core.display import HTML
style=open('notebooks.css', "r").read()
HTML(style)