基于样条曲线插值天然气偏差因子计算方法

姬安召，王玉风，张光生

（陇东学院，甘肃庆阳 745100）

天然气偏差因子是在一定温度压力下天然气体积与该气体在相同条件下按理想气体定律计算的气体体积的比值。天然气的偏差因子是天然气性质研究过程中一个很重要的参数，例如在密度、黏度、等温压缩系数以及体积系数计算过程中必须考虑的参数之一。该参数主要决定于气体的组成、温度以及压力。天然气的偏差因子虽然可以通过实验进行测试，但实验的周期长，成本高，对实验仪器的精度要求高。但也可以通过相关的计算方法得到。目前，天然气偏差因子获取大致有三类：图版插值法、基于状态方程法和图版拟合法的数值计算方法。

在基于状态方程计算天然气偏差因子的研究方面，对状态方程的维里方程法中的AGA8 方法进行了研究，给出了该方法的适应性评价[1]。根据天然气气体组成[2]，研究了天然气压缩因子的计算，得出对低中压含碳气体，AGA8 方程计算误差最小，而含硫气体，Piper-Mahmoud 方法误差最小。根据实验研究结果[3]，对于超高压气藏而言，现有的图版法和经验公式对于天然气偏差因子计算的误差很大，甚至有些超压气藏天然气偏差因子误差超过20%。对图版拟合的经验公式和基于状态方程两参数法计算天然气偏差因子的公式进行评价，给出适合超高压气藏天然气偏差因子的计算模型。根据气体状态方程的BWRS 维里方程[4]，采用非线性回归分析，对该方程指数参数进行修正，提出了更高精度天然气偏差因子的计算方法。基于Strobridge 方法[5]，对BWR（Benedict-Webb-Rubin）状态方程进行了修正，修正后的方法计算天然气偏差因子的精度控制在1%以内。根据Nishiumi-Saito 状态方程[6]，采用多元回归方程，给出了天然气偏差因子计算的解析表达式。但这些方法都要基于不同的状态方程，还要进行复杂的数学计算，计算的结果在精度评定时，基本要参照Standing-Katz 图版的实验数据。

基于图版拟合方法，在天然气偏差因子的研究方面，对DAK、HTP、DRP、LXF、HY、Sarem、AGA8-92DC 七种方法计算天然气偏差因子的结果进行精度评估[7-10]，给出了每种方法的适用条件。对DAK 方法进行了修正[11]，修正后的方法提高了天然气偏差因子的计算精度。根据DPR 方法计算天然气偏差因子的方法[12]，采用一阶中心差分方法计算了超高压气藏天然气等温拟对比压缩系数，拓展了Trube 使用范围。利用低压传统天然气偏差因子图版数据和实验测试的高压天然气数据进行拟合[13]，得到了适合于中低压和高压的天然气偏差因子计算GXH 模型。这些图版拟合结果，基本都是基于Standing-Katz 图版数字化后的结果，或者有些参照了状态方程的基本形式，其不同的模型都用一定的适用条件。

除了基于状态方程法和图版拟合法之外，还有采用BP 神经网络进行预测方法，以测试压力、温度、天然气的组成以及其临界参数、测试的气体的偏差因子作为输入[14]，采用BP 神经网络方法对249 个实验数据点进行了样本学习，建立基于BP 神经网络的天然气偏差因子的预测模型，但该模型对输入数据的依赖性强，实用性不强。基于气藏物质平衡原理进行天然气偏差因子计算的模型[15]，根据生产动态资料，结合干气气藏的物质平衡方程，提出了基于生产资料计算天然气偏差因子的新方法，但计算的精度依赖于气藏特性和生产动态资料的准确性。在商业软件进行天然气偏差因子的评价方面，对目前流行的四种商业软件[16]（PIPESIM 2012、PROPSER 13.0、PVTp 10.0、Multiflash 4.1）对天然气偏差因子计算模块进行了评价，其中PROPSER 13.0 计算精度最好，对比不同的状态方程计算的结果，SRK 方程计算结果最理想。

天然气中一般还有酸性气体（CO2、H2S），对于还有酸性气体的天然气计算的偏差因子需要做矫正。应用酸性气体偏差因子的计算方法[17]，根据Keys、SBV 以及RK 混合规则，酸性气体的临界参数校正方法，优选出了CO2气藏的偏差因子计算的高精度RK+WA 模型。研究富含酸性气CO2天然气偏差因子的计算方法[18-20]，通过实验结果与理论模型的对比，优选出来基于GXO校正的DAK 模型计算误差最小，并给出了富含CO2天然气偏差因子随温度与压力变化规律以及CO2含量对偏差因子计算结果的影响。通过对比目前天然气中含有酸性气体（CO2、H2S）时偏差因子的矫正方法[21]，得出DPR 模型和DAK 模型对含有酸性气体的矫正后天然气偏差因子计算结果最好。其他模型（GXQ、CKB）结果不理想。

根据前人的研究结果，以Poettman、Katz 和Smith对Standing-Katz 图版低压范围（0.2 ≤Ppr≤15.0、1.05≤Tpr≤3.00）数值化的5 940 个数据点[22，23]以及Poettman 对高压范围（15.0≤Ppr≤30.0、1.4≤Tpr≤2.8）数值化Katz 图版的1 208 个数据点为基础，采用样条函数差值的原理，对同一视对比温度Tpr下的天然气视对比压力与偏差因子进行样条函数差值，既保证了实验数据点不偏离差值函数，又保证了差值曲线的连续光滑性。若视对比温度Tpr不在实验数据给定的曲线上，则根据实验数据给定的视对比温度Tpr值在Standing-Katz 图版上选择距离被差值的视对比温度Tpr最近的两条曲线，进行样条曲线间的二次差值，则可以得到满足Standing-Katz 图版的任意视对比温度Tpr和视对比压力Ppr下的天然气偏差因子的值。如果天然气中含有酸性气体，则进行矫正。

1 基本原理

1.1 曲线差值方法

根据Standing-Katz 图版数值化的数据点，选取任意一条以视对比温度Tpr为参数数值化离散数据点，将视对比压力Ppr作为自变量，Standing-Katz 图版上偏差因子Z 为因变量。在每个相邻的视对比压力之间的子区间进行三次多项式差值，并且要求在整个区域子区间上的差值多项式函数具有二阶连续导数，这样得到差值函数即为三次样条差值函数。这样既保证了区间内端点的连续性，也保证了内端点处的光滑性，并且被差值的数据点落在差值曲线上。

假设在［Pprmin，Pprmax］区间上有n+1 个离散数据点Ppr（0）＜Ppr（1）＜…＜Ppr（n-1）＜Ppr（n），Pprmin为差值区间视对比压力的最小值，Pprmax为差值区间视对比压力的最大值。天然气偏差因子Zi（Ppr）在区间［Ppr（i-1），Ppr（i）］上是三次多项式函数，这个三次多项式函数可表示为，即有4 个待定系数。因为n+1 个离散数据点有n 个子区间，每个子区间要确定一个三次多项式，因此共计有4n 个待定系数需要计算。根据差值的基本条件，即共计n+1 个条件。在区间内端点处，Zi（Ppr）、Zi（Ppr）的一阶导函数以及二阶导函数为连续函数，即可得3n-3个条件。还需要两个条件，即边界条件，这里取差值区间的两个端点处的一阶导数即可，其值由区间起点与终点的相邻两点的直线的斜率确定。根据上述分析，在区间x∈［Ppr（i-1），Ppr（i）］上的三次样条差值多项式的数学模型可由（1）式表示：

其中：i=1，2，3，…，n。通过（2）式，可以计算位于差值区间任意视对比压力下的天然气偏差因子。

1.2 曲线间差值方法

2 计算结果

根据（2）式的分析，在进行曲线差值之前，需要计算Standing-Katz 图版数值化后的每条曲线两个端点的一阶导函数。这里分别取曲线每条数值化两个相邻端点的点进行一阶导数的计算。为了分析一阶导数对差值结果的影响以及偏差因子的计算精度，这里分别选取两个方案，方案一：以曲线端点相邻两点之间的直线的斜率代替一阶导数。方案二：以曲线端点相邻三点进行一阶导数加权计算，其中方案二的一阶导数计算公式由（4）式表示。

其中：当k=0 时，则j=1；当k=n 时，则j=n-1。

根据Standing-Katz 图版数值化的结果，结合上述一阶导数的计算方法，曲线的两个端点的一阶导数的结果（见表1）。从表1 可以看出，方案一与方案二计算的结果差值曲线两个端点的一阶导数在左端点和右端点的相对误差很小，基本控制在5%以内，其中视对比温度为1.8 和1.9 时，方案一与方案二计算的左端点的一阶导数值误差略大于5%。

根据表1 中方案一和方案二一阶导数的计算结果，这里取二者的平均值，将二者平均值代入公式（1）求解差值内节点的一阶导数值，将内节点的一阶导数值代入到（2）式，得到对应视对比温度下的样条函数；然后在区间［0.2，30］上离散视对比压力，这里视对比压力离散的间距为0.001，然后将离散的视对比压力代入到（2）式，得到样条函数差值后结果（见图1），图1中的点表示Standing-Katz 图版的数值化的结果。这里需要说明的是Standing-Katz 图版数值化的数据太多，在绘图时按照一定间隔进行了抽稀处理。通过图1 的对比，Standing-Katz 图版的数值化的结果与样条差值的计算结果匹配很好，并且样条差值法可以计算满足Standing-Katz 图版差值范围的任意给定视对比温度与视对比压力下的天然气偏差因子的值。

图1 Standing-Katz 图版与样条曲线拟合结果对比图

表1 方案一与方案二计算拟合曲线的左右端点的一阶导数值与相对误差

根据曲线间的差值原理，这里选取了视对比温度分别为1.08、1.55 和2.10 的三组参数。将这三组参数与Standing-Katz 图版的视对比温度对比可知，这三组参数均不是Standing-Katz 图版所给定的参数，因此这里只能通过曲线间差值的方法得到不同视对比压力下的偏差因子值。这三组视对比温度参数分别介于视对比温度为［1.05，1.10］、［1.50，1.60］和［2.00，2.20］三个区间，根据（3）式的计算方法，得到这三组参数差值的结果（见图2），Standing-Katz 图版数值化的数据见图2的离散点所示。从图2 可以看出，差值的实线的结果与Standing-Katz 图版数值化的变化趋势一致。

图2 Standing-Katz 图版与样条曲线拟合结果对比图

3 实例分析

在垂直管流计算过程中，天然气偏差系数的计算模型有很多种，根据前人的研究结果，本次研究主要采用了四种计算方法，即DAK（Dranchk-Abu-Kasse）模型，HTP（Hankinson-Thomas-Phillips）模型，李相方（LXF）模型和本文提出的样条函数差值法。

根据文献提出的气井井底静压的计算原理[25]，采用微元分析法的思想，将气井的井筒从井口到井底离散为多个单元，每一个计算微元体内的性质是一致的，其属性值可以看作常数处理。基于以上考虑，可以把井筒从井口到产层中深划分为n 个微元体，在每个微元体内，天然气的温度、偏差因子、黏度、压力以及与压力有关的量可以看作常数（都是基于前一个节点的压力与温度计算出的值）。

由于969 井区的天然气中含有酸性气体，考虑到酸性气体对天然气偏差因子的影响，采用了GXQ 方法和Wichert-Aziz 方法对天然气的临界温度和压力分别做了校正，具体方法参见文献[19-21]。根据试气资料的整理，统计了969 井区共计76 口井生产前的压力恢复数据，这些井也下入了井下压力计，测量了实际的静压。结合垂直管流微元分析法的计算思想，根据关井恢复后的井口油压和套压的测量值计算了这76 口井的静压。将静压计算的结果与实测的进行对比，得到了969井区分别采用Dranchk-Abu-Kasse 模型[6]、Hankinson-Thomas-Phillips 模型[8]、李相方（LXF）模型[7]以及样条函数差值模型计算的天然气偏差因子，然后根据微元分析进行井底静压计算，井底静压的误差分布（见图3）。从图3 可以看出，本文给出的样条函数差值法计算的偏差因子在计算井底静压时的相对误差最小；其次为LXF 法，因为该方法在计算偏差因子时，采用了分段函数进行处理，适用性较好；然后为DAK 和HTP法。其中井底静压相对误差大于5%的这些井，基本为气水同产，甚至有些井存在着井底积液，导致井底静压计算的相对误差较大。

图3 969 井区井底静压误差统计图

4 结论

（1）根据Standing-Katz 图版数值化的结果，采用三次样条函数的差值原理，给出了三次样条差值计算天然气偏差因子的基本原理。

（2）对于视对比温度Tpr不在Standing-Katz 图版上的情况，在Standing-Katz 图版上查找距离给定视对比温度Tpr最近的两条曲线上的视对比温度，然后根据视对比温度Tpr的值进行反距离加权二次差值。根据反距离加权二次差值，可以得到任意给定视对比温度Tpr（1.05≤Tpr≤3.00）和视对比压力Ppr（0.2≤Ppr≤15.0）下的天然气偏差因子的值。

（3）通过969 井区76 口气井的试气资料，本文提出的三次样条函数差值计算偏差因子方法具有较高的精度，并且计算简单，适用性强。