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

(1.大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024;2.中國(guó)航天科工集團(tuán)公司第九總體設(shè)計(jì)部，武漢 430040;3.西北工業(yè)大學(xué)燃燒、熱結(jié)構(gòu)與內(nèi)流場(chǎng)重點(diǎn)實(shí)驗(yàn)室，西安 710072)

微型固體姿控發(fā)動(dòng)機(jī)微噴管內(nèi)氣粒兩相流動(dòng)規(guī)律的CFD-DSMC研究①

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

(1.大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024;2.中國(guó)航天科工集團(tuán)公司第九總體設(shè)計(jì)部，武漢 430040;3.西北工業(yè)大學(xué)燃燒、熱結(jié)構(gòu)與內(nèi)流場(chǎng)重點(diǎn)實(shí)驗(yàn)室，西安 710072)

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

微型固體姿控發(fā)動(dòng)機(jī);微噴管;CFD;DSMC

0 引言

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

微型固體姿控發(fā)動(dòng)機(jī)微噴管中的氣粒兩相流動(dòng)一方面引發(fā)了比沖損失，同時(shí)又加劇了微噴管壁面的燒蝕。另一方面，對(duì)微噴管內(nèi)氣粒兩相流動(dòng)規(guī)律的研究是開展微型固體姿控發(fā)動(dòng)機(jī)氣粒兩相羽流污染、氣動(dòng)力和氣動(dòng)熱效應(yīng)研究的首要前提［7-9］。由此可見，開展微型固體姿控發(fā)動(dòng)機(jī)微噴管內(nèi)氣粒兩相流動(dòng)規(guī)律的研究工作，不僅對(duì)其優(yōu)化設(shè)計(jì)具有十分重要的理論價(jià)值，同時(shí)對(duì)微型航天器的結(jié)構(gòu)優(yōu)化也具有十分重要的現(xiàn)實(shí)意義。

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

1 數(shù)學(xué)模型

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

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

式中Vcell為網(wǎng)格單元體積;下標(biāo)k表示第k束顆粒;Nk為第k束顆粒的真實(shí)顆粒數(shù);Dk為第k束顆粒的直徑;為第k束顆粒的滑移速度;CDk為第k束顆粒的滯止系數(shù);Cp和C分別為顆粒和氣體的比熱容;Nuk為第k束顆粒的努謝爾特?cái)?shù);Pr為普朗特?cái)?shù)。

2 數(shù)值方法

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

對(duì)顆粒相方程組的求解，也采用時(shí)間推進(jìn)法。在每一個(gè)時(shí)刻步開始時(shí)，類似經(jīng)典的DSMC方法，在入口邊界處加入新的顆粒，其速度和溫度與氣相一致，位置隨機(jī)給定。在該時(shí)刻步內(nèi)，計(jì)算顆粒的位置、速度和溫度，并計(jì)算顆粒相對(duì)氣相作用的動(dòng)量和能量源項(xiàng)，同時(shí)對(duì)超出計(jì)算區(qū)域的顆粒施加相應(yīng)的邊界條件。在一個(gè)時(shí)刻步結(jié)束后，對(duì)顆粒的位置進(jìn)行重新編號(hào)。在流場(chǎng)穩(wěn)定后，統(tǒng)計(jì)流場(chǎng)信息。

對(duì)氣相和顆粒相之間的耦合計(jì)算，采用PSIC算法。計(jì)算對(duì)象為方形截面的硅質(zhì)微噴管，采用濕法化學(xué)腐蝕技術(shù)加工(腐蝕夾角約為54.7°)，其結(jié)構(gòu)示意圖如圖1所示(單位μm)。計(jì)算中，采用結(jié)構(gòu)網(wǎng)格離散計(jì)算區(qū)域，其示意圖如圖2所示。

圖1 微噴管幾何結(jié)構(gòu)示意圖Fig.1 Structure diagram of micro nozzle

圖2 計(jì)算網(wǎng)格示意圖Fig.2 Grid diagram for computation

氣相入口給定來流的總溫、總壓、馬赫數(shù)和方向角，氣相出口處各流動(dòng)參數(shù)均由內(nèi)場(chǎng)按二價(jià)外插獲得，氣相壁面采用絕熱、無滑移固壁邊界條件。

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

3 計(jì)算結(jié)果與分析

3.1 微噴管中純氣相流動(dòng)規(guī)律研究

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

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

圖4 不同來流總壓下微噴管中流動(dòng)狀態(tài)Fig.4 Flow state in the micro nozzle under different total pressures

3.2 微噴管中氣粒兩相流動(dòng)規(guī)律研究

3.2.1 氣相和顆粒相主要特性參數(shù)的分布規(guī)律

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

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

顆粒相主要特性參數(shù)的分布規(guī)律如圖6所示。由圖6可知，沿流動(dòng)方向顆粒相的速度逐漸增大，溫度逐漸降低。在噴管的擴(kuò)張段中，大部分區(qū)域?yàn)闊o粒子區(qū)(圖6中參數(shù)無梯度變化的區(qū)域?yàn)闊o粒子區(qū))。

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

3.2.2 顆粒質(zhì)量分?jǐn)?shù)對(duì)氣相馬赫數(shù)和溫度的影響

在不同顆粒質(zhì)量分?jǐn)?shù)下，氣相流場(chǎng)中心線處的馬赫數(shù)和溫度變化分別如圖7和圖8所示。由圖7和圖8可知，隨顆粒質(zhì)量分?jǐn)?shù)的增加，顆粒對(duì)氣相運(yùn)動(dòng)的阻礙效應(yīng)增強(qiáng)，導(dǎo)致氣相的馬赫數(shù)呈下降趨勢(shì)。同時(shí)，由于顆粒質(zhì)量分?jǐn)?shù)的增加，也直接導(dǎo)致了顆粒相與氣相間換熱效應(yīng)的增強(qiáng)，從而導(dǎo)致氣相溫度呈上升趨勢(shì)。

3.2.3 顆粒粒徑對(duì)氣相馬赫數(shù)和溫度的影響

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

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

圖8 不同顆粒質(zhì)量分?jǐn)?shù)下微噴管中心線上氣相溫度變化Fig.8 Temperature distribution of gas-phase along the micro nozzle centerline under different particle mass fractions

圖9 不同顆粒粒徑下微噴管中心線上氣相馬赫數(shù)變化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 結(jié)論

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

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

(3)隨來流顆粒質(zhì)量分?jǐn)?shù)的增加，氣相流場(chǎng)中心線處的馬赫數(shù)和速度呈下降趨勢(shì)，溫度呈上升趨勢(shì)。隨著來流顆粒粒徑的增大，氣相流場(chǎng)中心線處的馬赫數(shù)和速度均呈上升趨勢(shì)，而溫度呈下降趨勢(shì)。

［1］邢繼發(fā)，劉國(guó)球，黃堅(jiān)定，等.世界導(dǎo)彈與航天發(fā)動(dòng)機(jī)大全［M］.北京:軍事科學(xué)出版社，1999.

［2］張鋼錘，葉定友.一種新概念固體姿控發(fā)動(dòng)機(jī)簡(jiǎn)介［J］.固體火箭技術(shù)，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］周海清，張高飛，尤政.固體微推力器設(shè)計(jì)與數(shù)值分析［J］.推進(jìn)技術(shù)，2007，28(3).

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

［7］張建華，賀碧蛟，蔡國(guó)飆，等.衛(wèi)星姿控發(fā)動(dòng)機(jī)噴管羽流撞擊效應(yīng)試驗(yàn)［J］.空氣動(dòng)力學(xué)報(bào)，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］程曉麗，王強(qiáng)，閻喜勤.衛(wèi)星姿控發(fā)動(dòng)機(jī)高空羽流場(chǎng)工程分析［J］.空間科學(xué)學(xué)報(bào)，2006，26(1).

［10］陳偉芳，常雨.三維管道超聲速氣固兩相流動(dòng)的CFD/DSMC 仿真［J］.推進(jìn)技術(shù)，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。

國(guó)家自然科學(xué)基金(11105023);新世紀(jì)優(yōu)秀人才支持計(jì)劃(NCET-11-0054);工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室開放基金(GZ1101);中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金(DUT11ZD(G)02)。

夏廣慶(1979—)，男，講師，研究方向?yàn)榘l(fā)動(dòng)機(jī)總體設(shè)計(jì)。E-mail:gq.xia@dlut.edu.cn

(編輯:崔賢彬)