散熱損失對微推進器性能的影響

胡松啟，劉　歡，武冠杰，陳　靜，高勝靈

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

?

散熱損失對微推進器性能的影響

胡松啟，劉歡，武冠杰，陳靜，高勝靈

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

微推進器尺寸小，散熱損失對微推進器性能有較大影響。本文研究微推進器陣列燃燒室壁散熱情況；分析熱損失對單元發(fā)動機性能的影響；試驗測量微型發(fā)動機推力，并與理論計算值對比分析。研究結(jié)果表明，對于微推進器陣列，燃燒室壁面處溫度最高，熱應(yīng)力和熱變形也最大；硅材料的室壁熱應(yīng)力以及熱變形均比鋼材料的小。微推進器噴管喉部散熱損失最為嚴重。微推進器散熱損失使噴管出口溫度降低，出口速度減小，推力減小11%。對不同單元發(fā)動機實測推力值分析，發(fā)現(xiàn)各單元實測推力有差異，但都小于絕熱理論計算推力；實測推力的平均值比散熱損失計算值小。

微推進器；散熱損失；微推力；溫度分布；熱應(yīng)力

0　引言

隨著微衛(wèi)星技術(shù)的發(fā)展和其應(yīng)用領(lǐng)域的擴大，對微推進系統(tǒng)的需求越來越迫切[1]。其中，基于MEMS的固體微推進器應(yīng)用較廣，它具有尺寸質(zhì)量小、高度集成、工作時間短、費用低等特點[2]。微推進器由于尺寸的微型化，其工作特點和常規(guī)大尺寸火箭發(fā)動機大相徑庭。微推進器的小尺寸效應(yīng)，使得熱損失等表面效應(yīng)成為影響微推進器性能的主要因素之一。Zhang K L和Chou S K[3]利用CFD建模研究發(fā)現(xiàn)，燃燒室壁面熱損失對發(fā)動機性能有較大影響，可使總沖損失5%～7%。Leach等[4]在研究壁面軸向熱傳導(dǎo)和環(huán)境熱損失對微管內(nèi)氫氣/空氣燃燒效率的影響中認為，隨著燃燒室尺寸減小，逐漸加強的壁面軸向熱傳導(dǎo)，有利于燃燒的穩(wěn)定和能量密度的提高。劉茂省、楊衛(wèi)娟等[5]模擬圓管內(nèi)氫氣/空氣的預(yù)混燃燒，發(fā)現(xiàn)微燃燒器散熱量占反應(yīng)放熱的10%。其中，圓管外壁面輻射散熱量最大，占整個熱損失的80%～90%。楊海威、朱衛(wèi)兵等[6]對微噴管內(nèi)流場進行了模擬仿真。結(jié)果表明，隨入口壓力升高，微噴管壁面的熱損失也隨之增加。微推進器殼體材料的傳熱特性對工作性能有較大影響，周海清等[7]對幾種不同傳熱特性的殼體材料進行數(shù)值計算。研究表明，常用的硅材料將造成較大的比沖損失，玻璃及陶瓷材料具有良好的絕熱性能，能夠有效改善微推進器熱損失。Alexander B Kiskin等[8]研究發(fā)現(xiàn)，喉徑為850～900 μm的微推進器的熱損失達到20%。Moríigo José等[9]用溫度云圖分析了微推進器壁面熱損失與燃氣熱通量的關(guān)系。Grujicie[10]研究了微電子封裝中界面熱阻對傳熱的影響。他提出在界面兩種材料之間加入熱膨脹較小、且熱導(dǎo)率較高的界面材料，通過試驗仿真證實了在界面間加入該種界面材料，可較大程度地減小熱阻。

熱損失會使微推進器總沖和推力減小，從而降低微推進器性能。因此，研究熱損失對微型發(fā)動機性能的影響十分必要。本文研究微推進器陣列燃燒室散熱情況，分析散熱損失對單元發(fā)動機性能的影響；最后，試驗測量微型發(fā)動機推力，并與理論計算值對比分析。

1　微推進器陣列燃燒室散熱分析

1.1仿真模型與條件

選取硅和鋼2種材料作為燃燒室壁面材料，分別計算溫度場，熱應(yīng)力分布和形變量分布，分析2種材料的散熱性能。用混合比為5∶5的斯蒂芬酸鉛和硝化棉作為固體推進劑。利用最小吉布斯自由能法計算燃氣的熱力學(xué)參數(shù)，燃氣溫度設(shè)為2 000 K。硅的密度ρ= 2 300 kg/m3，熱導(dǎo)率λ= 160 W/(m·K)，比定壓熱容cp= 700 J/(kg·K)，線膨脹系數(shù)a= 4.15×10-6K-1，泊松比μ= 0.27，彈性模量E= 131 GPa。鋼的密度ρ= 7 850 kg/m3，熱導(dǎo)率λ= 48 W/(m·K)，比定壓熱容cp= 480 J/(kg·K)，線膨脹系數(shù)a= 1.2×10-5K-1，泊松比μ= 0.28，彈性模量E= 210 GPa。

選用的固體微推進器燃燒室橫截面為圓形。燃燒室裝填推進劑燃燒柱后，將燃燒室和噴管鍵合起來，構(gòu)成完整的微型固體發(fā)動機。微推進器陣列中，燃燒室數(shù)量較多，兩燃燒室之間的間隔距離很近。研究對象取為燃燒室層，燃燒室直徑1 000 μm，兩單元中心間距1 600 μm。

材料的熱傳導(dǎo)方程為

式中λ為熱導(dǎo)率；ρ為材料密度；cp為比定壓熱容；τ為時間步長。

材料的應(yīng)力-應(yīng)變的本構(gòu)方程選用廣義胡克定律。由于燃燒室層結(jié)構(gòu)對稱，只對其兩個燃燒室的1/2進行計算分析，如圖1所示。

采用結(jié)構(gòu)化網(wǎng)格計算區(qū)域進行劃分。模型邊界條件為對稱邊界，加固定約束，1、2為燃燒室壁，燃氣與燃燒室壁面發(fā)生熱作用，在1、2上加熱載荷，邊界條件設(shè)定為第三類邊界條件，燃燒室層初溫設(shè)為300 K。

圖1　微推進器陣列燃燒室網(wǎng)格圖

1.2仿真結(jié)果分析

壁面初始溫度設(shè)為300 K，計算時間設(shè)為10 ms，圖2是10 ms時燃燒室溫度分布。此時，硅材料燃燒室沿徑向溫度梯度比鋼材料燃燒室壁小。這是由于鋼的熱導(dǎo)率較小，單位時間傳導(dǎo)的熱量較少。由圖2可知，2種材料的最高溫度均出現(xiàn)在燃燒室壁處，硅材料燃燒室壁的最高溫度為328.272 K。鋼材料燃燒室壁的最高溫度320.772 K。因此，分別取燃燒室壁處溫度變化作曲線，結(jié)果如圖3所示。由圖3發(fā)現(xiàn)，兩種材料壁面處溫度隨時間的增加而升高，硅材料燃燒室壁面的升溫速率比鋼快。

(a) 硅材料

(b) 鋼材料

圖3　燃燒室壁面處溫度隨時間變化曲線

圖4為10 ms時2種材料燃燒室壁等效熱應(yīng)力分布圖，圖5為10 ms時2種材料燃燒室壁熱變形量分布圖。由此組圖可知，最大等效熱應(yīng)力與最大熱變形均出現(xiàn)在燃燒室壁面處。

(a) 硅材料

(b) 鋼材料

(b) 鋼材料

圖6(a)為不同燃燒室壁材料下，燃燒室最大等效熱應(yīng)力隨時間變化曲線；圖6(b)為不同燃燒室壁材料下，燃燒室最大形變量隨時間變化曲線。由圖6可知，2種材料的最大等效熱應(yīng)力和最大形變量均隨時間增加而增大。硅材料微推進器室壁最大等效熱應(yīng)力及熱變形量分別為39 MPa和51.8 nm，鋼材料微推進器室壁最大等效熱應(yīng)力及熱變形量分別為101 MPa和67.5 nm，大于硅。這是由于硅的導(dǎo)熱性更好，更易建立熱平衡，彈性模量和線膨脹系數(shù)小于鋼，受熱作用時，產(chǎn)生較小的熱應(yīng)力與形變量。

由數(shù)值模擬結(jié)果可知，微推進器工作時，燃燒室壁處溫度最高，壁面處溫度隨時間的增加而升高，硅材料燃燒室壁的升溫速率比鋼快。硅材料燃燒室的最高溫度為328.272 K，鋼材料燃燒室的最高溫度為320.772 K，它們的最高溫度都遠小于各自的熔點。因此，推進器工作時，對陣列材料影響不大。同時，燃燒室壁處產(chǎn)生的等效熱應(yīng)力和熱形變量也最大；以硅為材料的微推進器室壁最大等效熱應(yīng)力及熱形變量分別為39 MPa和51.8 nm，以鋼為材料的微推進器室壁最大等效熱應(yīng)力及熱形變量分別為101 MPa和67.5 nm。由此可知，硅比鋼更適合作微推進器陣列材料。

(a)熱應(yīng)力

(b)形變量

2　散熱損失對單元發(fā)動機性能影響分析

2.1計算模型與邊界條件

數(shù)值計算所用固體微推進器的殼體材料為硅，分析熱損失對微推進器單元發(fā)動機工作性能的影響。本節(jié)計算的微推進器單元發(fā)動機結(jié)構(gòu)如圖7所示，微噴管的收斂和擴張半角均為35.3°，具體尺寸如圖7所示，單位為μm。

圖7　單元發(fā)動機結(jié)構(gòu)尺寸

為簡化計算，作如下假設(shè)：

(1)裝藥的化學(xué)反應(yīng)產(chǎn)物的熱力學(xué)性能是均一的、不變的；

(2)所有做功的工質(zhì)都是氣相的，忽略凝聚相的做功；

流體動力學(xué)控制方程選用非定常N-S方程。流體熱傳導(dǎo)控制方程為三維常物性、無內(nèi)熱源、非穩(wěn)態(tài)的導(dǎo)熱微分方程：

式中λ為熱導(dǎo)率；ρ為材料密度；cp為比定壓熱容；τ為時間步長。

湍流模型選用低Re數(shù)κ-ε模型，其控制方程為

式中ρ為流體密度；μ、μt為層流和湍流粘性系數(shù)；Cμ、C1、C2為模型常數(shù)；fμ、f1、f2為衰減函數(shù)；σk、σε為普朗特數(shù)；Ui、Uj為x、y方向上速度分量的時均分布值；k為湍動能；ε為湍動能耗散率；D、E分別為引入的附加項，表達式見文獻[11-12]。

燃燒室壓強設(shè)為1 MPa，溫度為3 000 K，噴管出口壓強設(shè)為10 Pa，溫度為300 K。燃氣的熱導(dǎo)率為0.710 35 W/(m·K)，密度為1.253 kg/m3，比定壓熱容為2 715 J/(kg·K)，粘性系數(shù)為1.008 6×10-4Pa·s，質(zhì)量流率設(shè)為0.1 g/s。固體壁面密度為2 400 kg/m3，比定壓熱容為700 J/(kg·K)，熱導(dǎo)率為157 W/(m·K)。

用結(jié)構(gòu)網(wǎng)格對計算區(qū)域進行劃分， 模型邊界條件主要有固體壁面、質(zhì)量流率入口邊界、壓強出口邊界條件，壁面均采用無滑移邊界，網(wǎng)格示意圖如圖8所示。

2.2仿真結(jié)果分析

首先，以文獻[13]中的物理模型為例進行推力計算。其中，燃燒室尺寸選取r= 0.5 cm，其他參數(shù)見原文。由文中數(shù)據(jù)可求得質(zhì)量流率為0.1、0.2 g/s所對應(yīng)推力分別為0.074、0.157 N。該物理模型使用本文計算模型時，質(zhì)量流率為0.1、0.2 g/s所求得推力分別為0.075、0.159 N。推力數(shù)據(jù)大小分別相差1.4%、1.3%，誤差相對較小。由此可見，本文計算模型可靠。

壁面初始溫度設(shè)為300 K，計算時間設(shè)為2 ms，圖9為2 ms時考慮壁面?zhèn)鳠釙r微推進器溫度云圖。圖10為假設(shè)壁面絕熱的理想情況下2 ms時微推進器溫度云圖。與考慮傳熱的情況相比，噴管出口溫度更高。

圖8　計算網(wǎng)格示意圖

圖9　2 ms時壁面散熱情況下微推進器溫度云圖

圖10　2 ms 時壁面絕熱情況下微推進器溫度云圖

圖11為2 ms時軸線上壁面絕熱和傳熱情況下溫度變化曲線。由圖11可看出，2條曲線變化趨勢大致相同，在傳熱情況下，靠近噴管出口一側(cè)，溫度下降梯度較大，噴管出口溫度較低，壁面的導(dǎo)熱效應(yīng)較為顯著。圖12為微推進器壁面上2點溫度隨時間變化規(guī)律，點1在燃燒室壁面中點附近，點2在噴管喉部附近。由圖12可看出，壁面上2點溫度都隨時間的增加而升高，點1的溫度升高速率較小，且?guī)缀醪蛔儯c2的溫度先快速升高，后緩慢升高。點2溫度一直高于點1溫度。

圖11　2 ms時絕熱與散熱時軸線溫度對比曲線

圖12　微推進器壁面上2點溫度隨時間變化

圖13為2 ms時壁面?zhèn)鳠岷捅诿娼^熱情況下，微推進器單元發(fā)動機速度云圖。圖14為2 ms時2種情況下軸線上的速度變化曲線。由圖14可知，兩者的速度變化趨勢基本一致，燃燒室內(nèi)速度較小，進入噴管后，速度快速增大；壁面?zhèn)鳠崤c壁面絕熱相比，噴管出口軸向速度較小。

根據(jù)發(fā)動機推力公式，求出2種情況下單元發(fā)動機的推力值并作圖，如圖15所示。由圖15可知，在壁面絕熱的理想情況下，推進器能夠較快地到達穩(wěn)定平衡工作段，工作段推力為0.206 2 N。考慮散熱損失時，工作段推力約為0.183 6 N。由圖15中數(shù)據(jù)計算發(fā)現(xiàn)，考慮壁面熱損失，其推力要比絕熱壁面推力減小11.0%。

由數(shù)值模擬結(jié)果可知，熱損失使噴管出口溫度降低，出口軸向速度減小。殼體硅材料的散熱損失，使微推進器推力減小了11%，其熱損失已達到不可忽視的程度。若采用硅作為推進器殼體材料時，設(shè)計過程中，必須考慮熱損失，采取適當方法進行彌補。

(a) 壁面?zhèn)鳠?/p>

(b) 壁面絕熱

圖14　2種情況下軸線上的速度曲線

圖15　2種不同工況的推力曲線

3　單元發(fā)動機推力測試試驗

利用試驗測量單元發(fā)動機推力，對比分析實測推力與理論計算值之間的差異。

3.1試驗方案

采用懸擺法對固體微推進器進行沖量測試。沖量是力對時間的積累，其大小等于受力載體的動量增量。懸擺是一種將待測沖量轉(zhuǎn)化為擺動角度的間接沖量測量裝置。該系統(tǒng)主要由懸擺、激光干涉儀、阻尼器和減震平臺等設(shè)備構(gòu)成。

3.2推力實測結(jié)果

實際測試的單元發(fā)動機與數(shù)值模擬計算所用發(fā)動機型號相同。試驗選擇LS/NC為5∶5的推進劑配方，對多個單元發(fā)動機進行了點火測試，所得推力大小見圖16。各單元發(fā)動機實測推力變化范圍為0.147 6～ 0.204 3 N，單元發(fā)動機考慮散熱時理論推力為0.183 6 N，絕熱時理論推力為0.206 2 N?？梢姡?/p>

(1)實測推力數(shù)據(jù)呈散點分布，這是因為實測時推力大小與裝藥量、裝藥密度、點火絲與推進劑接觸緊密程度等因素有密切關(guān)系，導(dǎo)致不同單元發(fā)動機實測推力有差異。加熱電阻絲的阻值不同，點火功率不同，導(dǎo)致推進劑點火燃燒的劇烈程度不同，影響了推進劑的燃燒時間和能量，也對發(fā)動機實測推力不同有影響。

(2)實測推力都小于絕熱推力。發(fā)動機實際工作時，存在散熱損失等影響推力的因素。因此，絕熱理論計算推力比實測推力大。

(3)實測推力與散熱損失推力相比，有高有低，但實測推力的平均值比散熱損失計算值小。發(fā)動機工作時，除散熱損失外，還存在粘性損失，燃燒不完全損失，摩擦損失等，這是實測推力的平均值比散熱損失計算值小的原因。

通過理論計算推力與實測推力的比較分析，可證明理論計算所用模型是正確的，實測值與理論計算值一致。對于單元發(fā)動機推力存在差異，可進行如下改進：首先，應(yīng)規(guī)范燃燒方法，統(tǒng)一燃燒質(zhì)量與燃燒密度；其次，針對電阻阻值不同，應(yīng)改善加工工藝，以保證電阻阻值相同，加熱功率一致。

圖16　單元發(fā)動機推力實測數(shù)據(jù)

4　結(jié)論

(1)微推進器工作時，燃燒室壁受熱升溫，但燃燒室壁最高溫度都遠小于殼體材料的熔點，因而燃燒室壁散熱損失對推進器陣列影響不大。推進器工作10 ms時，以硅為材料的微推進器室壁最大等效熱應(yīng)力以及熱變形分別為39 MPa和51.8 nm，鋼為材料的微推進器室壁最大等效熱應(yīng)力以及熱變形分別為101 MPa和67.5 nm。

(2)對于微推進器，噴管喉部散熱損失最為嚴重。

(3)根據(jù)數(shù)值模擬結(jié)果可知，微推進器散熱損失，使噴管出口溫度降低，出口軸向速度減??；殼體和噴管的散熱損失，使推進器推力減小了11%。

(4)不同單元發(fā)動機實測推力范圍為0.147 6 ～0.204 3 N，差異較大，這是因為不同單元發(fā)動機的裝藥量、裝藥密度、加熱電阻的阻值、點火絲與推進劑接觸緊密程度不同，使得單元發(fā)動機內(nèi)推進劑燃燒的劇烈程度不同，影響了推進劑的燃燒時間和能量。

(5)實測推力都小于絕熱推力，實測推力與散熱損失推力相比，有高有低，但實測推力的平均值比散熱損失計算值小。

[1]張高飛, 尤政, 胡松啟, 等. 基于MEMS的固體推進器陣列[J]. 清華大學(xué)學(xué)報(自然科學(xué)版), 2004, 44(11): 1489-1492.

[2]洪延姬, 金星, 崔村燕, 等. 先進航天推進技術(shù)[M]. 北京：國防工業(yè)出版社, 2012: 140-142.

[3]Zhang K L, Chou S K. Performance prediction of a novel solid-propellant micro thruster[J]. Journal of Propulsion and Power, 2006, 22(1): 56-63.

[4]Timothy T Leach, Chistopher P Cadou. The role of structural heat exchange and heat loss in the design of efficient silicon micro-combustors[J]. Proceedings of the Combustion Institute, 2005, 30(2): 2437-2444.

[5]劉茂省, 楊衛(wèi)娟, 周俊虎, 等. HZ/空氣預(yù)混燃燒模擬計算和散熱的研究[J]. 燃燒科學(xué)與技術(shù), 2008, 14(3): 259-264.

[6]楊海威, 朱衛(wèi)兵, 趙陽. 微型推進器流場數(shù)值模擬[J]. 導(dǎo)彈與航天運載技術(shù), 2009(5): 53-56.

[7]周海清, 張高飛, 尤政. 固體微推力器設(shè)計與數(shù)值分析[J]. 推進技術(shù), 2007, 28(3): 230-234.

[8]Alexander B Kiskin, Vladimir N Simonenko, Lev K Gusachenko, et al. Experimental study of microthruster heat loss[J]. International Journal of Energetic Materials and Chemical Propulsion, 2009, 8(4): 329-344.

[11]Launder B E, Sharma B I. Application of the energy-dissi-pation model of turbulence to the calculation of flow neara spinning disc[J]. Letters in Heat and Mass Transfer, 1974, 1(74): 131-138.

[12]Chang K C, Hsieh W D, Chen C S. A modified low-Reynolds-number turbulence model applicable to recirculating flow in pipe expansion[J]. Transactions of the ASME，Journal of Fluids Engineering, 1995, 117(3): 417-423.

[13]林博穎, 張根烜, 劉明候, 等. 微小型化學(xué)推進器的性能分析[J]. 中國科學(xué)技術(shù)大學(xué)學(xué)報, 2008, 38(1): 94-104.

(編輯：崔賢彬)

Influences of heat loss on micro thruster

HU Song-qi, LIU Huan, WU Guan-jie, CHEN Jing, GAO Sheng-ling

(Science and Technology on Combustion, Internal Flow and Thermo-Structure Laboratory, Northwestern Polytechnical University, Xi'an710072, China)

The performance of micro thruster is influenced greatly by heat loss due to the small size of micro thruster. The heat loss of the combustion chamber wall in micro thruster array and the performance of unit engine affected by the heat loss were studied. Thrust produced by micro thruster was measured and compared with the theoretical value. The results indicate that temperature, thermal stress and thermal deformation are at their maximum on the surface of combustion chamber wall in micro thruster array. The thermal stress and thermal deformation of silicon wall are smaller than the steel ones. The heat loss at nozzle throat is the most serious. The heat loss makes temperature of nozzle outlet and exit velocity decrease, and the thrust is reduced by 11%. For different engine unit, the measured thrust is diverse, but it is smaller than the adiabatic thrust. What's more, mean of the measured thrust is smaller than the value obtained from the heat loss calculation.

micro thruster；heat loss；micro thrust；temperature distribution；thermal stress

2014-07-08；

2015-10-28。

西北工業(yè)大學(xué)基礎(chǔ)研究基金(JC201221)。

胡松啟(1976—)，教授，研究領(lǐng)域為固體推進劑及燃燒。E-mail:pinecore@nwpu.edu.cn

V435

A

1006-2793(2016)01-0050-06

10.7673/j.issn.1006-2793.2016.01.009