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

(1.北京石油化工學(xué)院機械工程學(xué)院深水油氣管線關(guān)鍵技術(shù)與裝備北京市重點實驗室，北京 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īng)用新興技術(shù)提升油氣管道的安全水平和運營效益[2]，智慧型管網(wǎng)應(yīng)運而生?？梢姡笮椭腔厶烊粴夤芫W(wǎng)是未來發(fā)展的必然趨勢。

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

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

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

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

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

1.1 水力方程

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

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

1.2 熱力方程

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

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

1.3 組分輸運方程

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

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

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

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

2.1 特征線法

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

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

(17)

(18)

(19)

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

2.2 有限差分法

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

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

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

(20)

隱式法中各項離散格式可寫為

(21)

(22)

(23)

2.3 有限體積法

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

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

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

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

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

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

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

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

步驟八求解密度。

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

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

2.4 有限元法

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

2.5 等效電路法

等效電路法的基本思想是將管路與電路進行等效[37]，如表2所示。壓強和電勢分別是驅(qū)動氣體和電流流動的動力，氣流和電流是兩種網(wǎng)絡(luò)中傳送的物質(zhì)的量化。流量隨時間的變化是由于流體的可壓縮性，可等效為電容效應(yīng)。管道的壓降是由于阻力效應(yīng)引起的壓降和通過管道的電感效應(yīng)引起的壓降之和。此時，瞬態(tài)流動的典型分支管網(wǎng)結(jié)構(gòu)的可以等效為圖3的電路。隨后，采用電路的相關(guān)理論分析管網(wǎng)，將偏微分管流控制方程轉(zhuǎn)化為一階常系數(shù)微分方程。這樣的形式避免了有限差分法和有體積法中的方程收斂問題，計算速度快，常應(yīng)用于管道的控制工程中。但該方法所求解的方程是電路相關(guān)，電路元件并不能完全代替管路元件，導(dǎo)致管網(wǎng)中多種現(xiàn)象無法被仿真，如壓縮機的喘振過程、溫度變化等，其工程應(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á)式。該方法計算精度可滿足工程需要，計算速度較有限差分法具有明顯的優(yōu)勢[39]，很適合應(yīng)用于管道控制和泄漏檢測等時效性強的工程問題[40]。但是該方法只求解管道入口和出口節(jié)點壓力或流量隨時間的變化趨勢，無法描述管道某一時刻沿線的壓力或流量參數(shù)。

(24)

(25)

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

2.7 低階模型法

低階模型法的基本思想是基于正交分解原理，采用只與時間有關(guān)的譜系數(shù)和只與空間有關(guān)的基函數(shù)相乘并線性疊加，以近似管道內(nèi)流動狀態(tài)的數(shù)值解[41]?；瘮?shù)與時間無關(guān)，可事先確定，管網(wǎng)瞬態(tài)仿真求解時只需求解數(shù)目極少的譜系數(shù)，達(dá)到降階的目的，可進一步參考文獻(xiàn)[41]。該方法具有較高仿真效率，但其計算精度受基函數(shù)的影響。文獻(xiàn)[42]采用16組和32組基函數(shù)所得的壓力偏差達(dá)5%。另外，低階模型構(gòu)造過程復(fù)雜，實施難度大，而且沒有壓縮機和閥門等非管道元件處理的相關(guān)報道，工程應(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ò)進行訓(xùn)練[45]。圖4為一個具有三層的人工神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖。數(shù)值結(jié)果表明，智能算法的計算精度較好，且具有高效的仿真效率，對于管道的控制工程具有較好的適用性[44-45]。但是，該算法目前只求解水力參數(shù)，未見有求解熱力方程和組分輸運方程的相關(guān)報道。

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

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

表3 不同求解方法的比較

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

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

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

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

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

水力熱力去耦的方式只弱化了水力熱力系統(tǒng)的耦合關(guān)系，并不改變數(shù)學(xué)模型方程形式的非線性，隨后所需求解的離散方程組仍為非線性方程。為此，張淼[48]針對管道離散方程組中的非線性項，基于Taylor級數(shù)展開，進行線性化處理。然而，離散方程非線性項數(shù)目眾多，該方法實施過程較為煩瑣。針對這一不足，鄭建國等[22]直接對天然氣管道數(shù)學(xué)模型[式(26)]，基于矩陣Taylor級數(shù)進行線性化，得線性化的方程[式(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ù)學(xué)模型方程形式簡潔程度差異較大，尤其是非線性項?；谶@一因素，為進一步提高管仿真效率，研究者首先對單一求解變量的能量方程進行研究，指出直接以溫度作為求解變量，求解速度最快[51]。隨后，Wang等[52]查明了壓力、流量、密度、流速等變量不同組合時對水力方程求解效率的影響機制，指出以密度和速度組合作為基本水力求解變量，可提高仿真效率約1.5倍?？梢?，建議天然氣管網(wǎng)仿真采用密度、速度和溫度作為基礎(chǔ)求解變量。

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

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

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

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

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

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所示。求解該方程組是一個大規(guī)?？茖W(xué)計算問題，其求解速度直接決定了天然氣管網(wǎng)仿真效率。對于天然氣管網(wǎng)的仿真中不規(guī)則大型稀疏矩陣，早期研究一般采用高斯消元和LU分解等直接求解法[58-59]。眾所周知，矩陣直接解法的計算量正比于待求變量數(shù)目的三次方。因此，采用矩陣直接解法時，仿真效率嚴(yán)重受限于管網(wǎng)規(guī)模[58]。

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

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

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

3.4 并行仿真的加速方法

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

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

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

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

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

4 結(jié)論

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