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

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

0 引言

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

1 致密氣流動(dòng)特征

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

式中λ表示分子平均自由程，m；H表示流動(dòng)通道特征長(zhǎng)度，m。

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

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

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

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

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

2 微尺度流動(dòng)模型的建立

2.1 物理模型

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

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

2.2 格子Boltzmann模型

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

圖3 描述致密氣微觀流動(dòng)的格子Boltzmann模型示意圖

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

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

其中

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

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

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

在格子Boltzmann模型中，運(yùn)動(dòng)黏度作如下處理可使得模型在計(jì)算時(shí)具有二階精度[15]，即

式中 表示格子空間運(yùn)動(dòng)黏度。

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

2.3 邊界條件

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

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

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

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

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

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

前人的研究中組合系數(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的計(jì)算式為：

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

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

2.4 無因次松弛時(shí)間的選取

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

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

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

反映了孔隙壁面對(duì)氣體分子運(yùn)動(dòng)過程的截?cái)嘧饔?，其表達(dá)式為：

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

式中下標(biāo)out代表通道出口中心位置。

3 模型驗(yàn)證

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

3.1 解析解驗(yàn)證

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

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

3.2 數(shù)值解驗(yàn)證

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

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

4 結(jié)果分析

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

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

4.1 Kn分布

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

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

4.2 壓力分布

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

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

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

4.3 速度分布

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

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

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

5 結(jié)論

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

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

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

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

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