翼型氣動優化設計問題多極值特性研究

何萌，白俊強，昌敏，張煜

(1.西北工業大學 航空學院，西安 710072) (2.西北工業大學 無人系統技術研究院，西安 710072)

0 引 言

飛行器氣動外形優化研究已取得了顯著的成果，通過采用優化算法，能夠設計出高性能的翼型，甚至能進行復雜全機構型的設計。影響飛行器氣動外形優化效率[1]的兩個主要因素為計算流體力學(Computational Fluid Dynamics，簡稱CFD)和優化算法，目前CFD作為一項成熟的技術已廣泛的應用于工程應用當中，而優化算法仍然有很大的發展空間。優化算法的合理選擇依賴于優化問題，必須考慮設計變量的類型(例如離散和連續)、約束的數量、設計空間的特點(一般來說，對于大部分優化問題，設計空間是光滑的)。

目前常用的優化算法[2]可以分為兩類：一類是基于梯度的優化算法，這類算法需要目標函數的梯度信息，例如最速下降法、序列二次規劃算法(Sequential Quadratic Programming，簡稱SQP)等；另一類是全局優化算法，這類算法只需要目標函數值，而不需要目標函數的梯度信息，例如粒子群算法、遺傳算法以及直接搜索算法等。全局算法耗時長，理論上其優化過程的極限可以收斂到全局最優解，優化效果比梯度算法好，但基于梯度的優化算法通常能迅速收斂到局部最優解。若氣動優化問題是一個高度非線性、多峰值的問題，基于梯度的優化算法在求解這類問題時對初始設計點的選擇有很強的依賴性，難以獲得全局最優解。如果氣動優化問題是一個單峰值問題，則梯度算法能夠得到相對滿意的最優解，并且有更高的優化效率。因此，對氣動外形優化問題設計空間的多極值特性有所了解是十分必要的，有助于我們在翼型設計階段合理選擇優化算法，提高飛行器氣動外形優化效率，縮短飛行器設計周期。

在一些氣動外形優化問題中，多極值特性是存在的，例如H.P.Buckley等[3]研究表明在翼型多點優化時出現兩個局部最優點；J.E.Hicken等[4]優化結果表明機翼展向垂直方向外形優化至少有兩個局部最優點，即上翹小翼和下翹小翼；而Lü Zhoujie等[5]探索了CRM(Common Research Model)機翼優化問題的多極值特性，其研究結果證明CRM機翼單點優化問題是單極值問題；O.Chernukhin等[1]探究了翼型和矩形機翼在定升力減阻優化問題中的多極值特性，結果表明，翼型優化在其給定的設計狀態和約束下是一個單極值問題(優化問題沒有考慮力矩約束的影響)，而機翼外形優化存在多極值特性；N.P.Bons等[6]研究了ADODG(Aerodynamic Design Optimization Discussion Group) case 6的多極值特性，發現機翼氣動外形優化結果取決于誘導阻力和粘性阻力之間的權衡。但上述研究對翼型在跨聲速氣動優化過程中的多極值特性探究較少，對優化中變量設計空間的認識并不明晰，缺少梯度算法和全局算法的優化效率和結果的詳細對比。

本文主要使用自由變形方法(Free-Form Deform，簡稱FFD)、基于伴隨的梯度優化算法、全局類優化算法探索ADODG case 2中 RAE2822翼型單點優化問題的多極值特性，并分析如果在設計空間中單極值存在的條件下，全局算法是否還有其優勢。

1 數值求解方法

本文使用RANS方程求解器進行定常求解。三維非定常雷諾平均N-S方程的控制方程為

(1)

式中：Q為守恒向量；Fc為對流矢通量；Fv為粘性通量。

空間離散格式為中心格式，優化時采用S-A(Spalart-Allmaras)模型作為湍流模型。

2 參數化方法及網格變形方法

采用FFD參數化方法[7]，在需要被參數化的幾何外形周圍建立FFD控制體，通過對FFD控制體頂點的移動，實現對FFD控制體中的目標幾何體上的變形。FFD參數化方法能以較少的設計變量光滑地描述曲線、曲面、三維幾何體的幾何外形，并能方便地應用于局部外形修型設計。對于本文的翼型優化，上下翼面各布置10個控制點來改變翼型型面，如圖1所示。

圖1 FFD控制框Fig.1 FFD control box

優化網格采用C網格，如圖2所示。

圖2 翼型優化網格Fig.2 Airfoil optimization mesh

采用基于逆距離權重插值(Inverse Distance Weighting Interpolation，簡稱IDW)方法[8]來實現空間網格更新，該方法魯棒性較高。翼型FFD框擾動及網格變形能力如圖3(a)所示，可以看出：擾動后的網格仍有較好的正交性；流場計算獲得收斂解，壓力分布如圖3(b)所示。

(a) FFD框擾動

(b) 變形后翼型壓力云圖圖3 網格變形Fig.3 Mesh deformation

3 離散伴隨方程法

基于梯度的優化算法主要優點是收斂快速，但是其最大的困難是需要高效準確的梯度計算。本文采用基于伴隨的梯度優化方法[2]，伴隨方程法通過求解原始方程的伴隨方程，可以一次求解出目標函數針對所有設計變量的導數。其計算量基本與設計變量個數無關，計算效率極高，并且能夠保持很好的精度。

離散伴隨方程法直接從已經離散的目標函數及流場控制方程出發，推導出離散形式的伴隨方程。如果伴隨方程的形式和原始流場控制方程的離散形式完全對應，則可以求解出對應于目標函數和控制方程離散形式的數值精確導數。

在氣動外形優化設計中，目標函數一般為阻力系數，構建目標函數表達式為：

FA=f[G(x),Q(x)]

(2)

式中：x為外形設計變量；G(x)為設計變量x確定的CFD計算網格；Q為流場解向量。

流場計算收斂表達式為

RA[Q(x),G(x)]=0

(3)

式中：RA為流場殘差。

構建伴隨方程，進行伴隨算子的求解。

(4)

求解目標函數對設計變量的導數：

(5)

4 優化算法

4.1 梯度優化算法

梯度優化的頂層控制算法采用SNOPT(Sparse Nonlinear Optimizer)算法[9]，能夠高效快速處理大規模非線性約束，同時對求解器的調用次數較少。SNOPT的基本思想是：將有約束的非線性優化問題在每一個迭代步上轉化為二次規劃子問題，在每一步上的二次規劃子問題中求解線性化約束條件下的修正拉格朗日函數二次模型的最優化問題。

4.2 全局優化算法

全局優化算法采用并行增廣拉格朗日乘子粒子群優化算法[10](Augmented Lagrange Multiplier Particle Swarm Optimization，簡稱ALPSO)。粒子群優化算法[11-12](Particle Swarm Optimization，簡稱PSO)是一種模擬鳥群和魚群社會行為的、基于群體智能的進化算法。PSO最初是由 Eberhart 和 Kennedy 設計研發的。在PSO中，優化問題的求解方案是通過追尋當前的最優值來探索設計空間。粒子通過跟隨當前稱為導向的最佳粒子在問題空間中探索。

ALPSO方法利用了PSO方法的優點，可解決不光滑目標函數優化問題，更有可能找到全局最優值。ALPSO方法使用增廣拉格朗日乘子處理約束，可以被用于解決非線性、不可微、非凸問題。

5 基于梯度優化算法的優化設計框架

本文采用的優化設計框架包括幾何參數化模塊、動網格模塊、流場求解模塊、梯度信息求解模塊及優化算法模塊，如圖4所示。首先進行流場求解，優化算法根據目標函數信息和梯度信息產生一組設計變量，更新翼型表面網格及空間網格；然后進行目標函數及目標函數對設計變量梯度的求解；最后將目標函數信息及梯度信息傳遞給優化算法，優化算法結合優化問題的約束要求，產生下一組設計變量，如此循環直至收斂。

圖4 梯度優化算法流程圖Fig.4 Gradient-based optimization algorithm flow chart

6 RAE2822翼型的減阻優化

ADODG定義的算例2：基于RANS方程的RAE2822翼型的跨聲速減阻優化。關于ADODG case2的研究，國內外已做過很多工作[13]。

自由來流馬赫數為0.734，雷諾數為650萬，升力系數為0.824(風洞試驗采集到的攻角大約為2.79°)，對REA2822翼型進行減阻，優化問題受翼型面積和俯仰力矩的約束。優化問題如下：

min:CD

subject to:CL=0.824

CMy≥-0.092

Area≥Areainitial

其中，CD表示阻力系數；CL表示升力系數；CMy表示俯仰力矩系數；Area表示翼型面積；Areainitial表示翼型初始面積。

優化過程中，為便于優化和直觀理解，將面積約束轉化為厚度約束，使優化后翼型的厚度不少于翼型的初始厚度，即t≥tinitial，翼型截面共布置了25個厚度約束。翼型優化也布置了前后緣約束，使前后緣控制點變化量大小相等、方向相反。前后緣約束的布置是為了避免因翼型剪切變形導致翼型攻角改變，造成最終的計算攻角不準確。所有幾何約束示意圖如圖5所示。

圖5 翼型優化幾何約束示意圖Fig.5 Airfoil optimization geometric constraints

翼型采用FFD參數化方法，上下翼面各10個控制點。設計變量包括FFD控制點及翼型攻角，共計21個。

6.1 初始翼型網格收斂性研究

初始翼型的網格收斂性研究結果如表1所示，阻力系數隨網格數的-2/3次冪變化的曲線如圖6所示。

隨著網格量的增加，阻力系數值逐漸收斂至1 count(1 count=0.000 1)以內，表明所使用的CFD方法及網格具有較好的收斂性。優化采用54 848的網格量。

表1 RAE2822翼型不同網格量計算結果Table 1 RAE2822 airfoil calculation results under different gridsize

圖6 RAE2822翼型網格收斂性研究Fig.6 Grid convergence study for the RAE2822 airfoil

6.2 基于梯度算法的單點優化多極值特性研究

采用拉丁超立方抽樣方法對翼型加入初始擾動，翼型控制點的變化范圍為(-0.05，0.05)。

生成七組初始擾動變量，初始翼型如圖7所示。壓力分布對比如圖8所示。

圖7 加入不同隨機初始擾動后的翼型Fig.7 Airfoil after different random initial disturbances

圖8 不同隨機初始擾動后的翼型壓力分布對比Fig.8 Comparison of airfoil pressure distribution after different random initial disturbances

優化后的翼型對比及優化后壓力分布對比分別如圖9～圖10所示。對7個有隨機初始擾動的翼型進行相同的優化，從圖9～圖10可以看出：它們之間只存在細微的差別，阻力系數大小幾乎沒有差別。

圖9 不同初始擾動下優化后翼型對比Fig.9 Comparison of optimized airfoils under different initial disturbance

圖10 不同初始擾動下優化后壓力分布對比Fig.10 Comparison of pressure distribution after optimization under different initial disturbance

梯度優化算法阻力收斂歷史和最優度(Optimality)[9,14]收斂歷史如圖11所示，可以看出：七次隨機初始外形優化迭代次數不同，但是結果均收斂，且目標函數值一樣。

SNOPT優化算法中最優度(Optimality)的含義是考慮了約束的 KKT 條件數值，用以判斷優化是否收斂。不同初始擾動下梯度算法優化結果如表2所示。

(a) 阻力收斂歷史

(b) 最優度收斂歷史圖11 不同初始擾動下梯度優化算法阻力和最優度收斂歷史Fig.11 Gradient-based optimization algorithm convergence history and optimality history under different initial disturbance

翼 型CLCD/countsCMyα/(°)RAE28220.824224.03-0.095463.0100RandomⅠ0.824138.20-0.092003.0704RandomⅡ0.824138.20-0.092003.0696RandomⅢ0.824138.20-0.092003.0700RandomⅣ0.824138.20-0.092003.0704RandomⅤ0.824138.20-0.092003.0701RandomⅥ0.824138.20-0.092003.0701RandomⅦ0.824138.20-0.092003.0699

為了更好地使設計空間可視化，任選三個優化結果，計算該三個優化設計點的設計變量的歐氏距離，如圖12所示。

圖12 多個局部極小值之間的歐氏距離Fig.12 Euclidean distance between multiple local minimums

當前設計空間的最大歐氏距離為0.45c，這三個優化設計點的設計變量之間的距離很短，最小歐氏距離為0.001 4 %c，最大為0.003 0 %c，而優化出的設計點的變量間的歐氏距離遠遠小于0.45c。基于計算數據，可以得出此氣動外形優化問題的設計空間是凸狀的。

為了驗證此氣動外形優化問題的設計空間可能是單峰值的，又選取以上三個優化點的設計變量的平均值，如圖12中圓點所示，其氣動力計算結果如表3所示。

表3 平均變量翼型的氣動力Table 3 Aerodynamics of mean variables airfoil

從表3可以看出：阻力已大幅減小，其數值與Random Ⅰ，Random Ⅱ，Random Ⅲ優化后的阻力數值一樣，并沒有發生突變，驗證了此優化空間很可能是單峰值的。

這次優化也表明氣動優化方法的魯棒性以及梯度伴隨方法的優越性。即使給翼型任意的初始擾動，也能優化出效果較好的翼型，設計點阻力大幅減小。

6.3 優化翼型網格收斂性研究

對Random Ⅰ的優化翼型進行網格收斂性研究，優化翼型的網格收斂性研究結果如表4所示，阻力系數隨網格數的-2/3次冪變化的曲線如圖13所示，可以看出：隨著網格量的增加，阻力系數值逐漸收斂至1 count以內。

表4 優化翼型不同網格量計算結果Table 4 Optimization airfoil calculation results under different gridsize

圖13 優化翼型網格收斂性研究Fig.13 Grid convergence study for the optimization airfoil

6.4 設計變量個數對優化設計結果的影響

進一步探究設計變量維數的變化對優化結果的影響，上下翼面各布置10、18和25個設計變量，梯度算法優化結果如表5和圖14所示。

表5 不同設計變量個數梯度算法優化結果Table 5 Gradient-based algorithm optimization results under different number of design variable

圖14 不同設計變量的優化結果Fig.14 Optimization result with different number of design variables

從表5和圖14可以看出：隨著設計變量的增多，減阻效果更好，但是采用10×2個設計變量和采用25×2個設計變量的梯度優化結果阻力系數差別在1 count以內。

優化后的翼型對比及優化后壓力分布對比如圖15～圖16所示，可以看出：優化結果基本一樣。

圖15 不同設計變量個數優化后翼型對比Fig.15 Comparison of optimized airfoils with different number of design variables

圖16 不同設計變量個數優化后壓力分布對比Fig.16 Comparison of pressure distribution after optimization with different number of design variables

不同設計變量個數梯度優化算法阻力收斂歷史和最優度(Optimality)收斂歷史如圖17所示，可以看出：設計變量增多會導致優化迭代次數增加，但是結果均收斂，且目標函數值相差不大。

(a) 阻力收斂歷史

(b) 最優度收斂歷史圖17 不同設計變量個數梯度優化算法阻力和最優度收斂歷史Fig.17 Gradient-based optimization algorithm convergence history and optimality history under different number of design variable

6.5 全局算法和梯度算法優化結果對比

采用全局優化算法中的粒子群算法求解上述優化問題，在使用相同核數的情況下(8核)，全局算法耗時約兩周，而梯度優化算法僅耗時約5小時即可獲得相對滿意的結果。

全局算法和梯度算法優化結果對比分別如表6、圖18～圖19所示。

表6 不同優化算法優化結果Table 6 The results of different optimization algorithms

圖18 不同優化算法優化前后翼型對比Fig.18 Initial and optimized airfoil comparison using different optimization algorithms

圖19 不同優化算法優化前后壓力分布對比Fig.19 Initial and optimized pressure distribution using different optimization algorithms

從表6可以看出：粒子群優化算法阻力減小85.84 counts，梯度優化算法阻力減小了85.83 counts，優化結果基本一樣。

從圖18～圖19可以看出：翼型與壓力分布優化結果也一致。但是基于伴隨的梯度優化方法計算時間短，效率更高，更適用于實際氣動外形優化時有大規模設計變量的問題。

7 結 論

(1) 在滿足設計約束的前提下，分別采用全局、梯度優化設計平臺對初始翼型進行單目標優化，改善翼型的性能。全局優化算法耗時久，優化效果最好；但梯度優化算法優化效率極高，最終的優化結果也有很好的收益。全局優化算法更適用于設計變量離散、目標函數不連續以及有多個局部最優解的設計問題，但是在大規模設計變量問題中，其計算代價是不可接受的。

(2) RAE2822翼型減阻設計的設計空間很可能是單極值的，但該結論還需要進一步證明。

(3) ADODG的case2可能是一個單峰值氣動設計問題，但是對于涉及大規模幾何設計變量的氣動外形優化問題或者非常規布局的優化設計問題，可能是多峰值的，其峰值數量可能與幾何變形的能力有關，仍需進一步探索。

(4) 設計變量增多對優化結果有一定的改善作用，但是設計變量增多會導致優化迭代次數增加，ADODG的case2中10×2個設計變量已能夠得到好的優化結果。