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

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

0 引言

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

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

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

1 滲流機(jī)理

1.1 裂縫分形特征

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

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

(1)

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

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

(2)

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

裂縫孔隙度主要受孔隙空間的聚集方式(可用分形維數(shù)表征)影響，與有效上覆壓力呈指數(shù)關(guān)系。引入裂縫壓縮系數(shù)Cf，考慮應(yīng)力敏感性的分形孔隙度φf(shuō)(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ì)非線性滲流

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

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

(5)

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

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

(6)

(7)

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

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

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

(8)

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

2 滲流模型

2.1 物理模型

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

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

(2)外邊界封閉或定壓，儲(chǔ)層厚度為h，地層擬初始?jí)毫閙i，水平井半長(zhǎng)為L(zhǎng)，氣藏半徑為re；

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

(4)考慮儲(chǔ)層的表皮因數(shù)、井筒儲(chǔ)集系數(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)系下，運(yùn)動(dòng)方程式(7)結(jié)合質(zhì)量守恒方程及狀態(tài)方程，考慮儲(chǔ)層應(yīng)力敏感性，得到雙重介質(zhì)致密分形氣藏水平井非線性滲流模型[18-19]：

(9)

(10)

限制條件：當(dāng)mfD-mmD<(c1D+c2D)/2時(shí)，δLVD=0；當(dāng)mfD-mmD≥(c1D+c2D)/2時(shí)，δ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 模型求解

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

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

(17)

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

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

(18)

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

(19)

對(duì)矩陣方程化簡(jiǎn)得

(20)

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

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

(21)

(22)

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

(23)

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

(24)

(25)

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

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

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

(26)

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

(27)

式中：S為表皮因數(shù)；s為L(zhǎng)aplace變量；CD為井筒儲(chǔ)集系數(shù)。

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

(28)

(29)

3 數(shù)值模擬

3.1 模型對(duì)比

雙重介質(zhì)致密氣藏水平井的達(dá)西模型、非線性模型及分形非線性模型的擬壓力及壓力導(dǎo)數(shù)曲線見(jiàn)圖2。其中：CD=100，S=0，ω=0.02，λ=10-5；對(duì)于達(dá)西模型，γmD=0.02，cD=0，θ=0，Df=2.0；對(duì)于非線性模型，γmD=0.02，cD=0.5，θ=0，Df=2.0；對(duì)于分形非線性模型，γ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個(gè)流動(dòng)階段，不同模型在各流動(dòng)階段的試井曲線呈不同形態(tài)，流動(dòng)階段劃分及模型對(duì)比見(jiàn)表1。

表1 流動(dòng)階段劃分及模型對(duì)比

3.2 敏感性分析

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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