繆佶, 龔春林, 李春娜
(西北工業大學 航天學院, 陜西 西安 710072)
小型無人飛行器由于其具有尺寸小、使用靈活、任務執行能力強等多項優勢,近些年來在各個領域都受到了廣泛關注。為了使飛行器具有更寬的飛行包線,更廣泛的應用范圍,其翼型需要擁有更優良的氣動性能。因此,對小型無人飛行器翼型進行氣動優化設計具有重要的意義。
氣動優化對于飛行器設計是一種有效的輔助方法。早期的氣動優化設計采用工程估算或面元法等低精度模型,單次分析耗時非常短,可以直接進行優化。隨著CFD理論和高性能計算機的發展,CFD數值模擬由于具有較高精度而越來越多地被應用于氣動優化問題求解。然而,由于CFD計算耗時很長,難以應用于大設計空間的全局優化。
基于代理模型的優化方法,由于能夠大大縮減CFD計算成本,逐漸成為氣動優化研究的重要分支和關鍵技術[1]。代理模型的概念最早在20世紀80年代被提出,其作用之一是在分析和優化設計過程中替代那些比較復雜和費時的分析[2]。代理模型方法不僅可以提高優化設計效率,而且有利于濾除數值噪聲和實現并行優化設計[3]。常用的代理模型方法主要包括:徑向基神經網絡模型[4]、多項式響應面[5]、Kriging模型等。然而直接構建代理模型也存在一定的局限性,當樣本點少時難以保證模型精度,樣本點多時計算耗時很長。對于外形復雜、設計變量較多的氣動優化問題,為了確保代理模型的精度,往往需要大量的樣本點,從而導致計算成本過大。為了解決這一問題,變復雜度建模方法(variable-fidelity modeling,VFM)被提出[6],其主要原理是:同時引入高精度分析模型和低精度分析模型,通過計算高、低精度分析模型在樣本點處的差異構建出變復雜度模型來描述優化問題,從而在保證優化精度的同時有效降低計算成本[7]。
Chang等[8]提出使用乘法標度函數得到低精度與高精度模型的差異,并據此建立變復雜度模型。Alexandrov等[9]提出了一階加法標度和一階乘法標度函數,基于標度函數建立了變復雜度模型并分別用于三維機翼和二維翼型的優化問題。Gano等[10]分別通過進行Taylor一階和二階展開的方式建立標度函數。該方法雖然能夠降低模型的構建成本,但該建模方法屬于局部近似建模,并不適用于全局優化問題。Leifsson等[11]提出了一種基于VFM的快速優化方法,其中大量的低精度氣動數據通過結合半工程方法快速計算得到,少量高精度氣動數據通過CFD計算得到,從而實現對機翼氣動外形的快速優化。Zhang等[12]將VFM和期望改進(expected improvement,EI)加點準則相結合,提出VF-EI方法,并結合分層Kriging代理模型,高效求解全局優化問題。Song等[13]提出了一種新的變復雜度優化策略,初始只構建基于低精度樣本點的Kriging代理模型,高精度樣本點則通過EI等加點策略添加到樣本點集中,不斷對低精度代理模型進行修正,構建變復雜度代理模型,該方法能夠有效降低由初始高精度樣本點帶來的計算成本。Han等[14]通過構建多級變復雜度分層Kriging代理模型的方法對二維翼型和三維機翼進行優化,相比單精度代理模型和兩級變復雜度代理模型方法,其進一步降低了計算成本。
當前的變復雜度模型通過少量計算成本大的高精度樣本數據修正大量計算成本小的低精度樣本數據。為了縮減計算成本,氣動優化問題中的低精度樣本可以由以下方式得到[15]:對復雜外形進行簡化,CFD數值模擬過程中采用網格數量少的粗糙網格。Leifur等[15]通過半工程方法計算低精度樣本的氣動參數,誘導阻力通過渦格法求得,零升阻力和激波阻力由片條理論得到,計算一個低精度樣本的氣動數據大約需要15 s,其結果與高精度的偏差在33%~50%。Alexandrov等[9]通過構建粗糙CFD網格的方式得到低精度樣本響應值,計算成本僅為高精度樣本的四分之一。Jiang等[7]將變復雜度模型應用于小水線面雙體船的設計中,低精度樣本響應值也通過粗糙CFD網格計算得到,其網格量為高精度樣本的五分之一。Mifsud等[16]為得到低精度樣本數據,不僅構建了粗糙CFD網格,且僅采用了一階精度進行迭代計算,從而進一步降低了低精度模型的計算耗時,高精度樣本數據則采用細網格和三階精度迭代計算得到。
針對小型無人飛行器的翼型氣動優化問題,為了使計算成本維持在較低水平,并保證優化結果精度足夠,本文提出一種基于CFD收斂提前終止和變復雜度模型的兩級氣動優化方法:首先,通過CFD收斂提前終止的方法獲得低精度樣本響應值,通過CFD完全收斂獲得高精度樣本響應值,并利用加法橋接函數構建變復雜度模型。為得到精確的最優解,先后利用多島遺傳算法和Hooke-Jeeves算法分別基于變復雜度模型和高精度CFD分析進行全局-局部兩級優化。最后,采用該優化方法對小型無人飛行器的翼型優化問題進行求解,表明了該優化方法的可行性和有效性。
Kriging模型最早由南非地質學家 Krige提出,并成功應用于地質勘查中。至今,Kriging模型已成為優化設計領域最具代表性、應用最廣泛的代理模型構建方法之一[17]。Kriging模型是一種預測方差最小的無偏估計模型,具有較好的擬合能力,適用于強非線性問題的擬合[18]。
(1)
式中:f(x)=[f1(x),f2(x),…,fj(x)],j=1,2,…,k為回歸函數;β={β1,…,βj}T為回歸系數;z(x)為滿足N(0,σ2)正態分布的隨機變量,方程右側兩項分別代表預測值的總體趨勢和局部偏差。定義n×k的矩陣
F=[fT(x(1)),fT(x(2)),…,fT(x(n))]T
(2)
在整個設計空間內,不同樣本點之間的關聯程度通過相關函數來描述。常用相關函數有高斯函數、指數函數、線性函數、球函數、三次函數和樣條函數等[19]。本文選用高斯函數作為相關函數,則樣本點的協方矩陣表示為
Cov[z(xi),z(xj)]=σ2R[z(xj),z(xj)]
(3)
相關函數可以表示為
(4)
式中:d為自變量的維數;θk為相關函數的參數。通過隨機過程理論,得到未知點x處Kriging模型的預測值如下
(5)

(6)
(7)
相關參數θ可以通過最大化似然函數求得,其中似然函數[19]的表達式為
(8)
進而可以得到未知點x處的預測方差表達式如下
(9)
單獨基于CFD數據建模會導致建模代價過大,而單獨基于低精度數據建模會難以保證模型精度。為了兼顧計算代價大的高精度模型和精度低的低精度模型,建立變復雜度模型是一種十分有效的方法。變復雜度方法可通過少量高、低精度樣本的差異構建差異響應模型,從而修正其他低精度響應值以建立代理模型或直接針對修正后的響應值進行優化,變復雜度模型通常可以表示為
(10)

首先,對低精度模型和高精度模型選取m個相同的采樣點并計算相應的高、低精度響應值,然后計算m個采樣點處高低響應值間的差異,并建立差異響應模型。最后將差異響應模型作為低精度與高精度響應的橋接,修正低精度響應值,從而建立變復雜度模型。修正方法可采用乘法標度或加法標度,表示如下
在加法修正和乘法修正中,C(x,a)是差異響應模型,預測高、低精度模型的差異。在實際應用中,加法標度函數在氣動優化問題中適用性更廣,而乘法標度函數可能出現奇異點,故本文采用加法標度函數構建變復雜度模型。由于低精度樣本的計算成本很低,不需要選取大量低精度樣本構建代理模型,可以直接基于修正后的響應值進行優化。
構建變復雜度模型所基于的低精度樣本與高精度樣本數據存在一定誤差。文獻[15]通過半工程方法得到低精度樣本零升阻力、誘導阻力以及激波阻力,與高精度樣本結果的相對誤差范圍為33%~50%,這種過大的差異會導致需要更多的高精度樣本,使得計算成本增加。根據經驗,高、低精度數據的相對差異在20%以內往往可以獲得精度較高的變復雜度模型,且低精度樣本的計算耗時應維持在較低水平,以確保優化成本較低。
文獻[8]通過使用粗糙CFD網格的方法獲得低精度數據,使得低精度模型構建成本僅為高精度模型的1/8。然而,粗糙網格受幾何模型的影響較大,外形越復雜的模型所需要的網格數越多。若對外形復雜的模型生成粗糙網格,由于優化過程中幾何外形在不斷地改變,有可能因網格量少或網格尺寸大,導致網格質量過差,從而使得CFD結果精度過低或網格生成過程失敗。因此,采用粗糙網格的方法存在不穩定因素,會影響外形復雜模型優化的可靠性。為了在保證結果精度足夠的同時盡可能縮減計算成本,本文采用CFD收斂提前終止的方法獲得低精度數據。

圖1 CFD收斂提前終止示意圖
CFD的收斂過程往往初始收斂很快,之后隨著迭代步數的增加其收斂速度逐漸減緩。由圖1所示,氣動力系數僅需很少的迭代步數即可接近完全收斂的解,若提前終止收斂,則計算成本遠遠小于完全收斂,且結果相差不大。而且,對于高、低精度樣本如果采用一樣的細網格,網格質量能得到保證,使得氣動優化過程具有很好的魯棒性。
由于殘差反映了當前流場與收斂的流場之間的差距,可以采用殘差來作為流場提前終止收斂的判定準則。對于網格質量高的細網格模型,殘差往往需要收斂到1×10-6才能使升阻力系數完全平穩,但殘差越小氣動力系數變化越小且收斂越緩慢,因此不等殘差完全收斂即提前終止的話可以大大減少計算耗時。
本文采用RAE2822翼型作為小型無人飛行器的基準翼型。為了獲得收斂提前終止準則,分別對殘差收斂到1×10-2,1×10-3和1×10-4時的計算耗時和計算精度進行對比。采用RAE2822翼型的標準算例來對比不同殘差收斂標準,該翼型的計算狀態為:Ma=0.729,Re=6.5×106,α=2.31°。
首先通過網格收斂性驗證,生成RAE2822翼型的細網格A1。網格收斂性驗證結果如表1所示:

表1 網格收斂性驗證
收斂性驗證后確定的細網格A1的網格單元數量為82 643,壁面附面層滿足Y+≈1,細網格A1的情況如圖2所示。
基于該網格,對不同收斂準則下的計算耗時和計算精度進行對比,情況如表2所示。
綜合考慮計算耗時和計算精度,選擇1×10-3作為提前收斂的殘差,此時計算耗時很少,僅為80 s,且與完全收斂結果的相對誤差不大,均在20%以內,能夠滿足變復雜度模型所需的低精度數據要求。因此,本文將殘差作為CFD收斂提前終止的條件:當殘差達到1×10-3時即終止收斂。

圖2 細網格

表2 不同收斂準則對比情況
為了確保細網格A1的計算結果準確,將其與相同計算狀態的實驗結果進行對比,對比情況如表3所示。為了將收斂提前終止方法與粗糙網格方法進行對比,除細網格A1外,再生成3個粗糙程度不一的網格M1,M2,M3,維持壁面Y+和邊界層網格層數不變,適當降低壁面網格密度,M1,M2和M3的網格單元量分別為54 841,26 409和13 195。分別將這4個不同精度的網格用于CFD分析,并與CFD收斂提前終止方法得到的結果進行對比,結果如表3所示。

表3 采用不同粗糙度網格的CFD結果對比
可以看出,細網格A1得到的結果與實驗結果十分接近,升阻力的相對誤差可由(13)至(14)式計算得到
(13)
(14)
升力系數誤差δCl約為1.2%,阻力系數誤差δCd僅為0.8%,故可將A1網格視為基準網格。M1網格的計算結果與基準網格差別不大,與實驗結果也比較接近。其運行時間相比基準網格減小了接近一半,但計算一個樣本點仍然需要430 s的時間。M2的網格數為M1的一半,由于網格數的減小,網格質量有少量下降,依然需要運行360 s才能使結果收斂,M2得到的氣動力系數仍較為接近真實結果。M3的網格數進一步減小了一半,但此時得到的結果已遠遠偏離了真實結果,且由于網格數量少,網格質量下降較多,仍需要300 s才能基本收斂,在優化過程的其他狀態下存在發散的可能。
然而,細網格A1收斂提前終止得到的結果精確度雖然不如M1和M2網格,但與真實結果誤差在20%以內,其升力系數Cl的誤差為16.8%,阻力系數Cd的誤差為4.9%,可以滿足構建變復雜度模型的低精度樣本需求。此外,其計算時間僅80 s,遠遠低于各粗糙網格的計算成本,更有利于全局優化問題的效率提高。
上述算例表明,隨著網格量的減小,CFD的結果誤差會逐漸加大。當網格數量減小到一定量時,相對誤差會急劇增大,由5%增至45%。因此,當網格粗糙到一定程度時,其結果與精確結果的誤差會驟增,從而無法作為低精度樣本參與優化過程。為了生成合適的粗糙網格,需要做不同粗糙度網格的驗證,否則隨著氣動外形優化的進行,模型的幾何外形一直在變化,生成的粗糙網格非常容易在優化過程中引起迭代發散或精度過低。相比而言,CFD收斂提前終止方法的計算成本低于粗糙網格,相對誤差同樣可以滿足低精度樣本的要求;更重要的是,CFD收斂提前終止的方法不需要生成低精度網格,低精度樣本和高精度樣本都采用同一套細網格,在降低工作量的同時提高了優化過程的魯棒性。
本文提出的基于CFD收斂提前終止和變復雜度模型的兩級氣動優化方法的流程如圖3所示,具體介紹如下:

圖3 優化框架流程圖
第一步 對優化問題建模,包括目標函數,設計變量和約束條件。
第二步 通過拉丁超立方試驗設計方法[21]生成少量用于構建差異響應代理模型的樣本點。
第三步 在樣本點處進行CFD模擬,獲得相應的高、低精度響應值。
第四步 在這些樣本點處計算高、低精度響應值的差異,利用Kriging建模方法,構造差異響應代理模型
C(xi)=F(xi)-fl(xi)
(15)
F(xi)為樣本點xi處高精度響應值,fl(xi)為相同樣本點處的低精度響應值。
第五步 將差異響應代理模型疊加在低精度樣本點上,即可建立變復雜度模型,從而獲得樣本點處變復雜度模型的響應值。隨機選取測試點,對比測試點處的高精度響應值與變復雜度響應值的差異,以驗證變復雜度模型精度是否滿足要求。若其精度足夠,則繼續第六步;若其精度不滿足要求,則在代理模型誤差最大處添加樣本點,并返回第三步。
第六步 基于建立的變復雜度模型進行第一級優化。采用多島遺傳算法直接基于變復雜度模型進行全局優化,獲得最優解。
第七步 以第一級全局優化得到的最優解為初值,進行第二級局部優化。采用Hooke-Jeeves算法直接基于高精度CFD分析進行局部優化,從而獲得更精確的最優解。
類別形狀函數變換方法(class and shape transformation,CST)[22]是一種有效的幾何外形參數表達方法,可以用較少的設計變量描述復雜外形[23],常應用于飛行器翼型和外形的設計。CST方法通過類型函數和形狀函數來構建二維或三維的外形。根據文獻,5階伯恩斯坦多項式即可將原翼型擬合得十分精確,故本文采用5階伯恩斯坦多項式描述RAE2822翼型的幾何外形。因此,氣動優化問題的設計變量一共12個,分別為描述上下翼面幾何外形的12個伯恩斯坦多項式系數增量。
由于高亞聲速狀態下飛行器周圍的流場變化更為復雜,對飛行器的飛行性能要求更高,故應重點考慮高亞聲速下的翼型氣動特性。參考多型高亞聲速無人飛行器的飛行狀態,本優化問題的設計狀態確定如下:Ma=0.75,H=7 500 m,α=1°。優化問題的數學模型可表示為:
(16)
式中:s(x),Cl(x)分別為翼型面積和升力系數;S0(x)為優化前的翼型面積,其值為0.078,Cl0(t)為優化前的升力系數,其值為0.51。對于該翼型優化問題,由于飛行狀態為跨聲速,翼型表面會存在局部激波,故在該飛行工況下優化方向以降低激波阻力為主。而翼型面積大小主要影響摩擦阻力,故該優化問題將翼型面積作為約束條件,以體現該優化方法對減小激波阻力起到的效果。優化后的翼型面積不應小于初始翼型面積,優化后的升力系數不應小于初始升力系數。伯恩斯坦多項式系數增量變化范圍應在-0.02和0.02之間。
小型無人飛行器初始翼型的氣動網格采用前面分析過的細網格A1。CFD數值求解過程中采用有限體積法求解N-S方程,湍流模型為SSTk-ω兩方程模型,空間離散采用二階迎風格式,遠場邊界條件為壓力遠場,物面邊界條件為無滑移壁面。
在構建初始VFM時,采用拉丁超立方取樣方法生成35個樣本點。為了證明低精度樣本數據的有效性,對35個樣本的高、低精度響應分布進行了對比,如圖4所示??梢园l現,低精度響應值與高精度響應值的相對誤差不大,且變化趨勢基本相同,升力系數Cl的平均誤差約為13.69%,阻力系數Cd的平均誤差約為5.04%,說明殘差收斂到1×10-3時,低精度樣本與高精度樣本的相對誤差保持在20%以內。此時的誤差可以滿足低精度數據的要求,即高、低精度樣本可用于構建差異響應代理模型。

圖4 高、低精度響應值對比
為驗證變復雜度模型的精度,隨機選取15個試驗點,分別在試驗點處進行CFD模擬得到高、低精度響應分布,通過差異響應代理模型修正低精度響應值,并將修正后的響應值與高精度響應值進行對比,結果如圖5所示。對比圖4和圖5,變復雜度模型的精度相比低精度模型有了明顯提升,修正后的響應值接近于高精度響應值。修正后的升力系數Cl與高精度結果平均誤差為1.63%,修正后的阻力系數Cd與高精度結果平均誤差為2.49%。因此,所建立的變復雜度模型可以代替高精度模型,并用于氣動優化問題的求解。

圖5 變復雜度模型精度驗證
首先,直接基于變復雜度模型進行第一級優化,優化算法采用多島遺傳算法。相比于傳統遺傳算法,多島遺傳算法具有更好的全局求解能力和更高的計算效率。根據設計變量的個數和優化問題復雜度,多島遺傳算法的控制參數設置如下:島數為5,子種群規模np=20,進化代數g=20,交叉概率pc=0.8,變異概率pm=0.01,島間遷移率pt=0.1。
第一級全局優化收斂歷程如圖6所示,其中圖6a)為目標函數阻力系數Cd的收斂歷程,圖6b)~6c)分別為2個約束條件升力系數Cl和翼型面積S的收斂歷程。圖6表明,目標函數Cd在收斂過程中有了明顯降低,約束條件隨著迭代進行在約束邊界附近波動,故第一級全局優化效果比較明顯。
用高精度CFD模型計算最優外形的氣動性能,并計算相應的目標函數值和約束函數值。表4對全局優化得到的最優解與CFD計算的結果進行了對比。

圖6 第一級全局優化迭代歷程
由表4可以看出,變復雜度模型精度較好,優化結果與真實結果相差很小。其中,優化得到的阻力系數Cd與CFD計算結果相對誤差為1.82%,升力系數Cl與CFD計算結果的相對誤差僅為1.46%。上述表明,第一級基于變復雜度模型的最優解與高精度結果相對誤差在2%以內,變復雜度模型具備較好的精度。
對比RAE2822翼型,阻力系數Cd經第一級全局優化有了顯著改善,降低了約13.3 counts;由約束條件分析,升力系數Cl變化較小,依然滿足約束條件,翼型面積S也與基準翼型相同,均為0.078,滿足約束條件。

表4 基于變復雜度模型優化結果與該點對應的CFD分析結果對比
為了提高最優解的質量,以第一級優化得到的最優解為初值,采用Hooke-Jeeves算法直接基于高精度模型進行局部搜索。Hooke-Jeeves算法的原理是從基點出發,通過軸向移動依次沿著不同坐標軸搜索函數值更小的點,作為新的基點,之后沿著相鄰兩個基點的連線方向繼續進行模式搜索,試圖使函數值更快地減小。該方法程序簡單,無需計算梯度,適應性較強,且局部搜索收斂快、效果好,適用于該優化問題。第二級局部優化迭代歷程如圖7所示,根據該算法的原理,選擇目標函數在連續25個迭代步均無法進一步改善作為優化終止判斷標準。當目標函數無法進一步改善時,即認為優化已達到局部最優,故可取優化過程中的目標函數最優值作為該優化問題的局部最優解。

圖7 第二級局部優化迭代歷程

圖8 優化前后翼型對比

圖9 優化前后壓力分布對比
第二級局部優化共調用50次高精度模型,計算耗時約為10 h。可以看出,相比第一輪優化的結果,目標函數Cd在第一輪優化的基礎上繼續降低了2 counts,這說明第二輪優化可以進一步提升優化結果的精度,體現了兩級氣動優化的有效性。約束條件Cl相比第一輪優化結果也僅有微小變化,仍然滿足約束要求。小型無人飛行器的初始RAE2822翼型與優化后翼型的氣動性能對比情況如表5所示。

表5 初始翼型與最優翼型氣動性能對比
小型無人飛行器的初始翼型與優化后翼型的外形和壓力分布對比如圖8和9所示。由圖8外形對比可以看出,翼型中前段彎度變小使壓力峰值前移,壓力恢復變緩,激波變弱。由圖9壓強分布對比可以看出,翼型下表面的壓強分布在優化前后變化不大,主要差異在于翼型上表面。初始翼型在上表面前部壓強緩慢減小,在中部達到峰值后壓強激增,從而產生強激波,而優化后由于壓強峰值前移,整個上表面壓強始終在緩慢增加,因此避免了壓強的激增,激波大大減弱。
圖10a)至10b)分別為小型無人飛行器的初始翼型和優化后翼型的流場壓強云圖??梢钥闯?初始RAE2822翼型在上表面中部存在非常明顯的強激波,而優化后消除了強激波,顯著地減小了激波阻力。

圖10 優化前后壓強云圖對比
采用本文基于CFD收斂提前終止和變復雜度模型的兩級氣動優化方法對飛行器翼型進行優化設計后,目標函數Cd總共降低了15.1 counts,且Cl和S均滿足約束條件,因此提出的氣動優化方法效果良好,可以獲得理想的最優解。
為了評估本文提出的兩級氣動優化方法的性能,采用基于單一精度Kriging代理模型的EGO方法對該翼型優化問題進行了求解,并將優化結果與3.2節中結果進行了對比。其中,CFD計算從網格到流場求解器均與3.1節中設置相同。文獻[24]指出,為了使初始代理模型精度足夠,通過試驗設計獲得的初始樣本點數應該大于m(m+1)/2,m為設計變量個數,故將構建VFM用到的34個樣本點和額外生成的44個樣本點共同作為初始樣本點,用CFD完全收斂的高精度模型求解相應的響應值來構建初始Kriging代理模型。此后,采用EGO方法對該翼型優化問題進行求解。最終,采用EGO方法獲得的最優Cd為0.008 56,相比基準狀態降低了15.7 counts,而本文提出方法獲得的最優Cd為0.008 62,2種方法的優化結果差異不到0.7%,故2種方法獲得的最優解很接近。EGO方法的整個優化過程耗時約為84 h,而本文提出方法的總耗時約為50 h,故本文方法的優化耗時更低。優化結果對比如表6所示。

表6 2種方法結果對比
由此可知,本文提出方法獲得的優化結果精度與EGO方法相近,但本文提出方法的耗時更少,相比EGO方法具有更高的優化效率。由圖11優化后翼型的外形對比也可以看出,2種方法獲得的最優翼型比較相似。

圖11 初始翼型與最優解對比
1) 為提高小型無人飛行器的氣動特性,針對氣動優化設計具有強非線性、多極值的特點,為進一步提高優化效率,本文提出了一種基于CFD收斂提前終止和變復雜度模型的兩級氣動優化方法。
2) 變復雜度模型所需要的低精度樣本響應值通過CFD收斂提前終止的方法獲得,該方法求解精度大于工程和半工程方法;計算耗時和穩定性均優于粗糙網格方法,更適合用于獲得低精度樣本響應值。
3) 采用本文氣動優化方法對小型無人飛行器的翼型進行氣動優化設計。相比于基準翼型,最優翼型的阻力系數下降了15.1 counts,且約束條件均滿足要求。與基于單一精度Kriging代理模型的EGO方法對比后,表明本文提出的方法在優化效率上更具優勢。