水下爆源位置诱导舰艇总强度鞭状响应规律

晏卫东，黄亚舒，刘智博，杨娜娜，王志凯

（1. 92941部队，辽宁葫芦岛， 125001；2.哈尔滨工程大学，黑龙江哈尔滨， 150001）

舰船的细长结构特征决定了其一阶固有频率与常见兵器爆炸后产生的气泡脉动频率相接近，因此，在气泡脉动载荷诱导下，舰船产生鞭状运动，将会给舰船的总纵强度带来严重的威胁。

舰艇的主要结构形式为由耐压外壳和加强筋组成的加筋壳。当受到武器攻击时，冲击波载荷和气泡脉动载荷先后会对舰艇结构造成破坏，使舰艇的生命力受到威胁，因此，加强对加筋壳结构毁伤特性研究的重要性不言而喻。

针对舰艇船体在水下爆炸作用下的鞭状响应，Keil A H指出了气泡脉动载荷对船体梁鞭状运动影响显著的主要原因，并基于球形气泡理论和船体梁理论，给出了结构在水下爆炸气泡载荷下动态响应的计算方法；Jeffrey W S对鞭状运动中的阻尼机制进行了讨论，并基于球形气泡理论和Morison公式，提出了鞭状运动的计算方法，结果显示，其方法与试验吻合良好；Zong Z对气泡载荷下自由梁的动态塑性响应进行研究；张振华使用有限元软件MSC.DYTRAN对冲击波载荷下加强壳动响应情况模拟仿真，其结果为环肋加强壳结构设计提供了参考；李玉节等研究了气泡脉动激起细长船模做鞭状响应的现象，得出了气泡作用下船体梁响应由刚体运动和弹性变形2 部分组成；张玮通过模态叠加法，计算了水下爆炸由气泡引起的舰艇鞭状运动规律；宫国田等采用有限元仿真分析软件ABAQUS 对水下爆炸载荷作用下加强壳的总体进行了模拟，为实际工程提供了有价值的参考资料；沈晓乐等采用试验的方法研究了在水下爆炸冲击载荷作用下的加肋加强壳结构模型的动态响应，证明了应变测量的方法对结构进行模态分析的可行性；刘云龙等采用双渐进法对水下爆炸冲击载荷作用下加强壳结构损伤问题进行了研究，发现爆距较近时，气泡脉动载荷占主要成分，较远时，冲击波占主要成分；王海坤等采用二阶DAA 流固解耦技术，揭示加强壳梁模型在水下爆炸载荷作用下的总体低频运动特性及鞭状响应规律，为研究细长体加强壳结构在水下爆炸载荷作用下的鞭状效应总体损伤，提供了试验基础及理论方法；黄晓明等利用试验的方法对加强壳在水下爆炸作用下的鞭状响应试验进行了研究；Zhang N在考虑气泡迁移、自由面效应和拖曳力的基础上，建立气泡动力学方程，并使用龙格库塔方法计算实际船模动态弹塑性响应。

由于采用有限元数值方法计算船体梁鞭状运动的计算成本通常较大，且试验多受现实条件制约，故需要大量工况计算的规律总结性工程应用并不适用。本文提出1 种在水下爆炸作用下，一维船体梁鞭状运动响应计算的解析方法，该解析法利用速度势函数和Morsion 相对速度公式计算气泡载荷水动力，利用Newmark 数值求解法求解水下爆炸气泡载荷作用下船体梁动态响应，建立了一维船体梁模型与气泡载荷之间的耦合计算模型，并对船体梁模型的鞭状响应进行了计算，计算效率大幅提高，可实现在舰艇设计初期总强度的快速校核。为此，在验证了该船体梁解析计算方法的正确性后，引入无量纲系数，研究了爆炸参数对鞭状运动响应的影响规律，旨在为舰艇结构总体损伤评估和结构设计研究提供参考。

1 鞭状运动响应解析法

1.1 气泡载荷水动力计算

Geers（2002）运用二阶DAA法对水下爆炸气泡以及气泡周围流体的相互作用进行研究，在线性可压缩流体介质的假设下，提出基于球形气泡的水下自由场动力学气泡模型，即Geers-Hunter气泡模型。

根据Geers-Hunter 气泡模型进行求解，可以得到半径及气泡上浮速度随时间变化的情况，由此可以求得气泡的速度势。水下爆炸气泡坐标系，如图1所示。

图1 水下爆炸气泡坐标系Fig.1 Coordinate system of underwater explosion bubble

速度势函数可以表示为：

式（1）中：为点源强度系数；为偶极强度系数；为流场中任一点到气泡中心径向距离；为任一点到偶极中心径向距离；和为点源及偶极径向与垂向之间的夹角。

全局使用笛卡尔坐标系，速度势函数分别关于、和方向求偏导数，可以求得气泡载荷作用下流场中任一点流体的速度。计算公式如下：

对于细长结构在流场中的振动问题，通常采用切片理论的思想，简化为二维问题分析。本文将流场假设为不可压缩，考虑流体惯性及拖曳力，使用Morison相对速度公式，计算作用于结构的水动力。Morison相对速度公式如下：

式（3）中：为作用于结构外表面单元的水动力；为流体密度；为表面单元的横截面积；为表面单元直径；为单元长度；和̇分别为流体质点在径向的速度和加速度；̇和̈分别为结构点在径向的速度和加速度；为流体惯性相对运动系数；为拖曳力相对运动系数，本文取1.0。利用Geers-Hunter气泡模型得到船体梁各位置的流体质点速度以及加速度，进而得到任意时刻船体梁所受的水下爆炸载荷。

1.2 船体梁响应求解

在由水下爆炸载荷引起的船体梁鞭状运动的研究中，通常将船体梁简化为包含集中质量的一维变截面梁。由于引起船体梁鞭状运动响应的主要原因是低阶频率共振，且一维梁模型在计算低阶频率响应方面较为准确，因此使用此模型具有合理性。

当受到外界激振力作用时，含有阻尼的多自由度振动系统会产生强迫振动，根据牛顿运动定律以及达朗贝尔原理建立的强迫振动微分方程如下：

式（4）中：为船体梁质量矩阵；为船体梁阻尼矩阵；为船体梁刚度矩阵；̈为加速度矢量；̇为速度矢量；为位移矢量；为激振力矢量。

在解析解法中，为了保证计算的准确性且能够降低计算难度，选用的阻尼类型为瑞利阻尼。瑞利阻尼矩阵可以表示为质量矩阵和刚度矩阵的线性叠加形式：

对于船体梁在气泡载荷阶段的动响应运动方程的求解，由于船体梁在水下爆炸载荷作用下的刚体位移成分占有一定比例，因此，运动方程更适合采用时域显式积分方法求解。本文采用Newmark 方法求解运动方程，Newmark方法的基本方程是对位移和速度进行Taylor 展开建立的，位移和速度的Taylor 展开式如下：

Newmark 指出，当取值0.25，取值0.5 时，求解过程无条件稳定。船体梁的运动方程可以写成如下形式：

根据上式可以求得，+时刻的平衡方程为：

Suzuki H指出，将Newmark 方法应用于非线性载荷作用下结构动力学方程求解问题时,建议采用Newmark-Raphson 迭代方法，故本文采用此迭代法对船体梁在气泡载荷阶段的动响应运动方程进行求解，得到任意时刻船体梁的位移、速度和加速度响应。

2 解析法验证

本研究中模型参数为：加强壳总长度6 m，壳体厚度2 mm，壳体总质量51.7 kg，总配重76.6 kg，总排水量128.3 kg。以某型号舰艇横舱壁位置为参考，确定加强筋位置。试验模型具体尺寸，如图2所示。

图2 模型尺寸参数及测点布置图Fig.2 Model size parameters and measuring point layout

爆心位置为加筋加强壳中部正下方，即方位角为90°。药量的选取原则为由炸药引起的气泡载荷的脉动周期与模型低阶固有周期接近。最终确定的试验工况，如表1所示。

表1 试验工况设置Tab.1 Setting of test conditions

在试验研究的基础之上建立加筋加强壳数值模型，基于声固耦合方法，对加筋加强壳鞭状运动响应进行数值仿真计算，工况设置与试验相同。模型使用的材料是Q235 钢，密度为7 850 kg/m，弹性模量为2.1×10Pa ，钢 的 泊 松 比 为0.3；水 的 密 度 为1 000 kg/m，水的体积模量为2.140 4×10Pa 。钢材本构模型选用Plastic-Kinematic 模型。加筋加强壳采用S4R四节点壳单元建立，外流场采用AC3D4四节点声学介质单元建立。加筋加强壳及外流场的三维有限元模型如图3、4所示。

图3 加筋加强壳有限元模型Fig.3 Finite element model of the reinforced cylindrical shell

图4 外部流场有限元模型Fig.4 External flow field finite element model

使用有限元仿真软件得到的前三阶湿模态固有频率数值解的相对误差分别为5.12%、1.99%和0.90%，可以认为，所建有限元模型的总体刚度分布、质量分布、外流场分布与真实模型较为接近。

为了验证解析法计算结果的正确性，将工况1、2的船体梁中部截面应变的数值解、试验解和解析解绘制在同一张图中，并从应变峰值和峰值出现时间2 方面进行分析。

图5 中部截面应变时历曲线数值解、试验解和解析解对比Fig.5 Comparison of numerical solution,experimental solution and analytical solution of the time-history curve of strain in the middle section

表2 中部截面应变最大值数值解、试验解和解析解对比Tab.2 Comparison of numerical solution,experimental solution and analytical solution for the maximum strain in the middle section

从表3 可以看出，数值解和解析解的应变幅值和首次气泡脉动产生鞭状运动峰值时间都在一定的误差允许范围之内。由此验证了鞭状运动数值仿真解法和解析解法的准确性。在鞭状运动应变幅值与首次气泡脉动产生应变峰值时间预报方面，解析解法和数值解法准确性相差不大，在计算成本方面，使用解析解法预报鞭状运动无须建立模型，计算时间也远远小于数值仿真解法，但解析解法对输入数据的准确性要求较高。

表3 首次鞭状运动峰值时间数值解、试验解和解析解对比Tab.3 Comparison of numerical solution,experimental solution and analytical solution of peak time of the first whip motion

3 船体梁的鞭状响应

基于解析法，计算该船体梁模型的鞭状运动响应。爆点与加筋加强壳的轴线保持在同一水平面，且位于加强壳长度方向的中部（1/2处，为加强壳总长度）。采用加强壳中部截面应变幅值来评估炸药深度对鞭状运动响应的影响，分别讨论了炸药深度、爆距和爆点轴向位置对鞭状运动响应的影响规律。

3.1 炸药深度对鞭状运动响应影响

炸药深度影响了气泡脉动周期，进而影响了鞭状运动周期和应变幅值，为了减少这部分误差，得到适用于不同模型的普遍规律，因此引入无量纲系数：

式（11）中：为气泡脉动周期；为结构一阶湿模态固有周期。

图6 显示了不同炸药深度条件下，加筋加强壳鞭状运动对比情况。从图中可以明显看出，随着炸药深度不断增加，加筋加强壳鞭状运动的周期缩短，气泡载荷脉动周期逐渐缩短。因此可以得出结论：气泡载荷脉动周期同鞭状运动周期呈正相关。

图6 不同炸药深度下加筋加强壳鞭状运动应变时历曲线Fig.6 Time-history curves of whip motion strain of the reinforced cylindrical shell at different explosive depths

具体工况设置，如表4 所示。表5 显示了不同工况下鞭状运动最大应变幅值同气泡脉动周期和加强壳一阶湿模态周期之比之间的关系。

表4 工况设置Tab.4 Setting of working conditions

表5 不同工况下鞭状运动最大应变幅值Tab.5 Maximum strain amplitude of whip motion under different working conditions

为了能够更加直观地分析规律，使用3 次样条插值的方法，将数据绘制成图7。

图7 鞭状运动最大应变幅值随δ 变化曲线Fig.7 Curve of maximum strain amplitude of whip motion with δ

从图7可以看出，随着周期比增加，加筋加强壳鞭状运动的应变幅值逐渐增加。当为1.34时，鞭状运动的应变幅值达到最大值，随后逐渐衰减。

综合不同工况分析可知，不同炸药深度对应工况的气泡脉动周期不同，气泡脉动时对应加强壳鞭状运动的波峰、波谷情况不同。当089 ≤≤134 时，多个气泡脉动的峰值正好作用于鞭状运动的上升阶段，气泡膨胀对结构做功为正，因此，对应的鞭状运动响应的幅值较大；当134 ≤≤145 时，多个气泡脉动的峰值处于结构鞭状运动的下落阶段，气泡压力方向与加强壳运动方向相反，气泡膨胀对结构做功为负，因此，加筋加强壳本身动能减少，鞭状运动响应变小。

3.2 爆距对鞭状运动响应影响

除了炸药深度会对鞭状运动产生影响外，爆距对于鞭状运动的幅值和变化趋势也会产生显著影响。不同长度的模型对于爆距的敏感程度不同，为了减少这部分误差，以得到适用于不同模型的普遍规律，因此引入无量纲系数：

式（12）中：为爆距；为结构特征长度。

具体工况设置，如表6所示。

表6 工况设置Tab.6 Setting of working conditions

不同爆距下，加筋加强壳鞭状运动响应应变幅值情况，如图8所示。

图8 不同爆距下加筋加强壳鞭状运动应变时历曲线Fig.8 Time-history curves of whip motion strain of the reinforced cylindrical shell at different detonation distances

从图中可以看出，爆距对于鞭状运动的周期影响较小，但对鞭状运动幅值影响却很大。随着爆距增大，鞭状运动的应变幅值逐渐减小，爆距和鞭状运动应变幅值呈负相关。

表7显示了不同工况下鞭状运动最大应变幅值同爆距和加强壳特征长度比之间的关系。

表7 不同工况下鞭状运动最大应变幅值Tab.7 Maximum strain amplitude of whip motion under different working conditions

为了能够更加直观地分析规律，使用3 次样条插值的方法，将数据绘制成图9。图中显示了炸药深度为15 m时，鞭状运动最大应变幅值随长度比的变化趋势。从图可以看出：当05 ≤≤07 时，应变幅值随长度比变化速率较快；当07 ≤≤167 时，应变幅值变化速率逐渐减慢。从整体上看，应变幅值和长度比呈反比例函数关系，其原因是：随着爆距增加，经过流场作用于结构体上的压力幅值逐渐减小；结构体特征长度相对于爆距比例逐渐减小，作用于结构体上的压力载荷沿着长度方向分布逐渐均匀，结构体响应中刚体位移所占比例逐渐升高，而鞭状运动就是由于作用在结构体上的载荷不均而造成的。最终得出：结构体鞭状运动应变幅值随长度比升高而降低。

图9 鞭状运动最大应变幅值随η 变化曲线Fig.9 Curve of maximum strain amplitude of whipping motion with η

3.3 爆点轴向位置对鞭状运动响应影响

爆点沿船长方向的位置也会对船体梁的响应产生影响，仍以试验加筋加强壳模型为考察对象，分析爆点位置偏移对鞭状运动响应的影响规律。同理设置无量纲系数：

式（13）中：为爆点的轴向位置；为结构特征长度。

由于加筋加强壳为对称结构，因此，只选取模型的一半进行考察。具体工况设置，如表8所示。

表8 工况设置Tab.8 Setting of working conditions

每种工况下，分别对=0 m、=1.5 m、=3 m、=4.5 m 和=6 m 的截面应变状况进行考察。

由于加强壳为对称结构，爆点轴向位置为0 m、1.5 m 时典型截面应变时历曲线与4.5 m、6 m 时应变时历曲线相同，因此只展示工况1、2 和3 的应变时历曲线。从图10可以看出，3种工况下,爆点轴向位置改变，但对于每个工况而言，船舯截面的响应应变幅值仍然要比两侧截面的应变幅值大。爆点轴向位置为3 m 时，不同截面的应变时历曲线波形较为简单；爆点轴向位置为0 m和1.5 m时，不同截面的应变时历曲线波形比较复杂，在原有基础波形上有其他小波叠加。

图10 不同轴向位置爆点下典型截面应变时历曲线Fig.10 Time-history curves of typical cross-sectional strain at different axial positons of the burst point

分别对不同工况下不同截面的应变时历曲线进行傅里叶变换，得到不同工况下的典型截面应变响应傅里叶谱,如图11所示。

图11 不同轴向位置爆点下典型截面应变响应傅里叶谱Fig.11 Fourier spectrum of typical cross-sectional strain response at different axial positions of the burst point

使用工程计算方法得到的加筋加强壳前三阶湿模态频率与图中傅里叶谱峰值对应频率相吻合。从图11 可以看出：爆点轴向位置3 m 时，爆炸载荷主要引起的是结构一阶模态共振，对应的应变曲线波形较为简单，只由一阶模态响应组成；爆点位置偏移至轴向1.5 m 时，加强壳的二阶、三阶模态比例升高，一阶模态响应比例下降，但对于总响应而言，仍占很大比重，对应的应变时历曲线由一阶、二阶和三阶模态响应共同组成，因此波形较为复杂；当爆点轴向位置为0 m 时，总体响应相对于爆点轴向位置1.5 m 时的二阶、三阶模态响应成分下降，一阶模态响应成分上升。

为探究爆点轴向位置的改变对加强壳鞭状运动响应的影响规律，现将考察截面位置确定，只改变爆点轴向位置，对比分析同一截面不同工况下鞭状运动响应情况。不同工况下考察截面的最大应变幅值，如表9所示。

表9显示了不同工况下鞭状运动最大应变幅值同爆点轴向位置和加强壳特征长度比之间的关系。为了能够更加直观地分析规律，使用3 次样条插值的方法，将数据绘制成图12，由于模型对称，因此只分析考察截面位置为0 m、1.5 m和3 m时的情况。

表9 不同工况下典型截面最大应变幅值Tab.9 Maximum strain amplitude of typical sections under different working conditions 10-6

图12 典型截面下鞭状运动最大应变幅值随λ 变化曲线Fig.12 Curve of maximum strain amplitude of whip motion under typical sections with λ

图12显示了不同考察截面位置下，船体梁鞭状运动响应应变幅值同爆点轴向位置和加强壳特征长度比之间的关系。总体上看，在不同考察截面位置条件下，船体梁鞭状运动应变幅值都存在随着先增大后减小的趋势，当=0.5 时，应变幅值达到峰值，因此，对于均匀对称分布的船体梁模型，爆点轴向位置在船舯时鞭状运动最为剧烈。

当考察截面位置=0 m 时，0 ≤≤0.5 部分比0.5 ≤≤1部分应变幅值要大，考察截面位置=1.5 m时也有类似规律，原因是0 ≤≤0.5 时考察截面与爆点位置之间距离比0.5 ≤≤1 时要大，鞭状运动响应随着爆距增大而减小，这与3.2 节得到的爆距对鞭状运动响应影响规律相吻合。

4 结论

本文基于水下爆炸载荷下加强壳结构总纵强度预报的工程解析解法，以试验模型为研究对象，计算试验工况下加强壳鞭状运动应变响应。并将通过试验、解析解法和数值仿真解法得到的应变响应曲线从最大应变幅值和首次鞭状运动峰值时间2方面进行对比分析，验证数值解法在预报鞭状运动应变响应方面的有效性。使用解析解法分别讨论了炸药深度、爆距和爆点轴向位置对鞭状运动响应的影响规律，并对典型工况下船体梁振型图进行分析，得到主要结论如下：

1）在鞭状运动应变幅值与首次气泡脉动产生应变峰值时间预报方面，解析解法和数值解法准确性相差不大，均能很好地计算船体梁在水下爆炸载荷下的响应；

2）结合文中所示不同工况条件下的模拟计算结果，当不考虑自由面效应时，船体梁结构鞭状运动应变幅值随着周期比先增加后减小，当周期比取值为1.34时，应变幅值达到峰值；

3）结合文中所示不同工况条件下的模拟计算结果，当不考虑自由面效应时，鞭状运动应变幅值和长度比呈反比例函数关系；

4）结合文中所示不同工况条件下的模拟计算结果，当不考虑自由面效应时，船体梁鞭状运动响应应变幅值随着长度比先增大后减小，当长度比取值为0.5时，应变幅值达到峰值。