爆轰是极其剧烈的物理化学反应，炸药爆炸过程扩展的速度高达每秒数千米到万米之间，所形成的温度约为3 000～5 000 ℃，压力高达10²～10⁴ MPa，发生在极短的时间（10⁻⁹ s量级）和极窄的空间（10⁻⁴ m 量级），释放极高的能量（10² GW/cm²量级）^[1-5]。爆轰过程的极端特征给试验精确测量和理论表征带来了极大挑战。由于爆轰试验成本高和化学反应的不稳定性，只有经过多次试验才能标定感兴趣量（quantity of interest, QoI）的取值。即使是常规炸药，也只有有限的试验数据。与试验相比，建模与模拟（modeling and simulation, M&S）是一种经济、安全的途径。为降低运算成本，提高计算可行性，通常还使用假设、简化和近似等手段，以降低模型的复杂性^{[1, 4-6]}。Chapman-Jouguet（C-J）理论假设流动是一维的，不考虑热传导、热辐射以及黏滞摩擦等耗散效应，忽略冲击起爆过渡和不定常效应，将复杂非线性多物理爆轰过程简化为一个具有冲击压缩间断面的一维问题。C-J理论能建立波后QoIs和爆速、初始密度以及唯象参数的函数关系式，能预测爆压和爆热等不易测量或测量成本较高的物理量，具有实用性强、适用范围广的优点，是高能炸药爆轰研究的有效工具^[5-6]。

C-J理论是学术界公认的成功理论，其模拟过程中使用的物理量和参数一直被认为是确定的。事实上，由于炸药组分的异质性以及物理量测量的随机特征，物理量和唯象参数都会受到不确定性的干扰。例如，爆速是描述炸药力学行为的基本物理量，某些HMX基高能炸药的爆速超过9×10³ m/s，速度的瞬时极大变化率导致缺乏精密的试验装置进行测试。传统的试验装置如圆筒试验，所用金属的纯度和密度不均匀性、加工公差的存在以及试验数据的拟合误差，都会导致测量结果的不确定性。同时，含能材料组分的细微变化、几何边界的不精确、极端荷载下金属延展性变化以及黏结剂的存在也会导致试验结果的波动^[7-9]。

密度是决定爆轰性能的另一个基本物理量。密度的不确定度来自于炸药加工过程中产生的空腔、孔隙、裂纹、扭结和颗粒结晶的随机性以及杂质的混入^[10]。此外，存储温度的变化也会引起密度的波动^[11]。密度的微小变化会激发感度和热点的变化，甚至会引起爆轰系统输出结果剧烈的、奇异的甚至目前理论无法解释的变动^{[7, 12]}。

综上所述，爆轰实际过程中，如爆速和初始密度等物理量都会受到噪声扰动。这类不确定度是炸药的固有属性，即使增加样本容量，也很难减少或降低这种不确定性。另一方面，爆轰产物状态方程、反应率函数等爆轰唯象模型中的参数，无法通过试验直接标定，需要凭借工程技术人员的经验^{[3, 13-14]}，这类不确定性也会严重影响炸药爆轰数值模拟的精度和可信性。

事实上，工程技术人员和学者们都已经意识到爆轰中的不确定因素^{[1, 5-6, 9, 15]}。文献[1, 6]甚至将物理量的容许不确定度（误差）定为 $\pm 5$ %。然而，很少有学者专门研究不确定因素对QoI的影响。实际上，不确定因素的存在使得研究人员一直面临着试验成本太高和数值模拟不稳定的矛盾心态。导致研究人员既想利用C-J理论这种经济、安全的手段，又对C-J理论的使用范围和预测结果心存疑虑。因此，爆轰不确定度量化（uncertainty quantification, UQ）研究是增强炸药爆轰唯象模型的可靠性、标定模型的使用范围、缓和试验与数值模拟矛盾的重要途径，具有重要的学术意义和实用价值。

近10年来，欧美UQ技术在水文、地理、预报、经济、自动控制和结构力学分析等领域蓬勃发展。但爆轰UQ研究结果较少，面临众多挑战。首先，在刻画炸药爆轰模型中输入不确定度方面，数理统计理论完备，技术成熟，是量化不确定度的最自然的工具。根据大数定律和中心极限定理，正态分布是描述输入不确定度的有效工具，且参数易于通过经典假设检验法标定。然而，正态分布无法保证密度和爆速的非负性。其次，在如何准确得到炸药产物状态方程（equation of state, EOS）和反应率函数等唯象模型中参数的可信取值范围方面，通常根据工程技术人员的经验给定取值范围，存在很大的人为因素。若假设参数服从此区间上的、无信息的（un-informative）均匀分布，很多重要信息会丢失。独立同分布是数理统计中的重要假设，但爆轰系统的输入不确定度未必满足此假设。除此之外，针对如何评估输入不确定度对系统响应量的影响方面，即使通过大规模工程计算，也很难给出QoIs的概率密度函数（probability density function, PDF）的显式函数关系。幸运的是，计算技术的快速发展和数值分析方法的进步提供了计算PDF的有效方法。与Monte Carlo方法相比，多项式混沌（polynomial chaos, PC）收敛速度快、运算成本低，是目前工程技术领域处理不确定度传播的主流工具^{[3, 16-18]}。特别是非嵌入式PC，不需要更改程序，且能给出Gauss统计量的显式表达，是处理爆轰UQ的潜在工具。但是在爆轰UQ研究中，会面临输入不确定度概率类型不一致、随机变量不独立的问题，导致PC理论无法直接应用。

面对众多挑战，本文中，结合真实的试验数据，通过参数估计和Anderson-Darling 假设检验法标定物理量的PDF，采用Beta分布定量刻画没有物理意义的、唯象参数的不确定度，评估输入不确定度对QoIs的影响，以增强C-J模型的可靠性和预测精度，以期为武器型号设计、装甲防护和土木施工管理等提供决策依据。

1.   Chapman-Jouguet理论

假设爆轰波运动方向和波阵面垂直，忽略波阵面的厚度和化学反应时间，忽略热传导和黏性效应以及非定常效应，则爆轰过程可简化为图1所示的C-J模型。图1中，ρ₀、p₀、U₀和e₀分别代表未反应物的初始密度、压力、速度和内能，D为爆速。爆炸物以D–U₀的速度进入波阵面，以D–U₁的速度流出波阵面。

图  1  爆轰波在高能炸药中的传播示意图

Figure  1.  Sketch of propagation of detonation wave into high explosive

下载: 全尺寸图片 幻灯片

根据Hugoniot守恒定律和C-J假设，波阵面两侧物理量满足如下关系。

(1) 质量守恒：

(2) 动量守恒：

(3) 能量方程：

(4) Chapman-Jouguet假设：

(5) 声速定义：

式中： $v_{1}={1}/{\rho_{1}}$ ，为爆轰产物的比容； $v_{0}={1}/{\rho_{0}}$ ，为未反应炸药的比容； $\rho_{1}$ 为爆轰产物的密度； ${p_1}$ 为波后压力； U₁为波后介质速度；e₁为爆轰产物的内能；c为声速；S代表等熵状态。

对于高能炸药，p₀ $\ll$ p₁，因而假设p₀=0。由于爆轰波未到达时，图1右侧高能炸药近似认为处于静止状态，进而假设U₀=0。且高能炸药反应物和产物的EOS均使用γ律：

由式(1)～(6)导出波后介质QoIs的函数表达式：

式中：ρ_J、p_J、U_J和c_J分别代表C-J理论导出的爆轰产物的初始密度、压力、速度和声速。式(7)建立了C-J状态下波后介质QoIs和物理量D、ρ₀以及唯象参数γ的函数关系式。由于波后介质物理量试验标定困难且成本较高，因此式(7)成为预测波后介质QoIs的实用工具。

2.   爆轰不确定因素的定量刻画

受测量技术、含能材料自身固有属性和人类认知的局限性等因素的影响，爆轰M&S中的物理量和唯象参数均受到噪声污染^[12-14]。ρ₀的不确定度主要来源于压制成型时炸药颗粒的扭结、随机结晶、空洞和裂纹的存在以及杂质的掺入；爆速D的不确定度主要来源于测量技术。即D和ρ₀的不确定度与材料的固有属性有关，即使增加样本容量，不确定度也不会消失。

2.1   不确定物理量的表征

根据Wilkins等^[19]的理论，在爆轰计算中，合适的概率分布能合理描述含能材料密度的非均匀性。本文中，假设PBX-9502的初始密度服从对数正态分布 $\mathfrak{L}\left(\hat\mu,\, \hat{\sigma}^{2}\right)$ ，其中 $\hat{\mu}$ 和 $\hat{\sigma}$ 分别为对数期望和对数标准差。与正态分布相比，对数正态分布的优势在于符合统计规律的同时能保证物理量的非负性。首先，利用浸液法测量炸药的密度^[20]，给出PBX-9502初始密度的统计直方图，如图2所示。根据矩估计法导出试验结果的期望μ=1.894 g/cm³和标准差σ=0.002 5 g/cm³。进而，使用Anderson-Darling假设检验法，在5%显著性水平下，ρ₀接受服从对数正态分布的原假设。最后，通过曲线拟合技术，得到PBX-9502的概率密度函数。如图2所示，初始密度的概率密度函数曲线与试验数据吻合很好。

图  2  PBX-9502初始密度的概率密度函数

Figure  2.  PDF of initial density of PBX-9502

下载: 全尺寸图片 幻灯片

利用公式：

导出对数期望和对数标准差分别为 $\hat{\mu}$ =0.6108， $\hat{\sigma}$ =0.0043。

爆速是爆轰研究中唯一能使用较简单试验装置和流程即可标定的物理量，这里的试验结果来自于测时仪法^[20]。假设 $D \sim \mathfrak{L}\left(\hat{\mu}_{\rm{D}},\, \hat{\sigma}_{\rm{D}}^{2}\right)$ ，使用Anderson-Darling假设检验，验证得到：在5%显著性水平下，D也接受服从对数正态分布的原假设。结合试验数据，使用矩估计法，导出爆速的期望μ_D=7 710 m/s，标准差σ_D=321 m/s。代入式(8)，计算得到 $\hat{\mu}_{\rm{D}}$ =8.9503， $\hat{\sigma}_{\rm{D}}$ =0.0025。使用曲线拟合技术，得到PBX-9502爆速的概率密度函数曲线如图3所示。

图  3  PBX-9502爆速的概率密度函数

Figure  3.  PDF of detonation velocity of PBX-9502

下载: 全尺寸图片 幻灯片

2.2   不确定的唯象参数

EOS中的参数γ无法通过试验直接标定。由工程经验和专家意见，假设γ～B[2.93, 2.99, 6, 6]，其中B[a, b, m, n]代表4参数Beta分布，[a, b]表示参数的取值范围，m、n表示形状参数。根据经验，形状参数取m=6，n=6。与正态分布相比，Beta分布的优点在于取值范围有界。与均匀分布相比，Beta分布的PDF曲线光滑而不间断。γ的PDF的详细信息如图4所示。

图  4  多方指数γ的概率密度函数

Figure  4.  PDF of polytrophic exponent γ

下载: 全尺寸图片 幻灯片

3.   不确定度传播方法

3.1   高维非嵌入多项式混沌

设 $\boldsymbol{ \xi }= {({\xi _1},{\xi _2}, \cdots ,{\xi _d})^{\text{T}}}$ 是概率空间 $(\varOmega, {F}, {P})$ 中的 $d$ 维独立标准正态随机向量，其中Ω为样本空间， ${F}$ 为定义在Ω上的σ-代数， ${P}$ 为概率测度。L₂(Ω)为Ω上的平方可积函数空间，设 $u \in L_{2}(\varOmega)$ 代表爆轰系统的QoI。若赋予内积：

式中： $\omega \in \varOmega$ 代表基本事件， $\rho (x) = {(2{\text{π}} )^{ -{d}/{2}}}{{\rm{e}}^{ - {x^{\text{T}}}x/2}}$ 为随机向量 $\boldsymbol{\xi}$ 的联合概率密度函数，则 $L_{2}(\varOmega)$ 构成Gauss-Hilbert空间。利用Cameron-Martin定理^[17]，u可以通过 $\boldsymbol \xi$ 的正交多项式展开：

式中：指标(indice) $\alpha = {({\alpha _1},{\alpha _2}, \cdots ,{\alpha _d})^{\text{T}}}$ ，指标集 ${\mathcal{J}_q} = \{ \alpha \big||\alpha | {\text{≤}} q\}$ ， $|\alpha | = \displaystyle\sum\limits_{k = 1}^d {{\alpha _k}}$ ， ${u_\alpha }$ 为待定系数。

随机基函数：

式中： $\alpha ! =\displaystyle \prod\limits_{i = 1}^d {{\alpha _i}} !$ ； ${H_\alpha }(\xi )$ 为 $d$ 元Hermite多项式，定义为 ${H_\alpha }(\xi ) = {H_{{\alpha _1}}}({\xi _1}) \times {H_{{\alpha _2}}}({\xi _2}) \times \cdots \times {H_{{\alpha _d}}}({\xi _d})$ ； $H_{a_{k}}\left(\xi_{k}\right)$ 代表 $\alpha_{k}$ 次一元Hermite多项式。令 $\tilde{H}_{0}(x)=1, \tilde{H}_{1}(x)=x, \cdots, \tilde{H}_{n}(x)=x^{n}, \cdots$ 为Gauss-Hilbert空间中的一列线性无关函数，由Gauss-Hilbert空间中的内积定义式(9)和Gram-Schmidt正交化方法，计算得正交函数H_k(x), k=0，1，…，n，具体步骤如下：

同时，式(11)可表示为 ${\psi _\alpha } = \displaystyle\prod\limits_{i = 1}^d {{\psi _{{\alpha _i}}}} ({\xi _i}) = \prod\limits_{i = 1}^d {\frac{{{H_{{\alpha _i}}}({\xi _i})}}{{\sqrt {({\alpha _i})!} }}}$ 。且由Gram-Schmidt正交化方法知， $< {\psi _{{\alpha _i}}}{\psi _{{\alpha _j}}} > = \displaystyle\int_{{{{R}}^{{d}}}} {{\psi _{{\alpha _i}}}} (x){\psi _{{\alpha _j}}}(x)\rho (x){\text{d}}x = {\delta _{{\alpha _i}{\alpha _j}}}$ ， $\delta$ 代表Delta算子。

由 $\left\{\psi_{a}\right\}$ 的正交性和式(10)，得：

当d=1时，由式(10)、(13)和Hermite-Gauss求积公式，得到：

式中：ξ_r、w_r分别为求积点和权重。本文中取S=6，此时，Hermite-Gauss积分的求积点和权重如表1所示。

表  1  6节点Hermite-Gauss积分的求积点和权重^[21]

Table  1.  Quadrature points and weights for Hermite-Gauss integration with six points^[21]

$\xi_{r}$	$w_{i}$		$\xi_{r}$	$w_{i}$		$\xi_{r}$	$w_{i}$
–3.324 26	0.002 56		–0.616 71	0.408 83		1.889 18	0.088 62
–1.889 18	0.088 62		0.616 71	0.408 83		3.324 26	0.002 56

下载: 导出CSV 

| 显示表格

当d>1时，使用全张量积Hermite-Gauss求积公式，即：

实际应用中，式(10)通常截断为：

式中： $K = \dfrac{{(d + Q)!}}{{d!Q!}}$ ， $Q$ 为阶数，且在概率意义下， $\mathop {\lim }\limits_{K \to \infty } {u_K}(\xi ) = u(\xi )$ 。

3.2   Rosenblatt变换

由第2节知，输入不确定度不服从标准正态分布且概率类型不一致，因此无法直接使用3.1节所述PC理论。本文中使用Rosenblatt变换将一列相关随机变量转化成一列服从标准正态分布的独立随机变量。具体步骤为：

设 ${\boldsymbol{X}} = {({X_1},{X_2}, \cdots ,{X_d})^{\text{T}}}$ 为 $d$ 维非高斯随机向量，令 ${F_X}({x_1},{x_2}, \cdots ,{x_d})$ 为 $\boldsymbol X$ 的联合累积分布函数(cumulative density function, CDF)。构造如下映射关系Y=G(X)^[22]：

式中： $F_{x_{d}}\left(x_{d} \mid x_{1}, x_{2}, \cdots, x_{d-1}\right)$ 为条件CDF， $\varPhi$ 为标准正态变量的CDF，则 ${{Y}}$ 服从独立标准正态分布。

4.   结果的UQ分析

利用3.2节Rosenblatt变换技术将第2节标定的输入不确定度D、ρ₀和γ转化为相互独立的标准正态分布。然后利用3.1节非嵌入PC结合式(7)计算波后QoIs的概率密度函数。具体流程如图5所示。

图  5  不确定分析算法流程图

Figure  5.  Flow chart of uncertainty analysis

下载: 全尺寸图片 幻灯片

本文中，d=3，Q=3，因此K=20。接着利用表1，构造二维全张量Hermite-Gauss积分求积点，如图6所示。本文中计算将会用到的前6阶Hermite多项式的具体形式如图7所示。结合式(1)～(7)和(15)～(16)，利用图5所示的不确定度分析流程计算出QoIs的概率密度函数，更多信息如图8所示。可以看出，波后QoIs的概率密度函数没有明显的对称性，且峰度指标均大于3，处于尖峰状态，不服从Gauss分布。特别是波后密度，峰度指标最大，曲线最尖峭，且具有明显的偏差，与正态分布差别很大。波后压力的峰度系数最小，对称性最好，形状最接近正态分布，可用正态分布近似替代。

图  6  二维全张量积Gauss-Hermite求积点的位置

Figure  6.  Locations of Gauss-Hermite quadrature points used for two-dimensional tensor product

下载: 全尺寸图片 幻灯片

图  7  前6阶单变量Hermite多项式

Figure  7.  The first six orders of univariate Hermite polynomials

下载: 全尺寸图片 幻灯片

图  8  系统响应量的概率密度函数

Figure  8.  PDF of system response quantities

下载: 全尺寸图片 幻灯片

此外，QoIs的统计量，如期望和方差，可通过 $u_{\alpha}$ 显性表达，即：

利用式(18)计算出波后QoIs的期望，标准差，具体内容见表2。从表2可以看出，波后介质密度的标准差与期望的比值小于5%，而波后介质压力的标准差与期望的比值大于5%。波后介质压力波动较大，取值较宽，与文献[6, 23]的“爆轰压力测量值分散性较大”的结论相吻合。

表  2  波后感兴趣量的试验和统计结果

Table  2.  Statistical and experimental result of QoIs at the rear of the shock wave

波后感兴趣量	ρJ/(g·cm−3)	UJ/(m·s−1)	pJ/GPa
平均值	2.529	1 943	28.410
置信区间下限	2.524	1 789	24.010
置信区间上限	2.534	2 108	33.390
试验数据[24]	2.525	1 922	28.301
标准差	0.003	81	2.376

下载: 导出CSV 

| 显示表格

根据图8信息，可计算出95%置信水平下波后QoIs的置信区间。同时，波后介质的压力、密度和速度等物理量也能通过试验标定^{[6, 23]}。数值预测结果与洛斯阿拉莫斯国家实验室（Los Alamos National Laboratory, LANL）的试验结果^[24]作比对。比较发现，波后介质的压力、密度和速度的试验数据落入95%置信水平下的置信区间内，进一步确认了方法的有效性。

5.   结果与讨论

C-J理论建立了波后介质压力、速度、声速和密度等系统响应量和初始密度、爆速以及多方指数的显式函数关系。初始密度易于测量，爆速也容易通过较简单的试验标定。C-J理论成为试验之外的波后物理量标定手段，成本低，实用性强。

本文中，在概率框架下研究C-J理论，结果以PDF表征，允许QoIs在一定范围内变动。首先给出了输入不确定度的定量表征，使用Anderson-Darling假设检验法验证概率假设的正确性。运用收敛速度更快的非嵌入PC结合Rosenblatt变换研究了输入不确定度对波后QoIs的影响，发现试验结果落入QoIs的置信区间内。这说明密度、爆速等物理量在生产、加工和测试过程中的正常扰动不影响C-J理论的使用。同时，模型也允许多方指数在一定范围内变动。基于C-J理论的爆轰UQ结果能为装药设计、装甲防护等提供决策依据。方法具有普遍性，能推广到其他状态方程描述的爆轰UQ研究。

下一步拟减少假设、简化和近似的使用次数，将C-J理论UQ研究推广到更加复杂的情况，如非γ律EOS，甚至推广到ZND理论的UQ研究。对于非理想爆轰，声速面后可能存在放能。引入放能后，预测结果的差异，即不同模型导致的不确定度，也叫模型形式不确定度量化。这是国际热点，也是不确定度量化领域的难点之一。

另外，本文中，研究沿用Oberkampf的思路^[25]，即将试验结果看作是基准的，用以评估模拟结果的可信度。事实上，QoIs的试验结果也受到不确定因素的扰动，下一步拟综合试验不确定度和数值模拟不确定度开展研究。同时，从爆轰研究角度看，更希望从包含不确定度的某一参数实测结果中通过统计分析得到尽量多相关参数的均值和标准差等信息。若能指导实验设计方法，用尽量少的实验次数得到QoIs的结果，是另一个值得探讨的课题。