Matlab编程天然气压缩因子计算模型
MATLAB辅助《燃气输配》教学
由于各 地 的工况 状态不 同 ， 需要 将流 量转 化为 统一 仅使 同学 们省 去复 杂 的计 算 ， 深 了对 理 论 的认 识 ， 加 还 标准 下的气 体流量 。 目前在 国内 ， 一般 按 照 国际通用 准 可 以模 拟燃 气 管道 水力工 况 ， 对教 学起 到 了事半 功倍 的 则将 工况 状态下 的气 体 体 积 折算 为标 准 状 况 下 的气 体 作 用 。下 面 以美 国燃 气协 会 （ Ｇ 推 荐 压缩 系数 计算 Ａ Ａ） 体积 ， 后者 简称 为标 方 。在 压 力 不 太 高 、 度 不太 低 的 方法 为例 进行 阐述 。 温
一种新型天然气压缩因子数值计算方法
0. 99
1. 58
2. 59
4. 30
6. 00
9. 33
由表 3 可以看出，在对比温度为 1. 05 时，本文
计算方法与其他计算方法一样误差较大，最大误差
达到了 60% 。从压缩因子三维图可以看出，对比
温度在 1. 05 ～ 1. 10 之间时，Z 值曲面表现出了较
强的扭曲性，这也是造成各方法预测精度均较低的
近几年国内相继发现了一批高温高压天然气 田［8］。原有压缩因子计算方法适用压力范围低的 弊端逐渐暴露。石油大学郭绪强教授针对这一问 题进行了相关试验，取得了丰富的高压天然气实验 数据［9］。
将郭绪强教授发表的高压天然气实验数据与 传统天然气 压 缩 因 子 图 版［10］ 叠 加，发 现 天 然 气 压 缩因子在高压阶段具有较强的延展性，表现出了较 好的规律。利用三维绘图软件将数据进行处理，可 以发现天然气压缩因子曲面较为复杂。因此，本文 利用传统压缩因子图版与郭绪强教授发表的高压 天然气实验数据进行拟合，尝试找到高精度的能够
引言
目前使用较多的天然气压缩因子计算方法，包 括 Dranchk － Abu － Kassem 方 法 （ DAK） ［1 － 2］，Hankinson － Thomas － Phillips 方 法 （ HTP ） ［3］，Dranchuk － Purvis － Robinson 方法（ DPR） ［4］，以及由石 油大学李相方教授根据天然气压缩因子图版拟合 的李 相 方 方 法 （ LXF） 。这 些 计 算 公 式 均 是 根 据 Standing 和 KatZ 1942 年提出的压缩因子图版［5］采 用不同拟合方法拟合得到的［6］。在不同的对比压 力及对比温度下，误差均较大。根据李相方教授的 统计，各方法的最大误差均超过了 55%［7］。
压缩因子
lst2.additem " 75℃,Zn/Zg=1.0000
lst2.additem " "
lst2.additem "***(绝压)压力等级0.2Mpa**
lst2.additem "-20℃,Zn/Zg=1.0034
lst2.additem "-15℃,Zn/Zg=1.0032
lst2.additem " 15℃,Zn/Zg=1.0023
lst2.additem " 20℃,Zn/Zg=1.0021
lst2.additem " 25℃,Zn/Zg=1.0020
lst2.additem " 30℃,Zn/Zg=1.0019
lst2.additem " 35℃,Zn/Zg=1.0018
lst2.additem " 45℃,Zn/Zg=1.0000
lst2.additem " 50℃,Zn/Zg=1.0000
lst2.additem " 55℃,Zn/Zg=1.0000
lst2.additem " 60℃,Zn/Zg=1.0000
lst2.additem " 65℃,Zn/Zg=1.0000
lst2.additem "-10℃,Zn/Zg=1.0455
lst2.additem " -5℃,Zn/Zg=1.0426
lst2.additem " 0℃,Zn/Zg=1.0400
lst2.additem " 5℃,Zn/Zg=1.0375
管输天然气贸易计量中压缩因子的计算
管输天然气贸易计量中压缩因子的计算肖迪;巩大利【摘要】管输天然气的贸易结算经常采用体积计量和能量计量两种方式,压缩因子作为计算参数直接影响到计量准确度.国家标准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种：用天然气的摩尔组成进行计算和用天然气的物性值进行计算。
Matlab编程天然气压缩因子计算模型
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 ∗。
基于Matlab的燃气锅炉尾气气体状态计算
件。因为压力和温度较便于测量 , 所 以程序 是以 三种气 体的 质量 分数和 混合 气体的 总压 力和温 度为 输入 参数,
计算气体的比体积, 并且计算混 合气体两种状态间的焓差。
关键词: 锅炉尾气; 混合气体; 实际气体状态方程 ; 第二维 里系数
中图分类号: O 522 2
文献标识码: A
文章编号: 1007 7804( 2010) 03 0047 04
1 混合气体的实际气体状态方程和实 际气体焓变的算法
1 1 混合气体成分的确定
现假设锅炉内为完全燃烧, 且尾气中含氧量为
2% , 根据甲烷、一氧化碳、氢气燃烧的化学方程
式得到燃烧产物中氧气、氮气、二氧化碳、水蒸气
的摩尔分数分别为 x1、 x2、 x3、x4。
燃烧
CH 4 + 2O 2
CO2 + 2H 2O
( 2)
其中, Z 为压缩因子; B 为第二维里系数; p 为气体压 力, bar; T 为气体温度, K; v 为气体比体积, cm3 /m o;l
R 为气体常数。
收稿日期: 2010 04 06 基金项目: 上海市教委科技资助项目 ( 09ZZ161)
48
低温与特气
第 28卷
本文研究的混合气体由四种气体组成, 且为了 提高计算速度, 采用的维里方程截断至第二项。即
表 3 混合气体各组分的摩尔百分数
T able 3 M o le percent of the components of a gas m ix tures
m 氢气 : m一 氧化碳 : m 甲烷 x 1
x2
x3
x4
20% 20% 60 0 02 0 6947 0 058 0 2273 20% 30% 50 0 02 0 6916 0 0592 0 2292 20% 40% 40 0 02 0 6881 0 0606 0 2313 20% 50% 30 0 02 0 6837 0 0623 0 234 30% 40% 30 0 02 0 6796 0 045 0 2554 40% 30% 30 0 02 0 6768 0 0335 0 2697
基于ASPENHYSYS和MATLAB天然气液化流程的优化
基于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 的仿真模拟能⼒结合起来，实现了设计流程的全局的优化，并降低了流程的液化率和⽐功耗，实现了节能。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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 二元交互作用参数值表。