0.   引言

煤层气等非常规天然气是一种可持续开采的清洁能源，在国内天然气占比已经超过30%^[1]，成为中国能源系统中无法替代的一部分。在煤层气开采过程中，由于煤基质系统的渗透性能远远小于裂隙系统，导致基质和裂隙中孔隙压力将长期处于非平衡状态，因此，研究此状态下煤层的气体流动对于深入认识渗透率的演化以及与此密切相关的煤层气开采具有重要意义。

围绕气体在煤中的流动，国内外学者通常采用应力控制（恒定围压、恒定有效应力和恒定孔隙压力）和位移控制（单轴应变条件和常体积条件）两种边界来研究煤渗透率的演化规律^[2-3]。应力控制试验一般用于室内试验研究，而位移控制试验通常认为与现场边界条件一致。在围压恒定条件下，通过降低孔隙压力不仅可以研究试样在不同有效应力（围压和孔隙压力之差）下的渗透率演变，还能研究不同孔隙压力下气体吸附对岩芯渗透率的影响。有效应力恒定则可以通过保持围压与孔隙压力增量一致来实现，这种试验条件可以消除裂隙压缩效应的影响，从而简化对岩芯渗透率影响的因素。当保持孔隙压力恒定时，试验中能够消除气体吸附效应的影响，使得煤的渗透率只受力学作用的影响，同时这些试验也可以用于计算岩芯的初始孔隙度和裂隙的可压缩常数等物理参数。对于位移控制的试验研究，单轴应变条件和常体积条件均能在不同程度上模拟现场应力条件。由于在室内试验难以实现这两种条件，目前只开展了少数的几个相关试验研究。

理论分析是从机理上分析气体在煤中流动规律的研究方法，一般通过构建渗透率模型来预测渗透率的演化规律。根据Liu等^[4]将这些模型归为两类：在单轴应变条件下建立的渗透率模型和在变应力条件下建立的渗透率模型。这些模型主要研究裂隙渗透率的演化，认为裂隙开度的变化由基质吸附作用和裂隙在有效应力下的变形构成。然而这种假设却忽略了基质变形对渗透率的影响，进而低估了裂隙的渗流能力^[5]。Liu等^[6]引入内部膨胀应力的概念解释基质–裂隙之间的相互作用，较好地拟合了圣胡安盆地的现场数据和Robertson的试验数据。李小春等^[7]建立了考虑微裂缝开度分布的表观渗透率模型，对气体渗流的影响因素进行了分析。Zhou等^[8]耦合Klinkenberg效应和裂隙开度提出了一个新的渗透率模型，并对恒定有效应力下的试验数据进行了较好的解释。为了更全面地覆盖现有模型，Liu等^[9]则提出了一个适合各种边界条件的渗透率模型，并考虑了基质–裂隙之间的压差来减少评估渗透率时产生的误差。

虽然上述试验和模型对于储层渗透率演化有较好的认识，但其侧重点主要考虑平衡状态下边界条件和应力状态对煤渗透率的影响。然而，由于基质和裂隙的导流能力不同，导致裂隙中气体压力的衰减速度远远快于在基质内的衰减速度，致使煤层气开采时基质和裂隙中的气体压力以及由此造成的变形将不再同步^[10]，从而造成储层长期处于非平衡状态。由于基质的渗透率远小于裂隙的渗透率，并且基质中气体压力难以测量，因此在室内试验和模型中往往忽略基质对气体运移性能的影响^[11-12]。然而，作为煤层气的主要载体，基质的气体运移性能对煤层气长期的生产起着至关重要的影响。因此需要分别研究基质和裂隙在开采过程中气体压力和渗透率在非平衡状态下的演化规律。

对此，本文将基质–裂隙视为双重孔隙介质，耦合固体变形、瓦斯解吸和气体扩散，提出了一个理论模型研究储层开采过程中基质–裂隙孔隙压力和相应的渗透率变化。然后采用现场数据对模型进行验证，并用多场耦合的有限元方法分析了非平衡状态下基质–裂隙的孔隙压力和岩芯渗透率随时空的演化规律。本研究为正确认识开采过程中储层的压力和渗透率变化提供了理论依据。

1.   孔隙度和渗透率模型

首先定义一组控制煤基质–裂隙变形和描述流体在煤层中流动的场方程，利用本文提出的煤基质–裂隙的渗透率和孔隙度模型对这些场方程进行耦合。这些推导基于以下4点假设：①煤是包含基质和裂隙的双重孔隙弹性连续介质，每种介质都是各向同性和均质的；②煤层中仅含饱和的单相气体，且气体在基质中的扩散符合菲克定律，在裂隙中的流动符合达西定律；③储层是等温的，气体在基质中扩散和裂隙中流动均不影响温度变化；④应变远小于长度尺度。

1.1   裂隙系统孔隙度和渗透率模型的建立

由立方定律，裂隙的孔隙度${\phi _{\text{f}}}$可定义为

式中，b为裂隙开度，s为基质宽度。在开采过程中，裂隙开度和基质宽度会随着气体压力衰减发生变形，从而引起孔隙度的变化，这一变化过程可通过对裂隙的孔隙度求导表示：

同理，裂隙的渗透率k_f也可由立方定律表示^[9]：

在气体衰减过程中，通过对裂隙渗透率进行求导可得到变化的渗透率为

由于煤层气储层异于常规储层的物理特性，煤层的变形主要受到有效应力和气体吸附/解吸产生的影响。在有效应力作用下，基质和裂隙产生的变形可由胡克定律计算：

式中下标e表示有效应力；下标m和f分别代表基质和裂隙；K为煤体的体积模量，K=E/3(1-2$\nu$)，$\nu$为泊松比，E为煤体的弹性模量；$\sigma$为外部应力；p为孔隙压力；$\alpha$为基质的Biot系数，$\alpha$=1-K_m/K_s，K_s为基质骨架的体积模量；$\beta$为裂隙的Biot系数，$\beta$=1- K/K_m。

假设吸附只发生在煤基质中，且气体的吸附/解吸符合Langmuir等温原理，则煤层由于吸附/解吸产生的基质块变形$\Delta {s_{\text{s}}}$为

式中，${\varepsilon _{{\text{ms}}}}$为煤基质吸附引起的应变，${\varepsilon _{\text{L}}}$为Langmuir体积应变常数，${p_\mathrm{L}}$为Langmuir压力常数。基质块的膨胀和收缩直接影响到裂隙开度的变化，通常认为吸附产生的应变完全等于裂隙开度的变形。则由上述分析可知，基质和裂隙在有效应力和气体吸附/解吸作用下产生的变形为

将式（8），（9）分别代入式（2），（4），可得裂隙孔隙度和渗透率的导数形式：

将式（10），（11）做进一步处理可得变应力作用下裂隙的孔隙度和渗透率模型：

式中，下标0表示初始状态。单轴应变条件下水平方向的有效应力增量为^[13]

由于基质的渗透率远小于裂隙的渗透率，因此流经基质和裂隙中的孔隙压力会存在一定的区别，对于双重孔隙度竖直方向的有效应力增量为^[14]

式中，${\sigma _{zz}}$为竖直方向的外部应力。根据式（14），（15），可得煤层气储层的平均有效应力为

将式（16）代入式（12），（13），可得单轴应变条件下裂隙的孔隙度和渗透率模型为

式（17）和（18）表明裂隙孔隙度和渗透率可以作为不同孔隙压力、煤的性质、气体吸附/解吸的函数来进行评价分析。

1.2   基质系统渗透率和孔隙度的建立

基于之前的研究，基质的孔隙度与基质的有效体积应变有关^[15]：

式中，Δε_me为基质的有效体积应变，定义为基质全局应变Δε_mg和基质块内孔隙局部应变Δε_pl之差，即

而基质的全局应变与煤的全局应变$\Delta {\varepsilon _{\text{v}}}$和基质裂隙间孔隙压力差引起的基质压缩应变有关，其关系可表示为

其中，基质块内孔隙局部应变受到基质孔隙压力对骨架的压缩作用和基质内孔隙表面吸附应变影响，即

将式（20）~（22）代入式（19），可得基质的孔隙度表达式：

根据立方定律，基质渗透率为

上述分析表明：基质的孔隙度和渗透率可以表示为煤的体应变、基质的吸附应变、基质的物理性质和孔隙内气体压力的函数进行分析研究。

2.   多物理场耦合方程

2.1   煤体变形控制方程

煤体变形的控制方程可表示为^[5]

式中，G为剪切模量，G=E/2(1+ν)，${u_i}$为位移分量，F为体积力。式（25）表明煤体变形受到重力（第一、二项）、基质、裂隙的孔隙压力（第三、四项）、吸附/解吸作用（第五项）的影响，主要与基质–裂隙的Biot系数和煤体的体积模量有关。吸附/解吸作用可用Langmuir等温原理求解，基质、裂隙压力可用随后描述的气体流动方程处理。

2.2   气体流动控制方程

通常认为煤层裂隙系统中只包含游离态的气体，而基质系统中同时存在吸附态和游离态气体，应用质量守恒定律可得煤层气储层中基质和裂隙的气体流动方程为^[16]

式中，m_m和m_f分别为基质和裂隙中的气体含量：

式中ρ_mg为基质中的气体密度；ρ_ga为标准大气压力下的气体密度；ρ_s为煤体密度，ρ_fg为裂隙中的气体密度。根据气体的理想定律可以计算出不同介质中的气体密度：

忽略重力对煤层中气体流动的影响，则达西速度可以定义为

式中，$\overrightarrow {{v_{{\text{mg}}}}}$和$\overrightarrow {{v_{{\text{fg}}}}}$分别为基质和裂隙中气体的达西速度。

将式（27），（29），（30）代入式（26），可得基质中气体流动方程为

同理，将式（28）~（30）代入式（26），可得裂隙中气体流动方程为

式（30）~（33）中，V_L为Langmuir体积常数；${p_{\text{a}}}$为标准大气压力（101.3 kPa），ω为式（33）所定义的气体质量交换系数，μ为气体的动力黏度系数。

式（25），（31），（32）为双重孔隙介质中耦合的固体变形和气体流动控制方程。

2.3   边界和初始条件

对于纳维叶类型方程，边界上位移和应力条件为

式中，$\widetilde {{u_i}}(t)$和$\widetilde {{F_i}}(t)$分别为在边界$\partial \mathit{\Omega}$上已知的位移和应力。${n_j}$为垂直于边界的余弦方向向量。在域Ω内的初始位移和应力可由下式描述：

式中，u₀和${\sigma _0}$为域Ω内初始的位移和应力。式（31）所示的基质气体流动方程，狄利克雷和诺伊曼边界条件定义为

同理，式（32）所示的裂隙气体流动方程，其狄利克雷和诺伊曼边界条件为

式中，$\widetilde p(t)$和${\widetilde Q_{\text{s}}}(t)$为边界$\partial \mathit{\Omega}$上规定的气体压力和气体流量。气体流动的初始条件为

3.   模型验证

自从Gray^[17]将煤基质的膨胀和收缩考虑进煤储层渗透率模型后，关于储层渗透率模型又继续发展了30余年，其中应用范围最广的模型为Shi-Durucan模型^[13]、Palmer-Mansoori模型^[18]和Cui-Bustin模型^[19]。

为了验证本文所建模型的正确性，选取了圣胡安盆地的部分测井数据进行验证比较。根据式（18）和表 1中的基本参数，通过调整裂隙的体积模量和基质骨架的体积模量对表 2的现场数据进行拟合。对于其他储层渗透率模型，杨氏模量和泊松比等基本物理参数直接从Shi等^[20]得到，并且将Shi-Durucan模型和Cui-Bustin模型中的裂隙压缩系数、Palmer-Mansoori模型中的裂隙孔隙度视为变量来拟合现场数据。由于无法在现场分别测试出基质和裂隙的气体压力，因此对模型做进一步简化，认为基质和裂隙的气体压力相同^{[3, 21]}。新模型和几种经典模型的拟合结果见图 1。

表  1  拟合圣胡安盆地煤层所用的参数

Table  1.  Parameters used to fit coal seams in San Juan Basin

符号	值	含义	单位
E	2902[14]	煤的杨氏模量	MPa
Em	8143[20]	基质的杨氏模量	MPa
$\nu$	0.35[14]	煤的泊松比	—
pL	4.3[14]	朗缪尔压力	MPa
εL	0.01266[14]	朗缪尔体积应变	—
b0	1×10-5[20]	裂隙开度	m
s0	0.01[20]	基质宽度	m

下载: 导出CSV 

| 显示表格

表  2  圣胡安盆地Fruitland煤层测量的现场数据^[23]

Table  2.  Field data measured in Fruitland coal seam in San Juan Basin ^[23]

井号	k0/(10-15·m2)	储层压力/MPa				渗透率/(10-15·m2)
		测点1	测点2	测点3		测点1	测点2	测点3
井1	1.8	4.85	3.12	2.43		7.2	16.5	20.6
井2	6.2	4.61	3.52	3.04		1.7	2.7	2.9
井3	5.3	4.52	3.03	2.07		4.0	18.5	28.3
井4	2.1	3.38	2.88	2.21		9.3	9.6	16.2
井5	9.1	3.20	2.48	2.05		10.6	14.5	23.5
井6	2.1	2.98	2.11	1.96		11.5	24.3	28.8

下载: 导出CSV 

| 显示表格

图  1  现场数据与模型的拟合图

Figure  1.  Fitting diagrams of field data with model

下载: 全尺寸图片 幻灯片

由图 1所示的拟合结果可知，4种模型均能不同程度拟合现场数据，然而本文所建立的新模型对现场数据的匹配程度最好。这是由于新模型明确考虑了基质在有效应力下产生的变形，在渗透率模型中包含了基质变形和裂隙变形，从而建立可以反映煤层气开采过程中动态变化的渗透率模型。因此，如果未能考虑上述因素，这些模型将无法准确地预测储层渗透率。

4.   数值模拟与分析

4.1   模型的描述

第3节通过现场数据验证了模型的有效性，本节将采用数值模拟的方法对模型做进一步分析，具体为：采用COMSOL Multiphysics对岩芯内的孔隙压力和渗透率的演化规律进行研究，在单轴应变条件下充分耦合煤的变形、气体流动和瓦斯开采过程中引起的气压衰减过程。采用的模型几何宽度为50 mm，高为100 mm，该模型是基于试验室测量岩芯渗透率而来，如图 2（a）所示，图中C1，C2，C3为岩芯中轴线的四等分点。图 2（b）为数值模拟的示意图，对于位移约束条件，左、右、下端均设置为辊轮支撑，上端不设置约束条件。对于气体流动条件，除了上端与大气压力接触外，其余3个边界均为无流动边界。基质和裂隙的初始压力均设置为1 MPa，模拟所用的参数如表 3所示。所用参数都是直接从已有文献选取，计算总时长为1000000 s。

图  2  试验原理和模拟计算图

Figure  2.  Diagram schematic of experimental and numeral calculation

下载: 全尺寸图片 幻灯片

表  3  输入数值模型的参数值

Table  3.  Input parameters in numerical model

符号	取值和来源	物理意义	单位
E	2902[20]	煤的弹性模量	MPa
Em	8143[24]	基质的弹性模量	MPa
ν	0.35[20]	泊松比	—
ρs	1250[15]	煤的密度	kg/m3
φf0	2[15]	初始裂隙孔隙度	%
φm0	2[15]	初始裂隙孔隙度	%
kf0	1×10-17	初始裂隙渗透率	m2
km0	1×10-19	初始基质渗透率	m2
b	0.01[24]	裂隙开度	m
s	1×10-5[24]	基质宽度	m
μ	1.84×10-5[15]	气体黏度	Pa·s
${p_{\text{L}}}$	4.3[20]	Langmuir压力常数	MPa
α	2/3[15]	基质的Biot系数	—
εL	0.01266[20]	Langmuir体积应变常数	—
p0	0.103	标准大气压	MPa

下载: 导出CSV 

| 显示表格

4.2   岩芯气体压力的时空分布规律

岩芯中裂隙气体压力在不同时间处的分布如图 3所示。在气体压力衰减的初始阶段（t < 3000 s），随着时间增长，受扰动的岩芯范围越来越大，直至蔓延整个岩芯。当压力衰减继续进行（20000 < t < 50000 s），气体压力沿岩芯长度不再像初始阶段的急剧下降，而是绝大部分呈现出平滑的下降趋势。当气体压力衰减到后期（t > 100000 s），沿整个岩芯，气体压力减小的极为缓慢，最终达到大气压力（t=1000000 s）并保持不变。

图  3  裂隙气体压力在不同时间处的分布

Figure  3.  Distribution of fracture gas pressure at different time

下载: 全尺寸图片 幻灯片

图 4为岩芯解吸过程中不同位置处裂隙气体压力随时间的分布。解吸时间接近100 s时，靠近出口端的C1处的气体压力首先受到扰动，随着时间的增加，扰动范围逐渐扩展到岩芯中心处的C2点和靠近入口端的C3点。当解吸时间达到200000 s时解吸基本完成，3点处的气体压力均降低至大气压。

图  4  不同位置处裂隙气体压力随时间的分布

Figure  4.  Distribution of fracture gas pressure with time at different locations

下载: 全尺寸图片 幻灯片

基质中气体压力在不同时间的分布见图 5。与图 3相比，可以明显看出基质的压力变化滞后于裂隙的压力变化，这是由于基质的导流能力远低于裂隙的导流能力，与之前的研究结果一致^[22]。在解吸的初始阶段（t < 3000 s），基质中受压力扰动的范围更小，压力变化的幅度更低。随着解吸时间增长，基质孔隙压力和裂隙孔隙压力的变化趋势趋于相同，并在解吸时间达到1000000 s时降至大气压力。对于气体流动的出口端，与裂隙孔隙压力的时空分布不同，基质孔隙压力在解吸初始阶段并没有完全降至大气压，而是保持一定的压力，随着解吸的进行才逐渐降低至大气压。

图  5  基质气体压力在不同时间处的分布

Figure  5.  Distribution of matrix gas pressure at different time

下载: 全尺寸图片 幻灯片

图 6为岩芯中不同位置处的基质气体压力随时间的分布图。与图 4所示裂隙气体压力分布类似，随着解吸时间增长，靠近出口端的C1点气体压力首先开始下降，受扰动范围逐渐扩展到远离出口端的C3点。不同之处在于，C1处基质气体压力下降的时间明显迟滞于裂隙气体压力，这种迟滞性在C2处和C3处表现不太明显，体现了解吸过程中气体流动状态的不平衡。

图  6  不同位置处基质气体压力随时间的分布

Figure  6.  Distribution of matrix gas pressure with time at different locations

下载: 全尺寸图片 幻灯片

4.3   岩芯渗透率的时空分布规律

气体压力的变化将影响煤层中气体流动并导致煤体产生相应的变形，进而影响到煤层的渗透性能。图 7为非平衡气体压力下裂隙渗透率在不同时间处的分布。

图  7  裂隙渗透率在不同时间处的分布

Figure  7.  Distribution of fracture permeability at different time

下载: 全尺寸图片 幻灯片

由图 7可知，在压力衰减的初始阶段，靠近气体出口端的岩芯渗透率急剧增加，靠近气体进口端岩芯的渗透率几乎没有发生变化。随着解吸的进行，靠近气体进口端的岩芯渗透率也缓慢增加，并在解吸完全发生后（t=1000000 s），整个岩芯的渗透率保持相同并稳定不变。这是由于在初始阶段裂隙压力在靠近出口端急剧减小，此时基质中的压力还来不及变化，为了保持系统的平衡，基质内气体将向外扩散从而增加裂隙的渗透率。随着基质解吸的持续进行，基质和裂隙的孔隙压差逐渐减小，裂隙内渗透率增加的速度也平缓下来。到解吸完全发生，基质–裂隙内气体压力达到平衡，与之相关的渗透率也最终达到稳定。

图 8为非平衡状态下基质渗透率在不同时间的分布。与图 7相比，基质和裂隙的渗透率变化趋势基本相同，然而基质渗透率的变化幅度却明显低于裂隙。这是由于裂隙的应力敏感性高于基质的应力敏感性，所以在相同的孔隙压力变化下，裂隙的渗透率变化将更加明显。

图  8  基质渗透率在不同时间处的分布

Figure  8.  Distribution of matrix permeability at different time

下载: 全尺寸图片 幻灯片

图 9，10分别为岩芯不同位置处裂隙渗透率和基质渗透率随时间的分布图。随着解吸的发生，靠近气体出口端的C1点渗透率最先开始增加，而靠近出口端的C3点则最后受到扰动。与裂隙渗透率相比，相同位置处基质渗透率发生变化的时间均较晚。此外，当解吸完成后，基质和裂隙渗透率均增加到一定值保持不变，但是基质渗透率的增加幅度远小于裂隙渗透率的增加幅度。这是因为解吸造成的基质收缩虽然也会在一定程度上使得基质孔隙变大，但是对裂隙开度的扩宽作用则更加显著。

图  9  不同位置处基质渗透率随时间的分布

Figure  9.  Distribution of matrix permeability with time at different locations

下载: 全尺寸图片 幻灯片

图  10  不同位置处裂隙渗透率随时间的分布

Figure  10.  Distribution of fracture permeability with time at different locations

下载: 全尺寸图片 幻灯片

5.   结论

本文将煤层气储层视为双重孔隙介质，耦合基质变形，解吸变形和气体流动提出了一种预测变应力状态下的储层渗透率模型。通过现场数据与传统的渗透率模型进行比较，数值分析了非平衡状态下基质–裂隙气体压力和渗透率的时空分布，得到以下4点结论。

（1）在气体解吸的初始阶段，岩芯中裂隙气体压力受扰动的范围大于基质气体压力受扰动的范围；随着解吸的进行，基质气体压力和裂隙气体压力变化将趋于相同，并在解吸完成后两者压力相等并保持不变。

（2）解吸过程中煤基质和裂隙的气体压力沿岩芯长度呈现非线性变化，靠近气体出口端处，裂隙气体压力迅速减小至大气压并蔓延至整个岩芯，而基质气体压力则先略微降低至某一数值，随后向岩芯蔓延，并随着解吸的进行最终趋于大气压。

（3）解吸过程中煤基质和裂隙渗透率的变化趋势大体相似。即在解吸的初始阶段，靠近气体出口端处基质和裂隙的渗透率均急剧增加，随着解吸进行，渗透率增加的趋势向整个岩芯蔓延，并在解吸完成后达到稳定。

（4）岩芯中基质–裂隙的渗透率在解吸过程中沿岩芯长度也呈现非线性分布，且裂隙渗透率受气体压力的影响高于基质，表现为裂隙渗透率的变化幅度高于基质渗透率的变化幅度。