LAD2D程序^[1]是一个非结构网格拉氏自适应的二维爆轰流体动力学软件, 此程序已通过软件质量保证(SQA)和大量测试模型的考核, 验证了程序的正确性^[2-4], 但由于爆轰流体力学物理过程的复杂性和人们认识的缺陷, 在建模过程中含有抽象、简化和近似, 逼真建模很难, 有时只能唯象建模和逐渐逼近真实情况, 建模过程含有不确定性。加之描述其过程的数学物理模型是高度非线性的偏微分方程组, 很难解析求解^[5-6]。在数值求解过程中, 由于连续到离散, 存在计算模型误差、离散误差、计算机舍入误差和分析误差等, 数值模拟过程始终是一种近似, 含有不确定性。为此, 数值模拟结果的置信度一直缺乏科学的论述^[7-8], 不确定度量化是置信度评估的核心, 不确定度量化(uncertainty quantification, UQ)方法分为概率框架和非概率框架下的多种量化方法^[9]。本文中采用概率框架下的误差累积分布函数(cumulative distribution function, CDF)和互补累积分布函数(complementary cumulative distribution function, CCDF), 进行统计分析, 给出了CDF和CCDF曲线, 形成爆轰计算中输入不确定性参数、网格尺度与输出响应量的误差马尾状图, 从而可量化分析爆轰模拟结果的置信度, 对多物理爆轰过程发展高置信度数值模拟软件有很好的借鉴作用。

1.   数学模型

1.1   CDF方法及马尾图

假定在模型计算区间[a, b]上有一组标准解：X=x₁, x₂, …, x_n, Y=y₁, y₂, …, y_n。输入参数为λ, 通过数值计算得到一组数值解:$\tilde X = \left( {{{\tilde x}_1}, {{\tilde x}_2}, \cdots , {{\tilde x}_\mathit{m}}} \right), \mathit{\tilde Y = }\left( {{{\tilde y}_1}, {{\tilde y}_2}, \cdots , {{\tilde y}_\mathit{m}}} \right)$。CDF方法的步骤为：

(1) 计算数值点上的标准解。在区间a, b上, 由标准解插值得到$\tilde X = \left( {{{\tilde x}_1}, {{\tilde x}_2}, \cdots , {{\tilde x}_\mathit{m}}} \right)$点上的标准解$Y' = \left( {{{y'}_1}, {{y'}_2}, \cdots , {{y'}_m}} \right)$：

这样通过式(1), 就可得到$\tilde X = \left( {{{\tilde x}_1}, {{\tilde x}_2}, \cdots , {{\tilde x}_\mathit{m}}} \right)$对应的标准解$Y' = \left( {{{y'}_1}, {{y'}_2}, \cdots , {{y'}_m}} \right)$。

(2) 计算误差。由${\tilde X}$的标准解$Y' = \left( {{{y'}_1}, {{y'}_2}, \cdots , {{y'}_m}} \right)$和数值解$\mathit{\tilde Y = }\left( {{{\tilde y}_1}, {{\tilde y}_2}, \cdots , {{\tilde y}_\mathit{m}}} \right)$, 就可以得到计算误差。即

(3) 统计分析。假设误差分为K段, 统计小于每段误差的个数, 然后计算出平均值, 即概率:

由式(3)就可以得到马尾图, 建立输入参数λ对输出响应量误差的关系图, 以量化输入参数λ对输出响应量的置信度, 判断输入参数计算条件的合理性。

1.2   爆轰流体力学模型

数值模拟使用的基本方程是不定常可压缩理想流体力学方程和化学动力学方程的藕合方程组:

式中:E=e+12u², 其中ρ、u、E、e、P分别表示密度、速度、单位质量的总能量、单位质量的内能与压力。方程组(4)~(7)的主要问题是状态方程(7)与一般的流体力学计算不同, 其中包含了燃烧函数F, 燃烧函数F要能正确反映化学反应的特性。爆炸产物的状态方程可以是JWL形式的状态方程, 也可以是理想气体状态方程P=(γ-1)ρe, 这里选用JWL形式状态方程：

式中:$V = \frac{v}{{{v_0}}}$, E为比热容力学能, A、B、R₁、R₂和w是常数, v与v₀是比容。

呈静止状态的凝固炸药P=0。为了进行数值计算, 用一条光滑曲线将它们连接起来, 通常引进所谓的燃烧函数F来表征炸药反应程度, 这里考虑如下燃烧函数:Wilkins函数(时间燃烧函数+C-J比容燃烧函数)：

式中：V_J=γV₀/(γ+1)是C-J比容, V₀=1/ρ₀, t_b是爆轰波刚到达计算网格的时刻(开始燃烧), 即起爆时间^[10], t是当前计算时刻, ΔL=r_bΔR/D_J, ΔR是网格宽度, D_J是C-J爆轰速度, n_b、r_b是可调参数。计算时, 考虑人为黏性。

(1) von Neumann-Richtmyer人为黏性(二次黏性):

式中:l_NR为具有长度量纲的量, $l_{{\rm{NR}}}^2 = a_{{\rm{NR}}}^2A$；a_NR为N-R人为黏性系数；A为计算网格面积。

(2) Landshoff人为黏性(一次黏性或线性黏性):

式中:l_L为具有长度量纲的量, ${l_{\rm{L}}} = {a_\rm{L}}\sqrt A$, a_L为Landshoff人为黏性系数；A为计算网格面积；c为当地声速。

1.3   CDF法与爆轰流体力学数值模拟结果置信度量化

1.2节中爆轰模型计算涉及很多输入参数, 包括产物JWL状态方程中A、B、R₁、R₂、w, 反应程度及燃烧函数中n_b、r_b和流体计算时人工黏性系数a_NR、a_L这些输入参数对输出计算结果有多大影响, 需要进行分析评估。

采用CDF法时, 首先凭经验给出一套计算参数, 针对其中一个参数进行扰动, 其他参数固定不变, 通过确定性程序, 进行系列模拟计算, 得到数值解。其次根据基准解, 可以是解析解, 也可以是高精度数值解, 计算数值解与基准解的误差。然后形成误差的CDF及马尾图, 以量化数值模拟结果的置信度。

2.   计算结果与分析

2.1   一维黎曼问题

SOD问题常称为激波管问题。许多间断解方法的设计和构造, 都利用这种激波管问题进行可靠性和准确度的数值实验, 从而判断和检验方法、格式和程序的优劣。SOD问题是一个非常柔和的例子, 精确解包含了一个左稀疏波、接触间断和一个激波。下面采用SOD问题作为数值例子测试考核LAD2D程序的正确性和数值格式的精度。

SOD问题一般计算区域为[-1, 1], 初值为:

式中：ρ_L、U_L、P_L分别指左端初始状态的密度、速度与压力, ρ_R、U_R、P_R分别指右端初始状态的密度、速度与压力。左右边界采用连续边界条件, 状态方程采用理想气体:P=(γ-1)ρe, γ=1.4。计算采用CFL=0.5, 一次黏性系数取0.06, 二次黏性系数取0.2。图 1给出了20~800个网格单元间隔10的16套网格计算到t=0.5时计算结果与解析解的比较。图 2给出了16套网格计算到t=0.5时计算结果与解析解的比较误差的累积分布函数CDF的统计分析马尾图。从图中可以清楚地看出计算误差与网格尺度的关系。

图  1  网格尺度0.1~0.025(16套)时计算结果与解析解的对比

Figure  1.  Comparison between numerical results and analytic solutionsfor different mesh scales of 0.1~0.025(16 sets mesh)

下载: 全尺寸图片 幻灯片

图  2  网格尺度0.1~0.025(16套)计算结果与解析解的对比误差统计分析的马尾图

Figure  2.  Horsetail figure of error between numerical results and analytic solutions for different mesh scales of 0.1~0.025(16 sets mesh)

下载: 全尺寸图片 幻灯片

从图 2可以清楚地查出LAD2D程序计算SOD问题的置信度情况。假设要求计算误差在0.01(1%), 那么LAD2D程序计算SOD问题只有在网格尺度为0.01(200个网格单元)下, 置信度才能超过90%。假设要求计算误差在0.02(2%), 那么LAD2D程序计算SOD问题只有在网格尺度为0.025(80个网格单元)下, 置信度就能超过90%。

爆炸波问题也是一个一维黎曼问题, 由于有一边状态的压力达到爆轰的压力, 所以称为爆炸波问题。同SOD问题相比, 此问题在激波处有一个很窄的强激波区, 对程序格式健壮性考核有重要意义。爆炸波问题一般计算区域为[-1, 1], 初值为:

式中：ρ_L、U_L、P_L、γ_L分别指左端状态初始的密度、速度、压力和理想气体状态方程系数, ρ_R、U_R、P_R、γ_R分别指右端状态初始的密度、速度、压力和理想气体状态方程系数。左右边界采用连续边界条件, 状态方程采用理想气体:P=(γ-1)ρe。计算采用CFL=0.5, 一次黏性系数取0.06, 二次黏性系数取0.2。图 3给出了14套网格计算到t=0.1时计算结果与解析解的比较误差的累积分布函数CDF的统计分析马尾图。从图中可以清楚地看出计算误差与网格尺度的关系。假设要求计算误差在0.01(1%), 那么LAD2D程序计算爆炸波问题只有在网格尺度为0.003 3(600个网格单元)下, 置信度才能超过90%。假设要求计算误差在0.02(2%), 那么LAD2D程序计算爆炸波问题只有在网格尺度为0.005(400个网格单元)下, 置信度能超过90%。从SOD和爆炸波的统计分析看, 网格计算尺度在0.04~0.025之间, 计算误差如果要求在0.01~0.05之间时, SOD问题的置信度为0.6~0.99, 爆炸波问题的置信度为0.4~0.7, 置信范围太大, 很难把握, 需缩小置信范围。

图  3  网格尺度0.04~0.025(14套)时计算结果与解析解的对比误差统计分析的马尾图

Figure  3.  Horsetail figure of error between numerical results and analytic solutions for different mesh scales of 0.1~0.025(14 sets mesh)

下载: 全尺寸图片 幻灯片

2.2   平面爆轰问题

众所周知, 爆轰波达到定态时, 数值计算一定要符合Chapman-Jouguet理论。也就是说, 在声速点的波后流场与CJ模型相一致, 即可达到CJ状态。炸药取PBX-9404炸药, 参数为:K＝2.827、D_J=8.88 km/s, ρ₀=1.842 g/cm³, 由CJ理论可以解析推出PBX-9404炸药CJ状态^[11]为:p_CJ=37.27 GPa。计算区域为[0, 10], 初始在左边起爆。此问题起爆后经过一段时间可以达到稳定的爆轰波及CJ状态, 以检验程序能否达到CJ状态, 是爆轰数值模拟软件适应性考核的最基本问题。

(1) 网格尺度置信度分析

在爆轰计算条件(主要是输入参数)固定的情况下, 改变网格尺度, 统计分析程序的置信度。起爆采用压缩比σ=1.03, 燃烧函数采用Wilkins函数, 参数n_b=1.0, r_b=2.1。状态方程采用JWL形式, 输入参数A=765.788 GPa、B=14.249 GPa、R₁=4.3、R₂=1.45、ω=0.28。计算CFL=0.7, 一次黏性系数为0.06, 二次黏性系数为2.0。图 4给出了1~10 μs时10套(网格尺度分别为0.1、0.09、0.08、0.07、0.06、0.05、0.04、0.03、0.02、0.01)网格计算结果不同时刻与基准解(选取10 000个网格的数值解作为基准解)对比误差的累积分布函数CDF。

图  4  网格尺度0.1~0.01(10套)时不同时刻计算结果误差统计分析的马尾图

Figure  4.  Horsetail figure of error between numerical results and analytic solutions at different times

下载: 全尺寸图片 幻灯片

针对一维平面爆轰, 假设要求计算误差在0.01(1%), 从图 4中的t=1 μs时计算结果与网格尺度变化来看, 网格尺度从0.1~0.01之间时, 置信度都可超过90%。但从t=4、8 μs的计算结果可以看出, 假设要求计算误差在0.01(1%), 网格尺度在0.1~0.01之间时, t=4 μs置信度为0.6~0.8, t=8 μs仅0.3~0.4。从这个分析可以看出, 随着计算时间的推进, 置信度逐渐下降, 甚至有的置信度很差。说明目前爆轰计算的建模很不适应, 需要大大改进计算模型。但无论从t=1、4 μs, 还是从t=8 μs计算结果看, 网格尺度达到0.01时, 置信度大于0.6, 所以在目前建模情况下, 爆轰计算的网格尺度要尽量小, 网格尺度最好小于0.01, 甚至更小, 才能大大提高爆轰计算的置信度。

(2) JWL状态方程中输入参数敏感性分析

爆轰计算中采用JWL状态方程式(8)时, 其参数A、B、R₁、R₂和w敏感性分析对爆轰模型的使用有重要意义, 其参数取值大小对实际问题模拟结果的影响是爆轰模型及模拟结果用于决策的重要依据。针对爆轰产物JWL状态方程参数随机选取了4套(见表 1)参数, 其他输入参数同(1)的计算条件, 在网格尺度分别为0.1和0.01情况下进行了数值模拟, 将网格尺度为0.001、JWL状态方程参数取第一组时的计算结果作为基准解。图 5为4套参数计算结果的误差统计的累积分布函数CDF, 其中基准解参数采用PAP1, 计算分点采用10 000个点。

表  1  JWL状态方程4套参数

Table  1.  Four sets of parameters for JWL state eqution

参数代号	A/GPa	B/GPa	R1	R2	ω
PAP1	765.788	14.249	4.3	1.45	0.28
PAP2	947.774	59.206	5.0	2.00	0.28
PAP2	1 502.851	134.913	6.3	2.40	0.28
PAP3	3 250.433	233.271	8.2	2.80	0.28

下载: 导出CSV 

| 显示表格

图  5  4套参数在不同网格尺度下的计算结果误差统计分析的马尾图(t=8 μs)

Figure  5.  Horsetail figure of error with four sets parameters for different mesh scales(t=8 μs)

下载: 全尺寸图片 幻灯片

假设基准解选取精度很高时, 从图 5可以看出, 当误差要求为0.01(1%)时, 4套参数对计算结果影响很大, 在网格尺度为0.01时, 密度分布置信度从0.9下降到0.3, 压力分布置信度从0.45下降到0.28。当误差要求为0.04(4%)时, 4套参数对计算结果密度影响不大, 置信度均大于0.9, 压力仍影响很大, 置信度仍从0.45下降到0.28。从这个计算结果说明, 选取爆轰产物JWL状态方程参数应该引起重视。但选取哪一套参数合理, 置信度高, 需要建立合理、精度高的基准解, 才能应用本文方法合理选取参数。

3.   结论

(1) 通过标准解与数值解之间的误差累积分布函数(CDF)及马尾图表征方法, 建立数值计算输入参数与网格尺度对输出响应量影响的置信度的表征方法。

(2) 基于自主开发的爆轰弹塑性流体力学软件LAD2D, 针对一维黎曼问题, 给出了网格尺度变化下数值模拟结果与解析解之间误差的马尾表征图。对于经典的SOD问题, 假设要求模拟误差达到0.01(1%), 那么计算网格尺度必须小于0.01, 置信度才能超过90%。假设要求计算误差达到0.02(2%), 即降低要求时, 计算网格尺度只要小于0.025, 置信度就会超过90%。对于强激波的爆炸波问题, 若要求计算误差达到0.01(1%), 那么计算网格尺度必须小于0.003 3, 置信度才能超过90%。假设要求计算误差达到0.02(2%), 即降低要求时, 计算网格尺度只要小于0.005, 置信度就会超过90%。但从爆炸波与SOD问题相比, 计算网格尺度与激波强弱关系很大, 激波越强, 需要计算网格尺度越小, 这与理论分析是一致的。

(3) 针对炸药爆轰产物JWL状态方程参数取值对计算结果的影响, 给出了4套参数在两种网格尺度下误差的马尾表征图。当误差要达到0.01(1%)时, 在网格尺度为0.01时, 密度分布置信度从0.9下降到0.3, 压力分布置信度从0.45下降到0.28。当误差要求达到0.04(4%)时, 即降低要求时, 密度分布置信度均大于0.9, 压力分布置信度仍从0.45下降到0.28。从这个分析结果说明, 选取爆轰产物JWL状态方程参数应该引起重视。

(4) 研究结果表明, 基于误差马尾图量化参数敏感性和置信度的思想是可行的, 为非线性、强间断、多物理耦合问题的数值模拟结果置信度评估提供了一种行之有效的方法。