摘 要：隨著城鎮(zhèn)燃?xì)鈿庠吹奶烊粴饣凸芫W(wǎng)的大型化，燃?xì)夤芫W(wǎng)水熱力計(jì)算時(shí)為滿足計(jì)算精度要求必須考慮燃?xì)獾膲嚎s性。本文采用基于SHBWR狀態(tài)方程求解城鎮(zhèn)燃?xì)獾膲嚎s因子，數(shù)值計(jì)算方法求解以燃?xì)饷芏葹樽宰兞康姆蔷€性方程。本文通過分析燃?xì)饷芏群瘮?shù)在鄰近有效解的范圍內(nèi)曲線的變化趨勢，選用牛頓迭代法、弦截法及二分法求解該方程，并編寫相應(yīng)的計(jì)算程序。通過算例分析，牛頓迭代法和弦截法是較好的求解方法，選取適當(dāng)?shù)某踔禃r(shí)，計(jì)算很少的次數(shù)（算例中約為2-6次）就可以達(dá)到同樣的精度。

關(guān)鍵詞：SHBWR狀態(tài)方程；城鎮(zhèn)燃?xì)?；燃?xì)饷芏?；壓縮因子；數(shù)值計(jì)算

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

據(jù)《中國城鄉(xiāng)建設(shè)統(tǒng)計(jì)年鑒2013》數(shù)據(jù)顯示，2013年底我國城鎮(zhèn)用氣人口數(shù)達(dá)4.08億，與2005年底相比增長了38.5%，而其中使用天然氣作為氣源的人口比例由24.1%增加到58.3%。隨著城市燃?xì)鈿庠刺烊粴饣内厔菀约安扇「邏和猸h(huán)儲(chǔ)氣調(diào)峰，很多大中城市出現(xiàn)了高壓、超高壓的燃?xì)夤芫W(wǎng)[1]。在高壓、超高壓運(yùn)行條件下進(jìn)行燃?xì)夤芫W(wǎng)水熱力計(jì)算時(shí)，若仍不考慮燃?xì)獾膲嚎s性，將不能滿足計(jì)算精度的要求。本文選用對天然氣狀態(tài)計(jì)算精度較高的SHBWR（BWRS）狀態(tài)方程確定城市燃?xì)獾膲嚎s因子Z[2][3]。

1 使用SHBWR狀態(tài)方程確定城鎮(zhèn)燃?xì)鈮嚎s因子Z的計(jì)算公式

城鎮(zhèn)燃?xì)馐且钥扇冀M分為主的混合氣體，可燃組分一般有碳?xì)浠衔铩浜鸵谎趸?，不可燃組分有二氧化碳、氮和氧等[4]。若燃?xì)庀到y(tǒng)的工作壓力為P和工作溫度為T，燃?xì)馐怯蒼種組分組成的混合氣體，其中第i組分的摩爾分?jǐn)?shù)為xi，則使用SHBWR狀態(tài)方程[5]確定燃?xì)鈮嚎s因子Z的計(jì)算公式見式（1）：

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

2 求解混合燃?xì)饷芏圈训臄?shù)值方法的選擇分析

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

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

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

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

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

根據(jù)上面的對比分析，本文選用二分法、牛頓法、弦截法迭代求解混合燃?xì)饷芏圈眩?jì)算終止條件選用，ε選用10-6。

3 數(shù)值方法求解混合燃?xì)饷芏鹊挠?jì)算算法

3.1 二分法求解燃?xì)饷芏鹊乃惴?/p>

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

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

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

④計(jì)算以及Fm=f（ρm）。

⑤判斷是否成立，若成立，則燃?xì)獾拿芏圈?=ρm；若不成立，轉(zhuǎn)⑥。

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

二分法求解燃?xì)饷芏鹊挠?jì)算框圖，見圖5。

3.2 牛頓迭代法求解燃?xì)饷芏鹊乃惴?/p>

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

②計(jì)算。

③判斷是否成立，若成立，則燃?xì)獾拿芏龋蝗舨怀闪?，則轉(zhuǎn)到④。

④ 判斷k+1

牛頓迭代法求解燃?xì)饷芏鹊挠?jì)算框圖，見圖6。

3.3 使用割線法求解混合燃?xì)饷芏鹊挠?jì)算算法

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

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

③計(jì)算。

④判斷是否成立，若成立，則燃?xì)獾拿芏?；若不成立，則轉(zhuǎn)到⑤。

⑤判斷k

弦截法求解燃?xì)饷芏鹊挠?jì)算框圖，見圖7。

4 算例計(jì)算結(jié)果

4.1 算例1

如表1給出的燃?xì)饨M分，當(dāng)工作壓力和溫度分別為1MPa和10℃、1MPa和30℃、1MPa和50℃、5MPa和35℃、10MPa和50℃時(shí)，燃?xì)饷芏群蛪嚎s因子的計(jì)算結(jié)果、與文獻(xiàn)計(jì)算結(jié)果對比結(jié)果見表2。

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

4.2 算例2

對應(yīng)于表1和表3中的燃?xì)饨M分，分別用二分法、牛頓迭代法、弦截法計(jì)算不同工作壓力和工作溫度時(shí)燃?xì)饷芏戎导皦嚎s因子，計(jì)算結(jié)果及迭代次數(shù)見表4。

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

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

5 結(jié)論

（1）基于SHBWR狀態(tài)方程求解混合燃?xì)鈮嚎s因子的過程中，關(guān)鍵的步驟是求解以混合燃?xì)饷芏圈褳樽宰兞康姆蔷€性方程，該方程需要使用數(shù)值方法進(jìn)行求解，本文通過對比分析，選用二分法、牛頓迭代法及弦截法分別編寫計(jì)算混合燃?xì)饷芏群蛪嚎s因子的計(jì)算程序；（2）通過算例分析，就對比迭代次數(shù)而言，牛頓迭代法和弦截法優(yōu)于二分法，牛頓迭代法優(yōu)于弦截法。但在每次迭代中牛頓迭代法需要計(jì)算兩個(gè)函數(shù)值和，所以就對比計(jì)算函數(shù)值次數(shù)而言，牛頓迭代法和弦截法優(yōu)于二分法，弦截法優(yōu)于牛頓迭代法；（3）通過算例分析，對于同種組分的燃?xì)?，為達(dá)到同樣的計(jì)算精度，壓力越高、溫度越低時(shí)需要的迭代次數(shù)越多。

參考文獻(xiàn)：

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

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

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

[4]段常貴.燃?xì)廨斉洌ǖ谖灏妫M].北京：中國建筑工業(yè)出版社，2011.

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

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

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

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

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

作者簡介：邊娟（1984-），女，江蘇淮安人，研究生，助教，研究方向：城市燃?xì)鈨?chǔ)存與運(yùn)輸技術(shù)。