楊希祥，周 張，彭 科

(1.國防科技大學(xué)航天科學(xué)與工程學(xué)院，長沙 410073;2.北京理工大學(xué) 宇航學(xué)院，北京 100081)

0 引言

氣動外形設(shè)計(jì)在運(yùn)載火箭總體設(shè)計(jì)中占有極其重要的地位。運(yùn)載火箭外形一般由頭部、箭身和尾舵(小型運(yùn)載火箭)等部分組成。頭部外形設(shè)計(jì)即整流罩外形設(shè)計(jì)對飛行性能影響至關(guān)重要，是外形設(shè)計(jì)的關(guān)鍵［1］。運(yùn)載火箭總體設(shè)計(jì)需重點(diǎn)考慮的氣動阻力、表面?zhèn)鳠崽匦院陀行лd荷容積等因素，都與整流罩外形息息相關(guān)。整流罩外形設(shè)計(jì)在國外飛行器設(shè)計(jì)中歷來受到高度重視［2－4］，國內(nèi)卻鮮有文獻(xiàn)對這一問題進(jìn)行研究。

氣動外形設(shè)計(jì)優(yōu)化過程往往需要多次調(diào)用氣動參數(shù)計(jì)算模型，直接采用高精度計(jì)算模型將面臨嚴(yán)重計(jì)算復(fù)雜性問題，近似建模方法應(yīng)運(yùn)而生。近似建模方法基本思想是用一個簡單逼近函數(shù)近似替代復(fù)雜高精度實(shí)際性能分析模型。根據(jù)近似函數(shù)所能模擬設(shè)計(jì)空間的大小，近似建模方法可分為局部近似方法、全局近似方法和中等范圍近似方法。全局近似方法中的Kriging函數(shù)法以具有樣本點(diǎn)處無偏估計(jì)、良好的高階非線性擬合能力和靈活的近似模型參數(shù)選擇等優(yōu)點(diǎn)，在近似建模領(lǐng)域得到廣泛應(yīng)用。Meunier和Laurenceau 等采用其建立機(jī)翼氣動力近似模型［5－6］;Han等提出一種基于Kriging函數(shù)的變復(fù)雜度近似建模方法，并應(yīng)用于氣動力近似建模［7］;Ahmed等分別研究了Kriging函數(shù)在高超聲速鈍錐體氣動力和氣動熱近似建模中的應(yīng)用［8］。

本文研究小型運(yùn)載火箭整流罩氣動外形設(shè)計(jì)優(yōu)化問題，近似模型構(gòu)造選用Kriging函數(shù)法。

1 Kriging函數(shù)法及其仿真實(shí)現(xiàn)

1.1 參數(shù)估計(jì)與預(yù)測

假定n維輸入變量X與輸出響應(yīng)y之間的函數(shù)關(guān)系y(X)未知，給定組s訓(xùn)練樣本數(shù)據(jù)(X，y)=［(X1，y1)，…，(Xs，ys)］T，結(jié)合線性回歸參數(shù)模型和隨機(jī)過程非參數(shù)模型，采用Kriging函數(shù)可構(gòu)造y(X)的近似模型 y^(X)如下［9］:

式中 fT(X)=［f1(X)，…，fe(X)］為回歸函數(shù);β=［β1，…，βe］T為回歸系數(shù);z(X)為隨機(jī)過程。

工程應(yīng)用中，線性回歸參數(shù)模型常采用二階多項(xiàng)式回歸函數(shù)，

隨機(jī)過程非參數(shù)模型考慮為真實(shí)響應(yīng)與假設(shè)線性模型的偏差，常采用均值為零、高斯相關(guān)函數(shù)的隨機(jī)過程模型，其協(xié)方差函數(shù)為

隨機(jī)變量 y=y(X)和 Ys=［y(X1)，…，y(Xs)］T服從如下多元正態(tài)分布:

其中:

式中 fT為e維回歸函數(shù)行向量;rT為s維相關(guān)函數(shù)行向量;F為s×e維回歸函數(shù)矩陣;R為s×s維相關(guān)函數(shù)矩陣。

y Ys服從如下條件分布:給定訓(xùn)練樣本數(shù)據(jù)Ys下，y的最小均方誤差估計(jì)為yYs的條件均值:

要得到式(6)的Kriging函數(shù)近似模型，需基于訓(xùn)練樣本數(shù)據(jù)估計(jì)近似模型參數(shù)β、θ，通常采用兩步估計(jì)法，β與θ相互獨(dú)立，即滿足下式:

第一步估計(jì)是假定已知相關(guān)函數(shù)模型及參數(shù)θ，得到β的廣義最小二乘估計(jì)(θ):

第二步估計(jì)是得到θ的某個估計(jì)θ^，然后將θ^代替θ，得到β^(θ^)作為β的估計(jì)。相關(guān)參數(shù)θ的估計(jì)方法有交互驗(yàn)證估計(jì)方法和極大似然估計(jì)方法，訓(xùn)練樣本點(diǎn)數(shù)較少時，可選擇交互驗(yàn)證估計(jì)方法，否則選擇極大似然估計(jì)方法。

相關(guān)參數(shù)θ的交互驗(yàn)證估計(jì)值θ^為所有樣本訓(xùn)練點(diǎn)處的均方預(yù)估誤差的最小值點(diǎn):

采用極大似然估計(jì)方法時，假設(shè)訓(xùn)練樣本Ys=［y(X1)，…，y(Xs)］T服從多元正態(tài)分布:

取似然函數(shù)FMLE為訓(xùn)練樣本數(shù)據(jù)Ys概率密度f(Ys):

不考慮常數(shù)項(xiàng)的對數(shù)似然函數(shù)為

將式(15)代入式(14)，得到只與相關(guān)參數(shù)θ相關(guān)的對數(shù)似然函數(shù):

式(18)所示的近似模型不僅可以預(yù)估輸入變量X和輸出響應(yīng)關(guān)系y(X)，而且還可以預(yù)估輸入變量X處一階偏導(dǎo)數(shù):

1．2 仿真實(shí)現(xiàn)

實(shí)際應(yīng)用中，由于訓(xùn)練樣本數(shù)量大，致使高維相關(guān)函數(shù)矩陣的求逆運(yùn)算計(jì)算復(fù)雜，同時訓(xùn)練樣本點(diǎn)選擇不合理會導(dǎo)致相關(guān)函數(shù)矩陣病態(tài)，無法計(jì)算逆矩陣。由于相關(guān)函數(shù)矩陣對稱，屬稀疏矩陣，為避免直接計(jì)算產(chǎn)生較大計(jì)算誤差，對相關(guān)函數(shù)矩陣進(jìn)行Cholesky分解;為避免不同的設(shè)計(jì)變量量級不同產(chǎn)生計(jì)算誤差，對訓(xùn)練樣本數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，具體步驟如下:

(1)對訓(xùn)練樣本數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，得到標(biāo)準(zhǔn)空間［－1，1］s上的訓(xùn)練樣本數(shù)據(jù)

(4)對回歸函數(shù)矩陣F進(jìn)行正交轉(zhuǎn)化并做QR分解:QG=C－1F。

(5)計(jì)算基于標(biāo)準(zhǔn)化訓(xùn)練樣本數(shù)據(jù)的Kriging函數(shù)近似模型:

(6)計(jì)算式(11)均方預(yù)估誤差值fCV(θ)或式(17)對數(shù)似然函數(shù)值LMLE(θ)。

(7)更新相關(guān)參數(shù)比例參數(shù)θ，重復(fù)執(zhí)行步驟(3)～(6)，直到找到目標(biāo)函數(shù)極小值。

(8)基于訓(xùn)練樣本數(shù)據(jù)的Kriging近似模型:

2 整流罩氣動外形設(shè)計(jì)優(yōu)化問題

2．1 設(shè)計(jì)變量選取

研究的運(yùn)載火箭整流罩前罩基線外形如圖1所示。為靈活表示整流罩氣動外形，便于選取控制變量，采用易修改性和局部調(diào)整特性好的非均勻有理B樣條曲線(Non-Uniform Rational B-Spline，NURBS)對整流罩外形進(jìn)行參數(shù)化表示。

控制點(diǎn)選取如圖2所示，除首末2點(diǎn)外，另取4個控制點(diǎn)，共6個控制點(diǎn)。其中，控制點(diǎn)2和控制點(diǎn)1具有相同橫坐標(biāo)(均為0)，目的是保持NURBS曲線光滑性，控制點(diǎn)3距離頭部頂點(diǎn)很近，目的是增強(qiáng)外形多樣性，控制點(diǎn)5與控制點(diǎn)6具有相同縱坐標(biāo)，目的是保持后端光滑性［4］。根據(jù)NURBS曲線性質(zhì)，改變控制點(diǎn)坐標(biāo)即可方便調(diào)整整流罩外形，為此，選取控制點(diǎn)坐標(biāo)為設(shè)計(jì)變量，共4個，具體如下:控制點(diǎn)3的橫、縱坐標(biāo)分別為x1和x2，控制點(diǎn)4的縱坐標(biāo)為x3，控制點(diǎn)5的橫坐標(biāo)為x4，即控制變量:

圖1 整流罩基線外形Fig.1 Primary shape of fairing

圖2 控制頂點(diǎn)選取Fig.2 Selection of control points

設(shè)計(jì)變量邊界值的確定對優(yōu)化設(shè)計(jì)速度和結(jié)果有重要影響。為加快收斂速度，同時確保得到全局最優(yōu)解，首先采用DATCOM氣動特性工程估算程序與Fay和Riddell熱流密度計(jì)算公式，結(jié)合Powell法，求得粗略最優(yōu)解，用于收縮設(shè)計(jì)空間，然后采用基于Kriging函數(shù)近似模型的方法進(jìn)一步精確求解。

DATCOM是美國空軍實(shí)驗(yàn)室研究的一套用于導(dǎo)彈氣動力工程估算的程序，對整流罩這樣的標(biāo)準(zhǔn)細(xì)長體預(yù)估誤差在10%～15%，對DATCOM源程序進(jìn)行校核和修正，形成核心計(jì)算程序，編制輸入輸出文件，嵌入優(yōu)化流程對迭代過程氣動力進(jìn)行計(jì)算。Fay和Riddell熱流近似計(jì)算公式用于求解軸對稱體駐點(diǎn)區(qū)熱流密度:

式中 下標(biāo)w、e、s分別表示壁面條件、邊界層外緣條件、駐點(diǎn)條件;qws為駐點(diǎn)熱流密度;ρw、ρs為密度;μw、μs為粘性系數(shù);hs、hw為比焓;hd為離解比焓;RN為鈍頭半徑。

為便于計(jì)算，Le和hd用下式近似計(jì)算［10］:

2.2 目標(biāo)函數(shù)和約束條件

運(yùn)載火箭大氣層飛行段阻力對運(yùn)載能力有重要影響，而氣動阻力主要集中在頭部。為此，以一級飛行段平均阻力最小為目標(biāo)函數(shù)，選取基線方案3個典型彈道特征點(diǎn)，阻力系數(shù)最大點(diǎn)、動壓最大點(diǎn)和一級發(fā)動機(jī)關(guān)機(jī)點(diǎn)，3點(diǎn)相關(guān)參數(shù)如表1所示。

表1 彈道特征點(diǎn)參數(shù)Table 1 Trajectory characteristic parameters

引入加權(quán)因子對目標(biāo)函數(shù)進(jìn)行處理，約束條件包括3個設(shè)計(jì)點(diǎn)駐點(diǎn)熱流密度最大值、整流罩內(nèi)部容積及設(shè)計(jì)變量約束:

式中 Di、qwi(i=1，2，3)分別表示設(shè)計(jì)點(diǎn)i阻力和駐點(diǎn)熱流密度;qw，max為允許熱流密度最大值;V表示整流罩內(nèi)部容積;Vb表示基線方案整流罩容積值，設(shè)計(jì)變量邊界根據(jù)初步優(yōu)化值和設(shè)計(jì)變量實(shí)際允許的邊界值得到。

3 設(shè)計(jì)優(yōu)化結(jié)果與分析

采用Kriging函數(shù)建立運(yùn)載火箭整流罩氣動性能計(jì)算近似模型，采用正交表型均勻拉丁超立方試驗(yàn)設(shè)計(jì)方法［11］選取樣本點(diǎn)，樣本點(diǎn)選取36個。

樣本訓(xùn)練采用非結(jié)構(gòu)網(wǎng)格和基于N-S方程的數(shù)值計(jì)算方法，湍流模型采用k-ε兩方程模型，空間離散采用上風(fēng)差分格式中的通量差分分裂(FDS)格式，時間分裂采用LU-SGS方法，由于3個特征設(shè)計(jì)點(diǎn)訓(xùn)練樣本總數(shù)較多，采用并行計(jì)算產(chǎn)生氣動力參數(shù)。優(yōu)化方法選取Powell法，初值設(shè)定采用2.1節(jié)所述方法。

式(25)中，w1、w2和 w3分別取為 0.2、0.5、0.3。經(jīng)4次迭代優(yōu)化過程收斂，優(yōu)化方案目標(biāo)函數(shù)和約束條件如表2所示，表中同時給出了基線方案和在3個特征點(diǎn)分別以阻力最小為目標(biāo)函數(shù)設(shè)計(jì)結(jié)果對應(yīng)的目標(biāo)函數(shù)和約束條件值。

優(yōu)化前后整流罩外形(含計(jì)算網(wǎng)格)對比如圖3所示。

特征點(diǎn)2處，優(yōu)化方案和基線方案表面壓力和熱流分布對比如圖4所示。

表2 整流罩外形優(yōu)化設(shè)計(jì)結(jié)果Table 2 Optimization design results of fairing shape

圖3 基線方案與優(yōu)化方案外形對比(含計(jì)算網(wǎng)格)Fig.3 Comparison of shape between basic scheme and optimal scheme(computation grid)

圖4 基線方案和優(yōu)化方案軸向參數(shù)分布對比Fig.4 Comparison of axial parameters between basic scheme and optimal scheme

由圖4可見，由于頭部鈍度減小，優(yōu)化方案在頭部具有更小表面壓力，但在距頭部約45～650 mm處，優(yōu)化方案表面壓力大于基線方案;優(yōu)化方案母線光滑過渡，表面壓力分布脈動較小，有利于整流罩結(jié)構(gòu)設(shè)計(jì);由圖4還可看出，頭部鈍度減小使得優(yōu)化方案駐點(diǎn)熱流密度較高，但滿足設(shè)計(jì)約束。

4 結(jié)論

(1)采用Kriging函數(shù)建立整流罩氣動參數(shù)計(jì)算近似模型，可有效避免整流罩氣動外形設(shè)計(jì)優(yōu)化計(jì)算復(fù)雜性問題，在保證優(yōu)化精度的前提下，實(shí)現(xiàn)較高的優(yōu)化效率。

(2)采用建立的整流罩氣動參數(shù)計(jì)算近似模型，實(shí)現(xiàn)了運(yùn)載火箭整流罩氣動外形設(shè)計(jì)優(yōu)化問題，優(yōu)化方案平均阻力比基線方案減小22.2%，各項(xiàng)約束條件得到良好滿足。

(3)研究成果對于方案論證與方案設(shè)計(jì)階段飛行器氣動外形設(shè)計(jì)優(yōu)化研究具有重要借鑒意義。

［1］龍樂豪.總體設(shè)計(jì)(上)［M］.北京:宇航出版社.1989.

［2］Lee Young Ki，Lee Jae-Woo，Byun Yung-Hwan.The design of space launch vehicle using numerical optimization and inverse method［C］//17th AIAA Applied Aerodynamics Conference.Norfolk，VA，1998:757-767.

［3］Deeppark N R，Ray T，Boyce R R.Evolutionary algorithm shape optimization of a hypersonic flight experiment nose cone［J］.Journal of Spacecraft and Rockets，2008，45(3):428-436.

［4］Lee Jae-Woo，Min Byung-Young，Byun Yung-Hwan，et al.Multipoint nose shape optimization of space launcher using response surface method［J］.Journal of Spacecraft and Rockets，2006，43(1):137-146.

［5］Meunier M.Simulation and optimization of flow control strategies for novel high-lift configurations［J］.AIAA Journal，2009，47(5):1145-1157.

［6］Laurenceau J，Sagaut P.Buliding efficient reponse surfaces of aerodynamic functions with Kriging and Cokriging ［J］.AIAA Journal，2008，46(2):498-507.

［7］Han Z H，Ralf Z，Gortz S.A New Cokriging method for variable-fidelity surrogate modeling of aerodynamic data［C］//48th AIAA Aerospace Sciences Meeting.Orlando，F(xiàn)lorida，2010:1-17.

［8］Ahmed M Y M，Qin N.Metamodels for aerothermodynamic design optimization of hypersonic spiked blunt bodies［C］//48th AIAA Aerospace Sciences Meeting.Orlando，F(xiàn)lorida，2010:1-17.

［9］Matheron G.Principles of geostatistics［J］.Economic Geology，1963(58):1246-1266.

［10］王國雄.彈頭技術(shù)(上)［M?.北京:宇航出版社，1993.

［11］張潤楚，馬長興.正交表型均勻LH設(shè)計(jì)和抽樣［J］.應(yīng)用概率統(tǒng)計(jì)，2001，17(3):243-254.