基于SHBWR狀態(tài)方程求解城鎮(zhèn)燃氣壓縮因子的計算方法

摘 要：隨著城鎮(zhèn)燃氣氣源的天然氣化和管網(wǎng)的大型化，燃氣管網(wǎng)水熱力計算時為滿足計算精度要求必須考慮燃氣的壓縮性。本文采用基于SHBWR狀態(tài)方程求解城鎮(zhèn)燃氣的壓縮因子，數(shù)值計算方法求解以燃氣密度為自變量的非線性方程。本文通過分析燃氣密度函數(shù)在鄰近有效解的范圍內(nèi)曲線的變化趨勢，選用牛頓迭代法、弦截法及二分法求解該方程，并編寫相應(yīng)的計算程序。通過算例分析，牛頓迭代法和弦截法是較好的求解方法，選取適當?shù)某踔禃r，計算很少的次數(shù)（算例中約為2-6次）就可以達到同樣的精度。

關(guān)鍵詞：SHBWR狀態(tài)方程；城鎮(zhèn)燃氣；燃氣密度；壓縮因子；數(shù)值計算

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

據(jù)《中國城鄉(xiāng)建設(shè)統(tǒng)計年鑒2013》數(shù)據(jù)顯示，2013年底我國城鎮(zhèn)用氣人口數(shù)達4.08億，與2005年底相比增長了38.5%，而其中使用天然氣作為氣源的人口比例由24.1%增加到58.3%。隨著城市燃氣氣源天然氣化的趨勢以及采取高壓外環(huán)儲氣調(diào)峰，很多大中城市出現(xiàn)了高壓、超高壓的燃氣管網(wǎng)[1]。在高壓、超高壓運行條件下進行燃氣管網(wǎng)水熱力計算時，若仍不考慮燃氣的壓縮性，將不能滿足計算精度的要求。本文選用對天然氣狀態(tài)計算精度較高的SHBWR（BWRS）狀態(tài)方程確定城市燃氣的壓縮因子Z[2][3]。

1 使用SHBWR狀態(tài)方程確定城鎮(zhèn)燃氣壓縮因子Z的計算公式

城鎮(zhèn)燃氣是以可燃組分為主的混合氣體，可燃組分一般有碳氫化合物、氫和一氧化碳，不可燃組分有二氧化碳、氮和氧等[4]。若燃氣系統(tǒng)的工作壓力為P和工作溫度為T，燃氣是由n種組分組成的混合氣體，其中第i組分的摩爾分數(shù)為xi，則使用SHBWR狀態(tài)方程[5]確定燃氣壓縮因子Z的計算公式見式（1）：

根據(jù)上述求解過程，得到使用SHBWR狀態(tài)方程求解混合燃氣壓縮因子Z的計算框圖，見圖1。

2 求解混合燃氣密度ρ的數(shù)值方法的選擇分析

若有兩類燃氣，其組分分別見表1和表3，當燃氣的工作壓力P為4MPa、工作溫度T為-20℃時，密度函數(shù)曲線f（ρ）的變化趨勢，見圖2。從圖2中可看出該工況下兩類燃氣的密度函數(shù)方程f（ρ）=0有兩個實數(shù)解，但對于燃氣的密度而言，正實數(shù)軸上的解才可能是該方程的有效解。圖3給出了對應(yīng)于表1中的燃氣組分，工作壓力P分別為4MPa和1atm，工作溫度T分別為-20℃、0℃、20℃時，密度函數(shù)曲線f（ρ）在靠近有效根的區(qū)間范圍內(nèi)的變化趨勢。從圖中可以看出，在靠近有效根的區(qū)域范圍內(nèi)，f（ρ）曲線是單調(diào)上升的。為求解該非線性方程，下面分析二分法、不動點迭代法、牛頓法、弦截法等數(shù)值方法的可行性。

（1）二分法[7][8]：雖然二分法有迭代速度慢，但它在有根區(qū)間[ρ0，ρ1]內(nèi)一定可以求出密度的有效根，所以可以用二分法求解密度值，特別是在其它的求解方法失效的情況下。二分法的另一個優(yōu)點是可以事先知道求出滿足精度要求的近似根所需要的迭代次數(shù)，n次迭代后的誤差小于。但為開始二分法算法，必須先找到使f（ρ0）·f（ρ1）<0的有根區(qū)間[ρ0，ρ1]。

（2）不動點迭代法[7][8]：式（2）比較簡單的等價形式，見式（5）。因不動點迭代法局部收斂的條件為：在不動點ρ*的某鄰域內(nèi)連續(xù)，且，則迭代公式（6）局部收斂。圖4給出對應(yīng)于表1、表3中的燃氣組分，其工作溫度均為-20℃時，的函數(shù)曲線。從圖中看出，當燃氣密度ρ*比較大時，出現(xiàn)，即迭代序列不收斂。

（3）牛頓迭代法[7][8]：牛頓迭代法的迭代格式見式（8），密度函數(shù)表達式（2）的一階導(dǎo)數(shù)形式，見式（7）。根據(jù)牛頓迭代法的收斂條件，可得出密度函數(shù)f（ρ）在有效根ρ*鄰近區(qū)間是收斂的。牛頓迭代法若初值選取恰當（本文迭代初值ρ0取P/RT），收斂速度快，但每次迭代時不僅需要計算密度函數(shù)值f（ρk），還需要計算密度函數(shù)的一階導(dǎo)數(shù)值f`（ρk）。

（4）弦截法[7][8]：弦截法的迭代格式見式（9），在求解ρk+1時需要用到前面兩步的結(jié)果ρk和ρk-1，所以使用該方法需要先給出兩個初始值ρ0和ρ1，本文分別取0和P/RT[5]。根據(jù)弦截法的收斂條件，密度函數(shù)f（ρ）在有效根ρ*鄰近區(qū)間是收斂的。若初值選取恰當，弦截法具有超線性的收斂性。

根據(jù)上面的對比分析，本文選用二分法、牛頓法、弦截法迭代求解混合燃氣密度ρ，迭代計算終止條件選用，ε選用10-6。

3 數(shù)值方法求解混合燃氣密度的計算算法

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

①選取密度初值ρok，搜索步長取Δρ，搜索次數(shù)k=1。

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

③判斷f（ρok）·f（ρ1k）<0是否成立，若成立，則有根區(qū)間為[ρ0，ρ1]，轉(zhuǎn)④；若不成立，，轉(zhuǎn)②。

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

⑤判斷是否成立，若成立，則燃氣的密度ρ*=ρm；若不成立，轉(zhuǎn)⑥。

⑥判斷F0·Fm<0是否成立，若成立，則和；若不成立，則和，然后轉(zhuǎn)④。

二分法求解燃氣密度的計算框圖，見圖5。

3.2 牛頓迭代法求解燃氣密度的算法

①選取迭代初值為ρ0，迭代次數(shù)k=0。

②計算。

③判斷是否成立，若成立，則燃氣的密度；若不成立，則轉(zhuǎn)到④。

④ 判斷k+1

牛頓迭代法求解燃氣密度的計算框圖，見圖6。

3.3 使用割線法求解混合燃氣密度的計算算法

①選取迭代初值為ρ0和ρ1，迭代次數(shù)k=1。

②計算f（ρ0），f（ρ1），判斷，若成立，則交換ρ0，ρ1的順序。

③計算。

④判斷是否成立，若成立，則燃氣的密度；若不成立，則轉(zhuǎn)到⑤。

⑤判斷k

弦截法求解燃氣密度的計算框圖，見圖7。

4 算例計算結(jié)果

4.1 算例1

如表1給出的燃氣組分，當工作壓力和溫度分別為1MPa和10℃、1MPa和30℃、1MPa和50℃、5MPa和35℃、10MPa和50℃時，燃氣密度和壓縮因子的計算結(jié)果、與文獻計算結(jié)果對比結(jié)果見表2。

從表2中可以看出本文編寫的程序計算結(jié)果與文獻計算結(jié)果偏差不大，相對偏差小于0.1%。

4.2 算例2

對應(yīng)于表1和表3中的燃氣組分，分別用二分法、牛頓迭代法、弦截法計算不同工作壓力和工作溫度時燃氣密度值及壓縮因子，計算結(jié)果及迭代次數(shù)見表4。

算例2中，確定二分法的有根區(qū)間[ρ0，ρ1]時，搜索步長Δρ取值為1kg/m3，由（ε取10-6），得出迭代次數(shù)n大于19.9時才能滿足預(yù)定的精度要求，與程序計算次數(shù)20是吻合的。減少二分法的迭代次數(shù)可通過縮小有根區(qū)間的寬度，但這是以增加搜索有根區(qū)間的計算次數(shù)為代價的。

從表5和表6中可以看出，當三種方法預(yù)定的計算精度ε都取10-6時，各種工況下計算結(jié)果是一致的，但迭代次數(shù)是有差別的，牛頓迭代法和弦截法收斂速度快，2-6次即可達到精度要求。

5 結(jié)論

（1）基于SHBWR狀態(tài)方程求解混合燃氣壓縮因子的過程中，關(guān)鍵的步驟是求解以混合燃氣密度ρ為自變量的非線性方程，該方程需要使用數(shù)值方法進行求解，本文通過對比分析，選用二分法、牛頓迭代法及弦截法分別編寫計算混合燃氣密度和壓縮因子的計算程序；（2）通過算例分析，就對比迭代次數(shù)而言，牛頓迭代法和弦截法優(yōu)于二分法，牛頓迭代法優(yōu)于弦截法。但在每次迭代中牛頓迭代法需要計算兩個函數(shù)值和，所以就對比計算函數(shù)值次數(shù)而言，牛頓迭代法和弦截法優(yōu)于二分法，弦截法優(yōu)于牛頓迭代法；（3）通過算例分析，對于同種組分的燃氣，為達到同樣的計算精度，壓力越高、溫度越低時需要的迭代次數(shù)越多。

參考文獻：

[1]彭繼軍，楊昭，田貫三.高壓天然氣管網(wǎng)水力計算精度的研究[J].天然氣工業(yè)，2005，25（07）：96-98.

[2]劉燕.燃氣管網(wǎng)計算理論分析與應(yīng)用的研究（博士學(xué)位論文）[D].天津：天津大學(xué)，2004：26-29.

[3]王興畏.城市天然氣輸配管網(wǎng)水力模擬研究與實踐（博士學(xué)位論文）[D].重慶：重慶大學(xué)，2013：11-14，129-131.

[4]段常貴.燃氣輸配（第五版）[M].北京：中國建筑工業(yè)出版社，2011.

[5]苑偉民，賀三，袁宗明等.求解BWRS方程中密度根的數(shù)值方法[J].天然氣與石油，2009，27（01）：4-6.

[6]馬沛生.化工熱力學(xué)（通用型）[M].北京：化學(xué)工業(yè)出版社，2009.

[7]李慶揚，王能超，易大義.數(shù)值分析（第五版）[M].北京：清華大學(xué)出版社，2008.

[8]Curtis F.Gerald，PACTRICK O.WHEATLEY.應(yīng)用數(shù)值分析[M].呂淑娟，譯.北京：機械工業(yè)出版社，2006.

[9]吳玉國，陳保東.BWRS方程在天然氣物性計算中的應(yīng)用[J].油氣儲運，2003，22（10）：16-21.

作者簡介：邊娟（1984-），女，江蘇淮安人，研究生，助教，研究方向：城市燃氣儲存與運輸技術(shù)。