楊希祥，周 張，彭 科

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

0 引言

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

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

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

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

1.1 參數(shù)估計與預(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ù)模型和隨機過程非參數(shù)模型，采用Kriging函數(shù)可構(gòu)造y(X)的近似模型 y^(X)如下［9］:

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

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

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

隨機變量 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的最小均方誤差估計為yYs的條件均值:

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

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

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

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

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

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

不考慮常數(shù)項的對數(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 仿真實現(xiàn)

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

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

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

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

(6)計算式(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è)計優(yōu)化問題

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

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

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

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

圖2 控制頂點選取Fig.2 Selection of control points

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

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

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

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

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

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

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

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

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

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

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

樣本訓(xùn)練采用非結(jié)構(gòu)網(wǎng)格和基于N-S方程的數(shù)值計算方法，湍流模型采用k-ε兩方程模型，空間離散采用上風(fēng)差分格式中的通量差分分裂(FDS)格式，時間分裂采用LU-SGS方法，由于3個特征設(shè)計點訓(xùn)練樣本總數(shù)較多，采用并行計算產(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個特征點分別以阻力最小為目標(biāo)函數(shù)設(shè)計結(jié)果對應(yīng)的目標(biāo)函數(shù)和約束條件值。

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

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

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

圖3 基線方案與優(yōu)化方案外形對比(含計算網(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è)計;由圖4還可看出，頭部鈍度減小使得優(yōu)化方案駐點熱流密度較高，但滿足設(shè)計約束。

4 結(jié)論

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

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

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

［1］龍樂豪.總體設(shè)計(上)［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］.應(yīng)用概率統(tǒng)計，2001，17(3):243-254.