刘 俊，宋文萍，韩忠华，许建华，樊艳红

（西北工业大学 翼型叶栅空气动力学国防科技重点实验室，陕西 西安 710072）

0 引 言

传统的基于CFD数值模拟的气动优化设计方法主要分为梯度法和非梯度法。非梯度方法使用遗传算法（GA）等直接调用流场分析程序进行优化，它们往往具有良好的全局性，但计算量较大。基于梯度的方法的关键在于梯度信息的计算，传统的有限差分法虽然简单，但在设计变量较多时，计算量较大。

从20世纪80年代末开始，由Jameson［1］等人发展起来的基于控制理论的优化方法（Adjoint方法）以其对多变量问题的适用性和高效的特点迅速成为气动设计的主流方法，并已大量应用于气动优化设计和反设计中［2－4］。直至当今，该方法仍然在气动设计中占据重要地位［5－6］。

虽然Adjoint方法具有高效的特点，但其仍然属于局部优化方法的范畴，且该方法通用性不强，针对不同类型的优化问题，需要重新推导伴随方程，重新编制伴随方程求解程序；在特殊情况下，如考虑自由转捩的NS方程求解器，Adjoint方法应用起来有难度。为了兼顾全局性与高效性，同时考虑到通用性的要求，基于代理模型的优化方法逐步流行起来，并且越来越受到重视。在众多的代理模型优化方法中，基于多项式响应面的方法（P－RSM）最为成熟，并已成功应用于气动优化设计与反设计中［7－8］。

近十年来，由于Kriging模型具有可以模拟复杂响应并可提供误差信息的优势，基于Kriging模型的优化方法［9］被受到重视并越来越多地应用于气动及多学科优化设计中［10－17］，以大大减少优化设计所需的计算时间并得到优良的设计结果。

尽管基于Kriging模型的优化方法在气动优化设计中取得了很大的成功，但一直未在反设计中得到应用。虽然Toal等［18－19］以翼型反设计为算例，比较了使用不同优化算法来优化Kriging超参数对最终反设计结果的影响，然而他们的出发点并非基于Kriging模型的优化方法在翼型反设计中的应用，也未得到真正有效的反设计结果。其原因可能是使用了不恰当的翼型参数化方法，且使用的样本点加点方法不当。可能正是受他们的影响，使得基于Kriging模型的优化方法未能应用于反设计中。因此，研究Kriging模型在翼型反设计中的可行性，发展有效、实用的基于Kriging模型的翼型反设计方法是十分必要的。

本文利用 Kriging模型，结合EI方法［9］及直接优化代理模型最小值等样本点加点准则、优化算法、CST参数化方法，进行了翼型单目标、多目标反设计。本文方法中，按照模型加点准则又分为三种：单独使用EI方法（SIC1）、单独使用直接优化代理模型最小值（SIC2）、同时使用上述两种加点准则（SIC3）。

首先在本文优化方法中分别使用三种不同加点方法进行了翼型单目标反设计，验证了本文方法的可行性并比较了不同加点方法对反设计结果的影响，同时与传统的基于二次响应面模型的方法及Adjoint方法进行了比较，结果表明，分别使用三种加点方法的基于Kriging的方法均适用于翼型反设计，且都优于二次响应面方法，且使用SIC2效果最好；与Adjoint方法相比，反设计效果及计算效率相当。

其次，在本文优化方法中应用SIC2加点方法进行了翼型多目标反设计，验证了本文方法在翼型多目标反设计中的适用性。

最后，将本文方法应用于接近工程实际的在已有翼型基础上修改压力分布进行了单目标、对目标反设计，验证了本文方法的可行性。

1 Kriging模型

Kriging模型是由南非地质学家Krige［20］提出来的一种统计学插值模型，并由Sacks［21］最早应用于计算机实验设计与分析，随后作为一种代理模型被广泛用于数值模拟分析与优化中响应值的预测。本文采用的是普通Kriging模型。

1.1 Kriging预测值与均方误差

Kriging模型用一个常数和一个随机过程的和来表示系统的响应值：

随机过程Z（·）的均值为0，协方差如下：

其中σ2是Z（·）的方差（假定对任意的x有（x）＝σ2（x）≡σ2），R是只与两点x和x′的距离有关的相关函数。

假定设计空间内任意一点的响应值可由已知样本点的响应值ys的线性组合来近似，则在未知点x，Kriging近似值有如下形式：

其中，w＝（w（1），...，w（ns））T为权系数。

用YS＝（Y（1），...，Y（ns））T替换上式中的yS＝（y（1），...，y（ns））T，最小化均方误 差 （Mean Squared Error）可得到Kriging的预测值如下：

其中，1是单位列向量，

任意未知点x处的Kriging预测值的均方误差（MSE）为：

其中，是方差σ2的估计值：

1.2 相关模型

相关矩阵R和相关向量r的建立与相关函数的选取有关，假定相关函数的函数值只与两点x（i）、x（j）之间的空间距离有关。相关模型的选取应该具有以下特征：a）当距离趋于0时，函数值趋于1；b）当距离增加时，函数值光滑的减小；c）当距离趋于无穷时，函数值趋于0；d）至少一阶可导。此处只考虑以下形式的相关模型：

本文采用三次样条函数作为相关函数：

其中

1.3 Kriging参数的获得

式（9）中的 Kriging超参数θ（hyper parameter）可以通过求解下列最大似然估计（MLE）这一优化问题获得：

2 加点准则

尽管Kriging模型具有较好的全局近似能力，但是仅依靠初始样本信息并不能找到设计空间内的最优值，因而需要通过增加新的样本点来寻找真实最优值并提高Kriging模型精度。本文采用了以下两种加点准则［22－23］。

2.1 最大化Expected Improvement（EI）

设样本集中的最优目标函数值为ymin，假设随机变量Y（x）服从均值为（x），标准差为s（x）的正太分布，Y（x）∈N［（x），s2］，其 中（x）为 Kriging 预 测值，见式（4），s2为均方误差（MSE），见式（7）。在x处目标函数相对的改进量为I＝ymin－Y（x）＞0，于是I的期望值由下式给出：

求出设计空间内EI最大的点作为新的样本点。

2.2 最小化代理模型 （MP）

用优化算法寻找代理模型上目标函数的最小值，将最小值点作为新的样本点加入样本集。

3 二次响应面模型

对于完全二阶多项式模型的一般公式为：

其中k为设计变量个数，c0、ci、cij是待定参数。

对于k元完全二阶多项式，待定参数的个数为：

利用样本点的响应值可求得所有待定参数。

4 流场求解

本文采用课题组自主开发的PMNS2D程序通过求解雷诺平均Navier－Stockes方程对翼型进行气动分析，采用中心格式离散，湍流模型采用Spalart－Allmaras模型，并使用了隐式残值光顺、当地时间步长、多重网格等加速收敛措施。计算网格为C型结构化网格。

5 几何参数化

5.1 CST参数化方法

CST参数化方法是由波音公司的工程师Kulfan［24］提出来的。该方法用一个类函数C（x／c）和一个型函数S（x／c）的乘积加上一个描述后缘厚度的函数表示一个翼型的几何形状：

其中c为翼型弦长，类函数C（x／c）具有以下的一般形式：

型函数S（x／c）由伯恩斯坦多项式各项的加权求和得到：其中Ai是加权系数，Kr，n为多项式的系数：

5.2 设计变量范围的确定

利用一个较“厚”的翼型和一个较“薄”的翼型，并使得这两个翼型之间的范围足够大，可以包含目标翼型，将两者之间作为真实的设计空间，如图1中阴影部分所示；将这两个翼型所对应的CST参数作为设计变量的上下限。

图1 翼型反设计的设计空间示意图Fig.1 Schematic diagram of the design space for airfoil inverse design

6 翼型反设计目标函数

以设计结果的压力分布和目标压力分布的差为目标函数，对于第k个设计状态：

其中，Cpi为计算翼型表面第i个点的压力系数，＾Cpi为目标翼型第i个点的压力系数。对于m个目标的反设计，加权目标函数为：

ωk为第k个目标的权重系数。对于更重要的目标压力分布，其权系数应取得更大。

7 代理模型优化方法

本文以当前翼型的压力分布和目标压力分布的差为目标函数，将反设计问题转化为优化问题，用基于代理模型的优化方法来求解。

基于代理模型的优化方法首先采用试验设计方法（Design of Experiment，DoE）生成设计空间内的样本点，并获得样本点的响应值，由样本信息建立代理模型，以近似描述空间内任意位置的真实响应值，采用优化算法（本文采用遗传算法）寻找新的样本点，将新的样本点加入原样本集，重新建立代理模型进行优化直至满足停止条件。

1）基于二次响应面模型的方法（P－RSM）

采用二次响应面模型，直接寻找二次响应面目标函数的最小值点作为新的样本点。

2）Kriging—SIC1

采用Kriging模型，寻找设计空间内EI最大值点作为新的样本点，即EI方法。

3）Kriging—SIC2

采用Kriging模型，直接寻找设计空间内目标函数的最小值点作为新的样本点，即MP方法。

4）Kriging—SIC3

采用Kriging模型，同时使用EI、MP两种加点方法，一次增加两个样本点。

流程图如图2所示。

图2 代理模型优化方法流程图Fig.2 Flowchart of the surrogate－based optimizer

8 方法验证

8.1 本文方法在反设计中的适应性研究以及与PRSM、Adjoint方法的比较

以跨声速翼型RAE2822单目标反设计为例进行测试。设计状态如下：

使用19个设计变量，使用代理模型方法进行反设计时，取相同的设计空间。

首先采用P－RSM方法进行反设计，反设计的压力分布及翼型几何形状与目标压力分布及目标翼型的比较见图3和图4，从图中看出，虽然压力分布及翼型几何形状与目标翼型的比较接近，但存在明显差距。

其次，用Adjoint方法进行反设计，初始翼型为NACA0012。反设计结果见图5、图6，从图中看出，反设计的压力分布及翼型几何形状与目标压力分布及目标翼型均吻合良好。

图3 P－RSM方法反设计压力分布与目标压力分布比较Fig.3 Comparison of the designed pressure distribution by P－RSM and the target one

图4 P－RSM方法反设计翼型与目标翼型比较Fig.4 Comparison of designed airfoil by P－RSM and the target one

图5 Adjoint方法反设计压力分布与目标压力分布比较Fig.5 Comparison of the designed pressure distribution by Adjoint and the target one

图6 Adjoin方法反设计翼型与目标翼型比较Fig.6 Comparison of designed airfoil by Adjoint and the target one

再次，使用本文基于Kriging模型的方法分别采用三种加点方法进行反设计。反设计的压力分布及翼型几何形状与目标压力分布及目标翼型的比较见图7～图12，从图可以看出，分别使用三种加点方法的基于Kriging模型的方法反设计的翼型压力分布与几何形状与目标翼型的均吻合良好，且均优于PRSM方法，从而验证了本文基于Kriging的方法在翼型反设计中的适应性。此外，从三种加点方法看，使用SIC2即仅采用直接优化代理模型最小值的加点方法所得的反设计效果最好。

图7 采用SIC1的Kriging方法反设计压力分布与目标压力分布比较Fig.7 Comparison of the designed pressure distribution by Kriging－SIC1and the target one

图8 采用SIC1的Kriging方法反设计翼型与目标翼型比较Fig.8 Comparison of the designed airfoil by Kriging－SIC1and the target one

图9 采用SIC2的Kriging方法反设计压力分布与目标压力分布比较Fig.9 Comparison of the designed pressure distribution by Kriging－SIC2and the target one

图10 采用SIC2的Kriging方法反设计翼型与目标翼型比较Fig.10 Comparison of the designed airfoil by Kriging－SIC2and the target one

图11 采用SIC3的Kriging方法反设计压力分布与目标压力分布比较Fig.11 Comparison of the designed pressure distribution by Kriging－SIC3and the target one

图12 采用SIC3的Kriging方法反设计翼型与目标翼型比较Fig.12 Comparison of the designed airfoil by Kriging－SIC3and the target one

表1给出了本文方法及P－RSM方法最终的目标值的比较。注意此表中的目标值是实际目标值与Kriging初始样本点中最大目标值的比值。从表中也可以看出，采用SIC2所得的目标值最优，SIC3次之，SIC1再次，但都优于P－RSM方法。

表1 四种方法目标值与CFD计算次数比较Table 1 Comparison of objective value and runs of flow solver for the 4methods

图13给出了本文方法及Adjoint方法反设计目标函数的收敛曲线，其中前面的方形符号表示建立Kriging模型初始样本点的目标值。注意，此图中的使用Kriging模型的目标值是真实目标值与初始样本点中最大目标值的比值；而Adjoint的目标值是实际目标值与初始目标值的比值，且把求解一次伴随方程的计算量视为一次流场求解的计算量。从该图可以看出，本文方法中采用SIC2的目标函数值下降最快，且最终结果最优。此外，使用SIC2的基于Kriging的方法与Adjoint方法相比，收敛速度相当且目标值下降量级相当。

图13 使用3种加点准则的Kriging方法及Adjoint方法反设计目标值收敛曲线比较（对于Adjoint方法，一次Adjoint方程求解视为一次流场求解）Fig.13 Comparison of convergence histories for the 3SIC and the Adjoint method（For Adjoint method，the cost of Adjoint solver is treated as equal to the cost of flowsolver）

8.2 翼型多目标反设计

为了验证本文方法在翼型多目标反设计中的适应性，以RAE2822翼型两个目标的反设计为例进行测试，即设计一个翼型同时满足两个状态的压力分布。设计状态如下：

设计变量及其范围与单目标算例相同。由8.1可知，使用SIC2所得的反设计效果最优，故此处采用SIC2加点方法。作为算例，此处第一个状态的权重系数取0.6。

反设计翼型的压力分布及几何形状见图14、图15，从图看出，压力分布与及几何形状与目标压力分布及几何形状都吻合得很好。图16给出了目标值的收敛历程，其中纵坐标为真实目标值与初始样本点中最大目标值的比值。

9 工程应用的可行性验证

反设计方法是工程实际中进行翼型设计的主要方法之一。翼型反设计一般首先在已有翼型基础上修改其压力分布，再进行反设计，以获得更优性能的翼型。为了验证本文方法在工程应用中的可行性，本节先对已有翼型的压力分布进行修改，再进行反设计。

图14 基于Kriging的方法反设计压力分布与目标压力分布比较Fig.14 Comparison the designed pressure distributions and the target ones

图15 基于Kriging的方法反设计翼型与目标翼型比较Fig.15 Comparison of the designed airfoil by Kriging－SIC2and the target one

图16 翼型多目标反设计目标值收敛历程Fig.16 Convergence history of the multiobjective inverse design

9.1 单目标反设计

以RAE2822为初始翼型，其在状态Ma＝0.73，Re＝6.5e6，α＝2.3°下，翼型上表面中部存在强激波，故其阻力较大，此处将上表面压力分布加以修改，以期在同样的流动状态下无激波，从而大大降低阻力。在反设计过程中，只对翼型上表面进行参数化，下表面外形保持不变。在修改上表面压力分布时，前缘附近（x／c＜0.1）通过手动调节，之后通过两段直线连接，且使其具有较大范围的超声速区并保证升力系数不变。

反设计所得压力分布、修改后的压力分布（即目标压力分布）以及原压力分布的对比如图17所示。从该图看出，反设计压力分布与目标压力分布吻合地很好，这是由于此处给出的目标压力分布比较合理；若目标压力分布不合理，则可能无法获得与之吻合很好的翼型，同时也说明当目标压力分布合理时，本文方法能获得满足目标压力分布的翼型。

图17 设计压力分布、修改后压力分布及原始压力分布比较Fig.17 Comparison of the designed、modified and initial pressure distributions

9.2 翼型多点反设计

仍以RAE2822为初始翼型，取两个设计状态

其中，第一个状态与上例中的相同，并将其定为主设计点。第二个状态为非设计点，希望在保证第一设计状态性能的同时降低该状态的阻力，从而改善翼型的阻力发散特性。与上例相同，仅对翼型上表面进行参数化，第一个状态的目标压力分布同上例。在修改第二状态的压力分布时，保证了升力系数不变，并假设翼型中部存在弱激波。

此处使用两组加权系数（0.6，0.4）、（0.4，0.6）分别进行反设计。采用两组加权系数进行反设计所得的压力分布与目标压力分布及原压力分布的对比分别见图18和图19。可以看出，两点设计所得的压力分布与目标压力分布差别较大，这是由于给出的目标压力分布并不一定合理，同时满足两个目标压力分布的翼型并不存在。此外，对比图18和图19，当第一个状态的权重系数减小时，设计的翼型在第一个状态的压力分布与目标压力分布吻合度降低，而在第二个状态的压力分布与目标压力分布吻合度提高。

图18 权重系数取（0.6，0.4）时，设计压力分布、修改后压力分布及原始压力分布比较Fig.18 Comparison of the designed、modified and initial pressure distributions

图19 权重系数取（0.4，0.6）时，设计压力分布、修改后压力分布及原始压力分布比较Fig.19 Comparison of the designed、modified and initial pressure distributions

图20给出了原始翼型、单目标设计、多目标设计翼型在迎角2.3°下的阻力系数随马赫数变化曲线，从该图看出，单点设计的阻力特性明显优于原翼型，而多点设计尽管没能很好吻合目标压力分布，但其阻力特性仍优于单点设计。

图20 阻力系数随马赫数变化对比Fig.20 Comparison of the drag coefficients at different Mach numbers

从前面两个算例看出，进行翼型反设计，首先需要给出合理的目标压力分布，尤其是多目标反设计。实际应用中，往往需要在反设计与修改压力分布之间进行多次迭代。

10 结论与分析

（1）在基于Kriging模型的气动优化设计中，通常使用EI方法作为寻优及提高模型精度的准则。而当应用于气动反设计时，应当使用MP加点准则，即直接优化目标函数最小值的方法效果更优。

（2）在进行翼型反设计时，Kriging模型明显优于多项式响应面模型。

（3）本文方法与Adjoint方法相比所需流场分析次数相当，但是其通用性好，可以方便的更换流场分析程序。

（4）Kriging模型不仅适用于气动优化设计，同时也适用于气动反设计。而合理的目标压力分布是反设计的关键，实际应用中往往需要在反设计与修改压力分布之间进行多次迭代。

致谢：感谢王乐提供Adjoint翼型反设计的结果。

［1］JAMESON A.Aerodynamic design via control theory［J］.Journal ofScientificComputation，1988，3（3）：233－260.

［2］JAMESON A.Control theory based airfoil design using the Euler equations［R］.AIAA Paper 94－4272.

［3］JAMESON A，PIERCE N A，MARTINELLI L.Optimum aerodynamic design using the Navier－Stockes equations［R］.AIAA Paper 97－0101.

［4］YANG X D，QIAO Z D，ZHU B.Control theory based aerodynamic shape inverse design for subsonic and transonic wings［J］.ACTAAerodynamicaSinica，2003，21（1）：12－19.（in Chinese）杨旭东，乔志德，朱兵.亚、跨音速三维机翼气动外形反设计的控制理论方法［J］.空气动力学学报，2003，21（1）：12－19.

［5］WALTHER B，NADARAJAH S.Constrained adjoint－based aerodynamic shape optimization in a multistage turbomachinery environment［R］.AIAA Paper 2012－0062.

［6］WINTZER M，KROO I.Optimization and adjoint－based CFD for the conceptual design of low sonic boom aircraft［R］.AIAA Paper 2012－0963.

［7］XIONG J T，QIAO Z D，HAN Z H.Response surface based aerodynamic design of transonic wings［J］.ACTAAeronaotica etAstronauticaSinica，2006，27（3）：399－402.（in Chinese）熊俊涛，乔志德，韩忠华.基于响应面法的跨声速机翼气动优化设计［J］.航空学报，2006，27（3）：399－402.

［8］DENG L，QIAO Z D，XIONG J T.et al.A multi－objective inverse design method for natural laminar airfoil［J］.ACTA AeronaoticaetAstronauticaSinica，2010，31（7）：1373－1378.（in Chinese）邓磊，乔志德，熊俊涛，等.多目标自然层流翼型反设计方法［J］.航空学报，2010，31（7）：1373－1378.

［9］JONES D R，SCHONLAU M，WELCH W J.Efficient global optimization of expensive black－box functions［J］.Journalof GlobalOptimization，1998，13：455－492.

［10］XU R F，SONG W P，HAN Z H.Study on the improving Kriging－based airfoil optimization design method［J］.Journalof NorthwesternPolytechnicalUniversity，2010，28（4）：503－510.（in Chinese）许瑞飞，宋文萍，韩忠华.改进Kriging模型在翼型气动优化设计中的应用研究［J］.西北工业大学学报，2010，28（4）：503－510.

［11］SU W，BAI J Q.A surrogate－based optimization algorithm and its application to aerodynamic optimization design［J］.Journal ofProjectiles，Rockets，MissilesandGuidance，2008，28（3）：199－202.（in Chinese）苏伟，白俊强.一种代理模型方法及其在气动优化设计中的应用［J］.弹箭与制导学报，2008，28（3）：199－202.

［12］WANG X F，XI G.Aerodynamic optimization design for airfoil based on Kriging model［J］.ACTAAeronaoticaetAstronautica Sinica，2005，26（5）：545－549.（in Chinese）王晓锋，席光.基于Kriging模型的翼型气动性能优化设计［J］.航空学报，2005，26（5）：545－549.

［13］JEONG S，MURAYAMA M，YAMAMOTO K..Efficient optimization design method using Kriging model［J］.Journalof Aircraft，2005，42（2）：413－420.

［14］GOTO Y，JEONG S，OBAYASHI S，et al.Design space exploration of supersonic formation flying focusing on drag minimization［J］.JournalofAircraft，2008，45（2）：430－439.

［15］KANAZAKI M，TANAKA K，JEONG S，et al.Multi－objective aerodynamic exploration of elements＇setting for high－lift airfoil using Kriging model［J］.JournalofAircraft，2007，44（3）：858－864.

［16］HOYLE N，BRESSLOFF N W，KEANE AJ.Design optimization of a two－dimensional subsonic engine air intake［J］.AIAA Journal，2006，44（11）：2672－2681.

［17］SONG W，KEANE A J.Surrogate－based aerodynamic shape optimization of a civil aircraft engine nacelle［J］.AIAAJour－nal，2007，45（10）：2565－2574.

［18］TOAL D J J，BRESSLOFF N W，KEANE A J.Kriging hyperparameter tuning strategy［J］.AIAAJournal，2008，46（5）：1240－1252.

［19］TOAL D J J，BRESSLOFF N W，KEANE A J，et al.The development of a hybridized particle swarm for Kriging hyperparameter tuning［J］.EngineeringOptimization，2011，43（6）：1－28.

［20］KRIGE D G.A statistical approach to some basic mine valuations problems on the witwatersrand［J］.JournaloftheChemical，MetallurgicalandMiningEngineeringSocietyofSouth Africa，1951，52（6）：119－139.

［21］SACKS J，WELCH W J，MITCHELL，et al.Design and analysis of computer experiments［J］.StatisticalScience，1989，4：409－423.

［22］JONES D R.A taxonomy of global optimization methods based on response surfaces［J］.JournalofGlobalOptimization，2001，21：345－383.

［23］LIU J，HAN Z H，SONG W P.Comparison of infill sampling criteria in Kriging－based aerodynamic optimization［A］.28thInternational Congress of the Aeronautical Sciences，Brisbane，Australia［C］.2012，9.

［24］KULFAN B M.Universal parametric geometry representation method［J］.JournalofAircraft，2008，45（1）：142－158.