唐林卡，李曉軒，孫朝翔，嚴(yán)東升，萬(wàn)金杰

（北京航天長(zhǎng)征飛行器研究所，北京，100076）

0 引 言

衛(wèi)星等空間在軌飛行器在軌運(yùn)行期間需要利用自帶動(dòng)力系統(tǒng)進(jìn)行軌道維護(hù)；某些再入飛行器，在被動(dòng)段飛行期間需要利用小發(fā)動(dòng)機(jī)維持其姿態(tài)穩(wěn)定。此類發(fā)動(dòng)機(jī)在真空或接近真空環(huán)境下工作時(shí)，都會(huì)在噴口處產(chǎn)生膨脹羽流，該羽流流場(chǎng)特性與發(fā)動(dòng)機(jī)在大氣環(huán)境下工作時(shí)的噴流流場(chǎng)特性顯著不同，造成了對(duì)飛行器力、熱環(huán)境影響的天地差異，對(duì)飛行器的環(huán)境適應(yīng)性及飛行結(jié)果分析產(chǎn)生顯著影響[1]。本文主要結(jié)合羽流流場(chǎng)特性數(shù)值仿真計(jì)算和地面試驗(yàn)結(jié)果，分析羽流對(duì)飛行器造成的影響，給出相應(yīng)的工程解決方案。

1 羽流對(duì)飛行器影響

1.1 羽流對(duì)飛行器彈道影響

某飛行器在底部中心位置安裝了小火箭發(fā)動(dòng)機(jī)，通過(guò)其兩側(cè)噴管產(chǎn)生推力使其旋轉(zhuǎn)，如圖1所示。

在飛行中，通過(guò)遙測(cè)數(shù)據(jù)發(fā)現(xiàn)起旋發(fā)動(dòng)機(jī)工作過(guò)程中除了產(chǎn)生滾轉(zhuǎn)角速率外，還會(huì)影響飛行器的軸向速度，如圖2所示。數(shù)值仿真結(jié)果表明，起旋發(fā)動(dòng)機(jī)工作時(shí)，真空環(huán)境下迅速膨脹的燃?xì)庥鹆髯饔迷陲w行器底部，在部分區(qū)域產(chǎn)生壓力，如圖3所示，該壓力形成一定軸向過(guò)載，使飛行器的軸向速度發(fā)生變化，造成飛行軌跡偏差。

圖2 起旋發(fā)動(dòng)機(jī)工作時(shí)的軸向過(guò)載Fig.2 Axial Overload of Spin Engine during Operation

圖3 數(shù)值仿真得到的羽流壓力分布Fig.3 Plume Pressure Distribution Obtained by Numerical Simulation

1.2 羽流對(duì)飛行器熱環(huán)境影響

某飛行器在其底部安裝了小型固體火箭，以便在空間飛行時(shí)產(chǎn)生軸向速度和旋轉(zhuǎn)速度。在大氣環(huán)境中進(jìn)行點(diǎn)火試驗(yàn)時(shí)，發(fā)動(dòng)機(jī)噴流未產(chǎn)生膨脹羽流，試驗(yàn)未發(fā)生故障。后續(xù)在不同容積的真空罐中又分別進(jìn)行了點(diǎn)火試驗(yàn)。在大容積（約1000 m3）真空罐中，發(fā)動(dòng)機(jī)羽流自由膨脹，使得發(fā)動(dòng)機(jī)噴口附近非金屬蒙皮發(fā)生燒穿故障，如圖4所示。在小容積（約100 m3）真空罐中，發(fā)動(dòng)機(jī)工作產(chǎn)生的燃?xì)馐拐婵展迌?nèi)壓力迅速上升，限制了羽流的自由膨脹，也沒(méi)有發(fā)生蒙皮被燒穿現(xiàn)象。該現(xiàn)象說(shuō)明發(fā)動(dòng)機(jī)羽流現(xiàn)象對(duì)環(huán)境壓力影響十分敏感，對(duì)地面試驗(yàn)條件有較為嚴(yán)格的要求。

圖4 發(fā)動(dòng)機(jī)真空羽流導(dǎo)致的表面燒穿示意Fig.4 Surface Burn-through Caused by Engine Vacuum Plume

2 羽流流場(chǎng)特性仿真

2.1 DSMC方法

DSMC方法[2]是目前用于高空環(huán)境羽流干擾流動(dòng)的主要分析方法。

DSMC方法的基本思想是用有限個(gè)仿真分子代替真實(shí)氣體分子，并在計(jì)算機(jī)中存儲(chǔ)仿真分子的位置坐標(biāo)、速度分量和內(nèi)能，其值隨仿真分子的運(yùn)動(dòng)、與邊界的作用以及仿真分子之間的碰撞而改變，最后通過(guò)統(tǒng)計(jì)網(wǎng)格內(nèi)仿真分子的運(yùn)動(dòng)狀態(tài)實(shí)現(xiàn)對(duì)真實(shí)氣體流動(dòng)問(wèn)題的模擬。計(jì)算時(shí)除考慮平動(dòng)能外, 還考慮了內(nèi)部能量。分子間相互作用模型采用VHS模型，物面邊界條件采用完全漫反射條件，平動(dòng)能與內(nèi)能的能量交換采用L-B模型，按照Bird的能量按自由度分配原則采用取舍法進(jìn)行抽樣分配[3]。

平動(dòng)能和轉(zhuǎn)動(dòng)能之間的抽樣符合P.S.Larsen和C.Borgnakke提出的L-B模型分布：

式中H為兩分子總轉(zhuǎn)動(dòng)能rot,abε在總能量Mε中所占比例，εM=εrota,b+εrel，εrel=0.5μg2為相對(duì)平動(dòng)能。采用取舍法抽樣得到H值后，可得碰撞后的分子總轉(zhuǎn)動(dòng)能和相對(duì)平動(dòng)能為

分子的轉(zhuǎn)動(dòng)能按照等概率原則進(jìn)行分配，則碰撞后的分子轉(zhuǎn)動(dòng)能為

2.2 發(fā)動(dòng)機(jī)噴管流場(chǎng)仿真

首先使用求解納維-斯托克斯方程（Navier-Stokes Equations，NS）方法獲得發(fā)動(dòng)機(jī)噴管內(nèi)流場(chǎng)[4]，然后取發(fā)動(dòng)機(jī)噴管出口截面流場(chǎng)參數(shù)，作為發(fā)動(dòng)機(jī)羽流流場(chǎng)的入口條件[5,6]，如圖5所示?？梢钥吹接捎趪姽艹隹?喉道面積比較小，出口馬赫數(shù)在3.5左右，壓力為60 000~70 000 Pa，速度約為2000 m/s。相比外部低壓靜止環(huán)境，噴管流動(dòng)屬于高壓高速流動(dòng)。

圖5 NS方程求解噴管內(nèi)流場(chǎng)Fig.5 Solution of Flow Field in the Nozzle by NS Equation

續(xù)圖5

圖6給出環(huán)境壓力為0時(shí)的噴管流場(chǎng)，圖7給出環(huán)境壓力為80 Pa時(shí)噴管流場(chǎng)，可以看到，外部壓力為0時(shí)，噴流流出噴管后發(fā)生顯著的膨脹和散射效應(yīng)。而當(dāng)外部環(huán)境壓力較高時(shí)，噴流從噴管流出后存在明顯的流動(dòng)壓縮，并存在明顯的核心流動(dòng)區(qū)域。

圖6 噴管流場(chǎng)云圖（P=0 Pa）Fig.6 Cloud Chart of Nozzle Flow Field (P=0 Pa)

圖7 噴管流場(chǎng)云圖（P=80 Pa）Fig.7 Cloud Chart of Nozzle Flow Field (P=80 Pa)

續(xù)圖7

取H0=0.5ρu2，H1=0.5u2兩個(gè)參數(shù)分別表征流場(chǎng)總能和氣體動(dòng)能，其中，ρ為流體密度，u為流動(dòng)速度。圖8給出兩個(gè)壓力條件下的噴流流場(chǎng)云圖，可以清楚看到在流場(chǎng)核心區(qū)內(nèi)表征總能的H0在不同壓力環(huán)境下基本一致，而動(dòng)能則存在明顯差異。這說(shuō)明當(dāng)環(huán)境壓力差異時(shí)，流場(chǎng)結(jié)構(gòu)差異顯著，但是對(duì)于流場(chǎng)核心區(qū)，總體差異不大。其中豎線處為距噴口300 mm位置，對(duì)于該位置，不同外部環(huán)境壓力時(shí)，噴流核心區(qū)（Z：-150~150 mm）流速和壓力基本相當(dāng)，只在核心區(qū)外存在明顯差異，這種差異主要由于外部環(huán)境壓力不同導(dǎo)致。

圖8 噴流流場(chǎng)參數(shù)云圖Fig.8 Cloud Chart of Jet Flow Field Parameters

3 羽流流場(chǎng)特性的試驗(yàn)驗(yàn)證

發(fā)動(dòng)機(jī)羽流流場(chǎng)特性試驗(yàn)可以在具備一定容積和一定真空度的真空罐內(nèi)進(jìn)行。采用高精度測(cè)力天平動(dòng)態(tài)測(cè)量羽流作用在飛行器表面產(chǎn)生的力和力矩，采用測(cè)壓傳感器測(cè)量環(huán)境壓強(qiáng)、壁面壓力分布、發(fā)動(dòng)機(jī)燃燒室總壓等數(shù)據(jù)，采用紅外相機(jī)記錄發(fā)動(dòng)機(jī)羽流流場(chǎng)，如圖9、圖10所示。試驗(yàn)結(jié)果可用于校核直接模擬蒙特卡羅法（Direct Simulation Monte Carlo，DSMC）計(jì)算結(jié)果或積累發(fā)動(dòng)機(jī)總沖及羽流效應(yīng)散差子樣。

圖9 力、力矩、壓強(qiáng)測(cè)量結(jié)果Fig.9 Measurement Results of Force, Moment and Pressure

圖10 羽流紅外圖像測(cè)量結(jié)果Fig.10 Measurement Results of Plume Infrared Image

4 羽流影響的工程解決方法

4.1 羽流影響的系統(tǒng)級(jí)預(yù)示與修正

4.1.1 預(yù) 示

如前文所述，發(fā)動(dòng)機(jī)羽流在真空中的膨脹對(duì)飛行器的力學(xué)和熱環(huán)境都會(huì)帶來(lái)較大影響，必須予以高度關(guān)注。因此需要準(zhǔn)確預(yù)示真空羽流的流場(chǎng)特性，確定作用于飛行器上的熱流、加熱量、壓強(qiáng)、過(guò)載等參數(shù)，然后根據(jù)實(shí)際情況予以防護(hù)或修正。目前可以采用以下3種預(yù)示方法。

a）數(shù)值仿真。

真空環(huán)境的羽流流場(chǎng)數(shù)值仿真一般采用第2節(jié)所述的DSMC方法，需要相對(duì)復(fù)雜的專業(yè)程序并花費(fèi)較多的計(jì)算機(jī)時(shí)。也可以基于一些假定條件，采用求解NS方程的連續(xù)流方法進(jìn)行近似估算，但誤差相對(duì)較大。

b）地面試驗(yàn)。

可以采用第3節(jié)所述的試驗(yàn)方法對(duì)真空環(huán)境的羽流流場(chǎng)和作用效果進(jìn)行測(cè)量，進(jìn)而評(píng)估羽流影響。開(kāi)展試驗(yàn)時(shí)，試驗(yàn)環(huán)境、試驗(yàn)產(chǎn)品、試驗(yàn)工裝、測(cè)量手段要盡可能模擬真實(shí)飛行環(huán)境，采集足夠的有效試驗(yàn)子樣，以保證統(tǒng)計(jì)結(jié)果可用于羽流影響的準(zhǔn)確評(píng)估。

c）飛行試驗(yàn)測(cè)量統(tǒng)計(jì)。

飛行試驗(yàn)與地面試驗(yàn)相比，具有真實(shí)的環(huán)境條件，但測(cè)量手段相對(duì)有限，且試驗(yàn)成本更高。目前主要可以通過(guò)高精度過(guò)載、角速率測(cè)量和壁面及燃燒室的壓力測(cè)量來(lái)評(píng)估羽流作用效果，需要重點(diǎn)關(guān)注測(cè)量設(shè)備的精度和采樣率。和地面試驗(yàn)一樣，飛行試驗(yàn)同樣需要一定數(shù)量的有效子樣，才能通過(guò)統(tǒng)計(jì)給出對(duì)羽流影響的準(zhǔn)確評(píng)估。

4.1.2 修 正

采用上述3種方法確定發(fā)動(dòng)機(jī)真空羽流流場(chǎng)的力、熱參數(shù)系統(tǒng)值后，可以系統(tǒng)修正其對(duì)飛行器飛行環(huán)境和飛行器飛行軌道的影響。例如對(duì)于圖5所示的羽流附加熱環(huán)境對(duì)蒙皮的破壞，可以采取提高蒙皮整體抗熱燒蝕性能，或在確定真空膨脹羽流流場(chǎng)的溫度場(chǎng)參數(shù)后，采取局部防熱燒蝕加強(qiáng)措施解決蒙皮的抗發(fā)動(dòng)機(jī)羽流防熱問(wèn)題。對(duì)于圖4所示的發(fā)動(dòng)機(jī)羽流對(duì)飛行器產(chǎn)生軸向力沖量和速度增量，進(jìn)而產(chǎn)生飛行軌跡偏差的解決方法一般有如下兩種。

a）基于零軌跡線理論的方法。

軌跡曲線L與分離點(diǎn)位置x、y、z，分離點(diǎn)速度Vx、Vy、Vz，有以下關(guān)系：

小偏差線性化：

軌跡的偏導(dǎo)數(shù)簡(jiǎn)化為

令LΔ=0，則零軌跡線對(duì)應(yīng)的俯仰角為

飛行器的起旋發(fā)動(dòng)機(jī)在該俯仰角條件下工作，便可消除發(fā)動(dòng)機(jī)羽流膨脹產(chǎn)生的軸向速度增量對(duì)飛行器軌道的影響。同理，對(duì)于橫向速度增量也可參照此解決。

b）基于修正諸元裝訂參數(shù)的方法。

根據(jù)飛行器結(jié)構(gòu)設(shè)計(jì)狀態(tài)和發(fā)動(dòng)機(jī)工作時(shí)的燃燒室理論參數(shù)，按本節(jié)前述的數(shù)學(xué)仿真法、地面試驗(yàn)測(cè)量統(tǒng)計(jì)法或飛行試驗(yàn)測(cè)量統(tǒng)計(jì)法求得飛行器軸向力或軸向沖量，由下式得到軸向速度增量：

式中xF為發(fā)動(dòng)機(jī)羽流作用于飛行器的軸向力均值；T為發(fā)動(dòng)機(jī)工作時(shí)間均值；M為飛行器質(zhì)量均值；I為發(fā)動(dòng)機(jī)總沖均值；F為發(fā)動(dòng)機(jī)推力均值。

在運(yùn)載器諸元計(jì)算準(zhǔn)備時(shí)，考慮此軸向速度增量值，即可修正其帶來(lái)的軌跡變化量。同理，對(duì)于橫向速度增量也可參照此解決。

4.2 羽流影響的隨機(jī)誤差評(píng)估

按照上節(jié)所述的方法可對(duì)羽流流場(chǎng)影響進(jìn)行預(yù)示和系統(tǒng)修正，但在實(shí)際飛行中，由于發(fā)動(dòng)機(jī)結(jié)構(gòu)制造和安裝、藥柱特性、飛行器質(zhì)量特性等因素均是存在隨機(jī)誤差的，此外還有流場(chǎng)數(shù)值計(jì)算、地面試驗(yàn)或飛行試驗(yàn)測(cè)量誤差等。因此采用上述方法進(jìn)行系統(tǒng)修正后，各項(xiàng)隨機(jī)誤差帶來(lái)的影響依然存在且不可忽視，其綜合各影響因素后的總量值約為系統(tǒng)影響量值的10%。因此進(jìn)行羽流影響隨機(jī)誤差的辨識(shí)和控制是非常重要復(fù)雜的工作。限于篇幅，本文僅就依據(jù)DSMC數(shù)值計(jì)算方法進(jìn)行預(yù)示和系統(tǒng)修正后，相關(guān)隨機(jī)誤差的辨識(shí)進(jìn)行簡(jiǎn)要描述。

當(dāng)采用DSMC數(shù)值計(jì)算方法，對(duì)實(shí)際飛行羽流影響造成的速度增量Vx進(jìn)行預(yù)示時(shí)，認(rèn)為引起隨機(jī)偏差ΔVx的因素主要有4項(xiàng)：發(fā)動(dòng)機(jī)總沖偏差σI，發(fā)動(dòng)機(jī)制造和安裝偏差σxyz，DSMC仿真方法誤差σDSMC，飛行器質(zhì)量特性偏差Mσ。認(rèn)為該4項(xiàng)偏差相互獨(dú)立，則有：

在具體工程實(shí)踐中，根據(jù)發(fā)動(dòng)機(jī)及飛行器的具體設(shè)計(jì)制造情況可對(duì)式（6）~（7）中的1、2、4項(xiàng)誤差進(jìn)行控制優(yōu)化，同時(shí)通過(guò)地面試驗(yàn)和飛行試驗(yàn)驗(yàn)證來(lái)修正DSMC計(jì)算基準(zhǔn)以減小其方法誤差，從而盡量減小總隨機(jī)誤差。

5 結(jié) 論

真空環(huán)境中，火箭發(fā)動(dòng)機(jī)噴流顯著膨脹形成羽流，會(huì)對(duì)飛行器力熱環(huán)境造成影響，需要通過(guò)DSMC仿真、地面試驗(yàn)等方法，評(píng)估其對(duì)彈道精度、熱防護(hù)方案等方面的影響，必要時(shí)進(jìn)行彈道系統(tǒng)修正和熱防護(hù)方案改進(jìn)。在彈道系統(tǒng)修正時(shí)必須辨識(shí)修正方法對(duì)應(yīng)的隨機(jī)誤差，分析其偏差構(gòu)成和影響程度，采取對(duì)應(yīng)的控制措施，以保證飛行器的性能指標(biāo)符合設(shè)計(jì)要求。