基于物质点法的深部煤层气水力割缝卸压解吸增透规律数值模拟研究

2022-10-23卢义玉夏彬伟

煤炭学报 2022年9期

周雷，彭雨，卢义玉，夏彬伟

(1.重庆大学煤炭灾害动力学与控制国家重点实验室，重庆 400030；2.重庆大学资源与安全学院，重庆 400030)

我国煤层气资源量约占全球的13.7%，其地质资源量在36万亿m左右，潜力可观，其中埋深1 000～1 500 m和1 500～2 000 m的深部煤层气资源量分别占30%和33%。开发深部煤层气对于保障我国能源安全、优化能源结构，加快建设清洁低碳、安全高效的现代能源体系尤为重要。深部煤层气储层具有渗透率低、储层压力大和地应力高等特点，导致煤层气解吸、运移困难。如何有效开发深部煤层气是目前面临的主要难题。

水力压裂是当前应用最为广泛的煤层气井人工增产改造技术。然而使用该技术开发临兴地区深部煤层气并未获得成功。其单井产量低于2 000 m/d,且无法形成稳产。通过分析，水力压裂失效的原因主要有以下2点：① 随着煤层深度增加，应力差增大，煤岩塑性增强，导致水力压裂难以在深部煤层气储层形成较大范围且复杂的体积裂缝网络，从而无法有效提高煤层的渗透率。② 水力压裂注入的高压流体使裂缝周围煤岩体应力提高，其诱发的流固耦合效应将进一步降低储层基质渗透率和煤层气解吸速率。水力割缝是一种高效的煤层增渗措施，通过水射流在煤层中形成卸压空间，能有效增加煤岩基质渗透率，促进瓦斯解吸。该技术已在井下煤层气开采中得到成功应用。基于此，重庆大学水射流研究团队提出地面定向井+水力割缝开发深部煤层气的新思路。通过定向井在煤层中切割产生多组盘状缝槽，诱导产生人工裂缝，沟通天然裂缝系统，提高煤层渗透率，同时增大煤层卸压范围，强化煤层气解吸与运移，提高煤层气产量。

水力割缝的增产效果主要受到卸压解吸增渗区域范围大小和其空间分布特征的影响。大量学者针对水力割缝的卸压增透规律开展了系统研究。ZHANG等基于动量守恒定律和摩尔-库仑准则等，建立了高压水射流割缝缝槽槽深的预测模型，对确定割缝卸压的最大影响范围具有重要意义。张永将等通过理论建模获得了高压水射流环切割缝自卸压技术改善煤层瓦斯流动机制，并分析了割缝后钻孔周边煤体应力演化规律。ZOU等通过室内实验采用不同倾角和孔槽比的煤样，研究了倾角和孔槽比对煤力学性能的弱化作用，还对开槽煤样的裂纹模式进行了识别，以揭示开槽煤样的弱化机理。由于水力割缝卸压解吸增透涉及弹塑性形变、损伤破坏等复杂力学过程，数值仿真成为研究该问题的主要手段。LIU等利用PFC研究了水力割缝后煤体的力学性质和损伤演化并建立了割缝煤样的损伤模型，该模型可以定量描述水力割缝对煤力学性质的影响。SI等使用FLAC考虑割缝的几何形状和割缝间的间距，量化这些关键参数与割缝引起的失效区体积之间的关系，为优化割缝参数提供理论依据。刘生龙等采用PFC对多割缝煤体开展了单轴压缩数值模拟试验，探究了割缝空间分布对煤层卸压增透的规律，确定了有利于煤层卸压增透的割缝最优空间分布模式。ZHAO等通过建立了非均质煤体中应力、损伤、瓦斯扩散和瓦斯流动的多场耦合数值模型，分析了均质系数、Langmuir体积应变常数和上覆应力对煤层损伤和瓦斯抽采的影响。

前述数值模拟研究主要采用有限元、有限差分等基于连续介质理论的数值方法和离散元等基于非连续介质理论的数值方法。由于煤岩强度较低、深部地应力高，割缝形成的卸压空间会使割缝周围煤岩发生大面积破坏，是一个涉及大变形、大位移、内边界界面持续变化和接触的力学过程。有限元(如Abaqus、FLAC)等基于连续介质理论的数值方法在计算大变形问题时会出现网格畸变和内嵌，导致刚度矩阵奇异，计算不收敛。针对该问题，需要对网格进行重划分和参数插值。此外，由于大变形和大位移导致边界界面变化复杂，增加了接触计算的难度。部分学者采用固定单元网格进行简化计算，忽略割缝内边界移动接触。这种处理方式在煤岩破坏面积较大时会出现灾难性结果，计算也不收敛。颗粒离散元(如PFC)可用于计算大位移和边界界面持续变化接触的问题，但由于其单元为刚体的假定，导致其难以计算大变形的问题。此外，该方法难以对应力和渗透率演化，以及瓦斯解吸进行量化分析和评价。块体离散元(如UDEC)以及有限元+离散元的混合数值方法通过力学本构描述单元的变形行为，但在计算大变形、大位移问题时存在和有限元同样的问题。因此，考虑大变形、大位移和内边界动态接触的水力割缝卸压解吸增渗规律研究鲜见报道。

物质点法是一种基于拉格朗日和欧拉混合描述的数值方法。该方法摒弃物理网格，通过携带材料信息的物质点离散材料区域，表征材料区域的运动和变形状态，适用于分析大变形和复杂接触的力学问题。笔者基于物质点法建立了适用于煤岩水力割缝卸压增透的数值模型。借助该模型研究了临兴地区深部煤层水力割缝诱导的煤岩应力、渗透率演化规律以及煤层气解吸规律。进一步分析了水力割缝缝槽尺寸和方位对煤层卸压增透的影响规律。研究结果为定向井水力割缝高效开发深部煤层气提供了理论依据。

1 基于物质点法的弹塑性损伤数值模型

1.1 控制方程

煤岩水力割缝是一个涉及弹塑性形变、损伤破坏的力学过程。笔者采用弹塑性理论描述其弹塑性形变过程。其中，平衡方程(式(1))、几何方程(式(2))和弹性本构方程(式(3))用于弹性试算，求取弹性条件下的应力、应变、速度和位移。

(1)

(2)

(3)

式中,Δ为应力张量增量，MPa；Δ为应变张量增量；Δ为体积应变增量；Δ为偏应变增量；为煤岩密度，kg/m；为容重，N/m；为速度，m/s；为时间；Δ为位移增量，m；和分别为体积模量和剪切模量，MPa。

当应力达到一定条件时，煤岩将发生塑性流变。笔者采用摩尔-库仑(式(4))和张拉塑性屈服准则(式(5))描述剪切和张拉主导的塑性流变发生的条件：

(4)

(5)

式中，，分别为剪切和张拉屈服强度，MPa；，分别最大和最小主应力，MPa；为黏聚力，MPa；=(1+sin)(1-sin),为内摩擦角，(°)；为张拉强度，MPa。

岩石材料达到屈服条件后，将产生塑性变形。塑性变形的大小和方向可由塑性流动法则确定。根据塑性位势原理可得塑性应变增量：

(6)

式中，为非负的塑性因子，其值可由一致性条件来确定；为塑性势函数，笔者采用关联流动法则，因此针对剪切和拉伸屈服的塑性势函数可由式(4)，(5)分别确定。

塑性应变可用于描述煤岩的损伤。笔者采用塑性应变的第二偏应变不变量和塑性张拉主应变来分别描述煤岩的剪切和张拉损伤：

(7)

(8)

式中，，分别为剪切和张拉塑性损伤；为塑性变形，下角1，2，3为主应力标示。

塑性损伤可进一步描述材料强度的弱化。笔者采用线性模型描述材料强度和塑性损伤之间的关系，设置塑性损伤临界值。当塑性损伤达到临界值时，材料强度不再弱化，达到残余强度。其中，黏聚力和张拉强度的残余值设置为0。

1.2 数值求解

笔者采用物质点法对控制方程进行数值求解。在物质点法中，研究对象通过物质点进行数值离散。图1简要地展示了其计算原理，包含5步计算过程。

图1 物质点法原理

(1)基于等效原理，将物质点携带的质量、动量和应力信息通过形函数以及形函数的导数映射到节点网格。

(9)

(10)

(11)

式中，为质量，kg；为时步；为动量，kg·m/s；,int()为物质点应力映射到节点的等效力，N；为形函数；为体积，m；下标和分别为背景网格节点和物质点；,为空间坐标。

(2)结合外力荷载条件计算网格节点所受到的合力：

(+1)=,load()+()-,int()-()

(12)

式中，,load()为外部施加荷载，N；为重力加速度，m/s；为阻尼系数，用于耗散系统动能。

(3)基于动量定律和牛顿第二定律计算节点的速度(式(13))和加速度(式(14)):

(13)

(14)

式中，为加速度，m/s；Δ为时间变化量，s。

(4)利用更新后的节点速度和加速度以及形函数插值更新物质点的位移(式(15))、速度(式(16))以及应变(式(17))。利用更新后的位移移动物质点到相应位置。

(15)

(16)

(17)

(5)通过应变增量和弹性本构(式(3))进行弹性试算，更新应力。进一步通过塑性屈服准则(式(4)，(5))判定是否发生屈服。如发生屈服，通过塑性流变法则(式(6))计算塑性应变，同时修正应力以满足屈服准则。

1.3 模型验证

现今还没有适用于岩土力学且成熟的物质点法商业计算软件。笔者基于物质点法的控制方程(2.1节)和计算原理(2.2节)，采用C++编程语言开发了相应的计算程序，同时利用Nvidia的CUDA开源库实现GPU并行加速。建立的数值模型和计算程序需要进行有效验证才能用于开展后续研究。如果把物质点看作有限元中的高斯点，物质点法的计算原理和动力显示有限元极为相似。其区别在于物质点法中增加了物质点迁移这一个步骤。因此，为验证模型弹塑性计算的正确性，在计算验证程序中关闭了物质点的迁移功能，使其可和传统的有限元计算进行有效对比。

用于程序验证的计算案例如图2所示。其几何形貌为长、宽各20 m的方形平板，平板外边界法线方向设置固定位移边界条件。初始应力为29.8 MPa和为35.6 MPa。其他物性参数见表1。平板中心开挖一长度为4 m、宽度为0.5 m的空槽。在卸压空间的作用下，平板应力将重新分布。笔者通过研发的物质点法程序对该力学过程进行弹塑性计算，并和FLAC软件的计算结果进行对比，对比结果如图3所示。物质点法和FLAC计算的弹塑性应力和塑性区分布基本一直，表明笔者所开发的计算程序能正确计算弹塑性力学过程。

图2 验证计算案例模型

图3 物质点法和FLAC3D弹塑性计算结果对比

表1 验证计算参数

物质点法的优势在于可以模拟大变形、大位移过程，自动识别接触。但验证算例的位移量较小，无法体现。因此，降低了材料强度参数，开启物质点迁移功能再次进行计算。其结果如图4所示。由于空槽周围产生大量塑性破坏，整个空槽边界向内收缩，直至空槽边界发生接触产生承载支撑力。同样采用FLAC对上述过程进行计算，采用2种方式：① 节点位置不随位移更新(不考虑大位移过程)；② 节点位置随位移更新，并在空槽边界上设置接触单元。当节点位置不随位移更新时，空槽边界不会发生接触产生接触力，导致整个平板都发生破坏，显然不合理。当节点位置随位移更新时，由于大变形的作用，导致单元网格畸形扭曲，刚度产生奇异性，且接触判定困难，最终出现计算错误。因此，通过上述结果表明，物质点法能更为准确地模拟空槽大变形卸压过程。

图4 大变形和大位移条件下物质点法和FLAC3D塑性区计算结果对比

2 临兴深部煤层水力割缝卸压解吸增透规律数值模拟

2.1 地质概况和数值模拟条件

依托临兴地区本溪组煤层气储层为对象开展水力割缝卸压增透研究。临兴位于鄂尔多斯盆地东缘，山西省吕梁市。本溪组煤层埋深1 500～2 000 m，属于典型的深部煤层气储层。以临兴LX26井本溪组煤层为参照，建立地质几何模型，并使用该储层的相关参数。该煤层高5 m，顶深2 000.8 m，底深2 005.8 m。基于平面应变原理，将计算模型简化为二维，如图5所示。几何模型在和方向分别长24，27 m。中心设置一个与最小主应力平行的水力割缝，缝槽长为5 m、宽为0.15 m。背景网格尺寸为0.05 m，每个背景网格含有4个物质点。因此模型共划分为26万个背景网格，含超过100万个物质点。模型四周采用固定位移边界条件。根据现场数据，初始地应力和材料物性参数见表2。

图5 水力割缝卸压模型

表2 数值计算参数

2.2 水力割缝作用下应力分布特征和卸压解吸增透量化分析

通过本文研发的物质点数值程序对2.1节中的工程情况进行了计算。由于割缝卸压解吸增透效果(卸压瓦斯解吸和渗透率增加)受孔隙体积变化的影响，因此，笔者主要分析割缝后平均应力的变化，其分布如图6(a)所示。平均应力变化整体成对称分布，在水力割缝的两端出现小范围应力集中，而垂直于割缝方向呈现大范围应力下降，应力下降区域的面积达到252 m。应力下降区明显分成两大区域，即近缝槽的应力骤降区(应力下降超过5 MPa)和远离缝槽的应力缓慢下降区。应力缓慢下降区的范围(240 m)远大于应力骤降区(12 m)。通过对比塑性区域发现(图6(b))，应力骤降区的范围和塑性区域的范围较为吻合，表明应力骤降主要由塑性大变形引起，而应力缓慢下降区则对应弹性小变形区域。

图6 割缝卸压后平均应力变化分布和塑性区分布

渗透率是表征煤层气在煤岩中流动难易程度的重要指标。通过LX26井测井数据表明，本溪组煤层渗透率在10～10m，属于典型的低渗储层。多孔介质渗透率和孔隙度密切相关。通过大量统计数据表明，渗透率和孔隙度成3次方的正比关系：

(18)

式中，为孔喉水力半径，m；为无量纲常数；为孔隙度。

煤层孔隙度变化主要由弹塑性变形引起。通过数值计算可获取其煤层弹塑性总变形。由于煤岩骨架的刚度比含孔隙煤岩大2～3个数量级，因此煤岩体积变形可看做其孔隙度的变化：

(19)

因此，结合式(18)和(19)可获得割缝卸压后渗透率变化率分布(图7(a))。在弹性区，渗透率为原始渗透率的1～5倍，该结果和实验测量值具有可比性。

图7 割缝卸压渗透率变化率分和单位体积解吸量分布

在塑性区，渗透率的变化幅度较大，大部分塑性区域的渗透率变化率介于10～20倍，部分区域(深红色)超过100倍。总体而言，整个应力下降区范围内的渗透率平均增加了0.64倍。

促进煤层气解吸是水力割缝的另一大优势。由于孔隙体积增加、气体压力下降(卸压)，部分吸附气将转化为游离气。通过质量守恒可知，卸压后孔隙内气体质量的增量等于解吸气的质量(式(20))。结合理想气体状态方程和Langmuir方程，可将式(20)改写为式(21)。因此，可通过求解式(21)得到孔隙体积变化条件下的解吸气量。

Δ=Δ

(20)

式中，为孔隙内气体质量，kg；为解吸气的质量，kg。

(21)

式中，左侧第1项为卸压后孔隙含气量，左侧第2项为初始孔隙含气量，右侧第1项为初始吸附量，右侧第2项为卸压后吸附量;为煤层气摩尔质量，kg/mol；为卸压后孔隙压力，MPa；为初始孔隙压力，MPa；为初始孔隙度；为煤层气密度，kg/m；为摩尔气体常数，J/(K·mol)；为温度，K；为 Langmuir体积，m/kg；为 Langmuir 压力，MPa。

卸压后单位体积的解吸气量分布如图7(b)所示，其分布规律与渗透率变化率一致，主要受到孔隙体积变化影响。在塑性区，单位体积解吸气量普遍超过4 m，总解吸量达到172 m。在弹性区，单位体积解吸量相对较小，但由于压降范围较大，其总解吸量达到307 m。综上所述，仅通过单一割缝就能快速解吸479 m煤层气，起到强化解吸的作用。

2.3 水力割缝几何参数对卸压解吸增透的影响规律

水力割缝的几何参数是影响煤层气解吸增透的重要因素。笔者从缝槽宽度、缝槽长度和缝槽方位角3个方面讨论其对煤层气解吸增透的影响规律。具体参数设置见表3。

表3 水力割缝几何参数

2.3.1 缝槽宽度对卸压解吸增透的影响规律

图8展示了不同缝槽宽度下应力变化的分布特征。总体而言，应力下降区域面积、塑性破坏区域面积和最大应力降随缝槽宽度的增加而增加。其次，应力下降区域的几何特征亦发生变化。塑性应力下降区在割缝方向的长度变化不大，但垂直于割缝方向的长度随缝槽宽度的增加而增加，使塑性区从扁平椭圆状向圆状变化。弹性应力下降区在缝槽宽度较小时呈现类椭圆状分布，其长轴垂直于割缝方向，且随缝槽宽度的增加，逐渐转化为圆状分布。从图9可以看出，弹塑性应力下降区面积、解吸气总量和缝槽宽度几乎呈线性相关。弹性应力下降区面积随缝槽宽度变化的增速明显大于塑性应力下降区。相反，塑性区解吸气体积的增速大于弹性应力下降区。其原因在于缝槽宽度的增加加剧了缝槽周围煤岩塑性变形，不仅增大了塑性区范围同时加剧了孔隙体积变化，而弹性应力下降区的面积虽增大，但弹性应力未有明显变化，因此，孔隙体积变化不大。

图8 不同缝槽宽度对煤层应力变化分布的影响规律

图9 不同缝槽宽度对煤层弹塑性应力下降区面积和煤层气解吸量的影响规律

2.3.2 缝槽长度对卸压解吸增透的影响规律

图10展示了不同缝槽长度下应力变化的分布特征。应力下降区域面积和塑性破坏区域面积随缝槽长度的增加而增加，和缝槽宽度的影响一致。但在几何特征的影响规律上和缝槽宽度有差异，随着缝槽长度的增加，弹性和塑性应力下降区逐渐从圆状分布转换为椭圆状分布。从图11可以看出，弹塑性应力下降区面积、解吸气总量和缝槽长度仍然呈线性相关。同样，弹性应力下降区面积随缝槽宽度变化的增速明显大于塑性应力下降区。但和缝槽宽度影响不同，塑性区解吸气量的增速小于弹性应力下降区。其原因在于缝槽长度的增加，虽然增加了塑性区面积，但是对塑性应变影响不大，而弹性区的面积显著增大，导致弹性区解吸气的增速高于塑性区。

图10 不同缝槽长度对煤层应力变化分布的影响规律

图11 不同缝槽长度对煤层弹塑性应力下降区面积和煤层气解吸量的影响规律

2.3.3 缝槽方位对卸压解吸增透的影响规律

图12为不同缝槽方位下应力变化的分布特征。应力下降区域随着割缝角度发生偏转，在0°和90°仍沿割缝对称分布，而在其他角度呈中心对称。随着缝槽和最小水平应力角度的增加，弹塑性应力下降区面积和煤层气解吸量呈下降趋势。由于缝槽体积不变，塑性区面积变化幅度较小，但其解吸气量降低明显，降速超过弹性区。其原因在于当最大水平应力和缝槽垂直时，缝槽周围煤岩应力强度比更高，产生的塑性应变更大，促进吸附气解吸(图13)。

图12 不同缝槽角度对煤层应力变化分布的影响规律

图13 不同缝槽角度对煤层弹塑性应力下降区面积和煤层气解吸量的影响规律

3 讨论

通过数值研究结果表明，单个割缝能在其长度的4～6倍内较为均匀地提高深部煤层气储层渗透率，达到均衡增透的目的。其次，近井区域由于塑性变形的影响，渗透率得到显著改善，能有效避免开采降压过程中由于井周附近有效应力显著增加导致渗透率骤减稳产困难的问题。另一方面，仅通过单一割缝就能快速解吸几百到上千方煤层气，起到了强化解吸的作用。因此，水力割缝是一种兼具增透和强化解吸双重作用的煤层气储层改造技术。

水力割缝卸压解吸和增透区域大小、几何特征与割缝参数密切相关。增加割缝宽度、割缝长度、以及让割缝方向平行于最小水平主应力有助于强化卸压解吸效果和增加增透范围。其中，割缝宽度对增透区域面积和解吸气量的影响最大，其次是割缝长度，割缝角度影响最小。适当增加割缝宽度更有利于煤层气开采。其次，提高割缝宽度和割缝长度的比值可获得偏圆状的卸压解吸增透区域，储层改造更均匀。因此，基于水力割缝卸压增透几何特征和影响规律，同时为克服水力割缝相较水力压裂影响范围小的问题，本文在定向井+水力割缝基础之上提出了平行水平井错位多割缝卸压解吸增透技术(图14)。通过多组平行于最大水平主应力的水平井+错位布置多组水力割缝，各水力割缝形成椭圆形的改造区域能够交错分布，填充卸压空白区域，最终形成更大范围的改造区域。而径向井多水力割缝(图15)只能在径向井方向形成范围有限的改造带区域，而且径向井之间存在大范围的改造空白区域。因此平行水平井错位多割缝卸压解吸增透技术能有效避免径向井水力割缝和水力压裂遗留的改造空白带，实现深部煤层气较大范围更为均衡的卸压解吸增透。