通用飞行器气动优化设计数字化集成平台——DIPasda

2020-06-03孙俊峰周铸黄勇庞宇飞卢风顺许勇

航空学报 2020年5期

孙俊峰，周铸，黄勇，庞宇飞，卢风顺，许勇

中国空气动力研究与发展中心计算空气动力研究所，绵阳 621000

随着计算机能力和数值计算技术的不断发展，数值优化设计方法在工业设计部门得到越来越广泛的应用[1-4]，改变了传统设计过程中，工程师用试凑法来进行设计的不足,提高了设计效率。随着对飞行器设计要求的不断提升，多学科优化设计(MDO)开始流行，设计师希望用精确、高效、鲁棒的多学科设计平台环境，来减少设计周期和设计成本，同时改善飞行器的设计质量。如美国波音公司开发的多学科优化软件MDOPT[5]，采用图形用户界面方式，集成了各类鲁棒数值优化算法、高阶CFD分析软件，可进行飞机气动、重量、性能及稳定与控制等学科的综合优化设计。多学科优化设计软件iSIGHT通过集成各种仿真软件和并行优化算法，提供工作流模式的设计流程定制方案，目前已在航空、航天、船舶、电子、汽车等多个工业部门获得广泛使用。空客公司采用iSIGHT并集成结构力学和空气动力学软件，在保证气动弹性约束条件下使A320机翼重量减小3%。此外还有NASA格伦研究中心的OpenMDAO[6]以及欧洲的开源框架项目FP7和Horizon 2020等，都极大提高了优化设计的易用性和优化效率。

引入优化设计框架的目的是为多学科优化设计应用提供支持。理想的多学科优化设计架构应满足以下设计需求：① 可伸缩性，用户能够根据现有的计算资源构建优化设计应用；② 灵活性，用户能够针对优化问题、优化策略选择不同的优化算法、解算器等设计工具；③ 可扩充性，系统能够方便、灵活地进行功能扩充。

本文第1节介绍了系统架构设计和系统的工作流程。第2节详细介绍了系统的功能模块，包括优化算法模块、几何管理模块、代理模型模块、学科分析模块等。第3节通过优化设计算例来展示系统进行多目标优化设计的能力。最后对本文进行了总结。

1 架构设计

作为一个数字化多学科综合设计和评估的集成平台，DIPasda 平台系统采用软件工程的设计方法，根据高内聚、低耦合的原则分成层次结构，包含图形用户界面层、MDO功能模块管理层、数据传输(TCP/IP,SSH)层以及学科管理层。图1给出了系统架构示意图。用户通过图形用户界面与系统进行交互操作，功能模块管理层是系统的主要部分。用于建立优化设计问题，对不同学科提供初始的输入参数。此外还控制着系统的优化设计流程以及提供优化设计过程中所需的工具。MDO功能模块管理层包含问题定义、代理模型、几何管理、优化算法以及输入/输出等功能模块。数据传输(TCP/IP,SSH)层维护底层的数据传输控制。用户通过用户界面输入的参数信息、优化信息以及相关的工程信息等都被打包进脚本文件，通过脚本文件在系统中进行传递。脚本文件用XML写成。XML文档可以跨平台、跨应用进行数据的交换，非常适合作为数据交换的标准。学科管理层提供对给定外形进行不同学科性能分析的功能，当前的DIPasda系统包含了空气动力学和气弹2个学科分析功能，可以完成气动多目标综合设计的相关功能，同时系统留有学科扩展功能接口，用户可以借助接口方便地进行学科的功能扩展，用以完成多目标/多学科的优化设计。

图1 系统架构示意图Fig.1 Diagram of system architecture

DIPasda系统既可以进行基于进化算法的全局多目标优化设计，也可以采用伴随/耦合伴随方法进行基于梯度信息的优化设计。图2给出了DIPasda系统进行全局优化设计流程的示意图，展示了使用DIPasda进行全局优化设计的关键步骤。按照顺时针方向，首先原始的几何模型进入系统，系统根据需要选择合适的模型参数化方法，进行几何模型的参数化，提取设计变量，定义设计空间。参数化完成后用户可以根据优化需要选取合适的优化工作流，如梯度类方法或是全局类的优化算法，系统会分别进入不同的工作流程。如果选择全局类方法，要根据设计变量和设计空间确定实验设计方法，定义采样点集合，对每个采样点进行不同学科的网格生成和学科分析，用学科分析得到的数据构建代理模型，代理模型主要用来替代耗时的学科分析数据，提高优化效率。用构建的代理模型分析优化问题的优化目标和约束，然后进入优化流程，经过多轮迭代后得到最终的优化结果。

图2 DIPasda系统全局优化设计流程示意图Fig.2 Diagram of DIPasda system global optimization process flow

目前，DIPasda系统已经基于PC和集群系统完成了2.0版本的开发。系统采用C++语言进行架构开发，学科解算器采用C++语言/FORTRAN语言进行编程，XML文件用于交换信息。基于Windows的用户界面采用C++语言和wxPython库开发完成。用户可以在PC机上完成优化问题的描述和参数的输入后，通过集群系统完成优化设计的流程。

2 系统模块

DIPasda 系统主要由以下功能模块构成，分别是：优化模块、几何管理模块、代理模型模块、学科分析模块、问题描述模块、图形用户接口(GUI)模块以及数据库管理模块。

2.1 优化模块

优化模块由优化算法库组成。优化算法库包含有多种局部和全局的寻优工具，例如伴随方法[7-9]、进化算法[10-11]、模拟退火算法、序列二次规划法等，用户可以通过应用界面来选择不同的优化算法完成鲁棒设计、多学科/多目标设计以及基于主成分分析(Principle Component Analysis，PCA)的多目标降维处理等功能。

许多真实世界的优化问题都涉及同时优化多个目标，在多目标优化问题中，设计者通常希望找到Pareto最优解。在当前的DIPasda系统中，采用进化算法结合代理模型是通常的多目标优化模式。进化算法是基于种群的全局搜索策略，遵从达尔文的自然选择理论。按照进化理论，种群中的个体会随着进化过程逐渐适应环境。进化算法从随机种群开始，对种群中的每个个体根据目标函数和约束函数计算适应值。一般进化算法涉及3种进化算子操作：选择算子、交叉算子和变异算子。选择算子用来选择参与产生后代的父代个体，常用的选择方法有锦标选择、按比例选择以及排序选择等方法。对多目标优化问题，排序通常基于非受控排序。种群排序后，通常采用精英策略保证优良的个体直接进入下一代种群。交叉算子和变异算子都用来产生后代并保持种群的多样性。通过进化算子的作用，种群不断进化直至得到最终的Pareto前沿。图3给出了进化算法的流程图。

基于灵敏度的求解主要包括复变量、有限差分法、符号微分等方法。Sobieszczanski-Sobiesk[12]系统总结了灵敏度分析方法在飞行器优化设计领域的应用。Jameson等[7]提出的基于变分思想的伴随方法，由于其与设计变量个数无关的优势，在灵敏度分析中扮演着重要角色，因此在气动优化设计中得到广泛应用。

对于气动优化设计的最小化问题[13]：

minI(W,X)

(1)

式中：I为目标函数；W为流场守恒量；X为设计变量。考虑残差约束R(W,X)=0，引入伴随算子Λ,可得目标函数表达式为

L=I+ΛTR

(2)

图3 进化算法流程图Fig.3 Flow chart of evolutionary algorithm

(3)

式(3)即流场伴随方程，采用式(2)和式(3)的处理方式以及多个学科残差约束，可以推导出相应的多学科伴随方程，求解Λ之后，进行目标函数梯度信息快速求解。

(4)

(5)

式中：ΔX为设计变量梯度。图4给出了基于伴随/耦合伴随的优化设计流程图。

2.2 几何管理模块

几何管理模块对几何外形操作功能进行管理，包括外形的输入/输出、约束评估、模型参数化和网格重构等。几何模型参数化过程应当尽可能自动化并对所有的学科分析提供一致的参数化表示。常用的几何模型参数化方法包括：NURBS(Non-Uniform Rational B-Splines) 曲线/曲面[13]、CST(Class function/Shape function Transformation)[14]、自由变形(Free-Form Deformation,FFD)[15-16]等。用户可以使用这些参数化方法表示飞行器翼型、机翼以及全机外形。参数化给出了优化设计变量和设计空间。

1) NURBS 曲线/曲面

NURBS 曲线可以表示为

(6)

图5 基于NURBS的参数化方法和控制点示意图Fig.5 Diagram of variables definition and control point based on NURBS

2) 自由变形(FFD)

FFD技术与计算机动画中的变形技术非常相似。可以模拟平面形状、扭转、上反、厚度和弯度变化等外形改变。FFD技术适合于飞行器设计的初步设计阶段。

FFD 技术可以表示为

(7)

FFD方法具体实现的基本步骤如下：

步骤1 选择适当的变形技术与变形体。在变形控制体坐标(ξ,η,ζ)与基本外形网格坐标(x,y,z)间建立映射关系。变形控制体可以是平行六面体、非平行六面体、Bezier体或NURBS体。这一步将待设计的部件用控制体包住，给出控制体的描述。

步骤2 建立基本外形网格坐标(x,y,z)到变形控制体坐标(ξ,η,ζ)的反映射。(ξ,η,ζ)这些映射参数保持固定，独立于形状变形，作为预处理工作，只需要求一次。也就是说，求出待设计部件表面网格点(x,y,z)在变形控制体里的对应坐标(ξ,η,ζ)。

步骤3 扰动变形控制体的控制参数，即设计变量，得到新的变形体。

步骤4 计算待设计部件表面网格点(ξ,η,ζ)在新变形体中的坐标(x,y,z)实现待设计部件的变形控制。

3) CST 方法

CST 方法由波音公司在2006年提出,几何外形表示为类别函数和型函数的乘积,即

(8)

类别函数表示为

(9)

指数N1、N2从 0～1可以产生任意外形。对于NACA对称翼型，指数N1、N2通常取0.5和1。

型函数S(ψ) 表示为

(10)

式中：Si(ψ)为i阶Bernstein多项式；Ai为多项式系数。

4) 计算网格重构

优化设计过程中，外形的改变通常需要重新生成计算网格，在自动化的优化设计流程中，网格重构是一项挑战性的工作，需要保证网格的质量，以满足性能评估的需要。

为了与几何外形的改变相适应。系统采用基于框架的方法[17]实现了对多块结构网格的网格重构。基于框架的网格重构方法首先从初始网格中提取网格的框架关系，保存多块结构网格的网格连线和界面网格的信息。这些信息可以保证网格重构时的拓扑关系。物面外形改变后，与改变物面相关联的框架线也随之变化，对变化后的网格框架可以采用超限插值(TransFinite Interpolation，TFI)的方法重新生成体网格，从而完成网格的重构。图6给出了AGARD445.6 机翼设计中面网格和框架变形的示意图。

图6 AGARD445.6机翼表面网格和框架变形示意图Fig.6 Framework of wing AGARD445.6 and its transformation diagran

2.3 代理模型模块

优化设计过程中，采用高精度学科分析解算器进行目标函数和约束函数的评估，需要耗费大量时间和资源成本。一种替代方法就是构造高精度学科分析解算器的代理模型，用代理模型来表征设计变量和目标函数之间的关系。使用代理模型可以快速评估目标函数和约束函数，从而可以节省时间提高优化效率。

DIPasda系统中，代理模型模块包含2部分内容：试验设计和代理模型方法。为了保证采样点在设计空间的均匀分布，系统可以采用拉丁超立方采样[18]和均匀采样[19]的方法进行试验设计选取采样点。Kriging模型[20]则用来作为目标函数和约束函数的近似模型。

Kriging模型是从地理统计学科发展起来的一种用来预测未知点函数分布的模型。通过Kriging模型可以得到未知点的函数值和不确定性，利用这些值，可以实现局部和全局搜索的平衡，这也就是期望改进(Expected Improvement，EI)准则概念，通过选择EI最大值的点加入到代理模型，可以改进模型的精度。图7给出了基于代理模型的优化流程示意图。

步骤1 利用拉丁超立方采样或者均匀采样构建初值均匀分布的采样点。

步骤2 利用高精度学科分析模块对采样点进行性能分析，构建针对目标函数和约束函数的Kriging代理模型。

图7 基于代理模型的优化流程示意图Fig.7 Flow chart of surrogate-based multi-objective optimization

步骤3 基于Kriging代理模型进行优化流程。

步骤4 基于EI准则添加采样点，重新构造代理模型。

回到步骤3重新进行优化，直至得到最终的最优解。

2.4 学科分析模块

DIPasda 系统气动分析模块和气弹分析模块，可以完成通用飞行器气动外形的综合设计，同时系统保留了学科的扩展接口，用户可以根据需要添加不同的学科分析模块，添加的接口可以实现与系统数据流的无缝对接。

气动分析模块主要用来提供气动力计算数据。DIPasda系统目前包含的气动解算器有： MBNS2D、Cart3D[21]、PMB3D[22]。这些解算器可以对飞行器外形完成基于欧拉方程和Navier-Stokes方程的流动分析。可以处理结构网格、非结构网格、笛卡尔网格以及重叠网格等多种网格形式。

DIPasda系统支持静气弹分析，针对复杂外形建立了基于Navier-Stokes方程和多块结构网格的高效、精确的静气弹分析方法。位移传递技术和载荷传递技术保证了耦合气动和结构计算的精度。发展了基于推进图的三维动态网格生成技术，实现了网格大幅变形后的计算网格的快速重构，同时保持网格质量满足计算精度的需求。

2.5 问题描述模块

问题描述模块用来定义优化问题的详细参数。指定设计变量及其范围，确定优化问题的目标函数和约束函数，设置学科分析模块的状态参数以及优化流程的控制参数。系统根据用户输入的参数，完成工作流参数的定义以及优化流程数据流的准备。图8给出了优化目标的定义和指定设计变量范围的操作示意图。用户通过图形界面完成优化问题的描述，相关信息通过XML文件在优化系统各功能模块之间进行传递。

图8 问题定义模块示意图Fig.8 Diagram of optimization problem definition

2.6 图形用户接口模块

图形用户界面是用户与系统交互的接口。DIPasda系统的各项功能都可以通过用户界面操作完成，系统提供了基于Windows的图形操作界面，适应用户的使用习惯。用户在前台通过操作界面完成优化问题的定义以及各项参数的设置，系统底层封装了与后台集群系统的连接以及信息的交换，可以实现本地计算与集群系统计算的无缝切换，用户可以通过输入/输出系统监视优化进程，方便处理各类操作问题。图9给出了系统应用界面的示意图。

图9 系统应用界面示意图Fig.9 Diagram of graphical user interface

2.7 数据库管理模块

数据库管理模块主要用于管理系统生成的数据，包括优化问题的定义文件、设计变量的定义、设计空间的范围、用户通过用户界面输入的学科分析的状态参数文件、优化参数的设置文件以及优化结果等数据。通过数据库管理模块，可以掌握优化流程中的关键信息，便于优化应用的交互管理和数据的保存。

3 应用算例

3.1 NACA0012翼型减阻优化

以AIAA气动优化设计讨论组(ADODG)给出的NACA0012翼型跨声速无黏流减阻优化算例来考核优化设计系统功能模块的性能。

优化设计问题描述如下：

以后缘修型后的NACA0012翼型作为初始翼型进行减阻优化设计，后缘修型为

0.351 6x2+0.284 3x3-0.103 6x4)

(11)

设计状态为：马赫数Ma=0.85无黏流，迎角α=0°，优化目标为阻力系数最小，几何约束条件为翼型厚度不减。优化问题可以描述为

(12)

式中:ybaseline为初始翼型的纵坐标值。

采用CST方法对翼型表面进行参数化，上下表面各取12个控制点作为设计变量。计算网格规模为769×129，如图10所示。采用进化算法进行翼型的减阻设计研究，种群规模240，进化100代，图11给出了优化的收敛历程，在50代左右已达到收敛。

表1给出了优化前后翼型的性能比较，在保持升力系数不变和厚度不减的约束下，阻力系数CD从0.047 2降到了0.004 2。图12给出了优化前后翼型的形状和表面压力分布的比较，图13给出了优化前后翼型马赫数云图比较。可以看出，优化翼型前缘半径增大，前后缘变钝，前后缘外形更加饱满。由于厚度约束的关系，翼型变厚，上表面更加平坦，激波位置从3/4弦长处移至翼型后缘，激波强度减弱，压力恢复更加缓和，降低了激波阻力。该算例展示了优化方法和参数化技术处理局部大变形和解算器的求解精度能力。

图10 计算网格Fig.10 Computing grid

图11 目标函数优化收敛历程(NACA0012)Fig.11 Convergence history of objective function optimization (NACA0012)

表1 基本翼型和优化翼型性能比较(NACA0012)Table 1 Comparison of characteristics between baseline and optimized airfoils(NACA0012)

图12 优化前后翼型形状和表面压力分布比较(NACA0012)Fig.12 Comparison of configuration and pressure distribution between initial and optimized airfoils (NACA0012)

图13 优化前后翼型马赫数云图比较 (NACA0012)Fig.13 Comparison of Mach number contour between initial and optimized airfoils (NACA0012)

3.2 RAE2822翼型减阻优化

考虑RAE2822翼型跨声速减阻优化算例。优化设计问题描述如下：

以RAE2822为初始翼型，进行减阻优化设计，设计状态为：Ma=0.734,CL=0.824，Re=6.5×106，优化目标为阻力系数最小，几何约束条件为翼型面积不减。优化问题可以描述为

(13)

采用CST方法对翼型表面进行参数化，上下表面各取13个控制点作为设计变量。多目标进化算法的参数设定为：种群规模200，交叉概率0.85，变异概率0.2，进化100代，图14给出了优化的收敛历程。计算网格规模为769×129。

表2给出了优化前后翼型的性能比较，在保持升力系数不变、力矩约束和面积不减的约束下，阻力系数降低了约87 counts。图15分别给出了优化前后翼型的形状和表面压力分布的比较，图16 给出了优化前后翼型压力云图的比较。可以看出，优化翼型消除了原始翼型在55%弦长位置附近的强激波，等值线变得平顺光滑。

图14 目标函数优化收敛历程(RAE2822)Fig.14 Convergence history of objective function optimization (RAE2822)

表2 基本翼型和优化翼型性能比较 (RAE2822)Table 2 Comparison of characteristics between baseline and optimized airfoils (RAE2822)

图15 优化前后翼型形状和表面压力分布比较 (RAE2822)Fig.15 Comparison of configuration and airfoil pressure distribution between initial and optimized airfoils (RAE2822)

图16 优化前后翼型压力云图比较 (RAE2822)Fig.16 Comparison of pressure contour between initial and optimized airfoils (RAE2822)

3.3 某机翼布局优化设计

对某翼身组合体外形进行机翼的优化设计，提高设计状态的升阻比。

1) 设计状态：

Ma= 0.35，Re=8.0×106

2) 优化目标：

MaximizeCL/CDatCL=0.98

3) 约束条件：

① 机翼最大厚度不减小。

② 保持俯仰力矩系数基本不变。

用CST方法进行机翼的参数化，采用Bernstein基函数的形式沿流向和展向统一构造机翼外形。沿机翼展向选取6个剖面构造CST的型函数，选取24个设计变量参与优化。图17给出了采用CST方法描述机翼变形的示意图。

机翼表面形状发生改变，与其相对应的空间结构网格块采用多块线对接结构网格的重构方法重新构造。半模的网格规模在1 000万左右，采用PMBNS3D软件进行气动力性能计算，湍流模型采用k-ωSST方法，利用多重网格技术进行收敛加速，在IBM工作站集群上进行采样点的并行计算。

采用基于Kriging代理模型的优化方法进行机翼的减阻优化设计，Kriging模型初始样本数目108，采用EI准则加点80次，图18给出了优化的收敛历程。表3给出了优化外形与基本外形气动力计算结果的比较，优化外形在保持机翼厚度不减和低头力矩基本不变的条件下，升阻比提高了2.2%，全机阻力系数减少了11 counts。

图19给出了优化外形和基本外形全机表面压力分布云图的比较，图20给出了优化外形的机翼和基本外形机翼4个剖面的几何外形和压力分布的对比，从几何对比上可以看出，优化翼型上表面前缘稍微变厚，下表面前缘变平，增加了机翼剖面的弯度。从压力分布的对比来看优化外形前缘负压降低,使得低头力矩增加。

图17 机翼CST参数化示意图Fig.17 Diagram of CST parameterization of wing

图18 目标函数优化收敛历程Fig.18 Convergence history of objective function optimization

表3 基本外形和优化性能比较Table 3 Comparison of characteristics between baseline and optimized airfoils

图19 优化外形和原始外形全机表面压力云图比较Fig.19 Comparison of pressure contour between initial and optimized wing

3.4 CRM宽体飞机气动设计

对宽体飞机CRM标模[23]巡航状态升阻比以及力矩特性进行多点综合优化，设计状态为：Ma=0.85，Re=5.0×106。

1) 优化目标：

图20 机翼沿展向4个剖面几何形状和压力分布的比较Fig.20 Comparison of shape and pressure distribution of four profiles along wing

2) 约束条件：

采用基于NURBS基函数的FFD方法进行参数化，图21给出了参数化示意图，共采用200个控制点实现了机翼的参数化建模。参数化需要在翼根、拐折位置以及翼尖处保持厚度几何约束。

基于序列二次规划(SQP)算法对多点开展伴随方法加权优化。作为对比，对巡航状态同时进行单点优化，两种方法均经过20代优化，图22给出了不同优化设计收敛历程。表4给出了单点优化和多点优化气动特性的比较，在保持几何约束条件下，多点优化设计的升阻比和阻力发散特性有明显改善，单点优化全机阻力系数减少了13 counts，升阻比提高了约4.7%。图23给出了展向Y=5,10,15 m站位压力分布的对比，可以看出激波强度均大幅减弱，单点设计与多点设计压力分布形态区别主要在外翼段。图24给出了单点优化外形和多点优化外形与初始外形压力系数云图的对比，可以看出单点优化消除了表面激波，多点优化呈现弱激波状态。文献[23]对该标模进行了更深入的优化和分析。

图21 CRM机翼FFD参数化示意图Fig.21 Diagram of FFD parameterization of CRM wing

图22 目标函数优化收敛历程比较Fig.22 Comparison of convergence history of objective function optimization

表4 不同优化进程气动特性比较Table 4 Comparison of aerodynamic characteristics between different design process

图24 单点和多点优化压力系数云图比较Fig.24 Comparison of pressure coefficient contours before and after sing-point and multi-point optimized configuration

3.5 飞翼布局机体/推进一体化设计

随着CFD技术和高性能计算机的发展，考虑进排气影响的飞行器机体-推进系统一体化设计开始成为可能。对于采用背负式动力形式的新概念民用飞机气动布局，推进系统的进排气效应对升力面流动形态的影响更为强烈，对设计参数的灵敏度产生重要影响，从而影响精细化设计效果。

对某飞翼布局外形[24]，进行考虑进排气影响下的减阻优化设计。优化目标为

MinimizeCDatMa=0.85,CL=0.36

采用序列二次规划算法进行带约束伴随优化，图25给出了飞翼布局外形的参数化示意图，其中1～8分别表示设计剖面控制点的位置。灵敏度采用伴随梯度计算方法，图26给出了有/无动力条件下设计变量灵敏度对比，可以看出动力效应对灵敏度产生了一定的影响，图27给出了一体化优化前后的压力云图的比较。可以看出优化后机翼激波强度明显减弱，验证了一体化设计的有效性，要充分挖掘飞行器气动设计的潜力，必须考虑推进系统动力影响，进一步发挥气动/推进一体化设计的优势。