基于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） 　　③计算。 　　④判断是否成立，若成立，则燃气的密度；若不成立，则转到⑤。 　　⑤判断kNmax是否成立，若成立，则，转到③；若不成立，则调用二分法计算密度的程序。 　　弦截法求解燃气密度的计算框图，见图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状态方程求解混合燃气压缩因子的过程中，关键的步骤是求解以混合燃气密度ρ为自变量的非线性方程，该方程需要使用数值方法进行求解，本文通过对比分析，选用二分法、牛顿迭