天然氣管道瞬態(tài)仿真研究綜述

王 鵬，童?？担Y若蘭，周冠行，高智文，張引弟

(1.北京石油化工學(xué)院機(jī)械工程學(xué)院深水油氣管線關(guān)鍵技術(shù)與裝備北京市重點(diǎn)實(shí)驗(yàn)室，北京 102617；2.長江大學(xué)石油工程學(xué)院，武漢 430100)

天然氣作為一種優(yōu)質(zhì)、清潔、高效的化石能源，是21世紀(jì)能源結(jié)構(gòu)優(yōu)化的重要發(fā)展對象。管道是天然氣大規(guī)模輸送的主要方式。近年來中國天然氣管道業(yè)務(wù)得到大規(guī)模發(fā)展，取得了相當(dāng)豐碩的成果。截至2019年底，中國天然氣主干線管道和城市管道的里程分別超過7.7萬km和70萬km[1]，總里程可繞地球赤道19圈半。一張由橫跨東西、縱貫?zāi)媳?、覆蓋全國、連通海外的大型干線管網(wǎng)與多源、多向、多區(qū)、多級的復(fù)雜城市管網(wǎng)相結(jié)合的全國性天然氣管網(wǎng)已基本形成。此外，隨著物聯(lián)網(wǎng)、大數(shù)據(jù)、人工智能等新興技術(shù)的發(fā)展及其廣泛應(yīng)用，管道運(yùn)營商也在探索應(yīng)用新興技術(shù)提升油氣管道的安全水平和運(yùn)營效益[2]，智慧型管網(wǎng)應(yīng)運(yùn)而生?？梢?，大型智慧天然氣管網(wǎng)是未來發(fā)展的必然趨勢。

隨著天然氣管網(wǎng)的大型化、復(fù)雜化、智慧化和多氣源供應(yīng)的發(fā)展，供氣可靠性[3]、多氣源置換[4]、實(shí)時(shí)監(jiān)測分析[5]、智能優(yōu)化調(diào)度[6]等一系列新興課題相繼出現(xiàn)，對管網(wǎng)的設(shè)計(jì)、運(yùn)行、管理提出了諸多新的挑戰(zhàn)。管道仿真技術(shù)正是為了滿足現(xiàn)代化管網(wǎng)的需求而發(fā)展起來的[7]，已成為管網(wǎng)建設(shè)、運(yùn)行管理及科學(xué)研究等不可或缺的核心基礎(chǔ)技術(shù)。天然氣管網(wǎng)仿真通過計(jì)算機(jī)計(jì)算能夠準(zhǔn)確地再現(xiàn)管道內(nèi)部氣體壓力、流量和溫度等參數(shù)[8]，指導(dǎo)管網(wǎng)的設(shè)計(jì)規(guī)劃、運(yùn)營管理和智慧化發(fā)展。因此，在大型智慧管網(wǎng)發(fā)展的大背景下，梳理中外天然氣管道瞬態(tài)仿真的研究歷程，總結(jié)當(dāng)前的研究成果，分析其今后的發(fā)展趨勢具有重要意義。

由于實(shí)際輸氣管道內(nèi)的氣流總是受到供氣與用氣的波動、設(shè)備故障、控制系統(tǒng)的調(diào)節(jié)、壓縮機(jī)的啟停等多種因素影響，沿線氣體參數(shù)(壓力、流量、溫度等)時(shí)刻發(fā)生變化，屬于瞬態(tài)流動[9]。所以，現(xiàn)針對天然氣管道瞬態(tài)仿真，通過系統(tǒng)地文獻(xiàn)調(diào)研，從數(shù)學(xué)模型、求解方法和加速仿真3個(gè)方面綜述其研究進(jìn)展。

1 天然氣管道仿真的數(shù)學(xué)方程

數(shù)學(xué)模型是數(shù)值仿真的基礎(chǔ)。天然氣管道數(shù)學(xué)模型包括描述天然氣在管道內(nèi)部流動的連續(xù)性方程、動量方程和能量方程，其中連續(xù)性方程、動量方程主要求解流速、流量及壓力等表征水力性質(zhì)的參數(shù)，稱為水力方程；能量方程主要求解溫度和焓等表征熱力性質(zhì)的參數(shù)，稱為熱力方程。此外，隨著多氣源供應(yīng)的發(fā)展，氣體組分追蹤逐漸成為研究熱點(diǎn)，組分輸運(yùn)方程也成為天然氣管道仿真數(shù)學(xué)模型的重要組成部分，天然氣管道仿真數(shù)學(xué)模型的歸納總結(jié)如表1所示。以下從水力方程、熱力方程和組分輸運(yùn)方程3個(gè)方面分別介紹。

表1 天然氣管道仿真數(shù)學(xué)模型

1.1 水力方程

20世紀(jì)90年代以前的研究中，由于對管道流動規(guī)律認(rèn)識不足，常假設(shè)溫度變化對氣體流動的影響很小，可以忽略不計(jì)，天然氣管道仿真中只求解水力方程。同時(shí)，受計(jì)算機(jī)計(jì)算能力的限制，水力方程仍需進(jìn)行一定的簡化。水力方程的簡化主要針對動量方程。簡化的方法是，根據(jù)不同流動過程假設(shè)，直接忽略掉動量方程中瞬態(tài)項(xiàng)或慣性項(xiàng)。簡化后水力方程主要有拋物線模型[10]和雙曲線模型[11]兩種形式，即式(1)～式(3)。Osiadacz[12]通過數(shù)值計(jì)算全面分析了動量方程中各項(xiàng)隨氣體流動狀態(tài)的變化規(guī)律，指出：輸氣管道負(fù)荷變化較小時(shí)，瞬態(tài)項(xiàng)和慣性項(xiàng)對動量方程的貢獻(xiàn)程度很小，忽略這兩項(xiàng)對計(jì)算結(jié)果影響不大；輸氣管道負(fù)荷變化很大時(shí)，忽略這兩項(xiàng)會導(dǎo)致計(jì)算結(jié)果的精度較差。Abbaspour等[13]進(jìn)一步研究表明慣性項(xiàng)的影響隨瞬變劇烈程度增加而變大，忽略該項(xiàng)會導(dǎo)致數(shù)值結(jié)果將無法準(zhǔn)確地描述快瞬變流動中水力突變過程。

20世紀(jì)90年后，隨著計(jì)算機(jī)性能的提高，天然氣管道仿真所采用的動量方程不需要再進(jìn)行簡化，此時(shí)水力方程可寫為式(4)、式(5)[14-15]。這種數(shù)學(xué)模型延續(xù)了引入波速變量的做法，一定程度上簡化了求解過程，而且也能應(yīng)用到液體管網(wǎng)求解中。但是波速的求解是假設(shè)管道處于絕熱或等溫的基礎(chǔ)上進(jìn)行的，因此，這種數(shù)學(xué)模型的適用范圍受到一定限制。中國早期的天然氣管道仿真研究，如李長俊等[15]的研究，就采用這種方程形式。為了避免這種限制，研究者們直接求解原來的連續(xù)性方程和動量方程，并通過氣體狀態(tài)方程來保證水力方程封閉。此時(shí)，水力方程可寫為式(6)～式(8)[16-17]。狀態(tài)方程可通過PK、PR和BWRS等公式求解，具體公式可參考文獻(xiàn)[18]。這種水力方程不需對原方程進(jìn)行任何變換，適用范圍廣，而且方程中各項(xiàng)的物理意義明確，有利于數(shù)值計(jì)算程序的編寫。唐建峰等[16]的研究表明，這種方程形式可較精確地模擬大管徑、高壓力、長距離管道動態(tài)輸氣過程，并在哈依燃?xì)忾L輸管道末段儲氣研究中得到較好的應(yīng)用。這也是目前天然氣管道仿真中標(biāo)準(zhǔn)的水力方程形式。

1.2 熱力方程

早期的天然氣管道仿真研究中，常將管道流動假設(shè)為等溫或絕熱兩種狀態(tài)[12]。波速的求解也是基于是這兩個(gè)過程，即式(9)和式(10)。實(shí)際上，天然氣為可壓縮氣體，氣體狀態(tài)由壓力和溫度共同決定，忽略溫度對管道流動的影響很顯然與工程實(shí)際不符。為此，研究者們開始探討溫度對天然氣管道流動過程的影響。Keena[19]在天然氣管網(wǎng)仿真中增加能量方程的求解，指出溫度對管流水力參數(shù)有一定影響。Osiadacz等[20]的研究也表明當(dāng)天然氣溫度不穩(wěn)定時(shí)，采用等溫假設(shè)將產(chǎn)生較大的誤差。2000年后，熱力方程成為天然氣管道仿真數(shù)學(xué)模型的重要組成部分。根據(jù)所求解變量的不同，熱力方程可分為以焓為求解變量[21]和以溫度為求解變量[20]兩種形式，即式(11)和式(12)。這兩種方程只是數(shù)學(xué)表達(dá)式不同，本質(zhì)上是可以相互變換的。前一種方程形式是基于能量守恒原理直接建立得到的，方程中各項(xiàng)的物理意義明確，并在小型國產(chǎn)天然氣管道仿真軟件得到較廣泛的應(yīng)用，如文獻(xiàn)[21]中所開發(fā)的燃?xì)夤芫W(wǎng)瞬態(tài)模擬及泄漏分析軟件。后一種方程是基于前一種方程進(jìn)一步推導(dǎo)得到的，直接以溫度作為求解變量，可與實(shí)際工程中溫度邊界條件對應(yīng)，便于工程應(yīng)用。例如，大型國產(chǎn)天然氣管道仿真軟件RealPipe-Gas[22]的內(nèi)核就是采用這種熱力方程。

管道內(nèi)氣體與周圍環(huán)境間的換熱量的計(jì)算是熱力仿真中重要部分。一般認(rèn)為管道周圍環(huán)境(土壤、海水等)處于穩(wěn)態(tài)，即周圍環(huán)境溫度恒定。此時(shí)，管材與周圍環(huán)境介質(zhì)的導(dǎo)熱系數(shù)當(dāng)量成固定的集總傳熱系數(shù)，氣體與周圍環(huán)境間的換熱量采用傅里葉定律計(jì)算[23]，即式(13)。實(shí)際上，由于管道周圍環(huán)境的蓄熱效應(yīng)，周圍環(huán)境與管道中氣體進(jìn)行熱交換后，其溫度不是恒定值，而是以管道為軸線的徑向變化的溫度場。針對這一問題，Chaczykowski[24]提出采用一維軸對稱的徑向傳熱方程組來描述瞬態(tài)換熱過程，即式(14)，數(shù)值結(jié)果表明采用瞬態(tài)傳熱模型較穩(wěn)態(tài)傳熱模型的計(jì)算結(jié)果精度高。丁延鵬等[25]則進(jìn)一步指出穩(wěn)態(tài)傳熱模型高估了管線中流量的波動幅度。但是，瞬態(tài)熱力方程大大地增加管道仿真計(jì)算量，而且實(shí)際工程問題中管線沿線附近地層構(gòu)造以及地層性質(zhì)(溫度、傳熱特性等)的數(shù)據(jù)也很難確定。因此，實(shí)際工程應(yīng)用中土壤穩(wěn)態(tài)傳熱模型仍然更為常用。

1.3 組分輸運(yùn)方程

近年來，天然氣用量迅速增加，多氣源供應(yīng)的格局正在逐步形成。液化天然氣和沼氣、氫氣等可再生氣體開始并網(wǎng)輸送，管網(wǎng)中氣體組分發(fā)生重要變化，這對管網(wǎng)的運(yùn)營產(chǎn)生許多挑戰(zhàn)。研究瞬態(tài)流動中天然氣組分對流動過程影響的需求日益增加。Guandalini等[26]和Hafsi等[27]研究了天然氣管網(wǎng)注入氫氣對瞬態(tài)流動過程的影響，結(jié)果表明注氫對管網(wǎng)壓降的影響較小，但對流體的密度和速度產(chǎn)生較大影響，導(dǎo)致原來以體積或質(zhì)量的計(jì)量方式產(chǎn)生較大誤差。這一研究結(jié)果得到的驗(yàn)證。Chaczykowski等[28]提出在天然氣管道仿真中求解組分輸運(yùn)方程，以預(yù)測氣體組分隨流動的變化過程，同時(shí)指出組分輸運(yùn)方程本來應(yīng)采用對流擴(kuò)散方程描述，即式(15)；但是，對于輸氣管道中的流動問題，軸向和徑向的擴(kuò)散項(xiàng)很小，并且沒有源項(xiàng)和匯項(xiàng)，可進(jìn)一步簡化為對流方程，即式(16)。上述這些研究為天然氣管道瞬態(tài)流動中組分輸運(yùn)奠定了良好的開端。但是，在這些研究中忽略多組分氣體在軸向的擴(kuò)散效應(yīng)，與實(shí)際工程不符。此外，目前的研究還只針對單根管道開展，缺乏管網(wǎng)中連接點(diǎn)處組分輸運(yùn)方程，導(dǎo)致無法直接應(yīng)用于管網(wǎng)。

總結(jié)天然氣管道瞬態(tài)仿真的數(shù)學(xué)方程的研究進(jìn)展可知，水力和熱力方程的研究起步較早，研究成果較為豐富，其數(shù)學(xué)方程形式已經(jīng)得到公認(rèn)，壓力、流量和溫度等參數(shù)的預(yù)測結(jié)果較為準(zhǔn)確。但是，天然氣管道瞬態(tài)仿真中組分輸運(yùn)方程的研究仍處于起步階段，相應(yīng)的數(shù)學(xué)方程形式還有待進(jìn)一步完善。

2 天然氣管道仿真的數(shù)值求解方法

自19世紀(jì)60年代Wilkinson 等[29]首次實(shí)現(xiàn)瞬態(tài)管流控制方程的數(shù)值求解后，許多不同數(shù)值方法被提出用于天然氣管網(wǎng)仿真，目前可用于天然氣管網(wǎng)仿真的數(shù)值算法主要為以下8種：特征線法、有限差分法、有限容積法、有限元法、等效電路法、狀態(tài)空間模型法、低階模型法和智能算法。

2.1 特征線法

特征線法是一種基于特征理論的求解雙曲型偏微分方程組的方法，將偏微分方程化為沿特征線上的全微分問題，具有物理概念明確和穩(wěn)定性好等優(yōu)點(diǎn)。此外，由于求解過程為顯式求解，復(fù)雜管網(wǎng)的內(nèi)部連接點(diǎn)及外部氣源的處理較為簡單，易于計(jì)算程序編制，而且計(jì)算性能對管網(wǎng)規(guī)模的適應(yīng)性較好。特征線法離散網(wǎng)格示意圖如圖1所示。離散方程簡寫為式(17)～式(19)，其中點(diǎn)G、M和H處的值由插值得到。宋濤[30]將其應(yīng)用于中國某長輸管道工程中，采用SCADA系統(tǒng)實(shí)測數(shù)據(jù)驗(yàn)證其正確性。但是，該方法仿真過程中時(shí)間步長受相容性條件限制，一般取值較小，常用持續(xù)時(shí)間短的快瞬變流動研究。

i和j分別為空間節(jié)點(diǎn)和時(shí)層的編號；Δx和Δt分別為空間步長和時(shí)間步長；A為j時(shí)層的空間節(jié)點(diǎn)i；G、M、H分別為j-1時(shí)層的3個(gè)空間點(diǎn)；AG、AM、AH分別為順特征線、材料特征線和逆特征線

(17)

(18)

(19)

式中：A、B和C為離散方程系數(shù)；上標(biāo)+、0和-分別為順特征線、材料特征線和逆特征線。

2.2 有限差分法

有限差分法的基本過程是將管道離散為一系列微小管段，針對每個(gè)微小管段，基于Taylor級數(shù)展開，采用差分項(xiàng)近似代替微分項(xiàng)，從而將偏微分形式的管流控制方程轉(zhuǎn)化為代數(shù)方程組，求解代數(shù)方程組得到每個(gè)微小管段上的水熱力參數(shù)。差分近似代替微分的過程稱為離散。根據(jù)離散過程中采用的不同時(shí)層上的水熱力參數(shù)，有限差分法分為顯式[24]和隱式[31]兩種過程。顯式有限差法在離散時(shí)采用已經(jīng)上一時(shí)層的參數(shù)，其計(jì)算特性和適用特點(diǎn)與特征線法類似。與顯式法相反，隱式法在離散時(shí)采用待求時(shí)層的參數(shù)。因此，各個(gè)節(jié)點(diǎn)上的待求參數(shù)都相互關(guān)聯(lián)，需要聯(lián)立所有離散代數(shù)方程進(jìn)行統(tǒng)一求解，導(dǎo)致求解速度較慢。但是，這種方法的時(shí)間步長不受空間步長的限制，取值范圍很廣，對持續(xù)時(shí)間很長的瞬變過程有很好的適用性。需要說明的是，隱式有限差分法是目前天然氣管道仿真最為常用的求解方法。眾多國產(chǎn)天然氣管道仿真軟件，如EGPNS[15]、RealPipe-Gas[22]和川氣東送管道仿真軟件[32]等，都是基于這種求解方法開發(fā)內(nèi)核求解器。隱式有限差分法中，管道離散網(wǎng)格示意圖如圖2所示，天然氣管道控制方程可寫為式(20)，控制方程中各項(xiàng)可采用式(21)～式(23)的離散格式。式(20)～式(23)中A、B和C的具體表達(dá)式可參見文獻(xiàn)[21]。

圖2 有限差分法網(wǎng)格離散示意圖

天然氣管道控制方程簡寫為

(20)

隱式法中各項(xiàng)離散格式可寫為

(21)

(22)

(23)

2.3 有限體積法

有限體積法的基本過程是將管道離散為一系列微小管段，然后基于積分的方法，將偏微分形的管流控制方程轉(zhuǎn)化為離散的代數(shù)方程組并求解。該方法求解過程中需基于連續(xù)性方程推導(dǎo)壓力修正方程來實(shí)現(xiàn)速度與壓力的耦合，并采用計(jì)算流動力學(xué)中SIMPLE算法逐步求解動量方程和壓力修正方程[33]。該方法的主要步驟如下：

步驟一離散管道控制方程中動量方程[式(7)]。

步驟二基于連續(xù)性方程推導(dǎo)壓力修正方程。

步驟三上一時(shí)層速度，記為un-1，計(jì)算動量離散方程中的系數(shù)及常數(shù)項(xiàng)。

步驟四上一時(shí)層壓力分布，記為p*=pn-1。

步驟五求解動量方程，得u*。

步驟六求解壓力修正方程，得到修正壓力p′。

步驟七以p′修正速度和壓力，u*=u*+dep′，p*=p*+p′，其中de為壓力修正的系數(shù)。

步驟八求解密度。

步驟九判斷是否收斂。是un=u*,pn=p*，進(jìn)入下一時(shí)層計(jì)算；若否，利用改進(jìn)后的速度和壓力分布作為下一次迭代的計(jì)算初始值，重復(fù)步驟三～步驟八。

有限體積法采用隱式求解法，時(shí)間步長不受限制，離散過程中易采用絕對穩(wěn)定的高精度對流項(xiàng)離散格式，如TVD格式，具有穩(wěn)定好和計(jì)算精度高的優(yōu)點(diǎn)[33]，對于快慢瞬變流動都有很好適用性。但是，該方法求解過程復(fù)雜，尤其是邊界條件的處理。此外，該方法的壓力值和速度值不在同一套網(wǎng)格上，不利于工程應(yīng)用。

2.4 有限元法

有限元法的基本過程是將管道離散為一系列微小單元，對每一微小單元假定一個(gè)合適的形函數(shù)，即微小單元上流動參數(shù)的近似分布函數(shù)，然后在整個(gè)管道長度上采用Galerkin加權(quán)殘差法來最小化近似誤差，從而將原來的偏微分形式的管流控制方程轉(zhuǎn)化為代數(shù)方程組并求解。有限元方法的實(shí)施過程很復(fù)雜，有興趣的讀者可參見文獻(xiàn)[34-36]。文獻(xiàn)[34]對比了隱式有限差分法和有限元法的計(jì)算精度和計(jì)算速度，發(fā)現(xiàn)有限元法在計(jì)算精度上具有一定優(yōu)勢，但計(jì)算速度較慢。近年來，該方法的穩(wěn)定性和計(jì)算耗時(shí)得到較大的改善，對于大型復(fù)雜管網(wǎng)也有較好的適用性[35-36]。

2.5 等效電路法

等效電路法的基本思想是將管路與電路進(jìn)行等效[37]，如表2所示。壓強(qiáng)和電勢分別是驅(qū)動氣體和電流流動的動力，氣流和電流是兩種網(wǎng)絡(luò)中傳送的物質(zhì)的量化。流量隨時(shí)間的變化是由于流體的可壓縮性，可等效為電容效應(yīng)。管道的壓降是由于阻力效應(yīng)引起的壓降和通過管道的電感效應(yīng)引起的壓降之和。此時(shí)，瞬態(tài)流動的典型分支管網(wǎng)結(jié)構(gòu)的可以等效為圖3的電路。隨后，采用電路的相關(guān)理論分析管網(wǎng)，將偏微分管流控制方程轉(zhuǎn)化為一階常系數(shù)微分方程。這樣的形式避免了有限差分法和有體積法中的方程收斂問題，計(jì)算速度快，常應(yīng)用于管道的控制工程中。但該方法所求解的方程是電路相關(guān)，電路元件并不能完全代替管路元件，導(dǎo)致管網(wǎng)中多種現(xiàn)象無法被仿真，如壓縮機(jī)的喘振過程、溫度變化等，其工程應(yīng)用受到一定限制。

表2 電路與管路的類比

圖3 典型分支管網(wǎng)瞬態(tài)流動的等效電路圖

2.6 狀態(tài)空間模型法

狀態(tài)空間模型法的基本思想是：基于Laplace變換，將偏微分形的管流控制方程轉(zhuǎn)變與復(fù)數(shù)z相關(guān)的狀態(tài)函數(shù)，通過傳遞函數(shù)直接構(gòu)建管道入口和出口之間流動參數(shù)的關(guān)系式[38]。式(24)為微分表達(dá)式，式(25)為離散表達(dá)式。該方法計(jì)算精度可滿足工程需要，計(jì)算速度較有限差分法具有明顯的優(yōu)勢[39]，很適合應(yīng)用于管道控制和泄漏檢測等時(shí)效性強(qiáng)的工程問題[40]。但是該方法只求解管道入口和出口節(jié)點(diǎn)壓力或流量隨時(shí)間的變化趨勢，無法描述管道某一時(shí)刻沿線的壓力或流量參數(shù)。

(24)

(25)

式中：ξ=x/L；L為管道長度；F1,1～F2,2為相應(yīng)的系數(shù)；上劃線-為上時(shí)層。

2.7 低階模型法

低階模型法的基本思想是基于正交分解原理，采用只與時(shí)間有關(guān)的譜系數(shù)和只與空間有關(guān)的基函數(shù)相乘并線性疊加，以近似管道內(nèi)流動狀態(tài)的數(shù)值解[41]?；瘮?shù)與時(shí)間無關(guān)，可事先確定，管網(wǎng)瞬態(tài)仿真求解時(shí)只需求解數(shù)目極少的譜系數(shù)，達(dá)到降階的目的，可進(jìn)一步參考文獻(xiàn)[41]。該方法具有較高仿真效率，但其計(jì)算精度受基函數(shù)的影響。文獻(xiàn)[42]采用16組和32組基函數(shù)所得的壓力偏差達(dá)5%。另外，低階模型構(gòu)造過程復(fù)雜，實(shí)施難度大，而且沒有壓縮機(jī)和閥門等非管道元件處理的相關(guān)報(bào)道，工程應(yīng)用較少。

2.8 人工智能法

人工智能法是一種新興的管道瞬態(tài)仿真方法，目前仍處于起步階段[43]?，F(xiàn)有智能算法可以分為兩種：第一種是先采用管道數(shù)學(xué)模型及離散方法得到相應(yīng)的代數(shù)方程組，再采用智能優(yōu)化算法[如粒子群優(yōu)化(PSO)[44]]求解代數(shù)方程組；第二種是舍棄管流控制方程，直接搭建人工神經(jīng)網(wǎng)絡(luò)模型，以管道入口工況(流量、入口壓力、平均溫度等)和管道參數(shù)(管徑、管長、摩阻系數(shù)等)作為輸入量，以管道出口參數(shù)作為輸出量，基于歷史數(shù)據(jù)對該網(wǎng)絡(luò)進(jìn)行訓(xùn)練[45]。圖4為一個(gè)具有三層的人工神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖。數(shù)值結(jié)果表明，智能算法的計(jì)算精度較好，且具有高效的仿真效率，對于管道的控制工程具有較好的適用性[44-45]。但是，該算法目前只求解水力參數(shù)，未見有求解熱力方程和組分輸運(yùn)方程的相關(guān)報(bào)道。

圖4 天然氣管網(wǎng)的人工神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖

上述介紹的8種天然氣管網(wǎng)仿真方法的特點(diǎn)總結(jié)如表3所示，可以看出，目前工程應(yīng)用最為廣泛的方法是特性線法和有限差分法；等效電路法、狀態(tài)空間和低階模型法等對管道實(shí)時(shí)控制和泄露檢測工程具有較好的適用性。智能算法是一種新興的算法，仍處于起步階段。但是，該算法具有高仿真效率的優(yōu)點(diǎn)，將是以后數(shù)值求解方法的一個(gè)重要發(fā)展方向。

表3 不同求解方法的比較

3 天然氣管網(wǎng)仿真的加速方法

21世紀(jì)后，隨著“天然氣黃金時(shí)代”的到來[46]，管網(wǎng)規(guī)模迅速增長，天然氣管道仿真的求解效率越來越不能滿足管網(wǎng)高效節(jié)能優(yōu)化運(yùn)行等的需求。近二十年來，眾多學(xué)者開展一系列研究來加速仿真效率，其中效果較為明顯的方法主要可歸納為數(shù)學(xué)模型處理、自適應(yīng)仿真、離散方程求解、并行仿真等4個(gè)方面，以下分別介紹。

3.1 數(shù)學(xué)模型處理方面的加速方法

天然氣管道的數(shù)學(xué)為典型的非線性偏微分方程，是導(dǎo)致求解效率低的原因之一。早期天然氣管網(wǎng)數(shù)值仿真通過忽略動量方程中的慣性項(xiàng)[10]來避免該方程的非線性問題。此時(shí)天然氣管道數(shù)學(xué)模型可寫為式(1)。隨著能量方程的加入，水力熱力系統(tǒng)相互耦合，上述方法的加速效果大打折扣。為削弱這種耦合，研究者從大量的數(shù)值實(shí)例中總結(jié)出一種水力熱力去耦求解方法[47]。該方法將水力系統(tǒng)和熱力系統(tǒng)交替求解，通過外插的方式進(jìn)行數(shù)據(jù)傳遞，其求解流程如圖5所示。水力熱力去耦求解方法可在幾乎不影響計(jì)算精度的前提下提高求解效率。

圖5 水熱力去耦求解方法流程圖

水力熱力去耦的方式只弱化了水力熱力系統(tǒng)的耦合關(guān)系，并不改變數(shù)學(xué)模型方程形式的非線性，隨后所需求解的離散方程組仍為非線性方程。為此，張淼[48]針對管道離散方程組中的非線性項(xiàng)，基于Taylor級數(shù)展開，進(jìn)行線性化處理。然而，離散方程非線性項(xiàng)數(shù)目眾多，該方法實(shí)施過程較為煩瑣。針對這一不足，鄭建國等[22]直接對天然氣管道數(shù)學(xué)模型[式(26)]，基于矩陣Taylor級數(shù)進(jìn)行線性化，得線性化的方程[式(27)]。這種形式的數(shù)學(xué)方程可極大地簡化數(shù)學(xué)模型的離散及求解過程，并取得較好的加速效果。數(shù)值結(jié)果表明該方法可在幾乎不降低精度的前提下，提高求解效率5倍以上。數(shù)學(xué)模型線性化處理和水熱力去耦求解技術(shù)在大型國產(chǎn)天然氣管道仿真軟件RealPipe-Gas[22]得到較好的應(yīng)用，較好地保證其仿真高效性。

非線性的天然氣管道數(shù)學(xué)方程為

(26)

線性化的天然氣管道數(shù)學(xué)方程為

(27)

此外，數(shù)學(xué)模型選取不同的基礎(chǔ)求解變量也會影響求解效率。這一現(xiàn)象最早由Mahgerefteh等[49]和丁延鵬等[50]在特征線法管道仿真研究中發(fā)現(xiàn)。其原因在于，天然氣管道仿真所需求解的變量眾多，選取不同基礎(chǔ)求解變量時(shí)，數(shù)學(xué)模型方程形式簡潔程度差異較大，尤其是非線性項(xiàng)?；谶@一因素，為進(jìn)一步提高管仿真效率，研究者首先對單一求解變量的能量方程進(jìn)行研究，指出直接以溫度作為求解變量，求解速度最快[51]。隨后，Wang等[52]查明了壓力、流量、密度、流速等變量不同組合時(shí)對水力方程求解效率的影響機(jī)制，指出以密度和速度組合作為基本水力求解變量，可提高仿真效率約1.5倍。可見，建議天然氣管網(wǎng)仿真采用密度、速度和溫度作為基礎(chǔ)求解變量。

3.2 自適應(yīng)仿真方面的加速方法

天然氣管道仿真中一般要將偏微分形管流控制方程離散成代數(shù)方程組。傳統(tǒng)的做法是基于某個(gè)固定的數(shù)學(xué)模型，采用固定時(shí)間步長和空間網(wǎng)格系統(tǒng)進(jìn)行離散。這樣的做法往往出現(xiàn)如下情況：本可采用簡化數(shù)學(xué)模型和較大時(shí)間步長來快速計(jì)算的慢瞬變工況，卻采用了復(fù)雜數(shù)學(xué)模型和較小時(shí)間步長，導(dǎo)致計(jì)算資源被極大地浪費(fèi)。自適應(yīng)仿真技術(shù)正好可以解決上述問題。該技術(shù)可以根據(jù)管網(wǎng)內(nèi)天然氣流動變化的劇烈程度，自適應(yīng)地配置時(shí)空間網(wǎng)格系統(tǒng)和選優(yōu)數(shù)學(xué)模型，實(shí)現(xiàn)計(jì)算資源的“按需分配”，從而提高仿真效率。例如，在圖6所示天然氣管道瞬態(tài)流動案例中，可以在紅色標(biāo)識的流動變化劇烈時(shí)間段(0～12 h和60～72 h)采用小時(shí)步和密網(wǎng)格設(shè)置，以保證計(jì)算精度；而在藍(lán)色標(biāo)識的流動變化平緩時(shí)間段(24～48 h)采用大時(shí)步和稀網(wǎng)格設(shè)置，以提高仿真效率。

圖6 自適應(yīng)仿真過程的設(shè)置

自適應(yīng)的配置時(shí)空間網(wǎng)格系統(tǒng)方面，Tentis等[53]和梁炎明等[54]首次將自適應(yīng)網(wǎng)格控制策略引入天然氣管網(wǎng)仿真中，開發(fā)出一種自適應(yīng)的顯式有限差分法。計(jì)算結(jié)果表明，自適應(yīng)法在相同的計(jì)算精度下，相比固定空間步長法具有較高的計(jì)算速度。但是這種方法只針對空間網(wǎng)格系統(tǒng)進(jìn)行自適應(yīng)，時(shí)間步長受空間網(wǎng)格系統(tǒng)限制，在天然氣管網(wǎng)仿真過程采用隱式方法求解。為了充分自適應(yīng)方法的優(yōu)勢，Wang等[55]將“當(dāng)?shù)貢r(shí)層誤差控制”和“多層次網(wǎng)格系統(tǒng)”原理相結(jié)合，提出了一種適用于隱式求解過程的天然氣管網(wǎng)自適應(yīng)仿真方法，該方法實(shí)現(xiàn)了時(shí)間步長和空間網(wǎng)格系統(tǒng)同時(shí)且互不影響的自適應(yīng)控制。天然氣管道末端關(guān)閥算例表明，在相同的精度條件下，該方法可提高仿真效率約6倍。

自適應(yīng)的選擇數(shù)學(xué)模型方面，Domschke等[56]依據(jù)管道中瞬變程度的大小，合理簡化管道的數(shù)學(xué)模型，并通過控制管網(wǎng)元件連接處的誤差，來實(shí)現(xiàn)管道數(shù)學(xué)模型的自適應(yīng)選擇。計(jì)算結(jié)果表明，該方法可提高求解效率60%以上。趙昆鵬[57]提出一種可依據(jù)實(shí)際測量的運(yùn)行狀態(tài)參數(shù)，對管道數(shù)據(jù)(長度、內(nèi)徑、粗糙度、高程等)進(jìn)行自適應(yīng)修正的自適應(yīng)仿真方法，有效地提高仿真精度?？梢姡烊粴夤艿雷赃m應(yīng)仿真方法是一種精度高、計(jì)算速度快的方法。

3.3 離散方程求解方面的加速方法

天然氣管網(wǎng)由于拓?fù)浣Y(jié)構(gòu)復(fù)雜化、管網(wǎng)元件種類多樣化，仿真過程中離散方程組的系數(shù)矩陣為非主對角占優(yōu)、非對稱、不規(guī)則結(jié)構(gòu)的大型稀疏矩陣，呈現(xiàn)為分散的塊狀和帶狀結(jié)構(gòu)，如圖7所示。求解該方程組是一個(gè)大規(guī)模科學(xué)計(jì)算問題，其求解速度直接決定了天然氣管網(wǎng)仿真效率。對于天然氣管網(wǎng)的仿真中不規(guī)則大型稀疏矩陣，早期研究一般采用高斯消元和LU分解等直接求解法[58-59]。眾所周知，矩陣直接解法的計(jì)算量正比于待求變量數(shù)目的三次方。因此，采用矩陣直接解法時(shí)，仿真效率嚴(yán)重受限于管網(wǎng)規(guī)模[58]。

圖7 某天然氣管網(wǎng)的離散方程系數(shù)矩陣結(jié)構(gòu)示意圖

基于天然氣管網(wǎng)拓?fù)浣Y(jié)構(gòu)，設(shè)計(jì)較優(yōu)的矩陣消元策略，可一定程度地提高矩陣直接求解法效率，也可提高對管網(wǎng)拓?fù)浣Y(jié)構(gòu)復(fù)雜度的適應(yīng)性[60]。但是，當(dāng)管網(wǎng)規(guī)模較大時(shí)，該方法的復(fù)雜度迅速增加，導(dǎo)致收效甚微?；谙∈杈仃嚰夹g(shù)的離散方程組求解方法是另一種有效的加速手段，實(shí)例表明，這類方法計(jì)算效率較高，較傳統(tǒng)解法可提高仿真效率約800倍[61-63]。然而，稀疏矩陣技術(shù)中影響加快效果的不可控因素很多，已有數(shù)值實(shí)例表明，該技術(shù)在管網(wǎng)拓?fù)浣Y(jié)構(gòu)極端復(fù)雜的情況下，未必能獲得較好的加速效果[63]。

近年來，研究者們嘗試開發(fā)適用于上述矩陣結(jié)構(gòu)的迭代求解方法來提高管網(wǎng)仿真效率。王鵬[64]基于方程余量歸一化原則對系數(shù)矩陣進(jìn)行預(yù)處理，使其變?yōu)闂l件數(shù)較低的正定對稱結(jié)構(gòu)，進(jìn)而可采用高效共軛梯度法(CG)快速求解管網(wǎng)內(nèi)部的連接點(diǎn)方程。Madoliat等[43-44]將粒子群優(yōu)化(particle swarm optimization，PSO)算法和CA算法等智能迭代求解算法引入天然氣管網(wǎng)仿真中，較傳統(tǒng)矩陣直接解法可提高仿真效率200倍以上。但是由于離散方程矩陣結(jié)構(gòu)的特點(diǎn)及計(jì)算機(jī)舍入誤差等原因，實(shí)際計(jì)算過程中迭代法的求解效率隨求解規(guī)模增大而下降[65]。

3.4 并行仿真的加速方法

隨著天然氣管網(wǎng)規(guī)模越來越大型化和復(fù)雜化，管網(wǎng)瞬態(tài)仿真效率仍不能滿足操作培訓(xùn)、優(yōu)化運(yùn)行需求[7]。并行計(jì)算技術(shù)是解決該問題的最有效方法[7,66]。由于天然氣管網(wǎng)并行仿真技術(shù)對算法的并行性能有很高的要求，目前這方面的研究仍處于起步階段。

Fleming等[67]直接針對大型復(fù)雜天然氣管網(wǎng)仿真中的離散代數(shù)方程組，采用大型稀疏矩陣并行算法MUMPs 算法進(jìn)行求解。在24核CPU的并行求解過程中，較串行計(jì)算可加快約10倍。以上研究是代數(shù)方程組求解層面的細(xì)粒度并行，其加速效率對方程組并行求解算法的依賴性強(qiáng)，且實(shí)施及調(diào)優(yōu)難度大?？紤]管網(wǎng)中多元件結(jié)構(gòu)，以元件為求解任務(wù)的粗粒度并行，不僅可以降低并行難度，而且可以成倍的提高其加速效果。謝黛茜等[68]基于特征線方法求解管網(wǎng)中管道邊界的流量和壓力，使各相連的管道不再相互影響，再基于動態(tài)規(guī)劃的思想對天然氣管道的流量進(jìn)行分配，從而創(chuàng)新性地實(shí)現(xiàn)管網(wǎng)中各元件并行求解。但這種解耦方式引入顯式求解過程，時(shí)間步長受限。

宇波等[69]、Wang等[70-71]提出一種基于分而治之思想的天然氣管網(wǎng)仿真方法。該方法先一次性求解管網(wǎng)中所有元件連接點(diǎn)，將管網(wǎng)解耦成若干個(gè)可獨(dú)立求解的元件。該方法首次實(shí)現(xiàn)了隱式差分算法中管網(wǎng)中各元件的解耦。向月[72]進(jìn)一步在GPU并行平臺上進(jìn)行實(shí)施這一方法，原理如圖8所示，其仿真效率較國際著名商業(yè)管道仿真軟件快50倍以上。但是，目前的分而治之方法的天然氣管網(wǎng)并行仿真技術(shù)中均只求解了水力方程，并未涉及溫度方程的求解。此外，目前天然氣管網(wǎng)并行仿真中任務(wù)的優(yōu)化分配等調(diào)優(yōu)研究未見報(bào)道。

圖8 天然氣管網(wǎng)GPU并行仿真示意圖

總結(jié)天然氣管道瞬態(tài)仿真的加速方法可知，目前已從數(shù)學(xué)模型線性化處理、自適應(yīng)仿真、離散代數(shù)方程組求解和并行仿真4個(gè)方面得到很好的研究，但求解效率仍不能滿足當(dāng)前的大型復(fù)雜管網(wǎng)的發(fā)展需求，并行仿真技術(shù)解決該問題的最有效方法，將是將以后高效管網(wǎng)仿真的一個(gè)發(fā)展方向。

4 結(jié)論

天然氣管網(wǎng)仿真技術(shù)是管網(wǎng)建設(shè)、運(yùn)行管理及科學(xué)研究等不可缺的核心基礎(chǔ)技術(shù)。調(diào)研了中外天然氣管道瞬態(tài)仿真的研究進(jìn)展，總結(jié)發(fā)現(xiàn)目前的仿真技術(shù)在數(shù)學(xué)模型、求解方法和加快計(jì)算等方面已取得了較多成果，可較準(zhǔn)確地預(yù)測管道內(nèi)天然氣壓力、流速、流量、溫度等參數(shù)的變化，對中小型管網(wǎng)的仿真效率較高，在供氣調(diào)度、優(yōu)化運(yùn)行、實(shí)時(shí)監(jiān)控等方面都有較好的應(yīng)用。隨著天然氣管網(wǎng)大型化、復(fù)雜化和智慧化發(fā)展，多氣源供應(yīng)格局的逐步形成，目前的管道仿真將面臨一系列新興工況，多氣源的多組分仿真、人工智能仿真和高效并行仿真技術(shù)將是以后的重要發(fā)展方向。