0 引言

飞机/发动机集成(下文简称：飞发集成)是翼吊布局飞机研制中的关键技术。随着涵道比和尺寸的增大，发动机短舱相对机翼的近距耦合安装不可避免，发动机和机翼间的气动干扰问题加剧。研究表明^[1]：发动机短舱/吊挂对内、外两侧机翼前缘附近的“上洗”气流分别产生增强和削弱的效果。国内多位学者也对飞发集成构型的复杂流场特征开展了分析研究^[2-4]。实际应用中，安装效应对机翼压力分布和展向载荷分布的影响范围超过50%翼展^{[1, 5]}。因此，在机翼设计中直接考虑短舱/吊挂的干扰作用变得很有必要。德国宇航院Hohciscl和巴西航空工业公司Guilherme等分别在文献[5]和文献[6]中详尽阐述了在飞发集成的环境中开展部件气动设计的重要性，Bousquet^[7]总结了ONERA在发动机集成方面的风洞试验研究。

在设计工具方法方面，由于跨声速阶段阻力对几何外形变化的响应过于敏感，依靠人工修形的经验设计方法难以最大限度地挖掘减阻潜力。搭建自动的优化系统开展气动外形优化设计是近年来研究的热点领域，主要涉及几何外形参数化、计算网格自动生成、气动性能评估、优化搜索算法等方面。东京大学Hyoung-Jin等^[8-10]使用伴随梯度优化方法开展了针对机翼/机身/短舱/吊挂构型(下文简称：WBNP构型)的多点优化设计工作，CFD计算用的是Euler求解器。在WBNP复杂外形几何参数化和网格变形方法方面，Takanaka^[11]、Truong^[12]和左英桃^[13]等都做了大量有意义的研究。胡晓东^[14]等在民机短舱安装位置优化中得出了安装位置对全机气动特性的影响规律。清华大学张宇飞等^[15]尝试了运用多目标优化算法在短舱干扰下开展机翼一体化设计，有较高的工程应用价值，但由于计算量巨大，设计周期是个问题。

为了实现在飞发集成构型中直接优化机翼，本文建立了基于直接RANS计算分析和遗传算法寻优的多目标优化系统，针对NASA CRM的WBNP构型开展了机翼多目标优化设计工作。在翼型CST参数化方法的基础上发展了一种简捷高效的三维机翼参数化描述方法，并将几何参数化与网格变形集成在一起，避免了调用CAD软件建模和网格生成软件重新生成计算网格的复杂CFD前处理过程，从而达到了简化优化流程、提高鲁棒性的目的。此外，实现了在“天河2号”超级计算机上同时对超过100个方案实施分布式并行计算分析，再加上适当降低收敛标准和一系列CFD加速收敛措施的合理使用，较好地解决了优化周期问题。

1 优化方法 1.1 几何外形参数化

民用飞机机翼气动设计一般针对图 1所示的控制剖面展开，最终的机翼几何模型则在CATIA或其他CAD软件中基于这些控制剖面放样生成。本文以文献[16]中二维CST参数化方法为基础，发展了一种简捷高效的机翼参数化方法。

对于每个控制翼型，上表面或下表面的CST参数化表达式为：

$ \xi \left( \psi \right) = C_{1.0}^{0.5}\left( \psi \right)S\left( \psi \right) + \psi \left[ {{\xi _T} - \Delta {\alpha _T}} \right] $

(1)

其中弦向相对坐标：

$ \psi = \left( {x - {x_{LE}}} \right)/c $

(2)

法向相对坐标：

$ \xi = \left( {y - {y_{LE}}} \right)/c $

(3)

后缘开口相对大小：

$ {\xi _T} = \Delta {y_{TE}}/c $

(4)

类型函数：

$ C_{1.0}^{0.5}\left( \psi \right) = \sqrt \psi \left( {1 - \psi } \right) $

(5)

形状函数由n阶Bernstein多项式叠加构成：

$ S\left( \psi \right) = \sum\limits_{r = 0}^n {{w_r}\frac{{n!}}{{r!\left( {n - r} \right)!}}} {x^r}{\left( {1 - x} \right)^{n - r}} $

(6)

w_r为第r项的权重系数，不同的权重系数组合描述不同形状的翼型，也就是翼型的CST参数。已知翼型的m个离散点坐标(m>n)，使用最小二乘逼近，求解n+1阶线性代数方程组即可得到描述该翼型的n阶CST参数。

这里将机翼表面视为由沿展向的无限多个剖面翼型构成，任一展向位置z(y指向上，x指向后)处剖面参数由控制剖面参数插值得到。采用如下插值方法：机翼剖面翼型的前缘点y坐标y_LE、当地扭转角Δα_T和翼型n阶CST参数(w₀, w₁，…，w_n)由控制剖面以三次样条插值得到；前缘点x_LE、当地弦长c、后缘开度ζ_T则由附近两个控制剖面线性插值得到。这就是本文提出的机翼三维参数化方法。给定机翼平面形状内任一位置(x, z)，可以唯一地算出上、下表面y坐标y_U、y_L。

本文的优化设计中，设计变量是控制剖面的扭转角或CST参数。以上下表面同步连动的方式，实现了厚度不变、仅改变剖面中弧线的优化设计。设计变量以扰动量的形式叠加到描述初始外形(以下称Baseline外形)的基本参数上，基本参数需在优化开始前对所有控制剖面经过参数反算获得。

1.2 计算网格与网格变形方法

Baseline外形的计算网格在ANSYS ICEM-CFD软件中生成(如图 2所示)，由239块结构化网格构成，总网格单元数约800万。相应的机翼表面网格由110片构成。Baseline外形的网格将作为网格变形的输入。本文的网格变形方法是与1.1节的机翼参数化方法集成在一起的。在给定一组机翼参数后，对于Baseline机翼表面的每一个网格点，根据其在平面形状中的展向、弦向位置，应用本文的参数化方法可以算出其在新机翼外形中的坐标。这样就可以把Baseline机翼表面网格“投影”到新机翼表面上。接下来以指数衰减规律将表面网格的位移，传递到每个网格块的角点。最后用无限插值(TFI)方法插值得到内部网格点新的空间坐标。与此同时，机翼厚度、容积等几何信息也都计算出来，这些量用来判断是否满足几何约束条件。

在变形过程中，网格拓扑结构不变，空间网格点的分布规律也变化不大。实际应用中，该网格变形方法对机翼扭转分布和剖面形状的改变具有足够的鲁棒性，最重要的是，整个变形过程可在10 s左右完成。此外，表面和空间网格的变形可以分开进行，变形后的表面网格只相当于空间网格变形的“边界条件”。在优化过程中，只需向“天河2号”传输表面网格，文件大小不到空间网格的1%。这样可以显著降低存储系统读写操作强度，利于提高并行计算规模。

1.3 CFD计算分析和自动后处理

使用内部CFD程序SFlow对外形方案计算分析，该CFD程序在多块点对点结构化网格上，以有限体积法求解雷诺平均N-S方程。无黏通量项用Roe平均迎风通量差分分裂格式(FDS)离散，黏性通量项用中心差分格式离散，时间推进计算采用隐式LU-SGS方法。该程序有SA一方程和SST两方程湍流模型可选用。本文的CFD计算用SST两方程湍流模型。采用MPI方式并行和基于贪婪算法的并行任务分配策略具有相当高的并行效率。该CFD程序具有定升力计算功能，在迭代计算过程中，通过不断调整来流迎角，逼近设定的升力系数值。

网格分辨率和迭代计算收敛标准严重影响气动力计算时间和可信度。使用多大网格量，如何设定收敛标准，是优化设计前必须解决的问题。图 3为针对WBNP构型，网格单元数分别为200万、800万、1600万和6400万的4套网格在马赫数0.85、升力系数0.5下的网格收敛性曲线。单元数1600万的网格是作者在气动性能计算工作中经常使用的，在这里与6400万网格计算的阻力相差不到3 counts。单元数200万的网格由1600万的网格经过半粗化得到。单元数800万网格是在1600万的网格基础上，适当降低机身、短舱/吊挂以及远场网格分辨率后得到，机翼附近网格分辨率基本不变。

图 4为使用单元数800万的网格，以收敛好的Baseline流场解为初场，扰动后外形CFD计算的收敛历史。在“天河2号”上使用48核并行，同时使用多重网格和当地时间步长等加速收敛措施，13 min左右即可完成600次迭代，阻力系数和力矩系数的收敛程度可以满足要求。经过计算精度和时间的权衡，下文的优化设计中使用单元数800万网格，而验证计算则使用1600万的网格。

CFD计算结束后，为优化设计而发展的后处理程序从流场计算结果中自动提取流场特征量和气动力系数。这里主要实现了指定截面压力分布特征量(吸力峰、激波位置等)和展向载荷分布、升力系数分布(图 5)的提取分析。这些特征量和力系数都可作为优化设计的目标或约束。

1.4 超级计算机上的优化设计流程

本文的优化设计系统以Windows/Linux跨平台的方式运行于“天河2号”超级计算机(下文简称“超算”)上。具体来说是以操作系统为Windows Server的胖节点(32核心，64线程，256GB内存)为整个优化设计系统的控制台，完成计算任务预处理，连接登录节点提交计算任务，以及优化搜索计算。登陆节点和计算节点操作系统为Linux，分别负责计算任务排队提交与计算结果回传、CFD计算与后处理。经过反复测试，这样的系统架构运行稳定，超算系统本身承受的压力较小。

具体优化流程如图 6所示。每个新个体都是在Baseline外形的基础上，叠加由寻优算法给出的扰动量得到，相应的表面网格也会做相应的扰动变形。与此同时，新外形的几何特征也会被分析出来，并判断是否满足几何约束。如果不满足几何约束，该个体将被淘汰，不再进行后续的CFD计算分析。这样做的好处在于尽早剔除不满足约束的个体，避免无效的CFD计算并加快优化进程。对于满足几何约束的个体，其表面网格将被传输给“天河二号”，通过登录节点提交计算作业，在计算节点上依次完成空间网格变形、CFD计算和自动后处理过程。最后，计算结果传回胖节点，判断是否满足气动约束，根据目标函数值判定个体优劣。完成一代种群的分析后，优化算法会生成新一代个体。本文的优化案例采用具备精英策略的非支配排序遗传算法NSGA-Ⅱ实现多目标寻优。

本文开展的基于大规模并行计算的民机复杂外形优化设计在超算应用领域属于高通量类型计算，即单个计算作业占用资源较少(本文中使用48核持续不到15 min)，但同时进行的计算作业数量庞大(数百甚至数千个)。每个计算作业的提交和监控都必须通过登录节点完成，每个作业都需要拷贝文件、新建文件和读写文件。高通量应用最为考验超级计算机的综合性能，对计算资源、I/O系统、网络系统、文件存储和文件管理等系统资源消耗巨大。超算的登录节点是一个公共平台，主要用于管理员对超算环境维护管理，及所有用户提交计算作业和上传下载数据，有时还作为各种自编程序的编译环境。因此登录节点资源不允许单个用户长时间或大量占用。

研究发现，高通量计算应用的主要瓶颈在于登录节点的线程处理能力和整个超算平台的I/O系统。若要提高计算规模，必须设法降低单个计算作业消耗的系统资源及通讯量。采用Python语言编写脚本，在胖节点上实现远程ssh连接、计算作业提交及相应文件上传。计算任务提交完毕之后，关闭对应ssh连接。计算完成后，再由Python脚本调用登录节点本地的shell脚本将计算结果下载回胖节点。在缓解I/O系统压力方面，通过代码优化、流程改进等措施尽可能减少文件创建以及读写操作，采用这种“登录节点和I/O系统减负”方案后，登录节点CPU占用率和内存消耗明显减少，优化系统运行更加稳定。

2 优化案例 2.1 CRM安装效应分析

NASA CRM是以巡航马赫数0.85的宽体客机为背景设计的标准模型，设计过程是先在翼身组合体构型(下文简称：WB构型)中设计一副性能优良并且足够鲁棒的机翼，以保证集成短舱、吊挂后性能损失不大。为了降低网格生成难度，动力装置用一个单圈通气短舱模拟，其流量系数与典型商用飞机巡航时相当。WB构型巡航设计点为马赫数Ma=0.85，基于平均气动弦长的雷诺数Re=4×10⁷，升力系数C_L=0.5。此工况下WB构型和相同迎角(1.875°)下WBNP构型表面等压线、机翼剖面压力分布的CFD计算结果对比分别如图 7、图 8所示。图 8中前两处剖面位于吊挂内侧，后两处位于外侧。集成短舱吊挂后，内翼上表面在附加上洗作用下吸力峰抬高，激波前移增强，内翼下表面在机翼、吊挂和短舱构成的通道内显著加速。吊挂外侧机翼上表面在下洗作用下吸力峰降低。根据设计经验，内翼激波增强会对阻力发散特性不利。

图 9对比了短舱、吊挂集成前后机翼展向载荷分布，图中横坐标Eta为机翼展向相对位置，纵坐标Load为当地升力系数与弦长的乘积。在迎角不变的情况下，WBNP构型机翼升力显著损失的范围超过60%翼展。而且在吊挂的阻断作用下，机翼载荷分布在吊挂位置存在间断，这是在飞发集成构型设计中需要考虑吊挂的一个重要原因。CRM的设计初衷除了用于验证数值计算工具和标定风洞外，还将用作考察优化设计工具的标准模型。WB构型的机翼载荷分布被设计成十分接近椭圆，因此就单独WB构型来说，很难在机翼厚度不减的情况下同时在多个状态点上取得显著减阻效果。而集成了短舱吊挂后的CRM在激波阻力和诱导阻力方面都应该有可降低的潜力。

2.2 机翼弯扭分布优化设计

平面形状确定后，机翼高速气动性能在很大程度上取决于厚度分布和弯扭分布的组合。厚度主要影响激波阻力，在内部容积、结构布置等条件的约束下，厚度分布往往没有太多缩减的余地。通过改变机翼弯扭分布可以调整展向升力分布和载荷分布，进而影响诱导阻力、俯仰力矩以及低速失速特性。使用本文的优化系统，在厚度不减的前提下对CRM WBNP构型的机翼扭转分布和控制剖面中弧线形状进行优化设计。优化问题定义如下：

$ \begin{gathered} \min \;\;{C_{D\_Ma = 0.83}}\;\& \;{C_{D\_Ma = 0.85}}\;\& \;{C_{D\_Ma = 0.87}} \hfill \\ \;\;\;\;\;\;\;\;{C_L} = 0.5 \hfill \\ {\text{s}}{\text{.t}}{\text{.}}\;\;\;\;{C_{m\_Ma = 0.85}} \geqslant - 0.145 \hfill \\ \;\;\;\;\;\;\;\;{\alpha _{\_Ma = 0.85}} \leqslant 2.2° \hfill \\ \end{gathered} $

(7)

本次实施的是3点3目标优化，3个设计点分别为马赫数0.83、0.85和0.87，每个设计点的阻力系数作为独立的目标函数。设计变量包括9个控制剖面的扭转角和表达中弧线的8阶CST参数，共90个。CST参数的变化范围设定为在Baseline值基础上扰动±0.05。翼根剖面的扭转角变化范围±0.5°，为保证扭转曲线单调下降，其余控制剖面的扭转角定义为在内侧相邻剖面扭角绝对值基础上的下降量，变化范围设定为在Baseline基础上扰动±0.5°。这些量值的设定依赖于大量应用后的经验总结，而且在下文的优化过程中，没有发现变量向边界收敛的现象。

设计约束包括：升力系数C_L=0.5；设计点Ma=0.85的俯仰力矩系数和迎角。在CFD程序定升力计算功能的帮助下，可直接实现升力系数的等式约束。计算评估后，不满足力矩系数和迎角约束的个体将会被列为不可行方案，在排序中不可行方案的等级低于所有可行方案。随着种群逐渐进化成熟，这种不可行个体所占比例也会越来越少。

本次优化种群规模为256，每批次同时对128个个体进行CFD计算，每个个体串行完成3个设计点的计算。在60 h内，完成超过10 000个个体的计算分析，遗传优化40代。图 10为优化过程中3个目标函数的演化历史，从下沿轮廓来看，优化已经趋于收敛。

多目标优化问题以找到Pareto前缘作为结束。最终在多个目标间权衡选择还需要由设计者根据应用需要完成。分析优化结果发现，马赫数0.83和0.85两个目标间基本不构成竞争关系。因此最优个体的遴选主要在马赫数0.85和0.87两个目标间的Pareto前缘上进行。图 11展示了优化过程中成功完成计算的个体的马赫数0.85和0.87两个目标函数值分布，蓝色方形标志为Baseline方案所在位置，在Pareto前缘上选定三个红色菱形标志的个体为此次优化结果的备选方案。图 12对比了他们和Baseline方案的表面等压线图，三个备选方案的内翼激波强度都明显降低。

表 1列出了他们与Baseline方案的阻力系数对比。ID号9740的个体在马赫数0.85时阻力最低，比Baseline方案降低3.1 counts，但在马赫数0.87时阻力仅降低3.3 counts。9415号个体在马赫数0.87时阻力下降最多，达到14.2 counts，但在马赫数0.85时阻力仅降低0.6 counts。7777号个体在马赫数0.85时阻力降低2.3 counts，在马赫数0.87时阻力降低10.4 counts，各设计点阻力下降比较均衡，下面重点分析该个体的细节特征。图 13对比了优化前后机翼扭转分布的差异，除了翼梢扭转角有所增加，大部分区域扭转角基本保持不变。图 14则对比了优化前后9个控制剖面的形状变化。

图 15和图 16分别对比了马赫数0.85和0.87优化前后四处剖面压力分布的变化。优化后，在马赫数0.85时，内翼部分激波强度明显减弱，吊挂外侧则表现为吸力峰提高。马赫数0.87时，各剖面激波后移减弱。

2.3 优化结果验证

本文的机翼三维参数化方法与CATIA中的放样成型不可能完全相同，同时优化过程中网格分辨率和收敛判据也有所降低，因此优化取得阻力下降的效果还需要进一步验证。将优化后的控制剖面在CATIA软件中放样成型，得到优化后机翼的CAD模型。在ICEM-CFD软件中重新生成单元数1600万的计算网格开展验证计算。计算从自由来流初场开始迭代，残值下降4个量级后提取计算结果。图 17为优化前后阻力发散曲线的对比，可见优化过程中所取得的减阻成果基本得到验证。

3 结论

探索了应用超级计算机开展民机飞发集成构型精细优化设计。直接使用RANS计算评估，同时在近距耦合部件间的干扰作用下开展全局寻优的多目标优化设计，得到的减阻效果可靠。主要工作简述如下：

1) 搭建了可用于民机WBNP构型的优化设计系统，采用遗传算法进行全局寻优，包括CFD前处理、RANS计算的CFD分析和自动后处理等功能模块。该系统在超级计算机上可同时对数百个设计案例开展计算分析。并行计算、加速收敛和适当降低收敛标准等技术的综合运用使优化过程能在可接受的时间周期内完成。

2) 在厚度不变的情况下，对CRM WBNP构型中的机翼弯扭分布开展了多目标优化设计。60 h内完成了40代遗传优化，各设计点阻力系数有2~10 counts的下降，阻力发散特性显著改善。验证计算结果表明，优化过程中取得的阻力下降量是可靠的。

致谢: 感谢广州超算中心对本文工作的大力支持。