中空结构物在深水环境下工作时，表面将会承受高静水压力，当中空结构物突然被压溃，高压水流向压溃中心汇聚相撞，形成内爆（Implosion）^[1-2]，产生内爆冲击波，对周围物体造成破坏。2014年美国深海科研潜艇在深海工作时，中空结构的外部工作设备在静水压下被压溃引发内爆，对潜艇壳体造成损伤，使其“遇难”。在中微子探测领域^[3-4]，由玻璃制成的具有中空结构的光电倍增管（Photomultiplier tube，PMT）是中微子探测的核心部件，工作在深水环境时容易发生内爆。2001年，日本“超级神冈”中微子试验站发生PMT连锁爆炸事故，炸毁近8 000个PMT，直接经济损失达3 000万美元^[5]。除此之外，其他水下工作设备如管道、照明灯、摄像机等都有可能发生内爆，因此内爆机理研究具有重要意义。

国内外学者针对水下内爆问题开展了诸多研究。美国海军作战中心（NSWC）^[6]在装有水的压力罐中进行了圆柱壳的水下内爆实验，得到了铝制圆柱壳的内爆冲击波压力曲线，并结合仿真研究了圆柱壳的内爆冲击波压力曲线的主要特征；Diwan等^[7]在压力罐中对PMT进行了2次水下内爆试验，得到了PMT内爆发生过程及内爆冲击波压力曲线，并利用任意拉格朗日欧拉（ALE）方法对PMT内爆进行了仿真，得到的仿真结果与试验结果吻合较好；杜志鹏等^[8]将水下爆炸气泡动力学与水下内爆相结合，基于能量守恒关系，推导出不可压缩流体中球形容器内爆理论模型；黄治新等^[9]根据应力波原理，提出了利用大型压力罐进行内爆试验的方法，并通过此方法测得了PMT内爆冲击波。

上述研究主要针对特定结构的圆柱壳、PMT进行试验和数值模拟研究，没有考虑静水压、真空体积对中空结构物内爆的影响。本研究拟开展PMT水下内爆试验，以验证ABAQUS软件CEL耦合计算方法模拟水下内爆的准确性，并通过有限元模拟方式探究外界静水压力、真空体积对内爆冲击波的影响规律，为深入分析中空结构物的水下内爆提供参考。

1.   PMT内爆试验与数值模拟

1.1   PMT水下内爆试验

PMT实物见图1，其形状如灯泡，总高度0.700 m，头部球体直径0.508 m，尾部圆柱直径0.100 m；其外壳材料为玻璃，厚度为5 mm，内部真空度为10⁻⁴ Pa。PMT水下内爆试验装置为直径3 m、长度5 m、壁厚30 mm的圆柱形钢制压力罐，其示意图见图2。压力罐侧面设有透明视窗，用于照明并放置高速摄影设备，可清晰地记录内爆全过程；压力罐内部设有工作平台，工作平台上放置液压式挤压装置，用于安装和引爆PMT。受试验装置尺寸限制，将压力罐装入水时，无法模拟PMT在50 m水深的工作环境，为此利用外置空压机对压力罐内施加0.5 MPa压力进行模拟。通常情况下，采用水中压力值代表冲击波压力，故在试验过程中主要测量水中压力值。试验前，压力罐内共布置4个PCB动水压传感器，测点F₁位于PMT正上方，距PMT中心0.41 m；测点F₂位于PMT赤道平面靠近压力罐入口处，距PMT中心0.55 m；测点F₃、测点F₄位于PMT赤道平面远离压力罐入口处，分别距PMT中心0.55和1.35 m。试验时，液压装置挤压PMT发生内爆，试验成功测得PMT水下内爆过程和压力数据。

图  1  PMT实物

Figure  1.  Picture of actual PMT

下载: 全尺寸图片 幻灯片

图  2  压力罐示意图

Figure  2.  Schematic of the pressure tank

下载: 全尺寸图片 幻灯片

1.2   水下内爆试验结果分析

完成PMT内爆试验后，根据高速摄影所记录的影像，对内爆过程进行详细分析。将PMT水下内爆发生时间点记为零时刻（t = 0 ms），图3为PMT发生内爆典型特征时刻的图像。可以看出，PMT内爆过程可分为4个典型阶段：（1）挤压阶段（−16.00～−14.00 ms），挤压装置挤压PMT使玻璃外壳产生微小变形，挤压装置附近玻璃外壳亮点缩小；（2）裂纹传播阶段（−14.00～−8.00 ms），PMT玻璃外壳在−12.33 ms产生初始裂纹，−8.00 ms裂纹迅速传播至整个球面；（3）整体压溃阶段（−8.00～0 ms），PMT玻璃外壳在0 ms完全被压溃，碎片迅速向中心运动，各方向水流发生碰撞，产生内爆冲击波；（4）冲击波传播阶段（0 ms以后），玻璃碎片迅速向外飞散，冲击波以近似声速向外传播，图像出现黑色斑点，这是因为冲击波到达罐壁时发生反射，反射波与入射波发生叠加产生了空化气泡。上述阶段划分充分展现了PMT水下内爆发生过程及内爆冲击波的产生过程，对水下内爆数值模拟具有一定的指导意义。

图  3  PMT内爆过程

Figure  3.  Implosion process of PMT

下载: 全尺寸图片 幻灯片

1.3   PMT水下内爆数值模拟

基于PMT水下内爆的试验结果，对PMT水下内爆进行数值模拟。由试验结果可知，在PMT内爆过程中，从裂纹产生到内爆发生共经历了14 ms，具有瞬时性；此外，考虑到PMT所处的环境，在内爆发生过程中还会存在流体与玻璃之间的耦合作用。根据上述分析，可利用有限元软件ABAQUS中的CEL算法对PMT水下内爆进行数值模拟。

图4为PMT内爆有限元模型中部截面图，整个球体水域采用欧拉单元（EC3D8R），中间网格密度大，外部网格密度小，以减小计算时间。红色部分的欧拉单元赋予水材料属性，水采用Grüneisen状态方程描述，密度为1.0 g/cm³，波速为1 480 m/s，黏度为1×10^–3 Pa·s，在水域中利用初始应力场施加0.5 MPa的初始静压力，且在水域外边界设置冲击波无反射、流体自由流入和流出的边界条件，以模拟无限水域。蓝色部分的欧拉网格赋予真空材料属性。红色部分与蓝色部分交界处为PMT玻璃外壳，采用六面体拉格朗日单元（C3D8R），材料参数如下：密度为2 230 kg/m³，弹性模量为2.77 GPa，泊松比为0.376，断裂极限为40 MPa，玻璃失效位移为0.1 µm。在模型中输出与试验相同测点的压力值，用于后期数据的相互验证及对比分析。数值模拟开始时，对PMT外表面局部区域施加瞬时压力，以模拟液压装置对PMT的挤压，使得PMT破碎，从而发生水下内爆。

图  4  PMT内爆仿真模型

Figure  4.  Implosion simulation model of PMT

下载: 全尺寸图片 幻灯片

欧拉网格具有与材料完全隔离的特性，根据网格内材料的占比不同，会呈现不同的网格颜色。当网格呈红色时，欧拉网格中水的占比为100%，即水填满整个网格；当网格呈蓝色时，水的占比为零，即没有水填充到网格。图5为数值模拟得到的流场运动过程。当PMT受到挤压时，局部发生破碎，水流从挤压位置涌入PMT内，导致PMT的玻璃外壳全部压溃，各方向水流迅速向PMT中心处汇聚并发生碰撞，产生冲击波。整个PMT内爆过程共经历13.6 ms，与高速摄像拍摄的内爆过程及发生内爆的时间相近，说明采用CEL耦合计算方法能够准确地模拟PMT发生内爆时水域流场变化过程。

图  5  水域流场变化过程

Figure  5.  Change process of flow field in the water

下载: 全尺寸图片 幻灯片

1.4   试验与数值模拟结果对比分析

为进一步验证PMT内爆数值模拟方法的正确性，对数值模拟和试验的压力测点数据进行对比分析。由于试验中使用的是动压传感器，无法测出水中的初始静压力，所以将试验中的压力值均增加0.5 MPa，以模拟所受的静水压作用。进一步调整试验压力值的起始位置，使试验与数值模拟压力峰值相对应。图6为各测点的内爆冲击波压力对比曲线，可以看出：无论是试验还是数值模拟，内爆冲击波压力峰值均较高，脉宽较小。数值模拟的压力曲线在零时刻附近出现振荡，这是由数值模拟过程中测点处的能量平衡引起的。在未达到压力峰值的曲线段，数值模拟与试验中均出现微弱压力降，这是由PMT玻璃外壳完全破碎后水流内涌存在速度势导致的。达到压力峰值后，数值模拟曲线快速恢复到平衡位置，而试验曲线出现振荡式衰减，这是由于数值模拟中设置了无反射边界的水域条件，而试验中由于压力罐尺寸的限制，罐壁对冲击波具有一定的反射作用。总体而言，数值模拟得到的内爆冲击波脉宽、峰值与试验结果吻合较好。

图  6  内爆压力对比曲线

Figure  6.  Comparison curves of implosion pressure

下载: 全尺寸图片 幻灯片

进一步对试验与数值模拟得到的压力峰值进行对比，结果列于表1。各测点峰值中，模拟峰值普遍小于试验峰值，这是由于数值模拟中玻璃材料本构方程与实际存在差异，致使玻璃的破碎速率不同，从而对内爆冲击波压力峰值产生一定的影响；此外，模拟时模型网格大小对结果精度也有一定的影响，使得冲击波压力峰值降低，但整体平均误差为20.8%，满足实际工程需求。选取内爆试验压力测点F₁的压力数据进行积分，得到比冲量峰值为2 700 Pa·s，该能量相当于50 g标准TNT炸药按Geers & Hunter水下爆炸理论公式计算所得的能量。

表  1  试验和模拟得到的冲击波压力峰值对比

Table  1.  Comparison of tested and simulated shock wave pressure peaks

Measuring point	Pressure peak/MPa		Error/%
	Test	Simulation
F1	14.63	11.28	22.90
F2	7.68	5.53	27.90
F3	6.18	5.53	10.51
F4	3.19	2.49	21.90

下载: 导出CSV 

| 显示表格

2.   内爆冲击波影响因素分析

在PMT水下内爆数值模拟方法正确建立的基础上，利用数值模拟对内爆冲击波的影响因素进行研究。为了方便建模，采用体积等效的方式将PMT等效为直径0.5 m的真空球体。

2.1   静水压力的影响

水下内爆发生时PMT玻璃外壳瞬间破碎，水流急速碰撞产生冲击波。为了简化计算，仅建立水域及真空球体区域模型。图7为建立的有限元简化模型。外部水域直径为5 m，内部真空球体直径为0.5 m；共设置5个压力测点，编号为S₁～S₅，分别距球心0.25、0.35、0.45、0.55、0.65 m；简化模型的材料本构、边界条件、载荷等计算参数保持不变；根据PMT的实际工作环境，共设置了静水压分别为0.3、0.4、0.5、0.6及0.7 MPa的5种内爆工况进行数值模拟。

图  7  内爆简化模型

Figure  7.  Simplified model of implosion

下载: 全尺寸图片 幻灯片

图8显示了0.5 MPa静水压作用下水域流场的变化过程。可以看出：计算开始时，水流迅速向中心流动，初始阶段，水流速度由零开始进行加速，水流总体运动较慢；随着时间的延长，水流获得了较大的运动速度，向中心加快运动，经历10.25 ms后，水流发生碰撞产生内爆冲击波。

图  8  内爆水流涌入过程

Figure  8.  Water influx of implosion

下载: 全尺寸图片 幻灯片

进一步对测点的压力情况进行分析，图9为真空球体在5种工况下，水流高速碰撞产生冲击波，各测点的冲击波压力时程曲线。可以看出：各工况下压力曲线的变化趋势大致相同，首先在短时间内达到内爆冲击波峰值，随后压力值出现了不同程度的振荡，最后均稳定为零压力，且各测点的振荡形式相似；同一静水压下，随着测点距离的增大，内爆冲击波压力峰值显著减小。比较各工况下同一测点测得的冲击波压力峰值对应的时间可知，随着静水压力的增加，内爆发生的时刻提前。

图  9  压力时程曲线

Figure  9.  Curves of pressure varied with time

下载: 全尺寸图片 幻灯片

对图9中的冲击波压力峰值进行分析，以测点距离为横坐标，各测点冲击波压力峰值为纵坐标，绘制出各静水压下冲击波压力峰值随测点距离的变化规律，如图10（a）所示。分析发现，同一静水压下，冲击波压力峰值p_m随测点距离x的变化规律可用指数函数p_m = ax^b进行拟合，拟合结果如表2所示。冲击波峰值衰减指数b在−1.067～−1.107范围，同比炸药的水下爆炸冲击波传播衰减指数（1.18）偏小。以冲击波压力峰值的10%对应的两个时刻之间的差值定义为内爆冲击波脉宽，选取0.3 MPa静水压下各测点脉宽进行讨论，脉宽随测点距离变化曲线如图10（b）所示，脉宽随距离的增大缓慢增大。

表  2  拟合系数值

Table  2.  Fitting coefficient values

Hydrostatic pressure/MPa	a	b
0.3	4.216	−1.067
0.4	4.957	−1.096
0.5	5.761	−1.084
0.6	6.488	−1.079
0.7	7.001	−1.107

下载: 导出CSV 

| 显示表格

图  10  冲击波压力峰值和脉宽随测点距离变化曲线

Figure  10.  Curves of shock wave peak and pulse width with measuring point distance

下载: 全尺寸图片 幻灯片

为准确分析静水压对内爆冲击波的影响规律，以静水压为横坐标，测点的冲击波压力峰值为纵坐标，对图9中各工况下不同测点冲击波压力峰值重新进行统计，如图11（a）所示。可以看出，对于同一测点，冲击波峰值压力p_m随外界静水压力p_h呈线性增长趋势，可用函数p_m = cp_h + d对其进行线性拟合，所得拟合系数值如表3所示，斜率c随测点距离的增大而减小，即距离内爆中心越远，冲击波压力峰值变化越平缓。选取各静水压下0.25 m处测点脉宽进行讨论，脉宽随静水压变化曲线如图11（b）所示，脉宽随静水压力的增加基本保持不变。

表  3  拟合系数值

Table  3.  Fitting coefficient values

Measuring point distance/m	c	d
0.25	34.075	8.631
0.35	23.771	6.004
0.45	17.434	4.770
0.55	14.026	3.932
0.65	11.409	3.338

下载: 导出CSV 

| 显示表格

图  11  冲击波压力峰值和脉宽随静水压变化曲线

Figure  11.  Curves of shock wave peak and pulse width with hydrostatic pressure

下载: 全尺寸图片 幻灯片

2.2   真空体积的影响

将静水压设为0.5 MPa，选取与2.1节相同的压力测点，真空球体半径依次修改为0.15、0.20、0.25、0.30及0.35 m，进行5种工况下的内爆数值模拟，研究真空体积对内爆冲击波压力峰值的影响规律。图12为不同工况各测点的内爆冲击波压力曲线。从图12可以看出：真空球体半径越大，内爆发生的时刻延后，冲击波峰值越大。

图  12  不同真空球体半径的冲击波压力时域曲线

Figure  12.  Curves of shock wave pressure with time under different vacuum sphere radii

下载: 全尺寸图片 幻灯片

对图12中的冲击波压力峰值数据进行统计，以真空球体半径为横坐标，冲击波压力峰值为纵坐标，绘制图像如图13(a)所示。

图  13  冲击波压力峰值和脉宽随真空半径变化曲线

Figure  13.  Curves of shock wave peak and pulse width with vacuum radius

下载: 全尺寸图片 幻灯片

从图13(a)可以看出，冲击波压力峰值p_m随真空球体半径r呈线性变化。利用p_m = mr + n进行线性拟合，拟合系数值列于表4。分析可得，当距内爆中心位置一定时，冲击波压力峰值随真空体积的增大呈线性增大，且距离内爆中心越近的位置增长越快。选取各工况下0.25 m处测点进行脉宽讨论，脉宽随真空半径变化曲线如图13(b)所示，脉宽随真空半径的增大缓慢降低。

表  4  拟合系数值

Table  4.  Fitting coefficient values

Measuring point distance/m	m	n
0.25	147.97	−11.125
0.35	100.68	−7.623
0.45	74.26	−5.351
0.55	58.72	−3.998
0.65	47.55	−3.045

下载: 导出CSV 

| 显示表格

3.   结　论

（1）对比PMT内爆试验与数值模拟结果可知，ABAQUS的CEL数值模拟方法在模拟水下内爆发生过程、水域流场变化过程、冲击波压力峰值和脉宽等方面具有较高的模拟精度，能够准确、全面地模拟PMT内爆的物理过程。

（2）通过对不同静水压力下的内爆进行数值模拟，得到了静水压力对内爆的影响规律：内爆发生的时刻随静水压力的增加而提前；内爆产生的冲击波压力峰值随测点距离增大呈指数衰减，且衰减系数比水下炸药爆炸冲击波衰减系数小；同一位置的内爆冲击波压力峰值随静水压力增加呈线性增长趋势，且距离内爆中心越远，冲击波压力峰值增长越缓慢；冲击波脉宽随静水压力的增加基本保持不变。

（3）通过对不同真空体积的内爆进行模拟，得到了真空体积对内爆的影响规律：内爆发生时刻随着真空体积的增大而延后；同一位置的内爆冲击波峰值随真空球体体积的增大呈线性增长，且距离内爆中心越近，冲击波峰值变化越快；冲击波脉宽随真空半径的增大缓慢降低。