何冬輝

(遼寧東科電力有限公司，沈陽 110006)

0 引 言

天然氣壓縮因子，又稱天然氣偏差系數(shù)，指在環(huán)境溫度、壓力下，真實氣體在相同條件下所占體積與相同量理想氣體所占體積的比值[1-2]。該參數(shù)是天然氣流量計算的最關(guān)鍵中間物性參數(shù)之一，其計算精度直接影響天然氣流量的計算精度。近年來，隨著燃?xì)廨啓C(jī)的廣泛應(yīng)用，天然氣常在非理想氣體工況下(壓力不大于1 MPa，溫度范圍不超過10～20 ℃)傳輸。在非理想氣體工況下，天然氣的流量計算需要在理想氣體狀態(tài)方程的基礎(chǔ)上引入壓縮因子。計算機(jī)高新技術(shù)的發(fā)展極大地提高了壓縮因子計算的精度和效率，其中基于AGA8-92DC模型的計算方法得到了國際標(biāo)準(zhǔn)化組織天然氣技術(shù)委員會的推薦。AGA8-92DC模型以狀態(tài)方程為基礎(chǔ)，因其計算精度較其他方法高而得到了廣泛應(yīng)用。該文以AGA8-92DC為計算模型，優(yōu)化壓縮因子迭代算法，采用面向?qū)ο缶幊痰腃#語言，設(shè)計并實現(xiàn)了天然氣壓縮因子計算軟件開發(fā)。

1 AGA8-92DC計算模型

1.1 計算模型介紹

美國標(biāo)準(zhǔn)學(xué)會以孔板流量計標(biāo)準(zhǔn)AGA No.3報告為基礎(chǔ)，對天然氣壓縮因子開展了深度研究，并發(fā)表了基于狀態(tài)方程的壓縮因子計算方法的AGA No.8報告，即AGA8-92DC方程。該方程的出發(fā)點是：任何天然氣容量性質(zhì)均可由一組合適的可測定物性值來表征和計算[1]。通過分析試驗氣體摩爾成分超過0.000 05的各種氣體組分，輸入實際工況下的壓力和溫度來計算壓縮因子。該方程的實際應(yīng)用條件為：輸氣管內(nèi)氣體絕對壓力為0～12 MPa，溫度范圍為263～338 K，計算結(jié)果不確定度約為0.1%。AGA No.8號報告中，以氣體組分成分為基礎(chǔ)進(jìn)行計算的方程如式(1)所示，以氣體物理特性數(shù)據(jù)為基礎(chǔ)進(jìn)行計算的方程如式(2)所示，前者屬于外延式二維維里方程，后者屬于三維維里方程。

(1)

(2)

1.2 計算過程

1.2.1 已知條件

該計算模型通過分析取樣氣體詳細(xì)的摩爾分?jǐn)?shù)組成，并根據(jù)實際工況壓力、溫度來計算氣體壓縮因子。

已知參數(shù)包括氣體實際工況下的氣體摩爾組分?jǐn)?shù)N、絕對壓力p和熱力學(xué)溫度T。

氣體摩爾組成以摩爾百分?jǐn)?shù)表示各氣體組分，氣體包括CO2、N2、H2、CO、CH4、C2H6、C3H8、i-C4H10、n-C4H10、i-C5H10、n-C5H12、n-C6H14、n-C7H16、n-C8H18。

氣體各組分的摩爾分?jǐn)?shù)分別表示為Xi，i=1～N。

1.2.2 計算步驟

1)計算第二維里系數(shù)B。

(3)

(4)

(WiWj+1-wn)wn

二元參數(shù)Eij和Gij，由以下兩式計算得出：

(5)

(6)

(7)

根據(jù)式(8)～(11)分別求解混合方程，計算混合物參數(shù)U，G，Q，F(xiàn)。

(8)

(9)

(10)

(11)

3)計算混合物體積參數(shù)K。

(12)

4)計算對比密度。

摩爾密度為：

ρm=p/(ZRT)

(13)

式中：T為熱力學(xué)溫度，K；R為摩爾氣體常數(shù)；p為絕對壓力，MPa。

對比密度ρr與摩爾密度ρm相關(guān)：

ρr=K3ρm

(14)

5)根據(jù)以上計算系數(shù)，代入式(1)所示的AGA8-92DC方程，用迭代方法求解天然氣壓縮因子。

1.3 算法優(yōu)化

根據(jù)利用AGA8-92DC模型計算天然氣壓縮因子的過程可看出，運算需要多次運用連加、二重求和甚至三重求和等數(shù)學(xué)方法，依次計算58種物質(zhì)的變量狀態(tài)方程參數(shù)、21種識別組分的特征參數(shù)、21種識別組分的二元交互作用參數(shù)組[2]，求解方程時需多次進(jìn)行二次插值和迭代計算。計算過程不僅復(fù)雜，而且計算誤差大、結(jié)果不準(zhǔn)確。為解決上述問題，對迭代算法和流程進(jìn)行全面優(yōu)化，用弦截法代替牛頓迭代法。弦截法是以牛頓迭代法為基礎(chǔ)，利用差商代替牛頓法中的導(dǎo)數(shù)值，通過迭代實現(xiàn)變種及優(yōu)化的算法，不僅收斂速度快，還避免了復(fù)雜計算函數(shù)導(dǎo)數(shù)值問題。弦截法的核心思想是用切線斜率近似割線斜率[3]，與牛頓迭代法相比，其計算量少且條理清晰，便于算法程序開發(fā)。

根據(jù)式(3)～(12)，計算各種輔助參數(shù)后，通過迭代計算摩爾密度以得到壓力p2，根據(jù)標(biāo)準(zhǔn)中定義將計算流程歸納如下：

1)初始化摩爾迭代計算初值，p1、p2；

2)計算摩爾密度ρn和低階摩爾密度ρn-1；

3)當(dāng)計算數(shù)值與初始值的誤差過大時，迭代計算新的密度ρn+1；

4)重新計算最新密度，然后重復(fù)地第2)至4)步，直到計算所得壓力p計算與實測壓力p輸入之差的絕對值小于10-6。

當(dāng)計算所得壓力與測量壓力之間的誤差小于一定范圍后，即得到壓縮因子的計算值，根據(jù)式(1)進(jìn)行壓縮因子計算，其計算流程如圖1所示。采用AGA8-92DC標(biāo)準(zhǔn)計算天然氣壓縮因子，算法流程主要分為4個步驟：調(diào)用已經(jīng)驗證過的待計算數(shù)據(jù)并進(jìn)行初始化操作、輔助參數(shù)計算、迭代計算、壓縮因子和密度計算。壓縮因子迭代計算流程中定義狀態(tài)參數(shù)函數(shù)、特征參數(shù)函數(shù)、二元交互作用參數(shù)函數(shù)三組結(jié)構(gòu)，分別對應(yīng)計算標(biāo)準(zhǔn)中變量狀態(tài)方程參數(shù)、特征參數(shù)、二元交互作用參數(shù)組。壓縮因子計算初始化過程需校驗相關(guān)計算參數(shù)數(shù)值的合法性，隨后根據(jù)AGA8-92DC標(biāo)準(zhǔn)以鍵值數(shù)據(jù)形式格式化計算參數(shù)寄存器以備調(diào)用，其詳細(xì)工作流程如圖1所示。

圖1 迭代計算流程Fig.1 Flow chart of iterative calculation

壓縮因子迭代計算占用了該算法的主要計算資源，因此消減壓縮因子迭代計算可顯著提高算法的計算效率。主要方式是優(yōu)化壓縮因子算法中迭代計算流程，通過快速迭代得出滿足|p2-p1|<10-6條件的摩爾密度ρm，使壓縮因子計算算法性能得以提高，其中，p1為實測壓力，p2為計算壓力，由式(15)計算得出。

p2=ρnTRZ

(15)

根據(jù)(15)可知，p2是1個關(guān)于ρn的函數(shù)，可表達(dá)為式(16)：

p2=p(ρn)

(16)

設(shè)f(ρn)是一個關(guān)于ρn的函數(shù)，且滿足式(17)：

f(ρn)=P(ρn)-p1

(17)

根據(jù)式(17)可以看出，f(ρn)是1個關(guān)于ρn的單調(diào)函數(shù)，當(dāng)P(ρ0)=p1時，f(ρ0)=0，那么P(ρn)=p1就是f(ρ0)的根，即式(17)可使用弦截法求解。當(dāng)P(ρn)=0時，f(ρn)=-p1，-p1<0，那么f(ρn)是一個單調(diào)遞增函數(shù)。結(jié)合式(15)～(17)可知求解|p2-p1|<10-6的過程即等同于求滿足f(ρn)<10-6的pn值。

使用快速收斂的弦截法(見圖2)求式(17)的詳細(xì)步驟如下所述。

圖2 弦截法示意圖Fig.2 schematic diagram of secant method

1)根據(jù)GBT 17747.2—2011標(biāo)準(zhǔn)可知,0≤ρn≤8，根據(jù)ρn范圍，設(shè)P(ρ0)=p1即f(ρ0)=0，選取ρ1=0，ρ2=15為迭代初始值，控制精度為10-6。根據(jù)式(17)計算f(p1)=f(ρn-1)、f(ρ2)=f(ρn)的值，此時n=2,n-1=1，n+1=3。

2)根據(jù)式(17)、(18)，分別求解ρn+1及f(ρn+1)的值。

(18)

3)對f(ρn-1)、f(ρn)、f(ρn+1)的值進(jìn)行判斷。由于f(ρn)是一個單調(diào)遞增函數(shù)，那么f(ρn-1)

①若f(ρn+1)<0且|f(ρn-1)|>|f(ρn+1)|>10-6，則賦值令n=n、n-1=n+1，重復(fù)2)與3)；

②若f(ρn+1)>0且|f(ρn)|>|f(ρn+1)|>10-6，則令n=n+1、n-1=n-1，重復(fù)2)與3)；

③若|f(ρn+1)|<10-6，則停止計算，此時的ρn+1即需要計算的ρm，迭代計算結(jié)束。

2 軟件開發(fā)

2.1 軟件設(shè)計

基于軟件功能需求，采用面向?qū)ο蟪绦虻腃#語言在Visual Studio.NET的WPF框架下對天然氣壓縮因子計算軟件進(jìn)行設(shè)計，實現(xiàn)人機(jī)交互可視化界面、輸入輸出數(shù)據(jù)處理、核心計算動態(tài)鏈接庫(Dynamic Link Library，DLL)等模塊的構(gòu)建。天然氣壓縮因子的計算主要針對實際工況下天然氣的詳細(xì)組分進(jìn)行計算，因此該程序共分為2個部分：人機(jī)交互界面部分和核心計算DLL部分。系統(tǒng)界面部分主要負(fù)責(zé)參數(shù)輸入和校驗、核心計算模塊的調(diào)用和計算結(jié)果的輸出顯示。核心計算DLL部分負(fù)責(zé)對輸入?yún)?shù)進(jìn)行計算。系統(tǒng)的計算功能在用戶界面窗口實現(xiàn)，用戶通過界面窗口與系統(tǒng)交互，界面組件接受數(shù)據(jù)輸入并檢測用戶輸入合法性，輸入非法數(shù)據(jù)將提示用戶。待數(shù)據(jù)輸入完整且正確后，點擊“計算”按鈕，系統(tǒng)調(diào)用核心計算DLL計算模塊，將最終計算結(jié)果綁定到界面組件并輸出顯示。

2.2 軟件實現(xiàn)

考慮到壓力和溫度有多種計量單位，為體現(xiàn)交互界面的直觀友好，系統(tǒng)提供多種單位數(shù)據(jù)輸入模式，用戶可根據(jù)實際應(yīng)用情況選擇。此外，由于天然氣氣體的成分較復(fù)雜，為方便用戶使用，把所有組分參數(shù)輸入設(shè)置在同一個界面上，組分默認(rèn)值都設(shè)為0，用戶可根據(jù)計算需要，填入相應(yīng)的氣體組分。由于每次用戶輸入量較多，系統(tǒng)設(shè)置了將用戶的歷次輸入進(jìn)行自動保存的功能，每次運行時系統(tǒng)會自動調(diào)出最近一次的輸入值，也可找出歷史記錄，同時，用戶可根據(jù)實際輸入情況進(jìn)行修改。

3 實例應(yīng)用分析

以某天然氣的取樣數(shù)據(jù)為例(見表1)，將表中數(shù)據(jù)輸入軟件界面，點擊“計算”按鈕，便可得出試驗工況和基準(zhǔn)工況下的壓縮因子和氣體真實密度，結(jié)果如圖3所示。

表1 天然氣取樣組分Table 1 Components of natural gassampling

圖3 實例演算界面Fig.3 Interface of example calculation

以GB/T17747.2附錄C中表C1中的1～6號氣體取樣數(shù)據(jù)為例[4]，應(yīng)用該軟件對不同工況下多樣本取樣數(shù)據(jù)進(jìn)行計算，計算結(jié)果如表2所示。與附錄C中表C2所給結(jié)果進(jìn)行對比，結(jié)果表明，該軟件計算結(jié)果平均誤差幾乎為0，完全能滿足工程實際應(yīng)用的需要。

表2 計算結(jié)果與附錄C中表C2結(jié)果的對比分析Table 2 Comparative analysis of calculation results and table C2 in appendix C

4 結(jié) 語

1)基于AGA8-92DC計算模型，將求解步驟進(jìn)行模塊分解，并探討了算法流程，用弦截法對牛頓迭代算法進(jìn)行優(yōu)化，全面優(yōu)化天然氣壓縮因子算法流程，提高了計算效率。

2)根據(jù)數(shù)學(xué)模型和算法流程，基于Visual Studio.NET開發(fā)平臺，采用C#語言設(shè)計并實現(xiàn)了壓縮因子計算軟件開發(fā)，將核心計算模塊封裝成動態(tài)鏈接庫形式，便于軟件直接調(diào)用，同時將不變量進(jìn)行緩存并可直接調(diào)用，具有良好的擴(kuò)展性及安全性，便于二次開發(fā)，提升了開發(fā)效率和計算精度。

3)采用GB/T 17747.2附錄C中的數(shù)據(jù)對所開發(fā)的軟件進(jìn)行驗證計算，與GB/T 17747.2表C2結(jié)果進(jìn)行對比分析，結(jié)果表明，該軟件計算平均相對誤差幾乎為0，能滿足工程實際應(yīng)用要求。