水平井压裂技术已经成为开发低渗透油气藏的关键技术[1], 由于单簇射孔储层改造体积的局限性, 发展形成了水平井多段分簇压裂技术。多簇水力裂缝改变了原地层应力场分布, 产生的诱导应力使水力裂缝偏离最大主应力方向扩展, 裂缝形态复杂[2]。常规的数值计算方法无法实现水力压裂过程中的流-固、孔压-岩石骨架和应力场等多场耦合问题的模拟, 而ABAQUS能够很好地解决这一问题。对于水力压裂多场耦合, 张汝生等[3]采用ABAQUS的cohesive单元模拟了水力压裂扩展过程, 但不涉及缝间干扰问题; 对于诱导应力场引起的缝间干扰, 潘林华等[4]用ABAQUS得出了水平井分簇压裂裂缝起裂和扩展规律, 但不涉及裂缝转向等动态扩展问题。笔者基于ABAQUS扩展有限元法[5-6], 根据流-固耦合方程和岩石断裂力学原理[7-8], 在前人研究的基础上通过用户子程序开发端口, 模拟水平井多段多簇压裂过程中储层岩石孔隙压力、渗透率、孔隙度、裂缝面滤失等动态变化过程, 分析水平应力差、簇间距、裂缝长度和压裂次序下缝间干扰对裂缝扩展规律及裂缝形态的影响, 并利用数值模拟结果对现场实例进行优化, 验证模型的适用性。

从力学方面分析, 水力压裂扩展有限元数学模型[9]的建立包括:①流-固耦合，压裂液及孔隙流体对多孔介质及裂缝面应力应变的影响; ②裂缝起裂, 水力压力克服地应力达到岩石的抗张强度使岩石破裂; ③裂缝扩展, 岩石破裂后水力压力达到岩石扩展所需最小压力, 裂缝向前延伸。

ABAQUS应用岩石的虚功与作用在岩石上的体力和面力产生的虚功相等, 根据有效应力及虚功原理简化得到描述水力压裂过程中随储层断裂岩石骨架变形的岩石应力平衡方程[10], 表达式为

式中, $ \bar \sigma $为储层基质的有效应力, Pa; pw为储层孔隙压力, Pa; m为[1, 1, 1, 0, 0, 0]; δε为储层基质的虚应变率, s-1; t为单元面载荷, N/m2; δv为储层基质虚运动速度, m/s; f为体力载荷, N/m3; dS为载荷作用面积, m2; dV为计算元的体积, m3。

假定水力压裂过程中渗流过程符合达西定律, 根据流体质量守恒原理, 一定时间内流入基质的流体质量总和等于基质内流体增量与流出质量之和, 推导得出流体连续性方程[10]为

式中, J为储层孔隙的体积改变率; nw为孔隙中液体体积与储层总体积的比值; ρw为储层液体密度, kg/m3; x为液体流动的方向, m; dt为时间步, s; νw为储层液体的流动速度。

水力压裂是裂缝动态扩展过程, 随着裂缝起裂并不断向前延伸, 储层孔隙压力及渗透率发生改变。孔隙度和渗透率动态变化方程[11]为

式中, φ为储层动态孔隙度; φ0为储层初始孔隙度; εv为计算的体积应变; k为储层动态渗透率, 10-3 μm2; k0为储层初始孔隙度, 10-3 μm2。

模型采用扩展有限单元模拟岩石断裂及裂缝扩展过程, 通过ABAQUS隐式求解器, 将式(1)~(4)离散为关于储层岩石应力、应变、位移、孔隙度、渗透率等相关参数的矩阵方程, 通过求解位移变量, 得出相关参数变化值。其中单元节点位移由连续位移场相关函数通过扩展有限元方法求得。扩展有限元法采用裂尖和裂缝面两个水平集函数确定裂缝的位置, 并在形函数的基础上总结形成能够描述裂缝扩展有限元模型的连续位移场函数[12], 表达式为

式中, N(x)为有限元形函数; F(x)为裂尖形函数; ui为标准自由度; I为所有有限元节点集合; J为裂缝穿过单元的节点集; H(x)为跃阶扩展形函数; aj、bk, l1和bk, l2分别为裂缝面、裂缝尖端一型裂缝和二型裂缝的自由度。

目前扩展有限元分析中判断裂缝起裂的准则主要有最大主应力准则、最大主应变准则、二次应力准则等, 经过大量数值模拟研究, 笔者认为最大主应力准则适合于水力压裂数值模拟中裂缝动态扩展研究。当水力压力大于储层岩石的临界最大主应力时, 损伤开始, 岩石破裂, 最大主应力准则[13]表达式为

式中, σmaxc为储层岩石承受的临界最大主应力, MPa, 超过此临界值, 即f > 1时岩石断裂, 损伤开始; σmax为储层岩石所受的最大主应力, MPa。

水力裂缝扩展是一个复杂的过程, 无法提前确定裂缝生成的类型, 在最大主应力准则的基础上引入能量释放率判断复合裂缝生成后裂缝的扩展。B-K准则[14]表达式为

式中, Gs和GT分别为两个切向的断裂能量释放率, N/mm, B-K准则认为两个切向的能量释放率相等; η为与材料特性相关的常数; Gc为复合型裂缝断裂临界能量释放率, N/mm; Gcs和Gcn分别为切向和法向临界能量释放率, N/mm。

由于室内试验无法模拟真实条件下的水平井分段分簇压裂的裂缝扩展, 所以使用室内大型真三轴水力压裂物模试验验证ABAQUS扩展有限元数值模型在水力压裂实践中的正确性。分别建立水力压裂数值模型和真三轴试验模型, 模型尺寸为30 cm×30 cm×30 cm, 井径为20 mm, 射孔方位角为45°, 两个模型采用相同的试验参数, 其中弹性模量为16.14 GPa, 泊松比为0.18, 最大和最小水平主应力分别为15和9 MPa, 岩石渗透率为15×10-3 μm2, 孔隙度为0.12, 压裂液黏度为40 mPa·s, 试验排量为100mL/min。模拟结果表明, 数值模型和真三轴模型在有应力差条件下, 裂缝均发生了转向, 且转向角及裂缝扩展形态趋于一致。这说明采用扩展有限元模拟裂缝转向及其动态扩展规律具有一定的可靠性, 能够应用于水平井多簇分段压裂中的裂缝扩展动态研究, 结果对比见图 1。

根据胜利油田某区块的地质资料建立水平井多段多簇压裂二维数值模型, 基本模型参数为:岩石弹性模量为12.79 GPa, 泊松比为0.25, 水平主应力为10~18 MPa, 储层上覆应力为20 MPa, 初始孔隙压力为10 MPa, 岩石抗拉强度为1.82 MPa, 岩石虑失系数为5.879×10-5 m3/min, 压裂液排量及黏度为6 m3/min和40 mPa·s, 孔隙度为0.25, 射孔长度为0.5 m, 储层渗透率为20×10-3 μm2。基于ABAQUS扩展有限元法, 使用Geostatic和Soil模块模拟真实地层条件下的多场耦合问题, 同时内置孔压单元模拟地层岩石孔隙及线弹性力学性质。模型建立过程中为了克服边界及远场地应力影响, 采用长度为160 m、宽度为80 m的计算模型, 同时根据对称原理选取1/4模型以减小计算代价。为了保证模型的计算精度, 对压裂段网格进行局部加密, 并引进线性减缩积分单元进行计算。实际模型中用预置裂缝模拟射孔簇, 其中水平井方向沿最小主应力方向, 射孔方向沿最大主应力方向，如图 2所示。

水平井多段分簇压裂中, 由于采用多簇射孔, 水力压裂过程中多条水力裂缝产生缝间应力干扰, 导致地应力场发生改变, 水力裂缝不再沿最大主应力方向扩展, 而是随着裂缝间距、裂缝长度、水平应力差等因素的改变出现不同程度的转向[15-16]。如果施工参数不合理, 可能会使裂缝急剧转向[17], 导致裂缝局部摩阻过大, 裂缝尖端所需的破裂能量增大, 同等排量下主力裂缝的有效长度变短。

在水力压裂初期, 水力裂缝沿着射孔方向向前延伸, 裂缝尖端出现应力集中, 当最大拉应力大于岩石的抗张强度时, 裂缝开启。随着压裂过程的延续, 裂缝间的干扰逐渐增强, 裂缝尖端主应力及方向急剧改变, 裂缝出现不同程度的转向, 表现为水平诱导应力差值越大, 裂缝偏转程度越大, 裂缝尖端应力集中越明显, 起裂所需水力能量越大。如图 3(泵注时间100 s, 簇间距15 m)所示。

水力裂缝长度是影响压裂成功与否的关键因素, 也是多段分簇压裂设计和评价的核心部分。在射孔簇数和间距一定的情况下, 裂缝长度增加, 诱导应力场的控制面积增加, 裂缝形态复杂程度增加。图 4为簇间距15 m、地应力差2 MPa时裂缝扩展规律。定义裂缝与最大主应力方向的夹角为裂缝偏转角, 裂缝偏转位置距射孔根部的距离为裂缝转向位置。由图 4可以看出:在水力压裂初始阶段(图 4(a))诱导应力场影响范围较小, 裂缝沿射孔方向扩展, 裂缝所需延伸压力低, 易形成有效裂缝; 随着裂缝向前延伸, 诱导应力场控制面积增加, 而裂缝尖端由于应力场干扰, 发生轻微扰动, 两条主裂缝开始向相反方向偏转, 偏转角为7°(图 4(b)); 裂缝长度增加, 诱导应力影响范围增加, 主裂缝偏转程度也增加(图 4(c)、(d))。这说明随着水力裂缝的扩展, 裂缝偏转程度越大, 裂缝形状越复杂, 裂缝扩展所需能量越大, 相同泵注程序下形成的有效裂缝长度越小。

水平井多簇压裂的簇间距D是水力压裂优化的重要参数, 直接影响压裂效果, 决定压裂增产效果。簇间距越小, 裂缝间的干扰越严重, 裂缝越容易转向, 裂缝形态越复杂。图 5为簇间距对裂缝扩展影响。其中压裂液注入时间为100 s, 水平应力差为2 MPa。由图 5可以看出:簇间距为10 m时, 水力裂缝在压裂初期明显偏转, 随着裂缝长度增加，偏转程度变大; 随着簇间距增大, 裂缝偏转程度减小, 裂缝偏转位置向后推移。4种簇间距情况下, 裂缝转向位置分别为3.2、4.2、5.8和16 m, 裂缝偏转角分别为22°、19°、10°和3°, 说明随着射孔簇间距的增加, 诱导应力的影响减小。图 6为相同泵注程序下两簇水力裂缝偏转位置和缝长随簇间距的变化。由图 6可以看出, 裂缝转向位置和裂缝长度随簇间距的增加呈递增关系, 且转向位置随间距的增加变化明显, 而裂缝长度随簇间距的增加变化程度较小, 说明在大簇间距情况下, 诱导应力影响递减明显, 裂缝扩展主要受原地应力场控制。

水平应力差Dp是裂缝扩展难易的决定性因素, 也是多段分簇压裂中必须考虑的问题。水平应力差越大, 诱导应力对裂缝扩展产生的影响越小, 在高应力差情况下, 裂缝趋于最大主应力方向扩展, 转向不明显, 裂缝形态单一。图 7为两簇裂缝随地应力差变化, 其中簇间距为10 m, 泵注时间为100 s。由图 7可以看出:在低水平应力差情况下, 裂缝扩展产生的诱导应力场影响显著, 水力裂缝转向明显, 且转向位置离射孔根部较近, 转向位置为2.2 m; 随着水平应力差增加, 裂缝偏转程度减小, 偏转位置距射孔根部距离增加, 诱导应力场影响程度减弱。图 8表明, 裂缝转向位置和裂缝长度随水平应力差的增加呈递增关系, 但在高应力差情况下, 诱导应力对裂缝形态的影响明显减弱, 裂缝趋于最大主应力方向扩展, 偏转程度较弱。水平应力差为2~8 MPa, 裂缝转向位置增加了13.8 m, 最长裂缝长度增加了10.2 m。

随着压裂次序的改变, 裂缝周围应力场出现不同程度的改变, 合理利用能有效优化裂缝长度。图 9为不同压裂次序下裂缝扩展情况, 其中射孔簇间距为15 m, 水平应力差为2 MPa, 泵注时间为300 s。将射孔簇从右到左编号(图 9(a)), 由图 9可以看出:将①、②、③射孔簇同时压裂, 由于①和③射孔簇的应力场叠加, 中间射孔簇裂缝沿最大主应力方向延伸, 不发生偏转, 而①和③射孔簇由于中间射孔簇诱导应力场的影响, 发生不同程度的偏转, 并且中间射孔簇的裂缝长度比①和③的裂缝长度增加了2.5 m; 将①、②、③依次分步压裂, 射孔簇①没有诱导应力场影响, 沿最大主应力方向扩展, 不发生偏转; 射孔簇②受到①所产生的诱导应力, 在裂缝长度达到4.2 m时发生偏转, 而射孔簇③同时受到①和②所产生诱导应力影响, 在3.8 m处发生偏转, 且裂缝弯曲程度明显高于②所产生的裂缝。在相同泵注程序下裂缝长度的关系为① > ② > ③, 分别为21.7、21.1和20.2 m。说明随着射孔簇数的增加, 裂缝产生的诱导应力场相互叠加, 缝间诱导应力场的影响范围扩大。图 9(d)为两步法压裂裂缝扩展情况, 首先将①和③射孔簇压裂, 再压裂射孔簇②, 3条主裂缝均沿最大主应力方向扩展, 偏转程度及裂缝复杂程度均小于同步压裂和顺序压裂, 裂缝总长度关系为两步法压裂 > 顺序压裂 > 同步压裂, 分别为72.2、63和59.8 m。说明两步法压裂能减小诱导应力场对水力裂缝的影响, 增加有效裂缝形成的概率。

应用数值模拟结果对胜利油田某区块的2口水平井的各段射孔和射孔簇进行优化。水平井所在目标储层厚度为32 m, 井深1 428 m, 上下都有隔层。该水平井邻近3口井都进行了多段分簇压裂(顺序压裂), 其中水平井长度为600 m, 分5段射孔, 优化簇间距为20 m, 每段射孔含3个射孔簇, 压裂施工后平均单井产量为15.8 m3/d。对每簇射孔的流量监测表明:每段射孔簇产量从右到左分别占总产量的比例为50%、30%和20%, 即80%以上的产量均由顺序压裂的前两段射孔提供。产量监测结果表明, 由于诱导应力场的叠加, 压裂次序最后的裂缝形成有效裂缝的概率降低, 较先压裂的裂缝产能降低。

根据储层物性参数, 对要优化的井进行多段分簇压裂数值模拟, 射孔参数与前3口压裂井一致, 改变压裂方法, 即使用两步法压裂对每个射孔段进行压裂。压裂施工后产量监测发现, 每个射孔簇产量差别不大, 但是优化后产量比邻井油气总产量更高, 如图 10所示。其中1、2、3井为顺序压裂井, 4和5两口井为两步法压裂井, 采用两步法压裂单井平均产量增加了2.4 m3/d, 增产效果明显。这也说明了优化参数形成的有效裂缝长度增加, 间接证明了数值模型的正确性, 对现场应用有一定的指导意义。

(1) 裂缝长度是影响诱导应力场控制面积的重要因素。随着裂缝长度的增加, 诱导应力场控制面积增加, 裂缝背离彼此偏转的程度增加, 扩展所需能量变大。

(2) 簇间距对诱导应力场的影响明显。射孔簇间距增加, 裂缝之间的诱导应力影响减弱, 裂缝转向位置增加, 偏转程度减小。在射孔簇间距为30 m时, 裂缝偏角仅为3°, 对有效裂缝长度的影响较弱。

(3) 水平应力差是诱导应力场影响程度的决定性因素。裂缝转向位置和裂缝长度随水平应力差的增加呈递增趋势, 但在高应力差情况下, 诱导应力对裂缝形态的影响明显减弱, 裂缝趋于最大主应力方向扩展, 偏转程度较小。水平应力差为2~8 MPa时, 裂缝转向距离增加了13.8 m, 最长裂缝长度增加了10.2 m。

(4) 改变簇间压裂次序是合理控制诱导应力影响程度的重要途径。在每段3簇射孔水平井压裂中, 两步法压裂裂缝复杂程度较低, 形成的有效裂缝长度大于同步压裂和顺序压裂。