# Introduction

This time, I would like to do a fluid simulation like the above gif using the lattice Boltzmann method.

** What is the lattice Boltzmann method **
> The lattice Boltzmann method (Lattice Boltzmann method) approximates a fluid with an aggregate of a large number of virtual particles (lattice gas model), and sequentially develops the collision and translation of each particle using the velocity distribution function of the particles. This is a simulation method for finding the heat flow field of a fluid by calculating the moment. -Wikipedia

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.

# Method

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$, \vec e_7 \cdot \vec u = (1,-1)\cdot (u_x,u_y) = u_x-u_y It will be.

Finally|u|^2 、 This is the square of the length of $\ vec u$, so |u|^2=u_x^2+u_y^2is.

# Obstacle handling

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.

# Implementation

### Function definition

#### 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])


### Parameter definition and initialization

#### 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)


### Obstacle

#### 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)


### plot

#### 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()


# result

Same as the one shown at the beginning,

# Summary

You'll be impressed if you can reproduce natural phenomena with your own code. Please give it a try.

All chords 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)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()

## reference [Fluid Dynamics Simulation](http://physics.weber.edu/schroeder/fluids/)