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

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

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

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

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

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)分布狀況隨時(shí)間的變化過程,下面以地層上界面為原點(diǎn),豎直向下的方向?yàn)閦軸(圖1),建立二維柱坐標(biāo)系,求解二維軸對(duì)稱非穩(wěn)態(tài)單相滲流壓強(qiáng)場(chǎng)方程(亦稱輸運(yùn)方程)[10,15-16]:

(3)

其中,r、z分別是極徑坐標(biāo)和高度坐標(biāo),t是時(shí)間,κ=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)

地層中初始?jí)簭?qiáng)值:p|t=0=p0,

(7)

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

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

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

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

令p=I(r,z)+K(r,z,t),其中I(r,z)是一個(gè)滿足邊界條件式(4)和(5)的調(diào)和函數(shù),則K(r,z,t)滿足式(3)．然而,經(jīng)過這種簡(jiǎn)單變型后雖然可使函數(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)均滿足二維軸對(duì)稱拉普拉斯方程,并保證z方向至少有一個(gè)邊界條件是齊次的;H(r,z,t)滿足輸運(yùn)方程,且邊界條件均為齊次的．這樣,上述3個(gè)函數(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)的邊界條件等價(jià)．此外,由于這3個(gè)函數(shù)值均有限且與時(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)等價(jià)．

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)可簡(jiǎn)化為

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

2.3　求解方法的推廣

根據(jù)以上推導(dǎo)過程,本研究提出了一種徑向第一類非齊次邊界條件和縱向第二類齊次邊界條件下的二維軸對(duì)稱輸運(yùn)方程的解析解法．為不失一般性,現(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)的意義在于,若同時(shí)考慮供水邊緣和生產(chǎn)井井壁處的z方向的壓強(qiáng)梯度,則利用上述方法同樣可以求出地層中的注水壓強(qiáng)分布函數(shù)．

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

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

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

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

本研究根據(jù)程序計(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個(gè)月時(shí)地層中的注水壓強(qiáng)分布圖．其中,圖2(a)為r> 670 m范圍內(nèi)的注水壓強(qiáng)分布情況．由圖可見,注水壓強(qiáng)隨r的減小而減?。畧D2(a)中較高的壓強(qiáng)主要集中于r> 685 m的區(qū)域,而圖2(b)中集中于r> 64 m的區(qū)域．說明注水后,注水壓強(qiáng)較高的區(qū)域逐漸接近生產(chǎn)井．圖2(c)為不同時(shí)刻的注水壓強(qiáng)與r的關(guān)系曲線．可以看出,注水第15天時(shí)注水壓強(qiáng)不為零的區(qū)域主要在r> 500 m的區(qū)域;第90天時(shí),注水壓強(qiáng)分布狀況趨于穩(wěn)定且與式(1)的計(jì)算結(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)．然而,隨著時(shí)間的推移,r> 685 m,z< 4 m的區(qū)域有一個(gè)注水壓強(qiáng)較低的區(qū)域(稱壓強(qiáng)異常區(qū))(圖3(b))．不難發(fā)現(xiàn),該區(qū)域的注水壓強(qiáng)隨r的減小而增大．這與實(shí)際情況不符．產(chǎn)生該問題的原因是本研究?jī)H在供水邊緣的邊界條件(式7)中考慮了重力的影響,而方程(3)本身卻并未考慮重力．因此,當(dāng)注水壓強(qiáng)分布逐漸趨于穩(wěn)定時(shí),供水邊緣底部高壓強(qiáng)區(qū)域的注入水會(huì)流至供水邊緣頂部,從而形成了壓強(qiáng)異常區(qū)．盡管如此,這對(duì)式(8)的實(shí)際應(yīng)用效果影響不大．因?yàn)?式(8)主要應(yīng)用于計(jì)算生產(chǎn)井中的產(chǎn)量．根據(jù)圖3(c)和(d),壓強(qiáng)異常區(qū)主要影響供水邊緣附近的注水壓強(qiáng),對(duì)生產(chǎn)井附近的注水壓強(qiáng)分布狀況影響較?。?/p>

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

圖3　考慮供水邊緣縱向壓強(qiáng)梯度時(shí)的注水壓強(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)可以計(jì)算生產(chǎn)井中不同時(shí)刻的日產(chǎn)量．計(jì)算日產(chǎn)量的步驟如下:

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

2) 將v代入產(chǎn)量的計(jì)算式[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層段的測(cè)井結(jié)論為油層,巖性為泥質(zhì)砂巖,屬于低滲地層需進(jìn)行注水開采．注水后該井段的日產(chǎn)量隨時(shí)間變化,其各時(shí)刻的試油資料為:注水后第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．注水層段中的滲流場(chǎng)可近似看作單相滲流場(chǎng),與式(8)的適用條件基本吻合．詳細(xì)數(shù)據(jù)見表1．

本研究應(yīng)用式(8)繪制了該層段的注水壓強(qiáng)分布圖及井內(nèi)流速曲線(圖4),并將式(2)和(8)計(jì)算的生產(chǎn)井中的日產(chǎn)量與表1中的試油資料進(jìn)行了對(duì)比(圖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)的計(jì)算結(jié)果,實(shí)線是不同深度條件下式(8)的計(jì)算結(jié)果．圖4(d)是生產(chǎn)井中不同深度不同時(shí)刻的滲流速度曲線．從圖4(c)和(d)中可以看出,距離生產(chǎn)井井軸距離相同的位置,壓強(qiáng)隨深度增大而增大,且在r=rw附近壓強(qiáng)梯度也隨深度增大而增大．說明只要地層滲透率均勻分布,則地層中滲流速度將隨深度增大而增大．

這與圖3(d)的結(jié)論一致,同時(shí)還可以看出,式(1)計(jì)算的注水壓強(qiáng)值偏?。?/p>

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

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

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

4　結(jié)　論

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

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

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

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

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

參考文獻(xiàn)：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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