[Python] Observons le vortex de Kalman par la méthode de Boltzmann en treillis. [Simulation de fluide]
introduction
Cette fois, je voudrais faire une simulation de fluide comme le gif ci-dessus en utilisant la méthode de Boltzmann en treillis.
** Qu'est-ce que la méthode Boltzmann en treillis **
> La méthode de Boltzmann en treillis (méthode de Koshi Boltzmann) approche un fluide avec un agrégat d'un grand nombre de particules virtuelles (modèle de gaz en treillis), et développe séquentiellement la collision et la translation de chaque particule en utilisant la fonction de distribution de vitesse des particules. Il s'agit d'une méthode de simulation permettant d'obtenir le champ de flux thermique d'un fluide en calculant le moment. -Wikipédia
Aussi, quand je l'ai écrit pour être aussi court que possible cette fois, c'était plus lent que l'exemple de code sur le site de référence (lien à la fin de l'article). Même Numpy, qui a un calcul matriciel rapide, semble être inutile.
Méthode
Ici, je vais omettre la théorie détaillée et expliquer le contenu nécessaire à la mise en œuvre.
La formule à calculer s'affiche en premier.
n_i'(r,t+1)=n_i(r+e_i,t)\tag{1}\\
n_i(r,t+1) = (1-\omega)n_i'(r,t+1)+\omega\bigg[\rho w_i(1+3\vec{e_i}\cdot\vec{u}+\frac{9}{2}(\vec{e_i}\cdot\vec{u})^2-\frac{3}{2}|\vec{u}|^2)\bigg] \tag{2}
ici,\\
n_i est le i composant de la fonction de distribution des particules,\\
\vec{e_i}Est le i-ème vecteur de position,\\
w_i=\begin{cases}
\frac{4}{9} & ( i = 0 ) \\
\frac{1}{9} & ( 1\le i \le 4 )\\
\frac{1}{36} & (5 \le i \le 8)
\end{cases}\\
\rho =\sum_{i} n_i(t)\\
\vec{u}Est en deux dimensions x,y vecteur de vitesse à 2 voies,\\
\omega = 1/\tau ， \tau est un facteur de relaxation unique.
Calculez ceci pour chaque $ i $, chaque point.
Tout d'abord, $ \ vec e_i $, qui est le vecteur de position vers votre propre carré + 8 carrés adjacents lorsque votre carré est l'origine. $ i = 0 $ est votre masse $ \ vec e_0 = (0, 0) $, $ 1 \ le i \ le 4 $ est haut, bas, gauche et droite $ \ vec e_ {1 ... 4} = (0,1), (0, -1), (1,0), (-1,0) $, $ 5 \ le i \ le 8 $ est adjacent au nom $ \ vec e_ {5 ... 8} = (1,1), (-1, -1), (1, -1), ( -1,1) $ Représente. Par conséquent, la formule $ (1) $ peut être interprétée comme représentant l'afflux de la masse adjacente.
Vient ensuite $ \ vec e_i \ cdot \ vec u $.
Puisque $ \ vec u = (u_x, u_y) $, par exemple, si $ i = 7 $,
finalement
Gestion des obstacles
Ce n'est pas très intéressant juste de couler, et le vortex de Kalman du sujet principal ne se produit pas, alors pensons au traitement lorsqu'un obstacle est installé. Cependant, c'est facile avec la méthode de Boltzmann en treillis (je pense que c'était un peu gênant lorsque j'ai approché l'équation de Navier-Stokes). À l'origine, il est calculé séparément pour chaque direction d'entrée, donc s'il rencontre un obstacle, il peut être acheminé vers le composant opposé à la direction d'entrée.
Passons maintenant à la mise en œuvre.
la mise en oeuvre
Définition des fonctions
Lattice_Boltzmann_Method.py
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from ipywidgets import interact
w = np.array([4 / 9, 1 / 9, 1 / 9, 1 / 9, 1 / 9, 1 / 36, 1 / 36, 1 / 36, 1 / 36])
e = np.array([[0, 0], [0, 1], [0, -1], [1, 0], [-1, 0], [1, 1], [-1, -1],[1, -1], [-1, 1]])
def stream(): #(1)Équation + traitement des collisions avec obstacles
global n
for i, axis in zip(range(1, 5),[[1, 1, 0, 1], [1, -1, 0, -1], [0, 1, 1, -1], [0, -1, 1, 1]]): #(1)formule
n[:, :, i] = np.roll(n[:, :, i], axis[1], axis=axis[0])
n[:, :, i + 4] = np.roll(np.roll(n[:, :, i + 4], axis[1],axis=axis[0]),axis[3],axis=axis[2])
for i in range(1, 9):
n[:, :, i][barrier[:, :, i]] = n[:, :, i - (-1)**i][barrier[:, :, 0]] #La quantité de collision est envoyée au composant dans la direction opposée.
def collide(): #(2)Équation + afflux de l'extrémité gauche
global n, u, u9
rho = n.sum(axis=2).reshape(H, W, 1) #Calcul de ρ
u = (n @ e) / rho #(ux,uy)Calculs de
n = (1 - om) * n + om * w * rho * (1 - 1.5 *(u**2).sum(axis=2).reshape(H, W, 1) +3 * (u @ e.T) + 4.5 * (u @ e.T)**2) #(2)formule
n[:, 0] = w * (1 - 1.5 * (u_o**2).sum(axis=1).reshape(H, 1) + 3 *(u_o @ e.T) + 4.5 * (u_o @ e.T)**2) #Flux entrant de l'extrémité gauche
def curl(): #Fonctions liées à la couleur lors du traçage
return np.roll(u[:, :, 0], -1, axis=1) - np.roll(u[:, :, 0], 1, axis=1) - np.roll(u[:, :, 1], -1, axis=0) + np.roll(u[:, :, 1], 1, axis=0)
def mk_barrier(arr): #Un enfant qui calcule la correspondance de la destination du rebond au moment de la collision lors du passage d'un tableau représentant des obstacles
global barrier
barrier = np.zeros((H, W, 9), bool)
barrier[:, :, 0] = arr
for i, axis in zip(range(1, 5),[[1, 1, 0, 1], [1, -1, 0, -1], [0, 1, 1, -1], [0, -1, 1, 1]]):
barrier[:, :, i] = np.roll(barrier[:, :, 0], axis[1], axis=axis[0])
barrier[:, :, i + 4] = np.roll(barrier[:, :, i], axis[3], axis=axis[2])
Définition et initialisation des paramètres
Lattice_Boltzmann_Method.py
H, W = 40, 100 #Vertical et horizontal
viscosity = .005 #viscosité
om = 1 / (3 * viscosity + 0.5) #Ce calcul suit le site de référence. Cela semble être l'expression relationnelle entre viscosité et ω.
u0 = [0, .12] #Vitesse initiale et vecteur de vitesse d'entrée de l'extrémité gauche(-uy,ux)Correspond à
u = np.ones((H, W, 2)) * u0
u_o = np.ones((H, 2)) * u0
n = w * (1 - 1.5 * (u**2).sum(axis=2).reshape(H, W, 1) + 3 * (u @ e.T) + 4.5 * (u @ e.T)**2)
Obstacle
Lattice_Boltzmann_Method.py
r = 3 #La longueur du mur ou le rayon du cercle
#mur
wall = np.zeros((H, W), bool)
wall[H // 2 - r:H // 2 + r, H // 2] = True
#Cercle
x, y = np.meshgrid(np.arange(W), np.arange(H))
circle = ((x - H / 2)**2 + (y - H / 2)**2 <= r**2)
mk_barrier(circle)
terrain
Lattice_Boltzmann_Method.py
def update(i):
for t in range(10):
stream()
collide()
im.set_array(curl() - barrier[:, :, 0])
return im,
"""Devenez interactif pour Jupyter
@interact(Viscosity=(0.05, 2.0, 0.01),v0=(0, 0.12, 0.01),Barrier=["circle", "wall"])
def a(Viscosity, v0, Barrier):
global viscosity, om, u0, u_o
viscosity = Viscosity * .1
om = 1 / (3 * viscosity + 0.5)
u0 = [0, v0]
u_o = np.ones((H, 2)) * u0
mk_barrier(eval(Barrier))
"""
fig = plt.figure()
stream()
collide()
im = plt.imshow(curl() - barrier[:, :, 0],cmap="jet",norm=plt.Normalize(-.1, .1))
anim = FuncAnimation(fig, update, interval=1, blit=True)
plt.show()
résultat
Identique à celui montré au début,
Résumé
Je suis impressionné que vous puissiez reproduire des phénomènes naturels avec votre propre code. Veuillez essayer.
<détails> Lattice_Boltzmann_Method.py
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from ipywidgets import interact
w = np.array([4 / 9, 1 / 9, 1 / 9, 1 / 9, 1 / 9, 1 / 36, 1 / 36, 1 / 36, 1 / 36])
e = np.array([[0, 0], [0, 1], [0, -1], [1, 0], [-1, 0], [1, 1], [-1, -1],[1, -1], [-1, 1]])
def stream(): #(1)Équation + traitement des collisions avec obstacles
global n
for i, axis in zip(range(1, 5),[[1, 1, 0, 1], [1, -1, 0, -1], [0, 1, 1, -1], [0, -1, 1, 1]]): #(1)formule
n[:, :, i] = np.roll(n[:, :, i], axis[1], axis=axis[0])
n[:, :, i + 4] = np.roll(np.roll(n[:, :, i + 4], axis[1],axis=axis[0]),axis[3],axis=axis[2])
for i in range(1, 9):
n[:, :, i][barrier[:, :, i]] = n[:, :, i - (-1)**i][barrier[:, :, 0]] #La quantité de collision est envoyée au composant dans la direction opposée.
def collide(): #(2)Équation + afflux de l'extrémité gauche
global n, u, u9
rho = n.sum(axis=2).reshape(H, W, 1) #Calcul de ρ
u = (n @ e) / rho #(ux,uy)Calculs de
n = (1 - om) * n + om * w * rho * (1 - 1.5 *(u**2).sum(axis=2).reshape(H, W, 1) +3 * (u @ e.T) + 4.5 * (u @ e.T)**2) #(2)formule
n[:, 0] = w * (1 - 1.5 * (u_o**2).sum(axis=1).reshape(H, 1) + 3 *(u_o @ e.T) + 4.5 * (u_o @ e.T)**2) #Flux entrant de l'extrémité gauche
def curl(): #Fonctions liées à la couleur lors du traçage
return np.roll(u[:, :, 0], -1, axis=1) - np.roll(u[:, :, 0], 1, axis=1) - np.roll(u[:, :, 1], -1, axis=0) + np.roll(u[:, :, 1], 1, axis=0)
def mk_barrier(arr): #Un enfant qui calcule la correspondance de la destination du rebond au moment de la collision lors du passage d'un tableau représentant des obstacles
global barrier
barrier = np.zeros((H, W, 9), bool)
barrier[:, :, 0] = arr
for i, axis in zip(range(1, 5),[[1, 1, 0, 1], [1, -1, 0, -1], [0, 1, 1, -1], [0, -1, 1, 1]]):
barrier[:, :, i] = np.roll(barrier[:, :, 0], axis[1], axis=axis[0])
barrier[:, :, i + 4] = np.roll(barrier[:, :, i], axis[3], axis=axis[2])
H, W = 40, 100 #Vertical et horizontal
viscosity = .005 #viscosité
om = 1 / (3 * viscosity + 0.5) #Ce calcul suit le site de référence. Cela semble être l'expression relationnelle entre viscosité et ω.
u0 = [0, .12] #Vitesse initiale et vecteur de vitesse d'entrée de l'extrémité gauche(-uy,ux)Correspond à
u = np.ones((H, W, 2)) * u0
u_o = np.ones((H, 2)) * u0
n = w * (1 - 1.5 * (u**2).sum(axis=2).reshape(H, W, 1) + 3 * (u @ e.T) + 4.5 * (u @ e.T)**2)
r = 3 #La longueur du mur ou le rayon du cercle
#mur
wall = np.zeros((H, W), bool)
wall[H // 2 - r:H // 2 + r, H // 2] = True
#Cercle
x, y = np.meshgrid(np.arange(W), np.arange(H))
circle = ((x - H / 2)**2 + (y - H / 2)**2 <= r**2)
mk_barrier(circle)
def update(i):
for t in range(10):
stream()
collide()
im.set_array(curl() - barrier[:, :, 0])
return im,
"""Devenez interactif pour Jupyter
@interact(Viscosity=(0.05, 2.0, 0.01),v0=(0, 0.12, 0.01),Barrier=["circle", "wall"])
def a(Viscosity, v0, Barrier):
global viscosity, om, u0, u_o
viscosity = Viscosity * .1
om = 1 / (3 * viscosity + 0.5)
u0 = [0, v0]
u_o = np.ones((H, 2)) * u0
mk_barrier(eval(Barrier))
"""
fig = plt.figure()
stream()
collide()
im = plt.imshow(curl() - barrier[:, :, 0],cmap="jet",norm=plt.Normalize(-.1, .1))
anim = FuncAnimation(fig, update, interval=1, blit=True)
plt.show()