固體微推力器氣粒兩相羽流場的數值模擬①

張 斌，毛根旺，胡松啟，陳茂林

(西北工業大學燃燒、熱結構與內流場重點實驗室，西安 710072)

0 引言

固體微推力器是基于MEMS技術制造的微型固體火箭發動機，屬于化學推進，主要用于微型衛星的姿態控制和軌道機動［1－3］。微型衛星的主要任務是組成分布式星載載波雷達、實現三維立體成像和高分辨率合成孔徑對地遙感等，其上搭載許多精密光學儀器，器件集成度高［4－5］。因此，對固體微推力器的羽流研究不容忽視。

固體微推力器工作條件下，外界環境接近真空。微噴管噴出的燃氣流動已明顯帶有稀薄流動的特性。同時，燃氣中還含有少量金屬氧化物。因此，該流動本身是典型的氣粒兩相稀薄流動。對于氣粒兩相稀薄流動的實驗研究較為困難，微尺度更增加了其實驗研究的難度。目前，在該領域已開展的研究工作多集中在對氣粒兩相稀薄流的數值研究。2001年，美國Sandia國家實驗室的Gallis等［6］提出了稀薄流中氣相分子對顆粒力和熱的作用模型，開辟了氣粒兩相稀薄流動數值研究的新領域。2004年，美國密歇根大學的Burt等［7］利用唯象論原理求解顆粒相對氣相分子的作用，同時結合Gallis的模型，較好地模擬了兩相間動量和能量的傳遞。2005年，美國南加州大學的Gimelshein等［8］充分發展了Burt的數學模型，結合唯象論原理，成功進行了航天器姿控發動機的氣粒兩相羽流場模擬。國內對該領域的實驗和數值研究報道較少，目前有見報道的只有國防科大的李潔等［9］，采用能量均分原理處理顆粒相對氣相作用導致氣相平動能和轉動能的變化，即采用熱力學模型處理顆粒相對氣相的作用。

數值求解氣粒兩相稀薄流包含3部分內容:(1)計算氣體分子自身的運動和碰撞;(2)計算顆粒自身的運動;(3)計算顆粒相與氣相間的動量和能量交換。

文中采用DSMC/EPSM混合算法計算多組分稀薄氣體自身流動，采用Gimelshein的模型求解氣相對顆粒相力和熱的作用，采用唯象論原理求解顆粒相對氣相的作用。通過以上數值方法，模擬了固體微推力器羽流場中氣粒兩相參數的分布規律，得到了令人滿意的結果。

1 數值計算方法

1.1 多組分氣體的DSMC/EPSM算法

Bird［10］提出的DSMC方法是一個用于求解玻爾茲曼方程的經典方法，也是目前求解稀薄氣體流動的主要手段。在DSMC方法中，大量的真實氣體分子被少數仿真分子代替，仿真分子按常速運動直到與另一個仿真分子碰撞。碰撞后，仿真分子的速度和內能隨之改變，但并不改變位置。其實質是將仿真分子的運動和碰撞有效解耦。計算中，大量計算機時被用于模擬仿真分子的碰撞過程。在多組分氣體的DSMC算法中，需分別跟蹤各組分氣體仿真分子的運動和碰撞。在碰撞計算中，除考慮各組分氣體仿真分子自身的碰撞外，還需考慮不同種類氣體仿真分子間的碰撞，其控制方程為混合氣體的Boltzman方程。

EPSM算法是目前求解過渡區流動最有效的數值手段。在1個時間步長下，如果網格內仿真分子的碰撞數足夠大，則經過1個時間步長后，網格中仿真分子的速度將達到局部平衡態，仿真分子對碰撞前后滿足動量和能量守恒，其速度抽樣分布即為Maxwell-Boltzman平衡態分布。EPSM算法不需處理網格中仿真分子的碰撞過程，所有仿真分子在一個時間步長后的狀態僅決定于其初始狀態，從而極大提高了程序的計算效率。對于多組分氣體的EPSM算法，需先求解網格內仿真分子的質量平均速度和熱能，然后對每一組分每一自由度重新建立能量分布。算法中有2項假設:(1)每一組分的平均速度等于網格中仿真分子的質量平均速度;(2)儲存在每一組分每一自由度中的熱能相等。

文中根據Bird的理論，認為1個時間步長內仿真分子碰撞數大于10的網格單元，即達到局部平衡，采用EPSM算法處理仿真分子的碰撞過程，從而構造出多組分氣體的DSMC/EPSM算法。

1.2 顆粒運動的描述

對顆粒運動的描述與對氣體分子運動的描述類似。計算中，只需將顆粒視為一種仿真分子，通過顆粒的質量、直徑、位置、速度和內能等參數表征顆粒屬性，即可用DSMC方法模擬顆粒的隨機運動。然而，與模擬氣體分子的不同之處在于顆粒速度和內能的變化，不是通過顆粒間碰撞得到的，而是通過求解周圍氣體分子對顆粒力和熱的作用，得到下1個時間步后顆粒的速度和溫度。

1.3 分子對顆粒的作用

Gimelshein的數學模型是目前有見報道的，并能最為有效描述稀薄流中氣相對顆粒相力和熱作用的數學模型。模型中做了以下假設:(1)顆粒是球形的;(2)顆粒內部沒有溫度梯度;(3)顆粒在氣相中是稀疏的;(4)忽略氣相和顆粒相間的傳質。理論上，分子對顆粒力和熱的作用，是由與顆粒發生碰撞的分子(包括入射分子、鏡面反射分子、等溫反射分子和絕熱反射分子)的動量和能量的變化率決定的。另一方面，單個顆粒所受到周圍氣體分子的力和熱，可通過相應網格中所有分子對該顆粒作用的統計得到。基于以上理論，Gimelshein給出了網格內單個仿真氣體分子對單個顆粒力和熱的作用公式分別為［8］

式中 m為氣體相對分子質量;Fnum為氣體分子權因子;為氣體分子相對于顆粒的相對速度矢量;Vcell為網格體積;g=||;Tp為顆粒溫度;ζrot為氣體分子的轉動能自由度;αm為動量調節系數;αt為平動能調節系數;αr為轉動能調節系數;αv為振動能調節系數;Erot為單個氣體分子的轉動能;Evib為單個氣體分子的振動能。

在每一時間步后，根據式(1)和式(2)計算得到顆

粒的速度和溫度表達式如下［8］:

1.4 顆粒對分子的作用

顆粒對分子的作用體現在碰撞后分子速度矢量的變化，其物理本質基于分子與顆粒表面碰撞的唯象論原理。

網格內分子與顆粒的碰撞數 ns由式(5)確定［7，11］。

式中 Np為顆粒權因子;ng為仿真分子數;(cr)max為分子與顆粒的最大相對速度。

式中 ur、vr和wr為分子與顆粒碰撞前的相對速度分量;δ為碰撞前后分子入射方向與反射方向的夾角，滿足分布函數:

其中，ε為反射速度矢量的方位角，服從［0，2π］的均勻分布;cr和分別為碰撞前后分子與顆粒間的相對速度，其中按“取舍法”由分布f()=2β43exp抽樣得到，β=［m/(2kBTp)］1/2。

2 計算條件

研究對象為固體微推力器陣列中的單元固體微推力器羽流場，采用軸對稱模型計算。噴管出口半徑3×10－4m。氣相為CO、H2和CO2組成的混合氣體，其摩爾比為0.5∶0.3∶0.2。通過求解N-S方程，得到微噴管出口截面上兩相的相關參數，采用相鄰網格節點值平均后作為相應入口處的邊界條件。顆粒相為Al2O3，直徑 1 ×10－6m，流量為2.5 ×10－1kg/(s·m2)，密度為3 970 kg/m3，定容比熱為765 J/(kg·K)。推力器壁面溫度300 K，采用漫反射邊界條件。分子碰撞模型采用VHS模型。

3 結果與分析

3.1 驗證算例

計算區域長20 mm，寬0.1 mm，壁面為鏡面反射邊界，如圖1所示。來流氣體為H2、CO和N2的混合氣體，來流分子數密度分別為 2×1023、1×1023、1×1023/m3，來流氣體速度 2 000 m/s，溫度 1 000 K。顆粒相為直徑3×10－6m和6×10－6m各一半的球形Al2O3，速度為 1 200 m/s，溫度為 2 200 K，流量為13.33 kg/(s·m2)。算例的計算結果見圖2。

圖1 算例的幾何尺寸和邊界類型示意圖Fig.1 Geometry dimensions and boundary types of numerical example

圖2 算例的計算結果Fig.2 Results of numerical example

算例結果與文獻［7，11］報道的結果(見圖3)基本一致。流動過程中，兩相間的動量傳遞導致軸線上沿流動方向氣相速度降低、顆粒速度增加;兩相間的能量傳遞導致軸線上沿流動方向的氣相溫度升高、顆粒溫度降低。

圖3 Burt 的計算結果Fig.3 Results of numerical example reported by Burt

3.2 固體微推力器羽流場中氣粒兩相參數分布規律

圖4(a)～(c)分別為氣相的壓強場、溫度場和軸向速度場，圖5(a)和(b)分別為顆粒相的溫度場和軸向速度場。各圖中，壓強、溫度和速度的單位分別為Pa、K 和 m/s。

由圖4可知，由于固體微推力器推進劑燃氣中顆粒含量少，顆粒相對氣相壓強場、溫度場和軸向速度場的影響并不顯著，各參數的分布規律基本與純氣相羽流場一致，羽流軸線附近區域氣相軸向速度在徑向方向變化幅度不大，沒有明顯的兩相界面效應存在，且未出現Gimelshein等［8］報道的V型氣相軸向速度場。

由圖5可知，氣相對顆粒相的作用顯著，沿羽流區軸線方向，顆粒相的軸向速度增加，同時存在于兩相間的能量傳遞導致顆粒相溫度降低。

3.3 氣相組分分離效應

圖6(a)為羽流區軸線上氣相各組分真實分子數百分數變化曲線，圖6(b)為羽流區壁面附近各組分真實分子數百分數變化曲線。

由圖6(a)可知，大分子量組分CO和CO2的分子數密度沿羽流區軸線逐漸增加，而小分子量組分H2的分子數密度沿羽流區軸線逐漸降低，說明H2在羽流場中的擴散效應明顯大于CO和CO2。

由圖6(b)可知，在微推力器壁面與噴管出口相鄰的極小區域內，氣相由少量CO和大量H2組成;隨著距噴管出口處距離增加，壁面附近的氣相中僅有H2存在;在噴管出口以上的壁面附近，幾乎沒有大分子量組分CO2存在。可見，推進劑燃氣中H2組分是造成羽流對固體微推力器壁面侵蝕和附著的主要氣體。

圖4 羽流場中氣相的壓強、溫度和軸向速度等值線Fig.4 Pressure contours，temperature contours and axial velocity of gas phase in plume field

圖5 羽流場中顆粒相的溫度和軸向速度等值線Fig.5 Temperature contours and axial velocity of gas phase in plume field

圖6 氣相(沿軸線和沿壁面)真實分子數百分數曲線Fig.6 Real molecular number percentage in gas phase(along axial line and along wall)

4 結論

通過Gimelshein的模型和唯象論原理描述稀薄流中氣相和顆粒相的雙向耦合作用，成功模擬了固體微推力器羽流場中的氣粒兩相稀薄流動。在羽流場中，氣相與顆粒相間的動量和能量傳遞能顯著影響顆粒相的速度和溫度，沿羽流流動軸線方向上顆粒相的速度增加、溫度降低;同時，由于顆粒相的相對含量較少，使其對氣相壓強場、溫度場和軸向速度場的影響并不顯著。羽流場中，氣相呈現出顯著的組分分離效應，小分子量氣體組分在羽流場中的高度擴散是造成羽流壁面侵蝕和附著效應的主要原因。

［1］Rossi C，Rouhani M D，Estève D.Prediction of the performance of a Si-micromachined microthruster by computing the subsonic gas flow inside the thruster［J］.Sensors and Actuators，2000，87:96-104.

［2］Rossi C，Orieux S，Larangot B，et al.Design，fabrication and modeling of solid propellant microrocket-application to micropropulsion［J］.Sensors and Actuators A，2002，99:125-133.

［3］Zhang K L，Chou S K，Ang S S，et al.A MEMS-based solid propellant microthruster with Au/Ti igniter［J］.Sensors and Actuators A，2005，122:113-123.

［4］張高飛，尤政，胡松啟，等.基于MEMS的固體推進器陣列［J］.清華大學學報(自然科學版)，2004，44(11):1489-1452.

［5］周海清，張高飛，尤政.固體微推力器設計與數值研究［J］.推進技術，2007，28(3):230-234

［6］Gallis M A，Torczynski J R，Rader D J.An approach for simulating the transport of spherical particles in a rarefied gas flow via the direct simulation Monte Carlo method［J］.Physics of Fluids，2001，13(11):3482-3492.

［7］Burt J M，Boyd I D.Development of a two-way coupled model for two phase rarefied flows［R］.AIAA 2004-1351.

［8］Gimelshein S F，Alexeenko A A，Wadsworth D C，et al.The influent of particulates on thruster plume/shock layer interaction at high altitudes［R］.AIAA 2005-766.

［9］李潔，任兵，陳偉芳.稀薄流過渡區氣固兩相噴流的建模與數值模擬［J］.空氣動力學學報，2005，23(4):484-489.

［10］Bird G A.Molecular gas dynamics and the direct simulation of gas flows［M］.Oxford:Clarendon Press，1994.

［11］Burt J M.Monte carlo simulation of solid rocket exhaust plumes at high altitude［D］.University of Michigan，2006.