碎软煤层韧性破坏-渗流耦合本构关系及其间接压裂工程验证

2021-04-11梁卫国李国富白建平王建美武鹏飞

煤炭学报 2021年3期

李浩，梁卫国,2,，李国富，白建平，王建美，武鹏飞

(1.太原理工大学原位改性采矿教育部重点实验室，山西太原 030024; 2.山西能源学院，山西晋中 030600; 3.太原理工大学矿业工程学院，山西太原 030024; 4.山西蓝焰煤层气集团有限责任公司，山西晋城 048026)

碎软煤层若直接水力压裂，将导致水力裂缝很短、煤层钻孔井壁坍塌，工程失败[1]。为了解决煤层气地面开发这一技术难题，张群等[2]提出了间接压裂煤层(即顶板水平井分段压裂)的新思路，并在淮北芦岭矿取得日产气量1×104m3的记录。此外，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符号。

(2)

式中，σi和ni(i=1,2,3)分别为主应力和主方向。

设拉、压损伤变量分别为d(+)和d(-)。假设在等温绝热条件下，弹性和塑性Helmholtz 自由能不耦合[15]。那么Helmholtz 自由能表达式为

Hp(-)(d(-),κs)]

(3)

(4)

将式(3)，(4)代入Clausius-Duhem不等式[16]，有

(5)

(6)

d=d(+)N(+)+d(-)(I-N(+))

(7)

α为有效应力系数，损伤前后表达式为

(8)

式中,Kb为多孔介质的排水体积模量；Ks为固体骨架的体积模量。

(9)

图2 偏平面上的屈服面以及无结构面煤的压、拉应力-应变曲线Fig.2 Yield surface and the compressive/tensile stress-strain curves of the intact coal

塑性势函数的表达式为

(10)

式中，当参数δ=0.1时，可使偏平面上qeff-peff曲线与l直线交点处的尖点退化为一个圆弧，有利于数值计算的收敛[17];σUTS为煤的拉伸强度;ψ为剪胀角。

(11)

(12)

在弹性变形阶段，煤基质的渗透系数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)

(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)

在数值模拟过程中，σ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]。此外，采用立方定律模拟裂缝中的切向流，同时考虑煤基质塑性损伤引起压裂液的滤失效应，推导得到结构面中流体的质量守恒方程:

(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

表2 煤结构面的材料参数Table 2 Material parameters of coal discontinuities

在此基础上，根据式(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

分析模拟实验结果(图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

基于上述结果，采用间接压裂煤层技术，通过优化水平井位置(图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