碎软煤层若直接水力压裂，将导致水力裂缝很短、煤层钻孔井壁坍塌，工程失败[1]。为了解决煤层气地面开发这一技术难题，张群等[2]提出了间接压裂煤层(即顶板水平井分段压裂)的新思路，并在淮北芦岭矿取得日产气量1×104 m3的记录。此外，OLSEN等[3]、贾慧敏等[4]、韩保山[5]将间接压裂技术分别用于美国Piceance盆地、山西郑庄矿区、湖南洪山殿矿区，利用顶板岩层坚硬、脆性性质，在保护钻井的同时实现顶板中较长的压裂裂缝，增加煤层气单井控制面积。然而，该技术往往受到碎软煤层特殊力学性质、复杂工程条件的制约，使其成功与否的关键，即水力裂缝有效穿透煤岩界面并进入煤层的临界条件与机理并不清楚。
建立适用于碎软低渗煤层的韧性破坏-渗流耦合本构关系，是预测上述临界条件、揭示煤岩界面阻碍水力裂缝扩展机理的重要手段。目前，多数学者基于线弹性断裂力学理论开展相关研究，发现水力裂缝在靠近煤、岩层界面时会产生复杂的相互作用，而低应力差、低界面强度、低弹性模量之比等因素会造成裂尖应力强度因子钝化，使得水力裂缝停止于界面[6-13]。
但是，前人研究并未考虑间接压裂煤层技术中的一项关键设计参数，即水平井与煤岩界面的间距。更重要的是，水力压裂会导致煤结构面出现韧性断裂与渗流、煤基质出现塑性损伤并加剧压裂液滤失，由此形成复杂的煤结构面与基质的韧性破坏-渗流耦合响应过程。这一过程深刻地影响水力裂缝穿透煤岩界面的机理，进而改变水力裂缝穿透煤岩界面的临界条件。这一问题仅由线弹性断裂力学、弹性损伤力学难以解决。
笔者针对这一根源问题，围绕碎软低渗煤层的韧性破坏-渗流特征，依次研究煤结构面韧性断裂-渗流规律、煤基质塑性损伤-渗流规律，以及2者之间的应力-渗流耦合关系，建立了碎软煤层的韧性破坏-渗流耦合本构方程组。在此基础上，模拟研究了间接压裂煤层过程中地应力、煤岩界面力学性能、水平井与煤岩界面间距等3个主要因素作用下水力裂缝的扩展形式，揭示了煤岩界面阻碍水力裂缝扩展的机理，得到了水力裂缝穿透煤岩界面的临界条件表达式，并在山西晋城矿区赵庄矿得到初步应用和检验。
1 煤的韧性破坏-渗流耦合本构方程
水力压裂过程中煤的韧性破坏-渗流耦合响应包括3个方面:煤结构面的韧性断裂与渗流、煤基质的塑性损伤与渗流(或滤失)，煤基质与结构面之间的应力-渗流耦合作用(图1)。
图1 水力压裂煤层过程中的流固耦合响应示意
Fig.1 Schematic diagram of fluid-solid coupling response during hydraulic fracturing of coal
1.1 煤基质的双标量型弹塑性损伤-渗流耦合本构关系
在数值模拟中，煤基质可视作多孔介质，可赋予实体单元模拟其流固耦合特性。煤基质的塑性变形和微裂隙的扩展是导致其出现韧性破坏的2个根本原因[14]。其中，微裂隙的基本破坏模式包括拉、压两类。张拉裂隙主要由拉应力引起，是引起煤基质渗透率改变、加剧压裂液滤失的主要原因。因此，建立双标量型弹塑性损伤-渗流本构模型之前需对应力进行分解。
规定拉应力为正，压应力为负，引入函数:
(1)
式中，σ为有效应力张量中的元素;〈·〉为Macauley符号。
由此，可将有效应力张量σeff分解为拉、压有效应力张量，分别记为
者的表达式为
(2)
式中，σi 和ni(i=1,2,3)分别为主应力和主方向。
设拉、压损伤变量分别为d(+)和d(-)。假设在等温绝热条件下，弹性和塑性Helmholtz 自由能不耦合[15]。那么Helmholtz 自由能表达式为
Hp(-)(d(-),κs)]
(3)
其中,
为弹性应变矢量;κs为硬化参数矢量，其内元素分别表示等效张拉塑性应变和等效剪切塑性应变，
根据实验结果，笔者假设煤在张拉变形过程中只出现张拉损伤，无塑性应变，即
和Hp(±)分别为拉、剪应力引起的弹性和塑性Helmholtz 自由能，其表达式分别为
(4)
将式(3)，(4)代入Clausius-Duhem不等式[16]，有
(5)
结合有效应力原理，并考虑式(5)中的任意性，可以得到煤基质的双标量型弹塑性损伤本构模型。表达式为
(6)
其中，
为总应力矢量;I为单位矩阵;d为四阶损伤张量；E0 为实体单元的刚度矩阵;εs 和
分别为实体单元的总应变矢量和压缩塑性应变矢量;pw 为孔隙水压力。四阶损伤张量d为
d=d(+)N(+)+d(-)(I-N(+))
(7)
式中，
为Heaviside 函数。
α为有效应力系数，损伤前后表达式为
(8)
式中,Kb为多孔介质的排水体积模量；Ks 为固体骨架的体积模量。
为了得到节点位移和孔隙水压力这两个未知变量，还需要对式(6)中的3个内变量，即
和 d(-) 的演化规律进行研究。对于
可由连续介质空间中的塑性理论求解。加载函数F[15]为
(9)
式中，qeff为有效偏应力；peff为有效球应力;
为最大有效主应力，σb0/σc0为双轴与单轴屈服应力之比;
和
为等效拉、压应力，且有
和Sc,i(i=1,2,3)分别为张拉和剪切偏应力分量;Kc控制π 平面上的屈服面形状，当Kc = 0.667 时与Mohr-Coulomb屈服准则相近(图2(a))。
图2 偏平面上的屈服面以及无结构面煤的压、拉应力-应变曲线
Fig.2 Yield surface and the compressive/tensile stress-strain curves of the intact coal
塑性势函数的表达式为
(10)
式中，当参数δ=0.1时，可使偏平面上qeff-peff曲线与l直线交点处的尖点退化为一个圆弧，有利于数值计算的收敛[17];σUTS 为煤的拉伸强度;ψ 为剪胀角。
基于式(9),(10)，结合塑性力学及其相应的数值解法，可以得到
值。而剩余的内变量，即d(+) 和 d(-)，可通过如下方式得到。
如上所述，拉和压是煤基质韧性破坏的两种基本模式。对于压缩损伤模式，非弹性应变
和
的关系式可由图2(b)所示的几何关系确定，表达式为
(11)
其中，
和εs,c可通过三轴压缩实验获得，数据处理方法可参考文献[17]。由此得到d(-) 的演化规律。
对于拉伸损伤模式，拉伸应变εs,t,d(+) 和
的关系可由图2(c)确定，表达式为
(12)
其中，εs,t和可通过无结构面半圆盘试件的紧凑拉伸实验得到[18]。由此可以得到d(+) 的演化规律。
在弹性变形阶段，煤基质的渗透系数knd主要由孔隙度npor和弹性体积应变εs,Ve控制;而在塑性损伤阶段，其渗透系数主要由张拉损伤变量d(+)和塑性体积应变εs,Vp控制。结合渗流过程的立方定律以及渗透系数的定义[19]，可得到煤基质渗透系数的表达式:
(13)
此外，煤基质及其内的流体还需分别满足平衡方程、几何方程，以及质量守恒方程与动量守恒方程，详见文献[20]。
1.2 煤结构面的韧性断裂-渗流本构关系
在数值模拟过程中，煤结构面可由黏聚力单元表示[21]。不同断裂模式下煤的断裂过程可以分为弹性变形和韧性断裂2个阶段(见1.1节)。弹性阶段各要素本构关系表达式为
σc=D0,cεc
(14)
式中，εc 为弹性应变矢量;D0,c为黏聚力单元的弹性刚度矩阵。
当黏聚力单元的应力达到以下条件时，煤结构面进入韧性断裂过程:
(15)
其中，
为法向和两个切向的(最大)应力。
为了得到煤结构面在韧性断裂阶段的力学本构关系，将适用于二维问题的Park-Paulino-Roesler(PPR)势能函数[22]扩展到三维水力压裂问题，势能函数Ψ(Sn,SS)的表达式为
(16)
其中,Sn 为结构面法向方向的位移，
和St为两个切向方向的位移，3者均为变量;sn,sS为结构面法向、切向方向的断裂位移，2者为常数，可通过断裂力学实验得到;Γn,ΓS 为断裂能常数，且有
(17)
其中，Gn和GS=Gs+Gt分别为法向、切向方向总的断裂能;β与γ 控制纯张拉型和纯剪切型断裂牵引力-位移曲线的形状，可由曲线拟合得到;参数m，l 与β，γ有关，且有
(18)
其中，χn=sn,p/sn,χS=sS,p/sS，sn,p和sS,p 为拉、剪断裂过程中的峰值位移。
通过计算Ψ(Sn,SS)对位移的一阶导数，得到混合型韧性断裂过程中应力与位移的关系式:
(19)
(20)
(21)
需要注意的是，式(19)～(21)只有在总位移
时才成立。一旦s 超过总的断裂位移sF，则黏聚力单元完全断裂，其内应力将由共节点的实体单元承担。
在数值模拟过程中，σc,n, σc,s 和 σc,t值可通过实体单元与黏聚力单元的共享节点处的节点力确定。由式(19)～(21)求解未知变量Sn, Ss和St的值。由此进一步得到3个方向的断裂能Gn,Gs,Gt，弹性能Ge与非弹性能Gin的表达式依次为
(22)
(23)
Gin=Gn+Gs+Gt-Ge
(24)
裂缝中不可压缩流体需要满足动量守恒方程[20]。此外，采用立方定律模拟裂缝中的切向流，同时考虑煤基质塑性损伤引起压裂液的滤失效应，推导得到结构面中流体的质量守恒方程:
(w3pt)+Qc
(25)
式中，w为裂缝开度;pn,cen和pn,boun分别为裂缝中部和裂缝边缘处的流体压力;ν 为动力黏度;pt为流体压力梯度;Qc 为裂缝入口处的总流量。
基于以上各式，建立了煤的韧性破坏-渗流(Ductile failure-seepage 或 DF-S)本构方程组，进而可开展数值计算(图3)。
图3 DF-S本构方程的数值计算流程
Fig.3 Numerical calculation process of DF-S constitutive equations
图4 断裂力学实验以及拉压损伤演化规律
Fig.4 Fracture mechanics experiments and the evolution law of tension/compression damage variables
2 参数识别与本构方程验证
DF-S本构方程中的材料参数通过断裂力学实验、全应力-应变过程中的渗流实验获得。本构方程组的合理性可通过对比分析数值模拟与实验结果给予评估。
2.1 材料参数识别
煤基质的渗透系数与压缩损伤参考文献[17]，结果分别如图4(d)和表1，2蓝色曲线所示。断裂参数和张拉损伤变量可通过如下方法获取:通过紧凑拉伸试验和贯穿剪切实验，得到张拉和剪切荷载下煤的荷载-位移曲线(图4(c))。
表1 完整煤岩的材料参数
Table 1 Material parameters of intact coal and rock
材料Es/GPaμψ/(°)δknd/10-9ms-1kd/(10-6m·s-1)σ^Tensions,eff/MPaσ^Compressions,eff/MPa煤基质1.610.35200.12.004.801.629.26顶板6.580.19280.10.430.223.8531.50
注:Es和μ为煤岩基质的弹性模量、泊松比。
表2 煤结构面的材料参数
Table 2 Material parameters of coal discontinuities
参量取值参量取值Enn/GPa0.98β3.56Ess/GPa1.83γ3.36Ett/GPa1.83χn0.19σ0c,n/MPa0.67χS0.29σ0c,s/MPa1.06Gn/(N·m-1)57.74σ0c,t/MPa1.06GS/(N·m-1)343.05T01m0.3sn/10-4m3.2l0.94sS/10-4m5.9Γn2.72sn,p/10-4m0.61ΓS-1 132.5sS,p/10-4m1.78
注:Enn,Ess 和 Ett 为弹性刚度矩阵D0,c主对角线上的元素，分别表示结构面法向和2个切向方向上的弹性模量。
在此基础上，根据式(12)反推得到煤基质的张拉损伤变量(图4(d)红色曲线)。对于断裂力学参数，Gn和GS为荷载-位移曲线下的面积;β与γ通过拟合峰后曲线得到;Γn,ΓS 与m，l通过结合图4(c)与式(16)～(17)计算得到(表1～2)。
2.2 本构方程验证
DF-S本构方程的合理性可通过3点弯曲断裂力学实验以及水力压裂实验来验证。
建立3点弯曲数值计算模型，模型尺寸及边界条件如图5(c)所示。采用泰森多面体模拟煤中的结构面空间网络[23]，并在此位置处嵌入零厚度的黏聚力单元，其力学属性由式(14)～(21)确定。其余部分为实体单元，其力学属性由式(6)～(12)确定。实验与模拟结果如图5所示。
如图5(b),(c)所示，尽管混合断裂模式下的荷载-位移关系是通过理论推导得到的，但是，由数值模拟所得的裂缝开裂模式以及荷载位移曲线与实验结果高度一致，这表明DF-S本构方程组中韧性断裂理论具有合理性。
图5 3点弯曲数值模拟与实验结果对比
Fig.5 Comparison of 3-point bending numerical simulation and experiment results
DF-S本构方程可由水力压裂实验进一步验证。为反映顶板中的水力裂缝穿透煤岩界面并进入煤层的过程，实验对象为煤与水泥的组合体试件(图6(b))，压裂段位于水泥块中，距离煤岩界面25 mm。实验采用太原理工大学自主研制的TCHFSM-I型三轴压裂渗流试验机，压裂液为蒸馏水，注入流量为20 mL/min，煤岩界面干燥且无黏结剂。水力压裂实验中竖向应力σv、最大水平主应力σH、最小水平主应力σh条件见表3。
图6 水力压裂试验机以及试件尺寸
Fig.6 Hydraulic fracturing test machine and specimen size
水力压裂数值计算模型根据图6(b)建立。其中，煤结构面的分布由文献[23]确定;而水泥中的结构面布置在其中部且与σh平行。在上述结构面及界面位置处嵌入零厚度的黏聚力单元，其余部分为实体单元。煤采用DF-S本构方程组，其力学参数见表1，2和图4所示;对于弹脆性水泥而言，其结构面的本构关系由式(14)～(15)和(26)确定，其实体单元的本构关系由式(6),(9),(10),(13)确定(其中d(+)=d(-)=0)。水泥的弹性模量为煤的3倍，断裂位移为煤的1/2。对于干燥界面，其摩擦因数及剪切模量可由直剪试验确定，且有fc,r=0.72，Ess=Ett = 0.032σn-0.018(单位:GPa)，为增加数值模拟的收敛性，Enn 取值为煤的 1/3。数值模型中超孔隙水压力设置为0，水泥与煤的孔隙比分别设置为0.2和0.01，饱和度设置为1。其余边界条件与水力压裂实验一致。
表3 水力压裂实验的应力条件
Table 3 Stress condition of hydraulic fracturing experiments MPa
试件编号应力σvσHσh1号9532号9533号8534号8535号753
分析模拟实验结果(图7)，可得:
(1)在干燥煤岩界面条件下，水力压裂实验与数值模拟结果均显示水力裂缝穿透煤岩界面的应力差阈值为6 MPa。
图7 水力压裂实验与模拟结果对比
Fig.7 Comparison of hydraulic fracturing experiment and simulation results
(2)水力压裂实验与数值模拟所得的压裂段处的水压力-试件曲线相似，均可划分为3个阶段，以图7(d)为例进行说明。在a—b段，水压力随着注水体积的增加而快速增加;b—c段，当水泥在b时刻破裂后水压力迅速下降，而后在弹脆性水泥中，水力裂缝在很短时间内到达界面处(c点);c—d—e阶段，若煤结构面在煤试件的中部且与水泥中的水力裂缝共面，则水压力会快速增加(1号试件)，而当煤结构面与水泥中的水力裂缝不共面时，会导致水压力缓慢增加至d点(2号试件)，而后随着煤结构面的周期性的断裂，使得水压力随时间呈现波动变化趋势。最终，随着煤或者水泥的完全断裂，水压力降为0。
需要注意的是，水力压裂实验所得的水压力-时间曲线，其起裂时刻(b点)滞后于数值模拟结果。这是因为实验所用的试件为干燥条件，而数值模拟假设试件为水饱和的。因此，在模拟过程中，一旦在压裂段注水，压裂段及其附近区域在瞬间产生由水压力引起的附加应力;而在实验过程中，只有当压裂段周围水泥逐渐饱和之后才会有较大的附加应力响应。
3 水力裂缝穿透煤岩界面的临界条件
3.1 工程概况与数值计算模型
赵庄矿位于晋城矿区东部，主采煤层为3号煤，平均厚度5 m，开采深度约450 m，煤层含气量在0.24～18.79 m3/t，其中CH4体积分数介于76.9%～99.6%，储层压力为1.97～6.25 MPa。井田北部无大的断层和陷落柱，但是发育数个次级褶皱，导致地应力有所变化，同时也可能会导致煤岩界面的黏结状态(即剪切强度)发生变化。
图8 间接压裂煤层技术的数值计算模型
Fig.8 Numerical calculation model of indirect coal fracturing technology
为研究间接压裂煤层过程中水力裂缝穿透煤岩界面的临界条件，建立了数值计算模型(图8)。考虑DF-S本构方程组的高度非线性，模型x，y，z三个方向上的尺寸控制在10,14和20 m。注水点位于xz平面的中心且与煤岩界面的间距D=0.5,1.0,1.5,2.0,2.5 m。在注水点下方0.5 m范围内设置初始断裂区以模拟间接压裂煤层工程中的射孔长度。数值模拟与工程条件一致，压裂液为水，注水流量为8 m3/min并持续100 min。σv,σH,σh分别沿着y,z和x方向施加，岩层应力分别设计为15.8/12.2/8.5，15.8/13.2/9.5，15.8/14.2/10.5，15.8/15.2/11.5，15.8/15.2/12.5 MPa，煤层应力值设置为15.9/8.9/7.9，15.9/9.9/8.9，15.9/10.9/9.9，15.9/11.9/10.9，15.9/12.9/11.9 MPa，即煤层中应力差为Δσ=σv-σh =8～4 MPa。在此条件下，改变煤岩界面的摩擦因数fc,r，直至水力裂缝穿透界面为止。由此，得到水力裂缝穿透煤岩界面的临界D，Δσ，fc,r这3个重要的参数。煤中结构面网络采用泰森多面体近似表示[23]，材料参数与图7一致。
3.2 模拟结果
不同D，Δσ，fc,r条件下的数值模拟结果共计115个，下面仅展示部分水力裂缝成功穿透煤岩界面的结果(图9)。
图9 不同D，Δσ，fc,r下水力裂缝扩展形式
Fig.9 Hydraulic fracture propagation form in different D，Δσ，fc,r conditions
由图9可知，当Δσ=4 MPa时，随着D从0.5 m增加至2.5 m，临界fc,r会从0.22快速增加至0.96;当D=0.5 m时，Δσ从4 MPa增加至8 MPa，临界fc,r会从0.22降低至0.022。在其他Δσ与D值条件下，临界fc,r也有类似的变化规律。以上规律说明:① 较大的Δσ(或较小的σh)会使水力裂缝更容易起裂与扩展，这与文献[6]所得的结果是一致的;② 较大的fc,r值会使压裂顶板的张拉荷载低损耗地通过界面传递至煤层中，有助于水力裂缝穿透煤岩界面;③ 在Δσ和fc,r较低的条件下，减小D值有助于实现水力裂缝穿透煤岩界面。
上述现象发生的机理可通过水力裂缝穿透煤岩界面时非弹性能占水力能量的比例解释。裂隙表面能(即弹性能)的增加是驱动水力裂缝扩展的根本原因[25]。而在间接压裂煤层工程中，水力裂缝穿透煤岩界面的临界条件比较严格，根本原因在于松软煤结构面的韧性断裂、煤基质的塑性损伤以及由此引发压裂液的滤失，会导致绝大部分水力能量消耗在无助于水力裂缝扩展的非弹性能上。如图10所示，当Δσ=6 MPa且D=1.0 m时，随着fc,r从0.05增加至0.58，煤岩界面附近的煤结构面的非弹性能占水力能量的比例从98%降低至84%;当D=1.0 m且fc,r=0.37时，随着Δσ从4 MPa增加至8 MPa，煤岩界面附近的煤结构面的非弹性能占水力能量的比例从86%降低至74%;而当Δσ=6 MPa且fc,r=0.37时，随着D从0.5 m增加至2.5 m，煤岩界面附近的煤结构面的非弹性能占水力能量的比例从70%增加至98%。
以上分析表明，在过小Δσ，fc,r和过大D值条件下，裂缝表面能(即弹性能)所占水力能量的比例最低会降低至2%。同时，水力压裂过程中煤基质塑性损伤会导致其渗透系数增加3～5个数量级，大大增加了水力裂缝中压裂液的滤失，导致水力能量耗散。为此，上述两方面是导致水力裂缝难以穿透煤岩界面的根本原因。
如图10所示，增加Δσ，fc,r并降低D可以将弹性能占水力裂缝的比例增加至30%，从而提高水力裂缝穿透煤岩界面的概率。为了进一步得到适用于工程实际的水力裂缝穿透煤岩界面的临界条件，采用幂函数对数值模拟得到的临界条件D，Δσ，fc,r值进行拟合，得到如图11所示的蓝色拟合曲面。
由此，得到水力裂缝穿透煤岩界面的临界条件数学表达式为
0.689-0.569D-0.536+3.515Δσ-1.162-
1.238D-0.536Δσ-1.162-fc,r>0
(26)
式中，D为注水点与煤岩界面的间距；Δσ为竖向应力与最小水平主应力之差。
需要注意的是，上式只适用于与赵庄煤层气地质与材料参数类似的区域，且限制条件为4 MPa<Δσ<8 MPa，0.5 m < D < 2.5 m。
基于上述结果，采用间接压裂煤层技术，通过优化水平井位置(图1红线)从而增加Δσ至6 MPa，同时确保煤层完整使得fc,r稳定在0.72，再将D控制在1.0 m以下，取得煤层气日产量3 000 m3以上的结果(图12)。这表明水力裂缝成功穿透煤岩界面，从而验证了上述结果的合理性。
图10 不同D，Δσ，fc,r条件下非弹性能占水力能量的比例以 及渗透系数的演化规律
Fig.10 Ratio of inelastic energy to hydraulic energy and the evolution law of permeability coefficient under different conditions of D，Δσ，fc,r
图11 水力裂缝穿透碎软煤与顶板界面的临界D，Δσ，fc,r 条件
Fig.11 Critical D，Δσ，fc,r conditions for hydraulic fractures to penetrate the interface between coal and roof
图12 间接压裂碎软低渗煤层后日产气量曲线
Fig.12 Daily CBM production curve under indirect fracturing coal technology
4 结 论
(1)基于Helmholtz 自由能、Park-Paulino-Roesler势能函数以及相应渗流理论推导得到的DF-S本构方程组，可以较好地反映碎软低渗煤层结构面的韧性断裂-渗流、煤基质的塑性损伤-渗流耦合响应，以及2者之间的应力-渗流相互作用。
(2)煤结构面的韧性断裂、煤基质的塑性损伤以及由此引发的压裂液滤失效应使得水力能量大量耗散，导致有助于水力裂缝扩展的弹性能最低只占到水力能量的2%，这是导致水力裂缝难以穿透煤岩界面的根本原因。
(3)水力裂缝穿透煤岩界面的临界fc,r与D正相关，而与Δσ负相关。D对临界fc,r的影响更大，将D控制在1 m之内有助于水力裂缝穿透煤岩界面，保证间接压裂煤层技术的成功率。
[1] 孙四清,张群,闫志铭,等.碎软低渗高突煤层井下长钻孔整体水力压裂增透工程实践[J].煤炭学报,2017,42(9):2337-2344.
SUN Siqing,ZHANG Qun,YAN Zhiming,et al.Practice of permeability enhancement through overall hydraulic fracturing of long hole in outburst-prone soft crushed coal seam with low permeability[J].Journal of China Coal Society,2017,42(9):2337-2344.
[2] 张群,葛春贵,李伟,等.碎软低渗煤层顶板水平井分段压裂煤层气高效抽采模式[J].煤炭学报,2018,43(1):150-159.
ZHANG Qun,GE Chungui,LI Wei,et al.A new model and application of coalbed methane high efficiency production from broken soft and low permeable coal seam by roof strata-in horizontal well and staged hydraulic fracture[J].Journal of China Coal Society,2018,43(1):150-159.
[3] OLSEN T N,BRATTON T R,DONALD A,et al.Application of indirect fracture for efficient stimulation of coalbed methane[A].SPE Rocky Mountain Oil & Gas Technology Symposium[C].Denver,Colorado,2007:16-18.
[4] 贾慧敏,胡秋嘉,祁空军,等.高阶煤煤层气直井低产原因分析及增产措施[J].煤田地质与勘探,2019,47(5):104-110.
JIA Huimin,HU Qiujia,QI Kongjun,et al.Reasons of low yield and stimulation measures for vertical CBM wells in high-rank coal[J].Coal Geology & Exploration,2019,47(5):104-110.
[5] 韩保山.构造煤煤层气压裂方式及机制探讨[J].煤矿安全,2019,50(7):211-214.
HAN Baoshan.Discussion on fracturing mode and mechanism of coalbed methane in tectonic coal[J].Safety in Coal Mines,2019,50(7):211-214.
[6] AIMENE Y,HAMMERQUIST C,OUENES A.Anisotropic damage mechanics for asymmetric hydraulic fracture height propagation in a layered unconventional gas reservoir[J].Journal of Natural Gas Science and Engineering,2019,67:1-13.
[7] ZHANG F S,DONTSOV E.Modeling hydraulic fracture propagation and proppant transport in a two-layer formation with stress drop[J].Engineering Fracture Mechanics,2018,199:705-720.
[8] POLUDASU S,AWOLEKE O,AHMADI M,et al.Using experimental design and response surface methodology to model induced fracture geometry in Shublik shale[J].Journal of Unconventional Oil Gas Resources,2016,15:43-55.
[9] WANG W,OLSON J 
M,et al.Interaction between cemented natural fractures and hydraulic fractures assessed by experiments and numerical simulations[J].Journal of Petroleum Science and Engineering,2018,167:506-516.
[10] LU C,LU Y X,GUO J C,et al.Stability of the formation interface under the impact of hydraulic fracture propagation in the vicinity of the formation interface[J].Petroleum Science,2020:1-18.
[11] OYEDERE M,GRAY K,MCCLURE M W.Numerical investigation of a novel hypothesis for fracture termination and crossing,with applications to lost circulation mitigation and hydraulic fracturing[J].Journal of Unconventional Oil Gas Resources,2016,15:122-133.
[12] GUO J C,LUO B,LU C,et al.Numerical investigation of hydraulic fracture propagation in a layered reservoir using the cohesive zone method[J].Engineering Fracture Mechanics,2017,186:195-207.
[13] LAN Z K,GONG B.Uncertainty analysis of key factors affecting fracture height based on box-behnken method[J].Engineering Fracture Mechanics,2020,228:106902.
[14] MENG Q B,ZHANG M W,HAN L J,et al.Effects of acoustic emission and energy evolution of rock specimens under the uniaxial cyclic loading and unloading compression[J].Rock Mechanics and Rock Engineering,2016,49(10):3873-3886.
[15] RUBIN M B.A new approach to modeling the thermomechanical,orthotropic,elastic-inelastic response of soft materials[J].Mechanics of Soft Material,2019,1(1):3.
[16] KOMAR L A,SVISTKOV A L.Thermodynamics of elastic material with relaxing heat flux[J].Mechanics of Solids,2020,55(4):584-588.
[17] 李浩,白海波,武建军,等.D-P 随机损伤本构模型及其在预防陷落柱突水中的应用[J].岩土力学,2018,39(12):574-582.
LI Hao,BAI Haibo,WU Jianjun,et al.D-P stochastic damage constitutive model and its application in preventing water inrush of karst collapsed column[J].Rock and Soil Mechanics,2018,39(12):574-582.
[18] 杨健锋.煤体黏聚裂纹本构方程研究及其在压裂工程中的应用[D].太原:太原理工大学,2019:63-75.
YANG Jianfeng.Study on the constitutive equation of the cohesive crack in coals and its application in fracturing engineering[D].Taiyuan:Taiyuan University of technology,2019:63-75.
[19] WU C,CHU J,WU S,et al.3D characterization of microbially induced carbonate precipitation in rock fracture and the resulted permeability reduction[J].Engineering Geology,2019,249:23-30.
[20] 赵阳升.多孔介质多场耦合作用及其工程响应[M].北京:科学出版社,2017:66-72.
[21] ORTIZ M,PANDOLFI A.Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis[J].International Journal for Numerical Methods in Engineering,2018,44(9):1267-1282.
[22] PARK K,PAULINO G H,ROESLER J R.Virtual internal pair-bond model for quasi-brittle materials[J].Journal of Engineering Mechanics,2008,134(10):856-866.
[23] 祝一搏.多孔介质孔隙结构特征及其输运属性分析[D].焦作:河南理工大学,2015:47-53.
ZHU Yibo.Pore structure characteristic and transport property analysis of porous medium[D].Jiaozuo:Henan University of Technology,2015:47-53.
[24] JIANG Y L,LIAN H J,NGUYEN V P,et al.Propagation behavior of hydraulic fracture across the coal-rock interface under different interfacial friction coefficients and a new prediction model[J].Journal of Natural Gas Science and Engineering,2019,68:102894.
[25] 李世愚,和泰名,尹祥础.岩石断裂力学[M].北京:科学出版社,2015:111-147.
