何冬輝

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

0 引 言

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

1 AGA8-92DC計算模型

1.1 計算模型介紹

美國標準學(xué)會以孔板流量計標準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ǔ)進行計算的方程如式(1)所示，以氣體物理特性數(shù)據(jù)為基礎(chǔ)進行計算的方程如式(2)所示，前者屬于外延式二維維里方程，后者屬于三維維里方程。

(1)

(2)

1.2 計算過程

1.2.1 已知條件

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

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

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

氣體各組分的摩爾分數(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]，求解方程時需多次進行二次插值和迭代計算。計算過程不僅復(fù)雜，而且計算誤差大、結(jié)果不準確。為解決上述問題，對迭代算法和流程進行全面優(yōu)化，用弦截法代替牛頓迭代法。弦截法是以牛頓迭代法為基礎(chǔ)，利用差商代替牛頓法中的導(dǎo)數(shù)值，通過迭代實現(xiàn)變種及優(yōu)化的算法，不僅收斂速度快，還避免了復(fù)雜計算函數(shù)導(dǎo)數(shù)值問題。弦截法的核心思想是用切線斜率近似割線斜率[3]，與牛頓迭代法相比，其計算量少且條理清晰，便于算法程序開發(fā)。

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

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

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

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

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

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

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

1)根據(jù)GBT 17747.2—2011標準可知,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)的值進行判斷。由于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框架下對天然氣壓縮因子計算軟件進行設(shè)計，實現(xiàn)人機交互可視化界面、輸入輸出數(shù)據(jù)處理、核心計算動態(tài)鏈接庫(Dynamic Link Library，DLL)等模塊的構(gòu)建。天然氣壓縮因子的計算主要針對實際工況下天然氣的詳細組分進行計算，因此該程序共分為2個部分：人機交互界面部分和核心計算DLL部分。系統(tǒng)界面部分主要負責參數(shù)輸入和校驗、核心計算模塊的調(diào)用和計算結(jié)果的輸出顯示。核心計算DLL部分負責對輸入?yún)?shù)進行計算。系統(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è)置在同一個界面上，組分默認值都設(shè)為0，用戶可根據(jù)計算需要，填入相應(yīng)的氣體組分。由于每次用戶輸入量較多，系統(tǒng)設(shè)置了將用戶的歷次輸入進行自動保存的功能，每次運行時系統(tǒng)會自動調(diào)出最近一次的輸入值，也可找出歷史記錄，同時，用戶可根據(jù)實際輸入情況進行修改。

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

以某天然氣的取樣數(shù)據(jù)為例(見表1)，將表中數(shù)據(jù)輸入軟件界面，點擊“計算”按鈕，便可得出試驗工況和基準工況下的壓縮因子和氣體真實密度，結(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ù)進行計算，計算結(jié)果如表2所示。與附錄C中表C2所給結(jié)果進行對比，結(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計算模型，將求解步驟進行模塊分解，并探討了算法流程，用弦截法對牛頓迭代算法進行優(yōu)化，全面優(yōu)化天然氣壓縮因子算法流程，提高了計算效率。

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

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