钢箱内部爆炸破坏的SPH 数值模拟*

强洪夫，孙新亚，王 广，黄拳章

（火箭军工程大学，陕西 西安 710025）

桥梁是公路、铁路等交通网的咽喉。为了减小恒载，对斜拉桥、悬索桥等大跨度桥广泛采用钢箱梁作其结构支撑的主梁。随着恐怖袭击的不断演化，战略性桥梁成为恐怖主义爆炸破坏的重点目标，因此，有必要对钢箱结构内部爆炸过程中钢箱的变形规律进行研究[1-4]。Tang 等[5]运用LS-DYNA 软件研究了汽车炸弹爆炸冲击作用下桥塔、桥墩、桥面结构的局部破坏模式；姚术健等[2,6]通过实验，研究了钢箱内部爆炸过程中内壁失效变形的模式, 并利用ANSYS 软件分析了箱体的损伤特性。钢箱内部爆炸相比外部爆炸而言，爆炸冲击波在约束空间内多次反射、叠加、汇聚，该物理过程相当复杂，目前的研究大多局限于爆炸过程中钢箱变形的宏观现象。光滑粒子流体动力学(smoothed particle hydrodynamics,SPH)作为一种无网格的粒子方法，因其在计算过程中可避免网格重分和算法耦合，被广泛用于研究爆炸冲击问题。自Monaghan 等[7]修正的人工黏度解决了SPH 计算结果的非物理震荡问题后，强洪夫等[8]、Qiang 等[9]和Liu 等[10]利用该方法对爆炸冲击波冲击钢板的过程成功地进行了数值模拟。本文中，拟利用自编的SPH 程序对钢箱内部爆炸过程进行数值模拟，并与实验对比，验证SPH 算法在模拟钢箱内部爆炸问题上的有效性；分析爆炸过程中钢箱表面压力、应力波的传播规律以及特定点处速度变化趋势；得到钢箱表面的von Mises 应力分布；进一步对钢箱内不同位置爆炸的过程进行数值模拟，得到不同工况下钢箱的失效变形模式，以期通过研究为钢箱类结构的爆炸毁伤提供一种评估方法。

1 模型与算法

1.1 控制方程和状态方程

1.1.1 控制方程

钢箱的内部爆炸是绝热过程，外力相对于内部驱动力可忽略不计，即可不考虑重力的影响，在流体控制方程中忽略黏性项。因此，其控制方程组表述为：

式中：ρ、v、p、e、x 分别为密度、速度、压强、能量、位移。

1.1.2 状态方程

数值模拟过程中，TNT 炸药的爆轰气体压力通过标准的JWL 状态方程[11]获得：

对于钢箱材料，选用Q235 钢。爆炸过程中钢箱的各向同性压力p 利用Mie-Grüneisen 状态方程[12]求得:

1.2 本构方程

钢箱材料模型采用Børvik 等[13]、Nilson[14]提出的修正Johnson-Cook 本构模型，该模型通过流动应力和失效应变描述金属材料的大变形、高应变率等动态力学特性。其流动应力的表达式为：

1.3 控制方程的SPH 离散

SPH 算法是模拟流体流动的一种拉格朗日型粒子方法，通过使用一系列任意分布的粒子来求解具有各种边界条件的积分方程或偏微分方程。SPH 方法通常通过核函数插值实现场变量或场变量梯度的插值，通过粒子近似实现对核函数估计积分表达式的粒子离散[9]。

本文中利用SPH 方法模拟钢箱内部爆炸过程属于该方法在具有大变形和高度非均匀的动力学极端情形下的应用。为了解决爆炸域内求解结果的非物理振荡，且将冲击面内的动能以热能的形式耗散，引入人工黏性[7]，同时为了更好地处理爆炸过程中钢箱大变形、大扭曲等引起的密度和光滑长度变化剧烈的问题，采用Qiang 等[9]提出的完全变光滑长度SPH 方法对钢箱内部爆炸过程进行数值模拟，结合后的方程组如下：

式中： fi为 修正系数； Wij=W(xi-xj,h)为核函数，一般选用三次样条函数；Πij为人工黏度；N 为求解区域内粒子总数；α、β 表示空间坐标轴方向； h 为光滑长度，在连续性方程中将其对时间求导，得：

式中：d 为空间维数。

2 钢箱内部爆炸过程的数值模拟

2.1 数值模型

为了证明SPH 方法在数值模拟钢箱结构内部爆炸过程时的有效性，本文中设定了与姚术健等[2]的实验一致的工况条件，如图1(a)所示，钢箱模型为立方体形状，边长为600 mm，且采用了120 mm 的延伸结构，Q235 钢板厚度为4 mm；TNT 炸药为立方体形状，边长为40 mm。具体粒子配置如图1（b）所示，粒子间距为4 mm，炸药粒子数为1 331，钢箱粒子数为265 600。

图 1 钢箱和TNT 炸药试件模型图Fig.1 Model diagrams for the steel box and TNT dynamite specimen

为保证箱体各壁面受力均衡，TNT 炸药位于钢箱内部正中心，且各边与箱体对应边平行，如图2 所示。计算中，TNT 炸药密度为1 500 kg/m3, 采用点起爆方式，起爆点设置在炸药上表面中心点（-73.9，5.6，50.7），爆速为6 930 m/s，初始内能为4 290 kJ/kg，状态方程选用JWL 状态方程；钢箱密度为7 800 kg/m3,使用Johnson-Cook(J-C)本构模型，具体材料参数A=369 MPa, B=684 MPa, n=0.73, C=0.083, m=1.7,Tm=775 K,=0.1 s-1, CV=875 J/(kg·K), D1=0.13, D2=0.13, D3=-1.5, D4=0.011, D5=0。在组成钢箱的各板之间施加焊接锚固作用，由文献[15]可知焊缝金属强度不应低于母材。为简化计算，忽略焊缝附近区域母板材料性能的变化，将焊缝与母材视为同一结构，即认为钢箱是一体的。

图 2 钢箱内部爆炸算例模型Fig.2 The case model for the explosion inside the steel box

2.2 钢箱表面变形情况

数值模拟中，核函数选用三次样条函数，光滑长度取1.5 倍粒子间距，时间积分采用蛙跳格式，时间步长取0.1 μs。钢箱内部爆炸的数值模拟中，钢箱表面的变形情况如图3 所示。起始时刻，钢箱表面无变形，1.97 ms 时中心处产生明显的向外鼓起，并向四周传播变形；3.34 ms 时变形首次到达表面边沿处，且表面角隅处发生变形；12.03 ms时，表面变形达到最大，并保持不变。

图4 给出了钢箱内部爆炸实验和数值模拟的对比，图中A 处钢板发生卷边破坏，B 处钢板向外鼓起，C 处钢板发生扭曲变形，通过对比发现，数值模拟和实验所得钢箱壁面失效变形模式吻合较好。图5 给出了钢箱内部爆炸过程中箱体壁面挠度的变化情况。从图5 可以看出，85 μs 时钢箱壁面开始变形，挠度值随时间逐渐增大，12.03 ms 后基本保持不变。数值模拟得到的钢箱壁面最终挠度值为31.80 mm，而实验值为32.93 mm，实验和数值模拟的相对误差为3.5%。产生这一误差的主要原因可能是：(1)炸药理论设计的质量为96 g，由于密度不均匀等原因实验中炸药的真实质量却为98.443 g；(2)本文数值模拟过程中Q235 钢材料密度是均匀的，没有考虑制作中的工艺误差。综上所述，可得SPH方法能够数值模拟钢箱内部爆炸过程中钢箱壁面的变形情况。

图 3 钢箱内部爆炸过程中不同时刻箱体内壁的变形情况Fig.3 Its inner wall deformation at different times for the steel box subjected to inner explosion loading

图 4 钢箱壁面破坏的数值模拟和实验对比Fig.4 Comparison of steel box wall damage between numerical simulation and experiment

图 5 钢箱壁面挠度变化趋势Fig.5 Deflection in the center of the steel box wall varying with time

2.3 钢箱表面压力和von Mises 应力分布

图 6 钢箱内部爆炸过程中，壁面von Mises 应力分布Fig.6 Von Mises stress inside the steel wall during the explosion in the steel box

图 7 钢箱内部爆炸过程中，壁面的压力分布Fig.7 Pressure distribution of the wall during the explosion in the steel box

装药在完全封闭的钢箱内部爆炸时，爆炸冲击波在钢箱四周壁面上发生多次反射，冲击波强度逐渐降低。图6 给出了钢箱内部爆炸时，钢箱表面的von Mises 应力分布；图7 给出了钢箱内部爆炸过程中，钢箱表面压力分布情况。从图6 和图7 可得：钢箱内部炸药爆炸后，很快产生爆炸冲击波并迅速（85.0 μs）到达钢箱表面，爆心投影点附近出现压强和应力最大值；随后冲击波沿着钢箱的内表面向板与板的连接处（角隅位置）运动，同时爆心投影点处超压峰值逐渐降低，436.1 μs 时各板向角隅运动的冲击波在角隅处汇聚，形成较强的冲击波压力，并出现应力集中和应力增强现象；随后各板反射的冲击波与角隅处汇聚加强了的冲击波向箱室中心汇聚，使箱室中心形成较高超压，箱室中心超压又会向四周传播，此过程与起爆后传播过程类似，钢箱表面承受多次类似的冲击载荷作用后，钢箱表面压强和von Mises 应力逐渐消失，且箱内维持较长时间的余压。

2.4 钢箱表面中心点处速度和压强

图8 给出了钢箱内部爆炸过程中，钢箱表面中心点处速度的变化趋势；图9 给出了钢箱内部爆炸过程中，钢箱表面中心点处压强的变化趋势。从图8 可以看出：101.2 μs 时（如点A 所示），表面中心点的速度跳跃式增加至22.3 m/s，其后速度值又迅速减小，最终减小至0，且速度值在B 区域内出现上下微小浮动。从图9 可以看出：85.0 μs 时，钢箱表面中心点处的压强迅速增大，达到峰值487 MPa； 131.6 μs 时表面压强达到第二次峰值394 MPa；随后压强值上下扰动并逐渐减弱，出现扰动的主要原因是爆炸冲击波在钢箱表面上的来回反射；2 000.0 μs 时压强值逐渐减弱直至消失。由文献[16]可得：压力和速度时程曲线形状与实验及经典冲击波理论相吻合，即初次峰值较大，然后发生多次反射，速度和压强逐渐减弱。

图 8 钢箱壁面中心点处的速度变化趋势Fig.8 Velocity varying with time at the wall center of the steel box

图 9 钢箱壁面中心点处的压强变化趋势Fig.9 Pressure varying with time at the wall center of the steel box

3 钢箱内部不同位置爆炸失效模式

为了进一步探讨同一爆炸当量下，钢箱内部不同位置爆炸时，钢箱结构的破坏程度。本文设定了以下3 种计算工况（如图10 所示）：工况1，炸药位于钢箱底部正中心；工况2，炸药在底部中心正上方60 mm 处；工况3，炸药位于底板与侧板之间的角隅处，距离各板60 mm。

图11 给出了3 种工况下钢箱表面的失效变形模式。图11(a)为工况1 钢箱的变形情况，钢箱底板局部特征破坏明显，产生圆形的冲切破坏，断裂口处相对比较平滑，且底板产生裂纹（A）；同时底板与侧板的角隅处产生了长条枣核形状的裂口（B）。图11(b)为工况2 钢箱的变形情况，钢箱底板剧烈破坏，产生不规则的撕裂口，断裂口呈锯齿状，且底板没有产生比较规则的裂纹（A）；同时底板与侧板的角隅处也产生了长条状枣核形状的裂口，裂口周围向外凸起（B）。图11(c)为工况3 钢箱的变形情况，距离炸药较近的角隅处发生严重破坏，断裂口周围钢板发生扭曲变形（A）；同时发现，在远离炸药的3 个钢板的角隅处发生了撕裂破坏（B）。结合各工况下钢箱失效变形的情况，可得炸药在钢箱角隅处爆炸时破坏最为严重。

图 10 钢箱中炸药位置示意图(a=60 mm)Fig.10 The explosive location diagram in the steel box(a=60 mm)

图 11 3 种工况下，钢箱表面局部破坏示意图Fig.11 Surface partial destruction diagrams of the steel box under three different conditions

4 结 论

采用光滑粒子流体动力学(smoothed particle hydrodynamics, SPH)方法对钢箱内部爆炸的过程进行数值模拟，研究了爆炸过程中钢箱表面的变形，分析了相同当量炸药在钢箱内部不同位置爆炸时钢箱表面的失效模式。

（1）通过比较实验与SPH 方法数值模拟获得的结果，数值模拟与实验所得的失效变形模式吻合较好，实验和数值模拟得到的钢箱表面挠度值相对误差为3.5%，验证了本文SPH 方法在模拟钢箱内部爆炸问题上的有效性；

（2）采用SPH 方法数值模拟钢箱内部爆炸的过程，得出与实验一致的变形结果。同时也发现，在钢箱内部爆炸过程中，钢箱角隅处出现超压和应力集中现象，这与实验吻合较好。

（3）进一步研究了相同当量炸药在钢箱内部不同位置爆炸时钢箱失效破坏情况，发现炸药在钢箱角隅处爆炸时钢箱破坏最严重。

以上研究结果可为工事结构的防护设计提供参考。