天然氣管道瞬態仿真研究綜述

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

(1.北京石油化工學院機械工程學院深水油氣管線關鍵技術與裝備北京市重點實驗室，北京 102617；2.長江大學石油工程學院，武漢 430100)

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

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

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

1 天然氣管道仿真的數學方程

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

表1 天然氣管道仿真數學模型

1.1 水力方程

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

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

1.2 熱力方程

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

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

1.3 組分輸運方程

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

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

2 天然氣管道仿真的數值求解方法

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

2.1 特征線法

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

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

(17)

(18)

(19)

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

2.2 有限差分法

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

圖2 有限差分法網格離散示意圖

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

(20)

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

(21)

(22)

(23)

2.3 有限體積法

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

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

步驟二基于連續性方程推導壓力修正方程。

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

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

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

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

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

步驟八求解密度。

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

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

2.4 有限元法

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

2.5 等效電路法

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

表2 電路與管路的類比

圖3 典型分支管網瞬態流動的等效電路圖

2.6 狀態空間模型法

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

(24)

(25)

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

2.7 低階模型法

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

2.8 人工智能法

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

圖4 天然氣管網的人工神經網絡結構圖

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

表3 不同求解方法的比較

3 天然氣管網仿真的加速方法

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

3.1 數學模型處理方面的加速方法

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

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

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

非線性的天然氣管道數學方程為

(26)

線性化的天然氣管道數學方程為

(27)

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

3.2 自適應仿真方面的加速方法

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

圖6 自適應仿真過程的設置

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

自適應的選擇數學模型方面，Domschke等[56]依據管道中瞬變程度的大小，合理簡化管道的數學模型，并通過控制管網元件連接處的誤差，來實現管道數學模型的自適應選擇。計算結果表明，該方法可提高求解效率60%以上。趙昆鵬[57]提出一種可依據實際測量的運行狀態參數，對管道數據(長度、內徑、粗糙度、高程等)進行自適應修正的自適應仿真方法，有效地提高仿真精度?？梢?，天然氣管道自適應仿真方法是一種精度高、計算速度快的方法。

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

天然氣管網由于拓撲結構復雜化、管網元件種類多樣化，仿真過程中離散方程組的系數矩陣為非主對角占優、非對稱、不規則結構的大型稀疏矩陣，呈現為分散的塊狀和帶狀結構，如圖7所示。求解該方程組是一個大規?？茖W計算問題，其求解速度直接決定了天然氣管網仿真效率。對于天然氣管網的仿真中不規則大型稀疏矩陣，早期研究一般采用高斯消元和LU分解等直接求解法[58-59]。眾所周知，矩陣直接解法的計算量正比于待求變量數目的三次方。因此，采用矩陣直接解法時，仿真效率嚴重受限于管網規模[58]。

圖7 某天然氣管網的離散方程系數矩陣結構示意圖

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

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

3.4 并行仿真的加速方法

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

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

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

圖8 天然氣管網GPU并行仿真示意圖

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

4 結論

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