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

(1.西北工業(yè)大學 航空學院，西安 710072) (2.西北工業(yè)大學 無人系統(tǒng)技術(shù)研究院，西安 710072)

0 引 言

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

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

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

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

1 數(shù)值求解方法

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

(1)

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

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

2 參數(shù)化方法及網(wǎng)格變形方法

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

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

優(yōu)化網(wǎng)格采用C網(wǎng)格，如圖2所示。

圖2 翼型優(yōu)化網(wǎng)格Fig.2 Airfoil optimization mesh

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

(a) FFD框擾動

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

3 離散伴隨方程法

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

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

在氣動外形優(yōu)化設計中，目標函數(shù)一般為阻力系數(shù)，構(gòu)建目標函數(shù)表達式為：

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

(2)

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

流場計算收斂表達式為

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

(3)

式中：RA為流場殘差。

構(gòu)建伴隨方程，進行伴隨算子的求解。

(4)

求解目標函數(shù)對設計變量的導數(shù)：

(5)

4 優(yōu)化算法

4.1 梯度優(yōu)化算法

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

4.2 全局優(yōu)化算法

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

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

5 基于梯度優(yōu)化算法的優(yōu)化設計框架

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

圖4 梯度優(yōu)化算法流程圖Fig.4 Gradient-based optimization algorithm flow chart

6 RAE2822翼型的減阻優(yōu)化

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

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

min:CD

subject to:CL=0.824

CMy≥-0.092

Area≥Areainitial

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

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

圖5 翼型優(yōu)化幾何約束示意圖Fig.5 Airfoil optimization geometric constraints

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

6.1 初始翼型網(wǎng)格收斂性研究

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

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

表1 RAE2822翼型不同網(wǎng)格量計算結(jié)果Table 1 RAE2822 airfoil calculation results under different gridsize

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

6.2 基于梯度算法的單點優(yōu)化多極值特性研究

采用拉丁超立方抽樣方法對翼型加入初始擾動，翼型控制點的變化范圍為(-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

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

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

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

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

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

(a) 阻力收斂歷史

(b) 最優(yōu)度收斂歷史圖11 不同初始擾動下梯度優(yōu)化算法阻力和最優(yōu)度收斂歷史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

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

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

當前設計空間的最大歐氏距離為0.45c，這三個優(yōu)化設計點的設計變量之間的距離很短，最小歐氏距離為0.001 4 %c，最大為0.003 0 %c，而優(yōu)化出的設計點的變量間的歐氏距離遠遠小于0.45c?；谟嬎銛?shù)據(jù)，可以得出此氣動外形優(yōu)化問題的設計空間是凸狀的。

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

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

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

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

6.3 優(yōu)化翼型網(wǎng)格收斂性研究

對Random Ⅰ的優(yōu)化翼型進行網(wǎng)格收斂性研究，優(yōu)化翼型的網(wǎng)格收斂性研究結(jié)果如表4所示，阻力系數(shù)隨網(wǎng)格數(shù)的-2/3次冪變化的曲線如圖13所示，可以看出：隨著網(wǎng)格量的增加，阻力系數(shù)值逐漸收斂至1 count以內(nèi)。

表4 優(yōu)化翼型不同網(wǎng)格量計算結(jié)果Table 4 Optimization airfoil calculation results under different gridsize

圖13 優(yōu)化翼型網(wǎng)格收斂性研究Fig.13 Grid convergence study for the optimization airfoil

6.4 設計變量個數(shù)對優(yōu)化設計結(jié)果的影響

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

表5 不同設計變量個數(shù)梯度算法優(yōu)化結(jié)果Table 5 Gradient-based algorithm optimization results under different number of design variable

圖14 不同設計變量的優(yōu)化結(jié)果Fig.14 Optimization result with different number of design variables

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

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

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

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

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

(a) 阻力收斂歷史

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

6.5 全局算法和梯度算法優(yōu)化結(jié)果對比

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

全局算法和梯度算法優(yōu)化結(jié)果對比分別如表6、圖18～圖19所示。

表6 不同優(yōu)化算法優(yōu)化結(jié)果Table 6 The results of different optimization algorithms

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

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

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

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

7 結(jié) 論

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

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

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

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