天然氣偏差因子計(jì)算新方法

張立俠 郭春秋

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

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

計(jì)算法主要是利用能夠代表偏差因子標(biāo)準(zhǔn)圖版的關(guān)系式進(jìn)行編程計(jì)算，其中顯式Z因子關(guān)系式多是通過(guò)回歸分析總結(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)方程在偏差因子計(jì)算方法的發(fā)展過(guò)程中具有重要意義，例如：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ù)偏差因子計(jì)算式即DAK方法[38]。這3種方法是應(yīng)用狀態(tài)方程計(jì)算天然氣偏差因子的經(jīng)典方法[39-44]，在一定范圍內(nèi)它們均能比較準(zhǔn)確地代表Standing-Katz圖版[37]，其中以DAK方法的計(jì)算精度最高[45-48]，但其在高溫高壓下的誤差也稍大。

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

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

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

然而，Nishiumi和Saito 指出BWRS方程仍不能有效地預(yù)測(cè)低對(duì)比溫度下的熱力學(xué)參數(shù)[49]，他們提出了另一種廣義的BWR方程，稱(chēng)之為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ù)。

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

(2)

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

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

(3)

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

(4)

(5)

式中：Tpr為(擬)對(duì)比溫度；Tpc為(擬)臨界溫度,K；ppr為(擬)對(duì)比壓力；ppc為(擬)臨界壓力,Pa；ρpr為(擬)對(duì)比密度;ρpc為(擬)臨界密度,kg/m3；Zc為臨界偏差因子，計(jì)算時(shí)一般取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 偏差因子計(jì)算方法

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

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

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

表1 5個(gè)偏差因子數(shù)據(jù)點(diǎn)的修改值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個(gè)偏差因子數(shù)據(jù)點(diǎn)，對(duì)式(6)作多元非線性回歸分析，得到各系數(shù)A1～A15的取值見(jiàn)表2。

式(6)和式(7)反映了偏差因子和對(duì)比密度ρ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為上次計(jì)算的對(duì)比密度；(ρpr)1為本次得到的對(duì)比密度;B1～B6為系數(shù)。

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

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

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

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

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

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

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

3 方法對(duì)比和討論

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

(13)

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

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

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

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

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

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

表3 誤差分析 Table 3 Error analysis %序號(hào)方法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)各等溫線的平均絕對(duì)誤差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)各等溫線的平均絕對(duì)誤差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)各等溫線的平均絕對(duì)誤差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) 對(duì)比油氣藏工程常用的DPR、HY、DAK方法，該方法的計(jì)算效果更好。除了Tpr=1.05等溫線上1.00≤ppr≤1.45范圍內(nèi)的10個(gè)點(diǎn)的計(jì)算誤差較大外，該方法既適用于常用壓力范圍(0.2≤ppr≤15.0 & 1.05≤Tpr≤3.00)，亦可應(yīng)用于高壓范圍(15.0≤ppr≤30.0 & 1.4≤Tpr≤2.8)，在這兩個(gè)范圍內(nèi)的平均絕對(duì)誤差分別為0.357%、0.066%，優(yōu)于前3種常用方法。