基于多项式混沌方法对C-J 爆轰参数不确定度的分析*

梁 霄，王瑞利，胡星志，陈江涛

（1. 山东科技大学数学与系统科学学院，山东 青岛 266590；2. 北京应用物理与计算数学研究所，北京 100094；3. 中国空气动力研究与发展中心，四川 绵阳 621000）

爆轰是极其剧烈的物理化学反应，炸药爆炸过程扩展的速度高达每秒数千米到万米之间，所形成的温度约为3 000～5 000 ℃，压力高达102～104MPa，发生在极短的时间（10-9s 量级）和极窄的空间（10-4m量级），释放极高的能量（102GW/cm2量级）[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×103m/s，速度的瞬时极大变化率导致缺乏精密的试验装置进行测试。传统的试验装置如圆筒试验，所用金属的纯度和密度不均匀性、加工公差的存在以及试验数据的拟合误差，都会导致测量结果的不确定性。同时，含能材料组分的细微变化、几何边界的不精确、极端荷载下金属延展性变化以及黏结剂的存在也会导致试验结果的波动[7-9]。

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

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

事实上，工程技术人员和学者们都已经意识到爆轰中的不确定因素[1,5-6,9,15]。文献[1, 6]甚至将物理量的容许不确定度（误差）定为 ±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 中，ρ0、p0、U0和e0分别代表未反应物的初始密度、压力、速度和内能，D为爆速。爆炸物以D–U0的速度进入波阵面，以D–U1的速度流出波阵面。

图1 爆轰波在高能炸药中的传播示意图Fig. 1 Sketch of propagation of detonation wave into high explosive

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

(1) 质量守恒：

(2) 动量守恒：

(3) 能量方程：

(4) Chapman-Jouguet 假设：

(5) 声速定义：

式中：v1=1/ρ1，为爆轰产物的比容；v0=1/ρ0，为未反应炸药的比容； ρ1为爆轰产物的密度；p1为波后压力；U1为波后介质速度；e1为爆轰产物的内能；c为声速；S代表等熵状态。

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

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

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

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

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

2.1 不确定物理量的表征

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

图2 PBX-9502 初始密度的概率密度函数Fig. 2 PDF of initial density of PBX-9502

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

图3 PBX-9502 爆速的概率密度函数Fig. 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 多方指数γ 的概率密度函数Fig. 4 PDF of polytrophic exponent γ

3 不确定度传播方法

3.1 高维非嵌入多项式混沌

设 ξ=(ξ1,ξ2,···,ξd)T是概率空间 ( Ω,F,P) 中的d维独立标准正态随机向量，其中Ω为样本空间，F为定义在Ω上的σ-代数，P为概率测度。L2(Ω)为Ω上的平方可积函数空间，设u∈L2(Ω) 代表爆轰系统的QoI。若赋予内积：

式中： ω ∈Ω 代表基本事件， ρ(x)=(2π)-d/2e-xTx/2为随机向量 ξ 的联合概率密度函数，则L2(Ω) 构成Gauss-Hilbert 空间。利用Cameron-Martin 定理[17]，u可以通过 ξ 的正交多项式展开：

随机基函数：

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

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

表1 6 节点Hermite-Gauss 积分的求积点和权重[21]Table 1 Quadrature points and weights for Hermite-Gauss integration with six points[21]

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

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

3.2 Rosenblatt 变换

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

设X=(X1,X2,···,Xd)T为d维非高斯随机向量，令FX(x1,x2,···,xd)为X的联合累积分布函数(cumulative density function, CDF)。构造如下映射关系Y=G(X)[22]：

式中：Fxd(xd|x1,x2,···,xd-1) 为条件CDF， Φ 为标准正态变量的CDF，则Y服从独立标准正态分布。

4 结果的UQ 分析

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

图5 不确定分析算法流程图Fig. 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 求积点的位置Fig. 6 Locations of Gauss-Hermite quadrature points used for two-dimensional tensor product

图7 前6 阶单变量Hermite 多项式Fig. 7 The first six orders of univariate Hermite polynomials

图8 系统响应量的概率密度函数Fig. 8 PDF of system response quantities

此外，QoIs 的统计量，如期望和方差，可通过uα显性表达，即：

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

表2 波后感兴趣量的试验和统计结果Table 2 Statistical and experimental result of QoIs at the rear of the shock wave

根据图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 的结果，是另一个值得探讨的课题。