基于樣條曲線插值天然氣偏差因子計算方法

姬安召，王玉風，張光生

（隴東學院，甘肅慶陽 745100）

天然氣偏差因子是在一定溫度壓力下天然氣體積與該氣體在相同條件下按理想氣體定律計算的氣體體積的比值。天然氣的偏差因子是天然氣性質(zhì)研究過程中一個很重要的參數(shù)，例如在密度、黏度、等溫壓縮系數(shù)以及體積系數(shù)計算過程中必須考慮的參數(shù)之一。該參數(shù)主要決定于氣體的組成、溫度以及壓力。天然氣的偏差因子雖然可以通過實驗進行測試，但實驗的周期長，成本高，對實驗儀器的精度要求高。但也可以通過相關(guān)的計算方法得到。目前，天然氣偏差因子獲取大致有三類：圖版插值法、基于狀態(tài)方程法和圖版擬合法的數(shù)值計算方法。

在基于狀態(tài)方程計算天然氣偏差因子的研究方面，對狀態(tài)方程的維里方程法中的AGA8 方法進行了研究，給出了該方法的適應(yīng)性評價[1]。根據(jù)天然氣氣體組成[2]，研究了天然氣壓縮因子的計算，得出對低中壓含碳氣體，AGA8 方程計算誤差最小，而含硫氣體，Piper-Mahmoud 方法誤差最小。根據(jù)實驗研究結(jié)果[3]，對于超高壓氣藏而言，現(xiàn)有的圖版法和經(jīng)驗公式對于天然氣偏差因子計算的誤差很大，甚至有些超壓氣藏天然氣偏差因子誤差超過20%。對圖版擬合的經(jīng)驗公式和基于狀態(tài)方程兩參數(shù)法計算天然氣偏差因子的公式進行評價，給出適合超高壓氣藏天然氣偏差因子的計算模型。根據(jù)氣體狀態(tài)方程的BWRS 維里方程[4]，采用非線性回歸分析，對該方程指數(shù)參數(shù)進行修正，提出了更高精度天然氣偏差因子的計算方法?；赟trobridge 方法[5]，對BWR（Benedict-Webb-Rubin）狀態(tài)方程進行了修正，修正后的方法計算天然氣偏差因子的精度控制在1%以內(nèi)。根據(jù)Nishiumi-Saito 狀態(tài)方程[6]，采用多元回歸方程，給出了天然氣偏差因子計算的解析表達式。但這些方法都要基于不同的狀態(tài)方程，還要進行復(fù)雜的數(shù)學計算，計算的結(jié)果在精度評定時，基本要參照Standing-Katz 圖版的實驗數(shù)據(jù)。

基于圖版擬合方法，在天然氣偏差因子的研究方面，對DAK、HTP、DRP、LXF、HY、Sarem、AGA8-92DC 七種方法計算天然氣偏差因子的結(jié)果進行精度評估[7-10]，給出了每種方法的適用條件。對DAK 方法進行了修正[11]，修正后的方法提高了天然氣偏差因子的計算精度。根據(jù)DPR 方法計算天然氣偏差因子的方法[12]，采用一階中心差分方法計算了超高壓氣藏天然氣等溫擬對比壓縮系數(shù)，拓展了Trube 使用范圍。利用低壓傳統(tǒng)天然氣偏差因子圖版數(shù)據(jù)和實驗測試的高壓天然氣數(shù)據(jù)進行擬合[13]，得到了適合于中低壓和高壓的天然氣偏差因子計算GXH 模型。這些圖版擬合結(jié)果，基本都是基于Standing-Katz 圖版數(shù)字化后的結(jié)果，或者有些參照了狀態(tài)方程的基本形式，其不同的模型都用一定的適用條件。

除了基于狀態(tài)方程法和圖版擬合法之外，還有采用BP 神經(jīng)網(wǎng)絡(luò)進行預(yù)測方法，以測試壓力、溫度、天然氣的組成以及其臨界參數(shù)、測試的氣體的偏差因子作為輸入[14]，采用BP 神經(jīng)網(wǎng)絡(luò)方法對249 個實驗數(shù)據(jù)點進行了樣本學習，建立基于BP 神經(jīng)網(wǎng)絡(luò)的天然氣偏差因子的預(yù)測模型，但該模型對輸入數(shù)據(jù)的依賴性強，實用性不強?；跉獠匚镔|(zhì)平衡原理進行天然氣偏差因子計算的模型[15]，根據(jù)生產(chǎn)動態(tài)資料，結(jié)合干氣氣藏的物質(zhì)平衡方程，提出了基于生產(chǎn)資料計算天然氣偏差因子的新方法，但計算的精度依賴于氣藏特性和生產(chǎn)動態(tài)資料的準確性。在商業(yè)軟件進行天然氣偏差因子的評價方面，對目前流行的四種商業(yè)軟件[16]（PIPESIM 2012、PROPSER 13.0、PVTp 10.0、Multiflash 4.1）對天然氣偏差因子計算模塊進行了評價，其中PROPSER 13.0 計算精度最好，對比不同的狀態(tài)方程計算的結(jié)果，SRK 方程計算結(jié)果最理想。

天然氣中一般還有酸性氣體（CO2、H2S），對于還有酸性氣體的天然氣計算的偏差因子需要做矯正。應(yīng)用酸性氣體偏差因子的計算方法[17]，根據(jù)Keys、SBV 以及RK 混合規(guī)則，酸性氣體的臨界參數(shù)校正方法，優(yōu)選出了CO2氣藏的偏差因子計算的高精度RK+WA 模型。研究富含酸性氣CO2天然氣偏差因子的計算方法[18-20]，通過實驗結(jié)果與理論模型的對比，優(yōu)選出來基于GXO校正的DAK 模型計算誤差最小，并給出了富含CO2天然氣偏差因子隨溫度與壓力變化規(guī)律以及CO2含量對偏差因子計算結(jié)果的影響。通過對比目前天然氣中含有酸性氣體（CO2、H2S）時偏差因子的矯正方法[21]，得出DPR 模型和DAK 模型對含有酸性氣體的矯正后天然氣偏差因子計算結(jié)果最好。其他模型（GXQ、CKB）結(jié)果不理想。

根據(jù)前人的研究結(jié)果，以Poettman、Katz 和Smith對Standing-Katz 圖版低壓范圍（0.2 ≤Ppr≤15.0、1.05≤Tpr≤3.00）數(shù)值化的5 940 個數(shù)據(jù)點[22，23]以及Poettman 對高壓范圍（15.0≤Ppr≤30.0、1.4≤Tpr≤2.8）數(shù)值化Katz 圖版的1 208 個數(shù)據(jù)點為基礎(chǔ)，采用樣條函數(shù)差值的原理，對同一視對比溫度Tpr下的天然氣視對比壓力與偏差因子進行樣條函數(shù)差值，既保證了實驗數(shù)據(jù)點不偏離差值函數(shù)，又保證了差值曲線的連續(xù)光滑性。若視對比溫度Tpr不在實驗數(shù)據(jù)給定的曲線上，則根據(jù)實驗數(shù)據(jù)給定的視對比溫度Tpr值在Standing-Katz 圖版上選擇距離被差值的視對比溫度Tpr最近的兩條曲線，進行樣條曲線間的二次差值，則可以得到滿足Standing-Katz 圖版的任意視對比溫度Tpr和視對比壓力Ppr下的天然氣偏差因子的值。如果天然氣中含有酸性氣體，則進行矯正。

1 基本原理

1.1 曲線差值方法

根據(jù)Standing-Katz 圖版數(shù)值化的數(shù)據(jù)點，選取任意一條以視對比溫度Tpr為參數(shù)數(shù)值化離散數(shù)據(jù)點，將視對比壓力Ppr作為自變量，Standing-Katz 圖版上偏差因子Z 為因變量。在每個相鄰的視對比壓力之間的子區(qū)間進行三次多項式差值，并且要求在整個區(qū)域子區(qū)間上的差值多項式函數(shù)具有二階連續(xù)導(dǎo)數(shù)，這樣得到差值函數(shù)即為三次樣條差值函數(shù)。這樣既保證了區(qū)間內(nèi)端點的連續(xù)性，也保證了內(nèi)端點處的光滑性，并且被差值的數(shù)據(jù)點落在差值曲線上。

假設(shè)在［Pprmin，Pprmax］區(qū)間上有n+1 個離散數(shù)據(jù)點Ppr（0）＜Ppr（1）＜…＜Ppr（n-1）＜Ppr（n），Pprmin為差值區(qū)間視對比壓力的最小值，Pprmax為差值區(qū)間視對比壓力的最大值。天然氣偏差因子Zi（Ppr）在區(qū)間［Ppr（i-1），Ppr（i）］上是三次多項式函數(shù)，這個三次多項式函數(shù)可表示為，即有4 個待定系數(shù)。因為n+1 個離散數(shù)據(jù)點有n 個子區(qū)間，每個子區(qū)間要確定一個三次多項式，因此共計有4n 個待定系數(shù)需要計算。根據(jù)差值的基本條件，即共計n+1 個條件。在區(qū)間內(nèi)端點處，Zi（Ppr）、Zi（Ppr）的一階導(dǎo)函數(shù)以及二階導(dǎo)函數(shù)為連續(xù)函數(shù)，即可得3n-3個條件。還需要兩個條件，即邊界條件，這里取差值區(qū)間的兩個端點處的一階導(dǎo)數(shù)即可，其值由區(qū)間起點與終點的相鄰兩點的直線的斜率確定。根據(jù)上述分析，在區(qū)間x∈［Ppr（i-1），Ppr（i）］上的三次樣條差值多項式的數(shù)學模型可由（1）式表示：

其中：i=1，2，3，…，n。通過（2）式，可以計算位于差值區(qū)間任意視對比壓力下的天然氣偏差因子。

1.2 曲線間差值方法

2 計算結(jié)果

根據(jù)（2）式的分析，在進行曲線差值之前，需要計算Standing-Katz 圖版數(shù)值化后的每條曲線兩個端點的一階導(dǎo)函數(shù)。這里分別取曲線每條數(shù)值化兩個相鄰端點的點進行一階導(dǎo)數(shù)的計算。為了分析一階導(dǎo)數(shù)對差值結(jié)果的影響以及偏差因子的計算精度，這里分別選取兩個方案，方案一：以曲線端點相鄰兩點之間的直線的斜率代替一階導(dǎo)數(shù)。方案二：以曲線端點相鄰三點進行一階導(dǎo)數(shù)加權(quán)計算，其中方案二的一階導(dǎo)數(shù)計算公式由（4）式表示。

其中：當k=0 時，則j=1；當k=n 時，則j=n-1。

根據(jù)Standing-Katz 圖版數(shù)值化的結(jié)果，結(jié)合上述一階導(dǎo)數(shù)的計算方法，曲線的兩個端點的一階導(dǎo)數(shù)的結(jié)果（見表1）。從表1 可以看出，方案一與方案二計算的結(jié)果差值曲線兩個端點的一階導(dǎo)數(shù)在左端點和右端點的相對誤差很小，基本控制在5%以內(nèi)，其中視對比溫度為1.8 和1.9 時，方案一與方案二計算的左端點的一階導(dǎo)數(shù)值誤差略大于5%。

根據(jù)表1 中方案一和方案二一階導(dǎo)數(shù)的計算結(jié)果，這里取二者的平均值，將二者平均值代入公式（1）求解差值內(nèi)節(jié)點的一階導(dǎo)數(shù)值，將內(nèi)節(jié)點的一階導(dǎo)數(shù)值代入到（2）式，得到對應(yīng)視對比溫度下的樣條函數(shù)；然后在區(qū)間［0.2，30］上離散視對比壓力，這里視對比壓力離散的間距為0.001，然后將離散的視對比壓力代入到（2）式，得到樣條函數(shù)差值后結(jié)果（見圖1），圖1中的點表示Standing-Katz 圖版的數(shù)值化的結(jié)果。這里需要說明的是Standing-Katz 圖版數(shù)值化的數(shù)據(jù)太多，在繪圖時按照一定間隔進行了抽稀處理。通過圖1 的對比，Standing-Katz 圖版的數(shù)值化的結(jié)果與樣條差值的計算結(jié)果匹配很好，并且樣條差值法可以計算滿足Standing-Katz 圖版差值范圍的任意給定視對比溫度與視對比壓力下的天然氣偏差因子的值。

圖1 Standing-Katz 圖版與樣條曲線擬合結(jié)果對比圖

表1 方案一與方案二計算擬合曲線的左右端點的一階導(dǎo)數(shù)值與相對誤差

根據(jù)曲線間的差值原理，這里選取了視對比溫度分別為1.08、1.55 和2.10 的三組參數(shù)。將這三組參數(shù)與Standing-Katz 圖版的視對比溫度對比可知，這三組參數(shù)均不是Standing-Katz 圖版所給定的參數(shù)，因此這里只能通過曲線間差值的方法得到不同視對比壓力下的偏差因子值。這三組視對比溫度參數(shù)分別介于視對比溫度為［1.05，1.10］、［1.50，1.60］和［2.00，2.20］三個區(qū)間，根據(jù)（3）式的計算方法，得到這三組參數(shù)差值的結(jié)果（見圖2），Standing-Katz 圖版數(shù)值化的數(shù)據(jù)見圖2的離散點所示。從圖2 可以看出，差值的實線的結(jié)果與Standing-Katz 圖版數(shù)值化的變化趨勢一致。

圖2 Standing-Katz 圖版與樣條曲線擬合結(jié)果對比圖

3 實例分析

在垂直管流計算過程中，天然氣偏差系數(shù)的計算模型有很多種，根據(jù)前人的研究結(jié)果，本次研究主要采用了四種計算方法，即DAK（Dranchk-Abu-Kasse）模型，HTP（Hankinson-Thomas-Phillips）模型，李相方（LXF）模型和本文提出的樣條函數(shù)差值法。

根據(jù)文獻提出的氣井井底靜壓的計算原理[25]，采用微元分析法的思想，將氣井的井筒從井口到井底離散為多個單元，每一個計算微元體內(nèi)的性質(zhì)是一致的，其屬性值可以看作常數(shù)處理?；谝陨峡紤]，可以把井筒從井口到產(chǎn)層中深劃分為n 個微元體，在每個微元體內(nèi)，天然氣的溫度、偏差因子、黏度、壓力以及與壓力有關(guān)的量可以看作常數(shù)（都是基于前一個節(jié)點的壓力與溫度計算出的值）。

由于969 井區(qū)的天然氣中含有酸性氣體，考慮到酸性氣體對天然氣偏差因子的影響，采用了GXQ 方法和Wichert-Aziz 方法對天然氣的臨界溫度和壓力分別做了校正，具體方法參見文獻[19-21]。根據(jù)試氣資料的整理，統(tǒng)計了969 井區(qū)共計76 口井生產(chǎn)前的壓力恢復(fù)數(shù)據(jù)，這些井也下入了井下壓力計，測量了實際的靜壓。結(jié)合垂直管流微元分析法的計算思想，根據(jù)關(guān)井恢復(fù)后的井口油壓和套壓的測量值計算了這76 口井的靜壓。將靜壓計算的結(jié)果與實測的進行對比，得到了969井區(qū)分別采用Dranchk-Abu-Kasse 模型[6]、Hankinson-Thomas-Phillips 模型[8]、李相方（LXF）模型[7]以及樣條函數(shù)差值模型計算的天然氣偏差因子，然后根據(jù)微元分析進行井底靜壓計算，井底靜壓的誤差分布（見圖3）。從圖3 可以看出，本文給出的樣條函數(shù)差值法計算的偏差因子在計算井底靜壓時的相對誤差最?。黄浯螢長XF 法，因為該方法在計算偏差因子時，采用了分段函數(shù)進行處理，適用性較好；然后為DAK 和HTP法。其中井底靜壓相對誤差大于5%的這些井，基本為氣水同產(chǎn)，甚至有些井存在著井底積液，導(dǎo)致井底靜壓計算的相對誤差較大。

圖3 969 井區(qū)井底靜壓誤差統(tǒng)計圖

4 結(jié)論

（1）根據(jù)Standing-Katz 圖版數(shù)值化的結(jié)果，采用三次樣條函數(shù)的差值原理，給出了三次樣條差值計算天然氣偏差因子的基本原理。

（2）對于視對比溫度Tpr不在Standing-Katz 圖版上的情況，在Standing-Katz 圖版上查找距離給定視對比溫度Tpr最近的兩條曲線上的視對比溫度，然后根據(jù)視對比溫度Tpr的值進行反距離加權(quán)二次差值。根據(jù)反距離加權(quán)二次差值，可以得到任意給定視對比溫度Tpr（1.05≤Tpr≤3.00）和視對比壓力Ppr（0.2≤Ppr≤15.0）下的天然氣偏差因子的值。

（3）通過969 井區(qū)76 口氣井的試氣資料，本文提出的三次樣條函數(shù)差值計算偏差因子方法具有較高的精度，并且計算簡單，適用性強。