基于PCA-BPNN的桥梁爆炸荷载时程预测

杜晓庆; 何益平; 邱涛; 程帅; 张德志

doi:10.11883/bzycj-2023-0343

基于PCA-BPNN的桥梁爆炸荷载时程预测

doi: 10.11883/bzycj-2023-0343

杜晓庆^1,,
何益平¹,
邱涛^1, ,,
程帅²,
张德志²

1.
上海大学力学与工程科学学院，上海 200444
2.
西北核技术研究所，陕西西安 710024

详细信息

作者简介:
杜晓庆（1973－　），男，博士，教授，dxq@shu.edu.cn

通讯作者:
邱　涛（1994－　），男，博士，qiutao@shu.edu.cn

中图分类号: O383
计量
- 文章访问数: 402
- HTML全文浏览量: 95
- PDF下载量: 134
- 被引次数: 4
出版历程
- 收稿日期: 2023-09-26
- 修回日期: 2024-04-19
- 网络出版日期: 2024-04-19
- 刊出日期: 2025-03-05

Prediction of blast loads on bridge girders based on PCA-BPNN

1.
School of Mechanics and Engineering Science, Shanghai University, Shanghai 200444, China
2.
Northwest Institute of Nuclear Technology, Xi’an 710024, Shaanxi, China

摘要

摘要: 人工智能方法是预测爆炸荷载的新手段，但现有方法主要用于预测爆炸冲击波的超压峰值或冲量，而用于预测反射超压时程的研究不多。针对这一问题，以平面冲击波绕射桥梁主梁为对象，提出了一种基于主成分分析（principal components analysis, PCA）和误差反向传播神经网络（backpropagation neural network, BPNN）的桥梁爆炸冲击波反射超压时程预测模型。该预测模型利用PCA降维处理时程数据，基于多任务学习的BPNN算法，提出了考虑超压峰值和冲量峰值影响的损失函数，使模型能有效预测不同入射超压下的桥梁冲击波荷载时程。通过分析多任务学习模型、多输入单输出模型和多输入多输出模型等3种BPNN模型，发现多任务学习模型的预测精度最高，而多输入多输出模型难以有效适应当前预测任务需求。采用多任务学习模型预测得到的桥梁表面各测点位置的反射超压时程、超压峰值精度较高，决定系数R²分别为0.792和0.987，作用在箱梁上的合力时程和扭矩时程预测值也与数值模拟值较为吻合。同时，该模型对内插值预测的表现优于外推值预测，但其在预测外推值方面同样展现出了一定的能力。
- 爆炸荷载预测 /
- 反射超压时程 /
- 误差反向传播神经网络 /
- 主成分分析 /
- 多任务学习
Abstract: Facing the challenges on the accurate and effective prediction under extreme loads, machine learning has gradually demonstrated its potential to replace traditional methods. Existing approaches primarily focus on predicting the peak overpressure or impulse of explosive shock waves, with limited research on predicting the reflected overpressure time history. Load-time history prediction encompasses not only the peak overpressure but also embraces various multi-dimensional information including duration, waveform, and impulse, thereby offering a more comprehensive depiction of the dynamic temporal and spatial characteristics of shock waves. To address this issue, a prediction model for bridge surface reflected overpressure time history is proposed, targeting a planar shock wave diffracting around a bridge section. This model is based on principal component analysis (PCA) and back propagation neural network (BPNN) algorithm with multi-task learning. A loss function considering the impact of peak overpressure and maximum impulse is introduced to fully consider the potential correlations between different modes after PCA dimension reduction. This enables the model to effectively predict bridge shock wave load time histories under varying incident overpressure. Through the analysis of three types of BPNN models, multi-task learning model, multi-input single-output model, and multi-input multi-output model. It was found that the multitask learning model has the highest prediction accuracy, while the multi-input multi-output model struggles to effectively adapt to the current predictive task requirements. The multitask learning model, used for predicting, achieves high precision in forecasting the time history of reflected overpressure at various measurement points on the bridge surface and the peak overpressure values, with R² values of 0.792 and 0.987. It also closely matches the simulation values in predicting the time history of combined forces and torque acting on the box girder. Additionally, this model performs better in interpolative value prediction than in extrapolative value prediction, but it also demonstrates a certain capability in predicting extrapolative values.
- explosion load prediction /
- reflected overpressure time history /
- back propagation neural network /
- principal components analysis /
- multi-task learning

HTML全文

预测岩石类材料在动态条件下的力学响应及破坏行为一直是数值计算中的重难点，建立合适的动态力学模型并选取适当的材料模型参数十分重要。目前常用的岩石类材料动态模型包括：Holmquist-Johnson-Cook模型（HJC模型）^[1]、Karagozian-and-Case模型（K&C模型）^[2]、Riedel-Hiermaier-Thoma模型（RHT模型）^[3]、Johnson-Holmquist 模型（JH-2 模型）^[4]、Talyor-Chen-Kuszmaul 模型（TCK 模型）^[5]等，并被应用于如冲击、爆炸、侵彻等多种动态过程的仿真分析中^[6-8]。研究表明，这些模型在运用时并不都能得到较理想的效果^[9-10]，原因在于模型本身存在一定的不足，无法很好地反映岩石材料的某些特征^[11-12]。一些学者对现有模型进行了部分修正，并将改进后的模型应用于相关的动态分析中，获得了较好的结果^{[10, 13-16]}。但建立更为合理的、完善的岩石类材料动态力学模型仍然十分必要。

Kong等^[11]基于塑性损伤理论建立了一种新的混凝土动态力学模型（Kong-Fang流体弹塑性损伤材料模型，KF模型），该模型定义了基于3个应力不变量的屈服面，并对压缩和拉伸损伤计算进行了区分。模型中引入了部分关联流动法则，可以较好地描述混凝土材料的剪胀行为。通过多种单元数值试验以及爆炸、冲击等动力学实例的模拟分析，对比HJC、K&C及RHT模型的模拟结果，验证了KF模型的适用性。Huang等^[17]进一步对KF模型进行了修正以适用于岩石材料，他们利用大量硬岩实验数据标定了模型中的部分参数，同时采用简明的半经验公式来描述动态增长因子随应变率的变化规律。通过多单元拉压数值试验及岩石劈裂、冲击实例的仿真计算，验证了修正模型对岩石的动态响应具有较好的预测能力。

考虑到天然岩石的力学性质差异较大，本文中将以山东五莲地区的花岗岩对研究对象，采用KF模型，根据模型参数定义，一方面，通过准静态单轴压缩、劈裂、常规三轴及动态分离式霍普金森杆压缩实验对模型中的强度参数进行确定，同时利用动态巴西圆盘实验结果验证应变率相关参数的有效性。另一方面，通过平板撞击实验拟合模型中的状态方程参数；最后结合实测参数值，利用KF模型对花岗岩侵彻试验过程进行仿真，通过对比计算得到的靶体侵彻深度和弹坑尺寸与实际试验的结果，验证材料模型及参数值的适用性。

1. KF模型简介及参数分类

1.1 KF模型简介

KF模型是一种可以描述岩石类材料在高压和高应变率条件下动态响应的计算本构模型^[11]。模型主要对屈服面、应变率效应和损伤累积进行了定义，此外，还采用了部分关联流动法则和多项式状态方程来分别描述岩石在低压压缩下的剪胀性和高压压缩下的变形特征。

1.1.1 基于三个应力不变量的屈服面

KF模型采用了压缩子午线的最大强度面σ_m和残余强度面σ_r，其表达式分别为^[11]：

${\sigma _{\text{m}}} = \left\{ {\begin{array}{*{20}{l}} 3\left[ {p/\left( {1 - D} \right) + T} \right] & p \text{≤} 0 \\ 1.5\left( {p + T} \right)/\psi & 0 \text{＜} p \text{＜} {f_{\text{c}}}/3 \\ {f_{\text{c}}} + \left( {p - {f_{\text{c}}}/3} \right)/\left( {{a_1} + {a_2}p} \right) & p \text{≥} {f_{\text{c}}}/3 \end{array}} \right.$

(1)

${\sigma _{\text{r}}} = p/\left( {{a_1} + {a_2}p} \right)$

(2)

式中：f_c和T分别为无侧限单轴压缩强度和拉伸强度；p为静水压力；D为总损伤；a₁、a₂为与强度面相关的材料常数；ψ为拉压子午线的比值，表示与压力相关的函数。最初，KF模型中定义ψ为p的分段函数，具体根据混凝土的三轴压缩和拉伸实验数据确定^[14]，后来Huang等^[17]考虑到硬岩的压缩性质与混凝土存在差异，对ψ进行了校准以适用于硬岩。

当材料并未达到最大强度面σ_m或残余强度面σ_r时，当前破坏面Y则依据当前岩石损伤情况对σ_m和σ_r进行插值确定^[11]：

$Y\left( {{\sigma _{ij}},D} \right) = \sqrt {3{J_2}} = r'\left[ {D\left( {{\sigma _{\text{r}}} - {\sigma _{\text{m}}}} \right) + {\sigma _{\text{m}}}} \right]$

(3)

式中：J₂为偏应力第二不变量；r′为当前子午线与压缩子午线的比值，其值可通过ψ(p)与Lode角计算获得^[11]，函数形式与RHT模型中的参数R₃(θ)（θ为Lode角）计算公式一致。

1.1.2 拉压分离的应变率效应

岩石材料动态强度具有明显的率效应^[18]，因此KF模型将当前破坏面Y进一步改进为：

$Y = \eta Y\left( {p/\eta } \right)$

(4)

式中：η为动态增长因子，其值为岩石动态强度与准静态强度的比值。

考虑到拉压条件下岩石将表现出不同的应变率效应，KF模型中分别定义了压缩强度动态增长因子η_c与拉伸强度动态增长因子η_t，由于原始计算公式较复杂，文献[17]中采用一种半经验的简洁公式对其进行了改进：

${\eta _{\text{c}}} = 1.0 + \alpha {{\text{e}}^{\beta \left( \lg {\dot \varepsilon } - \theta \right) }}$

(5)

${\eta _{\text{t}}} = 1.0 + ({\eta _{\text{c}}} - 1){{{f_{\text{c}}}} \mathord{\left/ {\vphantom {{{f_{\text{c}}}} T}} \right. } T}$

(6)

式中：α、β、θ为材料参数， $\dot \varepsilon$ 为应变率。

1.1.3 状态方程

当压力较高时，岩石材料的强度和结构的影响一般可以忽略^[19]，此时仅需考虑材料的体积变形，因此模型采用了单独的多项式状态方程来表征岩石在高压条件下的体积压缩行为：

$p = {k_1}\mu + {k_2}{\mu ^2} + {k_3}{\mu ^3}$

(7)

式中：k₁、k₂和k₃为状态方程材料参数，μ为体应变。

此外，模型中定义了材料损伤、变形及单元删除。（1）考虑到岩石材料在拉压荷载下的破坏机理不同，分别定义了拉伸损伤D_t和压缩损伤D_c，并根据D_t和D_c计算了岩石材料的总损伤D。在计算损伤时共定义了6个损伤累积相关参数，对于硬岩，文献[17]中给出了这类参数的建议值：d₁=0.04、d₂=2.3、γ=1、c₁=3、c₂=6.93、ε_frac=0.0015，因此，在本研究中也将这类参数作为常数处理。（2）为了表征岩石在低围压压缩荷载下的剪胀性，基于K&C模型中的塑性势函数，引入了部分关联的塑性流动法则，得到了相应的破坏面与应力更新计算公式。（3）模型中还采纳了两种单元删除准则，分别是基于损伤的拉伸删除准则和基于等效塑性应变的压缩删除准则，以表征岩石在遭受如爆炸等动载荷时材料前端和后端出现的成坑和剥落等破坏特征。由于损伤累积部分的参数为定值，而部分关联塑性流动法则与单元删除准则中没有涉及新的参数，因此本文中不再介绍。

1.2 模型参数分类及确定方法

为了获得待研究花岗岩材料的KF模型参数值，需确定各参数的实验方法。首先，根据模型定义可将待确定参数分为材料基本强度参数、强度相关参数和状态方程相关参数，见表1。材料基本强度参数f_c、E和ν的确定可通过无侧限单轴压缩实验进行；对于参数T的确定，由于直接拉伸实验费用高、成功率较低，而巴西劈裂实验制样方便、操作简单、成功率高^[20]，因此将采用劈裂实验获得的材料劈裂抗拉强度f_t来代替参数T；强度面相关的材料常数a₁、a₂的确定则需要大量三轴压缩实验数据；由式（5）～（6）可知，应变率相关的参数α、β、θ可通过SHPB实验确定，同时可利用SHPB-BD实验结果对拟合参数值的有效性进行验证；对状态方程相关参数的确定，可先通过平板撞击实验获得岩石的冲击雨果纽数据，然后拟合材料的压力-体应变关系曲线（p-μ曲线）获得。

表 1 KF模型参数的分类及实验确定方法

Table 1. Parameter classification and experimental determination method of KF model

分类	参数	物理意义	确定方法
材料基本强度参数	f_t	劈裂拉伸强度	准静态劈裂实验
	f_c	单轴压缩强度	准静态单轴压缩实验
	E	弹性模量
	υ	泊松比
强度面相关参数	a₁、a₂	强度面相关材料常数	准静态常规三轴实验
强度面相关参数	α、β、θ	应变率相关材料参数	SHPB/ SHPB-BD实验
状态方程相关参数	k₁、k₂、k₃	状态方程参数	平板撞击实验

下载: 导出CSV

| 显示表格

2. 准静态实验及相关参数确定

2.1 实验样品

选用产自山东五莲地区的粗中粒角闪二长花岗岩，岩石为二长半自形粒状结构，块状构造，主要矿物为斜长石(Pl)、钾长石(Kfs)、石英(Qtz)，共占矿物总量的90%。次要矿物为角闪石(Hbl)、黑云母(Bt)、辉石(Px)；蚀变及充填矿物为绢云母、绿帘石。岩石的实物照片和显微图像如图1所示。

图 1 岩石试件的实物照片和显微图像

Figure 1. Photographs and micrographs of rock samples

下载: 全尺寸图片幻灯片

根据相关的实验规范^[21]及实验机对试件的要求，制作各准静态实验的岩石试样：无侧限单轴压缩及常规三轴实验的试件为直径50 mm、高100 mm的标准圆柱形；劈裂实验试件为直径50 mm、高25 mm的圆盘。为了保证实验结果的可靠性，每种实验均开展了重复实验。

2.2 准静态单轴压缩及劈裂实验

分别利用TAW2000岩石三轴伺服压缩机及WDW-500E微机控制电子式万能试验机对花岗岩试件进行准静态单轴压缩和劈裂实验，如图2所示。实验中均采用变形控制加载，设置变形速率分别为0.06、0.03 mm/min，使得实验过程中试件的应变率为10⁻⁵ s⁻¹。

图 2 花岗岩单轴压缩及劈裂实验

Figure 2. Typical view of uniaxial compression test (UCT) and splitting test (ST)

下载: 全尺寸图片幻灯片

对单轴压缩实验获得的轴向压缩应力-应变曲线进行求导，其中弹性段导数即为试件的弹性模量E，同时采用f_c/2处的环向应变和轴向应变来计算各试件的泊松比υ^[22]。结合各试件的单轴抗压强度值f_c及劈裂拉伸强度f_t，利用两倍标准差法剔除异常值，并计算剩余有效数据的平均值作为花岗岩试件的KF模型参数值，结果见表2。

表 2 花岗岩的基本强度参数值

Table 2. Basic strength parameters of granite

试件	f_c /MPa	E/GPa	υ	试件	f_t /MPa
UCT-1	125.088	38.076	0.310	ST-1	7.714
UCT-2	133.206	41.82	0.352	ST-2	6.9
UCT-3	128.969	38.843	0.318	ST-3	8.12
UCT-4	133.802	39.534	0.498	ST-4	6.887
UCT-5	133.772	40.759	0.259	ST-5	7.409
KF模型参数值	130.967±3.856	39.806±1.497	0.347±0.091	KF模型参数值	7.406±0.532

下载: 导出CSV

| 显示表格

2.3 准静态常规三轴实验

利用MTS815型电液伺服岩石力学试验机对花岗岩试件进行准静态常规三轴实验，为了满足参数确定需求，设置了5种不同的围压水平：10、20、30、40、50 MPa。其中10、20 MPa工况各进行1次实验，30、40、50 MPa工况重复3次实验。如图3所示，将试件置于实验平台上，先施加一定的轴向压力防止试件滑动，然后再以1～3 MPa/min的加载速度施加围压至设定值。保持围压约5～10 min，确定不漏油后，以0.06 mm/min的加载速度施加轴压，直至试件破坏。

图 3 常规三轴实验

Figure 3. A view of triaxial compression test (TXC, σ₂=σ₃)

下载: 全尺寸图片幻灯片

花岗岩试件的常规三轴压缩应力-应变曲线如图4所示。由图4可知，围压实验的重复性较好，岩石的峰值强度随围压的增加而增加。由式（1）可知，确定参数a₁、a₂只需利用强度面应力值。因此，利用各试件的围压及相应的轴向峰值强度分别计算出各工况下的偏压力（σ_D=σ₁−σ₂）和静水压力（p=1/(3(σ₁+2σ₂))），结果见表3。进一步对a₁、a₂进行拟合，获得如图5所示的参数拟合值：a₁=0.36，a₂=0.09/f_c。

图 4 常规三轴实验花岗岩试件应力-应变曲线

Figure 4. Stress-strain curves of rock samples in TXC tests

下载: 全尺寸图片幻灯片

表 3 不同围压条件下试件的轴向峰值强度、偏压力与静水压力

Table 3. The strength, deviatoric pressure and hydrostatic pressure of rock samples under different confining pressures

试件	围压/MPa	轴向峰值强度/MPa	静水压力/MPa	偏应力/MPa
10-1	10	256.220	92.073	246.220
20-1	20	327.341	122.447	307.341
30-1	30	385.384	148.461	355.384
30-2		393.095	151.032	363.095
30-3		390.809	150.270	360.809
40-1	40	437.691	172.564	397.691
40-2		440.422	173.474	400.422
40-3		444.189	174.730	404.189
50-1	50	482.936	194.312	432.936
50-2		485.975	195.325	435.975
50-3		495.950	198.650	445.950

下载: 导出CSV

| 显示表格

图 5 花岗岩a₁、a₂参数拟合结果

Figure 5. Fitting results of parameters a₁ and a₂ of granite

下载: 全尺寸图片幻灯片

3. 动态实验及相关参数确定

3.1 SHPB及SHPB-BD实验

根据相关的实验标准^[23-24] ，基于50 mm杆径的分离式霍普金森杆实验系统开展SHPB及SHPB-BD实验。入射杆、透射杆分别长3.0、1.8 m，材料均为马氏体钢，两种实验采用的试件均为 $\varnothing$ 38 mm×19 mm的圆盘形试样。经过力平衡校验后，计算各有效实验的试件应变率 $\dot \varepsilon$ 、应变ε及应力σ：

$\dot \varepsilon {\text{ = }}{c_0}({\varepsilon _{\text{i}}} - {\varepsilon _{\text{r}}} - {\varepsilon _{\text{t}}})/{l_{\text{s}}}$

(8)

$\varepsilon {\text{ = (}}{c_0}{\text{/}}{l_{\text{s}}}{\text{)}}\int_0^{\text{t}} {({\varepsilon _{\text{i}}} - {\varepsilon _{\text{r}}} - {\varepsilon _{\text{t}}})} {\text{d}}t$

(9)

${\sigma _{\text{c}}}{\text{ = }}{A_0}{E_0}({\varepsilon _{\text{i}}} + {\varepsilon _{\text{r}}} + {\varepsilon _{\text{t}}})/(2{A_{\text{s}}})$

(10)

${\sigma _{\text{t}}}{\text{ = }}{A_0}{E_0}({\varepsilon _{\text{i}}} + {\varepsilon _{\text{r}}} + {\varepsilon _{\text{t}}})/({\text{π}}D_{\rm{s}}B)$

(11)

式中：c₀为实验杆的一维弹性波速，ε_i 、ε_r 、ε_t 分别为入射应变、反射应变和透射应变，l_s为试样长度，σ_c为SHPB实验中试样的压缩应力，σ_t为SHPB-BD实验中试样的拉伸应力，E₀为杆的弹性模量，A₀为杆的横截面积，A_s为试样的横截面积，D_s为试样直径，B为试样厚度。

为了确定并检验KF模型中的应变率相关参数α、β、θ，分别进行了14次动态压缩实验（剔除误差较大的一组数据）和15次动态劈裂实验，试件平均应变率为83～360 s⁻¹（动态压缩）和133～245 s⁻¹（动态劈裂）。实验所得的应力-应变曲线如图6所示，从图中可以发现，花岗岩具有一定的应变率效应。提取各实验中平均应变率及峰值强度并计算相应的η_c值与η_t值，结果见表4。利用表4中数据，结合式（5）～（6）对α、β、θ进行拟合及验证。图7中展示了拟合得到的参数值：α=0.06、β=3.47、θ=1.83。图8为利用SHPB-BD实验数据对参数进行验证的结果，由图8可知，拟合得到的参数值有较好的适用性。

图 6 花岗岩SHPB、SHPB-BD实验应力-应变曲线

Figure 6. Stress-strain curves of SHPB and SHPB-BD tests of granite samples

下载: 全尺寸图片幻灯片

表 4 SHPB、SHPB-BD实验的平均应变率及相应的

$\eta_{\rm c}$ 、

$\eta_{\rm t}$ 值

Table 4. Average strain rates and corresponding

$\eta_{\rm c}$ and

$\eta_{\rm t}$ values obtained from SHPB and SHPB-BD experiments

试件	动态压缩强度/MPa	η_c	平均应变率/s⁻¹	试件	动态劈裂拉伸强度/MPa	η_t	平均应变率/s⁻¹
SHPB-1	159.324	1.217	83	SHPB-BD-1	41.395	5.589	138
SHPB-2	142.925	1.091	105	SHPB-BD-2	37.363	5.045	135
SHPB-3	160.994	1.229	102	SHPB-BD-3	45.456	6.138	133
SHPB-4	156.455	1.195	85	SHPB-BD-4	45.548	6.150	135
SHPB-5	164.434	1.256	122	SHPB-BD-5	52.404	7.076	187
SHPB-6	180.403	1.377	144	SHPB-BD-6	52.654	7.110	183
SHPB-7	145.891	1.114	168	SHPB-BD-7	49.934	6.742	187
SHPB-8	140.933	1.076	187	SHPB-BD-8	49.225	6.647	188
SHPB-9	212.790	1.625	230	SHPB-BD-9	65.320	8.820	241
SHPB-10	260.696	1.991	325	SHPB-BD-10	54.190	7.317	245
SHPB-11	158.481	1.210	258	SHPB-BD-11	70.300	9.492	239
SHPB-12	177.977	1.359	292	SHPB-BD-12	60.436	8.160	243
SHPB-13	215.694	1.647	360	SHPB-BD-13	48.513	6.550	155
				SHPB-BD-14	61.372	8.287	199
				SHPB-BD-15	71.112	9.602	210

下载: 导出CSV

| 显示表格

图 7 基于SHPB实验数据拟合得到的α、β、θ参数值

Figure 7. Parameter values of α, β and θ fitted by SHPB tests data

下载: 全尺寸图片幻灯片

图 8 利用SHPB-BD实验数据对α、β、θ值进行验证的结果

Figure 8. Validation of α, β and θ values by SHPB-BD tests data

下载: 全尺寸图片幻灯片

3.2 平板撞击实验

为了获得更干净稳定的信号，基于57 mm口径一级轻气炮，采用反向撞击法开展实验^[25]。考虑到获取信号不能受到追赶稀疏波及边侧稀疏波的影响，飞片与靶板尺寸均设计为 $\varnothing$ 50 mm×10 mm，其中靶体材料为无氧铜。具体的实验原理为：当岩石飞片撞击无氧铜靶时，利用电探针测得飞片的速度w，利用DISAR（displacement interferometer system for any reflector）系统测得靶板背面的自由面质点速度u_p,Cu。根据应力波原理，结合冲击绝热线方程及Rankine-Hugoniot方程可知，岩石材料的波后质点速度u_p,rock、波后压力p_rock、岩石中的冲击波速度u_s,rock及岩石的体应变μ_rock^[26]分别为：

${u_{{\text{p}},{\text{rock}}}} = w - 0.5{u_{{\text{p}},}}_{{\text{Cu}}}$

(12)

${p_{{\text{rock}}}}{\text{ = }}0.5{u_{{\text{p}},{\text{Cu}}}}{\rho _{0,{\text{Cu}}}}\left( {{c_{0,{\text{Cu}}}} + 0.5{u_{{\text{p}},{\text{Cu}}}}{s_{{\text{Cu}}}}} \right)$

(13)

${u_{{\text{s}},{\text{rock}}}}{\text{ = }}{p_{{\text{rock}}}}/\left( {{\rho _{0,{\text{rock}}}}{u_{{\text{p}},{\text{rock}}}}} \right)$

(14)

${\mu _{{\text{rock}}}}{\text{ = }}{u_{{\text{p}},{\text{rock}}}}/{u_{{\text{s}},{\text{rock}}}}$

(15)

式中：ρ_0,Cu为铜的初始密度，ρ_0,Cu=8930 kg/m³；c_0,Cu、s_Cu为铜的材料参数， c_0,Cu=3940 m/s，s_Cu=1.49^[27]；ρ_0,rock为花岗岩的初始密度，其值通过密度实验获得：ρ_0,rock=2686 kg/m³。

实验时准备并安装好弹丸、电探针、靶板和DISAR系统。利用压缩空气（设计冲击速度小于400 m/s）或压缩氮气（设计冲击速度大于400 m/s）推动弹丸高速撞击靶板。为了拟合花岗岩的p_rock-μ_rock曲线，共设计了4种撞击力水平，并各重复3次实验。由式（12）～（15）获得各状态量的结果，见表5。基于实验数据拟合花岗岩试件的p_rock-μ_rock曲线，其形式如式（7），拟合结果如图9所示，拟合得到的p_rock-μ_rock方程为：

表 5 花岗岩平板撞击实验结果

Table 5. Plate impact test results of granite samples

实验	w/ (m∙s⁻¹)	u_p,Cu (m∙s⁻¹)	u_p,rock /(m∙s⁻¹)	p_rock /GPa	u_s,rock/ (m∙s⁻¹)	μ_rock
G-1	298.249	170.887	212.806	3.103	5430.281	0.039
G-2	298.586	171.614	212.779	3.117	5454.783	0.039
G-3	296.200	175.154	208.623	3.183	5681.906	0.037
G-4	364.353	210.768	258.969	3.856	5543.886	0.047
G-5	362.860	210.444	257.638	3.849	5563.631	0.046
G-6	370.399	217.718	261.540	3.988	5677.582	0.046
G-7	420.639	247.273	297.002	4.553	5708.849	0.052
G-8	425.834	247.173	302.248	4.552	5607.388	0.054
G-9	424.120	246.701	300.769	4.542	5623.706	0.053
G-10	536.010	308.458	381.771	5.743	5601.403	0.068
G-11	527.022	302.785	375.608	5.632	5582.955	0.067
G-12	536.027	312.243	370.879	5.817	5840.622	0.063

下载: 导出CSV

| 显示表格

图 9 拟合的花岗岩p_rock-μ_rock曲线

Figure 9. Fitting curve of plate impact data of granite samples

下载: 全尺寸图片幻灯片

${p_{{\text{rock}}}} = 45.002{\mu _{{\text{rock}}}} + 1413.751{\mu _{{\text{rock}}}}^2 - 12037.857{\mu _{{\text{rock}}}}^3$

(16)

式中：p_rock的单位为GPa。因此状态方程相关参数k₁、k₂、k₃的值分别为：k₁=45.002，k₂=1413.751，k₃= −12037.857。

综上所述，花岗岩的KF模型参数已被确定，见表6。

表 6 花岗岩KF模型参数拟合结果

Table 6. Fitting results of KF model parameters of granite

强度参数									状态方程参数
f_c/MPa	f_t/MPa	E/GPa	a₁	a₂/MPa ⁻¹	υ	α	β	θ	k₁/GPa	k₂/GPa	k₃/GPa
130.967	7.406	39.806	0.36	0.09/f_c	0.347	0.06	3.47	1.83	45.002	1413.751	−12037.857

下载: 导出CSV

| 显示表格

4. 花岗岩侵彻实验及仿真验证

为了验证实验确定的花岗岩KF模型参数值适用性，开展了花岗岩侵彻实验，并利用代入了实测参数值的KF模型对侵彻实验进行模拟，进一步对比仿真与实验得到的最终侵彻深度与成坑结果。

4.1 侵彻实验

花岗岩侵彻实验基于 $\varnothing$ 30 mm口径滑膛火炮开展，实验原理是利用火药气体燃烧产生的能量推动弹体侵彻岩石靶体。设计工况为子弹以约670 m/s的速度侵彻半无限厚靶体。

考虑到弹体发射的稳定性和岩石靶体取样大小的可操作性，本次实验采用次口径发射技术，利用直径20 mm、CRH为2.5的尖卵形头部杆形弹，外置直径30 mm的三瓣式弹托进行实验，弹体材料为超高强度合金钢30CrMnSiNi2A，力学性能指标及Johnson-Cook本构参数^[28]见表7。为了模拟半无限条件，需要合理设计靶体直径以忽略边界效应。当侵彻速度为670 m/s时，由已有的研究成果可知^[29]，若靶弹径比不小于30，则可以忽略有限边界效应。将满足需求的600 mm×600 mm×800 mm的花岗岩块放置在 $\varnothing$ 1200 mm×800 mm的钢套筒正中位置，并在空隙处浇筑C40混凝土振捣填实养护28 d制成靶体。

表 7 弹体材料参数^[28]

Table 7. Parameters of the projectile’s material^[28]

基本力学指标					JC本构模型参数
密度/ (kg∙m⁻³)	洛氏硬度 HRC	抗拉强度/ MPa	延伸率/ %	截面收缩率/ %	屈服强度A/ MPa	硬化系数B/ MPa	硬化指数n	应变率敏感系数c	温度软化系数m
7830	45	1380	12	58	1269	810	0.479	0.04	1

下载: 导出CSV

| 显示表格

实验系统现场布局如图10所示，整个系统主要包括3部分：发射平台、测速系统和高速摄影系统。其中发射平台为口径30 mm的滑膛火炮；测速系统用于获取子弹的侵彻速度，方法是通过在炮口和目标靶之间设置两道已知距离的测速靶纸，当弹体击穿靶纸时将分别输出2个电信号至测速仪，根据信号的时间差即可获得弹体的平均飞行速度。根据测得的飞行速度可进一步调整火药的填装量，以实现设计的发射初速。高速摄像系统用于观察子弹的脱壳情况，并判断弹体入射姿态是否满足垂直入射要求。

图 10 侵彻实验装置布局图

Figure 10. Device layout diagram of penetration test

下载: 全尺寸图片幻灯片

进行多次侵彻实验，根据高速摄影图片选取子弹侵彻着角小于5°的有效侵彻结果，图11为3个有效侵彻实验的弹坑照片与三维扫描结果。利用直尺和三维扫描仪对花岗岩侵彻后的靶体进行现场测量和光学扫描，测量及扫描得到的靶体侵彻深度及成坑参数见表8。

图 11 花岗岩靶体弹坑照片及三维扫描图

Figure 11. Photos and 3D scanning results of cratering failure on the impact surfaces of granite targets

下载: 全尺寸图片幻灯片

表 8 花岗岩靶成坑参数现场测量及三维扫描结果

Table 8. Direct measurement and 3D scanning results of crater parameters of granite targets

试验	子弹着角/(°)	侵彻速度/(m∙s⁻¹)	侵彻深度/mm		侵彻弹坑最大直径/mm		侵彻弹坑最小直径/mm
试验	子弹着角/(°)	侵彻速度/(m∙s⁻¹)	现场测量	三维扫描	现场测量	三维扫描	现场测量	三维扫描
靶1	−3.5	638	75.05	76.29	310	327.98	245	248.36
靶2	−3.5	667	72.19	72.48	430	408.24	275	283.00
靶3	−1.5	673	96.00	91.73	410	402.60	265	261.67

下载: 导出CSV

| 显示表格

4.2 仿真与实验结果对比

利用LS-DYNA对花岗岩侵彻实验进行仿真。考虑到实验的对称性，同时，为了减少计算成本，采用1/4对称模型进行仿真实验，网格设计如图12所示，将花岗岩靶体简化成了外套一圈薄层铁箍的 $\varnothing$ 1200 mm×800 mm的圆柱体，弹体则按实验中实用子弹进行建模，模型整体尺寸与实际实验中弹靶尺寸一致。

图 12 侵彻实验数值模型

Figure 12. Numerical model

下载: 全尺寸图片幻灯片

由于3组实验测得的子弹侵彻速度接近，因此模拟时设置弹体在靶体中心处以670 m/s的初速度垂直向下侵彻，弹体与靶体均采用实体单元网格，网格尺寸均设置为1 mm×1 mm×1 mm（弹头处网格做过渡处理）。靶体上表面设置为自由边界，下表面设置为无反射边界，靶体对称面上施加了对称边界。弹体和靶体间的接触采用*CONTACT_ERODING_SURFACE_TO_SURFACE接触算法定义。通过二次开发将KF模型嵌入LS-DYNA软件中，并作为靶体的材料模型，模型参数见表6。使用Johnson-Cook模型作为弹体的材料模型，具体参数见表7。

提取靶体单元的损伤图像并进行对称处理，得到的靶体损伤分布如图13所示。当损伤值为零时，表示单元未损伤，当损伤值为1时，表示单元完全损伤。由图13可知，计算得到的弹坑最大直径为400 mm，最小直径为300 mm。提取侵彻过程中子弹侵彻深度及加速度的历史曲线，如图14所示，可以发现，侵彻深度达到约80 mm时即保持不变，同时弹体的加速度时程曲线也呈现典型的三段式。

图 13 靶体损伤分布情况及成坑参数

Figure 13. Damage distribution and cratering parameters of target

下载: 全尺寸图片幻灯片

图 14 弹体的侵彻深度和加速度时程曲线

Figure 14. History curves of penetration depth and acceleration of projectile

下载: 全尺寸图片幻灯片

由于侵彻实验随机性较大，因此取成坑参数的平均值进行分析，结果见表9。对比实验和数值模拟得到的侵彻弹坑和侵彻深度发现，对于侵彻深度和弹坑最大直径，数值模拟与实验结果的误差小于5%；侵彻弹坑最小直径的模拟结果与实际结果误差较大，但仍小于15%。综上所述，采用KF材料模型，并利用实验获取的参数值得到的数值模拟结果与实验结果吻合较好，验证了材料模型及实测参数值具有较好的适用性。

表 9 实验和数值模拟得到的花岗岩靶侵彻结果比较

Table 9. Comparison of penetration results obtained from penetration test and numerical simulation

方法	侵彻深度/mm	侵彻弹坑最大直径/mm	侵彻弹坑最小直径/mm
实验	80.62±10.46	381.47±49.60	263.01±14.75
数值模拟	80.02	400.00	300.00
误差/%	−0.75	4.86	14.07

下载: 导出CSV

| 显示表格

5. 结　论

通过多种准静态及动态实验确定了山东五莲地区花岗岩的KF模型参数，并利用KF模型，结合实测参数值，对花岗岩侵彻实验进行了仿真计算，得到以下主要结论。

（1）通过准静态单轴压缩、劈裂及常规三轴实验确定了模型中的部分强度相关参数：f_c=130.967 MPa、f_t=7.406 MPa、E=39.806 GPa、υ=0.347、a₁=0.36、a₂=0.09/f_c。

（2）通过SHPB实验确定了模型中的应变率相关参数：α=0.06、β=3.47、θ=1.83，并利用SHPB-BD实验结果验证了参数的有效性。

（3）通过平板撞击实验确定了模型中的状态方程相关参数：k₁=45.002 GPa、k₂=1413.751 GPa、k₃=−12037.857 GPa。

（4）利用KF模型及实测参数值对花岗岩侵彻实验进行了模拟，结果表明，数值模拟获得的侵彻深度及成坑最大直径与实验结果的误差小于5%，侵彻弹坑最小直径的模拟结果与实际结果误差小于15%。总体上来说仿真结果与实际结果吻合较好，验证了材料模型及参数值的适用性。

图 1 主梁表面测点布置示意图

Figure 1. Schematic diagram of surface measurement points layout on the main beam

下载: 全尺寸图片幻灯片

图 2 计算域与边界条件

Figure 2. Computational domain and boundary conditions

下载: 全尺寸图片幻灯片

图 3 自由场超压时程曲线

Figure 3. Free-field overpressure time history curves

下载: 全尺寸图片幻灯片

图 4 PCA-BPNN训练和预测流程图

Figure 4. Process of training and predicting for PCA-BPNN

下载: 全尺寸图片幻灯片

图 5 累积方差解释率随主成分数量的变化

Figure 5. Variation of cumulative variance contribution with the number of principal components

下载: 全尺寸图片幻灯片

图 6 R²随主成分数量的变化

Figure 6. Variation of R² with the number of principal components

下载: 全尺寸图片幻灯片

图 7 多输出反向传播神经网络

Figure 7. Multi-output BPNN

下载: 全尺寸图片幻灯片

图 8 单输出反向传播神经网络

Figure 8. Single-output BPNN

下载: 全尺寸图片幻灯片

图 9 输入输出变量之间的相关矩阵

Figure 9. Correlation matrix between input and output

下载: 全尺寸图片幻灯片

图 10 K折交叉验证示意图

Figure 10. Schematic diagram of K-fold cross-validation

下载: 全尺寸图片幻灯片

图 11 入射超压0.3 MPa下测点的反射超压和冲量时程曲线

Figure 11. Reflected overpressure and impulse time histories at measurement points under 0.3 MPa incident overpressure

下载: 全尺寸图片幻灯片

图 12 入射超压0.8 MPa下测点的反射超压和冲量时程曲线

Figure 12. Reflected overpressure and impulse time histories at measurement points under 0.8 MPa incident overpressure

下载: 全尺寸图片幻灯片

图 13 入射超压1.5 MPa下测点的反射超压和冲量时程曲线

Figure 13. Reflected overpressure and impulse time histories at measurement points under 1.5 MPa incident overpressure

下载: 全尺寸图片幻灯片

图 14 入射超压0.3 MPa下箱梁的合力和扭矩时程曲线

Figure 14. Resultant force and torque time histories of box girder under 0.3 MPa incident overpressure

下载: 全尺寸图片幻灯片

图 15 入射超压0.8 MPa下箱梁的合力和扭矩时程曲线

Figure 15. Resultant force and torque time histories of box girder under 0.8 MPa incident overpressure

下载: 全尺寸图片幻灯片

图 16 入射超压1.5 MPa下箱梁的合力和扭矩时程曲线

Figure 16. Resultant force and torque time histories of box girder under 1.5 MPa incident overpressure

下载: 全尺寸图片幻灯片

图 17 不同入射超压下模型超压峰值预测拟合曲线

Figure 17. Model overpressure peak prediction fit graph under different incident overpressures

下载: 全尺寸图片幻灯片

图 18 超压峰值的数值模拟值与预测值对比

Figure 18. Comparison of peak overpressure between numerical simulation and predictive values

下载: 全尺寸图片幻灯片

表 1 超压峰值的各模型预测评价指标

Table 1. Evaluation metrics for peak overpressure prediction models

预测模型	RMSE	MAE	MAPE	R²
多输出模型	1.089	0.788	3.105	−0.386
单输出模型	0.152	0.092	0.220	0.976
多任务模型	0.115	0.074	0.180	0.987

下载: 导出CSV

表 2 冲量峰值的各模型预测评价指标

Table 2. Evaluation metrics for peak impulse prediction models

预测模型	RMSE	MAE	MAPE	R²
多输出模型	3.970	2.822	0.554	−0.403
单输出模型	0.866	0.646	0.133	0.934
多任务模型	0.693	0.542	0.109	0.958

下载: 导出CSV

表 3 超压时程曲线的各模型预测评价指标

Table 3. Evaluation metrics of overpressure time-history prediction models

预测模型	RMSE	MAE	MAPE	R²
多输出模型	0.173	0.068	−	−5.466
单输出模型	0.069	0.015	−	0.753
多任务模型	0.065	0.014	−	0.792

下载: 导出CSV

表 4 模型预测误差及经验公式计算误差对比

Table 4. Comparison of model prediction error and empirical formula calculation error

方法	p_i,max=0.3 MPa		p_i,max=0.8 MPa		p_i,max=1.5 MPa
方法	超压时程的R²	最大冲量的MAPE	超压时程的R²	最大冲量的MAPE	超压时程的R²	最大冲量的MAPE
PCA-BPNN	0.81	0.124	0.92	0.052	0.46	0.262
经验公式	−0.33	0.410	−2.05	0.573	−1.42	0.396

下载: 导出CSV

参考文献(27)

[1]	CLUBLEY S K. Non-linear long duration blast loading of cylindrical shell structures [J]. Engineering Structures, 2014, 59: 113–126. DOI: 10.1016/j.engstruct.2013.10.030.
[2]	DENNY J W, CLUBLEY S K. Long-duration blast loading & response of steel column sections at different angles of incidence [J]. Engineering Structures, 2019, 178: 331–342. DOI: 10.1016/j.engstruct.2018.10.019.
[3]	JIANG Y X, ZHANG B Y, WANG L, et al. Dynamic response of polyurea coated thin steel storage tank to long duration blast loadings [J]. Thin-Walled Structures, 2021, 163: 107747. DOI: 10.1016/j.tws.2021.107747.
[4]	SAKURAI A. Blast wave theory [M]. Mathematics Research Center, United States Army, University of Wisconsin, 1964.
[5]	RIGBY S E, AKINTARO O I, FULLER B J, et al. Predicting the response of plates subjected to near-field explosions using an energy equivalent impulse [J]. International Journal of Impact Engineering, 2019, 128(1): 24–36. DOI: 10.1016/j.ijimpeng.2019.01.014.
[6]	KINGERY C N, BULMASH G. Airblast parameters from TNT spherical air burst and hemispherical surface burst [R]. US Army Armament and Development Center, Ballistic Research Laboratory, 1984.
[7]	RIGBY S E, TYAS A, FAY S D, et al. Validation of semiempirical blast pressure predictions for far field explosions: is there inherent variability in blast wave parameters? [C]//Proceedings of the 6th International Conference on Protection of Structures Against Hazards. Sheffield, 2014.
[8]	LARCHER M, CASADEI F. Explosions in complex geometries: a comparison of several approaches [J]. International Journal of Protective Structures, 2010, 1(2): 169–195. DOI: 10.1260/2041-4196.1.2.169.
[9]	DENNIS A A, PANNELL J J, SMYL D J, et al. Prediction of blast loading in an internal environment using artificial neural networks [J]. International Journal of Protective Structures, 2021, 12(3): 287–314. DOI: 10.1177/2041419620970570.
[10]	HANSEN O R, HINZE P, ENGEL D, et al. Using computational fluid dynamics (CFD) for blast wave predictions [J]. Journal of Loss Prevention in the Process Industries, 2010, 23(6): 885–906. DOI: 10.1016/j.jlp.2010.07.005.
[11]	REMENNIKOV A M, ROSE T A. Modelling blast loads on buildings in complex city geometries [J]. Computers & Structures, 2005, 83(27): 2197–2205. DOI: 10.1016/j.compstruc.2005.04.003.
[12]	SOHAIMI A S M, RISBY M S. Using computational fluid dynamics (CFD) for blast wave propagation under structure [J]. Procedia Computer Science, 2016, 80: 1202–1211. DOI: 10.1016/j.procs.2016.05.463.
[13]	HASTIE T, TIBSHIRANI R, FRIEDMAN J H, et al. The elements of statistical learning: data mining, inference, and prediction [M]. New York: Springer, 2009.
[14]	FLOOD I, BEWICK B T, DINAN R J, et al. Modeling blast wave propagation using artificial neural network methods [J]. Advanced Engineering Informatics, 2009, 23(4): 418–423. DOI: 10.1016/j.aei.2009.06.005.
[15]	REMENNIKOV A M, ROSE T A. Predicting the effectiveness of blast wall barriers using neural networks [J]. International Journal of Impact Engineering, 2007, 34(12): 1907–1923. DOI: 10.1016/j.ijimpeng.2006.11.003.
[16]	BEWICK B, FLOOD I, CHEN Z. A neural-network model-based engineering tool for blast wall protection of structures [J]. International Journal of Protective Structures, 2011, 2(2): 159–176. DOI: 10.1260/2041-4196.2.2.159.
[17]	ALSHAMMARI O G, ISAAC O S, CLARKE S D, et al. Mitigation of blast loading through blast-obstacle interaction [J]. International Journal of Protective Structures, 2023, 14(3): 357–389. DOI: 10.1177/20414196221115869.
[18]	PANNELL J J, RIGBY S E, PANOUTSOS G. Physics-informed regularisation procedure in neural networks: an application in blast protection engineering [J]. International Journal of Protective Structures, 2022, 13(3): 555–578. DOI: 10.1177/20414196211073501.
[19]	PANNELL J J, RIGBY S E, PANOUTSOS G. Application of transfer learning for the prediction of blast impulse [J]. International Journal of Protective Structures, 2023, 14(2): 242–262. DOI: 10.177/20414196221096699.
[20]	LI Q L, WANG Y, SHAO Y D, et al. A comparative study on the most effective machine learning model for blast loading prediction: from GBDT to transformer [J]. Engineering Structures, 2023, 276: 115310. DOI: 10.1016/j.engstruct.2022.115310.
[21]	LI Q L, WANG Y, LI L, et al. Prediction of BLEVE loads on structures using machine learning and CFD [J]. Process Safety and Environmental Protection, 2023, 171: 914–925. DOI: 10.1016/j.psep.2023.02.008.
[22]	HUANG Y, ZHU S J, CHEN S W. Deep learning-driven super-resolution reconstruction of two-dimensional explosion pressure fields [J]. Journal of Building Engineering, 2023, 78: 107620. DOI: 10.1016/j.jobe.2023.107620.
[23]	孙传猛, 裴东兴, 陈嘉欣, 等. 基于深度学习的爆炸冲击波信号重构模型 [J]. 计测技术, 2022, 42(2): 57–67. DOI: 10.11823/j.issn.1674-5795.2022.02.07. SUN C M, PEI D X, CHEN J X, et al. Model for reconstruction of explosion shock wave signals based on deep learning [J]. Measurement Technology, 2022, 42(2): 57–67. DOI: 10.11823/j.issn.1674-5795.2022.02.07.
[24]	QIU T, CHENG S, DU X Q, et al. Spacing effects on blast loading characteristics of two tandem square columns under planar shock waves [J]. Physics of Fluids, 2023, 35(12): 127116. DOI: 10.1063/5.0177869.
[25]	KARLOS V, LARCHER M, SOLOMOS G. Analysis of the blast wave decay coefficient in the Friedlander equation using the Kingery-Bulmash data [R]. Joint Research Center, European Commission, 2015.
[26]	ZHOU Z H. Machine learning [M]. Springer Nature, 2021. DOI: 10.1007/978-981-15-1967-3.
[27]	MISRA D. Mish: a self-regularized non-monotonic activation function[R]. British Machine Vision Conference, 2019. DOI: 10.48550/arXiv.1908.08681.