CN112100890B - 一种压裂水平缝起裂扩展全三维数学模型的计算方法 - Google Patents

Abstract

本发明涉及采油工程技术领域，特别是属于一种压裂水平缝起裂扩展全三维数学模型的计算方法。本发明包括：1)建立弹性岩石力学方程；2)建立物质流量连续性方程；3)采用三节点三角形等参单元，对裂缝扩展区域进行网格剖分；4)采用Galerkin即伽辽金方法有限元法并结合网格剖分的节点数，对弹性岩石力学方程和物质流量连续性方程进行离散；5)采用新的皮卡迭代方法，对离散的方程进行求解；6)利用求解结果计算裂缝延伸判断准则，判断裂缝扩展前缘是否稳定；7)如果步骤6)中裂缝前缘不稳定则重新进行步骤3)到6)；如果步骤6)中裂缝前缘稳定则停止计算。本发明具有求解速度快，迭代次数少，计算精度高等优点。

Description

一种压裂水平缝起裂扩展全三维数学模型的计算方法

技术领域：

本发明涉及采油工程技术领域，特别是属于一种压裂水平缝起裂扩展全三维数学模型的计算方法。

背景技术：

上世纪80年代以前，国内外在进行水力压裂设计计算上均采用二维模型，代表性的模型包括PKN模型和KGD模型，但这些模型均用于模拟垂直裂缝情况，大量实验和矿场研究表明，裂缝的高度不是定值，因此限制了模型的使用，同时没有水平裂缝二维的经典计算模型。对于大型压裂过程，需要采用三维模型进行合理的模拟，三维模型包括拟三维和全三维模型，拟三维模型假设裂缝三维延伸和流体一维流动，而全三维模型则考虑三维裂缝扩展以及裂缝内流体的二维流动，并且将滤失因素考虑进去，压裂液垂直于裂缝壁面渗滤进入地层，适用于各种储层条件，更加真实地模拟水力压裂的物理过程。

近年来，对裂缝的起裂扩展机理、压裂液的滤失以及温度场变化的研究大部分采用拟三维模型，假设条件过多，不适用于现场工程应用。

因此，由于传统全三维压裂裂缝起裂扩展三维数学模型求解计算量大，影响采油工程矿场实际应用，急需对其数学模型进行设计，对其计算解法进行改进，目的是提高求解速度和精度。

发明内容：

本发明针对传统全三维压裂裂缝起裂扩展三维数学求解计算量大的不足，提出了一种压裂水平缝起裂扩展全三维数学模型的计算方法，达到求解速度快，迭代次数少，计算精度高的目的。

为实现上述目的，本发明采用以下技术方案：一种压裂水平缝起裂扩展全三维数学模型的计算方法，其特征在于，包括以下步骤：

步骤一：建立弹性岩石力学方程；

步骤二：建立物质流量连续性方程；

步骤三：采用三节点三角形等参单元，对裂缝扩展区域进行网格剖分；

步骤四：采用Galerkin即伽辽金方法有限元法并结合网格剖分的节点数，对弹性岩石力学方程和物质流量连续性方程进行离散；

步骤五：采用新的皮卡迭代方法，对离散的方程进行求解；

步骤六：利用求解结果计算裂缝延伸判断准则，判断裂缝扩展前缘是否稳定；

步骤七：如果步骤六中裂缝前缘不稳定则重新进行步骤三到六；如果步骤六中裂缝前缘稳定则停止计算。

所述步骤一的建立弹性岩石力学方程为：

式中：∫为积分符号；A为水平裂缝扩展积分平面面积，m²；x,y为网格交点即节点的横坐标值和纵坐标值，m；p(x,y)为裂缝内部压裂液的液体压力，MPa；σ(x,y)为垂直于裂缝壁面的最小主应力，Mpa；V(x,y)为势函数；dx和dy为对x和y求微分；E为岩石的等效弹性模量，MPa；R为被积函数积分点(x',y')与压力作用点(x,y)之间的距离，

m；x',y'为被积函数积分点的横坐标值和纵坐标值，m；w(x,y)为水平裂缝高度，m；

为求导运算；dx'和dy'为对x'和y'求微分。

所述步骤二的建立物质流量连续性方程为：

式中：μ为压裂液的粘度，mPa·s；q_I为压裂液的注入量，m/s；t为压裂施工总时间，s；c_l为压裂液综合滤失系数，

τ(x,y)为裂缝壁面上某点(x，y)开始滤失的时刻，s。

所述步骤三的采用三节点三角形等参单元，对裂缝扩展区域进行网格剖分为：

采用三节点三角形单元进行网格划分；考虑井筒存在的情况，对井筒附近进行网格加密处理；将裂缝区域主要分成两大块：流体流动区域和裂缝前缘端部区域，其中流体流动区域网格比裂缝前缘端部区域的网格稀疏，以此对裂缝前缘处进行细致的模拟，提高模拟的精度并降低计算的数据规模，剖分后得到网格节点的总数为N。

所述步骤四的采用Galerkin即伽辽金方法有限元法并结合网格剖分的节点数，对弹性岩石力学方程和物质流量连续性方程进行离散为：

对方程(1)和(2)进行Galerkin有限元法离散，采用线性二元多项式(形函数φ(x,y))模式进行插值计算，并且裂缝高度和裂缝内的压力的形函数形式一致，同时令势函数等于形函数即V(x,y)＝φ(x,y)；

方程(1)离散形成以下矩阵方程：

[B]{w}＝{f} (3)

式中：B为N×N阶矩阵，矩阵元素用B_ij表示；w为N×1阶矩阵，矩阵元素用w_i表示；f为N×1阶矩阵，矩阵元素用f_i表示；各矩阵元素表达式为：

f_i＝-∫_A[p(x,y)-σ(x,y)]φ_idxdy (5)

式中：i＝1,2,...N；j＝1,2,...N；N为网格剖分的节点总数；

设其中某个三角形单元的三个节点的编号为I，J，K，(I＝1,2,...N，J＝1,2,...N，K＝1,2,...N)，坐标分别为(x_I,y_I)、(x_J,y_J)、(x_K,y_K)，则

φ_i＝φ_I＝1-φ_J-φ_K (6)

|J|＝(x_J-x_I)(y_J-y_I)-(x_K-x_I)(y_J-y_I) (9)

则式(4)中的

方程(2)离散形成以下矩阵方程：

[K]{p}＝-{f_l}-{f_w}+{f_q} (10)

式中：K为N×N阶矩阵，矩阵元素用K_ij表示；p为N×1阶矩阵，矩阵元素用p_i表示；f_l为N×1阶矩阵，矩阵元素用f_li表示；f_w为N×1阶矩阵，矩阵元素用f_wi表示；f_q为N×1阶矩阵，矩阵元素用f_qi表示；各矩阵元素表达式为：

式中：

为压裂液的注入积分平面面积，m²；ds为对面积求微分；i＝1,2,...N；j＝1,2,...N；N为网格剖分的节点总数；

其他形式的求导运算结合公式(6)-公式(9)得出。

所述步骤五的采用新的皮卡迭代方法，对离散的方程进行求解为:

当施工进行到时间t时刻时，节点i处的裂缝高度w_i和压力p_i的具体形式分别用

和

表示。联立求解方程(3)和方程(10)，采用如下迭代方法：

将迭代参数顺序调换，首先假设裂缝高度的分布值

再根据方程(3)计算压力分布值

再根据得到的压力分布值

带入方程(10)得到新的裂缝各点高度分布值

再根据方程(3)计算新的压力分布值

如此完成了一次迭代过程。迭代收敛所满足的条件为：

式中：

为t时刻在第n次迭代时第i节点的裂缝高度值，m；

为t时刻在第n次迭代时第i节点的压力值，MPa；

为t时刻在第n+1次迭代时第i节点的裂缝高度值，m；

为t时刻在第n+1次迭代时第i节点的压力值，MPa；

为求和运算；N为网格剖分的节点总数；ε为指定的公差，1×10^-4；

初始裂缝高度

的假设：根据斯内登即Sneddon关于裂缝高度与作用在缝壁面上的正应力的关系的研究，假设在无限大均质单一的弹性介质中，压开一周域为圆的水平裂缝，如果作用在缝壁上的正应力，为函数σ(r)，那么裂缝的高度由下式确定：

式中：σ(r)为半径为r处的正应力值，MPa；w(r)为半径为r处的裂缝高度，m；r为半径，

m；v为岩石泊松比；R'为裂缝最大半径，m；E为岩石弹性模量，Mpa；μ为压裂液的粘度，mPa·s；ρ为比值

小于1；λμ为表示裂缝半径分数的积分变量；π为圆周率，3.1415926；

当r＝0时裂缝的高度最大，令方程(16)中ρ＝0并积分简化得到裂缝的最大高度与作用在缝壁上的正应力的关系：

式中：w_max为坐标为(0,0)点处的裂缝高度值，为裂缝的最大高度值，m；arccos()反三角函数中的反余弦；

只要知道了缝壁上的正应力分布函数σ(r)，就可以求出最大高度值w_max；

假设初始裂缝纵断面为抛物线，裂缝沿径向高度的变化由下式确定：

根据步骤三网格剖分结果，得到了某节点的编号i和坐标值(x,y),此时

则

通过公式(19)就得到假设裂缝高度的初始分布值

至此，通过提出的新的皮卡迭代方法，耦合方程(3)和方程(10)，迭代收敛后，就可以得到时间t后裂缝内各点的流体压力值

和裂缝高度值

接下来根据裂缝前缘处的应力强度因子判断裂缝是否继续延伸。

所述步骤六的利用求解结果计算裂缝延伸判断准则，判断裂缝扩展前缘是否稳定为：

采用线弹性断裂力学的断裂准则作为裂缝扩展的判断依据，即扩展点处的应力强度因子K_I保持为近似等于临界应力强度因子K_IC；裂缝边界上任意一点的应力强度因子为：

式中：G为岩石的剪切模量，MPa；

为步骤五中求得的所有裂缝高度的最小值，m；X，Y为裂缝高度最小值的节点的横坐标值和纵坐标值，m。

所述步骤七如果步骤六中裂缝前缘不稳定则重新进行步骤三到六；如果步骤六中裂缝前缘稳定则停止计算为：

利用公式(20)计算得到扩展点处的应力强度因子K_I值，如果K_I小于等于临界应力强度因子K_IC值，则计算停止；如果K_I大于临界应力强度因子K_IC值，则计算裂缝前缘各点向前延伸的距离并重新进行步骤三到六计算，直到K_I小于等于临界应力强度因子K_IC值。

本发明与已有技术相比具有如下有益效果：

1)采用Galerkin有限元法对弹性岩石力学方程和物质流量连续性方程进行离散，采用微分方程对应的弱形式，通过选取有限多项试函数，又称基函数或形函数，将它们叠加，要求结果在求解域内及边界上的加权积分满足原方程，便可以得到一组易于求解的线性代数方程，且自然边界条件能够自动满足,以上处理极大的提高了数值求解的精度和速度；2)采用三节点三角形等参单元对裂缝起裂扩展区域进行网格剖分，由于三角形单元对于复杂边界有较强的适应能力，很容易将一个二维区域离散成有限个三角形单元，在边界上以若干段直线近似原来的曲线边界，随着单元的增多，拟合的越来越精确，因此提高了数值求解的精度；3)提出了新的皮卡迭代方法，避免了传统计算方法设定初始压力值均等，导致后期迭代次数过多，增加的计算繁琐性，因此提高了计算的速度。

附图说明：

图1为计算物理模型图；图2为三角形单元剖分结果图。

具体实施方式：

下面结合附图及实施例子将对本发明作进一步说明：

一种压裂水平缝起裂扩展全三维数学模型的计算方法，包括以下步骤：

步骤一：建立弹性岩石力学方程；

步骤二：建立物质流量连续性方程；

步骤三：采用三节点三角形等参单元，对裂缝扩展区域进行网格剖分；

步骤四：采用Galerkin即伽辽金方法有限元法并结合网格剖分的节点数，对弹性岩石力学方程和物质流量连续性方程进行离散；

步骤五：采用新的皮卡迭代方法，对离散的方程进行求解；

步骤六：利用求解结果计算裂缝延伸判断准则，判断裂缝扩展前缘是否稳定；

步骤七：如果步骤六中裂缝前缘不稳定则重新进行步骤三到六；如果步骤六中裂缝前缘稳定则停止计算。

以压裂施工F1井为例，具体如下：

1、所述建立弹性岩石力学方程；具体如下：

式中：∫为积分符号；A为水平裂缝扩展积分平面面积，m²；x,y为如图2所示的网格交点即节点的横坐标值和纵坐标值，m；p(x,y)为裂缝内部压裂液的液体压力，MPa；σ(x,y)为垂直于裂缝壁面的最小主应力，Mpa；V(x,y)为势函数；dx和dy为对x和y求微分；E为岩石的等效弹性模量，MPa；R为被积函数积分点(x',y')与压力作用点(x,y)之间的距离，

m；x',y'为被积函数积分点的横坐标值和纵坐标值，m；w(x,y)为水平裂缝高度，m；

为求导运算；dx'和dy'为对x'和y'求微分。

2、所述建立物质流量连续性方程；具体如下：

式中：μ为压裂液的粘度，mPa·s；q_I为压裂液的注入量，m/s；t为压裂施工总时间，s；c_l为压裂液综合滤失系数，

τ(x,y)为裂缝壁面上某点(x，y)开始滤失的时刻，s。

3、所述采用三节点三角形等参单元，对裂缝扩展区域进行网格剖分；具体如下：

假设裂缝的扩展半径为30m，采用三节点三角形单元进行网格划分，考虑井筒存在的情况，将裂缝区域主要分成两大块，流体流动区域和裂缝前缘端部区域。其中流体流动区域网格较近裂缝端部区域网格稀疏，以此对裂缝前缘处进行细致的模拟。在本实例中，共计对裂缝扩展区域剖分成256个三角形单元，单元的总结点数为144。计算物理模型如图1所示，其中井筒半径实际为10cm，显示放大了30倍，利用三节点三角形单元对物理模型进行网格剖分的结果如图2所示。

4、所述采用Galerkin即伽辽金方法有限元法并结合网格剖分的节点数，对弹性岩石力学方程和物质流量连续性方程进行离散；具体如下：

对方程(1)和(2)进行Galerkin有限元法离散，采用线性二元多项式即形函数φ(x,y)模式进行插值计算，并且裂缝高度和裂缝内的压力的形函数形式一致，同时令势函数等于形函数即V(x,y)＝φ(x,y)。

方程(1)离散形成以下矩阵方程：

[B]{w}＝{f} (3)

式中：B为N×N阶矩阵，矩阵元素用B_ij表示；w为N×1阶矩阵，矩阵元素用w_i表示；f为N×1阶矩阵，矩阵元素用f_i表示；本实施例子中N为144；各矩阵元素表达式为：

f_i＝-∫_A[p(x,y)-σ(x,y)]φ_idxdy (5)

式中：i＝1,2,...N；j＝1,2,...N；N为网格剖分的节点总数，本实施例子中N为144。

设其中某个三角形单元的三个节点的编号为I，J，K，(I＝1,2,...N，J＝1,2,...N，K＝1,2,...N)，坐标分别为(x_I,y_I)、(x_J,y_J)、(x_K,y_K)，则：

φ_i＝φ_I＝1-φ_J-φ_K (6)

|J|＝(x_J-x_I)(y_J-y_I)-(x_K-x_I)(y_J-y_I) (9)

则式(4)中的

以图2中的某个单元为例，此三角形单元的三个节点的编号为I＝1、J＝2、K＝3，坐标分别为(17,0)、(20,4)、(27,0)，即x_I＝17,y_I＝0、x_J＝20,y_J＝4、x_K＝27,y_K＝0，则：

|J|＝(x_J-x_I)(y_J-y_I)-(x_K-x_I)(y_J-y_I)＝-28；

方程(2)离散形成以下矩阵方程：

[K]{p}＝-{f_l}-{f_w}+{f_q} (10)

式中：K为N×N阶矩阵，矩阵元素用K_ij表示；p为N×1阶矩阵，矩阵元素用p_i表示；f_l为N×1阶矩阵，矩阵元素用f_li表示；f_w为N×1阶矩阵，矩阵元素用f_wi表示；f_q为N×1阶矩阵，矩阵元素用f_qi表示；本实施例子中N为144；各矩阵元素表达式为：

式中：

为压裂液的注入积分平面面积，m²；ds为对面积求微分；i＝1,2,...N；j＝1,2,...N；N为网格剖分的节点总数，144；

其他形式的求导运算结合公式(6)-公式(9)得出。

以上离散方程中的数学符号的物理意义、单位及数值如表1所示。

5、所述采用新的皮卡迭代方法，对离散的方程进行求解；具体如下：

离散的方程包括方程(3)和方程(10)，对两个方程组成的方程组进行联立，并采用迭代法求解。

首先假设裂缝高度的分布值

其具体计算过程如下：

1)根据缝壁上的正应力分布函数σ(r)，求出最大高度值w_max；

式中：σ(r)为半径为r处的正应力值，MPa；r为半径，

m；v为岩石泊松比；R'为裂缝最大半径，m；E为岩石弹性模量，Mpa；w_max为坐标为(0,0)点处的裂缝高度值，为裂缝的最大高度，m；arccos()反三角函数中的反余弦。

将表1中的数据代入到公式(15)中，其中初始裂缝的最大半径R'为30m，计算得到最大高度值w_max为0.01424m。

2)假设初始裂缝纵断面为抛物线，裂缝沿径向高度的变化由下式确定：

式中：w(r)为半径为r处的裂缝高度，m。

例如以图2中的某个三角形单元，此三角形单元的三个节点的编号为I＝1、J＝2、K＝3，坐标分别为(17,0)、(20,4)、(27,0)，则：

根据节点的坐标，得到节点的半径值r，代入公式(16)中进行计算，计算得到节点1的初始裂缝假设的高度数值为0.009374m；节点2的初始裂缝假设的高度数值为0.008057m；节点3的初始裂缝假设的高度数值为0.004503m；其他三角形单元节点的初始裂缝的假设值同理得到。

通过以上计算得到了施工时间为6000s时刻的假设裂缝高度的初始分布值

将各节点在时间t＝6000s时刻的裂缝高度的分布值

带入到公式(3)中，对形成的方程组进行求解，计算得到各节点的压力分布值

再根据得到的压力分布值

带入公式(10)得到新的裂缝各点高度分布值

再根据公式(3)计算新的压力分布值

如此完成了一次迭代过程。以上矩阵方程的矩阵元素的积分和求解过程通过Excel表格实现。

迭代计算过程对计算结果进行收敛判断，依据如下：

式中：

为t时刻在第n次迭代时第i节点的裂缝高度值，m；

为t时刻在第n次迭代时第i节点的压力值，MPa；

为t时刻在第n+1次迭代时第i节点的裂缝高度值，m；

为t时刻在第n+1次迭代时第i节点的压力值，MPa；

为求和运算；N为网格剖分的节点总数，本例中为144；ε为指定的公差，1×10^-4。

通过一定次数的迭代，此时得到的

值，代入公式(17)中，计算得到：

此时迭代收敛。图2中的部分节点的半径和裂缝高度的迭代计算结果如表2所示。

6、所述利用求解结果计算裂缝延伸判断准则，判断裂缝扩展前缘是否稳定；具体如下：

采用线弹性断裂力学的断裂准则作为裂缝扩展的判断依据，即扩展点处的应力强度因子K_I保持为近似等于临界应力强度因子K_IC。裂缝边界上任意一点的应力强度因子为：

式中：G为岩石的剪切模量，MPa；

为步骤五中求得的所有裂缝高度的最小值，m；X，Y为裂缝高度最小值的节点的横坐标值和纵坐标值，m。

本例中，计算得到的节点的裂缝高度最小值为2.78mm，节点编号为7，其横坐标值和纵坐标值分别为27.5、5.4。将以上数据带入到公式(18)中，得到扩展点处的应力强度因子K_I：

7、所述如果步骤六中裂缝前缘不稳定则重新进行步骤三到六；如果步骤六中裂缝前缘稳定则停止计算。具体如下：

1)如果本例中临界应力强度因子K_IC值为

步骤6中计算得到扩展点处的应力强度因子K_I值为

K_I值大于K_IC，则计算裂缝前缘各点向前延伸的距离。

根据Mastorjannis，得到裂缝前缘各点向前延伸的距离为：

式中：Δd为裂缝前缘各点向前延伸的距离，m；σ为裂缝前缘的地应力值，MPa；h为裂缝的深度，m。

计算得到裂缝前缘各点向前延伸的距离Δd为：

此时假设新的裂缝半径值为30.3212m，并重新进行步骤3到6，直到K_I小于等于临界应力强度因子K_IC值，计算停止。

2)如果本例中的临界应力强度因子K_IC值为

K_I值小于K_IC，则裂缝不扩展，即Δd＝0，计算停止。

表1F1井地质和工程参数数据表

表2第一次迭代收敛后的图2中部分节点的半径和裂缝高度值

节点编号	1	3	4	5	6	7
							半径(m)	17	27	3	10	30	28
高度(mm)	4.92	3.26	5.76	5.53	0	2.78

Claims (8)

1.一种压裂水平缝起裂扩展全三维数学模型的计算方法，其特征在于，包括以下步骤：

步骤一：建立弹性岩石力学方程；

步骤二：建立物质流量连续性方程；

步骤三：采用三节点三角形等参单元，对裂缝扩展区域进行网格剖分；

步骤四：采用Galerkin即伽辽金方法有限元法并结合网格剖分的节点数，对弹性岩石力学方程和物质流量连续性方程进行离散；

步骤五：采用新的皮卡迭代方法，对离散的方程进行求解；

步骤六：利用求解结果计算裂缝延伸判断准则，判断裂缝扩展前缘是否稳定；

步骤七：如果步骤六中裂缝前缘不稳定则重新进行步骤三到六；如果步骤六中裂缝前缘稳定则停止计算。

2.根据权利要求1所述的一种压裂水平缝起裂扩展全三维数学模型的计算方法，其特征在于，所述步骤一的建立弹性岩石力学方程为：

式中：∫为积分符号；A为水平裂缝扩展积分平面面积，m²；x,y为网格交点即节点的横坐标值和纵坐标值，m；p(x,y)为裂缝内部压裂液的液体压力，MPa；σ(x,y)为垂直于裂缝壁面的最小主应力，Mpa；V(x,y)为势函数；dx和dy为对x和y求微分；E为岩石的等效弹性模量，MPa；R为被积函数积分点(x',y')与压力作用点(x,y)之间的距离，

m；x',y'为被积函数积分点的横坐标值和纵坐标值，m；w(x,y)为水平裂缝高度，m；

为求导运算；dx'和dy'为对x'和y'求微分。

3.根据权利要求1所述的一种压裂水平缝起裂扩展全三维数学模型的计算方法，其特征在于，所述步骤二的建立物质流量连续性方程为：

式中：μ为压裂液的粘度，mPa·s；q_I为压裂液的注入量，m/s；t为压裂施工总时间，s；c_l为压裂液综合滤失系数，

τ(x,y)为裂缝壁面上某点(x，y)开始滤失的时刻，s。

4.根据权利要求1所述的一种压裂水平缝起裂扩展全三维数学模型的计算方法，其特征在于，所述步骤三的采用三节点三角形等参单元，对裂缝扩展区域进行网格剖分为：

采用三节点三角形单元进行网格划分；考虑井筒存在的情况，对井筒附近进行网格加密处理；将裂缝区域主要分成两大块：流体流动区域和裂缝前缘端部区域，其中流体流动区域网格比裂缝前缘端部区域的网格稀疏，以此对裂缝前缘处进行细致的模拟，提高模拟的精度并降低计算的数据规模，剖分后得到网格节点的总数为N。

5.根据权利要求1所述的一种压裂水平缝起裂扩展全三维数学模型的计算方法，其特征在于，所述步骤四的采用Galerkin即伽辽金方法有限元法并结合网格剖分的节点数，对弹性岩石力学方程和物质流量连续性方程进行离散为：

对方程(1)和(2)进行Galerkin有限元法离散，采用线性二元多项式(形函数φ(x,y))模式进行插值计算，并且裂缝高度和裂缝内的压力的形函数形式一致，同时令势函数等于形函数即V(x,y)＝φ(x,y)；

方程(1)离散形成以下矩阵方程：

[B]{w}＝{f} (3)

式中：B为N×N阶矩阵，矩阵元素用B_ij表示；w为N×1阶矩阵，矩阵元素用w_i表示；f为N×1阶矩阵，矩阵元素用f_i表示；各矩阵元素表达式为：

f_i＝-∫_A[p(x,y)-σ(x,y)]φ_idxdy (5)

式中：i＝1,2,...N；j＝1,2,...N；N为网格剖分的节点总数；

设其中某个三角形单元的三个节点的编号为I，J，K，(I＝1,2,...N，J＝1,2,...N，K＝1,2,...N)，坐标分别为(x_I,y_I)、(x_J,y_J)、(x_K,y_K)，则

φ_i＝φ_I＝1-φ_J-φ_K (6)

|J|＝(x_J-x_I)(y_J-y_I)-(x_K-x_I)(y_J-y_I) (9)

则式(4)中的

方程(2)离散形成以下矩阵方程：

[K]{p}＝-{f_l}-{f_w}+{f_q} (10)

式中：K为N×N阶矩阵，矩阵元素用K_ij表示；p为N×1阶矩阵，矩阵元素用p_i表示；f_l为N×1阶矩阵，矩阵元素用f_li表示；f_w为N×1阶矩阵，矩阵元素用f_wi表示；f_q为N×1阶矩阵，矩阵元素用f_qi表示；各矩阵元素表达式为：

式中：

为压裂液的注入积分平面面积，m²；ds为对面积求微分；i＝1,2,...N；j＝1,2,...N；N为网格剖分的节点总数；

其他形式的求导运算结合公式(6)-公式(9)得出。

6.根据权利要求1所述的一种压裂水平缝起裂扩展全三维数学模型的计算方法，其特征在于，所述步骤五的采用新的皮卡迭代方法，对离散的方程进行求解为：

当施工进行到时间t时刻时，节点i处的裂缝高度w_i和压力p_i的具体形式分别用

和

表示， 联立求解方程(3)和方程(10)，采用如下迭代方法：

将迭代参数顺序调换，首先假设裂缝高度的分布值

再根据方程(3)计算压力分布值

再根据得到的压力分布值

带入方程(10)得到新的裂缝各点高度分布值

再根据方程(3)计算新的压力分布值

如此完成了一次迭代过程， 迭代收敛所满足的条件为：

式中：

为t时刻在第n次迭代时第i节点的裂缝高度值，m；

为t时刻在第n次迭代时第i节点的压力值，MPa；

为t时刻在第n+1次迭代时第i节点的裂缝高度值，m；

为t时刻在第n+1次迭代时第i节点的压力值，MPa；

为求和运算；N为网格剖分的节点总数；ε为指定的公差，1×10^-4；

初始裂缝高度

的假设：根据斯内登即Sneddon关于裂缝高度与作用在缝壁面上的正应力的关系的研究，假设在无限大均质单一的弹性介质中，压开一周域为圆的水平裂缝，如果作用在缝壁上的正应力，为函数σ(r)，那么裂缝的高度由下式确定：

式中：σ(r)为半径为r处的正应力值，MPa；w(r)为半径为r处的裂缝高度，m；r为半径，

m；v为岩石泊松比；R'为裂缝最大半径，m；E为岩石弹性模量，Mpa；μ为压裂液的粘度，mPa·s；ρ为比值

小于1；λμ为表示裂缝半径分数的积分变量；π为圆周率，3.1415926；

当r＝0时裂缝的高度最大，令方程(16)中ρ＝0并积分简化得到裂缝的最大高度与作用在缝壁上的正应力的关系：

式中：w_max为坐标为(0,0)点处的裂缝高度值，为裂缝的最大高度值，m；arccos()反三角函数中的反余弦；

只要知道了缝壁上的正应力分布函数σ(r)，就可以求出最大高度值w_max；

假设初始裂缝纵断面为抛物线，裂缝沿径向高度的变化由下式确定：

根据步骤三网格剖分结果，得到了某节点的编号i和坐标值(x,y),此时

则

通过公式(19)就得到假设裂缝高度的初始分布值

至此，通过提出的新的皮卡迭代方法，耦合方程(3)和方程(10)，迭代收敛后，就可以得到时间t后裂缝内各点的流体压力值

和裂缝高度值

接下来根据裂缝前缘处的应力强度因子判断裂缝是否继续延伸。

7.根据权利要求1所述的一种压裂水平缝起裂扩展全三维数学模型的计算方法，其特征在于，所述步骤六的利用求解结果计算裂缝延伸判断准则，判断裂缝扩展前缘是否稳定为：

采用线弹性断裂力学的断裂准则作为裂缝扩展的判断依据，即扩展点处的应力强度因子K_I保持为近似等于临界应力强度因子K_IC；裂缝边界上任意一点的应力强度因子为：

式中：G为岩石的剪切模量，MPa；

为步骤五中求得的所有裂缝高度的最小值，m；X，Y为裂缝高度最小值的节点的横坐标值和纵坐标值，m。

8.根据权利要求1所述的一种压裂水平缝起裂扩展全三维数学模型的计算方法，其特征在于，所述步骤七如果步骤六中裂缝前缘不稳定则重新进行步骤三到六；如果步骤六中裂缝前缘稳定则停止计算为：

利用公式(20)计算得到扩展点处的应力强度因子K_I值，如果K_I小于等于临界应力强度因子K_IC值，则计算停止；如果K_I大于临界应力强度因子K_IC值，则计算裂缝前缘各点向前延伸的距离并重新进行步骤三到六计算，直到K_I小于等于临界应力强度因子K_IC值。

Images