实验尺度下的煤岩水力压裂数值模拟研究*

崔俊飞

（1.瓦斯灾害监控与应急技术国家重点实验室；2.中煤科工集团重庆研究院有限公司）

随着煤矿资源开采技术的不断进步，煤层开采深度越来越深。对煤矿资源进行合理利用，进一步提高煤矿资源及煤矿内伴生资源的利用效率受到越来越多的关注［1-4］。煤矿内伴生的煤层气等资源由于地下的高压作用，又容易产生冲击地压等次生灾害。水力压裂技术作为非常规油气藏的重要开采技术，能够有效增加储层渗透率。该技术在煤层增渗和瓦斯泄压方面应用的可能性，使其在煤矿资源开采和煤矿安全防治等方面受到越来越多的关注。

近年来，既有学者提出引进页岩气资源开采的水力压裂技术增渗煤层，从而开采煤层气资源［3-4］，也有学者提出在煤层巷道内通过压裂贯通深部煤层气，降低煤层气压，从而降低煤矿灾害的发生［5-6］。水力压裂技术最早主要被用于低渗透、超低渗透页岩油气资源的开采中，压裂可形成多尺度裂缝网络，增加流体流通通道的连通性以及可流通区域体积，提高储层渗透率，进而提高开采效率［7］。对煤层水力压裂的研究方法主要分为3类：室内实验研究、数值模拟研究、工程试验研究［8］。其中，工程试验中容易受到监测技术、仪器精度等综合影响，难以掌握精确而全面的信息研究水力裂缝扩展机理。因此，对煤矿中水力裂缝扩展机理的研究，多通过室内实验和数值模拟研究进行讨论。Deng等［1-2］学者通过重庆大学自主研制的多功能真三轴试验系统对压裂进行了实验研究，讨论了煤岩水力压裂与页岩水力压裂的差异。黄赛鹏等［9］学者利用ANSYS有限元软件建立了煤层水力压裂模型，计算了煤岩破裂压力的影响因素，其模型中假设储层为全部均质储层，模拟了煤层气井压裂过程中降低煤层破裂压力的各种途径和方法，研究建立了套管射孔完井的应力计算模型，并计算了模拟煤岩破裂压力，分析了射孔长度、孔径、煤岩的泊松比与杨氏模量、地层深度等对煤岩破裂压力的影响。刘大锰等［10］学者从煤层埋深时所处的地应力可能存在差异的角度考虑，研究了煤储层地应力状态对煤层压裂增渗规律和裂缝扩展机理的影响。事实上，水力压裂已经被大量运用于低渗透页岩气开采中［11-12］，而在煤层中多处于实验和少量现场试验阶段。

考虑到水力压裂技术在页岩中已有大量现场运用，而在煤储层开发中大致处于论证阶段。本研究根据煤岩区别于页岩的天然裂缝网络特征，采用ABAQUS软件中的Cohesive单元，建立二维天然裂缝网络模型，通过该模型研究具备天然裂缝网络特征的煤岩试件的水力裂缝形成过程。同时，讨论初始破裂单元和流体注入速度对煤岩中水力裂缝扩展的影响。

1 有限元模型

1.1 模型假设

在确立几何模型和数学物理模型时，需要几何模型利于网格剖分，煤岩本构模型和水流流动的物理方程符合客观规律，并能保证数值计算的收敛性，做出如下假定。

（1）煤岩试件由煤岩颗粒组成的小尺度煤岩块体和裂隙组成，由于煤岩岩石具有明显的裂缝特征，假设完整的煤岩试件由可能发生断裂弱胶结裂隙面构成，该裂隙面称为煤岩岩石的天然裂隙面。

（2）小尺度煤岩块体由微观颗粒黏结在一起，为各向同性孔隙材料，具备低渗特性，其内的流体流动为达西流。

（3）假设所有天然裂隙面厚度为零，未破坏前为均质线弹性，弹性参数与试件内部岩块的弹性参数一致。

（4）压裂过程中，煤岩试件完全饱和，孔隙中的流体为不可压缩流体。

1.2 有限元模型

本研究选取了相关文献［1-2］中的一个煤岩表面裂缝图像，对压裂后的煤岩表面裂隙的素描图进行图像处理、裂缝数据的提取、几何模型的导入、以及全局cohesive单元的嵌入。经过此过程即形成具有一定网络特征的裂缝网络系统，有限元模型的处理过程如图1所示。

图1中，扫描图像的黑色线条部分为通过压裂后的岩石表面裂隙图像重构的煤岩天然裂隙。该有限元模型尺寸为0.1 m×0.1 m。分别固定图1中水平线上两边的纵向位移和纵向线两边的横向位移。

2 数学物理模型

2.1 流固耦合控制方程

水力压裂扩展过程是压裂液流体流动与煤岩变形的动态耦合过程，在数值模拟中可通过耦合求解应力平衡方程和流体连续性方程实现［13］。应力平衡方程为

式中，σ为有效应力矩阵，Pa；pw为孔隙压力，Pa；t为表面力矩阵，N/m2；δV为虚速度矩阵，m/s；f为体力矩阵，N/m3。

流体连续性方程为

式中，ρw为流体密度，kg/m3；nw为孔隙比，无因次；n为控制体方向向量；vw为流体渗流速度，m/s。

流体在岩石中的流动服从Darcy定律：

式中，k为渗透率矩阵，m/s；g为重力加速度向量，m/s2。

2.2 裂缝界面模型

裂隙界面通过零厚度cohesive单元模拟，可用于模拟裂缝的起裂、扩展等研究［14］。单元损伤之前，满足线弹性关系［15］

式中，σcoh为cohesive单元承受的应力矢量；σcoh-h、σcoh-s、σcoh-t分别为cohesive单元法向（垂直于cohesive单元上下表面的方向）、第一切向和第二切向（在二维情况不存在）承受的应力矢量；Kcoh为cohesive单元的刚度矩阵；εcoh为cohesive单元产生的应变矩阵；εcoh-n、εcoh-s、εcoh-t分别为cohesive单元法向、第一切向、第二切向产生的应变。

cohesive单元的损伤启始准则采用最大应力准则［16］：

式中，σn为法向应力，MPa；σs、σt为切向应力，MPa；为岩石的抗拉损伤的阈值应力，MPa；、为切向损伤的阈值应力，MPa；符号＜＞表示cohesive单元只抗拉不抗压。

对于切向流动，假设流体为不可压缩的牛顿流体，则有

式中，q为黏聚单元中切向上单位长度的体积流量向量；t为黏聚单元张开的厚度，m；μ为黏聚单元中压裂液黏性系数；p为黏聚单元中的流体压力，MPa。

压裂液在裂缝上下表面的滤失速率由下式计算［17］。

式中，qt和qb分别为流体流出cohesive单元上表面、下表面的体积流率，二维为面流率；ct和cb分别为上表面、下表面的滤失系数，m/min1/2；pt和pb为上表面、下表面的孔隙压力，MPa；pi为cohesive单元中面的流体压力，MPa。

3 数值模拟结果与讨论

根据文献中的参数取值范围［18-21］，煤岩试件数值模拟模型中主要参数的取值范围如表1所示。表1中包含了煤岩的岩石力学特征参数以及煤岩压裂模拟过程中使用的边界条件参数，其中，水平最小地应力为8 MPa，水平最大地应力为18 MPa，垂直地应力为28 MPa。

?

3.1 煤岩水力裂缝形成过程

通过将参考的取值赋予数值模拟模型后，将初始破裂单元个数设置为2个进行数值模拟计算，获得的流体压力曲线和入射点裂缝最大宽度变化曲线，如图2所示。

图2（a）表示流体压力曲线的演化过程，从图2（a）中可以看到：对于实验尺度的煤岩水力压裂而言，压裂过程可以分为4个阶段：①流体压力快速上升阶段，随着压裂流体开始注入煤岩试件，注入流体由于缺乏天然或人工的流体流动通道，流体压力快速上升，逐渐达到煤岩的起裂压力；②流体压力缓慢下降阶段，在流体压力达到煤岩水力压裂的起裂压力后，煤岩内部裂缝迅速扩展。此时，流体持续稳定地注入不会导致流体压力的进一步上升，显然是由于裂缝扩展形成了流体的泄压通道所致；③流体压力准稳态阶段，流体继续注入时，随着裂缝通道面积扩张，流动通道增大，流体与岩石接触面增大。④流体压力快速下降阶段，当流体压力准静态阶段持续一段时间后，水力裂缝逐渐触及甚至贯穿煤岩试件边界，导致流体快速倾泻，流体压力迅速降低至0 MPa附近。

图2（b）表示入射点处裂缝最大宽度的演化过程，从图2（b）可以看到压裂过程中入射点处裂缝最大宽度的演化过程主要分为3个阶段：裂缝最大宽度迅速扩张阶段、裂缝最大宽度缓慢扩展阶段、以及裂缝最大宽度迅速下降阶段。对比图2（a）和图2（b）：

①裂缝宽度迅速扩张阶段刚好对应于流体压力快速上升阶段，说明流体注入初期入射点附近可能存在微裂缝的发育和扩展；②最大裂缝宽度演化的第二阶段对应了流体压力的缓慢下降阶段和流体压力准稳态阶段，说明即使流体压力波动不大时，煤岩中的水力裂缝依然在发育。③与流体压力迅速下降阶段类似的，当裂缝逐渐触及甚至贯穿煤岩试件边界后，流体迅速倾泻，入射点处裂缝最大宽度迅速下降。同时，从图2（b）中可以看到，裂缝最大宽度并没有下降到0 mm附近，说明流体注入速度与流体滤失形成了一个准平衡。这说明由于水力压裂过程受到流体注入速度、地应力、煤岩内部天然裂缝形态、流体滤失、煤岩力学性质等综合作用，即使水力裂缝贯穿煤岩，流体从贯穿口倾泻后，入射点裂缝最大宽度可能不会下降至0 mm。特别地，从以上分析可以看出：实验尺度和现场尺度的水力压裂研究很容易受到边界范围的影响，实验尺度可看作有限边界，而现场压裂实质是无限边界，因此，亟需考虑实验尺度的煤岩水力压裂研究应该如何降低边界效应的影响。

压裂过程中，分别在压裂开始阶段t=10.9 s时、水力裂缝即将贯穿煤岩试件前t=37.9 s时、水力裂缝贯穿煤岩试件后t=48.43 s时的煤岩孔隙压力云图、位移云图、裂缝张开程度云图（PFOPEN），如图3所示。

图3中，t=10.9 s对应于入射点裂缝宽度稳定上升阶段和流体压力缓慢下降阶段后的一个时间点；t=37.9s对应于水力裂缝贯穿煤岩试件前的一个时间点；t=48.43 s对应于水力裂缝贯穿煤岩时间后的一个时间点。从图3中可以看到：①煤岩水力压裂过程中，在流体注入初期，由于水力裂缝尚未充分发育，压裂流体未能通过水力裂缝面大量流入煤岩内部，其孔隙压力云图未观察到明显的分区（图3（a）、（d））。随着流体持续注入，在水力裂缝贯穿煤岩试件后，孔隙压力云图出现孔隙压力的分区（图3（g））。说明在煤岩水力压裂过程中，随着水力裂缝的充分发育，注入的压裂流体不仅仅用于水力裂缝的形成，还导致了地层孔隙压力的升高和分区。同时，从孔隙压力云图中还可以看到，在水力裂缝贯穿岩石后，压裂过程中流入煤岩内部的压裂流体也会逐渐沿水力裂缝发生倾泻。②从位移云图3（b）、（e）、（h）可以看到，整个压裂过程中，由于煤岩本身的非均质特征，位移变换是不对称的，也是非均匀的。由于煤岩试件的天然裂缝特征、非均质特征等存在差异，因此，针对于实验尺度的煤岩水力压裂的实验和数值模拟研究中，建议尽可能多地对煤岩试件的非均质性进行区分和表征，从而更准确地分析煤岩中水力压裂扩展机理的影响因素。③从裂缝张开度云图图3（c）、（f）、（i）可以看到，煤岩试件中的水力裂缝宽度最为发育的时间出现在压裂流体贯穿煤岩试件前的一个时间点。

3.2 初始破裂单元对裂缝扩展的影响

水力压裂过程的影响因素很多，而诸如最大水平地应力、最小水平地应力等参数的研究在很多文献中已有，本研究选择初始破裂单元的个数作为裂缝扩展的影响因素进行分析。假设初始破裂单元的个数分别为2个、4个、6个、8个、10个，其流体压力曲线如图4所示。

图4描述了在煤岩水力压裂数值模拟模型建立过程中，其初始破裂单元设置对于流体压裂曲线和入射点裂缝最大宽度曲线的影响。图4（a）表示初始破裂单元设置对于流体压力曲线的影响，图4（b）表示初始破裂单元设置对于入射点裂缝最大宽度的影响。从图4（a）可以看到，初始破裂单元个数越多的流体压力曲线倾向于位于初始破裂单元个数少的流体压力曲线的下方。说明在数值模拟研究中，随着初始破裂单元个数的增加，煤岩水力压裂的流体压力曲线上的起裂压力逐渐降低，甚至可能低于扩展压力。一般而言，数值模拟中的初始破裂单元设置主要是为了求解的收敛性设置，通过图4（a），认为初始破裂单元设置可为煤岩的天然张开裂缝的属性特征提供支撑。从图4（b）可以看到，随着初始破裂单元个数增大，裂缝最大宽度增大，可能是由于初始破裂单元描绘了已开裂的煤岩天然裂缝。初始破裂单元个数越多，压裂过程中水力裂缝扩展所需的能量减少，从而促使更多的水力做功用于裂缝通道的扩宽所致。该分析可能为煤岩水力压裂区别于页岩水力压裂的模型构建提供参考。

3.3 注入速度对裂缝扩展的影响

将数值模型中的流体注入速度进行修改后，模拟结果如图5所示。图5中A、B、C、D、E5条曲线分别表示流体注入速度为1×10-7，8×10-8，6×10-8，4×10-8，2×10-8m/s时的曲线。

图5描述了压裂过程中，注入流体速度对入射口压力曲线和入射点裂缝最大宽度曲线的影响。从图5（a）可以看到，随着流体注入速度的降低，流体压力曲线的破裂压力值逐渐降低，说明流体注入速度对流体压力影响显著，且流体注入速度越大，煤岩水力压裂的起裂压力和破裂压力越大。这可能是由于流体速度不同时，流体压裂煤岩的过程中，裂缝扩展速度不同，裂缝形成的形态特征也有差异导致的。从图5（b）可以看到，大致上流体注入速度越大的裂缝最大宽度曲线位于流体注入速度越小的裂缝最大宽度曲线的上方，但这种趋势存在一定的波动性。说明流体注入速度的改变对水力裂缝的形态特征和发育状况的影响应该是一个综合性的影响，而不仅仅局限于裂缝最大宽度。因此，讨论流体注入速度对煤岩水力裂缝扩展的影响时，应该尽可能获取更多的数据，以判断其对水力裂缝扩展的综合性影响。

4 结 论

（1）煤岩水力压裂数值模拟过程中，流体压力曲线分为流体压力快速上升阶段、流体压力缓慢下降阶段、流体压力准稳态阶段、以及水力裂缝贯穿煤岩试件形成的流体压力快速下降阶段。入射点处裂缝最大宽度的演化过程主要分为3个阶段：裂缝最大宽度迅速扩张阶段、裂缝最大宽度缓慢扩展阶段以及裂缝最大宽度迅速下降阶段。煤岩水力压裂过程中，裂缝发育在流体压力曲线的前3个阶段均存在。

（2）初始破裂单元个数的增加，煤岩水力压裂的流体压力曲线上的起裂压力逐渐降低，甚至可能低于扩展压力。一般而言，随着初始破裂单元个数增大，裂缝最大宽度增大。

（3）流体注入速度越大，煤岩水力压裂的起裂压力和破裂压力越大。通过小尺度的实验或模拟进行煤岩水力压裂研究时，其尺度效应下的边界效应应该被考虑。