This time, I would like to do a fluid simulation like the above gif using the lattice Boltzmann method.
Also, when I wrote it as short as possible this time, it was slower than the sample code on the reference site (link at the end of the article). Even Numpy, which has a fast matrix calculation, doesn't seem to rely too much on it.
The detailed theory is omitted here, and the contents required for implementation are explained.
The formula to calculate is shown first.
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}
here,\\
n_i is the i component of the particle distribution function,\\
\vec{e_i}Is the i-th position vector,\\
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}Is two-dimensional x,y Two-way velocity vector,\\
\omega = 1/\tau , \tau is a single time relaxation factor.
Calculate this for each $ i $, each point.
First, $ \ vec e_i $, which is the position vector toward your own square + the adjacent 8 squares when your square is the origin. $ i = 0 $ is your cell $ \ vec e_0 = (0, 0) $, $ 1 \ le i \ le 4 $ is up / down / left / right $ \ vec e_ {1 ... 4} = (0,1), (0, -1), (1,0), (-1,0) $, $ 5 \ le i \ le 8 $ is adjacent to the name $ \ vec e_ {5 ... 8} = (1,1), (-1, -1), (1, -1), ( -1,1) $ Represents. Therefore, the $ (1) $ formula can be interpreted as representing the inflow from the adjacent square.
Next is $ \ vec e_i \ cdot \ vec u $.
Since $ \ vec u = (u_x, u_y) $, for example, if $ i = 7 $,
Finally
It's not very interesting just to flow, and the Karman vortex of the main subject does not occur, so let's think about the processing when an obstacle is installed. However, this is easy with the Lattice Boltzmann method (I think it was a little troublesome when approximating the Navier-Stokes equation). Originally, the calculation is done separately for each inflow direction, so if you hit an obstacle, you can just flow it to the component opposite to the inflow direction.
Now let's go with the implementation.
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)Formula + collision processing with 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)formula
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]] #The part that collided is sent to the component in the opposite direction.
def collide(): #(2)Equation + inflow from the left end
global n, u, u9
rho = n.sum(axis=2).reshape(H, W, 1) #Calculation of ρ
u = (n @ e) / rho #(ux,uy)Calculation
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)formula
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) #Inflow from the left end
def curl(): #Functions related to color when plotting
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): #A child who calculates the correspondence of the bounce destination at the time of collision when passing an array expressing 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])
Lattice_Boltzmann_Method.py
H, W = 40, 100 #Vertical and horizontal
viscosity = .005 #viscosity
om = 1 / (3 * viscosity + 0.5) #This calculation follows the reference site. It seems to be the relational expression between viscosity and ω.
u0 = [0, .12] #Initial velocity and inflow velocity vector from the left end(-uy,ux)Corresponds to
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)
Lattice_Boltzmann_Method.py
r = 3 #The length of the wall, the radius of the circle,
#wall
wall = np.zeros((H, W), bool)
wall[H // 2 - r:H // 2 + r, H // 2] = True
#Circle
x, y = np.meshgrid(np.arange(W), np.arange(H))
circle = ((x - H / 2)**2 + (y - H / 2)**2 <= r**2)
mk_barrier(circle)
Lattice_Boltzmann_Method.py
def update(i):
for t in range(10):
stream()
collide()
im.set_array(curl() - barrier[:, :, 0])
return im,
"""Interactive for 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()
Same as the one shown at the beginning,
You'll be impressed if you can reproduce natural phenomena with your own code. Please give it a try.
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)Formula + collision processing with 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)formula
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]] #The part that collided is sent to the component in the opposite direction.
def collide(): #(2)Equation + inflow from the left end
global n, u, u9
rho = n.sum(axis=2).reshape(H, W, 1) #Calculation of ρ
u = (n @ e) / rho #(ux,uy)Calculation
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)formula
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) #Inflow from the left end
def curl(): #Functions related to color when plotting
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): #A child who calculates the correspondence of the bounce destination at the time of collision when passing an array expressing 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 and horizontal
viscosity = .005 #viscosity
om = 1 / (3 * viscosity + 0.5) #This calculation follows the reference site. It seems to be the relational expression between viscosity and ω.
u0 = [0, .12] #Initial velocity and inflow velocity vector from the left end(-uy,ux)Corresponds to
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 #The length of the wall, the radius of the circle,
#wall
wall = np.zeros((H, W), bool)
wall[H // 2 - r:H // 2 + r, H // 2] = True
#Circle
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,
"""Interactive for 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()