0 引　言

为了防御鱼雷接触爆炸对大型舰船内部重要舱室造成破坏，在大型舰船水下舷侧设有多舱防护结构，其中防护液舱是防御鱼雷爆炸形成的高速破片的重要结构。大型舰船水下舷侧防护结构遭受鱼雷接触爆炸后，在舷侧外板上首先形成直径与鱼雷直径相当的剪切冲塞破片，此破片往舷侧空舱内部高速运动，直至撞击液舱前板。液舱前板在高速破片的撞击下发生剪切冲塞破坏，舷侧外板和液舱前板的冲塞破片以相同入水速度在液舱中运动，在水中形成冲击波，液舱后板在冲击波作用下发生大变形^[1]。如果防护液舱结构设计不合理，那么液舱后板轻则变形过大，重则被破片击穿。为了更好地设计防护液舱结构，弄清楚高速破片侵彻防护液舱过程中液舱后板的载荷特性具有重要意义。

国内外研究人员对不同工程背景下弹体或高速破片与液舱的相互作用问题进行了系列研究。国外研究主要针对弹体与航天器燃料舱的相互作用^{[2 – 6]}；国内研究主要针对高速破片与大型舰船防护液舱的相互作用^{[7 – 11]}，侧重于高速破片在液体中运动时产生的冲击波、压力、气穴、破片速度衰减特性及防护液舱的变形破坏等问题。尽管如此，高速破片侵彻防护液舱过程中液舱后板的载荷特性尚未研究清楚，因此本文对此进行数值仿真研究，旨在为更好地设计防护液舱结构提供参考。

1 有限元数值计算 1.1 有限元模型

为了研究高速破片侵彻防护液舱过程中液舱后板的载荷特性，利用LS_DYNA软件，建立如图1所示的有限元模型。

防护液舱的尺寸为18 m×1.8 m×9 m，采用Lagrange壳单元；欧拉域的尺寸为19 m×4 m×10 m，采用Euler单元。有限元网格的划分在中间较密、在四周较疏，Lagrange壳单元的边长约为60～300 mm，Euler单元的边长约为100～300 mm。液舱前板厚14 mm，液舱后板厚36 mm，液舱上下左右的板厚12 mm。防护液舱内为80%的水和20%的空气，防护液舱外为空气。欧拉域边界是1层厚0.1 m的环境单元，其内部压力恒为1.01×10⁵ Pa。

考虑到液舱前板在产自舷侧外板的高速破片撞击下的剪切冲塞过程较难模拟，为便于研究，不考虑液舱前板的剪切冲塞过程，而是直接在液舱前板中央开挖一个直径为533 mm（MK48鱼雷直径）的圆孔，并在此处放置一个直径同为533 mm的圆形破片，此破片的厚度是舷侧外板和液舱前板的厚度之和，在不同工况下破片的速度和厚度见表1。

采用ALE算法模拟空气、水和破片及防护液舱之间的流固耦合作用，采用自动单面接触算法模拟高速破片与液舱后板之间可能的接触作用。

1.2 材料模型及参数

此有限元模型涉及空气、水和钢这3种材料。

空气的初始密度为ρ₀=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\text{。}$

(1)

式中：p为压力；μ=ρ/ρ₀-1，ρ为密度；E为单位体积内能；C₀，C₁，C₂，C₃，C₄，C₅，C₆均为常数。各参数取值^[12]为：C₁=C₂=C₃=C₆=0，C₄=C₅=0.4，E=0.252 5 MPa。

水的初始密度ρ₀=1 000 kg/m³，采用Gruneisen状态方程，即

$p \!=\! \frac{{{\rho _0}{C^2}\mu \left[ {1 \!+\! \left( {1 \!-\! \frac{{{\gamma _0}}}{2}} \right)\mu \!-\! \frac{\alpha }{2}{\mu ^2}} \right]}}{{{{\left[ {1 \!-\! \left( {{S_1} \!-\! 1} \right)\mu \!-\! {S_2}\frac{{{\mu ^2}}}{{\mu \!+\! 1}} \!-\! {S_3}\frac{{{\mu ^3}}}{{{{\left( {\mu \!+\! 1} \right)}^2}}}} \right]}^2}}} \!+\! \left( {{\gamma _0} \!+\! \alpha \mu } \right)E\text{。}$

(2)

式中：p为压力；C为声速；μ=ρ/ρ₀-1，ρ为密度；E为单位体积内能；S₁，S₂，S₃，γ₀，α均为常数。各参数取值^[12]为：C=1 484 m/s，S₁=1.979，S₂=S₃=0，γ₀=0.11，α=3，E=0。

高速破片和防护液舱采用钢材，密度为7 800 kg/m³，采用双线性弹塑性本构模型和最大塑性应变失效模型，材料的应变率效应由Cowper-Symonds模型描述，动态屈服强度

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

(3)

式中：σ₀为屈服应力；E为杨氏模量；E_h为硬化模量；ε_p为等效塑性应变； $\dot \varepsilon $ 为等效塑性应变率；D和n为常数。各参数取值^[13]为：σ₀=685 MPa，E=210 GPa，E_h=1.218 GPa，D=8 000 s^–1，n=8，泊松比为0.3，最大塑性应变为0.3。

1.3 压力测点及计算工况

在防护液舱后板上选取如图2所示的25个典型单元作为压力测点，以研究防护液舱后板的载荷传播及衰减特性。

计算工况如表1所示。工况1～工况5的破片厚度δ=30 mm，破片速度v在1 000～1 400 m/s范围内变化，以研究破片速度v对防护液舱后板载荷特性的影响。工况1及工况6～工况9的破片速度v=1 200 m/s，破片厚度δ在26～34 mm范围内变化，以研究破片厚度δ对防护液舱后板载荷特性的影响。

2 计算结果及分析 2.1 工况1下防护液舱后板的载荷分析

在工况1下，测点A，J～Q（沿液舱后板的长对称轴分布）的压力时程曲线如图3所示，各测点的压力峰值p_max和比冲量i与测点到液舱后板中心距离r的关系分别如图4和图5所示。

由图3可见，在工况1下，测点A的压力峰值最大，约为61 MPa；随着测点到液舱后板中心的距离r增加，测点的压力峰值逐渐衰减。在以液舱后板中心为起点的其他射线方向上，液舱后板的载荷呈现出类似图3的传播和衰减特性，到液舱后板中心距离相同的几个测点的压力曲线基本重合。由图4和图5可见，在工况1下，各测点的压力峰值p_max和比冲量i随测点到液舱后板中心距离r变化的趋势线呈指数函数形式，经拟合，p_max=66.9E^–0.776r MPa，i=17 775E^–0.641r Pa·s，其中r以m为单位。

2.2 破片速度和厚度对防护液舱后板载荷的影响

各测点的压力峰值p_max和比冲量i随测点到液舱后板中心距离r变化的趋势线可取为p_max=K₁·e^αr和i=K₂·e^βr形式，其中，p_max，i和r分别以MPa，Pa·s和m为单位。

在各工况下，K₁和α值如表2所示，K₂和β值如表3所示。由表2和表3可见，破片厚度δ和速度v对K₁和K₂有显著影响，而对α和β的影响可以忽略不计，α和β的平均值分别为–0.767和–0.641。

由图6和图8可见，当破片厚度δ不变时，K₁和K₂均随着破片速度v增加而增加。由图7和图9可见，当破片速度v不变时，K₁和K₂均随着破片厚度δ增加而增加。当δ=30 mm时，K₁和K₂与v（以m/s为单位）的拟合关系式分别为K₁=0.058 8v^{0.992 7}和K₂=33.467v^{0.884 7}。当v=1 200 m/s时，K₁和K₂与δ（以mm为单位）的拟合关系式分别为K₁=6.823 8δ^{0.671 4}和K₂=2 167.8δ^{0.618 5}。

2.3 防护液舱后板载荷拟合公式及验证

基于上述分析，可以将液舱后板上任一点处的压力峰值p_max（以MPa为单位）和比冲量i（以Pa·s为单位）拟合成式（4）和式（5），即

${p_{\max }} = 0.006 \cdot {v^{0.9927}} \cdot {\delta ^{0.6714}} \cdot {e^{ - 0.767r}}\text{，}$

(4)

$i = 4.088 \cdot {v^{0.8847}} \cdot {\delta ^{0.6185}} \cdot {e^{ - 0.641r}}\text{。}$

(5)

式中：v为破片速度，m/s；δ为破片厚度，mm；r为液舱后板上任一点到其中心点的距离，m。

当v=1 050 m/s，δ=27 mm时，液舱后板上各测点的压力峰值p_max和比冲量i的有限元计算值和式（4）、式（5）计算值分别如图10和图11所示，两者吻合较好，从而对式（4）、式（5）进行验证。由式（4）、式（5）可看出，破片的速度v比厚度δ对液舱后板载荷的影响更显著。由于式（4）、式（5）是基于有限元计算结果得出的，并且未考虑液舱的宽度、液舱前后板的厚度及装载水的水位等因素对液舱后板载荷的影响，因此在工程应用上式（4）、式（5）有一定的局限性。尽管如此，对文中的防护液舱结构，当破片的速度v和厚度δ在一定范围内变化时，可以利用式（4）、式（5）对其液舱后板的载荷进行计算。

3 结　语

利用LS_DYNA软件，本文对不同速度和厚度的高速破片侵彻防护液舱过程中液舱后板的载荷特性进行数值仿真研究。基于仿真计算结果，首先分析了液舱后板载荷的空间分布特性，然后分析了破片速度和厚度对液舱后板载荷的影响，最后得到了液舱后板上任一点的压力峰值和比冲量的拟合计算公式。主要研究结果如下：

1）在空间分布上，液舱后板的载荷在板中心（即破片中心在液舱后板上的投影点）最大，随着到中心点的距离增加而呈指数衰减。

2）增加破片的速度或厚度（或者说增加破片的初始动量），将使液舱后板上任一点的压力峰值和比冲量均增大。

3）文中得到的液舱后板上任一点的压力峰值和比冲量的拟合计算公式可以为防护液舱结构的设计提供参考。