页岩气在纳微孔道中的流动模拟研究

曾彦，宁正福*，齐荣荣，黄亮，吕朝辉

1 中国石油大学(北京)石油工程教育部重点实验室, 北京 102249

2 中国石油大学(北京)油气资源与工程国家重点实验室, 北京 102249

页岩气在纳微孔道中的流动模拟研究

曾彦1,2，宁正福1,2*，齐荣荣1,2，黄亮1,2，吕朝辉1,2

1 中国石油大学(北京)石油工程教育部重点实验室, 北京 102249

2 中国石油大学(北京)油气资源与工程国家重点实验室, 北京 102249

富有机质页岩中存在大量的纳米孔隙，连续介质流的Navier-Stokes方程难以准确描述纳米孔隙中的气体流动行为。为此，采用格子Boltzmann方法模拟了气体在页岩纳微孔道中的流动。结果表明：页岩气在纳微孔道中的流动存在显著的纳米尺度现象，表现为气体的平均流速大于Poiseuille公式计算值；纳微孔道中气体的稀薄效应显著，孔壁处的气体速度不为零；且孔道直径越小，孔壁的滑移速度越大，甚至大于孔道中的气体速度、流速剖面也由抛物线型转变为柱塞状。页岩气在纳米孔中流动时发生“双滑脱效应”，流体粒子与孔壁发生碰撞，获得动能后反弹进入孔道中，增强整个纳微孔道气体的流动。“双滑脱效应”导致页岩气的渗透率大于等效液体渗透率，并引起克氏渗透率图中，气体渗透率与平均压力的倒数偏离直线关系。

页岩；纳米孔；格子Boltzmann方法；流动模拟；纳米尺度；双滑脱效应

富有机质页岩中聚集了大量的天然气，这使得页岩气成为全球油气资源勘探开发的热点[1-2]。前人已经采用多种技术手段来研究页岩储层的孔隙结构特征[3-6]。采用离子抛光制样和高分辨率的场发射扫描电镜实验显示页岩分散的有机质中存在大量的纳米孔隙，孔径大小为10～800 nm。尽管页岩中存在一些大孔隙，高压压汞和低温氮气吸附实验结果表明页岩主体孔径集中在50 nm以内。气体在纳米尺度空间中的流动规律不同于宏观尺度流动。在纳米尺度下，气固界面之间的分子间作用力更强，动量交换更加剧烈，甚至成为气体流动的主要影响因素[7]。

国内外学者对气体在页岩储层中的渗流行为进行了大量的研究。Javadpour[8-9]建立了页岩气在纳米孔隙中的传质模型，考虑了页岩气的扩散、滑脱和吸附解吸行为，通过引入视渗透率参数，修正了页岩气渗流方程。Freeman[10]利用尘气模型，研究了Knudsen扩散对低渗透致密介质流动的影响。王瑞等[11]考虑页岩气扩散类型的差异，分析了扩散对页岩气质量流量的贡献。这些都侧重于宏观渗流的研究，而对于页岩气在纳米孔隙中的微观渗流行为还鲜有报道。本文引入适合宏观和介观尺度的格子Boltzmann方法，考虑流体粒子之间、流体粒子与纳米孔壁之间的相互作用，从微观角度研究页岩气在纳微孔道中的流动行为。

1 页岩气在纳米孔隙中的流动状态

1.1 页岩中的纳米孔隙

图1为页岩有机质中的纳米孔隙。页岩样品来自四川盆地下寒武统牛蹄塘组以及西北地区六盘山盆地下白垩统乃家河组，有机碳含量为0.5%～9.15%，成熟度为0.6%～2.7%。页岩样品首先采用氩离子抛光技术进行处理，然后在FEITMQuantaTM200F场发射扫描电镜下观察孔隙结构，实验仪器的相关参数和实验方法见文献[5]。扫描电镜结果显示，页岩有机质孔隙直径为10～450 nm。

为了和可视化的电镜图片结果进行对比，对页岩样品进行了高压压汞实验。压汞实验的驱替压力最高达到200 MPa，测试孔隙直径下限为7.35 nm。压汞孔径分布曲线表明，页岩样品的孔隙直径主要为8～100 nm(图2)。

图1 页岩有机质中的纳米孔隙Fig. 1 Nanopores in organic matters

图2 页岩压汞法孔径分布曲线Fig. 2 Pore size distribution of shales

1.2 流动状态

图3 Knudsen数与气体流动状态Fig. 3 The relationship between Knudsen number and fl ow state

根据Knudsen数(气体平均分子自由程与流动空间特征尺度的比值)[12]的大小，将流体在多孔介质中的流动状态划分为(图3)：(1)当Kn＜0.001时，为连续介质流。流体可以认为是连续介质，气体渗流服从达西定律。(2)当0.001≤Kn＜0.1时，为滑脱流。气体分子之间的碰撞频率仍然较高，但是稀薄气体效应已不能再忽略，这时流动控制方程仍然是经典流体力学的Navier-Stokes方程，但必须考虑边界滑移速度和非等温条件下的温度跳跃，气体流动服从Klinkenberg方程。(3)当0.1≤Kn＜10时，为过渡流(或Knudsen流)。此时气体分子之间的碰撞频率和气体分子与固体之间的碰撞频率差不多，连续介质假设已不再成立，气体流动服从Burnett方程。(4)当Kn≥10时，为自由分子流。此时气体分子与固体之间的碰撞频率比分子之间的碰撞频率高得多，气体分子之间的碰撞可以忽略，主要采用分子动力学方法进行研究。

表1为典型页岩气藏的气体组成和气体性质参数。页岩气主要由CH4、C2H6和CO2组成，其摩尔分数分别为80%、10%和10%。图4为温度为300 K时，随着储层压力的变化(0.1～40 MPa)，页岩气在不同孔径中的Knudsen数。页岩的孔径分布主要为8～150 nm，因此从页岩储层到井底，页岩气在纳米孔隙中的流动主要为滑脱流和弱过渡流。在这两个流动阶段，气体分子的平均自由程大、密度小，稀薄效应明显，气体流动违背了控制方程中的连续介质假设，不存在滑脱边界的Navier-Stokes方程是无效的。因此，基于达西定律的经典渗流理论不适用于页岩气储层流动模拟。

2 格子Boltzmann方法

基于动力学理论的Boltzmann方程是连接宏观和微观的桥梁，可用于从连续流到自由分子流的跨尺度流动。求解Boltzmann方程非常困难，从Boltzmann方程和介观层次上的格子气方法发展而来的格子Boltzmann方法(LBM)成为一个较好的解决方案。格子Boltzmann方法将流体颗粒离散成一系列的流体粒子。这些粒子比分子级别大，但在宏观上又无限小，流体粒子在2D离散格子网格组成的毛细管模型中碰撞和流动。在一定的初始和边界条件下，这些流体粒子之间、以及流体粒子和毛细管孔壁之间相互作用，但它们共同地遵循Boltzman流动方程。以气体动力学理论为基础的格子Boltzmann方法已经被认为是模拟微观尺度介质的有效数值方法[12-13]。

表1 典型页岩气藏的气体组成和气体性质参数Table 1 Components properties of shale gas

图4 不同孔径下，页岩气在储层中的Knudsen数Fig. 4 The relationship between pressure and Knudsen number under different pore sizes

2.1 计算模型

Freeman等运用尘气模型(DGM)研究页岩储层中非吸附相的运移。DGM可以准确表征滑脱效应和扩散效应[14]。将DGM控制方程离散，引入格子Boltzmann方法，可以预测页岩气在多孔介质中的视渗透率。基于Maxewll-Stefan扩散方程，页岩气的流动方程可以表示为[15]：

式中，P是压力(Pa)；R是普氏气体常数，8.314 J/(mol·K)；T是温度(K)；μ是动力黏度(kg/(m·s))；Dk是Knudsen扩散系数(m2/s)。K0是固有渗透率，K0=d232，m2

Knudsen扩散系数可表示为：

式中，d是孔道直径，m；M是气体分子量，kg/mol。对于三维多孔介质模型，沿流动方向的有效Knudsen扩散系数可表示为：

式中，Lx，Ly，Lz分别为多孔介质x，y，z三个方向的长度，其中x方向为流动方向。ρin，ρout分别为进出口端密度。引入视渗透率Kapp，则有：

2.2 格子Boltzmann模型

本文采用格子Boltzmann-BGK方程D2Q9模型(二维九速模型)对纳微孔道(如图5)中的页岩气流动进行了模拟。

图5 纳米孔道物理模型Fig. 5 Physical model of nanopores

通过BGK近似、时空离散和速度离散，得到完全离散化的格子Boltzmann-BGK方程[12-13]，即LBM演化方程：

式中，δt表示时间步长；fα(r,t)表示t时刻r处速度为eα的粒子密度分布函数；为对应时刻点的局部平衡态分布函数；τ为松弛时间；Fα代表由孔道和外力项引起的体积力，可表示为：

其中，G为外力扰动项，Kapp为式(5)求得的视渗透率。Guo和Zhao将平衡态分布函数和外力项分布函数分别表示为：

式中，ρ、u表示流体宏观密度和流动速度；ωα为与离散速度方向矢量模量有关的权系数；cs为格子声速；eα表示α方向上的粒子速度。

模型的速度配置如下：

由式(7)可知，外力项F中也包含了速度u，将式(7)代入式(12)即可求得速度u。

在给定的离散空间和离散时间下，式(6)中的松弛时间计算公式为[17]：

式中，Ny表示y方向上的网格数，Kn为给定空间的Knudsen数。式(14)中Knudsen数与松弛时间的关系为[18-19]。

式中，p0为出口端压力；Kn0为出口端Knudsen数。

2.3 边界条件

在LBM模拟中，对于固体边界的处理方法主要有两种：无滑移边界的反弹格式和自由滑移边界的镜面反射格式。对于气体在页岩纳微孔道中的流动，既不能用简单的反弹格式处理，也不能用镜面反射来描述气体与固体壁面的相互作用与动量交换。将反弹和镜面反射相结合能比较准确地反映真实气体与固体之间的相互作用。在混合格式中，引入弹回比例系数r，表示粒子与边界作用后沿原路返回所占的比例。研究表明[20]，r取值为0.7时，能最佳地描述气体与固体壁面的相互作用。本次模拟，设定上下边界为封闭边界，采用反弹和镜面反射的组合格式，且弹回比例系数为0.7；左右边界为定压边界，左边界为入口，右边界为出口，孔道中流体在两端压差作用下进行流动。

定压边界采用的是外推边界条件。假设物理边界外还有一层虚拟节点，物理边界节点(i,1)和这层虚拟节点(i,0)上的分布函数都参与碰撞迁移过程。借鉴有限差分的概念，在边界节点(i,1)采用中心差分，确定物理边界节点上的未知分布函数。

外推计算之后，对所有节点执行迁移步骤。通过分布函数计算进出口端压力，并对其进行修正。使进出口端平均压力分别等于给定的压力值Pin，Pout。在随后的碰撞过程中，边界节点上的平衡态分布函数需要由修正后的压力值来计算。对压力进行修正，不仅要保证进出口端压力满足给定的压力条件，还使得其压力剖面符合内部流场特征。

(3)语法复习目标：掌握各项基础语法的基本概念及重要用法；熟练理解、掌握并能够在不同的语境中辨别、正确、灵活运用各种重要语法功能。

3 模型验证

应用LBM之前，必须对该方法和模型进行验证。LBM适用于连续介质流、滑脱流和弱过渡流。本文首先在连续介质流条件下，通过LBM的模拟结果与Poiseuille流动的解析解进行对比，说明LBM在连续介质流下的正确性。在连续性假设失效时，通过与他人的解析结果和数值结果进行对比，验证在纳米尺度下LBM的正确性。

3.1 连续介质区的流体流动验证

物理模型如图5，孔道上下为封闭边界，相距H，孔道中充满了黏度为µ的流体，在压力梯度驱动下，流体沿着x方向进行流动。流场空间划分为正方形均匀网格，网格数为4000×400，空间步长Δx=Δy，时间步长Δt=1，松弛时间τ=1(连续介质流不存在压缩效应，孔道两端压力一定时，沿孔道的Kn数、松弛时间为定值)，各单位均采用格子单位。在孔道两端给定一定的密度，使流体在x方向产生压力梯度。本次模拟中，进口密度为1.0012 kg/m3，出口密度为1 kg/m3。

图6 连续介质流的格子Boltzmann模拟结果与Poiseuille公式的解析解对比Fig. 6 Streamwise velocity at outlet

图6为连续介质流的出口端速度剖面。连续介质流条件下，由于流体的黏性耦合特征，孔道中流体的流动速度剖面为抛物线，即孔道中间的流体速度最大，孔壁处流体速度为零。从图中可以看出，格子Boltzmann方法的模拟结果与Poiseuille公式的解析解吻合非常好，这表明格子Boltzmann方法可以准确地描述连续介质流流体的动力学行为。

3.2 滑脱流和弱过渡流的流体流动验证

对于连续性假设失效或宏观方程难以描述的现象，采用微观分子动力学或介观气体动力学描述更加恰当。为了和他人的结果[21-24]进行对比，计算了孔道出口端的流速分布剖面。模型的基本参数选取如下：x方向网格数Nx=1100；y方向网格数Ny=11；出口端密度ρout=1；图7(a)为滑脱流阶段的模拟结果，其中出口端Knudsen数Kn0= 0.0194，入口端密度ρin=1.4；图7(b)为弱过渡流阶段的模拟结果，其中出口端Knudsen数Kn0= 0. 194，入口端密度ρin=2。

通过比较发现，笔者的模拟结果与Frederik Verhaeghe的模拟结果、Arkilic理论解以及DSMC方法的模拟结果重合性良好，验证了模型及程序的正确性，说明在滑移区和弱过渡区，采用LBM方法可以有效模拟气体在页岩纳微孔道中的气体流动行为。

4 模拟结果及讨论

下面模拟了页岩气在不同孔径纳微孔道中的流动。纳米孔孔径分别为10 nm、20 nm、50 nm、100 nm和1 μm，孔道的径长比为0.01，孔道的平均压力为0.1、0.5、1、2、5和10 MPa，储层温度300 K，页岩气藏的气体组成和气体性质参数见表1，计算得到页岩气在纳微孔道中流动时对应的Knudsen数为0.000 58～5.826 0。

4.1 平均流速

为了便于和经典流体力学的Poiseuille流动结果进行对比，定义纳微孔道中气体流动的无因次平均流速为：式中，为采用LBM方法模拟的纳微孔道的平均流速；采用Poiseuille公式求得流动速度的解析解。

模拟得到的不同直径孔道的无因次平均流速与平均压力的关系如图8。从图中可以看出，高压下，无因次平均流速趋于1，纳米孔的平均流速与经典Poiseuille公式的解析解一致。但当压力＜5 MPa，模拟得到的纳微孔道流速大于Poiseuille公式的解析解。压力低于2 MPa时，纳微孔道的平均流速比经典值高1～2个数量级。纳微孔道直径越小，流体流动速度和Poiseuille公式解析解的差别越大。低压下气体在纳米孔的平均流速的模拟结果显著高于Poiseuille流动解，尽管我们可以将这归因于“滑脱效应”的影响，但如此高的平均流速是非常少见的，分析认为这是由于纳微孔道中的“双滑脱效应”的促进作用，将在下文中进行解释。

图7 LBM方法模拟滑脱流和过渡流阶段的孔道出口端速度剖面与文献结果对比Fig. 7 Normalized streamwise velocity at different Knudsen number

图8 不同直径纳米孔的气体无因次平均流速Fig. 8 Dimensionless average velocity in different nanopore sizes

4.2 空间流速分布

随着压力降低，气体分子的平均自由程增大，Knudsen数增大，纳微孔道中气体的流动与Poiseuille流的差别更明显。图9和图10为平均压力为0.1 MPa时，页岩气分别在1 μm和50 nm纳米孔中的空间流速分布(无因次流速为流动速度与孔道出口端中心处流速的比值)。沿着流动方向(x方向)，纳米孔中心处(y/H=0)的流动速度越来越大。这是因为越靠近出口端，孔道内的压力越低，气体密度越小，Knudsen数越大，气体稀薄效应越严重。为了满足质量守恒，气体的流动速度增大。图中孔道的上下边界处(y/H=±0.5)表现出明显的滑脱效应，边界流速不为零。孔径为1 μm的(图9)边界滑移速度小于孔道内气体流速，但孔径为50 nm的纳米孔(图10)边界处的流速与孔道中气体的速度相当，这是由于孔径减小，气体与孔壁之间的作用力增强的结果。

图9 页岩气在1 μm孔中流动的空间速度分布Fig. 9 Dimensionless velocity distribution in 1 μm pore

图10 页岩气在50 nm孔中流动的空间速度分布Fig. 10 Dimensionless velocity distribution in 50 nm pore

4.3 流速剖面

图11为模拟得到的不同孔径纳米孔的出口端流速剖面。纳微孔道直径为1 μm(或者大于1 μm)时，尽管边界处的流速不为零，但出口端的流速剖面仍然为Poiseuille流的经典抛物线形状。随着孔道直径的降低，由于气体的滑脱效应增强，孔壁附近的滑移流速逐渐增大。当孔道直径减小至50 nm时，尽管孔道中心速度存在极值，但孔壁滑移速度与孔道中心速度相当，孔道出口端的流速剖面趋于均匀。孔径小于20 nm后，孔壁的滑移速度显著大于孔道中心速度，流动速度剖面变为柱塞状。低压下，气体在小尺寸纳米孔(＜50 nm)的流动行为与在大孔(1 μm或更大)中的流动出现根本区别。

为了进一步说明纳米孔的气体流动与大孔道中气体流动的区别，我们计算了出口端气体流动的动能剖面(图12和图13)。无因次动能为LBM模拟的动能与Poiseuille公式计算的出口端孔道中心(x/L=1，y/H=0)动能的比值。从图中可以看出，孔道为1 μm时，LBM的模拟结果与Poiseuille公式计算的动能(解析法)非常接近，孔道中心处的动能几乎一致；只是在孔壁处，考虑滑脱效应的LBM方法模拟的动能略大于解析结果(图12)。但当孔道直径为50 nm时，不仅孔壁处的动能大于解析法结果，而且孔道中心的动能也显著高于解析法结果(图13)。因此，随着孔径的降低，纳米孔孔壁处流体粒子的动能增大，并且孔道中心处气体的动能也增大。这是由于流体粒子与孔壁之间的相互作用引起的。气体在孔道中流动时主要发生两种基本相互作用：流体粒子间的碰撞，流体粒子与孔壁的碰撞[25]。当流体粒子的平均自由程小于孔隙直径时，常发生流体粒子间的碰撞；而对于那些平均自由程大于孔隙直径的流体粒子，则常发生流体粒子与孔隙壁间的碰撞。随着孔道直径的降低，流体粒子平均自由程大于孔隙直径，流体粒子与孔壁的碰撞越剧烈，滑脱现象愈显著。孔径小于20 nm时，流体粒子与孔壁的相互作用逐渐占主导地位，成为气体流动的主要影响因素，边界滑脱速度显著高于孔道中心的气体流速，流动速度剖面因此转变为“两端高中间平缓”的柱塞状。

图11 纳米孔的出口端流速剖面Fig. 11 Streamwise velocity at different nanopore sizes

图12 直径为1 μm孔道的出口端动能剖面Fig. 12 Kinetic energy distribution in 1 μm pore

由于孔壁处的流体粒子流速不为零，孔壁处的流体粒子与孔壁发生非弹性碰撞后获得了动能，并且反弹后的流体粒子在垂直于孔壁方向上的速度不为零。获得动能的流体粒子进入孔道中，使孔道中气体动能增大，流动速度增大，增强了孔道中心气体的流动(图14)。正是由于这个原因，从出口端气体流动的动能剖面上可以看出，随着孔径的减小，孔道中心(y/H=0)的气体动能大于Poiseuille公式计算的动能。因此，不仅管壁处的流体粒子发生滑脱，整个孔道中的流体也受到这种滑脱效应的促进作用，即所谓的“双滑脱”现象[26-28]。1941年Klinkenberg[29]首次提出滑脱效应，传统的滑脱效应主要用于表征滑脱区的流动特征，只考虑了孔壁处的气体流速不为零，没有考虑到孔壁对流体粒子在整个纳微孔道中流动的促进作用。随着流动尺度不断缩小，传统的滑脱效应越来越难以准确表征流体视渗透率。很多学者[30-31]选择更高阶的Klinkenberg方程来表征滑脱效应。页岩气储层主要孔径为几十个纳米，脱离了传统滑脱效应表征范围，“双滑脱”现象对于纳米孔中页岩气的流动促进作用非常重要。

图13 直径为50 nm孔道的出口端动能剖面Fig. 13 Kinetic energy distribution in 50 nm pore

图14 “双滑脱效应”中页岩气体分子与纳米孔孔壁的相互作用Fig. 14 The interaction between gas molecules and nanopore

4.4 孔道渗透率

页岩气体在纳微孔道中流动时，流体粒子与孔道之间存在强烈的相互作用。纳米尺度下气体的流动行为与常规尺度渗流有所不同。为了研究纳微孔道中，纳米尺度流动对页岩气渗流的影响，定义无因次渗透率为：

式中，KLBM为采用LBM模拟的不同平均压力下的纳米孔的气体渗透率；KPoiseuille为采用Poiseuille公式和达西定律计算的等效液体渗透率。

在经典的滑脱效应理论中，考虑气体滑脱的气测渗透率与孔道平均压力的倒数呈直线关系，或者气测渗透率和等效液体渗透率的比值与平均压力的倒数为直线。图15为页岩气在不同孔径纳米孔孔的无因次渗透率与平均压力倒数的关系(克氏渗透率校正图)。从图中可以看出，高压下(10 MPa)，大孔(孔径为1 000 nm)中的气测渗透率与等效液体渗透率相等，表现为无因次渗透率为1的直线。但在低压下，“双滑脱效应”对纳米孔的气体渗流产生了显著的影响：纳微孔道的气体渗透率大于等效液体渗透率，并且在克氏渗透率校正图中，孔道的气体渗透率与平均压力的倒数不再呈直线关系。随着孔道平均压力的降低，无因次渗透率与平均压力倒数的关系曲线发生上翘；纳微孔道的直径越小，无因次渗透率的上翘现象越明显。平均压力为2 MPa时，直径为50 nm孔的气体渗透率是等效液体渗透率的2.1倍，直径为20 nm孔的气体渗透率是液体渗透率的3.7倍。平均压力为1 MPa时，直径为20 nm孔的气体渗透率则是液体渗透率的9.8倍。“双滑脱”现象增强了页岩气在纳米孔隙中的渗流，并且孔径越小、压力越低，这种增强作用越明显。在页岩气开发时，从地层到井底，气体的流动压力逐渐降低，“双滑脱”现象在井底附近影响最大。因此，在室内渗透率的实验测试和实际页岩气开发中，应当考虑“双滑脱效应”对气体渗流的影响。

图15 页岩气在纳米孔中的无因次渗透率与平均压力倒数的关系Fig. 15 The relationship between dimensionless permeability and the reciprocal of the average pressure

为了研究页岩气在多孔介质中的渗流，基于扫描电镜图片，用MCMC方法构建出数字岩心。数字岩心尺寸为200×200×200，孔隙度为0.1，上方为入口端，下方为出口端。z方向上采用压力边界条件，x方向和y方向采用周期性边界条件。运用格子Boltzmann方法进行流动模拟，研究滑脱效应对视渗透率的影响。当流动达到稳定时，运用式(5)计算出数字岩心视渗透率和固有渗透率。

图16Kapp/K0和Kapp与有效流动半径的关系Fig. 16 Variation ofKappandKapp/K0with mean pore radius

图16表征了视渗透率(Kapp)和视渗透率同固有渗透率之比(Kapp/K0)与有效流动半径(r)之间的关系。模拟过程中，多孔介质平均压力分别为5，10，15 MPa，气体动力黏度为1.447×105Pa·s，有效流动半径变化区间为1～1000 nm。结果表明，当压力一定时，随着r的增大，Kapp/K0不断减小，Kapp不断增大。这是因为随着r的增大，气体流动逐渐变为连续流，滑脱效应逐渐减弱，视渗透率无限趋近于固有渗透率，Kapp/K0逐渐趋近于1。此外，随着r的增大Kapp的增幅逐渐减小。这是因为r较小时，Knudsen数较大，气体分子自由程与r处于相同数量级，r微小的变化，就能显著影响滑脱效应强度，导致Kapp较大的增幅。当r一定时，视渗透率随着压力的减小不断增大。但是当r足够大时，不同压力下的Kapp近似相等，这是因为气体流动空间足够大时，气体流动变为连续流，视渗透率近似等于固有渗透率。

为了进一步表征流动空间平均压力与视渗透率之间的关系，分别模拟了有效流动半径为5 nm，50 nm，500 nm时，不同平均压力下的视渗透率，压力变化区间为0.1～50 MPa，如图17所示。由模拟结果可知，当有效流动半径一定时，固有渗透率保持不变，平均压力越大，视渗透率越小。低压下的视渗透率远大于固有渗透率，并且随着压力的增大，视渗透率逐渐接近固有渗透率。当P=0.1 MPa，小孔径(5 nm)下气体视渗透率为6.13×10-3mD，固有渗透率为6.13×10-4mD，视渗透率是固有渗透率的10倍；大孔径(500 nm)下气体视渗透率为1.25 mD，固有渗透率为0.61 mD，视渗透率是固有渗透率的2倍。由此可知，随着有效流动空间的增大，视渗透率与固有渗透率之间的差异在不断缩小。

图17Kapp和K0与平均压力的关系Fig. 17 Variation ofKappandK0with average pressure

5 结论

(1)气体在页岩纳微孔道中的流动属于滑移区和弱过渡区范围，采用格子Boltzmann方法可以有效模拟页岩气的微观流动机理。

(2)页岩气在纳微孔道中的平均流速大于Poiseuille公式的解析解。压力越低，纳微孔道直径越小，气体流动速度与Poiseuille公式解析解的差别越大。

(3)随着页岩气体流动通道直径的降低，孔壁处的气体滑脱速度逐渐增大。孔径小于20 nm时，流体粒子的平均自由程远大于孔隙直径，流体粒子与孔壁的相互作用占主导地位，边界滑脱速度显著高于孔道中心的气体流速，孔道出口端的流速剖面也由Poiseuille流的经典抛物线形状转变为“两端高中间平缓”的柱塞状。

(4)由于“双滑脱效应”的影响，页岩气在纳微孔道中的渗透率大于等效液体渗透率；并且孔径越小、压力越低，这种增强作用越显著。同时，在克氏渗透率校正图中，气体渗透率与平均压力的倒数偏离了直线关系，发生上翘。

(5)压力一定时，随着流动空间的增大，视渗透率与固有渗透率的比值不断减小并逐渐趋于1。当流动空间一定时，随着压力的增大，视渗透率不断减小并无限趋近于固有渗透率。

[1] 邹才能, 杨智, 陶士振, 等. 纳米油气与源储共生型油气聚集[J]. 石油勘探与开发, 2012, 39(1): 13-26.[ZOU C N, YANG Z, TAO S Z, et al. Nano-hydrocarbon and the accumulation in coexisting source and reservoir[J]. Petroleum Exploration and Development, 2012, 39(1): 13-26.]

[2] 梁超, 姜在兴, 杨镱婷, 等. 四川盆地五峰组—龙马溪组页岩岩相及储集空间特征[J]. 石油勘探与开发, 2012, 39(6): 691-697. [LIANG C, JIANG Z X, YANG Y T, et al. Characteristics of shale lithofacies and reservoir space of the Wufeng-Longmaxi Formation, Sichuan Basin[J]. Petroleum Exploration and Development, 2012, 39(6): 691-697.]

[3] LOUCKS R G, REED R M, RUPPEL S C, et al. Morphology, genesis, and distribution of nanometer-scale pores in siliceous mudstones of the Mississippian Barnett shale[J]. Journal of Sedimentary Research, 2009, 79(12): 848-861.

[4] SONDERGELD C H, AMBROSE R J, RAI C S, et al. Micro-structural studies of gas shales[R]. SPE 131771, 2010.

[5] 杨峰, 宁正福, 胡昌蓬, 等. 页岩储层微观孔隙结构特征[J]. 石油学报, 2013, 34(2): 301-311.[YANG F, NING Z F, HU C P, et al. Characterization of microscopic pore structures in shale reservoirs[J]. Acta Petrolei Sinica, 2013, 34(2): 301-311.]

[6] 杨峰, 宁正福, 张世栋, 等. 基于氮气吸附实验的页岩孔隙结构表征[J]. 天然气工业, 2013, 33(4): 135-140.[YANG F, NING Z F, ZHANG S D, et al. Characterization of pore structures in shales through nitrogen adsorption experiment[J]. Natural Gas Industry, 2013, 33(4): 135-140.]

[7] 陶然, 权晓波, 徐建中. 微尺度流动研究中的几个问题[J]. 工程热物理学报, 2001, 22(5): 576-577.[TAO R, QUAN X B, XU J Z. Several questions in research of micro scale fl ow[J]. Journal of Engineering Thermophysics, 2001, 22(5): 576-577.]

[8] JAVADPOUR F, FISHER D, UNSWORTH M. Nanoscale gas flow in shale gas sediments and siltstones[J]. Journal of Canadian Petroleum Technology, 2007, 46(10): 55-61.

[9] JAVADPOUR F. Nanopores and apparent permeability of gas flow in mudstones (shales and siltstones)[J]. Journal of Canadian Petroleum Technology, 2009, 48(8): 16-21.

[10] FREEMAN C.M. A numerical study of microscale fl ow behave in tight gas and shale gas reservoir systems[R]. SPE 141125, 2010.

[11] 王瑞, 张宁生，刘晓娟, 等. 页岩气扩散系数和视渗透率的计算与分析[J]. 西北大学学报(自然科学版), 2013, 43(1): 75-80. [WANG R, ZHANG N S, LIU X J, et al. The calculation and analysis of diffusion coef fi cient and apparent permeability of shale gas[J]. Journal of Northwest University (Natural Science Edition), 2013, 43(1): 75-80.]

[12] 何雅玲, 王勇, 李庆. 格子Boltzmann方法的理论及应用[M]. 北京: 科学出版社, 2008.[HE Y L, WANG Y, LI Q. Lattice Boltzmann Method: Theory and applications[M]. Beijing: Science Press, 2008.]

[13] 郭照立, 郑楚光. 格子Boltzmann方法的原理及应用[M]. 北京：科学出版社, 2008. [GUO Z L, ZHENG C G. Theory and applications of Lattice Boltzmann Method[M]. Beijing: Science Press, 2008.]

[14] 曾彦,宁正福,王庆,等. 考虑微尺度效应的页岩气藏压裂水平井压力动态分析[J]. 中国科技论文, 2016, (23): 2679-2687.[ZENG Y, NING Z F, WANG Q, et al. Transient fl ow analysis of multiple fractured horizontal wells in shal gas reservoirs with consideration of microscale effect[J]. China Sciencepaper, 2016, (23): 2679-2687.]

[15] MASON E A, MALINAUSKAS A. Gas transport in porous media: The dusty-gas model[M]. Elsevier Science Ltd, 1983.

[16] GUO Z, ZHAO T S. Lattice Boltzmann model for incompressible fl ows through porous media[J]. Physical Review E, 2002, 66(3): 036304.

[17] ZHANG Y H, QIN R S, DAVID R E. Lattice Boltzmann simulation of rare fi ed gas fl ow in micro channels[J], Physical Review E, 2005, 71(4): 1-4.

[18] DONGARI N, AGRAWAL A, AGRAWAL A. Analytical solution of gaseous slip fl ow in long microchannels[J]. International Journal of Heat and Mass Transfer, 2007, 50(17-18): 3411-3421.

[19] LIM C Y, SHU C, NIU X D, et al. Application of Lattice Boltzmann Method to simulation microchannel fl ows[J]. Physics of Fluids, 2002, 14(7): 2299-2308.

[20] TANG G H, TAO W Q, HE Y L. Lattice Boltzmann Method for simulating gas flow in micro-channels[J]. International Journal of Modern Physics C, 2001, 15(2): 335-347.

[21] VERHAEGHE F, LUO L S, BLANPAIN B. Lattice Boltzmann modeling of microchannel flow in slip flow regime[J]. Journal of Computational Physics, 2009, 228(1): 147-157.

[22] ARKILIC E B. Measurement of the mass flow and tangential momentum accommodation coefficient in silicon micromachined channels[D]. Cambridge: Massachusetts Institue of Technology, 1997.

[23] ARKILIC E B, SCHMIDT M A. BREUER K S. Gaseous slip flow in long microchannels[J]. Journal of Microelectromechanical Systems, 1997, 6(2): 167-178.

[24] SHEN C, TIAN D B, XIE C, et al. Examination of the LBM in simulation of microchannel fl ow in transitional regime[J]. Microscale Thermophysical Engineering, 2004, 8(4): 423-432.

[25] 陈代珣. 渗流气体滑脱现象与渗透率变化的关系[J]. 力学学报, 2002, 34(1): 96-100. [CHEN D X. Gas slippage phenomenon and change of permeability when gas fl ows in tight porous media[J]. Acta Mechanica Sinica, 2002, 34(1): 96-100.]

[26] AKKUTLU I Y, FATHI E. Multiscale gas transport in shales with local kerogen heterogeneities[J]. SPE J, 2012, 17(4): 1001-1011.

[27] SWAMI V, CLARKSON C R, SETTARI A. Non-darcy flow in shale nanopores: Do we have a final answer? [R]. SPE Canadian Unconventional Resources Conference. Society of Petroleum Engineers, 2012.

[28] FATHI E, TINNI A, AKKUTLU I Y. Correction to Klinkenberg slip theory for gas fl ow in nano-capillaries[J]. International Journal of Coal Geology, 2012, 103(1): 51-59.

[29] KLINKENBERG L J. The permeability of porous media to liquid and gases[C]. New York: API Drilling and Production Practices: 200-231.

[30] HELLER R, ZOBACK M. Laboratory measurements of matrix permeability and slippage enhanced permeability in gas shales[C]. Unconventional Resources Technology Conference (URTEC), 2013.

[31] WANG J, CHEN L, KANG Q, et al. Apparent permeability prediction of organic shale with generalized Lattice Boltzmann model considering surface diffusion effect[J]. Fuel, 2016, 181: 478-490.

Simulation of transport of shale gas through the nanopores of shales

ZENG Yan1,2, NING Zhengfu1,2, QI Rongrong1,2, HUANG Liang1,2, LV Chaohui1,2
1 MOE Key Laboratory of Petroleum Engineering, China University of Petroleum-Beijing, Beijing 102249, China
2 State Key Laboratory of Petroleum Resources and Engineering, China University of Petroleum-Beijing, Beijing 102249, China

Abundant nanopores were observed in gas shales using field emission scanning electron microscopy. The Navier-Stokes equation based on the continuum assumption breaks down at the nanometer scale. The transport mechanisms of gas in nanopores of shales have been simulated using the Lattice Boltzmann Method. The simulation result shows that the microscale phenomenon of gas in nanopores of shales leads to bigger average fl ow rates than that from the classic Poiseuille equation. The rarefaction effect is very remarkable, and the slip velocity on the wall boundary is not equal to zero. The smaller the pore diameter is, the larger the slip velocity is, so much so that the slip velocity becomes bigger than gas velocity in the nanopores, and the parabola fl ow velocity pro fi le transforms from parabola to plunger. The “double skip effect” which results from the rebounding gas molecules with kinetic energy entering into the nanopores enhances the fl ow of gas in the whole nanopores, which bring about increased actual permeability, and the permeability of nanopores is bigger than that from liquid permeability. In addition, the “double skip effect” leads to nonlinear deviation from the classic Klinkenberg slip theory, especially in low fl ow pressure and smaller pore diameter.

shale; nanopore; Lattice Boltzmann Method; fl ow simulation; nanoscale; double skip effect

10.3969/j.issn.2096-1693.2017.01.007

(编辑 马桂霞)

*通信作者, nzf@cup.edu.cn

2016-08-19

国家自然科学基金“页岩气储层吸附解吸机理研究”(51274214)、教育部科学技术研究重大计划“页岩气流动机理与产能预测模型研究”(311008)、油气资源与探测国家重点实验室自主研究课题“页岩气吸附与解吸机理研究”(PRP/indep-3-1108)联合资助。

曾彦, 宁正福, 齐荣荣, 黄亮, 吕朝辉. 页岩气在纳微孔道中的流动模拟研究. 石油科学通报, 2017, 01: 64-75

ZENG Yan, NING Zhengfu, QI Rongrong, HUANG Liang, LV Chaohui. Simulation of transport of shale gas through the nanopores of shales. Petroleum Science Bulletin, 2017, 01: 64-75. doi: 10.3969/j.issn.2096-1693.2017.01.007