席柯，閻超*，黃宇，王文，袁武

(1.北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京100191;2.中國科學(xué)院 計(jì)算機(jī)網(wǎng)絡(luò)信息中心，北京100190)

飛行器動(dòng)態(tài)穩(wěn)定性導(dǎo)數(shù)(工程上常稱為“動(dòng)導(dǎo)數(shù)”)是飛行器動(dòng)態(tài)氣動(dòng)特性設(shè)計(jì)、彈道設(shè)計(jì)及動(dòng)態(tài)穩(wěn)定性分析中的關(guān)鍵參數(shù).動(dòng)導(dǎo)數(shù)是分析飛行器動(dòng)態(tài)品質(zhì)和設(shè)計(jì)操控系統(tǒng)的原始?xì)鈩?dòng)數(shù)據(jù)，動(dòng)導(dǎo)數(shù)的準(zhǔn)確評(píng)估對(duì)飛行器設(shè)計(jì)和飛行都有重要的意義.

傳統(tǒng)上，采用繞定軸俯仰振動(dòng)(強(qiáng)迫振動(dòng)或自由振動(dòng))方法得到的是俯仰阻尼導(dǎo)數(shù)的組合形式Cmq+，即直接阻尼導(dǎo)數(shù)(或稱旋轉(zhuǎn)導(dǎo)數(shù))Cmq和洗流時(shí)差導(dǎo)數(shù)(或稱加速度導(dǎo)數(shù))的組合.傳統(tǒng)位勢流理論認(rèn)為:時(shí)差導(dǎo)數(shù)在組合導(dǎo)數(shù)中所占的比例較小，采用Cmq+代替Cmq誤差不大.但是，特別是隨著現(xiàn)代科技的進(jìn)步和戰(zhàn)爭需求的升級(jí)，出現(xiàn)了具有大迎角、超機(jī)動(dòng)飛行能力的飛行器，其外形設(shè)計(jì)及運(yùn)動(dòng)方式都比傳統(tǒng)飛行器復(fù)雜，流動(dòng)的非線性及非定常性也更強(qiáng).上述結(jié)論不一定成立，比如大升力體飛行器的時(shí)差導(dǎo)數(shù)在組合導(dǎo)數(shù)中所占比例較大［1］(可達(dá)40% ～50%)，并且實(shí)際飛行中，飛行器的俯仰角速度和迎角變化率并不總是相等，因此需要分別確定組合導(dǎo)數(shù)的各個(gè)分量.

目前風(fēng)洞試驗(yàn)及計(jì)算對(duì)俯仰阻尼導(dǎo)數(shù)分量的研究并不多.一般風(fēng)洞提供的是直勻流，不能進(jìn)行直接阻尼導(dǎo)數(shù)Cmq的試驗(yàn).國內(nèi)，任玉新［2］求解流動(dòng)控制方程的敏感性方程，直接得到阻尼導(dǎo)數(shù)的各個(gè)分量，并研究了減縮頻率的影響，劉偉等人［3］采用強(qiáng)迫沉浮運(yùn)動(dòng)形式，數(shù)值模擬了高超聲速返回艙及彈道外形(HBS)的加速度導(dǎo)數(shù)，研究指出在組合導(dǎo)數(shù)小于零的情況下，加速度導(dǎo)數(shù)的符號(hào)可能大于零，起負(fù)阻尼作用，屬于動(dòng)不穩(wěn)定情況.孫海生［4］在低速風(fēng)洞中進(jìn)行了加速度導(dǎo)數(shù)的試驗(yàn)，指出在α＜30°范圍內(nèi)，加速度導(dǎo)數(shù)大約是組合導(dǎo)數(shù)的30% ～50%.國外，美國ARL(Army Research Laboratory)的 Weinacht［5-6］及英國的 Qin等人［7］采用定常方法求解錐轉(zhuǎn)運(yùn)動(dòng)及螺旋運(yùn)動(dòng)，獲得了ANSR(Army-Navy Spinner Rocket)和裙?fàn)畎l(fā)射體的直接阻尼導(dǎo)數(shù)Cmq及時(shí)差導(dǎo)數(shù)，但該方法只適用于軸對(duì)稱外形.

本文結(jié)合計(jì)算流體力學(xué)(CFD)技術(shù)開展了飛行器俯仰阻尼導(dǎo)數(shù)各個(gè)分量的數(shù)值模擬研究.采用準(zhǔn)定常方法直接獲得直接阻尼導(dǎo)數(shù)Cmq，求解強(qiáng)迫沉浮運(yùn)動(dòng)獲得洗流時(shí)差導(dǎo)數(shù)，求解強(qiáng)迫俯仰振動(dòng)獲得俯仰阻尼組合導(dǎo)數(shù)Cmq+.通過HBS及基本帶翼導(dǎo)彈(Basic Finner)標(biāo)模進(jìn)行了數(shù)值方法驗(yàn)證及研究，探討了質(zhì)心位置對(duì)各阻尼導(dǎo)數(shù)的影響以及大攻角下各阻尼導(dǎo)數(shù)的變化規(guī)律，在此基礎(chǔ)上計(jì)算了日本的再入飛行器Hyflex的各阻尼導(dǎo)數(shù).

1 計(jì)算方法

1.1 流動(dòng)控制方程

流動(dòng)控制方程為三維非定常可壓縮Navier-Stokes方程.一般曲線坐標(biāo)系中，無量綱化的方程守恒形式為［8］

式中，Q為守恒變矢量;F，G和H為對(duì)流通量;Fv，Gv，Hv為黏性通量;ξ，η 和 ζ為 3 個(gè)貼體坐標(biāo)系方向;t為時(shí)間;Re∞為自由來流雷諾數(shù).

流場求解采用基于結(jié)構(gòu)網(wǎng)格的有限體積法，空間格式采用Roe格式和Minmod限制器，時(shí)間離散方法為LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法，非定常時(shí)間推進(jìn)采用雙時(shí)間步方法.采用SA湍流模型模擬大攻角下可能存在的分離流動(dòng)，使用Open-MP并行技術(shù)提高計(jì)算效率.

1.2 直接阻尼導(dǎo)數(shù)計(jì)算方法

直接阻尼導(dǎo)數(shù)的物理意義是由于飛行器的姿態(tài)角速度的變化引起各處局部迎角的變化，結(jié)果引起質(zhì)心前后升力之差而產(chǎn)生附加的阻尼力矩，該阻尼力矩對(duì)角速度的導(dǎo)數(shù)即為直接阻尼導(dǎo)數(shù)，它本質(zhì)上是準(zhǔn)定常量.

設(shè)體軸坐標(biāo)系x軸朝前，y軸朝上，z軸滿足右手法則.本文采用準(zhǔn)定常方法，假設(shè)飛行器繞z軸以恒定俯仰角速度q拉起，某瞬時(shí)攻角α下有

式中，Cm為飛行器以勻角速率q拉起時(shí)瞬時(shí)攻角α處的俯仰力矩;Cm0為俯仰角速度為零時(shí)攻角α處的俯仰力矩;l為特征長度;u∞為來流速度;

Cm0通過定常計(jì)算獲得;Cm的計(jì)算需要考慮俯仰角速度引起的牽連速度的影響以反映局部迎角的變化，本文在空間格式及壁面邊界條件中加入牽連速度的貢獻(xiàn)，同時(shí)網(wǎng)格保持不動(dòng)，采用準(zhǔn)定常方法進(jìn)行時(shí)間推進(jìn)，經(jīng)過一段暫態(tài)效應(yīng)，最終力和力矩收斂到定常結(jié)果.

需要說明的是，本文采用的方法不同于文獻(xiàn)［5-7］中采用的方法，文獻(xiàn)中Weinacht等人在非慣性坐標(biāo)系下通過在流動(dòng)控制方程中引入源項(xiàng)來考慮牽連速度的影響，而本文方法是在慣性坐標(biāo)系下，求解帶運(yùn)動(dòng)速度物體的瞬態(tài)情況，只需將現(xiàn)有求解強(qiáng)迫運(yùn)動(dòng)的程序稍作修改即可實(shí)現(xiàn).

1.3 洗流時(shí)差導(dǎo)數(shù)計(jì)算方法

洗流時(shí)差導(dǎo)數(shù)本質(zhì)上是非定常量，反映了洗流時(shí)差對(duì)飛行器的阻尼特性.洗流時(shí)差導(dǎo)數(shù)的作用體現(xiàn)在飛行器受陣風(fēng)或直接力作用時(shí)飛行器動(dòng)態(tài)穩(wěn)定性的變化，它是飛行器設(shè)計(jì)時(shí)一個(gè)重要的設(shè)計(jì)參數(shù).

在體軸坐標(biāo)系，給定沉浮運(yùn)動(dòng)形式如下［3］:

在式(4)描述的運(yùn)動(dòng)形式下，可得瞬時(shí)攻角的運(yùn)動(dòng)形式為

式中，α0為起始攻角;αm為攻角振幅;θ為俯仰角;k為減縮頻率.

由式(4)、式(5)確定的沉浮運(yùn)動(dòng)中，確定運(yùn)動(dòng)的狀態(tài)變量只有攻角及其各階導(dǎo)數(shù)，根據(jù)Etkin非定常氣動(dòng)力模型，Cm可寫成

在小振幅運(yùn)動(dòng)情況下，將式(6)泰勒展開并略去高階小量得

俯仰阻尼導(dǎo)數(shù)Cmq+的計(jì)算采用強(qiáng)迫振動(dòng)方法，具體數(shù)值計(jì)算方法參見文獻(xiàn)［9-10］.

2 數(shù)值模擬結(jié)果及分析

本節(jié)采用上述數(shù)值計(jì)算方法，對(duì)HBS標(biāo)模外形及基本帶翼導(dǎo)彈外形進(jìn)行計(jì)算驗(yàn)證和研究分析，并在此基礎(chǔ)上對(duì)Hyflex飛行器進(jìn)行研究.

2.1 HBS 外形研究

HBS為一個(gè)半球柱、帶有兩段擴(kuò)張裙部的高超聲速彈道外形標(biāo)模，如圖1所示.其動(dòng)態(tài)特性有較為精確的試驗(yàn)結(jié)果［11］，常被用來驗(yàn)證計(jì)算結(jié)果的準(zhǔn)確程度.

圖1 彈道外形示意圖Fig.1 Schematic of hyperballistic shape(HBS)

計(jì)算條件為:Ma=6.85，以頭部直徑為參考長度的 Red=0.72 ×106，質(zhì)心位置 Xcg/L=0.72，α0=0°，5°，10°，15°，20°.強(qiáng)迫沉浮運(yùn)動(dòng)及角振動(dòng)的減縮頻率k定義為:k=ωd/2u∞，d為頭部直徑.本節(jié)計(jì)算時(shí)取 k=0.01，振幅取 αm=1°.

圖2給出了計(jì)算得到的俯仰力矩系數(shù)曲線.其中圖2(a)為HBS外形俯仰角速度取q=0，100，300，500，700，1 000(°)/s 時(shí)的俯仰力矩系數(shù)變化曲線，可以看出，俯仰力矩系數(shù)隨俯仰角速度呈完全線性變化，其斜率即為直接阻尼導(dǎo)數(shù)，斜率都為負(fù)，說明直接阻尼導(dǎo)數(shù)為負(fù)，起正阻尼作用.圖2(b)為強(qiáng)迫沉浮運(yùn)動(dòng)時(shí)俯仰力矩系數(shù)的遲滯環(huán)，圖2(c)為強(qiáng)迫角振動(dòng)時(shí)俯仰力矩系數(shù)的遲滯環(huán)，對(duì)比兩圖可以發(fā)現(xiàn)，沉浮運(yùn)動(dòng)的遲滯環(huán)比較瘦，并且0°～15°的遲滯環(huán)為順時(shí)針方向旋轉(zhuǎn)，由于遲滯環(huán)的旋轉(zhuǎn)方向決定動(dòng)導(dǎo)數(shù)的正負(fù)，面積代表動(dòng)導(dǎo)數(shù)的量值大小，這說明HBS的洗流時(shí)差導(dǎo)數(shù)在俯仰阻尼導(dǎo)數(shù)中所占比例較小并且0°～15°的洗流時(shí)差導(dǎo)數(shù)為正值，其洗流特性是負(fù)阻尼的，屬于動(dòng)不穩(wěn)定情況.而相應(yīng)的強(qiáng)迫振動(dòng)的遲滯環(huán)比較飽滿且各攻角下均為逆時(shí)針旋轉(zhuǎn).

圖3為計(jì)算所得俯仰阻尼導(dǎo)數(shù)與試驗(yàn)值［11］對(duì)比，表1為各攻角下直接阻尼導(dǎo)數(shù)、洗流時(shí)差導(dǎo)數(shù)及俯仰阻尼導(dǎo)數(shù)的計(jì)算結(jié)果，其中Cmq+Cmα·為Cmq和直接相加的結(jié)果為通過強(qiáng)迫振動(dòng)直接得到的俯仰阻尼導(dǎo)數(shù)，兩者誤差不超過2%，說明本文發(fā)展的計(jì)算方法是正確的，并且計(jì)算結(jié)果是精確的.

圖2 彈道外形計(jì)算結(jié)果Fig.2 Calculation results of hyperballistic shape(HBS)

圖3 彈道俯仰阻尼導(dǎo)數(shù)計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比Fig.3 Calculation and experiment results of pitch-damping derivatives for hyperballistic shape(HBS)

表1 彈道阻尼導(dǎo)數(shù)計(jì)算結(jié)果Table1 Caculation results of pitch-damping derivatives for hyperballistic shape(HBS)

2.2 基本帶翼導(dǎo)彈外形研究

基本帶翼導(dǎo)彈為一個(gè)尖錐形頭部、圓柱形彈身并帶有4個(gè)矩形尾翼的外形，如圖4所示.

圖4 基本帶翼導(dǎo)彈外形Fig.4 Schematic of Basic Finner

計(jì)算條件為:Ma=1.96，以底部直徑為參考長度的ReD=0.187×106.強(qiáng)迫沉浮運(yùn)動(dòng)及角振動(dòng)的減縮頻率k定義為:k=ωD/2u∞，D為底部直徑.本節(jié)計(jì)算時(shí)取 k=0.01，振幅取 αm=1°.

圖5展示了固定質(zhì)心位置Xcg/D=6.1時(shí)，導(dǎo)彈的直接阻尼導(dǎo)數(shù)Cmq、洗流時(shí)差導(dǎo)數(shù)及俯仰阻尼導(dǎo)數(shù)Cmq+隨攻角的變化規(guī)律(圖中曲線采用樣條曲線擬合而成，下同).計(jì)算攻角為α=0°，2.5°，4°，5°，10°，15°，20°，25°，30°.本文計(jì)算的雷諾數(shù)與風(fēng)洞試驗(yàn)［12］中的高雷諾數(shù)狀態(tài)相同.從圖5中可以看出，本文計(jì)算的俯仰阻尼導(dǎo)數(shù)在大攻角時(shí)(10°以上)與試驗(yàn)值吻合較好，而小攻角時(shí)有一定的偏離.因?yàn)樵囼?yàn)采用尾支撐方式，Uselton等人［13］注意到小攻角時(shí)(6°以下)有較強(qiáng)的支架干擾而導(dǎo)致試驗(yàn)結(jié)果重復(fù)度較差，直到攻角增大到一定程度才會(huì)減弱支架干擾的影響.直接阻尼導(dǎo)數(shù)一直為負(fù)值，且其絕對(duì)值隨攻角增大而近似線性增大;時(shí)差導(dǎo)數(shù)在0°攻角附近為正值，在4°攻角處過零點(diǎn)，然后向負(fù)向增大，在10°攻角附近達(dá)到負(fù)向最大，然后緩慢回落，并一直保持負(fù)值.時(shí)差導(dǎo)數(shù)的非線性變化導(dǎo)致了俯仰阻尼導(dǎo)數(shù)的非線性變化，但總體來看，導(dǎo)彈的縱向動(dòng)穩(wěn)定性隨攻角增大而增強(qiáng)，時(shí)差導(dǎo)數(shù)符號(hào)會(huì)有變化，但量值在組合導(dǎo)數(shù)中所占比例較小(10%以下).

圖5 各阻尼導(dǎo)數(shù)隨攻角變化曲線Fig.5 Variations of pitch-damping derivatives with angles of attack

圖6 各阻尼導(dǎo)數(shù)隨質(zhì)心位置變化曲線Fig.6 Variations of pitch-damping derivatives with center of gravity(CG)location

圖6為0°攻角時(shí)計(jì)算的導(dǎo)彈各阻尼導(dǎo)數(shù)隨質(zhì)心位置的變化曲線圖.計(jì)算的質(zhì)心位置為Xcg/D=5，6.1，7.從圖6中可以看出，時(shí)差導(dǎo)數(shù)為正值，在質(zhì)心Xcg/D=5處幾乎為零且隨質(zhì)心位置后移而線性增加;直接阻尼導(dǎo)數(shù)及組合導(dǎo)數(shù)則呈拋物分布，這與文獻(xiàn)［14］給出的質(zhì)心平移關(guān)系式(8)～式(10)吻合很好.

式中scg=(－xcg)/D，坐標(biāo)定義在體軸系.上式只適用于軸對(duì)稱物體.

雖然沒有在此列出，但計(jì)算中CNα，CNq，均為正值，因此時(shí)差導(dǎo)數(shù)數(shù)值隨質(zhì)心后移增加，即朝動(dòng)不穩(wěn)定性增大的方向發(fā)展.而直接阻尼導(dǎo)數(shù)與質(zhì)心偏移量是復(fù)雜的二次函數(shù)關(guān)系.從數(shù)值上來看，存在一個(gè)質(zhì)心位置使俯仰阻尼導(dǎo)數(shù)值最大，即動(dòng)穩(wěn)定性最差.

2.3 Hyflex外形研究

Hyflex(Hypersonic Flight Experiment)是日本HOPE-X計(jì)劃中有關(guān)大氣層再入項(xiàng)目的帶翼升力體外形高超聲速飛行器，用于驗(yàn)證制導(dǎo)和控制以及熱防護(hù)材料和結(jié)構(gòu)等技術(shù).它于1996年2月進(jìn)行飛行試驗(yàn)，本文選取其飛行末端彈道點(diǎn)(馬赫2.0，高度30 km)進(jìn)行研究［15］，該彈道點(diǎn)處于飛行器減速傘開傘前，其動(dòng)態(tài)特性對(duì)于減速傘的開啟有著至關(guān)重要的影響.飛行器外形三視圖及網(wǎng)格圖見圖7.

圖7 Hyflex外形及網(wǎng)格圖Fig.7 Schematic and mesh of Hyflex

計(jì)算條件為:Ma=2.0，H=30 km，Sref=4.27 m2，Lref=4 m.計(jì)算攻角為:α =0°，5°，10°，15°，20°.減縮頻率定義為:k= ωLref/2u∞，本節(jié)計(jì)算時(shí)取 k=0.01，振幅取 αm=1°.

表2給出了各攻角下的直接阻尼導(dǎo)數(shù)、洗流時(shí)差導(dǎo)數(shù)及俯仰阻尼導(dǎo)數(shù)計(jì)算結(jié)果.可以看出，直接計(jì)算得到的俯仰阻尼導(dǎo)數(shù)與兩個(gè)分量相加得到的俯仰阻尼導(dǎo)數(shù)Cmq+最大相對(duì)誤差為6%，雖然沒有試驗(yàn)數(shù)據(jù)作為對(duì)比，但通過兩種不同方法得到的俯仰阻尼導(dǎo)數(shù)相差不大，使得本文發(fā)展的阻尼導(dǎo)數(shù)計(jì)算方法的正確性得到了很好的自證.

表2 Hyflex阻尼導(dǎo)數(shù)計(jì)算結(jié)果Table2 Calculation results of pitch-damping derivatives for Hyflex

另外，從表2中可以看出，直接阻尼導(dǎo)數(shù)為負(fù)值，且絕對(duì)值隨攻角增加而單調(diào)增加;洗流時(shí)差導(dǎo)數(shù)也全為負(fù)值，且絕對(duì)值先增加后減小，隨攻角呈非線性變化，在－15°達(dá)到最大，在俯仰阻尼導(dǎo)數(shù)中所占比例不再是小量，反而起主導(dǎo)作用(最大比重達(dá)78%).與2.1節(jié)和2.2節(jié)研究的外形相比，Hyflex是帶翼升力體外形飛行器，尾部翼面的洗流時(shí)差效果更顯著.

3 結(jié)論

本文研究直接求解俯仰阻尼導(dǎo)數(shù)分量的方法，通過對(duì)HBS和基本帶翼導(dǎo)彈兩個(gè)標(biāo)模外形及Hyflex飛行器進(jìn)行數(shù)值模擬研究，結(jié)果表明:

1)采用的方法能夠準(zhǔn)確預(yù)測飛行器外形的俯仰阻尼導(dǎo)數(shù)的各個(gè)分量，即使對(duì)于大攻角狀態(tài)也具有較好的預(yù)測精度，具備一定的工程實(shí)用價(jià)值;

2)對(duì)于彈丸類彈道物體，其在超聲速及高超聲速狀態(tài)下，洗流時(shí)差導(dǎo)數(shù)在俯仰阻尼導(dǎo)數(shù)中所占比例較小，但其符號(hào)可能大于零，起動(dòng)不穩(wěn)定作用，并且數(shù)值隨質(zhì)心后移而增大;對(duì)于帶翼飛行器，超聲速狀態(tài)下，洗流時(shí)差導(dǎo)數(shù)在俯仰阻尼導(dǎo)數(shù)中所占比例較大.

直接阻尼導(dǎo)數(shù)不論從物理意義上來看還是數(shù)值計(jì)算都是負(fù)值，總起動(dòng)穩(wěn)定作用.

References)

［1］ 李周復(fù).風(fēng)洞特種試驗(yàn)技術(shù)［M］.北京:航空工業(yè)出版社，2010:210-250.Li Z F.Special wind tunnel testing technology［M］.Beijing:Aviation Industry Press，2010:210-250(in Chinese).

［2］ Ren Y X.Evaluation of the stability derivatives using the sensitivity equations［J］.AIAA Journal，2008，46(4):912-917.

［3］ 劉偉，楊小亮，趙云飛.高超聲速飛行器加速度導(dǎo)數(shù)數(shù)值模擬［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2010，28(4):426-429.Liu W，Yang X L，Zhao Y F.Numerical simulation of acceleration derivative of hypersonic aircraft［J］.Acta Aerodynamica Sinica，2010，28(4):426-429(in Chinese).

［4］ 孫海生.飛行器大攻角升沉平移加速度導(dǎo)數(shù)測量技術(shù)［J］.流體力學(xué)實(shí)驗(yàn)與測量，2001，15(4):15-19.Sun H S.A measurement technique for derivatives of aircraft due to acceleration in heave and sideslip at high angle of attack［J］.Experiments and Measurements in Fluid Mechanics，2001，15(4):15-19(in Chinese).

［5］ Weinacht P.Navier-Stokes predictions of the individual components of the pitch-damping coefficient sum，ARL-TR-3169［R］.Adelphi:Army Research Laboratory，2004.

［6］ Weinacht P.Projectile performance，stability and free flight motion prediction using computational fluid dynamics［J］.Journal of Spacecraft and Rockets，2004，41(2):257-263.

［7］ Qin N，Ludlow D K，Shaw S T，et al.Calculation of pitch damping coefficients for projectiles，AIAA-1997-0405［R］.Reston:AIAA，1997.

［8］ 閻超.計(jì)算流體力學(xué)方法及應(yīng)用［M］.北京:北京航空航天大學(xué)出版社，2006:18-25.Yan C.The methodology and application of computational fluid dynamics［M］.Beijing:Beihang University Press，2006:18-25(in Chinese).

［9］ McGowan G Z，Kurzen M J，Nance R P，et al.High fidelity approaches for pitch damping prediction at high angles of attack，AIAA-2012-2903［R］.Reston:AIAA，2012.

［10］ Hashimoto A，Hashizume M，Sunada S，et al.Unsteady analysis of aerodynamic derivatives on standard dynamics model，AIAA-2013-0343［R］.Reston:AIAA，2013.

［11］ East R A，Hutt G R.Comparison of predictions and experimental data for hypersonic pitching motion stability［J］.Journal of Spacecraft and Rockets，1988，25(3):225-233.

［12］ Uselton B L，Uselton J C.Test mechanism for measuring pitch damping derivatives of missile configurations at high angles of attack，AEDC-TR-75-43［R］.Tennessee:AEDC，1975.

［13］ Uselton B L，Jenke L M.Experimental missile pitch-and rolldamping characteristics at large angles of attack［J］.Journal of Spacecraft and Rockets，1977，14(4):241-247.

［14］ Murphy C H.Free flight motion of symmetric missiles，NO.1216［R］.Aberdeen:Army Ballistic Research Lab，1963.

［15］ Shirouzu M，Yamamoto M.Overview of the hyflex project，AIAA-1996-4524［R］.Reston:AIAA，1996.