張春光， 姜瑞忠， 喬 欣， 崔永正， 沈澤陽， 張福蕾， 原建偉

( 1. 中國石油大學(xué)(華東) 石油工程學(xué)院，山東 青島 266580； 2. 中國石油北京油氣調(diào)控中心，北京 100007 )

0 引言

研究雙重介質(zhì)致密分形氣藏的滲流特征，需要充分考慮致密基質(zhì)塊的非線性滲流、啟動壓力梯度及裂縫的應(yīng)力敏感性、分形特征等因素，應(yīng)用達(dá)西滲流理論描述該類儲層存在局限性[1-3]。

基于Warrant-Root模型，何國良[4]、同登科等[5]建立雙重介質(zhì)滲流模型，分析壓力對具有分形特征的孔隙度和滲透率的影響。宋付權(quán)[6]考慮啟動壓力梯度和介質(zhì)變形，推導(dǎo)低滲透油藏的產(chǎn)能公式并研究產(chǎn)量遞減規(guī)律。姜漢橋等[7]建立不完全可逆變形介質(zhì)的滲流模型，釆用全隱式和Richtmyer線性化方法對模型進(jìn)行數(shù)值求解。李允等[8]認(rèn)為低滲透氣藏存在啟動壓力梯度，在特定壓力范圍內(nèi)，滲流速度與壓力梯度呈非線性關(guān)系。孫賀東等[9]根據(jù)裂縫性應(yīng)力敏感氣藏的實驗數(shù)據(jù)和試井理論，引入滲透率模量，建立滲透率呈指數(shù)形式變化的雙重介質(zhì)試井模型。Zhang L等[10]考慮壓力變化對孔隙度和滲透率的動態(tài)影響，建立第二類邊界條件下可變形雙重介質(zhì)分形氣藏的非穩(wěn)態(tài)流數(shù)學(xué)模型，基于分形理論，證明模型存在收斂的有限元離散數(shù)值解。姜瑞忠等[11]分析低滲透儲層非線性滲流的微觀機(jī)理，引入屈服應(yīng)力和邊界層流控制的非線性參數(shù)，以Hagen-Poisseuille定律為基礎(chǔ)，推導(dǎo)得到低速非線性滲流模型?；谫|(zhì)量守恒原理和達(dá)西定律，Rao L等[12]建立雙孔介質(zhì)—三線性流動模型，分析雙重孔隙參數(shù)、裂縫流動能力、二次梯度項和表皮因數(shù)對量綱一的壓力的影響。Razminia K等[13]應(yīng)用分形幾何和分?jǐn)?shù)階微分的概念，研究雙區(qū)復(fù)合儲層壓力的數(shù)學(xué)解析方法。Amin D等[14]利用瞬態(tài)產(chǎn)量遞減分析估算，提出計算天然裂縫性儲層平均壓力的方法，與物質(zhì)平衡法得到的結(jié)果一致。劉建軍等[15]根據(jù)復(fù)雜學(xué)科的多尺度關(guān)聯(lián)方法及多尺度裂縫建模的相關(guān)文獻(xiàn)，提出儲層中多尺度裂縫的分級標(biāo)準(zhǔn),分析連續(xù)介質(zhì)和離散介質(zhì)構(gòu)建的代表性模型。

人們研究雙重介質(zhì)致密氣藏的壓力特點，簡化儲層的滲流模型，但不能充分體現(xiàn)雙重介質(zhì)致密氣藏的特點。筆者考慮該類儲層的裂縫分形特征、致密基質(zhì)塊的非線性滲流機(jī)理及滲透率各向異性，建立雙重介質(zhì)致密分形氣藏水平井非線性滲流模型，利用等價壓力點處理模型內(nèi)邊界條件，基于有限元原理求解井底壓力，繪制井底壓力動態(tài)曲線，對比不同模型滲流規(guī)律并對相關(guān)參數(shù)進(jìn)行敏感性分析。

1 滲流機(jī)理

1.1 裂縫分形特征

達(dá)西滲流理論基于歐式空間，描述儲層裂縫的復(fù)雜性與非均質(zhì)性具有局限性。引入分形理論，利用分形維數(shù)及異常擴(kuò)散系數(shù)表征分形裂縫的滲透率和孔隙度，使模型更精確地描述實際儲層的裂縫形態(tài)[16]。

令分形體內(nèi)流體儲集在某節(jié)點處且每處節(jié)點體積相同，由節(jié)點的密度得出滲流節(jié)點的數(shù)量α為

(1)

式中：Df為裂縫分形維數(shù)；ds為譜維數(shù)；μ為流體黏度；rw為井筒半徑；kfw為生產(chǎn)井處的滲透率；φfw為生產(chǎn)井處的孔隙度；Cft為裂縫綜合壓縮系數(shù)。

用d表示分形體嵌入巖塊的歐式維數(shù)，則裂縫孔隙度φf為

(2)

式中：B為描述對稱性的幾何常量；Vs為儲集流體的節(jié)點體積；Vf為裂縫孔隙體積；V為儲層總體積；r為距生產(chǎn)井的距離。

裂縫孔隙度主要受孔隙空間的聚集方式(可用分形維數(shù)表征)影響，與有效上覆壓力呈指數(shù)關(guān)系。引入裂縫壓縮系數(shù)Cf，考慮應(yīng)力敏感性的分形孔隙度φf(r)為

(3)

式中：mi為原始地層擬壓力；mf為裂縫系統(tǒng)擬壓力。

裂縫之間的連通性(可用異常擴(kuò)散系數(shù)表征)影響裂縫滲透率。Pedrosa O A研究應(yīng)力敏感地層不穩(wěn)定試井，證明滲透率與基巖上覆壓力呈指數(shù)關(guān)系[17]。引入異常擴(kuò)散系數(shù)θ，考慮應(yīng)力敏感性的分形滲透率kf(r)為

(4)

式中：γm為滲透率模數(shù)。

1.2 基質(zhì)非線性滲流

低滲透儲層研究常用模型為擬啟動壓力梯度模型、分段模型及連續(xù)模型。3種模型基于滲流實驗，由不同函數(shù)擬合實驗數(shù)據(jù)所得，由于物理背景不豐富，很難解釋低滲透儲層啟動壓力梯度和非線性滲流產(chǎn)生的根本原因。筆者在低滲透儲層微觀滲流機(jī)理及毛管滲流模型的基礎(chǔ)上，改進(jìn)并簡化低速非線性滲流模型。

考慮啟動壓力梯度及非線性滲流，在Hagen-Poisseuille定律的基礎(chǔ)上進(jìn)行修正可得

(5)

式中：ν為通過巖心的流速；K為巖心滲透率；r0為毛細(xì)管半徑；δ為表邊界層厚度；τ0為表流體屈服應(yīng)力；p為壓力梯度。

同一根毛管，壓力梯度越大，邊界層越薄，即δ/r0與壓力梯度成反比關(guān)系，令δ/r0=a1/p；同一種流體屈服應(yīng)力不變，令8τ0/3r0=a2(a1、a2為實驗數(shù)據(jù)擬合所得常數(shù))。式(5)變形并簡化為

(6)

(7)

式中：c1、c2為實驗數(shù)據(jù)擬合所得非線性參數(shù),c1反映流體的屈服應(yīng)力及邊界層對滲流的影響，c2主要反映邊界層對滲流的影響。

當(dāng)c1=0時，式(7)為達(dá)西模型；當(dāng)c2=0時，式(7)為擬啟動壓力梯度模型，c1/(p-c2)為啟動壓力梯度項。

令ν≥0，滿足p≥c1+c2，則最小啟動壓力梯度為pmin=c1+c2。當(dāng)c1+c2=0時，pmin=0，即滲流不存在啟動壓力梯度；當(dāng)c1+c2>0時，pmin>0，滲流存在啟動壓力梯度，式(7)變形為

(8)

此時，滲流曲線不通過原點，隨壓力梯度增大曲線趨于直線，即壓力梯度增大到一定程度，所有孔喉參與流動，邊界層厚度不再發(fā)生變化開始出現(xiàn)擬線性流。擬線性流段的反向延長線與橫軸的交點即啟動壓力梯度。

2 滲流模型

2.1 物理模型

雙重介質(zhì)致密分形氣藏水平井物理模型見圖1。假設(shè)條件：

(1)模型由基質(zhì)系統(tǒng)和裂縫系統(tǒng)組成，是基質(zhì)系統(tǒng)向裂縫系統(tǒng)發(fā)生擬穩(wěn)態(tài)竄流的雙孔單滲的圓形氣藏；

(2)外邊界封閉或定壓，儲層厚度為h，地層擬初始壓力為mi，水平井半長為L，氣藏半徑為re；

(3)水平井位于氣藏平面中心，地面產(chǎn)氣量為qsc；

(4)考慮儲層的表皮因數(shù)、井筒儲集系數(shù)及滲透率的各向異性，裂縫水平滲透率為Kfh，裂縫垂直滲透率為Kfv，基質(zhì)系統(tǒng)滲透率為Km；

(5)考慮裂縫系統(tǒng)的應(yīng)力敏感性和分形特征，以及基質(zhì)向裂縫系統(tǒng)的非線性滲流；

(6)忽略重力與毛管壓力的影響。

圖1 雙重介質(zhì)致密分形氣藏水平井物理模型

2.2 數(shù)學(xué)模型

2.2.1 模型建立

在圓柱坐標(biāo)系下，運動方程式(7)結(jié)合質(zhì)量守恒方程及狀態(tài)方程，考慮儲層應(yīng)力敏感性，得到雙重介質(zhì)致密分形氣藏水平井非線性滲流模型[18-19]：

(9)

(10)

限制條件：當(dāng)mfD-mmD<(c1D+c2D)/2時，δLVD=0；當(dāng)mfD-mmD≥(c1D+c2D)/2時，δLVD=1-c1D/[2(mfD-mmD)-c2D]。

初始條件為

mfD|tD=0=mmD|tD=0=0。

(11)

內(nèi)邊界條件為

(12)

式中：zwD、εD分別為量綱一的水平井縱向坐標(biāo)和微變量。

封閉外邊界條件為

(13)

(14)

式中：reD為量綱一的泄流半徑。

定壓外邊界條件為

mmD(rD=reD)=mmD(zD=0,1)=0;

(15)

mfD(rD=reD)=mfD(zD=0,1)=0。

(16)

2.2.2 模型求解

對雙重介質(zhì)致密分形氣藏水平井模型，采用隱式方法求解裂縫系統(tǒng)、顯式方法求解基質(zhì)系統(tǒng)。

應(yīng)用Galerkin法得到裂縫系統(tǒng)的有限元方程(不存在源匯項時)為

(17)

式中：Ni為形函數(shù)，i=1,2,…n。

以格林公式為基礎(chǔ)，通過分部積分得到內(nèi)部單元和封閉條件外邊界單元的有限元方程為

(18)

裂縫系統(tǒng)單元有限元方程的矩陣形式為

(19)

對矩陣方程化簡得

(20)

式中：Ke為非線性系數(shù)矩陣；Fe為單元載荷向量；n為時刻變量。

式(20)中的Ke、Fe表達(dá)式為

(21)

(22)

同理，基質(zhì)的有限元矩陣方程化簡得

(23)

式(23)中的Ke、Fe表達(dá)式為

(24)

(25)

式(20)和式(23)分別為裂縫系統(tǒng)和基質(zhì)系統(tǒng)的有限元單元平衡方程。根據(jù)基質(zhì)場的上一時刻值求解裂縫場的壓力解，再將該時刻裂縫場壓力解引入基質(zhì)場求基質(zhì)壓力解。

模型求解過程中需要處理水平井的內(nèi)邊界條件，無限導(dǎo)流能力模型符合水平井的實際情況，考慮無限導(dǎo)流能力模型比均勻流量模型求解困難，選取0.7L處為水平井在兩種模型下的等價壓力點，利用均勻流量模型評價水平井的井底壓力[20]。

將水平井考慮為單元內(nèi)源匯項并積分，應(yīng)用Delta函數(shù)將水平井劃分n個節(jié)點，得到量剛一的單元源匯項有限元方程為

(26)

考慮井筒儲集效應(yīng)與表皮因數(shù)，利用Duhamel原理并對壓力解進(jìn)行Laplace變換得

(27)

式中：S為表皮因數(shù)；s為Laplace變量；CD為井筒儲集系數(shù)。

采用Stehfest數(shù)值反演方法，得到考慮井筒儲集效應(yīng)與表皮因數(shù)的水平井井底擬壓力解為

(28)

(29)

3 數(shù)值模擬

3.1 模型對比

雙重介質(zhì)致密氣藏水平井的達(dá)西模型、非線性模型及分形非線性模型的擬壓力及壓力導(dǎo)數(shù)曲線見圖2。其中：CD=100，S=0，ω=0.02，λ=10-5；對于達(dá)西模型，γmD=0.02，cD=0，θ=0，Df=2.0；對于非線性模型，γmD=0.02，cD=0.5，θ=0，Df=2.0；對于分形非線性模型，γmD=0.02，cD=0.3，θ=0.2，Df=1.9。設(shè)非線性參數(shù)cD=c1D=c2D，分形指數(shù)β=2+θ-Df?？紤]應(yīng)力敏感性越強(qiáng)，曲線末端上翹幅度越大[21-22]。

根據(jù)曲線形態(tài)，將雙重介質(zhì)致密氣藏水平井滲流劃分7個流動階段，不同模型在各流動階段的試井曲線呈不同形態(tài)，流動階段劃分及模型對比見表1。

表1 流動階段劃分及模型對比

3.2 敏感性分析

對非線性參數(shù)cD、水平井長度LD及分形指數(shù)β進(jìn)行敏感性分析。

非線性參數(shù)對壓力動態(tài)曲線的影響見圖3。由圖3可見，非線性參數(shù)主要影響基質(zhì)向裂縫系統(tǒng)的竄流強(qiáng)度及出現(xiàn)時間。隨非線性參數(shù)的增大，竄流的非線性增強(qiáng)，竄流強(qiáng)度減弱且竄流段出現(xiàn)時間延遲，在擬壓力導(dǎo)數(shù)曲線上表現(xiàn)為下凹段右移且下凹程度減小，非線性參數(shù)不影響基質(zhì)向裂縫竄流的結(jié)束時間。

水平井長度對壓力動態(tài)曲線的影響見圖4。由圖4可見，隨水平井長度的增加，早期徑向流階段越明顯，早期線性流階段曲線的斜率越大，且不同長度水平井的擬壓力導(dǎo)數(shù)曲線在中期徑向流起點處相交。水平井增長延遲竄流出現(xiàn)的時間，不影響竄流持續(xù)的時間及強(qiáng)度。考慮儲層應(yīng)力敏感性與非線性滲流，不同長度水平井的擬壓力導(dǎo)數(shù)曲線在末端不是匯聚成一條直線，而是上翹形成一簇平行線。

分形指數(shù)對壓力動態(tài)曲線的影響見圖5。由圖5可見，初始階段分形指數(shù)對曲線形態(tài)影響不明顯，隨時間的增加，各導(dǎo)數(shù)曲線在中期徑向流的起點處相交。分形指數(shù)越大，中期徑向流階段后擬壓力導(dǎo)數(shù)曲線上翹幅度越大。分形特征使壓力波在儲層中逐漸傳播并形成動態(tài)邊界，導(dǎo)致系統(tǒng)在徑向流階段后擬壓力導(dǎo)數(shù)曲線上翹幅度逐漸增大，體現(xiàn)分形氣藏的典型特征[23]。

圖2 不同模型壓力動態(tài)曲線對比Fig.2 Comparison of different models pressure dynamic curves

圖3 非線性參數(shù)對壓力動態(tài)曲線的影響Fig.3 Influence of nonlinear parameter on pressure dynamic curves

圖4 水平井長度對壓力動態(tài)曲線的影響Fig.4 Influence of horizontal well length on pressure dynamic curves

圖5 分形指數(shù)對壓力動態(tài)曲線的影響Fig.5 Influence of fractal index on pressure dynamic curves

4 結(jié)論

(1)引入分形理論描述裂縫系統(tǒng)的應(yīng)力敏感性及分形特征，應(yīng)用非線性滲流模型描述致密基質(zhì)塊的滲流特性，考慮啟動壓力梯度，建立雙重介質(zhì)致密分形氣藏水平井非線性滲流模型。

(2)采用等價壓力點方法處理內(nèi)邊界條件，利用有限元原理求解模型；繪制并對比3種模型的井底壓力動態(tài)曲線，劃分7個滲流階段，即井筒儲存、表皮效應(yīng)過渡、早期徑向流、早期線性流、中期徑向流、基質(zhì)向裂縫竄流及基質(zhì)向裂縫擬徑向流。

(3)非線性參數(shù)增大，擬壓力導(dǎo)數(shù)曲線下凹段右移且下凹程度減??；水平井長度增大，竄流出現(xiàn)的時間延遲；分形指數(shù)越大，中期徑向流段后擬壓力導(dǎo)數(shù)曲線上翹幅度越大，各擬壓力導(dǎo)數(shù)曲線在中期徑向流的起點處相交。