王世斌，王 刚，，陈雪畅，范酒源，迟利辉

（1.山东科技大学能源与矿业工程学院，山东青岛 266590；2.山东科技大学安全与环境工程学院，山东青岛 266590）

水力压裂作为一种有效的强渗-增润手段被广泛应用于煤层气开采过程中，通过向煤层中打入压裂钻孔，注入高压水对煤体进行压裂，提高煤层的渗透率，从而对煤层中的瓦斯进行抽采[1-6]。水力压裂及瓦斯抽采的效果往往通过煤层的裂隙扩展、瓦斯抽采流量和抽采作用范围等指标进行评价，因此，为了达到理想的压裂和抽采效果，在实际应用前通常会采用数值模拟软件进行模拟，从而确定实际应用的工艺参数。

针对煤层水力压裂和瓦斯抽采的数值模拟研究，国内外研究学者开展了较多的研究工作。从常奎[7]利用FLAC3D数值模拟分析了顶板水力压裂裂缝扩展演化特征，讨论了压裂孔直径和压裂压力对裂缝发育的影响规律；闫炎等[8]基于损伤力学与流固耦合理论，建立了水力压裂过程中射孔段水泥环裂缝扩展的数值模型，对水平井段压裂工况下固井水泥环径向裂缝沿井筒轴向的扩展过程进行了数值模拟；付海峰等[9]通过三维离散格子模拟方法建立了等尺度射孔水力裂缝起裂扩展全三维数值模型，揭示了不同射孔模式下的裂缝起裂规律；张玉等[10]以多孔介质流体渗流和围岩应力耦合理论为基础，提出一种基于有限容积法（FVM）的水力耦合作用下射孔围岩水力压裂破裂数值模拟方法，该方法丰富了水力压裂破裂机理的研究，对实际工程设计提供重要的参考；刘嘉等[11]基于多孔弹性理论和能量最小化原理建立水力耦合的相场模型，采用交错策略的分离式耦合方法进行求解，并基于三维水力压裂模拟结果进行分析，验证了该方法对于模拟不同地应力条件下的水力裂纹扩展问题的适用性；闫晓等[12]提出了统一管道-界面单元法（UP-IEM），模拟可渗透性孔隙裂隙介质中的水力压裂，通过3 个经典数值模型验证了UP-IEM 模拟流体流动、裂纹扩展以及水力压裂的准确性及适用性；杨晨旭等[13]针对水力压裂过程中所涉及到的岩石变形与裂缝扩展、缝内流体的流动及滤失等力学问题，基于有限元法建立了全三维全耦合水力压裂数值模型，并采用数值模拟方法对模型进行验证；郭欣等[14]以瓦斯渗流理论为基础，根据达西定律和质量守恒定律，利用COMSOL Multiphysics 数值模拟软件建立考虑煤岩变形控制方程的气固耦合数学模型，研究相关因素对瓦斯有效抽采半径大小的影响；郝长胜等[15]运用FLUENT 软件模拟分析布置于平行回风巷且孔间距为10 m 的顶板走向长钻孔束在不同参数条件下的瓦斯抽采效果，为寺河矿5309 工作面上隅角瓦斯超限问题的解决提供了一种方法；贾毅超等[16]采用FLAC 数值模拟软件对高位钻孔布置优化进行了研究，并在实际工程应用中提高了瓦斯抽采浓度，减小了上隅角和回风流中的瓦斯浓度；汪腾蛟等[17]应用有限体积法分别离散了考虑温度变化下的采空区瓦斯抽采多场耦合二维数学模型，基于Microsoft Visual Basic 编制了计算机解算程序，利用Tecplot软件对求解结果进行可视化，研究了瓦斯抽采前后的采空区瓦斯压力、瓦斯体积分数、氧体积分数、固体温度和气体温度分布情况及前后变化趋势；王刚等[18]建立急倾斜煤层钻孔周围气体渗流单元模型，利用数值模拟分析了急倾斜煤层瓦斯钻孔有效抽采范围随时间的变化规律；王广宏[19]利用COMSOL Multiphysics 数值模拟软件对榆树田煤矿的定向长钻孔瓦斯抽采工艺技术进行了研究，以提高瓦斯抽采效率，确保矿井安全生产；张育磊[20]以万峰煤矿1201 工作面为工程背景，采用FLAC3D数值模拟软件，研究煤层开采后对邻近煤层的影响及合理的瓦斯抽采钻孔布置位置；张智渊等[21]基于建立的复合储层和单一储层瓦斯抽采流固耦合模型，模拟研究了不同倾角的穿层钻孔对复合储层和单一储层的瓦斯抽采效果，以及多钻孔耦合下的穿层钻孔瓦斯抽采效果。

从上述研究中可以看出，绝大多数学者仅针对煤层水力压裂和瓦斯抽采中1 种工艺方式进行了数值模拟，未对水力压裂作用后的煤体开展瓦斯抽采数值模拟，也即未对2 种工艺数值模拟间的“链接”进行构建，从而无法达到先压裂后抽采的模拟要求。针对这一问题，提出了一种基于PFC2D水力压裂和COMSOL 瓦斯抽采结合的数值模拟方法，将PFC2D水力压裂后裂隙扩展的图像导入到网格绘制的图像处理软件中，将图像进行像素级优化，按照像素点将图像中的裂隙进行绘制，然后将绘制的图像导入COMSOL 中进行瓦斯抽采数值模拟，从而实现了煤层先压裂后抽采的数值模拟研究，研究内容对煤层气开采具有重要的指导意义。

1 煤层水力压裂PFC2D 数值模拟研究

1.1 PFC2D 数值模拟模型及参数

模型构建采用PFC2D软件建立200 mm × 200 mm 的二维模型，模型的中心位置设定半径7.0 mm的水力压裂孔，模型4 个边界为固定边界，即对模型边界的颗粒进行固定。水力压裂模型示意图如图1。

PFC2D模拟时的宏观和细观参数如下[22-23]：

1）宏观力学参数。单轴抗压强度σmc为32.22 MPa，弹性模量E 为1.81 GPa，泊松比ν 为0.25，本体抗拉强度σmt为2.83 MPa。

2）颗粒力学参数。颗粒最小半径Rmin为2.50 mm，粒径比Rmax/Rmin为1.66，颗粒密度ρ 为1 280 kg/m3，孔隙率为9%，阻尼β 为0.50 N/（s·m）。

3）接触模型。细观弹性模量Ec为1.00 GPa，刚度比kn/ks为2.30，摩擦因数μ 为0.50，细观抗拉强度σc为11 MPa，黏聚力c 为20 MPa，摩擦角φ 为37.50°，法向临界阻尼比βn0.57。水力压裂模型示意图如图1。

图1 水力压裂模型示意图Fig.1 Schematic diagram of hydraulic fracturing model

除水力压裂的注水压力外，其余参数保持不变，模拟20 MPa 和40 MPa 的不同注水压力下煤体水力压裂，模拟时仅设置注水钻孔边界处的压力为20 MPa 和40 MPa，且在注水过程中注水压力不发生变化，不存在钻孔注水压力上升的过程，始终保持20 MPa 和40 MPa，且裂纹扩展后钻孔注水压力也为20 MPa 和40 MPa，直到模拟时注水压力不发生改变。计算步数为10 000，每一步的时间间隔为0.001 s，对应总时间为10 s。

1.2 水力压裂实现方式

PFC2D进行水力压裂时，为了定量计算流体交换时的流量，假设流体通道是由相邻两颗粒的接触点组成的1 个平行板通道，长度为L，开度为a，厚度为单位厚度，垂直于x-y 平面[24]。则2 个孔隙之间的流量可以通过Hagen-Poiseuille 方程表示：

式中：q 为流量，m3/s；k 为渗透系数，m/（Pa·s）；△p 为2 个孔隙网格之间的压力差，Pa。

在颗粒之间接触处，假设存在1 个初始开度a0，初始开度的存在允许即使两颗粒之间紧密接触的情况下其构成的流体通道也有流体流动，从而保证了材料的基质渗透率。流体通道的开度a 依赖于颗粒之间的接触力，当两颗粒之间的法向接触力为压应力时，流体通道开度a 为：

式中：F 为现在两颗粒间的压缩力，N；F0为管道开度降低为初始开度一半时的压缩力，N。

当两胶结颗粒处于张拉状态，或者两颗粒之间的黏结已破坏，颗粒在接触点断开时，开度a 为：

式中：d 为两颗粒之间的距离，m；R1、R2分别为两颗粒半径，m；λ 为1 个无量纲乘子，对于大部分的模型而言，颗粒粒径都比实际的颗粒粒径大得多，计算得到的开度也会偏大，因此，λ 常常取1 个小于1 的常数来获得1 个合理的开度。

在△t 的时间步里面，由于流体流动导致的孔隙流体压力变化由流体的体积压缩模量计算。考虑某个孔隙，其有N 条流体通道，在△t 的时间步里面，其流体总流量为∑q，则孔隙流体压力的变化△p 可表示为：

式中：Kf为流体的压缩模量，Pa；Vd为“域”，即孔隙体积，m3；△Vd为孔隙体积变化，m3。

在PFC2D中流体与颗粒间的耦合作用方式主要有以下3 种：①通过接触的张开与闭合或接触力的变化实现通道孔隙的变化；②通过改变研究区域的力学特性来改变其中的压力；③区域孔隙压力对其内部颗粒有推移作用。

采用显式方法求解，在所有域内交替应用流动方程和压力方程。保证模型运行稳定的条件就是水流引起的压力变化必须小于扰动压力，当两者相等时可求出临界时间步长△t 为：

式中：N 为1 个“域”所连接的管道数；r 为1 个“域”周围颗粒的平均半径，m。

1.3 水力压裂结果

不同注水压力的压裂范围可通过压裂后颗粒的位移云图来反映，不同注水压力的裂隙扩展规律可通过裂隙扩展图表现。2 种注水压力下颗粒位移云图如图2，2 种注水压裂状态下的裂隙扩展如图3。

图3 2 种注水压力下注水压裂裂隙扩展图Fig.3 Fracture propagation diagrams under two kinds of water injection pressure

从图2 中可以看出：2 种注水压力下的压裂范围有所不同；20 MPa 的注水压力时，仅有钻孔周围的小部分颗粒具有位移，说明20 MPa 注水压力下，注水压裂的范围较小；40 MPa 的注水压力时，钻孔周围的许多颗粒均产生了位移，且颗粒的位移量沿钻孔径向逐渐减小，表明40 MPa 注水压力下，注水压裂的范围较大，同时说明增加注水钻孔的注水压力可以提高水力压裂的压裂效果。

图2 2 种注水压力下颗粒位移云图Fig.2 Cloud diagrams of particle displacement under two kinds of water injection pressure

从图3 中可以看出：裂隙从钻孔附近沿径向逐渐扩展，不同注水压力下生成的裂隙数量和范围有所不同；20 MPa 注水压力生成的裂隙数量较少，仅有8 条裂隙；而40 MPa 注水压力生成的裂隙数量较多，有78 条裂隙；20 MPa 注水压力下裂隙扩展延伸的范围较小，压裂效果较差；40 MPa 注水压力下裂隙扩展延伸的范围较大，压裂效果较好；同时也表明增加钻孔的注水压力可以增加水力压裂的裂隙扩展延伸，提高水力压裂的效果。

2 煤层瓦斯抽采COMSOL 数值模拟

2.1 COMSOL 模型

将PFC2D水力压裂后裂隙扩展的图像导入到网格绘制的图像处理软件中，将图像进行像素级优化，按照像素点对图像中的裂隙进行绘制，将像素长度与实际长度进行匹配，然后将绘制的图像导出为COMSOL 能够识别的网格文件，最后将网格文件导入到COMSOL 数值模拟软件中，瓦斯抽采几何模型构建图如图4。

图4 瓦斯抽采几何模型构建图Fig.4 Geometric model construction diagrams of gas extraction

2.2 控制方程

模拟采用COMSOL 软件中的达西定律模块，控制方程为：

式中：εp为孔隙率；ρ 为瓦斯密度，kg/m3；u 为瓦斯的达西渗流速度，m/s；Qm为源汇项，kg/（m3·s）；k为渗透率，m2；μ 为瓦斯动力黏度，Pa·s；△pg为压力梯度，Pa/m；t 为时间变量，s。

将处于渗流运动中的瓦斯视为理想气体，满足理想气体状态方程。该方程表示平衡状态下，理想气体的压强、体积、物质的量以及温度之间的关系为：

式中：Mg为瓦斯的摩尔质量，kg/mol；pg为瓦斯压力，Pa；R 为气体摩尔常数，J/（mol·K）；T 为绝对温度，K。

2.3 边界条件及参数

瓦斯抽采数值模拟时的边界条件如图5。

图5 瓦斯抽采边界条件示意图Fig.5 Schematic diagram of boundary conditions of gas extraction

模型四周边界无瓦斯流动，煤层的初始瓦斯压力为p0，钻孔的瓦斯抽采压力为pout，由于PFC2D构建注水钻孔时采用的是删除颗粒的命令，当颗粒有不存在删除区域中的部分时，颗粒就不会被删除，因此会使构建的钻孔直径有所变化，因此在抽采时将钻孔直径设置为5 mm。各参数具体数值如下：①煤体孔隙率φ0为9.00%；②煤体孔隙渗透率k0为9×10-16m2；③煤体裂隙渗透率kl为9×10-13m2；④初始瓦斯压力p0为1 MPa；⑤抽采负压pout为0.2 MPa；⑥瓦斯动力黏度μ 为1×10-5Pa·s；⑦瓦斯密度ρ 为0.716 kg/m3；⑧气体摩尔常数R 为8.314 5 J/（mol·K）；⑨参考温度T 为293.15 K；⑩瓦斯摩尔质量Mg为16 g/mol。

2.4 瓦斯抽采结果

2 种压裂条件下瓦斯抽采10 s 的压力分布云图如图6。20 MPa 注水压裂煤层瓦斯抽采压力分布云图如图7，40 MPa 注水压裂煤层瓦斯抽采压力分布云图如图8。

图6 抽采10 s 瓦斯压力分布云图Fig.6 Gas pressure distribution diagrams of 10 s

图7 20 MPa 注水压裂煤层瓦斯抽采压力分布云图Fig.7 Diagrams of gas drainage pressure distribution in 20 MPa water injection fracturing coal seam

从图6 中可以看出：瓦斯压力沿着钻孔径向向中心逐渐减小，这是因为煤层中的瓦斯被抽采，靠近钻孔周围的煤层最先被抽采出，瓦斯压力减小，煤层中瓦斯压力分布不均匀，从而使得高压瓦斯不断向低压处运移。

瓦斯抽采范围可以通过瓦斯压力分布云图反映出，瓦斯压力产生变化的地方表明钻孔的瓦斯抽采范围能够达到此处。从图7 和图8 可以看出：随着抽采时间的增加，瓦斯抽采的范围逐渐增大；当抽采时间相同时，40 MPa 与20 MPa 水力压裂的煤层相比，前者的瓦斯抽采范围比后者大，表明增加水力压裂的注水压力不仅增加了煤层中裂隙的扩展范围，同时也提高了煤层的瓦斯抽采范围。

采用COMSOL 中的线积分对钻孔的瓦斯抽采流量进行计算，钻孔瓦斯抽采流量变化如图9。由于本文COMSOL 瓦斯抽采数值模拟采用的是二维几何模型，对于抽采流量的计算是采用线积分的方式，因此抽采流量单位为m2/s。

图9 钻孔瓦斯抽采流量变化Fig.9 Change of gas extraction flow rate from borehole

从图9 中可以看出，随着抽采时间的增加，抽采流量逐渐降低，在抽采开始后的10 s 内，钻孔瓦斯抽采流量下降的最快，20 MPa 水力压裂下瓦斯抽采流量从0.006 6 m2/s 降低为0.001 6 m2/s；40 MPa 水力压裂下瓦斯抽采流量从0.010 1 m2/s 降低为0.002 0 m2/s，这是因为此时煤层的瓦斯压力与钻孔抽采压力之间形成的压力梯度最大，瓦斯运移的速度快，从而使得钻孔的瓦斯抽采流量大。随着抽采的进行，煤层中的瓦斯含量和瓦斯压力均有所下降，瓦斯压力与抽采钻孔之间的压差逐渐减小，因此会造成瓦斯流量的降低。同时，2 种压裂条件下的瓦斯抽采流量也不相同，40 MPa 水力压裂后的煤层瓦斯抽采流量比20 MPa 水力压裂后的煤层大，表明增加水力压裂的注水压力在增加煤层裂隙扩展的同时也增加了瓦斯抽采的流量。

文献［25］的实际工程应用中也显示水力压裂后钻孔的瓦斯抽采流量有所增加，表明水力压裂可以提高煤层的透气性，提高煤层的瓦斯抽采效果，进一步验证所提出的数值模拟方法具有可行性。

3 结 语

对煤层开展了PFC2D水力压裂数值模拟，分析了不同注水压力下的压裂范围与裂隙扩展规律，并采用一种图像处理方法将压裂后的模型导入到COMSOL 软件中进行压裂后的瓦斯抽采数值模拟研究，分析了不同压裂条件下的瓦斯抽采范围与瓦斯抽采流量。

1）提出的几何模型图像像素级处理方法成功将PFC2D水力压裂后的裂隙几何模型导入到COMSOL软件中，完成了2 个数值模拟软件之间的结合，该方法成功实现了对水力压裂煤层进行瓦斯抽采的数值模拟。

2）40 MPa 与20 MPa 注水压力相比，注水压裂的范围较大，产生的裂隙数量较多，裂隙扩展范围较大，压裂效果好。随着抽采时间的增加，瓦斯抽采的范围逐渐增大，抽采流量逐渐降低。且相同抽采时间下，40 MPa 注水压裂的煤层瓦斯抽采范围和抽采流量大，水力压裂可以提高煤层的瓦斯抽采效果，与实际工程应用效果相符，模拟效果较好，该模拟方法具有可行性，能够为现场应用提供指导。