孔喉結構對致密氣微尺度滲流特征的影響

張烈輝 劉香禺 趙玉龍 周 源 單保超

1.“ 油氣藏地質及開發工程”國家重點實驗室?西南石油大學 2. 中國石油西南油氣田公司勘探開發研究院3. 煤燃燒國家重點實驗室?華中科技大學

0 引言

目前，對于致密氣藏的滲流機理通常根據低滲透氣藏的開發理論或實驗方法來加以模擬和研究，但傳統的數值模擬方法和宏觀實驗手段無法捕捉氣體在致密儲層中的微尺度流動特征；同時，微納米級孔喉是致密儲層儲集和運移流體的主要空間與通道，而關于孔喉結構對致密氣微尺度滲流特征影響的研究卻較少。為了克服宏觀方法的局限性，采用介觀尺度的格子Boltzmann方法（LBM）對致密氣的微觀滲流進行模擬，進而研究孔喉結構對滲流特征的影響，不失為一種有效的方法。由于格子Boltzmann方法具有介觀特性且并行性好，在致密砂巖氣、頁巖氣和煤層氣等非常規天然氣的微尺度滲流機理研究中得到了應用。吳子森等[1]結合四參數隨機生長法，采用LBGK-D2Q9模型研究了致密氣的微尺度效應。王華龍等[2]考慮氣體在多孔介質中的微尺度效應，將多孔介質等效為相互平行的圓柱體，將單通道LBM模型拓展并應用到孔隙群里氣體滲流的數值模擬，為LBM深入研究氣體在多孔介質中的滲流問題奠定了基礎。任嵐等[3]、姚軍等[4]采用LBGK-D2Q9模型對頁巖氣藏微納米孔隙中甲烷氣體的流動進行了模擬研究，得到了不同溫度、壓力和孔道直徑對頁巖氣微尺度流動的影響規律。張烈輝等[5]基于表征體元尺度的LBM模型，考慮滑脫效應，對頁巖氣在儲層中的流動進行了數值模擬研究，研究結果表明頁巖氣在有機質中的流速略大于無機質中的流速。Fathi等[6]采用LBM模擬了二維Poiseuille流動，對Klinkenberg滑移理論進行了修正，并提出了雙重滑脫模型。現有的文獻報道通常針對平直管道或隨機生成的多孔介質展開研究，未能體現出孔喉結構差異對氣體微尺度滲流的影響。為此，筆者從致密氣藏實際的溫度、壓力及儲層孔喉特征尺寸出發，研究致密氣的流態及采用LBM模擬致密氣流動的合理性，然后考慮微尺度效應、滑脫效應等因素的影響，基于LBGK-D2Q9模型建立了致密氣流動模型，進而針對平直通道開展流動模擬計算，將模擬結果與解析解、數值解的計算結果進行對比以驗證模型的正確性，在此基礎上，研究了孔喉結構對致密氣微尺度滲流特征的影響。

1 致密氣流動特征

Nelson[7]指出致密砂巖儲層孔喉直徑介于30～2 000 nm；鄒才能等[8]、楊智等[9]通過大量實驗數據總結出致密砂巖儲層孔喉直徑介于40～700 nm。在微納米尺度，克努森數（Kn）是氣體流動的特征參數，定義為分子平均自由程與流動通道特征長度的比值[10]，即

式中λ表示分子平均自由程，m；H表示流動通道特征長度，m。

分子平均自由程的計算式[11]為：

式中m表示分子質量，kg，考慮致密氣為單組分甲烷，取值為2.658×10-26kg；ρg表示氣體密度，kg/m3，本文采用美國國家標準與技術研究院（National Institute of Standards and Technology，簡稱為NIST）化學數據庫中數據；d表示分子直徑，m，取值為0.414×10-9m。

根據Kn可以將氣體流態劃分為以下4種：連續流、滑脫流、過渡流和分子自由流[12]。當流態為連續流和滑脫流時，流體流動均遵循Navier-Stokes方程，其中流態為滑脫流時需要對邊界進行滑移修正處理；流態為過渡流時，連續介質假設不再成立，相應的基于該假設的傳統計算流體力學方法（CFD）無法模擬對應尺度下流體的流動；對于分子自由流，通常需要采用分子動力學方法（MD）進行模擬。

為明確致密氣在地層條件下的流態，繪制了不同溫度（293.15 K、373.15 K）、壓力（介于1～70 MPa）和儲層孔喉直徑（30 nm、1 000 nm、2 000 nm）下致密氣的Kn變化曲線（圖1）。從圖1可以看出，當壓力大于3 MPa時，Kn均小于10-1，氣體流態為滑脫流和弱連續流，出現過渡流的概率較低。由于LBM是求解Navier-Stokes方程組的一種特殊離散方法。因此，采用LBM模擬致密氣的微尺度流動是可行的。

圖1 致密氣流態劃分圖版

2 微尺度流動模型的建立

2.1 物理模型

常見的喉道類型概括為以下4種：孔隙的縮小部分、可變斷面的收縮部分、片狀或彎片狀喉道及管束狀喉道[13]。為了便于研究孔喉結構對致密氣流動的影響，對物理模型進行如下簡化，考慮孔隙為二維平直通道，將喉道抽象為規則方形凸起，如圖2所示。

圖2 孔喉結構簡化示意圖

2.2 格子Boltzmann模型

1992年，Qian等[14]提出了格子Boltzmann方法的基本模型——DdQm模型，得到了廣泛應用。在模擬單相單組分的二維流動時，通常采用九速離散的LBGK-D2Q9模型（圖3）。不考慮外加作用力時，LBGK-D2Q9模型的演化方程如式（3）所示。在LBM中所有物理量的單位為格子空間單位，均無量綱。

圖3 描述致密氣微觀流動的格子Boltzmann模型示意圖

式中f表示粒子分布函數；下標i表示離散速度方向，i=0,1,…,8； 表示格子空間中某位置；表示不同離散速度方向格子速度；δt表示格子空間時間步長；t表示格子空間時間；τ表示無因次松弛時間；上標eq表示平衡態。

在LBGK-D2Q9模型中，平衡態分布函數表示為：

其中

式中ρ表示格子空間宏觀密度；wi表示權重系數；c表示格子速度； 表示格子空間宏觀速度；cs表示格子聲速；δx表示格子空間中的網格間距，一般其取值與δt數值相等，此時c=1。

格子空間宏觀密度ρ和宏觀速度 分別表示為：

對于LBGK-D2Q9模型，格子速度 和權重系數wi分別取值為：

在格子Boltzmann模型中，運動黏度作如下處理可使得模型在計算時具有二階精度[15]，即

式中 表示格子空間運動黏度。

同時，對于單組分單相流體，格子空間中壓力p與密度、格子聲速的關系式為：

2.3 邊界條件

對于微尺度流動問題，邊界條件起著極為重要的作用。與宏觀模擬方法不同，在LBM中邊界條件設置的對象為分布函數。對微尺度通道入口和出口采用Guo等[16]提出的非平衡外推格式，上下壁面采用反彈與鏡面反射組合格式（BSR）[17]以實現對滑移邊界條件的處理。

對于二維平直通道內流動（圖3），氣體在進出口壓差作用下，沿x軸正方向流動。以左側x=0位置處網格節點為例，在遷移步后未知的分布函數為f1、f5、f8，采用非平衡外推格式，其邊界處分布函數可以表示為平衡態和非平衡態兩部分之和，即

式中上標neq表示非平衡態；下標i=1,5,8。

對于壓力驅動的微尺度流動，進口位置處壓力已知，而速度未知。在非平衡外推格式中，未知的平衡態分布函數采用入口處壓力以及相鄰節點的速度進行求解，即

對于非平衡態部分，使用相鄰流體節點處的非平衡態分布函數近似處理，即

為了實現對壁面滑移的模擬，上下壁面采用BSR邊界條件處理格式。由于采用Half-Way類型的反彈格式可使粒子運動圖像清晰，計算精度高于標準反彈格式[10]。因此采用Half-Way類型的BSR邊界條件處理格式。以下邊界為例，遷移步后未知的分布函數為f2、f5、f6，此時可得

前人的研究中組合系數r多為定值[1,18-19]，但是根據Guo等[20]的研究結論，為消除數值離散效應，組合系數r應遵循如下選取規則，即

式中A1、A2分別表示與氣固相互作用有關的系數；N表示網格數。

Guo等[21]在廣義二階滑移邊界條件中給出A1、A2的計算式為：

式中α表示切向動量調節系數，一般取值為1[4,22]。

顯然，組合系數r不僅依賴于氣體與壁面之間的相互作用系數A1、A2，也與Kn和網格數相關。

2.4 無因次松弛時間的選取

無因次松弛時間τ與Kn之間的關系確定是LBM模擬氣體微尺度流動的一個基本問題，其中Kn是流動特征控制參數，τ與Kn之間的關系式[10]為：

實際上，氣體在微納米尺度空間中流動時還受到Knudsen層的影響。由于固體壁面的存在使得靠近壁面的氣體分子平均自由程被截斷，表現出不同于遠離壁面的氣體分子的運動特征，受到影響的這部分氣體層則稱為Knudsen層。Knudsen層的存在使得其氣體分子平均自由程減小，當流動通道的特征長度遠大于Knudsen層厚度時，Knudsen層對流體流動的影響不明顯；但隨著流動通道孔徑減小，Knudsen層的厚度與通道特征長度的比值增加，其對流體流動的影響增大，表現出了不同于在常規大尺度通道中流動的特征，即微尺度效應。當氣體流態為滑脫流時，分子之間的碰撞概率仍大于分子與壁面碰撞的概率，但Knudsen層的存在使得氣體在邊界的流動出現滑脫現象[23]。本文采用Guo等[24]提出的方法，將Knudsen層的影響考慮為使整個流動區域的氣體分子平均自由程減小，此時無因次松弛時間τ和Kn的關系式為：

式中 表示平均分子自由程修正函數。

反映了孔隙壁面對氣體分子運動過程的截斷作用，其表達式為：

考慮到本文研究的微尺度滲流通道為變截面通道，流動通道的特征長度及Kn在不同位置會存在差異，因此式（19）中的網格數N和Kn均為局部參數。相應Kn根據下式進行計算，即

式中下標out代表通道出口中心位置。

3 模型驗證

基于二維平板流動，假設入口處壓力為pin，出口處壓力為pout，將不同工況下本文模型的模擬結果分別與解析解及本文參考文獻[25]中數值解的計算結果進行對比，以驗證本文模型的正確性。

3.1 解析解驗證

進出口壓力分別設定為pin=1.01、pout=1.00，在運動黏度 分別取0.20、0.30的情況下，不考慮邊界Knudsen層的影響，采用無滑移LBGK-D2Q9模型模擬流體流動。模擬中網格劃分為Nx×Ny=300×100，計算得到出口處y方向的速度（Ux）剖面，并與解析解計算結果進行對比，結果顯示兩者具有很好的一致性（圖4）。

圖4 不同工況下本文模型數值解與解析解計算的速度剖面圖

3.2 數值解驗證

為進一步驗證模型的正確性，采用本文模型模擬了Kn為0.019 4和0.194 0兩種情況下出口處y方向的無因次速度剖面以及沿程無因次壓力偏移量，并與本文參考文獻[25]中采用DSMC、IP方法計算得到的結果進行對比。結果顯示在不同Kn條件下，本文模型計算得到的出口處y方向的無因次速度剖面與DSMC、IP方法的計算結果吻合程度均較高（圖5）；沿程無因次壓力偏移量與采用DSMC、IP方法計算的結果雖存在少許差異，但整體趨勢相符（圖6）。

圖5 不同Kn下DSMC、IP方法及本文模型計算的無因次速度剖面對比圖

4 結果分析

圖6 不同Kn下DSMC、IP方法及本文模型計算的沿程無因次壓力偏移量結果對比圖

在微尺度流動中，Kn是主控因素[22]，因此本文直接給定出口處Kn為0.019 4；采用壓差驅動，定進口壓力pin為1.1，出口壓力pout為1.0，模擬時間步為100 000步以保證流動達到穩定狀態。模擬計算中，網格數Nx為400、Ny為102，對直孔隙（對應孔喉比為5∶5）及孔喉比分別為5∶4、5∶3、5∶2和5∶1情況下的速度和壓力分布以及相關參數進行分析。

4.1 Kn分布

由于沿程壓力及孔喉比的變化，在不同位置處Kn的大小也存在差異。如圖7所示，相比于沿程壓力變化對Kn的影響，孔喉比的影響程度更大，在喉道處Kn將成倍增加；當孔喉比一定時，Kn沿喉道呈緩慢上升的趨勢，且孔喉比越大，Kn上升的趨勢越明顯。

圖7 不同孔喉比下沿程Kn分布圖

4.2 壓力分布

如圖8所示，在孔喉比為5∶5時氣體壓力呈近線性分布；而存在喉道時，近線性壓力分布規律被打破，氣體在接近喉道進口處形成相對高壓區，在離開喉道后形成相對低壓區，且孔喉比越大，壓力越低；在喉道位置處氣體壓力下降明顯，在喉道兩側氣體壓力降低于喉道內，并且隨著孔喉比的增加，喉道兩側氣體壓力變化越來越平緩。由于喉道的存在所導致的強非線性壓力分布在一定程度上降低了氣體流動動力，從而阻礙了氣體的流動。

圖8 不同孔喉比下孔道中心位置處壓力沿程變化圖

圖9 不同孔喉比下喉道內壓力降占比結果統計圖

4.3 速度分布

如圖10所示，當孔喉比大于1時，由于喉道結構的存在，氣體流動的橫截面積減小，導致喉道內外氣體的流速存在明顯差異，氣體流速在喉道末端位置處出現最大值，這是遵循質量守恒定律的必然結果。如圖11所示，在不同孔喉比下出口處y方向的速度分布剖面依然呈現為中間高、兩邊低的趨勢，并且邊界滑移速度隨著孔喉比的增加而下降。總的來說，隨著孔喉比增加，不同位置處的流速均較大程度地降低，從而導致流量降低。

圖10 不同孔喉比下孔道中心位置處沿程速度變化曲線圖

圖11 不同孔喉比下出口處速度剖面圖

5 結論

1）當壓力介于3～70 MPa、溫度介于293.15～373.15 K時，致密氣藏中Kn小于0.1，氣體流態為滑脫流和弱連續流，采用考慮微尺度效應和邊界滑移效應的LBGK-D2Q9模型進行致密氣流動的模擬是合理的。

2）本文模型的模擬結果與解析解及文獻中DSMC、IP方法等數值解的計算結果吻合程度均較高，證實本文模型可靠。

3）流動通道的特征長度對Kn的影響遠大于壓力變化對其產生的影響，當孔喉比一定時，Kn沿喉道呈緩慢上升的趨勢，且孔喉比越大，Kn上升的趨勢越顯著。

4）喉道的存在使得孔隙中壓力的非線性分布特征顯著，且壓力降主要位于喉道內，同時孔喉比越大，喉道內的壓降幅度越大。

5）壓力的非線性分布使得氣體的流動速度顯著降低，從而降低了流動通道內氣體的質量流量。