程海龙

（河南省机场集团有限公司，河南 郑州 450000）

0 引言

格子Boltzmann方法（Lattice Boltzmann Method，LBM）被学者提出来已有近30 年的时间［1］。经过几十年的发展，LBM 理论及工程应用研究等方面都有新的进展，并取得了卓越的成果。近年来，LBM 理论及应用研究成为研究热点之一，受到国内外众多学者的重视。其基本思想是将系统中流体离散成大量的流体粒子，在简单划定的网格上，有规律地进行反复碰撞和迁移，给出一个简化的动力模型。然后通过统计来求取平均值，得到宏观物理量（如密度、速度）的流体动力学方程［2］。在实际应用中，有很多以压力进行驱动的不可压缩流体有周期性流动的现象［3］，特别是在工业生产及医学研究中，泊肃叶流的应用效果显著。最初泊肃叶流动是为了对血管中的血液流动进行描述的，在生物力学中有着重要的地位，且许多微尺度换热器中的流体流动也是用泊肃叶定律进行描述的。

1 物理模型及相关假设

如图1 所示，两个平行的平板间充满流体，以沿x方向的压强梯度为动力，推动整个流体流动，这就叫泊肃叶流［4］。该流动是流体力学中十分经典的流动问题，且可用Navier-Stokes 直接算出其解析解［5］，因而对基于格子Boltzmann方法所开发的程序，经常使用该算例进行验证［6］。假定所述研究对象满足以下假设：平行板之间的流动为不可压流动；流动为二维流动，且只考虑XY方向；满足刚性固体壁的边界条件。

图1 二维平板间的流动模型图

2 控制方程

2.1 D2Q9模型

二维定常流动一般采用二维九速的D2Q9模型进行描述，该模型中的每个节点上都有9 个不同的速度方向。其中，每个节点上允许有一个静止粒子存在，加上与其相邻的8 个节点，记为D2Q9 模型。模型采用的是D2Q9模型平衡态分布函数，见式（1）。

式中：f eqα为平衡态分布函数；eα为速度向量；cs为格子声速；ωα为权系数；ρ为宏观密度；u为宏观速度。在D2Q9模型中，cs=c=δx/δt，δx为网格步长，δt为时间步长，假定x和y方向的网格步长均相等，且等于时间步长，单位均为1，即δx=δy=δt，可得c= 1。

2.2 宏观方程及格子Boltzmann-BGK 方程

为恢复宏观方程，平衡态分布函数需满足的矩方程见式（2）至式（4）。

式中：δij=ei·ej，且p=ρc2s；u为宏观速度；δij为一个矢量；p为宏观压力。因此，模型的宏观密度、速度、压力定义见式（5）至式（8）。

式（8）是考虑外力项影响时的格子Boltzmann-BGK 方程。式中：Fα(x,t) 是外力分布函数；τ=τ0/δt为无量纲松弛时间，它直接影响流体的黏性和宏观属性。对上述方程进行求解就是格子Boltzmann方法的核心。

3 计算流程

本研究以visual studio 平台为整个程序的运行环境，编程语言选用C++。用LBM 模型对泊肃叶流进行数值模拟计算。该方法在商业软件数值模拟中具有优势，可有针对性地进行计算，通过简单程序可循环进行碰撞和迁移两个过程，直到符合设置的精度，因而程序的编写和实施容易实现，且对计算机系统的要求不高。本研究的程序计算流程如下所示。

①初始化整个流场，确定各个节点的宏观量，并同时初始化分布函数，见式（9）。

②在时刻t执行碰撞，见式（10）。

③执行迁移，见式（11）。

④计算宏观量，见式（12）。

⑤重复步骤②③④，直到满足终止条件。

4 数值方法

本研究编写了基于LBM 模型的程序，通过运行程序，将模拟值、相应的解析值与已有的相关文献给出的结果进行分析对比，以验证所选模型和开发出的程序的正确性。用格子Boltzmann 方法对模拟二维速度驱动下的泊肃叶流动和用商业软件有限容积法模拟得出的结果进行对比，二者的结果非常吻合。本研究采用D2Q9 模型，分别对速度驱动和压力驱动进行模拟，并对模拟结果进行对比分析［7］。与现有的研究相比［8］，模拟中均采用格子单位。在速度驱动下，一是进口处给定恒定速度uin= 0.01，而压力驱动则没有初速度；二是给定一个沿x方向的力，大小为1.0×10-6，y轴方向的力为0，两种情况均采用均匀网格，网格数为Nx×Ny=500×50，pr数为0.71，粒子速度分布函数对应的松弛时间、进口边界上的宏观参数速度已知，并假定进口处的分布函数满足相应的平衡态分布，右端出口为充分发展边界，粒子速度分布函数和温度等于临近的流体节点对应的值，上下壁面为无滑移壁面，并采用非平衡外推格式对壁面进行处理。

5 结果分析

5.1 速度驱动

按照上述设定进行仿真模拟试验，试验结果见图2、图3。

图2 X=40截面的速度分布

图2是X=40截面处的速度分布图，经过比较可得，在X=40 截面处速度取得最大值，并且可以明显看出，在流动达到稳定后，沿着平板方向的截面速度呈抛物线分布。由流体力学知识可知，在此条件下的最大速度是平均速度的1.5倍［9］。由图2可知，模拟结果为0.014 698，与解析解最大误差为2.01%，可以看出该方法比文献解更接近解析解，说明格子Boltzmann 方法在处理二维平板间的不可压缩流动问题具有一定的精确度和可行性［10］。中心线上的速度大小在x方向是逐渐减小的，这是因为压强是线性变化的，压强梯度是流动的动力，当达到稳定后，速度大小保持不变。在以X/Lx为横坐标、（p-pout）/（pin-pout）为纵坐标的曲线上，压力变化如图3所示，入口和出口处稍有偏差，这是因为采用的边界处理格式不同，即假定入口处分布函数完全满足给定条件下的平衡态，而出口处采用的是充分发展格式，和实际还是存在偏差的。即在x方向上存在压力梯度。压力梯度也是导致平板间流体流动的原因，压力梯度是平板间流体流动的驱动力。

图3 截面Y/Ly=0.5的压力曲线

5.2 压力驱动

压力驱动和速度驱动的最大区别在于，压力驱动没有给定初速度，只给了一个沿x方向的力，当流动达到稳定后，其流动情况类似于速度驱动。在y方向上几乎没有速度，而x方向上的速度呈对称分布，并且从边缘处到中间是逐渐增大的，这是因为驱动力较小时，要考虑流体的黏度和边界层因素的影响。在边界层的影响下，同一垂直界面上的流体的不同位置的速度不同，中间位置离边界层最远，受边界层的影响最小，因此速度较大（见图4）。图5 给出了流动稳定后所获得的速度场矢量图，从图5 可以看出，当流动达到稳定后，速度矢量均指向x方向，速度大小随着y值的变化而变化，且呈抛物线分布，但这并不是在入口处就形成的。这是因为当流体以某一较小速度均匀进入平板间时，受黏性的影响，紧挨着壁面的流体会完全粘在壁面上，流体的速度逐渐降为0。根据流量连续原理，流过各个截面上的流量是相等的，因而越靠近平板的中心轴线，流速就越大，当流动达到稳定后，各截面上的流速分布规律就不再变化，呈抛物线形状。

图4 截面X/Lx=0.5的速度分布

图5 速度矢量图

6 结语

本研究采用格子Boltzmann 方法，选取基本的D2Q9 模型、标准的碰撞迁移，对速度驱动和压力驱动的泊肃叶流进行数值模拟，根据泊肃叶流的特征，边界处理采用非平衡态外推格式和充分发展格式两种方式。试验结果表明，基于LBM 模型得出的模拟解与泊肃叶流速度的解析解吻合紧密。同时，不同纵截面上的速度分布以及整体变化趋势同现有的研究成果也保持一致，这说明格子Boltzmann 方法完全适合处理压力驱动与速度驱动类的流动问题，并且具有很高的数值精度和数值稳定性，与其他方法相比，格子Boltzmann 方法具有很多优势。通过分析速度云图和矢量图可知，压力梯度是流体流动的驱动力，并且通过对比速度驱动和压力驱动，对泊肃叶流动机理进行分析，为进一步研究流体的流动与传热奠定理论基础。