微型固體推力器瞬態工作過程數值模擬①

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

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

0 引言

基于MEMS技術制造的固體微推力器具有結構簡單、功耗低、無可動部件及其帶來的摩擦損失小、可靠性高等諸多優勢。固體微推力器最大的優勢是集成度高，可集成在一塊基片上組成陣列。目前，陣列式推進系統已在動能攔截器等方面投入應用，如美國的增程攔截彈ERINT-1。基于上述特點，固體微推力器成為微推進系統方面新的研究熱點。

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

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

1 物理模型

1.1 推力器幾何模型

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

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

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

1.2 控制方程

1.2.1 流場區域控制方程

對微尺度下的流場，克努森數(Kn)是決定氣體流動狀態的關鍵無量綱參數，其定義為分子平均自由程(λ)與特征幾何長度(L)的比值，可表示為當地雷諾數(Re＇)、氣體比熱容比γ和當地馬赫數(Ma)的函數，其表示方式如下:

對于Kn＜0.14的微尺度流動，一階速度滑移和溫度突躍邊界條件能準確模擬微尺度效應對流場的影響，表示如下［10］:

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

壁面溫度突躍:

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

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

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

該通用方程包含了連續方程、動量方程和能量方程。

燃氣滿足理想氣體狀態方程:

燃氣粘度采用Sutherland方程描述:

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

微尺度下的層流到非充分發展湍流的轉捩雷諾數與常規尺度下表現出一定的差異，會發生轉捩提前的現象。根據國內外微流動元器件的部分實驗結果，可近似以式(7)概括［12］:

式中 L為流動方向上的長度;Dh為水力直徑。對本文的模型，推力器的，即 Recr＜446.7。

作為初步計算，本文采用標準k-ε模型進行流場仿真。求解器為隱式耦合非穩態求解器，計算方法采用SIMPLE算法。

首先，對不考慮微尺度效應的推力器流場進行計算，得到克努森數分布后，再考慮是否對壁面添加滑移邊界條件。

1.2.2 固體區域控制方程

對固體區域僅考慮熱傳導，其控制方程為二維常物性、無內熱源、非穩態的導熱微分方程:

1.3 物性參數

推進劑采用平臺雙基推進劑SD，其主要參數如表1所示［13］。其中，燃速采用添加了少量黑火藥后的推進劑在常壓環境下于微推力器內通過實驗完全燃燒后獲得的平均燃速［8］。藥柱參數參考同類推進劑相關參數，如表2所示。

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

表1 推進劑參數Table 1 Parameters of the propellant

表2 藥柱參數Table 2 Parameters of the grain

表3 壁面材料參數Table 3 Parameters of the wall

1.4 計算區域、網格劃分及邊界條件

本文考慮的計算模型區域和推力器區域如圖2所示。計算網格采用結構網格，圖3為計算網格(局部)。

圖2 計算區域Fig.2 Computational area

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

模型出口條件設定為壓力出口，壓強為101 325 Pa，計算域初始溫度為300 K。燃燒室外壁面設為絕熱壁面，根據第三類邊界條件:

設置燃燒室、噴管內壁面及推進劑表面為熱耦合壁面。

假定推進劑在燃燒表面完全燃燒，因此將燃氣作為單一氣體組分考慮，對推進劑表面相鄰流場區域單元通過FLUENT用戶自定義函數UDF加載燃氣質量源項、動量源項與能量源項。其公式表示如下:

式中 r為推進劑燃速;A為推進劑燃面面積;ρg為燃氣密度;cp為燃氣比定壓熱容;Tf為燃溫。

燃面退移過程通過動網格進行模擬，假定燃面按實測平均燃速4.6 mm/s平行退移，網格更新方法為動態分層法。時間步長為0.001 s，為保證計算準確性，內迭代步數為1 000。

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

2 計算結果

2.1 流動微尺度效應

首先，不考慮滑移邊界條件對微推力器工作過程進行計算，得到推力器內部基于喉部雷諾數的克努森數分布，推力器壁面材料選擇Macor Corning Ceramic。圖4為0.2 s和1.0 s時推力器內流場克努森數分布情況。

圖4 0.2 s和1.0 s時刻推力器內流場努森數分布Fig.4 Kudsen number distribution inside the thruster at 0.2 s and 1.0 s

可看出，推力器噴管內部噴喉處克努森數最大，這與文獻［14］計算得到的噴管克努森數分布有所不同。其原因在于本文計算的推力器噴喉處流速未達到聲速，噴喉未處于壅塞狀態，噴管擴張段壓力沿流動方向升高導致克努森數下降。根據克努森數云圖，噴管內克努森數均大于0.001，應對噴管內部流場考慮速度滑移邊界條件［11］。

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

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

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

2.2 推力器瞬態溫度場

忽略流動微尺度效應后分別對壁面材料為Macor Corning Ceramic和Si的推力器進行非穩態仿真。

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

根據式(9)，在初始時刻，對于噴管壁面，2種材料推力器的壁面邊界層換熱量相差較小，而Si的熱導率遠大于Macor Corning Ceramic的熱導率，因而Si壁面的表面溫度梯度較小。同時，Si的熱擴散率遠大于Macor Corning Ceramic的熱擴散率，固體區域內部溫度趨于一致的能力較強，這是導致上述現象的主要原因。

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

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

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

圖9 不同時刻推力器內壁面溫度分布曲線Fig.9 Temperature distribution curves at different moments

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

2.3 流動損失

在推進劑穩定燃燒的情況下，固體微推力器的推力損失主要由壁面熱損失和粘性損失構成。圖10為2種壁面材料推力器的推力-時間曲線。

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

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

2.3.1 熱損失

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

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

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

2.3.2 粘性損失

通常由粘性阻力的變化趨勢衡量粘性損失的變化趨勢。圖12為2種材料推力器的粘性阻力-時間曲線。由圖12可知，壁面采用Si材料的微推力器噴管壁面粘性阻力數值相對較小，根據文獻［5］，傳熱壁面的粘性阻力小于絕熱壁面，結合圖11熱損失的計算結果，發現Si壁面熱損失較大，粘性阻力較小;Macor Corning Ceramic壁面熱損失較小，粘性阻力較大。可進一步認為:粘性阻力幅值與壁面熱損失呈負相關性。

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

噴管壁面粘性阻力首先經歷短暫的上升，其主要原因是噴管內部速度場和溫度場的建立;隨后，在剩余工作時間內粘性阻力處于緩慢下降階段，其主要原因是流場穩定后燃氣往壁面的傳熱。

燃燒室壁面粘性阻力緩慢上升，這主要是由于燃面退移造成更大面積的壁面暴露于燃氣流場中所導致的，但由于燃燒室內燃氣流速較低，因此其幅值不大。可看到，總體粘性阻力的變化趨勢主要由噴管壁面粘性阻力的變化趨勢決定。

2.3.3 流動損失綜合分析

總結粘性損失、熱損失和推力的變化趨勢。可得出:對本文研究的推力器，在推進劑按實測燃速穩定燃燒的前提下，推力隨時間變化趨勢的主要影響因素是總體熱損失隨時間的變化趨勢，即推力器噴管壁面的熱損失與燃燒室壁面熱損失之和。對Macor Corning Ceramic材料，由于其熱擴散率小，壁面附近升溫明顯，有效降低了與燃氣間的溫度梯度，減小了噴管壁面熱損失，其降幅超過了燃燒室壁面熱損失的上升幅值，使推力可上升并趨于穩定。對Si材料，由于其熱擴散率大，壁面附近升溫緩慢，噴管壁面熱損失的下降無法抵消燃燒室壁面熱損失的上升，造成推力持續下降，無法維持。壁面采用Macor Corning Ceramic材料不僅從數值上減小了壁面熱損失，而且使壁面熱損失隨時間減小，抵消了燃燒室壁面熱損失的上升，推力可上升并維持。文獻［15］中的推力器實驗結果表明:在提高壁面等效熱導率前，推力器工作過程無法維持，而提高壁面等效熱導率后，推力可上升至平衡段，這與本文的計算結果是相似的。

可看到，微型固體推力器的瞬態流動損失是一個復雜的過程，其影響因素主要有壁面材料、燃氣溫度、粘性和推進劑燃燒模型。本文作為初步研究，假定推進劑完全燃燒，忽略了推進劑配方改變導致的燃燒模型變化以及輻射傳熱等因素，該問題仍有待進一步研究。

3 結論

(1)通過所建立的端燃裝藥固體微推力器瞬態工作過程的流固耦合仿真模型，實現了對推力器內流場、壁面及推進劑內部溫度場的一體化仿真。

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

(3)采用Si材料的壁面熱導率較高，因而在相近的初始傳熱量下壁面溫度梯度低，且其熱擴散率高，導致與燃氣接觸的噴管壁面升溫緩慢，溫度梯度減小緩慢。而采用Macor Corning Ceramic材料的壁面熱導率較低，因而在相近的初始傳熱量下壁面溫度梯度高，且其熱擴散率較低，導致壁面附近區域溫度升高較快，有效降低了溫度梯度。

(4)粘性阻力的幅值與壁面散熱量呈負相關性，總體粘性阻力的變化趨勢由噴管壁面粘性阻力變化趨勢主導。

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

［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］張斌，毛根旺，胡松啟，等.固體微推力器氣粒兩相羽流場的數值模擬［J］.固體火箭技術，2011，34(3).

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

［6］夏廣慶，張斌，孫得川，等.微型固體姿控發動機微噴管內氣粒兩相流動規律的CFD-DSMC研究［J］.固體火箭技術，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］楊海威，朱衛兵，趙陽.拉伐爾型微噴管流場的三維模擬［J］.半導體技術，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］計光華，計洪苗.微流動及其元器件［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］孫建威，劉明侯，張根烜，等.滑移邊界條件對二維微噴管數值模擬結果的影響［J］.宇航學報，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).