曹 粟，蔣 鋒，李 易

(1.西北工業(yè)大學(xué) 航天學(xué)院，陜西 西安 710072； 2.陜西省空天飛行器設(shè)計(jì)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710072； 3.上海航天技術(shù)研究院，上海 201109)

0 引言

與傳統(tǒng)飛行器相比，高超聲速滑翔飛行器突破了常規(guī)的彈道再入飛行模式，在再入大氣層后能以極高的速度飛行，具有很強(qiáng)的突防能力。因此，未來臨近空間防御問題也將愈加突出。通常情況下，為有效攔截臨近空間目標(biāo)，攔截器需獲得更高的飛行速度，而再入式飛行是獲得較高飛行速度的有效方式。然而，再入式飛行器通常面臨速域?qū)?、航程遠(yuǎn)和飛行環(huán)境惡劣的任務(wù)特點(diǎn)。在設(shè)計(jì)時(shí)，要考慮攔截過程中的總熱流量、飛行時(shí)間、航程等問題[1]。因這些指標(biāo)一般具有沖突性，故需采用多目標(biāo)優(yōu)化方法進(jìn)行設(shè)計(jì)和研究。

求解多目標(biāo)問題時(shí)，通常將原問題分解為多個(gè)單目標(biāo)子問題。文獻(xiàn)[2]基于邊界交叉法,對(duì)再入軌跡優(yōu)化問題進(jìn)行分解和優(yōu)化。文獻(xiàn)[3]提出一種三維流場(chǎng)的乘波體外形快速生成方法，通過調(diào)整其前緣曲線，實(shí)現(xiàn)飛行器升阻比和容積率的兩目標(biāo)問題求解。本文主要采用基于分解的多目標(biāo)進(jìn)化算法(MOEA/D)和標(biāo)量分解方法，將相互耦合的多目標(biāo)優(yōu)化問題轉(zhuǎn)化為單目標(biāo)問題，從而實(shí)現(xiàn)優(yōu)化[4]。相較于傳統(tǒng)的NSGA-Ⅱ等方法，本文方法具有更快的收斂速度，能得到更好的優(yōu)化結(jié)果[5]。

本文首先建立攔截器的氣動(dòng)力、氣動(dòng)熱和攔截彈道模型；然后以航程、總熱流量和飛行時(shí)間作為多目標(biāo)優(yōu)化的目標(biāo)函數(shù)，利用基于高斯偽譜法的GPOPS程序[6]對(duì)其進(jìn)行求解；最后利用MOEA/D計(jì)算得到Pareto前沿與優(yōu)化外形。優(yōu)化結(jié)果驗(yàn)證了該方法的有效性。

1 攔截器分析模型的建立

1.1 攔截器幾何外形

根據(jù)任務(wù)時(shí)序，三錐體攔截器會(huì)呈現(xiàn)2種形態(tài)，分別為安裝有整流罩的三錐體外形和具有較高升阻比的雙錐-尾翼組合體外形。在攔截飛行前段，剛進(jìn)入并穿越稀薄大氣區(qū)時(shí)，攔截器為三錐體初始形態(tài)；在進(jìn)入稠密大氣后，尾錐分離，呈現(xiàn)雙錐-尾翼外形，從而為飛行器提供較大升力。

基于參數(shù)化建模的思想，為降低建模復(fù)雜度并減小計(jì)算量，將該三錐體外形參數(shù)設(shè)定為球頭半徑RN，3段長(zhǎng)度L1，L2，L3，以及3段半徑R1，R2，R3。三錐體半徑的增加，可有效降低飛行器的熱流率，但也會(huì)因此而增大飛行器的阻力系數(shù)，降低其氣動(dòng)性能。在分離后的雙錐-尾翼組合體氣動(dòng)外形設(shè)計(jì)中，為簡(jiǎn)化計(jì)算，取尾翼前端位置為第3段錐體的中點(diǎn)。上述2種形態(tài)的飛行器的外形剖面如圖1，2所示。

圖1 三錐體氣動(dòng)外形剖面Fig.1 Aerodynamic profile of triple-cone interceptor

圖2 雙錐-尾翼組合體氣動(dòng)外形剖面Fig.2 Aerodynamic profile of combination

1.2 空氣動(dòng)力計(jì)算模型與驗(yàn)證

高超聲速飛行器的空氣動(dòng)力計(jì)算對(duì)于確定其氣動(dòng)外形、飛行軌跡和性能極為重要。因當(dāng)前的風(fēng)洞試驗(yàn)無法完全模擬高超聲速飛行的真實(shí)條件，故在概念設(shè)計(jì)中，需采用快捷有效的工程估算算法。

由于攔截器在飛行中主要處于連續(xù)流和自由分子流之間的稀薄過渡流區(qū)域，因此可用當(dāng)?shù)貥蚝瘮?shù)等加權(quán)函數(shù)來綜合考慮其動(dòng)力特性。在應(yīng)用當(dāng)?shù)貥蚝瘮?shù)前，需分別確定連續(xù)流和自由分子流區(qū)域的氣動(dòng)特性。

對(duì)于連續(xù)流區(qū)，采用Dahlem-Buck方法[7]估算物面壓力系數(shù)，并進(jìn)行馬赫數(shù)、背風(fēng)區(qū)的經(jīng)驗(yàn)修正，計(jì)算公式為

(1)

式中：Ccontinue為連續(xù)流區(qū)壓力系數(shù)；Ma∞為來流馬赫數(shù)；θ為來流與物面夾角；K為常數(shù)。

對(duì)于自由分子流區(qū)，通常采用經(jīng)驗(yàn)公式擬合方法[8]，具體采用文獻(xiàn)[9]中物面壓力系數(shù)計(jì)算公式

(2)

式中：Crare為自由分子流區(qū)壓力系數(shù)；fn為粘性系數(shù)；S為參考面積；TW為物面溫度；T∞為來流溫度。

在得到飛行器在連續(xù)流區(qū)的氣動(dòng)系數(shù)Ccontinue和非連續(xù)流區(qū)的氣動(dòng)系數(shù)Crare后，利用特定的橋函數(shù)方法將兩者進(jìn)行組合，得到稀薄氣體過渡流區(qū)的氣動(dòng)特性系數(shù)

Ctran=PbKn∞Crare+(1-PbKn∞)Ccontinue

(3)

式中：Pb=sin2φ，其中，φ=π(a1+a2lgKn∞)，此處，a1，a2為待選系數(shù)，Kn∞為來流努森數(shù)。

為驗(yàn)證上文介紹的快速估算算法的可靠性，需將其與計(jì)算流體動(dòng)力學(xué)(CFD)方法的計(jì)算結(jié)果進(jìn)行比較。CFD方法主要包括基于自由分子流理論的直接模擬蒙特卡羅(DSMC)法和基于連續(xù)流理論的納維-斯托克斯(NS)方程法。針對(duì)飛行器稀薄效應(yīng)較為顯著的50 km的飛行高度，本文分別采用DSMC法與NS方程法對(duì)三錐體外形進(jìn)行氣動(dòng)計(jì)算，并與橋函數(shù)工程算法進(jìn)行比較，結(jié)果如圖3所示。

圖3 三錐體氣動(dòng)系數(shù)計(jì)算對(duì)比Fig.3 Comparison of aerodynamic coefficients of triple-cone interceptor

由圖可知：橋函數(shù)法與2種CFD方法的計(jì)算結(jié)果并無較大差別，升力系數(shù)相對(duì)誤差最大不超過20%，阻力系數(shù)相對(duì)誤差最大不超過10%。兩者相對(duì)誤差最大值均在大攻角處，且在攻角小于30°時(shí)相對(duì)誤差均較小。因此，橋函數(shù)法的計(jì)算精度符合要求，下文在估算飛行器氣動(dòng)參數(shù)時(shí)將采用該方法。

1.3 攔截彈道計(jì)算模型

不考慮地球旋轉(zhuǎn)，在側(cè)滑角為0的假設(shè)條件下，攔截器無動(dòng)力3自由度質(zhì)點(diǎn)動(dòng)力學(xué)和運(yùn)動(dòng)學(xué)方程為

(4)

式中：r為飛行器地心距；θ為經(jīng)度；φ為緯度；v為飛行器相對(duì)地球速度；γ為速度傾角；ψ為速度方位角，ψ=0°表示正東方向；σ為傾斜角；D和L分別為阻力加速度和升力加速度；g為地球重力加速度，g=9.8 m/s2。攔截彈道計(jì)算方法采用高斯-偽譜法[10]。

1.4 熱流計(jì)算模型

采用Fay-Riddle公式估算攔截器的駐點(diǎn)熱流，將高溫氣體邊界層偏微分方程轉(zhuǎn)為常微分方程[11]。假設(shè)與氣體熱力學(xué)特性、輸運(yùn)特性有關(guān)的無因次參數(shù)為一系列常數(shù)：普朗特?cái)?shù)Pr=0.71；路易斯數(shù)Le=1.0～2.0；ρsμs/(ρwμw)=1.0，其中，ρs為駐點(diǎn)空氣密度，μs為駐點(diǎn)空氣粘性系數(shù)，ρw為壁面空氣密度，μw為壁面空氣粘性系數(shù)。當(dāng)飛行速度為1.77～7.00 km/s，飛行高度為7.6～37.0 km，壁溫為300～3 000 K時(shí)，利用有限差分法，對(duì)平衡、非平衡和凍結(jié)流進(jìn)行大量的數(shù)值計(jì)算。最終，駐點(diǎn)的熱流密度qs可表示為

(5)

(6)

(7)

式中：Ts為駐點(diǎn)溫度。

2 MOEA/D介紹

在此簡(jiǎn)要介紹MOEA/D的思路和流程。MOEA/D主要將多目標(biāo)的優(yōu)化問題分解為多個(gè)標(biāo)量子目標(biāo)進(jìn)行優(yōu)化，然后使用遺傳算法求解這些子問題。因相鄰兩個(gè)子問題的優(yōu)化解非常相似[2]，故在MOEA/D中，每個(gè)子問題都可使用其相鄰子問題的優(yōu)化信息，從而提高優(yōu)化效率，減小計(jì)算量。其流程如下[12]：

1) 初始化解集與種群。首先，令Pareto解集為?。然后，計(jì)算權(quán)重索引集，將其記為B(i)={i1,i1,…,iT}，其中，λi1,λi2,…，λiT為距離權(quán)重向量λi的歐幾里得距離最小的T個(gè)向量。之后，建立初始種群x1,x2,…,xN，并設(shè)置種群的支配變量Vi=F(xi)。根據(jù)具體問題，初始化最優(yōu)值的解集z=[z1z2…zi]T。

2) 更新解集。從B(i)中隨機(jī)選擇2個(gè)下標(biāo)k和l，利用遺傳算法的特點(diǎn)，從xk和xl中產(chǎn)生新解y，然后更新z。對(duì)j=1,2,…,m，若zj

3) 若滿足循環(huán)終止條件，則停止循環(huán)。

3 攔截器的多目標(biāo)優(yōu)化模型

3.1 優(yōu)化問題設(shè)置

3.1.1 設(shè)計(jì)變量

設(shè)球頭半徑為定值30 mm，設(shè)計(jì)變量為圖1中長(zhǎng)度變量L1，L2，L3，R1，R2，R3，決策變量x為

x=[L1L2L3R1R2R3]T

(8)

3.1.2 優(yōu)化目標(biāo)

該飛行器飛行時(shí)間長(zhǎng)，且要穿越稀薄大氣與稠密大氣區(qū)域，飛行器熱環(huán)境嚴(yán)峻，這要求攔截器攔截時(shí)間盡可能短，并保證攔截范圍盡可能大。因此，設(shè)置優(yōu)化目標(biāo)為最大化航程R、最小化飛行時(shí)間t和最小化飛行總熱流量Q，并將GPOPS的彈道優(yōu)化結(jié)果作為該氣動(dòng)外形的目標(biāo)函數(shù)值，即

minF=[-RtQ]T

(9)

3.1.3 約束條件

過程約束、控制量約束、初始狀態(tài)量和終端條件設(shè)置如下：最大熱流密度為1 600 kW/m2；法向過載nz=Lcosα+Dsinα<5；彈道傾角區(qū)間為[-5°，5°]，末端傾角為0°；攻角區(qū)間為[-10°,30°]；傾側(cè)角區(qū)間為[-70°，70°]；始端速度為7 864 m/s，末端速度為1 812 m/s；始端高度為100 km，末端高度為35 km；始端航向角為90°。

3.2 優(yōu)化變量

飛行器氣動(dòng)外形的設(shè)計(jì)變量的基準(zhǔn)值和上下限見表1。優(yōu)化初始種群規(guī)模為200，最大迭代步數(shù)為100，鄰近個(gè)體的數(shù)量為30，不滿足約束的懲罰因子為1 000。

表1 外形參數(shù)變量取值

4 仿真結(jié)果

圖4為優(yōu)化后的Pareto前沿三維示意圖。針對(duì)攔截器外形進(jìn)行單目標(biāo)優(yōu)化。單目標(biāo)優(yōu)化為多目標(biāo)優(yōu)化的一種特殊情況，即在分解目標(biāo)時(shí)將分解向量中的元素λi設(shè)為1，其余分量設(shè)為0。因此，航程、飛行時(shí)間和總熱流量單目標(biāo)優(yōu)化結(jié)果可近似認(rèn)為是單目標(biāo)最優(yōu)解。得到單目標(biāo)最優(yōu)解后，代入彈道程序進(jìn)行計(jì)算，可得飛行器的飛行軌跡，結(jié)果如圖5所示。圖中：實(shí)線部分為飛行器三錐體形態(tài)的軌跡；虛線部分為組合體形態(tài)的軌跡?；谠搩?yōu)化結(jié)果的飛行器外形如圖6所示。圖中：虛線部分為飛行器的基準(zhǔn)外形。

選擇航程作為單目標(biāo)優(yōu)化問題，航程最遠(yuǎn)彈道優(yōu)化結(jié)果如圖5(a)所示。圖中：實(shí)線部分為單目標(biāo)航程最遠(yuǎn)的氣動(dòng)外形對(duì)應(yīng)的彈道優(yōu)化結(jié)果；虛線部分為基準(zhǔn)外形的最遠(yuǎn)彈道優(yōu)化結(jié)果。因分離后的兩錐-尾翼組合體擁有更大的升阻比，故根據(jù)“理想狀況下升阻比越高，航程越遠(yuǎn)”的規(guī)律，飛行器分離時(shí)間較早。就優(yōu)化結(jié)果而言，其航程高于基準(zhǔn)值約17%，總熱流量增加15%，飛行時(shí)間增加19%。圖6(a)為優(yōu)化后的氣動(dòng)外形與基準(zhǔn)外形的對(duì)比圖。由圖可見：優(yōu)化后的外形長(zhǎng)細(xì)比更大，飛行器的升阻比等氣動(dòng)性能得到有效提高。

圖4 優(yōu)化后的Pareto前沿示意圖Fig.4 Schematic diagram of Pareto frontier after optimization

圖5 攔截彈道計(jì)算結(jié)果Fig.5 Results of intercept trajectory simulation

圖6 氣動(dòng)外形優(yōu)化結(jié)果Fig.6 Results of aerodynamic shape optimization

選擇飛行時(shí)間作為單目標(biāo)優(yōu)化問題，飛行時(shí)間最短彈道優(yōu)化結(jié)果如圖5(b)所示。圖中：實(shí)線部分為單目標(biāo)再入時(shí)間最短的氣動(dòng)外形對(duì)應(yīng)的彈道優(yōu)化結(jié)果；虛線部分為基準(zhǔn)外形的最短飛行時(shí)間優(yōu)化結(jié)果。因三錐體形態(tài)較分離后雙錐-尾翼組合體有更大的阻力系數(shù)，故飛行器在絕大多數(shù)時(shí)候都處于分離前的狀態(tài)，從而實(shí)現(xiàn)速度與高度的快速下降。與基準(zhǔn)外形相比，優(yōu)化后，飛行器航程縮短了7.3%，總熱流量減少了35.0%，飛行時(shí)間減少了9.9%。基于最短飛行時(shí)間的外形優(yōu)化結(jié)果如圖6(b)所示。由圖可見：優(yōu)化后的外形底部半徑略有增大，阻力系數(shù)也因此而有所增大。

選擇總熱流量為單目標(biāo)優(yōu)化問題，總熱流量最小彈道優(yōu)化結(jié)果如圖5(c)所示。圖中：實(shí)線部分為單目標(biāo)總熱流量最小的氣動(dòng)外形對(duì)應(yīng)的彈道優(yōu)化結(jié)果；虛線部分為基準(zhǔn)外形的最小總熱流量?jī)?yōu)化結(jié)果。與圖5(b)類似，飛行器直到最后時(shí)刻才進(jìn)行分離。類似地，其航程縮短了9.5%，總吸熱量減少了15.4%，飛行時(shí)間減少了8.7%。其外形設(shè)計(jì)結(jié)果如圖6(c)所示。更鈍的頭部可有效減小駐點(diǎn)熱流密度，從而減輕飛行器的熱防護(hù)壓力。此外，更大的底部半徑可增大飛行器的阻力系數(shù)，有利于減少飛行器的飛行時(shí)間。

通過對(duì)比3個(gè)單目標(biāo)的優(yōu)化結(jié)果可知，較大的頭部鈍度或較長(zhǎng)的底部半徑可減小飛行器的熱流密度，但會(huì)降低升阻比，進(jìn)而縮短飛行器的航程。由此可見，在單目標(biāo)優(yōu)化中，飛行時(shí)間與總吸熱量成正相關(guān)關(guān)系，但與航程存在沖突，并不存在同時(shí)滿足每個(gè)目標(biāo)函數(shù)均最優(yōu)的解法，因而需要在多目標(biāo)優(yōu)化設(shè)計(jì)中尋找平衡點(diǎn)。

基于Pareto前沿結(jié)果與單目標(biāo)優(yōu)化結(jié)果，選擇一個(gè)3目標(biāo)函數(shù)均優(yōu)于基準(zhǔn)值的外形。3目標(biāo)同時(shí)優(yōu)化的計(jì)算結(jié)果如圖5(d)所示。就優(yōu)化結(jié)果而言，航程增加了12.5%，飛行時(shí)間減少了9.5%，總吸熱量減少了10.0%。然而，由于盡可能簡(jiǎn)化了飛行器的幾何外形參數(shù)，同時(shí)在氣動(dòng)數(shù)據(jù)和飛行氣動(dòng)熱的計(jì)算中采用了快速估算算法，因此，優(yōu)化結(jié)果在一定程度上存在誤差，但仍具一定的參考價(jià)值。

5 結(jié)束語

本文針對(duì)一種可分離的三錐體臨近空間攔截器進(jìn)行了多目標(biāo)優(yōu)化。首先，結(jié)合基于高斯-偽譜法的彈道計(jì)算模型、氣動(dòng)力快速估算算法和氣動(dòng)熱計(jì)算方法，建立飛行器分析模型。之后，采用MOEA/D計(jì)算Pareto前沿，得到3個(gè)單目標(biāo)最優(yōu)解和一個(gè)3目標(biāo)同時(shí)改進(jìn)的較優(yōu)解。優(yōu)化結(jié)果表明：該目標(biāo)優(yōu)化算法具有一定的參考價(jià)值，可為未來的攔截器設(shè)計(jì)提供參考。由于采用了快速估算算法，計(jì)算結(jié)果存在一定誤差，因此，后續(xù)需研究更精確高效的氣動(dòng)模型計(jì)算方法，以及MOEA/D在飛行器優(yōu)化問題上的改進(jìn)方案。