微型固體姿控發動機微噴管內氣粒兩相流動規律的CFD-DSMC研究①

夏廣慶，張 斌，孫得川，陳茂林，3

(1.大連理工大學工業裝備結構分析國家重點實驗室，大連 116024;2.中國航天科工集團公司第九總體設計部，武漢 430040;3.西北工業大學燃燒、熱結構與內流場重點實驗室，西安 710072)

微型固體姿控發動機微噴管內氣粒兩相流動規律的CFD-DSMC研究①

夏廣慶1，張 斌2，孫得川1，陳茂林1，3

(1.大連理工大學工業裝備結構分析國家重點實驗室，大連 116024;2.中國航天科工集團公司第九總體設計部，武漢 430040;3.西北工業大學燃燒、熱結構與內流場重點實驗室，西安 710072)

微型固體姿控發動機在航天領域具有廣泛的應用前景。以基于MEMS技術的微噴管為研究對象，首先通過計算微噴管中的克努森數，得到了微噴管中的氣相流動狀態;然后，采用CFD-DSMC方法，模擬了微噴管中的氣粒兩相流動，并研究了顆粒相質量分數和粒徑對氣相流動的影響。結果表明，在所研究的來流條件下，微噴管中的連續介質假設是成立的;氣相與顆粒相間的動量和能量交換，導致氣相馬赫數降低、溫度升高，同時也導致顆粒相速度增加、溫度降低;顆粒相質量分數和粒徑均能顯著影響氣相的馬赫數和溫度。

微型固體姿控發動機;微噴管;CFD;DSMC

0 引言

固體姿控發動機結構簡單、可靠性高，在有效提高飛行器機動性能的同時，還能縮短控制系統的響應時間。因此，在國內外已得到了廣泛的研究和使用［1-2］，其典型代表如使用180個固體姿控發動機作為末端控制裝置的美國增程攔截彈ERINT-1。微型固體姿控發動機除具備以上優勢外，以其在小型航天器上所具有的廣泛應用前景，更成為當前的一個研究熱點，其代表如基于MEMS技術的固體微推力器及其陣列［3-6］。

微型固體姿控發動機微噴管中的氣粒兩相流動一方面引發了比沖損失，同時又加劇了微噴管壁面的燒蝕。另一方面，對微噴管內氣粒兩相流動規律的研究是開展微型固體姿控發動機氣粒兩相羽流污染、氣動力和氣動熱效應研究的首要前提［7-9］。由此可見，開展微型固體姿控發動機微噴管內氣粒兩相流動規律的研究工作，不僅對其優化設計具有十分重要的理論價值，同時對微型航天器的結構優化也具有十分重要的現實意義。

本文中，以基于MEMS技術加工的微噴管為研究對象。首先，開展了微噴管中純氣相流動規律研究，在得到氣相主要特性參數分布規律的同時，通過計算流場中的克努森數，驗證了氣相計算方法的有效性;然后，采用CFD-DSMC方法，開展了微噴管中的氣粒兩相流動規律研究，得到了氣相和顆粒相主要特征參數的分布規律，并研究了顆粒相濃度和粒徑對氣相相關參數的影響。

1 數學模型

氣相控制方程組為不計體積力和無內熱源的三維非定常N-S方程組。在笛卡爾坐標系下，其積分形式由式(1)表示為

顆粒相的控制方程組由式(4)～式(6)表示。

式中Vcell為網格單元體積;下標k表示第k束顆粒;Nk為第k束顆粒的真實顆粒數;Dk為第k束顆粒的直徑;為第k束顆粒的滑移速度;CDk為第k束顆粒的滯止系數;Cp和C分別為顆粒和氣體的比熱容;Nuk為第k束顆粒的努謝爾特數;Pr為普朗特數。

2 數值方法

氣相控制方程組采用中心有限體積法求解。其中，對流項采用添加各向異性人工粘性的中心差分格式求解，擴散項采用中心差分格式求解，時間項采用簡化五步Runge-Kutta法求解，并采用隱式殘值光順法加速收斂。

對顆粒相方程組的求解，也采用時間推進法。在每一個時刻步開始時，類似經典的DSMC方法，在入口邊界處加入新的顆粒，其速度和溫度與氣相一致，位置隨機給定。在該時刻步內，計算顆粒的位置、速度和溫度，并計算顆粒相對氣相作用的動量和能量源項，同時對超出計算區域的顆粒施加相應的邊界條件。在一個時刻步結束后，對顆粒的位置進行重新編號。在流場穩定后，統計流場信息。

對氣相和顆粒相之間的耦合計算，采用PSIC算法。計算對象為方形截面的硅質微噴管，采用濕法化學腐蝕技術加工(腐蝕夾角約為54.7°)，其結構示意圖如圖1所示(單位μm)。計算中，采用結構網格離散計算區域，其示意圖如圖2所示。

圖1 微噴管幾何結構示意圖Fig.1 Structure diagram of micro nozzle

圖2 計算網格示意圖Fig.2 Grid diagram for computation

氣相入口給定來流的總溫、總壓、馬赫數和方向角，氣相出口處各流動參數均由內場按二價外插獲得，氣相壁面采用絕熱、無滑移固壁邊界條件。

顆粒相入口給定顆粒的質量分數、溫度、速度和相關物性參數(密度為 4 004.6 kg/m3，比熱容為1 380 J/(kg·K))，顆粒相出口處施加吸收邊界，顆粒相壁面邊界按Tabakoff［11］的經驗公式，確定顆粒與壁面碰撞后的速度大小和方向。

3 計算結果與分析

3.1 微噴管中純氣相流動規律研究

在來流總溫2 982 K、總壓2 MPa條件下，微噴管中氣相主要特性參數的分布規律如圖3所示。由圖3可知，由于微噴管加工工藝的限制，微噴管擴張段中壁面附近氣流的膨脹并不充分，尤其是邊角處氣相的能量轉換效率最低，導致氣相的馬赫數較低、溫度較高。在微噴管的出口截面上，氣相亞聲速區域面積與該截面總面積的比值在擴張段內所有截面中位居最高，約為微噴管出口截面總面積的1/6。同時，氣相在微噴管出口截面邊角處的溫度也在該截面中最高，其值約為2 600 K。圖4為來流總溫2 982 K、總壓分別為1、2、3 MPa 3種情況下微噴管中的克努森數分布情況。由圖4可知，來流總壓越大，噴管中的克努森數越小。克努森數較大的區域位于噴管出口處的壁面附近。當來流總壓1 MPa時，其最大值約0.008 5，微噴管中的流動為符合連續介質假設的無滑移流。可見，以上3種工況下，在全場采用無滑移邊界的N-S方程進行流場求解是完全合理的。

圖3 微噴管中氣相主要特性參數分布規律Fig.3 Distribution contours of gas-phase characteristic parameters in the micro nozzle

圖4 不同來流總壓下微噴管中流動狀態Fig.4 Flow state in the micro nozzle under different total pressures

3.2 微噴管中氣粒兩相流動規律研究

3.2.1 氣相和顆粒相主要特性參數的分布規律

當顆粒質量分數為9.4%、粒徑為1 μm時，氣相主要特性參數的分布規律如圖5所示(氣相來流條件同3.1節)。由圖5可知，顆粒相的存在對氣相壓強場和密度場的影響不大，但對氣相溫度場和馬赫數場的影響較大，尤其是在擴張段中，在顆粒濃度較大的位置，氣相場的溫度升高，馬赫數降低。與3.1節中純氣相流動計算結果對比可知，在微噴管出口截面的中心線處，顆粒相的存在導致氣相在該處的馬赫數降低約6%，溫度升高約8%。

圖5 氣相主要特性參數分布規律Fig.5 Distribution contours of gas-phase characteristic parameters

顆粒相主要特性參數的分布規律如圖6所示。由圖6可知，沿流動方向顆粒相的速度逐漸增大，溫度逐漸降低。在噴管的擴張段中，大部分區域為無粒子區(圖6中參數無梯度變化的區域為無粒子區)。

圖6 顆粒相主要特性參數分布規律Fig.6 Distribution contours of particle-phase characteristic parameters

3.2.2 顆粒質量分數對氣相馬赫數和溫度的影響

在不同顆粒質量分數下，氣相流場中心線處的馬赫數和溫度變化分別如圖7和圖8所示。由圖7和圖8可知，隨顆粒質量分數的增加，顆粒對氣相運動的阻礙效應增強，導致氣相的馬赫數呈下降趨勢。同時，由于顆粒質量分數的增加，也直接導致了顆粒相與氣相間換熱效應的增強，從而導致氣相溫度呈上升趨勢。

3.2.3 顆粒粒徑對氣相馬赫數和溫度的影響

在不同顆粒粒徑下，氣相流場中心線處的馬赫數和溫度變化分別如圖9和圖10所示。由圖9和圖10可知，隨顆粒粒徑的增大，顆粒相與氣相接觸的表面積減小，對氣相的阻力下降，從而導致中心線上氣相的馬赫數和速度均呈上升趨勢。同時，相間接觸面積的減小，導致了顆粒相與氣相間換熱效應的減弱，從而導致氣相溫度呈下降趨勢。

圖7 不同顆粒質量分數下微噴管中心線上氣相馬赫數變化Fig.7 Mach number distribution of gas-phase along the micro nozzle centerline under different particle mass fractions

圖8 不同顆粒質量分數下微噴管中心線上氣相溫度變化Fig.8 Temperature distribution of gas-phase along the micro nozzle centerline under different particle mass fractions

圖9 不同顆粒粒徑下微噴管中心線上氣相馬赫數變化Fig.9 Mach number distribution of gas-phase along the micro nozzle centerline under different particle diameters

圖10 不同顆粒粒徑下微噴管中心線上氣相溫度變化Fig.10 Temperature distribution of gas-phase along the micro nozzle centerline under different particle diameters

4 結論

(1)在氣相來流總溫2 982 K、總壓1 MPa的條件下，微噴管中氣相克努森數的最大值約為0.008 5，氣相流動為符合連續介質假設的無滑移流。

(2)顆粒相的存在對氣相溫度場和馬赫數場影響較大，尤其是在氣相特性參數大幅變化的擴張段中，顆粒相的存在導致氣相馬赫數降低、溫度升高。在氣相來流總溫2 982 K、總壓2 MPa、顆粒質量分數9.4%、粒徑1 μm的條件下，顆粒相的存在導致微噴管出口截面中心線處的氣相馬赫數降低約6%，溫度升高約8%。

(3)隨來流顆粒質量分數的增加，氣相流場中心線處的馬赫數和速度呈下降趨勢，溫度呈上升趨勢。隨著來流顆粒粒徑的增大，氣相流場中心線處的馬赫數和速度均呈上升趨勢，而溫度呈下降趨勢。

［1］邢繼發，劉國球，黃堅定，等.世界導彈與航天發動機大全［M］.北京:軍事科學出版社，1999.

［2］張鋼錘，葉定友.一種新概念固體姿控發動機簡介［J］.固體火箭技術，2001，24(3).

［3］Muller J，Marrese C，Polk J，et al.An overview of MEMS-based micropropulsion developments at JPL［J］.Acta Astronautica，2003，52(9).

［4］Rossi C，Orieux S，Larangot B，et al.Design fabrication and modeling of solid propellant microrocket-application to micropropulsion［J］.Sensors and Actuators，2002，99(1).

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

［6］劉旭輝，方蜀州.微型固體推力器陣列尋址點火控制系統研究［J］.固體火箭技術，2010，33(6).

［7］張建華，賀碧蛟，蔡國飆，等.衛星姿控發動機噴管羽流撞擊效應試驗［J］.空氣動力學報，2007，25(2).

［8］Park J H，Baek S W，Kim J S.DSMC analysis of the interaction between thruster plume and satellite components［R］.AIAA 2002-0794.

［9］程曉麗，王強，閻喜勤.衛星姿控發動機高空羽流場工程分析［J］.空間科學學報，2006，26(1).

［10］陳偉芳，常雨.三維管道超聲速氣固兩相流動的CFD/DSMC 仿真［J］.推進技術，2005，26(3).

［11］Tabakoff W，Sugiyama Y.Experimental method of determining particle restitution coefficients［C］//New York:Proc.Symp.on Polyphase Flow and Transport，ASME，1980，8:2032210.

Research on the gas-particle two-phase flow in the micro nozzle of attitude control micro solid rocket motor

XIA Guang-qing1，ZHANG Bin2，SUN De-chuan1，CHEN Mao-lin1，3
(1.State Key Laboratory of Structural Analysis for Industrial Equipment，Dalian University of Technology，Dalian 116024，China;2.The Ninth System Design Department of China Aerospace Science and Industry Corporation，Wuhan 430040，China;3.Science and Technology on Combustion，Internal Flow and Thermal-Structure Laboratory，Northwestern Polytechnical University，Xi＇an 710072，China)

Attitude control micro solid rocket motor has wide application potential in the aerospace field.The gas-particle twophase flow in the micro nozzle based on the MEMS technology was investigated.Firstly，through calculating the Knudsen number of the micro nozzle，the gas phase flow state in the micro nozzle was obtained.Then the gas-particle two-phase flow in the micro nozzle was simulated by using the method of CFD-DSMC.The influence of particle mass fraction and particle diameter on the gas phase flow was studied.The result shows that the continuum assumption in the micro nozzle is established under the conditions of the defined flow in the study.The exchange of momentum and energy between the gas phase and the particle can not noly reduce the gas phase Mach number and raise the temperature，but also increase the particle phase velocity and decrease the temperature.The particle phase mass fraction and particle diameter can significantly influence the Mach number and temperature of gas phase.

attitude control micro solid rocket motor;micro nozzle;CFD;DSMC

V435

A

1006-2793(2012)03-0356-05

2012-01-09;

2012-04-08。

國家自然科學基金(11105023);新世紀優秀人才支持計劃(NCET-11-0054);工業裝備結構分析國家重點實驗室開放基金(GZ1101);中央高校基本科研業務費專項資金(DUT11ZD(G)02)。

夏廣慶(1979—)，男，講師，研究方向為發動機總體設計。E-mail:gq.xia@dlut.edu.cn

(編輯:崔賢彬)