[Python] Let's observe the Karman vortex by the lattice Boltzmann method. [Fluid simulation]

Introduction

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.

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.

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