基于格子Boltzmann方法的3D數(shù)字巖心滲流特征分析

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

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

0 引 言

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

1 格子Boltzmann方法

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

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

(1)

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

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

(2)

(3)

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

(4)

圖1 D3Q19格子模型示意圖

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

(5)

(6)

(7)

(8)

2 數(shù)字巖心重建

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

2.1 頁巖氣儲(chǔ)層信息提取

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

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

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

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

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

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

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

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

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

圖6 孔隙尺寸劃分方法

2.2 重建方法

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

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

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

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

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

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

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

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

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

Π=Π1+Π2

(9)

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

(10)

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

3 頁巖氣儲(chǔ)層微觀流動(dòng)模擬

3.1 模型驗(yàn)證

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

(11)

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

(12)

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

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

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

(13)

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

孔隙度的計(jì)算式為

(14)

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

4 結(jié) 論

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

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

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

參考文獻(xiàn):

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

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

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

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

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

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

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

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

[9] 李仁民, 劉松玉, 方磊, 等. 采用隨機(jī)生長四參數(shù)生成法構(gòu)造黏土微觀結(jié)構(gòu) [J]. 浙江大學(xué)學(xué)報(bào): 工學(xué)版, 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]. 測(cè)井技術(shù), 2008, 32(1): 25-28.