0.   引　　言

随着近年来国家能源结构的优化调整，煤层气作为一种清洁高效的非常规资源，其勘探与开发研究已得到高度重视^[1-2]。据国家统计局数据显示，2019年我国煤层气查明资源储量为3 040.7亿m³，截至2021年12月，煤层气累计产量已达到104.7亿m³。目前，国内已建成沁水盆地南部、鄂尔多斯盆地东缘2个煤层气商业性开发基地^[3]。但与国外煤层气储层相比，我国煤层气储层低渗特征显著^[4]，阻碍了煤层气的高效开发。

在煤层气开发过程中，孟艳军等^[5]将煤层气排采过程分为排水阶段、不稳定排气阶段、稳定排气阶段和衰减阶段。排水阶段结束后煤层气储层压力会降低至临界解吸压力以下，从而使甲烷发生解吸，水和甲烷发生气水两相渗流^[6]。初次开采阶段或提高采收率阶段时，煤储层中都存在多相流动现象^[7]。其中流体流动行为的比例以气水相对渗透率为主，而非绝对渗透率。相对渗透率被定义为多孔介质中有效渗透率与绝对渗透率之比^[7]，是影响产能的重要参数之一^[8-9]。单一流体在储层中的相对渗透率越大，代表其在储层中的阻力越小，流动性越强。尽管已有学者^[10-11]针对储层压力和流体的变化提出了考虑基质收缩、有效应力效应和气体滑脱的绝对渗透率动态模型，但在煤层气的真实开发过程中，绝对渗透率并不能充分表征储层流体的真实流动状态，需要借助相对渗透率才能够判断煤层气藏是否具有工业价值。

气水相对渗透率曲线是现场历史模拟研究、拟合以及气水产量预测的关键，但煤岩特殊的双孔隙结构与低孔低渗特性导致两相流体在储层中渗流困难。渗透率的试验精度问题也限制了相对渗透率的测定。煤层气储层不同于常规储层的排采模式，其压力传递特征决定了适用于常规储层的相对渗透率预测模型不能被直接采用。尽管20世纪80年代已经开展了煤层气储层的相对渗透率研究，但相关试验测定方法和数值模拟技术的发展较为缓慢，在一定程度上影响了煤层气商业开发的进展^[12-13]。笔者对煤层气储层相对渗透率的测定方法和数值模拟技术研究进行综述，以期为煤层气储层相对渗透率未来的研究工作提供思路。

1.   煤储层气水相对渗透率试验设计

1.1   试验方法

国内外典型地区煤层气储层渗透率分布如图1所示，与国外相比，我国煤储层渗透率普遍较低。目前稳态法和非稳态法是实验室直接获取煤储层相对渗透率常用的2种方法^[14]，其示意图如图2所示。这2种方法在单相渗透率测试中存在多方面差异（表1）。笔者将进一步介绍稳态法和非稳态法在相对渗透率测试中的应用。

图  1  典型地区煤层气储层渗透率^[4]

Figure  1.  The permeability of CBM reservoirsin typical areas ^[4]

下载: 全尺寸图片 幻灯片

1.1.1   稳态法

稳态法是将气体和水按一定流量比例同时以恒速注入岩样，当进出口压力及气水流量稳定时（即达到稳态），测定进出口压力及气水流量，再计算获得气水相对渗透率以及岩样含水饱和度。且可以通过改变注入岩心的气水流量比例获取相对渗透率曲线。稳态法已被广泛应用于相对渗透率的测试中。REZNIK等^[15]利用稳态法测定了Pittsburgh和Pocahontas煤样的相对渗透率，发现保持孔隙压力恒定的同时增加上覆压力会使含水饱和度更低。GASH^[16]以地层水作为水相流体测定了经过9天时间达到稳定状态的煤样的相对渗透率。吕祥锋等^[17]通过假设两相流体不互溶且不可压缩，探讨了稳态时围压对相对渗透率的影响。ALEXIS等^[18]在稳态条件下探究了不同有效应力下有效渗透率的变化规律，并建立了两者之间的函数关系。稳态法的优点在于可以获得范围更广的相对渗透率曲线，其测定结果可信度高且不需要单独获取岩心孔隙度。然而，稳态法也有局限性，其试验过程耗时长并要求岩样本身具有较高的渗透率。因此稳态法更适用于砂岩和碳酸盐岩储层^[12]，在煤层气储层中的应用较少。

图  2  稳态法与非稳态法测定方法示意^[7]

Figure  2.  Schematic diagrams of the steady-state method and the non-steady-state method ^[7]

下载: 全尺寸图片 幻灯片

表  1  稳态法与非稳态法试验原理对比（改自文献[9]）

Table  1.  Comparison of experimental principal between the steady-state method and the unsteady-state method (modified from [9])

参数	稳态法		非稳态法
	恒压法	恒流法
测定范围	>10−19 m2	>10−19 m2	<10−19 m2
基础理论	达西定律
核心方程	$K=\dfrac{qL}{A\mathrm{\Delta }H}$	$K=\dfrac{qL}{A\mathrm{\Delta }H}$	$K=\dfrac{\omega \mu {L{C}_{{\rm{f}}}V}_{{\rm{u}}}{V}_{{\rm{d}}} }{A\left({V}_{{\rm{u}}}+{V}_{{\rm{d}}}\right)}$
试验原理	流体在压差作用下，煤样两端的流量与压力梯度达到稳态时呈现线性关系，再利用流量测定获取渗透率	以恒定流速在煤样中注入流体，通过计算煤样两端压力差随时间的变化来获取渗透率	在仪器整体达到压力稳定后，使上游室的压力瞬间增大，由此产生可以传递至下游室的压力脉冲。在此过程中记录上下游室的压力差，从而计算得到煤样渗透率
试验时间	较长		较短
初始条件	$ H\left(z,0\right),0 < z < L $	$ H\left(z,0\right),0 < z < L $	$ H\left(z,0\right),0 < z < L $
边界条件	$ \begin{array}{c}H=0,z=0\\ H=\Delta H,z=L\end{array} $	$\begin{array}{c}H=0,z=0\\ \dfrac{\partial H}{\partial z}=\dfrac{I}{KA}\left(q-{C}_{{\rm{e}}}\dfrac{\partial H}{\partial t}\right),z=L\end{array}$	$\begin{array}{c}\dfrac{ {\mathrm{S} }_{{\rm{d}}} }{KA}\dfrac{\partial H}{\partial t}-\dfrac{\partial H}{\partial z}=0,z=0\\ \dfrac{ {\mathrm{S} }_{{\rm{u}}} }{KA}\dfrac{\partial H}{\partial t}+\dfrac{\partial H}{\partial z}=0,z=L\end{array}$
数据形态图
注：K—渗透率，mD；q—流泵产生的流量，cm3/t；L—煤样高度，cm；A—横截面积，cm2；ΔH—定水压力差，MPa；H—水压，MPa；z—与煤样端面间距，cm；Ce—流泵压缩储水系数，无量纲；Su，Sd—上、下游渗透系统压缩储水系数，无量纲；t—时间，s；ω—衰减曲线斜率；μ—流体黏度，mPa·s；Cf—流体可压缩系数，无量纲；Vu，Vd—上、下游容器体积，cm3；P—压力，MPa; Ps—稳态法初始注入压力，MPa；ΔP—压力变化值，MPa；P's—非稳态法初始两侧注入的压力，MPa；P'e—非稳态法上下游稳定压力，MPa；ΔP'—施加脉冲压力，MPa

下载: 导出CSV 

| 显示表格

稳态法忽略了毛细管压力和重力作用，针对气相和水相流体的有效渗透率可以根据达西定律进行计算，如下所示。

水相流体：

气相流体：

式中：$ {K}_{\mathrm{e}\mathrm{w}} $和$ {K}_{\mathrm{e}\mathrm{g}} $分别为水相和气相有效渗透率，mD；Q_w为水相流量，cm³/s；Q_g为气相流量, cm³/s；P₀为下游压力，MPa；P₁为上游压力，MPa；$ {\mu }_{\mathrm{g}} $为气体黏度，mPa·s。

相对渗透率计算方法如下。

水相相对渗透率：

气相相对渗透率：

其中：$ {K}_{{\rm{a}}} $为绝对渗透率，mD；$ {K}_{\mathrm{r}\mathrm{g}} $和$ {K}_{\mathrm{r}\mathrm{w}} $分别为气相和水相相对渗透率，无量纲。

1.1.2   非稳态法

非稳态法最初由RAPOPORT等^[19]提出。在非稳态法试验过程中，岩石中的一种流动相被另一种以恒定流速注入的不混溶流动相驱替^[20]。记录从上游至下游过程中脉冲压力的压差即可通过进一步计算来获取相对渗透率曲线^[21]。有学者进行煤岩非稳态法试验时先后测定了气驱水和水驱气的相对渗透率曲线^[22]，发现在不同气体驱替水后会使束缚水饱和度发生变化。此外，已有大量学者使用非稳态法进行了气驱水测试^[23-25]，证实了非稳态法可以为两相流体渗流行为的深入探究提供坚实的试验基础。非稳态法整体上测试程序相对简单，测定时间短，相较于稳态法而言对绝对渗透率要求低。但非稳态法获取相对渗透率的数据范围较窄，曲线端点往往需要推测获得。

非稳态法相对渗透率计算主要是由Johnson-Bosslex-Naimann方法（简称JBN法）计算。

式中：$ {f}_{\mathrm{w}2} $为流体中水所占的比例因子，无量纲；$ {\mu }_{\mathrm{w}} $为水的黏度，mPa·s。

可以进一步变形得到

其中：

式中：${I}_{{\rm{r}}}$为相对注入量，无量纲；$ {N}_{\mathrm{p}} $为孔隙体积中累计出水量，cm³；$ {Q}_{\mathrm{g}1} $为孔隙体积中累计气体注入量，cm³；i为累计次数。

其中：$ {f}_{\mathrm{g}2} $为流体中气体所占比例因子，无量纲。

但JBN法对储层中的气−水−煤相互作用机制欠缺考虑，会在一定程度上导致相对渗透率的计算结果存在误差。应用JBN法需要满足3个假设^[26]：① 忽略毛细管末端效应；② 两相流体的密度保持不变；③ 整个岩心存在活塞式位移。煤储层气水两相相对渗透率的测试过程很难满足后2个假设，因此不能涉及到相对渗透率曲线的全部内在因素（如气体吸附）。

总体而言，煤层气储层低渗、低孔及非均质性强的特点，导致稳态法测量相对渗透率相对困难且准确性较低，因此使用非稳态法进行相关研究逐渐占据了优势地位^[27]（表2和图3）。

表  2  稳态法与非稳态法测定渗透率部分结果统计

Table  2.  Statistics of partial results of permeability determined by the steady-state method and the unsteady-state method

方法	文献来源	测试流体		样品尺寸（直径×高度/mm×mm）	采样地区	最大相对渗透率
		气相	水相			气相	水相
非 稳 态 法	PURL等[28]	氦气	水/地层水	—	圣胡安河盆地	0.8～0.9	1
					黑勇士盆地	0.58～0.68	1
	GASH等[16]	氦气	普通水	～90×80	圣胡安河盆地	0.51～0.71	0.24～0.4
	RAHMAN等[29]	空气	普通水	—	澳大利亚	0.33～0.58	0.26～0.8
	SHEDID等[30]	甲烷/CO2	地层水	63.5×114.3	澳大利亚	0.12～0.57	0.8～0.96
	SHEN等[31]	甲烷	蒸馏水/地层水	50×100	沁水盆地	0.08～0.33	0.1～0.25
	ZHANG等[32]	氦气	普通水	48.5×100	沁水盆地	0.68～0.71	0.78～0.83
	SUN等[14]	氦气	地层水	25×60	沁水盆地	0.23	0.05
					鄂尔多斯盆地	0.68	0.2
	ZHANG等[33]	甲烷	地层水	25×50	沁水盆地	0.18～0.26	—
	SHAW等[34]	氦气	地层水	6.1×7.9	格洛斯特盆地	0.12～0.62	1
	SHEN等[35]	甲烷	普通水	直径：25	二连盆地	0.40～0.86	0.5～0.93
稳 态 法	GASH[16]	氦气	普通水	88.9×76.2	圣胡安河盆地	1	0.29
	ALEXIS等[18]	氦气/CO2	地层水	～50×100	西弗吉尼亚州马歇尔	0.02～0.11	0.03～0.9

下载: 导出CSV 

| 显示表格

图  3  试验流体使用频率及稳态法和非稳态法使用频次对比

Figure  3.  Frequency of selection of experimental fluid and comparison of the frequency of use of the steady-state method and the non-steady-state method

下载: 全尺寸图片 幻灯片

1.1.3   微流控试验

常规的室内试验测试无法定量研究流体在煤岩割理结构中的流动行为，也很难将相对渗透率的测定与孔裂隙系统的几何拓扑形状联系起来。近年来兴起的微流控系统在可视化多孔介质中多相流体的流动特性方面有巨大优势。JOSEPH等^[36]使用微流控技术测量了多孔介质中的渗透率，并研究了渗透率与有效孔隙度之间的关系。XU等^[37]发现了在均匀随机网络中的渗透率测试值与Carmen-Kozeny方程估计值保持一致。然而，这些研究存在的最大问题是没有考虑到真实储层的实际孔裂隙几何结构。MAHONEY等^[38]提出了一种使用激光刻蚀和反应离子刻蚀制备岩石内部微通道的新方法，并研究了煤润湿性的影响因素。GUNDA等^[39]通过真实样品的扫描电子显微镜（SEM）图像对储层进行了空间重构，进而将其制作成了微流控芯片。SONG等^[40]基于方解石的微流控芯片模拟了二次和三次采油的注酸过程，可视化了基质与流体的相互作用。目前针对二维的微流控试验已经有较多的研究成果，但由于驱替速度快、空间结构复杂，可视化三维裂隙体系中的流体行为较为困难，并且很少模拟多相流的情况。基于以上问题，GERAMI等^[41]利用微米CT和SEM成像试验获取并量化了裂隙网络，构建了一套可以高效进行多相微流控试验的流程（图4）。并计算了相对渗透率曲线，预测了水相最大相对渗透率为0.15。SONG等^[42]引入深度学习算法建立了相对渗透率智能模型，研究了两相相对渗透率的关键影响因素。但微流控试验在煤储层中的相对渗透率研究较少，也缺乏与地质控制模式之间的联系。

图  4  微流控试验流程示意图（改自文献[41]）

Figure  4.  Schematic diagram of the microfluidic experiment flow (revisd from [41])

下载: 全尺寸图片 幻灯片

1.2   试验流体

煤层气储层相对渗透率试验通常采用气水两相流，其中由于试验目的不同，水相相对渗透率测试所用流体也有所差异。在早期发展阶段，REZNIK等^[15]利用蒸馏水分别进行了稳态法和非稳态法试验的煤岩相对渗透率测定，但后续研究中蒸馏水利用的较少。为真实模拟储层条件，许多学者通常利用地层水进行试验^{[31, 34, 43]}。在无法获取煤层气储层原地地层水条件下，可利用KCl等溶液进行地层水盐度和PH值的配置以达到模拟地层水的效果^[33]。近几年来在水相流体的选择中，地层水或模拟地层水的利用率逐渐增加（图3），但由于试验需求不同或受试验条件限制多数研究仍使用普通水（未指定的情况下通常为自来水）作为测试流体。

相对于水相流体的选择，气体类型的选择更为多样（图3）。根据已有文献显示，REZNIK^[15]最早使用空气作为测试流体，但实际在后续的试验研究中，多数学者采用纯净气体进行测试。为避免煤基质在压力作用下发生气体吸附进而影响相对渗透率的测定，惰性气体常常被认为是控制变量的理想气体。在煤岩相对渗透率的测试研究中氦气的使用最为常见^{[23, 34, 44]}。为了探究煤层气储层中气体吸附/解吸影响下的渗透率变化机制，或模拟真实储层条件的需求下，甲烷在测试流体的使用率逐渐增多（图3）^[35]。此外，由于N₂和CO₂同样能被煤基质吸附，加之CH₄在室内试验存在一定安全风险，因此这2种气体也常作为CH₄的替代气体。

1.3   试验装置

稳态法和非稳态相对渗透率试验测试装置基本相同，主要由流体注入系统、驱替系统、高精度流量计、岩心夹持器和气水分离装置组成^[13]。由于煤岩的低渗特点，因此相对渗透率测试中气分离和流量精确测量是试验结果准确的关键。GASH等^[16]在气水分离器下安装电子秤，从而通过重量变化来估算煤样的含水饱和度（图5B1）。DURUCAN等^[44]为减少误差并准确记录流量，设计了高达1.5 m且足够细的气水分离试管（图5B2）。近几年来不同学者通过采用新兴技术手段对试验装置进行改进来实现含水饱和度和流量的准确测量，包括核磁共振技术、CT成像技术等。核磁共振技术的应用可以实时观察煤样驱替内在现象，并通过核磁共振数据计算含水饱和度解决了驱替过程中死体积带来的误差。潘一山等^[24]将核磁共振成像技术（NMRI）应用于非稳态法试验中，在三轴渗透仪周围安装NMRI成像探头，可获取煤样水驱气的实时图像特征，结果表明水驱前缘具有不一致性，并且驱替过程中存在保持不变的优势通道。WANG等^[45]没有直接将核磁共振设计到渗流仪中，而是在岩心加持器周围加入加热设备，旨在对不同温度下不同驱替阶段的煤样进行成像来观察驱替现象。发现温度的升高会使气体突破水相的速度加快，从而导致气水界面变得不规则更加复杂。SUN等^[14]通过特殊设计来实现非磁性和非金属材料的岩心夹持器与NMR分光计的连接（图5A3），研究表明T₂谱振幅与煤样中含水量有很强的线性关系；进而利用T₂谱振幅计算驱替过程中的含水饱和度并分析不同孔径在驱替过程中含水量的变化。但由于试验设备的限制，目前只能停留在尺寸为5 cm的样品。

图  5  煤层气储层相对渗透率测试的不同类型试验装置

Figure  5.  Different experimental setups for relative permeability testing

下载: 全尺寸图片 幻灯片

CT成像技术也可以被直接应用于煤岩气水两相流流动现象的观察。ALEXIS等^[18]分别进行He-地层水和CO₂-地层水条件下的稳态试验，并在岩心加持器周围布置CT扫描探头（图5A1），试验结果表明水分的存在对气体的流动有一定的阻碍作用。

2.   煤储层气水相对渗透率基本模型

1）Brook-Corey模型

Brook-Corey模型是目前多孔介质中多相流体相对渗透率计算最为广泛运用的模型，由BROOKS 和COREY^[46]在毛细管束模型的基础上（图6b）考虑孔径分布差异化提出的。其基本表达式如下：

图  6  煤层气储层理想几何模型

Figure  6.  Ideal geometric model of CBM

下载: 全尺寸图片 幻灯片

式中：$ {K'}_{\mathrm{r}\mathrm{w}} $和$ {K'}_{\mathrm{r}\mathrm{n}\mathrm{w}} $分别为润湿相和非润湿相相对渗透率，无量纲；λ为孔隙分布指数，无量纲；${S}_{{\rm{w}}}^{*}$为标准化处理的润湿相流体饱和度，可表示为

式中：$ {S}_{\mathrm{w}} $为含水饱和度，%；$ {S}_{\mathrm{w}\mathrm{r}} $和$ {S}_{\mathrm{n}\mathrm{w}\mathrm{r}} $分别为润湿相和非润湿相的残余饱和度，%。

不同学者以Brook-Corey模型为基础，结合绝对渗透率动态预测模型，来建立有效渗透率的预测新模型。XU等^[49]通过引入Palmer-Mansoori（PM）模型对有效应力和基质收缩效应进行了探究，建立了气水有效渗透率的预测模型。SUN等^[14]利用气水两相流偏微分方程，考虑生产过程中气体解吸和相对渗透率变化，利用迭代法求解推导出无限小煤层中压力和饱和度的关系，将Corey模型与PM模型相结合得到新的有效气/水渗透率的半解析模型。由于煤层气储层同时受有效应力和基质收缩正负效应的影响，在两模型耦合过程中通常需要校正端点有效渗透率。

该模型在砂岩等固结多孔介质中的适用性已被广泛验证。但由于其是基于毛细管模型推导而来，考虑到煤岩的割理系统中并不具有类似毛细管的几何形状，并且未将流体间的界面作用纳入其中，因此该模型对于煤岩渗透率计算的适用性仍存在争议。

2）多元线性回归方法

SHEN等^[31]通过试验数据发现气相和水相渗透率与含气饱和度之间关系可以用以下公式表示:

式中：S_g为气体饱和度，%；a，b，m，n为依赖于煤岩性质的系数，无量纲；通过非线性多元回归的最小二乘法拟合，相对渗透率曲线形状的变化仅受参数b，n的影响；a，m反映煤岩的绝对渗透率，且参数与煤阶和显微组分有一定的关系。通过多元回归分析可以获得一个与煤岩性质有关的综合模型：

式中： R_{o, max}为最大镜质体反射率，%；V、I、M和AC分别为镜质组、惰质组、矿物和灰分含量，%。S_ng为气相饱和度，%；

其中，S_ng可表示为：

其中， S_wo为束缚水饱和度，%；S_go为束缚气饱和度，%。

此外，ZHU等^[50]选用最大的水相和气相渗透率确定参数a、m，并用标准化的含水饱和度作为自变量。但该模型参数要求多，仅依赖有限的试验数据获取参数对结果影响较大。

3）CHEN模型^[51]

煤层气储层具有特殊的双重孔隙结构，因此基于毛细管束提出的相对渗透率模型会不适用于煤层气储层。CHEN等^[51]基于火柴棍模型（图6a），引入残余相饱和度和形状因子并考虑了气水流动过程中孔隙结构改变而带来的影响。该改进模型相对渗透率的计算依赖于饱和度和孔隙的变化，其基本表达如下：

式中：ξ为割理尺寸分布指数，无量纲；η为割理弯曲度指数，无量纲；J为用于校正因孔隙变化导致ξ变化的形状因子，无量纲；$ {K}_{\mathrm{r}\mathrm{w}}^{\mathrm{*}} $为润湿相终止点相对渗透率，无量纲；$ {K}_{\mathrm{r}\mathrm{n}\mathrm{w}}^{\mathrm{*}} $为非润湿相终止点相对渗透率，无量纲。

CHEN^[48]等利用中国、澳大利亚、美国和欧洲等地32组煤样相对渗透率进行了测试，结果表明ξ的值总体在0.3～8.8，在ξ>3时对相对渗透率曲线的影响不明显；且随着ξ的减小，气体渗透率增加，水相渗透率增大。η的值在0～2，随η值的增加气水相对渗透率下降，同时相对渗透率曲线由凸变凹。ξ随固定碳含量和最大镜质组反射率的增加成U型变化，而η与煤阶没有明显关系。

该模型是目前唯一充分考虑了煤储层自身特征建立的气水相对渗透率计算模型，因此对煤储层适用效果较好。

4）LIU模型^[8]

LIU等^[8]分别基于毛细管模型和裂隙模型被润湿相流体润湿前后的的几何参数和饱和度变化推导。得出相对渗透率与饱和度的函数关系后，在考虑气体滑脱效应对相对渗透率的影响基础上，发现了不同模型下函数形式仅有指数项存在差异，归纳出统一模型：

式中：α反映不考虑气体滑脱效应下，给定含水饱和度下气体的相对渗透率；β反映含水饱和度对气体滑脱渗透率的影响；$\sigma ={b}_{0}/{P}_{{\rm{m}}}$反映平均孔隙压力对相对渗透率的影响；b₀为气体滑脱因子，无量纲；P_m为毛细管平均孔径压力，MPa。

该模型是根据不同理想几何模型推导下的统一表达式，可以广泛应用多种非常规储层，仅需要确定不同非常规储层条件下的指数项。由于β决定了气体滑脱渗透率对气体相对渗透率的影响，随β的增大，气体的相对渗透率减小。随α的减小，气体的绝对相对渗透率增大，从而使气体的相对渗透率增大。此外，图7对比了LIU模型^[8]和SHEN模型^[31]，结果显示在多孔介质中考虑气体滑脱效应时，LIU模型的适用性更强。

图  7  SHEN模型和LIU模型预测效果对比图（改自文献[8]）

Figure  7.  Comparison of prediction effects of Shen model and Liu model (modified from [8])

下载: 全尺寸图片 幻灯片

但该模型仅限于气相渗透率的计算，水相渗透率的计算过程未能给出，无法获取完整的相对渗透率变化曲线。

5）相对渗透率分形模型

PURCELL等^[52]利用毛细管压力曲线逆推得到孔径，再带入泊肃叶定律并结合达西定律推导出相对渗透率计算公式；BURDINE等^[53]引入毛细管弯曲度对相对渗透率的计算公式进行了改进。压汞法是获取静态毛细管压力数据的方法之一，但在高压下煤基质会被压缩甚至破裂从而破坏原有的孔隙结构，导致利用分形维数计算的相对渗透率出现偏差，因此需要对高压段的压汞体积进行校正^[32]。前人研究表明储层岩石孔隙具有分形结构性质^[54-55]，贺承祖等^[56]率先结合分形理论对相对渗透率计算公式积分得到：

式中：${S}^{\mathrm{{'}}}$为有效润湿相饱和度，%；D为分形维数。

但贺承祖和华明琪的模型没有考虑不同围压条件下相对渗透率的变化特征^[56]。ZHANG等^[32]进一步对比了不同围压下相对渗透率测试数据，发现相对渗透率曲线形状相似但发生了水平方向上的位移，因此对该模型改进为：

式中：c为常数，仅取决于煤岩样品和围压。

PENG等^[57]基于Hagen-Poiseuille二阶滑动连续方程，引入权重系数和分形理论，考虑气体扩散和气体滑脱效应、孔喉结构参数、含水饱和度分布、气水界面压降等参数建立了可以反映不同温压条件的计算模型。

总体上来说，分形维数模型通过分形理论对煤岩孔隙分布特征进行了考虑，但仍基于毛细管束模型进行推导，对煤岩的适用性仍需进一步讨论。

由于煤岩独特的孔裂隙结构，目前对相对渗透率计算函数的理解还远远不够。尽管基于不同理论的相对渗透率计算模型都存在一定的优点（表3），但针对不同的适用条件，仍存在较大差异并难以实现完全统一。此外，以上数学模型大部分是将岩石物理性质集中到与饱和度有关的相对渗透率函数中。这样所带来的问题是无法探究内在孔裂隙结构与有效渗透率的联系。

表  3  不同基本计算模型优缺点统计

Table  3.  Statistics of advantages and disadvantages of different basic computing models

名称	基本表达式	几何模型	优点	缺点
Brook-Croey 模型[37]	${K{'} }_{\mathrm{r}\mathrm{w} }={({S}_{\mathrm{w} }^{\mathrm{*} })}^{\tfrac{2+3\lambda }{\lambda } }$ ${K{'} }_{\mathrm{r}\mathrm{n}\mathrm{w} }={({1-S}_{\mathrm{w} }^{\mathrm{*} })}^{2}{({1-S}_{\mathrm{w} }^{\mathrm{*} })}^{\tfrac{2+\mathrm{\lambda } }{\mathrm{\lambda } } }$	毛细管束模型	参数少，计算简单快捷	毛细管束模型对煤割理系统不完全适用，其次模型未考虑流体间界面作用
SHEN模型[31]	$ \begin{array}{c}{K}_{\mathrm{r}\mathrm{g}}=a*{S}_{\mathrm{g}}^{b}\\ {K}_{\mathrm{r}\mathrm{w}}={m}^{*}{\left(1-{S}_{\mathrm{g}}\right)}^{n}\end{array} $	—	参数通过含气饱和度和相对渗透率拟合获取，无需额外计算；考虑了煤岩自身物质组成特性	参数要求多，且仅根据数量和类别有限的岩石测试数据来计算参数会导致结果有较大差异，普适性较差
CHEN模型[51]	${K{'} }_{\mathrm{r}\mathrm{w} }={ {K}_{\mathrm{r}\mathrm{w} }^{\mathrm{*} }({S}_{\mathrm{w} }^{*})}^{\eta +1+\tfrac{2}{J*\xi } }$ ${K\mathrm{ {'} } }_{\mathrm{r}\mathrm{n}\mathrm{w} }={ {K}_{\mathrm{r}\mathrm{n}\mathrm{w} }^{\mathrm{*} }(1-{S}_{\mathrm{w} }^{\mathrm{*} })}^{1+2/\left(J*\xi \right)}$	火柴棍模型	基于火柴棍模型，更适用于煤割理结构，并考虑了孔隙结构变化	J参数为不代表物理性质的拟合参数，为多项式函数，增加里模型的复杂性却无实际物理意义
LIU模型[8]	${K}_{\mathrm{r}\mathrm{g} }={\left(1-{S}_{\mathrm{w} }\right)}^{\alpha }+\sigma {\left(1-{S}_{\mathrm{w} }\right)}^{\beta }$	毛细管束模型、裂隙模型	考虑了气体滑脱效应对相对渗透率的影响，并将孔隙和裂隙模型进行结合	该模型只进行了气相相对渗透率的计算，无法计算水相相对渗透率变化
分形维数 模型[43]	$\begin{array}{c}{\mathrm{K} }_{\mathrm{r}\mathrm{w} }={ {S}^{{'} } }^{\tfrac{11-3D}{3-D} }\\ {K}_{\mathrm{r}\mathrm{n}\mathrm{w} }={\left(1-{S}^{{'} }\right)}^{2}\left(1-{S{'} }^{\tfrac{5-D}{3-D} }\right)\end{array}$	毛细管模型	引入分形维数表征孔隙分布 特征	模型基于毛细管束模型，对煤岩的适用性仍需进一步讨论

下载: 导出CSV 

| 显示表格

3.   煤储层气水相对渗透率数值模拟

相对渗透率的试验测定通常会受控于样品的尺寸，无法避免毛细管末端效应^[34]。数值模拟技术可以利用生产数据通过数学方法进行拟合，并依据区块特征对参数进行调整以最优化求取相对渗透率曲线，从而弥补在特殊情况下实验室无法测定储层相对渗透率的缺点。

1）数值反演方法

数值反演方法最大的优势在于可以解释毛细管效应和储层的非均质性^[58]。数值反演模拟获取相对渗透率的曲线主要包括目标函数的确立、相对渗透率基本模型的确定和模型控制参数的调整3个步骤^[59]。相对渗透率模型已经在第2节中进行了探讨。模型控制参数调整主要通过优化算法来实现，目前常用的优化算法有梯度算法、随机算法和混合算法等。梯度算法通过计算目标函数的Hessian矩阵或梯度来确定目标函数搜索方向，其最大的优点在于高效和良好的收敛性，但耗时长，因此快速获取梯度是提高收敛速度的关键。需要注意的是，如果最初没有进行正确的初始化，基于梯度的方法很容易陷入局部最小值^[26]。在目标函数无法明确表达时，随机算法作为启发式算法便可作为解决这类复杂问题的途径，常用的随机算法包括遗传算法、粒子群优化算法，但其计算量较大，通常需要数百次甚至以上的模拟运行找到最优解。混合算法是几种算法的组合，相比单一算法，其计算精度和收敛速度都可以得到保证。

在数值反演计算过程中通常将端点饱和度作为已知条件，但初始条件对端点饱和度影响非常敏感，因此在计算中准确估算端点饱和度是仍需解决的问题，另外由于数值反演过程中数据规模和计算过程复杂，兼顾更高计算精度和更快收敛速度的优化算法需进一步研究。

2）历史拟合法

历史拟合法是迄今为止研究原位条件下计算相对渗透率的主要方法。历史拟合法是在地层深度、储层有效厚度、原始压力、初始含气量等储层性质参数基础上建立数值储层模型，通过调整参数使模拟结果与生产数据相拟合，直到模拟输出与历史数据控制在合理误差之内（图8）。

图  8  煤层气储层相对渗透率数值模拟流程

Figure  8.  Flowchart of numerical simulation on relative permeability in CBM reservoirs

下载: 全尺寸图片 幻灯片

在气−水−煤复杂系统中，室内测试并依据数学模型计算的相对渗透率精度不能得到保证，而历史拟合方法考虑到了气体可压缩性，吸附性和煤的渗透性等机理，实用性较好。SHAW等^[34]以修正的Brooks-Corey模型为模拟软件ECLIPSE的输入模型，初始化模拟阶段输入渗透率、孔隙度、毛细管压力等相关数据，通过调整拟合指数与非稳态法测试数据拟合。目前已有包括基于梯度方法、基于集成方法等多种算法辅助历史拟合，但历史拟合之前需假定气水两相饱和端点已知。ZHANG等^[26]基于混合贝叶斯优化和网格自适应直接搜索算法提出非稳态岩心驱替相对渗透率自动拟合生产数据的历史拟合方法，其最大优势在于无需确定端点饱和度，并且具有较高的收敛速度，模拟效果如图9b所示。历史拟合方法的准确性取决于基本输入参数，其输出的拟合结果需要检验，AMINIAN^[62]认为历史拟合是获取真实渗透率的唯一实用方法。由于试验测试的岩心柱和实际储层差异，历史拟合法获取的相对渗透率曲线往往与试验结果有较大差异。另外，历史拟合方法的准确性还依赖于输入模拟参数的准确性和全面性，因此在历史拟合中需要对输入参数的代表性和准确性进行讨论。

图  9  煤层气储层相对渗透率数值模拟结果对比

注：DRP为动态相对渗透率

Figure  9.  Comparison of numerical simulation results

下载: 全尺寸图片 幻灯片

3）生产数据分析

CONWAY等^[63]认为煤层气储层的相对渗透率可以通过试井分析进行估算。CLARKSON^[64]详细阐述了利用生产动态数据获取相对渗透率的过程：首先获取生产井的流动压力和关井压力、以及气水生产数据；利用适用于煤层气井的物质平衡方程匹配现有关井压力，再通过径向流动方程，求解各相的有效渗透率作为储层压力的函数；然后用校正过的物质平衡方程计算，计算归一化处理的含水饱和度并作为累计产水量的函数，以含水饱和度作为横坐标绘制气相和水相的有效渗透率曲线；最后估算基础渗透率以求取相对渗透率（图8b），生产数据分析的结果如图9a所示。XU^[49]等通过获取累计产气和产水曲线和利用物质平衡方程分析两相流阶段储层压力和含水饱和度变化，并用室内试验数据标定端点相对渗透率，结合Crory模型求取相对渗透率曲线，此外，由于煤储层和非均质性，试验测试很难获得代表性样品，因此生产数据分析和岩芯分析获取的相对渗透率曲线存在很大差异。CLARKSON等^[60]基于此系统研究了动态裂缝渗透率和孔隙度、多层效应和瞬变渗流效应等因素对煤层气相对渗透率曲线形态的影响。

目前历史拟合法在煤层气相对渗透率模拟中应用相对较为广泛，但总体上煤层气储层相对渗透率的数值模拟研究发展缓慢。早期源于常规储层的数值模拟技术，但由于煤层气储层特殊的双孔隙结构与吸附为主的赋存方式，不能直接应用排水降压的生产模式。

总之，煤储层相对渗透率在目前在试验开展，理论模型建立，数值模拟等方面均取得了一定的发展。不同研究方法手段的结合使得相对渗透率的分析更加准确科学（图10）。

图  10  煤层气储层相对渗透率试验及数值模拟技术研究框架

Figure  10.  Framework of relative permeability experiment and numerical simulation techniques of CBM reservoir

下载: 全尺寸图片 幻灯片

4.   结　　论

1）由于煤层气储层的超低孔、低渗特征，主要采用非稳态法对其相对渗透率进行试验测试，但其未能考虑储层介质空间的几何形态和复杂程度，微流控试验可以作为可视化和定量化拓扑结构约束下的多相流体渗流行为的工具，但很少与地质控因相联系；此外近年来测试流体中地层水和甲烷使用频率的升高表明未来研究趋向于还原真实的储层条件。

2）合理设计气水分离装置并且应用核磁共振技术、CT扫描技术等新兴试验手段进行改进可以减小渗透率等参数的计算误差，并且实时捕获两相流的驱替过程。

3）目前基于毛细管压力和理想几何模型已建立多种应用于煤层气储层的相对渗透率计算模型，但不同模型的假设条件和适用范围仍存在许多差异，需进一步验证比较，并建立考虑相饱和度润湿顺序和真实孔裂隙结构的模型。

4）数值模拟技术可以突破试验样品的尺寸限制，结合煤层气现场数据与数值反演、历史拟合、生产数据分析等方法可以动态刻画煤层气生产过程中的渗透率变化，但由于岩心尺度问题，存在室内测试结果与数值模拟不统一的情况。

5）未来煤层气储层相对渗透率发展趋势应至少包括以下3方面：① 尽管不同压力条件下的有效渗透率试验测试已经被广泛研究，但原位温度条件下的相对渗透率变化规律仍有待进一步揭示。② 相对渗透率试验测试流体多为单一气体和水相的渗透率测定，多种混合气体和水相渗透率测定尚未得到充分研究。③ 虽然微流控试验已经应用于相对渗透率曲线的测定，但由于其成本较高、操作相对复杂，目前未建立统一的测试工序流程，仍需深入拓展其应用场景。 ④数值模拟在煤层气储层相对渗透率分析中研究应用相对较少且尚未成熟，相对渗透率曲线往往与实验室获取的结果有较大差异，因此数值模拟和物理试验的联动分析与应用将是未来的研究重点。