MF-1彈道設(shè)計(jì)與蒙特卡羅飛行仿真

和爭春, 車 競, 肖涵山, 袁先旭

(1. 空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室, 綿陽 621000; 2. 中國空氣動(dòng)力研究與發(fā)展中心, 綿陽 621000)

0 引 言

模型飛行試驗(yàn)[1]是利用飛行器縮尺或簡化模型在真實(shí)大氣環(huán)境中飛行，進(jìn)行氣動(dòng)力/熱、飛行力學(xué)等科學(xué)問題研究的一種試驗(yàn)手段，是空氣動(dòng)力學(xué)研究三大手段之一。它具有模型飛行不受外界約束，試驗(yàn)高度和速度范圍大，能綜合研究氣動(dòng)力、氣動(dòng)熱、結(jié)構(gòu)、動(dòng)力、飛控等問題的特點(diǎn)，所得結(jié)果更真實(shí)可信。模型飛行試驗(yàn)的獨(dú)特優(yōu)勢和重要價(jià)值歷來都很受西方各航空航天發(fā)達(dá)國家的重視[2-3]，特別是近年來在高超聲速領(lǐng)域陸續(xù)開展了HIFiRE[4-7]、HyShot[8]、HyBoLT[9]、HyCAUSE[10]、EXPERT[11-12]、SHEFEX[13]、Pre-X[14]等飛行試驗(yàn)項(xiàng)目，為邊界層轉(zhuǎn)捩、激波邊界層干擾等高超聲速流動(dòng)機(jī)理深入研究和地面氣動(dòng)力熱預(yù)測方法修正提供了寶貴的真實(shí)飛行數(shù)據(jù)。

MF-1是中國空氣動(dòng)力研究與發(fā)展中心開展的單級火箭助推不分離無控高超聲速空氣動(dòng)力學(xué)飛行試驗(yàn)研究項(xiàng)目，主要目的是通過在試驗(yàn)?zāi)Ｐ吞囟ú课话惭b的溫度、壓力等傳感器，對飛行試驗(yàn)全過程參數(shù)變化歷程實(shí)現(xiàn)測量，為邊界層轉(zhuǎn)捩、激波邊界層干擾機(jī)理研究提供真實(shí)飛行數(shù)據(jù)。MF-1試驗(yàn)飛行器全長6.23 m，最大直徑0.6 m，由試驗(yàn)?zāi)Ｐ汀⒓夐g段、助推器和尾段等組成(圖1)。

本文研究建立了MF-1飛行試驗(yàn)基準(zhǔn)彈道設(shè)計(jì)及拉偏仿真方法模型，開展了MF-1飛行彈道設(shè)計(jì)和蒙特卡羅拉偏仿真，并對設(shè)計(jì)仿真結(jié)果和滾轉(zhuǎn)共振等問題進(jìn)行了分析。

圖1 MF-1試驗(yàn)飛行器外形圖Fig.1 Sketch of MF-1 flight vehicle

1 彈道設(shè)計(jì)及仿真模型

地球模型采用WGS84旋轉(zhuǎn)橢球模型，利用六自由度彈道積分，通過調(diào)整發(fā)射傾角和發(fā)射方位角滿足試驗(yàn)窗口和落點(diǎn)位置要求。

1.1 質(zhì)心運(yùn)動(dòng)模型

(sinφsinγ+cosφsinψcosγ)+gx+axe+axc

(-cosφsinγ+sinφsinψcosγ)+gy+aye+ayc

式中:xd、yd、zd為飛行器在地球固連系中的位置坐標(biāo)；Vxd、Vyd、Vzd為地球固連系速度分量；Nx、Ny、Nz為體軸系氣動(dòng)過載分量；Fx、Fy、Fz為發(fā)動(dòng)機(jī)推力在體軸系中的分量；φ、ψ、γ為俯仰、偏航、滾轉(zhuǎn)姿態(tài)角；gx、gy、gz為重力加速度分量；axe、aye、aze為向心牽連加速度分量；axc、ayc、azc為柯氏加速度分量；g0為常數(shù)9.81。

重力加速度、向心牽連加速度和柯氏加速度是位置坐標(biāo)的函數(shù)，需利用地球物理關(guān)系進(jìn)行推導(dǎo)和計(jì)算，具體公式和推導(dǎo)過程參見文獻(xiàn)[15]。

1.2 繞質(zhì)心轉(zhuǎn)動(dòng)模型

質(zhì)心轉(zhuǎn)動(dòng)的方程是建立在體軸系中的：

其中，ω1=(ωx1ωy1ωz1)T

M1=(Mx1a+Mx1TMy1a+My1TMz1a+Mz1T)T

Mx1a、My1a、Mz1a為繞質(zhì)心的氣動(dòng)力矩分量，Mx1T、My1T、Mz1T為發(fā)動(dòng)機(jī)推力偏斜及偏心產(chǎn)生的繞質(zhì)心的力矩分量。

姿態(tài)角運(yùn)動(dòng)方程為：

式中ωx、ωy、ωz為體軸系角速率分量。

1.3 助推發(fā)動(dòng)機(jī)推力模型

發(fā)動(dòng)機(jī)推力是基于發(fā)動(dòng)機(jī)生產(chǎn)廠家提供的內(nèi)彈道數(shù)據(jù)進(jìn)行高度修正得到：

F=F0+STref(p0-p)

式中:F0為從內(nèi)彈道數(shù)據(jù)差值得到的推力，STref為發(fā)動(dòng)機(jī)噴管出口面積，p0為內(nèi)彈道數(shù)據(jù)對應(yīng)的地面大氣壓，p為飛行器當(dāng)前高度對應(yīng)的大氣壓。

發(fā)動(dòng)機(jī)推力在飛行器體軸系中的分量為：

Fx=FcosθT

Fy=FsinθTcosφT

Fz=FsinθTsinφT(4)

式中:θT為推力線相對于飛行器縱軸的偏斜角度，φT為推力偏斜周向角，即推力偏斜方向。

發(fā)動(dòng)機(jī)推力產(chǎn)生的繞飛行器質(zhì)心的力矩為：

Mx1T=FzyT-FyzT

My1T=FxzT-FzxT

Mz1T=FxyT+FyxT(5)

式中xT、yT、zT為發(fā)動(dòng)機(jī)推力作用點(diǎn)在飛行器體軸系中的坐標(biāo)。

1.4 大氣模型和風(fēng)修正模型

彈道設(shè)計(jì)采用的大氣密度模型為靶場月平均大氣密度統(tǒng)計(jì)模型，每個(gè)月的模型均以與一系列不同高度對應(yīng)的大氣密度當(dāng)月平均值給出。彈道設(shè)計(jì)時(shí)，可以根據(jù)擬實(shí)施發(fā)射的月份選擇對應(yīng)的大氣密度模型。

彈道設(shè)計(jì)采用的風(fēng)場模型為靶場月平均風(fēng)場統(tǒng)計(jì)模型，每個(gè)月的月平均模型均以與一系列不同高度對應(yīng)的風(fēng)速、風(fēng)向當(dāng)月平均值給出。彈道設(shè)計(jì)時(shí)，可以根據(jù)擬實(shí)施發(fā)射的月份選擇對應(yīng)的風(fēng)場模型。

飛行器相對于大氣的速度矢量為：

V=Vd-Vw(6)

其中:Vd為地面固連坐標(biāo)系中的飛行器速度矢量，Vw為地面固連坐標(biāo)系中的風(fēng)速矢量。

在彈道設(shè)計(jì)和仿真中，飛行動(dòng)壓、迎角、側(cè)滑角等均是基于風(fēng)修正后的數(shù)據(jù)得到。

1.5 空氣動(dòng)力模型

彈道設(shè)計(jì)和仿真中的空氣動(dòng)力學(xué)模型是由MF-1氣動(dòng)力團(tuán)隊(duì)提供的，是對于地面風(fēng)洞試驗(yàn)和數(shù)值計(jì)算結(jié)果綜合分析并經(jīng)質(zhì)心換算等工作后得到的，包括無底阻六分量數(shù)據(jù)、底阻修正數(shù)據(jù)、不同高度摩阻修正數(shù)據(jù)、動(dòng)導(dǎo)數(shù)數(shù)據(jù)等。

1.6 發(fā)射架模型

MF-1采用品字型三點(diǎn)同時(shí)離軌下掛式傾斜發(fā)射，滑軌有效長度6.5 m。在離架前采用固定彈道傾角、彈道偏角和飛行器姿態(tài)角的單自由度彈道積分模型，離架后采用1.1節(jié)中的正常六自由度彈道積分模型。

2 基準(zhǔn)飛行彈道設(shè)計(jì)結(jié)果

MF-1基準(zhǔn)彈道設(shè)計(jì)，必須與助推發(fā)動(dòng)機(jī)能力、發(fā)射場和落區(qū)選擇等相結(jié)合來開展，其基本要求是能夠滿足對10～40 km高度試驗(yàn)窗口的覆蓋。而且，由于試驗(yàn)飛行器無控，為了避免彈道太高、無氣動(dòng)阻尼的真空段太長導(dǎo)致的再入大氣層迎角不確定度太大甚至尾部朝前的情況發(fā)生，彈道頂點(diǎn)高度在60～90 km之間為宜。在發(fā)射點(diǎn)和落區(qū)選定后，基準(zhǔn)彈道設(shè)計(jì)的任務(wù)，就是通過設(shè)計(jì)發(fā)射架俯仰角和發(fā)射方位角，使得試驗(yàn)飛行器落于所選落點(diǎn)附近。

本文基準(zhǔn)彈道設(shè)計(jì)算例采用的大氣密度模型和風(fēng)場模型來自于靶場12月份的月平均統(tǒng)計(jì)數(shù)據(jù)；固體火箭發(fā)動(dòng)機(jī)內(nèi)彈道數(shù)據(jù)采用生產(chǎn)廠家提供的-10 ℃發(fā)動(dòng)機(jī)內(nèi)彈道數(shù)據(jù)，發(fā)動(dòng)機(jī)推力作用線沿發(fā)動(dòng)機(jī)縱軸，無橫移無偏斜。

彈道設(shè)計(jì)采用的試驗(yàn)飛行器滿載和空載質(zhì)心位置及質(zhì)量慣量見表1，其中的結(jié)構(gòu)坐標(biāo)系定義為：錐理論頂點(diǎn)為原點(diǎn)，X軸向后，Y軸向上，Z軸與X、Y構(gòu)成右手系；體軸系定義為：質(zhì)心為原點(diǎn)，X軸向前，Y軸向上，Z軸與X、Y構(gòu)成右手系。

基于以上條件，對MF-1試驗(yàn)飛行器進(jìn)行基準(zhǔn)彈道設(shè)計(jì)得到的發(fā)射架仰角為62.148°、發(fā)射方位角為267.410°。圖2給出了MF-1基準(zhǔn)彈道設(shè)計(jì)結(jié)果的飛行剖面示意圖，表2給出了基準(zhǔn)彈道特征點(diǎn)參數(shù)，圖3給出了基準(zhǔn)彈道參數(shù)曲線。

表1 MF-1質(zhì)心位置和質(zhì)量慣量Table 1 MF-1 gravity center, mass and moment of inertia

圖2 MF-1飛行剖面示意圖Fig.2 Sketch of the MF-1 flight profile

特征點(diǎn)時(shí)間/s高度/km縱向航程/km速度/(m·s-1)馬赫數(shù)動(dòng)壓/kPa點(diǎn)火0.01.0350.00.00.00.0離架0.4321.0410.00337.140.110.84進(jìn)試驗(yàn)窗口15.5210.006.74312614.21332.2關(guān)機(jī)點(diǎn)19.5714.5810.4816455.58279.1出試驗(yàn)窗口44.6840.0033.8412433.923.098彈道頂點(diǎn)133.777.91112.8894.53.130.009落地點(diǎn)269.01.553224.1410.21.2394.30

從設(shè)計(jì)結(jié)果可見，發(fā)射離架后由于彈道下沉導(dǎo)致的最大飛行迎角不超過2°，此后飛行迎角迅速收斂到0°附近。隨著高度迅速升高大氣密度越來越稀薄，迎角振幅逐漸發(fā)散，最大振幅不超過0.2°，進(jìn)入下降段后迎角振幅又向0°收斂。整個(gè)飛行過程中最大動(dòng)壓約330 kPa，最大軸向過載14左右。經(jīng)MF-1氣動(dòng)力/熱研究團(tuán)隊(duì)評估，該彈道滿足MF-1氣動(dòng)研究試驗(yàn)窗口需求，氣動(dòng)熱環(huán)境也屬于可承受范圍。

圖3 MF-1基準(zhǔn)彈道參數(shù)曲線Fig.3 Curves of MF-1 standard trajectory parameters

3 蒙特卡羅拉偏仿真結(jié)果

MF-1蒙特卡羅拉偏仿真所采用的偏差因素及其分布，是由包括總體、氣動(dòng)、動(dòng)力、發(fā)射架等專業(yè)在內(nèi)的來自總體單位、型號部門、院校等單位的專家們經(jīng)過多次會(huì)議討論確定的，動(dòng)力、發(fā)射架等拉偏依據(jù)主要來自產(chǎn)品給出的精度數(shù)據(jù)，所有的小不對稱氣動(dòng)參數(shù)主要來自于尾翼安裝偏差影響的氣動(dòng)力計(jì)算結(jié)果，其它拉偏參數(shù)也借鑒了多種類似飛行器型號經(jīng)驗(yàn)。表3給出了基于設(shè)計(jì)計(jì)算預(yù)估數(shù)據(jù)的MF-1六自由度彈道拉偏仿真偏差因素及取值，其中的小不對稱參數(shù)按每片尾翼都存在7.5′安裝角偏差估算得到，比如，估算滾轉(zhuǎn)小不對稱氣動(dòng)力時(shí)，假設(shè)四片尾翼安裝角的方向均可產(chǎn)生同一個(gè)方向的滾轉(zhuǎn)力矩；估算俯仰小不對稱氣動(dòng)力時(shí)，假設(shè)四片尾翼安裝角的方向均可產(chǎn)生同一個(gè)方向的俯仰力矩和法向力；偏航方向依此類推。圖4給出了10 000條彈道的蒙特卡羅拉偏結(jié)果彈道參數(shù)散布圖，其中的彈道參數(shù)歷程散布是按照每條彈道起始時(shí)刻在0～10 s間隨機(jī)均勻分布、每10 s一個(gè)彈道點(diǎn)的方式進(jìn)行輸出和作圖的。從結(jié)果可見，在各種隨機(jī)拉偏因素的影響下，彈道參數(shù)和落點(diǎn)參數(shù)均有較大散布。主要參數(shù)散布情況統(tǒng)計(jì)結(jié)果為：落點(diǎn)縱向航程均值223.75 km，均方差15.93 km；落點(diǎn)側(cè)向航程均值8.14 km，均方差6.93 km；彈道頂點(diǎn)高度均值75.84 km，均方差6.09 km；彈道最大滾轉(zhuǎn)角速率均值9.34°/s，均方差473.9°/s。

圖4 基于預(yù)估數(shù)據(jù)的MF-1蒙特卡羅飛行仿真結(jié)果Fig.4 Monte Carlo flight simulation results for MF-1 based on predicted data

偏差類別偏差因素偏差分布?xì)鈩?dòng)力偏差小不對稱法向力系數(shù)正態(tài),均值0,均方差0.0003小不對稱側(cè)向力系數(shù)正態(tài),均值0,均方差0.0003小不對稱滾轉(zhuǎn)力矩系數(shù)正態(tài),均值0,均方差0.001小不對稱俯仰力矩系數(shù)正態(tài),均值0,均方差0.0003小不對稱偏航力矩系數(shù)正態(tài),均值0,均方差0.0003軸向力系數(shù)百分比偏差正態(tài),均值0,均方差5%法向力系數(shù)百分比偏差正態(tài),均值0,均方差3%側(cè)向力系數(shù)百分比偏差正態(tài),均值0,均方差3%滾轉(zhuǎn)力矩系數(shù)百分比偏差正態(tài),均值0,均方差10%偏航力矩系數(shù)百分比偏差正態(tài),均值0,均方差7%俯仰力矩系數(shù)百分比偏差正態(tài),均值0,均方差7%三個(gè)方向動(dòng)導(dǎo)數(shù)50%至200%均勻分布質(zhì)量慣量特性偏差起飛質(zhì)量百分比偏差正態(tài),均值0,均方差0.1%三個(gè)方向轉(zhuǎn)動(dòng)慣量百分比偏差正態(tài),均值0,均方差7%質(zhì)心位置x向偏差正態(tài),均值0,均方差6 mm質(zhì)心位置y、z向偏差正態(tài),均值0,均方差5/3 mm發(fā)動(dòng)機(jī)偏差發(fā)動(dòng)機(jī)推力線偏斜偏斜角的量值為均值0,均方差(5/3)′正態(tài)分布偏斜方向?yàn)?～180°均勻分布發(fā)動(dòng)機(jī)推力線橫移橫移量為均值0,均方差0.5 mm正態(tài)分布橫移方向?yàn)?～180°均勻分布發(fā)動(dòng)機(jī)總沖百分比偏差正態(tài),均值0,均方差0.2%發(fā)動(dòng)機(jī)高低溫影響發(fā)動(dòng)機(jī)內(nèi)彈道隨機(jī)在-15 ℃、-10 ℃和-5 ℃三套數(shù)據(jù)中抽取氣象參數(shù)偏差和發(fā)射角偏差發(fā)射傾角偏差正態(tài),均值0,均方差(0.05/3)°發(fā)射偏角偏差正態(tài),均值0,均方差(0.1/3)°大氣密度模型在靶場12月份平均大氣密度與標(biāo)準(zhǔn)大氣的百分比偏差的基礎(chǔ)上,對偏差乘以0～200%的均勻分布風(fēng)場采用靶場12月份的月平均風(fēng)速和風(fēng)向、對風(fēng)速乘以80%至120%均勻分布,再疊加上均值0、均方差1 m/s的正態(tài)分布、風(fēng)向0°～180°的均勻分布的隨機(jī)風(fēng)

另外，從圖4的高度-航程散布可見，有幾條彈道的高度明顯低于其它彈道，航程也要近得多，經(jīng)分析，這主要是由于在飛行過程中發(fā)生了滾轉(zhuǎn)共振[16]連鎖，氣動(dòng)阻力顯著增加而導(dǎo)致的。飛行器在飛行過程中發(fā)生滾轉(zhuǎn)共振連鎖的原因，是由于在某些條件下其滾轉(zhuǎn)頻率與氣動(dòng)俯仰偏航頻率在較長時(shí)間段內(nèi)一直比較接近，從而引發(fā)共振，使得飛行器的飛行總迎角顯著增大，橫法向過載顯著增加，甚至可能造成飛行器解體。以圖4中彈道高度最低的彈道(蒙特卡羅仿真的第9264條彈道)為例，圖5給出了該條彈道的滾轉(zhuǎn)頻率、俯仰頻率、滾轉(zhuǎn)與俯仰頻率之比，以及迎角側(cè)滑角的時(shí)間歷程。從圖可見，在關(guān)機(jī)點(diǎn)附近，俯仰頻率與滾轉(zhuǎn)頻率趨于接近，二者之比迅速達(dá)到1附近，隨著高度的增加和動(dòng)壓的減小，這兩個(gè)頻率的變化趨勢也非常接近，導(dǎo)致二者之比的量級在較長時(shí)間內(nèi)維持在1附近。這也使得側(cè)滑角在關(guān)機(jī)點(diǎn)前后迅速增大到-30°～-40°，并維持了很長時(shí)間，屬于典型的較長時(shí)間的滾轉(zhuǎn)共振連鎖。后來，隨著MF-1試驗(yàn)飛行器重新進(jìn)入稠密大氣內(nèi)飛行，俯仰頻率迅速增加，與此同時(shí)滾轉(zhuǎn)頻率逐漸減小，滾轉(zhuǎn)共振的條件不復(fù)存在，側(cè)滑角也就迅速恢復(fù)到0°附近并振蕩收斂。由此可見，關(guān)機(jī)點(diǎn)附近滾轉(zhuǎn)頻率與氣動(dòng)俯仰偏航頻率是否接近1，可以近似作為MF-1可能滾轉(zhuǎn)共振出現(xiàn)的判據(jù)。近似以滾轉(zhuǎn)頻率/俯仰頻率的比值在0.7～1.5之間作為可能發(fā)生滾轉(zhuǎn)共振的條件，從上述的拉偏仿真結(jié)果得到的可能發(fā)生滾轉(zhuǎn)共振的概率為6.3%。

圖5 第9264條彈道滾轉(zhuǎn)共振分析曲線Fig.5 Curves of simulated No. 9264 trajectory for MF-1 roll resonance analysis

在MF-1完成制造總裝出廠前，由于有了更為可信的氣動(dòng)力風(fēng)洞試驗(yàn)數(shù)據(jù)以及合練產(chǎn)品的安裝精度以及質(zhì)量特性實(shí)測數(shù)據(jù)，使得用于飛行分析的數(shù)據(jù)可信度有了明顯提高。特別是尾翼安裝后的實(shí)測結(jié)果表明，安裝精度遠(yuǎn)遠(yuǎn)高于之前的估計(jì)，四片尾翼安裝角代數(shù)和僅為5′量級，這也必將導(dǎo)致拉偏后試飛器的滾轉(zhuǎn)角速率散布顯降低。在此基礎(chǔ)上，進(jìn)行的10 000條拉偏彈道仿真結(jié)果表明，沒有一條彈道滿足上面的可能發(fā)生滾轉(zhuǎn)共振的條件，說明基于總裝后實(shí)測數(shù)據(jù)條件下，MF-1試飛器在上升段發(fā)生滾轉(zhuǎn)共振的概率接近0，試飛器不會(huì)發(fā)生滾轉(zhuǎn)共振而導(dǎo)致解體風(fēng)險(xiǎn)。圖6給出了基于總裝后實(shí)測數(shù)據(jù)的10 000條蒙特卡羅拉偏彈道的高度-時(shí)間分布。

圖6 基于出廠測量數(shù)據(jù)的MF-1蒙特卡羅飛行仿真結(jié)果Fig.6 Monte Carlo flight simulation results for MF-1 based on final factory test data

4 結(jié) 論

研究建立了MF-1試驗(yàn)飛行器彈道設(shè)計(jì)與拉偏仿真數(shù)學(xué)模型，并開展了飛行試驗(yàn)彈道設(shè)計(jì)與拉偏仿真分析。結(jié)果表明：

(1) 本文建立的彈道設(shè)計(jì)和仿真模型，可用于傾斜發(fā)射無控飛行器的彈道設(shè)計(jì)和偏差條件下的彈道及落點(diǎn)散布特性研究；

(2) 對于MF-1這類無控飛行器，在偏差條件下可能在飛行過程中發(fā)生滾轉(zhuǎn)共振，導(dǎo)致飛行風(fēng)險(xiǎn)增加；

(3) MF-1嚴(yán)格控制了尾翼加工安裝精度后，飛行過程中不會(huì)發(fā)生滾轉(zhuǎn)共振。