水下近场爆炸冲击波与射流载荷联合作用分析方法

郭 娅，唐文勇，郑绍文，李德聪，苏怡然

（1.上海交通大学 海洋工程国家重点实验室，上海 200240；2.中国舰船研究设计中心，武汉 430064）

水下近场爆炸载荷主要包括爆炸前期冲击波载荷与爆炸后期气泡载荷。局部结构响应主要由二者所致。为此，冲击波载荷、气泡载荷及与结构相互作用已成研究热点。数值模拟方法有基于双渐近近似DAA方法［1］及水下冲击波分析程序［2］，后者对水下爆炸领域影响深远。文献［3－4］对水下爆炸冲击波载荷及对结构毁伤进行研究。而对气泡的研究，自Reyleigh建立不可压缩流场中球形气泡运动方程以来气泡溃灭理论不断完善与发展，张阿漫等［5－7］基于边界积分法对气泡脉动过程进行过深入研究。

理想的水下近场爆炸数值模拟需同时考虑结构受冲击波作用的流固耦合效应及受射流冲击的流固耦合效应；但因软件功能限制尚不能模拟冲击波与气泡射流两阶段完整的流固耦合作用。通常为单独考虑冲击波载荷及射流载荷分别进行流固耦合模拟，再将两计算结果叠加［8］。而实际结构在冲击波后阶段已发生变形或破损，射流载荷实际作用于已变形或破损结构上，因而对结构造成进一步冲击。因此，简单叠加并不能真实反映结构响应结果。

为真实模拟整个水下近场爆炸过程，本文用MSC.Dytran软件据给定爆炸初始参数，利用基于欧拉法的气泡动力学程序模拟气泡近壁脉动过程［9］，获得射流基本参数；基于能量等效原则将射流对结构冲击简化为等效射流载荷。通过编制MSC.Dytran软件子程序，在冲击波数值模拟阶段后自动添加等效射流载荷，从而实现结构受冲击波－射流联合作用的动态响应数值模拟。在整个计算过程中冲击波、射流载荷对结构影响连贯统一，计算结果会更符合实际情况。

1 水下近场爆炸冲击波及射流数值模拟

1.1 水下近场爆炸载荷作用过程

水下近场爆炸过程中按时间顺序结构先后会受三种载荷，即冲击波载荷、气泡脉动压力及气泡射流冲击载荷。冲击波载荷压力峰值高、作用时间短，结构受冲击波作用后会在短时间内产生巨大速度、加速度，导致受损。冲击波作用后气泡周围水体受气泡内压力影响产生压力场，形成气泡脉动压力。气泡脉动后期受壁面影响，在气泡溃灭时会形成高速射流，对已损伤结构造成进一步破坏。据研究［10］，气泡脉动载荷的最大压力小于冲击波压力的10%～20%；据试验、数值仿真研究知，气泡脉动载荷主要引起中远场水下爆炸时结构整体鞭状运动，对近场爆炸的局部损伤影响较小，故本文数值模拟水下近场爆炸时仅考虑冲击波载荷及射流载荷。

1.2 冲击波阶段数值模拟

MSC.Dytran求解流固耦合问题时需定义固体或流体压力及密度与比内能间之关系，即状态方程。常规水下爆炸情况用TAIT状态方程描述水体为

水下爆炸，水与爆炸产物气体将在欧拉网格中共存。为描述该现象，应使用多物质欧拉单元。炸药及流体的初始分布、条件参数可通过多物质欧拉单元初始状态定义，见表1。流场边界均定义为水介质可流出的无反射边界。

表1 冲击波数值模拟模型欧拉网格初始状态Tab.1 Initial conditions for the Euler mesh in blast wave simulation model

1.3 射流阶段数值模拟

求解水下气泡脉非球状脉动过程的数值算法主要有基于势流理论的边界积分法、基于NS方程的欧拉法。边界积分法计算效率高，但适用范围有限；欧拉法适用范围广，但计算效率低。本文采用改进的欧拉法模拟气泡近壁脉动过程［9］。该方法在欧拉法基础上提出，可有效提高计算速度的近场压力边界条件。将近场压力边界外部流体设为理想流体，通过势流理论求解。外部流体为内部流体提供动态压力边界条件，实现欧拉域减小、有效提高欧拉法模拟水下气泡动态特性的计算效率。

图1 射流时刻气泡形状简化描述Fig.1 The simplification of the bubble shape at jet time

据给定的爆炸初始参数，采用改进的欧拉法编制气泡动力学程序模拟气泡近壁脉动过程，获得气泡脉动过程中流场压力分布及气泡射流基本参数。与射流冲击载荷计算相关参数为射流形状、射流速度及射流冲击时刻气泡离结构距离。考虑射流形状、尺寸参数时可将射流简化为凸台［8］，用3个尺寸描述射流形状，见图1。图中H为凸台高度；d1，d2分别为凸台上、下表面直径；Dist为射流距壁面距离；Vel为射流产生冲击时刻速度。

2 结构受冲击波－射流联合作用数值模拟方法

考虑冲击波载荷与射流载荷联合作用时进行冲击波作用于结构的流固耦合模拟。据1.3节所得射流基本参数进行单独射流冲击流固耦合数值模拟，获得结构不同位置受冲击载荷的时历曲线及结构变形能W0。由于冲击压力时历曲线波动性较大，故需适当简化方可用于Msc.Dytran二次开发加载。压力曲线简化方法为

设在射流冲击时对特定时间步及板格其等效射流冲击压力为恒定值，式（4）中压力简化主要思想为在时间、面积上对压力取平均，对给定时间区间Ti及面积区域（环形）Sj，等效射流压力记Pi，j。式（4）积分计算为

获得等效射流载荷后施加于结构进行流固耦合分析，获得结构变形能W1。将W1与射流水柱流固耦合模型所得变形能W0进行对比，两次所得变形能应保持一致。若两者相差较大，则修正等效射流载荷大小，直至W1与W0相等，此时等效射流载荷对结构作用效果与射流水柱直接作用效果应一致。

计算模型以冲击波流固耦合模型为基础，通过调用Msc.Dytran用户子程序完成结构等效射流冲击载荷施加。计算过程中可通过用户子程序标记失效单元标号，统计失效区域面积。计算流程见图2。

图2 结构受冲击波－射流联合作用的数值模拟流程Fig.2 Process for the numerical simulation of structure subjected to both explosive wave and jet impact

Msc.Dytran用户子程序用FORTRAN语言编写，其中EXFUNC函数允许用户定义随时间变化的外部函数，可用于描述随时间变化载荷。Msc.Dytran通过TABLEEX选项卡在每个时步内多次调用EXFUNC函数对结构加载。同一模型可对结构不同受载区域定义不同TABLEEX选项卡，该选项卡通过不同标识名区别彼此。在用户子程序中通过判断结构单元应变是否超过许用应变筛选破损单元，并统计破损区域及塑性区域面积。

3 冲击波－射流联合作用数值模拟实例

3.1 冲击波－射流联合作用结构响应数值模拟

图3 板架结构有限元模型Fig.3 The frame structure finite element model

据上节计算流程计算典型舷侧板架结构遭受水下接触爆炸作用后的动态响应。在水深5 m、距板架5 m处装500 kg当量TNT炸药，距离参数（爆距与气泡最大半径比值）为0.5。板架单层结构长宽均6 m，板厚12 mm，长宽方向各设用5根T材⊥加固，见图3。结构单元网格尺寸100 mm×100 mm。板结构弹性模量210 GPa，泊松比0.3，密度7.8 g／cm3，屈服强度440 MPa；材料破坏准则采用最大等效应变失效准则，即等效塑性应变超过0.3时结构发生破坏［3］，或取 0.25［12］，本文用 0.25作为临界失效应变。应变率强化模型用CowperSymond模型［4］，即

据爆炸初始参数用气泡动力学程序对气泡脉动过程进行数值模拟见图4，获得射流基本参数见表2。由表2知，射流时刻气泡贴在结构壁面上直接对结构造成冲击，射流速度120 m／s。经比较，射流参数与经验、试验数据规律相符［13－16］。

表2 水下爆炸近壁射流数值模拟计算结果Tab.2 Numerical results of the nearwall jet caused by underwater explosion

图4 结构受射流冲击作用的数值模型Fig.4 Numerical model for the simulation of the structure under jet impact load

射流流固耦合模型耦合面由结构单元与哑元构成。用多物质欧拉单元，欧拉区域初始状态见表3。

区域 状态方程密度／（kg·m－3）比内能／（J·kg－1）优先级全欧拉域 TAIT 1000 8.4×1041射流水柱 TAIT 1000 8.4×1042

在板架结构距中心3 m内由中心向外围均匀取30个监测点进行压力监测，记为P1～P30。部分监测点压力时历曲线见图5，可见压力变化规律与经验相符［8，15］。

图5 板架结构监测点压力时历曲线Fig.5 Time history graph of the pressure in several monitoring points

图6 等效射流载荷加载区域示意图Fig.6 Equivalent jet load loading area

图7 水下接触爆炸作用下板架结构吸能曲线Fig.7 Time evolution of structural energy for the shell structure

射流冲击时间取0.2 ms，以0.5 m为半径，最终得板架结构由中心向外5层圆环区域等效射流载荷分别为：第一层81 MPa，第二层73 MPa，第三层65 MPa，第四层4.8 MPa，第五层 0.07 MPa，其它区域载荷较小，可忽略，见图6。

获得等效射流冲击压力随时间、空间分布后，用上节用户子程序方法将该载荷数据添加到冲击波流固耦合模型中，实现结构受冲击波－射流联合作用的数值模拟。结构受冲击波作用的流固耦合数值模型与射流作用的流固耦合数值模型类似，耦合面亦为由外壳拉格朗日单元及哑元组合的闭合面，哑元不提供刚度。

3.2 数值模拟计算结果

用MSC.Dytran软件计算3.1节数值模型，获得板架结构在水下近场爆炸载荷作用下的动力响应结果。板架结构在冲击波－射流联合作用下的结构变形能及动能随时间变化曲线见图7，爆炸过程中结构等效塑性应变见图8。由图7看出，结构变形能有两次大幅度提升，分别为冲击波作用阶段与射流冲击阶段，且结构在冲击波作用阶段吸收的变形能累积到射流阶段；结构在射流阶段吸收的变形能小于冲击波作用阶段；结构动能经历两次先增大后减小过程：第一次结构动能增加因冲击波作用所致，且在阻尼作用下逐渐恢复静止；第二次动能增加因射流作用所致，且射流载荷对结构冲击小于冲击波作用。由图8可知，冲击波作用的板架结构等效塑性应变最大为0.157，射流冲击作用的板架等效塑性应变为0.175，因此若不考虑射流载荷影响，对结构评估结果偏危险。

由计算结果知，板架结构受冲击波瞬间作用而获得较大动能。虽冲击波作用时间极短，但冲击作用后10 ms内板架动能持续转化为结构变形能。图8（a）为结构受冲击波作用后板架的等效塑性应变云图，板架结构以小板格为单位发生局部弯曲，塑性变形在边界附近略高，此由结构边界约束偏刚性导致。射流冲击载荷的实际作用时间约在爆炸后1 s。由计算知冲击波作用后约15 ms时结构塑性变形能已趋于稳定，动能趋于零，可认为结构破坏与变形已充分发展。为节约计算时间可将射流冲击作用提前至冲击波作用后15 ms。板架结构受射流冲击作用后响应过程与冲击波阶段相似，即板架获得较大动能后逐渐转化为变形能。板架受射流冲击中心变形、应变进一步增大。图8（b）～（d）为板架受射流冲击作用后的变形与破坏形式。

图8 水下近场爆炸作用板架结构等效塑性应变云图Fig.8 Fringe figure of the effective plastic strain for the shell structure of underwater nearfield explosion

4 结 论

本文利用基于欧拉法的气泡动力学程序与Msc.Dytran软件对结构受冲击波－射流联合作用的动态过程进行数值模拟，通过算例对计算结果进行分析，结论如下：

（1）近场水下爆炸时射流冲击载荷不可忽略，射流冲击会使结构变形、应变进一步加剧，进而导致受冲击波作用未损坏结构在射流作用下发生破坏。

（2）本文综合考虑冲击波载荷、射流载荷作用，简化处理射流载荷，使其在冲击波流固耦合模拟后期自动添加等效射流冲击载荷，能较真实模拟水下近场爆炸过程，适合工程应用。

［1］Geers T L.Doubly asymptotic approximations for transient motions of submerged structures［J］.The Journal of the Acoustical Society of America，1978，64（5）：1500－1508.

［2］DeRuntz Jr J A.The underwater shock analysis code and its applications［C］.／／60th Shock and Vibration Symposium，SAVIAC，1989：89－107.

［3］陈永念，谭家华.固支方形板在水下爆炸冲击波载荷作用下动态响应数值仿真研究［J］.舰船科学技术，2008，30（1）：131－135.CHEN Yongnian，TAN Jiahua.Dynamical response of square plates under shock load［J］.Journal of Ship Science and Technology，2008，30（1）：131－135.

［4］袁建红，朱锡，张振华.水下爆炸载荷数值模拟方法［J］.舰船科学技术，2011，33（9）：18－24.YUAN Jianhong，ZHU Xi，ZHANG Zhenhua.Numerical simulation method study of underwater explosion load［J］.Journal of Ship Science and Technology，2011，33（9）：18－24.

［5］张阿漫.水下爆炸气泡三维动态特性研究［D］.哈尔滨：哈尔滨工程大学，2006.

［6］姚熊亮，汪玉，张阿漫.水下爆炸气泡动力学［M］.哈尔滨：哈尔滨工程大学出版社，2012.

［7］宗智，何亮，孙龙泉.水下爆炸气泡对水面舰船载荷的数值研究［J］.船舶力学，2008，12（5）：733－739.ZONG Zhi，HE Liang，SUN Longquan.Numerical study of loading on the ship near an underwater explosion bubble［J］.Journal of Ship Mechanics，2008，12（5）：733－739.

［8］崔杰，李世铭，黄超，等.射流冲击模型在水下爆炸实验室中的应用研究［J］.高压物理学报，2012，26（5）：523－530.CUI Jie，LI Shiming，HUANG Chao，et al.Applied research of jet impact model in underwater explosion experiment［J］.Chinese Journal of High Pressure Physics，2012，26（5）：523－530.

［9］苏怡然.水下接触爆炸作用下双层防护结构瞬态动力分析方法研究［D］.上海：上海交通大学，2013.

［10］姚熊亮，陈建平，任慧龙.水下爆炸气泡脉动压力下舰船动态响应分析［J］.哈尔滨工程大学学报，2000，21（1）：1－5.YAO Xiongliang，CHEN Jianping，REN Huilong.The analysis of dynamic response of ship hull subjected to gas bubble impulsive pressure of underwater explosions［J］.Journal of Harbin Engineering University，2000，21（1）：1－5.

［11］Popinet S，Zaleski S.Bubble collapse near a solid boundary：a numerical study of the influence of viscosity［J］.Journal of Fluid Mechanics，2002，464：137－163.

［12］李海涛，朱锡，张振华，等.弹塑性边界下水下爆炸气泡动态特性的数值仿真研究［J］.系统仿真学报，2009，21（14）：4532－4535.LI Haitao，ZHU Xi，ZHANG Zhenhua，et al.Numerical simulation study of pulsation dynamics of underwater explosion bubble located under an elastoplastic plane［J］.Journal of System Simulation，2009，21（14）：4532－4535.

［13］张阿漫，王诗平，白兆宏，等.不同环境下气泡脉动特性实验研究［J］.力学学报，2011，43（1）：71－82.ZHANG Aman，WANG Shiping，BAI Zhaohong，et al.Experimental study on bubble pulse features under different circumstances［J］.Chinese journal of Theoretical and Applied Mechanics，2011，43（1）：71－82.

［14］汪斌，张远平，王彦平.水中爆炸气泡与水底边界相互作用的水射流现象［J］.爆炸与冲击，2011，31（3）：251－256.WANG Bin，ZHANG Yuanping，WANG Yanping.Water jet phenomena induced by the interaction between underwater explosion bubbles and water bottom boundaries［J］.Explosion and Shock Waves，2011，31（3）：251－256.

［15］张伟，岳永威，张阿漫，等.基于AUTODYN的气泡与固定壁面相互作用数值模拟［J］.中国舰船研究，2012，7（6）：24－32.ZHANG Wei，YUE Yongwei， ZHANG Aman， et al.Numerical simulation of underwater explosion bubble interactions with a solid boundary based on Autodyn［J］.Chinese Journal of Ship Research，2012，7（6）：24－32.

［16］张阿漫，姚熊亮.水深和药量的变化对水下爆炸气泡射流的影响研究［J］.工程力学，2008，25（3）：222－230.ZHANG Aman，YAO Xiongliang.The effect of charge and water depth on the underwater explosion bubble［J］.Engineering Mechanics，2008，25（3）：222－230.

［17］梅志远，朱锡，侯海量.双层防爆舱壁结构动态吸能特性分析［A］.无锡：船舶力学学术会议暨《船舶力学》创刊十周年纪念学术会议论文集，2007.