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

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

0 引言

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

1 致密氣流動特征

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

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

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

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

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

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

圖1 致密氣流態(tài)劃分圖版

2 微尺度流動模型的建立

2.1 物理模型

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

圖2 孔喉結(jié)構(gòu)簡化示意圖

2.2 格子Boltzmann模型

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

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

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

在LBGK-D2Q9模型中，平衡態(tài)分布函數(shù)表示為：

其中

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

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

對于LBGK-D2Q9模型，格子速度 和權(quán)重系數(shù)wi分別取值為：

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

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

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

2.3 邊界條件

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

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

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

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

對于非平衡態(tài)部分，使用相鄰流體節(jié)點處的非平衡態(tài)分布函數(shù)近似處理，即

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

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

式中A1、A2分別表示與氣固相互作用有關(guān)的系數(shù)；N表示網(wǎng)格數(shù)。

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

式中α表示切向動量調(diào)節(jié)系數(shù)，一般取值為1[4,22]。

顯然，組合系數(shù)r不僅依賴于氣體與壁面之間的相互作用系數(shù)A1、A2，也與Kn和網(wǎng)格數(shù)相關(guān)。

2.4 無因次松弛時間的選取

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

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

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

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

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

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

3 模型驗證

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

3.1 解析解驗證

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

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

3.2 數(shù)值解驗證

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

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

4 結(jié)果分析

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

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

4.1 Kn分布

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

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

4.2 壓力分布

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

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

圖9 不同孔喉比下喉道內(nèi)壓力降占比結(jié)果統(tǒng)計圖

4.3 速度分布

如圖10所示，當(dāng)孔喉比大于1時，由于喉道結(jié)構(gòu)的存在，氣體流動的橫截面積減小，導(dǎo)致喉道內(nèi)外氣體的流速存在明顯差異，氣體流速在喉道末端位置處出現(xiàn)最大值，這是遵循質(zhì)量守恒定律的必然結(jié)果。如圖11所示，在不同孔喉比下出口處y方向的速度分布剖面依然呈現(xiàn)為中間高、兩邊低的趨勢，并且邊界滑移速度隨著孔喉比的增加而下降?？偟膩碚f，隨著孔喉比增加，不同位置處的流速均較大程度地降低，從而導(dǎo)致流量降低。

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

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

5 結(jié)論

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

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

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

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

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