基于SHBWR状态方程求解城镇燃气压缩因子的计算方法

摘 要：随着城镇燃气气源的天然气化和管网的大型化，燃气管网水热力计算时为满足计算精度要求必须考虑燃气的压缩性。本文采用基于SHBWR状态方程求解城镇燃气的压缩因子，数值计算方法求解以燃气密度为自变量的非线性方程。本文通过分析燃气密度函数在邻近有效解的范围内曲线的变化趋势，选用牛顿迭代法、弦截法及二分法求解该方程，并编写相应的计算程序。通过算例分析，牛顿迭代法和弦截法是较好的求解方法，选取适当的初值时，计算很少的次数（算例中约为2-6次）就可以达到同样的精度。

关键词：SHBWR状态方程；城镇燃气；燃气密度；压缩因子；数值计算

DOI：10.16640/j.cnki.37-1222/t.2016.06.200

据《中国城乡建设统计年鉴2013》数据显示，2013年底我国城镇用气人口数达4.08亿，与2005年底相比增长了38.5%，而其中使用天然气作为气源的人口比例由24.1%增加到58.3%。随着城市燃气气源天然气化的趋势以及采取高压外环储气调峰，很多大中城市出现了高压、超高压的燃气管网[1]。在高压、超高压运行条件下进行燃气管网水热力计算时，若仍不考虑燃气的压缩性，将不能满足计算精度的要求。本文选用对天然气状态计算精度较高的SHBWR（BWRS）状态方程确定城市燃气的压缩因子Z[2][3]。

1 使用SHBWR状态方程确定城镇燃气压缩因子Z的计算公式

城镇燃气是以可燃组分为主的混合气体，可燃组分一般有碳氢化合物、氢和一氧化碳，不可燃组分有二氧化碳、氮和氧等[4]。若燃气系统的工作压力为P和工作温度为T，燃气是由n种组分组成的混合气体，其中第i组分的摩尔分数为xi，则使用SHBWR状态方程[5]确定燃气压缩因子Z的计算公式见式（1）：

根据上述求解过程，得到使用SHBWR状态方程求解混合燃气压缩因子Z的计算框图，见图1。

2 求解混合燃气密度ρ的数值方法的选择分析

若有两类燃气，其组分分别见表1和表3，当燃气的工作压力P为4MPa、工作温度T为-20℃时，密度函数曲线f（ρ）的变化趋势，见图2。从图2中可看出该工况下两类燃气的密度函数方程f（ρ）=0有两个实数解，但对于燃气的密度而言，正实数轴上的解才可能是该方程的有效解。图3给出了对应于表1中的燃气组分，工作压力P分别为4MPa和1atm，工作温度T分别为-20℃、0℃、20℃时，密度函数曲线f（ρ）在靠近有效根的区间范围内的变化趋势。从图中可以看出，在靠近有效根的区域范围内，f（ρ）曲线是单调上升的。为求解该非线性方程，下面分析二分法、不动点迭代法、牛顿法、弦截法等数值方法的可行性。

（1）二分法[7][8]：虽然二分法有迭代速度慢，但它在有根区间[ρ0，ρ1]内一定可以求出密度的有效根，所以可以用二分法求解密度值，特别是在其它的求解方法失效的情况下。二分法的另一个优点是可以事先知道求出满足精度要求的近似根所需要的迭代次数，n次迭代后的误差小于。但为开始二分法算法，必须先找到使f（ρ0）·f（ρ1）<0的有根区间[ρ0，ρ1]。

（2）不动点迭代法[7][8]：式（2）比较简单的等价形式，见式（5）。因不动点迭代法局部收敛的条件为：在不动点ρ*的某邻域内连续，且，则迭代公式（6）局部收敛。图4给出对应于表1、表3中的燃气组分，其工作温度均为-20℃时，的函数曲线。从图中看出，当燃气密度ρ*比较大时，出现，即迭代序列不收敛。

（3）牛顿迭代法[7][8]：牛顿迭代法的迭代格式见式（8），密度函数表达式（2）的一阶导数形式，见式（7）。根据牛顿迭代法的收敛条件，可得出密度函数f（ρ）在有效根ρ*邻近区间是收敛的。牛顿迭代法若初值选取恰当（本文迭代初值ρ0取P/RT），收敛速度快，但每次迭代时不仅需要计算密度函数值f（ρk），还需要计算密度函数的一阶导数值f`（ρk）。

（4）弦截法[7][8]：弦截法的迭代格式见式（9），在求解ρk+1时需要用到前面两步的结果ρk和ρk-1，所以使用该方法需要先给出两个初始值ρ0和ρ1，本文分别取0和P/RT[5]。根据弦截法的收敛条件，密度函数f（ρ）在有效根ρ*邻近区间是收敛的。若初值选取恰当，弦截法具有超线性的收敛性。

根据上面的对比分析，本文选用二分法、牛顿法、弦截法迭代求解混合燃气密度ρ，迭代计算终止条件选用，ε选用10-6。

3 数值方法求解混合燃气密度的计算算法

3.1 二分法求解燃气密度的算法

①选取密度初值ρok，搜索步长取Δρ，搜索次数k=1。

②令ρ1k=ρok+Δρ，计算f（ρok），f（ρ1k）。

③判断f（ρok）·f（ρ1k）<0是否成立，若成立，则有根区间为[ρ0，ρ1]，转④；若不成立，，转②。

④计算以及Fm=f（ρm）。

⑤判断是否成立，若成立，则燃气的密度ρ*=ρm；若不成立，转⑥。

⑥判断F0·Fm<0是否成立，若成立，则和；若不成立，则和，然后转④。

二分法求解燃气密度的计算框图，见图5。

3.2 牛顿迭代法求解燃气密度的算法

①选取迭代初值为ρ0，迭代次数k=0。

②计算。

③判断是否成立，若成立，则燃气的密度；若不成立，则转到④。

④ 判断k+1

牛顿迭代法求解燃气密度的计算框图，见图6。

3.3 使用割线法求解混合燃气密度的计算算法

①选取迭代初值为ρ0和ρ1，迭代次数k=1。

②计算f（ρ0），f（ρ1），判断，若成立，则交换ρ0，ρ1的顺序。

③计算。

④判断是否成立，若成立，则燃气的密度；若不成立，则转到⑤。

⑤判断k

弦截法求解燃气密度的计算框图，见图7。

4 算例计算结果

4.1 算例1

如表1给出的燃气组分，当工作压力和温度分别为1MPa和10℃、1MPa和30℃、1MPa和50℃、5MPa和35℃、10MPa和50℃时，燃气密度和压缩因子的计算结果、与文献计算结果对比结果见表2。

从表2中可以看出本文编写的程序计算结果与文献计算结果偏差不大，相对偏差小于0.1%。

4.2 算例2

对应于表1和表3中的燃气组分，分别用二分法、牛顿迭代法、弦截法计算不同工作压力和工作温度时燃气密度值及压缩因子，计算结果及迭代次数见表4。

算例2中，确定二分法的有根区间[ρ0，ρ1]时，搜索步长Δρ取值为1kg/m3，由（ε取10-6），得出迭代次数n大于19.9时才能满足预定的精度要求，与程序计算次数20是吻合的。减少二分法的迭代次数可通过缩小有根区间的宽度，但这是以增加搜索有根区间的计算次数为代价的。

从表5和表6中可以看出，当三种方法预定的计算精度ε都取10-6时，各种工况下计算结果是一致的，但迭代次数是有差别的，牛顿迭代法和弦截法收敛速度快，2-6次即可达到精度要求。

5 结论

（1）基于SHBWR状态方程求解混合燃气压缩因子的过程中，关键的步骤是求解以混合燃气密度ρ为自变量的非线性方程，该方程需要使用数值方法进行求解，本文通过对比分析，选用二分法、牛顿迭代法及弦截法分别编写计算混合燃气密度和压缩因子的计算程序；（2）通过算例分析，就对比迭代次数而言，牛顿迭代法和弦截法优于二分法，牛顿迭代法优于弦截法。但在每次迭代中牛顿迭代法需要计算两个函数值和，所以就对比计算函数值次数而言，牛顿迭代法和弦截法优于二分法，弦截法优于牛顿迭代法；（3）通过算例分析，对于同种组分的燃气，为达到同样的计算精度，压力越高、温度越低时需要的迭代次数越多。

参考文献：

[1]彭继军，杨昭，田贯三.高压天然气管网水力计算精度的研究[J].天然气工业，2005，25（07）：96-98.

[2]刘燕.燃气管网计算理论分析与应用的研究（博士学位论文）[D].天津：天津大学，2004：26-29.

[3]王兴畏.城市天然气输配管网水力模拟研究与实践（博士学位论文）[D].重庆：重庆大学，2013：11-14，129-131.

[4]段常贵.燃气输配（第五版）[M].北京：中国建筑工业出版社，2011.

[5]苑伟民，贺三，袁宗明等.求解BWRS方程中密度根的数值方法[J].天然气与石油，2009，27（01）：4-6.

[6]马沛生.化工热力学（通用型）[M].北京：化学工业出版社，2009.

[7]李庆扬，王能超，易大义.数值分析（第五版）[M].北京：清华大学出版社，2008.

[8]Curtis F.Gerald，PACTRICK O.WHEATLEY.应用数值分析[M].吕淑娟，译.北京：机械工业出版社，2006.

[9]吴玉国，陈保东.BWRS方程在天然气物性计算中的应用[J].油气储运，2003，22（10）：16-21.

作者简介：边娟（1984-），女，江苏淮安人，研究生，助教，研究方向：城市燃气储存与运输技术。