相间力对微孔通道内两相驱替渗流影响规律的LBM研究

张一夫，梁 冰，张遵国

（1.辽宁工程技术大学 安全科学与工程学院，辽宁 葫芦岛 125105；2.矿山热动力灾害与防治教育部重点实验室,辽宁 葫芦岛 125105；3.辽宁工程技术大学，辽宁 阜新 123000）

孔隙介质尺度下的细观多相渗流的研究具有十分广泛的工程背景[1-2]，如煤层注水、石油天然气开采、核废料处理、有机物对土壤的污染范围等都与细观多相渗流的特点和机理有密切关系。对于孔隙尺度的研究，因其涉及到复杂的物理力学机制，对其内部现象的描述较为困难。既要考虑到不同流体之间的相互作用，还要考虑到流体与孔隙骨架之间的流固耦合，特别是当渗流通道达到了微孔尺度，许多宏观尺度下的作用力，如重力、惯性力等不再是影响流动特点的主要因素，而像表面张力、毛细力及孔隙介质表面的润湿机理等作用更加凸显，成为影响两相渗流规律的主要因素。对两相流动的数值研究，传统的方法是依托Navier-stokes方程并结合对界面追踪方程的求解来追踪两相界面运动，常见的方法有流体体积函数法（VOF）[3]、水平集方法（Level set）[4]、Front tracking[5]和 phase-field[6]等。这些方法通常是用来描述少数大的界面运动，但在孔隙尺度下，由于物理模型的复杂性，往往很难刻画出流动的细节和大量细小、分散界面的准确位置。格子Bohzmann方法(LBM)是近二十几年发展起来的一种基于分子动力学的模拟方法[7]，具有清晰的物理背景。由于其微观本质和介观特点，可以方便地描述出流动中不同相间的相互作用，特别是边界条件易于实施的特点，在模拟复杂流动问题上具有天然的优势。近年来提出的基于动理学理论的多组分流动LBM模型主要有颜色模型[8]、伪势模型[9]、自由能模型[10]等，其中伪势模型因其能方便地描述相间作用力、自动追踪相界面且计算程序简单易行，在多相流领域得到充分发展、应用最为广泛。利用LBM方法在孔隙尺度上对微孔介质内的流动进行数值模拟，拟通过改变相间力的大小和作用机制，找到反映微孔通道内两相渗流的运动规律和特点。

1 基于多组分伪势模型的LBM方法

取经典的二维九速正方格子-玻尔兹曼模型，简称D2Q9模型，对格子速度eα进行离散。采用单松弛时间的LBM多组分伪势模型对多孔介质内的两相驱替过程进行研究。D2Q9离散速度模型如图1。

图1 D2Q9离散速度模型

式中：α为9个格子速度分量的方向；c为格子速度。

1.1 考虑相间力的多组分伪势模型

在多组分伪势模型中，对多孔介质中两相介质分别采用不同的分布函数进行描述。考虑外力项且碰撞项采用线性化（LBGK近似）的第k种组分粒子分布函数f的演化方程（格子-Boltzmann方程）为：

其表达式为：

式中：u、ρ分别为混合流体的宏观速度和密度；cs为格子声速，取；wα为权重系数，在D2Q9模型中其值分别为w0=4/9，w1-4=1/9和w5-8=1/36。

每种组分流体的宏观密度和速度由下式确定：

混合流体的速度为：

作用在第k种组分上的总作用力Fk包括3个部分：流体组分之间的相互作用力、流体组分与固体壁面之间的作用力以及流体组分所受的重力。

伪势模型通过建立组分流体间的相互作用势，得到来自于相邻组分流体的作用力：

式中：ψk为第k种组分“有效密度”；Gkk为格林（Green）函数。

式中：g为第k种组分和第种组分之间相互作用的强度。

当存在固-液界面时，为了考虑固体表面润湿性对流体流动的影响时，固体表面对第k种组分流体的作用力可写成：

式中：s为固相标识函数，其值取为1为固相，取0为流体相；Gkw为控制流体和固体表面相互作用力强度的参数，通过调节Gkw值，可以得到不同的表面润湿性条件，其值取负代表为润湿，取正代表为非润湿。

在D2Q9模型中，Gkw值为：

当有重力存在时，可用下式表示：

式中：g为重力加速度。

当所有外力都给定时，混合流体的压强p可由系统的状态给定：

1.2 复杂曲线边界处理

当研究孔隙尺度下的多相渗流问题时，必须要考虑固体骨架的具体形态对流动的影响。由于LBM方法是通过碰撞和迁移过程来描述流动粒子运动，因此对于流体粒子而言，在固体骨架附近的粒子的迁移和碰撞过程结果的尤为重要。曲线边界处节点类型示意图如图2。

图2 曲线边界处节点类型示意图

粒子碰撞仅发生在流体节点自身所在的位置，而迁移过程需沿着网格线方向以格子速度将相应的分布函数传递到相邻节点。曲线边界处粒子的迁移过程如图3。

图3 曲线边界处粒子的迁移过程

图3中xf位置的分布函数在碰撞后需沿着1到8 个格子方向分别以 e1、e2、e3、...、e8的速度进行迁移。以发生在5方向xf点与xs点之间的迁移过程为例，由边界节点xw上的分布函数（xw，t）运动到流体节点xf的表达式为：

因此，对于直角正交网格下的LBM模型，在每一个碰撞-迁移循环中，迁移前均需利用插值获得与物理边界相交的网格线上的分布函数值（xw，t），插值结果即为此处的曲线边界条件，据此插值分布函数值（xw，t）进入式（12）的迁移过程。

为了获得从边界进入流场内部的分布函数值，采用Bouzidi等[11]提出的一种结合反弹格式和空间插值的边界处理方法。Bouzidi格式如图4。

图4 Bouzidi格式

该方法通过其引入表征边界位置与相邻固体节点、流体节点间相对距离的线性插值因子q，从而可以更加精准地描述固体边界的位置。其表达式为：

利用二阶线性插值格式获得固体边界附近流体节点的分布函数f值：

当q<1/2时：

当q≥1/2时：

1.3 数值方法验证

1.3.1 液滴静态接触角的验证

为了验证流固作用参数Gw对壁面润湿性情况以及液滴静态接触角的影响，建立了液滴-固壁物理模型。将初始半径为20的液滴置于100×100的二维格子系统中，液滴初始圆心置于（50,0），上下边界为固体壁面，采用半步长反弹边界条件，左右为周期性边界条件。2 种组分流体的参数分别为 ρ1，in＝ρ2，out＝1.0，v1，in＝v2，out＝0.16，组分间作用力参数 Gf为 2.4。壁面润湿性及稳态壁面接触角模拟结果如图5。

图5 壁面润湿性及稳态壁面接触角模拟结果

从图5的模拟结果来看：流固作用强度参数Gw数值的大小决定了液滴的稳定形态。随着Gw由0.5变化到－0.5，壁面由疏水性转变为亲水性，液滴与壁面间的静态接触角θ由大逐渐变小，且呈近似线性规律变化（图6）。从而验证了伪势模型对模拟两相流固耦合作用的敏感性和准确性。

图6 壁面接触角θ与流固作用强度参数Gw之间的关系

1.3.2 毛细力液桥现象的验证

为进一步验证流固作用力及两相间表面张力模型施加的正确性，建立不同流固作用参数Gw下液桥的稳定性模型。将初始宽度为10的液柱置于100×100的二维格子系统的中心位置，上下边界为固体壁面，采用半步长反弹边界条件，左右为周期性边界条件。2 种组分流体的参数分别设置为，ρ1，in＝ρ2，out＝1.0，v1，in＝v2，out＝0.16，组分间作用力参数 Gf为 2.7。不同流固作用强度参数Gw时，液桥或液滴形态的演变过程如图7。

图7 不同流固作用强度参数Gw时，液桥或液滴形态的演变过程（粗线为最终稳定状态）

从图7的模拟结果来看，流固作用参数强度Gw对稳态液柱形态有显著影响，相间作用力及流固作用力是决定能否形成稳定液桥的主要影响因素。当相间作用力（相间力作用参数Gf）不变时，流固作用力过大（图 7（a））或者过小（图 7（d））都很难形成稳定的液桥，只有在合适的流固作用力下（图7（b）和图7（c））才能形成稳定的液桥。

2 相间力对多孔介质内的两相驱替规律的影响

为探讨相间力对孔隙介质两相渗流的影响，建立500×100个格子单位的规则孔隙结构，上下边界为固体壁面，采用半步长反弹边界条件，左右Zuo-He边界条件，入口流速为0.1。2种组分流体的参数分别设置为 ρ1，in＝ρ2，out＝1.0，v1，in＝v2，out＝0.16，组分间作用力参数Gf为2.1。分别取不同的Gw值，对多孔介质内两相驱替过程的进行了数值模拟，得到了两相界面演变规律（图8）。

从图8中可以看出，随着两相流体的驱替过程向右推进到孔隙介质区域时，相界面开始发生扭曲、破碎和变形。随着时间的推移，原本连续的相界面变成了大量细小、分散的界面。当孔隙介质表面设置为中性润湿性界面（Gw取0）时，驱替流体前缘界面向前推进时，后续流体由于受孔隙介质骨架的阻碍，推进较慢，导致两相流体的相界面被拉得很长、黏性指进现象明显、驱替效率较低。

当孔隙介质表面的亲水性逐渐增强时（Gw=-0.5和 Gw=-1.0）（图 9 和图 10），随着驱替流体相界面的前缘向前推进，后续驱替流体依靠相间力和流固作用力，加快其在孔隙骨架表面上的润湿速度，驱替流体的整体推进速度较快、黏性指进现象减弱、驱替效果提高。

图9 Gw取-0.5时两相界面的时空演变规律

图10 Gw取-1.0时两相界面的时空演变规律

另外，从两相界面的总体驱进速度对比来看（图11），随着煤体润湿性的增强（Gw减小），两相界面的总体驱进速度也得到了提高。

由此可见，相间力对微孔通道内两相渗流驱替过程影响较为明显，随着孔隙介质表面润湿性的增强，两相驱替的总体效果增强，驱替效果较好。因此，为了提高煤层的注水效果，可以从改善煤体表面润湿性角度出发，通过注入液体改性来提高煤体表面的润湿性，特别是低渗透煤层，由于煤体渗透性差，孔隙盲孔、角隅较多，容易形成封闭气，通过改性能使在在孔隙通道内两相驱替过程中表面张力、毛细力的作用更加凸显，明显改善注水效果和驱替速度。

图11 Gw取不同值时，两相界面的总体驱进速度对比

3 结论

1）基于伪势模型的多相流格子波尔滋蔓方法可以准确地反映出相间力对微孔通道内两相渗流躯替过程力学特性的影响规律。

2）从模拟结果的分析来看，两相驱替的总体效果与相间力有较大关系，多孔骨架的亲水性越强，两相渗流的总体驱替效果就越好，从而有利于驱替速度的提高，反之亦然。

3）可以通过改进注入液体在煤体孔隙表面的润湿性，减少低渗透煤层注水过程中封闭气体的形成，从而有效提高注水效果。