魏闖,楊龍,李春鵬,張鐵軍
中國航空工業(yè)空氣動力研究院,高速高雷諾數(shù)氣動力航空科技重點實驗室,沈陽 110034
傳統(tǒng)的依賴經(jīng)驗、以試湊(Cut and Try)為主的人工修型氣動設(shè)計方法,設(shè)計效率低,很難獲得最佳的氣動外形。隨著計算流體力學(xué)(Computational Fluid Dynamics,CFD)技術(shù)、計算機輔助幾何造型設(shè)計技術(shù)(Computer-Aided Design,CAD)以及數(shù)值優(yōu)化等技術(shù)的發(fā)展和日趨成熟,基于高可信度計算流體力的氣動分析與優(yōu)化設(shè)計技術(shù),目前已廣泛應(yīng)用于航空、航天飛行器的設(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)外研究者均開展了廣泛的研究。如斯坦福大學(xué)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è)大學(xué)[17-19]、清華大學(xué)[20]、南京航空航天大學(xué)[21-23]、中國航空工業(yè)空氣動力研究院[24]等也開展了廣泛深入的研究,分別建立了各自的氣動優(yōu)化設(shè)計工具。
本文簡要介紹了中國航空工業(yè)空氣動力研究院氣動優(yōu)化軟件(Aerodynamics Research Institute Optimization Code,ARI_OPT))各個模塊的基本原理和方法,給出了基于數(shù)值優(yōu)化算例的功能驗證,詳細闡述了其在翼型、多段翼型和飛翼布局機翼等典型單目標(biāo)和多目標(biāo)氣動外形優(yōu)化問題的應(yīng)用算例,旨在促進氣動外形優(yōu)化設(shè)計領(lǐng)域的交流。
ARI_OPT是中國航空工業(yè)空氣動力研究院針對飛行器氣動外形優(yōu)化的氣動優(yōu)化軟件。它可以求解任意多約束的單目標(biāo)、多目標(biāo)優(yōu)化問題,主要應(yīng)用于基于高精度CFD的氣動優(yōu)化設(shè)計,也可用于氣動/結(jié)構(gòu)、氣動/隱身等多學(xué)科優(yōu)化及其他工程優(yōu)化設(shè)計問題,在軍、民用飛行器氣動外形優(yōu)化設(shè)計上得到了較為廣泛的應(yīng)用。優(yōu)化軟件集成于中國航空工業(yè)空氣動力研究院大型計算集群,可進行大規(guī)模并行優(yōu)化。該軟件包含氣動外形參數(shù)化、網(wǎng)格自動變形、高逼真度CFD數(shù)值模擬、代理模型和高效優(yōu)化算法等模塊。
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ù)化方法,該方法能準(zhǔn)確地生成飛行器三維外形,同時可以方便地輸出表面面積、內(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
優(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
優(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個模塊都采用有限體積法,包含多種常用空間差分格式和湍流模型,在長期應(yīng)用中得到了大量算例的廣泛驗證[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等求解器計算以及風(fēng)洞試驗值的對比結(jié)果[42],可見ARI_ENSMB與試驗值符合較好,具有較高的計算準(zhǔn)度。
圖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)
代理模型技術(shù)采用的算法和程序模塊主要來自西北工業(yè)大學(xué)韓忠華教授團隊[3,19,43-47]。試驗設(shè)計方法包括拉丁超立方(LHS)、均勻設(shè)計(UD)、蒙特卡洛抽樣(MC)等,代理模型包含二次響應(yīng)面(PRSM)、Kriging模型、梯度增強Kriging模型(GEK)、分層Kriging模型(HK)、徑向基函數(shù)(RBFs)等多種代理模型,改善期望(maximizing Expected Improvement,EI)、目標(biāo)函數(shù)值最小(Minimizing Surrogate Prediction,MSP)、目標(biāo)函數(shù)誤差最大 (maximizing Meansquared Error,MSE)、低置信邊界方法(minimizing Lower-Confidence Bounding,LCB)以及概率提升(Probability of Improvement,PI)等多種優(yōu)化加點準(zhǔn)則[19],既可單獨使用也可根據(jù)優(yōu)化任務(wù)需求組合使用多種加點方法。
優(yōu)化策略包含遺傳算法優(yōu)化和基于代理模型的優(yōu)化方法(簡稱代理優(yōu)化(Surrogate-Based Optimization,SBO))2種,其中代理優(yōu)化方法具有效率高、魯棒性好的優(yōu)點,得到了較為廣泛的應(yīng)用。典型代理優(yōu)化方法流程如圖6所示,典型優(yōu)化流程如下:① 對設(shè)計空間進行試驗設(shè)計,獲得初始樣本點并調(diào)用數(shù)值求解模塊獲得響應(yīng)值,構(gòu)建初始代理模型;② 基于代理模型,采用傳統(tǒng)優(yōu)化算法求解相應(yīng)的子優(yōu)化問題,以很小的計算代價,對最優(yōu)解進行預(yù)測,按照一定的優(yōu)化加點準(zhǔn)則獲得新樣本點;③ 調(diào)用數(shù)值求解模塊計算得到新樣本點響應(yīng)值,并將結(jié)果添加到現(xiàn)有數(shù)據(jù)集中,不斷更新代理模型,直到所產(chǎn)生的樣本點序列收斂于局部或全局最優(yōu)解。子優(yōu)化算法包括Hooke-Jeeves模式搜索、擬牛頓梯度優(yōu)化、序列二次規(guī)劃法(SQP)、單/多目標(biāo)遺傳算法等多種成熟優(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
六峰值駝背測試函數(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
ZDT3問題為多目標(biāo)優(yōu)化領(lǐng)域普遍使用的標(biāo)準(zhǔn)測試函數(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
針對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
針對某寬速域翼型,為提高其跨聲速和超聲速氣動性能進行多目標(biāo)優(yōu)化,優(yōu)化目的是提高Ma=0.72和Ma=3.15時的最大升阻比K。
優(yōu)化目標(biāo):
(5)
約束條件:
(6)
參數(shù)化方法采用改進的Hicks-Henne翼型參數(shù)化方法優(yōu)化,上下翼面總計16個設(shè)計變量,通過求解定??蓧篟ANS方程獲得優(yōu)化過程中翼型氣動性能,湍流模型為k-ωSST兩方程模型。優(yōu)化算法采用NSGA-2多目標(biāo)優(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 寬速域翼型多目標(biāo)優(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
針對帶擾流板下偏先進三段翼型著陸構(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 多段翼型多目標(biāo)優(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
針對典型飛翼布局,開展跨聲速減阻多目標(biāo)優(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)化目標(biāo)。
圖17 基于CATIA飛翼布局參數(shù)化建模Fig.17 Parametric modeling of flying wing based on CATIA automation
圖18 飛翼布局多目標(biāo)優(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個,進行代理模型訓(xùn)練時耗時約2 h,超過了單個樣本CFD數(shù)值求解的時間,制約著優(yōu)化效率的提高,如何提高大樣本點下代理模型訓(xùn)練效率是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
本文簡要介紹了ARI_OPT氣動優(yōu)化軟件各個模塊的基本原理和方法及單個模塊的驗證結(jié)果,給出了ARI_OPT針對數(shù)值函數(shù)算例的功能驗證結(jié)果和寬速域翼型、多段翼型和飛翼布局機翼等典型單目標(biāo)和多目標(biāo)氣動外形優(yōu)化問題的應(yīng)用算例,表明了其可靠性和適用性。得出以下結(jié)論:
1) 數(shù)值函數(shù)算例驗證結(jié)果表明,相比基于遺傳算法優(yōu)化,代理優(yōu)化方法計算較少的樣本點就能獲得較優(yōu)優(yōu)化結(jié)果,效率較高。
2) 經(jīng)過數(shù)值函數(shù)驗證和寬速域翼型、多段翼型、飛翼布局機翼等典型氣動外形優(yōu)化問題應(yīng)用檢驗,ARI_OPT氣動優(yōu)化設(shè)計工具是有效的且具有較高的優(yōu)化效率。
3) 大樣本點(1 000以上)情況下進行代理模型訓(xùn)練時耗時較長,制約著優(yōu)化效率的提高,如何提高大樣本點下代理模型訓(xùn)練效率是ARI_OPT軟件后續(xù)的一個改進方向。
4) 本文ARI_OPT軟件的多目標(biāo)優(yōu)化算例,僅僅是氣動性能的多目標(biāo)優(yōu)化,理論完全可以推廣至氣動/結(jié)構(gòu)、氣動/隱身等多學(xué)科優(yōu)化問題,如何進一步應(yīng)用于多學(xué)科優(yōu)化問題將是下一步工作的主要內(nèi)容。