基于流量分配優(yōu)化的天然氣管網(wǎng)模擬計算方法

林 敏

(中國石油化工股份有限公司天然氣榆濟管道分公司)

0 引 言

我國已成為天然氣市場的主要國家，截至2022年底，我國已成為天然氣進口大國。天然氣燃燒產(chǎn)生的熱值是其價值的主要標志，相同體積不同氣源的天然氣燃燒產(chǎn)生的熱值差別很大，而我國天然氣的體積熱值最小為33.9 MJ/m3，最大為45 MJ/m3，因此不同熱值的天然氣以同樣的體積價格進行交易是不合理的。我國天然氣計量、計價和貿易交接仍以體積計量為主，在管網(wǎng)互聯(lián)互通、氣質復雜多元的情況下，體積計量已不適應目前天然氣行業(yè)發(fā)展的需要。為此，國家相關部門和天然氣生產(chǎn)、輸送及銷售企業(yè)均在努力推進天然氣能量計量方式的實施。

在能量計量模式下，根據(jù)GB/T 22723—2008《天然氣能量的測定》[1]，精確地模擬天熱氣管網(wǎng)中氣體的流動狀態(tài)，是節(jié)約設備投資和計量站現(xiàn)場改造、施工費用的必要前提。同時GB/T 18603—2014《天然氣計量系統(tǒng)技術要求》中也規(guī)定了B級和C級計量站需要采用賦值方式獲取組成和發(fā)熱量數(shù)據(jù)的要求[2]。

本文將水力方程、能量方程、組分守恒方程及真實氣體狀態(tài)方程耦合，結合流量分配算法、組分混合算法，計算天然氣在管道中任意位置的流動參數(shù)與氣質組成，依據(jù)GB/T 11062—2020《天然氣發(fā)熱量、密度、相對密度和沃泊指數(shù)的計算方法》，實現(xiàn)標準狀態(tài)下天然氣發(fā)熱量的精確計算[3]，以期為能量計量的實施提供幫助，節(jié)約設備投資和計量站改造費用。

1 算法與數(shù)學模型

天然氣在管網(wǎng)內的流動狀態(tài)模擬以水力計算為基礎，通過求解水力方程組獲得管道內氣體的壓力和流速，使用能量守恒方程和狀態(tài)方程求解管道內氣體的密度及溫度，依據(jù)組分輸運方程求解天然氣組成，結合流量分配與組分混合算法，計算天然氣的發(fā)熱量。

1.1 控制方程組

天然氣流動狀態(tài)水力方程組包含連續(xù)性方程和動量守恒方程：

(1)

能量守恒方程為：

(2)

式中：ρ為密度，kg/m3；ω為流速，m/s；t為時間步長，s；x為空間步長，m；p為壓力，Pa；g為重力加速度，m/s2；θ為管道傾角，rad；λ為管道摩擦因數(shù)；D為管道內徑，m；A為管道截面積，m2；T為管內天然氣溫度，K；h為焓，J；T0為環(huán)境溫度，K；K為綜合換熱系數(shù)，W/(m2·K)。

方程(1)中的管道摩擦因數(shù)λ由 Colebrook-White 方程計算[4]，其表達式如下：

(3)

式中：e為管壁表面粗糙度，m；Re為雷諾數(shù)。

在天然氣運輸過程中，常有多個氣源混合的情況，混合后天然氣管網(wǎng)中的氣質組成情況可通過數(shù)值計算方法獲得，計算過程中需要耦合組分輸運方程：

(4)

式中：c為摩爾分數(shù)；Dx為軸向擴散系數(shù)，m2/s；Dr為徑向擴散系數(shù)，m2/s；r為半徑，m；Rc為源項，s-1。

在求解過程中，首先根據(jù)水力方程、能量守恒方程和GERG-2008 方程迭代計算出天然氣流動參數(shù)；然后利用組分輸運方程求解管道組分變化，并更新熱力學參數(shù)，求解節(jié)點的發(fā)熱量。

由于軸向擴散損失和徑向擴散損失相較于整個管網(wǎng)系統(tǒng)的損失較小，在計算時將擴散項和源項忽略，以節(jié)省計算資源，提升計算效率。此時組分輸運方程為[5-6]：

(5)

因為計算介質為天然氣，所以需要考慮氣體狀態(tài)方程才能使求解方程組封閉。在ISO標準和技術應用中，AGA-8 和 GERG-2008是2種常用的計算方程[7]。其中，GERG-2008方程是德國波鴻大學熱力學研究所提出，適用于求解天然氣和其他混合物在較寬壓力范圍內的體積性質和熱性質，改善了AGA-8方程對混合氣體性質的計算性能[8-10]。GERG-2008 方程可以通過無量綱亥姆霍茲自由能及其導數(shù)顯式表達天然氣的所有熱力學性質，天然氣的亥姆霍茲自由能計算與組分有關[11]?？梢允褂脽o量綱亥姆霍茲自由能及其導數(shù)簡化模擬方程組，并耦合氣體組分輸運方程實現(xiàn)天然氣流動速度、氣體流量分配及氣體組分混合結果的精確模擬。亥姆霍茲自由能公式如下：

(6)

無量綱亥姆霍茲自由能可以表示為：

(7)

(8)

ρr和Tr是混合物無量綱密度和無量綱溫度的函數(shù)，只取決于混合物的摩爾組成。對于能量守恒方程，使用無量綱亥姆霍茲自由能對內能和焓進行替換并進行處理，得到：

(9)

方程(9)中，無量綱亥姆霍茲自由能α的上標r表示剩余無量綱亥姆霍茲自由能，下標δ、τ表示無量綱亥姆霍茲自由能對無量綱密度及無量綱時間的偏導數(shù)。

由于無量綱亥姆霍茲自由能及其偏導數(shù)僅是氣體密度(無量綱密度)、溫度(無量綱溫度)及組分的函數(shù)，且無量綱亥姆霍茲自由能及其導數(shù)在上述方程中顯式表達，即在已知氣質組成的情況下可以根據(jù)初始條件直接計算。所以采用GERG-2008方程作為狀態(tài)方程時，不僅提升了混合氣體性質的計算性能，還減少了不必要的微分或積分表達式。由上述分析可知，能量守恒方程可以表示為無量綱密度、無量綱溫度、天然氣流速和無量綱亥姆霍茲自由能及其偏導數(shù)的函數(shù)。

天然氣管網(wǎng)中氣體的流動狀態(tài)實時變化，但單一管道的材料、狀態(tài)及使用年限等可以認為一致。對于長輸管道來說，管道長度足夠長，且管段曲率半徑遠大于直徑，因此，天然氣在管道中的流動可以視為一維非穩(wěn)態(tài)流動，管道的離散結構如圖1所示。采用如圖2所示的分段線性型線作為函數(shù)及其導數(shù)對時間及空間的局部型線，對守恒型方程組積分，得到離散后的控制方程，如式(10)所示。

圖1 管網(wǎng)空間離散示意圖Fig.1 Schematic diagram of spatial discretization of the pipeline network

圖2 分段線性型線Fig.2 Segmented linear lines

采用基于壓力的有限體積法進行計算，計算流程如圖3所示。

圖3 計算流程圖Fig.3 Calculation workflow

(10)

1.2 流量分配算法

在流量分配過程中，由于管道連接點處天然氣流動狀態(tài)和損失的復雜性，難以通過控制方程準確求解下游管道初始參數(shù)。在實際工況下，下游管道往往安裝壓力和溫度計量設備。所以可以在將儀表測量值作為限制條件的前提下，假設下游管道初始壓力、溫度和流速，進行迭代求解。這個過程可以視為多目標優(yōu)化問題，其中決策變量為下游管道初始壓力、溫度、流速，優(yōu)化目標如下：

(11)

式中：e1、e2為誤差；p為初始壓力，Pa；T為溫度，K；下標c代表計算值；下標r代表實際值。且2個誤差均小于相應的誤差要求；限制條件則需要對整個天然氣流動模擬過程進行處理分析。

在優(yōu)化過程中，需要多次使用前文章模型進行計算，這會耗費大量時間。代理模型可以將通過模擬獲得的樣本點進行擬合，從而得到每個目標函數(shù)與輸入變量間的近似預測模型，無需進行大量的數(shù)值求解工作[12]。本文通過抽樣方法選取樣本點后，建立Kriging代理預測模型，使用NSGA-Ⅱ遺傳算法進行優(yōu)化求解。

Kriging代理模型的建立需要在初始參數(shù)組成的樣本空間內進行抽樣，拉丁超立方抽樣(LHS)方法是一種隨機多維分層抽樣方法，根據(jù)初始參數(shù)的取值范圍，將樣本空間等分為N個相互獨立的子區(qū)域，并在每個子區(qū)域中進行等概率抽樣，得到具有良好空間均勻性和填充性的樣本集，因此本文采用LHS方法進行抽樣操作。

為了滿足Kriging代理模型的計算精度，初始樣本的數(shù)量不少于決策變量的10倍，本文流量分配優(yōu)化的決策變量為3個，為保證代理模型的通用性，環(huán)境溫度也應作為初始變量，因此初始樣本點的數(shù)量設計為50個[13-14]。

在代理模型建成后，采用R2誤差分析法驗證代理模型精度。R2誤差分析法反映了預測值與實測值的偏差程度，R2越接近1，代理預測模型的全局近似精度越高。其計算方法如下：

(12)

NSGA-Ⅱ作為一種采用快速非支配排序過程、精英保留策略的排序遺傳算法，與傳統(tǒng)遺傳算法相比，其算法運行效率高，收斂速度快，適用于3個及以下的多個目標優(yōu)化設計[15]。其主要計算流程為[16]：

(1)隨機產(chǎn)生個體數(shù)為N的初始種群Fg，按照非支配關系排序后進行交叉、變異操作，產(chǎn)生子代種群Sg。

(2)將父代種群Fg與子代種群Sg混合，按照非支配關系排序后從最優(yōu)個體開始，選取N個個體組成新的初始種群Fg+1。

(3)對種群Fg+1進行交叉、變異操作后，重復流程(1)和(2)過程，直到達到迭代步數(shù)。

基于上述方法，可以在最終進化結果中，構建出滿足下式的 Pareto前沿：

(13)

對Pareto前沿中的非劣解的目標函數(shù)值歸一化，選取距離理想結果最近的非劣解作為最優(yōu)解[17]。

優(yōu)化設計得到結果的優(yōu)劣很大程度上取決于代理模型的預測精度，合理的加點準則可以提高模型精度。把最優(yōu)解對應的決策變量代入控制方程進行模擬計算，將獲得的模擬結果與溫度計、壓力表示值進行對比，若誤差大于規(guī)定，說明代理模型在該決策變量取值下的預測結果不滿足要求。此時，將該點加入初始樣本點對代理模型進行更新，重新進行優(yōu)化計算，直到誤差滿足要求。將優(yōu)化結果作為初始條件，利用邊界條件隨時間的變化規(guī)律，使用控制方程迭代計算每個時刻管道各處的參數(shù)，然后根據(jù)計算誤差來確定是否需要做修正。流量分配計算流程如圖4所示。

圖4 流量分配計算流程圖Fig.4 Workflow of flow distribution calculation

1.3 組分混合算法

在天然氣管網(wǎng)結構中，往往會有多個天然氣氣源。當不同氣源的天然氣匯集到同一管道時，需要對混合后的天然氣氣質重新進行分析。FAN D.等[18]提出了一種管網(wǎng)結點結構的特殊處理方法，即將動量方程單獨進行離散，其結點離散方式如圖5所示。

圖5 結點單元示意圖Fig.5 Schematic diagram of node elements

如圖5a所示，動量方程在動量控制體進行離散，其余方程在圖5b所示的質量控制體內離散。其中，在質量控制體進行離散的方程均可用輸運方程的形式表達。

對于整個結點結構，其控制方程組可以表示為如下形式：

(14)

通過求解控制方程組(14)，可以獲得不同管道天然氣混合后的組分、壓力、流量及密度等參數(shù)。需要注意的是：天然氣管道長度較長，混合氣體在相對較短的距離內達到均勻狀態(tài)，所以這里沒有考慮氣體混合后的均勻性問題。

天然氣管網(wǎng)各處氣體組分含量已知后，使用ISO 20762—2中的方法對壓縮因子、聲速進行計算。依據(jù) GB/T 11062—2020 中的方法計算發(fā)熱量和沃泊指數(shù)。

2 工程實例對比分析

基于前文所述計算方法，編制了計算程序，分別對流量分配、組分追蹤和發(fā)熱量進行求解。以山東管網(wǎng)系統(tǒng)中3個計量站的實際運營數(shù)據(jù)進行驗證，其結構示意如圖6所示。

圖6 計量站示意圖Fig.6 Schematic diagram of the measuring station

在12：00-21：00，計量站A作為輸氣站，計量站B、C作為儲氣站，此時可以用來驗證流量分配計算模型。在1：00-10：00，計量站B、C作為輸氣站，計量站A作為儲氣站，此時可以用來驗證組分混合模型。管道的參數(shù)如表1所示。

表1 管道尺寸參數(shù)Table 1 Pipeline dimensions

2.1 控制方程組驗證

選擇圖6計量站A中的直管道對控制方程的計算精度進行驗證。管道的長度為1 500 m，內徑為585 mm。6個時刻的計算結果如表2所示，計算誤差如圖7所示。

表2 不同時刻參數(shù)Table 2 Parameters at different time

圖7 控制方程組計算誤差Fig.7 Calculation errors of governing equations

由圖7可以看出，使用控制方程組(式(10))計算得到的溫度最大差值為0.22 ℃，壓力最大誤差為0.17%，工況流量最大誤差為0.33%，均滿足GB/T 18603—2014對A級計量站的計量要求。

2.2 流量分配算法驗證

各計量站不同時刻的參數(shù)值如圖8所示，在這段時間內，氣質不發(fā)生變化。

圖8 不同時刻計量站計量數(shù)據(jù)Fig.8 Measuring data of the measuring stations at different time

以管道2為例，使用遺傳算法進行流量分配計算。首先，依據(jù)計量站B的歷史運營數(shù)據(jù)確定代理模型各物理量的取值范圍，采用LHS方法抽取50組初始樣本。將樣本分別代入第1章中的計算模型，利用計算結果分別構建儀表測點壓力、溫度的Kriging代理模型。獲得代理模型后，另取10組初始值進行驗證。經(jīng)驗證，R2值分別為0.97和0.95，說明該代理模型計算精度足夠高。

使用NSGA—Ⅱ遺傳算法進行優(yōu)化求解，其中，種群數(shù)設為200，變異率為0.1，交叉率為0.9，最大迭代次數(shù)設為500。管道2的初始壓力、溫度及流速的取值范圍如下：

0.95pb≤p0≤1.05pb

(15)

0.95Tb≤T0≤1.05Tb

(16)

0.95ωc≤ω0≤1.05ωc

(17)

式中：下標b代表管道1末端的物理量；下標0代表管道2初始時刻的物理量；ωc為以管道1末端溫度、壓力作為管道2初始值，使用二分法計算得到初始流速的估計值。管道2初始時刻的取值范圍如表3所示。

表3 管道2初始參數(shù)取值范圍Table 3 Ranges of initial parameters of Pipeline No.2

由于目標優(yōu)化的最終結果是壓力、溫度誤差取最小，所以只要Pareto最優(yōu)解滿足式(13)，就可以認為對應的決策變量為最優(yōu)。利用第一個時刻的決策變量最優(yōu)值，計算后續(xù)時刻計量儀表安置點的結果，誤差如圖9所示。

圖9 計算誤差圖Fig.9 Calculation errors

由圖9可知，最大溫度差值為0.26 K，最大壓力誤差為0.19%，滿足標準要求，證明該算法可以精確計算流量分配模型。需要注意的是：如果上游管道參數(shù)發(fā)生較大變化，模擬計算結果不一定滿足精度要求，此時會重新抽樣進行遺傳算法計算，從而保證每個時刻的模擬結果在誤差范圍內。

2.3 氣質組成和發(fā)熱量計算驗證

在1：00-10：00，計量站B、C作為輸氣站，計量站A作為儲氣站。此時，管道2、3有著不同的氣質組成，天然氣在管道1發(fā)生組分混合。

表4為不同時刻各計量站天然氣中CH4的實際摩爾分數(shù)。

表4 不同時刻各計量站CH4摩爾分數(shù)Table 4 CH4 molar fraction of each measuring station at different time

圖10為計量站B、C中除CH4外天然氣各主要組分實際摩爾分數(shù)。從表4和圖10中可以看出，隨著時間的變化，不同氣源的天然氣氣質組成會發(fā)生無規(guī)律波動。這是因為天然氣在開采時氣質組成會發(fā)生波動，同時在天然氣管網(wǎng)系統(tǒng)的工作過程中，管道中天然氣流向的改變、閥門的啟閉等也會對管道中的氣質組成造成影響。

圖10 不同時刻計量站B、C各組分摩爾分數(shù)Fig.10 Molar fractions of each component in the measuring stations B and C at different time

使用前文所述方法，根據(jù)現(xiàn)場實際數(shù)據(jù)，將模擬計算的時間步長取為1 h，進行氣質跟蹤模擬計算。圖11為計量站A 中天然氣主成分(甲烷、乙烷、丙烷)及其他成分的計量與模擬數(shù)據(jù)。

圖11 不同時刻計量站A各組分摩爾分數(shù)模擬值與計量值Fig.11 Simulated and measured molar fractions of each component in the measuring station A at different time

由圖11可以看出：甲烷含量在5：00、6：00時發(fā)生較大波動，計量誤差在5：00時達到最大，為0.178%，其余時刻的計量誤差在0.1%以下；乙烷和丙烷含量在5：00時也發(fā)生了明顯的變化，這是因為該時刻計量站C中天然氣正庚烷含量發(fā)生了較大變化，計量站B、C中乙烷含量也有波動；乙烷含量在1：00、2：00時發(fā)生較大變化，計量誤差在2：00時達到最大，為0.698%，這是因為在1：00、2：00時，計量站C中乙烷含量波動變大；除乙烷外氣體成分在7：00、8：00時刻均發(fā)生了波動，且丁烷含量波動較為劇烈，這是因為計量站B、C中正庚烷含量發(fā)生了明顯波動。除此之外，其余組分的計量誤差均在1%以下。

從上述結果可知，可以通過混合前管道氣質的波動程度判斷混合后天然氣氣質會發(fā)生變化，但這種變化無明顯規(guī)律。這是因為天然氣混合是一個十分復雜的過程，與流量、壓力、溫度及管道夾角等參數(shù)有關。

在氣質波動較大時，混合后氣體計量誤差與模擬誤差偏大，這是因為氣質波動大，管道內氣體流動不穩(wěn)定，計量設備的計量結果存在誤差。雖然各組分的計量值與模擬值存在誤差，但是只要最終計算得到的發(fā)熱量誤差在標準規(guī)定的范圍內，該模型的精確度便可以得到驗證。

根據(jù)模型計算得到的狀態(tài)參數(shù)和氣質組成，計算得到標準狀況下不同時刻的發(fā)熱量數(shù)據(jù)，分別將高位體積發(fā)熱量和低位體積發(fā)熱量數(shù)據(jù)繪圖，結果如圖12所示。

圖12 發(fā)熱量計量和模擬誤差Fig.12 Errors between measured and simulated calorific values

從圖12可以看出，發(fā)熱量最大偏差出現(xiàn)在1：00，其最大偏差分別為0.305%和0.313%，均小于GB/T 18603—2014 中規(guī)定的0.5%，滿足A級計量場站精度要求。發(fā)熱量在2：00、3：00、6：00時升高，在4：00、5：00時降低，這與氣質的變化是同步的。

3 結論及認識

(1)提出了一種基于GERG-2008 方程的天然氣發(fā)熱量賦值計算方法。使用壓力適用范圍廣、計算性能好的 GERG-2008 方程作為狀態(tài)方程，采用無量綱亥姆霍茲自由能及其導數(shù)表達熱力學參數(shù)，耦合組分輸運方程，對天然氣進行流動模擬和發(fā)熱量賦值計算，結果表明，使用該方法計算得到的天然氣溫度、流量及壓力的最大誤差(差值)分別為0.22 ℃、0.17%、0.33%。

(2)提出了基于分流管道下游計量數(shù)據(jù)的自適應流量分配方法。通過構建Kriging代理方程簡化計算過程，以溫度、壓力誤差最小為優(yōu)化目標，使用 NSGA-Ⅱ 遺傳算法計算下游管道初始參數(shù)，基于計算結果更新Kriging代理模型，并對整根管道進行自適應模擬計算。經(jīng)驗證，流量分配計算中溫度的最大差值為0.26℃，壓力的最大誤差為0.19%。

(3)采用交錯網(wǎng)格處理管道結點結構，動量方程在動量控制體進行離散，其余方程在質量控制體內離散，進而進行氣體混合計算，從而提升氣質跟蹤計算的精確性。經(jīng)驗證，天然氣組分混合計算中，甲烷的最大誤差為0.178%，乙烷的最大誤差為0.698%，丙烷的最大誤差為0.732%，其余組分誤差均在1%以下；高位和低位體積發(fā)熱量最大偏差分別為0.305%、0.313%，計算結果均滿足GB/T 18603—2014中A級計量站的計量要求。

(4)本文的方法可以為能量計量的實施提供幫助，節(jié)約設備投資和計量站現(xiàn)場改造、施工費用。