基于AGA8和AGA10的天然氣物性參數計算的軟件設計

陳玉潔 沈昱明

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

關鍵詞： MATLAB;天然氣;物性參數;AGA8-92DC方程;AGA10;優化算法

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

本文著錄格式：陳玉潔，沈昱明. 基于AGA8和AGA10的天然氣物性參數計算的軟件設計[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? 引言

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

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

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

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

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

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

2? 優化計算

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

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

——狀態方程參數，n=1～58

——天然氣組分特征參數

——二元交互作用參數

——比熱、焓和熵的計算參數

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

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

AGA8-92DC方程的計算關系如式（1）所示：

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

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

當滿足時，迭代結束，對應的即為所求解。

2.2? 向量化求解偏導

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

第二維利系數B，僅與天然氣溫度與組分有關，由式（6）計算。

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

第二維里系數偏導計算循環次數不多，優化效果并不明顯。但后續壓縮性系數偏導計算以及積分計算循環量大，引入向量化思想則計算速度提升較為明顯。

2.3? 機械求積法計算積分

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

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

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

3? 軟件設計

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

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

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

4? 測試

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

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

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

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

5? 結語

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

（2）計算結果與AGA報告數據相比，二者完全吻合。

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

參考文獻

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

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

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

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

[5]陳旻灝， 黃希揚， 沈昱明. Matlab與LabVIEW混合編程的天然氣壓縮因子軟件設計[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的車牌自動識別系統設計[J]. 軟件， 2016， 37（3）： 44-46.

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

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

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

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