接触爆炸荷载作用下带孔防护结构内冲击波的传播规律＊

范 进，徐大立，任新见

（1.南京理工大学土木工程系，江苏 南京 210094；2.总参工程兵科研三所，河南 洛阳471023）

防护结构是指能够抵抗预定杀伤武器破坏作用的工程结构，目前常见的野战阵地防护结构大多由钢筋、混凝土材料筑成，并设置了各类尺寸的孔口供射击、观察、通风和出入使用。当常规武器命中带孔防护结构后，接触爆炸产生的威力强大的冲击波经孔口进入防护结构内部，可能对结构内部的人员及设备造成一定程度的毁伤，甚至导致人员伤亡、设备瘫痪，使野战阵地丧失作战能力。而在现代战争中，带孔防护结构往往是精确制导武器打击的重点，其生存能力直接关系着战争的胜负。因此，关于常规武器对带孔防护结构毁伤效应的研究对防护工程的建设、综合作战效率的提高具有重大意义。

目前，已对自由场空爆、坑道内爆炸等条件下冲击波的传播规律进行了深入的研究，提出了诸如Henrych公式［1］、Sadovsk公式［2］、TM5－855－1公式［3］等具有广泛影响力的计算公式。C.R.Welch［4］在模型坑道中进行化爆实验，得到了堵口爆炸情况下坑道内空气冲击波冲量的经验公式；LS2000设计规范［5］中提出了堵口爆炸时坑道内空气冲击波冲量的确定方法；庞伟宾等［6］对爆炸冲击波在直坑道、T型坑道内的传播规律及其影响因素进行了原型坑道爆炸实验研究；黄广炎等［7］对炸药装药在腔室中心、墙角等不同约束情况下爆炸冲击波的传播特性进行了数值计算。

虽然对冲击波传播规律问题的研究已取得了一定成果，但既有研究均有其特定的适用范围［8］，尚不能准确反映接触爆炸荷载作用下带孔防护结构内冲击波的传播规律。随着精确制导技术的快速发展，战斗部直接命中防护结构并发生接触爆炸的可能性越来越大，这就要求对带孔防护结构内冲击波的传播规律有清晰的认识，以便准确评估爆炸冲击波对防护结构的毁伤效应。本文中，采用ATUODYN软件构建有限元模型，对接触爆炸荷载作用下带孔防护结构内冲击波的传播过程进行计算分析，以期获得冲击波各特征参数的变化规律。

1 计算参数标定

为了确认数值计算所采用软件的有效性、所使用材料状态方程及本构模型的适用性，对某直坑道化爆实验［9］进行相同工况下的数值分析，通过计算结果与实测结果的比对，判断误差是否在容许范围内。模型坑道总长21m，由21节可拼装式钢结构单元组成。每节坑道单元长1m，净截面尺寸为600mm×600mm。爆室壁厚80mm，其余单元壁厚20mm。坑道实物如图1所示，为保持模型坑道在爆炸实验过程中的气密性，坑道单元通过高强螺栓连接，单元与单元之间设有密封圈，坑道两端开口。

直通道测点布置如图2所示。直坑道口部化爆实验测点1（测点距炸药262cm）和测点4（测点距炸药562cm）的冲击波超压Δp计算波形和实测波形［9］的比较如图3所示。

图1 模型坑道实物Fig.1 Model tunnel for experiments

图2 直坑道化爆实验测点布置示意Fig.2 The layout of gaging points for chemical explosion experiments in the straight tunnel

图3 典型监测点的冲击波超压－时间曲线Fig.3 Overpressure－time curves of shock waves at typitcal gaging points

从直坑道化爆实验的情况来看，典型监测点计算波形与实验波形基本吻合，峰值压力、到达时间、正压作用时间误差均不超过20%。因此，本文计算所采用的方法和材料参数及求解设置均为合理的。

2 计算概述

为了便于计算，认为TNT炸药位于带孔防护结构中间孔口的正上方，利用对称性，建立1／2模型，如图4所示。带孔防护结构呈筒状，组成材料为混凝土和钢筋，混凝土标号为C35；结构高3 100mm，壁厚d＝500mm，内部半径r＝2 000mm；结构正面开有3个孔口，尺寸均为500mm×300mm，背面开有一尺寸为500mm×1 500mm的出入口；钢筋采用HRB335钢筋，直径为12mm，分3层布置，如图5所示。模型中，混凝土采用Lagrange网格进行建模，网格尺寸为50mm，底部固接；钢筋采用beam单元，采用join方法将钢筋与钢筋、钢筋与混凝土之间的节点进行连接；空气采用Euler算法，网格尺寸为25mm，边界面设为物质流出。由于带孔防护结构形状较复杂，在AUTODYN中建立模型难度大，对于该问题采用如下方法加以解决：在ANSYS主程序中建立带孔防护结构模型并对其进行网格划分，将所建模型保存成.k文件的形式，并导入AUTODYN。

图4 空气、炸药及结构模型的实体Fig.4 The model for air，TNT and structure

图5 钢筋模型Fig.5 The model for reinforcing bars

炸药爆炸的反应区压力模型采用高能炸药燃烧函数因子模型，状态方程采用JWL状态方程［10］，炸药爆轰过程中压力和比容的关系为：

式中：A、B、R1、R2、ω 为材料常数，p 为压力，V 为相对体积，e0为初始比内能。

计算中使用的TNT材料参数为：密度ρ＝1 630kg／m3，爆速D＝6 930m／s，A＝381.2GPa，B＝3.822GPa，R1＝4.15，R2＝0.9，ω＝0.35，e0＝6.12GJ／m3，pCJ＝21.4GPa。

空气采用理想气体状态方程［11］，初始密度为1.225kg／m3，比热比为1.4；混凝土采用修正后的RHT本构模型［12］；钢筋采用Johnson－Cook热黏塑性材料模型［13］。

TNT炸药位于结构正面中间孔口的正上方，到孔口上边缘的距离为x，为了研究爆炸位置对空气冲击波毁伤效应的影响，分别取x＝0.1，0.2，0.4，0.8m。TNT炸药质量为m，为了研究炸药当量对毁伤效应的影响，分别取m＝1，2，4，8kg。因此，模拟的工况共有4×4＝16个。

结构内布置2个观测点，测点1位于结构内部中心处，测点2位于孔口处，如图4所示。

3 带孔防护结构内冲击波传播规律分析

3.1 冲击波传播过程分析

分别对不同药量、不同爆炸位置的16个工况进行计算，得到了结构内部中心处（测点1）和孔口处（测点2）爆炸冲击波压力p随时间的变化过程。根据Δp＝pt－p0即可得到对应的爆炸冲击波超压。

图6 x＝0.4m，m＝4kg工况时，孔口处和结构内部中心处的冲击波超压－时间曲线Fig.6 Overpressure－time curves of the shock waves around the orifice and in the center of the structure in the case of x＝0.4m，m＝4kg

虽然在不同工况下炸药的质量和位置有所区别，但测点处爆炸冲击波超压－时间曲线的大致线形都较相似。以x＝0.4m，m＝4kg的工况为例，孔口处和结构内部中心处冲击波超压－时间曲线如图6所示。爆炸冲击波超压Δp最初为零，此时结构外TNT炸药虽然已经发生爆炸（炸药爆炸时为t＝0时刻），但由于冲击波尚未传播至测点处，故孔口处和结构内部中心处的空气压力仍与大气压力p0相同。当t＝0.616ms时，爆炸冲击波波阵面到达孔口处，测点2处压力急剧上升，并在t＝0.716ms时达到正压峰值，此时孔口处峰值超压Δpmax＝627.0kPa，如图6（a）所示。随后，t＝2.138ms时，冲击波到达结构内部中心处，此处超压由最初的零迅速达到正超压峰值Δpmax＝30.7kPa，如图6（b）所示。在爆炸冲击波波阵面经过结构内部中心处后，测点1的超压值并没有迅速地稳定为零，而是快速减小，直至出现了负压。随后，结构内部中心处出现了较长时间的超压震荡，并且在超压震荡期间，正压和负压的绝对值均较大，如图6（b）所示。

图7 炸药在带孔结构外表面爆炸的冲击波传播压力云图Fig.7 Pressure nephograms of the shock waves induced by TNT explosion outside the protective structure with holes

结合图7所示的冲击波传播压力云图，经分析可知，造成上述超压震荡现象的原因是：爆炸冲击波波阵面通过孔口进入结构内部（如图7（c）所示）后，迅速到达测点1处（如图7（d）所示），此时结构内部中心处出现峰值超压，随后，此处压力由正压逐渐变为负压；冲击波继续向前运动，并与结构内壁发生碰撞，产生的反射冲击波再次到达结构内部中心（如图7（e）所示），导致测点1处再次出现正超压；之后，冲击波与结构内壁发生了多次碰撞，反射冲击波多次经过结构内部中心，造成此处超压出现震荡，并且正压作用时间也大大增加。

当冲击波在结构内壁产生多次反射后，流场变得十分紊乱，呈现复杂波特性，部分位置压力将会增强，如图7（f）所示。

由图6可知，在相同的工况（x＝0.4m，m＝4kg）下，孔口处（测点2）爆炸冲击波波阵面的到达时间和峰值超压所对应的时刻显然要先于结构内部中心处（测点1），并且孔口处的峰值超压远大于结构内部中心处（627.0kPa≫30.7kPa）。在爆炸冲击波波阵面第1次经过孔口处后，测点2处的超压值虽然也会发生震荡，但震荡幅值相对于超压峰值要小得多。

在研究空气冲击波毁伤判据时，通常会考虑描述空气冲击波强弱的3个主要参数：峰值超压、正压区作用时间和正压冲量。冲击波峰值超压表示冲击波瞬间作用的量，而正压冲量表示在正压区时间范围内超压的持续作用量，在一定条件下两者均可作为目标毁伤程度的判定依据。在已获得冲击波超压－时间曲线的基础上，通过对超压Δp对正压作用时间t＋进行积分就可求得正压冲量I。由于图中冲击波超压在达到峰值后会出现一定时间的震荡，同时考虑到第1股爆炸冲击波对目标的毁伤最严重，所以在计算正压冲量时仅对第1段正压区进行积分。不同工况时带孔防护结构中心处冲击波的峰值超压和正压冲量如表1所示。

表1 不同工况下，结构内部中心处（测点1）冲击波的峰值超压和正压冲量Table 1 The maximum overpressures and the positive pressure impulse of the shock waves in the center of the structure（the gaging point 1）

在爆炸位置固定的前提下，比较不同TNT药量情况下的爆炸冲击波超压－时间曲线。如图8（a）所示，在工况1（x＝0.1m，m＝1kg）和工况2（x＝0.1m，m＝2kg）下，炸药的爆炸位置相同，均为x＝0.1m，工况2的药量是工况1的2倍，两者的冲击波超压－时间曲线线形相似，但工况2的峰值超压Δpmax，2＝74.803kPa高于工况1的峰值超压Δpmax，1＝37.335kPa，并且工况2的爆炸冲击波波阵面的到达时间和峰值超压所对应的时刻均要先于工况1。由此可知，爆炸位置固定的情况下，带孔防护结构内部中心处的峰值超压随着药量的增加而上升，爆炸冲击波波阵面的到达时间和峰值超压所对应的时刻随着药量的增加而提前。

在相同药量的前提下，研究TNT炸药爆炸位置对冲击波传播规律的影响。如图8（b）所示，在工况3（x＝0.1m，m＝4kg）和工况7（x＝0.2m，m＝4kg）下，炸药的药量相同，均为m＝4kg，工况7的炸药到正面中间孔口上边缘的距离x是工况3的2倍，两者的冲击波超压－时间曲线线形相似，但工况7的峰值超压Δpmax，7＝72.159kPa低于工况3的峰值超压 Δpmax，3＝143.524kPa，并且工况7的爆炸冲击波波阵面的到达时间和峰值超压所对应的时刻均要晚于工况2。由此可知，在相同药量的前提下，带孔防护结构内部中心处的峰值超压随着炸药至孔口距离的增加而下降，爆炸冲击波波阵面的到达时间和峰值超压所对应的时刻随着炸药至孔口距离的增大而推迟。

图8 不同工况下，测点1处的冲击波超压－时间曲线Fig.8 Overpressure－time curves of the shock waves at the gaging point 1in the different cases

3.2 冲击波特征参数计算公式

为了建立接触爆炸荷载作用下带孔防护结构内部中心处冲击波特征参数的计算公式，用量纲理论分析冲击波超压、正压冲量与时间的一般函数关系。

图9为弹药在混凝土结构表面爆炸的示意图，设装药爆炸能为Q，混凝土结构的主要几何参数分别为正面孔口面积S、孔口下边缘至结构内地面的距离h、背面入口面积Sm、内部半径r、净高H 和壁厚d。混凝土结构内部体积Vc＝πr2H，可由r和H 表示，不作为主定参量。爆心到混凝土结构孔口上边缘的距离为x，测点到混凝土结构墙壁内侧面的距离为l，周围大气压力和密度分别为p0和ρ0，测点处的峰值超压为Δpmax，从爆炸时刻算起，正压作用时间为t＋，正压冲量为I。

图9 弹药在混凝土结构表面爆炸示意图Fig.9 The layout of the explosive outside the concrete structure

于是该问题的主定参量组为：Q，S，Sm，r，h，H，d，x，l，p0，ρ0；待定参量组为：Δpmax，t＋，I。首先确定测点处的峰值超压Δpmax，如果用f表示爆炸冲击波参数关系，那么有以下函数：

采用L－M－T度量单位系统作为基本量纲，构建量纲矩阵：

并计算：

解齐次方程：

该齐次方程有12－3＝9个基本解，可取为：

峰值超压Δpmax的显式表达式为：

即：

式中：f1表示一个未定函数。

由上面的计算结果可以得知峰值超压Δpmax的主要影响因素为比例距，带孔防护结构的几何尺寸以及结构内径与爆心到孔口上边缘的距离之比

运用同样的方法，确定正压冲量：

3.3 数据拟合

在用实验方法整理函数（f·）的形式时，一般能在自变量的某一范围内采用指数关系式［14］，即π＝k0πk11·πk22·...·πknn。为了得到最终的拟合公式，将式（8）、（9）采用指数关系式表示，即：

对数值计算数据进行多元线性回归，得到式（10）、（11）中的待定参数αi，βi（i＝0，1，…，8），即可最终确定各待定函数（f·）的表达式，从而可以分别求得峰值超压Δpmax、正压冲量I的计算公式。

式（10）、（11）中待定变量过多，若在其基础上总结经验公式是比较困难的。对于本文所研究的防护结构，其几何尺寸、开孔大小都是确定的，此处可只保留结构内部半径r以反映结构尺寸对爆炸冲击波的影响。这样，函数拟合的难度就大大降低。结合数值计算所得结果（表1），针对本文所研究的带孔防护结构，结构内部中心处爆炸冲击波特征参数的拟合公式为：

与数值计算结果相比，由拟合公式计算所得的峰值超压和正压冲量的最大误差为25%（小于30%），满足冲击波特征参数预估的精度要求。

4 结 论

（1）当爆炸冲击波经孔口进入结构内部后，结构内部中心处的超压会快速到达正超压的峰值，并且该峰值大多出现在第1个波峰处。由于冲击波与结构内壁的反复碰撞，在峰值超压出现后，结构内部中心处的超压值会出现较长时间、较大幅值的震荡。

（2）孔口处爆炸冲击波波阵面的到达时间和峰值超压所对应的时刻先于结构内部中心处，并且孔口处峰值超压约为结构内部中心处的20倍。

（3）爆炸位置固定的情况下，带孔防护结构内部中心处的峰值超压随着药量的增加而上升，爆炸冲击波波阵面的到达时间和峰值超压所对应的时刻随着药量的增加而提前；在相同药量的前提下，带孔防护结构内部中心处的峰值超压随着炸药至孔口距离的增加而下降，爆炸冲击波波阵面的到达时刻和峰值超压所对应的时刻随着炸药至孔口距离的增大而后移。

（4）利用数值计算的结果，针对带孔防护结构，拟合得到了结构内部中心处爆炸冲击波特征参数的预估公式。与数值计算结果相比，由拟合公式计算所得的峰值超压和正压冲量的最大误差为25%。

［1］Henrych J.爆炸动力学及其应用［M］.熊建国，译.北京：科学出版社，1994：74－77.

［2］叶序双.爆炸作用理论基础［M］.南京：解放军理工大学，2001：64－66.

［3］Twisdale L A，Sues R H，Lavelle F M.Reliability－based analysis and design methods for reinforced concrete protective structures［M］.Tyndall Air Force Base：Air Force Civil Engineering Support Agency，1993：195－199.

［4］Welch C R.In－tunnel airblast engineering model for internal and external detonations［C］∥Proceedings of the 8th International Symposium on Interaction of the Effects of Munitions with Structures，1997：195－208.

［5］Anet B，Binggeli E D.The use of the LS 2000－design charts for predictions of air blast loading in tunnels due to HE－detonations at the tunnel entrance［C］∥Proceedings of the 9th International Symposium on Interaction of the Effects of Munitions with Structures.1999：225－234.

［6］庞伟宾，何翔，李茂生，等.空气冲击波在坑道内走时规律的实验研究［J］.爆炸与冲击，2003，23（6）：573－576.Pang Wei－bin，He Xiang，Li Mao－sheng，et al.The formula for airblast time of arrival in tunnel［J］.Explosion and Shock Waves，2003，23（6）：573－576.

［7］黄广炎，冯顺山，王云峰.装药约束对腔室爆炸冲击波传播特性影响的数值分析［J］.系统仿真学报，2009，21（1）：58－61.Huang Guang－yan，Feng Shun－shan，Wang Yun－feng.Numerical analysis for influence of restriction of charges on shock wave propagation inside cavity［J］.Journal of System Simulation，2009，21（1）：58－61.

［8］刘瑞朝，任新见.常规武器袭击下目标毁伤判据研究进展［J］.防护工程，2012，37（5）：60－66.Liu Rui－chao，Ren Xin－jian.Research progress of the damage criterion by conventional weapon［J］.Protective Engineering，2012，37（5）：60－66.

［9］任新见，张庆明，薛一江，等.坑道口部B炸药爆炸冲击波传播速度模型试验研究［J］.振动与冲击，2012，31（7）：71－73.Ren Xin－jian，Zhang Qing－ming，Xu Yi－jiang，et al.Scale model tests to determine in－tunnel blast impact wave velocity for B－charges at tunnel entrance［J］.Journal of Vibration and Shock，2012，31（7）：71－73.

［10］Welch C R.In－tunnel airblast engineering model for internal and external detonations［C］∥Proceedings of the 8th International Symposium on Interaction of the Effects of Munitions with Structures.1997：195－208.

［11］Century Dynamics Inc.Autodyn theory manual［M］.Concord，CA：Century Dynamics Inc，2006：133－141.

［12］Riedel W，Thoma K，Hiermaier S，et al.Penetration of reinforced concrete by BETA－B－500numerical analysis using a new macroscopic concrete model for hydrocodes［C］∥Proceedings of the 9th International Symposium on Interaction of the Effects of Munitions with Structures.1999：315－322.

［13］Johnson G R，Cook W H.A constitutive model and data for metals subjected to large stains，high stain rates and high temperatures［C］∥Proceedings of the 7th International Symposium on Ballistics，1983：541－547.

［14］谈庆明.量纲分析［M］.合肥：中国科学技术大学出版社，2005：120－123.