微納多孔結構中稀薄氣體流動滲透率的解析型預測模型

作者：楊光程鑫王崢王曄張良俊吳靜怡來源：《化工學報》日期：2022-10-28人氣：2283

多孔結構中的稀薄氣體流動機理在石油化工、能源利用、航空航天等領域具有廣闊的應用前景[1-3]，如航天推進劑的增壓輸送、非常規(guī)油氣開采、質(zhì)子交換膜燃料電池、CO2封存、核廢料處理等。近年來，3D打印、MEMS技術的發(fā)展也推動著對微納尺度氣體流動特性的深入研究[4-5]。當多孔結構的孔隙尺度足夠小，或者氣體工質(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效應。而目前常用的多孔結構滲透率模型，如Carman-Kozeny模型、Rumpf-Gupte模型等均未考慮氣體稀薄效應的影響。

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

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

1 理論分析

1.1 多孔結構固有滲透率表達式

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

$μ ?^{2} u = ? p$ （1） $? ? u = 0$ （2）

假設多孔結構孔隙中的任何一點都隸屬于某條連通了多孔區(qū)域入口和出口的流線S。令Ω表示所有流線所占據(jù)空間的集合，則孔隙率可定義為Ω和多孔介質(zhì)總體積V的比值，即

$ε = \frac{Ω}{V}$ （3）

多孔介質(zhì)中牛頓流體的緩慢流動特性可以由Darcy定律描述

$K_{\infty} = - \frac{Q μ L}{A Δ p}$ （4）

式中，Q為流體的流量；A為多孔介質(zhì)截面積；L為多孔區(qū)域長度； $Δ$ p為流過多孔區(qū)域的前后壓差。Darcy定律也可寫成沿流體區(qū)域的積分形式：

$K_{\infty} = - \frac{Q μ L}{A Δ p} = ε \frac{1}{Ω} \int_{Ω}^{} - μ ? p ? u {(\frac{L}{Δ p})}^{2} d V$ （5）

定義微元的滲透率因子 $k$ ：

$k = \frac{K_{\infty}}{ε} = \frac{1}{Ω} \int_{Ω}^{} - μ ? p ? u {(\frac{L}{Δ p})}^{2} d V$ （6）

按照前面的假設，將孔隙空間Ω離散為無數(shù)個流線S的集合，并定義單條流線上的滲透率因子：

$k (S) = - \frac{μ L^{2}}{Δ p}$ （7）

結合式（6）、式（7），可以得出[20]

$k = \frac{1}{Ω} \int_{?}^{} k (S) d Q_{S} = \frac{k (S)}{Ω} \int_{?}^{} d Q_{S} = \frac{Q}{Ω} k (S)$ （8）

式中，?為兩個標量函數(shù)x和y構成的二維面域。由于多孔結構孔隙中的流線特征與孔隙彎曲度和收縮-擴張?zhí)匦杂嘘P。下面分別對其進行定義。

孔隙的彎曲度可定義為多孔介質(zhì)長度L和實際流線長度Ls的比值[圖1(a)]：

$τ (S) = \frac{L}{L_{s}}$ （9）

圖1

圖1 理想多孔介質(zhì)中的微元流管

Fig.1 Micro-element flow tube in idealized porous media

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

對于變截面流動通道[圖1(b)]，通道長度為L，通道長度方向位置x處的流動截面積為 $A (x)$ ，其截面的收縮-擴張因子可以描述為：

$C = \frac{1}{L^{2}} \int_{0}^{L} A {(x)}^{2} d x \int_{0}^{L} \frac{1}{A {(x)}^{2}} d x$ （10）

根據(jù)面積與壓差的關系[21]，進一步可得

$C = \frac{1}{L^{2}} \int_{0}^{L} \frac{1}{? p (x)} d x \int_{0}^{L} ? p (x) d x$ （11）

將上式應用于流線S上，并利用沿流線的壓力梯度 $\frac{δ p}{δ S} = \frac{? p ? u}{u}$ 替換式（11）中的壓力梯度，可得到流線上的收縮-擴張因子表達式：

$C (S) = \frac{1}{L_{s}^{2}} \int_{S}^{} \frac{u}{? p ? u} d S \int_{S}^{} \frac{? p ? u}{u} d S = \frac{1}{L_{s}^{2}} Δ p \int_{S}^{} \frac{u}{? p ? u} d S$ （12）

根據(jù)水力傳導系數(shù)的定義[22]：

$B = - \frac{μ u^{2}}{? p ? u}$ （13）

則單條流線上的水力傳導系數(shù)為：

$B (S) = \int_{S}^{} B \frac{1}{u} d S = \int_{S}^{} - \frac{μ u}{? p ? u} d S$ （14）

結合式（7）、式（9）、式（12）、式（14），可得，在單條流線上滿足如下關系式：

$k (S) = \frac{B (S) τ {(S)}^{2}}{C (S)}$ （15）

將式（15）中的各變量在孔隙空間Ω內(nèi)進行體積積分，并結合式（6）可得

$K_{\infty} = \frac{τ^{2} B ε}{C}$ （16）

式（16）給出了固有滲透率、水力傳導系數(shù)、孔隙率、孔隙彎曲度和收縮-擴張因子之間的定量關系。

結合式（5）和式（16）可以得出，對于平直圓管中的流動

$K_{\infty} = - \frac{Q μ L}{A Δ p} = B$ （17）

另外，平直圓管流動可以利用Hagen-Poisseuille方程來描述

$Q = - \frac{π R^{4} Δ p}{8 μ L}$ （18）

式中，R為圓管的管道半徑。結合式（17）和式（18）可以得出管道半徑 $R = \sqrt[]{8 B}$ 。因此，通過上述類比，對于由復雜孔隙組成的多孔結構，可將表征孔隙有效尺寸的特征半徑取 $R_{c} = \sqrt[]{8 B}$ ，并結合式（16），即

$R_{c} = \sqrt[]{\frac{8 K_{\infty} C}{τ^{2} ε}}$ （19）

由于 $K_{\infty}$ 、C、 $τ$ 和 $ε$ 均只與多孔結構孔隙的幾何信息有關，按照式（18）定義的Rc也只取決于多孔結構的幾何特性。

1.2 多孔結構中的稀薄氣體滲透率模型

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

$K n = \frac{λ}{R}$ （20）

式中，λ為氣體分子平均自由程。根據(jù)Bird[23]提出的可變硬球模型（variable hard sphere, VHS），實際氣體的平均分子自由程可以表示為：

$λ = \frac{2 (5 - 2 ω) (7 - 2 ω)}{15} {(\frac{m}{2 π b T})}^{\frac{1}{2}} \frac{μ}{ρ}$ （21）

式中，μ為氣體動力黏度；ρ為氣體密度；m為氣體分子質(zhì)量；b為Boltzmann常數(shù)，取值1.38×10-23； $ω$ 為溫度系數(shù)，反映了黏度與溫度間的關系。

早在1909年，Knudsen發(fā)現(xiàn)隨著平均壓力降低氣體密度減小，在給定壓降下通過毛細管的氣體質(zhì)量流量會先減小后增大，并提出了如下公式來計算高Knudsen數(shù)下的氣體流量[24]：

$\frac{K_{a}}{K_{\infty}} = 1 + \frac{64}{3 π} \times \frac{1 + c_{1}^{K} p}{1 + c_{2}^{K} p} K n$ （22） $c_{1}^{K} p = \sqrt[]{\frac{π}{2}} \times \frac{2}{K n}, c_{2}^{K} p = \sqrt[]{\frac{π}{2}} \times \frac{2.47}{K n}$ （23）

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

$K n_{p} = \frac{λ τ \sqrt[]{ε}}{\sqrt[]{8 K_{\infty} C}}$ （24） $\frac{K_{a}}{K_{\infty}} = 1 + \frac{64}{3 π} \times \frac{1 + c_{1}^{K} p}{1 + c_{2}^{K} p} K n_{p}$ （25） $c_{1}^{K} p = \sqrt[]{\frac{π}{2}} \times \frac{2}{K n_{p}}, c_{2}^{K} p = \sqrt[]{\frac{π}{2}} \times \frac{2.47}{K n_{p}}$ （26）

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

2 DSMC模擬驗證

2.1 數(shù)學模型

為了驗證本文所提出的表觀滲透率模型的準確性，基于直接模擬Monte Carlo方法（DSMC）開展了不同條件下的多孔結構稀薄氣體流動特性仿真。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 物理模型與流動特性

在本文的驗證性模擬中，采用了三種形式的孔隙結構。如圖2所示，組成多孔結構的固體微元分別設置為三維立方體[圖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時不同形式的多孔結構孔隙中的速度分布云圖，不同結構內(nèi)的氣體速度沿流動方向均逐漸增大，這是因為隨著氣流向出口發(fā)展，局部壓力減小，導致氣體密度也逐漸減小。

圖2

圖2 多孔結構的物理模型及在Knp =1條件下的典型流場分布

Fig.2 Physical modes of the porous structures and flow fields at Knp =1

表1 氣體工質(zhì)物性［23］

Table 1 Physical properties of gas species［23］

氣體種類	分子直徑(d)/ 10-10 m	分子質(zhì)量(m)/ 10-23 kg	自由度(ξ)	溫度系數(shù)(ω)
氬氣 (Ar)	4.17	66.3	3	0.81
氦氣 (He)	2.33	6.65	3	0.66
甲烷 (CH4)	4.83	26.6	6.4	0.84
氮氣 (N2)	4.17	46.5	5	0.74

新窗口打開| 下載CSV

2.3 多孔結構中的稀薄氣體滲透率特性

圖3（a）為在固體微元為矩形[圖2（c）]，D1=D2=100 nm，孔隙率ε=0.75，入口壓力0.25 MPa，出口壓力0.05 MPa條件下，采用不同氣體工質(zhì)時多孔結構中心線（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時，多孔結構中心線（y=0）上的速度分布。可以看出，不同單元尺度下的速度分布是相似的，且在微元的尾流區(qū)域未出現(xiàn)回流現(xiàn)象，這是因為擴散作用相比對流作用更占主導地位。圖4分別給出了與上述相同邊界條件下微元尺寸為20 nm時的氬氣、氮氣、甲烷和氦氣在單元陣列中速度大小及流線的分布；可以看出，在相同條件下氦氣的最大流速遠高于氬氣和氮氣，這是因為氦氣具有更大的平均分子自由程，從而導致其Knudsen數(shù)比其他分子高得多，同時也證明了氣體種類對多孔介質(zhì)中滲流特性的重要影響。

圖3

圖3 多孔結構中心線（y=0）上的速度分布

Fig.3 Distributions of velocity magnititudes at y=0

圖4

圖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)計算：

$K_{a} = \frac{2 μ L p_{o u t} {\overset{ˉ}{u}}_{x, o u t}}{p_{i n}^{2} - p_{o u t}^{2}}$ （27）

式中， ${\overset{ˉ}{u}}_{x, o u t}$ 為多孔介質(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ù)點均可近似分布于同一條曲線上。該曲線恰好與本文所提出的理論滲透率模型一致。模擬結果與理論模型的平均相對誤差小于10%。

圖5

圖5 不同Knp 下的表觀滲透率與固有滲透率比值

Fig.5 Variation of Ka/K∞ with different Knp

3 結論

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

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

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

符號說明

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	流線

關鍵字：優(yōu)秀論文

上一篇：光熱驅(qū)動多孔氧化鈰熱化學循環(huán)解水制氫非熱質(zhì)平衡模型
下一篇：受限空間內(nèi)浮升氣泡與液體間傳質(zhì)行為實驗研究

欄目分類

熱門排行

推薦信息

期刊知識