微通道內(nèi)氣體流動(dòng)的DSMC方法

張文雯， 王學(xué)德

(南京理工大學(xué)能源與動(dòng)力工程學(xué)院， 江蘇南京 210094)

引 言

隨著微機(jī)電系統(tǒng)以及納米技術(shù)的飛速發(fā)展， 各式各樣的微機(jī)械設(shè)備與元件已應(yīng)用于自動(dòng)控制、 醫(yī)療儀器、 電子冷卻和航空航天等諸多領(lǐng)域中[1-4]， 其中涉及了許多有關(guān)微通道的流體流動(dòng)， 使得微通道內(nèi)部氣體的流動(dòng)特性受到研究人員的廣泛關(guān)注。近年來(lái)， 研究人員基于不同形狀的直微通道氣體流動(dòng)做了很多工作[5-7]， 而在微通道的實(shí)際應(yīng)用中普遍存在著內(nèi)部流道橫截面變化的復(fù)雜情況， 關(guān)于此類情況的相關(guān)研究較少?？卓谀Ｐ妥鳛槲C(jī)電系統(tǒng)中典型的流道截面變化結(jié)構(gòu)， 其流道截面出現(xiàn)收縮/擴(kuò)張現(xiàn)象， 導(dǎo)致流場(chǎng)產(chǎn)生壓降并伴隨流動(dòng)分離現(xiàn)象[8]。因此， 揭示氣體在具有孔口結(jié)構(gòu)微通道下的流動(dòng)機(jī)理對(duì)于微機(jī)電系統(tǒng)實(shí)際裝置的設(shè)計(jì)、 制造和運(yùn)行具有重要的意義。

本文采用直接模擬Monte Carlo(direct simul-ation Monte Carlo， DSMC)方法對(duì)微尺度氣體流動(dòng)進(jìn)行模擬， 網(wǎng)格尺寸和時(shí)間步長(zhǎng)等模擬參數(shù)會(huì)直接影響DSMC模擬結(jié)果的準(zhǔn)確性。在以往的微流體研究中， 網(wǎng)格尺寸與時(shí)間步長(zhǎng)研究較少， 國(guó)內(nèi)外研究人員通常做法是以Bird[9]提出的網(wǎng)格最大尺寸約為當(dāng)?shù)胤肿悠骄杂沙痰?/3作為劃分依據(jù)。但微通道流動(dòng)平均自由程較小， 根據(jù)此原則會(huì)導(dǎo)致網(wǎng)格數(shù)量過(guò)多以及模擬分子數(shù)急劇增加， 造成極大的運(yùn)算成本。

因此， 為解決以上問(wèn)題， 本文首先針對(duì)壓力梯度驅(qū)動(dòng)的微通道流動(dòng)發(fā)展了基于壓力邊界條件的DSMC方法。隨后定義CSL數(shù)和DCFL數(shù)兩個(gè)無(wú)量綱參數(shù)作為網(wǎng)格尺寸和時(shí)間步長(zhǎng)的約束條件， 通過(guò)對(duì)微通道流的模擬， 總結(jié)了在DSMC方法中微尺度條件下網(wǎng)格尺寸和時(shí)間步長(zhǎng)的一般選取原則。并在該選取原則的基礎(chǔ)上， 探究了單/雙孔口結(jié)構(gòu)微通道內(nèi)氣體流動(dòng)結(jié)構(gòu)， 并通過(guò)改變?nèi)肟趬毫Ψ治隽鞯澜孛孀兓瘜?duì)氣體流動(dòng)的影響。

1 微通道DSMC壓力邊界條件

微通道中的氣體流動(dòng)往往是由壓力梯度驅(qū)動(dòng)的低速流動(dòng)， 此時(shí)， 分子數(shù)密度、 來(lái)流速度等參數(shù)難以獲得， 壓力和溫度是邊界上唯一已知的物理量[10]， 因此， 采用DSMC壓力邊界條件。已知進(jìn)口壓力pin, 溫度Tin以及出口壓力pe, 流入或流出微通道的粒子數(shù)通量為

(1)

式中，n表示分子數(shù)密度，θ為速度方向與邊界法向的夾角， erf表示的是誤差函數(shù)，s為分子速率比，β表示分子最可幾速度cmp的倒數(shù)， 即

式中，T和u分別表示邊界溫度以及來(lái)流平均速度分量,m為分子質(zhì)量，k為Boltzmann常數(shù)。

采用基于特征理論的壓力邊界條件[11]進(jìn)行計(jì)算。對(duì)于進(jìn)出口邊界單元(in表示入口邊界， e表示出口邊界)， 每個(gè)單元j的宏觀速度u, 出口密度ρe以及出口溫度Te由下式與式(1)迭代求解

ρj=njm，pj=njkTj

Tj=(3Ttr+ζTrot)/(3+ζ)

式中，Nj為網(wǎng)格內(nèi)的采樣粒子數(shù)，Ttr為平動(dòng)溫度，Trot為轉(zhuǎn)動(dòng)溫度，ζ為轉(zhuǎn)動(dòng)自由度。

2 微通道DSMC方法驗(yàn)證

采用經(jīng)典微Poiseuille流驗(yàn)證微尺度DSMC方法的有效性。二維微通道幾何模型如圖 1所示， 計(jì)算條件取自文獻(xiàn)[12]， 見表 1。其中，H為 0.4 μm，L為2 μm。左右邊界為自由來(lái)流， 模擬氮?dú)庠谖⑼ǖ纼?nèi)流動(dòng)特性， 上下邊界為采用漫反射模型的固體壁面。無(wú)時(shí)間計(jì)時(shí)器技術(shù)(no time counter， NTC)和可變硬球(variable hard sphere， VHS)模型被用來(lái)模擬粒子之間的碰撞。

圖1 微通道基本結(jié)構(gòu)Fig. 1 Basic structure of the microchannel

表1 微通道計(jì)算參數(shù)Table 1 Numerical parameters for microchannel

圖2給出了兩種壓力比的中心線壓力分布與壁面滑移速度us變化的比較結(jié)果。圖2(a)首先將本文使用特征壓力邊界條件的DSMC計(jì)算結(jié)果與Arkilic等[13]通過(guò)1階解析壓力分布作比較， 結(jié)果表明本文計(jì)算結(jié)果與解析解具有較高的一致性， 驗(yàn)證了本文程序的準(zhǔn)確性。同時(shí)將壓力與壁面滑移速度與Fang[12]使用隱式邊界條件的結(jié)果進(jìn)行對(duì)比， 兩種分布量的模擬結(jié)果均符合較好， 進(jìn)一步證明本文方法的有效性。

(a) Centerline pressure distribution

(b) Wall slip velocity distribution圖2 微通道模擬值與文獻(xiàn)值的對(duì)比Fig. 2 Comparison between simulation and reference results

3 網(wǎng)格尺寸與時(shí)間步長(zhǎng)約束條件

DSMC方法的本質(zhì)是在一段時(shí)間步長(zhǎng)Δt內(nèi)， 使得分子運(yùn)動(dòng)和碰撞能夠解耦。但為了確保能達(dá)到正確解耦的目的， 一般需要保證以下兩個(gè)條件[14]: (1) 網(wǎng)格尺寸小于分子平均自由程λ, 一般取為平均自由程的1/3； (2) 時(shí)間步長(zhǎng)Δt小于平均碰撞時(shí)間τ。下面研究其在微通道內(nèi)的適用性和可擴(kuò)展性。為此， 首先定義一個(gè)無(wú)量綱參數(shù)CSL表示網(wǎng)格尺寸與平均自由程之間的關(guān)系。同時(shí)借用計(jì)算流體力學(xué)中CFL數(shù)來(lái)定義一個(gè)新的參數(shù)DCFL數(shù)限制時(shí)間步長(zhǎng)， 表示如下

以壓力比為2.5的微通道來(lái)驗(yàn)證網(wǎng)格尺寸和時(shí)間步長(zhǎng)對(duì)計(jì)算結(jié)果的影響。計(jì)算了6種不同CSL數(shù)和DCFL數(shù)組合工況， 如表 2所示。其中， 工況3為第2節(jié)驗(yàn)證算例中壓力比為2.5的計(jì)算參數(shù)， 以該參數(shù)下的計(jì)算結(jié)果作為比較的基準(zhǔn)。由表中數(shù)據(jù)可知， 除了工況6的DCFL數(shù)等于1外， 其余情況DCFL數(shù)均小于1。

表2 網(wǎng)格尺寸與時(shí)間步長(zhǎng)約束參數(shù)Table 2 Constraint parameters for grid size and time step

圖3給出了各種工況下中心線壓力分布與解析壓力值對(duì)比情況， 可以看出除了工況6之外， 其余工況的壓力分布均與解析值符合良好。為了進(jìn)一步分析時(shí)間步長(zhǎng)與網(wǎng)格參數(shù)對(duì)其他物理量的影響， 圖 4， 5分別給出了中心截面處速度分布與沿通道中心線溫度分布， 可以看出DCFL數(shù)小于1的工況1， 2， 4， 5的速度與溫度分布結(jié)果與經(jīng)過(guò)驗(yàn)證的工況3的結(jié)果趨于一致， 在工況6情況下， 中心截面速度與溫度的計(jì)算結(jié)果誤差均出現(xiàn)較大的偏離。綜合比較圖3， 4， 5的計(jì)算結(jié)果， 可以得到在選取網(wǎng)格尺寸和時(shí)間步長(zhǎng)時(shí)， 首先應(yīng)當(dāng)保證DCFL數(shù)的取值小于1作為時(shí)間步長(zhǎng)的約束條件。

圖3 中心線壓力分布對(duì)比Fig. 3 Comparison of centerline pressure distribution

在滿足DCFL數(shù)小于1的情況下， 進(jìn)一步分析CSL數(shù)對(duì)中心截面速度與溫度計(jì)算結(jié)果的影響。從圖 4看出工況1中心截面的流向速度明顯偏小， 而工況2， 4和5均能得到與工況工況3較為一致的模擬結(jié)果， 從而表明在CSL數(shù)小于1.865時(shí)也可以得到較為準(zhǔn)確的計(jì)算結(jié)果。從圖 5(b)的局部放大圖中可以看出工況5的中心線溫度較工況 3 明顯偏低。由此可見， 在選擇網(wǎng)格尺寸時(shí)不需要過(guò)分追求小CSL數(shù)， 當(dāng)網(wǎng)格尺寸很小時(shí)， 是同一個(gè)網(wǎng)格內(nèi)模擬分子碰撞的次數(shù)減少， 造成一定統(tǒng)計(jì)誤差， 并且由于網(wǎng)格數(shù)量變多， 流場(chǎng)收斂需要的時(shí)間變長(zhǎng)， 增加了運(yùn)算成本。因此， 在確保DCFL數(shù)小于1的前提下， 選擇CSL數(shù)接近于1可作為本文網(wǎng)格尺寸選取的原則。

圖4 中心截面處流向速度分布對(duì)比Fig. 4 Comparison of velocity distribution at the channel midspan

(a) Centerline temperature distribution

(b) Enlarged red frame圖5 中心線溫度分布對(duì)比Fig. 5 Comparison of centerline temperature distribution

4 單/雙孔口微通道流動(dòng)特性研究

4.1 幾何條件設(shè)置與物理模型

研究模型為具有單/雙孔口結(jié)構(gòu)的微通道， 如圖 6所示， 兩種不同孔口結(jié)構(gòu)的尺寸參數(shù)見表 3， 壁面采用漫反射模型， 模擬氣體為氮?dú)猓?不同入口壓力條件下[9]的來(lái)流條件見表 4， 進(jìn)出口流態(tài)均處于滑移流狀態(tài)。

(a) Single orifice model

(b) Double orifice model圖6 微通道幾何形狀示意圖Fig. 6 Schematic of the microchannel geometry

表3 不同形狀模型的幾何條件Table 3 Geometric conditions of different shapes

表4 孔口流通道參數(shù)Table 4 Orifice flow channel parameters

4.2 入口壓力對(duì)微通道流場(chǎng)的影響

圖 7分別給出了不同入口壓力下， 單/雙孔口模型通道內(nèi)的流線與壓力等值線分布圖。在單孔口模型中， 通道上游處流線沿拐角進(jìn)入孔， 在孔口內(nèi)部流線與通道平行， 主流經(jīng)孔口加速， 在孔口出口形成一股收縮流， 由于管道突然擴(kuò)張， 流動(dòng)減速以充滿整個(gè)通道?？卓诔隹谟捎诹鲃?dòng)速度降低， 壓力對(duì)應(yīng)升高， 在孔口下游的拐角處形成反向壓力梯度， 從而出現(xiàn)回流區(qū)。在雙孔口模型中， 具有與單孔口模型相似的流動(dòng)結(jié)構(gòu)， 且每一個(gè)孔口后均存在回流區(qū)。此外， 當(dāng)入口壓力增大之后， 由于出口壓力始終保持不變， 逆壓區(qū)的壓差增大， 回流區(qū)尺寸隨之變大。在雙孔口模型中， 重點(diǎn)關(guān)注兩孔口之間的過(guò)渡區(qū)域， 入口壓力為1.5×105Pa時(shí)， 在過(guò)渡區(qū)拐角處出現(xiàn)了兩處分離， 孔口2的上游也出現(xiàn)一個(gè)很小的分離區(qū)。這是由于流體在過(guò)渡區(qū)流程較短， 流體經(jīng)過(guò)孔口壓縮加速， 無(wú)法及時(shí)調(diào)整流向， 出現(xiàn)分離區(qū)。而隨著入口壓力的增加， 分離區(qū)的尺寸增大， 兩處分離逐漸融合， 形成一個(gè)大的分離渦。

(a) pin=1.5×105 Pa

(b) pin=2.0×105 Pa

(c) pin=2.5×105 Pa

(d) pin=3.0×105 Pa圖7 流線與壓力等值線云圖Fig. 7 Streamlines and pressure contours

圖8給出了單/雙孔口的中心線壓力分布曲線， 從中可以看出主要的壓力變化均集中于孔口位置， 兩種模型均在孔口位置出現(xiàn)了較大的壓降， 其中， 雙孔口模型經(jīng)歷了兩次壓降， 第2次壓降比第1次大很多。這是由于當(dāng)流道收縮使得流體加速時(shí)， 流體的靜壓能轉(zhuǎn)換為動(dòng)能， 孔口上游的壓力快速下降。此外， 最小壓力不是在截面最小的孔口處， 而是出現(xiàn)在孔口下游， 隨著入口壓力的增大， 最小壓力的位置也隨之后移， 并且孔口后的壓力甚至?xí)∮诔隹趬毫?。在最小壓力位置的下游?隨著動(dòng)能轉(zhuǎn)換回靜壓能， 流動(dòng)開始減速， 壓力逐漸增加至局部最大值， 即出口壓力。

(a) Single orifice model

(b) Double orifice model圖8 中心線壓力分布Fig. 8 Centerline pressure distribution

圖 9， 10分別為不同入口壓力下通道中心線的速度分布和溫度分布， 可以看出速度與溫度的變化規(guī)律相反。孔口結(jié)構(gòu)類似于一個(gè)收縮/擴(kuò)張的裝置， 因此氣體在經(jīng)過(guò)孔口時(shí)被壓縮， 導(dǎo)致孔口處的速度增加， 其后， 由于通道擴(kuò)張， 速度減小， 流場(chǎng)中溫度變化與速度變化相反。在單孔結(jié)構(gòu)中， 可以發(fā)現(xiàn)氣流經(jīng)過(guò)孔口時(shí)速度上升非?？?， 由于慣性作用， 在離開孔口后速度達(dá)到最大， 再緩慢降低， 且流場(chǎng)的最大速度隨著入口壓力的增大而增大。在雙孔模型中， 發(fā)現(xiàn)孔口1的速度變化幅度較小， 一方面孔口1的位置處于流場(chǎng)的前半部分， 本身速度就較低； 另一方面， 氣體速度與壓降大小有關(guān)。當(dāng)入口壓力為1.5×105Pa時(shí)， 由于兩次壓降大小幾乎相同， 因此經(jīng)過(guò)兩個(gè)孔口速度增加量相似； 而當(dāng)入口壓力逐漸增加， 經(jīng)孔口2的壓降比孔口1要大得多， 因此， 孔口2的速度變化更劇烈。對(duì)比同一入口壓力下兩個(gè)模型的速度大小， 單孔口的最大速度要大于雙孔口模型， 由于雙孔口模型經(jīng)兩次壓降， 在相同壓差情況下， 雙孔口的單次壓降小， 使得流體速度每次增加的幅度也小。另外， 由于雙孔口流場(chǎng)產(chǎn)生兩處分離區(qū)， 流場(chǎng)擾動(dòng)作用更加明顯， 對(duì)流體的速度也產(chǎn)生了一定影響。

(a) Single orifice model

(b) Double orifice model圖9 中心線速度分布Fig. 9 Centerline velocity distribution

(a) Single orifice model

(b) Double orifice model圖10 中心線溫度分布Fig. 10 Centerline temperature distribution

5 結(jié)論

本文基于壓力邊界條件的DSMC方法， 研究了不同網(wǎng)格尺寸和時(shí)間步長(zhǎng)對(duì)直微通道模型結(jié)果精度的影響， 得到了在微尺度下的網(wǎng)格尺寸和時(shí)間步長(zhǎng)選取的一般準(zhǔn)則。在此基礎(chǔ)上， 研究了不同入口壓力條件對(duì)單/雙孔口結(jié)構(gòu)微通道內(nèi)流動(dòng)過(guò)程的影響。主要得出以下結(jié)論:

(1) 通過(guò)CSL數(shù)和DCFL數(shù)兩個(gè)無(wú)量綱參數(shù)給出了網(wǎng)格尺寸和時(shí)間步長(zhǎng)的約束條件， 研究結(jié)果表面在滿足DCFL數(shù)嚴(yán)格小于1的情況下， 選擇CSL數(shù)接近或小于1， 可以作為網(wǎng)格尺寸與時(shí)間步長(zhǎng)的選取原則。

(2)由于壓降的存在使得孔口后形成低壓區(qū)， 因此氣體經(jīng)過(guò)孔口后會(huì)在拐角處形成分離區(qū)， 并且分離區(qū)大小隨入口壓力的增大而增大。雙孔口模型中， 兩孔口之間的過(guò)渡區(qū)域流程較短， 但也出現(xiàn)了分離現(xiàn)象， 并隨入口壓力增大， 過(guò)渡區(qū)兩處分離渦逐漸融合。在相同壓差情況下， 通過(guò)雙孔口模型流速會(huì)低于單個(gè)孔口通道的流速。

(3)在孔口模型中， 氣流在通過(guò)孔口時(shí)產(chǎn)生較大壓降。經(jīng)孔口處的氣流速度增加， 并隨入口壓力增大而增大。最大速度的位置由于慣性作用逐漸后移， 出現(xiàn)在孔口下游， 其后， 速度逐漸降低。溫度與速度變化相反。