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

西北工業(yè)大學(xué) 航空學(xué)院 翼型葉柵空氣動(dòng)力學(xué)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室，西安 710072

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

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

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

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

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

1 音爆預(yù)測(cè)方法

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

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

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

(1)

(2)

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

近場(chǎng)壓強(qiáng)信號(hào)在向遠(yuǎn)場(chǎng)傳播的過(guò)程中，由于激波的匯聚，擾動(dòng)壓強(qiáng)在一定條件下會(huì)在遠(yuǎn)場(chǎng)形成N型波。遠(yuǎn)場(chǎng)N型波強(qiáng)度的計(jì)算公式為

(3)

式(3)為自由空間中過(guò)壓值的計(jì)算式，當(dāng)觀測(cè)點(diǎn)位于地面時(shí)，考慮地面反射效應(yīng)因子KR，有

(4)

1.2 算例驗(yàn)證

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

r=

(5)

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

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

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

圖2 NASA圓錐體模型近場(chǎng)壓強(qiáng)信號(hào)試驗(yàn)觀測(cè)位置Fig.2 Experimentally observed position of near-field pressure signature for NASA cone model

圖3 近場(chǎng)壓強(qiáng)信號(hào)計(jì)算值與試驗(yàn)值的對(duì)比Fig.3 Comparison between calculation and experimental data of near-field pressure signature

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

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

2.1 試驗(yàn)設(shè)計(jì)

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

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

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

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

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

(6)

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

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

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

(7)

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

優(yōu)化加點(diǎn)準(zhǔn)則是指根據(jù)所建立的代理模型增加新樣本點(diǎn)的準(zhǔn)則。加點(diǎn)的目的是擴(kuò)充樣本集合，提高代理模型的預(yù)測(cè)精度。代理優(yōu)化的本質(zhì)是采用某種加點(diǎn)準(zhǔn)則產(chǎn)生新樣本點(diǎn)，并使這些樣本點(diǎn)逐步落入最優(yōu)點(diǎn)或其附近。目前已經(jīng)有多種加點(diǎn)準(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”的組合加點(diǎn)準(zhǔn)則。

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

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

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

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

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

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

本節(jié)將代理優(yōu)化方法應(yīng)用于NASA圓錐體模型和2014年音爆預(yù)測(cè)大會(huì)翼身組合體標(biāo)準(zhǔn)模型的低音爆優(yōu)化設(shè)計(jì)。以NASA圓錐體優(yōu)化算例來(lái)驗(yàn)證代理優(yōu)化應(yīng)用的適用性和優(yōu)越性，以翼身組合體優(yōu)化算例來(lái)說(shuō)明代理優(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)外形進(jìn)行低音爆優(yōu)化設(shè)計(jì)。外形參數(shù)化及優(yōu)化問(wèn)題的數(shù)學(xué)模型建立過(guò)程如下。

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

(8)

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

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

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

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è)計(jì)的目標(biāo)函數(shù)為遠(yuǎn)場(chǎng)(h/l=1 000，見(jiàn)圖 6)音爆N型波的峰值，3種設(shè)計(jì)狀態(tài)的來(lái)流馬赫數(shù)分別為1.26、1.60和2.00。相應(yīng)的數(shù)學(xué)模型為

(9)

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

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

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

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

圖7 3種設(shè)計(jì)狀態(tài)下NASA圓錐體模型低音爆優(yōu)化設(shè)計(jì)中目標(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è)計(jì)變量及遠(yuǎn)場(chǎng)聲壓級(jí)對(duì)比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時(shí)，表示該約束為主動(dòng)約束。

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

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

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

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

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

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

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

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

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

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

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

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

圖13 翼身組合體近場(chǎng)壓強(qiáng)信號(hào)計(jì)算驗(yàn)證Fig.13 Calculation verification of near-field pressure signatures for wing-body configuration

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

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

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

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

(10)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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

參 考 文 獻(xiàn)

[1] 朱自強(qiáng), 蘭世隆. 超聲速民機(jī)和降低音爆研究[J]. 航空學(xué)報(bào), 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ù)測(cè)法[J]. 航空動(dòng)力學(xué)報(bào), 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ī)聲爆計(jì)算與布局討論[J]. 航空工程進(jìn)展, 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] 馮曉強(qiáng). 聲爆計(jì)算方法研究及在超聲速客機(jī)設(shè)計(jì)的應(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] 馮曉強(qiáng), 宋筆鋒, 李占科, 等. 超聲速飛機(jī)低聲爆布局混合優(yōu)化方法研究[J]. 航空學(xué)報(bào), 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ìn)展[J]. 航空學(xué)報(bào), 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.