宋超，周鑄，李偉斌，羅驍

（中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000）

旋翼是直升機(jī)的核心氣動(dòng)部件，旋翼氣動(dòng)設(shè)計(jì)技術(shù)為直升機(jī)平臺(tái)性能提供堅(jiān)實(shí)的支撐［1］。旋翼翼型是旋翼氣動(dòng)設(shè)計(jì)的基礎(chǔ)，直接影響直升機(jī)的飛行包線、載荷、噪聲等關(guān)鍵指標(biāo)，其復(fù)雜的工作狀態(tài)決定了翼型設(shè)計(jì)是一項(xiàng)十分復(fù)雜的工作。旋翼槳葉之間存在強(qiáng)烈的氣動(dòng)干擾，且展向剖面翼型相對(duì)氣流速度受高速旋轉(zhuǎn)而差別巨大。在直升機(jī)飛行的不同工況下，前飛速度與旋轉(zhuǎn)速度疊加，使得旋翼在不同方位角的工況也劇烈變化。例如，在前飛時(shí)，前行槳葉槳尖區(qū)域上存在跨聲速流動(dòng)，伴隨激波的產(chǎn)生；在機(jī)動(dòng)飛行時(shí)，旋翼后行槳葉上存在分離或動(dòng)態(tài)失速現(xiàn)象。因此，先進(jìn)的旋翼翼型要求兼顧高速低阻、低速高升力、低零升力矩等特性［2］，其設(shè)計(jì)是典型的多設(shè)計(jì)點(diǎn)、多目標(biāo)和強(qiáng)約束優(yōu)化問(wèn)題。

基于數(shù)值優(yōu)化方法的翼型設(shè)計(jì)技術(shù)已經(jīng)取得了很大進(jìn)步，但如何滿(mǎn)足多設(shè)計(jì)點(diǎn)、多目標(biāo)、多約束的設(shè)計(jì)要求，仍是亟需突破的難點(diǎn)。對(duì)于多目標(biāo)優(yōu)化問(wèn)題，常見(jiàn)的做法是采用目標(biāo)加權(quán)方法將多目標(biāo)統(tǒng)一為單目標(biāo)。例如，楊慧等［3］在開(kāi)展旋翼翼型多目標(biāo)、多約束設(shè)計(jì)時(shí)，將懸停狀態(tài)阻力、機(jī)動(dòng)狀態(tài)升力、阻力發(fā)散馬赫數(shù)特性等多設(shè)計(jì)目標(biāo)加權(quán)統(tǒng)一為單目標(biāo)問(wèn)題。王清和招啟軍［4］提出了適用于中型運(yùn)輸直升機(jī)旋翼翼型設(shè)計(jì)的目標(biāo)及約束條件，利用遺傳算法對(duì)SC1095翼型進(jìn)行了優(yōu)化，也將多個(gè)設(shè)計(jì)目標(biāo)組合為單一目標(biāo)。

目標(biāo)加權(quán)方法的優(yōu)化結(jié)果依賴(lài)于權(quán)因子的選取，當(dāng)優(yōu)化目標(biāo)較多時(shí)，選取合適的權(quán)因子是十分困難的。采用基于Pareto非支配關(guān)系的多目標(biāo)進(jìn)化算法，如NSGA-Ⅱ，可以得到整個(gè)Pareto最優(yōu)前沿，以供設(shè)計(jì)者決策。Wang和Zhao［5］采用NSGA-Ⅱ算法進(jìn)行了考慮2個(gè)設(shè)計(jì)點(diǎn)的旋翼翼型設(shè)計(jì)，旨在提高翼型升阻特性。Massaro和Benini［6］以SC1095翼型為基準(zhǔn)，利用進(jìn)化算法與代理模型進(jìn)行了旋翼翼型兩目標(biāo)的設(shè)計(jì)，得到了最優(yōu)Pareto前沿。然而，進(jìn)化多目標(biāo)算法解決高維多目標(biāo)設(shè)計(jì)問(wèn)題（目標(biāo)數(shù)大于3）時(shí)面臨難以收斂的挑戰(zhàn)。具體而言，隨著目標(biāo)數(shù)的增多，Pareto占優(yōu)關(guān)系難以區(qū)分進(jìn)化個(gè)體，進(jìn)化動(dòng)力不足，描述Pareto前沿的解集數(shù)目也呈指數(shù)級(jí)增長(zhǎng)，算法收斂緩慢，甚至難以收斂。為了避免高維多目標(biāo)優(yōu)化問(wèn)題面臨的收斂困難問(wèn)題，Zhao等［7］利用主成分分析（principal component analysis，PCA）降維方法，將6個(gè)目標(biāo)的旋翼設(shè)計(jì)問(wèn)題轉(zhuǎn)化為2個(gè)目標(biāo)設(shè)計(jì)問(wèn)題。然而，多目標(biāo)自身相關(guān)的嚴(yán)苛要求和降維對(duì)結(jié)果的不確定影響，使得降維方法的應(yīng)用受到了一定限制。

在進(jìn)化算法研究領(lǐng)域，解決高維多目標(biāo)優(yōu)化算法難以收斂的問(wèn)題也是熱點(diǎn)之一。Deb和Jain［8］利用聚類(lèi)算子替換了NSGA-Ⅱ中的擁擠距離算子，提出了NSGA-Ⅲ算法，能夠有效應(yīng)用于高維多目標(biāo)優(yōu)化問(wèn)題。Zhang和Li［9］提出了一種基于分解的多目標(biāo)優(yōu)化算法（multi-objective evolutionary algorithm based on decomposition，MOEA/D），其優(yōu)越的性能受到了學(xué)者們的廣泛關(guān)注。利用高維多目標(biāo)優(yōu)化算法，能夠充分考慮復(fù)雜氣動(dòng)優(yōu)化所涉及的多個(gè)目標(biāo)，有助于提高綜合氣動(dòng)性能。然而，該類(lèi)算法對(duì)復(fù)雜氣動(dòng)優(yōu)化設(shè)計(jì)問(wèn)題的適應(yīng)性還未得到驗(yàn)證。

高維多目標(biāo)優(yōu)化設(shè)計(jì)面臨的另一個(gè)挑戰(zhàn)是目標(biāo)空間的可視化問(wèn)題。由于目標(biāo)空間維數(shù)的增加，Pareto前沿難以在二維或三維坐標(biāo)系中直觀展現(xiàn)。高維目標(biāo)空間的可視化也是亟待解決的問(wèn)題與研究的熱點(diǎn)。高維目標(biāo)空間可視化技術(shù)可分為2類(lèi)：①在平行坐標(biāo)系中表示目標(biāo)的解，如熱圖［10］和平行坐標(biāo)系［11］；②構(gòu)建高維到二維的映射關(guān)系，并保留解之間的距離信息，如自組織圖映射（self-organizing mapping，SOM）［12］、徑向坐標(biāo)可視化［13］、雷達(dá)圖［14］等。SOM方法憑借其直觀的解集性質(zhì)表達(dá)形式，成為了目前流行的可視化方法。Kumano等［15］采用多目標(biāo)遺傳算法進(jìn)行了涉及4個(gè)目標(biāo)的公務(wù)機(jī)機(jī)翼氣動(dòng)/結(jié)構(gòu)多學(xué)科優(yōu)化，并利用SOM 對(duì)最優(yōu)Pareto前沿進(jìn)行了可視化分析。王超［16］利用SOM對(duì)旋翼翼型與戰(zhàn)斗機(jī)翼型最優(yōu)Pareto前沿進(jìn)行了聚類(lèi)分析，為設(shè)計(jì)決策提供了全局視角。

為實(shí)現(xiàn)旋翼高維多目標(biāo)氣動(dòng)優(yōu)化設(shè)計(jì)，提高旋翼翼型綜合氣動(dòng)特性，基于MOEA/D算法建立了考慮高低速升阻特性、力矩特性、阻力發(fā)散特性等的旋翼翼型高維多目標(biāo)優(yōu)化設(shè)計(jì)方法，并采用高精度kriging模型以提高優(yōu)化設(shè)計(jì)效率。采用SOM方法對(duì)最優(yōu)解集進(jìn)行了聚類(lèi)分析，并對(duì)典型最優(yōu)翼型進(jìn)行了驗(yàn)證與分析。

1 高維多目標(biāo)優(yōu)化設(shè)計(jì)框架

針對(duì)旋翼翼型的高維多目標(biāo)優(yōu)化設(shè)計(jì)問(wèn)題，本文采用MOEA/D算法進(jìn)行問(wèn)題求解。算法在優(yōu)化過(guò)程中需要計(jì)算大量的種群，對(duì)利用高精度CFD分析的優(yōu)化問(wèn)題，需要耗費(fèi)巨大的計(jì)算資源。本文利用kriging模型，建立設(shè)計(jì)變量與氣動(dòng)特性之間較為精確的模型，以代替耗時(shí)的高精度CFD分析。本節(jié)簡(jiǎn)略介紹MOEA/D算法與kriging代理模型。

1.1 MOEA/D算法

MOEA/D算法并不直接近似Pareto前沿，而是將多目標(biāo)優(yōu)化問(wèn)題分解成一定數(shù)量的單目標(biāo)優(yōu)化子問(wèn)題，然后利用進(jìn)化算法同時(shí)求解這些單目標(biāo)子問(wèn)題。這些單目標(biāo)最優(yōu)解集的集合是真實(shí)Pareto前沿面的一個(gè)良好近似。每個(gè)單目標(biāo)子問(wèn)題通過(guò)初始權(quán)向量之間的距離組成鄰域，子問(wèn)題在對(duì)應(yīng)的鄰域內(nèi)與其他子問(wèn)題進(jìn)行協(xié)同優(yōu)化。一般可以采用加權(quán)函數(shù)、Tchebycheff函數(shù)、基于懲罰的邊界交叉（penalty-based boundary intersection，PBI）等方法將多目標(biāo)轉(zhuǎn)換為單目標(biāo)求解。本文采用Tchebycheff函數(shù)聚合方法，其數(shù)學(xué)表達(dá)式為

式中：x為決策變量；z*=為參考點(diǎn)；為權(quán)重向量，共有N組權(quán)重向量。

采用Das-Dennis方法［17］生成指定數(shù)量的相對(duì)均勻的權(quán)重向量，它們同時(shí)滿(mǎn)足如下條件：

式中：H為用戶(hù)定義的正整數(shù)；權(quán)重向量個(gè)數(shù)（種群個(gè)數(shù)）滿(mǎn)足。

下面簡(jiǎn)單描述MOEA/D算法步驟：

步驟1算法初始化。

構(gòu)造權(quán)重向量λ1，λ2，…，λN，并計(jì)算任意2個(gè)權(quán)重向量之間的歐氏距離，為每個(gè)權(quán)重向量選出最近的T個(gè)向量作為它的鄰域B（i）=｛i1，i2，…，iT｝。

初始化種群x1，x2，…，xN，設(shè)FVi=F（xi），i=1，2，…，N。

初始化z=（z1，z2，…，zk），zj=），其中j=1，2，…，k。

步驟2更新。

繁殖：從B（i）中選擇2個(gè)個(gè)體r1、r2，利用xr1、xr2交叉生成新的個(gè)體ˉy，再以一定的變異概率生成y。

修正：在y的基礎(chǔ)上應(yīng)用特定的修正或啟發(fā)式改進(jìn)策略，得到新的個(gè)體y′；若fj（y′）＜zj，則zj=fj（y′），j=1，2，…，k。

步驟3停止判斷。若滿(mǎn)足停止條件，算法停止，輸出Pareto前沿，否則執(zhí)行步驟2。

1.2 kriging模型

kriging模型是一種基于隨機(jī)過(guò)程的估計(jì)方差最小無(wú)偏估計(jì)模型，適用于擬合具有高度非線性、多峰值問(wèn)題，在氣動(dòng)優(yōu)化設(shè)計(jì)領(lǐng)域中有廣泛的應(yīng)用。隨機(jī)函數(shù)Y（x）定義為

式中：fj（x）為基函數(shù)，一般可選為簡(jiǎn)單的多項(xiàng)式；βj為基函數(shù)對(duì)應(yīng)的系數(shù)；代表Y（x）的數(shù)學(xué)期望值；Z（x）為均值為0、方差為的靜態(tài)隨機(jī)過(guò)程。對(duì)于設(shè)計(jì)空間中的任意2個(gè)位置x與x′，隨機(jī)量的協(xié)方差可表述為

式中：相關(guān)函數(shù)R（x，x′）只與空間距離有關(guān)，代表了不同位置隨機(jī)變量之間的相關(guān)性，相關(guān)函數(shù)在兩位置距離無(wú)窮大時(shí)R=0，距離為零時(shí)R=1；R隨距離的增大而減小。

本文選用常用的高斯函數(shù)：

式中：θ為kriging超參數(shù)，本文采用最大似然估計(jì)方法訓(xùn)練得到θ。

采用交叉驗(yàn)證（leave-one-out）方法來(lái)驗(yàn)證kriging模型精度，定義為

式中：y（xi）為第i個(gè)樣本點(diǎn)的響應(yīng)值；y-（i）（xi）為xi的預(yù)測(cè)值，特殊之處在于該預(yù)測(cè)值是由利用排除第i個(gè)樣本點(diǎn)之后的樣本集所建立模型的預(yù)測(cè)值。

2 求解器精度驗(yàn)證

進(jìn)行旋翼翼型氣動(dòng)優(yōu)化設(shè)計(jì)研究前，先驗(yàn)證求解器的精度。氣動(dòng)特性求解采用PMB3D［18］。求解RANS方程所采用的時(shí)間推進(jìn)方法為L(zhǎng)U-SGS，空間離散采用JST格式。湍流模型為Spalart-Allmaras模型。計(jì)算網(wǎng)格為C型結(jié)構(gòu)網(wǎng)格，網(wǎng)格示意圖如圖1所示。

圖1 OA309翼型計(jì)算網(wǎng)格Fig.1 Computational grids for OA309 airfoil

為排除計(jì)算網(wǎng)格對(duì)結(jié)果的影響，開(kāi)展網(wǎng)格無(wú)關(guān)性檢驗(yàn)，選取3套網(wǎng)格作為檢驗(yàn)對(duì)象。表1為不同網(wǎng)格量對(duì)應(yīng)的升力系數(shù)Cl和阻力系數(shù)Cd。計(jì)算狀態(tài)為Ma=0.6，Re=2.5×106，Cl=0.50。從表中可以看出，3組不同網(wǎng)格量計(jì)算得到的升力系數(shù)幾乎一致，阻力系數(shù)隨網(wǎng)格量的增大而減小，其中網(wǎng)格2與網(wǎng)格3的阻力系數(shù)僅相差0.67%。因此，兼顧網(wǎng)格規(guī)模和計(jì)算精度考慮，選用網(wǎng)格2進(jìn)行流場(chǎng)分析。

表1 OA309翼型氣動(dòng)特性隨網(wǎng)格量的變化Table 1 Aerodynamic performance of OA309 airfoil with different amounts of computational grids

為驗(yàn)證求解器的可靠性，進(jìn)行了OA309翼型不同狀態(tài)下計(jì)算值與實(shí)驗(yàn)值的對(duì)比。實(shí)驗(yàn)值來(lái)自西北工業(yè)大學(xué)翼型葉柵空氣動(dòng)力學(xué)國(guó)防科技重點(diǎn)實(shí)驗(yàn)室NF-6風(fēng)洞，實(shí)驗(yàn)雷諾數(shù)表示為Re=Ma×4.2×106。圖2為OA309翼型表面壓力分布計(jì)算值與實(shí)驗(yàn)值對(duì)比（Ma=0.6，Re=2.5×106，Cl=0.6），計(jì)算值在壓力峰值處略小于實(shí)驗(yàn)值，兩者總體吻合良好（Cp為表面壓力系數(shù)）。圖3為相同狀態(tài)下OA309翼型極曲線的計(jì)算值和實(shí)驗(yàn)值對(duì)比，可以看出兩者的阻力值吻合良好。在Cl=0.5時(shí)，阻力計(jì)算值相比實(shí)驗(yàn)值的計(jì)算誤差約為6.4%。由此可以看出，PMB3D計(jì)算結(jié)果能夠很好反映翼型的氣動(dòng)特性。

圖2 OA309翼型表面壓力系數(shù)分布計(jì)算值與實(shí)驗(yàn)值對(duì)比（Ma=0.6，Re=2.5×106，Cl=0.6）Fig.2 Comparison of pressure coefficient at surface between numerical results and experimental data for OA309 airfoil（Ma=0.6，Re=2.5×106，Cl=0.6）

圖3 OA309翼型極曲線計(jì)算值與實(shí)驗(yàn)值對(duì)比（Ma=0.6，Re=2.5×106）Fig.3 Comparison of polars between numerical results and experimental data for OA309 airfoil（Ma=0.6，Re=2.5×106）

3 旋翼翼型高維多目標(biāo)算例

旋翼不同站位的翼型流場(chǎng)環(huán)境隨旋翼工作狀態(tài)和旋翼方位角的變化而劇烈變化，因此，旋翼翼型的設(shè)計(jì)需要充分考慮其站位的不同。一般將旋翼翼型劃分為外段翼型（90%R～100%R，R為旋翼半徑）、中段翼型（80%R～90%R）和內(nèi)段翼型（75%R～80%R）。外段翼型強(qiáng)調(diào)高速性能；內(nèi)段翼型強(qiáng)調(diào)高升力性能；中段翼型則兼顧兩方面性能。然而，旋翼翼型設(shè)計(jì)指標(biāo)之間往往相互矛盾、相互制約，需要根據(jù)直升機(jī)性能需求確定設(shè)計(jì)目標(biāo)，通過(guò)仔細(xì)地權(quán)衡和折中，在上述要求中取得一個(gè)較好的平衡。

本文充分考慮旋翼翼型的各項(xiàng)性能指標(biāo)，提煉出高維多目標(biāo)設(shè)計(jì)優(yōu)化問(wèn)題，避免了設(shè)計(jì)目標(biāo)難以權(quán)衡的問(wèn)題。下面分別以布置在中段的9%厚度翼型和內(nèi)段的12%厚度翼型為例，進(jìn)行旋翼翼型高維多目標(biāo)優(yōu)化。

3.1 9%厚度旋翼翼型高維多目標(biāo)算例

以O(shè)A309翼型為初始翼型，翼型參數(shù)化采用CST（class function/shape function transformation）方法［19］，翼型上下表面型函數(shù)階數(shù)均為6，設(shè)計(jì)變量數(shù)目為14。以O(shè)A309翼型的CST參數(shù)為基準(zhǔn)，設(shè)計(jì)變量變化范圍為±20%。

旋翼設(shè)計(jì)同時(shí)考慮高升力能力、懸停效率和阻力發(fā)散特性，具體的設(shè)計(jì)目標(biāo)和約束條件如表2所示。Ma=0.4時(shí)的最大升力系數(shù)是衡量旋翼翼型性能的重要參數(shù)，同時(shí)需要考慮Ma=0.4時(shí)的力矩特性，對(duì)應(yīng)于表中的前2個(gè)設(shè)計(jì)目標(biāo)。Ma=0.6時(shí)，機(jī)動(dòng)狀態(tài)需要較大的最大升力系數(shù)，對(duì)應(yīng)表2中設(shè)計(jì)目標(biāo)3；在Ma=0.6時(shí)，為懸停狀態(tài)，為了提高懸停效率，需要減小翼型阻力，對(duì)應(yīng)設(shè)計(jì)目標(biāo)4。最后，考慮阻力發(fā)散特性，同時(shí)約束力矩系數(shù)的絕對(duì)值。本設(shè)計(jì)問(wèn)題的設(shè)計(jì)目標(biāo)個(gè)數(shù)為5，幾何約束與力矩約束共5個(gè)，涉及的計(jì)算狀態(tài)共7個(gè)。定義多目標(biāo)設(shè)計(jì)問(wèn)題為

表2 9%厚度旋翼翼型設(shè)計(jì)目標(biāo)與約束Table 2 Design objectives and constraints for a r otor airfoil with 9% thickness

式中：對(duì)5個(gè)設(shè)計(jì)目標(biāo)進(jìn)行了歸一化處理，參考值為OA309翼型的性能；下標(biāo)數(shù)字對(duì)應(yīng)設(shè)計(jì)狀態(tài)；設(shè)計(jì)狀態(tài)5、6、7對(duì)應(yīng)的馬赫數(shù)分別為0.80、0.81、0.82，Re=Ma×8×106；上標(biāo)0表示基準(zhǔn)翼型；A為翼型面積；t為翼型最大相對(duì)厚度；Cm為力矩系數(shù)。

根據(jù)設(shè)計(jì)問(wèn)題，采用進(jìn)化算法開(kāi)展優(yōu)化設(shè)計(jì)，需要大量評(píng)估翼型的氣動(dòng)性能，而高精度CFD流場(chǎng)分析耗時(shí)耗力。為提高優(yōu)化效率，采用kriging模型建立氣動(dòng)特性的響應(yīng)模型，替代優(yōu)化過(guò)程中的流場(chǎng)計(jì)算。首先，采用拉丁超立方試驗(yàn)設(shè)計(jì)方法（latin hypercube sampling，LHS）［20］生成700個(gè)樣本點(diǎn)，在7個(gè)給定的設(shè)計(jì)狀態(tài)下進(jìn)行翼型氣動(dòng)特性分析。進(jìn)而采用交叉驗(yàn)證方法檢驗(yàn)?zāi)Ｐ途?。以設(shè)計(jì)狀態(tài)1為例，在700個(gè)樣本點(diǎn)中隨機(jī)抽取1個(gè)樣本點(diǎn)進(jìn)行交叉驗(yàn)證，重復(fù)進(jìn)行10次，得到的最大預(yù)測(cè)誤差與平均預(yù)測(cè)誤差如表3所示。從表中可以看出力矩系數(shù)的預(yù)測(cè)誤差相對(duì)較大，但最大預(yù)測(cè)誤差僅為2.596%，阻力系數(shù)預(yù)測(cè)誤差保持在0.2%以下，升力系數(shù)預(yù)測(cè)更為準(zhǔn)確。交叉驗(yàn)證結(jié)果表明，kriging模型具有較高的預(yù)測(cè)精度，可以替代高精度的CFD流場(chǎng)分析。

表3 設(shè)計(jì)狀態(tài)1的kriging模型交叉驗(yàn)證結(jié)果Table 3 Results of cross validation for kriging model under design condition 1

在建立了高精度kriging模型之后，利用kriging模型進(jìn)行5目標(biāo)的全局優(yōu)化設(shè)計(jì)。首先，采用MOEA/D算法，參數(shù)設(shè)置采用Das-Dennis方法：H=9，得到權(quán)重向量個(gè)數(shù)為715，鄰居子問(wèn)題規(guī)模為20，交叉概率為0.5，最大迭代步數(shù)為200步。

圖4給出了優(yōu)化迭代200步后的Pareto前沿。從圖中可以看出，f1與f4、f1與f5存在較為明顯的非支配關(guān)系，然而f1與f2、f3的支配關(guān)系不明顯。類(lèi)似圖4的二維圖像不能很好表達(dá)設(shè)計(jì)目標(biāo)之間的關(guān)系，很難從中選擇出綜合性能良好的設(shè)計(jì)結(jié)果。

圖4 旋翼翼型優(yōu)化最優(yōu)Pareto前沿Fig.4 Optimal Pareto front for optimization of rotor airfoil

本文采用SOM處理最優(yōu)Pareto前沿，選取的SOM單元個(gè)數(shù)為6×6。SOM 聚類(lèi)分析將每組目標(biāo)視為一個(gè)個(gè)體，根據(jù)個(gè)體之間的距離將715組5維目標(biāo)映射在36個(gè)單元上，處理后的結(jié)果如圖5所示。比較圖5（a）與圖5（d），圖5（a）中深色區(qū)域基本對(duì)應(yīng)圖5（b）中淺色區(qū)域，反映出目標(biāo)1與目標(biāo)4是相互矛盾的，這也與圖4中反映的趨勢(shì)一致，說(shuō)明了SOM 聚類(lèi)分析的正確性。圖5（f）為每個(gè)SOM 單元包含的個(gè)體數(shù)，綜合圖5（a）～（e），折中選取圖5（f）亮黃色單元中的個(gè)體進(jìn)行分析。表4中給出了亮黃色單元中的8組目標(biāo)，可以看出8組最優(yōu)目標(biāo)具有高度相似性，其中目標(biāo)1和目標(biāo)4與基本翼型基本保持一致，目標(biāo)3略有下降，目標(biāo)2和目標(biāo)5提升明顯。下面選取其中第1個(gè)個(gè)體，并利用CFD進(jìn)行氣動(dòng)特性分析。

圖5 最優(yōu)Pareto前沿SOM可視化結(jié)果（9%厚度翼型）Fig.5 Visualization results of optimal Pareto front using SOM（airfoil with 9% thickness）

表4 Pareto前沿中選取的8組最優(yōu)目標(biāo)Table 4 Eight groups of optimal objectives selected from Pareto front

圖6為選取的優(yōu)化翼型形狀，表5為優(yōu)化翼型與OA309翼型的幾何信息，包括最大相對(duì)厚度與翼型面積，從表中可以看出，優(yōu)化翼型的相對(duì)厚度與無(wú)量綱面積均滿(mǎn)足約束條件。

表5 OA309翼型與優(yōu)化翼型幾何信息Table 5 Geometry information of OA309 and optimized airfoils

圖6 OA309翼型與優(yōu)化翼型形狀Fig.6 Geometry shape of OA309 and optimized airfoils

優(yōu)化翼型與OA309翼型的升阻特性對(duì)比如圖7所示，可以看出，優(yōu)化翼型與OA309翼型在Ma=0.4時(shí)的升阻特性基本一致。圖8給出了Ma=0.4時(shí)的力矩特性?？梢钥闯?，優(yōu)化翼型的零升力矩絕對(duì)值相比OA309翼型基本不變，在更大的升力系數(shù)范圍內(nèi)，優(yōu)化翼型相比OA309翼型的力矩系數(shù)更接近0，力矩系數(shù)幅值（Cl≈0.8）減小約50.7%，力矩特性得到改善。在高馬赫數(shù)Ma=0.6時(shí)，氣動(dòng)特性可以由圖9、圖10中看出，高馬赫數(shù)下的最大升力系數(shù)和最大升阻比均有明顯提高，其中最大升力系數(shù)提高約6.5%，最大升阻比提高約7.7%。對(duì)于力矩特性，優(yōu)化翼型在較大范圍具有更小的力矩幅值，力矩特性得到提升。

圖7 OA309翼型與優(yōu)化翼型升阻特性對(duì)比（Ma=0.4，Re=3.2×106）Fig.7 Comparison of lift drag ratio curves between OA309 and optimized airfoils（Ma=0.4，Re=3.2×106）

圖8 OA309翼型與優(yōu)化翼型力矩特性對(duì)比（Ma=0.4，Re=3.2×106）Fig.8 Comparison of moment coefficient curves between OA309 and optimized airfoils（Ma=0.4，Re=3.2×106）

圖9 OA309翼型與優(yōu)化翼型升阻特性對(duì)比（Ma=0.6，Re=4.8×106）Fig.9 Comparison of lift drag ratio curves between OA309 and optimized airfoils（Ma=0.6，Re=4.8×106）

圖10 OA309翼型與優(yōu)化翼型力矩特性對(duì)比（Ma=0.6，Re=4.8×106）Fig.10 Comparison of moment coefficient curves between OA309 and optimized airfoil（Ma=0.6，Re=4.8×106）

圖11給出了優(yōu)化翼型與OA309翼型的阻力發(fā)散特性。從圖中可以看出，在一定速度范圍內(nèi)，優(yōu)化翼型具有更低的阻力水平，在Ma=0.83之前，阻力系數(shù)增長(zhǎng)更加緩慢。盡管之后阻力系數(shù)增長(zhǎng)較快，但翼型很少工作在更高的馬赫數(shù)下。因此，優(yōu)化翼型具有更好的阻力發(fā)散特性。綜上所述，優(yōu)化翼型保持了低速的高升力特性，提高了高速時(shí)的阻力性能，同時(shí)提升了力矩特性。

圖11 OA309翼型與優(yōu)化翼型阻力系數(shù)隨馬赫數(shù)變化特性對(duì)比（Cl=0，Re=Ma×8×106）Fig.11 Comparison of drag coefficient varying with Mach number between OA309 and optimized airfoils（Cl=0，Re=Ma×8×106）

3.2 12%厚度旋翼翼型高維多目標(biāo)算例

以O(shè)A312翼型為初始翼型，進(jìn)行旋翼內(nèi)段翼型優(yōu)化，設(shè)計(jì)狀態(tài)如表6所示。設(shè)計(jì)考慮低速時(shí)的最大升力特性、高速時(shí)的翼型阻力及阻力發(fā)散馬赫數(shù)。定義設(shè)計(jì)問(wèn)題如下：

表6 12%厚度旋翼翼型設(shè)計(jì)目標(biāo)與約束Table 6 Design objectives and constraints for a rotor airfoil with 12% thickness

式中：下標(biāo)數(shù)字對(duì)應(yīng)設(shè)計(jì)狀態(tài)；設(shè)計(jì)狀態(tài)5、6、7對(duì)應(yīng)的馬赫數(shù)分別為0.74、0.75、0.76；Re=Ma×8×106；其余符號(hào)含義同3.1節(jié)。

翼型參數(shù)化方法、設(shè)計(jì)變量選取和優(yōu)化算法參數(shù)設(shè)置與3.1節(jié)相同。采用SOM方法進(jìn)行優(yōu)化結(jié)果后的處理。SOM 可視化結(jié)果如圖12所示。從圖中明顯看出，目標(biāo)1與目標(biāo)3、4是相互矛盾的，而目標(biāo)1與目標(biāo)5則有較好的正相關(guān)關(guān)系。OA312翼型具有較高的低速最大升力系數(shù)，因此，在確保升力系數(shù)不損失的情況下，盡可能降低高速工況阻力。最終，選取的一組解在圖12（f）中標(biāo)記出來(lái)，相應(yīng)的目標(biāo)值如表7所示。選取表中第2個(gè)最優(yōu)解，利用CFD進(jìn)行詳細(xì)性能分析。

圖12 最優(yōu)Pareto前沿SOM可視化結(jié)果（12%厚度翼型）Fig.12 Visualization results of optimal Pareto front using SOM（airfoil with 12% thickness）

表7 Pareto前沿中選取的5組最優(yōu)目標(biāo)Table 7 Five gr oups of optimal objectives selected from Pareto front

圖13為優(yōu)化翼型與OA312翼型的幾何形狀對(duì)比，表8中為翼型的幾何信息，包括最大相對(duì)厚度與翼型面積，從表中看出，優(yōu)化翼型的相對(duì)厚度與無(wú)量綱面積均滿(mǎn)足約束條件。

表8 OA312翼型與優(yōu)化翼型幾何信息Table 8 Geometr y information of OA312 and optimized airfoils

圖13 OA312翼型與優(yōu)化翼型形狀Fig.13 Geometry shape of OA312 and optimized airfoils

圖14給出了翼型在Ma=0.4時(shí)的升阻特性。從圖中看出，在升力線性段，兩者無(wú)明顯差別，優(yōu)化翼型的最大升力系數(shù)略有下降。圖15給出了Ma=0.4時(shí)的翼型力矩特性，優(yōu)化翼型相比OA312翼型的力矩特性有明顯改善。在Ma=0.6時(shí)，優(yōu)化翼型的升阻特性略?xún)?yōu)于OA312翼型（見(jiàn)圖16）。同時(shí)力矩特性改善明顯，如圖17所示，的最大值減小約11%。

圖14 OA312翼型與優(yōu)化翼型升阻比對(duì)比（Ma=0.4，Re=3.2×106）Fig.14 Comparison of lift-to-drag ratio curves between OA312 and optimized airfoils（Ma=0.4，Re=3.2×106）

圖15 OA312翼型與優(yōu)化翼型力矩特性對(duì)比（Ma=0.4，Re=3.2×106）Fig.15 Comparison of moment coefficient curves between OA312 and optimized airfoils（Ma=0.4，Re=3.2×106）

圖16 OA312翼型與優(yōu)化翼型升阻特性對(duì)比（Ma=0.6，Re=4.8×106）Fig.16 Comparison of lift-to-drag ratio curves between OA312 and optimized airfoils（Ma=0.6，Re=4.8×106）

圖17 OA312翼型與優(yōu)化翼型力矩特性對(duì)比（Ma=0.6，Re=4.8×106）Fig.17 Comparison of moment coefficient curves between OA312 and optimized airfoils（Ma=0.6，Re=4.8×106）

圖18給出了優(yōu)化翼型與OA312翼型的阻力發(fā)散特性。從圖中可以明顯看出，優(yōu)化翼型在較寬的馬赫數(shù)范圍內(nèi)保持著更低的阻力系數(shù)，其中在Ma=0.76處減阻效果為6.1%。

圖18 OA312翼型與優(yōu)化翼型阻力系數(shù)隨馬赫數(shù)變化特性對(duì)比（Cl=0，Re=Ma×8×106）Fig.18 Comparison of drag coefficient varying with Mach number between OA312 and optimized airfoils（Cl=0，Re=Ma×8×106）

綜上，旋翼內(nèi)段翼型的優(yōu)化取得了良好效果，在保持低速高升力特性的同時(shí)，明顯提高了高速時(shí)的阻力性能，同時(shí)提升了力矩特性。

4 結(jié)論

基于MOEA/D算法，本文提出了旋翼翼型高維多目標(biāo)優(yōu)化設(shè)計(jì)方法，采用kriging代理模型提高了求解效率，并應(yīng)用SOM方法對(duì)最優(yōu)Pareto解集進(jìn)行了可視化分析。旋翼翼型5目標(biāo)優(yōu)化設(shè)計(jì)的結(jié)果表明：

1）MOEA/D算法能夠有效求解高維多目標(biāo)氣動(dòng)優(yōu)化設(shè)計(jì)問(wèn)題；結(jié)合高精度代理模型，能夠高效解決旋翼翼型復(fù)雜氣動(dòng)優(yōu)化設(shè)計(jì)問(wèn)題。

2）高維多目標(biāo)優(yōu)化設(shè)計(jì)結(jié)果難以通過(guò)二維圖表直觀地展示最優(yōu)Pareto前沿。SOM 聚類(lèi)分析能夠?yàn)樵O(shè)計(jì)者的決策提供較好的全局可視化結(jié)果，有助于從Pareto前沿中選取綜合性能良好的設(shè)計(jì)結(jié)果。

3）旋翼翼型高維多目標(biāo)設(shè)計(jì)避免了設(shè)計(jì)目標(biāo)的加權(quán)或折中處理，利用高維多目標(biāo)進(jìn)化算法實(shí)現(xiàn)了翼型阻力特性、力矩特性和阻力發(fā)散特性的綜合提升，能夠提高旋翼懸停、機(jī)動(dòng)性能和操縱特性。