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

姬安召，王玉風，張光生

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

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

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

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

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

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

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

1 基本原理

1.1 曲線差值方法

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

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

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

1.2 曲線間差值方法

2 計算結果

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

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

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

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

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

表1 方案一與方案二計算擬合曲線的左右端點的一階導數值與相對誤差

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

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

3 實例分析

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

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

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

圖3 969 井區井底靜壓誤差統計圖

4 結論

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

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

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