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

1.國家管網(wǎng)集團科學技術(shù)研究總院分公司，河北廊坊 065000

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

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

1 管網(wǎng)數(shù)學模型

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

1.1 管道元件數(shù)學模型

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

連續(xù)性方程：

動量方程：

能量方程：

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

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

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

式中各參數(shù)的具體形式如表1所示。

表1 管道數(shù)學模型參數(shù)

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

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

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

1.2 非管道元件數(shù)學模型

天然氣管網(wǎng)元件可以分為管道元件和非管道元件，其中非管道元件包括壓縮機、閥門、氣源、節(jié)點等。非管道元件具體的數(shù)學模型可參考相關(guān)文獻[18]。

2 管網(wǎng)求解策略

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

圖1 天然氣管網(wǎng)水熱力解耦算法

2.1 分層算法計算簡介

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

圖2 分層算法原理

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

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

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

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

圖3 具有多層結(jié)構(gòu)的天然氣管網(wǎng)

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

圖4 與外界連接的子管網(wǎng)

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

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

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

子管網(wǎng)α內(nèi)部水力變量數(shù)目較水力方程數(shù)目多m，方程不封閉。為了使水力方程組達到封閉的條件，需要補充m個方程。

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

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

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

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

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

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

于是可得子管網(wǎng)α內(nèi)部所有水力變量關(guān)于選定的n個自由變量（v1，v2，…，vn）表達式：

根據(jù)上式可得連接點上選定的自由變量v與非自由變量之間關(guān)系式：

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

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

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

圖5 分層算法示意

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

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

圖6 具有多個調(diào)壓閥的天然氣管網(wǎng)

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

圖7 分層后的天然氣管網(wǎng)拓撲圖

子管網(wǎng)3-1 共有4 條管道，1 個閥門和2 個匯點。依舊假設(shè)管道沒有離散，只有起點和終點兩個節(jié)點。子管網(wǎng)3-1內(nèi)部的水力變量和方程數(shù)目如表2所示。

表2 子管網(wǎng)3-1水力變量與方程數(shù)目

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

其中，A51和B51代表解向量A1和B1中變量m5對應(yīng)的系數(shù)。

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

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

表3 不同層的子管網(wǎng)應(yīng)用分層算法時自由變量與方程

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

3 案例分析

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

圖8 具有復雜拓撲結(jié)構(gòu)的某天然氣管網(wǎng)

表4 管道長度

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

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

表6 壓縮機的流量壓比數(shù)據(jù)

管網(wǎng)t=0時刻已經(jīng)處于穩(wěn)態(tài)。當t=2 h時，閥門V1關(guān)閉，管網(wǎng)停止向出氣氣源S2供氣。

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

圖9 分區(qū)后的天然氣管網(wǎng)拓撲結(jié)構(gòu)

3.1 準確性對比

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

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

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

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

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

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

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

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

3.2 求解速度對比

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

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

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

4 結(jié)論

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