约束空间内爆超压载荷影响因素二维数值模拟

徐维铮，吴卫国

1武汉理工大学高性能舰船技术教育部重点实验室，湖北武汉430063

2武汉理工大学交通学院，湖北武汉430063

0 引 言

海洋是支撑未来发展的战略制高点，已成为世界各国竞相发展的领域，导致各国间的海洋权益摩擦、领海和岛屿主权争端不断。大型水面舰船是海上国防力量建设最为关键的部分之一，然而其在执行作战任务时，不可避免地面临着日益先进的反舰武器的打击，包括水下威胁（鱼雷、水雷等）和水上威胁（反舰导弹、自杀式飞机等）。

目前，在反舰导弹的研发和制造中，半穿甲爆破型战斗部被各国广泛采用，如美国的“鱼叉”、法国的“飞鱼”、德国的“鸬鹚”以及俄罗斯的“白蛉”等反舰导弹。从攻击模式和毁伤机理来看，采用半穿甲爆破型战斗部的反舰导弹距海面3～8 m掠海飞行，凭借弹体的动能侵彻舰船舷侧外板，同时触发延迟引信，穿入舰船舱室内部爆炸［1］。舱内爆炸由于局限在约束空间内，战斗部爆炸形成的冲击波经多次壁面反射及叠加，会形成持续时间较长的准静态压力，将对舱室内人员、设备及结构产生较大的毁伤。因此，深入研究约束空间内爆炸的超压载荷及其影响因素具有重要意义。

针对约束空间内爆炸的准静态超压载荷问题，国内外学者均有所研究。例如，Anderson等［2］采用量纲分析理论，在大量实验数据的基础上，提出了约束空间内爆炸的准静态超压载荷计算公式；Feldgun等［3］基于伯努利方程，推导出了泄压过程中的准静态超压公式；徐维铮等［4-6］基于自主开发的三维内爆炸波高精度数值计算程序，探讨了泄压口、WENO格式（Weighted Essentially Non-oscillation Scheme）的精度、炸药量等因素影响约束空间内爆炸准静态超压载荷的规律；Xu等［7］在文献［4，6］研究的基础上，进一步通过量纲分析理论给出了计算准静态超压载荷的简化公式。然而，上述研究尚存在一定的缺陷，研究结果与真实的爆炸过程也存在某些差异。

从真实的爆炸过程来看，在载体搭载下战斗部处于运动状态，即使其通过舰船舷侧穿甲进入舱内，仍保留了一定的剩余速度。从近期国内相关文献来看，有学者针对空中动态爆炸（动爆）超压场进行了研究。例如，聂源等［8］研究指出，由于炸药随载体运动，故对爆炸后产生的冲击波超压进行测试较为复杂，目前主要以地面静态爆炸（静爆）实验测量为主，在动爆与静爆冲击波场有着明显差异的情况下，静爆实验无法精确测量出战斗部的实际威力，所以研究炸药动爆冲击波超压场及其计算模型非常有必要。鉴于此，为了获得球形炸药动爆冲击波超压场的计算模型，在计算静爆冲击波超压的Baker公式中加入修正因子进行修正，建立了包含炸药运动速度、对比距离和方位角的修正因子函数的方法；蒋海燕等［9］利用AUTODYN商用软件，对炸药在不同运动速度下的空中爆炸（空爆）冲击波场进行数值计算，以定量研究运动状态下炸药空爆冲击波场的特性，分析了动爆与静爆冲击波场的关联特性，建立了动爆冲击波超压工程模型，其计算结果与动爆实验及仿真结果吻合较好。然而，对于炸药在运动状态下舱内爆炸的准静态超压载荷的影响规律鲜有研究。

值得注意的是，导弹战斗部壳体破裂后，会对爆轰产物阵面产生一定的初始扰动，因此在进行等效裸药数值模拟过程中，需要给予初始界面一定的随机扰动，文献［10-11］在球形炸药爆炸数值模拟过程中，均在初始界面处添加了一定的随机扰动，使数值计算更符合实际，但有关添加的初始扰动对内爆炸超压载荷的影响规律仍需进一步探讨。此外，在数值模拟过程中，有关网格尺寸对内爆炸超压载荷的影响规律也需考察。

综上所述，在当前的数值模拟研究中，对于炸药运动、爆炸初期界面的初始扰动及网格尺寸这些因素对内爆炸超压载荷的影响规律关注尚少，因此，本文拟通过自主开发的高精度二维数值计算程序对上述3种因素进行初步探讨，为真实爆炸过程中的抗爆设计和毁伤评估提供一定的指导及依据。

1 控制方程及数值求解

基于自主开发的约束空间内爆炸波高精度数值计算程序，通过求解二维可压缩欧拉方程［12］，分别采用五阶WENO有限差分格式和三阶TVD龙格—库塔（TVD-RK）法，对其空间项及时间项进行数值离散［13］。

2 内爆炸载荷数值计算

本节主要采用自主开发的计算程序开展不同内爆炸工况下的数值模拟研究，探讨炸药运动、炸药界面初始扰动及网格尺寸对舱内爆炸超压载荷的影响规律。

2.1 炸药运动的影响

为了探讨炸药运动对舱内爆炸超压载荷的影响，假设炸药自左向右沿舱室长度方向运动，并在舱内中部瞬时起爆。在数值模拟中，炸药运动速度v0分别取为0，600和1 200 m/s。

2.1.1 舱室尺寸及爆炸初场

图1所示为数值模拟计算用的长方形内爆炸舱室示意图。本文在瞬时爆轰假定［14］的基础上，将初始的方形炸药等效为均匀的高压爆轰产物，其边长为160 mm，密度为1 630 kg/m3，压力为3.057 9×109Pa。周围的空气密度为 1.0 kg/m3，压力为1.0×105Pa。舱室壁面边界条件设置为刚性边界（不考虑舱室结构变形），泄压口处的边界设置为透射边界［13］。计算网格数为180×80。在舱室右壁面中心处设置一个典型的测点1，以输出内爆炸超压载荷的时程曲线。

2.1.2 爆炸压力场分布

图2所示为典型时刻下3种运动速度的炸药在密闭舱内爆炸的压力场分布云图。由图可以看出，静爆压力场分布沿舱室长度方向呈对称分布，炸药运动状态改变了舱内压力场的分布形态，破坏了对称性；而对于动爆压力场，炸药运动状态改变了舱内压力分布形态。

2.1.3 爆炸载荷分析

图3所示为不同运动速度的炸药在舱室（密闭舱和泄压舱）内爆炸后壁面测点1处的超压载荷时程曲线。由图可以看出，炸药运动速度对冲击波峰值及到达时刻的影响较大，随着炸药运动速度的增大，冲击波峰值明显增大，初始冲击波峰值到达时刻提前。

为定量分析图3中静爆和动爆冲击波的差异性，并验证所开发数值计算程序对动爆冲击波计算的可靠性，参考文献［9］，利用AUTODYN商用程序对不同运动速度下炸药空中爆炸的冲击波场进行数值计算，在综合分析动爆与静爆冲击波场的关联特性的基础上，建立了动爆冲击波超压的工程计算模型。动爆冲击波超压峰值pp，d与炸药运动速度v0、静爆冲击波超压峰值pp，s的关系表达式为

式中：R为冲击波阵面至爆心的距离；c0为波阵面前的空气声速；α为冲击波阵面与炸药爆心的夹角。

由式（1）可知，增大炸药运动速度，将增大向右方向运动的冲击波强度。针对本节算例，在3种炸药运动速度下，分别测得内爆炸壁面测点1的初始反射冲击波峰值为184.4，327.8，485.1 MPa。根据强冲击波壁面正反射关系pr/pi=(3γ-1) /(γ-1)（其中pi为入射冲击波峰值，pr为反射冲击波峰值，γ为气体绝热指数），分别得到入射冲击波峰值为 23.05，40.98，60.64 MPa。由上式，分别得到600，1 200 m/s运动速度下炸药的动爆初始入射冲击波峰值为37.68，55.90 MPa，其数值解与拟合公式的相对误差分别为-8.05%，-7.82%。这初步验证了使用本文数值计算程序计算动爆冲击波的可靠性。

根据能量守恒定律，推导出考虑炸药运动初速的密闭空间内爆准静态超压峰值Δps的计算公式为［6］

式中：m为炸药质量；V为舱室体积；，为爆炸过程中释放的总能量（含运动状态下炸药动能），其中QTNT为炸药爆热，QTNT=4.69×106J/kg；p0为初始大气压力；ρE为炸药密度，ρE=1 630 kg/m3。

由式（2）可知，增大炸药运动速度将使舱内最终形成的内爆炸准静态压力增加，炸药初始动能最终转化为舱内气体内能。由式（2），得3种运动速度下对应的内爆炸准静态超压峰值分别为60，62，69 MPa，最大增加幅值约15%。然而，由于图3中爆炸初期前几个冲击波峰值很大，使得准静态超压峰值间的差异不能很好地体现。为此，本文单独给出了3种运动速度下炸药在密闭舱内爆炸的准静态超压峰值与理论解的对比，如图4所示，从图中可以明显看出两者的差异。

2.2 界面初始扰动的影响

为了探讨界面初始扰动对舱内爆炸超压载荷的影响，本文选取了圆形炸药进行数值计算研究。在保证炸药质量一定的前提下，圆形炸药界面初始扰动采用如下方式添加［15］：

式中：r为圆形炸药的半径；r0为无扰动时炸药的半径；ξ为扰动幅值；N为扰动波数；θ为界面上一点与圆心的夹角。选取2种扰动波数（8和16）进行数值计算。本算例中，ξ的取值为10 mm。

2.2.1 爆炸初场

本文在瞬时爆轰假定［14］的基础上，将初始圆形炸药等效为均匀高压爆轰产物，其半径为80 mm，密度为 1 630 kg/m3，压力为 3.057 9×109Pa。计算的舱室尺寸、炸药位置、周围空气参数、边界条件及测点设置与3.1.1节相同。对于爆炸初场计算，网格数取为360×160。在该算例中，针对泄压舱，泄压口设置在舱室上壁面中心处，尺寸取为200 mm。图5所示为不同界面扰动下的炸药初始形态。

2.2.2 爆炸压力场分布

图6所示为典型时刻下不同初始界面扰动时炸药在泄压舱内爆炸的压力场分布云图。从图可以看出，界面初始扰动主要改变了舱内压力分布的局部形态，对舱内整体压力分布形态的影响较小。

2.2.3 爆炸载荷分析

图7和图8分别给出了不同初始界面扰动下密闭、泄压舱壁面典型测点1处的超压载荷时程曲线。从图中可以看出，不同初始界面扰动下超压时程曲线具有相似的变化规律，主要差异在于冲击波峰值的大小及相位。由于求解的方程为双曲型强非线性方程，即使是微小的初始界面扰动，均可使得爆炸冲击波的峰值与相位产生差别。这说明界面初始扰动主要影响瞬态冲击波载荷，而对于持续时间较长的准静态超压载荷的影响很小。

从图中还可以看出，当初始扰动波数N=8时，超压峰值数值较大；当初始扰动波数N=16时，超压峰值与无界面扰动结果的接近程度增加。这说明，随着扰动波数的增大，界面上的总体扰动趋于均匀化，使得计算结果与无扰动时的更接近。若增加初始扰动波数至N=32，计算得到的局部放大图，如图7（b）和图 8（b）所示。从这两个图来看，同样可说明上述影响规律。

由上述分析可知，对于某爆炸工况，存在某一特定的扰动波数N可给出较高冲击波峰值的计算结果。

2.3 网格尺寸的影响

这里需要说明的是，在数值计算程序中，网格对炸药填充质量的影响较大，在分析网格对爆炸载荷的影响规律时，需要保证炸药质量一定。鉴此，本文在保证炸药质量一定的前提下，选取了3种网格数（80×80，160×160，320×320）进行数值计算。

2.3.1 舱室尺寸及爆炸初场

由于网格数增大会使计算时间迅速增加，所以在计算舱室的爆炸初场时，将舱室选为正方形，如图9所示。在瞬时爆轰假定［14］的基础上，本文将正方形炸药等效为均匀高压气团，其边长为20 mm，密度为1 630 kg/m3，压力为3.057 9×109Pa。炸药位置、周围空气参数、边界条件及测点设置与3.1.1节相同。

2.3.2 爆炸载荷分析

图10所示为不同网格尺寸下密闭舱和泄压舱壁面典型测点1处的超压载荷时程曲线。

从图10可以看出，随着网格尺寸的减小，冲击波峰值增大（尤其是初始冲击波峰值和第2次冲击波峰值），到达时间提前，不过网格尺寸对形成的准静态压力载荷影响很小。

为了定量分析冲击波峰值随网格数变化的规律，图11给出了壁面测点1处的超压载荷时程局部放大图和初始及第2次冲击波峰值随网格数的变化曲线。

这里需要说明的是，当网格数取为640×640时，计算非常缓慢，若计算到6 ms，需消耗大量时间。因此，图11中只在对比初始和第2次冲击波峰值时给出了网格数为640×640的计算结果，且仅计算到t=0.20 ms。

分析图11可知，舱壁面反射冲击波载荷的计算与网格尺寸相关性很大，随着网格加密，反射冲击波载荷迅速增加，且收敛速度较慢，不仅如此，计算时间迅速增加。因此，在数值计算时不能将网格划分得过细。由于冲击波阵面非常陡峭，在数值计算中会跨越2～3个网格节点，至于采用多大的网格尺寸既能捕捉冲击波阵面，又能保证峰值不失真，目前还没有合适的准则来判断，作者也将在后续研究中对此进一步探讨。

3 结论

通过本文的研究，得到如下主要结论：

1）炸药运动改变了舱内爆炸压力场的分布形态，增大炸药运动速度能较明显地提高冲击波峰值，但对内爆炸准静态超压载荷的影响较小。

2）界面初始扰动主要影响瞬态冲击波峰值，而对持续时间较长的准静态超压载荷的影响很小；对于内爆炸工况，存在某一特定的扰动波数可给出较高冲击波峰值的计算结果。

3）网格尺寸对冲击波峰值及到达时间的影响较大，但对形成准静态压力载荷的影响很小。随着网格尺寸的减小，冲击波峰值增大，到达时间提前。尤其是在进行结构与冲击波耦合响应计算时，需要专门对网格进行对比分析，以选取合适的网格尺寸。