Python-LBM(格子玻尔兹曼)程序源码实例分析—圆柱绕流篇

初次学习LBM计算方法，找到一个比较优秀的用python语言编写的圆柱绕流的实例，对每段代码详细添加了注释，帮助自己总结，也为初学者提供一点帮助。先放一张结果图像：

weixin_58913471

11466人浏览 · 2021-06-12 18:21:58

weixin_58913471 · 2021-06-12 18:21:58 发布

初次学习LBM计算方法，找到一个比较优秀的用python语言编写的圆柱绕流的实例，对每段代码详细添加了注释，帮助自己总结，也为初学的朋友们提供一点帮助（全部代码在文章最后）。

先放一张结果图像：

1、导入模块，只需要两个

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

2、定义LBM程序所需要的变量：

maxIter = 200000 #总计算次数
Re      = 100.0  #模型要模拟的雷诺数
nx = 520         #模型x轴长度
ny = 180         #模型y轴长度
ly=ny-1.0        #稍后用于计算入口速度，为模型增加扰动，因为当雷诺数较小时，计算缺少扰动。
q = 9            #本文选择的时D2Q9模型
cx = nx/4        #圆柱圆心x坐标
cy=ny/2          #圆柱圆心y坐标
r=ny/9           #圆柱半径长度
uLB = 0.04       #模型的入口速度
nulb=uLB*r/Re    #模型流体的粘性系数，（这里用的是半径，严格来说应该用直径）
omega = 1.0 / (3.*nulb+0.5)  #LBGK(单松弛模型)的弛豫时间

需要有一点注意的就是弛豫时间（有的说松弛时间，我也不确定）的计算，在LBM中模拟流动问题，只需要控制雷诺数相同即可。

雷诺数为：

$Re=\frac{UD}{v}$

在圆柱绕流中U表示入口平均速度，D表示圆柱的直径，v（代码中的nulb）表示流体粘性系数

在LBGK中流体粘性系数和松弛时间的关系则是：

$\tau =\frac{c^{2}}{3}\left ( v-0.5 \right )dt$

其中 $\tau$ 代表松弛时间，c代表格子速度（dx/dt），v代表粘性系数，dt表示时间步长

通过这两个公式就可以将速度、雷诺数、粘性系数和松弛时间联系起来，方便以后模型参数灵活调控。

3、定义分布函数、格子方向等变量

#c代表格子速度，形状为（9，2），9代表有9个方向，2分别代表x,y方向，具体参加示意图
c = array([(x,y) for x in [0,-1,1] for y in [0,-1,1]])

#t代表权重系数，N-S方程的权重系数有三个取值：1/36，1/9，4/9
t = 1./36. * ones(q)            
t[asarray([norm(ci)<1.1 for ci in c])] = 1./9.
t[0] = 4./9.

#noslip为格子速度的反方向
noslip = [c.tolist().index((-c[i]).tolist()) for i in range(q)] 

#i1为6，7，8方向，用于出口边界（右边界）边界处理
i1 = arange(q)[asarray([ci[0]<0  for ci in c])] # Unknown on right wall.

#i2为0，1，2方向，用于入口边界（左边界）边界处理
i2 = arange(q)[asarray([ci[0]==0 for ci in c])] # Vertical middle.

#i3为3，4，5方向，用于入口边界（左边界）边界处理
i3 = arange(q)[asarray([ci[0]>0  for ci in c])] # Unknown on left wall.

本算例中的格子方向的编号和一般情况有所不同：

其实根据c的取值可以判断其方向如c[1]取值为[0,-1]，其代表向量方向为[0,-1]，即沿着y轴负方向。

4、定义平衡态方程

碰撞和扩散对所有LBM问题基本都是相似的，只有平衡态方程决定着问题的属性，即问题的控制方程，N-S方程下的LBM平衡态函数表示为：

$f^{eq}=w_{i}\rho \left [ 1+\frac{c_{i}\cdot u}{c_{s}^{2}}+\frac{\left ( c_{i}\cdot u \right )^{2}}{2c_{s}^{2}}-\frac{u^{2}}{2c_{s}^{2}} \right ]$

使用代码表示为：

sumpop = lambda fin: sum(fin,axis=0) #定义计算宏观密度的函数，fin代表粒子分布函数

def equilibrium(rho,u):              #平衡态计算

    cu   = 3.0 * dot(c,u.transpose(1,0,2))  #

    usqr = 3./2.*(u[0]**2+u[1]**2)

    feq = zeros((q,nx,ny))

    for i in range(q): feq[i,:,:] = rho*t[i]*(1.+cu[i]+0.5*cu[i]**2-usqr)

    return feq

其中较为难理解的是cu的计算，dot表示numpy中的点积，c的形状为（9，2），而u的形状为（2，nx，ny），u.transpose(1,0,2)，记为up，表示将第一轴与第二轴互换，其形状变为（nx,2,ny）。

这是一个二维矩阵和三维矩阵的乘法，运算的结果为cu的形状为(9,nx,ny)

运算的规则比较复杂，可以结合格点的物理意义来理解：

首先c保持不动,从c[0]开始，取出up的第一个维度的第一列，就是up[0,:,0]

c[0]*up[0,:,0]等于cu[0,0,0]第一个维度第一行的第一个取值，也就是第一个格点0方向的取值

然后取出up的第一个维度的第二列，就是up[0,:,1]，c[0]*up[0,:,1]等于cu[0,0,1]第一个维度第一行的第二个取值，也就是第二个格点0方向的取值

依次类推，可以求得up的全部取值。虽然计算比较抽象，需要多理解、多运算，但是比其写循环来计算要快很多。

5、定义圆柱相关

obstacle = fromfunction(lambda x,y: (x-cx)**2+(y-cy)**2<r**2, (nx,ny))

#vel用于初始化速度，d的取值是0，1，代表x方向和y方向，ulb后面的代表微小偏量
#其实直接全部赋值0.04也没关系
vel = fromfunction(lambda d,x,y: (1-d)*uLB*(1.0+1e-4*sin(y/ly*2*pi)),(2,nx,ny))

#计算所有节点的平衡态
feq = equilibrium(1.0,vel)
#初始化分布函数
fin = feq.copy()

需要注意的是numpy.fromfunction（f,shape）这个函数的使用，有两个参数f函数是用于构建数组的函数，shape是数组输出的形状，其中f函数的输入和shape有关，具体用法就不写了。

obstacle是和fin形状一样的数组，其中圆柱外的节点是false，圆柱内的节点是true

vel是初始化速度的变量

6、开始计算

for time in range(maxIter):

    #出口边界条件，充分发展的流动格式即出口未知的分布函数与上一个节点相同
    fin[i1,-1,:] = fin[i1,-2,:] 
    
    rho = sumpop(fin)           #计算宏观密度
    
    u = dot(c.transpose(), fin.transpose((1,0,2)))/rho  #计算宏观速度

    #入口边界条件，zou-he速度边界格式
    u[:,0,:] =vel[:,0,:] 
    rho[0,:] = 1./(1.-u[0,0,:]) * (sumpop(fin[i2,0,:])+2.*sumpop(fin[i1,0,:]))
    
    #计算平衡态
    feq = equilibrium(rho,u)
    
    #计算入口边界的6，7，8方向的分布函数
    fin[i3,0,:] = fin[i1,0,:] + feq[i3,0,:] - fin[i1,0,:]
    
    #碰撞过程
    fout = fin - omega * (fin - feq)
    
    #对圆柱内的节点应用反弹格式
    for i in range(q): fout[i,obstacle] = fin[noslip[i],obstacle]
    
    #扩散过程
    for i in range(q):

        fin[i,:,:] = roll(roll(fout[i,:,:],c[i,0],axis=0),c[i,1],axis=1)
    
    #绘图
    if (time%10000==0):

        plt.clf(); plt.imshow(sqrt(u[0]**2+u[1]**2).transpose(),cmap=cm.Reds)

        plt.savefig("vel."+str(time/10000).zfill(4)+".png")

程序的重头戏，在LBM程序设计中，原则上只要将流场中所有节点的值都更新一遍就可以，没有严格的先后顺序，但是一般大家都按照先碰撞、再扩散、再边界处理的方式。不过本例子就是灵活设计的，下面详细讲解一下：

根据LGBK的控制方程：

粒子在当前时间步当前节处完成一次碰撞后，8个方向的粒子会继续沿着8个方向到达下一节点（0方向不变），如节点（33，33）处6方向（沿X轴正方向）的粒子会到达（33，34）处，

即

$f_{i}\left ( x+e_{i}\delta t,t \right )=f_{i}\left ( x,t \right )$

当所有节点都完成迁移运动之后，LBM一次计算就完成了。

在本例中使用的是numpy.roll（a,shift,axis）函数实现迁移功能，这个函数的作用将数组a就是沿axis轴（0表示沿着列（y方向），1表示沿着行（x方向））向前滚动，滚动的间隔是shift。

也就是说axis=0时，将第0行变为第1行，第1行变为第2行，最后一行变为第0行，当axis=1时，将第0列变为第1列，第1列变为第2列，最后一列变为第0列

在本例中，roll运算的过程是fout先沿着asix=0向前滚动一行，此时相当于其x坐标加shift，滚动的值应该是c[i,0]。此时若i=6，7，8,c为正值，x坐标加1,若c为3，4，5,c为负值，x坐标减1.其他方向类似，第二层roll也类似。

这样计算对上下边界是没问题的，因为上边界出去的粒子（2，5，8方向）滚动以后变成了下边界进入的粒子（1，4，7方向），这实际上是LBM中的周期性边界条件。这种边界条件认为上下边界是开放的，对与圆柱绕流计算是适用的（也有认为上下边界是固体边界，那就需要使用反弹边界条件）

但是对左右边界存在一个问题，我们知道左右边界是速度边界条件，就是左边界的格点的分布函数是我们通过zou-he边界条件计算出来的。所以我们就需要在下一次循环刚开始的时候就重新计算左右边界的值

还有关于圆柱边界的处理，本例中采用的是直接反弹格式，我们知道，在扩散过程中会有流体节点的粒子以i方向进入固体节点，同时会有固体节点的粒子以-i的方向进入流体节点，反弹格式就是将进入固体节点的粒子方向进行反转，反转后该方向的粒子就会在下一个扩散计算中重新进入流体节点中。

本例中的实现方式就是fout[i]=fin[noslip[i]]

所以我们的一次循环的计算过程就是：

1-计算右边界分布函数，
2-计算宏观密度，宏观速度，
注意此时左边界处的密度、速度是不对的，因为左边界的分布函数是右边界迁移过来的，但是没关系因为，因为左边界的密度（1.0）、速度（0.04）是已知量，所以，
3-重新计算左边界分布函数，
这个时候我们已经把所有节点都进行了一次运算，然后，
4-对所有节点计算碰撞。
5-处理圆柱边界
6-扩散计算
7-可视化

7、模型验证

通过阻力系数和升力系数可以验证计算结果，其计算公式为：

在这里纠结了很久，主要原因是fx的计算，有些文献中的公式是：

这里的Fd表示总的曳力，Fx表示x方向的曳力，两者正好相差一个A（二维情况是圆的直径），在LBM计算中fx比较容易取得，所以我们采用的是第一个公式。

而Fx的计算是采用的动量交换法，具体的计算公式是：

$F_{x}=\sum e_{i}(fout(s,t)+fout(w,t))$

其中s是固体内的节点，w是与s相邻的流体的节点。

python代码为：

fx=0
fy=0

for i in range(q):

    for x in range(100,170,1): 

        for y in range(50,130,1):

            if obstacle[x,y]:

                x1=x+c[i,0]
                y1=y+c[i,1]

                if not obstacle[x1,y1]:

                    fx+=c[i,0]*(fout[i,x,y]+fout[noslip[i],x1,y1])
                    fy+=c[i,1]*(fout[i,x,y]+fout[noslip[i],x1,y1])


fx_1=abs(2*fx/(r*0.04**2))
fy_1=abs(2*fy/(r*0.04**2))

print(fx_1)
print(fy_1)

经过20万次的计算，计算耗时2小时左右，cd的值在3.2左右波动，符合文献中的计算结果。至此，圆柱绕流的计算圆满成功！附上所有代码：

# 2D flow around a cylinder
from numpy import *; from numpy.linalg import *
import matplotlib.pyplot as plt; from matplotlib import cm
###### Flow definition #########################################################
maxIter = 200000 # Total number of time iterations.
Re      = 220.0  # Reynolds number.
nx = 520; ny = 180; ly=ny-1.0; q = 9 # Lattice dimensions and populations.
cx = nx/4; cy=ny/2; r=ny/9;          # Coordinates of the cylinder.
uLB     = 0.04                       # Velocity in lattice units.
nulb    = uLB*r/Re; omega = 1.0 / (3.*nulb+0.5); # Relaxation parameter.

###### Lattice Constants #######################################################
c = array([(x,y) for x in [0,-1,1] for y in [0,-1,1]]) # Lattice velocities.
t = 1./36. * ones(q)                                   # Lattice weights.
t[asarray([norm(ci)<1.1 for ci in c])] = 1./9.; t[0] = 4./9.
noslip = [c.tolist().index((-c[i]).tolist()) for i in range(q)] 
i1 = arange(q)[asarray([ci[0]<0  for ci in c])] # Unknown on right wall.
i2 = arange(q)[asarray([ci[0]==0 for ci in c])] # Vertical middle.
i3 = arange(q)[asarray([ci[0]>0  for ci in c])] # Unknown on left wall.

###### Function Definitions ####################################################
sumpop = lambda fin: sum(fin,axis=0) # Helper function for density computation.
def equilibrium(rho,u):              # Equilibrium distribution function.
    cu   = 3.0 * dot(c,u.transpose(1,0,2))
    usqr = 3./2.*(u[0]**2+u[1]**2)
    feq = zeros((q,nx,ny))
    for i in range(q): feq[i,:,:] = rho*t[i]*(1.+cu[i]+0.5*cu[i]**2-usqr)
    return feq

###### Setup: cylindrical obstacle and velocity inlet with perturbation ########
obstacle = fromfunction(lambda x,y: (x-cx)**2+(y-cy)**2<r**2, (nx,ny))

#vel用于初始化速度，d的取值是0，1，代表x方向和y方向，ulb后面的代表微小偏量
#其实直接全部赋值0.04也没关系
#vel = fromfunction(lambda d,x,y: (1-d)*uLB*(1.0+1e-4*sin(y/ly*2*pi)),(2,nx,ny))
vel = fromfunction(lambda d,x,y: (1-d)*uLB,(2,nx,ny))

feq = equilibrium(1.0,vel); fin = feq.copy()

w=zeros((9,nx,ny))
w[:,obstacle]=1


###### Main time loop ##########################################################
for time in range(maxIter):
    
    fin[i1,-1,:] = fin[i1,-2,:] #出口边界条件，充分发展的流动格式即出口未知的分布函数与上一个节点相同
    
    rho = sumpop(fin)           #计算宏观密度
    
    u = dot(c.transpose(), fin.transpose((1,0,2)))/rho  #

    u[:,0,:] =vel[:,0,:] # 
    rho[0,:] = 1./(1.-u[0,0,:]) * (sumpop(fin[i2,0,:])+2.*sumpop(fin[i1,0,:]))

    feq = equilibrium(rho,u) # 
    fin[i3,0,:] = fin[i1,0,:] + feq[i3,0,:] - fin[i1,0,:]
    fout = fin - omega * (fin - feq)  # 碰撞过程
    
    for i in range(q): fout[i,obstacle] = fin[noslip[i],obstacle]
    
    for i in range(q): #扩散过程
        
        fin[i,:,:] = roll(roll(fout[i,:,:],c[i,0],axis=0),c[i,1],axis=1)
    
    if (time%2000==0): # Visualization

        fd=0
        ff=0

        for i in range(q):

            for x in range(100,170,1): 

                for y in range(50,130,1):

                    if obstacle[x,y]:

                        x1=x+c[i,0]
                        y1=y+c[i,1]

                        if not obstacle[x1,y1]:

                            fd+=c[i,0]*(fout[i,x,y]+fout[noslip[i],x1,y1])
                            ff+=c[i,1]*(fout[i,x,y]+fout[noslip[i],x1,y1])


        fd_1=abs(2*fd/(r*0.04**2))
        ff_1=abs(2*ff/(r*0.04**2))

        print(fd_1)
        print(ff_1)
            
    plt.clf(); plt.imshow(sqrt(u[0]**2+u[1]**2).transpose(),cmap=cm.Reds)
    plt.savefig("vel."+str(time/10000).zfill(4)+".png")