張立俠 郭春秋

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

天然氣偏差因子(Z)是油氣藏工程、天然氣化工等領(lǐng)域常用的重要參數(shù)，它表征了某一溫度、壓力條件下相同數(shù)量真實氣體與理想氣體體積的比值，通常又稱壓縮因子或偏差系數(shù)。確定偏差因子的方法主要有3類，即：實驗測量法、圖版法(查表法)和計算法。實驗測量法雖然直接可靠，但耗時、成本高；圖版法(查表法)難以得到連續(xù)的值，因此計算法以其簡便性和實用性被廣泛應(yīng)用。

計算法主要是利用能夠代表偏差因子標(biāo)準(zhǔn)圖版的關(guān)系式進(jìn)行編程計算，其中顯式Z因子關(guān)系式多是通過回歸分析總結(jié)得到，如Beggs、Gopal、Kumar、Azizi、Heidaryan和Salarabadi、Heidaryan和Moghadasi等提出的表達(dá)式[1-8]；而隱式關(guān)系式則一般源于狀態(tài)方程。狀態(tài)方程有Van der Waals方程、RK方程、SRK方程、PR方程等立方型狀態(tài)方程[9-19]，以及Kamerlingh Onnes方程、Beattie方程、BWR方程、Strobridge方程、Carnahan方程、Morsy方程、BWRS方程、LeE-Kesler方程等維里狀態(tài)方程[20-35]。維里狀態(tài)方程在偏差因子計算方法的發(fā)展過程中具有重要意義，例如：Dranchuk等利用BWR方程擬合1500組偏差因子數(shù)據(jù)提出了8系數(shù)偏差因子表達(dá)式即DPR方法[36-37];Hall和Yarborough在Carnahan-Starling硬球方程的基礎(chǔ)，提出了一種Z因子關(guān)系式即HY方法[27-28]；Dranchuk和Abou-Kassem 則應(yīng)用BWRS方程提出了一種11系數(shù)偏差因子計算式即DAK方法[38]。這3種方法是應(yīng)用狀態(tài)方程計算天然氣偏差因子的經(jīng)典方法[39-44]，在一定范圍內(nèi)它們均能比較準(zhǔn)確地代表Standing-Katz圖版[37]，其中以DAK方法的計算精度最高[45-48]，但其在高溫高壓下的誤差也稍大。

為更加準(zhǔn)確地預(yù)測整個壓力范圍(0.2≤ppr≤30.0)內(nèi)的天然氣偏差因子，本文基于Nishiumi-Saito方程[49]，提出一種新的偏差因子計算方法，以期為油氣藏工程相關(guān)研究者提供參考。

1 氣體狀態(tài)方程分析

利用狀態(tài)方程擬合偏差因子數(shù)據(jù)進(jìn)而總結(jié)經(jīng)驗關(guān)系式的方法是計算天然氣偏差因子的有效工具。Benedict等 提出的BWR方程頗具代表性[22-24]，因其能比較準(zhǔn)確地計算氣體的熱力學(xué)性質(zhì)而迅速得到應(yīng)用。Dranchuk 基于BWR方程提出了DPR方法[36],隨后許多學(xué)者對其進(jìn)行了修正和推廣，以便更加準(zhǔn)確、有效地預(yù)測更大溫度與壓力范圍內(nèi)的純物質(zhì)和混合物的熱力學(xué)性質(zhì)參數(shù)。其中，Starling修正的BWR方程(即BWRS方程)應(yīng)用更為廣泛[30-34]，DAK方法便是基于此。

然而，Nishiumi和Saito 指出BWRS方程仍不能有效地預(yù)測低對比溫度下的熱力學(xué)參數(shù)[49]，他們提出了另一種廣義的BWR方程，稱之為Nishiumi-Saito方程，其表達(dá)式如式(1)：

(1)

式中：p為壓力，Pa；R為摩爾氣體常數(shù)[50-51]，8.314 459 8 J/(mol·K)；T為溫度，K；ρ為密度，kg/m3；M為摩爾質(zhì)量，kg/mol；B0、A0、C0、D0、E0、b、a、d、e、f、α、c、g、h、γ為與狀態(tài)方程相關(guān)的系數(shù)。

真實氣體狀態(tài)方程為：

(2)

式中：Z為氣體偏差因子。

在式(1)兩端同時乘以M/(ρRT)，得到：

(3)

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

(4)

(5)

式中：Tpr為(擬)對比溫度；Tpc為(擬)臨界溫度,K；ppr為(擬)對比壓力；ppc為(擬)臨界壓力,Pa；ρpr為(擬)對比密度;ρpc為(擬)臨界密度,kg/m3；Zc為臨界偏差因子，計算時一般取0.27。

將式(4)和式(5)代入式(3)，整理得到：

Z= 1+(A1-A2Tpr-1-A3Tpr-3+A4Tpr-4-

A5Tpr-5)ρpr+(A6-A7Tpr-1-A8Tpr-2-

A9Tpr-5-A10Tpr-24)ρpr2+A11(A7Tpr-1+

A8Tpr-2+A9Tpr-5+A10Tpr-24)ρpr5+

(A12Tpr-3+A13Tpr-9+A14Tpr-18)ρpr2(1+

A15ρpr2)exp(-A15ρpr2)

(6)

(7)

(8)

式(6)即為基于Nishiumi-Saito狀態(tài)方程導(dǎo)出的新的隱式關(guān)系式。式(8)中A1～A15為系數(shù)。

2 偏差因子計算方法

Poettman 率先發(fā)表了從Standing-Katz圖版(見圖1)上讀取的偏差因子數(shù)據(jù)表(0.2≤ppr≤15.0 & 1.05≤Tpr≤3.00)[37,52]，Katz 和Smith 沿用了Poettman的數(shù)值化成果并對個別數(shù)據(jù)點進(jìn)行了更正[53-54]。

然而，發(fā)現(xiàn)Poettman數(shù)據(jù)表(共297×20=5940個數(shù)據(jù)點)中的5個數(shù)據(jù)點疑似有誤，因為在這些點處，Z值發(fā)生了突變。參考這些點附近壓力、溫度范圍內(nèi)Z的取值，相應(yīng)作了更正(見表1)。表1中，Z1、Z2、Z3所在列分別為Poettman、Katz、Smith數(shù)據(jù)的偏差因子取值，Zsta所在列為本文更正的偏差因子值。

對于常用的壓力范圍(即中低壓力范圍0.2≤ppr≤15.0 & 1.05≤Tpr≤3.00)，以更正后的Poettman數(shù)值化的偏差因子數(shù)據(jù)作為標(biāo)準(zhǔn)。類似地，對于高壓范圍15.0≤ppr≤30.0 & 1.40≤Tpr≤2.80)內(nèi)的Katz圖版(見圖2)[53]，同樣利用數(shù)值化處理結(jié)果作為標(biāo)準(zhǔn)數(shù)據(jù)(此處每條等溫線上ppr每隔0.1取一個點，共151×8=1208個數(shù)據(jù)點)。

表1 5個偏差因子數(shù)據(jù)點的修改值Table 1 Five corrected values for Z-factor dataTprpprZ1Z2Z3Zsta1.202.500.4190.5190.5190.5191.800.250.9980.9980.9880.9881.057.650.9870.9870.9870.9771.158.500.0461.0461.0461.0461.1014.351.6501.6501.6501.654

利用上述常用壓力范圍和高壓范圍的7148個偏差因子數(shù)據(jù)點，對式(6)作多元非線性回歸分析，得到各系數(shù)A1～A15的取值見表2。

式(6)和式(7)反映了偏差因子和對比密度ρpr的隱式關(guān)系。采用牛頓迭代法可快速解得Z值[55]。首先構(gòu)造如下函數(shù)f(ρpr)：

表2 Nishiumi-Saito狀態(tài)方程的系數(shù)Table 2 Coefficients of the Nishiumi-Saito Equation of State系數(shù)0.2≤ppr≤15.015.0

f(ρpr)= 1+B1·ρpr+B2·ρpr2+B3·ρpr5+

B4(ρpr2+B5ρpr4)exp(-B5ρpr2)-

B6·ρpr-1

(9)

B1=A1-A2Tpr-1-A3Tpr-3-A4Tpr-4-A5Tpr-5B2=A6-A7Tpr-1-A8Tpr-2-A9Tpr-5-A10Tpr-24B3=A11(A7Tpr-1+A8Tpr-2+A9Tpr-5+A10Tpr-24)B4=A12Tpr-3+A13Tpr-9+A14Tpr-18B5=A15B6=0.27pprTpr-1

(10)

求f(ρpr)關(guān)于ρpr的一階導(dǎo)數(shù)：

f′(ρpr)=B6·ρpr-2+B1+2B2·ρpr+5B3·ρpr4+

(11)

建立牛頓迭代公式[55]：

(12)

式中：(ρpr)0為上次計算的對比密度；(ρpr)1為本次得到的對比密度;B1～B6為系數(shù)。

結(jié)合式(9)～式(12)，設(shè)置如下迭代步驟求解偏差因子：

(1) 給定ppr、Tpr，設(shè)定計算精度eps和最大迭代次數(shù)M。

(2) 假定偏差因子Z0=1,令(ρpr)1=B6/Z0，迭代次數(shù)m=0。

(3) 將(ρpr)1賦給(ρpr)0。

(4) 利用式(12)求得對比密度(ρpr)1；將m+1賦值給m。

(5) 比較(ρpr)1和(ρpr)0的差距，若到達(dá)精度要求或迭代次數(shù)超過限制(即abs[(ρpr)1-(ρpr)0]M)則終止迭代；否則，繼續(xù)步驟(3)～(5)。

(6) 計算偏差因子，Z=B6/(ρpr)1。

3 方法對比和討論

利用Standing-Katz圖版數(shù)值化的5940組數(shù)據(jù)和Katz圖版數(shù)值化的1208組數(shù)據(jù)[37，52-54]，對DPR、HY、DAK和本方法進(jìn)行誤差評價，誤差計算式如式(13)：

(13)

式中：AAE為平均絕對誤差，%；Zcal為計算的Z值；Zsta為Z的圖版數(shù)值化值。

中低壓范圍(0.2≤ppr≤15.0 & 1.05≤Tpr≤3.00)內(nèi)的計算誤差記為AAE1，相對高壓范圍(15.0≤ppr≤30.0 & 1.4≤Tpr≤2.8)內(nèi)的誤差記為AAE2，總共7148個數(shù)據(jù)點的平均絕對誤差記為AAE3。

表3列出了4種偏差因子計算方法的平均絕對誤差；表4～表6列出了“0.2≤ppr≤15.0 & 1.05≤Tpr≤3.00”以及“15.0≤ppr≤30.0 & 1.4≤Tpr≤2.8”范圍內(nèi)各等溫線的計算誤差。由表3～表6可看出，本方法的AAE1為0.357%，AAE2為0.066%，其計算精度高于其他3種方法，對比低溫(尤其是Tpr=1.1)和高壓(ppr>15.0)條件下的計算結(jié)果，優(yōu)勢更為明顯。

圖3～圖10統(tǒng)計了DPR、HY、DAK和本方法在中低壓和高壓范圍內(nèi)的偏差因子計算誤差，圖中顏色從藍(lán)到紅體現(xiàn)誤差由小變大。在中低壓(0.2≤ppr≤15.0 & 1.05≤Tpr≤3.00)范圍內(nèi)，4種方法在Tpr=1.05等溫線上均存在誤差較大的點，圖3～圖6中缺失部分表示誤差大于10%的區(qū)域。在高壓(15.0≤ppr≤30.0 & 1.4≤Tpr≤2.8)范圍內(nèi)，圖7和圖9中缺失部分表示誤差大于2%的區(qū)域。

對于中低壓區(qū)，4種方法中HY方法的缺失區(qū)域最大，在“1.30≤ppr≤2.75 &Tpr=1.05”和“1.50≤ppr≤1.65 &Tpr=1.1”范圍內(nèi)的共計34個點的誤差大于10%(DPR有31個點，DAK有29個點)；DPR方法的平均絕對誤差最大；而本方法的計算精度最高，且誤差大于10%的區(qū)域最小(只有“1.00≤ppr≤1.45 &Tpr=1.05”范圍內(nèi)的10個點)，整體誤差也最小。

對于高壓區(qū)，DPR和DAK方法在部分高溫區(qū)域的計算誤差大于2%，DAK方法的平均絕對誤差最大(AAE2為0.979%)，HY方法的誤差相對較小但誤差分布不平衡，而本方法的誤差分布十分平滑且AAE2僅為0.066%。

表3 誤差分析 Table 3 Error analysis %序號方法AAE1AAE2AAE31DPR0.5140.9230.5832HY0.4410.4180.4373DAK0.4350.9790.5274本方法0.3570.0660.307

表4 0.2≤ppr≤15.0 & 1.05≤Tpr≤1.50范圍內(nèi)各等溫線的平均絕對誤差Table 4 Isotherm average absolute errors in the range of 0.2≤ppr≤15.0 & 1.05≤Tpr≤1.50%方法Tpr=1.05Tpr=1.10Tpr=1.15Tpr=1.20Tpr=1.25Tpr=1.30Tpr=1.35Tpr=1.40Tpr=1.45Tpr=1.50DPR2.8671.0740.3510.3540.3800.3980.3080.3320.2620.228HY2.9321.2600.3690.2870.2930.2860.2750.2980.2150.217DAK2.6721.1280.3860.2630.2920.3420.2110.2200.1220.142本方法2.4690.4870.3330.3030.2310.2350.1640.1710.1760.187

表5 0.2≤ppr≤15.0 & 1.60≤Tpr≤3.00范圍內(nèi)各等溫線的平均絕對誤差Table 5 Isotherm average absolute errors in the range of 0.2≤ppr≤15.0 & 1.60≤Tpr≤3.00%方法Tpr=1.60Tpr=1.70Tpr=1.80Tpr=1.90Tpr=2.00Tpr=2.20Tpr=2.40Tpr=2.60Tpr=2.80Tpr=3.00DPR 0.3740.4000.2780.2830.2390.3690.4040.4170.4580.510HY0.1790.1650.2070.2160.2250.2450.2190.1980.2620.471DAK0.2950.3060.2930.2530.1880.2050.2510.2920.3520.488本方法0.1920.1800.1850.1760.1730.2280.2280.2140.2990.501

表6 15.0≤ppr≤30.0 & 1.4≤Tpr≤2.8范圍內(nèi)各等溫線的平均絕對誤差Table 6 Isotherm average absolute errors in the range of 15.0≤ppr≤30.0 & 1.4≤Tpr≤2.8%方法Tpr=1.4Tpr=1.6Tpr=1.8Tpr=2.0Tpr=2.2Tpr=2.4Tpr=2.6Tpr=2.8DPR0.4940.6870.5180.5310.8031.0641.4101.880HY0.5961.1640.5340.2180.0440.0720.2300.489DAK0.2851.0110.7870.7040.8601.0451.3461.795本方法0.0580.0510.0610.1140.0580.0710.0370.076

4 結(jié) 論

(1) 基于Nishiumi修正的BWR狀態(tài)方程(Nishiumi-Saito方程)導(dǎo)出了15系數(shù)偏差因子關(guān)系式，進(jìn)而提出了一種確定氣體偏差因子的新方法。

(2) 對比油氣藏工程常用的DPR、HY、DAK方法，該方法的計算效果更好。除了Tpr=1.05等溫線上1.00≤ppr≤1.45范圍內(nèi)的10個點的計算誤差較大外，該方法既適用于常用壓力范圍(0.2≤ppr≤15.0 & 1.05≤Tpr≤3.00)，亦可應(yīng)用于高壓范圍(15.0≤ppr≤30.0 & 1.4≤Tpr≤2.8)，在這兩個范圍內(nèi)的平均絕對誤差分別為0.357%、0.066%，優(yōu)于前3種常用方法。