LBM学习记录3 Python实现D2Q9圆柱绕流

在https://blog.csdn.net/weixin_58913471/article/details/117561995?spm=1001.2101.3001.6650.1&utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7Edefault-1.no_search_link&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7Edefault-1.no_search_link 基础上进行理解，重写，补充。

from numpy import *
import matplotlib.pyplot as plt
from matplotlib import cm

雷诺数（Re）是反应粘性和惯性之间平衡的无量纲数。
$u\frac{L}{\nu}$
流体速度u，特征长度 $L$ ， $ν\nu$ 流体粘度
低速、高粘度和封闭的流体条件导致低Re，粘性力占主导地位。如果Re<<1，这种流体被称为Stokes或者Creeping flow. 由于孔径小，这种流动在许多多孔介质中的液体中是很常见的。
omega 松弛频率
如果 $ω\omega$ 很小，流体只能缓慢地收敛到其平衡状态：高粘性。
流体的粘性与松弛参数 $ω\omega$ 成反比。
$\nu=\delta t c_{s}^{2}\left(\frac{1}{\omega}-\frac{1}{2}\right)$
其中， $cs2=13c_{s}^{2}=\frac{1}{3}$

maxIter = 200000  # 迭代次数
Re = 220.0         # 雷诺数
uLB = 0.04        # 模型的入口速度
nx, ny = 420, 180 # x, y轴长度
ly = ny-1         # 用于计算入口，为模型增加扰动，当雷诺数较小时，计算缺少扰动
cx, cy, r = nx//4, ny//2, ny//9  # 圆柱坐标
nulb = uLB*r/Re   # 粘性系数
omega = 1/(3*nulb+0.5)  # 松弛频率

$t_i$ 用于补偿不同长度的速度

流入条件

密度

2 5 8
\ | /
1–>4<–7
/ | \
0 3 6

因为

$ρ(x,t)=∑i=08fiin(x,t)\rho(\boldsymbol{x}, t)=\sum_{i=0}^{8} f_{i}^{\mathrm{in}}(\boldsymbol{x}, t)$

$u(x,t)=1ρ(x,t)δxδt∑i=08vifiin(x,t)\boldsymbol{u}(\boldsymbol{x}, t)=\frac{1}{\rho(\boldsymbol{x}, t)} \frac{\delta x}{\delta t} \sum_{i=0}^{8} \boldsymbol{v}_{i} f_{i}^{\mathrm{in}}(\boldsymbol{x}, t)$

定义

$\begin{aligned} &\rho_{1}=f_{0}+f_{1}+f_{2}(\text { unknown }) \\ &\rho_{2}=f_{3}+f_{4}+f_{5} \text { (known) } \\ &\rho_{3}=f_{6}+f_{7}+f_{8}(\text { known }) \end{aligned}$

同时

$\begin{aligned} &\rho=\rho_{1}+\rho_{2}+\rho_{3} \\ &\rho u_{x}=\rho_{1}-\rho_{3} \end{aligned}$

$u_x$ 速度 $x$ 方向的分量

因此

$ρ=ρ2+2ρ31−ux\rho=\frac{\rho_{2}+2 \rho_{3}}{1-u_{x}}$

rho[0,:] = 1/(1-u[0,0,:])*(sum(fin[col2,0,:],axis=0)+2*sum(fin[col3,0,:],axis=0))

流入侧f0 f1 f2的密度分布函数

$\rho, \boldsymbol{u})=\rho t_{i}\left(1+\frac{\boldsymbol{v}_{\boldsymbol{i}} \cdot \boldsymbol{u}}{c_{s}^{2}}+\frac{1}{2 c_{s}^{4}}\left(\boldsymbol{v}_{\boldsymbol{i}} \cdot \boldsymbol{u}\right)^{2}-\frac{1}{2 c_{s}^{2}}|\boldsymbol{u}|^{2}\right)$

密度总是接近于他们的平衡状态。
首先将未知密度分布函数初始化为其平衡值。
随后，检查相反分布函数偏离平衡的程度，再加上这个值作为修正。

$\begin{aligned} f_{0}^{\mathrm{in}} &=E(0, \rho, \boldsymbol{u})+\left(f_{8}^{\mathrm{in}}-E(8, \rho, \boldsymbol{u})\right) \\ f_{1}^{\mathrm{in}} &=E(1, \rho, \boldsymbol{u})+\left(f_{7}^{\mathrm{in}}-E(7, \rho, \boldsymbol{u})\right) \\ f_{2}^{\mathrm{in}} &=E(2, \rho, \boldsymbol{u})+\left(f_{6}^{\mathrm{in}}-E(6, \rho, \boldsymbol{u})\right) \end{aligned}$

fin[[0,1,2],0,:] = feq[[0,1,2],0,:] + fin[[8,7,6],0,:] - feq[[8,7,6],0,:]

v = array([[1,1], [1,0], [1,-1],
           [0,1], [0,0], [0,-1],
           [-1,1], [-1,0], [-1,-1]])

t = array([1/36, 1/9, 1/36,
          1/9, 4/9, 1/9,
          1/36, 1/9, 1/36])

col1 = array([0, 1, 2])
col2 = array([3, 4, 5])
col3 = array([6, 7, 8])

密度

$ρ(x,t)=∑i=08fiin(x,t)\rho(\boldsymbol{x}, t)=\sum_{i=0}^{8} f_{i}^{\mathrm{in}}(\boldsymbol{x}, t)$

rho = zeros((nx, ny))
for ix in range(nx):
    for iy in range(ny):
        rho[ix, iy] = 0
        for i in range(9):
            rho[ix, iy] += fin[i, ix, iy]

压力
压力正比于密度，根据理想气体状态方程，在等温气体中
$p=cs2ρp=c_{s}^{2} \rho$
在D2Q9模型中
$cs2=13δx2δt2c_{s}^{2}=\frac{1}{3} \frac{\delta x^{2}}{\delta t^{2}}$

速度

6 3 0
\ | /
7<–4–>1
/ | \
8 5 2

$v_0$ (1,1), $v_1$ (1,0), $v_2$ (1,-1)
$v_3$ (0,1), $v_4$ (0,0), $v_5$ (0,-1)
$v_6$ (-1,1), $v_7$ (-1,0), $v_8$ (-1,-1)

v = np.array([[1,1], [1,0], [1,-1], [0,1], [0,0], [0,-1], [-1,1], [-1,0], [-1,-1]])
u = zeros((2, nx, ny))
for ix in range(nx):
    for iy in range(ny):
        u[0, ix, iy] = 0
        u[1, ix, iy] = 0
        for i in range(9):
            u[0,ix,iy] += v[i,0]*fin[i,ix,iy]
            u[1,ix,iy] += v[i,1]*fin[i,ix,iy]

def macroscopic(fin):
    rho = sum(fin, axis=0)
    u = zeros((2, nx, ny))
    for i in range(9):
        u[0,:,:] += v[i,0] * fin[i,:,:]
        u[1,:,:] += v[i,1] * fin[i,:,:]
    u /= rho
    return rho, u

平衡方程
$\rho, \boldsymbol{u})=\rho t_{i}\left(1+\frac{\boldsymbol{v}_{\boldsymbol{i}} \cdot \boldsymbol{u}}{c_{s}^{2}}+\frac{1}{2 c_{s}^{4}}\left(\boldsymbol{v}_{\boldsymbol{i}} \cdot \boldsymbol{u}\right)^{2}-\frac{1}{2 c_{s}^{2}}|\boldsymbol{u}|^{2}\right)$

# 平衡态计算
def equilibrium(rho, u):
    usqr = 3/2 * (u[0]**2 + u[1]**2)
    feq = zeros((9, nx, ny))
    for i in range(9):
        cu = 3 * (v[i,0]*u[0,:,:] + v[i,1]*u[1,:,:])
        feq[i,:,:] = rho*t[i] * (1 + cu + 0.5*cu**2 - usqr)
    return feq

碰撞

$f_{i}^{\text {out }}-f_{i}^{\text {in }}=-\omega\left(f_{i}^{\text {in }}-E(i, \rho, \vec{u})\right)$

fout = fin-omega*(fin-eq)

迁移

for ix in range(nx):
    for iy in range(ny):
        for i in range(9):
            next_x = ix + v[i,0]
            if next_x<0:
                next_x = nx-1
            if next_x>=nx:
                next_y = 0
            
            next_y = iy + v[i,1]
            if next_y<0:
                next_y = ny-1
            if next_y>=ny:
                next_y = 0    
            
            fin[i,next_x,next_y] = fout[i,ix,iy]

np.roll(a,shift,axis=None) 将数组a，沿着axis方向，滚动shift长度
可改写成

for i in range(9):
    fin[i, :, :] = roll(roll(fout[i,:,:],v[i,0], axis=0), v[i,1], axis=1)

边界条件
Bounce-back BCs
反弹边界条件是处理静止无滑移壁面的一类常用格式，是指当离散分布函数到达边界节点时，将沿着其进入的方向散射回流体，包括on-grid bounce-back（边界与晶格点对齐）和 mid-grid bounce-back（边界在界外节点和界内节点之间的中心）。

on-grid bounce-back

on-grid

mid-grid bounce-back

mid-grid

$(x,t)f_{i}^{\text {in }}(x, t+1)=f_{j}^{\text {out }}(x, t)$

$v_{i}=-v_{j}$

for i in range(9):
    fout[i,obstacle] = fin[8-i, obstacle]

def obstacle_fun(x, y):
    return (x-cx)**2 + (y-cy)**2 < r**2

fromfunction从函数中创建数组，返回数组，符合条件值为True，不符合为False。

obstacle = fromfunction(obstacle_fun, (nx, ny))

# 初始速度曲线：几乎为零，有一个轻微的扰动来触发不稳定。
def inivel(d, x, y):
    return (1-d) * uLB * (1+1e-4*sin(y/ly*2*pi))

vel = fromfunction(inivel, (2, nx, ny))

# 以给定的速度初始化处于平衡状态的种群
fin = equilibrium(1, vel)

for time in range(maxIter):
    # 右边界分布函数
    fin[col3, -1, :] = fin[col3, -2, :]
    
    #计算宏观密度和速度
    rho, u = macroscopic(fin)
    
    # 重新计算左边界的分布函数
    u[:, 0, :] = vel[:, 0, :]
    rho[0,:] = 1/(1-u[0,0,:]) * (sum(fin[col2,0,:], axis=0) + 2*sum(fin[col3,0,:], axis=0))
    
    # 计算平衡态
    feq = equilibrium(rho, u)
    fin[[0,1,2],0,:] = feq[[0,1,2],0,:]+fin[[8,7,6],0,:]-feq[[8,7,6],0,:]
    
    # 碰撞过程
    fout = fin - omega * (fin - feq)
    
    # 对圆柱内节点进行反弹
    for i in range(9):
        fout[i, obstacle] = fin[8-i, obstacle]
    
    # 扩散过程
    for i in range(9):
        fin[i, :, :] = roll(roll(fout[i,:,:],v[i,0], axis=0), v[i,1], axis=1)
        
    # 可视化
    if (time%100==0):
        plt.clf()
        plt.imshow(sqrt(u[0]**2+u[1]**2).transpose(), cmap=cm.Reds)
        plt.savefig("/share/home/yszhang/PBS/LBM/pic/"+str(time//100)+".jpg")

import cv2

file_dir = '/pic/'
video = cv2.VideoWriter('video.avi',cv2.VideoWriter_fourcc(*'MJPG'),5,(1280,720))

for i in range(2000):
    img = cv2.imread(file_dir+str(i)+'.jpg')     
    img = cv2.resize(img,(1280,720))
    video.write(img)

video.release()