基于高斯賽德爾迭代的輸氣管道狀態(tài)靜態(tài)模擬

文/王煒斌，何璐楠

0 引言

近年來，我國(guó)開建的天然氣管道數(shù)量日益增多、管道的長(zhǎng)度也急劇增加。如何保證輸送管道的穩(wěn)定可靠運(yùn)行是一項(xiàng)重要而艱巨的任務(wù)。當(dāng)輸送管道的流動(dòng)參數(shù)不隨時(shí)間變化時(shí)，意味著管道進(jìn)入穩(wěn)定運(yùn)行狀態(tài)。管道穩(wěn)態(tài)仿真原理是在已知的管道結(jié)構(gòu)，根據(jù)初始節(jié)點(diǎn)的壓力和載荷求解各個(gè)節(jié)點(diǎn)的壓力和流量。而針對(duì)天然氣管道穩(wěn)態(tài)的仿真，有助于輸送管道系統(tǒng)的優(yōu)化設(shè)計(jì)和管理。穩(wěn)態(tài)仿真不但可以減少投資浪費(fèi)，優(yōu)化系統(tǒng)設(shè)計(jì)，而且對(duì)于后期生產(chǎn)也有保障作用。

管道仿真技術(shù)就是為了發(fā)展現(xiàn)代管道技術(shù)應(yīng)運(yùn)而生的。計(jì)算機(jī)仿真技術(shù)就是采用計(jì)算機(jī)求解系統(tǒng)模型的數(shù)值解并對(duì)結(jié)果進(jìn)行分析論證的一門綜合性技術(shù)。仿真結(jié)果的優(yōu)劣取決于模型是否足夠精確（即理論模型與實(shí)際系統(tǒng)的誤差是否達(dá)到要求的精度）、是否能夠全面反映實(shí)際系統(tǒng)的特征。嚴(yán)格說來，對(duì)于氣體輸送管道的建模，流動(dòng)氣體的熱力學(xué)參數(shù)與運(yùn)動(dòng)學(xué)參數(shù)是相互關(guān)聯(lián)的。如果壓力和溫度較低，可以忽略壓力和溫度變化，將其當(dāng)作等溫氣體處理，簡(jiǎn)化計(jì)算復(fù)雜度。但對(duì)于高溫氣體必須采用特殊方法。此外，流動(dòng)的氣體常常處于不穩(wěn)定的狀態(tài)，若不穩(wěn)定狀態(tài)程度較小，可以將其蛻化為穩(wěn)定狀態(tài)處理。但當(dāng)存在較強(qiáng)不穩(wěn)定狀態(tài)時(shí)，必須采用動(dòng)態(tài)方法來處理。

目前在工業(yè)中應(yīng)用的天然氣管道計(jì)算模型，大多采用的穩(wěn)定流、簡(jiǎn)化的水力熱力公式，傳統(tǒng)的氣體狀態(tài)方程。由于這些公式是在很多理想化的條件下建立的，因而所得到的計(jì)算結(jié)果精度和使用性較差。根據(jù)大量工程實(shí)例，當(dāng)輸送管道的氣體壓力較小且溫度較低時(shí)，采用傳統(tǒng)的計(jì)算方法得到的結(jié)果與實(shí)際數(shù)據(jù)吻合度較高。但對(duì)于高溫高壓輸送管道，傳統(tǒng)的設(shè)計(jì)方法就失效了。不僅計(jì)算過程較為煩瑣，在計(jì)算精度方面也同樣無法滿足要求。這種情況在壓力超過6 Mpa時(shí)尤為明顯。因此急需尋找一種新的求解氣體管流偏微分方程組的有效方法。

本文選用高斯-賽德爾迭代（Gauss–Seidel method）求解線性方程組。該方法與雅可比法都采用了矩陣分解原理，但比雅克比法的求解速度快，而且省去了對(duì)原函數(shù)進(jìn)行微積分運(yùn)算，計(jì)算簡(jiǎn)便，在求解微分方程領(lǐng)域是一個(gè)重要的工具。

1 理論模型

當(dāng)管道內(nèi)的天然氣穩(wěn)定流動(dòng)時(shí)，管道內(nèi)的溫度T、氣體的壓力P和流速三者之間滿足連續(xù)性方程、能量方程、熱平穩(wěn)方程，具體如下所述。

氣體內(nèi)勢(shì)能的減少等于動(dòng)能、位勢(shì)能的增加與摩阻損失之和，因此滿足能量方程dP+pTds+pgdz+pdw2/2+0。其中，dP為氣體的內(nèi)勢(shì)能，Tds為摩阻損失，pgdz為氣體的位勢(shì)能，pdw2/2為氣體的動(dòng)能。P為氣體的壓強(qiáng)，單位為Pa。z為管道高度，單位為m。

天然氣在管道內(nèi)流動(dòng)滿足熱平穩(wěn)方程。dh+d(w2/2)/dx=dq/dx。其中，h為氣體的焓，單位為J/kg。q為單位質(zhì)量的氣體吸熱量，單位為J/kg。

沿線高程差小于200 m的輸氣管可以作為水平輸氣管來計(jì)算。定D義為管道內(nèi)徑，單位為m；λ為摩阻系，為無量綱參數(shù)。摩阻損失Tds可以簡(jiǎn)化為：由于h=f(P,T)，式（1）可以寫成：

將連續(xù)性方程、能量方程、熱平穩(wěn)方程及方程（2）進(jìn)行整理，寫成式（3）所示的矩陣形式。常微分方程方程組式（3）即為天然氣管道時(shí)的動(dòng)態(tài)方程。

求解式（3）所描述的常微分方程必須已知邊界條件。邊界條件分為初值問題和多點(diǎn)邊界兩類。本文所討論的是初值問題，即已知輸氣管道起始點(diǎn)位置的狀態(tài)參數(shù)值，由微分方程求管道其他位置的流量、壓力、溫度值。求解初始值問題的方法有很多，在此選用高斯賽德爾迭代法。該方法屬于顯式方法，是在雅克比矩陣迭代法的基礎(chǔ)上發(fā)展來的。采用高斯賽德爾迭代法進(jìn)行輸氣管靜態(tài)模擬的計(jì)算過程是依次計(jì)算公式（4）、公式（5）、公式（6）。其中公式（4）、公式（5）、公式（6）均是由公式（1）至公式（3）并結(jié)合高斯賽德爾迭代經(jīng)過推導(dǎo)演變獲得的，其參數(shù)含義、單位和之前的相同。

根據(jù)上述計(jì)算步驟，使用數(shù)值分析理論和Matlab編程，給定天然氣管道的初始數(shù)據(jù)和已知?dú)怏w、環(huán)境參數(shù)，即可求得管道沿線各截面的溫度、壓力和流速分布情況，同時(shí)可以利用Matlab的繪圖功能。將沿途壓力、溫度、流速的變化以折線圖的形式進(jìn)行表達(dá)，從而更加直觀，也便于在以后的工程應(yīng)用中，快速掌握管道情況，進(jìn)行實(shí)時(shí)監(jiān)測(cè)控制。

2 計(jì)算示例

某輸氣管道全長(zhǎng)957 km。設(shè)輸氣量為Q＝8x108 m3／年。起點(diǎn)壓力為8.9 MPa。管徑為635×8.74(24英寸)。天然氣的組成中甲烷為94.81%，乙炔為2.45%，二氧化碳為2.16%，氮?dú)鉃?.58%。土壤的導(dǎo)熱系數(shù)為1.2 W/(m·k)。管道埋深1.3 m。管道粗糙度0.016 38 mm。土壤月平均最低溫度5 ℃。

為了說明本文所提出的高斯賽德爾迭代方法的有效性，與雅克比迭代相比較求解計(jì)算實(shí)例。仿真結(jié)果給出了輸氣管沿線的溫度如圖1所示，輸氣管道沿線的壓力如圖2所示，輸氣管道沿線的流速圖3所示。

在圖1中實(shí)線代表實(shí)際溫度，點(diǎn)劃線代表采用雅克比方法求解的溫度，虛線代表采用本文方法求解的溫度。從圖1可知輸氣管道沿線的溫度隨輸送距離的增加不斷減小，本文所提出的方法求解的溫度更接近實(shí)際數(shù)據(jù)。在圖2中實(shí)線代表實(shí)際壓力，點(diǎn)劃線代表采用雅克比方法求解的壓力，虛線代表采用本文方法求解的壓力。從圖2可知輸氣管道沿線的壓力隨輸送距離的增加不斷減小，相比于雅克比方法，本文所提出的方法求解的壓力也更接近實(shí)際數(shù)據(jù)。在圖3中實(shí)線代表實(shí)際氣體流速，點(diǎn)劃線代表采用雅克比方法求解的氣體流速。虛線代表采用本文方法求解的氣體流速。從圖3可知輸氣管道沿線的流速隨輸送距離的增加而增加，相比于雅克比方法，本文所提出的方法求解的流速也更接近實(shí)際數(shù)據(jù)。

圖1 輸氣管道沿線的溫度

圖2 輸氣管道沿線的壓力

圖3 輸氣管道沿線的流速

綜上所述，輸氣管沿線壓力逐漸降低，流速逐漸減小，溫度逐漸降低，與雅克比求解方法相比，更接近實(shí)際情況。

3 結(jié)語

本文通過Matlab編程，采用高斯賽德迭代法對(duì)穩(wěn)態(tài)天然氣管道狀態(tài)進(jìn)行模擬，比四階龍格庫(kù)塔法計(jì)算復(fù)雜度低、計(jì)算時(shí)間短，比雅克比方法精度更高，為工程應(yīng)用提供了理論參考。同時(shí)，將該方法做進(jìn)一步改進(jìn)和拓展，使之適應(yīng)現(xiàn)有的自動(dòng)化設(shè)備，并將其應(yīng)用于長(zhǎng)距離天然氣管道輸送，為工業(yè)自動(dòng)化、現(xiàn)代化和管理的系統(tǒng)化、智能化保駕護(hù)航。