港池环境近水面水下爆炸特性及其毁伤效应

董琪; 韦灼彬; 唐廷; 李凌锋; 刘靖晗

doi:10.11858/gywlxb.20180638

摘要: 为研究港池环境近水面水下爆炸载荷及其对码头结构的损伤特性，设计了一种典型码头结构，并构建港池环境，运用LS-DYNA程序开展水下爆炸数值模拟研究，对爆炸现象、荷载特性、结构动态响应和能量吸收特性4个方面进行了详细研究，分析了边界、比例爆距等参数的影响规律。结果表明：爆炸气泡脉动主要受到码头结构边界和水面的影响，水底和有限港池内的流体运动对其亦有一定影响；冲击波荷载以比例爆深为中心呈垂向对称分布，气泡脉动荷载主要分布于比例爆深以下位置；结构变形和毁伤主要在冲击波传播阶段形成，气泡脉动和射流的二次毁伤效果较弱；混凝土和沉箱内填土是主要能量吸收部分。

关键词:

Abstract: In order to study the load of underwater explosion near the surface of the harbor basin and the damage effect on the wharf, we designed a typical wharf structure and built a harbor basin environment. Then through a series of numerical analysis which was accomplished based on finite element program LS-DYNA, the explosion phenomenon, loading characteristics, structural dynamic response and energy absorption characteristics were studied in details, the influence rules and action mechanism of the boundary, scaled collapse distance and other parameters were analyzed. The results show that: the explosion bubble pulsation is mainly affected by wharf structure boundary and free surface, the bottom and the movement of fluid in the limited space also have some impacts. The shock wave load is symmetrically distributed in the vertical direction with the scaled explosion depth as the center, and the bubble pulsation load is mainly distributed below the scaled explosion depth. The structural deformation and damage are mainly formed on the propagation of shock wave, and the secondary damage effect of bubble pulsation and jet is weak. Concrete and caisson fill absorb most of the explosive energy.

Key words:

冲击载荷下延性金属的层裂破坏是一个复杂的损伤演化过程，目前被广泛认同的延性金属损伤发展模式是孔洞首先在材料内部成核，随后孔洞长大，当孔洞长大到一定程度时，孔洞之间开始聚集，最终在材料内部形成宏观断裂面^[1-2]，整个过程和物理机制十分复杂。孔洞聚集是层裂过程中的关键环节，因此，研究孔洞聚集行为对揭示材料层裂的整个物理过程和本质具有重要意义。

近年来，研究人员针对影响层裂行为的多种因素进行了研究，如外载荷条件^[3-4]、材料微观结构^[5]等。Curran等^[6]综合各种实验结果，总结出两种微孔洞的聚集方式：一种是孔洞之间的直接聚集，由塑性应变场相互作用造成，一般发生在延性金属中；另一种是孔洞之间材料的局部塑性失稳，孔洞之间出现一个很窄的局部化区域，孔洞通过该区域连接，从而发生聚集。同时，也有学者为了描述孔洞成核和长大过程，建立了一些微观层裂模型，如NAG模型^[1]、VG模型^[7]等。目前，孔洞聚集的物理机制研究相对较少，部分学者^[8-9]通过逾渗理论描述微孔团簇之间的相互作用，但这也是基于平均场理论建模的。

随着计算能力的不断提高，多尺度数值模拟方法越来越多地运用于延性金属中的孔洞聚集研究中。Tu等^[10]通过分子动力学方法模拟了单晶铝中两个孔洞的聚集过程，结果表明，随着压缩应变的增大，屈服强度呈现先降低后升高的趋势，压缩阶段孔洞表面产生的位错阻碍了其在拉伸阶段的聚集。Yang等^[11]对铝中的孔洞聚集进行了分子动力学模拟，结果表明，孔洞聚集的持续时间随温度的升高而延长。Zhao等^[12]对高应变率拉伸载荷下单晶铜中孔洞长大和聚集进行了分子动力学模拟，结果表明，应力三轴度的峰值先随初始韧带距离（initial ligament distance，ILD）的增大而增大，直到达到孔洞聚集的临界点，随后减小，应力三轴度可能是表征孔洞是否发生聚集的一个重要指标。王永刚等^[13]采用有限元分析方法，研究了延性金属动态拉伸断裂过程中孔洞的贯通行为，发现孔洞之间开始贯通的临界韧带距离约为0.5倍孔洞的实时直径，并指出临界孔洞间韧带距离判据和临界孔洞体积分数判据在物理本质上是一致的。Hure^[14]采用有限元方法，利用有效各向同性屈服应力，提出了一个多孔单晶材料的聚集准则。Holte等^[15]通过对应变梯度增强材料基体中嵌入离散空隙的单胞模型进行有限元分析，结果表明，在没有应变梯度的影响下，减小孔洞间韧带距离会增强孔洞的聚集倾向。

由于实验过程中难以实时观测材料的内部信息，人们对孔洞聚集行为的认识仍不够深刻，数值模拟方法便成为研究层裂过程中孔洞聚集行为的有效途径之一。考虑到分子动力学模拟技术所涉及的尺度相对较小，无法真实还原层裂实验中孔洞聚集至靶板断裂的整个损伤演化过程，为此，我们早期通过有限元数值模拟，以厚度比为1∶2建立了一维应变下的平板撞击模型^[16]，进行了层裂过程中的损伤演化分析。但是，对于随机孔洞的损伤演化而言，其聚集行为仍很难用定量的方式剖析孔洞之间的相互作用。为定量分析孔洞在层裂过程中的聚集行为，本研究拟建立平板撞击条件下高纯铜层裂损伤演化的有限元模型，重点探讨孔洞间的初始韧带距离、孔洞直径以及孔洞分布位置对孔洞聚集行为的影响规律，为建立合理的孔洞聚集模型奠定基础，进一步揭示延性金属层裂的物理机制。

1. 实验现象

在延性金属的初始层裂实验中，通过软回收靶板，观察到孔洞之间的聚集现象。图1(a)、图1(b)、图1(c)是早期实验照片^[17]，可以观察到相同孔洞（孔径相差不大）、不同孔洞以及三孔洞聚集后的现象；从图1 (d)^[18]中可以观察到不同位置的孔洞聚集后的现象（图中均用红圈标出）。基于这些实验现象，本研究将通过有限元方法进行对比分析。

图 1 软回收靶板中孔洞的光学照片：(a) 相同孔洞，(b) 不同孔洞，(c) 3个孔洞，(d) 不同位置孔洞

Figure 1. Optical photo of voids in soft recycling target: (a) same voids, (b) different voids, (c) three voids, (d) voids in different position

下载: 全尺寸图片幻灯片

2. 有限元模型

在ABAQUS有限元软件中建立高纯铜飞片和靶板模型，为减少计算时间、节约计算成本，这里建立的飞片和靶板模型均为二维模型，所建模型在厚度方向上与实验保持一致，如图2所示，其中飞片的宽度为2.0 mm，厚度为0.4 mm，靶板的尺寸为4.0 mm×0.4 mm。为了模拟一维应变条件，对飞片和靶板施加轴对称周期边界条件。在对称碰撞条件下，由波传播理论分析可知，当飞片以一定的速度撞击靶板时，接触面会向飞片和靶板分别传入冲击压缩波，冲击压缩波分别在飞片和靶板的自由面上反射形成卸载波，相向的卸载波在靶板中相遇，产生拉伸应力区。当飞片和靶板的厚度比为1∶2时，靶板中心及其附近区域拉伸应力的持续时间最长，损伤主要集中在该区域内，因此在靶板中心位置设置边长为0.4 mm的正方形作为损伤区域。在损伤区域的中心布置微米级的填充孔洞作为初始成核源。

图 2 高纯铜层裂的几何模型

Figure 2. Geometric model of high-purity copper spallation

下载: 全尺寸图片幻灯片

在压缩过程中，孔洞具有与周围靶板相同的材料性能，孔洞的成核应力 ${p_{\rm c}}$ ^[19]可以估算如下

${p\rm_c} = \frac{2}{3}\left( {{\sigma _0}+\frac{{{k\rm_y}}}{{\sqrt {{d\rm_G}} }}} \right)\left[ {1 - \ln \frac{3}{2}\left( {\frac{{{\sigma _0}}}{E}+\frac{{{k\rm_y}}}{{E\sqrt {{d\rm_G}} }}} \right)} \right]$

(1)

式中： $E$ 为杨氏模量； ${d\rm_G}$ 为晶粒度，这里设为78 μm^[16]； ${\sigma _0}$ 和 ${k_{\rm{y}}}$ 为模型参数，对于高纯铜， ${\sigma _0}$ =200 MPa， ${k\rm_y}$ =0.14 MPa·m^1/2。

计算中，高应变率加载下高纯铜材料的本构关系采用理想流体弹塑性本构关系，即

${\sigma _{ij}} = - p+{S_{ij}}$

(2)

式中： ${\sigma _{ij}}$ 为应力张量分量， $p$ 为静水压力， ${S_{ij}}$ 为t时刻下的偏应力张量分量。若von Mises屈服函数满足

$\phi = \frac{1}{2}{S_{ij}}{S_{ij}} - \frac{{\sigma _{\rm{y}}^2}}{2} \leqslant 0$

(3)

则偏应力张量为

$S_{ij}^{t+1} = 2G\left( {{\varepsilon _{ij}} - {\varepsilon _v}} \right)$

(4)

反之，若von Mises屈服函数满足

$\phi = \frac{1}{2}{S_{ij}}{S_{ij}} - \frac{{\sigma _{\rm{y}}^2}}{2} > 0$

(5)

则偏应力张量为

$S_{ij}^{t+1} = \frac{{{\sigma _{\rm{y}}}}}{{{S^*}}}{S_{ij}},\;\;\;\;\;\;{S^*} = {\left( {\frac{2}{3}{S_{ij}}{S_{ij}}} \right)^{1/2}}$

(6)

式中： $S_{ij}^{t+1}$ 为t+1时刻下的偏应力张量分量， $G$ 为剪切模量， ${\varepsilon _v}$ 为体积应变， ${\sigma \rm_y}$ 为屈服应力。 ${\sigma \rm_y}$ 满足Johnson-Cook本构，即

${\sigma }_{\text{y}}=\left(A+B{\varepsilon }_{\rm p}^{n}\right)\left(1+C\mathrm{ln}\,{\dot{\varepsilon }}_{\rm p}^{*}\right)\left(1-{T}^{*}{}^{m}\right)$

(7)

式中： ${\varepsilon \rm_p}$ 为等效塑性应变； ${\dot{\varepsilon }}_{\rm p}^{*}\text=({\dot{\varepsilon }}_{\rm p}/{\dot{\varepsilon }}_{0})$ ， ${\dot \varepsilon \rm_p}$ 为等效塑性应变率， ${\dot \varepsilon _0}$ 为参考应变率； $A$ 为参考应变率和参考温度下的初始屈服应力； $B$ 为材料应变硬化模量； $C$ 为材料应变率强化参数； $n$ 为硬化指数； $m$ 为材料热软化指数； ${T}^{*}=(T-{T}_{\text{room }})/({T}_{\text{melt }}-{T}_{\text{room }})$ ，T_room为室温，T_melt为熔点。采用Johnson-Cook损伤模型模拟材料的损伤和失效，即

${\varepsilon _{\rm f}} = \left[ {{d_1}+{d_2}\exp \left( {{d_3}{\sigma ^*}} \right)} \right]\left( {1+{d_4}\ln {{\dot \varepsilon }\rm_p}^*} \right)\left( {1+{d_5}{T^*}} \right)$

(8)

式中： ${\sigma ^*} = - p/{\sigma _{\rm eq}}$ ，其中 ${\sigma _{\rm eq}}$ 为Mises等效应力； ${d_1}$ 、 ${d_2}$ 、 ${d_3}$ 、 ${d_4}$ 、 ${d_5}$ 为材料损伤系数。

静水压力 $p$ 由Mie-Grüneisen状态方程给出

$p = \frac{{K\eta }}{{{{(1 - s\eta )}^2}}}\left( {1 - \frac{{{\varGamma _0}\eta }}{2}} \right)+{\varGamma _0}{\rho _0}{e\rm_m}\eta = \frac{\rho }{{{\rho _0}}} - 1$

(9)

式中： $K$ 为体积模量， $\rho$ 为当前的密度， ${\rho _0}$ 为材料初始密度， $s$ 为材料常数, ${\varGamma _0}$ 为材料的Grüneisen系数， ${e\rm_m}$ 为材料的比内能。高纯铜的J-C本构、损伤及状态方程参数见表1。

表 1 高纯铜的Johnson-Cook本构、损伤及状态方程参数^[20-21]

Table 1. Johnson-Cook constitutive, damage and equation of state parameters of high-purity copper^[20-21]

ρ/(kg·m^－3)	A/MPa	B/MPa	C	m	n	${\dot \varepsilon}_0$ /s^–1	T_room/K	T_melt/K
8910	90	292	0.025	1.09	0.31	1	300	1356

C₀/(m·s^－1)	s	Γ₀	G/GPa	d₁	d₂	d₃	d₄	d₅
3910	1.51	2	46.6	0.54	4.89	3.03	0.014	1.12

下载: 导出CSV

| 显示表格

3. 计算结果与分析

3.1 双孔的聚集行为

3.1.1 初始韧带距离的影响

采用上述有限元模型对两个相同孔洞的聚集行为进行数值模拟。初始成核源孔洞直径设置为5 μm，如图2(a)所示，飞片的撞击速度设置为100 m/s，对这两个孔洞在不同孔洞间初始韧带距离条件下的聚集行为进行分析，在一系列模拟中，孔洞间的初始韧带距离d_IL0分别设置为20、30、40和50 μm。王永刚等^[13]的研究表明：在孔洞独立长大阶段，孔洞为“球对称”长大，即图1(a)中孔洞上A、B两点的位移是一致的；A、B两点的位移发生明显改变，表明孔洞之间的塑性应变场交叠影响了B点的位移增长速率，判断此时孔洞开始聚集。其中，B点的位移增长速率减小直接影响孔洞直径的变化，因此通过实时跟踪孔洞长大、聚集过程中孔洞直径的变化可以判断孔洞何时发生聚集。

以孔洞间初始韧带距离d_IL0=20 μm的两孔为例，图3给出了孔洞上A点和B点的位移时程曲线以及孔洞直径增长速率时程曲线。从A、B两点的位移时程曲线可以发现：t=1.57 μs左右时，A、B点的位移确实存在一个偏差点，此时孔洞开始发生聚集。从孔洞直径增长速率曲线可以发现：孔洞直径增长速率起初缓慢增长，随后进入近似线性增长阶段，此阶段的直径增长加速度K_g约为2.995 Gm/s²；当A、B两点位移出现偏差时，对应的孔洞直径增长速率出现非线性变化，此阶段的直径增长加速度K_c约为1.717 Gm/s²。孔洞间塑性应变场的交叠会影响孔洞直径的增长，孔洞在长大和聚集过程中的直径增长加速度存在明显变化，以此可判定孔洞的聚集情况；在t=1.63 μs后，孔洞直径出现了减速增长的变化趋势，对应图4中两孔发生聚集的一侧应变局部化严重，受两孔之间的应力联结作用，最终导致两孔贯通。

图 3 点A和点B的位移时程曲线以及孔洞直径增长速率时程曲线

Figure 3. Time history curves of displacement of points A and B and void diameter growth rate

下载: 全尺寸图片幻灯片

图 4 不同时刻下d_IL0=20 μm的孔洞聚集演化云图

Figure 4. Evolution distribution of voids coalescence in d_IL0=20 μm at different time

下载: 全尺寸图片幻灯片

图4给出了不同时刻d_IL0=20 μm的孔洞聚集的等效塑性应变演化云图。可见，孔洞的长大和聚集强烈依赖塑性应变。图4中两孔洞各自形成几乎均匀的球对称塑性应变场；当t=1.57 μs时，孔洞开始发生聚集，孔洞周围的塑性应变场相互交叠，孔洞仍为规则的圆形；当t=1.63 μs时，对应图3中孔洞直径的增长速率呈非线性增长，孔洞之间的塑性应变场严重交叠，阻碍了孔洞的增长，孔洞形状也发生明显变化；当t=1.65 μs时，孔洞之间发生贯穿。孔洞周围产生的塑性应变场相互作用形成新的塑性流动路径，这一路径最终发展为宏观的断裂面。图4显示的结果与图1(a)的实验结果相似。

图5(a)给出了不同孔洞间初始韧带距离下孔洞直径增长速率时程曲线，根据上述判据判断孔洞之间开始聚集的起始时刻（图5(a)中用虚线标出）。当初始韧带距离d_IL0为20、30、40、50 μm时，孔洞长大阶段的直径增长加速度相差不大，分别为2.995、2.770、2.315、2.117 Gm/s²，发生聚集后的加速度均出现骤降，分别为1.717、1.483、0.984、0.602 Gm/s²。随着d_IL0的增大，孔洞直径增长速率峰值持续时间也越来越长，出现减速增长趋势的时间点也越来越延后，意味着孔洞贯通时间也延后。图5(b)将图5(a)中各组孔洞开始聚集的起始时间t_n与d_IL0=20 μm的孔洞开始聚集的起始时间t₁进行归一化，结果显示：随着d_IL0的增加，孔洞之间开始聚集的起始时间逐渐延后。图6给出了不同d_IL0条件下实时孔洞间的韧带距离d_IL随时间的演化曲线，根据孔洞之间开始聚集的起始时间，容易确定不同d_IL0条件下孔洞开始聚集的临界韧带距离（d_ILc），可见，d_ILc随d_IL0的增大而增大。

图 5 不同d_IL0下孔洞直径增长速率时程曲线以及归一化聚集时间曲线

Figure 5. Time history curve of void diameter growth rate with different d_IL0 and curve of normalized coalescence time

下载: 全尺寸图片幻灯片

图 6 实时韧带距离演化时程曲线

Figure 6. Time history curves of evolution of ILD

下载: 全尺寸图片幻灯片

3.1.2 孔洞位置分布的影响

在层裂实验中，成核孔洞随机分布在材料内部，本研究通过保持孔洞直径不变、改变孔洞位置，探究孔洞位置分布对孔洞聚集的影响。设置孔洞分别为水平（0°）分布、垂直（90°）分布和呈45°角分布，如图2(a)所示，孔洞间初始韧带距离d_IL0均设置为30 μm，飞片的撞击速度为100 m/s。统计3种位置分布方式的孔洞直径增长速率，结果如图7所示。从图7中可以看到：0°分布和90°分布的两个孔洞直径增长速率无明显差异，而45°分布的两个孔洞在聚集后期存在差异；在孔洞长大阶段，45°分布的孔洞直径增长加速度最大，约为3.179 Gm/s²，90°分布方式次之，约为2.900 Gm/s²， 0°分布方式最慢，约为2.770 Gm/s²；当孔洞直径增长速率出现非线性变化时，45°分布、90°分布、0°分布的孔洞直径增长加速度分别降至2.112、1.500和1.483 Gm/s²。通过上述孔洞发生聚集的判据，可以发现，45°分布、90°分布、0°分布孔洞出现聚集的时间分别为1.57、1.60和1.58 μs，即在这3种孔洞位置分布中，45°分布方式的孔洞更容易发生聚集。

图 7 不同位置分布的孔洞直径增长速率时程曲线

Figure 7. Time history curves of voids diameter growth rate with different arrangement

下载: 全尺寸图片幻灯片

图8为3种几何构型孔洞聚集过程的演化云图。0°分布的两个孔洞左右两侧的塑性应变场几乎呈对称分布，孔洞间的拉伸破坏使孔洞间出现贯通的损伤带，最终导致两孔洞贯通；90°分布的两个孔洞的塑性应变场则几乎呈上下对称，最终在竖直方向上出现孔洞贯通；45°分布的两个孔洞的塑性应变场在垂直于45°方向上几乎呈对称分布，最终两孔洞沿45°方向贯通。图8显示的结果与图1(d)中90°和45°分布孔洞的实验结果相似。

图 8 不同时刻下不同位置分布的孔洞聚集演化云图

Figure 8. Evolution distribution of voids coalescence in different arrangement at different time

下载: 全尺寸图片幻灯片

3.1.3 孔洞直径的影响

层裂实验中，成核孔洞的大小是随机的，为了探究随机大、小两个孔洞之间的聚集行为，模拟了孔径不同的两个孔洞之间的聚集行为，保持其中一个初始成核孔洞的直径d₁为5 μm不变，另一个孔洞直径d₂分别设置为1、3和7 μm，则孔径比分别为1/5、3/5、7/5，孔洞间的d_IL0设置为30 μm，飞片的撞击速度设置为100 m/s。

图9(a)给出了不同孔径比条件下各孔洞直径增长速率时程曲线。可以发现：孔洞直径增长速率均表现为先缓慢增长，随后进入近似线性增长阶段，接着出现非线性增长阶段，最后孔洞直径减速增长。当孔径比为1/5时，小孔和大孔的直径增长速率差异明显，随着孔径比的逐渐增大，这种差异越来越不明显，说明大孔抑制小孔的生长，但是这种抑制效果随孔径比的增大而变弱。根据上述孔洞发生聚集的判据发现：当孔径比为1/5时，小孔出现聚集的时间约为1.57 μs，直径增长加速度约为2.228 Gm/s²，大孔出现聚集的时间约为1.59 μs，直径增长加速度约为0.869 Gm/s²；当孔径比为3/5时，小孔出现聚集的时间约为1.58 μs，直径增长加速度约为1.166 Gm/s²，大孔出现聚集的时间约为1.59 μs，直径增长加速度约为1.141 Gm/s²；当孔径比为7/5时，大孔和小孔出现聚集的时间均约为1.58 μs，直径增长加速度约为1.441 Gm/s²。统计孔径比为1/5的孔洞在不同时刻下水平长度与竖直长度的比值（x/y，简称横纵比），如图9(b)所示，结果显示：小孔的横纵比变化明显大于大孔的横纵比变化，即小孔更早出现破坏，这一结果在图10中表现为靠近小孔一侧的韧带破坏更加严重。

图 9 不同孔洞的直径增长速率时程曲线及孔洞横纵比曲线

Figure 9. Time history curves of different voids diameter growth rate and curves of voids horizontal to vertical ratio

下载: 全尺寸图片幻灯片

图 10 不同时刻不同孔径比的孔洞聚集演化云图

Figure 10. Evolution distribution of voids coalescence in different aperture ratios at different time

下载: 全尺寸图片幻灯片

图10为不同孔径比条件下孔洞聚集的演化云图。从图10中可以看到：当孔径比为1/5时，随着时间的推移，大孔几乎呈球形变化，小孔出现了严重的破坏变形。在初期阶段，孔洞周围形成几乎呈球形的塑性应变场，但是这两个孔洞的塑性应变场不对称，与相同孔洞之间韧带的均匀破坏不同，不同孔径两孔之间的韧带出现小孔被大孔拖拽，发生非对称破坏的现象。图10显示的结果与图1(b)中大小孔洞聚集的实验结果相似。

3.2 三孔聚集机制

基于上述两孔聚集行为的模拟和分析，采用三孔几何模型探讨孔洞聚集时的竞争机制。当两孔发生聚集时，孔洞只有一侧与另一孔洞相互作用；当3个孔洞聚集时，中间孔洞的左右两侧均会受到来自其他两个孔洞的作用，导致中间孔洞与附近孔洞发生聚集时存在一定的先后顺序和竞争关系，以下将通过改变孔洞间初始韧带距离和孔洞直径来研究三孔聚集时的竞争机制。

3.2.1 初始韧带距离的影响

三孔的几何模型相比于上述两孔的几何模型做了一些调整，在损伤区域内设置3个直径相同的孔洞，直径大小为5 μm，如图2(b)所示。通过改变孔洞间初始韧带距离来研究孔洞之间聚集的竞争机制，这里选取的3组孔洞中，孔洞1与孔洞2之间的韧带距离统一设置为30 μm，孔洞2与孔洞3之间的韧带距离分别设置为30、50、70 μm，飞片的撞击速度设置为100 m/s。

图11给出了不同孔洞间初始韧带距离下各孔洞的直径增长速率曲线。从图11可以看出：当3个孔洞之间的初始韧带距离均为30 μm时，孔洞之间同时发生聚集，且孔洞之间聚集的起始时间约为1.57 μs，直径增长加速度约为2.484 Gm/s²；当孔洞2与孔洞3之间的韧带距离增加至50 μm时，孔洞1与孔洞2之间聚集的起始时间约为1.57 μs，直径增长加速度约为2.118 Gm/s²，孔洞3向孔洞2聚集的起始时间稍有延迟，约为1.58 μs，直径增长加速度约为1.168 Gm/s²；当孔洞2与孔洞3之间的韧带距离增加至70 μm时，也可观察到与50 μm时相似的情况。通过对比发现，孔洞之间聚集的起始时间随着初始韧带距离的增加而延长。孔洞2由于受两侧塑性应变场的影响，其直径增长速率与孔洞1和孔洞3存在明显区别，但是随着初始韧带距离的增加，孔洞2受孔洞3的影响变小，其直径增长速率越来越接近其他两个孔洞。

图 11 不同ILD下各孔洞直径增长速率时程曲线

Figure 11. Time history curves of void diameter growth rate with different ILD

下载: 全尺寸图片幻灯片

图12给出了不同时刻、不同初始韧带距离下3个孔洞聚集的演化云图。当t=1.54 μs时，损伤区域的拉伸应力达到孔洞成核应力（阈值），孔洞成核。当d_IL0=30 μm时，孔洞周围的塑性应变场几乎呈对称变化，3个孔洞之间同时发生贯穿，其中孔洞1和孔洞3的最终形状十分相似，而孔洞2则受两侧塑性应变场的影响，在水平方向上被拉长，且最终的形状也与孔洞1和孔洞3存在区别；随着孔洞2与孔洞3之间初始韧带距离的增加，可以看到孔洞1与孔洞2之间的塑性应变场优先相互作用，即初始韧带距离越短，孔洞越容易发生聚集。相同初始韧带距离下的孔洞1和孔洞2发生贯通的时间均为1.69 μs，而孔洞 2与孔洞3之间贯通的时间随着初始韧带距离的增加而延迟。

图 12 不同时刻不同初始韧带距离下三孔聚集演化云图

Figure 12. Evolution distribution of three voids coalescence in different ILD at different time

下载: 全尺寸图片幻灯片

3.2.2 孔洞直径的影响

通过改变成核孔洞的大小研究不同孔径对孔洞聚集时竞争机制的影响。对损伤区域内成核孔洞的直径进行调整，孔洞1和孔洞2的直径相等，设置为5 μm，孔洞3的直径d₃分别设置为5、10、15 μm，孔洞间初始韧带距离设置为50 μm，飞片的撞击速度设置为100 m/s。

图13给出了三孔聚集时各孔洞直径增长速率时程曲线，结果显示：在d₃为5、10、15 μm的条件下，孔洞2与孔洞1几乎同一时刻发生聚集，分别为1.58、1.57、1.57 μs，直径增长加速度分别约为1.916、2.022和2.887 Gm/s²；当d₃=5 μm时，孔洞2与孔洞3也几乎同一时刻发生聚集，随着d₃的增大，孔洞3的聚集时间为1.58 μs，稍有推迟，直径增长加速度约为0.695 Gm/s²。通过对比发现，孔洞2优先与孔洞1发生聚集，即孔洞优先与相同直径的孔洞聚集，不同直径孔洞之间发生聚集的时刻稍有延迟。

图 13 不同d₃下孔洞直径增长速率时程曲线

Figure 13. Time history curves of void diameter growth rate with different d₃

下载: 全尺寸图片幻灯片

图14给出了不同时刻、不同d₃下3个孔洞的长大、聚集演化云图。当t=1.54 μs时，各孔洞成核；随着时间的推移，当d₃=5 μm时，孔洞周围的塑性应变场几乎呈对称变化，并且3个孔洞同时贯通；当d₃逐渐增大时，孔洞1与孔洞2的塑性应变场优先交叠，并且优先发生贯通，孔洞1与孔洞3之间的塑性应变低于前者。图14显示的结果与图1(c)的实验结果相似。

图 14 不同时刻不同d₃条件下三孔聚集演化云图

Figure 14. Evolution distribution of three voids coalescence with different d₃ at different time

下载: 全尺寸图片幻灯片

4. 结　论

采用ABAQUS动力学有限元程序，对二维模型中双孔和三孔的聚集行为进行了数值模拟，讨论了孔洞间初始韧带距离、孔洞直径和孔洞位置分布的影响，得到了如下4点认识。

(1) 可以通过孔洞直径增长速率判断孔洞之间何时开始聚集，数值模拟结果显示，该方法具有一定的合理性。

(2) 相同孔径孔洞之间开始聚集的临界韧带距离随孔洞间初始韧带距离的增加而增大，相同孔径孔洞之间聚集时，孔洞各自均匀长大，不同孔径孔洞之间聚集时，小孔洞优先开始向大孔洞聚集。

(3) 相同孔洞条件下，在孔洞长大阶段，45°分布方式孔洞的增长加速度相比于0°和90°分布时孔洞的增长加速度更大，且45°分布方式的孔洞之间聚集的起始时间最早。

(4) 三孔聚集时，相同韧带距离下的相同孔洞之间几乎同时聚集且同时贯通，孔洞间聚集的起始时间随初始韧带距离的增加而延后；孔径不同时，大孔向附近小洞聚集的起始时间较迟，相同孔洞之间优先发生聚集并贯通，不同孔洞之间贯通较慢。

图 1 港池内水下爆炸环境示意图

Figure 1. Diagram of the environment of underwater explosion in harbor basin

下载: 全尺寸图片幻灯片

图 2 有限元计算模型(单位：cm)

Figure 2. FEM calculation model (Unit: cm)

下载: 全尺寸图片幻灯片

图 3 码头结构三视图和剖面图（单位：cm）

Figure 3. Three views and sectional views of wharf（Unit: cm）

下载: 全尺寸图片幻灯片

图 4 炸药及测点位置

Figure 4. Positions of explosive and measure point

下载: 全尺寸图片幻灯片

图 5 工况1气泡形态图

Figure 5. Configuration changing process of bubble in Condition 1

下载: 全尺寸图片幻灯片

图 9 工况5气泡形态图

Figure 9. Configuration changing process of bubble in Condition 5

下载: 全尺寸图片幻灯片

图 6 工况2气泡形态图

Figure 6. Configuration changing process of bubble in Condition 2

下载: 全尺寸图片幻灯片

图 7 工况3气泡形态图

Figure 7. Configuration changing process of bubble in Condition 3

下载: 全尺寸图片幻灯片

图 8 工况4气泡形态图

Figure 8. Configuration changing process of bubble in Condition 4

下载: 全尺寸图片幻灯片

图 10 荷载沿水深分布

Figure 10. Load distribution at different depth

下载: 全尺寸图片幻灯片

图 11 工况1中结构变形与毁伤过程

Figure 11. Process of the deformation and damage of structure in Condition 1

下载: 全尺寸图片幻灯片

图 12 工况3中结构变形与毁伤过程

Figure 12. Process of the deformation and damage of structure in Condition 3

下载: 全尺寸图片幻灯片

表 1 主要部位混凝土厚度及配筋情况

Table 1. Concrete thickness and matching bar condition of main members

Position	Concrete thickness/cm	Reinforcement situation	Cover thickness/cm
Cabin ex-wall	60	Double two way, $\varnothing$ 2.2 cm，@ 60 cm	20
Breast wall	60	Double two way, $\varnothing$ 2.2 cm，@ 60 cm	20
Partition	60	No reinforcement
Cabin floor	60	No reinforcement
Sealed plate	30	No reinforcement
Face plate	30	No reinforcement

下载: 导出CSV

表 2 材料参数^[11-12]

Table 2. Material parameters^[11-12]

Material	$\rho$ _a0/(kg·m^–3)	E_a/(MJ·kg^–1)	C₀	C₁	C₂	C₃	C₄	C₅	C₆
Air	1.293	0.25	0	0	0	0	0.4	0.4	0

Material	$\rho$ _w0/(kg·m^–3)	C/(m·s^–1)	S₁	S₂	S₃	$\gamma$ ₀
Water	1000	1480	2.56	–1.986	1.2268	0.5

Material	$\rho$ _e0/(kg·m^–3)	A/GPa	B/GPa	$\omega$	R₁	R₂	D/(m·s^–1)	p_CJ/GPa
Explosive	1654	374	3.23	0.3	4.15	0.95	6390	27

Material	$\rho$ _s0/(kg·m^–3)	E_s/MPa	G_s/MPa
Soil	1860	22.4	8

下载: 导出CSV

表 3 填土参数^[14]

Table 3. Parameters of backfill^[14]

$\rho$ ₀/(kg·m^–3)	E/MPa	G/MPa	Yield stress/MPa	Cutoff pressure/MPa	Failure strain
1800	47.38	16.01	7.70	–0.70	1.2

下载: 导出CSV

表 4 钢筋材料参数

Table 4. Parameters of steel bar

Density/ （g·cm^–3）	Poisson’s ratio	Initial yield stress/MPa	Elastic modulus/GPa	Tangent modulus/GPa	Strain rate/s^–1	Strain rate parameter	Failure strain	Reinforcement parameter
7.85	0.3	335	210	1.2	40	5	0.12	0

下载: 导出CSV

表 5 C40混凝土HJC模型参数

Table 5. Parameters of C40 concrete used in HJC model

$\rho$ /(kg·m^–3)	G/GPa	F_c'/MPa	A	B	C	N
2440	11.01	31.60	0.79	1.6	0.007	0.61

S_max	D₁	D₂	EFMIN	T/MPa	P_crush/MPa	$\mu$ _crush
7.0	0.036	1.0	0.0080	3.49	10.53	0.0013

P_lock/GPa	${\mu _{{\rm{lock}}}}$	k₁/GPa	k₂/GPa	k₃/GPa	${{\dot \varepsilon }_0}/{{\rm{s}}^{ - 1}}$	f_s
0.80	0.11	85	–171	208	1	0.004

下载: 导出CSV

表 6 C45混凝土HJC模型参数

Table 6. Parameters of C45 concrete used in HJC model

${\rho _0}$ /(kg·m^-3)	G/GPa	F_c'/MPa	A	B	C	N
2440	11.68	35.55	0.79	1.6	0.007	0.61

S_max	D₁	D₂	EFMIN	T/MPa	P_crush/MPa	${\mu _{\rm crush}}$
7.0	0.037	1.0	0.0085	3.70	11.85	0.0014

P_lock/GPa	${\mu _{{\rm{lock}}}}$	k₁/GPa	k₂/GPa	k₃/GPa	${{\dot \varepsilon }_0}/{{\rm{s}}^{ - 1}}$	f_s
0.80	0.11	85	–171	208	1	0.004

下载: 导出CSV

表 7 工况设置

Table 7. Simluation conditions

Condition	Explosive position	W/kg	R/m	H/m	D/m	$\overline R$ /(m·kg^–1/3)	$\overline H$ /(m·kg^–1/3)	$\overline D$ /(m·kg^–1/3)
1	C₁	100	2	3	0	0.43	0.65	0
2	C₂	100	4	3	0	0.86	0.65	0
3	C₃	100	6	3	0	1.29	0.65	0
4	C₄	100	8	3	0	1.72	0.65	0
5	C₅	100	10	3	0	2.16	0.65	0

下载: 导出CSV

表 8 气泡脉动规律

Table 8. Pattern of bubble impulse

Condition	${\bar R}$ /(m·kg^–1/3)	t_re/s	Maximum shape	D_re/m	t₁/s	Minimum shape	t₂/s	t_b/s
1	0.43	0.24	Semi pyriform cavity	7.8		Irregular discrete bubble		0.18
2	0.86	0.34	Semi pyriform cavity	10.0		Irregularly multiconnected domain		0.42
3	1.29	0.32	Offside platy pyriform cavity	11.4	0.78	Cyclic multiconnected domain	1.50	0.58
4	1.72	0.32	Vertical symmetry pyriform cavity	12.0	0.78	Cyclic multiconnected domain	1.38	1.48
5	2.16	0.32	Vertical symmetry pyriform cavity	12.2	0.78	Cyclic multiconnected domain	1.42	2.52

下载: 导出CSV

表 9 码头结构各部分吸收能量

Table 9. Energy absorption of different parts of the harbor basin

$\overline R$ /(m·kg^–1/3)	E_c/MJ	${\eta _{ {\rm{c} }} }$ /%	E_cs/MJ	${\eta _{{\rm{cs}}}}$ /%	E_s/MJ	${\eta _{{\rm{s}}}}$ /%	E_ts/MJ	${\eta _{{\rm{ts}}}}$ /%	E_sum/MJ
0.43	5.08	57.44	3.56	40.29	0.1306	1.48	0.0703	0.79	8.84
0.86	1.97	44.30	2.41	54.16	0.0069	0.15	0.0618	1.39	4.46
1.29	0.81	29.98	1.85	68.19	0.0014	0.05	0.0485	1.78	2.71
1.72	0.36	19.03	1.49	78.97	0.0003	0.01	0.0374	1.99	1.88
2.16	0.19	12.72	1.26	85.41	0.0002	0.02	0.0272	1.85	1.47

下载: 导出CSV

[1]	ZAMYSHLYAYEV B V. Dynamic loads in underwater explosion: AD-757183 [R]. Suitland: Naval Intelligence Support Center, 1972.
[2]	SWISDAK M M. Explosion effects and properties (Part II) —explosion effects in water: AD-A056694 [R]. Sliver Spring: Naval Surface Weapons Center, 1978.
[3]	LI J, RONG J L. Bubble and free surface dynamics in shallow underwater explosion [J]. Ocean Engineering, 2011, 38(17/18): 1861–1868.
[4]	CUI P, ZHANG A M, WANG S P. Small-charge underwater explosion bubble experiments under various boundary conditions [J]. Physics of Fluids, 2016, 28(11): 1–25.
[5]	WANG G H, ZHANG S R, YU M, et al. Investigation of the shock wave propagation characteristics and cavitation effects of underwater explosion near boundaries [J]. Applied Ocean Research, 2014, 46(2): 40–53.
[6]	张显丕, 刘建湖, 潘建强, 等. 基于效应靶的装药水下近场爆炸威力评估方法 [J]. 兵工学报, 2016, 37(8): 1430–1435. doi: 10.3969/j.issn.1000-1093.2016.08.013 ZHANG X P, LIU J H, PAN J Q, et al. An evaluation method for near-field underwater explosion power based on effect target [J]. Acta Armamentarii, 2016, 37(8): 1430–1435. doi: 10.3969/j.issn.1000-1093.2016.08.013
[7]	王振雄, 顾文彬, 陈江海, 等. 浅水中爆炸水底介质对水中冲击波峰值压力影响的试验研究 [J]. 振动与冲击, 2017, 36(4): 243–248. WANG Z X, GU W B, CHEN J H, et al. Experimental study on the influence of the bottom medium on the peak pressure of explosion shock waves in shallow water [J]. Journal of Vibration and Shock, 2017, 36(4): 243–248.
[8]	韦灼彬. 钢筋混凝土桩基梁板码头爆炸毁伤及抢修技术研究[D]. 天津: 天津大学, 2005.
[9]	侯晓峰, 王全胜, 钱展芃, 等. 水中爆炸条件下结构毁伤评估方法研究 [J]. 防护工程, 2014, 36(1): 40–45. HOU X F, WANG Q S, QIAN Z P, et al. Study on the assessment methods of structures subjected to underwater explosion [J]. Protective Engineering, 2014, 36(1): 40–45.
[10]	韦灼彬, 唐廷, 王立军. 港口水下爆炸荷载冲击特性研究 [J]. 振动与冲击, 2014, 33(6): 18–22. WEI Z B, TANG T, WANG L J. Shock characteristics of underwater explosion in port [J]. Journal of Vibration and Shock, 2014, 33(6): 18–22.
[11]	董琪, 韦灼彬, 唐廷, 等. 爆炸深度对浅水爆炸气泡脉动的影响 [J]. 高压物理学报, 2018, 32(2): 024102. doi: 10.11858/gywlxb.20170580 DONG Q, WEI Z B, TANG T, et al. Influence of explosion depth on bubble pulsation in shallow water explosion [J]. Chinese Journal of High Pressure Physics, 2018, 32(2): 024102. doi: 10.11858/gywlxb.20170580
[12]	董琪, 韦灼彬, 唐廷, 等. 浅水爆炸特性数值模拟研究 [J]. 海军工程大学学报, 2018, 30(3): 76–81. DONG Q, WEI Z B, TANG T, et al. Numerical simulation of characters of shallow water explosion [J]. Journal of Naval University of Engineering, 2018, 30(3): 76–81.
[13]	中国水电顾问集团西北勘测设计研究院. 水工混凝土结构设计规范: DLT 5057-2009 [S]. 北京: 中国建筑工业出版社, 2015.
[14]	高洪泉, 卢芳云, 赵宏伟. 不同土壤介质中爆炸的数值模拟[C]//第六届全国工程结构安全防护学术会议论文集. 洛阳, 2007: 85-90.
[15]	HOLMQUIST T J, JOHNSON G R. A computational constitutive model for glass subjected to large strains, high strain rates and high pressures [J]. Journal of Applied Mechanics, 2011, 78(5): 051003. doi: 10.1115/1.4004326
[16]	张凤国, 李恩征. 大应变、高应变率及高压强条件下混凝土的计算模型 [J]. 爆炸与冲击, 2002, 22(3): 198–202. doi: 10.3321/j.issn:1001-1455.2002.03.002 ZHANG F G, LI E Z. A computational model for concrete subjected to large strains, high strain rates, and high pressures [J]. Explosion and Shock Waves, 2002, 22(3): 198–202. doi: 10.3321/j.issn:1001-1455.2002.03.002