基于BWRS狀態(tài)方程的天然氣偏差因子計(jì)算方法

張立俠 郭春秋

中國(guó)石油勘探開發(fā)研究院

天然氣偏差因子（Z）通常又稱壓縮因子或偏差系數(shù)，它表征了某一壓力、溫度條件下相同數(shù)量真實(shí)氣體和理想氣體的體積之比，是油氣藏工程計(jì)算常用的重要參數(shù)；為避免與氣體壓縮系數(shù)（Cg）混淆，后文均稱之為偏差因子。確定偏差因子的方法主要有實(shí)驗(yàn)測(cè)量、圖版法和計(jì)算法3類。實(shí)驗(yàn)法雖然直接可靠，但耗時(shí)、成本高，而圖版法或查表則不能得到連續(xù)的值，因此，計(jì)算法以其簡(jiǎn)便性和實(shí)用性得到了廣泛應(yīng)用。計(jì)算法主要是利用能夠代表偏差因子標(biāo)準(zhǔn)數(shù)據(jù)的經(jīng)驗(yàn)關(guān)系進(jìn)行編程計(jì)算，這些偏差因子關(guān)系式多是源于狀態(tài)方程，或直接通過(guò)回歸分析總結(jié)得到。

利用狀態(tài)方程擬合偏差因子數(shù)據(jù)進(jìn)而總結(jié)經(jīng)驗(yàn)公式的方法是計(jì)算天然氣偏差因子的有效工具。氣體狀態(tài)方程有許多種，如Van der Waals方程［1-3］、RK方程［4-6］、SRK方程［7-9］、PR方程［10-12］等立方型狀態(tài)方程，以及Starling方程［13-14］、BWR方程［15-16］、BWRS方程［17-21］等維里狀態(tài)方程。維里狀態(tài)方程在氣體混合物偏差因子計(jì)算方法的發(fā)展過(guò)程中起到了重要的推動(dòng)作用，較有代表性的有：Hall（1973）［13］在Carnahan-Starling（1969）硬球方程［14］的基礎(chǔ)上提出了一種迭代計(jì)算偏差因子的關(guān)系式，即HY方法；Dranchuk等人（1973）［22］利用BWR狀態(tài)方程擬合偏差因子數(shù)據(jù)［23］提出了一種8系數(shù)表達(dá)式，即DPR方法；Dranchuk和Abou-Kassem（1975）［24］則應(yīng)用BWRS方程提出了一種11系數(shù)偏差因子計(jì)算式，即DAK方法。這3種方法是應(yīng)用狀態(tài)方程計(jì)算天然氣偏差因子的經(jīng)典方法，在一定范圍內(nèi)它們均能比較準(zhǔn)確地代表Standing-Katz圖版數(shù)據(jù)［23］，其中DAK方法的計(jì)算精度最高，但其在高溫高壓下的計(jì)算誤差也稍大，存在改進(jìn)的空間。為此，筆者基于改進(jìn)的BWRS狀態(tài)方程提出一種新的偏差因子計(jì)算方法，以期更精確地確定天然氣偏差因子。

1 改進(jìn)的BWRS狀態(tài)方程

Benedict等人（1940，1942）［15-16］提出了一種剩余功（residual work content）與密度、絕對(duì)溫度之間的經(jīng)驗(yàn)方程，給出了與之對(duì)應(yīng)的狀態(tài)方程，即Benedict-Webb-Rubin方程，其表達(dá)式為方程（1）。

BWR狀態(tài)方程因其能比較準(zhǔn)確地計(jì)算氣體的熱力學(xué)性質(zhì)而得到了廣泛應(yīng)用，例如，Dranchuk等人（1973）［22］利用BWR狀態(tài)方程擬合平滑處理后的Standing-Katz圖版［23］上的1 500個(gè)數(shù)據(jù)點(diǎn)，提出了DPR方法。后來(lái)，許多學(xué)者對(duì)BWR方程進(jìn)行了修正和推廣，以期更加準(zhǔn)確、有效地預(yù)測(cè)更大壓力和溫度范圍內(nèi)的純物質(zhì)及混合物的熱力學(xué)性質(zhì)參數(shù)。其中，Starling（1970）［17-18］修正的BWR方程（即BWRS方程）頗具代表性，應(yīng)用較為廣泛［19-21］，其表達(dá)式為方程（2）。

式中，p為壓力，Pa；R為摩爾氣體常數(shù)，8.3144598 J/（mol·K）；T為溫度，K；ρ為密度，kg/m3；M為摩爾質(zhì)量，kg/mol；B0、A0、C0、b、a、α、c、γ均為與狀態(tài)方程相關(guān)的系數(shù)；D0、E0和d為BWRS方程相較于BWR方程增加的3個(gè)系數(shù)。

BWRS狀態(tài)方程是基于統(tǒng)計(jì)力學(xué)攝動(dòng)理論對(duì)BWR方程中的系數(shù)C0和a相關(guān)項(xiàng)作了修改而得到的。增加的3個(gè)系數(shù)D0、E0和d有助于擬合并預(yù)測(cè)流體在低溫及臨界區(qū)的熱力學(xué)參數(shù)（焓、熵、蒸氣壓等）。利用BWRS方程，Dranchuk（1975）［24］同樣應(yīng)用DPR關(guān)系式所用的1 500個(gè)偏差因子數(shù)據(jù)點(diǎn)進(jìn)行擬合，提出了DAK方法。由于狀態(tài)方程的改進(jìn)，DAK方法比DPR（1973）方法具有更高的計(jì)算精度，但其在15≤ppr≤30 & 1.4≤Tpr≤2.8高壓范圍內(nèi)的計(jì)算結(jié)果與Katz（1959）高壓偏差因子圖版之間的誤差稍大，仍可改進(jìn)。胡建國(guó)（2013）［25］利用6 548組偏差因子數(shù)據(jù)對(duì)DAK方法的11個(gè)系數(shù)進(jìn)行了修正，改善了高壓下的偏差因子計(jì)算精度，但這種修正在0.2≤ppr≤15 & 1.05≤Tpr≤3.0范圍內(nèi)的計(jì)算結(jié)果卻與Standing-Katz圖版數(shù)據(jù)不符。

研究發(fā)現(xiàn)，對(duì)式（2）中的系數(shù)γ進(jìn)行調(diào)整可明顯提高基于BWRS狀態(tài)方程的偏差因子計(jì)算方法的計(jì)算精度。為此，將式（2）改寫為

由真實(shí)氣體狀態(tài)方程得

在式（3）兩端同時(shí)乘以M/(ρRT)，得

引入對(duì)比密度ρpr、對(duì)比溫度Tpr和對(duì)比壓力ppr

則式（5）可化為

其中

式中，γ1、γ2為改進(jìn)的BWRS方程的系數(shù)；Z為氣體偏差因子；Tpr為（擬）對(duì)比溫度；ppr為（擬）對(duì)比壓力；Tpc為（擬）臨界溫度；ppc為（擬）臨界壓力；ρpr為（擬）對(duì)比密度；ρpc為（擬）臨界密度；Zc為臨界偏差因子，計(jì)算時(shí)一般取0.27。

式（7）和式（8）反映了偏差因子Z與Tpr、ppr和ρpr的隱式關(guān)系，相較于DAK方法，該關(guān)系式增加了一個(gè)系數(shù)，同樣可利用迭代法求解Z。

2 偏差因子計(jì)算方法

使用式（7）計(jì)算Z需要先確定A1～A12等12個(gè)系數(shù)的值。Poettmann（1952）［26］、Katz（1959）［27］及Smith（1990）［28］均發(fā)表了0.2≤ppr≤15 & 1.05≤Tpr≤3.0范圍內(nèi)Standing-Katz圖版（圖1）數(shù)值化的標(biāo)準(zhǔn)偏差因子數(shù)據(jù)表（共297×20=5 940個(gè)數(shù)據(jù)點(diǎn)）；Poettmann［26］數(shù)據(jù)表中的5個(gè)數(shù)據(jù)點(diǎn)疑似有誤（因?yàn)樵谶@些點(diǎn)處，Z值發(fā)生了突變），作如表1修改。

圖1 Standing-Katz 偏差因子圖版［23］Fig.1 Standing-Katz Z-factor chart［23］

表1 5個(gè)偏差因子數(shù)據(jù)點(diǎn)的修改值Table 1 Modified values of 5 Z-factor data points

對(duì)于常用的壓力范圍，以Poettmann（1952）數(shù)值化的偏差因子數(shù)據(jù)作為標(biāo)準(zhǔn)；對(duì)于“15≤ppr≤30、1.4≤Tpr≤2.8”范圍內(nèi)的Katz（1959）［27］圖版（圖2），同樣利用數(shù)值化處理結(jié)果作為參照依據(jù)（每條等溫線上對(duì)比壓力ppr每隔0.1取一個(gè)點(diǎn)，共151×8=1 208個(gè)數(shù)據(jù)點(diǎn)）。

圖2 15 ≤ ppr ≤ 30 范圍內(nèi)的Katz 圖版［27］Fig.2 Katz chart for 15 ≤ ppr ≤ 30［27］

為方便數(shù)據(jù)擬合，式（7）可寫為：

函數(shù)g（Tpr，ρpr，Z）以Tpr、ρpr、Z為自變量，且其函數(shù)值應(yīng)恒為0；利用上述7 148組數(shù)據(jù)點(diǎn)分2個(gè)區(qū)擬合該函數(shù)，求滿足式（10）時(shí)各系數(shù)A1～A12的值，此過(guò)程可用Matlab中的非線性最小二乘回歸（nonlinear least squares regression）函數(shù)編程實(shí)現(xiàn)。各系數(shù)A1～A12的取值見(jiàn)表2。

表2 改進(jìn)BWRS狀態(tài)方程的系數(shù)Table 2 Coefficients of modified BWRS Equation

利用牛頓迭代法［29］求解偏差因子Z較為方便，可令函數(shù)f(ρpr)為

其中

f(ρpr)關(guān)于ρpr的導(dǎo)數(shù)為

牛頓迭代公式為

式中，(ρpr)0為上次計(jì)算的對(duì)比密度；(ρpr)1為本次得到的對(duì)比密度。

具體求解步驟如下：

（1）給定ppr、Tpr，設(shè)定最大迭代次數(shù)N和計(jì)算精度eps；

（2）令(ρpr)1=0.27ppr/Tpr，k=0；

（3）將(ρpr)1賦給(ρpr)0；

（4）利用式（13）得到對(duì)比密度(ρpr)1；k+1→k；

（5）比較(ρpr)1和(ρpr)0的絕對(duì)差，若k＞N或abs［(ρpr)1-(ρpr)0］＜eps，則迭代終止；否則，繼續(xù)步驟（3）～（5）。

（6）計(jì)算偏差因子，Z=0.27ppr/［Tpr（ρpr)1］。

3 計(jì)算結(jié)果對(duì)比

利用Standing-Katz圖版數(shù)值化的5 940組數(shù)據(jù)［26-28］和Katz（1959）圖版數(shù)值化的1 208組數(shù)據(jù)對(duì)表2中3種取值方法進(jìn)行評(píng)價(jià)（平均絕對(duì)誤差記為AAE），0.2≤ppr≤15 & 1.05≤Tpr≤3.0范圍（中低壓區(qū)）內(nèi)的計(jì)算誤差記為AAE1，15≤ppr≤30& 1.4≤Tpr≤2.8范圍（高壓區(qū)）內(nèi)的計(jì)算誤差記為AAE2，總共7 148個(gè)數(shù)據(jù)點(diǎn)的計(jì)算誤差記為AAE3。

式中，Zcal為計(jì)算的Z值。

表3列出了這3種相關(guān)的偏差因子計(jì)算方法的平均絕對(duì)誤差；本文提出的12系數(shù)關(guān)系式的AAE1為0.382%，相較于原DAK方法，其計(jì)算精度提高了12.084%；而AAE2為0.205%，比DAK方法的計(jì)算精度高了79.041%；0.2≤ppr≤15 & 1.05≤Tpr≤3.0以及15≤ppr≤30 & 1.4≤Tpr≤2.8范圍內(nèi)各等溫線的計(jì)算誤差見(jiàn)表4、5。

表3 誤差分析Table 3 Error analysis

表4 0.2≤ppr≤15 & 1.05≤Tpr≤3.0范圍內(nèi)各等溫線的平均絕對(duì)誤差 %Table 4 Average absolute error of each isothermal line for 0.2≤ppr≤15 & 1.05≤Tpr≤3.0

表5 15≤ppr≤30 & 1.4≤Tpr≤2.8范圍內(nèi)各等溫線的平均絕對(duì)誤差 %Table 5 Average absolute error of each isothermal line for 15≤ppr≤30 & 1.4≤Tpr≤2.8

圖3、圖4統(tǒng)計(jì)了DAK（1975）［24］、DAK-HJG（2013）［25］和本文提出的新方法在中低壓和高壓范圍內(nèi)的偏差因子計(jì)算誤差，其大小由顏色從藍(lán)到紅來(lái)體現(xiàn)。對(duì)于“0.2≤ppr≤15 & 1.05≤Tpr≤3.0”的中低壓情形，計(jì)算誤差以10%為界，圖3中缺失部分表示誤差大于10%的區(qū)域。同理，對(duì)于“15≤ppr≤30& 1.4≤Tpr≤2.8”的高壓情形，由于計(jì)算誤差整體較小，圖4只展示誤差小于2%的區(qū)域。

圖3 0.2≤ppr≤15 & 1.05≤Tpr≤3.0范圍內(nèi)3種方法的誤差分布Fig.3 Error distribution for 0.2≤ppr≤15 & 1.05≤Tpr≤3.0 by three methods

圖4 15≤ppr≤30 & 1.4≤Tpr≤2.8范圍內(nèi)3種方法的誤差分布Fig.4 Error distribution for 15≤ppr≤30 & 1.4≤Tpr≤2.8 by three methods

對(duì)于中低壓區(qū)，DAK-HJG（2013）方法的缺失區(qū)域較大，藍(lán)色區(qū)域占比較低，其AAE1=12.733%，由此可見(jiàn)，該方法在0.2≤ppr≤15、1.05≤Tpr≤3.0范圍內(nèi)并不適用；而DAK（1975）和新方法的計(jì)算精度較高（AAE1分別為0.435%、0.382%），但新方法誤差大于10%的區(qū)域更小，整體誤差也更小。

對(duì)于高壓區(qū)，DAK（1975）方法在高溫區(qū)域部分點(diǎn)的計(jì)算誤差大于2%（部分區(qū)域缺失），其平均絕對(duì)誤差也最大（AAE2為0.979%）；而DAK-HJG方法和新方法的誤差較?。ˋAE2分別為0.219%和0.205%），但新方法的最大誤差更?。▋H為1.091%）。

筆者認(rèn)為BWRS方程中指數(shù)項(xiàng)中的γ與前一項(xiàng)的γ取值不相同時(shí)該方程能更準(zhǔn)確地代表氣體（尤其是在高壓下）的性質(zhì)；由于增加了一個(gè)系數(shù)，修正后的方程對(duì)于天然氣偏差因子數(shù)據(jù)的代表能力應(yīng)更強(qiáng)。此外，本文使用了0.2≤ppr≤15 & 1.05≤Tpr≤3.0范圍內(nèi)的5 940個(gè)數(shù)據(jù)點(diǎn)（對(duì)Poettmann發(fā)布的偏差因子標(biāo)準(zhǔn)數(shù)據(jù)［26］進(jìn)行了更正）以及15≤ppr≤30 & 1.4≤Tpr≤2.8范圍內(nèi)的1 208個(gè)數(shù)據(jù)點(diǎn)，偏差因子Z數(shù)據(jù)體更大，故新方法適用范圍更廣，計(jì)算精度更高。

4 結(jié)論

（1）基于改進(jìn)的BWRS狀態(tài)方程和DAK方法導(dǎo)出了一種12系數(shù)偏差因子關(guān)系式，利用7 148組偏差因子數(shù)據(jù)結(jié)合非線性回歸分析確定了這些系數(shù)，提出了一種新的Z值確定方法。

（2）由于BWRS方程中指數(shù)項(xiàng)的改進(jìn)，該方程能更準(zhǔn)確地代表及預(yù)測(cè)氣體（尤其是在高壓下）的物理性質(zhì)。

（3）由于偏差因子數(shù)據(jù)的合理性以及狀態(tài)方程的改進(jìn)，新方法的計(jì)算精度更高。新方法在中低壓和高壓下的平均絕對(duì)誤差分別為0.382%、0.205%，這比原DAK方法的計(jì)算精度分別高出12%和79%；推薦使用范圍為：0.2≤ppr≤15 & 1.05＜Tpr≤3.0以及15≤ppr≤30 & 1.4≤Tpr≤2.8。而胡建國(guó)修正的DAK方法只適用于高壓（15≤ppr≤30）的情形。