湍流模型系数不确定度对翼型绕流模拟的影响

赵辉，胡星志，张健，陈江涛，马明生

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

近些年来，随着网格生成、流场求解和高性能计算机等学科的快速发展，计算流体力学（Computational Fluid Dynamics，CFD）在工业领域发挥了越来越重要的作用。但是在复杂工程外形的数值模拟过程中，存在着多种来源的不确定性因素，包括几何外形、网格分布、物理模型及数值计算方法等［1－2］，这使得CFD给出的最终结果也存在着复杂的不确定性。在工程外形的评估或者优化问题中，如果不考虑CFD结果的不确定性，往往不能准确反映复杂物理条件下真实外形的气动性能。因此CFD的不确定度量化成为近期研究的热点。

对于工程中的湍流问题，通常都是通过湍流模型来模化雷诺应力，封闭雷诺平均方程。湍流模型是典型的认知不确定性问题。由于对湍流物理本质的理解不完备，导致了湍流模型存在明显的不确定性。比如模型中的系数很多都是通过试探方法和经验给出，或者应用渐进性原则通过简单流动的直接数值模拟或试验结果来确定的［3］。但是这样得到的湍流模型并不能保证对所有流动情况都适用。因此评估湍流模型系数的不确定性对数值模拟结果的影响对于CFD的确认问题或者工程外形的设计优化都十分重要。

在常用的不确定度量化方法中，蒙特卡罗（Monte Carlo，MC）方法［4］由于构造简单、易于实现被广泛使用。但是MC方法收敛速度慢，需要大量的采样点才能得到满足精度要求的结果，因此限制了其在大维度多变量不确定性问题中的应用。随机谱方法，包括本文要讨论的多项式混沌（Polynomial Chaos，PC）方法［5－6］，是一种更有效率的不确定度量化方法。PC方法将整个数值模拟过程认定为一个随机过程，根据输入参数的概率密度分布函数，选择合适的正交基函数，然后将该随机过程的输出在基函数构成的谱空间中展开。通过嵌入式或者非嵌入式方法，确定展开式中的自由度，从而得到输出的完整的谱空间展开和各种统计信息，比如平均值、均方根以及每个不确定输入变量的全局敏感度指标。使用PC方法的不确定度量化已经有了很多成功的应用［7－14］。

本文使用非嵌入式PC方法研究了Spalart-Allmaras（SA）模型系数的不确定度对RAE2822翼型数值模拟的影响。首先，介绍了非嵌入式PC方法和使用的流场解算工具，研究从单变量入手，考察了卡门常数的不确定度对升阻力系数、压力分布、摩擦系数分布和马赫数的影响。然后，同时考虑了SA模型中9个参数的不确定度对数值模拟的影响。最后，给出了研究结论。

1 PC方法展开

PC方法源于Wiener提出的各向同性混沌理论［15］，用于湍流问题中随机变量的谱空间展开。Ghanem和Spanos将该方法应用到有限元计算中，用来分析固体力学中的不确定性问题［16］。Xiu和Karniadakis将适合高斯随机过程的 Hermite PC推广到 Wiener-Askey PC，适用于满足更一般的概率密度分布的随机过程［6］。PC方法基于不确定量的谱空间展开，将随机量分解为确定性的和随机的两部分。下面首先从单输入变量入手，介绍PC方法。

假定随机的输入变量ξ，满足某种概率密度分布f（ξ）。对于任意的随机输出变量y（可以是升力、阻力系数等积分量，也可以是流场中某处的压力系数、摩擦力系数或者马赫数，或者分离区长度、再附点位置等流场特征量），可以在输入变量ξ的正交基函数序列构成的谱空间中展开，即

式中：αj为确定性分量，也可称为自由度；ψj为随机变量的基函数。理论上式（1）的展开有无穷多项，但是在实际应用中通常将展开在某阶次处截断，即

正交基函数的选择由随机输入变量的概率密度函数决定，满足：

即 ｛ψj，j≥0｝是带权函数为f（ξ）的正交多项式序列。对于均匀分布的随机变量，最优的基函数序列是Legendre多项式；而对于正态分布的随机变量，最优的基函数序列是Hermite多项式。这里谈到的最优，是指系统输出的统计信息随着多项式阶数增加而收敛的速度与理论一致。

根据与求解器耦合方式的不同，可以将PC分为嵌入式多项式混沌（Intrusive Polynomial Chaos，IPC）法和非嵌入多项式混沌（Non-Intrusive Polynomial Chaos，NIPC）法。在IPC中，控制方程中的不确定变量和参数都用其PC展开替代，经过投影之后，得到与原控制方程相似的随机控制方程，自变量为PC展开中的自由度。IPC需要对现有代码做较大的改动，不仅工作量大，而且很多时候现有代码都是封闭的，这限制了IPC在工程问题中的应用。而NIPC无需修改控制方程，将现有的求解器视为一个黑匣子。首先在输入变量的随机空间里通过恰当的抽样方法获得一系列的采样点，传递给求解器，进行确定性的计算，得到各采样点对应的输出响应值，然后通过下文介绍的投影法或回归法，得到完整的PC展开。本文主要讨论NIPC方法。

在NIPC中，主要有两种方法求解PC展开式中的自由度，投影法和回归法。投影法利用基函数的正交特性，在式（2）的左右两侧同时乘以ψj，并在ξ的支撑集内积分，利用基函数的正交性质，可以得到

式中：右侧积分〈yψj〉可以通过蒙特卡罗积分或者高斯积分得到。当输入变量的维数较大时，投影法需要大量的采样点才能满足精度要求，即使使用稀疏网格的高斯积分，计算量也是巨大的。

在回归法中，通过确定性的CFD计算得到采样点 ｛ξ1，ξ2，…，ξN｝对 应 的 响 应 值 ｛y1，y2，…，yN｝。通常情况下，采样点的数目大于PC展开中自由度的个数。因此形成一个超定方程组Ψα＝Y，其中Ψ是测量矩阵，Ψij＝ψj（ξi）。通过最小二乘法可以得到自由度组成的向量α。本文中自由度计算都是采用回归法得到。在单变量的不确定度量化中，使用高斯积分法进行结果的对比。

采样点的数目N跟PC展开中的自由度个数有关，即

式中：n为输入随机变量的维数；p为混沌多项式的阶数；np定义为过采样率。Hosder等［17］在计算中发现，np＝2时，PC展开对输出变量统计信息的近似较好。因此在本文的计算中np均取为2。

得到输出的PC展开之后，可以根据基函数的正交特性直接得到输出的均值μ、均方根σ等统计信息，即

对于多输入变量的情况，如果各个变量之间是相互独立的，多维基函数可以直接取为单维基函数的乘积形式，联合概率密度函数为各变量概率密度函数的乘积。

在多输入变量的情况下，每个变量的不确定度对系统输出的不确定度的贡献大小不一，很多时候需要量化每个变量的相对贡献大小，分析输入变量的全局敏感性。Sobol指标［18］提供了工具来量化每个输入变量对系统输出方差的贡献程度。Sobol指标通过Sobol分解得到，是一种基于方差分解的全局敏感性分析方法，其优点是考虑了随机输入参数在整个取值范围内对输出响应的贡献以及随机输入参数的交互作用［17］。Sobol指标的详细推导过程和具体形式可以参考文献［18-19］。

2 计算外形和方法

本文考虑的算例是RAE2822翼型［20］，这是一个典型的超临界翼型，也是测试CFD程序跨声速流动模拟能力的经典算例，国内外研究机构针对该翼型开展了大量的风洞试验和数值模拟研究。本文的计算状态是：马赫数Ma＝0.729，雷诺数Rec＝6.5×106，迎角α＝2.31°。使用的网格如图1所示。

计算使用的程序是课题组自行研发的MFlow［21］。该程序采用基于格心的有限体积法，能够处理多种网格单元类型（六面体、四面体、三棱柱、金字塔以及使用几何多重过程中融合产生的多面体）。单元内使用线性重构达到二阶精度。梯度计算使用基于格点的格林高斯方法［22］。使用Venkatakrishnan限制器［23］来抑制间断附近的振荡。使用Roe格式来计算无黏通量。计算使用一阶欧拉后差来推进到收敛状态，使用局部时间步长加速收敛过程。Jacobian通量使用一阶迎风格式得到。

图1 RAE2822翼型网格分布Fig.1 Grid distribution for RAE2822airfoil

计算使用SA一方程模型［24］，这是在工程界应用非常广泛的湍流模型。在使用中，假定流动为全湍流，忽略了原始模型中的转捩项（Trip Term）。因此模型中有9个常数，分别为：cb1、σ、cb2、κ、cw2、cw3、cv1、ct3、ct4。SA 模型的输运方程形式和模型系数是通过量纲分析、伽利略不变性和部分试验（包括二维混合层、尾流和平板边界层等）结果给出的［25］。模型参数的取值存在着一定的不确定度。比较典型的是卡门常数κ，κ描述了边界层内平均流向速度U＋随离开壁面距离y＋变化的对数关系［26］。在大气表层试验中，κ的值最大为Sheppard测得的0.46，最小为Kansas的结果0.35［27－28］。也有学者指出κ跟雷诺数和压力梯度有关［29－30］。佘振苏给出了规范壁湍流中普适卡门常数（κ＝0.45）的定义［31－32］。本文参考了文献［33］中8个变量的取值范围，如表1所示，假定每个参数在各自的支撑集内为均匀分布。

表1 SA模型的标准参数和支撑Table 1 SA model closure coefficients and associated support

在进行不确定度量化之前需要指出的是本文侧重于研究SA模型系数的不确定度对该算例模拟结果的不确定度影响，没有考虑网格分布、数值格式等方面的不确定度带来的影响。

3 不确定度分析

3.1 单变量

从单变量入手，首先研究了卡门常数的不确定度对RAE2822翼型模拟的影响。由于事先无法知晓系统输出对卡门常数的响应关系，所以假定PC展开的截断阶次分别为p＝1，2，3，4。图2和图3给出了翼型升力、阻力系数CL和CD的平均值和均方差值结果。可以清楚地看到，随着PC展开的阶次增加，升力、阻力系数的统计信息变化很小，特别是平均值只有0.01%以内的差别。图4给出了所有采样点的升力、阻力系数随着卡门常数变化的曲线，据此可以推断升力、阻力系数对卡门常数是近似线性的响应关系。图5给出了升力、阻力系数的PC展开中各项自由度随多项式阶次的变化，二阶和二阶以上基函数对应的自由度比一阶基函数对应的自由度量值小两个量级以上，进一步证实了上述结论。

图2 升力系数的平均值和均方根值Fig.2 Mean values and standard deviations of lift coefficients

图3 阻力系数的平均值和均方根值Fig.3 Mean values and standard deviations of drag coefficients

图4 所有采样点的升力、阻力系数Fig.4 Lift and drag coefficients for all samples

图5 升力、阻力系数的自由度Fig.5 Freedoms of lift and drag coefficients

得到输出的PC展开后，可以预测任意一点的输入对应的输出响应。表2给出了κ为模型标准取值0.41时，PC预测值和CFD计算的结果，两者的误差几乎可以忽略不计。说明建立的PC展开能够准确地描述输出对输入的响应。

在单变量情况下，利用高斯积分法和回归法计算PC展开的自由度，计算量相当。图6给出了两种方法得到的升力、阻力系数PC展开中各项自由度，两者结果也基本一致。

在上面的分析中，系统输出关注的是升力、阻力系数等积分量，下面进一步讨论流场中某处的压力系数、摩擦力系数或者马赫数。图7和图8给出了壁面压力分布的不确定度分析结果。其中包括压力分布的试验值、平均值、最大值、最小值和方差值。κ的不确定度主要影响翼型上表面，特别是激波附近的压力分布。减小κ会使得激波位置后移。

图9给出了壁面摩擦系数Cf预测的均值和极值。κ的不确定度影响着翼型上、下表面的摩擦系数预测。当κ在其支撑集内取最小值时，预测的壁面摩擦系数最小；当κ取最大值时，预测的壁面摩擦系数最大。对于激波位置和壁面摩擦系数随着κ变化的原因，可以从湍流涡黏系数υt的变化来解释。在简单平面剪切湍流的等雷诺应力层，涡黏系数υt与κ成正比［3］，即υt＝κuτd，其中uτ为壁面摩擦速度，d为到壁面的距离，因此减小κ会得到较小的涡黏系数。图10给出了κ取最大值时计算的涡黏系数与κ取最小值时的涡黏系数的差值云图。从图可以清楚地看到，在流场的大部分区域，减小κ会减小涡黏系数，从而导致预测的壁面摩擦系数减小，同时边界层内流体微团动能损失较少，抵御逆压梯度的能力更强，预测的激波位置因此后移。

表2 PC展开与CFD得到的升力、阻力系数比较Table 2 Comparison of lift and drag coefficients obtained with PC and CFD

图6 升力、阻力系数的自由度比较（回归法和高斯积分法）Fig.6 Comparison of freedoms of lift and drag coefficients（regression method and Gauss integration method）

图7 单变量下的压力分布Fig.7 Pressure distribution under single variable

图8 压力分布的方差Fig.8 Variance of pressure distribution

图11是马赫数方差的云图，κ的不确定度主要影响上表面激波位置的预测，从而也会影响激波附近马赫数的预测。

图9 单变量下壁面摩擦系数的平均值和极值Fig.9 Mean and extreme values of skin friction coefficient under single variable

图10 涡黏系数差值云图（υt｜κ＿max－υt｜κ＿min）

图11 马赫数方差云图Fig.11 Contours of variance of Mach number

3.2 多变量

通过上面的分析可以知道κ的变化对激波位置和壁面摩擦预测的影响，也得到了系统输出对κ的响应。下面同时考虑SA模型中9个参数的不确定度对数值模拟的影响。计算中PC展开的阶次为p＝2，考虑了任意两个参数的联合影响。在9维空间中，采样点通过拉丁超立方抽样得到。

表3给出了升力、阻力系数预测中9个参数各自的Sobol指标，按照贡献大小排序。与文献［33］的结果存在明显的差异，主要表现为在本文的计算中ct3、ct4这两个参数对升力、阻力系数方差的贡献很小。这是因为本文使用的SA模型忽略了原始模型中的转捩项，由此很多学者指出包含ct3、ct4的ft2也是可以忽略的。同时cb2、cw3这两个参数的贡献也是可以忽略不计的。cb1、σ、κ、cw2、cv1这5个参数对于升力、阻力系数的不确定度影响都相对比较大，κ是其中贡献最大的参数，其他4个参数的贡献大小排序并不完全一致。

表3 升力、阻力系数预测中9个参数的Sobol指标Table 3 Sobol indices of 9closure coefficients for lift and drag coefficients

通过升力、阻力系数的响应，可以得到其平均值、均方差等统计信息，如表4所示。其中升力系数 的 最 大 值 取 在 ｛cb1＿max、σ＿max，κ＿min，cw2＿max，cv1＿max｝，最 小 值 取 在 ｛cb1＿min，σ＿min，κ＿max，cw2＿min，cv1＿min｝。阻力 系 数 的 最 大 值 取 在 ｛cb1＿max，σ＿max，κ＿max，cw2＿max，cv1＿min｝，最 小 值 取 在 ｛cb1＿min，σ＿min，κ＿min，cw2＿min，cv1＿max｝。ct3、ct4、cb2、cw3这 4 个 参 数 在支撑集内变化对升力、阻力系数影响很小。

表4 升力、阻力系数的均值和均方差值Table 4 Mean value and standard deviation of lift and drag coefficients

得到输出的PC展开后，可以预测任意一点的输入对应的输出。表5给出了模型参数为标准取值时，PC预测值和CFD计算的对比，其中，CDp为压差阻力，CDf为摩擦阻力，两者误差很小。说明建立的PC展开能够准确地描述输出对多变量输入的响应。

表5 PC展开与CFD得到的升力、阻力系数比较Table 5 Comparison of lift and drag coefficients obtained through PC and CFD

通过确定的CFD计算可以得到壁面压力分布响应。图12给出了压力的试验值、平均值和极值。在机翼的上表面，特别是激波附近，压力的不确定度表现显著。上表面压力极小值取值在｛cb1＿max，σ＿max，κ＿min，cw2＿max，cv1＿max｝，极 大 值 取 值 在｛cb1＿min，σ＿min，κ＿max，cw2＿min，cv1＿min｝，与 升 力 系 数 的情况类似。图13给出了9个输入参数对压力方差的贡献度。与升力、阻力系数贡献结论一致，ct3、ct4、cb2、cw3这4个参数对压力的不确定度贡献相对较小，其他5个参数，特别是σ、κ贡献较大。

图14给出了壁面摩擦系数的平均值和极值。与压力分布不同的是，ct3、ct4、cb2、cw3这4个参数在机翼的前缘对摩擦系数的不确定度也有着不可忽视的贡献，如图15所示。

Fig.12 Pressure distribution under multi-variables

图13 压力分布的Sobol指标分布Fig.13 Sobol index distributions for pressure distribution

图14 多变量下壁面摩擦系数的平均值和极值Fig.14 Mean and extreme values of skin friction coefficient under multi-variables

图15 壁面摩擦系数分布的Sobol指标分布Fig.15 Sobol index distributions for skin friction coefficient distribution

SA模型中的9个参数综合影响了湍流涡黏系数的生成、扩散和耗散。比如参数cb1与涡黏系数的生成项和耗散项（部分）成正比，该参数的改变会同时增大或者减小生成和耗散，对涡黏系数的综合作用尚且未知。从理论上很难得到参数变化对涡黏系数影响的规律，而且很可能在流场不同区域产生相反的结论。因此本文只对κ、ct3、ct4这3个参数引起的不确定度进行简单分析，其他参数的影响需要继续深入研究。

参数ct3、ct4对模型的影响体现在生成衰减项ft2里。本文使用的SA模型忽略了原始模型中的转捩项，由此ft2只在翼型前缘附近涡黏系数很小的区域发挥作用，在其他区域的边界层内基本没有作用。因此ct3、ct4的不确定度对模拟的影响主要体现在前缘附近，特别是对摩擦系数的预测。

图16给出了马赫数的平均值和方差的云图，在激波附近马赫数表现出较大的不确定度，与设想的一致。

4 结 论

本文使用非嵌入式多项式混沌方法研究了湍流模型系数的不确定度对RAE2822跨声速翼型绕流模拟的影响。

Fig.16 Contours of mean values and variances of Mach number

1）从单输入变量入手，研究了卡门常数κ的不确定度对数值模拟的影响。研究发现升力、阻力系数对κ是近似线性的响应关系。在单变量情况下，利用高斯积分法和回归法计算PC展开的自由度，两者结果也基本一致。κ的不确定度主要影响翼型上表面，特别是激波附近的压力分布，同时还影响翼型上、下表面的摩擦预测结果和激波附近马赫数的预测结果。κ减小时，预测的激波位置后移，摩擦系数预测值变小。

2）同时考虑模型中9个参数的不确定性带来的影响。通过升力、阻力系数预测中9个参数各自的Sobol指标可以发现，ct3、ct4、cb2、cw3这4个参数对升力、阻力系数方差的贡献很小，可以忽略不计。其他5个参数对于结果的不确定度影响都相对比较大，κ是其中贡献最大的参数。建立的PC能够给出壁面压力分布、摩擦系数分布和马赫数分布的不确定度。

3）通过PC展开和CFD计算的对比证实了建立的PC展开能够准确地描述输出对输入的响应。这为量化计算结果的不确定度提供了切实可行、有较高效率的解决方案。需要指出的是本文的计算只考虑了RAE2822跨声速翼型模拟的单一计算状态，影响规律是否可以推及其他工况和算例需要进一步检验。