PMT水下内爆冲击波强度与传播特性研究

金 键， 侯海量， 陈鹏宇， 朱 锡， 吴林杰， 何 苗, 王天穹

(1.海军工程大学 舰船工程系，武汉 430033；2.中国科学院 高能物理研究所，北京 100049)

中微子是组成自然界的基本粒子，它不仅是物理学家们探索宇宙形成和演化的途径，同时，其不带电、微质量、光速传递等特点使其有望作为理想信息传递载体在通讯、军事等领域带来革命性的突破。目前，世界上共有27个主要的中微子探测试验，我国在大亚湾反应堆中微子实验中，发现了中微子的第三种振荡模式，首次精确地测定了中微子的混合角，被誉为基础物理学的一大突破[1]。大亚湾实验二期工程即江门中微子实验(在建)旨在解决中微子质量谱排序问题[2]。

光电倍增管(Photo-Mutiplier Tube,PMT)是用于探测宇宙中微子的核心部件，2万只PMT依附在水中的中心探测器上构成了光信号的探测构件。2001年日本超级神冈中微子探测器发生PMT内爆事故[3]，其中心探测器底部一个PMT发生内爆后，冲击波导致7 877个PMT发生连环爆炸，爆炸直接损失超过3 000万美元，延误实验时间4年。这次事故使各国大型水下中微子探测工程开始重视水下PMT的防爆必要性。

PMT是由透明玻璃吹制而成，内部接近真空，工作状态下承受静水压力。PMT的意外破坏有多种原因：①玻璃表面在静水压下发生蠕变；②工作周期内水池反复的排水、注水致使PMT承压状态变化，从而导致玻壳发生疲劳；③水中微生物对玻壳的侵蚀；④不可避免的初始缺陷以及安装过程中的磨损、擦伤。当PMT在静水压力下发生内爆后，产生的冲击波在向四周传播的过程中有可能击破相邻的PMT，从而导致水中的所有PMT发生连锁内爆。Ling等[4]利用有限元仿真的方法模拟了多个PMT水下殉爆过程，说明了初始内爆引发水下PMT连环内爆的可能性。

本文针对江门中微子实验中的PMT防爆工作，对PMT水下内爆冲击波的产生、气泡压溃时间、冲击波的强度以及冲击波的传播规律开展了理论分析和数值仿真研究。

1 工程背景

江门中微子中心探测器安装于直径43 m、深45 m的圆柱型水池中(如图1所示)；中心探测器外层与水池之间注满超纯水，密度为1 000 kg/m3；20 000只PMT密布在图1所示的不锈钢网架和水池内壁上，水池内的PMT承受0.1～0.54 MPa的静水压力。

图1 中心探测器结构

PMT实物图及外形尺寸如图2图所示。PMT由玻璃外壳(简称玻壳)及内部电气原件组成。玻壳厚度约为4～6 mm。内部电气结构为金属材料，工作在真空之中，内部真空度为10-5Pa以上。

(a) PMT示意图

(b) PMT外形尺寸图(mm)

图2 PMT示意图及外形尺寸图

Fig.2 Diagram and dimension of PMT

2 球形容器水下内爆理论分析

球形容器水下内爆过程类似于水下爆炸后的气泡脉动过程和水下空泡破裂过程。外部水压高于气泡内部气体压力，一旦平衡状态被打破，外部高压水迅速向中心汇涌，压缩气泡内部气体，直到气泡达到最小半径之后，气泡开始膨胀并向外辐射冲击波。有所区别的是气泡形成的方式，水下爆炸产生的气泡是由爆轰产物的混合物形成的；典型水下空泡由相变、成核与空化产生[5]；而本例中则由容器整体破坏产生。

在对气泡动力学特性进行分析时，通常以无限均匀温度液体区域内的单个球形对称气泡为研究对象，而热效应、浮力、液体的黏性、液体的可压缩性和空泡表面张力等则是影响因素，在具体研究中为了简化问题往往忽略一些因素。Rayleigh[6]首先推导并应用了不考虑表面张力和黏性的空泡动力学方程；Plesset[7]在Rayleigh工作的基础上得到了广义的空泡动力学R-P方程(式(1))。

(1)

式中:R为球形气泡半径；Pi是气泡内初始压力;p∞为距离气泡无穷远处的静压力;μ,ρ分别为液体的黏性和密度；σ为表面张力。

(2)

当气体含量很小时，可近似认为气泡半径被压缩为0，从而得到了气泡由R=R0到R=0完全压缩所需的时间为

(3)

在对气泡回弹阶段冲击波形成的研究中，Ivany等[8]通过研究证实表面张力和黏性对冲击波没有重要影响。Hicking等[9]利用可压缩流动方程的数值解探究了气泡回弹阶段冲击波的形成过程，并指出液体的可压缩性对气泡动力学特性的影响是轻微的，其主要作用是在气泡破裂后回弹阶段中对冲击波的形成，得到了压力脉冲远离气泡传播时呈现出近似几何形状衰减的结论。

研究球形容器水下内爆可以利用气泡动力学的相关理论，但也有其特殊性。从能量的角度分析，由于容器薄壁的存在，容器空腔蓄积的势能一部分需要用来打破容器壁并加速碎片向内运动，这必然会使容器水下内爆产生的冲击波峰值和比冲量比相同尺寸、相同气体压力的气泡破裂所产生的冲击波峰值和比冲量要小。另外，Turner[10]通过将玻璃球水下内爆实验与数值仿真结果进行对比，得出了容器结构的破坏时间对压力峰值具有很大影响的结论。

3 数值模拟

由于水池深45 m，对整个水池进行全尺寸建模将会带来巨大的网格数及计算量，因此只选取45 m水深处PMT的工作环境进行建模(如图3所示)；又因为模型是轴对称的，为进一步减少计算量，仅建立四分之一模型。有限元计算模型使用Patran进行几何模型建模，采用kg-mm-ms单位制。

图3 有限元计算模型示意图

建模时，以PMT上半球最大半径处几何中心作为原点，在X轴方向上，0～8 000 mm划分100个发散型网格；Y轴方向上，-8 000～0 mm、0～8 000 mm分别划分100个发散型网格；Z轴方向上的网格划分与X轴相同；水域外围边界设置为固支边界，中间设置为对称边界。

用稀薄空气来代替PMT内部真空度，假设空气为理想气体，其状态方程如下

P=(γ-1)ρe

(4)

式中:P为空气压力，设为10-5Pa；ρ为空气密度；e为内能，设为1.972 7×105J/kg；γ为比热，设为1.4。

由式(4)可知，要使PMT内部真空度达到10-5Pa，对应的内部稀薄空气密度应为1.28×10-10kg/m3。

水的状态方程采用Mie-Gruneisen方程

(5)

(6)

式中:p为水的压力；μ为水的压缩度；V0为水的相对初始体积；ρ0为水的密度，取1 000 kg/m3；c为vs-vp曲线的斜率(声速)，取1 484 m/s；γ0为格吕内森参数，取0.11；a为γ0和μ的一阶体积修正量，取3；S1、S2、S3分别为曲线拟合参数，取S1=1.979、S2=S3=0。E为水的初始体积内能，由于水几乎不可压缩，取值设为0。

基于安全性考虑，选取水池最深45 m处为初始水压，即0.54 MPa。通过对式(5)、(6)的计算可以得到v0=0.999 754 6。

PMT玻璃外壳采用壳单元建立(如图4所示)，利用这些壳单元组成的Part对内部空间进行初始化，使用LS-DYNA中的*INITIAL_VOLUME_FRAC-TION_GEOMETRY关键字将壳单元内部包裹空间初始化成稀薄空气。

图4 PMT玻壳有限元模型

为了在计算开始时删除图4所示的壳单元以便水涌入空穴中，在LS-DYNA的材料模型中添加*MAT_ADD_EROSION关键字，并在FAILTM中加入失效时间。该卡片可以让不易失效的材料在计算开始后的某个时刻强制失效，本例中将失效时间设为0.001 ms。

有限元模型中的流体全部为欧拉单元，利用LS-DYNA程序的ALE算法对模型进行计算。

4 算法验证

为了验证该算法的准确性，使用该算法对Turner做的玻璃球水下内爆实进行数值模拟计算。实验工况如图5所示，玻璃球直径762 mm、内部压力为0.101 3 MPa的玻璃球在6.996 MPa的静水压力下内爆。三个夹角120°布置的自由场压力传感器距球心都为101.6 mm。

图5 玻璃球水下内爆实验布置图

Fig.5 Experimental arrangement of glass implosion on underwater

建立有限元模型时，水的状态方程中的初始体积V0设为0.996 843，空气的密度ρ设为1.29×10-9kg/m3，球形空气半径设为37.3 mm(玻璃球厚度0.8 mm)。数值仿真的计算结果(距离球心101.6 mm处某单元压力、超压比冲量时程曲线)如图6、7所示。

图6 压力时程曲线

图7 超压比冲量时程曲线

计算开始后，从气泡表面径向向外辐射稀疏波，稀疏波经过的区域，其初始静压会突然下降，如图6所示，大约0.04 ms时刻压力大幅下降了3.039 MPa。随着高压水进一步压缩气泡，压力值又缓慢上升，直到气泡回弹阶段，出现了较大的压力峰值。数值仿真结果和实验所测均值列于表1。

表1实验与数值仿真结果对比

Tab.1Comparisionbetweenexperimentresultandnumericalsimulationresult

压力峰值/MPa超压比冲量/(Pa·s)仿真值24.47672.14实验值26.10601.33误差-6.24%11.75%

对比表1所示结果，就超压比冲量而言，仿真结果要略大于实验值，而压力峰值的仿真结果却小于实验值。这是因为波阵面的厚度非常小，约为10-5～10-6cm[11]，而有限元软件在输出压力结果时是以一个单元上的平均压力输出，单元长度又远远大于波阵面厚度，从而导致压力峰值的仿真值小于实际值。另外，由于该实验中无法准确获得容器压溃时间，仅对气泡压溃时间的仿真值与理论值进行对比。气泡压溃时间的仿真值为0.406 ms，根据方程4得到的半径37.3 mm气泡压溃时间为0.411 ms，两者误差1.2%。总体而言，数值仿真计算结果与实验值的误差控制在一定范围之内，两者吻合较好。

5 PMT水下内爆数值仿真结果及分析

对PMT在0.54 MPa静水压力下内爆过程进行数值仿真计算时，以单元中心到原点的距离作为衡量值，选取X+,Y+,Y-三个方向上不同距离处的6个单元作为压力输出点，各点压力时程曲线及比冲量时程曲线如图8～图13所示。

图8 X+方向上各点压力时程曲线

Fig.8 Pressure time history for orientation ofX+

图9 X+方向上各点比冲量时程曲线

图10 Y+方向上各点压力时程曲线

Fig.10 Pressure time history for orientation ofY+

从压力时程曲线来看，Y+方向上155 mm处单元峰值压力发生在9.222 ms，Y-方向上136 mm处单元峰值压力发生在9.237 ms，假设冲击波在水中的以声速1 484 m/s向外均匀传播，根据两者压力峰值到达时间差得到PMT内爆中心应在PMT几何中心上方20.56 mm处。这时，X+方向上141 mm处单元距内爆中心约142.5 mm，该单元压力峰值发生在9.227 5 ms，相比Y方向上冲击波的传播速度，X方向上冲击波的传播速度比其仅大0.4%，因此，可以认为PMT内爆中心坐标为(0，20.56，0)。

根据式(4)可以计算出半径250 mm球状气泡溃灭的理论时长为9.844 ms；在数值仿真计算中，距内爆中心142.5 mm处压力峰值发生在9.227 5 ms，以声速1 484 m/s作为冲击波传播速度，则内爆中心的向外传播冲击波的时间为9.13 ms；数值仿真得到的气泡溃灭时长比理论时长大7.8%，两者基本吻合。

图11 Y+方向上各点比冲量时程曲线

Fig.11 Impulse time history for orientation ofY+

图12 Y-方向上各点压力时程曲线

Fig.12 Pressure time history for orientation ofY-

图13 Y-方向上各点比冲量时程曲线

Fig.13 Impulse time history for orientation ofY-

由于PMT外形轮廓不是球形，导致其内爆产生的冲击波并不是以球形波阵面向外传播，冲击波在X+、Y+、Y-三个方向上的传播存在一定的差异。通过对三个方向上距内爆中心不同距离处的压力峰值、超压比冲量进行曲线拟合(如图14、15所示)，可以较为直观对其进行比较分析。

从图14、15可以看出：①对于压力峰值拟合曲线，Y方向上的两个幂函数拟合曲线与测点吻合较好，且两个拟合曲线的整体趋势基本相同；而X+方向上的幂函数拟合曲线与测量点的偏差较大，这是因为X+方向所取的测点单元并不与PMT内爆中心在一条直线上，且冲击波向外传播的能量在空间上是不均匀的；②对于超压比冲量拟合曲线，各数据点与拟合曲线吻合较好；③距内爆中心相同距离处，Y-方向上压力峰值、超压比冲量明显比另外两个方向上的压力峰值要大；④PMT内爆后产生的冲击波在远离气泡传播时呈现出幂函数形式的衰减。

图14 压力峰值拟合曲线

图15 超压比冲量拟合曲线

6 结 论

通过对PMT水下内爆现象的分析，结合气泡动力学理论、数值仿真和算例验证，得到以下结论：

(1) 数值仿真得到的冲击波压力峰值和超压比冲量与实验算例结果吻合较好，数值仿真得到的气泡压溃时间与理论值吻合较好，说明本文提供的数值算法对于容器水下内爆模型具有适用性和准确性。

(2) PMT外形轮廓对冲击波的传播具有一定的影响，PMT水下内爆后向下方传播的冲击波能最大，冲击波在远离气泡传播时呈现幂函数形式的衰减。