張思勤, 汪志明, 洪凱, 程青松, 曾泉樹, 趙巖龍, 郭肖

(1.中石化石油工程機械有限公司第四機械廠, 湖北 荊州 434022; 2.中國石油大學(北京)油氣資源與探測國家重點實驗室, 北京 102249; 3.中國石油南方石油勘探開發(fā)有限責任公司, 廣東 廣州 510240;4.中國石油天然氣股份有限公司塔里木油田分公司, 新疆 庫爾勒 841000)

0 引 言

與國外相比,中國頁巖氣儲層物性差、各向異性突出,孔隙結(jié)構(gòu)存在多尺度特性。對頁巖氣儲層進行滲流行為模擬時達西定律不再適用,常規(guī)的解析方法也較難對其滲流性質(zhì)進行準確模擬[1-4]。此外,由于其孔隙結(jié)構(gòu)的多尺度性和組分的復雜性,較難對其物性參數(shù)進行定量描述,提出一種微觀圖像構(gòu)建數(shù)字巖心與Image J軟件[5-6]相結(jié)合的方法定量提取孔隙信息。采用MATLAB圖像處理技術(shù)和Image J軟件對頁巖二維切片圖像進行處理,根據(jù)所提取的孔隙結(jié)構(gòu)信息重建頁巖氣儲層的三維數(shù)字巖心,然后采格子Boltzmann方法[7-8]對所建的三維數(shù)字巖心進行滲流行為的模擬,并分析孔隙結(jié)構(gòu)對滲透率的影響。

1 格子Boltzmann方法

格子Boltzmann方法是一種新的模擬流體的方法,其模型主要包含流體粒子的離散速度集合、格子結(jié)構(gòu)和演化方程3個要素。它描述了具有離散速度的流體粒子分布函數(shù)在1個固定格點上的運動過程。

粒子分布函數(shù)fi(x,t)的時間演化方程為

(1)

式中,fi(x,t)表示在時刻t、位置x處的節(jié)點上第i個方向的粒子分布函數(shù);x是格子上的一個格點;{ci:i=1,2,…,b}是流體粒子的離散速度集合;δt是離散時間步長;t為當前時間步;τ為松弛時間。分布函數(shù)的演化過程主要分為碰撞和傳播2步,在碰撞過程中,通過每個格點的質(zhì)量和動量重新分布;在傳播過程中,粒子向鄰近點流動。

本文采用的格子模型為D3Q19模型,構(gòu)造格子Boltzmann模型的關(guān)鍵在于選取合適的平衡態(tài)分布函數(shù),離散速度的個數(shù)太少可能導致某些物理量不滿足守恒定律,離散速度的個數(shù)太多又會造成計算上的浪費,D3Q19模型其平衡態(tài)分布函數(shù)

(2)

(3)

D3Q19模型(見圖1)的速度配置為

(4)

圖1 D3Q19格子模型示意圖

根據(jù)質(zhì)量守恒定律和動量守恒定律原則,可以得到宏觀物理參數(shù)與微觀粒子分布函數(shù)的關(guān)系式

(5)

(6)

(7)

(8)

2 數(shù)字巖心重建

通過MATLAB圖像處理技術(shù)對頁巖氣儲層二維巖心切片進行閾值分割和二值化處理,得到巖心切片的二值化圖像;采用Image J軟件對二值圖像進行信息提取,得到孔隙分布個數(shù)及其半徑大小。根據(jù)所提取的孔隙信息將按孔隙半徑劃分為2個孔隙,小于1 μm的為小孔隙,大于1 μm的為大孔隙。以此為基礎(chǔ)在切片圖像中選擇具有代表性的一部分大孔隙和小孔隙進行三維數(shù)字巖心的重建。

2.1 頁巖氣儲層信息提取

重建三維數(shù)字巖心需要3張相互垂直的獨立的二維圖像。采用四參數(shù)隨機生長方法,故只需要1張頁巖氣儲層二維切片圖像即可;頁巖儲層二維巖心切片取自數(shù)字巖心物理實驗室(見圖2)。

圖2 頁巖SEM圖像與二值圖像

在圖2中選取具有代表性的一部分大孔隙和小孔隙用于重建三維數(shù)字巖心的基礎(chǔ)圖像,應(yīng)用MATLAB計算機圖像處理技術(shù)對圖像進行處理,并對圖像孔隙結(jié)構(gòu)特征進行分析和統(tǒng)計,其結(jié)果見表1。

表1 頁巖孔隙結(jié)構(gòu)參數(shù)

將表征單元體作為選擇巖體介質(zhì)類型的依據(jù),表征單元體是指多孔介質(zhì)的一個控制體積,它包含足夠多的微小孔隙,尺寸遠大于單個孔的尺寸,但又遠小于宏觀流動區(qū)域的尺度,故相比于宏觀滲流區(qū)域,REV可以看成是一個點。其概念類似于流體連續(xù)介質(zhì)假定中流體質(zhì)點的概念。以大孔隙頁巖切片圖像對頁巖REV進行說明,圖3(右)為采用四參數(shù)隨機方法重建的頁巖切片圖像。圖4為頁巖切片重建圖像與原圖自相關(guān)函數(shù)和變差函數(shù)對比圖。重建圖像與原圖的自相關(guān)函數(shù)和變差函數(shù)吻合性較好,驗證了四參數(shù)隨機方法的可行性。自相關(guān)函數(shù)圖中,當測量單元尺寸L>20時就能滿足條件。故在后期重建過程中取L=50。

圖3 大孔隙頁巖巖心切片重構(gòu)圖像與原圖對比圖

圖4 頁巖切片自相關(guān)函數(shù)與變差函數(shù)對比圖

圖5 頁巖數(shù)字巖心的孔隙半徑頻率分布

圖5為頁巖切片圖像的孔隙半徑頻率分布圖。圖6是2種對孔隙大小的劃分標準,一種是Choquette and Pray于1970將孔隙劃分為5種基本類型;另一種是Loucks等2012年將孔隙劃分為3種基本類型?？倛D像中孔隙主要以納米孔隙和微孔隙為主,中孔和大孔較少。以2種不同孔隙重建數(shù)字巖心。

圖6 孔隙尺寸劃分方法

2.2 重建方法

建立三維數(shù)字巖心的方法分為物理實驗方法和數(shù)值重建方法兩大類,物理實驗方法重構(gòu)的三維圖像效果較好,但其價格昂貴,不便于大規(guī)模使用,數(shù)值重建方法一般是以1張或數(shù)張二維切片圖像為基礎(chǔ)結(jié)合一定的數(shù)學方法重建數(shù)字巖心。要對所建的頁巖氣儲層數(shù)字巖心采用格子Boltzmann方法進行流動模擬,采用與格子Boltzmann方法緊密結(jié)合的四參數(shù)隨機生長方法建立頁巖氣儲層的三維數(shù)字巖心。

圖7 壓力分布偏差與出口端速度剖面(Kn=0.194)

采用四參數(shù)隨機生長法(Quartet Structure Generation Set, QSGS)構(gòu)造多孔介質(zhì)微觀REV結(jié)構(gòu)[9-10]。首先建立孔隙和基質(zhì)的三維多孔介質(zhì),令基質(zhì)為生長相,孔隙為非生長相,初始相全為孔隙。

(1) 在構(gòu)造網(wǎng)格上以概率pc1隨機布置第1個生長相的生長核,pc1要小于該相的體積分數(shù)。

(2) 對于每一個生長核,按照不同方向上給定的生長概率pi生長,即第i個方向的相鄰單元的狀態(tài)以概率pi變?yōu)榈?相態(tài)(i=1,2,…,18)。

(3) 重復步驟2,直至非生長相(孔隙)達到給定的孔隙率n。

(4) 最后,網(wǎng)格中沒有被生長相占據(jù)的位置即為非生長相單元。

2.3 重建的3D數(shù)字巖心

采用隨機生長方法得到小孔隙三維數(shù)字巖心,再用四參數(shù)隨機生長方法重建大孔隙三維數(shù)字巖心,基于所建立的大孔隙和小孔隙三維數(shù)字巖心為基礎(chǔ),進行模型的組合。組合表達式為

Π=Π1+Π2

(9)

式中,Π為總模型中參數(shù);Π1為大孔隙模型中參數(shù);Π2為小孔隙模型中參數(shù)。大孔隙和小孔隙的基礎(chǔ)模型中,均只包含基質(zhì)和孔隙2部分,其組合類型有4種。

(10)

在組合過程中,先放進模型的組分會有部分被后放進去的所覆蓋,在實際建模過程中,結(jié)合各組分的特點,根據(jù)各種組分的重要性,設(shè)定優(yōu)先級。在組合過程中,設(shè)定孔隙度的優(yōu)先級高于基質(zhì),即孔隙度后組合到總模型中,組合完成后,得到的三維數(shù)字巖心。

3 頁巖氣儲層微觀流動模擬

3.1 模型驗證

選取單孔隙通道對LBM模型進行驗證,采用D2Q9模型模擬二維圓管Poiseuille流動。頁巖氣儲層孔隙小,其滲流行為屬于微尺度流動[11-12],對于微尺度流動要解決2個基本問題:一是確定松弛時間和Knudsen數(shù)之間的關(guān)系;另一個是微尺度流動的邊界條件處理。引入Knudsen數(shù)與松弛時間的關(guān)系式

(11)

通過引入Knudsen數(shù),使格子Boltzmann方法模擬微尺度流動成為可能?？紤]到氣體的壓縮性,對松弛參數(shù)進一步修正,采用Nie等提出的修正松弛參數(shù)τ,其修正表達式

(12)

圖7給出了在Kn=0.388情況下微通道中心處壓力分布偏差和出口端速度分布。通過與Frederik、IP-DSMC、Slip NS、DSMC等的模擬結(jié)果進行對比分析,發(fā)現(xiàn)模擬結(jié)果與Frederik模擬結(jié)果和Slip NS的結(jié)果符合性較好,但與IP-DSMC和DSMC等的模擬結(jié)果有一定的偏差。分析原因主要是當Kn=0.388時,此時的流動為弱過渡流,采用的是單精度的模型,需要引進更加復雜精度的模型。在研究的Kn數(shù)范圍內(nèi),該精度的模型模擬的結(jié)果還是有效的,在一定程度上驗證了模型和方法的準確性。

3.2 頁巖氣儲層數(shù)字巖心滲流特性模擬

采用格子Boltzmann方法對3種不同類型的頁巖氣儲層數(shù)字巖心的孔隙度和滲透率進行模擬計算,選取D3Q19模型,出入口邊界采用壓力邊界,初始格子速度都為0,初始格子密度都為0,可用格子Boltzmann方法計算所建數(shù)字巖心的固有滲透率[13-15],重構(gòu)頁巖儲層三維數(shù)字巖心的孔隙度和滲透率結(jié)果見表2。隨著孔隙度增大,所重建數(shù)字巖心的滲透率明顯增大,組合模型的滲透率不僅分別大于大孔隙數(shù)字巖心的滲透率和小孔隙數(shù)字巖心的滲透率,而且遠大于兩者總和,這說明小孔隙的存在可以大大地提高數(shù)字巖心的連通性,可以極大地提高整體數(shù)字巖心的滲透率。

(13)

式中,K為格子滲透率;N為孔隙空間的格子總數(shù);n為格子方向,此時n=18。

孔隙度的計算式為

(14)

表2 數(shù)字巖心物性參數(shù)與模擬結(jié)果對比

4 結(jié) 論

(1) 隨機生長方法重建數(shù)字巖心過程簡單、成本低,只需要1張二維切片圖像,且與其他重建方法相比,采用隨機生長方法和組合法配合重建數(shù)字巖心能夠極大地提高計算效率。

(2) 通過MATLAB圖像處理技術(shù)和Image J軟件配合,統(tǒng)計出頁巖氣儲層孔隙主要以納米孔和微米孔為主,通過對自相關(guān)函數(shù)和變差函數(shù)分析,發(fā)現(xiàn)四參數(shù)隨機方法重建圖像能較好反映原圖性質(zhì)。

(3) 滲透率隨孔隙度的增大而增大,且組合數(shù)字巖心的滲透率分別大于大孔隙數(shù)字巖心和小孔隙數(shù)字巖心的滲透率,且遠大于二者的總和,小孔隙的存在可以極大地提高儲層的連通性。

參考文獻:

[1] 李建忠, 董大忠, 陳更生, 等. 中國頁巖氣資源前景與戰(zhàn)略地位 [J]. 天然氣工業(yè), 2009(5): 11-16, 134.

[2] 陳尚斌, 朱炎銘, 王紅巖, 等. 中國頁巖氣研究現(xiàn)狀與發(fā)展趨勢 [J]. 石油學報, 2010(4): 689-694.

[3] 張小濤, 吳建發(fā), 馮曦, 等. 頁巖氣藏水平井分段壓裂滲流特征數(shù)值模擬 [J]. 天然氣工業(yè), 2013(3): 47-52.

[4] 段永剛, 李建秋. 頁巖氣無限導流壓裂井壓力動態(tài)分析 [J]. 天然氣工業(yè), 2010(10): 26-29, 116.

[5] 王銀改. Image J軟件在檢驗醫(yī)學圖像分析處理中的應(yīng)用 [J]. 中華檢驗醫(yī)學雜志2005, 28(7), 747-748.

[6] 徐祖新, 郭少斌. 基于氬離子拋光-SEM和ImageJ軟件的頁巖儲層孔隙結(jié)構(gòu)分析——以中揚子地區(qū)陡山沱組為例 [J]. 大慶石油學院學報, 2014(4): 45-51.

[7] 何雅玲, 王勇, 李慶. 格子Boltzmann方法的理論及應(yīng)用 [M]. 北京: 科學出版社, 2009.

[8] 郭照立, 鄭楚光. 格子Boltzmann方法的原理及應(yīng)用 [M]. 北京: 科學出版社, 2008.

[9] 李仁民, 劉松玉, 方磊, 等. 采用隨機生長四參數(shù)生成法構(gòu)造黏土微觀結(jié)構(gòu) [J]. 浙江大學學報: 工學版, 2010, 44(10): 1897-1901.

[10] Wang M, Wang J, Pan N, et al. Mesoscopic Predictions of the Effective Thermal Conductivity for Micro-scale Random Porous Media [J]. Phys Rev E 75, 036702, 2007.

[13] Chai Zhenhua, Guo Zhaoli, Shi Baochang, et al. Predictionof Permeability in Porous Media with Multi-relaxation-time Lattice Boltzmann Method [J]. Journal of Engineering Thermo Physics, 2010, 31(1): 107-109.

[14] 朱益華, 陶果, 方偉, 等. 3D多孔介質(zhì)滲透率的格子Boltzmann模擬 [J]. 測井技術(shù), 2008, 32(1): 25-28.