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

梁卫国1,2,杨健锋1,2,廉浩杰1,2,王志勇3,沈文豪3

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

:线弹性断裂力学作为一种十分成功的断裂理论框架,已被广泛地应用于表征固体材料中裂纹扩展行为。对于线弹性岩石断裂力学来说,岩石一般被简化为脆性材料,相对于裂纹尺寸及试件尺寸,其裂纹尖端前断裂过程区(Fracture process zone,FPZ)范围很小可以被忽略。而另一方面,煤的破坏形式通常表现为韧性破坏,即其应力峰值后存在明显的应变软化区。对于这种韧性材料,其断裂过程区尺寸范围相对较大且会对材料的断裂行为产生很大的影响,因此线弹性断裂理论不再适用于描述煤体中裂纹扩展。而黏聚型模型(Cohesive zone model,CZM)被证明是一种有效的理论工具,能够描述韧性材料断裂过程区中的断裂行为。在该黏聚型本构模型理论中,裂纹尖端前的断裂过程区被简化为一条闭合的裂纹或闭合的裂纹面(分别对应二维及三维情况),其中断裂过程区内非线性断裂行为通过黏聚力与相对位移之间的本构关系进行表征。通过对煤进行圆盘形紧凑拉伸试验建立了不同煤阶煤(其中包括弱黏煤、气煤、肥煤、贫瘦煤及无烟煤)的黏聚型裂纹本构关系,试验结果表明,随着煤试件煤阶的升高,其初始刚度及峰值载荷逐渐升高,最大张开位移逐渐降低,试验峰后软化阶段载荷-CTOD曲线趋于线性变化且破坏形式逐渐趋于脆性破坏。采用Karihaloo多项式黏聚型本构方程对5种煤阶煤软化曲线进行拟合,得到煤体中黏聚型裂纹模型本构关系的一般形式。针对煤层松软的力学特性和韧性破坏特征,建立了基于黏聚型裂纹本构关系的煤岩水力压裂多场耦合方程组,包括多孔介质变形方程、孔隙渗流方程、裂隙渗流方程及Karihaloo多项式本构关系方程。并采用包含裂隙流水压自由度的黏聚型界面单元法进行数值模拟。结合大型真三轴水力压裂实验,验证所得煤岩水力压裂模型的正确性;根据数值模拟和物理实验结果,讨论了煤岩松软的力学特性及其裂纹尖端过程区对水力压裂的影响。

关键词:煤岩;水力压裂;韧性破坏;黏聚型裂纹

中图分类号:P618.11;TQ531

文献标志码:A

文章编号:0253-9993(2019)01-0263-08

移动阅读

梁卫国,杨健锋,廉浩杰,等.基于黏聚型裂纹本构关系的煤岩水力压裂韧性破坏模型[J].煤炭学报,2019,44(1):263-270.doi:10.13225/j.cnki.jccs.2019.0040

LIANG Weiguo,YANG Jianfeng,LIAN Haojie,et al.Modelling ductile failure of coals in hydraulic fracturing based on the constitutive equations of cohesive cracks[J].Journal of China Coal Society,2019,44(1):263-270.doi:10.13225/j.cnki.jccs.2019.0040

收稿日期:20181204

修回日期:20190108

责任编辑:郭晓炜

基金项目:国家自然科学基金面上资助项目(51874206);国家杰出青年科学基金资助项目(51225404);“三晋学者”支持计划资助项目(2013)

作者简介:梁卫国(1972—),男,山西盂县人,教授,博士生导师。E-mail:liangweiguo@tyut.edu.cn

Modelling ductile failure of coals in hydraulic fracturing based on the constitutive equations of cohesive cracks

LIANG Weiguo1,2,YANG Jianfeng1,2,LIAN Haojie1,2,WANG Zhiyong3,SHEN Wenhao3

(1.College of Mining Engineering,Taiyuan University of Technology,Taiyuan 030024,China; 2.Key Laboratory of Insitu Prooperty-improving Under Mining of Ministry of Education,Taiyuan University of Technology,Taiyuan 030024,China; 3.College of Mechanics,Taiyuan University of Technology,Taiyuan 030024,China)

Abstract:Linear elastic fracture mechanics (LEFM) has become an enormously successful theory framework in characterizing the crack propagation in a wide range of solid materials.For linear elastic rock fracture mechanics,rocks is generally simplified to brittle materials,in which the fracture process zone (FPZ),i.e.the region ahead of the crack tip where micro-cracks initiate and coalesce,is small enough to be ignored compared to the size of the crack and the size of the specimen.On the other hand,coals usually exhibit quasi-brittle failure behaviors,char-acterized by a strain softening regime after the peak stress.For the ductile material,the size of FPZ is considerably large and has strong impact on the fracture behavior,and thus the theory of LEFM is no longer suitable for characterizing the crack propagation in coals.The cohesive zone model (CZM) has proven a useful theoretical tool to describe the fracture behavior in the FPZ of ductile materials.In the theory of CZM,the FPZ in front of the real crack tip is lumped hypothetically into a discrete line or plane corresponding to two-dimensional or three-dimensional cases,respectively,and the nonlinear behavior occurrence in the FPZ is represented by a con-stitutive equation that relates the cohesive stresses to the displacement jump across this line or plane.The consti-tutive relationships of CZMs for the different rank coals,including weakly caking coals,gas coals,fat coals,meager-lean coals and anthracite,have been determined by the disk-shaped compact tension (DC(T)) tests.The results show that the initial stiffness and peak loads increase with the coal rank increase.As the coal rank increases,the critical crack separation displacement reduces,and meanwhile the shape of the post-peak softening curves tends to be linear and the failure mode becomes more brittle.To arrive at a general form of constitutive relationships of CZMs in coals,the Karihaloo’s polynomial cohesion-separation law is applied to fit the softening curves of the five different rank coals.By considering the mechanical properties and the ductile failure of coals,the authors establish the multi-physics coupling equations for the hydraulic fracturing in coals,including the deformation equation of porous media,fluid flow equation in the pores,fluid flow equation in the fractures and Karihaloo’s polynomial constitutive equation,and simulate it with cohesive interface elements.Combined with a large-scale hydraulic fracturing experiment,the obtained mathematical model of hydraulic fracturing in coals is verified,and the effect of coal properties and the fracture process zone on hydraulic fracturing is discussed.

Key words:coals;hydraulic fracturing;ductile failure;cohesive crack

煤层气是指赋存于煤层中以甲烷(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)

其中,σij为多孔介质总应力;为有效应力,此处假定受到拉应力状态为正;p为孔隙水压力;δij为kronecker-delta符号;α为比奥特系数,其计算公式为

(3)

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

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

(4)

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

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

(5)

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

1.2 孔隙渗流方程

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

(6)

式中,φ为孔隙率; ρsρw分别为固体和流体的密度,因此多孔介质的密度是为固体骨架的速度;为流体和固体间的相对速度。

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

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

式中,ρs0为多孔介质固体骨架初始密度;为有效应力;为初始有效应力;ρ0为流体初始密度;cf为流体滤失系数。

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

p+ρwg)(8)

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

1.3 裂隙流方程

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

(9)

式中,q为裂隙流的流量;为裂纹的开度;Q(t)为注射点流量;qtqb分别为压裂液从裂纹上下表面的滤失量;δ(x,y)为Dirac-delta函数;裂隙流的流量和压力梯度满足立方定律:

(10)

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

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 实验用煤样工业分析及产地
Table 1 Proximate analyses and the deposit location of these coals used for testing

煤试件类型弱黏煤气煤肥煤贫瘦煤无烟煤含水率/%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 煤岩黏聚型裂纹本构关系的建立

将以上实验过程中得到的荷载与裂纹张开位移分别除以其最大值,即归一化后,可得到不同煤阶煤的峰后软化曲线,其中,P为拉伸载荷,Pcr为载荷峰值;其中,w为裂纹张开位移,wc为最大裂纹张开位移。图4中给出了5种不同煤阶煤的曲线。

图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

为了将不同煤阶煤的黏聚型裂纹本构关系统一于同一个数学表达式,采用Karihaoo多项式本构方程对上述不同煤阶煤的曲线进行拟合,Karihaloo多项式本构方程为KARIHALOO 与XIAO[20-21]通过黏聚型裂纹应力场与位移场关系推导得到的解析本构方程,即

(12)

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

表2 Karihaloo多项式本构方程拟合不同煤阶煤的结果与不同煤阶煤断裂能
Table 2 Fitting results of the different rank coals by Karihaloo polynomial cohesion-separation laws and the fracture energy of the different rank coals

煤试件类型弱黏煤气煤肥煤贫瘦煤无烟煤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)中a1a2a3a4a5五个拟合参数。5种煤阶煤DC(T)实验测定的曲线进一步反映出煤岩的韧性破坏特征,并表明不同煤阶煤的软化曲线有别于常见的线性模型、双线性模型、指数曲线模型等。

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

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

参考文献

[1] ADACHI J,SIEBRITS E,PEIRCE A,et al.Computer simulation of hydraulic fractures[J].International Journal of Rock Mechanics and Mining Sciences,2007,44(5):739-757.

[2] NORDGREN R P.Propagation of a vertical hydraulic fracture[J].Society of Petroleum Engineers Journal,1972,12(4):306-314.

[3] PERKINS T K,KERN L R.Widths of hydraulic fractures[J].Journal of Petroleum Engineering,1961,13(9):937-949.

[4] GEERTSMA J,DE KLERK F.A rapid method of predicting width and extent of hydraulically induced fractures[J].Journal of Petroleum Technology,1969,21(12):1571-1581.

[5] KHRISTIANOVIC S A,ZHELTOV Y P.Formation of vertical fractures by means of highly viscous liquid[A].Proceedings of the fourth world petroleum congress[C].Rome,1955:579-86.

[6] ZHANG X,DETOURNAY E,JEFFREY R.Propagation of a penny-shaped hydraulic fracture parallel to the free-surface of an elastic half-space[J].International Journal of Fracture,2002,115(2):125-158.

[7] SETTARI A,CLEARY M P.Development and testing of a pseudo-three-dimensional model of hydraulic fracture geometry[J].SPE Production Engineering,1986,1(6):449-466.

[8] CARTER B J,DESROCHES J,INGRAFFEA A R,et al.Simulating fully 3D hydraulic fracturing[J].Modeling in Geomechanics,2000,200:525-557.

[9] BOONE T J,INGRAFFEA A R.A numerical procedure for simulation of hydraulically-driven fracture propagation in poroelastic media[J].International Journal for Numerical and Analytical Methods in Geomechanics,1990,14(1):27-47.

[10] 杨健锋,梁卫国,陈跃都,等.不同水损伤程度下泥岩断裂力学特性试验研究[J].岩石力学与工程学报,2017,36(10):95-104.

YANG Jianfeng,LIANG Weiguo,CHEN Yuedu,et al.Experiment research on the fracturing characteristics of mudstone with different degrees of water damage[J].Chinese Journal of Rock Mechanics and Engineering,2017,36(10):95-104.

[11] YANG J,LIAN H,LIANG W,et al.Experimental investigation of the effects of supercritical carbon dioxide on fracture toughness of bituminous coals[J].International Journal of Rock Mechanics and Mining Sciences,2018,107:233-242.

[12] HILLERBORG A,M MODER,PETERSSON P E.Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[J].Cement & Concrete Research,1976,6(6):773-781.

[13] XU X,NEEDLEMAN A.Numerical simulation of fast crack growth in brittle solids[J].Journal of the Mechanics and Physics of Solids,1994,42(9):1397-1434.

[14] CARRIER B,GRANET S.Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model[J].Engineering Fracture Mechanics,2012,79(Pt8):312-328.

[15] NGUYEN V P,LIAN H,RABCZUK T,et al.Modelling hydraulic fractures in porous media using flow cohesive interface elements[J].Engineering Geology,2017,225:68-82.

[16] SHEN W,ZHAO Y P.Combined effect of pressure and shear stress on penny-shaped fluid-driven cracks[J].Journal of Applied Mechanics,2018,85(3):031003.

[17] SHEN W,ZHAO Y P.Quasi-static crack growth under symmetrical loads in hydraulic fracturing[J].Journal of Applied Mechanics,2017,84(8):081009.

[18] WAGONER M P.Disk-shaped compact tension test for asphalt concrete fracture[J].Experimental Mechanics,2005,45(3):270-277.

[19] SONG S H,WAGONER M P,PAULINO G H,et al.δ25 Crack opening displacement parameter in cohesive zone models:Experiments and simulations in asphalt concrete[J].Fatigue & Fracture of Engineering Materials & Structures,2010,31(10):850-856.

[20] KARIHALOO B L,XIAO Q Z.Asymptotic fields at the tip of a cohesive crack[J].International Journal of Fracture,2008,150(1-2):55-74.

[21] KARIHALOO B L,XIAO Q Z.Asymptotic fields ahead of mixed mode frictional cohesive cracks[J].ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift f1/4r Angewandte Mathematik und Mechanik,2010,90(9):710-720.