CN103246820B - 一种油气藏数值模拟计算方法 - Google Patents

Abstract

本发明提供一种油气藏数值模拟计算方法，该方法包括：求解空间离散后模拟过程各时刻的油相压力、水相压力、气相压力；根据所述空间离散后模拟过程各时刻的油相压力、水相压力、气相压力计算流体流动方向的总体积流量；根据所述流体流动方向的总体积流量计算含油饱和度、含气饱和度、含水饱和度；根据物质守恒原则，对所述含油饱和度、含气饱和度、含水饱和度进行校正。本发明采用“一步计算压力，一步计算速度，多步计算分流量和饱和度”的方法，引入速度项计算，避免了非均质性、流体相态变化、生产制度不合理等因素而引起的计算不收敛现象，不仅使数值模拟稳定性更好，而且收敛条件更易判断，适合在油气藏数值模拟常规方法陷入不收敛时应用。

Description

一种油气藏数值模拟计算方法

技术领域

本发明涉及数值模拟计算领域，具体地，涉及一种油气藏数值模拟计算方法。

背景技术

数值模拟计算在油田开发调整方案编制、开发效果预测方面得到广泛应用。常用的数值模拟计算方法主要包括隐式计算压力、显式计算饱和度（IMPES）方法，隐式交替求解（IMPIMS）方法，半隐式和全隐式方法，其中，IMPES和IMPIMS方法属于顺序求解。

当前，我国东部各大油田大多进入高含水期，聚合物驱、可动凝胶、复合驱等化学驱技术开始在矿场广泛应用，同时，由于低渗透、超低渗透油气藏注水困难，现场开展了一系列CO₂驱、CO₂吞吐等气驱现场试验，并取得理想效果。但低渗透油气藏非均质性，注化学剂或气体时导致的流体剧烈相变、各相流动能力差异巨大，模拟计算过程中生产制度不合理等因素往往使对应化学驱或气驱数值模拟计算工作量大、不收敛现象严重。

针对此问题，前人提出了“一步计算压力，多步计算分流量和饱和度”的方法，但是当各相流动能力差异巨大时，仅缩小步长仍很难使饱和度计算收敛，目前的化学驱或气驱数值模型模拟计算常需要10天以上的时间，严重影响了数值模拟历史拟合的工作效率和精度。所以，改进油气藏数值模拟的计算方法，使其在气驱或化学驱模拟计算时更加稳定，具有重要的理论和实际意义。

发明内容

本发明实施例的主要目的在于提供一种油气藏数值模拟计算方法，以解决因油藏地质模型非均质性、流体剧烈相变、各相流动能力差异大、模拟计算过程中生产制度不合理等因素导致的计算不收敛现象严重、计算稳定性差、计算时间长等问题。

为了实现上述目的，本发明实施例提供一种油气藏数值模拟计算方法，包括：

求解空间离散后模拟过程各时刻的油相压力、水相压力、气相压力；

利用以下公式，根据所述空间离散后模拟过程各时刻的油相压力、水相压力、气相压力计算流体流动方向的总体积流量：

$\left\{\begin{array}{c}{v}_o^{n+1}=-\frac{{\mathrm{kk}}_{\mathrm{ro}}^{n}A}{{\mathrm{\μ}}_o^n}[{\▿p}_o^{n+1}-({\mathrm{\γ}}_O^n+{\mathrm{\γ}}_{\mathrm{gd}}^n)\▿D]\\ v_w^{n+1}=-\frac{{\mathrm{kk}}_{\mathrm{rw}}^{n}A}{{\mathrm{\μ}}_w^n}({\▿p}_w^{n+1}-{\mathrm{\γ}}_W^n\▿D)\\ v_g^{n+1}=-\frac{{\mathrm{kk}}_{\mathrm{rg}}^{n}A}{{\mathrm{\μ}}_g^n}({\▿p}_g^{n+1}-{\mathrm{\γ}}_G^n\▿D)\\ v_t^{n+1}=v_o^{n+1}+v_w^{n+1}+v_g^{n+1}\\ {\▿p}_o^{n+1}=\frac{p_o^{n+1}-p_o^n}{\mathrm{\Δu}}\\ {\▿p}_w^{n+1}=\frac{p_w^{n+1}-p_w^n}{\mathrm{\Δu}}\\ {\▿p}_g^{n+1}=\frac{p_g^{n+1}-p_g^n}{\mathrm{\Δu}}\end{array}\right.$

其中，n为空间离散后的模拟过程各时刻，且n为自然数；分别为模拟过程第n时刻流体流动方向上的油相有效渗透率、水相有效渗透率、气相有效渗透率；分别为模拟过程第n+1时刻流体流动方向的油相体积流量、水相体积流量、气相体积流量、总体积流量；分别为模拟过程第n时刻油相粘度、水相粘度、气相粘度；分别为模拟过程第n+1时刻流体流动方向的油相压力梯度、水相压力梯度、气相压力梯度；分别为模拟过程第n时刻油组分的重度、水组分的重度、气组分的重度、溶解在油中的气组分的重度；D为流体流动方向的深度；A为流体流动方向流体通过的网格横截面积，所述网格由空间离散得到；▽为哈密顿算子；Δu为流体流动方向的空间步长；分别为模拟过程第n+1时刻、第n时刻的油相压力；分别为模拟过程第n+1时刻、第n时刻的水相压力；分别为模拟过程第n+1时刻、第n时刻的气相压力；

利用以下公式，根据所述流体流动方向的总体积流量计算含油饱和度、含气饱和度、含水饱和度：

$\left\{\begin{array}{c}{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}=\frac{{\mathrm{kk}}_{\mathrm{ro}}^{n+\frac{i}{m}}}{{\mathrm{\μ}}_o^n}\\ {\mathrm{\λ}}_w^{n+\frac{i+1}{m}}=\frac{{\mathrm{kk}}_{\mathrm{rw}}^{n+\frac{i}{m}}}{{\mathrm{\μ}}_w^n}\\ {\mathrm{\λ}}_g^{n+\frac{i+1}{m}}=\frac{{\mathrm{kk}}_{\mathrm{rg}}^{n+\frac{i}{m}}}{{\mathrm{\μ}}_g^n}\\ f_o^{n+\frac{i+1}{m}}=\frac{{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}}{{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_w^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_g^{n+\frac{i+1}{m}}}\\ f_w^{n+\frac{i+1}{m}}=\frac{{\mathrm{\λ}}_w^{n+\frac{i+1}{m}}}{{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_w^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_g^{n+\frac{i+1}{m}}}\\ f_g^{n+\frac{i+1}{m}}=\frac{{\mathrm{\λ}}_g^{n+\frac{i+1}{m}}}{{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_w^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_g^{n+\frac{i+1}{m}}}\\ v_o^{n+\frac{i+1}{m}}=f_o^{n+\frac{i+1}{m}}{v}_t^{n+1}\\ v_w^{n+\frac{i+1}{m}}=f_w^{n+\frac{i+1}{m}}{v}_t^{n+1}\\ v_g^{n+\frac{i+1}{m}}=f_g^{n+\frac{i+1}{m}}{v}_t^{n+1}\\ S_o^{n+\frac{i+1}{m}}=\frac{S_o^{n+\frac{i}{m}}{V}^{n+\frac{i+1}{m}}{\mathrm{\ρ}}_o^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_o^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}+\underset{l}{\mathrm{\Σ}}\frac{v_o^{n+\frac{i+1}{m}}{\mathrm{\ρ}}_o^{n+\frac{i+1}{m}}-v_o^{n+\frac{i}{m}}{\mathrm{\ρ}}_o^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_o^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}\\ S_w^{n+\frac{i+1}{m}}=\frac{S_w^{n+\frac{i}{m}}{V}^{n+\frac{i}{m}}{\mathrm{\ρ}}_w^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_w^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}+\underset{l}{\mathrm{\Σ}}\frac{v_w^{n+\frac{i+1}{m}}{\mathrm{\ρ}}_w^{n+\frac{i+1}{m}}-v_w^{n+\frac{i}{m}}{\mathrm{\ρ}}_w^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_w^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}\\ S_g^{n+\frac{i+1}{m}}=\frac{S_g^{n+\frac{i}{m}}{V}^{n+\frac{i}{m}}{\mathrm{\ρ}}_g^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_g^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}+\underset{l}{\mathrm{\Σ}}\frac{v_g^{n+\frac{i+1}{m}}{\mathrm{\ρ}}_g^{n+\frac{i+1}{m}}-v_g^{n+\frac{i}{m}}{\mathrm{\ρ}}_g^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_g^{n+\frac{i+1}{m}}+V^{n+\frac{i+1}{m}}}\end{array}\right.$

其中，分别为模拟过程第时刻油相流度、水相流度、气相流度；i为整数，且i＝0,1,2，...m-1；m为多步计算饱和度时，将第n+1时刻与第n时刻的时间间距Δt划分的小段数量；分别为模拟过程第时刻的油相分流量方程、水相分流量方程、气相分流量方程；分别为模拟过程第时刻流体流动方向的油相体积流量、水相体积流量、气相体积流量；分别为模拟过程第时刻流体流动方向的油相体积流量、水相体积流量、气相体积流量；分别为模拟过程第时刻的含油饱和度、含气饱和度、含水饱和度；分别为模拟过程第时刻的含油饱和度、含气饱和度、含水饱和度；为模拟过程第时刻网格单元的孔隙体积；为模拟过程第时刻网格单元的孔隙体积；分别为模拟过程第时刻的油相密度、水相密度、气相密度；分别为模拟过程第时刻的油相密度、水相密度、气相密度。

根据物质守恒原则，对所述含油饱和度、含气饱和度、含水饱和度进行校正。

借助于上述技术方案，本发明以IMPES方法为基础，改进现有的“一步计算压力，多步计算分流量和饱和度”方法，形成“一步计算压力，一步计算速度，多步计算分流量和饱和度”方法，该方法在油气藏数值模拟计算中引入速度项计算，并根据速度项计算饱和度，可避免因油藏地质模型非均质性、流体剧烈相变、各相流动能力差异大、模拟计算过程中生产制度不合理等因素导致的计算不稳定现象，本发明不仅能使数值模拟计算稳定性更好，而且收敛条件更好判断，尤其适合在油气藏化学驱或气驱数值模拟常规方法陷入不收敛时应用。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例一提供的油气藏数值模拟计算方法流程示意图；

图2是本发明实施例二提供的POLY软件的计算流程示意图；

图3是本发明实施例二提供的POLY软件与ECLIPSE软件的计算步长对比示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例一

本实施例提供一种油气藏数值模拟计算方法，如图1所示，该方法包括：

步骤S101，求解空间离散后模拟过程各时刻的油相压力、水相压力、气相压力。

该步骤分为以下四个过程：

过程1，建立线性差分方程中的压力方程，即按照现有的IMPES方法对常规线性差分方程组进行消元化简变形，消去空间离散后的含水饱和度差δS_w和含气饱和度差δS_g后，得到只包含变量为空间离散后的油相压力差δp_o的方程；

其中，δS_w＝S_w ⁿ⁺¹-S_w ⁿ；S_w ⁿ⁺¹、S_w ⁿ分别为空间离散后模拟过程第n+1时刻和第n时刻的含水饱和度；n为空间离散后的模拟过程各时刻，且n为自然数；

δS_g＝S_g ⁿ⁺¹-S_g ⁿ；S_g ⁿ⁺¹、S_g ⁿ分别为空间离散后模拟过程第n+1时刻和第n时刻的含气饱和度；

δp_o＝p_o ⁿ⁺¹-p_o ⁿ；p_o ⁿ⁺¹、p_o ⁿ分别为空间离散后模拟过程第n+1时刻、第n时刻的油相压力；

本过程中，线性差分方程由建立的油气藏流动数学模型差分得到，具体差分方法可采用现有的常规方法；n取值的上限根据实际的油气藏模拟计算过程具体决定。

过程2，利用不完全LU分解预处理共轭梯度法求解上述只包含变量为空间离散后的油相压力差δp_o的方程。

过程3，按照常规方法，根据毛管压力方程，得到空间离散后的水相压力差δp_w和空间离散后的气相压力差δp_g；

所述毛管压力方程指：p_o-p_w=p_cow、p_g-p_o=p_cog；其中，p_o为油相压力，p_w为水相压力，p_g为气相压力，p_cow为油相、水相毛管压力，p_cog为油相、气相毛管压力。空间离散后的水相压力差δp_w和空间离散后的气相压力差δp_g计算公式如下：

δp_o＝δ(p_cow+p_w)＝(p_w ⁿ⁺¹-p_cow ⁿ)-(p_w ⁿ-p_cow ⁿ)＝δp_w

δp_o=δ(p_g-p_cog)=(p_g ⁿ⁺¹-p_cog ⁿ)-(p_g ⁿ-p_cog ⁿ)=δp_g

其中，δp_w=p_w ⁿ⁺¹-p_w ⁿ；p_w ⁿ⁺¹、p_w ⁿ分别为空间离散后模拟过程第n+1时刻、第n时刻的水相压力；

δp_g＝p_g ⁿ⁺¹-p_g ⁿ；p_g ⁿ⁺¹、p_g ⁿ分别为空间离散后模拟过程第n+1时刻、第n时刻的气相压力；

过程4，根据过程2、3得到的δp_o、δp_w、δp_g计算下一时刻的油相、水相、气相压力。具体计算公式如下：p_o ⁿ⁺¹=p_o ⁿ+δp_o、p_w ⁿ⁺¹=p_w ⁿ+δp_w、p_g ⁿ⁺¹=p_g ⁿ+δp_g。

步骤S102，利用公式1，根据所述空间离散后模拟过程各时刻的油相压力、水相压力、气相压力计算流体流动方向的总体积流量；

$\left\{\begin{array}{c}{v}_o^{n+1}=-\frac{{\mathrm{kk}}_{\mathrm{ro}}^{n}A}{{\mathrm{\μ}}_o^n}[{\▿p}_o^{n+1}-({\mathrm{\γ}}_O^n+{\mathrm{\γ}}_{\mathrm{gd}}^n)\▿D]\\ v_w^{n+1}=-\frac{{\mathrm{kk}}_{\mathrm{rw}}^{n}A}{{\mathrm{\μ}}_w^n}({\▿p}_w^{n+1}-{\mathrm{\γ}}_W^n\▿D)\\ v_g^{n+1}=-\frac{{\mathrm{kk}}_{\mathrm{rg}}^{n}A}{{\mathrm{\μ}}_g^n}({\▿p}_g^{n+1}-{\mathrm{\γ}}_G^n\▿D)\\ v_t^{n+1}=v_o^{n+1}+v_w^{n+1}+v_g^{n+1}\\ {\▿p}_o^{n+1}=\frac{p_o^{n+1}-p_o^n}{\mathrm{\Δu}}\\ {\▿p}_w^{n+1}=\frac{p_w^{n+1}-p_w^n}{\mathrm{\Δu}}\\ {\▿p}_g^{n+1}=\frac{p_g^{n+1}-p_g^n}{\mathrm{\Δu}}\end{array}\right.$ （公式1）

上述公式1中，n为空间离散后模拟过程的各时刻，且n为自然数；分别为模拟过程第n时刻流体流动方向上的油相有效渗透率、水相有效渗透率、气相有效渗透率；分别为模拟过程第n+1时刻流体流动方向的油相体积流量、水相体积流量、气相体积流量、总体积流量；分别为模拟过程第n时刻油相粘度、水相粘度、气相粘度；分别为模拟过程第n+1时刻流体流动方向的油相压力梯度、水相压力梯度、气相压力梯度；分别为模拟过程第n时刻油组分的重度、水组分的重度、气组分的重度、溶解在油中的气组分的重度；D为流体流动方向的深度；A为流体流动方向流体通过的网格横截面积，所述网格由空间离散得到；▽为哈密顿算子；Δu为流体流动方向的空间步长；分别为模拟过程第n+1时刻、

第n时刻的油相压力；分别为模拟过程第n+1时刻、第n时刻的水相压力； 分别为模拟过程第n+1时刻、第n时刻的气相压力。

本实施例中，流体流动方向是二维的或三维的，对于三维问题，l∈{x,y,z}，其中l为流体流动方向，x、y、z分别为油气藏流动数学模型的三维坐标。

步骤S103，利用公式2，根据所述流体流动方向的总体积流量计算含油饱和度、含气饱和度、含水饱和度：

$\left\{\begin{array}{c}{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}=\frac{{\mathrm{kk}}_{\mathrm{ro}}^{n+\frac{i}{m}}}{{\mathrm{\μ}}_o^n}\\ {\mathrm{\λ}}_w^{n+\frac{i+1}{m}}=\frac{{\mathrm{kk}}_{\mathrm{rw}}^{n+\frac{i}{m}}}{{\mathrm{\μ}}_w^n}\\ {\mathrm{\λ}}_g^{n+\frac{i+1}{m}}=\frac{{\mathrm{kk}}_{\mathrm{rg}}^{n+\frac{i}{m}}}{{\mathrm{\μ}}_g^n}\\ f_o^{n+\frac{i+1}{m}}=\frac{{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}}{{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_w^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_g^{n+\frac{i+1}{m}}}\\ f_w^{n+\frac{i+1}{m}}=\frac{{\mathrm{\λ}}_w^{n+\frac{i+1}{m}}}{{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_w^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_g^{n+\frac{i+1}{m}}}\\ f_g^{n+\frac{i+1}{m}}=\frac{{\mathrm{\λ}}_g^{n+\frac{i+1}{m}}}{{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_w^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_g^{n+\frac{i+1}{m}}}\\ v_o^{n+\frac{i+1}{m}}=f_o^{n+\frac{i+1}{m}}{v}_t^{n+1}\\ v_w^{n+\frac{i+1}{m}}=f_w^{n+\frac{i+1}{m}}{v}_t^{n+1}\\ v_g^{n+\frac{i+1}{m}}=f_g^{n+\frac{i+1}{m}}{v}_t^{n+1}\\ S_o^{n+\frac{i+1}{m}}=\frac{S_o^{n+\frac{i}{m}}{V}^{n+\frac{i+1}{m}}{\mathrm{\ρ}}_o^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_o^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}+\underset{l}{\mathrm{\Σ}}\frac{v_o^{n+\frac{i+1}{m}}{\mathrm{\ρ}}_o^{n+\frac{i+1}{m}}-v_o^{n+\frac{i}{m}}{\mathrm{\ρ}}_o^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_o^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}\\ S_w^{n+\frac{i+1}{m}}=\frac{S_w^{n+\frac{i}{m}}{V}^{n+\frac{i}{m}}{\mathrm{\ρ}}_w^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_w^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}+\underset{l}{\mathrm{\Σ}}\frac{v_w^{n+\frac{i+1}{m}}{\mathrm{\ρ}}_w^{n+\frac{i+1}{m}}-v_w^{n+\frac{i}{m}}{\mathrm{\ρ}}_w^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_w^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}\\ S_g^{n+\frac{i+1}{m}}=\frac{S_g^{n+\frac{i}{m}}{V}^{n+\frac{i}{m}}{\mathrm{\ρ}}_g^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_g^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}+\underset{l}{\mathrm{\Σ}}\frac{v_g^{n+\frac{i+1}{m}}{\mathrm{\ρ}}_g^{n+\frac{i+1}{m}}-v_g^{n+\frac{i}{m}}{\mathrm{\ρ}}_g^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_g^{n+\frac{i+1}{m}}+V^{n+\frac{i+1}{m}}}\end{array}\right.$ （公式2）

上述公式2中，分别为模拟过程第时刻油相流度、水相流度、气相流度；i为整数，且i＝0,1,2，...m-1；m为多步计算饱和度时，将第n+1时刻与第n时刻的时间间距Δt划分的小段数量；分别为模拟过程第时刻的油相分流量方程、水相分流量方程、气相分流量方程；分别为模拟过程第时刻流体流动方向的油相体积流量、水相体积流量、气相体积流量； 分别为模拟过程第时刻流体流动方向的油相体积流量、水相体积流量、气相体积流量；分别为模拟过程第时刻的含油饱和度、含气饱和度、含水饱和度；分别为模拟过程第时刻的含油饱和度、含气饱和度、含水饱和度；为模拟过程第时刻网格单元的孔隙体积；为模拟过程第时刻网格单元的孔隙体积；分别为模拟过程第时刻的油相密度、水相密度、气相密度；分别为模拟过程第时刻的油相密度、水相密度、气相密度。

该步骤针对各相流动能力差异较大情况下，每个时间步长的饱和度变化可能非常剧烈，很难收敛的问题，采用缩小步长，将第n+1时刻与第n时刻的间距Δt分为m小段，然后再利用步骤S102中计算得到的模拟过程第n+1时刻流体流动方向的总体积流量计算模拟过程第n+1时刻的含油饱和度含气饱和度含水饱和度该步骤避免了油藏地质模型非均质性、流体相态变化、模拟计算过程中生产制度不合理等因素而引起的计算不收敛现象。

该步骤中，计算模拟过程第时刻网格单元的孔隙体积可以采用现有油气藏数值模拟计算中常用的求解网格单元孔隙体积的方法，也可以采用公式3求解：

$\left\{\begin{array}{c}{V}^{n+\frac{i+1}{m}}=V^{n+\frac{i}{m}}+{\mathrm{\δ\φ}}^{n+\frac{i+1}{m}}X\\ {\mathrm{\δ\φ}}^{n+\frac{i+1}{m}}={\mathrm{\φ}}^0{C}_R{\mathrm{\δp}}^{n+\frac{i+1}{m}}={\mathrm{\φ}}^0{C}_R(p^{n+\frac{i+1}{m}}-p^{n+\frac{i}{m}})\end{array}\right.$ （公式3）

上述公式3中，X为岩石总体积；φ⁰为参考压力下的孔隙度；C_R为孔隙压缩系数；为模拟过程第时刻的压力值；为模拟过程第时刻的压力值。

该步骤中，间距Δt的确定可以采用现有油气藏数值模拟计算中常用的时间步长计算方法，也可以采用公式4进行确定：

$\mathrm{\Δt}\≤\frac{V^n(1+S_{\mathrm{or}}^n-S_{\mathrm{wc}}^n-S_{\mathrm{gc}}^n)}{\sqrt{3}\underset{l}{\mathrm{\Σ}}{v}_t^n}$ （公式4）

上述公式4中，Vⁿ为模拟过程第n时刻油气藏总孔隙体积；分别为模拟过程第n时刻残余油饱和度、束缚水饱和度、束缚气饱和度；为模拟过程第n时刻流体流动方向的总体积流量。

步骤S104，根据物质守恒原则，对所述含油饱和度、含气饱和度、含水饱和度进行校正。

该步骤是指针对油气藏数值模拟计算过程中发生的相态变化，先假设一定的相态进行计算，再对其中物理不合理之处的计算进行校正，使之在保证物质守恒的前提下符合新旧相态交替的实际情况，然后再按新的相态计算下去。

该步骤具体为当油气藏由油水两相变为油气水三相时，按物质守恒关系式对步骤S103计算得到的含油饱和度、含气饱和度、含水饱和度参数进行校正；

校正前的饱和度计算结果为公式5：

$\left\{\begin{array}{c}{S}_o^n=1-S_w^n\\ S_g^n=0\end{array}\right.$ （公式5）

上述公式5中，分别为模拟过程第n时刻的含油饱和度、含气饱和度、含水饱和度；

校正后的饱和度计算结果为公式6：

$\left\{\begin{array}{c}{\stackrel{\‾}{p}}_b^n={\stackrel{\‾}{p}}^n=p^n+\mathrm{\Δp}\\ {\stackrel{\‾}{S}}_o^n=S_o^n+{\mathrm{\ΔS}}_o\\ {\stackrel{\‾}{S}}_w^n=S_w^n+{\mathrm{\ΔS}}_w\\ {\stackrel{\‾}{S}}_g^n={\mathrm{\ΔS}}_g=-({\mathrm{\ΔS}}_o+{\mathrm{\ΔS}}_w)\end{array}\right.$ （公式6）

上述公式6中，分别为模拟过程第n时刻经过修正后的泡点压力、油相压力、含油饱和度、含水饱和度、含气饱和度；Δp、ΔS_o、ΔS_w、ΔS_g分别为油相压力校正值、含油饱和度校正值、含水饱和度校正值、含气饱和度校正值。

该步骤中采用的物质守恒关系式为公式7：

$\left\{\begin{array}{c}\mathrm{\Δ}\left(S_o{\mathrm{\ρ}}_o\right)=0\\ \mathrm{\Δ}(S_g{\mathrm{\ρ}}_g+S_o{\mathrm{\ρ}}_{\mathrm{gd}})=0\\ \mathrm{\Δ}\left(S_w{\mathrm{\ρ}}_w\right)=0\end{array}\right.$ （公式7）

上述公式7中，Δ(S_oρ_o)为模拟过程第n+1时刻至第n时刻所有网格内油组分质量的变化；Δ(S_gρ_g+S_oρ_gd)为模拟过程第n+1时刻至第n时刻所有网格内气组分质量的变化；Δ(S_wρ_w)为模拟过程第n+1时刻至第n时刻所有网格内水组分质量的变化；ρ_gd为溶解在油相中的气组分密度。ρ_o、ρ_w、ρ_g、ρ_gd为油相、水相、气相压力的函数。

本实施例以IMPES方法为基础，改进现有的“一步计算压力，多步计算分流量和饱和度”方法，形成“一步计算压力，一步计算速度，多步计算分流量和饱和度”方法，该方法在油气藏数值模拟计算中引入速度项计算，并根据速度项计算饱和度，可避免因因油藏地质模型非均质性、流体剧烈相变、各相流动能力差异大、模拟计算过程中生产制度不合理等因素引起的计算不收敛现象，本发明不仅能使数值模拟计算稳定性更好，而且收敛条件更好判断，尤其适合在油气藏化学驱或气驱数值模拟常规方法陷入不收敛时应用。

实施例二

本实施例以一聚合物-微凝胶驱数值模拟软件（以下简称为POLY软件）的计算流程为例对实施例一提供的油气藏数值模拟计算方法的应用进行详细介绍，POLY软件基于实施例一提供的油气藏数值模拟计算方法为计算原理。

如图2所示，POLY软件的计算流程如下：

步骤S201，程序启动，建立程序运行所需的内存环境；

步骤S202，输入数值模拟所需的不同方向网格数目、各个网格的宽度、厚度、渗透率、孔隙度数据；

步骤S203，输入数值模拟所需的油气水粘度、溶解气油比、体积系数、相渗参数；

步骤S204，输入数值模拟所需每个网格的初始压力、油气水饱和度；

步骤S205，判断是否开始聚合物-微凝胶驱；若是，则执行步骤S206；若否，则执行步骤S207。

步骤S206，输入聚合物粘度、吸附、滞留、不可及孔隙体积等相关物性参数；

步骤S207，数值模拟开始计时；

步骤S208，输入数值模拟所需单井井号、射孔数据、注采数据、计算步长控制参数；

步骤S209，求解线性差分方程的压力方程，得到空间离散后的压力差δp；该步骤对应实施例一的步骤S101；

步骤S210，计算模拟过程第n+1时刻流体流动方向的总体积流量该步骤对应实施例一的步骤S102；

步骤S211，计算模拟过程第n+1时刻的含油饱和度含气饱和度含水饱和度该步骤对应实施例一的步骤S103；

步骤S212，判断计算结果是否满足计算步长控制参数；若是，则执行步骤S213；若否，则执行步骤S214。

步骤S213，根据物质守恒原则，对模拟过程第n+1时刻的含油饱和度含气饱和度含水饱和度进行校正；该步骤对应实施例一的步骤S104；

步骤S214，调整时间步长；

步骤S215，输出本时间段的计算结果；

步骤S216，判断计算时间是否满足退出要求；若是，则执行步骤S217；若否，则执行步骤S207。

步骤S217，清理内存环境，结束模拟计算。

为验证本发明提供的油气藏数值模拟计算方法的正确性和优越性，建立一个油气藏典型模型，分别使用本实施例的POLY软件和成熟商业软件ECLIPSE进行模拟计算，该油气藏典型模型基本参数如表1所列：

表1

韵律	反韵律	原油粘度/mPa.s	70
				渗透率范围/10-3μm2	0.01～3000	孔隙度范围/%	6～35
Kz/Kx	0.01	井网	反九点井网
				采液速度/PV/a	0.025	井距/m	350
油层厚度/m	20	小层数	4

为更好地测试本发明提供的油气藏数值模拟计算方法在强非均质性下的计算稳定性，典型模型渗透率分布在0.01～3000×10^-3μm²，孔隙度分布在6%～35%，最大计算步长限制为1天。如图3所示，模拟计算4年（1440天）后，POLY软件与商业软件ECLIPSE计算结果可比情况下，通过对比它们的计算步长，发现商业软件计算不收敛线性较严重，出现了长达6000次左右的不收敛，计算步长波动较大，而采用本发明算法的POLY软件计算步长稳定在0.5天。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种油气藏数值模拟计算方法，其特征在于，包括：

求解空间离散后模拟过程各时刻的油相压力、水相压力、气相压力；

利用以下公式，根据所述空间离散后模拟过程各时刻的油相压力、水相压力、气相压力计算流体流动方向的总体积流量：

$\left\{\begin{array}{c}{v}_o^{n+1}=-\frac{{\mathrm{kk}}_{ro}^{n}A}{{\mathrm{\μ}}_o^n}\[\▿p_o^{n+1}-({\mathrm{\γ}}_O^n+{\mathrm{\γ}}_{gd}^n)\▿D\]\\ v_w^{n+1}=-\frac{{\mathrm{kk}}_{rw}^{n}A}{{\mathrm{\μ}}_w^n}(\▿p_w^{n+1}-{\mathrm{\γ}}_W^n\▿D)\\ v_g^{n+1}=-\frac{{\mathrm{kk}}_{rg}^{n}A}{{\mathrm{\μ}}_g^n}(\▿p_g^{n+1}-{\mathrm{\γ}}_G^n\▿D)\\ v_t^{n+1}=v_o^{n+1}+v_w^{n+1}+v_g^{n+1}\\ \▿p_o^{n+1}=\frac{p_o^{n+1}-p_o^n}{\mathrm{\Δ}u}\\ \▿p_w^{n+1}=\frac{p_w^{n+1}-p_w^n}{\mathrm{\Δ}u}\\ \▿p_g^{n+1}=\frac{p_g^{n+1}-p_g^n}{\mathrm{\Δ}u}\end{array}\right.$

其中，n为空间离散后的模拟过程各时刻，且n为自然数；分别为模拟过程第n时刻流体流动方向上的油相有效渗透率、水相有效渗透率、气相有效渗透率；分别为模拟过程第n+1时刻流体流动方向的油相体积流量、水相体积流量、气相体积流量、总体积流量；分别为模拟过程第n时刻油相粘度、水相粘度、气相粘度；分别为模拟过程第n+1时刻流体流动方向的油相压力梯度、水相压力梯度、气相压力梯度；分别为模拟过程第n时刻油组分的重度、水组分的重度、气组分的重度、溶解在油中的气组分的重度；D为流体流动方向的深度；A为流体流动方向流体通过的网格横截面积，所述网格由空间离散得到；▽为哈密顿算子；Δu为流体流动方向的空间步长；分别为模拟过程第n+1时刻、第n时刻的油相压力；分别为模拟过程第n+1时刻、第n时刻的水相压力；分别为模拟过程第n+1时刻、第n时刻的气相压力；

利用以下公式，根据所述流体流动方向的总体积流量计算含油饱和度、含气饱和度、含水饱和度：

$\left\{\begin{array}{c}{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}=\frac{{\mathrm{kk}}_{ro}^{n+\frac{i}{m}}}{{\mathrm{\μ}}_o^n}\\ {\mathrm{\λ}}_w^{n+\frac{i+1}{m}}=\frac{{\mathrm{kk}}_{rw}^{n+\frac{i}{m}}}{{\mathrm{\μ}}_w^n}\\ {\mathrm{\λ}}_g^{n+\frac{i+1}{m}}=\frac{{\mathrm{kk}}_{rg}^{n+\frac{i}{m}}}{{\mathrm{\μ}}_g^n}\\ f_o^{n+\frac{i+1}{m}}=\frac{{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}}{{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_w^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_g^{n+\frac{i+1}{m}}}\\ f_w^{n+\frac{i+1}{m}}=\frac{{\mathrm{\λ}}_w^{n+\frac{i+1}{m}}}{{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_w^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_g^{n+\frac{i+1}{m}}}\\ f_g^{n+\frac{i+1}{m}}=\frac{{\mathrm{\λ}}_g^{n+\frac{i+1}{m}}}{{\mathrm{\λ}}_o^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_w^{n+\frac{i+1}{m}}+{\mathrm{\λ}}_g^{n+\frac{i+1}{m}}}\\ v_o^{n+\frac{i+1}{m}}=f_o^{n+\frac{i+1}{m}}{v}_t^{n+1}\\ v_w^{n+\frac{i+1}{m}}=f_w^{n+\frac{i+1}{m}}{v}_t^{n+1}\\ v_g^{n+\frac{i+1}{m}}=f_g^{n+\frac{i+1}{m}}{v}_t^{n+1}\\ S_o^{n+\frac{i+1}{m}}=\frac{S_o^{n+\frac{i}{m}}{V}^{n+\frac{i}{m}}{\mathrm{\ρ}}_o^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_o^{n+\frac{i}{m}}{V}^{n+\frac{i+1}{m}}}+\underset{l}{\Σ}\frac{v_o^{n+\frac{i+1}{m}}{\mathrm{\ρ}}_o^{n+\frac{i+1}{m}}-v_o^{n+\frac{i}{m}}{\mathrm{\ρ}}_o^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_o^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}\\ S_w^{n+\frac{i+1}{m}}=\frac{S_w^{n+\frac{i}{m}}{V}^{n+\frac{i}{m}}{\mathrm{\ρ}}_w^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_w^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}+\underset{l}{\Σ}\frac{v_w^{n+\frac{i+1}{m}}{\mathrm{\ρ}}_w^{n+\frac{i+1}{m}}-v_w^{n+\frac{i}{m}}{\mathrm{\ρ}}_w^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_w^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}\\ S_g^{n+\frac{i+1}{m}}=\frac{S_g^{n+\frac{i}{m}}{V}^{n+\frac{i}{m}}{\mathrm{\ρ}}_g^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_g^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}+\underset{l}{\Σ}\frac{v_g^{n+\frac{i+1}{m}}{\mathrm{\ρ}}_g^{n+\frac{i+1}{m}}-v_g^{n+\frac{i}{m}}{\mathrm{\ρ}}_g^{n+\frac{i}{m}}}{{\mathrm{\ρ}}_g^{n+\frac{i+1}{m}}{V}^{n+\frac{i+1}{m}}}\end{array}\right.$

其中，分别为模拟过程第时刻油相流度、水相流度、气相流度；i为整数，且i＝0,1,2，...m-1；m为多步计算饱和度时，将第n+1时刻与第n时刻的时间间距Δt划分的小段数量；分别为模拟过程第时刻的油相分流量方程、水相分流量方程、气相分流量方程；分别为模拟过程第时刻流体流动方向的油相体积流量、水相体积流量、气相体积流量；分别为模拟过程第时刻流体流动方向的油相体积流量、水相体积流量、气相体积流量；分别为模拟过程第时刻的含油饱和度、含气饱和度、含水饱和度；分别为模拟过程第时刻的含油饱和度、含气饱和度、含水饱和度；为模拟过程第时刻网格单元的孔隙体积；为模拟过程第时刻网格单元的孔隙体积；分别为模拟过程第时刻的油相密度、水相密度、气相密度；分别为模拟过程第时刻的油相密度、水相密度、气相密度；l为流体流动方向；

根据物质守恒原则，对所述含油饱和度、含气饱和度、含水饱和度进行校正。

2.根据权利要求1所述的油气藏数值模拟计算方法，其特征在于，还包括：

利用以下公式求解模拟过程第时刻网格单元的孔隙体积

$\left\{\begin{array}{c}{V}^{n+\frac{i+1}{m}}=V^{n+\frac{i}{m}}+{\mathrm{\δ\φ}}^{n+\frac{i+1}{m}}X\\ {\mathrm{\δ\φ}}^{n+\frac{i+1}{m}}={\mathrm{\φ}}^0{C}_R{\mathrm{\δp}}^{n+\frac{i+1}{m}}={\mathrm{\φ}}^0{C}_R(p^{n+\frac{i+1}{m}}-p^{n+\frac{i}{m}})\end{array}\right.$

其中，X为岩石总体积；φ⁰为参考压力下的孔隙度；C_R为孔隙压缩系数；为模拟过程第时刻的压力值；为模拟过程第时刻的压力值。

3.根据权利要求1所述的油气藏数值模拟计算方法，其特征在于，还包括：

利用以下公式，确定模拟过程第n+1时刻与第n时刻的间距Δt：

$\mathrm{\Δ}t\≤\frac{V^n(1-S_{or}^n-S_{wc}^n-S_{gc}^n)}{\sqrt{3}\underset{l}{\Σ}{v}_t^n}$

其中，Vⁿ为模拟过程第n时刻油气藏总孔隙体积；分别为模拟过程第n时刻残余油饱和度、束缚水饱和度、束缚气饱和度；为模拟过程第n时刻流体流动方向的总体积流量。