下载文档原格式
管输天然气贸易计量中压缩因子的计算肖迪;巩大利【摘要】管输天然气的贸易结算经常采用体积计量和能量计量两种方式,压缩因子作为计算参数直接影响到计量准确度.国家标准GB/T 17747提供了天然气压缩因子的两种计算方法:摩尔组成法和物性值法.目前国内管输天然气压力普遍在6 MPa 以上、12 MPa以下,在这种工况条件下,物性值法计算压缩因子与摩尔组成法计算结果偏差比较大,尤其是非烃含量高(高含N2或CO2)的气体,采用物性值法更需慎重.在管输天然气贸易计量中,应采用适用范围更广,计算精度更高的摩尔组成法；物性值法是在现场增设在线物性参数测量仪器而采用的简单方法,此方法适用于无法得到气体组成且对计量准确度要求不高的情况.【期刊名称】《油气田地面工程》【年(卷),期】2011(030)009【总页数】1页(P24)【关键词】天然气;压缩因子;摩尔组成法;物性值法【作者】肖迪;巩大利【作者单位】国家石油天然气大流量计量站;国家石油天然气大流量计量站【正文语种】中文近年来，我国天然气工业迅速发展，建设了一批管道工程项目，引进了多条跨国管道。
管输天然气的贸易结算经常采用体积计量和能量计量两种方式，压缩因子作为计算参数直接影响到计量准确度。
国家标准《天然气压缩因子的计算（GB/T17747-1999）》规定了天然气压缩因子的两种计算方法，通过对两种方法比较，可明确各自的适用范围，确保国家和企业的合法权益。
国家标准《天然气压缩因子的计算GB/T 17747）》提供了天然气的压缩因子的两种计算方法：摩尔组成法和物性值法。
摩尔组成法也叫详细特征法（源自AGA8-92DC），采用已知天然气的详细摩尔组成和相关压力、温度计算压缩因子；物性值法，又称为总体特征法（源自SGERG-88），通过获取天然气的高位发热量、相对密度、CO2含量和N2含量中任意3个变量作为输入变量的压缩因子计算方法。
利用物性值计算压缩因子时，GB/T 17747不推荐采用N2含量作为输入变量之一，只给出了前3个变量作为输入变量时的压缩因子计算方法。
天然气压缩因子的分析及其计算谢莉莉;刘劲松【摘要】根据天然气压缩因子的2种计算方法:用摩尔组成进行计算和用物怀值进行计算编制计算机程序,并运用此程序研究天然气压缩因子与温度、压力之间的关系.【期刊名称】《上海计量测试》【年(卷),期】2011(000)005【总页数】4页(P27-30)【关键词】天然气;压缩因子;计算方法【作者】谢莉莉;刘劲松【作者单位】上海公正燃气计量站;上海公正燃气计量站【正文语种】中文0 引言天然气是重要的能源之一，随着天然气贸易量的增加，其流量计量越来越被人们重视。
在天然气流量计量中，天然气压缩因子是决定其准确与否的关键因素之一。
天然气压缩因子是实际气体状态采用理想气态方程时引入的偏差修正系数。
实际上，符合理想气态方程的理想气体是不存在的，实验表明，只有在低压高温下实际气体才可以近似被看作理想气体。
由于实际气体与理想气体的差异，使得对气体流量测量的准确性和可靠性难以评价，特别是低温、高压管道气体流量的测量，在这种情况下，管道中的被测介质就不能用理想气体状态方程进行描述。
在高压、低温下，任何气体理想状态方程都会出现明显的偏差，而且压力越高，温度越低，这种偏差就越大，因而需要引入一个压缩校正因子Z来修正气体的状态方程，如式（1）所示。
因此，天然气压缩因子Z在天然气这一重要能源计量中起着举足轻重的作用。
虽然GB/T 17747-1999《天然气压缩因子的计算》对天然气压缩因子进行了详细的描述，但是国内大部分是使用超压缩因子来计算天然气流量，对于压缩因子大多是文献上查得的或是通过图表获得。
若是用图表方式，则整个计算过程不仅费时费力，而且计算误差大，结果不准确。
而国外的进口流量计，像压缩因子等技术核心不公开，因此有必要编制一套计算程序来计算天然气压缩因子，确保天然气流量计量的准确性。
本文将介绍程序的编制简要以及运用该程序研究压缩因子与温度、压力之间的关系，并对两种方法进行比较。
1 计算程序编制天然气压缩因子的计算方法有2种：用天然气的摩尔组成进行计算和用天然气的物性值进行计算。
1程序目的利用AGA8-92DC模型计算天然气的压缩因子，该程序主要应用于在输气和配气正常进行的压力P和温度T范围内的管输气的压缩因子计算2数学模型：AGA8-92DC模型2.1模型介绍此模型是已知气体详细的摩尔分数组成和相关压力、温度来计算气体压缩因子。
输入变量包括绝对压力、热力学温度和摩尔组成。
摩尔组成是以摩尔分数表示下列组分：CO2、N2、H2、CO、CH4、C2H6、C3H8、i-C4H10、n-C4H10、i-C5H12、n-C5H12、n-C6H14、n-C7H16、n-C8H18。
2.2 模型适用条件绝对压力：0MPa＜P＜12MPa热力学温度：263K≤T≤338K高位发热量：30MJ·m-3≤HS≤45 MJ·m-3 相对密度：0.55≤d≤0.80天然气中各组分的摩尔分数应在以下范围内：CH4：0.7≤xCH4≤1.0N2：0≤xN2≤0.20CO2：0≤xCO2≤0.20C2H6：0≤xC2H6≤0.10C3H8：0≤xC3H8≤0.035C4H10：0≤xC4H10≤0.015C5H12：0≤xC5H12≤0.005C6H14：0≤xC6H14≤0.001C7H16：0≤xC7H16≤0.0005C8H18和更高碳数烃类：C8H18：0≤xC8H18≤0.0005H2：0≤xH2≤0.10CO ：0≤x CO ≤0.03如果已知体积分数组成，则应将其换算成摩尔分数组成。
所有摩尔分数大于0.00005的组分都不可忽略。
2.3 模型描述2.3.1 已知条件绝对压力P 、热力学温度T 、组分数N ； 各组分的摩尔分数X i ，i = 1～N ； 查附表1、2、3得到的以下数据：58种物质的状态方程参数a n ，b n ， c n ，k n ，u n ，g n ，q n ，f n ，s n ，w n ； 14种识别组分的特征参数M i ，E i ，K i ，G i ，Q i ，F i ，S i ，W i ；14种识别组分的二元交互作用参数E ij ∗，U ij ，K ij ，G ij ∗。
基于ASPENHYSYS和MATLAB天然⽓液化流程的优化- 50 -技术交流⽯油和化⼯设备2014年第17卷基于ASPEN HYSYS和MATLAB天然⽓液化流程的优化⿅院卫，刘丽华，吴⽟庭，马重芳（北京⼯业⼤学环境与能源⼯程学院传热强化与过程节能教育部重点实验室，北京 100124）[摘要] 为设计⼀种更节能的⼩型天然⽓液化装置，本⽂通过MATLAB中的ActiveX组件将ASPEN HYSYS与MATLAB连接起来，MATLAB利⽤ASPEN HYSYS中的spreadsheet对流程的参数进⾏读取，利⽤MATLAB的计算能⼒对液化流程中的参数进⾏相关的计算并返回ASPEN HYSYS中进⾏验证，以液化率和⽐功耗为流程性能评价指标找到了参数的最优值，实现了流程的优化。
[关键词] 天然⽓液化；MATLAB；ASPEN HYSYS；液化率；⽐功耗；优化天然⽓是当今世界能源消耗中的重要组成部分，它与煤炭、⽯油并称为世界能源的三⼤⽀柱[1]。
我国存在⼤量的天然⽓资源，但由于我国天然⽓⽥具有零、散、⼩的特点，许多偏远和⼩⽓⽥的⾮常规天然⽓都没有得到有效的开采[2-6]。
据统计，我国⾮常规⼩⽓⽥⽐例占所有⽓⽥总量的86%，⽽这些⽓⽥⽬前不具备开采的条件，导致这部分天然⽓不能合理有效利⽤[7]。
⼩型液化装置和LNG ⾮管道运输从技术上打破了零散⽓⽥和边际⽓⽥进⼊天然⽓终端市场的屏障。
如今设计的⼩型天然⽓液化流程尚有进⼀步优化的空间，因其采⽤的⽅法⼤多是利⽤ASPEN HYSYS ⾃带的优化器或者设置步长法进⾏优化[8]，具有⼀定的局限性。
MATLAB 以COM 技术为基础[9]，⽀持ActiveX 组件，它具备强⼤的计算能⼒，我们通过ActiveX 组件将MATLAB 和ASPEN HYSYS 连接起来，在MATLAB 平台环境下实现对ASPEN HYSYS 流程的读写和程序控制[10]，将MATLAB 的计算能⼒和ASPEN HYSYS 的仿真模拟能⼒结合起来，实现了设计流程的全局的优化，并降低了流程的液化率和⽐功耗，实现了节能。
标准研究/StandardResearch天然气计量中物性计算方法适用性探讨连子超1杨妮2李学成3许佳4代晓雨5吴萍4（1.华北油田公司华港燃气集团；2.中国石油西南油气田公司华油公司重庆凯源石油天然气有限责任公司；3.国家石油天然气管网集团北方管道大庆输油气分公司；4.国家管网集团山东省分公司德州作业区；5.中国石油吐哈油田分公司工程技术研究院地面工程设计所）摘要：针对目前天然气体积计量中面临的物性参数计算问题，在GB/T 17747.2—2011和ISO 20765-2：2015的基础上，采用Matlab 软件建立AGA8-92DC 和GERG-2008状态方程天然气物性求解程序，以相对偏差（RD）和平均相对偏差（ARD）为评价指标，评估了两种方程在计算不同种类天然气物性上的准确性。
结果表明，在管输天然气压力0～10MPa、温度280～320K 的范围内，AGA8-92DC 和GERG-2008状态方程的计算结果准确度一致，ARD 均为0.03%；对于含重烃天然气，压力小于30MPa、温度250～500K 的范围内，GERG-2008状态方程的计算表现更优，压力大于30MPa，部分温度范围内AGA8-92DC 状态方程的计算表现更优；AGA8-92DC 状态方程和GERG-2008状态方程分别在计算高含硫天然气和液化天然气物性上具有优越性，但当含硫量和重烃含量较大时，偏差会显著增大。
研究结果可为天然气计量工作的持续推进提供实际参考。
关键词：天然气计量；AGA8-92DC 方程；GERG-2008方程；压缩因子；物性计算方法DOI :10.3969/j.issn.2095-1493.2024.01.014Research on the adaptability of physical property calculation method in natural gas measurementLIAN Zichao 1，YANG Ni 2，LI Xuecheng 3，XU Jia 4，DAI Xiaoyu 5，WU Ping 41Huagang Gas Group of Huabei Oilfield Company2Huayou Company Chongqing Kaiyuan Oil &Gas Co.，Ltd.，Southwest Oil and Gas Field Company，CNPC3Daqing Oil and Gas Transmission Company of North Pipeline Co.，Ltd.，PipeChina 4Dezhou Operation Area of Shandong Company，PipeChina5Surface Engineering Design of Engineer Technology Research Institute of Tuha Oilfield，CNPCAbstract：At present，based on GB/T 17747.2—2011and ISO 20765-2：2015，faced with the cal-culation problem of physical property in the volumetric measurement of natural gas，the Matlab soft-ware is used to establish the natural gas physical property solving programs for AGA8-92DC and GERG-2008equation of ing relative deviation （RD）and average relative deviation （ARD）as evaluation indexes，the accuracy of the two equations in calculating the physical properties of differ-ent kinds of natural gas is evaluated.The results show that when the pressure of pipeline natural gas ranges from 0MPa to 10MPa and the temperature ranges from 280K to 320K，the accuracy of AGA8-92DC and GERG-2008equation of state is consistent and ARD is 0.03%.For natural gas con-taining heavy hydrocarbon，the GERG-2008equation of state is performed better when the pressure is less than 30MPa and the temperature is ranges from 250K to 500K while the calculation performance of AGA8-92DC equation is better when the pressure is greater than 30MPa and some temperature第一作者简介：连子超，2018年毕业于河北工业大学（工商管理专业）省任丘市万丰佳园小区，062550。
1程序目的利用AGA8-92DC模型计算天然气的压缩因子，该程序主要应用于在输气和配气正常进行的压力P和温度T围的管输气的压缩因子计算2数学模型：AGA8-92DC模型2.1模型介绍此模型是已知气体详细的摩尔分数组成和相关压力、温度来计算气体压缩因子。
输入变量包括绝对压力、热力学温度和摩尔组成。
摩尔组成是以摩尔分数表示下列组分：CO2、N2、H2、CO、CH4、C2H6、C3H8、i-C4H10、n-C4H10、i-C5H12、n-C5H12、n-C6H14、n-C7H16、n-C8H18。
2.2 模型适用条件绝对压力：0MPa＜P＜12MPa热力学温度：263K≤T≤338K高位发热量：30MJ·m-3≤HS≤45 MJ·m-3 相对密度：0.55≤d≤0.80天然气中各组分的摩尔分数应在以下围：CH4：0.7≤xCH4≤1.0N2：0≤xN2≤0.20CO2：0≤xCO2≤0.20C2H6：0≤xC2H6≤0.10C3H8：0≤xC3H8≤0.035C4H10：0≤xC4H10≤0.015C5H12：0≤xC5H12≤0.005C6H14：0≤xC6H14≤0.001C7H16：0≤xC7H16≤0.0005C8H18和更高碳数烃类：C8H18：0≤xC8H18≤0.0005H2：0≤xH2≤0.10CO ：0≤x CO ≤0.03如果已知体积分数组成，则应将其换算成摩尔分数组成。
所有摩尔分数大于0.00005的组分都不可忽略。
2.3 模型描述2.3.1 已知条件绝对压力P 、热力学温度T 、组分数N ； 各组分的摩尔分数，i = 1～N ； 查附表1、2、3得到的以下数据：58种物质的状态方程参数，， ，，，，，，， ； 14种识别组分的特征参数，，，，，，， ； 14种识别组分的二元交互作用参数，，，。
2.3.2 待求量压缩因子 Z 2.3.3 计算步骤a) 第二维利系数B 的计算：318*2111B (K K )nN Nu n i j ijijn i j a Tx x B -====∑∑∑11*22(G 1g )(1)(F F 1f )(S S 1s )(WW 1w )nnn n ng q f s w nijij n i j n i jn i j n i j nB QQ q =+-+-+-+-+-二元参数E ij 和G ij ，由以下两式计算：1*2(E E )ij iji j E E =*()/2ij ij i j G G G G =+b) 计算系数，n = 13～58*2(1)()(1)n n n n n g q f u u n n n n n C a G g Q Q q F f U T -=+-+-+-用以下方程求解混合方程，计算混合物参数U ，G ，Q 。
55525221111(2(1)())i iijNN Ni i j i i j U x E U E E -===+=+-∑∑∑1*1112(1)()N N Ni i i jiji j i i j i G x G x x GG G -===+=+-+∑∑∑1Ni i i Q x Q ==∑21Ni i i F x F ==∑c) 计算混合物体积参数K ；55152522i 111[]2(K 1)(K K )NN Ni ii j ijj i i j i K x K x y -===+=+-∑∑∑d) 计算对比密度摩尔密度为：/(ZRT)m P ρ=式中，P 为绝对压力，Mpa ；R 为摩尔气体常数；T 为热力学温度，K 。
对比密度ρr 同摩尔密度ρm 相关：3r m K ρρ=e) 利用AGA8-92DC 方程，对压缩因子进行迭代计算1858**n n n n 13131(b c k )exp(c )n n n k b k m r nn r r r n n Z B C C ρρρρρ===+-+--∑∑迭代过程：给出Z0的初始值为1，先计算出ρm ，将ρm 、K 和已知量带入AGA8-92DC 方程方程，得到新的Z 值，当（Z-Z0）的绝对值小于0.000001时，停止迭代，得到Z 值。
3 程序代码function [ Z ] = YSYZ( T,p,x) %计算天然气给定组分的压缩因子% x 为天然气组分，按照CO2 N2 H2 CO CH4 C2H6 C3H8 i-C4H10 n-C4H10 i-C5H12 n-C5H12% n-C6H14 n-C7H16 n-C8H18的顺序输入 %T 为温度，单位为K %P 为压力，单位为兆帕 N=14; R=8.314;%状态参数值b=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,3,3,3,3, 3,3,3,3,3,3,4,4,4,4,4,4,4,5,5,5,5,5,6,6,7,7,8,8,8,9,9];k=[0,0,0,0,0,0,0,0,0,0,0,0,3,2,2,2,4,4,0,0,2,2,2,4,4,4,4,0,1,1,2, 2,3,3,4,4,4,0,0,2,2,2,4,4,0,2,2,4,4,0,2,0,2,1,2,2,2,2];c=[0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,0,1,1,1, 1,1,1,1,1,1,0,0,1,1,1,1,1,0,1,1,1,1,0,1,0,1,1,1,1,1,1];g=[0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0, 1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,0,0];f=[0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0, 0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];q=[0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0, 0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,1,0,0,1,0,0,0,0,0,1];s=[0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];w=[0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];a=[0.153832600,1.341953000,-2.998583000,-0.048312280,0.375796500,-1.589575000,...-0.053588470,0.886594630,-0.710237040,-1.471722000,1.321850350,-0.786659250,...2.291290e-9,0.157672400,-0.436386400,-0.044081590,-0.003433888,0.032059050,...0.024873550,0.073322790,-0.001600573,0.642470600,-0.416260100,-0.066899570,...0.279179500,-0.696605100,-0.002860589,-0.008098836,3.150547000,0.007224479,...-0.705752900,0.534979200,-0.079314910,-1.418465000,-5.99905e-17,0.105840200,...0.034317290,-0.007022847,0.024955870,0.042968180,0.746545300,-0.291961300,...7.294616000,-9.936757000,-0.005399808,-0.243256700,0.049870160,0.033733797,...1.874951000,0.002168144,-0.658716400,0.000205518,0.009776195,-0.020487080,...0.015573220,0.006862415,-0.001226752,0.002850908];u=[0.0,0.5,1.0,3.5,-0.5,4.5,0.5,7.5,9.5,6.0,12.0,12.5,-6.0,2.0,3.0,2.0,2.0,11.0,-0.5,...0.5,0.0,4.0,6.0,21.0,23.0,22.0,-1.0,-0.5,7.0,-1.0,6.0,4.0,1.0,9.0,-13.0,21.0,8.0,...-0.5,0.0,2.0,7.0,9.0,22.0,23.0,1.0,9.0,3.0,8.0,23.0,1.5,5.0,-0.5,4.0,7.0,3.0,0.0,1.0,0.0];%特征参数值M=[44.0100,28.0135,2.0159,28.0100,16.0430,30.0700,44.0970,58.1230 ,58.1230,72.1500,72.1500,86.1770,100.2024,114.2310];E=[241.960600,99.737780,26.957940,105.534800,151.318300,244.16670 0,298.118300,324.068900,337.638900,365.599900,370.682300,402.636293,4 27.722630,450.325022];G=[0.189065,0.027815,0.034369,0.038953,0.0,0.079300,0.141239,0.25 6692,0.281835,0.332267,0.366911,0.289731,0.337542,0.383381];Q=[0.690,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0];K=[0.4557489,0.4479153,0.3514916,0.4533894,0.4619255,0.5279209,0. 5837490,0.6406937,0.6341423,0.6738577,0.6798307,0.7175118,0.7525189,0 .7849550];F=zeros(N);S=zeros(N);W=zeros(N);F(3)=1;%交互作用参数值Ex=[1.0,1.022740,1.281790,1.5,0.960644,0.925053,0.960237,0.906849 ,0.897362,0.726255,0.859764,0.855134,0.831229,0.808310;...1.022740,1.0,1.086320,1.005710,0.971640,0.970120,0.945939,0.946914,0. 973384,0.959340,0.945520,1.0,1.0,1.0;...1.281790,1.086320,1.0,1.1,1.170520,1.164460,1.034787,1.3,1.3,1. 0,1.0,1.0,1.0,1.0;...1.5,1.005710,1.1,1.0,0.990126,1.0,1.0,1.0,1.0049,1.0,1.0,1.0,1. 0,1.0;...0.960644,0.971640,1.170520,0.990126,1.0,1.0,0.994635,1.019530,0 .989844,1.00235,0.999268,1.107274,0.88080,0.880973;...0.925053,0.970120,1.164460,1.0,1.0,1.0,1.022560,1.0,1.013060,1. 0,1.005320,1.0,1.0,1.0;...0.960237,0.945939,1.034787,1.0,0.994635,1.022560,1.0,1.0,1.0049 ,1.0,1.0,1.0,1.0,1.0;...0.906849,0.946914,1.3,1.0,1.019530,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;...0.897362,0.973384,1.3,1.0,0.989844,1.01306,1.0049,1.0,1.0,1.0,1 .0,1.0,1.0,1.0;...0.726255,0.959340,1.0,1.0,1.00235,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1 .0,1.0;...0.859764,0.945520,1.0,1.0,0.999268,1.00532,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;...0.855134,1.0,1.0,1.0,1.107274,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1 .0;...0.831229,1.0,1.0,1.0,0.88088,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1. 0;...0.808310,1.0,1.0,1.0,0.880973,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1 .0];Gx=[1.0,0.982746,1.0,1.0,0.807653,0.370296,1.0,1.0,1.0,1.0,1.0,1. 0,1.0,1.0;0.982746,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;1.0,1.0,1.0,1.0,1.95731,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;0.807653,1.0,1.95731,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1. 0;0.370296,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0];Ux=[1.0,0.835058,1.0,0.9,0.963827,0.969870,1.0,1.0,1.0,1.0,1.0,1. 066638,1.077634,1.088178;...0.835058,1.0,0.408838,1.0,0.886106,0.816431,0.915502,1.0,0.9935 56,1.0,1.0,1.0,1.0,1.0;...1.0,0.408838,1.0,1.0,1.156390,1.616660,1.0,1.0,1.0,1.0,1.0,1.0, 1.0,1.0;...0.9,1.0,1.0,1.0,1.156390,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;...0.963827,0.886106,1.156390,1.0,1.0,1.0,0.990877,1.0,0.992291,1. 0,1.00367,1.302576,1.191904,1.205769;...0.969870,0.816431,1.616660,1.0,1.0,1.0,1.065173,1.25,1.25,1.25,1.25,1.0,1.0,1.0;...1.0,0.915502,1.0,1.0,0.990877,1.065173,1.0,1.0,1.0,1.0,1.0,1.0, 1.0,1.0;...1.0,1.0,1.0,1.0,1.0,1.25,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;...1.0;...1.0,1.0,1.0,1.0,1.0,1.25,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;...1.0,1.0,1.0,1.0,1.00367,1.25,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;...1.066638,1.0,1.0,1.0,1.302576,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1 .0;...1.077634,1.0,1.0,1.0,1.191904,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1 .0;...1.088178,1.0,1.0,1.0,1.205769,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1 .0];Kx=[1.0,0.982361,1.0,1.0,0.995933,1.00851,1.0,1.0,1.0,1.0,1.0,0.9 10183,0.895362,0.881152;0.982361,1.0,1.03227,1.0,1.00363,1.00796,1.0,1.0,1.0,1.0,1.0,1. 0,1.0,1.0;1.0,1.03227,1.0,1.0,1.02326,1.02034,1.0,1.0,1.0,1.0,1.0,1.0,1.0 ,1.0;1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;0.995933,1.00363,1.02326,1.0,1.0,1.0,1.007619,1.0,0.997596,1.0,1.002529,0.982962,0.983565,0.982707;1.00851,1.00796,1.02034,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0 ,1.0;1.0,1.0,1.0,1.0,1.007619,0.986893,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1 .0;1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;1.0,1.0,1.0,1.0,0.997596,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;1.0,1.0,1.0,1.0,1.002529,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;0.910183,1.0,1.0,1.0,0.982962,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1 .0;.0;0.881152,1.0,1.0,1.0,0.982707,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1 .0];Z0=1;M0=0;B=0;%参数Eij、Gij以及第二维利系数的计算sum1=0;sum2=0;for n=1:18sum=0;ZJCS=(T^-u(n));for i=1:14for j=1:14Eij(i,j)=Ex(i,j)*(E(i)*E(j)^0.5);Gij(i,j)=Gx(i,j)*(G(i)+G(j))/2;Bij(i,j)=((Gij(i,j)+1-g(n))^g(n))*((Q(i)*Q(j)+1-q(n))^q(n))*(((F(i)*F(j)^0.5)+1-f(n))^f(n))*((S(i)*S(j)+1-s(n))^s(n))*((W(i)*W(j)+1-w(n))^w(n));sum=sum+x(i)*x(j)*Bij(i,j)*(K(i)*K(j)^1.5);endendB=B+a(n)*ZJCS*sum;end%K值的计算F0=0;G0=0;Q0=0;U0=0;for i=1:NF0=F0+(x(i)^2)*F(i);Q0=Q0+x(i)*Q(i);sum1=sum1+x(i)*G(i);sum2=sum2+x(i)*(E(i)^2.5);endfor i=1:(N-1)for j=2:NG0=G0+x(i)*x(j)*(Gx(i,j)-1)*(G(i)+G(j));U0=U0+((Ux(i,j)^5)-1)*((E(i)*E(j))^2.5);endendG0=sum1+2*G0;U0=((sum2^2)+2*U0)^0.2;sum1=0;for i=1:Nsum1=sum1+x(i)*(K(i)^2.5);endsum2=0;for i=1:(N-1)for j=2:Nsum2=sum2+x(i)*x(j)*((Kx(i,j)^5)-1)*((K(i)*K(j))^2.5); endendK0=(((sum1^2)+2*sum2)^0.2);%对比密度的计算pr=p*(K0^3)/(Z0*R*T);SUM1=0;SUM2=0;%计算Cn的值，共有46个值for n=13:18Cn=a(n)*((G0+1-g(n))^g(n))*(((Q0^2)+Q0-q(n))^q(n))*((F0+1-f(n))^f(n))*(U0^u(n))*(T^-u(n));SUM1=SUM1+Cn;endfor i=13:58Cn=a(i)*((G0+1-g(i))^g(i))*((Q0*Q0+Q0-q(i))^q(i))*((F0+1-f(i))^f(i))*(U0^u(i))*(T^-u(i));SUM2=SUM2+Cn*(b(i)-c(i)*k(i)*(pr^k(i)))*(pr^b(i))*exp(-c(i)*(pr^k(i)));end% Z的计算Z=1+B*(pr/(K0^3))-pr*SUM1+SUM2;%迭代过程while abs(Z-Z0)>=0.000001Z0=(Z+Z0)/2;pr=p*(K0^3)/(Z0*R*T);SUM2=0;for n=13:58Cn=a(i)*((G0+1-g(i))^g(i))*((Q0*Q0+Q0-q(i))^q(i))*((F0+1-f(i))^f(i))*(U0^u(i))*(T^-u(i));SUM2=SUM2+Cn*(b(i)-c(i)*k(i)*(pr^k(i)))*(pr^b(i))*exp(-c(i)*(pr^k(i)));endZ=1+B*(pr/(K0^3))-pr*SUM1+SUM2;endZend4 运行结果对气样如下表的天然气计算在300K、0.1MPa下的压缩因子表4.1 气样组分气体组成摩尔分数CO2 0.006N2 0.003H2 0CO 0CH4 0.965C2H6 0.018C3H8 0.0045i-C4H10 0.001n-C4H10 0.001i-C5H12 0.0005n-C5H12 0.0003C6H14 0.0007C7H16 0C8H18 0附录附表1 状态参数方程表附表2 特征参数表附表3 二元交互作用参数值表。
