張文雯, 王學(xué)德
(南京理工大學(xué)能源與動力工程學(xué)院, 江蘇南京 210094)
隨著微機(jī)電系統(tǒng)以及納米技術(shù)的飛速發(fā)展, 各式各樣的微機(jī)械設(shè)備與元件已應(yīng)用于自動控制、 醫(yī)療儀器、 電子冷卻和航空航天等諸多領(lǐng)域中[1-4], 其中涉及了許多有關(guān)微通道的流體流動, 使得微通道內(nèi)部氣體的流動特性受到研究人員的廣泛關(guān)注。近年來, 研究人員基于不同形狀的直微通道氣體流動做了很多工作[5-7], 而在微通道的實(shí)際應(yīng)用中普遍存在著內(nèi)部流道橫截面變化的復(fù)雜情況, 關(guān)于此類情況的相關(guān)研究較少。孔口模型作為微機(jī)電系統(tǒng)中典型的流道截面變化結(jié)構(gòu), 其流道截面出現(xiàn)收縮/擴(kuò)張現(xiàn)象, 導(dǎo)致流場產(chǎn)生壓降并伴隨流動分離現(xiàn)象[8]。因此, 揭示氣體在具有孔口結(jié)構(gòu)微通道下的流動機(jī)理對于微機(jī)電系統(tǒng)實(shí)際裝置的設(shè)計、 制造和運(yùn)行具有重要的意義。
本文采用直接模擬Monte Carlo(direct simul-ation Monte Carlo, DSMC)方法對微尺度氣體流動進(jìn)行模擬, 網(wǎng)格尺寸和時間步長等模擬參數(shù)會直接影響DSMC模擬結(jié)果的準(zhǔn)確性。在以往的微流體研究中, 網(wǎng)格尺寸與時間步長研究較少, 國內(nèi)外研究人員通常做法是以Bird[9]提出的網(wǎng)格最大尺寸約為當(dāng)?shù)胤肿悠骄杂沙痰?/3作為劃分依據(jù)。但微通道流動平均自由程較小, 根據(jù)此原則會導(dǎo)致網(wǎng)格數(shù)量過多以及模擬分子數(shù)急劇增加, 造成極大的運(yùn)算成本。
因此, 為解決以上問題, 本文首先針對壓力梯度驅(qū)動的微通道流動發(fā)展了基于壓力邊界條件的DSMC方法。隨后定義CSL數(shù)和DCFL數(shù)兩個無量綱參數(shù)作為網(wǎng)格尺寸和時間步長的約束條件, 通過對微通道流的模擬, 總結(jié)了在DSMC方法中微尺度條件下網(wǎng)格尺寸和時間步長的一般選取原則。并在該選取原則的基礎(chǔ)上, 探究了單/雙孔口結(jié)構(gòu)微通道內(nèi)氣體流動結(jié)構(gòu), 并通過改變?nèi)肟趬毫Ψ治隽鞯澜孛孀兓瘜怏w流動的影響。
微通道中的氣體流動往往是由壓力梯度驅(qū)動的低速流動, 此時, 分子數(shù)密度、 來流速度等參數(shù)難以獲得, 壓力和溫度是邊界上唯一已知的物理量[10], 因此, 采用DSMC壓力邊界條件。已知進(jìn)口壓力pin, 溫度Tin以及出口壓力pe, 流入或流出微通道的粒子數(shù)通量為
(1)
式中,n表示分子數(shù)密度,θ為速度方向與邊界法向的夾角, erf表示的是誤差函數(shù),s為分子速率比,β表示分子最可幾速度cmp的倒數(shù), 即
式中,T和u分別表示邊界溫度以及來流平均速度分量,m為分子質(zhì)量,k為Boltzmann常數(shù)。
采用基于特征理論的壓力邊界條件[11]進(jìn)行計算。對于進(jìn)出口邊界單元(in表示入口邊界, e表示出口邊界), 每個單元j的宏觀速度u, 出口密度ρe以及出口溫度Te由下式與式(1)迭代求解
ρj=njm,pj=njkTj
Tj=(3Ttr+ζTrot)/(3+ζ)
式中,Nj為網(wǎng)格內(nèi)的采樣粒子數(shù),Ttr為平動溫度,Trot為轉(zhuǎn)動溫度,ζ為轉(zhuǎn)動自由度。
采用經(jīng)典微Poiseuille流驗(yàn)證微尺度DSMC方法的有效性。二維微通道幾何模型如圖 1所示, 計算條件取自文獻(xiàn)[12], 見表 1。其中,H為 0.4 μm,L為2 μm。左右邊界為自由來流, 模擬氮?dú)庠谖⑼ǖ纼?nèi)流動特性, 上下邊界為采用漫反射模型的固體壁面。無時間計時器技術(shù)(no time counter, NTC)和可變硬球(variable hard sphere, VHS)模型被用來模擬粒子之間的碰撞。
圖1 微通道基本結(jié)構(gòu)Fig. 1 Basic structure of the microchannel
表1 微通道計算參數(shù)Table 1 Numerical parameters for microchannel
圖2給出了兩種壓力比的中心線壓力分布與壁面滑移速度us變化的比較結(jié)果。圖2(a)首先將本文使用特征壓力邊界條件的DSMC計算結(jié)果與Arkilic等[13]通過1階解析壓力分布作比較, 結(jié)果表明本文計算結(jié)果與解析解具有較高的一致性, 驗(yàn)證了本文程序的準(zhǔn)確性。同時將壓力與壁面滑移速度與Fang[12]使用隱式邊界條件的結(jié)果進(jìn)行對比, 兩種分布量的模擬結(jié)果均符合較好, 進(jìn)一步證明本文方法的有效性。
(a) Centerline pressure distribution
(b) Wall slip velocity distribution圖2 微通道模擬值與文獻(xiàn)值的對比Fig. 2 Comparison between simulation and reference results
DSMC方法的本質(zhì)是在一段時間步長Δt內(nèi), 使得分子運(yùn)動和碰撞能夠解耦。但為了確保能達(dá)到正確解耦的目的, 一般需要保證以下兩個條件[14]: (1) 網(wǎng)格尺寸小于分子平均自由程λ, 一般取為平均自由程的1/3; (2) 時間步長Δt小于平均碰撞時間τ。下面研究其在微通道內(nèi)的適用性和可擴(kuò)展性。為此, 首先定義一個無量綱參數(shù)CSL表示網(wǎng)格尺寸與平均自由程之間的關(guān)系。同時借用計算流體力學(xué)中CFL數(shù)來定義一個新的參數(shù)DCFL數(shù)限制時間步長, 表示如下
以壓力比為2.5的微通道來驗(yàn)證網(wǎng)格尺寸和時間步長對計算結(jié)果的影響。計算了6種不同CSL數(shù)和DCFL數(shù)組合工況, 如表 2所示。其中, 工況3為第2節(jié)驗(yàn)證算例中壓力比為2.5的計算參數(shù), 以該參數(shù)下的計算結(jié)果作為比較的基準(zhǔn)。由表中數(shù)據(jù)可知, 除了工況6的DCFL數(shù)等于1外, 其余情況DCFL數(shù)均小于1。
表2 網(wǎng)格尺寸與時間步長約束參數(shù)Table 2 Constraint parameters for grid size and time step
圖3給出了各種工況下中心線壓力分布與解析壓力值對比情況, 可以看出除了工況6之外, 其余工況的壓力分布均與解析值符合良好。為了進(jìn)一步分析時間步長與網(wǎng)格參數(shù)對其他物理量的影響, 圖 4, 5分別給出了中心截面處速度分布與沿通道中心線溫度分布, 可以看出DCFL數(shù)小于1的工況1, 2, 4, 5的速度與溫度分布結(jié)果與經(jīng)過驗(yàn)證的工況3的結(jié)果趨于一致, 在工況6情況下, 中心截面速度與溫度的計算結(jié)果誤差均出現(xiàn)較大的偏離。綜合比較圖3, 4, 5的計算結(jié)果, 可以得到在選取網(wǎng)格尺寸和時間步長時, 首先應(yīng)當(dāng)保證DCFL數(shù)的取值小于1作為時間步長的約束條件。
圖3 中心線壓力分布對比Fig. 3 Comparison of centerline pressure distribution
在滿足DCFL數(shù)小于1的情況下, 進(jìn)一步分析CSL數(shù)對中心截面速度與溫度計算結(jié)果的影響。從圖 4看出工況1中心截面的流向速度明顯偏小, 而工況2, 4和5均能得到與工況工況3較為一致的模擬結(jié)果, 從而表明在CSL數(shù)小于1.865時也可以得到較為準(zhǔn)確的計算結(jié)果。從圖 5(b)的局部放大圖中可以看出工況5的中心線溫度較工況 3 明顯偏低。由此可見, 在選擇網(wǎng)格尺寸時不需要過分追求小CSL數(shù), 當(dāng)網(wǎng)格尺寸很小時, 是同一個網(wǎng)格內(nèi)模擬分子碰撞的次數(shù)減少, 造成一定統(tǒng)計誤差, 并且由于網(wǎng)格數(shù)量變多, 流場收斂需要的時間變長, 增加了運(yùn)算成本。因此, 在確保DCFL數(shù)小于1的前提下, 選擇CSL數(shù)接近于1可作為本文網(wǎng)格尺寸選取的原則。
圖4 中心截面處流向速度分布對比Fig. 4 Comparison of velocity distribution at the channel midspan
(a) Centerline temperature distribution
(b) Enlarged red frame圖5 中心線溫度分布對比Fig. 5 Comparison of centerline temperature distribution
研究模型為具有單/雙孔口結(jié)構(gòu)的微通道, 如圖 6所示, 兩種不同孔口結(jié)構(gòu)的尺寸參數(shù)見表 3, 壁面采用漫反射模型, 模擬氣體為氮?dú)猓?不同入口壓力條件下[9]的來流條件見表 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
圖 7分別給出了不同入口壓力下, 單/雙孔口模型通道內(nèi)的流線與壓力等值線分布圖。在單孔口模型中, 通道上游處流線沿拐角進(jìn)入孔, 在孔口內(nèi)部流線與通道平行, 主流經(jīng)孔口加速, 在孔口出口形成一股收縮流, 由于管道突然擴(kuò)張, 流動減速以充滿整個通道??卓诔隹谟捎诹鲃铀俣冉档停?壓力對應(yīng)升高, 在孔口下游的拐角處形成反向壓力梯度, 從而出現(xiàn)回流區(qū)。在雙孔口模型中, 具有與單孔口模型相似的流動結(jié)構(gòu), 且每一個孔口后均存在回流區(qū)。此外, 當(dāng)入口壓力增大之后, 由于出口壓力始終保持不變, 逆壓區(qū)的壓差增大, 回流區(qū)尺寸隨之變大。在雙孔口模型中, 重點(diǎn)關(guān)注兩孔口之間的過渡區(qū)域, 入口壓力為1.5×105Pa時, 在過渡區(qū)拐角處出現(xiàn)了兩處分離, 孔口2的上游也出現(xiàn)一個很小的分離區(qū)。這是由于流體在過渡區(qū)流程較短, 流體經(jīng)過孔口壓縮加速, 無法及時調(diào)整流向, 出現(xiàn)分離區(qū)。而隨著入口壓力的增加, 分離區(qū)的尺寸增大, 兩處分離逐漸融合, 形成一個大的分離渦。
(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)流道收縮使得流體加速時, 流體的靜壓能轉(zhuǎn)換為動能, 孔口上游的壓力快速下降。此外, 最小壓力不是在截面最小的孔口處, 而是出現(xiàn)在孔口下游, 隨著入口壓力的增大, 最小壓力的位置也隨之后移, 并且孔口后的壓力甚至?xí)∮诔隹趬毫?。在最小壓力位置的下游?隨著動能轉(zhuǎn)換回靜壓能, 流動開始減速, 壓力逐漸增加至局部最大值, 即出口壓力。
(a) Single orifice model
(b) Double orifice model圖8 中心線壓力分布Fig. 8 Centerline pressure distribution
圖 9, 10分別為不同入口壓力下通道中心線的速度分布和溫度分布, 可以看出速度與溫度的變化規(guī)律相反??卓诮Y(jié)構(gòu)類似于一個收縮/擴(kuò)張的裝置, 因此氣體在經(jīng)過孔口時被壓縮, 導(dǎo)致孔口處的速度增加, 其后, 由于通道擴(kuò)張, 速度減小, 流場中溫度變化與速度變化相反。在單孔結(jié)構(gòu)中, 可以發(fā)現(xiàn)氣流經(jīng)過孔口時速度上升非???, 由于慣性作用, 在離開孔口后速度達(dá)到最大, 再緩慢降低, 且流場的最大速度隨著入口壓力的增大而增大。在雙孔模型中, 發(fā)現(xiàn)孔口1的速度變化幅度較小, 一方面孔口1的位置處于流場的前半部分, 本身速度就較低; 另一方面, 氣體速度與壓降大小有關(guān)。當(dāng)入口壓力為1.5×105Pa時, 由于兩次壓降大小幾乎相同, 因此經(jīng)過兩個孔口速度增加量相似; 而當(dāng)入口壓力逐漸增加, 經(jīng)孔口2的壓降比孔口1要大得多, 因此, 孔口2的速度變化更劇烈。對比同一入口壓力下兩個模型的速度大小, 單孔口的最大速度要大于雙孔口模型, 由于雙孔口模型經(jīng)兩次壓降, 在相同壓差情況下, 雙孔口的單次壓降小, 使得流體速度每次增加的幅度也小。另外, 由于雙孔口流場產(chǎn)生兩處分離區(qū), 流場擾動作用更加明顯, 對流體的速度也產(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
本文基于壓力邊界條件的DSMC方法, 研究了不同網(wǎng)格尺寸和時間步長對直微通道模型結(jié)果精度的影響, 得到了在微尺度下的網(wǎng)格尺寸和時間步長選取的一般準(zhǔn)則。在此基礎(chǔ)上, 研究了不同入口壓力條件對單/雙孔口結(jié)構(gòu)微通道內(nèi)流動過程的影響。主要得出以下結(jié)論:
(1) 通過CSL數(shù)和DCFL數(shù)兩個無量綱參數(shù)給出了網(wǎng)格尺寸和時間步長的約束條件, 研究結(jié)果表面在滿足DCFL數(shù)嚴(yán)格小于1的情況下, 選擇CSL數(shù)接近或小于1, 可以作為網(wǎng)格尺寸與時間步長的選取原則。
(2)由于壓降的存在使得孔口后形成低壓區(qū), 因此氣體經(jīng)過孔口后會在拐角處形成分離區(qū), 并且分離區(qū)大小隨入口壓力的增大而增大。雙孔口模型中, 兩孔口之間的過渡區(qū)域流程較短, 但也出現(xiàn)了分離現(xiàn)象, 并隨入口壓力增大, 過渡區(qū)兩處分離渦逐漸融合。在相同壓差情況下, 通過雙孔口模型流速會低于單個孔口通道的流速。
(3)在孔口模型中, 氣流在通過孔口時產(chǎn)生較大壓降。經(jīng)孔口處的氣流速度增加, 并隨入口壓力增大而增大。最大速度的位置由于慣性作用逐漸后移, 出現(xiàn)在孔口下游, 其后, 速度逐漸降低。溫度與速度變化相反。