微型固體推力器瞬態(tài)工作過(guò)程數(shù)值模擬①

李 騰，方蜀州，劉旭輝，馬紅鵬，湯 旭，趙 婷

(1.北京理工大學(xué)宇航學(xué)院，北京 100081;2.中國(guó)空間技術(shù)研究院北京控制工程研究所，北京 100190)

0 引言

基于MEMS技術(shù)制造的固體微推力器具有結(jié)構(gòu)簡(jiǎn)單、功耗低、無(wú)可動(dòng)部件及其帶來(lái)的摩擦損失小、可靠性高等諸多優(yōu)勢(shì)。固體微推力器最大的優(yōu)勢(shì)是集成度高，可集成在一塊基片上組成陣列。目前，陣列式推進(jìn)系統(tǒng)已在動(dòng)能攔截器等方面投入應(yīng)用，如美國(guó)的增程攔截彈ERINT-1?；谏鲜鎏攸c(diǎn)，固體微推力器成為微推進(jìn)系統(tǒng)方面新的研究熱點(diǎn)。

目前，對(duì)固體微推力器工作特性的研究主要集中在其微尺度噴管上，國(guó)內(nèi)外學(xué)者對(duì)微尺度下噴管內(nèi)粘性和熱損失效應(yīng)，稀薄效應(yīng)對(duì)噴管內(nèi)流場(chǎng)的影響、微噴管內(nèi)兩相流方面的研究開(kāi)展了較多的工作［1－6］，但對(duì)推力器燃燒室-噴管一體化結(jié)構(gòu)的推力器非穩(wěn)態(tài)工作過(guò)程研究相對(duì)較少。2011年，西班牙國(guó)家宇航中心的José A Morínigo等對(duì)其設(shè)計(jì)的一種新型固體微推力器結(jié)構(gòu)在近真空環(huán)境下進(jìn)行了非穩(wěn)態(tài)一體化仿真。結(jié)果表明，這種結(jié)構(gòu)可很好地減小熱損失，二階滑移邊界條件對(duì)推力幅值有45%的影響［7］。目前，對(duì)固體微推力器非穩(wěn)態(tài)工作情況的研究仍多以實(shí)驗(yàn)手段為主，但微型推力器工作過(guò)程時(shí)間通常較短(不到2 s)，且尺寸為亞毫米級(jí)，體積微小，測(cè)試難度大，難以直接獲得內(nèi)流場(chǎng)各物理量數(shù)據(jù)。因此，通過(guò)數(shù)值模擬手段對(duì)推力器一體化模型及其瞬態(tài)工作過(guò)程的研究是有必要的。

由于藥柱成型及裝填難度的限制，當(dāng)前國(guó)內(nèi)外的固體微推力器多采用端燃裝藥，由靠近推進(jìn)劑表面的點(diǎn)火電阻點(diǎn)燃推進(jìn)劑。其中，以法國(guó)國(guó)家科學(xué)研究中心(LAAS-CNRS)開(kāi)展的相關(guān)研究最具有代表性，本文主要基于LAAS-CNRS的平板結(jié)構(gòu)推力器模型進(jìn)行初步研究。

1 物理模型

1.1 推力器幾何模型

推力器蝕刻深度545 μm，燃燒室寬度1 500 μm，長(zhǎng)度5 300 μm，噴管入口寬度與燃燒室寬度相等，喉部寬度為150 μm，出口寬度為600 μm，燃燒室壁厚625 μm，擴(kuò)張半角 10°［8］，假定推力器壁面完全光滑。

本文研究的推力器蝕刻深度大于喉部寬度的3倍，因此可忽略三維壁面效應(yīng)［9］。同時(shí)，出于簡(jiǎn)化分析和計(jì)算成本的考慮，設(shè)置推力器在蝕刻方向上的壁面為絕熱壁面，采用二維CFD模型進(jìn)行分析。微型固體推力器幾何模型如圖1所示。

圖1 微型固體推力器幾何模型Fig.1 Geometric model of solid propellant micro-thruster

1.2 控制方程

1.2.1 流場(chǎng)區(qū)域控制方程

對(duì)微尺度下的流場(chǎng)，克努森數(shù)(Kn)是決定氣體流動(dòng)狀態(tài)的關(guān)鍵無(wú)量綱參數(shù)，其定義為分子平均自由程(λ)與特征幾何長(zhǎng)度(L)的比值，可表示為當(dāng)?shù)乩字Z數(shù)(Re＇)、氣體比熱容比γ和當(dāng)?shù)伛R赫數(shù)(Ma)的函數(shù)，其表示方式如下:

對(duì)于Kn＜0.14的微尺度流動(dòng)，一階速度滑移和溫度突躍邊界條件能準(zhǔn)確模擬微尺度效應(yīng)對(duì)流場(chǎng)的影響，表示如下［10］:

考慮熱蠕變的一階速度滑移邊界條件:

壁面溫度突躍:

式中 ug和Tg分別為緊鄰壁面處流體沿壁面的切向速度和溫度;uw和Tw分別為壁面移動(dòng)速度(本文中壁面靜止，uw=0)和溫度;σv為壁面動(dòng)量調(diào)節(jié)系數(shù);σT為熱調(diào)節(jié)系數(shù)。

對(duì)本文的情況，溫度突躍對(duì)流場(chǎng)的影響可以忽略［11］。因此，只考慮一階速度滑移，按完全漫反射模型取σv=1。

主流區(qū)控制方程仍然為通用N-S方程:

該通用方程包含了連續(xù)方程、動(dòng)量方程和能量方程。

燃?xì)鉂M足理想氣體狀態(tài)方程:

燃?xì)庹扯炔捎肧utherland方程描述:

式中 μ0為參考溫度T0=273.11 K時(shí)的動(dòng)力粘度，在此對(duì)應(yīng)的值為1.663×10－5kg/(m·s);S1為常數(shù)，其值為106.67 K。

微尺度下的層流到非充分發(fā)展湍流的轉(zhuǎn)捩雷諾數(shù)與常規(guī)尺度下表現(xiàn)出一定的差異，會(huì)發(fā)生轉(zhuǎn)捩提前的現(xiàn)象。根據(jù)國(guó)內(nèi)外微流動(dòng)元器件的部分實(shí)驗(yàn)結(jié)果，可近似以式(7)概括［12］:

式中 L為流動(dòng)方向上的長(zhǎng)度;Dh為水力直徑。對(duì)本文的模型，推力器的，即 Recr＜446.7。

作為初步計(jì)算，本文采用標(biāo)準(zhǔn)k-ε模型進(jìn)行流場(chǎng)仿真。求解器為隱式耦合非穩(wěn)態(tài)求解器，計(jì)算方法采用SIMPLE算法。

首先，對(duì)不考慮微尺度效應(yīng)的推力器流場(chǎng)進(jìn)行計(jì)算，得到克努森數(shù)分布后，再考慮是否對(duì)壁面添加滑移邊界條件。

1.2.2 固體區(qū)域控制方程

對(duì)固體區(qū)域僅考慮熱傳導(dǎo)，其控制方程為二維常物性、無(wú)內(nèi)熱源、非穩(wěn)態(tài)的導(dǎo)熱微分方程:

1.3 物性參數(shù)

推進(jìn)劑采用平臺(tái)雙基推進(jìn)劑SD，其主要參數(shù)如表1所示［13］。其中，燃速采用添加了少量黑火藥后的推進(jìn)劑在常壓環(huán)境下于微推力器內(nèi)通過(guò)實(shí)驗(yàn)完全燃燒后獲得的平均燃速［8］。藥柱參數(shù)參考同類(lèi)推進(jìn)劑相關(guān)參數(shù)，如表2所示。

壁面分別設(shè)置為文獻(xiàn)［8］中組成推力器殼體的2種材料進(jìn)行計(jì)算。由表3結(jié)合熱擴(kuò)散率的定義式計(jì)算得到二者的熱擴(kuò)散率分別為 8.6572 ×10－5、7.33 ×10－7m2/s。

表1 推進(jìn)劑參數(shù)Table 1 Parameters of the propellant

表2 藥柱參數(shù)Table 2 Parameters of the grain

表3 壁面材料參數(shù)Table 3 Parameters of the wall

1.4 計(jì)算區(qū)域、網(wǎng)格劃分及邊界條件

本文考慮的計(jì)算模型區(qū)域和推力器區(qū)域如圖2所示。計(jì)算網(wǎng)格采用結(jié)構(gòu)網(wǎng)格，圖3為計(jì)算網(wǎng)格(局部)。

圖2 計(jì)算區(qū)域Fig.2 Computational area

圖3 局部計(jì)算網(wǎng)格(推力器噴管及附近區(qū)域)Fig.3 Partial computational mesh(nozzle and the area nearby)

模型出口條件設(shè)定為壓力出口，壓強(qiáng)為101 325 Pa，計(jì)算域初始溫度為300 K。燃燒室外壁面設(shè)為絕熱壁面，根據(jù)第三類(lèi)邊界條件:

設(shè)置燃燒室、噴管內(nèi)壁面及推進(jìn)劑表面為熱耦合壁面。

假定推進(jìn)劑在燃燒表面完全燃燒，因此將燃?xì)庾鳛閱我粴怏w組分考慮，對(duì)推進(jìn)劑表面相鄰流場(chǎng)區(qū)域單元通過(guò)FLUENT用戶自定義函數(shù)UDF加載燃?xì)赓|(zhì)量源項(xiàng)、動(dòng)量源項(xiàng)與能量源項(xiàng)。其公式表示如下:

式中 r為推進(jìn)劑燃速;A為推進(jìn)劑燃面面積;ρg為燃?xì)饷芏?cp為燃?xì)獗榷▔簾崛?Tf為燃溫。

燃面退移過(guò)程通過(guò)動(dòng)網(wǎng)格進(jìn)行模擬，假定燃面按實(shí)測(cè)平均燃速4.6 mm/s平行退移，網(wǎng)格更新方法為動(dòng)態(tài)分層法。時(shí)間步長(zhǎng)為0.001 s，為保證計(jì)算準(zhǔn)確性，內(nèi)迭代步數(shù)為1 000。

本文推力的計(jì)算采用將壁面壓力對(duì)面積積分，再與粘性阻力矢量求和，最后取軸向分量的方法獲得。

2 計(jì)算結(jié)果

2.1 流動(dòng)微尺度效應(yīng)

首先，不考慮滑移邊界條件對(duì)微推力器工作過(guò)程進(jìn)行計(jì)算，得到推力器內(nèi)部基于喉部雷諾數(shù)的克努森數(shù)分布，推力器壁面材料選擇Macor Corning Ceramic。圖4為0.2 s和1.0 s時(shí)推力器內(nèi)流場(chǎng)克努森數(shù)分布情況。

圖4 0.2 s和1.0 s時(shí)刻推力器內(nèi)流場(chǎng)努森數(shù)分布Fig.4 Kudsen number distribution inside the thruster at 0.2 s and 1.0 s

可看出，推力器噴管內(nèi)部噴喉處克努森數(shù)最大，這與文獻(xiàn)［14］計(jì)算得到的噴管克努森數(shù)分布有所不同。其原因在于本文計(jì)算的推力器噴喉處流速未達(dá)到聲速，噴喉未處于壅塞狀態(tài)，噴管擴(kuò)張段壓力沿流動(dòng)方向升高導(dǎo)致克努森數(shù)下降。根據(jù)克努森數(shù)云圖，噴管內(nèi)克努森數(shù)均大于0.001，應(yīng)對(duì)噴管內(nèi)部流場(chǎng)考慮速度滑移邊界條件［11］。

圖5為0.1～1 s間噴喉處壁面和噴管出口處壁面滑移速度隨時(shí)間的變化規(guī)律?？煽吹?，噴喉壁面處的速度滑移現(xiàn)象較為明顯，其值分布于23～28 m/s之間;而噴管出口處速度滑移現(xiàn)象不明顯，其值分布于2～4 m/s之間。這與噴管內(nèi)克努森數(shù)的分布情況是相一致的，即克努森數(shù)越大，稀薄效應(yīng)越強(qiáng)，壁面滑移速度越大。從圖6所示的推力-時(shí)間曲線可看出，對(duì)本文的推力器及其工作條件，滑移邊界條件對(duì)推力的影響很小，對(duì)平衡段推力最大只有1.5%的差異，這與Ivanov等的結(jié)論相一致［2］。因此，在下面計(jì)算中忽略滑移邊界條件。

圖5 0.1～1 s間噴喉和噴管出口壁面滑移速度-時(shí)間曲線Fig.5 Wall slip velocity-time curves for nozzle throat and nozzle outlet between 0.1 s and 1s

圖6 考慮一階速度滑移邊界條件和無(wú)滑移邊界條件后的推力-時(shí)間曲線Fig.6 Thrust-time curves with and without 1st order velocity slip boundary condition

2.2 推力器瞬態(tài)溫度場(chǎng)

忽略流動(dòng)微尺度效應(yīng)后分別對(duì)壁面材料為Macor Corning Ceramic和Si的推力器進(jìn)行非穩(wěn)態(tài)仿真。

根據(jù)圖7、圖8可定性地看出，推進(jìn)劑表面出現(xiàn)了一定的溫度梯度，噴管內(nèi)流場(chǎng)區(qū)域溫度變化很小。采用Macor Corning Ceramic材料的推力器噴喉壁面附近存在明顯的高溫區(qū)。隨著燃燒過(guò)程的進(jìn)行，高溫區(qū)逐漸沿軸向、徑向擴(kuò)散，造成噴管收斂段和擴(kuò)張段壁面的溫度也有所上升;而采用Si材料的推力器噴喉壁面及附近區(qū)域溫升不明顯。

根據(jù)式(9)，在初始時(shí)刻，對(duì)于噴管壁面，2種材料推力器的壁面邊界層換熱量相差較小，而Si的熱導(dǎo)率遠(yuǎn)大于Macor Corning Ceramic的熱導(dǎo)率，因而Si壁面的表面溫度梯度較小。同時(shí)，Si的熱擴(kuò)散率遠(yuǎn)大于Macor Corning Ceramic的熱擴(kuò)散率，固體區(qū)域內(nèi)部溫度趨于一致的能力較強(qiáng)，這是導(dǎo)致上述現(xiàn)象的主要原因。

圖7 采用Macor Corning Ceramic材料的推力器內(nèi)部溫度云圖Fig.7 Temprature profile for thruster with Macor Corning Ceramic material

圖8 采用Si材料的推力器內(nèi)部溫度云圖Fig.8 Temprature profile for thruster with silicon material

在圖9中，上述壁面溫度分布曲線都有2個(gè)峰值，橫坐標(biāo)較小處的峰值表示與推進(jìn)劑燃面接觸的壁面溫度，橫坐標(biāo)較大處的峰值表示噴喉處壁面溫度。兩峰值間的拐點(diǎn)表示噴管入口處壁面溫度。

圖9 不同時(shí)刻推力器內(nèi)壁面溫度分布曲線Fig.9 Temperature distribution curves at different moments

由圖9可見(jiàn)，對(duì)于壁面采用Macor Corning Ceramic材料的推力器，噴管內(nèi)壁面升溫很快，在0.05 s內(nèi)溫度就由300 K上升至約1 000 K。隨后，升溫趨勢(shì)隨時(shí)間減小，最后溫度趨于穩(wěn)定。燃燒室內(nèi)溫度起初在燃面與噴管入口間存在較大的軸向梯度，隨著燃燒的進(jìn)行，壁面被不斷加熱，溫度曲線趨于平緩。而對(duì)于壁面采用Si材料的推力器，噴管壁面升溫較慢，高溫區(qū)不明顯。由于其較大的熱擴(kuò)散率，推進(jìn)劑燃面位置之前的壁面溫度也有所上升。

2.3 流動(dòng)損失

在推進(jìn)劑穩(wěn)定燃燒的情況下，固體微推力器的推力損失主要由壁面熱損失和粘性損失構(gòu)成。圖10為2種壁面材料推力器的推力-時(shí)間曲線。

圖10 2種不同壁面材料推力器的推力-時(shí)間曲線Fig.10 Thrust-time curves for thrusters with different wall material

由圖10可見(jiàn)，采用Macor Corning Ceramic材料的推力器能夠維持推力上升至平衡段，而采用Si材料的推力器無(wú)法維持推力上升。

2.3.1 熱損失

推力器內(nèi)的總體熱損失由通過(guò)燃燒室壁面的熱損失和通過(guò)噴管壁面的熱損失組成。圖11為2種不同壁面材料推力器的熱損失-時(shí)間曲線?？梢?jiàn)，在數(shù)值上，采用Macor Corning Ceramic材料的推力器熱損失明顯小于采用Si材料的推力器熱損失。同時(shí)，可看到采用Macor Corning Ceramic材料和采用Si材料的推力器，其熱損失隨時(shí)間的變化趨勢(shì)具有一定共性:通過(guò)噴管壁面的熱損失在工作過(guò)程中主要表現(xiàn)為下降趨勢(shì)，而通過(guò)燃燒室壁面的熱損失在工作過(guò)程中主要表現(xiàn)為上升趨勢(shì)。結(jié)合上節(jié)的溫度場(chǎng)分析，可認(rèn)為其原因在于:流體到固體壁面的換熱主要在表面邊界層中進(jìn)行，噴喉處粘性邊界層較薄，表面溫度梯度相對(duì)較大，流速最快，單位面積換熱最劇烈，因而噴喉附近壁面溫度上升較快。隨著表面溫度的上升，流場(chǎng)與壁面間溫度梯度下降，因而通過(guò)噴管壁面的熱損失表現(xiàn)為下降趨勢(shì)。而隨著推進(jìn)劑燃面的不斷退移，更多的推力器內(nèi)壁面暴露于燃?xì)猸h(huán)境中，且由圖9可看出，燃燒室表面溫度上升緩慢。因此，通過(guò)燃燒室壁面的熱損失不斷上升。

圖11 2種不同壁面材料推力器的熱損失-時(shí)間曲線Fig.11 Heat loss-time curves for thrusters with different wall materials

但是，采用2種不同壁面材料的固體推力器內(nèi)的總體熱損失卻表現(xiàn)出相反的變化趨勢(shì):采用Macor Corning Ceramic材料的推力器總體熱損失在下降了約0.25 s后趨于穩(wěn)定，降幅約為25%;采用Si材料的推力器總體熱損失卻一直上升，增幅約為20%。這說(shuō)明前者噴管壁面熱損失的下降大于燃燒室壁面熱損失的升高，而后者噴管壁面熱損失的下降小于燃燒室壁面熱損失的升高。二者總體熱損失的變化趨勢(shì)恰好與二者的推力曲線變化趨勢(shì)相反。

2.3.2 粘性損失

通常由粘性阻力的變化趨勢(shì)衡量粘性損失的變化趨勢(shì)。圖12為2種材料推力器的粘性阻力-時(shí)間曲線。由圖12可知，壁面采用Si材料的微推力器噴管壁面粘性阻力數(shù)值相對(duì)較小，根據(jù)文獻(xiàn)［5］，傳熱壁面的粘性阻力小于絕熱壁面，結(jié)合圖11熱損失的計(jì)算結(jié)果，發(fā)現(xiàn)Si壁面熱損失較大，粘性阻力較小;Macor Corning Ceramic壁面熱損失較小，粘性阻力較大?？蛇M(jìn)一步認(rèn)為:粘性阻力幅值與壁面熱損失呈負(fù)相關(guān)性。

圖12 2種不同壁面材料推力器的粘性阻力-時(shí)間曲線Fig.12 Viscous force-time curves for thrusters with different wall materials

噴管壁面粘性阻力首先經(jīng)歷短暫的上升，其主要原因是噴管內(nèi)部速度場(chǎng)和溫度場(chǎng)的建立;隨后，在剩余工作時(shí)間內(nèi)粘性阻力處于緩慢下降階段，其主要原因是流場(chǎng)穩(wěn)定后燃?xì)馔诿娴膫鳠帷?/p>

燃燒室壁面粘性阻力緩慢上升，這主要是由于燃面退移造成更大面積的壁面暴露于燃?xì)饬鲌?chǎng)中所導(dǎo)致的，但由于燃燒室內(nèi)燃?xì)饬魉佥^低，因此其幅值不大?？煽吹剑傮w粘性阻力的變化趨勢(shì)主要由噴管壁面粘性阻力的變化趨勢(shì)決定。

2.3.3 流動(dòng)損失綜合分析

總結(jié)粘性損失、熱損失和推力的變化趨勢(shì)?？傻贸?對(duì)本文研究的推力器，在推進(jìn)劑按實(shí)測(cè)燃速穩(wěn)定燃燒的前提下，推力隨時(shí)間變化趨勢(shì)的主要影響因素是總體熱損失隨時(shí)間的變化趨勢(shì)，即推力器噴管壁面的熱損失與燃燒室壁面熱損失之和。對(duì)Macor Corning Ceramic材料，由于其熱擴(kuò)散率小，壁面附近升溫明顯，有效降低了與燃?xì)忾g的溫度梯度，減小了噴管壁面熱損失，其降幅超過(guò)了燃燒室壁面熱損失的上升幅值，使推力可上升并趨于穩(wěn)定。對(duì)Si材料，由于其熱擴(kuò)散率大，壁面附近升溫緩慢，噴管壁面熱損失的下降無(wú)法抵消燃燒室壁面熱損失的上升，造成推力持續(xù)下降，無(wú)法維持。壁面采用Macor Corning Ceramic材料不僅從數(shù)值上減小了壁面熱損失，而且使壁面熱損失隨時(shí)間減小，抵消了燃燒室壁面熱損失的上升，推力可上升并維持。文獻(xiàn)［15］中的推力器實(shí)驗(yàn)結(jié)果表明:在提高壁面等效熱導(dǎo)率前，推力器工作過(guò)程無(wú)法維持，而提高壁面等效熱導(dǎo)率后，推力可上升至平衡段，這與本文的計(jì)算結(jié)果是相似的。

可看到，微型固體推力器的瞬態(tài)流動(dòng)損失是一個(gè)復(fù)雜的過(guò)程，其影響因素主要有壁面材料、燃?xì)鉁囟?、粘性和推進(jìn)劑燃燒模型。本文作為初步研究，假定推進(jìn)劑完全燃燒，忽略了推進(jìn)劑配方改變導(dǎo)致的燃燒模型變化以及輻射傳熱等因素，該問(wèn)題仍有待進(jìn)一步研究。

3 結(jié)論

(1)通過(guò)所建立的端燃裝藥固體微推力器瞬態(tài)工作過(guò)程的流固耦合仿真模型，實(shí)現(xiàn)了對(duì)推力器內(nèi)流場(chǎng)、壁面及推進(jìn)劑內(nèi)部溫度場(chǎng)的一體化仿真。

(2)對(duì)本文所涉及的工作條件，推力器噴管內(nèi)部克努森數(shù)雖然使噴喉壁面流場(chǎng)出現(xiàn)較明顯的滑移現(xiàn)象，但噴喉處未達(dá)到壅塞狀態(tài)，擴(kuò)張段內(nèi)壓力升高，使得擴(kuò)張段內(nèi)克努森數(shù)降低。因此，擴(kuò)張段內(nèi)滑移現(xiàn)象不明顯，對(duì)推力的影響很小，可忽略不計(jì)。

(3)采用Si材料的壁面熱導(dǎo)率較高，因而在相近的初始傳熱量下壁面溫度梯度低，且其熱擴(kuò)散率高，導(dǎo)致與燃?xì)饨佑|的噴管壁面升溫緩慢，溫度梯度減小緩慢。而采用Macor Corning Ceramic材料的壁面熱導(dǎo)率較低，因而在相近的初始傳熱量下壁面溫度梯度高，且其熱擴(kuò)散率較低，導(dǎo)致壁面附近區(qū)域溫度升高較快，有效降低了溫度梯度。

(4)粘性阻力的幅值與壁面散熱量呈負(fù)相關(guān)性，總體粘性阻力的變化趨勢(shì)由噴管壁面粘性阻力變化趨勢(shì)主導(dǎo)。

(5)對(duì)本文研究的固體微推力器，在推進(jìn)劑按實(shí)測(cè)燃速穩(wěn)定燃燒的前提下，影響推力隨時(shí)間變化趨勢(shì)的主要因素是壁面熱損失隨時(shí)間的變化趨勢(shì)，推力能否上升并維持的關(guān)鍵在于噴管壁面熱損失的減小能否抵消由于燃面退移造成的燃燒室壁面熱損失的增大。

［1］Markelov G N，Ivanov M S.Numerical study of 2D/3D micronozzle flows［C］//Proceedings of 22nd Symposium on Rarefied Gas Dynamics.2001.

［2］Ivanov M S，Markelov G N，Ketsdever A D，et al.Numerical study of cold gas micronozzle flow ［C］//37th Aerospace Sciences Meeting and Exhibit.1999.

［3］Alexeenko A A，Levin D A.Numerical simulation of hightemperature gas flows in a millimeter-scale thruster［J］.Journal of Thermophysics and Heat Transfer，2002，16(1).

［4］張斌，毛根旺，胡松啟，等.固體微推力器氣粒兩相羽流場(chǎng)的數(shù)值模擬［J］.固體火箭技術(shù)，2011，34(3).

［5］吳素麗，胡松啟，張斌，等.MEMS-SPMT噴管內(nèi)的粘性和熱損失研究［J］.固體火箭技術(shù)，2012，35(1).

［6］夏廣慶，張斌，孫得川，等.微型固體姿控發(fā)動(dòng)機(jī)微噴管內(nèi)氣粒兩相流動(dòng)規(guī)律的CFD-DSMC研究［J］.固體火箭技術(shù)，2012，35(3).

［8］Chaalane A，Rossi C，Esteve D.The formulation and testing of new solid propellant mixture(DB+x%BP)for a new MEMS-based microthruster［J］.Sensors and Actuators，A 138(2007).

［9］楊海威，朱衛(wèi)兵，趙陽(yáng).拉伐爾型微噴管流場(chǎng)的三維模擬［J］.半導(dǎo)體技術(shù)，2009，34(10).

［10］Von Smoluchowski M.Ueber Würmeleitung in verdunnten gasen［J］.Annalen der Physik und Chemie 1898，64.

［11］Zhang Genxuan，Liu Minghou，Zhang Xianfeng，et al.Continuum based slip model and its validity for micro-channel flows［J］.Chinese Science Bulletin，2006，51(9).

［12］計(jì)光華，計(jì)洪苗.微流動(dòng)及其元器件［M］.北京:高等教育出版社，2009.

［13］Gilles Fonblanc，Bruno Herran.A new insensitive，non-toxic double base propellant for rocket motors［C］//30th AIAA/ASME/SAE/ASEE Joint Propulsion Conference.2004.

［14］孫建威，劉明侯，張根烜，等.滑移邊界條件對(duì)二維微噴管數(shù)值模擬結(jié)果的影響［J］.宇航學(xué)報(bào)，2005，26(6).

［15］Rossi C，Benoit Larangot，Denis Lagrange，et al.Final characterizations of MEMS-based pyrotechnical microthrusters［J］.Sensors and Actuators，2005，A 121.

［16］Kaili Zhang，Chou S K，Simon S Ang.MEMS-based solid propellant microthruster design，simulation，fabrication，and testing［J］.Journal of Microelectro-Mechanical Systems，2004，13(2).