基于Kriging 模型和網(wǎng)格變形的翼型優(yōu)化

胡勇軍 高一星

（貴州貴飛飛機(jī)設(shè)計(jì)研究院有限公司，貴州 安順 561000）

機(jī)翼通常是飛行器產(chǎn)生升力和阻力的主要部件，而在機(jī)翼設(shè)計(jì)中作為基本要素之一的翼型對飛行器氣動(dòng)特性和飛行性能都有很大的影響，因此選用擁有突出的氣動(dòng)特性的翼型對飛行器氣動(dòng)設(shè)計(jì)非常重要。早期，傳統(tǒng)的翼型設(shè)計(jì)方式都是首先取得已有翼型的試驗(yàn)資料，然后再結(jié)合翼型理論和一定的設(shè)計(jì)經(jīng)驗(yàn)等，加以剪裁和修形，最后才通過風(fēng)洞試驗(yàn)進(jìn)行驗(yàn)證。該設(shè)計(jì)過程往往要求工程設(shè)計(jì)人員必須具備較高水平的專業(yè)技術(shù)知識(shí)和較豐富的工程設(shè)計(jì)經(jīng)驗(yàn)，隨著現(xiàn)代飛機(jī)設(shè)計(jì)中設(shè)計(jì)周期的急劇降低，這種方法使用的越來越少。二十世紀(jì)七十年代起，由于計(jì)算機(jī)技術(shù)的快速發(fā)展，計(jì)算流體力學(xué)（CFD）同時(shí)獲得了飛速發(fā)展，并成為空氣動(dòng)力學(xué)中發(fā)展最快的方向之一，也成為了飛行器設(shè)計(jì)過程中必不可少的技術(shù)手段。CFD 為飛機(jī)翼型設(shè)計(jì)帶來了全新的技術(shù)手段，將原由風(fēng)洞完成的大量試驗(yàn)變?yōu)榱嘶贑FD 的大量仿真計(jì)算再加少量風(fēng)洞試驗(yàn)驗(yàn)證，從而大大縮短了周期，同時(shí)提高了獲得具有優(yōu)秀氣動(dòng)特性的翼型的可能性。

在基于CFD 的翼型設(shè)計(jì)方法中，數(shù)值優(yōu)化設(shè)計(jì)方法是常用的方法之一[1-2]。該方法可以選取目標(biāo)翼型的指定迎角的升力系數(shù)、升阻比、阻力、俯仰力矩等氣動(dòng)特性作為設(shè)計(jì)目標(biāo)，疊加諸如相對厚度分布等設(shè)計(jì)約束，通過遺傳算法等方法，快速獲得滿足要求的翼型。國內(nèi)外均有許多學(xué)者進(jìn)行了該方向的研究。氣動(dòng)外形優(yōu)化設(shè)計(jì)的基本過程就是不斷調(diào)用參數(shù)計(jì)算模型進(jìn)行計(jì)算的過程，如果該過程使用高精度模型用于計(jì)算，計(jì)算過程將會(huì)變得非常復(fù)雜，為降低優(yōu)化過程計(jì)算的復(fù)雜程度，使用的模型最好使用近似方法建立，即用一個(gè)簡單近似函數(shù)代替高精度的復(fù)雜的計(jì)算模型，比如采用局部近似、全局近似等。在全局近似法中，Kriging 模型因其所具有的非線性擬合能力強(qiáng)、近似模型參數(shù)選取靈活性高等優(yōu)勢，在翼型優(yōu)化等氣動(dòng)優(yōu)化設(shè)計(jì)實(shí)踐中得到了廣泛應(yīng)用，如Meunier 就通過該方法建立了機(jī)翼氣動(dòng)力模型[3]等。

本文在基準(zhǔn)翼型上，采用網(wǎng)格變形和CFD 方法獲得試驗(yàn)樣本點(diǎn)數(shù)據(jù)，通過Kriging 模型優(yōu)化設(shè)計(jì)[4-5]出具有要求的目標(biāo)特性的翼型，采用該翼型設(shè)計(jì)的機(jī)翼隨全機(jī)一起進(jìn)行風(fēng)洞試驗(yàn)，飛機(jī)氣動(dòng)特性滿足全機(jī)設(shè)計(jì)要求，成功驗(yàn)證了翼型設(shè)計(jì)。

1 Kriging 模型

響應(yīng)面法的基本思想是假設(shè)一個(gè)包含一些未知參數(shù)的相對簡單的顯式函數(shù)表達(dá)式，通過試驗(yàn)設(shè)計(jì)對給定的設(shè)計(jì)集合進(jìn)行試驗(yàn)，然后用插值等方式確定表達(dá)式中的未知參數(shù)，得到顯式的響應(yīng)面模型以逼近目標(biāo)函數(shù)，用來估計(jì)非試驗(yàn)點(diǎn)的目標(biāo)響應(yīng)值。Kriging 模型方法屬于一種響應(yīng)面法，在1951 年由南非的Krige 首次提出，該方法由于有對試驗(yàn)點(diǎn)無偏估計(jì)、非線性逼近能力較好等優(yōu)點(diǎn)，在高度非線性的優(yōu)化問題中得以大量應(yīng)用。該方法主要利用對未知點(diǎn)附近的信息的加權(quán)線性組合來估計(jì)該點(diǎn)的值，其本質(zhì)就是通過將估計(jì)值的誤差方差最小化來確定加權(quán)系數(shù)。

Kriging 模型[6]由線性回歸部分（多項(xiàng)式部分）和非參數(shù)部分（多項(xiàng)式的偏離部分）組成。假設(shè)在一個(gè)d 維設(shè)計(jì)空間中，有n 個(gè)樣本點(diǎn)的集合X=[x1，x2，…，xn]T，單個(gè)樣本點(diǎn)為xi=[x，x…，x]，樣本點(diǎn)處的響應(yīng)值為Y=[y1，y2，…，yn]T，則采用kriging 函數(shù)構(gòu)造如下形式近似模型。

其中，f (X，β) 為用線性組合的多項(xiàng)式表示的回歸函數(shù)，是線性回歸部分，β 為回歸模型系數(shù)，整個(gè)部分也可理解為在設(shè)計(jì)空間中建立的全局近似，該多項(xiàng)式的階次通常可取零階、一階、二階等，一般不會(huì)取高于二階，而據(jù)尹大偉等研究[7]，通?；貧w模型的階次對精度的影響并不是很明顯；z(X)是均值為零的非參數(shù)模型，屬于隨機(jī)函數(shù)部分，它也可理解為平均數(shù)為零的隨機(jī)的服從正態(tài)分布但協(xié)方差不為0 的局部逼近。隨機(jī)函數(shù)z(X)代表了真實(shí)響應(yīng)與假設(shè)的全局逼近的局部偏差，其統(tǒng)計(jì)學(xué)特性及樣本點(diǎn)之間協(xié)方差為：

其中，σ2為隨機(jī)過程方差；R 是相關(guān)函數(shù)矩陣；r 則是與二個(gè)樣本點(diǎn)xi，xj之間的距離相關(guān)的空間相關(guān)函數(shù)；θ為空間相關(guān)函數(shù)的比例參數(shù)向量，代表了兩個(gè)樣本點(diǎn)之間的相關(guān)性隨兩點(diǎn)距離增大而衰減的程度?？臻g相關(guān)函數(shù)可以是高斯函數(shù)、指數(shù)函數(shù)、樣條函數(shù)等。在使用高斯函數(shù)時(shí)，輸入向量xi、xj的第k 個(gè)分量的空間相關(guān)函數(shù)為：

其中，θk是組成向量θ 中的第k 個(gè)量，dk為：

相關(guān)函數(shù)的選取既具體決定了代理模型中如何擬合數(shù)據(jù),也反映了與全局逼近的局部偏差,一般選取相對比較簡單的相關(guān)函數(shù)。在本文中，選擇的是ANSYS WORKBENCH 自帶的高斯函數(shù)。

2 網(wǎng)格變形

在翼型優(yōu)化中，對翼型的處理和描述有很多種方法，比如采用CST 參數(shù)化描述翼型[8]，采用網(wǎng)格變形獲得翼型等。采用CST 參數(shù)化方法描述翼型時(shí)，參數(shù)化表示方法直接影響翼型設(shè)計(jì)方法的合理性，進(jìn)而使得翼型表示方法成為采用該方法設(shè)計(jì)翼型的重要領(lǐng)域之一。而網(wǎng)格變形是通過對網(wǎng)格模型的形狀進(jìn)行修改的幾何處理技術(shù)，網(wǎng)格變形技術(shù)避免了反復(fù)進(jìn)行CAD 設(shè)計(jì)更改和再生成網(wǎng)格的環(huán)節(jié)，是近年來幾何處理領(lǐng)域的熱點(diǎn)研究問題之一。用于網(wǎng)格變形的方法通常有自由變形技術(shù)、骨架驅(qū)動(dòng)變形技術(shù)、基于網(wǎng)格曲面變形技術(shù)、基于徑向基函數(shù)的網(wǎng)格變形技術(shù)等。

本文在基準(zhǔn)翼型網(wǎng)格基礎(chǔ)上，使用ANSYS 自帶的網(wǎng)格快速變形技術(shù)，建立局部網(wǎng)格變形區(qū)域的拓?fù)淇刂茀?shù)，通過參數(shù)值變化控制網(wǎng)格的隨移變形，改變網(wǎng)格節(jié)點(diǎn)位置，實(shí)現(xiàn)基于純網(wǎng)格的設(shè)計(jì)變形，獲得不同的氣動(dòng)外形網(wǎng)格。

3 流動(dòng)控制方程及邊界條件

本文CFD 計(jì)算使用的控制方程組為基于Reynolds時(shí)均法的方程組，湍流模型為SST K-ω 模型，邊界條件為壓力遠(yuǎn)場邊界條件。方程組的求解方法則選擇了基于壓力求解的方法，通過耦合隱式方法進(jìn)行有限體積離散，而梯度離散方法采用基于網(wǎng)格點(diǎn)的Green-Gauss 方法。本文采用NASA-LANGLEY LS(1)-0413 作為初始翼型，進(jìn)行多目標(biāo)和約束的翼型優(yōu)化設(shè)計(jì)。

4 翼型氣動(dòng)優(yōu)化設(shè)計(jì)

4.1 設(shè)計(jì)變量選取

初始翼型網(wǎng)格見圖1。圖中同時(shí)給出了翼型上部網(wǎng)格的變形區(qū)域，變形區(qū)域從翼型前緣到后緣，可用于上部網(wǎng)格變形的7 個(gè)控制點(diǎn)均勻分布在弦線上且不含前后緣點(diǎn)，變形方向?yàn)樯舷伦冃?，從前往后分別為控制點(diǎn)1至控制點(diǎn)7。同理，用于翼型下部網(wǎng)格變形的7 個(gè)控制點(diǎn)也均勻分布在弦線上且不含前后緣點(diǎn)，變形方向也為上下方向，從前往后分別為控制點(diǎn)8 至控制點(diǎn)14。控制點(diǎn)合計(jì)14 個(gè)，控制點(diǎn)變形參數(shù)作為輸入變量。為了分別控制翼型上表面和下表面變形，雖然控制點(diǎn)分布相同，但區(qū)分用于上下表面控制的控制點(diǎn)。網(wǎng)格變形過程中前后緣點(diǎn)位置保持不變。為避免極端情況下出現(xiàn)負(fù)體積和上下表面交錯(cuò)的情況，各允許變形控制點(diǎn)的變形參數(shù)范圍需在優(yōu)化前通過網(wǎng)格變形測試獲得，控制點(diǎn)參數(shù)設(shè)置時(shí)作為限制條件?？刂泣c(diǎn)1 至控制點(diǎn)12 的輸入?yún)?shù)均為0.05、控制點(diǎn)13 輸入?yún)?shù)為0.03、控制點(diǎn)14 輸入?yún)?shù)為0 時(shí)，網(wǎng)格變形結(jié)果見圖2。在輸出變量中，包括設(shè)計(jì)目標(biāo)所需的2 個(gè)迎角的升力系數(shù)、阻力系數(shù)。此外，為了控制厚度分布，控制點(diǎn)2 與控制點(diǎn)9 的差量，控制點(diǎn)3 與控制點(diǎn)10 的差量，控制點(diǎn)4 與控制點(diǎn)11 的差量，這3 個(gè)量也作為輸出變量。輸出變量合計(jì)7 個(gè)。

圖1 基準(zhǔn)翼型網(wǎng)格

圖2 某輸入?yún)?shù)下網(wǎng)格變形結(jié)果

4.2 目標(biāo)和設(shè)計(jì)約束

目標(biāo)翼型在Ma=0.4 巡航狀態(tài)需具有較高的升阻比、較小的阻力以具有較長的航時(shí)，在Ma=0.15 大迎角狀態(tài)還需具有較高的升力系數(shù)。由于飛機(jī)目標(biāo)為按照機(jī)翼安裝角-3 度設(shè)計(jì)并在0 度迎角附近巡航，翼型設(shè)計(jì)目標(biāo)為：（1）Ma=0.4，在迎角-3 度升力系數(shù)不低于0.15 且盡量大，阻力系數(shù)不高于0.0053；（2）Ma=0.15，在迎角10 度時(shí)升力系數(shù)不低于1.7 且盡量大。設(shè)計(jì)約束條件主要是厚度條件，相對厚度可以減小但不能低于12%。

4.3 樣本分布生成

在建立Kriging 模型前，需要完成試驗(yàn)設(shè)計(jì)（DOE）以獲得樣本點(diǎn)。試驗(yàn)設(shè)計(jì)中，通過取樣策略確定樣本點(diǎn)的數(shù)量和樣本空間分布。常用的樣本設(shè)計(jì)方法一般有正交設(shè)計(jì)、中心復(fù)合設(shè)計(jì)、均勻設(shè)計(jì)、隨機(jī)取點(diǎn)設(shè)計(jì)、拉丁超立方設(shè)計(jì)等方法,也可以自己創(chuàng)建樣本設(shè)計(jì)方法。抽取的樣本點(diǎn)希望覆蓋比較均勻，又希望顯著減少樣本數(shù)量以降低計(jì)算規(guī)模和計(jì)算周期。本文中，采用的是ANSYS WORKBENCH 中的Custom+Sampling 方法先建立總數(shù)為300 的初始試驗(yàn)樣本，在必要時(shí)將后續(xù)優(yōu)化驗(yàn)證結(jié)果也添加到試驗(yàn)樣本中以有效填充設(shè)計(jì)空間。

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

在ANSYS WORKBENCH 中采用上述方法建立優(yōu)化流程，并采用CFD 方法獲得試驗(yàn)樣本作為訓(xùn)練樣本。在獲取了各樣本點(diǎn)的具體數(shù)據(jù)之后，采用Kriging 模型構(gòu)建通過翼型設(shè)計(jì)參數(shù)計(jì)算氣動(dòng)特性的近似模型，在上述約束條件下采用多目標(biāo)優(yōu)化方法進(jìn)行翼型優(yōu)化，不進(jìn)行目標(biāo)函數(shù)歸一化處理以便直接獲得設(shè)計(jì)點(diǎn)處的設(shè)計(jì)值。優(yōu)化設(shè)計(jì)前后翼型的氣動(dòng)特性結(jié)果對比見表1，優(yōu)化設(shè)計(jì)前后的翼型形狀對比見圖3。優(yōu)化后，翼型相對厚度從13%減小到12.5%。從表1 可見，經(jīng)過優(yōu)化設(shè)計(jì)后，在Ma=0.4 時(shí)，阻力系數(shù)增加較小，但升力系數(shù)顯著提高，升阻比提高較多；低速大迎角時(shí)的升力系數(shù)也有所增加。

表1 優(yōu)化設(shè)計(jì)前后對比

圖3 優(yōu)化前后翼型

5 全機(jī)試驗(yàn)驗(yàn)證

采用上述優(yōu)化后的翼型生成機(jī)翼，裝配于機(jī)身，按照1:7.5 的比例縮比設(shè)計(jì)和制造了風(fēng)洞試驗(yàn)?zāi)Ｐ停贔L-8風(fēng)洞中完成了全機(jī)風(fēng)洞試驗(yàn)。全機(jī)試驗(yàn)?zāi)Ｐ鸵妶D4，全機(jī)典型狀態(tài)試驗(yàn)結(jié)果見圖5、圖6。采用優(yōu)化后的翼型設(shè)計(jì)生成的機(jī)翼及全機(jī)，在典型狀態(tài)的氣動(dòng)力試驗(yàn)結(jié)果滿足了設(shè)計(jì)要求，達(dá)到了預(yù)期目的。

圖4 試驗(yàn)?zāi)Ｐ蛨D

圖5 風(fēng)洞試驗(yàn)結(jié)果（升力系數(shù)～阻力系數(shù)）

圖6 風(fēng)洞試驗(yàn)結(jié)果（升力系數(shù)～迎角）

6 結(jié)論

本文通過使用網(wǎng)格變形技術(shù)和CFD 方法獲得翼型優(yōu)化的試驗(yàn)樣本，采用Kriging 模型進(jìn)行了多目標(biāo)多約束條件下的低速翼型優(yōu)化設(shè)計(jì)研究，獲得了滿足設(shè)計(jì)目標(biāo)的翼型。將優(yōu)化獲得的翼型隨全機(jī)進(jìn)行了風(fēng)洞試驗(yàn)驗(yàn)證，全機(jī)氣動(dòng)特性滿足設(shè)計(jì)目標(biāo)要求。上述方法設(shè)計(jì)的翼型通過全機(jī)試驗(yàn)驗(yàn)證，表明該翼型優(yōu)化設(shè)計(jì)方法可以快速有效地用于飛機(jī)翼型設(shè)計(jì)，極大地推動(dòng)了目標(biāo)飛機(jī)進(jìn)一步的氣動(dòng)設(shè)計(jì)工作。