頁巖氣在納微孔道中的流動模擬研究

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

1 中國石油大學(北京)石油工程教育部重點實驗室, 北京 102249

2 中國石油大學(北京)油氣資源與工程國家重點實驗室, 北京 102249

頁巖氣在納微孔道中的流動模擬研究

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

1 中國石油大學(北京)石油工程教育部重點實驗室, 北京 102249

2 中國石油大學(北京)油氣資源與工程國家重點實驗室, 北京 102249

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

頁巖；納米孔；格子Boltzmann方法；流動模擬；納米尺度；雙滑脫效應

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

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

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

1.1 頁巖中的納米孔隙

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

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

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

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

1.2 流動狀態(tài)

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

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

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

2 格子Boltzmann方法

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

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

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

2.1 計算模型

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

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

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

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

式中，Lx，Ly，Lz分別為多孔介質(zhì)x，y，z三個方向的長度，其中x方向為流動方向。ρin，ρout分別為進出口端密度。引入視滲透率Kapp，則有：

2.2 格子Boltzmann模型

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

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

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

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

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

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

模型的速度配置如下：

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

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

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

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

2.3 邊界條件

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

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

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

(3)語法復習目標：掌握各項基礎語法的基本概念及重要用法；熟練理解、掌握并能夠在不同的語境中辨別、正確、靈活運用各種重要語法功能。

3 模型驗證

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

3.1 連續(xù)介質(zhì)區(qū)的流體流動驗證

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

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

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

3.2 滑脫流和弱過渡流的流體流動驗證

對于連續(xù)性假設失效或宏觀方程難以描述的現(xiàn)象，采用微觀分子動力學或介觀氣體動力學描述更加恰當。為了和他人的結(jié)果[21-24]進行對比，計算了孔道出口端的流速分布剖面。模型的基本參數(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)為弱過渡流階段的模擬結(jié)果，其中出口端Knudsen數(shù)Kn0= 0. 194，入口端密度ρin=2。

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

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

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

4.1 平均流速

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

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

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

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

4.2 空間流速分布

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

圖9 頁巖氣在1 μm孔中流動的空間速度分布Fig. 9 Dimensionless velocity distribution in 1 μm pore

圖10 頁巖氣在50 nm孔中流動的空間速度分布Fig. 10 Dimensionless velocity distribution in 50 nm pore

4.3 流速剖面

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

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

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

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

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

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

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

4.4 孔道滲透率

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

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

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

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

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

圖16Kapp/K0和Kapp與有效流動半徑的關系Fig. 16 Variation ofKappandKapp/K0with mean pore radius

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

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

圖17Kapp和K0與平均壓力的關系Fig. 17 Variation ofKappandK0with average pressure

5 結(jié)論

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

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

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

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

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

[1] 鄒才能, 楊智, 陶士振, 等. 納米油氣與源儲共生型油氣聚集[J]. 石油勘探與開發(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] 梁超, 姜在興, 楊鐿婷, 等. 四川盆地五峰組—龍馬溪組頁巖巖相及儲集空間特征[J]. 石油勘探與開發(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] 楊峰, 寧正福, 胡昌蓬, 等. 頁巖儲層微觀孔隙結(jié)構(gòu)特征[J]. 石油學報, 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] 楊峰, 寧正福, 張世棟, 等. 基于氮氣吸附實驗的頁巖孔隙結(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] 陶然, 權曉波, 徐建中. 微尺度流動研究中的幾個問題[J]. 工程熱物理學報, 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] 王瑞, 張寧生，劉曉娟, 等. 頁巖氣擴散系數(shù)和視滲透率的計算與分析[J]. 西北大學學報(自然科學版), 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方法的理論及應用[M]. 北京: 科學出版社, 2008.[HE Y L, WANG Y, LI Q. Lattice Boltzmann Method: Theory and applications[M]. Beijing: Science Press, 2008.]

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

[14] 曾彥,寧正福,王慶,等. 考慮微尺度效應的頁巖氣藏壓裂水平井壓力動態(tài)分析[J]. 中國科技論文, 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)象與滲透率變化的關系[J]. 力學學報, 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

國家自然科學基金“頁巖氣儲層吸附解吸機理研究”(51274214)、教育部科學技術研究重大計劃“頁巖氣流動機理與產(chǎn)能預測模型研究”(311008)、油氣資源與探測國家重點實驗室自主研究課題“頁巖氣吸附與解吸機理研究”(PRP/indep-3-1108)聯(lián)合資助。

曾彥, 寧正福, 齊榮榮, 黃亮, 呂朝輝. 頁巖氣在納微孔道中的流動模擬研究. 石油科學通報, 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