# The Fluid-Solid Interaction Dynamics between Underwater Explosion Bubble and Corrugated Sandwich Plate

Authors: Hao Wang; Yuan Sheng Cheng; Jun Liu; Lin Gan
Journal: Shock and Vibration
 (2016)

---

## Abstract

Lightweight sandwich structures with highly porous 2D cores or 3D (three-dimensional) periodic cores can effectively withstand underwater explosion load. In most of the previous studies of sandwich structure antiblast dynamics, the underwater explosion (UNDEX) bubble phase was neglected. As the UNDEX bubble load is one of the severest damage sources that may lead to structure large plastic deformation and crevasses failure, the failure mechanisms of sandwich structures might not be accurate if only shock wave is considered. In this paper, detailed 3D finite element (FE) numerical models of UNDEX bubble-LCSP (lightweight corrugated sandwich plates) interaction are developed by using MSC.Dytran. Upon the validated FE model, the bubble shape, impact pressure, and fluid field velocities for different stand-off distances are studied. Based on numerical results, the failure modes of LCSP and the whole damage process are obtained. It is demonstrated that the UNDEX bubble collapse jet local load plays a more significant role than the UNDEX shock wave load especially in near-field underwater explosion.


---

## Body

## 1. Introduction

Underwater blast is much more destructive than free field explosion in air. Since World War 1 and World War 2, a lot of naval vessels were attacked and destroyed by underwater explosion weapons [1]. On March 26, 2010, the Cheonan Naval Ship (PCC-772) was attacked and sank into West sea of Korea. The whole ship was cut into two pieces (Figure 1). According to the official Joint Investigation Report (JIR), this attack was underwater explosion caused by a torpedo with 300 kg TNT at a depth of about 6~9 m [2]. Thus, the UNDEX-resistant structures that can withstand these extreme loading conditions and sustain their functionality are critical in modern military setting. Due to structural efficiency and high energy absorption capability, sandwich systems have been extensively used in a variety of applications for many years. They may be the potential structure types of naval vessel hulls. The failure mechanism of this type structure subject to UNDEX is quite different from the ones of traditional marine constructions of stiffened plates.Figure 1
Example of Cheonan Naval Ship damage after close-in UNDEX [2].To develop and test innovative lightweight structural concepts for the US Navy, Wiernicki et al. [3] firstly investigated the elastic and plastic dynamic behaviors of LCSP subjected to air blast loading. A set of relatively simple closed analytical expressions were also given to quickly identify the effects of geometric parameters on dynamic behaviors. The optimization problem of LCSP under blast load was firstly discussed by Liang et al. [4]. The Feasible Direction Method (FDM) coupled with the Backtrack Program Method (BPM) was used in the optimum design algorithm, and the main structural parameters including corrugation angle, face sheet thickness, core thickness, and corrugation pitch were selected as design variables in optimum mathematical model. The calculations by Wiernicki et al. [3] and Liang et al. [4] are based on the semianalytical empirical formulations which cannot capture the details of air blast phenomenon.By adopting the Taylor plate assumptions, the fluid-solid interaction (FSI) analysis of LCSP under UNDEX shock wave loading was firstly considered in the investigation proposed by Xue and Hutchinson [5]. If the blast medium is water, it was found that FSI effecting can reduce the momentum imparted to a sandwich plate by almost a factor of two relative to that imparted to a solid plate of the same weight. Vaziri and Hutchinson [6] further complemented the previous studies on the role of FSI by accounting for the nonlinear compressibility and finite shock behavior of air medium. The results showed that the FSI enhances the performance of LCSP relative to monolithic plates under intense air explosion, but not as significantly as for UNDEX. In practical cases, if the impulse of impact loading is sufficiently large, the damage of LCSP may occur. Based on the nonlinear finite element computations by adopting ABAQUS/Explicit code, the failure mode maps of LCSP under intense uniform impulsive pressure loads were obtained from the investigation of Vaziri et al. [7]. Here, it must be pointed out that the FSI is not considered in the numerical simulations proposed by Vaziri et al. [7] (if the effects of FSI are included, the failure mechanisms of LCSP may be different). To investigate structural designs of vessels against collisions, Rubino et al. [8] measured and analyzed the dynamic performance of sandwich beams with the Y-frame and corrugated cores. Both the experimental and finite element (FE) results revealed that these two topology type cores sandwich beams with equal mass have similar dynamic behavior. Qin et al. [9] studied low velocity impact resistance of a LCSP struck by a heavy mass. In addition, new analytical predictions were obtained and the predicted results agreed well with FE results. Rimoli et al. [10] utilized experimental tests and numerical methods to investigate the dynamic responses of edge-clamped LCSP subject to the shock of explosively driven wet sand. In the modeling, a decoupled wet sand loading curve was incorporated into FE simulation.Recently, Wadley et al. [11] further investigated this sand-water-structure coupling phenomenon in mine explosion accident. The dynamic deformation and fracture processes were both included in this analysis by employing a particle-based method. The high speed fragment penetration problem of aluminum alloy LCSP with empty and alumina filled core has been done experimentally by Wadley et al. [12]. Zhang et al. [13] also conducted a study on the dynamic response of LCSP with unfilled and foam-filled sinusoidal plate cores. A novel analytical procedure was built to evaluate the dynamic response of LCSP. Very few experimental investigations have been conducted on the blast resistance of corrugated sandwich panels. Li et al. [14] experimentally investigated failure mechanisms of two configurations of the LCSP specimen by using a ballistic pendulum system. The deflection modes demonstrated that only global deformation and small tearing crack occurred. In the previous investigation of our group, a three-dimensional fully coupled simulation is conducted to analyze the dynamic response of sandwich panels comprising equal thicknesses face sheets sandwiching a corrugated core when subjected to localized impulse created by the detonation of cylindrical explosive by Zhang et al. [15]. The numerical simulation results showed that the core configuration has a negligible influence on the peak reflected pressure, but it has an effect on the deflection of a panel.Most of these previous investigations are concerned with the “pure” shock wave impact loading, while the multidamage sources of load environment such as bubble are often neglected. Regarding the knowledge of authors, the detailed failure mechanisms of LCSP subjected to close-in UNDEX bubble loading are not very clear. Thus, the present paper is primarily concerned with the detailed dynamic behavior and the failure mechanism of LCSP subjected to close-in UNDEX shock wave load and following bubble pulse. An outline of this paper is as follows. In the initial part, the geometry characteristics of LCSP and the theoretical background of UNDEX bubble are introduced. To verify the numerical simulation model, the 3D near free surface UNDEX bubble problem is firstly analyzed and discussed by adopting multimaterial Euler-Lagrange coupling method in MSC.Dytran code. The detailed integrated response and actual deformation characteristics of LCSP are presented for three different cases. The bubble shape, fluid field velocity, and bubble collapse jet are investigated. In addition, some new remarks about failure mode of LCSP are discussed in the last part. The objective of this paper is to study the nonlinear inelastic responses and damage of LCSP under UNDEX shock wave and bubble loading and to develop more accurate predictions of bubble-LCSP interaction behavior.


## 2. Theoretical Background of UNDEX Bubble

### 2.1. Migration of Explosion Bubble

An early investigation of considerable significance in the field of UNDEX bubble research was conducted by Lamb [17] for analyzing the collapse of a spherical transient bubble in an infinite fluid. In that study, it was assumed that the pulsation pressure within the explosion bubble varies as(1)
p
=
p
0
V
0
V
γ
,where V is volume of bubble, γ is a constant, and the subscript denotes initial values. Lamb carried through the analysis for the cases γ
=
1 and γ
=
4
/
3.Based on the work of Lamb, including the buoyancy forces, the influence of viscous damping, and the ocean surface interaction, the complex relation of bubble is given by Vernon [18]: (2a)
x
˙
=
σ
,
(2b)
ζ
˙
=
λ
,
(2c)
σ
˙
=
-
3
δ
2
δ
-
β
x
σ
2
x
1
-
2
β
x
3
δ
-
λ
2
6
x
+
ζ
x
ζ
0
-
γ
-
1
k
x
3
γ
+
1
+
β
x
4
δ
2
C
d
λ
2
4
x
+
σ
λ
3
-
x
ζ
0
,
(2d)
λ
˙
=
-
3
α
1
ζ
0
+
σ
λ
x
-
C
d
λ
2
4
x
+
β
x
4
δ
2
3
σ
2
+
x
σ
˙
,where (3)
x
=
a
L
;
x
0
=
a
0
L
;
ζ
=
Z
L
;
ζ
0
=
Z
0
L
;
δ
=
d
L
;
τ
=
t
T
;
x
˙
=
d
x
d
τ
;
ζ
˙
=
d
ζ
d
τ
;
σ
˙
=
d
σ
d
τ
;
λ
˙
=
d
λ
d
τ
;and z is the water pressure head and the subscript denotes initial values,d is the depth of charge at initial time, γ is the adiabatic gas constant being equal to 1.25, and α is the migration control coefficient which can be set to the value of 1.0 if the migration is considered.  β is a free surface effect control coefficient which will always be equal to 0 or 1. β
=
0 means that free surface effect is not considered. Here, α
=
1 and β
=
1 are adopted in the analysis. The other three unknown parameters in (2a)–(2d) are the length scale factorL, time scale factor T, and the nondimensional energy factor k which can be defined as following relations [18]:(4)
L
=
3
E
0
4
π
ρ
g
Z
0
1
/
3
,
T
=
3
2
g
Z
0
1
/
2
×
L
,
k
=
ρ
g
Z
0
γ
-
1
γ
-
1
k
1
W
E
0
γ
-
1
,where E
0 is the total energy of explosion, ρ is the density of fluid, g is acceleration due to gravity,W is the charge weight, and k
1 is the parameter based on the charge type. For the nondimensional energy charge parameter k
1 for TNT, the simplified expression is specified:(5)
k
1
≅
0.0743
Z
0
0.25
.Based on the previous bubble dynamics equations (2a)–(5), the characteristics of bubble radius velocity and bubble vertical velocity can be solved using the fourth-order Runge-Kutta method once the initial conditions are given. And the initial conditions of UNDEX bubble can be defined as x
=
x
0, σ
=
σ
0, ζ
=
ζ
0, and λ
=
λ
0 and the initial radius of bubble R
0 can be obtained by solving the following energy conservation equation [19]: (6)
1.39
×
1
0
5
Δ
P
1
γ
-
1
3
W
4
π
R
m
3
1
-
R
0
R
m
-
3
γ
-
1
=
R
0
R
m
3
-
1
,where Δ
P
=
P
∞
-
P
c is the condensing steam pressure of water and P
∞ is the static water pressure at infinite distance. R
m (in m) is the first maximum radius of UNDEX bubble and is presented as follows [20]:(7)
R
m
=
K
1
W
1
/
3
d
+
10.34
1
/
3
.The first period of the bubble pulse wave T
1 (in s) can be also expressed using the empirical formula:(8)
T
1
=
K
2
W
1
/
3
d
+
10.34
5
/
6
,where K
1 and K
2 are the constants depending on explosive charge types (for TNT charge, K
1
=
3.5, K
2
=
2.11), W is the mass of the charge in kilograms, and d is the depth of charge in meters.


### 2.2. Kelvin Impulse and Blake Criteria

The Kelvin impulse is a particularly valuable concept in unsteady fluid dynamics. Benjamin and Ellis [21] seem to be the first to have realized its value in bubble dynamics. The Kelvin impulse may be used to determine aspects of the gross bubble motion, and it is defined as follows [22]:(9)
I
=
ρ
∫
S
ϕ
n
d
S
=
∫
0
t
F
t
d
t
,
(10)
F
t
=
ρ
g
V
e
+
ρ
∫
∑
b
1
2
∇
ϕ
2
n
-
∂
ϕ
∂
n
∇
ϕ
d
S
,where ρ is the fluid density, V is the bubble volume, and e is the unit vector of buoyancy force direction. ϕ is the velocity potential, S is surface of the bubble, ∑ is the boundary, and n is the outward normal to the fluid.The first part in (10) indicates the buoyancy force of bubble, and the second part in (10) is the Bjerknes force caused by the variation of fluid field near the boundary. Thus, the sum of these two forces leads to the bubbles migrating. For rigid boundary condition, the second part in (10) can be also written as [22](11)
ρ
∫
∑
b
1
2
∇
ϕ
2
n
-
∂
ϕ
∂
n
∇
ϕ
d
S
=
-
ρ
m
2
t
16
π
h
2
t
.According to the Rayleigh spherical bubble model [23] (which can be simplified from (2a)–(2d)), one may have (12)
R
d
2
R
d
t
2
+
3
2
d
R
d
t
+
Δ
P
ρ
=
0
.And suppose that h
(
t
) is constant throughout the motion and equals its initial value h
0. So one may obtain(13)
m
t
=
4
π
R
2
R
˙
=
±
4
π
R
2
2
3
Δ
P
ρ
R
m
3
R
3
-
1
1
/
2
.By using these results and integrating (9) over the lifetime of bubble, the Kelvin impulse at the end of collapse is obtained [22]:(14)
I
=
2
6
π
R
m
5
ρ
Δ
P
1
/
2
9
h
0
2
2
ε
2
δ
2
B
11
6
,
1
2
-
B
7
6
,
3
2
,where B
(
x
,
y
) is beta function, and ε and δ are defined as follows:(15)
ε
=
h
R
m
,
δ
=
ρ
g
R
m
Δ
P
1
/
2
.Blake and Cerone [22] pointed out that if I
>
0 the bubble will migrate away from the rigid boundary. For I
<
0, the bubble will migrate towards the rigid boundary and for I
=
0 the following relation between ε and δ can be observed: (16)
ε
δ
=
B
7
/
6,3
/
2
2
B
11
/
6,1
/
2
1
/
2
=
0.442
.This is well known as Blake criteria [22].


## 2.1. Migration of Explosion Bubble

An early investigation of considerable significance in the field of UNDEX bubble research was conducted by Lamb [17] for analyzing the collapse of a spherical transient bubble in an infinite fluid. In that study, it was assumed that the pulsation pressure within the explosion bubble varies as(1)
p
=
p
0
V
0
V
γ
,where V is volume of bubble, γ is a constant, and the subscript denotes initial values. Lamb carried through the analysis for the cases γ
=
1 and γ
=
4
/
3.Based on the work of Lamb, including the buoyancy forces, the influence of viscous damping, and the ocean surface interaction, the complex relation of bubble is given by Vernon [18]: (2a)
x
˙
=
σ
,
(2b)
ζ
˙
=
λ
,
(2c)
σ
˙
=
-
3
δ
2
δ
-
β
x
σ
2
x
1
-
2
β
x
3
δ
-
λ
2
6
x
+
ζ
x
ζ
0
-
γ
-
1
k
x
3
γ
+
1
+
β
x
4
δ
2
C
d
λ
2
4
x
+
σ
λ
3
-
x
ζ
0
,
(2d)
λ
˙
=
-
3
α
1
ζ
0
+
σ
λ
x
-
C
d
λ
2
4
x
+
β
x
4
δ
2
3
σ
2
+
x
σ
˙
,where (3)
x
=
a
L
;
x
0
=
a
0
L
;
ζ
=
Z
L
;
ζ
0
=
Z
0
L
;
δ
=
d
L
;
τ
=
t
T
;
x
˙
=
d
x
d
τ
;
ζ
˙
=
d
ζ
d
τ
;
σ
˙
=
d
σ
d
τ
;
λ
˙
=
d
λ
d
τ
;and z is the water pressure head and the subscript denotes initial values,d is the depth of charge at initial time, γ is the adiabatic gas constant being equal to 1.25, and α is the migration control coefficient which can be set to the value of 1.0 if the migration is considered.  β is a free surface effect control coefficient which will always be equal to 0 or 1. β
=
0 means that free surface effect is not considered. Here, α
=
1 and β
=
1 are adopted in the analysis. The other three unknown parameters in (2a)–(2d) are the length scale factorL, time scale factor T, and the nondimensional energy factor k which can be defined as following relations [18]:(4)
L
=
3
E
0
4
π
ρ
g
Z
0
1
/
3
,
T
=
3
2
g
Z
0
1
/
2
×
L
,
k
=
ρ
g
Z
0
γ
-
1
γ
-
1
k
1
W
E
0
γ
-
1
,where E
0 is the total energy of explosion, ρ is the density of fluid, g is acceleration due to gravity,W is the charge weight, and k
1 is the parameter based on the charge type. For the nondimensional energy charge parameter k
1 for TNT, the simplified expression is specified:(5)
k
1
≅
0.0743
Z
0
0.25
.Based on the previous bubble dynamics equations (2a)–(5), the characteristics of bubble radius velocity and bubble vertical velocity can be solved using the fourth-order Runge-Kutta method once the initial conditions are given. And the initial conditions of UNDEX bubble can be defined as x
=
x
0, σ
=
σ
0, ζ
=
ζ
0, and λ
=
λ
0 and the initial radius of bubble R
0 can be obtained by solving the following energy conservation equation [19]: (6)
1.39
×
1
0
5
Δ
P
1
γ
-
1
3
W
4
π
R
m
3
1
-
R
0
R
m
-
3
γ
-
1
=
R
0
R
m
3
-
1
,where Δ
P
=
P
∞
-
P
c is the condensing steam pressure of water and P
∞ is the static water pressure at infinite distance. R
m (in m) is the first maximum radius of UNDEX bubble and is presented as follows [20]:(7)
R
m
=
K
1
W
1
/
3
d
+
10.34
1
/
3
.The first period of the bubble pulse wave T
1 (in s) can be also expressed using the empirical formula:(8)
T
1
=
K
2
W
1
/
3
d
+
10.34
5
/
6
,where K
1 and K
2 are the constants depending on explosive charge types (for TNT charge, K
1
=
3.5, K
2
=
2.11), W is the mass of the charge in kilograms, and d is the depth of charge in meters.


## 2.2. Kelvin Impulse and Blake Criteria

The Kelvin impulse is a particularly valuable concept in unsteady fluid dynamics. Benjamin and Ellis [21] seem to be the first to have realized its value in bubble dynamics. The Kelvin impulse may be used to determine aspects of the gross bubble motion, and it is defined as follows [22]:(9)
I
=
ρ
∫
S
ϕ
n
d
S
=
∫
0
t
F
t
d
t
,
(10)
F
t
=
ρ
g
V
e
+
ρ
∫
∑
b
1
2
∇
ϕ
2
n
-
∂
ϕ
∂
n
∇
ϕ
d
S
,where ρ is the fluid density, V is the bubble volume, and e is the unit vector of buoyancy force direction. ϕ is the velocity potential, S is surface of the bubble, ∑ is the boundary, and n is the outward normal to the fluid.The first part in (10) indicates the buoyancy force of bubble, and the second part in (10) is the Bjerknes force caused by the variation of fluid field near the boundary. Thus, the sum of these two forces leads to the bubbles migrating. For rigid boundary condition, the second part in (10) can be also written as [22](11)
ρ
∫
∑
b
1
2
∇
ϕ
2
n
-
∂
ϕ
∂
n
∇
ϕ
d
S
=
-
ρ
m
2
t
16
π
h
2
t
.According to the Rayleigh spherical bubble model [23] (which can be simplified from (2a)–(2d)), one may have (12)
R
d
2
R
d
t
2
+
3
2
d
R
d
t
+
Δ
P
ρ
=
0
.And suppose that h
(
t
) is constant throughout the motion and equals its initial value h
0. So one may obtain(13)
m
t
=
4
π
R
2
R
˙
=
±
4
π
R
2
2
3
Δ
P
ρ
R
m
3
R
3
-
1
1
/
2
.By using these results and integrating (9) over the lifetime of bubble, the Kelvin impulse at the end of collapse is obtained [22]:(14)
I
=
2
6
π
R
m
5
ρ
Δ
P
1
/
2
9
h
0
2
2
ε
2
δ
2
B
11
6
,
1
2
-
B
7
6
,
3
2
,where B
(
x
,
y
) is beta function, and ε and δ are defined as follows:(15)
ε
=
h
R
m
,
δ
=
ρ
g
R
m
Δ
P
1
/
2
.Blake and Cerone [22] pointed out that if I
>
0 the bubble will migrate away from the rigid boundary. For I
<
0, the bubble will migrate towards the rigid boundary and for I
=
0 the following relation between ε and δ can be observed: (16)
ε
δ
=
B
7
/
6,3
/
2
2
B
11
/
6,1
/
2
1
/
2
=
0.442
.This is well known as Blake criteria [22].


## 3. Numerical Modeling and Simulation

Considering the practical application of LCSP in ship building industry, the LCSP with an exposed areaa
∗
b
=
1200.0
m
m
∗
1000.0
m
m (the total mass of LCSP M
s = 67.6 kg) and its geometric description are plotted in Figure 2. The geometric parameters of LCSP are φ
=
π
/
4, d
c = 120.0 mm, t
f = 2.5 mm, t
b = 2.5 mm, t
c = 1.5 mm, and H
c = 38.250 mm. The calculated parameters for all 3 cases are shown in Table 1.Table 1
Various computational cases considered.
Case number
a (mm)
b (mm)
W (kg)
d (depth of charge, m)
T
1 (s)
R
m
a
x (m)
d/R
m
a
x
1
1200.00
1000.00
0.05
0.6
0.106
0.581
~1.0
2
1200.00
1000.00
0.05
0.3
0.108
0.586
~0.5
3
1200.00
1000.00
5.00
0.6
0.491
2.695
~0.2Figure 2
Schematic diagram of a lightweight corrugated sandwich plate (LCSP).


### 3.1. Lagrange Finite Element Model

The face sheets and corrugated core are modeled as a plane plate using 52000 quadrilateral shell elements (CQUAD4, KEYHOFF formulation, hourglass control, five degrees of freedom per nodeu
x, u
y, u
z, θ
x, θ
y, and finite membrane strains elements, with 5 integration points). The detailed Lagrange FE model of LCSP can be seen in Figure 3. In the simulation, the Lagrange material of the LCSP is modeled to be Q235 steel. In order to consider the strain rate effect, the Cowper-Symonds model [24] is adopted in the analysis. And the material property constants are described in previous investigation [24]. Here, the FAILMP Sentry was used to represent the element failure model. The value of 0.24 was used from MSC/MVISION database according to the test results [25].Figure 3
Lagrange FE model of LCSP and location of charge.


### 3.2. Euler Finite Element Model

Boundary integral method (BIM) was widely used in the bubble dynamics early in the 1960s, which has been validated by many experiments. Gong and Khoo [26] analyzed the transient response of stiffened composite submersible hull subjected to underwater explosion bubble by adopting the coupled BEM-FEM method to handle the interaction of the glass-epoxy composite structures and the underwater explosion bubble. And the effect of bubble locations on the composite submersible hull is also studied and analyzed. By using boundary integral method (BIM) and multiple vortex rings model, Zhang et al. found some new phenomena such as more splits after the first split of the toroidal bubble in the splitting of a toroidal bubble near a rigid boundary according to the numerical simulation and two experiments [27]. Based on the vortex ring for arbitrary location in 3D model and a new density potential method (DMP), higher accuracy and stability are obtained to capture detailed features of bubble deformation especially for the large deforming problem and the toroidal bubble phase [28]. Though boundary integral method (BIM) is a traditional technique for underwater explosion bubble problem, some very complex phenomena such as the water splash (especially when the free surface is considered), the fracture of structures, and mesh distortion are hard to overcome. Moreover, some remeshing techniques introduced to decrease the error due to mesh distortion are studied in [27, 28]. However, the artificial numerical noise cannot be easily avoided during the whole simulation process. And these complicated numerical algorithms will increase the complexity of algorithm and the amount of calculation. Some other different numerical techniques, for example, finite volume method (FVM) and smoothed particle hydrodynamics (SPH), are adopted to deal with underwater explosion bubble problem. Thus, the CFD solver of MSC.Dytran by using an Eulerian approach and a finite volume method is adopted in the present simulation. And the fluid governing equations are the conservation laws and are integrated in time by a first-order explicit dynamic procedure [25].To model the fluid inside and outside the LCSP, two Euler domains are used. The outer domain has the LCSP surface (including top face sheet, bottom face sheet, and out-off rigid wall) as part of the fluid boundary. Euler material is outside the LCSP surface and there is no material inside the LCSP surface. The contents inside the LCSP are modeled in the inner domain and this domain is also enclosed by the LCSP surface. Therefore, both Euler domains use the LCSP surface as part of their enclosure (Figures4 and 5). The outer Euler domain and its enclosing surface are shown in Figure 4.Figure 4
The outer Euler domain.Figure 5
The inner Euler domain of LCSP.The outer boundary of the outer domain is given by a sufficiently large fixed box. Pressure at the outer boundary is set to the hydrostatic pressure by using HYDRSTAT and FLOWDEF keyword cards. This behaves as the open boundary. The Euler mesh contains the water and the air on the top of the water. The fluid mesh used for this problem consists of a block of elements, with the dimensions2.5
m
∗
2.5
m
∗
4.0
m. This fluid block of water and air was meshed with 100
∗
100
∗
60 hexahedron elements, and the total number of fluid elements is 900000 (Figure 6). The gird size of inner Euler domain is 0.02 m in this simulation. All boundary conditions for the fluid mesh shown in Figures 4 and 5 were given a “flow” boundary condition by adopting TICEUL and TICVAL keyword cards.Figure 6
Euler domain mesh.In order to model the fresh water, a polynomial equation of state was conducted. This state equation (EOS) of fresh water relates the pressure in the fluid to the acoustic condensationμ and the specific internal energy by(17a)
P
=
a
1
μ
+
a
2
μ
2
+
a
3
μ
3
+
b
0
+
b
1
μ
+
b
2
μ
2
ρ
0
E
μ
>
0
,
(17b)
P
=
a
1
μ
+
b
0
+
b
1
μ
2
ρ
0
E
μ
<
0
,where μ
=
(
ρ
-
ρ
0
)
/
ρ
0, ρ
0 is the initial density of fresh water. E is the specific internal energy per unit mass, and a
1, a
2, a
3, b
0, b
1, and b
2 are the constants of the fluid, respectively. And (17a) applies to a fluid in a compressed state, while (17b) applies to a fluid in an expanded state. And the constants in (17a) and (17b) are provided in Table 2 [29].Table 2
Material properties (the fresh water and the air).
Material
MSC.Dytran model
Input parameters
Water
Polynomial equation of state
a
1 = 2.314 × 109 Pa, a
2 = 6.561 × 109 Pa, a
3 = 1.126 × 1010 Pa, b
0 = 0.4934, b
1 = 1.3937,
b
2 = 0.0000, ρ
0 = 1025 kg/m3, and E = 3750.4 J/kg
Air
Gamma law
γ = 1.40, ρ = 1.185 kg/m3, and E = 287 J/kgThe gamma law gas model is adopted for the EOS of air:(18)
P
=
γ
-
1
ρ
E
,where ρ is the density of air,  γ is the heat capacities of the gas, andE is the specific internal energy of air. The initial pressure of air is set to 1.0 × 105 Pa.The TNT explosive can be modeled by a JWL EOS in MSC.Dytran. However, if the explosive is a spherical ball, the radius of this ball is only 0.04 m when the mass of TNT charge is set to 0.5 kg. A finer gird has to be modeled to simulate this small ball. In this analysis, the TNT explosive is defined as a compressed hot gas (γ = 1.25; see (2a)–(2d)). The mass and the specific internal energy are those of the TNT charge. The radius of this hot gas ball is calculated using (6) and the density of air is adjusted based on the equivalent mass of the TNT explosive.


### 3.3. Coupled Fluid-Structure Interaction

Multimaterial Euler solver in MSC.Dytran allows for up to 9 different Eulerian materials to be presented in a given investigation. The two different models (General Coupling method and Arbitrary Lagrange-Euler method) are available to calculate the FSI problem between Eulerian and Lagrangian materials. The detailed descriptions of these two FSI algorithms can be seen in MSC.Dytran user’s manual [25].In the present analysis, the “General Coupling” algorithm is used. In this algorithm, the Lagrangian and Eulerian meshes are geometrically independent and interact via closed coupling surface attached to LCSP (also see Figures4 and 5). The deformation coupling surface “cuts across” Eulerian elements which contain multimaterial including air and water, changing their volume and surface areas. As the LCSP FE mesh deforms under the action of the impact pressure from the Eulerian mesh, the resulting FE deflection then influences subsequent material flow and pressure forces in the Eulerian mesh, resulting in automatic and precise coupling of FSI (Figure 7).Figure 7
General Coupling function and LCSP-UNDEX bubble simulation model.When the pressure force of bubble is sufficiently large, the breach damage may occur in the simulation. Under this condition, the flow transfer of fluid from the breach between outer Euler domain and inner Euler domain will be important. To simulate the flow transfer of breach, the adaptive multiple Euler domains technology and failure of coupling surface technology are adopted in the simulation. The PARAM, FASTCOUP, INPLANE, and FAIL keyword cards are used to consider failure of the coupling surface, and the PORFLCPL and COUPOR keyword cards are adopted here to model the transfer of different Euler domains.


## 3.1. Lagrange Finite Element Model

The face sheets and corrugated core are modeled as a plane plate using 52000 quadrilateral shell elements (CQUAD4, KEYHOFF formulation, hourglass control, five degrees of freedom per nodeu
x, u
y, u
z, θ
x, θ
y, and finite membrane strains elements, with 5 integration points). The detailed Lagrange FE model of LCSP can be seen in Figure 3. In the simulation, the Lagrange material of the LCSP is modeled to be Q235 steel. In order to consider the strain rate effect, the Cowper-Symonds model [24] is adopted in the analysis. And the material property constants are described in previous investigation [24]. Here, the FAILMP Sentry was used to represent the element failure model. The value of 0.24 was used from MSC/MVISION database according to the test results [25].Figure 3
Lagrange FE model of LCSP and location of charge.


## 3.2. Euler Finite Element Model

Boundary integral method (BIM) was widely used in the bubble dynamics early in the 1960s, which has been validated by many experiments. Gong and Khoo [26] analyzed the transient response of stiffened composite submersible hull subjected to underwater explosion bubble by adopting the coupled BEM-FEM method to handle the interaction of the glass-epoxy composite structures and the underwater explosion bubble. And the effect of bubble locations on the composite submersible hull is also studied and analyzed. By using boundary integral method (BIM) and multiple vortex rings model, Zhang et al. found some new phenomena such as more splits after the first split of the toroidal bubble in the splitting of a toroidal bubble near a rigid boundary according to the numerical simulation and two experiments [27]. Based on the vortex ring for arbitrary location in 3D model and a new density potential method (DMP), higher accuracy and stability are obtained to capture detailed features of bubble deformation especially for the large deforming problem and the toroidal bubble phase [28]. Though boundary integral method (BIM) is a traditional technique for underwater explosion bubble problem, some very complex phenomena such as the water splash (especially when the free surface is considered), the fracture of structures, and mesh distortion are hard to overcome. Moreover, some remeshing techniques introduced to decrease the error due to mesh distortion are studied in [27, 28]. However, the artificial numerical noise cannot be easily avoided during the whole simulation process. And these complicated numerical algorithms will increase the complexity of algorithm and the amount of calculation. Some other different numerical techniques, for example, finite volume method (FVM) and smoothed particle hydrodynamics (SPH), are adopted to deal with underwater explosion bubble problem. Thus, the CFD solver of MSC.Dytran by using an Eulerian approach and a finite volume method is adopted in the present simulation. And the fluid governing equations are the conservation laws and are integrated in time by a first-order explicit dynamic procedure [25].To model the fluid inside and outside the LCSP, two Euler domains are used. The outer domain has the LCSP surface (including top face sheet, bottom face sheet, and out-off rigid wall) as part of the fluid boundary. Euler material is outside the LCSP surface and there is no material inside the LCSP surface. The contents inside the LCSP are modeled in the inner domain and this domain is also enclosed by the LCSP surface. Therefore, both Euler domains use the LCSP surface as part of their enclosure (Figures4 and 5). The outer Euler domain and its enclosing surface are shown in Figure 4.Figure 4
The outer Euler domain.Figure 5
The inner Euler domain of LCSP.The outer boundary of the outer domain is given by a sufficiently large fixed box. Pressure at the outer boundary is set to the hydrostatic pressure by using HYDRSTAT and FLOWDEF keyword cards. This behaves as the open boundary. The Euler mesh contains the water and the air on the top of the water. The fluid mesh used for this problem consists of a block of elements, with the dimensions2.5
m
∗
2.5
m
∗
4.0
m. This fluid block of water and air was meshed with 100
∗
100
∗
60 hexahedron elements, and the total number of fluid elements is 900000 (Figure 6). The gird size of inner Euler domain is 0.02 m in this simulation. All boundary conditions for the fluid mesh shown in Figures 4 and 5 were given a “flow” boundary condition by adopting TICEUL and TICVAL keyword cards.Figure 6
Euler domain mesh.In order to model the fresh water, a polynomial equation of state was conducted. This state equation (EOS) of fresh water relates the pressure in the fluid to the acoustic condensationμ and the specific internal energy by(17a)
P
=
a
1
μ
+
a
2
μ
2
+
a
3
μ
3
+
b
0
+
b
1
μ
+
b
2
μ
2
ρ
0
E
μ
>
0
,
(17b)
P
=
a
1
μ
+
b
0
+
b
1
μ
2
ρ
0
E
μ
<
0
,where μ
=
(
ρ
-
ρ
0
)
/
ρ
0, ρ
0 is the initial density of fresh water. E is the specific internal energy per unit mass, and a
1, a
2, a
3, b
0, b
1, and b
2 are the constants of the fluid, respectively. And (17a) applies to a fluid in a compressed state, while (17b) applies to a fluid in an expanded state. And the constants in (17a) and (17b) are provided in Table 2 [29].Table 2
Material properties (the fresh water and the air).
Material
MSC.Dytran model
Input parameters
Water
Polynomial equation of state
a
1 = 2.314 × 109 Pa, a
2 = 6.561 × 109 Pa, a
3 = 1.126 × 1010 Pa, b
0 = 0.4934, b
1 = 1.3937,
b
2 = 0.0000, ρ
0 = 1025 kg/m3, and E = 3750.4 J/kg
Air
Gamma law
γ = 1.40, ρ = 1.185 kg/m3, and E = 287 J/kgThe gamma law gas model is adopted for the EOS of air:(18)
P
=
γ
-
1
ρ
E
,where ρ is the density of air,  γ is the heat capacities of the gas, andE is the specific internal energy of air. The initial pressure of air is set to 1.0 × 105 Pa.The TNT explosive can be modeled by a JWL EOS in MSC.Dytran. However, if the explosive is a spherical ball, the radius of this ball is only 0.04 m when the mass of TNT charge is set to 0.5 kg. A finer gird has to be modeled to simulate this small ball. In this analysis, the TNT explosive is defined as a compressed hot gas (γ = 1.25; see (2a)–(2d)). The mass and the specific internal energy are those of the TNT charge. The radius of this hot gas ball is calculated using (6) and the density of air is adjusted based on the equivalent mass of the TNT explosive.


## 3.3. Coupled Fluid-Structure Interaction

Multimaterial Euler solver in MSC.Dytran allows for up to 9 different Eulerian materials to be presented in a given investigation. The two different models (General Coupling method and Arbitrary Lagrange-Euler method) are available to calculate the FSI problem between Eulerian and Lagrangian materials. The detailed descriptions of these two FSI algorithms can be seen in MSC.Dytran user’s manual [25].In the present analysis, the “General Coupling” algorithm is used. In this algorithm, the Lagrangian and Eulerian meshes are geometrically independent and interact via closed coupling surface attached to LCSP (also see Figures4 and 5). The deformation coupling surface “cuts across” Eulerian elements which contain multimaterial including air and water, changing their volume and surface areas. As the LCSP FE mesh deforms under the action of the impact pressure from the Eulerian mesh, the resulting FE deflection then influences subsequent material flow and pressure forces in the Eulerian mesh, resulting in automatic and precise coupling of FSI (Figure 7).Figure 7
General Coupling function and LCSP-UNDEX bubble simulation model.When the pressure force of bubble is sufficiently large, the breach damage may occur in the simulation. Under this condition, the flow transfer of fluid from the breach between outer Euler domain and inner Euler domain will be important. To simulate the flow transfer of breach, the adaptive multiple Euler domains technology and failure of coupling surface technology are adopted in the simulation. The PARAM, FASTCOUP, INPLANE, and FAIL keyword cards are used to consider failure of the coupling surface, and the PORFLCPL and COUPOR keyword cards are adopted here to model the transfer of different Euler domains.


## 4. Numerical Verification: 3D Near Free Surface UNDEX Bubble

To verify the reliability of the developed FE model, a three-dimensional (3D) near free surface UNDEX bubble model is first built for calibration purpose. And the computational parameters of validation case are as follows:W = 0.3 kg (mass of TNT charge), ρ
T
N
T = 1630 kg/m3 (density of TNT charge), γ = 1.25 (for bubble gas), γ = 1.40 (for air),d = 1.0 m, R
0 = 0.075 m (initial bubble radius, from (6)), ρ
0 = 168.976 kg/m3 (initial bubble density), P
0 = 84.68 MPa (initial bubble pressure), R
m
a
x = 1.04 m (maximum bubble radius), and T
1 = 0.151 s (first period of the bubble pulse).The entire volume of the bubble in the first bubble circulation, with use of the developed FE model as simulated by DYTRAN, is presented in Figure8. Figure 8(a) shows the initial conditions of the bubble, at which time the bubble rapidly expands outward with high internal pressure. Att = 0.07 s, the bubble expands to its maximum size. Att = 0.09 s, it begins to shrink again. The first bubble circulation time obtained from the simulated results is 0.145 s, while the result calculated from (8) is 0.151 s. The simulated motion of bubble is in excellent agreement with that in the empirical formula. The relative error is 0.75%. The numerical computational maximum bubble radius is 1.067 m while the result of (7) is 1.04 m. The relative error is 2.5%. According to the Blake criteria (see (16)), the value εδ of verification case is equal to 0.297 which means the bubble will move from the free surface. And the bubble collapse jet does not occur. As the effect of free surface boundary condition, the bubble shape is not spherical in the migration process. Thus, the Blake criteria in this case will have some error, especially the value εδ which is near 0.442. In the previous analysis which was given by Zhang et al. [30], the similar conclusion is also given.Figure 8
Calculated evolution of the bubble at 0.3 kg detonation at a depthd = 1.0 m.
(a)
t = 0.01 s
(b)
t = 0.02 s
(c)
t = 0.03 s
(d)
t = 0.05 s
(e)
t = 0.07 s
(f)
t = 0.09 s
(g)
t = 0.11 s
(h)
t = 0.13 s
(i)
t = 0.15 s


## 5. UNDEX Simulation Results and Discussions

### 5.1. Bubble Shape and Fluid Dynamics

#### 5.1.1. Case  1:d
/
R
m
a
x
~
1.0

The bubble dynamics of case  1 (W = 0.05 kg andd = 0.6 m) are firstly calculated. Figure 9 shows the whole process of bubble shape and free surface characteristics at different time instants. As illustrated in Figure 9, the bubble shapes remain spherical in the firstt = 0.0 ms~0.8 ms. The local cavitation behind the LCSP occurs due to the reflection of shock wave front whent = 0.5 ms, and this cavitation region becomes larger att = 0.8 ms. At this time instant, the bubble shape also changes to nonspherical significantly because of the complex interaction of bubble-LSCP-free surface. The distorted bubble shape is similar to a peach which can be called “peach-bubble” here (t = 0.8 ms~1.0 ms). It should be noted that, as shown in Figure 9, the “peach-bubble” changes to the “hill” (here, which is called “hill-bubble”) at about 2.0 ms. As a result, the hill-bubble moves from the LCSP.Figure 9
Bubble shape state of case  1 (W = 0.05 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msTo illustrate the bubble impact pressure under case 1, the contour of its characteristics is shown in Figure10. The bubble pressure is dissipated rapidly during its propagating process (t = 0.0 ms~0.5 ms) in the water. And the bubble shock wave pressure front arrives at the LCSP att = 0.5 ms (the reflected pressure wave is clear in Figure 10). This reflected wave is tensile in nature, as opposed to other compressive wave effects. And it is produced from the rarefaction of the shock wave from the free surface. Since water cannot sustain a significant amount of tension, local cavitation occurs below LCSP att = 0.8 ms. Figures 10(d)–10(h) show the typical bulk cavitation zone. In the figure, the cavitation zone can be seen to be symmetric about the vertical axis.Figure 10
Bubble impact pressure contour (W = 0.05 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msFigure11 depicts the fluid particles velocity behavior of case  1. During the initial state, the gas bubble expands outward rapidly with the high-pressure and high-temperature gas in it. Att = 0.5 ms, as the reflected wave passes, the particles are acted on by gravity and atmospheric pressure. Surface effects also occur as a result of an underwater explosion bubble att = 0.8 ms. It is obvious that a spray dome jet is formed as a result of the bubble pressure pulse, whose velocity V
m
a
x is nearly equal to 80 m/s. Later, as presented in Figures 10(e)–10(h), the range of spray dome jet is becoming larger. Furthermore, the jet penetrates the upper surface of the bubble att = 2.5 ms as shown in Figure 10(h) because of the change of bubble shape. But this jet does not cause the significant damage to LCSP.Figure 11
Fluid particle velocity contour (W = 0.05 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 ms


#### 5.1.2. Case  2:d
/
R
m
a
x
~
0.5

Figure12 depicts the bubble shape characteristics of case  2 (W = 0.05 kg andd = 0.5 m). Like case  1, in the first time periodt = 0.0 ms~0.5 ms, the bubble generated by the explosion is almost spherical during its initial stage of expansion and contraction. And the incident shock wave, which is compressive, reflects from the free surface and results in a tensile reflected wave. Notice that the cavitation occurs (seet = 0.8 ms) when the absolute pressure in the water drops below the cavitation pressure, which is about a negative pressure.Figure 12
Bubble shape state of case  2 (W = 0.05 kg,d = 0.3 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msThe pressure contour characteristic is presented in Figure13. From the figure, it is obvious that the pressure distribution characteristic has some differences with that of case  1. The first difference is the cavitation area. The one of case  2 is larger during the pressure wave front propagation than that of case  1. And the second is the peak pressure value. The peak pressure value in case  2 is about two times the one in case  1.Figure 13
Bubble impact pressure contour (W = 0.05 kg,d = 0.3 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msIt is noted that the bulk cavitation area at timet = 0.8 ms as shown in Figure 12 is not similar to the one in case  1. As the distance d in case  2 is closer than that of case  1, the bubble shape is significantly different in the stage oft = 0.8 ms~2.5 ms. Particularly, the characteristics of bubble shape for case  2 are similar to a spindle (which is called “spindle bubble” here). The length of this spindle bubble becomes longer when the bubble becomes larger. From the figure, it is clear that a water hammer (similar phenomenon was found in experiment study [16]) is formed. The average velocity of this water hammer is about 40 m/s in the initial stage (t = 1.5 ms~2.0 ms; see Figure 14), and it increases quickly in the next stage (about 65 m/s,t = 2.0 ms~2.5 ms; see Figure 14). The shape of this water hammer obtained by numerical analysis is similar to that of previous experimental image results as shown in Figures 15 and 16 [16].Figure 14
Fluid particle velocity contour (W = 0.05 kg,d = 0.3 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msFigure 15
UNDEX bubble shape and water hammer [16].Figure 16
Geometric description of water hammer from present numerical investigation.


### 5.2. Structural Deformation Patterns

The structural deformation of cases  1 and 2 are presented in Figures17 and 18, respectively. In the initial stage (t = 0.0 ms~1.0 ms for case  1,t = 0.0 ms~0.8 ms for case  2), the characteristics of deflection of these two cases are similar. And the time interval (~0.2 ms) is caused by the bubble pressure propagation difference in the water. During the following stage (t = 1.0 ms~2.0 ms for case  1,t = 0.8 ms~1.5 ms for case  2), due to the local cavitation phenomenon behind the back face sheet of LCSP, the maximum structural deformation value is much smaller than that of initial stage. But the closing impact pressure of local cavitation phenomenon does not cause significant damage to LCSP. It is noted that the local large deformation near boundary is formed att = 1.5 ms in both cases  1 and 2. Compared with the bubble shape evolution contour (see Figures 9 and 12) and fluid velocity distribution (see Figures 11 and 14) at different times, the local high speed of water plume loaded on the wet face sheet of LCSP is mainly reasonable. At the last stage (t = 1.0 ms~2.0 ms for case  1 andt = 0.8 ms~1.5 ms for case  2), it is clearly observed that the region of this local large deflection is becoming larger as the velocity of water plume jet is becoming faster.Figure 17
Deformation of case  1 (W = 0.05 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msFigure 18
Deformation of case  2 (W = 0.05 kg,d = 0.3 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 ms


### 5.3. Failure Mode and Structural Damage

As stated in the previous section, serious damage and failure do not occur for cases  1 and 2. Thus, the mass of TNT charge adopted for case  3 is increased in the simulation as shown in Table1. As the bubble shape is not spherical during the interaction process between bubble and LCSP, (6) and (7) do not fit for this case. But the equation can give a reasonable prediction for bubble initial evolution.The bubble shape evolution process is shown in Figure19. Like the ones in cases  1 and 2, at the initial staget = 0.0 ms~0.8 ms, the bubble shape still remains spherical as the reflection pressure wave does not affect the bubble migration. At the stage oft = 0.5 ms~0.8 ms, the bubble shape becomes nonspherical due to the propagation of reflection pressure wave. Unlike the results of cases  1 and 2, the top-half part of UNDEX bubble connects the cavitation area (t = 0.8 ms~1.5 ms). So the “peach-bubble” (for case  1) and “spindle bubble” (for case  2) do not occur here. It should be pointed out that the simulation of case 3 is terminated att = 1.5 ms due to too small time interval of the computational stability step.Figure 19
Bubble shape state of case 3 (W = 5.00 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.3 ms
(d)
t = 0.5 ms
(e)
t = 0.6 ms
(f)
t = 0.8 ms
(g)
t = 1.4 ms
(h)
t = 1.5 msTo illustrate the damage process of UNDEX bubble collapse during the time periodt = 1.0 ms~1.5 ms, the fluid velocity contour is presented in Figure 20. It is noted that the local high velocity region is formed att = 1.0 ms (see the shape of red dashed point) when the UNDEX bubble connected with the cavitation area. And the local high velocity region changes to two parts att = 1.1 ms~1.2 ms. In this period, the shape of local high velocity region is similar to a butterfly which can be called “butterfly bubble.” This “butterfly bubble” has been reported by few studies [31, 32]. And the shape of this bubble collapse jet changes during this process due to the interaction between LCSP and bubble (see the red dashed line of Figure 20). The fluid velocity of “butterfly bubble” top part (the detailed shape of which can be referred to in Figure 21) is about 150 m/s which is very high local impact loading. Under this loading, the significant damage (large deflection) is formulated at aboutt = 1.3 ms~1.4 ms. At aboutt = 1.4 ms, the large deflection reaches its limit where a crevasse can be found as shown in Figures 19 and 20.Figure 20
Bubble jet velocity shape at time = 1.0 ms, 1.1 ms, 1.2 ms, 1.3 ms, 1.4 ms, and 1.5 ms.Figure 21
Deformation of case 3 (W = 5.00 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.3 ms
(g)
t = 1.4 ms
(h)
t = 1.5 msIn order to compare the damage characteristics of these three cases, the curves of equivalent plastic strain of representative element (central element of water contacted surface) are presented in Figure22. As plotted in Figure 22(a), the plastic occurs at time t
=
0.8 ms simultaneously with the time when reflected shock wave passed over LCSP (see Figure 9). Furthermore, the bubble jet is formulated around time t
d
a
m = 1.0 ms; the plastic damage level increased very quickly as shown in Figure 22(a). However, the largest plastic value of case  1 for central element is equal to 6.7 × 10−3, which means that the bubble induced structural damage is minor. The equivalent plastic strain results of case  2, as illustrated in Figure 22(b), give a clear description that the plastic damage occurs at time t
d
a
m
=
0.3 ms, which is much earlier than that of case  1. And this is mainly caused by the fact that reflected shock wave front arrives at the LCSP earlier for case  2 (this phenomenon also exists for case 3, in which t
d
a
m is nearly equal to 0.18 ms though the initial damage is very minor). For both cases  1 and 2, it also shows that the plastic damage increases very fast during the bubble-LCSP interaction process during which the bubble jet is formulated and can be proved in Figures 11 and 14. Extremely, the high speed local bubble jet loading on the LCSP water contacted surface makes the breach phenomenon of LCSP as depicted in Figure 21 (and this also can be validated in Figure 22(c) as fast increase of plastic strain).Figure 22
Plastic damage of three cases.
(a)
Case  1 (W = 0.05 kg,d = 0.6 m)
(b)
Case  2 (W = 0.05 kg,d = 0.3 m)
(c)
Case  3 (W = 5.00 kg,d = 0.6 m)


## 5.1. Bubble Shape and Fluid Dynamics

### 5.1.1. Case  1:d
/
R
m
a
x
~
1.0

The bubble dynamics of case  1 (W = 0.05 kg andd = 0.6 m) are firstly calculated. Figure 9 shows the whole process of bubble shape and free surface characteristics at different time instants. As illustrated in Figure 9, the bubble shapes remain spherical in the firstt = 0.0 ms~0.8 ms. The local cavitation behind the LCSP occurs due to the reflection of shock wave front whent = 0.5 ms, and this cavitation region becomes larger att = 0.8 ms. At this time instant, the bubble shape also changes to nonspherical significantly because of the complex interaction of bubble-LSCP-free surface. The distorted bubble shape is similar to a peach which can be called “peach-bubble” here (t = 0.8 ms~1.0 ms). It should be noted that, as shown in Figure 9, the “peach-bubble” changes to the “hill” (here, which is called “hill-bubble”) at about 2.0 ms. As a result, the hill-bubble moves from the LCSP.Figure 9
Bubble shape state of case  1 (W = 0.05 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msTo illustrate the bubble impact pressure under case 1, the contour of its characteristics is shown in Figure10. The bubble pressure is dissipated rapidly during its propagating process (t = 0.0 ms~0.5 ms) in the water. And the bubble shock wave pressure front arrives at the LCSP att = 0.5 ms (the reflected pressure wave is clear in Figure 10). This reflected wave is tensile in nature, as opposed to other compressive wave effects. And it is produced from the rarefaction of the shock wave from the free surface. Since water cannot sustain a significant amount of tension, local cavitation occurs below LCSP att = 0.8 ms. Figures 10(d)–10(h) show the typical bulk cavitation zone. In the figure, the cavitation zone can be seen to be symmetric about the vertical axis.Figure 10
Bubble impact pressure contour (W = 0.05 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msFigure11 depicts the fluid particles velocity behavior of case  1. During the initial state, the gas bubble expands outward rapidly with the high-pressure and high-temperature gas in it. Att = 0.5 ms, as the reflected wave passes, the particles are acted on by gravity and atmospheric pressure. Surface effects also occur as a result of an underwater explosion bubble att = 0.8 ms. It is obvious that a spray dome jet is formed as a result of the bubble pressure pulse, whose velocity V
m
a
x is nearly equal to 80 m/s. Later, as presented in Figures 10(e)–10(h), the range of spray dome jet is becoming larger. Furthermore, the jet penetrates the upper surface of the bubble att = 2.5 ms as shown in Figure 10(h) because of the change of bubble shape. But this jet does not cause the significant damage to LCSP.Figure 11
Fluid particle velocity contour (W = 0.05 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 ms


### 5.1.2. Case  2:d
/
R
m
a
x
~
0.5

Figure12 depicts the bubble shape characteristics of case  2 (W = 0.05 kg andd = 0.5 m). Like case  1, in the first time periodt = 0.0 ms~0.5 ms, the bubble generated by the explosion is almost spherical during its initial stage of expansion and contraction. And the incident shock wave, which is compressive, reflects from the free surface and results in a tensile reflected wave. Notice that the cavitation occurs (seet = 0.8 ms) when the absolute pressure in the water drops below the cavitation pressure, which is about a negative pressure.Figure 12
Bubble shape state of case  2 (W = 0.05 kg,d = 0.3 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msThe pressure contour characteristic is presented in Figure13. From the figure, it is obvious that the pressure distribution characteristic has some differences with that of case  1. The first difference is the cavitation area. The one of case  2 is larger during the pressure wave front propagation than that of case  1. And the second is the peak pressure value. The peak pressure value in case  2 is about two times the one in case  1.Figure 13
Bubble impact pressure contour (W = 0.05 kg,d = 0.3 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msIt is noted that the bulk cavitation area at timet = 0.8 ms as shown in Figure 12 is not similar to the one in case  1. As the distance d in case  2 is closer than that of case  1, the bubble shape is significantly different in the stage oft = 0.8 ms~2.5 ms. Particularly, the characteristics of bubble shape for case  2 are similar to a spindle (which is called “spindle bubble” here). The length of this spindle bubble becomes longer when the bubble becomes larger. From the figure, it is clear that a water hammer (similar phenomenon was found in experiment study [16]) is formed. The average velocity of this water hammer is about 40 m/s in the initial stage (t = 1.5 ms~2.0 ms; see Figure 14), and it increases quickly in the next stage (about 65 m/s,t = 2.0 ms~2.5 ms; see Figure 14). The shape of this water hammer obtained by numerical analysis is similar to that of previous experimental image results as shown in Figures 15 and 16 [16].Figure 14
Fluid particle velocity contour (W = 0.05 kg,d = 0.3 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msFigure 15
UNDEX bubble shape and water hammer [16].Figure 16
Geometric description of water hammer from present numerical investigation.


## 5.1.1. Case  1:d
/
R
m
a
x
~
1.0

The bubble dynamics of case  1 (W = 0.05 kg andd = 0.6 m) are firstly calculated. Figure 9 shows the whole process of bubble shape and free surface characteristics at different time instants. As illustrated in Figure 9, the bubble shapes remain spherical in the firstt = 0.0 ms~0.8 ms. The local cavitation behind the LCSP occurs due to the reflection of shock wave front whent = 0.5 ms, and this cavitation region becomes larger att = 0.8 ms. At this time instant, the bubble shape also changes to nonspherical significantly because of the complex interaction of bubble-LSCP-free surface. The distorted bubble shape is similar to a peach which can be called “peach-bubble” here (t = 0.8 ms~1.0 ms). It should be noted that, as shown in Figure 9, the “peach-bubble” changes to the “hill” (here, which is called “hill-bubble”) at about 2.0 ms. As a result, the hill-bubble moves from the LCSP.Figure 9
Bubble shape state of case  1 (W = 0.05 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msTo illustrate the bubble impact pressure under case 1, the contour of its characteristics is shown in Figure10. The bubble pressure is dissipated rapidly during its propagating process (t = 0.0 ms~0.5 ms) in the water. And the bubble shock wave pressure front arrives at the LCSP att = 0.5 ms (the reflected pressure wave is clear in Figure 10). This reflected wave is tensile in nature, as opposed to other compressive wave effects. And it is produced from the rarefaction of the shock wave from the free surface. Since water cannot sustain a significant amount of tension, local cavitation occurs below LCSP att = 0.8 ms. Figures 10(d)–10(h) show the typical bulk cavitation zone. In the figure, the cavitation zone can be seen to be symmetric about the vertical axis.Figure 10
Bubble impact pressure contour (W = 0.05 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msFigure11 depicts the fluid particles velocity behavior of case  1. During the initial state, the gas bubble expands outward rapidly with the high-pressure and high-temperature gas in it. Att = 0.5 ms, as the reflected wave passes, the particles are acted on by gravity and atmospheric pressure. Surface effects also occur as a result of an underwater explosion bubble att = 0.8 ms. It is obvious that a spray dome jet is formed as a result of the bubble pressure pulse, whose velocity V
m
a
x is nearly equal to 80 m/s. Later, as presented in Figures 10(e)–10(h), the range of spray dome jet is becoming larger. Furthermore, the jet penetrates the upper surface of the bubble att = 2.5 ms as shown in Figure 10(h) because of the change of bubble shape. But this jet does not cause the significant damage to LCSP.Figure 11
Fluid particle velocity contour (W = 0.05 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 ms


## 5.1.2. Case  2:d
/
R
m
a
x
~
0.5

Figure12 depicts the bubble shape characteristics of case  2 (W = 0.05 kg andd = 0.5 m). Like case  1, in the first time periodt = 0.0 ms~0.5 ms, the bubble generated by the explosion is almost spherical during its initial stage of expansion and contraction. And the incident shock wave, which is compressive, reflects from the free surface and results in a tensile reflected wave. Notice that the cavitation occurs (seet = 0.8 ms) when the absolute pressure in the water drops below the cavitation pressure, which is about a negative pressure.Figure 12
Bubble shape state of case  2 (W = 0.05 kg,d = 0.3 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msThe pressure contour characteristic is presented in Figure13. From the figure, it is obvious that the pressure distribution characteristic has some differences with that of case  1. The first difference is the cavitation area. The one of case  2 is larger during the pressure wave front propagation than that of case  1. And the second is the peak pressure value. The peak pressure value in case  2 is about two times the one in case  1.Figure 13
Bubble impact pressure contour (W = 0.05 kg,d = 0.3 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msIt is noted that the bulk cavitation area at timet = 0.8 ms as shown in Figure 12 is not similar to the one in case  1. As the distance d in case  2 is closer than that of case  1, the bubble shape is significantly different in the stage oft = 0.8 ms~2.5 ms. Particularly, the characteristics of bubble shape for case  2 are similar to a spindle (which is called “spindle bubble” here). The length of this spindle bubble becomes longer when the bubble becomes larger. From the figure, it is clear that a water hammer (similar phenomenon was found in experiment study [16]) is formed. The average velocity of this water hammer is about 40 m/s in the initial stage (t = 1.5 ms~2.0 ms; see Figure 14), and it increases quickly in the next stage (about 65 m/s,t = 2.0 ms~2.5 ms; see Figure 14). The shape of this water hammer obtained by numerical analysis is similar to that of previous experimental image results as shown in Figures 15 and 16 [16].Figure 14
Fluid particle velocity contour (W = 0.05 kg,d = 0.3 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msFigure 15
UNDEX bubble shape and water hammer [16].Figure 16
Geometric description of water hammer from present numerical investigation.


## 5.2. Structural Deformation Patterns

The structural deformation of cases  1 and 2 are presented in Figures17 and 18, respectively. In the initial stage (t = 0.0 ms~1.0 ms for case  1,t = 0.0 ms~0.8 ms for case  2), the characteristics of deflection of these two cases are similar. And the time interval (~0.2 ms) is caused by the bubble pressure propagation difference in the water. During the following stage (t = 1.0 ms~2.0 ms for case  1,t = 0.8 ms~1.5 ms for case  2), due to the local cavitation phenomenon behind the back face sheet of LCSP, the maximum structural deformation value is much smaller than that of initial stage. But the closing impact pressure of local cavitation phenomenon does not cause significant damage to LCSP. It is noted that the local large deformation near boundary is formed att = 1.5 ms in both cases  1 and 2. Compared with the bubble shape evolution contour (see Figures 9 and 12) and fluid velocity distribution (see Figures 11 and 14) at different times, the local high speed of water plume loaded on the wet face sheet of LCSP is mainly reasonable. At the last stage (t = 1.0 ms~2.0 ms for case  1 andt = 0.8 ms~1.5 ms for case  2), it is clearly observed that the region of this local large deflection is becoming larger as the velocity of water plume jet is becoming faster.Figure 17
Deformation of case  1 (W = 0.05 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 msFigure 18
Deformation of case  2 (W = 0.05 kg,d = 0.3 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.5 ms
(g)
t = 2.0 ms
(h)
t = 2.5 ms


## 5.3. Failure Mode and Structural Damage

As stated in the previous section, serious damage and failure do not occur for cases  1 and 2. Thus, the mass of TNT charge adopted for case  3 is increased in the simulation as shown in Table1. As the bubble shape is not spherical during the interaction process between bubble and LCSP, (6) and (7) do not fit for this case. But the equation can give a reasonable prediction for bubble initial evolution.The bubble shape evolution process is shown in Figure19. Like the ones in cases  1 and 2, at the initial staget = 0.0 ms~0.8 ms, the bubble shape still remains spherical as the reflection pressure wave does not affect the bubble migration. At the stage oft = 0.5 ms~0.8 ms, the bubble shape becomes nonspherical due to the propagation of reflection pressure wave. Unlike the results of cases  1 and 2, the top-half part of UNDEX bubble connects the cavitation area (t = 0.8 ms~1.5 ms). So the “peach-bubble” (for case  1) and “spindle bubble” (for case  2) do not occur here. It should be pointed out that the simulation of case 3 is terminated att = 1.5 ms due to too small time interval of the computational stability step.Figure 19
Bubble shape state of case 3 (W = 5.00 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.3 ms
(d)
t = 0.5 ms
(e)
t = 0.6 ms
(f)
t = 0.8 ms
(g)
t = 1.4 ms
(h)
t = 1.5 msTo illustrate the damage process of UNDEX bubble collapse during the time periodt = 1.0 ms~1.5 ms, the fluid velocity contour is presented in Figure 20. It is noted that the local high velocity region is formed att = 1.0 ms (see the shape of red dashed point) when the UNDEX bubble connected with the cavitation area. And the local high velocity region changes to two parts att = 1.1 ms~1.2 ms. In this period, the shape of local high velocity region is similar to a butterfly which can be called “butterfly bubble.” This “butterfly bubble” has been reported by few studies [31, 32]. And the shape of this bubble collapse jet changes during this process due to the interaction between LCSP and bubble (see the red dashed line of Figure 20). The fluid velocity of “butterfly bubble” top part (the detailed shape of which can be referred to in Figure 21) is about 150 m/s which is very high local impact loading. Under this loading, the significant damage (large deflection) is formulated at aboutt = 1.3 ms~1.4 ms. At aboutt = 1.4 ms, the large deflection reaches its limit where a crevasse can be found as shown in Figures 19 and 20.Figure 20
Bubble jet velocity shape at time = 1.0 ms, 1.1 ms, 1.2 ms, 1.3 ms, 1.4 ms, and 1.5 ms.Figure 21
Deformation of case 3 (W = 5.00 kg,d = 0.6 m).
(a)
t = 0.0 ms
(b)
t = 0.2 ms
(c)
t = 0.5 ms
(d)
t = 0.8 ms
(e)
t = 1.0 ms
(f)
t = 1.3 ms
(g)
t = 1.4 ms
(h)
t = 1.5 msIn order to compare the damage characteristics of these three cases, the curves of equivalent plastic strain of representative element (central element of water contacted surface) are presented in Figure22. As plotted in Figure 22(a), the plastic occurs at time t
=
0.8 ms simultaneously with the time when reflected shock wave passed over LCSP (see Figure 9). Furthermore, the bubble jet is formulated around time t
d
a
m = 1.0 ms; the plastic damage level increased very quickly as shown in Figure 22(a). However, the largest plastic value of case  1 for central element is equal to 6.7 × 10−3, which means that the bubble induced structural damage is minor. The equivalent plastic strain results of case  2, as illustrated in Figure 22(b), give a clear description that the plastic damage occurs at time t
d
a
m
=
0.3 ms, which is much earlier than that of case  1. And this is mainly caused by the fact that reflected shock wave front arrives at the LCSP earlier for case  2 (this phenomenon also exists for case 3, in which t
d
a
m is nearly equal to 0.18 ms though the initial damage is very minor). For both cases  1 and 2, it also shows that the plastic damage increases very fast during the bubble-LCSP interaction process during which the bubble jet is formulated and can be proved in Figures 11 and 14. Extremely, the high speed local bubble jet loading on the LCSP water contacted surface makes the breach phenomenon of LCSP as depicted in Figure 21 (and this also can be validated in Figure 22(c) as fast increase of plastic strain).Figure 22
Plastic damage of three cases.
(a)
Case  1 (W = 0.05 kg,d = 0.6 m)
(b)
Case  2 (W = 0.05 kg,d = 0.3 m)
(c)
Case  3 (W = 5.00 kg,d = 0.6 m)


## 6. Conclusions

In the present investigation, the high-resolution numerical analysis model of underwater explosion bubble and corrugated sandwich plate are built and developed using the fluid-solid interaction nonlinear software MSC.Dytran. The multi-Euler algorithm is adapted to describe the detailed dynamic behavior between underwater explosion bubble and corrugated sandwich plate. It is found that the response and damage of corrugated sandwich plate in the whole interaction process are significantly affected by the bubble response and such effect should not be neglected. In summary, the following conclusions from the viewpoint of gas bubble pulse dominating failure can be drawn:(i)
The underwater explosion bubble shape changes to nonspherical significantly caused by the interaction of bubble-LSCP-free surface. This distorted bubble shape is different when the ratio of charge depth to maximum bubble radius (d
/
R
m
a
x) decreases. Here, whend/R
m
a
x~1.0, the nonspherical bubble is similar to a peach. But ifd/R
m
a
x~0.5, the nonspherical bubble is similar to a spindle. Furthermore, the nonspherical bubble is similar to a butterfly when d
/
R
m
a
x
<
0.2.
(ii)
It is found that the bulk cavitation area connected with the underwater explosion bubble whend
/
R
m
a
x
≲
0.5. So, the estimation formula of maximum radius of underwater explosion bubble from previous study (e.g., (7)) cannot be used in these cases.
(iii)
Based on the nonlinear finite element analysis results, the failure modes of corrugated sandwich plate subjected to underwater explosion bubble are also clarified. The major damage of corrugated sandwich plate is large plastic deformation caused by UNDEX shock wave whend/R
m
a
x
≳
0.5. And the crevasse damage caused by butterfly bubble local high velocity jet is major failure type when d
/
R
m
a
x
≲
0.2. It is demonstrated that the UNDEX bubble collapse jet load plays a more significant role than the UNDEX shock wave load in near-field underwater explosion.



---

*Source: 6057437-2016-09-20.xml*