0 引　言

水下爆炸近场条件下高瞬态、强冲击、强辐射的流场特征给试验研究带来很大困难，测试设备不仅需要满足量程高、延迟小、频率高的基本要求，还要在强干扰、强破坏的环境下保证测试精度，并尽量保证可靠工作时间，这意味着由试验独立对水下爆炸下冲击载荷进行综合全面的研究探索比较困难。水下爆炸数值仿真研究随着计算机技术的飞跃发展不断深入，弥补了试验和理论计算的缺点，其计算快捷省力，有效避免试验高危险性、高耗资等弊端。

目前，Ansys/Ls-dyna，Dytran，Abaqus等有限元动力学数值软件被用于水下爆炸问题的求解。Geers^[1]提出一个用于自由场中水下爆炸载荷的模型，在该模型中，冲击波、气泡脉动为水下爆炸的两阶段，冲击波阶段为气泡脉动的发生创造了条件。Webster^[2]基于 Ls-dyna 中 ALE 算法仿真分析了类舰船结构在近场水下爆炸载荷作用下的冲击响应，并通过对比试验结果有效验证了仿真结果。Arami等^[3]在Dytran中使用双重渐近近似方法（DAA）研究水下爆炸作用下平板的动态响应，仿真结果表明平板由气泡脉动作用下的响应远大于冲击波作用，与试验结果相吻合。Gupta等^[4]通过 Abaqus 软件，考虑结构应变率影响，模拟在不同炸距离和材料下加筋板和不加筋板的损伤，总结出2种平板失效模型，并证实冲击波载荷峰值对结构毁伤具有影响意义。Ramajeyathilagam 等^[5]通过实验总结了4种水下爆炸的薄矩形板失效模式：I型（大变形），II*型（部分拉伸撕裂），IIa型（完全拉伸撕裂，中心点挠度增加）和III型（剪切破坏）。基于Ramajeyathilagam 等的工作，Zong 等^[6]用声结构耦合法模拟了矩形板和加强板的变形和破裂，数值结果与Ramajeyathilagam 等的实验结果接近。

本文用耦合欧拉-拉格朗日(CEL)方法在 Abaqus 求解器中对水下平板爆炸响应进行数值仿真模拟。保持爆炸距离不变，分析不同爆炸药量下平板的动态响应，为结构毁伤评估和安全预警提供依据。

1 CEL方法

水下爆炸载荷对舰船结构作用过程是一个复杂的流固耦合问题，炸药爆炸会在极短时间内产生高达GPa压力的冲击波并造成结构的巨大变形，基于网格的数值算法—拉格朗日方法和欧拉方法，都会在模拟时因网格畸变而无法较好地独立模拟水下爆炸问题。拉格朗日法在计算过程中材料粘附在网格上，材料移动则单元网格也会随之运动和变形，单元质量不变，是离散化的一种数值计算方法。欧拉法在计算过程中网格固定，材料在网格间运输，欧拉单元网格与物体不一定完全重合，每一个时间步需重新计算欧拉物质边界^[7]。为缩短计算时间需在计算前设置合适尺寸的欧拉网格。一般欧拉网格为四边形，在计算时可很好与欧拉物质边界重合，允许欧拉物质发生大变形。Abaqus中的CEL法用欧拉体积分数（Eulerian volume fraction，EVF）描述材料在欧拉域中的运动，允许在同一模型中欧拉场和拉格朗日结构之间的相互作用。该方法结合了2种网格算法的优点，在计算中拉格朗日法求解结构运动，欧拉法求解流体运动，在接触面上进行耦合计算。因此，采用CEL方法计算平板结构与水下爆炸冲击之间的相互作用。

本文CEL方法将采用欧拉算法计算在爆炸中发生大变形的水、炸药及空气域的动态响应，计算过程中，材料在固定的欧拉网格间运动，网格形状、质量大小不变；采用拉格朗日法计算平板结构的变形，计算时，网格随物质运动而变形。在接触面上，通过罚函数的施加实现欧拉体和拉格朗日体的耦合。

2 数值模型

在数值计算中，用状态方程确定CEL域中流体介质材料密度、温度、能量等和压力的关系。共涉及到水、空气、炸药爆轰产物共3种物质，状态方程的选择对模拟计算结果的准确性影响较大。

2.1 水介质状态方程

假设水的流动是黏性的、不可压缩的层流流动，采用线性Us-Up状态方程对其进行模拟，通过N-S方程精确模拟水下爆炸作用下水的运动。水的状态方程为^[8]

$ P = \frac{{{\rho _0}{\text{c}}_{\text{0}}^{\text{2}}\eta }}{{{{\left( {1 - s\eta } \right)}^2}}}\left( {1 - \frac{{\varGamma {}_0\eta }}{2}} \right) + \varGamma {}_0{\rho _0}{E_m} \text{。}$

(1)

在模拟中参数设置为水的密度 $\rho = 1\;000$ kg·m⁻³； ${c_0} = 1\;450$ m· s⁻¹；s = 0。

2.2 空气介质状态方程

空气采用理想气体状态方程：

$ P = \left( {\gamma - 1} \right)\rho {E_m} - {P_a} \text{。}$

(2)

式中： $ {P_a} $ 为外界压力，取101300 Pa； $ \gamma $ 为绝热指数，取1.4； $ {E_m} $ 为单位质量内能； $ \rho $ 为空气密度，取1.225 kg·m⁻³；气体常数R = 278；比热c_v = 716 J /(kg·K)。

2.3 炸药爆轰状态方程

炸药的状态方程应用较多的有BKW状态方程^[9]、KHT状态方程^[10]、VLW状态方程^[11]、JWL状态方程^[12]，使用Abaqus内的JWL状态方程描述TNT爆炸后的状态响应，方程如下：

$ P = A\left( {1 - \frac{{\omega \eta }}{{{R_1}}}} \right){e^{ - \frac{{{R_1}}}{\eta }}} + B\left( {1 - \frac{{\omega \eta }}{{{R_2}}}} \right){e^{ - \frac{{{R_2}}}{\eta }}} + \omega \eta {\rho _0}e \text{。}$

(3)

式中： $ \eta $ 为爆轰产物与初始炸药的密度比，即 $ \eta = \rho /{\rho _0} $ ；A，B，R_1，R_2， $ \omega $ 为JWL方程中描述炸药状态的常数；e为高含能炸药的单位质量内能。

在模拟中参数设置为：ρ=1630 kg·m⁻³，爆速6930 m· s⁻¹，A = 3.7377×10¹¹ Pa，B = 3.7414×10⁹ Pa，R₁ = 4.15，R₂ = 0.9， $ \omega $ = 0.35，e = 4.26×10⁶ J·kg⁻¹。

2.4 自由场气泡脉动验证

为进一步测试 CEL 方法处理近场水下爆炸的可靠性，研究水下爆炸对舰船及简单结构的毁伤特性，进行了自由场气泡脉动验证。

爆炸载荷来源于冲击波、气泡和空化，这3种非接触爆炸载荷都会对结构产生毁伤。从力学的角度分析，水下爆炸大致可分为4个主要过程：炸药的爆轰、冲击波的形成和传播、气泡的脉动和上浮以及冲击波在自由水面和结构的相互作用下产生的空化^[13]。图1，模拟了40g TNT爆炸形成的气泡脉动过程。

气泡在开始膨胀阶段体积变化较快，而当气泡快膨胀到最大时，气泡表面径向速度较低且持续时间较长；然后气泡开始收缩，半径迅速变小；当气泡收缩到最小时，射流将气泡撕裂。由于此时气泡内的压力高于周围水的压力，气泡再次开始膨胀。图2为射流形成阶段仿真气泡运动纵截面图，对比文献[14]中气泡脉动的试验结果，仿真结果与试验结果相吻合，验证了该方法可以很好地对水下爆炸的载荷冲击进行模拟。

3 数值结果与讨论

为模拟平板结构遭受近场水下爆炸冲击的动态响应，建立计算仿真模型如下：欧拉计算域尺寸为2 m×2 m×2 m，其中水域深度为1.5 m，水域上方空气域高度为0.5 m。将平板四周刚性固定水平沉浸在水域中，平板端面在自由液面处，欧拉的边界条件设为无反射边界。球形PEK炸药中心距离水域底部1.35 m，距离钢制平板中心0.15 m。平板采用壳体网格单元进行建模，单元尺寸为0.01 m×0.01 m。由于近场水下爆炸的特性属性，将欧拉域网格单元尺寸设置为0.025 m×0.025 m×0.025 m，以保证仿真计算的精度。

水域、空气域以及爆轰炸药采用上述状态方程及相关参数。平板材料采用Q235钢，其密度为7 800 kg·m⁻³，弹性模量取210 GPa，泊松比取0.3。钢制平板的材料本构关系由Johnson－Cook模型定义，J-C本构模型使用广泛，考虑了大变形、高应变率和高压等影响因子，具体表达形式如下^[15]:

$ {\sigma _Y} = \left( {A + B\varepsilon _p^n} \right)\left[ {1 + C\ln \left( {\frac{{\dot \varepsilon }}{{{{\dot \varepsilon }_0}}}} \right)} \right]\left( {1 - {T^{*m}}} \right) \text{。}$

(4)

式中： $ {\sigma _Y} $ 为动态屈服应力；A为静态屈服应力取2.492×10⁸ Pa；B为硬化参数取8.89×10⁸ Pa；n为硬化指数取0.746；C取0.058为应变率参数， $ \dot \varepsilon $ 为有效塑性应变率； $ {\dot \varepsilon _0} $ 为参考应变率（一般取 $ {\dot \varepsilon _0} $ =1s⁻¹）；m为温度指数取0.94。采用Johnson-Cook失效模型来模拟材料的失效，失效参数D设置为D₁=0.38；D₂=1.47；D₃=2.58；D₄=−0.0015；D₅=8.07。

为研究不同炸药量对钢制平板在同一爆距下毁伤情况的影响，仿真中固定炸药中心到平板距离为0.15 m，使用耦合欧拉－拉格朗日方法探究在近场水下爆炸下钢制平板的破裂状况。为了研究相同爆距不同爆炸药量下钢制平板的毁伤模式，仿真计算中设定炸药距离平板的距离是0.15 m，为探究水下爆炸的近场情况，利用有限元仿真方法得到不同爆炸药量下钢板破裂情况，共设置了4组工况，如表1所示。表中用冲击因子 $ 0.45\sqrt {k/d} $ 来衡量冲击波能量密度；k为药量；d为爆距。试验用PEK型炸药，是1.17倍的TNT。

图3为各工况钢制平板的失效状况。左图是试验结果^[5]，右图是仿真结果图。表2为计算与仿真结果误差对比，其误差范围均在10%以内。

可以清楚地从图3的试验结果与仿真结果看出钢制固支平板在不同炸药量下破损的情况。工况1爆炸药量为0.01 kg和工况2爆炸药量为0.02 kg时，平板仅在中心区域出现一定程度的挠度变形，并没有出现破损情况且四周也未与固支边界发生撕裂；工况3爆炸药量为0.04 kg时，试验结果中，平板结构产生一定的塑性变形，并在平板结构边缘处出现部分撕裂现象，在仿真结果中虽也出现边缘撕裂现象，但由于数值模拟没有考虑其他非实验因素，是在理想工况下进行仿真分析，故其边缘撕裂呈对称形式；工况4中平板在0.08 kg炸药量作用下已发生严重的剪切破坏，钢制平板主体从固支边界脱离，板中心出现撕裂。对比试验、与仿真结果，脱落钢板中心凸起撕裂，在水下爆炸巨大冲击作用下，钢板出现中心撕裂，并在剪切力作用下边缘撕裂整体脱落。爆炸药量的逐渐增大，钢制固支平板的最终破裂方式会发生改变，由开始的板中心发生塑性大变形到固支边界处撕裂再到板整体从固支边界脱落。

4 结　语

1）采用耦合欧拉朗格朗日方法对自由场下爆炸气泡进行了仿真，仿真结果与试验吻合良好，并对不同爆炸药量下平板的动态响应进行数值模拟，分析了不同炸药量下板的毁伤形式。

2）平板在近距水下爆炸作用下的破坏响应会随爆炸药量的增加而改变。最先发生损坏的是平板的边缘，属于易损结构。当爆炸药量较高时，平板中心会形成大破口并脱落，结构失效造成安全威胁及经济损失。在对简单舰船结构进行抗冲击性能设计时可参考本文的仿真结果。