Stable Fluids

Stable Fluids

Jos Stam
Alias wavefront
Abstract
articles have been published in various areas on how to compute
these equations numerically. Which solver to use in practice de-
Building animation tools for ﬂuid-like motions is an important and
pends largely on the problem at hand and on the computing power
challenging problem with many applications in computer graphics.
available. Most engineering tasks require that the simulation pro-
The use of physics-based models for ﬂuid ﬂow can greatly assist
vide accurate bounds on the physical quantities involved to answer
in creating such tools. Physical models, unlike key frame or pro-
questions related to safety, performance, etc. The visual appearance
cedural based techniques, permit an animator to almost effortlessly
(shape) of the ﬂow is of secondary importance in these applications.
create interesting, swirling ﬂuid-like behaviors. Also, the interac-
In computer graphics, on the other hand, the shape and the behav-
tion of ﬂows with objects and virtual forces is handled elegantly.
ior of the ﬂuid are of primary interest, while physical accuracy is
Until recently, it was believed that physical ﬂuid models were too
secondary or in some cases irrelevant. Fluid solvers, for computer
expensive to allow real-time interaction. This was largely due to the
graphics, should ideally provide a user with a tool that enables her
fact that previous models used unstable schemes to solve the phys-
to achieve ﬂuid-like effects in real-time. These factors are more im-
ical equations governing a ﬂuid. In this paper, for the ﬁrst time,
portant than strict physical accuracy, which would require too much
we propose an unconditionally stable model which still produces
computational power.
complex ﬂuid-like ﬂows. As well, our method is very easy to im-
In fact, most previous models in computer graphics were driven
plement. The stability of our model allows us to take larger time
by visual appearance and not by physical accuracy. Early ﬂow
steps and therefore achieve faster simulations. We have used our
models were built from simple primitives. Various combinations of
model in conjuction with advecting solid textures to create many
these primitives allowed the animation of particles systems [15, 17]
ﬂuid-like animations interactively in two- and three-dimensions.
or simple geometries such as leaves [23]. The complexity of the
ﬂows was greatly improved with the introduction of random tur-
CR Categories: I.3.7 [Computer Graphics]: Three-Dimensional
bulences [16, 20]. These turbulences are mass conserving and,
Graphics and Realism—Animation
therefore, automatically exhibit rotational motion. Also the tur-
Keywords: animation of ﬂuids, Navier-Stokes, stable solvers, im-
bulence is periodic in space and time, which is ideal for motion
plicit elliptic PDE solvers, interactive modeling, gaseous phenom-
“texture mapping” [19]. Flows built up from a superposition of
ena, advected textures
ﬂow primitives all have the disadvantage that they do not respond
dynamically to user-applied external forces. Dynamical models
of ﬂuids based on the Navier-Stokes equations were ﬁrst imple-
1
Introduction
mented in two-dimensions. Both Yaeger and Upson and Gamito
et al. used a vortex method coupled with a Poisson solver to cre-
One of the most intriguing problems in computer graphics is the
ate two-dimensional animations of ﬂuids [24, 8]. Later, Chen et
simulation of ﬂuid-like behavior. A good ﬂuid solver is of great
al. animated water surfaces from the pressure term given by a two-
importance in many different areas. In the special effects industry
dimensional simulation of the Navier-Stokes equations [2]. Their
there is a high demand to convincingly mimic the appearance and
method unlike ours is both limited to two-dimensions and is un-
behavior of ﬂuids such as smoke, water and ﬁre. Paint programs
stable. Kass and Miller linearize the shallow water equations to
can also beneﬁt from ﬂuid solvers to emulate traditional techniques
simulate liquids [12]. The simpliﬁcations do not, however, cap-
such as watercolor and oil paint. Texture synthesis is another pos-
ture the interesting rotational motions characteristic of ﬂuids. More
sible application. Indeed, many textures result from ﬂuid-like pro-
recently, Foster and Metaxas clearly show the advantages of us-
cesses, such as erosion. The modeling and simulation of ﬂuids is,
ing the full three-dimensional Navier-Stokes equations in creating
of course, also of prime importance in most scientiﬁc disciplines
ﬂuid-like animations [7]. Many effects which are hard to key frame
and in engineering. Fluid mechanics is used as the standard math-
manually such as swirling motion and ﬂows past objects are ob-
ematical framework on which these simulations are based. There
tained automatically. Their algorithm is based mainly on the work
is a consensus among scientists that the Navier-Stokes equations
of Harlow and Welch in computational ﬂuid dynamics, which dates
are a very good model for ﬂuid ﬂow. Thousands of books and
back to 1965 [11]. Since then many other techniques which Fos-
ter and Metaxas could have used have been developed. However,
¡
Alias wavefront, 1218 Third Ave, 8th Floor, Seattle, WA 98101, U.S.A.
their model has the advantage of being simple to code, since it is
jstam@aw.sgi.com
based on a ﬁnite differencing of the Navier-Stokes equations and
an explicit time solver. Similar solvers and their source code are
also available from the book of Griebel et al. [9]. The main prob-
lem with explicit solvers is that the numerical scheme can become
unstable for large time-steps. Instability leads to numerical sim-
ulations that “blow-up” and therefore have to be restarted with a
smaller time-step. The instability of these explicit algorithms sets
serious limits on speed and interactivity. Ideally, a user should be
able to interact in real-time with a ﬂuid solver without having to
worry about possible “blow ups”.
In this paper, for the ﬁrst time, we propose a stable algorithm
that solves the full Navier-Stokes equations. Our algorithm is very

easy to implement and allows a user to interact in real-time with
time
, then the evolution of these quantities over time is given
£!
three-dimensional ﬂuids on a graphics workstation. We achieve
by the Navier-Stokes equations [3]:
this by using time-steps much larger than the ones used by Fos-
"$#

ter and Metaxas. To obtain a stable solver we depart from Foster
(1)
£%
&
1
and Metaxas’ method of solution. Instead of their explicit Eulerian

#0"
"
"87

(2)
&
schemes, we use both Lagrangian and implicit methods to solve the
£(')¥

'
¡4365
3@9A©
2

Navier-Stokes equations. Our method cannot be found in the com-
where
is the kinematic viscosity of the ﬂuid, 2 is its density and
putational ﬂuids literature, since it is custom made for computer
5
is an external force. Some readers might be unfamiliar with this
graphics applications. The model would not be accurate enough
9
compact version of the Navier-Stokes equations. Eq. 2 is a vec-
for most engineering applications. Indeed, it suffers from too much
tor equation for the three (two in two-dimensions) components of
“numerical dissipation”, i.e., the ﬂow tends to dampen too rapidly
#
the velocity ﬁeld. The “ ” denotes a dot product between vec-
as compared to actual experiments. In a computer graphical appli-
"
tors, while the symbol
is the vector of spatial partial deriva-
cation, on the other hand, this is not so bad, especially in an interac-
&DC0&
&DC0&
"
tives. More precisely,
in two-dimensions and
tive system where the ﬂow is “kept alive” by an animator applying
£B¥
§©

&DC0&
&DC0&
&DC0&
"
in three-dimensions. We have also used
external forces. In fact, a ﬂow which does not dampen at all might
£E¥
§©
©

"
"E#F"
the shorthand notation
7
. The Navier-Stokes equations
be too chaotic and difﬁcult to control. As our results demonstrate
£
are obtained by imposing that the ﬂuid conserves both mass (Eq. 1)
we are able to produce nice swirling ﬂows despite the numerical
and momentum (Eq. 2). We refer the reader to any standard text
dissipation.
on ﬂuid mechanics for the actual derivation. These equations also
In this paper we apply our ﬂows mainly to the simulation of
have to be supplemented with boundary conditions. In this paper
gaseous-like phenomena. We employ our solver to update both
we will consider two types of boundary conditions which are use-
the ﬂow and the motion of densities within the ﬂow. To further
ful in practical applications: periodic boundary conditions and ﬁxed
increase the complexity of our animations we advect texture co-
boundary conditions. In the case of periodic boundaries the ﬂuid is
ordinates along with the density [13]. In this manner we are able
deﬁned on an
-dimensional torus (
). In this case there
to synthesize highly detailed “wispy” gaseous ﬂows even with low
G
G$£IH©QP
are no walls, just a ﬂuid which wraps around. Although such ﬂu-
resolution grids. We believe that the combination of physics-based
ids are not encountered in practice, they are very useful in creating
ﬂuid solvers and solid textures is the most promising method of
evolving texture maps. Also, these boundary conditions lead to a
achieving highly complex ﬂows in computer graphics.
very elegant implementation that uses the fast Fourier transform as
The next section presents the Navier-Stokes equations and the
shown below. The second type of boundary condition that we con-
derivation which leads to our method of solution. That section con-
sider is when the ﬂuid lies in some bounded domain
. In that case,
tains all the fundamental ideas and techniques needed to obtain a
R
the boundary conditions are given by a function TS
deﬁned on the
stable ﬂuids solver. Since it relies on sophisticated mathematical
&
boundary
of the domain. See Foster and Metaxas’ work for a
techniques, it is written in a mathematical physics jargon which
R
good discussion of these boundary conditions in the case of a hot
should be familiar to most computer graphics researchers working
ﬂuid [7]. In any case, the boundary conditions should be such that
in physics-based modeling. The application oriented reader who
the normal component of the velocity ﬁeld is zero at the boundary;
wishes only to implement our solver can skip Section 2 entirely and
no matter should traverse walls.
go straight to Section 3. There we describe our implementation of
The pressure and the velocity ﬁelds which appear in the Navier-
the solver, providing sufﬁcient information to code our technique.
Stokes equations are in fact related. A single equation for the ve-
Section 4 is devoted to several applications that demonstrate the
locity can be obtained by combining Eq. 1 and Eq. 2. We brieﬂy
power of our new solver. Finally, in Section 5 we conclude and
outline the steps that lead to that equation, since it is fundamen-
discuss future research. To keep this within the conﬁnes of a short
tal to our algorithm. We follow Chorin and Marsden’s treatment
paper, we have decided not to include a “tutorial-type” section on
of the subject (p. 36ff, [3]). A mathematical result, known as the
ﬂuid dynamics, since there are many excellent textbooks which pro-
Helmholtz-Hodge Decomposition, states that any vector ﬁeld
can
vide the necessary background and intuition. Readers who do not
U
uniquely be decomposed into the form:
have a background in ﬂuid dynamics and who wish to fully under-
stand the method in this paper should therefore consult such a text.
"4W

(3)
UV£
3
©
Mathematically inclined readers may wish to start with the excel-
"X#
W

lent book by Chorin and Marsden [3]. Readers with an engineering
where has zero divergence:
and
is a scalar ﬁeld. Any
£Y
bent on the other hand can consult the didactic book by Abbott [1].
vector ﬁeld is the sum of a mass conserving ﬁeld and a gradient
Also, Foster and Metaxas’ paper does a good job of introducing the
ﬁeld. This result allows us to deﬁne an operator
which projects
`
concepts from ﬂuid dynamics to the computer graphics community.
any vector ﬁeld
onto its divergence free part
. The
U
£a`bU
operator is in fact deﬁned implicitly by multiplying both sides of
"
Eq. 3 by “ ”:
2
Stable Navier-Stokes
"$#
"87QWdc
(4)
UV£
W
This is a Poisson equation for the scalar ﬁeld
with the Neumann
2.1
Basic Equations
&
boundary condition
on
. A solution to this equation is
e0f
£i
R
ehg
In this section we present the Navier-Stokes equations along with
used to compute the projection :
the manipulations that lead to our stable solver. A ﬂuid whose den-
"4Wdc

£!`bUV£pUq'
sity and temperature are nearly constant is described by a velocity
ﬁeld and a pressure ﬁeld . These quantities generally vary both
If we apply this projection operator on both sides of Eq. 2 we obtain
¡
in space and in time and depend on the boundaries surrounding the
a single equation for the velocity:
ﬂuid. We will denote the spatial coordinate by
, which for two-
&
¢

dimensional ﬂuids is
and three-dimensional ﬂuids is
#0"
"87

(5)
&
¢¤£¦¥¨§©
£!`$rs')¥

365
3t9huv©
equal to
. We have decided not to specialize our results

¥¨§©©
"
for a particular dimension. All results are thus valid for both two-
where we have used the fact that

and
. This is
`
£
`
¡w£i
dimensional and three-dimensional ﬂows unless stated otherwise.
our fundamental equation from which we will develop a stable ﬂuid
Given that the velocity and the pressure are known for some initial
solver.

q
C
UT
T

, where
is the spacing of their computational
£
)HPQ£SR
£SR
grid. Therefore, for small separations and/or large velocities, very
w2
w
small time steps have to be taken. On the other hand, we use a to-
1
w
3
tally different approach which results in an unconditionally stable
solver. No matter how big the time step is, our simulations will
¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡
u
never “blow up”. Our method is based on a technique to solve par-
¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡
¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡
.u=0
tial differential equations known as the method of characteristics.
¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡
w
Since this method is of crucial importance in obtaining our stable
¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡
w0
4
solver, we provide all the mathematical details in Appendix A. The
¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡ ¡ ¡ ¡ ¢ ¡
method, however, can be understood intuitively. At each time step
all the ﬂuid particles are moved by the velocity of the ﬂuid itself.
Figure 1: One simulation step of our solver is composed of steps.
Therefore, to obtain the velocity at a point
at the new time
,
¢
F3¡£)
The ﬁrst three steps may take the ﬁeld out of the space of divergent
we backtrace the point
through the velocity ﬁeld
over a time
¢
U
V)
free ﬁelds. The last projection step ensures that the ﬁeld is divergent
. This deﬁnes a path ¤
corresponding to a partial stream-
£
)
¥¨¢T©3W0
free after the entire simulation step.
line of the velocity ﬁeld. The new velocity at the point
is then
¢
set to the velocity that the particle, now at
, had at its previous
¢
location a time
ago:
£
)
c
¤
x
7
U
¥¨¢
T£pU0)0¥
¥¨¢T©
'
¥£)s
p(x,s)
Figure 2 illustrates the above. This method has several advantages.
Most importantly it is unconditionally stable. Indeed, from the
above equation we observe that the maximum value of the new
p(x,−∆t)
ﬁeld is never larger than the largest value of the previous ﬁeld.
Secondly, the method is very easy to implement. All that is re-
s
quired in practice is a particle tracer and a linear interpolator (see
0
−∆t
next Section). This method is therefore both stable and simple to
implement, two highly desirable properties of any computer graph-
ics ﬂuid solver. We employed a similar scheme to move densities
Figure 2: To solve for the advection part, we trace each point of
through user-deﬁned velocity ﬁelds [19]. Versions of the method of
the ﬁeld backward in time. The new velocity at
is therefore
characteristics were also used by other researchers. The application
¢
the velocity that the particle had a time
ago at the old location
was either employed in visualizing ﬂow ﬁelds [13, 18] or improv-
£
)
¤
.
ing the rendering of gas simulations [21, 5]. Our application of
¥¨¢T©
'
¥£)
the technique is fundamentally different, since we use it to update
the velocity ﬁeld, which previous researchers did not dynamically
2.2
Method of Solution
animate.
The third step solves for the effect of viscosity and is equivalent
Eq. 5 is solved from an initial state §¦

by marching
£
¥¨¢T©s A
to a diffusion equation:
through time with a time step
. Let us assume that the ﬁeld has
£
)
been resolved at a time and that we wish to compute the ﬁeld at a
&
7

"87
c
U
later time
. We resolve Eq. 5 over the time span
in four
7
&
£!5
U

3
¨£)
£
)

steps. We start from the solution
¦

of the previous
U
¥¨¢
v£
¥¨¢T©
s
time step and then sequentially resolve each term on the right hand
This is a standard equation for which many numerical procedures
side of Eq. 5, followed by a projection onto the divergent free ﬁelds.
have been developed. The most straightforward way of solving this
"
7
The general procedure is illustrated in Figure 1. The steps are:
equation is to discretize the diffusion operator
and then to do
an explicit time step as Foster and Metaxas did [7]. However, this
FE
54
D
©

©213 54
7
8@9A
method is unstable when the viscosity is large. We prefer, therefore,
!#"%$#&
!#"%$#&
!#"%$#&
!#"%$#&
c
¦
to use an implicit method:
7
U
¥¨¢

'
('
U
0)0¥¨¢

'
6'
U
¥¨¢

'
6'
U
CBA¥¨¢

'
('
U
CGF¥¨¢

The solution at time
is then given by the last velocity ﬁeld:
"87
3¨£)
7
't5(£)
U
CBF¥¨¢
T£pU
¥¨¢
Q©
r
YX
u

. A simulation is obtained by iterating these
¥¨¢T©
3H£)£pUCGF¥¨¢

steps. We now explain how each step is computed in more detail.
where
is the identity operator. When the diffusion operator is
X
The easiest term to solve is the addition of the external force .
discretized, this leads to a sparse linear system for the unknown
9
If we assume that the force does not vary considerably during the
ﬁeld
. Solving such a system can be done efﬁciently, however
U
CB
time step, then
(see below).
The fourth step involves the projection step, which makes the
¦
resulting ﬁeld divergence free. As pointed out in the previous sub-
U
0)
¥¨¢
T£pU
¥¨¢

3
I£)
9¥¨¢T©
s
section this involves the resolution of the Poisson problem deﬁned
is a good approximation of the effect of the force on the ﬁeld over
by Eq. 4:
the time step
. In an interactive system this is a good approxi-
"87QW
"$#
"4Wdc
£
)
mation, since forces are only applied at the beginning of each time
£
U
CB
U
CG
£pUCB
'
step.
The projection step, therefore, requires a good Poisson solver.
The next step accounts for the effect of advection (or convec-
Foster and Metaxas solved a similar equation using a relaxation
tion) of the ﬂuid on itself. A disturbance somewhere in the ﬂuid
scheme. Relaxation schemes, though, have poor convergence and
#"
propagates according to the expression

. This term
usually require many iterations. Foster and Metaxas reported that
')¥

makes the Navier-Stokes equations non-linear. Foster and Metaxas
they obtained good results even after a very small number of re-
resolved this component using ﬁnite differencing. Their method
laxation steps. However, since we are using a different method to
is stable only when the time step is sufﬁciently small such that
resolve for the advection step, we must use a more accurate method.

Indeed, the method of characteristics is more precise when the ﬁeld
is close to divergent free. More importantly from a visual point of
view, the projection step forces the ﬁelds to have vortices which re-
sult in more swirling-like motions. For these reasons we have used
a more accurate solver for the projection step.
The Poisson equation, when spatially discretized, becomes a
sparse linear system. Therefore, both the projection and the viscos-
ity steps involve the solution of a large sparse system of equations.
Multigrid methods, for example, can solve sparse linear systems in
linear time [10]. Since our advection solver is also linear in time,
Figure 3: The values of the discretized ﬁelds are deﬁned at the cen-
the complexity of our proposed algorithm is of complexity
.

8¥¢¡

ter of the grid cells.
Foster and Metaxas’ solver has the same complexity. This perfor-
mance is theoretically optimal since for a complicated ﬂuid, any
algorithm has to consult at least each cell of the computational grid.
ﬂuid and a texture coordinate. The evolution of this scalar ﬁeld is
conveniently described by an advection diffusion type equation:
2.3
Periodic Boundaries and the FFT
&
#h"
"87
8

When we consider a domain with periodic boundary conditions, our
&
£$'
8
b3@9BA
'
DCEA68b3GFHA
©

algorithm takes a particularly simple form. The periodicity allows
us to transform the velocity into the Fourier domain:
where
is a diffusion constant,
is a dissipation rate and
is
9BA
C"A
FIA
a source term. This equation is very similar in form to the Navier-
c
Stokes equation. Indeed, it includes an advection term, a diffusion

¥¨¢T©
s
'
6'¤£
¥
¢¥
©
s
term and a “force term”
. All these terms can be resolved exactly
FHA
"
In the Fourier domain the gradient operator “ ” is equivalent to
in the same way as the velocity of the ﬂuid. The dissipation term
the multiplication by
, where
1
. Consequently, both the
not present in the Navier-Stokes equation is solved as follows over
¦§¥
¦
v£©¨
'
diffusion step and the projection step are much simpler to solve.
a time-step:
Indeed the diffusion operator and the projection operators in the
c
1
Fourier domain are
¥
3
I£)CEAP8
¥¨¢T©
3I£)T£8
¥¨¢T©
s
Similar equations were used by Stam and Fiume to simulate ﬁre
"
7
1
7

'
5
(£)
'
6'
365(£)
X
and other gaseous phenomena [21]. However, their velocity ﬁelds
1
#
were not computed dynamically.
£
`bUB'6'
`
£
U
¥
¢¥
£
£
U
¥
¢¥'
¥
¢¥
U
£
¥
¢¥¥
©
7

We hope that the material in this section has convinced the reader
that our stable solver is indeed based on the full Navier-Stokes
T
T
where
. The operator £ projects the vector
onto the

£
¥
`
U
£
¥
¢¥
equations.
Also, we have pointed to the numerical techniques
plane which is normal to the wave number . The Fourier transform
¥
which should be used at each step of our solver. We now proceed
of the velocity of a divergent free ﬁeld is therefore always perpen-
to describe the implementation of our model in more detail.
dicular to its wavenumbers. The diffusion can be interpreted as a
low pass ﬁlter whose decay is proportional to both the time step
and the viscosity. These simple results demonstrate the power of
3
Our Solver
the Fourier transform. Indeed, we are able to completely transcribe
our solver in only a couple of lines. All that is required is a particle
3.1
Setup
tracer and a fast Fourier transform (FFT).
Our implementation handles both the motion of ﬂuids and the prop-
FourierStep( ¦ ,
,
):
agation by the ﬂuid of any number of substances like mass-density,
U
U
CG
£
)
add force:
¦
temperature or texture coordinates. Each quantity is deﬁned on ei-
U
0)v£pU
3
I£)
9
advect:
¤
ther a two-dimensional (NDIM=2) or three-dimensional (NDIM=3)
7
U
¥¨¢

£pU0)0¥
¥¨¢T©
'
¥£)
transform:
grid, depending on the application. The grid is deﬁned by its phys-
7
7
'&
U
£
£
! " "#%$
U
C
diffuse:
1
7
ical dimensions: origin O[NDIM] and length L[NDIM] of each
7
U
£
CBF¥¢¥T££
U
¥
¢¥
¥
365(£)

project:
side, and by its number of cells N[NDIM] in each coordinate. This
£
U
£
CG
£
`
(£
U
CB
in turn determines the size of each voxel D[i]=L[i]/N[i]. The
transform:
)
&
U
CG
£
! " "#0)
$1£
U
CG
deﬁnition of the grid is an input to our program which is speci-
ﬁed by the animator. The velocity ﬁeld is deﬁned at the center of
Since the Fourier transform is of complexity
, this

8¥¢¡32
4657¡

each cell as shown in Figure 3. Notice that previous researchers,
method is theoretically slightly more expensive than a method of
e.g., [7], deﬁned the velocity at the boundaries of the cells. We
solution relying on multi-grid solvers. However, this method is
prefer the cell-centered grid since it is more straightforward to im-
very easy to implement. We have used this algorithm to generate
plement. We allocate two grids for each component of the velocity:
the “liquid textures” of Section 4.
U0[NDIM] and U1[NDIM]. At each time step of our simulation
one grid corresponds to the solution obtained in the previous step.
2.4
Moving Substances through the Fluid
We store the new solution in the second grid. After each step, the
grids are swapped. We also allocate two grids to hold a scalar ﬁeld
A non-reactive substance which is injected into the ﬂuid will be ad-
corresponding to a substance transported by the ﬂow. Although our
vected by it while diffusing at the same time. Common examples of
implementation can handle any number of substances, for the sake
this phenomenon include the patterns created by milk stirred in cof-
of clarity we present only the algorithm for one ﬁeld in this section.
fee or the smoke rising from a cigarette. Let
be any scalar quantity
This scalar quantity is stored in the grids S0 and S1. The speed of
8
which is moved through the ﬂuid. Examples of this quantity include
interactivity is controlled by a single time step dt, which can be as
the density of dust, smoke or cloud droplets, the temperature of a
large as the animator wishes, since our algorithm is stable.

The physical properties of the ﬂuid are a function of its viscosity
for(i=0;i<NDIM;i++)
visc alone. By varying the viscosity, an animator can simulate a
addForce ( U0[i], F[i], dt );
wide range of substances ranging from glue-like matter to highly
for(i=0;i<NDIM;i++)
turbulent ﬂows. The properties of the substance are modeled by a
Transport ( U1[i], U0[i], U0, dt );
diffusion constant kS and a dissipation rate aS. Along with these
for(i=0;i<NDIM;i++)
parameters, the animator also must specify the values of these ﬁelds
Diffuse ( U0[i], U1[i], visc, dt );
on the boundary of the grid. There are basically two types: peri-
Project ( U1, U0, dt );
odic or ﬁxed. The boundary conditions can be of a different type
for each coordinate. When periodic boundary conditions are cho-
sen, the ﬂuid wraps around. This means that a piece of ﬂuid which
The general structure of the scalar ﬁeld solver is very similar to the
leaves the grid on one side reenters the grid on the opposite side.
above. It involves four steps: add the source, transport the ﬁeld by
In the case of ﬁxed boundaries, the value of each physical quantity
the velocity, diffuse and ﬁnally dissipate the ﬁeld. The scalar ﬁeld
must be speciﬁed at the boundary of the grid. The simplest method
solver shares some of the routines called by the velocity solver:
is to set the ﬁeld to zero at the boundary. We refer the reader to
Foster and Metaxas’ paper for an excellent description of different
boundary conditions and their resulting effects [7]. In the results
Sstep ( S1, S0, k, a, U, source, dt )
section we describe the boundary conditions chosen for each an-
addForce ( S0, source, dt );
imation. For the special case when the boundary conditions are
Transport ( S1, S0, U, dt );
periodic in each coordinate, a very elegant solver based on the fast
Diffuse ( S0, S1, k, dt );
Fourier transform can be employed. This algorithm is described in
Dissipate ( S1, S0, a, dt );
Section 2.3. We do not repeat it here since the solver in this section
is more general and can handle both types of boundary conditions.
The ﬂuid is set into motion by applying external forces to it.
The addForce routine adds the force ﬁeld multiplied by the time
We have written an animation system in which an animator with
step to each value of the ﬁeld. The dissipation routine Dissipate
a mouse can apply directional forces to the ﬂuid. The forces can
divides each element of the ﬁrst array by 1+dt*a and stores it
also be a function of other substances in the ﬂuid. For example,
in the new array. The Transport routine is a key step in our
a temperature ﬁeld moving through the ﬂuid can produce buoyant
simulation. It accounts for the movement of the substance due to
and turbulent forces. In our system we allow the user to create
the velocity ﬁeld. More importantly it is used to resolve the non-
all sorts of dependencies between the various ﬁelds, some of which
linearity of the Navier-Stokes equations. The general structure of
are described in the results section of this paper. We do not describe
this routine (in three-dimensions) is
our animation system in great detail since its functionality should
be evident from the examples of the next section. Instead we focus
Transport ( S1, S0, U, dt )
on our simulator, which takes the forces and parameters set by the
for each cell (i,j,k) do
animator as an input.
X = O+(i+0.5,j+0.5,k+0.5)*D;
TraceParticle ( X, U, -dt, X0 );
3.2
The Simulator
S1[i,j,k] = LinInterp ( X0, S0 );
end
Once we worked out the mathematics underlying the Navier-Stokes
equations in Section 2, our implementation became straightfor-
ward. We wish to emphasize that the theoretical developments of
The routine TraceParticle traces a path starting at X through
Section 2 are in no way gratuitous but are immensely useful in cod-
the ﬁeld U over a time -dt. The endpoint of this path is the new
ing compact solvers. In particular, casting the problem into a math-
point X0. We use both a simple second order Runge-Kutta (RK2)
ematical setting has allowed us to take advantage of the large body
method for the particle trace [14] and an adaptive particle tracer,
of work done in the numerical analysis of partial differential equa-
which subsamples the time step only in regions of high velocity gra-
tions. We have written the solver as a separate library of routines
dients, such as near object boundaries. The routine LinInterp
that are called by the interactive animation system. The entire li-
linearly interpolates the value of the scalar ﬁeld S at the location
brary consists of only roughly 500 lines of C code. The two main
X0. We note that we did not use a higher order interpolation, since
routines of this library update either the velocity ﬁeld Vstep or
this might lead to instabilities due to the oscillations and overshoots
a scalar ﬁeld Sstep over a given time step. We assume that the
inherent in such interpolants. On the other hand, higher order spline
external force is given by an array of vectors F[NDIM] and that
approximants may be used, though these tend to smooth out the re-
the source is given by an array Ssource for the scalar ﬁeld. The
sulting ﬂows.
general structure of our simulator looks like
To solve for the diffusion (Diffuse) and to perform the projec-
tion (Project) we need a sparse linear solver SolveLin. The
while ( simulating ) $
best theoretical choice is the multi-grid algorithm [10]. However,
/* handle display and user interaction */
we used a solver from the FISHPAK library since it was very easy to
/* get forces F and sources Ssource from the UI */
incorporate into our code and gave good results [22]1. In practice,
Swap(U1,U0); Swap(S1,S0);
it turned out to be faster than our implementation of the multi-grid
Vstep ( U1, U0, visc, F, dt );
algorithm. In Appendix B, we show exactly how these routines are
Sstep ( S1, S0, kS, aS, U1, Ssource, dt );
used to perform both the Diffuse step and the Project step.
&
These routines are ideal for domains with no internal boundaries.
The velocity solver is composed of four steps: the forces are added
When complex boundaries or objects are within the ﬂow, one can
to the ﬁeld, the ﬁeld is advected by itself, the ﬁeld diffuses due to
either use a sophisticated multi-grid solver or a good relaxation rou-
viscous friction within the ﬂuid, and in the ﬁnal step the velocity is
tine [9]. In any case, our simulator can easily accomodate new
forced to conserve mass. The general structure of this routine is:
solvers.
Vstep ( U1, U0, visc, F, dt )
1FISHPAK is available from http://www.netlib.org.

4
Results
Acknowledgments
Our Navier-Stokes solver can be used in many applications requir-
I would like to thank Marcus Grote for his informed input on ﬂuid
ing ﬂuid-like motions. We have implemented both the two- and the
dynamics and for pointing me to reference [3]. Thanks to Duncan
three-dimensional solvers in an interactive modeler that allows a
Brinsmead for his constructive criticisms all along. Thanks also to
user to interact with the ﬂuids in real-time. The motion is modeled
Pamela Jackson for carefully proofreading the paper and to Brad
by either adding density into the ﬂuid or by applying forces. The
Clarkson for his help with creating the videos.
evolution of the velocity and the density is then computed using our
solver. To further increase the visual complexity of the ﬂows, we
add textural detail to the density. By moving the texture coordinates
A
Method of Characteristics
using the scalar solver as well, we achieve highly detailed ﬂows. To
compensate for the high distortions that the texture maps undergo,
The method of characteristics can be used to solve advection equa-
we use three sets of texture coordinates which are periodically reset
tions of the type
to their initial (unperturbed) values. At every moment the resulting
&
texture map is the superposition of these three texture maps. This
#0"

8
¥¨¢T©
s
¦
&
£$'£¢
¥¨¢

8
¥¨¢T©
s
8
¥¨¢T© A£8
¥¨¢
Q©
idea was ﬁrst suggested by Max et al. [13].

Figure 4.(a) shows a sequence of frames from an animation
where
is a scalar ﬁeld,
is a steady vector ﬁeld and ¦ is the ﬁeld
where the user interacts with one of our liquid textures. The ﬁgure
8
¢
8
at time
. Let ¤
¦
denote the characteristics of the vector
on the backcover of the SIGGRAPH proceedings is another frame
v£Y
¥¨¢
©
s
ﬁeld
which ﬂow through the point
¦
at
:
of a similar sequence with a larger grid size (1
¢
¢
£!
7
).
F
Figures 4.(b) through 4.(g) show frames from various animations
¤
c

¤
¦
¤
¦
¤
¦
¦
that we generated using our three-dimensional solver. In each case
¤
¥¨¢
©
sT£¥¢
¥
¥¨¢
©
s
¥¨¢
©s A£p¢

the animations were created by allowing the animator to place den-
sity and apply forces in real-time. The gases are volume rendered
Now let
¦
¤
¦
be the value of the ﬁeld along the
8
¦
¥¨¢
©
sT£8
¥
¥¨¢
©
sQ©
s
using the three-dimensional hardware texture mapping capabilities
characteristic passing through the point
¦
at
. The variation
¢

£$
of our SGI Octane workstation. We also added a single pass that
of this quantity over time can be computed using the chain rule of
computes self-shadowing effects from a directional light source in
differentiation:
&
¤
a ﬁxed position. It should be evident that the quality of the render-
#0"
c
8
¦
8
ings could be further improved using more sophisticated rendering
&
¤
£
3
§¢
8
8£!

hardware or software. Our grid sizes ranged from 1¡ B to
B
with
PF
This shows that the value of the scalar does not vary along the
frame rates fast enough to monitor the animations while being able
streamlines. In particular, we have
¦
¦
¦
¦
.
to control their behavior. In most of these animations we added a
8
¦
¥¨¢
©
s
£
¨¦
8
¥¨¢
©s Av£
8
¥¨¢

Therefore, the initial ﬁeld and the characteristics entirely deﬁne the
“noise” term which is proportional to the amount of density (the
solution to the advection problem. The ﬁeld for a given time
and
factor of proportionality being a user deﬁned parameter). This pro-

location
is computed by ﬁrst tracing the location
back in time
duced nice billowing motions in some of our animations. In Fig-
¢
¢
along the characteristic to get the point
¦
, and then evaluating the
ures 4.(d)-(e) we used a fractal texture map, while in Figure 4.(g)
¢
initial ﬁeld at that point:
we used a texture map consisting of evenly spaced lines.
All of our animations were created on an SGI Octane worksta-
c
¤
¦
¦
¦
8
¥
¥¨¢
©sQ©
sT£8
¥¨¢

tion with a R10K processor and 192 Mbytes of memory.
We use this method to solve the advection equation over a time
interval
for the ﬂuid. In this case,

and ¦ is
©

©
3
£
)
¢
w£
¥¨¢T©
s
8
5
Conclusions
any of the components of the ﬂuid’s velocity at time .

The motivation of this paper was to create a general software sys-
tem that allows an animator to design ﬂuid-like motions in real time.
B
FISHPAK Routines
Our initial intention was to base our system on Foster and Metaxas’
work. However, the instabilities inherent in their method forced us
The linear solver POIS3D from FISHPAK is designed to solve a
to develop a new algorithm. Our solver has the property of being
general system of ﬁnite difference equations of the type:
unconditionally stable and it can handle a wide variety of ﬂuids in
both two- and three-dimensions. The results that accompany this
K1*(S[i-1,j,k]-2*S[i,j,k]+S[i+1,j,k]) +
paper clearly demonstrate that our solver is powerful enough to al-
K2*(S[i,j-1,k]-2*S[i,j,k]+S[i,j+1,k]) +
low an animator to achieve many ﬂuid-like effects. We therefore
A[k]*S[i,j,k-1]+B[k]*S[i,j,k]+
.
believe that our solver is a substantial improvement over previous
For the diffusion solver, the values of the constants on the left hand
work in this area. The work presented here does not, however, dis-
side are:
credit previous, more visually oriented models. In particular, we
believe that the combination of our ﬂuid solvers with solid textures,
K1 = -dt*k/(D[0]*D[0]),
for example, may be a promising area of future research [4]. Our
K2 = -dt*k/(D[1]*D[1]),
ﬂuid solvers can be used to generate the overall motion, while the
A[k] = C[k] = -dt*k/(D[2]*D[2]) and
solid texture can add additional detail for higher quality animations.
B[k] = 1+2*dt*k/(D[2]*D[2]),
Also we have not addressed the problem of simulating ﬂuids with
free boundaries, such as water [6]. This problem is considerably
while the right hand side is equal to the grid containing the previous
more difﬁcult, since the geometry of the boundary evolves dynam-
solution: F=S0. In the projection step these constants are equal to
ically over time. We hope, however, that our stable solvers may
be applied to this problem as well. Also, we wish to extend our
K1 = 1/(D[0]*D[0]), K2 = 1/(D[1]*D[1]),
solver to ﬁnite element boundary-ﬁtted meshes. We are currently
A[k] = C[k] = 1/(D[2]*D[2]) and
investigating such extensions.
B[k] = -2/(D[2]*D[2]),

while the right hand side is equal to the divergence of the velocity
[12] M. Kass and G. Miller. Rapid, Stable Fluid Dynamics for
ﬁeld:
Computer Graphics. ACM Computer Graphics (SIGGRAPH
’90), 24(4):49–57, August 1990.
F[i,j,k] = 0.5*((U[i+1,j,k]-U[i-1,j,k])/D[0]+
(U[i,j+1,k]-U[i,j-1,k])/D[1]+
[13] N. Max, R. Crawﬁs, and D. Williams. Visualizing Wind Ve-
(U[i,j,k+1]-U[i,j,k-1])/D[2]).
locities by Advecting Cloud Textures. In Proceedings of Vi-
sualization ’92, pages 179–183, Los Alamitos CA, October
The gradient of the solution is then subtracted from the previous
1992. IEEE CS Press.
solution:
[14] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetter-
U1[0][i,j,k] = U0[0][i,j,k] -
ling. Numerical Recipes in C. The Art of Scientiﬁc Computing.
0.5*(S[i+1,j,k]-S[i-1,j,k])/D[0],
Cambridge University Press, Cambridge, 1988.
U1[1][i,j,k] = U0[1][i,j,k] -
0.5*(S[i,j+1,k]-S[i,j-1,k])/D[1],
[15] W. T. Reeves. Particle Systems. A Technique for Modeling
U1[2][i,j,k] = U0[2][i,j,k] -
a Class of Fuzzy Objects. ACM Computer Graphics (SIG-
0.5*(S[i,j,k+1]-S[i,j,k-1])/D[2].
GRAPH ’83), 17(3):359–376, July 1983.
[16] M. Shinya and A. Fournier. Stochastic Motion - Motion Under
The FISHPAK routine is also able to handle different types of
the Inﬂuence of Wind. In Proceedings of Eurographics ‘92,
boundary conditions, both periodic and ﬁxed.
pages 119–128, September 1992.
References
[17] K. Sims. Particle Animation and Rendering Using Data Paral-
lel Computation. ACM Computer Graphics (SIGGRAPH ’90),
24(4):405–413, August 1990.
[1] M. B. Abbott. Computational Fluid Dynamics: An Introduc-
tion for Engineers. Wiley, New York, 1989.
[18] K. Sims. Choreographed Image Flow. The Journal Of Visual-
ization And Computer Animation, 3:31–43, 1992.
[2] J. X. Chen, N. da Vittoria Lobo, C. E. Hughes, and J. M.
Moshell. Real-Time Fluid Simulation in a Dynamic Virtual
[19] J. Stam. A General Animation Framework for Gaseous Phe-
Environment. IEEE Computer Graphics and Applications,
nomena.
ERCIM Research Report, R047, January 1997.
pages 52–61, May-June 1997.
http://www.ercim.org/publications/technical reports/047-abstract.html.
[3] A. J. Chorin and J. E. Marsden. A Mathematical Introduc-
[20] J. Stam and E. Fiume. Turbulent Wind Fields for Gaseous
tion to Fluid Mechanics. Springer-Verlag. Texts in Applied
Phenomena. In Proceedings of SIGGRAPH ’93, pages 369–
Mathematics 4. Second Edition., New York, 1990.
376. Addison-Wesley Publishing Company, August 1993.
[4] D. Ebert, K. Musgrave, D. Peachy, K. Perlin, and S. Worley.
[21] J. Stam and E. Fiume. Depicting Fire and Other Gaseous
Texturing and Modeling: A Procedural Approach. AP Profes-
Phenomena Using Diffusion Processes. In Proceedings of
sional, 1994.
SIGGRAPH ’95, pages 129–136. Addison-Wesley Publishing
Company, August 1995.
[5] D. S. Ebert, W. E. Carlson, and R. E. Parent. Solid Spaces
and Inverse Particle Systems for Controlling the Animation of
[22] P. N. Swarztrauber and R. A. Sweet. Efﬁcient Fortran Subpro-
Gases and Fluids. The Visual Computer, 10:471–483, 1994.
grams for the Solution of Separable Elliptic Partial Differen-
tial Equations. ACM Transactions on Mathematical Software,
[6] N. Foster and D. Metaxas.
Realistic Animation of Liq-
5(3):352–364, September 1979.
uids. Graphical Models and Image Processing, 58(5):471–
[23] J. Wejchert and D. Haumann.
Animation Aerodynamics.
483, 1996.
ACM Computer Graphics (SIGGRAPH ’91), 25(4):19–22,
[7] N. Foster and D. Metaxas. Modeling the Motion of a Hot,
July 1991.
Turbulent Gas. In Computer Graphics Proceedings, Annual
[24] L. Yaeger and C. Upson. Combining Physical and Visual Sim-
Conference Series, 1997, pages 181–188, August 1997.
ulation. Creation of the Planet Jupiter for the Film 2010. ACM
Computer Graphics (SIGGRAPH ’86), 20(4):85–93, August
[8] M. N. Gamito, P. F. Lopes, and M. R. Gomes.
Two-
1986.
dimensional Simulation of Gaseous Phenomena Using Vor-
tex Particles. In Proceedings of the 6th Eurographics Work-
shop on Computer Animation and Simulation, pages 3–15.
Springer-Verlag, 1995.
[9] M. Griebel, T. Dornseifer, and T. Neunhoeffer.
Numeri-
cal Simulation in Fluid Dynamics: A Practical Introduction.
SIAM, Philadelphia, 1998.
[10] W. Hackbusch.
Multi-grid Methods and Applications.
Springer Verlag, Berlin, 1985.
[11] F. H. Harlow and J. E. Welch.
Numerical Calculation of
Time-Dependent Viscous Incompressible Flow of Fluid with
Free Surface. The Physics of Fluids, 8:2182–2189, December
1965.

(a)
(b)
(c)
(d)
(e)
(b)
(f)
(g)
Figure 4: Snapshots from our interactive ﬂuid solver.