7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Compressibility and Unsteady Effects on Particles in Explosions
Y. Ling* A. Haselbacher* and S. Balachandar *
Department of Mechanical and Aerospace Engineering, University of Florida, Gainesville, FL 32611, USA
yueling@ufl.edu, balals@ufl.edu and haselbac@ufl.edu
Keywords: Particle dispersal, Explosion, LagrangianEulerian Simulation, Compressibility, Unsteady effect
Abstract
An unsteady compressible multiphase flow model is proposed in the present study to solve the extended problem of
Brode (Brode, H. L., "Numerical Solutions of Spherical Blast Waves," J. Appl. Phys., Vol. 26, 1955, pp. 766775),
namely an explosion generated by a spherical charge containing a compressed gasparticle mixture. Oneway coupled
simulations are carried out for different particle properties. The main goal of the present study is to investigate the
importance of compressibility and unsteady effects on particles' behavior in explosion.
Nomenclature
Roman symbols
ReP Particle Reynolds number
Mp Particle Mach number
NuP Particle Nusselt number
Pr Gas Prandtl number
Subscripts
i Lagrangian particle index
qs Quasisteady component
iu Inviscidunsteady component
Superscripts
g Gas
p Particle
Introduction
Different from the conventional single phase (gas) ex
plosions, multiphase explosions contain solid particles
as an additional phase. The problem of multiphase
explosions with solid particles has drawn more atten
tion recently, due to some important applications, such
as EBX(Enhanced Blast Explosive), (see Zhang et al.
(2001); Tanguay et al. (2007)), and enhanced mixing
and afterburn of detonation products (Balakrishnan &
Menon (2010)). Generally, the material of these solid
particles is reactive, such as Aluminum; and the motiva
tion of adding particles in explosives is to increase the
total energy release. The design of such explosives re
quires a good understanding of particle behavior in ex
plosions. Accurate prediction of particle motion is cru
cial, because the burning of a particle is directly related
to its location in the flow. For example, the region be
tween the main shock and the gas contact is good for
particle burning, because the gas there is heated by the
main shock, and is rich of oxygen. On the other extreme,
if a particle locates in the region inside the gas contact,
its surrounding gas is expanded detonation products, and
the reactions of particles may be completely suppressed.
However, accurate prediction of particle behavior in
explosions is a very difficult task. The unsteady com
pressible nature of the explosion flow coupled with the
motion of dispersing particles makes this problem chal
lenging. In order to build a better understanding of the
physics of multiphase explosion, we employ a classic
explosion model problem considered by Brode (1955),
and extend it to multiphase case. In his study, Brode sim
ply modeled the explosion process by a sudden release
of a spherical charge of highly pressurized gas. This
problem has since been intensively studied in the single
phase flow regime, see Brode (1955), Boyer (1960), and
Friedman (1961). It can also be treated as an extension
of the classic point source explosion theory developed
by Taylor (1950) and Sedov (1959). The added compli
cation considered in the present study is to replace the
pressurized gas inside the spherical charge by the pres
surized gasparticle mixture.
The problem of explosion of a spherical charge con
taining solid inert particles is far more complicated than
the case without particles. Unlike planar geometry, in
the case of a spherical multiphase explosion, often the
particle contact (the front of the particle cloud) can move
ahead of the gas contact, thus exposing particles to am
bient air. In fact, it has been observed that under some
initial conditions, the particle contact can even penetrate
ahead of the main shock, see Lanovets et al. (1993);
Zhang et al. (2001). The complicated interactions be
tween particles and different waves in compressible flow
need to be captured to guarantee accurate prediction on
the particle dispersal process.
Based on the previous works in both the single phase
and multiphase regimes, in the present study, we carry
out a systematic study of the extended Brode's problem
through LagrangianEulerian (LE) simulation. The gas
is treated as a continuum and computed using the Eu
lerian framework; while the particles are treated as a
discrete phase and described under Lagrangian frame
work. The strong compressibility and unsteady effects
associated with the rapidly expanding flow introduce im
portant difficulties in multiphase flow modeling. The
commonly used standard drag law and heat transfer law,
such as the correlation by Clift & Gauvin (1970) and
Whitaker (1972), ignore both the compressibility and
unsteady effects, and thus lead to large errors in the cur
rent problem. The recent advancements on particulate
force and heat transfer by Parmar et al. (2009, 2010)
and Balachandar & Ha (2001) are employed here to
develop a multiphase model for unsteady compressible
flow. Simulations were performed in the limit of one
way coupling in the spherical geometry. A wide range
of particle density and size are considered to systemati
cally investigate the compressibility and unsteady effects
on different particles in explosion flows.
Mathematical Model and Numerical Method
The physical model of multiphase flow presented here is
based on the following assumptions: (1) The fluid may
be represented by a perfect gas. (2) The fluid motion
may be regarded as inviscid, so its viscosity and con
ductivity are neglected except in the interaction with the
particles. (3) The particles are inert, rigid, and spherical.
(4) The particles have constant heat capacity and uni
form temperature distribution. (5) The volume fraction
of the particles is negligible. (6) The particles do not
collide with each other. (7) The multiphase flow is con
sidered to be oneway coupled, namely particles have no
influence on gas.
The EulerianLagrangian approach is employed to
model the multiphase flow. The gas phase is governed
by the Euler equations with source terms. Based on the
assumption of spherical symmetry, the mass, momentum
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
and energy equations of gas can be written as:
9(p9) d(p9u9) 2p9u9
at Or r
d(p9gv) d[p(u)2 +pg] 2p (U )2
d(pgEg) d(p9H9u9) 2p9H9u9
at i3r r
0, (1)
0, (2)
0, (3)
where p, u,p, E, and H represent the density, velocity,
pressure, total energy, and total enthalpy. The super
scripts g and p represent the gas and particle phase, re
spectively. The subscript i denotes quantities pertaining
to the ith Lagrangian particle (for example .., denotes
the velocity of the ith Lagrangian particle).
Assuming that the gas follows the perfect gas law and
the particles not to contribute to the pressure, the equa
tion of state can be written as
1' = pSRST.
The corresponding position, momentum, and energy
equations of the particles in the Lagrangian framework
can be written as:
where x., ..p and Tf represent the position, velocity, and
temperature of the ith particle.
In the present study, both the particle force and heat
addition are calculated by the sum of quasisteady and
inviscidunsteady components. The force exerted on a
particle per unit mass is given as:
fi = fqs,i + f,i i (8)
where the subscripts qs and iu denote quasisteady and
inviscidunsteady components, respectively. The quasi
steady force ff,, used here can be expressed as
S CDiReP
fqs,i TP 24 (9)
The quantities with superscript g and subscript i denote
gas quantities at the ith Lagrangian particle location (for
example 'represents the gas velocity at the ith particle
location). ,Ti represents the mechanical response time of
the particle, which is defined as:
tP 1(di)2
Ti ,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
where ;, ', p, and df represent the gas viscosity, parti
cle density and diameter, respectively. ;, is computed
according to the Sutherland's law. CD,i in Eq. (9) is the
quasisteady drag coefficient, which is a function of par
ticle Reynolds number Ref and Mach number I. The
expression for CD,i CD,i(ReP, .L ) is given by the
empirical correlation in Parmar et al. (2010).
The definitions of Ref and I are:
Rep 1(11)
,,;
where al = vl^RT[ is the speed of sound of the gas.
Different from the integral form suggested by Parmar
et al. (2009), the inviscidunsteady force ff, used here
is in the incompressible limit I 0) for simplicity.
ftu, Cam,ef p( Dt D dt Dt
(13)
where C,,,.ef is the effective addedmass coefficient,
which is a function of I Here we simply used
C,,a,ef 1/2, namely the value of effective added
mass coefficient at I 0.
Similarly, the heat transferred to a particle per unit
mass can also be expressed as the sum of quasisteady
and inviscidunsteady components:
qi qP ', + qP, (14)
The quasisteady heat addition to a particle can be ex
pressed as:
T9 T ''11i
T^= 2 (15)
O,i
where 7P is the thermal response time of the particle,
which is defined as:
p fpfP(di)2
P (16)
',i 12_', (16)
where Cf and are the specific heat of the particle and
the gas thermal conductivity. '.1I, is the particle Nus
selt number given by the heattransfer law of Whitaker
(1972):
': 2 + [0.4(E. + 0.06(ReP)2/3]Pr04 (17)
where Pr is the Prandtl number of gas,
/K9
and Cfl represents the specific heat at constant pressure
of the gas. In the present study, Pr is considered to be
constant.
The inviscidunsteady heat transfer in the unsteady
compressible flow is not found in the literature. How
ever, the expression in the incompressible limit given by
Balachandar & Ha (2001) can be employed as an ap
proximation, which is consistent with the treatment for
the inviscidunsteady force.
(DT'\
... _1 DT) (19)
The numerical approach for the gas equations is based
on the cellcentered finitevolume method. The inviscid
flux for the gas is calculated by the approximate Rie
mann solver of Roe (1981). The facestates are obtained
by a simplified secondorder accurate weighted essen
tially nonoscillatory (WENO) scheme (see, e.g., Jiang
& Shu (1996)), which modifies the gradients computed
using the leastsquares reconstruction method of Barth
(1991). The temporal integration method for both gas
and particle equations is the fourthorder RungeKutta
method. The solution method has been extensively ver
ified and validated for gas and gasparticle flows with
and without shock waves. For brevity, the results from
these studies are not reproduced here. For more de
tails, see Haselbacher (2005); Haselbacher et al. (2007).
The gas variables appearing in the particle equations are
evaluated through Hermite interpolation from gas solu
tions stored in neighboring cell centers, see Ling et al.
(2010a). The material derivatives of gas velocity and
temperature in Eqs. (13) and (19) are calculated by the
method in Ling et al. (2010b).
Results
In the present problem, the explosion is generated from
a spherical charge containing pressurized gasparticle
mixture. The particles are initially uniformly distributed
inside the charge, which are also considered to be in
equilibrium with the surrounding gas. After the bound
ary of the spherical charge is removed, a shock wave (the
socalled main shock) propagates outward, followed by
the gas contact separating the gases originally contained
inside and outside the charge. Because of the finite in
ertia of the particles, the particle contact will initially
lag behind the gas contact. Simultaneously, an expan
sion fan propagates inward. In contrast to the planar
shocktube problem, a second shock is generated due
to the radial effect. This second shock first travels out
ward, but eventually reverses direction. Figure 1 shows
a schematic describing the problem, where the numbers
4, 3, 2, 1, and 0 represent the regions inside of the ex
pansion fan, between the expansion fan and the second
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Compressed
gasparticle mixture
Air
Contact
Shock
s Gas contact
Second shock
Head of Expansion Fan
Figure 1: Schematic illustration of explosive dispersal
from a spherical charge containing a gasparticle mix
ture.
shock, between the second shock and the gas contact,
between the gas contact and the main shock, and out
side of the main shock, respectively. The shaded area
indicates the presence of particles.
In the present problem, for simplicity the gas is taken
to be air (y9 1.4, R = 287 J/(kg K), Pr 0.72).
As we only consider oneway coupled simulations, the
gas solution is fully characterized by the initial condi
tions inside and outside of the charge i' /,, / ; and
,,,), and the radius of the spherical charge Rse. In
Brode's study, he used 14/1.,, 121, p /p = 121,
and Rs 0.0254 m, which are maintained the same in
the present study. If Rs and the ambient speed of sound
a are taken to be the typical dimensional length and ve
locity scales, the governing equations can be rewritten in
dimensionless form. The particle phase can be charac
terized by the particle diameter dp, density pP, specific
heat CP. In order to investigate the compressibility and
unsteady effects on particles of different size and den
sity, multiple computations were carried out for dP/Rsc
varying from 4 x 104 to 4 x 10 2, and pP/, varying
from 1 to 100. As the thermal behavior is not the fo
cus of the present paper, CP is fixed as 1600 J/(kg K).
All the cases of different particles were computed using
both the present model and also the conventional model
(namely the standard drag law (Clift & Gauvin (1970))
and heat transfer law (Whitaker (1972)). So that, we can
compare the results of the two models, and the improve
ment by including compressibility and unsteady effects
in the present model can be evaluated.
Figure 2 shows the r t diagrams of the solutions
for different pP and dp. Generally, a particle initially lo
cated closer to the charge boundary can be accelerated
by the stronger initial period of the expansion fan. The
width of expansion fan increases when it propagates in
ward. The particle exposed to stronger expansion fan
will gain larger momentum, and propelled to larger dis
tance. During the interaction with the strong expansion
fan, the acceleration of the gas around the particles is
very large, therefore, the inviscidunsteady force is ex
S'   Main shock
5 6  Gas contact
2nd shock
5 t / Gas contact
t/ (present)
4e4 (present)
=4e6(CG)
=4e5(CG)
S  i ="4e4(CG)
6 5 10 15 20
r/Rsc
(b) d /R 4 x 103
10
Figure 2: rt diagram of main shock, gas contact, second
shock wave, and particle contacts for different particle
density and diameter. "CG" and 'ilc'icii denote results
computed by the drag law of lift & dPauvin (1970) and
entthe present force model.
pected to be large, too. Figure 2(a) shows the trajec
tories of the particle contacts for different particle den
sity, when diameter is fixed at dp/R8s = 4 x 10o3. It
can be observed that, for large particle density, such as
p/pf4 > 10, the difference between results computed by
the present force model and those by the standard drag
increase it rtie araadensity. Similarly, Fig. 2(b)t shows
the trajectories of particle contact with different diame
density and diameter. "CG" and Ili d i/R denote results
ter but with the same density law of Clif, & Gauvin 10). The devi
ation be observed that, for large particle size is large, such
as / > 4 x The difference between results computed bysent
the present force model and those by the standard drag law increases with particle di
law is pretty large. The deviation of results computed
by standard drag law is due to the neglect of compress
ibility and unsteady effects on particles, which seems to
increase with particle density. Similarly, Fig. 2(b) shows
the trajectories of particle contact with different diame
ter but with the same density (, .'/, . 10). The devi
ation is also ,ignili,,.mi when particle size is large, such
as dP/Rsc > 4 x 10 The difference between present
model and standard drag law increases with particle di
ameter.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
12
present
O 8 Chlft & Gauvin
06
04
02
0 20 40 60 80 100 120
t/(Rsclag)
Figure 3: Evolution of particle velocity for the leading
particle (particle contact) of dP/R, 4 x 103 and
pP/p = 10.
The impact of the compressibility and unsteady ef
fects on particles can be explained more clearly through
the time evolution of particle velocity and force break
down, which are shown in Figs. 3 and 4. It can be ob
served in Fig. 3 that the particle velocity by the present
model increases faster and reaches a larger maximum
value than by the standard drag law. This rapid accelera
tion happens when the particles interact with the expan
sion fan. The results show that if the compressibility and
unsteady effects are ignored, the particle acceleration
can be underestimated. After the particles exit the ex
pansion fan (t/(Rs/a ) > 3, particles start to deceler
ate), then move in a region where the acceleration of the
surrounding gas is not large, then it can be observed that,
the particle velocity evolutions between the two force
models tend to merge together. In this time period, the
compressibility and unsteady effects are no important,
standard drag law is sufficient to describe the particle
motion. Figure 4 shows the quasisteady and inviscid
unsteady components of particle force, from which it
can be seen that the inviscidunsteady force is absolutely
dominant at the early time (t/(Rs/a ) < 3), when the
particles accelerate inside the expansion fan. After that,
the quasisteady force becomes the dominant force, ex
cept when the particles interact with the second shock
wave and the gas contact. It can be observed that the
inviscidunsteady force becomes dominant again for a
very short time when particle cross the second shock.
However, the influence of inviscidunsteady force on the
particle motion seems to be not that important since its
duration is very short.
The compressibility and unsteady effects also have
similar influence on particle thermal behavior. Figure 5
shows the temperature evolution of the leading particle.
Similar trend as the particle velocity can be observed.
1 5E+07
1E+07
z
S5E+06
LL
102 10 100 101 1C
t/(Rj/ag)
Figure 4: Evolution of particle force per unit mass for
the leading particle (particle contact) of dP/R, 4 x
103 and pP/, 10.
present
Whitaker
098)
096
So
E
094
092
"0 20 40 60 80
t/(Rslao)
100 120
Figure 5: Evolution of particle temperature for the lead
ing particle (particle contact) of dP /Rs 4 x 103 and
pP/p= 10.
The particle temperature decreases faster and reaches to
a lower minimum value with the present model than with
the standard heattransfer law (Whitaker (1972)). As the
gas temperature drops very fast inside the strong expan
sion fan (large negative DT9/Dt), the inviscidunsteady
component of heat transfer is dominant. After particle
exits the expansion fan, the quasisteady component be
comes dominant.
Finally, Fig. 6 shows the deviation between particle
contact locations computed by the present model and the
standard drag law, (xzresent xzta)/R at a very long
time after the explosion starts (t/(Rs,/a ) 136) as a
function of particle density and diameter. This figure can
be viewed as a summary of the significance of the com
pressibility and unsteady effects on particle motion over
a wide range of particle properties. Several important
35 ,
4 .
1 05 0 05 1
L :,.,ii,i' 1 ,
Figure 6: Deviation of particle contact location after
a long time (t/(Rsc/ac) 136) computed by stan
dard drag law from the results computed by the present
model, as a function of pP and dp.
observations can be made here. First, it can be observed
that the deviation generally increases with diameter and
density, though when particle diameter is large, the vari
ation of the deviation may be nonmonotonic with parti
cle density. Second, for the range of particles considered
here, the error of particle contact location by ignoring
compressibility and unsteady effects can reach 30 times
of the charge radius. Third, when both particle density
and diameter are very small (left bottom comer of the
contour), the deviation between the two models is neg
ligibly small (less than 0.01 charge radius). In this case,
the standard drag law is sufficient.
Conclusions
The classic explosion problem of Brode (1955) is ex
tended to multiphase case in the present study to in
vestigate the compressibility and unsteady effects on
particles in explosion flows. In this problem, particles
are accelerated by the strong expansion fan generated
when the pressurized gasparticle mixture is suddenly
released. When particle interacts with the strong ex
pansion fan, the inviscidunsteady components of both
particle force and heat transfer are important for a wide
range of particle diameter and density, as Dug /D and
DT9/Dt are both very large inside the expansion fan.
Ignoring the compressibility and unsteady effects will
lead to i iniik.iia underestimate on particle acceleration
and temperature decrease. Though the inviscidunsteady
force and heat transfer are active only in a short time pe
riod, they are important to the particle behavior. The rea
son is the inviscidunsteady components have strong im
pact on how much momentum and energy the particles
gain from the explosion. After acceleration inside the
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
expansion fan, the particles decelerate and are mainly
governed by the quasisteady force and heat transfer.
Acknowledgments
The authors gratefully acknowledge support by the
National Science Foundation under grant number
EAR0609712.
References
Balachandar S. and Ha M. Y., Unsteady heat transfer
from a sphere in a uniform crossflow, Phys. Fluids, Vol.
13, pp. 37143728, 2001
Balakrishnan K. and Menon S., On the Role of Ambi
ent Reactive Particles in the Mixing and Afterburn be
hind Explosive Blast Waves, Combust. Sci. Technol.,
Vol. 182, pp. 186214, 2010
Barth T.J., A 3D Upwind Euler Solver for Unstructured
Meshes, AIAA Paper 911548, 1991
Boyer D.W., An Experimental Study of the Explosion
Generated by a Pressurized Sphere, J. Fluid Mech., Vol.
9, pp. 401429, 1960
Brode H.L., Numerical Solutions of Spherical Blast
Waves, J. Appl. Phys., Vol. 26, pp. 766775, 1955
Clift R. and Gauvin W. H., The motion of particles in
turbulent gas streams, Proc. Chemeca '70, Vol. 1, pp.
1428, 1970
Friedman M.P, A Simplified Analysis of Spherical and
Cylindrical Blast Waves, J. Fluid Mech., Vol. 11, pp. 1
15, 1961
Haselbacher A., A WENO Reconstruction Algorithm
for Unstructured Grids Based on Explicit Stencil Con
struction, AIAA Paper 20050879, 2005
Haselbacher A., Najjar F.M., Balachandar S. and Ling
Y., Lagrangian Simulations of ShockWave Diffraction
at a RightAngled Comer in a ParticleLaden Gas, 6th
International Conference on Multiphase Flow, Leipzig,
Germany. 2007
Jiang G.S. and Shu C.W., Efficient Implementation of
Weighted ENO Schemes, J. Comput. Phys., Vol. 126,
pp. 202228, 1996
Lanovets V.S., Levich V.A., Rogov N.K., Tunik Yu.V.
and Shamshev K.N., Dispersion of the detonation prod
ucts of a condensed explosive with solid inclusions,
Combustion, Explosion, and Shock Waves, Vol. 29, pp.
638641, 1993
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Ling Y., Haselbacher A. and Balachandar S., A numeri
cal source of smallscale numberdensity fluctuations in
EulerianLagrangian simulations of multiphase flows, J.
Comput. Phys., Vol. 229, pp. 18281851, 2010a
Ling Y., Haselbacher A. and Balachandar S., Numerical
investigation of particle dispersal in multiphase explo
sions, AIAA Paper 20100768, 2010b
Parmar M. and Haselbacher A. and Balachandar S., On
the unsteady inviscid force on cylinders and spheres in
subcritical compressible flow, Phil. Trans. R. Soc. A,
Vol. 366, pp. 21612175, 2009
Parmar M. and Haselbacher A. and Balachandar S., An
improved drag correlation for spheres and application to
shocktube experiments, AIAA J., under review, 2010
Roe PL., Approximate Riemann Solver, Parameter Vec
tors, and Difference Schemes, J. Comput. Phys., Vol. 43,
pp. 357372, 1981
Sedov L.I., Similarity and Dimensional Methods in Me
chanics, Academic Press, New York, 1959
Tanguay V., Higgins A.J. and Zhang F., A simple ana
lytical model for reactive particle ignition in explosives,
Propell. Explos. Pyrot., Vol. 32, pp. 371384, 2007
Taylor G.I., The Formation of a Blast Wave by a Very
Intense Explosion. I. Theoretical Discussion, Proc. Roy
Soc A, Vol. 201, pp. 159174, 1950
Whitaker S., Forced Convection Heat Transfer Correla
tions for Flow in Pipes, Past Flat Plates, Single Spheres,
and for Flow in Packed Beds and Tube Bundles, AIChE
J., Vol. 18, pp. 361371, 1972
Zhang F, Frost D.L., Thibault PA. and Murray S.B.,
Explosive dispersal of solid particles, Shock Waves, Vol.
10, pp. 431443, 2001
