Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (2011)
A. Bargteil and M. van de Panne (Editors)

Procedural Fluid Modeling of Explosion Phenomena
Based on Physical Properties
Genichi Kawada and Takashi Kanai
University of Tokyo, Graduate School of Arts and Sciences, Japan

Abstract
We propose a method to procedurally model the fluid flows of explosion phenomena by taking physical properties
into account. Explosion flows are always quite difficult to control, because they easily disturb each other andchange rapidly. With this method, the target flows are described by control paths, and the propagation flows
are controlled by following these paths. We consider the physical properties, which are the propagations of thepressure generated by the ignition, the detonation state caused by the pressure and the fuel combustions. Velocity,
density, temperature and pressure fields are generated procedurally, and the fluid flows are computed from thesefour fields based on grid-based fluid simulations. Using this method, we can achieve a fluid motion that closely
resembles one generated solely through simulation. This method realizes the modeling of flows controlled frameby frame and follows the flow’s physical properties.
Categories and Subject Descriptors (according to ACM CCS): I.3.7 [Computer Graphics]: Three-Dimensional
Graphics and Realism—Animation; I.6.8 [Simulation and Modeling]: Types of Simulation—Animation

1. Introduction
Explosion phenomena are defined as the propagation of
fluid, the flowing density of flame with compressible expansion behavior caused by the chemical process consisting oflow- and high-speed fluid aspects (see [Bak73] and [HS94]).
The high-speed aspect is where ignition occurs and the fluids propagate rapidly by expansion. The low-speed aspect
is where the fluids continue to propagate with a gradualloss of velocity, and then converge at the destination. Explosion flows show properties such as whirling and streaming motions throughout the entire process. Ignition causes
the propagation of pressure from the explosion source, andthis strong pressure is called a shock wave. This pressure
determines the physical explosion properties. We focus ona model of the explosion phase prior to the flow reachingits destination, while previously-proposed approaches havedealt with what occurs after the flow reaches its destination.
Explosion behaviors are extremely complex, in particular
high-speed flows and large turbulence have to be consideredwhen attempting to express such behaviors. Therefore, explosion flows are difficult to define using simple geometricfeatures or only physical parameter settings at the ignition.
Adding to that, a grid size for the simulation typically has tobe small enough to realize the whirling and streaming flows.
c The Eurographics Association 2011.
⃝

However, with the small grid size it is difficult to realize
high-speed flows, because the time steps need to be small
according to a CFL condition and that leads to huge computational costs.
In related work, Kim et al. proposed a method to guide
flows by generating velocity fields along a path to control
smoke as a general fluid [KMD06]. Simulations for explosions are, however, different from simulations for other fluid
phenomena, and it is difficult to lead the behaviors withlarge vortices and high-speed flows using only the methodof Kim et al.. Dobashi et al. [DSYA09] proposed a control method for explosion by directly adjusting the velocity fields in grid-based simulation. It is, however, quite difficult for this method to realize the detailed control for highspeed flows as a unique explosion property, or to generate
the highly curved flows with occlusions from the source. We
believe that explosions must be characterized step by step
based on the propagation behaviors of the flows during thesimulation. This is because the model should have enoughinformation about the flows during each step to describe thewhole propagation in a relatively short time before the flow
reaches its destination.
We propose a procedural method to model explosion phenomena by considering the physical properties. In addition

G. Kawada & T. Kanai / Procedural Fluid Modeling of Explosion Phenomena Based on Physical Properties

to controlling the explosion flows based on the paths specified by the user, the flow fields are procedurally generated
and the grid-based simulations are applied to the fields. In
this way, a smooth transition from high-speed flows to lowspeed flows can be realized. We generate density, velocity,temperature and pressure fields based on explosion curves
derived from the physical explosion models, and realize thecomplex explosion behaviors including large vortices andpressure change in flows. Our method also takes into consideration the detonation phenomena to cause the drasticchange of pressure and the fuel combustion to cause large
vortices.

Figure 1: Explosion curves. Left: Pressure magnitude curve.
Right: Density opacity curve, pressure propagation curveand density propagation curve.

We believe our method has the following advantages:• While our mathematical equations are directly and theoretically derived from the physics, we first present practical ways to apply those equations to compressible andexplosive fluids for computer graphics simulation.
• Although the flow behaviors are less physically accurate
than full grid-based simulation, each flow following acontrol path exhibits the appropriate physical properties.
• The user inputs are the specifications of non-uniform rational B-spline (NURBS) curves and the parameters of the
explosion curve. The explosion flows are automaticallycalculated once the inputs are complete. Therefore, ourmethod realizes modeling of the flows that is free, easy to
predict and intuitive for the user.
• In our method, the transition between compressible andincompressible states as a unique explosion property is
consistently dealt with by grid-based simulation. Afterthe transition to low-speed flows, various previouslyproposed approaches for grid-based simulations can beapplied.
The rest of the paper is organized as follows. In Section 2 webriefly review the relevant related work. Section 3 presents
the details of our model and the algorithm. In Section 4 weshow our experimental results and discuss the advantages
and limitations of our approach. Section 5 presents conclusions and suggestions for future work.
2. Related Work
Modeling explosions. Yngve et al. [YOH00] proposed amodel in which compressible Navier-Stokes equations withsmall time steps are used to generate blast waves of explosions. Sewall et al. [SGTL09] introduced a finite volume
method for resolving shock waves that requires small time
steps to describe the compressible phenomena. In a later
study, a divergence modification was added to [Sta99] tohandle the reactive ignition described in [FOA03]. While[FOA03] can realize the drastic effects of the ignition, thehigh-speed flows generated by this method are hard to predict.
Takeshita et al. [TOT∗ 03] modeled flames of explosions
using the Lagrangian method. Ihm et al. [IKC04] and Kang

et al. [KJI07] extended the method of Takeshita et al. tocreate a model for explosive fluid by using chemical reaction processes. Gao et al. [GLHB09] proposed a hybrid
solver of Lagrangian and Eulerian methods to deal effectively with high- and low-speed fluids in one model. Kwatra
et al. [KGF10] proposed a method to physically calculate an
explosion’s compressible-to-incompressible transitions withlarger time steps in a semi-implicit method. Since those previous works determine the explosion behaviors only by thephysical parameters, it is difficult to control the complex explosion behaviors.

Control for fluid simulations. Treuille et al. [TMPS03]
proposed the density keyframing algorithm whereby thedriving forces to generate specific keyframes are optimized.
To follow up this research, Fattal and Lischeiski [FL04]
modeled the algorithm to avoid optimization. To get smoke
fluids to follow specified paths, Kim et al. [KMD06] proposed a path-based control method by generating velocity fields. Shi et al. [SY05] controlled target liquids, whichchange their shapes rapidly by external forces. Thuerey et al.
[TKRP09] proposed a control method by decomposing velocity fields into low- and high-frequency components withvarious control forces.
Pighin et al. [PCS04] edited buoyancy-driven fluids by
specifying the density at various keyframes. Beaudoin etal. [BPP01] and Lamorlette et al. [LF02] realized the behaviors of fire by using primitives or particle elements based on
physical properties without numerical simulations.
Although various types of fluid control models have been
proposed, unfortunately only a few methods for modeling
explosion fluids have been developed. The step size needs tobe small enough for high-speed explosions, as in the methodof Shi et al. [SY05]. This leads to huge computational costs
not suitable for user turnaround. Dobashi et al. [DSYA09]
first proposed a fluid control model specialized for explosion. However, this method’s primary purpose is to effectively generate the end flow as it reaches the destination,rather than to simulate the entire propagation of the flow
frame by frame.
c The Eurographics Association 2011.
⃝

G. Kawada & T. Kanai / Procedural Fluid Modeling of Explosion Phenomena Based on Physical Properties

Figure 2: Our algorithm. Pre-processing: flow control path converted from user-specified NURBS curve; Process 1 and Process
2: density and temperature fields and velocity fields calculated based on density propagation and opacity curves, respectively;
Process 3: pressure field determined based on the pressure propagation curve; Process 4: detonation state for drastic pressurechange and fuel combustion for large vortices; Process 5: grid-based simulation of the field calculations.
3. Model and Algorithm
3.1. Model
We define a model based on physical properties to generate the behaviors of flows in explosion phenomena. We consider a one-dimensional model as the physical model. This
one-dimensional model represents the propagation properties over time in a certain direction from the source to thedestination. Explosion flows are described by applying thisone-dimensional model to each of the control paths specified
by the user.
For the model to calculate the flows along the control
paths, we utilize the pressure magnitude curve, pressure
propagation curve, density propagation curve, and density opacity curve defined by [Bak73] and [HS94] (or see[GD77]) (Figure 1). The vertical axes are different for eachcurve, while the horizontal axes of all those curves describe
time. The pressure magnitude curve exhibits an important
property, and the other curves are derived from this curve.
In other words, the pressure obtained by the pressure magnitude curve determines other physical properties (the density
opacity, the velocity of density, the magnitude of pressure,and the velocity of pressure).
The magnitude of propagating pressure front p at the explosion can be approximated by the Friedlander equation[Bak73] as described below. We define p as the standard
pressure. The pressure magnitude curve (Figure 1 left) formularizes the following equation:

− tbt
t

d ,
i f t ≤ td ,
 p0 + Pp (1 − td )em p0 − 2P(t−t),
i f t > tddtn& t ≤ td + t2n ,p(t) =
2Pmtn p0 − tn (td + tn − t), i f t > td + 2& t ≤ td + tn ,p0 ,i f t > td + tn ,(1)
where p0 is the ambient pressure, Pp is the peak over presc The Eurographics Association 2011.
⃝

sure, Pm is the minimum negative pressure, td is the time
to reach p0 , tn is the time to reach p0 from td , and b is thedecreasing coefficient. This curve indicates that a huge pressure arises around the ignition, and that pressure gradually
decreases and converges before the flow reaches its destination. On the other hand, the density propagation curve describes the propagation velocity of flow (density), and thedensity opacity curve describes the opacity (density value)of flow. The pressure propagation curve describes the propagation velocity of standard pressure. The latter three curves
are formulated in later sections.
In our method, each explosion flow propagates on eachcontrol path until the destination time based on the propagation velocity magnitude of the flow is reached. This velocity magnitude is determined by the explosion curves represented by the one-dimensional model described above. The
distance from the ignition source to the destination, that is,
the length of each control path, is different, while we also
consider all flows to reach the ending points of the paths at
the same time. To this end, the propagation velocity magnitude of each flow is scaled based on the actual length of eachpath. Therefore, the longer the control path is, the larger theactual velocity magnitude becomes, while the relative velocity magnitude on each path between times is maintained.
In contrast, the explosion behaviors for the entire scene
domain are realized by the grid-based simulation. Namely,
our method combines the grid-based simulation with thegeneration of velocity, density, temperature and pressure
fields from multiple control paths. Therefore, consistency for
all explosion behaviors is maintained, while each explosion
flow is controlled.

3.2. Algorithm
Figure 2 shows our algorithm. There is one pre-processingin our algorithm. There are also five modeling processes per
step. The five processes are executed until the flows reachG. Kawada & T. Kanai / Procedural Fluid Modeling of Explosion Phenomena Based on Physical Properties

the destination, while only Process 5 is executed after theflows reach the destination.
• As Pre-processing, NURBS curves representing the target flows are specified by the user. These curves are automatically converted into flow control paths and pressure
control paths (Section 3.3).
• As Process 1, to realize the compressibility of the flow,
the density opacity during the current step is calculated.
The density and temperature fields at the region along theflow control paths are also generated (Section 3.4.1). Furthermore, fuel particles are randomly generated based on
the temperature (Section 3.4.6).
• As Process 2, the propagation velocity magnitude of thedensity during the current step is calculated, and velocity
fields at the region along the paths are generated. Vortex
particles are also randomly generated (Section 3.4.2).
• As Process 3, the velocity and the magnitude of the pressure front during the current step are calculated, and thepressure fields at the region along the pressure control
paths are generated (Section 3.4.3).
• As Process 4, by detecting the detonation state that willcause the drastic change in the pressure, the pressure fields
are amplified (Section 3.4.5). Large vortices are generated
from fuel particles by considering the combustion (Section 3.4.6).
• As Process 5, the whole grid domain is calculated by executing grid-based simulation [Sta99] for one step (Section
3.4.4).
Details of each process in our algorithm are described in thefollowing sections.

Figure 3: Left: Density Field. Right: Pressure field.
NURBS curve is automatically converted into a flow control
path (blue line in Figure 2) defined as a discretized curve
composed of a sequence of center points on the grid.
The pressure propagates faster than the flow, and the propagation distance of the pressure is also longer than that of
the flow. Therefore, the pressure control path (green line in
Figure 2) is generated by extending the curve along the averaged unit tangent vector of the curve’s last several points(black points in path extension in Figure 2).
T is discretized by n sets of time interval ∆t, and the computation for T is described as n steps. The density flows startfrom the explosion source at t = 0, propagate along the specified curves in each step, and reach the ends of the curves at
the n th step (time t = T ).
3.4. Process per Step

3.3. User Input and Control Path
The parameters for the pressure magnitude curve (Equation(1)) are determined by the user. Also, NURBS curves (browncurves in Figure 2) are specified to prescribe the directionsof the desired target flows. The beginning locations of those
NURBS curves describe the explosion sources, and the end
locations describe the destination points. Our method considers that the density to represent the flow proceeds along
this curve of width w. The propagation time of the flow fromthe source to the destination is T , which is the actual timecorresponding to the time tv when the propagation magnitude becomes 0 at the density propagation curve (Figure 1
right). Here, tv is automatically determined once td is obtained, and td < tv is always true. The user specifies w andT.

In this subsection, the process at step i (i = 0 . . . n), in otherwords, when time t = ti , is described. Five processes per
step, which are the density and temperature generation, velocity generation, pressure generation, pressure amplification, fuel combustion, and grid-based simulation, are executed in order.
3.4.1. Density and Temperature Field Generation
First of all, the density fields along the flow control paths
during step i are generated (Process 1 in Figure 2). Those
fields are calculated based on both the density propagation
curve to describe the density propagation velocity and thedensity opacity curve to describe the density value.
The velocity magnitude, v(t), at the density propagation
curve is described [HS94] as follows:1
cp(t)γ + 1 p(t) − 2v(t) =
(2)
(1 +
) ,
γp0
2γ p0

The pressure magnitude curve in Equation (1) determines
three curves in Equation (2), (3), and (6) described below.
The horizontal axis, t, is scaled by considering that the time
when density propagation magnitude becomes 0 in Equation(2) is T .

where c is the speed of sound in the air, and γ is the ratio of
the specific heats (1.4 at the ambient pressure).

Actual calculation is executed in the grids covering the entire domain to describe the target explosion. Therefore, each

The front-most location of the flow on the path during thecurrent step i is calculated based on this velocity magnitude,c The Eurographics Association 2011.
⃝

G. Kawada & T. Kanai / Procedural Fluid Modeling of Explosion Phenomena Based on Physical Properties

The direction velocity field (Figure 4 left) during the step i
is generated at each grid inside the sweep region along Si
in the same way as the density field. Namely, the direction
velocity ud is allocated at the each grid, and the velocity is
described as:
ud (G(g)) = V (ti )t(g),(5)

where g (Figure 4 left, light blue grid) is a grid on Si , G(g)is a set of grids (Figure 4 left, pink region) inside the circle
region with width w surrounding g, and t(g) is the unit tangent vector derived from g on the flow control path. V (ti ) is
the same value for all grids inside the sweep region.
Figure 4: Velocity field.
which describes only the relative value. Therefore, this velocity magnitude needs to be scaled to become the magnitude along each actual control path. Namely, we consider
that the flow propagates based on the velocity magnitude described by Equation (2), and then the propagation distancefrom the source until time ti is l(ti ) = ∑ik=1 v(tk )∆t. Therefore, the propagation distance until time T is described asl(T ). We also let the corresponding length of a control pathbe L, and then the propagation velocity magnitude on thev(t )Lactual control path becomes V (ti ) = l(Ti ) .
The density value ρ is determined based on the density
opacity curve and ρ at this curve is described [HS94] as follows:
2γp0 + (γ + 1)p(t)ρ(t) = ρ0,
2γp0 + (γ − 1)p(t)(3)

where ρ0 is the ambient density specified by the user.
The propagation distance of the flow during step i is
V (ti )∆t. We define the control path segment corresponding
to this distance as Si , and then the density fields are generated on the grids corresponding to the sweep region made by
the width w that surrounds Si (Figure 3 left, red region). Indeed, the uniform value ρ(ti ) is allocated at each grid inside
this sweep region.
The temperature fields are also generated in the same way
as the density fields. The uniform temperature value (specified by the user) is allocated at each grid inside the same
sweep region. This temperature value decreases at a constantrate as time passes after the ignition based on the principle ofthermal radiation. Those temperature fields computed procedurally are mainly utilized for generating large vortices by
combustion (described in Section 3.4.6).
3.4.2. Velocity Field Generation
The velocity field is computed by combining direction velocity field Ud and vortex velocity field Uv as follows (Process
2 in Figure 2):
U = Ud + Uv .
c The Eurographics Association 2011.
⃝(4)

The vortex velocity field Uv (Figure 4 right) is an important explosion property. We apply the vortex particle methodby [SRF05] to add fine details of whirling motions. The vortex particles with uniform magnitude are randomly generated at the grids inside the sweep region Si in the same way
as the direction velocity field. Namely, the vortex velocity
field whose center is a particle is added to this region. Particles also advect along the field as the step proceeds. Also,
vortex vectors of the field are rotated per step without changing its magnitude [SRF05]. This prevents the total magnitudes on the partial region from increasing as the step proceeds. More dynamic and complex behaviors of explosion
phenomena are realized by such vortex velocity fields.
3.4.3. Pressure Field Generation
The pressure field propagating along the control path is calculated according to the propagation velocity determined by
the pressure propagation curve (Process 3 in Figure 2). As
described in Section 3.3, since the pressure propagates fasterthan the density flow, the pressure field propagates more
forward than the density field. We consider the pressure topropagate until time td . The velocity magnitude, v p (t), at thepressure propagation curve is described [HS94] as follows:1v p (t) = c(1 +γ + 1 p(t) 2
) .
2γ p0

(6)

The actual velocity magnitude on the pressure control path,V p (ti ), is calculated in the same way as V (ti ), the actual velocity magnitude on the flow control path (Figure 3 right).
Namely, the scale adjustment is applied to v p based on boththe control path length and the propagation distance at thepressure propagation curve.
It is very important to note that the method to calculate thepressure field is totally different from the methods to calculate the density and velocity fields, as described previously.
The density and velocity fields are generated in the region
where the flow propagates only at the current step based on
the corresponding explosion curves. In contrast, the pressure distribution during the current step between the pressure
front and the pressure end (the source) is independent of theprevious frame’s distribution, which is calculated based on
the previous step. In other words, each frame’s distribution

G. Kawada & T. Kanai / Procedural Fluid Modeling of Explosion Phenomena Based on Physical Properties

is calculated based only on that frame’s information. Therefore, the pressure fields need to be calculated step by step byconsidering the entire distribution between the front and theend of the current pressure.
Our method to calculate the pressure field is describedp
as follows [WBB99]. First of all, we define S̃i , the pressure control path segment ranging from its top to its end,p
as shown in Figure 3 right. The distance for S̃i is thenpip
L (ti ) = ∑k=1 V (tk )∆t. The pressure value 0.4P(ti ) is generated at each grid inside the sweep region along the segment ranging from the end to the point corresponding to
0.5L p (ti ). P(ti ) is the adjusted value determined by considering the standard pressure, p(ti ), to be linearly† proportional to the velocity magnitude of the pressure. The scaled
pressure magnitude, P(ti ), is then obtained. For the rest ofp
the segment S̃i , we consider 0.4P(ti ) for the location ofp
0.5L (ti ) and P(ti ) for the front. In this case the pressure
field values inside the corresponding region are linearly interpolated.
3.4.4. Grid-Based Simulation for Fields
The density, temperature, velocity, and pressure are determined around the control paths as our algorithm proceeds toProcess 3. The grid fields for the entire domain are, however,not determined yet. Therefore, by applying grid-based simulation [Sta99] for one step, the density, temperature, velocity,and pressure fields are updated for the entire domain.
Our method updates each field by two-stage calculations,as follows. For the first stage, advection, pressure (projection), and diffusion terms in grid-based simulation [Sta99]
are solved for one step without utilizing the pressure fields
calculated in Section 3.4.3. The reason is that p(t) is included in the right term of Equation (2) and (3) for the density and velocity fields in Section 3.4.1 and 3.4.2. Therefore,
the density and velocity fields take the effects by the pressure
calculated in Section 3.4.3 into account.
For the second stage, first the pressure fields are calculated as described in Section 3.4.3. Those calculated fields
are used to determine the detonation state in Section 3.4.5.
The pressure field for the whole domain is generated by solving the pressure term of the grid-based simulation. This field
is mainly utilized for the visualization of pressure (shockwave).
In our method, Process 5 continues after time T while
all other processes are terminated. All fields are stably andconsistently solved from procedurally determined pieces of† According to Equation (6), the pressure magnitude needs to be
proportional to a square of the velocity magnitude of the pressure.
However, we found by our experiments that the pressure magnitudedifference between the paths becomes quite huge, so we apply thislinear relation equation.

Figure 5: Rankine-Hugoniot Relation Equations.
fields, especially for advection, and pressure terms. Therefore, the pressure calculated at the compressible state is
not added to the incompressible solver as a huge numerical burden. This is the advantage of our method in terms
of preventing the velocity field from diverging, as discussed
in [KGF10].
3.4.5. Considering Drastic Change of Pressure
The pressure is changed drastically in the domain where theshock wave pressure propagates. The state that yields thedrastic pressure change is called the detonation state. In thissubsection, we explain our method to realize the pressurechange caused by a strong detonation state. By consideringsuch a change in state, large vortices and changes over theentire flows are realized.
We use the Rankine-Hugoniot relation equations [Law06](see graphs in Figure 5) to express the drastic pressurechange, as follows:( p̂ +

γ−1
γ−1
γ−1
4γ
+ 2q(
)(η̂ −) =), (7)γ+1γ+1γ+1(γ + 1)2M2 = −p̂ − 1,γ(η̂ − 1)(8)

where p̂ denotes the pressure change rate between before
and after the pressure propagation, and η̂ denotes the change
rate for the inverse of the density (= 1/ρ). Both of theseequations describe the relationship between p̂ and η̂. Here,q denotes the amount of heat release per unit mass flux, andM denotes the Mach number (= v p (t)/c). Equation (7) represents two different curves (each is the Hugoniot curve)
depending on the q value. Equation (8) is also called theRayleight curve.
The detonation state when q = 0 is the strongest detonation state that yields the largest pressure change between before and after the propagation. We consider the intersection
point, (X0 ), of the Rayleight curve and the Hugoniot curvewhen q = 0 (solid red line in Figure 5).
In contrast, the detonation state when q > 0 yields lesspressure change, because the energy is partially used for heatgeneration. This state is considered to be the state of a usual
explosion. There are two intersection points (X1 , X2 ) for theRayleight curve and the Hugoniot curve when q > 0 (dottedc The Eurographics Association 2011.
⃝

G. Kawada & T. Kanai / Procedural Fluid Modeling of Explosion Phenomena Based on Physical Properties

red line in Figure 5). In general, the change rate for the pressure on X0 ( p̂X0 ) has a larger value than the change rates on
X1 , X2 ( p̂X1 , p̂X2 ).
In our method, the strong detonation state is realized by
the amplification of the pressure as follows: the pressure
2 p̂change rate r p̂ = p̂X +Xp̂0 X is multiplied to P(ti ) if M > 1 (the12condition to become the detonation state by considering theMach number). P(ti )r p̂ is replaced by the new pressure value
in Process 3 (Section 3.4.3). This routine is added as Process
4, as shown in Figure 2.
For the temperature field, the temperature value is amplified by using the same pressure change rate, r p̂ .
3.4.6. Fuel and Combustion
In this subsection, we consider the relationship between
fuel and combustion and realize large vortices, including theflows called turbulence. As an actual routine, fuel particles
are generated depending on fuel quantity. We define fuel particles as the seeds to generate the large vortices with the duration time determined by the temperature. In addition, thoseparticles advect themselves with the density flows.
We consider two sphere-shaped explosions, A and B. Let
the ignition fuel quantity for A and B be WA and WB , respectively, and the distance in a certain direction from the source
to the destination for A and B be dA and dB , respectively. The
( )1
dA 3
A
relationship W=
exists in those values [HS94]. We
WB
dB
now regard W as the generation probability of fuel particles,
in other words, the quantity to determine the percentage ofparticles for a certain region. Therefore, W is determined according to the control path length, d, by utilizing the above
relationship. Our method for generating fuel particles is described as follows. In the density field generation in Section
3.4.1, fuel particles are at the same time generated at theprobability rate of W , which is much smaller than the rate
for the vortex particles here.
Next, the duration time for the vortex is determined basedon the temperature of the grid where each fuel particle isgenerated. According to the Arrehenius equation, the energy,β2

K (= β1 e− Tem ), is calculated, where Tem denotes the temperature of the grid, β1 denotes the frequency factor coefficient, and β2 denotes the energy coefficient. The calculated
energy determines the magnitude of the vortex. This vortex
decreases at a constant rate as the time passes, and the duration time is the time when the vortex magnitude becomes 0
from the generation of the vortex.
Similarly this algorithm to determine the duration time
can also be applied to the vortex particles described in Section 3.4.2. We also simply realize the density transition offlame to smoke at each grid by considering the temperature
below a certain threshold specified by the user. To generatesmoke more specifically, density values can be added to thegrid with the temperature below the threshold.
c The Eurographics Association 2011.
⃝

Figure 6: Our method (one path). First row (without the
drastic pressure change or combustion). Second row (withonly the drastic pressure change). Third row (with only the
combustion). Fourth row (with both the drastic pressurechange and combustion).
4. Results and Discussion
Experimental results obtained by using our method areshown in this section. The results are generated using a PC
with a 2.8-GHz Intel Core i7-930.
Figure 6 shows the comparison with and without considering the drastic pressure change and combustion by using
one vertical control path specified from the bottom to the
top of the domain (1003 grids, Left: 10th frame (destinationtime), Right: 20th frame). The example shows that the quantity of the flow movements is increased only with the drasticpressure change (second rows). Also, the flow disturbancesby large vortices can only be seen in the example with the
combustion (third rows). Since we apply a 1D model to a
3D domain, as Figure 6 shows, the velocity fields include
some discontinuous artifacts, especially in the case of one
path compared to the case of multiple paths.
Figure 7 shows the comparison between the previous
method [FOA03] that uses only the grid-based method withvortex particles [SRF05] and our method (1003 grids). As

G. Kawada & T. Kanai / Procedural Fluid Modeling of Explosion Phenomena Based on Physical Properties

Figure 7: Comparison of our method with the previous
method (ground explosion No. 1 at the 50th frame, which isafter the destination). Left: our method (input curves). Right:
previous method [FOA03].

the figure shows, our method realizes similar propagation
behaviors from the source as those described using the previous method, while it is difficult to compare the methods
under exactly the same conditions. In addition, the example shows that our method models the explosions along thecontrol paths. One of our main contributions is high controllability in flow directions and specific target shapes, which
purely grid-based simulation cannot realize.
Figure 9 demonstrates the results using two different setsof control paths by our method. Both results are generated
with 2003 grids. The first two rows in the figure are examples
with the paths specified to all directions uniformly, to represent a hemisphere-shaped explosion, while last two rows are
the examples with the paths specified to only several specificdirections. All examples show our method to freely model
the explosions along the paths. The generation of explosions
along such path directions have not been addressed in a previous method, and this is the main advantage of our method.
Figure 8 shows a comparison on a camel model with andwithout the drastic pressure change. Both results are generated with 2003 grids and the same control paths. Regarding
the control paths, the NURBS curves to a head, front legs,
and back legs are specified by the user, while the curves
to vertices of the body are automatically generated. The top
row in the figure shows the example with the drastic pressurechange, and the bottom row in the figure shows the example without the drastic pressure change. The largely curved
paths in the top row show that the increased quantity of theflow makes the flows larger than the target shape, especiallyaround the neck. For these kinds of shapes, the algorithm toconsider geometric features of the paths has to be developed
as a future work.
The computational cost per frame is about 60 seconds inboth Figure 8 and Figure 9 with 2003 grids (10 seconds on
average for Process 1 - 4). Approximately 30 seconds haveto be added when the pressure term in the grid-based simulation is solved for the purpose of visualizing the pressure(shock wave). The computational cost to complete the processes of our method is quite low, while the computation isspent mostly on the pressure term for the visualization.

Figure 8: Our method (camel). Top row (with the drasticpressure change) left: 10th frame (destination time) and input curves. Top row right: 55th frame. Bottom row (without
the drastic pressure change) left: 10th frame. Bottom row
right: 55th frame.
5. Conclusions and Future WorkIn this paper, we proposed a method to procedurally model
explosion phenomena by considering physical properties
along the control paths specified by the user. The user’s intension can be taken into account by our method, and at thesame time explosion flows with complex behaviors can berealized by considering the propagations of the pressure anddensity flow, the fuel combustion and the detonation state torepresent the drastic pressure change.
As our future work, more refined algorithms for generating control paths will be considered. For example, the pathwidths can be automatically calculated based on the relationship between the target explosion shape and the specified curves, and tree-structural paths can also be considered
to describe the flow details. Automatic systems to optimize
the number of paths can also be developed by considering
the velocity fields among the paths. In addition, interactions
with objects have to be taken into account by considering thepressure reflection effects during the interactions. Finally,
we would like to reduce computational time by developing amore interactive system for users.
References[Bak73] BAKER W.: Explosions in Air. University of Texas Press,
1973. 1, 3
[BPP01] B EAUDOIN P., PAQUET S., P OULIN P.: Realistic andcontrollable fire simulation. In Proc. Graphics Interface (2001),
pp. 159–166. 2[DSYA09] D OBASHI Y., S ATO S., YAMAMOTO T., A NJYO K.:
Controlling explosion simulation. In SIGGRAPH ASIA 2009
Sketches (Yokohama, Japan, December 2009). 1, 2
[FL04] FATTAL R., L ISCHINSKI D.: Target-driven smoke animation. ACM Transactions on Graphics 23, 3 (2004), 441–448.
2
c The Eurographics Association 2011.
⃝

G. Kawada & T. Kanai / Procedural Fluid Modeling of Explosion Phenomena Based on Physical Properties[FOA03] F ELDMAN B. E., O’B RIEN J. F., A RIKAN O.: Animating suspended particle explosions. ACM Transactions on Graphics 22, 3 (2003), 708–715. 2, 7, 8
[GD77] G LASSTONE S., D OLAN P. J.: The Effects of Nuclear
Weapons. United States Department of Defense and the Energy
Research and Development Administration, 1977. 3[GLHB09] G AO Y., L I C.-F., H U S.-M., BARSKY B. A.: Simulating gaseous fluids with low and high speeds. Computer Graphics Forum (Proc. Pacific Graphics) 28, 7 (2009), 1845–1852. 2[HS94] H ETHERINGTON J. G., S MITH P. D.: Blast and Ballistic
Loading of Structures. Butterworth and Heinemann Ltd., Oxford,
1994. 1, 3, 4, 5, 7
[IKC04] I HM I., K ANG B., C HA D.: Animation of reactive
gaseous fluids through chemical kinetics. In Proc. Symposiumon Computer Animation (2004), pp. 203–212. 2[KGF10] K WATRA N., G RÉTARSSON J. T., F EDKIW R.: Practical animation of compressible flow for shock waves and related phenomena. In Proc. Symposium on Computer Animation(2010), pp. 207–215. 2, 6
[KJI07] K ANG B., JANG Y., I HM I.: Animation of chemicallyreactive fluids using a hybrid simulation method. In Proc. Symposium on Computer Animation (2007), pp. 199–208. 2[KMD06] K IM Y., M ACHIRAJU R., DAVID T.: Path-based control of smoke simulations. In Proc. Symposium on Computer
Animation (2006), pp. 33–42. 1, 2
[Law06] L AW C.: Comubustion Physics. Cambridge University
Press, 2006. 6[LF02] L AMORLETTE A., F OSTER N.: Structural modeling offlames for a production environment. ACM Transactions on
Graphics 21, 3 (2002), 729–735. 2[PCS04] P IGHIN F., C OHEN J. M., S HAH M.: Modeling andediting flows using advected radial basis functions. In Proc. Symposium on Computer Animation (2004), pp. 223–232. 2[SGTL09] S EWALL J., G ALOPPO N., T SANKOV G., L IN M.: Visual simulation of shockwaves. Graphical Models 71, 4 (2009),
126–138. 2[SRF05] S ELLE A., R ASMUSSEN N., F EDKIW R.: A vortex particle method for smoke, water and explosions. ACM Transactions
on Graphics 24, 3 (2005), 910–914. 5, 7
[Sta99] S TAM J.: Stable fluids. In Proc. ACM SIGGRAPH ’99(1999), pp. 121–128. 2, 4, 6
[SY05] S HI L., Y U Y.: Taming liquids for rapidly changing
targets. In Proc. Symposium on Computer Animation (2005),
pp. 229–236. 2[TKRP09] T HÜREY N., K EISER R., RÜDE U., PAULY M.:
Detail-preserving fluid control. Graphical Models 71, 6 (2009),
221–228. 2[TMPS03] T REUILLE A., M C NAMARA A., P OPOVI Ć Z., S TAM
J.: Keyframe control of smoke simulations. ACM Transactions
on Graphics 22, 3 (2003), 716–723. 2[TOT∗ 03] TAKESHITA D., OTA S., TAMURA M., F UJIMOTO T.,
M URAOKA K., C HIBA N.: Particle-based visual simulation of
explosive flames. In Proc. Pacific Graphics (2003), pp. 482–486.
2[WBB99] W INGERDEN K., B JERKETVEDT D., BAKKE J. R.:
Detonations in pipes and in the open. In Proc. the Petro-Chemical
Congress (1999). 6[YOH00] Y NGVE G. D., O’B RIEN J. F., H ODGINS J. K.: Animating explosions. In Proc. ACM SIGGRAPH 2000 (2000),
pp. 29–36. 2
c The Eurographics Association 2011.
⃝

Figure 9: Our method with both the drastic pressure change
and combustion. First row: ground explosion No 2, 10thframe (destination time) and input curves. Second row:
ground explosion No 2, 30th frame. Third row: ground explosion No 3, 25th frame (10th frame is destination time)and input curves. Fourth row: ground explosion No 3, 70thframe.