页岩气藏水平井分段压裂缝间应力干扰全三维模拟

唐慧莹，张东旭，刘环竭,梁海鹏，张烈辉

(1.西南石油大学“油气藏地质及开发工程”国家重点实验室，四川 成都610500；2.中石油西南油气田 开发事业部，四川 成都 610017)

引 言

由于页岩等非常规储层基质渗透率极低，减少裂缝间距能够增大井筒与储层的接触面积，提高后期的油气产量。但是裂缝间距不能无限制减小。一方面，当裂缝间距足够小时，增大裂缝条数对生产影响不明显；另一方面，裂缝间距减小，裂缝间相互干扰增强、相互挤压，造成部分裂缝欠发育或者不发育[1-2]。

近十年来，为了定量分析多裂缝间的相互作用，考虑裂缝间相互影响的多裂缝扩展压裂模型研究不断增加。Cheng等[3]采用二维位移不连续方法对比了不同裂缝间距下裂缝开度的差异，但该模型并未考虑裂缝高度的影响。Olson等[4]引入高度修正公式，计算了多条裂缝以亚临界速度扩展时的裂缝形态。Wu等[5]通过对Olson高度修正公式的进一步改进，得到了更为准确的基于二维位移不连续方法的多裂缝应力干扰计算结果。但地层条件下裂缝均以三维形态存在，二维模型在刻画裂缝高度生长以及复杂裂缝形态等方面仍有欠缺。Tang等[6]构建了基于平面三维位移不连续方法的多裂缝扩展模型，计算结果表明纵向地应力分布也会影响裂缝间相互干扰的强弱。Peirce等[7]采用平面三维有限元与水平集方法优化了同压裂段内各射孔簇的排布方式。Cheng等[8]通过能量的观点推导了多条径向型裂缝同时扩展时，裂缝半径与时间的解析解。平面三维模型是对二维模型的一大改进，平面三维模型具备描述整个裂缝平面生长的能力，对裂缝高度生长的刻画及对裂缝内流体流动过程的描述都明显优于二维模型。但是，当初始裂缝面不与最小地应力垂直或地应力差较弱时，裂缝平面有可能发生明显偏转。在这类情况下，平面三维模型不再适用，此时需要采用全三维的压裂模型。

Castonguay等[9]采用四边形网格刻画了压裂过程中在不同应力状态下多裂缝的扩展形态。Kumar等[10]同样使用四边形网格模拟了压裂过程中多裂缝同时扩展等物理过程。Maerten等[11]采用三角形网格剖分裂缝，模拟地层中断层等不连续面在外加载荷下的变形，着重对III型撕拉型裂缝的扩展进行了研究。上述压裂模型均采用位移不连续方法求解全三维裂缝的扩展问题，这是由于边界元方法只需对裂缝面进行求解，与求解全空间的有限元等方法相比，一方面避免了剖分三维网格的困难，另一方面能够有效提高计算效率。虽然已有部分基于位移不连续方法的全三维压裂模型，但文献中并未与二维计算结果进行对比，也未对多裂缝扩展时三维空间应力状态的变化进行分析。

本文首先介绍基于位移不连续方法的全三维多裂缝扩展模型，重点介绍了近似奇异、超奇异影响系数的数值计算，裂缝前沿单元的加密与合并，自适应积分点选取方式等。在对模型计算结果进行验证后，先后对二维与三维模型的多裂缝应力干扰强度、诱导应力分布、多裂缝扩展过程进行了对比。

1 计算方法

1.1 基于三角形单元的三维位移不连续方法

三维与二维压裂模型的显著区别在于裂缝维度的增加。在边界元方法中，二维空间内裂缝可以用线段表征，在三维空间内，裂缝为二维平面，需要用相应的二维网格单元进行刻画。假设裂缝平面由一系列三角形单元网格构成，每个长方形单元有3个方向的位移不连续量：与裂缝面平行且相互正交的切向位移不连续量Dx与Dy以及垂直于裂缝平面的法向位移不连续量Dz。位移不连续量在单元局部坐标系(x-y-z)定义为：

Dx=ux(x,y,0-)-ux(x,y,0+)，

Dy=uy(x,y,0-)-uy(x,y,0+)，

Dz=uz(x,y,0-)-uz(x,y,0+)。

(1)

裂缝单元单位位移不连续量在空间任一点(x,y,z)诱导的位移及应力大小为：

ux=C{[2(1-v)I,z-zI,xx]Dx-zI,xyDy-[(1-2v)I,x+zI,xz]Dz}，

uy=C{-zI,xyDx+[2(1-v)I,z-zI,yy]Dy-[(1-2v)I,y+zI,yz]Dz}，

uz=C{[(1-2v)I,x-zI,xz]Dx+[(1-2v)I,y-zI.yz]Dy+[2(1-v)I,z-zI,zz]Dz}；

σxx=2C{[2I,xz-zI,xxx]Dx+[2vI,yz-zI,xxy]Dy+[I,zz+(1-2v)I,yy-zI,xxz]Dz}，

σyy=2C{[2vI,xz-zI,xyy]Dx+[2I,yz-zI,yyy]Dy+[I,zz+(1-2v)I,xx-zI,yyz]Dz}，

σzz=2C{-zI,xzzDx-zI,yzzDy+[I,zz-zI,zzz]Dz}；

τxy=2C{[(1-v)I,yz-zI,xxy]Dx+[(1-v)I,xz-zI,xyy]Dy-[(1-2v)I,xy+zI,xyz]Dz}，

τyz=2C{-[vI,xy+zI,xyz]Dx+[I,zz+vI,xx-zI,yyz]Dy-zI,yzzDz}，

τzx=2C{[I,zz+vI,yy-zI,xxz]Dx-[vI,xy+zI,xyz]Dy-zI,xzzDz} 。

(2)

图1 加密与合并操作Fig.1 Densifying and merging

式(2)中：C为与岩石力学性质相关的常数，

(3)

v为岩石基质泊松比；函数I为半无穷空间格林函数，

(4)

式(4)中S为当前三角形单元面积。假设单元i满足应力边界条件，所受法向应力为σzi，沿局部坐标系x、y方向的切向应力分别为τxi和τyi。此时所有裂缝单元作用在单元i的应力总和加上就地应力，应等于作用在单元表面的应力边界条件：

(5)

1.2 三角形单元影响系数计算方法

三维位移不连续方法的一大难点在于影响系数的计算。三角形单元存在特殊的解析求解方法[12-13]，但推导过程及表达式都比较繁琐，同时对于某些位置点的应力计算无效， 因此对于不规则几何单元，更为便捷的方式是寻求影响系数的数值解。

当函数类型非奇异时，数值积分可以取得较好的准确度，但对于位移不连续方法，裂缝单元对自身的积分不只存在奇异性，还存在超奇异性。目前已有多种求解奇异积分的方法[14]，本文采用基于有限部分的Hadamard积分[15]与三角形分解方法进行计算，该方法表达形式简洁，同时有较高的计算精度。本文同时采用Li等[16]提出的平方根形式尖端单元进一步提高模型计算精度。

三角形位移不连续单元数值积分的精度与待计算点位置、单元形态、积分点数量与位置的选取都有密切关系。待计算点距离单元越近，则所需积分点越多。因此，本文引入接近度这一概念，用于刻画待计算点与单元的距离

e=D/Lmax。

(6)

式(6)中，D为待计算点与单元中心距离，Lmax为三角形单元最长边长度。当接近度e<1.5时，本文模型将选取更多的高斯积分点(本文中选取9，其余位置选择6)。同时，本文模型对裂缝前沿单元采用了网格加密与合并(图1)。当相邻前沿节点距离较远时，采用加密操作，在相邻节点中增加新的单元节点。反之，当相邻前沿节点过于接近时，对这对节点进行合并，作为一个节点。图2(a)是采用加密与合并操作后，PKN型裂缝扩展过程的模拟结果(11个生长步)， 图2(b)是未对前沿节点进行调整时得到的裂缝形态(5个生长步)。可以看出，当裂缝在各个方向的生长不均时，若不对前沿节点进行调整，将会出现畸形网格，进而影响模型的计算精度。

图2 加密与合并操作前后PKN型裂缝扩展模拟结果Fig.2 Simulation results of PKN-type crack propagation before and after densifying and merging

全三维裂缝扩展时，可采用Schöllmann等[17]的最大主应力准则，该准则认为裂缝沿主应力最大的平面扩展。本文采用的网格延伸方式与Meng 等[18]的相似。以尖端单元的中点为起点，根据计算得到的扩展长度Δd和扩展方向θ0，得到向量终点，再将所有向量终点连接起来，得到新的裂缝边界。

2 计算结果与分析

2.1 三维裂缝应力干扰计算

计算不同三维裂缝形态下多裂缝间的应力干扰强度，并将基于三维位移不连续方法的计算结果与二维计算结果进行对比。模型输入基本参数见表1。

表1 算例基本参数Tab.1 Basic parameters of an example

对比了半径为50 m的径向裂缝以及长度和高度分别为100 m、20 m，100 m、40 m，200 m、40 m的PKN型裂缝，在不同裂缝间距下外侧裂缝与内侧裂缝的开度差异。图3为上述4种裂缝在裂缝间距为30 m时的裂缝开度分布。从图中可以看出，径向裂缝间应力干扰最为明显。对于PKN型裂缝，增加裂缝高度和裂缝长度会增强裂缝间的挤压作用，裂缝间应力干扰强度对高度的敏感性强于对长度的敏感性。

图3 不同裂缝形态在裂缝间距为30 m时的开度分布Fig.3 Opening distribution of different fracture patterns at fracture spacing of 30 m

图4中展示了采用纯二维模型与Olson的高度修正二维模型3个裂缝高度对应的模拟结果。对比图4与图3看到，二维模型计算得到的外侧裂缝开度更大，内侧裂缝开度更小。从图5的对比图中可以更加清晰地看到二维与三维模型计算结果的差异，纯二维模型计算得到的裂缝间挤压作用显著(黑色虚线)。PKN型裂缝间挤压作用的强弱与裂缝高度有更强的相关性，即裂缝间应力干扰更易受裂缝中较短边控制。

通过以上分析可知二维压裂模型预测内外侧裂缝开度差异大于三维模型，可以推测二维模型预测的裂缝诱导地应力变化应强于三维模型。图6是3条PKN型裂缝(长度100 m，高度30 m，间距30 m)并排分布时二维应力计算结果(定义拉应力为正)，每条裂缝内流体净压力为2 MPa，图中显示应力仅为裂缝变形诱导应力。由于压裂缝沿Y轴生长，因此在该例子中X方向为最小主应力方向。图6中，裂缝内侧空间由于受到挤压作用呈现较大的压应力。同时，裂缝内侧部分的应力差减小，地应力差减小，压裂后裂缝更有可能形成复杂缝网，同时已有天然裂缝也更容易被激活[19]。

图4 二维模型与裂缝的高度修正二维模型(裂缝高度为20、40、60 m)的模拟结果Fig.4 Simulation results by using two-dimensional model and modified two-dimensional model(crack height 20,40,60 m)

图5 不同裂缝间距下外侧与内侧裂缝中心开度比值Fig.5 Center opening ratio of lateral crack and medial crack under different crack spacing

图6 二维模型应力计算结果Fig.6 Stress calculated using 2D model

图7—图9对应半径50 m、间距30 m的3条径向裂缝变形诱导的应力场，裂缝内净压力同样为2 MPa。从图中可以看出，三维径向裂缝诱导压应力小于二维裂缝场。三维空间应力差值DS(Y方向应力与X方向应力之差)小于二维裂缝模拟结果，即三维模型计算得到的应力差减小程度小于二维模型。

图7 3条径向裂缝诱导不同截面X方向应力分布图Fig.7 X-direction stress distribution map of different sections induced by three radial cracks

图8 3条径向裂缝诱导不同截面Y方向应力分布图Fig.8 Y-direction stress distribution map of different sections induced by three radial cracks

图9 3条径向裂缝诱导不同截面Y方向应力与X方向应力差分布图Fig.9 Y-X direction dtress difference distribution map of different sections induced by three radial cracks

图10对比了3条径向裂缝与3条PKN型裂缝诱导的应力差分布。从图中可以看到，与径向裂缝相比，PKN型裂缝诱导地应力差减小范围更小，且减小幅度更小。

图10 径向裂缝和PKN型裂缝地应力差等值云图Fig.10 Equivalent nephograms of in-situ stress difference of radial fracture and PKN fracture

2.2 三维裂缝扩展过程模拟

以上结果均为静态裂缝分析，以下将对裂缝动态扩展过程进行比较。模拟使用参数见表2。假设地层应力分布均匀，裂缝条数为3，裂缝初始半径为10 m，裂缝间距分别为20、30、40 m，扩展15个生长步后裂缝形态如图11所示。可以看到，随着裂缝间距的增加，内侧裂缝与外侧裂缝面积与开度的差异逐渐减小。同时，由于裂缝中心部分所受挤压较大，外侧裂缝外沿开度反而大于中心部位开度。

表2 裂缝动态扩展使用参数Tab.2 Parameters for crack propagation

图12为裂缝高度为20 m与40 m、裂缝初始长度20 m时采用高度修正二维模型的模拟结果。可以看出，当假设裂缝高度为20 m时，二维模型得到的裂缝绝对开度与内外侧裂缝开度差异小于三维模型(图1)，而当裂缝高度为30 m时，二维模型计算结果得到的外侧裂缝开度更大，同时外侧裂缝对内侧裂缝挤压更为明显。同时，二维模型计算得到的裂缝长度总大于三维径向裂缝长度。

图11 三条裂缝在间距为20 m 、30 m 、40 m时扩展15个生长步后裂缝开度分布Fig.11 Opening distribution after three cracks propagating 15 growth steps under fracture spacing of 20 m,30 m and 40 m

图12 裂缝高度为20 m和30 m，间距为20 m、30 m、40 m时二维模型模拟结果Fig.12 Simulation results of two-dimensional model under crack height of 20 m and 30 m and crack spacing of 20 m,30 m and 40 m

3 结 论

(1)PKN型裂缝间应力干扰强度与裂缝高度、裂缝长度均有关系，应力干扰对短边长度的敏感度更高。径向裂缝间的应力干扰强度强于PKN型裂缝。

(2)采用高度未修正的二维模型预测的裂缝间应力干扰最强，高度修正二维模型预测的裂缝间应力干扰强度次之，三维模型预测的裂缝间应力干扰最弱。

(3)用二维模型预测的裂缝变形诱导应力变化大于三维模型，基于二维模型计算得到的裂缝间应力差减小程度亦大于三维模型，即根据二维模型结果判断的应力转向区域将大于三维模型的结果。径向裂缝诱导应力变化大于PKN型裂缝。

(4)用二维模型预测的裂缝扩展速度大于三维模型。即使采用高度修正的二维模型，也很难获得与三维模型相似的结果，尤其难以刻画径向裂缝的生长规律。

(5)用二维模型预测裂缝间应力干扰的准确度不高。为了更准确地刻画裂缝间应力干扰，达到优化裂缝间距及压裂方案的目的，仍需要采用三维模型。