利用格子Boltzmann方法计算页岩渗透率

张 磊,姚 军,孙 海,孙致学

(中国石油大学石油工程学院,山东青岛266580)

页岩气作为非常规油气资源日益受到重视，成为常规能源的替代能源，Javadpour等[1]对北美9个油藏152块岩心进行实验分析，发现页岩的基岩渗透率主要在5.4×10-8μm2左右，90%的页岩渗透率小于1.5×10-7μm2。由于页岩的低孔低渗等特性，实验室研究其物理特性难度较大。格子Boltzmann(简称LB)方法是近十几年来发展起来的介观模拟方法，可以对岩心的SEM图像或CT图像经过简单二值化[2]处理后直接进行模拟；学者们[3-5]利用LB方法计算多孔介质的渗透率，并与实验结果和经验公式进行对比，拟合结果比较好。同时，Javadpour等[1]指出页岩气的主要流动孔隙直径在4～200 nm，页岩的孔隙直径从纳米级到微米级，气体在纳米孔隙中的流动不同于达西流动，气体分子可以自由地在孔隙表面滑动以及与壁面或者其他分子发生碰撞，因此利用LB方法模拟页岩中的流动必须考虑页岩孔隙的微尺度效应。笔者利用Knudsen数修正的LB模型进行模拟，计算页岩的绝对渗透率。

1 格子Boltzmann方法

格子Boltzmann模型一般由格子模型、平衡态分布函数和演化方程组成，Qian等[6]提出的DdQm(d为空间维数，m为离散速度个数)模型是格子Boltzmann方法的基本模型。格子Boltzmann方法的计算过程由碰撞和迁移两部分组成，常用的Bhatnagar-Gross-Krook(BGK)模型(不包含外力项)的方程为

式中，α=0，1，2，…，N，表示共有N个速度方向；r为粒子的空间位置；eα为α方向的速度；t为时间；δ为时间步长；fα为离散速度空间α方向上的分布函数；为离散速度空间的局部平衡态分布函数；ρ和u分别为宏观的粒子密度和速度；cs为格子声速；ωα为权系数。

在每个格子点上满足以下关系：

1.1 基本模型

应用最广泛的模型是D2Q9和D3Q19模型，其模型的速度方向示意图见图1。在格子边长和时间步长均取1的情况下，D2Q9模型各方向速度分布为

D3Q19模型的速度分布：

格子声速和权系数分别为

图1 D2Q9和D3Q19模型示意图Fig.1 Schematic diagrams of D2Q9 and D3Q19 models

1.2 松弛时间和边界处理

利用格子Boltzmann方法模拟微尺度流动时，须解决两个基本问题，松弛时间τ与Knudsen数(Kn)的关系和边界条件的处理。对松弛时间τ与Kn的关系进行修正[7]：

为了描述边界上的滑移流动，边界条件不再采用标准的反弹格式，而是采用标准反弹与镜面反弹混合的反射格式，见图2。

式(7)中的α是标准反弹和镜面反弹的比例系数，当α=1时，是标准反弹格式；当α=0时，是完全镜面反弹格式；当α=0.5时，是理想漫反射格式。边界采取这种简单的处理方式是为了在三维实际数字岩心中便于实现。

图2 粒子边界反弹示意图Fig.2 Schematic diagram of particle boundary bounceback

1.3 Kn的计算

Knudsen数定义为Kn=λ/H(λ为气体的分子平均自由程，H为流场的特征长度)，由于利用格子Boltzmann方法模拟流动时涉及到实际物理单位与格子单位之间的转化，同时Kn为无量纲数值，在两种单位下应该保持不变。

在格子单位下，Kn和渗透率等参数主要依赖于格子的划分，与格子分辨率有关，同时存在一些不变量，如格子声速cs在D2Q9和D3Q19模型中都为1/，即。由动理学理论可知，气体的平均自由程λ与动力黏性系数μ、压力p和温度T之间关系为

式中，v为运动黏性系数；N为特征尺度方向上划分的格子数；R为摩尔气体常数。

Kn的表达式为

体外癫痫模型海马神经元中线粒体衔接蛋白Miro1表达的改变及其意义 … 王英，连亚军，谢南昌，等 130

当选定松弛时间τ和格子速度c，并且格子划分确定后，就可以得到Kn，但是在实际中，Kn与温度、压力等有关，对于硬球分子，其分子平均自由程为

式中，m为分子质量，g；σ为分子直径，m；M为摩尔分子质量，g/mol-1。

Kn的表达式为

式中，Na为阿伏伽德罗常数，Na=6.022×1023mol-1。

对于页岩气，主要成分为甲烷，分子直径σ约为0.414×10-9m，这样给定温度和压力以及流动通道的特征尺寸，就可以计算出Kn。图3给出在温度为350 K，不同特征长度下Kn与压力p的关系。在利用格子Boltzmann方法模拟页岩中流体流动时，计算得到Kn后，根据上式对格子模型中的参数N、δx、δ和τ进行设置，也可以利用式(6)对τ直接进行修正。

图3 温度T=350 K和不同特征长度H下Kn与压力p的关系Fig.3 Relation between Kn and p with different H at T=350 K

2 渗透率计算

2.1 平板模型

模型采用二维平板间的缝隙流动，其渗透率理论值为l2/12，l为两平板之间的距离，见图4。

图4 平板模型示意图Fig.4 Schematic diagram of flat model

利用格子Boltzmann方法计算其渗透率，出、入口边界条件不采用周期性边界，而是采用更符合实际情况的压力边界，不考虑微尺度流动的影响，固体边界采用标准反弹格式，二维平板模型的理论渗透率与LB计算渗透率的对比见图5。由图5可以看出，LB模拟的渗透率与理论渗透率拟合很好。

图5 LB渗透率与理论渗透率对比Fig.5 Comparison between LB permeability and theoretical permeability

图6 考虑Kn影响的渗透率曲线Fig.6 Permeability curve considering Kn effect

在考察新模型中Kn对计算结果的影响时，模型中采取在x和y方向格子数分别取200和20，在不考虑Kn影响的一般模型中，通过计算得到格子单位下渗透率K=33.08509。图6给出的是在不同Kn值影响下计算得到修正后的渗透率K′及两者的相对误差。从图6中可以看出，随着Kn增大，修正后的渗透率K′越小，以相对误差5%为界限，在相对误差超过5%时必须考虑Kn的影响，此时对应的Kn约在0.14；假如给定温度为350 K，压力为2 MPa，此时可以算出对应的孔隙尺寸约为22.6 nm，也就是对于实际的岩心，在该温度、压力条件下，当孔隙喉道尺寸小于22.6 nm时，必须考虑Kn对渗透率的影响。

采用5对平板并排的模型(图7)研究Kn对通道内流动的影响。平板之间的宽度l分别为10、20、30、40 和50 nm，温度T=250 K，压力p=2 MPa，利用式(11)计算得到其Kn依次为：0.317 1、0.158 55、0.1057、0.079 275和0.063 42，流体流动方向由左向右。图8为在格子时间t=800时，在x=L/2处，x方向上的速度分布，其中L为平板长度，图中黑线为考虑Kn修正的模型计算得到的x方向速度U′x，红线为不考虑Kn的一般模型计算得到的x方向速度Ux。从图8中可以看出，随着宽度增加，这5条通道中流体流动的速度最大值显然是逐渐增大的。对于同一通道，考虑Knudsen数Kn的影响后，通道中间部分格子点x方向的速度变大，而且随着Kn增大(宽度减小)，U′x和Ux之间的差别逐渐增大。图中右上角小图是把虚线框中50 nm通道中的曲线放大，可以看出在靠近边界的部分，考虑Kn影响的x方向的速度是变小的，在数据上是靠近边界的第一个格子点的速度变小，其他格子点的速度都是变大的，这在其他几个通道中也是同样的结果。Kn的增大意味着分子平均自由程增大，在相同压差条件下流动更快。

图7 多平板模型示意图Fig.7 Schematic diagram of multi-flat model

图8 x=L/2处的x方向速度分布Fig.8 x-direction velocity distribution at x=L/2

2.2 页岩渗透率计算

图9 页岩SEM扫描图像Fig.9 SEM image of shale

对四川盆地某地区页岩岩心进行SEM电镜扫描，得到分辨率为2 μm的图像(图9)，由于页岩的致密性，孔隙连通性较差，因此不对其做二维的平面LB流动模拟，通过数值重构方法得到平面岩心切片三维图像的方法主要有高斯模拟法[8]、模拟退火法[9]、过程模拟法[10]、多点统计法[11]和基于马尔科夫链的蒙特卡洛重建法[12]。利用文献[13]中随机建模的方法，将页岩二维的SEM图像重构成三维图像，然后利用经过Knudsen数修正的格子Boltzmann模型对三维图像进行模拟。由于进行SEM扫描时，选取的是页岩中孔隙比较发育的部分，得到数字岩心的孔隙度为14.999 7%，计算出的渗透率代表孔隙比较发育部分的渗透率。图10为流体在数字岩心中流动达到平衡状态时的速度分布。由图10可以看出，流体流速在绝大部分孔隙中处于相对低速，只有在很少的大孔隙中流速相对较高，图中速度为格子单位下的速度，通过计算得到该页岩的绝对渗透率为0.0185 μm2。

图10 页岩数字岩心流体速度分布场Fig.10 Velocity distribution field in shale digital core

3 结束语

利用格子Boltzmann模拟页岩中的流体流动，Knudsen越大，得到的修正渗透率与原渗透率差别越大。在已知压力和温度，一定相对误差标准下，可以计算出模型必须考虑Knudsen效应的孔隙最大宽度。Knudsen效应使得通道中间部分流体速度增大，而在边界附近流体速度会减小。用格子Boltzmann模拟方法计算渗透率相对于试验室物理实验方法有很大的便捷性。

[1] JAVADPOUR F，FISHER D.UNSWORTH M.Nanoscale gas flow in shale gas sediments[J].Journal of Canadian Petroleum Technology，2007，46(10)：55-61.

[2] 赵秀才，姚军，房克荣.合理分割岩心微观结构图像的新方法[J].中国石油大学学报：自然科学版，2009，33(1)：64-67.ZHAO Xiu-cai，YAO Jun，FANG Ke-rong.A new reasonable segmentation method for microstructure image of reservoir rock[J].Journal of China University of Petroleum(Edition of Natural Science)，2009，33(1)：64-67.

[3] 朱益华，陶果，方伟，等.3D多孔介质渗透率的格子Boltzmann模拟[J].测井技术，2008，32(1)：25-28.ZHU Yi-hua，TAO Guo，FANG Wei，et al.Lattice Boltzmann simulation of permeability in 3D porous medium[J].Well Logging Technology，2008，32(1)：25-28.

[4] DEGRUYTER W，BURGISSER A，BACHMANNET O，et al.Synchrotron X-ray microtomography and lattice Boltzmann simulations of gas flow through volcanic pumices[J].Geosphere，2010，6(5)：470-481.

[5] JIN G，PATZEK T W，SILIN D B.Direct prediction of the absolute permeability of unconsolidated and consolidated reservoir rock[R].SPE 90084，2004.

[6] QIAN Y，DHUMIERES D，LALLEMAND P.Lattice BGK models for Navier-Stokes equation[J].Europhysics Letters，1992，17(6)：479.

[7] GUO Z，ZHAO T，SHI Y.Physical symmetry，spatial accuracy，and relaxation time of the lattice Boltzmann equation for microgas flows[J].Journal of Applied Physics，2006，99，074903.

[8] JOSHI M.A class of stochastic models for porous media[D].Lawrence：Universtiy of Kansas，1974.

[9] HAZLETT R D.Statistical characterization and stochastic modeling of pore networks in relation to fluid flow[J].Mathematical Geology，1997，29(6)：801-822.

[10] BRYANT S，BLUNT M.Prediction of relative permeability in simple porous media[J].Physical Review A，1992，46(4)：2004-2012.

[11] OKABE H，BLUNT M J.Prediction of permeability for porous media reconstructed using multiple-point statistics[J].Physical Review E，1997，55：1959-1978.

[12] WU K，VAN DIJKE M I J，COUPLES G D，et al.3D stochastic modelling of heterogeneous porous media-applications to reservoir rocks[J].Transport in Porous Media，2006，65(3)：443-467.

[13] 赵秀才.数字岩心及孔隙网络模型重构方法研究[D].青岛：中国石油大学石油工程学院，2009.ZHAO Xiu-cai.Numerical rock construction and pore network extraction[D].Qingdao：School of Petroleum Engineering in China University of Petroleum，2009.