数值模拟方法对NASA CRM模型阻力预测的影响

陈江涛，赵娇，章超，刘深深，张耀冰，吴晓军

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

气动力和力矩的精细预测是商用运输机气动性能评估和外形优化设计的关键因素之一。近些年来随着网格生成技术、流场求解和高性能计算机等因素的持续发展，CFD在预测复杂外形气动性能方面取得了长足的进步，评估三维真实运输机外形气动性能成为可能。为了评估现阶段CFD预测商用运输机气动性能的能力，AIAA应用空气动力学委员会在2001年发起了AIAA阻力预测会议，来自世界各地的研究机构和公司使用不同的程序、网格、湍流模型和数值格式对翼身组合体外形进行了模拟研究[1-6]。

在复杂工程外形的数值模拟中，网格类型、规模和分布、湍流模型、数值格式等，都会不同程度影响模拟结果。如何评估这些模拟方法对模拟结果的影响，并识别对模拟结果影响较显著的因素，对于后续CFD发展努力的方向有积极的借鉴指导意义。Rumsey和Allison[7]在抖振边界条件下对某民用运输机进行了数值模拟研究，考虑了不同软件、湍流模型和机翼后缘处的网格处理等影响因素，发现湍流模型对模拟结果的影响最大。国内也有很多学者从不同的角度，比如离散精度[8-9]、网格拓扑[10]、数值方法[11-12]等，研究了不同模拟方法对阻力会议外形模拟的影响。但是这些研究大多没有考虑多种因素同时作用下对预测的影响，比如Rumsey和Allison[7]在研究湍流模型对模拟影响的时候，保持使用的网格和程序不变，这样得到的结论有一定的局限性。因此如何在多种影响因素共同作用下评估模拟结果的变化规律，识别出对模拟影响较大的因素是本文试图解决的问题。

为了综合研究不同因素对商用运输机外形阻力预测的影响，本文以AIAA第六届阻力会议外形NASA CRM[13]为研究对象，同时考虑了网格、湍流模型、无黏通量格式和体心梯度求解方法等因素，使用枚举法和正交试验设计两种方法研究了不同因素对阻力预测的影响大小，识别出了关键因素，为下一步方法的发展指出了方向。

1 计算外形、网格和方法

本文考虑的外形是AIAA第六届阻力会议的NASA CRM模型，如图1所示。该模型代表了典型的现代跨声速运输机，设计巡航状态是马赫数Ma=0.85，升力系数CL=0.5。该外形在NASA NTF风洞、NASA Ames 11 ft风洞[14-15]和欧洲ETW风洞[16]都进行了测试，有大量数据可供CFD确认用。但是由于不同风洞运行环境的差异等因素，试验数据也存在着明显的不确定性。图2给出了NTF风洞和Ames风洞多个车次的试验结果对比，通过插值可以得到该模型在定CL=0.5条件下在NTF风洞测得的阻力系数CD为0.024 90～0.024 97，在Ames风洞测得的阻力系数CD为0.024 14～0.024 34。本文更多关注的是模拟方法对阻力预测的影响，没有通过试验结果评判模拟方法的优劣。

图1 NASA CRM外形

图2 CRM不同风洞试验数据对比

这里关注的是定升力系数CL=0.5条件下风洞试验状态的阻力预测，来流雷诺数为5×106。研究中综合考虑了网格规模、湍流模型、无黏通量格式和体心梯度求解方法对阻力预测的影响。下面分别介绍不同因素的情况。

首先使用了依次加密的3套非结构网格来进行模拟，网格严格按照组委会提供的规范来生成，加密过程中保持网格整体加密，最粗的一套网格表面如图3所示，网格信息列于表1。

计算使用课题组发展的MFlow程序[17], 这是基于二阶有限体积法的非结构网格求解器，能够处理多种网格单元类型。单元内使用线性重构使得空间精度达到二阶。研究中使用了两种梯度求解方法来计算体心的梯度，分别为基于节点的格林高斯方法[18](GN)和基于节点的最小二乘方法[19](LS)。考虑了3种不同的无黏通量格式，分别为Roe[20]、AUSMPW+[21](图片标识中简写为AUS)、HLLC+[22](图片标识中简写为HLL)。计算中假定流动为全湍流，使用SA一方程模型[23]和剪切应力输运(SST)k-ω两方程模型[24]。两种模型都是使用了QCR修正[25-26]之后的版本。研究发现[8,27-28]，对于该外形，QCR修正能够提高翼身结合处流动的模拟精度，计算结果与试验结果更加吻合。

图3 粗网格表面网格

表1 计算使用的网格信息

本文的计算从自由来流开始，所有的计算密度残差都下降了两个量级以上，气动力的振荡在0.2% 以内，如图4所示，因此计算收敛性引入的影响可以忽略不计。

图4 密度残差和气动力系数收敛曲线

2 阻力预测影响分析

2.1 枚举法结果

本文计算考虑了3种网格、3种无黏通量格式、2种湍流模型、2种体心梯度求解方法对阻力预测的影响。本节将所有因素可能的搭配都进行了模拟，计算安排如表2所示, 其中网格C代表粗网格，M代表中等网格、F代表密网格。

对于总阻力，使用粗网格、SA模型和GN方法的预测值高于上限值。对于压差阻力，使用粗网格、SSTk-ω模型、AUSMPW+和GN方法的预测值高于上限值，使用密网格、SA模型和LS方法的预测值低于下限值。对于摩擦阻力，均在上下限之间。需要指出的是，计算结果在统计的上下限之外不代表该结果是有问题的。

表2 计算安排设置

图5 阻力系数预测值的统计结果

图6给出了不同模拟方法搭配下的网格收敛性曲线。从图中可以看出随着网格加密，总阻力和压差阻力减小，摩擦阻力增大。摩擦阻力收敛图上SA和SSTk-ω的结果区分的很明显。压差阻力图上湍流模型和梯度求解方法影响都很大，特别是湍流模型。总阻力图上梯度求解方法影响很大，随着网格加密不同湍流模型的解之间的差别变小。

图6 阻力预测的网格收敛特性曲线

图7和图8分别给出了定升力状态下的迎角α和俯仰力矩系数Cm随网格加密的变化趋势。除了使用SA模型得到的迎角会有非单调变化,大部分结果都呈现单调的变化。总的来说，迎角和俯仰力矩随着网格加密变化都不大，SA模型计算得到的迎角因网格因素导致的变化不超过0.004°，SST模型为0.02°，SA模型得到的俯仰力矩系数变化不超过0.000 3，SST模型为0.002。网格收敛性研究基于Taylor展开，需要指定网格尺度h。在三维复杂外形计算中，经常使用N-1/3来近似h。假定粗、中、细3套网格计算得到的目标变量分别为f3、f2、f1, 3套网格的网格尺度分别为h3、h2、h1。通过简单的证明可以知道：只有当0<(f1-f2)/(f2-f3)ln(h2/h1)/ln(h3/h2), 这时计算结果随着网格加密呈现单调发散的趋势，基于Taylor展开的网格收敛性分析不再适用。这里需要指出的是，用N-1/3来近似h可能会给网格收敛性分析带来较大的误差，特别是在非结构混合网格的框架下，如何构造相容的网格序列开展网格收敛性研究也是验证与确认领域研究的难点。

图8 俯仰力矩系数预测的网格收敛特性曲线

对于不同模拟方法对阻力预测的影响，可以定性地做出判断。首先考虑湍流模型的影响，在固定网格、无黏通量格式和体心梯度求解方法的前提下，对于压差阻力，SSTk-ω模型比SA大，随着网格加密差量增大；对于摩擦阻力，SSTk-ω模型比SA小，但随着网格加密差距减小，总的量级大于压差阻力，因此总阻力SSTk-ω模型比SA小，这点在粗网格上更加明显，随着网格加密差距减小。

考虑梯度求解方法的影响：GN比LS压差阻力大，摩擦阻力稍大，总阻力大，这点在粗网格上和使用SSTk-ω模型时更加明显。

而无黏通量格式对阻力预测的影响几乎可以忽略不计。需要指出的是这一结论是基于使用的这3种格式而言的，如果使用其他与这几种格式性能有明显差别的格式，比如Steger-Warming、Van-Leer等，很有可能得到不一致的结论。

2.2 基于聚类算法的敏感性分析

4种模拟因素同时作用下，需要对每种因素的影响大小进行量化和判断，识别出对模拟结果影响较大的因素，为深入研究提供努力的方向。

本节使用了Kmeans聚类算法[29]来分析目标变量和输入因素之间的潜在联系。Kmeans聚类算法通过衡量不同点间的欧式距离来表征样本点之间的相似度，被聚为一类的样本点具有较高的相似度，通过分析同类样本中输入因素的特征，来辨别不同因素的影响大小。

分别针对总阻力、压差阻力和摩擦阻力，将所有样本点的目标变量分为两类，如图9～图11所示。黑色实线为大体的两类之间的分割线，黑色虚线分别为两类样本目标变量的中心值，不同颜色的标志代表各因素的不同方法。

以总阻力为例，图9给出了所有样本通过4个因素分类的情况，并用黑色实线将各个因素的结果分为两类。在图9(a)中第2类包括了全部细网格、大量中网格和少量粗网格样本，不能有效区分网格因素；图9(b)中SA样本几乎均等的分布在第1类和第2类，不能有效区分湍流模型的影响；图9(c)中3种无黏通量格式样本同样几乎均等分布在第1类和第2类，不能有效区分无黏通量格式的影响，在图9(d)中，第1类大部分为GN结果，第2类大部分为LS结果，起到了有效区分，所以梯度求解方法对于总阻力的影响最大。类似地，对于压差阻力，湍流模型发挥了最大的影响，对于摩擦阻力，湍流模型发挥了最大的影响。

图9 总阻力系数的聚类分析结果(分2类)

本文计算中，有的模拟因素比如网格和无黏通量格式，都存在3种选项，是否应该将所有样本归为3类，结果如图12～图14所示。对于总阻力，梯度求解方法发挥了最大的影响，第1类全部为GN的结果，第3类全部为LS的结果，而其他因素在3个 类里都存在不可忽视的数量标识。对于压差阻力，湍流模型发挥了最大的影响，第1类全部为SSTk-ω结果，第3类全部为SA结果，而其他因素在3个类里都存在不可忽视的数量标识。对于摩擦阻力，同样是湍流模型发挥了最大的影响，第1类 全部为SA结果，第2类和第3类全部为SSTk-ω结果，而其他因素在3个类里都存在不可忽视的数量标识。

在本算例中，将所有样本点不管是分为2类还是3类，都可以得到一致的敏感性结论。即对于总阻力，梯度求解方法发挥了最大的影响。对于压差阻力和摩擦阻力，湍流模型发挥了最大的影响。

2.3 基于McKay的敏感性分析

聚类算法只是定性地给出了因素的敏感性分析，下面使用McKay主影响分析来量化4种因素各自的敏感性，定义为

图10 压差阻力系数的聚类分析结果(分2类)

(1)

式中：V(Y)表示目标变量的方差估计；V[E(Y|Xi)]表示目标变量在因素Xi作用下条件期望的方差估计。该指标量化了因素Xi导致的目标方差的变化。表3分别给出了压差阻力、摩擦阻力和总阻力预测中各因素的敏感性指标。

从表3可以分析出，对于压差阻力，湍流模型对预测的影响最大，其次是网格和梯度求解方法，无黏通量格式的影响很小；对于摩擦阻力，湍流模型对预测的影响最大，网格、无黏通量格式和梯度求解方法对结果的影响很小，均可以忽略不计；对于总阻力，梯度求解方法的影响最大，其次是网格和湍流模型，无黏通量格式的影响可以忽略不计。

图12 总阻力系数的聚类分析结果(分3类)

图13 压差阻力系数的聚类分析结果(分3类)

虽然在压差阻力和摩擦阻力预测中，湍流模型都是影响最大的因素，但是在总阻力预测中，湍流模型并不是影响最大的因素。这里考虑是因为湍流模型对压差阻力和摩擦阻力的影响趋势相反。具体来说，SA模型预测的压差阻力较SSTk-ω模型预测偏小，但摩擦阻力预测的偏大，这样导致湍流模型导致的总阻力变化较小。而梯度求解方法对于压差阻力和摩擦阻力的影响是一致的，GN预测的压阻和摩阻都比LS偏大，所以总阻力预测中，这一趋势被放大。

图14 摩擦阻力系数的聚类分析结果(分3类)

表3 McKay主影响指标的参数敏感性

需要指出的是对模拟因素敏感性的分析都是基于使用的这几种网格、模型和格式，对这一特定算例分析得到的。如果使用不同的模拟方法，比如其他的网格、模型和格式，或者针对新的算例，可能会得到不一致的结论。但本文给出了一般的因素敏感性分析方法，针对具体问题，可以仿照这一流程进行分析。

2.4 正交试验设计分析

在上文的分析中，考虑了4种影响因素对阻力预测的影响，每种因素考虑了2～3种不同的选项，所有可能选项遍历下来计算量不小，因此在本节中尝试使用更有效率的正交试验设计[30]来完成上述分析。在正交试验设计中，试验因素为网格、湍流模型、无黏通量格式和梯度求解方法，各因素划分为表4所示的不同水平因子。

首先不考虑各个因素之间的交互作用，使用极差来辨别因素的主次，选择正交表L9(34)，采用拟水平法进行正交试验设计，试验设计方案如表5所示。

(2)

式中：pj为第j个因素的水平数目。极差越大，则说明该因素对结果的影响越大，即该因素越重要。

表4 正交试验设计的因素与水平因子

表5 拟水平法下的L9(34)正交表

通过极差分析可得不同因素对结果影响的主次顺序：① 对于压差阻力，湍流模型的影响>网格的影响>梯度求解方法的影响>无黏通量格式的影响；② 对摩擦阻力，湍流模型的影响>网格的影响>梯度求解方法的影响>无黏通量格式的影响；③ 对总阻力，梯度求解方法的影响>网格的影响>湍流模型的影响>无黏通量格式的影响。这里得到的结论与枚举法的敏感性分析结果一致，但是计算量只是枚举法的1/4。

除此之外，还可以看出在不同水平下网格、湍流模型对压差阻力和摩擦阻力的影响趋势是相反的，而梯度求解方法对压差阻力和摩擦阻力的影响趋势是相同的，这也解释了在总阻力预测中，梯度求解方法影响最大的原因。

图15 不同因素不同水平下阻力均值的变化曲线

上面的分析中没有考虑各个因素之间的交互作用。从枚举法的结果中可以看到，两个因素之间存在一定的交互作用，如图16所示，可以明显看出网格和湍流模型之间存在一定的交互作用，这在统计学上称为正交互。

为了进一步考查各因素之间的交互作用贡献，根据混合水平正交表L18(2×37)安排试验，如表6所示，第3、6、7、8列作为误差列。需要说明的是，混合水平正交表L18(2×37)可以考察置于第1列的二水平因子与置于第2列的三水平因子的交互作用，但不能考察其他列间的交互作用，故

图16 湍流模型和网格因素的交互作用

表6 考虑两因素交互作用的拟水平法正交试验设计

此处只求解湍流模型和网格的交互作用，梯度求解方法与网格、湍流模型间的交互作用可用相同方法进行分析。

对此可以通过衡量各因素的贡献率来评估因素的重要性，贡献率的详细推导过程可以查阅文献[31]，各因素对阻力预测的贡献率分析结果如表7所示。

通过正交表结果分析，可以得到以下结论：

1) 贡献率分析的结果与极差分析保持一致：对压差阻力，湍流模型对指标影响最大，其次是网格；对摩擦阻力，湍流模型对指标影响最大；对总阻力，梯度求解法对指标影响最大，这也与枚举法的敏感性分析结果一致。

2) 梯度求解方法与网格、梯度求解方法与湍流模型的交互作用对压差阻力、摩擦阻力以及总阻力的贡献率均可忽略；湍流模型与网格的交互作用引起的数据波动只对总阻力有较显著的贡献，但也与误差引起的数据波动的贡献率差不多，故也可忽略。

表7 阻力预测的贡献率分析

3 结 论

使用枚举法来考察所有因素的综合作用，通过基于聚类分析的定性敏感性分析方法和基于Mckay主影响分析的定量方法，识别出了对总阻力预测影响较大的因素为梯度求解方法。而湍流模型是影响压差阻力和摩擦阻力预测的最关键因素。

为了减轻计算负担，使用正交试验设计方法完成上述研究，敏感性分析结果与枚举法基本一致，但计算量大大减少，这位解决类似的因素影响分析提供了切实可行的思路。

需要指出的是，本文得到的模拟因素敏感性结论都是基于使用的这几种网格、模型和格式，对这一特定算例分析得到的。如果使用不同的模拟方法，比如其他的网格、模型和格式，或者针对新的算例，可能会得到不一致的结论。但本文给出了一般的因素敏感性分析方法。针对具体问题，可以仿照这一流程进行分析。