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 d'étudier l'évolution d'un pulse de Pression à l'aide d'un schéma LBM $D2Q4$.
Le modèle à 4 vitesses de la méthode de Boltzmann sur réseau est un modèle linéaire permettant de résoudre les équations de propagation d’ondes - type d’Alembert ou Helmholtz - et son algorithme s’exprime de façon très simple. Les paramètres du modèle sont résumés dans le tableau suivant:
| $c_\alpha$ | $\omega_\alpha$ | $c_0$ | $\tau_g$ |
| (1,0) (0,1) (-1 0) (0 -1) | $\dfrac{1}{4},\dfrac{1}{4},\dfrac{1}{4},\dfrac{1}{4}$ | $\dfrac{1}{\sqrt{2}}$ | $\dfrac{1}{2}$ |
La fonction d'équilibre du modèle s'écrit:
$$ \displaystyle{g_{\alpha}^{eq}=\rho\omega_\alpha +\omega_\alpha\dfrac{(g_1-g_3).c_{\alpha,x}+(g_2-g_4).c_{\alpha,y}}{c_0^2}} $$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}^4 g_{\alpha} $$
$$ \rho u_x=\sum_{\alpha=1}^4 c_{\alpha,x}g_{\alpha} $$$$ \rho u_y=\sum_{\alpha=1}^4 c_{\alpha,y}g_{\alpha} $$Le pulse de pression est un cas test très utilisé en simulation numérique car on connait une solution analytique. Ainsi il devient très utile lorsque l'on souhaite évaluer quantitativement les capacités d'un schéma numérique. Ce cas test consiste à initialiser le champ de pression avec une gaussienne centrée autour d'un point $(x_0,y_0)$ et le champ de vitesse à 0:
$$ \rho(x,y,t_0)=1.+\varepsilon\exp\left(-\alpha\eta^2 \right) $$et
$$ u(x,y,t_0)=0 $$Avec $\eta = \sqrt{(x-x_0)^2 +(y-y_0)^2}$, $\alpha=\dfrac{ln(2)}{b^2}$ où $b$ est la demie-largeur de la gaussienne. La solution analytique de ce cas test est donnée en fonction du temps par:
$$ \displaystyle{\rho (x,y,t) = 1.+\dfrac{\varepsilon}{2\alpha}\int_0^{\infty} \xi exp\left[-\dfrac{\xi^2}{4\alpha}\right] cos(t\xi)J_0(\xi \eta) d\xi} $$et
$$ \displaystyle{u_x(x,y,t) = \dfrac{(x-x_0)\varepsilon}{2\alpha\eta}\int_0^{\infty} \xi exp\left[-\dfrac{\xi^2}{4\alpha}\right] sin(t\xi)J_1(\xi \eta) d\xi} $$avec les fonction de Bessel sphérique de type $1$, $J_0$ d'ordre $0$ et $J_1$ d'ordre 1.
Pour les paramètres du calcul, on prendra un domaine de $(400,400)$ mailles avec $\tau_g=0.5$, $N_t=240$, $b=12$ et les paramètres du modèle D2Q4 décrits plus haut.
Completer les cellules suivantes pour pouvoir effectuer ce calcul.
#Import des modules utiles
import numpy as np
import matplotlib.pyplot as plt
import time
#Pour pouvoir afficher les figures dans le Notebook:
%matplotlib inline
# Parametres du modele D2Q4
ca=...# vitesses discrètes
w=...# coefficients omega_alpha
c0=...# coefficient de vitesse du son
nx,ny,taug,b,nt=...# a définir
Compléter les fonctions suivantes correspondants aux étapes principales de l'algorithme LBM:
def init():
# fonction d'initialisation des variables macroscopiques
# et de la fonction d'équilibre
geq=np.zeros((nx,ny,4))
ux,uy=np.zeros((nx,ny)),np.zeros((nx,ny))
alpha=np.log(2.)/b**2
x,y=np.arange(nx),np.arange(ny)
rho=... # On pourra ici utiliser la fonction einsum
eq(geq,rho,ux,uy)
return geq,rho,ux,uy
def collide(gcoll,g,geq,taug):
# Etape de collision
gcol[:,:,:]=...
def propagate(g,gcoll):
# Etape de propagation
g[?,?,0]=gcoll[?,?,0]
g[?,?,1]=gcoll[?,?,1]
g[?,?,2]=gcoll[?,?,2]
g[?,?,3]=gcoll[?,?,3]
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]=...
geq[:,:,1]=...
geq[:,:,2]=...
geq[:,:,3]=...
def boundary(g,gcoll):
# Conditions aux limites periodiques
# Entree
g[0,:,?]=...
#Sortie
g[-1,:,?]=...
#Bas
g[:,0,?]=...
#Haut
g[:,-1,?]=...
Une fois les paramètres et les fonctions définis, on procède au calcul à proprement parler:
geq,rho,ux,uy=init()
g,gcoll=geq.copy(),geq.copy()
prof=np.zeros((nt,nx,2))#profil de pression et vitesse en fonction du temps
for t in range(nt):
collide(gcoll,g,geq,taug) # Collision
propagate(g,gcoll) # Propagation
boundary(g,gcoll) # Conditions aux limites
macro(g,rho,ux,uy) # Calcul des moments
eq(geq,rho,ux,uy) # Calcul de feq
prof[t,:,0],prof[t,:,1]=... # Ecriture des résultats
Représenter la solution numérique obtenus. On pourra par exemple tracer $\rho$ et $u_x$ à l'aide de la fonction pcolormesh
#Affichage du resultat
fig=plt.figure()
fig.add_subplot(121)
plt.pcolormesh(rho.T)
fig.add_subplot(122)
plt.pcolormesh(ux.T)
Calculer et tracer la solution théorique donnée plus haut à l'aide des fonctions integrande1 et integrande2.
from scipy.integrate import quad
from scipy.special import j0,j1
def integrande1(xi):
global alpha,eta,t
y=xi*np.exp(-xi**2/(4.*alpha))*np.cos(t*xi)*j0(xi*eta)
return y
def integrande2(xi):
global alpha,eta,t
y=xi*np.exp(-xi**2/(4.*alpha))*np.sin(t*xi)*j1(xi*eta)
return y
# Utiliser la fonction quad(f,a,b) pour intégrer la fonction f entre a et b
rho_th=...
ux_th=...
fig=plt.figure()
fig.add_subplot(121)
plt.plot(rho_th,'k',label='Theorique')
#plt.plot(rho_num,'r',label='D2Q4')
plt.legend()
fig.add_subplot(122)
plt.plot(ux_th,'k',label='Theorique')
#plt.plot(ux_num,'r',label='D2Q4')
plt.legend()
Comparer les profils de pression et de vitesse numérique et théorique. Commenter les résultats obtenus
En remarquant que la fonction d'équilibre est une version linéaire, on remarque qu'il devient possible d'écrire l'algorithme général sous une forme très simple. Proposer une écriture explicite en simplifiant les termes linéaire et montrer que l'on peut écrire l'algorithme en utilisant uniquement des sommes de $f$.
from IPython.core.display import HTML
style=open('notebooks.css', "r").read()
HTML(style)