基于代理模型的高效全局低音爆優(yōu)化設(shè)計方法

喬建領(lǐng)，韓忠華，宋文萍

西北工業(yè)大學(xué) 航空學(xué)院 翼型葉柵空氣動力學(xué)國家級重點實驗室，西安 710072

第一代超聲速客機商業(yè)運營的失敗，很大程度上是由音爆噪聲等環(huán)境問題造成的。隨著航空科學(xué)技術(shù)的進步，20世紀(jì)90年代中后期又掀起了一輪低音爆超聲速客機研究熱潮。發(fā)展新一代低音爆超聲速客機主要面臨兩個關(guān)鍵問題：超聲速飛行過程中地面音爆強度計算和低音爆氣動布局優(yōu)化設(shè)計[1]。

音爆是飛行器在超聲速飛行時產(chǎn)生的一種特有的聲學(xué)現(xiàn)象。在超聲速飛行條件下，飛機對氣流擾動會形成激波和膨脹波；激波和膨脹波在傳播過程中經(jīng)過演化、匯聚、耗散等作用后傳至地面時，人們會聽到類似“爆炸”的聲音。在對音爆進行研究時，一般將音爆從飛機到地面的傳播過程分為近場、中場和遠(yuǎn)場。其中，遠(yuǎn)場(地面)音爆強度水平是超聲速飛機設(shè)計中最須關(guān)注的問題之一。國外從20世紀(jì)50年代開始，對音爆現(xiàn)象作了系統(tǒng)性研究，發(fā)展了基于線化理論的音爆預(yù)測方法[2-6]、簡化音爆預(yù)測方法[7]、波形參數(shù)法[8]和Burgers方程法[9]。近年來，隨著計算流體力學(xué)(CFD)技術(shù)的發(fā)展，CFD被應(yīng)用于音爆預(yù)測。目前國外學(xué)者發(fā)展了CFD與波形參數(shù)法[10]、全流場CFD求解法[11]和CFD與Burgers方程匹配法[12-13]。針對音爆預(yù)測方法，國內(nèi)研究人員也開展了大量研究[1, 14-16]。

低音爆氣動布局優(yōu)化設(shè)計是通過優(yōu)化算法得到音爆響應(yīng)較低的外形。馮曉強等[17]運用遺傳算法開展了機體容積和音爆響應(yīng)的多目標(biāo)優(yōu)化研究，并得到考慮音爆響應(yīng)的多目標(biāo)問題的Pareto解。Choi等[18]在多學(xué)科飛機綜合分析工具(PASS)基礎(chǔ)上通過建立飛行性能和音爆的多可信度響應(yīng)面模型，并運用單純形方法(Nelder-Mead Simplex)進行了低音爆超聲速公務(wù)機的優(yōu)化。優(yōu)化過程中以建立的響應(yīng)面模型代替直接的飛行性能分析和音爆分析，從而減少費時的數(shù)值分析次數(shù)。Chan[19]在其博士論文中運用多目標(biāo)遺傳算法等優(yōu)化方法，對兩款超聲速公務(wù)機進行了低阻力和低音爆的優(yōu)化設(shè)計。Farhat等[20]將Whitham修正線化音爆預(yù)測理論與梯度優(yōu)化算法相結(jié)合，發(fā)展了適用于超聲速客機的快速低音爆優(yōu)化設(shè)計方法。Aftosmis等[21]將基于Adjoint的梯度優(yōu)化方法應(yīng)用于低音爆優(yōu)化設(shè)計，完成了灣流公司超聲速公務(wù)機的優(yōu)化設(shè)計，地面音爆響應(yīng)降低了約10 dB。Wintzer和Kroo[22]也將Adjoint的梯度優(yōu)化方法應(yīng)用于低音爆優(yōu)化設(shè)計，完成了低音爆公務(wù)機的概念設(shè)計。Rallabhandi[23-24]推導(dǎo)出耦合CFD近場計算與基于Burgers方程遠(yuǎn)場傳播模型的伴隨方程，并應(yīng)用到超聲速客機的低音爆優(yōu)化設(shè)計中。

由上述文獻調(diào)研可知，在低音爆優(yōu)化設(shè)計領(lǐng)域，遺傳算法等啟發(fā)式優(yōu)化和基于Adjoint的梯度優(yōu)化是目前應(yīng)用較多的方法[17-27]。遺傳算法等啟發(fā)式優(yōu)化算法雖然具有全局優(yōu)化能力，但是優(yōu)化效率低，在引入高精度的音爆預(yù)測方法(需要求解CFD)時，其計算量難以承受。而梯度優(yōu)化方法雖然效率高，但屬于局部優(yōu)化，優(yōu)化結(jié)果依賴于起始點的選擇。

為了解決低音爆優(yōu)化中全局性和高效性難以兼顧的問題，本文將最新發(fā)展的代理優(yōu)化(Surrogate-Based Optimization, SBO)算法應(yīng)用于低音爆優(yōu)化設(shè)計，發(fā)展了一種高效全局的低音爆優(yōu)化設(shè)計方法。其中，遠(yuǎn)場音爆分析可采用高精度的方法，但為了方便研究代理優(yōu)化在低音爆領(lǐng)域的應(yīng)用潛力，本文采用Whitham方法對遠(yuǎn)場音爆強度進行分析。選取NASA圓錐體模型對低音爆優(yōu)化方法進行驗證，并與梯度優(yōu)化和遺傳算法的結(jié)果進行比較，之后將該方法應(yīng)用于2014年音爆預(yù)測大會翼身組合體標(biāo)模的低音爆優(yōu)化設(shè)計，展示了該方法在復(fù)雜外形低音爆優(yōu)化設(shè)計方面的潛力。

1 音爆預(yù)測方法

Whitham發(fā)展的線化音爆預(yù)測理論是根據(jù)超聲速旋成體理論修正得到的(本文稱為“Whitham理論”)，不用求解流場，計算簡單、速度快，計算結(jié)果精度比較合理，很適合于飛機初步設(shè)計階段的低音爆布局優(yōu)化和低音爆優(yōu)化方法研究。Whitham理論詳細(xì)推導(dǎo)過程參見文獻[2]。

1.1 線化音爆預(yù)測理論

Whitham將線化超聲速細(xì)長旋成體理論中特征線是相互平行的假設(shè)進行非線性修正，給出精度較高的特征線描述，獲得用于計算近場壓強信號和遠(yuǎn)場音爆強度的公式。其中，計算近場壓強信號的公式為

(1)

(2)

其中：S″(ξ)為軸向站位ξ處馬赫錐所截有效截面積對ξ的二階導(dǎo)數(shù)；y為特征線相關(guān)參數(shù)[2]。

近場壓強信號在向遠(yuǎn)場傳播的過程中，由于激波的匯聚，擾動壓強在一定條件下會在遠(yuǎn)場形成N型波。遠(yuǎn)場N型波強度的計算公式為

(3)

式(3)為自由空間中過壓值的計算式，當(dāng)觀測點位于地面時，考慮地面反射效應(yīng)因子KR，有

(4)

1.2 算例驗證

根據(jù)線化音爆預(yù)測方法開發(fā)了相應(yīng)程序，并選取NASA多段圓錐體模型的試驗數(shù)據(jù)[28]進行驗證。圖 1為模型示意圖，其半徑公式為

r=

(5)

式中：l=2 in(1 in≈2.54 cm)；x為圓錐體軸向坐標(biāo)。

計算時來流馬赫數(shù)分別為1.26、1.41和2.01，近場音爆信號的觀測位置為模型正下方離軸線h=10l處(見圖 2)。圖 3為應(yīng)用線化音爆預(yù)測方法計算的近場壓強信號與試驗值的對比。由計算結(jié)果與試驗值的對比可知，雖然較高馬赫數(shù)下(圖 3(c))的計算結(jié)果與試驗值還存在一定誤差，但總體而言，3種狀態(tài)下計算值與試驗值均吻合較好，驗證了本文開發(fā)的音爆預(yù)測程序的正確性。同時可以看到，由于該音爆預(yù)測方法沒有考慮大氣黏性效應(yīng)，所以圖 3中計算所得近場壓強信號的后激波過強。

圖1 NASA圓錐體模型Fig.1 NASA cone model

圖2 NASA圓錐體模型近場壓強信號試驗觀測位置Fig.2 Experimentally observed position of near-field pressure signature for NASA cone model

圖3 近場壓強信號計算值與試驗值的對比Fig.3 Comparison between calculation and experimental data of near-field pressure signature

2 代理優(yōu)化方法及低音爆優(yōu)化設(shè)計

本節(jié)概述代理優(yōu)化方法的幾個關(guān)鍵要素和基于代理模型的低音爆優(yōu)化設(shè)計流程。代理優(yōu)化的關(guān)鍵要素包括試驗設(shè)計、代理模型及模型訓(xùn)練、優(yōu)化加點及約束處理、優(yōu)化收斂判斷等，其詳細(xì)內(nèi)容見文獻[29]。

2.1 試驗設(shè)計

樣本點集是建立代理模型的基礎(chǔ)?，F(xiàn)代試驗設(shè)計方法主要有擬蒙特卡洛方法、準(zhǔn)蒙特卡洛方法、拉丁超立方抽樣(Latin Hypercube Sampling, LHS)、正交試驗設(shè)計和均勻設(shè)計等[30]，其中較為流行的是拉丁超立方和均勻設(shè)計方法。本文在進行優(yōu)化時選擇拉丁超立方抽樣方法。

2.2 代理模型及模型訓(xùn)練

代理模型是指在分析和優(yōu)化過程中可“代替”相對復(fù)雜和費時的物理分析模型的一種近似數(shù)學(xué)模型。目前已發(fā)展出多種模型，本文在優(yōu)化過程中選擇Kriging模型[31-33]。

Kriging模型將未知函數(shù)看作某一高斯隨機過程的具體表現(xiàn)，相應(yīng)隨機函數(shù)Y(x)的表達(dá)式為

Y(x)=β0+Z(x)

(6)

式中：β0為未知函數(shù)，是Y(x)的數(shù)學(xué)期望；Z(x)為均值為0、方差為σ2的靜態(tài)隨機過程。

Kriging模型的具體建模過程詳見文獻[29]。建模過程涉及相關(guān)函數(shù)的選擇和計算，本文選擇的相關(guān)函數(shù)為“三次樣條函數(shù)”，其表達(dá)式為

Rk(θk,x(i)-x(j))=

(7)

2.3 優(yōu)化加點準(zhǔn)則

優(yōu)化加點準(zhǔn)則是指根據(jù)所建立的代理模型增加新樣本點的準(zhǔn)則。加點的目的是擴充樣本集合，提高代理模型的預(yù)測精度。代理優(yōu)化的本質(zhì)是采用某種加點準(zhǔn)則產(chǎn)生新樣本點，并使這些樣本點逐步落入最優(yōu)點或其附近。目前已經(jīng)有多種加點準(zhǔn)則，較為流行的有EI(Expected Improvement)準(zhǔn)則、MSP(Minimum of Surrogate Prediction)準(zhǔn)則、MSE(Mean Squared Error)準(zhǔn)則、LCB(Lower Confidence Bounding)準(zhǔn)則和PI(Probability of Improvement)準(zhǔn)則等[34-36]。本文選擇“EI+MSP”的組合加點準(zhǔn)則。

2.4 優(yōu)化收斂判斷條件

本文綜合采用文獻[29]中所述的4種判斷方法進行優(yōu)化收斂判斷。

2.5 低音爆優(yōu)化設(shè)計流程

將代理優(yōu)化方法應(yīng)用于低音爆優(yōu)化設(shè)計，以提高低音爆優(yōu)化的效果和效率，特別是當(dāng)遠(yuǎn)場音爆計算方法涉及流場求解時其效果更為顯著。圖4為運用代理優(yōu)化方法進行低音爆優(yōu)化設(shè)計的流程，低音爆優(yōu)化方法不單指前面介紹的Whitham理論，也可以是其他精度較高的遠(yuǎn)場音爆預(yù)測算法。本文為方便驗證代理優(yōu)化應(yīng)用于低音爆優(yōu)化的可行性和應(yīng)用于復(fù)雜外形低音爆優(yōu)化的潛力，選擇計算速度高、精度可接受的Whitham理論進行遠(yuǎn)場音爆分析。在低音爆優(yōu)化過程中，可將遠(yuǎn)場N型波峰值或遠(yuǎn)場感覺聲壓級作為優(yōu)化目標(biāo)，對飛機外形進行優(yōu)化設(shè)計[22, 24-26]。

圖4 基于代理模型的低音爆優(yōu)化設(shè)計流程Fig.4 Framework of surrogate-based low boom optimization design

3 代理優(yōu)化應(yīng)用于低音爆優(yōu)化設(shè)計的算例

本節(jié)將代理優(yōu)化方法應(yīng)用于NASA圓錐體模型和2014年音爆預(yù)測大會翼身組合體標(biāo)準(zhǔn)模型的低音爆優(yōu)化設(shè)計。以NASA圓錐體優(yōu)化算例來驗證代理優(yōu)化應(yīng)用的適用性和優(yōu)越性，以翼身組合體優(yōu)化算例來說明代理優(yōu)化在低音爆優(yōu)化領(lǐng)域的應(yīng)用潛力。

3.1 NASA圓錐體模型的低音爆優(yōu)化

3.1.1 外形參數(shù)化及優(yōu)化數(shù)學(xué)模型

以1.2節(jié)中的NASA圓錐體模型(圖 1)為基準(zhǔn)外形進行低音爆優(yōu)化設(shè)計。外形參數(shù)化及優(yōu)化問題的數(shù)學(xué)模型建立過程如下。

圓錐體模型參數(shù)化如圖 5所示。r1、r2分別表示頭部圓錐體的最大半徑和整個模型的最大半徑，l1、l2、l3分別表示NASA圓錐體模型頭部圓錐體長度、中間圓柱長度和后部圓臺長度。模型參數(shù)化后，其半徑沿軸向的分布為

(8)

表1為表征NASA圓錐體模型的5個參數(shù)的取值范圍。

圖5 NASA多段圓錐體模型參數(shù)化示意圖Fig.5 Sketch of parameterization of NASA stepped cone model

表1NASA圓錐體模型低音爆優(yōu)化設(shè)計變量的取值范圍

Table1RangeofdesignvariablesforlowboomoptimizationdesignofNASAconemodel

Typel1/inl2/inl3/inr1/inr2/inBaseline0．51．00．50．50．08/π0．04/πUpperbound0．20．50．20．050．09Lowerbound1．01．51．00．090．13

優(yōu)化設(shè)計的目標(biāo)函數(shù)為遠(yuǎn)場(h/l=1 000，見圖 6)音爆N型波的峰值，3種設(shè)計狀態(tài)的來流馬赫數(shù)分別為1.26、1.60和2.00。相應(yīng)的數(shù)學(xué)模型為

(9)

式中：V0=0.038 047 in3為基準(zhǔn)模型的體積。前2個約束(C1、C2)為多段圓錐體的長度約束；最后1個約束(C3)為體積約束，以保證l1、l2、l3三段總體積不減。

圖6 NASA圓錐體模型低音爆優(yōu)化中遠(yuǎn)場N型波的計算位置Fig.6 Calculated position of far-field N-wave for low-boom optimization of NASA cone model

3.1.2 低音爆優(yōu)化設(shè)計結(jié)果

低音爆優(yōu)化設(shè)計流程是在課題組開發(fā)的代理優(yōu)化軟件SurroOpt[37]上，結(jié)合1.1節(jié)中介紹的線性音爆預(yù)測方法建立的。優(yōu)化過程中，采用LHS試驗設(shè)計得到20個初始樣本點，并建立Kriging模型，加點準(zhǔn)則為“EI+MSP”組合加點，最大樣本點個數(shù)為300。3種設(shè)計狀態(tài)的收斂歷程如圖 7所示，圖中縱坐標(biāo)“Obj”為優(yōu)化目標(biāo)響應(yīng)值，橫坐標(biāo)“Evaluation”為音爆分析次數(shù),“DoE”表示通過試驗設(shè)計獲得的初始樣本點。

圖7 3種設(shè)計狀態(tài)下NASA圓錐體模型低音爆優(yōu)化設(shè)計中目標(biāo)函數(shù)的收斂歷程Fig.7 Convergence histories of objective function values at three design states for low boom optimization design of NASA cone model

表2 優(yōu)化前后NASA圓錐體模型的設(shè)計變量及遠(yuǎn)場聲壓級對比Table 2 Comparison of design variables and far-field sound pressure levels between baseline and optimal NASA cone models

注：表中數(shù)據(jù)為四舍五入后的值；約束≥0表示滿足約束，當(dāng)約束=0時，表示該約束為主動約束。

3.1.3 代理優(yōu)化、梯度優(yōu)化和遺傳算法的結(jié)果對比

為說明代理優(yōu)化算法在低音爆優(yōu)化設(shè)計領(lǐng)域中應(yīng)用的優(yōu)越性，分別采用遺傳算法和梯度優(yōu)化算法(序列二次規(guī)劃(SQP)[38])對Ma=1.26狀態(tài)下的NASA圓錐體模型進行低音爆優(yōu)化設(shè)計，并與代理優(yōu)化結(jié)果進行對比。遺傳算法的種群規(guī)模為500，最大進化代數(shù)為100，交叉概率為0.9，變異概率為0.05。由于梯度優(yōu)化算法的優(yōu)化結(jié)果與起始點的選擇密切相關(guān)，在運用SQP算法進行低音爆優(yōu)化時，考慮基準(zhǔn)外形、最優(yōu)外形附近及其他位置選取了4個起始點(分別記為x1、x2、x3、x4)，其具體取值及相應(yīng)點處的約束和目標(biāo)函數(shù)值如表 3所示。表 4為代理優(yōu)化結(jié)果與遺傳算法和SQP算法的優(yōu)化結(jié)果對比，“SQP+xi”表示以xi(i=1,2,3,4)為起始點的SQP算法，最右邊一列中“147f”表示147次的目標(biāo)函數(shù)計算，“127g”表示127次的目標(biāo)函數(shù)梯度計算。

圖8 3種設(shè)計狀態(tài)下NASA圓錐體模型優(yōu)化前后的外形對比Fig.8 Comparison of geometries between baseline and optimal NASA cone models at three design states

圖9 優(yōu)化前后NASA圓錐體模型的遠(yuǎn)場N型波和近場壓強信號對比(Ma=1.26)Fig.9 Comparison of far-field N-wave and near-field pressure signatures between baseline and optimal NASA cone models (Ma=1.26)

圖10 優(yōu)化前后NASA圓錐體模型的遠(yuǎn)場N型波和近場壓強信號對比(Ma=1.60)Fig.10 Comparison of far-field N-wave and near-field pressure signatures between baseline and optimal NASA cone models (Ma=1.60)

圖11 優(yōu)化前后NASA圓錐體模型的遠(yuǎn)場N型波和近場壓強信號對比(Ma=2.00)Fig.11 Comparison of far-field N-wave and near-field pressure signatures between baseline and optimal NASA cone models (Ma=2.00)

由表 3和表 4可知，針對該低音爆優(yōu)化設(shè)計問題，SQP算法的優(yōu)化結(jié)果在一定程度上依賴于起始點的選擇，起始點不同，優(yōu)化解不同。分別以基準(zhǔn)外形(x1)和最優(yōu)外形附近點(x4)為起始點時，遠(yuǎn)場音爆N型波強度基本降低到了最低水平；而以其他位置(x2和x3)為起始點時，優(yōu)化外形相對于基準(zhǔn)外形的遠(yuǎn)場音爆N型波強度雖然都有不同程度的減弱，但都沒有達(dá)到全局最優(yōu)。出現(xiàn)這一現(xiàn)象的原因是設(shè)計空間內(nèi)存在多極值，梯度優(yōu)化易陷入局部最優(yōu)。此外，由表 4可知，遺傳算法雖然具有全局優(yōu)化能力，但在優(yōu)化過程中需要大量的目標(biāo)函數(shù)計算，造成巨大的信息浪費，因此其優(yōu)化效率很低。

綜上所述，與梯度優(yōu)化相比，代理優(yōu)化算法的結(jié)果不依賴于起始點的選擇，能夠獲得全局最優(yōu)解；而與遺傳算法相比，代理優(yōu)化算法的效率高了2個量級以上。結(jié)果表明，代理優(yōu)化算法在低音爆優(yōu)化設(shè)計領(lǐng)域中具有明顯的優(yōu)勢。

表3 SQP算法起始點位置的約束和目標(biāo)函數(shù)值Table 3 Constraints and objective function values at start points of SQP algorithm

注：x1為基準(zhǔn)外形點，x4為最優(yōu)外形附近點，x2、x3為其他位置點。

表4 3種優(yōu)化算法對相同低音爆優(yōu)化問題的優(yōu)化結(jié)果對比Table 4 Comparison of results of three optimization algorithms for same low boom optimization problem

3.2 翼身組合體外形的低音爆優(yōu)化

2014年1月，AIAA第1屆音爆預(yù)測大會在美國馬里蘭州舉行，旨在評估和探討近場音爆水平的預(yù)測技術(shù)[39]。為方便研究不同音爆預(yù)測方法之間的差異，大會給出3個標(biāo)準(zhǔn)模型，分別為NASA圓錐體、三角翼翼身組合體(69° DWB)以及洛克希德·馬丁的超聲速客機模型(LM1021)，其中前2個模型為無升力模型[40]。69° DWB外形是Hunton等[41]研究不同翼面形狀對音爆影響所用的模型之一，圖 12為模型示意圖，模型全長17.52 cm，機頭長7.01 cm，機頭母線沿飛機軸向符合二次函數(shù)分布，機身最大半徑為0.54 cm，三角翼的前緣后掠角為69°，三角翼翼型為菱形翼型，其最大厚度為弦長(c)的5%。

圖12 69°后掠三角翼翼身組合體標(biāo)準(zhǔn)模型示意圖Fig.12 Sketch of 69° sweepback delta wing-body configuration standard model

大會給出的69°后掠三角翼翼身組合體標(biāo)模的狀態(tài)為Ma=1.7,迎角α=0°，來流雷諾數(shù)Re=2.43×106，采用Whitham理論對模型中心線正下方h=62.99 cm處的壓強信號進行計算，計算結(jié)果與試驗值[42]的對比如圖 13所示。

圖13 翼身組合體近場壓強信號計算驗證Fig.13 Calculation verification of near-field pressure signatures for wing-body configuration

3.2.1 外形參數(shù)化與優(yōu)化數(shù)學(xué)模型

機翼參數(shù)化如圖 14所示。為保證三角翼和菱形翼型特征，機翼控制參數(shù)為：機翼1/2弦線后掠角Λ1/2、翼根翼型的弦長c、翼根翼型最大厚度與弦長的比值t、機翼半展長b和機翼位置lw共5個參數(shù)。

機身參數(shù)化如圖 15所示，通過對機身母線的控制來實現(xiàn)。機身母線采用Bezier曲線進行參數(shù)化，共選取10個控制點，坐標(biāo)為(xi,Bi)。xi為第i個控制點的軸向站位(圖15(a)),在母線參數(shù)化和低音爆優(yōu)化過程中保持不變。Bi為第i個控制點的縱坐標(biāo)(圖15(b))，是母線參數(shù)化和低音爆優(yōu)化的變量。

采用Whitham理論計算遠(yuǎn)場音爆強度，優(yōu)化時目標(biāo)為Ma=1.7狀態(tài)下飛機軸線正下方1 000倍機長處N型波的壓強峰值，由于模型為對稱模型，因此在優(yōu)化過程中不考慮升力約束。相應(yīng)優(yōu)化問題的數(shù)學(xué)模型如下

(10)

圖14 翼身組合體的機翼參數(shù)化示意圖Fig.14 Sketch of wing parameterization of wing-body configuration

圖15 翼身組合體的機身母線參數(shù)化示意圖Fig.15 Sketch of fuselage generatrix parameterization of wing-body configuration

式中：Vf、Vf,0分別為外形優(yōu)化后機身容積和基準(zhǔn)機身容積；Vw、Vw,0分別為外形優(yōu)化后機翼容積和基準(zhǔn)機翼容積；lf,0為基準(zhǔn)機身長度；Sw、Sw,0分別為優(yōu)化后機翼面積和基準(zhǔn)機翼面積。

控制翼身組合體外形的15個設(shè)計變量的取值范圍如表 5所示?？刂茩C翼的5個參數(shù)的取值范圍為基準(zhǔn)參數(shù)的(1±0.2)倍?？刂茩C身的10個參數(shù)的取值范圍為基準(zhǔn)參數(shù)的(1±0.5)倍，圖15(b)中控制站位處箭頭長短代表了相應(yīng)控制點可移動的范圍大小，圖 16給出了機身母線能夠變化的范圍。

3.2.2 低音爆優(yōu)化結(jié)果

表5 翼身組合體設(shè)計變量及設(shè)計空間Table 5 Design variables and design space of wing-body configuration

圖16 翼身組合體機身母線變化的范圍以及范圍內(nèi)的兩條曲線Fig.16 Variation range of fuselage generatrix of wing-body configuration and two curves within the range

圖19和圖 20分別給出優(yōu)化前后的外形對比和優(yōu)化前后的機身母線對比。可以看出，優(yōu)化后機翼的后掠角變大，翼展略有減小，在機身與機翼結(jié)合部位，機身的半徑變小，同時為保證機身容積約束，機頭處半徑變大。這樣的布局形式可以使飛機在巡航狀態(tài)(Ma=1.7)下，具有低音爆的有效截面積分布。圖 21給出了優(yōu)化前后的有效截面積(Seff)分布，可知優(yōu)化外形的有效截面積分布在機頭處斜率增大，隨后變化平緩且最大有效截面積減小，使得近場壓強信號表現(xiàn)為第1道激波略有增強，第2道激波明顯減弱，且兩道激波之間距離增大，圖 22為優(yōu)化后的近場(h=62.99 cm)壓強信號。當(dāng)這樣的近場壓強信號傳播到遠(yuǎn)場時，第1道激波與第2道激波合并，使N型波的峰值減小，從而達(dá)到降低音爆強度的目的。

圖17 翼身組合體低音爆優(yōu)化的收斂歷程Fig.17 Convergence history of low boom optimization for wing-body configuration

圖18 翼身組合體優(yōu)化前后遠(yuǎn)場N型波對比Fig.18 Comparison of far-field N-wave between baseline and optimal wing-body configurations

表6 翼身組合體優(yōu)化前后參數(shù)對比Table 6 Comparison of parameters between baseline and optimal wing-body configurations

圖19 翼身組合體優(yōu)化前后外形對比Fig.19 Comparison of geometries between baseline and optimal wing-body configurations

圖20 翼身組合體優(yōu)化前后的機身母線對比Fig.20 Comparison of fuselage generatrix between baseline and optimal wing-body configurations

圖21 翼身組合體優(yōu)化前后有效截面積分布對比Fig.21 Comparison of effective cross-section distribution between baseline and optimal wing-body configurations

圖22 翼身組合體優(yōu)化前后近場壓強信號對比Fig.22 Comparison of near-field pressure signature between baseline and optimal wing-body configurations

4 結(jié) 論

本文將代理優(yōu)化應(yīng)用于飛行器的低音爆優(yōu)化設(shè)計，發(fā)展了一種基于代理模型的低音爆優(yōu)化設(shè)計方法。

1) 用本文方法進行了NASA圓錐體模型的低音爆優(yōu)化設(shè)計，3種優(yōu)化設(shè)計狀態(tài)下，遠(yuǎn)場音爆強度相比于基準(zhǔn)外形降低了近1 dB，表明該方法應(yīng)用于低音爆優(yōu)化設(shè)計的可行性。

2) 針對NASA圓錐體模型的低音爆優(yōu)化設(shè)計，代理優(yōu)化方法相比于梯度優(yōu)化和遺傳算法，具有很好的優(yōu)越性。與梯度優(yōu)化結(jié)果相比具有全局性，與遺傳算法相比具有高效性。

3) 將代理優(yōu)化方法應(yīng)用于翼身組合體外形的低音爆優(yōu)化設(shè)計，優(yōu)化后的模型更具有低音爆的特征，遠(yuǎn)場N型波峰值減少了27.4%，表明本文方法在復(fù)雜外形低音爆優(yōu)化設(shè)計中具有很大的應(yīng)用潛力。

本文的重點是展示基于代理模型的全局優(yōu)化方法在低音爆優(yōu)化設(shè)計中的適用性，并為進一步研究打下基礎(chǔ)。一方面，本文選擇了計算效率高、計算精度合理的Whitham方法來預(yù)測遠(yuǎn)場音爆。這種方法對于復(fù)雜外形的適應(yīng)性和音爆預(yù)測精度都存在一定的不足，下一步擬發(fā)展耦合高可信度CFD和非線性Burgers方程的高精度音爆預(yù)測方法。另一方面，本文優(yōu)化沒有考慮氣動性能的目標(biāo)和約束。為了增強設(shè)計結(jié)果的工程適用性，下一步需要在優(yōu)化中考慮氣動性能的目標(biāo)(如升阻比或阻力)和約束(如升力和力矩)。

參 考 文 獻

[1] 朱自強, 蘭世隆. 超聲速民機和降低音爆研究[J]. 航空學(xué)報, 2015, 36(8): 2507-2528.

ZHU Z Q, LAN S L. Study of supersonic commercial transport and reduction of sonic boom[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(8): 2507-2528 (in Chinese).

[2] WHITHAM G. The flow pattern of a supersonic project[J]. Communications on Pure and Applied Mathematics, 1952, 5(3): 301-347.

[3] HAYES W D. Brief review of basic theory: Sonic boom research: NASA SP-147[R]. Washington, D.C.: NASA, 1967.

[4] SEEBASS R. Sonic boom theory[J]. Journal of Aircraft, 1969, 6(3):177-184.

[5] SEEBASS R, GEORGE A R. Sonic boom minimization[J]. Journal of the Acoustical Society of America, 1972, 51(2): 686-694.

[6] SEEBASS R, ARGROWB. Sonic boom minimization revisited: AIAA-1998-2956[R]. Reston, VA: AIAA, 1998.

[7] CARLSON H W. Simplified sonic boom prediction: NASA TP-1122[R]. Washington, D.C.: NASA, 1978.

[8] THOMAS L C. Extrapolation of sonic boom pressure signatures by the waveform parameter method: NASA TN D-6832[R]. Washington, D.C.: NASA, 1972.

[9] CLEVELAND O R. Propagation of sonic booms through a real, stratified atmosphere[D]. Austin: University of Texas at Austin, 1995: 6-37.

[10] POTAPKIN A V, KOROTAEVA T A, MOSKVICHEV D Y, et al. Anadvanced approach for far-field sonic boom prediction: AIAA-2009-1056[R]. Reston, VA: AIAA, 2009.

[11] YAMASHITA R, SUZUKI K. Full-field sonic boom simulation in real atmosphere: AIAA-2014-2269[R]. Reston, VA: AIAA, 2014.

[12] WINTZER M, NEMEC M, AFTOSMIS J M. Adjoint-based adaptive mesh refinement for sonic boom prediction: AIAA-2008-6593[R]. Reston, VA: AIAA, 2008.

[13] RALLABHANDI S K. Advanced sonic boom prediction using the augmented Burgers equation[J]. Journal of Aircraft, 2011, 48(4): 1245-1253.

[14] 陳鵬, 李曉東. 基于Khokhlov-Zabolotskaya-Kuznetsov方程的聲爆頻域預(yù)測法[J]. 航空動力學(xué)報, 2010, 25(2): 359-365.

CHEN P, LI X D. Frequency domain method for predicting sonic boom propagation based on Khokhlov-Zabolotskaya-Kuznetsov equation[J]. Journal of Aerospace Power, 2010, 25(2): 359-365 (in Chinese).

[15] 但聃, 楊偉. 超音速公務(wù)機聲爆計算與布局討論[J]. 航空工程進展, 2012, 3(1): 8-15.

DAN D, YANG W. Supersonic business jet sonic boom computation and layout discussion[J]. Advances in Aeronautical Science and Engineering, 2012, 3(1): 8-15 (in Chinese).

[16] 馮曉強. 聲爆計算方法研究及在超聲速客機設(shè)計的應(yīng)用[D]. 西安: 西北工業(yè)大學(xué), 2012: 9-24.

FENG X Q. The research of sonic boom prediction method and application in supersonic aircraft design[D]. Xi’an: Northwestern Polytechnical University, 2012: 9-24 (in Chinese).

[17] 馮曉強, 宋筆鋒, 李占科, 等. 超聲速飛機低聲爆布局混合優(yōu)化方法研究[J]. 航空學(xué)報, 2013, 34(8): 1768-1777.

FENG X Q, SONG B F, LI Z K, et al. Hybrid optimization approach research for low sonic boom supersonic aircraft configuration[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(8): 1768-1777 (in Chinese).

[18] CHOI S, ALONSO J J, KROO M I. Multi-fidelity design optimization of low-boom supersonic business jets: AIAA-2004-4371[R]. Reston, VA: AIAA, 2004.

[19] CHAN K M. Supersonicaircraft optimization for minimizing drag and sonic boom[D]. Palo Alto: Stanford University, 2003: 48-96.

[20] FARHAT C, MAUTE K, ARGROW B, et al. Shape optimization methodology for reducing the sonic boom initial pressure rise[J]. AIAA Journal, 2007, 45(5): 1007-1018.

[21] AFTOSMIS M J, NEMEC M, CLIFF S E. Adjoint-based low-boom design with Cart3D: AIAA-2011-3500[R]. Reston, VA: AIAA, 2011.

[22] WINTZER M, KROO I. Optimization and adjoint-based CFD for the conceptual design of low sonic boom aircraft: AIAA-2012-0963[R]. Reston, VA: AIAA, 2012.

[23] RALLABHANDI K S. Sonic boom adjoint methodology and its applications: AIAA-2011-3497[R]. Reston, VA: AIAA, 2011.

[24] RALLABHANDI K S. Application of adjoint methodology to supersonic aircraft design using reversed equivalent areas[J]. Journal of Aircraft, 2014, 51(6): 1873-1882.

[25] MINELLI A, SALAH E D I, CARRIER G. Advancedoptimization approach for supersonic low-boom design: AIAA-2012-2168[R]. Reston, VA: AIAA, 2012.

[26] LI J H, TIM W, RAMESH K A. Shapeoptimization of supersonic bodies to reduce sonic boom signature: AIAA-2016-3432[R]. Reston, VA: AIAA, 2016.

[27] LI J H, KRAMPF J, MITCHELL J, et al. Prediction,minimization and propagation of sonic boom from supersonic bodies: AIAA-2017-0277[R]. Reston, VA: AIAA, 2017.

[28] CARLSON H W, MACK R J, MORRIS O A. A wind tunnel investigation of the effect of body shape on sonic boom pressure distributions: NASA TN D-3106[R]. Washington, D.C.: NASA, 1965.

[29] 韓忠華. Kriging模型及代理優(yōu)化算法研究新進展[J]. 航空學(xué)報, 2016, 37(11): 3197-3225.

HAN Z H. Kriging surrogate model and its application to design optimization: A review of recent progress[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(11): 3197-3225 (in Chinese).

[30] GIUNTA A A, WOJTKIEWICZ J S F, ELDRED M S. Overview of modern design of experiments methods for computational simulations: AIAA-2003-649[R]. Reston, VA: AIAA, 2003.

[31] KRIGE D G. A statistical approach to some basic mine valuation problems on the Witwatersrand[J]. Journal of the Chemical Metallurgical & Mining Society of South Africa, 1951, 52(6): 119-139.

[32] HAN Z H, GOERTZ S. Hierarchical kriging model for variable-fidelity surrogate modeling[J]. AIAA Journal, 2012, 50(5): 1285-1296.

[33] HAN Z H, ZHANG Y, SONG C X, et al. Weighted gradient-enhanced kriging for high-dimensional surrogate modeling and design optimization[J]. AIAA Journal, 2017, 55(12): 4330-4346.

[34] FORRESTER A I J, KEANE A J. Recentadvances in surrogate-based optimization[J]. Progress in Aerospace Sciences, 2009, 45(1): 50-79.

[35] JONES D R. Ataxonomy of global optimization methods based on response surfaces[J]. Journal of Global Optimization, 2001, 21(4): 345-383.

[36] LIU J, SONG W P, HAN Z H, et al. Efficient aerodynamic shape optimization of transonic wings using a parallel infilling strategy and surrogate models[J]. Structural and Multidisciplinary Optimization, 2016, 55(3):925-943.

[37] HAN Z H. SurroOpt a generic surrogate-based optimization code for aerodynamic and multidisciplinary design: ICAS-2016-0281[R]. Bonn: ICAS, 2016.

[38] GILL P E, MURRAY W. NPSOL: AFORTRAN package for nonlinear programming[EB/OL]. (2013-07-15) [2016-07-20]. http:∥www.sbsi-sol-optimize.com.

[39] AFTOSMIS M J, NEMEC M. Cart3D simulations for the first AIAA sonic boom prediction workshop: AIAA-2014-0558[R]. Reston, VA: AIAA, 2014.

[40] DAGRAU F, LOSEILLE A, SALAH E D I. Computational and experimental assessment of models for the first AIAA sonic boom prediction workshop using high fidelity CFD methods with mesh adaptation: AIAA-2014-2009[R]. Reston, VA: AIAA,2014.

[41] HUNTON W L, HICKS M R, MENDOZA P J. Some effects of wing planform on sonic boom: NASA TN D-7160[R]. Washington, D.C.: NASA, 1973.

[42] MA P B, WANG G, REN J, et al. Near field sonic boom analysis with HUNS3D solver: AIAA-2017-0038[R]. Reston, VA, AIAA, 2017.