吳素麗，胡松啟，張 斌，劉迎吉，劉 凱，王鵬飛

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

0 引言

基于MEMS的固體微推力器(MEMS－based solid propellant microthruster，MEMS－SPMT)尺寸小，微沖量精確度高，集成性好，在微小衛(wèi)星姿軌控技術中具有廣闊的應用前景［1－2］。隨著微推力器的火箭發(fā)動機尺寸的縮小，其面積容積比增大，微噴管內(nèi)的雷諾數(shù)特別低(一般低于 1 000)，粘性耗散特別嚴重［3－4］;同時，噴管殼體材料通常選用硅(熱導率較高)，高溫氣體向壁面的熱傳導嚴重［5－6］。因此，需研究微噴管內(nèi)的換熱和粘性效應對微噴管推力性能的影響。

國內(nèi)外關于微推力器噴管內(nèi)的粘性耗散以及熱量損失做了很多的研究工作。Bayt和Breuer最早研究了超聲速微噴管內(nèi)的流動狀態(tài)，發(fā)現(xiàn)粘性亞聲速層對流場的作用較大，推力損失較大［7－8］。2003年，Kujawa J和Hitt D L等用CFD軟件研究了NASA/GSFC的雙組元發(fā)動機噴管的性能，得出的結論是降低固壁溫度，能降低噴管12%的推力，較大的噴管擴張角能減小粘性作用［9］。2007年，Louisos W F和 Hitt D L研究了2D和3D超聲速微噴管內(nèi)的粘性和熱效應對噴管性能的影響。結果表明，在三維模擬中，粘性作用更加明顯，流動產(chǎn)生的熱損失能降低粘性作用，提高微噴管的性能［3，5］。2002年，國內(nèi)李德桃等研究發(fā)現(xiàn)，通過燃燒室壁面的熱損失很大，出口溫度降低了13%［10］。2009年，楊海威等采用直接蒙特卡羅法，建立微噴管流場和壁面溫度場耦合仿真模型，研究了以壁面固定溫度加載和以時間更新的熱流密度加載2種工況下微噴管的熱損失［11］。

本文從湍流模型、噴管壁面初始溫度和噴管構型3方面研究噴管內(nèi)的粘性和換熱損失，利用多參數(shù)分別評價亞聲速邊界層造成的粘性損失和高溫燃氣向固體壁面的傳熱損失對微噴管性能的影響。

1 計算模型

1.1 物理模型

3種微型發(fā)動機噴管結構及網(wǎng)格劃分見圖1。圖1尺寸單位為μm。

圖1中，微噴管截面為正方形結構，Lt代表微噴管喉部邊長，微噴管的收斂和擴張半角均為35.2°。燃燒室壓強設為1 MPa，溫度為3 000 K，噴管出口壓強為1 kPa。噴管材料選用硅，導熱系數(shù)為150 W/(m·K)，密度為 2 330 kg/m3，比熱容為700 J/(kg·K)。推進劑采用1∶1混合的斯蒂芬酸鉛和硝化棉，利用最小吉布斯自由能法計算燃氣的熱力學參數(shù)。

此方法是利用體系在等壓條件下達到平衡狀態(tài)時，體系吉布斯自由能具有最小值，通過求解平衡方程組得到平衡組分的熱力學參數(shù)［12］。流場計算中，假定燃氣為理想氣體，燃氣的導熱系數(shù)為0.710 35 W/(m·K)，密度為 1.253 kg/m3，比熱容為2 715 J/(kg·K)，粘性系數(shù)為1.008 6 × 10－4Pa·s。

1.2 數(shù)學模型

1.2.1 湍流模型

S－A模型是針對大網(wǎng)格的簡單湍流模型，適于求解有壁面影響流動及有逆壓力梯度下的邊界層流動。低Re數(shù)k－ε模型是標準k－ε模型修正的方案，可自動適應不同Re數(shù)的區(qū)域，適于求解微尺度下噴管內(nèi)的流動。

三維微噴管內(nèi)的流動是由燃氣的質(zhì)量、動量和能量守恒方程控制的，具體可見文獻［3］。低Re數(shù)k－ε模型的輸運方程:

式中 n代表壁面法向坐標;u為與壁面平行的流速;系數(shù) G1ε、G2ε、Cμ、σk、σε及產(chǎn)生項 Gk與標準 k－ε 模型的方程相同。

三維常物性、無內(nèi)熱源、非穩(wěn)態(tài)的導熱微分方程(即拉普拉斯方程)為

微噴管內(nèi)流場計算時，已經(jīng)計算出噴管壁面處溫度分布，故微噴管的內(nèi)壁邊界條件設定為第三類邊界條件

式中 λ為導熱系數(shù);ρ為材料密度;c為定壓比熱容;h為邊界上物體與周圍流體的表面換熱系數(shù);tw為壁面溫度;tf為周圍流體溫度。

根據(jù)熱力計算結果，估算Lt=106 μm的微噴管喉部和出口的Re數(shù)分別是792和135。因此，本文選用這兩種低Re數(shù)湍流模型，研究不同的湍流模型在計算噴管內(nèi)粘性效應的差異。

1.2.2 邊界條件與初始化

采用隱式耦合非穩(wěn)態(tài)求解器，入口處給定來流總溫、總壓和方向角，出口采用外推邊界，壁面無滑移，流體向壁面的熱傳導采用流－固耦合邊界條件。為節(jié)省計算時間，只模擬噴管1/4型面流場。2個對稱面，2個傳熱壁面，邊界條件見圖2。

2 結果分析

采用高速攝影系統(tǒng)，測得微型固體火箭發(fā)動機的工作過程約1 ms，故本文分析發(fā)動機開始工作到工作結束(總共1 ms)壁面溫度的變化。

2.1 無傳熱時噴管的粘性損失

選用Lt=106 μm的噴管構型，壁面絕熱時，壁面初始溫度為300 K。S－A、低Re數(shù)k－ε模型下微噴管內(nèi)亞聲速輪廓見圖3。

圖4和圖5分別是2種模型下噴管中心線上沿z軸的速度分布和垂直對稱壁面上噴管內(nèi)輪廓線沿z軸的溫度分布。

由圖3計算出S－A模型下噴管出口截面的亞聲速邊界層(Ma＜1)面積占噴管出口面積的30.3%，亞聲速邊界層厚度遠大于常規(guī)尺寸噴管的邊界層，燃氣的粘性效應不能忽略。由圖5看出，2種模型計算出的噴管擴張段壁面溫度不同，S－A模型高于低Re數(shù)k－ε模型。其原因是2種湍流模型采用不同的湍動粘度的計算式，S－A模型中的輸運變量在近壁處的梯度要比低Re數(shù)k－ε模型中的小，這使得該模型對網(wǎng)格粗糙帶來數(shù)值誤差不太敏感，故S－A模型計算出的燃氣粘性較大。比較圖3中(a)、(b)看出，S－A模型的亞聲速層厚度大于低Re數(shù)k－ε模型。因此，S－A模型的計算結果是在擴張段速度滯止較大，壁面溫度升高。兩模型計算的主流速度是吻合的，差別僅在于靠近噴管壁面附近的粘性計算上。

2.2 存在傳熱時噴管的粘性和熱損失

選用Lt=106 μm的噴管構型，設定噴管壁面為流－固耦合壁面。1 ms時，低Re數(shù)k－ε模型下噴管選定截面的壁面溫度輪廓見圖6。

圖7是2種模型下微噴管水平對稱面與出口平面交線沿x軸馬赫數(shù)曲線。

圖8是2種模型下噴管垂直對稱面上噴管內(nèi)輪廓壁面沿z軸的溫度曲線。

圖9是2種模型下固體壁面的垂直對稱面上入口、喉部及出口這3條輪廓線(圖6中標注的3條輪廓線:入口、喉部和出口)沿y軸方向的溫度曲線。

圖10是低Re數(shù)k－ε模型下，固體壁面的入口和喉部上平分這2條線的6個點(圖6中的點1～6)隨時間的變化曲線。

結合圖6和圖9可看出，高溫燃氣向壁面的熱傳導過程，S－A模型計算的噴管壁面沿y軸的溫度在喉部之前都小于低Re數(shù)k－ε模型，喉部之后，二者相差不大，圖7可解釋這一現(xiàn)象。低Re數(shù)k－ε模型下噴管入口處壁面最低溫度為662 K，喉部壁面最低溫度為609 K，出口壁面最低溫度為410 K，3處溫度分別是初始壁溫的2.21、2.03、1.37 倍，溫度損失較大。

由圖10看出，壁面某點處溫度隨時間的變化基本是線性的，點2和點5在1 ms的時間內(nèi)溫度分別上升了278、240 K，溫升明顯。點2與點3、點4與點5的溫度變化曲線幾乎重合，但與壁面處點的溫度曲線有段距離，說明靠近壁面處溫度梯度較大，遠離壁面時溫差減小，傳熱減弱。

2.3 粘性和熱損失分析

表1列出了在無傳熱和有傳熱時，2種湍流模型下噴管的內(nèi)推力Fp、粘性阻力 Fμ、噴管比沖Isp、粘性與總的內(nèi)推力之比Fμ/FA(FA=Fp－Fμ)、噴管出口亞聲速輪廓與噴管出口面積比ε及燃氣與壁面間的導熱量q參數(shù)。由表1中看出，無熱損失時，噴管出口截面的亞聲速層面積占整個出口截面面積的23.1% ～30.6%，噴 管 的 粘 性 阻 力 占 總 推 力 的22.5% ～25.7%，說明粘性作用會顯著影響微噴管推力性能。

通過比較有無傳熱時的ε值可發(fā)現(xiàn)，對于同一種模型計算的ε值，有傳熱較無傳熱降低，說明氣體向固壁的傳熱造成了噴管壁面的粘性亞聲速層尺寸減小，噴管的粘性損失降低，在一定程度上能提高微噴管性能。因此，降低粘性損失和散熱損失的一個途徑是預加熱噴管壁面。

表1 4種情況下噴管流場參數(shù)(1/4噴管)Table 1 Calculation results of the four situations(1/4 of the nozzle)

有傳熱時，噴管的內(nèi)推力和比沖較無傳熱時小，這是由于高溫燃氣向壁面導熱損失了部分熱量。在無傳熱時，這部分熱量用來膨脹做功產(chǎn)生推力，因此微噴管性能降低。總體來說，高溫燃氣向壁面的熱傳導既減小亞聲速邊界層尺寸，一定程度上能提高微噴管性能，又損失部分能量降低微噴管性能。因此，S－A模型計算的由于傳熱造成比沖的損失為17.88%，標準k－ε模型計算的傳熱損失為17.04%，傳熱損失顯著。2種模型計算的噴管壁面的導熱量相差不大。

2.4 壁面初始溫度對傳熱態(tài)下噴管的粘性和熱損失影響

根據(jù)納衛(wèi)星的工作高度［13］，設定噴管壁面溫度為179 K和239 K，選用低Re數(shù)k－ε模型模擬不同壁面溫度條件下微噴管的粘性和換熱損失。

圖11是微噴管出口截面與水平對稱面的交線沿x軸的壓強曲線。

圖12是微噴管中心軸線沿z軸的馬赫數(shù)曲線。

圖13是噴管壁面的內(nèi)輪廓線沿z軸的溫度曲線。

圖14是圖6中標注的3條輪廓線(入口、喉部和出口)沿y軸方向的溫度曲線。由各圖可看出，壁面溫度為239 K和179 K時噴管內(nèi)流場基本相同，與壁面為常溫時有微小區(qū)別。

據(jù)表2，噴管壁面溫度升高時，由于壁面與燃氣的溫差降低，故燃氣向壁面的傳熱量減小，熱損失降低，對提高微噴管比沖是正效應，而壁面溫度的升高使得粘性力占總推力的比例增大，粘性亞聲速邊界層面積增大，粘性損失增加，對提高微噴管比沖是負效應，熱損失和粘性力二者正負作用結果是比沖小幅增大(最多提高64 N·s/kg)。

表2 3種初始溫度下的噴管流場參數(shù)(1/4噴管)Table 2 Calculation results of three initial temperture(1/4 of the nozzle)

2.5 3種噴管構型傳熱態(tài)下粘性和熱損失分析

選定3種本研究所用且已加工出的噴管構型(具體尺寸如圖1中(a)～(c)所示)，研究噴管內(nèi)流場，預測其能量水平。

圖15是微噴管出口截面與水平對稱面的交線沿x軸的馬赫數(shù)曲線。

圖16是微噴管中心軸線沿z軸的馬赫數(shù)曲線。

圖17是噴管壁面的內(nèi)輪廓線沿z軸方向的溫度曲線。

圖18是圖6中標注的3條輪廓線(入口、喉部和出口)沿y軸方向的溫度曲線。

表3列出了噴管的流場參數(shù)，其中A代表微噴管2個傳熱壁面的面積和。

表3 3種尺寸的噴管流場參數(shù)(1/4噴管)Table 3 Calculation results of three sizes of the nozzle(1/4 of the nozzle)

結合圖和表數(shù)據(jù)看出，噴管構型的變化使得微噴管的傳熱壁面面積發(fā)生變化，噴管喉部截面邊長從106 μm到178 μm，傳熱面積逐漸增大，高溫燃氣向壁面的傳熱量增大了57.4%，是負效應。微噴管壁面溫度升高，壁面粘性力小幅增大，但噴管的內(nèi)推力顯著增大。因此，粘性力占總推力的比例由22.6%降至11.8%，是正效應。正負效應綜合作用下，微噴管比沖提高64 N·s/kg。因此，增大微噴管尺寸能降低粘性效應，提高微噴管性能指標。

3 結論

(1)2種模型模擬微噴管流場結果相差不大，S－A模型計算的粘性效應比低Re數(shù)k－ε模型大。

(2)不考慮熱損失時，采用低Re數(shù)k－ε模型計算的噴管出口截面的亞聲速層面積占噴管出口面積的23.1%，粘性力占噴管推力25.7%，粘性損失明顯。

(3)熱傳導減小粘性亞聲速層尺寸，一定程度上通過降低粘性損失提高微噴管性能。另外，向壁面導熱會損失高溫燃氣能量，總體上微噴管比沖降低17%左右。

(4)提高微噴管壁面初始溫度，熱損失降低，粘性力增加，總體上比沖小幅增大。

(5)增大微噴管尺寸，粘性效應減弱，微噴管性能提高。

［1］Darren L Hittl，Charles M Zakrzwski，Michael A Thomas.MEMS－based satellite micropropulsion via catalyzed hydrogen peroxide decomposition［J］.Smart Materials and Structures，2001，10(2):1163－1175.

［2］Norman Chigier，Tevfik Gemci.A review of micro propulsion technology［J］.American Institute of Aeronautics and Astronautics，2003，670(2):1－12.

［3］Louisos W F，Hitt D L.Heat transfer＆ viscous effects in 2D＆ 3D supersonic micro－nozzle flows［C］//37th AIAA Fluid Dynamics Conference，Miami，2007.

［4］Louisos W F，Alexeenko A A，Hitt D L，et al.Design considerations for supersonic micronozzles［J］.International Journal of Manufacturing Research，2008，3(1):80－113.

［5］Louisos W F，Hitt D L.Viscous effects in supersonic micronozzle flows:transient analysis［C］//36th AIAA Fluid Dynamics Conference and Exhibit，AIAA 2006－2874.

［6］Alexeenko A A，Gimelshein S F.Numerical modeling of axisymmetric and three－dimensional flows in MEMS nozzles［C］//36th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit，AIAA 2000－3668.

［7］Bayt R L，Breuer K S.Analysis and testing of a silicon intrinsic－point heater in a micropropulsion application［J］.Sensors and Actuators A:Physical，2001，91(3):249－255.

［8］Bayt R L，Breuer K S.Fabrication and testing of micronsized cold－gas thrusters［J］.Progress in Astronautics and Aeronautics，2000，187(1):381－398.

［9］Kujawa J，Hitt D L，Cretu G.Numerical simulation of supersonic flow in realistic nozzle geometries with heat loss［C］//33rd AIAA Fluid Dynamics Conference，2003.

［10］李德桃，鄧軍，等.微型發(fā)動機燃燒室的模擬研究［J］.機械工程學報，2002，38(10):59－61.

［11］楊海威，朱衛(wèi)兵，等.微型推進器流場數(shù)值模擬［J］.導航與航天運載技術，2009，5(4):53－56.

［12］李宜敏，張中飲，等.固體火箭發(fā)動機原理［M］.北京:北京航空航天大學出版社，1991:69－74.

［13］李運澤，寧獻文，等.納衛(wèi)星軌道熱環(huán)境仿真分析［J］.系統(tǒng)仿真學報，2007，19(14):3166－3175.