頁巖納米級孔隙甲烷滲透特性模擬

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

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

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

本工作通過四參數隨機增長法構造二維多孔介質模型，采用格子Boltzmann方法模擬甲烷在多孔介質中的流動，研究了多孔介質孔隙率（ε）、進出口壓差和吸附解吸對于甲烷氣體滲透率（K）的影響。

1 模型的建立及邊界條件

1.1 基本模型

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

圖1 D2Q9模型Fig.1 D2Q9 model.

通過BGK近似、時空離散和速度離散，得到完全離散化的格子Boltzmann-BGK方程[4]，見式（1）。

局部平衡態分布函數見式（2）。

模型的速度配置見式（3）。

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

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

1.2 松弛時間與Kn

Kn為頁巖中的氣體流動特征參數，需要通過Kn確定松弛時間[5]，見式（7）。

式中，L為計算域的特征長度。

Kn是模擬的關鍵數據[6]，由式（8）計算得出。

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

ψ（Kn）為矯正函數，由Guo等[7]提出，見式（9）。

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

1.3 邊界條件

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

本工作采用假定速度（us）來表示微觀尺度下頁巖基質中氣流的滑移效應[10]，見式（13）。

α為吸附氣體分子占比，見式（14）。

氣體擴散時，吸附氣在壁面發生表面擴散，自由氣在空隙中產生分子擴散，在此模型中，用Fick定律計算吸附氣和自由氣的流量，見式（15）。

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

β為Langmuir分配系數，由式（17）計算。

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

式中，fi'為標準反彈分布函數；fi''為鏡面反彈分布函數。

反彈系數（Rb）[5]，由式（19）計算。

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

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

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

壁面處密度分布函數見式（21）。

式中，i=0～8。

而進出口邊界采用Zou-He邊界模型[12]，流體在通道內靠兩端壓差驅動，孔道入口端和出口端的分布函數見式（22）～式（27）。

入口端：

出口端：

2 模擬結果及討論

2.1 四參數隨機增長法構造多孔介質

采用Wang等[13]提出的隨機生長四參數生成法（QSGS）構造多孔介質微觀結構。該方法可通過參數調整控制生長介質的形貌特征。構造二維多孔介質時，令固體顆粒為生長相，孔隙為非生長相。實驗過程如下：1）在構造區域內的網格上以中心核生長概率（cdd）隨機配置生長相的生長核，cdd要小于ε；2）按照不同方向上給定的生長概率（dj）（j=1，2，…，8，j為生長方向），固相生長核向周圍鄰點生長；3）重復步驟2），直到非生長相達到給定的ε。

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

2.2 頁巖氣微尺度流動模擬及影響因素

氣體在納米孔隙中的流動受很多因素影響，不同頁巖氣藏的ε、進出口壓差和吸附作用差異較大，這些差異造成了不同頁巖氣藏中氣體滲流規律的不同，理清這些因素對微尺度氣體流動的影響機制具有重要意義。K是表征多孔介質流的重要參數及衡量標準，根據達西定律得到K的計算公式，見式（28）。

2.2.1ε對頁巖氣流動的影響

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

通過四參數隨機增長法構造了不同ε的多孔介質模型，由于不同ε下流通通道的不同，導致了Kn及Rb等的差異，ε對頁巖氣流動有很大的影響。圖4為不同ε下的K。由圖4可知，隨著ε的降低，K也隨之降低。ε高時，孔道連通效果好，孔隙半徑大，有利于流體流動，因此流速會高，K也大；而ε低時，孔道情況會變得復雜，孔隙半徑變小，平均流速會減小，導致K減小。

2.2.2 壓差對頁巖氣流動的影響

頁巖氣開采過程中的驅動力主要是地層的壓力，而地層環境復雜，導致壓差分布不均，因此需考慮壓力對頁巖氣流動的影響。本工作采用壓力驅動模型，進出口為Zou-He壓力邊界，出入口壓力由密度確定，當ε=0.7時，模擬頁巖氣在不同壓差下的K，結果見圖5。

圖2 多孔介質模型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可知，隨著進出口壓差的增大，K也增大。進出口壓差增大，流體流動的推動力增大，流體流速提高，K會升高，但在進出口壓差變得很大時K才會有顯著的提高。

2.2.3 吸附/解吸、擴散對頁巖氣流動的影響

頁巖氣開采過程中吸附/解吸、擴散等微尺度效應對油氣運移具有一定影響，在模擬時通過邊界條件體現，需綜合考慮PL、Ds、吸附氣密度等諸多因素的影響。分別模擬了吸附氣密度確定時PL與Ds變化及PL確定時Ds與吸附氣密度變化對K的影響，考察微納米孔隙內甲烷氣體吸附/解吸、表面擴散對流動的影響，結果見圖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可知，在吸附氣濃度不變時，隨著Ds與分子擴散系數（Dk）比值的提高，K會下降，而PL提高時，K也降低；在PL不變的情況下，吸附氣濃度提高時，K下降，但兩種情況下K下降幅度都不大。在頁巖氣流動時，PL和Ds只影響頁巖氣的壁面速度，而壁面速度的微小改變不會對整體K造成太大的影響。

3 結論

1）頁巖氣開采過程為微納米孔徑內的流動，使用格子Boltzmann方法可很好地模擬甲烷在多孔介質中的流動。

2）模擬頁巖氣在多孔介質中的流動時需考慮Kn及吸附/解吸的影響。

3）多孔介質ε對K的影響最為明顯，ε越高越有利于甲烷流動；壓差也會對K產生影響，壓差大時驅動力大，有利于甲烷流動。甲烷在頁巖表面的吸附/解吸會對壁面速度產生一定影響，但對整體K的影響不大。

符 號 說 明