一种具有泡沫油现象的稠油油藏溶解气驱数值模拟方法与流程

本发明涉及一种具有泡沫油现象的稠油油藏溶解气驱数值模拟方法，属于稠油油藏开发的模拟技术领域。

背景技术：

世界范围内稠油资源十分丰富，占油气资源总量的三分之一以上。我国探明和控制储量19亿吨以上，主要分布于辽河、新疆、胜利等十几个油田。此外，随着我国油气勘探开发的国际化，大量国外稠油资源亟待开发。因此，在常规油气勘探开发难度越来越大，能源矛盾日益突出的当下，大力发展稠油提高采收率技术符合国家经济发展战略。

与蒸汽吞吐、蒸汽重力辅助泄油等稠油热采方法相比，稠油油藏溶解气驱开发过程是经济、高效开发稠油资源的重要方法之一。在加拿大、委内瑞拉、阿尔巴尼亚以及中国部分稠油油藏溶解气驱过程中，出现了异于常规稠油油藏的生产特征，主要集中在缓慢的气油比上升速度和较高的油藏采收率上。研究表明，泡沫油现象的存在是上述特殊稠油油藏高产的主要原因。因此，如何有效模拟泡沫油现象，准确模拟该类特殊稠油溶解气驱过程，预测油田产能，制定合理开发对策成为提高该类特殊稠油油藏采收率的关键。

泡沫油现象的产生是由于地层压力降至泡点压力后，溶解气从稠油中析出，气泡开始成核，并随着压力进一步下降而生长，由于特殊稠油较高 的粘滞力和压力梯度，较低的扩散速度，气泡分散在油相中，而不是聚并形成自由气。目前，国内外学者提出了以下理论和方法来模拟上述过程：①常规稠油油藏溶解气驱模型。该模型以常规稠油油藏溶解气驱模型为基础，通过降低气相相对渗透率曲线、调高孔隙度等特性来模拟泡沫油现象。②“假泡点”模型：将假泡点压力设为可调参数，反映泡沫油藏原始压力保持水平高、生产气油比低和采收率高的异常生产特性。②分相流动模型：认为气体分相流动随气体饱和度呈直线上升，超过弥散气体积分数极限后产生游离气。因此，通过改进气相相对渗透率和组分特性来模拟泡沫油现象。③粘度降低模型，通过调低稠油粘度来模拟泡沫油现象产生的粘度降低现象。

由于泡沫油现象的生成和消失是一个动态过程，因此，上述理论和方法只通过单次调节相对渗透率曲线或物性特征无法描述气泡成核、生长等微观过程，难以准确地表征泡沫油现象对稠油溶解气驱开发效果的影响，导致模拟计算结果误差较大，因此，无法用于模拟具有泡沫油现象的稠油油藏溶解气驱过程。

技术实现要素：

针对现有理论和方法的不足，本发明提供一种具有泡沫油现象的稠油油藏溶解气驱数值模拟方法。

发明概述：

本发明所述方法从微观角度出发，建立数学模型描述气泡成核、生长规律。在此基础上，推导气体从溶解气中向气泡中的扩散速度方程，建立三维二相多组分溶解气驱数学模型，并结合有限差分方法和隐压显饱法(IMPES)对其求解，确定三维稠油油藏溶解气驱过程中油藏压力分布、 油相与气相的饱和度分布、溶解气油比、产油量和产气量。之后，以具有泡沫油现象的稠油油藏溶解气驱实验为依据，通过实验拟合的方法确定模型的不确定参数。最后，通过拟合后的模型揭示开发效果参数影响规律，理解泡沫油现象，预测油田产能，制定开发对策提高具有泡沫油现象的稠油油藏采收率。

本发明的具体技术方案如下：

一种具有泡沫油现象的稠油油藏溶解气驱数值模拟方法，包括如下步骤：

步骤1.根据油藏或者实验岩心物性参数，读入原始压力、饱和度、孔隙度和渗透率数据，初始化模型；根据模拟精度要求和计算机计算能力，划分模型网格；根据油藏或者实验岩心孔隙分布特征，选择密度分布函数；

步骤2.根据当前时间的平衡压力p_e1计算下一时间的平衡压力p_e2：油井定产液量Q生产时，根据物质平衡原理，在平衡状态下，油藏的产量等于因油藏压力下降引起的油藏孔隙体积的变化，即

$p_{e2}=p_{e1}-\frac{Q\·\mathrm{\Δ}t}{C_{t}V}---\left(I\right)$

式(I)中，p_e1和p_e2分别为当前时间与下一时间的平衡压力，MPa；初始计算时，p_e1为初始油藏压力；C_t为油藏综合压缩系数，1/MPa；V为油藏孔隙总体积，cm³；Δt为时间步长，s；

步骤3.由于气泡分散在油相中，泡沫油的油气相渗关系不仅需要考虑气相饱和度对油气相对渗透率的影响，而且需要考虑稠油的黏滞力对油气相对渗透率的影响；因此，本发明通过以下公式计算油气相对渗透率：

$k_{ro}=k_{ro}^0{\left(\frac{1-S_g}{1-S_{gc}}\right)}^{n_o}---\left(II\right)$

$k_{rg}=k_{rg}^0{\left(\frac{S_g-S_{gc}}{1-S_{gc}}\right)}^{n_g}---\left(III\right)$

其中：

$k_{rg}^0=k_{rgi}\exp (-{{\mathrm{\αv\μ}}_o}^{1/2})---\left(IV\right)$

S_gc＝S_gci+βvμ_o^1/2 (V)

式(II)-(V)中，S_g为气相饱和度；S_gc为临界气相饱和度；n_o为油相Corey指数；n_g为气相Corey指数；v为流体流动速度，cm/s；μ_o为油相粘度，mPa·s；k_ro为油相相对渗透率；k_rg为气相相对渗透率；k_rg⁰为气相相对渗透率的端点值；k_ro⁰为油相相对渗透率的端点值；k_rgi为气相相对渗透率的端点值的最大值；α和β为相关系数；S_gci为最小临界气相饱和度；

步骤4.当油藏压力p低于泡点压力以后，原油处于过饱和状态，通过以下公式判断原油的气泡是否成核：

$\mathrm{\Δ}p=p_e-p\≥\frac{2\mathrm{\γ}}{r_p}---\left(VI\right)$

$\mathrm{\&Delta;}p=p_e-p<\frac{2\mathrm{\&gamma;}}{r_p}---\left(VII\right)$

式(VI)-(VII)中，Δp为过饱和压力，MPa；p为油藏压力，MPa；p_e为平衡压力，MPa；γ为油气界面张力，dyne/cm；r_p为多孔介质孔隙中孔隙半径，m；

若则气泡开始成核，泡沫油现象产生，此时需要考虑 气泡生长的影响，因此，至步骤5；

若则气泡尚未成核，油藏中稠油仍以单相形式流动，没有产生泡沫油现象；因此，不需要计算气泡生长以及气体从溶解气中向气泡中的扩散速度，直接进行步骤7；

步骤5.由于气体浓度在气泡边界与油相中存在差异，因此，油相中的气体会不断扩散到气泡中使得气泡生长，气泡半径的计算公式如下：

$\frac{dr}{dt}=\frac{DRT(R_s-R_e){\mathrm{\&rho;}}_{gsc}}{{\mathrm{prB}}_o}---\left(VIII\right)$

式(VIII)中，为气泡半径随时间的变化率；r为气泡半径，m；D为气体在油相中的扩散系数，m²/s；R为通用气体常数；T为温度，K；R_s为泡沫油的溶解气油比，m³/m³；R_e为平衡状态下的溶解气油比，m³/m³；ρ_gsc为气体摩尔密度，mol/m³；B_o为油相体积系数，m³/m³；

步骤6.通过以下方程计算气体从溶解气中向气泡中的扩散速度q_s→b：

$q_{s\&RightArrow;b}=\frac{8{\mathrm{\&pi;RTD}}^2{(R_s-R_e)}^2{\mathrm{\&rho;}}_{gsc}}{{\mathrm{pB}}_o^2}\underset{j}{\&Sigma;}{N}_j^{4/3}{t}_j---\left(IX\right)$

式(IX)中，N_j为气泡数目，个；t_j为气泡的生长时间，s；

步骤7.根据质量守恒原则，考虑气体从溶解气中向气泡中的扩散，建立稠油、溶解气和气泡组成的三维两相数学模型如下：

$\&dtri;\&CenterDot;\&lsqb;{\mathrm{\&lambda;}}_o\&dtri;(p-{\mathrm{\&rho;}}_{o}gH)\&rsqb;+q_{vo}=\frac{\&part;}{\&part;t}\left(\frac{{\mathrm{\&phi;S}}_o}{B_o}\right)$

$\&dtri;\&CenterDot;\&lsqb;R_s{\mathrm{\&lambda;}}_o\&dtri;(p-{\mathrm{\&rho;}}_{o}gH)\&rsqb;-q_{s\&RightArrow;b}+R_s{q}_{vo}=\frac{\&part;}{\&part;t}\left(\frac{{\mathrm{\&phi;R}}_s{S}_o}{B_o}\right)$

$\&dtri;\&CenterDot;\&lsqb;{\mathrm{\&lambda;}}_g\&dtri;(p-{\mathrm{\&rho;}}_{g}gH)\&rsqb;+q_{s\&RightArrow;b}+q_{vg}=\frac{\&part;}{\&part;t}\left(\frac{{\mathrm{\&phi;S}}_g}{B_g}\right)---\left(X\right)$

式(X)中，H为标高，由某一基准面算起的垂直方向深度，m；λ_o和λ_g分别是油相和气相的流度；S_o为油相饱和度；q_vo和q_vg分别为产油量和产气量，cm³/s；φ为孔隙度；ρ_o为原油密度，kg/m³；ρ_g为天然气密度，kg/m³；B_g为气相体积系数，m³/m³；

步骤8.已知油藏的初始压力和初始饱和度，油藏外边界为封闭边界，内边界为定液生产，计算油藏的压力、油相与气相的饱和度、溶解气油比、产油量和产气量；

步骤9.根据以填砂管或者岩心为基础的稠油油藏溶解气驱实验，调整不确定参数：密度分布函数参数、油相Corey指数n_o和气相Corey指数n_g、气体在油相中的扩散系数D拟合实验数据，为后续实际油藏研究提供可靠的模型。

以拟合后的模型为基础，揭示密度分布函数参数、衰竭速率、粘度、油气界面张力、临界含气饱和度、渗透率和扩散系数对具有泡沫油现象的稠油油藏溶解气驱开发过程的影响，预测油田产能，制定开发对策提高油藏采收率。

根据本发明优选的，所述步骤1中，密度分布函数采用正态、瑞利和莱斯密度分布函数，通过实验拟合最终确定合适的密度分布函数及其参数：

正态密度分布：

瑞利密度分布：

莱斯密度分布：

式(XI)中，f(r_p)为密度分布函数，σ为标准差，μ、R_m和A均是常数， I₀()是修正的0阶第一类贝塞尔函数。

根据本发明优选的，所述步骤2中，公式(I)中综合压缩系数C_t用以下公式计算：

$C_t=S_o{\&lsqb;-\frac{1}{B_o}\frac{\&part;B_o}{\&part;p}+\frac{B_g}{B_o}\frac{\&part;R_s}{\&part;p}\&rsqb;}_T-\frac{S_g}{B_g}{\left(\frac{\&part;B_g}{\&part;p}\right)}_T+C_f---\left(XII\right)$

式(XII)中，C_f为油藏或实验岩心中岩石骨架的压缩系数，1/MPa。

根据本发明优选的，所述步骤6中，气泡数目通过以下公式计算：

$N_j=\frac{\mathrm{\&delta;}\mathrm{\&phi;}}{V_s}=\frac{\underset{r_{p1}}{\overset{r_{p2}}{\&Integral;}}f\left(r_p\right){\mathrm{dr}}_p\&CenterDot;\mathrm{\&phi;}}{V_s}---\left(XIII\right)$

式(XIII)中，δ为含有激活气泡的孔隙占总孔隙的分数；V_s为孔隙平均体积，cm³；r_p1为当前时刻t_i对应的孔隙半径，m；r_p2为下一时刻t_i+1对应的孔隙半径，m；

气泡的生长时间t_j＝t_i+1-t_i。

根据本发明优选的，所述步骤7中，联立稠油、溶解气和气泡组成的三维两相数学模型，得以下公式：

$B_g\&lsqb;\&dtri;\&CenterDot;({\mathrm{\&lambda;}}_g\&dtri;p)+{\mathrm{CG}}_g+q_{s\&RightArrow;b}+q_{vg}\&rsqb;+B_o\&lsqb;\&dtri;\&CenterDot;({\mathrm{\&lambda;}}_o\&dtri;p)+{\mathrm{CG}}_o+q_{vo}\&rsqb;={\mathrm{\&phi;C}}_t\frac{\&part;p}{\&part;t}---\left(XIV\right)$

$\&dtri;\&CenterDot;(R_s{\mathrm{\&lambda;}}_o\&dtri;p)+{\mathrm{CG}}_{sg}-q_{s\&RightArrow;b}-R_s\&lsqb;\&dtri;\&CenterDot;({\mathrm{\&lambda;}}_o\&dtri;p)+{\mathrm{CG}}_o\&rsqb;=\frac{{\mathrm{\&phi;S}}_o}{B_o}\frac{\&part;R_s}{\&part;t}---\left(XV\right)$

式(XIV)、(XV)中，CG_o＝-▽·[λ_o▽(ρ_ogH)]；CG_g＝-▽·[λ_g▽(ρ_g gH)]；

CG_sg＝-▽·[R_sλ_o▽(ρ_o gH)]。

根据本发明优选的，所述步骤8中，油藏外边界为封闭边界，通过以下公式实现：

$\frac{\&part;P}{\&part;n}{|}_{\mathrm{\&Gamma;}}=0---\left(XVI\right)$

内边界定液生产，产液量为Q，cm³/s；

计算油藏的压力、油相与气相的饱和度、溶解气油比、产油量和产气量的方法为：首先，对两个公式(XIV和XV)进行有限差分。由于公式(XIV)中只含有p一个未知数，求解公式(XIV)得到油藏压力；然后，求解公式(XV)得到泡沫油的溶解气油比R_s；利用隐压显饱法(IMPES)求解得到油相的饱和度S_o，根据归一化原理S_o+S_g＝1，得到气相的饱和度S_g，最后，利用以下公式计算产油量和产气量：

$q_{vo}=\frac{Q}{B_o(1+\frac{k_{rg}\&CenterDot;{\mathrm{\&mu;}}_o}{k_{ro}\&CenterDot;{\mathrm{\&mu;}}_g})}$

$q_{vg}=\frac{Q\left(\frac{k_{rg}\&CenterDot;{\mathrm{\&mu;}}_o}{k_{ro}\&CenterDot;{\mathrm{\&mu;}}_g}\right)}{B_g(1+\frac{k_{rg}\&CenterDot;{\mathrm{\&mu;}}_o}{k_{ro}\&CenterDot;{\mathrm{\&mu;}}_g})}---\left(XVII\right)$

式(XVI)、(XVII)中，Γ为油藏外边界，n为油藏外边界的法线方向，μ_g为气相粘度，mPa·s。

本发明的优势在于：

本发明提供了一种具有泡沫油现象的稠油油藏溶解气驱数值模拟方法，能够从微观角度出发，描述特殊稠油溶解气驱过程中气泡成核、生长规律，计算气体从溶解气中向气泡中的扩散速度，确定三维稠油油藏溶解气驱过程中的压力分布、油相与气相的饱和度分布、溶解气油比、产油量和产气量。从而能够较好地模拟泡沫油现象对稠油溶解气驱开发过程的影响，大大提高了预测的准确性，对揭示参数影响规律，理解泡沫油现象，预测油田产能，制定开发对策提高该类特殊稠油油藏采收率具有十分重要 的意义。

附图说明

图1是本发明所述模拟方法的流程图；

图2是本发明实施例1中孔隙体积分布示意图；

图3a是本发明实施例1中计算结束时填砂管中饱和度分布图，其中，稠油相饱和度分布图；

图3b是本发明实施例1中计算结束时填砂管中饱和度分布图，其中，气相饱和度分布图；

图4是本发明实施例1中实验结果与数模计算结果对比图；

图5a是本发明实施例1中油藏平均压力参数影响规律图，其中，密度分布函数中标准差σ对油藏平均压力的影响；

图5b是本发明实施例1中油藏平均压力参数影响规律图，其中，油气界面张力γ对油藏平均压力的影响；

图6a是本发明实施例2中具有泡沫油现象的稠油溶解气驱计算结束后油藏压力和油相饱和度分布示意图，其中，油藏压力分布图；

图6b是本发明实施例2中具有泡沫油现象的稠油溶解气驱计算结束后油藏压力和油相饱和度分布示意图，其中，油相饱和度分布图；

图7是本发明实施例2中具有泡沫油现象的稠油油藏溶解气驱与常规稠油油藏溶解气驱采收率计算结果对比图。

具体实施方式

下面结合实施例和说明书附图对本发明做详细的说明，但不限于此。

实施例1、

本实施例1为拟合一维具有泡沫油现象的稠油溶解气驱实验数据，确 定模型不可靠参数，为实施例2中的三维实际油藏预测提高可靠的模型参数。实验选用填砂管岩心长为60cm，横截面积为20cm²，模型渗透率为1.25D，孔隙度为0.386，岩石的压缩系数为3.2×10^-3MPa^-1。实验温度为25℃，模型初始压力为4.19MPa，原油泡点压力为3.96MPa，油气界面张力γ＝30dyne/cm，泡点压力下的原油粘度、原油体积系数、溶解气油比分别为1100cp、1.05、15.57m³/m³。填砂管左端封闭，右端定产液量Q＝8.3333×10^-4cm³/s生产，测试得到实验过程中的填砂管压力、产油量和产气量变化。

图1是本发明实施例1的计算流程图。

一种具有泡沫油现象的稠油油藏溶解气驱数值模拟方法，包括如下步骤：

步骤1.根据油藏或者实验岩心物性参数，读入原始压力、饱和度、孔隙度和渗透率数据，初始化模型；根据模拟精度要求和计算机计算能力；根据油藏或者实验岩心孔隙分布特征，选择密度分布函数；

所述步骤1中，密度分布函数采用正态、瑞利和莱斯密度分布函数，通过实验拟合最终确定合适的密度分布函数及其参数：

正态密度分布：

瑞利密度分布：

莱斯密度分布：

式(XI)中，f(r_p)为密度分布函数，σ为标准差，μ、R_m和A均是常数，I₀()是修正的0阶第一类贝塞尔函数。

本实施例1中，根据上述填砂管实际原始压力、饱和度、孔隙度和渗透率数据初始化模型；根据模拟精度要求和计算机计算能力，划分模型网格为120×1×1；根据填砂管岩心孔隙分布特征，采用瑞利密度分布函数，如图2所示，其它相关参数见表1所示。

表1模型相关参数

步骤2.根据当前时间的平衡压力p_e1计算下一时间的平衡压力p_e2：油井以产量Q＝8.3333×10^-4cm³/s生产，根据物质平衡原理，在平衡状态下，油藏的产量等于因油藏压力下降引起的油藏孔隙体积的变化，即

$p_{e2}=p_{e1}-\frac{Q\&CenterDot;\mathrm{\&Delta;}t}{C_{t}V}---\left(I\right)$

式(I)中，p_e1和p_e2分别为当前时间与下一时间的平衡压力，MPa；初始计算时，p_e1为初始油藏压力；C_t为油藏综合压缩系数，1/MPa；V为油藏孔隙总体积，cm³；Δt为时间步长，s；

所述步骤2中，公式(I)中综合压缩系数C_t用以下公式计算：

$C_t=S_o{\&lsqb;-\frac{1}{B_o}\frac{\&part;B_o}{\&part;p}+\frac{B_g}{B_o}\frac{\&part;R_s}{\&part;p}\&rsqb;}_T-\frac{S_g}{B_g}{\left(\frac{\&part;B_g}{\&part;p}\right)}_T+C_f---\left(XII\right)$

式(XII)中，C_f为油藏或实验岩心中岩石骨架的压缩系数，1/MPa。

本实施例1中，时间步长为Δt＝50s，初始油藏压力为p_e1＝4.19MPa，当前时间计算所得综合压缩系数为C_t＝3.2×10^-3MPa^-1，则下一时间t＝50s时的平衡压力p_e2＝4.16MPa；

步骤3.由于气泡分散在油相中，泡沫油的油气相渗关系不仅需要考虑气相饱和度对油气相对渗透率的影响，而且需要考虑稠油的黏滞力对油气相对渗透率的影响；因此，本发明通过以下公式计算油气相对渗透率：

$k_{ro}=k_{ro}^0{\left(\frac{1-S_g}{1-S_{gc}}\right)}^{n_o}---\left(II\right)$

$k_{rg}=k_{rg}^0{\left(\frac{S_g-S_{gc}}{1-S_{gc}}\right)}^{n_g}---\left(III\right)$

其中：

S_gc＝S_gci+βvμ_o^1/2 (V)

式(II)-(V)中，S_g为气相饱和度；S_gc为临界气相饱和度；n_o为油相Corey指数；n_g为气相Corey指数；v为流体流动速度，cm/s；μ_o为油相粘度，mPa·s；k_ro为油相相对渗透率；k_rg为气相相对渗透率；k_rg⁰为气相相对渗透率的端点值；k_ro⁰为油相相对渗透率的端点值；k_rgi为气相相对渗透率的端点值的最大值；α和β为相关系数；S_gci为最小临界气相饱和度；由于填砂管饱和油，初始油相饱和度S_o＝1，初始气相饱和度S_g＝0，油相的初始相对渗透率k_ro＝1，气相的初始相对渗透率为k_rg＝0，气泡成核之前，油相和气相的饱和 度保持不变，因此，油相和气相的相对渗透也保持不变；

步骤4.当油藏压力p低于泡点压力以后，原油处于过饱和状态，通过以下公式判断原油的气泡是否成核：

$\mathrm{\&Delta;}p=p_e-p\&GreaterEqual;\frac{2\mathrm{\&gamma;}}{r_p}---\left(VI\right)$

$\mathrm{\&Delta;}p=p_e-p<\frac{2\mathrm{\&gamma;}}{r_p}---\left(VII\right)$

式(VI)-(VII)中，Δp为过饱和压力，MPa；p为油藏压力，MPa；p_e为平衡压力，MPa；γ为油气界面张力，dyne/cm；r_p为多孔介质孔隙中孔隙半径，m；

若则气泡开始成核，泡沫油现象产生，此时需要考虑气泡生长的影响，因此，至步骤5；

若则气泡尚未成核，油藏中稠油仍以单相形式流动，没有产生泡沫油现象；因此，不需要计算气泡生长以及气体从溶解气中向气泡中的扩散速度，直接进行步骤7；

本实施例1中，油气界面张力γ＝30dyne/cm，孔隙最大半径r_pmax＝1.0e-7，该时间步，油藏压力等于平衡压力，Δp＝p_e-p＝0。因此，气泡尚未成核，油藏中稠油仍以单相形式流动，没有产生泡沫油现象。因此，不需要计算气泡生长以及气体从溶解气中向气泡中的扩散速度，即q_s→b为0，省略步骤5和步骤6；

步骤5.由于气体浓度在气泡边界与油相中存在差异，因此，油相中的气体会不断扩散到气泡中使得气泡生长，气泡半径的计算公式如下：

$\frac{dr}{dt}=\frac{DRT(R_s-R_e){\mathrm{\&rho;}}_{gsc}}{{\mathrm{prB}}_o}---\left(VIII\right)$

式(VIII)中，为气泡半径随时间的变化率；r为气泡半径，m；D为气体在油相中的扩散系数，m²/s；R为通用气体常数；T为温度，K；R_s为泡沫油的溶解气油比，m³/m³；R_e为平衡状态下的溶解气油比，m³/m³；ρ_gsc为气体摩尔密度，mol/m³；B_o为油相体积系数，m³/m³；

步骤6.通过以下方程计算气体从溶解气中向气泡中的扩散速度q_s→b：

$q_{s\&RightArrow;b}=\frac{8{\mathrm{\&pi;RTD}}^2{(R_s-R_e)}^2{\mathrm{\&rho;}}_{gsc}}{{\mathrm{pB}}_o^2}\underset{j}{\&Sigma;}{N}_j^{4/3}{t}_j---\left(IX\right)$

式(IX)中，N_j为气泡数目，个；t_j为气泡的生长时间，s；

所述步骤6中，气泡数目通过以下公式计算：

$N_j=\frac{\mathrm{\&delta;}\mathrm{\&phi;}}{V_s}=\frac{\underset{r_{p1}}{\overset{r_{p2}}{\&Integral;}}f\left(r_p\right){\mathrm{dr}}_p\&CenterDot;\mathrm{\&phi;}}{V_s}---\left(XIII\right)$

式(XIII)中，δ为含有激活气泡的孔隙占总孔隙的分数；V_s为孔隙平均体积，cm³；r_p1为当前时刻t_i对应的孔隙半径，m；r_p2为下一时刻t_i+1对应的孔隙半径，m；

气泡的生长时间t_j＝t_i+1-t_i；

步骤7.根据质量守恒原则，考虑气体从溶解气中向气泡中的扩散，建立稠油、溶解气和气泡组成的三维两相数学模型如下：

$\&dtri;\&CenterDot;\&lsqb;{\mathrm{\&lambda;}}_o\&dtri;(p-{\mathrm{\&rho;}}_{o}gH)\&rsqb;+q_{vo}=\frac{\&part;}{\&part;t}\left(\frac{{\mathrm{\&phi;S}}_o}{B_o}\right)$

$\&dtri;\&CenterDot;\&lsqb;R_s{\mathrm{\&lambda;}}_o\&dtri;(p-{\mathrm{\&rho;}}_{o}gH)\&rsqb;-q_{s\&RightArrow;b}+R_s{q}_{vo}=\frac{\&part;}{\&part;t}\left(\frac{{\mathrm{\&phi;R}}_s{S}_o}{B_o}\right)$

$\&dtri;\&CenterDot;\&lsqb;{\mathrm{\&lambda;}}_g\&dtri;(p-{\mathrm{\&rho;}}_{g}gH)\&rsqb;+q_{s\&RightArrow;b}+q_{vg}=\frac{\&part;}{\&part;t}\left(\frac{{\mathrm{\&phi;S}}_g}{B_g}\right)---\left(X\right)$

式(X)中，H为标高，由某一基准面算起的垂直方向深度，m；λ_o和λ_g分别是油相和气相的流度；S_o为油相饱和度；q_vo和q_vg分别为产油量和产气量，cm³/s；φ为孔隙度；ρ_o为原油密度，kg/m³；ρ_g为天然气密度，kg/m³；B_g为气相体积系数，m³/m³；

所述步骤7中，联立稠油、溶解气和气泡组成的三维两相数学模型，得以下公式：

$B_g\&lsqb;\&dtri;\&CenterDot;({\mathrm{\&lambda;}}_g\&dtri;p)+{\mathrm{CG}}_g+q_{s\&RightArrow;b}+q_{vg}\&rsqb;+B_o\&lsqb;\&dtri;\&CenterDot;({\mathrm{\&lambda;}}_o\&dtri;p)+{\mathrm{CG}}_o+q_{vo}\&rsqb;={\mathrm{\&phi;C}}_t\frac{\&part;p}{\&part;t}---\left(XIV\right)$

$\&dtri;\&CenterDot;(R_s{\mathrm{\&lambda;}}_o\&dtri;p)+{\mathrm{CG}}_{sg}-q_{s\&RightArrow;b}-R_s\&lsqb;\&dtri;\&CenterDot;({\mathrm{\&lambda;}}_o\&dtri;p)+{\mathrm{CG}}_o\&rsqb;=\frac{{\mathrm{\&phi;S}}_o}{B_o}\frac{\&part;R_s}{\&part;t}---\left(XV\right)$

式(XIV)、(XV)中，CG_o＝-▽·[λ_o▽(ρ_ogH)]；CG_g＝-▽·[λ_g▽(ρ_g gH)]；

CG_sg＝-▽·[R_sλ_o▽(ρ_o gH)]；

本实施例1中，由于填砂管压降实验为一维模型，因此，将上述三维两相数学模型化简为以下一维模型：

$\frac{\&part;}{\&part;x}\left({\mathrm{\&lambda;}}_o\frac{\&part;p}{\&part;x}\right)+q_{vo}=\frac{\&part;}{\&part;t}\left(\frac{{\mathrm{\&phi;S}}_o}{B_o}\right)$

$\frac{\&part;}{\&part;x}\left(R_s{\mathrm{\&lambda;}}_o\frac{\&part;p}{\&part;x}\right)-q_{s\&RightArrow;b}+R_s{q}_{vo}=\frac{\&part;}{\&part;t}\left(\frac{{\mathrm{\&phi;R}}_s{S}_o}{B_o}\right)$

$\frac{\&part;}{\&part;x}\left({\mathrm{\&lambda;}}_g\frac{\&part;p}{\&part;x}\right)+q_{s\&RightArrow;b}+q_{vg}=\frac{\&part;}{\&part;t}\left(\frac{{\mathrm{\&phi;S}}_g}{B_g}\right)$

步骤8.已知油藏的初始压力和初始饱和度，油藏外边界为封闭边界，内边界为定液生产，计算油藏的压力、油相与气相的饱和度、溶解气油比、产油量和产气量；

所述步骤8中，油藏外边界为封闭边界，通过以下公式实现：

$\frac{\&part;P}{\&part;n}{|}_{\mathrm{\&Gamma;}}=0---\left(XVI\right)$

内边界定液生产，产液量为Q，cm³/s；

计算油藏的压力、油相与气相的饱和度、溶解气油比、产油量和产气量的方法为：首先，对两个公式(XIV和XV)进行有限差分。由于公式(XIV)中只含有p一个未知数，求解公式(XIV)得到油藏压力；然后，求解公式(XV)得到泡沫油的溶解气油比R_s；利用隐压显饱法(IMPES)求解得到油相的饱和度S_o，根据归一化原理S_o+S_g＝1，得到气相的饱和度S_g，最后，利用以下公式计算产油量和产气量：

$q_{vo}=\frac{Q}{B_o(1+\frac{k_{rg}\&CenterDot;{\mathrm{\&mu;}}_o}{k_{ro}\&CenterDot;{\mathrm{\&mu;}}_g})}$

$q_{vg}=\frac{Q\left(\frac{k_{rg}\&CenterDot;{\mathrm{\&mu;}}_o}{k_{ro}\&CenterDot;{\mathrm{\&mu;}}_g}\right)}{B_g(1+\frac{k_{rg}\&CenterDot;{\mathrm{\&mu;}}_o}{k_{ro}\&CenterDot;{\mathrm{\&mu;}}_g})}---\left(XVII\right)$

式(XVI)、(XVII)中，Γ为油藏外边界，n为油藏外边界的法线方向，μ_g为气相粘度，mPa·s；

设置下一时间步，重复步骤1～8，计算得到气泡成核之前各时间步下的油藏压力、油相与气相饱和度、溶解气油比、产油量和产气量，计算时公式中q_s→b等于0。

当过饱和度时，气泡开始成核，通过步骤5和步骤6计算气泡半径以及气体从溶解气中向气泡中扩散的速度q_s→b；通过上述一维模型计算该时间步油藏压力、油相与气相饱和度、溶解气油比、产油量和产气量，公式中q_s→b不等于0；

重复步骤1～8，计算气泡成核之后各时间步下的油藏压力、油相与气相饱和度(图3a和图3b)、溶解气油比、产油量和产气量，公式中q_s→b不等于0，直到时间t等于最大时间步t_max，计算结束；

步骤9.根据以填砂管或者岩心为基础的稠油油藏溶解气驱实验，调整不确定参数：密度分布函数参数、油相Corey指数n_o和气相Corey指数n_g、气体在油相中的扩散系数D拟合实验数据，为后续实际油藏研究提供可靠的模型。

本实施例1中，调整后的密度分布函数参数、油相Corey指数n_o和气相Corey指数n_g、气体在油相中的扩散系数D见表2。由图4所示油藏平均压力的实验结果与模拟计算结果对比图可知，计算结果与实验结果基本一致，表明计算模型中调整后的孔隙分布参数、油相Corey指数n_o和气相Corey指数n_g、气体在油相中的扩散系数D有较高的可靠性，为后续实际油藏研究提高可靠的模型参数。

表2拟合的模型相关参数

以拟合后的模型为基础，揭示密度分布函数参数、衰竭速率、粘度、油气界面张力、临界含气饱和度、渗透率和扩散系数对具有泡沫油现象的 稠油油藏溶解气驱开发过程的影响，预测油田产能，制定开发对策提高油藏采收率。

本实施例1研究了密度分布函数中标准差σ和油气界面张力γ对油藏平均压力的影响。由图5a可知，标准差σ越大，最大孔隙半径r_pmax越大，气泡成核所需要的过饱和压力越小，泡沫油现象越容易形成；由图5b可知，油气界面张力γ越小，气泡成核所需要的过饱和压力也越小，气泡成核越容易。因此，油气界面张力和密度分布函数参数是影响气泡成核的重要因素，可通过降低油气界面张力促进泡沫油现象的形成，进而提高具有泡沫油现象的稠油油藏采收率。

实施例2、

如实施例1所述的一种具有泡沫油现象的稠油油藏溶解气驱数值模拟方法，其区别在于，油藏面积为300m×300m，油藏厚度为50m，油藏中部一口井以定产量15m³/d生产，其他参数采用实施例1中拟合得到的参数。划分模型网格为15×15×5，其它参数使用实施例1中实验拟合后的参数，计算步骤与实施例1相同。

模拟计算后的油藏压力分布见图6a，油相饱和度分布见图6b。由图6a可以看出，由于生产井的生产，从油藏边界到油藏中部，油藏压力逐渐减小并且减小的速度越来越快，呈现漏斗式下降。由图6b可以看出，从油藏边界到油藏中部，油相饱和度逐渐减小并且减小的速度越来越快。上述研究结果符合具有泡沫油现象的稠油油藏溶解气驱实际开发过程，表明该数值模拟方法预测结果的可靠性。

如图7可知，与常规稠油油藏溶解气驱相比，具有泡沫油现象的稠油油藏溶解气驱采收率提高3.64％，说明泡沫油现象可以有效提高稠油油藏的 采收率。由此可知，本发明一种具有泡沫油现象的稠油油藏溶解气驱数值模拟方法可以有效地模拟泡沫油现象，对于研究具有泡沫油现象的稠油油藏溶解气驱过程，预测油田产能，制定开发对策提高该类特殊油藏采收率具有十分重要的意义。