張海瑞，秦 夢，周國峰，王 浩

(中國運(yùn)載火箭技術(shù)研究院，北京 100076)

新型高超聲速飛行器具有飛行速度快、突防能力強(qiáng)、機(jī)動(dòng)范圍大等優(yōu)勢，是世界主要軍事大國的研究熱點(diǎn)之一。外形優(yōu)化設(shè)計(jì)是影響這一類新型飛行器總體性能的關(guān)鍵技術(shù)之一，旨在尋找一定參數(shù)約束條件下飛行器總體性能最優(yōu)的氣動(dòng)外形，由于其本身的基礎(chǔ)性和全局性，在總體設(shè)計(jì)中占有重要地位[1]。結(jié)合新型飛行器大空域、寬速域特點(diǎn)，有必要考慮飛行剖面的影響，開展飛行器氣動(dòng)-彈道一體化外形優(yōu)化設(shè)計(jì)。

傳統(tǒng)飛行器外形優(yōu)化設(shè)計(jì)方法通常針對典型工況開展氣動(dòng)外形設(shè)計(jì)，給出該工況條件下的外形設(shè)計(jì)方案，在工程中有著廣泛的應(yīng)用。馬洋[2]基于類型函數(shù)/形狀函數(shù)轉(zhuǎn)換(Class/Shape function Transformation,CST)技術(shù)，以升阻比、容積率和壓心位置為目標(biāo)函數(shù)，開展了氣動(dòng)外形的多目標(biāo)優(yōu)化設(shè)計(jì)；Deng等[3]基于二次曲線參數(shù)化模型，以升阻比為目標(biāo)函數(shù)，并以其他總體性能參數(shù)為約束條件，開展了總體氣動(dòng)布局優(yōu)化設(shè)計(jì)。然而，隨著新型飛行器總體性能的不斷提升，傳統(tǒng)方法難以綜合考慮大空域、寬速域飛行剖面特征對氣動(dòng)外形適應(yīng)性的影響，嚴(yán)重制約了飛行器綜合性能的有效提升。針對此問題，在傳統(tǒng)氣動(dòng)外形優(yōu)化模型基礎(chǔ)上，進(jìn)一步耦合了彈道的影響。

本文以新型飛行器為研究對象，以外形參數(shù)為設(shè)計(jì)變量，以最大縱程為優(yōu)化目標(biāo)，綜合利用CST參數(shù)化方法、氣動(dòng)工程估算方法和Radau偽譜法，建立了飛行器氣動(dòng)-彈道一體化模型。由于耦合了彈道的影響，極大增加了一體化模型的復(fù)雜程度。為進(jìn)一步提升優(yōu)化求解效率，引入了Kriging代理模型，提出了一種基于氣動(dòng)-彈道一體化模型的飛行器外形高效全局優(yōu)化設(shè)計(jì)方法。

1 氣動(dòng)-彈道一體化模型

1.1 氣動(dòng)-彈道學(xué)科分析模型

在新型飛行器的方案設(shè)計(jì)階段，外形設(shè)計(jì)參數(shù)決定飛行器的氣動(dòng)特性，并最終影響彈道特性。為了描述這些設(shè)計(jì)參數(shù)對彈道的影響，需要建立一個(gè)既能反映系統(tǒng)特點(diǎn)，又能快速分析和優(yōu)化的學(xué)科分析模型。本文的學(xué)科分析模型主要包括：外形參數(shù)化模型、氣動(dòng)分析模型和彈道分析模型。

1)外形參數(shù)化模型通過一組參數(shù)表達(dá)和約束飛行器幾何外形[2]。采用CST方法[4-6]描述飛行器底部截面曲線η(ψ)，通過類型函數(shù)C(ψ) 定義底部截面曲線的種類，利用形狀函數(shù)S(ψ) 修正類型函數(shù)定義的幾何外形，采用冪函數(shù)[7]控制飛行器俯視面輪廓曲線，進(jìn)而生成參數(shù)化幾何外形，涉及的曲線控制方程如下：

(1)

(2)

其中，(ξ,ψ,η)為飛行器物理坐標(biāo)(x,y,z)歸一化后的參數(shù)坐標(biāo)。式(1)表示升力體底部截面輪廓控制曲線，Nc為類型函數(shù)指數(shù)；式(2)表示飛行器俯視面輪廓控制曲線，xn為俯視面輪廓控制指數(shù)。

2)氣動(dòng)分析模型根據(jù)飛行器外形以及飛行狀態(tài)參數(shù)等計(jì)算氣動(dòng)力系數(shù)。目前常用的氣動(dòng)分析方法主要有風(fēng)洞試驗(yàn)、數(shù)值計(jì)算方法和快速工程估算方法。前兩者能夠更為精確地獲取氣動(dòng)參數(shù)，缺點(diǎn)在于成本高、計(jì)算周期長。為兼顧精度和快速性的計(jì)算要求，本文采用快速工程估算方法求解氣動(dòng)特性，將飛行器表面劃分面元網(wǎng)格，計(jì)算各面元壓力系數(shù)，對所有面元上的氣動(dòng)力進(jìn)行求和，給出飛行器氣動(dòng)特性。其中迎風(fēng)面的計(jì)算使用修正牛頓理論，背風(fēng)面的計(jì)算使用Prandtl-Meyer法[8]。

3)彈道分析模型是在氣動(dòng)分析模型的基礎(chǔ)上，根據(jù)飛行器動(dòng)力學(xué)方程實(shí)現(xiàn)飛行過程仿真。為引入全空、速域變化對氣動(dòng)外形優(yōu)化設(shè)計(jì)的影響，本文將彈道學(xué)科模型處理為軌跡優(yōu)化問題，在一定的約束條件下，以飛行器縱向射程最大為優(yōu)化目標(biāo)，利用Radau偽譜法實(shí)現(xiàn)軌跡優(yōu)化。采用簡化的縱向運(yùn)動(dòng)方程作為彈道計(jì)算模型，飛行器在此階段只受氣動(dòng)力和重力影響，不考慮地球的自轉(zhuǎn)與非球形攝動(dòng)項(xiàng)等因素帶來的影響。用飛行速度V、速度傾角θ、地心距r和射程角φ描述速度系下的飛行器運(yùn)動(dòng)方程：

(3)

式中，CL、CD分別為升力系數(shù)和阻力系數(shù)，m為飛行器質(zhì)量，q為飛行動(dòng)壓，SM為飛行器參考面積。針對縱向動(dòng)力學(xué)方程組，采用Radau偽譜法開展快速軌跡優(yōu)化。Radau偽譜法[9-11]的基本原理是先將狀態(tài)變量和控制變量在(-1,1]上的N個(gè)LGR (Legendre-Gauss-Radau)配點(diǎn)上進(jìn)行離散，然后以N個(gè)配點(diǎn)和初始點(diǎn)τ0=-1為節(jié)點(diǎn)構(gòu)造Lagrange插值多項(xiàng)式來逼近狀態(tài)和控制變量，把最優(yōu)控制問題轉(zhuǎn)化成非線性規(guī)劃問題，再使用某種非線性規(guī)劃算法求得使性能指標(biāo)最優(yōu)的參數(shù)，最后代回原問題進(jìn)行求解。

1.2 一體化優(yōu)化模型

1)目標(biāo)函數(shù)。外形設(shè)計(jì)參數(shù)對飛行器容積率ηv和氣動(dòng)性能具有顯著影響[2]，為了結(jié)果的可比較性，將飛行器的容積率固定，以最大縱程L為優(yōu)化目標(biāo)開展氣動(dòng)-彈道一體化優(yōu)化設(shè)計(jì)。

2)設(shè)計(jì)變量。CST參數(shù)化外形的設(shè)計(jì)參數(shù)中，對升阻比影響最大的5個(gè)參數(shù)依次為底部截面高度xHl、xHu，俯視面輪廓控制指數(shù)xn，底部截面類型函數(shù)指數(shù)xNcu、xNcl[2]。文獻(xiàn)[10]指出，射程最大的飛行控制方式是最大升阻比，因此設(shè)計(jì)變量選為X=(xHu,xHl,xNcu,xNcl,xn)?？紤]飛行器總體規(guī)模及外形尺寸限制，設(shè)飛行器最大寬度和最大長度為定值。

(4)

式中：Kn為飛行器相關(guān)經(jīng)驗(yàn)參數(shù)，與飛行器頭部駐點(diǎn)曲率半徑、飛行器外形及材料相關(guān)；ρ為大氣密度。

綜上，優(yōu)化模型可表示為：

(5)

式中，Re為地球平均半徑。

2 高效全局優(yōu)化方法

由于引入了Radau偽譜法軌跡優(yōu)化模塊，極大地增加了一體化模型的復(fù)雜度，增大了單次迭代的計(jì)算時(shí)間。針對這一問題，引入Kriging代理模型代替設(shè)計(jì)空間內(nèi)的響應(yīng)值。Kriging代理模型[12-13]是一種基于隨機(jī)過程的無偏估計(jì)模型，其預(yù)測的響應(yīng)值服從正態(tài)分布，不僅能給出對未知函數(shù)的預(yù)估值，還能給出預(yù)估值的誤差估計(jì)，在工程中有著廣泛的應(yīng)用。

采用高效全局優(yōu)化(Efficient Global Optimization,EGO)方法[14-15]進(jìn)行求解與優(yōu)化。EGO方法作為一種全局優(yōu)化算法，具有很高的效率和精度[16]，其基本原理是通過期望改善函數(shù)(Expected Improvement Function,EIF)描述代理模型預(yù)測值優(yōu)于真實(shí)值的可能性，并在EIF最大處加入樣本點(diǎn)以更新代理模型。

圖1 高效全局優(yōu)化方法流程Fig.1 Flow chart of efficient global optimization

然后，以初始樣本點(diǎn)為輸入、對應(yīng)的最大縱程L為輸出，構(gòu)建Kriging代理模型為：

(6)

(7)

式中的均值和方差都可以由式(8)[14]計(jì)算：

(8)

(9)

EIF(X)=E(P(X))

(10)

(11)

3 算例

以升力體構(gòu)型飛行器開展優(yōu)化設(shè)計(jì)，飛行器質(zhì)量為m=900 kg，參考面積SM考慮飛行器底部最大截面積，容積率ηv不小于0.292。飛行器最大寬度和長度分別為4 m和2.5 m，外形設(shè)計(jì)變量取值范圍如表1所示。

表1 外形設(shè)計(jì)變量及變化范圍Tab.1 Design variables and ranges of variation

表2 偽譜法狀態(tài)變量初值和終值設(shè)置Tab.2 Initial and final value settings of state variables in pseudospectral method

利用高效全局優(yōu)化方法進(jìn)行氣動(dòng)-彈道一體化優(yōu)化設(shè)計(jì)。在初始階段采用全因子采樣方案對初始參數(shù)設(shè)計(jì)空間進(jìn)行等密度采樣，構(gòu)建初始Kriging代理模型。為對比驗(yàn)證，利用高效全局優(yōu)化方法，開展Ma=10、α=8°工況下的氣動(dòng)外形優(yōu)化設(shè)計(jì)，構(gòu)建初始代理模型的方法同上，進(jìn)而將優(yōu)化結(jié)果代入彈道模型，采用Radau偽譜法以最大射程為目標(biāo)開展快速軌跡優(yōu)化。優(yōu)化結(jié)果對比如表3所示，外形對比如圖2所示。

表3 氣動(dòng)優(yōu)化與氣動(dòng)-彈道一體化優(yōu)化結(jié)果對比Tab.3 Comparison of results of aerodynamic-only optimization with aerodynamic and trajectory integrated optimization

由表3和圖2可以看出，兩種方法給出的俯視面輪廓控制指數(shù)xn相同，俯視面輪廓外形完全一致；一體化優(yōu)化方法給出的底部截面高度相對較低，通過優(yōu)化底部截面類型函數(shù)，在保證容積率相同條件下，獲取較優(yōu)氣動(dòng)外形。

為了進(jìn)一步說明氣動(dòng)-彈道一體化模型的優(yōu)勢，給出AeroOpt-8和JointOpt升力體的升阻比隨攻角變化曲線，如圖3所示?？梢钥吹剑?dāng)攻角α=7.71°時(shí)，AeroOpt-8升力體的升阻比最大，為4.604 8；當(dāng)攻角α=7.39°時(shí)，JointOpt升力體的升阻比最大，為4.602 3。AeroOpt-8和JointOpt兩者的最大升阻比相同，然而JointOpt升力體的攻角-升阻比曲線更平緩，這是因?yàn)锳eroOpt-8只考慮了攻角α=8°時(shí)的優(yōu)化外形，而氣動(dòng)-彈道一體化模型耦合了彈道，考慮了飛行過程攻角擾動(dòng)對外形設(shè)計(jì)的影響，得到全空、速域下的優(yōu)化外形。

(a) 氣動(dòng)優(yōu)化外形(a) Aerodynamic optimization shape

圖3 升阻比隨攻角變化曲線對比Fig.3 Comparison of lift-to-drag ratio curves

分別求解AeroOpt-8和JointOpt升力體的氣動(dòng)特性，最大射程對比如表4所示，偽譜法控制飛行攻角的變化曲線對比如圖4所示。JointOpt的最大升阻比略小于AeroOpt-8的最大升阻比，但前者的最大縱向射程比后者提升了28.04 km；JointOpt升力體的飛行攻角整體小于AeroOpt-8，二者的飛行攻角分別在各自最大升阻比對應(yīng)攻角附近小幅度振蕩。結(jié)合圖表可以說明，當(dāng)飛行攻角在最大升阻比對應(yīng)攻角附近發(fā)生擾動(dòng)時(shí)，JointOpt升阻比的降低量比AeroOpt-8更小，前者受到攻角擾動(dòng)的影響比AeroOpt-8更小，說明采用高效全局優(yōu)化算法針對氣動(dòng)-彈道一體化模型開展優(yōu)化設(shè)計(jì)可以提升飛行器對工況變化的適應(yīng)能力，避免因攻角擾動(dòng)而影響升阻比最終導(dǎo)致射程的降低，進(jìn)一步驗(yàn)證了氣動(dòng)-彈道一體化模型能夠綜合考慮全空、速域下優(yōu)化外形的結(jié)論。

圖4 飛行攻角變化曲線對比Fig.4 Comparison of the angle of attack

表4 最大升阻比及最大縱向射程對比Tab.4 Comparison of the maximum lift-to-drag ratio and maximum trajectory

4 結(jié)論

1)以新型飛行器為研究對象，為考慮大空域、寬速域飛行剖面特征對氣動(dòng)外形適應(yīng)性的影響，采用CST參數(shù)化建模方法結(jié)合快速工程估算方法獲取氣動(dòng)特性，通過偽譜法給出彈道特性，建立飛行器氣動(dòng)-彈道一體化模型，提升了外形優(yōu)化模型精度。

2)針對這一模型的特點(diǎn)，以最大縱向射程為優(yōu)化目標(biāo)，以外形參數(shù)為設(shè)計(jì)變量，并明確約束條件，給出了基于Kriging的氣動(dòng)外形優(yōu)化方法，實(shí)現(xiàn)在多參數(shù)約束條件下的飛行器氣動(dòng)外形高效全局優(yōu)化，提升了復(fù)雜外形優(yōu)化模型的求解效率。

3)以升力體構(gòu)型飛行器為例，開展氣動(dòng)外形優(yōu)化設(shè)計(jì)，結(jié)果表明，氣動(dòng)-彈道一體化模型能夠較好描述大空域、寬速域的氣動(dòng)和彈道特征，有效提升新型飛行器氣動(dòng)外形設(shè)計(jì)精度和水平，具有一定的工程應(yīng)用價(jià)值。