陳玉潔 沈昱明

摘? 要： 根據(jù)美國燃氣協(xié)會（AGA）發(fā)表的AGA No.8及AGA No.10報告中提出的計算方程，采用向量法、機械求積法等對AGA推薦算法進行了優(yōu)化;并通過MATLAB編程，計算了天然氣物性參數(shù)，包括：壓縮性系數(shù)、比熱、聲速、等熵指數(shù)、焓和熵;分析了AGA8中天然氣壓縮性系數(shù)隨壓力與溫度的變化關(guān)系。測試結(jié)果表明，計算誤差低于0.0003‰，優(yōu)化后程序運行平均耗時僅150 ms左右。

關(guān)鍵詞： MATLAB;天然氣;物性參數(shù);AGA8-92DC方程;AGA10;優(yōu)化算法

中圖分類號： TP391.75? ? 文獻標識碼： A? ? DOI：10.3969/j.issn.1003-6970.2020.09.027

本文著錄格式：陳玉潔，沈昱明. 基于AGA8和AGA10的天然氣物性參數(shù)計算的軟件設(shè)計[J]. 軟件，2020，41（09）：9699+132

【Abstract】： According to the calculation equations proposed in the AGA No. 8 and AGA No. 10 reports issued by the American Gas Association （AGA）， the AGA recommendation algorithm was optimized using the vector method and the mechanical quadrature method; and the physical properties of natural gas was calculated through MATLAB， including compressibility factor， heat capacity， speed of sound， isentropic index， enthalpy and entropy; the relationship between the compressibility factor of natural gas in AGA8 and the pressure and temperature was analyzed. The test results show that the calculation error is less than 0.0003‰， and the average running time of the optimized program is only about 150 ms.

【Key words】： MATLAB; Natural gas; Physical properties; AGA8-92DC equation; AGA10; Optimization

0? 引言

天然氣物性參數(shù)，特別是壓縮性系數(shù)的精確度，對天然氣流量計量影響巨大[1]。對于實際氣體，常采用PV=ZnRT來代替理想氣體狀態(tài)方程。通過引入壓縮性系數(shù)Z，以衡量實際氣體相對理想氣體的偏離程度。

目前，壓縮性系數(shù)的計算方法有RK公式、Redlieh- Kwong對比態(tài)方程[2]、SGERC-88方程、AGA8-92DC方程等。已有大量實驗證實了AGA8-92DC方程的精確性，因此在行業(yè)中被廣泛采用。

文獻3、4、5實現(xiàn)了對壓縮性系數(shù)的編程計算，但計算結(jié)果與AGA8所給結(jié)果存在較大誤差;文獻6探討了AGA8-92DC方程計算值與實驗值的差異;文獻7基于AGA8-92DC的狀態(tài)方程，提出了一種簡便、準確的計算天然氣臨界流質(zhì)量流量的方法。

AGA10雖然給出了計算其他物性參數(shù)的算法及源代碼，但其算法復(fù)雜、代碼量大，不利于芯片化運行。本文在其基礎(chǔ)上進行了拓展，增加求解了天然氣其他重要的物性參數(shù);優(yōu)化了AGA8、AGA10中的幾個偏微分的算法。計算結(jié)果與AGA8、AGA10中的結(jié)果完全吻合，也大大縮短了計算時間。經(jīng)測試，平均計算時間約為150 ms。

1? AGA8-92DC方程及其計算精度

在1992年6月召開的國際標準化組織（ISO）上，美國燃氣協(xié)會（AGA）提出了計算天然氣壓縮性系數(shù)的AGA8-92DC方程（以下簡稱AGA8方程）。AGA8方程是參照天然氣組分摩爾分數(shù)來計算的擴展維里方程。在甲烷含量不低于70%，溫度263 K～338 K，絕對壓力低于20 MPa的工作條件下，計算不確定度優(yōu)于± 0.2%，能滿足絕大部分的工業(yè)條件。

2? 優(yōu)化計算

輸入變量為21種天然氣組分的摩爾分數(shù)、溫度及絕對壓力，MATLAB軟件編程計算后，輸出天然氣物性參數(shù)值，包括壓縮因子、比熱、聲速、等熵指數(shù)、焓和熵。其中，比熱、聲速、等熵指數(shù)、焓、熵的計算公式參照AGA No.10報告[8]。

計算輸入常量說明如下，詳見GB/T17747.2-1999[9]表B.1、表B.2、表B.3和表4：

——狀態(tài)方程參數(shù)，n=1～58

——天然氣組分特征參數(shù)

——二元交互作用參數(shù)

——比熱、焓和熵的計算參數(shù)

圖1給出了天然氣物性參數(shù)的計算流程，下文將說明計算過程采用的優(yōu)化方法。

2.1? 密度檢索法迭代求解壓縮因子

AGA8-92DC方程的計算關(guān)系如式（1）所示：

式（1）～（3）表明壓縮因子、摩爾密度、對比密度、絕對壓力一一對應(yīng)，計算壓縮因子即為求解方程。該方程為非線性方程，無法使用直接求解法求解。為方便說明，令：

本文采用對密度參量進行迭代對壓縮因子求解。圖2為描點曲線圖，可見曲線單調(diào)，滿足二分法使用條件。一般工況條件下，Z處于0.4至1.2之間[3]，故迭代起始條件為：

當滿足時，迭代結(jié)束，對應(yīng)的即為所求解。

2.2? 向量化求解偏導(dǎo)

天然氣物性參數(shù)計算公式較為復(fù)雜，偏導(dǎo)計算更為冗雜。向量化計算是一種特殊的并行化計算方式。由于計算機并行的硬件結(jié)構(gòu)特點，本文發(fā)現(xiàn)，對于偏導(dǎo)計算采用向量化方法，不但可以減少“for循環(huán)”，還可避免大量重復(fù)運算，極大提高了計算效率。以下以第二維里系數(shù)偏導(dǎo)計算予以說明。

第二維利系數(shù)B，僅與天然氣溫度與組分有關(guān)，由式（6）計算。

由此可見，僅需對（n=1～18）構(gòu)成的三個向量進行點積運算即可。

第二維里系數(shù)偏導(dǎo)計算循環(huán)次數(shù)不多，優(yōu)化效果并不明顯。但后續(xù)壓縮性系數(shù)偏導(dǎo)計算以及積分計算循環(huán)量大，引入向量化思想則計算速度提升較為明顯。

2.3? 機械求積法計算積分

比熱、焓和熵的計算涉及求解關(guān)于摩爾密度的定值積分。被積函數(shù)復(fù)雜，包含累加和偏導(dǎo)計算，本文使用機械求積法求解。與定義法相比，機械求積法能夠劃分更少的區(qū)間而達到相同的計算精度，算法時間復(fù)雜度得以降低。物性參數(shù)計算方程為式（10）～（12）所示：

將被積函數(shù)參照AGA No.10報告[8]中式（3.12）和式（3.14）展開，其中常數(shù)項積分直接按照積分表求解，并記被積函數(shù)中關(guān)于摩爾密度的非常數(shù)積分項分別為Q1、Q2、Q3：

積分計算時，本文選用梯形公式代入。表1顯示當選則不同區(qū)間長度，計算相對誤差與平均耗時的變化情況。根據(jù)測試結(jié)果，取區(qū)間長度為0.000002最為合適。

3? 軟件設(shè)計

GUI（Graphical User Interface），指利用圖形控件進行人機交互操作。為了方便用戶的使用，本文引入MATLAB的GUI技術(shù)[10]，只需用戶在人機交互界面上輸入天然氣溫度、壓力和組分摩爾數(shù)，即可自行計算出物性參數(shù)。

圖形界面設(shè)計完成后，通過Matlab的自動導(dǎo)出功能可將計算程序*.m文件與GUI圖形界面*.fig文件合并為一個.m文件。隨后C/C++編譯器生成*.exe可執(zhí)行文件。未安裝MATLAB的機器僅需下載編譯環(huán)境庫MATLAB Runtime（MCR）或者將MATLAB安裝路徑下的/toolbox/compiler/Deploy/win32的MCRInstaller.exe與編譯的執(zhí)行文件打包發(fā)送即可，執(zhí)行文件運行時會自動調(diào)用下載的動態(tài)鏈接庫DLL。

但MCRInstaller.exe所占內(nèi)存較大，可移植性不高，對應(yīng)用范圍有一定限制。既需MATLAB強大的計算能力，又需生成便捷的計算機軟件，也可采用其他編程語言調(diào)用MATLAB程序。如Mideva工具編譯*.m文件后可生成能被Visual C++[12]、Visual Basic以及Delphi調(diào)用的DLL動態(tài)鏈接庫。Microsoft VB在軟件開發(fā)界面上具有獨到之處，但在計算方面略顯不足，這方面恰可與MATLAB互補[13]。圖3是本文軟件界面設(shè)計。

4? 測試

本文將軟件計算結(jié)果與GB/T 17747-1999與AGA No.10報告中算例進行比較，計算結(jié)果中壓縮性系數(shù)完全一致，其他參數(shù)如表1所示，計算相對誤差小于0.0003‰，說明軟件計算正確，滿足精度要求。

MATLAB源程序中，本文利用tic， toc語句對程序進行時間測試。不采用任何優(yōu)化方法計算（默認使用循環(huán)結(jié)構(gòu)和MATLAB自帶積分函數(shù)），平均耗時為1.5 s;優(yōu)化后平均耗時為150 ms，降為原來的10%，優(yōu)化效果明顯。

表2中顯示的是AGA10更寬的應(yīng)用范圍，本文在AGA8原有數(shù)據(jù)的基礎(chǔ)上，計算了更高溫度、更大壓力條件的壓縮性系數(shù)。對于AGA8所展示的6種天然氣，其壓縮性系數(shù)與溫度、壓力的關(guān)系如圖4所示。

可以看出：對同一種類天然氣，壓力為6 MPa的壓縮性系數(shù)變化曲線與12 MPa情況下相差較大，而12 MPa與20 MPa曲線幾乎重合。

5? 結(jié)語

（1）本文采用向量化方法、機械求積法、密度檢索法迭代等優(yōu)化算法，極大提高了計算效率。

（2）計算結(jié)果與AGA報告數(shù)據(jù)相比，二者完全吻合。

（3）對同一種類天然氣，壓力為6 MPa的壓縮性系數(shù)變化曲線與12 MPa情況下相差較大;而12 MPa以上曲線幾乎重合;且隨溫度升高，不同壓力下天然氣壓縮性系數(shù)差值減小。

參考文獻

[1]孟祥適， 姜印平， 劉玉杰， 等. 基于天然氣壓縮系數(shù)Z提高天然氣計量準確度的方法[J]. 測控技術(shù)， 2004， 23（6）： 16-17， 28.

[2]黃希揚， 沈昱明. 基于Redlieh-Kwong對比態(tài)方程的天然氣壓縮因子計算[J]. 化工自動化及儀表， 2019， 46（02）： 107-110.

[3]李丹華， 姜東琪. AGA8-92DC計算方法天然氣壓縮因子計算[J]. 煤氣與熱力， 2011， 31（03）： 43-46.

[4]梁光川， 左果， 文雯， 等. 應(yīng)用VB實現(xiàn)AGA8-92DC氣體壓縮因子計算方法[J]. 遼寧化工， 2014， 43（05）： 609-612.

[5]陳旻灝， 黃希揚， 沈昱明. Matlab與LabVIEW混合編程的天然氣壓縮因子軟件設(shè)計[J]. 化工自動化及儀表， 2019， 46（10）： 844-848.

[6]蘇中良， 劉曰武， 張俊清. AGA8-92DC方法計算天然氣偏差因子的研究[J]. 油氣井測試， 2010， 19（06）： 29-36+83.

[7]D.G. Stewart， J.T.R. Watson， A.M. Vaidya. A new correlation for the critical mass flux of natural gas mixtures[J]. Flow Measurement and Instrumentation， 2000， 11（4）.

[8]AGA Report No. 10 Speed of Sound in Natural Gas and Other Related Hydrocarbon Gases[S]. Transmission Measurement Committee. 2002.11.15.

[9]GB/T 17747.2-2011， 天然氣壓縮因子計算—第2部分： 用摩爾組成進行計算（ISO 12213-2：2006）[S]. 2011， 12， 30.

[10]姚楠， 耿奇. 基于MATLAB GUI的車牌自動識別系統(tǒng)設(shè)計[J]. 軟件， 2016， 37（3）： 44-46.

[11]曹妍， 陳偉. 基于MATLAB GUI的指紋圖像增強技術(shù)研究[J]. 軟件， 2018， 39（11）： 30-32.

[12]徐治. Visual C++調(diào)用MATLAB函數(shù)庫的混合編程技術(shù)[J]. 軟件， 2015， 36（2）： 55-58.

[13]曹妍， 陳偉. 基于MATLAB圖像處理的藥片計數(shù)方法研究[J]. 軟件， 2018， 39（9）： 13-15.

[14]劉艷華. 數(shù)字圖像置亂加密解密技術(shù)的研究與實現(xiàn)[J]. 軟件， 2018， 39（3）： 141-144.