基于進(jìn)化算法和復(fù)合參數(shù)化方法的低速翼型優(yōu)化

王順順, 郭 正

(國(guó)防科技大學(xué)空天科學(xué)學(xué)院， 湖南長(zhǎng)沙 410000)

引 言

隨著CFD的發(fā)展、 優(yōu)化算法的日趨完善與計(jì)算能力的提升， 特定工況下氣動(dòng)外形優(yōu)化成本越來越低， 使該技術(shù)的工程應(yīng)用越發(fā)成為可能。氣動(dòng)優(yōu)化可以顯著提升飛行器的飛行性能和飛行品質(zhì)， 對(duì)飛行安全、 飛行效率與經(jīng)濟(jì)性都有決定性影響[1]。翼型優(yōu)化則是其中的重要內(nèi)容， 波音公司的Garner等[2]認(rèn)為， 對(duì)一架大型運(yùn)輸機(jī)， 在相同迎角下升力系數(shù)增加0.1， 即可因增大擦尾角縮短起落架而使飛機(jī)質(zhì)量減少635 kg,而同一速度下， 最大升力系數(shù)增大1.5%則意味著載重增加3 000 kg。

翼型優(yōu)化綜合了實(shí)驗(yàn)方法、 參數(shù)化方法、 氣動(dòng)計(jì)算方法、 優(yōu)化算法等多個(gè)學(xué)科。20世紀(jì)70年代Hicks等[3]首次將數(shù)值優(yōu)化思想引入到翼型設(shè)計(jì)中， Yamamoto等[4]第一次將遺傳算法引入到氣動(dòng)外形設(shè)計(jì)中， Quagliarella等[5]等利用遺傳算法對(duì)翼型進(jìn)行了系列優(yōu)化設(shè)計(jì)， Jeong[6]將Krige模型結(jié)合遺傳算法應(yīng)用到翼型優(yōu)化中， 國(guó)內(nèi)如西北工業(yè)大學(xué)等單位也長(zhǎng)期開展翼型優(yōu)化相關(guān)研究， 周旺儀等[7]基于Krige代理模型與多目標(biāo)遺傳算法對(duì)增升裝置進(jìn)行多目標(biāo)優(yōu)化設(shè)計(jì)， 白俊強(qiáng)等[8]利用RBF動(dòng)網(wǎng)格和改進(jìn)粒子群算法進(jìn)行多段翼優(yōu)化,韓忠華等[9]總結(jié)了代理模型在氣動(dòng)優(yōu)化方面的研究進(jìn)展， 李潤(rùn)澤等[10]對(duì)“人在回路”思想在氣動(dòng)優(yōu)化中的發(fā)展進(jìn)行了綜述， 高正紅等[1]、 余雄慶[11]分別對(duì)飛機(jī)總體氣動(dòng)設(shè)計(jì)優(yōu)化進(jìn)行了回顧總結(jié)。總體來講當(dāng)前對(duì)于翼型氣動(dòng)優(yōu)化的研究工作以驗(yàn)證優(yōu)化方法為主， 面向工程應(yīng)用的較少。

無人飛行器是目前航空領(lǐng)域研究的熱點(diǎn)， 如美軍的RQ-7影子無人機(jī)， 是美軍旅團(tuán)級(jí)作戰(zhàn)單元單獨(dú)獲取地面情報(bào)的主力， 為100～500 kg量級(jí)的戰(zhàn)術(shù)無人機(jī)， 該類無人機(jī)攜帶載荷較全、 質(zhì)量較小， 應(yīng)用最為廣泛。但隨著任務(wù)的不斷拓展， 戰(zhàn)術(shù)無人機(jī)對(duì)氣動(dòng)性能也提出了更高要求， 尤其對(duì)滯空時(shí)間與短距起降的需要日趨增加， 而當(dāng)前無人機(jī)翼型研究主要針對(duì)高空長(zhǎng)航時(shí)無人機(jī)， 對(duì)于戰(zhàn)術(shù)無人機(jī)翼型研究較少， 本文研究的目的就是為該量級(jí)無人機(jī)設(shè)計(jì)高質(zhì)量翼型， 優(yōu)化的主要目標(biāo)是綜合提高其起降與巡航性能。

1 復(fù)合參數(shù)化方法

翼型參數(shù)化是用代數(shù)變量表征翼型幾何形狀的過程， 好的參數(shù)化方法可以用盡可能少的參數(shù)表征范圍更廣的翼型， 是翼型優(yōu)化的基礎(chǔ)。目前常用的翼型參數(shù)化方法有CST， PARSEC， Hicks-Henne等。CST方法的主要思想是使用由Bezier多項(xiàng)式構(gòu)造的型函數(shù)修正類別函數(shù)以描述外形變化特征； PARSEC通過基函數(shù)的加權(quán)疊加分別表示翼型上下表面曲線； Hicks-Henne方法是將連續(xù)光滑的Hicks-Henne型函數(shù)加權(quán)疊加到初始翼型上去得到新的翼型輪廓。 廖炎平等[12]通過對(duì)比研究發(fā)現(xiàn)不同的參數(shù)化方法表現(xiàn)翼型的能力有較大的區(qū)別。此外當(dāng)前的參數(shù)化方法已經(jīng)可以用來對(duì)機(jī)翼甚至整機(jī)進(jìn)行優(yōu)化， 逐漸成為了氣動(dòng)優(yōu)化領(lǐng)域的熱點(diǎn)[13]。但單一的參數(shù)化方法往往存在不足， 比如光滑性與細(xì)節(jié)表現(xiàn)力是一對(duì)基本的矛盾， 王迅等[14]通過采用B樣條疊加小波分解進(jìn)行局部光順的方法解決這個(gè)問題， 取得了較好效果。本文通過對(duì)在光滑性與細(xì)節(jié)表現(xiàn)性上各有優(yōu)劣的兩種參數(shù)化方法進(jìn)行數(shù)值疊加， 吸收不同參數(shù)化方法的優(yōu)勢(shì)， 拓展翼型表現(xiàn)范圍。

本文采用NACA 4位數(shù)翼型參數(shù)化方法與Hicks-Henne型函數(shù)法進(jìn)行組合。其中NACA 4位數(shù)翼型是基于NACA變密度風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)提出的系列翼型， 它的特點(diǎn)是特征參數(shù)較少、 可解析描述、 翼型庫豐富、 應(yīng)用廣泛可靠性好[15]， 同時(shí)特征參數(shù)直接反映翼型幾何特征， 便于結(jié)合飛機(jī)設(shè)計(jì)的結(jié)構(gòu)要求對(duì)翼型的幾何特征進(jìn)行約束， 因此極適合作為優(yōu)化的初始翼型庫， 但由于翼型解析公式較為簡(jiǎn)單, 細(xì)節(jié)表現(xiàn)能力較弱。

為了彌補(bǔ)NACA 4位數(shù)翼型細(xì)節(jié)表現(xiàn)能力不足的缺點(diǎn)， 在其坐標(biāo)基礎(chǔ)上采用型函數(shù)線性擾動(dòng)法疊加擾動(dòng)， 形成復(fù)合參數(shù)化方法。翼型的形狀由基準(zhǔn)翼型和擾動(dòng)型函數(shù)的線性疊加共同決定， 一般擾動(dòng)型函數(shù)選用的是Hicks-Henne型函數(shù)[16]。原方法在翼型后緣表現(xiàn)方面有缺陷， 王建軍等[17]對(duì)Hicks-Henne型函數(shù)進(jìn)行了相關(guān)改進(jìn),彌補(bǔ)了這一不足。由于理論上可以采用任意高階型函數(shù)， 具有很強(qiáng)的細(xì)節(jié)表現(xiàn)能力， 但選用參數(shù)范圍過大時(shí)則可能會(huì)出現(xiàn)不適合工程應(yīng)用的奇異翼型。同時(shí)該方法也無法對(duì)翼型幾何特征進(jìn)行有效控制， 所以一般用來進(jìn)行小范圍修型。其表達(dá)式如下， 其中ytop0、ylow0為原翼型上下翼面y坐標(biāo)，xk為控制變化的變量。

最終復(fù)合方法為將上兩種方法進(jìn)行數(shù)值疊加， 如下所示， 其中ylowf與ytopf為最終翼型的縱坐標(biāo)，ytop1與ylow1分別為NACA 4位數(shù)翼型坐標(biāo)， 后面的項(xiàng)為Hicks-Henne型函數(shù)數(shù)值擾動(dòng)。

圖1是3種參數(shù)化方法實(shí)現(xiàn)效果的對(duì)比， 可以看到單純的NACA 4位數(shù)翼型細(xì)節(jié)表現(xiàn)不足， 但可靠性好并容易施加約束， Hicks-Henne型函數(shù)法容易出現(xiàn)奇異翼型,在某些情況下甚至出現(xiàn)不光順的情況， 而通過兩種參數(shù)化方法的疊加組合， 既能夠保證翼型曲線的表現(xiàn)能力和翼型非奇異， 同時(shí)還可以約束彎度、 厚度等幾何特征參數(shù)。

(a) Airfoil represented by NACA four-digit parameterization method

(b) Airfoil represented by Hicks-Henne shape functions

(c) Hicks-Henne functions superimposed on NACA 4412 airfoil圖1 不同翼型參數(shù)化方法之間的比較Fig. 1 Comparison between different airfoil parameterization methods

2 分步計(jì)算方法

為了節(jié)省前期計(jì)算量， 盡快縮小優(yōu)化范圍， 初步優(yōu)化采用Xfoil軟件進(jìn)行氣動(dòng)估算， 使用高階的面元法和完全耦合的黏流/無黏流迭代方法計(jì)算阻力、 邊界層轉(zhuǎn)捩和分離[18]。盡管從理論上講， 利用面元法進(jìn)行流場(chǎng)模擬的近似等級(jí)最低， 但由于Xfoil發(fā)展早， 修正方法比較完善， 工程應(yīng)用性較強(qiáng)， 成為目前使用廣泛的氣動(dòng)估算軟件。經(jīng)檢驗(yàn), 在低Reynolds數(shù)和小迎角工況下Xfoil與CFD計(jì)算結(jié)果具有較好的一致性[19]。精細(xì)優(yōu)化采用求解Reynolds平均 Navier-Stokes (RANS)方程進(jìn)行流動(dòng)數(shù)值計(jì)算， 利用軟件腳本自動(dòng)生成結(jié)構(gòu)網(wǎng)格， 圖2為生成的網(wǎng)格， 可以看到對(duì)于不同翼型自動(dòng)生成的網(wǎng)格均有較好質(zhì)量。

圖2 生成的網(wǎng)格Fig. 2 Generated mesh

求解器采用Fluent軟件， 湍流模型選擇k-ωSST模型， 這是一種用來結(jié)合RANS方程可以有效計(jì)算翼型外流場(chǎng)的湍流模型， 計(jì)算對(duì)象為低速不可壓流。圖3為Re=1×106， 14°迎角條件下NACA 4412翼型的算例驗(yàn)證， 數(shù)據(jù)來自伊利諾伊大學(xué)低速風(fēng)洞實(shí)驗(yàn)結(jié)果， 其中Cp與x分別為無量綱壓力系數(shù)和長(zhǎng)度。

圖3 算例驗(yàn)證Fig. 3 Example verification

實(shí)際優(yōu)化過程中發(fā)現(xiàn)采用分步計(jì)算可以較快地使翼型參數(shù)空間收斂到原范圍的1/5， 小空間有利于后續(xù)優(yōu)化算法進(jìn)行精確尋優(yōu)?？紤]到氣動(dòng)估算基本不考慮時(shí)間成本， 所以采用估算方法與精確求解相結(jié)合的方法較全部采用精確求解方法進(jìn)行優(yōu)化耗費(fèi)計(jì)算時(shí)間大大縮短， 同時(shí)優(yōu)化精度不變。

3 優(yōu)化算法與流程

翼型優(yōu)化往往是典型的多目標(biāo)優(yōu)化問題,本文采用非歸一化的NSGA-Ⅱ 算法。非歸一化算法不用通過設(shè)置權(quán)重將多目標(biāo)轉(zhuǎn)化為單目標(biāo)， 而是使所求解集的前沿與Pareto前沿盡量接近。NSGA-Ⅱ 算法是對(duì)NSGA算法的改良版， 優(yōu)點(diǎn)在于探索性能更好， Pareto前沿前進(jìn)能力增強(qiáng)， 收斂速度加快[20]。本文中群體規(guī)模取為40， 最大進(jìn)化代數(shù)為200， 雜交率為0.9， 變異概率為0.01。圖4為精細(xì)優(yōu)化過程中的2°迎角升力系數(shù)收斂進(jìn)程， 橫坐標(biāo)代表進(jìn)化算法推進(jìn)序列，可以看到在50代以后算法就已經(jīng)收斂到最優(yōu)值附近， 說明該算法的有效性。

圖4 算法收斂進(jìn)程Fig. 4 Algorithm convergence process

本文的優(yōu)化流程為： ①根據(jù)適用工況確定優(yōu)化目標(biāo)、 約束、 初始翼型、 參數(shù)變化范圍； ②利用Xfoil結(jié)合NACA 4位數(shù)翼型參數(shù)化方法進(jìn)行初步優(yōu)化， 選擇一個(gè)合適的優(yōu)化解作為下一步優(yōu)化的起始值； ③確定優(yōu)化參數(shù)范圍， 采用求解Reynolds平均 Navier-Stokes 方程的氣動(dòng)解算方法進(jìn)行精細(xì)優(yōu)化； ④從最優(yōu)解集中比較選擇出合適的最優(yōu)解； ⑤在上一步最優(yōu)解上疊加Hicks-Henne擾動(dòng)型函數(shù)進(jìn)行細(xì)節(jié)優(yōu)化； ⑥通過分析比較得到最優(yōu)解， 檢驗(yàn)優(yōu)化結(jié)果， 流程結(jié)束。

4 優(yōu)化案例

4.1 約束條件與初始條件

本節(jié)采用的初始翼型為NACA 5314， 彎度為5%、 彎度位置30%、 厚度為14%， 如圖5所示， 是NACA 4位數(shù)字翼型中最大升力系數(shù)較高的低速翼型， 同時(shí)具有一定程度的緩失速性能。

圖5 初始翼型Fig. 5 Initial airfoil

本節(jié)設(shè)計(jì)工況為Re=1×106， 以下計(jì)算中Rey-nolds數(shù)保持不變， 為實(shí)現(xiàn)氣動(dòng)性能綜合優(yōu)化的目的， 考慮飛機(jī)起降與巡航階段對(duì)于氣動(dòng)性能的不同要求及機(jī)翼結(jié)構(gòu)的強(qiáng)度要求， 選擇較高迎角下8°的升力系數(shù)CLα=8°與較低迎角下2°的升力系數(shù)CLα=2°和升阻比Kα=2°共3個(gè)參數(shù)作為優(yōu)化目標(biāo)， 以最大厚度t、 最大彎度位置xcamber、 最大彎度fmax為設(shè)計(jì)變量并確保翼型最大厚度t不小于5%， 可以表示為： 尋找一組設(shè)計(jì)變量x， 使

4.2 優(yōu)化過程

初步優(yōu)化中各參數(shù)上下各變化50%作為參數(shù)范圍。最終圖6為翼型初步優(yōu)化前后對(duì)比， 可以發(fā)現(xiàn)翼型厚度明顯降低， 同時(shí)彎度增大， 最大彎度位置后移。

圖6 初步優(yōu)化前后翼型對(duì)比Fig. 6 Comparison of airfoils before and after preliminary optimization

圖7為優(yōu)化前后2°迎角下壓力分布對(duì)比(單位為Pa)， 圖7(b)為優(yōu)化后翼型,可以明顯看到隨著彎度增加后加載明顯增大， 同時(shí)前緣厚度減小導(dǎo)致吸力峰強(qiáng)度略有降低。

(a) Before optimization

(b) After optimization圖7 壓力分布對(duì)比Fig. 7 Pressure distribution comparison

初步優(yōu)化后選取翼型的參數(shù)： 彎度為5.7%、 彎度位置在40.3%、 厚度為10.6%， 以得到的新翼型參數(shù)上下變化20%為新的參數(shù)范圍， 通過求解RANS方程解算氣動(dòng)特性進(jìn)行精確優(yōu)化。精確優(yōu)化得到的部分解集如圖8所示， 得到的最優(yōu)解有3個(gè)， 分別側(cè)重于CLα=8°,CLα=2°,Kα=2°， 如圖8所示， 3個(gè)結(jié)果分別標(biāo)示為Op1， Op2， Op3， 但Op3擁有較均勻的氣動(dòng)特性優(yōu)勢(shì)故被選擇為精確優(yōu)化結(jié)果。

圖8 精確優(yōu)化的解集Fig. 8 Solutions of precise optimization

精確優(yōu)化后翼型與原翼型對(duì)比見圖9， 可以看到翼型厚度減小彎度增大， 最大彎度位置向后移。

圖9 精確優(yōu)化前后翼型對(duì)比Fig. 9 Comparison of airfoils before and after precise optimization

圖10為翼型優(yōu)化前后2°迎角下壓力對(duì)比(單位為Pa)， 圖10(b)為優(yōu)化后翼型， 上部曲率分布更加均勻、 前緣略平緩導(dǎo)致精確優(yōu)化后翼型上部低壓區(qū)范圍明顯增大， 同時(shí)彎度增大引起翼型下方壓力增大。

(a) Before optimization

(b) After optimization圖10 壓力分布對(duì)比Fig.10 Pressure distribution comparison

結(jié)合大量計(jì)算數(shù)據(jù)還可以進(jìn)一步分析翼型幾何參數(shù)與氣動(dòng)性能之間的關(guān)系， 見圖11。

(a) KLα=2°

(b) CLα=2°

(c) CLα=8°圖11 幾何參數(shù)與氣動(dòng)性能的關(guān)系Fig. 11 Relationship between geometric parameters and aerodynamic performance

可以看到翼型的彎度對(duì)升力系數(shù)影響尤為明顯， 這也是大量實(shí)驗(yàn)結(jié)果的結(jié)論， 而厚度的減小則會(huì)減少阻力系數(shù)， 提高升阻比。

對(duì)精確優(yōu)化結(jié)果Op2采用改進(jìn)型Hicks-Henne型函數(shù)攝動(dòng)進(jìn)行細(xì)節(jié)優(yōu)化， 本文共設(shè)置12個(gè)攝動(dòng)變量,上下翼面各6個(gè)， 每個(gè)控制變量的范圍為-0.005～0.005， 各點(diǎn)在翼面上均勻分布。具體細(xì)節(jié)如圖12所示。

圖12 Hicks-Henne型函數(shù)攝動(dòng)Fig.12 Hicks-Henne function perturbation

細(xì)節(jié)優(yōu)化的部分優(yōu)化解集如圖13， 其最優(yōu)解也可以分為3類， 通過比較選擇Op1作為最終優(yōu)化結(jié)果。細(xì)節(jié)優(yōu)化前后翼型對(duì)比見圖14， 可以發(fā)現(xiàn)翼型彎度略有提升至8%， 最大彎度位置略有前移至40%， 翼型厚度略有減小至8%， 最大厚度位置前移至24%， 其2°迎角下壓力分布對(duì)比(單位為Pa)見圖15， 圖15(b)為優(yōu)化后翼型， 吸力峰強(qiáng)度與范圍明顯變大， 且翼型下方高壓區(qū)擴(kuò)大。

圖13 細(xì)節(jié)優(yōu)化的解集Fig. 13 Solutions of detailed optimization

圖14 細(xì)節(jié)優(yōu)化前后翼型對(duì)比Fig. 14 Comparison of airfoils before and after detailed optimization

(a) Before optimization

(b) After optimization圖15 壓力分布對(duì)比Fig.15 Pressure distribution comparison

4.3 優(yōu)化結(jié)果

三輪優(yōu)化結(jié)果對(duì)比如表1所示， 可以看到最終所有目標(biāo)氣動(dòng)性能均有大幅度提升， 以2°迎角下的升阻比為例， 假設(shè)無人機(jī)在該迎角下巡航飛行， 根據(jù)布雷蓋航程公式飛機(jī)航程與升阻比成正比， 則其最大航程也會(huì)比使用NACA 5314顯著提高。優(yōu)化結(jié)果見表1。

表1 優(yōu)化結(jié)果對(duì)比Table 1 Comparison of optimized results

圖16反映了優(yōu)化前后兩種翼型用Xfoil軟件計(jì)算的氣動(dòng)性能之間的對(duì)比， 圖16(a)可以看到經(jīng)過優(yōu)化的翼型在10°迎角以前的升力系數(shù)得到提升， 圖16(b)中反映升力系數(shù)與升阻比的關(guān)系， 可以看到優(yōu)化后的翼型在升力系數(shù)小于1.5時(shí)升阻大于原翼型， 這可以有效提升無人機(jī)的巡航性能。

(a) CL& α

(b) K & CL圖16 優(yōu)化前后氣動(dòng)性能對(duì)比Fig. 16 Aerodynamic performance comparison

5 結(jié)論

本文基于復(fù)合參數(shù)化方法結(jié)合進(jìn)化算法， 對(duì)翼型在2°和8°兩種迎角下的氣動(dòng)性能進(jìn)行3輪優(yōu)化， 優(yōu)化后的翼型氣動(dòng)性能得到較大提升， 完全符合預(yù)定優(yōu)化目標(biāo)。

(1)實(shí)現(xiàn)了優(yōu)化過程中對(duì)翼型幾何形狀的控制， 厚度保持在5%以上， 最大彎度低于9%， 滿足了應(yīng)用要求。

(2)采用多步翼型優(yōu)化流程進(jìn)行高效全局多目標(biāo)求解， 翼型在兩種工況下的升力系數(shù)、 升阻比都提高了20%以上。