王 倩 何小英 蔡國飆

（北京航空航天大學(xué)宇航學(xué)院）

1 引言

航天器的姿態(tài)控制往往利用姿態(tài)控制發(fā)動(dòng)機(jī)來進(jìn)行調(diào)節(jié)，發(fā)動(dòng)機(jī)噴出的羽流會(huì)對流場中的航天器帶來羽流污染、氣動(dòng)力、氣動(dòng)熱效應(yīng)等影響。羽流污染可能會(huì)導(dǎo)致航天器上敏感原件無法正常工作，甚至影響航天器的壽命，所以研究羽流的撞擊效應(yīng)就顯得非常重要。

國外對羽流撞擊效應(yīng)的研究開展較早，1986年德國的Dettleff試驗(yàn)研究了模型發(fā)動(dòng)機(jī)羽流場[1]；1984年法國的Allegre，在SR3風(fēng)洞中測量了羽流撞擊平板的氣動(dòng)力效應(yīng)[2]；1990年Legge測量了羽流撞擊傾斜擋板的氣動(dòng)力[3]；Boyd等人在1992年對低密度噴管流動(dòng)和羽流場進(jìn)行了試驗(yàn)研究和數(shù)值模擬[4]。在實(shí)驗(yàn)研究的基礎(chǔ)上，國外還進(jìn)行了相關(guān)數(shù)值模擬研究，1985年Legge等人對衛(wèi)星模型發(fā)動(dòng)機(jī)的羽流場和羽流撞擊效應(yīng)進(jìn)行了數(shù)值模擬[5]；2002年P(guān)ark數(shù)值模擬了羽流和衛(wèi)星的相互作用[6]。國內(nèi)由于試驗(yàn)條件的限制，研究工作主要集中在數(shù)值計(jì)算方面。

本文采用自開發(fā)的PWS軟件，對推力0.5N的模型發(fā)動(dòng)機(jī)進(jìn)行了羽流撞擊效應(yīng)DSMC模擬。工質(zhì)為CO2。首先，對流場的壓強(qiáng)、溫度、馬赫數(shù)、流線圖等進(jìn)行研究；然后，對擋板壓強(qiáng)及羽流污染情況進(jìn)行了分析。

2 DSMC方法及PWS軟件介紹

DSMC方法的基本思想是：利用少量的模擬分子代替真實(shí)流場內(nèi)數(shù)目眾多的氣體分子，用計(jì)算模擬實(shí)際物理過程，即對模擬分子進(jìn)行跟蹤，記錄它們的位置和速度的變化，計(jì)算分子之間，分子和物面間的碰撞，再經(jīng)統(tǒng)計(jì)平均，獲得所需的各種流動(dòng)參數(shù)。DSMC的基本假設(shè)是：模擬每個(gè)網(wǎng)格中的分子間的碰撞時(shí)，忽略分子的具體位置。

基于DSMC的PWS軟件的主程序包括網(wǎng)絡(luò)模塊、粒子參數(shù)模塊、粒子進(jìn)入模塊、粒子運(yùn)動(dòng)及邊界模塊、粒子碰撞模塊、并行模塊、統(tǒng)計(jì)輸出模塊共七大模塊。在這七大基本模塊的基礎(chǔ)上，針對各種實(shí)際問題，PWS羽流計(jì)算軟件都可以方便的搭建數(shù)值模擬應(yīng)用平臺(tái)[7]。

PWS軟件是根據(jù)模塊化和面向?qū)ο蟮脑瓌t對DSMC計(jì)算程序進(jìn)行改進(jìn)，具體如下：

（1）使用松散的軟件架構(gòu)保證各個(gè)計(jì)算模塊的相對獨(dú)立性，即改進(jìn)其中一個(gè)模塊時(shí)，基本上不需要改動(dòng)其他模塊；

（2）對各個(gè)模塊的功能實(shí)現(xiàn)采用通用性設(shè)計(jì)，即在使用單/多組分分子、不同分子模型、不同固體邊界條件時(shí)不需要改動(dòng)源代碼。

使用圖1的步驟生成所需邊界條件，這樣可以將粒子-固體邊界處理歸結(jié)到幾種基本類型，可以方便的實(shí)現(xiàn)通用性設(shè)計(jì)。

（1）定義幾種常用的平面和曲面（直線）構(gòu)成基本元素；

（2）使用多個(gè)平面或曲面相交構(gòu)成凸多面體；

（3）使用多個(gè)凸多面體來構(gòu)成最終的復(fù)雜邊界。

圖1 邊界生成過程說明

3 計(jì)算模型

本文采用推力0.5N模型發(fā)動(dòng)機(jī)，邊長300mm方形擋板，工質(zhì)為CO2，總壓Pc=8000Pa，總溫Tc=900K，環(huán)境溫度為300K。

3.1 物理模型

發(fā)動(dòng)機(jī)與擋板空間位置如圖2所示，設(shè)定X軸為對稱軸。發(fā)動(dòng)機(jī)軸線經(jīng)過擋板中心，a=0.1m，角度ω=0°，n=m=300mm，h=0.1mm，Dt=0.6mm，De=4.7mm。

圖2 發(fā)動(dòng)機(jī)與擋板空間位置

3.2 邊界條件

DSMC粒子入口條件：以0.5N發(fā)動(dòng)機(jī)出口截面為粒子入口截面，模擬粒子選取CO2。真空邊界：認(rèn)為粒子逃逸，即粒子經(jīng)過邊界后注銷。粒子入口截面參數(shù)分布如圖3至圖6所示。圖3為粒子入口密度分布圖，圖4為粒子入口溫度分布圖，圖5為粒子入口軸向速度分布圖，圖6為粒子入口徑向速度分布圖。

圖3 粒子入口密度分布圖

圖4 粒子入口溫度分布圖

圖5 粒子入口軸向速度分布圖

圖6 粒子入口徑向速度分布圖

3.3 網(wǎng)格劃分

X，Y，Z方向網(wǎng)格數(shù)目分別為60,120,120，且網(wǎng)格寬度為等比數(shù)列，X方向比例因子為1.035，Y,Z方向網(wǎng)格寬度從中心向兩側(cè)等比增加，比例因子為1.046，Y-Z平面網(wǎng)格視圖如圖7所示。

圖7 網(wǎng)格劃分

4 計(jì)算結(jié)果及分析

使用PWS軟件對羽流撞擊效應(yīng)進(jìn)行計(jì)算。采用64臺(tái)節(jié)點(diǎn)機(jī)并行計(jì)算約2.9h，總共計(jì)算模擬分子數(shù)達(dá)到950萬個(gè)以上，迭代到10,000步時(shí)開始采樣，繼續(xù)統(tǒng)計(jì)20,000步得出羽流流場。流場參數(shù)分布如圖8-圖 11。

圖8 流場壓強(qiáng)分布圖

圖9 流場溫度分布圖

圖10 流場馬赫數(shù)分布圖

由圖8至圖10可以看出，發(fā)動(dòng)機(jī)羽流場的壓強(qiáng)在粒子入口處最高，最大壓強(qiáng)約為950Pa，在軸線和擋板附近一定區(qū)域內(nèi)比較高，隨著氣流向真空中擴(kuò)散而降低。在軸線附近一橢球形區(qū)域內(nèi)，溫度較低，馬赫數(shù)較高；而在這個(gè)橢球形邊緣的附近，由于從發(fā)動(dòng)機(jī)噴出的氣流和反射氣流相遇，造成局部溫度的升高，馬赫數(shù)減小；在橢球形以外的流場中，則以反射氣流為主導(dǎo)，壓強(qiáng)減小，溫度下降，馬赫數(shù)增加。最高溫度約為450K，最大馬赫數(shù)約為12。根據(jù)圖11可以看出，氣流沿著發(fā)動(dòng)機(jī)出口呈放射狀流向擋板，噴流撞擊到擋板上發(fā)生反射，擋板表面附近的氣流向四周反射流走。

擋板壓強(qiáng)分布如圖12所示，由于核心氣流對擋板的撞擊，在發(fā)動(dòng)機(jī)軸線與擋板的交點(diǎn)處壓強(qiáng)較高，以該交點(diǎn)為圓心，板上的壓強(qiáng)分布具有對稱性，且沿徑向逐漸減小。平板中心最高壓強(qiáng)可達(dá)110Pa左右。由此可見，擋板中心處羽流效應(yīng)影響最嚴(yán)重，影響程度由中心沿徑向逐漸減小。

圖12 擋板壓強(qiáng)分布圖

5 結(jié)論

本文采用自開發(fā)的PWS軟件，對推力0.5N的模型發(fā)動(dòng)機(jī)進(jìn)行了羽流撞擊效應(yīng)進(jìn)行DSMS模擬，工質(zhì)為CO2，得到了流場的壓強(qiáng)、溫度、馬赫數(shù)、流線圖；然后，對擋板壓強(qiáng)和羽流效應(yīng)情況進(jìn)行了分析，得到以下結(jié)論：

（1）發(fā)動(dòng)機(jī)羽流場壓強(qiáng)在粒子入口處最高，在軸線和擋板附近一定區(qū)域內(nèi)比較高，隨著氣流向真空中擴(kuò)散而降低。

（2）在軸線附近一個(gè)橢球形區(qū)域內(nèi)，溫度較低，馬赫數(shù)較高；而在這個(gè)橢球形邊緣的附近，由于從發(fā)動(dòng)機(jī)噴出的氣流和反射氣流相遇，造成局部溫度的升高，馬赫數(shù)減小；在橢球形以外的流場中，則以反射氣流為主導(dǎo)，壓強(qiáng)減小，溫度下降，馬赫數(shù)增加。

（3）氣流沿著發(fā)動(dòng)機(jī)出口呈放射狀流向擋板，噴流撞擊到擋板上發(fā)生反射，擋板表面附近的氣流向四周反射流走。

（4）由于核心氣流對擋板的撞擊，在發(fā)動(dòng)機(jī)軸線與擋板的交點(diǎn)處壓強(qiáng)較高，以該交點(diǎn)為圓心，板上的壓強(qiáng)分布具有對稱性，且沿徑向逐漸減小。擋板中心處羽流效應(yīng)影響最嚴(yán)重，影響程度由中心沿徑向逐漸減小。 ◇

［1］ DETTLEFF G,BOETTCHER R D,DANKERT C,KPPPENWALLNER G,LEGGE H.Attitude Control Thruster Plume Flow Modeling and Experiments［J］.Journal of Spacecraft and Rockets,1986,23（5）:476-481.

［2］ALLEGRE J,RAFFIN M,LENGRAND J C.Forces Induced by a Simulated Rocket Exhaust Plume Impinging upon a Flat Plate［A］.Proceedings of the 14th International Symposium on Rarefied Gas Dynamics［C］,1984:287-294.

［3］LEGGE H.Plume impingement forces on Inclined Flat Plates［A］.Proceedings of the 17th International Symposium on Rarefied Gas Dynamics［C］,1990:955-962.

［4］BOYD I D,PENKO P F,MEISSNER D L,DEWITT K J.Experimental and Numerical Investigations of Low-density Nozzle and Plume Flows of Nitrogen ［J］.AIAA Journal,1992,30（10）:2453-2461.

［5］LEGGE H,BOETTCHER R D.Modelling Control Thruster Plume Flow and Impingement［A］.Proceedings of the 13thInternational Symposium on Rarefied Gas Dynamics［C］,1985:983-992.

［6］PARK J H,BAEK S W,KIM J S.DSMC Analysis of the Interaction between Thruster Plume and Satellite Components［R］.AIAA 2002-0794.

［7］蔡國飆，賀碧蛟，PWS軟件應(yīng)用于探月著陸器羽流效應(yīng)數(shù)值模擬研究［A］.航天器環(huán)境工程［J］.2010:18-23.