基于黏聚型裂纹本构关系的煤岩水力压裂韧性破坏模型

梁卫国，杨健锋，廉浩杰，王志勇，沈文豪

(1.太原理工大学 矿业工程学院，山西 太原 030024; 2.太原理工大学 原位改性采矿教育部重点实验室,山西 太原 030024; 3.太原理工大学 力学学院，山西 太原 030024)

煤层气是指赋存于煤层中以甲烷(CH4)为主要成分的非常规天然气，是一种重要的优质清洁能源。我国煤层气资源丰富，埋深浅于 2 000 m 的煤层气有 36.81 万亿m3(据 2006年国土资源部全国煤层气资源评价)，居世界第三位，是我国重要的能源储备。实现煤层气的高效开采和利用，能够弥补我国常规油气资源的不足，并且对保护大气环境、促进煤矿安全生产有重大意义。煤层渗透率通常较低，导致煤层气开发难度大，效率低。水力压裂是目前国内外开发非常规天然气的重要技术手段。

水力压裂属于复杂的多场耦合问题[1]，学术界对水力压裂进行了大量的研究，并发展出了一些经典的水力压裂力学模型来指导工程实践，例如 PKN模型[2-3]、KGD模型[4-5]、Penny-shape模型[6]、 拟三维模型[7]和全三维模型[8]等。在经典水力压裂模型中，岩石固体骨架被简化为不可渗透的线弹性体。对于孔隙充分发育的材料，如煤岩等，岩石内部的渗流现象明显,因而孔隙压的作用不可忽略。为考虑孔隙压的影响，BOONE等[9]在水力压裂模型中将岩石抽象为多孔介质，耦合了孔隙内部渗流场和固体骨架的变形。

但多孔介质压裂模型中仍然将岩石假设为多孔弹性体，一般采用线弹性断裂力学(Linear Elastic Fracture Mechanics,LEFM)理论对其裂纹扩展行为进行研究[10-11]，而没有考虑塑性变形和损伤的影响。在现实情况中，煤岩力学特性较软，是一种典型的准脆性材料，其断裂的重要特征是韧性破坏，即起裂之后能量并非瞬间释放，而是呈现一个损伤逐渐积累直至完全破坏的渐进式过程，因而线弹性断裂力学并不适用。解释材料渐进式破坏的有效理论工具是黏聚型裂纹模型[12-13]。黏聚型裂纹模型是针对裂尖附近断裂过程区(Fracture Process Zone,FPZ)的破坏规律建立的本构关系。断裂过程区是裂纹起裂扩展过程中主要的能量耗散区和非线性响应区。在多场作用下，断裂过程区的微观孔裂隙会不断萌生、扩张和融合，由此决定了裂纹的扩展规律和缝网形态。对于水力压裂来说，断裂过程区的损伤还会导致局部渗流特征的显著变化，形成更为复杂的流固耦合机制[14-15]。黏聚型裂纹模型的原理是把断裂过程区假设为真实裂纹前端的虚拟裂纹，虚拟裂纹的上下表面之间存在黏聚力的作用(图1)。黏聚力是裂纹上下表面相对位移的函数(黏聚型本构关系)，并且受界面损伤因子控制。在外界作用下，黏聚区域的黏聚力不断损伤，当其完全损伤时，意味着黏聚力完全消失，同时也标志着新的宏观裂纹面的产生，而这个损伤过程所消耗的能量大小就是断裂能。黏聚型裂纹本构关系通常是非线性的，并且依赖于加载、卸载路径和界面损伤变量的演化历史。常见的黏聚型裂纹本构关系有双线型，指数型，梯型等，但具体形式往往需要根据实验和理论共同确定，而且还需要考虑复合加载[16]、摩擦[17]、加载率等诸多因素可能造成的影响。由此可见，黏聚型裂纹的本构模型反映了断裂过程区的损伤破坏规律，对准脆性材料的断裂特征有关键影响，但是，目前还尚未有学者针对煤岩断裂过程区的特点对其黏聚型裂纹本构模型进行研究，因而也就不能建立起符合实际的煤岩水力压力控制方程。

基于以上分析，笔者拟通过实验建立煤岩断裂过程区的黏聚型裂纹本构关系,并在其基础上构建煤岩水力压裂多场耦合模型。模型结果的正确性将通过真三轴水力压裂实验验证，并结合数值模拟和实验数据加以讨论。此项工作有助于揭示煤岩水力压裂韧性破坏的规律，完善煤岩水力压裂的理论模型，促进煤层气资源的高效安全绿色开采。

1 煤岩水力压裂控制方程

如引言所述，煤岩水力压裂的力学机制是一组复杂的非线性耦合方程。本文将耦合方程分成固体骨架变形方程、孔隙渗流方程、裂隙流方程以及裂尖塑性区(断裂过程区)4组分别进行分析和解释。需要特别指出的是，这4组方程属于强耦合，需要同时求解。方程组最终的未知量为固体位移ui,孔隙流水压p,裂隙流水压pf。在求出这些未知量后，即可代回相应的本构方程计算应力场及损伤区。

1.1 固体骨架变形方程

多孔介质平衡方程:

σij,j+bi=0(1)

式中，σij,j为总应力的偏导数求和;bi为体力;下角标i,j为空间维度分量;下角标“,”为求偏导数运算。总应力由固体骨架的有效应力和空隙压组合而成，即

(2)

(3)

式中，K为多孔介质的体积模量;Km为排水体积模量。

固体骨架的有效应力和应变之间的关系为

(4)

式中，λ和μ为拉梅常数;εkk为正应变之和;εij为无穷小应变。

基于小应变假设，应变和位移满足如下几何方程:

(5)

式中，εij为应变张量;ui为位移向量;ui,j及uj,i为对位移求偏微导运算。

1.2 孔隙渗流方程

流体和固体骨架构成的多孔介质满足如下质量守恒方程:

(6)

此外，固体密度和流体密度满足如下的的状态方程:

ρw=ρ0[1+cf(p-p0)](7)

本模型中假设孔隙渗流满足达西定律:

式中，uw为流体速度;us为多孔介质固体骨架速度;μw为流体动力黏度系数;p为流体压力梯度;g为重力加速度;为梯度算子。

1.3 裂隙流方程

模型中裂隙流的质量守恒方程为

(9)

(10)

其中，pf为裂隙流水压;μf为压裂液黏度;s为裂纹扩展路径的曲线坐标。立方定律假设裂纹上下表面光滑，且雷诺数较低，满足层流。除此之外，由于煤岩孔裂隙较为发育，压裂液沿裂纹上下表面向岩石内部的滤失量qt和qb不可忽略。值得注意的是，在经典水力压裂模型中，裂隙流体在裂尖的压强会出现奇异性，这并不符合物理实际。在本文中，由于黏聚型裂纹的引入，不仅去除了裂尖奇异型，也消除了压强奇异性。

1.4 裂尖塑性区演化方程

建立黏聚型裂纹的关键是确定黏聚力-裂纹开度曲线:

σy=g(w)(11)

如图1所示，黏聚力σy是裂纹张开位移w的函数，随着裂纹开度增大，黏聚力逐渐减小。黏聚力的减小意味着裂纹尖端区域开始出现损伤。当裂纹张开位移达到其最大值wc时，黏聚力减小为零, 这标志着裂纹之间的黏聚力消失，也就代表着这一区域完全破坏，新的裂纹面产生。在水力压裂过程中，裂隙间的流体压力决定了裂纹塑性区的演化，从而决定了裂纹的起裂扩展。黏聚力-裂纹开度之间的函数关系是材料属性，也是裂尖塑性区的一种本构方程，需要物理实验测定。但目前黏聚型裂纹的测定主要集中于混凝土的研究工作，而煤岩的黏聚型本构模型尚未见到相关论文。

图1 黏聚型裂纹力学模型Fig.1 Mechanical model of cohesive crack

2 煤岩黏聚型裂纹本构关系测试

2.1 圆盘形紧凑拉伸(DC(T))试验

为了将第2节介绍的多孔介质水力压裂控制方程应用于煤岩水力压裂，关键一步就是实验测定1.4节中的煤岩黏聚型裂纹本构关系。本文采用 WAGONER[18]提出的圆盘形紧凑拉伸DC(T)实验方法，其样本形状如图2所示。DC(T)样本直径为D,中间预制长度为a的裂纹，裂纹最大扩展路径长度为W。通过等距离分布在预制裂纹两侧的加载孔对试样施加上下对称的拉力为P，加载孔的直径为φ。

图2 圆盘形紧凑拉伸(DC(T))试件Fig.2 Disk-shaped compact tension (DC(T)) specimen

本次研究中，将对弱黏煤、气煤、肥煤、贫瘦煤和无烟煤进行DC(T)实验，上述5种不同煤阶煤样的工业分析及取样地点详见表1。其中，样本直径D=95 mm;厚度B为35 mm;预制裂纹和扩展长度的比值a/W设置为0.25;加载孔直径设置为13.5 mm;其距离裂纹尖端和开口的距离分别设置为d=8 mm 和c=13 mm。 最终，不同煤阶煤的DC(T)测试通过在加载孔上施加相对拉伸载荷实现。实验过程中，煤DC(T)试件的裂尖开度(Crack Tip Opening Displacement，CTOD)[19]通过一对精度为0.15 μm的位移传感器(Linear Variable Differential Transformer,LVDT)测量得到,本次实验采用位移控制方式加载，加载速度设置为在0.05 mm/min。

表1实验用煤样工业分析及产地
Table1Proximateanalysesandthedepositlocationofthesecoalsusedfortesting

煤试件类型弱黏煤气煤肥煤贫瘦煤无烟煤含水率/%6.232.080.631.632.54灰度/%16.789.237.9012.3219.37挥发分/%42.3441.2530.3514.198.02固定碳/%56.2057.2861.4378.1771.33镜质体反射率/%0.850.751.121.822.53煤样产地大同忻州窑煤矿吕梁斜沟煤矿临汾长风煤矿左权阜生矿区大同河寺煤矿

通过以上的DC(T)实验，可以测得不同煤阶煤的荷载-CTOD 曲线,如图3所示，随着煤试件煤阶的升高，煤试件的初始刚度与峰值载荷大小均逐渐上升;同时，煤DC(T)试件的临界裂纹张开位移wc随着煤阶升高逐渐减小，弱黏煤裂纹最大张开位移wc为1.120 mm，而无烟煤的最大裂纹张开位移wc已减小到0.254 mm。此外，很明显可以看出荷载峰值后会出现软化现象:随着裂纹开度缓慢增加，荷载逐渐下降。这是一个明显的准脆性破坏特征曲线。且不同煤阶的煤的软化曲线形状差别很大，对于较低阶煤弱黏煤及气煤，其载荷初始下降较为缓慢;而对于煤阶较高的煤，其载荷初始下降幅度较大，随后载荷平缓地降低;该种应力软化形式对于煤体的断裂过程有重要影响。确定峰后软化曲线变化规律的数学方程式，即黏聚型裂纹本构关系，是确定煤体黏聚型裂纹模型的关键步骤。

2.2 煤岩黏聚型裂纹本构关系的建立

图3 不同煤阶煤DC(T)试件典型载荷-CTOD曲线Fig.3 Representative load-CTOD curves for different rank coal DC(T) specimens

图4 不同煤阶煤试件的曲线 curves of the different rank coal specimens

(12)

其中，L=5;ai为相应的常数系数，i=1,…,5。以上的Karihaloo 多项式是从黏聚型裂纹尖端附近的应力场及位移场渐进解推导出来的，因而具有严格的力学基础。并且由于采用多项式形式，参数非常容易拟合。将以上的 Karihaloo多项式与实验数据相对照，可以发现Karihaloo多项式可以很好地拟合实验数据:对应着弱黏煤，气煤，肥煤，贫瘦煤、无烟煤这5种不同阶煤，其拟合度R2值分别达到 0.998,0.997,0.999,0.992，0.989。表2中给出Karihaloo多项式黏聚型本构方程对于不同煤阶煤试验结果的拟合结果。

表2Karihaloo多项式本构方程拟合不同煤阶煤的结果与不同煤阶煤断裂能
Table2FittingresultsofthedifferentrankcoalsbyKarihaloopolynomialcohesion-separationlawsandthefractureenergyofthedifferentrankcoals

煤试件类型弱黏煤气煤肥煤贫瘦煤无烟煤a1-0.888-0.716-1.7352.5310.266a212.92513.5291.763-31.259-10.728a3-76.692-96.989-14.08275.32239.230a4157.437216.96039.084-68.913-67.506a5-138.613-200.474-40.17917.77453.744拟合度R20.9980.9970.9990.9920.989GF /(N·m-1)67.3857.6849.3736.1932.76

此外,软化曲线包围着的面积即为断裂能,GF。断裂能标志着生成新的裂纹表面所需要能量的大小，计算公式为

(13)

其中，P为施加的外荷载大小;W1为裂纹扩展长度;B为DC(T)试件厚度，不同煤阶煤的断裂能详见表2。根据试验结果可知，随着煤试件煤阶的提升，其断裂能不断降低。

3 煤岩水力压裂数值模拟及验证

在建立了黏聚型裂纹本构关系之后，通过考虑黏聚型裂纹与固体骨架变形、孔裂隙渗流、压裂液滤失等物理机制的相互作用，就可以建立煤岩水力压裂多场耦合方程组。黏聚型界面单元法是一种模拟黏聚型裂纹的常用数值模拟方法。该数值模拟方法中涉及的编程过程相对简便，通过该模拟方法可以对各种黏聚型裂纹本构关系进行表达。因此，本文的数值模拟采用黏聚型单元法进行数值模拟。传统的黏聚型界面单元只适用于无表面力裂纹，本文采用一种改进的黏聚型界面单元，即在传统的黏聚型界面单元中引入水压自由度，从而能够反映裂隙流水压的作用(图5)。此外，本文中所采用的黏聚型界面单元引入了第2节中实验测定的黏聚型裂纹本构关系，反映了煤岩塑性区对水力压裂的影响。

图5 孔压黏聚型界面单元Fig.5 Pore pressure interface cohesive element

第1种数值模拟算例针对的是煤体内部水压裂纹穿越煤岩界面的物理实验。如图6所示，实验采用30 cm×30 cm×30 cm的正方体，内部为15 cm×15 cm×15 cm的煤样，外部包裹混凝土作为岩石的相似材料。作为对比，在本次模拟中虽然考虑了多孔介质耦合方程，但暂不考虑裂纹尖端断裂过程区，模拟结果如图7所示。裂纹开始沿最大主应力方向扩展，之后在煤岩界面处分叉并沿着界面继续扩展。在煤岩界面的应力集中区域穿越界面进入岩层。本次模拟大致符合水力压裂的裂纹扩展规律，说明了水力压裂多孔介质耦合方程的有效性，但是由于没有考虑煤的裂尖断裂过程区，与实验仍然有一定偏差，尤其体现在压强变化曲线和裂纹开度等。

图6 大型真三轴水力压裂试验设备Fig.6 True triaxial simulation equipment for hydraulic fracturing

图7 水力压裂数值模拟(不考虑裂尖断裂过程区)Fig.7 Numerical simulation of hydraulic fracturing (the fracture process zone at the crack tip is not considered)

图8 水力压裂数值模拟(考虑裂尖断裂过程区)Fig.8 Numerical simulation of hydraulic fracturing (the fracture process zone at the crack tip is considered)

第2种数值模拟算例是基于黏聚型裂纹的多孔介质水力压裂模拟,考虑了裂纹尖端断裂过程区的影响(图8)。由图8可以看出相较于图7，裂纹的扩展路径更为曲折，并且在裂纹扩展过程中，出现了较为密集的细小分叉，这些分叉都有利于增加岩石的渗透性。但另一方面，宏观裂纹开度较小，抑制了煤层气的渗透。

4 结 论

(1)对于不同煤阶煤的DC(T)试件紧凑拉伸试验，随着煤试件煤阶的升高，其初始刚度及峰值载荷逐渐升高;且裂纹最大张开位移逐渐降低，例如临界张开位移由低阶弱粘煤的1.120 mm降低到高阶无烟煤的0.254 mm。同时黏聚型裂纹断裂能不断下降。

(2)通过DC(T)实验测定了不同煤阶煤的黏聚型裂纹本构关系，采用Karihaloo多项式本构关系对其进行表征，同时确定了5种煤阶煤的Karihaloo多项式本构方程(11)中a1，a2，a3，a4与a5五个拟合参数。5种煤阶煤DC(T)实验测定的曲线进一步反映出煤岩的韧性破坏特征，并表明不同煤阶煤的软化曲线有别于常见的线性模型、双线性模型、指数曲线模型等。

(3)由于黏聚型裂纹模型的引入，水力压裂多孔介质模型考虑了裂纹尖端塑性区的影响。不仅消除了裂纹尖端的应力奇异型，也消除了裂隙流压强在裂尖附近的奇异性，更为符合物理实际。

(4)在忽略煤岩的裂纹尖端断裂过程区时，裂纹只有在遇到界面处或者在应力集中区域，才会出现裂纹分叉。在考虑了裂纹尖端断裂过程区的影响后，裂纹的分叉更为密集，但同时主裂纹开度较小。裂纹的密集分叉有利于增加煤岩的渗透性，但裂纹开度较小又不利于开采煤层气。因此，在水力压裂过程中，有必要根据水力压裂数学模型来优化各项工程参数。