CN113158425B - 一种基于边界元方法的全三维裂缝相交过程模拟方法 - Google Patents

Abstract

本发明涉及裂缝模拟技术领域，具体地说，涉及一种基于边界元方法的全三维裂缝相交过程模拟方法，其包括以下步骤：(1)时间步与注入压力初始化；(2)、裂缝开度计算；(3)、质量守恒检查；(4)、判断是否满足裂缝扩展条件；(5)、裂缝扩展计算；(6)、自适应调整天然裂缝网格；(7)、自适应调整人工裂缝网格；(8)、进入下一个时间步的计算。本发明可实现人工裂缝与天然裂缝相交过程中，人工与天然裂缝内网格的动态、自适应调整，进而为利用边界元方法进行三维裂缝相交模拟提供重要的技术支撑，推动全三维裂缝扩展模拟的进一步发展。

Description

一种基于边界元方法的全三维裂缝相交过程模拟方法

技术领域

本发明涉及裂缝模拟技术领域，具体地说，涉及一种基于边界元方法的全三维裂缝相交过程模拟方法。

背景技术

利用边界元方法进行裂缝扩展模拟具有待求解未知数少、计算精度高等优点，适用于模拟大型水力压裂过程。但是目前基于边界元方法的全三维裂缝扩展模型尚未考虑裂缝相交的情形，导致对水力压裂缝网真实形态的刻画不足，影响了压裂施工方案的制定与优化。网格问题是阻碍全三维压裂缝网模型发展的一个主要原因。三维裂缝平面相交时，网格处理非常复杂，尤其是当裂缝动态相交时，需要实时对网格进行调整和更新，既要保证网格质量又要避免过度的网格加密带来的计算负担。现有处理裂缝相交网格的技术均是针对裂缝面固定的情形，无法处理裂缝扩展过程中裂缝形态的动态变化。同时，三维裂缝面相交过程还需对人工裂缝扩展网格进行特殊处理，以实现对相交过程的完整刻画。

发明内容

本发明的内容是提供一种基于边界元方法的全三维裂缝相交过程模拟方法，其能够克服现有技术的某种或某些缺陷。

根据本发明的一种基于边界元方法的全三维裂缝相交过程模拟方法，其包括以下步骤：

(1)时间步与注入压力初始化；

将上个时间步n的时间步长Δtⁿ和注入点压力

作为当前时间步n+1的初始猜测，下标inj表示注入；

(2)、裂缝开度计算；

使所有裂缝单元压力均与注入点压力相等，采用边界元算法计算当前时间步长及注入点压力下裂缝开度

下标i为裂缝单元编号，k为更新注入点压力的的迭代步编号，j为更新时间步长的迭代步编号；

(3)、质量守恒检查；

检查当前压力迭代步对应注入压力下，当前时间步长内注液总量与裂缝体积总变化量是否相等，若相等则进入下一步，否则采用二分法对注入点压力进行更新，并返回步骤(2)；

(4)、判断是否满足裂缝扩展条件；

计算当前时间步长下，人工裂缝尖端的最大等效应力集中因子

与天然裂缝尖端单元的最大F值

下标HF与NF分别代表人工裂缝与天然裂缝；检查注入当前时间步长流体后，人工裂缝或天然裂缝是否能发生拟稳态扩展，若不能，则采用二分法对时间步长进行更新，并返回步骤(2)，直到人工裂缝或天然裂缝满足拟稳态扩展条件，则进入下一步；

(5)、裂缝扩展计算；

判断裂缝是否扩展及计算裂缝扩展长度；

(6)、自适应调整天然裂缝网格；

采用逐次节点捕获技术对于人工裂缝相交后的天然裂缝平面网格进行自适应调整与优化；

如果人工裂缝相邻前沿网格点均恰落在天然裂缝平面的边界位置；

初始时刻，天然裂缝网格并未沿着交线分布；

天然裂缝网格点的调整分为两类：a.移动网格点与交线端点重合；b.调整网格点使之落在交线内部，对于上述两类情形，天然裂缝网格调整均需经历以下3个步骤：

a.1、找到网格平面内的捕获点；

a.2、寻找与捕获点相邻的图距离为N_dis的节点，定义为松弛节点；

a.3、等间距移动捕获点至目标点位置；在松弛过程中，采用优化算法对被松弛网格总体质量进行优化；如果总体网格质量达不到要求，则增加松弛节点的搜索图距离；

重复a.2与a.3，直至达到最大图距离N_relax；

对于a类网格调整，a.1中的捕获点为与端点最近的节点，a.3中的目标点为交线端点；

对于b类网格，捕获点搜索方法如下：

网格平面空间被交线上最后一个节点与交线分为三个区域，三个区域内分别对应S＝1、0、-1；与交线上最后一个节点相连的单元若有边界跨过交线，则从该边界上的两个节点中寻找捕获点；将与交线较为靠近的节点选为捕获点；目标点为捕获点在交线上的投影；

若目标位置过于靠近交线端点或落在交线外侧，则按照a类网格调整方法，移动捕获点与端点重合；

(7)、自适应调整人工裂缝网格；

(7.1)、对人工裂缝扩展前沿进行平滑处理，通过3阶Savitzky-Golay滤波器插值方法，将裂缝扩展向量转移到裂缝前沿网格点上，裂缝扩展向量与裂缝前沿网格点所在切线相垂直；

(7.2)、构建新的裂缝单元；

当人工裂缝与天然裂缝相交时，需重构多个裂缝单元；

当相邻两人工裂缝前沿扩展向量均为不0时，主要由较长扩展向量对应的前沿点与交线上天然裂缝的网格点连接，构成新的人工裂缝单元；

若有一侧人工裂缝前沿已与天然裂缝相交，则所有新增单元均与未相交的人工裂缝前沿网格点连接；

当人工裂缝与天然裂缝边界相交时，需对人工裂缝扩展向量方向进行调整，若相邻网格点中有一个网格点的扩展向量相交与天然裂缝平面外时，需要对相交与平面外的扩展向量或相交于平面内的扩展向量进行移动；使新增网格统一处于平面内或平面外；

(8)、进入下一个时间步的计算。

作为优选，步骤(2)中，裂缝开度计算的方法为：

设裂缝平面由一系列三角形单元网格构成，每个长方形单元有3个方向的位移不连续量：

与裂缝面平行且相互正交的切向位移不连续量D_x与D_y以及垂直于裂缝平面的法向位移不连续量D_z；

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

u_x＝C{[2(1-v)I_,z-zI_,xx]D_x-zI_,xyD_y-[(1-2v)I_,x+zI_,xz]D_z}

u_y＝C{-zI_,xyD_x+[2(1-v)I_,z-zI_,yy]D_y-[(1-2v)I_,y+zI_,yz]D_z}

u_x＝C{[(1-2v)I_,x-zI_,xz]D_x+[(1-2v)I_,y-zI_.yz]D_y+[2(1-v)I_,z-zI_,zz]D_z}

σ_xx＝2C{[2I_,xz-zI_,xxx]D_x+[2vI_,yz-zI_,xxy]D_y+[I_,zz+(1-2v)I_,yy-zI_,xxz]D_z}

σ_yy＝2C{[2vI_,xz-zI_,xyy]D_x+[2I_,yz-zI_,yyy]D_y+[I_,zz+(1-2v)I_,xx-zI_,yyz]D_z}

σ_zz＝2C{-zI_,xzzD_x-zI_,yzzD_y+[I_,zz-zI_,zzz]D_z}

τ_xy＝2C{[(1-v)I_,yz-zI_,xxy]D_x+[(1-v)I_,xz-zI_,xyy]D_y-[(1-2v)I_,xy+zI_,xyz]D_z}

τ_yz＝2C{-[vI_,xy+zI_,xyz]D_x+[I_,zz+vI_,xx-zI_,yyz]D_y-zI_,yzzD_z}

τ_zx＝2C{[I_,zz+vI_,yy-zI_,xxz]D_x-[vI_,xy+zI_,xyz]D_y-zI_,xzzD_z}；

其中，C为与岩石力学性质相关的常数：

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

S为当前三角形单元面积；

设单元i满足应力边界条件，所受法向应力为σ_z ⁱ，对应裂缝内压强大小，沿局部坐标系x、y方向的切向应力分别为τ_x ⁱ和τ_y ⁱ；

此时所有裂缝单元作用在单元i的应力总和加上就地应力，应等于作用在单元表面的应力边界条件：

上式中A为影响系数；

裂缝单元i的开度w_i对应单元z方向位移不连续量

的负数，即

作为优选，步骤(4)中，人工裂缝尖端的最大等效应力集中因子

与天然裂缝尖端单元的最大F值

分别计算如下：

采用最大主应力原理判断人工裂缝是否发生扩展，最大主应力计算公式为：

式中σ_θ、σ_z、τ_θz分别为径向坐标系下的轴向应力、垂向应力与切向应力；

σ_θ、σ_z、τ_θz可基于I型、II型与III型断裂应力集中因子K_I，K_II和K_III进行计算：

式中r与θ分别为裂缝尖端径向坐标系下的轴向距离与上倾角；应力集中因子可根据位移不连续量D_i计算：

式中E和v分别为储层岩石的杨氏模量与泊松比；使得最大主应力取最大值的倾角θ₀为裂缝扩展方向，需要满足下列条件：

倾角θ₀下裂缝尖端等效应力集中因子K_eq用于判断裂缝是否扩展及相应扩展距离：

对于天然裂缝尖端单元，采用F判据判断天然裂缝单元是否激活：

对于天然裂缝面内的裂缝单元激活有θ＝0⁰；式中K_IC与K_IIC分别为天然裂缝的I型与II型断裂的临界应力集中因子；当F值大于1时，天然裂缝单元被激活。

作为优选，步骤(4)中，判断人工裂缝或天然裂缝是否能发生拟稳态扩展的方法为：

最大等效应力集中因子

或天然裂缝尖端单元的最大F值

在该时间步内取值在区间[(1-α,1+α)]内，其中α为接近0的小数。

作为优选，步骤(5)中人工裂缝或天然裂缝扩展距离及是否扩展的判断依据如下：

式中下标HF与NF分别代表人工裂缝与天然裂缝；d_max是单步最长扩展距离。

作为优选，步骤(6)，a.3步骤中的网格优化采用Matlab自带的带约束非线性优化函数fmincon实现，对于位于平面内部的节点，约束为移动后的节点不得脱离裂缝平面；对于位于平面边界的节点，约束条件为移动后的节点需仍位于平面边界处；优化函数为：

下标R代表被松弛的节点或网格组合，v代表节点，κ代表网格；B是满足所有约束条件的点的集合，网格质量Q定义如下：

式中

δ等于0.3，A是网格面积，L_j,j＝1,2,3是网格κ的边长度；Q越大，网格越接近等边三角形；当所有被松弛网格中有网格质量Q小于0.4，则认为网格质量未达到标准，需扩大搜索捕获点的图距离，重新更新节点位置。

作为优选，步骤(7.2)中，采用下述3类方法对新增人工裂缝网格形态进行调整：劈分、合并和扩展，上述操作进行的前提条件如下：

式中Δd为扩展向量的长度，L_new为新单元相邻前沿点的距离，当扩展后扩展向量大于网格边界的2倍，即再补充进行劈分操作。

作为优选，步骤(7.2)中，当人工裂缝越过天然裂缝后，天然裂缝的“背面”再次与天然裂缝相交时，若交点与原始交线较为接近，则在原始交线上寻找与交点最近的点，并将扩展向量移动至该处；若交点与原始交线距离较远，则需构建一条新的交线。

本发明可实现人工裂缝与天然裂缝相交过程中，天然裂缝内网格的动态、自适应调整，该方法不需要进行网格加密，只需在原有天然裂缝网格基础上进行网格点坐标的微小调整，既可以有效实现相交平面网格内的网格调整，又可以在不进行网格加密的前提下生成质量较好的新网格系统。同时提出相交过程中人工裂缝网格的自适应调整方案，进而为利用边界元方法进行三维裂缝相交模拟提供重要的技术支撑，推动全三维裂缝扩展模拟的进一步发展。

附图说明

图1为实施例1中一种基于边界元方法的全三维裂缝相交过程模拟方法的流程图；

图2为实施例1中三维边界元裂缝单元示意图；

图3为实施例1中扩展向量及扩展前沿计算方法示意图；

图4为实施例1中对新增网格形态进行调整的示意图；

图5为实施例1中相邻人工裂缝扩展向量均不为0和一侧人工裂缝前沿已与天然裂缝相交的示意图；

图6为实施例1中交叉点移动示意图；

图7为实施例1中越过天然裂缝后人工裂缝扩展向量调整的示意图；

图8为实施例1中人工裂缝与天然裂缝相交过程中天然裂缝网格单元调整过程示意图；

图9为实施例1中天然裂缝网格自适应调整过程示意图；

图10为实施例1中交线内部捕获点的搜索方法示意图；

图11为实施例1中0～2.5s不同时间裂缝开度分布示意图；

图12为实施例1中5.9～30s不同时间裂缝开度分布示意图。

具体实施方式

为进一步了解本发明的内容，结合附图和实施例对本发明作详细描述。应当理解的是，实施例仅仅是对本发明进行解释而并非限定。

实施例1

本实施例提供了一种基于边界元方法的全三维裂缝相交过程模拟方法，模型假设边界条件为定注液流量边界，其包括以下步骤：

(1)时间步与注入压力初始化；将上个时间步n的时间步长Δtⁿ和注入点压力

作为当前时间步n+1的初始猜测，下标inj表示注入；

(2)、裂缝开度计算；使所有裂缝单元压力均与注入点压力相等，采用边界元算法计算当前时间步长及注入点压力下裂缝开度

下标i为裂缝单元编号，k为更新注入点压力的的迭代步编号，j为更新时间步长的迭代步编号；

裂缝开度计算的方法为：在边界元方法中，二维空间内裂缝可以用线段表征，在三维空间内，裂缝为二维平面，需要用相应的二维网格单元进行刻画。假设裂缝平面由一系列三角形单元网格构成，每个长方形单元有3个方向的位移不连续量：与裂缝面平行且相互正交的切向位移不连续量D_x与D_y以及垂直于裂缝平面的法向位移不连续量D_z(图2)；

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

u_x＝C{[2(1-v)I_,z-zI_,xx]D_x-zI_,xyD_y-[(1-2v)I_,x+zI_,xz]D_z}

u_y＝C{-zI_,xyD_x+[2(1-v)I_,z-zI_,yy]D_y-[(1-2v)I_,y+zI_,yz]D_z}

u_x＝C{[(1-2v)I_,x-zI_,xz]D_x+[(1-2v)I_,y-zI_.yz]D_y+[2(1-v)I_,z-zI_,zz]D_z}

σ_xx＝2C{[2I_,xz-zI_,xxx]D_x+[2vI_,yz-zI_,xxy]D_y+[I_,zz+(1-2v)I_,yy-zI_,xxz]D_z}

σ_yy＝2C{[2vI_,xz-zI_,xyy]D_x+[2I_,yz-zI_,yyy]D_y+[I_,zz+(1-2v)I_,xx-zI_,yyz]D_z}

σ_zz＝2C{-zI_,xzzD_x-zI_,yzzD_y+[I_,zz-zI_,zzz]D_z}

τ_xy＝2C{[(1-v)I_,yz-zI_,xxy]D_x+[(1-v)I_,xz-zI_,xyy]D_y-[(1-2v)I_,xy+zI_,xyz]D_z}

τ_yz＝2C{-[vI_,xy+zI_,xyz]D_x+[I_,zz+vI_,xx-zI_,yyz]D_y-zI_,yzzD_z}

τ_zx＝2C{[I_,zz+vI_,yy-zI_,xxz]D_x-[vI_,xy+zI_,xyz]D_y-zI_,xzzD_z}；

其中，C为与岩石力学性质相关的常数：

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

S为当前三角形单元面积；设单元i满足应力边界条件，所受法向应力为σ_z ⁱ，对应裂缝内压强大小，沿局部坐标系x、y方向的切向应力分别为τ_x ⁱ和τ_y ⁱ；此时所有裂缝单元(假设总单元数为N)作用在单元i的应力总和加上就地应力(下标为0)，应等于作用在单元表面的应力边界条件(缝内流体压力与剪切应力)：

上式中A为影响系数；裂缝单元i的开度w_i对应单元z方向位移不连续量

的负数，即

本实施例采用基于有限部分的Hadamard积分与三角形分解方法进行计算，该方法表达形式简洁，同时有较高的计算精度。本文同时采用平方根形式尖端单元进一步提高模型计算精度。

(3)、质量守恒检查；

检查当前压力迭代步对应注入压力下，当前时间步长内注液总量与裂缝体积总变化量是否相等，若相等则进入下一步，否则采用二分法对注入点压力进行更新，并返回步骤(2)；

(4)、判断是否满足裂缝扩展条件；

计算当前时间步长下，人工裂缝尖端的最大等效应力集中因子

与天然裂缝尖端单元的最大F值

下标HF与NF分别代表人工裂缝与天然裂缝；检查注入当前时间步长流体后，人工裂缝或天然裂缝是否能发生拟稳态扩展，若不能，则采用二分法对时间步长进行更新，并返回步骤(2)，直到人工裂缝或天然裂缝满足拟稳态扩展条件，则进入下一步；

人工裂缝尖端的最大等效应力集中因子

与天然裂缝尖端单元的最大F值

分别计算如下：

采用最大主应力原理判断人工裂缝是否发生扩展，最大主应力计算公式为：

式中σθ、σ_z、τ_θz分别为径向坐标系下的轴向应力、垂向应力与切向应力；σθ、σ_z、τ_θz可基于I型、II型与III型断裂应力集中因子K_I，K_II和K_III进行计算：

式中r与θ分别为裂缝尖端径向坐标系下的轴向距离与上倾角；应力集中因子可根据位移不连续量D_i计算：

式中E和v分别为储层岩石的杨氏模量与泊松比；使得最大主应力取最大值的倾角θ₀为裂缝扩展方向，需要满足下列条件：

倾角θ₀下裂缝尖端等效应力集中因子K_eq用于判断裂缝是否扩展及相应扩展距离：

对于天然裂缝尖端单元，采用F判据判断天然裂缝单元是否激活：

对于天然裂缝面内的裂缝单元激活有θ＝0⁰；式中K_IC与K_IIC分别为天然裂缝的I型与II型断裂的临界应力集中因子；当F值大于1时，天然裂缝单元被激活。

判断人工裂缝或天然裂缝是否能发生拟稳态扩展的方法为：

最大等效应力集中因子

或天然裂缝尖端单元的最大F值

在该时间步内取值在区间[(1-α,1+α)]内，其中α为接近0的小数，在模型中取0.1。

(5)、裂缝扩展计算；

判断裂缝是否扩展及计算裂缝扩展长度；

人工裂缝或天然裂缝扩展距离及是否扩展的判断依据如下：

式中下标HF与NF分别代表人工裂缝与天然裂缝；d_max是单步最长扩展距离。

(6)、自适应调整天然裂缝网格；

采用逐次节点捕获(successive node snapping，简称SNS)技术对于人工裂缝相交后的天然裂缝平面网格进行自适应调整与优化；SNS处理该问题有以下优点:(1)网格节点的数量或连通性不需要做任何改变；(2)只有曲线周边少数节点需要移动；(3)可以考虑多条曲线存在的情形。

在本实施例涉及问题中，曲线为相邻人工裂缝扩展向量与天然裂缝的交点之间的连线，交点为曲线的端点。人工裂缝与天然裂缝相交过程中天然裂缝网格单元调整过程如图8所示。图8中，(a)初次相交；(b)初次相交后网格调整结果；(c)再次相交；(d)再次相交后网格调整结果。

如果人工裂缝相邻前沿网格点均恰落在天然裂缝平面的边界位置(图9)；初始时刻，天然裂缝网格并未沿着交线分布(图9(a))；天然裂缝网格点的调整分为两类：a.移动网格点与交线端点重合；b.调整网格点使之落在交线内部，对于上述两类情形，天然裂缝网格调整均需经历以下3个步骤：

a.1、找到网格平面内的捕获点；

a.2、寻找与捕获点相邻的图距离为N_dis的节点(图距离N_dis不大于预设的最远图距离N_relax)，定义为松弛节点；

a.3、等间距移动捕获点至目标点位置；在松弛过程中，采用优化算法对被松弛网格总体质量进行优化；如果总体网格质量达不到要求，则增加松弛节点的搜索图距离，重复a.2与a.3，直至达到最大图距离N_relax；

对于a类网格调整，a.1中的捕获点为与端点最近的节点，a.3中的目标点为交线端点；

对于b类网格，b类网格调整与a类网格调整接近，区别主要在于捕获点的搜索方法。图10展示了如何搜索交线内部的捕获点。捕获点搜索方法如下：网格平面空间被交线上最后一个节点与交线分为三个区域，三个区域内分别对应S＝1、0、-1；与交线上最后一个节点相连的单元若有边界跨过交线，则从该边界上的两个节点中寻找捕获点；将与交线较为靠近的节点选为捕获点；目标点为捕获点在交线上的投影；若目标位置过于靠近交线端点或落在交线外侧，则按照a类网格调整方法，移动捕获点与端点重合；图10中，单元边界节点分别位于(a)S＝-1与S＝1区域或(b)S＝0与S＝1区域。

a.3步骤中的网格优化采用Matlab自带的带约束非线性优化函数fmincon实现，对于位于平面内部的节点，约束为移动后的节点不得脱离裂缝平面；对于位于平面边界的节点，约束条件为移动后的节点需仍位于平面边界处；优化函数为：

下标R代表被松弛的节点或网格组合，v代表节点，κ代表网格；B是满足所有约束条件的点的集合，网格质量Q定义如下：

式中

δ等于0.3，A是网格面积，L_j,j＝1,2,3是网格κ的边长度；Q越大，网格越接近等边三角形；当所有被松弛网格中有网格质量Q小于0.4，则认为网格质量未达到标准，需扩大搜索捕获点的图距离，重新更新节点位置。

(7)、自适应调整人工裂缝网格；

(7.1)、首先根据前述算法计算得到沿裂缝边界中点先出发的裂缝扩展向量的长度与方向(图3中右方箭头)，对人工裂缝扩展前沿进行平滑处理，通过3阶Savitzky-Golay滤波器插值方法，将裂缝扩展向量转移到裂缝前沿网格点上(左侧箭头)，裂缝扩展向量与裂缝前沿网格点所在切线相垂直；

(7.2)、构建新的裂缝单元；

扩展前沿确定之后，需要构建新的裂缝单元。为了提升新增网格的质量，采用下述3类方法对新增人工裂缝网格形态进行调整：劈分、合并和扩展(图4)。图4中，新增单元特殊操作：(a)劈分，(b)合并，(c)扩展：(c-1)直接扩展；(c-2)扩展后重新劈分。

上述操作进行的前提条件如下：

式中Δd为扩展向量的长度，L_new为新单元相邻前沿点的距离，当扩展后扩展向量大于网格边界的2倍，即再补充进行劈分操作。

当人工裂缝与天然裂缝相交时，由于相邻扩展向量与天然裂缝的交点间可能存在多个天然裂缝网格点，因此在重构人工裂缝前沿新网格时需重构多个裂缝单元，图5展示了两种可能出现的情形。图5中，(a)相邻人工裂缝扩展向量均不为0；(b)一侧人工裂缝前沿已与天然裂缝相交。当相邻两人工裂缝前沿扩展向量均为不0时，主要由较长扩展向量对应的前沿点与交线上天然裂缝的网格点连接，构成新的人工裂缝单元；若有一侧人工裂缝前沿已与天然裂缝相交，则所有新增单元均与未相交的人工裂缝前沿网格点连接；

当人工裂缝与天然裂缝边界相交时，为避免出现网格交叉的情形，需对人工裂缝扩展向量方向进行调整。如图6所示，若相邻网格点中有一个网格点的扩展向量相交与天然裂缝平面外时，需要对相交与平面外的扩展向量或相交于平面内的扩展向量进行移动；使新增网格统一处于平面内或平面外；移动的目标点C为交点A2，B2连线与天然裂缝边界的交点。图6交点的移动：(a)A2-C>B2-C；(b)B2-C>A2-C。

同时，当人工裂缝越过天然裂缝后，有可能在天然裂缝的“背面”再次与天然裂缝相交(点A)，或与其他的人工裂缝壁面重合(点B与点C)(图7)。对于点A型裂缝前沿，若与天然裂缝的相交与原始交线较为接近，则在原始交线上寻找与交点最近的点，并将扩展向量移动至该处。若交点与原始交线距离较远，则需构建一条新的交线。

(8)、进入下一个时间步的计算。

计算结果展示

利用本实施例所列方法，对人工裂缝与天然裂缝相交过程进行模拟，假设该参数条件下人工裂缝无法直接穿过天然裂缝，模型参数如表1所示。

表1演示算例模型参数

图11与图12分别展示了不同注液时间对应的裂缝形态及开度分布。从模拟结果可以看出，本专利可高效、准确的对相交裂缝网格进行动态的、自适应的优化调整，可实现对人工裂缝与天然裂缝相交、人工裂缝越过天然裂缝、人工裂缝与天然裂缝再相交及人工裂缝合并等复杂过程的模拟。同时，该方法可实现对人工裂缝与天然裂缝开度及扩展方向的准确判断，提高边界元算法对水力压裂裂缝形态的模拟能力，尤其是对全三维空间裂缝形态的捕捉。该专利所述方法为利用边界元方法进行三维裂缝相交模拟提供重要的技术支撑，可推动全三维裂缝扩展模拟的进一步发展。

以上示意性的对本发明及其实施方式进行了描述，该描述没有限制性，附图中所示的也只是本发明的实施方式之一，实际的结构并不局限于此。所以，如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结构方式及实施例，均应属于本发明的保护范围。

Claims (6)

1.一种基于边界元方法的全三维裂缝相交过程模拟方法，其特征在于：包括以下步骤：

(1)时间步与注入压力初始化；将上个时间步n的时间步长Δtⁿ和注入点压力

作为当前时间步n+1的初始猜测，下标inj表示注入；

(2)裂缝开度计算；使所有裂缝单元压力均与注入点压力相等，采用边界元算法计算当前时间步长及注入点压力下裂缝开度

下标i为裂缝单元编号，k为更新注入点压力的的迭代步编号，j为更新时间步长的迭代步编号；

(3)质量守恒检查；检查当前压力迭代步对应注入压力下，当前时间步长内注液总量与裂缝体积总变化量是否相等，若相等则进入下一步，否则采用二分法对注入点压力进行更新，并返回步骤(2)；

(4)判断是否满足裂缝扩展条件；计算当前时间步长下，人工裂缝尖端的最大等效应力集中因子

与天然裂缝尖端单元的最大F值

下标HF与NF分别代表人工裂缝与天然裂缝；检查注入当前时间步长流体后，人工裂缝或天然裂缝是否能发生拟稳态扩展，若不能，则采用二分法对时间步长进行更新，并返回步骤(2)，直到人工裂缝或天然裂缝满足拟稳态扩展条件，则进入下一步；

人工裂缝尖端的最大等效应力集中因子

与天然裂缝尖端单元的最大F值

分别计算如下：

采用最大主应力原理判断人工裂缝是否发生扩展，最大主应力计算公式为：

式中σ_θ、σ_z、τ_θz分别为径向坐标系下的轴向应力、垂向应力与切向应力；σ_θ、σ_z、τ_θz可基于I型、II型与III型断裂应力集中因子K_I，K_II和K_III进行计算：

式中r与θ分别为裂缝尖端径向坐标系下的轴向距离与上倾角；应力集中因子可根据位移不连续量D_i计算：

式中E和v分别为储层岩石的杨氏模量与泊松比；使得最大主应力取最大值的倾角θ₀为裂缝扩展方向，需要满足下列条件：

倾角θ₀下裂缝尖端等效应力集中因子K_eq用于判断裂缝是否扩展及相应扩展距离：

对于天然裂缝尖端单元，采用F判据判断天然裂缝单元是否激活：

对于天然裂缝面内的裂缝单元激活有θ＝0⁰；式中K_IC与K_IIC分别为天然裂缝的I型与II型断裂的临界应力集中因子；当F值大于1时，天然裂缝单元被激活；

(5)裂缝扩展计算；判断裂缝是否扩展及计算裂缝扩展长度；

(6)自适应调整天然裂缝网格；采用逐次节点捕获技术对于人工裂缝相交后的天然裂缝平面网格进行自适应调整与优化；

如果人工裂缝相邻前沿网格点均恰落在天然裂缝平面的边界位置；初始时刻，天然裂缝网格并未沿着交线分布；天然裂缝网格点的调整分为两类：a.移动网格点与交线端点重合；b.调整网格点使之落在交线内部，对于上述两类情形，天然裂缝网格调整均需经历以下3个步骤：

a.1、找到网格平面内的捕获点；

a.2、寻找与捕获点相邻的图距离为N_dis的节点，定义为松弛节点；

a.3、等间距移动捕获点至目标点位置；在松弛过程中，采用优化算法对被松弛网格总体质量进行优化；如果总体网格质量达不到要求，则增加松弛节点的搜索图距离，重复a.2与a.3，直至达到最大图距离N_relax；

a.3步骤中的网格优化采用Matlab自带的带约束非线性优化函数fmincon实现，对于位于平面内部的节点，约束为移动后的节点不得脱离裂缝平面；对于位于平面边界的节点，约束条件为移动后的节点需仍位于平面边界处；优化函数为：

下标R代表被松弛的节点或网格组合，v代表节点，κ代表网格；B是满足所有约束条件的点的集合，网格质量Q定义如下：

式中

δ等于0.3，A是网格面积，L_j,j＝1,2,3是网格κ的边长度；Q越大，网格越接近等边三角形；当所有被松弛网格中有网格质量Q小于0.4，则认为网格质量未达到标准，需扩大搜索捕获点的图距离，重新更新节点位置；

对于a类网格调整，a.1中的捕获点为与端点最近的节点，a.3中的目标点为交线端点；

对于b类网格，捕获点搜索方法如下：网格平面空间被交线上最后一个节点与交线分为三个区域，三个区域内分别对应S＝1、0、-1；与交线上最后一个节点相连的单元若有边界跨过交线，则从该边界上的两个节点中寻找捕获点；将与交线较为靠近的节点选为捕获点；目标点为捕获点在交线上的投影；若目标位置过于靠近交线端点或落在交线外侧，则按照a类网格调整方法，移动捕获点与端点重合；

(7)自适应调整人工裂缝网格；

(7.1)对人工裂缝扩展前沿进行平滑处理，通过3阶Savitzky-Golay滤波器插值方法，将裂缝扩展向量转移到裂缝前沿网格点上，裂缝扩展向量与裂缝前沿网格点所在切线相垂直；

(7.2)构建新的裂缝单元；

当人工裂缝与天然裂缝相交时，需重构多个裂缝单元，当相邻两人工裂缝前沿扩展向量均为不0时，由较长扩展向量对应的前沿点与交线上天然裂缝的网格点连接，构成新的人工裂缝单元；若有一侧人工裂缝前沿已与天然裂缝相交，则所有新增单元均与未相交的人工裂缝前沿网格点连接；

当人工裂缝与天然裂缝边界相交时，需对人工裂缝扩展向量方向进行调整，若相邻网格点中有一个网格点的扩展向量相交与天然裂缝平面外时，需要对相交与平面外的扩展向量或相交于平面内的扩展向量进行移动；使新增网格统一处于平面内或平面外；

(8)进入下一个时间步的计算。

2.根据权利要求1所述的一种基于边界元方法的全三维裂缝相交过程模拟方法，其特征在于：步骤(2)中，裂缝开度计算的方法为：

设裂缝平面由一系列三角形单元网格构成，每个长方形单元有3个方向的位移不连续量：与裂缝面平行且相互正交的切向位移不连续量D_x与D_y以及垂直于裂缝平面的法向位移不连续量D_z；

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

u_x＝C{[2(1-v)I_,z-zI_,xx]D_x-zI_,xyD_y-[(1-2v)I_,x+zI_,xz]D_z}

u_y＝C{-zI_,xyD_x+[2(1-v)I_,z-zI_,yy]D_y-[(1-2v)I_,y+zI_,yz]D_z}

u_x＝C{[(1-2v)I_,x-zI_,xz]D_x+[(1-2v)I_,y-zI_.yz]D_y+[2(1-v)I_,z-zI_,zz]D_z}

σ_xx＝2C{[2I_,xz-zI_,xxx]D_x+[2vI_,yz-zI_,xxy]D_y+[I_,zz+(1-2v)I_,yy-zI_,xxz]D_z}

σ_yy＝2C{[2vI_,xz-zI_,xyy]D_x+[2I_,yz-zI_,yyy]D_y+[I_,zz+(1-2v)I_,xx-zI_,yyz]D_z}

σ_zz＝2C{-zI_,xzzD_x-zI_,yzzD_y+[I_,zz-zI_,zzz]D_z}

τ_xy＝2C{[(1-v)I_,yz-zI_,xxy]D_x+[(1-v)I_,xz-zI_,xyy]D_y-[(1-2v)I_,xy+zI_,xyz]D_z}

τ_yz＝2C{-[vI_,xy+zI_,xyz]D_x+[I_,zz+vI_,xx-zI_,yyz]D_y-zI_,yzzD_z}

τ_zx＝2C{[I_,zz+vI_,yy-zI_,xxz]D_x-[vI_,xy+zI_,xyz]D_y-zI_,xzzD_z}；

其中，C为与岩石力学性质相关的常数：

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

S为当前三角形单元面积；单元i满足应力边界条件，所受法向应力为

对应裂缝内压强大小，沿局部坐标系x、y方向的切向应力分别为τ_x ⁱ和τ_y ⁱ；此时所有裂缝单元作用在单元i的应力总和加上就地应力，应等于作用在单元表面的应力边界条件：

上式中A为影响系数；裂缝单元i的开度w_i对应单元z方向位移不连续量

的负数，即

3.根据权利要求2所述的一种基于边界元方法的全三维裂缝相交过程模拟方法，其特征在于：判断人工裂缝或天然裂缝是否能发生拟稳态扩展的方法为：

最大等效应力集中因子

或天然裂缝尖端单元的最大F值

在该时间步内取值在区间[(1-α,1+α)]内，其中α为接近0的小数。

4.根据权利要求3所述的一种基于边界元方法的全三维裂缝相交过程模拟方法，其特征在于：步骤(5)中人工裂缝或天然裂缝扩展距离及是否扩展的判断依据如下：

式中下标HF与NF分别代表人工裂缝与天然裂缝；d_max是单步最长扩展距离。

5.根据权利要求1所述的一种基于边界元方法的全三维裂缝相交过程模拟方法，其特征在于：步骤(7.2)中，采用下述3类方法对新增人工裂缝网格形态进行调整：劈分、合并和扩展，上述操作进行的前提条件如下：

式中Δd为扩展向量的长度，L_new为新单元相邻前沿点的距离，当扩展后扩展向量大于网格边界的2倍，即再补充进行劈分操作。

6.根据权利要求1所述的一种基于边界元方法的全三维裂缝相交过程模拟方法，其特征在于：步骤(7.2)中，当人工裂缝越过天然裂缝后，天然裂缝的“背面”再次与天然裂缝相交时，若交点与原始交线近，则在原始交线上寻找与交点最近的点，并将扩展向量移动至该处；若交点与原始交线距离远，则需构建一条新的交线。

Images