水下爆炸冲击波作用下船体刚体动态响应特性

王军，郭君，孙丰，杨棣

(哈尔滨工程大学 船舶工程学院，黑龙江 哈尔滨 150001)

舰船在水下爆炸环境中主要发生2种动态响应:刚体运动及自身结构弹塑性变形破坏。若船体的刚体运动不能进行正确分析，则后续舰船的弹塑性动态响应计算及强度评估等工作便无法准确进行［1-2］。关于水下爆炸冲击波作用下船体动态响应的研究开展了很多，但大部分是结构的弹塑性响应分析［3-5］，船体的刚体运动响应的影响一般被忽略。李海涛等［6］提出了一种计算船体梁在球面冲击波作用下发生整体运动响应的理论方法，是针对试验模型的船体梁开展的研究，炸药在船体结构中部正下方，将球面波看作是移动的平面波，认为船体梁横向截面受到的冲击波压力为均匀分布。张弩等［7］研究了船体在水下爆炸气泡作用下的刚体运动对其鞭状运动的影响，并未考虑冲击波阶段的影响。

针对目前研究进展，考虑到实船可能受到任意位置水下武器的攻击，同时冲击波是以球面冲击波的形式作用在船体上，船体不同位置受到作用力的时间不同。为此在文中，首先分析船体在水下爆炸冲击波作用下的刚体运动，即船体在水中爆炸发生后的起伏运动，然后将刚体运动和弹塑性振动分开研究。建立的理论模型以球面冲击波的形式进行加载，并考虑冲击波的反射、流场和船体的流固耦合作用，分析药包在水下任意位置船体的刚体动态响应。之后通过数值模拟的方法计算出冲击波入射能，得出船体刚体运动的动能及船体结构弹塑性变形等吸收的能量。此过程应用动量定理进行求解，从而间接求出船体振动变形等吸收的能量，进而对船体的弹塑性振动变形等进行分析。

1 船体垂向刚体运动计算模型

在计算舰船受到水下爆炸冲击波的作用时，由于不同的船体外壳其计算模型不同，为理论计算带来很大困难。文献［9］中给出了以流固耦合的Taylor平板理论计算船体动态响应的合理性，因在已有理论基础上，应用Taylor平板来代替复杂形状的船体外壳，分析船体在冲击波作用下的整体垂向运动。

对于铸装TNT球形药包，空间任意一点爆炸冲击波的压力随时间和位置的变化为

式中:Pm为峰值压力;θ为衰减时间常数，对于板上不同的位置而言，其距药包的距离不同，所以爆距不同。设药包位置为O(xW，yW，zW)，水中声速为c，板上离药包最近点为A(x0，y0，z0)，如图1所示，则冲击波到达平板的时间:

以冲击波开始作用平板时记为初始时刻，以图1所示坐标系进行计算，则对于板上任意一点B(x，y，z)，相对于A点受到压力有一时间延迟:

在td时刻之前，B点并未受到水下爆炸冲击波压力，在此定义一个便于计算的辅助函数:

此时平板任意一点受到的球面冲击波入射压力可以表示为

得到冲击波入射压力后，计算冲击波反射压力时采用流固耦合的Taylor平板理论，假定入射波全部反射，并且这种假定相当近似于紧接冲击后的实际情况［9］。考虑波以与板面法向成α角度冲到钢板上，如图2所示，对于空气背衬板，流固耦合界面处的压力可表示为

式中:Pr(x，y，t)为反射压力，应用Taylor平板理论求解流固耦合界面处的压力为

图1 球面冲击波作用于刚性平板Fig.1 The spherical shock wave acting on a rigid plate

图2 平板和流体相互作用Fig.2 Interaction between plate and fluid

整个平板受到球面冲击波的作用力随时间的变化可表示为

入射角α(x，y)也是随平板上作用点变化而发生变化的量，用坐标可表示为

Taylor平板理论没有正确模拟流体惯性(水动力质量)的影响，在此研究平板的垂向运动，需考虑附加质量的影响。为得到较准确的理论计算结果，不使用空气背衬单面浸水的平板附加质量计算公式，而直接利用流固耦合动力学中的边界元法［10］计算船体的垂向平动附加质量，计算公式为

式中:φ为速度势，S为所有边界面，ρf为流场密度。则由动量定理，平板所受垂向力与自身垂向速度关系为

式中:m为船体自身质量。联立式(8)～(11)，则平板刚体运动速度随时间变化为

上式含有随位置和时间变化无积分性质的辅助函数δ(x，y，t)，同时等式两边都有要求解的未知函数v(t)，因此无法得到解析解。但将平板面积进行微分，计算时间离散，通过迭代方法，可以方便进行求解。因平板在受到爆炸冲击波之前，垂向速度为零，初始冲量为零，所以求解思路可表示为

式中:大括号内表示平板整体在某时刻受到的压力，Ii表示冲量对时间的积分值，Δx、Δy、Δt分别为板长度、宽度和计算时间的微元，根据求解精度及计算时间的需要可进行设置。

2 理论模型的试验验证

2.1 舰船实尺度舱段试验验证

为验证本文理论计算方法的合理性，通过某舰船舱段水下爆炸试验的高速摄影试验资料，得到试验舱段在水下爆炸冲击波作用下的最大速度，将试验值与理论计算值进行比较。试验舱段模型长度方向取某舰船的5倍跨长L，宽度方向取船半宽B，高度方向取四甲板高度H，具体结构如图3所示。

图4为试验时舱段在水下爆炸响应时的图片，试验过程中设置3个工况，如图5所示，都采用8 kg的TNT炸药，药包位于舱段中部下方，具体位置如表1所示。通过高速摄影技术，可测出舱段在整个水下爆炸过程中的刚体位移响应，将位移微分便可得到舱段的速度时历响应。

图3 舱段结构尺寸Fig.3 Structural size of section model

图4 水下爆炸试验Fig.4 The test of underwater exp losion

图5 工况设置Fig.5 Structural size of section model

表1 试验工况Table 1 Test conditions

不同工况得到的试验值和理论值的速度最大值如表2所示，爆距越小，舱段的刚体运动速度越大，总体上理论值和试验值比较接近，误差控制在10%以内，充分说明了本文理论计算方法在预报船体水下爆炸刚性运动的有效性。其中工况1的误差最大，可能由于爆距较近时在试验过程中模型发生了一定程度的旋转，增大了试验值测量的误差。

以工况1为例，计算出的舱段整体受到的冲击波垂向压力如图6，在求解空气背衬板的反射压力时应用了Taylor平板理论，此理论没有考虑气泡脉动、气泡射流等影响，但已经能成功模拟经受水下爆炸时船体结构在早期运动的动态响应［9］。本文在此不考虑气泡脉动对舰船刚体运动的影响，以此求出舱段刚体运动速度如图7所示，出现最大刚体运动速度的时间大约为1.5 ms，在8 ms时速度趋近于0。

表2 试验值与理论值对比Table 2 Comparison of experimental and theoretical results

图6 舱段受到冲击波压力随时间变化曲线Fig.6 Variation curve of the section receiving the shock wave pressure with time

图7 舱段刚体运动速度随时间变化曲线Fig.7 Variation curve of the section velocity of rigid movem ent with time

2.2 实船试验验证

为验证本文方法应用于实船在水下爆炸冲击环境中计算刚体运动速度的有效性，利用文献［11］中的实船水下爆炸试验数据进行对比。试验船体长54 m，宽8.5 m，排水量570 t。应用本文方法计算得到的理论值为2.07 m/s，而实船试验值为2.17 m/s，相对误差5%。

实船的受力和刚体速度响应变化趋势和舱段类似，但实船受到的冲击波压力比舱段大很多，作用时间也长，同时最大速度出现的时间比舱段晚，以冲击波开始作用船体为初始时刻，舱段试验的工况1出现最大刚体运动速度的时间大约为1.5 ms，而实船约为4 ms。

由舱段和实船的水下爆炸试验可知，不论是舰船舱段还是整船试验的刚性运动速度，以上理论分析得到的理论值与试验值的相对误差，都控制在10%以内，充分说明了本文理论计算方法在预报船体水下爆炸刚性运动的有效性。

3 刚体运动对舰船总体响应影响

上文实尺度舱段试验的模型取材于某舰船，拟通过局部结构试验来考核全船的局部强度，研究舰船在水下爆炸冲击环境中的响应。下面以此舰船为研究对象，分析在不同工况的水下爆炸环境中舰船刚体速度响应对舰船总体响应的影响，从而为局部结构试验设计提供理论指导。

3.1 刚体运动速度对舰船总体速度影响

舰船水下爆炸节点速度响应可认为是刚体运动速度和弹塑性变形响应速度的叠加，为把刚体运动速度从节点速度响应的结果分离出去，需要得到在这个过程中刚体运动速度。

现根据舰船各站的剖面惯性矩和分布质量，把整个舰船等效成变截面梁，船体共分20站，船体中部较重，两端质量相对较小，将水下爆炸的冲击波载荷根据不同分段受到不同的冲击波作用力分别加载在变截面梁对应节点上，这样分析的问题就成为变截面梁在冲击载荷作用下的动态响应求解问题。此时可用多体系统的传递矩阵法或有限元法进行分析，传递矩阵法编程较复杂，且求解效率上不如有限元方法，同时有限元在求解变截面梁的动响应问题已相当成熟，求解精度也较高，所以在此选用有限元方法对变截面梁在冲击载荷作用下的动态响应问题进行分析。在对实船进行考核时，采用1 t的TNT炸药，使实船的壳体冲击因子与舱段试验的3个工况对应相等，药包的位置设置在舰船试验舱段中部的正下方，相应的爆距设置如表3所示。

表3 舰船试验工况Table 3 Ship test conditions

在将爆炸冲击波压力转换为节点力输入时，不同位置力的作用时间不同，相对于船体中部正下方爆炸工况，船体中部先受到力的作用，之后随着冲击波的传播，两侧开始受力。上面已验证了理论方法在实船计算中的有效性，以工况Ⅲ为例，应用式(8)可编程得到船体受到的冲击作用力，如图8，第10站为船中位置，其相对于一侧的2站在受力时，约提前18 ms，而位于船长1/4位置的第5站，受力时间则在第10站和第2站之间，约比第10站延迟8 ms。从图中还可看出离药包越近的位置，受力越大，但作用时间也相对于爆距远的位置作用时间短。在船体受力达到峰值时，爆距越远，则力变化的越缓慢。

图8 沿船长不同位置受力随时间变化Fig.8 Variation curve of the different locations along the ship receiving the force with time

沿船长方向取不同节点，通过有限元计算得到的不同工况节点最大速度如图9所示。工况Ⅰ的速度响应最大，工况Ⅲ最小。3种工况下船体两端的速度响应较其他位置偏大，这与船体首尾重量较轻有关。同时船体中间部分沿船长有较小的起伏变化，但整体在某一稳定值附近变化。通过本文计算船体刚体运动速度的方法可得到舰船的刚体运动速度，用有限算法得到节点整体速度响应减去刚体运动速度则可认为是船体结构弹塑性变形的速度。选取船体中部节点，以工况Ⅰ为例，绘制出整体变形速度、刚体运动速度及弹塑性变形速度随时间的变化如图10所示，在冲击波阶段，3种速度响应都在8 ms左右出现最大值，之后由于刚体速度逐渐趋于零，整体速度和弹塑性变形速度变化曲线基本重合。

船体在3种工况下沿船长分布的各个节点最大速度响应的平均值如表4所示，从表中可以看出对于冲击因子最大的工况Ⅰ整体速度响应最大，刚体运动速度也最大，同时弹塑性变形也是最大的，根据刚体速度占整体速度比重，得出全船实际得到的节点速度响应，包含总速度1/4左右的刚体运动速度，如不考虑刚体运动速度的响应，会对整体结构动响应问题分析带来很大的误差。

图9 不同工况船体节点速度沿船长变化规律Fig.9 The change law of hull nodes'speed along the ship length in different cases

图10 船体中部节点速度响应Fig.10 Velocity response of the hull central node

表4 船体最大速度响应Tab le 4 M axim um velocity of the hull

3.2 刚体运动动能对舰船响应影响

通过上文分析，刚体运动速度在船体整体速度响应分析中不可忽略，由刚体运动产生的刚体运动动能会削弱冲击波入射能量作用给船体结构弹塑性变形的变形能，下面将讨论刚体运动动能对船体结构吸收冲击波入射能量的影响。

冲击波对结构的入射能量Es可通过结构遮挡的冲击波能量求得［12］，在无限域流场中可认为冲击波能量均布于整个球形冲击波波阵面，冲击波的入射能量为

式中:E为冲击波总能量，R为爆距，Se为结构在垂直于冲击波波阵面上的投影面积，可通过文献［12］中的方法进行求解，图11为数值模拟冲击波波阵面求解全船在波阵面上投影面积Se示意图。

图11 全船数值模拟Fig.11 The schematic diagram of ship numerical simulation

在一次爆炸中，冲击波携带的能量为

式中:Ef为能流密度，指冲击波在传播方向上经过某一固定的单位面积内的能流，其表达式［1］为

由式(15)、(16)即可求得结构受到的冲击波能量Es，通过式(12)可得到船体刚体运动最大速度，即可求得船体刚体运动消耗的能量Er，则船体弹塑性变形等吸收的变形能为(包括船体结构的弹塑性变形能和流体的耗散能，在此不研究流体耗散能量，将变形能定义为冲击波入射能减去刚体运动动能):

在对船体舱段进行试验考核全船的水下抗爆性能时，只是按照冲击因子相等来设置试验，而实际要达到的目的是要弹塑性变形的能量相同，舱段和实船实现相同的弹塑性变形响应，从而考核实船的抗爆性能。舱段试验的工况1、2、3分别对应全船水下爆炸工况的Ⅰ、Ⅱ、Ⅲ，他们的冲击因子分别相等。但根据舱段和药包的相对位置及舱段试验设计时大量的数值仿真计算结果，舱段试验的3个工况是按照要考核的药包在攻角30°时全船舱段对应位置响应进行设置的。工况Ⅳ、Ⅴ、Ⅵ分别对应实船工况Ⅰ、Ⅱ、Ⅲ，药量和爆距均对应相同，但设置攻角为30°，如图12 所示。

根据试验舱段在全船水下爆炸工况中的投影面积，可求出舱段的冲击波入射能;已知全船刚体运动速度，根据舱段占全船重量的比重，可得到舱段在全船工况中的刚体运动动能。以此得到的全船上局部舱段受到的结构变形能，大小应与舱段试验的变形能相等。在按照舱段和全船冲击因子相等设置工况时，舱段试验和全船不同工况的冲击波能量如表5～6所示。

图12 全船工况设置Fig.12 Test conditions of ship

表5 舱段试验能量计算Table 5 Energy calculation of cabin test

表6 舱段在全船工况下的能量计算Table 6 Energy calculation of cabin in ship conditions

从表5和表6中数据可以看出，无论是舱段试验还是全船水下爆炸工况，船体刚体运动的动能相对于冲击波入射能不能忽略，尤其舱段试验的刚体运动动能占冲击波入射能的比重更大，在工况2时高达24%。

最终计算出的舱段试验各工况的变形能与全船各个工况的变形能并不相等，尽管舱段和全船各工况的冲击因子对应相等。舱段试验工况1的变形能为1 230 kJ，而全船上在对应冲击因子相同的工况Ⅳ时对应舱段位置的变形能为1 720 kJ，比舱段变形能高出40%。工况Ⅴ实船上对应舱段位置变形能和舱段试验工况2的变形能比较接近，认为此试验舱段工况可以考核对应实船冲击因子的工况，而实船工况Ⅵ的对应舱段位置变形能比舱段试验工况3的变形能小，可见在攻角30°时，大冲击因子的实船工况，舱段试验需要通过适当加大冲击因子来达到实船相等的变形能，在冲击因子较小时舱段变形能可能大于实船的变形能。

通过上面分析，若要考核全船某一工况的爆炸冲击动响应问题，只是根据冲击因子相同考核局部舱段的动响应，并不能反映实船的真实响应。较为准确的方法是按照本文思路计算全船上对应舱段位置的变形能，根据此变形能和本文计算模型刚体运动速度的方法得出舱段刚体运动动能，得到试验的冲击波入射能，最终确定试验的药量和爆距。具体求解流程如图13所示。

图13 局部结构试验工况设置求解框图Fig.13 The solving block diagram of the test conditions set of local structure

通过此求解思路，计算出本文舱段试验的工况只能考核如表7所示实船攻角为30°的工况。

表7 舱段试验能考核的实船工况Table 7 Real ship conditions of the cabin test assessment

若要考核实船要求的工况Ⅳ、Ⅴ、Ⅵ冲击因子，则根据同样的算法，求出舱段试验需要设置的工况冲击因子如表8所示。

表8 要考核实船工况需设置的舱段试验工况Table 8 Cabin test conditions to be set for assessment of ship conditions

4 结论

本文将复杂船体外壳等效为刚性平板，研究了船体在水下爆炸环境中球面冲击波作用下的刚体动态响应特性，计算出船体受到的冲击波作用力及刚体运动速度，研究了船体整体速度响应和刚体运动速度之间及冲击波入射能、刚体运动动能之间的关系，得出如下结论:

1)船体在球面冲击波作用下，计及冲击波的反射及船体与流场流固耦合的附加质量，对于舱段模型试验及实船试验，当药包在水下任意位置时，应用本方法可得出与试验结果较符合的刚体运动速度;

2)在本文设定的实船工况中，船体整体速度响应、刚体运动速度及弹塑性变形速度的最大值基本出现在同一时刻，刚体速度响应占总速度响应的1/4左右，在水下爆炸冲击波作用下进行结构动响应问题分析时，需考虑结构刚体运动对其的影响;

3)舱段模型试验中刚体运动动能占冲击波入射能的比重大于全船中舱段刚体运动动能占冲击波入射能的比重，为使舱段模型试验局部结构吸收的变形能与全船对应位置相等，可通过改变冲击因子的方法来设置局部结构试验的工况，考核实船对应冲击因子的结构动态响应。

［1］COLE R H.Underwater explosion［M］.Princeton:Princeton University Press，1948:270-352.

［2］李良军，蒋伟康，艾艳辉，等.水下爆炸载荷下小型目标变形及冲击损伤试验［J］.上海交通大学学报，2010，44(10):1450-1455.LILiangjun，JIANG Weikang，AI Yanhui，et al.Experimental study on deformation and shock damage of small target subjected to underwater explosions［J］.Journal of Shanghai Jiaotong University，2010，44(10):1450-1455.

［3］LIJ，RONG JL.Experimental and numerical investigation of the dynamic response of structures subjected to underwater explosion［J］.European Journal of Mechanics，2012，32:59-69.

［4］JIN Q K，DING G Y.A finite element analysis of ship sections subjected to underwater explosion［J］.International Journal of Impact Engineering，2011，38(7):558-566.

［5］JEN C Y，TAI Y S.Deformation behavior of a stiffened panel subjected to underwater shock loading using the non-linear finite elementmethod［J］.Materials ＆ Design，2010，31(1):325-335.

［6］李海涛，朱锡，张振华.水下爆炸球面冲击波作用下船体梁的刚塑性动响应特性［J］.工程力学，2010，27(10):202-207.LIHaitao，ZHUXi，ZHANG Zhenhua.Dynamic rigid-plastic response of ship-like beam subjected to underwater spherical shockwaves［J］.Engineering Mechanics，2010，27(10):202-207.

［7］ZHANG N，ZONG Z.The effect of rigid-body motions on the whipping response of a ship hull subjected to an underwater bubble［J］.Journal of Fluids and Structures，2011，27:1326-1336.

［8］姚熊亮.舰船结构振动冲击与噪声［M］.北京:国防工业出版社，2007:144-146.

［9］汪玉，华宏星.舰船现代冲击理论及应用［M］.北京:科学出版社，2005:220-223.

［10］张阿漫，戴绍仕.流固耦合动力学［M］.北京:国防工业出版社，2011:35-47.

［11］张文鹏，宗智，汪玉，等.冲击波作用下舰船刚体运动响应［J］. 噪声与振动控制，2012，32(6):134-136.ZHANGWenpeng，ZONG Zhi，WANG Yu，et al.Motion response of rigid ship under shock wave［J］.Noise and Vibration Control，2012，32(6):134-136.

［12］曹宇.水下爆炸冲击环境相似性研究［D］.哈尔滨:哈尔滨工程大学，2006:27-30.CAO Yu.Research on environment comparability of underwater explosion［D］.Harbin:Harbin Engineering University，2006:27-30.