頁(yè)巖納米級(jí)孔隙甲烷滲透特性模擬

張朋旗，吳家全，王桂珠，張永康，郭麗梅

（天津科技大學(xué) 化工與材料學(xué)院，天津 300457）

隨著常規(guī)油氣資源的減少，頁(yè)巖氣等非常規(guī)資源的勘探和開發(fā)成為了新的研究熱點(diǎn)。頁(yè)巖氣藏主要由基質(zhì)（有機(jī)質(zhì)、無(wú)機(jī)質(zhì)） 和天然裂縫組成，且絕大部分有機(jī)質(zhì)和無(wú)機(jī)質(zhì)中的孔隙均為納米級(jí)孔隙[1]。而氣體在納米級(jí)尺度中流動(dòng)時(shí)，往往會(huì)出現(xiàn)微尺度效應(yīng)。Knudsen數(shù)（Kn）是微尺度流動(dòng)的特征參數(shù)，根據(jù)Kn可將流動(dòng)劃分為4個(gè)區(qū)域[2]：連續(xù)介質(zhì)區(qū)Kn≤0.001，流體可視為連續(xù)介質(zhì)，采用無(wú)滑移的Navier-Stokes（N-S）方程描述；滑移區(qū)0.001≤Kn≤0.1，流動(dòng)存在微尺度效應(yīng)，仍可視為連續(xù)介質(zhì)，采用滑移邊界的N-S方程描述；過(guò)渡區(qū)0.1≤Kn≤10，連續(xù)介質(zhì)假設(shè)不再成立，可采用Burnett方程描述；自由分子流區(qū)Kn≥10，主要采用分子動(dòng)力學(xué)方法進(jìn)行研究。對(duì)于微尺度氣體流動(dòng)的研究，目前常用的方法有分子動(dòng)力學(xué)方法、直接模擬蒙特卡洛方法和格子Boltzmann方法等[3]。由于格子Boltzmann方法具有編程簡(jiǎn)單，易于并行計(jì)算，且能夠處理復(fù)雜邊界條件的優(yōu)點(diǎn)，近幾年在頁(yè)巖氣數(shù)值模擬研究領(lǐng)域得到廣泛的應(yīng)用。

本工作通過(guò)四參數(shù)隨機(jī)增長(zhǎng)法構(gòu)造二維多孔介質(zhì)模型，采用格子Boltzmann方法模擬甲烷在多孔介質(zhì)中的流動(dòng)，研究了多孔介質(zhì)孔隙率（ε）、進(jìn)出口壓差和吸附解吸對(duì)于甲烷氣體滲透率（K）的影響。

1 模型的建立及邊界條件

1.1 基本模型

采用D2Q9格子Boltzmann-BGK[4]模型，如圖1所示。

圖1 D2Q9模型Fig.1 D2Q9 model.

通過(guò)BGK近似、時(shí)空離散和速度離散，得到完全離散化的格子Boltzmann-BGK方程[4]，見(jiàn)式（1）。

局部平衡態(tài)分布函數(shù)見(jiàn)式（2）。

模型的速度配置見(jiàn)式（3）。

式中，c為格子速度，c=δx/δt。

模型的宏觀密度、壓力和速度定義見(jiàn)式（4）～式（6）。

1.2 松弛時(shí)間與Kn

Kn為頁(yè)巖中的氣體流動(dòng)特征參數(shù)，需要通過(guò)Kn確定松弛時(shí)間[5]，見(jiàn)式（7）。

式中，L為計(jì)算域的特征長(zhǎng)度。

Kn是模擬的關(guān)鍵數(shù)據(jù)[6]，由式（8）計(jì)算得出。

式中，d為分子的直徑，m；m為分子的質(zhì)量；H表示計(jì)算域高度；χ為矯正因子。

ψ（Kn）為矯正函數(shù)，由Guo等[7]提出，見(jiàn)式（9）。

Yao等[8]考慮到氣體致密性的影響，將松弛時(shí)間進(jìn)行更進(jìn)一步修正，見(jiàn)式（10）～式（12）。

1.3 邊界條件

邊界條件在頁(yè)巖氣模擬中起著重要作用，頁(yè)巖氣在通道中流動(dòng)時(shí)，必須考慮在壁面上的吸附/解吸及擴(kuò)散等現(xiàn)象，通過(guò)格子Boltzmann方程邊界條件體現(xiàn)微尺度效應(yīng)的影響。在一些研究中，在格子Boltzmann模型的邊界條件中引入Langmuir吸附等溫線，以呈現(xiàn)頁(yè)巖中的吸附/解吸情況[9]。

本工作采用假定速度（us）來(lái)表示微觀尺度下頁(yè)巖基質(zhì)中氣流的滑移效應(yīng)[10]，見(jiàn)式（13）。

α為吸附氣體分子占比，見(jiàn)式（14）。

氣體擴(kuò)散時(shí)，吸附氣在壁面發(fā)生表面擴(kuò)散，自由氣在空隙中產(chǎn)生分子擴(kuò)散，在此模型中，用Fick定律計(jì)算吸附氣和自由氣的流量，見(jiàn)式（15）。

吸附氣引起的壁面速度[11]，見(jiàn)式（16）。

β為L(zhǎng)angmuir分配系數(shù)，由式（17）計(jì)算。

考慮氣體在壁面處的滑脫情況，引入混合反彈邊界條件[12]，見(jiàn)式（18）。

式中，fi'為標(biāo)準(zhǔn)反彈分布函數(shù)；fi''為鏡面反彈分布函數(shù)。

反彈系數(shù)（Rb）[5]，由式（19）計(jì)算。

式中，σ=1，B1=1-0.181 7σ，B=0.8，σv=(2-σ)/σ。

將Langmuir等溫吸附方程引入本模型的邊界條件中，需對(duì)平衡態(tài)分布函數(shù)以及非平衡態(tài)分布函數(shù)進(jìn)行組合，得到新的分布函數(shù)，根據(jù)Chapman-Enskog方法，定義壁面處的分布函數(shù)[12]，見(jiàn)式（20）。

式中，η=1/（τ+0.5）。

壁面處密度分布函數(shù)見(jiàn)式（21）。

式中，i=0～8。

而進(jìn)出口邊界采用Zou-He邊界模型[12]，流體在通道內(nèi)靠?jī)啥藟翰铗?qū)動(dòng)，孔道入口端和出口端的分布函數(shù)見(jiàn)式（22）～式（27）。

入口端：

出口端：

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

2.1 四參數(shù)隨機(jī)增長(zhǎng)法構(gòu)造多孔介質(zhì)

采用Wang等[13]提出的隨機(jī)生長(zhǎng)四參數(shù)生成法（QSGS）構(gòu)造多孔介質(zhì)微觀結(jié)構(gòu)。該方法可通過(guò)參數(shù)調(diào)整控制生長(zhǎng)介質(zhì)的形貌特征。構(gòu)造二維多孔介質(zhì)時(shí)，令固體顆粒為生長(zhǎng)相，孔隙為非生長(zhǎng)相。實(shí)驗(yàn)過(guò)程如下：1）在構(gòu)造區(qū)域內(nèi)的網(wǎng)格上以中心核生長(zhǎng)概率（cdd）隨機(jī)配置生長(zhǎng)相的生長(zhǎng)核，cdd要小于ε；2）按照不同方向上給定的生長(zhǎng)概率（dj）（j=1，2，…，8，j為生長(zhǎng)方向），固相生長(zhǎng)核向周圍鄰點(diǎn)生長(zhǎng)；3）重復(fù)步驟2），直到非生長(zhǎng)相達(dá)到給定的ε。

模擬儲(chǔ)層區(qū)域壓力約30 MPa，溫度約80 ℃，頁(yè)巖孔徑分布范圍在5～150 nm之間[14]，油層氣主要成分為甲烷，在此條件下的密度為166.782 2 kg/m3，聲速為622.823 5 m/s，黏度為1.298 43×10-7m2/s。網(wǎng)格步長(zhǎng)（δx）在實(shí)際空間的長(zhǎng)度為0.72 nm，設(shè)置實(shí)際模擬區(qū)域長(zhǎng)約200 nm，寬約100 nm，建立幾種不同ε的多孔介質(zhì)模型，見(jiàn)圖2。

2.2 頁(yè)巖氣微尺度流動(dòng)模擬及影響因素

氣體在納米孔隙中的流動(dòng)受很多因素影響，不同頁(yè)巖氣藏的ε、進(jìn)出口壓差和吸附作用差異較大，這些差異造成了不同頁(yè)巖氣藏中氣體滲流規(guī)律的不同，理清這些因素對(duì)微尺度氣體流動(dòng)的影響機(jī)制具有重要意義。K是表征多孔介質(zhì)流的重要參數(shù)及衡量標(biāo)準(zhǔn)，根據(jù)達(dá)西定律得到K的計(jì)算公式，見(jiàn)式（28）。

2.2.1ε對(duì)頁(yè)巖氣流動(dòng)的影響

采用甲烷為模擬頁(yè)巖氣，將模擬進(jìn)出口壓差設(shè)為 0.4 MPa，Langmuir壓力（PL）為 10 MPa，Langmuir體積為1.7 m3/t，表面擴(kuò)散系數(shù)（Ds）為1×10-6m3/s，吸附氣密度為游離氣密度的2倍。圖3為不同ε下流體的流速分布。

通過(guò)四參數(shù)隨機(jī)增長(zhǎng)法構(gòu)造了不同ε的多孔介質(zhì)模型，由于不同ε下流通通道的不同，導(dǎo)致了Kn及Rb等的差異，ε對(duì)頁(yè)巖氣流動(dòng)有很大的影響。圖4為不同ε下的K。由圖4可知，隨著ε的降低，K也隨之降低。ε高時(shí)，孔道連通效果好，孔隙半徑大，有利于流體流動(dòng)，因此流速會(huì)高，K也大；而ε低時(shí)，孔道情況會(huì)變得復(fù)雜，孔隙半徑變小，平均流速會(huì)減小，導(dǎo)致K減小。

2.2.2 壓差對(duì)頁(yè)巖氣流動(dòng)的影響

頁(yè)巖氣開采過(guò)程中的驅(qū)動(dòng)力主要是地層的壓力，而地層環(huán)境復(fù)雜，導(dǎo)致壓差分布不均，因此需考慮壓力對(duì)頁(yè)巖氣流動(dòng)的影響。本工作采用壓力驅(qū)動(dòng)模型，進(jìn)出口為Zou-He壓力邊界，出入口壓力由密度確定，當(dāng)ε=0.7時(shí)，模擬頁(yè)巖氣在不同壓差下的K，結(jié)果見(jiàn)圖5。

圖2 多孔介質(zhì)模型Fig.2 Porous media model.Porosity(ε)：a 0.9；b 0.8；c 0.7；d 0.6

圖3 流速分布Fig.3 Flow velocity distribution.ε：a 0.9；b 0.8；c 0.7；d 0.6

由圖5可知，隨著進(jìn)出口壓差的增大，K也增大。進(jìn)出口壓差增大，流體流動(dòng)的推動(dòng)力增大，流體流速提高，K會(huì)升高，但在進(jìn)出口壓差變得很大時(shí)K才會(huì)有顯著的提高。

2.2.3 吸附/解吸、擴(kuò)散對(duì)頁(yè)巖氣流動(dòng)的影響

頁(yè)巖氣開采過(guò)程中吸附/解吸、擴(kuò)散等微尺度效應(yīng)對(duì)油氣運(yùn)移具有一定影響，在模擬時(shí)通過(guò)邊界條件體現(xiàn)，需綜合考慮PL、Ds、吸附氣密度等諸多因素的影響。分別模擬了吸附氣密度確定時(shí)PL與Ds變化及PL確定時(shí)Ds與吸附氣密度變化對(duì)K的影響，考察微納米孔隙內(nèi)甲烷氣體吸附/解吸、表面擴(kuò)散對(duì)流動(dòng)的影響，結(jié)果見(jiàn)圖6。

圖4 不同ε下甲烷的KFig.4 Methane permeability(K) at different ε.

圖5 不同壓差下甲烷的KFig.5 K of methane at different pressure differences.

圖6 不同Ds/Dk下甲烷的KFig.6 K of methane at different Ds/Dk.Ds：surface diffusion coefficient；Dk：molecular diffusion coefficient；Cμ：adsorbed gas concentration；C：free gas concentration；PL：Langmuir pressure.a Cμ/C=5；b PL=10 MPa

由圖6可知，在吸附氣濃度不變時(shí)，隨著Ds與分子擴(kuò)散系數(shù)（Dk）比值的提高，K會(huì)下降，而PL提高時(shí)，K也降低；在PL不變的情況下，吸附氣濃度提高時(shí)，K下降，但兩種情況下K下降幅度都不大。在頁(yè)巖氣流動(dòng)時(shí)，PL和Ds只影響頁(yè)巖氣的壁面速度，而壁面速度的微小改變不會(huì)對(duì)整體K造成太大的影響。

3 結(jié)論

1）頁(yè)巖氣開采過(guò)程為微納米孔徑內(nèi)的流動(dòng)，使用格子Boltzmann方法可很好地模擬甲烷在多孔介質(zhì)中的流動(dòng)。

2）模擬頁(yè)巖氣在多孔介質(zhì)中的流動(dòng)時(shí)需考慮Kn及吸附/解吸的影響。

3）多孔介質(zhì)ε對(duì)K的影響最為明顯，ε越高越有利于甲烷流動(dòng)；壓差也會(huì)對(duì)K產(chǎn)生影響，壓差大時(shí)驅(qū)動(dòng)力大，有利于甲烷流動(dòng)。甲烷在頁(yè)巖表面的吸附/解吸會(huì)對(duì)壁面速度產(chǎn)生一定影響，但對(duì)整體K的影響不大。

符 號(hào) 說(shuō) 明