在https://blog.csdn.net/weixin_58913471/article/details/117561995?spm=1001.2101.3001.6650.1&utm_medium=distribute.pc_relevant.none-task-blog-2~default~CTRLIST~default-1.no_search_link&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2~default~CTRLIST~default-1.no_search_link 在理解的基础上，重写，补充。
from numpy import * import matplotlib.pyplot as plt from matplotlib import cm R e = u L ν Re = u\frac{L}{\nu} Re=uνL 流体速度u，特征长度 L L L， ν \nu ν流体粘度 低速、高粘度和封闭流体条件导致流体条件低Re，如果粘性力占主导地位。Re<<这种流体被称为Stokes或者Creeping flow. 由于孔径小，这种流动在许多孔介质中的液体中很常见。 如果 ω \omega ω非常小，流体只能慢慢收敛到其平衡：高粘度。 流体的粘度和松弛参数 ω \omega ω成反比。 ν = δ t c s 2 ( 1 ω ? 1 2 ) \nu=\delta t c_{s}^{2}\left(\frac{1}{\omega}-\frac{1}{2}\right) ν=δtc
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 t_i ti 用于补偿不同长度的速度
5 8 \ | / –>4<–7 / | \ 3 6
因为
ρ ( x , t ) = ∑ i = 0 8 f i i n ( x , t ) \rho(\boldsymbol{x}, t)=\sum_{i=0}^{8} f_{i}^{\mathrm{in}}(\boldsymbol{x}, t) ρ(x,t)=i=0∑8fiin(x,t)
u ( x , t ) = 1 ρ ( x , t ) δ x δ t ∑ i = 0 8 v i f i i n ( 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) u(x,t)=ρ(x,t)1δtδxi=0∑8vifiin(x,t)
定义
ρ 1 = f 0 + f 1 + f 2 ( unknown ) ρ 2 = f 3 + f 4 + f 5 (known) ρ 3 = f 6 + f 7 + f 8 ( known ) \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} ρ1=f0+f1+f2( unknown )ρ2=f3+f4+f5 (known) ρ3=f6+f7+f8( known )
同时
ρ = ρ 1 + ρ 2 + ρ 3 ρ u x = ρ 1 − ρ 3 \begin{aligned} &\rho=\rho_{1}+\rho_{2}+\rho_{3} \\ &\rho u_{x}=\rho_{1}-\rho_{3} \end{aligned} ρ=ρ1+ρ2+ρ3ρux=ρ1−ρ3
u x u_x ux速度 x x x方向的分量
因此
ρ = ρ 2 + 2 ρ 3 1 − u x \rho=\frac{\rho_{2}+2 \rho_{3}}{1-u_{x}} ρ=1−uxρ2+2ρ3
rho[0,:] = 1/(1-u[0,0,:])*(sum(fin[col2,0,:],axis=0)+2*sum(fin[col3,0,:],axis=0))
E ( i , ρ , u ) = ρ t i ( 1 + v i ⋅ u c s 2 + 1 2 c s 4 ( v i ⋅ u ) 2 − 1 2 c s 2 ∣ u ∣ 2 ) E(i, \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) E(i,ρ,u)=ρti(1+cs2vi⋅u+2cs41(vi⋅u)2−2cs21∣u∣2)
密度总是接近于他们的平衡状态。 首先将未知密度分布函数初始化为其平衡值。 随后，检查相反分布函数偏离平衡的程度，再加上这个值作为修正。
f 0 i n = E ( 0 , ρ , u ) + ( f 8 i n − E ( 8 , ρ , u ) ) f 1 i n = E ( 1 , ρ , u ) + ( f 7 i n − E ( 7 , ρ , u ) ) f 2 i n = E ( 2 , ρ , u ) + ( f 6 i n − E ( 6 , ρ , u ) ) \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} f0inf1inf2in=E(0,ρ,u)+(f8in−E(8,ρ,u))=E(1,ρ,u)+(f7in−E(7,ρ,u))=E(2,ρ,u)+(f6in−E(6,ρ,u))
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 = 0 8 f i i n ( x , t ) \rho(\boldsymbol{x}, t)=\sum_{i=0}^{8} f_{i}^{\mathrm{in}}(\boldsymbol{x}, t) ρ(x,t)=i=0∑8fiin(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 = c s 2 ρ p=c_{s}^{2} \rho p=cs2ρ 在D2Q9模型中 c s 2 = 1 3 δ x 2 δ t 2 c_{s}^{2}=\frac{1}{3} \frac{\delta x^{2}}{\delta t^{2}} cs2=31δt2δx2
6 3 0 \ | / 7<–4–>1 / | \ 8 5 2
v 0 v_0 v0(1,1), v 1 v_1 v1(1,0), v 2 v_2 v2(1,-1) v 3 v_3 v3(0,1), v 4 v_4 v4(0,0), v 5 v_5 v5(0,-1) v 6 v_6 v6(-1,1), v 7 v_7 v7(-1,0), v 8 v_8 v8(-1,-1)
u ( x , t ) = 1 ρ ( x , t ) δ x δ t ∑ i = 0 8 v i f i i n ( 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) u(x,t)=ρ(x,t)1δtδxi=0∑8vifiin(x,t)
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
E ( i , ρ , u ) = ρ t i ( 1 + v i ⋅ u c s 2 + 1 2 c s 4 ( v i ⋅ u ) 2 − 1 2 c s 2 ∣ u ∣ 2 ) E(i, \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) E(i,ρ,u)=ρti(1+cs2vi⋅u+2cs41(vi⋅u)2−2cs21∣u∣2)
# 平衡态计算
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 out − f i in = − ω ( f i in − E ( i , ρ , u ⃗ ) ) f_{i}^{\text {out }}-f_{i}^{\text {in }}=-\omega\left(f_{i}^{\text {in }}-E(i, \rho, \vec{u})\right) fiout −fiin =−ω(fiin −E(i,ρ,u ))
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
mid-grid bounce-back
f i in ( x , t + 1 ) = f j out ( x , t ) f_{i}^{\text {in }}(x, t+1)=f_{j}^{\text {out }}(x, t) fiin (x,t+1)=fjout (x,t)
v i = − v j v_{i}=-v_{j} vi=−vj
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()