水下爆炸冲击波作用下的空化效应数值仿真研究

宗 智，陈 岗，叶 帆，李海涛，赵延杰

（1大连理工大学 工业装备结构分析国家重点实验室 船舶与海洋工程学院，辽宁 大连 116024；2中国船舶与设计研究院，上海 200011）

1 引 言

水下爆炸主要包含两个阶段：冲击波阶段和气泡脉动阶段，但实际上除了上述两个主要爆炸载荷外，还会有空化效应引起的二次加载对结构的影响。冲击波在接触到自由液面及结构湿表面时会发生反射，使附近水域压力骤降，当水域压力下降至水的空化极限压力时，就会产生空化效应。空化效应具体可以分为区域空化和局部空化。区域空化主要是由于自由表面反射引起的反射稀疏波与入射波叠加后形成的负压区引起的，而局部空化是由于结构物的反射波与入射波叠加后形成的负压区引起的。

针对水下爆炸引起的空化现象，Aron[1]研究了区域空化的形成，理论推导了区域空化的上下边界计算公式。Makinen[2]采用平面波近似法，对比分析了4种不同的空化模型，并对简单例子分别利用近似方法和经典理论求解，比较了计算结果。Rajendran[3]分析了金属平板在遭受非接触爆炸时受到空化效应和气泡脉动的再加载现象，将试验得到的总体损伤与理论计算得到的初次冲击波损伤进行了对比分析。李海涛[4-5]利用平面波理论对局部空化和区域空化的形成特性进行了研究，并通过具体实验对局部空化理论进行了验证；理论分析了结构目标尺度的变化对空化区域形成的影响，两者符合较好；通过试验和数值仿真方法研究了气泡脉动引起的锥形空化的形成特点，初步分析了锥形空化的形成原因。许斐[6]运用了ABAQUS中内嵌的声—固耦合法对三块不同尺寸的气背圆板进行了水下冲击波载荷作用下的数值模拟。考虑了空背板材料的应变率，模型中网格的密度，水域中的空化和炸药参数对数值结果的影响。叶帆[7]模拟了水下爆炸冲击波产生的空化效应对船舶等水面结构的响应造成的影响，特别是片空化溃灭时会对结构造成较大的二次加载。

随着计算机技术的发展，现在许多著名的商业有限元计算软件都包含有类似的计及空化效应的声学有限元方法。例如NASTRAN-CFA-DAA代码[8]就采用了类似的气穴声学有限元方法。而ABAQUS软件也融入了一种类似的计及空化效应的声学有限元方法[9]，其中流体域是采用曲面波近似边界进行截断，以流体的压力作为因变量来判断是否启动空化效应的计算。

本文运用大型商业有限元程序ABAQUS，数值模拟了水下爆炸产生的空化区域的膨胀扩大、收缩溃灭，得到了水中测点的压力曲线，与实验进行了对比。然后通过空背方板的数值模拟，捕捉到了局部空化对结构的二次加载现象，并与实验值进行了对比。

2 ABAQUS水下爆炸空化效应模拟

Abaqus/Explicit计算水下爆炸时可以同时考虑空化效应的影响。Abaqus中，水域由声学单元来模拟，水节点的绝对压力为静水压力与计算所得的动压之和，如果绝对压力下降到了空化极限压力，则开始计算空化效应，该部分水介质假定开始自由膨胀而动压没有相应的下降。一旦发生了自由膨胀，水域中的压力分布将会重新计算。考虑空化效应后，水介质的压力可由下式来表述：

其中：pc为空化极限压力，p0为静水压力，pv为水节点计算所得动压，Kf为水的体积模量，x为水节点的空间坐标，uf为水节点的位移向量，θi表示i个相互独立的场变量，例如温度、空气湿度、水的盐度等等影响水介质密度和体积拖曳系数的因素。

要计及空化效应的影响，必须运用总波公式来计算水下爆炸，并且需添加两个关键字。首先需要对水域设置静水压力初始条件，此静水压力只为计算空化效应提供判断依据，而不会影响水域中动压的计算。选择所有水域节点创建一个nodeset，使用关键字*INITIAL CONDITIONS,TYPE=ACOUSTIC STATIC PRESSURE，对此nodeset施加静水压力。格式如下：

其中：Setname是水域节点nodeset的名称，P1是第一个参考点的静水压力，X1、Y1、Z1分别是第一个参考点的x、y、z坐标，P2是第二个参考点的静水压力，X2、Y2、Z2分别是第二个参考点的x、y、z坐标，根据两个参考点的静水压力，由程序自动线性计算水域节点nodeset中的所有节点的静水压力。两个参考点一般设在沿着重力方向的一条直线上。其次，空化极限压力由关键字*ACOUSTIC MEDIUM,CAVITATION LIMIT来设置，格式如下：

其中：PC—指定的空化极限压力（Pa）。

3 区域空化效应

当水下爆炸产生的冲击波到达水面之后会反射形成反射波，反射波与入射波叠加之后会使水域的绝对压力骤降，当压力下降到水的空化压力时就会出现空化现象。

图1 水下爆炸冲击波载荷几何示意图Fig.1 Underwater shock attack geometry

图2 冲击波压力的截断Fig.2 Shock wave pressure profile with cut-off

图3 区域空化示意图Fig.3 Bulk cavitation zone

为了更好地反应冲击波的反射现象，如图1所示，实爆点在自由表面的另一侧形成一个虚爆点，我们认为反射波压力由虚爆点产生，由于入射波比反射波先前到达测点，所以入射波会由于压力衰减而小于随后到达的反射波压力，引起水域绝对压力的急剧下降，表现出截断现象，一旦绝对压力到达空化极限压力，空化就会形成，见图2。图3是区域空化示意图，区域空化是以爆点为中心对称的，因此图中只显示了一半。其存在上下边界，此边界不是某时刻空化区域的边界，而表示的是整个空化过程中会发生空化的水的区域。

如图1所示，爆点与测点之间的距离为

虚爆点与测点之间的距离为

根据爆炸冲击理论中经典的经验公式，冲击波压力可以表示成

其中：t0为冲击波到达测点处的初始时刻，t≥t0；Pmax为测点处的峰值压力，单位MPa；W为炸药当量，单位kg；R为爆距，单位m；K1、K2、A、B为冲击波常数，与所使用的炸药相关；θ为冲击波衰减常数。

考虑静水压力，包括大气压和水压，记为Ph，Ph=Patm+ρgy。则任意时刻任意测点处的总压力为

可以根据Arons[1]的方法来确定区域空化的上下边界。对于上边界，取t=R2/c，令（8）式等于零；对于下边界，令测点正压力（包括入射波压力和静水压力）的衰减速度等于反射波压力的衰减速度。经过推导，得到两个方程：

以上两个方程可以使用Matlab编译程序编码，可以计算任意当量的各种类型炸药在任意爆炸深度下形成的区域空化的边界。

Marcus[10]在他的报告中做了关于区域空化的实验，如图4所示，此实验采用彭托利特炸药，药量30.8 kg，埋深21.3 m，压力传感器放置在水面下3.05 m，与爆点的水平距离51.8 m的地方，实验证明水下爆炸产生的区域空化在闭合溃灭时会产生较大的冲击压力，并提供了相关实验数据。

图4 实验示意图Fig.4 Sketch map of the test

图5 区域空化的范围Fig.5 Bulk cavitation zone

首先应有上述提到的编译代码可以计算得到此实验条件下区域空化的大致范围，如图5所示。可以看到空化区域半径约290 m，深度约7.5 m，这个数值对数值模拟时的建模具有一定的参考价值。由于空化区域是旋转体形的，关于爆点中心对称的，在Abaqus中我们里采用二维模型来模拟，选取有限元模型水域宽度300 m，深10 m，几何模型及边界条件设置如图6所示。模型A，最小单元尺寸0.1 m，共300 000个单元，模型B，最小单元尺寸0.05 m，共1 200 000个单元，为了得到空化区域的半径值，取计算时间0.2 s。由于计算时间较短，气泡脉动载荷可以忽略。本文取水的空化极限为零，即认为当水压为负压时就开始计算空化。两个模型计算得到的空化区域如图7和图8所示，深色区域是绝对压力为零的区域，即是空化区域，可以看到空化区域的扩大和缩小，直至溃灭。

从模型A和模型B计算得到的空化区域随时间变化的图示可以看出两者的趋势基本一致，但由于两者的网格密度不同，进而影响到水域压力分布，所以计算得到的空化区域稍微有点不同，同时可以看出模型A的空化区域在154 ms全部溃灭，此时溃灭处的半径大概249 m，与理论计算值相比，误差为14.1%，这是由于网格较大，水域中压力衰减较快，而模型B中的空化区域在181 ms全部溃灭。此时溃灭处的半径大概281 m，与理论计算值相比，误差为3.1%，这是由于图网格更密，模型B计算得到的空化区域显示得更细腻，衰减也较慢。

图6 边界条件示意图Fig.6 Schematic of the fluid boundary conditions

图7 模型A计算得到的空化区域Fig.7 Bulk cavitation region calculated by Model A

图8 模型B计算得到的空化区域Fig.8 Bulk cavitation region calculated by Model B

图9 模型A与模型B计算得到的测点压力时程曲线与实验对比图Fig.9 Comparison of the absolute pressure

表1 冲击波压力峰值、二次冲击时刻及压力峰值统计表Tab.1 Statistics of the shock peak value,reloading peak value and moment

模型A与模型B计算得到的测点压力时程曲线与实验测得的压力时程曲线如图9所示，模拟得到的测点处压力时程曲线与实验曲线相同，都表现出了明显的二次冲击压力，二次冲击压力出现的时间与实验基本一致，但峰值有些差别。模型A计算得到的冲击波压力峰值明显小于实验值，模型B计算得到的冲击波压力峰值接近于实验值，说明网格越密，冲击波在传播过程中衰减得越少。模型A计算得到的二次压力峰值大于模型B的计算值，但都小于实验值。得到的持续时间比实验值略长。模型A和模型B的二次冲击持续时间均比冲击波持续时间长，约为后者的2倍。

模型A中测点处进入空化发生在24.5 ms，二次冲击发生在38.9 ms时，模型B中测点处进入空化发生在24.6 ms时，二次冲击发生在38.2 ms时，观察模型A和模型B在24 ms和38 ms左右测点附近空化区域变化情况，如图10-13所示。可以看到，压力测点附近空化区域24 ms开始经过测点，38 ms发生了闭合，闭合时刻与二次冲击时刻吻合良好。

图10 模型A测点附近区域进入空化示意图Fig.10 Schematic of cavitation occur near the measuring point in Model A

图11 模型B测点附近区域进入空化示意图Fig.11 Schematic of cavitation occur near the measuring point in Model B

图12 模型A测点附近空化区域闭合示意图Fig.12 Schematic of cavitation closure near the measuring point in Model A

图13 模型B测点附近空化区域闭合示意图Fig.13 Schematic of cavitation closure near the measuring point in Model B

4 局部空化效应

水下爆炸冲击波接触结构物表面后会发生反射，入射波和反射波叠加会使流固耦合面附近的水压力下降，当达到水的空化压力时，就会出现空化现象，这样的空化称为局部空化。若结构位于浅水中，则容易发生局部空化。如果结构位于足够深度的水中，由于静水压力较大，局部空化一般不会发生。

Hammond[11]进行了一系列水下爆炸实验，以研究水中空气背衬钢板的结构响应。其将方形钢板密封固定在坚固的空心钢制方箱上，方箱由绳索连接到浮体而悬浮在水中，如图14所示。Hammond让钢板中心位于水面下2.3 m处，并将炸药置于远距离爆炸，则爆炸载荷主要由冲击波及局部空化效应构成，而避免了自由表面、海底反射、区域空化效应及气泡脉动的影响。

其论文中给出了某工况下钢板中心的壁压曲线及位移曲线，此工况采用的炸药为0.164 kg彭托利特，置于水下2.3 m，正对钢板中心，爆距3.63 m，如图15所示。目标正方形钢板边长0.78 m，板厚 3 mm，钢板密度 ρ=7 850 kg/m3，杨氏模量 E=2.11×1011Pa，泊松比 μ=0.3。

采用计算区域空化上下边界程序，计算得到此工况下的区域空化范围，如图16所示。区域空化最深不超过2 m，且空化区域距方箱有一定距离，因此确实可以忽略区域空化效应的影响。

图14 方箱布置图Fig.14 Schematic of the chain arrangement used to secure the box

图15 炸药位置示意图Fig.15 Schematic of the charge position

图16 此工况下区域空化范围示意图Fig.16 Bulk cavitation zone

此问题具有对称性，因此可以建立二分之一有限元模型。目标钢板和方箱用共节点壳单元建立，单元大小为1.5 cm，共8 112个单元。根据实验情况，近似认为目标钢板四周是刚性固定的，因此在目标钢板四周施加全位移约束。采用圆柱形水域，水域半径1.56 m，钢板附近水单元大小为1.5 cm，向外逐渐变大，共706 992个单元。以目标钢板中心为坐标原点，x轴水平沿着钢板表面，y轴指向钢板的外法向，z轴竖直向上。在对称面上施加对称条件，水域其他表面施加无反射边界条件。有限元模型如图17所示。爆炸载荷加载时间0.05 s，水的空化极限设为零。

计算得到的钢板中心位移与实验测得曲线比较如图18所示，可见板受到水下爆炸作用后发生了往复振动，而未产生塑性变形，计算得到的位移变化情况与实验基本一致。实验测得的板中心最大凹陷23 mm，计算最大值为15 mm；实验测得的板中心最大外凸11 mm，计算最大值也为11 mm。

图17 有限元模型示意图Fig.17 Schematic of the finite element model

图18 计算得到的钢板中心位移与实验比较图Fig.18 Comparison of experimental and calculated results for the displacement of the plate centre

图19 实验测得的钢板中心壁压曲线Fig.19 Experimental pressure curve

图20 计算得到的钢板中心壁压曲线Fig.20 Calculated pressure curve

图19、20分别为钢板中心壁压的实验测量曲线与计算曲线，两者的变化趋势一致，冲击波到达时壁压迅速增大，然后很快降为零，说明附近水域发生了空化现象，然后又受到了二次加载，二次加载壁压峰值约为冲击波壁压峰值的30%，且衰减速度较慢。由表2可以看到，计算得到的壁压值略小于实验值，二次加载的间隔时间也略小于实验值。总体来看，数值模拟的结果与实验是比较接近的。

表2 冲击波壁压峰值、二次加载壁压峰值及间隔时间的比较Tab.2 Comparison of the shock peak pressure,reloading peak pressure and their interval

图21显示的是计算得到的目标钢板附近水域发生局部空化的情况，深色区域为空化区域。可以看到钢板迎爆面附近水域发生了空化现象，空化区域迅速扩大，然后缩小、溃灭。在6.6 ms左右，板附近空化区域几乎全部溃灭，与二次加载时刻一致。而由图18位移曲线可以看到6.6 ms左右，无论是实验曲线还是计算曲线都显示钢板中心仍往y轴反方向移动。因此可以说明二次加载是由于局部空化在流固耦合面附近闭合溃灭引起的，而不是由于钢板回弹与水碰撞引起的。

图21 钢板附近水域局部空化情况Fig.21 Local cavitation region near the plate

实验没有测量钢板中心的速度，计算得到的钢板中心速度曲线如图22所示。从图中可以看到受到冲击波作用后钢板中心速度迅速增大，而后减小、振荡，6 ms后由于受到二次加载作用而再次大幅增大，然后又减小、振荡。这说明二次加载会对结构响应产生强烈的影响，二次加载引起的钢板中心速度峰值甚至大于冲击波引起的速度峰值。

图22 计算得到的钢板中心y向速度曲线Fig.22 Calculated y direction velocity curve of the plate centre

5 结 论

（1）Abaqus可以数值模拟水下爆炸产生的空化区域的膨胀，缩小直至溃灭的整个过程，实验证明在空化区域闭合溃灭时附近水质点会对结构造成较大的二次冲击，数值模拟也验证了这一结果，并与实验数据基本一致。

（2）水下爆炸冲击波到达自由液面后会反射引起稀疏波，稀疏波与入射波叠加在水面下形成较大的负压区，出现空化区域，根据Arons[1]的方法利用matlab编程可以预报区域空化的范围，对流体有限元模型尺寸的确立有一定的指导意义。

（3）通过对局部空化模拟得知二次加载不是由于钢板回弹与水碰撞引起的，而是因为局部空化区域在流固耦合面附近闭合溃灭，从而对周围水质点及结构施加较大的二次加载。

（4）区域空化引起的二次加载相对于冲击波虽然峰值较小，但是持续时间大约是冲击波的两倍，作用效果显著。其引起的结构速度响应可以达到冲击波相当的量级。所以水下爆炸考虑空化效应是必须的。

（5）网格划分的差异对数值模拟结果有一定的影响。较密的网格更能捕捉到空化区域的边界，更精确地模拟空化效应，所以对于水下爆炸这一瞬时高度非线性问题加密水域网格是必要的。

[1]Arons A B,et al.Long range shock propagation in underwater explosion phenomena II[C].Underwater Explosions Compendium,1943,I.

[2]Makinen.Cavitation models for structures excited by a plane shock wave[J].Journal of Fluids and Structures,1998,12:85-101.

[3]Rajendran.Reloading effects on plane plates subjected to non-contact underwater explosion[J].Journal of Materials Processing Technology,2008,206:275-281.

[4]李海涛,朱 锡,黄晓明,牟金磊.水下爆炸冲击波作用下空化区域形成的特性研究[J].高压物理学报,2008,22(2):181-186.

[5]李海涛,朱 锡,牟金磊,黄晓明.水下近距爆炸作用下弹性钢板处的空化特性研究[J].海军工程大学学报,2008,20(1):21-24.

[6]许 斐,周 力,宗 智.铝合金舰艇在水下冲击波作用下动态响应的数值研究[J].舰船科学技术,2011,33(1):31-40.

[7]叶 帆,宗 智,李海涛.水下爆炸冲击波产生的区域空化效应二次加载的并行数值研究[J].船科学技术,2011,33(9):11-17.

[8]Shin Y S,Santiago L D.Surface ship shock modeling and simulation:Two-dimensional analysis[J].Shock and Vibration,1998,5:129-137.

[9]Coupled acoustic-structural medium analysis[K].ABAQUS Theory Manual.

[10]Marcus M H.The response of a cylindrical shell to bulk cavitation loading[R].NSWC TR 81-295.1983.

[11]Lloyd H,Raphael G.Structural response of submerged air-backed plates by experimental and numerical analyses[J].Shock and Vibration,2000,7:333-341.