在重大复杂工程领域, 众多工程事故已经引发社会各界的关注。人们总是抱着一个良好的愿望, 希望能逼真建模, 研制高可信度的数值模拟软件, 精确模拟各类事故的发生过程，科学分析其产生的机理。因此, 数值模拟软件必须能真实反映实际系统的外部表征、内在特性及其动作行为等现象, 其数值模拟结果才能使人置信。但由于物理过程的复杂性及人们的认知缺陷, 物理建模与数值模拟始终存在误差和不确定性。为此, 要利用数值模拟结果分析机理和模拟过程, 必须量化数值模拟结果的不确定性, 评估建模与模拟的可信度。但长期以来, 物理建模与数值模拟软件的研制者重点关注的是软件自身的研制和对实际问题的顺利模拟, 即算下去的问题, 对软件自身的正确性、模拟结果的不确定性和建模与模拟的可信度, 即算多准的问题, 一直没有给予足够的重视。尤其是对于理论、实验无法解释或难于解决的科学问题, 或做实验耗资巨大的工程设计, 很想利用数值模拟真实、全过程, 全系统、反复精密地进行模拟的研究者持一种矛盾的心态：既想利用数值模拟这条快捷经济的重要研究途径, 又对程序的正确性和数值模拟提供的结果心存疑虑。为了给复杂装置设计提供高效可靠的数值模拟工具, 必须展开计算模型正确性和可信度的评估研究。

流体力学方程及求解方法是众多复杂工程多物理过程建模与模拟(modeling & simulation，M&S)中的最基本模型, 此计算模型的验证与确认(verification and validation, V&V)受到高度重视^[1-7]。W.L.Oberkampf等^[8]综述了机械工程领域中科学计算V&V的发展, 系统论述了科学计算中V&V的基本概念、原理及步骤和发展过程。王瑞利等^[5]综述了V&V的相关概念、术语、规范、可信度评估方法和应用等的研究进展, 就M&S的V&V的几个关键问题进行了概括, 构建了复杂工程M&S的V&V的知识指南。确认是建立在定量准确评价的概念上, 是评价模型描述实际物理过程/实验(试验)的准确程度、模型形式与参数的适应范围, 确认刻画实际问题/实验(试验)区域的形状, 预测可信度评估等重要活动。对于流体力学确认的核心是随机输入参数对模拟输出响应量不确定传播进行量化的方法。

不确定度量化(uncertainty quantification, UQ)方法是利用各种数学和统计学方法, 对随机因素或认知缺陷尽可能给出分布或区间描述, 建立建模和模拟中随机或认知不确定性因素和建模与模拟描述系统响应量之间的关系, 进而量化影响系统性能的不确定性, 给出系统响应量性能分析的不确定度度量的方法。对于不确定度量化有多种方法, 其中敏感度分析方法、抽样方法、代理模型方法、随机微分方程方法和参数率定与优化方法已在复杂工程中得到很好的应用。

在上述几种方法中, 随机微分方程是20世纪中叶发展起来的一种不确定度量化方法。广义的随机微分方程包括随机过程驱动的微分系统和系数为随机量的微分方程, 狭义上的随机微分方程往往专指前者。对于随机过程驱动的微分系统, 方程的解往往是一个随机过程函数。虽然求解常微分方程的方法不适用于求解随机微分方程, 然而, 有些复杂过程可以直接使随机过程建模、模型的不确定性, 直接表现在模型本身上而不是输入参数的不确定性, 随机微分方程是研究不确定度量化的很好途径。另一方面, 系数为随机量的微分方程可以直接处理为随机参数的微分方程, 常见的分析方法是随机谱方法。

随机谱方法是使用谱逼近随机微分方程或参量, 然后将其分解为独立的确定性分量和随机分量来量化不确定性。其数学思想是数学模型中的每一个参量利用正交多项式如(如Hermite多项式)展开成无穷级数项, 在实际应用中取有限项。无穷级数第1项表示参量的期望, 第2项表示高斯随机波动, 第3项及高阶项表示非高斯随机波动。随机谱方法已在建模和模拟不确定度量化中得到了广泛应用。常见的2种方法:KL展开(Karhunen-Loeve expansion)和多项式混沌方法(polynomial chaos, PC)。本文中，将着重介绍PC在Riemann问题不确定度量化中的应用。

1.   流体力学方程组与迎风格式

一维拉氏流体力学方程组^[9]:

式中：ρ、u、E、e和p分别为流体密度、流体速度、总能、比内能和流体压强, v为比容, γ为理想气体多方指数。一维拉氏流体力学方程组迎风型格式见文献[9]

2.   多项式混沌法

多项式混沌法其数学思想是数学模型中的每一个参量利用正交多项式(如Hermite多项式)展开成无穷级数项。设系统的解具有如下形式：

式中：${\mathit{\boldsymbol{\widetilde U}}_i}(x, t) = \left( {{{\tilde v}_i}, {\rm{ }}{{\tilde u}_i}, {{\tilde E}_i}} \right)$为确定性向量函数, ψ_i(ξ)为随机变量ξ的多项式函数, 且{ψ_i(ξ)}_i=0^x为正交多项式族, ξ是概率空间(Ω, F, P)上的随机变量, U(x, t, ξ)为概率空间(Ω, F, P)上的随机函数, Ω为样本空间, F是Ω上的σ代数, P:F→[0, 1]是概率测度, θ为事件, θ∈Ω。不同的随机变量类型对应不同的多项式ψ_i^[10]。由于炸药颗粒服从正态分布, 因此本文中选取Hermite多项式。在实际应用过程中, 使用有限项展开, 设：

式中：${N_{{\rm{PC}}}} = \frac{{(d + n)!}}{{d!n!}}-1$为阶段长度，n为随机变量个数, d为多项式阶数。本文中, n=1, 则N_PC=d。

从式(2)~(3)可以看出, 只要用某种方式得到了级数中的确定性向量函数$\mathit{\boldsymbol{\tilde U}}$_i(x, t), 则可得到随机函数U(x, t, ξ(θ))或U(x, t, ξ(θ))的近似值。根据求解$\mathit{\boldsymbol{\tilde U}}$_i(x, t)方法的不同, 将PC分为嵌入多项式混沌(intrusive polynomial chaos, IPC)和非嵌入多项式混沌(non-intrusive polynomial chaos, NIPC)。

IPC是将数学模型中的每个参量进行正交多项式展开, 然后利用Galerkin映射, 得到相应的随机控制方程, 但是每一次展开都需要重新书写控制方程。因此，IPC无法利用已有程序, 需要对已有的数值求解程序进行大量修改或者重新研制计算程序。

NIPC是把已有的数值求解程序作为一个黑匣子, 在随机空间里通过一定的抽样方法, 获得若干个样本点, 将样本点输入到确定性程序求解, 然后对确定性输出结果进行统计分析, 以获得相关数值求解结果的统计特征, 来评估输入参数或计算条件的不确定性在计算过程中传播的影响。NIPC采用已有的数值模拟程序, 不需要对控制方程进行修改, 不需要重新编写程序。NIPC不确定度评估的计算流程见文献[7]。

3.   P-box的构造方法

若问题中同时存在偶然不确定和认知不确定, 需要将其结合起来考虑, 此时可采用概率盒方法进行描述，即将各参数分为内外2层循环：外层循环考虑认知不确定参数, 使用拉丁超立方取样(Latin hypercube sampling, LHS)方法抽样；内层循环考虑偶然不确定参数, 取值由随机抽样给出。

P-box计算步骤:

(1) 外层循环考虑认知不确定参数, 从其区间中, 采用LHS方法得到N个样本；内层循环考虑偶然不确定参数, 利用Monte Carlo(MC)抽样得到M个样本, 作为左端密度初值;

(2) 固定一个外层认知不确定变量γ, 计算出相应的M个结果;

(3) 寻找相应物理量的最大、最小值, 然后分成100份, 计算每一份对应的样本容量, 绘出相应的累积分布曲线;

(4) 对剩余的N-1个认知不确定数施行同样运算, 划出N条累积分布曲线;

(5) 选取最内层和最外层的那一条作为概率边界，构成一个P-box;

(6) 读出概率盒的意义。

4.   混合不确定量化方法在Sod问题中的应用

Sod激波管问题是初始间断的分解问题, 在一个无限长管道中, 在x=0处的左端是高压气体, 右边是低压气体, 在初始时刻处于静止状态, 两者之间用薄膜隔离。t=0时刻薄膜破裂, 形成向左稀疏波和向右激波。Sod问题计算区域为[-1, 1], 初值为:

x≤0时，ρ_L~N(1.0, 2.5×10^-7), u_L=0.0, p_L=1.0, γ_L∈[1.4, 3.1];

x＞0时, ρ_R=0.125, u_R=0.0, p_R=0.1, γ_R=1.4。

采用理想气体状态方程:

p=(γ-1)ρe

此问题的特点是左右2种状态的γ不一样, 是一个多介质激波管问题, 有解析解。但其中左右密度量含有不确定因素, 对其量化才能实施计算模型的确认。对于左端炸药密度, 由于加工过程中炸药颗粒凝结的不均匀性, 导致密度存在随机性。

假设其服从正态分布, 期望为μ=1.4 g/cm³, 标准差σ=0.5 mg/cm³, 即炸药密度服从正态分布：

对炸药密度采用MC方法, 取样1 000次, 分布直方图如图 1所示。

图  1  炸药密度分布直方图

Figure  1.  Histogram showing distribution of explosive density

下载: 全尺寸图片 幻灯片

对于多方指数γ，不知道它的任何信息, 只知道它位于[1.4，3.1], 在实际设计或计算中采用哪个数值由设计者决定, 这种选取的过程与认识水平有关, 属于认知不确定度。利用LHS采样10次, 如图 2所示。

图  2  多方指数分布散点

Figure  2.  Scatter of polytrophic exponent

下载: 全尺寸图片 幻灯片

当固定γ时, 会得到1 000组状态确定性Sod问题的样本点, 将其代入迎风格式的流体力学程序和精确解法器, 得到1 000组响应量的样本状态, γ=2.29时, 计算到t=0.22时, 密度、压力、速度、能量的状态构成的马尾图, 如图 3所示。从图 3可以看出, 在某些地方解的分散度很明显。当γ变化时, 会得到10张类似的马尾图。

图  3  Sod问题1000组计算结果

Figure  3.  1000 sets of computational results for the Sod problem

下载: 全尺寸图片 幻灯片

4.1   可信范围

对于上述计算结果, 信息量巨大, 现在提取有用的信息。首先, 将10张马尾图合并, 得到10×1 000=10 000组响应量样本状态, 采用NIPC方法开展不确定度量化。

图 4给出了计算得到的密度、压力、速度、能量的期望和方差, 从图可以看出有些地方方差很大, 精确求解黎曼解法器的方差比迎风格式的高。

图  4  t=0.22时Sod问题的期望和方差

Figure  4.  Expectation and variance for the Sod problem at t=0.22

下载: 全尺寸图片 幻灯片

图 5给出了迎风格式可信范围和精确黎曼解法器可信范围, 从图可以看出有80%的地方完全重叠, 其他局部重叠或者不重叠, 说明采用迎风格式求解Sod问题是可信的。

图  5  2种方法计算Sod可信区间的比较

Figure  5.  Comparison of confidential intervals from two computational methods in solving the Sod problem

下载: 全尺寸图片 幻灯片

4.2   Sod问题的P-box

利用第3节中的方法, 构造Sod问题中迎风格式的P-box, 做出在x=0.1, t=0.22处, 密度函数的累积分布函数, 如图 6所示, 然后选取最内层和最外层概率边界曲线构成一个P-box, 如图 7所示。

图  6  密度函数的累积分布函数族

Figure  6.  Cumulative distribution function ensemble of density function

下载: 全尺寸图片 幻灯片

图  7  利用迎风格式构造的密度函数的P-box

Figure  7.  P-box constructed by the upwind scheme for density function

下载: 全尺寸图片 幻灯片

从图 7可以看出，在t=0.22, x=0.1时, 密度ρ(x, t)≤0.555 1的概率0.1＜P＜0.6, 0.95置信区间位于[0.590, 0.609], 并且在不确定初边值条件下, ρ(x, t)≤0.54的概率为0, ρ(x, t)＜0.61依概率1成立, 为不确定评估找到了依据。

同时, 做出黎曼解法器的P-box, 并与迎风格式做比较, 如图 8所示, 不难看出迎风格式构成的P-box基本包含黎曼解法器的P-box, 两者差异不大, 这说明采用迎风格式求解Sod问题是有效的。

图  8  用2种方法构成的密度函数的P-box

Figure  8.  P-box constructed by two methods for density function

下载: 全尺寸图片 幻灯片

5.   结论

(1) 流体力学数值模拟过程中初始状态或输入条件常含有偶然和认知混合型的不确定性，对模拟结果的响应量有很大影响。

(2) 针对流体力学初始状态中含有偶然不确定性，给出了非嵌入式多项式混沌法量化方法。针对初始状态中含有认知不确定度，给出了P-box量化方法，提出了一种流体力学输入参数对模拟输出响应量不确定性传播量化的方法。

(3) 以迎风格式和黎曼解法器求解Sod问题为例，给出了Sod问题左状态密度ρ_L偶然不确定度和多方指数γ_L认知不确定度，分别采用非嵌入式多项式混沌法和P-box量化结果，给出了流体力学方程组迎风格式数值求解可信度与输入不确定度的范围关系，为流体力学数值求解过程中不确定性传播量化和可信度评估提供了一种有效方法。