旋翼翼型高維多目標氣動優化設計

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

（中國空氣動力研究與發展中心 計算空氣動力研究所，綿陽 621000）

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

基于數值優化方法的翼型設計技術已經取得了很大進步，但如何滿足多設計點、多目標、多約束的設計要求，仍是亟需突破的難點。對于多目標優化問題，常見的做法是采用目標加權方法將多目標統一為單目標。例如，楊慧等［3］在開展旋翼翼型多目標、多約束設計時，將懸停狀態阻力、機動狀態升力、阻力發散馬赫數特性等多設計目標加權統一為單目標問題。王清和招啟軍［4］提出了適用于中型運輸直升機旋翼翼型設計的目標及約束條件，利用遺傳算法對SC1095翼型進行了優化，也將多個設計目標組合為單一目標。

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

在進化算法研究領域，解決高維多目標優化算法難以收斂的問題也是熱點之一。Deb和Jain［8］利用聚類算子替換了NSGA-Ⅱ中的擁擠距離算子，提出了NSGA-Ⅲ算法，能夠有效應用于高維多目標優化問題。Zhang和Li［9］提出了一種基于分解的多目標優化算法（multi-objective evolutionary algorithm based on decomposition，MOEA/D），其優越的性能受到了學者們的廣泛關注。利用高維多目標優化算法，能夠充分考慮復雜氣動優化所涉及的多個目標，有助于提高綜合氣動性能。然而，該類算法對復雜氣動優化設計問題的適應性還未得到驗證。

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

為實現旋翼高維多目標氣動優化設計，提高旋翼翼型綜合氣動特性，基于MOEA/D算法建立了考慮高低速升阻特性、力矩特性、阻力發散特性等的旋翼翼型高維多目標優化設計方法，并采用高精度kriging模型以提高優化設計效率。采用SOM方法對最優解集進行了聚類分析，并對典型最優翼型進行了驗證與分析。

1 高維多目標優化設計框架

針對旋翼翼型的高維多目標優化設計問題，本文采用MOEA/D算法進行問題求解。算法在優化過程中需要計算大量的種群，對利用高精度CFD分析的優化問題，需要耗費巨大的計算資源。本文利用kriging模型，建立設計變量與氣動特性之間較為精確的模型，以代替耗時的高精度CFD分析。本節簡略介紹MOEA/D算法與kriging代理模型。

1.1 MOEA/D算法

MOEA/D算法并不直接近似Pareto前沿，而是將多目標優化問題分解成一定數量的單目標優化子問題，然后利用進化算法同時求解這些單目標子問題。這些單目標最優解集的集合是真實Pareto前沿面的一個良好近似。每個單目標子問題通過初始權向量之間的距離組成鄰域，子問題在對應的鄰域內與其他子問題進行協同優化。一般可以采用加權函數、Tchebycheff函數、基于懲罰的邊界交叉（penalty-based boundary intersection，PBI）等方法將多目標轉換為單目標求解。本文采用Tchebycheff函數聚合方法，其數學表達式為

式中：x為決策變量；z*=為參考點；為權重向量，共有N組權重向量。

采用Das-Dennis方法［17］生成指定數量的相對均勻的權重向量，它們同時滿足如下條件：

式中：H為用戶定義的正整數；權重向量個數（種群個數）滿足。

下面簡單描述MOEA/D算法步驟：

步驟1算法初始化。

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

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

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

步驟2更新。

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

修正：在y的基礎上應用特定的修正或啟發式改進策略，得到新的個體y′；若fj（y′）＜zj，則zj=fj（y′），j=1，2，…，k。

步驟3停止判斷。若滿足停止條件，算法停止，輸出Pareto前沿，否則執行步驟2。

1.2 kriging模型

kriging模型是一種基于隨機過程的估計方差最小無偏估計模型，適用于擬合具有高度非線性、多峰值問題，在氣動優化設計領域中有廣泛的應用。隨機函數Y（x）定義為

式中：fj（x）為基函數，一般可選為簡單的多項式；βj為基函數對應的系數；代表Y（x）的數學期望值；Z（x）為均值為0、方差為的靜態隨機過程。對于設計空間中的任意2個位置x與x′，隨機量的協方差可表述為

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

本文選用常用的高斯函數：

式中：θ為kriging超參數，本文采用最大似然估計方法訓練得到θ。

采用交叉驗證（leave-one-out）方法來驗證kriging模型精度，定義為

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

2 求解器精度驗證

進行旋翼翼型氣動優化設計研究前，先驗證求解器的精度。氣動特性求解采用PMB3D［18］。求解RANS方程所采用的時間推進方法為LU-SGS，空間離散采用JST格式。湍流模型為Spalart-Allmaras模型。計算網格為C型結構網格，網格示意圖如圖1所示。

圖1 OA309翼型計算網格Fig.1 Computational grids for OA309 airfoil

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

表1 OA309翼型氣動特性隨網格量的變化Table 1 Aerodynamic performance of OA309 airfoil with different amounts of computational grids

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

圖2 OA309翼型表面壓力系數分布計算值與實驗值對比（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翼型極曲線計算值與實驗值對比（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 旋翼翼型高維多目標算例

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

本文充分考慮旋翼翼型的各項性能指標，提煉出高維多目標設計優化問題，避免了設計目標難以權衡的問題。下面分別以布置在中段的9%厚度翼型和內段的12%厚度翼型為例，進行旋翼翼型高維多目標優化。

3.1 9%厚度旋翼翼型高維多目標算例

以OA309翼型為初始翼型，翼型參數化采用CST（class function/shape function transformation）方法［19］，翼型上下表面型函數階數均為6，設計變量數目為14。以OA309翼型的CST參數為基準，設計變量變化范圍為±20%。

旋翼設計同時考慮高升力能力、懸停效率和阻力發散特性，具體的設計目標和約束條件如表2所示。Ma=0.4時的最大升力系數是衡量旋翼翼型性能的重要參數，同時需要考慮Ma=0.4時的力矩特性，對應于表中的前2個設計目標。Ma=0.6時，機動狀態需要較大的最大升力系數，對應表2中設計目標3；在Ma=0.6時，為懸停狀態，為了提高懸停效率，需要減小翼型阻力，對應設計目標4。最后，考慮阻力發散特性，同時約束力矩系數的絕對值。本設計問題的設計目標個數為5，幾何約束與力矩約束共5個，涉及的計算狀態共7個。定義多目標設計問題為

表2 9%厚度旋翼翼型設計目標與約束Table 2 Design objectives and constraints for a r otor airfoil with 9% thickness

式中：對5個設計目標進行了歸一化處理，參考值為OA309翼型的性能；下標數字對應設計狀態；設計狀態5、6、7對應的馬赫數分別為0.80、0.81、0.82，Re=Ma×8×106；上標0表示基準翼型；A為翼型面積；t為翼型最大相對厚度；Cm為力矩系數。

根據設計問題，采用進化算法開展優化設計，需要大量評估翼型的氣動性能，而高精度CFD流場分析耗時耗力。為提高優化效率，采用kriging模型建立氣動特性的響應模型，替代優化過程中的流場計算。首先，采用拉丁超立方試驗設計方法（latin hypercube sampling，LHS）［20］生成700個樣本點，在7個給定的設計狀態下進行翼型氣動特性分析。進而采用交叉驗證方法檢驗模型精度。以設計狀態1為例，在700個樣本點中隨機抽取1個樣本點進行交叉驗證，重復進行10次，得到的最大預測誤差與平均預測誤差如表3所示。從表中可以看出力矩系數的預測誤差相對較大，但最大預測誤差僅為2.596%，阻力系數預測誤差保持在0.2%以下，升力系數預測更為準確。交叉驗證結果表明，kriging模型具有較高的預測精度，可以替代高精度的CFD流場分析。

表3 設計狀態1的kriging模型交叉驗證結果Table 3 Results of cross validation for kriging model under design condition 1

在建立了高精度kriging模型之后，利用kriging模型進行5目標的全局優化設計。首先，采用MOEA/D算法，參數設置采用Das-Dennis方法：H=9，得到權重向量個數為715，鄰居子問題規模為20，交叉概率為0.5，最大迭代步數為200步。

圖4給出了優化迭代200步后的Pareto前沿。從圖中可以看出，f1與f4、f1與f5存在較為明顯的非支配關系，然而f1與f2、f3的支配關系不明顯。類似圖4的二維圖像不能很好表達設計目標之間的關系，很難從中選擇出綜合性能良好的設計結果。

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

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

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

表4 Pareto前沿中選取的8組最優目標Table 4 Eight groups of optimal objectives selected from Pareto front

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

表5 OA309翼型與優化翼型幾何信息Table 5 Geometry information of OA309 and optimized airfoils

圖6 OA309翼型與優化翼型形狀Fig.6 Geometry shape of OA309 and optimized airfoils

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

圖7 OA309翼型與優化翼型升阻特性對比（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翼型與優化翼型力矩特性對比（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翼型與優化翼型升阻特性對比（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翼型與優化翼型力矩特性對比（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給出了優化翼型與OA309翼型的阻力發散特性。從圖中可以看出，在一定速度范圍內，優化翼型具有更低的阻力水平，在Ma=0.83之前，阻力系數增長更加緩慢。盡管之后阻力系數增長較快，但翼型很少工作在更高的馬赫數下。因此，優化翼型具有更好的阻力發散特性。綜上所述，優化翼型保持了低速的高升力特性，提高了高速時的阻力性能，同時提升了力矩特性。

圖11 OA309翼型與優化翼型阻力系數隨馬赫數變化特性對比（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%厚度旋翼翼型高維多目標算例

以OA312翼型為初始翼型，進行旋翼內段翼型優化，設計狀態如表6所示。設計考慮低速時的最大升力特性、高速時的翼型阻力及阻力發散馬赫數。定義設計問題如下：

表6 12%厚度旋翼翼型設計目標與約束Table 6 Design objectives and constraints for a rotor airfoil with 12% thickness

式中：下標數字對應設計狀態；設計狀態5、6、7對應的馬赫數分別為0.74、0.75、0.76；Re=Ma×8×106；其余符號含義同3.1節。

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

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

表7 Pareto前沿中選取的5組最優目標Table 7 Five gr oups of optimal objectives selected from Pareto front

圖13為優化翼型與OA312翼型的幾何形狀對比，表8中為翼型的幾何信息，包括最大相對厚度與翼型面積，從表中看出，優化翼型的相對厚度與無量綱面積均滿足約束條件。

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

圖13 OA312翼型與優化翼型形狀Fig.13 Geometry shape of OA312 and optimized airfoils

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

圖14 OA312翼型與優化翼型升阻比對比（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翼型與優化翼型力矩特性對比（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翼型與優化翼型升阻特性對比（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翼型與優化翼型力矩特性對比（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給出了優化翼型與OA312翼型的阻力發散特性。從圖中可以明顯看出，優化翼型在較寬的馬赫數范圍內保持著更低的阻力系數，其中在Ma=0.76處減阻效果為6.1%。

圖18 OA312翼型與優化翼型阻力系數隨馬赫數變化特性對比（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）

綜上，旋翼內段翼型的優化取得了良好效果，在保持低速高升力特性的同時，明顯提高了高速時的阻力性能，同時提升了力矩特性。

4 結論

基于MOEA/D算法，本文提出了旋翼翼型高維多目標優化設計方法，采用kriging代理模型提高了求解效率，并應用SOM方法對最優Pareto解集進行了可視化分析。旋翼翼型5目標優化設計的結果表明：

1）MOEA/D算法能夠有效求解高維多目標氣動優化設計問題；結合高精度代理模型，能夠高效解決旋翼翼型復雜氣動優化設計問題。

2）高維多目標優化設計結果難以通過二維圖表直觀地展示最優Pareto前沿。SOM 聚類分析能夠為設計者的決策提供較好的全局可視化結果，有助于從Pareto前沿中選取綜合性能良好的設計結果。

3）旋翼翼型高維多目標設計避免了設計目標的加權或折中處理，利用高維多目標進化算法實現了翼型阻力特性、力矩特性和阻力發散特性的綜合提升，能夠提高旋翼懸停、機動性能和操縱特性。