[Python] Beobachten wir den Kalman-Wirbel nach der Gitter-Boltzmann-Methode. [Flüssigkeitssimulation]
Einführung
Dieses Mal möchte ich eine Flüssigkeitssimulation wie das obige GIF mit der Gitter-Boltzmann-Methode durchführen.
** Was ist die Gitter-Boltzmann-Methode **
> Die Gitter-Boltzmann-Methode (Koshi-Boltzmann-Methode) approximiert eine Flüssigkeit mit einem Aggregat einer großen Anzahl virtueller Partikel (Gittergasmodell) und entwickelt nacheinander die Kollision und Translation jedes Partikels unter Verwendung der Geschwindigkeitsverteilungsfunktion der Partikel. Dies ist eine Simulationsmethode, um das Wärmestromfeld eines Fluids durch Berechnung des Moments zu erhalten. -Wikipedia
Als ich es diesmal so kurz wie möglich schrieb, war es langsamer als der Beispielcode auf der Referenzseite (Link am Ende des Artikels). Sogar Numpy, der eine schnelle Matrixberechnung hat, scheint nutzlos zu sein.
Methode
Ich werde hier auf die detaillierte Theorie verzichten und die für die Implementierung erforderlichen Inhalte erläutern.
Die zu berechnende Formel wird zuerst angezeigt.
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}
Hier,\\
n_i ist die i-Komponente der Partikelverteilungsfunktion,\\
\vec{e_i}Ist der i-te Positionsvektor,\\
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}Ist zweidimensional x,y 2-Wege-Geschwindigkeitsvektor,\\
\omega = 1/\tau ， \Tau ist ein einmaliger Entspannungsfaktor.
Berechnen Sie dies für jedes $ i $, jeden Punkt.
Zuerst $ \ vec e_i $, der Positionsvektor zu Ihrem eigenen Quadrat + 8 benachbarten Quadraten, wenn Ihr Quadrat der Ursprung ist. $ i = 0 $ ist deine Masse $ \ vec e_0 = (0, 0) $, $ 1 \ le i \ le 4 $ ist oben, unten, links und rechts $ \ vec e_ {1 ... 4} = (0,1), (0, -1), (1,0), (-1,0) $, $ 5 \ le i \ le 8 $ steht neben dem Namen $ \ vec e_ {5 ... 8} = (1,1), (-1, -1), (1, -1), ( -1,1) $ Repräsentiert. Daher kann die $ (1) $ -Formel so interpretiert werden, dass sie den Zufluss aus der benachbarten Masse darstellt.
Als nächstes kommt $ \ vec e_i \ cdot \ vec u $.
Da $ \ vec u = (u_x, u_y) $ zum Beispiel, wenn $ i = 7 $,
Schließlich
Hindernisbehandlung
Es ist nicht sehr interessant, nur zu fließen, und der Kalman-Wirbel des Hauptmotivs tritt nicht auf. Denken wir also über die Verarbeitung nach, wenn ein Hindernis installiert ist. Dies ist jedoch mit der Gitter-Boltzmann-Methode einfach (ich denke, es war ein wenig mühsam, als ich die Navier-Stokes-Gleichung approximierte). Ursprünglich wird es für jede Zuflussrichtung separat berechnet. Wenn es also auf ein Hindernis trifft, kann es zu der Komponente geleitet werden, die der Zuflussrichtung entgegengesetzt ist.
Lassen Sie uns nun mit der Implementierung fortfahren.
Implementierung
Funktionsdefinition
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)Gleichung + Kollisionsverarbeitung mit Hindernissen
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)Formel
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]] #Der kollidierende Teil wird in entgegengesetzter Richtung an die Komponente gesendet.
def collide(): #(2)Gleichung + Zufluss vom linken Ende
global n, u, u9
rho = n.sum(axis=2).reshape(H, W, 1) #Berechnung von ρ
u = (n @ e) / rho #(ux,uy)Berechnung von
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)Formel
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) #Zufluss vom linken Ende
def curl(): #Farbbezogene Funktionen beim Plotten
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): #Ein Kind, das die Entsprechung des Sprungziels zum Zeitpunkt der Kollision berechnet, wenn es ein Array passiert, das Hindernisse darstellt
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])
Parameterdefinition und Initialisierung
Lattice_Boltzmann_Method.py
H, W = 40, 100 #Vertikal und horizontal
viscosity = .005 #Viskosität
om = 1 / (3 * viscosity + 0.5) #Diese Berechnung folgt der Referenzstelle. Es scheint der relationale Ausdruck zwischen Viskosität und ω zu sein.
u0 = [0, .12] #Anfangsgeschwindigkeit und Zuflussgeschwindigkeitsvektor vom linken Ende(-uy,ux)Entspricht
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)
Hindernis
Lattice_Boltzmann_Method.py
r = 3 #Die Länge der Wand oder der Radius des Kreises
#Wand
wall = np.zeros((H, W), bool)
wall[H // 2 - r:H // 2 + r, H // 2] = True
#Kreis
x, y = np.meshgrid(np.arange(W), np.arange(H))
circle = ((x - H / 2)**2 + (y - H / 2)**2 <= r**2)
mk_barrier(circle)
Handlung
Lattice_Boltzmann_Method.py
def update(i):
for t in range(10):
stream()
collide()
im.set_array(curl() - barrier[:, :, 0])
return im,
"""Werden Sie interaktiv für 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()
Ergebnis
Gleich wie am Anfang gezeigt,
Zusammenfassung
Sie werden beeindruckt sein, wenn Sie das natürliche Phänomen mit dem selbst geschriebenen Code reproduzieren können. Bitte probieren Sie es aus.
Alle Codes Tunagetano
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)Gleichung + Kollisionsverarbeitung mit Hindernissen
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)Formel
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]] #Der kollidierende Teil wird in entgegengesetzter Richtung an die Komponente gesendet.
def collide(): #(2)Gleichung + Zufluss vom linken Ende
global n, u, u9
rho = n.sum(axis=2).reshape(H, W, 1) #Berechnung von ρ
u = (n @ e) / rho #(ux,uy)Berechnung von
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)Formel
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) #Zufluss vom linken Ende
def curl(): #Farbbezogene Funktionen beim Plotten
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): #Ein Kind, das die Entsprechung des Sprungziels zum Zeitpunkt der Kollision berechnet, wenn es ein Array passiert, das Hindernisse darstellt
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 #Vertikal und horizontal
viscosity = .005 #Viskosität
om = 1 / (3 * viscosity + 0.5) #Diese Berechnung folgt der Referenzstelle. Es scheint der relationale Ausdruck zwischen Viskosität und ω zu sein.
u0 = [0, .12] #Anfangsgeschwindigkeit und Zuflussgeschwindigkeitsvektor vom linken Ende(-uy,ux)Entspricht
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 #Die Länge der Wand oder der Radius des Kreises
#Wand
wall = np.zeros((H, W), bool)
wall[H // 2 - r:H // 2 + r, H // 2] = True
#Kreis
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,
"""Werden Sie interaktiv für 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()
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)Gleichung + Kollisionsverarbeitung mit Hindernissen
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)Formel
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]] #Der kollidierende Teil wird in entgegengesetzter Richtung an die Komponente gesendet.
def collide(): #(2)Gleichung + Zufluss vom linken Ende
global n, u, u9
rho = n.sum(axis=2).reshape(H, W, 1) #Berechnung von ρ
u = (n @ e) / rho #(ux,uy)Berechnung von
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)Formel
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) #Zufluss vom linken Ende
def curl(): #Farbbezogene Funktionen beim Plotten
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): #Ein Kind, das die Entsprechung des Sprungziels zum Zeitpunkt der Kollision berechnet, wenn es ein Array passiert, das Hindernisse darstellt
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 #Vertikal und horizontal
viscosity = .005 #Viskosität
om = 1 / (3 * viscosity + 0.5) #Diese Berechnung folgt der Referenzstelle. Es scheint der relationale Ausdruck zwischen Viskosität und ω zu sein.
u0 = [0, .12] #Anfangsgeschwindigkeit und Zuflussgeschwindigkeitsvektor vom linken Ende(-uy,ux)Entspricht
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 #Die Länge der Wand oder der Radius des Kreises
#Wand
wall = np.zeros((H, W), bool)
wall[H // 2 - r:H // 2 + r, H // 2] = True
#Kreis
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,
"""Werden Sie interaktiv für 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()