考虑动态滑脱的页岩微裂隙表观渗透率演化机制

丁云娜 ，李波波,2,3 ，成巧耘 ，段淑蕾 ，宋浩晟

（1.贵州大学 矿业学院, 贵州 贵阳 550025；2.贵州大学 喀斯特地区优势矿产资源高效利用国家地方联合工程实验室, 贵州 贵阳 550025；3.贵州省非金属矿产资源综合利用重点实验室, 贵州 贵阳 550025）

0 引 言

我国页岩气储量十分丰富，预计储量达36.1 万亿m3[1]。页岩气属于非常规天然气，较常规天然气储层、成藏条件及赋存特征具有明显的特殊性。在原位环境下，储层结构复杂且致密，具有低渗[2]、气体传输机制多样和渗流相关参数难以精确测量等特性[3]。其中，渗透率是决定页岩气采收率的关键因素之一，而渗透率的大小主要由渗流通道(孔隙和裂隙)决定。在页岩气开采过程中，储层压力减小，有效应力增大，微裂隙被压缩导致渗透率减小。此外，随孔隙压力逐渐减小，滑脱效应显著，对渗透率影响较大。因此，探究多因素耦合下的微裂隙渗流机制对于页岩气开采具有重要意义。

页岩气开采过程中，有效应力及气体吸附行为对固有渗透率影响较为突出，是国内外学者广为研究的热点。孔隙压力增大导致气体吸附量增大，其诱发的基质膨胀变形逐步增大，导致渗流通道减小，进一步对渗透率造成影响。孔隙压力对固有渗透率影响方面：张宏学等[4]通过气体吸附/解吸等多物理场耦合作用，得出随气体解吸，页岩基质尺寸逐渐减小，进一步诱导裂隙孔隙率增大，渗透率增大。LI 等[5]发现渗透率随着孔隙压力增大而减小，其中气体(CO2、CH4、N2)吸附膨胀主导了该变化规律。页岩气持续生产过程中，孔隙压力减小，从而改变页岩储层有效应力。BHANDARI 等[6]通过各向异性渗透率实验证明页岩渗透率与有效应力有关。BUSTIN等[7]研究表明页岩固有渗透率对有效应力表现出较强依赖性，渗透率随有效应力增大而显著减小。FAN 等[8]发现增压过程中，裂隙被压密，有效裂隙宽度减小导致渗透率减小。

大量试验表明表观渗透率与固有渗透率具有一定的差异，主要由于致密岩石内存在滑脱、扩散等多种流态，导致表观渗透率大于固有渗透率[9]。因此，部分研究学者通过修正滑脱边界条件，将表观渗透率视为固有渗透率和滑脱系数的函数。KLINKENBERG 等[10]认为在微孔和中孔中，气体分子与壁面的碰撞导致边界处的气体滑脱速度不为0。因此，引入滑脱系数修正固有渗透率，建立表观渗透率模型。但当克努森数大于1 时，传统一阶边界条件下的Klinkenberg 模型高估了气体滑脱速度，导致渗透率精确度不高[9]。鉴于此，HSIA 等[11]基于Maxwell 模型，给出精度更高的二阶滑脱模型。HATAMI 等[12]考虑气体在孔隙壁面流动时非恒定速度梯度的影响，进一步修正二阶滑脱表观渗透率模型，并发现裂隙宽度较小时，表观渗透率随气体压力的变化规律与经典的Klinkenberg 理论不一致。LYU 等[13]为研究渗透率与应力的相关性，基于双流速动力学的Klinkenberg 二阶修正方程，提出一种改进的弹性圆管模型。此外，在现有研究基础上，相关学者基于分形理论对页岩表观渗透率进行研究。分形维数是表征孔裂隙结构变形的有效手段，可描述多孔介质裂隙分布的非均质性，进而有效评估并量化孔隙结构特征与页岩储层渗透特性的相关机制[14]。LI 等[15]耦合孔隙力学、气体吸附变形和滑脱效应，采用孔隙分形维数表征裂隙宽度，建立二阶滑脱边界条件下的表观渗透率模型。

综上所述，以往研究大多将一阶、二阶滑脱系数视为恒定值，鲜有对页岩滑脱效应的动态分析。然而，滑脱系数通常随着微裂隙宽度、孔隙压力和固有渗透率的改变而动态变化，导致页岩渗透率演化机理极其复杂。因此，笔者基于立方体模型，考虑动态滑脱系数与固有渗透率间的函数关系，定量分析了一阶、二阶滑脱系数的演化规律，并探讨有效应力和气体吸附对页岩气体滑脱系数的作用机制。同时，考虑微裂隙结构变化及滑脱效应的影响，结合分形理论建立微裂隙变形和滑脱效应综合作用的表观渗透率模型，并利用公开发表的试验数据进行验证。此外，分析了微裂隙压缩性系数、分形维数和内膨胀系数对一阶、二阶滑脱系数的影响，以期加深对页岩气开采过程中渗流特性的认识。

1 页岩表观渗透率模型

1.1 页岩微裂隙分形特征

通常，可将页岩视为由基质和裂隙构成的双孔介质[16]。其中，页岩微观结构以纳米尺度为主，微裂隙形状、类型及其分布特征复杂无序，导致气体流动状态多样、传输机制复杂，其传输机制主要包括滑脱流动、克努森扩散、表面扩散[15]，如图1 所示。而传统几何计算方法较为简化，未能量化微裂隙尺寸、形状及分布规律，进而影响微裂隙中气体运移规律研究的准确性[17]。而多孔介质具有统计自相似性分形特征，分形理论是研究微裂隙结构及气体传输特性的有效手段[15]，因此，笔者基于分形理论探讨微裂隙渗流机制的演变。

图1 页岩微裂隙内部气体传输机理Fig.1 Mechanism of gas transport in shale microfractures

页岩微裂隙中微裂隙总数N与微裂隙宽度b间的关系可用分形标度表征[15]为

式中：L为微裂隙长度，m；bmax为微裂隙最大宽度，m；bmin为微裂隙最小宽度，m；Df为微裂隙分形维数，1≤Df≤2，其定义[17]为

式中：dE为欧几里得维数，二维空间里dE=2； ϕ为孔隙率，%。

由于微裂隙数目较大，可认为式(1)是连续可微函数，对式(1)求b的微分得

式中：-dN>0 为微裂隙宽度增大导致微裂隙数目减少，与实际相符。

1.2 页岩微裂隙动态演化

1) 气体吸附诱导的微裂隙变形。气体吸附对渗透率的影响不可忽略。随孔隙压力逐渐增大，气体吸附量增大而诱导基质膨胀变形增大，导致微裂隙宽度减小，进一步诱发页岩微裂隙渗透率呈减小的变化趋势。气体吸附诱导的基质体积应变变化量ΔεS满足Langmuir 型曲线，可表征[15]为

式中：ΔεS为气体吸附诱导的基质体积应变变化量；εL为Langmuir 体积应变；PL为Langmuir 压力，MPa；P为孔隙压力，MPa；P0为初始孔隙压力，MPa。

假设页岩储层表现为各向同性，此时单个方向所产生的吸附变形为

研究表明，页岩储层中微裂隙并未完全贯穿基质，气体吸附诱导的基质变形仅部分作用于微裂隙[19]。因此，引入内膨胀系数f量化基质膨胀对微裂隙变形的贡献。假设压缩为负，则气体吸附诱导的微裂隙变形为

将式(5)代入式(6)，可得气体吸附诱导微裂隙应变变化量为

式中： ϕ0为初始孔隙率，%，表示页岩微裂隙体积占单元总体积的百分比，ϕ0=3b0/a0[19]。

2) 有效应力诱导的微裂隙变形。对于多孔介质，有效应力诱导的微裂隙应变[20]可表示为

3)有效应力与气体吸附综合作用的微裂隙宽度变化。综合考虑有效应力与气体吸附的影响，页岩微裂隙有效宽度beff表示为

式中：b0为初始状态下微裂隙宽度，m；Δεf为有效应力与气体吸附引起的应变量。

其中，微裂隙应变为有效应力与气体吸附诱导[18]的应变之和

将式(7)、式(8)和式(10)代入式(9)，可得有效应力和气体吸附综合作用下页岩微裂隙宽度

式中：Kf=1/Cf[14]，Cf为微裂隙压缩性系数，MPa-1，Cf表示力作用于微裂隙时的难易程度。

由式(11)可得，分形集合内考虑有效应力和气体吸附影响时微裂隙最大宽度和最小宽度分别为

此外，微裂隙弯曲长度随有效应力的增大而增大。微裂隙在有效应力影响下实际长度[21]为

式中：v为泊松比；E为弹性模量，MPa；L0为微裂隙初始长度，m。

其中，有效线性应变[22]可表示为

式中：Δεe为有效应力作用下的应变变化量。

实际上，岩石中的有效体积应变可由气体吸附变形和应力导致的弹性变形组成[23]。结合页岩的内膨胀机制，有效线性应变可表示为体积应变与吸附应变的函数：

式中：Δεv为体积应变变化量，Δεv=-Δσeff/K[24]；K为体积模量，MPa，K=E/[3(1-2v)]。

将式(14)和式(15)代入式(13)，可得有效应力和气体吸附影响下有效微裂隙长度

1.3 基于页岩微裂隙动态演化的表观渗透率模型

微裂隙中气体的流动状态与克努森数的大小相关，克努森数为气体分子平均自由程λ与微裂隙宽度beff的比值[25]

式中：Kn为克努森数；λ为气体分子平均自由程。根据克努森数的取值范围，可分为4 个流动区域：Kn≤0.001 适用于连续介质理论，连续介质区；0.001≤Kn≤0.1，气体分子与微裂隙壁面的碰撞不可忽视，滑移区；0.1≤Kn≤10 时连续介质不适用，过渡区；Kn>10 为自由分子区，如图2 所示。

图2 页岩不同传输流态示意[25]Fig.2 Schematic diagram of different transport flow patterns of shale[25]

同时，考虑真实气体效应，微裂隙中气体平均分子自由程λ[15]可表示为

式中：R为气体常数，J/(mol·K)；T为温度，K；M为气体摩尔质量，kg/mol；μg为气体黏度系数，Pa·s；Z为气体压缩性。

其中，μg、Z[26]可表示为

式中：μg0为理想气体黏度系数，Pa·s；Tr为气体对比温度；Pr为气体对比压力；Pg为气体临界压力，MPa；Tg为气体临界温度，K；Y1、Y2、Y3均为拟合常数。

对Hagen-Poisseuille 公式[26]进行二阶修正，并考虑滑脱效应和介质变形，得到单条微裂隙中的气体体积流量q为

式中：C1为一阶滑脱系数；C2为二阶滑脱系数。

由于微裂隙动态变化，式(3)可重写为

根据质量守恒定律，对式(22)在最小微裂隙宽度到最大微裂隙度范围内进行积分，得到页岩微裂隙总体体积流量[15]为

将式(16)—式(18)以及式(22)、(23)代入式(24)中，可得有效应力与气体吸附影响下页岩微裂隙的有效气体流量为

式中：δ=beff,min/beff,max，为微裂隙最小宽度与最大宽度的比值。

由达西定律可得页岩表观渗透率为

式中：A微裂隙横截面积，且A=hb0，m2；h为微裂隙高度，且h=ξb0，m，ξ为纵横比(无因次)。

表观渗透率是滑脱系数与固有渗透率的函数。当不考虑滑脱效应时，得到的表观渗透率值即为固有渗透率，是页岩固有属性。在该条件下，不考虑有效应力及气体吸附的影响时渗透率为初始固有渗透率。根据式(26)，固有渗透率K∞和初始固有渗透率K0∞为

式中： ψ=bmin/bmax， ψ=0.01[15]。

孔隙压力较低时，页岩渗透率主要受滑脱效应影响。而一阶滑脱系数、二阶滑脱系数与固有渗透率的大小密切相关，呈负相关的变化关系。其中一阶、二阶滑脱系数[27]可表示为

式中：B1、B2为与气体黏度、温度、弯曲度和多孔介质孔隙度等相关的系数。

由式(29)可得，C1和C10、C2和C20与固有渗透率的关系可表示为

式中：C10为初始状态下的一阶滑脱系数；C20为初始状态下的二阶滑脱系数。

将式(30)代入式(26)可得考虑微裂隙变形和滑脱效应综合作用的表观渗透率模型

2 模型验证与分析

2.1 页岩表观渗透率演化机制

1) 恒定外应力下表观渗透率演化分析。GHANIZADEH 等[28]以德国北部Toarcian Posidonia 地区的页岩样品为研究对象，其有机质成熟度为0.53%～1.45%。试验气体为CH4。对于相同页岩样品，该试验进行恒定外应力条件下的渗透率测定。2.1 节均采用GHANIZADEH 等的试验数据进行验证，模型计算结果与试验数据对比如图3 所示。

图3 恒定外应力条件下表观渗透率随孔隙压力的变化Fig.3 Variation of apparent permeability with pore pressure under constant external stress

由图3 可知，模型值与试验值具有相似的变化趋势，表明所建渗透率模型可较好地表征微裂隙表观渗透率随孔隙压力的变化关系。随孔隙压力增大，表观渗透率逐渐减小并趋于平缓。其原因为：随孔隙压力增大，气体吸附量增大，导致基质吸附膨胀变形量增大，微裂隙宽度减小，气体流动阻力增大，进而气体流量减小，使表观渗透率逐渐减小。

表1 为对应试验条件下的模型引用参数。为验证新建微裂隙变形和滑脱效应综合作用表观渗透率模型的可靠性，将表1 中模型引用参数与GHANIZADEH 等[28]的试验数据代入式(31)，可得恒定外应力条件下页岩微裂隙表观渗透率随孔隙压力的变化规律，如图3 所示。此外，通过数据拟合可得初始一阶、二阶滑脱系数、内膨胀系数和微裂隙压缩性系数，见表2。

表1 恒定外应力条件下的模型引用参数Table 1 Model reference parameters under constant external stress

表2 恒定外应力条件下的模型拟合参数Table 2 Model fitting parameters under constant external stress conditions

2) 恒定有效应力下表观渗透率演化分析。LI 等[33]采用脉冲衰减渗透仪测定Barnett 地区页岩样品渗透率，条件为恒定有效应力，改变孔隙压力。选取来自Barnett 地层有效应力为3.4 MPa 和6.8 MPa 的2 组试验数据进行比较。模型计算结果与试验数据对比如图4 所示。

图4 恒定有效应力条件下表观渗透率随孔隙压力的变化Fig.4 Variation of apparent permeability with pore pressure under constant effective stress

表3 为对应试验条件下的模型引用参数。因此将表3 中模型引用参数与LI 等[33]的试验数据代入式(31)，可得恒定有效应力条件下页岩微裂隙表观渗透率随孔隙压力的变化规律，如图4 所示。此外，通过数据拟合可得初始一阶、二阶滑脱系数和内膨胀系数，见表4。

表3 恒定有效应力条件下的模型拟合参数Table 3 Model fitting parameters under ponstant pffective ptress ponditions

表4 恒定有效应力条件下的模型引用参数Table 4 Model reference parameters under constant effective stress

由图4 可知：模型值与试验值具有相似的变化趋势，表明所建渗透率模型可较好地表征恒定有效应力条件下微裂隙表观渗透率随孔隙压力的变化关系。其原因为：有效应力恒定，随孔隙压力增大，气体分子吸附于微裂隙壁面，使得气体流动通道减小，流动阻力增大，表观渗透率减小。

2.2 气体吸附与有效应力综合作用下的滑脱系数动态演变

致密多孔介质中，气体分子与微裂隙壁面频繁碰撞，产生大于初始流动速度的“滑动”现象，称为滑脱效应。在气体渗流过程中滑脱效应产生的附加流量使页岩固有渗透率与表观渗透率间存在一定偏差。滑脱系数(C1、C2、B)是决定滑脱效应强弱的关键参数。B为整体滑脱系数(也称Klinkenberg 系数、滑脱系数)，其与一阶、二阶滑脱系数的关系[27]为

1) 恒定外应力条件下滑脱系数的演变。为进一步分析气体吸附及有效应力影响下滑脱系数的演化机制，将GHANIZADEH 等[28]的试验数据及表1 中模型引用参数和表2 中拟合参数代入式(30)、式(32)，计算得出恒定外应力条件下一阶、二阶与整体滑脱系数，如图5 所示。

图5 恒定外应力条件下滑脱系数随孔隙压力的变化Fig.5 Variation of slip coefficient with pore pressure under constant external stress

由图5 可知，当外应力恒定时，随孔隙压力增大，C1、C2逐渐减小，而B逐渐增大，且C1大于C2为主导参数。这与CHEN 等[27]研究结果一致。其原因为：随孔隙压力增大，气体分子平均自由程减小，气体分子与微裂隙壁面碰撞次数减小，滑脱效应减弱导致滑脱系数减小，因此C1、C2减小。而由式(32)可知，而整体滑脱系数B与孔隙压力成正相关，随孔隙压力的增大而增大[28]。

2) 恒定有效应力条件下滑脱系数的演变。将LI 等[33]的试验数据及表3 中模型引用参数和表4中拟合参数代入式(30)、式(32)，计算得出恒定有效应力条件下一阶、二阶及整体滑脱系数，如图6所示。

图6 恒定有效应力条件下滑脱系数随孔隙压力的变化Fig.6 Variation of slip coefficient with pore pressure under constant effective stress

由图6 可知，当有效应力恒定时，随孔隙压力增大，C1、C2和B逐渐增大，这与ZHOU 等[18]研究结果一致。其原因为：随孔隙压力增大，气体吸附变形量增大而诱导基质吸附膨胀变形量增大，进而导致微裂隙渗流通道变窄。因此，气体分子与微裂隙壁面的碰撞频率增大，滑脱效应增大，C1、C2和B逐渐增大。另外，由于其他因素影响导致C1、C2随孔隙压力增加其增加速率不太一致。

3 结果与讨论

3.1 参数敏感性分析

滑脱系数是影响页岩渗透率的主要因素之一。因此，为进一步探讨各因素影响下滑脱系数的演化规律，3.1 节选取GHANIZADEH 等[28]的试验数据进行详细讨论不同参数(微裂隙压缩性系数、分形维数和内膨胀系数)对滑脱系数的影响变化规律。

将表1 中模型引用参数、表2 中σm=16 MPa 条件下的模型拟合参数与GHANIZADEH 等[28]恒定外应力条件下σm=16 MPa 的一组试验数据代入式(30)、式(32)，可得微裂隙压缩性系数、分形维数和内膨胀系数对气体滑脱效应的影响。

1) 微裂隙压缩性系数对气体滑脱效应的影响。图7 为不同微裂隙压缩性系数下一阶和二阶滑脱系数随孔隙压力的变化规律。

图7 微裂隙压缩性系数对滑脱系数的影响Fig.7 Influence of microcrack compressibility coefficient on slipping coefficient

由图7 可知：当孔隙压力恒定时，随微裂隙压缩性系数的增大，C1、C2、B减小。其原因为：微裂隙压缩性系数表示在应力作用下微裂隙压缩的难易程度。因此微裂隙压缩性系数越大，表明页岩微裂隙可压缩空间越大，固有渗透率越大。结合式(30)可知，滑脱系数与固有渗透率成反比，即滑脱效应减弱，C1、C2、B减小。

2) 分形维数对气体滑脱效应的影响。图8 为不同分形维数下滑脱系数随孔隙压力的变化规律。

图8 微裂隙分形维数对滑脱系数的影响Fig.8 Influence of fractal dimension of micro-cracks on slipping coefficient

由图8 可知，相同孔隙压力条件下，C1、C2、B均随分形维数的增大而减小。原因为：分形维数表征微裂隙发育程度，孔隙率随分形维数的增大而增大，固有渗透率随之增大。而滑脱系数与固有渗透率成反比，因此一阶、二阶滑脱系数C1、C2、B减小。另外，相同分形维数条件下，C1、C2随孔隙压力的增大而减小。原因为：随孔隙压力增大，气体分子平均自由程减小，分子与壁面碰撞频率减小，分子邻近壁面的滑动速度逐渐减小，因此滑脱效应减弱，C1、C2减小。

3) 内膨胀系数对气体滑脱效应的影响。图9 为不同内膨胀系数下滑脱系数随孔隙压力的变化规律。

由图9a 可知：相同孔隙压力条件下，随内膨胀系数增大，C1、C2、B逐渐增大。其原因为：内膨胀系数表示基质变形对微裂隙变形的贡献。内膨胀系数越大，表明微裂隙收缩变形量越大，相应的固有渗透率越小，因此C1、C2、B逐渐增大。

由图9b 可知：相同内膨胀系数条件下，随孔隙压力增大，B呈先增大后减小的变化趋势。其原因为：当P＜0.9 MPa 时，随孔隙压力增大，基质吸附膨胀变形增大，导致微裂隙宽度减小。同时，气体分子平均自由程减小，但其减小速率小于微裂隙宽度的减小速率，因此，气体分子与壁面的碰撞频率呈增大趋势，气体滑脱效应增强，整体滑脱效应B相应增大。当P>0.9 MPa 时，气体分子平均自由程减小速率大于微裂隙宽度减小速率，导致滑脱效应相对减弱，因此B减小。

3.2 页岩表观渗透率与滑脱效应的响应机制分析

当孔隙压力较低时，滑脱效应显著，对页岩表观渗透率起积极作用。而页岩气开采过程中，滑脱效应可增大页岩表观渗透率，进一步提高页岩气采收率。因此，探究滑脱效应影响下页岩表观渗透率的动态演化机制对页岩气开采具有一定实际意义。

将GHANIZADEH 等[28]的试验数据和表1 中模型引用参数，以及表2 中拟合参数代入式(31)、式(32)，可得表观渗透率与滑脱系数、孔隙压力间的变化关系，如图10 所示。

图10 表观渗透率与整体滑脱系数、孔隙压力间的变化关系Fig.10 Relationship between apparent permeability, slip coefficient and pore pressure

由图10 可知：恒定外应力条件下，随孔隙压力增大，页岩表观渗透率均逐渐减小并趋于平缓。其原因为：孔隙压力增大前期，以滑脱效应为主，随孔隙压力逐渐增大，滑脱效应减弱，页岩表观渗透率迅速减小。而孔隙压力增大后期，气体吸附量逐渐趋于饱和，对基质变形影响可忽略，表观渗透率减小速率趋于平缓。

4 结 论

1) 不同边界条件下(恒定外应力、恒定有效应力)，随孔隙压力逐渐增大到一定程度时，在气体吸附与滑脱效应的竞争过程中，以气体吸附为主，气体吸附量增大而诱导基质膨胀变形量增大，微裂隙渗流通道减小，导致页岩表观渗透率呈减小的变化趋势。

2) 当外应力恒定时，随孔隙压力逐渐增大，气体分子平均自由程逐渐减小，分子与微裂隙壁面碰撞次数减小，滑脱效应减弱，C1、C2减小。当有效应力恒定，随孔隙压力增大，吸附变形量增大而导致基质膨胀变形量增大，固有渗透率减小，因此C1、C2、B增大。

3) 基于所建页岩表观渗透率模型进行敏感性分析，探讨分形维数、微裂隙压缩性系数、内膨胀系数与滑脱系数之间的关系。微裂隙压缩性系数、分形维数增大导致固有渗透率增大，C1、C2、B减小。内膨胀系数的增大表明微裂隙收缩变形量越大，其微裂隙宽度减小，分子与微裂隙壁面碰撞次数增大，因此C1、C2、B逐渐增大。