復雜微通道氣體流動的格子Boltzmann模擬

姚 杰， 董 波， 李維仲

(大連理工大學 海洋能源利用與節(jié)能教育部重點實驗室，大連 116024)

1 引 言

目前，在油氣開采[1-3]和微機電裝置[4]等領域普遍存在的微通道內氣體流動現象受到廣泛關注。在微通道中，由于特征長度H很小，能達到微米級甚至納米級，氣體分子平均自由程λ通常為10-7m～10-8m(以甲烷分子為例，273 K和0.1 MPa時，λ=0.058 μm)，此時努森數Kn=λ/H>10-3，流動處于滑移區(qū)，在研究中需要考慮努森層對流動的影響。努森層是氣體流經壁面時在近壁面處存在的一個局部非平衡區(qū)。在努森層內，氣體密度較小，分子間碰撞相對較少，流體的動力學行為主要受氣體分子和固壁之間的碰撞影響，流體粘度顯著降低。由于固壁作用，近壁面氣體速度不同于壁面速度，這意味著會出現流動不連續(xù)或速度滑移。此時流體不能按照連續(xù)介質處理，基于連續(xù)介質假設的N-S(Navier-Stokes)方程組不能描述氣體流動狀態(tài)。

通常，當流動處于滑移區(qū)(10-3≤Kn≤0.1)時，在努森層影響下，氣體流動狀態(tài)可以用N-S方程加滑移邊界來描述，N-S方程組依舊可以很好地描述主流區(qū)內流體流動。當流動處于過渡區(qū)(0.1

格子Boltzmann方法(LBM)是一種基于分子動理論的介觀數值模擬方法，計算效率高，且易于處理復雜邊界，尤其適用于孔隙尺度內流動模擬[10]。目前已有學者利用LBM研究微氣體流動問題。Zhao等[1]建立了一個微尺度多松弛時間LB模型，用于模擬氣體在各向異性多孔介質中的流動。姚軍等[3]采用考慮正則化過程的微尺度LB模型研究微裂縫對致密多孔介質中氣體滲流的影響機制。上述研究采用的動力學邊界Maxwell反射邊界是一個完全漫反射邊界[11]，模擬時得到的滑移速度會比實際滑移速度大。Chai等[12]將Maxwell漫反射邊界與半步長反彈結合提出了新的曲面滑移邊界，可用于復雜微通道氣體流動。Chen等[13]提出了考慮Klinkenberg效應的表征單元體尺度LB模型，并采用該模型研究了微裂縫對頁巖氣滲流的影響。Suga等[10,14]討論了微尺度LB模型在Poiseuille流、Couette流和組合微通道流等算例中的應用，并對比了采用不同離散速度，邊界條件模型的計算結果。以上采用的模型沒有考慮孔隙截面變化帶來的Kn數的變化。

盡管前人已經利用LBM對氣體微流動進行了大量研究，但是采用模型仍具有一定局限性。實際多孔介質中，孔隙結構蜿蜒變化，交錯縱橫，因此孔隙內部各點的特征尺寸也不一樣，且對于含有微縫隙的介質，特征尺寸跨度很大。由于孔隙的特征尺寸與Kn數息息相關，而Kn數又決定著流體流動狀態(tài)，因而采用單一特征Kn數來代表整個多孔介質內每一點的Kn數勢必會帶來較大的誤差。

本文建立了適用于微通道氣體流動的多松弛格子Boltzmann模型，考慮了努森層內分子平均自由程的修正，采用半步長反彈和Maxwell漫反射邊界組合的曲面邊界格式，將局部Kn數引入演化方程。該模型可以模擬滑移區(qū)和過渡區(qū)的滲流問題。

2 數值方法

2.1 適用于微通道流動模擬的格子Boltzmann模型

格子Boltzmann方程(LBE)可以看作Boltzmann方程的一種離散格式[15,16]，且已證明是研究微尺度流動的有效工具[17-20]。演化方程為

fi(x+ciδt,t+δt)-fi(x,t)=Ωi(f)+δtFi

(i=0,1,…,b-1)(1)

式中fi(x,t)為速度為ci的氣體分子在時刻t空間位置x處的分布函數；Fi為系統(tǒng)受到的體積力的源項，δt為時間步長，b為空間的離散方向數，對于D2Q9模型，b=9。離散速度為

(2)

式中c=δx/δt為格子速度，本文c=1。Ωi(f)為離散的碰撞算子，且對于LBGK模型[21]為

(3)

(4)

(5)

LBGK模型在很多領域得到廣泛的應用，然而最近的一些研究指出LBGK模型在計算氣體微流動具有一定局限性[22]，主要在于應用動力學邊界時會引入非物理因素。根據氣體動理論，滑移速度us只取決于Kn數以及氣體與壁面的相互作用系數，但是在LBGK模型中，滑移速度還和網格分辨率相關。而多松弛模型(MRT)則可以通過選擇合適的松弛時間避免這種非物理因素的干擾。因此本文采用多松弛(MRT)模型碰撞算子[23]為

(6)

式中M為一個b×b階轉換矩陣，用于將速度空間的分布函數fi投影到矩空間m=Mf，其中f=(f0,f1,…,fb - 1)T，S=diag(τ0,τ1,…,τb - 1)-1為非負對角矩陣，τi為第i個矩的松弛時間。轉換矩陣為

(7)

對應的松弛矩陣

(8)

式(1)的演化過程可以寫作兩步為

碰撞

(9)

遷移

(10)

離散的外力源項Fi可表示[24]為

(11)

(12)

式中G為系統(tǒng)受到的體積力的合力。

2.2 松弛時間及曲面滑移邊界

2.2.1 松弛時間

在LBM中，宏觀尺度下流動特征參數是雷諾數，通常由Re數來確定松弛時間；而對于微尺度下的氣體流動，特征參數是努森數，需要通過Kn數來確定松弛時間。由分子動理論以及格子Boltzmann方法中粘度與松弛時間的關系可推得松弛時間τ與Kn數的基本關系式為

(13)

式中N為模型中特征尺寸所占的網格數。由于努森層的存在，微通道系統(tǒng)中近壁面處氣體的運動是受限的，分子平均自由程小于遠離壁面處分子的平均自由程。開放系統(tǒng)和無邊界系統(tǒng)中分子平均自由程的關系為[25]

(14)

式(14)表明，受限程度隨Kn數的變大而顯著。當Kn數趨于0，也就是努森層占比很小時，此時有效平均自由程又恢復到開放系統(tǒng)的平均自由程。

將有效平均自由程代替平均自由程以修正Kn數較大時努森層對主流層的影響，可以在Kn數變化更大的范圍內模擬氣體的微流動現象。因此，本文微通道中流體流動的MRT模型的松弛時間為

(15)

其他參數可取τe=1.0，τe=τε=1/1.4，τd=1.0，τq的取值和滑移邊界相關。

2.2.2 曲面滑移邊界

將滑移邊界引入到LBM的方法主要有兩種。第一種方法是用滑移速度來反映邊界處滑移長度的影響[26]。第二種方法是采用組合滑移格式。對于處于滑移區(qū)的氣體流動，可采用LBM邊界包括鏡面反射(SR)、標準反彈(BB)和麥克斯韋漫反射(MD)等格式組合的混合格式[12]。根據計算需要，本文采用MD和BB的組合格式。該格式具有清晰的物理意義，即運動到固壁的氣體分子，一部分直接反彈回去，另一部分與固壁相互作用并達到Maxwell平衡分布再反射回去。在原格式的基礎上，考慮努森層以及局部Kn數的影響。曲面滑移邊界可表示為

(16)

(17)

根據具有解析解的平板Poiseuille流求取組合系數r。MRT-LBM結合式(16)可推得近壁面相鄰流層速度之間的關系[12]為

u2=Au1+Bδta

(18)

式中 系數A和B為

A=(3-2r+rτs)/(1-r+rτs)

(19)

B=[(11/2-2τq-8τs+4τqτs)r+4τq+4τs-

8τqτs-5]/[(2τs-1)(1-r+rτs)]

(20)

Poiseuille流的解析解為

u=(4U0y/H)(1-y/H)+us

(21)

式中y=(j-Δ)δx。對于半步長反彈，Δ=1/2，可得兩流層之間的速度關系式。

u1=(2u0δx/H)[(2H-δx)/(2H)]+us

u2=(6u0δx/H)[(2H-3δx)/(2H)]+us

(22)

代入式(18)，并令ε=δx/H，可得

(23)

聯立式(19,23)，可得

us/u0=[(2r-2)(1-2τs)/(1+r)]ε+

{[16(τq-0.5)(τs-0.5)-3]/3}ε2

(24)

此速度是根據滑移邊界條件求得的滑移速度，與實際滑移速度比較即可得到組合系數r。實際平板Poiseuille流滑移速度可表示為[23]

Us=4C1Kn+8C2Kn2

(25)

(26)

2.3 局部Kn數

多孔介質是由固體骨架和相互連通的孔隙所構成，孔隙內部各點的特征尺寸也不一樣，且對于含有微縫隙的材料，特征尺寸跨度很大。Zhao等[1]提出了一個獲取孔隙特征尺寸的方法，通過邊緣細化算法[27,28]提取孔隙中軸，然后根據獲取的孔隙中軸計算孔隙的局部特征尺寸。經細化算法處理后的多孔介質如圖1所示。對于每一個中軸點，該點處的特征尺寸為中軸點到最近固壁距離的2倍。這樣便得到了孔隙內任意的特征尺寸H(r)。

選取合適的參考特征尺寸Href、參考壓力Pref和參考溫度Tref，可確定參考努森數Knref，則局部努森數Kn(r,t)可表示為孔隙內某一點溫度T(r,t)、壓力P(r,t)以及特征尺寸H(r)的函數，可表示為

(27)

代入式(15)，可得局部松弛時間

(28)

圖1 邊緣細化算法獲取中軸Fig.1 Central axis obtained by edge thinning algorithm

3 模型驗證

為了驗證所構建模型的正確性，本文采用Poiseuille流對模型在過渡區(qū)以及滑移區(qū)的計算滑移速度準確性進行驗證。

平板長度為L，寬為H，區(qū)域網格為100×15，上下邊界為固體節(jié)點，左右邊界為周期性邊界，固體邊界采用滑移邊界，令σ=1.0，即固壁是全反射的。水平方向添加外力場Fx=5×10-5，由式(11)引入模型中，收斂條件為

(29)

式中ux和uy分別為通道中氣體沿x方向和y方向的流動速度。

本文對不同Kn數下的氣體流動情況進行模擬，將模擬結果與直接模擬蒙特卡洛方法(DSMC)、信息保存-直接模擬蒙特卡洛方法(IP -DSMC)[29]以及未進行有效平均自由程修正和邊界采用DBB格式[30]的計算結果進行比較。當模擬達到穩(wěn)態(tài)后，統(tǒng)計x=L/2處的無因次速度分布如圖2所示。

圖2 不同Kn數下Poiseuille流無因次速度分布Fig.2 Dimensionless velocity profile comparisons of Poiseuille flows under different Kn numbers

本文在滑移區(qū)和過渡區(qū)模擬結果與平板滑移模型解析解吻合較好。當Kn數較小時，不采用有效平均自由程修正的LBM結果與解析解和DSMC方法的解吻合較好，但是在Kn數較大的過渡區(qū)，其滑移速度偏小，這是因為努森層占主流區(qū)比例變大，壁面對流動影響愈加明顯，此時表觀滑移速度將明顯小于實際滑移速度。這與本文分析結果一致，且進一步證明有效平均自由程修正的必要性。

4 微尺度氣體滲流規(guī)律研究

應用本文提出的模型研究氣體在各向異性多孔介質中的流動，如圖3所示。模擬區(qū)域為400×300，其中孔隙率ε=0.597，孔隙實際物理尺寸為7.5 μm～0.5 μm?？紫吨辛鲃拥臍怏w為273 K和0.1 MPa狀態(tài)下的甲烷，對應局部Kn數范圍為0.0077～0.116，平均努森數Knm=0.0192。本文通過改變網格分辨率來調整Kn的數值。由達西定律

k=μ·Q·ΔL/(A·ΔP)

(30)

圖3 多孔介質Fig.3 Porous media

式中Q為體積流量，A為通道的截面積，對于二維模型，A=Ly，μ為氣體的動力粘度。對于本算例，體積流量Q的計算式為

Q=∑uxδy

(31)

Knm平均努森數可計算為

Knm=∑Kn(r)、N

(32)

式中N為孔隙點的個數。選取模擬的Knm數范圍為0.0001～3.8。

4.1 不同Knm數下多孔介質內氣體速度分布

圖4為穩(wěn)態(tài)時多孔介質內氣體無因次流速的縱向分布情況。無因次流速為孔隙中某一點的流速和孔隙最大流速的比值，可以反映在Knm數變化時孔隙內流速分布的變化情況。

圖4 孔隙內無因次速度分布情況Fig.4 Dimensionless velocity distribution in pores

可以看出，當Knm數較小時(圖4(a))，多孔介質中主流通道較少，此時氣體主要流經相對大的孔隙，多孔介質的非均質性對氣體的流動影響明顯；而隨著Kn數的增加，小孔隙中的氣體流動受努森層的影響變大，相對大孔隙而言，小孔隙內的氣體流速增加，主流通道明顯變多(圖4(f))，多孔介質的非均質性影響減弱。為了進一步研究Knm對孔隙中速度的影響，提取了不同Knm下截面x=115和x=240上的速度分布情況，如圖5所示。氣體的流速分布在各個截面上仍呈現近似的拋物線分布，這與Poiseuille流是一致的，隨著多孔介質內Knm的增加，通道內各點的無因次速度相應變大，且可以明顯看到小孔隙內氣體速度增加幅度更大，孔隙愈小，增加幅度越大。這是由于小孔隙的局部Kn數更大，微尺度效應更明顯。隨著Knm數的增加，大孔隙和小孔隙流動阻力均減弱，小孔隙中流動阻力減小得更多。這與Klinkenberg對多孔介質內氣側滲透率大于固有滲透率且展現出壓力依賴性的解釋一致[31]。

圖5 不同截面上流體無因次速度分布情況Fig.5 Dimensionless velocity distribution of the fluid at different cross -sections

4.2 不同Knm數下多孔介質滲流特性

為了研究Knm數對孔隙內滲流特性的影響，首先，將孔隙內氣體的流速無因次化，然后，計算多孔介質的等效流速v=Q總/Ly，其中Q總為出口截面的體積流量，Ly為出口截面的長度。本文為了直觀觀察孔隙中氣體的流動情況，令孔隙內速度大于等效流速的孔隙為高滲點，將孔隙分成高滲區(qū)域和低滲區(qū)域。定義高滲區(qū)域與孔隙面積的比值為多孔介質的高滲面積比R，高滲面積比隨Knm數的變化如圖6所示。結果表明隨著Knm的增加，高滲面積比逐漸增大。

當Knm數增大時，小孔隙中高滲區(qū)域相比大孔隙增加更加明顯，且小孔隙的高滲區(qū)域相比于大孔隙更靠近壁面。這是由于小孔隙中滑移速度更大，速度分布更均勻所致。

圖6 高滲面積比隨Knm的變化Fig.6 Variation of high permeability area ratio with Knm

在此基礎上，本文分別計算在不同Knm數下孔隙表觀滲透率ka，并和Klinkenberg的宏觀滲透率模型的計算值進行比較，ka的表達式為

ka=k∞(1+6c·Kn)

(33)

式中k∞為固有滲透率，c為相稱性因子。為了便于比較，令kr=ka/k∞，用kr來表征表觀滲透率的大小。對于相稱性因子，取值與多孔介質有關，平板模型一般取c=1.0，對于其他多孔介質，一般c<1.0，當c=0.4，與本文算例計算結果吻合較好，這也證明了本文構建模型用于計算多孔介質的正確性。本文計算結果與Klinkenberg模型計算值的比較列入表1。

表1 不同Knm下相對表觀滲透率比較Tab.1 Comparison between current results and result obtained by Klinkenberg model in terms of relative apparent permeability at different Knm

5 結 論

(1) 考慮努森層的影響，本文構建了適用于滑移區(qū)和過渡區(qū)的MRT-LBM模型，通過算例Poiseuille流對模型進行了驗證，所得結果與用直接蒙特卡洛方法和分子模擬獲得的結果吻合較好。本文的模型能夠較好模擬滑移區(qū)以及過渡區(qū)的氣體滑移現象，并且在Kn數較大情況下(Kn>0.1)，比不考慮有效分子平均自由程修正的滑移速度更加準確。

(2) 本文利用所建立的模型研究了微尺度下非均質多孔介質內氣體滲流特性。結果表明多孔介質的表觀滲透率隨著平均努森數Knm增加而增加，這與Klinkenberg宏觀滲流模型理論一致。隨著Knm數的增加，小孔隙中氣體滲流速度相對大孔隙內氣體滲流速度有所增加，主滲流通道數增加。通過分析不同Knm數下多孔介質內高滲區(qū)域分布發(fā)現，隨著Knm數增加，孔隙中高滲區(qū)域占總孔隙比例增加，且高滲區(qū)優(yōu)先從小孔隙近壁面處增加，這是因為小孔隙中努森層占比更大，氣體分子和壁面間碰撞所占比例增大，滑移流速增大，氣體流動阻力減小，截面速度分布更加均勻。