TBCC飛行器發(fā)動(dòng)機(jī)尺寸選型及爬升策略設(shè)計(jì)

李 樂(lè), 姜光泰, 褚顯應(yīng), 王 健

(北京機(jī)電工程研究所，北京 100074)

0 引 言

高度在20～100 km的空間被稱為臨近空間，伴隨著動(dòng)力技術(shù)和控制技術(shù)的發(fā)展，臨近空間高速飛行器成為大國(guó)爭(zhēng)相發(fā)展的熱點(diǎn)裝備。該類飛行器要實(shí)現(xiàn)遠(yuǎn)程飛行，需盡量降低爬升段油耗，提高巡航段升阻比和比沖。為此，氣動(dòng)外形與動(dòng)力系統(tǒng)匹配設(shè)計(jì)以及爬升策略選擇是關(guān)鍵。一方面，臨近空間高速飛行器動(dòng)力系統(tǒng)通常選用渦輪沖壓組合(Turbine Based Combined Cycle, TBCC)發(fā)動(dòng)機(jī)[1]，內(nèi)外流耦合作用明顯，研究發(fā)動(dòng)機(jī)尺寸選型對(duì)巡航航程的影響非常重要。另一方面，臨近空間高速飛行器爬升段高度變化大，加速性能要求高，采用TBCC發(fā)動(dòng)機(jī)易出現(xiàn)高空推力不足問(wèn)題，可考慮采用固體火箭助推器在適當(dāng)時(shí)機(jī)輔助加速爬升，為巡航段剩余更多燃油。

近年來(lái)，國(guó)內(nèi)利用偽譜法開(kāi)展臨近空間飛行器航跡優(yōu)化的研究有很多[2-4]?；趥巫V法研究TBCC發(fā)動(dòng)機(jī)尺寸與飛行器適配性和助推器與飛行器的適配性問(wèn)題的研究較少。

本文擬通過(guò)hp自適應(yīng)Radau偽譜法，以消耗燃料質(zhì)量最小為目標(biāo)優(yōu)化航跡，研究TBCC發(fā)動(dòng)機(jī)尺寸選型和固體火箭助推器質(zhì)量、推力、點(diǎn)火時(shí)刻、工作時(shí)間等因素對(duì)臨近空間高速飛行器爬升航跡和巡航航程的影響。以此為基礎(chǔ)，提出一種TBCC發(fā)動(dòng)機(jī)尺寸選型的依據(jù)和固體火箭助推器輔助加速的爬升策略。為未來(lái)臨近空間高速飛行器發(fā)動(dòng)機(jī)尺寸選型和選用固體火箭助推器提供參考。

1 臨近空間高速飛行器模型

本文設(shè)計(jì)了一型臨近空間高速飛行器。其設(shè)計(jì)條件為：巡航高度25 km，巡航馬赫數(shù)3.5。動(dòng)力系統(tǒng)采用串聯(lián)式渦輪沖壓組合(TBCC)發(fā)動(dòng)機(jī)。經(jīng)過(guò)設(shè)計(jì)計(jì)算，飛行器氣動(dòng)外形如圖1所示。飛行器基本外形和質(zhì)量參數(shù)如表1所示。

2 航跡優(yōu)化數(shù)學(xué)模型

研究臨近空間高速飛行器航跡優(yōu)化問(wèn)題，需要建立動(dòng)力學(xué)模型、大氣模型、發(fā)動(dòng)機(jī)模型和氣動(dòng)力模型。

2.1 動(dòng)力學(xué)模型

忽略地球自轉(zhuǎn)影響，考慮地球曲率，只研究飛行器縱向平面內(nèi)運(yùn)動(dòng)，建立的動(dòng)力學(xué)模型如式(1)。

(1)

式中r=Re+h，h為飛行高度，Re為地球半徑，v為飛行速度，m為飛行器質(zhì)量，θ為航跡傾角，T為發(fā)動(dòng)機(jī)推力，ms為發(fā)動(dòng)機(jī)耗油率，k為偽控制量，D為阻力，L為升力，α為迎角，μ為地球引力常數(shù)。

2.2 發(fā)動(dòng)機(jī)模型

本文設(shè)計(jì)的臨近空間高速飛行器采用串聯(lián)軸對(duì)稱渦輪沖壓組合發(fā)動(dòng)機(jī)(TBCC)提供推力，組合發(fā)動(dòng)機(jī)原型的推力和耗油率數(shù)據(jù)參考J-58發(fā)動(dòng)機(jī)[5-7]得到，與高度、來(lái)流馬赫數(shù)和迎角相關(guān)。0度迎角下，發(fā)動(dòng)機(jī)推力和耗油率隨飛行高度和速度的變化規(guī)律如圖2所示。來(lái)流Ma3.5時(shí)，不同高度發(fā)動(dòng)機(jī)推力和耗油率隨迎角變化規(guī)律如圖3所示。

2.3 大氣模型

大氣模型采用1976年美國(guó)標(biāo)準(zhǔn)大氣模型，依據(jù)高度插值得到聲速和空氣密度。

2.4 氣動(dòng)力模型

飛行器飛行過(guò)程中受到的升力、阻力計(jì)算公式如式(2)和式(3)。

(2)

(3)

式中:S為參考面積，即為機(jī)翼投影面積，ρ為大氣密度，v為飛行速度，CL和CD分別為升力系數(shù)和阻力系數(shù)，是迎角α和馬赫數(shù)M的函數(shù)，其利用Digital DATCOM[8]計(jì)算得到，如圖4所示。

2.5 約束條件

(1)爬升段初末狀態(tài)約束

表2 爬升段初末狀態(tài)列表Table 2 The initial and final states of ascent trajectory

(2)爬升段路徑約束

在吸氣式組合動(dòng)力飛行器爬升過(guò)程中，約束航跡傾角以保證足夠的推阻余量；約束迎角以保證發(fā)動(dòng)機(jī)進(jìn)氣量和進(jìn)氣道啟動(dòng)；迎角變化率約束主要考慮飛行器法向過(guò)載約束；飛行器發(fā)動(dòng)機(jī)可靠工作、翼面載荷和舵面鉸鏈力矩要求對(duì)最大動(dòng)壓予以約束；頭部駐點(diǎn)熱流密度也需要約束以滿足飛行器熱防護(hù)系統(tǒng)要求。

航跡傾角約束：-40°≤θ≤40°,

迎角約束：-5°≤α≤10°,

偽控制量約束：0≤k≤1,

法向過(guò)載約束：

(4)

頭部駐點(diǎn)熱流約束：

(5)

最大動(dòng)壓約束：

(6)

分析仿真結(jié)果得到主要約束條件是航跡傾角約束、迎角變化率約束和偽控制量約束。

(3)巡航段初末狀態(tài)及路徑約束

巡航段初始條件即是爬升段末端條件，巡航段的末端條件設(shè)為飛行器剩余5%油量(用于降落)。巡航段路徑約束條件為：

h=25 km，Ma=3.5，θ=0°。

3 航跡優(yōu)化問(wèn)題求解

針對(duì)最優(yōu)控制問(wèn)題，本文采用Radau偽譜法求解。該方法在一組Legendre-Gauss-Radau(LGR)點(diǎn)上構(gòu)造全局Lagrange插值多項(xiàng)式，近似狀態(tài)變量和控制變量，用插值多項(xiàng)式導(dǎo)數(shù)近似微分方程中狀態(tài)變量的導(dǎo)數(shù)，將最優(yōu)控制問(wèn)題轉(zhuǎn)化為非線性規(guī)劃問(wèn)題[9-10]，再通過(guò)序列二次規(guī)劃(Sequential Quadratic Programming,SQP)算法[11]對(duì)其進(jìn)行求解。

基于hp自適應(yīng)Radau偽譜法汲取了全局偽譜法和網(wǎng)格細(xì)分法的優(yōu)點(diǎn)，采用雙層優(yōu)化策略，具有更合理的LGR點(diǎn)分布，在保證優(yōu)化精度的前提下有效提高了算法的優(yōu)化效率。其基本思想是檢驗(yàn)各區(qū)間內(nèi)的相對(duì)誤差，若超出規(guī)定值，則采用細(xì)化區(qū)間或增加區(qū)間內(nèi)LGR點(diǎn)個(gè)數(shù)來(lái)提高優(yōu)化精度；而相對(duì)誤差小于規(guī)定值的區(qū)間則保持不變[12-13]。這樣只需在個(gè)別區(qū)間細(xì)化網(wǎng)格，增加采樣密度，從而在滿足精度要求的情況下，有效提高了優(yōu)化速度?；趆p自適應(yīng)Radau偽譜法的流程框圖如圖5所示。

4 TBCC發(fā)動(dòng)機(jī)尺寸選型

TBCC發(fā)動(dòng)機(jī)作為動(dòng)力系統(tǒng)對(duì)飛行器的飛行性能有重要影響。本節(jié)討論TBCC發(fā)動(dòng)機(jī)尺寸選型對(duì)臨近空間高速飛行器爬升航跡和巡航航程的影響。

首先，本文參考J-58發(fā)動(dòng)機(jī)建立的基準(zhǔn)TBCC發(fā)動(dòng)機(jī)主要參數(shù)如表3所示。

表3 基準(zhǔn)TBCC發(fā)動(dòng)機(jī)主要參數(shù)列表Table 3 Parameters of basic TBCC engine

在方案論證階段，發(fā)動(dòng)機(jī)性能參數(shù)主要依靠估算得到。由于可供參考的TBCC發(fā)動(dòng)機(jī)型號(hào)有限，因此本文通過(guò)應(yīng)用經(jīng)驗(yàn)和參考相關(guān)型號(hào)參數(shù)，得到如下近似擬合公式：

T=T0·K2

(7)

(8)

發(fā)動(dòng)機(jī)質(zhì)量遵從如下變化規(guī)律

m=m0·K1.35

(9)

不同尺寸發(fā)動(dòng)機(jī)主要通過(guò)三種途徑影響爬升航跡。首先，不同尺寸發(fā)動(dòng)機(jī)的推力不同；其次，不同尺寸發(fā)動(dòng)機(jī)的質(zhì)量不同，飛行器起飛總重不變條件下，發(fā)動(dòng)機(jī)質(zhì)量變化會(huì)影響飛行器載油質(zhì)量，如表4所示；最后，發(fā)動(dòng)機(jī)的直徑不同，引起的氣動(dòng)力也不同。為考察爬升段內(nèi)發(fā)動(dòng)機(jī)尺寸對(duì)飛行器氣動(dòng)力的影響，本文進(jìn)行了大量仿真計(jì)算。這里選取兩個(gè)典型狀態(tài)展示其影響，分別是模態(tài)轉(zhuǎn)換起始時(shí)(Ma2)和迎角為0.5°時(shí)的情況，討論不同尺寸發(fā)動(dòng)機(jī)對(duì)飛行器氣動(dòng)力影響隨迎角和馬赫數(shù)的變化規(guī)律，結(jié)果如圖6～7所示。

觀察圖6～7，一方面，可以發(fā)現(xiàn)發(fā)動(dòng)機(jī)尺寸對(duì)飛行器升力系數(shù)的影響較小，對(duì)阻力系數(shù)影響較大；另一方面，發(fā)動(dòng)機(jī)尺寸對(duì)飛行器氣動(dòng)力的影響隨馬赫數(shù)變化較為敏感，且對(duì)爬升過(guò)程中的阻力峰值影響較大。在不同迎角下，發(fā)動(dòng)機(jī)尺寸增大導(dǎo)致飛行器的升力和阻力同時(shí)增加。但發(fā)動(dòng)機(jī)尺寸對(duì)飛行器氣動(dòng)力的影響隨馬赫數(shù)的變化則不規(guī)律。

表4 不同尺寸發(fā)動(dòng)機(jī)的直徑和質(zhì)量列表Table 4 Diameter and mass of different size of TBCCengines

本文以所設(shè)計(jì)的臨近空間高速飛行器為原型，基于第2節(jié)建立的航跡優(yōu)化數(shù)學(xué)模型，代入不同尺寸TBCC發(fā)動(dòng)機(jī)，依據(jù)優(yōu)化得到的控制變量變化規(guī)律，仿真得到爬升航跡參數(shù)變化曲線如圖8～11所示，其中K為發(fā)動(dòng)機(jī)直徑比。

通過(guò)圖8～9可以看出，當(dāng)采用較小尺寸發(fā)動(dòng)機(jī)時(shí)，飛行器高度爬升有遲滯，甚至有俯沖加速的趨勢(shì)；當(dāng)采用較大尺寸發(fā)動(dòng)機(jī)時(shí)，飛行器可以直接爬升到巡航高度，但速度爬升末段存在遲滯，甚至有減速爬升的趨勢(shì)。其原因是為減小爬升段消耗燃料，如圖10所示飛行器在相應(yīng)階段燃料消耗率明顯降低。圖11顯示配備大尺寸發(fā)動(dòng)機(jī)的飛行器在轉(zhuǎn)平階段短時(shí)間內(nèi)有負(fù)迎角，且發(fā)動(dòng)機(jī)尺寸越大負(fù)迎角度數(shù)和持續(xù)時(shí)間也越長(zhǎng)。隨著發(fā)動(dòng)機(jī)尺寸變大，飛行器爬升消耗燃料質(zhì)量減小。但飛行器爬升航跡趨于一致，反映在爬升消耗燃料質(zhì)量趨同，如圖12所示。

圖13比較了配備相同尺寸發(fā)動(dòng)機(jī)的飛行器，分別以爬升消耗燃料質(zhì)量最小和爬升所用時(shí)間最小為目標(biāo)優(yōu)化得到的結(jié)果，兩種目標(biāo)優(yōu)化得到的航跡，其所用時(shí)間和爬升消耗燃料質(zhì)量差值隨發(fā)動(dòng)機(jī)尺寸增大而減小。這說(shuō)明較大尺寸發(fā)動(dòng)機(jī)可以近似同時(shí)滿足爬升最快和最省燃料的要求。

接下來(lái)討論不同尺寸發(fā)動(dòng)機(jī)對(duì)巡航航程的影響。首先，相同尺寸發(fā)動(dòng)機(jī)和巡航飛行條件一定情況下，爬升消耗燃料質(zhì)量越小，意味著有更多燃料可供巡航飛行，因而爬升消耗燃料質(zhì)量越小，飛行器巡航航程越大。以此為基礎(chǔ)開(kāi)展飛行器巡航段仿真。

本文將飛行器需用推力與此環(huán)境下飛行器可用推力之比定義為發(fā)動(dòng)機(jī)推力利用率。仿真結(jié)果圖14顯示，隨著發(fā)動(dòng)機(jī)尺寸增大，巡航段推力利用率逐步降低，即過(guò)大尺寸發(fā)動(dòng)機(jī)會(huì)導(dǎo)致巡航推力的較大浪費(fèi)，且會(huì)導(dǎo)致較大的飛行阻力和減小飛行器所載燃料質(zhì)量，不利于提升航程。

由圖15可看出，巡航起始點(diǎn)升阻比隨發(fā)動(dòng)機(jī)尺寸增大先增大后減小，巡航航程與巡航起始點(diǎn)升阻比的變化趨勢(shì)基本一致，說(shuō)明影響巡航航程的主要因素仍是巡航升阻比。綜上所述，TBCC飛行器的動(dòng)力系統(tǒng)尺寸匹配問(wèn)題應(yīng)綜合考慮爬升和巡航段要求，雖然大尺寸發(fā)動(dòng)機(jī)在爬升階段可縮短時(shí)間和節(jié)省燃料，但會(huì)造成巡航段推力利用率低。對(duì)于遠(yuǎn)程飛行器，發(fā)動(dòng)機(jī)尺寸選型應(yīng)主要考慮其對(duì)飛行器巡航升阻比的影響。

吸氣式組合動(dòng)力飛行器的爬升方法常采用等動(dòng)壓爬升和等熱流密度爬升方法。本文基于1.291倍尺寸發(fā)動(dòng)機(jī)，通過(guò)將動(dòng)壓關(guān)于高度積分后得到爬升平均動(dòng)壓qaver，如式(10)。

(10)

式中qe為飛行過(guò)程中來(lái)流動(dòng)壓，h0為起始點(diǎn)高度，hf為巡航點(diǎn)高度。基于參考文獻(xiàn)[14]方法生成qaver的等動(dòng)壓爬升航跡如圖16所示。進(jìn)一步研究發(fā)動(dòng)機(jī)尺寸選型對(duì)等動(dòng)壓爬升航跡與優(yōu)化航跡偏差的影響，結(jié)果如圖17所示。

通過(guò)觀察圖16和圖17可以發(fā)現(xiàn)存在最優(yōu)尺寸發(fā)動(dòng)機(jī)，其以爬升消耗燃料質(zhì)量最少為目標(biāo)的優(yōu)化航跡可以部分近似表示為等動(dòng)壓爬升航跡，這為解析表示最優(yōu)航跡創(chuàng)造了條件。

5 爬升段助推器選型研究

通過(guò)分析上節(jié)結(jié)果可知，TBCC發(fā)動(dòng)機(jī)尺寸選型在滿足最優(yōu)巡航升阻比條件下存在巡航狀態(tài)推力余量較大的問(wèn)題。為獲取最大航程，需要盡量提高巡航起始點(diǎn)推力利用率。這是一個(gè)多點(diǎn)設(shè)計(jì)問(wèn)題，其實(shí)質(zhì)是受限于渦扇和沖壓發(fā)動(dòng)機(jī)工作特性，飛行器與TBCC發(fā)動(dòng)機(jī)在起飛、爬升和高空巡航階段難以完全匹配。為此，本文提出一種折衷方案，即考慮在爬升段采用固體火箭助推器輔助加速，以提升飛行器爬升后剩余質(zhì)量，提高巡航段起始點(diǎn)推力利用率。

經(jīng)過(guò)分析，固體火箭助推器潛在使用階段包括：起飛階段，爬升階段和轉(zhuǎn)平階段。在起飛階段使用固體火箭助推器可快速脫離稠密大氣，減小飛行阻力；在爬升階段使用固體助推器有利于飛行器快速通過(guò)跨聲速阻力峰值；在轉(zhuǎn)平階段使用固體助推器有助于緩解TBCC發(fā)動(dòng)機(jī)高空爬升動(dòng)力不足的問(wèn)題。

針對(duì)如上助推器的三種使用方案，本文以爬升段最省燃料為目標(biāo)，選取助推器質(zhì)量、裝藥質(zhì)量流量和點(diǎn)火時(shí)刻作為優(yōu)化變量。動(dòng)力系統(tǒng)采用1.155倍尺寸TBCC發(fā)動(dòng)機(jī)。采用內(nèi)外兩層優(yōu)化循環(huán)，內(nèi)層循環(huán)采用Radau偽譜法對(duì)爬升航跡優(yōu)化，相關(guān)狀態(tài)變量、控制變量約束和路徑約束同之前仿真條件。外層循環(huán)依據(jù)內(nèi)層優(yōu)化得到的爬升消耗燃料質(zhì)量，采用梯度算法優(yōu)化得到優(yōu)化變量取值。具體優(yōu)化框圖如圖18所示。

采用兩枚助推器，在飛行器上掛載示意圖如圖19所示。飛行器起飛總重保持不變，爬升消耗燃料質(zhì)量包括助推器質(zhì)量。并且考慮加裝助推器對(duì)飛行器氣動(dòng)性能的影響。

固體火箭助推器的其余參數(shù)由式(11-13)推導(dǎo)得到：

(11)

(12)

(13)

表5 固體火箭發(fā)動(dòng)機(jī)不同高度性能列表Table 5 Performance of solid rocket booster in variable altitudes

助推器裝藥燃燒速率取為40 mm/s，裝藥密度取為1.78 g/cm3。單枚助推器總質(zhì)量為2.5噸情況下，不同裝藥質(zhì)量流量對(duì)應(yīng)的助推器外形參數(shù)如表6所示。

表6 助推器不同裝藥質(zhì)量流量對(duì)應(yīng)的助推器外形Table 6 Appearance feature of solid rocket booster of different fuel flow rate

助推器燃后拋離，因而在助推段與下一段的質(zhì)量連接條件如式(14)所示。

m2f-m30=0.12mz

(14)

本文設(shè)計(jì)飛行器在Ma1.5左右阻力較大，圖21顯示優(yōu)化得到的助推器工作馬赫數(shù)主要覆蓋這一區(qū)間，因而助推器的主要作用是克服爬升過(guò)程中的阻力峰值。圖22顯示配備越大推力的固體火箭助推器對(duì)巡航航程的提升作用越大。

接下來(lái)討論三種使用階段固體火箭助推器的作用效果，結(jié)果如表7所示。

通過(guò)觀察表7可得出如下結(jié)論：在爬升不同階段均傾向于使用大推力固體火箭助推器。在爬升克服阻力峰值時(shí)使用助推器節(jié)省燃料效果最明顯，在起飛階段使用助推器的效果稍差，在轉(zhuǎn)平階段使用助推器會(huì)增加爬升所需燃料。其原因是在相同裝藥質(zhì)量下，大推力固體火箭助推器能在較短時(shí)間內(nèi)將飛行器助推至較高高度，有利于減小氣動(dòng)阻力，節(jié)省燃料；大推力助推器可助推飛行器快速通過(guò)阻力峰值，減小阻力沖量；在轉(zhuǎn)平段使用固體助推器由于點(diǎn)火時(shí)刻較晚，助推器帶來(lái)的阻力增量作用于全部爬升段航跡，阻力增量做功抵消了助推器做功，因而在轉(zhuǎn)平段使用助推器反而會(huì)增加爬升所需燃料。

表7 不同飛行階段優(yōu)化得到的助推器參數(shù)及其影響Table 7 Optimized performance of solid rocket booster in different ascent phases

基于以上分析結(jié)果，本文提出一種采用固體火箭助推器輔助加速的爬升策略，首先采用固體火箭助推器輔助加速可以提升飛行器巡航段推力利用率。在對(duì)飛行器氣動(dòng)性能影響較小的情況下，選用大質(zhì)量、大推力的固體火箭助推器，在飛行器阻力峰值附近點(diǎn)火加速的效果最好。

6 結(jié) 論

本文設(shè)計(jì)了一型臨近空間高速飛行器，基于hp自適應(yīng)Radau偽譜法開(kāi)展爬升段航跡優(yōu)化，分析了TBCC動(dòng)力系統(tǒng)尺寸選型和選用固體火箭助推器輔助加速策略對(duì)飛行器爬升段航跡和巡航航程的影響?；诜治鼋Y(jié)果，首先提出TBCC飛行器動(dòng)力系統(tǒng)尺寸選型的依據(jù)，應(yīng)綜合考慮動(dòng)力需求和對(duì)巡航配平升阻比的影響，對(duì)于遠(yuǎn)程飛行器應(yīng)重點(diǎn)考慮尺寸選型對(duì)配平升阻比的影響；其次證明了爬升段采用固體火箭助推器加速爬升策略可以提高飛行器巡航段推力利用率，增大航程。在飛行器阻力峰值附近采用大質(zhì)量、大推力的固體火箭助推器輔助加速的效果最好。

[1] 王巍巍, 李丹,曾軍.國(guó)外水平起降臨近空間高速飛機(jī)動(dòng)力的發(fā)展[J]. 燃?xì)鉁u輪試驗(yàn)與研究, 2014, 27(1):57-62. [Wang Wei-wei, Li Dan, Zeng Jun. Development of propulsion system for near-space HTHL high-speed vehicles abroad [J].Gas Turbine Experiment and Research, 2014, 27(1):57-62.]

[2] 龔春林,韓璐,谷良賢.適應(yīng)于RBCC運(yùn)載器的軌跡優(yōu)化建模研究[J]. 宇航學(xué)報(bào)，2013，34(12):1592-1598.[Gong Chun-lin, Han Lu, Gu Liang-xian. Research on modeling of trajectory optimization for RBCC-powered RLV [J]. Journal of Astronautics, 2013, 34(12):1592-1598]

[3] 明超,孫瑞勝,白弘陽(yáng),等.吸氣式超聲速導(dǎo)彈爬升段多約束軌跡優(yōu)化[J].宇航學(xué)報(bào),2016,37(9):1063-1070.[Ming Chao, Sun Rui-sheng, Bai Hong-yang, et al. Climb trajectory optimization with multiple constraints for air-breathing supersonic missile[J]. Journal of Astronautics,2016,37(9):1063-1070.]

[4] 劉瑞帆,于云峰,閆斌斌.基于改進(jìn)hp自適應(yīng)偽譜法的高超聲速飛行器上升段軌跡規(guī)[J], 西北工業(yè)大學(xué)學(xué)報(bào),2016,34(5):790-797.[Liu Rui-fan, Yu Yun-feng, Yan Bin-bin. Ascent phase trajectory optimization for hypersonic vehicle based on hp-adaptive pseudospectral method[J].Journal of Northwes-tern Polytechnical University, 2016,34(5):790-797]

[5] James D. H. Exhaust emission calibration of two J-58 afterburning turbojet engines at simulated high-altitude, supersonic flight conditions[R]. Washington D. C.: National Aeronautics and Space Administration, February 1976.

[6] James D. H. Measurement of exhaust emissions from two J-58 engines at simulated supersonic cruise flight conditions[R]. Washington D. C.:National Aeronautics and Space Administration, March 1976.

[7] David H C. F-12 series aircraft propulsion system performance and development[C], AIAA 5thAircraft Design, Flight Test and Operations Meeting, ST. Louis, USA, August 6-8, 1973.

[8] Mc Donnell Douglas Astronautics Company. The USAF stability and control digital DATCOM, volume I, user’s manual[R]. Missouri: McDonnell Douglas Astronautics Company, April 1979.

[9] Divya G, Michael P, William W. H, et al. A unified framework for the numerical solution of optimal control problems using pseudospectral methods[J], Automatica, 2010,46(11): 1843-1851.

[10] Han P, Shan J Y. Re-entry trajectory optimization using radau pseudospectral method [J]. Control Theory & Applications, 2013, 30(8): 1027-1032.

[11] Philip E.G, Walter M, Michael A.S. SNOPT: an SQP algorithm for large-scale constrained optimization[J], SIAM Journal on Optimization, 2002, 12(4):979-1006.

[12] Michael A P, Anil V R.Exploiting sparsity in direct collocation pseudospectral methods for solving optimal control problems[J], Journal of Spacecraft and Rockets, 2012, 49(2): 364-377.

[13] 雷虎民, 劉滔，李炯,等. 網(wǎng)格縮減的自適應(yīng)hp偽譜法 [J]. 控制理論與應(yīng)用 2016, 33(8): 1061-1067.[Lei Hu-min, Liu Tao, Li Jiong, et al. Adaptive mesh refinement of hp pseudospectral method using mesh size reduction [J].Control Theory & Applications, 2016, 33(8)：1061-1067.]

[14] 賈曉娟, 閆曉東. 吸氣式組合動(dòng)力飛行器爬升軌跡設(shè)計(jì)方法研究 [J]. 西北工業(yè)大學(xué)學(xué)報(bào), 2015, 33(1): 104-109. [ Jia Xiao-juan, Yan Xiao-dong. Ascent trajectory design method for air-breathing powered propulsion system [J]. Journal of Northwestern Polytechnical University, 2015, 33(1)：104-109.]