CN104112057A - 一种大尺度裂缝性油藏数值模拟方法 - Google Patents

Abstract

本发明提供了一种大尺度裂缝性油藏数值模拟方法，属于油气田开发领域。所述方法将油藏中的大尺度裂缝作为油藏的内边界进行降维处理，建立油藏几何模型，然后采用三角形网格对油藏几何模型进行几何剖分形成离散单元，同时形成油藏几何信息描述文件、几何离散信息文件及井信息描述文件，在离散单元基础上构建数值计算单元，建立表征此类油藏流体流动的数学模型，采用积分法在每个数值计算单元上建立其数值计算格式，形成大型线性代数方程组，并进行求解，从而实现大尺度裂缝性油藏数值模拟。本发明能够精确表征大尺度裂缝性油藏，物理意义明确，对大尺度裂缝油藏数值模拟具有很好的适应性。

Description

一种大尺度裂缝性油藏数值模拟方法

技术领域

本发明属于油气田开发领域，具体涉及一种大尺度裂缝性油藏数值模拟方法。

背景技术

大尺度裂缝性油藏是指储层中存在开度大于1毫米裂缝的油藏。大尺度裂缝在油藏中起着特殊导流作用，其存在导致了储层的强烈非均质性，造成了油藏复杂的渗流机理和相对较低的采收率，为了高效开发此类油藏，需要采用数值模拟手段对其注水开发过程进行精确模拟，设计合理的开发方案，避免水沿大尺度裂缝窜进，达到提高采收率的目的。

目前，在大尺度裂缝性油藏数值模拟中，常采用常规连续性介质模型表征此类油藏，如：单孔隙介质模型，双重孔隙介质模型以及等效连续介质模型等。单孔隙介质模型能够较精确的模拟流体流动过程，但由于基岩大小与裂缝开度在尺度上的巨大差异，在进行数值计算时需要在裂缝周围进行大量的网格加密，大大降低了计算效率；双重孔隙介质模型将裂缝性油藏的基岩系统与裂缝系统看作两个相互交错的连续性多孔介质系统，两系统间通过窜流量系数进行联系，该模型适用于裂缝发育程度较高，连通性较好的网状裂缝油藏，对于裂缝发育程度不高或存在主导流体流动的大尺度裂缝时，数值模拟结果与实际相差较大；等效连续介质模型本身并不复杂，但该模型的有效性难以判定，且裂缝性油藏的等效渗透率张量难以求解。采用离散介质模型表征大尺度裂缝性油藏，能够真实地代表裂缝的地质分布特征，基于有限元法可以较好地进行该类油藏的数值模拟，但对于严重非均质性油藏中水驱油等多相流问题不能保证局部质量守恒，不适合描述基岩-裂缝交界面处饱和度的不连续性。

综上所述，当前各种常规大尺度裂缝油藏的数值模拟方法都存在不同的问题：基于单孔隙介质模型的方法，存在计算效率低的问题；基于双重孔隙介质模型的方法，仅适用于裂缝发育程度较高，连通性较好的网状裂缝油藏；基于等效连续介质模型的方法，其有效性难以判定，且裂缝性油藏的等效渗透率张量难以求解；基于离散裂缝模型能够较真实地刻画大尺度裂缝的地质分布特征，但采用有限元数值计算格式进行模拟，存在着计算量大，物理意义不明确，在含大尺度裂缝的油藏中不能严格保证局部质量守恒等问题。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种大尺度裂缝性油藏数值模拟方法，高效精确地进行大尺度裂缝性油藏的数值模拟，为分析裂缝对注水开发效果的影响，优选布井方式，避免注入水沿裂缝窜进，预测提高采收率措施的效果提供可靠依据，同时物理意义明确，能够严格保证局部质量守恒。

本发明是通过以下技术方案实现的：

一种大尺度裂缝性油藏数值模拟方法，所述方法将油藏中的大尺度裂缝作为油藏的内边界进行降维处理，建立油藏几何模型，然后采用三角形网格对油藏几何模型进行几何剖分形成离散单元，同时形成油藏几何信息描述文件、几何离散信息文件及井信息描述文件，在离散单元基础上构建数值计算单元，建立表征此类油藏流体流动的数学模型，采用积分法在每个数值计算单元上建立其数值计算格式，形成大型线性代数方程组，并进行求解，从而实现大尺度裂缝性油藏数值模拟。

所述方法具体包括以下步骤：

(1)根据实际油藏特征及岩心、测井及流体的测试资料确定油藏储层物性参数，建立表征大尺度裂缝性油藏流体流动的数学模型；

(2)依据油藏实际情况及油藏描述成果，确定油藏几何边界范围，建立油藏几何模型，形成油藏几何信息描述文件；

(3)对步骤(2)建立的油藏模型进行几何离散，生成油藏模型的几何离散信息文件：平面上基岩部分采用Delaunay三角形单元进行划分，大尺度裂缝用线单元代表，纵向上采用储层厚度代表单元厚度，生成油藏模型的几何离散信息文件；

(4)基于步骤(3)的几何离散信息文件，连接所有三角形的重心与三边的中点形成联络线，以三角形的顶点为中心，其周围的所有联络线组成的多边形即为数值计算单元；

(5)计算数值计算单元节点参数值：设定饱和度变量(S_w，S_o)在每个数值计算单元内恒定不变，其数值通过测井数据取得，而数值计算单元的流动势变量(φ_w，φ_o，φ_c)由以该数值计算单元的节点为顶点的相邻三角网格的值通过线性逼近计算得到：

$\mathrm{\Φ}\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{N}_i\left(x\right){\mathrm{\Φ}}_i---\left(7\right)$

其中，n为三角形单元的顶点数；φ_i代表节点i的任意流动势变量；x为节点i的坐标；N_i(x)为形函数，为坐标的独立函数；

由式(7)可得在一个三角形内任意变量的梯度恒为：

$\▿\mathrm{\Φ}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\Φ}}_i\▿N_i\left(x\right)---\left(8\right)$

(6)对步骤(1)建立的数学模型在步骤(4)建立的数值计算单元上进行积分，建立流动控制方程的数值计算方程，数值计算单元节点参数值是由步骤(5)计算得到的；

(7)确定油藏的边界条件与初始条件，然后根据步骤(6)建立的每一个数值计算单元上的数值计算方程，构建线性代数方程组，并进行求解，完成大尺度裂缝性油藏的数值模拟。

所述步骤(1)中建立的表征大尺度裂缝性油藏流体流动的数学模型如下：

${\∫\∫}_{\mathrm{\Ω}}\mathrm{FEQd\Ω}={\∫\∫}_{{\mathrm{\Ω}}_m}\mathrm{FEQd}{\mathrm{\Ω}}_m+e\×{\∫}_{{\mathrm{\Ω}}_f^{\′}}\mathrm{FEQd}{\mathrm{\Ω}}_f^{\′}---\left(6\right)$

其中FEQ代表流动控制方程式(1)和式(5)：

$\frac{\∂\left(\mathrm{\φ}{S}_w\right)}{\∂t}-\▿\·({\mathrm{\λ}}_w\▿{\mathrm{\Φ}}_w)-q_w=0---\left(1\right)$

$\▿\·(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\▿{\mathrm{\Φ}}_w)+\▿\·({\mathrm{\λ}}_o\▿{\mathrm{\Φ}}_c)+(q_o+q_w)=0---\left(5\right)$

式中：油藏整体区域用Ω表示，其中基岩部分用Ω_m代表，简化降维的大尺度裂缝部分用Ω′_f表示，e表示裂缝的开度；w，o分别表示水相和油相；φ为油藏孔隙度；S_w为油藏含水饱和度；S_o为油藏含油饱和度；为哈密顿算子；q_w，q_o为源汇项；λ_w为油藏水流度；λ_o为油藏原油流度；Φ_w为油藏水流动势；Φ_o为油藏原油流动势；φ_c为毛管压力势；p_c为毛管压力；ρ_w为油藏水的密度；ρ_o为油藏原油的密度；g为重力加速度；z为油藏垂向坐标，向上为正。

所述步骤(2)包括：

(21)建立外边界轨迹文件，该文件的内容包括：描述油藏边界轨迹的点的X、Y坐标值及储层深度值；

(22)将大尺度裂缝作为油藏内边界进行简化降维处理，根据油藏描述得到的大尺度裂缝具体位置，建立大尺度裂缝的具体信息文件，具体信息文件的内容包括；裂缝轨迹点对应的裂缝编号、X、Y坐标值、裂缝高度、裂缝开度；

(23)根据油藏中井的情况，建立井信息文件，该井信息文件的内容包括井号、X、Y坐标值、井所在的组别及井类别(生产井、注水井等)。

所述步骤(22)中的所述简化降维处理是指将三维油藏中的大尺度裂缝作为二维面单元，将二维油藏面上的裂缝作为一维线单元。

所述步骤(3)中的所述几何离散信息文件包括：基岩中三角形单元的序号、相邻三角形序号、三个顶点序号、三条边序号以及重心的坐标，裂缝中线单元的序号、两个顶点序号、中点坐标值。

与现有技术相比，本发明的有益效果是：

(1)本发明将油藏中的大尺度裂缝作为内边界进行简化降维处理，减少了网格划分的数量，采用积分的方式建立的数值计算格式，能够严格保证局部质量守恒，计算高效且稳定性强。

(2)本发明能够精确表征大尺度裂缝性油藏，物理意义明确，对大尺度裂缝油藏数值模拟具有很好的适应性，与常规数值模拟计算结果一致，计算效率提高2.5倍以上，而且可操作性强，有效实用。

附图说明

图1是本发明中大尺度裂缝油藏平面几何离散示意图。

图2是本发明中的数值计算单元的构成示意图。

图3是实施例1中大尺度裂缝油藏的物理模型。

图4(a)是依据本发明方法进行的实施例1中大尺度裂缝油藏的剖分示意图。

图4(b)是依据常规方法I进行的实施例1中大尺度裂缝油藏的剖分示意图。

图4(c)是依据常规方法II进行的实施例1中大尺度裂缝油藏的剖分示意图。

图5-1(a)是当注水总体积分别为0.5PV时，采用本发明方法计算得到的实施例1油藏物理模型的含水饱和度分布平面图。

图5-1(b)是当注水总体积分别为1.5PV时，采用本发明方法计算得到的实施例1油藏物理模型的含水饱和度分布平面图。

图5-1(c)是当注水总体积分别为3.0PV时，采用本发明方法计算得到的实施例1油藏物理模型的含水饱和度分布平面图。

图5-2(a)是当注水总体积分别为0.5PV时，采用常规方法I计算得到的实施例1油藏物理模型的含水饱和度分布平面图。

图5-2(b)是当注水总体积分别为1.5PV时，采用常规方法I计算得到的实施例1油藏物理模型的含水饱和度分布平面图。

图5-2(c)是当注水总体积分别为3.0PV时，采用常规方法I计算得到的实施例1油藏物理模型的含水饱和度分布平面图。

图6是基于本方法和常规方法I计算的实施例1油藏物理模型中不同裂缝倾角下采收程度随时间的变化对比图。

图7是实施例2中大尺度裂缝油藏的物理模型。

图8是依据本方法进行的实施例2中大尺度裂缝油藏的剖分示意图。

图9-1(a)是当注水时间为150天时，采用本发明方法计算得到的实施例2油藏物理模型的含水饱和度分布平面图。

图9-1(b)是当注水时间为300天时，采用本发明方法计算得到的实施例2油藏物理模型的含水饱和度分布平面图。

图9-1(c)是当注水时间为600天时，采用本发明方法计算得到的实施例2油藏物理模型的含水饱和度分布平面图。

图9-2(a)是当注水时间为1200天时，采用本发明方法计算得到的实施例2油藏物理模型的含水饱和度分布平面图。

图9-2(b)是当注水时间为2400天时，采用本发明方法计算得到的实施例2油藏物理模型的含水饱和度分布平面图。

图9-2(c)是当注水时间为4800天时，采用本发明方法计算得到的实施例2油藏物理模型的含水饱和度分布平面图。

图10是基于本发明方法计算的实施例2油藏物理模型中含大尺度裂缝情况下和均质油藏情况下采收程度随时间的变化对比图。

图11是基于本发明方法计算的实施例2油藏物理模型中含大尺度裂缝情况下和均质油藏情况下含水率随时间的变化对比图。

图12是本发明大尺度裂缝性油藏数值模拟方法的步骤框图。

具体实施方式

下面结合附图对本发明作进一步详细描述：

如图1所示，一种大尺度裂缝性油藏数值模拟方法，对油藏中的大尺度裂缝进行简化降维处理，将大尺度裂缝看作油藏的内边界，利用非结构化网格对油藏进行网格划分(即几何离散)，并在此基础上构建数值计算单元。建立大尺度裂缝性油藏的数学模型，在每个数值计算单元上采用积分的方式构建流动控制方程的数值计算格式，进行大尺度裂缝性油藏的数值模拟。在保证高效计算的前提下，精确考虑了大尺度裂缝性油藏的非均质性，同时不必计算基岩和裂缝间的窜流函数，具有严格保证局部质量守恒，物理意义明确，计算效率高等优点。

图1中，三角形代表基岩单元，粗线段代表裂缝单元，二者纵向都具有一定厚度，三角形厚度为储层单层厚度，线段厚度为裂缝高度。

本发明是按照以下方式实现的：

对油藏中大尺度裂缝采取降维的处理方式，建立大尺度裂缝油藏的数学模型，采用Delaunay三角形对油藏进行几何离散，并以此为基础构建数值计算单元，在数值计算单元对流动控制方程上进行积分，形成其数值计算格式，实现大尺度裂缝性油藏的数值模拟(即基于积分的方法建立大尺度裂缝性油藏中油水两相流动控制方程的数值计算格式，进行大尺度裂缝油藏的数值模拟)。

具体实施方式如下：

选择需要进行数值模拟的含有大尺度裂缝的油藏区块，参照附图对本发明方法作详细描述，具体步骤如图12所示，包括以下步骤：

(1)根据实际油藏特征及岩心、流体等的测试资料确定油藏储层物性参数，建立表征大尺度裂缝性油藏油水两相流体流动的数学模型。假设经测试是不可压缩油藏，其油水两相流动控制方程为：

$\frac{\∂\left(\mathrm{\φ}{S}_w\right)}{\∂t}-\▿\·({\mathrm{\λ}}_w\▿{\mathrm{\Φ}}_w)-q_w=0---\left(1\right)$

$\frac{\∂\left(\mathrm{\φ}{S}_o\right)}{\∂t}-\▿\·({\mathrm{\λ}}_o\▿{\mathrm{\Φ}}_o)-q_o=0---\left(2\right)$

式中：w，o分别表示水相和油相；φ为油藏孔隙度；S_w为油藏含水饱和度；S_o为油藏含油饱和度；为哈密顿算子；q_w，q_o为源汇项；λ_w为油藏水流度；λ_o为油藏原油流度；Φ_w为油藏水流动势；Φ_o为油藏原油流动势。

定义油藏中：

S_w+S_o＝1    (3)

φ_c＝φ_o-φ_w＝p_c+(ρ_o-ρ_w)gz    (4)

式中：φ_c为毛管压力势；p_c为毛管压力；ρ_w为油藏水的密度；ρ_o为油藏原油的密度；g为重力加速度；z为油藏垂向坐标，向上为正。

将式(1)(2)相加，结合式(3)(4)，得出式(5)：

$\▿\·(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\▿{\mathrm{\Φ}}_w)+\▿\·({\mathrm{\λ}}_o\▿{\mathrm{\Φ}}_c)+(q_o+q_w)=0---\left(5\right)$

油藏中油水两相流动的控制方程可由式(1)和(5)表示，式(1)称为油藏含水饱和度方程，式(5)称为油藏流体流动势方程。

采用平行板模型和立方定律描述裂缝中的流体渗流，由N-S方程给出裂缝中渗流的一般形式，根据立方定律求得沿裂缝开度的流量以及速度分布，依据Darcy定律求取裂缝的等效渗透率，其沿裂缝开度方向是不变的，根据等效后流体流动参数及相关物理量沿裂缝开度方向的不变性，对裂缝进行简化降维处理建立大尺度裂缝性油藏流体流动数学模型。

本方法的大尺度裂缝性油藏，油藏整体区域用Ω表示，其中基岩部分用Ω_m代表，简化降维的大尺度裂缝部分用Ω′_f表示，e表示裂缝的开度，则大尺度裂缝性油藏流体流动数学模型可以表示为：

${\∫\∫}_{\mathrm{\Ω}}\mathrm{FEQd\Ω}={\∫\∫}_{{\mathrm{\Ω}}_m}\mathrm{FEQd}{\mathrm{\Ω}}_m+e\×{\∫}_{{\mathrm{\Ω}}_f^{\′}}\mathrm{FEQd}{\mathrm{\Ω}}_f^{\′}---\left(6\right)$

其中FEQ代表流动控制方程式(1)和式(5)。

(2)依据油藏实际情况及油藏描述成果，确定油藏的几何边界范围，建立如表1所示的外边界轨迹文件(即油藏边界范围文件)，该文件具体信息包括：描述油藏边界轨迹的点的X、Y坐标值与储层深度值；

表1

将大尺度裂缝作为油藏内边界进行简化降维处理，根据油藏描述得到的大尺度裂缝具体位置，建立如表2所示的大尺度裂缝的描述文件，文件具体内容包括；裂缝轨迹点对应的裂缝编号、X、Y坐标值、裂缝高度、裂缝开度；

表2

根据油藏中井的情况，建立如表3所示的井信息文件，文件内容包括井号、X、Y坐标值、井所在的组别及井类别(生产井、注水井等)。

表3

(3)对要进行数值模拟的大尺度裂缝油藏进行几何离散，基岩部分平面采用Delaunay三角形划分，三角形具有一定厚度，纵向上所有三角形厚度相加即为储层厚度，大尺度裂缝平面上用线单元代表，线单元具有一定宽度，代表裂缝真实高度。离散结果如图1所示，并生成油藏的离散信息文件。文件内容包括：基岩中三角形单元的序号、相邻三角形序号、三个顶点序号、三条边序号、重心的坐标以及三角形单元厚度，裂缝中线单元的序号、两个顶点序号、中点坐标值、裂缝高度及开度等；

(4)在上述网格划分的基础上构建数值计算单元，将三角形单元重心与三个边中点相连将单元划分为三个子面积，相应于每一个节点的表征区域是由此节点相连的所有子面积构成，是一个多边形面积。如图2所示，节点a的相邻节点为{b₁，b₂…b₆}，以节点a为顶点的各三角形为{T₁，T₂…T₆}，各三角形的重心为{G₁，G₂…G₆}，以节点a为顶点的各边中点分别为连接各三角形重心以及相应边中点，可得节点a的表征区域(即数值计算单元)为多边形其中，代表裂缝。

(5)定义饱和度变量(S_w，S_o)在每个数值计算单元内恒定不变，而流动势变量(φ_w，φ_o，φ_c)由每个Delaunay三角网格的值通过线性逼近估计：

$\mathrm{\Φ}\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{N}_i\left(x\right){\mathrm{\Φ}}_i---\left(7\right)$

其中，n为三角形单元的顶点数；φ_i代表节点i的任意流动势变量；x为节点i的坐标；N_i(x)为形函数，为坐标的独立函数。

由式(7)可得在一个三角形内任意变量的梯度恒为：

$\▿\mathrm{\Φ}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\Φ}}_i\▿N_i\left(x\right)---\left(8\right)$

(6)对大尺度裂缝性油藏油水两相流动控制方程在数值计算单元上进行积分，构建其数值计算格式。首先对含水饱和度方程(1)在任意数值计算单元CV_i上进行积分，即：

${\∫\∫}_{\mathrm{\Ω}}\frac{\∂\left(\mathrm{\φ}{S}_w\right)}{\∂t}\mathrm{dA}-{\∫\∫}_{\mathrm{\Ω}}\▿\·({\mathrm{\λ}}_w\▿{\mathrm{\Φ}}_w)\mathrm{dA}-{\∫\∫}_{\mathrm{\Ω}}{q}_w\mathrm{dA}=0---\left(9\right)$

设定孔隙度仅在空间位置上变化，储层厚度为单位厚度，并对式(9)第二项应用高斯散度定理将面积分转化为线积分，可得：

${\∫\∫}_{\mathrm{\Ω}}\mathrm{\φ}\frac{\∂S_w}{\∂t}\mathrm{dA}-{\∫}_{\mathrm{\Γ}}({\mathrm{\λ}}_w\▿{\mathrm{\Φ}}_w)\·\mathrm{nd\Γ}-{\∫\∫}_{\mathrm{\Ω}}{q}_w\mathrm{dA}=0---\left(10\right)$

式中，Γ指的是数值计算单元CV_i的边界；n是Γ上的外法向量。

方程(10)左边第一项积分可以近似表示为：

${\∫\∫}_{\mathrm{\Ω}}\mathrm{\φ}\frac{{\∂S}_w}{\∂t}\mathrm{dA}\≈A_{\mathrm{\φi}}\frac{\∂S_w^m}{\∂t}---\left(11\right)$

其中，

$A_{\mathrm{\φi}}=\underset{k=1}{\overset{t}{\mathrm{\Σ}}}{\mathrm{\ψ}}_k{A}_k{\mathrm{\φ}}_k^m+\underset{l=1}{\overset{s}{\mathrm{\Σ}}}\frac{d{S}_w^f}{d{S}_w^m}{e}_l\left|L_l\right|{\mathrm{\φ}}_l^f---\left(12\right)$

式(12)中A_φi为CV_i孔隙体积；右边第一项表示CV_i内基岩孔隙体积；第二项表示CV_i内裂缝孔隙体积。式中，t表示以节点i为顶点的Delaunay三角形单元总数；ψ_k为数值计算单元CV_i内三角形k的面积占所在Delaunay三角形K的面积比例；A_k为Delaunay三角形K的面积；指三角形K内基岩部分的孔隙度；s表示CV_i内裂缝的总条数；e_l和|L_l|分别为数值计算单元CV_i内第l条裂缝的孔隙度、开度与长度。

方程(10)左边第二项积分可以表示为：

${\∫}_{\mathrm{\Γ}}({\mathrm{\λ}}_w\▿{\mathrm{\Φ}}_w)\·\mathrm{nd\Γ}\≈\underset{k=1}{\overset{t}{\mathrm{\Σ}}}\left|s_k\right|{[{\mathrm{\λ}}_w^m\left(S_w^{m,\mathrm{up}}\right)\▿{\mathrm{\Φ}}_w]}_k\·n_k+\underset{l=1}{\overset{s}{\mathrm{\Σ}}}{e}_l{\mathrm{\λ}}_w^f\left(S_m^{f,\mathrm{up}}\right)\frac{\∂{\mathrm{\Φ}}_w^f}{\∂\mathrm{\ξ}}---\left(13\right)$

式(13)中右边第一项代表通过数值计算单元CV_i边界的流量；右边第二项表示通过CV_i内每条裂缝的流量。其中上角标up指上游饱和度；|s_k|为CV_i在三角形K内的边界，具有向外的单位法向量n_k；为|s_k|处的水相流动势梯度，可由式(8)近似计算；ξ表示沿裂缝方向的局部坐标；项代表裂缝l上的流动势梯度，裂缝上的流动可以看成一维的，因此可以通过式(14)近似计算。

$\frac{d{\mathrm{\Φ}}_w^f}{\mathrm{d\ξ}}=\frac{{\mathrm{\Φ}}_j-{\mathrm{\Φ}}_i}{2\left|L_l\right|}---\left(14\right)$

方程(10)的第三项积分可以近似表达为：

∫∫_Ωq_wdA≈q_wiA_i    (15)

其中，A_i指的是数值计算单元CV_i的面积，通过式(16)计算。

$q_{\mathrm{wi}}{A}_i=q_{\mathrm{wi}}^m\underset{k=1}{\overset{t}{\mathrm{\Σ}}}{\mathrm{\ψ}}_k{A}_k+\underset{l=1}{\overset{s}{\mathrm{\Σ}}}{e}_l\left|L_l\right|q_{w,l}^f---\left(16\right)$

因此饱和度方程(1)在数值计算单元CV_i上的数值计算格式为：

$A_{\mathrm{\φi}}\frac{\∂S_w^m}{\∂t}-\left[\underset{k=1}{\overset{t}{\mathrm{\Σ}}}\right|s_k\left|\right[{\mathrm{\λ}}_w^m\left(S_w^{m,\mathrm{up}}\right)\▿{\mathrm{\Φ}}_w{]}_k\·n_k+\underset{l=1}{\overset{s}{\mathrm{\Σ}}}{e}_l{\mathrm{\λ}}_w^f\left(S_m^{f,\mathrm{up}}\right)\frac{\∂{\mathrm{\Φ}}_w^f}{\∂\mathrm{\ξ}}]-q_{\mathrm{wi}}{A}_i=0---\left(17\right)$

为了避免产生无意义的解，饱和度采用上游权取值原则，具有明确的物理意义，取决于两相邻数值计算单元CV_i与CV_j之间流体流动方向，对于具有边界s_k(k＝1…t)的数值计算单元CV_i，在基岩区域利用式(18)确定，裂缝区域应用式(19)确定。

$S_w^{m,\mathrm{up}}=\left\{\begin{array}{cc}{S}_{w_i}^m& \mathrm{if}{(-\▿{\mathrm{\Φ}}_w\·n)}_k>0\\ S_{w_j}^m& \mathrm{else}.\end{array}\right.----\left(18\right)$

$S_w^{f,\mathrm{up}}=\left\{\begin{array}{cc}{S}_{w_i}^f& \mathrm{if}({\mathrm{\Φ}}_{w_i}>{\mathrm{\Φ}}_{w_j})\\ S_{w_i}^f& \mathrm{else}.\end{array}---\left(19\right)\right.$

同样地，应用上述方法得到数值计算单元CV_i的流动势方程(5)的数值计算格式为：

$\left[\underset{k=1}{\overset{t}{\mathrm{\Σ}}}\right|s_k\left|\right[{\mathrm{\λ}}^m\▿{\mathrm{\Φ}}_w+{\mathrm{\λ}}_o^m\▿{\mathrm{\Φ}}_c{]}_k\·n_k+\underset{l=1}{\overset{s}{\mathrm{\Σ}}}[{\mathrm{\λ}}^f\frac{\∂{\mathrm{\Φ}}_w}{\∂\mathrm{\ξ}}+{\mathrm{\λ}}_o^f\frac{\∂{\mathrm{\Φ}}_c}{\∂\mathrm{\ξ}}]e_l]$

$+(q_{\mathrm{wi}}+q_{\mathrm{oi}})A_i=0---\left(20\right)$

式(20)中，λ＝λ_w+λ_o表示总流度，所有流度取值同样采取如式(18)(19)所示的上游权标准来确定，此处不再赘述。裂缝中的毛管力势梯度由式(21)确定。

$\frac{\∂{\mathrm{\Φ}}_c^f}{\∂\mathrm{\ξ}}=\frac{{\mathrm{\Φ}}_{\mathrm{cj}}-{\mathrm{\Φ}}_{\mathrm{ci}}}{2\left|L_l\right|}---\left(21\right)$

(7)确定油藏的边界条件和初始条件，定压边界采用式(23)的格式，定流量边界采用式(24)的格式，初始条件为式(25)的格式。

φ_α(x，t)＝φ_α，S_α(x，t)＝S_α，onΓ_D    (23)

${\mathrm{\υ}}_{\mathrm{\α}}\·n=(-{\mathrm{\λ}}_{\mathrm{\α}}\▿{\mathrm{\Φ}}_{\mathrm{\α}})\·n=0,\▿S_{\mathrm{\α}}\·n=0,\mathrm{on}{\mathrm{\Γ}}_N---\left(24\right)$

式中，n为外边界的外法线方向向量。

φ_α(x，0)＝φ_α(x)，S_α(x，0)＝S_α(x)，att＝0    (25)

采用隐式压力显式饱和度的方法顺序求解流动势和饱和度，所有与饱和度相关变量都由前一时间步长计算得出，方程(20)用来隐式求解φ_w，控制体CV_i上计算格式如式(26)所示：

$-\left[\underset{k=1}{\overset{t}{\mathrm{\Σ}}}\right|s_k|({\mathrm{\λ}}^{m,n}{\▿\mathrm{\Φ}}_w^{n+1}{)}_k\·n_k+\underset{l=1}{\overset{s}{\mathrm{\Σ}}}{e}_l{\mathrm{\λ}}^{f,n}\frac{\∂{\mathrm{\Φ}}_w^{n+1}}{\∂\mathrm{\ξ}}]$

$=\left[\underset{k=1}{\overset{t}{\mathrm{\Σ}}}\right|s_k|{({\mathrm{\λ}}_o^{m,n}\▿{\mathrm{\Φ}}_c^n)}_k\·n_k+\underset{l=1}{\overset{s}{\mathrm{\Σ}}}{e}_l{\mathrm{\λ}}_o^{f,n}\frac{\∂{\mathrm{\Φ}}_c^n}{\∂\mathrm{\ξ}}]+(q_{\mathrm{wi}}+q_{\mathrm{oi}})A_i---\left(26\right)$

上角标n代表时间步长。

利用方程(17)显式求解含水饱和度，数值计算单元CV_i上计算格式如式(27)：

$A_{\mathrm{\φi}}\frac{S_{\mathrm{wi}}^{n+1}-S_{\mathrm{wi}}^n}{\mathrm{\Δt}}=\left[\underset{k=1}{\overset{t}{\mathrm{\Σ}}}\right|s_k\left|\right[{\mathrm{\λ}}_w^{m,n}\▿{\mathrm{\Φ}}_w^{n+1}{]}_k\·n_k+\underset{l=1}{\overset{s}{\mathrm{\Σ}}}{e}_l{\mathrm{\λ}}_w^{f,n}\frac{{\∂\mathrm{\Φ}}_w^{n+1}}{\∂\mathrm{\ξ}}]+q_{\mathrm{wi}}{A}_i---\left(27\right)$

将式(26)加以整理，写成一般形式，在整个数值计算单元系数上形成能够求解的线性方程组，求出流动势φⁿ⁺¹后，即可以利用饱和度方程(27)组成的方程组显式求解

这样可以得到每个数值计算单元节点在每个时间步长的流动势和饱和度的值，即完成了油藏的数值模拟。

下面通过几个实施例来说明：

实施例1：

研究如图3所示含大尺度裂缝的油藏模型，其尺寸大小为W(1m)×H(1m)，油藏厚度10cm，假设模型基岩为均质各向同性，基岩孔隙度为φ＝0.2，渗透率K_m＝10^-3μm²，模型左下角(0，0)m处为一口注水井，以0.01PV/d的速度注水，右上角(1，1)m处有一口生产井，以0.01PV/d的速度产液。假设模型中初始饱和油，束缚水饱和度和残余油饱和度均为0，其中水的黏度μ_w＝1mPa·s，油的黏度μ_o＝5mPa·s，水的密度ρ_w＝1000kg/m³，油的密度ρ_o＝800kg/m³。模型中存在一条方位角为θ的裂缝，其中心与模型中心点重合，其开度e＝1mm，裂缝渗透率K_f＝e²/12＝8.33×10⁴μm²。在离散裂缝模型中，裂缝由一条线表示，裂缝两端点坐标分别为(0.2，0.2)m和(0.8，0.8)m；在单孔隙介质模型中，裂缝为长度与离散裂缝模型中线段长度相等，宽度等于裂缝开度的矩形区域。基岩和裂缝的水相相对渗透率k_rw＝S_w，油相相对渗透率k_ro＝1-S_w，假设模型初始压力为20MPa，初始S_w为0。据此建立油藏的边界信息文件、裂缝信息文件及井信息文件。

分别采用本方法、常规方法I(裂缝周围网格加密的单孔隙介质模型)常规方法II(裂缝周围网格不加密的单孔隙介质模型)进行水驱油计算模拟。以常规方法I的计算结果作为参考解，验证本方法的正确性；通过比较本方法与常规方法II相当条件下的计算时间，验证本方法的高效性。

以含θ＝45°裂缝的物理模型为例，对研究区域平面进行Delaunay剖分，然后构建所需的数值计算单元。如图4为基于本方法、常规方法I和常规方法II的网格划分。图4(a)为基于本发明方法得到的划分结果，平面上用一维线单元代表裂缝，纵向裂缝高度为储层厚度，共得到878个数值计算单元节点；图4(b)为基于单常规方法I的划分结果，平面上用二维矩形区域代表裂缝，矩形长为裂缝长度，矩形宽为裂缝开度，高度为油藏厚度，其大小与研究区域长度相差3个数量级，为了计算得到较精确的结果，在裂缝周围进行了大量网格加密，共划分9343个数值计算单元节点，数量远远大于本方法；图4(c)为基于常规方法II的网格剖分结果，裂缝周围没有进行网格加密，共得到1088个数值计算单元节点，与本方法网格剖分数量相差不大。

图5列出了当注水总体积分别为0.5PV、1.5PV和3.0PV时，分别基于本方法与常规方法I的含水饱和度分布平面图，其中，图5-1(a)是当注水总体积分别为0.5PV时，采用本发明方法计算得到的实施例1油藏物理模型的含水饱和度分布平面图，图5-1(b)是当注水总体积分别为1.5PV时，采用本发明方法计算得到的含水饱和度分布平面图，图5-1(c)是当注水总体积分别为3.0PV时，采用本发明方法计算得到的含水饱和度分布平面图，图5-2(a)是当注水总体积分别为0.5PV时，采用常规方法I计算得到的含水饱和度分布平面图，图5-2(b)是当注水总体积分别为1.5PV时，采用常规方法I计算得到的的含水饱和度分布平面图，图5-2(c)是当注水总体积分别为3.0PV时，采用常规方法I计算得到的含水饱和度分布平面图。

从图中可以看出本方法能够较好地模拟大尺度裂缝的特殊导流能力，水沿裂缝快速窜进到生产井，基于常规方法I也得到同样模拟效果，二者计算结果具有一致性(从图中水沿裂缝的快速窜进，含水饱和度的值大小即可看出)。

按本方法计算后，列出了θ为0°、45°和135°方位角的单裂缝多孔介质模型注水开发过程中采收率随注水时间的变化情况，从图6中可以看到基于本方法和常规方法I得到的采收率随时间变化相同，验证了本方法的正确性。

将基于三种模型的CPU计算时间做对比，其中所用计算机CPU主频为双核2.66GHz，不同裂缝倾角情况下的多孔介质模型计算时间如表4所示：

表4

从表4中可以看出基于本方法与常规方法I精度相当，但前者计算时间远小于后者；本方法与常规方法II剖分网格数相差不多，但后者计算时间为前者的2.5-3.0倍，验证了本方法的高效性。

实施例2

图7所示的含复杂大尺度裂缝的油藏模型，油藏大小为300m×300m，厚度为10m，油藏中含有3条裂缝，裂缝起始坐标分别为(30，150)、(240，225)、(125，30)、(200，275)和(50，100)、(250，50)。均质各向同性基岩的孔隙度φ＝0.2，渗透率K_m＝1μm²，油水物性参数和裂缝参数与实施例1相同。基岩和裂缝的水相相对渗透率油相相对渗透率k_ro＝(1-S_e)²，其中S_e＝(S_w-S_wr)/(1-S_wr-S_or)，残余油饱和度为0.2，油藏初始压力为10MPa，注水井和生产井速度均为10m³/day。

据此建立油藏的边界信息文件、裂缝信息文件及井信息文件。网格剖分共得到如图8所示的1377个数值计算单元结点，并建立油藏离散文件。

图9为不同注水时刻复杂大尺度裂缝性油藏含水饱和度剖面图，其中，图9-1(a)是当注水时间为150天时，采用本发明方法计算得到的含水饱和度分布平面图，图9-1(b)是当注水时间为300天时，采用本发明方法计算得到的含水饱和度分布平面图，图9-1(c)是当注水时间为600天时，采用本发明方法计算得到的含水饱和度分布平面图，图9-2(a)是当注水时间为1200天时，采用本发明方法计算得到的含水饱和度分布平面图，图9-2(b)是当注水时间为2400天时，采用本发明方法计算得到的含水饱和度分布平面图，图9-2(c)是当注水时间为4800天时，采用本发明方法计算得到的含水饱和度分布平面图。

从图中可以看出本方法能够较好地模拟随着注水开发的进行，由于裂缝具有较大的渗透率，注入水沿裂缝大量窜进，沿着裂缝方向的水驱油速度相对于基岩大大加快，裂缝中的含水饱和度迅速增大的现象。

图10、图11分别为基于该方法模拟的复杂大尺度裂缝性油藏与均质油藏的采出程度与油井含水率随时间的变化曲线，从图中可以看出：该方法模拟结果与实际油田开发相符，裂缝的存在使得油藏采出程度变小，油井见水时间提前，且采油速度变慢，不利于油藏的开采。

从上面的实施例可以看出本方法能够很好地进行大尺度裂缝性油藏的数值模拟。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (6)

1.一种大尺度裂缝性油藏数值模拟方法，其特征在于：所述方法将油藏中的大尺度裂缝作为油藏的内边界进行降维处理，建立油藏几何模型，然后采用三角形网格对油藏几何模型进行几何剖分形成离散单元，同时形成油藏几何信息描述文件、几何离散信息文件及井信息描述文件，在离散单元基础上构建数值计算单元，建立表征此类油藏流体流动的数学模型，采用积分法在每个数值计算单元上建立其数值计算格式，形成大型线性代数方程组，并进行求解，从而实现大尺度裂缝性油藏数值模拟。

2.根据权利要求1所述的大尺度裂缝性油藏数值模拟方法，其特征在于：所述方法具体包括以下步骤：

(1)根据实际油藏特征及岩心、测井及流体的测试资料确定油藏储层物性参数，建立表征大尺度裂缝性油藏流体流动的数学模型；

(2)依据油藏实际情况及油藏描述成果，确定油藏几何边界范围，建立油藏几何模型，形成油藏几何信息描述文件；

(3)对步骤(2)建立的油藏模型进行几何离散，生成油藏模型的几何离散信息文件：平面上基岩部分采用Delaunay三角形单元进行划分，大尺度裂缝用线单元代表，纵向上采用储层厚度代表单元厚度，生成油藏模型的几何离散信息文件；

(4)基于步骤(3)的几何离散信息文件，连接所有三角形的重心与三边的中点形成联络线，以三角形的顶点为中心，其周围的所有联络线组成的多边形即为数值计算单元；

(5)计算数值计算单元节点参数值：设定饱和度变量(S_w，S_o)在每个数值计算单元内恒定不变，其数值通过测井数据取得，而数值计算单元的流动势变量(φ_w，φ_o，φ_c)由以该数值计算单元的节点为顶点的相邻三角网格的值通过线性逼近计算得到：

$\mathrm{\Φ}\left(x\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{N}_i\left(x\right){\mathrm{\Φ}}_i---\left(7\right)$

其中，n为三角形单元的顶点数；φ_i代表节点i的任意流动势变量；x为节点i的坐标；N_i(x)为形函数，为坐标的独立函数；

由式(7)可得在一个三角形内任意变量的梯度恒为：

$\▿\mathrm{\Φ}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\Φ}}_i\▿N_i\left(x\right)---\left(8\right)$

(6)对步骤(1)建立的数学模型在步骤(4)建立的数值计算单元上进行积分，建立流动控制方程的数值计算方程，数值计算单元节点参数值是由步骤(5)计算得到的；

(7)确定油藏的边界条件与初始条件，然后根据步骤(6)建立的每一个数值计算单元上的数值计算方程，构建线性代数方程组，并进行求解，完成大尺度裂缝性油藏的数值模拟。

3.根据权利要求2所述的大尺度裂缝性油藏数值模拟方法，其特征在于：所述步骤(1)中建立的表征大尺度裂缝性油藏流体流动的数学模型如下：

${\∫\∫}_{\mathrm{\Ω}}\mathrm{FEQd\Ω}={\∫\∫}_{{\mathrm{\Ω}}_m}\mathrm{FEQd}{\mathrm{\Ω}}_m+e\×{\∫}_{{\mathrm{\Ω}}_f^{\′}}\mathrm{FEQd}{\mathrm{\Ω}}_f^{\′}---\left(6\right)$

其中FEQ代表流动控制方程式(1)和式(5)：

$\frac{\∂\left(\mathrm{\φ}{S}_w\right)}{\∂t}-\▿\·({\mathrm{\λ}}_w\▿{\mathrm{\Φ}}_w)-q_w=0---\left(1\right)$

$\▿\·(({\mathrm{\λ}}_o+{\mathrm{\λ}}_w)\▿{\mathrm{\Φ}}_w)+\▿\·({\mathrm{\λ}}_o\▿{\mathrm{\Φ}}_c)+(q_o+q_w)=0---\left(5\right)$

式中：油藏整体区域用Ω表示，其中基岩部分用Ω_m代表，简化降维的大尺度裂缝部分用Ω′_f表示，e表示裂缝的开度；w，o分别表示水相和油相；φ为油藏孔隙度；S_w为油藏含水饱和度；S_o为油藏含油饱和度；为哈密顿算子；q_w，q_o为源汇项；λ_w为油藏水流度；λ_o为油藏原油流度；Φ_w为油藏水流动势；Φ_o为油藏原油流动势；φ_c为毛管压力势；p_c为毛管压力；ρ_w为油藏水的密度；ρ_o为油藏原油的密度；g为重力加速度；z为油藏垂向坐标，向上为正。

4.根据权利要求3所述的大尺度裂缝性油藏数值模拟方法，其特征在于：所述步骤(2)包括：

(21)建立外边界轨迹文件，该文件的内容包括：描述油藏边界轨迹的点的X、Y坐标值及储层深度值；

(22)将大尺度裂缝作为油藏内边界进行简化降维处理，根据油藏描述得到的大尺度裂缝具体位置，建立大尺度裂缝的具体信息文件，具体信息文件的内容包括；裂缝轨迹点对应的裂缝编号、X、Y坐标值、裂缝高度、裂缝开度；

(23)根据油藏中井的情况，建立井信息文件，该井信息文件的内容包括井号、X、Y坐标值、井所在的组别及井类别。

5.根据权利要求4所述的大尺度裂缝性油藏数值模拟方法，其特征在于：所述步骤(22)中的所述简化降维处理是指将三维油藏中的大尺度裂缝作为二维面单元，将二维油藏面上的裂缝作为一维线单元。

6.根据权利要求5所述的大尺度裂缝性油藏数值模拟方法，其特征在于：所述步骤(3)中的所述几何离散信息文件包括：基岩中三角形单元的序号、相邻三角形序号、三个顶点序号、三条边序号以及重心的坐标，裂缝中线单元的序号、两个顶点序号、中点坐标值。