鄭 琴 董傳杰 高淑芳

（中國(guó)石油大港油田勘探開發(fā)研究院）

求解異常高壓氣藏單相非線性滲流模型的一種新方法

鄭 琴 董傳杰 高淑芳

（中國(guó)石油大港油田勘探開發(fā)研究院）

針對(duì)傳統(tǒng)異常高壓定參數(shù)滲流模型不能滿足工程需要的現(xiàn)狀，建立滲透率、孔隙度等參數(shù)隨有效應(yīng)力變化的異常高壓氣藏單相非線性滲流模型。為使井底附近壓力分布的計(jì)算更加準(zhǔn)確，采用對(duì)數(shù)等距網(wǎng)格進(jìn)行網(wǎng)格剖分，即越靠近井底處，網(wǎng)格分布越密，沿徑向向外，網(wǎng)格逐漸變稀。在求解差分方程時(shí)，創(chuàng)新采用自適應(yīng)延拓法，克服了傳統(tǒng)Newton法、Brown法等數(shù)值方法要求迭代初值與真解充分接近的局限性。編程計(jì)算結(jié)果表明：考慮了滲透率、孔隙度等參數(shù)隨有效應(yīng)力變化規(guī)律的異常高壓氣藏單相非線性滲流模型比常用的常系數(shù)模型更加貼近工程實(shí)際，自適應(yīng)延拓法的應(yīng)用對(duì)解決初值問(wèn)題及奇異問(wèn)題也是簡(jiǎn)單有效的。

異常高壓 非線性 滲流 有限元 數(shù)值模擬

國(guó)內(nèi)現(xiàn)已發(fā)現(xiàn)多個(gè)異常高壓氣田。氣藏衰竭式開采過(guò)程中，流體壓力逐漸下降，儲(chǔ)層巖石骨架受到的有效應(yīng)力會(huì)隨之增大，導(dǎo)致巖石發(fā)生顯著的彈塑性形變，并減小巖石孔隙度、滲透率、壓縮系數(shù)等物性參數(shù)的值，進(jìn)而影響到氣藏的最終開發(fā)效果。因此，在進(jìn)行異常高壓氣藏研究時(shí)應(yīng)充分考慮儲(chǔ)層應(yīng)力敏感性的影響，掌握儲(chǔ)層物性參數(shù)隨有效應(yīng)力改變的變化規(guī)律，建立能更準(zhǔn)確描述異常高壓氣藏滲流機(jī)理的非線性滲流模型，進(jìn)而更加精確地對(duì)該類氣藏進(jìn)行開發(fā)動(dòng)態(tài)模擬和最優(yōu)控制。

1 儲(chǔ)層物性參數(shù)的應(yīng)力敏感性

國(guó)內(nèi)外專家從上世紀(jì)70年代起，就開始研究孔隙度隨凈圍壓變化的規(guī)律，他們得出了孔隙度隨有效壓力呈指數(shù)式變化的結(jié)論［1－4］。這種指數(shù)式關(guān)系可用數(shù)學(xué)表達(dá)式寫為：

類似，滲透率與有效壓力之間的關(guān)系也可用指數(shù)函數(shù)式來(lái)表達(dá)：

2 異常高壓?jiǎn)蜗喾蔷€性滲流微分方程的建立

氣體的等溫壓縮系數(shù)：

由（1）、（2）知地層滲透率、孔隙度隨有效應(yīng)力呈指數(shù)式變化。由狀態(tài)方程、運(yùn)動(dòng)方程和連續(xù)性方程最終可得：

在式（4）中，由于K、Z都是壓力的函數(shù)，因此不能提到算子之外，引入擬壓力函數(shù)的概念［5］：

經(jīng)變換可得到以擬壓力形式表示的異常高壓氣藏不穩(wěn)定滲流微分方程：

因αΦ的值與氣體壓縮系數(shù)C（p）比較起來(lái)非常的小，可以忽略不計(jì)。令c＝C（p）＋αΦ≈C（p），則得到與常見形式相同的異常高壓氣藏不穩(wěn)定滲流微分方程：

在圓柱坐標(biāo)系（r、θ、z）中，由于地層均質(zhì)，壓力p與θ角度無(wú)關(guān)，可得到柱坐標(biāo)形式下的滲流微分方程：

3 異常高壓裂縫性氣藏變系數(shù)滲流模型數(shù)值方法研究

前面我們已經(jīng)推導(dǎo)出了異常高壓氣藏單相非線性滲流方程，結(jié)合相應(yīng)的初邊值條件，如果僅考慮二維平面圓形封閉邊界問(wèn)題，模型可簡(jiǎn)化為：

3.1 網(wǎng)格的選擇

根據(jù)滲流力學(xué)理論可知在井底附近存在壓降漏斗，井筒附近地層壓力變化非常大，遠(yuǎn)離井筒區(qū)域壓力變化小。因此，采用非等距網(wǎng)格進(jìn)行定義域剖分，即沿徑向向外，網(wǎng)格由密逐漸變稀，以達(dá)到所需的計(jì)算精度。為使計(jì)算結(jié)果具有更佳的穩(wěn)定性，我們?cè)谶@里采用對(duì)數(shù)等距點(diǎn)中心網(wǎng)格進(jìn)行劃分。

3.2 差分方程的形成

根據(jù)前文分析，孔隙度、滲透率隨有效壓力呈指數(shù)式變化規(guī)律，而氣體在開發(fā)過(guò)程中隨壓力減小而發(fā)生膨脹，其壓縮系數(shù)變化趨勢(shì)與巖石壓縮系數(shù)相反，是隨有效壓力增大而增大的，假設(shè)其變化也遵循指數(shù)式變化的規(guī)律：

引入差分記號(hào)，再代入傳導(dǎo)系數(shù)的邊界，最終可得：

簡(jiǎn)記Pn＋1為Pi，引入記號(hào)后有：

式中：

3.3 定解條件的處理

·初始條件：

作為初始層上的沿空間各節(jié)點(diǎn)處的函數(shù)值，即

·內(nèi)邊界條件：

·外邊界條件：

可得：

3.4 非線性方程組

綜合式（13）、（18）、（20）即可得到由N個(gè)方程構(gòu)成的N元非線性方程組：

于是，對(duì)（14）式模型的求解，就轉(zhuǎn)化成了對(duì)非線性方程組（21）的逐層求解。

3.5 自適應(yīng)延拓法

下面我們討論對(duì)（21）式表示的方程組進(jìn)行求解的算法。設(shè)F：D?Rn→Rn，含有n個(gè)自變量的N元非線性方程組問(wèn)題F（x）＝0的求解可選取常用的Newton法、Brown法、割線法、擬牛頓法等。而這些方法都具有局部收斂性，在使用它們進(jìn)行計(jì)算時(shí)要求迭代初值x0與解x*充分接近，因此實(shí)際使用時(shí)往往很難滿足。延拓法是一種可以擴(kuò)大局部收斂區(qū)域的算法，它在計(jì)算非線性方程組的迭代初值時(shí)是一種很有效的方法。在運(yùn)用延拓法尋求得到充分靠近真解的初值后，還可以和其他局部收斂的迭代算法復(fù)合起來(lái)使用，以求得更加精確的近似數(shù)值解。

在使用延拓法求F（x）＝0的解x*∈D時(shí)，先要人為地引入一個(gè)參數(shù)t，構(gòu)造函數(shù)族H（x，t），使之滿足同倫條件：

其中F0（x）＝0的解已知。

構(gòu)造H（x，t）的方法有很多種，但并非任意構(gòu)造H（x，t）的都有效。在此，根據(jù)方程的具體形式，我們?nèi)?/p>

引進(jìn)一個(gè)控制參數(shù)對(duì)普通延拓法進(jìn)行改進(jìn)，使其除了保持一般延拓法的大范圍收斂性外，還可通過(guò)改變參數(shù)，重新定義同倫，以避免Hx（x（t），a（t），t）出現(xiàn)奇異。在進(jìn)行具體的計(jì)算機(jī)數(shù)值求解時(shí)，可自動(dòng)控制非奇異參數(shù)解決奇異問(wèn)題，并能自動(dòng)選擇步長(zhǎng)以及判斷是否進(jìn)入Newton迭代收斂域，這種改進(jìn)的延拓法就是自適應(yīng)延拓法。

將（23）式中的常量a參數(shù)化為a（t），從而構(gòu)造新的同倫：

這里，a（t）是一個(gè)含參的已知函數(shù)，在tε［0，1］上連續(xù)可微，比如a（t）＝a（1＋t）。

顯然，對(duì)于（24），當(dāng)時(shí)t＝0，時(shí)H（x0，0，a（0））＝F（x0），當(dāng)t＝1時(shí)，H（x，1，a（1））＝F（x），故H（x（t），a（t），t）滿足同倫條件（22）。

為了求解x（1），在（23）的兩端對(duì)t求導(dǎo)并整理得：

其中g(shù)（t）＝a（t）（1－t）。在t?［0，1］上，由于a（t）的連續(xù)可微性，g（t）在上t?［0，1］也是連續(xù)可微的。

于是，與（24）等價(jià)的Cauchy問(wèn)題是

容易得出，只要g（t）取得足夠大，總能使矩［F′（x）＋g（t）I］陣非奇異。而改變a（t），就可以重新定義H（x（t），a（t），t）。

首先采用Euler法進(jìn)行迭代初值的近似：

進(jìn)一步精確化則采用Newton迭代：

在計(jì)算F′（x）時(shí)采用差商近似J（x，δ），精確化則采用Newton迭代法，在Euler公式中結(jié)合使用1次2階的中點(diǎn)求積公式，于是編程進(jìn)行數(shù)值計(jì)算的迭代算法可以寫成：

式中

ej＝（0，…1，…，0）T為第j個(gè)分量為1的單位坐標(biāo)向量，i，j＝1，2，……n。

4 異常高壓氣藏單相非線性滲流模型算例

本文用理論數(shù)據(jù)模擬了某異常高壓氣藏。設(shè)原始地層壓力pi＝70 MPa，原始孔隙度10%，初始滲透率0.03μm2，原始?xì)怏w壓縮系數(shù)，4×10－3MPa－1，αΦ＝0.0008，αK＝0.06，αc＝0.015，黏度為0.35 mPa·s，氣層有效厚度20 m，產(chǎn)量q＝5×104m3，，井筒半徑rw＝0.1 m，地層半徑re＝1000 m，根據(jù)換算得初始擬壓力Pi＝1120。構(gòu)造非等距網(wǎng)格如圖1所示：

圖1 有限差分非等距網(wǎng)格圖

得到有限差分方程組后，直接采用Newton法進(jìn)行迭代求解，取初值，3……N－1，無(wú)法得到收斂值。采用自適應(yīng)延拓法，求得當(dāng)生產(chǎn)進(jìn)行300天后，擬壓力及壓力沿徑向距離的分布。

圖2 300天后擬壓力隨徑向距離變化趨勢(shì)圖

圖3 300天后壓力隨徑向距離變化趨勢(shì)圖

由圖中可以看出，在同一時(shí)間內(nèi)，距離氣井遠(yuǎn)近不同的點(diǎn)，其壓力變化是不一樣的，離井越近，壓力下降越快，離井越遠(yuǎn)，變化相對(duì)越平緩。

圖4 擬壓力隨時(shí)間變化趨勢(shì)圖

圖5 壓力隨時(shí)間變化趨勢(shì)圖

圖4、圖5是根據(jù)程序求得的井底附近擬壓力及對(duì)應(yīng)的壓力隨時(shí)間變化的趨勢(shì)圖，由圖中可以看出，擬壓力和壓力變化的趨勢(shì)是一致的。與低壓和常壓氣藏不同，在異常高壓氣藏開發(fā)初期，壓力下降較中后期緩慢，這可以解釋為由于儲(chǔ)層巖石具有較高的壓縮性，隨地層壓力的下降，巖石顆粒膨脹使壓力得以保持。到了氣藏開發(fā)后期，巖石壓縮系數(shù)降低，氣體壓縮系數(shù)增大，此時(shí)巖石壓縮系數(shù)相對(duì)于氣體壓縮系數(shù)已很小，氣產(chǎn)量主要由天然氣膨脹貢獻(xiàn)，所以氣藏的動(dòng)態(tài)非常接近于定容氣藏。異常高壓氣藏的這個(gè)特點(diǎn)，在開發(fā)動(dòng)態(tài)預(yù)測(cè)及原始儲(chǔ)量計(jì)算時(shí)應(yīng)給予高度重視。

5 結(jié)論

在進(jìn)行異常高壓氣藏研究時(shí)應(yīng)充分考慮儲(chǔ)層應(yīng)力敏感性的影響，掌握儲(chǔ)層物性參數(shù)隨有效應(yīng)力改變的變化規(guī)律。本文建立了滲透率、孔隙度等參數(shù)均隨有效應(yīng)力變化的異常高壓氣藏單相非線性滲流模型；為能準(zhǔn)確地計(jì)算井底附近壓力分布，采用了對(duì)數(shù)等距網(wǎng)格進(jìn)行網(wǎng)格剖分，即沿徑向向外，網(wǎng)格由密逐漸變稀；采用有限差分法對(duì)模型進(jìn)行離散求解，并創(chuàng)新采用自適應(yīng)延拓法，克服了傳統(tǒng)數(shù)值方法要求迭代初值與真解充分接近的局限性。編程計(jì)算結(jié)果表明，考慮了滲透率、孔隙度等參數(shù)隨有效應(yīng)力變化規(guī)律的異常高壓氣藏單相非線性滲流模型是符合工程實(shí)際的，自適應(yīng)延拓法的應(yīng)用對(duì)解決初值問(wèn)題及奇異問(wèn)題也是簡(jiǎn)單有效的。

符號(hào)注釋：

Φ—當(dāng)前壓力下的孔隙度，%；

Φ0—初始孔隙度，%；

αΦ—孔隙度變化系數(shù)；

Δp—有效覆壓，MPa；

K—目前壓力下的滲透率，μm2；

K0—初始滲透率，；

αK—滲透率變化系數(shù)；

μ—天然氣黏度，mPags；

p—壓力，MPa；

T—絕對(duì)溫度，K；

V—?dú)怏w體積，；

R—?dú)怏w常數(shù)，0.008314 MPa·m3／kmol·()K；

ρ—?dú)怏w密度，kg／m3；

Z—?dú)怏w偏差因子；

t—時(shí)間，d；

P—距離井r處在t時(shí)刻的擬壓力；

P0—初始擬壓力；

q—?dú)饩a(chǎn)量，104m3／d；

h—?dú)鈱雍穸龋琺；

rw—?dú)饩霃剑琺；

re—供給半徑，m。

1 Geertsma J.Theeffectof pressure decline on volumetric changes of porous rocks［J］.Petroleum Transport AME，1957，210：331－3402 Vairogs J，Hearn CL，Dareing DW，Rhoades VW.Effectof rock stress on gas production from low－permeability from reservoirs［J］.Journal of Petroleum Technology，1971，9：1161－1167

3 王江，王玉英.異常高壓、特低滲透油藏儲(chǔ)層壓力敏感性研究［J］.大慶石油地質(zhì)與開發(fā)，2003，22（5）：28－31.

4 黃繼新，彭仕宓，黃述旺，等.異常高壓氣藏儲(chǔ)層應(yīng)力敏感性研究［J］.西安石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2005，20（4）：21－25.

5 葛家理.油氣層滲流力學(xué)［M］.北京：石油工業(yè)出版社，1982.

（修改回稿日期 2013－07－10 編輯 文敏）

鄭琴，女，1981年出生，2011年獲石油工程計(jì)算技術(shù)專業(yè)博士學(xué)位；現(xiàn)在大港油田從事博士后研究。主要研究方向?yàn)橛蜌獠財(cái)?shù)值模擬及優(yōu)化控制。地址：（300280）天津市大港區(qū)大港油田勘探開發(fā)研究院歧口開發(fā)一室。電話：18622042029。E－mail：zhengqin710＠126.com