ARI_OPT气动优化软件研究进展及应用

2020-06-03魏闯杨龙李春鹏张铁军

航空学报 2020年5期

魏闯，杨龙，李春鹏，张铁军

中国航空工业空气动力研究院，高速高雷诺数气动力航空科技重点实验室，沈阳 110034

传统的依赖经验、以试凑(Cut and Try)为主的人工修型气动设计方法，设计效率低，很难获得最佳的气动外形。随着计算流体力学(Computational Fluid Dynamics，CFD)技术、计算机辅助几何造型设计技术(Computer-Aided Design，CAD)以及数值优化等技术的发展和日趋成熟，基于高可信度计算流体力的气动分析与优化设计技术，目前已广泛应用于航空、航天飞行器的设计与研制，在提高飞行器气动与综合性能方面正发挥着越来越重要的作用[1-2]。

完整的飞行器气动外形优化流程包括气动外形参数化、网格自动生成、气动性能评估和高效优化算法等。根据优化算法的不同现有的优化设计方法主要分为3类[3]：①遗传算法、粒子群算法等启发式优化算法；②梯度优化算法；③代理优化算法。3种方法各有优劣，国内外研究者均开展了广泛的研究。如斯坦福大学Jameson等在伴随方法方面开展了卓有成效的先驱性研究工作[4]，而后该校Alonso团队开发了基于伴随方法的SU2[5-6]开源气动优化程序，波音公司研制了MDOPT[7]优化软件用于气动外形优化，美国NASA埃姆斯研究中心CFD求解器FUN3D[8-9]、德国宇航院非结构求解器TAU[10-11]以及法宇航的求解器e1sA[12-13]均集成发展了气动外形优化设计功能。国内中国空气动力研究与发展中心[14-16]、西北工业大学[17-19]、清华大学[20]、南京航空航天大学[21-23]、中国航空工业空气动力研究院[24]等也开展了广泛深入的研究，分别建立了各自的气动优化设计工具。

本文简要介绍了中国航空工业空气动力研究院气动优化软件(Aerodynamics Research Institute Optimization Code,ARI_OPT))各个模块的基本原理和方法，给出了基于数值优化算例的功能验证，详细阐述了其在翼型、多段翼型和飞翼布局机翼等典型单目标和多目标气动外形优化问题的应用算例，旨在促进气动外形优化设计领域的交流。

1 ARI_OPT气动优化软件

ARI_OPT是中国航空工业空气动力研究院针对飞行器气动外形优化的气动优化软件。它可以求解任意多约束的单目标、多目标优化问题，主要应用于基于高精度CFD的气动优化设计，也可用于气动/结构、气动/隐身等多学科优化及其他工程优化设计问题，在军、民用飞行器气动外形优化设计上得到了较为广泛的应用。优化软件集成于中国航空工业空气动力研究院大型计算集群，可进行大规模并行优化。该软件包含气动外形参数化、网格自动变形、高逼真度CFD数值模拟、代理模型和高效优化算法等模块。

1.1 气动外形参数化方法

ARI_OPT包含多种气动外形参数化方法，如Hicks-Henne[25]、类别形状函数变换(Class-Shape Transformation，CST)方法[26-27]、自由曲面变形(Free Form Deform，FFD)[28]、直接操作自由曲面变形(Directly manipulated FFD，DFFD)[29]等参数化方法。同时发展了基于CATIA二次开发的气动外形参数化方法，该方法能准确地生成飞行器三维外形，同时可以方便地输出表面面积、内部容积、前后梁高度和位置等几何信息。另外可直接输出可供上下游设计阶段直接利用的CAD模型，相比其他参数化方法如直接控制网格点的方法等，无需再根据优化结果重新建立外形CAD模型，避免了模型精度损失，提高了工程实用性。典型的机翼参数化流程如图1所示。

图1 基于CATIA民机机翼参数化建模Fig.1 Parametric modeling of wing based on CATIA automation

1.2 网格自动变形技术

优化过程中利用径向基函数(Radial Basis Function，RBF)网格变形技术实现网格自动生成。径向基函数的基本形式为[30]

(1)

本文径向基函数采用计算效率与网格变形质量都较好的Wendland’s C2函数[30-31]，并通过贪心算法[32]缩减变形矩阵维数和基于MPI(Message Passing Interface)的并行求解技术提高变形效率。典型带流板下偏的三段翼构型，8.1万网格条件下，30次网格变形平均耗时不超过1 s，网格变形结果如图2所示；NASA的CRM翼身组合体构型，1 017万网格条件下，30次网格变形耗时平均约9.2 s，如图3所示，具有较高的网格变形质量。

图2 多段翼型网格变形结果Fig.2 Grid deformation of multi-element airfoil

图3 翼型组合体网格变形结果Fig.3 Grid deformation of wing-body

1.3 高逼真度CFD数值模拟技术

优化过程中气动性能评估依托中国航空工业空气动力研究院航空大规模CFD平台ARI_CFD(如图4所示)，主要使用求解雷诺平均Naiver-Stoke(RANS)方程的结构网格模块ARI_ENSMB和非结构动态重叠网格模块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万计算网格下DLR-F6翼身组合体构型ARI_ENSMB计算马赫数Ma=0.75、雷诺数Re=3.0×106时不同迎角α下升力系数CL、阻力系数CD与德宇航TAU、NASA FUN3D、商业软件FLUENT等求解器计算以及风洞试验值的对比结果[42]，可见ARI_ENSMB与试验值符合较好，具有较高的计算准度。

图4 ARI_CFD航空大规模CFD平台Fig.4 Aeronautical large-scale CFD platform of ARI_CFD

图5 DLR-F6翼身组合体构型ARI_ENSMB计算与试验值及其他求解器结果对比(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 代理模型技术

代理模型技术采用的算法和程序模块主要来自西北工业大学韩忠华教授团队[3,19,43-47]。试验设计方法包括拉丁超立方(LHS)、均匀设计(UD)、蒙特卡洛抽样(MC)等，代理模型包含二次响应面(PRSM)、Kriging模型、梯度增强Kriging模型(GEK)、分层Kriging模型(HK)、径向基函数(RBFs)等多种代理模型，改善期望(maximizing Expected Improvement，EI)、目标函数值最小(Minimizing Surrogate Prediction，MSP)、目标函数误差最大 (maximizing Meansquared Error，MSE)、低置信边界方法(minimizing Lower-Confidence Bounding，LCB)以及概率提升(Probability of Improvement，PI)等多种优化加点准则[19]，既可单独使用也可根据优化任务需求组合使用多种加点方法。

1.5 高效优化算法

优化策略包含遗传算法优化和基于代理模型的优化方法(简称代理优化(Surrogate-Based Optimization，SBO))2种，其中代理优化方法具有效率高、鲁棒性好的优点，得到了较为广泛的应用。典型代理优化方法流程如图6所示，典型优化流程如下：① 对设计空间进行试验设计，获得初始样本点并调用数值求解模块获得响应值，构建初始代理模型；② 基于代理模型，采用传统优化算法求解相应的子优化问题，以很小的计算代价，对最优解进行预测，按照一定的优化加点准则获得新样本点；③ 调用数值求解模块计算得到新样本点响应值，并将结果添加到现有数据集中，不断更新代理模型，直到所产生的样本点序列收敛于局部或全局最优解。子优化算法包括Hooke-Jeeves模式搜索、拟牛顿梯度优化、序列二次规划法(SQP)、单/多目标遗传算法等多种成熟优化算法等。对于典型三维优化问题，在一般情况下，利用500 CPU核完成一轮优化耗时不超过24 h。

图6 ARI _OPT软件基于代理模型优化流程Fig.6 Flow chart of surrogate-based optimization of ARI _OPT software

2 基于数值优化算例的功能验证

2.1 六峰值驼背测试函数单目标优化

六峰值驼背测试函数(Aix-hump Camel Function)具有6个局部极小点,其中有2个为全局极小点，优化问题描述为

(2)

分别进行了遗传算法和代理优化方法测试，其中遗传算法采用多岛遗传算法(Multi-Island GA，MIGA)，每代种群总数为32，分4个岛，优化50代，共调用1 600 次计算。代理优化方法初始采样点20个，试验设计方法为LHS，代理模型为Kriging模型，EI+MSP+LCB+PI混合加点方式，每次加16点，共加16轮，共调用计算276次。表1给出了30次测试的平均结果，可见代理优化方法在效率和优化结果上都优于多岛遗传算法。

表1 多岛遗传算法和代理优化方法测试结果Table 1 Comparison of test results of SBO with MIGA

2.2 ZDT3多目标函数

ZDT3问题为多目标优化领域普遍使用的标准测试函数，真实解集Pareto凸、非连续。优化问题描述为

(3)

分别进行了遗传算法和代理优化方法测试，其中遗传算法采用快速非支配排序遗传算法(NSGA-2)，每代种群总数为160，经过70代优化达到收敛要求，共调用11 200次计算。代理优化方法初始采样点20个，试验设计方法为LHS，代理模型为Kriging模型，Pareto优化时的MSP加点方法，总计调用计算1 260次。图7给了2种方法的优化结果,可见2种方法都得到的Pareto前沿与真实Pareto前沿基本重合，都获得了较优的优化结果，但代理优化方法调用计算次数较少，效率较高。

图7 NSGA-2和代理优化方法测试结果Fig.7 Test results of NSGA-2 and SBO

3 典型应用案例

3.1 翼型减阻单目标优化

针对RAE2822翼型进行跨声速减阻优化设计，设计状态如下：

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

优化问题描述为

(4)

式中:Cm为俯仰力矩系数。

参数化方法采用CST翼型参数化方法，上下翼面总计18个设计变量，通过求解定常可压RANS方程获得优化过程中翼型气动性能，湍流模型为k-ωSST两方程模型。优化采用代理优化方法，初始采样点36个，试验设计方法为LHS，代理模型为Kriging模型，EI+MSP+LCB+PI混合加点方式，每次加点4个，共加43次，总计调用CFD计算208次。表2给出优化翼型和RAE2822气动力结果对比，优化构型减阻约31.74%。

表2 RAE2822和优化翼型气动性能对比Table 2 Comparison of aerodynamic performance of RAE2822 and optimal airfoil

图8和图9分别给出了优化翼型和RAE2822外形压力系数Cp分布对比曲线，c为弦长。可见优化翼型显著地减弱了激波强度，使阻力大幅降低。

图8 RAE2822和优化翼型外形对比Fig.8 Comparison of shapes for RAE2822 and optimum airfoils

图9 RAE2822翼型优化前后的表面压力系数分布对比Fig.9 Comparison of surface pressure coefficient distributions for RAE2822 airfoil before and after optimum

3.2 宽速域翼型减阻多目标优化

针对某宽速域翼型，为提高其跨声速和超声速气动性能进行多目标优化，优化目的是提高Ma=0.72和Ma=3.15时的最大升阻比K。

优化目标：

(5)

约束条件：

(6)

参数化方法采用改进的Hicks-Henne翼型参数化方法优化，上下翼面总计16个设计变量，通过求解定常可压RANS方程获得优化过程中翼型气动性能，湍流模型为k-ωSST两方程模型。优化算法采用NSGA-2多目标优化算法，每代种群有60个个体，总代数为50，总共完成3 000个翼型构型计算评估。

图10给出优化得到的Pareto前沿，共计89个优化构型，相比初始构型，Pareto前沿优化构型Ma=0.72、α=3.0°(KDP1)和Ma=3.15、α=6.0°(KDP2)下升阻比都有所增加，其中OPT1和OPT3分别为KDP2和KDP1最优构型。沿着Pareto前沿由构型OPT1到OPT3，优化构型Ma=0.72、α=3.0°升阻比越来越大，相比初始构型最大增加37.17%(OPT3)；反之，Ma=3.15,α=6.0° 升阻比数越来越大，相比初始构型最大增加11.39%(OPT3)，OPT2为KDP1和KDP2都增加较多优化构型，相比初始构型分别增加21.48% 和8.01%。

图11给出了优化翼型外形对比结果对比，图12 和图13给出了优化构型和初始构型Ma=0.72 和Ma=3.15下升力和升阻比对比结果。可见Ma=0.72下优化构型升力性能略有增加，Ma= 3.15下升力性能基本不变，满足约束，Ma= 0.72和Ma=3.15优化构型升阻比都获得了较大提升，达到了优化目的。

图10 宽速域翼型多目标优化Pareto前沿Fig.10 Pareto multi-optimal front of wide-speed range airfoil

图11 优化前后翼型外形对比Fig.11 Comparison of shapes for initial and optimum airfoils

图12 Ma=0.72时初始和优化翼型气动性能对比Fig.12 Comparison of aerodynamic performance of initial and optimal airfoils at Ma=0.72

图13 Ma=3.15时初始和优化翼型气动性能对比Fig.13 Comparison of aerodynamic performance of initial and optimal airfoils at Ma=3.15

3.3 多段翼型多目标优化

针对带扰流板下偏先进三段翼型着陆构型，开展气动外型和缝道参数优化，以提高其升力性能，计算条件为Ma=0.2、Re=45.6×106，优化问题描述如下：

(7)

通过椭圆方程控制生成多段翼型外型，包括6个外型设计参数，前缘缝翼和后缘襟翼缝道宽度、搭接量和偏角6个设计参数，以及扰流板下偏角，总计13个设计变量。优化采用NSGA-2每代种群有70个个体，总代数为50，总共完成3 500个构型计算评估。图14给出优化结束后，优化得到的Pareto前沿，相比初始构型，Pareto前沿优化构型迎角8°和27°下升力系数都有所增加，其中OPT1和OPT3分别为迎角27°和8°升力性能最优构型。沿着Pareto前沿由构型OPT1到OPT3，迎角8°时升力系数越来越大，反之，迎角27°升力系数越来越大，OPT2为综合迎角8°和迎角27°升力性能的优化构型。

图15给出了多段翼型典型优化构型外形对比，表3给出了典型优化构型设计点升力系数对比结果。相比基本构型，优化构型升力系数有明显增加，其中8°迎角升力系数最大增加 9.24%(OPT3)，27°迎角升力系数最大增加1.93%(OPT1)；从图16多段翼型典型优化构型升力系数随迎角变化曲线可以看出，构型OPT1有最大的最大升力系数，线性段增加较小，而构型OPT3线性段升力系数增加较大，最大升力系数增加较小，OPT2构型8°迎角和27°升力系数都有一定的增加，分别增加7.74%和1.43%。

图14 多段翼型多目标优化Pareto前沿Fig.14 Pareto multi-optimal front of multi-element airfoil

图15 优化前后多段翼型外形对比Fig.15 Comparison of shapes for baseline and optimum multi-element airfoils

表3 基础和优化多段翼型升力系数对比Table 3 Comparison of lift coefficient of baseline and optimal multi-element airfoils

图16 基础和优化多段翼型升力系数曲线Fig.16 Curves of lift coefficient of baseline and optimal multi-element airfoils

3.4 飞翼布局多目标优化

针对典型飞翼布局，开展跨声速减阻多目标优化设计，优化目的是提高跨声速设计点升阻比和阻力发散马赫数，优化问题描述如下：

(8)

式中:y为机翼展向站位；B为机翼展长。

参数化方法采用翼型剖面CST+CATIA二次开发三维成型方法，共计5个优化剖面(图17红色曲线),每个剖面18个设计变量，总计90个设计变量。采用代理优化方法，初始采样点270个，试验设计方法为LHS，代理模型为Kriging模型，采用Pareto优化时的EI加点方法，总计调用CFD计算1 200次。

图18给出了优化Pareto前沿，表4给出了典型优化构型与基础构型阻力系数对比结果，相比基础构型，Ma=0.72、CL=0.35时优化构型减阻13.17%，Ma=0.78、CL=0.30时优化构型减阻19.56%。从图19和图20不同马赫数下优化构型和基础构型上翼面压力云图对比结果可以看出，优化构型显著地减小了上翼面激波强度，从而使阻力大幅降低。从图21不同升力系数下阻力系数随马赫数变化结果可以看出，优化构型阻力发散马赫数获得了一定的增加，达到了优化目标。