爆炸水池是开展武器弹药毁伤效能评估的基础条件，是完成战斗部破碎性实验，进行水下威力参数测定以及开展水下毁伤技术研究的重要场所^[1-2]。美国西南研究所对水槽、抑制挡板、填充沙隔墙、织物隔板的抗冲击侵彻能力进行研究，设计建造了内径9.15 m、深9.15 m的钢质圆顶、水槽式筒形钢筋混凝土爆轰室，可承受45 kg TNT当量的威力实验^[1]。目前国内同等尺寸爆炸水池的设计承试能力还远远达不到45 kg TNT当量。国内爆炸水池的工程设计通常是基于静态力学理论，按照等效静载的方式评估结构强度，但是对于高强度、高应变率的爆炸载荷，静态力学设计方法显然不能反映爆炸水池的实际动力学响应，也难以评估其结构的薄弱环节、危险区域以及作用机理。另外，破碎性实验^[2]是统计战斗部自然破片质量分布的重要手段，对杀伤力评估具有重要意义，为了隔绝水介质使壳体充分破碎，该实验通常将战斗部放置在空气桶里面引爆，为了削弱水中冲击载荷，池壁前方通常会设置一层气泡帷幕。空气桶和气泡帷幕将空气因素引入单一水介质，对水中冲击载荷的强度会产生不可忽视的显著影响，也大大增加了结构设计的难度。关于水中冲击波的传播规律和衰减规律已有大量研究^[3-4]，而对气泡帷幕防护作用的相关研究却比较少^[5-7]，关于空气桶的作用目前还停留在定性认识阶段，并没有相关论文发表。对于集成了空气桶、气泡帷幕、合金钢内衬、混凝土围堰、无限介质土壤和水介质的复杂抗爆体系，目前还没有成熟可靠的设计方法可以遵循，采用理论解析的方法更是难以实现。

因此本文中针对复合结构爆炸水池强度设计问题，利用非线性动力学程序LS-DYNA对10 kg TNT爆炸后的水中冲击波传播规律及爆炸水池结构动态响应情况进行数值模拟，对空气桶和气泡帷幕削弱水中冲击波的能力进行定量计算，指出内衬钢板和混凝土围堰的危险区域，分析冲击载荷作用以及结构动态响应机理，为相关爆炸水池的结构设计提供参考。

1.   计算模型与参数

1.1   模型的建立

爆炸实验水池如图 1所示，其基本参数为：内径10 m、深8 m；钢板内衬为15MnV合金钢材料，壁厚0.02 m；围堰为混凝土材料，厚1 m；设定土壤为无限介质，采用透射边界。装药为10 kg TNT球形炸药，爆心设置在距底面4 m的中心位置。空气柱直径0.5 m、高1 m。气泡帷幕厚度设为0.08 m。因此，本文中计算模型中包含6种材料：15MnV合金钢内衬、混凝土围堰、无限介质土壤、水、TNT炸药、以及空气。

图  1  爆炸实验水池示意图

Figure  1.  Sketch of explosion testing pool

下载: 全尺寸图片 幻灯片

为了研究不同实验条件下水池的冲击动力学性能，设计了4种不同的实验模型：模型Ⅰ为初始模型；模型Ⅱ在模型Ⅰ的基础上设置了空气桶；模型Ⅲ在模型Ⅰ的基础上设置气泡帷幕；模型Ⅳ同时设置了空气桶和气泡帷幕。设计模型Ⅰ的目的主要是为了参照P.Cole经验公式进行可行性分析，同时作为另外3种模型的对比模型。另外为了便于分析，取A、B、C等3个典型位置，如图 2所示：位置A为壁面中心，位置B为池壁与底面的交界处，位置C为底面中心。

图  2  爆炸水池模型及其网格划分情况

Figure  2.  Sketch of finite element structure and meshing condition

下载: 全尺寸图片 幻灯片

采用TrueGrid软件建立1/4三维数值计算模型, 如图 2所示，其中炸药区域采用在k文件中填充的方式“*INITIAL_VOLUME_FRACTION_GEOMETRY”加入，模型中水、炸药和空气采用ALE算法，编为一个ALE多物质组，无限介质土壤、混凝土围堰和钢板内衬采用Langrange算法，两者之间应用流-固耦合，空气上表面和土壤外侧施加透射边界。为了简化气泡帷幕模型，根据空气压缩机单位时间气流量以及气泡从池底升到水面的时间，将气泡帷幕等效为0.08 m厚的空气层，气泡帷幕也采用填充方式生成。

1.2   材料模型及参数

TNT炸药采用Mat_High_Explosive_Burn材料模型和JWL状态方程，具体参数如表 1所示^[8]，其中ρ为材料密度；D₀为炸药爆速；P_CJ为炸药爆炸产生的爆轰波C-J面压力；E₀为初始内能；A、B、R₁、R₂、ω为实验确定的常数。

表  1  TNT炸药材料参数

Table  1.  Material parameters of TNT explosive

ρ/(g·cm-3)	D0/(km·s-1)	PCJ/GPa	A/GPa	B/GPa	R1	R2	ω	E0/MJ
1.63	6.93	27	371	7.43	4.15	0.95	0.3	0.007

下载: 导出CSV 

| 显示表格

内衬钢板采用Plastic-Kinematic模型，该模型考虑了材料的弹塑性特性、强化效应以及应变率效应，材料参数见表 2^[9]，其中E为弹性模量；λ为泊松比；E_t为剪切弹性模量；σ_s为材料的弹性屈服强度；C、P为材料应变率相关常数；β为随动硬化系数。

表  2  内衬钢板材料参数

Table  2.  Material parameters of steel inner-lined plate

ρ/(g·cm-3)	E/GPa	λ	σs/GPa	Et/GPa	C	P	β
7 850	206	0.3	0.337	9.4	40.4	5	1

下载: 导出CSV 

| 显示表格

由于混凝土不直接承受冲击加载，不会产生破碎、裂纹等大变形，因此采用Plastic-Kinematic模型足以描述应力应变状态及风险区域，具体材料参数列于表 3^[10]。

表  3  混凝土材料参数

Table  3.  Material parameters of concrete

ρ/(g·cm-3)	E/GPa	λ	σs/GPa	Et/GPa	C	P	β
2 400	30	0.18	0.035	12.7	0	0	0.8

下载: 导出CSV 

| 显示表格

土壤采用Mat_Soil_And_Foam模型，密度ρ=1.8 g/cm³；泊松比λ=0.48；弹性模量E=47.38 MPa；体积载模量K=394.8 MPa，设土壤单元的初始体积为V₀，受力压缩后的体积为V₁，则体应变表示为：-ln(V₁/V₀)，土壤的压力与体应变关系如表 4^[10]所示。

表  4  土壤压力与体应变关系

Table  4.  Relation between pressure and volumetric strain

体应变	0.0	0.05	0.09	0.15	0.19	0.21	0.22	0.25
压力/MPa	0.0	0.34	0.45	1.27	2.08	2.71	3.92	5.66

下载: 导出CSV 

| 显示表格

空气采用空物质材料本构模型和线性多项式状态方程。水采用空物质材料本构模型和Grüneisen状态方程，具体参数列于表 5^[10]，其中c是u_s-u_p(冲击波速度-质点速度^[11])曲线的截距；S₁、S₂和S₃是u_s-u_p曲线斜率的系数；γ₀为Grüneisen常数；α是对γ₀的一阶体积修正；E₀为单位体积初始内能。

表  5  水的状态参数

Table  5.  State parameters of water

ρ/(g·cm-3)	c/(m·s-1)	S1	S2	S3	γ0	α	E0/MJ
1 020	1 650	1.92	-1.98	0.22	0.5	0	0

下载: 导出CSV 

| 显示表格

2.   水中爆炸冲击波压力模拟结果及分析

在大量理论和实验研究基础上，前人总结了一些水中爆炸冲击波参数的计算模型。应用最广的是P.Cole关于无限水中爆炸冲击波峰值压力计算公式^[12]：

式中:r为装药中心距测点的距离，m; p_m为测点r处得冲击波峰值压力，MPa; C为炸药当量, kg。公式(1)的适用范围为(20~200)r₀，r₀为装药半径。

将与装药中心等高、距离装药中心分别2.0、2.5、3.0、3.5、4.0、4.5和5.0 m的位置设为观测点。表 6是不同位置冲击波峰值压力数值计算结果与P.Cole公式计算值对比，除了2 m处数值模拟结果偏高10.6%，其他位置数值模拟结果与P. Cole公式相对值误差均小于±10%。表明水中爆炸冲击波作用的计算模型和参数可信，计算结果正确，满足计算精度要求。另外，池壁位置(r=5 m)的反射冲击波压力目前还没有经验公式可以计算，因为反射面的材质、结构特征、波阻抗等力学性质对水中冲击波的反射系数都有显著影响。水是可压缩性较小的流体介质，对于刚性壁面，当入射压力为1~100 MPa时，其反射系数K_r≈2^[13]。对于实际工程问题，则有1＜K_r＜2。

表  6  水中冲击波峰值压力计算值与理论值对比

Table  6.  Contrast of numerical calculation and empirical formulas

距离/m	峰值压力/MPa		计算与理论值的相对误差/%
	P. Cole理论计算值	数值计算结果平均值
2.0	57.97	64.1	10.6
2.5	45.05	47.0	4.3
3.0	36.66	37.4	2.0
3.5	30.80	29.8	-3.2
4.0	26.49	25.3	-4.5
4.5	23.19	21.6	-6.8
5.0	20.59	31.2*
注：由于池壁位置入射压力峰值被反射压力峰值覆盖，因此此值仅为反射压力峰值。

下载: 导出CSV 

| 显示表格

如图 3所示，对比拟合曲线与P.Cole压力曲线可以看出：初始阶段相对误差较大，但随着比例距离的增加，相对误差可以控制在10%以内。由于比例距离越小，冲击波高频成分越多，对传感器的响应频率要求越高，冲击波近场压力测试结果一般都会偏低，这也是P.Cole经验公式在比例距离较小时计算结果会偏低的原因。在比例距离较大时计算结果还是比较准确的。参照P.Cole经验公式，可以认为在本研究条件下，所建立的模型对于水中爆炸载荷的计算结果真实可信。

图  3  数值计算拟合曲线与P.Cole压力曲线的对比

Figure  3.  Contrast curves of numerical calculation and P. Cole empirical formulas

下载: 全尺寸图片 幻灯片

3.   防护性能及动力学响应分析

3.1   压力与比冲量载荷的变化规律

以模型Ⅰ为例，模拟结果如图 4所示：爆炸发生后，很快在水中形成冲击波，t=1.760 ms时，水中冲击波的波阵面清晰而狭窄，同时在爆心位置会形成由内向外传播的稀疏波。t=3.139 ms时冲击波到达左右边界，并发生正反射。冲击波压力在开口自由端迅速卸载，并形成向下传播的稀疏波，从而导致底面边界反射冲击波的作用时间减少、衰减速度较壁面反射冲击波明显提升，由图 4(c)可以看出：t=3.979 ms时，底面边界反射冲击波的强度明显低于壁面反射冲击波，波阵面也较为模糊。t=4.338 ms时，左右边界反射冲击波上端也同样受到了自由面传过来的稀疏波的影响，从底端到自由水面，反射冲击波的强度逐渐降低。

图  4  模型Ⅰ水中冲击波传播过程

Figure  4.  Propagation process of chock wave in water (model Ⅰ)

下载: 全尺寸图片 幻灯片

冲击波对边界钢板的作用效果不仅取决于冲击强度，还受到作用时间的显著影响。将A、C两个典型位置的冲击波强度曲线对时间积分，可以得到比冲量(单位面积冲量)变化曲线，如图 5所示，由于水中爆炸产生的气泡脉动以及半封闭空间冲击波反射叠加等效应，位置A比冲量从t=2.9 ms时刻起，经历了2次明显的阶跃上升，位置B比冲量从t=2.3 ms起，经历了2次阶跃上升和1段平缓上升。比冲量的第1次阶跃上升是首次冲击波反射压力作用的结果，位置A的作用时间大约为1.3 ms，要明显长于位置B的0.5 ms。首次冲击，位置A获得的比冲量为0.131 MPa·s，位置B获得的比冲量为0.054 MPa·s，位置A获得的比冲量高出位置B一倍以上。第2次阶跃上升开始于气泡二次脉动冲击波到达的时刻。t=15 ms时爆炸作用基本结束，位置A获得的冲量高出B位置30%。这种现象产生的主要原因是水池开口自由端形成的稀疏波主要向下传播，降低了底面反射冲击波的作用时间，并加快了其衰减速率，因此水池壁面中心位置最容易产生塑性变形，需要进行重点防护。

图  5  观察点A、C两点比冲量时程曲线

Figure  5.  Contrast curves of specific impulse in positions A, C

下载: 全尺寸图片 幻灯片

3.2   空气桶和气泡帷幕对冲击波的削弱作用

以装药中心同一水平高度、距离装药中心2.0 m位置为主要观测点，得到不同模型的冲击波压力对比曲线。图 6为模型Ⅱ与模型Ⅰ的对比曲线，可见设置空气桶对峰值压力有明显的削弱作用，由局部视图可以看出，无空气桶时，冲击波衰减曲线类似于直角斜三角形，陡然上升，平缓下降。有空气桶时，冲击波衰减曲线出现明显的波动，类似于锯齿形，这是因为爆炸产物首先驱动空气运动，形成空气冲击波，并迅速衰减，然后空气冲击波作用于水介质，形成水中冲击波，由于空气和水的波阻抗不同，会在两个介质之间形成透射波与反射波，如此反复几次，就会形成锯齿状波动的衰减曲线。随着距离的增加削弱效果略有增加，在距离爆心4.0 m处，本模型的削弱作用为48.8%。

图  6  模型Ⅱ与模型Ⅰ冲击波压力对比曲线

Figure  6.  Contrast curves of shock wave pressure (models Ⅰ & Ⅱ)

下载: 全尺寸图片 幻灯片

从图 7的对比曲线中可以看出，模型Ⅰ的冲击波压力时程曲线在5、10 ms左右出现了2次峰值，这是由于水中爆炸气泡脉动和冲击波反射叠加等效应引起的。但模型Ⅲ设置了气泡帷幕，从而抹平了该时刻的2个峰值，可见气泡帷幕对壁面反射冲击波的削弱作用十分显著，避免了反射波在中心汇聚产生二次脉动。由于气泡帷幕的作用，冲击波从水中进入空气，相当于从高阻抗介质进入一种很低阻抗的介质，会首先反射一个很强的卸载波，将冲击波压力迅速卸载掉。由图 7中气泡帷幕在冲击波作用下的压力云图可以看出：气泡帷幕类似于自由面，冲击波压力在该位置几乎被卸载殆尽，冲击波波阵面也逐渐模糊；最大压力出现在壁面中心位置，达到4.3 MPa，而模型Ⅰ无气泡帷幕情况下，反射压力高达31.2 MPa，气泡帷幕对冲击波反射压力的削弱能力高达86.2%。伍俊等^[14]在直径10 m、深10 m的爆炸水池中(混凝土围堰+钢衬内壁结构，气幕发生管道和空气压缩机与本文研究对象为同一规格产品)，测试了气泡帷幕的防护能力：1 kg TNT进行了两次对比实验，气泡帷幕对反射冲击波的削弱作用为89%；5 kg TNT仅进行了气幕防护实验，实验结果相对理论值削减了85.1%。本文计算的TNT当量为10 kg，得到气泡帷幕对反射冲击波的削弱作用为86.2%。随着入射冲击波强度的增加，气泡帷幕的削弱能力逐渐降低，对于5~20 MPa的入射压力，气泡帷幕的削弱作用在80%以上。

图  7  模型Ⅲ与模型Ⅰ冲击波压力对比曲线

Figure  7.  Contrast curves of shock wave pressure (models Ⅰ & Ⅲ)

下载: 全尺寸图片 幻灯片

3.3   钢板内衬动态响应分析

爆炸水池钢板内衬作为大型圆柱薄壳，其内部爆炸冲击载荷具有轴向对称性，因此剪切应力不会是主要屈服条件，可采用Mises屈服准则对其进行应力分析^[15]。本文中主要通过分析钢板内衬单元的Mises等效应力来评估塑性应变风险。数值计算结果显示，模型Ⅰ的钢板内衬没有产生任何塑性形变，位置B和C的应力都远未到达钢板屈服极限，只有位置A(壁面中心)有较大区域等效应力达到328.1 MPa，接近钢板屈服极限337 MPa，因此主要通过钢板内衬位置A的应力应变状态进行分析。

图 8给出了不同模型在位置A处所承受的比冲量时程曲线，空气桶的设置使位置A的比冲量削弱了16.2%，但是基本没有改变比冲量时程曲线阶跃上升的形态。气泡帷幕的设置使位置A的比冲量削弱了75.6%，同时消除了比冲量时程曲线陡峭的上升沿，使其上升更加平缓。图 9给出了不同模型在位置A处的等效应力变化曲线，各模型的应力峰值并非在冲击加载的瞬间就能到达，首次冲击从t=2.9 ms开始最多持续1 ms，但是应力峰值通常出现在t=6.5 ms，此时第1次冲击加载已经结束，钢板内衬在惯性作用下持续沿径向移动，达到最大位移，此时周向拉伸应力达到最大，也是在Mises准则下最容易发生屈服的时刻。

图  8  不同模型位置A比冲量时程曲线

Figure  8.  Contrast curves of specific impulse of steel in position A

下载: 全尺寸图片 幻灯片

图  9  不同模型位置A等效应力时程曲线

Figure  9.  Contrast curves of effective stress of steel in position A

下载: 全尺寸图片 幻灯片

表 7给出了空气桶和气泡帷幕对钢板内衬典型位置应力峰值的影响情况，气泡帷幕对位置A和B的衰减作用较为明显，空气桶对位置C的衰减作用较为明显。对于最容易产生塑性变形的位置A，空气桶和气泡帷幕的联合作用，使其峰值应力衰减了92.6%。气泡帷幕虽然不直接对位置C进行防护，但是大大减弱了壁面冲击波反射叠加从而形成二次脉动的能力，因此气泡帷幕对位置C的峰值应力有34.5%的削弱能力。

表  7  空气桶和气泡帷幕对钢板内衬典型位置应力峰值的影响

Table  7.  Effect of air tank and bubble curtain on effective stress in typical positions of steel lining

典型位置	应力峰值/MPa					衰减率/%
	模型Ⅰ	模型Ⅱ	模型Ⅲ	模型Ⅳ		模型Ⅱ	模型Ⅲ	模型Ⅳ
A	328.1	231.5	53.5	24.2		29.4	83.6	92.6
B	76.9	44.9	31.7	25.2		43.9	58.8	67.2
C	201.5	108.2	132.0	83.1		46.3	34.5	58.8

下载: 导出CSV 

| 显示表格

3.4   混凝土动态响应分析

混凝土作为一种波阻抗较高的介质，可以最大限度透射钢板内部应力，其自身阻抗越高，透射应力就越高。图 10给出了不同模型混凝土介质的塑性应变区域，模型Ⅰ底面中心位置C处产生的塑性变形最大，达到了0.24%，其次是底边折角B位置，大约在0.15%~0.17%的水平，最后是壁面中心点位置A，只有不到0.1%的塑性形变。在钢板内衬应力分析中，3个位置中位置A最危险，但是混凝土的应力情况却是位置A最安全。其原因主要有两点：一方面，钢板的惯性移动受到了混凝土的阻滞，而混凝土的惯性移动却不能被土壤有效阻滞，导致底面折角和底面中心处有应力集中，产生较大塑性形变，而池壁混凝土的惯性移动不会产生应力集中，所以塑性形变较小；另一方面，根据一维弹性应力波理论，土壤密度虽然和混凝土差距不大，但是弹性模量远远低于混凝土，从而导致土壤介质的声速和波阻抗都很小，透射应力波能力有限，导致混凝土背面产生明显的反射拉伸塑性应变区，加重了混凝土在位置A产生塑性变形的风险。在混凝土较容易发生塑性变形的区域，可以采用重晶石和铁矿石配制而成的重混凝土进行浇筑，以提高结构抗爆能力。

图  10  不同模型混凝土围堰的塑性应变区域

Figure  10.  Plastic strain zone of concrete cofferdam in each model

下载: 全尺寸图片 幻灯片

针对混凝土风险区域位置B和C，图 11~12给出了不同模型该位置的等效应力时程曲线。位置B作为应力集中最明显的区域，其应力在t=13 ms之前，维持了较高的平台，平台期大约为9 ms。随后进入弹性卸载阶段。由图 11可以看出，空气桶和气泡帷幕对位置B应力削弱明显，模型Ⅳ相对于模型Ⅰ削弱了56.7%。由图 12可以看出，位置C等效应力波动较为明显，与正反射冲击波压力作用和有关，该位置大约在t=9 ms时刻开始进入弹性卸载阶段，模型Ⅳ相对于模型Ⅰ应力峰值削弱了53.3%。

图  11  不同模型位置B的等效应力时程曲线

Figure  11.  Contrast curves of effective stress in position B

下载: 全尺寸图片 幻灯片

图  12  不同模型C位置的等效应力时程曲线

Figure  12.  Contrast curves of effective stress in position C

下载: 全尺寸图片 幻灯片

3.5   能量分析

在本计算模型中，爆炸能量主要以内能和动能的形式在各物质中体现，内能包括弹性变形能和塑性变形能，弹性变形能最终会通过其他形式释放出来，残存的能量为塑性变形能。图 13给出了不同模型钢板内衬的内能变化曲线，模型Ⅰ内能峰值为5.86 MJ，相当于1.3 kg TNT的能量，模型Ⅱ在此基础上削减了45.2%，模型Ⅲ在此基础上削减了96.3%，模型Ⅳ在此基础上削减了98.4%。气泡帷幕的作用要强于空气桶。内能变化曲线在经历第1个波动周期之后，逐渐趋于稳定，仅有小幅震荡，其中弹性变形能的成分要超过90%。图 14给出了不同模型混凝土的内能变化曲线，其波形与钢板内能变化曲线类似。都是在经历大的波动之后，逐渐趋于稳定，并有小幅震荡。模型Ⅰ混凝土获得的最大内能为0.6 MJ，相当于该模型中钢板内能峰值的10.2%，模型Ⅳ混凝土获得的最大内能为0.04 MJ，相当于该模型中钢板内能峰值的43.5%。可见混凝土吸收的能量比较小，其主要作用在于阻碍钢板的惯性移动，维持结构的稳固。

图  13  不同模型钢板内衬的内能变化曲线

Figure  13.  Contrast curves of internal energy in steel lining

下载: 全尺寸图片 幻灯片

图  14  不同模型混凝土的内能变化曲线

Figure  14.  Contrast curves of internal energy in concrete cofferdam

下载: 全尺寸图片 幻灯片

4.   结论

利用非线性动力学程序LS-DYNA对10 kg TNT爆炸后的水中冲击波传播规律及水介质威力池结构的动态响应情况进行了数值模拟研究，结果表明：

(1) 数值模拟能够较好地反映水中爆炸冲击波的分布规律和衰减规律，计算结果与经验公式的相对误差在10%以内，具有较高的可信度。

(2) 空气桶对冲击波峰值压力的削弱作用接近50%，对比冲量的削弱作用达到16.2%；气泡帷幕对壁面反射冲击波的削弱作用高达86.2%，对比冲量的削弱作用达到75.6%，气泡帷幕可有效减弱水中气泡脉动。

(3) 钢板的应力峰值出现在惯性移动最大时刻，弹性震荡周期为6.9ms，从钢板透射到混凝土中的能量仅占爆炸总能量的1.4%，混凝土吸能作用有限，其主要作用在于阻碍钢板的惯性移动，维持整体结构的稳固。混凝土围堰的危险区域和钢板截然相反，底面中心和底边折角处产生应力集中明显，且底面中心位置会有明显反射拉伸塑性应变区。