基于格子Boltzmann方法的致密氣藏微尺度效應(yīng)研究

吳子森，董平川，袁忠超，張雪嬌，曹耐，楊書

（1.中國石油大學(xué)（北京）石油工程教育部重點(diǎn)實(shí)驗(yàn)室，北京 102249；2.中海油研究總院開發(fā)研究院，北京 100027）

基于格子Boltzmann方法的致密氣藏微尺度效應(yīng)研究

吳子森1，董平川1，袁忠超2，張雪嬌1，曹耐1，楊書1

（1.中國石油大學(xué)（北京）石油工程教育部重點(diǎn)實(shí)驗(yàn)室，北京 102249；2.中海油研究總院開發(fā)研究院，北京 100027）

致密儲層孔隙結(jié)構(gòu)復(fù)雜，氣體在致密孔喉中的流動存在微尺度效應(yīng)，宏觀流動規(guī)律難以準(zhǔn)確描述其滲流特征。針對致密多孔介質(zhì)孔隙主要分布在微納米尺度的特點(diǎn)，建立修正克努森數(shù)（Kn）、固體邊界處考慮鏡面反彈（邊界滑移效應(yīng)）的格子Boltzmann模型。通過模擬壓差驅(qū)動下的二維平板流動，驗(yàn)證了模型的正確性，并分析了克努森數(shù)對流動速度的影響?；陔S機(jī)生長四參數(shù)生成法，對致密儲層的二維微觀孔隙結(jié)構(gòu)進(jìn)行了重構(gòu)，利用修正格子Boltzmann模型進(jìn)行流動模擬。結(jié)果表明：氣體在致密儲層中的流動存在微尺度效應(yīng)，滑脫效應(yīng)使得通道中間部分流體速度增大；在一定的壓力梯度下，滲透率隨著克努森數(shù)的增加而呈線性增加；克努森數(shù)不變時，滲透率隨著平均壓力倒數(shù)的增加而呈線性增加，即隨著平均壓力的增加，巖心的絕對滲透率減小。

致密砂巖；微尺度效應(yīng)；格子Boltzmann方法；克努森數(shù)；滲透率

隨著非常規(guī)油氣資源日益受到重視，致密氣藏的勘探開發(fā)逐漸成為熱點(diǎn)［1-3］。由于致密儲層巖性致密，滲透率極低，實(shí)驗(yàn)時對設(shè)備和技術(shù)人員要求較高，常規(guī)的實(shí)驗(yàn)方法難以保證實(shí)驗(yàn)數(shù)據(jù)的準(zhǔn)確性和精度。作為一種介觀模擬方法，格子Boltzmann方法可以對處理過的多孔介質(zhì)圖像進(jìn)行流動仿真模擬［4-6］，進(jìn)而研究流體在多孔介質(zhì)中的流動規(guī)律?；诟褡覤oltzmann方法，張?jiān)频龋?-10］計(jì)算得到多孔介質(zhì)滲透率，并與實(shí)驗(yàn)結(jié)果進(jìn)行擬合，取得較好的效果。致密儲層的孔徑尺度從納米級到微米級，存在多尺度效應(yīng)。王勇杰等［11］研究指出，氣體在致密多孔介質(zhì)中流動時，存在明顯的滑脫效應(yīng)。氣體分子與壁面或者其他分子發(fā)生碰撞，產(chǎn)生滑脫效應(yīng)，使得氣測滲透率大于多孔介質(zhì)的絕對滲透率，因此，在研究致密氣藏時，需要考慮孔隙的微尺度效應(yīng)。前人在研究致密氣藏中氣體的微尺度效應(yīng)時，多采用實(shí)驗(yàn)手段，數(shù)值模擬研究較少，且未考慮微尺度效應(yīng)。筆者基于格子Boltzmann方法，采用隨機(jī)生長四參數(shù)生成法，對儲層的二維微觀孔隙結(jié)構(gòu)進(jìn)行了重構(gòu)，通過修正克努森數(shù)（Kn），考慮到氣體的滑脫效應(yīng)，模擬單相氣體在致密儲層中的流動，并分析了滲透率與克努森數(shù)、平均壓力之間的關(guān)系。

1 格子Boltzmann模型

格子Boltzmann模型（LBM）一般由格子模型、平衡態(tài)分布函數(shù)和演化方程組成?；谠撃Ｐ?，筆者建立了致密氣藏單相格子Boltzmann模型，模擬致密氣藏的氣相流動。模型計(jì)算過程由碰撞步和遷移步兩部分組成，格子Boltzmann方程形式如下［12］：

Boltzmann方程中碰撞項(xiàng)采用單松弛時間形式：

對D2Q9模型［11］，平衡態(tài)分布函數(shù)可表示為

式中：fi（x，t）為t時刻x點(diǎn)處的密度分布函數(shù)；ei為離散速度向量，e0=（0，0），e1=（1，0），e2=（0，1），e3=（-1，0），e4=（0，-1），e5=（1，1），e6=（-1，1），e7=（-1，-1），e8=（1，-1）；t為時間；Δt為時間步長；τ為松弛時間，與流體黏度有關(guān)；ρ為流體的宏觀密度，g/cm3；wi為權(quán)系數(shù)；c為格子速度，一般取值為1；u為流體的宏觀速度，m/s；cs為格子聲速；下標(biāo)i表示離散速度方向 （i= 0，1，…，8）。

1.1 松弛時間和邊界處理

確定碰撞算子的τ，Kn以及微尺度流動中的邊界處理，是格子Boltzmann方法模擬微尺度流動時需要解決的關(guān)鍵問題。宏觀流動中，可用黏性系數(shù)確定τ；對于微尺度流動，特征參數(shù)是克努森數(shù)，在格子Boltzmann模型中，通過引入Kn，進(jìn)而可修正τ。

對D2Q9模型，修正的松弛時間計(jì)算公式為［13-17］

邊界滑移是微尺度流動的重要特征，采取有效準(zhǔn)確的邊界處理方法，才能正確地計(jì)算邊界滑移速度，實(shí)現(xiàn)真實(shí)氣體與固體之間的相互作用。本文采用將無滑移的反彈格式和鏡面反彈按照一定比例組合而成的混合邊界處理格式。

式中：α為標(biāo)準(zhǔn)反彈和鏡面反彈的彈回比例系數(shù)。

α∈［0，1］。α=1時，表示標(biāo)準(zhǔn)反彈格式；α=0時，表示鏡面反彈格式；α=0.5時，表示理想的漫反射。Tang等［16］研究發(fā)現(xiàn)，α=0.7時的計(jì)算結(jié)果與理論值吻合度較高。1.2 克努森數(shù)

在微尺度流動中，正確計(jì)算Kn是微尺度模擬流動的關(guān)鍵。

在實(shí)際中，Kn不僅與流動通道的特征尺寸有關(guān)，還與壓力、溫度有關(guān)。對于硬球分子：

式中：m為分子質(zhì)量，g；R為氣體常數(shù)，R=8.314 Pa·m3/（K·mol）；T為溫度，K；p為壓力，Pa；M為分子摩爾質(zhì)量，g/mol；d為分子直徑，m。

將式（7）代入式（6），整理可得：

式中：NA為阿伏伽德羅常數(shù)（NA=6.022×1023mol-1）。

致密氣的主要成分為甲烷，d約為0.414×10-9m，將給定溫度和壓力以及流動通道的特征尺寸帶入式（8），即可算出Kn值。在T為350 K、H分別為10，20，50，100 nm時，計(jì)算得到Kn與p的關(guān)系（見圖1）。計(jì)算得到Kn后，根據(jù)式（4）對τ進(jìn)行修正，進(jìn)而利用修正的格子Boltzmann模型模擬致密氣流動。

圖1 Kn-p關(guān)系

2 格子Boltzmann模型驗(yàn)證

為了驗(yàn)證格子Boltzmann方法和模型的正確性，基于格子Boltzmann方法，模擬壓力梯度驅(qū)動下二維平板間的泊肅葉流動，其滲透率理論值為l2/12（其中，l為兩平板間的寬度，m）。

計(jì)算模型采用D2Q9模型，格子步長Δx=Δy=1；時間步長Δt=1；不考慮微尺度流動的影響，松弛時間τ= 1；出口邊界和入口邊界采用定壓邊界，上下固體邊界采用標(biāo)準(zhǔn)反彈格式反彈處理。在x/l=0.5處截面上，無因次速度（u/umax）沿y方向的分布呈拋物線分布，可以看出模擬結(jié)果與解析解吻合較好（見圖2）。

圖2 x/l=0.5截面上的速度分布

微流動中，存在復(fù)雜的微尺度效應(yīng)，因此利用格子Boltzmann方法研究微流動時，必須使用修正模型和正確的邊界處理格式。圖3為不同Kn對通道內(nèi)流動的影響，流體流動方向由左往右，左右采用定壓邊界，上下采用反彈邊界與鏡面邊界的混合邊界格式，α=0.7。從圖3可以看出，Kn對氣體流速和邊界滑移速度有很大的影響。

圖3 不同Kn的無因次速度剖面

隨著Kn的增加，壁面邊界滑移速度逐漸增大，且通道流速剖面逐漸偏離拋物線分布的規(guī)律。當(dāng)Kn＜0.1時，此時流動區(qū)域處于滑移區(qū)，通道流速剖面呈拋物線型，壁面滑移速度增加幅度不大，這主要是因?yàn)楫?dāng)流動處于滑移區(qū)時，通道寬度遠(yuǎn)遠(yuǎn)大于分子平均自由程，此時分子間的碰撞占主導(dǎo)地位，流動以黏性流為主。當(dāng)0.1＜Kn≤10，隨著Kn的增大，邊界滑移速度增幅明顯，主要是因?yàn)楫?dāng)流動區(qū)域處于過渡區(qū)時，分子與壁面間的碰撞所占比例逐漸增大，流動以努森流為主，邊界滑移效應(yīng)顯著。

3 致密儲層微尺度效應(yīng)

由于致密儲層孔隙結(jié)構(gòu)復(fù)雜，連通性差，常規(guī)的電鏡和CT掃描獲得的巖心圖像難以用于微觀流動模擬。而通過數(shù)值重建方法，可以獲取連通性較好的數(shù)字巖心［18-20］。目前數(shù)值重建方法主要有隨機(jī)生長法和過程法。筆者基于Wang等［21］提出的隨機(jī)生長四參數(shù)生成法（Quartet structure generation set，QSGS），重構(gòu)致密儲層微觀孔喉結(jié)構(gòu)，初始相全為孔隙，固體顆粒為生長相，構(gòu)造的致密氣藏?cái)?shù)字巖心如圖4a所示（黑色表示固體顆粒，白色表示孔隙顆粒）。

采用修正的格子Boltzmann模型，模擬單相氣體在致密氣藏中的流動。在模擬中，左端為入口，右端為出口，左右邊界采用定壓邊界，入口壓力p1=1.01 MPa，出口壓力p2=1.00 MPa，上下采用混合邊界格式，α= 0.7，Kn=0.050 0。圖4b為氣體在數(shù)字巖心中流動達(dá)到平衡狀態(tài)時的速度分布?？梢钥闯觯瑲怏w流速在絕大部分連通孔隙中處于相對低速，只有在很少的孔隙中流速相對較高。

圖4 致密氣藏?cái)?shù)字巖心

計(jì)算不同Kn和不同壓力梯度下通過多孔介質(zhì)的流體流量，然后帶入達(dá)西公式，就可以計(jì)算出巖心的滲透率K［22］。K與Kn的關(guān)系如圖5所示。從圖5可以看出，在一定的壓力梯度下，滲透率隨著克努森數(shù)的增加而呈線性增加。這說明克努森數(shù)越大，氣體分子的平均自由程越大，導(dǎo)致在相同壓差下氣體分子流動更快，氣體滑脫效應(yīng)越明顯。

K與巖心進(jìn)出口平均壓力pˉ之間的關(guān)系見圖6（其中，pˉ=（p1+p2）/2）。從圖6可以看出，在相同的克努森數(shù)下，滲透率隨著平均壓力倒數(shù)的增加而呈線性增加，即隨著平均壓力的增加，巖心的絕對滲透率減小。這主要是因?yàn)槠骄鶋毫υ叫?，氣體密度越小，氣體分子間的碰撞就越少，氣體滑脫效應(yīng)越明顯，滲透率越大。

圖5 一定壓力梯度下K與Kn的關(guān)系（p1=1.01 MPa，p2=1.00 MPa）

圖6 Kn不變時K與pˉ的關(guān)系（Kn=0.050 0，p2=1.00 MPa）

4 結(jié)論

1）建立修正克努森數(shù)、固體邊界處考慮鏡面反彈（邊界滑移效應(yīng)）的格子Boltzmann模型，并驗(yàn)證了模型的正確性。

2）氣體在納米孔隙中的流動存在微尺度效應(yīng)，滑脫效應(yīng)使得通道中間部分流體速度增大。

3）在一定的壓力梯度下，滲透率隨著克努森數(shù)的增加而呈線性增加；克努森數(shù)不變時，滲透率隨著平均壓力倒數(shù)的增加而呈線性增加，即隨著平均壓力的增加，巖心的絕對滲透率減小。

4）利用修正的格子Boltzmann模擬模型，可便捷研究致密氣藏微尺度，作為實(shí)驗(yàn)室物理實(shí)驗(yàn)方法的有效補(bǔ)充。

［1］趙崇鎮(zhèn).新場氣田須五致密氣藏縫網(wǎng)壓裂技術(shù)［J］.石油鉆探技術(shù)，2015，43（6）：70-75.

［2］趙超.致密氣藏多段壓裂水平井非穩(wěn)態(tài)復(fù)合產(chǎn)能模型［J］.斷塊油氣田，2015，22（5）：651-655.

［3］孫元偉，程遠(yuǎn)方，張礦生，等.考慮非達(dá)西效應(yīng)的致密氣藏裂縫參數(shù)優(yōu)化設(shè)計(jì)［J］.石油鉆探技術(shù)，2014，42（6）：87-91.

［4］王晨晨，姚軍，楊永飛，等.基于格子玻爾茲曼方法的碳酸鹽巖數(shù)字巖心滲流特征分析［J］.中國石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2012，36（6）：95-98.

［5］劉學(xué)鋒，孫建孟，王海濤，等.順序指示模擬重建三維數(shù)字巖心的準(zhǔn)確性評價［J］.石油學(xué)報(bào)，2009，30（3）：391-395.

［6］吳子森，董平川，雷剛，等.基于格子Boltzmann方法的油水兩相流動規(guī)律［J］.斷塊油氣田，2016，23（3）：338-341.

［7］張?jiān)?多孔介質(zhì)中流動的格子Boltzmann模擬［D］.北京：中國石油大學(xué)（北京），2011.

［8］張磊，姚軍，孫海，等.利用格子Boltzmann方法計(jì)算頁巖滲透率［J］.中國石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2014，38（1）：87-91.

［9］寧正福，王波，楊峰，等.頁巖儲集層微觀滲流的微尺度效應(yīng)［J］.石油勘探與開發(fā)，2014，41（4）：445-452.

［10］ETREKIN T，KING G R，SCHWERER F C.Dynamic gas slip-page：a unique dual-mechanism approach to the flow of gas in tight formattions［J］.SPE Formation Evaluation，1986，1（1）：43-52.

［11］王勇杰，王昌杰，高家碧.低滲透多孔介質(zhì)中氣體滑脫行為研究［J］.石油學(xué)報(bào)，1995，16（3）：101-105.

［12］郭照立，鄭楚光.格子Boltzmann方法的原理及應(yīng)用［M］.北京：科學(xué)出版社，2009：218-223.

［13］DONGARI N，AGRAWAL A.Analytical solution of gaseous slip flow inlong microchannels［J］.International Journal of Heat and Mass Transfer，2007，50（17/18）：3411-3421.

［14］LIM C Y，SHU C，NIU X D，et al.Application of lattice Boltzmann methodto simulate microchannel flows［J］.Physics of Fluids，2002，14（7）：2299-2308.

［15］ZHANG Y H，QIN R S，SUN Y H，et al.Gas flow in micro channels：a lattice Boltzmann method approach［J］.Journal of Statistical Physics，2005，21（1/2）：257-267.

［16］TANG G H，TAO W Q，HE Y L.Lattice Boltzmann method for simulating gas flow in microchannels［J］.International Journal of Modern Physics C，2004，15（2）：335-347.

［17］PERUMALDA，KRISHNAV，SARVESHG，etal.Numericalsimulation of gaseous microflows by lattice Boltzmann method［J］.International Journal of Recent Trends in Engineering，2009，1（5）：15-20.

［18］趙秀才.數(shù)字巖心及孔隙網(wǎng)絡(luò)模型重構(gòu)方法研究［D］.青島：中國石油大學(xué)（華東），2009.

［19］劉偉，張德峰，劉海河，等.數(shù)字巖心技術(shù)在致密砂巖儲層含油飽和度評價中的應(yīng)用［J］.斷塊油氣田，2013，20（5）：593-596.

［20］李仁民，劉松玉，方磊，等.采用隨機(jī)生長四參數(shù)生成法構(gòu)造黏土微觀結(jié)構(gòu)［J］.浙江大學(xué)學(xué)報(bào)（工學(xué)版），2010，44（10）：1897-1901.

［21］WANG M，PAN N.Numerical analyses of effective dielectric constant of multiphase microporous media［J］.Journal of Applied Physics，2007，101（11）：114102.

［22］楊勝來，魏俊之.油層物理學(xué)［M］.北京：石油工業(yè)出版社，2007：136-138.

（編輯 趙衛(wèi)紅）

Micro-scale effect in tight gas reservoir based on lattice Boltzmann method

WU Zisen1,DONG Pingchuan1,YUAN Zhongchao2,ZHANG Xuejiao1,CAO Nai1,YANG Shu1
(1.MOE Key Laboratory of Petroleum Engineering,China University of Petroleum,Beijing 102249,China;2.CNOOC Research Institute,Beijing 100027,China)

The micro-pore structure of tight gas reservoir is complex,and gas flow in tight reservoir can lead to micro-scale effect which can′t be described by flowing law of macro fluid.Considering the micro-nano scale pores in tight sandstone,a lattice Boltzmann model(LBM)amended by Knudsen number has been established with a specular bounceback condition(slippage effect) on the solid boundary.The correctness of the new model is verified by the simulation of two parallel plates flow forced by differential pressure.The influence of Knudsen number on velocity was analyzed.With the multi-parameter random generation growth method (MRGGM),a 2D digital core of tight reservoir was reconstructured and the fluid flow behavior in this core was studied by the LBM. The results show that there exists micro-scale effect when gas flows in tight reservoir and the predicted flow velocity in the middle of the flat for the novel model will become larger because of the slippage effect.Under a certain pressure gradient,the permeability increases linearly with the increase of Knudsen number;when the Knudsen number is constant,the permeability increases linearly with the reciprocal of the average pressure,which means the absolute permeability of core decreases with the increase of average pressure.

tight sandstone;micro-scale effect;lattice Boltzmann method;Knudsen number;permeability

國家自然科學(xué)基金項(xiàng)目“基于孔隙網(wǎng)絡(luò)模型的多孔介質(zhì)近混相油氣水三相流動模擬研究”（11072268）

TE311

A

10.6056/dkyqt201606022

2016-04-14；改回日期：2016-08-26。

吳子森，男，1984年生，在讀博士研究生，從事油氣田開發(fā)工程研究工作。E-mail：wuzisen@163.com。

吳子森，董平川，袁忠超，等.基于格子Boltzmann方法的致密氣藏微尺度效應(yīng)研究［J］.斷塊油氣田，2016，23（6）：793-796.

WU Zisen，DONG Pingchuan，YUAN Zhongchao，et al.Micro-scale effect in tight gas reservoir based on lattice Boltzmann method［J］. Fault-Block Oil&Gas Field，2016，23（6）：793-796.