林 敏

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

0 引 言

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

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

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

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

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

1.1 控制方程組

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

(1)

能量守恒方程為：

(2)

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

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

(3)

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

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

(4)

式中：c為摩爾分?jǐn)?shù)；Dx為軸向擴(kuò)散系數(shù)，m2/s；Dr為徑向擴(kuò)散系數(shù)，m2/s；r為半徑，m；Rc為源項(xiàng)，s-1。

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

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

(5)

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

(6)

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

(7)

(8)

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

(9)

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

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

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

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

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

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

圖3 計(jì)算流程圖Fig.3 Calculation workflow

(10)

1.2 流量分配算法

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

(11)

式中：e1、e2為誤差；p為初始?jí)毫Γ琍a；T為溫度，K；下標(biāo)c代表計(jì)算值；下標(biāo)r代表實(shí)際值。且2個(gè)誤差均小于相應(yīng)的誤差要求；限制條件則需要對(duì)整個(gè)天然氣流動(dòng)模擬過(guò)程進(jìn)行處理分析。

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

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

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

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

(12)

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

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

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

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

基于上述方法，可以在最終進(jìn)化結(jié)果中，構(gòu)建出滿足下式的 Pareto前沿：

(13)

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

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

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

1.3 組分混合算法

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

圖5 結(jié)點(diǎn)單元示意圖Fig.5 Schematic diagram of node elements

如圖5a所示，動(dòng)量方程在動(dòng)量控制體進(jìn)行離散，其余方程在圖5b所示的質(zhì)量控制體內(nèi)離散。其中，在質(zhì)量控制體進(jìn)行離散的方程均可用輸運(yùn)方程的形式表達(dá)。

對(duì)于整個(gè)結(jié)點(diǎn)結(jié)構(gòu)，其控制方程組可以表示為如下形式：

(14)

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

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

2 工程實(shí)例對(duì)比分析

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

圖6 計(jì)量站示意圖Fig.6 Schematic diagram of the measuring station

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

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

2.1 控制方程組驗(yàn)證

選擇圖6計(jì)量站A中的直管道對(duì)控制方程的計(jì)算精度進(jìn)行驗(yàn)證。管道的長(zhǎng)度為1 500 m，內(nèi)徑為585 mm。6個(gè)時(shí)刻的計(jì)算結(jié)果如表2所示，計(jì)算誤差如圖7所示。

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

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

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

2.2 流量分配算法驗(yàn)證

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

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

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

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

0.95pb≤p0≤1.05pb

(15)

0.95Tb≤T0≤1.05Tb

(16)

0.95ωc≤ω0≤1.05ωc

(17)

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

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

由于目標(biāo)優(yōu)化的最終結(jié)果是壓力、溫度誤差取最小，所以只要Pareto最優(yōu)解滿足式(13)，就可以認(rèn)為對(duì)應(yīng)的決策變量為最優(yōu)。利用第一個(gè)時(shí)刻的決策變量最優(yōu)值，計(jì)算后續(xù)時(shí)刻計(jì)量?jī)x表安置點(diǎn)的結(jié)果，誤差如圖9所示。

圖9 計(jì)算誤差圖Fig.9 Calculation errors

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

2.3 氣質(zhì)組成和發(fā)熱量計(jì)算驗(yàn)證

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

表4為不同時(shí)刻各計(jì)量站天然氣中CH4的實(shí)際摩爾分?jǐn)?shù)。

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

圖10為計(jì)量站B、C中除CH4外天然氣各主要組分實(shí)際摩爾分?jǐn)?shù)。從表4和圖10中可以看出，隨著時(shí)間的變化，不同氣源的天然氣氣質(zhì)組成會(huì)發(fā)生無(wú)規(guī)律波動(dòng)。這是因?yàn)樘烊粴庠陂_(kāi)采時(shí)氣質(zhì)組成會(huì)發(fā)生波動(dòng)，同時(shí)在天然氣管網(wǎng)系統(tǒng)的工作過(guò)程中，管道中天然氣流向的改變、閥門的啟閉等也會(huì)對(duì)管道中的氣質(zhì)組成造成影響。

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

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

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

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

從上述結(jié)果可知，可以通過(guò)混合前管道氣質(zhì)的波動(dòng)程度判斷混合后天然氣氣質(zhì)會(huì)發(fā)生變化，但這種變化無(wú)明顯規(guī)律。這是因?yàn)樘烊粴饣旌鲜且粋€(gè)十分復(fù)雜的過(guò)程，與流量、壓力、溫度及管道夾角等參數(shù)有關(guān)。

在氣質(zhì)波動(dòng)較大時(shí)，混合后氣體計(jì)量誤差與模擬誤差偏大，這是因?yàn)闅赓|(zhì)波動(dòng)大，管道內(nèi)氣體流動(dòng)不穩(wěn)定，計(jì)量設(shè)備的計(jì)量結(jié)果存在誤差。雖然各組分的計(jì)量值與模擬值存在誤差，但是只要最終計(jì)算得到的發(fā)熱量誤差在標(biāo)準(zhǔn)規(guī)定的范圍內(nèi)，該模型的精確度便可以得到驗(yàn)證。

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

圖12 發(fā)熱量計(jì)量和模擬誤差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級(jí)計(jì)量場(chǎng)站精度要求。發(fā)熱量在2：00、3：00、6：00時(shí)升高，在4：00、5：00時(shí)降低，這與氣質(zhì)的變化是同步的。

3 結(jié)論及認(rèn)識(shí)

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

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

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

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