考虑数值离散误差的湍流模型选择引入的不确定度量化

陈江涛，章超，吴晓军，赵炜

1. 中国空气动力研究与发展中心 空气动力学国家重点实验室，绵阳 621000

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

近些年来，随着流场解算器、网格生成、前后置处理和高性能计算机等学科的快速发展，计算流体力学(Computational Fluid Dynamics，CFD)在工业领域得到了广泛应用，也发挥了日益重要的作用。但目前中国在CFD软件体系化、工程化及大规模推广应用等方面还与国外先进水平存在较大差距。因此中国2018年底启动了国家数值风洞(National Numerical Windtunnel，NNW)工程。该工程由中国空气动力研究与发展中心牵头建设，旨在建成拥有自主知识产权、中国开放共享的空气动力数值模拟平台，满足航空航天等领域对自主CFD软件产品的迫切需求。经过两年多的努力，NNW工程在模型算法、软件研发等一系列关键技术上取得了长足的进展[1]。

在复杂工程问题的数值模拟过程中，存在着多种来源的不确定性因素，包括网格、物理模型、数值离散方法等。在CFD中，数值、模型选择和模型预测偏差是3类重要的不确定因素。数值离散引入数值误差，但对于工程问题数学模型的真实解通常未知，因此不能准确估计数值离散误差，只能采用数值不确定度进行表征。CFD中模型选择引入的不确定性非常普遍。对同一种自然现象，通常存在多种可能的物理模型模化，如各种湍流模型、亚格子模型和化学反应模型等，不存在唯一的模型。人们试图在可能的备选模型中选择最好的模型，但由于工程问题的复杂性以及使用者认知的局限性，这种选择本身也有一定的不确定性。对于某一个特定的模型，即使能够准确求解，其预测结果也可能与真实现象存在偏差。很多情况下这种偏差不能准确地量化(因缺少试验数据或者试验存在显著的不确定度)，转化为模型预测偏差的不确定度。

不确定因素使CFD模拟结果也存在着不可忽视的不确定度。传统CFD应用过程中，从业者基于工业部门提供的数模，从经验或理论分析出发生成分布直觉上合理的网格，使用常用的数学模型和数值格式运行程序得到模拟结果。在整个过程中，没有考虑到不确定性因素导致的结果的不确定度。这也给工业产品设计优化、性能评估和模型确认等活动带来潜在的风险。科学、定量地描述这些因素对模拟结果的影响已成为科学计算领域的重点方向。NASA[2]认为，误差和不确定度管理是未来CFD需要具备的六大基本能力之一。由于CFD中存在显著的不确定因素，无法准确地给出关注的气动力结果(给出的是某种模型、参数和离散方法下的数值结果)，但通过系统的不确定度量化工作可以给出气动力的置信区间或概率分布。这也可能改变CFD给工业部门提供数据的传统方式。

针对模型选择引入的不确定度，通常可用两种方法进行量化：调整参数法[3-5]和贝叶斯模型平均方法[6-8]。调整参数法在缺少试验数据的情况下，通过专家经验给备选模型中的每一个模型分配相应的概率，通过平均得到目标量的均值和标准差，以此量化模型选择的不确定度；该方法虽然取得了部分成功应用，但是过分依赖于专家的经验，不同专家可能给出不同的概率分配，这进一步引入了不确定性。贝叶斯模型平均方法在获得试验结果的情况下，通过构建模型似然函数评价模型预测和试验结果的吻合程度，使用贝叶斯定理更新模型的概率；该方法能够比较客观地分配每个模型的概率，得到了广泛的应用。

贝叶斯模型平均方法依赖于模型预测结果分配模型后验概率，但在CFD领域通常无法获得模型的真实解，CFD得到的是该模型在某种特定数值离散(格式、网格)下的数值解。使用数值解代替模型真实解会给模型后验概率计算带来未知的不确定性，但现有应用算例均未考虑数值离散对模型后验概率的影响。为考查这一影响，本文发展考虑数值离散误差的贝叶斯模型平均方法，使用概率的上、下限表征数值离散和模型选择两者的总不确定度。首先对贝叶斯模型平均方法和数值离散误差估计方法进行介绍，然后提出考虑数值离散误差的模型选择不确定度量化方法，最后通过NACA0012翼型和CHN-T1模拟两个算例分析湍流模型选择中的不确定度。

1 贝叶斯模型平均方法

在存在多种可能模型的情况下引入模型概率的概念，定义为在备选模型集合中某一模型是最好模型的可信度。给定试验数据集D和包含n个备选模型的集合，使用贝叶斯定理将模型的先验概率更新为后验概率：

(1)

式中：P(Mk|D)为给定试验数据集D的情况下模型Mk的后验概率；P(D|Mk)为在给定试验数据D的情况下模型Mk的似然函数，即该模型预测结果为试验结果的概率，这项是贝叶斯模型平均方法的关键，需要建立计算数据和试验数据的概率关系，计为L(Mk|D)；P(Mk)为模型Mk的先验概率，即在获得试验数据之前根据专家意见或问题分析分配的概率，当对模型没有明显偏好时，可认为每个模型的先验概率相等，即P(Mk)=1/n。通常做法是假定模型预测值和试验值之间的差异为均值为0的高斯随机变量，即

(2)

(3)

(4)

此时

(5)

可以看到似然函数与模型预测值和试验的偏差成反比，与直觉想象一致。

考虑到试验数据有m个状态点的情况，假定每个状态都是独立无关的，有

(6)

现可以通过式(1)给出模型的后验概率。综合各个模型的后验概率得到一个新的融合模型，其预测值的概率描述为

(7)

由于每个模型输出的概率表示P(y|Mk,D)均为独立的高斯随机变量，因此总的输出也是高斯分布，其均值和方差为

(8)

在实际工程中，比较关注的是预测值在一定置信水平下的置信区间。由于输出变量满足高斯分布，根据其均值和方差得到一定置信水平下的置信区间非常方便，比如95%的置信区间为

(9)

值得注意的是，通过简单运算可知，输出变量的统计特性只与试验值和备选模型集合中每个模型的预测值有关。这也为第3节构建目标量的概率边界提供了方便。

2 数值离散误差估计方法

数值离散包括离散格式和计算域离散，必然引入数值误差。所有针对模型的验证与确认活动(如模型确认、模型平均等)都必须考虑数值误差的影响，将模型结果和模拟结果区分开。数值误差估计就是通过一系列的数值计算估计模型方程真实解的情况。但对于工程问题，准确的数值误差估计非常困难，此时数值误差转化成数值不确定度问题。数值不确定度属于认知不确定度，因为数值解和模型真实解的误差是定值，但是由于认知水平局限不知道模型真实解，所以也无法准确估计离散误差。通过一系列逐渐加密的自相似网格进行Richardson插值分析是数值离散误差估计的主要手段[9-11]。

Richardson插值是基于Taylor展开的，假定离散解满足

(10)

式中：fi为在该套网格上关注输出量的离散解；f0为关注输出量的真实解；α为未知常数；hi为网格的特征尺度，在二维计算中，通常使用N-1/2近似(N为单元总量)，在三维计算中，使用N-1/3近似；p为网格加密时数值误差的收敛速度。该展开有3个未知数(f0、α和p)，理论上只需要在3套满足自相似和全计算域同一比例加密的网格上进行计算即可。但是要求3套网格的离散解都位于真实解的渐进收敛区域内。假定粗、中、细3套网格计算得到的目标变量分别为f3、f2、f1，3套网格的网格尺度分别为h3、h2、h1。通过简单的证明可知，只有当0<(f1-f2)/(f2-f3)0，该展开才是合理的。

但即使实际收敛速度p>0，也要考虑与理论收敛速度的差异。当实际收敛速度与理论收敛速度差异较大时，会对误差估计产生显著的影响。CFD软件使用二阶的空间离散格式。密网格离散解的相对误差定义为

(11)

通过简单的运算可以得到相对误差，也可以写为

(12)

当实际收敛速度0.5≤p≤2.0时，通过式(10) 进行离散误差估计是合理的；当p>2.0时，数值误差估计过小；当0

因使用软件的理论精度为二阶，考虑使用基于一阶、二阶和一阶/二阶混合展开进行误差估计，分别如式(13)～式(15)所示。

一阶展开为

fi=f0+αhi

(13)

二阶展开为

(14)

一阶/二阶混合展开为

(15)

式中：α′、β和β′均为展开常数。

式(10)、式(13)～式(15)给出了4种离散误差展开形式，具体使用哪种展开需要根据离散解收敛性质和展开拟合误差决定。

首先定义4种展开的拟合误差。当使用ng(>3)套网格时，拟合误差标准差的无偏估计分别为

(16)

式中:下标1、2、12和p分别表示一阶、二阶、一阶/二阶混合展开及使用式(10)的p阶展开。

当只有3套网格时，σp和σ12均为0。

已知至少3套自相似和全计算域同一比例加密的网格上的离散解fi和相应网格的全局网格尺度hi，对空间离散理论精度为二阶的CFD软件进行数值离散误差估计的流程如下：

1) 根据式(10)拟合得到实际收敛阶p。

2) 确定拟合关系式：

① 当0.5≤p≤2.0时，确定式为最终的展开拟合关系式。

② 当p>2.0时,在一阶展开和二阶展开中，选择拟合误差最小的展开为最终的拟合关系式。

③ 当0

④ 当p≤0时,在一阶展开、二阶展开和一阶/二阶混合展开中，选择拟合误差最小的展开为最终的拟合关系式。

3) 根据选定的展开关系式进行最小二乘拟合，得到f0、拟合值fi,fit和拟合误差标准差

求解的最小二乘问题分别为

(17)

4) 根据拟合结果和实际收敛阶p判断误差估计过程是否可靠，并给出安全因子Fs取值。

5) 给出考虑数值不确定度的模型真解的95%置信区间：

[fi-Δ,fi+Δ]

(18)

当σ<Δ时，拟合误差可接受，此时定义

Δ=Fs|fi-f0|+σ+|fi-fi,fit|

(19)

当σ>Δ时，拟合误差较大，此时定义

(20)

虽然步骤5)可以针对网格序列中的所有网格离散解，但一般在细网格上进行，即模型真解的置信区间为

[f1-Δ,f1+Δ]

(21)

在实际应用中，直觉上应该是密网格的计算结果比粗网格更可信。因此在求解的最小二乘问题中考虑网格尺度的权重wi，修改为

(22)

3 考虑数值离散误差的模型选择不确定度量化方法

贝叶斯模型平均方法给出了每个模型输出为定值情况下的模型后验概率；数值离散误差估计给出了模型真实解的置信区间。模型真实解是不确定的，在这个区间内都有可能。此时模型的后验概率和输出的统计特性都是不确定的。

为考虑数值不确定度对模型概率更新的影响，借鉴认知随机混合不确定性因素量化的二阶概率框架[12]。该方法将认知和随机不确定性因素分离开，通过嵌套的迭代方法构建概率盒。其内层循环是随机不确定性因素，外层循环是认知不确定性因素。当固定认知不确定性因素为某一确定取值时，只有随机不确定性因素发挥作用，可以通过经典的概率方法量化关注输出量的随机响应，获得其概率密度函数和累计分布函数。当认知不确定性因素在给定区间内变化时，可以得到若干条可能的累积分布曲线。取该曲线簇的上、下限即为目标量的概率盒。

提出考虑数值离散误差的模型选择不确定度量化方法，具体流程如下：

1) 针对备选模型集合中的每个模型估计其数值离散误差，得到各自的置信区间。

2) 当每个模型在其置信区间内循环取值时，根据贝叶斯模型平均方法得到累积分布曲线簇，如图1中样本1～样本5所示。

3) 通过简单的代数运算得到目标变量累积分布函数的上、下限，如图1所示，即概率分布边界，将其定义为目标变量的概率盒。

图1 NACA0012翼型累积分布曲线簇和分布边界Fig.1 Clusters and bounds of cumulative distribution functions for NACA0012 airfoil

4) 通过概率盒可以分析一定置信水平下目标变量的置信区间，此时区间的端点值不是固定值，而是一个区间值。

4 计算结果与分析

在CFD工程应用中，通常使用湍流模型模化雷诺应力对平均流动的影响。学者们从涡黏性假设或雷诺应力输运方程出发，发展了数十种的湍流模型；如果加上各种湍流模型的修正形式，总共有上百种湍流模型。但现阶段没有对所有问题普适的湍流模型。许多CFD使用者还是基于以往的经验选择某种湍流模型完成计算。

针对湍流模型选择问题，通过NACA0012翼型低速绕流和CHN-T1翼身组合体跨声速绕流两个算例证明第3节提到的模型选择不确定度量化框架。需要指出的是，量化模型选择不确定度时，需要计算所有可能的备选模型。出于计算可行性考虑，选用两种使用最为广泛的湍流模型，即SA(Spalart-Allmaras)[13]和Menterk-ωSST(Shear Stress Transfer)[14]模型。

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

4.1 NACA0012翼型低速绕流

第1个算例是NACA0012二维翼型绕流问题，计算状态为Ma=0.15、临界雷诺数Rec=6.0×106，攻角为-3.99°～17.13°，试验结果可参考文献[20]。计算使用了4套网格，均为四边形单元。最密的一套网格使用了917 504个单元，C型网格拓扑，流向和径向单元数分别是1 792 和512，标示为Grid0。在密网格基础上通过相邻4个单元聚合形成下一级粗网格，网格尺度加倍。最粗的网格标识为Grid3。

4套网格计算的升、阻力系数如图2所示，直观上感觉最粗网格计算结果与其他网格计算结果差异较大，特别是在中到大攻角下。在进行数值离散误差估计时，由于要同时协调所有网格计算结果完成数据拟合，最粗网格的结果导致产生非常大的置信区间，失去了研究问题的意义。因此该算例使用较密的3套网格进行后续分析。此时通过第2节的流程方法可以得到升、阻力系数的置信区间，如图3所示。在小攻角时，置信区间很窄，说明此时数值不确定度较小；在大攻角时，置信区间变宽，说明此时数值离散引入的不确定度较大，离散格式或网格会对模拟结果产生显著影响。

图2 NACA0012翼型不同网格下SA模型升力和阻力预测结果Fig.2 Lift and drag prediction results with SA model using different grids for NACA0012 airfoil

图3 NACA0012翼型SA和SST模型升力和阻力预测的置信区间Fig.3 Confidence interval of lift and drag prediction with SA and SST models for NACA0012 airfoil

两种湍流模型计算产生比较显著的差异，如图4所示。对于升力系数预测，直观上无法判断哪种模型更接近试验结果；对于阻力系数预测，SA模型结果更接近试验值，特别是在中到大攻角下。

图4 NACA0012翼型密网格下SA和SST模型升力和阻力预测结果Fig.4 Lift and drag prediction results with SA and SST models using finest grid for NACA0012 airfoil

使用贝叶斯模型平均方法可以定量地给出两种湍流模型的后验概率，如图5所示。对于升力系数预测，大多数攻角下都是SST概率更大；对于阻力系数，大多数攻角下SA概率更大。

值得注意的是，两个模型的后验概率在负攻角和对应的正攻角区域不尽一致。这是因为贝叶斯模型平均方法与模型预测结果和试验结果密切相关，数据的偏差会显著影响后验概率。表1给出了在攻角-1.98°和2.00°的情况下使用最密网格时升力系数的计算结果和试验结果。此时试验结果或计算结果的微小误差都会影响模型概率，不过两种模型的计算结果差别不大，概率的偏差也不会过多影响输出的统计特性。

表1 计算和试验升力预测对比

从图5可以看出，网格对模型后验概率计算还是存在一定的影响。后验概率的差异必然会导致输出统计性质的不同。图6给出了不同网格下，考虑模型选择和模型预测偏差不确定度给出的关注输出量的累积分布函数，可以估算出阻力系数一定置信水平下的置信区间，如在攻角为11.13°时，根据密网格结果估计阻力系数的95%置信区间为[0.011 04, 0.016 61]，根据粗网格结果估计的置信区间为[0.010 52, 0.018 51]。可看出网格因素对置信区间估计影响显著。

图5 NACA0012翼型不同网格下SA模型升力和阻力预测的后验概率Fig.5 Posterior probabilities of lift and drag prediction using SA model with different grids for NACA0012 airfoil

图6 NACA0012翼型不同网格得到的累积分布曲线Fig.6 Cumulative distribution functions using different grids for NACA0012 airfoil

如要考查模型本身的后验概率，抛掉数值离散的影响，需使用第3节介绍的方法。此时，SA和SST模型的真实解无法确定，只能假设其在置信区间内都有可能，因此模型的后验概率也在一定的区间内都有可能，如图7中的灰色区域所示。

图7 NACA0012翼型考虑数值离散误差的SA模型升力和阻力预测的后验概率Fig.7 Posterior probabilities of lift and drag prediction using SA model considering discretization errors for NACA0012 airfoil

当只能确定模型后验概率的区间时，输出的统计特性也呈现一定的不确定度，可以通过概率的上、下限(如图8所示)给出置信区间端值的区间。如在攻角为13.31°时，升力系数的95%置信区间的左端值为区间[1.335 5, 1.344 0]，右端值为区间[1.407 0, 1.439 5]；在攻角为11.13°时，阻力系数的95%置信区间的左端值为区间[0.010 92, 0.011 25]，右端值为区间[0.015 86, 0.017 35]。

图8 NACA0012翼型累积分布函数边界Fig.8 Bounds of cumulative distribution function for NACA0012 airfoil

4.2 CHN-T1跨声速绕流

第2个算例是CHN-T1翼身组合体跨声速绕流[21-22]。该外形是中国空气动力研究与发展中心自行设计的单通道运输机标模，采用机翼/机身/平尾/立尾构型。计算状态为Ma=0.78，Rec=3.3×106，攻角为2.0°～4.0°，间隔0.5°。计算使用官方提供的非结构混合网格，单元数目从6 518 485(标示为Grid3)提升至162 253 613(标示为Grid0)。具体网格信息可以参考文献[23]。

通过第2节给出的离散误差估计方法可得两种湍流模型各自的95%置信区间，如图9所示。置信区间给出了模型真实解的估计情况，与离散数值解有不同的含义。

图9 CHN-T1外形SA和SST模型升力和阻力预测的置信区间Fig.9 Confidence interval of lift and drag prediction with SA and SST models for CHN-T1 configuration

通过贝叶斯模型平均方法可计算两种模型的后验概率，如图10所示。考虑数值离散的不确定度之后，无法直接判断哪个模型的后验概率更大。

图10 CHN-T1外形考虑数值离散误差的SA模型升力和阻力预测的后验概率Fig.10 Posterior probabilities of lift and drag prediction using SA model considering discretization errors for CHN-T1 configuration

考虑数值离散、模型选择、模型预测偏差3种不确定性因素后，得到升、阻力系数累积分布函数的上、下限，如图11所示。此时可以估计其置信区间的范围。比如对于升力系数，在2°攻角下，其95%的置信区间的左端值为区间[0.370 90, 0.372 75]，右端值为区间[0.433 0, 0.461 4]。

图11 CHN-T1外形累积分布曲线簇和边界Fig.11 Clusters and bounds of cumulative distribution curves for CHN-T1 configuration

5 结 论

将数值离散、模型选择、模型预测偏差3种不确定性因素统一考虑，发展了考虑数值离散误差的贝叶斯模型平均方法，给出关注目标量累积分布函数的上、下限。该方法抛掉了传统贝叶斯模型平均方法使用的数值离散解，通过嵌套方法构建概率盒。

针对湍流模型选择问题，通过NACA0012翼型低速绕流和CHN-T1翼身组合体跨声速绕流两个算例给出了该方法实施的具体示例。

关于“已知试验结果，通过计算结果估算其置信区间的意义”问题，需要指出的是模型确认是确定模型在预期用途内是否能表征真实流体系统的过程。当存在高可信度试验结果时，可以比较计算与试验结果对模型进行评价；当不存在试验结果时，需要根据已有的模型确认结论进行推论。本文发展的方法可以在有试验结果的情况下，计算关注目标量的置信区间。研究结论可以推广到没有试验结果的物理过程或状态，通过有限的数值计算估算关注量的置信区间。这给工业部门合理使用CFD模拟结果提供了新的思路。