微納多孔結(jié)構(gòu)中稀薄氣體流動(dòng)滲透率的解析型預(yù)測(cè)模型

楊光，程鑫，王崢，王曄，張良俊，吳靜怡

（1上海交通大學(xué)機(jī)械與動(dòng)力工程學(xué)院，上海 200240； 2上海宇航系統(tǒng)工程研究所，上海 201108）

引 言

多孔結(jié)構(gòu)中的稀薄氣體流動(dòng)機(jī)理在石油化工、能源利用、航空航天等領(lǐng)域具有廣闊的應(yīng)用前景[1-3]，如航天推進(jìn)劑的增壓輸送、非常規(guī)油氣開采、質(zhì)子交換膜燃料電池、CO2封存、核廢料處理等。近年來，3D 打印、MEMS 技術(shù)的發(fā)展也推動(dòng)著對(duì)微納尺度氣體流動(dòng)特性的深入研究[4-5]。當(dāng)多孔結(jié)構(gòu)的孔隙尺度足夠小，或者氣體工質(zhì)處于低壓狀態(tài)時(shí)，氣體分子的平均自由程與孔隙的特征長(zhǎng)度相當(dāng)，氣體分子與固體壁面的碰撞頻率和氣體分子間的碰撞頻率相近，即產(chǎn)生稀薄氣體效應(yīng)。該現(xiàn)象通常用Knudsen 數(shù)（Kn）來量化，它表征了氣體分子平均自由程（λ）與流動(dòng)特征長(zhǎng)度（Lc）的比值，并可將氣體的流動(dòng)大致分為四種狀態(tài)：連續(xù)流（Kn≤0.001），滑移流（0.001＜Kn≤0.1），過渡流（0.1＜Kn≤10）和自由分子流（Kn＞10）。對(duì)于非連續(xù)區(qū)的高Kn流動(dòng)，壁面流動(dòng)速度的“無滑移”假設(shè)不再成立，從而導(dǎo)致氣體的表觀滲透率高于材料的固體滲透率，這種現(xiàn)象也被稱為Klinkenberg 效應(yīng)。而目前常用的多孔結(jié)構(gòu)滲透率模型，如Carman-Kozeny 模型、Rumpf-Gupte 模型等均未考慮氣體稀薄效應(yīng)的影響。

Klinkenberg[6]提出了多孔結(jié)構(gòu)中稀薄氣體表觀滲透率的表達(dá)式Ka=K∞(1+n/p)。其中，K∞為材料的固有滲透率，只與多孔材料本身的幾何特征有關(guān)；n為修正因子，取決于孔隙尺度的流動(dòng)形態(tài)，也可以表示為Kn的函數(shù)。在過去的幾十年中，國內(nèi)外學(xué)者在上述表達(dá)式的基礎(chǔ)上提出了更多形式的表觀滲透率關(guān)聯(lián)式[7-13]，以適用于具有復(fù)雜孔隙形式的多孔結(jié)構(gòu)以及高壓等復(fù)雜條件。然而，現(xiàn)有稀薄氣體流動(dòng)滲透率的關(guān)聯(lián)式大多為半解析型或經(jīng)驗(yàn)型，其模型參數(shù)通過數(shù)值模擬或?qū)嶒?yàn)數(shù)據(jù)擬合的方式獲得。此外，由于在不同應(yīng)用中的多孔結(jié)構(gòu)形式多樣、孔隙的幾何拓?fù)鋸?fù)雜，不同關(guān)聯(lián)式通常只適用于特定形式的多孔結(jié)構(gòu)。隨著微納制造技術(shù)的發(fā)展，使得具有有序性孔隙的微納多孔結(jié)構(gòu)得以工程化應(yīng)用[14-19]。由于其孔隙的三維拓?fù)浣Y(jié)構(gòu)具有精確的數(shù)學(xué)描述，為從孔隙尺度流動(dòng)機(jī)理層面上獲得表觀滲透率的解析型模型提供了可能。

本文從多孔結(jié)構(gòu)孔隙尺度的流動(dòng)機(jī)理出發(fā)，確定了多孔結(jié)構(gòu)固有滲透率、孔隙率、彎曲度、收縮-擴(kuò)張因子和有效孔隙尺寸之間的定量關(guān)系。并以此為基礎(chǔ)，定義一種新的多孔結(jié)構(gòu)有效孔隙尺寸。該有效孔隙尺寸可以完全利用多孔結(jié)構(gòu)的幾何參數(shù)進(jìn)行確定，不依賴于仿真或?qū)嶒?yàn)數(shù)據(jù)的輸入。通過將該有效孔隙尺寸與稀薄氣體管道流動(dòng)模型相結(jié)合，推導(dǎo)出一種具有通用性的稀薄氣體表觀滲透率模型?；诟呔鹊闹苯幽MMonte Carlo 方法（direct simulation Monte Carlo,DSMC）對(duì)不同工況下的流動(dòng)狀態(tài)進(jìn)行數(shù)值模擬分析，對(duì)所提模型的準(zhǔn)確性進(jìn)行驗(yàn)證。

1 理論分析

1.1 多孔結(jié)構(gòu)固有滲透率表達(dá)式

固有滲透率，也稱為絕對(duì)滲透率，反映了多孔材料本身的屬性，其值為滿足Darcy 定律條件時(shí)的單相流動(dòng)滲透率。忽略多孔結(jié)構(gòu)中的閉合孔隙，當(dāng)孔隙空間Ω中充滿不可壓縮的低Reynolds數(shù)（Re＜1）牛頓流體時(shí)，孔隙中的穩(wěn)態(tài)流動(dòng)可由Stokes 方程和連續(xù)性方程控制：

其值越小代表孔隙的彎曲程度越大（也有文獻(xiàn)中采用相反比值的定義）。

圖1 理想多孔介質(zhì)中的微元流管Fig.1 Micro-element flow tube in idealized porous media

結(jié)合式（7）、式（9）、式（12）、式（14），可得，在單條流線上滿足如下關(guān)系式：

由于K∞、C、τ和ε均只與多孔結(jié)構(gòu)孔隙的幾何信息有關(guān)，按照式（18）定義的Rc也只取決于多孔結(jié)構(gòu)的幾何特性。

1.2 多孔結(jié)構(gòu)中的稀薄氣體滲透率模型

氣體的稀薄程度可以用Knudsen 數(shù)來表示，對(duì)于平直圓管流動(dòng)：

對(duì)于實(shí)際多孔材料，由于孔隙幾何特征的復(fù)雜性，特征長(zhǎng)度不能通過直徑或半徑來直接獲得。由于從本質(zhì)上來講，Knudsen 數(shù)的分母中的特征長(zhǎng)度為氣體分子與固體分子間的平均碰撞距離，可按照式（19）的定義對(duì)其進(jìn)行確定。因此，可將式（20）、式（22）、式（23）進(jìn)一步寫為

因此，根據(jù)式(21)、式(24)～式(26)，可以在孔隙幾何結(jié)構(gòu)和氣體的狀態(tài)參數(shù)已知的情況下，對(duì)多孔結(jié)構(gòu)中稀薄氣體的表觀滲透率進(jìn)行預(yù)測(cè)計(jì)算。

2 DSMC模擬驗(yàn)證

2.1 數(shù)學(xué)模型

為了驗(yàn)證本文所提出的表觀滲透率模型的準(zhǔn)確性，基于直接模擬Monte Carlo方法（DSMC）開展了不同條件下的多孔結(jié)構(gòu)稀薄氣體流動(dòng)特性仿真。DSMC是一種基于粒子運(yùn)動(dòng)的模擬方法，并收斂于Boltzmann方程[25]，因此可以對(duì)稀薄狀態(tài)下的氣體流動(dòng)過程進(jìn)行精確模擬。DSMC的模擬精度已通過不同Kn區(qū)間的實(shí)驗(yàn)數(shù)據(jù)所驗(yàn)證[26]。其在特定問題上的模擬精度高于格子Boltzmann方法，效率高于分子動(dòng)力學(xué)模擬。本文的DSMC代碼已在基于C++的開源平臺(tái)OpenFOAM上實(shí)現(xiàn)（dsmcFoam）[27]。流體區(qū)域中的計(jì)算網(wǎng)格尺寸設(shè)置為氣體最小平均自由程的1/3，計(jì)算時(shí)間步長(zhǎng)設(shè)置為分子的平均碰撞時(shí)間的1/5。為了減少統(tǒng)計(jì)誤差和提高模擬精度，每個(gè)計(jì)算單元內(nèi)設(shè)置不少于10個(gè)粒子[28]。采用無時(shí)間計(jì)數(shù)器（no time counter,NTC）方法選擇碰撞對(duì)，采用可變硬球模型（VHS）和Larsen-Borgnkke內(nèi)能重分布模型模擬粒子間的碰撞行為。DSMC模擬的典型步驟包括：（1）將計(jì)算粒子索引到計(jì)算單元；（2）計(jì)算追蹤粒子的運(yùn)動(dòng)；（3）進(jìn)行碰撞并計(jì)算碰撞后的狀態(tài)。重復(fù)這些流程直到宏觀變量的統(tǒng)計(jì)誤差足夠小。最后通過對(duì)每個(gè)單元的粒子特性進(jìn)行統(tǒng)計(jì)采樣和平均化處理，得到壓力、密度、速度、溫度等宏觀參數(shù)[29]。本文所采用的DSMC方法已在作者前期的研究中得到驗(yàn)證[29]。

2.2 物理模型與流動(dòng)特性

在本文的驗(yàn)證性模擬中，采用了三種形式的孔隙結(jié)構(gòu)。如圖2 所示，組成多孔結(jié)構(gòu)的固體微元分別設(shè)置為三維立方體[圖2（a）]、二維圓形[圖2（b）]和二維矩形[圖2（c）]。沿x方向設(shè)置壓力入口和壓力出口，y（z）方向采用周期性邊界條件。氣體入口溫度和固體微元表面溫度均恒定為273 K。固體微元的尺寸區(qū)間范圍為20～400 nm。通過改變微元間距，調(diào)節(jié)孔隙率范圍為0.12≤ε≤0.90。分別以氬氣（Ar）、氦氣（He）、氮?dú)猓∟2）和甲烷（CH4）為流動(dòng)工質(zhì)，其物性如表1 所示。此外，圖2 也給出了在Knp=1 時(shí)不同形式的多孔結(jié)構(gòu)孔隙中的速度分布云圖，不同結(jié)構(gòu)內(nèi)的氣體速度沿流動(dòng)方向均逐漸增大，這是因?yàn)殡S著氣流向出口發(fā)展，局部壓力減小，導(dǎo)致氣體密度也逐漸減小。

表1 氣體工質(zhì)物性［23］Table 1 Physical properties of gas species［23］

圖2 多孔結(jié)構(gòu)的物理模型及在Knp=1條件下的典型流場(chǎng)分布Fig.2 Physical modes of the porous structures and flow fields at Knp=1

2.3 多孔結(jié)構(gòu)中的稀薄氣體滲透率特性

圖3（a）為在固體微元為矩形[圖2（c）]，D1=D2=100 nm，孔隙率ε=0.75，入口壓力0.25 MPa，出口壓力0.05 MPa 條件下，采用不同氣體工質(zhì)時(shí)多孔結(jié)構(gòu)中心線（y=0）上的速度分布。在該邊界條件下，氦氣分子的平均自由程最大，因此氣體的稀薄性最強(qiáng)，相應(yīng)的Kn也最高。由于該中心線交替穿過固體和流體區(qū)域，所有工質(zhì)對(duì)應(yīng)的速度曲線都呈現(xiàn)出波動(dòng)形式；對(duì)于不同氣體種類，沿流程的速度波動(dòng)幅值：氦氣＞甲烷＞氮?dú)猓練鍤?。隨著氣體壓力從入口到出口的下降，氣體密度沿著流動(dòng)方向逐漸減小，速度則逐漸增大。對(duì)于不同的氣體工質(zhì)，在相同邊界條件下的體積流量相對(duì)大小為：Q(氦氣)＞Q(甲烷)＞Q(氮?dú)?＞Q(氬氣)。然而，由于氣體密度的差異，它們對(duì)應(yīng)質(zhì)量流量的相對(duì)大小則相反。在圖3（a）的條件下，氦氣、甲烷、氮?dú)夂蜌鍤庠趩挝婚L(zhǎng)度入口下的質(zhì)量流量分別為5.53、11.1、14.0和16.5 kg/(m·s)。圖3（b）為以氮?dú)鉃楣べ|(zhì)并將圖3（a）中的固體微元尺寸分別設(shè)定為20、100 和400 nm 時(shí)，多孔結(jié)構(gòu)中心線（y=0）上的速度分布?？梢钥闯?，不同單元尺度下的速度分布是相似的，且在微元的尾流區(qū)域未出現(xiàn)回流現(xiàn)象，這是因?yàn)閿U(kuò)散作用相比對(duì)流作用更占主導(dǎo)地位。圖4分別給出了與上述相同邊界條件下微元尺寸為20 nm 時(shí)的氬氣、氮?dú)狻⒓淄楹秃庠趩卧嚵兄兴俣却笮〖傲骶€的分布；可以看出，在相同條件下氦氣的最大流速遠(yuǎn)高于氬氣和氮?dú)?，這是因?yàn)楹饩哂懈蟮钠骄肿幼杂沙蹋瑥亩鴮?dǎo)致其Knudsen 數(shù)比其他分子高得多，同時(shí)也證明了氣體種類對(duì)多孔介質(zhì)中滲流特性的重要影響。

圖3 多孔結(jié)構(gòu)中心線（y=0）上的速度分布Fig.3 Distributions of velocity magnititudes at y=0

圖4 D=20 nm時(shí)氬氣、氮?dú)狻⒓淄楹秃獾牧魉俸土骶€分布Fig.4 Velocity profiles and streamlines for argon,nitrogen,methane,and helium flows at D=20 nm

實(shí)際上，上述流動(dòng)規(guī)律均可通過稀薄氣體表觀滲透率模型來進(jìn)行定量描述。根據(jù)Darcy 定律，在速度場(chǎng)已知的情況下，多孔介質(zhì)中可壓縮流體的表觀滲透率（Ka）可以用式(27)計(jì)算：

式中，uˉx,out為多孔介質(zhì)出口處的氣體宏觀平均速度。本文對(duì)不同孔隙形式（方形、圓形、矩形、立方體）、孔隙率（0.17～0.90）、氣體成分（氦氣、甲烷、氮?dú)狻鍤猓┮约安煌M(jìn)出口壓力下50 組不同工況的氣體流動(dòng)狀態(tài)進(jìn)行DSMC 模擬，并利用式（27）進(jìn)行表觀滲透率計(jì)算。通過圖5 可以發(fā)現(xiàn)，如果以式（24）所表示的Knp為橫軸，以相應(yīng)的Ka/K∞為縱軸，在整個(gè)Knp=0.01～10范圍內(nèi)，所有數(shù)據(jù)點(diǎn)均可近似分布于同一條曲線上。該曲線恰好與本文所提出的理論滲透率模型一致。模擬結(jié)果與理論模型的平均相對(duì)誤差小于10%。

圖5 不同Knp下的表觀滲透率與固有滲透率比值Fig.5 Variation of Ka/K∞with different Knp

3 結(jié) 論

多孔結(jié)構(gòu)中稀薄氣體流動(dòng)的表觀滲透率是航空航天、能源化工等領(lǐng)域相關(guān)應(yīng)用中的重要參數(shù)。本文從多孔結(jié)構(gòu)孔隙的幾何拓?fù)涮匦猿霭l(fā)，將固有滲透率、孔隙率、孔隙彎曲度和孔隙收縮-擴(kuò)張因子作為基本參數(shù)定義了一種新的多孔結(jié)構(gòu)有效孔隙尺寸。將該有效孔隙尺寸與現(xiàn)有的稀薄氣體管道流動(dòng)模型相結(jié)合，理論推導(dǎo)出了一種具有通用性的稀薄氣體滲透率模型。通過該模型，可以在孔隙三維幾何結(jié)構(gòu)和氣體的狀態(tài)參數(shù)已知的情況下對(duì)多孔結(jié)構(gòu)中稀薄氣體的表觀滲透率進(jìn)行預(yù)測(cè)計(jì)算。

為了驗(yàn)證本文所提滲透率模型的準(zhǔn)確性，在開源平臺(tái)OpenFOAM 上開展了直接Monte Carlo（DSMC）模擬。對(duì)不同孔隙微元形式（方形、圓形、矩形、立方體）、孔隙率（0.17～0.90）、氣體成分（氦氣、甲烷、氮?dú)?、氬氣）以及不同進(jìn)出口壓力下共50組不同工況的氣體流動(dòng)狀態(tài)進(jìn)行DSMC 模擬，并基于模擬結(jié)果計(jì)算出氣體的表觀滲透率。氣體成分對(duì)流動(dòng)狀態(tài)的影響顯著，在相同的壓力邊界條件下不同氣體的體積流量的相對(duì)大小為：氦氣＞甲烷＞氮?dú)猓練鍤?。而由于氣體密度的差異，它們質(zhì)量流量的大小順序則相反。在整個(gè)Knp=0.01～10 范圍內(nèi)，不同氣體、不同孔隙結(jié)構(gòu)以及不同壓力條件下的模擬結(jié)果與滲透率模型預(yù)測(cè)結(jié)果間的平均相對(duì)誤差小于10%。

因此，本文所提出的解析型模型可以充分描述多孔結(jié)構(gòu)中氣體流動(dòng)的Klinkenberg 效應(yīng)，且不依賴于仿真或?qū)嶒?yàn)數(shù)據(jù)作為模型輸入條件，為分析和評(píng)估微納尺度或低氣壓狀態(tài)下有序性多孔結(jié)構(gòu)的氣體滲透率特性提供了一種新的解決思路。

符 號(hào) 說 明

B——水力傳導(dǎo)系數(shù)

C——收縮-擴(kuò)張因子

D——模型中微單元的特征尺寸，m

k——微元滲透率，m2

L——多孔區(qū)域長(zhǎng)度，m

R——半徑，m

Rc——有效孔隙半徑，m

u——流體速度矢量，m/s

V——多孔介質(zhì)總體積，m3

ρ——?dú)怏w密度，kg/m3

τ——彎曲度

ω——溫度系數(shù)

下角標(biāo)

in——入口

out——出口

s——流線