裂縫-孔隙介質(zhì)儲(chǔ)層滲透率表征及其影響因素分析

熊繁升， 甘利燈， 孫衛(wèi)濤,2*， 楊昊

1 清華大學(xué)周培源應(yīng)用數(shù)學(xué)研究中心, 北京 100084 2 清華大學(xué)航天航空學(xué)院, 北京 100084 3 中國(guó)石油集團(tuán)勘探開發(fā)研究院， 北京 100083

0 引言

巖心切片電子顯微鏡實(shí)驗(yàn)和數(shù)字巖心技術(shù)等最新手段研究證實(shí)，在儲(chǔ)層巖石內(nèi)部普遍存在裂縫-孔隙形成的復(fù)雜連通網(wǎng)絡(luò).在天然巖石中，流體流動(dòng)發(fā)生的空間是全局性裂縫-孔隙網(wǎng)絡(luò)，而不是局部的單一孔隙或孔管束，因此三維裂縫-孔隙網(wǎng)絡(luò)模型是描述孔隙空間非均質(zhì)性的一種合理方法.Fatt(1955)較早提出了孔隙介質(zhì)簡(jiǎn)化網(wǎng)絡(luò)模型，通過規(guī)則網(wǎng)絡(luò)上的管道表示孔隙空間.Bryant等(1993)提出了一個(gè)多孔介質(zhì)的無序網(wǎng)絡(luò)模型，認(rèn)為真實(shí)多孔系統(tǒng)孔隙尺寸隨機(jī)分布假設(shè)可能是無效的.Seeburger和Nur(1984)提出了二維孔隙空間網(wǎng)絡(luò)模型，用于預(yù)測(cè)砂巖的靜滲透率，砂巖和花崗巖的計(jì)算結(jié)果表明，網(wǎng)絡(luò)模型可以根據(jù)圍壓函數(shù)預(yù)測(cè)滲透率大小.基于流體力學(xué)基本方程.Bernabé(2009)建立了牛頓流體在剛性和彈性管網(wǎng)中的水力傳導(dǎo)率和波頻散方程.Bernabé(1995)通過管道和管網(wǎng)對(duì)牛頓流體的輸運(yùn)特性進(jìn)行了數(shù)值模擬.Johnson等(1987)得到了含流體管道中滲透率對(duì)頻率的依賴性.Blunt和King(1991)、Blunt(2001)研究了由細(xì)管形成的網(wǎng)絡(luò)中相對(duì)滲透率與飽和度的關(guān)系，表明網(wǎng)絡(luò)模型能夠解釋微觀流動(dòng)如何影響厘米尺度上的平均特性.基于孔隙尺度的潤(rùn)濕性變化，三維孔隙網(wǎng)絡(luò)也被用于模擬砂巖內(nèi)部排水和注水過程(Nilsen et al., 1996; Bakke and ?ren, 1997).

大量研究表明，雖然局部非均勻孔隙/裂縫參數(shù)對(duì)孔隙流體會(huì)產(chǎn)生一定影響，但是在孔隙網(wǎng)絡(luò)節(jié)點(diǎn)相互連通的“均質(zhì)化”效應(yīng)下，其對(duì)儲(chǔ)層巖石整體滲透率的影響十分有限.平行排列的裂縫系統(tǒng)能夠較好的反映地層各向異性屬性，但是很難描述跨平行裂縫層的流體流動(dòng)現(xiàn)象.因此，相對(duì)于研究局部孔隙-裂縫相互作用的模型來說，更需要將跨微觀-宏觀尺度的三維裂縫/孔隙網(wǎng)絡(luò)系統(tǒng)作為研究對(duì)象，從整體角度考察裂縫網(wǎng)絡(luò)對(duì)巖石孔隙度、滲透率、地震波頻散衰減、提高儲(chǔ)層采收率等諸多方面的影響.

滲透率是多孔介質(zhì)的一種固有特性，它表征了多孔介質(zhì)中流體流動(dòng)的難易程度.達(dá)西定律描述了流體在多孔介質(zhì)中的穩(wěn)態(tài)流動(dòng)，在這種情況下滲透率是一個(gè)靜態(tài)值.然而，天然巖石孔隙結(jié)構(gòu)十分復(fù)雜，在波場(chǎng)作用下內(nèi)部流體-固體交界面振動(dòng)引起孔隙連通性改變.彎曲度是描述多孔介質(zhì)中孔隙連通性和擴(kuò)散路徑曲率的常用量，隨周圍壓力的變化而變化(Johnson et al., 1987; Smeulders et al., 2006).因此，振蕩流動(dòng)過程中存在“動(dòng)態(tài)滲透率”，換句話說，如果流體流動(dòng)是由振蕩壓力引起的，滲透率將隨頻率而變化.

實(shí)驗(yàn)測(cè)量發(fā)現(xiàn)，滲透率在低頻率下是一個(gè)恒定值，在高頻率下依賴于頻率的倒數(shù)(頻率范圍從0.1 Hz到1 kHz)(Charlaix et al., 1988)，一般認(rèn)為隨著頻率增大，存在著從黏性流態(tài)到慣性流態(tài)的過渡，實(shí)驗(yàn)觀測(cè)也提供了黏彈性流體滲透率增強(qiáng)的證據(jù)(Mena et al., 1979; Castrejón-Pita et al., 2003).Johnson等(1987)進(jìn)行了動(dòng)態(tài)滲透率模型的研究，采用分支函數(shù)法推導(dǎo)了隨機(jī)孔隙中與波的耗散相關(guān)的動(dòng)態(tài)滲透率一般表達(dá)式，并通過設(shè)計(jì)不同實(shí)驗(yàn)證實(shí)這個(gè)模型可以很好描述高低頻極限下彈性波的響應(yīng)規(guī)律，在中間頻率能給出較合理的近似.Blanc等(2013)從流體固體的本構(gòu)關(guān)系和動(dòng)量守恒方程式出發(fā)，將動(dòng)態(tài)滲透率結(jié)合到二維孔隙介質(zhì)波動(dòng)方程中，修正了耗散項(xiàng)，通過算子轉(zhuǎn)化實(shí)現(xiàn)在時(shí)間域進(jìn)行求解.Müller等(2007)研究了隨機(jī)分布孔隙空間中的動(dòng)態(tài)滲透率變化規(guī)律，并推導(dǎo)此類分布下動(dòng)態(tài)滲透率的數(shù)學(xué)表達(dá)式.Chapman等(2002)提出了從微觀尺度建立孔隙彈性介質(zhì)模型的方法，該模型假設(shè)孔隙網(wǎng)絡(luò)節(jié)點(diǎn)與臨近其他節(jié)點(diǎn)之間的流體流動(dòng)滿足達(dá)西定律，利用達(dá)西定律導(dǎo)出宏觀流動(dòng)-壓力方程，但是該模型并沒有給出宏觀裂縫系統(tǒng)滲透率表達(dá)公式，滲透率作為一個(gè)計(jì)算參數(shù)，需要通過實(shí)驗(yàn)手段測(cè)量得到.Xiong和Sun(2020) 分析了任意三維孔隙網(wǎng)絡(luò)中流體流變學(xué)特征和孔隙連通關(guān)系對(duì)巖石滲透率的影響，發(fā)現(xiàn)滲透率變化受脈沖壓力波場(chǎng)頻率、流體流變學(xué)屬性和孔隙連通性的影響.

在含有裂縫/軟孔隙的儲(chǔ)層巖石中，流體流動(dòng)路徑往往受到孔隙截面閉合的影響而發(fā)生改變.但是，目前關(guān)于裂縫-孔隙網(wǎng)絡(luò)的研究多集中于圓形截面的管道網(wǎng)絡(luò)，這種網(wǎng)絡(luò)空間可以用來表征硬孔隙對(duì)滲透率的影響，但無法解釋在壓力作用下軟孔隙壓縮甚至閉合帶來的影響.利用橢圓截面微管網(wǎng)絡(luò)模擬儲(chǔ)層巖石流體流動(dòng)空間的工作尚不多見，主要原因在于圓形截面的孔隙管道流體流動(dòng)問題相對(duì)簡(jiǎn)單，易于給出簡(jiǎn)潔的表達(dá)式；而橢圓截面孔隙管道的流場(chǎng)速度數(shù)學(xué)建模復(fù)雜，很難用簡(jiǎn)單函數(shù)得到解析表達(dá)式.

基于巖石內(nèi)部的裂縫連通網(wǎng)絡(luò)特征，本文提出了具有橢圓截面的裂縫/軟孔隙三維網(wǎng)絡(luò)模型(圖1).該模型通過在橢圓坐標(biāo)系下為單個(gè)橢圓截面裂縫/軟孔隙管道流場(chǎng)速度建模，利用無量綱流量與縱橫比無關(guān)的條件，將橢圓坐標(biāo)系下馬蒂厄函數(shù)轉(zhuǎn)換為圓柱坐標(biāo)系下貝塞爾函數(shù)，獲得裂縫/軟孔隙流量解析表達(dá)式，并利裂縫/軟孔隙網(wǎng)絡(luò)節(jié)點(diǎn)流通量守恒條件，得到振蕩壓力作用下儲(chǔ)層巖石的整體滲透率計(jì)算方法.

1 三維橢圓截面微管網(wǎng)絡(luò)模型

在三維橢圓截面微管網(wǎng)絡(luò)模型中，引入了裂縫體密度、裂縫縱橫比和裂縫長(zhǎng)度、孔隙度等參數(shù)，并在宏觀滲透率與微觀參數(shù)之間建立了聯(lián)系，研究發(fā)現(xiàn)這些參數(shù)對(duì)于孔隙介質(zhì)整體滲透率具有重要影響.

給定邊長(zhǎng)為l的巖石樣本單元立方體，沿X、Y、Z方向上分別排布有M、N、L個(gè)橢圓截面微管(圖1).在X-Y平面內(nèi)，分布有M×N個(gè)裂縫，并沿Z方向排布有(L-1)個(gè)間隔，因此沿Z方向上共有M×N×(L-1)個(gè)微裂縫.同理，在Y方向上分布有M×L×(N-1)個(gè)微裂縫，在X方向上分布有N×L×(M-1)個(gè)微裂縫.單元體內(nèi)的裂縫數(shù)密度定義為：

圖1 由橢圓截面毛細(xì)管構(gòu)成的三維裂縫網(wǎng)絡(luò)模型Fig.1 Schematic diagram of 3D fracture network model composed of capillary tubes with elliptic cross section

(1)

(2)

繼而得到孔隙度為：

(3)

在后續(xù)計(jì)算中，常常需要給定裂縫網(wǎng)絡(luò)體密度、孔隙度，以便考察滲透率隨裂縫縱橫比的變化情況.因此，這里給出已知裂縫體密度參數(shù)M、N、L和孔隙度φ，求解裂縫主半徑長(zhǎng)度的公式：

(4)

其中，巖石樣本邊長(zhǎng)l可以取做實(shí)驗(yàn)特征長(zhǎng)度，或者取做滲透率預(yù)測(cè)所需要尺度.在給定孔隙度情況下，當(dāng)裂縫的縱橫比很小，裂縫的寬度可能非常大，甚至超過了巖石樣本的邊長(zhǎng)l，此時(shí)把裂縫的寬度限制在一定范圍內(nèi)，例如在X方向上，裂縫最大寬度小于l/(M-1).這種情況下，裂縫網(wǎng)絡(luò)的孔隙度可能會(huì)小于預(yù)設(shè)值.因此，在縱橫比取非常小的數(shù)值時(shí)，單個(gè)裂縫/軟孔隙的空間受到限制，模型孔隙度存在一個(gè)上限.該模型具有以下特點(diǎn)：

(1)單個(gè)管道橫截面是縱橫比可變的橢圓形，可以模擬從狹窄裂縫/軟孔隙到圓形硬孔隙的不同類型裂縫-孔隙網(wǎng)絡(luò)空間；

(2)在微觀尺度上單個(gè)微管(半徑處于微米量級(jí))相互連接，形成宏觀網(wǎng)絡(luò)(整體單元處于厘米～米量級(jí))，實(shí)現(xiàn)孔隙介質(zhì)中跨尺度整體建模，能夠更好的反映孔隙系統(tǒng)整體效應(yīng)；

(3)將裂縫/軟孔隙近似成橢圓微管網(wǎng)絡(luò)，可以簡(jiǎn)化物理建模和數(shù)學(xué)推導(dǎo)，得到精度很高的解析解，同時(shí)又不失對(duì)裂縫/軟孔隙整體空間的特征刻畫；

(4)裂縫/軟孔隙網(wǎng)絡(luò)模型的滲透率依賴于壓力振蕩頻率，因此可以應(yīng)用于不同頻率下滲透率變化情況.

當(dāng)然，裂縫網(wǎng)絡(luò)是一種等效孔隙介質(zhì)模型，其中包含的裂縫寬度、長(zhǎng)度等參數(shù)與巖石切片中觀察到的實(shí)際尺寸并非一一對(duì)應(yīng)，而是對(duì)真實(shí)儲(chǔ)層巖石中裂縫/孔隙理想化建模的等效參數(shù).等效參數(shù)與實(shí)驗(yàn)觀測(cè)數(shù)值之間存在不一致，甚至是相差很大，但是對(duì)于等效巖石物理模型來說，這些參數(shù)數(shù)值是合適的，能夠從整體上實(shí)現(xiàn)對(duì)儲(chǔ)層巖石滲透率的預(yù)測(cè).通過對(duì)來自不同文獻(xiàn)的多種類型砂巖實(shí)驗(yàn)數(shù)據(jù)對(duì)比表明，三維裂縫/軟孔隙網(wǎng)絡(luò)模型確實(shí)可以定量預(yù)測(cè)出多種巖樣滲透率分布范圍，并估計(jì)出對(duì)應(yīng)巖石孔隙縱橫比.

2 橢圓截面微管中流體速度場(chǎng)和流量

首先采用橢圓柱坐標(biāo)系，給出含不可壓縮牛頓流體在橢圓截面微管中的流量表達(dá)式(圖2).

圖2 橢圓柱坐標(biāo)系下橢圓截面微管模型Fig.2 Microtubule model with elliptic section in elliptic cylindrical coordinate system

橢圓柱坐標(biāo)徑向、軸向坐標(biāo)分別以ξ、ψ和z表示，根據(jù)微管道對(duì)稱性條件，流體質(zhì)點(diǎn)軸向速度可以表示為u=u(r,z,t)，徑向速度表示為v=v(r,z,t)，流體壓力表示為p=p(r,z,t)，壓力振蕩頻率為ω，微管中流體的密度、粘度分別為ρf、η.

基于長(zhǎng)波長(zhǎng)假設(shè)，壓力梯度沿微管道軸向分布，徑向速度分量為零(Maslen,1958).微管中的流體在壓力作用下發(fā)生流動(dòng)，并使壓力在微管流體中發(fā)生傳播.孔隙流體流速很低，因此流體壓縮性可以忽略.考慮到流體速度非零分量沿z軸(軸向)方向，其大小隨徑向坐標(biāo)而變化，流體滿足方程：

(5)

其中ρf是流體密度，p是流體壓力.設(shè)橢圓截面兩個(gè)焦點(diǎn)之間的距離是2d，ψ=ψ0表示一族共焦點(diǎn)雙曲線，ξ=ξ0表示一族共焦點(diǎn)橢圓，這里ξ0和ψ0是常數(shù).利用坐標(biāo)變換關(guān)系式：

(6)

得到橢圓柱坐標(biāo)系下的Navier-Stokes方程：

(7)

u(ξ,ψ,z,t)=U(ξ,ψ)e-iω(t-z/c)，

(8)

p(ξ,ψ,z,t)=P(ξ,ψ)e-iω(t-z/c)，

(9)

代入Navier-Stokes方程得到：

(10)

式(10)是關(guān)于U的非齊次微分方程, 這里c是流體中聲波的速度(復(fù)數(shù)).

當(dāng)裂縫/孔隙壓力隨時(shí)間和空間變化時(shí)，流體的速度場(chǎng)是時(shí)間和空間的函數(shù)，可以表示為穩(wěn)態(tài)場(chǎng)和非穩(wěn)態(tài)場(chǎng)的疊加：

u(x,y,t)=u0(x,y,z)+u1(x,y,z,t)，

(11)

代入Navier-Stokes方程并利用邊界條件和分離變量方法，得到穩(wěn)態(tài)速度場(chǎng)的空間分布(Haslam, 1997)：

u0(ξ,ψ)=

(12)

通過在橢圓微管截面上積分，容易得到穩(wěn)態(tài)流量表達(dá)式：

(13)

將速度寫為壓力P(ξ,ψ)的函數(shù)，代入Navier-Stokes方程得到非穩(wěn)態(tài)解(Haslam, 1997)：

(14)

(15)

(16)

(17)

這里Qc(t)是圓形截面微管的流量.因此可以通過圓形截面微管流量和縱橫比來計(jì)算橢圓形截面微管流量.

3 裂縫/軟孔隙流體流量與壓力振動(dòng)頻率關(guān)系

Bernabé(Bernabé, 2009)研究了諧波振蕩壓力下圓柱形管道的流體流動(dòng)問題，給出了流體流量與頻率、流體密度、粘性、管道半徑和長(zhǎng)度等參數(shù)之間的關(guān)系式.本文將該方法推廣到了具有橢圓形截面的微管流量計(jì)算中.下面考慮橢圓裂縫微管長(zhǎng)度為L(zhǎng)0，管兩端壓力分別為PU、PD，流量分別為QU、QD，此時(shí)管內(nèi)壓力為沿相反方向傳播的兩種壓力波的疊加：

p(z,t)=B+e-iω(t-z/c)+B-e-iω(t+z/c)，

(18)

管內(nèi)總流量表示為

q(z,t)=Q+e-iω(t-z/c)+Q-e-iω(t+z/c).

(19)

利用邊界條件p(0,0)=PU和p(L0,0)=PD，得到：

(20)

利用初始邊界條件和無量綱流量關(guān)系式，得到壓力波作用下橢圓截面微管流量表達(dá)式：

(21)

其中：

(22)

結(jié)合達(dá)西定律得到橢圓柱體通道的流體傳導(dǎo)系數(shù)：

(23)

其中A=aπR2為截面面積，R是橢圓截面長(zhǎng)軸半徑.因此，流量與壓力之間滿足：

Qe=k(PU-PD)，

(24)

4 裂縫/軟孔隙流體流量與圍壓的關(guān)系

對(duì)于橢圓截面的裂縫，當(dāng)裂縫截面縱橫比a較大時(shí)，截面接近正圓，孔隙結(jié)構(gòu)抵抗外部圍壓的能力較強(qiáng)，不易變形；當(dāng)裂縫截面縱橫比a很小時(shí)，裂縫截面兩側(cè)邊緣互相接近，抵抗外部圍壓的能力較弱，容易發(fā)生較大變形，影響流體的流動(dòng)特征.因此，當(dāng)a較小時(shí)，在孔隙/裂縫滲透率計(jì)算中需要引入孔隙/裂縫的彈性形變理論.

考慮截面積相對(duì)較小的裂縫，其長(zhǎng)度方向延伸較長(zhǎng)，因此可以將裂縫的彈性形變看作平面應(yīng)變問題.在含裂縫巖石中，當(dāng)裂縫十分狹窄時(shí)，裂縫兩個(gè)邊緣十分接近，裂縫尖端呈楔形.設(shè)裂縫截面的長(zhǎng)軸和短軸半徑分別為R1、R2，Mavko和Nur(1978)研究了楔形尖端裂縫在圍壓作用下的變形問題，并給出了裂縫長(zhǎng)度和寬度方向半徑隨圍壓的變化關(guān)系.當(dāng)孔隙介質(zhì)外部受到外部圍壓pc作用時(shí)，半徑形變與壓力pc之間滿足：

(25)

(26)

其中ν、E是固體骨架礦物材料的泊松比和楊氏模量.因此，裂縫截面縱橫比可以表示為圍壓的函數(shù)：

(27)

將半徑-壓力關(guān)系式代入速度場(chǎng)流量表達(dá)式，得到裂縫流量隨圍壓的變化關(guān)系式為

(28)

其中:

(29)

5 含三維裂縫/軟孔隙網(wǎng)絡(luò)多孔介質(zhì)滲透率計(jì)算

上面得到了單裂縫流量隨縱橫比、壓力變化頻率、巖石圍壓等參數(shù)變化關(guān)系式，為了建立由微裂縫組成的三維裂縫/軟孔隙網(wǎng)絡(luò)孔隙介質(zhì)滲透率，首先采用流量計(jì)算公式得到每個(gè)微管的流道，然后利用微管兩端每個(gè)節(jié)點(diǎn)的流量守恒條件(克?；舴蚨?，得到全部網(wǎng)絡(luò)節(jié)點(diǎn)滿足的流量-壓力線性方程組.該方程組未知量為各節(jié)點(diǎn)的壓力，入口和出口端的壓力邊界條件是方程非齊次項(xiàng).求解該方程組得到各節(jié)點(diǎn)壓力后，計(jì)算巖石樣本整體三維網(wǎng)絡(luò)流量，進(jìn)而根據(jù)壓力邊界條件和流量得到巖石樣本滲透率.網(wǎng)絡(luò)節(jié)點(diǎn)上流量守恒滿足：

(30)

這里Z為與節(jié)點(diǎn)i相連接的其他節(jié)點(diǎn)總數(shù)目.記N為網(wǎng)格的節(jié)點(diǎn)個(gè)數(shù)，Γ為節(jié)點(diǎn)連接矩陣，?！蔙N×N，Γ的元素僅由0和1組成，Γ矩陣元素的賦值規(guī)定為：若節(jié)點(diǎn)i和節(jié)點(diǎn)j相連γij=1，否則γij=0；對(duì)任意節(jié)點(diǎn)i均有γii=0，即節(jié)點(diǎn)與其本身不相連；矩陣Γ是對(duì)稱的.

考慮連接網(wǎng)絡(luò)節(jié)點(diǎn)i和j的任意一段微管，節(jié)點(diǎn)處的壓力分別為pi和pj，微管流量是pi和pj的函數(shù)，采用線性關(guān)系可以表示為

Qij=cijpi+dijpj，

(31)

這里系數(shù)cij和dij是與所選這兩個(gè)節(jié)點(diǎn)之間微管參數(shù)(管長(zhǎng)、半徑)相關(guān)的參數(shù)，采用式(31)作為壓力差與流量的關(guān)系式，則有cij=ij，dij=δij，其中：

(32)

(33)

(34)

這里L(fēng)ij與Rij分別表示節(jié)點(diǎn)i與j所連接微管的長(zhǎng)度與主半徑，aij是截面縱橫比，這里有cij=cos(ωLij/c0)·dij，c0=1/(βρ)是聲波在流體中的傳播速度，β為流體壓縮系數(shù)，ω是脈沖壓力波場(chǎng)頻率.

對(duì)各節(jié)點(diǎn)列流量平衡方程式，得到以下形式方程組：

(35)

將系數(shù)矩陣分裂為兩個(gè)矩陣相加，A=M+N，其中對(duì)角陣M表示為

(36)

6 數(shù)值算例

6.1 砂巖滲透率預(yù)測(cè)與實(shí)驗(yàn)數(shù)據(jù)對(duì)比

為驗(yàn)證所建立的三維裂縫/軟孔隙多孔介質(zhì)滲透率模型可靠性，選取4類巖樣實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比分析：致密氣砂巖，固結(jié)砂巖，非固結(jié)砂巖，新沉積海灘砂.巖樣滲透率隨孔隙度變化的實(shí)驗(yàn)數(shù)據(jù)來自于PetroWiki(2015).

在本算例中，三維裂縫/軟孔隙模型計(jì)算參數(shù)為：空間三個(gè)方向上網(wǎng)絡(luò)維度為M=5,N=5,L=5，模型巖樣邊長(zhǎng)1 cm, 孔隙度范圍為φ=0～55%, 縱橫比范圍為a=0.001～1.砂巖礦物材料密度為2650 kg·m-3，體積模量為37 GPa，剪切模量為44 GPa，砂巖泊松比采用0.2.孔隙流體為水，其密度為1000 kg·m-3，體積模量為2.15 Gpa，黏性為0.001 Pa·s.

圖3 含水砂巖滲透率-孔隙度交會(huì)圖 采用新沉積的海灘砂和來自油氣田測(cè)井的三套樣品(Petrowiki, 2015)，TG-SS: 致密氣砂巖；C-SS:固結(jié)砂巖；PC-SS: 弱固結(jié)砂巖； ND-S: 新沉積海灘砂.Fig.3 Permeability-porosity cross plot of water-bearing sandstone Using newly deposited beach sand and three sets of samples from oil and gas field logging (Petrowiki, 2015).TG-SS: Tight gas sandstone; C-SS: Consolidated sandstone; PC-SS: Weakly consolidated sandstone; ND-S: Newly deposited beach sand.

這里計(jì)算了在圍壓為10 MPa件下含水砂巖滲透率-孔隙度交會(huì)圖(圖3)，孔隙流體承受的振蕩壓力頻率為20 Hz.圖3給出了不同縱橫比下得到的含三維裂縫/軟孔隙網(wǎng)絡(luò)砂巖滲透率-孔隙度關(guān)系曲線.計(jì)算結(jié)果顯示，理論模型預(yù)測(cè)的滲透率-孔隙度變化范圍與實(shí)驗(yàn)數(shù)據(jù)吻合很好(新沉積海灘砂除外).隨著縱橫比從0.04增加到1，滲透率從10-4mD增加至105mD(1 mD=0.987×10-3μm2)，對(duì)于同一個(gè)縱橫比，滲透率隨著孔隙度增加而增大，但是增加趨勢(shì)逐漸平緩.

在數(shù)據(jù)圖中，致密氣砂巖滲透率范圍較低(10-4～10-2mD).模型預(yù)測(cè)結(jié)果顯示，實(shí)驗(yàn)數(shù)據(jù)點(diǎn)分布在縱橫比分別為0.04和0.08的界限之間(圖3).這個(gè)結(jié)果表明，真實(shí)致密氣砂巖中的裂縫/孔隙對(duì)滲透率的影響，可以用三維裂縫網(wǎng)絡(luò)模型中縱橫比為0.04和0.08的等效裂縫/孔隙網(wǎng)絡(luò)表示.

在準(zhǔn)確預(yù)測(cè)出滲透率范圍的同時(shí)，理論模型預(yù)測(cè)的裂縫/軟孔隙縱橫比能否與真實(shí)巖石中縱橫比吻合，這是人們關(guān)心的重要問題.Smith等(2009)在研究低孔低滲致密氣砂巖的速度-孔隙度關(guān)系時(shí)發(fā)現(xiàn)，微裂縫對(duì)致密氣砂巖的P波速度起到至關(guān)重要的作用，需要引入較大的縱橫比分布范圍才能解釋所有的觀測(cè)數(shù)據(jù)點(diǎn).他們的結(jié)果顯示(Smith et al., 2009)，致密氣砂巖的實(shí)驗(yàn)數(shù)據(jù)點(diǎn)分布在以縱橫比0.05為中心的一個(gè)變化范圍內(nèi)，這與本文模型預(yù)測(cè)的縱橫比分布范圍(0.04～0.08)一致，從多個(gè)文獻(xiàn)中不同數(shù)據(jù)和分析方法的相互印證可以看出，三維裂縫網(wǎng)絡(luò)模型對(duì)致密氣砂巖滲透率和孔隙縱橫比的預(yù)測(cè)是可信的.

固結(jié)砂巖滲透率實(shí)驗(yàn)數(shù)據(jù)分布的范圍較為寬泛，覆蓋了從0.01 mD到100 mD的范圍(圖3).Burns等(1990)研究了利用直接觀測(cè)和速度反演手段獲得固結(jié)砂巖孔隙縱橫比的方法，對(duì)Navajo、Weber和Kayenta砂巖樣本的實(shí)驗(yàn)數(shù)據(jù)分析顯示，砂巖孔隙的縱橫比均大于0.01，其中含有黏土的固結(jié)砂巖孔隙縱橫比主要分布在以0.2為均值的范圍內(nèi).這個(gè)結(jié)果與本文對(duì)固結(jié)砂巖的預(yù)測(cè)范圍(0.08～0.4)十分吻合(圖3)，可以看出理論模型預(yù)測(cè)的縱橫比數(shù)值再一次得到了實(shí)驗(yàn)數(shù)據(jù)的驗(yàn)證.

弱固結(jié)砂巖的滲透率實(shí)驗(yàn)數(shù)據(jù)分布在約100 mD到1000 mD之間(圖3).Alam等 (2014a,b)通過實(shí)驗(yàn)手段測(cè)量了不同溫度和圍壓下經(jīng)過短時(shí)間(24 h)固結(jié)的Shikotsu熔結(jié)凝灰?guī)r和Kimachi砂巖孔隙縱橫比，與地下儲(chǔ)層巖石億萬(wàn)年的長(zhǎng)期固結(jié)作用相比，短時(shí)間固結(jié)的樣本是一種弱固結(jié)砂巖.他們的測(cè)量結(jié)果顯示(Alam et al., 2014a,b)，在1 MPa圍壓作用和353 K溫度下，巖石孔隙縱橫比主要分布在0.2至0.85之間，其中以0.45至0.65占優(yōu)；在15 MPa圍壓作用和295 K溫度下，縱橫比主要分布在0.3至0.8之間，其中以0.4至0.6占用.本文模型對(duì)弱固結(jié)砂巖縱橫比的預(yù)測(cè)范圍是0.4至0.7(圖3)，這與Alam等學(xué)者對(duì)弱固結(jié)砂巖的實(shí)驗(yàn)結(jié)果吻合的很好.

新沉積海灘砂的滲透率實(shí)驗(yàn)數(shù)據(jù)分布在約50～100 D之間(圖3)，由于砂粒完全沒有固結(jié)，顆粒之間的孔隙空間松散，本文滲透率模型預(yù)測(cè)的結(jié)果低于實(shí)驗(yàn)測(cè)量值.

6.2 裂縫/軟孔隙體密度對(duì)巖石滲透率的影響

巖石中裂縫/軟孔隙的連通程度對(duì)于滲透率有重要影響.當(dāng)裂縫體密度增加時(shí)，裂縫之間的間距降低，彼此連通度變得更加緊密；當(dāng)體密度降低時(shí)，裂縫之間的間距增加，裂縫變得更加稀疏，相互之間的連通關(guān)系降低.Lucia(1983)研究發(fā)現(xiàn)，碳酸鹽滲透率與裂縫寬度三次方成正比，與裂縫間距成反比.我們分別計(jì)算了裂縫網(wǎng)絡(luò)規(guī)模從4×4×4、5×5×5到6×6×6不斷增加的情況，其對(duì)應(yīng)裂縫體密度分別為144、300和540.為了對(duì)比分析裂縫密度、縱橫比的影響，我們還對(duì)不同孔隙度、不同縱橫比下的網(wǎng)絡(luò)進(jìn)行了計(jì)算(圖4，圖5).

由圖4可以看出，恒定孔隙度條件下，滲透率隨裂縫/軟孔隙密度增加而增大.我們發(fā)現(xiàn)，當(dāng)裂縫/軟孔隙網(wǎng)絡(luò)規(guī)模從4×4×4加密到5×5×5時(shí)，滲透率曲線約上升17%左右；當(dāng)裂縫/軟孔隙網(wǎng)絡(luò)規(guī)模從5×5×5加密到6×6×6時(shí)，滲透率曲線約上升25%左右.因此可以看出，隨著裂縫/軟孔隙密度的增加，裂縫之間的連通程度提高，巖石滲透率逐步增大；反之，當(dāng)裂縫之間連通程度降低時(shí)，滲透率下降.相比于Lucia(1983)對(duì)碳酸鹽的研究工作，本文的理論建模過程從不同角度出發(fā)得到了一致的結(jié)論；同時(shí)，Davudov和Moghanloo(2018)對(duì)頁(yè)巖滲透率隨孔隙可壓縮性和連通性變化規(guī)律的研究發(fā)現(xiàn)，隨著等效壓力的增加，裂縫/孔隙連接數(shù)量逐漸下降，滲透率隨著逐漸降低.這些工作從多個(gè)方面研究了滲透率-裂縫體密度之間的關(guān)系，通過對(duì)比前人工作發(fā)現(xiàn)，三維裂縫/軟孔隙網(wǎng)絡(luò)模型符合已有研究得到的規(guī)律，并且能夠?qū)崿F(xiàn)滲透率隨裂縫參數(shù)變化規(guī)律的定量預(yù)測(cè)計(jì)算.

從圖4、圖5中可以看出，滲透率隨著縱橫比的增加而顯著增大(跨越近8個(gè)數(shù)量級(jí))；我們發(fā)現(xiàn)，隨著縱橫比接近1，這種增大趨勢(shì)逐漸降低.同時(shí)，在給定縱橫比情況下，滲透率隨著孔隙度增加而增大，這種增加趨勢(shì)并非線性，而是逐漸趨于緩和.

在圖4a中，由于縱橫比很小(0.04)，其裂縫網(wǎng)絡(luò)存在一個(gè)孔隙度上限，因此孔隙度范圍的最大值只有約15%左右；當(dāng)縱橫比上增加到0.08時(shí)，裂縫網(wǎng)絡(luò)對(duì)應(yīng)的孔隙度上限也有所增大，達(dá)到約30%左右(圖4b).當(dāng)裂縫縱橫比繼續(xù)增大時(shí)，孔隙度可以達(dá)到預(yù)設(shè)范圍最大值55%左右(圖5).

7 結(jié)論

基于橢圓截面微管道在振蕩壓力作用下的流體速度場(chǎng)和流量建模，提出了含三維裂縫/軟孔隙網(wǎng)絡(luò)多孔介質(zhì)滲透率計(jì)算方法，得到了含多孔介質(zhì)宏觀滲透率與孔隙度、縱橫比、裂縫/軟孔隙體密度等參數(shù)之間的關(guān)系.數(shù)值計(jì)算結(jié)果表明，相同孔隙度下縱橫比對(duì)滲透率有顯著影響，滲透率的數(shù)值可以跨越幾個(gè)數(shù)量級(jí)；相同縱橫比條件下，滲透率隨著孔隙度增加而增大，同時(shí)，這種增長(zhǎng)趨勢(shì)依賴于孔隙度數(shù)值，低孔隙度下增長(zhǎng)趨勢(shì)更明顯.相同孔隙度下，滲透率隨著網(wǎng)絡(luò)密度的增加而增大，而且在更高的裂縫/軟孔隙體密度下，滲透率的增強(qiáng)效應(yīng)更明顯.理論模型數(shù)值計(jì)算與實(shí)驗(yàn)數(shù)據(jù)對(duì)比分析表明，三維裂縫/軟孔隙網(wǎng)絡(luò)模型得到的認(rèn)識(shí)與現(xiàn)有實(shí)驗(yàn)結(jié)果吻合很好，更重要的是該模型提供了一種根據(jù)裂縫幾何特征、壓力、頻率、巖性等實(shí)驗(yàn)可測(cè)量參數(shù)定量計(jì)算儲(chǔ)層巖石宏觀滲透率的方法，計(jì)算結(jié)果考慮到了裂縫/軟孔隙的大范圍連通效應(yīng)，相比局部孔隙、管道建模方法更具真實(shí)性.

圖4 滲透率隨孔隙度、裂縫/軟孔隙密度曲線 (a) 縱橫比為0.04； (b) 縱橫比為0.08.其中δκ1為裂縫網(wǎng)絡(luò)規(guī)模從4×4×4加密到5×5×5時(shí)引起的滲透率變化， δκ2為裂縫網(wǎng)絡(luò)規(guī)模從5×5×5加密到6×6×6時(shí)引起的滲透率變化.Fig.4 Permeability curves as function of porosity, fracture/soft pore density (a) Aspect ratio of 0.04; (b) Aspect ratio is 0.08.δκ1 is the permeability change caused by the size of the fracture network from 4×4×4 to 5×5×5, and δκ2 is the permeability change caused by the size of the fracture network from 5×5×5 to 6×6×6.

圖5 滲透率隨孔隙度、裂縫/軟孔隙體密度曲線 (a) 縱橫比為0.4； (b) 縱橫比為0.7. 其中δκ1為裂縫網(wǎng)絡(luò)規(guī)模從4×4×4加密到5×5×5時(shí)引起的滲透率變化， δκ2為裂縫網(wǎng)絡(luò)規(guī)模從5×5×5加密到6×6×6時(shí)引起的滲透率變化.Fig.5 Permeability curves as function of porosity, fracture/soft pore body density (a) Aspect ratio of 0.4; (b) Aspect ratio is 0.7. δκ1 is the permeability change caused by the size of the fracture network from 4×4×4 to 5×5×5, and δκ2 is the permeability change caused by the size of the fracture network from 5×5×5 to 6×6×6.