适用于任意形状倾斜裂缝面的三维嵌入式离散裂缝模型前处理算法

2020-04-08程林松曹仁义安小平雷征东

科学技术与工程 2020年3期

饶翔，程林松，曹仁义*，安小平，雷征东

(1.中国石油大学(北京)石油工程教育部重点实验室，北京 102249；2.中国石油大学(北京)石油工程学院，北京 102249；3.中国石油长庆油田分公司勘探开发研究院，西安 710018；4.中国石油勘探开发研究院，北京 100083)

在开展复杂裂缝网络传质的数值模拟时，为了能够准确刻画油藏中裂缝的形态，离散裂缝模型(DFM)需要生成高质量的非结构网格去匹配裂缝在油藏空间中的几何形态。对于在三维油藏中的大量倾斜裂缝及其形成的复杂缝网，与之匹配的高质量非结构网格的生成是十分困难的。如果需要考虑裂缝的动态行为，则需要在每一时间步重新生成非结构网格去匹配该时间步的缝网形态，会显著降低计算效率。而嵌入式离散裂缝模型(EDFM)仅需要利用结构化网格剖分基质系统，形成背景网格，将显式表征的裂缝嵌入到背景网格之中形成裂缝网格，再获取网格之间的连接及相应的传导系数，可以解决高质量非结构网格生成的问题，在考虑裂缝扩展或增删裂缝的情况下无需更新背景网格，模拟效率高，操作灵活。Lee等[1]首次提出将裂缝嵌入到基质网格中，采用与水平井相似的方式将裂缝视为基质网格中的源汇项来处理基质与裂缝之间的传质。随后，Li等[2]扩展了前述的初步想法首次提出了EDFM。在EDFM中，由于结构化网格的使用，一般采用既能满足局部物质守恒，又简单易行的有限体积、有限差分或者模拟有限差分等方法去获取渗流控制方程的离散格式。Hajibeygi等[3]利用循环多尺度有限体积方法建立了相应的EDFM，可以有效提高模拟计算效率。Moinfar等[4-5]首次将EDFM从二维扩展到三维，并简单分析了天然裂缝和水力压裂缝动态行为对渗流场及水平井产能的影响。EDFM已被用于开展多种渗流问题中复杂缝网传质的数值模拟研究[6-12]。Zhang等[13]、Yan等[14]利用模拟有限差分方法离散渗流控制方程，达到同时求解速度场和压力场的目的；Cao等[15]基于边界积分方程得到基质网格与裂缝网格之间传质量的离散格式；Rao等[16]基于格林元方法中对方程中含时项处理的思想提出了新的基质网格与裂缝网格传质量的近似格式；Rao等[17]利用EDFM研究水侵层对注水吞吐过程的影响。目前，EDFM已被广泛应用于各种复杂渗流问题的模拟研究中，然而Tene等[18]研究发现EDFM无法有效处理裂缝网格渗透率低于所在基质网格渗透率的情况，如油藏中的断层、隔夹层等。Jiang等[19]的研究也表明EDFM在多相流横穿裂缝时会出现明显误差。造成EDFM局限性的本质原因在于仅将裂缝处理为其所在基质网格中的源汇项，一般仅能适用于多段压裂水平井生产或注入过程的数值模拟。为了解决EDFM的上述局限性，Tene等[18]提出了基于投影的嵌入式离散裂缝模型(pEDFM)；Jiang等[19]对pEDFM做出了一些改进，改进之处主要是提出了裂缝单元投影面的选取原则以及修正了传导系数的计算公式。但pEDFM与EDFM的相似之处在于都是将裂缝嵌入到被结构化网格剖分的基质系统获取裂缝网格的分布，建立网格之间的连接，不同的是pEDF较EDFM添加了一类新的连接。

EDFM和pEDFM的前处理主体部分基本相同，从二维模型升级到三维模型的本质即是前处理算法的维度升级，由于三维油藏中裂缝面可以任意倾斜或者具有复杂的几何形状，这给前处理过程带来了挑战，并且在以往的三维EDFM中一般只有矩形平面缝，且均未给出相应的前处理算法细节。为此，提出一种可适用于由复杂形状倾斜缝构成的复杂裂缝网络的EDFM前处理算法，以期依据该前处理算法实现相应的三维EDFM。

1 裂缝面的参数化

在EDFM中，对裂缝采用降维处理，即在三维油藏中，用一个二维平面或曲面刻画一条裂缝(记作裂缝面)。一般情况下，三维油藏用三维欧式空间刻画(油藏空间尺度不大，无需考虑地面曲率的影响)，裂缝面是一个二维曲面。根据经典的曲面论，二维裂缝曲面可以用该曲面的第一、第二基本形式和该面上某点的向径唯一确定。然而，如果裂缝面用二维曲面表征，则难以获取裂缝面与基质网格之间的连接关系及相应裂缝网格的几何参数，并且在EDFM或者DFM中，裂缝一般是二维带边平面或者一维直线段，因此利用二维平面来刻画三维油藏中的裂缝，即一条裂缝对应着一个二维平面。而要确定一个三维欧式空间(油藏)中的二维平面，则需要该平面上某点A0的向径(记作基准向量β0)和该裂缝面上的两个线性无关(或正交)的向量β1和β2，这两个线性无关的向量实际上构成了该裂缝所在平面以A0为原点的局部坐标系的两个基向量，即该裂缝所在平面上的点可以表示为

r=β0+uβ1+vβ2

(1)

式(1)中：r为裂缝所在平面上点的矢径；u和v分别为β1和β2对应的参数，即偶对 (u,v)在局部坐标系{A0，β1,β2}中的坐标，与该裂缝所在平面上的点一一对应。

为了进一步确定裂缝在该平面上的形状，则需要给出u和v需要满足的方程函数关系f(u,v)≤0，其中，f(u,v)=0确定裂缝的边界形状，f(u,v)<0即是裂缝内部。

基于以上的裂缝面参数化表示，得出具体的前处理算法如下。

2 前处理具体算法

EDFM实质上是将DFM生成非结构网格匹配裂缝几何形态的困难转变为获取裂缝网格分布及网格之间连接关系的困难，这种转变使得EDFM避免了在复杂缝网情况下高质量非结构网格生成的困难，并且只需要一套预先设定的固定的背景网格，使得在求解大量问题时的计算效率更高。因此获取裂缝网格分布及网格之间连接关系是嵌入式离散裂缝模型前处理的关键，同时也说明从二维EDFM到三维EDFM维度升级的本质困难在于前处理算法的维度升级。

在获取裂缝网格分布及网格之间连接关系这个问题上，三维EDFM比二维EDFM要复杂得多，因为在二维EDFM中是二维的基质网格和一维直线段裂缝，故只需要求解一维直线段裂缝与基质网格各边的交点以及不同直线段裂缝之间的交点即可确定裂缝网格的分布以及网格之间的连接关系，但在三维EDFM中是三维的基质网格和二维的裂缝面，裂缝面与三维基质网格之间以及不同裂缝面之间的几何关系复杂。对此，三维EDFM前处理算法的基本思路分为六个步骤：①计算裂缝面(裂缝面指裂缝降维处理后的几何实体；裂缝所在平面是指裂缝面实体所在的三维欧式空间中的无限大平面)与基质网格线的交点；②计算裂缝面边界与基质网格面的交点；③确定预裂缝网格分布；④求解不同裂缝面之间的交线方程；⑤确定最终的裂缝网格分布；⑥确定网格之间连接关系。

按照上述的六个步骤，具体的算法细节如下。

2.1 计算裂缝面与基质网格线的交点

假设该裂缝面的参数化形式如下：

r=β0+uβ1+vβ2，f(u,v)≤0

(2)

式(2)中：β0=(β01,β02,β03)，β1=(β11,β12,β13)，β2=(β21,β22,β23)。

则计算裂缝面与基质网格线交点的具体算法步骤如下。

步骤1首先判断该裂缝面是否垂直于x-y平面Γxy、y-z平面Γyz或x-z平面Γxz，判断方法是：先计算该裂缝面法向量n=β1β2，若n·(0,0,1)=0，则该裂缝面垂直于Γxy，若n·(1,0,0)=0，则该裂缝面垂直于Γyz，若n·(0,1,0)=0，则该裂缝面垂直于Γxz。

步骤2计算网格线与裂缝所在平面的交点，以垂直于x-y平面的z方向网格线为例说明，对于某一条z方向网格线，其x和y是确定的，不妨假设为x0和y0，则将x0和y0代入式(2)，可以得到线性方程组：

(3)

若裂缝所在平面不垂直于x-y平面[图1(a)]，根据几何关系可知裂缝所在平面与该z方向网格线必有唯一交点，即上述线性方程组系数矩阵非奇异必有唯一解。解方程组得到(u0,v0)后，再判断(u0,v0)是否满足f(u0,v0)≤0，若不满足，则说明该z方向网格线与裂缝面无交点，若满足，则说明(u0,v0)对应的点是该z方向网格线与裂缝面的交点，则计算z0=β03+u0β13+v0β23，存储(x0,y0,z0)和(u0,v0)分别为该交点的三维欧式空间坐标和二维裂缝面上局部坐标。

图1 裂缝面与基质网格线之间的三种几何关系

若裂缝所在平面垂直于x-y平面，根据几何关系可知上述线性方程组系数矩阵奇异，或者无解,如图1(b)所示，对应于该网格线与裂缝所在平面没有交点；或者无穷多组解，如图1(c)所示，对应于该网格线包含在裂缝所在平面有无穷多个交点。此时需要判断是上述哪种情况，方法是任意取一个z0值，将y0和z0代入式(2)，可以得到一个线性方程组：

(4)

若式(4)中方程组的系数矩阵奇异，即说明原z方向网格线与裂缝所在平面没有交点，即无解；反之，则说明该z方向网格线包含在该裂缝平面内，是第二种情况。在无解情况下，由于没有交点，因此无需进行其他操作；在无穷多组解下，由于该z网格方向线包含在裂缝所在平面，因此通过二分法结合f(u0,v0)≤0的条件来计算z0的取值范围，从而将z0取值范围内的同时也属于基质网格顶点(格点)的(x0,y0,z0)以及(x0,y0,z0,max)、(x0,y0,z0,min)存储为该z方向网格线与裂缝面的交点，同时存储这些交点对应二维裂缝面上的局部坐标(u0,v0)[该局部坐标通过方程(3)求解得到]。

同理，对于x方向坐标线和y方向坐标线，采用类似的方法计算得到坐标线与裂缝面的交点。

2.2 计算裂缝面边界与基质网格面的交点

依然假设该裂缝面的参数化表示如式(2)所示，则计算裂缝面边界与基质网格面交点的具体算法步骤如下。

步骤1首先判断裂缝面Γf是否平行于x-y平面Γxy、y-z平面Γyz或x-z平面Γxz，其判断方法类似于2.1节的第一步；

步骤2接下来以某一x-y基质网格面Γ1(其z值是一个常数，记为z0)为例介绍计算裂缝面Γf与Γ1交点的算法。

图2 裂缝面与基质网格面之间的四种几何关系

若裂缝面Γf不与Γxy平行，若存在裂缝面某边界Γ1Γf在基质网格面Γ1内，如图2(a)所示，因为这两个面均是平面，故该边界是一条直线段，且该直线段一般与Γ1的网格线有交点，即这些交点在2.1节中已计算出来，此时无需任何操作；同时因为该边界是一条直线，故方程(5)是线性方程组，因此若该方程组的系数矩阵是否奇异则可以判断∂1Γf在基质网格面Γ1内，如图2(b)～图2(d)所示，在这种情况下由几何关系可知该裂缝面边界与Γ1至多有有限个交点，裂缝面Γf边界方程是f(u,v)=0，因此交点对应的裂缝面Γf上的局部坐标(u0,v0)满足下述方程组：

(5)

求解式(5)即可解出交点的二维裂缝面局部坐标(u0,v0)，再将(u0,v0)代入到裂缝面Γf参数化形式中，计算交点对应的三维欧式空间坐标(x0,y0,z0)。特殊的，当裂缝面Γf是矩形、平行四边形或者其他多边形时，裂缝面Γf边界是由直线段构成的，因此上述方程组实际上是线性方程组，易于求解。若裂缝面Γf边界是椭圆或者更高次代数曲线时，则需要求解二次或更高次的二元代数方程组，此时有两种方法，第一种是利用NR迭代计算，但需要注意若迭代不收敛，则说明裂缝面Γf边界与Γxy没有交点；第二种是先利用方程组中的第一个线性方程，用u0表示v0，或者用v0表示u0，再代入f(u0,v0)=0中得到一元二次或更高次方程求解实数根，若没有实数根，则说明裂缝面Γf边界与Γxy没有交点。此处是三维EDFM前处理中矩形缝与椭圆形等复杂形状缝在几何参数计算中的唯一不同之处。

若裂缝面平行于Γxy，则可知裂缝面边界与Γ1要么没有交点，此时无需任何操作；要么无穷多个交点(即该裂缝面边界全包含在Γ1内)，这种情况在EDFM中一般不会出现，因为此时裂缝面与基质网格面重合，即裂缝面在几何形态上成为了基质网格的边界。

2.3 确定预裂缝网格分布

在前述基质网格线与每个裂缝面之间以及基质网格面与每个裂缝面边界之间的所有交点计算得到后，判断同一裂缝面上的所有交点所属的基质网格编号，在同一基质网格中的所有交点构成的一个凸多边形，定义为该裂缝面上的一个预裂缝网格(此处是预裂缝网格，是因为最终的裂缝网格还需要考虑不同裂缝面之间的交线对预裂缝网格的进一步剖分)，一个裂缝面上的全体预裂缝网格构成该裂缝面上的预裂缝网格分布。如图3所示，点1、2、3、4、5、6、7是通过2.1节算法计算得到的该椭圆裂缝面与基质网格线的交点，点8、9、10、11、12、13、14、15、16、17是通过2.2节算法计算得到的椭圆形裂缝面与基质网格面的交点，然后属于同一基质网格的节点构成一个预裂缝网格，从而得到12个预裂缝网格，分别是：1-2-8、1-2-9、2-3-10-8、2-3-11-9、3-4-12-10、3-4-13-11、4-5-14-12、4-5-15-13、5-6-16-14、5-6-17-1、6-7-16、6-7-17。

图3 椭圆形裂缝面预裂缝网格示意图

需要指出以下问题。

(1)其中凸多边形的生成及相应的交点排列顺序是利用MATLAB中DelaunayTri函数和convexhull函数得到的，需要注意的是，根据这两个函数的使用规则，并不能够处理这些交点的三维坐标(x,y,z)，因此需要将每个交点对应的二维裂缝面上的局部坐标(u,v)代入函数中处理得到凸多边形，这也说明将二维裂缝面参数化从而只需要用两个参数u、v去刻画裂缝面的处理是十分必要的。

(2)如果一个基质网格内包含某裂缝面交点的数量小于3，则无法构成一个凸多边形，则认为该裂缝面没有属于这个基质网格中的预裂缝网格。

按照上述方法，可以求出所有裂缝面上的预裂缝网格分布，其中每个预裂缝网格是该凸多边形顶点编号按照顶点顺序组成的有序数组来表征和存储的。

2.4 求解不同裂缝面之间的交线方程

由上述算法可知，预裂缝网格分布仅是基质网格与每一个单独裂缝面相交形成的裂缝网格分布，而实际情况中很可能还存在多个裂缝面相交，继而将预裂缝网格分布继续剖分的情况，因此接下来首先是要判断各个裂缝面之间是否会相交，如果相交，则需要求解交线的方程。具体算法如下。

不妨假设是对裂缝面1和裂缝面2进行处理，其裂缝面参数化形式分别如下。

裂缝面1：

r=α0+uα1+vα2，f1(u,v)≤0

(6)

式(6)中：α0=(α01,α02,α03)；α1=(α11,α12,α13)；α2=(α21,α22,α23)。

裂缝面2：

r=β0+aβ1+bβ2，f2(u,v)≤0

(7)

式(7)中：β0=(β01,β02,β03)；β1=(β11,β12,β13)；β2=(β21,β22,β23)。

通过计算裂缝面1和裂缝面2的法向量n1=α1×α2和n2=β1×β2，若n1×n2=0，则裂缝面1与裂缝面2平行，即这两个裂缝面之间不会相交；若n1×n2≠0，则裂缝面1与裂缝面2不平行，即裂缝面1所在平面和裂缝面2所在平面(此处是裂缝面所在平面，而不是裂缝面)必会存在交线。然后计算这两个裂缝面所在平面的交线方程，提出的方法如下：

首先证明一个引理。

引理如图4所示，若裂缝面1和裂缝面2不平行，则直线r1=α0+uα1和r2=α0+vα2中至少存在一条直线与裂缝面2所在平面存在唯一交点。

证明若直线r1=α0+uα1和r2=α0+vα2均与裂缝面2所在平面没有交点或者有无穷多个交点，则说明这两条直线的方向向量α1和α2均平行于裂缝2所在平面，又因为α1和α2线性无关，故推出裂缝面1所在平面与裂缝面2所在平面平行，产生矛盾，因此假设不成立，即直线r1=α0+uα1和r2=α0+vα2中至少存在一条直线与裂缝面2所在平面存在唯一交点。

图4 引理证明过程示意图

基于上述引理，如图5所示，不妨设直线r1=α0+uα1与裂缝面2所在平面存在唯一交点，结合裂缝面2的参数化形式得到方程α0+uα1=β0+aβ1+bβ2存在唯一解，即下述线性方程组存在唯一解：

(8)

图5 交线方程计算示意图

求解上述线性方程组得到即可得到(u,a,b)，继而计算出裂缝面1所在平面上直线r1=α0+uα1与裂缝面2所在平面的交点，记作点A。将直线r1=α0+uα1沿α2方向平移得到直线r3=(α0+α2)+uα1，可知直线r3仍是裂缝面1所在平面上的直线并与直线r1平行，因此直线r3与裂缝面2所在平面也存在唯一交点，即方程(α0+α2)+uα1=β0+aβ1+bβ2也存在唯一解，类似得解如下线性方程组即可计算出直线r3与裂缝面2所在平面的交点，记作点B。

(9)

点A和点B都是裂缝面1所在平面和裂缝面2所在平面之间的交点，又因为两点确定一条直线，故AB连线所在直线即是裂缝面1所在平面与裂缝面2所在平面之间的交线，从而可由A、B这两点的坐标确定出交线的方程。

同理，通过上述算法可以求解出各裂缝面所在平面之间的交线方程。

2.5 确定最终的裂缝网格分布

以裂缝面1为例说明具体的过程，若某裂缝面(称为裂缝面2)所在平面与裂缝面1所在平面存在交线且交线方程为r=γ0+mγ1。接下来就是要判断该交线会不会继续剖分某些预裂缝网格，具体算法如下。

以裂缝面1上的某一个预裂缝网格(称为预裂缝网格1)为例说明，假设该预裂缝网格是一个凸四边形，记作ABCD，并且ABCD是其构成凸多边形的连点顺序。

步骤1首先判断边AB、BC、CD、DA是否与交线平行，若平行，则说明该边与交线没有交点或者该边包含在交线内，这种情况下无需任何操作，若不平行，则该边所在直线与该交线必有唯一交点E，求解出该交点E，并判断点E是否在该边上(因为该边是一条有限长线段)，若在该边上，则存储交线与该边的交点E。

步骤2若存储的交点数小于2，则说明该交线对该“预裂缝网格”没有影响，无需继续剖分；若存储的交点数等于2(不可能大于2，因为“预裂缝网格是凸多边形”)，则需要根据存储的交点所在的具体位置将该“预裂缝网格”继续剖分成两个子凸多边形裂缝网格。以图6为例，凸四边形ABCD是某裂缝面上的一个预裂缝网格，交线与该网格相交于E、F两点，故将ABCD剖分为AEFD和BCFE这两个裂缝网格。

图6 某预裂缝网格被交线剖分为两个裂缝网格

步骤3对裂缝面1上的所有“预裂缝网格”按照步骤1、步骤2循环，可以得到裂缝面2所在平面与裂缝面1所在平面的交线r=γ0+mγ1对裂缝面1上“预裂缝网格”继续剖分后的网格分布。

步骤4对与裂缝面1所在平面存在交线的裂缝面循环步骤3，则可以得到裂缝面1上最终的裂缝网格分布。

同理可以得到其他裂缝面上最终的裂缝网格分布。

如图7(a)所示，两条裂缝相交后，两条裂缝面上的“预裂缝网格分布”被交线进一步剖分，图7(b)、图7(c)分别是两条裂缝面上最终裂缝网格分布的二维视角图。

图7 两条裂缝相交后裂缝网格分布

2.6 确定网格之间的连接关系

在EDFM中，存在四类网格之间的连接，分别是：①两相邻基质网格之间的连接；②同一裂缝面上且位于不同基质网格内的两相邻裂缝网格之间的连接；③同一基质网格内裂缝网格之间的连接；④基质网格与其所包含的裂缝网格之间的连接。

2.6.1 第一类连接

仅涉及相邻基质网格之间的连接，是比较平凡的，剩余三类均与裂缝网格有关，2.6节针对这三类连接分别简要给出确定相应连接关系的方法。

2.6.2 第二类连接

基于每个裂缝面上的裂缝网格分布，其中，每个裂缝网格是由该凸多边形网格顶点编号按照顶点顺序组成的有序数组来表征的，记这个有序数组为该裂缝网格的表征数组。对于同一裂缝面上的两个裂缝网格(记为e1、e2)，首先计算这两个裂缝网格表征数组中相同顶点编号的数量(记为n)，若n<2，则e1和e2之间没有公共边，即e1和e2无连接；由于都是凸多边形裂缝网格，因此不可能出现n>2的情况；若n=2，不妨设e1和e2公共的两个顶点分别是A、B,则e1和e2之间具有公共边AB，同时可以计算出该公共边长度(记为lAB)，e1和e2中心到AB的距离(分别记为de1,AB，de2,AB)，并可以计算后续块中心有限体积离散格式中传导系数所需的一些几何参数值。

2.6.3 第三类及第四类连接

在第三类连接中，有多个裂缝网格在同一个基质网格中，这种情况是由于有不同的裂缝面在该基质网格中相交造成的。在三维基质网格和二维裂缝面情况下，该基质网格与包含的裂缝网格之间的几何关系远比二维基质网格和一维裂缝线段的情况复杂得多，缺乏一个绝对精确的确定裂缝网格之间连接关系及相应几何因子计算的方法。对此，假设在同一基质网格中的所有裂缝网格之间均存在连接(在实际地质条件下，这种假设也是合理的，因为基质网格尺寸一般不大)。基于上述假设，由于2.3、2.5节中的算法可以获取每个基质网格内包含的裂缝网格编号，因此便可以建立在同一基质网格中的这些裂缝网格之间的连接，同时也可以建立该基质网格与所包含的每个裂缝网格之间的连接。

3 数值实例

根据前处理算法，建立能够适用于任意形状倾斜缝的三维EDFM，并给出一个三维油藏复杂缝网情况下单相原油衰竭开发的数值实例。

如图8所示，油藏工区大小为1 000 m×700 m×20 m，基质系统网格的总数是100×70×2，每个基质网格的尺寸为10 m×10 m×10 m。有一口多段压裂水平井定井底流压(10 MPa)生产，体积压裂形成的有效缝网包含10条裂缝，其中5条人工主裂缝和5条次生裂缝，其中有两条椭圆形主裂缝。可以看到这10条裂缝具有不同程度的倾斜(红色网格表示形成的裂缝网格)，且在纵向上贯穿的层位各有不同。油藏及流体的基本物性参数如表1所示。