基于格子Boltzmann方法的非常规天然气微尺度流动基础模型

赵玉龙，刘香禺，张烈辉，单保超

（1. 西南石油大学油气藏地质及开发工程国家重点实验室，成都 610500；2. 华中科技大学煤燃烧国家重点实验室，武汉 430074）

0 引言

气体的微尺度流动模拟研究极大地促进了科技进步，影响着微电子机械系统、生物制药、仪表仪器、半导体材料以及航天航空等领域的发展[1-2]。近年来，页岩气、煤层气和致密砂岩气等非常规天然气资源的高效开发受到广泛关注。由于此类储集层岩石致密，气体流动空间多为微纳米级孔道且流动机理复杂[3-7]，常规物理实验和数值模拟方法无法揭示气体微观流动细节。为此，研究者开始利用微尺度模拟技术对非常规天然气进行流动模拟研究[7-11]。

在微纳米级别通道中，气体分子数目过少不足以支撑连续介质假设，另一方面，气体流动存在压缩效应和滑脱现象，且壁面附近的边界努森层对流动有重要影响，使得微尺度气体流动表现出不同于常规宏观尺度流动的特征。努森数（Kn）是气体流动特征控制参数，定义为分子平均自由程和流动特征长度的比值。依据Kn取值范围，可将气体流态划分为连续流（Kn≤0.001）、滑移流（0.001＜Kn≤0.1）、过渡流（0.1＜Kn≤10）和分子自由流（Kn＞10）[12]。在连续流区，连续介质假设成立，气体流动满足N-S（Navier-Stokes）方程；在滑移流区，连续介质假设基本成立，采用考虑边界滑移修正的 N-S方程进行模拟[13-15]；在过渡流区，连续介质假设失效，可采用Burnett方程进行描述；在分子自由流区，一般采用分子动力学方法进行模拟。在微尺度复杂流动体系中，流动通道尺寸变化范围大，气体可能同时存在上述几种流动状态。采用直接模拟蒙特卡洛方法（DSMC）和分子动力学方法可以对所有流态气体进行模拟，但是计算成本过于高昂；考虑不同模拟方法的适用条件，对混合流态流动采用不同方法进行耦合计算的处理方式难度大、计算复杂，且不能保证模拟结果的正确性。此时，可以选择用于模拟流体流动和输运现象的可行方法——格子 Boltzmann方法（LBM）[7-8,16]。不同于传统计算流体力学的方法，LBM基于微观粒子和介观动力学方程，具有清晰的物理背景，易于编程实现，并且具有天然的并行性，可以处理复杂边界，在研究非常规天然气微尺度流动时极具优势，进而得到广泛应用[11,16-27]。

采用LBM模拟微尺度气体流动时，无因次松弛时间和边界条件的确定至关重要[28]。目前通过分子动力学理论推导建立无因次松弛时间和Kn的关系式，从而利用Kn确定模拟条件下的无因次松弛时间[28-29]。由于前人在推导公式过程中采用不同的计算精度和物理参数计算表达式，关于Kn和无因次松弛时间的关系式非常多[14-15,28,30-35]，还考虑了不同因素的影响对无因次松弛时间进行修正[36-42]。前人均基于分子动力学理论推导无因次松弛时间表达式，采用其他思路求取无因次松弛时间的报道尚未出现。本文在前人研究的基础上，考虑相似准则和气体真实物性参数，推导建立了适用于高温高压条件的非常规天然气微尺度流动模拟新模型。

1 基于LBM的微尺度流动模拟新模型

1.1 LBM基本模型理论

1992年，Qian等[43]提出的 LBM 基本模型——DnQm模型得到广泛应用。模拟非常规天然气流动时，常采用BGK近似碰撞项[44]的单松弛格子Boltzmann模型（LBGK）。LBGK模型是连续Boltzmann-BGK方程的特殊离散形式，而连续Boltzmann-BGK方程是连续Boltzmann方程的简化。对于一个完整的LBGK模型，应该包括演化方程、平衡态分布函数和格子模型，其中演化方程表示如下：

在模拟二维空间流动时，常采用九速离散的D2Q9格子模型，此时，平衡态分布函数为：

格子空间动量表示为：

在D2Q9格子模型中，权重系数iω和离散格子速度c分别表示为：

此外，格子空间中运动黏度与无因次松弛时间满足关系式[45]：

上述黏度处理方式使得LBGK模型在计算时具有二阶精度。格子空间中，单组分单相气体压力表示为[28]：

外加作用力的离散表达式[28]：

1.2 无因次松弛时间的选取

表1 文献中无因次松弛时间和努森数关系

在微尺度流动模拟中，还需要考虑边界努森层带来的影响。固体壁面的存在使得靠近壁面位置的分子平均自由程被截断，导致一定范围内气体分子（即努森层）的运动规律受到影响[8,11]。当流动通道特征长度较大时，努森层在整个流动通道中所占比例小，其对流动的影响可以忽略不计；但随着特征长度的减小，努森层在流动通道内所占比例逐渐增大，其对流动产生的影响也越来越大[8]，产生微尺度流动特征。因此，研究者对无因次松弛时间进行了相应的修正。总体来说，主要通过修正分子平均自由程以表征边界努森层的影响。一些学者将努森层中气体分子平均自由程减小等价转换为整个流场中气体分子平均自由程的减小。Stops等[37]提出了二维平板流动下的分子平均自由程修正式λe=λψ（Kn） ，但修正系数ψ（Kn） 过于复杂，实际计算不方便；Guo等[32]在此基础上提出新的修正系数表达式；Li等[38]根据气体动力黏度与分子平均自由程的正比例关系，采用 Michalis等[39]提出的Bosanquet-type型修正系数。Zhang等[40]、Tang等[41]、Guo等[42]考虑流场中不同位置与壁面的距离对该位置处分子平均自由程进行修正，推导出相应的无因次松弛时间。

上述无因次松弛时间表达式均通过分子动力学理论推导而来，在微尺度流动模拟研究中应用广泛。如图 1所示，作出部分表达式中无因次松弛时间随努森数的变化曲线，可以看到在相同网格数及努森数条件下，无因次松弛时间值存在差异，并且在过渡流区尤其明显。无因次松弛时间对微尺度气体流动模拟研究影响重大，微小的取值差别会导致模拟结果出现差异，可能造成模拟结果偏离正确值。

图1 不同努森数下无因次松弛时间取值（N=100）

前人在推导无因次松弛时间时，分子平均自由程与松弛时间、粒子平均运动速度满足关系式[36]：

动力黏度和压力、松弛时间满足关系[32]：

考虑高温高压非理想气体稠密性修正时[36]，动力黏度可表示为：

不考虑壁面产生的影响，绘制出温度为350 K时（11）式对应的真实物理空间中甲烷气体分子动力黏度与压力（1～100 MPa）关系图（见图2），并与美国国家标准与技术研究院（NIST）化学数据库中的数据对比。从图中可以看出，压力高于20 MPa时，各表达式计算的甲烷气体分子动力黏度均与 NIST数据库中数据存在较大偏差。此时，上述无因次松弛时间表达式无法反映气体真实情况。

图2 不同表达式计算的的甲烷动力黏度随压力的变化

在微尺度流动模拟中，应使模拟的格子空间量和真实物理空间量满足相似准则[46-47]，如马赫数、努森数及雷诺数相等，才能真正反映实际情况，揭示客观规律。在宏观液体流动模拟中，考虑数值模拟条件与真实流动条件下的雷诺数相等，实际上是满足格子空间和真实物理空间中的惯性力与摩擦力之比相等。值得注意的是，对于宏观大尺度液体流动，可忽略液体的压缩性，但针对微观气体流动应该考虑气体的压缩性。此时，需要严格满足模拟条件下的马赫数和实际马赫数相等，也就是满足二者的惯性力和弹性力之比相等。

由格子空间和真实物理空间下的雷诺数相等可得：

格子空间和真实物理空间中的努森数相等，则：

满足格子空间和真实物理空间中的马赫数相等，则：

联立（6）式及（12）—（14）式可得：

真实物理空间中气体分子平均自由程计算式为：

前人的研究结果表明，可以通过修正无因次松弛时间表达式以表征微尺度条件下固体壁面附近努森层对流动的影响[32,38]：

图3 不同努森数下修正系数的取值及差异

1.3 边界条件及其关键参数的确定

一般微尺度模拟中，考虑气体在微通道内周期流动时，通常在流动方向进出口处采用周期边界条件[47]；考虑气体在微通道中受进出口压差驱动流动时，在进出口处采用Zou等[48]提出的压力边界条件或Guo等[49]提出的非平衡态外推边界条件。相比之下孔道上下壁面处边界处理更复杂，通过在壁面处采用合理的边界条件可以有效捕捉边界滑移等现象。

目前，处理壁面处微尺度气体流动时通常采用两种边界条件[7,35]：组合反弹/镜面反射边界条件（CBBSR）和离散 Maxwell边界条件（DM）。这两种边界条件在本质上是等价的，本文采用CBBSR边界条件[35]。以二维微通道下边界为例，网格划分如图 4所示，本文将网格节点设置在壁面边界位置处有利于直接捕捉边界滑移速度。经过碰撞步和迁移步后，未知分布函数为f2、f5和f6，采用CBBSR边界条件可以得到：

图4 网格划分示意图（数字表示离散速度方向）

在微尺度流动模拟中，组合系数r的取值直接影响模拟得到的边界滑移速度，若将r设为定值[50]则无法反映出模拟条件的差异。为消除数值离散效应并使模拟满足二阶滑移边界条件，根据Guo等[35]的相关理论，结合本文提出的无因次松弛时间推导得出组合系数表达式：Guo等[51]在广义二阶滑移边界条件中给出：

1.1 资料来源 选择2013年2月-2017年7月在本院住院和门诊就诊的妊娠孕妇3 920例，年龄21～40岁，平均年龄(27.23±2.31)岁；经产妇2 137例，初产妇1 783例；孕周16～41周，平均孕周(28.34±6.24)周。纳入标准：自愿参与本次研究，并签署知情同意书者；外院检查疑似或确诊有畸形者；早期有感冒病史；经引产或分娩证实为胎儿畸形。排除标准：合并妊娠糖尿病、妊娠期高血压疾病等疾病者；存在胎儿宫内生长受限、脐绕颈、羊水异常等。本研究经我院伦理委员会审核并批准。

2 模型验证

2.1 无限长微通道中体积力驱动模拟验证

首先，利用等温条件下甲烷气体在无限长微通道中的体积力驱动流动进行模型验证。通过在进出口设置周期流动边界条件[47]，可以将气体流动等效为无限长微通道中的流动。通常将压力梯度设置为很小的常量[38]，并且将整个流场中外力项统一为F=（Fx,Fy），本模拟取Fx为10−6，Fy为零。此时，沿流动方向气体压力变化微小，可以忽略气体宏观物性参数差异。将整个流场中的无因次松弛时间设为恒定，同时可以忽略压缩效应的影响，只需要考虑气体的滑脱效应和边界努森层的影响。网格划分为Nx×Ny=30×300，由于采用周期流动边界条件，流动方向（即x方向）网格数少一些。（18）式和（20）式中涉及到的所有甲烷气体物性参数均来自NIST标准参考数据库，通过改变模拟的微通道上下壁面间距Hr可以得到不同Kn及对应的无因次松弛时间和组合系数。由于需要模拟较高Kn下的流动，为保证通道特征长度为纳米级，模拟压力不能太高，因此取Tr为400 K，pr为0.1 MPa。

由图5可以看出，当ψ（Kn） =ψ（Kn）_1时，Kn取值不同本文模型计算精度不同。当Kn值为0.112 8时，本文模型计算结果与 Ohwada等直接求解线性化Boltzmann方程得到的结果具有较好的一致性，但是在中心线位置处本文模型计算结果大于Zhang等和Guo等的结果；随着Kn值的增大，本文模型计算结果逐渐偏离Ohwada等的结果，具体表现为边界处滑移速度偏大，而通道中间速度偏小，与 Zhang等的计算结果呈现出很好的一致性；当Kn＞2.5，模拟不收敛。当ψ（Kn） =ψ（Kn）_2时，本文模型 计 算 结 果与 Ohwada等的模拟结果始终保持很好的一致性，说明此时模型具有更高的计算精度，可以在较大的Kn变化范围内有效捕捉微尺度气体流动边界滑移现象。此外，本文模拟结果与Li等（Li等使用了相同的修正系数ψ（Kn）_2）的计算结果吻合度较高，说明ψ（Kn）对边界滑移速度的准确捕捉具有重要影响。本文模型计算结果与 Guo等计算结果始终存在偏差，并且在Kn值较大时偏差非常明显。可以发现，随着Kn值的增加，边界滑移速度显著增大，滑脱效应对流动的影响越来越明显。

图5 甲烷气体在微通道中周期流动时出口处无因次速度剖面

图6 甲烷在微通道中周期流动时无因次流量随努森数的变化

在上述两种修正系数下计算出不同Kn值条件下出口处的归一化边界滑移速度Uslip/U0。将计算结果同Arkilic等[55]考虑一阶滑移边界推导得到的解析解以及Nie等[30]的结果进行对比（见图 7）。可以看出，在滑移流区内采用不同修正系数得到的计算结果基本吻合，但随着Kn值的增加二者差距逐渐显现；在过渡流区内，二者差异进一步增加。同时，在整个Kn取值范围内，当ψ（Kn） =ψ（Kn）_1时，本文模型计算结果同Arkilic等的解析解具有更好的一致性。实际上，在Arkilic等的推导过程中仅考虑一阶滑移边界条件，因此采用ψ（Kn）_1时本文模型只具有一阶计算精度。在Nie等的研究中，采用的是简单的反弹边界条件，可以看到在滑移流区，Nie等的计算结果小于其余3种计算结果，在过渡流区则相反。当ψ（Kn） =ψ（Kn）_2时，过渡流区本文模型模拟计算结果小于其余 3种结果，说明采用ψ（Kn）_1导致计算得到的边界滑移速度在过渡流区偏大。

图7 甲烷在微尺度二维平板中周期流动时归一化边界滑移速度随努森数的变化

2.2 长直微通道中进出口压差驱动模拟验证

采用本文模型模拟甲烷在长直微通道中受进出口压差驱动的流动。对于微尺度气体流动，在研究范围内温度压力变化微小，气体流动速度缓慢，属于典型的低马赫数流动。虽然理论上可以将微尺度气体流动视为不可压缩流动，但是在整个宏观尺度上来看，压缩性是客观存在的，忽略压缩效应的模拟必然存在误差。在研究此类问题时，需要同时考虑压缩效应、气体滑脱效应和边界努森层的影响。由于气体的可压缩性，微通道内沿流动方向压力表现为非线性分布。将本文模拟结果同直接模拟蒙特卡洛方法（DSMC）计算结果、信息保存-直接模拟蒙特卡洛方法（IP-DSMC）计算结果、Li等[38]、Guo等[53]、Shen等[56]的LBM模拟结果以及Arkilic等[55]的解析解进行对比，以验证模型的准确性。在上述文献中，进出口压力比为1.4∶1.0或2.0∶1.0，考虑模拟的压力梯度不能太大，设定出口压力为0.01 MPa，模拟温度恒定为400 K，通过改变微通道上下壁面间距以实现不同Kn值条件下的流动模拟。气体沿流动通道方向压力发生变化，此时在整个流场内气体的宏观物性参数也会发生相应变化。根据（18）式和（20）式，可以看出不同位置处的无因次松弛时间、努森数以及组合系数不同，在模拟时需要考虑各参数值随位置的变化。流场任意位置处努森数表示为Kn（x,y） =Knoutρout/ρ（x,y） 。值得注意的是，在计算无因次松弛时间和组合系数时涉及到νr/csr，需要进行相应的预处理。通过 NIST数据库查询得到νrcsr和ρr在进出口压力范围内的一系列离散值，然后拟合出二者关系式νr/csr=3×10-8/ρr。实际模拟时，设置格子空间密度和真实物理空间密度在数值上相等，此时可结合上述关系式计算出任意位置处的无因次松弛时间和组合系数。

取ψ（Kn） =ψ（Kn）_2，采用本文模型计算出Knout分别为0.019 4、0.194 0和0.388 0时的沿程无因次压力偏移量 （p-pline）/pout。网格划分为NxNy=2 000×20，与参考文献[55]一致。

从图 8可以看出，采用本文模型计算得到的结果明显优于Arkilic等、Shen等和Guo等的结果，与Li等的计算结果相近；但相比之下本文模型计算结果与DSMC、IP-DSMC法的模拟结果更具有一致性，说明本文模型能更为准确地捕捉微尺度气体的压缩效应。

图8 本文模型计算的长直微通道中气体沿程无因次压力偏移量与文献数据对比

3 结论

与基于连续介质假设的传统计算流体力学方法相比，格子 Boltzmann方法由于其介观特性及天然的并行性等特点，在微尺度气体流动模拟研究方面更具优越性，但当前广泛采用的格子 Boltzmann模型恢复出的气体黏度在高温高压条件下与真实值存在较大差异，可能导致模拟结果偏离正确值。

本文模型可有效表征微尺度气体滑脱效应和压缩效应，准确捕捉边界努森层的影响。对甲烷在微通道中的周期流动进行模拟，结果表明采用 Michalis等提出的 Bosanquet-type型黏度修正系数使本文模型具有更高的计算精度和更广的适用范围。此时，采用本文模型计算得到的无因次边界滑脱速度和无因次流量与文献中数据具有很好的一致性。模拟甲烷在长直微通道中受进出口压差驱替的流动，结果表明本文模型可有效表征微尺度气体压缩效应的影响。与文献中常用微纳尺度模型相比，本文模型可更为精确地捕捉边界滑移速度和气体压缩性的影响。

本文模型提出的无因次松弛时间表达式中除努森数外还包含气体在模拟温度和压力条件下的真实运动黏度、声速及气体分子平均自由程等参数，能更全面地反应非常规天然气在储集层条件下微尺度流动空间中的真实状态。此外，结合吸附解吸、表面扩散等边界条件可将本文模型进行扩展应用。

符号注释：

A1，A2——二阶滑移边界条件中与气固相互作用有关的系数；a——常数；C——经验参数；c——格子速度；ci——格子空间i方向离散速度；cs——格子空间声速；csr——真实物理空间声速，m/s；——格子空间粒子平均速度；D——模拟空间维数；dr——真实物理空间气体分子直径，甲烷气体取dr= 0.414×10−9m；F——格子空间外力；Fx——格子空间x方向外力分量；Fy——格子空间y方向外力分量；Fi——格子空间i方向离散作用力；fi(X,t)——粒子分布函数；fi,eq(X,t)——平衡态粒子分布函数；——壁面节点处经过迁移步之后的粒子分布函数；H——格子空间流动特征长度；Hr——真实物理空间流动特征长度，m；i——离散速度方向，i=0，1，…，m-1；I——二阶单位张量；j——y方向的网格序号；Kn——努森数，无因次；Knout——微通道出口处努森数，无因次；m——离散速度数；mr——真实物理空间气体分子质量，甲烷气体取mr=2.664×10−26kg；N——特征长度所占网格数；Nx——x方向网格数；Ny——y方向网格数；n——空间维数；p——格子空间压力；pline——格子空间压力线性分布值；pin——格子空间微通道入口处压力；pout——格子空间出口处压力；pr——真实空间压力，MPa；Q——格子空间无因次流量；R——格子空间气体常数；r——组合系数，控制反弹与镜面反射间的比例，取 0＜r＜1；T——格子空间温度；Tr——真实物理空间温度，K；t——格子空间时间；U——格子空间微通道出口处沿流体流动方向的速度；Uave——格子空间微通道出口处沿流体流动方向的平均速度；Uslip——格子空间微通道出口处边界滑移速度；U0——格子空间微通道出口处中轴线上流动方向速度；u——格子空间宏观速度矢量；u——格子空间速度标量；ur——真实物理空间运动速度，m/s；X——格子空间位置；x——距离入口端的距离；y——距离微通道下壁面的距离；γ——格子空间气体比热容；δt——格子空间时间步；δx——格子空间网格间距，通常取 δx=δt；λ——格子空间气体分子平均自由程；λe——格子空间有效气体分子平均自由程；λr——真实物理空间气体分子平均自由程，m；ρ——格子空间宏观密度；ρout——格子空间微通道出口处宏观密度；ρr——真实物理空间气体密度，kg/m3；σ——流向动量调节系数；σv——切向动量调节系数，通常取σv=1.0；τ——无因次松弛时间；τs——格子空间松弛时间；μ——格子空间动力黏度；μr——真实物理空间动力黏度，Pa·s；ν——格子空间运动黏度；νr——真实物理空间运动黏度，m2/s；ωi——i方向权重系数；χ——气体稠密性修正系数，无因次；ψ(Kn)——平均分子自由程修正系数；Δψ(Kn)——修正系数差。