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

基于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是否成立，若成立，则有根