ARI_OPT氣動優(yōu)化軟件研究進展及應用

魏闖，楊龍，李春鵬，張鐵軍

中國航空工業(yè)空氣動力研究院，高速高雷諾數(shù)氣動力航空科技重點實驗室，沈陽 110034

傳統(tǒng)的依賴經(jīng)驗、以試湊(Cut and Try)為主的人工修型氣動設(shè)計方法，設(shè)計效率低，很難獲得最佳的氣動外形。隨著計算流體力學(Computational Fluid Dynamics，CFD)技術(shù)、計算機輔助幾何造型設(shè)計技術(shù)(Computer-Aided Design，CAD)以及數(shù)值優(yōu)化等技術(shù)的發(fā)展和日趨成熟，基于高可信度計算流體力的氣動分析與優(yōu)化設(shè)計技術(shù)，目前已廣泛應用于航空、航天飛行器的設(shè)計與研制，在提高飛行器氣動與綜合性能方面正發(fā)揮著越來越重要的作用[1-2]。

完整的飛行器氣動外形優(yōu)化流程包括氣動外形參數(shù)化、網(wǎng)格自動生成、氣動性能評估和高效優(yōu)化算法等。根據(jù)優(yōu)化算法的不同現(xiàn)有的優(yōu)化設(shè)計方法主要分為3類[3]：①遺傳算法、粒子群算法等啟發(fā)式優(yōu)化算法；②梯度優(yōu)化算法；③代理優(yōu)化算法。3種方法各有優(yōu)劣，國內(nèi)外研究者均開展了廣泛的研究。如斯坦福大學Jameson等在伴隨方法方面開展了卓有成效的先驅(qū)性研究工作[4]，而后該校Alonso團隊開發(fā)了基于伴隨方法的SU2[5-6]開源氣動優(yōu)化程序，波音公司研制了MDOPT[7]優(yōu)化軟件用于氣動外形優(yōu)化，美國NASA埃姆斯研究中心CFD求解器FUN3D[8-9]、德國宇航院非結(jié)構(gòu)求解器TAU[10-11]以及法宇航的求解器e1sA[12-13]均集成發(fā)展了氣動外形優(yōu)化設(shè)計功能。國內(nèi)中國空氣動力研究與發(fā)展中心[14-16]、西北工業(yè)大學[17-19]、清華大學[20]、南京航空航天大學[21-23]、中國航空工業(yè)空氣動力研究院[24]等也開展了廣泛深入的研究，分別建立了各自的氣動優(yōu)化設(shè)計工具。

本文簡要介紹了中國航空工業(yè)空氣動力研究院氣動優(yōu)化軟件(Aerodynamics Research Institute Optimization Code,ARI_OPT))各個模塊的基本原理和方法，給出了基于數(shù)值優(yōu)化算例的功能驗證，詳細闡述了其在翼型、多段翼型和飛翼布局機翼等典型單目標和多目標氣動外形優(yōu)化問題的應用算例，旨在促進氣動外形優(yōu)化設(shè)計領(lǐng)域的交流。

1 ARI_OPT氣動優(yōu)化軟件

ARI_OPT是中國航空工業(yè)空氣動力研究院針對飛行器氣動外形優(yōu)化的氣動優(yōu)化軟件。它可以求解任意多約束的單目標、多目標優(yōu)化問題，主要應用于基于高精度CFD的氣動優(yōu)化設(shè)計，也可用于氣動/結(jié)構(gòu)、氣動/隱身等多學科優(yōu)化及其他工程優(yōu)化設(shè)計問題，在軍、民用飛行器氣動外形優(yōu)化設(shè)計上得到了較為廣泛的應用。優(yōu)化軟件集成于中國航空工業(yè)空氣動力研究院大型計算集群，可進行大規(guī)模并行優(yōu)化。該軟件包含氣動外形參數(shù)化、網(wǎng)格自動變形、高逼真度CFD數(shù)值模擬、代理模型和高效優(yōu)化算法等模塊。

1.1 氣動外形參數(shù)化方法

ARI_OPT包含多種氣動外形參數(shù)化方法，如Hicks-Henne[25]、類別形狀函數(shù)變換(Class-Shape Transformation，CST)方法[26-27]、自由曲面變形(Free Form Deform，F(xiàn)FD)[28]、直接操作自由曲面變形(Directly manipulated FFD，DFFD)[29]等參數(shù)化方法。同時發(fā)展了基于CATIA二次開發(fā)的氣動外形參數(shù)化方法，該方法能準確地生成飛行器三維外形，同時可以方便地輸出表面面積、內(nèi)部容積、前后梁高度和位置等幾何信息。另外可直接輸出可供上下游設(shè)計階段直接利用的CAD模型，相比其他參數(shù)化方法如直接控制網(wǎng)格點的方法等，無需再根據(jù)優(yōu)化結(jié)果重新建立外形CAD模型，避免了模型精度損失，提高了工程實用性。典型的機翼參數(shù)化流程如圖1所示。

圖1 基于CATIA民機機翼參數(shù)化建模Fig.1 Parametric modeling of wing based on CATIA automation

1.2 網(wǎng)格自動變形技術(shù)

優(yōu)化過程中利用徑向基函數(shù)(Radial Basis Function，RBF)網(wǎng)格變形技術(shù)實現(xiàn)網(wǎng)格自動生成。徑向基函數(shù)的基本形式為[30]

(1)

本文徑向基函數(shù)采用計算效率與網(wǎng)格變形質(zhì)量都較好的Wendland’s C2函數(shù)[30-31]，并通過貪心算法[32]縮減變形矩陣維數(shù)和基于MPI(Message Passing Interface)的并行求解技術(shù)提高變形效率。典型帶流板下偏的三段翼構(gòu)型，8.1萬 網(wǎng)格條件下，30次網(wǎng)格變形平均耗時不超過1 s，網(wǎng)格變形結(jié)果如圖2所示；NASA的CRM翼身組合體構(gòu)型，1 017萬網(wǎng)格條件下，30次網(wǎng)格變形耗時平均約9.2 s，如圖3所示，具有較高的網(wǎng)格變形質(zhì)量。

圖2 多段翼型網(wǎng)格變形結(jié)果Fig.2 Grid deformation of multi-element airfoil

圖3 翼型組合體網(wǎng)格變形結(jié)果Fig.3 Grid deformation of wing-body

1.3 高逼真度CFD數(shù)值模擬技術(shù)

優(yōu)化過程中氣動性能評估依托中國航空工業(yè)空氣動力研究院航空大規(guī)模CFD平臺ARI_CFD(如圖4所示)，主要使用求解雷諾平均Naiver-Stoke(RANS)方程的結(jié)構(gòu)網(wǎng)格模塊ARI_ENSMB和非結(jié)構(gòu)動態(tài)重疊網(wǎng)格模塊ARI_Overset，2個模塊都采用有限體積法，包含多種常用空間差分格式和湍流模型，在長期應用中得到了大量算例的廣泛驗證[33-41]。對于ARI_ENSMB求解模塊，黏性項采用2階中心差分格式，空間離散無黏項采用Roe-FDS(Roe’s Flux Difference Splitting)分裂格式，隱式LU-SGS(Lower-Upper Symmetric-Gauss-Seidel)時間推進，物面給定絕熱壁、無滑移條件，遠場采用無反射壓力遠場邊界條件，湍流模型采用k-ωSST二方程模型。圖5給出了1 000萬計算網(wǎng)格下DLR-F6翼身組合體構(gòu)型ARI_ENSMB計算馬赫數(shù)Ma=0.75、雷諾數(shù)Re=3.0×106時不同迎角α下升力系數(shù)CL、阻力系數(shù)CD與德宇航TAU、NASA FUN3D、商業(yè)軟件FLUENT等求解器計算以及風洞試驗值的對比結(jié)果[42]，可見ARI_ENSMB與試驗值符合較好，具有較高的計算準度。

圖4 ARI_CFD航空大規(guī)模CFD平臺Fig.4 Aeronautical large-scale CFD platform of ARI_CFD

圖5 DLR-F6翼身組合體構(gòu)型ARI_ENSMB計算與試驗值及其他求解器結(jié)果對比(Ma=0.75、Re=3.0×106)Fig.5 Comparison of results between test data and ARI_ENSMB and other solvers for DLR-F6 wingbody(Ma=0.75,Re=3.0×106)

1.4 代理模型技術(shù)

代理模型技術(shù)采用的算法和程序模塊主要來自西北工業(yè)大學韓忠華教授團隊[3,19,43-47]。試驗設(shè)計方法包括拉丁超立方(LHS)、均勻設(shè)計(UD)、蒙特卡洛抽樣(MC)等，代理模型包含二次響應面(PRSM)、Kriging模型、梯度增強Kriging模型(GEK)、分層Kriging模型(HK)、徑向基函數(shù)(RBFs)等多種代理模型，改善期望(maximizing Expected Improvement，EI)、目標函數(shù)值最小(Minimizing Surrogate Prediction，MSP)、目標函數(shù)誤差最大 (maximizing Meansquared Error，MSE)、低置信邊界方法(minimizing Lower-Confidence Bounding，LCB)以及概率提升(Probability of Improvement，PI)等多種優(yōu)化加點準則[19]，既可單獨使用也可根據(jù)優(yōu)化任務需求組合使用多種加點方法。

1.5 高效優(yōu)化算法

優(yōu)化策略包含遺傳算法優(yōu)化和基于代理模型的優(yōu)化方法(簡稱代理優(yōu)化(Surrogate-Based Optimization，SBO))2種，其中代理優(yōu)化方法具有效率高、魯棒性好的優(yōu)點，得到了較為廣泛的應用。典型代理優(yōu)化方法流程如圖6所示，典型優(yōu)化流程如下：① 對設(shè)計空間進行試驗設(shè)計，獲得初始樣本點并調(diào)用數(shù)值求解模塊獲得響應值，構(gòu)建初始代理模型；② 基于代理模型，采用傳統(tǒng)優(yōu)化算法求解相應的子優(yōu)化問題，以很小的計算代價，對最優(yōu)解進行預測，按照一定的優(yōu)化加點準則獲得新樣本點；③ 調(diào)用數(shù)值求解模塊計算得到新樣本點響應值，并將結(jié)果添加到現(xiàn)有數(shù)據(jù)集中，不斷更新代理模型，直到所產(chǎn)生的樣本點序列收斂于局部或全局最優(yōu)解。子優(yōu)化算法包括Hooke-Jeeves模式搜索、擬牛頓梯度優(yōu)化、序列二次規(guī)劃法(SQP)、單/多目標遺傳算法等多種成熟優(yōu)化算法等。對于典型三維優(yōu)化問題，在一般情況下，利用500 CPU核完成一輪優(yōu)化耗時不超過24 h。

圖6 ARI _OPT軟件基于代理模型優(yōu)化流程Fig.6 Flow chart of surrogate-based optimization of ARI _OPT software

2 基于數(shù)值優(yōu)化算例的功能驗證

2.1 六峰值駝背測試函數(shù)單目標優(yōu)化

六峰值駝背測試函數(shù)(Aix-hump Camel Function)具有6個局部極小點,其中有2個為全局極小點，優(yōu)化問題描述為

(2)

分別進行了遺傳算法和代理優(yōu)化方法測試，其中遺傳算法采用多島遺傳算法(Multi-Island GA，MIGA)，每代種群總數(shù)為32，分4個島，優(yōu)化50代，共調(diào)用1 600 次計算。代理優(yōu)化方法初始采樣點20個，試驗設(shè)計方法為LHS，代理模型為Kriging模型，EI+MSP+LCB+PI混合加點方式，每次加16點，共加16輪，共調(diào)用計算276次。表1給出了30次測試的平均結(jié)果，可見代理優(yōu)化方法在效率和優(yōu)化結(jié)果上都優(yōu)于多島遺傳算法。

表1 多島遺傳算法和代理優(yōu)化方法測試結(jié)果Table 1 Comparison of test results of SBO with MIGA

2.2 ZDT3多目標函數(shù)

ZDT3問題為多目標優(yōu)化領(lǐng)域普遍使用的標準測試函數(shù)，真實解集Pareto凸、非連續(xù)。優(yōu)化問題描述為

(3)

分別進行了遺傳算法和代理優(yōu)化方法測試，其中遺傳算法采用快速非支配排序遺傳算法(NSGA-2)，每代種群總數(shù)為160，經(jīng)過70代優(yōu)化達到收斂要求，共調(diào)用11 200次計算。代理優(yōu)化方法初始采樣點20個，試驗設(shè)計方法為LHS，代理模型為Kriging模型，Pareto優(yōu)化時的MSP加點方法，總計調(diào)用計算1 260次。圖7給了2種方法的優(yōu)化結(jié)果,可見2種方法都得到的Pareto前沿與真實Pareto前沿基本重合，都獲得了較優(yōu)的優(yōu)化結(jié)果，但代理優(yōu)化方法調(diào)用計算次數(shù)較少，效率較高。

圖7 NSGA-2和代理優(yōu)化方法測試結(jié)果Fig.7 Test results of NSGA-2 and SBO

3 典型應用案例

3.1 翼型減阻單目標優(yōu)化

針對RAE2822翼型進行跨聲速減阻優(yōu)化設(shè)計，設(shè)計狀態(tài)如下：

Ma=0.73，Re=6.5×106，CL=0.8

優(yōu)化問題描述為

(4)

式中:Cm為俯仰力矩系數(shù)。

參數(shù)化方法采用CST翼型參數(shù)化方法，上下翼面總計18個設(shè)計變量，通過求解定?？蓧篟ANS方程獲得優(yōu)化過程中翼型氣動性能，湍流模型為k-ωSST兩方程模型。優(yōu)化采用代理優(yōu)化方法，初始采樣點36個，試驗設(shè)計方法為LHS，代理模型為Kriging模型，EI+MSP+LCB+PI混合加點方式，每次加點4個，共加43次，總計調(diào)用CFD計算208次。表2給出優(yōu)化翼型和RAE2822氣動力結(jié)果對比，優(yōu)化構(gòu)型減阻約31.74%。

表2 RAE2822和優(yōu)化翼型氣動性能對比Table 2 Comparison of aerodynamic performance of RAE2822 and optimal airfoil

圖8和圖9分別給出了優(yōu)化翼型和RAE2822外形壓力系數(shù)Cp分布對比曲線，c為弦長?？梢妰?yōu)化翼型顯著地減弱了激波強度，使阻力大幅降低。

圖8 RAE2822和優(yōu)化翼型外形對比Fig.8 Comparison of shapes for RAE2822 and optimum airfoils

圖9 RAE2822翼型優(yōu)化前后的表面壓力系數(shù)分布對比Fig.9 Comparison of surface pressure coefficient distributions for RAE2822 airfoil before and after optimum

3.2 寬速域翼型減阻多目標優(yōu)化

針對某寬速域翼型，為提高其跨聲速和超聲速氣動性能進行多目標優(yōu)化，優(yōu)化目的是提高Ma=0.72和Ma=3.15時的最大升阻比K。

優(yōu)化目標：

(5)

約束條件：

(6)

參數(shù)化方法采用改進的Hicks-Henne翼型參數(shù)化方法優(yōu)化，上下翼面總計16個設(shè)計變量，通過求解定?？蓧篟ANS方程獲得優(yōu)化過程中翼型氣動性能，湍流模型為k-ωSST兩方程模型。優(yōu)化算法采用NSGA-2多目標優(yōu)化算法，每代種群有60個個體，總代數(shù)為50，總共完成3 000個翼型構(gòu)型計算評估。

圖10給出優(yōu)化得到的Pareto前沿，共計89個 優(yōu)化構(gòu)型，相比初始構(gòu)型，Pareto前沿優(yōu)化構(gòu)型Ma=0.72、α=3.0°(KDP1)和Ma=3.15、α=6.0°(KDP2)下升阻比都有所增加，其中OPT1和OPT3分別為KDP2和KDP1最優(yōu)構(gòu)型。沿著Pareto前沿由構(gòu)型OPT1到OPT3，優(yōu)化構(gòu)型Ma=0.72、α=3.0°升阻比越來越大，相比初始構(gòu)型最大增加37.17%(OPT3)；反之，Ma=3.15,α=6.0° 升阻比數(shù)越來越大，相比初始構(gòu)型最大增加11.39%(OPT3)，OPT2為KDP1和KDP2都增加較多優(yōu)化構(gòu)型，相比初始構(gòu)型分別增加21.48% 和8.01%。

圖11給出了優(yōu)化翼型外形對比結(jié)果對比，圖12 和圖13給出了優(yōu)化構(gòu)型和初始構(gòu)型Ma=0.72 和Ma=3.15下升力和升阻比對比結(jié)果?？梢奙a=0.72下優(yōu)化構(gòu)型升力性能略有增加，Ma= 3.15下升力性能基本不變，滿足約束，Ma= 0.72和Ma=3.15優(yōu)化構(gòu)型升阻比都獲得了較大提升，達到了優(yōu)化目的。

圖10 寬速域翼型多目標優(yōu)化Pareto前沿Fig.10 Pareto multi-optimal front of wide-speed range airfoil

圖11 優(yōu)化前后翼型外形對比Fig.11 Comparison of shapes for initial and optimum airfoils

圖12 Ma=0.72時初始和優(yōu)化翼型氣動性能對比Fig.12 Comparison of aerodynamic performance of initial and optimal airfoils at Ma=0.72

圖13 Ma=3.15時初始和優(yōu)化翼型氣動性能對比Fig.13 Comparison of aerodynamic performance of initial and optimal airfoils at Ma=3.15

3.3 多段翼型多目標優(yōu)化

針對帶擾流板下偏先進三段翼型著陸構(gòu)型，開展氣動外型和縫道參數(shù)優(yōu)化，以提高其升力性能，計算條件為Ma=0.2、Re=45.6×106，優(yōu)化問題描述如下：

(7)

通過橢圓方程控制生成多段翼型外型，包括6個外型設(shè)計參數(shù)，前緣縫翼和后緣襟翼縫道寬度、搭接量和偏角6個設(shè)計參數(shù)，以及擾流板下偏角，總計13個設(shè)計變量。優(yōu)化采用NSGA-2每代種群有70個個體，總代數(shù)為50，總共完成3 500個 構(gòu)型計算評估。圖14給出優(yōu)化結(jié)束后，優(yōu)化得到的Pareto前沿，相比初始構(gòu)型，Pareto前沿優(yōu)化構(gòu)型迎角8°和27°下升力系數(shù)都有所增加，其中OPT1和OPT3分別為迎角27°和8°升力性能最優(yōu)構(gòu)型。沿著Pareto前沿由構(gòu)型OPT1到OPT3，迎角8°時升力系數(shù)越來越大，反之，迎角27°升力系數(shù)越來越大，OPT2為綜合迎角8°和迎角27°升力性能的優(yōu)化構(gòu)型。

圖15給出了多段翼型典型優(yōu)化構(gòu)型外形對比，表3給出了典型優(yōu)化構(gòu)型設(shè)計點升力系數(shù)對比結(jié)果。相比基本構(gòu)型，優(yōu)化構(gòu)型升力系數(shù)有明顯增加，其中8°迎角升力系數(shù)最大增加 9.24%(OPT3)，27°迎角升力系數(shù)最大增加1.93%(OPT1)；從圖16多段翼型典型優(yōu)化構(gòu)型升力系數(shù)隨迎角變化曲線可以看出，構(gòu)型OPT1有最大的最大升力系數(shù)，線性段增加較小，而構(gòu)型OPT3線性段升力系數(shù)增加較大，最大升力系數(shù)增加較小，OPT2構(gòu)型8°迎角和27°升力系數(shù)都有一定的增加，分別增加7.74%和1.43%。

圖14 多段翼型多目標優(yōu)化Pareto前沿Fig.14 Pareto multi-optimal front of multi-element airfoil

圖15 優(yōu)化前后多段翼型外形對比Fig.15 Comparison of shapes for baseline and optimum multi-element airfoils

表3 基礎(chǔ)和優(yōu)化多段翼型升力系數(shù)對比Table 3 Comparison of lift coefficient of baseline and optimal multi-element airfoils

圖16 基礎(chǔ)和優(yōu)化多段翼型升力系數(shù)曲線Fig.16 Curves of lift coefficient of baseline and optimal multi-element airfoils

3.4 飛翼布局多目標優(yōu)化

針對典型飛翼布局，開展跨聲速減阻多目標優(yōu)化設(shè)計，優(yōu)化目的是提高跨聲速設(shè)計點升阻比和阻力發(fā)散馬赫數(shù)，優(yōu)化問題描述如下：

(8)

式中:y為機翼展向站位；B為機翼展長。

參數(shù)化方法采用翼型剖面CST+CATIA二次開發(fā)三維成型方法，共計5個優(yōu)化剖面(圖17紅色曲線),每個剖面18個設(shè)計變量，總計90個設(shè)計變量。采用代理優(yōu)化方法，初始采樣點270個，試驗設(shè)計方法為LHS，代理模型為Kriging模型，采用Pareto優(yōu)化時的EI加點方法，總計調(diào)用CFD計算1 200次。

圖18給出了優(yōu)化Pareto前沿，表4給出了典型優(yōu)化構(gòu)型與基礎(chǔ)構(gòu)型阻力系數(shù)對比結(jié)果，相比基礎(chǔ)構(gòu)型，Ma=0.72、CL=0.35時優(yōu)化構(gòu)型減阻13.17%，Ma=0.78、CL=0.30時優(yōu)化構(gòu)型減阻19.56%。從圖19和圖20不同馬赫數(shù)下優(yōu)化構(gòu)型和基礎(chǔ)構(gòu)型上翼面壓力云圖對比結(jié)果可以看出，優(yōu)化構(gòu)型顯著地減小了上翼面激波強度，從而使阻力大幅降低。從圖21不同升力系數(shù)下阻力系數(shù)隨馬赫數(shù)變化結(jié)果可以看出，優(yōu)化構(gòu)型阻力發(fā)散馬赫數(shù)獲得了一定的增加，達到了優(yōu)化目標。

圖17 基于CATIA飛翼布局參數(shù)化建模Fig.17 Parametric modeling of flying wing based on CATIA automation

圖18 飛翼布局多目標優(yōu)化Pareto前沿Fig.18 Pareto multi-optimal front of flying wing

表4 基礎(chǔ)和優(yōu)化構(gòu)型氣動性能對比Table 4 Comparison of aerodynamic performance of baseline and optimal configuration

另外，此算例中樣本點總數(shù)達到了1 200個，進行代理模型訓練時耗時約2 h，超過了單個樣本CFD數(shù)值求解的時間，制約著優(yōu)化效率的提高，如何提高大樣本點下代理模型訓練效率是ARI_OPT軟件后續(xù)的一個改進方向。

圖19 CL=0.35基礎(chǔ)和優(yōu)化構(gòu)型壓力系數(shù)云圖對比Fig.19 Comparison of pressure coefficient contour of baseline and optimal configuration atCL=0.35

圖20 CL=0.30基礎(chǔ)和優(yōu)化構(gòu)型壓力系數(shù)云圖對比Fig.20 Comparison of pressure coefficient contour of baseline and optimal configuration atCL=0.30

圖21 基礎(chǔ)和優(yōu)化多段翼構(gòu)型阻力隨馬赫數(shù)變化曲線Fig.21 Curvers of CD-Ma of baseline and optimal multi-element airfoils

4 結(jié) 論

本文簡要介紹了ARI_OPT氣動優(yōu)化軟件各個模塊的基本原理和方法及單個模塊的驗證結(jié)果，給出了ARI_OPT針對數(shù)值函數(shù)算例的功能驗證結(jié)果和寬速域翼型、多段翼型和飛翼布局機翼等典型單目標和多目標氣動外形優(yōu)化問題的應用算例，表明了其可靠性和適用性。得出以下結(jié)論：

1) 數(shù)值函數(shù)算例驗證結(jié)果表明，相比基于遺傳算法優(yōu)化，代理優(yōu)化方法計算較少的樣本點就能獲得較優(yōu)優(yōu)化結(jié)果，效率較高。

2) 經(jīng)過數(shù)值函數(shù)驗證和寬速域翼型、多段翼型、飛翼布局機翼等典型氣動外形優(yōu)化問題應用檢驗，ARI_OPT氣動優(yōu)化設(shè)計工具是有效的且具有較高的優(yōu)化效率。

3) 大樣本點(1 000以上)情況下進行代理模型訓練時耗時較長，制約著優(yōu)化效率的提高，如何提高大樣本點下代理模型訓練效率是ARI_OPT軟件后續(xù)的一個改進方向。

4) 本文ARI_OPT軟件的多目標優(yōu)化算例，僅僅是氣動性能的多目標優(yōu)化，理論完全可以推廣至氣動/結(jié)構(gòu)、氣動/隱身等多學科優(yōu)化問題，如何進一步應用于多學科優(yōu)化問題將是下一步工作的主要內(nèi)容。