基于三维嵌入式离散裂缝模型的致密油藏体积压裂水平井数值模拟

赵国翔 姚约东 王 链 唐 康 吴文炜

（1.中国石油大学（北京）油气资源与探测国家重点实验室，北京 102249；2.宁波锋成先进能源材料研究院有限公司，浙江 宁波 315000）

0 引 言

致密油藏的高效开发是中国油气勘探开发过程中的重点工作，而水平井体积压裂则是动用这类油藏的关键技术［1-9］。体积压裂施工后形成的复杂缝网为实现致密油藏商业化开发提供了可能［10-12］。体积压裂致密油藏复杂裂缝网络数值模型的建立难度较大，虽然目前已经建立了多种针对致密油藏缝网表征及数值模拟的方法［13-16］，但依然存在模型维度低、裂缝形态简单、实际生产过程中三维油藏复杂裂缝网络的分布特征难以刻画、与真实油藏中的复杂地质条件适应性差等问题，导致致密油藏开发动态模拟及产能预测难度较大。

国内外学者为实现复杂裂缝性油藏的数值模拟所建立的模型包括连续介质模型、等效连续介质模型及离散裂缝模型［14-17］。其中，连续介质模型及等效连续介质模型将裂缝及基质看作连续体，将裂缝渗透率平均到整个基质。等效连续介质模型精度较低，难以准确表征压裂缝形态及流体在裂缝中的特殊流动状态，处理大尺度裂缝时会低估其导流能力，增加计算误差［18-21］。离散裂缝模型（DFN）将裂缝显式处理［22-25］，通过生成与基质相匹配的裂缝网格对复杂裂缝形态进行表征，具有较高的计算精度。但由于裂缝形态分布复杂，DFN需要通过非结构化网格进行网格划分，网格剖分过程复杂，计算量大，裂缝面附近难以生成高质量的网格，导致模型计算效率低，收敛性及稳定性较差。嵌入式离散裂缝模型（EDFM）采用结构化的基质网格，将裂缝嵌入到正交基质网格之中，将裂缝处理为基质网格中的源汇项，避免了生成高质量匹配型非结构化网格的困难，保证了前处理的效率，也具备高于连续介质模型的精度，在复杂裂缝表征应用中具有明显优势［26-31］。但在嵌入式离散裂缝模型应用过程中，仍难以对三维空间中具有复杂缝网分布形态的三维裂缝进行表征，对于三维油藏模型中裂缝的立体分布、裂缝斜角、裂缝方位角及复杂裂缝形态仍难以准确刻画。

本文基于三维笛卡尔空间中的裂缝面参数化表征方法，形成了适用于任意形态倾斜裂缝的三维嵌入式离散裂缝模型网格剖分方法，建立了考虑致密油藏应力敏感特征及非线性渗流特征的三维嵌入式离散裂缝模型。采用与Eclipse 计算的结果对比并通过拟合新疆油田实际生产数据的方法验证模型的准确性与可靠性，对致密油藏进行体积压裂的产能影响因素进行了分析。

1 基于三维笛卡尔网格的嵌入式离散裂缝网格剖分方法

1.1 裂缝面的参数化表征

在三维欧式空间中，任意裂缝面的大小和位置，均可由裂缝面中心点位置坐标A（x0，y0，z0），裂缝面边界约束条件及裂缝面法向量n唯一确定，裂缝面法向量决定了裂缝面的方位角和倾角，对应的裂缝面方程可表示为

式中：l——裂缝面法向量的横坐标，m；x——裂缝面上任一点的横坐标，m；x0——裂缝面中心点的横坐标，m；m——裂缝面法向量的纵坐标，m；y——裂缝面上任一点的纵坐标，m；y0——裂缝面中心点的纵坐标，m；n——裂缝面法向量的竖坐标，m；z——裂缝面上任一点的竖坐标，m；z0——裂缝面中心点的竖坐标，m。

应用向量形式，则裂缝面方程可表示为

式中：n——裂缝面法向量；P——裂缝面上任一点；A——裂缝面中心点。

法向量n可由裂缝面的倾角α及方位角β表征，其中倾角α为裂缝面与x-y平面所成夹角，方位角β为裂缝所在平面与x-y平面交线与x轴正方向所成夹角，如图1所示。裂缝面法向量n的坐标计算方法为：

图1 裂缝面倾角α及方位角β与法向量n的关系Fig.1 Relation among dip angle α,azimuth angle β,and normal vector n

式中：α——裂缝面倾角，（°）；β——裂缝面方位角，（°）。

裂缝面方程可表示为

在三维欧式空间中，任意平面可进一步由该平面上的一点A(x0，y0，z0)及平面上两个线性无关的基向量表征，即

式中：u——裂缝面上的基向量；v——裂缝面上与u正交的基向量；a——基向量u对应的系数；b——基向量v对应的系数。

在已知平面法向量n时，u和v可利用平面法向量与坐标向量的叉积求取。裂缝边界形状可由长轴半径、短轴半径及相应的边界形状方程表征。

1.2 网格剖分方法

在三维EDFM 模型中，其网格剖分增加了求解面与面之间交线的过程，并要求识别裂缝网格面的复杂连接关系。本文三维EDFM 前处理算法的基本思路分为3 个步骤：（1）基于基质网格分布的初次网格剖分；（2）基于裂缝面交线的网格二次剖分；（3）确定网格间的连接关系。

1.2.1 基于基质网格对裂缝面进行初次网格剖分

在三维嵌入式离散裂缝模型中，基质部分利用三维笛卡尔网格进行剖分，形成的基质网格将与其相交的裂缝面剖分，其剖分节点可根据基质网格线与裂缝面的相互关系确定。

以基质网格两相邻节点Pa、Pb为节点的直线方程可表示为

式中：γ——直线方程中的比例系数；Pa，Pb——基质网格中的两相邻节点。

取网格线上任一已知点代入裂缝面方程，可判断网格线是否在裂缝面上。若网格线在裂缝面上，则取基质网格相邻节点Pa、Pb及网格线与裂缝面边界交点Pm、Pn为网格面剖分节点。若网格线与裂缝面相交，其交点可表示为

以上方法计算得到的基质网格线与裂缝面的交点有可能位于裂缝面形状约束范围外，因此应用面积判别法对所得交点是否位于裂缝面形状约束范围内进行判断。根据剖分域内一点与各节点组成的向量在由u和v构成的局部坐标系中与方向向量的夹角对各节点进行排序（图2），形成凸多边形网格。

图2 裂缝面剖分域内各节点逆时针排序结果Fig.2 Counter-clockwise sort of nodes in fracture surface subdivision domain

对于裂缝面边界节点，需要对其是否在目标网格内进行判断，如图3（a）所示。基质网格的六面体包含3 组平行平面，当目标点P与2 个平行平面向外的法向量的夹角为一个锐角和一个钝角时，则点P位于2 个平行平面内，当目标点P与2个平行平面向外的法向量的夹角均为锐角或均为钝角时，则点P位于两平行平面外，若点P均在3组平行平面内，则点P位于基质网格内。当裂缝位于基质网格内或未与基质网格线形成有效交点时，需对此类裂缝进行识别（图3（b）、（c））并进行网格剖分。当裂缝面完全处于基质网格内时，其网格剖分结果则可直接应用裂缝面节点信息获得，如图3（c）所示。

图3 裂缝面网格剖分过程中的特殊情况Fig.3 Special cases in mesh generation of fracture surface

在获得目标裂缝面所有剖分节点后，将各基质网格内的剖分节点以逆时针方向进行连接形成凸多边形，可得到基于基质网格的裂缝面初步剖分结果，如图4所示。

图4 基于基质网格的裂缝面初次网格剖分结果Fig.4 Initial mesh generation of fracture surface based on matrix grid

1.2.2 基于裂缝面交线进行网格的二次剖分

前述网格剖分过程仅考虑了基质网格对裂缝面的剖分，但实际情况下，多条裂缝可能相交，初次裂缝面网格剖分结果会被相交的裂缝进一步剖分。在两个相交的裂缝面（法向量分别为n1和n2）上，过正交坐标向量u和v的直线，必有一条直线与另一裂缝相交。假设过平面1 上中心点的直线方程与平面2 相交，其交点P1可表示为

式中：P1——过平面1中心点的直线与平面2的交点；A1——平面1过中心点且与基向量u平行的点；n2——平面2 的法向量；A2——平面2 的中心点。

平面1 和平面2 的交线方程l1可表示为

式中：l1——平面1 和平面2 的交线；n1——平面1的法向量。

在求得交线方程l1的基础上，根据裂缝面网格初次剖分所得的凸多边形，首先判断交线是否与各边平行，若不平行，则计算交线与目标边的交点，并判断此点是否处于目标边上，若此交点处于目标边内，则保存此交点作为裂缝面二次剖分的节点。对同一凸多边形内，根据存储交点及初次剖分的凸多边形节点，对目标凸多边形进行二次剖分。考虑到可能存在多个裂缝面相交（图5），则对各个裂缝面的相互关系、网格二次剖分过程循环处理，从而得到最终的网格剖分结果。

图5 基于裂缝面交线进行的最终网格剖分结果Fig.5 Final mesh generation of fracture surface based on intersecting line

1.2.3 网格间连接关系的确定

以图5为例，各网格间存在4 种连接关系：（1）基质网格与基质网格之间的连接；（2）基质网格与基质网格内裂缝网格的连接；（3）同一裂缝面内各网格之间的连接；（4）同一基质网格内相交裂缝网格之间的连接。由于每个节点均包含三维笛卡尔坐标系内的坐标信息，因此各网格的面积、体积、公共边长度等几何参数均可计算得到。

2 三维嵌入式离散裂缝模型求解方法

2.1 渗流控制方程

单相原油基本渗流微分方程为

式中：K——渗透率，10-3μm2；μ——黏度，mPa·s；B——体积系数；∇p——流体压力梯度，MPa/ m；p——压力，MPa；g——重力加速度，m/ s2；ρ——地面标况下的密度，g/ cm3；q——地面标况下的体积流量（源汇项），m3/s；V——网格体积，m3；δ——狄拉克函数；ϕ——孔隙度，%；Z——垂向高度，m；t——时间，s。

致密油藏应力敏感特征，本文采用常用的指数关系模型［32］，其表达式为：

式中：K0——初始油藏压力下的渗透率，10-3μm2；pi——原始地层压力，MPa；CK——渗透率应力敏感系数，MPa-1；Cϕ——孔隙度应力敏感系数，MPa-1；ϕ0——初始油藏压力下的孔隙度，%。

基质非线性渗流模型所采用模型［33］为

式中cl——非线性参数，MPa/ m。

井模型与边界条件。当井在裂缝网格上时，考虑裂缝网格所在平面可能与压裂水平井所在直线夹角及裂缝形状的不规则性，井指数计算公式为：

式中：WI——井指数，10-3m3/（GPa·s）；bf——裂缝缝宽，m；θ——水平井所在直线与裂缝平面法向量的夹角，(°)；re——等效供给半径，m；rw——井筒半径，m；S——水平井穿过的裂缝网格的面积，m2。

2.2 差分方法及传导率计算

本文采用满足局部质量守恒且各参数物理意义明确的块中心有限体积法对控制方程进行离散求解，采用封闭边界条件，其离散结果为：

（1）基质系统，其方程式可表达为

（2）裂缝系统，其方程式可表达为

式中：T——网格间的传导率，m3；n——时间步序号；Δt——时间步长，s；下标m 和f 分别表示基质系统和裂缝系统；下标i和j为两相邻网格的编号。

模型离散方程组可以表示为F(x) = 0，应用牛顿迭代法［27］及自动微分方法［34］求解此方程组。

考虑网格剖分后出现的4 种网格连接关系，其传导率计算公式为［27，35］：

式中：Ag——相邻网格接触面积，m2；d——网格中心点到接触面中心点的距离，m。

3 模型验证

以Eclipse 模型及新疆油田JT2-H 水平井生产数据为例［36］，对本文三维嵌入式离散裂缝模型的准确性进行了验证，模拟参数如表1所示。Eclipse模型所用参数中油藏规模（长×宽×高）为400 m×400 m×20 m，裂缝数量为1 条，忽略基质非线性渗流特征，其余与表1所示参数相同。从图6可以看出，本文与Eclipse 模拟结果基本一致，并实现了与实际生产数据的拟合。

表1 三维嵌入式离散裂缝数值模拟参数Table 1 Parameters for numerical simulation by 3D EDFM

图6 本文模型与Eclipse模型验证Fig.6 Validation of this model and Eclipse

4 数值模拟实例

在模型验证结果的基础上，分别进行了裂缝形态、裂缝方位角、倾角、天然裂缝及基质非线性渗流特征影响的模拟计算。油藏规模（长×宽×高）为1 400 m×800 m×60 m，原始地层压力30 MPa，井底流压25 MPa，主压裂缝7 条，其余参数与模型验证中相同。

4.1 裂缝形态的影响

对于不同的裂缝形态，分别进行了常规压裂缝、非贯穿缝、不等长缝、交错缝、椭圆缝以及复杂缝网6 种压裂缝形态的对比。在传统模型中，复杂缝网可利用双重介质模型进行等效处理，本文采用7 条压裂缝，28 条次级裂缝，200 条形状、倾角、方位角、尺寸各异的天然缝进行复杂缝网的表征。

常规缝及复杂缝网的网格剖分结果如图7（a）、（b）所示。由图7（c）—（h）可知，随着压裂缝控制面积的增加，压裂缝区域受到的影响随之增大，对于交错分布的裂缝（图7（f）），虽然交错分布方式增加了压裂缝在纵向上的控制体积，但降低了压裂缝在横向上的控制体积，导致两条相邻压裂缝之间存在压力波未波及区域。因此，在致密油藏压裂施工过程中，增加裂缝控制面积的同时，也应对压裂缝间距进行优化。

图7 不同裂缝形态网格剖分及定压生产100 d的压力场Fig.7 Mesh generation and pressure field of 100 d constant pressure production of different fracture geometry

裂缝形态在影响开发过程中储层压力分布的同时，也导致压裂水平井产能随着裂缝形态的改变而受到显著影响（图8）。在100 d 的定压生产过程中，裂缝复杂程度越高，改造体积越大，其产能越大。并且裂缝形态对产能的影响主要体现在早期线性流阶段，在后期拟径向流阶段，裂缝形态对产能的影响逐渐减小。

图8 不同形态裂缝影响的产能随时间的变化Fig.8 Change of productivity with time affected by different fracture geometry

4.2 天然裂缝的影响

在天然裂缝性致密油藏中，天然裂缝发育，分布形式复杂。图9中共包含200 条方位角分别为45°或-45°的天然裂缝、28 条次级裂缝以及7 条主压裂缝；天然裂缝宽度为0.5 mm，渗透率为50×10-3μm2，其余参数与模型验证中参数相同，其裂缝分布、网格剖分及模拟结果如图9所示。由于大量天然裂缝的存在，油藏压力分布变化显著，在主裂缝与天然裂缝沟通的区域，流体经过天然裂缝所在基质网格流入裂缝，导致此基质压力与裂缝压力明显小于其余区域，而对于不与裂缝连通的天然裂缝，压力降低幅度较小。同时，由图9（b）、（d）可以看出，当存在天然裂缝时，累计产量更高，说明天然裂缝对于产能有积极作用，因此，致密油藏压裂过程中应选择天然裂缝发育的区域，并使主压裂缝尽可能与天然裂缝沟通，形成控制体积更大的复杂缝网。

图9 天然裂缝分布、网格剖分及模拟结果Fig.9 Natural fractures distribution,mesh generation and numerical simulation result

4.3 裂缝方位角及倾角的影响

为分析裂缝方位角及倾角的影响，分别进行了压裂缝方位角及倾角为15°、30°、45°、60°、75°、90°的模拟计算，各压裂缝面积相同。由图10可以看出，随着裂缝方位角的逐渐减小，压裂缝控制区域逐渐重叠并导致沿水平井长度方向的压力波及范围逐渐减小，压裂缝之间的干扰变大，使得产能逐渐降低。对于等间距分布的压裂缝，小倾角裂缝在重力作用下使得累计产量有所增加，但压裂改造体积并未随着裂缝倾角的变化而明显增加，因此裂缝倾角对累计产油量的影响较小可以忽略不计（图11）。

图10 不同方位角裂缝影响的产能随时间的变化Fig.10 Change of productivity with time affected by different fracture azimuth

图11 不同倾角裂缝影响的产能随时间的变化Fig.11 Change of productivity with time affected by different fractures dips

4.4 基质非线性渗流特征的影响

致密油藏基质渗透率低，孔隙结构复杂，包含大量的微纳米孔隙，具备显著的微尺度效应及非线性渗流特征。基于以上特点，分别选择非线性渗流参数为0.05、0.5、5 MPa/m 来表征弱非线性、中非线性及强非线性渗流特征对产能的影响，如图12 所示。在衰竭式开发过程中，由于非线性渗流特征产生额外的渗流阻力，压力波及范围随非线性程度的增加而降低，导致累计产油量显著下降。

图12 不同非线性参数影响的产能随时间的变化Fig.12 Change of productivity with time affected by different nonlinear parameters

5 结 论

（1）针对传统嵌入式离散裂缝模型难以对具有复杂形态分布的裂缝进行数值刻画的问题，建立了三维笛卡尔坐标系中裂缝面的参数化表征方法及可应用于任意形状倾斜裂缝面的三维嵌入式离散裂缝模型网格剖分方法。并考虑重力、基质非线性渗流特征及应力敏感特征建立了三维嵌入式离散裂缝模型，实现了具有复杂裂缝形态及分布特征的体积压裂致密油藏的开发动态模拟。

（2）裂缝形态对致密油藏产能有明显影响，压裂缝控制体积越大，产能越高。与主裂缝连通的天然裂缝对油藏开发过程中的压力分布有明显影响，非连通天然裂缝压力降低程度较低，对产能的影响较连通天然裂缝小。同时，裂缝倾角对产能影响较小，裂缝方位角对产能的影响较明显。因此，在体积压裂致密油藏开发动态模拟过程中，可利用本文所建立的三维嵌入式离散裂缝模型对具有复杂形态及分布特征的天然裂缝进行表征，并对人工压裂缝形态及分布特征进行对比优化，为压裂设计及产能优化提供参考。

（3）基质非线性渗流特征对致密油藏产能有明显影响，通过实验或模拟计算等方法合理选择基质非线性渗流参数对体积压裂致密油藏的开发动态模拟有重要意义。