Animating Sand as a Fluid
Robert Bridson∗
University of British Columbia

Yongning Zhu∗
University of British Columbia

Figure 1: The Stanford bunny is simulated as water and as sand.

Abstract
We present a physics-based simulation method for animating sand.
To allow for efficiently scaling up to large volumes of sand, weabstract away the individual grains and think of the sand as a continuum. In particular we show that an existing water simulator canbe turned into a sand simulator with only a few small additions toaccount for inter-grain and boundary friction.

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
to just a handful of sand is clearly infeasible; a sand dune containing
trillions of grains is out of the question. Thus we take a continuum
approach, in essence squinting our eyes and pretending there are
no separate grains but rather the sand is a continuous material. The
question is then what the constitutive laws of this continuum should
be: how should it respond to force?

We also propose an alternative method for simulating fluids. Our
core representation is a cloud of particles, which allows for accurate
and flexible surface tracking and advection, but we use an auxiliary
grid to efficiently enforce boundary conditions and incompressibility. We further address the issue of reconstructing a surface fromparticle data to render each frame.

The dynamics of granular materials have been studied extensively
in the field of soil mechanics, but for the purposes of plausible animation we are willing to simplify their models drastically to allow
efficient implementation. In section 3 we detail this simplification,
arriving at a model which may be implemented by adding just alittle code to an existing water simulation.

CR Categories: I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling—Physically based modeling

We also present a new fluid simulation method in section 4. As
previous papers on simulating fluids have noted, grids and particles
have complementary strengths and weaknesses. Here we combine
the two approaches, using particles for our basic representation andfor advection, but auxiliary grids to compute all the spatial interactions (e.g. boundary conditions, incompressibility, and in the caseof sand, friction forces).

Keywords: sand, water, animation, physical simulation

1

Introduction

The motion of sand—how it flows and also how it stabilizes—hastransfixed countless children at the beach or playground. More generally, granular materials such as sand, gravel, grains, soil, rubble,
flour, sugar, and many more play an important role in the world,
thus we are interested in animating them.

Our simulations only output particles that indicate where the fluid
is. To actually render the result, we need to reconstruct a smoothsurface that wraps around these particles. Section 5 explains ourwork in this direction, and some of the advantages this frameworkhas over existing approaches.

∗ email: {yzhu, rbridson}@cs.ubc.ca

2

Related Work

We briefly review relevant papers, primarily in graphics, that provide a context for our research. In later sections we will refer tothose that explicitly guided our work.

Several authors have worked on sand animation, beginning with
Miller and Pearce[1989] who demonstrated a simple particle system model with short-range interaction forces which could be tuned
to achieve (amongst other things) powder-like behavior. Later Luciani et al.[1995] developed a similar particle system model specifically for granular materials using damped nonlinear springs at acoarse length scale and a faster linear model at a fine length scale.
Li and Moshell[1993] presented a dynamic height-field simulationof soil based on the standard engineering Mohr-Coulomb constitutive model of granular materials. Sumner et al.[1998] took a similarheightfield approach with simple displacement and erosion rules tomodel footprints, tracks, etc. Onoue and Nishita[2003] recently
extended this to multi-valued heightfields, allowing for some 3D
effects.
As we mentioned above, one approach to granular materials
is directly simulating the grains as interacting rigid bodies.
Milenkovic[1996] demonstrated a simulation of 1000 rigid spheres
falling through an hour-glass using optimization methods for resolving contact. Milenkovic and Schmidl[2001] improved the capability of this algorithm, and Guendelman et al.[2003] further accelerated rigid body simulation, showing 500 nonconvex rigid bodies
falling through a funnel.
Particles have been a popular technique in graphics for water simulation. Desbrun and Cani[1996] introduced Smooth Particle Hydrodynamics (see Monaghan[1992] for the standard review of SPH) to
the animation literature, demonstrating a viscous fluid flow simulation using particles alone. Müller et al.[2003] developed SPH further for water simulation, achieving impressive interactive results.
Premoze et al.[2003] used a variation on SPH with an approximate
projection solve for their water simulations, but generated the final
rendering geometry with a grid-based level set. Particles were also
used by Takeshita et al.[2003] for fire.
There has been more work on animating water with grids, however. Foster and Metaxas[1996] developed the first grid-based fully
3D water simulation in graphics. Stam[1999] provided the semiLagrangian advection method for faster simulation, but whose excessive numerical dissipation was later mitigated by Fedkiw etal.[2001] with higher order interpolation and vorticity confinement.
Foster and Fedkiw[2001] incorporated this into a water solver, andadded a level set—augmented by marker particles to counter mass
loss—for high quality surface tracking. To this model Génevaux
et al.[2003] added two-way fluid-solid interaction forces. Carlson
et al.[2002] added variable viscosity to a water solver, providing
beautiful simulations of wet sand dripping (achieved simply by increased viscosity). Enright et al.[2002b] extended the Foster andFedkiw approach with particles on both sides of the interface andvelocity extrapolation into the air. Losasso et al.[2004] extended
this to dynamically adapted octree grids, providing much enhanced
resolution, and Rasmussen et al.[2004] improved the boundary conditions and velocity extrapolation while adding several other features important for visual effects production. Carlson et al.[2004]
added better coupling between fluid and rigid body simulations—which we parenthetically note had its origins in scientific work on
simulating flow with sand grains. Hong and Kim[2003] and Takahashi et al.[2003] introduced volume-of-fluid algorithms for animating multi-phase flow (water and air, as opposed to just the waterof the free-surface flows animated above).
For more general plastic flow, most of the graphics work has dealt
with meshes. Terzopoulos and Fleischer[1988] first introduced
plasticity to physics-based animation, with a simple 1D model applied to their springs. O’Brien et al.[2002] added plasticity to a
fracture-capable tetrahedral-mesh finite element simulation to support ductile fracture. Irving et al.[2004] introduced a more sophisticated plasticity model along with their invertible finite elements,

also for tetrahedral meshes. Real-time elastic simulation including plasticity has been the focus of several papers (e.g. [Müller and
Gross 2004]). Closest in spirit to the current paper, Goktekin etal.[2004] added elastic forces and associated plastic flow to a water solver, enabling animation of a wide variety of non-Newtonian
fluids.
For a scientific description of the physics of granular materials,
we refer the reader to Jaeger et al.[1996]. In the soil mechanics
literature there is a long history of numerically simulating granular materials using a elasto-plastic finite element formulation, with
Mohr-Coulomb or Drucker-Prager yield conditions and various
non-associated plastic flow rules. We highlight the standard workof Nayak and Zienkiewicz[1972], and the book by Smith and Griffiths[1998] which contains code and detailed comments on FEM
simulation of granular materials. Landslides and avalanches are two
granular phenomena of particular interest, and generally have beenstudied using depth-averaged 2D equations introduced by Savage
and Hutter[1989]. For alternatives that simulate the material at thelevel of individual grains, Herrmann and Luding’s review[1998] isa good place to start.
Our new fluid code traces its history back to the early particle-incell (PIC) work of Harlow[1963] for compressible flow. Shortly
thereafter Harlow and Welch[1965] invented the marker-and-cell
method for incompressible flow, with its characteristic staggered
grid discretization and marker particles for tracking a free surface—
the other fundamental elements of our algorithm. PIC suffered from excessive numerical dissipation, which was cured by
the Fluid-Implicit-Particle (FLIP) method[Brackbill and Ruppel
1986]; Kothe and Brackbill[1992] later worked on adapting FLIP
to incompressible flow. Compressible FLIP was also extended
to a elasto-plastic finite element formulation, the Material Point
Method[Sulsky et al. 1995], which has been used to model granular materials at the level of individual grains[Bardenhagen et al.
2000] amongst other things. MPM was used by Konagai and Johansson[2001] for simulating large-scale (continuum) granular materials, though only in 2D and at considerable computational expense.

3

Sand Modeling

3.1

Frictional Plasticity

The standard approach to defining the continuum behavior of sandand other granular (or “frictional”) materials is in terms of plastic
yielding. Suppose the stress tensor σ has mean stress σm √= tr(σ )/3and Von Mises equivalent shear stress σ̄ = |σ − σm δ |F / 2 (where| · |F is the Frobenius norm). The Mohr-Coulomb1 law says thematerial will not yield as long as√

3σ̄ < sin φ σm(1)

where φ is the friction angle. Heuristically, this is just the classic Coulomb static friction law generalized to tensors: the shear
stress σ − σm δ that tries to force particles to slide against one another is annulled by friction if the pressure σm forcing the particles
against each other is proportionally strong enough. The coefficientof friction µ commonly used in Coulomb friction is replaced hereby sin φ .
1 Technically Mohr-Coulomb uses the Tresca definition of equivalent
shear stress, but since it is not smooth and poses numerical difficulties, mostpeople use Von Mises.

If the yield condition is met and the sand can start to flow, a flow
rule is required. The simplest reasonable flow direction is σ − σm δ .
Heuristically this means the sand is allowed to slide against itself
in the direction that the shearing force is pushing. Note that thisis nonassociated, i.e. not equal to the gradient of the yield condition. While associated flow rules are simpler and are usually assumed for plastic materials (and have been used exclusively before in graphics[O’Brien et al. 2002; Goktekin et al. 2004; Irving
et al. 2004]), they are impossible here. The Mohr-Coulomb law
compares the shear part of the stress to the dilatational part, unlike
normal plastic materials which ignore the dilation part. Thus an associated flow rule would allow the material to shear if and only if
it dilates proportionally—the more the sand flows, the more it expands. Technically sand does dilate a little when it begins to flow—
the particles need some elbow room to move past each other—but
this dilation stops as soon as the sand is freely flowing. This is afairly subtle effect that is very difficult to observe with the eye, so
we confidently ignore it for animation.
Once plasticity has been defined, the standard soil mechanics approach to simulation would be to add a simple elasticity relation,and use a finite element solver to simulate the sand. However,
we would prefer to avoid the additional complexity and expense of
FEM (which would require numerical quadrature, an implicit solver
to avoid stiff time step restrictions due to elasticity, return mapping
for plasticity, etc.). For the purposes of plausible animation we will
make some further simplifying assumptions so we can retrofit sandmodeling into an existing water solver.

3.2

A Simplified Model

We first ignore the nearly imperceptible elastic deformation regime(due to rock grains deforming) and the tiny volume changes thatsand undergoes as it starts and stops flowing. Thus our domain
can be decomposed into regions which are rigidly moving (where
the Mohr-Coulomb yield condition has not been met) and the rest
where we have an incompressible shearing flow.
We then further assume that the pressure required to make the entire velocity field incompressible will be similar to the true pressurein the sand. This is certainly false: for example, it is well known
that a column of sand in a silo has a maximum pressure independent of the height of the column2 whereas in a column of water
the pressure increases linearly with height. However, we can only
think of one case where this might be implausible in animation: thepressure limit means sand in an hour glass flows at a constant rate
independent of the amount of sand remaining (which is the reason
hour glasses work). We suggest our inaccurate results will still beplausible in most other cases.
Since we are neglecting elastic effects, we assume that where thesand is flowing the frictional stress is simplyD
σ f = − sin φ p p
1/3|D|F(2)

where D = (∇u + ∇uT )/2 is the strain rate. That is, the frictional
stress is the point on the yield surface that most directly resists thesliding.

Consider one grid cell, of side length ∆x. The relative sliding velocities between opposite faces in the cell are Vr = ∆xD. The massof the cell is M = ρ ∆x3 . Ignoring other cells, the forces required
to stop all sliding motion in a time step ∆t are F = −Vr M/∆t =−(∆xD)(ρ ∆x3 )/∆t, derived from integrating Newton’s law to get
the update formula V new = V + ∆tF/M. Dividing F by the crosssectional area ∆x2 to get the stress gives:

σrigid = −ρ D∆x2∆t(3)

We can check if this satisfies our yield inequality: if it does (i.e. if
the sand can resist forces and inertia trying to make it slide), thenthe material in that cell should be rigid at the end of the time step.
This gives us the following algorithm. Every time step, after advection, we do the usual water solver steps of adding gravity, applying
boundary conditions, solving for pressure to enforce incompressibility, and subtracting ∆t/ρ ∇p from the intermediate velocity field
to make it incompressible. Then we evaluate the strain rate tensor
D in each grid cell with standard central differences. If the resulting σrigid from equation 3 satisfies the yield condition (using thepressure we computed for incompressibility) we mark the grid cellas rigid and store σrigid at that cell. Otherwise, we store the sliding frictional stress σ f from equation 2. We then find all connectedgroups of rigid cells, and project the velocity field in each separate
group to the space of rigid body motions (i.e. find a translationaland rotational velocity which preserve the total linear and angularmomentum of the group). The remaining velocity values we update with the frictional stress, using standard central differences for
u+ = ∆t/ρ ∇ · σ f .

3.3

Frictional Boundary Conditions

So far we have covered the friction within the sand, but there isalso the friction between the sand and other objects (e.g. the solid
wall boundaries) to worry about. Our simple solution is to use thefriction formula of Bridson et al.[2002] when we enforce the u ·n ≥ 0 boundary condition on the grid. When we project out thenormal component of velocity, we reduce the tangential componentof velocity proportionately, clamping it to zero for the case of static
friction:µ |u · n|uT = max 0, 1 −
uT
(4)
|uT |

where µ is the Coulomb friction coefficient between the sand andthe wall.

This is a crucial addition to a fluid solver, since the boundary conditions previously discussed in the animation literature either don’t
permit any kind of sliding (u = 0) which means sand sticks even tovertical surfaces, or always allow sliding (u · n ≥ 0, possibly witha viscous reduction of tangential velocity) which means sand piles
will never be stable. Compare figure 2, which includes friction, tofigure 1 which doesn’t.

3.4

Cohesion

We finally need a way of determining the yield condition (whenthe sand should be rigid) without tracking elastic stress and strain.

A common enhancement to the basic Mohr-Coulomb condition isthe addition of cohesion:√(5)
3σ̄ < sin φ σm + c

2 This is due to the force that resists the weight of the sand above beingtransferred to the walls of the silo by friction—sometimes causing silos to
unexpectedly collapse as a result.

where c > 0 is the cohesion coefficient. This is appropriate forsoils or other slightly sticky materials that can support some tension
before yielding.

Figure 2: The sand bunny with a boundary friction coefficient µ = 0.6: compare to figure 1 where zero boundary friction was used.

We have found that including a very small amount of cohesion improves our results for supposedly cohesion-less materials like sand.
In what should be a stable sand pile, small numerical errors can allow some slippage; a tiny amount of cohesion counters this error
and makes the sand stable again, without visibly effecting the flow
when the sand should be moving.

to represent the water surface, quickly rounding off any small features. The most attractive alternative to level sets is volume-of-fluid(VOF), which conserves water up to floating-point precision: however it has difficulties maintaining an accurate but sharply-definedsurface. Coupling level sets with VOF is possible[Sussman 2003]
but at a significant implementation and computational cost.

However, for modeling soils with their true cohesion coefficient,
our method is too stable: obviously unstable structures are implausibly held rigid. We will investigate this issue in the future.

On the other hand, particles trivially handle advection with excellent accuracy—simply by letting them flow through the velocity
field using ODE solvers—but have difficulties with pressure andthe incompressibility condition, often necessitating smaller than desired time steps for stability. SPH methods in particular can’t tolerate nonuniform particle spacing, which can be difficult to enforce
throughout the entire length of a simulation.

4

Fluid Simulation Revisited

4.1

Grids and Particles

There are currently two main approaches to simulating fully threedimensional water in computer graphics: Eulerian grids and Lagrangian particles. Grid-based methods (recent papers include[Carlson et al. 2004; Goktekin et al. 2004; Losasso et al. 2004])
store velocity, pressure, some sort of indicator of where the fluid
is or isn’t, and any additional variables on a fixed grid. Usually astaggered “MAC” grid is used, allowing simple and stable pressure
solves. Particle-based methods, exemplified by Smoothed Particle
Hydrodynamics (see [Monaghan 1992; Müller et al. 2003; Premoze
et al. 2003] for example), dispense with grids except perhaps for accelerating geometric searches. Instead fluid variables are stored onparticles which represent actual chunks of fluid, and the motion of
the fluid is achieved simply by moving the particles themselves.
The Lagrangian form of the Navier-Stokes equations (sometimes
using a compressible formulation with a fictitious equation of state)is used to calculate the interaction forces on the particles.
The primary strength of the grid-based methods is the simplicity ofthe discretization and solution of the incompressibility condition,
which forms the core of the fluid solver. Unfortunately, grid-based
methods have difficulties with the advection part of the equations.
The currently favored approach in graphics, semi-Lagrangian advection, suffers from excessive numerical dissipation due to accumulated interpolation errors. To counter this nonphysical effect,a nonphysical term such as vorticity confinement must be added(at the expense of conserving angular or even linear momentum).
While high resolution methods from scientific computing can accurately solve for the advection of a well-resolved velocity field
through the fluid, their implementation isn’t entirely straightforward: their wide stencils make life difficult on coarse animation
grids that routinely represent features only one or two grid cells
thick. Moreover, even fifth-order accurate methods fail to accurately advect level sets[Enright et al. 2002a] as are commonly used

Realizing the strengths and weaknesses of the two approaches complement each other Foster and Fedkiw[2001] and later Enright etal.[2002b] augmented a grid-based method with marker particles
to correct the errors in the grid-based level set advection. Thisparticle-level set method, most recently implemented on an adaptive octree[Losasso et al. 2004], has produced the highest fidelity
water animations in the field. However, we believe particles can be
exploited even further, simplifying and accelerating the implementation, and affording some new benefits in flexibility.

4.2

Particle-in-Cell Methods

Continuing the graphics tradition of recalling early computational
fluid dynamics research, we return to the Particle-in-Cell (PIC)
method of Harlow[1963]. This was an early approach to simulating
compressible flow that handled advection with particles, but everything else on a grid. At each time step, the fluid variables at a grid
point were initialized as a weighted average of nearby particle values, then updated on the grid with the non-advection part of theequations. The new particle values were then interpolated from the
updated grid values, and finally the particles moved according to
the grid velocity field.
The major problem with PIC was the excessive numerical diffusion caused by repeatedly averaging and interpolating the fluid
variables. Brackbill and Ruppel[1986] cured this with the FluidImplicit-Particle (FLIP) method, which achieved “an almost total
absence of numerical dissipation and the ability to represent large
variations in data.” The crucial change was to make the particles thefundamental representation of the fluid, and use the auxiliary grid
simply to increment the particle variables according to the change
computed on the grid.

We have adapted PIC and FLIP to incompressible flow3 as follows:• Initialize particle positions and velocities• For each time step:• At each staggered MAC grid node, compute a weighted average 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 particle.
• For PIC: Interpolate the new grid velocity to the particles.
• Move particles through the grid velocity field with an ODE
solver, making sure to push them outside of solid wallboundaries.
• Write the particle positions to output.

Observe there is no longer any need to implement grid-based advection, or the matching schemes such as vorticity confinement andparticle-level set to counter numerical dissipation.
4.2.1

Initializing Particles

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 flow is underresolved at thesimulation resolution. With less particles per grid cell, we tend to
run into occasional “gaps” in the flow, and more particles allow for
too much noise. To help with surface reconstruction later (section5) we reposition particles that lie near the initial water surface (saywithin one grid cell width) to be exactly half a grid cell away fromthe surface.
4.2.2

Transferring to the Grid

Each grid point takes a weighted average of the nearby particles.
We define “near” as being contained in the cube of twice the gridcell width centered on the grid point. (Note that on a staggered
MAC grid, the different components of velocity will have offset
neighborhoods.) The weight of a particle in this cube is the standard trilinear weighting. We also mark the grid cells that contain atleast one particle. In future work we will investigate using a moreaccurate indicator, e.g. reconstructing a signed distance field on thegrid from distance values stored on the particles. This would allow
us to easily implement a second-order accurate free surface boundary condition[Enright et al. 2003] which significantly reduces grid
artifacts.
We also note that the grid in our simulation is purely an auxiliary
structure to the fundamental particle representation. In particular,
we are not required to use the same grid every time step. An obvious optimization which we have not yet implemented—rendering
is currently our bottleneck—is to define the grid according to thebounding box of the particles every time step. This also means wehave no a priori limits on the simulation domain. We plan in the
future to detect when clusters of particles have separated and use
separate bounding grids for them, for significantly improved efficiency.
3 The one precedent for this that we could find is a reference to an unpublished manuscript[Kothe and Brackbill 1992].

Figure 3: FLIP vs. PIC velocity update for the same simulation.
Notice the small-scale velocities preserved by FLIP but smoothedaway by PIC.

4.2.3

Solving on the Grid

We first add the acceleration due to gravity to the grid velocities.
We then construct a distance field φ (x) in the unmarked (non-fluid)part of the grid using fast sweeping[Zhao 2005] and use that to extend the velocity field outside the fluid 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 field again using fast sweeping (as we found it tobe marginally faster and simpler to implement than fast marching[Adalsteinsson and Sethian 1999]).
4.2.4

Updating Particle Velocities

At each particle, we trilinearly interpolate either the velocity (PIC)or the change in velocity (FLIP) recorded at the surrounding eightgrid points. For viscous flows, such as sand, the additional numerical viscosity of PIC is beneficial; for inviscid flows, such as water,
FLIP is preferable. See figure 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, aviscosity step on the grid in concert with PIC would be required forhighly viscous fluids).
4.2.5

Moving Particles

Once we have a velocity field defined on the grid (extrapolated toexist everywhere) we can move the particles around in it. We use a
simple RK2 ODE solver with five substeps each of which is limitedby the CFL condition (so that a particle moves less than one grid cellin each substep), similar to Enright et al.[2002b]. We additionallydetect when particles penetrate solid wall boundaries due simply totruncation error in RK2 and the grid-based velocity field, and movethem in the normal direction back to just outside the solid, to avoid
the occasional stuck-particle artifact this would otherwise cause.

5

Surface Reconstruction from Particles

Our simulations output the positions of the particles that define thelocation of the fluid. For high quality rendering we need to reconstruct a surface that wraps around the particles. Of course we cangive up on direct reconstruction, e.g. running a level set simulation
guided by the particles[Premoze et al. 2003]. While this is an effective solution, we believe a fully particle-based reconstruction canhave significant advantages.

frame; in the absence of a fast method of computing these to the desired precision, we currently fix all the particle radii to the constantaverage particle spacing and simply adjust our initial partial positions so that the surface particles are exactly this distance from thesurface. A small amount of additional grid smoothing, restricted todecreasing φ to avoid destroying features, reduces bump artifacts atlater frames.

Figure 4: Comparison of our implicit function and blobbies formatching a perfect circle defined by the points inside. The blobby
is the outer curve, ours is the inner curve.

Naturally the first approach that comes to mind is blobbies[Blinn1982]. For irregularly shaped blobs containing only a few particles,
this works beautifully. Unfortunately, it seems unable to match asurface such as a flat plane, a cone, or a large sphere from a large
quantity of irregularly spaced particles—as we might well see in atleast the initial conditions of a simulation. Bumps relating to theparticle distribution are always visible. A small improvement tothis was suggested in [Müller et al. 2003], where the contribution
from a particle was divided by the SPH estimate of density. This
overcomes some scaling issues but does not significantly reducethe bumps on what should be flat surfaces.
We thus take a different approach, guided by the principle that weshould exactly reconstruct the signed distance field of an isolatedparticle: for a single particle x0 with radius r0 , our implicit function
must be:φ (x) = |x − x0 | − r0
(6)
To generalize this, we use the same formula with x0 replaced by aweighted average of the nearby particle positions and r0 replaced aweighted average of their radii:

φ (x)

=

x̄

=

|x − x̄| − r̄(7)∑ wi xi(8)∑ wi ri(9)k(|x − xi |/R)
∑ j k(|x − x j |/R)

(10)

i

r̄

=

i

wi

=

On the other hand, we do enjoy significant advantages over reconstruction methods that are tied to level set simulations. The cost perframe is quite low—even in our unoptimized version which calculates and writes out a full 2503 grid, we average 40–50 seconds aframe on a 2GHz G5, with the bulk of the time spent on I/O. Moreover, every frame is independent, so this is easily farmed out tomultiple CPU’s and machines running in parallel.
If we can eliminate the need for grid-based smoothing in the future,
perhaps adapting an MLS approach such as in Shen et al.[2004],
we could do the surface reconstruction on the fly during rendering.
Apart from speed, the biggest advantage this would bring would
be accurate motion blur for fluids. The current technique for generating intermediate surfaces between frames (for a Monte Carlomotion blur solution) is to simply interpolate between level set gridvalues[Enright et al. 2002b]. However, this approach destroys smallfeatures that move further than their width in one frame: the interpolated values may all be positive. Unfortunately it’s exactly thesesmall and fast-moving features that most demand motion blur. With
particle-based surface reconstruction, the positions of the particles
can be interpolated instead, so that features are preserved at intermediate times.

6

Examples

We used pbrt[Pharr and Humphreys 2004] for rendering. For the
textured sand shading, we blended together a volumetric texture
advected around by the simulation particles, similar to the approachof Lamorlette[2001] for fireballs.
The bunny example in figures 1 and 2 were simulated with 269,322
particles on a 1003 grid, taking approximately 6 seconds per frameon a 2Ghz G5 workstation. We believe optimizing the particle
transfer code and using BLAS routines for the linear solve would
substantially improve this performance. Unoptimized smoothing
and text I/O cause the surface reconstruction on a 2503 grid to take
40–50 seconds per frame.

where k is a suitable kernel function that smoothly drops to zero—we used k(s) = max(0, (1 − s2 )3 ) since it avoids the need for squareroots and is reasonably well-shaped—and where R is the radius of
the neighborhood we consider around x. Typically we choose R to
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 flat orsmooth surfaces. See figure 4 for a comparison with blobbies using
the same kernel function.

The column test in figures 5 and 6 was simulated with 433,479 particles on a 100 × 60 × 60 grid, taking approximately 12 seconds perframe. Our cost is essentially linearly proportional to the numberof particles (or equivalently, occupied grid cells).

The most significant problem with this definition is artifacts in concave regions: spurious blobs of surface can appear, since x̄ may erroneously end up outside the surface in concavities. However, sincethese artifacts are significantly smaller than the particle radii we caneasily remove them without destroying true features by sampling
φ (x) on a higher resolution grid and then doing a simple smoothingpass.

We have presented a method for converting an existing fluid solver
into one capable of plausibly animating granular materials such as
sand. In addition, we developed a new fluid solver that combines
the strengths of both particles and grids, offering enhanced flexibility and efficiency. We offered a new method for reconstructing implicit surfaces from particles. Looking toward the future, weplan to more aggressively exploit the optimizations available with
PIC/FLIP (e.g. using multiple bounding grids), increase the accuracy of our boundary conditions, and implement motion blur of thereconstructed surface.

A secondary problem is that we require the radii to be accurate estimates of distance to the surface. This is nontrivial after the first7

Conclusion

Figure 5: A column of regular liquid is released.

B LINN , J. 1982. A generalization of algebraic surface drawing.
ACM Trans. Graph. 1, 3, 235–256.
B RACKBILL , J. U., AND RUPPEL , H. M. 1986. FLIP: a method
for adaptively zoned, particle-in-cell calculuations of fluid flows
in two dimensions. J. Comp. Phys. 65, 314–343.
B RIDSON , R., F EDKIW, R., AND A NDERSON , J. 2002. Robust
treatment of collisions, contact and friction for cloth animation.
ACM Trans. Graph. (Proc. SIGGRAPH) 21, 594–603.
C ARLSON , M., M UCHA , P., VAN H ORN III, R., AND T URK ,G.
2002. Melting and flowing. In Proc. ACM SIGGRAPH/Eurographics Symp. Comp. Anim., 167–174.
C ARLSON , M., M UCHA , P. J., AND T URK , G. 2004. Rigid fluid:
animating the interplay between rigid bodies and fluid. ACM
Trans. Graph. (Proc. SIGGRAPH) 23, 377–384.
D ESBRUN , M., AND C ANI , M.-P. 1996. Smoothed particles: A
new paradigm for animating highly deformable bodies. In Comput. Anim. and Sim. ’96 (Proc. of EG Workshop on Anim. andSim.), Springer-Verlag, R. Boulic and G. Hegron, Eds., 61–76.
Published under the name Marie-Paule Gascuel.
E NRIGHT, D., F EDKIW, R., F ERZIGER , J., AND M ITCHELL , I.
2002. A hybrid particle level set method for improved interface
capturing. J. Comp. Phys. 183, 83–116.

Figure 6: A column of granular material is released.

8

Acknowledgements

This work was in part supported by a grant from the Natural Sciences and Engineering Research Council of Canada. We would liketo thank Dr. Dawn Shuttle for her help with granular constitutive
models.

References

E NRIGHT, D., M ARSCHNER , S., AND F EDKIW, R. 2002. Animation and rendering of complex water surfaces. ACM Trans.
Graph. (Proc. SIGGRAPH) 21, 3, 736–744.
E NRIGHT, D., N GUYEN , D., G IBOU , F., AND F EDKIW, R. 2003.
Using the particle level set method and a second order accuratepressure boundary condition for free surface flows. In Proc. 4th
ASME-JSME Joint Fluids Eng. Conf., no. FEDSM2003–45144,
ASME.
F EDKIW, R., S TAM , J., AND J ENSEN , H. 2001. Visual simulationof smoke. In Proc. SIGGRAPH, 15–22.
F OSTER , N., AND F EDKIW, R. 2001. Practical animation of liquids. In Proc. SIGGRAPH, 23–30.
F OSTER , N., AND M ETAXAS , D. 1996. Realistic animation ofliquids. Graph. Models and Image Processing 58, 471–483.

A DALSTEINSSON , D., AND S ETHIAN , J. 1999. The fast construction of extension velocities in level set methods. J. Comput.
Phys. 148, 2–22.

G ÉNEVAUX , O., H ABIBI , A., AND D ISCHLER , J.-M. 2003. Simulating fluid-solid interaction. In Graphics Interface, 31–38.

BARDENHAGEN , S. G., B RACKBILL , J. U., AND S ULSKY, D.
2000. The material-point method for granular materials. Comput. Methods Appl. Mech. Engrg. 187, 529–541.

G OKTEKIN , T. G., BARGTEIL , A. W., AND O’B RIEN , J. F. 2004.
A method for animating viscoelastic fluids. ACM Trans. Graph.
(Proc. SIGGRAPH) 23, 463–468.

G UENDELMAN , E., B RIDSON , R., AND F EDKIW, R. 2003. Nonconvex rigid bodies with stacking. ACM Trans. Graph. (Proc.
SIGGRAPH) 22, 3, 871–878.

NAYAK , G. C., AND Z IENKIEWICZ , O. C. 1972. Elasto-plastic
stress analysis. A generalization for various constitutive relations
including strain softening. Int. J. Num. Meth. Eng. 5, 113–135.

H ARLOW, F., AND W ELCH , J. 1965. Numerical Calculationof Time-Dependent Viscous Incompressible Flow of Fluid withFree Surface. Phys. Fluids 8, 2182–2189.

O’B RIEN , J., BARGTEIL , A., AND H ODGINS , J. 2002. Graphicalmodeling of ductile fracture. ACM Trans. Graph. (Proc. SIGGRAPH) 21, 291–294.

H ARLOW, F. H. 1963. The particle-in-cell method for numericalsolution of problems in fluid dynamics. In Experimental arithmetic, high-speed computations and mathematics.

O NOUE , K., AND N ISHITA , T. 2003. Virtual sandbox. In PacificGraphics, 252–262.

H ERRMANN , H. J., AND L UDING , S. 1998. Modeling granular
media on the computer. Continuum Mech. Therm. 10, 189–231.
H ONG , J.-M., AND K IM , C.-H. 2003. Animation of bubbles inliquid. Comp. Graph. Forum (Eurographics Proc.) 22, 3, 253–
262.
I RVING , G., T ERAN , J., AND F EDKIW, R. 2004. Invertible finite
elements for robust simulation of large deformation. In Proc.
ACM SIGGRAPH/Eurographics Symp. Comp. Anim., 131–140.
JAEGER , H. M., NAGEL , S. R., AND B EHRINGER , R. P. 1996.
Granular solids, liquids, and gases. Rev. Mod. Phys. 68, 4, 1259–
1273.
KONAGAI , K., AND J OHANSSON , J. 2001. Two dimensional
Lagrangian particle finite-difference method for modeling large
soil deformations. Structural Eng./Earthquake Eng., JSCE 18,
2, 105s–110s.
KOTHE , D. B., AND B RACKBILL , J. U. 1992. FLIP-INC: aparticle-in-cell method for incompressible flows. Unpublishedmanuscript.
L AMORLETTE , A. 2001. Shrek effects—flames and dragon fireballs. SIGGRAPH Course Notes, Course 19, 55–66.
L I , X., AND M OSHELL , J. M. 1993. Modeling soil: Realtime
dynamic models for soil slippage and manipulation. In Proc.
SIGGRAPH, 361–368.
L OSASSO , F., G IBOU , F., AND F EDKIW, R. 2004. Simulating water and smoke with an octree data structure. ACM Trans. Graph.
(Proc. SIGGRAPH) 23, 457–462.

P HARR , M., AND H UMPHREYS , G. 2004. Physically Based Rendering: From Theory to Implementation. Morgan Kaufmann.
P REMOZE , S., TASDIZEN , T., B IGLER , J., L EFOHN , A., ANDW HITAKER , R. 2003. Particle–based simulation of fluids. In
Comp. Graph. Forum (Eurographics Proc.), vol. 22, 401–410.
R ASMUSSEN , N., E NRIGHT, D., N GUYEN , D., M ARINO , S.,
S UMNER , N., G EIGER , W., H OON , S., AND F EDKIW, R.
2004. Directible photorealistic liquids. In Proc. ACM SIGGRAPH/Eurographics Symp. Comp. Anim., 193–202.
S AVAGE , S. B., AND H UTTER , K. 1989. The motion of a finite
mass of granular material down a rough incline. J. Flui Mech.
199, 177–215.
S HEN , C., O’B RIEN , J. F., AND S HEWCHUK , J. R. 2004. Interpolating and approximating implicit surfaces from polygon soup.
ACM Trans. Graph. (Proc. SIGGRAPH) 23, 896–904.
S MITH , I. M., AND G RIFFITHS , D. V. 1998. Programming the
Finite Element Method. J. Wiley & Sons.
S TAM , J. 1999. Stable fluids. In Proc. SIGGRAPH, 121–128.
S ULSKY, D., Z HOU , S.-J., AND S CHREYER , H. L. 1995. Application of particle-in-cell method to solid mechanics. Comp.
Phys. Comm. 87, 236–252.
S UMNER , R. W., O’B RIEN , J. F., AND H ODGINS , J. K. 1998.
Animating sand, mud, and snow. In Graphics Interface, 125–
132.
S USSMAN , M. 2003. A second order coupled level set and volumeof-fluid method for computing growth and collapse of vapor bubbles. J. Comp. Phys. 187, 110–136.

L UCIANI , A., H ABIBI , A., AND M ANZOTTI , E. 1995. A multiscale physical model of granular materials. In Graphics Interface, 136–146.

TAKAHASHI , T., F UJII , H., K UNIMATSU , A., H IWADA , K.,
S AITO , T., TANAKA , K., AND U EKI , H. 2003. Realistic animation of fluid with splash and foam. Comp. Graph. Forum(Eurographics Proc.) 22, 3, 391–400.

M ILENKOVIC , V. J., AND S CHMIDL , H. 2001. Optimizationbased animation. In Proc. SIGGRAPH, 37–46.

TAKESHITA , D., OTA , S., TAMURA , M., F UJIMOTO , T., AND
C HIBA , N. 2003. Visual simulation of explosive flames. In
Pacific Graphics, 482–486.

M ILENKOVIC , V. J. 1996. Position-based physics: simulating the
motion of many highly interacting spheres and polyhedra. In
Proc. SIGGRAPH, 129–136.

T ERZOPOULOS , D., AND F LEISCHER , K. 1988. Modeling inelastic deformation: viscoelasticity, plasticity, fracture. Proc.
SIGGRAPH, 269–278.

M ILLER , G., AND P EARCE , A. 1989. Globular dynamics: a connected particle system for animating viscous fluids. In Comput.
& Graphics, vol. 13, 305–309.

Z HAO , H. 2005. A fast sweeping method for Eikonal equations.
Math. Comp. 74, 603–627.

M ONAGHAN , J. J. 1992. Smoothed particle hydrodynamics. Annu.
Rev. Astron. Astrophys. 30, 543–574.
M ÜLLER , M., AND G ROSS , M. 2004. Interactive virtual materials.
In Graphics Interface.
M ÜLLER , M., C HARYPAR , D., AND G ROSS , M. 2003. Particlebased fluid simulation for interactive applications. In Proc. ACM
SIGGRAPH/Eurographics Symp. Comp. Anim., 154–159.