Kriging模型在翼型反設(shè)計中的應(yīng)用研究

劉 俊，宋文萍，韓忠華，許建華，樊艷紅

（西北工業(yè)大學(xué) 翼型葉柵空氣動力學(xué)國防科技重點實驗室，陜西 西安 710072）

0 引 言

傳統(tǒng)的基于CFD數(shù)值模擬的氣動優(yōu)化設(shè)計方法主要分為梯度法和非梯度法。非梯度方法使用遺傳算法（GA）等直接調(diào)用流場分析程序進行優(yōu)化，它們往往具有良好的全局性，但計算量較大?；谔荻鹊姆椒ǖ年P(guān)鍵在于梯度信息的計算，傳統(tǒng)的有限差分法雖然簡單，但在設(shè)計變量較多時，計算量較大。

從20世紀80年代末開始，由Jameson［1］等人發(fā)展起來的基于控制理論的優(yōu)化方法（Adjoint方法）以其對多變量問題的適用性和高效的特點迅速成為氣動設(shè)計的主流方法，并已大量應(yīng)用于氣動優(yōu)化設(shè)計和反設(shè)計中［2－4］。直至當(dāng)今，該方法仍然在氣動設(shè)計中占據(jù)重要地位［5－6］。

雖然Adjoint方法具有高效的特點，但其仍然屬于局部優(yōu)化方法的范疇，且該方法通用性不強，針對不同類型的優(yōu)化問題，需要重新推導(dǎo)伴隨方程，重新編制伴隨方程求解程序；在特殊情況下，如考慮自由轉(zhuǎn)捩的NS方程求解器，Adjoint方法應(yīng)用起來有難度。為了兼顧全局性與高效性，同時考慮到通用性的要求，基于代理模型的優(yōu)化方法逐步流行起來，并且越來越受到重視。在眾多的代理模型優(yōu)化方法中，基于多項式響應(yīng)面的方法（P－RSM）最為成熟，并已成功應(yīng)用于氣動優(yōu)化設(shè)計與反設(shè)計中［7－8］。

近十年來，由于Kriging模型具有可以模擬復(fù)雜響應(yīng)并可提供誤差信息的優(yōu)勢，基于Kriging模型的優(yōu)化方法［9］被受到重視并越來越多地應(yīng)用于氣動及多學(xué)科優(yōu)化設(shè)計中［10－17］，以大大減少優(yōu)化設(shè)計所需的計算時間并得到優(yōu)良的設(shè)計結(jié)果。

盡管基于Kriging模型的優(yōu)化方法在氣動優(yōu)化設(shè)計中取得了很大的成功，但一直未在反設(shè)計中得到應(yīng)用。雖然Toal等［18－19］以翼型反設(shè)計為算例，比較了使用不同優(yōu)化算法來優(yōu)化Kriging超參數(shù)對最終反設(shè)計結(jié)果的影響，然而他們的出發(fā)點并非基于Kriging模型的優(yōu)化方法在翼型反設(shè)計中的應(yīng)用，也未得到真正有效的反設(shè)計結(jié)果。其原因可能是使用了不恰當(dāng)?shù)囊硇蛥?shù)化方法，且使用的樣本點加點方法不當(dāng)。可能正是受他們的影響，使得基于Kriging模型的優(yōu)化方法未能應(yīng)用于反設(shè)計中。因此，研究Kriging模型在翼型反設(shè)計中的可行性，發(fā)展有效、實用的基于Kriging模型的翼型反設(shè)計方法是十分必要的。

本文利用 Kriging模型，結(jié)合EI方法［9］及直接優(yōu)化代理模型最小值等樣本點加點準(zhǔn)則、優(yōu)化算法、CST參數(shù)化方法，進行了翼型單目標(biāo)、多目標(biāo)反設(shè)計。本文方法中，按照模型加點準(zhǔn)則又分為三種：單獨使用EI方法（SIC1）、單獨使用直接優(yōu)化代理模型最小值（SIC2）、同時使用上述兩種加點準(zhǔn)則（SIC3）。

首先在本文優(yōu)化方法中分別使用三種不同加點方法進行了翼型單目標(biāo)反設(shè)計，驗證了本文方法的可行性并比較了不同加點方法對反設(shè)計結(jié)果的影響，同時與傳統(tǒng)的基于二次響應(yīng)面模型的方法及Adjoint方法進行了比較，結(jié)果表明，分別使用三種加點方法的基于Kriging的方法均適用于翼型反設(shè)計，且都優(yōu)于二次響應(yīng)面方法，且使用SIC2效果最好；與Adjoint方法相比，反設(shè)計效果及計算效率相當(dāng)。

其次，在本文優(yōu)化方法中應(yīng)用SIC2加點方法進行了翼型多目標(biāo)反設(shè)計，驗證了本文方法在翼型多目標(biāo)反設(shè)計中的適用性。

最后，將本文方法應(yīng)用于接近工程實際的在已有翼型基礎(chǔ)上修改壓力分布進行了單目標(biāo)、對目標(biāo)反設(shè)計，驗證了本文方法的可行性。

1 Kriging模型

Kriging模型是由南非地質(zhì)學(xué)家Krige［20］提出來的一種統(tǒng)計學(xué)插值模型，并由Sacks［21］最早應(yīng)用于計算機實驗設(shè)計與分析，隨后作為一種代理模型被廣泛用于數(shù)值模擬分析與優(yōu)化中響應(yīng)值的預(yù)測。本文采用的是普通Kriging模型。

1.1 Kriging預(yù)測值與均方誤差

Kriging模型用一個常數(shù)和一個隨機過程的和來表示系統(tǒng)的響應(yīng)值：

隨機過程Z（·）的均值為0，協(xié)方差如下：

其中σ2是Z（·）的方差（假定對任意的x有（x）＝σ2（x）≡σ2），R是只與兩點x和x′的距離有關(guān)的相關(guān)函數(shù)。

假定設(shè)計空間內(nèi)任意一點的響應(yīng)值可由已知樣本點的響應(yīng)值ys的線性組合來近似，則在未知點x，Kriging近似值有如下形式：

其中，w＝（w（1），...，w（ns））T為權(quán)系數(shù)。

用YS＝（Y（1），...，Y（ns））T替換上式中的yS＝（y（1），...，y（ns））T，最小化均方誤 差 （Mean Squared Error）可得到Kriging的預(yù)測值如下：

其中，1是單位列向量，

任意未知點x處的Kriging預(yù)測值的均方誤差（MSE）為：

其中，是方差σ2的估計值：

1.2 相關(guān)模型

相關(guān)矩陣R和相關(guān)向量r的建立與相關(guān)函數(shù)的選取有關(guān)，假定相關(guān)函數(shù)的函數(shù)值只與兩點x（i）、x（j）之間的空間距離有關(guān)。相關(guān)模型的選取應(yīng)該具有以下特征：a）當(dāng)距離趨于0時，函數(shù)值趨于1；b）當(dāng)距離增加時，函數(shù)值光滑的減?。籧）當(dāng)距離趨于無窮時，函數(shù)值趨于0；d）至少一階可導(dǎo)。此處只考慮以下形式的相關(guān)模型：

本文采用三次樣條函數(shù)作為相關(guān)函數(shù)：

其中

1.3 Kriging參數(shù)的獲得

式（9）中的 Kriging超參數(shù)θ（hyper parameter）可以通過求解下列最大似然估計（MLE）這一優(yōu)化問題獲得：

2 加點準(zhǔn)則

盡管Kriging模型具有較好的全局近似能力，但是僅依靠初始樣本信息并不能找到設(shè)計空間內(nèi)的最優(yōu)值，因而需要通過增加新的樣本點來尋找真實最優(yōu)值并提高Kriging模型精度。本文采用了以下兩種加點準(zhǔn)則［22－23］。

2.1 最大化Expected Improvement（EI）

設(shè)樣本集中的最優(yōu)目標(biāo)函數(shù)值為ymin，假設(shè)隨機變量Y（x）服從均值為（x），標(biāo)準(zhǔn)差為s（x）的正太分布，Y（x）∈N［（x），s2］，其 中（x）為 Kriging 預(yù) 測值，見式（4），s2為均方誤差（MSE），見式（7）。在x處目標(biāo)函數(shù)相對的改進量為I＝y(tǒng)min－Y（x）＞0，于是I的期望值由下式給出：

求出設(shè)計空間內(nèi)EI最大的點作為新的樣本點。

2.2 最小化代理模型 （MP）

用優(yōu)化算法尋找代理模型上目標(biāo)函數(shù)的最小值，將最小值點作為新的樣本點加入樣本集。

3 二次響應(yīng)面模型

對于完全二階多項式模型的一般公式為：

其中k為設(shè)計變量個數(shù)，c0、ci、cij是待定參數(shù)。

對于k元完全二階多項式，待定參數(shù)的個數(shù)為：

利用樣本點的響應(yīng)值可求得所有待定參數(shù)。

4 流場求解

本文采用課題組自主開發(fā)的PMNS2D程序通過求解雷諾平均Navier－Stockes方程對翼型進行氣動分析，采用中心格式離散，湍流模型采用Spalart－Allmaras模型，并使用了隱式殘值光順、當(dāng)?shù)貢r間步長、多重網(wǎng)格等加速收斂措施。計算網(wǎng)格為C型結(jié)構(gòu)化網(wǎng)格。

5 幾何參數(shù)化

5.1 CST參數(shù)化方法

CST參數(shù)化方法是由波音公司的工程師Kulfan［24］提出來的。該方法用一個類函數(shù)C（x／c）和一個型函數(shù)S（x／c）的乘積加上一個描述后緣厚度的函數(shù)表示一個翼型的幾何形狀：

其中c為翼型弦長，類函數(shù)C（x／c）具有以下的一般形式：

型函數(shù)S（x／c）由伯恩斯坦多項式各項的加權(quán)求和得到：其中Ai是加權(quán)系數(shù)，Kr，n為多項式的系數(shù)：

5.2 設(shè)計變量范圍的確定

利用一個較“厚”的翼型和一個較“薄”的翼型，并使得這兩個翼型之間的范圍足夠大，可以包含目標(biāo)翼型，將兩者之間作為真實的設(shè)計空間，如圖1中陰影部分所示；將這兩個翼型所對應(yīng)的CST參數(shù)作為設(shè)計變量的上下限。

圖1 翼型反設(shè)計的設(shè)計空間示意圖Fig.1 Schematic diagram of the design space for airfoil inverse design

6 翼型反設(shè)計目標(biāo)函數(shù)

以設(shè)計結(jié)果的壓力分布和目標(biāo)壓力分布的差為目標(biāo)函數(shù)，對于第k個設(shè)計狀態(tài)：

其中，Cpi為計算翼型表面第i個點的壓力系數(shù)，＾Cpi為目標(biāo)翼型第i個點的壓力系數(shù)。對于m個目標(biāo)的反設(shè)計，加權(quán)目標(biāo)函數(shù)為：

ωk為第k個目標(biāo)的權(quán)重系數(shù)。對于更重要的目標(biāo)壓力分布，其權(quán)系數(shù)應(yīng)取得更大。

7 代理模型優(yōu)化方法

本文以當(dāng)前翼型的壓力分布和目標(biāo)壓力分布的差為目標(biāo)函數(shù)，將反設(shè)計問題轉(zhuǎn)化為優(yōu)化問題，用基于代理模型的優(yōu)化方法來求解。

基于代理模型的優(yōu)化方法首先采用試驗設(shè)計方法（Design of Experiment，DoE）生成設(shè)計空間內(nèi)的樣本點，并獲得樣本點的響應(yīng)值，由樣本信息建立代理模型，以近似描述空間內(nèi)任意位置的真實響應(yīng)值，采用優(yōu)化算法（本文采用遺傳算法）尋找新的樣本點，將新的樣本點加入原樣本集，重新建立代理模型進行優(yōu)化直至滿足停止條件。

1）基于二次響應(yīng)面模型的方法（P－RSM）

采用二次響應(yīng)面模型，直接尋找二次響應(yīng)面目標(biāo)函數(shù)的最小值點作為新的樣本點。

2）Kriging—SIC1

采用Kriging模型，尋找設(shè)計空間內(nèi)EI最大值點作為新的樣本點，即EI方法。

3）Kriging—SIC2

采用Kriging模型，直接尋找設(shè)計空間內(nèi)目標(biāo)函數(shù)的最小值點作為新的樣本點，即MP方法。

4）Kriging—SIC3

采用Kriging模型，同時使用EI、MP兩種加點方法，一次增加兩個樣本點。

流程圖如圖2所示。

圖2 代理模型優(yōu)化方法流程圖Fig.2 Flowchart of the surrogate－based optimizer

8 方法驗證

8.1 本文方法在反設(shè)計中的適應(yīng)性研究以及與PRSM、Adjoint方法的比較

以跨聲速翼型RAE2822單目標(biāo)反設(shè)計為例進行測試。設(shè)計狀態(tài)如下：

使用19個設(shè)計變量，使用代理模型方法進行反設(shè)計時，取相同的設(shè)計空間。

首先采用P－RSM方法進行反設(shè)計，反設(shè)計的壓力分布及翼型幾何形狀與目標(biāo)壓力分布及目標(biāo)翼型的比較見圖3和圖4，從圖中看出，雖然壓力分布及翼型幾何形狀與目標(biāo)翼型的比較接近，但存在明顯差距。

其次，用Adjoint方法進行反設(shè)計，初始翼型為NACA0012。反設(shè)計結(jié)果見圖5、圖6，從圖中看出，反設(shè)計的壓力分布及翼型幾何形狀與目標(biāo)壓力分布及目標(biāo)翼型均吻合良好。

圖3 P－RSM方法反設(shè)計壓力分布與目標(biāo)壓力分布比較Fig.3 Comparison of the designed pressure distribution by P－RSM and the target one

圖4 P－RSM方法反設(shè)計翼型與目標(biāo)翼型比較Fig.4 Comparison of designed airfoil by P－RSM and the target one

圖5 Adjoint方法反設(shè)計壓力分布與目標(biāo)壓力分布比較Fig.5 Comparison of the designed pressure distribution by Adjoint and the target one

圖6 Adjoin方法反設(shè)計翼型與目標(biāo)翼型比較Fig.6 Comparison of designed airfoil by Adjoint and the target one

再次，使用本文基于Kriging模型的方法分別采用三種加點方法進行反設(shè)計。反設(shè)計的壓力分布及翼型幾何形狀與目標(biāo)壓力分布及目標(biāo)翼型的比較見圖7～圖12，從圖可以看出，分別使用三種加點方法的基于Kriging模型的方法反設(shè)計的翼型壓力分布與幾何形狀與目標(biāo)翼型的均吻合良好，且均優(yōu)于PRSM方法，從而驗證了本文基于Kriging的方法在翼型反設(shè)計中的適應(yīng)性。此外，從三種加點方法看，使用SIC2即僅采用直接優(yōu)化代理模型最小值的加點方法所得的反設(shè)計效果最好。

圖7 采用SIC1的Kriging方法反設(shè)計壓力分布與目標(biāo)壓力分布比較Fig.7 Comparison of the designed pressure distribution by Kriging－SIC1and the target one

圖8 采用SIC1的Kriging方法反設(shè)計翼型與目標(biāo)翼型比較Fig.8 Comparison of the designed airfoil by Kriging－SIC1and the target one

圖9 采用SIC2的Kriging方法反設(shè)計壓力分布與目標(biāo)壓力分布比較Fig.9 Comparison of the designed pressure distribution by Kriging－SIC2and the target one

圖10 采用SIC2的Kriging方法反設(shè)計翼型與目標(biāo)翼型比較Fig.10 Comparison of the designed airfoil by Kriging－SIC2and the target one

圖11 采用SIC3的Kriging方法反設(shè)計壓力分布與目標(biāo)壓力分布比較Fig.11 Comparison of the designed pressure distribution by Kriging－SIC3and the target one

圖12 采用SIC3的Kriging方法反設(shè)計翼型與目標(biāo)翼型比較Fig.12 Comparison of the designed airfoil by Kriging－SIC3and the target one

表1給出了本文方法及P－RSM方法最終的目標(biāo)值的比較。注意此表中的目標(biāo)值是實際目標(biāo)值與Kriging初始樣本點中最大目標(biāo)值的比值。從表中也可以看出，采用SIC2所得的目標(biāo)值最優(yōu)，SIC3次之，SIC1再次，但都優(yōu)于P－RSM方法。

表1 四種方法目標(biāo)值與CFD計算次數(shù)比較Table 1 Comparison of objective value and runs of flow solver for the 4methods

圖13給出了本文方法及Adjoint方法反設(shè)計目標(biāo)函數(shù)的收斂曲線，其中前面的方形符號表示建立Kriging模型初始樣本點的目標(biāo)值。注意，此圖中的使用Kriging模型的目標(biāo)值是真實目標(biāo)值與初始樣本點中最大目標(biāo)值的比值；而Adjoint的目標(biāo)值是實際目標(biāo)值與初始目標(biāo)值的比值，且把求解一次伴隨方程的計算量視為一次流場求解的計算量。從該圖可以看出，本文方法中采用SIC2的目標(biāo)函數(shù)值下降最快，且最終結(jié)果最優(yōu)。此外，使用SIC2的基于Kriging的方法與Adjoint方法相比，收斂速度相當(dāng)且目標(biāo)值下降量級相當(dāng)。

圖13 使用3種加點準(zhǔn)則的Kriging方法及Adjoint方法反設(shè)計目標(biāo)值收斂曲線比較（對于Adjoint方法，一次Adjoint方程求解視為一次流場求解）Fig.13 Comparison of convergence histories for the 3SIC and the Adjoint method（For Adjoint method，the cost of Adjoint solver is treated as equal to the cost of flowsolver）

8.2 翼型多目標(biāo)反設(shè)計

為了驗證本文方法在翼型多目標(biāo)反設(shè)計中的適應(yīng)性，以RAE2822翼型兩個目標(biāo)的反設(shè)計為例進行測試，即設(shè)計一個翼型同時滿足兩個狀態(tài)的壓力分布。設(shè)計狀態(tài)如下：

設(shè)計變量及其范圍與單目標(biāo)算例相同。由8.1可知，使用SIC2所得的反設(shè)計效果最優(yōu)，故此處采用SIC2加點方法。作為算例，此處第一個狀態(tài)的權(quán)重系數(shù)取0.6。

反設(shè)計翼型的壓力分布及幾何形狀見圖14、圖15，從圖看出，壓力分布與及幾何形狀與目標(biāo)壓力分布及幾何形狀都吻合得很好。圖16給出了目標(biāo)值的收斂歷程，其中縱坐標(biāo)為真實目標(biāo)值與初始樣本點中最大目標(biāo)值的比值。

9 工程應(yīng)用的可行性驗證

反設(shè)計方法是工程實際中進行翼型設(shè)計的主要方法之一。翼型反設(shè)計一般首先在已有翼型基礎(chǔ)上修改其壓力分布，再進行反設(shè)計，以獲得更優(yōu)性能的翼型。為了驗證本文方法在工程應(yīng)用中的可行性，本節(jié)先對已有翼型的壓力分布進行修改，再進行反設(shè)計。

圖14 基于Kriging的方法反設(shè)計壓力分布與目標(biāo)壓力分布比較Fig.14 Comparison the designed pressure distributions and the target ones

圖15 基于Kriging的方法反設(shè)計翼型與目標(biāo)翼型比較Fig.15 Comparison of the designed airfoil by Kriging－SIC2and the target one

圖16 翼型多目標(biāo)反設(shè)計目標(biāo)值收斂歷程Fig.16 Convergence history of the multiobjective inverse design

9.1 單目標(biāo)反設(shè)計

以RAE2822為初始翼型，其在狀態(tài)Ma＝0.73，Re＝6.5e6，α＝2.3°下，翼型上表面中部存在強激波，故其阻力較大，此處將上表面壓力分布加以修改，以期在同樣的流動狀態(tài)下無激波，從而大大降低阻力。在反設(shè)計過程中，只對翼型上表面進行參數(shù)化，下表面外形保持不變。在修改上表面壓力分布時，前緣附近（x／c＜0.1）通過手動調(diào)節(jié)，之后通過兩段直線連接，且使其具有較大范圍的超聲速區(qū)并保證升力系數(shù)不變。

反設(shè)計所得壓力分布、修改后的壓力分布（即目標(biāo)壓力分布）以及原壓力分布的對比如圖17所示。從該圖看出，反設(shè)計壓力分布與目標(biāo)壓力分布吻合地很好，這是由于此處給出的目標(biāo)壓力分布比較合理；若目標(biāo)壓力分布不合理，則可能無法獲得與之吻合很好的翼型，同時也說明當(dāng)目標(biāo)壓力分布合理時，本文方法能獲得滿足目標(biāo)壓力分布的翼型。

圖17 設(shè)計壓力分布、修改后壓力分布及原始壓力分布比較Fig.17 Comparison of the designed、modified and initial pressure distributions

9.2 翼型多點反設(shè)計

仍以RAE2822為初始翼型，取兩個設(shè)計狀態(tài)

其中，第一個狀態(tài)與上例中的相同，并將其定為主設(shè)計點。第二個狀態(tài)為非設(shè)計點，希望在保證第一設(shè)計狀態(tài)性能的同時降低該狀態(tài)的阻力，從而改善翼型的阻力發(fā)散特性。與上例相同，僅對翼型上表面進行參數(shù)化，第一個狀態(tài)的目標(biāo)壓力分布同上例。在修改第二狀態(tài)的壓力分布時，保證了升力系數(shù)不變，并假設(shè)翼型中部存在弱激波。

此處使用兩組加權(quán)系數(shù)（0.6，0.4）、（0.4，0.6）分別進行反設(shè)計。采用兩組加權(quán)系數(shù)進行反設(shè)計所得的壓力分布與目標(biāo)壓力分布及原壓力分布的對比分別見圖18和圖19?？梢钥闯觯瑑牲c設(shè)計所得的壓力分布與目標(biāo)壓力分布差別較大，這是由于給出的目標(biāo)壓力分布并不一定合理，同時滿足兩個目標(biāo)壓力分布的翼型并不存在。此外，對比圖18和圖19，當(dāng)?shù)谝粋€狀態(tài)的權(quán)重系數(shù)減小時，設(shè)計的翼型在第一個狀態(tài)的壓力分布與目標(biāo)壓力分布吻合度降低，而在第二個狀態(tài)的壓力分布與目標(biāo)壓力分布吻合度提高。

圖18 權(quán)重系數(shù)?。?.6，0.4）時，設(shè)計壓力分布、修改后壓力分布及原始壓力分布比較Fig.18 Comparison of the designed、modified and initial pressure distributions

圖19 權(quán)重系數(shù)取（0.4，0.6）時，設(shè)計壓力分布、修改后壓力分布及原始壓力分布比較Fig.19 Comparison of the designed、modified and initial pressure distributions

圖20給出了原始翼型、單目標(biāo)設(shè)計、多目標(biāo)設(shè)計翼型在迎角2.3°下的阻力系數(shù)隨馬赫數(shù)變化曲線，從該圖看出，單點設(shè)計的阻力特性明顯優(yōu)于原翼型，而多點設(shè)計盡管沒能很好吻合目標(biāo)壓力分布，但其阻力特性仍優(yōu)于單點設(shè)計。

圖20 阻力系數(shù)隨馬赫數(shù)變化對比Fig.20 Comparison of the drag coefficients at different Mach numbers

從前面兩個算例看出，進行翼型反設(shè)計，首先需要給出合理的目標(biāo)壓力分布，尤其是多目標(biāo)反設(shè)計。實際應(yīng)用中，往往需要在反設(shè)計與修改壓力分布之間進行多次迭代。

10 結(jié)論與分析

（1）在基于Kriging模型的氣動優(yōu)化設(shè)計中，通常使用EI方法作為尋優(yōu)及提高模型精度的準(zhǔn)則。而當(dāng)應(yīng)用于氣動反設(shè)計時，應(yīng)當(dāng)使用MP加點準(zhǔn)則，即直接優(yōu)化目標(biāo)函數(shù)最小值的方法效果更優(yōu)。

（2）在進行翼型反設(shè)計時，Kriging模型明顯優(yōu)于多項式響應(yīng)面模型。

（3）本文方法與Adjoint方法相比所需流場分析次數(shù)相當(dāng)，但是其通用性好，可以方便的更換流場分析程序。

（4）Kriging模型不僅適用于氣動優(yōu)化設(shè)計，同時也適用于氣動反設(shè)計。而合理的目標(biāo)壓力分布是反設(shè)計的關(guān)鍵，實際應(yīng)用中往往需要在反設(shè)計與修改壓力分布之間進行多次迭代。

致謝：感謝王樂提供Adjoint翼型反設(shè)計的結(jié)果。

［1］JAMESON A.Aerodynamic design via control theory［J］.Journal ofScientificComputation，1988，3（3）：233－260.

［2］JAMESON A.Control theory based airfoil design using the Euler equations［R］.AIAA Paper 94－4272.

［3］JAMESON A，PIERCE N A，MARTINELLI L.Optimum aerodynamic design using the Navier－Stockes equations［R］.AIAA Paper 97－0101.

［4］YANG X D，QIAO Z D，ZHU B.Control theory based aerodynamic shape inverse design for subsonic and transonic wings［J］.ACTAAerodynamicaSinica，2003，21（1）：12－19.（in Chinese）楊旭東，喬志德，朱兵.亞、跨音速三維機翼氣動外形反設(shè)計的控制理論方法［J］.空氣動力學(xué)學(xué)報，2003，21（1）：12－19.

［5］WALTHER B，NADARAJAH S.Constrained adjoint－based aerodynamic shape optimization in a multistage turbomachinery environment［R］.AIAA Paper 2012－0062.

［6］WINTZER M，KROO I.Optimization and adjoint－based CFD for the conceptual design of low sonic boom aircraft［R］.AIAA Paper 2012－0963.

［7］XIONG J T，QIAO Z D，HAN Z H.Response surface based aerodynamic design of transonic wings［J］.ACTAAeronaotica etAstronauticaSinica，2006，27（3）：399－402.（in Chinese）熊俊濤，喬志德，韓忠華.基于響應(yīng)面法的跨聲速機翼氣動優(yōu)化設(shè)計［J］.航空學(xué)報，2006，27（3）：399－402.

［8］DENG L，QIAO Z D，XIONG J T.et al.A multi－objective inverse design method for natural laminar airfoil［J］.ACTA AeronaoticaetAstronauticaSinica，2010，31（7）：1373－1378.（in Chinese）鄧磊，喬志德，熊俊濤，等.多目標(biāo)自然層流翼型反設(shè)計方法［J］.航空學(xué)報，2010，31（7）：1373－1378.

［9］JONES D R，SCHONLAU M，WELCH W J.Efficient global optimization of expensive black－box functions［J］.Journalof GlobalOptimization，1998，13：455－492.

［10］XU R F，SONG W P，HAN Z H.Study on the improving Kriging－based airfoil optimization design method［J］.Journalof NorthwesternPolytechnicalUniversity，2010，28（4）：503－510.（in Chinese）許瑞飛，宋文萍，韓忠華.改進Kriging模型在翼型氣動優(yōu)化設(shè)計中的應(yīng)用研究［J］.西北工業(yè)大學(xué)學(xué)報，2010，28（4）：503－510.

［11］SU W，BAI J Q.A surrogate－based optimization algorithm and its application to aerodynamic optimization design［J］.Journal ofProjectiles，Rockets，MissilesandGuidance，2008，28（3）：199－202.（in Chinese）蘇偉，白俊強.一種代理模型方法及其在氣動優(yōu)化設(shè)計中的應(yīng)用［J］.彈箭與制導(dǎo)學(xué)報，2008，28（3）：199－202.

［12］WANG X F，XI G.Aerodynamic optimization design for airfoil based on Kriging model［J］.ACTAAeronaoticaetAstronautica Sinica，2005，26（5）：545－549.（in Chinese）王曉鋒，席光.基于Kriging模型的翼型氣動性能優(yōu)化設(shè)計［J］.航空學(xué)報，2005，26（5）：545－549.

［13］JEONG S，MURAYAMA M，YAMAMOTO K..Efficient optimization design method using Kriging model［J］.Journalof Aircraft，2005，42（2）：413－420.

［14］GOTO Y，JEONG S，OBAYASHI S，et al.Design space exploration of supersonic formation flying focusing on drag minimization［J］.JournalofAircraft，2008，45（2）：430－439.

［15］KANAZAKI M，TANAKA K，JEONG S，et al.Multi－objective aerodynamic exploration of elements＇setting for high－lift airfoil using Kriging model［J］.JournalofAircraft，2007，44（3）：858－864.

［16］HOYLE N，BRESSLOFF N W，KEANE AJ.Design optimization of a two－dimensional subsonic engine air intake［J］.AIAA Journal，2006，44（11）：2672－2681.

［17］SONG W，KEANE A J.Surrogate－based aerodynamic shape optimization of a civil aircraft engine nacelle［J］.AIAAJour－nal，2007，45（10）：2565－2574.

［18］TOAL D J J，BRESSLOFF N W，KEANE A J.Kriging hyperparameter tuning strategy［J］.AIAAJournal，2008，46（5）：1240－1252.

［19］TOAL D J J，BRESSLOFF N W，KEANE A J，et al.The development of a hybridized particle swarm for Kriging hyperparameter tuning［J］.EngineeringOptimization，2011，43（6）：1－28.

［20］KRIGE D G.A statistical approach to some basic mine valuations problems on the witwatersrand［J］.JournaloftheChemical，MetallurgicalandMiningEngineeringSocietyofSouth Africa，1951，52（6）：119－139.

［21］SACKS J，WELCH W J，MITCHELL，et al.Design and analysis of computer experiments［J］.StatisticalScience，1989，4：409－423.

［22］JONES D R.A taxonomy of global optimization methods based on response surfaces［J］.JournalofGlobalOptimization，2001，21：345－383.

［23］LIU J，HAN Z H，SONG W P.Comparison of infill sampling criteria in Kriging－based aerodynamic optimization［A］.28thInternational Congress of the Aeronautical Sciences，Brisbane，Australia［C］.2012，9.

［24］KULFAN B M.Universal parametric geometry representation method［J］.JournalofAircraft，2008，45（1）：142－158.