CN112630090B - 一种考虑稠油体积膨胀的甲烷-二氧化碳-丙烷混合气扩散系数的计算方法 - Google Patents

Abstract

一种考虑稠油体积膨胀的甲烷‑二氧化碳‑丙烷混合气扩散系数的计算方法，本发明以菲克定律和质量守恒原理为基础，考虑了扩散过程中稠油体积的膨胀与气体压缩因子的变化，结合本发明提出的二元相互作用系数计算公式计算了稠油体系密度，并利用隐式差分法建立了一种新的混合气扩散系数计算方法，该方法能准确计算扩散实验过程中气体扩散系数、体系密度和气相浓度等参数。解决了长期以来对混合气在稠油中扩散系数研究不足的问题，并克服了现有模型未同时考虑稠油体积的膨胀与气体压缩因子的变化的缺陷。满足了生产现场工作人员简便掌握稠油注甲烷‑二氧化碳‑丙烷混合气体扩散研究关键参数的需求。本发明对于研究稠油注气提高采收率技术、油气传质理论及油田地面工程建设等具有重要意义，因此能够在稠油注气高效开发领域广泛应用。

Description

一种考虑稠油体积膨胀的甲烷-二氧化碳-丙烷混合气扩散系 数的计算方法

技术领域

本发明涉及一种考虑稠油体积膨胀的甲烷-二氧化碳-丙烷混合气扩散系数的计算方法，属于稠油油藏开发的技术领域。

背景技术

在新的能源战略形势下，稠油油藏提高采收率技术亟待发展。目前，稠油常用的开采方式是热力采油，然而对于薄层、低渗、含边底水等复杂油气藏，稠油热采技术存在热损大、效率低、伤害储层等问题，并且开采成本较高，不够环保。但是稠油注气技术可有效避免上述问题。稠油注气技术能促进稠油膨胀，补充地层能量，降低稠油密度和粘度，减小稠油流动阻力，具有良好的应用价值与经济效益。

在稠油注气过程中，油藏流体会产生浓度梯度，气体分子会从高浓度区域向低浓度区域迁移，而扩散系数就是定量描述气体在体系中扩散能力的重要参数。稠油注气技术的效率很大程度上取决于气体在稠油中的扩散能力以及气体溶解后原油物性的变化。因此，建立气体在稠油中的扩散模型对于稠油注气提高采收率研究，稠油油藏经济高效开发具有重要意义。稠油注气技术常用气体包括甲烷、二氧化碳与丙烷等，其中丙烷可以提高稠油的膨胀能力，降低稠油的粘度和界面张力；甲烷与二氧化碳为稠油油藏开采时产出气的主要成分，价格低廉，加入后可有效防止丙烷液化，从而减少丙烷使用量。因此确定上述混合气体的扩散系数、稠油-混合气体系密度和气相浓度对于稠油注混合气提高采收率技术具有重要意义。

近年来稠油注气开发技术得以迅速发展，国内外相关研究人员对气体在稠油中扩散规律进行了一些有针对性的研究，取得了许多成果，但仍然存在部分尚待解决的问题，总结如下：

①气体在稠油中扩散时，稠油体积逐渐膨胀，气相压力逐渐减小，气体压缩因子发生变化。但目前所建立的气体扩散计算方法往往没有同时考虑稠油体积膨胀和气体压缩因子的变化。

②目前在稠油注气扩散实验中大多数还是仅使用单一气体，对甲烷-二氧化碳-丙烷混合气体在稠油中扩散行为研究尚未涉及。

③稠油-甲烷-二氧化碳-丙烷体系属于非常规流体，存在超额体积，因此现有混合物密度计算公式对该体系不适用，难以精确描述气体扩散过程中的膨胀效应。

发明内容

针对现有技术的不足，本发明提供一种考虑稠油体积膨胀的甲烷-二氧化碳-丙烷混合气扩散系数的计算方法。

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

一种考虑稠油体积膨胀的甲烷-二氧化碳-丙烷混合气扩散系数的计算方法，其特征在于，计算方法包括如下步骤：

步骤1.确定油藏温度T与混合气的初始注入压力P(t₀)、混合气的初始气相高度h_g(t₀)、扩散实验容器高度H、油气界面面积S、稠油质量m_o；

各组分相对分子质量：M₁是丙烷的相对分子质量、M₂是二氧化碳的相对分子质量、M₃是甲烷的相对分子质量、M₄是稠油的相对分子质量；

混合气体的临界压缩因子Z_c；

混合气中各组分所占摩尔分数：y₁是丙烷所占摩尔分数、y₂是二氧化碳所占摩尔分数、y₃是甲烷所占摩尔分数；

气体组分性质参数：T_c1是丙烷的临界温度、T_c2是二氧化碳的临界温度、T_c3是甲烷的临界温度；P_c1是丙烷的临界压力、P_c2是二氧化碳的临界压力、P_c3是甲烷的临界压力；w₁是丙烷的偏心因子、w₂是二氧化碳的偏心因子、w₃是甲烷的偏心因子；

步骤2.进行稠油注甲烷-二氧化碳-丙烷混合气扩散实验：扩散容器内气体扩散示意图如图1所示，获取混合气压力P随实验时间t变化数据；

步骤3.从实验总时长出发选择一个时间步长Δt，确定多个时刻t∈(t₀＝0；t₁＝t₀+Δt；......t_x＝t₀+xΔt)，根据本发明所提出的气相压力公式(1)计算气体扩散过程中各时刻下的气相压力P(t)：

P(t)＝B₁exp(-B₂t)+B₃exp(-B₄t)+B₅exp(-B₆t)+B₇ (1)

式中，B₁，B₂，B₃，B₄，B₅，B₆，B₇为公式系数，通过拟合步骤2的实验数据得到；t为实验时间，h；P(t)为时刻t时的气相压力，MPa。

步骤4.根据气体压缩因子计算模型，即公式(2)～(4)计算在各时刻下的混合气体压缩因子Z(t)：

式中，Z(t)为各时刻下混合气压缩因子，无量纲；ρ_pr(t)为各时刻下的气体对应密度，无量纲；T为油藏温度，K；n为气体种数；P_cj为各气体组分的临界压力，kPa；T_cj为各气体组分的临界温度，K；Z_c为混合气临界压缩因子；A₁～A₁₁为公式常数，无量纲，其具体值在是实施例中给出。

步骤5.根据以下脱气稠油密度公式，即公式(5)～(8)计算各时刻下的脱气稠油密度ρ₄(t)：

ρ_s＝a₁+a₂T+a₃T² (6)

A＝c₁+c₂T (7)

B＝b₁+b₂T^-1+b₃T^-2 (8)

式中，ρ₄(t)为各时刻下的脱气稠油密度，kg/m³；T为油藏温度，K；ρ_s为大气压下脱气稠油密度，kg/m³；A、B、a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂均为公式常数，无量纲，其具体数值在实施例给出。

步骤6.根据变形的PR状态方程，即公式(9)计算各时刻下的各气体参数C_j(t)值，进而得到各气体摩尔体积V_j(t)，再通过公式(18)计算在各时刻下稠油-甲烷-二氧化碳-丙烷体系中各气体组分密度ρ_j(t)：

a_j＝a_cjα(T_rj,ω_j) (j＝1,2,3) (13)

T_rj＝T/T_cj (j＝1,2,3) (16)

式中，A_j(t)、B_j(t)和C_j(t)为PR状态方程变形时所引入的参量；a_j为吸引力系数，无量纲；b_j为范德华摩尔体积，m³/kmol；V_j(t)为时刻t下各气体摩尔体积，m³/kmol；R为通用气体常数，8.314kPa·m³/(K·kmol)；T_rj为混合气中各气体组分的折算温度，无量纲；ω_j为各气体组分的偏心因子，无量纲；M_j为各气体组分的相对分子质量，g/mol；ρ_j(t)为时刻t下各气体组分的密度，kg/m³。

步骤7.根据本发明拟合所得的公式(20)计算不同时刻下各气体组分间的二元相互作用系数，其中β₁₂为C₃H₈和CO₂二元相互作用系数、β₁₃为C₃H₈和CH₄二元相互作用系数、β₂₃为CO₂和CH₄二元相互作用系数；根据本发明拟合所得的公式(21)计算不同时刻下气体与稠油间的密度二元相互作用系数β_ij(t)，其中β₁₄为C₃H₈和稠油二元相互作用系数、β₂₄为CO₂和稠油二元相互作用系数、β₃₄为CH₄和稠油二元相互作用系数：

式中，ρ_i(t)为时刻t下组分i的密度，kg/m³；ρ_j(t)为时刻t下组分j的密度，kg/m³；v_Nij(t)为时刻t下组分i和组分j的归一化比容差，无量纲；式(20)中的β_ij(t)为气体组分间的二元相互作用系数，无量纲；式(21)中的β_ij(t)为气体与稠油间的二元相互作用系数，无量纲。

步骤8.估计一个t₁时刻下气相高度h_g(t₁)的值，结合质量守恒原理与真实气体状态方程，按公式(23)计算t₁时刻下的溶解气质量m_d(t₁)：

M＝y₁M₁+y₂M₂+y₃M₃ (22)

式中，m_d(t₁)为t₁时刻下的气体扩散质量，kg；M₁、M₂、M₃分别为甲烷-二氧化碳-丙烷混合气体的相对分子质量，kg/kmol；P(t₀)为初始气相压力，kPa；Z(t₀)为初始状态下气体压缩因子；h_g(t₀)为初始气相高度，m；P(t₁)为t₁时刻下的气体压缩因子；Z(t₁)为t₁时刻下的气体压缩因子；h_g(t₁)为t₁时刻下的气相高度，m；S为油气界面的面积，m²。

步骤9.根据超额体积混合准则，按照公式(30)计算在t₁时刻下的溶气稠油的密度ρ_mix(t₁)：

x₁(t₁)＝y₁n_g(t₁) (25)

x₂(t₁)＝y₂n_g(t₁) (26)

x₃(t₁)＝y₃n_g(t₁) (27)

x₄(t₁)＝1-n_g(t₁) (28)

式中，n_g(t₁)为t₁时刻下溶气稠油中混合气的物质的量，mol；m_d(t₁)为t₁时刻下的气体扩散质量，kg；x_i(t₁)为t₁时刻下溶气稠油中组分i的物质的量，mol；n为混合体系中组分总数；ρ_mix(t₁)时刻下为t₁时刻下溶气稠油密度，kg/m³；w_i(t₁)为t₁时刻下组分i的质量分数，无量纲；ρ_i(t₁)为t₁时刻下组分i的密度，kg/m³；w_j(t₁)为t₁时刻下为组分j的质量分数，无量纲；ρ_j(t₁)为t₁时刻下组分j的密度，kg/m³；β_ij(t₁)为t₁时刻下组分i与组分j的密度二元相互作用系数。

步骤10.计算t₁时刻下溶气稠油体积V(t₁)，进而求得t₁时刻下的气相高度h_g’(t₁)，比较步骤8的h_g(t₁)估算值与步骤10的h_g’(t₁)计算值：

若相对误差小于等于0.5×10^-2，则计算结果满足精度要求；

若相对误差大于0.5×10^-2，以步骤10所得h_g’(t₁)为新估算值，重复步骤8～10，直到相对误差满足精度要求，此时再利用h_g’(t₁)计算值按式(24)计算t₁时刻下的溶解气质量m_d(t₁)；

h_g ^,(t₁)＝H-h_t(t₁) (33)

式中，V₁(t₁)为t₁时刻下溶气稠油体积，m³；m_o为稠油质量，kg；h(t₁)为t₁时刻下的油相高度，m；h_g’(t₁)为t₁时刻下气相高度，m；H为扩散实验容器高度，m。

步骤11.代入不同时刻下的混合气压力P(t)与实验时间t数据，重复步骤8～10所述计算方法，利用式(24)求出各时刻下的溶解气质量m_d(t_i)、利用式(30)求取溶气稠油密度ρ(t_i)，利用式(33)求取气相高度h_g’(t_i)(i＝2,3,…,n)；

步骤12.根据亨利定律，计算各时刻下油气界面的气体饱和浓度c_sat(t)：

式中，P(t)为t时刻下的气相压力，kPa；K_j为第j种气体的亨利常数，kPa·m³/kg；c_j(t)为t时刻下第j种气体在油气界面的浓度，kg/m³；y_j为第j种气体的摩尔分数，mol％；n为混合气体中气体的种类。

步骤13.如公式(36)与公式(37)所示，通过菲克定律与质量守恒原理来求解各时刻下的气体扩散系数：

进一步优选的，扩散实验开始时稠油中不含气相，因此求解初始条件如公式(38)所示；扩散过程中油气界面饱和气体，两相平衡，且气体扩散至容器底部时，无法继续向下扩散，因此上下边界条件分别如公式(39)与(40)所示；

进一步优选的，为便于计算，对公式(36)与(37)进行离散化处理，因此根据扩散实验容器高度，将稠油油藏划分为N个网格，网格高度按公式(44)计算。根据公式(41)～(43)，结合边界条件使用隐式差分法求解各时刻下的混合气扩散系数D^j与各网格中的气相浓度c^j _i，最后根据公式(45)求取稠油中的平均气相浓度c^j _p；

c|_t＝0＝0(0≤x＜h_t) (38)

式中，c为气体浓度，kg/m³；t为实验时间，s；D为扩散系数，m²/s，x为坐标位置，m；D^j为第j时刻下的气体扩散系数，m²/s；λ为差分过程所引入的变量，无物理意义；Δt为时间步长，s；Δx为空间步长，m；c^j+1 _i-1为第j+1时刻下第i-1网格中的气相浓度，kg/m³；c^j+1 _i为第j+1时刻下第i网格中的气相浓度，kg/m³；c^j+1 _i+1为第j+1时刻下第i+1网格中的气相浓度，kg/m³；c^j _i为第j时刻下第i网格中的气相浓度，kg/m³；c^j _sat为第j时刻下的油气界面气相浓度，kg/m³；c^j ₁为第j时刻下第1个网格中的气相浓度，kg/m³；h_g ^j+1为第j+1时刻下混合气相高度，m；P^{j +1}为第j+1时刻下混合气压力，m；Z^j+1为第j+1时刻下混合气压缩因子，无量纲；h_g ^j为第j时刻下混合气相高度，m；P^j为第j时刻下混合气压力，m；Z^j为第j时刻下混合气压缩因子，无量纲；c^j _p为第j时刻下的平均气相浓度，kg/m³；N为网格数；h_t ^j为第j时刻下油相高度，m。

本发明的有益效果

1.目前稠油注气扩散系数研究中大多数研究人员还是使用单一气体，本发明所提出的扩散系数计算模型能够准确计算扩散实验中各时刻下甲烷-二氧化碳-丙烷混合气的扩散系数、稠油-混合气体系密度和气相浓度，对于稠油注混合气提高采收率技术具有重要意义。

2.本发明所述的计算方法解决了现有模型未同时考虑混合气扩散过程中气体压缩因子变化和稠油体积膨胀效应的缺陷，提高了模型计算的准确性。

3.本发明所采用的计算方法简便，所需基础物性参数少，便于现场技术人员在稠油油藏注混合气开发初期快速获取扩散系数、体系密度和气相浓度等关键参数的需求。

4.本发明所述计算方法得到的稠油注甲烷-二氧化碳-丙烷混合气扩散实验的关键参数对于研究稠油注气提高采收率技术、油气传质理论，及油田地面工程建设等具有重要意义，因此能够在稠油注气高效开发领域广泛应用。

附图说明

图1是本发明所述稠油注甲烷-二氧化碳-丙烷混合气扩散实验中扩散容器内气体扩散示意图。

图2是本发明提出的考虑稠油体积膨胀的甲烷-二氧化碳-丙烷混合气扩散系数计算流程图。

图3是本发明实施实例1中气相压力实验值与拟合值随时间的变化曲线。

图4是本发明实施实例1中混合气压缩因子随时间的变化曲线。

图5是本发明实施实例1中各时刻下的脱气稠油密度。

图6是本发明实施实例1中各时刻下的丙烷密度随时间的变化曲线。

图7是本发明实施实例1中各时刻下的二氧化碳密度随时间的变化曲线。

图8是本发明实施实例1中各时刻下的甲烷密度随时间的变化曲线。

图9是本发明实施实例1中各时刻下的稠油-丙烷二元相互作用系数β₄₁(t)随时间的变化曲线。

图10是本发明实施实例1中各时刻下的稠油-二氧化碳二元相互作用系数β₄₂(t)随时间的变化曲线。

图11是本发明实施实例1中各时刻下的稠油-甲烷二元相互作用系数β₄₃(t)随时间的变化曲线。

图12是本发明实施实例1中各时刻下的甲烷-丙烷二元相互作用系数β₃₁(t)随时间的变化曲线。

图13是本发明实施实例1中各时刻下的甲烷-二氧化碳二元相互作用系数β₃₂(t)随时间的变化曲线。

图14是本发明实施实例1中各时刻下的二氧化碳-丙烷二元相互作用系数β₂₁(t)随时间的变化曲线。

图15是本发明实施实例1中各时刻下油气界面的气体饱和浓度。

图16和图17分别是本发明实施实例1中各时刻下的混合气扩散系数与平均气相浓度。

图18和图19分别是本发明实施实例2中各时刻下的混合气扩散系数与平均气相浓度。

图20是本发明实施实例中气相浓度预测值与实验值对比结果图。

具体实施方式

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

在实施例中利用Orinoco油藏MPE3区块所产出的脱气稠油样品进行稠油注甲烷-二氧化碳-丙烷混合气扩散实验。通过扩散实验所得压力随实验时间变化数据，计算了各时刻下混合气扩散系数，体系密度等参数。并通过气相浓度实验值与气相浓度预测值进行对比，验证了本发明所述计算方法的可靠性。

实施例1

一种考虑稠油体积膨胀的甲烷-二氧化碳-丙烷混合气扩散系数的方法，计算模型原理如图1所示，计算流程图如图2所示，按如下步骤进行计算：

考虑稠油体积膨胀的甲烷-二氧化碳-丙烷混合气扩散系数的计算方法，包括如下步骤：

步骤1.确定油藏温度T与混合气的初始注入压力P(t₀)、混合气的初始气相高度h_g(t₀)、扩散实验容器高度H、油气界面面积S、稠油质量m_o；

各组分相对分子质量：M₁是丙烷的相对分子质量、M₂是二氧化碳的相对分子质量、M₃是甲烷的相对分子质量、M₄是稠油的相对分子质量；

混合气体的临界压缩因子Z_c；

混合气中各组分所占摩尔分数：y₁是丙烷所占摩尔分数、y₂是二氧化碳所占摩尔分数、y₃是甲烷所占摩尔分数；

气体组分性质参数：T_c1是丙烷的临界温度、T_c2是二氧化碳的临界温度、T_c3是甲烷的临界温度；P_c1是丙烷的临界压力、P_c2是二氧化碳的临界压力、P_c3是甲烷的临界压力；w₁是丙烷的偏心因子、w₂是二氧化碳的偏心因子、w₃是甲烷的偏心因子；

本实施例中，油藏温度T为54.2℃；初始注入压力P(t₀)为5MPa；初始气相高度h_g(t₀)为0.838m；扩散实验容器高度H为0.2538m；油气界面面积S为1×10^-3m²；稠油质量m_o为0.1668kg；混合气中丙烷的相对分子质量M₁为44.10；二氧化碳的相对分子质量M₂为44.01；甲烷的相对分子质量M₃为16.04；稠油的相对分子质量M₄为557.81；混合气体的临界压缩因子为0.284；混合气中丙烷的摩尔分数y₁为28mol％；二氧化碳的摩尔分数y₂为8mol％；甲烷的摩尔分数y₃为64mol％；混合气中丙烷的临界温度T_c1为369.80K；二氧化碳的临界温度T_c2为304.14K；甲烷的临界温度T_c3为190.55K；混合气中丙烷的临界压力T_c1为4246kPa；二氧化碳的临界压力T_c2为7378kPa；甲烷的临界压力T_c3为4600kPa；丙烷、二氧化碳和甲烷偏心因子w₁，w₂和w₃分别为0.008、0.224和0.152。

步骤2.进行稠油注甲烷-二氧化碳-丙烷混合气扩散实验：扩散容器内气体扩散示意图如图1所示，获取混合气压力P随实验时间t变化数据；

本实施例中实验过程中压力P随实验时间t变化如图3中的实验值所示。

步骤3.从实验总时长出发选择一个时间步长Δt，确定多个时刻t∈(t₀＝0；t₁＝t₀+Δt；......t_x＝t₀+xΔt)，从实验总时长出发选择一个时间步长Δt，确定多个时刻t，根据本发明所提出的气相压力公式(1)拟合实验数据，以计算气体扩散过程中各时刻下的气相压力P(t)；

P(t)＝B₁exp(-B₂t)+B₃exp(-B₄t)+B₅exp(-B₆t)+B₇ (1)

式中，B₁，B₂，B₃，B₄，B₅，B₆，B₇为公式系数，通过拟合实验数据得到；t为扩散时长，h；P(t)为时刻t时的气相压力，MPa。

本实施例中，实验总时长为360h，因此选择时间步长Δt为5h(1.8×10⁴s)，本次实验中，B₁＝0.2422，B₂＝27.6973，B₃＝0.2847，B₄＝0.0742，B₅＝0.9361，B₆＝0.0051，B₇＝3.5370，气相压力拟合值随时间的变化曲线如图3中的拟合值所示。

步骤4.根据气体压缩因子计算模型，即公式(2)～(4)计算在各时刻下的混合气体压缩因子Z(t)：

式中，Z(t)为各时刻下混合气压缩因子，无量纲；ρ_pr(t)为各时刻下的气体对应密度，无量纲；T为油藏温度，K；n为气体种数；P_cj为各气体组分的临界压力，kPa；T_cj为各气体组分的临界温度，K；Z_c为混合气临界压缩因子；A₁～A₁₁为公式常数，无量纲。

本实施例中，参数A₁～A₁₁如表1所示，通过牛顿迭代法计算各时刻下的混合气体压缩因子，结果如图4所示。

步骤5.根据以下脱气稠油密度公式，即公式(5)～(8)计算各时刻下的脱气稠油密度ρ₄(t)：

ρ_s＝a₁+a₂T+a₃T² (6)

A＝c₁+c₂T (7)

B＝b₁+b₂T^-1+b₃T^-2 (8)

式中，ρ₄(t)为各时刻下的脱气稠油密度，kg/m³；T为油藏温度，K；ρ_s为大气压下脱气稠油密度，kg/m³；A、B、a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂均为公式常数，无量纲，其具体数值在实施例给出。

表1公式(3)中参数值

本实施例中，参数A和B如表2所示，计算所得各时刻下的脱气稠油密度如图5所示。

步骤6.根据变形的PR状态方程，即公式(9)计算各时刻下的各气体参数C_j(t)值，进而得到各气体摩尔体积V_j(t)，再通过公式(18)计算在各时刻下稠油-甲烷-二氧化碳-丙烷体系中各气体组分密度ρ_j(t)：

a_j＝a_cjα(T_rj,ω_j)(j＝1,2,3) (13)

T_rj＝T/T_cj(j＝1,2,3) (16)

式中，A_j(t)、B_j(t)和C_j(t)为PR状态方程变形时所引入的参量；a_j为吸引力系数，无量纲；b_j为范德华摩尔体积，m³/kmol；V_j(t)为时刻t下各气体摩尔体积，m³/kmol；R为通用气体常数，8.314kPa·m³/(K·kmol)；T_rj为混合气中各气体组分的折算温度，无量纲；ω_j为各气体组分的偏心因子，无量纲；M_j为各气体组分的相对分子质量，g/mol；ρ_j(t)为时刻t下各气体组分的密度，kg/m³。

本实施例中，计算所得各时刻下的丙烷、二氧化碳和甲烷密度分别如图6-8所示。

步骤7.根据本发明拟合所得的公式(20)计算不同时刻下各气体组分间的二元相互作用系数，其中β₁₂为C₃H₈和CO₂二元相互作用系数、β₁₃为C₃H₈和CH₄二元相互作用系数、β₂₃为CO₂和CH₄二元相互作用系数；根据本发明拟合所得的公式(21)计算不同时刻下气体与稠油间的密度二元相互作用系数β_ij(t),其中β₁₄为C₃H₈和稠油二元相互作用系数、β₂₄为CO₂和稠油二元相互作用系数、β₃₄为CH₄和稠油二元相互作用系数：

式中，ρ_i(t)为时刻t下组分i的密度，kg/m³；ρ_j(t)为时刻t下组分j的密度，kg/m³；v_Nij(t)为时刻t下组分i和组分j的归一化比容差，无量纲；式(20)中的β_ij(t)为气体组分间的二元相互作用系数，无量纲；式(21)中的β_ij(t)为气体与稠油间的二元相互作用系数，无量纲。

表2公式(5)至(8)中公式参数

本实施例中，计算所得各时刻下的稠油-丙烷β₄₁(t)、稠油-二氧化碳β₄₂(t)、稠油-甲烷二元相互作用系数β₄₃(t)与甲烷-丙烷β₃₁(t)、甲烷-二氧化碳β₃₂(t)、二氧化碳-丙烷二元相互作用系数β₂₁(t)分别如图9-14所示。

步骤8.估计一个t₁时刻下气相高度h_g(t₁)的值，结合质量守恒原理与真实气体状态方程，按公式(24)计算t₁时刻下的溶解气质量m_d(t₁)：

M＝y₁M₁+y₂M₂+y₃M₃ (23)

式中，m_d(t₁)为t₁时刻下的气体扩散质量，kg；M为甲烷-二氧化碳-丙烷混合气体的相对分子质量；P(t₀)为初始气相压力，kPa；Z(t₀)为初始状态下气体压缩因子；h_g(t₀)为初始气相高度，m；P(t₁)为t₁时刻下的气体压缩因子；Z(t₁)为t₁时刻下的气体压缩因子；h_g(t₁)为t₁时刻下的气相高度，m；S为油气界面的面积，m²。

本实施例中，t₁时刻下气相高度h_g(t₁)估计值为0.0833m，按公式(24)计算t₁时刻下的溶解气质量m_d(t₁)为3.954×10^-4kg。

步骤9.根据超额体积混合准则，按照公式(30)计算在t₁时刻下的溶气稠油的密度ρ_mix(t₁)：

x₁(t₁)＝y₁n_g(t₁) (25)

x₂(t₁)＝y₂n_g(t₁) (26)

x₃(t₁)＝y₃n_g(t₁) (27)

x₄＝1-n_g(t₁) (28)

式中，n_g(t₁)为t₁时刻下溶气稠油中混合气的物质的量，mol；x_i(t₁)为t₁时刻下溶气稠油中组分i的物质的量，mol；n为混合体系中组分总数；ρ_mix(t₁)时刻下为t₁时刻下溶气稠油密度，kg/m³；w_i(t₁)为t₁时刻下组分i的质量分数，无量纲；ρ_i(t₁)为t₁时刻下组分i的密度，kg/m³；w_j(t₁)为t₁时刻下为组分j的质量分数，无量纲；ρ_j(t₁)为t₁时刻下组分j的密度，kg/m³；β_ij(t₁)为t₁时刻下组分i与组分j的密度二元相互作用系数。

本实施例中，按照公式(30)计算在t₁时刻下的溶气稠油的密度ρ_mix(t₁)的值为981.11kg/m³。

步骤10.计算t₁时刻下溶气稠油体积V(t₁)，进而求得t₁时刻下的气相高度h_g’(t₁)，比较步骤8的h_g(t₁)估算值与步骤10的h_g’(t₁)计算值：

若相对误差小于等于0.5×10^-2，则计算结果满足精度要求；

若相对误差大于0.5×10^-2，以步骤10所得h_g’(t₁)为新估算值，重复步骤8～10，直到相对误差满足精度要求，此时再利用h_g’(t₁)计算值按式(24)计算t₁时刻下的溶解气质量m_d(t₁)：

h_g,(t₁)＝H-h_t(t₁) (33)

式中，V₁(t₁)为t₁时刻下溶气稠油体积，m³；m_o为稠油质量，kg；h(t₁)为t₁时刻下的油相高度，m；h_g’(t₁)为t₁时刻下气相高度，m；H为扩散实验容器高度，m。

本实施例中，t₁时刻下气相高度h_g(t₁)计算值为0.0834m，因此相对误差为1.2×10^-3，满足精度要求。

步骤11.代入不同时刻下的混合气压力P(t)与实验时间t数据，重复步骤8～10所述计算方法，利用式(24)求出各时刻下的溶解气质量m_d(t_i)、利用式(30)求取溶气稠油密度ρ(t_i)，利用式(33)求取气相高度h_g’(t_i)(i＝2,3,…,n)；

步骤12.根据亨利定律，计算各时刻下油气界面的气体饱和浓度c_sat(t)：

式(34)中，P(t)为t时刻下的气相压力，kPa；K_j为第j种气体的亨利常数，kPa·m³/kg；c_j(t)为t时刻下第j种气体在油气界面的浓度，kg/m³；c_sat(t)为t时刻下混合气在油气界面的浓度，kg/m³；y_j为第j种气体的摩尔分数，mol％；式(35)中，n为混合气体中气体的种类。

本实施例中，K₁为34.8kPa·m³/kg，K₂为71.7kPa·m³/kg，K₃为607.6kPa·m³/kg，计算所得各时刻下油气界面的气体饱和浓度c_sat(t)结果如图15所示。

步骤13.如公式(36)与公式(37)所示，通过菲克定律与质量守恒原理来求解各时刻下的气体扩散系数；

进一步优选的，扩散实验开始时稠油中不含气相，因此求解初始条件如公式(38)所示；扩散过程中油气界面饱和气体，两相平衡，且气体扩散至容器底部时，无法继续向下扩散，因此上下边界条件分别如公式(39)与(40)所示；

进一步优选的，为便于计算，对公式(36)与(37)进行离散化处理，因此根据扩散实验容器高度，将稠油油藏划分为N个网格，网格高度按公式(44)计算。根据公式(41)～(43)，结合边界条件使用隐式差分法求解各时刻下的混合气扩散系数D^j与各网格中的气相浓度c^j _i，最后根据公式(45)求取稠油中的平均气相浓度c^j _p：

c|_t＝0＝0(0≤x＜h_t) (38)

式中，c为气体浓度，kg/m³；t为实验时间，s；D为扩散系数，m²/s，x为坐标位置，m；D^j为第j时刻下的气体扩散系数，m²/s；λ为差分过程所引入的变量，无物理意义；Δt为时间步长，s；Δx为空间步长，m；c^j+1 _i-1为第j+1时刻下第i-1网格中的气相浓度，kg/m³；c^j+1 _i为第j+1时刻下第i网格中的气相浓度，kg/m³；c^j+1 _i+1为第j+1时刻下第i+1网格中的气相浓度，kg/m³；c^j _i为第j时刻下第i网格中的气相浓度，kg/m³；c^j _sat为第j时刻下的油气界面气相浓度，kg/m³；c^j ₁为第j时刻下第1个网格中的气相浓度，kg/m³；h_g ^j+1为第j+1时刻下混合气相高度，m；P^{j +1}为第j+1时刻下混合气压力，m；Z^j+1为第j+1时刻下混合气压缩因子，无量纲；h_g ^j为第j时刻下混合气相高度，m；P^j为第j时刻下混合气压力，m；Z^j为第j时刻下混合气压缩因子，无量纲；c^j _p为第j时刻下的平均气相浓度，kg/m³；N为网格数；h_t ^j为第j时刻下油相高度，m。

本实施例中，计算所得的各时刻下的混合气扩散系数与平均气相浓度分别如图16和图17所示。

实施例2

一种考虑稠油体积膨胀的甲烷-二氧化碳-丙烷混合气扩散系数的计算方法。

如实施例1所述的一种考虑稠油体积膨胀的甲烷-二氧化碳-丙烷混合气扩散系数的计算方法，其不同之处在于改变了稠油注甲烷-二氧化碳-丙烷混合气扩散实验的初始注入压力，将初始注入压力分别设为8.65和10MPa进行实验，按照实施例1的计算步骤计算甲烷-二氧化碳-丙烷混合气扩散系数与平均气相浓度，结果如图18和图19所示。

效果例

为验证本发明所得计算结果的准确性，通过实验测定了中各时刻下稠油中平均气相浓度。平均气相浓度预测值与实验值如图20所示。从图20可以看出，气相浓度的预测值与实验值相差很小，从而验证了模型的准确性。

Claims (3)

1.一种考虑稠油体积膨胀的甲烷-二氧化碳-丙烷混合气扩散系数的计算方法，其特征在于，计算方法包括如下步骤：

步骤1.确定油藏温度T与混合气的初始注入压力P(t₀)、混合气的初始气相高度h_g(t₀)、扩散实验容器高度H、油气界面面积S、稠油质量m_o；

各组分相对分子质量：M₁是丙烷的相对分子质量、M₂是二氧化碳的相对分子质量、M₃是甲烷的相对分子质量、M₄是稠油的相对分子质量；

混合气体的临界压缩因子Z_c；

混合气中各组分所占摩尔分数：y₁是丙烷所占摩尔分数、y₂是二氧化碳所占摩尔分数、y₃是甲烷所占摩尔分数；

气体组分性质参数：T_c1是丙烷的临界温度、T_c2是二氧化碳的临界温度、T_c3是甲烷的临界温度；P_c1是丙烷的临界压力、P_c2是二氧化碳的临界压力、P_c3是甲烷的临界压力；w₁是丙烷的偏心因子、w₂是二氧化碳的偏心因子、w₃是甲烷的偏心因子；

步骤2.进行稠油注甲烷-二氧化碳-丙烷混合气扩散实验：获取混合气压力P随实验时间t变化数据；

步骤3.从实验总时长出发选择一个时间步长Δt，确定多个时刻t∈(t₀＝0；t₁＝t₀+Δt；......t_x＝t₀+xΔt)，根据气相压力公式(1)计算气体扩散过程中各时刻下的气相压力P(t)；

P(t)＝B₁exp(-B₂t)+B₃exp(-B₄t)+B₅exp(-B₆t)+B₇ (1)

式中，B₁，B₂，B₃，B₄，B₅，B₆，B₇为公式系数，通过拟合步骤2的实验数据得到；t为实验时间，h；P(t)为时刻t时的气相压力，MPa；

步骤4.根据气体压缩因子计算模型，即公式(2)～(4)计算在各时刻下的混合气体压缩因子Z(t)：

式中，Z(t)为各时刻下混合气压缩因子，无量纲；ρ_pr(t)为各时刻下的气体对应密度，无量纲；T为油藏温度，K；n为气体种数；P_cj为各气体组分的临界压力，kPa；T_cj为各气体组分的临界温度，K；Z_c为混合气临界压缩因子；A₁～A₁₁为公式常数，无量纲；

步骤5.根据以下脱气稠油密度公式，即公式(5)～(8)计算各时刻下的脱气稠油密度ρ₄(t)；

ρ_s＝a₁+a₂T+a₃T² (6)

A＝c₁+c₂T (7)

B＝b₁+b₂T^-1+b₃T^-2 (8)

式中，ρ₄(t)为各时刻下的脱气稠油密度，kg/m³；T为油藏温度，K；ρ_s为大气压下脱气稠油密度，kg/m³；A、B、a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂均为公式常数，无量纲；

步骤6.根据变形的PR状态方程，即公式(9)计算各时刻下的各气体参数C_j(t)值，进而得到各气体摩尔体积V_j(t)，再通过公式(18)计算在各时刻下稠油-甲烷-二氧化碳-丙烷体系中各气体组分密度ρ_j(t)：

a_j＝a_cjα(T_rj,ω_j) (j＝1,2,3) (13)

T_rj＝T/T_cj (j＝1,2,3) (16)

式中，A_j(t)、B_j(t)和C_j(t)为PR状态方程变形时所引入的参量；a_j为吸引力系数，无量纲；b_j为范德华摩尔体积，m³/kmol；V_j(t)为时刻t下各气体摩尔体积，m³/kmol；R为通用气体常数，8.314kPa·m³/(K·kmol)；T_rj为混合气中各气体组分的折算温度，无量纲；ω_j为各气体组分的偏心因子，无量纲；M_j为各气体组分的相对分子质量，g/mol；ρ_j(t)为时刻t下各气体组分的密度，kg/m³；

步骤7.根据公式(20)计算不同时刻下各气体组分间的二元相互作用系数β_ij(t)，其中β₁₂为C₃H₈和CO₂二元相互作用系数、β₁₃为C₃H₈和CH₄二元相互作用系数、β₂₃为CO₂和CH₄二元相互作用系数；根据公式(21)计算不同时刻下气体与稠油间的密度二元相互作用系数β_ij(t)，其中β₁₄为C₃H₈和稠油二元相互作用系数、β₂₄为CO₂和稠油二元相互作用系数、β₃₄为CH₄和稠油二元相互作用系数：

式中，ρ_i(t)为时刻t下组分i的密度，kg/m³；ρ_j(t)为时刻t下组分j的密度，kg/m³；v_Nij(t)为时刻t下组分i和组分j的归一化比容差，无量纲；式(20)中的β_ij(t)为气体组分间的二元相互作用系数，无量纲；式(21)中的β_ij(t)为气体与稠油间的二元相互作用系数，无量纲；

步骤8.估计一个t₁时刻下气相高度h_g(t₁)的值，结合质量守恒原理与真实气体状态方程，计算t₁时刻下的溶解气质量m_d(t₁)；

步骤9.根据超额体积混合准则，按照公式(30)计算在t₁时刻下的溶气稠油的密度ρ_mix(t₁)；

x₁(t₁)＝y₁n_g(t₁) (25)

x₂(t₁)＝y₂n_g(t₁) (26)

x₃(t₁)＝y₃n_g(t₁) (27)

x₄(t₁)＝1-n_g(t₁) (28)

式中，n_g(t₁)为t₁时刻下溶气稠油中混合气的物质的量，mol；m_d(t₁)为t₁时刻下的气体扩散质量，kg；x_i(t₁)为t₁时刻下溶气稠油中组分i的物质的量，mol；n为混合体系中组分总数；ρ_mix(t₁)时刻下为t₁时刻下溶气稠油密度，kg/m³；w_i(t₁)为t₁时刻下组分i的质量分数，无量纲；ρ_i(t₁)为t₁时刻下组分i的密度，kg/m³；w_j(t₁)为t₁时刻下为组分j的质量分数，无量纲；ρ_j(t₁)为t₁时刻下组分j的密度，kg/m³；β_ij(t₁)为t₁时刻下组分i与组分j的密度二元相互作用系数；

步骤10.计算t₁时刻下溶气稠油体积V(t₁)，进而求得t₁时刻下的气相高度h_g’(t₁)，比较步骤8的h_g(t₁)估算值与步骤10的h_g’(t₁)计算值：

若相对误差小于等于0.5×10^-2，则计算结果满足精度要求；

若相对误差大于0.5×10^-2，以步骤10所得h_g’(t₁)为新估算值，重复步骤8～10，直到相对误差满足精度要求，此时再利用h_g’(t₁)计算值按式(24)计算t₁时刻下的溶解气质量m_d(t₁)；

h_g,(t₁)＝H-h_t(t₁) (33)

式中，V₁(t₁)为t₁时刻下溶气稠油体积，m³；m_o为稠油质量，kg；h(t₁)为t₁时刻下的油相高度，m；h_g’(t₁)为t₁时刻下气相高度，m；H为扩散实验容器高度，m；

步骤11.代入不同时刻下的混合气压力P(t)与实验时间t数据，重复步骤8～10所述计算方法，利用式(24)求出各时刻下的溶解气质量m_d(t_i)、利用式(30)求取溶气稠油密度ρ(t_i)，利用式(33)求取气相高度h_g’(t_i)(i＝2,3,…,n)；

步骤12.根据亨利定律，计算各时刻下油气界面的气体饱和浓度c_sat(t)；

步骤13.扩散实验开始时求解初始条件如公式(38)所示；扩散过程中上下边界条件分别如公式(39)与(40)所示；

对菲克定律(即公式(36))与质量守恒原理(即公式(37))进行离散化处理，因此根据扩散实验容器高度，将稠油油藏划分为N个网格，网格高度按公式(44)计算；根据公式(41)～(43)，结合边界条件使用隐式差分法求解各时刻下的混合气扩散系数D^j与各网格中的气相浓度c^j _i，最后根据公式(45)求取稠油中的平均气相浓度

c|_t＝0＝0 (0≤x＜h_t) (38)

式中，c为气体浓度，kg/m³；t为实验时间，s；D为扩散系数，m²/s，x为坐标位置，m；D^j为第j时刻下的气体扩散系数，m²/s；λ为差分过程所引入的变量，无物理意义；Δt为时间步长，s；Δx为空间步长，m；c^j+1 _i-1为第j+1时刻下第i-1网格中的气相浓度，kg/m³；c^j+1 _i为第j+1时刻下第i网格中的气相浓度，kg/m³；c^j+1 _i+1为第j+1时刻下第i+1网格中的气相浓度，kg/m³；c^j _i为第j时刻下第i网格中的气相浓度，kg/m³；c^j _sat为第j时刻下的油气界面气相浓度，kg/m³；c^j ₁为第j时刻下第1个网格中的气相浓度，kg/m³；h_g ^j+1为第j+1时刻下混合气相高度，m；P^{j +1}为第j+1时刻下混合气压力，m；Z^j+1为第j+1时刻下混合气压缩因子，无量纲；h_g ^j为第j时刻下混合气相高度，m；P^j为第j时刻下混合气压力，m；Z^j为第j时刻下混合气压缩因子，无量纲；c^j _p为第j时刻下的平均气相浓度，kg/m³；N为网格数；h_t ^j为第j时刻下油相高度，m。

2.如权利要求1所述的一种考虑稠油体积膨胀的甲烷-二氧化碳-丙烷混合气扩散系数的计算方法，其特征在于，所述的步骤8按照以下公式计算t₁时刻下的溶解气质量m_d(t₁)；

M＝y₁M₁+y₂M₂+y₃M₃ (22)

式中，m_d(t₁)为t₁时刻下的气体扩散质量，kg；M₁、M₂、M₃分别为甲烷、二氧化碳、丙烷混合气体的相对分子质量，kg/kmol，y₁是丙烷所占摩尔分数、y₂是二氧化碳所占摩尔分数、y₃是甲烷所占摩尔分数；P(t₀)为初始气相压力，kPa；Z(t₀)为初始状态下气体压缩因子；h_g(t₀)为初始气相高度，m；P(t₁)为t₁时刻下的气体压缩因子；Z(t₁)为t₁时刻下的气体压缩因子；h_g(t₁)为t₁时刻下的气相高度，m；S为油气界面的面积，m²。

3.如权利要求1所述的一种考虑稠油体积膨胀的甲烷-二氧化碳-丙烷混合气扩散系数的计算方法，其特征在于，所述的步骤12中的c_sat(t)通过以下公式计算：

式中，P(t)为t时刻下的气相压力，kPa；K_j为第j种气体的亨利常数，kPa·m³/kg；c_j(t)为t时刻下第j种气体在油气界面的浓度，kg/m³；y_j为第j种气体的摩尔分数，mol％；式中，n为混合气体中气体的种类。

Images