基于REV 尺度格子Boltzmann 方法的頁巖氣藏滲流規律

符東宇，李勇明，趙金洲，江有適，陳曦宇，許文俊

（油氣藏地質及開發工程國家重點實驗室·西南石油大學，成都 610500）

0 引言

近年來，隨著常規油氣資源儲量的減少以及油氣勘探開發技術的提高，我國正持續加大以頁巖氣、致密氣為代表的非常規油氣資源的勘探開發力度［1-3］。不同于常規油氣藏，頁巖氣藏富含微納米級孔隙，孔滲條件差，基質孔隙度通常＜12%，滲透率通常＜1 mD［4］。目前，頁巖氣藏的商業開發主要依賴于通過大型水力壓裂來盡可能溝通儲層中發育的天然微裂縫。所以，開展氣體在頁巖儲層中的滲流機理研究，對于頁巖氣藏的滲流理論研究及工程實踐都具有重要意義。

由于頁巖氣藏獨特的儲層物性，頁巖氣的運移過程呈現出從分子尺度到宏觀尺度的多尺度的特性。頁巖氣的流動機理包括了自由氣在人工裂縫中的達西流動和在基質微孔隙中的連續流、滑脫流以及過渡流，吸附氣在微孔隙表面吸附/解吸以及擴散［5］。通常對非常規裂縫性氣藏滲流規律的研究手段主要有室內實驗以及數值模擬方法［6-7］。由于非常規氣藏超低孔、低滲的儲層特性，對實驗設備要求較高，流體的滲流過程也難以進行直觀地觀測，所以，通過數值方法模擬流體在多孔介質中的滲流方式得到廣泛的應用。目前已有大量的用于表征氣體在納米級孔隙中運移的模型。此類模型主要分為2 類：一類是基于Hagen-Poiseuille 模型對其邊界條件進行修正，從而表征頁巖氣的多種運移機制［8-10］；另一類是按一定的權重系數對多種運移機制進行耦合并考慮氣體分子間的相互作用建立的真實氣體傳輸模型［11］。此類模型僅能用于描述頁巖氣在單一納米孔隙中的運移規律，尚不能直接用于研究頁巖氣在天然微裂縫發育的復雜多孔介質中的運移特征。

格子Boltzmann 方法（LBM）以其易于處理復雜幾何邊界問題、編程較簡單等優勢在近年來被大量用于非常規氣藏微觀滲流模擬研究［12-14］。LBM模擬分為孔隙尺度以及表征單元（REV）尺度，其中孔隙尺度主要用于模擬流體在微納米尺度孔隙中的流動，而不能直接用于模擬流體在天然裂縫發育的多孔介質中的流動過程［15］。Guo 等［16］通過在廣義格子Boltzmann 方程中添加修正項來實現孔隙尺度向REV 尺度的升級。此后，該方法被大量運用于非常規油氣藏滲流模擬。Zhao 等［17］運用REV尺度的格子Boltzmann 方法，模擬了氣相在重構多孔介質中的滲流過程，并對多孔介質滲透率進行了計算。Zhao 等［15］應用格子Boltzmann 方法對頁巖氣在基質中的滲流過程進行了REV 尺度的模擬，并對頁巖氣滲流規律的影響因素進行了分析。張烈輝等［18］基于頁巖掃描電鏡圖像運用REV 尺度的格子Boltzmann 方法模擬了頁巖氣在單一天然微裂縫中的滲流。Zhang 等［19］基于頁巖掃描電鏡建立的二值圖像模擬了水相在富含有機質頁巖多孔介質中的滲流過程，并通過在二值圖像上規則排布天然裂縫的方式分析了天然裂縫的存在對頁巖氣藏滲流規律的影響。

本文首先重構天然微裂縫發育的頁巖儲層多孔介質模型，然后對Zhao 等［15］以及Liu 等［20］提出的格子Boltzmann 模型圍繞頁巖氣藏特性對演化方程以及邊界條件進行針對性修正，并將其升級到REV 尺度，使其能夠表征頁巖氣在裂縫性多孔介質中滑脫、吸附/解吸以及氣體擴散等運移機制；應用修正模型模擬氣體在天然微裂縫發育的多孔介質中的流動過程；最后基于模擬結果對氣體滲流特性以及滲透率影響因素進行分析，以期為頁巖氣藏的滲流規律研究提供新思路。

1 頁巖氣藏表征單元體尺度微觀滲流模型

1.1 物理模型

本文建立的頁巖氣藏滲流模擬流動域如圖1所示。整個滲流過程考慮二維流動，流場上、下均為不滲透邊界；整個滲流過程為從入口處通過壓力驅動從左至右的流動；氣體滲流過程的速度分布用格子Boltzmann 方法中的二維九速模型（D2 Q9）進行表征。由于模型的求解過程復雜，為了提高求解效率，對模型作出了如下假設：①流場中僅考慮甲烷氣體的單相流動；②滲流過程為等溫滲流過程；③對于天然微裂縫的物性參數僅考慮裂縫長度、裂縫寬度以及裂縫密度。

圖1 低滲致密裂縫性氣藏滲流場物理模型及氣體分子D2 Q9 模型速度分布Fig.1 Schematic of 2D physic model for simulating gas transport in naturally micro-fractured gas reservoir and the velocity directions of gas transport in D2 Q9 model

1.2 修正REV 尺度的格子Boltzmann 模型

1.2.1 外力項修正

本文使用LBM 來表征氣體在上述物理模型中的滲流過程。該方法是將Boltzmann 方程進行離散，把流體分子運動分為遷移和碰撞等2 個過程［21］。廣義格子Boltzmann 演化方程為

式中：x為格子空間格點位置；t為格子空間的某一時刻；fi（x，t）為t時刻格子空間中x格點沿i方向的氣體密度分布函數；Δt為格子空間的時間步長；τ為馳豫時間；Fi（x，t）為格子空間x格點沿i方向的外力項函數；ei為格子空間的離散速度。

ei在D2 Q9 模型中氣體分子的運動方向分為以下9 個方向。

式中：i為方向編號。

式中：ωi為加權系數，取值為ω0=4/9，ω1-4=1/9，ω5-8=1/36；cs為格子空間聲速，c一般取值1，故；I為單位矩陣；ρ為格子空間流體密度；u為格子空間流體速度；φ為多孔介質孔隙度，%。

馳豫時間對于LBM 模型至關重要，由于在頁巖氣藏中，氣體的流動特性由努森數（Kn）控制，所以對于馳豫時間必須進行相應的修正，其表達式［13］為

式中：L為流場特征長度，m；Δx格子尺寸，m；φ(Kn)為考慮努森層的修正函數，φ(Kn)=；χ為修正系數。

χ的計算式為

式中：b=2 πd3/(3m)。其中，d為甲烷分子直徑，m；m為甲烷分子質量，kg。

氣體的動力黏度也是影響滲流過程的重要參數，其在格子空間的計算式為

式（1）中的格子空間外力項Fi考慮孔隙表面對于氣體分子的作用，其離散形式為

式中：Φ(x,t)為勢函數；F為格子空間氣體分子所受總外力項。

F的計算式為

式中：G為格子空間氣體分子所受外部力，主要為基質及裂縫表面對氣體分子的吸附力；Fε為幾何形狀因子；K為固有滲透率，m2。

Fε與K的計算式［23］為

式中：dp為平均孔隙直徑，m。

氣體的宏觀密度以及速度可以通過下式進行計算

由于式（8）中也包含流體流速u，所以要對式（8）進行求解，必須聯立式（8）及式（11）

V為臨時變量，其計算式為

聯立式（12）及式（13）則可以得到氣體宏觀速度u的計算式為

d0及d1的完整表達式為

1.2.2 邊界條件修正

本文模型考慮了頁巖氣藏中氣體吸附/解吸及擴散效應，所以必須對傳統的格子Boltzmann 方法的邊界條件進行修正。巖石表面的氣體滑移流動速度的表達式［24］為

式中：ug為離邊界層最近的氣流速度，m/s；uw為假定的固體表面速度，m/s；α為比例系數。

表征的氣體在孔隙及裂縫表面的吸附情況，可由下式進行計算

式中：Cμ及Cμmax分別為吸附氣濃度及最大吸附濃度，g/m3；PL為蘭格繆爾壓力，MPa。

式（16）中的uw為假定的邊界滑移速度，在頁巖氣藏須要考慮氣體的吸附以及表面擴散作用，其表達式為

式中：DK為氣體分子擴散系數，m2/s；DS為氣體表面擴散系數，m2/s。

在格子空間s可由下式進行計算［25］：

2 結果與討論

2.1 模型驗證

Guo 等［23］提出了對REV 尺度的格子Boltzmann模型的驗證方法，后來該方法被大量應用于REV尺度格子Boltzmann 模型準確性的驗證［15，19-20］。本文也選取該方法對所建模型的準確性進行驗證。

首先設置模型上、下不滲透平板間的多孔介質孔隙度為20%，流體在模型入口處受到恒定的壓力驅動，其流動過程滿足Brinkman 達西方程

其解析解為

式中：a=；ν，νe分別是流體的黏度和有效黏度，本文取ν/νe=1，m2/s；L，H分別為流場的長度和寬度，m。

格子空間中的一些必要參數取值為：網格數Lx×Ly=100×100，流體密度ρ=1，黏度ν=1/6，驅動壓差ΔP=1。圖2 為不同達西數下的解析解與本文所建立模型的數值解的對比。從圖2 可看出，本文模型的計算結果與經典解析解吻合，證明了本文所建立的REV 尺度格子Boltzmann 模型（LB 模型）具有較高的準確性。

圖2 不同達西數下不滲透平板間的無因次速度剖面Fig.2 Velocity profiles of fluid flow between two plates under different Darcy number

2.2 算例及參數敏感性分析

2.2.1 基本算例

本文模型模擬所用的基本參數及其格子空間尺度向真實物理空間尺度的轉換關系如表1 所列。在LBM 中主要通過保持格子空間和物理空間流動特征數不變的準則來實現尺度的轉換。將表1 中的參數帶入本文所建立的LB 模型中，然后通過編程求解可以得到氣體在重構頁巖基質多孔介質中的滲流動態。圖3 展示了在不同時間步長下的氣體流速分布情況。從圖3 可以看出，當天然裂縫較發育時，大量氣體沿天然裂縫竄進，天然裂縫溝通程度越顯著時氣流速度越大。以上結論證實了天然裂縫發育及儲層中天然裂縫的溝通作用對于頁巖氣藏開發十分重要。

表1 模型基本參數Table 1 Parameters used for simulation

圖3 不同時間步長下的氣流速度分布情況Fig.3 Gas velocity distribution under different time steps

2.2.2 天然微裂縫物性參數分析

為了進一步明確天然微裂縫對于頁巖氣藏氣體滲流特性的影響，圖4 展示了不同裂縫密度下（ρf=0.02，0.04，0.08）的氣體流速分布情況。在3 種模擬條件下，僅改變裂縫的線密度、裂縫寬度Wf=（2～4）×10-6m，裂縫長度lf=（1～4）×10-4m、裂縫傾角θ=0°～60°，其余流體參數保持表1 中的數值不變。從圖4 可以看出，當裂縫密度增大時，氣體流動通道增加，滲流能力顯著增強。

圖4 不同裂縫密度下氣體流速分布云圖Fig.4 Distributions of gas velocity under different fracture density

圖5為不同天然裂縫密度條件下，所建立裂縫性多孔介質的表觀滲透率與基質固有滲透率比值的變化情況。從圖5 可以看出，隨著天然裂縫密度的增加，表觀滲透率將顯著增大。

同樣僅改變裂縫開度，保持天然裂縫ρf=0.04及lf=（1～4）×10-4m 不變，將裂縫參數以及流體參數帶入LB 模型，計算所得表觀滲透率的增大情況如圖6 所示。從圖6 可以看出，隨著裂縫開度的增加，表觀滲透率增量會隨之增大。這是由于天然裂縫開度增加后，氣體在天然裂縫的流動通道增大，從而引起氣體表觀滲透率增大。

圖5 不同裂縫密度下的表觀滲透率增大倍數Fig.5 Apparent permeability increase multiple under different fracture density

圖6 不同裂縫開度下的表觀滲透率增大倍數Fig.6 Apparent permeability increase multiple under different fracture aperture

僅改變裂縫長度，保持其他裂縫參數為ρf=0.04及Wf=（2～4）×10-6m 不變，同樣將上述裂縫參數及流體物性參數帶入本文LB 模型進行計算，得到表觀滲透率增大倍數情況如圖7 所示。從圖7 可以看出，隨著裂縫長度的增加，表觀滲透率的增大倍數也隨之增加。對比裂縫寬度對于表觀滲透率增大量的貢獻，裂縫長度對于表觀滲透率的增大作用更為顯著。這是由于隨著裂縫長度的增加，基質向裂縫的竄流會加劇，從而引起表觀滲透率的大幅上升。

2.2.3 氣體滲流特征參數分析

圖7 不同裂縫長度下的表觀滲透率增大倍數Fig.7 Apparent permeability increase multiple under different fracture length

在頁巖氣藏中，氣體運移規律不同于常規儲層。為了進一步明確頁巖儲層中氣體多重運移的規律，有必要對氣體滲流特征參數進行敏感性分析。

當儲層致密，氣體滲流通道狹窄時，氣體分子自由程將接近于其流動通道的特征長度，此時Kn對于氣體流動的影響將加劇。氣體在頁巖儲層中的流動可由Kn的大小將其劃分為4 種流態［26］。其中Kn≤0.001 為連續介質流；0.001＜Kn≤0.100 為滑脫流；0.1＜Kn≤10.0 為過渡流；Kn>10 為自由分子流。為了能夠定量分析Kn對于氣體流動特征的影響，筆者模擬了Kn為0.007，0.069，0.687 條件下的氣體滲流過程。

圖8 不同Kn 下的氣體流速分布云圖Fig.8 Distribution of gas velocity under differentKn

圖8 表明，隨著Kn的增大，氣體流動過程中克氏效應越顯著，氣體流速也隨之增大。為了能夠直觀觀測Kn對于氣體滲流速度的影響，圖9 展示了不同Kn下的氣體流速剖面。從圖9 可看出，當Kn持續增大后，裂縫壁面的滑脫效應更加顯著，使得氣體沿天然微裂縫竄進的速度大幅提升；特別是當Kn>0.1 時，流態進入過渡流后，氣體流速增幅顯著變大。圖10 為不同Kn下的裂縫性多孔介質表觀滲透率隨出口壓力的變化情況。從圖10 可以看出，隨著Kn的增大，表觀滲透率也隨之增大；特別是當Kn>0.1 時，表觀滲透率增幅明顯增大。

圖9 不同Kn 下的氣體流速剖面Fig.9 Profiles of gas velocity magnitude under differentKn

圖10 不同Kn 下的表觀滲透率隨出口壓力的變化情況Fig.10 Variance of apparent permeability with the rising outlet pressure under differentKn

頁巖氣藏中氣體的吸附/解吸對于氣體滲流有巨大的影響。為了能夠直觀表征氣體吸附/解吸作用對于頁巖氣藏表觀滲透率的影響，本文計算了不同吸附氣量下多孔介質的表觀滲透率。圖11 為不同吸附氣濃度下表觀滲透率的變化情況。從圖11可以看出，隨著吸附氣量的增加，多孔介質的表觀滲透率隨之下降。出現該結果的主要原因是由于在頁巖有機質中，吸附氣吸附在孔隙表面會占據部分游離氣的流動通道，使基質的有效孔隙減小，從而影響氣體的滲流［27］。

圖11 不同吸附氣濃度下的表觀滲透率變化情況Fig.11 Variance of apparent permeability under different adsorption strength

目前，大量研究表明［28-30］，頁巖氣藏有機質中吸附氣的表面擴散作用是影響頁巖氣滲流過程的關鍵因素。在本文模型中，氣體表面擴散作用的強弱是通過表面擴散系數Ds來進行表征。圖12 為不同表面擴散系數比值下（Ds/Dk=0，0.01，0.1，1），頁巖氣藏多孔介質表觀滲透率隨壓力的變化情況。從圖12 可以看出，在壓力較低的情況下，考慮氣體表面擴散效應的表觀滲透率數值較不考慮表面擴散效應的情況能夠提升2～5 倍。這樣的模擬結果表明氣體表面擴散效應對于表觀滲透率的貢獻作用是不可忽略的。

圖12 不同表面擴散系數下的表觀滲透率變化情況Fig.12 Variance of apparent permeability under different surface diffusivity

在頁巖有機質中，氣體表面擴散效應并不能無限變化，其強弱受制于吸附氣量的多少。所以，表面擴散系數以及吸附氣量是影響頁巖氣表面擴散效應的2 個關鍵參數。圖13 為考慮不同表面擴散效應下，多孔介質表觀滲透率隨吸附氣量的變化情況。從圖13 可以看出，表面擴散效應對表觀滲透率的貢獻會隨著吸附氣量的增大而增加，吸附氣量越大時，不同表面擴散系數對表觀滲透率的影響也愈發顯著。

圖13 不同表面擴散系數下的表觀滲透率隨吸附氣量增長的變化情況Fig.13 Variance of apparent permeability with ascending adsorbed gas concentration under different surface diffusivity

3 結論

（1）重構了天然微裂縫發育的頁巖儲層多孔介質，考慮頁巖氣滲流過程中的滑脫、吸附/解吸、表面擴散效應等運移機制對廣義格子Boltzmann 模型進行了修正，并基于經典REV 尺度格子Boltzmann 模型驗證方法對重構模型的正確性進行了驗證。

（2）應用重構模型對頁巖氣滲流過程進行了模擬。結果表明，當儲層天然微裂縫較發育時，絕大部分氣體沿天然微裂縫竄進，天然微裂縫成為頁巖氣滲流的主要通道。通過對天然微裂縫物性參數進行敏感性分析發現，裂縫密度對于頁巖氣滲流的促進作用最為顯著；隨著裂縫密度增大3～4 倍，儲層表觀滲透率可增大10 倍以上；裂縫長度及裂縫開度對表觀滲透率的貢獻則較小。

（3）參數敏感性分析表明，努森數（Kn）是影響頁巖氣滲流特征的主要參數，隨著Kn的增大，克氏效應顯著，氣體流速隨之增大，特別當Kn>0.1 時，表觀滲透率增幅明顯增大；表觀滲透率會隨著吸附氣量的增大而減小，特別當儲層壓力較低時，該現象更為顯著；表面擴散效應也是影響頁巖氣滲流的主要因素之一，同等條件下考慮氣體表面擴散效應的表觀滲透率較忽略該因素可提升2～5 倍，但提升作用受制于吸附氣量的多少。