利用格子Boltzmann方法計算頁巖滲透率

張 磊,姚 軍,孫 海,孫致學

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

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

1 格子Boltzmann方法

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

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

在每個格子點上滿足以下關系：

1.1 基本模型

應用最廣泛的模型是D2Q9和D3Q19模型，其模型的速度方向示意圖見圖1。在格子邊長和時間步長均取1的情況下，D2Q9模型各方向速度分布為

D3Q19模型的速度分布：

格子聲速和權系數分別為

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

1.2 松弛時間和邊界處理

利用格子Boltzmann方法模擬微尺度流動時，須解決兩個基本問題，松弛時間τ與Knudsen數(Kn)的關系和邊界條件的處理。對松弛時間τ與Kn的關系進行修正[7]：

為了描述邊界上的滑移流動，邊界條件不再采用標準的反彈格式，而是采用標準反彈與鏡面反彈混合的反射格式，見圖2。

式(7)中的α是標準反彈和鏡面反彈的比例系數，當α=1時，是標準反彈格式；當α=0時，是完全鏡面反彈格式；當α=0.5時，是理想漫反射格式。邊界采取這種簡單的處理方式是為了在三維實際數字巖心中便于實現。

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

1.3 Kn的計算

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

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

式中，v為運動黏性系數；N為特征尺度方向上劃分的格子數；R為摩爾氣體常數。

Kn的表達式為

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

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

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

Kn的表達式為

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

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

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

2 滲透率計算

2.1 平板模型

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

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

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

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

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

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

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

圖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 頁巖滲透率計算

圖9 頁巖SEM掃描圖像Fig.9 SEM image of shale

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

圖10 頁巖數字巖心流體速度分布場Fig.10 Velocity distribution field in shale digital core

3 結束語

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

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

[2] 趙秀才，姚軍，房克榮.合理分割巖心微觀結構圖像的新方法[J].中國石油大學學報：自然科學版，2009，33(1)：64-67.ZHAO Xiu-cai，YAO Jun，FANG 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多孔介質滲透率的格子Boltzmann模擬[J].測井技術，2008，32(1)：25-28.ZHU Yi-hua，TAO Guo，FANG 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] 趙秀才.數字巖心及孔隙網絡模型重構方法研究[D].青島：中國石油大學石油工程學院，2009.ZHAO Xiu-cai.Numerical rock construction and pore network extraction[D].Qingdao：School of Petroleum Engineering in China University of Petroleum，2009.