引言

CO$ _2 $驱开发过程中，需针对地下CO$ _2 $-原油体系相态进行研究^[1-3]。CO$ _2 $与原油之间的相态计算是注气机理研究和流体特征计算的关键，由实验数据给出的相图虽能直观给出油气烃类体系的相态特征变化，但仅依靠相态实验测试，会受实验仪器工作温度、压力范围限制，不能得到完整相图和全部参数^[4-6]。通过流体相平衡计算来描述和预测CO$ _2 $-原油体系的相态特征和变化规律，是现有常用的研究手段^[7]。

从1802年提出的理想气体状态方程至今，状态方程已有长足发展，广泛应用于体系相态计算中。目前，按应用形式和理论基础，现有状态方程可以分为van der Waals(vdW)型及其改进型^[8-11]、varial型^[12-13]、多参数型^[14-16](BWR、BWRS等)以及具有严格统计力学基础的类型(微扰硬链和转子链等)^[17-18]。

以上方程虽取得长足发展，但仍然存在不足。(1) vdW型方程对高温高压条件考虑不足，影响实用。(2) varial型方程由于参数过多且不易通过实验获取，只能常用二次或三次截断型，适用范围较低。(3)多参数型方程在非烃气体、重组分及较低温度条件下，计算精度不太令人满意。(4)具有严格统计力学基础的状态方程形式比较复杂，工程应用较少。

CO$ _2 $的分子为线性四偶极距分子^[19]，同烃类作用时具有特殊性，特别在临界状态^[20]，相间发生严重组分交换会放大预测方法本身的缺陷，导致状态方程对CO$ _2 $-烃体系预测精度偏低。因此，需要对现有状态方程参数进行修正，得到精度高且计算简便的方程形式。

1 修正的PR状态方程

PR是应用广泛的vdW型方程之一，能满足vdW型方程所具有的临界点条件，即一种物质在临界点时，临界等温线上的压力对体积的一阶和二阶导数等于零^[21-23]。

$ \left\{ \begin{array}{l} {\left( {\dfrac{{\partial p}}{{\partial v}}} \right)_{{T_{ \rm c }}}} = {\left( {\dfrac{{{\partial ^2}p}}{{\partial {v^2}}}} \right)_{{T_{ \rm c }}}}{ \rm = }0\\[10pt] p = \dfrac{{{ \rm R }T}}{{v - {b_i}}} - \dfrac{{{a_i} \alpha \left( T \right)}}{{v\left( {v + {b_i}} \right) + {b_i}\left( {v - {b_i}} \right)}} \end{array} \right. $

(1)

则可以得到式中引力修正系数$ a_i $和体积修正系数$ b_i $

$ \left\{ \begin{array}{l} {a_i} = {\varOmega _{ \rm a }}\dfrac{{{{ \rm R }^2}T_{{ \rm c }i}^2}}{{{p_{{ \rm c }i}}}}\\[10pt] {b_i} = {\varOmega _{ \rm b }}\dfrac{{{ \rm R }{T_{{ \rm c }i}}}}{{{p_{{ \rm c }i}}}} \end{array} \right. $

(2)

$ \alpha \left( T \right) = {\left[ {1 + {m_i}\left( {1 - T_{{ \rm r }i}^{0.5}} \right)} \right]^2} $

(3)

$ {m_i} = 0.37464 + 1.48503{\omega _i} - 0.26992\omega _i^2 $

(4)

$ {\varOmega _{ \rm a }} = 0.45724, {\varOmega _{ \rm b }} = 0.07780 $

(5)

PR方程中，无因次系数通过实验回归得到，回归系数应与组分具有相关性而非常数；同时，体积和引力与温度具有较强相关性，无因次系数也应与温度具有较强相关性。因此，对每个组分，将$ \varOmega _{\rm{a}} $和$ \varOmega _{\rm{b}} $处理成温度的函数，为了保持其普适性和无因次形式，利用物质对比态相似性原理和对比态无因次化特性，关联为对比温度的函数。

$ \begin{eqnarray} \left\{ \begin{array}{l} {\varOmega _{{\rm{a}}{T_{{\rm{r}}i}}}} = f\left( {{T_{{\rm{r}}i}}} \right)\\ {\varOmega _{{\rm{b}}{T_{{\rm{r}}i}}}} = f\left( {{T_{{\rm{r}}i}}} \right) \end{array} \right. \end{eqnarray} $

(6)

假设这些与温度有关的参数是对比温度函数的二次方多项式。将式(6)代入式(2)，则

$ \begin{eqnarray} \left\{ \begin{array}{l} {a_i} = {\varOmega _{{\rm{a}}{T_{{\rm{r}}i}}}}\dfrac{{{{\rm{R}}^2}T_{{\rm{c}}i}^2}}{{{p_{{\rm{c}}i}}}}\\[10pt] {b_i} = {\varOmega _{{\rm{b}}{T_{{\rm{r}}i}}}}\dfrac{{{\rm{R}}{T_{{\rm{c}}i}}}}{{{p_{{\rm{c}}i}}}} \end{array} \right. \end{eqnarray} $

(7)

根据实验数据，通过拟合液体密度、饱和蒸汽压的测量值(图 1，图 2)，利用二项式回归确定$ {\varOmega _{{\rm{a}}{T_{{\rm{r}}i}}}} $、$ {\varOmega _{{\rm{b}}{T_{{\rm{r}}i}}}} $。

根据对比温度的不同，可以得到$ {\varOmega _{{\rm{a}}{T_{{\rm{r}}i}}}} $的拟合关系式为

$ \begin{eqnarray} {\varOmega _{\rm{a}}}_{{T_{{\rm{r}}i}}}\! = \!\left\{ \begin{array}{l} 0.0305T_{{\rm{r}}i}^2\!-\!0.1017{T_{{\rm{r}}i}}\!+\!0.4923, \; \; {{T_{{\rm{r}}i}} < 1} \\ -0.0711T_{{\rm{r}}i}^2\!+\!0.0827{T_{{\rm{r}}i}}\!+\!0.3601, \; \; {{T_{{\rm{r}}i}} \geqslant 1} \end{array} \right. \end{eqnarray} $

(8)

根据对比温度的不同，可以得到$ {\varOmega _{{\rm{b}}{T_{{\rm{r}}i}}}} $的拟合关系式为

$ \begin{eqnarray} {\varOmega _{\rm{b}}}_{{T_{{\rm{r}}i}}}\! = \!\left\{ \begin{array}{l} 0.1012T_{{\rm{r}}i}^2 - 0.1479{T_{{\rm{r}}i}}\!+\!0.1321, \; \; {{T_{{\rm{r}}i}} < 1} \\[10pt] -0.0166T_{{\rm{r}}i}^2\!+\!0.0293{T_{{\rm{r}}i}}\!+\!0.0559, \; \; {{T_{{\rm{r}}i}} \geqslant 1} \end{array} \right. \end{eqnarray} $

(9)

2 状态方程在体系相平衡计算中的应用原理

设已知一个由$ n $个组分构成的烃类体系，取1 mol的物质量作为分析单元，那么当其在开发过程中处于任一气液相平衡状态时应有以下特征。

(1) 平衡气、液相的摩尔分数$ n_{\rm{V}} $和$ n_{\rm{L}} $分别在0$ \sim $1变化，且恒满足质量归一化条件

$ \begin{eqnarray} {n_{\rm{V}}} + {n_{\rm{L}}} = 1 \end{eqnarray} $

(10)

(2) 平衡气、液相的组成$ y_1 $，$ y_2 $，…，$ y_i $，…，$ y_n $及$ x_1 $，$ x_2 $，…，$ x_i $，…，$ x_n $应分别满足组成归一化条件

$ \begin{eqnarray} \left\{ \begin{array}{l} \sum {{x_i}} = 1\\[8pt] \sum {{y_i}} = 1\\[8pt] \sum {\left( {{y_i} - {x_i}} \right)} = 0 \end{array}\right. \end{eqnarray} $

(11)

(3) 平衡气、液相各组分的摩尔质量应满足物质平衡条件

$ \begin{eqnarray} {y_i}{n_{\rm{V}}} + {x_i}{n_{\rm{L}}} = {z_i} \end{eqnarray} $

(12)

(4) 任一组分在平衡气、液相中的分配比例可用平衡常数来描述，即

$ \begin{eqnarray} {K_i} = {{{y_i}} / {{x_i}}} \end{eqnarray} $

(13)

经数学处理，即可得到由平衡气、液相组成方程和物质平衡方程所构成的物质平衡方程

$ \begin{eqnarray} \sum {\left( {{y_i} - {x_i}} \right)} = {\mbox{$\sum$}} {\dfrac{{{z_i}\left( {{K_i} - 1} \right)}}{{1 + \left( {{K_i} - 1} \right){n_{\rm{V}}}}}} = 0 \end{eqnarray} $

(14)

当油气体系达到气液平衡时，体系中各组分在气液相中的逸度应相等^[24-25]，即

$ \begin{eqnarray} {y_i}\phi _i^{\rm{V}}p = f_i^{\rm{V}} = f_i^{\rm{L}} = {x_i}\phi _i^{\rm{L}}p \end{eqnarray} $

(15)

则，气液平衡常数可以表征为带有逸度的函数

$ \begin{eqnarray} {K_i} = \dfrac{{{y_i}}}{{{x_i}}} = \dfrac{{\phi _i^{\rm{L}}}}{{\phi _i^{\rm{V}}}} = \dfrac{{{{f_i^{\rm{L}}} / {{x_i}}}}}{{{{f_i^{\rm{V}}} / {{y_i}}}}} \end{eqnarray} $

(16)

根据热力学原理求解逸度的严格积分方程为

$ \begin{eqnarray} \left\{ \begin{array}{l} {\rm{R}}T\ln {\dfrac{{f_i^{\rm{V}}}}{{{y_i}p}}} = {\mbox{$\int$}}_V^\infty {\left[{{{\left( {\dfrac{{\partial p}}{{\partial {n_{{\rm{V}}i}}}}} \right)}_{{v_{\rm{V}}}, \; T, \; {n_{{\rm{V}}j}}}}- \dfrac{{{\rm{R}}T}}{{{v_{\rm{V}}}}}}\right]{\rm{d}}{v_{\rm{V}}}}-{\rm{R}}T\ln {z_{\rm{V}}}\\[15pt] {\rm{R}}T\ln{\dfrac{{f_i^{\rm{L}}}}{{{x_i}p}}} = {\mbox{$\int$}}_L^\infty {\left[{{{\left({\dfrac{{\partial p}}{{\partial {n_{{\rm{L}}i}}}}} \right)}_{{v_{\rm{L}}}, \; T, \; {n_{{\rm{L}}j}}}}- \dfrac{{{\rm{R}}T}}{{{v_{\rm{L}}}}}}\right]{\rm{d}}{v_{\rm{L}}}}-{\rm{R}}T\ln {z_{\rm{L}}} \end{array} \right. \end{eqnarray} $

(17)

相态计算的热力学平衡条件目标方程组为

$ \begin{eqnarray} \left\{ \begin{array}{l} {F_1}\left( {{x_i}, {y_i}, p, T} \right) = f_1^{\rm{L}} - f_1^{\rm{V}} = 0\\[9pt] {F_2}\left( {{x_i}, {y_i}, p, T} \right) = f_2^{\rm{L}} - f_2^{\rm{V}} = 0\\[9pt] \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \vdots \\[9pt] {F_i}\left( {{x_i}, {y_i}, p, T} \right) = f_i^{\rm{L}} - f_i^{\rm{V}} = 0\\[9pt] \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \vdots \\[9pt] {F_n}\left( {{x_i}, {y_i}, p, T} \right) = f_n^{\rm{L}} - f_n^{\rm{V}} = 0 \end{array} \right. \end{eqnarray} $

(18)

在实际应用过程中，由于研究对象是混合体系，需要引入混合法则对状态方程进行处理^[26-28]，相应的对常用混合法则修正，分析CO$ _2 $-原油的一般体系可知，各组分间极性和对称性存在差异，加入斥力作用系数和引力作用系数，可以更好地描述极性和对称性存在差异的组分体系，从而得到针对CO$ _2 $-原油体系的混合法则如下所示

$ \begin{eqnarray} \left\{ \begin{array}{l} {a_m} = \sum\limits_i {\sum\limits_j {{y_i}{y_i}{{\left( {{a_i}{a_j}} \right)}^{{1 / 2}}}\left( {1 - {k_{ij}}} \right)} } \\[13pt] {b_m} = \sum\limits_i {\sum\limits_j {{y_i}{y_j}\left[ {{{\left( {{b_i} + {b_j}} \right)} / 2}} \right]\left( {1 - {d_{ij}}} \right)} } \end{array} \right. \end{eqnarray} $

(19)

3 应用实例 3.1 密度计算

应用室内实验验证修正PR方程的计算结果。利用不同方法计算的密度见表 1，密度计算方法通过状态方程计算出的体积进行转换见式(20)，结果表明，采用RK和PR状态方程计算的密度偏低，存在较大误差，本研究建立的方程计算的密度相对误差为0.03%$ \sim $0.74%，与实验结果接近。表 2进一步列出了不同区块的油井油样在原始地层压力下密度的实验结果和计算结果的对比。

图 3$ \sim $图 6也给出了计算结果和误差的对比。可以看出，用本研究建立方程计算的密度相对误差均小于1%，表明该方程对CO$ _2 $-原油体系相态参数计算的适用性。

$ \begin{eqnarray} m = \rho v \end{eqnarray} $

(20)

3.2 最小混相压力(MMP)计算

根据实际油田油样组成(表 3)，用不同方法计算的MMP见表 4，结果表明，本研究建立方程计算的密度相对误差是2%$ \sim $8%，与实验结果接近。图 7$ \sim $图 8给出了计算结果和误差的对比。可以看出，用本研究建立方程计算的MMP误差均小于5%。

本文修正的状态方程不仅考虑了温度、压力和混合组分的影响，并且也修正了方程参数，使得非烃类组分也被考虑进来，因而可以获得比较满意的精度。利用本文修正后的状态方程计算所得到的CO$ _2 $驱最小混相压力与实验值比较接近，平均相对误差不超过5%；而常用估算最小混相压力经验关联式的计算结果与实际有较大的偏差，平均相对误差超过20%。同时，在计算量上，本方法虽然需要用到迭代求解，但是计算过程比较简洁，计算量不是很大，因而综合考虑精度和计算速度，本文所得到的方法比常用经验关联式具有更高的实际应用价值。

4 结论

(1) 常用的状态方程中，PR状态方程虽然相对误差20%左右，但由于形式简洁，有较好的应用效果。

(2) 通过对常用的PR状态方程引力项系数(温度系数)和斥力项系数(压力系数)的修正，可以进一步提高其精度，使其工程应用范围进一步扩展，可用于描述CO$ _2 $-原油体系等非理想性体系的相态变化特征。

(3) 修正后的状态方程对CO$ _2 $-原油体系的体积和密度参数计算相对误差1%、对MMP参数计算相对误差小于5%，在精度提高的基础上，计算速度和参数获取复杂程度与原方程简洁程度相似，具有较好应用价值。

符号说明

$ p $-压力，MPa；

$ v $-体积，m$ ^3 $；

R-普适气体常数，R = 8.314 J/(mol$ \cdot $K)；

$ T_{\rm{c}} $-临界温度，K；

$ a_i $- $ i $组分的引力修正系数，kPa$ \cdot $(dm$ ^3 $/mol)$ ^2 $；

$ b_i $- $ i $组分的体积修正系数，dm$ ^3 $；

$ \varOmega_{\rm{a}} $-引力修正系数中的无因次系数，无因次；

$ \varOmega_{\rm{b}} $-体积修正系数中的无因次系数，无因次；

$ T_{{\rm{c}}i} $- $ i $组分的临界温度，K；

$ p_{{\rm{c}}i} $- $ i $组分的临界压力，MPa；

$ \alpha(T) $-温度相关函数，无因次；

$ m_i $- $ i $组分的偏心因子相关系数，无因次；

$ T_{{\rm{r}}i} $- $ i $组分的对比温度，无因次；

$ \omega_i $- $ i $组分的偏心因子，无因次；

$ \varOmega_{{\rm{a}}T_{{\rm{r}}i}} $- $ \varOmega_{\rm{a}} $关联为对比温度函数后的参数形式，无因次；

$ \varOmega_{{\rm{b}}T_{{\rm{r}}i}} $- $ \varOmega_{\rm{b}} $关联为对比温度函数后的参数形式，无因次；

$ n_{\rm{V}} $-气相摩尔分数，无因次；

$ n_{\rm{L}} $-液相摩尔分数，无因次；

$ x_i $，$ y_i $，$ z_i $- $ i $组分液，气，总摩尔分数，无因次；

$ K_i $- $ i $组分气液平衡常数，无因次；

$ \phi _i^{\rm{V}} $- $ i $组分的气相逸度系数，无因次；

$ \phi _i^{\rm{L}} $- $ i $组分的液相逸度系数，无因次；

$ f_i^{\rm{V}} $- $ i $组分的气相逸度，无因次；

$ f_i^{\rm{L}} $- $ i $组分的液相逸度，无因次；

$ n_{{\rm{V}}i} $，$ n_{{\rm{V}}j} $- $ i $组分和$ j $组分气相摩尔分数，无因次；

$ n_{{\rm{L}}i} $，$ n_{{\rm{L}}j} $- $ i $组分和$ j $组分液相摩尔分数，无因次；

$ v_{\rm{V}} $，$ v_{\rm{L}} $-气相和液相体积，m$ ^3 $；

$ z_{\rm{V}} $，$ z_{\rm{L}} $-气相和液相偏差因子，无因次；

$ a_m $-体系混合引力修正系数，kPa$ \cdot $(dm$ ^3 $/mol)$ ^2 $；

$ b_m $-体系混合体积修正系数，dm$ ^3 $；

$ a_j $- $ j $组分的引力修正系数，kPa$ \cdot $(dm$ ^3 $/mol)$ ^2 $；

$ b_j $- $ j $组分的体积修正系数，dm$ ^3 $；

$ k_{ij} $，$ d_{ij} $- $ i $组分和$ j $组分的引力和斥力修正项，无因次；

$ m $-质量，$ \times $10$ ^3 $ kg；

$ \rho $-密度，$ \times $10$ ^3 $ kg/m$ ^3 $。