重力作用下軸對稱非穩(wěn)態(tài)滲流壓強(qiáng)場的解析解

劉智穎,彭文耀,章成廣*

(1.油氣資源與勘探技術(shù)教育部重點實驗室(長江大學(xué)),2.長江大學(xué)地球物理與石油資源學(xué)院,湖北 武漢 430100;3.江漢石油管理局有限公司測錄井工程公司,湖北 潛江 433123)

通過注水來增大儲層內(nèi)的壓強(qiáng)是各大油田常用的增產(chǎn)或穩(wěn)產(chǎn)的手段．注水進(jìn)入地層使地層中的注水壓強(qiáng)增大,從而驅(qū)替地層中的油氣進(jìn)入生產(chǎn)井以提高采收率．為了準(zhǔn)確描述地層中的注水壓強(qiáng)分布狀況,從而計算滲流速度并預(yù)測注水效果和生產(chǎn)井中的產(chǎn)量,需要建立滲流力學(xué)方程并求解出不同條件下的滲流場注水壓強(qiáng)分布函數(shù)．

前人進(jìn)行過許多相關(guān)研究,主要研究內(nèi)容及成果可分兩類．一類是滲流力學(xué)方程的建立,包括單相線性滲流方程、兩相線性滲流方程、稠油非線性滲流方程[1]等,主要是給定邊界條件和初始條件,求解這些方程即可得出地層中注水壓強(qiáng)分布函數(shù)[2-3]．另一類是滲流力學(xué)方程的求解方法,該方法可分為數(shù)值解法和解析解法．?dāng)?shù)值解法的研究成果較多,如用有限差分法計算公路隧道附近區(qū)域的滲流場[4]、水庫大壩中的滲流場[5]以及用有限元法求解低滲油藏開采過程中注水壓強(qiáng)分布狀況[6-7],人工壓裂后地層流體壓強(qiáng)分布狀況[8]等．由于絕大多數(shù)滲流力學(xué)方程的定解問題無法解析求解,且解析解的形式大都比較復(fù)雜,因而解析解法的研究成果較少．盡管如此,解析解不僅是數(shù)值解的重要理論基礎(chǔ)之一,而且計算程序比較容易實現(xiàn),因而也是常見的研究課題．主要成果有平面徑向注水壓強(qiáng)分布函數(shù)[9-10]、柱坐標(biāo)下非達(dá)西滲流力學(xué)方程的解析解[11]以及非線性三參數(shù)模型[12-14]等．這些解析解主要是基于平面徑向穩(wěn)定滲流模型求出的．該模型的基本假定為:徑向上供水邊緣的壓強(qiáng)最大,生產(chǎn)井內(nèi)壓強(qiáng)最小;縱向上壓強(qiáng)不變,且地層中注水壓強(qiáng)分布狀況穩(wěn)定,不隨時間變化．因此,該模型不僅未考慮重力作用下注水壓強(qiáng)隨深度增大而增大的特點,也未反映注水壓強(qiáng)隨時間變化的過程,因而難以準(zhǔn)確描述驅(qū)替過程中注水壓強(qiáng)的分布狀況[9-10]．事實上,由于注入水受到重力的作用,供水邊緣的注水壓強(qiáng)隨深度增大而增大;此外,滲流過程中注水壓強(qiáng)分布狀況會不斷發(fā)生變化．這種情況下,方程的供水邊緣的壓強(qiáng)邊界條件是一個與深度和縱向壓強(qiáng)梯度有關(guān)的函數(shù),屬于非齊次邊界條件,而且方程還包含時間項,使得求解過程更加復(fù)雜．針對該問題,本研究采用分離變量法將非齊次問題轉(zhuǎn)化為齊次問題后再求方程的解析解,并將其應(yīng)用于石油開采過程中注水后的產(chǎn)液量預(yù)測,驗證所求解析解的正確性．

1　注水壓強(qiáng)分布函數(shù)滿足的微分方程

根據(jù)文獻(xiàn)[10],平面徑向注水壓強(qiáng)分布函數(shù)p與井中秒產(chǎn)量Q分別為:

(1)

(2)

其中,pw是井壁的壓強(qiáng),pe是供水邊緣的壓強(qiáng),r0是供水邊緣到生產(chǎn)井井軸的距離,rw是生產(chǎn)井半徑,K是地層滲透率,h是地層厚度,μ是流體黏度．

為研究注水壓強(qiáng)分布狀況隨時間的變化過程,下面以地層上界面為原點,豎直向下的方向為z軸(圖1),建立二維柱坐標(biāo)系,求解二維軸對稱非穩(wěn)態(tài)單相滲流壓強(qiáng)場方程(亦稱輸運方程)[10,15-16]:

(3)

其中,r、z分別是極徑坐標(biāo)和高度坐標(biāo),t是時間,κ=K/μCt為地層導(dǎo)壓系數(shù)[10],Ct為地層綜合壓縮系數(shù)[10]．

假設(shè)目的層是滲透性地層,上下圍巖為非滲透地層,則層界面附近滲流速度垂直于層界面的分量為零．由于滲流速度往往很小,可以認(rèn)為供水邊緣的注水壓強(qiáng)等于靜水壓強(qiáng)．這樣,方程(3)的邊界條件和初始條件分別為

供水邊緣的壓強(qiáng)邊界:p|r=r0=pe+ρgz,

(4)

生產(chǎn)井處的邊界:p|r=rw=pw,

(5)

(6)

地層中初始壓強(qiáng)值:p|t=0=p0,

(7)

其中,ρ是注入水密度,g是重力加速度,p0為地層中的初始壓強(qiáng)值,可近似看作常數(shù)．

2　注水壓強(qiáng)分布函數(shù)的求解過程

方程(3)可用分離變量法求解,但若使用分離變量法則會因徑向非齊次邊界條件式(4)和(5)而無法求解．因此,需要先將供水邊緣的非齊次邊界條件化為齊次邊界條件再求解方程(3)．

2.1　供水邊緣非齊次邊界條件的處理方法

令p=I(r,z)+K(r,z,t),其中I(r,z)是一個滿足邊界條件式(4)和(5)的調(diào)和函數(shù),則K(r,z,t)滿足式(3)．然而,經(jīng)過這種簡單變型后雖然可使函數(shù)K(r,z,t)在r=r0處的邊界條件轉(zhuǎn)化為齊次的,但在z=±h處的邊界條件卻變?yōu)榉驱R次的,仍然難以求解．因此,需要進(jìn)一步將函數(shù)K(r,z,t)分解為U(r,z)、V(r,z)和H(r,z,t)之和,使得U(r,z)和V(r,z)均滿足二維軸對稱拉普拉斯方程,并保證z方向至少有一個邊界條件是齊次的;H(r,z,t)滿足輸運方程,且邊界條件均為齊次的．這樣,上述3個函數(shù)滿足的定解問題均可以用分離變量法求解析解．具體的變型結(jié)果如下:

p(r,z,t)=I(r,z)+U(r,z)+

V(r,z)+H(r,z,t),

(8)

這樣,U(r,z)、V(r,z)和H(r,z,t)分別滿足

(9)

(10)

(11)

不難證明,函數(shù)U(r,z)、V(r,z)和H(r,z,t)之和滿足的邊界條件與方程(3)的邊界條件等價．此外,由于這3個函數(shù)值均有限且與時間無關(guān),因此只要r滿足rw≤r≤r0,就有

p|t=0=H|t=0+U(r,z)+V(r,z)+

I(r,z)≡p0,

因此,經(jīng)上述變型后,式(11)的初始條件與式(7)等價．

2.2　方程的求解步驟

2.2.1求解U(r，z)與V(r，z)之和

U(r,z)和V(r,z)均滿足拉普拉斯方程．由分離變量法,可得到其通解[17]．這里僅寫出U(r,z)的求解過程,V(r,z)的求解過程與之類似,不做贅述．

由于U(r,z)在r方向上滿足第一類齊次邊界條件,因此其通解形式為

(12)

將U(r,z)在z方向上的邊界條件代入式(12)得

(13)

(14)

(15)

由式(13)可得Bn=An．因此,式(14)可簡化為

(16)

(17)

將式(17)代入式(12)即可確定U(r,z),同理可求出V(r,z)．

令W(r,z)=U(r,z)+V(r,z),則有

(18)

其中，

2.2.2求解H(r,z,t)

式(11)同樣可以采用分離變量法求解．令H(r,z,t)=R(r)Z(z)T(t),則式(11)可分解為關(guān)于R(r)、Z(z)和T(t)的常微分方程．不妨設(shè)α和β分別為函數(shù)T(t)和R(r)的本征值．通過分析α和β的取值范圍可知,當(dāng)且僅當(dāng)0<β≤α?xí)r,式(11)存在滿足邊界條件的解．這樣,式(11)的通解形式為

(19)

Bm0≡0,

(20)

(21)

將式(20)和(21)代入式(19)得

(22)

將式(11)的初始條件代入式(22),可得

H(r,z,0)=p0-I(r,z)-W(r,z),

(23)

(24)

(25)

(26)

這樣

(27)

其中,

將式(27)代入式(24)和(25)得

(28)

(29)

將式(28)和(29)代入式(22)即得H(r,z,t),故P(r,z,t)也隨即確定．

函數(shù)P(r,z,t)需要用程序進(jìn)行計算．下面給出計算程序中涉及的變量及相關(guān)算法的實現(xiàn)過程．程序?qū)方向和z方向的長度r0和h分別設(shè)定為700和20 m．W(r,z)函數(shù)的級數(shù)項數(shù)n取60,H(r,z,t)函數(shù)的級數(shù)項數(shù)m取130、q取50．式(11)中μ值采用精度較高的12階高斯-洛巴托公式進(jìn)行計算,該公式中勒讓德多項式的導(dǎo)數(shù)采用步長為5×10-3的4步理查森外推公式進(jìn)行計算．

2.3　求解方法的推廣

根據(jù)以上推導(dǎo)過程,本研究提出了一種徑向第一類非齊次邊界條件和縱向第二類齊次邊界條件下的二維軸對稱輸運方程的解析解法．為不失一般性,現(xiàn)將方程(3)中的p改作u．這樣,方程(3)及其邊界條件和初始條件分別為

(30)

u|r=r0=ξ(z)，u|r=r1=η(z)，

其中,r0>r1．

求解步驟如下:

3) 用分離變量法求出v1(r,z)、v2(r,z)和w(r,z,t)，即得式(30)的解．

需要指出,若ξ(z)和η(z)均為z的一次函數(shù),即ξ(z)=a0+b0z，η(z)=a1+b1z，則f(r,z)可選取

f(r,z)=a1+b1z+

(31)

式(31)的意義在于,若同時考慮供水邊緣和生產(chǎn)井井壁處的z方向的壓強(qiáng)梯度,則利用上述方法同樣可以求出地層中的注水壓強(qiáng)分布函數(shù)．

圖2　不考慮供水邊緣縱向壓強(qiáng)梯度時的注水壓強(qiáng)分布狀況Fig.2 The distribution of injection pressure while ignoring the vertical pressure gradient at the injection edge

3　求解結(jié)果的檢驗及其應(yīng)用

3.1　求解結(jié)果的檢驗

檢驗求解結(jié)果正確性的最簡單方法是證明當(dāng)g=0且t→∞時,式(8)與式(1)等價，即若不考慮重力及時間,檢驗式(8)是否與式(1)吻合．不難發(fā)現(xiàn),當(dāng)g=0時,I(r,z)退化為式(1),且W(r,z)≡0;當(dāng)t→∞時,有H(r,z,t)≡0．因此,當(dāng)g=0且t→∞時P(r,z,t)由較復(fù)雜的級數(shù)式退化為式(1),說明式(8)的求解方法及結(jié)果正確．

本研究根據(jù)程序計算結(jié)果繪制了兩種條件下注水壓強(qiáng)分布圖以及注水壓強(qiáng)曲線．供水邊緣壓強(qiáng)pe設(shè)為1.2 MPa;生產(chǎn)井井眼半徑rw設(shè)為0.101 6 m;井壁壓強(qiáng)pw設(shè)為1.0 MPa．圖2給出了不考慮供水邊緣的縱向壓強(qiáng)梯度,也就是g=0條件下的注水壓強(qiáng)分布圖．圖2(a)和(b)分別為注水后第15天和注水后第3個月時地層中的注水壓強(qiáng)分布圖．其中,圖2(a)為r> 670 m范圍內(nèi)的注水壓強(qiáng)分布情況．由圖可見,注水壓強(qiáng)隨r的減小而減小．圖2(a)中較高的壓強(qiáng)主要集中于r> 685 m的區(qū)域,而圖2(b)中集中于r> 64 m的區(qū)域．說明注水后,注水壓強(qiáng)較高的區(qū)域逐漸接近生產(chǎn)井．圖2(c)為不同時刻的注水壓強(qiáng)與r的關(guān)系曲線．可以看出,注水第15天時注水壓強(qiáng)不為零的區(qū)域主要在r> 500 m的區(qū)域;第90天時,注水壓強(qiáng)分布狀況趨于穩(wěn)定且與式(1)的計算結(jié)果基本吻合．

下部壓強(qiáng)隨z坐標(biāo)值增大,地層導(dǎo)壓系數(shù)κ設(shè)為100 cm2/s,注入水密度ρ設(shè)為1.0 t/m3,重力加速度g設(shè)為9.8 m/s2,地層厚度h設(shè)為20 m,其余參數(shù)與圖(2)中相同．圖3中給出了考慮供水邊緣的縱向壓強(qiáng)梯度．可以看出,地層底部的注水壓強(qiáng)明顯高于頂部,說明地層底部因重力作用而產(chǎn)生了較高的注水壓強(qiáng)．然而,隨著時間的推移,r> 685 m,z< 4 m的區(qū)域有一個注水壓強(qiáng)較低的區(qū)域(稱壓強(qiáng)異常區(qū))(圖3(b))．不難發(fā)現(xiàn),該區(qū)域的注水壓強(qiáng)隨r的減小而增大．這與實際情況不符．產(chǎn)生該問題的原因是本研究僅在供水邊緣的邊界條件(式7)中考慮了重力的影響,而方程(3)本身卻并未考慮重力．因此,當(dāng)注水壓強(qiáng)分布逐漸趨于穩(wěn)定時,供水邊緣底部高壓強(qiáng)區(qū)域的注入水會流至供水邊緣頂部,從而形成了壓強(qiáng)異常區(qū)．盡管如此,這對式(8)的實際應(yīng)用效果影響不大．因為,式(8)主要應(yīng)用于計算生產(chǎn)井中的產(chǎn)量．根據(jù)圖3(c)和(d),壓強(qiáng)異常區(qū)主要影響供水邊緣附近的注水壓強(qiáng),對生產(chǎn)井附近的注水壓強(qiáng)分布狀況影響較小．

通過圖2(c)與圖3(d)的對比看出,圖3(d)中注水壓強(qiáng)與柱坐標(biāo)系極徑的關(guān)系曲線形狀與圖2(c)中曲線相似．兩者的不同之處在于,圖3(d)表明當(dāng)z=20 m時,注水壓強(qiáng)曲線值以及r→0時的曲線斜率均高于z=10 m和z=0,亦即地層底部的注水壓強(qiáng)及壓強(qiáng)梯度均較大,因而當(dāng)?shù)貙又械臐B透率處處相等時,底部的滲流速度較快．

圖3　考慮供水邊緣縱向壓強(qiáng)梯度時的注水壓強(qiáng)分布狀況Fig.3 The distribution of injection pressure while considering the vertical pressure gradient at the injection edge

3.2　應(yīng)用方法及效果

應(yīng)用式(8)可以計算生產(chǎn)井中不同時刻的日產(chǎn)量．計算日產(chǎn)量的步驟如下:

1) 將式(8)代入達(dá)西公式求出生產(chǎn)井井壁處的滲流速度v．

2) 將v代入產(chǎn)量的計算式[10]求出秒產(chǎn)量Q,日產(chǎn)量也隨即確定．

表1　江漢油田Hxx3試油層段數(shù)據(jù)

Tab.1　Data of testing interval of well Hxx3 in Jianghan oilfield

K/(10-3μm2)μ/(mPa·s)pe/MPap0/MPapw/MParw/cmCt/(10-10Pa-1)r0/mρ/(t·m-3)13.3615.314.3711.711.710.161.777401.06

下面介紹式(8)在江漢油田Hxx3井中的應(yīng)用效果．該井的1 317.5～1 326.6 m層段的測井結(jié)論為油層,巖性為泥質(zhì)砂巖,屬于低滲地層需進(jìn)行注水開采．注水后該井段的日產(chǎn)量隨時間變化,其各時刻的試油資料為:注水后第15天,日產(chǎn)油0.52 m3/d;注水后第30天,日產(chǎn)油3.12 m3/d;注水后第60天,日產(chǎn)油6.36 m3/d;注水后第90天,日產(chǎn)油6.04 m3/d．注水層段中的滲流場可近似看作單相滲流場,與式(8)的適用條件基本吻合．詳細(xì)數(shù)據(jù)見表1．

本研究應(yīng)用式(8)繪制了該層段的注水壓強(qiáng)分布圖及井內(nèi)流速曲線(圖4),并將式(2)和(8)計算的生產(chǎn)井中的日產(chǎn)量與表1中的試油資料進(jìn)行了對比(圖5)．

從圖4(a)和(b)中不難看出,注水90 d后,注水壓強(qiáng)波傳播到生產(chǎn)井處,且注水壓強(qiáng)分布狀況已趨于穩(wěn)定．圖4(c)是注水90 d后地層中注水壓強(qiáng)p與柱坐標(biāo)系極徑r的關(guān)系曲線圖．圖中,虛線是將供水邊緣的頂部壓強(qiáng)pe代入式(2)的計算結(jié)果,實線是不同深度條件下式(8)的計算結(jié)果．圖4(d)是生產(chǎn)井中不同深度不同時刻的滲流速度曲線．從圖4(c)和(d)中可以看出,距離生產(chǎn)井井軸距離相同的位置,壓強(qiáng)隨深度增大而增大,且在r=rw附近壓強(qiáng)梯度也隨深度增大而增大．說明只要地層滲透率均勻分布,則地層中滲流速度將隨深度增大而增大．

這與圖3(d)的結(jié)論一致,同時還可以看出,式(1)計算的注水壓強(qiáng)值偏小．

圖4　Hxx3井的層內(nèi)注水壓力分布狀況及滲流速度Fig.4 Formation injection pressure and the seepage velocity in well Hxx3

圖5為兩種方法計算的日產(chǎn)量與實際試油結(jié)果的對比圖．通過對比不難看出,式(8)計算的日產(chǎn)量及其隨時間變化的規(guī)律與實際試油結(jié)果比較一致,說明式(8)比式(2)更能準(zhǔn)確計算出生產(chǎn)井中的日產(chǎn)量．此外,還能看出90 d之內(nèi),日產(chǎn)量逐漸增長并趨于穩(wěn)定,說明注水開始后,地層中壓強(qiáng)波的傳播過程需要時間,因而生產(chǎn)井中的日產(chǎn)量并不會迅速提高．

圖5　兩種方法計算的日產(chǎn)量與實際試油結(jié)果的對比圖Fig.5 The comparison between day production calculated by capillary pressure formula and by practical production testing

4　結(jié)　論

本研究通過求解重力作用下軸對稱非穩(wěn)態(tài)滲流壓強(qiáng)場方程的解析解,提出了一種徑向第一類非齊次邊界條件和縱向第二類齊次邊界條件下的二維軸對稱輸運方程的一般求解方法,并獲得了如下結(jié)論:

1) 注水開采時,地層中的注水壓強(qiáng)及注水壓強(qiáng)梯度均隨深度增大而增大．對于滲透率分布均勻的地層,其底部的滲流速度更快．

2) 隨著壓強(qiáng)波的傳播,地層中的注水壓強(qiáng)分布狀況逐漸趨于穩(wěn)定．在該過程中,生產(chǎn)井內(nèi)的日產(chǎn)量逐漸上升至最大值．

3) 由于地層底部的注水壓強(qiáng)梯度較大,因此,如果注水之后還需要通過壓裂、酸化等措施提高地層滲透率,則應(yīng)該優(yōu)先針對地層頂部實施．

4) 僅在邊界條件中考慮了重力的作用,而式(3)本身并未反映重力的影響,今后將針對該問題對式(3)進(jìn)行改進(jìn),并提出新的求解方法．

參考文獻(xiàn)：

[1]姚同玉,黃廷章,李繼山.孔隙介質(zhì)中稠油流體非線性滲流方程[J].力學(xué)學(xué)報,2012,44(1):106-110.

[2]林振生.兩類橢圓偏微分方程解的存在性問題[J].福建師范大學(xué)學(xué)報(自然科學(xué)版),2010,4(1):1-33.

[3]沈欽銳,程立新,周宗福.一類三階多時滯p-Laplacian方程周期解的存在性[J].武漢大學(xué)學(xué)報(理學(xué)版),2013,59(6):505-510.

[4]嚴(yán)紹洋,李亮輝,高燕希,等.公路隧道開挖滲流場的有限差分法分析[J].中外公路,2007,27(6):120-123.

[5]杜澤.基于有限差分計算的某水庫大壩滲流場穩(wěn)定性分析[J].安徽農(nóng)業(yè)科學(xué),2016,44(3):309-311.

[6]楊正明,于榮澤.特低滲透油藏非線性滲流數(shù)值模擬[J].石油勘探與開發(fā),2010，37(1):35-37.

[7]周志軍.低滲透儲層流固耦合滲流理論及應(yīng)用研究[D].大慶:大慶石油學(xué)院,2003:76-82.

[8]劉振宇.有限元法在油藏滲流中的理論和應(yīng)用[D].大慶:大慶石油學(xué)院,2003:44-65.

[9]王君連.用勢函數(shù)法計算滲流干擾區(qū)的地下水流態(tài)[J].水利水電技術(shù),1999,30(2):19-25.

[10]張建國,杜殿發(fā),侯健,等.油氣層滲流力學(xué)[M].2版.北京:中國石油大學(xué)出版社,2010:63-67.

[11]同登科,葛家理.分形油藏非達(dá)西低速滲流模型及其解[J].大慶石油地質(zhì)與開發(fā),1996,15(3):18-23.

[12]李宏偉.帶啟動壓力梯度的滲流方程求解[D].長春:東北石油大學(xué),2013:18-21.

[13]鄧英爾,劉慈群.低滲油藏非線性滲流規(guī)律數(shù)學(xué)模型及其應(yīng)用[J].石油學(xué)報,2001,22(4):72-77.

[14]廖作才.非線性滲流方程解析方程研究及應(yīng)用[D].北京:中國科學(xué)院滲流力學(xué)研究所,2015:15-22.

[15]王曉冬,侯曉春,郝明強(qiáng),等.低滲透介質(zhì)有啟動壓力梯度的非穩(wěn)態(tài)壓力分析[J].石油學(xué)報,2011,32(5):847-851.

[16]石達(dá)友.超稠油油藏水平井產(chǎn)量變化規(guī)律研究[D].青島:中國石油大學(xué)(華東),2009:33-42.

[17]梁昆淼.數(shù)學(xué)物理方法 [M].3版.北京:高等教育出版社,1998:261-263.