王永亮 ，刘娜娜 ，王 昊

（1.中国矿业大学(北京) 力学与建筑工程学院, 北京 100083；2.中国矿业大学(北京) 煤炭资源与安全开采国家重点实验室, 北京 100083）

0 引 言

水力压裂是深部致密岩体油气储层开发的核心技术，水平井多射孔簇分段压裂方法涉及储层和孔隙-裂隙内流体之间的热扩散、流体流动与岩体基质变形，热扩散效应和多物理场耦合作用是深部致密岩体压裂的典型特征。水力压裂过程中，缝网三维扩展之间的应力阴影效应和多裂缝扰动偏转行为成为影响空间缝网扩展形态和压裂效果的重要因素[1-4]。三维裂缝在水压裂缝扩展过程中伴随着裂缝的空间偏转和裂缝之间的相互挤压作用，导致裂缝的非稳定扩展[5-6]。研究发现分段压裂的射孔簇间距和起裂顺序会造成不同程度的裂缝偏转，裂缝的非稳定扩展影响裂缝的定向控制和缝网设计[7-8]。压裂裂缝扩展引起的地层应力干扰区域在多裂缝中存在叠加、覆盖行为，形成的应力阴影区将扰动裂缝扩展[9-11]。压裂裂缝三维扩展的热扩散效应影响、不同射孔簇间距以及典型压裂方案（顺序、同步、交替压裂）下裂缝网络的扩展扰动行为成为目前理论和工程实际研究的重点。

数值计算成为研究水力压裂缝网三维扩展扰动中的多物理场耦合、多裂缝非稳定扩展、工程尺度分析的重要手段。采用岩体压裂裂缝扩展的数值模型和复杂工况分析方案，压裂缝网相互作用和应力阴影效应得到了定量化分析，并对压裂裂缝起裂、扩展扰动和偏转的机理进行了探索[12-14]。通过结合具有简单几何形状裂缝的有限元模型和平面压裂裂缝解析解，TAGHICHIAN 等[15]发展了应力阴影和裂缝相互作用分析的数值模型，研究了致密岩体中压裂裂缝间的应力阴影区域和裂缝扩展高度。基于增强的二维位移不连续性方法，KRESSE 等[16]建立了复杂水压裂缝网络中分支裂缝应力阴影和裂缝扩展高度修正模型，研究裂缝的扰动偏转行为。基于有限元和边界元方法，KUMAR 等[17]模拟分析了水平井多射孔簇的裂缝扩展和地应力场重分布行为。自适应有限元方法在振动、稳定、损伤分析等复杂问题中展现出很好的求解效力[18-21]，传统有限元因在裂缝尖端区域网格划分的限制[22]，较难可靠、有效求解断裂问题，该自适应分析方法使其成为可能。笔者利用网格自适应重划分方法，结合有限元和离散元耦合技术，研究了考虑热扩散效应的超临界CO2压裂裂缝平面内动态扩展[23]、水力压裂裂缝三维扩展演化行为[24]。

笔者综合考虑热-流-固耦合效应，进一步研究水力压裂缝网三维扩展中的应力阴影效应和多裂缝扰动偏转行为，引入自适应有限元-离散元算法，通过三维裂缝尖端局部区域的自适应网格重划分获得高精度应力解答并得以有效描述裂缝动态扩展，形成分析策略和求解方案。该研究建立水平井分段压裂的工程尺度三维数值模型，利用典型工况计算分析了压裂裂缝三维扩展的热扩散效应影响、不同射孔簇间距以及不同分段压裂方案（顺序、同步、交替压裂）下裂缝网络的空间扩展扰动行为。

1 岩体压裂的热-流-固耦合控制方程

笔者研究岩体压裂过程涉及到的物理场包括温度场、流体场、固体场，下面分别对岩体基质的固体变形、孔隙-裂隙内的流体流动和热扩散的微分控制方程进行介绍。

1.1 固体变形

岩体基质考虑为多孔弹性介质，变形控制方程[23-24]为

1.2 流体流动

1.3 热扩散

岩体基质、孔隙-裂隙内流体之间的热扩散控制方程[23]为

式中：kb为 导热系数；Tf为流体温度； ρb为容积密度；cb为比热系数； ρf为 流体密度；cf为 流体比热系数；qf为Darcy 流体通量。

热-流-固耦合作用和参数传递中，流体与固体的温度传递通过热扩散实现，温度变化引起固体的应力场变化；孔隙流体对固体的作用通过有效应力原理实现，流体向固体传递孔隙流体压力场；裂隙流体对固体的作用通过裂隙壁面作用的水压力实现，流体向固体传递裂缝内流体压力场；固体对孔隙流体的作用通过体积应变实现，固体向流体传递应变场。

2 裂缝扩展过程热传递与网格自适应分析

2.1 压裂裂缝扩展准则

记单元的抗拉强度和断裂能分别为 σt、Gf，单元应力-应变关系如图1 所示，当其最大拉应力达到 σt，则单元开始损伤(d=0)，拉应力为0 时，单元损伤达到最大(d=1)，将发生断裂[24]。该应力-应变曲线与x轴所围面积即为断裂能Gf，应力-应变曲线在损伤阶段的斜率为H，损伤阶段弹性模量为E～，E为弹性阶段的弹性模量。裂缝扩展中实施网格重分[24]：在流体驱动裂缝扩展过程中，在裂缝尖端出现应力集中，形成损伤区(d=1)；损伤区表示单元已经破坏，裂缝沿垂直最大主应力方向扩展，裂缝长度即该方向在损伤区的直线距离；当裂缝预测长度达到提前给定的扩展长度时，采用常规离散元处理单元间破裂的方式进行裂缝扩展，扩展长度即裂缝预测长度，网格在裂缝尖端指定区域内进行重新划分，并计算得到新的损伤区。

图1 应力-应变关系及损伤分析Fig.1 Relationship of stress-stain and damage analysis

2.2 单元离散及节点间热传递

岩体基质的固体变形、孔隙-裂隙内的流体流动和热扩散的微分控制方程式(1)-式(3)采用常规有限元法进行离散，推导过程及离散后的常微分方程组形式可参见文献[24-25]。通过有限元单元离散后，岩石基质区域与裂缝区域利用裂缝面节点进行热传递[23]，单元节点间的热传递形式如图2 所示，节点处的温度和热通量[23]为

图2 单元节点间的热传递Fig.2 Heat transfer between nodes

式中：ΔT为岩石单元温度变化量；ΔV为体积变化量；V为初始体积； αT为岩石基质的线性热膨胀系数。其通过以上所述的储层岩体热传递诱发体积收缩或膨胀，从而引起岩体固体应变场、应力场变化。

2.3 网格局部自适应重划分方法

研究采用超收敛拼片法[26]计算得到的应力超收敛解比原有限元解具有至少高一阶的收敛阶，可以用其对原有限元解进行误差估计和控制。在每个单元K上定义误差估计参数[18,20-21]：

式中：hK为 单元K的当前网格尺寸。需要指出的是，为了更合理地生成网格、高效地进行计算，可以仅在裂缝尖端区域进行网格重划分。对裂缝尖端的单元进行误差估计和网格重划分，使三维问题的自适应求解大大降低运算时间，同时也能保证很好描述裂缝扩展路径。

3 水平井分段压裂三维数值模型

研究建立图3 所示深部致密岩体水平井多射孔簇分段压裂的工程尺度初始几何模型，模型尺寸为400 m×600 m×400 m。模型的基本物理参数见表1，部分选自文献[23]中致密页岩材料参数。压裂中设置5 个射孔簇位置，射孔簇编号依次为1～5。根据射孔簇不同压裂扩展的顺序，顺序压裂方案的顺序为1→2→3→4→5，交替压裂方案的顺序为1→3→2→5→4，同步压裂方案的顺序为1–2–3–4–5（即5 个射孔簇同时压裂扩展）。表2 所示为顺序和同步压裂的各阶段持续时间和总时间。5 个射孔簇沿水平井段均匀分布，坐标位置为yi（i为射孔簇编号值，取为1～5）：

表1 模型基本物理参数Table 1 Basic physical parameters of the model

表2 顺序和同步压裂的各阶段持续时间和总时间Table 2 Duration and total time of multiple fracturing stages for sequential and alternate fracturing

式中：a为射孔簇间距。

为研究射孔簇间距对裂缝扩展的影响，对表3中列出的100、75、50、25 m 等射孔簇间距的压裂工况进行计算分析。表4 为不同压裂液和岩石基质的温度梯度工况，用于研究储层温度梯度对裂缝扩展的影响。需要说明的是，研究中主要针对这两种考虑温度效应（热-流-固耦合）工况与不考虑温度效应（仅流-固耦合）工况进行计算与对比分析。本文集中在讨论有无温度效应、温度梯度不同情况下裂缝的扩展行为，并没有针对不同温度梯度下研究裂缝的不同扩展行为。为可靠描述初始阶段的裂缝起裂，在射孔区域进行初始网格细化，该区域为150 m≤x≤250 m、(y1– 50) m≤y≤(y5+50) m、150 m≤z≤250 m。研究采用线性四面体单元对模型进行离散化处理，图4 为模型网格划分截面图，在每个射孔区域采用了更为细密的初始网格。上述水平井分段压裂三维数值模型的计算通过岩石力学计算分析程序软件ELFEN[24]实施。

表3 各压裂方案射孔簇间距设置Table 3 Perforation cluster spaces for fracturing scenarios

图4 初始网格Fig.4 Initial mesh

4 水力压裂缝网三维扩展扰动模拟分析

4.1 不同射孔簇间距下裂缝网络扩展扰动

对数值模型压裂过程分析的网格自适应优化细分，并使用各优化网格对不同射孔簇间距下裂缝网络扩展扰动行为进行研究。

图5 为顺序压裂下裂缝区域局部网格重划分，分别为第一、第三、第五阶段的裂缝扩展和局部网格重划分结果。外部粗糙网格为储层模型外部使用的网格，裂缝周围区域使用了更为细密的网格。可以看出，随着裂缝的动态扩展，裂缝周围区域的网格进行自适应重划分加密，形成的高质量网格确保裂缝尖端区域应力场解答的精确性和裂缝扩展路径的可靠性。

图5 顺序压裂下裂缝区域局部网格重划分Fig.5 Local mesh refinement around fracture domains in sequential fracturing

图6 为顺序压裂方案中不同射孔簇间距下的最终缝网形态和应力场结果，可以看出裂缝出现不同程度的相互扰动偏转。图6a 为射孔簇间距a=100 m的裂缝扩展形态，由于各射孔簇间距较大，裂缝接近平行稳定扩展。图6b 为射孔簇间距a=75 m 的裂缝扩展形态，由于各射孔簇间距缩小，裂缝间扰动应力场出现相互干扰，从第2 条裂缝朝左侧（远离其它裂缝侧）开始出现轻微偏转，后续裂缝扩展扰动累计、偏转量增大。图6c 为射孔簇间距a=50 m 的裂缝扩展形态，由于各射孔簇间距继续缩小，裂缝间扰动应力场出现叠加和干扰，从第2 条裂缝开始偏转加剧。图6d 为射孔簇间距a=25 m 的裂缝扩展形态，由于各射孔簇距离骤减，裂缝间扰动应力场出现强烈干扰，从第2 条裂缝开始剧烈偏转，射孔簇间距成为影响三维裂缝非稳定扩展的重要因素。

图6 顺序压裂下最终缝网形态和应力场结果Fig.6 Final morphology of fracture network and stress results in sequential fracturing

4.2 顺序、同步、交替压裂下裂缝网络扩展扰动

在本研究模型物理力学参数和工况设置下，可以看出射孔簇间距a=75 m 时，裂缝出现空间扩展扰动现象；下面基于该射孔簇间距，讨论顺序、同步、交替压裂下裂缝网络扩展扰动行为。

图7 为顺序压裂下各阶段裂缝扩展和应力场演化结果。图7a 为第一阶段裂缝展布情况，可以看出由于没有其它裂缝的干扰，第1 条裂缝在空间接近平面扩展。图7b 为第二阶段裂缝形态，可以看出由于第1 条裂缝对地应力场的扰动，使得第2 条裂缝在空间出现轻微偏转。图7c 为第三阶段裂缝形态，可以看出由于第1 条、第2 条裂缝对地应力场的扰动，使得第3 条裂缝偏转加剧。图7d、图7e 为第四、第五阶段裂缝形态，可以看出剩余的射孔簇压裂过程中出现累积偏转，裂缝均朝同侧偏转，成为顺序压裂的典型裂缝扩展行为。

图7 顺序压裂下裂缝扩展和应力场演化结果(射孔簇间距a=75 m)Fig.7 Dynamic propagation of fracture network and evolution of stress results in sequential fracturing (Perforation cluster space a=75 m)

图8 为同步压裂下各阶段裂缝扩展和应力场演化结果。图8a 为中间阶段(t=1 210 s)裂缝展布情况，可以看出各裂缝接近平面稳定扩展；图8b 为最终时刻裂缝形态，可以看出位于射孔簇两侧的第1 条、第5 条裂缝分别向外侧偏转，成为同步压裂的典型裂缝扩展行为。

图8 同步顺序压裂下裂缝扩展和应力场演化结果Fig.8 Dynamic propagation of fracture network and evolution of stress results in simultaneous fracturing

图9 为交替压裂下各阶段裂缝扩展和应力场演化结果。图9a 为第一阶段裂缝展布情况，该阶段与图7a 所示的第一阶段压裂条件相同，取得了相同的压裂结果；即由于没有其它裂缝的干扰的原因，第1条裂缝在空间接近平面扩展。图9b 为第二阶段裂缝形态，可以看出由于第3 条裂缝开始交替扩展；此时，第1 条、第3 条裂缝间距为原始射孔簇间距的2 倍，降低了裂缝之间的相互扰动，使得第3 条裂缝仍呈平面扩展。图9c 为第三阶段裂缝形态，可以看出介于第1 条、第3 条裂缝中间的第2 条裂缝开始扩展，在两侧裂缝相当的影响下，第2 条裂缝最终呈平面展布。同理，第4 条、第5 条也进行交替式压裂扩展，如图9d、图9e 所示裂缝扩展过程仅出现轻微偏转。综上所述，在未改变射孔簇间距的基础上，通过控制起裂顺序，多裂缝得以稳定扩展；射孔簇起裂顺序的变化使得射孔簇起裂时的间距增大，进而减弱多裂缝的偏转程度和应力阴影效应。

图9 交替压裂下裂缝扩展和应力场演化结果(射孔簇间距a=75 m)Fig.9 Dynamic propagation of fracture network and evolution of stress results in alternate fracturing (Perforation cluster space a=75 m)

为定量分析裂缝扩展行为，本研究统计了表5所示的各压裂方案下压裂裂缝面积和体积结果，给出射孔簇间距100 m 和75 m 的结果。可以看出在不同压裂方案中，压裂裂缝最终总面积和体积随射孔簇间距的减小均有不同程度降低，射孔簇间距在起控制作用。对比各压裂方案，可以看出同步压裂比顺序压裂得出更小的裂缝面积和体积，交替压裂可以得出更大的裂缝面积和体积，射孔簇的起裂顺序在起控制作用。在研究裂缝之间显著的影响作用以及应力干扰效应时，本研究将射孔簇间距继续减小到a=50、25、12.5 m 时，裂缝面积和体积将出现非常显著的差异性，得出射孔簇间距对裂缝偏转、应力场干扰的影响[27]。综上，在相同射孔簇间距设置下，交替压裂成为缓解压裂缝网三维扩展扰动、优化空间缝网形态的有效方案。

表5 各压裂方案下裂缝面积和体积结果Table 5 Results of fracture areas and volumes in sequential fracturing

4.3 压裂裂缝三维扩展的热扩散效应影响

岩体基质、孔隙-裂隙内流体之间发生热扩散，图10 为压裂裂缝区域的热扩散行为，压裂液和岩石基质之间的温度梯度使得三维裂缝面环状区域发生热扩散。

图10 压裂裂缝区域热扩散Fig.10 Thermal diffusion around fracture domains

本文研究的储层深度约2 000 m，在该储层上开展水平井分段压裂，在该深度的储层温度约60 ℃，温度梯度为30 ℃/km。为分析深部致密岩体裂缝内的压裂液与岩体基质的热交换，本研究对比储层和压裂液不同温度梯度（60、20 ℃；60、35 ℃）下压裂裂缝面积和体积，定量分析热扩散影响压裂裂缝扩展的行为。为分析热扩散效应的影响，本研究对考虑温度场的热-流-固耦合作用、不考虑温度场的流-固耦合作用的压裂过程进行计算。图11 为不同压裂方案下压裂裂缝面积对比，可以看出考虑温度场作用的最终压裂面积均大于不考虑温度场作用的情况，忽略热扩散效应将低估压裂缝网扩展。对于图11a所示裂缝扩展扰动弱的工况（射孔簇间距a=100 m）和图11b 所示扩展扰动增强的工况（射孔簇间距a=75 m），裂缝面积出现一致的结果：温度梯度越大，热扩散效应越强；相比顺序和同步压裂，对于有利于裂缝稳定扩展的交替压裂方案，热扩散效应提高裂缝扩展的作用显著。

图11 不同压裂方案下压裂裂缝面积对比Fig.11 Comparison of fracture areas in different fracturing scenarios

图12 为不同压裂方案下压裂裂缝体积对比，可以看出考虑温度场作用的最终压裂体积均大于不考虑温度场作用的情况，忽略热扩散效应将低估压裂缝网扩展。对于图12a 所示裂缝扩展扰动弱的工况和图12b 所示扩展扰动增强的工况，裂缝体积出现一致的结果：温度梯度越大，热扩散效应越强；顺序、同步、交替压裂中，热扩散效应促进裂缝扩展。

5 结 论

1）压裂裂缝动态扩展引起的应力扰动区域在多裂缝中存在叠加、覆盖行为，形成应力阴影效应、造成裂缝空间偏转。

2）水平井多射孔簇间距的减小，将增大应力阴影区、加剧裂缝间相互干扰。

3）相比多射孔簇顺序压裂，同步压裂将增大应力阴影区、交替压裂可减小应力阴影区，交替压裂成为缓解压裂缝网三维动态扩展扰动、优化空间缝网形态的有效方案。

4）深部致密岩体裂缝内的压裂液与岩体基质进行热交换，对比储层和压裂液不同温度梯度（60、20 ℃；60、35 ℃）下压裂裂缝面积，定量分析了热扩散促进压裂裂缝扩展的行为。压裂液与岩石基质的温度梯度越大，热扩散效应越强；顺序、同步、交替压裂中，热扩散效应对裂缝扩展有促进作用，增大裂缝扩展面积和体积。