伴随压力分布反设计方法在大型客机气动优化中的初步探索

刘峰博，蒋城，马涂亮，梁益华

1. 中国航空工业西安航空计算技术研究所，西安 710065 2. 中国商飞上海飞机设计研究院，上海 201210

国产大型客机研制是关系国家经济产业升级和创新技术发展的一项重大战略举措。而大型民用客机在“四性三减”(四性：安全性、经济性、环保性、舒适性；三减：减重、减阻、减排)上的极致追求，是提高设计能力和市场竞争力的核心驱动力。在“四性三减”极致性能要求下，民用客机的气动设计思路、方法及工具都需要得到升级和改进。

当下基于计算流体力学(Computational Fluid Dynamics, CFD)的分析及设计技术已大量应用于大型客机研发过程中，并且已大大缩减了设计周期及成本，波音、空客等民机设计力量均拥有成熟的CFD设计程序。在CFD的发展带动下，目前基于伴随理论的梯度优化设计方法与基于遗传算法的非梯度优化设计两种优化设计方法均取得了长足的发展，且都一定程度推广应用到了型号任务设计中。而在大型民机气动优化设计问题中，由于大规模设计变量及约束的引入和优化设计效率的需求，使得基于伴随理论的梯度优化设计方法的工程应用需求更加迫切。

国内外在气动伴随优化设计方面开展了很多的研究，国外有斯坦福大学Jameson团队[1]、多伦多大学Zingg团队[2-3]、密歇根大学Martins的多学科优化设计(MDO)团队[4-6]等研究成果。国内有西北工业大学的杨旭东[7-8]、高正红团队[9-10]、白俊强团队[11-12]等，南京航空航天大学的唐智礼[13]、夏健[14]、刘学强[15]，中国空气动力研究与发展中心的邓有奇[16]、吴文华[17]、黄江涛[18-19]等研究成果。上述研究在处理优化问题时，均集中在以多设计点、多目标、多约束、鲁棒优化等优化方法的基础上。而在工程实际应用设计中尚存在些许困难。

相比于伴随优化设计方法，压力分布反设计优化方法由来已久，从早期国外Tranen[20]在二维跨声速翼型上的应用，Henne[21]基于全速势方程在三维跨声速机翼上的应用，到20世纪90年代国内乔志德[22]在自然层流超临界翼型设计上的应用和朱自强等[23]在跨声速翼型和机翼上的应用。而近几年将反设计方法与遗传类优化方法相结合的研究较多，如西北工业大学宋文萍、韩忠华[24-25]等，但大多集中在层流翼型的设计中，对三维机翼反设计的研究不多。而将压力分布反设计方法与伴随方法相结合的研究相对较少，由于目标压力分布装配面临的困难特别是在非结构混合网格上，使得直接应用于复杂三维超临界机翼设计的更少。通常经过反设计的二维翼型装配到机翼上之后压力分布会发生变化甚至恶化，这种变化和恶化往往很难从二维反设计中去消除，所以开展高效的三维超临界机翼压力分布反设计，具有很大的工程意义。

压力分布特征作为翼型和机翼设计的核心，是翼型或者机翼最明显的特征，也是设计师设计思想的集中表现。而当前的气动优化设计方法常以简单气动力指标为设计目标，导致优化设计得到的构型在设计点处气动性能异常优越，偏离设计点的气动性能却恶化严重，这种优化设计结果，往往在工程应用中很难被采用。对于此类问题，已有研究者将压力分布特征以约束形式引入优化设计中，来提高这种自动寻优方法的工程应用性，且目前主要集中在遗传类的非梯度优化方法中，如李焦赞等[26]在翼型设计中的应用，中国航天空气动力技术研究院的杨昆淼等[27]在C919机翼上的应用，清华大学的陈海昕团队[28-29]提出的“人在回路”思想等，均取得了很好的工程应用效果。这些研究从一定程度上解决了优化设计自动化寻优和设计师工程设计思想之间的矛盾。但当前尚未有将压力分布约束引入伴随气动优化设计方法的工程类研究。而凭借伴随气动优化在高效、大规模设计变量方面的优势，开展此项研究对于工程实际应用具有非常大的意义。

本文研究直接应用于三维机翼的伴随压力分布反设计方法和压力分布约束方法，通过尝试与探索，提出了三维构型伴随压力分布反设计与伴随优化相结合的优化设计思路，增强伴随梯度类优化方法中引入设计思想的可行性，提升伴随气动优化方法在工程型号设计中的参与程度。

1 离散伴随方程和敏感性导数

考虑流场空间离散残值R为设计变量D、流动变量Q和计算网格X的函数，对定常流动有R(D,Q,X)=0。同时线弹性体网格变形方程G(D,X)=0。引入拉格朗日函数：

L(D,Q,X,Λf,Λg)=f(D,Q,X)+

(1)

式中：f为优化目标函数，目前可选为各种气动力系数及其组合；Λf和Λg均为拉格朗日算子。将式(1)拉格朗日函数对设计标量变分整理得

(2)

由于Λf和Λg可为任意值，可令

(3)

以避免∂Q/∂D和∂X/∂D的计算。相应地，式(2)变为

(4)

方程组式(3)中的两式分别称为离散流场伴随方程和离散网格伴随方程，其求解规模与设计变量个数无关。对其中各方程采用隐式时间推进并行求解，离散所得线性方程组采用GMRES方法求解。求得Λf和Λg，代入式(4)便可得到目标函数对设计变量的导数。可见，采用伴随方法求敏感度，在流场求解基础上还需求解一次流场伴随方程，其计算量与流场求解相当。

进行压力分布反设计之前，在优化定义及设置中需要给定期望的目标压力分布，以初始翼型(二维)或某一剖面(三维)的压力分布与目标压力分布之差的平方和作为目标函数，对于第j个状态：

(5)

图1 压力分布插值Fig.1 Pressure distribution interpolation

对于有k个截面目标的反设计，加权的目标函数为

(6)

2 翼面压力分布约束控制方法

在针对三维机翼进行考虑压力分布约束的伴随气动优化设计时，压力分布目标/约束需要通过以下步骤进行添加：

1) 确定机翼展向控制剖面站位信息(展向Y值坐标)。

2) 对各剖面的压力分布进行上下翼面分割，分开上下翼面压力分布，确定分割信息。

3) 分别给出上下翼面目标/约束压力分布信息(可分段给出X和Cp信息)。

对于如何分段给出目标/约束压力分布，以典型超临界压力分布为例。如图2所示，可将上翼面压力分布分为4段，用来控制超临界机翼的压力分布特性，具体控制方法如下：

1)X1～X2段用于控制吸力平台区域。

2)X2～X3段用于控制激波位置及强度。

3)X3～X4段用于控制波后加速区。

4)X4～X5段用于控制压力恢复区。

而下翼面的压力分布可以以同样的方式给出，也可按整体处理。

图2 压力分布约束Fig.2 Pressure distribution constraints

3 面向工程应用的伴随优化及压力反设计

针对自研非结构流场解算器(WiseCFD-UG)发展了离散伴随求解模块。WiseCFD-UG软件采用基于格点的二阶有限体积法数值求解雷诺平均Navier-Stokes(RANS)方程和SA湍流模型方程来实现工程复杂湍流模拟。空间离散基于包含四面体、六面体、三棱柱、金字塔等单元类型的混合网格，以利用混合网格对复杂外形的适应能力。无黏通量离散采用Roe格式，并采用非结构重构来实现二阶精度。黏性通量离散采用中心格式。梯度计算采用Green-Gauss方法。

图3 考虑压力分布约束的气动外形优化设计流程Fig.3 Framework of aerodynamic optimization design with pressure distribution constraints

图3为气动优化设计时提出的面向工程应用的离散伴随优化设计流程，采用气动优化设计与压力分布反设计相结合的设计思路。气动优化设计进行自动寻优，对寻优结果的压力分布进行筛选和修正，再采用伴随压力分布反设计对自动寻优构型进行反设计，或者压力分布作为约束被引入到此轮优化设计中，从而实现离散伴随优化方向的可控性。解决了自动优化寻优与工程适用性之间的矛盾，使得设计师的设计思想在自动优化寻优的过程中得到体现。

4 超临界机翼气动压力反设计算例

对于三维伴随气动反设计优化的验证，选取具有典型特性的超临界构型CRM(Common Research Model)机翼，此模型作为美国航空航天学会组织空气动力优化讨论会(Aerodynamic Design Optimization Discussion Group, ADODG)发布的优化案例模型，具有一定的研究价值。在此构型上定义两个优化算例，分别命名为Case 1和Case 2，其具体设计要求如表1所示。

4.1 Case 1压力分布反设计优化

以随机扰动变形后的CRM机翼为初始构型，而以CRM机翼原始构型在巡航设计状态下8个不同站位的压力分布作为目标压力分布，完成在马赫数Ma=0.85,雷诺数Re=5.0×106巡航状态下压力分布反设计。反设计优化所用的计算网格为Pointwise软件生成的非结构混合网格，网格节点数量约为200万。远场距离为参考弦长的100倍，网格法向增长率1.2，并保证物面第1层网格y+小于1，表面网格如图4所示。

采用三维自由变形(Free-Form Deformation, FFD)参数化方法对CRM机翼进行控制，FFD控制框如图5所示，通过11 个不同站位下FFD控制点对其进行控制，控制点分布为10×11×2，并以控制点作为设计变量，共计220个。

表1 Case 1和Case 2定义Table 1 Definition of Case 1 and Case 2

图4 CRM机翼表面网格Fig.4 Surface grid of CRM wing

图5 CRM机翼FFD控制体Fig.5 FFD control box of CRM wing

此算例的目标函数(Object_function)定义为8个不同站位剖面压力分布与目标压力分布的偏差和，即

式中：n为截取的第i个站位剖面上的网格节点数，而目标压力分布在网格节点上的装配方式如第1节所示。

图6给出了反设计优化目标函数的收敛曲线，目标函数值在20轮迭代后几乎不再变化，目标函数值不断向“0”逼近，并最后达到收敛标准。反设计经历59轮迭代，在96核并行集群上运行24 h，可以得到反设计结果，可以实现“一日设计循环”的工程设计效率要求。

表2为反设计前后机翼气动力参数的变化对比，其中α为迎角，CL为升力系数,CD为阻力系数,Cmy为力矩系数反设计构型与目标构型的气动力参数值吻合较好，阻力系数相差0.2 counts，基本可以忽略不计，单从气动力系数方面也可以看出压力分布反设计的效果很好。

图6 Case 1目标函数收敛历程Fig.6 Convergence history of objective function for Case 1 study

表2 Case 1反设计前后气动力性能对比Table 2 Comparison of aerodynamic performance of base- line and inverse design result in Case 1 study

巡航设计状态下初始的CRM随机自由扰动机翼(Baseline)、CRM原始构型机翼(Target)、CRM压力反设计机翼(Inverse design)的表面及剖面压力分布对比如图7所示。图中：y为机翼展向位置，b为机翼半展长，可以看出，伴随压力反设计优化能有效地将初始的CRM随机自由扰动机翼糟糕的压力分布反设计为带有超临界压力分布特征的目标压力分布，在经过反设计后与给定的目标压力分布具有较好一致性，且在不同剖面上压力分布的前缘峰值、吸力平台区、前后加载、波后加速区、激波强度和激波位置等特征方面都与目标压力分布一致，由于翼梢位置没有给定目标压力分布，所以翼梢处表现稍有差别。图8给出了反设计后6个典型控制剖面的几何对比，反设计后的典型剖面的翼型与目标翼型几何对应较好。图中：z为机翼法向位置。上述结果表明，三维伴随压力分布反设计方法能够实现具有超临界压力分布特性的宽体客机机翼气动优化反设计。

4.2 Case 2压力分布反设计优化

对于Case 2，将CRM机翼原始构型作为反设计的初始构型，而目标压力分布选取离散伴随三点优化得到的巡航设计状态(Ma=0.85)下8个 不同站位的压力分布，完成马赫数Ma=0.85，雷诺数Re=5.0×106巡航状态下的压力分布反设计。这里的离散伴随三点优化得到的压力分布是指：考虑阻力发散特性的三点优化，即在马赫数Ma=0.83、0.85、0.87三个设计点，权重因子分别为0.25、0.5、0.25下的加权减阻优化，详细的优化设置及描述可参考文献[30]。FFD控制框及控制点的分布与4.1节中Case 1相同。由于三点优化得到的压力分布是含迎角α作为设计变量的结果，所以此算例进行反设计时也将攻角作为设计变量，设计变量个数共计221个。而目标函数(Object_function)的定义与4.1节Case1的定义方式相同。

图9给出了Case 2的收敛历史曲线，反设计经历48轮，同样在96核并行集群上运行19 h可达到收敛标准，得到反设计构型，同样满足工程设计效率要求。

表3给出了反设计构型与目标构型气动力系数的对比，可以看出经过压力分布反设计之后的反设计构型，升阻力系数及俯仰力矩系数值与目标构型值吻合较好，迎角α也增加到目标构型对应值接近的范围。从气动力系数的比较上，进一步验证了此算例压力分布反设计的效果。

图8 Case 1初始构型、目标构型与反设计构型6个典型站位翼型形状对比Fig.8 Comparision of sectional airfoil shapes for baseline, target and inverse design configuration at 6 span-wise positions in Case 1 study

图9 Case 2目标函数收敛历程Fig.9 Convergence history of objective function for Case 2 study

表3 Case 2反设计前后气动力性能对比Table 3 Comparison of aerodynamic performance of base- line and inverse design result in Case 2 study

图10给出巡航设计状态下初始的CRM原始构型机翼(Baseline)、CRM多点优化后目标构型机翼(Target)、CRM压力反设计机翼(Inverse design)的表面及剖面压力分布对比。可以看出反设计构型的压力分布与给定的经多点优化后的目标压力分布云图吻合较好，对比典型站位处剖面压力分布也与目标压力分布一致。图11给出了反设计后6个典型剖面的几何对比，反设计后的典型剖面的翼型与目标翼型几何对应较好。上述结果表明，本文的伴随压力分布反设计方法能够直接针对大型客机三维超临界机翼进行气动优化设计。

图10 Case 2反设计机翼的表面压力云图及典型站位压力分布Fig.10 Pressure contour and pressure distribution of typical span-wise stations of inverse design wing in Case 2 study

图11 Case 2初始构型、目标构型与反设计构型在6个典型站位翼型形状对比Fig.11 Comparision of sectional airfoil shapes for baseline, target and inverse design configuration at 6 span-wise positions in Case 2 study

5 宽体客机压力分布修正优化设计

考虑到在真实的工程设计应用中，直接给出满足设计要求的目标压力分布具有很大的困难，为此本节将伴随压力分布反设计方法中对压力分布的处理方式以压力分布约束的形式引入到宽体客机的伴随优化设计中，即将伴随压力分布反设计中的优化目标转换为伴随优化中的优化约束。针对某宽体客机，开展考虑压力分布修正的伴随多点优化设计。首先对宽体客机进行考虑阻力发散特性的两点优化(Ma=0.85和Ma=0.87)，并定义优化算例为Design 1。对Design 1在巡航设计点的优化压力分布进行修正，并结合工程设计要求，将修正后的压力分布引入优化设计系统中，对Design 1的结果进行二次优化设计。在二次设计中，主要对主设计点马赫数Ma=0.85进行优化设计，将压力分布作为约束，目标函数以及其余的几何、气动约束与多点优化一致，巡航攻角不再作为设计变量，定义此轮优化算例为Design 2。

5.1 Design 1不考虑压力分布约束的两点优化

采用常规离散伴随多点优化策略，针对宽体客机全机构型，开展考虑阻力发散特性的两点优化设计。具体的优化问题定义如表4所示，FFD控制体分布及约束布置如图12所示。

两点优化前后的巡航状态下的压力分布比较见图13。图中：虚线为原始机翼压力分布曲线，实线为两点优化后压力分布曲线。在巡航状态下，优化后机翼的各剖面压力分布激波强度减弱，且维持着超临界压力分布特性，激波位置不同程度的前移，这对阻力发散特性是有益的。但激波后存在不同程度的二次加速区，且内翼压力分布存在抖动，趋于无激波形态，机翼压力分布等压线沿展向的一致性变差，这往往反应在阻力蠕增特征曲线上，低马赫数范围的阻力蠕增特性变差，这对工程应用带来挑战。

表4 宽体客机两点气动优化问题说明Table 4 Wide-body aircraft multipoint aerodynamic shape optimization problem statement

图12 宽体客机FFD控制体分布及约束分布Fig.12 Constraints location and FFD control box of wide-body aircraft

5.2 Design 2考虑压力分布约束的优化

在5.1节Design 1的优化结果的基础上，选取优化构型的巡航状态为优化对象。对压力分布适当修正，并以约束的形式添加到第2轮优化Design 2中，具体的压力分布修正准则如下：

1) 展向直等压线分布控制

2) 吸力平台压力波动修正

对10个控制剖面上翼面10%～50%弦长位置进行压力系数采样，采样的压力系数采用最小二乘法进行线性拟合，得到

式中:xij为第i剖面上第j个采样点的x位置信息，定义采样点与线性拟合值的偏差为吸力平台约束函数：

3) 激波后二次加速区修正

将修正后的压力分布以压力分布约束的形式添加到Design 1优化得到的巡航点设计结果中，再进行一轮伴随单点优化设计。图14给出了初始构型及Design 1和Design 2的巡航马赫数下的优化压力分布对比。可以看出考虑压力分布约束优化后的Design 2构型，展向直等压线分布得到改善，内翼吸力平台波动得到改善，波后的二次加速区得到一定的抑制。对初始构型以及Design 1和Design 2的构型进行阻力发散特性评估，评估结果如图15所示。从图中可以看出：Design 1在高马赫数范围内获得了不错的机翼阻力收益，但低马赫数范围的阻力蠕增特性变差，而考虑压力分布约束的Design 2优化结果，使得在马赫数小于巡航设计马赫数时出现阻力蠕增现象得到改善，巡航设计点获得2.9 counts的收益。

图13 宽体客机机翼第1轮优化后表面压力云图及典型站位压力分布Fig.13 Pressure contour and pressure distribution of typical span-wise stations of optimal wing in the first round optimization design

图14 宽体客机机翼第2轮优化后表面压力云图及典型站位压力分布Fig.14 Pressure contour and pressure distribution of typical span-wise stations of optimal wing in the second round optimization design

由此说明，应用的伴随压力分布反设计方法在宽体客机超临界机翼设计中具有一定的工程实用性。为进一步开展宽体客机气动优化设计方法研究打下一定的基础。

图15 阻力发散曲线对比Fig.15 Comparison of drag divergence curves

6 结 论

1) 采用目标点邻域线性插值的方法解决了目标压力分布在非结构混合网格剖面上的装配难题。

2) 直接对三维CRM机翼进行压力分布反设计，验证本文的伴随压力分布反设计方法能有效的对宽体客机机翼进行反设计优化，反设计效果明显，效率高，结果可靠。

3) 利用伴随压力分布反设计方法直接对三维机翼进行反设计，相比于传统翼型反设计再装配方案，具有更高的效率和精度。

4) 在宽体全机构型离散伴随两点气动优化结果的基础上，引入压力分布反设计约束方法，改进了展向直等压线分布和吸力平台区压力分布波动等问题，使优化构型更加接近工程实用性。

5) 对于宽体客机超临界机翼设计，采用考虑压力分布约束的伴随多点优化是很有必要的。这种方法将数值优化与人工经验结合，可以促进伴随气动优化方法在工程应用中的实用性。