近场水下爆炸气泡射流载荷冲击船体外板的动响应分析

姜忠涛， 李 烨， 庞学佳， 王诗平

(1. 哈尔滨工程大学 船舶工程学院， 哈尔滨 150001； 2. 湖南磁浮技术研究中心有限公司， 长沙 410000；3. 中国船舶重工集团公司第七0三研究所, 哈尔滨 150078)

船体结构在遭受近场水下爆炸气泡射流的冲击作用时，首当其冲的应当是其船体外板，船体板架根据不同的水下工况自然会产生不同的弹性或者塑性响应。关于气泡载荷[1]及气泡射流[2]的研究通常从实验和理论两个方面进行， Brunton等[3-6]针对不同类型的材料进行了相应的研究，认为射流载荷下的结构变形是由平板上的冲击流体产生的“水锤”压力及高速横向射流对材料表面进行腐蚀和冲刷作用而产生，引起易变形材料结构表面的塑性凹陷及脆性材料的剪切破坏。

关于船体外板在水下爆炸作用下的动响应多从理论、试验及数值仿真等方面进行研究。Taylor等[7-9]从理论上研究了板承受水下爆炸载荷的结构响；吴成等[10](2007)从理论出发对固支方板刚塑性动态响应进行了分析，并给出相应的运动控制方程和结构变形解析解；张振华等[11]对水下爆炸冲击波作用下舰局部板架结构的动响应的相似规律进行研究；牟金磊等[12]提出了水下爆炸载荷作用下加筋板结构的计算方法。在试验研究方面，Ming等[13]对水下接触爆炸作用下加筋平板的毁伤特性开展了试验研究，并将试验结果与SPH (Smoothed Particle Hydrodynamics)计算结果进行对比；崔杰等[14]针对实尺度舱段模型进行了水下爆炸实验，并对气泡射流载荷的影响范围进行分析；Gifford等[15](1988)结合试验对存在裂纹的焊接厚板进行了近场爆炸分析，并分析了其动态响应。在数值模拟部分，近场爆炸数值模拟方法大多借助于任意欧拉-拉格朗日算法(Arbitrary Lagrangian-Eulerian， ALE) 和耦合欧拉-拉格朗日算法[16]( Coupling Eulerian-Lagrangian， CEL)。王龙侃等[17]采用LS-DGF-DG方法对船体板架结构在近场水下爆炸载荷作用下的毁伤特性进行研究。在具体到内部气泡载荷的分析过程而言，在计算量允许的前提下边界附近气泡的运动形态更多的是采用边界元法[18-19]进行模拟。将气泡射流进行几何等效后借助边界元法计算气泡射流的相关几何数值特征，最后结合CEL或ALE对水射流冲击结构的响应进行模拟计算[20]。

本文采用ABAQUS软件中的CEL算法对气泡射流冲击简单平板的过程进行数值模拟以规避材料的高应变率问题，单元大变形等很多常规问题。

1 基本理论和方法

1.1 耦合欧拉-拉格朗日算法

ABAQUS/Explicit可以实现将划分网格后的欧拉域和拉格朗日域装配到同一个计算模型中。其中射流对应的水域采用欧拉网格进行离散，而平板或者板架结构则遵循传统的拉格朗日网格离散方式。两种网格之间可以交叉，并且可以约束材料防止其进入欧拉域。在欧拉域和拉格朗日域的交界处，后者的网格为前者的网格提供了几何流动边界，而前者的网格为后者网格提供了压力边界。

在ABAQUS中，耦合欧拉拉格朗日算法(CEL)的原理基于中心差分法。求解过程中用位移的形式来表示相应的加速度和速度，其各自的表达式以及运动方程分别如式(1)和式(2)所示

(1)

(2)

将t时刻的速度和加速度结合式(1)与式(2)即可得到中心差分表达式

(3)

式中:M为质量阵;C为阻尼阵;K为阻尼阵；Q为载荷。

(4)

1.2 射流载荷针对CEL算法的等效原则

由于气泡射流速度非常快，且气泡凹陷形成的射流近似为圆台形状，为了便于工程应用，本文在采用CEL算法施加射流载荷时将气泡射流圆台形状近似为圆柱，如图1所示。当气泡以横剖图显示时，S表示射流形状横剖面面积，L表示气泡射流轮廓的高度，D表示射流形状近似后圆柱的直径。表达式为

(5)

图1 气泡射流形状等效方法

在给定初始条件的前提下，采用程序结合边界元算法可直接确定特定工况下圆柱形气泡射流的速度、宽度以及高度，为后续CEL计算射流冲击提供基本参数。

2 数值结果及讨论

2.1 有效性验证

本文采用ABAQUS/Explicit中的CEL算法对射流载荷的相关机理和载荷特性进行分析，势必需要考虑其计算的精度问题。本小节在此以Obara等[21]的射流试验数据以及射流载荷的经验公式为基础，通过在同样的工况条件下与CEL算法的计算结果进行对比验证数值算法的准确性。

Obara等的实验采用高速水射流冲击一定距离下的玻璃钢材料板，通过在特定位置设置压力感应装置获取压力载荷，同时结合高速摄影装置捕捉实验过程中射流冲击产生的冲击波、稀疏波现象。高速水射流通过装置可将射流速度控制在90～1 000 m/s内。试验中射流冲击速度为570 m/s。靶板板材规格为50 mm×50 mm×5.9 mm。

针对CEL算法的数值模拟而言，模型如图2所示。为了与上述试验条件保持一致，平板结构采用实体单元以拉格朗日网格进行离散，拉格朗日域大小为50 mm×50 mm×5.9 mm，水域以欧拉域代替，所建域底部与拉格朗日域重合，高度为20 mm。射流载荷用圆柱体代替，其高度为12 mm，射流区域的直径为3 mm，以570 m/s的初始速度冲向平板，并得到结构和射流流体的动响应。

图2 计算模型示意图

图3 射流域在耦合模型中的位置

图4 板架有限元模型

经验公式、CEL耦合算法以及试验结果的压力时历曲线对比如图5所示。其中试验结果峰值大小为850 MPa，CEL算法结果的峰值大小为817 MPa，三条曲线的峰值相差较小。同时，实验值在驻压段基本维持在180 MPa附近， CEL算法与经验值和试验值相差较小，能够满足相应工程使用精度。

图5 三种算法下射流载荷中心点压力曲线对比

2.2 基于CEL的射流载荷冲击钢板的动响应分析

2.2.1 数值模型

板架材料采用弹性钢板，屈服极限为235 MPa，密度为7 850 kg/m3，弹性模量为2.1×1011Pa，泊松比为0.3，采用拉格朗日网格进行离散。射流部分的密度为1 000 kg/m3，采用圆柱形水柱对射流进行替代，圆柱体直径10 cm，高度50 cm，水柱冲击速度为200 m/s，水的状态方程为线性Us-Up方程，水中声速为1 500 m/s。欧拉域尺寸为156 cm×156 cm×70 cm，整个水域采用154 560个欧拉单元进行离散。耦合模型及平板模型分别如图6和图7所示。

图6 欧拉域与平板耦合模型

图7 加筋平板模型

2.2.2 射流的载荷压力特性

当设定板的厚度为10 mm时，在200 m/s的速度冲击下板架仅发生弹性变形，板架大部分区域应力最大值不超过235 MPa，没有出现塑性变形，计算结果云图如图8所示，且加强筋未出现明显的应力分布，仅仅是加筋结构与面板交接中心区域出现应力集中。而当速度设定为1 000 m/s时在射流中心区域出现部分塑性变形，如图9所示，其局部应力已经超过235 MPa，但因射流半径较小并未达到塑性破坏所需应变，因此没有出现破口。

图8 未出现塑性应变时板架应力云图

图9 出现塑性应变时板架应力云图

射流载荷在水下对舰船结构进行冲击时射流区域的冲击压力的分布一直以来都是业界十分关心的问题。因此在射流半径范围内流体与板架接触的区域选定部分参考点提取其压力进行直观分析，测点位置示意图如图10所示，图中V0表示射流初始速度，Vs表示水流受到平板约束后横向流动后产生的横向速度。

图10 射流中心附近测点示意图

根据图10所示沿径向在射流半径范围内均匀布点，提取其压力时历曲线如图11所示，其中R表示射流半径。根据文献[22]计算射流载荷的经验公式可知，在该速度下“水锤”压力和水动压力经验公式值分别为1.5 GPa和0.5 GPa，从图11可知数值结果的“水锤”压力峰值在1.65 GPa左右，t=0.15 ms后出现压力衰减，随后进入水动压力阶段，图中压力值在0.5 GPa附近振动，与经验公式结果基本吻合。

同时从图11可以看出，当射流以圆柱形冲击平板时，射流半径范围内的外围流体出现压力爬升的时间比中心附近流体较早，即射流中心产生压力的时间较晚。同时，对比水动压力段可以发现，在射流作用的水平压力输出端，射流中心区域附近的压力明显高于外围接触部分，射流中心处的水动压力峰值近似为半径处(r=R)附近区域的2倍。另外，且当水柱对平板作用结束后整个半径区域内的压力迅速衰减到零。

图11 射流半径范围内典型径向接触面节点压力时历曲线

2.2.3 射流的流体速度分布特性

以上述水柱冲击平板为例，选定V0=200 m/s，在流体和结构接触的表面上的部分欧拉单元中设置速度输出参考点，分别以射流中心点为中心沿着径向分散，测点位置与图10相同，并以射流半径为无量纲标准量进行无量纲化，各个参考点的速度时间曲线如图12所示。

图12 射流中心附近欧拉域各径向测点速度时历曲线

从图12中曲线可以看出，在射流与结构接触后，测点速度在极短时间内达到最大值后开始衰减，速度衰减过程中会出现少量的振荡，这种现象可能与平板出现弹性振动有关，但振荡幅值不大。对比射流半径范围内的三个参考点可知，随着远离射流中心点，流体最终的稳定速度也会逐渐变大。而对于射流半径以外的流体，其能达到的速度峰值明显大于射流冲击速度。对于r=1.36R的欧拉节点，其最大速度达到663.78 m/s，甚至超过射流速度的3倍，这一现象与Brunton的研究结果相吻合。随着进一步远离射流中心点，流体速度峰值有所下降，但依然大于原始冲击速度，且随着远离射流中心点距离较大，流体最终趋于稳态流动，速度保持在180 m/s附近，因此距离射流中心较远的节点的最终稳定速度保持一致。射流半径以内的流体最终速度明显小于半径以外区域的速度。

为了直观表示不同射流冲击时刻射流区域及其附近的流体速度在径向的分布情况，在流体径向扩展的方向上设置一系列参考点，将同一时刻下的流体在径向的分布情况如图13所示给出。其中横坐标表示距离射流中心点的无量纲距离，以射流半径表示参考量，由上图可知，在t=0.3 ms时射流中心速度有所增加但中心附近半径范围内的流体速度出现衰减，这一现象与图12的结果向吻合，此时最大速度出现在r=1.5R的环形区域附近，但在最大速度以外区域速度则全部为零。随着时间进一步推进，当t=0.6～0.9 ms时最大速度区域已经位于r=3R附近，但随着远离射流中心，速度峰值出现明显的衰减，且这一阶段射流半径范围内的速度基本保持稳定状态。当t=1.2 ms以后虽然速度进一步向边缘传播，但在r≤3R范围内速度基本保持不变，这一结果也与图12所示结果相同。此后，因为载荷的逐步衰减以及流体的逐渐稳定，速度响应也基本不再出现变化。

图13 不同时刻射流中心附近径向流体速度分布

2.2.4 射流引起的结构剪应力分布特性

当高速射流流体在受到结构板架的阻碍作用而产生横向流动时会产生比射流原始冲击速度更大的速度，甚至一度能达到冲击速度的3～4倍，在结构上出现如此高速的流动势必会在结构表面上产生极大的剪切应力。图14所示为在1 000 m/s射流冲击下平板沿Y方向的剪应力分布云图，由图可知剪应力梯度分别沿从y轴旋转45°方向各自以加强筋为界分成四个区域扩展。其中y轴正向45°为负值，负向45°为正值。由于板架边界设定为简支，当载荷压力结束后，剪应力最大值及应力集中现象最终位于平板四个角落处。

(a) t=0.173 ms

(b) t=0.36 ms

(c) t=0.76 ms

(d) t=1.5 ms

为了更加细致地分析射流冲击平板结构时结构上的剪应力分布，如图15所示在Y轴正向和XY之间的45°方向以射流半径为单位选取一系列参考点，并对比参考点上的剪应力大小。由图16可知45°方向上的剪应力明显大于Y轴正向上应力，且在射流作用的初始阶段会出现一个比较大的瞬时剪应力峰值，且该峰值的出现时间紧随载荷的峰值压力而出现。在峰值过后剪应力则在一个较小的应力值附近振荡，而Y轴正方向剪应力则没有瞬时峰值的出现，即整个冲击过程剪应力都很小。

通过对比不同径向位置剪应力分布可知，r=2R和r=3R位置处45°方向剪应力峰值明显大于r=R处的峰值，且平稳振荡阶段的剪应力也大于半径范围处节点。由此可见，在射流冲击结构后不仅产生速度极高的横向射流，同时在射流径向的某些方位上会对结构产生极高的剪切应力，不同方位的应力之间会产生2～3倍的差别，且射流半径范围外的流体由于其速度更高，产生的剪应力也更高，其对结构表面的材料可能造成十分严重的剥蚀现象。因此，针对水下爆炸气泡射流产生的毁伤效应，不仅需要关注其垂向变形和破口，同时也需要重点关注横向剪切剥蚀毁伤。

图15 方位示意图

(a) r=R

(b) r=2R

(b) r=3R

3 结 论

本文以射流载荷冲击船体外板结构为背景，采用CEL算法对水下爆炸气泡射流载荷冲击简单平板的过程进行了数值模拟。将计算结果与试验进行了对比并分析了相关规律，本文结论具有一般性，具体如下：

(1) 采用CEL算法计算了射流载荷冲击PMMA平板的过程，结合试验和经验公式对比了其载荷曲线，结果表明CEL算法可有效地模拟水射流对平板的冲击过程。

(2) 采用CEL算法模拟了射流载荷冲击弹性钢板的过程，对射流半径范围内结构上的接触载荷压力进行了分析，发现载荷压力在0.6R(R为射流半径)范围内基本保持一定，在0.6R以外区域出现衰减，在r=0.8R附近压力出现较小值。

(3) 通过分析射流半径范围内速度分布发现，流体速度最大峰值出现在半径以外区域(约r=1.36R)，通常能达到射流冲击速度的3～4倍，且超过一定范围后速度峰值随远离射流中心而逐渐变小，但最终速度趋于一致。

(4) 对平板结构上的剪应力分布进行分析，发现在与平板骨材交叉方向存在较大剪应力，约为与骨材平行方向大小的2～3倍，即射流的存在对平板平面形成巨大的剪切应力和剥蚀效应。

参 考 文 献

[1] ZHANG A M, CUI P, CUI J, et al. Experimental study on bubble dynamics subject to buoyancy[J]. Journal of Fluid Mechanics, 2015, 776: 137-160.

[2] FIELD J E, LESSER M B. On the mechanics of high speed liquid jets[C]∥Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. London:The Royal Society, 1977: 143-162.

[3] BRUNTON J H. High speed liquid impact[J]. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 1966, 260(1110): 79-85.

[4] BOWDEN F P, BRUNTON J H. The deformation of solids by liquid impact at supersonic speeds[C]∥Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. London: The Royal Society, 1961: 433-450.

[5] BOWDEN F P, FIELD J E. The brittle fracture of solids by liquid impact, by solid impact, and by shock[C]∥Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. London: The Royal Society, 1964: 331-352.

[6] SMITH D G, KINSLOW R. Pressure due to high-velocity impact of a water jet[J]. Experimental Mechanics, 1976, 16(1): 21-25.

[7] TAYLOR G I. The pressure and impulse of submarine explosion waves on plates[J]. The Scientific Papers of GI Taylor, 1963, 3: 287-303.

[8] COX A D, MORLAND L W. Dynamic plastic deformations of simply-supported square plates[J]. Journal of the Mechanics and Physics of Solids, 1959, 7(4): 229-241.

[9] HUANG H. Transient bending of a large elastic plate by an incident spherical pressure wave[J]. Journal of Applied Mechanics, 1974, 41(3): 772-776.

[10] 吴成, 倪艳光, 郭磊, 等. 水下爆炸载荷作用下气背固支方板的动态响应分析[J]. 北京理工大学学报, 2007, 27(3): 205-209.

WU Cheng, NI Yanguang, GUO Lei, et al. Dynamic response of rectangular air-back plate to underwater explosion loading[J]. Transaction of Beijing Institute of Technology, 2007, 27(3): 205-209.

[11] 张振华，陈平毅，漆万鹏, 等.舰船局部板架结构在水下爆炸冲击波下动态响应的相似率研究[J]. 振动与冲击, 2008, 27(6): 81-86.

ZHNAG Zhenhua, CHEN Pingyi, QI Wanpeng, et al. Scaling law of dynamic response of stiffened plates for a ship subjected to underwater shock[J]. Journal of Vibration and Shock, 2008, 27(6): 81-86.

[12] 牟金磊，朱锡，张振华，等.爆炸冲击作用下加筋板结构变形研究[J]. 海军工程大学学报, 2007, 19(6): 12-19.

MU Jinlei, ZHU Xi, ZHANG Zhenhua, et al. A study on deformation of blast-loaded stiffened plates[J].Journal of Naval University of Engineering, 2007, 19(6): 12-19.

[13] MING F R, ZHANG A M, XUE Y Z, et al. Damage characteristics of ship structures subjected to shockwaves of underwater contact explosions[J]. Ocean Engineering, 2016, 117: 359-382.

[14] 崔杰，张阿漫，郭君, 等.舱段结构在气泡射流作用下的毁伤效果[J]. 爆炸与冲击, 2012, 32(4): 355-361.

CUI Jie, ZHANG Aman, GUO Jun, et al. Damage effect of cabin structure subjected to bubble jet[J]. Explosive and Shock Waves, 2012, 32(4): 355-361.

[15] GIFFORD L N, CARLBERG J R, WIGGS A J, et al. Explosive testing of full thickness precracked weldments[C]∥Fracture Mechanics Twenty First Symposium, ASTM-STP-1074. Philadelphia: ASTM， 1990: 157-77.

[16] 王耀辉, 陈海龙, 岳永威, 等. 水下接触爆炸作用下的船体板架结构毁伤研究[J]. 中国舰船研究, 2012, 7(4): 76-82.

WANG Yaohui, CHEN Hailong, YUE Yongwei, et al. Damage study on ship plate frame subjected to the underwater contact explosion[J]. Chinese Journal of Ship Research, 2012, 7(4): 76-82.

[17] 王龙侃，张之凡，郎济才.基于LS-DGF-DG方法的船体板架结构近场水下爆炸毁伤特性研究[J].振动与冲击,2016,35(16):64-71.

WANG Longkan,ZHANG Zhifan,LANG Jicai,et al.Damage characteristics of hull grillage subjected to near-field underwater explosion based on an LS-DGF-DG method[J].Journal of Vibration and Shock,2016,35(16):64-71.

[18] ZHANG A M, LIU Y L. Improved three-dimensional bubble dynamics model based on boundary element method[J]. Journal of Computational Physics, 2015, 294: 208-223.

[19] ZHANG A M, LI S. CUI J. Study on splitting of a toroidal bubble near a rigid boundary[J]. Physics of Fluids, 2015, 27(6): 809-822.

[20] KLASEBOER E, KHOO B C, HUNG K C. Dynamics of an oscillating bubble near a floating structure[J]. Journal of Fluids and Structures, 2005, 21(4): 395-412.

[21] OBARA T, BOURNE N K, FIELD J E. Liquid-jet impact on liquid and solid surfaces[J]. Wear, 1995, 186(95): 388-394.

[22] COOK S S. Erosion by water-hammer[C]∥Proceedings of the Royal Society of London A: Containing Papers of a Mathematical and Physical Character. London: The Royal Society, 1928: 481-488.