曾彥，寧正福*，齊榮榮，黃亮，呂朝輝

1 中國(guó)石油大學(xué)(北京)石油工程教育部重點(diǎn)實(shí)驗(yàn)室, 北京 102249

2 中國(guó)石油大學(xué)(北京)油氣資源與工程國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 102249

頁(yè)巖氣在納微孔道中的流動(dòng)模擬研究

曾彥1,2，寧正福1,2*，齊榮榮1,2，黃亮1,2，呂朝輝1,2

1 中國(guó)石油大學(xué)(北京)石油工程教育部重點(diǎn)實(shí)驗(yàn)室, 北京 102249

2 中國(guó)石油大學(xué)(北京)油氣資源與工程國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 102249

富有機(jī)質(zhì)頁(yè)巖中存在大量的納米孔隙，連續(xù)介質(zhì)流的Navier-Stokes方程難以準(zhǔn)確描述納米孔隙中的氣體流動(dòng)行為。為此，采用格子Boltzmann方法模擬了氣體在頁(yè)巖納微孔道中的流動(dòng)。結(jié)果表明：頁(yè)巖氣在納微孔道中的流動(dòng)存在顯著的納米尺度現(xiàn)象，表現(xiàn)為氣體的平均流速大于Poiseuille公式計(jì)算值；納微孔道中氣體的稀薄效應(yīng)顯著，孔壁處的氣體速度不為零；且孔道直徑越小，孔壁的滑移速度越大，甚至大于孔道中的氣體速度、流速剖面也由拋物線型轉(zhuǎn)變?yōu)橹麪?。?yè)巖氣在納米孔中流動(dòng)時(shí)發(fā)生“雙滑脫效應(yīng)”，流體粒子與孔壁發(fā)生碰撞，獲得動(dòng)能后反彈進(jìn)入孔道中，增強(qiáng)整個(gè)納微孔道氣體的流動(dòng)?！半p滑脫效應(yīng)”導(dǎo)致頁(yè)巖氣的滲透率大于等效液體滲透率，并引起克氏滲透率圖中，氣體滲透率與平均壓力的倒數(shù)偏離直線關(guān)系。

頁(yè)巖；納米孔；格子Boltzmann方法；流動(dòng)模擬；納米尺度；雙滑脫效應(yīng)

富有機(jī)質(zhì)頁(yè)巖中聚集了大量的天然氣，這使得頁(yè)巖氣成為全球油氣資源勘探開(kāi)發(fā)的熱點(diǎn)[1-2]。前人已經(jīng)采用多種技術(shù)手段來(lái)研究頁(yè)巖儲(chǔ)層的孔隙結(jié)構(gòu)特征[3-6]。采用離子拋光制樣和高分辨率的場(chǎng)發(fā)射掃描電鏡實(shí)驗(yàn)顯示頁(yè)巖分散的有機(jī)質(zhì)中存在大量的納米孔隙，孔徑大小為10～800 nm。盡管頁(yè)巖中存在一些大孔隙，高壓壓汞和低溫氮?dú)馕綄?shí)驗(yàn)結(jié)果表明頁(yè)巖主體孔徑集中在50 nm以?xún)?nèi)。氣體在納米尺度空間中的流動(dòng)規(guī)律不同于宏觀尺度流動(dòng)。在納米尺度下，氣固界面之間的分子間作用力更強(qiáng)，動(dòng)量交換更加劇烈，甚至成為氣體流動(dòng)的主要影響因素[7]。

國(guó)內(nèi)外學(xué)者對(duì)氣體在頁(yè)巖儲(chǔ)層中的滲流行為進(jìn)行了大量的研究。Javadpour[8-9]建立了頁(yè)巖氣在納米孔隙中的傳質(zhì)模型，考慮了頁(yè)巖氣的擴(kuò)散、滑脫和吸附解吸行為，通過(guò)引入視滲透率參數(shù)，修正了頁(yè)巖氣滲流方程。Freeman[10]利用塵氣模型，研究了Knudsen擴(kuò)散對(duì)低滲透致密介質(zhì)流動(dòng)的影響。王瑞等[11]考慮頁(yè)巖氣擴(kuò)散類(lèi)型的差異，分析了擴(kuò)散對(duì)頁(yè)巖氣質(zhì)量流量的貢獻(xiàn)。這些都側(cè)重于宏觀滲流的研究，而對(duì)于頁(yè)巖氣在納米孔隙中的微觀滲流行為還鮮有報(bào)道。本文引入適合宏觀和介觀尺度的格子Boltzmann方法，考慮流體粒子之間、流體粒子與納米孔壁之間的相互作用，從微觀角度研究頁(yè)巖氣在納微孔道中的流動(dòng)行為。

1 頁(yè)巖氣在納米孔隙中的流動(dòng)狀態(tài)

1.1 頁(yè)巖中的納米孔隙

圖1為頁(yè)巖有機(jī)質(zhì)中的納米孔隙。頁(yè)巖樣品來(lái)自四川盆地下寒武統(tǒng)牛蹄塘組以及西北地區(qū)六盤(pán)山盆地下白堊統(tǒng)乃家河組，有機(jī)碳含量為0.5%～9.15%，成熟度為0.6%～2.7%。頁(yè)巖樣品首先采用氬離子拋光技術(shù)進(jìn)行處理，然后在FEITMQuantaTM200F場(chǎng)發(fā)射掃描電鏡下觀察孔隙結(jié)構(gòu)，實(shí)驗(yàn)儀器的相關(guān)參數(shù)和實(shí)驗(yàn)方法見(jiàn)文獻(xiàn)[5]。掃描電鏡結(jié)果顯示，頁(yè)巖有機(jī)質(zhì)孔隙直徑為10～450 nm。

為了和可視化的電鏡圖片結(jié)果進(jìn)行對(duì)比，對(duì)頁(yè)巖樣品進(jìn)行了高壓壓汞實(shí)驗(yàn)。壓汞實(shí)驗(yàn)的驅(qū)替壓力最高達(dá)到200 MPa，測(cè)試孔隙直徑下限為7.35 nm。壓汞孔徑分布曲線表明，頁(yè)巖樣品的孔隙直徑主要為8～100 nm(圖2)。

圖1 頁(yè)巖有機(jī)質(zhì)中的納米孔隙Fig. 1 Nanopores in organic matters

圖2 頁(yè)巖壓汞法孔徑分布曲線Fig. 2 Pore size distribution of shales

1.2 流動(dòng)狀態(tài)

圖3 Knudsen數(shù)與氣體流動(dòng)狀態(tài)Fig. 3 The relationship between Knudsen number and fl ow state

根據(jù)Knudsen數(shù)(氣體平均分子自由程與流動(dòng)空間特征尺度的比值)[12]的大小，將流體在多孔介質(zhì)中的流動(dòng)狀態(tài)劃分為(圖3)：(1)當(dāng)Kn＜0.001時(shí)，為連續(xù)介質(zhì)流。流體可以認(rèn)為是連續(xù)介質(zhì)，氣體滲流服從達(dá)西定律。(2)當(dāng)0.001≤Kn＜0.1時(shí)，為滑脫流。氣體分子之間的碰撞頻率仍然較高，但是稀薄氣體效應(yīng)已不能再忽略，這時(shí)流動(dòng)控制方程仍然是經(jīng)典流體力學(xué)的Navier-Stokes方程，但必須考慮邊界滑移速度和非等溫條件下的溫度跳躍，氣體流動(dòng)服從Klinkenberg方程。(3)當(dāng)0.1≤Kn＜10時(shí)，為過(guò)渡流(或Knudsen流)。此時(shí)氣體分子之間的碰撞頻率和氣體分子與固體之間的碰撞頻率差不多，連續(xù)介質(zhì)假設(shè)已不再成立，氣體流動(dòng)服從Burnett方程。(4)當(dāng)Kn≥10時(shí)，為自由分子流。此時(shí)氣體分子與固體之間的碰撞頻率比分子之間的碰撞頻率高得多，氣體分子之間的碰撞可以忽略，主要采用分子動(dòng)力學(xué)方法進(jìn)行研究。

表1為典型頁(yè)巖氣藏的氣體組成和氣體性質(zhì)參數(shù)。頁(yè)巖氣主要由CH4、C2H6和CO2組成，其摩爾分?jǐn)?shù)分別為80%、10%和10%。圖4為溫度為300 K時(shí)，隨著儲(chǔ)層壓力的變化(0.1～40 MPa)，頁(yè)巖氣在不同孔徑中的Knudsen數(shù)。頁(yè)巖的孔徑分布主要為8～150 nm，因此從頁(yè)巖儲(chǔ)層到井底，頁(yè)巖氣在納米孔隙中的流動(dòng)主要為滑脫流和弱過(guò)渡流。在這兩個(gè)流動(dòng)階段，氣體分子的平均自由程大、密度小，稀薄效應(yīng)明顯，氣體流動(dòng)違背了控制方程中的連續(xù)介質(zhì)假設(shè)，不存在滑脫邊界的Navier-Stokes方程是無(wú)效的。因此，基于達(dá)西定律的經(jīng)典滲流理論不適用于頁(yè)巖氣儲(chǔ)層流動(dòng)模擬。

2 格子Boltzmann方法

基于動(dòng)力學(xué)理論的Boltzmann方程是連接宏觀和微觀的橋梁，可用于從連續(xù)流到自由分子流的跨尺度流動(dòng)。求解Boltzmann方程非常困難，從Boltzmann方程和介觀層次上的格子氣方法發(fā)展而來(lái)的格子Boltzmann方法(LBM)成為一個(gè)較好的解決方案。格子Boltzmann方法將流體顆粒離散成一系列的流體粒子。這些粒子比分子級(jí)別大，但在宏觀上又無(wú)限小，流體粒子在2D離散格子網(wǎng)格組成的毛細(xì)管模型中碰撞和流動(dòng)。在一定的初始和邊界條件下，這些流體粒子之間、以及流體粒子和毛細(xì)管孔壁之間相互作用，但它們共同地遵循Boltzman流動(dòng)方程。以氣體動(dòng)力學(xué)理論為基礎(chǔ)的格子Boltzmann方法已經(jīng)被認(rèn)為是模擬微觀尺度介質(zhì)的有效數(shù)值方法[12-13]。

表1 典型頁(yè)巖氣藏的氣體組成和氣體性質(zhì)參數(shù)Table 1 Components properties of shale gas

圖4 不同孔徑下，頁(yè)巖氣在儲(chǔ)層中的Knudsen數(shù)Fig. 4 The relationship between pressure and Knudsen number under different pore sizes

2.1 計(jì)算模型

Freeman等運(yùn)用塵氣模型(DGM)研究頁(yè)巖儲(chǔ)層中非吸附相的運(yùn)移。DGM可以準(zhǔn)確表征滑脫效應(yīng)和擴(kuò)散效應(yīng)[14]。將DGM控制方程離散，引入格子Boltzmann方法，可以預(yù)測(cè)頁(yè)巖氣在多孔介質(zhì)中的視滲透率?；贛axewll-Stefan擴(kuò)散方程，頁(yè)巖氣的流動(dòng)方程可以表示為[15]：

式中，P是壓力(Pa)；R是普氏氣體常數(shù)，8.314 J/(mol·K)；T是溫度(K)；μ是動(dòng)力黏度(kg/(m·s))；Dk是Knudsen擴(kuò)散系數(shù)(m2/s)。K0是固有滲透率，K0=d232，m2

Knudsen擴(kuò)散系數(shù)可表示為：

式中，d是孔道直徑，m；M是氣體分子量，kg/mol。對(duì)于三維多孔介質(zhì)模型，沿流動(dòng)方向的有效Knudsen擴(kuò)散系數(shù)可表示為：

式中，Lx，Ly，Lz分別為多孔介質(zhì)x，y，z三個(gè)方向的長(zhǎng)度，其中x方向?yàn)榱鲃?dòng)方向。ρin，ρout分別為進(jìn)出口端密度。引入視滲透率Kapp，則有：

2.2 格子Boltzmann模型

本文采用格子Boltzmann-BGK方程D2Q9模型(二維九速模型)對(duì)納微孔道(如圖5)中的頁(yè)巖氣流動(dòng)進(jìn)行了模擬。

圖5 納米孔道物理模型Fig. 5 Physical model of nanopores

通過(guò)BGK近似、時(shí)空離散和速度離散，得到完全離散化的格子Boltzmann-BGK方程[12-13]，即LBM演化方程：

式中，δt表示時(shí)間步長(zhǎng)；fα(r,t)表示t時(shí)刻r處速度為eα的粒子密度分布函數(shù)；為對(duì)應(yīng)時(shí)刻點(diǎn)的局部平衡態(tài)分布函數(shù)；τ為松弛時(shí)間；Fα代表由孔道和外力項(xiàng)引起的體積力，可表示為：

其中，G為外力擾動(dòng)項(xiàng)，Kapp為式(5)求得的視滲透率。Guo和Zhao將平衡態(tài)分布函數(shù)和外力項(xiàng)分布函數(shù)分別表示為：

式中，ρ、u表示流體宏觀密度和流動(dòng)速度；ωα為與離散速度方向矢量模量有關(guān)的權(quán)系數(shù)；cs為格子聲速；eα表示α方向上的粒子速度。

模型的速度配置如下：

由式(7)可知，外力項(xiàng)F中也包含了速度u，將式(7)代入式(12)即可求得速度u。

在給定的離散空間和離散時(shí)間下，式(6)中的松弛時(shí)間計(jì)算公式為[17]：

式中，Ny表示y方向上的網(wǎng)格數(shù)，Kn為給定空間的Knudsen數(shù)。式(14)中Knudsen數(shù)與松弛時(shí)間的關(guān)系為[18-19]。

式中，p0為出口端壓力；Kn0為出口端Knudsen數(shù)。

2.3 邊界條件

在LBM模擬中，對(duì)于固體邊界的處理方法主要有兩種：無(wú)滑移邊界的反彈格式和自由滑移邊界的鏡面反射格式。對(duì)于氣體在頁(yè)巖納微孔道中的流動(dòng)，既不能用簡(jiǎn)單的反彈格式處理，也不能用鏡面反射來(lái)描述氣體與固體壁面的相互作用與動(dòng)量交換。將反彈和鏡面反射相結(jié)合能比較準(zhǔn)確地反映真實(shí)氣體與固體之間的相互作用。在混合格式中，引入彈回比例系數(shù)r，表示粒子與邊界作用后沿原路返回所占的比例。研究表明[20]，r取值為0.7時(shí)，能最佳地描述氣體與固體壁面的相互作用。本次模擬，設(shè)定上下邊界為封閉邊界，采用反彈和鏡面反射的組合格式，且彈回比例系數(shù)為0.7；左右邊界為定壓邊界，左邊界為入口，右邊界為出口，孔道中流體在兩端壓差作用下進(jìn)行流動(dòng)。

定壓邊界采用的是外推邊界條件。假設(shè)物理邊界外還有一層虛擬節(jié)點(diǎn)，物理邊界節(jié)點(diǎn)(i,1)和這層虛擬節(jié)點(diǎn)(i,0)上的分布函數(shù)都參與碰撞遷移過(guò)程。借鑒有限差分的概念，在邊界節(jié)點(diǎn)(i,1)采用中心差分，確定物理邊界節(jié)點(diǎn)上的未知分布函數(shù)。

外推計(jì)算之后，對(duì)所有節(jié)點(diǎn)執(zhí)行遷移步驟。通過(guò)分布函數(shù)計(jì)算進(jìn)出口端壓力，并對(duì)其進(jìn)行修正。使進(jìn)出口端平均壓力分別等于給定的壓力值Pin，Pout。在隨后的碰撞過(guò)程中，邊界節(jié)點(diǎn)上的平衡態(tài)分布函數(shù)需要由修正后的壓力值來(lái)計(jì)算。對(duì)壓力進(jìn)行修正，不僅要保證進(jìn)出口端壓力滿(mǎn)足給定的壓力條件，還使得其壓力剖面符合內(nèi)部流場(chǎng)特征。

(3)語(yǔ)法復(fù)習(xí)目標(biāo)：掌握各項(xiàng)基礎(chǔ)語(yǔ)法的基本概念及重要用法；熟練理解、掌握并能夠在不同的語(yǔ)境中辨別、正確、靈活運(yùn)用各種重要語(yǔ)法功能。

3 模型驗(yàn)證

應(yīng)用LBM之前，必須對(duì)該方法和模型進(jìn)行驗(yàn)證。LBM適用于連續(xù)介質(zhì)流、滑脫流和弱過(guò)渡流。本文首先在連續(xù)介質(zhì)流條件下，通過(guò)LBM的模擬結(jié)果與Poiseuille流動(dòng)的解析解進(jìn)行對(duì)比，說(shuō)明LBM在連續(xù)介質(zhì)流下的正確性。在連續(xù)性假設(shè)失效時(shí)，通過(guò)與他人的解析結(jié)果和數(shù)值結(jié)果進(jìn)行對(duì)比，驗(yàn)證在納米尺度下LBM的正確性。

3.1 連續(xù)介質(zhì)區(qū)的流體流動(dòng)驗(yàn)證

物理模型如圖5，孔道上下為封閉邊界，相距H，孔道中充滿(mǎn)了黏度為μ的流體，在壓力梯度驅(qū)動(dòng)下，流體沿著x方向進(jìn)行流動(dòng)。流場(chǎng)空間劃分為正方形均勻網(wǎng)格，網(wǎng)格數(shù)為4000×400，空間步長(zhǎng)Δx=Δy，時(shí)間步長(zhǎng)Δt=1，松弛時(shí)間τ=1(連續(xù)介質(zhì)流不存在壓縮效應(yīng)，孔道兩端壓力一定時(shí)，沿孔道的Kn數(shù)、松弛時(shí)間為定值)，各單位均采用格子單位。在孔道兩端給定一定的密度，使流體在x方向產(chǎn)生壓力梯度。本次模擬中，進(jìn)口密度為1.0012 kg/m3，出口密度為1 kg/m3。

圖6 連續(xù)介質(zhì)流的格子Boltzmann模擬結(jié)果與Poiseuille公式的解析解對(duì)比Fig. 6 Streamwise velocity at outlet

圖6為連續(xù)介質(zhì)流的出口端速度剖面。連續(xù)介質(zhì)流條件下，由于流體的黏性耦合特征，孔道中流體的流動(dòng)速度剖面為拋物線，即孔道中間的流體速度最大，孔壁處流體速度為零。從圖中可以看出，格子Boltzmann方法的模擬結(jié)果與Poiseuille公式的解析解吻合非常好，這表明格子Boltzmann方法可以準(zhǔn)確地描述連續(xù)介質(zhì)流流體的動(dòng)力學(xué)行為。

3.2 滑脫流和弱過(guò)渡流的流體流動(dòng)驗(yàn)證

對(duì)于連續(xù)性假設(shè)失效或宏觀方程難以描述的現(xiàn)象，采用微觀分子動(dòng)力學(xué)或介觀氣體動(dòng)力學(xué)描述更加恰當(dāng)。為了和他人的結(jié)果[21-24]進(jìn)行對(duì)比，計(jì)算了孔道出口端的流速分布剖面。模型的基本參數(shù)選取如下：x方向網(wǎng)格數(shù)Nx=1100；y方向網(wǎng)格數(shù)Ny=11；出口端密度ρout=1；圖7(a)為滑脫流階段的模擬結(jié)果，其中出口端Knudsen數(shù)Kn0= 0.0194，入口端密度ρin=1.4；圖7(b)為弱過(guò)渡流階段的模擬結(jié)果，其中出口端Knudsen數(shù)Kn0= 0. 194，入口端密度ρin=2。

通過(guò)比較發(fā)現(xiàn)，筆者的模擬結(jié)果與Frederik Verhaeghe的模擬結(jié)果、Arkilic理論解以及DSMC方法的模擬結(jié)果重合性良好，驗(yàn)證了模型及程序的正確性，說(shuō)明在滑移區(qū)和弱過(guò)渡區(qū)，采用LBM方法可以有效模擬氣體在頁(yè)巖納微孔道中的氣體流動(dòng)行為。

4 模擬結(jié)果及討論

下面模擬了頁(yè)巖氣在不同孔徑納微孔道中的流動(dòng)。納米孔孔徑分別為10 nm、20 nm、50 nm、100 nm和1 μm，孔道的徑長(zhǎng)比為0.01，孔道的平均壓力為0.1、0.5、1、2、5和10 MPa，儲(chǔ)層溫度300 K，頁(yè)巖氣藏的氣體組成和氣體性質(zhì)參數(shù)見(jiàn)表1，計(jì)算得到頁(yè)巖氣在納微孔道中流動(dòng)時(shí)對(duì)應(yīng)的Knudsen數(shù)為0.000 58～5.826 0。

4.1 平均流速

為了便于和經(jīng)典流體力學(xué)的Poiseuille流動(dòng)結(jié)果進(jìn)行對(duì)比，定義納微孔道中氣體流動(dòng)的無(wú)因次平均流速為：式中，為采用LBM方法模擬的納微孔道的平均流速；采用Poiseuille公式求得流動(dòng)速度的解析解。

模擬得到的不同直徑孔道的無(wú)因次平均流速與平均壓力的關(guān)系如圖8。從圖中可以看出，高壓下，無(wú)因次平均流速趨于1，納米孔的平均流速與經(jīng)典Poiseuille公式的解析解一致。但當(dāng)壓力＜5 MPa，模擬得到的納微孔道流速大于Poiseuille公式的解析解。壓力低于2 MPa時(shí)，納微孔道的平均流速比經(jīng)典值高1～2個(gè)數(shù)量級(jí)。納微孔道直徑越小，流體流動(dòng)速度和Poiseuille公式解析解的差別越大。低壓下氣體在納米孔的平均流速的模擬結(jié)果顯著高于Poiseuille流動(dòng)解，盡管我們可以將這歸因于“滑脫效應(yīng)”的影響，但如此高的平均流速是非常少見(jiàn)的，分析認(rèn)為這是由于納微孔道中的“雙滑脫效應(yīng)”的促進(jìn)作用，將在下文中進(jìn)行解釋。

圖7 LBM方法模擬滑脫流和過(guò)渡流階段的孔道出口端速度剖面與文獻(xiàn)結(jié)果對(duì)比Fig. 7 Normalized streamwise velocity at different Knudsen number

圖8 不同直徑納米孔的氣體無(wú)因次平均流速Fig. 8 Dimensionless average velocity in different nanopore sizes

4.2 空間流速分布

隨著壓力降低，氣體分子的平均自由程增大，Knudsen數(shù)增大，納微孔道中氣體的流動(dòng)與Poiseuille流的差別更明顯。圖9和圖10為平均壓力為0.1 MPa時(shí)，頁(yè)巖氣分別在1 μm和50 nm納米孔中的空間流速分布(無(wú)因次流速為流動(dòng)速度與孔道出口端中心處流速的比值)。沿著流動(dòng)方向(x方向)，納米孔中心處(y/H=0)的流動(dòng)速度越來(lái)越大。這是因?yàn)樵娇拷隹诙?，孔道?nèi)的壓力越低，氣體密度越小，Knudsen數(shù)越大，氣體稀薄效應(yīng)越嚴(yán)重。為了滿(mǎn)足質(zhì)量守恒，氣體的流動(dòng)速度增大。圖中孔道的上下邊界處(y/H=±0.5)表現(xiàn)出明顯的滑脫效應(yīng)，邊界流速不為零?？讖綖? μm的(圖9)邊界滑移速度小于孔道內(nèi)氣體流速，但孔徑為50 nm的納米孔(圖10)邊界處的流速與孔道中氣體的速度相當(dāng)，這是由于孔徑減小，氣體與孔壁之間的作用力增強(qiáng)的結(jié)果。

圖9 頁(yè)巖氣在1 μm孔中流動(dòng)的空間速度分布Fig. 9 Dimensionless velocity distribution in 1 μm pore

圖10 頁(yè)巖氣在50 nm孔中流動(dòng)的空間速度分布Fig. 10 Dimensionless velocity distribution in 50 nm pore

4.3 流速剖面

圖11為模擬得到的不同孔徑納米孔的出口端流速剖面。納微孔道直徑為1 μm(或者大于1 μm)時(shí)，盡管邊界處的流速不為零，但出口端的流速剖面仍然為Poiseuille流的經(jīng)典拋物線形狀。隨著孔道直徑的降低，由于氣體的滑脫效應(yīng)增強(qiáng)，孔壁附近的滑移流速逐漸增大。當(dāng)孔道直徑減小至50 nm時(shí)，盡管孔道中心速度存在極值，但孔壁滑移速度與孔道中心速度相當(dāng)，孔道出口端的流速剖面趨于均勻?？讖叫∮?0 nm后，孔壁的滑移速度顯著大于孔道中心速度，流動(dòng)速度剖面變?yōu)橹麪?。低壓下，氣體在小尺寸納米孔(＜50 nm)的流動(dòng)行為與在大孔(1 μm或更大)中的流動(dòng)出現(xiàn)根本區(qū)別。

為了進(jìn)一步說(shuō)明納米孔的氣體流動(dòng)與大孔道中氣體流動(dòng)的區(qū)別，我們計(jì)算了出口端氣體流動(dòng)的動(dòng)能剖面(圖12和圖13)。無(wú)因次動(dòng)能為L(zhǎng)BM模擬的動(dòng)能與Poiseuille公式計(jì)算的出口端孔道中心(x/L=1，y/H=0)動(dòng)能的比值。從圖中可以看出，孔道為1 μm時(shí)，LBM的模擬結(jié)果與Poiseuille公式計(jì)算的動(dòng)能(解析法)非常接近，孔道中心處的動(dòng)能幾乎一致；只是在孔壁處，考慮滑脫效應(yīng)的LBM方法模擬的動(dòng)能略大于解析結(jié)果(圖12)。但當(dāng)孔道直徑為50 nm時(shí)，不僅孔壁處的動(dòng)能大于解析法結(jié)果，而且孔道中心的動(dòng)能也顯著高于解析法結(jié)果(圖13)。因此，隨著孔徑的降低，納米孔孔壁處流體粒子的動(dòng)能增大，并且孔道中心處氣體的動(dòng)能也增大。這是由于流體粒子與孔壁之間的相互作用引起的。氣體在孔道中流動(dòng)時(shí)主要發(fā)生兩種基本相互作用：流體粒子間的碰撞，流體粒子與孔壁的碰撞[25]。當(dāng)流體粒子的平均自由程小于孔隙直徑時(shí)，常發(fā)生流體粒子間的碰撞；而對(duì)于那些平均自由程大于孔隙直徑的流體粒子，則常發(fā)生流體粒子與孔隙壁間的碰撞。隨著孔道直徑的降低，流體粒子平均自由程大于孔隙直徑，流體粒子與孔壁的碰撞越劇烈，滑脫現(xiàn)象愈顯著?？讖叫∮?0 nm時(shí)，流體粒子與孔壁的相互作用逐漸占主導(dǎo)地位，成為氣體流動(dòng)的主要影響因素，邊界滑脫速度顯著高于孔道中心的氣體流速，流動(dòng)速度剖面因此轉(zhuǎn)變?yōu)椤皟啥烁咧虚g平緩”的柱塞狀。

圖11 納米孔的出口端流速剖面Fig. 11 Streamwise velocity at different nanopore sizes

圖12 直徑為1 μm孔道的出口端動(dòng)能剖面Fig. 12 Kinetic energy distribution in 1 μm pore

由于孔壁處的流體粒子流速不為零，孔壁處的流體粒子與孔壁發(fā)生非彈性碰撞后獲得了動(dòng)能，并且反彈后的流體粒子在垂直于孔壁方向上的速度不為零。獲得動(dòng)能的流體粒子進(jìn)入孔道中，使孔道中氣體動(dòng)能增大，流動(dòng)速度增大，增強(qiáng)了孔道中心氣體的流動(dòng)(圖14)。正是由于這個(gè)原因，從出口端氣體流動(dòng)的動(dòng)能剖面上可以看出，隨著孔徑的減小，孔道中心(y/H=0)的氣體動(dòng)能大于Poiseuille公式計(jì)算的動(dòng)能。因此，不僅管壁處的流體粒子發(fā)生滑脫，整個(gè)孔道中的流體也受到這種滑脫效應(yīng)的促進(jìn)作用，即所謂的“雙滑脫”現(xiàn)象[26-28]。1941年Klinkenberg[29]首次提出滑脫效應(yīng)，傳統(tǒng)的滑脫效應(yīng)主要用于表征滑脫區(qū)的流動(dòng)特征，只考慮了孔壁處的氣體流速不為零，沒(méi)有考慮到孔壁對(duì)流體粒子在整個(gè)納微孔道中流動(dòng)的促進(jìn)作用。隨著流動(dòng)尺度不斷縮小，傳統(tǒng)的滑脫效應(yīng)越來(lái)越難以準(zhǔn)確表征流體視滲透率。很多學(xué)者[30-31]選擇更高階的Klinkenberg方程來(lái)表征滑脫效應(yīng)。頁(yè)巖氣儲(chǔ)層主要孔徑為幾十個(gè)納米，脫離了傳統(tǒng)滑脫效應(yīng)表征范圍，“雙滑脫”現(xiàn)象對(duì)于納米孔中頁(yè)巖氣的流動(dòng)促進(jìn)作用非常重要。

圖13 直徑為50 nm孔道的出口端動(dòng)能剖面Fig. 13 Kinetic energy distribution in 50 nm pore

圖14 “雙滑脫效應(yīng)”中頁(yè)巖氣體分子與納米孔孔壁的相互作用Fig. 14 The interaction between gas molecules and nanopore

4.4 孔道滲透率

頁(yè)巖氣體在納微孔道中流動(dòng)時(shí)，流體粒子與孔道之間存在強(qiáng)烈的相互作用。納米尺度下氣體的流動(dòng)行為與常規(guī)尺度滲流有所不同。為了研究納微孔道中，納米尺度流動(dòng)對(duì)頁(yè)巖氣滲流的影響，定義無(wú)因次滲透率為：

式中，KLBM為采用LBM模擬的不同平均壓力下的納米孔的氣體滲透率；KPoiseuille為采用Poiseuille公式和達(dá)西定律計(jì)算的等效液體滲透率。

在經(jīng)典的滑脫效應(yīng)理論中，考慮氣體滑脫的氣測(cè)滲透率與孔道平均壓力的倒數(shù)呈直線關(guān)系，或者氣測(cè)滲透率和等效液體滲透率的比值與平均壓力的倒數(shù)為直線。圖15為頁(yè)巖氣在不同孔徑納米孔孔的無(wú)因次滲透率與平均壓力倒數(shù)的關(guān)系(克氏滲透率校正圖)。從圖中可以看出，高壓下(10 MPa)，大孔(孔徑為1 000 nm)中的氣測(cè)滲透率與等效液體滲透率相等，表現(xiàn)為無(wú)因次滲透率為1的直線。但在低壓下，“雙滑脫效應(yīng)”對(duì)納米孔的氣體滲流產(chǎn)生了顯著的影響：納微孔道的氣體滲透率大于等效液體滲透率，并且在克氏滲透率校正圖中，孔道的氣體滲透率與平均壓力的倒數(shù)不再呈直線關(guān)系。隨著孔道平均壓力的降低，無(wú)因次滲透率與平均壓力倒數(shù)的關(guān)系曲線發(fā)生上翹；納微孔道的直徑越小，無(wú)因次滲透率的上翹現(xiàn)象越明顯。平均壓力為2 MPa時(shí)，直徑為50 nm孔的氣體滲透率是等效液體滲透率的2.1倍，直徑為20 nm孔的氣體滲透率是液體滲透率的3.7倍。平均壓力為1 MPa時(shí)，直徑為20 nm孔的氣體滲透率則是液體滲透率的9.8倍。“雙滑脫”現(xiàn)象增強(qiáng)了頁(yè)巖氣在納米孔隙中的滲流，并且孔徑越小、壓力越低，這種增強(qiáng)作用越明顯。在頁(yè)巖氣開(kāi)發(fā)時(shí)，從地層到井底，氣體的流動(dòng)壓力逐漸降低，“雙滑脫”現(xiàn)象在井底附近影響最大。因此，在室內(nèi)滲透率的實(shí)驗(yàn)測(cè)試和實(shí)際頁(yè)巖氣開(kāi)發(fā)中，應(yīng)當(dāng)考慮“雙滑脫效應(yīng)”對(duì)氣體滲流的影響。

圖15 頁(yè)巖氣在納米孔中的無(wú)因次滲透率與平均壓力倒數(shù)的關(guān)系Fig. 15 The relationship between dimensionless permeability and the reciprocal of the average pressure

為了研究頁(yè)巖氣在多孔介質(zhì)中的滲流，基于掃描電鏡圖片，用MCMC方法構(gòu)建出數(shù)字巖心。數(shù)字巖心尺寸為200×200×200，孔隙度為0.1，上方為入口端，下方為出口端。z方向上采用壓力邊界條件，x方向和y方向采用周期性邊界條件。運(yùn)用格子Boltzmann方法進(jìn)行流動(dòng)模擬，研究滑脫效應(yīng)對(duì)視滲透率的影響。當(dāng)流動(dòng)達(dá)到穩(wěn)定時(shí)，運(yùn)用式(5)計(jì)算出數(shù)字巖心視滲透率和固有滲透率。

圖16Kapp/K0和Kapp與有效流動(dòng)半徑的關(guān)系Fig. 16 Variation ofKappandKapp/K0with mean pore radius

圖16表征了視滲透率(Kapp)和視滲透率同固有滲透率之比(Kapp/K0)與有效流動(dòng)半徑(r)之間的關(guān)系。模擬過(guò)程中，多孔介質(zhì)平均壓力分別為5，10，15 MPa，氣體動(dòng)力黏度為1.447×105Pa·s，有效流動(dòng)半徑變化區(qū)間為1～1000 nm。結(jié)果表明，當(dāng)壓力一定時(shí)，隨著r的增大，Kapp/K0不斷減小，Kapp不斷增大。這是因?yàn)殡S著r的增大，氣體流動(dòng)逐漸變?yōu)檫B續(xù)流，滑脫效應(yīng)逐漸減弱，視滲透率無(wú)限趨近于固有滲透率，Kapp/K0逐漸趨近于1。此外，隨著r的增大Kapp的增幅逐漸減小。這是因?yàn)閞較小時(shí)，Knudsen數(shù)較大，氣體分子自由程與r處于相同數(shù)量級(jí)，r微小的變化，就能顯著影響滑脫效應(yīng)強(qiáng)度，導(dǎo)致Kapp較大的增幅。當(dāng)r一定時(shí)，視滲透率隨著壓力的減小不斷增大。但是當(dāng)r足夠大時(shí)，不同壓力下的Kapp近似相等，這是因?yàn)闅怏w流動(dòng)空間足夠大時(shí)，氣體流動(dòng)變?yōu)檫B續(xù)流，視滲透率近似等于固有滲透率。

為了進(jìn)一步表征流動(dòng)空間平均壓力與視滲透率之間的關(guān)系，分別模擬了有效流動(dòng)半徑為5 nm，50 nm，500 nm時(shí)，不同平均壓力下的視滲透率，壓力變化區(qū)間為0.1～50 MPa，如圖17所示。由模擬結(jié)果可知，當(dāng)有效流動(dòng)半徑一定時(shí)，固有滲透率保持不變，平均壓力越大，視滲透率越小。低壓下的視滲透率遠(yuǎn)大于固有滲透率，并且隨著壓力的增大，視滲透率逐漸接近固有滲透率。當(dāng)P=0.1 MPa，小孔徑(5 nm)下氣體視滲透率為6.13×10-3mD，固有滲透率為6.13×10-4mD，視滲透率是固有滲透率的10倍；大孔徑(500 nm)下氣體視滲透率為1.25 mD，固有滲透率為0.61 mD，視滲透率是固有滲透率的2倍。由此可知，隨著有效流動(dòng)空間的增大，視滲透率與固有滲透率之間的差異在不斷縮小。

圖17Kapp和K0與平均壓力的關(guān)系Fig. 17 Variation ofKappandK0with average pressure

5 結(jié)論

(1)氣體在頁(yè)巖納微孔道中的流動(dòng)屬于滑移區(qū)和弱過(guò)渡區(qū)范圍，采用格子Boltzmann方法可以有效模擬頁(yè)巖氣的微觀流動(dòng)機(jī)理。

(2)頁(yè)巖氣在納微孔道中的平均流速大于Poiseuille公式的解析解。壓力越低，納微孔道直徑越小，氣體流動(dòng)速度與Poiseuille公式解析解的差別越大。

(3)隨著頁(yè)巖氣體流動(dòng)通道直徑的降低，孔壁處的氣體滑脫速度逐漸增大。孔徑小于20 nm時(shí)，流體粒子的平均自由程遠(yuǎn)大于孔隙直徑，流體粒子與孔壁的相互作用占主導(dǎo)地位，邊界滑脫速度顯著高于孔道中心的氣體流速，孔道出口端的流速剖面也由Poiseuille流的經(jīng)典拋物線形狀轉(zhuǎn)變?yōu)椤皟啥烁咧虚g平緩”的柱塞狀。

(4)由于“雙滑脫效應(yīng)”的影響，頁(yè)巖氣在納微孔道中的滲透率大于等效液體滲透率；并且孔徑越小、壓力越低，這種增強(qiáng)作用越顯著。同時(shí)，在克氏滲透率校正圖中，氣體滲透率與平均壓力的倒數(shù)偏離了直線關(guān)系，發(fā)生上翹。

(5)壓力一定時(shí)，隨著流動(dòng)空間的增大，視滲透率與固有滲透率的比值不斷減小并逐漸趨于1。當(dāng)流動(dòng)空間一定時(shí)，隨著壓力的增大，視滲透率不斷減小并無(wú)限趨近于固有滲透率。

[1] 鄒才能, 楊智, 陶士振, 等. 納米油氣與源儲(chǔ)共生型油氣聚集[J]. 石油勘探與開(kāi)發(fā), 2012, 39(1): 13-26.[ZOU C N, YANG Z, TAO S Z, et al. Nano-hydrocarbon and the accumulation in coexisting source and reservoir[J]. Petroleum Exploration and Development, 2012, 39(1): 13-26.]

[2] 梁超, 姜在興, 楊鐿婷, 等. 四川盆地五峰組—龍馬溪組頁(yè)巖巖相及儲(chǔ)集空間特征[J]. 石油勘探與開(kāi)發(fā), 2012, 39(6): 691-697. [LIANG C, JIANG Z X, YANG Y T, et al. Characteristics of shale lithofacies and reservoir space of the Wufeng-Longmaxi Formation, Sichuan Basin[J]. Petroleum Exploration and Development, 2012, 39(6): 691-697.]

[3] LOUCKS R G, REED R M, RUPPEL S C, et al. Morphology, genesis, and distribution of nanometer-scale pores in siliceous mudstones of the Mississippian Barnett shale[J]. Journal of Sedimentary Research, 2009, 79(12): 848-861.

[4] SONDERGELD C H, AMBROSE R J, RAI C S, et al. Micro-structural studies of gas shales[R]. SPE 131771, 2010.

[5] 楊峰, 寧正福, 胡昌蓬, 等. 頁(yè)巖儲(chǔ)層微觀孔隙結(jié)構(gòu)特征[J]. 石油學(xué)報(bào), 2013, 34(2): 301-311.[YANG F, NING Z F, HU C P, et al. Characterization of microscopic pore structures in shale reservoirs[J]. Acta Petrolei Sinica, 2013, 34(2): 301-311.]

[6] 楊峰, 寧正福, 張世棟, 等. 基于氮?dú)馕綄?shí)驗(yàn)的頁(yè)巖孔隙結(jié)構(gòu)表征[J]. 天然氣工業(yè), 2013, 33(4): 135-140.[YANG F, NING Z F, ZHANG S D, et al. Characterization of pore structures in shales through nitrogen adsorption experiment[J]. Natural Gas Industry, 2013, 33(4): 135-140.]

[7] 陶然, 權(quán)曉波, 徐建中. 微尺度流動(dòng)研究中的幾個(gè)問(wèn)題[J]. 工程熱物理學(xué)報(bào), 2001, 22(5): 576-577.[TAO R, QUAN X B, XU J Z. Several questions in research of micro scale fl ow[J]. Journal of Engineering Thermophysics, 2001, 22(5): 576-577.]

[8] JAVADPOUR F, FISHER D, UNSWORTH M. Nanoscale gas flow in shale gas sediments and siltstones[J]. Journal of Canadian Petroleum Technology, 2007, 46(10): 55-61.

[9] JAVADPOUR F. Nanopores and apparent permeability of gas flow in mudstones (shales and siltstones)[J]. Journal of Canadian Petroleum Technology, 2009, 48(8): 16-21.

[10] FREEMAN C.M. A numerical study of microscale fl ow behave in tight gas and shale gas reservoir systems[R]. SPE 141125, 2010.

[11] 王瑞, 張寧生，劉曉娟, 等. 頁(yè)巖氣擴(kuò)散系數(shù)和視滲透率的計(jì)算與分析[J]. 西北大學(xué)學(xué)報(bào)(自然科學(xué)版), 2013, 43(1): 75-80. [WANG R, ZHANG N S, LIU X J, et al. The calculation and analysis of diffusion coef fi cient and apparent permeability of shale gas[J]. Journal of Northwest University (Natural Science Edition), 2013, 43(1): 75-80.]

[12] 何雅玲, 王勇, 李慶. 格子Boltzmann方法的理論及應(yīng)用[M]. 北京: 科學(xué)出版社, 2008.[HE Y L, WANG Y, LI Q. Lattice Boltzmann Method: Theory and applications[M]. Beijing: Science Press, 2008.]

[13] 郭照立, 鄭楚光. 格子Boltzmann方法的原理及應(yīng)用[M]. 北京：科學(xué)出版社, 2008. [GUO Z L, ZHENG C G. Theory and applications of Lattice Boltzmann Method[M]. Beijing: Science Press, 2008.]

[14] 曾彥,寧正福,王慶,等. 考慮微尺度效應(yīng)的頁(yè)巖氣藏壓裂水平井壓力動(dòng)態(tài)分析[J]. 中國(guó)科技論文, 2016, (23): 2679-2687.[ZENG Y, NING Z F, WANG Q, et al. Transient fl ow analysis of multiple fractured horizontal wells in shal gas reservoirs with consideration of microscale effect[J]. China Sciencepaper, 2016, (23): 2679-2687.]

[15] MASON E A, MALINAUSKAS A. Gas transport in porous media: The dusty-gas model[M]. Elsevier Science Ltd, 1983.

[16] GUO Z, ZHAO T S. Lattice Boltzmann model for incompressible fl ows through porous media[J]. Physical Review E, 2002, 66(3): 036304.

[17] ZHANG Y H, QIN R S, DAVID R E. Lattice Boltzmann simulation of rare fi ed gas fl ow in micro channels[J], Physical Review E, 2005, 71(4): 1-4.

[18] DONGARI N, AGRAWAL A, AGRAWAL A. Analytical solution of gaseous slip fl ow in long microchannels[J]. International Journal of Heat and Mass Transfer, 2007, 50(17-18): 3411-3421.

[19] LIM C Y, SHU C, NIU X D, et al. Application of Lattice Boltzmann Method to simulation microchannel fl ows[J]. Physics of Fluids, 2002, 14(7): 2299-2308.

[20] TANG G H, TAO W Q, HE Y L. Lattice Boltzmann Method for simulating gas flow in micro-channels[J]. International Journal of Modern Physics C, 2001, 15(2): 335-347.

[21] VERHAEGHE F, LUO L S, BLANPAIN B. Lattice Boltzmann modeling of microchannel flow in slip flow regime[J]. Journal of Computational Physics, 2009, 228(1): 147-157.

[22] ARKILIC E B. Measurement of the mass flow and tangential momentum accommodation coefficient in silicon micromachined channels[D]. Cambridge: Massachusetts Institue of Technology, 1997.

[23] ARKILIC E B, SCHMIDT M A. BREUER K S. Gaseous slip flow in long microchannels[J]. Journal of Microelectromechanical Systems, 1997, 6(2): 167-178.

[24] SHEN C, TIAN D B, XIE C, et al. Examination of the LBM in simulation of microchannel fl ow in transitional regime[J]. Microscale Thermophysical Engineering, 2004, 8(4): 423-432.

[25] 陳代珣. 滲流氣體滑脫現(xiàn)象與滲透率變化的關(guān)系[J]. 力學(xué)學(xué)報(bào), 2002, 34(1): 96-100. [CHEN D X. Gas slippage phenomenon and change of permeability when gas fl ows in tight porous media[J]. Acta Mechanica Sinica, 2002, 34(1): 96-100.]

[26] AKKUTLU I Y, FATHI E. Multiscale gas transport in shales with local kerogen heterogeneities[J]. SPE J, 2012, 17(4): 1001-1011.

[27] SWAMI V, CLARKSON C R, SETTARI A. Non-darcy flow in shale nanopores: Do we have a final answer? [R]. SPE Canadian Unconventional Resources Conference. Society of Petroleum Engineers, 2012.

[28] FATHI E, TINNI A, AKKUTLU I Y. Correction to Klinkenberg slip theory for gas fl ow in nano-capillaries[J]. International Journal of Coal Geology, 2012, 103(1): 51-59.

[29] KLINKENBERG L J. The permeability of porous media to liquid and gases[C]. New York: API Drilling and Production Practices: 200-231.

[30] HELLER R, ZOBACK M. Laboratory measurements of matrix permeability and slippage enhanced permeability in gas shales[C]. Unconventional Resources Technology Conference (URTEC), 2013.

[31] WANG J, CHEN L, KANG Q, et al. Apparent permeability prediction of organic shale with generalized Lattice Boltzmann model considering surface diffusion effect[J]. Fuel, 2016, 181: 478-490.

Simulation of transport of shale gas through the nanopores of shales

ZENG Yan1,2, NING Zhengfu1,2, QI Rongrong1,2, HUANG Liang1,2, LV Chaohui1,2
1 MOE Key Laboratory of Petroleum Engineering, China University of Petroleum-Beijing, Beijing 102249, China
2 State Key Laboratory of Petroleum Resources and Engineering, China University of Petroleum-Beijing, Beijing 102249, China

Abundant nanopores were observed in gas shales using field emission scanning electron microscopy. The Navier-Stokes equation based on the continuum assumption breaks down at the nanometer scale. The transport mechanisms of gas in nanopores of shales have been simulated using the Lattice Boltzmann Method. The simulation result shows that the microscale phenomenon of gas in nanopores of shales leads to bigger average fl ow rates than that from the classic Poiseuille equation. The rarefaction effect is very remarkable, and the slip velocity on the wall boundary is not equal to zero. The smaller the pore diameter is, the larger the slip velocity is, so much so that the slip velocity becomes bigger than gas velocity in the nanopores, and the parabola fl ow velocity pro fi le transforms from parabola to plunger. The “double skip effect” which results from the rebounding gas molecules with kinetic energy entering into the nanopores enhances the fl ow of gas in the whole nanopores, which bring about increased actual permeability, and the permeability of nanopores is bigger than that from liquid permeability. In addition, the “double skip effect” leads to nonlinear deviation from the classic Klinkenberg slip theory, especially in low fl ow pressure and smaller pore diameter.

shale; nanopore; Lattice Boltzmann Method; fl ow simulation; nanoscale; double skip effect

10.3969/j.issn.2096-1693.2017.01.007

(編輯 馬桂霞)

*通信作者, nzf@cup.edu.cn

2016-08-19

國(guó)家自然科學(xué)基金“頁(yè)巖氣儲(chǔ)層吸附解吸機(jī)理研究”(51274214)、教育部科學(xué)技術(shù)研究重大計(jì)劃“頁(yè)巖氣流動(dòng)機(jī)理與產(chǎn)能預(yù)測(cè)模型研究”(311008)、油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室自主研究課題“頁(yè)巖氣吸附與解吸機(jī)理研究”(PRP/indep-3-1108)聯(lián)合資助。

曾彥, 寧正福, 齊榮榮, 黃亮, 呂朝輝. 頁(yè)巖氣在納微孔道中的流動(dòng)模擬研究. 石油科學(xué)通報(bào), 2017, 01: 64-75

ZENG Yan, NING Zhengfu, QI Rongrong, HUANG Liang, LV Chaohui. Simulation of transport of shale gas through the nanopores of shales. Petroleum Science Bulletin, 2017, 01: 64-75. doi: 10.3969/j.issn.2096-1693.2017.01.007