王曉偉，武 越，杜春林，韓 峰，張志軍*，張世偉

（1. 東北大學(xué) 機(jī)械工程與自動化學(xué)院，沈陽 110819； 2. 北京衛(wèi)星環(huán)境工程研究所，北京 100094）

0 引言

在沒有任何初始壓力梯度的稀薄氣體中，非均勻溫度場會誘導(dǎo)氣體流動，且浸沒在氣體中的結(jié)構(gòu)會受到輻射力[1-3]。該現(xiàn)象最早發(fā)現(xiàn)于19 世紀(jì)后期William Crookes 爵士發(fā)明的輻射計中[4]。該輻射計由安裝在抽真空的玻璃燈泡內(nèi)部的一組薄葉片組成，葉片一面涂成黑色，另一面為光亮的白色。將該裝置放置于陽光下，葉片開始旋轉(zhuǎn)，旋轉(zhuǎn)方向為從黑色一面指向白色一面。

隨著航空航天科學(xué)技術(shù)和微機(jī)電系統(tǒng)（MEMS）的發(fā)展，研究人員對應(yīng)用輻射力的興趣日益濃厚。Jones 異想天開地提出基于輻射力的直升機(jī)[5]和微型巡航導(dǎo)彈[6]?；谳椛淞Φ耐七M(jìn)系統(tǒng)還可用在臨近空間飛行器上[7]，由安裝在航天器表面的多個方形葉片陣列為航天器提供升力或者阻力補(bǔ)償力。研究發(fā)現(xiàn)，在40～80 km 的高度上，這種推進(jìn)系統(tǒng)葉片上的輻射力能夠達(dá)到25 mN/m2。此外，由于對環(huán)境壓力、溫度以及氣體成分較為敏感，輻射力也被用于MEMS 傳感器中的壓力測量[8]和氣體檢測[9-20]。

最近，Strongrich 等[9-10]提出一種基于輻射力驅(qū)動的設(shè)備（Microelectromechanical In-plane Knudsen Radiometric Actuator, MIKRA），可用于氣體檢測。研究表明，物理和幾何參數(shù)顯著影響該設(shè)備的性能[21-22]：當(dāng)冷熱微梁溫差從10 K 增大到100 K 時，最大輻射力將增大7 倍以上；輻射力隨環(huán)境壓力的增大呈現(xiàn)出“鐘形”分布規(guī)律；微梁寬度的增大將使設(shè)備整體所受輻射力減小，且最大輻射力在更小的環(huán)境壓力下取得；微梁高度增大3 倍時，最大輻射力將增大3.66 倍；當(dāng)梁的高度與寬度的比值增加9 倍時，最大輻射力能增大5.33 倍。需要指出的是，在這些研究中氣體與所有梁表面之間的相互作用規(guī)律都是基于完全漫反射假設(shè)；然而，在實際應(yīng)用中這種假設(shè)并非總是成立的。有關(guān)研究表明，氣-面相互作用規(guī)律（表面調(diào)節(jié)系數(shù)）對流場結(jié)構(gòu)、傳熱特性、壓力分布以及輻射力都有顯著的影響[23-25]。Zhou 等的研究表明[23]，法向能量調(diào)節(jié)系數(shù)（NEAC）和切向動量調(diào)節(jié)系數(shù)（TMAC）對微通道內(nèi)的壓力分布和流場結(jié)構(gòu)會產(chǎn)生很大的影響；并且，相比于NEAC，改變TMAC 對流動特性的影響更為顯著。此外，Kalempa 和Wang 等[24-25]研究了氣-面相互作用規(guī)律對稀薄氣體中球體和微梁周圍流體特性的影響，發(fā)現(xiàn)對于輻射力，改變NEAC 和TMAC 會產(chǎn)生相反的作用。值得注意的是，在文獻(xiàn)[25]的研究中，微梁（一對冷熱梁）上的NEAC 和TMAC 是同步變化的，并沒有考慮其中一根梁上表面調(diào)節(jié)系數(shù)單獨變化的影響情況。

為了補(bǔ)充氣-面相互作用規(guī)律對輻射力影響的理論研究，本文采用直接模擬蒙特卡羅（DSMC）方法研究微梁周圍的稀薄氣體流動特性，主要考查過渡流態(tài)下單梁表面調(diào)節(jié)系數(shù)變化對流場結(jié)構(gòu)、傳熱特性、壓力分布以及輻射力等的影響，以期為相關(guān)器件的設(shè)計制造和應(yīng)用提供理論參考與指導(dǎo)。

1 氣-面相互作用模型

常見的氣-面相互作用模型包括鏡面反射、漫反射、Maxwell 以及CL（Cercignani-Lampis）模型。對于不同的氣-面相互作用模型，氣體分子撞擊壁面并發(fā)生反射的速度分布規(guī)律[26]如圖1 所示。鏡面反射模型是最簡單的氣-面相互作用模型，它假設(shè)撞擊到表面的氣體分子類似于完全彈性的球體一樣被反射，即，氣體分子的法向速度分量反向而切向速度分量保持不變，如圖1(a)所示，這意味著在這種氣-面相互作用過程中，氣體分子與壁面之間不發(fā)生熱傳遞和摩擦行為。對于漫反射模型，撞擊到壁面的氣體分子首先與表面進(jìn)行能量交換，使氣體分子溫度和壁面溫度達(dá)到熱平衡；然后根據(jù)由壁面溫度確定的半范圍Maxwell 分布，隨機(jī)分配反射分子的速度，如圖1(b)所示。前面這2 種氣-面相互作用模型都不太切合實際，因此出現(xiàn)了Maxwell模型。實際上，Maxwell 模型是漫反射與鏡面反射模型的組合[27]，在模型中引入調(diào)節(jié)系數(shù)α（0≤α≤1）來指定入射分子在壁面上的反射形式，即，入射分子的α部分為漫反射，剩余的(1-α)部分為鏡面反射。相比于Maxwell 模型，CL 模型能更真實地反映物理現(xiàn)實中反射分子的速度和能量分布[28]，如圖1(c)所示。該模型包括2 個表面調(diào)節(jié)系數(shù)（NEAC和TMAC），二者相互獨立，對于一般的光滑表面，這2 個調(diào)節(jié)系數(shù)的數(shù)值均在0～1 之間變化。本研究的所有DSMC 模擬中，氣-面相互作用模型均選用CL 模型。

圖1 不同氣-面相互作用模型中反射分子的速度分布示意Fig. 1 Schematics diagram of velocity distributions of reflected molecules based on different gas-surface interaction models

2 問題陳述與數(shù)值方法

由于MIKRA 設(shè)備中冷熱梁的長度（1.4 mm）遠(yuǎn)大于其橫截面尺寸（50 μm×50 μm），所以只需進(jìn)行二維模擬以近似表征實際的氣體流動情況，模擬域如圖2 所示。左右兩表面為對稱邊界，上方是自由流邊界以表征H2-N2混合氣體的流入或流出，下方為基底壁面邊界。冷熱梁懸浮于距基底上方4 μm處。冷梁和熱梁具有相同的尺寸，且在模擬過程中假設(shè)兩者的間距g保持不變，為20 μm。此間距也用于計算Knudsen 數(shù)（Kn）。

圖2 模擬域和邊界條件Fig. 2 Simulation domain and boundary conditions

DSMC 方法基于粒子求解Boltzmann 方程[29]，使用大量的模擬粒子再現(xiàn)實際氣體分子的運(yùn)動情況，跟蹤并記錄這些粒子的位置、速度和內(nèi)能等信息，最后對模擬粒子的微觀值進(jìn)行平均從而獲得宏觀量。該方法為了降低統(tǒng)計誤差，需要較長的模擬時間。本文使用開源的DSMC 求解器dsmcFoam-Plus[30]來模擬微梁周圍的氣體流動；采用無時間計數(shù)器（NTC）方案選取粒子碰撞對；另外，因為模擬域氣體介質(zhì)是50%H2+50%N2的雙原子混合氣體[13,18,25]，所以采用變徑軟球（VSS）模型和Larsen-Borgnakke 模型分別模擬粒子間的碰撞和粒子本身動能與內(nèi)能的交換過程。

所有模擬均在環(huán)境壓力P=387 Pa（Kn=0.74）的氣氛中進(jìn)行；冷梁和均勻加熱的熱梁的溫度分別為Tc=304 K 和Th=341.5 K。值得一提的是，這2 個溫度是通過在純N2氣氛中將相應(yīng)壓力條件下測得的實驗數(shù)據(jù)[9-10]取平均值得到的。根據(jù)這些參數(shù)，使用Δx=Δy=2 μm 的網(wǎng)格均勻地劃分模擬域[10,13,18,25]，每個網(wǎng)格中的模擬粒子數(shù)不少于20 個。模擬計算的時間步長Δt=1 ns，總模擬時間不短于5 ms[10,13,18,25]；另外，對于表面調(diào)節(jié)系數(shù) NEAC（αn）和TMAC（σt）均獨立設(shè)為0.2、0.4、0.6、0.8 和1。本研究中只單獨改變冷梁的表面調(diào)節(jié)系數(shù)，而熱梁的表面調(diào)節(jié)系數(shù)恒為1（對于冷、熱梁同步變化的情況，請參見文獻(xiàn)[25]）。在物理實際中，可通過改變制造冷梁和熱梁所用的材料來實現(xiàn)其表面調(diào)節(jié)系數(shù)的調(diào)整。本文為區(qū)分冷梁和熱梁的表面調(diào)節(jié)系數(shù)，分別以上標(biāo)Cold 和Hot 予以標(biāo)記。

3 結(jié)果與討論

3.1 流場結(jié)構(gòu)

圖3 不同NEAC 下的速度流線和溫度分布云圖Fig. 3 Velocity streamlines and temperature contours for different NEACs (In each of the sub-figures, the left half shows the results of αnCold=αn, αnHot=1, the right half shows the results extracted from Ref. [25], αnCold=αnHot=αn)

TMAC 變化對微梁附近流場的影響如圖4 所示，同樣在每幅圖中均給出了本文和文獻(xiàn)[25]的模擬結(jié)果。不難看出，相比于NEAC 對流場結(jié)構(gòu)的影響，切向TMAC 的影響力較弱，且在僅改變冷梁TMAC 以及冷梁和熱梁TMAC 同步變化的情況中均是如此。這一結(jié)果對微梁制作材料的選取具有一定的指導(dǎo)意義。同樣由于熱梁TMAC 的降低導(dǎo)致氣體與壁面交換的能量減少，使得冷梁和熱梁TMAC 同步改變所引起的流場結(jié)構(gòu)的變化比僅改變冷梁TMAC 的情況更為明顯，且流場溫度更高。

圖4 不同TMAC 下的速度流線和溫度分布云圖Fig. 4 Velocity streamlines and temperature contours for different TMACs (In each of the sub-figures, the left half shows the results of =σt, =1, the right half shows the results extracted from Ref. [25], ==σt)

3.2 傳熱特性

不同NEAC 對冷梁各表面?zhèn)鳠崃康挠绊懭鐖D5所示，每幅圖中虛線數(shù)據(jù)為文獻(xiàn)[25]的模擬結(jié)果，下同。傳熱量數(shù)值為正表示該表面處于吸熱狀態(tài)，為負(fù)則表示其處于放熱狀態(tài)?？梢钥吹剑艿礁邷?zé)崃旱挠绊懀淞旱纳媳砻婧陀冶砻嫣幱谖鼰釥顟B(tài)；對于下表面，由于其溫度比底部基底的更高且距離較近，故處于放熱狀態(tài)；左表面處于放熱狀態(tài)的原因是，其遠(yuǎn)離熱梁且表面溫度比其周圍氣體的溫度高。正如預(yù)期的那樣，NEAC 的減小導(dǎo)致表面?zhèn)鳠崃康慕^對值減小，即，表面的吸熱或放熱能力減弱。另外，對于僅改變冷梁的NEAC 的情況，表面的吸熱或放熱能力比冷梁和熱梁的NEAC 同步變化時更強(qiáng)。由于冷梁右表面距離熱梁較近，這種現(xiàn)象尤為明顯，如圖5(d)所示。

圖5 不同NEAC 對冷梁各表面?zhèn)鳠崃康挠绊慒ig. 5 Heat transfer distributions along surfaces of cold beam with different NEACs (In each of the sub-figures, dashed lines show the results extracted from Ref. [25])

不同TMAC 對冷梁各表面?zhèn)鳠崃康挠绊懭鐖D6所示?？梢钥闯?，各表面的吸熱或放熱狀態(tài)以及傳熱分布規(guī)律與改變NEAC 時的一致。從傳熱量上來看，改變TMAC 比改變NEAC 能夠獲得更大的絕對值。

圖6 不同TMAC 對冷梁各表面?zhèn)鳠崃康挠绊慒ig. 6 Heat transfer distributions along surfaces of cold beam with different TMACs (In each of the sub-figures, dashed lines show the results extracted from Ref. [25])

這對于物理應(yīng)用中的熱管理具有一定的參考價值，即，針對材料不同的微梁需要配備不同換熱能力的散熱系統(tǒng)。另外，僅改變冷梁的TMAC 獲得的結(jié)果與冷梁和熱梁的TMAC 同步變化所得到的結(jié)果之間的差異在TMAC 較大的情況下并不明顯，只在TMAC 較小時才能凸顯。

3.3 壓力分布特性

圖7(a)和圖7(b)分別為NEAC 和TMAC 變化對冷梁左右表面壓力分布的影響。

圖7 不同表面調(diào)節(jié)系數(shù)下冷梁左右表面的壓力分布Fig. 7 Pressure distributions along left and right surfaces of the cold beam with different surface accommodation coefficients (In each of the sub-figures, dashed lines show the results extracted from Ref. [25])

由圖7(a)可以看到，僅改變冷梁的NEAC 在引起梁表面壓力變化方面與冷梁和熱梁的NEAC 同步變化的情況顯示出完全不同的結(jié)果。更確切地說，如果僅僅改變冷梁的NEAC，冷梁右表面壓力始終大于左表面壓力，且隨著NEAC 的減小，右表面壓力明顯增大。這意味著冷梁與熱梁兩者之間相互排斥。而對于冷梁和熱梁NEAC 同步變化的情況，則會出現(xiàn)冷梁左表面壓力大于右表面壓力的現(xiàn)象，且隨著NEAC 的降低右表面壓力明顯減小，這種現(xiàn)象更為顯著[25]。這意味著冷梁與熱梁兩者之間在NEAC 較大時是相互排斥的，而當(dāng)NEAC 較小時則表現(xiàn)為相互吸引。因此，在物理實際中可通過選用不同的材料來實現(xiàn)對微梁移動方向的控制。

由圖7(b)可以看到，無論是僅改變冷梁的TMAC還是同步改變冷梁和熱梁的TMAC，冷梁右表面壓力均始終大于左表面壓力，即，冷梁與熱梁兩者之間相互排斥。此外，隨著TMAC 的減小，右表面壓力明顯增大，且同步降低冷梁和熱梁的TMAC 在增大右表面壓力方面比僅降低冷梁的TMAC 的情況表現(xiàn)出更顯著的效果。

3.4 輻射力

在DSMC 中，施加在冷梁上的輻射力F可通過撞擊在冷梁表面的入射粒子和反射粒子之間的動量變化來計算，

式中：Nequ為一個DSMC 模擬粒子所代表的真實氣體分子數(shù)量；N為與表面相互作用的總粒子數(shù)；m為粒子的質(zhì)量；c為粒子的速度，上標(biāo)i 和r 分別表示入射粒子和反射粒子。

圖8(a)和圖8(b)分別為NEAC 和TMAC 變化對冷梁所受輻射力的影響，其中虛線為文獻(xiàn)[25]的結(jié)果。

由圖8(a)可以看出，對于僅改變冷梁的NEAC與冷梁和熱梁的NEAC 同步變化這2 種情況，冷梁所受輻射力對NEAC 的依賴性顯示出完全不同的情形。更確切地說，如果同步減小冷梁和熱梁的NEAC，輻射力隨著NEAC 的減小呈非線性減小，當(dāng)αn＜0.2 之后出現(xiàn)輻射力方向反轉(zhuǎn)的現(xiàn)象，且其絕對值具有隨著NEAC 的減小而增大的趨勢（輻射力曲線的斜率隨著NEAC 的減小而增大）；而對于僅減小冷梁的NEAC 的情況，輻射力隨著NEAC的減小呈近似線性增大。

由圖8(b)可以看出，在僅改變冷梁表面TMAC的情況中，輻射力隨著TMAC 的減小呈非線性減小；而對于冷梁和熱梁的TMAC 同步變化的情況，輻射力卻隨著TMAC 的減小呈非線性增大。同時，我們注意到，在這2 種情況中，隨著TMAC 的減小輻射力曲線的斜率（絕對值）均增大。

圖8 不同表面調(diào)節(jié)系數(shù)下冷梁所受輻射力Fig. 8 Radiometric force of the cold beam with different surface accommodation coefficients (In each of the sub-figures, dashed lines show the results extracted from Ref. [25])

4 結(jié)束語

本文采用DSMC 方法并基于CL 氣-面相互作用模型研究法向能量調(diào)節(jié)系數(shù)（NEAC）和切向動量調(diào)節(jié)系數(shù)（TMAC）對一組微梁周圍稀薄氣體流動的影響。結(jié)果表明，表面調(diào)節(jié)系數(shù)對流場結(jié)構(gòu)、傳熱特性、壓力分布和輻射力的影響是顯著的。為了在實際應(yīng)用中盡可能地獲得最大的輻射力，在設(shè)計和制造MEMS 器件時，如果冷梁和熱梁使用相同的材料，則應(yīng)該選取NEAC 較大而TMAC 較小的材料；如果冷梁和熱梁使用不同的材料，則冷梁應(yīng)該選具有較小NEAC 和較大TMAC 的材料，而熱梁應(yīng)該選NEAC 和TMAC 均較大（最好接近于1）的材料。