基于BWRS状态方程的天然气偏差因子计算方法

张立侠 郭春秋

中国石油勘探开发研究院

天然气偏差因子（Z）通常又称压缩因子或偏差系数，它表征了某一压力、温度条件下相同数量真实气体和理想气体的体积之比，是油气藏工程计算常用的重要参数；为避免与气体压缩系数（Cg）混淆，后文均称之为偏差因子。确定偏差因子的方法主要有实验测量、图版法和计算法3类。实验法虽然直接可靠，但耗时、成本高，而图版法或查表则不能得到连续的值，因此，计算法以其简便性和实用性得到了广泛应用。计算法主要是利用能够代表偏差因子标准数据的经验关系进行编程计算，这些偏差因子关系式多是源于状态方程，或直接通过回归分析总结得到。

利用状态方程拟合偏差因子数据进而总结经验公式的方法是计算天然气偏差因子的有效工具。气体状态方程有许多种，如Van der Waals方程［1-3］、RK方程［4-6］、SRK方程［7-9］、PR方程［10-12］等立方型状态方程，以及Starling方程［13-14］、BWR方程［15-16］、BWRS方程［17-21］等维里状态方程。维里状态方程在气体混合物偏差因子计算方法的发展过程中起到了重要的推动作用，较有代表性的有：Hall（1973）［13］在Carnahan-Starling（1969）硬球方程［14］的基础上提出了一种迭代计算偏差因子的关系式，即HY方法；Dranchuk等人（1973）［22］利用BWR状态方程拟合偏差因子数据［23］提出了一种8系数表达式，即DPR方法；Dranchuk和Abou-Kassem（1975）［24］则应用BWRS方程提出了一种11系数偏差因子计算式，即DAK方法。这3种方法是应用状态方程计算天然气偏差因子的经典方法，在一定范围内它们均能比较准确地代表Standing-Katz图版数据［23］，其中DAK方法的计算精度最高，但其在高温高压下的计算误差也稍大，存在改进的空间。为此，笔者基于改进的BWRS状态方程提出一种新的偏差因子计算方法，以期更精确地确定天然气偏差因子。

1 改进的BWRS状态方程

Benedict等人（1940，1942）［15-16］提出了一种剩余功（residual work content）与密度、绝对温度之间的经验方程，给出了与之对应的状态方程，即Benedict-Webb-Rubin方程，其表达式为方程（1）。

BWR状态方程因其能比较准确地计算气体的热力学性质而得到了广泛应用，例如，Dranchuk等人（1973）［22］利用BWR状态方程拟合平滑处理后的Standing-Katz图版［23］上的1 500个数据点，提出了DPR方法。后来，许多学者对BWR方程进行了修正和推广，以期更加准确、有效地预测更大压力和温度范围内的纯物质及混合物的热力学性质参数。其中，Starling（1970）［17-18］修正的BWR方程（即BWRS方程）颇具代表性，应用较为广泛［19-21］，其表达式为方程（2）。

式中，p为压力，Pa；R为摩尔气体常数，8.3144598 J/（mol·K）；T为温度，K；ρ为密度，kg/m3；M为摩尔质量，kg/mol；B0、A0、C0、b、a、α、c、γ均为与状态方程相关的系数；D0、E0和d为BWRS方程相较于BWR方程增加的3个系数。

BWRS状态方程是基于统计力学摄动理论对BWR方程中的系数C0和a相关项作了修改而得到的。增加的3个系数D0、E0和d有助于拟合并预测流体在低温及临界区的热力学参数（焓、熵、蒸气压等）。利用BWRS方程，Dranchuk（1975）［24］同样应用DPR关系式所用的1 500个偏差因子数据点进行拟合，提出了DAK方法。由于状态方程的改进，DAK方法比DPR（1973）方法具有更高的计算精度，但其在15≤ppr≤30 & 1.4≤Tpr≤2.8高压范围内的计算结果与Katz（1959）高压偏差因子图版之间的误差稍大，仍可改进。胡建国（2013）［25］利用6 548组偏差因子数据对DAK方法的11个系数进行了修正，改善了高压下的偏差因子计算精度，但这种修正在0.2≤ppr≤15 & 1.05≤Tpr≤3.0范围内的计算结果却与Standing-Katz图版数据不符。

研究发现，对式（2）中的系数γ进行调整可明显提高基于BWRS状态方程的偏差因子计算方法的计算精度。为此，将式（2）改写为

由真实气体状态方程得

在式（3）两端同时乘以M/(ρRT)，得

引入对比密度ρpr、对比温度Tpr和对比压力ppr

则式（5）可化为

其中

式中，γ1、γ2为改进的BWRS方程的系数；Z为气体偏差因子；Tpr为（拟）对比温度；ppr为（拟）对比压力；Tpc为（拟）临界温度；ppc为（拟）临界压力；ρpr为（拟）对比密度；ρpc为（拟）临界密度；Zc为临界偏差因子，计算时一般取0.27。

式（7）和式（8）反映了偏差因子Z与Tpr、ppr和ρpr的隐式关系，相较于DAK方法，该关系式增加了一个系数，同样可利用迭代法求解Z。

2 偏差因子计算方法

使用式（7）计算Z需要先确定A1～A12等12个系数的值。Poettmann（1952）［26］、Katz（1959）［27］及Smith（1990）［28］均发表了0.2≤ppr≤15 & 1.05≤Tpr≤3.0范围内Standing-Katz图版（图1）数值化的标准偏差因子数据表（共297×20=5 940个数据点）；Poettmann［26］数据表中的5个数据点疑似有误（因为在这些点处，Z值发生了突变），作如表1修改。

图1 Standing-Katz 偏差因子图版［23］Fig.1 Standing-Katz Z-factor chart［23］

表1 5个偏差因子数据点的修改值Table 1 Modified values of 5 Z-factor data points

对于常用的压力范围，以Poettmann（1952）数值化的偏差因子数据作为标准；对于“15≤ppr≤30、1.4≤Tpr≤2.8”范围内的Katz（1959）［27］图版（图2），同样利用数值化处理结果作为参照依据（每条等温线上对比压力ppr每隔0.1取一个点，共151×8=1 208个数据点）。

图2 15 ≤ ppr ≤ 30 范围内的Katz 图版［27］Fig.2 Katz chart for 15 ≤ ppr ≤ 30［27］

为方便数据拟合，式（7）可写为：

函数g（Tpr，ρpr，Z）以Tpr、ρpr、Z为自变量，且其函数值应恒为0；利用上述7 148组数据点分2个区拟合该函数，求满足式（10）时各系数A1～A12的值，此过程可用Matlab中的非线性最小二乘回归（nonlinear least squares regression）函数编程实现。各系数A1～A12的取值见表2。

表2 改进BWRS状态方程的系数Table 2 Coefficients of modified BWRS Equation

利用牛顿迭代法［29］求解偏差因子Z较为方便，可令函数f(ρpr)为

其中

f(ρpr)关于ρpr的导数为

牛顿迭代公式为

式中，(ρpr)0为上次计算的对比密度；(ρpr)1为本次得到的对比密度。

具体求解步骤如下：

（1）给定ppr、Tpr，设定最大迭代次数N和计算精度eps；

（2）令(ρpr)1=0.27ppr/Tpr，k=0；

（3）将(ρpr)1赋给(ρpr)0；

（4）利用式（13）得到对比密度(ρpr)1；k+1→k；

（5）比较(ρpr)1和(ρpr)0的绝对差，若k＞N或abs［(ρpr)1-(ρpr)0］＜eps，则迭代终止；否则，继续步骤（3）～（5）。

（6）计算偏差因子，Z=0.27ppr/［Tpr（ρpr)1］。

3 计算结果对比

利用Standing-Katz图版数值化的5 940组数据［26-28］和Katz（1959）图版数值化的1 208组数据对表2中3种取值方法进行评价（平均绝对误差记为AAE），0.2≤ppr≤15 & 1.05≤Tpr≤3.0范围（中低压区）内的计算误差记为AAE1，15≤ppr≤30& 1.4≤Tpr≤2.8范围（高压区）内的计算误差记为AAE2，总共7 148个数据点的计算误差记为AAE3。

式中，Zcal为计算的Z值。

表3列出了这3种相关的偏差因子计算方法的平均绝对误差；本文提出的12系数关系式的AAE1为0.382%，相较于原DAK方法，其计算精度提高了12.084%；而AAE2为0.205%，比DAK方法的计算精度高了79.041%；0.2≤ppr≤15 & 1.05≤Tpr≤3.0以及15≤ppr≤30 & 1.4≤Tpr≤2.8范围内各等温线的计算误差见表4、5。

表3 误差分析Table 3 Error analysis

表4 0.2≤ppr≤15 & 1.05≤Tpr≤3.0范围内各等温线的平均绝对误差 %Table 4 Average absolute error of each isothermal line for 0.2≤ppr≤15 & 1.05≤Tpr≤3.0

表5 15≤ppr≤30 & 1.4≤Tpr≤2.8范围内各等温线的平均绝对误差 %Table 5 Average absolute error of each isothermal line for 15≤ppr≤30 & 1.4≤Tpr≤2.8

图3、图4统计了DAK（1975）［24］、DAK-HJG（2013）［25］和本文提出的新方法在中低压和高压范围内的偏差因子计算误差，其大小由颜色从蓝到红来体现。对于“0.2≤ppr≤15 & 1.05≤Tpr≤3.0”的中低压情形，计算误差以10%为界，图3中缺失部分表示误差大于10%的区域。同理，对于“15≤ppr≤30& 1.4≤Tpr≤2.8”的高压情形，由于计算误差整体较小，图4只展示误差小于2%的区域。

图3 0.2≤ppr≤15 & 1.05≤Tpr≤3.0范围内3种方法的误差分布Fig.3 Error distribution for 0.2≤ppr≤15 & 1.05≤Tpr≤3.0 by three methods

图4 15≤ppr≤30 & 1.4≤Tpr≤2.8范围内3种方法的误差分布Fig.4 Error distribution for 15≤ppr≤30 & 1.4≤Tpr≤2.8 by three methods

对于中低压区，DAK-HJG（2013）方法的缺失区域较大，蓝色区域占比较低，其AAE1=12.733%，由此可见，该方法在0.2≤ppr≤15、1.05≤Tpr≤3.0范围内并不适用；而DAK（1975）和新方法的计算精度较高（AAE1分别为0.435%、0.382%），但新方法误差大于10%的区域更小，整体误差也更小。

对于高压区，DAK（1975）方法在高温区域部分点的计算误差大于2%（部分区域缺失），其平均绝对误差也最大（AAE2为0.979%）；而DAK-HJG方法和新方法的误差较小（AAE2分别为0.219%和0.205%），但新方法的最大误差更小（仅为1.091%）。

笔者认为BWRS方程中指数项中的γ与前一项的γ取值不相同时该方程能更准确地代表气体（尤其是在高压下）的性质；由于增加了一个系数，修正后的方程对于天然气偏差因子数据的代表能力应更强。此外，本文使用了0.2≤ppr≤15 & 1.05≤Tpr≤3.0范围内的5 940个数据点（对Poettmann发布的偏差因子标准数据［26］进行了更正）以及15≤ppr≤30 & 1.4≤Tpr≤2.8范围内的1 208个数据点，偏差因子Z数据体更大，故新方法适用范围更广，计算精度更高。

4 结论

（1）基于改进的BWRS状态方程和DAK方法导出了一种12系数偏差因子关系式，利用7 148组偏差因子数据结合非线性回归分析确定了这些系数，提出了一种新的Z值确定方法。

（2）由于BWRS方程中指数项的改进，该方程能更准确地代表及预测气体（尤其是在高压下）的物理性质。

（3）由于偏差因子数据的合理性以及状态方程的改进，新方法的计算精度更高。新方法在中低压和高压下的平均绝对误差分别为0.382%、0.205%，这比原DAK方法的计算精度分别高出12%和79%；推荐使用范围为：0.2≤ppr≤15 & 1.05＜Tpr≤3.0以及15≤ppr≤30 & 1.4≤Tpr≤2.8。而胡建国修正的DAK方法只适用于高压（15≤ppr≤30）的情形。