天然氣管網仿真分層計算方法

閻 濤，劉天堯，王雪莉，薛向東，2，康 陽，朱 峰

1.國家管網集團科學技術研究總院分公司，河北廊坊 065000

2.西安交通大學軟件學院，陜西西安 710049

隨著經濟的持續發展，我國對能源的需求越來越大[1]。由于天然氣具有高熱值、低碳、清潔等優點[2-3]，其消費量逐年增長。預計到2030 年，全國天然氣消費量將達到（5 500～6 000）×108m3[4]。出于經濟原因，天然氣通常通過管道進行大規模輸送[5-6]。通過對天然氣管網進行仿真，可以明確天然氣在管網內部流動過程中的水熱力變化情況，為天然氣管網的優化、管理、調度及運行提供技術支持[7-10]。截至2021年，我國主干天然氣管道總里程達到11.6×104km[11]。天然氣管網規模的不斷擴大對管網仿真提出了越來越高的要求。對于大規模天然氣管網，采用隱式差分法離散水熱力系統得到的代數方程組規模較大。受限于代數方程組巨大的規模，LU 分解等直接解法的[12]求解速度慢。童?？档萚13]采用迭代法替代高斯消元法來求解水熱力方程組，有效提高了求解速度。然而，使用迭代法仍要對整個方程組進行求解操作，這對大規模管網來說仍是個較大的挑戰。為提高天然氣管網的求解速度與效率，提出了一種管網分層求解的算法。該算法將管網劃分為不同層次的子管網，自下而上傳遞變量關系式，自上而下傳遞變量值，最終實現整個管網仿真過程的快速求解。

1 管網數學模型

天然氣管網由管道、閥門、壓縮機等元件組成，因此描述氣體流經這些元件行為的數學模型就組成天然氣管網的數學模型。

1.1 管道元件數學模型

天然氣管道的數學模型由連續性方程、動量方程和能量方程共同組成。其中，連續性方程和動量方程被稱為管道的水力方程，能量方程則被稱為管道的熱力方程[14]。采用質量流量m、壓力p、溫度T作為描述管道內部水力和熱力過程基礎變量[15]的方程形式如下式所示：

連續性方程：

動量方程：

能量方程：

式中：p為管道內氣體壓力，Pa；t為時間變量，s；T為管道內氣體的溫度，K；A為管道橫截面積，m2；ρ為天然氣密度，kg/m3；m為管道內氣體的質量流量，kg/s；z為軸向空間變量，m；λ為水力摩阻系數；g為重力加速度，m/s2；θ為管道傾斜角，rad；cv為氣體的定容比熱容，J/（kg·K）；u為管道內氣體流速，m/s；D為管道內徑，m；Kt為管道內氣體與環境之間的總傳熱系數，W/（m2·K）；Tg為管道周圍土壤的溫度，K。

其中，質量流量m和壓力p稱為水力變量，溫度T稱為熱力變量。

為了加快管網仿真過程求解速度，可以對上述方程進行線性化處理。線性化的統一表達式[16]如下：

式中各參數的具體形式如表1所示。

表1 管道數學模型參數

采用Implicit Cell Centered Method[17]方法對線性化后的水力方程進行離散，整理后可以得到如下方程：

線性化熱力方程的時間項采用向前差分格式，對流項采用一階迎風格式進行離散，離散結果整理后如下：

式中：UPi、CEi、DWi為根據上一時層計算得到的熱力系數。

1.2 非管道元件數學模型

天然氣管網元件可以分為管道元件和非管道元件，其中非管道元件包括壓縮機、閥門、氣源、節點等。非管道元件具體的數學模型可參考相關文獻[18]。

2 管網求解策略

為加快天然氣水熱力仿真的速度，采用水熱力解耦[19]的求解策略。如圖1所示，該方法將天然氣管網的水熱力過程解耦為水力過程和熱力過程再分別求解，能夠在保證仿真精度的前提下提高仿真速度。

圖1 天然氣管網水熱力解耦算法

2.1 分層算法計算簡介

如圖2所示，直接求解算法是直接采用LU分解算法對整個天然氣管網的水力（熱力）方程組進行求解。由于天然氣管網元件眾多且管道離散節點較多，該方程組規模較大，直接求解算法較為耗時。

圖2 分層算法原理

針對這一問題提出了一種分層算法。該算法主要針對不同的子管網集合進行計算，因此在應用分層算法之前需要將管網按照一定的標準分解為若干子管網的集合。對于常見的燃氣管網或者長輸天然氣管網，可以按照調壓站（調壓閥）或增壓站（壓縮機）所在的位置對管網進行分解操作，從而可以形成不同層級的子管網。需要指出的是，本文提出的分層算法包括圖2 的步驟2 和步驟3，但是步驟1并不屬于分層算法的內容。

將管網分為若干個子管網之后，將其組織為層級結構，之后便可以采用分層算法，開始通過LU 分解方法對各個子管網進行單獨求解，求解過程中采用CPU 并行技術[20-21]將同一層級的子管網同時進行計算，從而加快仿真求解速度。

2.2 分層算法求解水熱力過程

圖3 代表由a層b個子管網通過拓撲連接形成的天然氣管網。子管網內部元件之間的拓撲結構并不會影響分層算法的應用，因此圖3并未顯示各個子管網內部的拓撲結構。管網分層操作使得同一個物理節點同時屬于兩個子管網，通過虛線來表示這種子管網間的拓撲連接。

圖3 具有多層結構的天然氣管網

如圖4 所示，某層的子管網α 與m個子管網連接，其中有n個頂層子管網和m-n個底層子管網。需要說明的是，頂層子管網和底層子管網都是相對本層子管網而言的。

圖4 與外界連接的子管網

假設圖4中的子管網α內部有Np個管道、Nv個閥門、Nc個壓縮機、Ne個氣源以及Nd個節點。在水熱力求解過程中，為了保證管道內部變量的求解精度，通常將管道離散為較多的管段。此處為了簡化說明過程，子管網內部的管道沒有進行離散操作，但這并不會影響分層算法的具體實施流程。

管道、閥門和壓縮機都具有起點和終點兩個端點，氣源只有一個端點，而每個端點上都具有壓力和流量兩個水力變量，因此子管網α內部管道元件、閥門元件、壓縮機元件、氣源元件端點處的水力變量數目共計4Np+4Nv+4Nc+2Ne。

管道、閥門和壓縮機都具有2 個特性方程，而氣源具有1 個特性方程為已知的邊界條件值，因此，子管網內部管道水力方程，閥門、壓縮機和氣源的水力特性方程共計2Np+2Nv+2Nc+Ne個。根據圖論知識可知，與m個外界子管網連接的子管網內部節點上流量平衡方程和壓力平衡方程數目為2Np+2Nv+2Nc+Ne-m。至此，方程數目為4Np+4Nv+4Nc+2Ne-m。

子管網α內部水力變量數目較水力方程數目多m，方程不封閉。為了使水力方程組達到封閉的條件，需要補充m個方程。

如圖4所示，本層子管網與m-n個底層子管網連接。這些底層子管網向子管網α傳入m-n個補充方程，子管網中水力方程組數目變為4Np+4Nv+4Nc+2Ne-n。該補充方程為連接點上兩個水力變量之間的關系式，即連接點上質量流量m和壓力p之間的關系式。

子管網α 通過n個連接點與頂層子管網連接，選用該n個連接點的質量流量m（或壓力p）作為自由變量v，可將管網內部所有水力變量表示為關于該n個自由變量（v1，v2，…，vn）的關系式。

令[v1，v2，…，vn]=[1，0，…，0，…，0 ]，代表關于自由變量（v1，v2，…，vn）的n個方程，其中v1= 1，v2～n= 0。

將該n個方程作為4Np+4Nv+4Nc+2Ne-n個方程的補充方程，通過求解這些方程組可得解向量D1。求解方程組可以采用LU 分解等常見的方程組求解算法。

將關于自由變量（v1，v2，…，vn）的n個方程中的各個變量依次等于1，重復上述過程，于是得到一系列的解向量D2，D3，Dn。

最后令[v1,v2,…,vk]=[0,0,…,0,…,0 ]，求解得到Dn+1。

于是可得子管網α內部所有水力變量關于選定的n個自由變量（v1，v2，…，vn）表達式：

根據上式可得連接點上選定的自由變量v與非自由變量之間關系式：

上式中c表示與子管網α 連接的子管網編號。若自由變量v代表壓力p，非自由變量vˉ則代表流量m；若自由變量v代表流量m，非自由變量vˉ則代表壓力p。

將上述n個關系式作為補充方程傳入頂層子管網中，正如之前由底層子管網傳上來的m-n個補充方程。

如圖5所示，底層的子管網不斷向頂層的子管網傳遞補充方程，直至最頂層的子管網。最頂層的子管網待底層子管網傳入補充方程后，水力方程組閉合。求解該方程組便可得到最頂層的子管網內部所有變量的值，包括最頂層的子管網與其底層子管網連接點上自由變量的值。值得一提的是，分層算法使得每層之間的各個子網的求解互相解耦，因此可利用并行計算等方式實現計算速度的提升。

圖5 分層算法示意

將這些自由變量的值傳入底層子管網，根據前述子管網內部計算得到關系式，得到子管網內部所有變量的值。不斷重復該過程至最底層管網，從而得到整個管網水力變量值，完成管網水力過程的求解。下面以一個具有詳細拓撲結構的天然氣管網為例展示分層算法求解水力仿真的過程。

如圖6 所示的天然氣管網具有9 個氣源，S1 為進氣氣源（源點），S2～S9 為出氣氣源（匯點），其中進氣氣源可以代表儲氣庫或上游管網來氣。該天然氣管網通過閥門V1～V6 調節壓力。實際管網中存在多種環狀管網，不能一一列舉，因此圖中只采用三角形環狀管網，但這并不意味著分層算法只能應用于三角形環狀結構的管網。

圖6 具有多個調壓閥的天然氣管網

按照調壓站（閥）的位置將管網分為如圖7所示的具有3 個不同層級的子管網，其中最底層有4個子管網（3-1，3-2，3-3，3-4），中間層有2 個子管網（2-1，2-2），最頂層有1個子管網（1-1）。圖中紅色虛線代表分層之后不同子管網之間的拓撲連接關系。由于管網在閥門處進行分層，因此選用分界處的閥門變量作為自由變量。

圖7 分層后的天然氣管網拓撲圖

子管網3-1 共有4 條管道，1 個閥門和2 個匯點。依舊假設管道沒有離散，只有起點和終點兩個節點。子管網3-1內部的水力變量和方程數目如表2所示。

表2 子管網3-1水力變量與方程數目

子管網3-1 與外界有1 個連接點，從表2 可知其水力變量數目較方程數目多1。為了使方程組閉合，需補充1 個關于自由變量p5的方程p5=1，并且其余方程右邊的常數項設置為0。采用LU 分解算法求解閉合的方程組可得到解向量A1。同樣，在原先缺1 個方程的方程組基礎上再添加p5=0 的方程，其余方程右邊的常數項也設置為0，采用LU分解算法求解閉合的方程組可以得到解向量B1。因此可以得到如下關系式：

其中，A51和B51代表解向量A1和B1中變量m5對應的系數。

子管網3-1 和3-2 將分別將表達式m5=p5A51+B51和m9=p9A92+B92傳遞給頂層子管網2-1。子管網2-1 與外界有3 個連接點，因此其需要3 個方程才能達到方程組閉合的條件，底層子管網已經傳遞了2 個方程，因此其只需要補充1 個方程即可。接著重復子管網3-1的求解流程即可。

表3 是不同層的子管網應用分層算法的過程中各個連接點選取的自由變量以及產生的一些方程。

表3 不同層的子管網應用分層算法時自由變量與方程

采用分層算法對天然氣管網進行熱力求解與上述水力求解類似，此處不再贅述。

3 案例分析

某天然氣管網具有如圖8 所示的復雜拓撲結構，該管網具有105 條管道（見表4）、3 臺壓縮機、7個閥門以及32個氣源。

圖8 具有復雜拓撲結構的某天然氣管網

表4 管道長度

所有管道的規格均為D720 mm×10 mm，出氣氣源（S1～S14，S16～S32）都采用流量邊界（見表5），進氣氣源（S15）采用壓力邊界（3 MPa），壓縮機（C1～C3）具有相同的壓比-流量曲線（見表6）。

表5 各個出氣氣源（匯點）的流量

表6 壓縮機的流量壓比數據

管網t=0時刻已經處于穩態。當t=2 h時，閥門V1關閉，管網停止向出氣氣源S2供氣。

按照管網中調壓閥V2、V3、V6、V5和增壓站C2、C3 的位置將管網分為如圖9 所示的不同層級子管網集合。

圖9 分區后的天然氣管網拓撲結構

3.1 準確性對比

采用分層算法對上述管網算例進行求解。由于停止向氣源S2 供氣之后，管道P4、P5、P6 及氣源S1、S3、S4 這些元件離關閉的閥門V1 較近，受到的影響最為明顯，因此選擇將管網停止向氣源S2 供氣前2 h 和后24 h 總共26 h 的這些元件的壓力或流量進行對比。同時為了保證對比的充分性，還將對比管網再次達到穩態后采用分層算法與不采用分層算法（即直接求解算法）計算的不同位置處管道（P7、P8、P9、P10、P11、P12）、氣源（S5、S9、S19、S21、S26）、壓縮機（C1、C3、C4）的流量、壓力以及溫度。

需要說明的是，兩種算法對管網仿真過程中出現的一些小方程組都采用LU 分解法進行求解。直接求解算法與分層算法的水力和熱力仿真時間步長為20 s，水力和熱力計算的空間步長都為1 km。

如圖10～圖14 所示，無論是關閥事件發生之后一段時間內閥門附近受到影響最大元件的流量、壓力和溫度對比情況，還是關閉閥門之后天然氣管網再次達到穩態時不同管道、氣源和壓縮機的流量、壓力和溫度對比情況都表明，分層算法計算得到計算結果與直接求解的計算結果完全一致，這說明分層算法具有和直接求解算法同等的計算精度，分層算法的正確性和準確性得到驗證。

圖10 管道P4、P5、P6起點流量隨時間變化情況

圖11 管道P4、P5、P6起點壓力隨時間變化情況

圖12 管道P4、P5、P6起點溫度隨時間變化情況

圖13 氣源S1、S3、S4壓力隨時間變化情況

圖14 氣源S1、S3、S4溫度隨時間變化情況

3.2 求解速度對比

對同在一層的各個子網，可以利用CPU 并行計算實現子管網計算過程的加速。本實驗基于AMD Ryzen 72700x 的CPU 平臺進行測試，表7 是分層算法與直接求解算法在不同管道總長度的管網算例上仿真耗時的對比情況。

表7 兩種算法仿真耗時對比

由表7 可知，在不同規模管網的仿真過程中，分層算法都具有較快的仿真速度，加速比平均在2倍左右。值得一提的是，本案例中僅對比了同樣拓撲結構不同管網規模下的加速效果，可知對其他拓撲結構的管網，當每層管網的獨立子網數量越多時加速效果會更加明顯。

4 結論

針對大規模天然氣管網水熱力模型計算時直接求解算法求解速度慢的問題，提出了一種將管網分層的算法。該算法將管網分為若干層的子管網，通過自下而上的補充方程傳遞和自上而下的值傳遞過程完成管網方程組的求解。通過天然氣管網向匯點供氣的算例來說明分層算法應用于管網水熱力仿真過程的計算精度和速度。算例的計算結果表明：本算法的應用不會改變模型求解結果的精度，能夠確?；A模型的準確性；本算法可以結合并行計算實現計算加速，對本案例中不同規模的天然氣管網仿真計算，平均可達2倍左右的加速比。