引言

目前，中国二氧化碳排放量居世界第二位，并仍在快速增长，预计到2025年前后，中国的二氧化碳排放总量会超过美国，居世界第一位。另一方面，中国石油供求矛盾日益突出，2010年以后进口依赖度始终约在50%以上。中国大陆东部地区老油气田开采区是封闭条件非常好的地下储气库，向油层中注入二氧化碳可以实现温室气体的长期埋存，同时还能大幅度提高原油的采收率^[1-2]。数值模拟是油藏二氧化碳驱项目方案设计、预测和管理的主要技术手段。由于混相驱或非混相驱过程中二氧化碳同原油和地下水的复杂相互作用，常用的黑油模型已不适用，必须采用组分模型进行计算分析。

本文的目的是建立裂缝性油藏二氧化碳驱的多相多组分数值计算模型，给出相关的计算方法和主要公式；这一计算模型主要包括裂缝性油藏二氧化碳驱多相多组分的基本控制方程、二氧化碳的吸附及原油重组分的沉淀、多组分气相和油相物性的计算、沥青质沉积对渗透率的影响、油-水-气(二氧化碳)三相地下条件下的相平衡计算方法等。

裂缝性油藏数值模拟最常采用的理论是双孔隙度-单渗透率、双孔隙度-双渗透率等双重介质模型^[3-5]。从双重介质的概念出发，很多学者都提出了相应的三重介质模型，其中，基于多重相关连续介质的多重介质模型灵活有效、适用性广^[6-8]，本文以此多重介质模型为基础，研究裂缝性油藏数值模拟的组分模型。

多相多组分流体在地下条件下的相平衡计算是组分模型的核心问题之一，主要有平衡常数法、基于Rachford-Rice函数的闪蒸计算法和Gibbs自由能最小化法3种计算方法。本文分别给出了采用这3种方法计算油-水-气(二氧化碳)三相流体相平衡的主要计算公式，并简要分析了这3种相平衡计算方法的优缺点。

1 基本控制方程

注入的二氧化碳在油藏条件下通常呈超临界状态(以下简称为气相)。假设二氧化碳驱过程中等温、不发生化学反应、且孔隙中相平衡能够瞬时完成，则多相多组分流体的流动由质量守恒定律来控制。

1.1 质量守恒方程

裂缝性油藏可用裂缝-基质的双重介质理论来描述，裂缝和基质同为储集空间，裂缝为主要的流动通道。根据文献[6-8]的多重介质理论，双重介质中裂缝间的流动、基质间的流动、裂缝-基质间的多相流动，都可抽象地归纳为单元间的流动，从而采用统一的质量守恒表达式。同时，这一理论还能够考虑不同的裂缝分布模式，如图 1所示的裂缝分布，可以抽象为图 2所示的概念模型，并根据此概念模型确定裂缝系统和基质系统间的窜流系数。其他的分布模式可参考文献[8]。

本文考虑的流体包括气相(二氧化碳)、水相和油相，每相流体都包含有若干质量组分。多孔介质中单元n中组分κ的质量守恒方程为

$ \frac{{\rm{d}}}{{{\rm{d}}t}}\int\limits_{{V_n}} {{M^\kappa }} {\rm{d}}{V_n} = \int\limits_{{\mathit{\Gamma} _n}} {{{\vec F}^\kappa } \cdot \vec n} {\rm{d}}{\mathit{\Gamma} _n} + \int\limits_{{V_n}} {{q^\kappa }} {\rm{d}}{V_n} $

(1)

孔隙介质中各相流体组分κ的质量累积为

$ {M^\kappa } = \phi \sum\limits_\beta {{S_\beta }} {\rho _\beta }x_\beta ^\kappa + M_{{\rm{ads}}}^\kappa $

(2)

1.2 多相流动及组分扩散

多孔介质中多相流动用扩展的Darcy定律来描述，组分κ因渗流引起的质量传输项为

$ {{\vec F}^\kappa }{{\rm{|}}_{{\rm{adv}}}} = \sum\limits_\beta {x_\beta ^\kappa {{\vec F}_\beta }} $

(3)

$ {{\vec F}_\beta } = {\rho _\beta }{{\vec u}_\beta } = - K\frac{{{K_{{\rm{r}}\beta }}{\rho _\beta }}}{{{\mu _\beta }}}\left( {\nabla {p_\beta } - {\rho _\beta }{\rm{\tilde g}}} \right) $

(4)

组分κ(如二氧化碳等)在多孔介质中的扩散采用扩展的Fick定律来描述，相β中组分κ因扩散引起的质量传输项为

$ {\left. {F_\beta ^\kappa } \right|_{{\rm{dis}}}} =-\phi {\rho _\beta }{\rm{D}}_\beta ^\kappa \nabla x_\beta ^\kappa $

(5)

1.3 二氧化碳的吸附及原油重组分的沉淀

二氧化碳在岩石或黏土固体颗粒上的吸附是一个动态的过程，假设吸附完成的速度相对于油藏中其他物理过程来说足够快，则这一吸附过程可用瞬时、可逆、等温吸附来描述

$ M_{{\rm{ads}}}^\kappa = (1-\phi ){\rho _{\rm{s}}}K_{{\rm{d}}, \beta }^\kappa {\rho _\beta }x_\beta ^\kappa $

(6)

二氧化碳注入后，原油中重组分的沉淀会对储层渗透性能造成严重影响。重组分的沉淀会阻塞孔隙喉道、降低渗透率，并可能会使岩石的湿润性由亲水转变为亲油。为了考虑这种情况并避免进行复杂的四相混合物相平衡计算，本文假设重组分的沉淀可用式(6)的等温吸附来描述，采用非线性的Langmuir等温吸附时，吸附平衡系数K_{d, β}^κ是相β中二氧化碳摩尔分数的函数，可由试验来确定。

2 本构关系和流体的物性

本构关系用来确定多相多组分流体流动各个物理过程中的变量和参数之间的内在关系和约束。根据本构关系，可由主变量来计算各相流体的物理性质。

2.1 本构关系

多相多组分流动体系中，各相流体饱和度及各组分的摩尔分数应当满足如下的约束方程

$ \sum\limits_\beta {{{\rm{S}}_\beta }} = 1 $

(7)

$ \sum\limits_\kappa {x_\beta ^\kappa } = 1 $

(8)

气-水-油的三相体系中，毛细管力由流相的饱和度来确定。如果选择气相为参考相，则水相和油相的压力计算公式为

$ {p_{\rm{w}}} = {p_{\rm{g}}}-{p_{{\rm{cgw}}}}\left( {{S_{\rm{w}}}, {\mkern 1mu} {S_{\rm{g}}}} \right) $

(9)

$ {p_{\rm{o}}} = {p_{\rm{g}}}-{P_{{\rm{cgo}}}}\left( {{S_{\rm{w}}}, {\mkern 1mu} {S_{\rm{o}}}} \right) $

(10)

气水毛细管力p_cgw和气油毛细管力p_cgo是水相饱和度S_w和油相饱和度S_o的函数。油水毛细管力的计算公式为

$ {p_{{\rm{cow}}}} = {p_{{\rm{cgw}}}}-{p_{{\rm{cgo}}}} = {p_{\rm{o}}}-{p_{\rm{w}}} $

(11)

气-水-油三相体系的三相相对渗透率对于准确预测多相流体流动动态十分重要。以往的数值计算模型中，通常基于这样的假定：气相为完全非润湿相，三相体系中一液相为完全的润湿相，另一液相介于二者之间为中间润湿相；用完全非润湿相和完全润湿相的两相相对渗透率作为它们各自的三相相对渗透率，中间润湿相的三相相对渗透率则由一些简单模型计算而得，如Stone Ⅱ方法等^[9-11]。

然而，Dria D E和Pope G A的研究结果表明^[12]，对于二氧化碳驱条件下的气-水-油三相体系，各相的三相相对渗透率并不是上述计算方法中依赖于两相的饱和度，而是仅与该相的饱和度有关。据此，相β的三相相对渗透率可表示为该相饱和度的函数

$ {K_{{\rm{r}}{\kern 1pt} \beta }} = {\mkern 1mu} {\mkern 1mu} {K_{{\rm{r}}{\kern 1pt} \beta }}\left( {{S_\beta }} \right) $

(12)

2.2各相流体的物性

二氧化碳组分在不同温度和压力下的密度采用文献[13]给出的关系式，不同温度和压力下的黏度可采用文献[14]中的数据。气相的主要组分为二氧化碳和天然气中烃类气体。气相的PVT性质，采用Soave-Redlich-Kwong (SRK)状态方程来描述^[15-16]。纯组分的SRK状态方程为

$ {Z^3}- {Z^2} + \left[{{A^*}-{B^*}-{{({B^*})}^2}} \right]Z -{A^*}{B^*} = 0 $

(13)

其中，$ {A^*} = ap/({{\rm{R}}^2}{T^2}) $和$ {B^*} = bp{\rm{/R}}T $是取决于压力p和热力学温度T的两个参数，参数a和b的计算公式为^[17]

$ a = {\left\{ {1 + \left( {0.48 + 1.574\omega- 0.176{\omega ^2}} \right)\left[{1-{{\left( {T/{T_c}} \right)}^{1/2}}} \right]} \right\}^2} \cdot \frac{{0.42748{{\rm{R}}^2}{T_{\rm{c}}}^2}}{{{p_{\rm{c}}}}} $

(14)

$ b = \frac{{0.08664{\rm{R}}{T_{\rm{c}}}}}{{{p_{\rm{c}}}}} $

(15)

对于多组分的气相，则可采用文献[17]提出的混合规则进行计算，混合物的a、b系数设为a_m和b_m，可以由纯组分的参数求得，设纯组分κ的a和b系数分别是a_κ和b_κ，则计算公式为

$ {a_{\rm{m}}} = {\left( {\sum\limits_\kappa {x_{\rm{g}}^\kappa \sqrt {{a_\kappa }} } } \right)^2} $

(16)

$ {b_{\rm{m}}} = \sum\limits_\kappa {x_{\rm{g}}^\kappa } {b_\kappa } $

(17)

根据式(13)得到气相的压缩因子后，其摩尔密度为

$ {\rho _{\rm{g}}} = p{\rm{R}}TZ $

(18)

纯组分κ的黏度μ_κ是压力和温度的函数。根据Chapman-Enskog动力学理论，气相的动力黏度μ_g可由纯组分的动力黏度μ_κ和组成表述

$ {\mu _{\rm{g}}} = \sum\limits_{\kappa = 1}^N {\frac{{{x_{\rm{g}}}^\kappa {\mu _\kappa }}}{{\sum\nolimits_{\lambda = 1}^N {{x_{\rm{g}}}^\lambda {\mathit{\Phi} _{\kappa \lambda }}} }}} $

(19)

Φ_κλ是二元相互作用参数，表述为

$ {\mathit{\Phi} _{\kappa \lambda }} = \frac{{{{\left[ {1 + {{\left( {{\mu _\kappa }/{\mu _\kappa }/{\mu _\lambda }} \right)}^{1/2}}{{\left( {{M_\lambda }/{M_\lambda }} \right)}^{1/4}}} \right]}^2}}}{{{{\left[ {8\left( {1 + {M_\kappa }/{M_\lambda }} \right)} \right]}^{1/2}}}} $

(20)

对于油相，摩尔密度的计算公式为

$ {\rho _{\rm{o}}} = p{\rm{R}}T{Z_{\rm{o}}} $

(21)

压缩因子Z_o根据原油的组成成分按照式(13)计算。

原油纯组分的黏度计算公式为

$ \ln \mu = A'-\frac{{B'}}{T} + C'T + D'{T^2} $

(22)

油相的黏度可由其组成按照式(23)计算

$ {\mu _{\rm{o}}} = \prod\limits_{\kappa = 1}^{{N_{\rm{c}}}} {{\mu _\kappa }^{x_{\rm{o}}^\kappa }} $

(23)

对于水的摩尔密度和黏度，则根据文献[18]的关系式和数据来确定。由于水相中二氧化碳和烃类溶解度相对较小，因此假定溶解于水的二氧化碳和烃类对水相本身的摩尔密度和黏度影响不大。

2.3 沥青质沉积对渗透率的影响

沥青质沉积导致的渗透率降低不仅取决于孔隙减少量，还取决于孔隙结构和沉积物质在孔隙中的分布。如图 3所示^[19]，可采用文献[19]的管串联模型，将孔隙结构考虑为由孔隙和喉道串联构成的毛管，R是孔隙的半径，r是喉道的半径，沥青质在串联毛管内壁均匀沉积，渗透率主要由喉道控制。

考虑沥青质沉积后的渗透率按式(24)式(26)计算

$ \frac{K}{{{K_0}}} = {\theta ^2}\frac{{1- {\rm{ }}\mathit{\Gamma} + {\rm{ }}\mathit{\Gamma} /{w^2}}}{{1- \mathit{\Gamma} + \mathit{\Gamma} {{\left[{\theta /\left( {\theta + w-1} \right)} \right]}^2}}} $

(24)

$ \theta = \frac{{1-{S_{\rm{s}}}-{\phi _r}}}{{1-{\phi _r}}} $

(25)

$ w = 1 + \frac{{1/\mathit{\Gamma} }}{{1/{\phi _r}-1}} $

(26)

3 相平衡计算

无论是混相二氧化碳驱油还是非混相二氧化碳驱油，其形成过程都是随着油层温度、压力和原油组成的变化而产生的相态平衡变化过程。加滕完司^[20]对细管二氧化碳驱油的研究表明：在温度、相对渗透率和黏度相同、压力相差不是很大的(10.3 MPa和12.4 MPa)的情况下，二氧化碳同原油的相平衡特性有着十分显著的差异；原油组分的不同，二氧化碳平衡常数变化明显。可见，准确的相平衡计算对于油藏二氧化碳驱的数值模拟至关重要。已有的一些组分模型数值计算方法需要在相平衡计算方面进行改进^[21-25]。多相多组分流体在地下条件下的相平衡计算主要有平衡常数法、基于Rachford-Rice函数的闪蒸计算法和Gibbs自由能最小化法3种计算方法。

3.1 平衡常数法

对于多组分体系内的相α和相β，组分κ的平衡常数定义为

$ K_{\alpha \beta }^\kappa = \frac{{x_\alpha ^\kappa }}{{x_\beta ^\kappa }} $

(27)

这一平衡常数应该是多组分混合物温度、压力和组成的函数。对于烃类体系，文献[26]根据试验建议了关系式近似计算气-油平衡常数

$ K_{{\rm{go}}}^i = \left( {\frac{{{p_{{\rm{c}}, i}}}}{p}} \right)\exp \left[{5.373\left( {1 + {\omega _i}} \right)\left( {1-\frac{{{T_{{\rm{c}}, i}}}}{T}} \right)} \right] $

(28)

在此基础上，Almehaideb R A等[27]根据大量的试验结果，建议了一个适用于高压体系的平衡常数计算公式

$ K_{{\rm{go}}}^i = {\left( {\frac{{{p_{{\rm{c}}, i}}}}{{{p_{\rm{k}}}}}} \right)^{A-1}}\left( {\frac{{{p_{{\rm{c}}, i}}}}{p}} \right)\exp \left( {AK_i^*} \right) $

(29)

其中：$ K_i^* = \frac{{{a_{{\rm{T1}}}}}}{{{T^2}}} + \frac{{{a_{{\rm{T2}}}}}}{T} + {a_{{\rm{T3}}}} + {a_{{\rm{p1}}}}\ln p + \frac{{{a_{{\rm{p2}}}}}}{{{p^2}}} + \frac{{{a_{{\rm{p3}}}}}}{p} + \frac{{{a_\omega }}}{\omega } $；$ A = 1-{\left( {\frac{{p-14.7}}{{{p_{\rm{k}}}-14.7}}} \right)^{0.6}} $。

文献[28]基于多重介质模型，采用平衡常数法研究了二氧化碳驱多相多组分数值模拟方法，并编制了相应的计算程序。

平衡常数法采用基于试验的经验公式来计算平衡常数，能够避免复杂、耗时的相平衡计算，计算速度快、稳定性良好，但由于缺乏严密的热力学理论基础，难以很好地考虑原油不同组成的影响，计算结果可能同真实相平衡状态有一定的偏差。

3.2 基于Rachford-Rice函数的闪蒸计算法

以气油两相体系为例，基于Rachford-Rice方法^[29]的闪蒸计算的迭代步骤如下

(1) 估计组分i初始的油气平衡常数K_goⁱ和气相的摩尔分数n_v；

(2) 按照式(30)~式(32)计算Rachford-Rice函数和新的气相摩尔分数(n_v)_new；

$ f\left( {{n_{\rm{v}}}} \right) = \sum {\frac{{{z_i}\left( {K_{{\rm{go}}}^i-1} \right)}}{{1 + \left( {K_{{\rm{go}}}^i-1} \right){n_{\rm{v}}}}}} $

(30)

$ f'\left( {{n_{\rm{v}}}} \right) =- \sum {\frac{{{z_i}{{\left( {K_{{\rm{go}}}^i- 1} \right)}^2}}}{{{{\left[{1 + \left( {K_{{\rm{go}}}^i-1} \right){n_{\rm{v}}}} \right]}^2}}}} $

(31)

$ {\left( {{n_{\rm{v}}}} \right)_{{\rm{new}}}} = {n_{\rm{v}}}-\frac{{f\left( {{n_{\rm{v}}}} \right)}}{{f'\left( {{n_{\rm{v}}}} \right)}} $

(32)

(3) 迭代计算至n_v同(n_v)_new的误差满足允许的要求；

(4) 计算组分i在气相中的摩尔分数y_i和油相中的摩尔分数x_i；

$ {x_i} = \frac{{{z_i}}}{{1-{n_{\rm{v}}} + {n_{\rm{v}}}K_{{\rm{go}}}^i}} $

(33)

$ {y_i} = \frac{{{z_i}K_{{\rm{go}}}^i}}{{1-{n_{\rm{v}}} + {n_{\rm{v}}}K_{{\rm{go}}}^i}} $

(34)

(5) 根据气相和油相的摩尔组成分别计算气相中组分i的逸度系数φ_gⁱ和油相中组分i的逸度系数φ_oⁱ；

(6) 计算新的油气平衡常数，返回步骤1进行迭代，直至新的平衡常数满足误差要求。

$ K_{{\rm{go}}}^i = \frac{{\varphi _{\rm{o}}^{\rm{i}}}}{{\varphi _{\rm{g}}^i}} $

(35)

上述闪蒸计算方法同基于经验公式的平衡常数法相比较，有较高的精度，但是因为需要反复迭代，因此计算量较平衡常数法大。组分模拟迭代求解过程中，当式(30)右项分母接近于零时，则该算法出现不收敛的情况。此外，基于Rachford-Rice函数的闪蒸计算法主要用于两相体系的相平衡计算，适用于三相体系相应的算法还有待研究和完善^[30]。

3.3 Gibbs自由能最小化法

根据多相体系相平衡的热力学理论，多组分混合物处于稳定状态时应满足Gibbs自由能最小的条件。油-水-气三相体系总的Gibbs自由能为

$ {n^{\rm{T}}}G = \sum\limits_i {{n^i}G_0^i} + {\rm{R}}T\sum\limits_\beta {\sum\limits_i {\left[{{n^{\rm{T}}}L_\beta ^{}x_\beta ^i\ln \left( {\varphi _\beta ^ix_\beta ^ip} \right)} \right]} } $

(36)

相β中组分i的逸度系数可由状态方程及式(37)计算得到

$ \ln \varphi _\beta ^i = {\left\{ {\frac{1}{{{\rm{R}}T}}\int_V^\infty {\left[{{{\left( {\frac{{\partial P}}{{\partial {n^j}}}} \right)}_{T, V, {n^{j \ne i}}}}-\frac{{{\rm{R}}T}}{V}} \right]} {\rm{d}}V -\ln Z} \right\}_\beta } $

(37)

由于体系达到相平衡时，体系的摩尔Gibbs自由能应达到最小，因此，油-水-气三相的平衡可通过搜索式(36)中摩尔Gibbs自由能最小值来进行。Gibbs自由能最小化方法理论严密、适应性广、计算稳定、收敛性好，比较适合于多相平衡计算，但是计算量比基于Rachford-Rice函数的闪蒸计算法要大。

4 数值解法

对于式(1)的质量守恒方程，采用文献[6]，文献[7]和文献[28]中的积分差分法进行空间离散、采用一阶向后差分进行时间离散，可写成残差形式

$ R_n^{\kappa, t + \Delta t} = M_n^{\kappa, t + \Delta t}-M_n^{\kappa, t}-\frac{{\Delta t}}{{{V_n}}}\left\{ {\sum\limits_m {{A_{nm}}F_{nm}^{\kappa, t + \Delta t}} + {V_n}q_n^{\kappa, t + \Delta t}} \right\} = 0 $

(38)

可采用Newton-Raphson迭代法对方程(38)进行求解。

5 结语

裂缝性油藏二氧化碳驱的多相多组分数值计算方法主要涉及如下问题的计算模拟：裂缝性油藏的多重介质模型、多相流动及多组分的扩散、二氧化碳在岩石固体骨架上的吸附、原油重组分的沉淀、储层条件下不同组分流体物性的计算方法、沥青质沉积对渗透率的影响、多相多组分的相平衡计算方法和数值离散方法。

本文研究了裂缝性油藏二氧化碳驱的油-水-气(二氧化碳)三相、多组分的数值计算模型，较为全面地给出了上述相关问题的计算方法和主要公式。对于裂缝性油藏可采用文献[6-8]和[28]的多重介质理论来描述，以便考虑不同的裂缝分布模式。裂缝系统内、基质系统内、裂缝-基质之间的多相流动，采用扩展的Darcy定律来描述。二氧化碳在油相中的扩散可采用Fick定律来描述。采用非线性的等温吸附来描述二氧化碳在岩石颗粒上的吸附及油相中沥青质在岩石骨架上的沉淀。沥青质沉积阻塞孔道对储层绝对渗透率的影响，采用串状毛管模型来考虑。各组分的质量守恒方程采用积分差分法进行空间离散，可采用全隐式格式通过Newton迭代法进行求解。

相平衡是二氧化碳驱组分模型的核心问题之一。不同压力及组成下二氧化碳同原油相平衡特性的差异十分显著，准确的相平衡计算对于油藏二氧化碳驱的数值模拟至关重要。对于油-水-气(二氧化碳)三相的相平衡计算，分别给出了平衡常数法、基于Rachford-Rice函数的闪蒸计算法和Gibbs自由能最小化法的计算方法，并简述了3种方法各自的优点和不足。