Jtam.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

42, 1, pp. 125-138, Warsaw 2004

COMPUTER SIMULATION OF CATCHING OF

EXPLOSIVELY PROPELLING METAL FRAGMENTS BY

PROTECTIVE CASING

Karol Jach
Robert Świerczyński

Institute of Optoelectronics, Military University of Technology, Warsaw

e-mail: kjach@wat.edu.pl

Adam Wiśniewski

Military Institute of Armament Technology, Zielonka

e-mail: zak 18@witu.mil.pl

In the paper a mathematical-physical model, results of computer simu-
lation of the propelling of metal fragments by pressure of detonation
products and their catching by the protective casing made of a quasi-
aramide fabric havebeen presented.The impact on the protective casing
due to shock waves generated by the explosion has been examined. De-
pending on the nature of the issue, the 2D free particles or 1DLagrange
method has been used in the computer simulation.

Key words: computer simulation, mutal fragments, protective casing

1. Explosive propelling of metal fragments

In order to allow safe detonation of a terrorist bomb, located in a public
place, different kinds of the protective casingmade of steel, gum, etc. are used.

One of the important problems that must be considered before the con-
struction of such a protective layer is gaining theoretical knowledge about
behaviour of the debris during the explosion of the bomb. In this paper, si-
mulation of the explosion has been described as the impact of a fragment
propelled by explosion products of a 200g trinitrotoluene block on the pro-
tective casing. The influence of the air shock wave caused by the explosion
on the casing and the time of flight of the fragment and of this wave to the



126 K.Jach et al.

internal surface of the protective layer has been modelled as well. Its is espe-
cially important to evaluate the power of the destruction of the casing by the
air shock wave during the flight of the fragment to the internal surface of the
protective layer. It has been assumed that the protective casing was made of
a quasi-aramide fabric.
Because of the lack of literaturedata concerning thebehaviour of thequasi-

aramid fabricunderdynamic loading, computer simulationhasbeenperformed
in theway to obtain, on the base of the knowledge of static parameters of this
fabric andmultiple computer simulations, sufficientmatching of experimental
results with the results of this simulation. Such a method of computer analy-
sis permits one to choose constants of the equation of state and constitutive
relations of the quasi-aramid fabric, and then to proceed with computer si-
mulations with varying input data of the casing in order to optimise its mass,
etc.
To describe the behaviour of metals under high dynamic loading during

explosive propelling of fragments, amodel of elastoplastic bodyhas beenused.
The system of equations expressing the conservation laws and constitutive
relations for thismodel have the following form (axial symmetry), seeWilkins
(1984), Jach et al. (2001), Jach andWłodarczyk (1992)

dρ

dt
+ρ∇·w=0

ρ
du

dt
=−
∂p

∂r
+
∂Srr
∂r
+
∂Srz
∂z
+
Srz −Sϕϕ
r

ρ
dv

dt
=−
∂p

∂z
+
∂Srz
∂r
+
∂Szz
∂z
+
Srz
r

ρ
de

dt
=−p∇·w+Srr

∂u

∂r
+Sϕϕ

u

r
+Szz

∂v

∂z
+Srz

(∂u

∂z
+
∂v

∂r

)

(1.1)

dSrr
dt
=2µ

(∂u

∂r
−
1

3
∇·w

)

+
(∂u

∂z
−
∂v

∂r

)

Srz

dSϕϕ
dt
=2µ

(u

r
−
1

3
∇·w

)

dSzz
dt
=2µ

(∂v

∂z
−
1

3
∇·w

)

−
(∂u

∂z
−
∂v

∂r

)

Srz

dSrz
dt
=µ
(∂u

∂z
+
∂v

∂r

)

−
1

2

(∂u

∂z
−
∂v

∂r

)

(Srr−Szz)

The von Mises limit of elasticity is assumed in the form

SijSij ¬
2

3
Y 2 (1.2)



Computer simulation of catching of ... 127

where

SijSij =S
2
rr+S

2
zz+S

2
ϕϕ+2S

2
rz (1.3)

The equation of state for metals is accepted in the form

p= k1x+k2x
2+k3x

3+γ0ρ0e
(1.4)

x=1−
ρ0
ρS

k2 =0 if x< 0

The temperature of the metal can be calculated from the relation

T =300
e0−e
e00

e0 = e00+e01x+e02x
2+e03x

3+e04x
4 (1.5)

where: r,z are space coordinates (axial symmetry), t– time, ρ–density, u,v –
mass velocity components along the r,z coordinates, p –pressure, e – internal
energy, T – temperature, ρs – density of the solid phase, Sik – components
of the stress deviator, Y – yield strength, µ – shear modulus.

The equation of state in the form of (1.4)1 is very convenient from the
practical point of view because on one hand it is valid in a wide range of the
pressure and temperature and on the other hand one can find in literature,
see Barbee et al. (1972), exact values of the constant coefficients: k1, k2, k3,
e00, e01, e02, e03, e04 for several most frequently used metals in the related
research.

For the description of strength properties of the metals, a modifiedmodel
using elements of the Steinberg-Guinan and Johnson-Cook models (Johnson
and Cook, 1983; Johnson and Lindholm, 1983; Steinberg, 1991; Steinberg et
al., 1980; Steinberg and Lund, 1989) is used

Y = [A+B(εp)n](1+C ln ε̇p
∗
)(1+ bp−Tm

∗
)F(ρS)

[A+B(εp)n]¬Ymax Y =0 if T >Tm0
µ=µ0(1+ bp−Tm∗ )F(ρS) (1.6)

εp =

√
2

3

√

(ε
p
rr −εpzz)2+(εprr−εpϕϕ)2+(εpzz −εpϕϕ)2+

3

2
(ε
p
rz)2

F(ρS)=



















1 if ρS ­ ρS1
ρS −ρS2
ρS1−ρS2

if ρS2 ¬ ρS <ρS1
0 if ρS <ρS2



128 K.Jach et al.

where ε
p
ik
denotes the components of tensor of plastic strain, εp is the in-

tensity of plastic strain, ε̇
p
∗ = ε̇

p/ε̇
p
0 – plastic strain rate for ε̇

p
0 = 1.0s

−1,
T∗ = (T −T0)/(Tm0 −T0), T0 and Tm0 – initial temperature and the melting
point temperature, A,B,C, n,m – material constants.
The system of equations describing the dynamics of the volume increase of

microcracks (microvoids) is assumedas inmodifiedFortow’smodel (Agurejkin
et al., 1984; Barbee et al., 1972; Johnson, 1981; Sugak et al., 1983)

dVc
dt
=















−ksgnp
[

|p|−σ0
Vc1
Vc+Vc1

]

(Vc+Vc0) if |p| ­σ0
Vc1
Vc+Vc1

0 if |p|<σ0
Vc1
Vc+Vc1

(1.7)

1

ρ
=Vc+

1

ρs
σ0 =σ00F(ρS)H(ε

p)(1−Tm
∗
)

The limitation of the strength properties brought about by appearing micro-
cracks is modeled bymultiplying Y , µ and η by the suitable function G(Vc)

Y⊤ =YG(Vc) µ
⊤ =µG(Vc)

(1.8)

(k1,k2,k3)
⊤ =(k1,k2,k3)G(Vc)

The functions G(Vc) and H(ε
p) are assumed in the form

G(Vc)=1−ρVc H(εp)= exp
(

−1
2
εp
)

(1.9)

2. Equations describing explosive detonation

In the description of processes of an explosive detonation, classical equ-
ations of hydrodynamics have been used

dρ

dt
+ρ∇·w=0 ρ

du

dt
=−
∂p

∂r

ρ
dv

dt
=−∂p
∂z

ρ
de

dt
=−p∇·w

(2.1)

System of equations (2.1) was then completed by the equation of state of
detonation products. It was assumed in the form of JWL

pPD =A
(

1− δ
R1V

)−R1V

+B
(

1− δ
R2V

)−R2V

+ δρe (2.2)



Computer simulation of catching of ... 129

where V = ρ0/ρ; A,B,R1,R2, δ denote empirical constants.

The systemof equationsmentioned abovewas applied in all cases, inwhich
so-called ”detonative optics” approximationwas used, i.e. when the knowledge
of the detonation wave front shape and of parameters on its front (Chapman-
Jouguet’s parameters) was assumed.

3. Behaviour of protective casing under pressure of detonation

products and shock wave in the air

Todescribe the process of protective casing deformation one can use a sys-
temof equations for an elastic solid in the spherical symmetry, see Stanyukovic
(1975)

dρ

dt
+ρ∇·w=0

ρ
du

dt
=−
∂p

∂r
+
∂S1
∂r
+2
S1−S2
r

ρ
de

dt
=−p∇·w+S1

∂u

∂r
+2S2

u

r
(3.1)

dS1
dt
=2µ

(∂u

∂r
− 1
3
∇·w

)

dS2
dt
=2µ

(u

r
− 1
3
∇·w

)

where

∇·w= ∂u
∂r
+2
u

r
(3.2)

The system of equations of the problem was then completed by the equ-
ation of state in the form of

p=K
(

1−
ρ0
ρ

)

(3.3)

The stresses were calculated as follows

σ1 =−p+S1 σ2 =−p+S2 (3.4)

where: σ1 stands for the radial stress, σ2 – longitudinal stress, K – bulk
modulus, µ – shear modulus.



130 K.Jach et al.

4. Model of deceleration of steel fragments in the air

The equation describing motion of metal fragments in the air has been
assumed in the following form

m
dv

dt
=−
1

2
cρav

2s m=
4

3
π
(D

2

)3

ρs (4.1)

where: D is thediameter of the steel ball, s– surface of thebody (s=πD2/4),
v – velocity of the ball, ρa – density of the air, ρs – density of steel, m –mass
of the steel ball.

The coefficient cwas approximated as follows (Cerny̌ı, 1988)

c=−0.2
(Re

105

)

+0.7 (4.2)

where Re is the Reynold’s number

Re=
ρpvD

η

5. Constant coefficients of materials used in calculations

The coefficients determining material properties of the examined bodies
are presented in the Table 1 and Table 2.

Table 1.Coefficients of the equation of state JWL for trinitrotoluene

ρ0 [g/cm
3] 1.63 A [GPa] 373.8 R2 0.9

D [m/s] 6930 B [GPa] 3.747 δ 0.35

pCJ [GPa] 21.0 R1 4.15 ρ0e0 [GPa] 5.9

For the casing one has K = 8.8GPa, µ = 300GPa, Y = 0.1GPa,
pmin = −0.3GPa, ρcr = 0.5ρ0 (pmin and ρcr represent the critical pressu-
re and density values at tearing the quasi-aramide fabric). The air viscosity
η=3 ·10−4g/(cm·s) has been assumed in the calculation.
To simulate the penetration of the quasi-aramide casing by a propelled

metal fragment, equations (1.1)-(1.3) of elastoplastic body were used. In this
case, for the quasi-aramide fabric, a simplified strength model was applied:
Y = const1 and µ= const2, in the plasticity area.



Computer simulation of catching of ... 131

Table 2. Values of coefficients occurring in the equation of state, in the
model of microcracks formation and in the Johnson-Cookmodel for steel

ρ0 [g/cm
3] 7.9 C 0

γ0 2.17 m 0.55

k1 [GPa] 164.8 n 0.32

k2 [GPa] 312.4 Ymax [GPa] 2.0

k3 [GPa] 564.9 Tm0 [K] 1811

ε00 [J/g] −1.340 ·102 k [1/(Pa·s)] 0.25
ε01 [J/g] −2.908 ·102 σ0 [GPa] 2.5
ε02 [J/g] 1.012 ·104 µ0 [GPa] 77
ε03 [J/g] 2.051 ·104 VC0 [cm3/g] 1.27 ·10−5
ε04 [J/g] 2.901 ·104 VC1 [cm3/g] 6.33 ·10−4
A [GPa] 344 ρS1 [g/cm

3] 6.87

B [GPa] 680 ρS2 [g/cm
3] 5.84

6. Analysis of calculations results

Figure 1 presents a sequence of snapshots of the process of propelling of
the standard steel fragment of about 1.1g, located on a 90mm long, 200g
trinitrotoluene block.

Fig. 1. Snapshot sequence of the standard steel fragment propelling process. The
results were obtained using the free particle method (Jach et al., 2001). Initial
distance of trinitrotoluene block (1) frommetal fragment (2) is 90mm



132 K.Jach et al.

The fragment, propelled explosively, reaches the velocity of about 800m/s
(Fig.2). A change in the velocity caused by deceleration in the air, and the
increasing rangeof the fragment in time fordifferent initial velocities are shown
in Fig.3. For ranges of the order of 0.5-1.0m, the deceleration due to the air
is negligible.

Fig. 2. Velocity V and range L of the fragment propelled by products of detonation
(the initiation as in Fig.1) as a function of time t (since the initiation of detonation)

Fig. 3. Changes of velocity V and range L of fragments during motion in the air for
different initial velocities V0

As a result of detonation of the trinitrotoluene block, the shape of the
propelled fragment becomes a little deformed plastically in the first phase of
the propelling, when the pressure of detonation products exceeds the yield
strength (Fig.4).

Figure 5 shows additionally changes of the energy density of the fragment
for different initial velocities (as a result of deceleration in the air). In the
distance of several metres from the explosion, the energy density still signi-
ficantly exceeds the value which is assumed as dangerous for human life, i.e.
100-150J/cm2.

The following figures present basic characteristics concerning the process
of shock wave propagation in the air and of its influence on the protective



Computer simulation of catching of ... 133

Fig. 4. The shape of the propelled fragment. Results of the computer simulations
based on the free particles method

Fig. 5. Energy density E
w
of fragments as a function of their range L in the air

casing. The results have been obtained using the Lagrange method. Changes
of the pressure on the shock wave front in the air are shown in Fig.6. This
wave reaches the surface of the internal protective casing and, as a result of
reflection from it, significantly increases its own amplitude. The change of the
pressure on the edge of the protective casing in function of time is presented
in Fig.7. The second, moderate local maximum of the pressure results from
reaching the edge by the secondary shock wave reflected from the border:
products of detonation – the air.



134 K.Jach et al.

Fig. 6. The pressure at the shock wave front in the air as a function of its range

Fig. 7. The pressure on the internal surface of the protective casing as a function of
time

The maximum tensile stresses, stretching the casing as a result of the
pressure action (Fig.7), occur on the internal edge of the protective layer. The
changes of the tensile stresses in time are presented in Fig.8.

Fig. 8. Longitudinal tensile stress on the internal surface of the protective casing

Another important issue is also the problem of the time required by the
air shock wave and fragments to reach the protective casing. Changes of the
range of the shock wave and of the propagating fragment (with the velocity
of 800m/s) are illustrated in Fig.9. From this figure it results that the shock
wave in the air overtakes the fragment significantly. In the time of about



Computer simulation of catching of ... 135

160µs the shockwave reaches the edge of the protective layer (with the radius
of 400mm), and in this time the fragment travels a distance of only about
120mm.

Fig. 9. Comparison of the range of the shock wave in the air with the range of the
steel fragment

Themultilayered fabric is impacted bya randomside of the propelled frag-
ment. In order to evaluate the catching properties of the casing, the standard
cylindrical fragment was used in the calculations. Figure 10 shows results of
computer simulation, based on the free particles method, of penetration of
the protective fabric by a steel fragment propelled to the velocity of 900m/s
(the last picture refers to themomentwhen the fragment has been completely
stopped). As one can see, during the penetration of the protective casing, the
fragment is hardly deformed.

Fig. 10. Computer simulation of penetration of the protective multilayered fabric by
the steel fragment



136 K.Jach et al.

Figure 11 shows the decrease of the steel fragment velocity as a result of
its interaction with the protective fabric for three different values of velocities
of the impact of the fragment.

Fig. 11. The decrease of the velocity of the steel fragment as a result of deceleration
in the protective fabric for three different values of velocities of the impact of the
fragment – results of calculations obtained using the free particles method

7. Conclusions

The following remarks can be formulated based on the results of computer
simulations carried out throughout the work:

• The time that the shockwave needs to reach the inner surface of the pro-
tective casing is significantly shorter than the analogous time pertaining
to the fragment.

• The change of the fragment shape, i.e. mushrooming of the surface of
the fragment is caused by the pressure impact of the explosion products.
It does not occur during the penetration through the protective casing.

• It is possible to select appropriate input data in the calculations so that
the simulation results would agree with the results of the experiment
(e.g. the depth of the protective casing penetration ocurred to be in line
with the experimental results).

• The computer codes used can be applied to optimisation of parameters
of the protective casing, i.e. its thickness, sensitivity to the kind, shape
and mass of the explosive as well as the kind, shape and mass of the
propagating fragments.



Computer simulation of catching of ... 137

References

1. Agurěıkin V.A., Anisimov S.I., Bushman A.V., Kanel G.I., Kary-
aginV.P.,KonstantinovA.B.,KryukovB.P.,MininV.F.,Razorenov
S.V., Sagdeev R.Z., Sugak S.G., FortovV.E., 1984,Thermophisical and
gas dynamics problems of protection againstmeteors of the Vega space vehicle
[in Russian],Teplofizika Vysokikh Temperatur, 22, 5

2. Barbee T.W. Jr., Seaman L., CrewdsonR., CurranD.R., 1972,Dyna-
mic fracture criteria for ductile and brittle metals, J. Mater. Sci., 7, 393-401

3. Cierny̌ı G.G., 1988,Gas dynamics [in Russian], Science, Moscow

4. Jach K., Morka A., Mroczkowski M., Panowicz R., Sarzyński A.,
Stępniewski W., Świerczyński R., Tyl J., 2001, Computer modelling of
dynamic interaction of bodies by free particles method [in Polish], PWN

5. Jach K., Włodarczyk E., 1992, Solutions of the initial-value problems of
the viscoplastic – nonstationary theory for the description shaped charge jet
formation and targetpenetration,Proceedings of 13th International Symposium
on Ballistics, Ballisics’92, Stockholm, Sweden

6. JohnsonJ.N., 1981,Dynamic fractureand spallation inductile solids,J.Appl.
Phys., 52, 4, 2812-2825

7. JohnsonG.R., Cook W.M., 1983,A constitutivemodel and data formetals
subjected to large strains and high temperatures, in: 7th Int. Symposium on
Ballistics, Hague

8. Johnson G.R., Lindholm U.S., 1983, Strain-rate effects in metals at large
shear strains, Material behavior under high stress and ultrahigh loading rates,
Proc. 29th Sagamore Army Mater. Res. Conf., Lake Placid 1982, NewYork

9. Stanyukovic K., 1975,Physics of explosion [in Russian], Science, Moscow

10. Steinberg D.J., 1991, Equation of state and strength properties of selected
materials, Lawrence LivermoreNat. Lab. UCRL-MA-106439

11. SteinbergD.J., Cochran S.G.,GuinanM.W., 1980,A constitutivemodel
for metals applicable at high-strain rate, J. Appl. Phys., 51, 3, 1498-1504

12. SteinbergD.J., LundC.M., 1989,A constitutivemodel for strain rates from
104 to 106s−1, J. Appl. Phys., 65, 4, 1528-1533

13. Sugak S.G., Kanel G.I., Fortov V.E., Ni A.L., Stelmah B.G., 1983,
Numericalmodelling of the action of an explosion on an iron plate [inRussian],
FGV, 19, 2

14. Wilkins M.L., 1984,Modelling the behaviour of materials. Structural impact
and crashworthiness,Proc. Intern. Conf., 2, London, NewYork



138 K.Jach et al.

Komputerowa symulacja wychwytywania przez warstwę ochronną

metalowych odłamków napędzonych wybuchem

Streszczenie

Wpracyprzedstawionomatematyczno-fizycznymodel orazwyniki symulacji kom-
puterowych napędzania metalowych odłamków przez produkty detonacji i wychwy-
tywania ich przez warstwę ochronną wykonaną z tkaniny paradramidowej. Przeana-
lizowano również wpływ generowanej wybuchem fali uderzeniowej w powietrzu na
przebieg zachodzących procesów. W zależności od charakteru tych cząstkowych za-
gadnieńwykorzystywanodo symulacji komputerowej kody typu 2D (metoda punktów
swobodnych) lub 1D (metoda Lagrange’a).

Manuscript received July 10, 2003; accepted for print September 25, 2003