REV尺度多孔介质格子Boltzmann方法的数学模型及应用的研究进展.pdf

CHEMICAL INDUSTRY AND ENGINEERING PROGRESS 2016 年第 35 卷第 6 期 1698 化 工 进 展 REV 尺度多孔介质格子 Boltzmann 方法的数学模型及应用的 研究进展 张潇丹 1，2，雍玉梅2，李文军3，赵元生4，李媛媛2，杨巧文1，杨超2 （1中国矿业大学（北京）化学与环境工程学院，北京 100083；2中国科学院过程工程研究所绿色过程与工程重 点实验室，北京 100190；3华北科技学院环境工程学院，河北 廊坊 065201；4中国石油化工研究院渣油加氢实 验室，北京 102200） 摘要：综述了多孔介质表征体元尺度（REV）格子 Boltzmann 模型的研究进展，根据对多孔介质处理方式主要分 为部分反弹模型和阻力模型两类，分析归纳了各类模型的优缺点。由于阻力模型中渗流的广义格子 Boltzmann 方程（GLBE）的作用力是基于 GUO 等的作用力模型，可以准确得到宏观方程，不存在离散误差，且模型的平 衡分布函数和作用力项中都包含反应介质特性的孔隙率，因而应用最为广泛。本文还重点介绍了 REV 尺度多孔 介质 LBE 模型在流动、传热、传质、化学反应及相变等过程中的具体应用，认为 REV 尺度多孔介质内的三传一 反数学模型中需要加入孔隙尺度因素，在更大工程尺度上应该考虑过程参数的各向异性，展望了 REV 尺度多孔 介质 LBE 模型的发展和应用前景。 关键词：多孔介质；表征体元尺度；格子 Boltzmann 方法；流动；传热；传质 中图分类号：TQ021.9 文献标志码：A 文章编号：10006613（2016）06169815 DOI：10.16085/j.issn.1000-6613.2016.06.010 Models and application of lattice Boltzmann method at REV-scale in porous media ZHANG Xiaodan1 ，2， YONG Yumei2， LI Wenjun3， ZHAO Yuansheng4， LI Yuanyuan2， YANG Qiaowen1， YANG Chao2 （1School of Chemical & Environmental Engineering，China University of Mining & Technology (Beijing)，Beijing 100083， China； 2Key Laboratory of Green Process and Engineering， Institute of Process Engineering， Chinese Academy of Sciences，Beijing 100190，China；3School of Environmental Engineering，North China Institute of Science and Technology，Langfang 065201，Hebei，China；4Laboratory of Residue Hydrotreating，Research Institute of Petroleum Processing，PetroChina，Beijing 102200，China） Abstract： This paper discusses the lattice Boltzmann model at representative elementary volume (REV) scale for porous media. According to different treatments of porous media， the lattice Boltzmann model at REV-scale for porous media can be classified into two categories，the partially bouncing-back model and the resistance model. The advantages and disadvantages of various models are analyzed. The Generalized lattice Boltzmann equation (GLBE model) in the resistance model is most widely used. Firstly， the force item of the GLBM model is based on the method proposed by Guo et al， which can be 第一作者：张潇丹（1989） ，女，硕士研究生，主要从事化学工程 数值模拟。联系人：雍玉梅，研究员，主要从事化学工程数值模拟 E-mail ymyong 。 特约评述 收稿日期：2016-01-19；修改稿日期：2016-01-31。 基金项目：国家自然科学基金（21276256,21490584） 、国家重点基础研究发 展计划（2012CB224806）及中央高校基本科研业务费（3142013097）项目。 第 6 期 张潇丹等：REV 尺度多孔介质格子 Boltzmann 方法的数学模型及应用的研究进展 1699 accurately recovered to the macroscopic equation without discretization error. Secondly，the equilibrium distribution function and the force items involve porosity， which reflects the characteristics of porous media. This review introduces the application of the lattice Boltzmann model at REV-scale in porous media with flow， heat transfer， mass transfer， chemical reaction and phase transition. We should take the pore scale factors into consideration when studying momentum transport， heat transport， mass transport and reaction engineering in porous media at REV-scale. The anisotropic nature of process parameters should be considered at the larger engineering scale. It gives some predictions and perspectives of applications for the lattice Boltzmann model at REV-scale. Key words：porous media；representative elementary volume (REV)；lattice Boltzmann method (LBM)；flow；heat transfer；mass transfer 流体在多孔介质内的流动是许多领域都涉及的 一类复杂流动，在能源化工、生物医药、空气净化、 污水处理等领域都有广泛应用。多孔介质内的流动 是一种典型的多尺度问题，通常涉及 3 个尺度，即 孔隙尺度、表征体元尺度（REV）和宏观尺度。格 子 Boltzmann 方法（LBM）是介于宏观和微观两种 方法之间的介观数值模拟方法，它既克服了宏观方 法需要处理复杂边界的劣势，又克服了微观方法计 算尺度受限的局限性，在模拟多孔介质方面具有独 特优势1。 由于 LBM 算法简单， 计算网格容易生成， 可以比较方便处理流体与边界之间、不同流体组分 之间等复杂的相互作用， 且程序代码简单容易并行， 应用逐渐广泛。格子 Boltzmann 方法模拟多孔介质 渗流主要有孔隙尺度和 REV 尺度。孔隙尺度 LBM 的研究对象是单个或者若干个孔隙内的流体，把介 质的骨架作为流场的边界，用来探索渗流机理和基 本规律。REV 尺度 LBM 是以一个控制体为研究对 象， 求解 REV 体积平均的宏观流动量， 流体与介质 的相互作用通过适当的处理方法来描述。 REV 尺度 LBM 方法相对于孔隙尺度 LBM 方法，优点是前者 不需要介质的细节结构、仅仅依赖于介质的统计参 数，因而计算效率较高，比较适合大区域的工程渗 流计算。已有文献总结了 LBM 方法的基本要素2， 介绍了 LBM 的起源、连续 Boltzmann 方程、边界 条件和数值稳定性，讨论了 LBM 在流体流动、多 相流和多组分流、反应扩散中的应用。CHEN 等3 则对 LBM 在湍流、多孔介质流、燃烧、多孔介质 内燃烧的应用进行了综述， 李楠等4分析了 LBM 在 模拟这些复杂流体问题时存在的优缺点。SUKOP 等5通过比较两种 REV 尺度 LBE 模型， 提出 Dardis 和McCloskey的LBE模型不能用于各向异性的多孔 介质的观点，但没有详细论述两模型的具体差异。 目前还没有 REV 尺度多孔介质 LBM 方法的综述， 因而本文总结了 REV 尺度多孔介质 LBE 各类数学 模型，分析了各模型的优缺点，综述了该方法在 REV 尺度多孔介质内流动、传热、传质、化学反应 和相变等领域的应用，分析了未来应用前景，并对 其发展趋势提出建议。 1 REV 尺度 LBE 模型 REV 尺度 LBE 模型根据对多孔介质的处理方 式不同主要分为两类：部分反弹模型6-11，此种 模型对经典的固壁回弹模型进行修改来拟合达西公 式；阻力模型12-21，此模型通过施加外力公式， 得到 Stokes 或 Brinkman 方程。 1.1 部分反弹模型 部分反弹模型由 3 个步构成，如式(1)。前两个 步是碰撞步和流动步， 与标准 LBE 模型相同。第三 个步是部分反弹步，该过程中流体分布函数根据该 点的固体密度重新定义。部分反弹模型中引入了多 孔介质连续表示方法，并将反弹边界的作用耦合到 演化方程中7。 (eq) 1 (,)( , )( , )( , )( , ) ii ttiiii fctftftftt += + xxxxx (1) 式中，i 是速度离散方向；t 是离散时间；f 为 分布函数； 为松弛时间；i (x,t)是流体粒子与固 体骨架之间的碰撞，表示多孔介质对流体影响。 最早的部分反弹模型是 GAO 和 SHARMA6提 出的多孔介质格子气自动模型，可延伸到 LBM 中 （GS-LBE 模型)。该模型能够模拟非均质多孔介质 内的流动。该模型中 i (x,t)如式(2)。 s ( , )=( , )( , ) ii i tnftft xxx (2) 式中，ns是固体格点的平均散射密度，取值范 围为 01，代表多孔介质的固体结构；i 是 i 的相 化 工 进 展 2016 年第 35 卷 1700 反方向。 但是 GS-LBE 模型存在格子气的统计噪声。 1998 年， DARDIS 和 MCCLOSKEY7则提出了第一 个部分反弹 LBE 模型（DM-LBE 模型)，是对 Brinkman 方程的间接求解。与 GS-LBE 模型相比， DM-LBE模型用ns代替孔隙尺度多孔介质的结构参 数布尔变量，消除了格子气的不稳定性。渗透率依 赖于固体格点的平均散射密度 ns，并且对应的有效 黏性系数与流体的黏性系数相同。 该模型中为 i (x,t) 如式(3)。 s ( , )(, )( , ) iiti i tnftft=+ xxcx (3) THORNE 和 SUKOP8根据部分反弹形式不同 提出了第二类部分反弹格式（TS-LBE 模型） ，该模 型流出分布函数反映临近点流入分布函数碰撞过 程，该类反弹格式需计算每个格点的碰撞步来估算 流出分布函数密度。该模型中 i (x,t)如式(4)。 c*c* s ( , )(,)( ,) ii ti i tnf tft =+ xxcx (4) 式中 t*是碰撞前的时间，满足 tt*t+1，fic 为碰撞后的分布函数。 DM-LBE模型和TS-LBE模型的部分反弹格式 需要根据临近格点计算得到，计算复杂，不利于并 行计算。WALSH、BURWINKLE 和 SAAR9引入 了新的部分反弹方法（WBS-LBE 模型） ，该方法 在碰撞步前改变流入流体方向， 同时计算碰撞步和 反弹步，且计算中不包含临近格点的流体分布函 数，使计算更容易并行，从而节省计算时间。与 DM-LBE模型和TS-LBE模型的部分反弹方法相比 较， 该方法在非均质多孔介质流动过程中满足质量 守恒，且提高了 LB 方法数值模拟浮力驱动及应力 和浓度梯度引发的扩散过程的精度。该模型中 i (x,t)如式(5)。 cc* s ( , )( , )( ,) ii i tnftft = xxx (5) WBS-LBE 模型虽然在非均质多孔介质内满足 质量守恒，更加适用于浮力驱动的流动和扩散过 程， 但是该模型在不同渗透特性的多孔介质界面处 速度不连续，ns接近 1 时会导致数值不稳定。 YEHYA、NAJI 和 SUKOP10修正了该模型，提出 新的模型（YNS-LBE 模型）。该模型也在碰撞前 改变流入流体方向，但是公式中使用 i 方向的分布 函数代替i 方向的分布函数， 使其在界面处速度连 续，如式(6)。 * s ( , )( , )( ,) c iii tnftft = xxx (6) WBS-LBE 模型在 i 方向上的流体分布函数是 碰撞后的分布函数，而i 方向上的分布函数则是在 碰撞前的分布函数，导致反弹回的流体不仅沿其反 方向流体，也会向其他方向流动，引起数值不稳定 性。 ZHU 和 MA11重新分配同一时刻同一格点碰撞 后方向相反的一对流体质点，每个方向上碰撞后的 流体质点分为两部分，一部分沿该方向流向下一格 点，另一部分由于固体壁面的作用沿相反方向反弹 回去，提出新的部分反弹模型（ZM-LBE 模型） ，该 模型在同一时间的格点仅改变方向相反的一对流体 质点，因此该模型满足质量守恒，且数值稳定性比 WBS-LBE 模型好。该模型的松弛时间接近 1 时， 有效黏度系数与流体黏度系数相等。 该模型中 i (x,t) 如式(7)。 cc* s ( , )( , )( ,) ii i tnftft = xxx (7) 部分反弹模型区别主要在于如何处理反弹边 界的作用、有效黏度、平均散射密度及有效渗透 率，总结目前主要部分反弹模型见表 1。表中包 括 DM-LBE 模型、 TS-LBE 模型、 WBS-LBE 模型、 YNS-LBE 模型和 ZM-LBE 模型， 这 5 个模型的渗 透率根据每个格点的平均散射密度 ns计算得到。 前 4 个模型在 ns=0 时，均可简化为标准 LBE 模 型. 在 ns=0 时，ZM-LBE 模型不能自动满足反弹 边界，因为该模型在同一时间格点仅改变方向相 反的一对流体质点。WBS-LBE 模型、YNS-LBE 模型和 ZM-LBE 模型在 ns=1 时能够满足全反弹 LBE 模型，其余模型不满足全反弹 LBE 模型，因 为当 ns=1 时，DM-LBE 部分反弹模型在不透水边 界处有质量损失，而 TS-LBE 部分反弹格式模型 仅反映一个时间步的流体函数而不是原则上的两 个时间步。其次，GS-LBE 模型、WBS-LBE 模型、 YNS-LBE 模型和 ZM-LBE 模型能够满足质量守 恒，DM-LBE 模型和 TS-LBE 模型在高渗透性和 低渗透性多孔介质边界之间不对称而不满足质量 守恒，需要在格点的边缘设置 ns来修正模型，只 有在密度恒定或者平均散射密度不变时才满足质 量守恒。CHEN 等 22比较了 GS-LBE 模型、 DM-LBE 模型和 TS-LBE 模型，分析了在松弛时 间 =1 时， 为了保证流体黏度为正值， 3 个模型中 ns的上限值。GS-LBE 模型中 ns (1/23)0.289 DM-LBE模型中 ns (571/28)0.234，TS-LBE 第 6 期 张潇丹等：REV 尺度多孔介质格子 Boltzmann 方法的数学模型及应用的研究进展 1701 表 1 REV 尺度多孔介质内部分反弹 LBE 模型（f为孔隙流体黏度）11 模型 ( , ) i tx 有效黏度 ns（=1） 有效渗透率 GS-LBM s ( , )( , ) i i nftft xx () ss s 121(1)(21) 6 12 nn n + 0.289 f s 2n DM-LBM s (, )( , ) i ti i nftft+ xcx () ssss s 121(1)(21)(12) 6 12 nnnn n + + 0.234 f s 2n TS-LBM c*c* s (,)( ,) ti i nftft i + xcx ()()() () () s s s 212611 1 6 121 n n n + + 1 f s 2n WBS-LBM c* s ( , )( ,) i i nftft xx ()() ()() ss ss 21121 6 1211 nn nn + + 1 () sf s 1 2 n n YNS-LBM c* s ( , )( ,) ii nftft xx 1 - ZM-LBM cc* s ( , )( ,) ii nftft xx ()()() ()() 2 ss ss 268112121 6 12112 nn nn + + 0.5 () sf s 12 2 n n 模型和WBS-LBE模型中 ns在01之间。ZM-LBE 模型的 ns上限为0.5， 说明ZM-LBE模型最多有一 半的流体质点可以反弹回相反的方向。综上， YNS-LBE模型在 ns=0时，均可简化为标准LBE 模型，在 ns=1时能够满足全反弹LBE模型，能够 满足质量守恒，数值稳定性也较好，因此应用比 较好。 1.2 阻力模型 REV尺度多孔介质的另一种LBE模型是阻力 模型12-21，即多孔介质的影响通过施加外力得到 Stokes或Brinkman方 程 ， 这 个 思 想 是 基 于 NITHIARASU等23提出的求解不同孔隙度渗流的 Navier-Stokes方程，用其他CFD方法（有限差分、 有限体积）也可以方便求解，采用LBM方法是众 多CFD方法的一种。该类模型的实现途径有两类： 第一类是修正速度12-17，即将多孔介质的阻碍作用 通过修正速度来体现；另一种是流体与多孔介质构 造体-固体的相互作用力直接加入分布函数的演化 方程中18-21。 1.2.1 修正速度的阻力模型 修正速度的REV尺度多孔介质LBE阻力模型 是在碰撞步中修正瞬时速度来描述多孔介质对流体 的影响，所以其演化方程与标准BGK演化方程相 同12，如式(8)。 (eq) 1 (,)( , )( , )( , ) iittii ftftftft += xcxxx (8) 1997年SPAID和PHELAN12建立了直接求解 BRINKMAN方程的LBE模型（SP-LBE模型） 。该 模型将流体-固体的相互作用力的影响通过修正平 衡分布函数的速度来实现，这种作用力处理方式对 应的宏观流动方程中存在一些误差，忽略误差后， SP-LBE对应的宏观方程在稳态时就是Brinkman方 程。SP-LBE模型平衡态分布函数的计算采用修正 后的平衡态速度，修正的平衡态速度包含了多孔介 质的阻力作用，而阻力大小又与当地速度相关，如 式(9)式(11)。 () 2 (eq) (eq)(eq) 2 (eq) 242 sss () ( , )1 22 i i ii ft ccc =+ cu cuu x (9) (eq) F =+uu (10) = Fu (11) 式中， v/K，是渗透势能参数，是黏度 v 与 渗透率 K 的比值。 与SP-LBE模型类同，FREED13把多孔介质阻 力通过碰撞过程在平均速度中体现Freed-LBE模 型，式(13)，可以用于Mach数较大的流动过程。 与多孔介质流动控制方程相比，Freed-LBE模型对 应的宏观方程存在人工多余项，需要通过具体问题 来分析。该模型比较适合用于具有裂隙及导电介质 的岩石孔隙中流体及其他工程模拟。Freed-LBE模 型公式如式(12)式(14)。 ()2 2 (eq) 242 sss ( , )1 22 i i ii ft ccc =+ cucuu x (12) 11 = 1) + 22 F =+( uuuu (13) 化 工 进 展 2016 年第 35 卷 1702 Gu =u (14) () =Fuu (15) 式中，u为碰撞后速度， u 碰撞前速度，u 为 平均速度，G是速度张量， 取决于流体的运动黏度、 松弛因子和渗透系数。 Shan-Chen模型14-15最初是应用于多相流和多 组分流动，把体积力 F 在平衡速度中体现，体积力 F 主要包括液液流体间相互作用力、液固相间相互 作用力和外力。LI等16-17基于Shan-Chen模型提出 单相单组分LBM方法（Li-LBE模型） ，其平衡分 布函数可以简化成式(16)，平衡态速度为式(17)。该 模型与其他模型相比，不用计算体积力 F，节省计 算时间，提高计算效率。 () () (eq) eq 2 s ,1 i ii ft c =+ cu x (16) 8 (eq)0 2(12 ) ii i f = =+ c uu (17) 8 1 0.5 1 2 i it i t f K = + = + cG u (18) 式中， K 为渗透率， 为孔隙率， G 为外部体积力。 1.2.2 修正演化方程的阻力模型 修正演化方程的阻力模型是将多孔介质效应直 接加入到演化方程，其演化方程19为式(19)。 (eq) 1 (,)( , )( , )( , ) iittiitii ftftftftF += + xcxxx (19) 式中，Fi是 F 的函数，根据其处理方式不同， 可以分成不同的阻力模型。 MARTYS18对SP-LBE模 型 进 行 改 进 （Martys-LBE模型） ， 采用LUO24的作用力处理方 法描述介质阻力，在演化方程中直接增加了一个作 用力项来消除部分离散误差， 但是不能完全消除， 如式(20)。其平衡态分布函数与SP-LBE模型相同， 即式(9)，但是作用力不修正平衡态分布函数的速 度，而是直接修正演化方程。Martys-LBE模型比 SP-LBE模型的数值稳定性好，因为SP-LBE模型 修正平衡分布函数的速度，增加了非物理项。 Martys-LBE模型其作用力表达式如式(20)此时 c =Fa， i i f=。 22 ss () iii ii FF cc =+ cucu c (20) GUO等19在2002年提出了求解包含线性阻 力、 黏性项和非线性阻力项的通用渗流模型的REV 尺度渗流的广义LBE模型（GLBE模型） 。该模型 的平衡分布函数和作用力项中均中包含反映介质特 性的孔隙率 ，如式(21)和式(22)，而且对渗流速度 大小没有要求。此模型采用GUO等25提出的作用 力模型， 不存在离散误差， 可得到准确的宏观方程。 其中，F是包含介质阻力和外部体积力的总力，式 (23)右端第一和第二项是线性（Darcy）和非线性 （Forchheimer）介质阻力。F是结构参数，与孔隙 率有关，对由固体颗粒构成的多孔介质，如 公式(24)。 ()2 2 (eq) 242 sss ( ) ( , )1 22 i i ii ft ccc =+ cucu u x (21) 2 s 24 ss :()1 (1) 2 ii i ii c F cc =+ cFuFc cI (22) F KK +F =uu uG (23) 3 1.75 150 F = (24) 目前还有一种趋势，将LBM与其他数值方法 结合起来模拟多孔介质内流动和传递过程。 ZARGHAMI等20把GLBE模型和有限体积法相结 合，用有限体积方法离散对流扩散项，如式(25)， 并用校正因子来提高数值稳定性。该模型通过 Poiseuille流、Couette流和平板驱动流的验证，用 于模拟均质和非均质多孔介质管道流。 () eq neqneqneq ,1,1 , neqneqneq 1,11,1 neqneqneq 1,11,11,1 1 d 11 48 1 16 ii acbd iii x yxyx y x y iii xyx yxy iii xyxyxy ffA fff A fff fff + + + = + + + (25) 式中，Ax.y为 abcd 的面积，ac、ab、bd、bc 分 别为格子边界。 ZHANG等21提出了Navier-Stokes方程的显性格 式解的LBE模型（Zhang-LBE模型） ，该模型为通过 Chapman-Enskog和泰勒展开，省略高阶项之后，得 到平均体积Navier-Stokes方程，在演化方程式(26) 中加入了附加项 Pi，如式(28)。该模型只在平衡分布 第6期 张潇丹等：REV尺度多孔介质格子Boltzmann方法的数学模型及应用的研究进展 1703 函数中加入反应介质特性的孔隙率，如式(27)，但是 与GLBE模型中孔隙率加入方式不同， 作用力项处理 方法采用式(29)，F公式与GLBE模型一致， 见式(23)。 该模型能有效地解决多孔介质流和两相流。 (eq) 1 (,)( , )( , )( , ) ii ttiiit iti f xctfx tfx tfx tPF += + (26) ()2 2 (eq) 242 sss ( ) ( , )1 22 i i ii u ft ccc =+ cuc u x (27) 24 ss 1 (1) 2 ii iii Pp cc =+ cucu c (28) 24 ss 1 (1) 2 ii iii F cc =+ cucu c F (29) 两种类型多孔介质REV尺度LBE模型，应 用较为广泛的是Freed-LBE模型和GLBE模型。 如果考虑不可压流时，GLBE模型更加精确；对 于稳态不可压低速多孔介质流动问题，两者的区 别在于GLBE模型加入了非线性阻力项，在低速 情况下，该项的作用较弱26。GLBE模型的作用 力是基于GUO等的作用力模型， 可以准确得到宏 观方程，不存在离散误差；其次，该模型的平衡 分布函数中包含反应介质特性的孔隙率，而其余 REV尺度LBE模型的平衡分布函数与标准模型 类似，不能反映介质特性；最后，GLBE模型可 以用于高速流动。 2 REV 尺度 LBE 模型应用现状 2.1 REV尺度LBE模型应用于多孔介质内流动 各REV尺度LBE模型在多孔介质内流动中均 有运用，但是运用最广泛的是GLBE模型。 WBS-LBE模型、YNS-LBE模型和ZM-LBE模 型均被用于模拟不同渗透特性多孔介质内的流 动10-11，如图1。模拟结果显示WBS-LBE模型在 界面处不连续，YNS-LBE模型和ZM-LBE模型在 界面处速度连续，而且在 ns较小时的连续性比较好 （ns1=ns3=0.2，ns2=0.1） ，且多孔介质的孔隙度变化 不会影响中心处的速度。 REV尺度LBE模型不仅被用于分层多孔介质 内速度的连续性， 还被用于多孔介质内Klinkenberg 效应等基础研究。孔隙大小在纳米到微米范围内， 克努森数 （Kn 数， 分子平均自由程和流动的宏观特 征长度之比）相对较大时，气体不再满足连续介质 假设，导致气体滑移流、瞬变流，甚至真空分子流。 图1 不同渗透特性多孔介质示意图及其速度分布图10-11 由于存在滑移效应，测得表观渗透率高于固有渗透 率，两者之差随着 Kn 数增大而增大。CHEN等27 基于四参数随机生成方法（QSGS）并运用基本元 素建模模型（EBB模型）重新构造页岩，采用 BESKOK和KARNIADAKIS28提出的表观渗透率 公式，在 Kn 数为常数时根据固有渗透率计算表观 渗 透 率 ， 运 用GLBE模 型 研 究 了 二 维 含 有 Klinkenberg效应的多孔介质内流动。高孔隙率 （=0.8）时，Klinkenberg效应可以忽略，随着孔隙 率减小，Klinkenberg效应在低渗透多孔介质比较明 显， 并且随着 Kn 数增大、Klinkenberg效应越明显， 化 工 进 展 2016年第35卷 1704 随着压力的减小、Klinkenberg效应更加明显 （图2） 。 随后CHEN等29通过Darcy定律预测REV尺度的 整个计算区域内的有效渗透率， 研究多孔介质成分、 颗粒大小、 压力、 滑移及吸附作用对REV尺度渗透 率的影响。 在致密储气岩层或页岩层开采天然气时， 如果不考虑Klinkenberg效应， 用表观渗透率代替固 有渗透率，不仅会导致定量错误而且会导致定性 错误。 图2 含有Klinkenberg效应和不含Klinkenberg效应的 速度分布图27（1psi=6894.76Pa） 除了用于致密岩石中渗流， 土壤渗流也是REV 尺度多孔介质LBE模型的研究热点。 樊火30通过绘 制渗流场几何剖面图或者读入图形文件，格点化流 场、生成LBM格式的边界数据文件，基于GLBE 型开发了土壤渗流LBM_Seepage软件平台， 生成数 据文件，以直观可视的等值线或者云图实现渗流场 速度实时可视化，并对竖向具有不同渗透率的土体 进行了模拟。结果表明，最上层土的相对的渗透能 力在渗透初期对渗流速度的传播起控制作用，只有 经过一段时间后土体才会表现出等效应的渗透能力 （图3） 。 申林芳等31进一步运用GLBE模型探讨了 土体在压力作用下孔隙率、渗透率及渗透压力等影 响因素与渗流速度的相互关系。 Freed-LBE模型被推广到非均匀网格中32，采 用HE等33提出的插入格式描述非均匀网格，该模 型在流动和碰撞过程之间采用插值法，不仅能保证 数值准确性而且易于实施。该模型对多孔介质、孔 隙和边界的处理方式一样，避免了不同尺度采用不 同数学方法而产生的误差，与其他多尺度LB模型 相比，它消除了内部边界，而且采用阻力模型，因 此可以用于大尺度。KANG等32在裂隙周围细化网 格研究裂隙对流动的影响，远离裂隙的地方采用粗 网格（图4） ，模拟了复杂裂隙和简化裂隙的多孔介 质内的流动，通过 u l K p = 预测REV尺度多孔介质 内的有效渗透率（表2） ，并提出可以把该模型用曲 线坐标和复合网格的设想。 REV尺度LBE模型不仅可以模拟含裂隙的均 质岩石中的流动，还能根据通过高分辨率的扫描图 像确定空腔、黏土、长石、石英不同成分构成的非 均质砂岩结构，如图5(a)，基于GLBE模型获得砂 岩样品内的流体流动速度，然后根据达西定律求得 整体的渗透率34-35。由图5(b)可以看出，样品大于 图3 不同渗透率的的边界示意图及渗透率 K 由小到大和 K 由大到小中心线速度图30 第6期 张潇丹等：REV尺度多孔介质格子Boltzmann方法的数学模型及应用的研究进展 1705 图4 裂隙示意图、非均匀网格示意图、计算结果32 表2 渗透率值32（格子单位） 项目 上半部分 下半部分 裂隙 整体 传统 LBE 模型 0.64 0.64 35.31 1.94 Freed-LBE 模型 0.61 0.61 37.01 1.99 误差 4.7% 4.7% 4.8% 2.6% 10mm时整体的渗透率趋于常数，对于砂岩样品， 大于10mm时可以达到REV尺度的要求。 由图5(c) 和图5(d)可以看出，流场主要受无渗透性的石英和 空腔影响，由于各成分的孔隙率不同导致各成分边 界处的流线不闭合。但是CT扫描大孔时容易忽略 狭小的孔喉结构而造成误差，因而扫描图像的分辨 率必须要高36。 REV尺度多孔介质模型除了可以通过渗流速 度预测多孔介质的渗透率外，还可以分析多孔介质 图5 砂岩样品 几何结构特性及外部条件对渗流速度的影响。 GLBE模型被用于在压力驱动下的非牛顿流体的电 渗 流 过 程 37 ， 在 演 化 方 程 中 加 入 了 基 于 Herschel-Bulkley模型的非牛顿流体力项， 研究外部 电势驱动流、离子浓度和电场强度、固体颗粒直径 和孔隙度、幂律指数和屈服应力、德拜长度对渗透 速度的影响。多孔介质内的流速随着幂律指数和屈 化 工 进 展 2016年第35卷 1706 服应力的减小而增大，壁面处的电势主要影响壁面 处流速，多孔介质的电势对通道中心处流速影响较 大。低幂律指数时，壁面处的电势在对流速的影响 更加敏感，特别是在通道高度与德拜长度的比值较 小时。 HUANG等38对LBM两相流模型及其应用 做了详细描述， 目前LBM两相流模型已经应用于 REV尺度多孔介质内多相流的驱替过程中。 SPAID和PHELAN39对纤维多孔介质的流动进行 了计算，纤维内部使用SP-LBE模型，外部使用 标准LBE模型，模拟圆形纤维内树脂驱替空气多 组分流动过程。纤维间的流动比纤维内的流动速 度快，导致纤维内部含有空气空穴，空穴最初是 长轴垂直于流动方向的椭圆形，随后由于表面张 力的影响，从椭圆形变为圆形（图6） 。根据达西 定律预测的不饱和渗透率在定性上正确。但是需 要更多研究证明纤维内的空洞随时间消失模拟结 果是否准确。SCHAAP等 40应用阻力模型中 Shan-Chen格式多相流LBE模型，模拟了充满玻 璃珠的多孔介质内水-空气及水-石油精体系驱替 实验。该模型基于显微层析构造的多孔介质几何 结构，忽略重力、惯性力和黏性力，只考虑毛细 作用力，设置合理的表面张力和接触角，得到压 力与饱和度的关系。水-空气体系模拟值与测量值 相差不大，但是水-石油精体系模拟值与测量值差 异很大（图7） ，原因是模拟条件与实验条件不一 致，如接触角。该方法也可以用于微尺度界面现 象的模拟。PORTER等41也运用该LBE模型估算 润湿-非润湿相界面的面积，模拟值与实验值 吻合。 图6 不同时刻单行圆形纤维树脂驱替空气模拟结果图39 （黑色：空气；灰色：树脂） 图7 压力与饱和度演化图 40 阻力模型中的单松弛GLBE模型推广到多松弛 模型，LIU等42运用二维D2Q8 MRT-GLBE模型模 拟不可压多孔介质流，在低黏性时的数值稳定性比 单松弛模型好。随后该模型被用于解含有时间导数 和非线性对流项平均体积动量方程，模拟了在多孔 介质圆柱静止时上下壁面以恒定速度运动和上下壁 面静止时管道内圆柱多孔介质以恒定速度运动的两 种情况下的流动43。 结果表明R1和R2的本征相平 均速度吻合，二者的相平均速度在流体区吻合，但 是在多孔介质区域内有偏差图8(b)，并且随着孔 隙度的减小而增大。多孔介质阻力对本征相平均速 度和相平均速度有类似规律。只有含有本征相平均 速度的宏观方程满足伽利略不变性。多松弛模型与 单松弛模型相比，克服了许多限制条件，比如运动 黏度和体积黏度比不再是定值，但是多松弛模型增 加了计算量。 REV尺度LBM已经应用到工程渗流问题。何 莹松44将REV尺度GLBE模型应用于煤矿开采过 程中回采工作面瓦斯渗流问题，该研究中认为煤体 的孔隙度受应力重新分布的影响很小，孔隙度基本 不发生变化，渗透率沿流动方向呈现分段函数的不 均匀分布，得到了不同时刻非均质煤层中瓦斯压力 分布和速度分布。这些数据能够生动再现瓦斯流动 第6期 张潇丹等：REV尺度多孔介质格子Boltzmann方法的数学模型及应用的研究进展 1707 图8 多孔介质圆柱运动示意图及=0.7、Re=100通道中心 固有相平均速度、相平均速度图42 过程，有助于指导煤矿安全生产。 2.2 REV尺度LBE模型应用于多孔介质内流动与 传热耦合过程 流动、 传热和传质两个及两个以上过程耦合时， LBE模型可以分为3种；即多速模型、双分布函数 模型（DDF）与差分方法相结合的混合模型。目前 文献大多采用DDF模型来模拟多过程， 该模型使用 两个（及以上）分布函数演化方程，分别用于速度 场和温度场（浓度场） ，格子结构简单，且具有良好 的数值稳定性。 等温GLBE模型通过DDF模型推广到多孔介 质热流动中（DDF-GLBE模型）45，速度场用等温 GLBE模型模拟，而温度场使用一个新的分布函数 描述。DDF-GLBE模型被用于多孔介质方腔模拟孔 隙度恒定45-49、孔隙度变化50、含有热源51的速度 场和温度场，研究孔隙度、雷诺数（Re 数） 、达西 数（Da 数）等参数对速度场和温度场的影响。 DDF-GLBE模型还被用于研究边界充满多孔介质 和中心充满多孔介质两种情况下管道内的流动与传 热过程52-53。 REV尺度的GLBE模型可以直接将孔隙尺度模 拟得到的渗透率和有效热导率带入REV模型， 避免 了采用经验公式可能得到的错误结论54。模拟结果 表明，REV尺度得到的热壁上的平均Nusselt数与 孔隙尺度下的趋势一致， 证明REV尺度在定性上是 正确的。在瑞利数（Ra 数）低于107时，GLBE模 型能定量预测孔隙尺度下自然对流的宏观规律，Ra 数大的时候只能定性预测孔隙尺度的宏观规律（图 9） 。REV尺度GLBE模型对渗透率如何影响多孔介 质自然对流宏观现象刻画很好，与孔隙尺度模拟结 果一致，但是往往夸大了流固热导率带来的影响。 图9 不同Ra数下热壁上平均Nusselt数随方柱个数N的变 化结果对比54 虚线REV；实线孔隙尺度 REV尺度的LBE模型为太阳墙工程应用研究 提供参考， 通过含热源项的LBE模型对太阳墙系统 中复杂流场及热传递进行模拟55-56，可以预测太阳 墙在不同工作条件下的状况，并优化太阳墙系统的 设计。 针对工程中常见的多孔介质圆管内传热问题， 荣 伏 梅 等 57-58 构 造 了REV尺 度 轴 对 称 热 DDF-GLBE模型，包含弱可压轴对称LBE模型和 不可压轴对称LBE模型， 把刻画多孔介质的参数加 入到分布函数中并构造合适的分布函数，使之能恢 复正确的宏观方程，具有形式简单、外力项不含复 杂梯度形式的优点， 并应用于非线性渗流数值模拟。 分别研究了覆盖多孔介质的方柱绕流、充满多孔介 质内的变黏性Rayleigh-Benard对流、 以及圆管内插 入多孔介质强化对流传热非线性渗流问题。由于这 类问题的计算量都很大，单个CPU的计算效率较 低，因此采用GPU（图形处理器）加速计算，可支 化 工 进 展 2016年第35卷 1708 持复杂的计算任务59。 单松弛DDF-GLBE模型推广到多松弛模型中， 建立REV尺度二维多松弛（MRT）DDF-GLBE方 程，用D2Q9 MRT-GLBE模型模拟流场，同时用 D2Q5 MRT-LBE模型模拟温度场， 模拟结果和解 析解与文献计算结果吻合60，而且在低黏度时， MRT-LBE模型的数值稳定性比LBGK模型好，但 是计算量大。该直角坐标系下的MRT-DDF-GLBE 模型延伸到轴对称模型61。 2.3 REV尺度LBE模型应用于多孔介质内流动、 传质和化学反应耦合过程 与流动和传热过程类同， 仍使用DDF模型描述 流动和传质两个过程，即运用两个分布函数分别模 拟流场和浓度场，并在浓度演化方程中加入反应源 项， 用于模拟化学反应。TIAN等62基于DDF-GLBE 模型模拟了CO2地下填埋地球化学反应过程，研究 CO2注入过程中溶质离子反应及扩散、岩石骨架 CaCO3溶解相互作用， 引起孔隙率和渗透率变化， 反作用于流场。如图10，开始在一段时间内出口处 的H+、Ca2+和HCO3离子浓度为0，随着反应的进 行，3种离子浓度逐渐增加。由于反应消耗和对流 扩散的相互作用，在出口处反应物H+离子浓度在 1105时出现变化。Ca2+离子浓度到达最大后保持平 衡，是因为Ca2+离子是可逆反应。 该模型还用于研究有机废水溶液中化学反应、 燃料电池和质子交换膜燃料电池中56 ，63-66。LIAO 等63和杨艳霞64研究了有机废水中含扰流光合细 菌包埋颗粒的流动、传质及其内部的光生化反应过 程，分析光照强度、进口流速、包埋颗粒的渗透率 及孔隙率对流场、浓度场及产氢性能的影响。XU 等65-66把该模型应用于固体燃料电池多组分化学反 应，计算浓差极化电势，并研究了电池电极厚度、 电解质浓度、孔隙率、燃料组成等结构参数和工作 图10 出口离子浓度图62 温度操作参数对浓差极化、超电势、电化学效率、 电池电压和功率密度的影响，LBM模拟结果与其他 模拟结果相比更接近实验值。该模型也被用于模拟 质子交换膜