河南科技
Henan Science and Technology 信息技术
总783期第十三期
2022年7月
基于格子Boltzmann方法的泊肃叶流数值模拟
研究
程海龙
（河南省机场集团有限公司，河南郑州450000）
摘要：本研究采用格子Boltzmann方法分别对压力驱动和速度驱动的二维泊肃叶不可压缩流体的流动进行数值模拟，模型采用D2Q9模型及标准的碰撞迁移规则，压力驱动时采用的是含有力项的Boltzmann-BGK方程，边界采用非平衡态外推格式和充分发展格式，基于格子Boltzmann方法对流体流动的速度、压力、速度矢量进行分析对比。模拟试验结果表明，当流体流动达到稳定后，沿着平板方向截面上的速
度大小呈现抛物线分布，其中心线上的压强沿着平板方向线性减小，在速度驱动下，某截面速度的模拟解与解析解最大误差仅为2.01%，从而得出在较小压力驱动下，边界层和驱动力对流体产生双重影响。
关键词：格子Boltzmann；数值模拟；D2Q9；泊肃叶流
中图分类号：TG333文献标志码：A文章编号：1003-5168（2022）13-0019-04 DOI:10.19968/jki.hnkj.1003-5168.2022.13.004
Numerical Simulation of Poiuille Flow Bad on Lattice Boltzmann
Method
CHENG Hailong
(Henan Airport Group Co.,Ltd.,Zhengzhou450000,China)
Abstract:In this study,the lattice Boltzmann method is ud to numerically simulate the pressure-driven and velocity-driven two-dimensional Poiuille incompressible fluid flow,respectively.The D2Q9model is ud,the standard collision migration is ud,the Boltzmann-BGK equation with force
term is ud for pressure driving,and the boundary adopts the non-equilibrium extrapolation scheme and the fully devel⁃oped scheme are ud to analyze and compare the flow velocity,pressure and velocity vector bad on the lattice Boltzmann method.The simulation results show that:after the flow is stable,the velocity on the c⁃tion along the direction of the plate prents a parabolic distribution,and it is concluded that the pressure on the center line decreas linearly along the direction of the plate.The analytical solution is in clo agreement,and the maximum error with the analytical solution is only2.01%.It is concluded that under the driving of small pressure,the boundary layer and driving force have dual effects on the fluid. Keywords:lattice Boltzmann;numerical simulation;D2Q9;poiuille
0引言
格子Boltzmann方法（Lattice Boltzmann Method，LBM）被学者提出来已有近30年的时间［1］。经过几十年的发展，LBM理论及工程应用研究等方面都有新的进展，并取得了卓越的成果。近年来，LBM理论及应用研究成为研究热点之一，受到国内外众多学者的重视。其基本思想是将系统中流体离散成
收稿日期：2022-05-24
作者简介：程海龙（1992—），男，硕士，研究方向：建筑工程管理。
大量的流体粒子，在简单划定的网格上，有规律地进行反复碰撞和迁移，给出一个简化的动力模型。然后通过统计来求取平均值，得到宏观物理量（如密度、速度）的流体动力学方程［2］。在实际应用中，有很多以压力进行驱动的不可压缩流体有周期性流动的现象［3］，特别是在工业生产及医学研究中，泊肃叶流的应用效果显著。最初泊肃叶流动是为了对血管中的血液流动进行描述的，在生物力学中有着重要的地位，且许多微尺度换热器中的流体流动也是用泊肃叶定律进行描述的。
1
物理模型及相关假设
如图1所示，两个平行的平板间充满流体，以沿x 方向的压强梯度为动力，推动整个流体流动，这就叫泊肃叶流［4］
。该流动是流体力学中十分经
典的流动问题，且可用Navier-Stokes 直接算出其解析解［5］，因而对基于格子Boltzmann 方法所开发的程序，经常使用该算例进行验证［6］。假定所述研究对象满足以下假设：平行板之间的流动为不可压流动；流动为二维流动，且只考虑XY 方向；满足刚性固体壁的边界条件。
2
控制方程
2.1
D2Q9模型
二维定常流动一般采用二维九速的D2Q9模型进行描述，该模型中的每个节点上都有9个不同的速度方向。其中，每个节点上允许有一个静止粒子存在，加上与其相邻的8个节点，记为D2Q9模型。模型采用的是D2Q9模型平衡态分布函数，见式（1）。
f eq
α=ρωα[1+e α·u c 2s +(e α·u )22c 4s
-u 22c 2s ]（1）
式中：
f eq α为平衡态分布函数；e α为速度向量；c s 为格子声速；
ωα为权系数；ρ为宏观密度；u 为宏观速度。在D2Q9模型中，c s =c 3，ω0=49，ω1,2,3,4=19，
ω5,6,7,8=1
36
初二英语作文。c =δx /δt ，δx 为网格步长，δt 为时间步
长，假定x 和y 方向的网格步长均相等，且等于时间
步长，单位均为1，即δx =δy =δt ，可得c =1。
2.2
股票型基金宏观方程及格子Boltzmann-BGK 方程
为恢复宏观方程，平衡态分布函数需满足的矩方程见式（2）至式（4）。
∑αf
eq α
=ρ（2）∑α
f
eq α
e α=ρu
（3）∑α
f
eq α
e αi e αj =ρu i u j +p δij
（4）
式中：δij =e i ·e j ，且p =ρc 2s ；
qq授权u 为宏观速度；δij 为一个矢量；p 为宏观压力。因此，模型的宏观密度、速度、压力定义见式（5）至式（8）。
ρ=∑α
f
eq
α
（5）u =
1
ρ∑αf eq αe α（6）p =ρc 2s
（7）f α(x +e αδt ,t +δt )-f α(x ,t )=-
1
τ
[f α（x ,t ）-f eq
α(x ,t )]+δtF α(x ,t )
（8）
式（8）是考虑外力项影响时的格子Boltzmann-
BGK 方程。式中：
F α(x ,t )是外力分布函数；τ=τ0/δt 为无量纲松弛时间，它直接影响流体的黏性和
宏观属性。对上述方程进行求解就是格子
Boltzmann 方法的核心。
3
人生观价值观
计算流程
本研究以visual studio 平台为整个程序的运行环境，编程语言选用C++。用LBM 模型对泊肃叶流进行数值模拟计算。该方法在商业软件数值模拟中具有优势，可有针对性地进行计算，通过简单程序可循环进行碰撞和迁移两个过程，直到符合设置的精度，因而程序的编写和实施容易实现，且对计算机系
统的要求不高。本研究的程序计算流程如下所示。
①初始化整个流场，确定各个节点的宏观量，并同时初始化分布函数，见式（9）。
f (x ,0)(i =1,2,...,b )
（9）
②在时刻t 执行碰撞，见式（10）。f 'i (x ,t )=f i (x ,t )+Ωi (x ,t )(i =1,2,...,b )（10）
③执行迁移，见式（11）。
f (x +c i δt ,t +δt )=f 't (x ,t )(i =1,2,...,b )（11）
④计算宏观量，见式（12）。
ρ(x ,t +δt )=∑i
f i (x ,t +δt )
ρu (x ,t +δt )=∑i
c i f i (x ,t +δt )
（12）
⑤重复步骤②③④，直到满足终止条件。
图1二维平板间的流动模型图
进口
出口
Y
X
4
数值方法
本研究编写了基于LBM 模型的程序，通过运行程序，将模拟值、相应的解析值与已有的相关文献给出
的结果进行分析对比，以验证所选模型和开发出的程序的正确性。用格子Boltzmann 方法对模拟二维速度驱动下的泊肃叶流动和用商业软件有限容积法模拟得出的结果进行对比，二者的结果非常吻合。本研究采用D2Q9模型，分别对速度驱动和压力驱动进行模拟，并对模拟结果进行对比分析［7］。与现有的研究相比［8］，模拟中均采用格子单位。在速度驱动下，一是进口处给定恒定速度u in =0.01，而压力驱动则没有初速度；二是给定一
个沿x 方向的力，大小为1.0×10-6
，y 轴方向的力为0，两种情况均采用均匀网格，网格数为N x ×N y =500×50，p r 数为0.71，粒子速度分布函数对应的松弛时
间、进口边界上的宏观参数速度已知，并假定进口处的分布函数满足相应的平衡态分布，右端出口为
充分发展边界，粒子速度分布函数和温度等于临近的流体节点对应的值，上下壁面为无滑移壁面，并采用非平衡外推格式对壁面进行处理。
5
结果分析
5.1
速度驱动
按照上述设定进行仿真模拟试验，试验结果见图2、图3。
图2是X =40截面处的速度分布图，经过比较可得，在X =40截面处速度取得最大值，并且可以明显看出，在流动达到稳定后，沿着平板方向的截面速度呈抛物线分布。由流体力学知识可知，在此条件下的最大速度是平均速度的1.5倍［9］。由图2可知，模拟结果为0.014698，与解析解最大误差为
2.01%，可以看出该方法比文献解更接近解析解，说
明格子Boltzmann 方法在处理二维平板间的不可压缩流动问题具有一定的精确度和可行性［10］。中心线上的速度大小在x 方向是逐渐减小的，这是因为压强是线性变化的，压强梯度是流动的动力，当达到稳定后，速度大小保持不变。在以X /L x 为横坐标、（p -p out ）/（p in -p out ）为纵坐标的曲线上，压力变化
如图3所示，入口和出口处稍有偏差，这是因为采用的边界处理格式不同，即假定入口处分布函数完全满足给定条件下的平衡态，而出口处采用的是充分发展格式，和实际还是存在偏差的。即在x 方向上存在压力梯度。压力梯度也是导致平板间流体流动的原因，压力梯度是平板间流体流动的驱动力。
男女身高差5.2
压力驱动
压力驱动和速度驱动的最大区别在于，压力驱动没有给定初速度，只给了一个沿x 方向的力，当流动达到稳定后，其流动情况类似于速度驱动。在y 方向上几乎没有速度，而x 方向上的速度呈对称分布，并且从边缘处到中间是逐渐增大的，这是因为驱动力较小时，要考虑流体的黏度和边界层因素的影响。在边界层的影响下，同一垂直界面上的流体的不同位置的速度不同，中间位置离边界层最远，受边界层的影响最小，因此速度较大（见图4）。
图5给出了流动稳定后所获得的速度场矢量图，从图5可以看出，当流动达到稳定后，速度矢量均指向x 方向，速度大小随着y 值的变化而变化，且呈抛物线分布，但这并不是在入口处就形成的。这是因为当流体以某一较小速度均匀进入平板间时，受黏性的影响，紧挨着壁面的流体会完全粘在壁面上，流体的速度逐渐降为0。根据流量连续原理，流过各个截面上的流量是相等的，因而越靠近平板的中心轴线，流速就越大，当流动达到稳定后，各截面上
图2X =40截面的速度分布
文献解解析解模拟解reply
0.0160.0140.0120.0100.0080.0060.0040.0020.000
Y /L y
V
0.00.2
0.4
0.6
0.8
1.0
图3
截面Y /L y =0.5的压力曲线
X /L y
0.0
0.2
0.4
0.6
0.8
1.0
p
0.39250.39200.39150.39100.39050.3900
的流速分布规律就不再变化，呈抛物线形状。6结语
本研究采用格子Boltzmann 方法，选取基本的D2Q9模型、标准的碰撞迁移，对速度驱动和压力驱动的泊肃叶流进行数值模拟，根据泊肃叶流的特征，边界处理采用非平衡态外推格式和充分发
展格式两种方式。试验结果表明，基于LBM 模型得出的模拟解与泊肃叶流速度的解析解吻合紧密。同时，不同纵截面上的速度分布以及整体变化趋势同现有的研究成果也保持一致，这说明格
子Boltzmann 方法完全适合处理压力驱动与速度驱动类的流动问题，并且具有很高的数值精度和数值稳定性，与其他方法相比，格子Boltzmann 方
法具有很多优势。通过分析速度云图和矢量图可知，压力梯度是流体流动的驱动力，并且通过对比速度驱动和压力驱动，对泊肃叶流动机理进行分析，为进一步研究流体的流动与传热奠定理论基础。
参考文献：
［1］CHEN S ，DOOLEN G D.Lattice Boltzmann method for fluid flows ［J ］.Annual Review of Fluid Mechanics ，1998：329-364
满月请帖
［2］陶文铨.数值传热学［M ］.2版.西安：西安交通大学出版社，2001.
［3］Roach P J.Computationnal Fluid Dynamics ［M ］.Al ‐buquerque ：Hemosa Publishers ，1976.
［4］STUART J T.On the non-linear mechanics of wave disturbances in stable and unstable parallel flows Part1.The basic behaviour in plane Poiuille flow ［J ］.Journal of Fluid Mechanics ，1960（3）：353-370.
［5］CHEN S ，DAWSON S P ，DOOLEN G D.Lattice Boltzmann model for simulation of magnetohydrodynam ‐ics ［J ］.Physical Review Letters ，1991（27）：637-667.
［6］陶文铨.计算传热学的近代发展［M ］.北京：科学出版社，2000.
［7］郭照立，郑楚光.格子Boltzmann 方法的原理及应用［M ］.北京：科学出版社，2008.
［8］吴念.基于格子Boltzmann 方法的散热数值仿真研究［D ］.哈尔滨：哈尔滨工程大学，2010
［9］何雅玲，王勇，李庆.格子Boltzmann 方法的理论及应用［M ］.北京：科学出版社，2009.
［10］TEIXEIRA C M.Incorporating turbulence models into the lattice-Boltzmann method ［J ］.International Journal of Modern Physics C ，1998（8）：1159-1175.
图4
截面X /L x =0.5的速度分布
L y
V
0.015
0.010
0.005
101520
25
30354045
05图5速度矢量图
X /L y
10.80.60.40.20
怎样写请柬Y /L y
0.20.4
0.60.80 1.0
