韓少強(qiáng),宋文萍,韓忠華*,王樂

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

基于梯度增強(qiáng)型Kriging模型的氣動(dòng)反設(shè)計(jì)方法

韓少強(qiáng),宋文萍,韓忠華*,王樂

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

基于Kriging模型的代理優(yōu)化算法目前在氣動(dòng)優(yōu)化設(shè)計(jì)中得到了廣泛應(yīng)用。但在高維(設(shè)計(jì)變量大于30個(gè))氣動(dòng)優(yōu)化中,計(jì)算量過大的問題對其進(jìn)一步發(fā)展產(chǎn)生了嚴(yán)重制約。將翼型和機(jī)翼氣動(dòng)反設(shè)計(jì)問題轉(zhuǎn)化為優(yōu)化問題,采用Adjoint方法進(jìn)行快速梯度求解,利用基于梯度增強(qiáng)型Kriging(GEK)模型的代理優(yōu)化算法分別開展了18、36和108個(gè)設(shè)計(jì)變量的氣動(dòng)反設(shè)計(jì)。首先,通過采用在設(shè)計(jì)空間局部建立GEK模型的方法成功地將基于代理優(yōu)化算法的氣動(dòng)反設(shè)計(jì)問題的維度拓展到了100維以上。其次,研究了梯度計(jì)算精度對基于GEK模型的反設(shè)計(jì)的影響,發(fā)現(xiàn)梯度精度越高,反設(shè)計(jì)的最終效果越好,同時(shí)效率相當(dāng)。最后,通過不同維度的氣動(dòng)反設(shè)計(jì)算例,比較了改進(jìn)擬牛頓法(BFGS)、基于GEK模型和Kriging模型的代理氣動(dòng)反設(shè)計(jì)方法,結(jié)果表明基于GEK模型的代理優(yōu)化算法的效率大幅度高于基于Kriging模型的代理優(yōu)化算法,并且維度越高,效率優(yōu)勢越明顯;同時(shí),基于GEK模型的代理優(yōu)化算法在優(yōu)化效果及分析程序調(diào)用次數(shù)上相比于BFGS方法也略有優(yōu)勢。

設(shè)計(jì)優(yōu)化;Kriging;梯度增強(qiáng)型Kriging(GEK);代理模型;BFGS;翼型反設(shè)計(jì);機(jī)翼反設(shè)計(jì)

Kriging模型具有較好預(yù)測非線性、多峰值函數(shù)的能力,且能給出未知點(diǎn)預(yù)測值的誤差估計(jì)。因此,Kriging代理優(yōu)化算法[1-2]受到了廣泛的重視。近年來,Kriging模型被越來越多地應(yīng)用于氣動(dòng)及多學(xué)科優(yōu)化設(shè)計(jì)中[3-8],在得到優(yōu)良的設(shè)計(jì)結(jié)果的同時(shí)大大減少了優(yōu)化設(shè)計(jì)所需要的計(jì)算時(shí)間。另一方面,相對于優(yōu)化方法,反設(shè)計(jì)方法具有計(jì)算效率高、針對性強(qiáng)和可信度高等特點(diǎn)[9],在飛行器氣動(dòng)設(shè)計(jì)領(lǐng)域仍然是不可或缺的[10]。目前Kriging代理優(yōu)化算法在氣動(dòng)反設(shè)計(jì)中的應(yīng)用較少,并且以設(shè)計(jì)變量個(gè)數(shù)較少的翼型反設(shè)計(jì)為主。Toal等[11-12]通過翼型反設(shè)計(jì),著重比較了使用不同優(yōu)化算法優(yōu)化Kriging超參數(shù)對最終反設(shè)計(jì)結(jié)果的影響。但是他們并非著眼于Kriging代理優(yōu)化算法在翼型反設(shè)計(jì)中的應(yīng)用研究,也未能得到真正有效的反設(shè)計(jì)結(jié)果。劉俊等[13]將翼型氣動(dòng)反設(shè)計(jì)問題轉(zhuǎn)化為優(yōu)化設(shè)計(jì)問題,利用Kriging代理優(yōu)化算法進(jìn)行了翼型的單目標(biāo)、多目標(biāo)氣動(dòng)反設(shè)計(jì);研究了模型加點(diǎn)準(zhǔn)則對設(shè)計(jì)結(jié)果的影響;并和多項(xiàng)式響應(yīng)面模型優(yōu)化方法進(jìn)行了比較研究。結(jié)果表明,Kriging代理優(yōu)化算法的優(yōu)化效果和效率都明顯優(yōu)于基于多項(xiàng)式響應(yīng)面的優(yōu)化方法。

隨著人們對設(shè)計(jì)結(jié)果的性能要求以及可信度要求的提高,一方面更高可信度的數(shù)值模擬程序被應(yīng)用于飛行器的氣動(dòng)設(shè)計(jì)中,這就意味著單次數(shù)值模擬所需的計(jì)算花費(fèi)將急劇增長。另一方面隨著設(shè)計(jì)變量個(gè)數(shù)(設(shè)計(jì)空間維度)的增加,Krig-ing模型精確建模所需要的CFD調(diào)用次數(shù)呈指數(shù)增長[14],這種現(xiàn)象又被稱為“維度之咒”。對于機(jī)翼甚至更復(fù)雜的三維外形優(yōu)化,設(shè)計(jì)變量多,Kriging優(yōu)化過程需要大規(guī)模調(diào)用耗時(shí)的CFD分析。因此,針對高維問題亟待發(fā)展更加高效的氣動(dòng)優(yōu)化設(shè)計(jì)方法。

為了提高Kriging代理優(yōu)化效率,梯度增強(qiáng)型Kriging(Gradient-Enhanced Kriging,GEK)模型在Kriging模型中引入梯度信息,即在建立代理模型時(shí)除了利用樣本點(diǎn)處的函數(shù)值,同時(shí)還要利用樣本點(diǎn)處的梯度值,從而提高模型的精度。這樣,建立和Kriging模型相同精度的GEK模型,需要的樣本點(diǎn)數(shù)量將大大減少。在氣動(dòng)設(shè)計(jì)領(lǐng)域,當(dāng)前已具備高效計(jì)算梯度的技術(shù)——Adjoint方法[15-17],只需額外花費(fèi)一次與CFD相當(dāng)?shù)挠?jì)算量即可獲得目標(biāo)函數(shù)對所有設(shè)計(jì)變量的梯度。也就是說,當(dāng)使用高效梯度計(jì)算方法時(shí),建立相同樣本點(diǎn)數(shù)的GEK模型只需比Kriging模型增加較少的計(jì)算量,而建模精度卻可以得到大幅度提高,從而充分發(fā)揮GEK模型的優(yōu)勢。Liu[18]和Chung[19]等分別提出了一種間接GEK模型和一種直接GEK模型。間接GEK模型利用一階泰勒展開將梯度值變成附加樣本點(diǎn)的函數(shù)值,然后按常規(guī)Kriging模型方法建立近似模型。GEK模型直接將樣本點(diǎn)的梯度值和樣本點(diǎn)的函數(shù)值一起加權(quán)到模型的預(yù)測值中。Laurenceau等[20-21]比較了Kriging、直接GEK、間接GEK模型在氣動(dòng)問題中的預(yù)測精度,發(fā)現(xiàn)在相同樣本點(diǎn)數(shù)下兩種GEK模型的精度都高于Kriging模型。Han等[22]發(fā)展了一種更一般化的直接GEK模型,允許梯度信息在任意位置引入,且允許在某些樣本點(diǎn)處不使用梯度信息。Han等在此基礎(chǔ)上,又發(fā)展了一種新的GEK模型[23],通過建立一系列子GEK模型來代替原有GEK模型,從而大幅度降低相關(guān)矩陣的維度,提高GEK建模的效率。宋超等[24]利用梯度信息,選取擾動(dòng)步長得到初始樣本點(diǎn)附近的派生點(diǎn),以派生點(diǎn)擬合出能夠預(yù)測目標(biāo)函數(shù)趨勢的低精度Kriging模型。然后以初始樣本點(diǎn)修正該模型得到高精度的Kriging模型。顯然,該方法并未直接將梯度信息耦合進(jìn)Kriging建模過程中。

目前,GEK模型已在氣動(dòng)優(yōu)化設(shè)計(jì)中得到一些初步應(yīng)用。Chung[25]和Xuan[26]等將GEK模型與序列二次規(guī)劃(SQP)法相結(jié)合,進(jìn)行了初步氣動(dòng)優(yōu)化設(shè)計(jì),他們用有限差分法來計(jì)算目標(biāo)函數(shù)的梯度,但未對優(yōu)化效率及優(yōu)化效果進(jìn)行研究。劉俊等[27]進(jìn)行了基于Kriging模型和GEK模型的翼型氣動(dòng)反設(shè)計(jì)比較研究,使用Adjoint方法計(jì)算梯度,發(fā)現(xiàn)在計(jì)算量一致的前提下,GEK模型優(yōu)化結(jié)果更好且收斂更快。對于設(shè)計(jì)變量較多(大于30個(gè))的優(yōu)化問題,Kriging模型的“維度之咒”逐漸凸顯。GEK模型能否在一定程度上緩解“維度之咒”,在設(shè)計(jì)變量個(gè)數(shù)上有所突破,這成了當(dāng)前代理優(yōu)化算法在氣動(dòng)優(yōu)化應(yīng)用中鮮有涉及卻十分重要的研究方向。另外,由于GEK模型每次建模過程中用到的響應(yīng)值和梯度值遠(yuǎn)遠(yuǎn)多于梯度算法每次尋優(yōu)過程所用到的響應(yīng)值和梯度值,所以有必要將基于GEK模型的代理優(yōu)化算法和改進(jìn)擬牛頓(BFGS)方法進(jìn)行比較研究。

本文使用GEK代理優(yōu)化算法進(jìn)行了翼型及機(jī)翼的壓力分布反設(shè)計(jì),包括18維翼型氣動(dòng)反設(shè)計(jì)和36維、108維機(jī)翼氣動(dòng)反設(shè)計(jì),并和Kriging代理優(yōu)化方法、BFGS方法進(jìn)行了對比研究。通過在設(shè)計(jì)空間局部建立GEK模型,使得基于代理優(yōu)化算法的氣動(dòng)反設(shè)計(jì)的維度突破了100維,“維度之咒”得到明顯緩解。

1 優(yōu)化設(shè)計(jì)方法

1.1 幾何參數(shù)化與設(shè)計(jì)變量范圍確定

在利用優(yōu)化方法進(jìn)行翼型優(yōu)化設(shè)計(jì)之前先要用一些設(shè)計(jì)參數(shù)來代替翼型,即對翼型進(jìn)行參數(shù)化。CST參數(shù)化方法是由波音公司的工程師Kulfan提出來的。該方法用一個(gè)類函數(shù)C(ˉx)和一個(gè)型函數(shù)Sˉ()x 的乘積加上一個(gè)描述后緣厚度的函數(shù)Δyte來表示一個(gè)翼型的幾何形狀:

關(guān)于CST參數(shù)化求解方法,具體可參見文獻(xiàn)[28]。

同時(shí),還要給定優(yōu)化過程中設(shè)計(jì)參數(shù)的變化范圍。利用一個(gè)較“厚”的翼型和一個(gè)較“薄”的翼型,并使得這兩個(gè)翼型之間的范圍足夠大,可以包含目標(biāo)翼型。將兩個(gè)翼型之間的空間作為真實(shí)設(shè)計(jì)空間,如圖1中陰影部分所示,c為翼型弦長。兩個(gè)翼型所對應(yīng)的CST參數(shù)即為設(shè)計(jì)變量上下限。機(jī)翼反設(shè)計(jì)采用控制剖面翼型的CST參數(shù)作為設(shè)計(jì)變量。

1.2 優(yōu)化問題數(shù)學(xué)模型

氣動(dòng)反設(shè)計(jì)問題轉(zhuǎn)化為優(yōu)化設(shè)計(jì)問題后,目標(biāo)函數(shù)可以寫為

式中:pd為目標(biāo)壓力分布;p為反設(shè)計(jì)得到的壓力分布;B為沿物面積分;S為積分面。

1.3 試驗(yàn)設(shè)計(jì)方法

在構(gòu)建Kriging模型之前,需要采用試驗(yàn)設(shè)計(jì)方法對設(shè)計(jì)空間進(jìn)行抽樣,來生成一定數(shù)量的初始樣本點(diǎn)。

1.3.1 拉丁超立方試驗(yàn)設(shè)計(jì)

翼型反設(shè)計(jì)中一般初始樣本點(diǎn)個(gè)數(shù)大于設(shè)計(jì)變量個(gè)數(shù),采用拉丁超立方抽樣[29](Latin Hypercube Sampling,LHS)試驗(yàn)設(shè)計(jì)方法生成初始樣本點(diǎn)。它的基本原理為:如果進(jìn)行ns次抽樣(即總共ns個(gè)樣本點(diǎn)),那么就把nv個(gè)設(shè)計(jì)變量都分成等概率的ns個(gè)區(qū)間,于是整個(gè)抽樣空間就被分為等概率的(nv)ns個(gè)小區(qū)間。然后將ns個(gè)樣本{0,1,…,ns-1}隨機(jī)地放入設(shè)計(jì)空間中并滿足下列兩個(gè)條件:

1)每一個(gè)樣本點(diǎn)在小區(qū)間內(nèi)隨機(jī)地分布。

2)將所有樣本點(diǎn)投影到任意一維時(shí),每一個(gè)小區(qū)間內(nèi)有且僅有一個(gè)樣本點(diǎn)。

1.3.2 蒙特卡羅試驗(yàn)設(shè)計(jì)

當(dāng)初始樣本點(diǎn)數(shù)小于設(shè)計(jì)變量個(gè)數(shù)時(shí),很多試驗(yàn)設(shè)計(jì)方法已經(jīng)無法實(shí)現(xiàn)而蒙特卡羅(Monte Carlo,MC)抽樣[30]仍然適用。機(jī)翼反設(shè)計(jì)中由于設(shè)計(jì)變量個(gè)數(shù)較多,初始樣本點(diǎn)個(gè)數(shù)一般少于設(shè)計(jì)變量個(gè)數(shù),所以采用蒙特卡羅抽樣試驗(yàn)設(shè)計(jì)方法生成初始樣本點(diǎn)。其實(shí)現(xiàn)比較簡單,假設(shè)給定區(qū)間[xl, xu],xu和xl為上下限,在該區(qū)間內(nèi)生成一個(gè)隨機(jī)數(shù),該隨機(jī)數(shù)即所生成的樣本點(diǎn)。對于nv維空間,即設(shè)計(jì)變量為nv,在每一維空間內(nèi)生成一個(gè)隨機(jī)數(shù),這nv個(gè)隨機(jī)數(shù)即構(gòu)成了設(shè)計(jì)空間內(nèi)的一個(gè)樣本點(diǎn)。

1.4 梯度增強(qiáng)型Kriging模型

本文不對Kriging模型進(jìn)行詳細(xì)介紹,詳見文獻(xiàn)[31-33]。GEK模型是在Kriging模型基礎(chǔ)上發(fā)展起來的,其思想、原理和建模方法與Kriging模型基本一致。但是GEK模型不僅利用了樣本點(diǎn)處的函數(shù)值,還利用了樣本點(diǎn)處的梯度值。下面介紹本文所采用的由Han等[22]發(fā)展的一種直接GEK模型。

首先,在Kriging模型中,對于m維問題,假設(shè)對目標(biāo)函數(shù)(或狀態(tài)變量)抽樣獲得n個(gè)函數(shù)值,抽樣位置及相應(yīng)的樣本點(diǎn)響應(yīng)值可寫成

式中:S為樣本矩陣,每一行代表一個(gè)樣本點(diǎn);ys為列向量,每一個(gè)元素代表一個(gè)樣本點(diǎn)的函數(shù)值。建立近似模型的目的就是利用已知的樣本數(shù)據(jù)(S,ys),來獲得任意未知點(diǎn)處的預(yù)測值。

對于GEK模型,將上述樣本數(shù)據(jù)的響應(yīng)值列向量加以拓展,把樣本點(diǎn)處的梯度信息以偏導(dǎo)數(shù)的形式加入到響應(yīng)值中。于是得到含梯度信息的樣本及響應(yīng)值矩陣為式中:xg(j)(j∈{1,2,…,n′})為第j個(gè)提供偏導(dǎo)數(shù)的樣本點(diǎn)(每個(gè)偏導(dǎo)數(shù)看作一個(gè)樣本點(diǎn));向量ys?y(j)中的?y(j)= (?x)表示第j個(gè)偏導(dǎo)數(shù)信息(對

k第k個(gè)獨(dú)立變量的偏導(dǎo)數(shù),k∈ {1,2,…,m})。

需要說明的是,這里梯度信息的引入具有一般性:梯度信息的抽樣位置可以與函數(shù)值的抽樣位置相同,也可不同;可以引入所有偏導(dǎo)數(shù)信息,也可只引入部分偏導(dǎo)數(shù)信息。也就是說,總的偏導(dǎo)數(shù)個(gè)數(shù)n′可以小于n·m。如果沒有任何偏導(dǎo)數(shù)信息,GEK模型將退化為傳統(tǒng)的Kriging模型。因此,可以認(rèn)為,這種GEK模型是一種可考慮梯度信息的更一般化的Kriging模型。

隨著優(yōu)化設(shè)計(jì)問題維度的增加,如果在優(yōu)化過程中采用所有的樣本點(diǎn)進(jìn)行GEK建模,模型訓(xùn)練過程將面臨非常巨大的計(jì)算量(對于108維氣動(dòng)反設(shè)計(jì)算例,每個(gè)樣本點(diǎn)都有1個(gè)響應(yīng)值和108個(gè)梯度值,以100個(gè)樣本點(diǎn)為例,GEK模型需構(gòu)建的相關(guān)矩陣大小將是10 900×10 900)。在本文的108維算例中,將反設(shè)計(jì)問題作為局部優(yōu)化問題處理,GEK模型每次建模只使用當(dāng)前所有樣本點(diǎn)中距離當(dāng)前最優(yōu)點(diǎn)最近的10個(gè)樣本點(diǎn),這樣每次建立局部代理模型時(shí)相關(guān)矩陣的大小都可以控制為1 090×1 090,大幅減小了優(yōu)化建模部分的計(jì)算量。事實(shí)證明這種方法在高維機(jī)翼氣動(dòng)反設(shè)計(jì)中是高效可行的。

1.4.1 GEK預(yù)測值與均方誤差

GEK模型對未知函數(shù)的預(yù)測值定義為所有抽樣函數(shù)值和所有抽樣偏導(dǎo)數(shù)值的加權(quán),即

式中:w(i)為第i個(gè)函數(shù)值的加權(quán)系數(shù);λ(j)為第j個(gè)偏導(dǎo)數(shù)信息的加權(quán)系數(shù)。與Kriging模型類似,引入靜態(tài)隨機(jī)過程假設(shè):

式中:β0∈1為未知常數(shù);Z(·)為均值為0的靜態(tài)隨機(jī)過程,且協(xié)方差滿足

式中:σ2為隨機(jī)過程Z(·)的方差;R為僅與空間距 離 有 關(guān) 的 “相 關(guān) 函 數(shù) ”;和分別為相關(guān)函數(shù)相對于x(i)和x(j)處第k個(gè)分量的偏導(dǎo)數(shù)。

與Kriging模型推導(dǎo)過程類似,可得到GEK模型的預(yù)測值為

式中:R為樣本點(diǎn)相關(guān)矩陣;r為樣本點(diǎn)與未知點(diǎn)x的相關(guān)向量;F=[1 … 1 0 … 0]T;β0=(FTR-1F)-1FTR-1ys。

均方誤差MSE為

1.4.2 相關(guān)函數(shù)模型

相關(guān)矩陣R和相關(guān)向量r的建立需要計(jì)算相關(guān)函數(shù)及其一階導(dǎo)數(shù)、二階導(dǎo)數(shù)值。假設(shè)兩個(gè)樣本點(diǎn)x(i)和x(j)之間的相關(guān)函數(shù)僅與兩點(diǎn)之間的距離有關(guān)(即R(x(i),x(j))=R(x(i)-x(j))),于是,此處考慮以下形式的相關(guān)模型:

本文采用3次樣條相關(guān)函數(shù):

1.4.3 參數(shù)模型

由式(10)可知,GEK模型中存在超參數(shù)θ,

與Kriging模型相同,該超參數(shù)通過采用最大似然估計(jì)(MLE)法求解下列優(yōu)化問題獲得

1.5 加點(diǎn)準(zhǔn)則及其約束處理

代理優(yōu)化方法的一個(gè)關(guān)鍵技術(shù)是加點(diǎn)準(zhǔn)則。由于建立一個(gè)較精確的全局近似模型需要大量的樣本點(diǎn),特別是設(shè)計(jì)變量較多時(shí),因此,為了減少總的樣本點(diǎn)數(shù),首先使用較少的樣本點(diǎn)建立初始代理模型,再尋找一定準(zhǔn)則下的最優(yōu)點(diǎn)并將其作為新的樣本點(diǎn)來更新代理模型重新優(yōu)化。加點(diǎn)準(zhǔn)則的使用直接決定了優(yōu)化的成敗及優(yōu)化效率。

本文的氣動(dòng)反設(shè)計(jì)優(yōu)化算例中,加點(diǎn)方法采用MSP(Minimum of Surrogate Prediction)[34]準(zhǔn)則。該準(zhǔn)則認(rèn)為代理模型足夠精確,直接尋找代理模型上目標(biāo)函數(shù)的最小值,即建立代理模型后,求解下列子優(yōu)化問題:

式中:NG為約束個(gè)數(shù);^y(x)和^gi(x) 分別為目標(biāo)函數(shù)和約束函數(shù)的代理模型。本文使用遺傳算法和梯度優(yōu)化相結(jié)合的方法求解式(13)子優(yōu)化問題,可得到代理模型上預(yù)測的目標(biāo)函數(shù)最優(yōu)解x*。再對x*進(jìn)行精確數(shù)值模擬分析,并將結(jié)果作為新的樣本數(shù)據(jù),添加到現(xiàn)有樣本數(shù)據(jù)集中,重新建立代理模型,直至整個(gè)優(yōu)化過程收斂。該準(zhǔn)則是一種局部加點(diǎn)準(zhǔn)則,在氣動(dòng)反設(shè)計(jì)問題中有較好的效果。

1.6 優(yōu)化設(shè)計(jì)流程

本文氣動(dòng)反設(shè)計(jì)優(yōu)化采用自主開發(fā)的優(yōu)化軟件SurroOpt,是結(jié)合代理模型、試驗(yàn)設(shè)計(jì)方法、多點(diǎn)加點(diǎn)準(zhǔn)則、遺傳算法、梯度優(yōu)化算法等發(fā)展的代理優(yōu)化設(shè)計(jì)方法。圖2給出了代理模型優(yōu)化設(shè)計(jì)方法的流程圖,圖中DoE(Design of Experiment)表示試驗(yàn)設(shè)計(jì)。

2 流場求解程序驗(yàn)證

采用自主開發(fā)的二維和三維流場求解程序PMNS2D/3D,通過求解雷諾平均Navier-Stokes方程對翼型及機(jī)翼進(jìn)行氣動(dòng)分析?？臻g離散采用中心有限體積法,LU-SGS隱式時(shí)間推進(jìn),使用了隱式殘值光順、當(dāng)?shù)貢r(shí)間步長、多重網(wǎng)格等加速收斂措施。

2.1 RAE2822翼型繞流計(jì)算

翼型為RAE2822亞臨界翼型,計(jì)算網(wǎng)格如圖3所示。全湍流計(jì)算,湍流模型為剪切應(yīng)力輸運(yùn)(SST)模型,計(jì)算狀態(tài)為馬赫數(shù)Ma=0.72,迎角α=2.8°,雷諾數(shù)Re=6.5×106。根據(jù)圖4可以看出:計(jì)算壓力系數(shù)Cp分布與試驗(yàn)值吻合較好,驗(yàn)證了二維流場求解程序的可靠性。

2.2 M6機(jī)翼繞流計(jì)算

M6機(jī)翼的表面網(wǎng)格和空間網(wǎng)格如圖5所示,網(wǎng)格量為50萬。全湍流計(jì)算,采用Spalart-Allmaras湍流模型,計(jì)算狀態(tài)為Ma=0.839 5,α=3.06°,Re=1.17×107。由圖6可以看出,M6機(jī)翼不同站位的壓力分布計(jì)算值與試驗(yàn)值吻合較好,驗(yàn)證了三維流場求解程序的可靠性。圖中:z/b表示機(jī)翼展向位置坐標(biāo)與半展長的比值。

3 Adjoint梯度求解程序驗(yàn)證

采用自主開發(fā)的二維和三維梯度求解程序,通過求解連續(xù)Adjoint方程得到反設(shè)計(jì)目標(biāo)函數(shù)對所有設(shè)計(jì)變量的梯度。并和二階精度有限差分方法計(jì)算的梯度進(jìn)行比較,來驗(yàn)證Adjoint梯度求解程序的可靠性。

3.1 RAE2822翼型壓力分布反設(shè)計(jì)梯度

目標(biāo)壓力分布為RAE2822翼型在計(jì)算狀態(tài)下的壓力分布。計(jì)算翼型為NACA0012翼型,計(jì)算網(wǎng)格的類型和規(guī)模與圖3所示的網(wǎng)格一致。計(jì)算狀態(tài)為 Ma=0.72,α=2.8°,Re=6.5×106。根據(jù)圖7可以看出2種方法計(jì)算的梯度值吻合較好,驗(yàn)證了二維Adjoint梯度求解程序的可靠性。

3.2 M6機(jī)翼壓力分布反設(shè)計(jì)梯度

目標(biāo)壓力分布為M6機(jī)翼在計(jì)算狀態(tài)下的表面壓力分布。計(jì)算機(jī)翼的翼根和翼梢剖面形狀均為NACA0012翼型,計(jì)算網(wǎng)格的類型和規(guī)模與圖5所示網(wǎng)格一致。計(jì)算狀態(tài)為Ma=0.839 5,α=3.06°,Re=1.17×107。根據(jù)圖8可以看出2種方法計(jì)算的梯度值吻合較好,驗(yàn)證了三維Adjoint梯度求解程序的可靠性。

4 反設(shè)計(jì)算例與分析

使用基于GEK模型的代理優(yōu)化方法進(jìn)行了18、36、108維的氣動(dòng)反設(shè)計(jì)研究,同時(shí)和Kriging代理優(yōu)化設(shè)計(jì)方法及BFGS方法進(jìn)行了對比。需要說明的是,GEK模型建模和BFGS方法所需要的梯度信息均采用Adjoint方法求解,用Adjoint方法求解一次反設(shè)計(jì)目標(biāo)函數(shù)對所有設(shè)計(jì)變量的梯度,花費(fèi)的計(jì)算量相當(dāng)于一次CFD求解。對于三維外形,CFD分析的計(jì)算量遠(yuǎn)大于代理模型建模的計(jì)算量,所以僅以CFD分析次數(shù)作為計(jì)算量統(tǒng)計(jì)的依據(jù)。

4.1 RAE2822翼型18維反設(shè)計(jì)

以RAE2822翼型單目標(biāo)反設(shè)計(jì)為例。設(shè)計(jì)狀態(tài)為Ma=0.734,α=2.79°,Re=6.5×106。在求解Navier-Stokes方程時(shí),采用SST湍流模型進(jìn)行全湍流計(jì)算。

使用8階CST參數(shù)化方法,共18個(gè)設(shè)計(jì)變量。使用基于GEK模型的代理優(yōu)化方法開展了RAE2822翼型氣動(dòng)反設(shè)計(jì),并和Kriging代理優(yōu)化方法及BFGS方法進(jìn)行了對比研究。取相同的設(shè)計(jì)空間,代理優(yōu)化方法取20個(gè)初始樣本點(diǎn),加點(diǎn)準(zhǔn)則為MSP。

根據(jù)表1展示的優(yōu)化設(shè)計(jì)結(jié)果可以看出,GEK代理優(yōu)化方法的效率是Kriging代理優(yōu)化方法的3倍,最終得到的目標(biāo)函數(shù)值是Kriging代理優(yōu)化方法的1/5。GEK代理優(yōu)化方法的效率是BFGS的2.5倍,最終得到的目標(biāo)函數(shù)值略小于BFGS方法。

表1 RAE2822翼型反設(shè)計(jì)結(jié)果Table 1 Inverse design results of RAE2822 airfoil

根據(jù)圖9所示的目標(biāo)函數(shù)收斂歷程,可以看出相比Kriging模型,GEK模型使得翼型代理優(yōu)化的效率大幅提高。說明梯度信息可以提高代理模型精度,從而使用更少的樣本點(diǎn)即可鎖定最優(yōu)點(diǎn)的位置。本算例中,相對于BFGS方法,GEK代理優(yōu)化方法也具有較明顯的效率優(yōu)勢,最終優(yōu)化效果也更好。

根據(jù)圖10(a)和圖10(b)可以看出,使用3種優(yōu)化方法的反設(shè)計(jì)翼型外形及壓力分布與目標(biāo)翼型及壓力分布吻合較好,說明3種方法都達(dá)到了較好的反設(shè)計(jì)效果。

通過控制連續(xù)Adjoint方程求解殘差的下降量級(jí),可以間接控制Adjoint梯度的求解精度。分別將Adjoint求解殘差下降1、2、3個(gè)量級(jí)的梯度用于基于GEK模型的代理優(yōu)化。從表1的優(yōu)化結(jié)果和圖11給出的優(yōu)化收斂曲線可以看出,GEK模型所使用的梯度精度越高,最終的反設(shè)計(jì)目標(biāo)函數(shù)越小,反設(shè)計(jì)效果越好,同時(shí)計(jì)算量相當(dāng)。

4.2 M6機(jī)翼36維反設(shè)計(jì)

以M6機(jī)翼單目標(biāo)反設(shè)計(jì)為例。設(shè)計(jì)狀態(tài)為Ma=0.839 5,α=3.06°,Re=1.17×107。

剖面翼型使用8階CST參數(shù)化方法,翼根翼梢兩個(gè)剖面共36個(gè)設(shè)計(jì)變量,通過兩個(gè)剖面沿展向線性插值即可得到M6機(jī)翼。使用基于GEK模型的代理優(yōu)化方法開展了M6機(jī)翼氣動(dòng)反設(shè)計(jì)優(yōu)化,并和Kriging代理優(yōu)化方法及BFGS方法進(jìn)行了對比研究。取相同的設(shè)計(jì)空間,代理優(yōu)化方法取10個(gè)初始樣本點(diǎn),加點(diǎn)準(zhǔn)則為MSP。

根據(jù)表2展示的優(yōu)化設(shè)計(jì)結(jié)果可以看出,基于Kriging模型的代理優(yōu)化方法花費(fèi)了將近8倍于GEK模型的計(jì)算量,最終得到的目標(biāo)函數(shù)是GEK模型的2.5倍。GEK代理優(yōu)化方法的優(yōu)化效果略優(yōu)于BFGS方法,但BFGS調(diào)用的流場求解次數(shù)是GEK代理優(yōu)化方法的2.6倍。

根據(jù)圖12所示的目標(biāo)函數(shù)收斂歷程,可以看出在36維的氣動(dòng)反設(shè)計(jì)問題中,Kriging模型的“維度之咒”已經(jīng)開始顯現(xiàn)。此時(shí)采用GEK模型可以使機(jī)翼反設(shè)計(jì)的效率大幅度提高。

表2 M6機(jī)翼反設(shè)計(jì)結(jié)果Table 2 Inverse design results of M6 wing

圖13(a)和圖13(b)給出了反設(shè)計(jì)機(jī)翼和M6機(jī)翼的19%站位剖面翼型的外形及壓力分布比較。圖14(a)和圖14(b)給出了反設(shè)計(jì)機(jī)翼和M6機(jī)翼的79%站位剖面翼型的外形及壓力分布比較?？梢钥闯鍪褂?種方法的反設(shè)計(jì)機(jī)翼外形及壓力分布與目標(biāo)機(jī)翼及壓力分布均吻合較好。

根據(jù)圖15所示的目標(biāo)函數(shù)收斂歷程可以看出,GEK模型所使用的梯度精度越高,最終目標(biāo)函數(shù)越小,反設(shè)計(jì)效果越好。

4.3 CRM機(jī)翼108維反設(shè)計(jì)

對第五屆AIAA阻力預(yù)測研討會(huì)(DPW-5)所使用的由NASA開發(fā)的CRM機(jī)翼進(jìn)行氣動(dòng)反設(shè)計(jì)(Ma=0.85,α=2.16°,Re=5.0×106)。

CRM機(jī)翼平面形狀比較復(fù)雜,以圖16所示的9個(gè)剖面翼型形狀參數(shù)作為設(shè)計(jì)變量進(jìn)行壓力分布反設(shè)計(jì)。將9個(gè)剖面沿展向樣條插值后即可近似得到CRM機(jī)翼。采用5階CST參數(shù)化方法,對剖面翼型進(jìn)行參數(shù)化,每個(gè)剖面有12個(gè)設(shè)計(jì)變量,9個(gè)剖面共108個(gè)設(shè)計(jì)變量。

CRM機(jī)翼繞流計(jì)算使用的網(wǎng)格為C-H結(jié)構(gòu)化網(wǎng)格,如圖17所示。網(wǎng)格量約140萬,遠(yuǎn)場約為翼根弦長的30倍。

如前文所述,為了減少高維問題中GEK模型訓(xùn)練所花費(fèi)的時(shí)間,在設(shè)計(jì)空間局部建立GEK模型來進(jìn)行CRM機(jī)翼氣動(dòng)反設(shè)計(jì),并和Kriging代理優(yōu)化方法及BFGS方法進(jìn)行對比研究。取相同的設(shè)計(jì)空間,代理優(yōu)化取100個(gè)初始樣本點(diǎn),使用MSP加點(diǎn)準(zhǔn)則。

根據(jù)表3給出的CRM機(jī)翼代理氣動(dòng)反設(shè)計(jì)結(jié)果可以看出,使用GEK模型最終得到的目標(biāo)函數(shù)大小約為Kriging模型的1/65,而花費(fèi)的計(jì)算量不到Kriging模型的1/2。預(yù)計(jì)在本算例中,要達(dá)到相當(dāng)?shù)膬?yōu)化效果,基于Kriging模型的代理優(yōu)化方法花費(fèi)的計(jì)算量將是基于GEK模型的代理優(yōu)化方法的數(shù)十倍。同時(shí),在108維反設(shè)計(jì)算例中,GEK代理優(yōu)化方法的優(yōu)化效果略優(yōu)于BFGS方法,流場求解的調(diào)用次數(shù)也略少于BFGS方法。需要說明的是,在本算例當(dāng)中,由于只在設(shè)計(jì)空間局部建立GEK模型,GEK代理優(yōu)化算法的優(yōu)勢可能沒有得到充分發(fā)揮。

圖18給出了不同方法反設(shè)計(jì)收斂曲線的比較,圖19～圖21分別給出了反設(shè)計(jì)機(jī)翼0.5%站位(翼根)、40.8%站位及99.4%站位(翼梢)處的剖面翼型及壓力分布和CRM機(jī)翼的比較?？梢钥闯?基于GEK模型的代理優(yōu)化方法使用不到500次CFD計(jì)算得到的外形和CRM機(jī)翼吻合得很好,而基于Kriging模型的代理優(yōu)化方法花費(fèi)了1 000次CFD計(jì)算,得到的外形和CRM機(jī)翼仍有較明顯的差距。

5 結(jié) 論

1)對于翼型和機(jī)翼的氣動(dòng)反設(shè)計(jì),基于GEK模型的代理優(yōu)化方法比基于Kriging模型的代理優(yōu)化方法具有明顯的效率優(yōu)勢。相比于BFGS方法,基于GEK模型的代理優(yōu)化方法仍然有更高的效率,同時(shí)也可以得到更好的優(yōu)化效果。

2)Adjoint梯度求解精度越高,基于GEK模型的代理氣動(dòng)反設(shè)計(jì)效果越好。

3)對大于30個(gè)設(shè)計(jì)變量的高維優(yōu)化問題,基于Kriging模型的代理優(yōu)化開始顯現(xiàn)“維度之咒”。GEK模型由于在建模過程中使用了梯度信息,可以使用較少的樣本點(diǎn)來建立高精度的模型,在高維問題中表現(xiàn)出了巨大的效率優(yōu)勢?；贕EK模型的代理優(yōu)化方法可在一定程度上緩解“維度之咒”,使得代理優(yōu)化問題的維度成功拓展到100維以上。

本文所采用的GEK模型,在設(shè)計(jì)變量和樣本點(diǎn)數(shù)量較多的情況下,建模所花費(fèi)的時(shí)間仍然較多,所以在108維算例中采用了建立局部代理模型的方法來減少計(jì)算成本。今后需要對該模型做進(jìn)一步改進(jìn),提高模型訓(xùn)練效率,在高維問題中利用所有樣本點(diǎn)的響應(yīng)值和梯度值來建立高精度代理模型,從而更加充分地發(fā)揮GEK模型的優(yōu)勢。

致 謝

感謝劉俊博士在優(yōu)化設(shè)計(jì)平臺(tái)開發(fā)方面所做的工作。

[1] 韓忠華.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).

[2] HAN Z H.SurroOpt:A generic surrogate-based optimization code for aerodynamic and multidisciplinary design[C]//30th Congress of the International Council of the Aeronautical Aciences,2016.

[3] 王曉鋒,席光.基于Kriging模型的翼型氣動(dòng)性能優(yōu)化設(shè)計(jì)[J].航空學(xué)報(bào),2005,26(5):545-549.WANG X F,XI G.Aerodynamic optimization design for airfoil based on Kriging model[J].Acta Aeronautica et Astronautica Sinica,2005,26(5):545-549(in Chinese).

[4] 許瑞飛,宋文萍,韓忠華.改進(jìn)Kriging模型在翼型氣動(dòng)優(yōu)化設(shè)計(jì)中的應(yīng)用研究[J].西北工業(yè)大學(xué)學(xué)報(bào),2010,28(4):503-510.XU R F,SONG W P,HAN Z H.Application of improved Kriging-model-based optimization method in airfoil aerodynamic design[J].Journal of Northwestern Polytechnical University,2010,28(4):503-510(in Chinese).

[5] 白俊強(qiáng),王波,孫智偉,等.基于松散式代理模型管理框架的亞音速機(jī)翼優(yōu)化設(shè)計(jì)方法研究[J].西北工業(yè)大學(xué)學(xué)報(bào),2011,29(4):515-519.BAI J Q,WANG B,SUN Z W,et al.Developing optimization design of subsonic wing with losse type of agent model[J].Journal of Northwestern Polytechnical University,2011,29(4):515-519(in Chinese).

[6] 張科施,韓忠華,李為吉,等.基于近似技術(shù)的高亞聲速運(yùn)輸機(jī)機(jī)翼氣動(dòng)/結(jié)構(gòu)優(yōu)化設(shè)計(jì)[J].航空學(xué)報(bào),2006,27(5):810-815.ZHANG K S,HAN Z H,LI W J,et al.Multidisciplinary aerodynamic/structural design optimization for high subsonic transport wing using approximation technique[J].Acta Aeronautica et Astronautica Sinica,2006,27(5):810-815(in Chinese).

[7] 姚拴寶,郭迪龍,孫振旭,等.基于Kriging代理模型的高速列車頭型多目標(biāo)優(yōu)化設(shè)計(jì)[J].中國科學(xué):技術(shù)科學(xué),2013,43(2):186-200.YAO S B,GUO D L,SUN Z X,et al.Multi-objective optimization of the streamlined head of high-speed trains based on the Kriging model[J].Scientia Sinica:Techchnologica,2013,43(2):186-200(in Chinese).

[8] 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&Multidisciplinary Optimization,2016:1-19.

[9] KIM H J,RHO O H.Aerodynamic design of transonic wings using the target pressure optimization approach[J].Journal of Aircraft,1998,35(5):671-677.

[10] 趙國慶,招啟軍.基于目標(biāo)壓力分布的旋翼先進(jìn)氣動(dòng)外形反設(shè)計(jì)分析方法[J].航空學(xué)報(bào),2014,35(3):744-755.ZHAO G Q,ZHAO Q J.Inverse design analysis method on rotor with advanced aerodynamic configuration based upon target pressure distribution[J].Acta Aeronautica et Astronautica Sinica,2014,35(3):744-755(in Chinese).

[11] TOAL D J J,BRESSLOFF N W,KEANE A J.Kriging hyperparameter tuning strategies[J].AIAA Journal,2008,46(5):1240-1252.

[12] TOAL D J J,BRESSLOFF N W,KEANE A J,et al.The development of a hybridized particle swarm for kriging hyperparameter tuning[J].Engineering Optimization,2011,43(6):675-699.

[13] 劉俊,宋文萍,韓忠華,等.Kriging模型在翼型反設(shè)計(jì)中的應(yīng)用研究[J].空氣動(dòng)力學(xué)學(xué)報(bào),2014,32(4):518-526.LIU J,SONG W P,HAN Z H,et al.Kriging-based airfoil inverse design[J].Acta Aerodynamica Sinica,2014,32(4):518-526(in Chinese).

[14] SHAN S.Survey of modeling and optimization strategies for high-dimensional design problems[J].Chemical Engineering Science,2015,141(2):184-194.

[15] JAMESON A,REUTHER J.Control theory based airfoil design using the Euler equations[C]//5th Symposium on Multidisciplinary Analysis and Optimization.Reston:AIAA,1994.

[16] JAMESON A.Aerodynamic design via control theory[J].Journal of Scientific Computing,1988,3(3):233-260.

[17] JAMESON A.Optimum aerodynamic design using CFD and control theory:AIAA-1995-1729[R].Reston:AIAA,1995.

[18] LIU W Y,BATILL S M.Gradient-enhanced response surface approximations using kriging models[C]//9th AIAA/ISSMO Symposium and Exhibit on Multidisciplinary Analysis and Optimization.Reston:AIAA,2002.

[19] CHUNG H S,ALONSO J J.Using gradients to construct cokriging approximation models for high-dimensional design optimization problems:AIAA-2002-0317[R].Reston:AIAA,2002.

[20] LAURENCEAU J,SAGAUT P.Building efficient response surfaces of aerodynamic functions with kriging and cokriging[J].AIAA Journal,2008,46(2):498-507.

[21] LAURENCEAU J,MEAUX M,MONTAGNAC M,et al.Comparison of gradient-based and gradient-enhanced response-surface-based optimizers[J].AIAA Journal,2010,48(5):981-994.

[22] HAN Z H,G?RTZ S,ZIMMERMANN R.Improving variable-fidelity surrogate modeling via gradient-enhanced kriging and a generalized hybrid bridge function[J].Aerospace Science and Technology,2013,25(1):177-189.

[23] HAN Z H,SONG C X,ZHANG Y.A new formulation of gradient-enhanced surrogate model and application to aerodynamic design[C]//34th AIAA Applied Aerodynamics Conference.Reston:AIAA,2016.

[24] 宋超,楊旭東,宋文萍.耦合梯度與分級(jí)Kriging模型的高效氣動(dòng)優(yōu)化方法[J].航空學(xué)報(bào),2016,37(7):2144-2155.SONG C,YANG X D,SONG W P.Efficient aerodynamic optimization method using hierarchical Kriging combined with gradient[J].Acta Aeronautica et Astronautica Sinica,2016,37(7):2144-2155(in Chinese).

[25] CHUNG H S,ALONSO J J.Design of a low-boom supersonic business jet using cokriging approximation models:AIAA-2002-5598[R].Reston:AIAA,2002.

[26] XUAN Y,XIANG J H,ZHANG W H,et al.Gradientbased Kriging approximate model and its application research to optimization design[J].Science in China Series E:Technological Sciences,2009,52(4):1117-1124.

[27] 劉俊,宋文萍,韓忠華,等.梯度增強(qiáng)的Kriging模型與Kriging模型在優(yōu)化設(shè)計(jì)中的比較研究[J].西北工業(yè)大學(xué)學(xué)報(bào),2015,33(5):819-826.LIU J,SONG W P,HAN Z H,et al.Comparative study of GEK(Gradient-Enhanced Kriging)and Kriging when applied to design optimization[J].Journal of Northwestern Polytechnical University,2015,33(5):819-826(in Chinese).

[28] KULFAN B M.Universal parametric geometry representation method[J].Journal of Aircraft,2008,45(1):142-158.

[29] GIUNTA A A,WOJTKIEWICZ S F,ELDRED M S.Overview of modern design of experiments methods for computational simulations[C]//Proceedings of the 41st AIAA Aerospace Sciences Meeting and Exhibit.Reston:AIAA,2003.

[30] METROPOLIS N,ULAM S.The Monte Carlo method[J].Journal of the American Statistical Association,1949,44(247):335-341.

[31] HOYLE N,BRESSLOFF N W,KEANE A J.Design optimization of a two-dimensional subsonic engine air intake[J].AIAA Journal,2006,44(11):2672-2681.

[32] HAN Z H,ZHANG K S,SONG W P,et al.Optimization of active flow control over an airfoil using a surrogatemanagement framework[J].Journal of Aircraft,2010,47(2):603-612.

[33] HAN Z H,ZHANG K S,LIU J,et al.Surrogate-based aerodynamic shape optimization with application to wind turbine airfoils:AIAA-2013-1108[R].Reston:AIAA,2013.

[34] LIU J,HAN Z H,SONG W P.Comparison of infill sampling criteria in Kriging-based aerodynamic optimization[C]//28th International Congress of the Aeronautical Sciences,2012.

Aerodynamic inverse design method based on gradient-enhanced Kriging model

HAN Shaoqiang,SONG Wenping,HAN Zhonghua*,WANG Le

National Key Laboratory of Science and Technology on Aerodynamic Design and Research,School of Aeronautics,Northwestern Polytechnical University,Xi'an 710072,China

Up to date,the surrogate-based optimization methods using Kriging model are widely used in the aerodynamic optimization design.However,the huge time costs seriously restrict the in-depth development of Kriging model when it comes to cases with dimension course(more than 30 design variables).lnverse design problems of airfoils and wings are formulated as optimization problems.Surrogate-based optimization using gradient-enhanced Kriging(GEK)model is adopted to conduct aerodynamic inverse design of 18,36 and 108 variables,and gradient is solved via the efficient adjoint method.lt is worth mentioning that the dimension of aerodynamic inverse design using surrogate-based optimization is successfully expanded to more than 100 by constructing GEK model in part of the design space.Besides,the influence of gradient accuracy on surrogate-based optimization using GEK model is studied.The results show that the effect of optimization is improved with more accurate gradients.BFGS(Broyden,Fletcher,Goldfarb and Shanno)and surrogate-based optimizations using GEK and Kriging models are compared through aerodynamic inverse design of different dimensions.The results show that surrogate-based optimization using GEK model has significantly higher efficiency,and the advantage of efficiency is more obvious with higher dimension of design.At the same time,surrogate-based optimization using GEK model has better optimization effect and less number of calls of analysis programs than the BFGS method.

design optimization;Kriging;gradient-enhanced Kriging(GEK);surrogate model;BFGS;inverse design of airfoil;inverse design of wing

2016-09-26;Revised:2016-10-17;Accepted:2016-11-15;Published online:2016-11-24 10:34

URL:www.cnki.net/kcms/detail/11.1929.V.20161124.1034.006.html

National Natural Science Foundation of China(11272265)

V211.3

A

1000-6893(2017)07-120817-15

10.7527/S1000-6893.2016.0299

2016-09-26;退修日期:2016-10-17;錄用日期:2016-11-15;網(wǎng)絡(luò)出版時(shí)間:2016-11-24 10:34

www.cnki.net/kcms/detail/11.1929.V.20161124.1034.006.html

國家自然科學(xué)基金(11272265)

*通訊作者.E-mail:hanzh@nwpu.edu.cn

韓少強(qiáng),宋文萍,韓忠華,等.基于梯度增強(qiáng)型Kriging模型的氣動(dòng)反設(shè)計(jì)方法[J].航空學(xué)報(bào),2017,38(7):120817.HAN S Q,SONG W P,HAN Z H,et al.Aerodynamic inverse design method based on gradient-enhanced Kriging model[J].Acta Aeronautica et Astronautica Sinica,2017,38(7):120817.

(責(zé)任編輯:李明敏)

*Corresponding author.E-mail:hanzh@nwpu.edu.cn