0 引言
针对水力裂缝扩展这一物理行为，前人进行了较多研究。多数学者将储层考虑为均质地层，在不考虑天然裂缝影响下，研究了多簇裂缝同步或拉链式扩展的物理过程；少数学者虽考虑了储层弱面结构影响，但未研究拉链式压裂裂缝扩展及窜通行为。Roussel等基于位移不连续法，不考虑天然裂缝影响，模拟研究裂缝起裂顺序、簇间距、应力差、排量等因素对裂缝扩展的影响规律，发现交替压裂作业可使裂缝扩展更均匀，同步压裂则会导致两侧裂缝优势扩展，中间裂缝受抑制，大簇间距、高应力差、较小排量下，裂缝扩展整体更均匀[6⇓⇓-9]。Zou等研究了天然裂缝存在时天然裂缝胶结强度、逼近角、地应力差、施工排量等因素对单簇/多簇裂缝扩展规律以及水力裂缝与天然裂缝相互作用机制的影响[10⇓⇓⇓⇓-15]，但未涉及双水平井拉链式压裂时裂缝的扩展规律。Hou等建立了考虑天然裂缝影响的裂缝扩展模型，分析了多裂缝竞争起裂时裂缝的扩展过程，发现天然裂缝对裂缝扩展轨迹及复杂缝网的形成有重要影响，但未进一步探讨拉链式压裂裂缝的扩展过程[16⇓⇓⇓-20]。Wu等研究了均质地层井间单簇和多簇裂缝的扩展规律，探讨了压裂顺序、井间距、簇间距等因素影响下拉链式压裂裂缝的扩展规律[21⇓⇓⇓-25]。上述研究虽然考虑了拉链式压裂的影响，但受数值模拟技术的限制，并未探索天然裂缝、裂缝带对拉链式压裂裂缝延伸与窜扰规律的影响。
综上可知，目前主要探索了单井或多井的单簇、多簇裂缝扩展规律，且不考虑天然裂缝影响。当存在天然裂缝带时，拉链式压裂裂缝的扩展与窜扰规律研究还有待深入。因此，本文采用有限元-离散元法（FDEM），建立拉链式压裂的流-固耦合数值模型，根据川南W平台的地质特征参数，研究天然裂缝带对井间裂缝扩展规律的影响，探索存在天然裂缝带时井间裂缝的窜通机理，以期为井工厂拉链式压裂施工提供技术支持。
1 模型构建
1.1 概念模型
井工厂开发模式已被广泛应用于页岩气开发，页岩气储层通常天然裂缝发育。为探究天然裂缝带对拉链式压裂裂缝扩展与窜扰规律的影响，建立如图1 所示的拉链式压裂概念模型，其中天然裂缝带通过离散天然裂缝集群近似表征。基于该模型模拟现场真实拉链式压裂的3个过程：水平井a压裂、停泵、水平井b压裂。通过3个过程的顺序耦合模拟，可捕捉压裂井与响应井之间的干扰动态。
1.2 数学模型
本文模型在考虑流-固耦合基础上，克服常规模型尺度小、计算成本高等缺陷，采用过程顺序耦合方法，将响应井压裂、关井、激动井压裂过程中应力、裂缝轨迹、孔隙压力等参数，依据矿场作业时间顺序，进行自然耦合，不需要进行数据处理或者外部调用，实现了大尺度、多过程的连续模拟。
FDEM的基本思想是将模拟域连续体离散为有限数量的三角形单元网格，并在所有相邻单元公共边上插入节理单元（见图2a ）。由于节理单元具有黏结作用，基质与天然裂缝的起裂与扩展都可以通过节理单元的损伤破坏来表征。在节理单元未破坏时，基质的变形利用具有较强黏聚性的节理单元进行模拟；天然裂缝通过具有初始厚度、黏聚力的节理单元表征（见图2a ）；破裂后的节理单元则形成水力裂缝。
1.2.1 应力平衡方程
任意连续体的应力分布都可由位移计算得到，而位移计算可归结为节点位移的计算。根据FDEM，三角形单元的质量、载荷、位移、速度、加速度场可以用节点值表示。根据节点质量、接触力、三角形单元变形引起的反作用力和节点单元的黏结应力，结合牛顿第二定律，这些节点的位置和速度可以由运动方程确定[26]：
$m\overline{\overline{x}}\text{+}C\overline{x}={{F}_{\text{c}}}\left( x \right)+{{F}_{\text{d}}}\left( x \right)+{{F}_{\text{j}}}\left( x \right)+{{F}_{\text{ext}}}\left( x \right)$
1.2.2 接触力计算方程
${{f}_{\text{n}}}={{p}_{\text{n}}}\int_{{{\beta }_{\text{t}}}\cap {{\beta }_{\text{c}}}}{\left[ \text{grad}{{\varphi }_{\text{c}}}\left( {{P}_{\text{c}}} \right)-\text{grad}{{\varphi }_{\text{t}}}\left( {{P}_{\text{t}}} \right) \right]}\text{d}A$
Mahabadi等在二维FDEM中根据摩尔库伦摩擦准则，提出了切向接触力的计算方法[29]：
${{f}_{\text{t, }n+1}}={{f}_{\text{t, }n}}-{{p}_{\text{t}}}\Delta {{u}_{\text{t}}}$
1.2.3 三角单元应力计算
三角形单元的总应变可以通过三角形单元当前位移和初始位移计算得到，然后根据应力-应变本构关系，得到具有黏性阻尼的各向同性线弹性三角形单元的应力[28]：
$T\text{=}\frac{1}{\sqrt{\left| \det \mathbf{F} \right|}}\left( \frac{E}{1+\nu }{{E}_{\text{d}}}+\frac{E}{1-\nu }{{E}_{\text{s}}}+2\lambda D \right)$
1.2.4 流体流动模型
基于霍夫定律，井底压力可以采用井筒沿程摩阻、孔眼摩阻及缝口压力表示：
${{p}_{\text{w}}}={{f}_{\text{w}}}+{{f}_{\text{p}}}+{{p}_{\text{f,}\ \text{in}}}$
各簇裂缝流入流量满足质量守恒定律：
$Q=\sum\limits_{i=1}^{{{n}_{\text{hf}}}}{{{Q}_{i}}}$
对于裂缝内流体流动，如图2a 所示，p1和p2分别为节点P1和P2处的压力，P2与P1间的总压差为[30]：
$\Delta p={{p}_{1}}-{{p}_{2}}+{{\rho }_{\text{w}}}g\left( {{y}_{2}}-{{y}_{1}} \right)$
根据立方定律，P2与P1间的流量为：
${{q}_{_{21}}}=\frac{1}{12\mu }{{a}^{3}}\frac{\Delta p}{l}$
断裂节理单元左右两侧与岩石基质单位时间内的流体交换量为：
$\left\{ \begin{align} & q_{\text{e}}^{+}={{h}_{\text{c}}}\left( p_{\text{p}}^{+}-{{p}_{\text{c}}} \right)L \\ & q_{\text{e}}^{-}=h_{\text{c}}^{{}}\left( p_{\text{p}}^{-}-{{p}_{\text{c}}} \right)L \\ \end{align} \right.$
断裂节理单元内流体在高压作用下进入基质孔隙，在FDEM中，孔隙渗流发生在相邻三角形单元间的未断裂节理之间。根据达西定律，单位横截面积的流速可表示为：
$v=\frac{K}{\mu }\frac{\partial p}{\partial x}{{f}_{\text{s}}}$
对于任何给定单元体，流体压力变化可由下式求得：
$\frac{\partial p}{\partial t}=\frac{M}{V}\left( {{Q}_{\text{total}}}+\alpha \frac{\partial V}{\partial t} \right)$
1.2.5 节理单元断裂模型
本文模型定义了3种断裂模式节理单元，分别是Ⅰ型（拉伸破坏）、Ⅱ型（剪切破坏）和Ⅰ-Ⅱ型（拉剪混合破坏）。通过定义损伤因子d（取值为0~1）来定量表征节理单元的破坏程度，当d=0时，表示没有损伤；当d=1时，表示完全损伤，即节理单元破裂。
$\frac{o-{{o}_{\text{p}}}}{{{o}_{\text{r}}}-{{o}_{\text{p}}}}=d$
$\frac{s-{{s}_{\text{p}}}}{{{s}_{\text{r}}}-{{s}_{\text{p}}}}=d$
${{\left( \frac{o-{{o}_{\text{p}}}}{{{o}_{\text{r}}}-{{o}_{\text{p}}}} \right)}^{2}}+{{\left( \frac{s-{{s}_{\text{p}}}}{{{s}_{\text{r}}}-{{s}_{\text{p}}}} \right)}^{2}}=d$
1.2.6 模型耦合框架
2 模型验证
2.1 水力裂缝与天然裂缝相交行为验证
水力裂缝与天然裂缝相交后扩展行为表征的准确与否是准确模拟裂缝性地层水力压裂过程的关键。Blanton通过物理实验获取的水力裂缝与天然裂缝相交准则得到了广泛认可，本节采用Blanton实验准则[31]验证本文模型。基于FDEM建立水力裂缝与天然裂缝相交模型，尺寸为10 m×20 m×10 m，在水平主应力差2 MPa、注入排量1 m3/min、黏度3 mPa·s条件下，模拟得到了不同裂缝逼近角条件下水力裂缝与天然裂缝相交后的扩展模式（见图5 ）。可以看到：裂缝逼近角为30°、45°时，水力裂缝被天然裂缝捕获；裂缝逼近角为60°时，水力裂缝则穿过天然裂缝。12组模拟结果（见图6 ）显示：在不同水平主应力差及逼近角条件下，水力裂缝与天然裂缝相交后的扩展模式完全按Blanton扩展模式分界线所划分的区域分布，说明本文模型是准确可靠的。
2.2 现场实例验证
川南W平台深层页岩气储层天然裂缝发育但分布规律性较差，受天然裂缝带发育的影响，该平台6口水平井压裂过程中，邻井压力涨幅显著，窜扰现象严重。为进一步验证本文模型的准确性，基于该平台2 井第13段（响应井）与3井第12段（激动井）地质参数，结合蚂蚁体显示的裂缝带分布，建立了如图7a 所示的拉链式压裂模型。根据测井曲线及实验数据可知，目标层段水平最大主应力为95 MPa，水平最小主应力为85 MPa，弹性模量为34 GPa，泊松比0.2。该压裂段设计5簇裂缝，施工排量14 m3/min，井间距300 m，考虑到裂缝带发育易导致压窜，注入总液量降低至1 600 m3。基于现场施工参数及本文建立的数值模型，模拟得到了目标层段3个压裂阶段的裂缝展布及地面压力曲线（见图7b —图7d 、图8 ）。从裂缝展布可以看到，在天然裂缝带影响下，两口井压裂裂缝最终窜通。将现场监测激动井、响应井的地面压力变化曲线与模拟地面压力数据进行对比，可以看到监测压力曲线与模拟压力曲线基本吻合。模拟与现场监测响应井压力涨幅分别为6.5 MPa、6.0 MPa，准确率为91.6%，说明本文模型具有较高的精确度。
3 裂缝扩展、窜通规律影响因素
3.1 模型构建及输入参数
根据W平台5井第20段（响应井）、6井第21段（激动井）物性参数，建立如图9 所示双水平井拉链式压裂物理模型，模型尺寸为100 m×200 m×10 m，划分离散单元数量为42 680，每段射孔3簇。天然裂缝带由集中分布的天然裂缝构成，天然裂缝逼近角与裂缝带逼近角保持一致，裂缝带处于两水平井之间。为准确观察拉链式压裂过程中的窜扰行为，分别模拟了响应井压裂、响应井压后停泵、激动井压裂3个连续阶段。两口井所有施工参数保持一致，基础参数如表1 所示。主要研究裂缝带逼近角、裂缝带内天然裂缝长度、裂缝带宽、裂缝带与孔眼相对位置对拉链式压裂裂缝扩展、窜扰规律的影响。
表1 W区块主要基础参数 |
| 参数名称 | 数值 | 参数名称 | 数值 |
|---|---|---|---|
| 孔隙度 | 4.8% | 基质内摩擦角 | 30° |
| 垂向主应力 | 96 MPa | 基质抗张强度 | 15 MPa |
| 水平最小主应力 | 90 MPa | 天然裂缝逼近角 | 5°~75° |
| 两向主应力差 | 12 MPa | 天然裂缝内聚力 | 5 MPa |
| 地层压力系数 | 2.0 | 天然裂缝抗张强度 | 3 MPa |
| 弹性模量 | 40 GPa | 天然裂缝内摩擦角 | 12° |
| 泊松比 | 0.23 | 起裂簇数 | 3 |
| 压裂液黏度 | 3 mPa·s | 簇间距 | 10 m |
| 压裂液排量 | 8 m3/min | 5、6井压裂时间 | 1 h |
| 基质内聚力 | 20 MPa | 停泵时间 | 1 h |
3.2 裂缝带逼近角的影响
以图9 物理模型为基础，取天然裂缝长20 m、带宽20 m，改变裂缝带逼近角（θ），模拟拉链式压裂过程中天然裂缝带逼近角对压裂裂缝扩展轨迹（见图10 ）、裂缝窜通（见图11 ）及裂缝长度（见图12 ）的影响：①图10a 显示，受天然裂缝捕获水力裂缝的影响，天然裂缝带逼近角越大，响应井水力裂缝前沿距离激动井井筒越远，说明大逼近角裂缝带对水力裂缝扩展具有阻挡作用。②对比图10a 与图10b 发现，受页岩基质超低渗透性与低滤失影响，停泵期间水力裂缝在缝内高净压力驱动下存在继续扩展的现象。③结合图10c 、图11 、图12 可知，在激动井压裂结束后，裂缝带逼近角为15°时，由于天然裂缝更容易发生张剪破坏，激动井水力裂缝在压裂0.52 h时通过天然裂缝带与响应井裂缝快速沟通，使响应井压力上涨8.3 MPa；裂缝带逼近角为60°时，激动井压裂裂缝受天然裂缝带阻挡先沿着天然裂缝带扩展，但受高水平应力差及净压力增加的影响，激动井压裂到0.94 h时，裂缝突破裂缝带阻挡与响应井压裂裂缝相交形成窜通，使响应井压力上涨5.4 MPa。④裂缝带逼近角为35°、45°时，受天然裂缝阻隔影响，水力裂缝优先沿天然裂缝扩展未形成井间窜通。由此可见，随着逼近角的增加，响应井压力涨幅呈现先减小后增加的变化趋势。
为评价压窜对压裂效果的影响，统计了压裂裂缝总长度与裂缝带逼近角的关系（见图12 ）。对比图10c 、图11 、图12 可知：①随裂缝带逼近角增加，压裂裂缝总长度呈先增加再减小的趋势；②与小逼近角天然裂缝带相比，大逼近角裂缝带窜通所需时间更长，窜通后响应井压力涨幅相对更小，窜通难度相对更大；③受井间压裂裂缝窜通的影响，压裂效率降低，窜通井裂缝总长度减小。④在高应力差下，5°~15°逼近角裂缝带比60°~75°逼近角裂缝带更易窜通，而60°~75°逼近角裂缝带又比35°~45°逼近角裂缝带更易窜通。
3.3 裂缝带宽的影响
页岩储层中天然裂缝带带宽是影响压裂裂缝窜通的重要因素。同样以图9 物理模型为基础，设置裂缝带逼近角为30°，天然裂缝长20 m，改变裂缝带带宽（wn），模拟裂缝带带宽对拉链式压裂裂缝扩展、窜扰规律的影响。图13 —图15 分别为裂缝带宽对激动井压裂停泵时刻裂缝轨迹、压裂裂缝窜通情况、压力涨幅与裂缝长度的影响。分析可知：①裂缝带带宽为5 m与15 m时，由于裂缝带带宽较小，压裂裂缝沿最大主应力方向扩展的阻力较小，两井压裂裂缝窜通，响应井压力上涨7.5 MPa与4.2 MPa；②裂缝带带宽为25 m与35 m时，带宽较大，裂缝带阻挡作用显著，两井裂缝并未窜通，响应井压力涨幅不显著；③随着裂缝带带宽的增加，拉链式压裂裂缝窜通概率相对减小，两井窜通所需要时间逐渐增加；④裂缝带带宽越大，响应井压力涨幅越小，压裂裂缝越复杂，压裂裂缝总长度更大，储层改造效率更高。
3.4 裂缝带内天然裂缝长度影响
以图9 物理模型为基础，设置裂缝带逼近角为15°，带宽20 m，改变裂缝带内天然裂缝长度，模拟裂缝带内天然裂缝长度对拉链式压裂裂缝扩展及窜扰规律的影响。图16 为响应井、激动井压裂结束时的压裂裂缝展布，图17 与图18 分别代表不同天然裂缝长度下拉链式压裂裂缝窜通情况、对应窜通时刻以及窜扰程度。分析可知，由于裂缝带逼近角较小，不同天然裂缝长度下两井压裂裂缝都形成了窜通，且窜通所需时间随天然裂缝长度增大而减少；响应井压力涨幅随天然裂缝长度增大而增大；拉链式压裂结束时形成的裂缝总长度随天然裂缝长度增大而减小。说明受小逼近角天然裂缝易激活及强捕获水力裂缝能力的影响，随着裂缝带内天然裂缝长度的增加，拉链式压裂窜通风险增大，储层改造效率降低。针对天然裂缝尺度较大、连续性较好的裂缝带，压裂施工过程中，更需要适当控制液量、排量，增加液体黏度，尽量预防或延缓裂缝窜通。
3.5 裂缝带与孔眼相对位置的影响
对于天然裂缝发育的页岩储层，水平井筒穿越裂缝带的情况难以避免，为减少井间裂缝窜通，优化设计裂缝带与孔眼的相对位置成为关键。基于图9 所示物理模型，设置裂缝带逼近角15°，裂缝带带宽20 m，天然裂缝长度10 m，定义错位距离ds表征裂缝带与孔眼相对位置（见图19 ），模拟裂缝带与孔眼相对位置对拉链式压裂裂缝扩展及窜通行为的影响。图19 为错位距离为0，5，10，20 m时对应的裂缝扩展情况。可以看出，错位距离为0，5 m时，受天然裂缝带影响，井间拉链式压裂裂缝形成了窜通；当错位距离达到10，20 m时，压裂结束时两井压裂裂缝未形成井间窜通。可见错位射孔对裂缝扩展与窜通具有重要影响，且错位距离增加，裂缝窜通概率减小。分析图20 可知，与完全正对压裂（错位距离为0）相比，错位压裂（错位距离为5 m）井间窜通所需时间更长。图21 则表明错位距离与压裂形成的裂缝总长度及响应井压力涨幅没有显著的规律性关系，虽然错位距离为0与5 m时，形成了井间裂缝窜通，但压力涨幅并非错位距离为0时最大，裂缝总长度也并非错位距离为0时最短，反而错位距离为5 m时对应的压力涨幅更大，错位距离为20 m对应的裂缝总长度最小。究其原因，错位距离为0，5 m时压裂形成的裂缝较大比例处于裂缝带内，压裂裂缝展布相对复杂，而错位距离为20 m时压裂形成的裂缝都属于基质型裂缝，压裂裂缝展布相对简单。因此，对于裂缝带发育段，考虑其特殊性，针对性设计射孔位置对于预防或减缓井间裂缝窜通具有重要意义。
4 结论
深层页岩气储层拉链式压裂裂缝扩展及窜扰规律受天然裂缝带逼近角、天然裂缝长度、裂缝带带宽及射孔位置等因素影响。大逼近角裂缝带对压裂裂缝正向扩展及井间窜通具有阻挡作用，停泵期间水力裂缝在净压力驱动下存在继续扩展行为；高应力差下，随裂缝带逼近角增加，响应井压力涨幅、压裂裂缝总长度分别呈先减小再增加及先增加再减小的趋势；与小逼近角裂缝带相比，大逼近角天然裂缝带引发窜通所需时间更长，窜通难度更大。裂缝带带宽、天然裂缝长度分别与响应井压力涨幅呈负、正相关，分别与窜通所需时间、压裂裂缝总长及压裂效率呈正、负相关；随布井错位距离增加，裂缝窜通概率减小，但错位距离与响应井压力涨幅、裂缝总长度之间的规律性不明显。
符号注释：
a——断裂节理单元的平均张开度，m；A——三角单元重叠区域面积，m2；C——阻尼矩阵，N·s/m；d——损伤因子，无因次；ds——两井射孔簇错位距离，m；det——求方阵的行列式，无因次；D——应变速率张量，s−1；E——弹性模量，Pa；Ed——格林应变张量，无因次；Es——圣维南应变张量，无因次；fn——法向接触力，Pa；fp——孔眼摩阻，Pa；fs——饱和度修正系数，%；ft——切向接触力，Pa；fw——井筒沿程摩阻，Pa；F——变形梯度矩阵，无因次；Fc——节点上的接触力矢量，N；Fd——三角单元变形诱导产生的节点力矢量，N；Fj——节点上的黏结力，N；Fext——节点上的外力矢量，N；g——重力加速度，m/s2；hc——裂缝与岩石基质之间的流体交换系数，m/(Pa·s )；i——射孔簇编号；K——渗透率张量，m2；l——两个压力节点之间的距离，m；ln——天然裂缝长度，m；L——裂缝长度，m；m——质量矩阵，N·s2/m；M——Biot模量，Pa；n，n+1——上一时间步与当前时间步，无因次；nhf——总裂缝簇数；o——法向张开量，m；op——临界法向张开量，m；or——最大法向张开量，m；p——流体压力，Pa；p1——裂缝初始端节点压力，Pa；p2——裂缝末端节点压力，Pa；pc——裂缝内流体压力，Pa；pf,in——各簇裂缝入口压力，Pa；pn——法向罚参数，Pa/m；$p_{\text{p}}^{+}$，$p_{\text{p}}^{-}$——裂隙两侧岩石基质孔隙压力，Pa；pt——切向罚参数，Pa/m；pw——各压裂段入口压力，Pa；q21——节点P1至节点P2方向上的流量，m2/s；$q_{\text{e}}^{+}$，$q_{\text{e}}^{-}$——节理单元左右两侧与岩石基质单位时间内的流体交换量，m2/s；Qi——第i簇裂缝流入流量，m3/s；Q——注入总流量，m3/s；Qtotal——单位时间内进入单元体的净流量，m3/s；s——切向滑移量，m；sp——临界剪切滑移量，m；sr——最大剪切滑移量，m；t——时间，s；T——应力张量，Pa；V——任意给定单元体体积，m3；v——单位截面通过的流体流速，m/s；wn——裂缝带带宽，m；x，y——直角坐标系，m；x——节点坐标的向量，m；$\overline{x}$——节点速度向量，m/s；$\overline{\overline{x}}$——节点加速度向量，m/s2；y1，y2——节点1和2的y轴坐标，m；α——Biot系数，无因次；βt，βc——靶体，目标体；θ——逼近角，(°)；ν——泊松比，无因次；λ——黏滞阻尼系数，Pa·s；ρw——流体密度，kg/m3；μ——流体黏度，Pa·s；σ——抗拉强度，Pa；σH，σh——水平最大、最小主应力，Pa；σt——法向应力，Pa；τ——抗剪强度，Pa；τs——切向应力，Pa；Δp——节点P1与P2间的总压差，Pa；Δut——某一时刻接触点的相对位移，m；${{\varphi }_{\text{t}}}\left( {{P}_{\text{t}}} \right)$，${{\varphi }_{\text{c}}}\left( {{P}_{\text{c}}} \right)$——重叠区域内Pt、Pc处的势，无因次；${{\beta }_{\text{t}}}\cap {{\beta }_{\text{c}}}$——靶体βt和目标体βc重叠区域的边界，无因次。