CN112163379B - 一种基于孔隙网络模型的非稳态两相相对渗透率计算方法 - Google Patents

Abstract

本发明公开了一种基于孔隙网络模型的非稳态两相相对渗透率计算方法，包括以下步骤：在孔隙网络模型的基础上，通过的非稳态渗流方程计算流体流量及压力；根据计算得到的流体流量及压力，结合相对渗透率计算公式，计算得到所述非稳态两相相对渗透率；根据驱替过程中不同含水饱和度下的两相相对渗透率数据，绘制两相相对渗透率曲线。本发明能够计算得到更符合储层物性特征的相对渗透率，并以此绘制两相相对渗透率曲线，为油气藏开发提供准确的技术指导。

Description

一种基于孔隙网络模型的非稳态两相相对渗透率计算方法

技术领域

本发明涉及油气田开发技术领域，特别涉及一种基于孔隙网络模型的非稳态两相相对渗透率计算方法。

背景技术

油气资源是现全世界使用最多与最重要的能源之一，如何高效合理开发油气资源是每一个油藏工程师面对的难题。两相(油水两相、气水两相)相对渗透率曲线是油气田开发工程中非常重要的基础数据，获取该类数据的方法常常有实验法与数值模拟法，每种方法都各有优缺点。实验方法包括常用的稳态法与非稳态法，非稳态法实验枯燥且难度大，耗时长，非稳态法(常用的JBN实验方法)虽然耗时短，但数据常常难以处理；常用的以黑油模型为基础数值模拟法往往较少考虑储层岩石的物性特征。基于孔隙网络模型的数值模拟方法一般假设岩石内部孔喉具有一定的形状，且可以通过岩心分析实验准确获取这类资料，模型复用性高，相比实验，可变条件(流速、压力等)宽泛，且稳定的孔隙网络模拟方法可无限重复使用，具有强大的经济效益。

目前，国内外许多研究人员通常采用连续介质理论研究多孔介质的多相渗流，但由于粘滞力、毛管力在孔隙尺度上具有非连续性，且考虑裂缝介质后其渗流规律如何变化尚不明确；实际流体在储层内部通常具有可压缩性(尤其是气体)，这与流体能瞬时从入口传到出口的稳态渗流理论相悖；实际生产情况中，多相流体的渗流规律极其复杂，常规的渗流理论不能够准确的指导油气藏的开发，预测其生产动态，而目前对非稳态两相渗流模拟研究方法还存在一定的局限性，常规的商业数值模拟类软件能大致模拟两相渗流过程，但是绘制的相对渗透率曲线往往与实际产生的准确性有所偏颇，这极大地限制了油气资源的开采。

发明内容

针对上述问题，本发明旨在提供一种基于孔隙网络模型的非稳态两相相对渗透率计算方法。

本发明的技术方案如下：

一种基于孔隙网络模型的非稳态两相相对渗透率计算方法，包括以下步骤：

通过的非稳态渗流方程计算流体流量及压力，所述非稳态渗流方程为：

式中：

为哈密尔顿算子；G为传导率；p为压力，MPa；φ₀为初始压力下的孔隙度，无量纲；C_t为综合压缩系数，无量纲；t为时间，s；

根据计算得到的流体流量及压力，结合相对渗透率计算公式，计算得到所述非稳态两相相对渗透率。

作为优选，若考虑汇源项，则所述非稳态渗流方程为：

式中：q为汇或源的流量，m³/s。

作为优选，若考虑网络管束，则非稳态渗流方程为：

式中：Δ'为孔隙网络模型的所有方向，当所述孔隙网络模型为SC网络模型时为6个方向，当所述孔隙网络模型为BCC网络模型时为8个方向，当所述孔隙网络模型为FCC网络模型时为12个方向；Δ”p为每个方向的压差，MPa；V_b为网格体积，cm³；Δt为时间差，s。

作为优选，计算油水两相的相对渗透率时，两相流中水相和油相的传导率分别为：

式中：g_w、g_o分别为水相传导率和油相传导率；r_ij为节点i和节点j之间的孔道半径，cm；μ_w、μ_o分别为水相粘度、油相粘度，Pa·s；l_ij为节点i和节点j之间的孔道长度，cm；B_o为油相体积系数，无量纲；

计算气水两相的相对渗透率时，两相流中水相的传导率通过式(4)进行计算，两相流中气体的传导率为：

式中：g_g为气相传导率；B_g为气相体积系数，无量纲；μ_g为气相粘度，Pa·s；p_sc为地面大气压，MPa；Z_sc为地面气体偏差因子，无量纲；T_sc为地面温度，℃；Z为地下气体偏差因子，无量纲；T表示地下温度，℃；<p>表示地下气体压力，MPa；<p>＝(p_i+p_j)/2，pi和pj为管束两端节点i和节点j的压力，MPa。

作为优选，所述相对渗透率计算公式为：

式中：k_ro、k_rg、k_rw1、k_rw2分别为油相、气相、油水两相中的水相、气水两相中的水相的相对渗透率，％；q_o(t)、q_g(t)、q_w1(t)、q_w2(t)分别为t时刻时油相、气相、油水两相中的水相、气水两相中的水相的流量，m³/s；α₁、α₂均为与岩心末端及驱替时间末端相关的经验常数；k_{ro_min}、k_{rg_min}分别为油相、气相的最小渗透率，mD；q_omax、q_gmax分别为整个流动过程中的油相、气相的最大流量，m³/s。

作为优选，本发明所述的非稳态两相相对渗透率计算方法还包括根据驱替过程中不同含水饱和度下的两相相对渗透率数据，绘制两相相对渗透率曲线的步骤。

作为优选，绘制两相相对渗透率曲线的具体步骤如下：

首先，根据所述非稳态渗流方程求得岩心末端部分两相初始状态的流量及压力，根据所述初始状态的流量及压力，结合相对渗透率计算公式，计算得到两相初始状态的相对渗透率；

然后，进行两相驱替，根据两相流动时的流动控制方程及最小时间步长，获得驱替过程中每个时间步长的含水饱和度及传导率，将所述传导率代入所述非稳态渗流方程，迭代计算获得整个驱替过程中的流量及压力变化数据，以此计算得到驱替过程中的两相相对渗透率数据；

最后，以含水饱和度为横坐标，以相对渗透率为纵坐标，绘制得到所述两相相对渗透率曲线。

作为优选，所述流动控制方程为：

μ_eff＝μ_wl_ijx_ij+μ_o(l_ij-l_ijx_ij)，μ_eff＝B_gμ_gl_ijx_ij+μ_w(l_ij-l_ijx_ij) (12)

式中：q_ij为两相的流量，cm³/s；g_ij为两相的传导率；p_i和p_j为管束两端节点i和节点j的压力，MPa；p_cij为圆管内毛管压力，MPa；μ_eff为单个管道中两相有效粘度，Pa·s；μ_w、μ_o、μ_g分别为水相粘度、油相粘度、气相粘度，Pa·s；l_ij为节点i和节点j之间的孔道长度，cm；x_ij为与凹液面位置有关的无量纲数，0≤x_ij≤1，即凹液面所在位置横坐标除以孔隙吼道的长度；B_g为气相体积系数，无量纲；γ为界面张力，mN/m；θ_ij为两相界面中润湿流体的润湿角，°；r_ij为节点i和节点j之间的孔道半径，cm。

作为优选，所述最小时间步长的计算方法如下：

对于任一与孔隙相连通的喉道i，假定其中被驱替相部分完全被驱替相流体充满所需时间为t_i，则有：

式中：V_i,nw为被驱替相流体体积，m³；q_i,w为被驱替相流体流量，m³/s；m为总迭代次数；

则驱替相流体充填所有孔隙网络的节点单元所需的最小时间步长t_min可表示为：

下一次的时间步长更新，将已经整管占据的管束将其设置为1，存在两相界面的管束设置为0，其时间t计算为：

式中：l_ij为节点i和节点j之间的孔道长度，cm；x_ij为与凹液面位置有关的无量纲数，0≤x_ij≤1，即凹液面所在位置横坐标除以孔隙吼道的长度；v为有效粘度下的线速度；

然后再统计时间步长更新这一过程中的最小时间步长，作为下一次最小时间步长。

作为优选，所述含水饱和度的计算方法如下：

V_w＝∑₁V_i (18)

式中：S_w为含水饱和度，％；V_w为孔隙在时间步长中被水相全部占据的管束体积之和，cm³；V_p为总孔隙体积，其值为每个连通管束的体积和，cm³；V_i为第i个被水相全部占据的管束体积，cm³。

与现有技术相比，本发明具有如下优点：

本发明能够通过结合两相非稳态数理方程与孔隙网络模型，计算得到更符合储层物性特征的相对渗透率，并以此绘制两相相对渗透率曲线，为油气藏开发提供准确的技术指导。

附图说明

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

图1为孔隙网络的二维局部节点示意图；

图2为岩心末端部分流量示意图；

图3为水驱油时圆管内油水界面示意图；

图4为水驱气时圆管内气水界面示意图；

图5为本发明实施例两相相对渗透率曲线结果示意图。

具体实施方式

下面结合附图和实施例对本发明进一步说明。需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的技术特征可以相互结合。除非另外定义，本发明公开使用的技术术语或者科学术语应当为本公开所属领域内具有一般技能的人士所理解的通常意义。本发明公开使用的“包括”或者“包含”等类似的词语意指出现该词前面的元件或者物件涵盖出现在该词后面列举的元件或者物件及其等同，而不排除其他元件或者物件。

一种基于孔隙网络模型的非稳态两相相对渗透率计算方法，包括以下步骤：

S1：通过的非稳态渗流方程计算流体流量及压力，所述非稳态渗流方程为：

式中：

为哈密尔顿算子；G为传导率；p为压力，MPa；φ₀为初始压力下的孔隙度，无量纲；C_t为综合压缩系数，无量纲；t为时间，s；

若考虑汇源项，则所述非稳态渗流方程为：

式中：q为汇或源的流量，m³/s。

若考虑网络管束，则利用非稳态渗流积分方程计算流体流量及压力，所述非稳态渗流积分方程为：

式中：Δ'为孔隙网络模型的所有方向，当所述孔隙网络模型为SC网络模型时为6个方向，当所述孔隙网络模型为BCC网络模型时为8个方向，当所述孔隙网络模型为FCC网络模型时为12个方向；Δ”p为每个方向的压差，MPa；V_b为网格体积，cm³；Δt为时间差，s。

公式(1)-(3)中的传导率计算方法如下：

计算油水两相的相对渗透率时，两相流中水相和油相的传导率分别为：

式中：g_w、g_o分别为水相传导率和油相传导率；r_ij为节点i和节点j之间的孔道半径，cm；μ_w、μ_o分别为水相粘度、油相粘度，Pa·s；l_ij为节点i和节点j之间的孔道长度，cm；B_o为油相体积系数，无量纲；

计算气水两相的相对渗透率时，两相流中水相的传导率通过式(4)进行计算，两相流中气体的传导率为：

式中：g_g为气相传导率；B_g为气相体积系数，无量纲；μ_g为气相粘度，Pa·s；p_sc为地面大气压，MPa；Z_sc为地面气体偏差因子，无量纲；T_sc为地面温度，℃；Z为地下气体偏差因子，无量纲；T表示地下温度，℃；<p>表示地下气体压力，MPa；<p>＝(p_i+p_j)/2，pi和pj为管束两端节点i和节点j的压力，MPa。

下面以公式(3)为例，介绍所述非稳态渗流方程的求解方法：

对公式(3)右侧时间进行差分：

式中：n为当前时刻状态下的各参数值；Δt为时间步长，s；i、j、k分别表示三维空间的三个方向；pⁿ⁺¹、pⁿ分别为下一时刻和当前时刻的压力值，MPa。

基于质量守恒定律，以图1所示的二维模型对公式(3)左边进行求解，每个节点流入流出流量体积之和为0，每个节点的流量如下：

节点1：

节点2：

节点3：

节点4：

节点5：

节点6：

公式(23)-(28)中下角标代表节点或喉道，例如p₁代表节点1处的压力，g₁₀代表节点0和节点1之间的喉道中流体的传导率。

以节点2为例，考虑隐式时间推进过程，变形得到：

以隐式方法处理所有节点得到：

将上述各节点的守恒方程写成矩阵的形式可以得到：

得到如[A]ⁿ⁺¹[P]n⁺¹＝[B]ⁿ形式的矩阵方程，其中：

采用梯度下降法求解矩阵f(P)＝AP-B：

f(P)＝AP-B (32)

f′(P)＝A (33)

其中，Pⁱ⁺¹为当前时刻压力场迭代求解过程中下一步压力场的迭代值(由初始值迭代而来)，当迭代过程中残差值f(P)误差小于1E-7(10的负7次方)时，认为收敛，此时得到的压力场即为当前时刻压力场。

通过计算所述矩阵即可得到这一时刻所有节点的流量和压力，这些流量和压力会形成流量场和压力场(矩阵形式)。

S2：根据计算得到的流体流量及压力，结合相对渗透率计算公式，计算得到所述非稳态两相相对渗透率。所述相对渗透率计算公式为：

式中：k_ro、k_rg、k_rw1、k_rw2分别为油相、气相、油水两相中的水相、气水两相中的水相的相对渗透率，％；q_o(t)、q_g(t)、q_w1(t)、q_w2(t)分别为t时刻时油相、气相、油水两相中的水相、气水两相中的水相的流量，m³/s；α₁、α₂均为与岩心末端及驱替时间末端相关的经验常数；k_{ro_min}、k_{rg_min}分别为油相、气相的最小渗透率，mD；q_omax、q_gmax分别为整个流动过程中的油相、气相的最大流量，m³/s。

上述的相对渗透率计算公式适用于各种渗透率的岩石，且能获取精确的相对渗透率。公式推导过程如下：

对于常规达西流动，两相相对渗透率可以表示为：

式中：k_rw、k_ro分别为水相相对渗透率和油相相对渗透率，％；q_w、q_s、q_o分别为两相流中水相流量、相同条件下(压差，流量，界面移动情况均相等)流过岩石内部单相流量、两相流中油相流量，cm³/s；ΔP_s、ΔP_w、ΔP_o为岩心两端稳定压差，MPa；μ_w、μ_s、μ_o分别为水相粘度、单相流体粘度、油相粘度，Pa·s。

根据式(35)，只要获取单向流动和两相流动时对应状态下各相流体的压差和流量，则可以计算相对渗透率，但在这个过程中，很难保证同模型中单相流和两相流动时的情况一样(主要是压差)，因此较难获取相同条件下的单相渗透率，故该公式计算误差较大。

在中低渗岩石中，如图2所示，任意时刻q_in不等于q_out，不满足达西渗流无法直接用进出口流量计算相渗，取模型或岩石的末端部分，即近似稳态渗流状态下的流动区域内，得到q_out＝q_out'，计算该区域内含水饱和度，以及该区域内的压力、粘度、分流量等数据进行相渗计算。

以水驱油为例，在驱替过程中，岩心末端部分流体流动状态分为三个阶段：第一阶段，单相油流；第二阶段，油水两相流；第三阶段，单相水流，整个计算过程忽略毛管力影响，饱和度计量也只取所述的岩心末端部分。

在第一阶段与第三阶段，q_s等于不同时刻t时的油相流量q_o与水相流量q_w之和(此时管束内流体视为单相流)，油相相对渗透率计算公式为：

此时：

在两相流阶段，可以近似认为油相流量与水相流量之和占总流量的100％，油相流量占比40％～60％，水相流量占比40％～60％，此时单管内两相界面中水相入口处至水驱前缘总长l_w占管束总长度40％～60％，l_o也占管束总长度40％～60％，由泊肃叶定律，两相流动时，每个存在两相界面的毛管中，单一相流量为：

因此：

所以油相相对渗透率为：

又因为：

所以：

同理可得：

考虑到所取位置越靠近岩心末端，流动越接近达西流，令：

得到：

油驱水或水驱气或气驱水等理论相同，综上得到式(8)-式(10)所示的能够适用于各种储层且计算结果精确的相对渗透率计算公式。

在另一个具体的实施例中，所述非稳态两相相对渗透率计算方法还包括步骤S3：根据驱替过程中不同含水饱和度下的两相相对渗透率数据，绘制两相相对渗透率曲线，具体步骤如下：

S31：根据所述非稳态渗流方程求得两相初始状态的流量及压力，根据所述初始状态的流量及压力，结合相对渗透率计算公式，计算得到两相初始状态的相对渗透率。

S32：进行两相驱替，根据两相流动时的流动控制方程及最小时间步长，获得驱替过程中每个时间步长(所述时间步长大于等于所述最小时间步长)的含水饱和度及传导率，将所述传导率代入所述非稳态渗流方程，迭代计算获得整个驱替过程中的流量及压力变化数据，以此计算得到驱替过程中的两相相对渗透率数据。

两相驱替过程中，两相界面如图3和图4所示，所述流动控制方程为：

μ_eff＝μ_wl_ijx_ij+μ_o(l_ij-l_ijx_ij)，μ_eff＝B_gμ_gl_ijx_ij+μ_w(l_ij-l_ijx_ij) (12)

式中：q_ij为两相的流量，cm³/s；gij为两相的传导率；p_i和p_j为管束两端节点i和节点j的压力，MPa；p_cij为圆管内毛管压力，MPa；μ_eff为单个管道中两相有效粘度，Pa·s；μ_w、μ_o、μ_g分别为水相粘度、油相粘度、气相粘度，Pa·s；l_ij为节点i和节点j之间的孔道长度，cm；x_ij为与凹液面位置有关的无量纲数，0≤x_ij≤1，即凹液面所在位置横坐标除以孔隙吼道的长度；B_g为气相体积系数，无量纲；γ为界面张力，mN/m；θ_ij为两相界面中润湿流体的润湿角，°；r_ij为节点i和节点j之间的孔道半径，cm。

所述最小时间步长的计算方法如下：

对于任一与孔隙相连通的喉道i，假定其中被驱替相部分完全被驱替相流体充满所需时间为t_i，则有：

式中：V_i,nw为被驱替相流体体积，m³；q_i,w为被驱替相流体流量，m³/s；m为总迭代次数；

则驱替相流体充填所有孔隙网络的节点单元所需的最小时间步长t_min可表示为：

下一次的时间步长更新，将已经整管占据的管束将其设置为1，存在两相界面的管束设置为0，其时间t计算为：

式中：l_ij为节点i和节点j之间的孔道长度，cm；x_ij为与凹液面位置有关的无量纲数，0≤x_ij≤1，即凹液面所在位置横坐标除以孔隙吼道的长度；v为有效粘度下的线速度；

然后再统计时间步长更新这一过程中的最小时间步长，作为下一次最小时间步长。

所述含水饱和度的计算方法如下：

V_w＝∑₁V_i (18)

式中：S_w为含水饱和度，％；V_w为孔隙在时间步长中被水相全部占据的管束体积之和，cm³；V_p为总孔隙体积，其值为每个连通管束的体积和，cm³；V_i为第i个被水相全部占据的管束体积，cm³。

S33：以含水饱和度为横坐标，以相对渗透率为纵坐标，绘制得到所述两相相对渗透率曲线。

在一个具体的实施例中，进行气驱水、水驱油、水驱气三种驱替实验，采用本发明所述的非稳态两相相对渗透率计算方法计算驱替过程中的相对渗透率，并绘制两相相对渗透率曲线，结果如图5所示。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容做出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (9)

1.一种基于孔隙网络模型的非稳态两相相对渗透率计算方法，其特征在于，包括以下步骤：

通过的非稳态渗流方程计算流体流量及压力，所述非稳态渗流方程为：

式中：

为哈密尔顿算子；G为传导率；p为压力，MPa；φ₀为初始压力下的孔隙度，无量纲；C_t为综合压缩系数，无量纲；t为时间，s；

根据计算得到的流体流量及压力，结合相对渗透率计算公式，计算得到所述非稳态两相相对渗透率；所述相对渗透率计算公式为：

式中：k_ro、k_rg、k_rw1、k_rw2分别为油相、气相、油水两相中的水相、气水两相中的水相的相对渗透率，％；q_o(t)、q_g(t)、q_w1(t)、q_w2(t)分别为t时刻时油相、气相、油水两相中的水相、气水两相中的水相的流量，m³/s；α₁、α₂均为与岩心末端及驱替时间末端相关的经验常数；k_{ro_min}、k_{rg_min}分别为油相、气相的最小渗透率，mD；q_omax、q_gmax分别为整个流动过程中的油相、气相的最大流量，m³/s。

2.根据权利要求1所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法，其特征在于，若考虑汇源项，则所述非稳态渗流方程为：

式中：q为汇或源的流量，m³/s。

3.根据权利要求2所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法，其特征在于，若考虑网络管束，则所述非稳态渗流方程为：

式中：Δ'为孔隙网络模型的所有方向，当所述孔隙网络模型为SC网络模型时为6个方向，当所述孔隙网络模型为BCC网络模型时为8个方向，当所述孔隙网络模型为FCC网络模型时为12个方向；Δ”p为每个方向的压差，MPa；V_b为网格体积，cm³；Δt为时间差，s。

4.根据权利要求1-3中任意一项所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法，其特征在于，计算油水两相的相对渗透率时，两相流中水相和油相的传导率分别为：

式中：g_w、g_o分别为水相传导率和油相传导率；r_ij为节点i和节点j之间的孔道半径，cm；μ_w、μ_o分别为水相粘度、油相粘度，Pa·s；l_ij为节点i和节点j之间的孔道长度，cm；B_o为油相体积系数，无量纲；

计算气水两相的相对渗透率时，两相流中水相的传导率通过式(4)进行计算，两相流中气相的传导率为：

式中：g_g为气相传导率；B_g为气相体积系数，无量纲；μ_g为气相粘度，Pa·s；p_sc为地面大气压，MPa；Z_sc为地面气体偏差因子，无量纲；T_sc为地面温度，℃；Z为地下气体偏差因子，无量纲；T表示地下温度，℃；<p>表示地下气体压力，MPa；<p>＝(p_i+p_j)/2，p_i和p_j为管束两端节点i和节点j的压力，MPa。

5.根据权利要求1-3中任意一项所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法，其特征在于，还包括根据驱替过程中不同含水饱和度下的两相相对渗透率数据，绘制两相相对渗透率曲线的步骤。

6.根据权利要求5所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法，其特征在于，绘制两相相对渗透率曲线的具体步骤如下：

选取岩心末端部分，根据所述非稳态渗流方程求得该区域内两相初始状态的流量及压力，根据所述初始状态的流量及压力，结合相对渗透率计算公式，计算得到两相初始状态的相对渗透率；

然后，进行两相驱替，根据两相流动时的流动控制方程及最小时间步长，获得驱替过程中每个时间步长的含水饱和度及传导率，将所述传导率代入所述非稳态渗流方程，迭代计算获得整个驱替过程中的流量及压力变化数据，以此计算得到驱替过程中的两相相对渗透率数据；

最后，以含水饱和度为横坐标，以相对渗透率为纵坐标，绘制得到所述两相相对渗透率曲线。

7.根据权利要求6所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法，其特征在于，所述流动控制方程为：

μ_eff＝μ_wl_ijx_ij+μ_o(l_ij-l_ijx_ij)，μ_eff＝B_gμ_gl_ijx_ij+μ_w(l_ij-l_ijx_ij) (12)

式中：q_ij为两相的流量，cm³/s；g_ij为两相的传导率；p_i和p_j为管束两端节点i和节点j的压力，MPa；p_cij为圆管内毛管压力，MPa；μ_eff为单个管道中两相有效粘度，Pa·s；μ_w、μ_o、μ_g分别为水相粘度、油相粘度、气相粘度，Pa·s；l_ij为节点i和节点j之间的孔道长度，cm；x_ij为与凹液面位置有关的无量纲数，0≤x_ij≤1，即凹液面所在位置横坐标除以孔隙吼道的长度；B_g为气相体积系数，无量纲；γ为界面张力，mN/m；θ_ij为两相界面中润湿流体的润湿角，°；r_ij为节点i和节点j之间的孔道半径，cm。

8.根据权利要求6所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法，其特征在于，所述最小时间步长的计算方法如下：

对于任一与孔隙相连通的喉道i，假定其中被驱替相部分完全被驱替相流体充满所需时间为t_i，则有：

式中：V_i,nw为被驱替相流体体积，m³；q_i,w为被驱替相流体流量，m³/s；m为总迭代次数；

则驱替相流体充填所有孔隙网络的节点单元所需的最小时间步长t_min可表示为：

下一次的时间步长更新，将已经整管占据的管束将其设置为1，存在两相界面的管束设置为0，其时间t计算为：

式中：l_ij为节点i和节点j之间的孔道长度，cm；x_ij为与凹液面位置有关的无量纲数，0≤x_ij≤1，即凹液面所在位置横坐标除以孔隙吼道的长度；v为有效粘度下的线速度；

然后再统计时间步长更新这一过程中的最小时间步长，作为下一次最小时间步长。

9.根据权利要求6所述的基于孔隙网络模型的非稳态两相相对渗透率计算方法，其特征在于，所述含水饱和度的计算方法如下：

V_w＝∑₁V_i (18)

式中：S_w为含水饱和度，％；V_w为孔隙在时间步长中被水相全部占据的管束体积之和，cm³；V_p为总孔隙体积，其值为每个连通管束的体积和，cm³；V_i为第i个被水相全部占据的管束体积，cm³。

Images