基于粗糙离散裂隙网络模型的煤层气流动特性研究

李文彬 ，张书金 ，衡献伟 ,3

（1.贵州省煤矿设计研究院有限公司，贵州 贵阳 550025；2.贵州省矿山安全科学研究院有限公司，贵州 贵阳 550025；3.贵州大学 矿业学院，贵州 贵阳 550025）

煤层气是一种以甲烷为主要成分的清洁高效能源，对其开采不仅可以减少温室气体排放、接替常规天然气能源，同时也可减少煤矿灾害的发生，保障矿井安全生产等[1-2]。目前的研究认为[3-5]，煤体是一种多孔介质材料，其内部广泛分布的裂隙结构是气体运移的重要通道，对产气量有显著影响。煤体内的裂隙分布呈现明显的随机性和无序性，现有模拟储层裂隙的方法主要分为等效连续介质模型[6]和离散裂隙网络模型[7]，2 种方法各有优缺点，在实际应用时，等效连续介质模型对裂隙的假设过于理想化，无法描述真实裂隙的细节，而离散裂隙网络模型恰好弥补了这一不足，可以充分刻画裂隙的真实细节和分布特征，因而被广泛应用。其中，采用Monte Carlo 方法生成离散裂隙网络模型是目前最多的方法，许多学者利用该方法分析了煤层气在裂隙内的流动规律。然而传统离散裂隙网络模型仍然是假设裂隙为平滑裂隙，对真实裂隙细节的描述仍然存在不足。而真实煤体内的裂隙绝大多数是凹凸不平的，气体在其内的流动也不同于平滑裂隙网络模型。因此，采用一种新的粗糙离散裂隙网络建模方法建立裂隙网络几何模型，该方法相比Monte Carlo 方法可充分考虑裂隙粗糙特征，因而更加符合实际情况。为此，基于建立的裂隙网络模型，结合数值模拟分析煤层气在裂隙内的流动规律，探讨不同裂隙结构特征参数对煤层气流动特性的影响。

1 粗糙离散裂隙网络模型

实际煤体内的裂隙往往是弯曲的，具有一定的粗糙特性，且裂隙相互交叉形成裂隙网络，组合形式多为网状和不规则网状，渗透率相对较好。为了更真实地建立符合实际工程情况的粗糙离散裂隙网络模型，WU 等[8]和XIA 等[9]等提出了基于像元概率分解算法建立粗糙离散裂隙网络模型的方法，根据其理论，在像元空间中，任意1 条裂隙都可以由点集合的形式组成，即：

式中：GF为裂隙点集合系统；Di为裂隙点的成长方向；Mi为成长概率；ti为成长步长。

根据像元概率分解算法的基本理论，借助MATLAB 软件编制二维粗糙离散裂隙网络生成程序，以实现像元空间范围内多条裂隙的随机生成。由式（1）可知，如果确定了裂隙的成长方向、成长概率、成长步长，裂隙点集合系统也就确定了。为了更好地说明其具体实现过程，假定裂隙的成长概率M1={0.67, 0.33, 0, 0, 0, 0, 0, 0}，为二项随机，成长步长ti=15，则1 条在像元空间范围[100 100]中裂隙点集合，粗糙离散裂隙网络生成过程如图1。

图1 粗糙离散裂隙网络生成过程Fig.1 Generation process of rough discrete fracture network

当需要生成具有同样特征的150 条裂隙点集合时，原理与生成1 条裂隙一致，只需要在程序中设置裂隙总数，裂隙在像元空间中随机分布，结果如图1（b）。由于真实煤岩体中的裂隙相互交错，为生成更为符合真实情况的裂隙网络，这里裂隙的分布特征参考文献[10-11]获得的统计规律，即裂隙长度幂律符合分布特征，倾角分布符合Fisher分布特征。在图1（b）中增加150 条裂隙，成长概率M2={0, 0, 0, 0, 0, 0, 0.33, 0.67}，成长步长ti=15，因此，具有不同成长概率，数量为300 条的二维像元空间中的裂隙网络集合系统如图1（c）。

2 粗糙离散裂隙网络渗透率模型

煤体内的裂隙相互交叉形成裂隙网络，对于任意复杂的裂隙网络均是由一条条单裂隙所组成。目前描述流体在单裂隙内的流动多采用经典的平行板立方定律[12]，即

式中：q为流过裂隙的总流量，m3/s；α为裂隙开度，m；l为裂隙长度，m；θ1、θ2分别为裂隙方位角和倾角，（°）；μ为流体动力黏度，μPa·s；Δp为压差，Pa；L0为特征单元体长度，m。

MIAO 等[11]借鉴多孔介质模型中孔隙网络渗透率模型的推导方法，在上述平行板立方定律的基础上，结合分形结合推导得到了适用于随机分布的裂隙网络渗透率模型，即：

式中：Kf为裂隙网络渗透率，m2；β为比例系数，取值范围一般为10-3～10-1；Dl为裂隙网络的分形维数；ϕf为裂隙率；lmax为最大裂隙长度，m。

然而，上述模型是假设煤岩体内的裂隙为平滑裂隙，没有考虑裂隙的粗糙或迂曲度特征。针对上述平行板裂隙网络模型存在的不足，LUO 等[12]改进了上述模型，在上述模型中增加了裂隙的粗糙特征，即：

式中：DTf为裂隙迂曲度分形维数。

式（4）考虑了裂隙的粗糙特征，因而也更加符合实际情况。因此，在后续章节将主要结合该渗透率模型和粗糙离散裂隙网络建模方法分析气体在裂隙内的流动特性。

3 煤层气流动数值模拟

由于实际煤层气的流动是一个极其复杂的过程，涉及的影响因素众多，研究仅分析裂隙结构特征参数的变化对煤层渗透性的影响。因此，研究重点在于关注不同裂隙结构特征对渗流的影响，故不考虑其它因素的影响，即假设气体为单相流，遵循达西定律，忽略应力、温度等因素的变化对渗流的影响。

3.1 边界条件和网格划分及气体流动控制方程

为采用像元概率分解算法建立的粗糙离散裂隙网络几何模型，然后将建立的裂隙网络模型导入到COMSOL 多物理场耦合软件中进行煤层气流动数值计算，煤层气流动数值模型如图2。模型尺寸为100 mm×100 mm，在模型的左侧设置压力入口，入口压力为1 000 Pa，模型右侧为出口压力，压力为0，将模型下边界固定，上边界为不渗流边界。然后对模型进行网格划分，这里采用自由三角形网格对其进行划分，对特别关注的裂隙区域进行网格加密处理。三角形单元数量为57 617 个，边单元数量为373 个，最小单元质量0.393 1 mm，平均单元质量0.824 6 mm。

图2 煤层气流动数值模型Fig.2 Numerical modeling of coalbed methane flow

数值计算过程中，煤层气的流动控制方程采用COMSOL 内置的Navier-Stokes 方程[13]，即：

式中：V为气体流速，m/s； ∇为哈密顿算子。

3.2 模拟结果

设定边界条件及相关参数后，通过达西稳态流求解控制方程（式（6）），可得到气体压力分布云图以及气体在裂隙内的流速分布云图，气体压力分布云图如图3，气体在裂隙网络中的流速分布云图如图4。

图3 气体压力分布云图Fig.3 Cloud diagram of gas pressure distribution

图4 气体在裂隙网络中的流速分布云图Fig.4 Cloud image of gas velocity distribution in fracture network

由图3 可知，由于煤体内裂隙的存在，使得压力分布呈现非均匀分布特征。由图4 可知，气体的流动主要发生在煤体内相互连通的裂隙网络内，裂隙的分布特征影响着气体的流动特性。

4 参数分析及讨论

由式（4）可知，裂隙网络渗透率主要与裂隙率、分形维数、迂曲度分形维数、裂隙倾角、最大裂隙长度以及比例系数等参数有关。因此，为了进一步分析裂隙特征参数对煤层渗透性的影响规律，采用像元概率分解算法建立具有不同结构特征的粗糙离散裂隙网络模型，结合数模模拟方法分析不同裂隙特征参数对煤层渗透性的影响，裂隙特征参数见表1。

表1 裂隙特征参数Table 1 Parameter values of fracture characteristics

4.1 裂隙率和分形维数

裂隙率和分形维数是表征裂隙复杂程度的重要参数，其中裂隙率为裂隙的总面积与测量煤岩样品的比值，其值越大表明裂隙数量越多，越容易形成更多的流动通道。同样，分形维数也是表征裂隙复杂程度常用的定量表征参数，是裂隙占有度量空间的多少，裂隙分形维数越大，裂隙分布越密集，越容易形成更多的渗流通道。为具体分析裂隙率和分形维数的变化对煤层渗透性的影响，固定其它裂隙参数不变（其中每个案例中裂隙的倾角为45°、迂曲度分形维数1.015、比例系数为0.014 以及最大裂隙长度为17 mm），采用像元概率分解算法生成具有不同裂隙结构参数的粗糙离散裂隙网络几何模型，裂隙率和分形维数的变化范围分别为0.1～0.25、1.55～1.85。

裂隙率与裂隙网络模型出口流速的关系如图5，分形维数与裂隙网络模型出口流速的关系如图6。

图5 裂隙率与裂隙网络模型出口流速的关系Fig.5 Relationship between fracture rate and outlet velocity of fracture network model

图6 分形维数与裂隙网络模型出口流速的关系Fig.6 Relationship between fractal dimension and outlet velocity of fracture network model

从图5 可以看出：随着裂隙率的增大，模型出口流速也随之增大，通过对模拟数据拟合发现，二者基本呈现较好的幂律增长关系。

由图6 可知：裂隙网络分形维数的增加使得模型出口速度也随之增大。对其数据进行拟合发现，分形维数与出口流速之间同样展现了较好的幂律增长关系。

4.2 迂曲度分形维数和最大裂隙长度

分析裂隙迂曲度分形维数的变化对煤层渗透性的影响所需参数见表1。固定其它裂隙参数不变，迂曲度分形维数和最大裂隙长度的变化范围分别为1.0～1.02 和11～26 mm。迂曲度分形维数与裂隙网络模型出口流速的关系如图7，最大裂隙长度与裂隙网络模型出口流速的关系如图8。

图7 迂曲度分形维数与裂隙网络模型出口流速的关系Fig.7 Relationship between tortuosity fractal dimension and outlet velocity of fracture network model

图8 最大裂隙长度与裂隙网络模型出口流速的关系Fig.8 Relationship between maximum fracture length and outlet velocity of fracture network model

由图7 可知：迂曲度分形维数与流速之间呈现了负相关关系，即随着迂曲度分形维数的增大，裂隙网络模型出口流速逐渐减小，这是因为迂曲度分形维数表征的是裂隙的弯曲程度，迂曲度分形维数值越大，说明裂隙越弯曲（或粗糙），则气体在裂隙内的流动阻力也就越大，因此，出口流速逐渐减小。

由图8 可知：随着裂隙长度的逐渐增大，模型出口流速也随之增大。这是因为裂隙最大长度的增加意味着裂隙之间更易形成网络，使得裂隙的连通性整体性增强，进而影响着煤层渗透性。

4.3 比例系数和裂隙倾角

同样，为分析比例系数和裂隙倾角对煤层渗透性的影响，固定裂隙长度、迂曲度分形维数等参数不变，比例系数变化范围为0.006～0.08，倾角变化范围为10°～85°。比例系数与裂隙网络模型出口流速的关系如图9，裂隙倾角与裂隙网络模型出口流速的关系如图10。

图9 比例系数与裂隙网络模型出口流速的关系Fig.9 Relationship between proportional coefficient and outlet velocity of fracture network model

图10 裂隙倾角与裂隙网络模型出口流速的关系Fig.10 Relationship between fracture inclination angle and outlet velocity of fracture network model

由图9 可知：随着比例系数的增加，裂隙网络模型出口流速逐渐增大。其原因在于，当裂隙长度不变时，比例系数与裂隙的开度成正相关关系，比例系数的增加意味着裂隙开度的增大，因此气体在裂隙内的流动更为容易，导致裂隙网络模型出口流速增大。此外，由图中数据的变化可知，比例系数的增大对裂隙网络流速的影响最为显著。

由图10 可知：随着倾角的逐渐增大，裂隙网络模型出口流速逐渐减小，这是由于较小裂隙倾角将更有利于气体流动，而倾角越大，气体流动阻力越大，可见煤层中裂隙倾角分布方向决定裂隙网络的渗流方向。

综上，为了更好地观察裂隙特征参数与流速之间的关系，对模拟结果进行非线性拟合，拟合结果表明，裂隙率、分形维数、最大裂隙长度以及比例系数与裂隙网络出口流速之间呈现较好的幂律增长关系，而迂曲度分形维数、裂隙倾角与裂隙网络出口流速之间呈现幂律降低的关系。

5 结 语

1）详细介绍了粗糙裂隙离散裂隙网络的建模理论及其实现过程，并基于该方法建立了更符合工程实际的粗糙离散裂隙网络模型。

2）由于煤体内裂隙网络的存在，气体压力分布呈现非均匀分布特征，气体的流动主要发生在相互连通的裂隙网络中，煤层渗透率主要为裂隙网络所控制。

3）裂隙率、分形维数、最大裂隙长度以及比例系数与裂隙网络出口流速之间呈现较好的幂律增长关系，而迂曲度分形维数、裂隙倾角与裂隙网络出口流速之间呈现幂律降低的关系；在这些参数中，比例系数、裂隙倾角和最大裂隙长度对煤岩渗透率的影响最为显著。