0.   引言

高效开发非常规天然气是增强我国能源保障、落实绿色发展战略、助力“双碳”目标的重要途径^[1]。煤系气是指蕴藏于煤系中的非常规天然气，其中煤层瓦斯是主要组成部分^[2]。煤层瓦斯高效抽采是煤矿安全生产的重要前提^[3]，然而目前缺乏对煤层瓦斯运移机制的充分认识。

目前，学者们在研究煤体瓦斯微观运移特性时，多以等效孔隙网络模型为切入点，对微观孔隙结构进行渗流过程研究^[4]。刘学锋^[5]综合考虑浸润性、流体动力学黏度等因素，基于双相格子Boltzmann模型对灰岩数字岩心油水驱替过程进行模拟。崔冠哲等^[6]基于砂岩孔隙空间结构网络模型，使用格子波尔兹曼法对各孔径范围的孔隙结构渗流规律进行研究，发现微孔数量对煤岩体连通性的影响较大。王刚等^[7]基于CT三维重建模型重构了6种煤样的孔隙网络结构模型，探究非达西渗流中速度、压力、非达西系数等参数对渗流的影响。白若男^[8]模拟研究了微观尺度下煤岩体孔隙结构单相、多向低压水渗流过程。虽然诸多学者探讨了煤微观孔隙瓦斯运移特性^[9]，但忽略了瓦斯解吸−扩散过程。

本文以原生结构焦煤孔隙系统为研究对象，综合运用压汞测试、纳米级工业CT扫描和数值仿真，精准重构煤微观孔隙空间结构，并对其进行定量表征，研究瓦斯解吸−扩散−渗流的演化过程，初步探讨煤孔隙空间结构对瓦斯运移的影响，旨在为煤层瓦斯运移理论提供新的见解。

1.   样品及实验

1.1   煤样处理

本文实验样品取自华晋焦煤有限责任公司沙曲一号煤矿四采区5号煤和6号煤，均属焦煤。其中，5号煤赋存于山西组下部，煤层平均厚度为3.30 m，属稳定大部可采煤层，瓦斯含量为11.16 m³/t，顶板为泥岩、砂质泥岩及极少量中−细粒砂岩，底板为粉砂岩、砂质泥岩、泥岩。6号煤赋存于太原组上部，上距5号煤平均距离为16.53 m，煤层平均厚度为0.55 m，属不稳定的局部可采煤层，瓦斯含量为10.05 m³/t，顶板为石灰岩，底板为泥岩或粉砂岩。原始煤样采集自工作面新鲜暴露的煤壁，每个煤层分别采集3块10 cm×10 cm×10 cm块状样品，保鲜膜包裹后运至实验室。块状煤样经破碎后打磨加工成5 mm×5 mm×10 mm的长方体，用于先后完成X射线CT（X−ray μCT）扫描实验和压汞实验，实验样品如图1所示。

图  1  实验样品

Figure  1.  Experimental samples

下载: 全尺寸图片 幻灯片

1.2   实验方法

X−ray μCT扫描实验用于获取并重构煤体微观孔隙结构，进而开展瓦斯渗流−扩散数值模拟。实验采用Phoenix v|tome|xs型多功能高分辨率CT扫描仪，该设备配备有高功率纳米焦点X射线管，空间对比能力和图像对比分辨率较强，如图2（a）所示。实验过程中选用的工作电压为130 kV，工作电流为70 mA，分辨率为0.5 μm，曝光时间为500 ms，采用局部扫描方式。煤样经过X−ray μCT扫描后，进行压汞测试。煤样经过二次干燥和脱气后，完成压汞测试，从而获取煤样孔隙率，用于煤微观孔隙精准重构。实验仪器为AutoPore IV 9500型压汞仪，如图2（b）所示。分析孔径范围为5 nm～800 μm，最大进汞压力为228 MPa。实验前样品经过12 h烘干，并在膨胀仪中抽真空6 h。

图  2  实验仪器

Figure  2.  Experimental instruments

下载: 全尺寸图片 幻灯片

1.3   CT数据处理

基于X−ray μCT扫描获得的原始CT数据存在像素噪声，影响后续孔隙重构的准确性，有必要进行降噪处理^[10]。常见的消除噪声方法有高斯滤波、均值滤波和中值滤波^[11]。高斯滤波主要用于处理符合正态分布的噪声，具有一定局限性^[12]。均值滤波能扩大目标像素点数，以8个像素点的灰度平均值作为目标点的灰度值，误差较大^[13]。中值滤波可在有效去除噪点的同时保护图像边缘，数据重建更为真实^[14]。因此，本文采用中值滤波对煤样的原始CT图像进行预处理，处理结果如图3所示。可看出，处理后的CT图像更清晰、噪点更少。

图  3  降噪前后CT图像对比

Figure  3.  CT image comparison before and after noise reduction

下载: 全尺寸图片 幻灯片

CT图像中，煤孔隙的识别依赖于阈值分割。常用的阈值分割法主要包括视觉判别、经验选取和模型反演^[15]。本文采用文献[16]提出的Bi−PTI模型识别煤孔隙结构，Bi−PTI模型对缺陷位置拓扑特征具有较好的识别能力，可真实还原孔裂隙拓扑体系的孔喉丰度和连通特性，对沙曲焦煤的适应性良好。Bi−PTI模型数学公式为

式中：Φ（G_m）为煤样孔隙率；G_m为孔隙的CT灰度阈值，HU；$\omega $为阈值权重；γ₁，γ₂分别为孔裂隙和基质骨架的阈值均值，HU；ζ₁，ζ₂分别为孔裂隙和基质骨架的灰度敏感系数。

利用图形处理软件提取CT数据中不同孔隙率对应的灰度阈值，然后用Bi−PTI模型进行拟合。将压汞孔隙率代入模型并反演孔隙灰度阈值，拟合结果如图4所示，拟合参数见表1。

图  4  孔隙灰度阈值的Bi−PTI拟合结果

Figure  4.  Bi-PTI fitting resules of pore gray threshold

下载: 全尺寸图片 幻灯片

表  1  Bi−PTI模型拟合参数

Table  1.  Fitting parameters of Bi-PTI model

煤样	γ1/ 104 HU	γ2/ 104 HU	ζ1/10−4	ζ2/10−4	$\omega $	相关 系数	Gm/HU
5号煤	5.32	3.81	3.110	0.94	0.44	0.99	10475
6号煤	1.31	2.45	2.100	0.81	0.65	0.99	12923

下载: 导出CSV 

| 显示表格

结合图4和表1可知，煤体CT数据孔隙率与灰度阈值的Bi−PTI拟合相关系数均超过0.99，拟合效果良好。5号煤的灰度阈值为10 475 HU，6号煤的灰度阈值略大于5号煤，为12 923 HU。

2.   数值模拟

2.1   数学模型

瓦斯在煤微观孔隙中的连续渗流过程符合低雷诺数和低速流的条件^[17]，因此渗流过程可采用Navier−Stokes（N−S）方程描述：

式中：ρ为瓦斯密度，kg/m³；v为流体速度，m/s；t为瓦斯扩散时间，s；$ \nabla $为拉普拉斯算子；P为瓦斯压力，MPa；μ为动力黏度，Pa·s；F为体积力，N。

由质量守恒定律可知，瓦斯在任意时刻流入和流出代表性体积（Representative Elementary Volume，REV）单元体截面的质量差与单元体内部的瓦斯质量变化量相等^[18]，故瓦斯在REV单元体内流动的连续性方程可描述为

依据前述假设条件，瓦斯密度ρ为常数，故可将式（3）简化为

受煤岩体孔隙结构非均质性影响，瓦斯在孔隙结构内部流动过程中，任意一点的速度大小和方向均不相同，随着空间位置和时间的变化而变化^[19]，因此可将式（4）转换为

联立式（3）、式（5）、式（6）可得到N−S方程的一般形式：

煤基质内瓦斯的解吸−扩散过程遵循含气体源项的菲克第二定律：

式中：C为物质的量浓度，与气体所处的孔隙空间坐标位置（x, y, z）及扩散时间t有关，mol/L；D为气体扩散系数，m²/s；S为气体源项。

对于瓦斯气体解吸−扩散耦合方程，式（7）中的气体源项可用瓦斯浓度随时间的变化情况来描述^[20]：

式中C_ad为煤基质内瓦斯物质的量浓度，mol/L。

式中：n_ad为煤岩基质内瓦斯气体量，mol；V_e为数值模拟计算中网格单位的体积，m³；s为网格单位所包含孔隙表面积，nm²；M为单位孔隙表面积吸附的气体量，mol；V_m为标准状况下的气体摩尔体积，22.4 L/mol。

单位孔隙表面积吸附的气体量M为^[12]

式中：V_ad为煤基质吸附的瓦斯体积，m³；N_solid为固体体素总数量；V_voxel为单位体素的体积，m³；ρ_true为煤岩密度，kg/m³；f为计算单元内部孔隙总表面积，m²。

假设瓦斯只吸附在孔隙内表面，当煤体内部任意位置体素点的灰度g（x, y, z）满足下列方程中的一个时，即可识别出孔隙空间内表面^[15]。孔隙部分体素的灰度为1，其余部分的灰度为0。

煤基质吸附的瓦斯体积V_ad可用朗格缪尔扩展型方程表示：

式中：a为甲烷的朗格缪尔体积常数，m³/kg； b为甲烷的朗格缪尔压力常数，MPa⁻¹。

联立式（7）—式（17）可得煤基质内部瓦斯的解吸−扩散耦合方程：

由于瓦斯具有可压缩性，其状态方程为

式中：n为气体物质的量，mol；Z为压缩因子，与瓦斯温度和压力有关；R为气体常数，J/（K·mol）；T为温度，K。

在考虑瓦斯气体压缩效应的情况下，煤基质内瓦斯气体解吸−扩散耦合方程为

2.2   模型及参数

基于重构的煤微观孔隙结构，借助COMSOL仿真软件开展煤体瓦斯微观解吸−扩散−渗流数值计算。模型选取尺寸为60×60×60像素的立方体单元。煤孔隙几何模型重建后需要优化调试，主要包括消除重合线条、倒角、缺口等，经优化后的四面体网格可直接用于后续数值模拟计算中^[13]。

数值计算所需的参数见表2，其中甲烷气体相关参数均按照标准状况下的气体属性设置，密度为0.714 kg/m³，动力黏度为1.80×10⁻⁵ Pa·s，温度为293 K。对于瓦斯解吸−扩散数值模拟，初始条件下，各煤体孔隙几何模型中孔隙内表面饱和瓦斯气体的初始压力为1×10⁻² Pa，且孔隙内部压力设置为0。瓦斯边界条件如图5所示。对于瓦斯渗流数值模拟，瓦斯渗流过程为单向渗流。其中，在5号煤中沿着x方向进行，其中孔隙外表面x+、x−分别设置为渗流进出口，对应瓦斯初始压力分别为1.1，0.1 MPa；在6号煤中沿着z方向进行，其中孔隙外表面z−、z+分别设置为渗流进出口，对应瓦斯初始压力分别为1.1，0.1 MPa。将其余4个孔隙外表面设置为自由滑移壁面，孔隙内表面设置为无滑移壁面。

表  2  数值模拟参数

Table  2.  Numerical simulation parameters

参数	5号煤	6号煤
R/（J·K−1·mol−1）	8.314	8.314
T/K	303	303
s/m2	5.6×10−9	5.9×10−9
a/（m3·kg−1）	0.0112	0.0112
b/MPa−1	1.86×10−7	1.86×10−7
Vm/（m3·kg−1）	0.024	0.024
f/m2	1×10−9	1.2×10−9
Nsolid	2.2×1017	2.1×1017
Vvoxel/m3	1×10−18	1×10−18
ρ/（kg·m−3）	1260	1260
Ve/m3	3.1×10−14	3.0×10−14
D/（m2·s−1）	3.60×10−12	3.60×10−12
P/MPa	1.0×10−2	1.0×10−2

下载: 导出CSV 

| 显示表格

图  5  瓦斯边界条件

Figure  5.  Gas boundary condition

下载: 全尺寸图片 幻灯片

3.   结果及讨论

3.1   煤微观孔隙结构定量表征

借助Thermo Sientific Avizo软件对数字岩心进行重构，重构体尺寸均为500×500×1 000像素。CT数据体的REV单元是探索孔裂隙精细结构特征的一种有效手段，通常用于地球物理领域的数字煤岩结构研究^[20]。在重构的数字岩心结构中，选取感兴趣的4个目标点，选取的REV单元尺寸均设置为200×200×200像素，位置如图6所示。

图  6  数字岩心重构结果

Figure  6.  Digital core reconstruction results

下载: 全尺寸图片 幻灯片

基于各煤样最佳孔隙灰度阈值，通过Avizo构建REV表征单元内部的微观孔隙空间结构，如图7所示。分析可知，5号煤含有大量孔隙，A2内部有多条连通的裂隙，A4裂隙含量最少，均存在较多孔隙团。6号煤含有大量孔隙，4个REV表征单元内部的裂隙表现出不同形态，B1内部2条呈V字形的连通裂隙贯穿整个表征单元，B4下方有多条呈片状横向平行分布的裂隙。

图  7  REV表征单元孔隙空间结构

Figure  7.  Representative elementary volume(REV) characterisation of unit pore space structure

下载: 全尺寸图片 幻灯片

各煤样REV单元微观孔隙等效直径的分布规律如图8所示。分析可知，5号煤4个REV单元的孔隙等效直径分布范围较为杂乱，4 μm的孔隙均占比最高，其在A4中的占比达26%。A2的等效直径分布范围最窄，为1～8 μm，等效直径为4 μm的孔隙占比约为21%。与5号煤相比，6号煤微观孔隙等效直径分布范围、优势孔隙分布均非常相似，表明5号煤的非均质性较强，其中等效直径为3.5 μm的孔隙占比约为23%，相对应的孔隙累计占比为57%。可见，5号煤和6号煤的孔隙等效直径分布规律大致相同。

图  8  孔隙等效直径分布规律

Figure  8.  Distribution law of pore equivalent diameter

下载: 全尺寸图片 幻灯片

煤微观孔隙体积分布规律如图9所示。分析可知，5号煤的孔隙体积主要分布在1～800 μm³，平均值为419 μm³，A1中的优势孔最大体积为530 μm³，相对应孔隙体积在1～530 μm³的孔隙累计占比达63%。6号煤的孔隙体积分布规律与5号煤较为相似，不同之处在于6号煤中孔隙体积为260 μm³的孔隙占比较高，表明6号煤体积的小孔隙占主导地位。可见，2种煤样的孔隙等效体积分布范围不存在明显差异。

图  9  孔隙体积分布规律

Figure  9.  Distribution law of pore volume

下载: 全尺寸图片 幻灯片

3.2   瓦斯解吸−扩散过程

瓦斯在煤微观孔隙REV单元中的压力分布情况如图10所示。可看出，随着瓦斯解吸−扩散时间的增加，瓦斯压力自REV单元中心至边缘逐渐降低，同时在不同时间及位置上瓦斯压力分布规律差异明显。造成瓦斯压力分布差异性的原因在于各REV单元中孔隙与喉道的半径、长度、形状、连通性能不同。

图  10  瓦斯压力空间分布

Figure  10.  Space distribution of gas pressure

下载: 全尺寸图片 幻灯片

为进一步分析瓦斯解吸−扩散过程中瓦斯压力在不同微观孔隙系统中的分布情况，对切片中瓦斯压力场进行定量分析。分别截取各REV单元体在x轴5，10，15，20，25，30 μm位置的y−z平面，观察不同位置切片压力场随时间变化情况，如图11所示。

分析可知，在同一REV切片中，瓦斯在解吸−扩散时由压力大的中心位置向压力较小的边缘位置扩散。同一时间，在不同切片中，中心位置（10，15，20，25 μm）切片瓦斯压力相对较高，而边缘位置（5，30 μm）切片瓦斯压力相对较低。可见，瓦斯气体在孔隙系统中由内而外地逐步完成解吸−扩散。此外，尺寸较大的连通孔隙能够为气体分子提供不同方向上的多元化运动空间，进而提高瓦斯解吸−扩散速率。但连通孔隙的小尺寸结构限制了气体分子空间扩散的广度，极大阻碍瓦斯的解吸−扩散进程。

图  11  REV单元y−z截面孔隙压力分布

Figure  11.  Pore pressure distribution in y-z section of REV unit

下载: 全尺寸图片 幻灯片

3.3   瓦斯渗流过程

在煤三维孔隙空间结构中，由于孔隙压力分布、孔隙通道连通性的差异性，造成瓦斯气体在煤岩体孔隙空间中的渗流速度大小和方向发生明显变化。为研究瓦斯渗流速度变化规律，保持瓦斯初始压力和扩散系数不变，分别截取各REV单元不同时间节点（1，10，20，40 s）的速度场流线分布情况，如图12所示，图中不同颜色的箭头指示线条表示不同时间节点瓦斯在煤岩体中的渗流过程，颜色由深蓝色逐渐过渡到深红色表示瓦斯流动速度不断增大。分析发现，5号煤中瓦斯渗流速度流线的分布范围、密度较6号煤更大，表明5号煤的连通孔隙数量更多，瓦斯渗流通道更加丰富。5号煤和6号煤渗流过程中瓦斯在各煤体三维孔隙空间内的渗流速度大小和方向表现出显著的差异性。在诸如5号煤的强非均质连通孔隙结构中，瓦斯渗流分散而高效，能通过广泛沟通煤基质完成气体由扩散到渗流的转变，提升瓦斯传质效率。然而，在诸如6号煤的弱非均质连通孔隙结构中，气体渗流路径单一、流线集中，渗流传质阻力较大，气体分子由扩散到渗流的转变效率低，不利于瓦斯高效运移。

图  12  瓦斯渗流速度流线分布

Figure  12.  Velocity streamline distribution of gas seepage

下载: 全尺寸图片 幻灯片

4.   结论

1） 从微观角度看，瓦斯在解吸−扩散时由孔隙中心位置向边缘位置进行。中心位置瓦斯压力相对较高，而边缘位置瓦斯压力相对较低。大尺寸连通孔隙结构能够为气体分子提供不同方向上的多元化运动空间，削弱尺寸效应对气体分子扩散广度的影响，进而促进瓦斯解吸−扩散速率。

2） 瓦斯在微观孔隙空间结构内的渗流速度大小和方向表现出明显的差异性，这主要是由于不同煤体微观孔隙系统存在高度非均质性。强非均质连通孔隙结构中的渗流传质阻力较小，能通过广泛沟通煤基质完成气体由扩散到渗流的高效转变，提升瓦斯传质效率。研究结果能从微观角度丰富煤体瓦斯运移理论，为瓦斯抽采工程实践提供一定理论基础。