0 引　言

舰船结构在水下非接触爆炸下除可能产生局部破坏外，还可能发生总纵弯曲破坏。当水下爆炸气泡脉动周期与舰船结构的低阶固有频率接近或相当时，舰船将出现鞭状运动，严重时将导致舰船折断而沉没。国内外学者采用理论分析、模型试验和数值计算等方法对舰船结构的鞭状运动进行了研究^[1-9]，但基于Ls-dyna软件对水下非接触爆炸下舰船结构鞭状运动进行研究的文献并不多见。鉴于此，本文基于Ls-dyna软件的ALE算法，对排水量5890 t的某船体梁在315 kg TNT水下近距爆炸下的鞭状运动进行数值模拟，分析药深对船体梁鞭状运动的影响，并探究船体梁鞭状运动与水下爆炸气泡脉动过程在时间上的同步性，为舰船抗爆设计提供参考。

1 有限元计算模型 1.1 有限元建模

考虑对称性，建立如图1所示的1/2有限元模型。该模型由船体梁和欧拉域两大部分组成，欧拉域包括空气域、水域、TNT、水边界层和空气边界层。

图1所示1/2船体梁为等截面梁，长150 m，宽5 m，高10 m，吃水5 m，排水量2945 t；水线以上为方形截面，水线以下为1/4圆形截面。船体梁自上而下以2.5 m为间距设置4层甲板，沿梁长方向以10 m为间距设置14块横隔板，两端各设置1块封闭板。船体梁各板厚度统一设为16 mm。船体梁采用拉格朗日壳单元，网格划分如图2所示，单元尺度约为0.4～0.6 m，共15786个壳单元。为便于研究水下近距爆炸下船体梁的鞭状运动，选取中纵剖面、顶层甲板和横隔板及封闭板的交汇点作为测点，沿X正向依次编号为节点1～16，计算输出这16个节点的Z向位移时程曲线。

图1所示1/2欧拉域长180 m，宽30 m，高67 m。其中，空气域高30 m，水域高37 m。欧拉域网格划分如图3所示，主要为六面体单元，单元尺度约为0.5～1.5 m，共768496个体单元。将欧拉域边界的1层单元（不包括在1/2模型对称面上而实际在全模型内部的单元）作为边界层，用来定义水、空气等材料的流入流出边界。

315 kg TNT炸药分别在船体梁正中部水下7 m，9 m，11 m爆炸，依次编号为工况1、工况2、工况3。在1/2模型中，炸药质量为157.5 kg。分别选取船体梁正中部水下7 m，9 m，11 m处的1个单元（其体积为0.125 m³）作为炸药，炸药密度为1260 kg/m³，炸药单位体积内能为5.411×10⁹ J/m³。

1.2 材料模型及参数

空气的密度为1.28 kg/m³，采用Linear-Polynomial状态方程，即

$ p = {C_0} + {C_1}\mu + {C_2}{\mu ^2} + {C_3}{\mu ^3} + ({C_4} + {C_5}\mu + {C_6}{\mu ^2}){E_0} 。$

(1)

式中：p为压力；μ = 1/V – 1；E₀为单位体积内能；C₀，C₁，C₂，C₃，C₄，C₅，C₆为常数，各参数取值^[10]为：C₄ = C₅ = 0.4，E₀ = 2.525×10⁵ J/m³，V₀ = 1，其余为0。

水的密度为1000 kg/m³，采用Gruneisen状态方程，即

$ p = \dfrac{{{\rho _0}{c^2}\mu \left[ {1 + \left( {1 - \dfrac{{{\gamma _0}}}{2}} \right)\mu - \dfrac{\alpha }{2}{\mu ^2}} \right]}}{{{{\left[ {1 - \left( {{S_1} - 1} \right)\mu - {S_2}\dfrac{{{\mu ^2}}}{{\mu + 1}} - {S_3}\dfrac{{{\mu ^3}}}{{{{\left( {\mu + 1} \right)}^2}}}} \right]}^2}}} + \left( {{\gamma _0} + \alpha \mu } \right){E_0}。$

(2)

式中：p为压力；c为声速；μ = 1/V – 1；E₀为单位体积内能；S₁，S₂，S₃，γ₀，α为常数，各参数取值^[8]为：c = 1484 m/s，S₁ = 1.979，γ₀ = 0.11，α = 3，V₀ = 0.999954，其余为0。

炸药的密度为1260 kg/m³，CJ压力为21 GPa，爆速为6930 m/s；采用JWL状态方程，即

$ p = A\left( {1 - \frac{\omega }{{{R_1}V}}} \right){e^{ - {R_1}V}} + B\left( {1 - \frac{\omega }{{{R_2}V}}} \right){e^{ - {R_2}V}} + \frac{{\omega {E_0}}}{V}。$

(3)

式中：p为压力；V为相对体积；E₀为单位体积内能；A，B，R₁，R₂，ω为特征参数，各参数取值^[10]为：A = 371.2 GPa，B = 3.231 GPa，R₁ = 4.15，R₂ = 0.95，ω= 0.3，E₀ = 5.411×10⁹ J/m³，V₀ = 1。

船体梁材料为Q235钢，密度为7800 kg/m³，采用双线性弹塑性本构模型和最大塑性应变失效模型，用关键字*MAT_PLASTIC_KINEMATIC定义；应变率效应采用Cowper-Symonds模型描述，动态屈服强度σ_d为

$ {\sigma _d} = \left( {{\sigma _0} + \frac{{E{E_h}{\varepsilon _p}}}{{E - {E_h}}}} \right)\left[ {1 + {{\left( {\frac{{\dot \varepsilon }}{C}} \right)}^{\frac{1}{n}}}} \right]。$

(4)

式中：σ₀为静态屈服强度；E为杨氏模量；E_h为硬化模量；ε_p为等效塑性应变； $ \dot \varepsilon $ 为等效塑性应变率；C，n为常数；各材料参数取值^[10]为：E = 210 GPa，泊松比为0.3，σ₀ = 235 MPa，E_h =250 MPa，C = 40.4 s⁻¹，n = 5，失效应变为0.28。

由于建模时未考虑加强筋和设备，全部钢板密度取为7800 kg/m³时，船体梁总质量是小于排水量的。考虑到水下非接触爆炸下船体梁最上2层甲板之间的横隔板及封闭板的变形可以忽略不计，通过增大此部分板的密度（实质是将未考虑的加强筋和设备质量附加到此部分板上）以使船体梁总质量等于排水量5890 t。

1.3 重力场与边界层

为正确模拟水下爆炸气泡脉动和船体梁鞭状运动的过程，必须考虑船体梁和流体受到的重力作用，因此用关键字*LOAD_BODY_PARTS和*LOAD_BODY_Z定义重力场。

空气域和水域的边界层定义为环境单元，使其内部压力在数值计算过程中保持为初始静压力；定义流入流出边界，使空气和水等流体材料在边界上能够流入流出。

1.4 流固耦合与接触算法

采用关键字*CONSTRANED_LAGRANGE_IN_SOLID定义ALE（任意拉格朗日欧拉）算法，以计算空气、水及TNT炸药和船体梁结构之间的流固耦合作用。采用关键字*CONTACT_AUTOMATIC_SINGLE_SURFACE定义自动单面接触算法，以计算船体梁结构中各钢板之间的接触作用。

2 有限元计算结果及分析 2.1 水下爆炸气泡

为便于与文献[11]对比分析水下爆炸气泡脉动过程，忽略船体梁对水下爆炸气泡脉动的影响。在自由场中，315 kg TNT在水下9 m爆炸时的气泡脉动过程如图4所示。在约0.48 s时，气泡第1次膨胀至最大，呈梨形，水平半径约为8.7 m，垂向半径约为9.6 m；在约0.96 s时，气泡第1次收缩至最小，呈环形，水面形成高达11 m的射流，射流顶部速度由0.48 s时刻的9.85 m/s逐渐减至0.96 s时刻的4.81 m/s；在约1.40 s时，气泡第2次膨胀至最大，射流呈皇冠型；在约1.80 s时，气泡第2次收缩至最小，并基本上浮至水面。随后，气泡逐渐逸出水面并溃灭。

315 kg TNT在水下9 m爆炸时药包上方7 m处的压力时程曲线如图5所示。可见，气泡在约0～0.96 s时间段内完成第1次脉动过程，在约0.96～1.80 s时间段内完成第2次脉动过程，在约1.80～2.40 s时间段内完成第3次脉动过程。

根据水下爆炸基本理论^[11]，无限水域中水下爆炸气泡的最大半径和第一次脉动周期可分别由式(5)和式(6)进行估算：

$ {R_{\max }} = 3.5{\left( {\frac{w}{{s + 10.3}}} \right)^{1/3}}，$

(5)

$ {T_1} = 2.11\frac{{{w^{1/3}}}}{{{{\left( {s + 10.3} \right)}^{5/6}}}} 。$

(6)

式中：R_max为气泡最大半径，m；T₁为气泡第1次脉动周期，s；w为TNT炸药质量，kg；s为药深，m。

由式(5)和式(6)计算得到，315 kg TNT在水下9 m爆炸时，气泡的最大半径理论值为8.9 m，第1次脉动周期理论值为1.22 s。气泡的最大半径与第1次脉动周期的仿真值和理论值比较吻合，因此本文采用的数值模拟方法可行。

2.2 船体梁的鞭状运动

仿真计算输出船体梁上测点1～16的Z向位移时程曲线，以便分析水下近距爆炸下船体梁的鞭状运动，以工况2为例进行分析。

工况2下节点1、节点8及二者相对的Z向位移时程曲线如图6所示。可见，船体梁在约0～0.4 s时间段内第1次中拱运动，在约0.4～1.0 s时间段内第1次中垂运动，在约1.0～2.0 s时间段内第2次中拱运动，在2.0 s后第2次中垂运动，到3.0 s时仍未停止中垂运动。船体梁的第1次中拱、中垂运动与炸药的第1次膨胀、收缩过程在时间上基本同步，而船体梁的后续中拱、中垂运动与炸药的后续膨胀、收缩过程在时间上稍有延迟。节点8相对节点1的Z向位移在0.386 s时达到最大值1.11 m，在0.970 s时达到最小值−0.60 m，在2.014 s时达到极大值0.67 m。

工况2下，水下爆炸气泡与船体梁的相互作用如图7所示，船体梁的最大中拱和最大中垂如图8所示。可见：在0.48 s时，气泡已经自顶部塌陷收缩；船体梁呈中拱状态，底板中部10 m范围内向内凹陷变形达2.6 m。在0.96 s时，气泡几乎收缩至最小，呈环状；船体梁呈中垂状态，底板凹陷变形范围及大小增大。在2.00 s时，船体梁呈中拱状态，底板中部30 m范围内向内凹陷变形达5.3 m。

2.3 药深对船体梁鞭状运动的影响

不同工况下节点8相对节点1的Z向位移时程曲线如图9所示。可见：不同药深下，船体梁的第1次中拱、中垂运动具有较好的时间同步性；随着药深增大，船体梁的第2次中拱、中垂运动逐渐增大，并且在时间上逐渐延迟。

不同工况下船体梁的最大中拱（t ≈ 0.38 s）和最大中垂（t ≈ 0.96 s）如图10所示。可见：工况1和工况2下船体梁的最大中拱程度相当，工况3下船体梁的最大中拱程度相比工况1、工况2是偏小的；随着药深增大，船体梁的最大中垂程度增大。

3 结　语

本文基于Ls-dyna的ALE算法，对排水量为5890 t的船体梁在315 kg TNT水下近距爆炸下的鞭状运动进行了数值模拟研究。水下爆炸气泡的最大半径与第一次脉动周期的仿真值和理论值比较吻合，因此本文采用的数值模拟方法可行。研究表明：船体梁的第1次中拱、中垂运动与炸药的第1次膨胀、收缩过程在时间上基本同步，而船体梁的后续中拱、中垂运动与炸药的后续膨胀、收缩过程在时间上稍有延迟；随着药深增大，船体梁的第2次中拱、中垂运动逐渐增大，并在时间上逐渐延迟。