基于MOEA/D的三錐體攔截器氣動外形優化設計

曹 粟，蔣 鋒，李 易

(1.西北工業大學 航天學院，陜西 西安 710072； 2.陜西省空天飛行器設計重點實驗室，陜西 西安 710072； 3.上海航天技術研究院，上海 201109)

0 引言

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

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

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

1 攔截器分析模型的建立

1.1 攔截器幾何外形

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

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

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

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

1.2 空氣動力計算模型與驗證

高超聲速飛行器的空氣動力計算對于確定其氣動外形、飛行軌跡和性能極為重要。因當前的風洞試驗無法完全模擬高超聲速飛行的真實條件，故在概念設計中，需采用快捷有效的工程估算算法。

由于攔截器在飛行中主要處于連續流和自由分子流之間的稀薄過渡流區域，因此可用當地橋函數等加權函數來綜合考慮其動力特性。在應用當地橋函數前，需分別確定連續流和自由分子流區域的氣動特性。

對于連續流區，采用Dahlem-Buck方法[7]估算物面壓力系數，并進行馬赫數、背風區的經驗修正，計算公式為

(1)

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

對于自由分子流區，通常采用經驗公式擬合方法[8]，具體采用文獻[9]中物面壓力系數計算公式

(2)

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

在得到飛行器在連續流區的氣動系數Ccontinue和非連續流區的氣動系數Crare后，利用特定的橋函數方法將兩者進行組合，得到稀薄氣體過渡流區的氣動特性系數

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

(3)

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

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

圖3 三錐體氣動系數計算對比Fig.3 Comparison of aerodynamic coefficients of triple-cone interceptor

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

1.3 攔截彈道計算模型

不考慮地球旋轉，在側滑角為0的假設條件下，攔截器無動力3自由度質點動力學和運動學方程為

(4)

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

1.4 熱流計算模型

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

(5)

(6)

(7)

式中：Ts為駐點溫度。

2 MOEA/D介紹

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

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

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

3) 若滿足循環終止條件，則停止循環。

3 攔截器的多目標優化模型

3.1 優化問題設置

3.1.1 設計變量

設球頭半徑為定值30 mm，設計變量為圖1中長度變量L1，L2，L3，R1，R2，R3，決策變量x為

x=[L1L2L3R1R2R3]T

(8)

3.1.2 優化目標

該飛行器飛行時間長，且要穿越稀薄大氣與稠密大氣區域，飛行器熱環境嚴峻，這要求攔截器攔截時間盡可能短，并保證攔截范圍盡可能大。因此，設置優化目標為最大化航程R、最小化飛行時間t和最小化飛行總熱流量Q，并將GPOPS的彈道優化結果作為該氣動外形的目標函數值，即

minF=[-RtQ]T

(9)

3.1.3 約束條件

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

3.2 優化變量

飛行器氣動外形的設計變量的基準值和上下限見表1。優化初始種群規模為200，最大迭代步數為100，鄰近個體的數量為30，不滿足約束的懲罰因子為1 000。

表1 外形參數變量取值

4 仿真結果

圖4為優化后的Pareto前沿三維示意圖。針對攔截器外形進行單目標優化。單目標優化為多目標優化的一種特殊情況，即在分解目標時將分解向量中的元素λi設為1，其余分量設為0。因此，航程、飛行時間和總熱流量單目標優化結果可近似認為是單目標最優解。得到單目標最優解后，代入彈道程序進行計算，可得飛行器的飛行軌跡，結果如圖5所示。圖中：實線部分為飛行器三錐體形態的軌跡；虛線部分為組合體形態的軌跡。基于該優化結果的飛行器外形如圖6所示。圖中：虛線部分為飛行器的基準外形。

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

圖4 優化后的Pareto前沿示意圖Fig.4 Schematic diagram of Pareto frontier after optimization

圖5 攔截彈道計算結果Fig.5 Results of intercept trajectory simulation

圖6 氣動外形優化結果Fig.6 Results of aerodynamic shape optimization

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

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

通過對比3個單目標的優化結果可知，較大的頭部鈍度或較長的底部半徑可減小飛行器的熱流密度，但會降低升阻比，進而縮短飛行器的航程。由此可見，在單目標優化中，飛行時間與總吸熱量成正相關關系，但與航程存在沖突，并不存在同時滿足每個目標函數均最優的解法，因而需要在多目標優化設計中尋找平衡點。

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

5 結束語

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