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

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

引 言

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

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

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

1 理論分析

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

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

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

圖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ù)來表示，對于平直圓管流動：

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

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

2 DSMC模擬驗證

2.1 數(shù)學模型

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

2.2 物理模型與流動特性

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

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

圖2 多孔結(jié)構(gòu)的物理模型及在Knp=1條件下的典型流場分布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ì)時多孔結(jié)構(gòu)中心線（y=0）上的速度分布。在該邊界條件下，氦氣分子的平均自由程最大，因此氣體的稀薄性最強，相應的Kn也最高。由于該中心線交替穿過固體和流體區(qū)域，所有工質(zhì)對應的速度曲線都呈現(xiàn)出波動形式；對于不同氣體種類，沿流程的速度波動幅值：氦氣＞甲烷＞氮氣＞氬氣。隨著氣體壓力從入口到出口的下降，氣體密度沿著流動方向逐漸減小，速度則逐漸增大。對于不同的氣體工質(zhì)，在相同邊界條件下的體積流量相對大小為：Q(氦氣)＞Q(甲烷)＞Q(氮氣)＞Q(氬氣)。然而，由于氣體密度的差異，它們對應質(zhì)量流量的相對大小則相反。在圖3（a）的條件下，氦氣、甲烷、氮氣和氬氣在單位長度入口下的質(zhì)量流量分別為5.53、11.1、14.0和16.5 kg/(m·s)。圖3（b）為以氮氣為工質(zhì)并將圖3（a）中的固體微元尺寸分別設定為20、100 和400 nm 時，多孔結(jié)構(gòu)中心線（y=0）上的速度分布?？梢钥闯?，不同單元尺度下的速度分布是相似的，且在微元的尾流區(qū)域未出現(xiàn)回流現(xiàn)象，這是因為擴散作用相比對流作用更占主導地位。圖4分別給出了與上述相同邊界條件下微元尺寸為20 nm 時的氬氣、氮氣、甲烷和氦氣在單元陣列中速度大小及流線的分布；可以看出，在相同條件下氦氣的最大流速遠高于氬氣和氮氣，這是因為氦氣具有更大的平均分子自由程，從而導致其Knudsen 數(shù)比其他分子高得多，同時也證明了氣體種類對多孔介質(zhì)中滲流特性的重要影響。

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

圖4 D=20 nm時氬氣、氮氣、甲烷和氦氣的流速和流線分布Fig.4 Velocity profiles and streamlines for argon,nitrogen,methane,and helium flows at D=20 nm

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

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

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

3 結(jié) 論

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

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

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

符 號 說 明

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

C——收縮-擴張因子

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

k——微元滲透率，m2

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

R——半徑，m

Rc——有效孔隙半徑，m

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

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

ρ——氣體密度，kg/m3

τ——彎曲度

ω——溫度系數(shù)

下角標

in——入口

out——出口

s——流線