基于NACA 翼型参数化方法的气动优化设计

罗福星

（国网西藏电力有限公司超高压分公司，西藏 拉萨 850000）

0 引言

在飞行器机翼的设计过程中，高升阻比、低阻力等多方位的需求与约束，使得优化设计需要对机翼气动外形进行多变量、大范围迭代。在此过程中，机翼几何外形的参数化表达尤为重要，优秀的参数化方法既可以满足光滑的几何外形，增加设计的鲁棒性，又可以用较少的设计变量描述关键的翼型特性。常用的参数化方式有[1]：NACA 翼型、GA（W）翼型等由严格数学表达式给出翼型定义方式，由基准翼型和扰动形函数线性叠加的形函数扰动法[2]，由Sobiesky[3]提出的特征参数描述法，利用正交基函数描述外形的正交基函数法[4]。此外，还有较为常用的CST 参数化方法，该方法由波音公司的Kulfan 等[5，6]提出，适应性强。

而在气动计算中，得益于发展成熟的算力和CFD 算法较好的鲁棒性，RANS 计算方法在科研及工程过程中扮演着相当重要的角色，相较于采用直接数值模拟（DNS）和大涡模拟（LES），其占用计算资源大幅减少。

采用较为成熟的NACA 四位数翼型参数化方式，翼型最大相对厚度、最大相对弯度、最大弯度位置的3 个参数作为设计变量，用以描述翼型外形。在气动计算中采用RANS 计算方法，并结合SST-kw 湍流模型开展计算。

优化流程中采用拉丁超立方法进行实验设计（DOE），得到优化设计的响应面代理模型，在此基础上依据二次拉格朗日非线性规划（NLPQL）开展气动优化，由以上算法构筑了一套机翼优化设计方案。

1 翼型流场分析

1.1 NACA 翼型参数化方法

翼型参数化方法的选择对气动优化设计中的结果优劣及优化效率有相当关键的作用。采用典型的低速翼型——NACA 四位数翼型，该系列翼型为美国国家航空咨询委员会提出，其应用广泛，在翼型发展史中有着重要地位。其为经典的低速翼型，满足低速飞机的较大升阻比、较低的最小阻力系数且低阻范围宽以及失速过程缓和的要求。

NACA 四位数翼型的中弧线、厚度沿弦向的分布均由严格的数学表达式给出，翼型参数描述如图1所示[7]。

图1 翼型相关参数示意

中弧线表达式如下：

式中，h和hmax表示弯度及最大弯度，x为弦向位置，p为最大弯度位置，c为弦长。

厚度表达式如下：

式中t和tmax表示厚度和最大厚度。

上弧线和下弧线表达式为：

式中，（xu，yu）为上弧线坐标，（xd，yd）为下弧线坐标，θ为中弧线对应切线角度。

将中弧线表达式（1）及厚度表达式（2）代入（3）式中即可得到NACA 四位数翼型的参数表达式。

1.2 基于RANS 的数值模拟

采用基于求解雷诺平均的N-S 方程（RANS）进行流场的数值模拟，其利用Reynolds 时均应力简化连续方程和动量方程。

RANS 所求解的连续方程和动量方程如下：

1.3 标准算例验证

为验证采用的气动数值计算方法的准确性，采用马赫数为0.63、迎角为2 度工况下的NACA0012 翼型试验数据[8]作为验证依据。

流体域选择流向往前15 倍、往后25 倍的特征长度作为计算域范围，以减小远场边界的影响。采用RANS 求解，考虑计算工况选择SST-kw 湍流模型，详细计算格式见表1，网格及边界条件设置见表2。

表1 求解格式选择

表2 计算域网格设置

网格采用的是C-H 型拓扑构建的结构化网格，以使得网格离散方式更加贴合流场梯度的变化方向，且结构化的网格节点编码方式更有利于提高求解效率。网格具体如图2～3 所示。

图2 全计算域网格示意

图3 翼型表面边界层网格

用表2 的三套网格计算，结果见表3。表3 结果表明满足网格无关性要求，为此后续计算选择中等密度网格。获取该算例的翼型表面压力系数沿弦向的分布，如图4 所示，可以看到CFD 的计算结果与试验数据相差不大，表明采用的数值算法可靠。

表3 网格无关性验证

图4 翼型表面压力系数对比

2 优化设计方法

采用基于代理模型的优化设计流程框架，以翼型的气动特性参数作为目标函数开展优化分析，流程框架如图5 所示。

图5 翼型优化设计流程

2.1 响应面代理模型

在DOE 中通过计算足够数量的采样点构建样本空间，以此作为响应面模型的输入。样本采用拉丁超立方设计，该方法原理为：在m维的设计空间中，将每一维的坐标空间均分为n个子区间，在随机取n个点时保证每个子区间只被研究一次，从而构成样本数为n的拉丁超立方采样设计。在完成样本空间的目标函数计算后，可以据此构建代理模型。

响应面模型则是一种应用广泛、高鲁棒性的代理模型，采用多项式函数来拟合设计空间，反应设计变量和目标函数的关系[9]。采用二项式多项式作为响应面方程，其表达式如下：

此外，由于代理模型为近似模型，故需要取样本空间中的部分设计点用于交叉验证，进行代理模型的精度检测。

2.2 NLPQL 优化算法

在得到代理模型后，采用非线性序列二次规划法（NLPQL）求解其最优解。NLPQL 是应用二级泰勒级数来展开目标函数并将约束条件线性化，以此将非线性问题转化为了二次规划问题。该算法采用类牛顿矩阵Bk 定义Lagrange 函数的逼近和Hessian 矩阵的逼近，以此得到一个NLPQL 的子问题，表达式如下[10]：

式中：d为搜索方向，Bk为类牛顿矩阵，f（x）为目标函数，g（x）为约束函数，q为变量，xu和xl分别表示边界约束的上下限。

3 翼型的优化设计

3.1 二维翼型减阻优化

选择NACA0012 作为优化的基准翼型，设计工况为：马赫数取0.73，迎角取2.5°，雷诺数取6.5×106。

根据NACA 翼型参数化表达方式，以翼型最大相对厚度、最大相对弯度、最大弯度位置作为优化的设计变量；而优化设计的目标函数为设计气动工况下的阻力系数。此外，再考虑优化问题的约束条件：

（1）几何约束：最大相对厚度不小于原翼型；

（2）气动约束：升力系数不低于原始翼型。

最终得到如（8）所示的优化问题表达式。此外，补充设计变量的限制范围，控制优化翼型与原始翼型间的差异性，见表4。

表4 设计变量限制范围

针对该减阻优化问题，依据图5 所示的优化流程，首先采用拉丁超立方法进行DOE，构建样本空间，依次对设计点进行流体的数值仿真模拟，并依据约束条件筛选可行解集，依据该可行解集得到响应面代理模型，其中关于最大弯度及最大弯度位置的三维响应面如图6 所示。

图6 三维响应面示意

在得到响应面的基础上，依据二次拉格朗日非线性规划开展气动优化，该算例共计63 步迭代数，得到该代理模型下的最优解，并对该设计点进行单独的CFD 验证，最终输出优化后翼型。

3.2 优化结果对比

在该优化框架下翼型优化前后的参数见表5。从表5 可以看到，通过增加了向下的最大弯度，并将最大弯度位置后移，实现了在设计工况下翼型的阻力系数降低34%。

表5 优化前后参数对比

此处最大厚度优化前后无变化，是由于阻力基本与最大厚度正相关，而气动约束驱使迭代主要朝着减小阻力的方向开展优化，同时几何约束限制了厚度的减小，故在该优化算法下，最大厚度的初值即为优化后的可行解。

优化前后的外形及压力系数对比见图7 和图8。由外形图可以看到，因为NACA 翼型是根据中弧线叠加厚度的参数化表达方式，为此优化后的翼型依旧保有较高的光顺性，且优化后翼型前缘位置的下翼面厚度大于上翼面厚度[11]。此外，优化后的翼型具有较好的压力系数分布规律，并且最大相对弯度、最大弯度位置的微小改变便能引起较大的阻力系数变化，说明阻力系数对外形变化相当敏感。

图7 翼型优化外形对比

图8 翼型优化压力系数对比

对原始翼型和优化后翼型的流场计算结果进行基于Q 准则的涡识别，如图9 所示，其中Q＞0 表示存在涡结构。可以看到设计工况下的流场主要存在以下两方面区别：

图9 基于Q 准则的涡识别对比

1）翼型在13%弦长位置增加向下的弯度后，上曲面前部下移，更贴近设计工况的气流方向，使得优化后翼型上表面的涡区域更薄；

2）通过优化最大相对弯度及其位置，使翼型下表面后部涡尺度更小，且尾缘后方脱落涡位置延后，相较于原翼型，优化后翼型在尾缘较近区域内不存在明显的涡结构。

将优化结果与参考文献[12]中的结果进行对比见表6。从表6 结果可以看出，依据翼型参数化的表达方式能得到与CST 参数化方法效果相当（在文献中的对应设计变量下）的优化结果，而且仅采用3 个优化设计变量，说明该参数化表达方式的变量能够精准地体现出翼型设计的关键变量，这与飞机设计手册[13]中的表述规律一致，最大厚度在12%～17%时能得到较大的升力系数，而且最大弯度及最大弯度位置是影响阻力系数等气动特性的关键参数。

表6 不同参数化方法结果对比

4 结论

（1）NACA 翼型参数化方法的设计变量少且变量有明显的物理含义，参数的小幅变化能较大程度影响翼型的气动特性，此外基于该表示函数的翼型光顺性和连续性好，可以较为高效地构建用于优化设计的翼型，减少设计点数量，节省计算资源。

（2）结合基于拉丁超立方的DOE、响应面代理模型和NLPQL 优化算法的优化设计框架可以开展较为高效的优化设计。

（3）与较为通用的基于CST 参数化方法的翼型优化方法相比，翼型优化设计框架可以在部分工况下得到效果相当的优化结果，体现了该框架具有一定的实用价值。