0 引　言

潜艇遭受鱼水雷等水中兵器爆炸攻击，主要分为接触爆炸和非接触爆炸，非接触爆炸又分为近场非接触爆炸和中远场非接触爆炸^{[1 − 3]}。中远场非接触水下爆炸作用下，潜艇除受到冲击波载荷的瞬时冲击而造成艇体的剧烈振动和结构的弹塑性变形外，还会受到气泡脉动的影响，当气泡脉动的基频与潜艇的低阶固有频率接近时，则会造成艇体严重的鞭状运动，从而导致艇上各类重要设备广泛的冲击破坏及潜艇总体强度的破损^{[4 − 6]}。

目前，国内外学者针对水下爆炸作用下对舰艇毁伤特性和结构响应研究已经做了大量工作^{[7 − 11]}。张阿漫^[12]利用声固耦合方法研究了气泡脉动、波浪和静水弯矩等外载荷联合作用下的水面舰船总体强度问题，姚熊亮等^[13]利用ALE流固耦合法模拟了冲击波载荷在圆柱壳结构周围的绕射特性，给出了冲击波与圆柱壳相互作用下的压力场分布情况。刘云龙等^[14]通过自主开发的流固耦合FSLAB软件对水面舰船在近场和中远场水下爆炸作用下，冲击响应特性进行了数值模拟，结果表明在中、远场条件下水面舰船抗冲击能力的薄弱环节是船载设备和人员。张轶凡等^[15]采用任意拉格朗日－欧拉算法研究分析了近场水下爆炸下，不同装药距离和方位对圆柱壳结构的损伤特性。嵇春艳等^[16]对大型舰船设备在水下爆炸荷载作用下的冲击环境进行数值仿真，通过结果对比分析，给出大型舰船设备冲击环境的预报公式。

在中、远场水下爆炸对结构的流固耦合毁伤问题中，常用的方法包括声学有限元法^{[17 − 18]}、谱单元D法^{[19 − 20]}、双渐近法（Doubly Asymptotic Approximations，DAA法）^[21]和二阶双渐近法（Second-order Doubly Asymptotic Approximations，DAA₂法）^[22]等。Abaqus/Explicit软件采用声固耦合方法求解流体与结构之间的冲击耦合过程，主要包括计入空化效应的总波公式和线性化的散波公式2种求解方法。本文通过采用声学单元来描述流场，将流体动压力处理为声压，并将水下爆炸载荷作为入射波加载到流固耦合界面进行计算，并据此分析中远场水下爆炸作用下，双层加筋圆柱壳结构响应特性。

1 声固耦合方法理论

使用声固耦合法简化计算。其基本思想是将流场线性化处理为声学介质，将流体动压力处理为声压，则可通过简化的方法对流固耦合问题进行求解。对于小扰动可压缩绝热流体，其控制方程可写作：

$ \frac{{\partial p}}{{\partial x}} + \gamma \left( {x,{\theta _i}} \right){\dot u^f} + {\rho _f}\left( {x,{\theta _i}} \right){\ddot u^f} = 0 。$

(1)

式中：p为流体动超压；x为流体质点的空间位置矢量；$ {\dot u^f} $和$ {\ddot u^f} $分别为流体质点速度和加速度；$ {\rho _f} $为流体密度；$ \gamma $为体积阻尼系数；$ {\theta _i} $为i个独立场变量，例如温度、气体湿度、流体盐度等。式中的达朗贝尔项写作不包含对流项形式，这对于马赫数小于0.1的问题通常具有足够精度。

假设流体为无粘，线性可压缩，于是由Bernoulli方程可得：

$ p = - {K_f}\left( {x,{\theta _i}} \right)\frac{{\partial {u^f}}}{{\partial x}}。$

(2)

式中：$ {K_f} $为流体的体积模量，$ {u^f} $为流体质点位移。

对于考虑声学介质空化效应的声固耦合法采用总波公式来实现。将式（1）两侧对x求偏导，忽略$ {\gamma / {{\rho _f}}} $，并结合式（2）对t的偏导数方程，可得到流体压力变化方程：

$ \frac{1}{{{K_f}}}\ddot p + \frac{\gamma }{{{\rho _f}{K_f}}}\dot p - \frac{\partial }{{\partial x}} \cdot \left( {\frac{1}{{{\rho _f}}}\frac{{\partial p}}{{\partial x}}} \right) = 0 。$

(3)

其中，$ {\gamma / {{\rho _f}}} $为小量的假设对于大部分的连续介质是足够精确的，但当2个相邻单元的$ {\gamma /{{\rho _f}}} $不连续时，在两单元的交界面上需特殊处理。

通过引入声压变分$ \delta p $，将式（3）在流体区域内进行积分计算，再应用格林公式，经运算得到：

$\begin{split} &\int\limits_{{\Omega _f}} {\left[ {\delta p\left( {\frac{1}{{{K_f}}}\ddot p + \frac{\gamma }{{{\rho _f}{K_f}}}\dot p} \right) + \frac{1}{{{\rho _f}}}\frac{{\partial \delta p}}{{\partial x}} \cdot \frac{{\partial p}}{{\partial x}}} \right]} {\rm{d}}V + \\ & \int_S {\delta p\left( {\frac{1}{{{\rho _f}}}n - \frac{{\partial p}}{{\partial x}}} \right)} {\rm{d}}S = 0。\end{split}$

(4)

将在获得水下爆炸载荷后，将水下爆炸载荷作为入射波加载到流固耦合界面进行计算。将计算域离散为有限个单元，单元内任一节点的声压和质点的位移、速度和加速度通过节点值和形函数插值近似得到。变分约去声压，得到结构流体运动方程：

$\begin{split} & \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{M}}_S}}&0 \\ {{\rho _f}R}&{{{\boldsymbol{M}}_f}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {\ddot U} \\ {\ddot P} \end{array}} \right\} + \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{C}}_S}}&0 \\ 0&{{C_f}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {\dot U} \\ {\dot P} \end{array}} \right\} + \\ & \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{K}}_S}}&{ - {R^{\rm{T}}}} \\ 0&{{{\boldsymbol{K}}_f}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} U \\ P \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{F_S}} \\ 0 \end{array}} \right\}。\end{split}$

(5)

式中：M_S、C_S、K_S分别为结构质量矩阵、结构阻尼矩阵和结构刚度矩阵；M_f、C_f、K_f分别为流体质量矩阵、声阻尼矩阵和流体刚度矩阵；R为流体结构耦合矩阵；U为节点位移向量；P为声压向量；F_S为结构载荷向量。

2 有限元模型 2.1 计算模型

国内外潜艇主要分为单壳和双壳，其中双壳潜艇相较于单壳潜艇有着诸多优点，如有充足的舷间空间，易于总布置，且非耐压壳体可保护耐压壳等。本研究针对典型双壳潜艇结构为目标，将其简化为双层加筋圆柱壳结构，对其在水下爆炸载荷下的鞭状运动进行研究。计算模型长为74.4 m，非耐压壳厚为10 mm，非耐压壳直径为8.6 m，耐压壳体厚度为28 mm，耐压壳直径为为7 m，肋距为0.6 m，其中每4档肋距为一个实肋板，各舱段长度为9.6 m，建立有限元模型时主要用到的是壳单元和梁单元，壳单元主要用于模拟耐压壳、非耐压壳、舱壁板等，梁单元主要用于模拟舱壁板、圆柱壳结构中加强材、实肋板等构件。壳单元类型为S4R、S3，梁单元类型为B31，杆单元类型为T3D2。

把流场分为内外2个部分构建。外流场的外部边界设定为零阻抗模拟无限流场域。外部流场、内外壳之间的水均以声学单元模拟，声学单元类型为AC3D4。外部水域由居中的圆柱和两端的半球体组成，模型位于外部水域正中间。结构模型结构形式如图1所示。外流部场及计算域模型如图2所示。

2.2 材料模型及参数

结构材料属性：弹性模量E=210 GPa；泊松比λ=0.3；密度ρ=7800 kg/m³；屈服强度σ_p=355 MPa；极限强度σ_s=560 MPa；极限延伸率取0.18。流体材料属性：声速c=1460 m/s；密度ρ=1000 kg/m³。采用Johnson-Cook本构模型定义材料的应变率效应中方程形式如下：

$ \sigma = \left( {{\sigma _0} + B{\varepsilon ^n}} \right)\left( {1 + C\ln \frac{{\dot \varepsilon }}{{{{\dot \varepsilon }_0}}}} \right)\left[ {1 - {{\left( {\frac{{T - {T_\gamma }}}{{{T_m} - {T_\gamma }}}} \right)}^m}} \right]。$

(6)

式中：$ {\sigma _e} $为Von-Mises等效塑性应力；$ {\varepsilon _e} $为等效塑性应变；$ {\sigma _0} $为参考应变率下的屈服应力；$ \dot \varepsilon $为等效应变率；$ {\dot \varepsilon _0} $为参考应变率；$ T $为温度；$ {T_m} $为金属熔点；$ {T_\gamma } $为参考温度（一般取室温）；$ B $为应变硬化系数；$ C $为应变率硬化系数；$ m $为温度软化指数；$ n $为应变硬化指数。

3 计算结果 3.1 工况设计

采用药量、水深、爆距、爆炸方位作为自变量，研究不同药量（450 kg、225 kg）、不同爆距（100 m、50 m）、不同水深（100 m、200 m）、不同爆炸方位（0、0.3、0.5）等参数下，双层加筋圆柱壳各舱的冲击环境差异，共计24个工况，具体工况设计如表1所示。

3.2 水下爆炸载荷计算

水下爆炸载荷采用Geers-Hunter模型^{[23 − 24]}计算得到。其中，水下爆炸的冲击波压力P(t)计算式为：

$ P(t) = {P_m}{e^{ - \frac{t}{\theta }}}，$

(7)

$ {P_m} = 5.24 \times {10^7}{\left( {\frac{{{W^{^{0.3333}}}}}{R}} \right)^{1.13}}，$

(8)

$ \theta = {0.45{R_0} \cdot {{\left( {\frac{R}{{{R_0}}}} \right)}^{0.45}} \cdot {{10}^{ - 3}}}，$

(9)

$ {R_0} = 0.053{W^{0.3333}}。$

(10)

式中：t为时间，s；P_m为冲击波峰压，MPa；θ为冲击波衰减时间常数；W为等效TNT当量，kg；R为爆源与模型壳板的最小距离，m；R₀为装药的等效半径，m。

气泡脉动载荷计算：

$ \left\{ \begin{gathered} {P_b} = {P_{mb}} \cdot {e^{\frac{{t - {T_m}}}{{{\theta _b}}}}}，\begin{array}{*{20}{c}} {{T_m} - {t_2} < t \leqslant {T_m}} ，\end{array} \\ {P_b} = {P_{mb}} \cdot {e^{ - \frac{{t - {T_m}}}{{{\theta _b}}}}}，\begin{array}{*{20}{c}} {{T_m} < t < {T_m}} \end{array} + {t_2} 。\\ \end{gathered} \right. $

(11)

$ {P_{mb}} = \frac{{39 \times {{10}^6} + 24{P_0}}}{{{{\bar R}_{bc}}}}，$

(12)

$ {T_m} = 2.01\frac{{{W^{0.3333}}}}{{Z_0^{0.8333}}}，$

(13)

$ {Z_0} = H + 10.3 ，$

(14)

$ {t_2} = 1645\frac{{{R_0}}}{{P_0^{0.71}}}。$

(15)

式中：P_b为气泡脉动压力，MPa；θ_b为气泡衰减时间常数，s；t₂为气泡脉动压力正压作用时间的一半，s；t为气泡脉动压力作用时间，s；T_m为气泡脉动周期，s；Z₀为气泡等效深度，m；R_bc为气泡脉动中心距离到模型壳板的最小距离，m；H为潜深，m。

3.3 模态分析

对于冲击问题而言，影响其结果的一方面是输入载荷的特性，另一方面是结构本身的固有属性。其中，影响圆柱壳鞭状运动的固有属性则主要是结构低阶的模态频率。为此，计算得到的双层圆柱壳结构一二低阶模态频率如图3所示。

对应的各阶模态固有频率如表2所示，可看出双层圆柱壳模型的干模态和湿模态存在较大差距，其湿模态与常规模型结构一阶模态频率较接近，当气泡脉动周期同结构固有周期处于同一量级，会诱发鞭状运动，加剧结构破坏程度。

3.4 水下爆炸载荷诱导鞭状运动研究 3.4.1 不同药量影响

以100 m水深、50 m爆距的2组工况为例分析不同药量对双层加筋圆柱壳结构的影响。其中，225 kg当量TNT工况对应工况3，450 kg当量TNT工况对应工况9。从图4可看出，工况9产生的冲击波峰值载荷相较于工况3提升约30%，第一次的气泡脉动载荷增大约18%，气泡脉动周期增大约25%。

由图4可知，模型中部位移响应周期与气泡载荷曲线相对应，在冲击波载荷作用后，模型开始产生形变，并在0.1 s左右达到变形最大值。且大当量工况下，模型中部的最大位移远大于小当量情况。在模型到达最大变形后，变形逐渐恢复，由于模型自身的阻尼振荡，模型在相反方向的变形很小，只略大于平衡时的位置。在气泡达到第一个脉动周期时，模型完成一个完整鞭状运动周期，在第一个气泡脉动载荷和模型的回落运动叠加作用下，模型发生第二次鞭状运动。2种工况的圆柱壳到达各波峰波谷时的变形和应力云图，如图5和图6所示。为便于观察现象，将图中的结构变形放大了50倍显示。由图5可知，工况3的模型在水下爆炸前期的冲击波作用下就已被激起了一阶鞭状运动，模型的应变主要集中在圆柱的上下表面。由图6可知，大当量工况下圆柱总变形及局部应变远大于小当量工况，应变峰值增大30%左右。

3.4.2 不同爆距影响

保持药量、爆炸水深、爆炸方位参数不变，研究不同爆距对模型爆炸冲击响应的影响，图7和图8分别为225 kg TNT当量和450 kg TNT当量下，不同爆距的计算结果。从不同爆距下计算得到的载荷输入可知，随着爆距增加1倍，冲击波载荷、气泡脉动载荷下降了约1倍，但气泡脉动周期主要受爆源自身特性影响，因此气泡脉动周期几乎不变。

由图可知，在远场爆炸载荷减小的情况下，模型变形变小，同时发现模型中部的位移响应幅值均有明显缩减，但位移响应曲线的周期变化不大。

3.4.3 不同水深影响

在保证其余变量不变的情况下，调整爆炸水深，研究不同水深对模型爆炸冲击响应的影响。图9和图10分别为225kg TNT当量和450kg TNT当量下不同水深的计算结果。225kg TNT当量对应工况3和工况15，450kg TNT当量下对应的工况9和工况21。

随着模型爆炸时所处水深不同，模型所受到的预应力也不相同，此外，爆源所处深度不同时，气泡脉动周期也不相同，在二者的叠加下，会影响模型的爆炸冲击响应。由图9可知，随着水深增加，对水下爆炸载荷的幅值影响不大，但第一次气泡脉动周期从0.245 s缩减到0.15 s，此时爆炸载荷的频率从4.1 Hz增加到6.7 Hz，而双层加筋圆柱壳模型的一阶湿模态频率在2.65 Hz。工况3对应的载荷频率与模型一阶湿模态频率相近，因此在载荷幅值相近的情况下，工况3的模型中部变形位移更大，与气泡载荷周期相对应，呈明显的一阶振型，而工况15的载荷频率更高，模型模型的位移响应幅值小于工况3的位移响应幅值。由图10可知，气泡的第一次脉动周期从0.31 s缩减到0.18 s，此时爆炸载荷的频率从3.2 Hz增加到5.6 Hz，同样工况9更接近双层加筋圆柱壳模型的一阶湿模态频率，因此在共振的作用下，模型位移响应更大。

3.4.4 不同爆炸方位影响

为研究不同爆炸方位对模型爆炸冲击响应影响，保持爆炸药量不变为225 kg，爆炸水深不变为100 m，爆距不变为50 m，分别对应工况1、工况2和工况3。图11和图12分别为提取不同舱室的位移曲线。

由图11可知，在当前中远场爆炸工况下，各舱位移曲线变化周期相近。其中，爆炸方位为0时，爆点位于模型头部；爆炸方位为0.5时，爆点位于模型中部。这2种工况在1舱和4舱的位移幅值差距最大，峰值差距达30%。当爆炸方位为0.3时，爆点位于3舱之间，位移响应位于2种工况之间，离爆点较近的2舱位移响应略大于另外2个工况。图12给出了450 kg TNT当量下不同爆炸方位的计算结果，在大当量情况下，基本规律与225 kg TNT当量下工况相似。在爆点靠近模型中部时，模型的鞭状运动最为显著，从4舱和3舱的位移曲线可看出，爆点居中的工况下，模型前2个周期的位移幅值相近，衰减较小，具有明显地周期特征。而当变形载荷传递到模型两侧时，从1舱响应曲线可看出，模型位移响应曲线周期性减弱。当爆点在模型一端工况下，模型本身自身鞭状运动特点减弱，传递到模型中部的载荷也相对较小，位移响应曲线衰减很快，从爆点距模型最近的1舱的响应曲线可看出，此时模型位移响应在爆炸载荷直接作用下，具有相对明显地周期性特征，但响应幅值衰减明显。

3.5 模型冲击环境预报

在中远场爆炸下，模型的结构变形虽不大，但模型上装载的设备系统可能承受较大冲击加速度，而简单的位移变化并不能完全反映真实的冲击过程，因此需对不同工况下的模型冲击环境进行预报。

描述冲击运动一般有2种不同的方法：一种是冲击响应的时域描述；另一种是冲击响应的频域描述，即Fourier变换或冲击谱分析。冲击谱的优点在于可将冲击的剧烈程度进行量化，参数化描述冲击环境的恶劣程度。从德国BV标准到国军标舰载设备的考核，设备抗冲击的标准是以冲击谱参数的形式进行规定，可见冲击谱在冲击环境分析及设备考核中都有至关重要的作用。根据模型重要部位的加速度时历曲线，加速度测点如图13所示，计算得到典型冲击谱如图14所示。

由于计算工况比较多，典型分析了各工况的谱加速度幅值，如表2所示。从各工况的谱加速度幅值对比可以发现，当药量从225 kg增加到450 kg时，谱加速度幅值增加约50%；当爆距从50 m增大到100 m时，水下爆炸载荷幅值减小，谱加速度幅值减小100%；当爆炸水深从100 m增大到200 m时，水下爆炸载荷幅值变化不大，因此谱加速度幅值变化不大；当水下爆炸爆炸方位从模型一侧移向模型中部时，谱加速度变化规律较复杂，当爆点位于模型一侧和模型中部时，相对应舱室处的谱速度幅值明显大于其余舱室，当爆点位于模型一侧1/3处时，谱速度幅值变化较小。

4 结　语

1）当水下爆炸药量从225 kg增大到450 kg时，水下爆炸载荷方面：水下爆炸冲击载荷和气泡脉动载荷约增大30%和18%，气泡脉动周期增大25%。模型冲击响应方面：在中远场小变形形况下，模型变形仍属弹性应变范围内，最大应变增大30%，模型鞭状运动周期与气泡脉动周期基本保持一致。冲击环境方面：谱加速度幅值增大约50%。

2）当水下爆炸爆距从50 m增大到100 m时，水下爆炸载荷方面：随着爆距增加1倍，冲击波载荷、气泡脉动载荷下降了约1倍，但气泡脉动周期主要受爆源自身特性影响，因此气泡脉动周期几乎不变。模型冲击响应方面：模型中部的位移响应幅值均有明显缩减，但位移响应曲线的周期变化不大。冲击环境方面：谱加速度幅值减小100%。

3）当水下爆炸水深从100 m增大到200 m时，水下爆炸载荷方面：随着爆炸水深增大，水下爆炸载荷峰值变化不大，但气泡脉动周期减小约40%。模型冲击响应方面：尽管模型在静水压作用下，所受预应力增大1倍，但模型鞭状运动受到水下载荷周期性影响显著，由于100 m水深下，水下爆炸气泡脉动周期与模型一阶湿模态相近，共振的作用导致模型位移响应比200 m水深更大，且位移响应曲线周期性更显著。冲击环境方面：谱加速度幅值受冲击载荷幅值影响，变化不大。

4）当水下爆炸方位改变时，水下爆炸载荷保持不变。在爆点位于模型中部附近时，模型变形与模型湿模态一阶振型相同，此时模型鞭状运动最显著，模型中部的位移响应大于模型两侧位移响应，且曲线周期性显著，位移响应曲线幅值衰减较小。当爆点位于模型首部时，模型鞭状运动减小，在爆点较近的舱室，位移响应幅值大于其他工况，但位移响应曲线周期性减弱，位移响应曲线幅值衰减较大。冲击环境方面：当水下爆炸爆炸方位从模型一侧移向模型中部时，谱加速度变化规律较复杂，当爆点位于模型一侧和模型中部时，相对应舱室处的谱速度幅值明显大于其余舱室，当爆点位于模型一侧1/3处时，不同舱室的谱速度幅值变化较小。