陈 鹏， 孙可明*,2， 张 宇

(1.辽宁工程技术大学 力学与工程学院，阜新 123000； 2.青岛理工大学 理学院，青岛 266000)

1 引 言

我国煤层气资源丰富，具有巨大的开发潜力，但抽采量和利用量较低，增速缓慢，主要原因是煤储层具有低压力、低渗透率以及低含气饱和度的三低储层特征[1]，渗透率一般在(0.1～0.001)×10-3μm2范围内[2]，导致煤层气抽采难度大，时间长，成本高，因此煤储层渗透率低成为商业化开采的瓶颈。超临界CO2气爆技术是一种新型的致裂增透技术，相比深孔化学炸药爆破的增透方法，超临界CO2物理气爆增透方法具有安全、高效、环保、适用范围广和可控性好等优点，特别对于高瓦斯煤层气开发其优势更为明显，超临界CO2近似液体的密度使气爆具有较强的冲击效果，使煤体内生成大量气爆裂缝，进而提高了增透效果。

天然煤体经过长期地质作用，孕育了大量的割理裂隙和层理等结构弱面，这些结构弱面的存在使煤体表现较强的非均质性，并影响着煤体的力学特性。宋浩然等[3]基于煤体的各向异性和非均质性，考虑煤体应力变形场和瓦斯渗流场的交叉耦合作用，分析了煤层抽采中水力割缝钻孔周围瓦斯压力以及渗透率的时空演化规律；宋红华[4]等基于波速测试、CT扫描和三维重构技术，分析了煤岩内部原生裂隙、孔隙以及矿物夹杂分布的非均质特征，研究了单轴受压条件下，沿不同方向加载时，煤样内部结构分布的非均质性对煤岩破坏特征的影响；贾慧敏[5]研究发现，煤岩孔隙结构具有分段分形特征，煤岩分形维数越大，其孔隙度和渗透率越小，孔隙分选性越差，非均质性越强；赵瑜等[6]采用非均质固气耦合数值试验的方法，研究不同瓦斯压力作用下，非均质煤岩抗压强度尺寸效应的影响规律及其作用机理。孙可明等[7-10]进行了大量的超临界CO2气爆煤体致裂实验研究及均质煤体的气爆模拟研究，且目前关于模拟非均质煤体爆破的相关研究鲜有报道，因此有必要进一步开展超临界CO2气爆致裂非均质煤体的模拟研究，对提高超临界CO2气爆技术的应用水平具有重要作用。

2 超临界CO2气爆煤体力学模型

2.1 爆轰波守恒方程

超临界CO2气爆瞬间产生的爆轰波满足质量、动量和能量守恒方程[11]

(1,2)

(3)

2.2 煤体变形场方程

(1) 增量型平衡微分方程

(4)

(2) 增量型几何方程

(5)

式中Δεi j为应变的增量，Δui，j为位移对坐标偏导数的增量。

(3) 煤体损伤本构模型

由于煤体气爆是大应变和高应变率的动态过程，煤体材料采用Johnson-Cook模型，其屈服应力可表示为[12]

(6)

3 超临界CO2气爆非均质煤体有限元模型

为了建立非均质含弱面煤体的几何模型，基于数字图像方法开发了对实际煤体图像识别处理的Matlab程序，将实际煤体图像灰度化，利用灰度值将400×400像素的图像分隔为400×400个节点矩阵，利用此节点矩阵和像素代表的真实长度构建出非均质煤体的几何表征模型，将其导入Abaqus中生成非均质煤体有限元模型，实际煤体图像与有限元模型对比如图1所示。

图1 实际煤体图像与非均质煤体有限元模型

为减少计算工作量和便于对比分析，建立的煤体模型尺寸为400 mm×400 mm×20 mm，煤体中心处气爆孔直径为16 mm(与实验模型[9]尺寸一致)，与爆孔轴线平行的四个边界面设置固定约束，且各边界无反射，煤体x和y方向施加的初应力如图2所示，α为弱面倾角。煤体上的监测点布置方式如图3所示，x和y方向布置的监测点到爆孔中心的距离相同，分别为10 mm，30 mm，60 mm，100 mm，140 mm及180 mm。实验得到了煤体及弱面的物理力学参数[9]，列入表1。

图2 气爆煤体计算模型

图3 监测点布置

表1 煤体及弱面的物理力学参数

模拟中采用JWL状态方程确定超临界CO2气爆过程中爆生气体的压力变化，其表达式为[15]

(7)

模拟过程采用了SPH与FEM联合求解的计算方法，煤体与超临界CO2单元类型均为C3D8R，SPH粒子直径为1 mm，煤体模型中实体单元达到转化阈值时转变为SPH粒子，转化阈值采用应力标准，转变为SPH粒子的区域即为煤体破坏区；计算开始时,超临界CO2有限单元直接转变为SPH粒子，粒子可以进入非有限单元区域与有限单元发生相互作用；粒子与粒子、粒子与有限单元之间接触类型采用通用接触来满足质量、动量和能量守恒。

4 含单个弱面煤体气爆的模拟结果及分析

地下煤层往往受地应力作用，文献[16]给出地应力公式为

σv=γH

(8)

式中σv为铅锤应力，γ为上覆岩体的平均重力密度，H为岩体单元的深度。

本文选取地下800 m深度煤层所受地应力情况进行模拟，煤层顶板岩石的平均重力密度γ为 25 kN/m3，由式(8)可得σv为20 MPa，因此模拟中将煤体x和y方向的初应力工况设置为静水应力状态σx=σy=20 MPa及非静水应力状态σx=0 MPa，σy=20 MPa和σx=20 MPa，σy=10 MPa。

为了控制变量的统一，将一个含单个弱面煤体的图像经过图像的拉伸和旋转得到图4和图8所示的图像，保证了这些图像中煤体的弱面为相同弱面，弱面厚度约为10 mm。图4弱面倾角为0°，弱面到爆孔中心距离分别为40 mm,80 mm,120 mm和160 mm；图8弱面到爆孔中心距离为80 mm，弱面倾角分别为0°，15°，30°和45°。

图4 实际煤体图像

4.1 弱面到爆孔中心距离对气爆致裂效果的影响

图5(a)气爆产生的主裂缝沿x和y方向扩展，可以看出弱面对裂缝扩展具有阻碍作用，随着弱面到爆孔中心距离的增大，弱面对裂缝扩展的阻碍作用逐渐增大，穿过弱面扩展的裂缝密度逐渐减小；图5(b)气爆产生的主裂缝只沿y方向扩展，随着弱面到爆孔中心距离的增大，穿过弱面扩展的裂缝尺度逐渐减小，弱面到爆孔中心距离为160 mm时，没有裂缝穿过弱面扩展；图5(c)气爆产生的主裂缝只沿x方向扩展，主裂缝上生成的分支裂缝偏向y方向扩展。

图5 不同初应力下弱面到爆孔中心不同距离气爆结果

由模拟结果可知，弱面到爆孔中心距离较小时，有裂缝穿过弱面扩展，随着距离的增大，穿过弱面扩展的裂缝尺度和密度不断减小，当弱面距离爆孔中心足够远时，弱面完全阻断了裂缝的扩展，没有裂缝穿过弱面。

图6和图7给出了部分含单个弱面煤体气爆时裂缝扩展及应力波传播过程。在超临界CO2气爆初始时间段，爆孔周围介质在气爆应力波作用下发生粉碎性破坏；45 μs 后爆孔周围介质继续破坏，爆孔四周出现环向裂缝，裂缝随着应力波传播时间的增加不断扩展，裂缝扩展到弱面时，弱面阻碍了裂缝的扩展，仅有少量裂缝穿过弱面继续扩展；240 μs后主裂缝不再扩展，只有分支裂缝在扩展。初始时应力波呈现圆环状向外传播，75 μs时应力波传播到弱面，表现非环状特征，75 μs～120 μs 可以看出弱面对应力波的阻隔作用，应力波未穿透弱面，说明应力波在弱面上发生吸收和反射，135 μs后应力波传播到弱面上侧的煤体上，应力波发生衰减，165 μs后弱面下侧的应力波穿过计算区域，240 μs后应力波穿过整个煤体。

图6 裂缝扩展过程

结合应力分析可知，气爆初始时，应力波主要表现为压力波，对爆孔周围煤体进行挤压破坏形成粉碎区，接着爆孔周围产生环向裂缝，爆生气体进入初始裂缝驱动其继续扩展形成裂隙区；应力波传播到弱面时，应力值出现不连续，入射应力波分化为反射应力波和透射应力波，反射应力波表现为拉应力，使弱面内侧产生向爆孔方向扩展的反射裂缝，在入射压应力波和反射拉应力波的作用下，造成煤体拉压累积损伤，反射裂缝与爆孔处产生的径向裂缝交汇贯通，使弱面和爆孔之间形成裂缝密集区；透射压应力波继续在煤体内传播，由于弱面对裂缝扩展的阻碍，弱面外侧煤体上只生成少量裂缝。

图7 应力波传播过程

图8 实际煤体图像

4.2 弱面倾角对气爆致裂效果的影响

由图9(a，b)可以看出，随着弱面倾角的增大，弱面对裂缝扩展的阻碍作用逐渐增大，穿过弱面扩展的裂缝尺度和密度逐渐减小；由图9(b，c)可以看出，45°倾角弱面对裂缝扩展的阻碍作用最大，几乎没有裂缝穿过弱面。

图9 不同初应力下弱面不同倾角气爆结果

对比图5和图9的3种初应力工况下的模拟结果可得，气爆主裂缝沿最大初应力方向扩展，初应力对裂缝的萌生和扩展具有导向作用。

根据图3监测不同弱面倾角下y方向各个监测点处的振动速度峰值，监测数据列入表2。

表2 质点振动速度峰值

图10 质点振动速度峰值曲线

由图10可以看出，监测点1和监测点2与弱面距离较大，受弱面的影响较小，不同倾角时同一监测点的振动速度相差不大；弱面倾角为0°和15°时，监测点3接近弱面，倾角为30°和45°时，监测点4接近弱面，由曲线可看出，接近弱面处的监测点振动速度发生了突变，因为应力波传播到弱面时，在弱面上产生能量聚集，造成应力集中，导致弱面附近的质点振动速度大幅度增大；弱面后的质点振动速度峰值随弱面倾角的增大而减小。

应力波通过弱面后会发生衰减，弱面不同倾角时，监测点6均位于弱面之后，监测得到该点的Mises应力峰值列入表3。数据表明，弱面后Mises应力峰值随其倾角的增大而减小，由图11可知应力波通过弱面时，弱面倾角越大，应力波衰减越多。

表3 监测点6的Mises应力峰值

图11 监测点6的Mises应力峰值曲线

文献[17]给出爆炸应力波透射率为透射应力与入射应力的比值，即

T=σt/σi

(9)

式中T为应力波透射率,σi为入射应力,σt为透射应力。

按照应力波向弱面入射和透射的时间及其入射和透射方向，监测弱面上及接近弱面两端的煤体上的应力，得到应力波由煤体传入弱面及由弱面传入煤体时的入射应力和透射应力，根据式(8)计算出弱面不同倾角时的应力波透射率(表4)，并得到 图12 所示的曲线。

表4 应力波透射率

结果表明，应力波由煤体传入弱面时的透射率小于1，应力波由弱面传入煤体时的透射率大于1；随着弱面倾角的增大，应力波透射率不断减小，应力波透射率衰减度不断增大。

5 含多弱面煤体气爆的模拟结果及分析

图1给出了含多个弱面煤体的实际图像及其经过图像识别处理得到的有限元模型，给煤体x和y方向均施加20 MPa初应力，得到气爆裂缝扩展过程及应力波传播过程如图13和图14所示。

对比含单个弱面的煤体气爆结果可知，多个弱面的存在改变了煤体气爆的粉碎区和裂隙区形貌，应力波波形变得不规则；煤体中弱面越多，应力波衰减越多，弱面处生成向爆孔扩展的反射裂缝越多，爆孔与弱面之间煤体上的裂缝密度越大。

图12 应力波透射率曲线

图14 应力波传播过程

图15 x方向应力时程曲线

对比同一监测点上的速度和应力时程曲线可知，质点在应力波作用下产生振动，质点起振时刻与应力波到达时刻一致。x方向监测点7处介质在气爆冲击载荷作用下瞬间达到抗压强度发生压缩破坏，承载应力得以释放，应力变为0恒定不变；监测点8处介质受到较大冲击应力的作用并发生屈服，之后回弹，应力逐渐减小；应力波随时间的推移向外传播，传到弱面时发生反射和吸收，反射应力波为拉力波，煤体发生拉压累积破坏，应力波通过弱面后发生衰减，造成监测点处的速度和应力大幅度减小，曲线波动较小。y方向在爆孔边缘存在一个弱面，受弱面影响，监测点1的应力和速度发生突变，其后监测点的速度和应力时程曲线变化规律与x方向的规律比较接近。

图16 y方向应力时程曲线

图17 x方向速度时程曲线

图18 y方向速度时程曲线

6 结 论

本文基于数字图像方法开发了对实际煤体图像识别处理的Matlab程序，建立了非均质煤体有限元模型，模拟了不同初应力条件下超临界CO2气爆非均质煤体的致裂过程，分析了初应力及煤体的弱面对超临界CO2气爆裂隙演化和应力波传播的影响规律，得到以下主要结论。

(1) 弱面倾角相同时，弱面到爆孔中心距离距离越大，弱面对裂缝扩展的阻碍作用越大，穿过弱面的裂缝尺度和密度越小。

(2) 弱面到爆孔中心距离相同时，弱面倾角越大，弱面对裂缝扩展的阻碍作用越大(弱面倾角为45°时阻碍作用最大)，穿过弱面的裂缝尺度和密度越小，应力波透射率越小，应力波透射率衰减度越大。

(3) 煤体受初应力作用时，主裂缝会沿最大初应力方向扩展，初应力对气爆裂缝的萌生和扩展有导向作用。

(4) 天然煤体中的弱结构面对超临界CO2气爆致裂效果影响较大，弱面数量越多、形态越复杂，煤体上生成的气爆裂缝密度越大，可较好地提高煤体的致裂增透效果。