Animating Sand As A Fluid

Animating Sand as a Fluid
Yongning Zhu∗
Robert Bridson∗
University of British Columbia
University of British Columbia
Figure 1: The Stanford bunny is simulated as water and as sand.
Abstract
For small numbers of grains (say up to a few thousand) it is possible
to simulate each one as an individual rigid body, but scaling this up
We present a physics-based simulation method for animating sand.
to just a handful of sand is clearly infeasible; a sand dune containing
To allow for efﬁciently scaling up to large volumes of sand, we
trillions of grains is out of the question. Thus we take a continuum
abstract away the individual grains and think of the sand as a con-
approach, in essence squinting our eyes and pretending there are
tinuum. In particular we show that an existing water simulator can
no separate grains but rather the sand is a continuous material. The
be turned into a sand simulator with only a few small additions to
question is then what the constitutive laws of this continuum should
account for inter-grain and boundary friction.
be: how should it respond to force?
We also propose an alternative method for simulating ﬂuids. Our
The dynamics of granular materials have been studied extensively
core representation is a cloud of particles, which allows for accurate
in the ﬁeld of soil mechanics, but for the purposes of plausible ani-
and ﬂexible surface tracking and advection, but we use an auxiliary
mation we are willing to simplify their models drastically to allow
grid to efﬁciently enforce boundary conditions and incompressibil-
efﬁcient implementation. In section 3 we detail this simpliﬁcation,
ity. We further address the issue of reconstructing a surface from
arriving at a model which may be implemented by adding just a
particle data to render each frame.
little code to an existing water simulation.
CR Categories: I.3.5 [Computer Graphics]: Computational Ge-
We also present a new ﬂuid simulation method in section 4. As
ometry and Object Modeling—Physically based modeling
previous papers on simulating ﬂuids have noted, grids and particles
have complementary strengths and weaknesses. Here we combine
Keywords: sand, water, animation, physical simulation
the two approaches, using particles for our basic representation and
for advection, but auxiliary grids to compute all the spatial interac-
1
Introduction
tions (e.g. boundary conditions, incompressibility, and in the case
of sand, friction forces).
The motion of sand—how it ﬂows and also how it stabilizes—has
transﬁxed countless children at the beach or playground. More gen-
Our simulations only output particles that indicate where the ﬂuid
erally, granular materials such as sand, gravel, grains, soil, rubble,
is. To actually render the result, we need to reconstruct a smooth
ﬂour, sugar, and many more play an important role in the world,
surface that wraps around these particles. Section 5 explains our
thus we are interested in animating them.
work in this direction, and some of the advantages this framework
has over existing approaches.
∗email: {yzhu, rbridson}@cs.ubc.ca
2
Related Work
We brieﬂy review relevant papers, primarily in graphics, that pro-
vide a context for our research. In later sections we will refer to
those that explicitly guided our work.

Several authors have worked on sand animation, beginning with
also for tetrahedral meshes. Real-time elastic simulation includ-
Miller and Pearce[1989] who demonstrated a simple particle sys-
ing plasticity has been the focus of several papers (e.g. [M¨uller and
tem model with short-range interaction forces which could be tuned
Gross 2004]). Closest in spirit to the current paper, Goktekin et
to achieve (amongst other things) powder-like behavior. Later Lu-
al.[2004] added elastic forces and associated plastic ﬂow to a wa-
ciani et al.[1995] developed a similar particle system model specif-
ter solver, enabling animation of a wide variety of non-Newtonian
ically for granular materials using damped nonlinear springs at a
ﬂuids.
coarse length scale and a faster linear model at a ﬁne length scale.
Li and Moshell[1993] presented a dynamic height-ﬁeld simulation
For a scientiﬁc description of the physics of granular materials,
of soil based on the standard engineering Mohr-Coulomb constitu-
we refer the reader to Jaeger et al.[1996]. In the soil mechanics
tive model of granular materials. Sumner et al.[1998] took a similar
literature there is a long history of numerically simulating granu-
heightﬁeld approach with simple displacement and erosion rules to
lar materials using a elasto-plastic ﬁnite element formulation, with
model footprints, tracks, etc. Onoue and Nishita[2003] recently
Mohr-Coulomb or Drucker-Prager yield conditions and various
extended this to multi-valued heightﬁelds, allowing for some 3D
non-associated plastic ﬂow rules. We highlight the standard work
effects.
of Nayak and Zienkiewicz[1972], and the book by Smith and Grif-
ﬁths[1998] which contains code and detailed comments on FEM
As we mentioned above, one approach to granular materials
simulation of granular materials. Landslides and avalanches are two
is directly simulating the grains as interacting rigid bodies.
granular phenomena of particular interest, and generally have been
Milenkovic[1996] demonstrated a simulation of 1000 rigid spheres
studied using depth-averaged 2D equations introduced by Savage
falling through an hour-glass using optimization methods for re-
and Hutter[1989]. For alternatives that simulate the material at the
solving contact. Milenkovic and Schmidl[2001] improved the capa-
level of individual grains, Herrmann and Luding’s review[1998] is
bility of this algorithm, and Guendelman et al.[2003] further accel-
a good place to start.
erated rigid body simulation, showing 500 nonconvex rigid bodies
falling through a funnel.
Our new ﬂuid code traces its history back to the early particle-in-
cell (PIC) work of Harlow[1963] for compressible ﬂow. Shortly
Particles have been a popular technique in graphics for water simu-
thereafter Harlow and Welch[1965] invented the marker-and-cell
lation. Desbrun and Cani[1996] introduced Smooth Particle Hydro-
method for incompressible ﬂow, with its characteristic staggered
dynamics (see Monaghan[1992] for the standard review of SPH) to
grid discretization and marker particles for tracking a free surface—
the animation literature, demonstrating a viscous ﬂuid ﬂow simula-
the other fundamental elements of our algorithm.
PIC suf-
tion using particles alone. M¨uller et al.[2003] developed SPH fur-
fered from excessive numerical dissipation, which was cured by
ther for water simulation, achieving impressive interactive results.
the Fluid-Implicit-Particle (FLIP) method[Brackbill and Ruppel
Premoze et al.[2003] used a variation on SPH with an approximate
1986]; Kothe and Brackbill[1992] later worked on adapting FLIP
projection solve for their water simulations, but generated the ﬁnal
to incompressible ﬂow.
Compressible FLIP was also extended
rendering geometry with a grid-based level set. Particles were also
to a elasto-plastic ﬁnite element formulation, the Material Point
used by Takeshita et al.[2003] for ﬁre.
Method[Sulsky et al. 1995], which has been used to model gran-
ular materials at the level of individual grains[Bardenhagen et al.
There has been more work on animating water with grids, how-
2000] amongst other things. MPM was used by Konagai and Jo-
ever. Foster and Metaxas[1996] developed the ﬁrst grid-based fully
hansson[2001] for simulating large-scale (continuum) granular ma-
3D water simulation in graphics. Stam[1999] provided the semi-
terials, though only in 2D and at considerable computational ex-
Lagrangian advection method for faster simulation, but whose ex-
pense.
cessive numerical dissipation was later mitigated by Fedkiw et
al.[2001] with higher order interpolation and vorticity conﬁnement.
Foster and Fedkiw[2001] incorporated this into a water solver, and
added a level set—augmented by marker particles to counter mass
3
Sand Modeling
loss—for high quality surface tracking. To this model G´enevaux
et al.[2003] added two-way ﬂuid-solid interaction forces. Carlson
et al.[2002] added variable viscosity to a water solver, providing
3.1
Frictional Plasticity
beautiful simulations of wet sand dripping (achieved simply by in-
creased viscosity). Enright et al.[2002b] extended the Foster and
The standard approach to deﬁning the continuum behavior of sand
Fedkiw approach with particles on both sides of the interface and
and other granular (or “frictional”) materials is in terms of plastic
velocity extrapolation into the air. Losasso et al.[2004] extended
yielding. Suppose the stress tensor σ has mean stress σm = tr(σ )/3
this to dynamically adapted octree grids, providing much enhanced
and Von Mises equivalent shear stress ¯
σ =
resolution, and Rasmussen et al.[2004] improved the boundary con-
|σ −σmδ|F/√2 (where
ditions and velocity extrapolation while adding several other fea-
| · |F is the Frobenius norm). The Mohr-Coulomb1 law says the
tures important for visual effects production. Carlson et al.[2004]
material will not yield as long as
added better coupling between ﬂuid and rigid body simulations—
√
which we parenthetically note had its origins in scientiﬁc work on
3 ¯
σ < sinφσm
(1)
simulating ﬂow with sand grains. Hong and Kim[2003] and Taka-
hashi et al.[2003] introduced volume-of-ﬂuid algorithms for ani-
where φ is the friction angle. Heuristically, this is just the clas-
mating multi-phase ﬂow (water and air, as opposed to just the water
sic Coulomb static friction law generalized to tensors: the shear
of the free-surface ﬂows animated above).
stress σ −σmδ that tries to force particles to slide against one an-
other is annulled by friction if the pressure σm forcing the particles
For more general plastic ﬂow, most of the graphics work has dealt
against each other is proportionally strong enough. The coefﬁcient
with meshes.
Terzopoulos and Fleischer[1988] ﬁrst introduced
of friction µ commonly used in Coulomb friction is replaced here
plasticity to physics-based animation, with a simple 1D model ap-
by sin φ .
plied to their springs. O’Brien et al.[2002] added plasticity to a
fracture-capable tetrahedral-mesh ﬁnite element simulation to sup-
1Technically Mohr-Coulomb uses the Tresca deﬁnition of equivalent
port ductile fracture. Irving et al.[2004] introduced a more sophis-
shear stress, but since it is not smooth and poses numerical difﬁculties, most
ticated plasticity model along with their invertible ﬁnite elements,
people use Von Mises.

If the yield condition is met and the sand can start to ﬂow, a ﬂow
Consider one grid cell, of side length ∆x. The relative sliding ve-
rule is required. The simplest reasonable ﬂow direction is σ −σmδ.
locities between opposite faces in the cell are Vr = ∆xD. The mass
Heuristically this means the sand is allowed to slide against itself
of the cell is M = ρ∆x3. Ignoring other cells, the forces required
in the direction that the shearing force is pushing. Note that this
to stop all sliding motion in a time step ∆t are F = −VrM/∆t =
is nonassociated, i.e. not equal to the gradient of the yield condi-
−(∆xD)(ρ∆x3)/∆t, derived from integrating Newton’s law to get
tion. While associated ﬂow rules are simpler and are usually as-
the update formula V new = V + ∆tF/M. Dividing F by the cross-
sumed for plastic materials (and have been used exclusively be-
sectional area ∆x2 to get the stress gives:
fore in graphics[O’Brien et al. 2002; Goktekin et al. 2004; Irving
et al. 2004]), they are impossible here. The Mohr-Coulomb law
ρ
σ
D∆x2
compares the shear part of the stress to the dilatational part, unlike
rigid = − ∆
(3)
t
normal plastic materials which ignore the dilation part. Thus an as-
sociated ﬂow rule would allow the material to shear if and only if
We can check if this satisﬁes our yield inequality: if it does (i.e. if
it dilates proportionally—the more the sand ﬂows, the more it ex-
the sand can resist forces and inertia trying to make it slide), then
pands. Technically sand does dilate a little when it begins to ﬂow—
the material in that cell should be rigid at the end of the time step.
the particles need some elbow room to move past each other—but
This gives us the following algorithm. Every time step, after advec-
this dilation stops as soon as the sand is freely ﬂowing. This is a
tion, we do the usual water solver steps of adding gravity, applying
fairly subtle effect that is very difﬁcult to observe with the eye, so
boundary conditions, solving for pressure to enforce incompress-
we conﬁdently ignore it for animation.
ibility, and subtracting ∆t/ρ∇p from the intermediate velocity ﬁeld
Once plasticity has been deﬁned, the standard soil mechanics ap-
to make it incompressible. Then we evaluate the strain rate tensor
proach to simulation would be to add a simple elasticity relation,
D in each grid cell with standard central differences. If the result-
and use a ﬁnite element solver to simulate the sand. However,
ing σrigid from equation 3 satisﬁes the yield condition (using the
we would prefer to avoid the additional complexity and expense of
pressure we computed for incompressibility) we mark the grid cell
FEM (which would require numerical quadrature, an implicit solver
as rigid and store σrigid at that cell. Otherwise, we store the slid-
to avoid stiff time step restrictions due to elasticity, return mapping
ing frictional stress σ f from equation 2. We then ﬁnd all connected
for plasticity, etc.). For the purposes of plausible animation we will
groups of rigid cells, and project the velocity ﬁeld in each separate
make some further simplifying assumptions so we can retroﬁt sand
group to the space of rigid body motions (i.e. ﬁnd a translational
modeling into an existing water solver.
and rotational velocity which preserve the total linear and angular
momentum of the group). The remaining velocity values we up-
date with the frictional stress, using standard central differences for
3.2
A Simpliﬁed Model
u+ = ∆t/ρ∇ ·σf .
We ﬁrst ignore the nearly imperceptible elastic deformation regime
3.3
Frictional Boundary Conditions
(due to rock grains deforming) and the tiny volume changes that
sand undergoes as it starts and stops ﬂowing. Thus our domain
So far we have covered the friction within the sand, but there is
can be decomposed into regions which are rigidly moving (where
also the friction between the sand and other objects (e.g. the solid
the Mohr-Coulomb yield condition has not been met) and the rest
wall boundaries) to worry about. Our simple solution is to use the
where we have an incompressible shearing ﬂow.
friction formula of Bridson et al.[2002] when we enforce the u ·
We then further assume that the pressure required to make the en-
n ≥ 0 boundary condition on the grid. When we project out the
tire velocity ﬁeld incompressible will be similar to the true pressure
normal component of velocity, we reduce the tangential component
in the sand. This is certainly false: for example, it is well known
of velocity proportionately, clamping it to zero for the case of static
that a column of sand in a silo has a maximum pressure indepen-
friction:
µ
dent of the height of the column2 whereas in a column of water
u
|u·n|
T = max
0, 1 −
uT
(4)
the pressure increases linearly with height. However, we can only
|uT|
think of one case where this might be implausible in animation: the
where µ is the Coulomb friction coefﬁcient between the sand and
pressure limit means sand in an hour glass ﬂows at a constant rate
the wall.
independent of the amount of sand remaining (which is the reason
hour glasses work). We suggest our inaccurate results will still be
This is a crucial addition to a ﬂuid solver, since the boundary con-
plausible in most other cases.
ditions previously discussed in the animation literature either don’t
permit any kind of sliding (u = 0) which means sand sticks even to
Since we are neglecting elastic effects, we assume that where the
vertical surfaces, or always allow sliding (u · n ≥ 0, possibly with
sand is ﬂowing the frictional stress is simply
a viscous reduction of tangential velocity) which means sand piles
will never be stable. Compare ﬁgure 2, which includes friction, to
σ
D
ﬁgure 1 which doesn’t.
f = −sinφ p
(2)
1/3|D|F
where D = (∇u + ∇uT )/2 is the strain rate. That is, the frictional
3.4
Cohesion
stress is the point on the yield surface that most directly resists the
sliding.
A common enhancement to the basic Mohr-Coulomb condition is
the addition of cohesion:
We ﬁnally need a way of determining the yield condition (when
the sand should be rigid) without tracking elastic stress and strain.
√3 ¯σ < sinφσm+c
(5)
2This is due to the force that resists the weight of the sand above being
where c > 0 is the cohesion coefﬁcient. This is appropriate for
transferred to the walls of the silo by friction—sometimes causing silos to
soils or other slightly sticky materials that can support some tension
unexpectedly collapse as a result.
before yielding.

Figure 2: The sand bunny with a boundary friction coefﬁcient µ = 0.6: compare to ﬁgure 1 where zero boundary friction was used.
We have found that including a very small amount of cohesion im-
to represent the water surface, quickly rounding off any small fea-
proves our results for supposedly cohesion-less materials like sand.
tures. The most attractive alternative to level sets is volume-of-ﬂuid
In what should be a stable sand pile, small numerical errors can al-
(VOF), which conserves water up to ﬂoating-point precision: how-
low some slippage; a tiny amount of cohesion counters this error
ever it has difﬁculties maintaining an accurate but sharply-deﬁned
and makes the sand stable again, without visibly effecting the ﬂow
surface. Coupling level sets with VOF is possible[Sussman 2003]
when the sand should be moving.
but at a signiﬁcant implementation and computational cost.
However, for modeling soils with their true cohesion coefﬁcient,
On the other hand, particles trivially handle advection with excel-
our method is too stable: obviously unstable structures are implau-
lent accuracy—simply by letting them ﬂow through the velocity
sibly held rigid. We will investigate this issue in the future.
ﬁeld using ODE solvers—but have difﬁculties with pressure and
the incompressibility condition, often necessitating smaller than de-
sired time steps for stability. SPH methods in particular can’t tol-
4
Fluid Simulation Revisited
erate nonuniform particle spacing, which can be difﬁcult to enforce
throughout the entire length of a simulation.
4.1
Grids and Particles
Realizing the strengths and weaknesses of the two approaches com-
plement each other Foster and Fedkiw[2001] and later Enright et
There are currently two main approaches to simulating fully three-
al.[2002b] augmented a grid-based method with marker particles
dimensional water in computer graphics: Eulerian grids and La-
to correct the errors in the grid-based level set advection. This
grangian particles.
Grid-based methods (recent papers include
particle-level set method, most recently implemented on an adap-
[Carlson et al. 2004; Goktekin et al. 2004; Losasso et al. 2004])
tive octree[Losasso et al. 2004], has produced the highest ﬁdelity
store velocity, pressure, some sort of indicator of where the ﬂuid
water animations in the ﬁeld. However, we believe particles can be
is or isn’t, and any additional variables on a ﬁxed grid. Usually a
exploited even further, simplifying and accelerating the implemen-
staggered “MAC” grid is used, allowing simple and stable pressure
tation, and affording some new beneﬁts in ﬂexibility.
solves. Particle-based methods, exempliﬁed by Smoothed Particle
Hydrodynamics (see [Monaghan 1992; M¨uller et al. 2003; Premoze
et al. 2003] for example), dispense with grids except perhaps for ac-
celerating geometric searches. Instead ﬂuid variables are stored on
4.2
Particle-in-Cell Methods
particles which represent actual chunks of ﬂuid, and the motion of
the ﬂuid is achieved simply by moving the particles themselves.
The Lagrangian form of the Navier-Stokes equations (sometimes
Continuing the graphics tradition of recalling early computational
using a compressible formulation with a ﬁctitious equation of state)
ﬂuid dynamics research, we return to the Particle-in-Cell (PIC)
is used to calculate the interaction forces on the particles.
method of Harlow[1963]. This was an early approach to simulating
compressible ﬂow that handled advection with particles, but every-
The primary strength of the grid-based methods is the simplicity of
thing else on a grid. At each time step, the ﬂuid variables at a grid
the discretization and solution of the incompressibility condition,
point were initialized as a weighted average of nearby particle val-
which forms the core of the ﬂuid solver. Unfortunately, grid-based
ues, then updated on the grid with the non-advection part of the
methods have difﬁculties with the advection part of the equations.
equations. The new particle values were then interpolated from the
The currently favored approach in graphics, semi-Lagrangian ad-
updated grid values, and ﬁnally the particles moved according to
vection, suffers from excessive numerical dissipation due to accu-
the grid velocity ﬁeld.
mulated interpolation errors. To counter this nonphysical effect,
a nonphysical term such as vorticity conﬁnement must be added
The major problem with PIC was the excessive numerical diffu-
(at the expense of conserving angular or even linear momentum).
sion caused by repeatedly averaging and interpolating the ﬂuid
While high resolution methods from scientiﬁc computing can ac-
variables. Brackbill and Ruppel[1986] cured this with the Fluid-
curately solve for the advection of a well-resolved velocity ﬁeld
Implicit-Particle (FLIP) method, which achieved “an almost total
through the ﬂuid, their implementation isn’t entirely straightfor-
absence of numerical dissipation and the ability to represent large
ward: their wide stencils make life difﬁcult on coarse animation
variations in data.” The crucial change was to make the particles the
grids that routinely represent features only one or two grid cells
fundamental representation of the ﬂuid, and use the auxiliary grid
thick. Moreover, even ﬁfth-order accurate methods fail to accu-
simply to increment the particle variables according to the change
rately advect level sets[Enright et al. 2002a] as are commonly used
computed on the grid.

We have adapted PIC and FLIP to incompressible ﬂow3 as fol-
lows:
• Initialize particle positions and velocities
• For each time step:
• At each staggered MAC grid node, compute a weighted av-
erage of the nearby particle velocities
• For FLIP: Save the grid velocities.
• Do all the non-advection steps of a standard water simulator
on the grid.
• For FLIP: Subtract the new grid velocities from the saved
velocities, then add the interpolated difference to each par-
Figure 3: FLIP vs. PIC velocity update for the same simulation.
ticle.
Notice the small-scale velocities preserved by FLIP but smoothed
away by PIC.
• For PIC: Interpolate the new grid velocity to the particles.
• Move particles through the grid velocity ﬁeld with an ODE
solver, making sure to push them outside of solid wall
4.2.3
Solving on the Grid
boundaries.
• Write the particle positions to output.
We ﬁrst add the acceleration due to gravity to the grid velocities.
Observe there is no longer any need to implement grid-based ad-
We then construct a distance ﬁeld φ (x) in the unmarked (non-ﬂuid)
vection, or the matching schemes such as vorticity conﬁnement and
part of the grid using fast sweeping[Zhao 2005] and use that to ex-
particle-level set to counter numerical dissipation.
tend the velocity ﬁeld outside the ﬂuid with the PDE ∇u ·∇φ = 0
(also solved with fast sweeping). We enforce boundary conditions
and incompressibility as in Enright et al.[2002b], then extend the
new velocity ﬁeld again using fast sweeping (as we found it to
4.2.1
Initializing Particles
be marginally faster and simpler to implement than fast march-
ing[Adalsteinsson and Sethian 1999]).
We have found a reasonable effective practice is to create 8 particles
in every grid cell, randomly jittered from their 2 × 2 × 2 subgrid
positions to avoid aliasing when the ﬂow is underresolved at the
4.2.4
Updating Particle Velocities
simulation resolution. With less particles per grid cell, we tend to
run into occasional “gaps” in the ﬂow, and more particles allow for
At each particle, we trilinearly interpolate either the velocity (PIC)
too much noise. To help with surface reconstruction later (section
or the change in velocity (FLIP) recorded at the surrounding eight
5) we reposition particles that lie near the initial water surface (say
grid points. For viscous ﬂows, such as sand, the additional numer-
within one grid cell width) to be exactly half a grid cell away from
ical viscosity of PIC is beneﬁcial; for inviscid ﬂows, such as water,
the surface.
FLIP is preferable. See ﬁgure 3 for a 2D view of the differences
between PIC and FLIP. A weighted average of the two can be used
to tune just how much numerical viscosity is desired (of course, a
4.2.2
Transferring to the Grid
viscosity step on the grid in concert with PIC would be required for
highly viscous ﬂuids).
Each grid point takes a weighted average of the nearby particles.
We deﬁne “near” as being contained in the cube of twice the grid
cell width centered on the grid point. (Note that on a staggered
4.2.5
Moving Particles
MAC grid, the different components of velocity will have offset
neighborhoods.) The weight of a particle in this cube is the stan-
Once we have a velocity ﬁeld deﬁned on the grid (extrapolated to
dard trilinear weighting. We also mark the grid cells that contain at
exist everywhere) we can move the particles around in it. We use a
least one particle. In future work we will investigate using a more
simple RK2 ODE solver with ﬁve substeps each of which is limited
accurate indicator, e.g. reconstructing a signed distance ﬁeld on the
by the CFL condition (so that a particle moves less than one grid cell
grid from distance values stored on the particles. This would allow
in each substep), similar to Enright et al.[2002b]. We additionally
us to easily implement a second-order accurate free surface bound-
detect when particles penetrate solid wall boundaries due simply to
ary condition[Enright et al. 2003] which signiﬁcantly reduces grid
truncation error in RK2 and the grid-based velocity ﬁeld, and move
artifacts.
them in the normal direction back to just outside the solid, to avoid
We also note that the grid in our simulation is purely an auxiliary
the occasional stuck-particle artifact this would otherwise cause.
structure to the fundamental particle representation. In particular,
we are not required to use the same grid every time step. An obvi-
ous optimization which we have not yet implemented—rendering
5
Surface Reconstruction from Particles
is currently our bottleneck—is to deﬁne the grid according to the
bounding box of the particles every time step. This also means we
have no a priori limits on the simulation domain. We plan in the
Our simulations output the positions of the particles that deﬁne the
future to detect when clusters of particles have separated and use
location of the ﬂuid. For high quality rendering we need to recon-
separate bounding grids for them, for signiﬁcantly improved efﬁ-
struct a surface that wraps around the particles. Of course we can
ciency.
give up on direct reconstruction, e.g. running a level set simulation
guided by the particles[Premoze et al. 2003]. While this is an effec-
3The one precedent for this that we could ﬁnd is a reference to an un-
tive solution, we believe a fully particle-based reconstruction can
published manuscript[Kothe and Brackbill 1992].
have signiﬁcant advantages.

frame; in the absence of a fast method of computing these to the de-
sired precision, we currently ﬁx all the particle radii to the constant
average particle spacing and simply adjust our initial partial posi-
tions so that the surface particles are exactly this distance from the
surface. A small amount of additional grid smoothing, restricted to
decreasing φ to avoid destroying features, reduces bump artifacts at
later frames.
On the other hand, we do enjoy signiﬁcant advantages over recon-
struction methods that are tied to level set simulations. The cost per
frame is quite low—even in our unoptimized version which calcu-
lates and writes out a full 2503 grid, we average 40–50 seconds a
Figure 4: Comparison of our implicit function and blobbies for
frame on a 2GHz G5, with the bulk of the time spent on I/O. More-
matching a perfect circle deﬁned by the points inside. The blobby
over, every frame is independent, so this is easily farmed out to
is the outer curve, ours is the inner curve.
multiple CPU’s and machines running in parallel.
If we can eliminate the need for grid-based smoothing in the future,
Naturally the ﬁrst approach that comes to mind is blobbies[Blinn
perhaps adapting an MLS approach such as in Shen et al.[2004],
1982]. For irregularly shaped blobs containing only a few particles,
we could do the surface reconstruction on the ﬂy during rendering.
this works beautifully. Unfortunately, it seems unable to match a
Apart from speed, the biggest advantage this would bring would
surface such as a ﬂat plane, a cone, or a large sphere from a large
be accurate motion blur for ﬂuids. The current technique for gen-
quantity of irregularly spaced particles—as we might well see in at
erating intermediate surfaces between frames (for a Monte Carlo
least the initial conditions of a simulation. Bumps relating to the
motion blur solution) is to simply interpolate between level set grid
particle distribution are always visible. A small improvement to
values[Enright et al. 2002b]. However, this approach destroys small
this was suggested in [M¨uller et al. 2003], where the contribution
features that move further than their width in one frame: the inter-
from a particle was divided by the SPH estimate of density. This
polated values may all be positive. Unfortunately it’s exactly these
overcomes some scaling issues but does not signiﬁcantly reduce
small and fast-moving features that most demand motion blur. With
the bumps on what should be ﬂat surfaces.
particle-based surface reconstruction, the positions of the particles
can be interpolated instead, so that features are preserved at inter-
We thus take a different approach, guided by the principle that we
mediate times.
should exactly reconstruct the signed distance ﬁeld of an isolated
particle: for a single particle x0 with radius r0, our implicit function
must be:
φ(
6
Examples
x) = |x−x0|−r0
(6)
To generalize this, we use the same formula with x0 replaced by a
weighted average of the nearby particle positions and r
We used pbrt[Pharr and Humphreys 2004] for rendering. For the
0 replaced a
weighted average of their radii:
textured sand shading, we blended together a volumetric texture
advected around by the simulation particles, similar to the approach
φ(
of Lamorlette[2001] for ﬁreballs.
x)
=
|x− ¯x|− ¯r
(7)
¯
x
=
∑wixi
(8)
The bunny example in ﬁgures 1 and 2 were simulated with 269,322
i
particles on a 1003 grid, taking approximately 6 seconds per frame
¯r
=
∑w
on a 2Ghz G5 workstation. We believe optimizing the particle
iri
(9)
i
transfer code and using BLAS routines for the linear solve would
substantially improve this performance. Unoptimized smoothing
k(
w
|x−xi|/R)
i
=
and text I/O cause the surface reconstruction on a 2503 grid to take
∑
(10)
j k(|x − xj|/R)
40–50 seconds per frame.
where k is a suitable kernel function that smoothly drops to zero—
The column test in ﬁgures 5 and 6 was simulated with 433,479 par-
we used k(s) = max(0, (1 −s2)3) since it avoids the need for square
ticles on a 100 ×60×60 grid, taking approximately 12 seconds per
roots and is reasonably well-shaped—and where R is the radius of
frame. Our cost is essentially linearly proportional to the number
the neighborhood we consider around x. Typically we choose R to
of particles (or equivalently, occupied grid cells).
be twice the average particle spacing. As long as the particle radii
are bounded away from zero (say at least 1/2 the particle spacing)
we have found this formula gives excellent agreement with ﬂat or
smooth surfaces. See ﬁgure 4 for a comparison with blobbies using
7
Conclusion
the same kernel function.
The most signiﬁcant problem with this deﬁnition is artifacts in con-
We have presented a method for converting an existing ﬂuid solver
cave regions: spurious blobs of surface can appear, since ¯
x may er-
into one capable of plausibly animating granular materials such as
roneously end up outside the surface in concavities. However, since
sand. In addition, we developed a new ﬂuid solver that combines
these artifacts are signiﬁcantly smaller than the particle radii we can
the strengths of both particles and grids, offering enhanced ﬂex-
easily remove them without destroying true features by sampling
ibility and efﬁciency. We offered a new method for reconstruct-
φ(x) on a higher resolution grid and then doing a simple smoothing
ing implicit surfaces from particles. Looking toward the future, we
pass.
plan to more aggressively exploit the optimizations available with
PIC/FLIP (e.g. using multiple bounding grids), increase the accu-
A secondary problem is that we require the radii to be accurate es-
racy of our boundary conditions, and implement motion blur of the
timates of distance to the surface. This is nontrivial after the ﬁrst
reconstructed surface.

Figure 5: A column of regular liquid is released.
BLINN, J. 1982. A generalization of algebraic surface drawing.
ACM Trans. Graph. 1, 3, 235–256.
BRACKBILL, J. U., AND RUPPEL, H. M. 1986. FLIP: a method
for adaptively zoned, particle-in-cell calculuations of ﬂuid ﬂows
in two dimensions. J. Comp. Phys. 65, 314–343.
BRIDSON, R., FEDKIW, R., AND ANDERSON, J. 2002. Robust
treatment of collisions, contact and friction for cloth animation.
ACM Trans. Graph. (Proc. SIGGRAPH) 21, 594–603.
CARLSON, M., MUCHA, P., VAN HORN III, R., AND TURK,
G.
2002.
Melting and ﬂowing.
In Proc. ACM SIG-
GRAPH/Eurographics Symp. Comp. Anim., 167–174.
CARLSON, M., MUCHA, P. J., AND TURK, G. 2004. Rigid ﬂuid:
animating the interplay between rigid bodies and ﬂuid. ACM
Trans. Graph. (Proc. SIGGRAPH) 23, 377–384.
DESBRUN, M., AND CANI, M.-P. 1996. Smoothed particles: A
new paradigm for animating highly deformable bodies. In Com-
put. Anim. and Sim. ’96 (Proc. of EG Workshop on Anim. and
Sim.), Springer-Verlag, R. Boulic and G. Hegron, Eds., 61–76.
Published under the name Marie-Paule Gascuel.
ENRIGHT, D., FEDKIW, R., FERZIGER, J., AND MITCHELL, I.
2002. A hybrid particle level set method for improved interface
capturing. J. Comp. Phys. 183, 83–116.
ENRIGHT, D., MARSCHNER, S., AND FEDKIW, R. 2002. Ani-
mation and rendering of complex water surfaces. ACM Trans.
Figure 6: A column of granular material is released.
Graph. (Proc. SIGGRAPH) 21, 3, 736–744.
ENRIGHT, D., NGUYEN, D., GIBOU, F., AND FEDKIW, R. 2003.
8
Acknowledgements
Using the particle level set method and a second order accurate
pressure boundary condition for free surface ﬂows. In Proc. 4th
ASME-JSME Joint Fluids Eng. Conf., no. FEDSM2003–45144,
This work was in part supported by a grant from the Natural Sci-
ASME.
ences and Engineering Research Council of Canada. We would like
to thank Dr. Dawn Shuttle for her help with granular constitutive
FEDKIW, R., STAM, J., AND JENSEN, H. 2001. Visual simulation
models.
of smoke. In Proc. SIGGRAPH, 15–22.
FOSTER, N., AND FEDKIW, R. 2001. Practical animation of liq-
uids. In Proc. SIGGRAPH, 23–30.
References
FOSTER, N., AND METAXAS, D. 1996. Realistic animation of
liquids. Graph. Models and Image Processing 58, 471–483.
ADALSTEINSSON, D., AND SETHIAN, J. 1999. The fast con-
struction of extension velocities in level set methods. J. Comput.
G ´ENEVAUX, O., HABIBI, A., AND DISCHLER, J.-M. 2003. Sim-
Phys. 148, 2–22.
ulating ﬂuid-solid interaction. In Graphics Interface, 31–38.
BARDENHAGEN, S. G., BRACKBILL, J. U., AND SULSKY, D.
GOKTEKIN, T. G., BARGTEIL, A. W., AND O’BRIEN, J. F. 2004.
2000. The material-point method for granular materials. Com-
A method for animating viscoelastic ﬂuids. ACM Trans. Graph.
put. Methods Appl. Mech. Engrg. 187, 529–541.
(Proc. SIGGRAPH) 23, 463–468.

GUENDELMAN, E., BRIDSON, R., AND FEDKIW, R. 2003. Non-
NAYAK, G. C., AND ZIENKIEWICZ, O. C. 1972. Elasto-plastic
convex rigid bodies with stacking. ACM Trans. Graph. (Proc.
stress analysis. A generalization for various constitutive relations
SIGGRAPH) 22, 3, 871–878.
including strain softening. Int. J. Num. Meth. Eng. 5, 113–135.
HARLOW, F., AND WELCH, J.
1965.
Numerical Calculation
O’BRIEN, J., BARGTEIL, A., AND HODGINS, J. 2002. Graphical
of Time-Dependent Viscous Incompressible Flow of Fluid with
modeling of ductile fracture. ACM Trans. Graph. (Proc. SIG-
Free Surface. Phys. Fluids 8, 2182–2189.
GRAPH) 21, 291–294.
H
ONOUE, K., AND NISHITA, T. 2003. Virtual sandbox. In Paciﬁc
ARLOW, F. H. 1963. The particle-in-cell method for numerical
solution of problems in ﬂuid dynamics. In Experimental arith-
Graphics, 252–262.
metic, high-speed computations and mathematics.
PHARR, M., AND HUMPHREYS, G. 2004. Physically Based Ren-
dering: From Theory to Implementation. Morgan Kaufmann.
HERRMANN, H. J., AND LUDING, S. 1998. Modeling granular
media on the computer. Continuum Mech. Therm. 10, 189–231.
PREMOZE, S., TASDIZEN, T., BIGLER, J., LEFOHN, A., AND
WHITAKER, R. 2003. Particle–based simulation of ﬂuids. In
HONG, J.-M., AND KIM, C.-H. 2003. Animation of bubbles in
Comp. Graph. Forum (Eurographics Proc.), vol. 22, 401–410.
liquid. Comp. Graph. Forum (Eurographics Proc.) 22, 3, 253–
262.
RASMUSSEN, N., ENRIGHT, D., NGUYEN, D., MARINO, S.,
SUMNER, N., GEIGER, W., HOON, S., AND FEDKIW, R.
IRVING, G., TERAN, J., AND FEDKIW, R. 2004. Invertible ﬁnite
2004.
Directible photorealistic liquids.
In Proc. ACM SIG-
elements for robust simulation of large deformation. In Proc.
GRAPH/Eurographics Symp. Comp. Anim., 193–202.
ACM SIGGRAPH/Eurographics Symp. Comp. Anim., 131–140.
SAVAGE, S. B., AND HUTTER, K. 1989. The motion of a ﬁnite
JAEGER, H. M., NAGEL, S. R., AND BEHRINGER, R. P. 1996.
mass of granular material down a rough incline. J. Flui Mech.
Granular solids, liquids, and gases. Rev. Mod. Phys. 68, 4, 1259–
199, 177–215.
1273.
SHEN, C., O’BRIEN, J. F., AND SHEWCHUK, J. R. 2004. Inter-
KONAGAI, K., AND JOHANSSON, J.
2001. Two dimensional
polating and approximating implicit surfaces from polygon soup.
Lagrangian particle ﬁnite-difference method for modeling large
ACM Trans. Graph. (Proc. SIGGRAPH) 23, 896–904.
soil deformations. Structural Eng./Earthquake Eng., JSCE 18,
SMITH, I. M., AND GRIFFITHS, D. V. 1998. Programming the
2, 105s–110s.
Finite Element Method. J. Wiley & Sons.
KOTHE, D. B., AND BRACKBILL, J. U.
1992.
FLIP-INC: a
STAM, J. 1999. Stable ﬂuids. In Proc. SIGGRAPH, 121–128.
particle-in-cell method for incompressible ﬂows. Unpublished
manuscript.
SULSKY, D., ZHOU, S.-J., AND SCHREYER, H. L. 1995. Ap-
plication of particle-in-cell method to solid mechanics. Comp.
LAMORLETTE, A. 2001. Shrek effects—ﬂames and dragon ﬁre-
Phys. Comm. 87, 236–252.
balls. SIGGRAPH Course Notes, Course 19, 55–66.
SUMNER, R. W., O’BRIEN, J. F., AND HODGINS, J. K. 1998.
LI, X., AND MOSHELL, J. M. 1993. Modeling soil: Realtime
Animating sand, mud, and snow. In Graphics Interface, 125–
dynamic models for soil slippage and manipulation. In Proc.
132.
SIGGRAPH, 361–368.
SUSSMAN, M. 2003. A second order coupled level set and volume-
LOSASSO, F., GIBOU, F., AND FEDKIW, R. 2004. Simulating wa-
of-ﬂuid method for computing growth and collapse of vapor bub-
ter and smoke with an octree data structure. ACM Trans. Graph.
bles. J. Comp. Phys. 187, 110–136.
(Proc. SIGGRAPH) 23, 457–462.
TAKAHASHI, T., FUJII, H., KUNIMATSU, A., HIWADA, K.,
S
L
AITO, T., TANAKA, K., AND UEKI, H. 2003. Realistic an-
UCIANI, A., HABIBI, A., AND MANZOTTI, E. 1995. A multi-
imation of ﬂuid with splash and foam. Comp. Graph. Forum
scale physical model of granular materials. In Graphics Inter-
(Eurographics Proc.) 22, 3, 391–400.
face, 136–146.
TAKESHITA, D., OTA, S., TAMURA, M., FUJIMOTO, T., AND
MILENKOVIC, V. J., AND SCHMIDL, H. 2001. Optimization-
CHIBA, N. 2003. Visual simulation of explosive ﬂames. In
based animation. In Proc. SIGGRAPH, 37–46.
Paciﬁc Graphics, 482–486.
MILENKOVIC, V. J. 1996. Position-based physics: simulating the
TERZOPOULOS, D., AND FLEISCHER, K. 1988. Modeling in-
motion of many highly interacting spheres and polyhedra. In
elastic deformation: viscoelasticity, plasticity, fracture. Proc.
Proc. SIGGRAPH, 129–136.
SIGGRAPH, 269–278.
MILLER, G., AND PEARCE, A. 1989. Globular dynamics: a con-
ZHAO, H. 2005. A fast sweeping method for Eikonal equations.
nected particle system for animating viscous ﬂuids. In Comput.
Math. Comp. 74, 603–627.
& Graphics, vol. 13, 305–309.
MONAGHAN, J. J. 1992. Smoothed particle hydrodynamics. Annu.
Rev. Astron. Astrophys. 30, 543–574.
M ¨
ULLER, M., AND GROSS, M. 2004. Interactive virtual materials.
In Graphics Interface.
M ¨
ULLER, M., CHARYPAR, D., AND GROSS, M. 2003. Particle-
based ﬂuid simulation for interactive applications. In Proc. ACM
SIGGRAPH/Eurographics Symp. Comp. Anim., 154–159.