304L/Q235B大面积金属板爆炸焊接物质点法模拟分析

王宇新; 李晓杰; 杨国俊; 范述宁; 王小红; 闫鸿浩

doi:10.11883/bzycj-2021-0198

摘要: 大面积金属板材304L/Q235B的爆炸焊接过程涉及炸药爆轰、金属板材的高速碰撞和塑性变形等。采用有限元法计算模拟这个问题时，网格单元会发生扭曲畸变现象，导致计算精度下降，甚至出现单元负体积而使计算终止，并且炸药爆轰形成气体产物飞散过程也很难模拟。为了能模拟大面积金属板材的爆炸焊接整个过程并获得合理的技术工艺参数，采用物质点法进行三维数值模拟分析。物质点法作为一种无网格法，在模拟冲击动力学问题中主要采用显式积分算法。通过将拉格朗日质点单元与固定的欧拉背景网格相结合，可以实现爆炸焊接的复板与基板的高速碰撞、炸药滑移爆轰、金属板面的塑性变形过程的数值模拟，并给出爆炸复合板材的形变、有效塑性应变和复板与基板的碰撞速度的计算结果。采用物质点法模拟的复合板材变形与爆炸焊接实验结果基本一致。计算复板与基板的碰撞速度这个重要的物理参数时，物质点法与Richter理论公式的相对误差不超过13%。数值计算和实验结果表明，物质点法在数值精度和计算效率方面具有优势，物质点法是研究金属焊接爆炸问题的一种有效数值方法。

关键词:

Abstract: The explosive welding process of large-area metal plates 304L/Q235B involves explosive detonation, high-speed collision and plastic deformation of two metal plates, etc. If the finite element method is used to simulate this problem, twist or distortion of the elements will lead to a decrease in the calculation accuracy, even negative volume elements will appear, resulting in the termination of the calculation. Moreover, it is difficult to use the finite element method to calculate and simulate the dispersion process of the gas products formed by explosive detonation. In order to simulate the whole process of explosive welding of large-area metal plates and obtain reasonable technical parameters, the material point method is used for three-dimensional numerical simulation analysis. As one of the meshless methods, the material point method mainly uses the explicit integral algorithm in the simulation of impact dynamics problems. By combining the Lagrangian particle elements with the fixed Euler background grid, this method can realize numerical simulations of the high-speed collision of clad and base plates, the explosive sliding detonation and the plastic deformation process of the metal plate surface, and obtain the results of deformation, effective plastic strain and collision velocity between the clad plate and the base plate. The deformation of the composite plate simulated by material point method is basically consistent with the experimental results of explosive welding. The collision velocity between the clad plate and the base plate is an important physical parameter. The relative error between the results of the material point method and Richter formula is less than 13%. Based on the numerical calculation by the material point method and the experimental results, it is shown that the material point method has advantages in numerical accuracy and computational efficiency, and it is also verified that the material point method is an effective numerical method to study metal explosive welding.

Key words:

material point method /
explosive welding /
large-area metal plates /
collision velocity between clad plate and base plate

意外爆炸、恐怖袭击、战争打击等爆炸事件会对居民的生命安全和公共财产构成严重威胁，研究爆炸对复杂结构的毁伤从而采取合理的措施是近年来的一个热点问题^[1-4]。爆炸发生时，爆源的能量以高温、高压气体的形式快速释放，其与周围空气相互作用产生冲击波。冲击波会对其经过的复杂结构产生巨大的爆炸载荷，爆炸载荷可通过关键位置处的超压峰值、正压冲量等冲击波参数进行量化。对爆炸载荷进行快速估计，一方面是复杂结构抗爆设计时多场景和多变量优化的关键，另一方面也是爆炸发生后进行紧急救援的重要信息支撑。近年来，随着人工智能（artificial intelligence, AI）的发展，深度学习方法在爆炸载荷估计领域的研究逐渐兴起。利用爆炸试验数据或模拟数据训练得到的深度学习模型能够对超压峰值、正压冲量等表征爆炸载荷的冲击波参数进行快速、准确的估计，目前已经有学者进行了相关研究。

基于深度学习方法对爆炸载荷进行估计，国内外学者的研究可以分成2个阶段。在早期阶段，研究人员专注于探索深度学习在爆炸载荷估计领域的适用性，用简单的神经网络实现爆炸载荷的快速估计。这一阶段的研究问题相对简单，神经网络的结构和输入输出模式并不复杂，浅层的全连接神经网络（fully connected neural network, FCNN）在爆炸载荷估计中表现出了较好的性能。对于防爆墙场景下目标结构的爆炸载荷估计问题，Remennikov等^[5]在不同距离和不同高度处均匀设置了压力测点，结合实验数据，优化了所构建的FCNN的隐藏层配置，实现了对超压峰值和峰值比冲的较为准确的估计。针对城市场景下的爆炸载荷估计，Remennikov等^[6]设计了街道宽度、街道距离和建筑高度可变的不同城市场景，结合爆炸模拟器Air3D生成了组合不同街道特征的数值模拟数据，在此基础上训练所得的FCNN模型对于关键位置处超压峰值和脉冲估计的决定系数（R²）分别为0.993和0.995，证明了FCNN能够对城市场景爆炸载荷实现准确估计。Dennis等^[7]将FCNN应用在近场爆炸情况下的爆炸载荷估计，同样取得了较好的效果。在一般FCNN模型的基础上，Flood等^[8-9]和Bewick等^[10]构建了径向高斯增量网络（radial Gaussian incremental network, RGIN），用高斯单元作为隐藏层神经元，并在训练过程中通过增量的方式构建网络，基于经验建模系统PURWall^[8-9]或数值模拟软件DYSMAS^[10]生成的爆炸数据优化网络隐藏层配置，实现了对防爆墙后结构的超压峰值的较为准确的估计。

在第2阶段，研究人员专注于将冲击波传播时的物理特性以及空间特征与更先进的深度学习模型相结合，构建性能更强的基于深度学习的爆炸载荷估计方法。在这一阶段，物理信息与场景空间信息的融入使得神经网络的结构和输入输出模式更加复杂，深度学习模型对于不同场景下爆炸载荷估计的性能也进一步提升。

考虑到用于爆炸分析的神经网络通常需要提供特定研究场景的输入，这在一定程度上限制了深度学习方法在特定场景之外的爆炸载荷估计性能，Dennis等^[11]将冲击波传播过程中的物理信息与神经网络相结合，在输入特征中考虑了波反射和叠加效应，构造了一种新的方向编码神经网络（direction-encoded neural network, DENN）。DENN通过底层物理过程提供的信息学习冲击波的行为，以感兴趣位置点（point of interest, POI）与周围障碍物的距离以及冲击波到达POI的最短路径作为输入，以POI的超压峰值作为输出，实现对二维障碍物环境下爆炸载荷的快速、准确的估计。

为进一步提高深度学习方法在近场爆炸载荷估计中的准确性和泛化性，Pannell等^[12]将冲击波传播的物理原理与神经网络相结合，构造了一种物理引导神经网络（physics-guided neural network, PGNN）。近场冲击波传播过程中的物理原理被作为正则化项整合到神经网络的损失函数中，有效提高了神经网络的爆炸载荷估计性能。

考虑到不同爆源产生的冲击波有不同的特性，Pannell等^[13]构造了一种迁移神经网络（transfer neural network, TNN），基于已有先验爆源知识的神经网络估计不同爆源产生的爆炸载荷。将使用球形爆源的近场峰值比冲分布数据集训练得到的FCNN模型与一个附加网络（用于处理新的长径比参数）相结合作为TNN模型，用该TNN模型估计圆柱形爆源在相似比例距离范围内产生的峰值比冲，结果表明，相比于普通的神经网络，TNN在估计性能上有所提升。该方法通过将爆源特性与神经网络相结合，在提升深度学习模型性能的同时，提高了爆炸数据的利用率，减少了计算资源消耗。

为了将场景空间信息与神经网络相结合，Kang等^[14]构建了一种三维卷积神经网络（3D convolutional neural network, 3DCNN）模型，利用三维卷积核有效地提取三维空间信息，从而更好地捕捉建筑布局和测点超压峰值之间的关系。结果表明，对于三维建筑环境中的爆炸载荷预测，3DCNN模型的预测精度远高于FCNN模型，且计算时间更短。

此外，Huang等^[15]将超分辨率重构（super-resolution reconstruction, SRR）引入爆炸领域。SRR是一种基于深度学习的数据重构方法，通过学习大量低分辨率到高分辨率的数据对之间的映射关系构建神经网络。结果表明，SRR能够以较低的计算成本实现二维爆炸场的超分辨率重构，因此，能够降低获取高分辨率爆炸场数据的计算成本。Li等^[16-19]将集成学习算法^[16]、FCNN^[17]、图神经网络^[18]、基于自注意力机制的神经网络^[19]等方法运用在沸腾液体膨胀蒸气爆炸（boiling liquid expanding vapor explosions, BLEVEs）的爆炸载荷估计，结合商业软件FLACS获取的BLEVEs数值模拟数据，对所构建的模型进行训练和优化，取得了较好的爆炸载荷估计效果，证明了先进的深度学习模型在爆炸载荷估计领域的应用潜能。此外，深度注意网络^[20]、长短期记忆神经网络（long short-term memory, LSTM）^[21]、基于卷积自编码器和LSTM神经网络的降阶模型^[22]以及随机森林^[23]等方法已被应用于空气动力学领域的气动预测研究，均实现了较高的预测精度。基于深度学习的方法开发的气动和流动模型比理论模型的精度更高，比数值模拟方法的计算成本更低^[24]。

上述研究证明，普通的神经网络模型能够对简单场景下的爆炸载荷实现较为准确的估计，且比传统数值模拟方法的计算耗时更少。而与冲击波传播的物理知识及空间特征等相结合的更先进的深度学习模型，对不同场景下爆炸载荷的估计性能更强。先进的深度学习模型在爆炸载荷估计领域能够取得较好的效果。但对于复杂结构的爆炸载荷估计问题而言，现有研究还存在一些不足。目前，对于实际场景中爆炸载荷估计的研究多集中在城市建筑群环境以及简单防爆墙场景等，在场景设计阶段，目标结构常被简化为无形变的规则长方体，但现实中的复杂结构（如飞机、汽车、楼房等）往往是不规则的几何体，现有研究中建模的结构体不够复杂。另外，现有研究所考虑的爆源情况（如爆源位置、当量）较少，并很少考虑爆源的起爆速度，因此，不能满足现实场景中爆源情况复杂多变的要求。此外，现有的基于深度学习的爆炸载荷估计方法虽然能够取得较好的效果，但只能给出关键位置处爆炸载荷的点估计，没有量化估计结果的不确定性。不确定性量化的缺失意味着模型无法提供估计结果的置信度，在实际决策中不能帮助用户决定是否能信任该结果，而估计结果的可信性对爆炸载荷估计领域十分重要。

针对现有研究的不足，本文中，提出一种基于贝叶斯深度学习（Bayesian deep learning, BDL）的复杂结构爆炸载荷估计方法。在场景构造阶段，基于建模软件SpaceClaim构造包含圆柱体、长方体及其拼接结构的复杂结构，并在结构表面的关键位置处设置200个压力测点。构造的场景考虑30种不同的爆源位置、30种不同的当量以及10种不同的起爆速度，以模拟现实爆炸事件中爆源情况复杂多变的情况。基于开源软件blastFoam进行数值模拟，获得120组针对目标结构的不同爆炸情况下关键位置处的超压峰值数据集。在此基础上，构建BDL模型，以爆源位置、当量以及测点的位置和比例距离等作为输入，以超压峰值作为输出，对复杂结构的爆炸载荷进行快速估计，并给出该结果不确定性的量化表述，以检验估计结果的可信性。

1. 场景构建与数据生成

1.1 场景构建

爆炸发生时复杂结构体爆炸载荷的快速估计方法需要基于建模软件构造复杂结构体，并在此基础上设计符合现实爆炸事件的爆炸场景。发生在现实场景中的爆炸事件具有不可预知性，爆源的当量、起爆位置以及起爆速度复杂多变。因此，为了确保所构建的爆炸载荷快速估计方法的泛化性和普适性，在场景构建阶段需要考虑在较大范围内可变的爆源参数。

基于建模软件SpaceClaim构建复杂结构模型，如图1所示，其中目标结构中心为笛卡尔坐标系的坐标原点（星标所示位置）。结构包含放置方式不同的圆柱体、长方体及其拼接体，用以模拟现实结构体的多样形态。为增强所构建方法对不同爆源情况的泛化性，构建的场景（表1）考虑了30种不同的爆源位置（距离目标结构中心40～100 m范围内随机选取）、30种不同的当量（200～2 000 kg TNT）和10种不同的起爆速度（0～1 020 m/s），速度的方向均设置为指向坐标原点。实验中，将上述3种爆源相关参数的取值进行随机组合，共得到9 000种不同参数设置的算例。为减少计算资源的消耗，从中随机抽取120组算例进行数值模拟，生成训练深度学习模型所需的数据。

图 1 复杂结构模型及压力监测点分布

Figure 1. Complex structure model and distribution of pressure monitoring points

下载: 全尺寸图片幻灯片

表 1 爆源参数的范围

Table 1. The range of explosion source parameters

爆源参数	符号	单位	参数范围
当量	Q	kg	(200, 2000)
爆心距	R	m	(40, 100)
起爆速度	U	m/s	(0, 1020)

下载: 导出CSV

| 显示表格

为全面评估目标结构的爆炸载荷，每次数值模拟时，都在目标结构的200个关键位置处（图1中的红色圆点）设置压力测点，记录冲击波传播过程中的压力。模拟结束后，提取每个压力测点的超压峰值，以超压峰值这一重要的冲击波参数作为结构关键位置处爆炸载荷的量化指标。

综上所述，经过数值模拟后，共获得120×200条样本数据，即最终的爆炸数据集包含24 000条样本，用于后续所构建的深度学习模型的训练和测试。

1.2 数据生成

1.2.1 基于blastFoam的数值模拟

blastFoam是一款用于模拟冲击波传播的开源计算流体动力学（computational fluid dynamics, CFD）软件。基于blastFoam进行数值模拟，可以对爆炸事件的影响进行定量分析，从而有助于在爆炸事件发生前制定有效的安全防护措施和应急方案，以及在爆炸事件发生后进行迅速响应和灾后紧急救援。

基于数值模拟软件blastFoam，对背景空气域和球形爆源进行建模，并对冲击波的传播进行计算。空气采用理想气体状态方程，TNT采用Jones-Wilkins-Lee（JWL）状态方程^[25]：

$\rho=A\left(1-\frac{\lambda_1\eta}{R_1}\right)\mathrm{e}^{-\frac{R_1}{\eta}}+B\left(1-\frac{\lambda_1\eta}{R_2}\right)\mathrm{e}^{-\frac{R_2}{\eta}}+\lambda_1\rho e$

(1)

式中： $\rho$ 为密度； $e$ 为比内能； $\eta =\rho /{\rho }_{0}$ 为相对密度， ${\rho }_{0}$ 为初始密度； $A$ 、 $B$ 、 ${R}_{1}$ 、 ${R}_{2}$ 、 $\lambda_1$ 为经验常数。经验常数的取值如表2所示。

表 2 经验常数的取值

Table 2. The values of empirical constants

A/TPa	B/TPa	R₁	R₂	λ₁	e/(TJ·m⁻³)	ρ₀/(kg·m⁻³)
609.77	12.95	4.50	0.9	0.25	9.0	1 630

下载: 导出CSV

| 显示表格

假定复杂结构体的所有壁面和地面都是无形变的，即采用刚性边界条件。计算域的四周和顶部（距离感兴趣的区域足够远处）采用流出边界条件。计算域用空气填充，其初始参数与正常大气条件相同，即密度 $\rho$ =1.225 kg/m³，温度T=298.15 K，静压p_s=101.325 kPa，比定容热容c_p=1004 J/(kg·K)。

1.2.2 网格敏感性分析

数值模拟中，基于blastFoam提供的blockMesh和snappyHexMesh等网格绘制工具划分背景网格以及细化结构附近的网格，并使用dynamicMesh工具进行自适应网格细化（adaptive mesh refinement, AMR），对冲击波传播位置处的网格数量进行动态调整。为验证数值模拟数据的可靠性，进行网格敏感性分析，在不同的网格尺寸下，将数值模拟结果与经验公式计算结果进行比较，以获取兼顾计算速度和计算精度的网格尺寸。

设置当量分别为200和2000 kg TNT，在距离爆源10～150 m范围内均匀布置一系列测点（每隔10 m放置一个测点），记录不同背景网格设置下冲击波传播过程中测点的超压峰值。基于经验公式^[26]，计算测点的超压峰值：

$\mathrm{\Delta}p_0=\left\{1+0.53Ma\frac{1}{0.27+0.48d}\left[0.32+0.80\;\mathrm{sin}\left(\frac{\theta+2.42}{3.58}{\text{π}} \right)\right]\right\}\times\left(\frac{0.067}{d}+\frac{0.301}{d^2}+\frac{0.431}{d^3}\right)$

(2)

式中： $\mathrm{\Delta }{p}_{0}$ 为超压峰值，Pa；Ma为爆源的马赫数，即起爆速度U与音速（340 m/s）的比值；d为比例距离， $\mathrm{m}/\mathrm{kg}^{\frac{1}{3}}$ ； $\theta$ 为方位角，即爆心连线与爆源运动方向的夹角，rad。

数值模拟计算的测点的超压峰值受网格尺寸影响，将式(2)与数值模拟结果进行对比，如图2所示，其中红色实线代表位于不同位置处各测点的超压峰值经验值，2条黑色细实线代表超压峰值经验值的上下10%包络线，当数值模拟结果位于2条黑色细实线之间时，表明模拟效果较为准确。从图2可以看出：网格尺寸为10.0 m时，超压峰值的数值模拟结果（黄色实线）位于2条黑色细实线之外；网格尺寸为5.0 m时，超压峰值的数值模拟结果（蓝色实线）位于2条黑色细实线之间，且与网格尺寸为7.5 m时的数值模拟结果（绿色实线）相比，更接近超压峰值的经验值（红色实线）。因此，网格尺寸为5.0 m时，超压峰值的数值模拟结果最为准确。网格尺寸不同时数值模拟所需的时间如表3所示。可以看出，随着网格尺寸的减小，数值模拟与经验公式的计算结果越来越接近，模拟所需时间也越长。当背景网格尺寸为5.0 m，结构附近细化等级为5，AMR（adaptive mesh refinement）动态细化等级为1时，数值模拟与经验公式计算结果的误差不超过10%，因此，数值模拟时，背景网格尺寸设置为5.0 m，结构附近细化等级为5，AMR动态细化等级为1。

图 2 网格尺寸不同时数值模拟与经验公式计算的结果对比

Figure 2. Comparison of numerical simulation results with empirical formulas under different grid settings

下载: 全尺寸图片幻灯片

表 3 不同背景网格设置下数值模拟所需的计算时间

Table 3. Calculation time required for numerical simulation with different background grid settings

背景网格尺寸/m	5.0	7.5	10.0
计算时间/h	10	2	1

下载: 导出CSV

| 显示表格

2. 基于贝叶斯深度学习的爆炸载荷快速估计方法

与传统深度学习不同，贝叶斯深度学习将网络参数视为服从一定概率分布的随机变量，并采用贝叶斯推理来训练预测模型^[27-28]。基于贝叶斯深度学习构建的爆炸载荷快速估计模型能在不牺牲预测速度的前提下获得更好的预测精度，并使模型具备不确定性量化能力^[29-30]。以下先介绍贝叶斯深度学习原理和模型实现细节，最后评估所提出的模型的性能。

2.1 基于变分推断的贝叶斯深度学习

基于变分推断的贝叶斯神经网络是一种结合了贝叶斯推理和神经网络技术的模型。该模型通过在神经网络的权重上引入概率分布，并利用变分推断来近似这些分布的后验，提高了对不确定性的处理能力。定义贝叶斯神经网络 $y={f}_{{\varTheta }}\left({x}\right)$ ，其中 ${\varTheta }=\left\{{{\theta }}_{1},{{\theta }}_{2},\dots \right\}$ 为参数集合，x和y分别表示神经网络的输入和输出。考虑一个具有M个样本的数据集 $D=\left\{X,Y\right\}$ ，其中 $X={\left\{{{x}}_{{i}}\right\}}_{i=1}^{M}$ 表示输入， $Y={\left\{{y}_{i}\right\}}_{i=1}^{M}$ 表示标签。给定参数先验分布 $p\left({\varTheta }\right)$ ，似然函数可以表示为 $P\left(Y|X,{\varTheta }\right)=\mathop{\prod}\limits_{i=1}^{M}p\left({y}_{i}\right|{x}_{i},{\varTheta })$ 。根据贝叶斯定理，模型参数 ${\varTheta }$ 的后验分布为：

$p\left({\varTheta }|X,Y\right)=\frac{p\left({\varTheta }\right) \mathop{\prod}\limits_{i=1}^{M}p\left({y}_{i}|{{x}}_{{i}},{\varTheta }\right)}{{\int }_{{\varTheta }}^{}p\left({\varTheta }\right)\mathop{\prod}\limits_{i=1}^{M}p\left({y}_{i}|{{x}}_{{i}},{\varTheta }\right){\mathrm{d}}{\varTheta }}$

(3)

对于新的输入 $x'$ ，基于后验分布，贝叶斯神经网络给出的输出 $y'$ 的预测分布为：

$p\left({y}'|{x}',X,Y\right)=\int p\left({\varTheta }|X,Y\right) p\left({y}'|x',{\varTheta }\right){\mathrm{d}}{\varTheta }$

(4)

贝叶斯深度学习的关键在于获得参数的后验分布。由式(3)可知，参数后验分布 $p\left({\varTheta }|X,Y\right)$ 的计算涉及高维积分，难以直接解析获得。变分推断是一种近似推断概率模型中未知变量后验分布的方法，通过在某个预定义的分布类中搜索最接近真实后验分布的分布来近似计算后验分布。其基本思想是将推断问题转化为一个优化问题，通过优化变分分布的参数来逼近真实后验分布。定义一个简单的变分分布 ${q}_{\lambda }\left({\varTheta }\right)$ ，如多元正态分布 $\mathcal{N}\sim\left(\mu,\mathit{\Sigma}\right),\lambda=\left\{\mu,\mathit{\Sigma}\right\}$ ，其中μ为均值， $\mathit{\Sigma}$ 为标准差。变分推断通过调优变分参数 $\lambda$ ，尽可能使 ${q}_{\lambda }\left({\varTheta }\right)$ 接近真实的后验分布 $p\left({\varTheta }|X,Y\right)$ ，采用KL（Kullback-Leibler）散度来量化2个分布的相似度：

$\begin{split}KL\left[{q}_{\lambda }\left({\varTheta }\right)\left|\right|p\left({\varTheta }|D\right)\right] =&\int {q}_{\lambda }\left({\varTheta }\right)\mathrm{lg}\left( \frac{{q}_{\lambda }\left({\varTheta }\right)} {p\left({\varTheta }|D\right)} \right){\mathrm{d}}\varTheta = \int {q}_{\lambda }\left( {\varTheta }\right) \mathrm{lg} \left(\frac{{q}_{\lambda }\left({\varTheta } \right)}{p\left({\varTheta },D\right)/p\left(D\right)}\right){\mathrm{d}}\varTheta =\\ &\int {q}_{\lambda }\left({\varTheta }\right)\mathrm{lg}\,p\left(D\right){\mathrm{d}}\varTheta - \int {q}_{\lambda }\left({ \varTheta }\right)\mathrm{lg}\left(\frac{p\left({\varTheta },D\right)}{{q}_{\lambda }\left({\varTheta }\right)} \right) d\varTheta =\mathrm{lg}\,p\left(D\right)-\int {q}_{\lambda }\left({\varTheta }\right)\mathrm{lg}\left(\frac{p\left({\varTheta },D\right)}{{q}_{\lambda }\left({\varTheta }\right)}\right){\mathrm{d}}\varTheta \end{split}$

(5)

$\begin{split}{\lambda }^{*}= & \mathrm{arg}\underset{\lambda }{\mathrm{min}}\left\{-\int {q}_{\lambda }\left({\varTheta }\right){{\mathrm{lg}}}\left(\frac{p\left({\varTheta },D\right)}{{q}_{\lambda }\left({\varTheta }\right)}\right){\mathrm{d}}{\varTheta }\right\} =\mathrm{arg}\underset{\lambda }{\mathrm{min}}\left\{\int {q}_{\lambda }\left({\varTheta }\right){{\mathrm{lg}}}\left(\frac{{q}_{\lambda }\left({\varTheta }\right)}{p\left({\varTheta }\right)}\right){\mathrm{d}}\varTheta -\int {q}_{\lambda }\left({\varTheta }\right){{\mathrm{lg}}}\,p\left(D|{\varTheta }\right){\mathrm{d}}{\varTheta } \right\}=\\ & \mathrm{arg}\underset{\lambda }{\mathrm{min}}\left\{KL\left[{q}_{\lambda }\left({\varTheta }\right) \left|\right|p\left({\varTheta }\right)\right]-{E}_{{\varTheta }\sim {q}_{\lambda }}\left[{{\mathrm{lg}}}\,p\left(D|{\varTheta }\right)\right]\right\}\end{split}$

(6)

式中： $KL\left[{q}_{\lambda }\left({\varTheta }\right)\left|\right|p\left({\varTheta }\right)\right]$ 为简单变分分布 ${q}_{\lambda }\left({\varTheta }\right)$ 与模型参数先验分布 $p\left({\varTheta }\right)$ 之间的KL散度，λ^*为最优变分参数， ${E}_{{\varTheta }\sim{q}_{\lambda }}\left[{\mathrm{lg}}\,p\left(D|{\varTheta }\right)\right]$ 表示对数似然函数。

对于每一个测试样本，基于贝叶斯神经网络模型的输出是一组概率分布而不是确定的点估计结果。对于给定的测试样本 ${X}'$ ，输出预测目标 ${Y}'$ 的分布为：

$p\left({Y}'|{X}',D\right)=\int p\left(Y|X,{\varTheta }\right){\mathrm{lg}}\,p\left({\varTheta }|D\right){\mathrm{d}}{\varTheta }$

(7)

目前，常用的不确定性量化指标有预测方差、置信区间等。基于统计理论，置信区间通过给定置信水平，估计真实参数值在区间内的可能性。本研究采用预测结果的95%置信区间作为不确定性量化指标，95%置信区间越大，模型对预测结果的不确定性越高。

基于变分推断的BDL模型与普通深度学习网络对训练数据的要求并无明显区别。但相较于仅给出点估计的普通深度学习网络，BDL模型能够输出预测结果的概率分布，因此，它能量化预测的不确定性。对于爆炸载荷预测这样需要高可靠性预测的领域，量化不确定性十分重要。并且，通过在神经网络的权重上引入先验分布并考虑后验分布，BDL模型能够在训练过程中自然地实现正则化，从而在一定程度上防止过拟合。此外，深度神经网络被称为不可解释的黑盒模型，而通过对权重的后验分布进行分析，BDL模型能提供关于模型参数的额外信息（如置信区间），有助于增强模型的可解释性。然而，BDL模型存在一些局限性，如由于训练过程中对后验分布的近似和采样，BDL模型会比传统的深度学习网络训练更为复杂和耗时。计算资源有限时，不宜采用BDL模型。

2.2 实现细节

2.2.1 数据集划分与数据预处理

在模型训练阶段，将120组样本划分为训练集和测试集，并采用5次5折交叉验证，即每次按照4∶1的比例划分训练集和测试集，进行5次重复实验，取每次独立训练和测试结果的均值为最终结果。

输入特征的质量影响深度学习模型的性能。原始数据集的输入特征包括爆源位置（爆心位于 ( ${x}_{{\mathrm{b}}},{y}_{{\mathrm{b}}}, {{\textit{z}}}_{{\mathrm{b}}})$ ）、爆源当量Q、起爆速度U以及测点的位置坐标( ${x}_{{\mathrm{p}}},{y}_{{\mathrm{p}}},{{\textit{z}}}_{{\mathrm{p}}})$ 。为了提高模型的性能，在原有特征的基础上进一步构造新的特征。首先，将速度矢量分解为x、y、z方向的3个分量 ${v}_{x}、{v}_{y}、{v}_{{\textit{z}}}$ ，作为新的输入特征。此外，由式(2)可知，超压峰值主要与测点到爆源的比例距离d有关，将1/d、1/d²、1/d³作为新的特征。同样地，将测点与起爆速度方向的夹角 $\theta$ 也作为输入特征之一。因此，深度学习模型的输入特征向量x可表示为：

${{\boldsymbol{x}}}=\left[x_{\mathrm{b}},y_{\mathrm{b}},\text{z}_{\mathrm{b}},Q,\mathit{U},x_{\mathrm{p}},y_{\mathrm{p}},\text{z}_{\mathrm{p}},v_x,v_y,v_{\text{z}},d,\frac{1}{d},\frac{1}{d^2},\frac{1}{d^3},\theta\right]$

(8)

实验结果表明，增加了特征丰富性的输入向量比仅包含原始特征的输入向量获得的模型训练效果更好。

2.2.2 基于网格搜索确定最优超参数

构建基于变分推断的贝叶斯深度神经网络，估计复杂结构关键位置处的超压峰值，采用网格搜索的方法确定模型结构的最优超参数，最后通过训练学习获得模型的参数。

基于网格搜索确定最优超参数，需要先定义待调节的超参数及其取值范围，并根据每个超参数的候选值生成超参数的组合网格；然后，采用网格搜索算法遍历所有的超参数组合，对每个参数组合进行模型训练和评估；最后，将在验证集上表现最优的超参数组合作为最优超参数。本研究中，定义的待调节的超参数及其取值范围如表4所示，其中ReLU和Leaky ReLU为2种常用的激活函数，前者称为修正线性单元，后者称为泄露修正线性单元。基于最优超参数，在训练集上重新对模型进行训练，在测试集上对模型性能进行评估。

表 4 待调节的超参数及其取值范围

Table 4. The hyperparameter to be adjusted and the corresponding range

神经网络层数	每层神经元数	激活函数	学习率
2～8	16, 32, 64, 128, 256	ReLU, Leaky ReLU	0.001, 0.005, 0.010

下载: 导出CSV

| 显示表格

2.3 模型评估

本研究采用5次5折交叉验证评估模型的性能，采用平均绝对百分比误差（mean absolute percentage error, MAPE）和决定系数（R²）作为评价模型精度的指标。其中，δ_mape（MAPE值）是相对误差，完美模型的δ_mape为0%，劣质模型的δ_mape超过100%。决定系数（R²）是评价回归模型的拟合优劣程度的统计指标，R²的取值范围为[0, 1]，R²越大，模型的拟合效果越好。δ_mape和R²可表示为：

${\delta }_{\mathrm{mape}}=\frac{1}{M}\sum_{i=1}^{M}\left|\frac{{y}_{l,i}-\widehat{{y}_{i}}}{{y}_{l,i}}\right|\times 100{\text{%}}$

(9)

${R}^{2}=1-\frac{{\displaystyle\sum }_{i=1}^{M}{\left({y}_{l,i}-{\widehat{y}}_{l}\right)}^{2}}{{\displaystyle\sum }_{i=1}^{M}{\left({y}_{l,i}-\overline{y}\right)}^{2}}$

(10)

式中： ${y}_{l,i}$ 和 $\widehat{{y}_{i}}$ 为测试样本的真实值和预测值， $\overline{y}$ 为所有测试样本中目标的平均值，M为样本数量。

采用预测区间覆盖概率（prediction interval coverage probability, PICP）和归一化平均预测区间宽度（normalized mean prediction interval width, NMPIW）来衡量贝叶斯神经网络对预测不确定性的量化能力。p_picp（PICP值）反映了真实值落在预测区间内的概率，用来度量模型置信区间的准确性。p_picp的取值范围为0～1，p_picp越趋近于1，模型预测的置信区间越准确。w_nmpiw（NMPIW值）进行了归一化处理，它衡量了平均预测区间的宽度。w_nmpiw越小，预测区间越紧凑，模型对预测不确定性的量化能力越好。p_picp和w_nmpiw可表示为：

$p_{\mathrm{picp}}=\frac{1}{M}\sum_{i=1}^Mc_i,\qquad\quad c_i=\left\{\begin{array}{ll}1 & \qquad i_{\mathrm{l}}\text{＜}y_{\mathrm{r},i}\text{＜}i_{\mathrm{u}} \\ 0 & \qquad\mathrm{else}\end{array}\right.$

(11)

$w_{\mathrm{nmpiw}}=\frac{\dfrac{1}{M}\displaystyle\sum_{i=1}^M i_{\mathrm{u}}-i_{\mathrm{l}}}{G}$

(12)

式中： $y_{\mathrm{r},i}$ 为第i个样本的真实值， $i_{\mathrm{l}}$ 和 $i_{\mathrm{u}}$ 分别为第i个样本的预测置信区间的下限和上限，G为所有测试样本中预测最大值与最小值的差值。

3. 实验结果

3.1 实验平台

数值模拟分别在天河二号超级计算中心和本地台式电脑上进行，模拟平台的相关参数如表5所示。天河二号计算中心上，采用blastFoam6.2.0软件进行计算，每个算例设置144个并行度，模拟时间设置为0.3 s，每个算例的平均运行时间约为2 h。在本地台式电脑上，基于python3.7，构建和训练贝叶斯深度学习模型。

表 5 模拟平台的参数

Table 5. Parameters of the simulation platform

平台	处理器	内存	操作系统
天河二号超级计算中心	2×12 Intel Xeon E5-2692 v2/单节点	128 GB/单节点	Linux lon26 3.10.0-514.el7.x86_64
本地台式机	AMD Ryzen7 3700x8-Core八核	32 GB	Windows10

下载: 导出CSV

| 显示表格

3.2 模型结构

采用TensorFlow probability（TFP）库和负对数似然（negative log-likelihood, NLL）构建贝叶斯神经网络和损失函数。基于网格搜索方法，最终所构建的贝叶斯深度学习模型包含4个隐藏层，每个隐藏层中包含128个神经元，激活函数为ReLU，学习率为0.005。

3.3 模型估计结果及其评估

从估计精度、不确定性量化、估计速度3方面评估模型效果。基于特征输入和所构建的贝叶斯深度学习模型，对于120组样本进行5次5折交叉验证。训练过程中，模型损失与训练迭代次数的关系如图3所示，其中黄色和蓝色曲线分别为训练集和验证集上的损失。可见，随着训练次数的增多，模型损失逐渐下降。当模型损失趋于平缓时，模型训练完毕。最终，模型对超压峰值的估计精度和不确定性如表6所示。

图 3 模型训练过程中损失与训练迭代次数的关系

Figure 3. Relationship between iteration times and loss during model training

下载: 全尺寸图片幻灯片

表 6 超压峰值的估计精度和不确定性

Table 6. The estimation accuracy and uncertainty of the peak overpressure

估计精度				不确定性
δ_mape		R²		p_picp		w_nmpiw
训练集	测试集	训练集	测试集	训练集	测试集	训练集	测试集
0.104	0.122	0.834	0.830	0.921	0.816	0.041	0.057

下载: 导出CSV

| 显示表格

由表6可知，对于复杂结构上的200个关键位置，所构建的贝叶斯深度学习模型在训练集和测试集上的平均MAPE约为10%和12%，说明模型不存在过拟合现象，能较准确地估计复杂结构的爆炸载荷。在训练集和测试集上，模型的PICP值都稳定在81%以上，NMPIW值保持在5.7%以下，说明模型的估计值具有较高的准确性。

为更好地分析模型的预测效果，将样本按照真实超压峰值的大小进行升序排序和编号，模型在训练集和测试集上的超压峰值的估计结果如图4所示，其中，蓝色区域表示预测结果的95%置信区间。可以看出，蓝色区域面积较小，说明预测结果的不确定性较小，模型对预测值的置信度较高，预测结果可靠。

图 4 超压峰值的估计结果

Figure 4. The estimated result of the peak overpressure

下载: 全尺寸图片幻灯片

对于单组测试样本（200条样本数据），模型对超压峰值的估计时间在4 s以内；对于单个样本，模型对超压峰值的估计时间约为20 ms。与数值模拟方法相比，所构建的BDL模型对超压峰值估算有明显的速度优势。

3.4 模型预测分析

为更好地理解模型的行为，对BDL模型进行更详细和严格的性能分析，将120组样本按照1∶9的比例划分为测试集和训练集，即测试集包含12×200条样本，训练集包含108×200条样本。根据目标值范围和比例距离，在测试集上分析模型预测的误差。

3.4.1 散点图

图5对比了数值模拟得到的超压峰值（真实值）和BDL模型预测的超压峰值（预测值），其中红色实线为预测结果与目标一致的完美模型，两侧的黑色实线为相对误差为±30%的包络线，圆点代表2 400个样本点，圆点的颜色反映测点的数量密度，密度越大，数量占比越接近1。在整个测试集上，BDL模型预测的超压峰值较为准确，δ_mape为12.2%，R²为0.848，接近1，大部分样本点都在2条黑色实线内，相对误差不超过30%。按照真实值的大小将2 400个测试样本分成2组。图6(a)为模型对超压峰值在0～100 kPa范围内的测点的预测效果，δ_mape为11.5%，样本点基本位于2条黑色实线之间。图6(b)为模型对超压峰值在100～400 kPa范围内的测点的预测效果，相对于图6(a)，图6(b)的R²更小（样本量小），样本点分布更加分散，δ_mape为19.9%，说明超压峰值较大时，BDL模型对样本的预测性能不太理想。从图5～6可以看出，BDL模型能较准确地预测样本的超压峰值，当超压峰值较小时，BDL模型的预测性能较好。

图 5 超压峰值的真实值与估计值

Figure 5. The true and estimated values of the peak overpressure

下载: 全尺寸图片幻灯片

图 6 范围不同时超压峰值的真实值与估计值

Figure 6. The true and estimated values of the peak overpressure for different ranges

下载: 全尺寸图片幻灯片

3.4.2 基于目标范围的误差分解

将2 400个样本根据超压峰值的取值区间分为4组，表7给出了预测结果在不同的MAPE区间内的样本数量占该组测试子集样本总数的百分比（φ）。可以看出：超压峰值较小时，BDL模型的预测精度较高；超压峰值较大时（150 kPa < Δp₀ ≤ 400 kPa），BDL模型的预测精度较低。原因有二：(1) 从模型学习的角度来看，样本点超压峰值较大时，样本组数据较少，模型难以充分学习以获得良好的预测能力；(2) 从爆炸工程的角度来看，较大的超压峰值通常发生在爆源的近场，而近场冲击波的传播规律还相对匮乏。总体而言，BDL模型在测试集上预测的超压峰值的相对误差约为12%，95%以上测试样本的预测误差小于50%。

表 7 基于目标范围的误差分解

Table 7. Error decomposition based on target range

超压峰值/kPa	样本数量	δ_mape/%	φ/%
超压峰值/kPa	样本数量	δ_mape/%	0 < δ_mape ≤ 10%	10% < δ_mape ≤ 20%	20% < δ_mape ≤ 30%	30% < δ_mape ≤ 50%	δ_mape ≤ 50%
[0, 50]	1927	10.173	66.27	21.33	7.16	3.68	98.44
(50, 100]	282	20.492	37.94	24.82	15.60	16.67	95.03
(100, 150]	124	17.714	41.13	20.97	18.55	14.52	95.15
(150, 400]	83	25.008	18.07	18.07	26.51	32.53	95.18

下载: 导出CSV

| 显示表格

3.4.3 基于比例距离的误差分解

BDL模型对冲击波传播的预测性能受到爆源当量Q和爆心距R的影响。比例距离d与当量Q及爆心距R的关系为：

$d=\frac{R}{\sqrt[3]{Q}}$

(13)

将2400个测试样本根据d的取值区间分成8组，表8给出了BDL模型对超压峰值的预测性能与比例距离之间的关系，其中φ_d为不同的d区间内的样本数量占该组测试子集样本总数的百分比。比例距离较小时，BDL模型预测的超压峰值的精度较低，δ_mape约为20%；比例距离较大时 ( $7\;\mathrm{m}/\mathrm{ kg}^{\frac{1}{3}} {\text{＜}} d{\text{≤}} 10\;\mathrm{m}/\mathrm{ kg}^{\frac{1}{3}}$ )时，BDL模型预测的超压峰值的精度较高，δ_mape约为8%。其原因可能是，较小的比例距离意味着更近的爆心距和更大的爆源当量，而近场冲击波的传播规律还相对匮乏。总体而言，BDL模型在测试集上预测的峰值超压的相对误差约为12%，93%以上测试样本的预测误差小于50%。

表 8 基于比例距离的误差分解

Table 8. Error decomposition based on scaled distance

比例距离/ $(\mathrm{m}\mathrm{\cdot kg}^{-\frac{1}{3}})$	样本数量	δ_mape/%	φ_d/%
比例距离/ $(\mathrm{m}\mathrm{\cdot kg}^{-\frac{1}{3}})$	样本数量	δ_mape/%	0 < δ_mape ≤ 10%	10% < δ_mape ≤ 20%	20% < δ_mape ≤ 30%	30% < δ_mape ≤ 50%	δ_mape ≤ 50%
(3, 5]	388	21.114	34.79	23.97	16.75	18.04	93.55
(5, 7]	963	13.091	57.22	23.88	8.72	7.79	97.61
(7, 10]	949	8.389	71.97	18.12	6.43	2.85	99.37

下载: 导出CSV

| 显示表格

综上所述，超压峰值较小、比例距离较大时，BDL模型能给出较为准确的预测结果，有必要进行误差分解。

4. 结　论

结合贝叶斯理论和深度学习，构造了复杂的目标结构。构造的场景考虑了30种爆源位置、30种当量以及10种起爆速度。基于blastFoam软件，获得了120组针对目标结构的不同爆炸情况下关键位置处的超压峰值数据集。在此基础上，采用基于变分推断的贝叶斯深度学习方法，构建了爆炸载荷估计模型，以爆源位置、当量以及测点的位置和比例距离等作为输入，以超压峰值作为输出，对复杂结构的爆炸载荷进行了快速估计，量化了估计结果的不确定性，得出以下主要结论。

(1) 在测试集上，所构建的贝叶斯深度学习模型预测的MAPE约为12%，在训练集上，模型预测的MAPE约为10%，该模型能较准确地估计复杂目标结构上的爆炸载荷。

(2) 在训练集和测试集上，模型的PICP值都稳定在81%以上，NMPIW值保持在5.7%以下，模型的估计值具有较高的准确性。所构建的不确定性量化方法为评价估计结果的可信性提供了一种新途径。

(3) 对于单组测试样本，模型对超压峰值的估计时间约为4 s；对于单个样本点，模型对超压峰值的估计时间约为20 ms，与数值模拟方法相比，所构建的贝叶斯深度学习模型对超压峰值估算有明显的速度优势。

所构建的基于贝叶斯深度学习的复杂结构爆炸载荷的估计方法不仅可以较快地估计目标结构上关键位置处的超压峰值，还可以量化估计结果的不确定性，是实现复杂结构爆炸载荷快速、可信估计的新方法。

图 1 背景网格与粒子单元

Figure 1. Background meshes and particles

下载: 全尺寸图片幻灯片

图 2 爆炸焊接布置

Figure 2. Distribution pattern of explosive welding

下载: 全尺寸图片幻灯片

图 3 爆炸焊接前处理模型

Figure 3. The preprocessing models for explosive welding

下载: 全尺寸图片幻灯片

图 4 爆炸焊接全过程的模拟

Figure 4. Simulation of the whole process of explosive welding

下载: 全尺寸图片幻灯片

图 5 金属板材的有效塑性应变

Figure 5. Effective plastic strains of the metal plates

下载: 全尺寸图片幻灯片

图 6 金属板材的碰撞速度

Figure 6. Impact velocities of the metal plates

下载: 全尺寸图片幻灯片

图 7 爆炸焊接的304L/Q235B复合板

Figure 7. An explosive-welded 304L/Q235B plate

下载: 全尺寸图片幻灯片

图 8 爆炸焊接的304L/Q235B复合板边缘

Figure 8. Edge of the explosive-welded 304L/Q235B plate

下载: 全尺寸图片幻灯片

图 9 复板的碰撞速度曲线

Figure 9. Impact velocity curves of the clad plate

下载: 全尺寸图片幻灯片

表 1 304L/Q235B材料模型参数^[14-15]

Table 1. Material parameters for 304L/Q235B^[14-15]

材料	ρ/(kg·m⁻³)	E/GPa	v	A/MPa	B/MPa	n	C	$\dot{{\varepsilon }_{0}}$ /s⁻¹
Q235B	7 850	210	0.3	380	275	0.36	0.022	0.001
304L	7 930	200	0.3	270	350	0.65	0.070	0.001

下载: 导出CSV

[1]	邵丙璜, 张凯. 爆炸焊接原理及工程应用 [M]. 大连: 大连工学院出版社, 1987: 1−5, 22−33.
[2]	CORIGLIANO P, CRUPI V, GUGLIELMINO E. Non linear finite element simulation of explosive welded joints of dissimilar metals for shipbuilding applications [J]. Ocean Engineering, 2018, 160: 346–353. DOI: 10.1016/j.oceaneng.2018.04.070.
[3]	GRIGNON F, BENSON D, VECCHIO K S, et al. Explosive welding of aluminum to aluminum: analysis, computations and experiments [J]. International Journal of Impact Engineering, 2004, 30(10): 1333–1351. DOI: 10.1016/j.ijimpeng.2003.09.049.
[4]	王宇新, 杨文彬, 李晓杰. 爆炸焊接二维复板飞行姿态计算机仿真 [J]. 爆破器材, 1999, 28(5): 1–4. DOI: 10.3969/j.issn.1001-8352.1999.05.001. WANG Y X, YANG W B, LI X J. Computerized simulation of two dimensional flying attitude of fly plate in explosive welding [J]. Explosive Materials, 1999, 28(5): 1–4. DOI: 10.3969/j.issn.1001-8352.1999.05.001.
[5]	王宇新, 李晓杰, 闫鸿浩, 等. 爆炸焊接CAE软件开发及工程应用 [J]. 工程爆破, 2018, 24(1): 1–7, 26. DOI: 10.3969/j.issn.1006-7051.2018.01.001. WANG Y X, LI X J, YAN H H, et al. CAE software development and engineering application of explosive welding [J]. Engineering Blasting, 2018, 24(1): 1–7, 26. DOI: 10.3969/j.issn.1006-7051.2018.01.001.
[6]	ZHANG Z L, LIU M B. Numerical studies on explosive welding with ANFO by using a density adaptive SPH method [J]. Journal of Manufacturing Processes, 2019, 41: 208–220. DOI: 10.1016/j.jmapro.2019.03.039.
[7]	SULSKY D, CHEN Z, SCHREYER H L. A particle method for history-dependent materials [J]. Computer Methods in Applied Mechanics and Engineering, 1994, 118(1/2): 179–196. DOI: 10.1016/0045-7825(94)90112-0.
[8]	WIĘZCKOWSKI Z. The material point method in large strain engineering problems [J]. Computer Methods in Applied Mechanics and Engineering, 2004, 193(39/40/41): 4417–4438. DOI: 10.1016/j.cma.2004.01.035.
[9]	WANG Y X, BEOM H G, SUN M, et al. Numerical simulation of explosive welding using the material point method [J]. International Journal of Impact Engineering, 2011, 38(1): 51–60. DOI: 10.1016/j.ijimpeng.2010.08.003.
[10]	ZHANG D Z, ZOU Q S, VANDERHEYDEN W B, et al. Material point method applied to multiphase flows [J]. Journal of Computational Physics, 2008, 227(6): 3159–3173. DOI: 10.1016/j.jcp.2007.11.021.
[11]	CHEN Z, HU W, SHEN L, et al. An evaluation of the MPM for simulating dynamic failure with damage diffusion [J]. Engineering Fracture Mechanics, 2002, 69(17): 1873–1890. DOI: 10.1016/S0013-7944(02)00066-8.
[12]	HU W Q, CHEN Z. Model-based simulation of the synergistic effects of blast and fragmentation on a concrete wall using the MPM [J]. International Journal of Impact Engineering, 2006, 32(12): 2066–2096. DOI: 10.1016/j.ijimpeng.2005.05.004.
[13]	黄鹏, 张雄, 马上, 等. 基于OpenMP的三维显式物质点法并行化研究 [J]. 计算力学学报, 2010, 27(1): 21–27. DOI: 10.7511/jslx20101004. HUANG P, ZHANG X, MA S, et al. Parallelization of 3D explicit material point method using OpenMP [J]. Chinese Journal of Computational Mechanics, 2010, 27(1): 21–27. DOI: 10.7511/jslx20101004.
[14]	JOHNSON G R, COOK W H. Fracture Characteristics of three metals subjected to various strains, strain rates, temperatures and pressures [J]. Engineering Fracture Mechanics, 1985, 21(1): 31–48. DOI: 10.1016/0013-7944(85)90052-9.
[15]	王宇新, 李晓杰, 王小红, 等. 爆炸焊接界面波物质点法三维数值模拟 [J]. 爆炸与冲击, 2014, 34(6): 716–722. DOI: 10.11883/1001-1455(2014)06-0716-07. WANG Y X, LI X J, WANG X H, et al. Numerical simulation on interfacial wave formation in explosive welding using material point method [J]. Explosion and Shock Waves, 2014, 34(6): 716–722. DOI: 10.11883/1001-1455(2014)06-0716-07.