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

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

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

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

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

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

(1)

(2)

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

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

(3)

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

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

供水邊緣的壓強邊界:p|r=r0=pe+ρgz,

(4)

生產井處的邊界:p|r=rw=pw,

(5)

(6)

地層中初始壓強值:p|t=0=p0,

(7)

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

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

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

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

令p=I(r,z)+K(r,z,t),其中I(r,z)是一個滿足邊界條件式(4)和(5)的調和函數(shù),則K(r,z,t)滿足式(3)．然而,經(jīng)過這種簡單變型后雖然可使函數(shù)K(r,z,t)在r=r0處的邊界條件轉化為齊次的,但在z=±h處的邊界條件卻變?yōu)榉驱R次的,仍然難以求解．因此,需要進一步將函數(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ù)滿足的定解問題均可以用分離變量法求解析解．具體的變型結果如下:

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ù)值均有限且與時間無關,因此只要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)可分解為關于R(r)、Z(z)和T(t)的常微分方程．不妨設α和β分別為函數(shù)T(t)和R(r)的本征值．通過分析α和β的取值范圍可知,當且僅當0<β≤α時,式(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)需要用程序進行計算．下面給出計算程序中涉及的變量及相關算法的實現(xiàn)過程．程序將r方向和z方向的長度r0和h分別設定為700和20 m．W(r,z)函數(shù)的級數(shù)項數(shù)n取60,H(r,z,t)函數(shù)的級數(shù)項數(shù)m取130、q取50．式(11)中μ值采用精度較高的12階高斯-洛巴托公式進行計算,該公式中勒讓德多項式的導數(shù)采用步長為5×10-3的4步理查森外推公式進行計算．

2.3　求解方法的推廣

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

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

3　求解結果的檢驗及其應用

3.1　求解結果的檢驗

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

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

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

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

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

3.2　應用方法及效果

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

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

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

表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井中的應用效果．該井的1 317.5～1 326.6 m層段的測井結論為油層,巖性為泥質砂巖,屬于低滲地層需進行注水開采．注水后該井段的日產量隨時間變化,其各時刻的試油資料為:注水后第15天,日產油0.52 m3/d;注水后第30天,日產油3.12 m3/d;注水后第60天,日產油6.36 m3/d;注水后第90天,日產油6.04 m3/d．注水層段中的滲流場可近似看作單相滲流場,與式(8)的適用條件基本吻合．詳細數(shù)據(jù)見表1．

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

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

這與圖3(d)的結論一致,同時還可以看出,式(1)計算的注水壓強值偏?。?/p>

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

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

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

4　結　論

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

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

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

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

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

參考文獻：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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