天然气偏差因子(Z)是油气藏工程、天然气化工等领域常用的重要参数，它表征了某一温度、压力条件下相同数量真实气体与理想气体体积的比值，通常又称压缩因子或偏差系数。确定偏差因子的方法主要有3类，即：实验测量法、图版法(查表法)和计算法。实验测量法虽然直接可靠，但耗时、成本高；图版法(查表法)难以得到连续的值，因此计算法以其简便性和实用性被广泛应用。
计算法主要是利用能够代表偏差因子标准图版的关系式进行编程计算，其中显式Z因子关系式多是通过回归分析总结得到，如Beggs、Gopal、Kumar、Azizi、Heidaryan和Salarabadi、Heidaryan和Moghadasi等提出的表达式[1-8]；而隐式关系式则一般源于状态方程。状态方程有Van der Waals方程、RK方程、SRK方程、PR方程等立方型状态方程[9-19]，以及Kamerlingh Onnes方程、Beattie方程、BWR方程、Strobridge方程、Carnahan方程、Morsy方程、BWRS方程、LeE-Kesler方程等维里状态方程[20-35]。维里状态方程在偏差因子计算方法的发展过程中具有重要意义，例如：Dranchuk等利用BWR方程拟合1500组偏差因子数据提出了8系数偏差因子表达式即DPR方法[36-37]; Hall和Yarborough在Carnahan-Starling硬球方程的基础，提出了一种Z因子关系式即HY方法[27-28]；Dranchuk和Abou-Kassem则应用BWRS方程提出了一种11系数偏差因子计算式即DAK方法[38]。这3种方法是应用状态方程计算天然气偏差因子的经典方法[39-44]，在一定范围内它们均能比较准确地代表Standing-Katz图版[37]，其中以DAK方法的计算精度最高[45-48]，但其在高温高压下的误差也稍大。
为更加准确地预测整个压力范围(0.2≤ppr≤30.0)内的天然气偏差因子，本文基于Nishiumi-Saito方程[49]，提出一种新的偏差因子计算方法，以期为油气藏工程相关研究者提供参考。
利用状态方程拟合偏差因子数据进而总结经验关系式的方法是计算天然气偏差因子的有效工具。Benedict等提出的BWR方程颇具代表性[22-24]，因其能比较准确地计算气体的热力学性质而迅速得到应用。Dranchuk基于BWR方程提出了DPR方法[36], 随后许多学者对其进行了修正和推广，以便更加准确、有效地预测更大温度与压力范围内的纯物质和混合物的热力学性质参数。其中，Starling修正的BWR方程(即BWRS方程)应用更为广泛[30-34]，DAK方法便是基于此。
然而，Nishiumi和Saito指出BWRS方程仍不能有效地预测低对比温度下的热力学参数[49]，他们提出了另一种广义的BWR方程，称之为Nishiumi-Saito方程，其表达式如式(1)：
式中：p为压力，Pa；R为摩尔气体常数[50-51]，8.314 459 8 J/(mol·K)；T为温度，K；ρ为密度，kg/m3；M为摩尔质量，kg/mol；B0、A0、C0、D0、E0、b、a、d、e、f、α、c、g、h、γ为与状态方程相关的系数。
真实气体状态方程为：
式中：Z为气体偏差因子。
在式(1)两端同时乘以M/(ρRT)，得到：
引入(拟)对比密度ρpr、(拟)对比温度Tpr和(拟)对比压力ppr：
式中：Tpr为(拟)对比温度；Tpc为(拟)临界温度, K；ppr为(拟)对比压力；ppc为(拟)临界压力, Pa；ρpr为(拟)对比密度; ρpc为(拟)临界密度, kg/m3；Zc为临界偏差因子，计算时一般取0.27。
将式(4)和式(5)代入式(3)，整理得到：
式(6)即为基于Nishiumi-Saito状态方程导出的新的隐式关系式。式(8)中A1~A15为系数。
Poettman率先发表了从Standing-Katz图版(见图 1)上读取的偏差因子数据表(0.2≤ppr≤15.0 & 1.05≤Tpr≤3.00)[37, 52]，Katz和Smith沿用了Poettman的数值化成果并对个别数据点进行了更正[53-54]。
然而，发现Poettman数据表(共297×20=5940个数据点)中的5个数据点疑似有误，因为在这些点处，Z值发生了突变。参考这些点附近压力、温度范围内Z的取值，相应作了更正(见表 1)。表 1中，Z1、Z2、Z3所在列分别为Poettman、Katz、Smith数据的偏差因子取值，Zsta所在列为本文更正的偏差因子值。
对于常用的压力范围(即中低压力范围0.2≤ppr≤15.0 & 1.05≤Tpr≤3.00)，以更正后的Poettman数值化的偏差因子数据作为标准。类似地，对于高压范围15.0≤ppr≤30.0 & 1.40≤Tpr≤2.80)内的Katz图版(见图 2)[53]，同样利用数值化处理结果作为标准数据(此处每条等温线上ppr每隔0.1取一个点，共151×8=1208个数据点)。
利用上述常用压力范围和高压范围的7148个偏差因子数据点，对式(6)作多元非线性回归分析，得到各系数A1~A15的取值见表 2。
式(6)和式(7)反映了偏差因子和对比密度ρpr的隐式关系。采用牛顿迭代法可快速解得Z值[55]。首先构造如下函数f(ρpr)：
求f(ρpr)关于ρpr的一阶导数：
建立牛顿迭代公式[55]：
式中：(ρpr)0为上次计算的对比密度；(ρpr)1为本次得到的对比密度; B1~B6为系数。
结合式(9)~式(12)，设置如下迭代步骤求解偏差因子：
(1) 给定ppr、Tpr，设定计算精度eps和最大迭代次数M。
(2) 假定偏差因子Z0=1, 令(ρpr)1=B6/Z0，迭代次数m=0。
(3) 将(ρpr)1赋给(ρpr)0。
(4) 利用式(12)求得对比密度(ρpr)1；将m+1赋值给m。
(5) 比较(ρpr)1和(ρpr)0的差距，若到达精度要求或迭代次数超过限制(即abs[(ρpr)1-(ρpr)0] < eps或m>M)则终止迭代；否则，继续步骤(3)~(5)。
(6) 计算偏差因子，Z=B6/(ρpr)1。
利用Standing-Katz图版数值化的5940组数据和Katz图版数值化的1208组数据[37, 52-54]，对DPR、HY、DAK和本方法进行误差评价，误差计算式如式(13)：
式中：AAE为平均绝对误差，%；Zcal为计算的Z值；Zsta为Z的图版数值化值。
中低压范围(0.2≤ppr≤15.0 & 1.05≤Tpr≤3.00)内的计算误差记为AAE1，相对高压范围(15.0≤ppr≤30.0 & 1.4≤Tpr≤2.8)内的误差记为AAE2，总共7148个数据点的平均绝对误差记为AAE3。
表 3列出了4种偏差因子计算方法的平均绝对误差；表 4~表 6列出了“0.2≤ppr≤15.0 & 1.05≤Tpr≤3.00”以及“15.0≤ppr≤30.0 & 1.4≤Tpr≤2.8”范围内各等温线的计算误差。由表 3~表 6可看出，本方法的AAE1为0.357%，AAE2为0.066%，其计算精度高于其他3种方法，对比低温(尤其是Tpr=1.1)和高压(ppr>15.0)条件下的计算结果，优势更为明显。
图 3~图 10统计了DPR、HY、DAK和本方法在中低压和高压范围内的偏差因子计算误差，图中颜色从蓝到红体现误差由小变大。在中低压(0.2≤ppr≤15.0 & 1.05≤Tpr≤3.00)范围内，4种方法在Tpr=1.05等温线上均存在误差较大的点，图 3~图 6中缺失部分表示误差大于10%的区域。在高压(15.0≤ppr≤30.0 & 1.4≤Tpr≤2.8)范围内，图 7和图 9中缺失部分表示误差大于2%的区域。
对于中低压区，4种方法中HY方法的缺失区域最大，在“1.30≤ppr≤2.75 & Tpr=1.05”和“1.50≤ppr≤1.65 & Tpr=1.1”范围内的共计34个点的误差大于10%(DPR有31个点，DAK有29个点)；DPR方法的平均绝对误差最大；而本方法的计算精度最高，且误差大于10%的区域最小(只有“1.00≤ppr≤1.45 & Tpr=1.05”范围内的10个点)，整体误差也最小。
对于高压区，DPR和DAK方法在部分高温区域的计算误差大于2%，DAK方法的平均绝对误差最大(AAE2为0.979%)，HY方法的误差相对较小但误差分布不平衡，而本方法的误差分布十分平滑且AAE2仅为0.066%。
(1) 基于Nishiumi修正的BWR状态方程(Nishiumi-Saito方程)导出了15系数偏差因子关系式，进而提出了一种确定气体偏差因子的新方法。
(2) 对比油气藏工程常用的DPR、HY、DAK方法，该方法的计算效果更好。除了Tpr=1.05等温线上1.00≤ppr≤1.45范围内的10个点的计算误差较大外，该方法既适用于常用压力范围(0.2≤ppr≤15.0 & 1.05≤Tpr≤3.00)，亦可应用于高压范围(15.0≤ppr≤30.0 & 1.4≤Tpr≤2.8)，在这两个范围内的平均绝对误差分别为0.357%、0.066%，优于前3种常用方法。
2019, Vol. 48 Issue (1): 91-98
2019, Vol. 48 Issue (1): 91-98