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

（西北工業(yè)大學(xué) 翼型葉柵空氣動(dòng)力學(xué)國防科技重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710072）

0 引 言

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

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

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

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

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

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

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

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

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

1 Kriging模型

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

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

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

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

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

假定設(shè)計(jì)空間內(nèi)任意一點(diǎn)的響應(yīng)值可由已知樣本點(diǎn)的響應(yīng)值ys的線性組合來近似，則在未知點(diǎn)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是單位列向量，

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

其中，是方差σ2的估計(jì)值：

1.2 相關(guān)模型

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

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

其中

1.3 Kriging參數(shù)的獲得

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

2 加點(diǎn)準(zhǔn)則

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

2.1 最大化Expected Improvement（EI）

設(shè)樣本集中的最優(yōu)目標(biāo)函數(shù)值為ymin，假設(shè)隨機(jī)變量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ù)相對的改進(jìn)量為I＝y(tǒng)min－Y（x）＞0，于是I的期望值由下式給出：

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

2.2 最小化代理模型 （MP）

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

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

對于完全二階多項(xiàng)式模型的一般公式為：

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

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

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

4 流場求解

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

5 幾何參數(shù)化

5.1 CST參數(shù)化方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2）Kriging—SIC1

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

3）Kriging—SIC2

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

4）Kriging—SIC3

采用Kriging模型，同時(shí)使用EI、MP兩種加點(diǎn)方法，一次增加兩個(gè)樣本點(diǎn)。

流程圖如圖2所示。

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

8 方法驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖13 使用3種加點(diǎn)準(zhǔn)則的Kriging方法及Adjoint方法反設(shè)計(jì)目標(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è)計(jì)

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

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

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

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

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

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

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

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

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

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

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

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

9.2 翼型多點(diǎn)反設(shè)計(jì)

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

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

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

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

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

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

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

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

10 結(jié)論與分析

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

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

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

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

致謝：感謝王樂提供Adjoint翼型反設(shè)計(jì)的結(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）楊旭東，喬志德，朱兵.亞、跨音速三維機(jī)翼氣動(dòng)外形反設(shè)計(jì)的控制理論方法［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，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)面法的跨聲速機(jī)翼氣動(dòng)優(yōu)化設(shè)計(jì)［J］.航空學(xué)報(bào)，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）鄧?yán)?，喬志德，熊俊濤，?多目標(biāo)自然層流翼型反設(shè)計(jì)方法［J］.航空學(xué)報(bào)，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）許瑞飛，宋文萍，韓忠華.改進(jìn)Kriging模型在翼型氣動(dòng)優(yōu)化設(shè)計(jì)中的應(yīng)用研究［J］.西北工業(yè)大學(xué)學(xué)報(bào)，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）蘇偉，白俊強(qiáng).一種代理模型方法及其在氣動(dòng)優(yōu)化設(shè)計(jì)中的應(yīng)用［J］.彈箭與制導(dǎo)學(xué)報(bào)，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模型的翼型氣動(dòng)性能優(yōu)化設(shè)計(jì)［J］.航空學(xué)報(bào)，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.