# Damage Characteristics of Coated Cylindrical Shells Subjected to Underwater Contact Explosion

Authors: Zhi-fan Zhang; Fu-ren Ming; A-man Zhang
Journal: Shock and Vibration
 (2014)

---

## Abstract

It is of great significance for the protective design of submarine to study the influences of coverings on the damage characteristics of single and double cylindrical shells subjected to underwater contact explosions. The SPH models of single and double cylindrical shells coated with foam silicone rubber are established to analyze shockwave propagation, damage characteristics, and elastoplastic responses, which provides reasonable parameters of covering position and thickness. The results can be concluded as follows: the superposition of multiple waves may cause the inhomogeneity and discontinuity; for the single cylindrical shell with inner or outer coverings, the damage mode is mainly tensile and shear failure is caused by detonation waves and detonation products; compared with out-covering approach, the in-covering approach has better antishock performance; the best protective effect comes out when the thickness of covering is close to that of the shell; as for the double cylindrical shell without interlayer water, the destruction of inner shell mainly results from the puncture of high-speed fragments from the outer shell, so for the outer shell, out-covering is a better choice; however, since the interlayer water is very effective in protecting the inner shell, in-covering will be better for the inner shell.


---

## Body

## 1. Introduction

With the development of precision-guided weapons, the probability of naval structures being attacked/damaged by underwater contact explosion increases gradually. Therefore, special materials or structures are always adopted to resist the severe shock. As a submerged vehicle, submarine is always coated with multiple acoustic coverings [1] in order to improve the concealment, and thus it is of great significance to study the influences of the coverings and their antishock performance. Numerous researches on the covering performance were carried out and can be classified into mainly two categories: those on the acoustic performance such as the vibration and acoustic radiation of finite submerged cylindrical shell [2] and the acoustic radiation of laminated composite shells [3], and those on the antishock performance such as the dynamic response of coated double-hull subjected to underwater explosion [4], antishock characteristics of the sandwich structure subjected to underwater explosion [5, 6], the theoretical and experimental researches on antishock characteristics of a new cavity structure [7–9], the analyses of structural dynamic responses and failure modes based on finite element program [10], and the antishock performance of plate-frame covering [11, 12]. Owing to the high-cost and poor-portability of experimental research and the limitation of theoretical research, the numerical research is more popular. Due to the presence of large-deformation, moving interface, multiphase mixing, and other challenges, traditional methods based on mesh come across many difficulties in solving such problems [13–15]. However, as the oldest mesh-free method, characterized by Lagrangian particle property [16], SPH method has great advantages in dealing with these problems [17, 18].Foam silicone rubber is chosen as the covering material since it has no significant influence on acoustic performance. SPH method is applied to study the influences of the covering on the damage characteristics of single and double cylindrical shells. Firstly, the feasibility of SPH method to solve strong discontinuity problems will be discussed in terms of dynamic continuity and motion continuity; and then the coated single and double cylindrical shells will be modeled with SPH method, the results of which will be compared with those of AUTODYN to verify the effectiveness of the present SPH model; on the basis of these, the parameters of covering position and covering thickness will be considered to discuss the influences on the structural damage of single and double cylindrical shells; more reasonable parameters of the coverings will be presented in terms of antishock performances to give a reference for the structural design of submarines.


## 2. Theoretical Background

### 2.1. SPH Approximation Equations

In SPH method, the approximation of functionf
(
x
) and its derivative ∇
·
f
(
x
) can be discretized as [19]
(1)
f
(
x
i
)
≈
∑
j
m
j
ρ
j
f
(
x
j
)
W
i
j
,
∇
·
f
(
x
i
)
≈
∑
j
=
1
N
m
j
ρ
j
f
(
x
j
)
·
∇
i
W
i
j
.Thus, the conservation of mass, momentum and energy in SPH can be expressed as (without regard to body force) [19]
(2)
[
d
ρ
i
d
t
d
v
i
a
d
t
d
e
i
d
t
]
=
[
m
i
·
∑
j
=
1
N
v
i
j
b
-
∑
j
=
1
N
m
j
ρ
j
(
σ
i
a
b
+
σ
j
a
b
ρ
i
+
Π
i
j
)
1
2
∑
j
=
1
N
m
j
ρ
j
(
σ
i
a
b
+
σ
j
a
b
ρ
i
+
Π
i
j
)
v
i
j
b
]
·
∂
W
i
j
∂
x
i
b
,
where ρ, m, v, e, t, x, and σ denote density, mass, velocity, energy, time, coordinates, and stress, respectively; a and b indicate the direction along axis; W
i
j is the smoothed function of a pair of particles i and j, the cubic spline function is applied in present paper; Π
i
j represents artificial viscosity [19].


### 2.2. Constitutive Equation

As for fluid, the stressσ is composed of isotropic pressure P and viscous shear stress τ, namely, [19]
(3)
σ
a
b
=
-
P
δ
a
b
+
τ
a
b
=
-
P
δ
a
b
+
μ
ε
˙
a
b
,
where P can be obtained from equation of state (EOS), Mie-Gruneisen EOS is used for water [20]; τ is the viscous shear stress; ε
˙
a
b is the strain rate. As for Mie-Gruneisen EOS of water [20], ρ
0 is initial density of water; C is sound velocity; η is density ratio before and after the perturbation; a is volume correction coefficient; S
1, S
2, and S
3 are fitting coefficients; E
0 is initial energy. The parameters for EOS of water are given in Table 1.Table 1
Parameters in Mie-Gruneisen EOS for water [20].
ρ
0 (kg/m3)
C (m/s)
γ
0
a
S
1
S
2
S
3
E
0 (J/kg)
1000
1480
0.5
0
2.56
1.986
1.2268
357.1As for solid, the stressσ is composed of isotropic pressure P and deviatoric stress tensor S, namely, [19]
(4)
σ
a
b
=
-
P
δ
a
b
+
S
a
b
.The solutions ofS and P are explained as follows.(1) The Solution of
S. The deviatoric stress S can be drawn from the stress rate S
˙ by integration, and yet S
˙ can be gained from Jaumann stress rate [19], which can be noted as
(5)
S
˙
a
b
=
2
G
(
ε
˙
a
b
-
δ
a
b
·
ε
˙
c
c
3
)
·
R
˙
b
c
+
S
a
c
·
R
˙
b
c
+
S
c
b
·
R
˙
c
c
,
where G is shear modulus, ε
˙ is the strain rate, and R
˙ is the rotation tensor rate.Whether the stress should be updated or not is determined by Mises yield criterion [21], in which the Mises stress σ is compared with yield strength Y to judge the stress state of particle i. If σ
≤
Y, particle i is in elastic state, and there is no need to modify the stress component. On the contrary, if σ
>
Y, particle i is in plastic state, it is necessary to modify the stress component as [21]
(6)
S
a
b
=
S
a
b
Y
σ
Mises
,
where Y means the yield strength, Johnson-cook model is applied to calculate the yield strength of the steel [22]. Meanwhile, the yield strength of the foam silicone rubber can be drawn from the dynamic constitutive equation [23].In Johnson-Cook constitutive model of steel [22], σ
0 is the static yield strength; B and n are strain hardening coefficients; C is strain rate strengthening coefficient; M is thermal softening coefficient; ε
e is equivalent plastic strain. ε
e
=
(
2
/
3
)
3
I
2, where I
2 is the second invariant of deviator stress tensor; ε
˙
e is equivalent plastic strain rate; ε
˙
0 is reference strain rate; T
ϕ is dimensionless temperature which corresponds to T; and T
ϕ
=
(
T
-
T
r
)
/
(
T
m
-
T
r
), where T
r is room temperature and T
r is melting temperature. T
=
T
r
+
(
e
-
e
0
)
/
M
C
v, where e is specific internal energy, e
0 is initial specific internal energy, and C
v is specific heat; G is shear modulus of steel. Some of the parameters are listed in Table 2. In the dynamic constitutive equation of the foam silicone rubber [23], ε
˙
e is equivalent plastic strain rate; ε
˙
0 is reference strain rate; ε
e is equivalent plastic strain; A
1, B
0, B
1, B
2, B
3, B
4, and B
5 are fitting coefficients. The detailed parameters are shown in Table 3.Table 2
Parameters in Johnson-Cook constitutive model of steel [22].
σ (MPa)
B (MPa)
ε
˙
0 (s−1)
n
C
m
T
m (K)
T
r (K)
C
v (J/kg·K)
G (GPa)
350
275
1
0.36
0.22
1.0
1811
288
452
77Table 3
Parameters in the dynamic constitutive equation of the foam silicone rubber [23].
A
1
ε
˙
0
B
0
B
1
B
2
B
3
B
4
B
5
0.03813
11022
11.85
561
−4528
16301
−26044
15570(2) The Solution of
P. The pressure P of detonation products, steel and foam silicone rubber, can be obtained from Jones-Wilkins-Lee (JWL) EOS, Mie-Gruneisen EOS for solid, and Murnaghan EOS, respectively [24–27]. For JWL EOS for explosive gas [24], ρ
0 is initial density; D is detonation velocity; A, B, R
1, R
2, and ω are experimental fitting coefficients; E
0 is detonation energy per unit mass. Then in Mie-Gruneisen EOS for steel [25], ρ
0 is initial density; Γ is Gruneisen parameter; C
s is linear participation coefficient of impact velocity and particles velocity; S
s is slope; E
0 is initial energy of steel; ν is Poisson’s ratio. Finally, as for Murnaghan EOS for foam silicone rubber [26, 27], ρ
0 is initial density; C
0 is linear participation coefficient of impact velocity and particles velocity; b is slope of Hugoniot D
-
u curve [26]. And some parameters are listed in Tables 4, 5, and 6.Table 4
Parameters in Jones-Wilkins-Lee (JWL) EOS for explosive gas [24].
ρ
0 (kg/m3)
D (m/s)
A  (Pa)
B   (Pa)
R
1
R
2
ω
E
0 (J/kg)
1630
6930
3.712
×
1
0
11
3.21
×
1
0
9
4.15
0.95
0.30
4.20
×
1
0
6Table 5
Parameters in Mie-Gruneisen EOS for steel [25].
ρ
0 (kg/m3)
Γ
C
s (m/s)
S
s
E
0 (J/kg)
v
7890
1.587
3075
1.294
0
0.3Table 6
Parameters in Murnaghan EOS for foam silicone rubber [26, 27].
ρ
0 (kg/m3)
C
0 (m/s)
b
1037
1.863
1.596


### 2.3. Continuity Condition of Discontinuity Surface

According to the requirement of continuum on wavefront, the particle displacement is continuous, but its derivatives, such as velocity and strain, may be discontinuous for shock problems. In this section, the essential conditions of stress and strain on discontinuity surface, namely, the conditions of dynamic continuity and motion continuity [28], will be deduced with SPH method.Given a finite regionΩ, as shown in Figure 1 [28], V and S denote the volume and the surface. The region is divided into two parts by a propagating shockwave Π
(
t
), the volume and surface of the two parts are V
i, V
j and S
i, S
j, respectively; the wavefront velocity in the direction of outward normal is assumed to be D, the corresponding particle velocity and normal velocity of the two parts are v
i, v
j and v
i
n, v
j
n, respectively.Figure 1
Kinematic state of shock wave [28].The functionf
(
x
,
t
) is assumed to be continuous and differentiable in V
i and V
j, yet it may be discontinuous on Π
(
t
). The particle position x is a function of time, and the density ρ is a scalar function of x and t, namely, ρ
=
ρ
(
x
,
t
).


#### 2.3.1. Dynamic Continuity Condition

With the propagation of shock wave, not only the functionf but also V
i including S
i and Π
(
t
) would change. Therefore, considering the derivative relations [28]
(7)
d
f
i
d
t
=
∂
f
i
∂
t
+
∂
f
i
∂
x
i
·
d
x
S
i
d
t
+
∂
f
i
∂
x
i
·
d
x
Π
d
t
and the particle approximation of function and its derivative, the above equation in V
i can be discretized as
(8)
d
d
t
∫
v
i
f
i
d
V
=
∫
v
i
∂
f
i
∂
t
d
V
+
∫
s
i
v
i
→
·
n
→
f
i
d
s
+
∫
Π
D
→
·
n
→
f
i
d
s
=
∫
v
i
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
∂
x
∂
t
d
V
+
∫
s
i
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
v
i
n
d
s
+
∫
Π
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
D
d
s
.
Similarly, the equation in V
j can be obtained as follows:
(9)
d
d
t
∫
v
j
f
j
d
V
=
∫
v
j
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
∂
x
∂
t
d
V
+
∫
s
j
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
v
j
n
d
s
+
∫
Π
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
(
-
D
)
d
s
.(1) Replace
f
with Function
ρ
in (8) and (9), namely,
f
(
x
,
t
)
=
ρ
(
x
,
t
). Assume that V
i and V
j tend to the wavefront Π
(
t
) and approach zero, and thus, both the boundaries S
i and S
j tend to Π
(
t
). Since the outward normal and Π
(
t
) are opposites when S
i
→
Π
(
t
), v
i
n should be -
v
i
n in (8). Then, adding (8) to (9),
(10)
d
d
t
(
∫
v
i
ρ
i
d
V
+
∫
v
j
ρ
j
d
V
)
=
∫
Π
[
-
∑
j
m
j
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
+
∑
i
m
i
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
]
d
s
.Due to the mass conservation in any volumeV and the arbitrariness of the wave front Π, (10) can be written as
(11)
∑
j
m
j
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
=
∑
i
m
i
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
.(2) Replace
f
in (8) and (9) as
f
i
=
ρ
i
v
i
α
,
f
j
=
ρ
j
v
j
α. Similarly, (8) is added to (9) as follows:
(12)
d
d
t
(
∫
v
i
ρ
i
v
i
α
d
V
+
∫
v
j
ρ
j
v
j
α
d
V
)
=
∫
Π
[
-
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
+
∑
i
m
i
ρ
i
ρ
i
v
i
α
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
]
d
s
.As a result of the opposite directions betweenv
i
α and Π
(
t
) when S
i
→
Π
(
t
), by applying Green formula, the momentum conservation equation regardless of the external force [19]
(13)
d
v
α
d
t
=
1
ρ
·
∂
σ
α
β
∂
x
β
can be derived as [28]
(14)
d
d
t
∫
v
i
ρ
i
v
i
α
d
V
=
-
∫
v
i
(
σ
α
β
,
β
)
i
d
V
=
-
∫
Π
σ
i
α
β
·
n
β
d
s
,
d
d
t
∫
v
i
ρ
i
v
j
α
d
V
=
∫
v
i
(
σ
α
β
,
β
)
j
d
V
=
∫
Π
σ
j
α
β
·
n
β
d
s
.Substituting the above equations into (12), we have the following:
(15)
-
∫
Π
σ
i
α
β
·
n
β
d
s
+
∫
Π
σ
j
α
β
·
n
β
d
s
=
∫
Π
[
-
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
+
∑
i
m
i
ρ
i
ρ
i
v
i
α
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
]
d
s
.Owing to the arbitrariness of wavefrontΠ, we have
(16)
(
σ
j
α
β
-
σ
i
α
β
)
·
n
β
=
-
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
+
∑
i
m
i
ρ
i
ρ
i
v
i
α
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
.Substituting (11) into the above equation yields
(17)
(
σ
j
α
β
-
σ
i
α
β
)
·
n
β
=
(
v
i
n
-
D
)
[
-
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
+
∑
j
(
m
j
/
ρ
j
)
ρ
j
(
∂
W
i
j
/
∂
x
i
)
∑
i
(
m
i
/
ρ
i
)
ρ
i
(
∂
W
i
j
/
∂
x
j
)
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
·
∑
i
m
i
ρ
i
ρ
i
v
i
α
∂
W
i
j
∂
x
j
]
.Applying the particle approximation, the following equation can be drawn:(18)
(
σ
j
α
β
-
σ
i
α
β
)
·
n
β
=
(
v
i
n
-
D
)
(
-
ρ
i
v
i
α
+
ρ
i
ρ
j
·
ρ
j
v
j
α
)
=
ρ
i
(
v
i
n
-
D
)
(
v
j
α
-
v
i
α
)
,
which is the dynamic continuity condition [28].


#### 2.3.2. Motion Continuity Condition

The discontinuity surface turns to beΣ
(
t
+
Δ
t
) after Δ
t from Σ
(
t
). Therefore, attend t
+
Δ
t, the discontinuous values on the discontinuity surface can be given as [28]
(19)
[
f
]
=
f
j
-
f
i
,
[
f
′
]
=
f
j
′
-
f
i
′
.Then duringΔ
t, the rate of discontinuous value [
f
] with respect to time is [28] as follows:
(20)
Δ
[
f
]
Δ
t
=
[
f
′
]
Δ
t
-
[
f
]
Δ
t
=
Δ
f
j
Δ
t
-
Δ
f
i
Δ
t
.IfΔ
t
→
0, the equation above turns to be [28] as follows:
(21)
d
[
f
]
d
t
=
d
f
j
d
t
-
d
f
i
d
t
=
[
d
f
d
t
]
.In general, the propagation velocityD is much larger than the particle velocity v
i
n and v
j
n, and thus the particle velocity in (7) can be ignored [28] as follows:
(22)
(
d
f
d
t
)
i
=
(
∂
f
∂
t
)
i
+
(
∂
f
∂
x
)
i
·
d
x
Π
d
t
,
(23)
(
d
f
d
t
)
j
=
(
∂
f
∂
t
)
j
+
(
∂
f
∂
x
)
j
·
d
x
Π
d
t
.Discretizing the two equations above, subtract (23) from (22) as follows:
(24)
d
[
f
]
d
t
=
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
∂
x
i
∂
t
-
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
∂
x
j
∂
t
+
[
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
-
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
]
·
D
·
n
α
.Replacef
(
x
,
t
) with the displacement u
(
x
,
t
), due to the continuity of the particle displacement, namely, d
[
f
]
/
d
t
=
0, so we get
(25)
∑
j
m
j
ρ
j
u
j
∂
W
i
j
∂
X
i
∂
X
i
∂
t
-
∑
i
m
i
ρ
i
u
i
∂
W
i
j
∂
X
j
∂
X
j
∂
t
=
-
[
∑
j
m
j
ρ
j
u
j
∂
W
i
j
∂
X
i
-
∑
i
m
i
ρ
i
u
i
∂
W
i
j
∂
X
j
]
·
D
·
n
α
.Therefore, restoring the discretized equation (25), in Cartesian coordinate, it turns to be as follows:
(26)
[
∂
u
α
∂
t
]
=
-
D
n
β
[
u
α
,
β
]
.According to the Helmholtz velocity decomposing theorem [29], in the case of small-deformation, the deformation can be divided into linear deformation and angular deformation. Therefore, the motion continuity condition can be derived as [28]
(27)
[
v
α
]
=
-
D
n
β
[
ε
α
β
+
ω
α
β
]
.In addition, we use a slight penalty force of Lennard-Jones model to solve interface problem, and the molecular force is so slight that it just prevents particles’ penetration. When particles on both sides of an interface tend to penetrate, in the case wherer
0
≥
r
i
j, there would be a molecular force F
i
j
=
f
·
[
(
r
0
/
r
i
j
)
a
-
(
r
0
/
r
i
j
)
b
]
·
x
i
j
/
r
i
j
2 acting on two approaching particles, where r
0 is cutoff radius and it is generally close to the initial spacing of particles and r
i
j is the distance between particles i and j. The direction of the F
i
j is along the center line of the particles, preventing the particle from penetrating;f, a, and b are set parameters. Consequently, we use the ratio of r
0 to r
i
j to deal with the interface; this can guarantee good numerical stability without errors caused by the smoothing length.


## 2.1. SPH Approximation Equations

In SPH method, the approximation of functionf
(
x
) and its derivative ∇
·
f
(
x
) can be discretized as [19]
(1)
f
(
x
i
)
≈
∑
j
m
j
ρ
j
f
(
x
j
)
W
i
j
,
∇
·
f
(
x
i
)
≈
∑
j
=
1
N
m
j
ρ
j
f
(
x
j
)
·
∇
i
W
i
j
.Thus, the conservation of mass, momentum and energy in SPH can be expressed as (without regard to body force) [19]
(2)
[
d
ρ
i
d
t
d
v
i
a
d
t
d
e
i
d
t
]
=
[
m
i
·
∑
j
=
1
N
v
i
j
b
-
∑
j
=
1
N
m
j
ρ
j
(
σ
i
a
b
+
σ
j
a
b
ρ
i
+
Π
i
j
)
1
2
∑
j
=
1
N
m
j
ρ
j
(
σ
i
a
b
+
σ
j
a
b
ρ
i
+
Π
i
j
)
v
i
j
b
]
·
∂
W
i
j
∂
x
i
b
,
where ρ, m, v, e, t, x, and σ denote density, mass, velocity, energy, time, coordinates, and stress, respectively; a and b indicate the direction along axis; W
i
j is the smoothed function of a pair of particles i and j, the cubic spline function is applied in present paper; Π
i
j represents artificial viscosity [19].


## 2.2. Constitutive Equation

As for fluid, the stressσ is composed of isotropic pressure P and viscous shear stress τ, namely, [19]
(3)
σ
a
b
=
-
P
δ
a
b
+
τ
a
b
=
-
P
δ
a
b
+
μ
ε
˙
a
b
,
where P can be obtained from equation of state (EOS), Mie-Gruneisen EOS is used for water [20]; τ is the viscous shear stress; ε
˙
a
b is the strain rate. As for Mie-Gruneisen EOS of water [20], ρ
0 is initial density of water; C is sound velocity; η is density ratio before and after the perturbation; a is volume correction coefficient; S
1, S
2, and S
3 are fitting coefficients; E
0 is initial energy. The parameters for EOS of water are given in Table 1.Table 1
Parameters in Mie-Gruneisen EOS for water [20].
ρ
0 (kg/m3)
C (m/s)
γ
0
a
S
1
S
2
S
3
E
0 (J/kg)
1000
1480
0.5
0
2.56
1.986
1.2268
357.1As for solid, the stressσ is composed of isotropic pressure P and deviatoric stress tensor S, namely, [19]
(4)
σ
a
b
=
-
P
δ
a
b
+
S
a
b
.The solutions ofS and P are explained as follows.(1) The Solution of
S. The deviatoric stress S can be drawn from the stress rate S
˙ by integration, and yet S
˙ can be gained from Jaumann stress rate [19], which can be noted as
(5)
S
˙
a
b
=
2
G
(
ε
˙
a
b
-
δ
a
b
·
ε
˙
c
c
3
)
·
R
˙
b
c
+
S
a
c
·
R
˙
b
c
+
S
c
b
·
R
˙
c
c
,
where G is shear modulus, ε
˙ is the strain rate, and R
˙ is the rotation tensor rate.Whether the stress should be updated or not is determined by Mises yield criterion [21], in which the Mises stress σ is compared with yield strength Y to judge the stress state of particle i. If σ
≤
Y, particle i is in elastic state, and there is no need to modify the stress component. On the contrary, if σ
>
Y, particle i is in plastic state, it is necessary to modify the stress component as [21]
(6)
S
a
b
=
S
a
b
Y
σ
Mises
,
where Y means the yield strength, Johnson-cook model is applied to calculate the yield strength of the steel [22]. Meanwhile, the yield strength of the foam silicone rubber can be drawn from the dynamic constitutive equation [23].In Johnson-Cook constitutive model of steel [22], σ
0 is the static yield strength; B and n are strain hardening coefficients; C is strain rate strengthening coefficient; M is thermal softening coefficient; ε
e is equivalent plastic strain. ε
e
=
(
2
/
3
)
3
I
2, where I
2 is the second invariant of deviator stress tensor; ε
˙
e is equivalent plastic strain rate; ε
˙
0 is reference strain rate; T
ϕ is dimensionless temperature which corresponds to T; and T
ϕ
=
(
T
-
T
r
)
/
(
T
m
-
T
r
), where T
r is room temperature and T
r is melting temperature. T
=
T
r
+
(
e
-
e
0
)
/
M
C
v, where e is specific internal energy, e
0 is initial specific internal energy, and C
v is specific heat; G is shear modulus of steel. Some of the parameters are listed in Table 2. In the dynamic constitutive equation of the foam silicone rubber [23], ε
˙
e is equivalent plastic strain rate; ε
˙
0 is reference strain rate; ε
e is equivalent plastic strain; A
1, B
0, B
1, B
2, B
3, B
4, and B
5 are fitting coefficients. The detailed parameters are shown in Table 3.Table 2
Parameters in Johnson-Cook constitutive model of steel [22].
σ (MPa)
B (MPa)
ε
˙
0 (s−1)
n
C
m
T
m (K)
T
r (K)
C
v (J/kg·K)
G (GPa)
350
275
1
0.36
0.22
1.0
1811
288
452
77Table 3
Parameters in the dynamic constitutive equation of the foam silicone rubber [23].
A
1
ε
˙
0
B
0
B
1
B
2
B
3
B
4
B
5
0.03813
11022
11.85
561
−4528
16301
−26044
15570(2) The Solution of
P. The pressure P of detonation products, steel and foam silicone rubber, can be obtained from Jones-Wilkins-Lee (JWL) EOS, Mie-Gruneisen EOS for solid, and Murnaghan EOS, respectively [24–27]. For JWL EOS for explosive gas [24], ρ
0 is initial density; D is detonation velocity; A, B, R
1, R
2, and ω are experimental fitting coefficients; E
0 is detonation energy per unit mass. Then in Mie-Gruneisen EOS for steel [25], ρ
0 is initial density; Γ is Gruneisen parameter; C
s is linear participation coefficient of impact velocity and particles velocity; S
s is slope; E
0 is initial energy of steel; ν is Poisson’s ratio. Finally, as for Murnaghan EOS for foam silicone rubber [26, 27], ρ
0 is initial density; C
0 is linear participation coefficient of impact velocity and particles velocity; b is slope of Hugoniot D
-
u curve [26]. And some parameters are listed in Tables 4, 5, and 6.Table 4
Parameters in Jones-Wilkins-Lee (JWL) EOS for explosive gas [24].
ρ
0 (kg/m3)
D (m/s)
A  (Pa)
B   (Pa)
R
1
R
2
ω
E
0 (J/kg)
1630
6930
3.712
×
1
0
11
3.21
×
1
0
9
4.15
0.95
0.30
4.20
×
1
0
6Table 5
Parameters in Mie-Gruneisen EOS for steel [25].
ρ
0 (kg/m3)
Γ
C
s (m/s)
S
s
E
0 (J/kg)
v
7890
1.587
3075
1.294
0
0.3Table 6
Parameters in Murnaghan EOS for foam silicone rubber [26, 27].
ρ
0 (kg/m3)
C
0 (m/s)
b
1037
1.863
1.596


## 2.3. Continuity Condition of Discontinuity Surface

According to the requirement of continuum on wavefront, the particle displacement is continuous, but its derivatives, such as velocity and strain, may be discontinuous for shock problems. In this section, the essential conditions of stress and strain on discontinuity surface, namely, the conditions of dynamic continuity and motion continuity [28], will be deduced with SPH method.Given a finite regionΩ, as shown in Figure 1 [28], V and S denote the volume and the surface. The region is divided into two parts by a propagating shockwave Π
(
t
), the volume and surface of the two parts are V
i, V
j and S
i, S
j, respectively; the wavefront velocity in the direction of outward normal is assumed to be D, the corresponding particle velocity and normal velocity of the two parts are v
i, v
j and v
i
n, v
j
n, respectively.Figure 1
Kinematic state of shock wave [28].The functionf
(
x
,
t
) is assumed to be continuous and differentiable in V
i and V
j, yet it may be discontinuous on Π
(
t
). The particle position x is a function of time, and the density ρ is a scalar function of x and t, namely, ρ
=
ρ
(
x
,
t
).


### 2.3.1. Dynamic Continuity Condition

With the propagation of shock wave, not only the functionf but also V
i including S
i and Π
(
t
) would change. Therefore, considering the derivative relations [28]
(7)
d
f
i
d
t
=
∂
f
i
∂
t
+
∂
f
i
∂
x
i
·
d
x
S
i
d
t
+
∂
f
i
∂
x
i
·
d
x
Π
d
t
and the particle approximation of function and its derivative, the above equation in V
i can be discretized as
(8)
d
d
t
∫
v
i
f
i
d
V
=
∫
v
i
∂
f
i
∂
t
d
V
+
∫
s
i
v
i
→
·
n
→
f
i
d
s
+
∫
Π
D
→
·
n
→
f
i
d
s
=
∫
v
i
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
∂
x
∂
t
d
V
+
∫
s
i
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
v
i
n
d
s
+
∫
Π
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
D
d
s
.
Similarly, the equation in V
j can be obtained as follows:
(9)
d
d
t
∫
v
j
f
j
d
V
=
∫
v
j
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
∂
x
∂
t
d
V
+
∫
s
j
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
v
j
n
d
s
+
∫
Π
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
(
-
D
)
d
s
.(1) Replace
f
with Function
ρ
in (8) and (9), namely,
f
(
x
,
t
)
=
ρ
(
x
,
t
). Assume that V
i and V
j tend to the wavefront Π
(
t
) and approach zero, and thus, both the boundaries S
i and S
j tend to Π
(
t
). Since the outward normal and Π
(
t
) are opposites when S
i
→
Π
(
t
), v
i
n should be -
v
i
n in (8). Then, adding (8) to (9),
(10)
d
d
t
(
∫
v
i
ρ
i
d
V
+
∫
v
j
ρ
j
d
V
)
=
∫
Π
[
-
∑
j
m
j
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
+
∑
i
m
i
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
]
d
s
.Due to the mass conservation in any volumeV and the arbitrariness of the wave front Π, (10) can be written as
(11)
∑
j
m
j
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
=
∑
i
m
i
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
.(2) Replace
f
in (8) and (9) as
f
i
=
ρ
i
v
i
α
,
f
j
=
ρ
j
v
j
α. Similarly, (8) is added to (9) as follows:
(12)
d
d
t
(
∫
v
i
ρ
i
v
i
α
d
V
+
∫
v
j
ρ
j
v
j
α
d
V
)
=
∫
Π
[
-
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
+
∑
i
m
i
ρ
i
ρ
i
v
i
α
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
]
d
s
.As a result of the opposite directions betweenv
i
α and Π
(
t
) when S
i
→
Π
(
t
), by applying Green formula, the momentum conservation equation regardless of the external force [19]
(13)
d
v
α
d
t
=
1
ρ
·
∂
σ
α
β
∂
x
β
can be derived as [28]
(14)
d
d
t
∫
v
i
ρ
i
v
i
α
d
V
=
-
∫
v
i
(
σ
α
β
,
β
)
i
d
V
=
-
∫
Π
σ
i
α
β
·
n
β
d
s
,
d
d
t
∫
v
i
ρ
i
v
j
α
d
V
=
∫
v
i
(
σ
α
β
,
β
)
j
d
V
=
∫
Π
σ
j
α
β
·
n
β
d
s
.Substituting the above equations into (12), we have the following:
(15)
-
∫
Π
σ
i
α
β
·
n
β
d
s
+
∫
Π
σ
j
α
β
·
n
β
d
s
=
∫
Π
[
-
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
+
∑
i
m
i
ρ
i
ρ
i
v
i
α
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
]
d
s
.Owing to the arbitrariness of wavefrontΠ, we have
(16)
(
σ
j
α
β
-
σ
i
α
β
)
·
n
β
=
-
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
+
∑
i
m
i
ρ
i
ρ
i
v
i
α
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
.Substituting (11) into the above equation yields
(17)
(
σ
j
α
β
-
σ
i
α
β
)
·
n
β
=
(
v
i
n
-
D
)
[
-
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
+
∑
j
(
m
j
/
ρ
j
)
ρ
j
(
∂
W
i
j
/
∂
x
i
)
∑
i
(
m
i
/
ρ
i
)
ρ
i
(
∂
W
i
j
/
∂
x
j
)
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
·
∑
i
m
i
ρ
i
ρ
i
v
i
α
∂
W
i
j
∂
x
j
]
.Applying the particle approximation, the following equation can be drawn:(18)
(
σ
j
α
β
-
σ
i
α
β
)
·
n
β
=
(
v
i
n
-
D
)
(
-
ρ
i
v
i
α
+
ρ
i
ρ
j
·
ρ
j
v
j
α
)
=
ρ
i
(
v
i
n
-
D
)
(
v
j
α
-
v
i
α
)
,
which is the dynamic continuity condition [28].


### 2.3.2. Motion Continuity Condition

The discontinuity surface turns to beΣ
(
t
+
Δ
t
) after Δ
t from Σ
(
t
). Therefore, attend t
+
Δ
t, the discontinuous values on the discontinuity surface can be given as [28]
(19)
[
f
]
=
f
j
-
f
i
,
[
f
′
]
=
f
j
′
-
f
i
′
.Then duringΔ
t, the rate of discontinuous value [
f
] with respect to time is [28] as follows:
(20)
Δ
[
f
]
Δ
t
=
[
f
′
]
Δ
t
-
[
f
]
Δ
t
=
Δ
f
j
Δ
t
-
Δ
f
i
Δ
t
.IfΔ
t
→
0, the equation above turns to be [28] as follows:
(21)
d
[
f
]
d
t
=
d
f
j
d
t
-
d
f
i
d
t
=
[
d
f
d
t
]
.In general, the propagation velocityD is much larger than the particle velocity v
i
n and v
j
n, and thus the particle velocity in (7) can be ignored [28] as follows:
(22)
(
d
f
d
t
)
i
=
(
∂
f
∂
t
)
i
+
(
∂
f
∂
x
)
i
·
d
x
Π
d
t
,
(23)
(
d
f
d
t
)
j
=
(
∂
f
∂
t
)
j
+
(
∂
f
∂
x
)
j
·
d
x
Π
d
t
.Discretizing the two equations above, subtract (23) from (22) as follows:
(24)
d
[
f
]
d
t
=
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
∂
x
i
∂
t
-
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
∂
x
j
∂
t
+
[
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
-
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
]
·
D
·
n
α
.Replacef
(
x
,
t
) with the displacement u
(
x
,
t
), due to the continuity of the particle displacement, namely, d
[
f
]
/
d
t
=
0, so we get
(25)
∑
j
m
j
ρ
j
u
j
∂
W
i
j
∂
X
i
∂
X
i
∂
t
-
∑
i
m
i
ρ
i
u
i
∂
W
i
j
∂
X
j
∂
X
j
∂
t
=
-
[
∑
j
m
j
ρ
j
u
j
∂
W
i
j
∂
X
i
-
∑
i
m
i
ρ
i
u
i
∂
W
i
j
∂
X
j
]
·
D
·
n
α
.Therefore, restoring the discretized equation (25), in Cartesian coordinate, it turns to be as follows:
(26)
[
∂
u
α
∂
t
]
=
-
D
n
β
[
u
α
,
β
]
.According to the Helmholtz velocity decomposing theorem [29], in the case of small-deformation, the deformation can be divided into linear deformation and angular deformation. Therefore, the motion continuity condition can be derived as [28]
(27)
[
v
α
]
=
-
D
n
β
[
ε
α
β
+
ω
α
β
]
.In addition, we use a slight penalty force of Lennard-Jones model to solve interface problem, and the molecular force is so slight that it just prevents particles’ penetration. When particles on both sides of an interface tend to penetrate, in the case wherer
0
≥
r
i
j, there would be a molecular force F
i
j
=
f
·
[
(
r
0
/
r
i
j
)
a
-
(
r
0
/
r
i
j
)
b
]
·
x
i
j
/
r
i
j
2 acting on two approaching particles, where r
0 is cutoff radius and it is generally close to the initial spacing of particles and r
i
j is the distance between particles i and j. The direction of the F
i
j is along the center line of the particles, preventing the particle from penetrating;f, a, and b are set parameters. Consequently, we use the ratio of r
0 to r
i
j to deal with the interface; this can guarantee good numerical stability without errors caused by the smoothing length.


## 2.3.1. Dynamic Continuity Condition

With the propagation of shock wave, not only the functionf but also V
i including S
i and Π
(
t
) would change. Therefore, considering the derivative relations [28]
(7)
d
f
i
d
t
=
∂
f
i
∂
t
+
∂
f
i
∂
x
i
·
d
x
S
i
d
t
+
∂
f
i
∂
x
i
·
d
x
Π
d
t
and the particle approximation of function and its derivative, the above equation in V
i can be discretized as
(8)
d
d
t
∫
v
i
f
i
d
V
=
∫
v
i
∂
f
i
∂
t
d
V
+
∫
s
i
v
i
→
·
n
→
f
i
d
s
+
∫
Π
D
→
·
n
→
f
i
d
s
=
∫
v
i
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
∂
x
∂
t
d
V
+
∫
s
i
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
v
i
n
d
s
+
∫
Π
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
D
d
s
.
Similarly, the equation in V
j can be obtained as follows:
(9)
d
d
t
∫
v
j
f
j
d
V
=
∫
v
j
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
∂
x
∂
t
d
V
+
∫
s
j
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
v
j
n
d
s
+
∫
Π
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
(
-
D
)
d
s
.(1) Replace
f
with Function
ρ
in (8) and (9), namely,
f
(
x
,
t
)
=
ρ
(
x
,
t
). Assume that V
i and V
j tend to the wavefront Π
(
t
) and approach zero, and thus, both the boundaries S
i and S
j tend to Π
(
t
). Since the outward normal and Π
(
t
) are opposites when S
i
→
Π
(
t
), v
i
n should be -
v
i
n in (8). Then, adding (8) to (9),
(10)
d
d
t
(
∫
v
i
ρ
i
d
V
+
∫
v
j
ρ
j
d
V
)
=
∫
Π
[
-
∑
j
m
j
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
+
∑
i
m
i
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
]
d
s
.Due to the mass conservation in any volumeV and the arbitrariness of the wave front Π, (10) can be written as
(11)
∑
j
m
j
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
=
∑
i
m
i
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
.(2) Replace
f
in (8) and (9) as
f
i
=
ρ
i
v
i
α
,
f
j
=
ρ
j
v
j
α. Similarly, (8) is added to (9) as follows:
(12)
d
d
t
(
∫
v
i
ρ
i
v
i
α
d
V
+
∫
v
j
ρ
j
v
j
α
d
V
)
=
∫
Π
[
-
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
+
∑
i
m
i
ρ
i
ρ
i
v
i
α
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
]
d
s
.As a result of the opposite directions betweenv
i
α and Π
(
t
) when S
i
→
Π
(
t
), by applying Green formula, the momentum conservation equation regardless of the external force [19]
(13)
d
v
α
d
t
=
1
ρ
·
∂
σ
α
β
∂
x
β
can be derived as [28]
(14)
d
d
t
∫
v
i
ρ
i
v
i
α
d
V
=
-
∫
v
i
(
σ
α
β
,
β
)
i
d
V
=
-
∫
Π
σ
i
α
β
·
n
β
d
s
,
d
d
t
∫
v
i
ρ
i
v
j
α
d
V
=
∫
v
i
(
σ
α
β
,
β
)
j
d
V
=
∫
Π
σ
j
α
β
·
n
β
d
s
.Substituting the above equations into (12), we have the following:
(15)
-
∫
Π
σ
i
α
β
·
n
β
d
s
+
∫
Π
σ
j
α
β
·
n
β
d
s
=
∫
Π
[
-
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
+
∑
i
m
i
ρ
i
ρ
i
v
i
α
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
]
d
s
.Owing to the arbitrariness of wavefrontΠ, we have
(16)
(
σ
j
α
β
-
σ
i
α
β
)
·
n
β
=
-
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
(
v
i
n
-
D
)
+
∑
i
m
i
ρ
i
ρ
i
v
i
α
∂
W
i
j
∂
x
j
(
v
j
n
-
D
)
.Substituting (11) into the above equation yields
(17)
(
σ
j
α
β
-
σ
i
α
β
)
·
n
β
=
(
v
i
n
-
D
)
[
-
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
+
∑
j
(
m
j
/
ρ
j
)
ρ
j
(
∂
W
i
j
/
∂
x
i
)
∑
i
(
m
i
/
ρ
i
)
ρ
i
(
∂
W
i
j
/
∂
x
j
)
∑
j
m
j
ρ
j
ρ
j
v
j
α
∂
W
i
j
∂
x
i
·
∑
i
m
i
ρ
i
ρ
i
v
i
α
∂
W
i
j
∂
x
j
]
.Applying the particle approximation, the following equation can be drawn:(18)
(
σ
j
α
β
-
σ
i
α
β
)
·
n
β
=
(
v
i
n
-
D
)
(
-
ρ
i
v
i
α
+
ρ
i
ρ
j
·
ρ
j
v
j
α
)
=
ρ
i
(
v
i
n
-
D
)
(
v
j
α
-
v
i
α
)
,
which is the dynamic continuity condition [28].


## 2.3.2. Motion Continuity Condition

The discontinuity surface turns to beΣ
(
t
+
Δ
t
) after Δ
t from Σ
(
t
). Therefore, attend t
+
Δ
t, the discontinuous values on the discontinuity surface can be given as [28]
(19)
[
f
]
=
f
j
-
f
i
,
[
f
′
]
=
f
j
′
-
f
i
′
.Then duringΔ
t, the rate of discontinuous value [
f
] with respect to time is [28] as follows:
(20)
Δ
[
f
]
Δ
t
=
[
f
′
]
Δ
t
-
[
f
]
Δ
t
=
Δ
f
j
Δ
t
-
Δ
f
i
Δ
t
.IfΔ
t
→
0, the equation above turns to be [28] as follows:
(21)
d
[
f
]
d
t
=
d
f
j
d
t
-
d
f
i
d
t
=
[
d
f
d
t
]
.In general, the propagation velocityD is much larger than the particle velocity v
i
n and v
j
n, and thus the particle velocity in (7) can be ignored [28] as follows:
(22)
(
d
f
d
t
)
i
=
(
∂
f
∂
t
)
i
+
(
∂
f
∂
x
)
i
·
d
x
Π
d
t
,
(23)
(
d
f
d
t
)
j
=
(
∂
f
∂
t
)
j
+
(
∂
f
∂
x
)
j
·
d
x
Π
d
t
.Discretizing the two equations above, subtract (23) from (22) as follows:
(24)
d
[
f
]
d
t
=
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
∂
x
i
∂
t
-
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
∂
x
j
∂
t
+
[
∑
j
m
j
ρ
j
f
j
∂
W
i
j
∂
x
i
-
∑
i
m
i
ρ
i
f
i
∂
W
i
j
∂
x
j
]
·
D
·
n
α
.Replacef
(
x
,
t
) with the displacement u
(
x
,
t
), due to the continuity of the particle displacement, namely, d
[
f
]
/
d
t
=
0, so we get
(25)
∑
j
m
j
ρ
j
u
j
∂
W
i
j
∂
X
i
∂
X
i
∂
t
-
∑
i
m
i
ρ
i
u
i
∂
W
i
j
∂
X
j
∂
X
j
∂
t
=
-
[
∑
j
m
j
ρ
j
u
j
∂
W
i
j
∂
X
i
-
∑
i
m
i
ρ
i
u
i
∂
W
i
j
∂
X
j
]
·
D
·
n
α
.Therefore, restoring the discretized equation (25), in Cartesian coordinate, it turns to be as follows:
(26)
[
∂
u
α
∂
t
]
=
-
D
n
β
[
u
α
,
β
]
.According to the Helmholtz velocity decomposing theorem [29], in the case of small-deformation, the deformation can be divided into linear deformation and angular deformation. Therefore, the motion continuity condition can be derived as [28]
(27)
[
v
α
]
=
-
D
n
β
[
ε
α
β
+
ω
α
β
]
.In addition, we use a slight penalty force of Lennard-Jones model to solve interface problem, and the molecular force is so slight that it just prevents particles’ penetration. When particles on both sides of an interface tend to penetrate, in the case wherer
0
≥
r
i
j, there would be a molecular force F
i
j
=
f
·
[
(
r
0
/
r
i
j
)
a
-
(
r
0
/
r
i
j
)
b
]
·
x
i
j
/
r
i
j
2 acting on two approaching particles, where r
0 is cutoff radius and it is generally close to the initial spacing of particles and r
i
j is the distance between particles i and j. The direction of the F
i
j is along the center line of the particles, preventing the particle from penetrating;f, a, and b are set parameters. Consequently, we use the ratio of r
0 to r
i
j to deal with the interface; this can guarantee good numerical stability without errors caused by the smoothing length.


## 3. Contact Explosion of Single Cylindrical Shell

### 3.1. Numerical Model and Parameters

The model of a coated single-hull submarine subjected to the attack of torpedo can be simplified as Figure2, and (a) indicates no-covering, (b) denotes in-covering, and (c) represents out-covering. The relevant parameters and cases are listed in Tables 7 and 8. The origin of coordinates is positioned at the center of ring-shaped water. The column TNT is detonated from the center. The TNT dosage is about 12.8 kg, a popular amount in most model tests. The entire model is discretized with nonuniform particle spacing, 0.004 m for the cylindrical shell and 0.008 m for TNT and water, and some verification for the initial particle spacing is explained in Section 3.2. The numbers of different kinds of particles are listed as Table 8. The test point A is placed at (0, 2.1, 0).Table 7
Parameters of single cylindrical shell.
Outside radius of water (m)
Inside radius of water (m)
Inside radius of shell (m)
Thickness of shell (m)
Radius of TNT (m)
3.0
1.22
1.2
0.02
0.05Table 8
Different cases of single cylindrical shell and numbers of particles (different positions are listed as Cases  1–3, in the case where the rubber thickness is 0.020 m; different thicknesses are listed as Cases  4–7, in the case where the inner surface is coated with rubber).
Cases
Name
Numbers of particles
Water
Steel
TNT
Rubber
Total
Case  1
No-covering
364507
9300
121
0
373928
Case  2
In-covering
364507
9300
121
9153
383081
Case  3
Out-covering
355250
9300
121
9257
373928
Case  4
0.010 m
364507
9300
121
13680
387608
Case  5
0.015 m
364507
9300
121
11406
385334
Case  6
0.025 m
364507
9300
121
6886
380814
Case  7
0.030 m
364507
9300
121
4575
378503Models of single-hull.
(a)
(b)
(c)


### 3.2. Numerical Verification

In order to verify the effectiveness of the present SPH method, the no-covering model of underwater contact explosion is selected, the SPH result will be compared with that of AUTODYN.The propagations of shockwave obtained from different methods are shown in Figure3 and the pressure of test point A is shown in Figure 4. It is obvious that the propagation laws of shockwave agree well with each other; moreover, the curves obtained from different methods also show good agreements, which all prove the effectiveness and accuracy of the present SPH method to solve the problems of underwater contact explosion.Pressure nephograms of different methods: (a) SPH method, (b) AUTODYN.
(a)
(b)Figure 4
The pressure-time curves of test point A.The pressure of test point A in the cases of different initial particle spacing is shown in Figure5. We note that the ratio of initial particle spacing between water and shell is n and discuss the different cases (resp., n = 1 : 1, 2 : 1, 3 : 1, and 4 : 1). The following figures are pressure-time curves of test point A in different cases and the verified model. As Figure 5 shows, the pressure-time curves in cases where n = 1 : 1 and n = 2 : 1 are basically consistent, but the case where n = 1 : 1 requires a significantly larger amount of computation than the case where n = 2 : 1. However, there are significant pressure noises at the interface in cases where n = 3 : 1 and n = 4 : 1 and serious numerical oscillation occurs when n = 4 : 1. Therefore, n = 2 : 1 is chosen.Figure 5
The pressure-time curves of test point A in the cases with different initial particle spacing;n is the ratio of initial particle spacing between water and shell.


### 3.3. Results and Discussion

#### 3.3.1. Shockwave Propagation

After TNT detonation, the shockwaves propagates and overlaps in multilayer media, which will cause the inhomogeneity in time and space. According to the principle of impedance matching, the results above are in good accordance with the physics law, and the numerical results are fairing and smooth, which verifies that the dynamic and motion continuity conditions are satisfied in the SPH method.For the single cylindrical shell, the pressure of point A is shown in Figure6. It is obvious that the peak pressure of the out-covering is smaller than those of the other two cases. The reason may be that a rarefaction wave is reflected when shock wave transmits from TNT of high-impedance to rubber of low-impedance, which will result in low pressure.Figure 6
The pressure of test point A; Cases1
–
3 correspond to no-covering, in-covering, and out-covering.The pressure of the in-covering with different covering thickness is shown in Figure7. Define the dimensionless thickness C
C as the ratio of covering thickness C and the thickness of steel. When C
C
<
1, the so-called “double peak” or “multimodality” of pressure is observed. With the increase of C
C, the peak pressure goes up. When C
C
=
1, it is up to the maximum, and does not increase with the increase of covering thickness.Figure 7
The pressure of test point A in the cases of in-covering with different covering thicknesses;C is covering thickness.


#### 3.3.2. Damage Characteristics of Single-Layer Cylindrical Shell

Damage Characteristics of Different Covering Positions. The failure modes of shell subjected to underwater explosion can be divided into three modes [30]:Mode 1: plastic large-deformation;
Mode 2: tensile tearing in outer fibers at the support;
Mode 3: transverse shear failure at the support.The equivalent plastic strain of coated shell is shown in Figure8. As for (a), it is observed that the damage process experiences plug failure, dent, and rolling. As for (b), the shell firstly develops a pit with a greater plastic deformation and becomes ruptured at about 0.9 ms. The size of the fragment is close to the contact area of TNT. As for (c), rupture occurs at the center of the shell first. With the increase of deflection, the tensile failure arises. Consequently, it is concluded that the main failure mode of no-covering shell is plastic large-deformation and tensile failure, while the mode will be tension failure and shear failure if the shell has coverings. The size of fragment is related to the contact area. In addition, the differences between in-covering and out-covering are initial rupture positions and initial rupture time.The equivalent plastic strain of coated shell; (a), (b), and (c) correspond to no-covering, in-covering, and out-covering, respectively.
(a)
(b)
(c)Damage Characteristics of Different Covering Thicknesses. The equivalent plastic strain of different covering thicknesses is shown in Figure 9. For (a), (b), and (c) where C
C = 0, 0.5 and 0.75, respectively, the main failure mode is plastic large-deformation and tensile failure, and the crevasse tends to roll; as for (d) where C
C
=
1.0, the tension failure and shear failure should be the main failure modes, and the size of fragment is close to the contact area of TNT; as for (e) and (f) where C
C = 1.25 and 1.5, which are similar to C
C
=
1.0, the ruptures come out at different positions. In a word, when the covering thickness is close to that of steel, the protective effect is better.The equivalent plastic strain of different covering thicknesses at about 2 ms; (a)–(f) correspond to the covering thickness of 0, 1.0 cm, 1.5 cm, 2.0 cm, 2.5 cm, and 3.0 cm, respectively.
(a)
(b)
(c)
(d)
(e)
(f)Crevasse Radius. The crevasse radius of different coated shells is shown in Figure 10. The dimensionless radius K is defined as the ratio of crevasse radius to the TNT radius. As for Case  1, the crevasse radius increases through time, on the contrast, there is a little rebound process in Case  2 and Case  3. The turning point of Case  3 further proves that the main failure mode of out-covering is shear failure. In brief, the covering may not always take positive effects on the shell, better protective effect is found in the case of in-covering.Figure 10
The crevasse radius of different coated shells; Cases1
–
3 correspond to the no-covering, in-covering, and out-covering.The crevasse radius of different covering thicknesses is shown in Figure11. When C
C
≤
0.75, the crevasse radius increases with time first, and decreases with the increase of covering thickness. When C
C
=
1.0, the crevasse radius undergoes a little rebound after the peak. When C
C = 1.25, 1.5, a turning point turns out which also proves that shear failure is one of the main failure modes; in addition, the thicker the covering is, the larger K will be. To sum up, the better protective effect appears when the covering thickness is equal to shell thickness.Figure 11
The crevasse radius of different covering thicknesses;C is covering thickness.


## 3.1. Numerical Model and Parameters

The model of a coated single-hull submarine subjected to the attack of torpedo can be simplified as Figure2, and (a) indicates no-covering, (b) denotes in-covering, and (c) represents out-covering. The relevant parameters and cases are listed in Tables 7 and 8. The origin of coordinates is positioned at the center of ring-shaped water. The column TNT is detonated from the center. The TNT dosage is about 12.8 kg, a popular amount in most model tests. The entire model is discretized with nonuniform particle spacing, 0.004 m for the cylindrical shell and 0.008 m for TNT and water, and some verification for the initial particle spacing is explained in Section 3.2. The numbers of different kinds of particles are listed as Table 8. The test point A is placed at (0, 2.1, 0).Table 7
Parameters of single cylindrical shell.
Outside radius of water (m)
Inside radius of water (m)
Inside radius of shell (m)
Thickness of shell (m)
Radius of TNT (m)
3.0
1.22
1.2
0.02
0.05Table 8
Different cases of single cylindrical shell and numbers of particles (different positions are listed as Cases  1–3, in the case where the rubber thickness is 0.020 m; different thicknesses are listed as Cases  4–7, in the case where the inner surface is coated with rubber).
Cases
Name
Numbers of particles
Water
Steel
TNT
Rubber
Total
Case  1
No-covering
364507
9300
121
0
373928
Case  2
In-covering
364507
9300
121
9153
383081
Case  3
Out-covering
355250
9300
121
9257
373928
Case  4
0.010 m
364507
9300
121
13680
387608
Case  5
0.015 m
364507
9300
121
11406
385334
Case  6
0.025 m
364507
9300
121
6886
380814
Case  7
0.030 m
364507
9300
121
4575
378503Models of single-hull.
(a)
(b)
(c)


## 3.2. Numerical Verification

In order to verify the effectiveness of the present SPH method, the no-covering model of underwater contact explosion is selected, the SPH result will be compared with that of AUTODYN.The propagations of shockwave obtained from different methods are shown in Figure3 and the pressure of test point A is shown in Figure 4. It is obvious that the propagation laws of shockwave agree well with each other; moreover, the curves obtained from different methods also show good agreements, which all prove the effectiveness and accuracy of the present SPH method to solve the problems of underwater contact explosion.Pressure nephograms of different methods: (a) SPH method, (b) AUTODYN.
(a)
(b)Figure 4
The pressure-time curves of test point A.The pressure of test point A in the cases of different initial particle spacing is shown in Figure5. We note that the ratio of initial particle spacing between water and shell is n and discuss the different cases (resp., n = 1 : 1, 2 : 1, 3 : 1, and 4 : 1). The following figures are pressure-time curves of test point A in different cases and the verified model. As Figure 5 shows, the pressure-time curves in cases where n = 1 : 1 and n = 2 : 1 are basically consistent, but the case where n = 1 : 1 requires a significantly larger amount of computation than the case where n = 2 : 1. However, there are significant pressure noises at the interface in cases where n = 3 : 1 and n = 4 : 1 and serious numerical oscillation occurs when n = 4 : 1. Therefore, n = 2 : 1 is chosen.Figure 5
The pressure-time curves of test point A in the cases with different initial particle spacing;n is the ratio of initial particle spacing between water and shell.


## 3.3. Results and Discussion

### 3.3.1. Shockwave Propagation

After TNT detonation, the shockwaves propagates and overlaps in multilayer media, which will cause the inhomogeneity in time and space. According to the principle of impedance matching, the results above are in good accordance with the physics law, and the numerical results are fairing and smooth, which verifies that the dynamic and motion continuity conditions are satisfied in the SPH method.For the single cylindrical shell, the pressure of point A is shown in Figure6. It is obvious that the peak pressure of the out-covering is smaller than those of the other two cases. The reason may be that a rarefaction wave is reflected when shock wave transmits from TNT of high-impedance to rubber of low-impedance, which will result in low pressure.Figure 6
The pressure of test point A; Cases1
–
3 correspond to no-covering, in-covering, and out-covering.The pressure of the in-covering with different covering thickness is shown in Figure7. Define the dimensionless thickness C
C as the ratio of covering thickness C and the thickness of steel. When C
C
<
1, the so-called “double peak” or “multimodality” of pressure is observed. With the increase of C
C, the peak pressure goes up. When C
C
=
1, it is up to the maximum, and does not increase with the increase of covering thickness.Figure 7
The pressure of test point A in the cases of in-covering with different covering thicknesses;C is covering thickness.


### 3.3.2. Damage Characteristics of Single-Layer Cylindrical Shell

Damage Characteristics of Different Covering Positions. The failure modes of shell subjected to underwater explosion can be divided into three modes [30]:Mode 1: plastic large-deformation;
Mode 2: tensile tearing in outer fibers at the support;
Mode 3: transverse shear failure at the support.The equivalent plastic strain of coated shell is shown in Figure8. As for (a), it is observed that the damage process experiences plug failure, dent, and rolling. As for (b), the shell firstly develops a pit with a greater plastic deformation and becomes ruptured at about 0.9 ms. The size of the fragment is close to the contact area of TNT. As for (c), rupture occurs at the center of the shell first. With the increase of deflection, the tensile failure arises. Consequently, it is concluded that the main failure mode of no-covering shell is plastic large-deformation and tensile failure, while the mode will be tension failure and shear failure if the shell has coverings. The size of fragment is related to the contact area. In addition, the differences between in-covering and out-covering are initial rupture positions and initial rupture time.The equivalent plastic strain of coated shell; (a), (b), and (c) correspond to no-covering, in-covering, and out-covering, respectively.
(a)
(b)
(c)Damage Characteristics of Different Covering Thicknesses. The equivalent plastic strain of different covering thicknesses is shown in Figure 9. For (a), (b), and (c) where C
C = 0, 0.5 and 0.75, respectively, the main failure mode is plastic large-deformation and tensile failure, and the crevasse tends to roll; as for (d) where C
C
=
1.0, the tension failure and shear failure should be the main failure modes, and the size of fragment is close to the contact area of TNT; as for (e) and (f) where C
C = 1.25 and 1.5, which are similar to C
C
=
1.0, the ruptures come out at different positions. In a word, when the covering thickness is close to that of steel, the protective effect is better.The equivalent plastic strain of different covering thicknesses at about 2 ms; (a)–(f) correspond to the covering thickness of 0, 1.0 cm, 1.5 cm, 2.0 cm, 2.5 cm, and 3.0 cm, respectively.
(a)
(b)
(c)
(d)
(e)
(f)Crevasse Radius. The crevasse radius of different coated shells is shown in Figure 10. The dimensionless radius K is defined as the ratio of crevasse radius to the TNT radius. As for Case  1, the crevasse radius increases through time, on the contrast, there is a little rebound process in Case  2 and Case  3. The turning point of Case  3 further proves that the main failure mode of out-covering is shear failure. In brief, the covering may not always take positive effects on the shell, better protective effect is found in the case of in-covering.Figure 10
The crevasse radius of different coated shells; Cases1
–
3 correspond to the no-covering, in-covering, and out-covering.The crevasse radius of different covering thicknesses is shown in Figure11. When C
C
≤
0.75, the crevasse radius increases with time first, and decreases with the increase of covering thickness. When C
C
=
1.0, the crevasse radius undergoes a little rebound after the peak. When C
C = 1.25, 1.5, a turning point turns out which also proves that shear failure is one of the main failure modes; in addition, the thicker the covering is, the larger K will be. To sum up, the better protective effect appears when the covering thickness is equal to shell thickness.Figure 11
The crevasse radius of different covering thicknesses;C is covering thickness.


## 3.3.1. Shockwave Propagation

After TNT detonation, the shockwaves propagates and overlaps in multilayer media, which will cause the inhomogeneity in time and space. According to the principle of impedance matching, the results above are in good accordance with the physics law, and the numerical results are fairing and smooth, which verifies that the dynamic and motion continuity conditions are satisfied in the SPH method.For the single cylindrical shell, the pressure of point A is shown in Figure6. It is obvious that the peak pressure of the out-covering is smaller than those of the other two cases. The reason may be that a rarefaction wave is reflected when shock wave transmits from TNT of high-impedance to rubber of low-impedance, which will result in low pressure.Figure 6
The pressure of test point A; Cases1
–
3 correspond to no-covering, in-covering, and out-covering.The pressure of the in-covering with different covering thickness is shown in Figure7. Define the dimensionless thickness C
C as the ratio of covering thickness C and the thickness of steel. When C
C
<
1, the so-called “double peak” or “multimodality” of pressure is observed. With the increase of C
C, the peak pressure goes up. When C
C
=
1, it is up to the maximum, and does not increase with the increase of covering thickness.Figure 7
The pressure of test point A in the cases of in-covering with different covering thicknesses;C is covering thickness.


## 3.3.2. Damage Characteristics of Single-Layer Cylindrical Shell

Damage Characteristics of Different Covering Positions. The failure modes of shell subjected to underwater explosion can be divided into three modes [30]:Mode 1: plastic large-deformation;
Mode 2: tensile tearing in outer fibers at the support;
Mode 3: transverse shear failure at the support.The equivalent plastic strain of coated shell is shown in Figure8. As for (a), it is observed that the damage process experiences plug failure, dent, and rolling. As for (b), the shell firstly develops a pit with a greater plastic deformation and becomes ruptured at about 0.9 ms. The size of the fragment is close to the contact area of TNT. As for (c), rupture occurs at the center of the shell first. With the increase of deflection, the tensile failure arises. Consequently, it is concluded that the main failure mode of no-covering shell is plastic large-deformation and tensile failure, while the mode will be tension failure and shear failure if the shell has coverings. The size of fragment is related to the contact area. In addition, the differences between in-covering and out-covering are initial rupture positions and initial rupture time.The equivalent plastic strain of coated shell; (a), (b), and (c) correspond to no-covering, in-covering, and out-covering, respectively.
(a)
(b)
(c)Damage Characteristics of Different Covering Thicknesses. The equivalent plastic strain of different covering thicknesses is shown in Figure 9. For (a), (b), and (c) where C
C = 0, 0.5 and 0.75, respectively, the main failure mode is plastic large-deformation and tensile failure, and the crevasse tends to roll; as for (d) where C
C
=
1.0, the tension failure and shear failure should be the main failure modes, and the size of fragment is close to the contact area of TNT; as for (e) and (f) where C
C = 1.25 and 1.5, which are similar to C
C
=
1.0, the ruptures come out at different positions. In a word, when the covering thickness is close to that of steel, the protective effect is better.The equivalent plastic strain of different covering thicknesses at about 2 ms; (a)–(f) correspond to the covering thickness of 0, 1.0 cm, 1.5 cm, 2.0 cm, 2.5 cm, and 3.0 cm, respectively.
(a)
(b)
(c)
(d)
(e)
(f)Crevasse Radius. The crevasse radius of different coated shells is shown in Figure 10. The dimensionless radius K is defined as the ratio of crevasse radius to the TNT radius. As for Case  1, the crevasse radius increases through time, on the contrast, there is a little rebound process in Case  2 and Case  3. The turning point of Case  3 further proves that the main failure mode of out-covering is shear failure. In brief, the covering may not always take positive effects on the shell, better protective effect is found in the case of in-covering.Figure 10
The crevasse radius of different coated shells; Cases1
–
3 correspond to the no-covering, in-covering, and out-covering.The crevasse radius of different covering thicknesses is shown in Figure11. When C
C
≤
0.75, the crevasse radius increases with time first, and decreases with the increase of covering thickness. When C
C
=
1.0, the crevasse radius undergoes a little rebound after the peak. When C
C = 1.25, 1.5, a turning point turns out which also proves that shear failure is one of the main failure modes; in addition, the thicker the covering is, the larger K will be. To sum up, the better protective effect appears when the covering thickness is equal to shell thickness.Figure 11
The crevasse radius of different covering thicknesses;C is covering thickness.


## 4. Contact Explosion of Double Cylindrical Shell

### 4.1. Numerical Model and Parameters

Based on the researches of single cylindrical shell, the influences of covering position and interlayer water on the damage characteristics of the double cylindrical shell are studied. The basic model is shown in Figure12; the detailed parameters and setting cases are listed in Tables 9 and 10. The column TNT is detonated from the center. The TNT dosage is about 12.8 kg, a popular amount in most model tests. The origin of coordinates is also positioned at the center of the model. The two test points are A (0, 2.1, 0) and B (0, 0.73, 0). In case (a), there is no interlayer water while case (b) has full interlayer water. The circular region represents the main research object. The whole model is discretized with nonuniform particle spacing, 0.003 m for the cylindrical shell and 0.006 m for TNT and water. And the numbers of different kinds of particles are listed in Table 10.Table 9
Parameters of double cylindrical shell.
Outside radius of water (m)
Inside radius of water (m)
Outside radius of shell  (m)
Thickness of shell (m)
Radius of TNT (m)
Inner shell
Outer shell
Inner shell
Outer shell
3.0
1.22
0.865
1.22
0.02
0.015
0.05Table 10
Different cases of double cylindrical shell (Cases  1–5 belong to the case of double-hull without interlayer water; Cases  6–10 belong to the case of double-hull full of interlayer water).
Cases
Name
Numbers of particles
Water
Steel
TNT
Rubber
Total
Case  1
No-covering
648852
36298
212
0
685362
Case  2
In-covering on inner shell
648852
36298
212
11632
696994
Case  3
Out-covering on inner shell
648852
36298
212
11795
697157
Case  4
In-covering on outer shell
648852
36298
212
16269
701631
Case  5
Out-covering on outer shell
644625
36298
212
16821
697956
Case  6
No-covering
709479
36298
212
0
745989
Case  7
In-covering on inner shell
709479
36298
212
11632
757621
Case  8
Out-covering on inner shell
706401
36298
212
11795
754706
Case  9
In-covering on outer shell
705143
36298
212
16269
757922
Case  10
Out-covering on outer shell
705127
36298
212
16821
758458Models of double-hull.
(a)
(b)


### 4.2. Results and Discussion

#### 4.2.1. Propagation of Shockwaves

The propagation of shockwaves in the case of in-covering on inner surface is shown in Figure13. Shock wave A
s propagates from water to inner shell at t
≈
0.20 ms; at the same time, a reflected shock wave B
s is generated on the lower surface of inner shell. When arriving at outer shell, a reflected shock wave C
s is generated on the upper surface of outer shell, which is shown in the rectangular region. At t
≈
0.58 ms, C
s not only reflects a shock wave D
s on the upper surface of outer shell (shown in the rectangular region) but also generates a shock wave E
s (shown in the elliptical region) transmitting to the water; the pressure in the circular region is larger because of the superposition of A
s, B
s, C
s, and D
s. In the process, it is obvious that the shock wave velocity in steel is larger than that in water or rubber, yet the velocity in water is close to that in rubber because of similar impedances. The impedance mismatch causes the complexity of shock wave propagation in different media.Pressure nephograms of double-hull.
(a)
(b)
(c)The pressure at test point A is shown in Figure14. As for 1, due to the air-backed outer shell, a rarefaction wave is generated, which will decrease the pressure at test point A. However, as for 2, the rarefaction wave takes little effects on the test point. So the peak pressure of 2 is larger than that of 1. In addition, as for 2, the second peak pressure appears at 0.9 ms because there is a reflected shock wave transmitting from the inner shell to water.Figure 14
The pressure of test point A; 1 and 2 mean no interlayer water and full interlayer water; cases are shown as Table10.The pressure at test point B in the case of full interlayer water is shown in Figure15. It turns out that the most notable feature is multimodality and the peak pressure decays with time. It is not beyond imagination that the superposition of incident wave, reflected wave, and transmitted wave appears at inner shell, outer shell, and the water between them.Figure 15
The pressure of test point B; 2 means full interlayer water; cases are shown in Table10.


#### 4.2.2. Damage Characteristics of Single Cylindrical Shell

Crevasse Mode. Failure modes 1–3 are plastic large-deformation, tensile failure, and shear failure, respectively. Moreover, the shear failure can also be classified to three categories 3a–3c. The failure modes of double cylindrical shell coated on different surfaces are shown in Figures 16 and 17.Failure modes of different coated shells at about 2 ms; (a)–(d) mean no-covering (no interlayer water), no-covering (full interlayer water), in-covering of inner shell (full interlayer water), and out-covering of inner shell (full interlayer water).
(a)
(b)
(c)
(d)Figure 17
Failure modes of different coated shells; Cases  1–10 are shown in Table10.As for the outer shell, the damage process experiences plug failure, dent, and rolling. As for the inner shell, in the case of no interlayer water and no-covering, the failure of inner shell is caused by high-speed fragments from outer shell, which is called “3a,” as shown in Figure16(a); in the case of full interlayer water and no-covering, the crack appears at an angle of 45° from the normal direction, which is called “3c,” as shown in Figure 16(b); in the case of full interlayer water and in-covering of inner shell, only plastic large-deformation comes out, as shown in Figure 16(c); in the case of full interlayer water and out-covering of inner shell, the crack appears in the center of inner shell in normal direction, which is called “3b,” as shown in Figure 16(d). To sum up, when the interlayer water is full, the in-covering of the inner shell will better protect the inner shell and the interlayer water can reduce the damage caused by fragments.Crevasse Radius. The crevasse variation of the outer shell without interlayer water is shown in Figure 18. In case 1, without covering, the crevasse increases as time goes on. However, in other cases, there is a little rebound when up to a peak and tends to be steady finally. Obviously, the out-covering provides better protection to the outer shell.Figure 18
The crevasse variation of outer shell without interlayer water; cases are shown in Table10.The crevasse variation of the outer shell with full interlayer water is shown in Figure19. In cases  6, 7, and 8, the crevasses increase through time, and the coincidence indicates that the protective effect on outer shell is not obvious when inner shell is coated. Nonetheless, in cases 9 and 10, the crevasse is smaller, which shows that the covering on the outer shell is beneficial for the outer shell. In addition, the water has certain influences on the outer shell because the water would participate in energy absorption.Figure 19
The crevasse variation of outer shell with full interlayer water; cases are shown in Table10.The crevasse variation of the inner shell without interlayer water is shown in Figure20. The turning points account for the fracture of the inner shell. In case  9 with in-covering on the inner shell, the crevasse is larger than case  8 without covering, which indicates that the covering does not always take positive effects.Figure 20
The crevasse variation of inner shell without interlayer water; cases are shown in Table10.The crevasse variation of the inner shell with full interlayer water is shown in Figure21. Except cases  9 and 10, the curves appear to overlap, which indicates that the protection is not obvious. It seems that the effect of in-covering on the inner shell is better; however, considering the operation in practice, out-covering on inner shell would be a practical choice. In addition, comparing the failure modes and the sizes of crevasse between cases with/without interlayer water, it is found that the inner shell can be protected by the water from the damage caused by fragments.Figure 21
The crevasse variation of inner shell with full interlayer water; cases are shown in Table10.


## 4.1. Numerical Model and Parameters

Based on the researches of single cylindrical shell, the influences of covering position and interlayer water on the damage characteristics of the double cylindrical shell are studied. The basic model is shown in Figure12; the detailed parameters and setting cases are listed in Tables 9 and 10. The column TNT is detonated from the center. The TNT dosage is about 12.8 kg, a popular amount in most model tests. The origin of coordinates is also positioned at the center of the model. The two test points are A (0, 2.1, 0) and B (0, 0.73, 0). In case (a), there is no interlayer water while case (b) has full interlayer water. The circular region represents the main research object. The whole model is discretized with nonuniform particle spacing, 0.003 m for the cylindrical shell and 0.006 m for TNT and water. And the numbers of different kinds of particles are listed in Table 10.Table 9
Parameters of double cylindrical shell.
Outside radius of water (m)
Inside radius of water (m)
Outside radius of shell  (m)
Thickness of shell (m)
Radius of TNT (m)
Inner shell
Outer shell
Inner shell
Outer shell
3.0
1.22
0.865
1.22
0.02
0.015
0.05Table 10
Different cases of double cylindrical shell (Cases  1–5 belong to the case of double-hull without interlayer water; Cases  6–10 belong to the case of double-hull full of interlayer water).
Cases
Name
Numbers of particles
Water
Steel
TNT
Rubber
Total
Case  1
No-covering
648852
36298
212
0
685362
Case  2
In-covering on inner shell
648852
36298
212
11632
696994
Case  3
Out-covering on inner shell
648852
36298
212
11795
697157
Case  4
In-covering on outer shell
648852
36298
212
16269
701631
Case  5
Out-covering on outer shell
644625
36298
212
16821
697956
Case  6
No-covering
709479
36298
212
0
745989
Case  7
In-covering on inner shell
709479
36298
212
11632
757621
Case  8
Out-covering on inner shell
706401
36298
212
11795
754706
Case  9
In-covering on outer shell
705143
36298
212
16269
757922
Case  10
Out-covering on outer shell
705127
36298
212
16821
758458Models of double-hull.
(a)
(b)


## 4.2. Results and Discussion

### 4.2.1. Propagation of Shockwaves

The propagation of shockwaves in the case of in-covering on inner surface is shown in Figure13. Shock wave A
s propagates from water to inner shell at t
≈
0.20 ms; at the same time, a reflected shock wave B
s is generated on the lower surface of inner shell. When arriving at outer shell, a reflected shock wave C
s is generated on the upper surface of outer shell, which is shown in the rectangular region. At t
≈
0.58 ms, C
s not only reflects a shock wave D
s on the upper surface of outer shell (shown in the rectangular region) but also generates a shock wave E
s (shown in the elliptical region) transmitting to the water; the pressure in the circular region is larger because of the superposition of A
s, B
s, C
s, and D
s. In the process, it is obvious that the shock wave velocity in steel is larger than that in water or rubber, yet the velocity in water is close to that in rubber because of similar impedances. The impedance mismatch causes the complexity of shock wave propagation in different media.Pressure nephograms of double-hull.
(a)
(b)
(c)The pressure at test point A is shown in Figure14. As for 1, due to the air-backed outer shell, a rarefaction wave is generated, which will decrease the pressure at test point A. However, as for 2, the rarefaction wave takes little effects on the test point. So the peak pressure of 2 is larger than that of 1. In addition, as for 2, the second peak pressure appears at 0.9 ms because there is a reflected shock wave transmitting from the inner shell to water.Figure 14
The pressure of test point A; 1 and 2 mean no interlayer water and full interlayer water; cases are shown as Table10.The pressure at test point B in the case of full interlayer water is shown in Figure15. It turns out that the most notable feature is multimodality and the peak pressure decays with time. It is not beyond imagination that the superposition of incident wave, reflected wave, and transmitted wave appears at inner shell, outer shell, and the water between them.Figure 15
The pressure of test point B; 2 means full interlayer water; cases are shown in Table10.


### 4.2.2. Damage Characteristics of Single Cylindrical Shell

Crevasse Mode. Failure modes 1–3 are plastic large-deformation, tensile failure, and shear failure, respectively. Moreover, the shear failure can also be classified to three categories 3a–3c. The failure modes of double cylindrical shell coated on different surfaces are shown in Figures 16 and 17.Failure modes of different coated shells at about 2 ms; (a)–(d) mean no-covering (no interlayer water), no-covering (full interlayer water), in-covering of inner shell (full interlayer water), and out-covering of inner shell (full interlayer water).
(a)
(b)
(c)
(d)Figure 17
Failure modes of different coated shells; Cases  1–10 are shown in Table10.As for the outer shell, the damage process experiences plug failure, dent, and rolling. As for the inner shell, in the case of no interlayer water and no-covering, the failure of inner shell is caused by high-speed fragments from outer shell, which is called “3a,” as shown in Figure16(a); in the case of full interlayer water and no-covering, the crack appears at an angle of 45° from the normal direction, which is called “3c,” as shown in Figure 16(b); in the case of full interlayer water and in-covering of inner shell, only plastic large-deformation comes out, as shown in Figure 16(c); in the case of full interlayer water and out-covering of inner shell, the crack appears in the center of inner shell in normal direction, which is called “3b,” as shown in Figure 16(d). To sum up, when the interlayer water is full, the in-covering of the inner shell will better protect the inner shell and the interlayer water can reduce the damage caused by fragments.Crevasse Radius. The crevasse variation of the outer shell without interlayer water is shown in Figure 18. In case 1, without covering, the crevasse increases as time goes on. However, in other cases, there is a little rebound when up to a peak and tends to be steady finally. Obviously, the out-covering provides better protection to the outer shell.Figure 18
The crevasse variation of outer shell without interlayer water; cases are shown in Table10.The crevasse variation of the outer shell with full interlayer water is shown in Figure19. In cases  6, 7, and 8, the crevasses increase through time, and the coincidence indicates that the protective effect on outer shell is not obvious when inner shell is coated. Nonetheless, in cases 9 and 10, the crevasse is smaller, which shows that the covering on the outer shell is beneficial for the outer shell. In addition, the water has certain influences on the outer shell because the water would participate in energy absorption.Figure 19
The crevasse variation of outer shell with full interlayer water; cases are shown in Table10.The crevasse variation of the inner shell without interlayer water is shown in Figure20. The turning points account for the fracture of the inner shell. In case  9 with in-covering on the inner shell, the crevasse is larger than case  8 without covering, which indicates that the covering does not always take positive effects.Figure 20
The crevasse variation of inner shell without interlayer water; cases are shown in Table10.The crevasse variation of the inner shell with full interlayer water is shown in Figure21. Except cases  9 and 10, the curves appear to overlap, which indicates that the protection is not obvious. It seems that the effect of in-covering on the inner shell is better; however, considering the operation in practice, out-covering on inner shell would be a practical choice. In addition, comparing the failure modes and the sizes of crevasse between cases with/without interlayer water, it is found that the inner shell can be protected by the water from the damage caused by fragments.Figure 21
The crevasse variation of inner shell with full interlayer water; cases are shown in Table10.


## 4.2.1. Propagation of Shockwaves

The propagation of shockwaves in the case of in-covering on inner surface is shown in Figure13. Shock wave A
s propagates from water to inner shell at t
≈
0.20 ms; at the same time, a reflected shock wave B
s is generated on the lower surface of inner shell. When arriving at outer shell, a reflected shock wave C
s is generated on the upper surface of outer shell, which is shown in the rectangular region. At t
≈
0.58 ms, C
s not only reflects a shock wave D
s on the upper surface of outer shell (shown in the rectangular region) but also generates a shock wave E
s (shown in the elliptical region) transmitting to the water; the pressure in the circular region is larger because of the superposition of A
s, B
s, C
s, and D
s. In the process, it is obvious that the shock wave velocity in steel is larger than that in water or rubber, yet the velocity in water is close to that in rubber because of similar impedances. The impedance mismatch causes the complexity of shock wave propagation in different media.Pressure nephograms of double-hull.
(a)
(b)
(c)The pressure at test point A is shown in Figure14. As for 1, due to the air-backed outer shell, a rarefaction wave is generated, which will decrease the pressure at test point A. However, as for 2, the rarefaction wave takes little effects on the test point. So the peak pressure of 2 is larger than that of 1. In addition, as for 2, the second peak pressure appears at 0.9 ms because there is a reflected shock wave transmitting from the inner shell to water.Figure 14
The pressure of test point A; 1 and 2 mean no interlayer water and full interlayer water; cases are shown as Table10.The pressure at test point B in the case of full interlayer water is shown in Figure15. It turns out that the most notable feature is multimodality and the peak pressure decays with time. It is not beyond imagination that the superposition of incident wave, reflected wave, and transmitted wave appears at inner shell, outer shell, and the water between them.Figure 15
The pressure of test point B; 2 means full interlayer water; cases are shown in Table10.


## 4.2.2. Damage Characteristics of Single Cylindrical Shell

Crevasse Mode. Failure modes 1–3 are plastic large-deformation, tensile failure, and shear failure, respectively. Moreover, the shear failure can also be classified to three categories 3a–3c. The failure modes of double cylindrical shell coated on different surfaces are shown in Figures 16 and 17.Failure modes of different coated shells at about 2 ms; (a)–(d) mean no-covering (no interlayer water), no-covering (full interlayer water), in-covering of inner shell (full interlayer water), and out-covering of inner shell (full interlayer water).
(a)
(b)
(c)
(d)Figure 17
Failure modes of different coated shells; Cases  1–10 are shown in Table10.As for the outer shell, the damage process experiences plug failure, dent, and rolling. As for the inner shell, in the case of no interlayer water and no-covering, the failure of inner shell is caused by high-speed fragments from outer shell, which is called “3a,” as shown in Figure16(a); in the case of full interlayer water and no-covering, the crack appears at an angle of 45° from the normal direction, which is called “3c,” as shown in Figure 16(b); in the case of full interlayer water and in-covering of inner shell, only plastic large-deformation comes out, as shown in Figure 16(c); in the case of full interlayer water and out-covering of inner shell, the crack appears in the center of inner shell in normal direction, which is called “3b,” as shown in Figure 16(d). To sum up, when the interlayer water is full, the in-covering of the inner shell will better protect the inner shell and the interlayer water can reduce the damage caused by fragments.Crevasse Radius. The crevasse variation of the outer shell without interlayer water is shown in Figure 18. In case 1, without covering, the crevasse increases as time goes on. However, in other cases, there is a little rebound when up to a peak and tends to be steady finally. Obviously, the out-covering provides better protection to the outer shell.Figure 18
The crevasse variation of outer shell without interlayer water; cases are shown in Table10.The crevasse variation of the outer shell with full interlayer water is shown in Figure19. In cases  6, 7, and 8, the crevasses increase through time, and the coincidence indicates that the protective effect on outer shell is not obvious when inner shell is coated. Nonetheless, in cases 9 and 10, the crevasse is smaller, which shows that the covering on the outer shell is beneficial for the outer shell. In addition, the water has certain influences on the outer shell because the water would participate in energy absorption.Figure 19
The crevasse variation of outer shell with full interlayer water; cases are shown in Table10.The crevasse variation of the inner shell without interlayer water is shown in Figure20. The turning points account for the fracture of the inner shell. In case  9 with in-covering on the inner shell, the crevasse is larger than case  8 without covering, which indicates that the covering does not always take positive effects.Figure 20
The crevasse variation of inner shell without interlayer water; cases are shown in Table10.The crevasse variation of the inner shell with full interlayer water is shown in Figure21. Except cases  9 and 10, the curves appear to overlap, which indicates that the protection is not obvious. It seems that the effect of in-covering on the inner shell is better; however, considering the operation in practice, out-covering on inner shell would be a practical choice. In addition, comparing the failure modes and the sizes of crevasse between cases with/without interlayer water, it is found that the inner shell can be protected by the water from the damage caused by fragments.Figure 21
The crevasse variation of inner shell with full interlayer water; cases are shown in Table10.


## 5. Conclusions

An SPH method with mesh-free and Lagrange properties is applied to solve extremely dynamic problems of cylindrical shell subjected to underwater contact explosion in this paper. The influences of coverings on damage characteristics of cylindrical shells are investigated. Through the parametric analyses of covering position and covering thickness, the following conclusions can be drawn.(1)
The propagation laws of shockwaves in multilayer media correspond with the theory of impedance matching, and the continuity conditions of discontinuity surface are verified by the smooth SPH results; the results of SPH agree well with those of AUTODYN, which verifies the feasibility and effectiveness of the present SPH method.
(2)
For single cylindrical shell, the failure mode will be large-deformation and tensile failure when the covering is thinner than the shell, and tension failure and shear failure otherwise. The antishock performance does not necessarily increase with covering thickness; better protective effects appear when the covering thickness is close to that of the shell.
(3)
As for double cylindrical shell, the obvious protective effect is found when there is out-covering on the outer shell or in-covering on the inner shell. However, considering the operability, the out-covering on inner shell is also a good choice when the interlayer is full of water.
(4)
The damage of single shell is mainly caused by detonation wave and detonation products; as for double shell without interlayer water, the inner shell is mainly disrupted by high-speed fragments. Moreover, a certain protective effect will be provided by the interlayer water.



---

*Source: 763607-2014-02-24.xml*