張 磊,姚 軍,孫 海,孫致學(xué)

(中國(guó)石油大學(xué)石油工程學(xué)院,山東青島266580)

頁(yè)巖氣作為非常規(guī)油氣資源日益受到重視，成為常規(guī)能源的替代能源，Javadpour等[1]對(duì)北美9個(gè)油藏152塊巖心進(jìn)行實(shí)驗(yàn)分析，發(fā)現(xiàn)頁(yè)巖的基巖滲透率主要在5.4×10-8μm2左右，90%的頁(yè)巖滲透率小于1.5×10-7μm2。由于頁(yè)巖的低孔低滲等特性，實(shí)驗(yàn)室研究其物理特性難度較大。格子Boltzmann(簡(jiǎn)稱LB)方法是近十幾年來(lái)發(fā)展起來(lái)的介觀模擬方法，可以對(duì)巖心的SEM圖像或CT圖像經(jīng)過(guò)簡(jiǎn)單二值化[2]處理后直接進(jìn)行模擬；學(xué)者們[3-5]利用LB方法計(jì)算多孔介質(zhì)的滲透率，并與實(shí)驗(yàn)結(jié)果和經(jīng)驗(yàn)公式進(jìn)行對(duì)比，擬合結(jié)果比較好。同時(shí)，Javadpour等[1]指出頁(yè)巖氣的主要流動(dòng)孔隙直徑在4～200 nm，頁(yè)巖的孔隙直徑從納米級(jí)到微米級(jí)，氣體在納米孔隙中的流動(dòng)不同于達(dá)西流動(dòng)，氣體分子可以自由地在孔隙表面滑動(dòng)以及與壁面或者其他分子發(fā)生碰撞，因此利用LB方法模擬頁(yè)巖中的流動(dòng)必須考慮頁(yè)巖孔隙的微尺度效應(yīng)。筆者利用Knudsen數(shù)修正的LB模型進(jìn)行模擬，計(jì)算頁(yè)巖的絕對(duì)滲透率。

1 格子Boltzmann方法

格子Boltzmann模型一般由格子模型、平衡態(tài)分布函數(shù)和演化方程組成，Qian等[6]提出的DdQm(d為空間維數(shù)，m為離散速度個(gè)數(shù))模型是格子Boltzmann方法的基本模型。格子Boltzmann方法的計(jì)算過(guò)程由碰撞和遷移兩部分組成，常用的Bhatnagar-Gross-Krook(BGK)模型(不包含外力項(xiàng))的方程為

式中，α=0，1，2，…，N，表示共有N個(gè)速度方向；r為粒子的空間位置；eα為α方向的速度；t為時(shí)間；δ為時(shí)間步長(zhǎng)；fα為離散速度空間α方向上的分布函數(shù)；為離散速度空間的局部平衡態(tài)分布函數(shù)；ρ和u分別為宏觀的粒子密度和速度；cs為格子聲速；ωα為權(quán)系數(shù)。

在每個(gè)格子點(diǎn)上滿足以下關(guān)系：

1.1 基本模型

應(yīng)用最廣泛的模型是D2Q9和D3Q19模型，其模型的速度方向示意圖見(jiàn)圖1。在格子邊長(zhǎng)和時(shí)間步長(zhǎng)均取1的情況下，D2Q9模型各方向速度分布為

D3Q19模型的速度分布：

格子聲速和權(quán)系數(shù)分別為

圖1 D2Q9和D3Q19模型示意圖Fig.1 Schematic diagrams of D2Q9 and D3Q19 models

1.2 松弛時(shí)間和邊界處理

利用格子Boltzmann方法模擬微尺度流動(dòng)時(shí)，須解決兩個(gè)基本問(wèn)題，松弛時(shí)間τ與Knudsen數(shù)(Kn)的關(guān)系和邊界條件的處理。對(duì)松弛時(shí)間τ與Kn的關(guān)系進(jìn)行修正[7]：

為了描述邊界上的滑移流動(dòng)，邊界條件不再采用標(biāo)準(zhǔn)的反彈格式，而是采用標(biāo)準(zhǔn)反彈與鏡面反彈混合的反射格式，見(jiàn)圖2。

式(7)中的α是標(biāo)準(zhǔn)反彈和鏡面反彈的比例系數(shù)，當(dāng)α=1時(shí)，是標(biāo)準(zhǔn)反彈格式；當(dāng)α=0時(shí)，是完全鏡面反彈格式；當(dāng)α=0.5時(shí)，是理想漫反射格式。邊界采取這種簡(jiǎn)單的處理方式是為了在三維實(shí)際數(shù)字巖心中便于實(shí)現(xiàn)。

圖2 粒子邊界反彈示意圖Fig.2 Schematic diagram of particle boundary bounceback

1.3 Kn的計(jì)算

Knudsen數(shù)定義為Kn=λ/H(λ為氣體的分子平均自由程，H為流場(chǎng)的特征長(zhǎng)度)，由于利用格子Boltzmann方法模擬流動(dòng)時(shí)涉及到實(shí)際物理單位與格子單位之間的轉(zhuǎn)化，同時(shí)Kn為無(wú)量綱數(shù)值，在兩種單位下應(yīng)該保持不變。

在格子單位下，Kn和滲透率等參數(shù)主要依賴于格子的劃分，與格子分辨率有關(guān)，同時(shí)存在一些不變量，如格子聲速cs在D2Q9和D3Q19模型中都為1/，即。由動(dòng)理學(xué)理論可知，氣體的平均自由程λ與動(dòng)力黏性系數(shù)μ、壓力p和溫度T之間關(guān)系為

式中，v為運(yùn)動(dòng)黏性系數(shù)；N為特征尺度方向上劃分的格子數(shù)；R為摩爾氣體常數(shù)。

Kn的表達(dá)式為

體外癲癇模型海馬神經(jīng)元中線粒體銜接蛋白Miro1表達(dá)的改變及其意義 … 王英，連亞軍，謝南昌，等 130

當(dāng)選定松弛時(shí)間τ和格子速度c，并且格子劃分確定后，就可以得到Kn，但是在實(shí)際中，Kn與溫度、壓力等有關(guān)，對(duì)于硬球分子，其分子平均自由程為

式中，m為分子質(zhì)量，g；σ為分子直徑，m；M為摩爾分子質(zhì)量，g/mol-1。

Kn的表達(dá)式為

式中，Na為阿伏伽德羅常數(shù)，Na=6.022×1023mol-1。

對(duì)于頁(yè)巖氣，主要成分為甲烷，分子直徑σ約為0.414×10-9m，這樣給定溫度和壓力以及流動(dòng)通道的特征尺寸，就可以計(jì)算出Kn。圖3給出在溫度為350 K，不同特征長(zhǎng)度下Kn與壓力p的關(guān)系。在利用格子Boltzmann方法模擬頁(yè)巖中流體流動(dòng)時(shí)，計(jì)算得到Kn后，根據(jù)上式對(duì)格子模型中的參數(shù)N、δx、δ和τ進(jìn)行設(shè)置，也可以利用式(6)對(duì)τ直接進(jìn)行修正。

圖3 溫度T=350 K和不同特征長(zhǎng)度H下Kn與壓力p的關(guān)系Fig.3 Relation between Kn and p with different H at T=350 K

2 滲透率計(jì)算

2.1 平板模型

模型采用二維平板間的縫隙流動(dòng)，其滲透率理論值為l2/12，l為兩平板之間的距離，見(jiàn)圖4。

圖4 平板模型示意圖Fig.4 Schematic diagram of flat model

利用格子Boltzmann方法計(jì)算其滲透率，出、入口邊界條件不采用周期性邊界，而是采用更符合實(shí)際情況的壓力邊界，不考慮微尺度流動(dòng)的影響，固體邊界采用標(biāo)準(zhǔn)反彈格式，二維平板模型的理論滲透率與LB計(jì)算滲透率的對(duì)比見(jiàn)圖5。由圖5可以看出，LB模擬的滲透率與理論滲透率擬合很好。

圖5 LB滲透率與理論滲透率對(duì)比Fig.5 Comparison between LB permeability and theoretical permeability

圖6 考慮Kn影響的滲透率曲線Fig.6 Permeability curve considering Kn effect

在考察新模型中Kn對(duì)計(jì)算結(jié)果的影響時(shí)，模型中采取在x和y方向格子數(shù)分別取200和20，在不考慮Kn影響的一般模型中，通過(guò)計(jì)算得到格子單位下滲透率K=33.08509。圖6給出的是在不同Kn值影響下計(jì)算得到修正后的滲透率K′及兩者的相對(duì)誤差。從圖6中可以看出，隨著Kn增大，修正后的滲透率K′越小，以相對(duì)誤差5%為界限，在相對(duì)誤差超過(guò)5%時(shí)必須考慮Kn的影響，此時(shí)對(duì)應(yīng)的Kn約在0.14；假如給定溫度為350 K，壓力為2 MPa，此時(shí)可以算出對(duì)應(yīng)的孔隙尺寸約為22.6 nm，也就是對(duì)于實(shí)際的巖心，在該溫度、壓力條件下，當(dāng)孔隙喉道尺寸小于22.6 nm時(shí)，必須考慮Kn對(duì)滲透率的影響。

采用5對(duì)平板并排的模型(圖7)研究Kn對(duì)通道內(nèi)流動(dòng)的影響。平板之間的寬度l分別為10、20、30、40 和50 nm，溫度T=250 K，壓力p=2 MPa，利用式(11)計(jì)算得到其Kn依次為：0.317 1、0.158 55、0.1057、0.079 275和0.063 42，流體流動(dòng)方向由左向右。圖8為在格子時(shí)間t=800時(shí)，在x=L/2處，x方向上的速度分布，其中L為平板長(zhǎng)度，圖中黑線為考慮Kn修正的模型計(jì)算得到的x方向速度U′x，紅線為不考慮Kn的一般模型計(jì)算得到的x方向速度Ux。從圖8中可以看出，隨著寬度增加，這5條通道中流體流動(dòng)的速度最大值顯然是逐漸增大的。對(duì)于同一通道，考慮Knudsen數(shù)Kn的影響后，通道中間部分格子點(diǎn)x方向的速度變大，而且隨著Kn增大(寬度減小)，U′x和Ux之間的差別逐漸增大。圖中右上角小圖是把虛線框中50 nm通道中的曲線放大，可以看出在靠近邊界的部分，考慮Kn影響的x方向的速度是變小的，在數(shù)據(jù)上是靠近邊界的第一個(gè)格子點(diǎn)的速度變小，其他格子點(diǎn)的速度都是變大的，這在其他幾個(gè)通道中也是同樣的結(jié)果。Kn的增大意味著分子平均自由程增大，在相同壓差條件下流動(dòng)更快。

圖7 多平板模型示意圖Fig.7 Schematic diagram of multi-flat model

圖8 x=L/2處的x方向速度分布Fig.8 x-direction velocity distribution at x=L/2

2.2 頁(yè)巖滲透率計(jì)算

圖9 頁(yè)巖SEM掃描圖像Fig.9 SEM image of shale

對(duì)四川盆地某地區(qū)頁(yè)巖巖心進(jìn)行SEM電鏡掃描，得到分辨率為2 μm的圖像(圖9)，由于頁(yè)巖的致密性，孔隙連通性較差，因此不對(duì)其做二維的平面LB流動(dòng)模擬，通過(guò)數(shù)值重構(gòu)方法得到平面巖心切片三維圖像的方法主要有高斯模擬法[8]、模擬退火法[9]、過(guò)程模擬法[10]、多點(diǎn)統(tǒng)計(jì)法[11]和基于馬爾科夫鏈的蒙特卡洛重建法[12]。利用文獻(xiàn)[13]中隨機(jī)建模的方法，將頁(yè)巖二維的SEM圖像重構(gòu)成三維圖像，然后利用經(jīng)過(guò)Knudsen數(shù)修正的格子Boltzmann模型對(duì)三維圖像進(jìn)行模擬。由于進(jìn)行SEM掃描時(shí)，選取的是頁(yè)巖中孔隙比較發(fā)育的部分，得到數(shù)字巖心的孔隙度為14.999 7%，計(jì)算出的滲透率代表孔隙比較發(fā)育部分的滲透率。圖10為流體在數(shù)字巖心中流動(dòng)達(dá)到平衡狀態(tài)時(shí)的速度分布。由圖10可以看出，流體流速在絕大部分孔隙中處于相對(duì)低速，只有在很少的大孔隙中流速相對(duì)較高，圖中速度為格子單位下的速度，通過(guò)計(jì)算得到該頁(yè)巖的絕對(duì)滲透率為0.0185 μm2。

圖10 頁(yè)巖數(shù)字巖心流體速度分布場(chǎng)Fig.10 Velocity distribution field in shale digital core

3 結(jié)束語(yǔ)

利用格子Boltzmann模擬頁(yè)巖中的流體流動(dòng)，Knudsen越大，得到的修正滲透率與原滲透率差別越大。在已知壓力和溫度，一定相對(duì)誤差標(biāo)準(zhǔn)下，可以計(jì)算出模型必須考慮Knudsen效應(yīng)的孔隙最大寬度。Knudsen效應(yīng)使得通道中間部分流體速度增大，而在邊界附近流體速度會(huì)減小。用格子Boltzmann模擬方法計(jì)算滲透率相對(duì)于試驗(yàn)室物理實(shí)驗(yàn)方法有很大的便捷性。

[1] JAVADPOUR F，F(xiàn)ISHER D.UNSWORTH M.Nanoscale gas flow in shale gas sediments[J].Journal of Canadian Petroleum Technology，2007，46(10)：55-61.

[2] 趙秀才，姚軍，房克榮.合理分割巖心微觀結(jié)構(gòu)圖像的新方法[J].中國(guó)石油大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，33(1)：64-67.ZHAO Xiu-cai，YAO Jun，F(xiàn)ANG Ke-rong.A new reasonable segmentation method for microstructure image of reservoir rock[J].Journal of China University of Petroleum(Edition of Natural Science)，2009，33(1)：64-67.

[3] 朱益華，陶果，方偉，等.3D多孔介質(zhì)滲透率的格子Boltzmann模擬[J].測(cè)井技術(shù)，2008，32(1)：25-28.ZHU Yi-hua，TAO Guo，F(xiàn)ANG Wei，et al.Lattice Boltzmann simulation of permeability in 3D porous medium[J].Well Logging Technology，2008，32(1)：25-28.

[4] DEGRUYTER W，BURGISSER A，BACHMANNET O，et al.Synchrotron X-ray microtomography and lattice Boltzmann simulations of gas flow through volcanic pumices[J].Geosphere，2010，6(5)：470-481.

[5] JIN G，PATZEK T W，SILIN D B.Direct prediction of the absolute permeability of unconsolidated and consolidated reservoir rock[R].SPE 90084，2004.

[6] QIAN Y，DHUMIERES D，LALLEMAND P.Lattice BGK models for Navier-Stokes equation[J].Europhysics Letters，1992，17(6)：479.

[7] GUO Z，ZHAO T，SHI Y.Physical symmetry，spatial accuracy，and relaxation time of the lattice Boltzmann equation for microgas flows[J].Journal of Applied Physics，2006，99，074903.

[8] JOSHI M.A class of stochastic models for porous media[D].Lawrence：Universtiy of Kansas，1974.

[9] HAZLETT R D.Statistical characterization and stochastic modeling of pore networks in relation to fluid flow[J].Mathematical Geology，1997，29(6)：801-822.

[10] BRYANT S，BLUNT M.Prediction of relative permeability in simple porous media[J].Physical Review A，1992，46(4)：2004-2012.

[11] OKABE H，BLUNT M J.Prediction of permeability for porous media reconstructed using multiple-point statistics[J].Physical Review E，1997，55：1959-1978.

[12] WU K，VAN DIJKE M I J，COUPLES G D，et al.3D stochastic modelling of heterogeneous porous media-applications to reservoir rocks[J].Transport in Porous Media，2006，65(3)：443-467.

[13] 趙秀才.數(shù)字巖心及孔隙網(wǎng)絡(luò)模型重構(gòu)方法研究[D].青島：中國(guó)石油大學(xué)石油工程學(xué)院，2009.ZHAO Xiu-cai.Numerical rock construction and pore network extraction[D].Qingdao：School of Petroleum Engineering in China University of Petroleum，2009.