AGA8—92DC计算方法天然气压缩因子计算(最漂亮的)

AGA8—92DC计算方法天然气压缩因子计算摘要：按照GB／T 17747．2—1999《天然气压缩因子的计算第2部分：用摩尔组成进行计算》，采用AGA8—92DC计算方法，用VB编程计算了天然气压缩因子。

用二分法求解状态方程，精度满足工程需要。

关键词：压缩因子；AGA8—92DC计算方法；二分法1概述工作状态下的压缩因子是天然气最重要的物性参数之一，涉及到天然气的勘探、开发、输送、计量和利用等各个方面。

实测天然气压缩因子所需的仪器设备价格高，不易推广，因此计算方法发展很快，主要为经验公式和状态方程计算方法。

1992年6月26日，国际标准化组织(ISO)天然气技术委员会(TC193)及分析技术分委员会(TC193/SC1)在挪威斯泰万格(Stavanger)召开了第四次全体会议，会上推荐了两个精度较高的计算工作状态下天然气压缩因子的方程，目PAGA8-92DC方程、SGERG-88方程[1]。

随后，国际标准化组织于1994年形成了国际标准草案[2]。

AGA8-92DC方程来自美国煤气协会(AGA)。

美国煤气协会在天然气压缩因子和超压缩因子表的基础上，开展了大量研究，于1992年发表了以状态方程为基础计算压缩因子的AGA No.8报告及AGA8-92DC方程[2]。

1994年，四川石油管理局天然气研究所遵照中国石油天然气总公司技术监督局的指示，对国际标准化组织1992年挪威斯泰万格会议推荐的AGA8-92DC 方程、SGERG-88方程进行验证研究，于1996年底基本完成[2]。

1999年，四川石油管理局天然气研究院(前身为四川石油管理局天然气研究所)起草的《天然气压缩因子的计算》GB/T 17747.1～3—1999被批准、发布。

《天然气压缩因子的计算》GB/T 17747.1～3—1999包括3个部分：《天然气压缩因子的计算第1部分：导论和指南》GB/T 17747.1—1999，《天然气压缩因子的计算第2部分：用摩尔组成进行计算》GB/T 17747.2—1999，《天然气压缩因子的计算第3部分：用物性值进行计算》GB/T 17747.3—1999。

GB/T 17747.1等效采用ISO 12213—1：1997《天然气压缩因子的计算导论和指南》。

GB/T 17747.2等效采用ISO 12213-2：1997《天然气压缩因子的计算用摩尔组成进行计算》，给出了用已知的气体的详细的摩尔组成计算压缩因子的方法，目PAGA8—92DC计算方法。

GB/T 17747.3等效采用ISO 12213-3：1997《天然气压缩因子的计算用物性值进行计算》，给出了用包括可获得的高位发热量(体积基)、相对密度、C02含量和H2含量(若不为零)等非详细的分析数据计算压缩因子的方法，即SGERG-88计算方法。

笔者在输气管道和城镇高压燃气管道水力计算中，按照GB/T 17747.2采用AGA8-92DC计算方法进行天然气压缩因子计算，效果良好。

本文对其中的一些问题进行探讨，受篇幅所限，一些内容文中适当省略，详见GB/T 17747.2。

2AGA8—92DC方法的计算过程2.1已知条件、待求量、计算步骤2.1.1已知条件按照GB／T 17747．2的要求，以CH4、N2、CO2、C2H6、C3H8、H2O、H2S、H2、CO、O2、i−C4H10、n−C4H10、i−C5H12、n−C5H12、n−C6H14、n−C7H16、n−C8H18、n−C9H20、n−C10H22、He、Ar共21种组分的摩尔分数表示气体的组成。

如果C7H16、C8H18、C9H20、C10H22摩尔分数未知，允许用C6+表示总的摩尔分数。

应进行敏感度分析，以检验此近似法是否会使计算结果变差。

在输入摩尔组成时，将各组分按上述顺序排列(CH4、N2、...)，输入摩尔分数值，合值为 1 。

若不存在某组分，则其摩尔分数值为0 。

以最后一个摩尔分数不为0的组分来计算组分数。

例如，某天然气的已知摩尔组成见表1，则输入数据见表2，组分数为15。

表1 某天然气的已知摩尔组若已知体积分数组成，则将其换算成摩尔分数组成，具体换算方法见GB ／T 11062--1998《天然气发热量、密度、相对密度和沃泊指数的计算方法》。

标准名称中的“沃泊指数”现称为“华白数”。

需要注意的是，GB／T 11062中体积分数到摩尔分数的换算方法与一些技术专著不同。

相比而言，国家标准比技术专著的权威性更强，因此，若基础数据(如摩尔质量、体积发热量)、计算方法(如密度、相对密度的计算)在国家标准中已有规定，那么，应优先执行国家标准。

2.1.2 已知条件为：a) 绝对压力P 、热力学温度T 、组分数N ； b) 各组分的摩尔分数X i ，i = 1～N ；c) 可查GB ／T 17747．2附录B 中表B1、B2、B3得到的数据：● 58种物质的状态方程参数a n ，b n ， c n ，k n ，u n ，g n ，q n ，f n ，s n ，w n ；● 21种识别组分的特征参数M i ，E i ，K i ，G i ，Q i ，F i ，S i ，W i ；● 21种识别组分的二元交互作用参数E ij ∗，U ij ，K ij ，G ij ∗ 。

2.1.3 待求量：压缩因子 Z 。

2.1.4 计算步骤a) 计算第二维利系数B(1个值)；b) 计算系数C n ∗，n = 13～58，共46个值； c) 计算混合物体积参数K(1个值)； d) 形成压力的状态方程；e) 解压力的状态方程，求得压缩因子Z 。

2.2 计算第二维利系数B二元参数E ij ，G ij 的计算式为：E ij =E ij ∗(E i E j )12式(1) G ij =G ij ∗(G i +G j )2⁄式(2)除了GB ／T 17747．2附录B 中表B3中给出的E ij ∗，U ij ，K ij ，G ij ∗外，所有其他二元交互作用参数E ij ∗，U ij ，K ij ，G ij ∗的值都是1.0 。

B nij ∗的计算公式为：B nij ∗=(G ij +1−g n )g n∙(Q i Q j +1−q n )q n∙(F i 12F j 12+1−f n )f n∙(S i S j +1−S n )S n∙(W i W j +1−W n )W n式(3)B 按照式（4）计算：B =∑a n T−u n18n=1∑∑x i x j B nij ∗Nj=1N i=1E ij u n (K i K j )32 式(4)2.3 计算系数C *nC n ∗(n =13∼58)由式（5）计算： C n ∗=a n (G +1−g n )g n ∙(Q 2+1−q n )q n ∙(F +1−f n )f n ∙U u n T −u n 式(5)式（5）中的U,G,Q,F 只与天然气的组成有关，按式(6)～(9)计算：U 5=[∑x i E i 52Ni=1]2+2∑∑x i x j (U ij 5−1)(E i E j)52Nj=i+1N−1i=1式(6)G =∑x i G i +2∑∑x i x j (G ij ∗−1)(G i +G j )Nj=i+1N−1i=1N i=1式(7) Q =∑x i Q i Ni=1式(8) F =∑x i 2F i Ni=1式(9)2.4 计算混合物体积参数KK 按下式计算：K 5=[∑x i K i 52Ni=1]2+2∑∑x i x j (K ij 5−1)(K i K j)52Nj=i+1N−1i=1式(10)2.5 形成压力的状态方程GB ／T 17747．2中式(3)变形为：p m Z −p (RT )=0式(11) 式中p m ----摩尔密度，kmol m 3⁄R----摩尔气体常数，MJ (kmol ·K)⁄，取8.31451×104 MJ (kmol ·K)⁄GB ／T 17747.2附录B 中式(B10)变形为：ρm [1+Bρm −ρr ∑C n ∗18n=13+∑C n ∗(b n −c n k n ρr kn )ρr bn exp(−c n ρr kn )58n=13]−pRT=0 式(12)将GB ／T 17747．2中式(2)代入式(12)并展开 得：ρm+Bρm2−ρm2K3∑C n∗18n=13+∑C n∗(b n−c n k n K3k nρm k n)K3b nρm b n exp(−c n K3k nρm k n) 58n=13−pRT=0整理得：ρm+(B−K3∑C n∗18n=13)ρm2+∑C n∗K3b n(b n−c n k n K3k nρm k n)ρm b n exp(−c n K3k nρm k n) 58n=13−p RT=0式(13)令：A1=B−K3∑C n∗18n=13A2n=C n∗K3b nA3n=c n k n K3k nA4n=c n K3k n则式(13)变形为：ρm+A1ρm2+∑A2n(b n−A3nρm k n)ρm b n exp(−A4nρm k n)58n=13−pRT=0式(14)2.6解方程，求得压缩因子Z2.6.1方程曲线的形状将式（14）左边用f(ρm)表示，通过计算，得到若干( ρm ,f(ρm))数对，进行描点，得到方程曲线的形状，见下图。

图：方程曲线的形状2.6.2求解方法的确定式（14）是一个超越方程，求解方法有二分法、牛顿法、近似牛顿法等。

经分析，方程解的区间可以确定，采用二分法比较简捷。

2.6.3求解步骤1)给出方程解的区间（a,b）对一般的p 、T 条件，Z 必然处于0.45和1.2之间。

因此，取：a =p (1.2RT )⁄=0.833p (RT )⁄b =p (0.4RT )⁄=2.5p (RT )⁄2) 对有根区间取中值c =(a +b )2⁄,计算f (ρm )的值f (c )由式(11)、(14)可得出：cZ −p (RT )⁄=f (c ) 式(15)由式(15)可得：Z =[f (c )+p (RT )⁄]c ⁄ 式(16) 令：p cal =cZRT式中p cal ----由Z 得出的压力的计算值，MPa若| p −p cal |<1.0×10−6，则达到精度，输出Z ，计算结束；否则，继续下一步3)。

3) 若f (c )>0，则将c 赋值给b ；否则，将c 赋值给a 。

转向步骤2) 。

3 算例用VB 编程进行计算，气体组成采用GB/T17747.2附录C 中表C1中1～6号气样的摩尔组成，将计算结果与附录C 中表C2给出的结果进行对比，见表3。