Abstract
Building animation tools for fluid-like motions is an important and challenging problem with many applications in computer graphics. The use of physics-based models for fluid flow can greatly assist in creating such tools. Physical models, unlike key frame or pro- cedural based techniques, permit an animator to almost effortlessly create interesting, swirling fluid-like behaviors. Also, the interac- tion of flows with objects and virtual forces is handled elegantly. Until recently, it was believed that physical fluid models were too expensive to allow real-time interaction. This was largely due to the fact that previous models used unstable schemes to solve the phys- ical equations governing a fluid. In this paper, for the first time, we propose an unconditionally stable model which still produces complex fluid-like flows. As well, our method is very easy to im- plement. The stability of our model allows us to take larger time steps and therefore achieve faster simulations. We have used our model in conjuction with advecting solid textures to create many fluid-like animations interactively in two- and three-dimensions.
CR Categories: I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism—Animation
Keywords: animation of fluids, Navier-Stokes, stable solvers, im- plicit elliptic PDE solvers, interactive modeling, gaseous phenom- ena, advected textures
1 Introduction
One of the most intriguing problems in computer graphics is the simulation of fluid-like behavior. A good fluid solver is of great importance in many different areas. In the special effects industry there is a high demand to convincingly mimic the appearance and behavior of fluids such as smoke, water and fire. Paint programs can also benefit from fluid solvers to emulate traditional techniques such as watercolor and oil paint. Texture synthesis is another pos- sible application. Indeed, many textures result from fluid-like pro- cesses, such as erosion. The modeling and simulation of fluids is, of course, also of prime importance in most scientific disciplines and in engineering. Fluid mechanics is used as the standard math- ematical framework on which these simulations are based. There is a consensus among scientists that the Navier-Stokes equations are a very good model for fluid flow. Thousands of books and
Alias wavefront, 1218 Third Ave, 8th Floor, Seattle, WA 98101, U.S.A.
jstam@aw.sgi.com
articles have been published in various areas on how to compute these equations numerically. Which solver to use in practice de- pends largely on the problem at hand and on the computing power available. Most engineering tasks require that the simulation pro- vide accurate bounds on the physical quantities involved to answer questions related to safety, performance, etc. The visual appearance (shape) of the flow is of secondary importance in these applications. In computer graphics, on the other hand, the shape and the behav- ior of the fluid are of primary interest, while physical accuracy is secondary or in some cases irrelevant. Fluid solvers, for computer graphics, should ideally provide a user with a tool that enables her to achieve fluid-like effects in real-time. These factors are more im- portant than strict physical accuracy, which would require too much computational power.
In fact, most previous models in computer graphics were driven by visual appearance and not by physical accuracy. Early flow models were built from simple primitives. Various combinations of these primitives allowed the animation of particles systems [15, 17] or simple geometries such as leaves [23]. The complexity of the flows was greatly improved with the introduction of random tur- bulences [16, 20]. These turbulences are mass conserving and, therefore, automatically exhibit rotational motion. Also the tur- bulence is periodic in space and time, which is ideal for motion “texture mapping” [19]. Flows built up from a superposition of flow primitives all have the disadvantage that they do not respond dynamically to user-applied external forces. Dynamical models of fluids based on the Navier-Stokes equations were first imple- mented in two-dimensions. Both Yaeger and Upson and Gamito et al. used a vortex method coupled with a Poisson solver to cre- ate two-dimensional animations of fluids [24, 8]. Later, Chen et al. animated water surfaces from the pressure term given by a two- dimensional simulation of the Navier-Stokes equations [2]. Their method unlike ours is both limited to two-dimensions and is un- stable. Kass and Miller linearize the shallow water equations to simulate liquids [12]. The simplifications do not, however, cap- ture the interesting rotational motions characteristic of fluids. More recently, Foster and Metaxas clearly show the advantages of us- ing the full three-dimensional Navier-Stokes equations in creating fluid-like animations [7]. Many effects which are hard to key frame manually such as swirling motion and flows past objects are ob- tained automatically. Their algorithm is based mainly on the work of Harlow and Welch in computational fluid dynamics, which dates back to 1965 [11]. Since then many other techniques which Fos- ter and Metaxas could have used have been developed. However, their model has the advantage of being simple to code, since it is based on a finite 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 fluid solver without having to worry about possible “blow ups”.
In this paper, for the first time, we propose a stable algorithm that solves the full Navier-Stokes equations. Our algorithm is very
Stable Fluids
Jos Stam Alias wavefront
easy to implement and allows a user to interact in real-time with three-dimensional fluids on a graphics workstation. We achieve 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 and Metaxas’ method of solution. Instead of their explicit Eulerian schemes, we use both Lagrangian and implicit methods to solve the Navier-Stokes equations. Our method cannot be found in the com- putational fluids literature, since it is custom made for computer graphics applications. The model would not be accurate enough for most engineering applications. Indeed, it suffers from too much “numerical dissipation”, i.e., the flow tends to dampen too rapidly as compared to actual experiments. In a computer graphical appli- cation, on the other hand, this is not so bad, especially in an interac- tive system where the flow is “kept alive” by an animator applying external forces. In fact, a flow which does not dampen at all might be too chaotic and difficult to control. As our results demonstrate we are able to produce nice swirling flows despite the numerical dissipation.
In this paper we apply our flows mainly to the simulation of gaseous-like phenomena. We employ our solver to update both the flow and the motion of densities within the flow. To further increase the complexity of our animations we advect texture co- ordinates along with the density [13]. In this manner we are able to synthesize highly detailed “wispy” gaseous flows even with low resolution grids. We believe that the combination of physics-based fluid solvers and solid textures is the most promising method of achieving highly complex flows in computer graphics.
The next section presents the Navier-Stokes equations and the derivation which leads to our method of solution. That section con- tains all the fundamental ideas and techniques needed to obtain a stable fluids solver. Since it relies on sophisticated mathematical techniques, it is written in a mathematical physics jargon which should be familiar to most computer graphics researchers working in physics-based modeling. The application oriented reader who wishes only to implement our solver can skip Section 2 entirely and go straight to Section 3. There we describe our implementation of the solver, providing sufficient information to code our technique. Section 4 is devoted to several applications that demonstrate the power of our new solver. Finally, in Section 5 we conclude and discuss future research. To keep this within the confines of a short paper, we have decided not to include a “tutorial-type” section on fluid dynamics, since there are many excellent textbooks which pro- vide the necessary background and intuition. Readers who do not have a background in fluid dynamics and who wish to fully under- stand the method in this paper should therefore consult such a text. Mathematically inclined readers may wish to start with the excel- lent book by Chorin and Marsden [3]. Readers with an engineering bent on the other hand can consult the didactic book by Abbott [1]. Also, Foster and Metaxas’ paper does a good job of introducing the concepts from fluid dynamics to the computer graphics community.
time , then the evolution of these quantities over time is given by the Navier-Stokes equations [3]:
(1) (2)
2 2.1
Stable Navier-Stokes Basic Equations
is the kinematic viscosity of the fluid,
is an external force. Some readers might be unfamiliar with this
compact version of the Navier-Stokes equations. Eq. 2 is a vec- tor equation for the three (two in two-dimensions) components of the velocity field. The “ ” denotes a dot product between vec-
the shorthand notation
where
is its density and
is the vector of spatial partial deriva- in two-dimensions and in three-dimensions. We have also used . The Navier-Stokes equations are obtained by imposing that the fluid conserves both mass (Eq. 1) and momentum (Eq. 2). We refer the reader to any standard text on fluid mechanics for the actual derivation. These equations also have to be supplemented with boundary conditions. In this paper we will consider two types of boundary conditions which are use- ful in practical applications: periodic boundary conditions and fixed boundary conditions. In the case of periodic boundaries the fluid is defined on an -dimensional torus ( ). In this case there are no walls, just a fluid which wraps around. Although such flu- ids are not encountered in practice, they are very useful in creating evolving texture maps. Also, these boundary conditions lead to a very elegant implementation that uses the fast Fourier transform as shown below. The second type of boundary condition that we con- sider is when the fluid lies in some bounded domain . In that case, the boundary conditions are given by a function defined on the boundary of the domain. See Foster and Metaxas’ work for a good discussion of these boundary conditions in the case of a hot fluid [7]. In any case, the boundary conditions should be such that the normal component of the velocity field is zero at the boundary;
tors, while the symbol tives. More precisely,
no matter should traverse walls.
The pressure and the velocity fields which appear in the Navier-
Stokes equations are in fact related. A single equation for the ve- locity can be obtained by combining Eq. 1 and Eq. 2. We briefly outline the steps that lead to that equation, since it is fundamen- tal to our algorithm. We follow Chorin and Marsden’s treatment of the subject (p. 36ff, [3]). A mathematical result, known as the Helmholtz-Hodge Decomposition, states that any vector field can uniquely be decomposed into the form:
(3)
where has zero divergence:
vector field is the sum of a mass conserving field and a gradient field. This result allows us to define an operator which projects any vector field onto its divergence free part . The operator is in fact defined implicitly by multiplying both sides of Eq.3by“ ”:
(4)
This is a Poisson equation for the scalar field with the Neumann boundary condition on . A solution to this equation is
and is a scalar field. Any
In this section we present the Navier-Stokes equations along with the manipulations that lead to our stable solver. A fluid whose den- sity and temperature are nearly constant is described by a velocity field and a pressure field . These quantities generally vary both in space and in time and depend on the boundaries surrounding the fluid. We will denote the spatial coordinate by , which for two- dimensional fluids is and three-dimensional fluids is equal to . We have decided not to specialize our results for a particular dimension. All results are thus valid for both two- dimensional and three-dimensional flows unless stated otherwise. Given that the velocity and the pressure are known for some initial
used to compute the projection
:
If we apply this projection operator on both sides of Eq. 2 we obtain a single equation for the velocity:
(5)
where we have used the fact that and . This is our fundamental equation from which we will develop a stable fluid solver.
7
EM < < LKD'J!%&#$!IH<
7><= <
7
456
4
?@9" G7F87>< 4 EBCADA
9?@93%&87
9 ;:
9" ! 87
2301/
/
.- %
+,
+, +, -
+, +, *
)'(!%&#$!"
#
'
u
q
w2 w w1 3
w w0 4
, where is the spacing of their computational grid. Therefore, for small separations and/or large velocities, very small time steps have to be taken. On the other hand, we use a to- tally different approach which results in an unconditionally stable solver. No matter how big the time step is, our simulations will never “blow up”. Our method is based on a technique to solve par- tial differential equations known as the method of characteristics. Since this method is of crucial importance in obtaining our stable solver, we provide all the mathematical details in Appendix A. The method, however, can be understood intuitively. At each time step all the fluid particles are moved by the velocity of the fluid itself. Therefore, to obtain the velocity at a point at the new time , we backtrace the point through the velocity field over a time
. This defines a path corresponding to a partial stream- line of the velocity field. 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:
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 field is never larger than the largest value of the previous field. Secondly, the method is very easy to implement. All that is re- quired in practice is a particle tracer and a linear interpolator (see next Section). This method is therefore both stable and simple to implement, two highly desirable properties of any computer graph- ics fluid solver. We employed a similar scheme to move densities through user-defined velocity fields [19]. Versions of the method of characteristics were also used by other researchers. The application was either employed in visualizing flow fields [13, 18] or improv- ing the rendering of gas simulations [21, 5]. Our application of the technique is fundamentally different, since we use it to update the velocity field, which previous researchers did not dynamically animate.
The third step solves for the effect of viscosity and is equivalent to a diffusion equation:
This is a standard equation for which many numerical procedures have been developed. The most straightforward way of solving this equation is to discretize the diffusion operator and then to do an explicit time step as Foster and Metaxas did [7]. However, this method is unstable when the viscosity is large. We prefer, therefore, to use an implicit method:
where is the identity operator. When the diffusion operator is discretized, this leads to a sparse linear system for the unknown field . Solving such a system can be done efficiently, however (see below).
The fourth step involves the projection step, which makes the resulting field divergence free. As pointed out in the previous sub- section this involves the resolution of the Poisson problem defined by Eq. 4:
The projection step, therefore, requires a good Poisson solver. Foster and Metaxas solved a similar equation using a relaxation scheme. Relaxation schemes, though, have poor convergence and usually require many iterations. Foster and Metaxas reported that they obtained good results even after a very small number of re- laxation steps. However, since we are using a different method to resolve for the advection step, we must use a more accurate method.
x
0
p(x,s)
p(x,−∆t) s −∆t
.u=0
Figure 1: One simulation step of our solver is composed of steps. The first three steps may take the field out of the space of divergent free fields. The last projection step ensures that the field is divergent free after the entire simulation step.
Figure 2: To solve for the advection part, we trace each point of the field backward in time. The new velocity at is therefore
the velocity that the particle had a time .
2.2 Method of Solution
ago at the old location
Eq. 5 is solved from an initial state
through time with a time step . Let us assume that the field has been resolved at a time and that we wish to compute the field at a later time . We resolve Eq. 5 over the time span in four steps. We start from the solution of the previous time step and then sequentially resolve each term on the right hand side of Eq. 5, followed by a projection onto the divergent free fields. The general procedure is illustrated in Figure 1. The steps are:
The solution at time is then given by the last velocity field: . A simulation is obtained by iterating these
steps. We now explain how each step is computed in more detail. The easiest term to solve is the addition of the external force . If we assume that the force does not vary considerably during the
time step, then
is a good approximation of the effect of the force on the field over the time step . In an interactive system this is a good approxi- mation, since forces are only applied at the beginning of each time step.
The next step accounts for the effect of advection (or convec- tion) of the fluid on itself. A disturbance somewhere in the fluid propagates according to the expression . This term makes the Navier-Stokes equations non-linear. Foster and Metaxas resolved this component using finite differencing. Their method is stable only when the time step is sufficiently small such that
by marching
!.
87
9!6
?@9"*+7F/+7*+793%&
I6'1!7F67
% 7 #
I 6 L 7
!
? I 6 7F6 % 7
45 667+45320
*+7 :
'
3%7F6.*+7K%&#J:;H !
%
?%&%7 )I6
?
./+7 )*+7 $ % 7 $ 7 "#-.,)'(&%"#!
6
./+7F0! 6
7
Indeed, the method of characteristics is more precise when the field is close to divergent free. More importantly from a visual point of view, the projection step forces the fields 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, the complexity of our proposed algorithm is of complexity . Foster and Metaxas’ solver has the same complexity. This perfor- mance is theoretically optimal since for a complicated fluid, any algorithm has to consult at least each cell of the computational grid.
2.3 Periodic Boundaries and the FFT
When we consider a domain with periodic boundary conditions, our algorithm takes a particularly simple form. The periodicity allows us to transform the velocity into the Fourier domain:
Figure 3: The values of the discretized fields are defined at the cen- ter of the grid cells.
fluid and a texture coordinate. The evolution of this scalar field is conveniently described by an advection diffusion type equation:
where is a diffusion constant, is a dissipation rate and is a source term. This equation is very similar in form to the Navier- Stokes equation. Indeed, it includes an advection term, a diffusion term and a “force term” . All these terms can be resolved exactly in the same way as the velocity of the fluid. The dissipation term not present in the Navier-Stokes equation is solved as follows over a time-step:
Similar equations were used by Stam and Fiume to simulate fire and other gaseous phenomena [21]. However, their velocity fields were not computed dynamically.
We hope that the material in this section has convinced the reader that our stable solver is indeed based on the full Navier-Stokes equations. Also, we have pointed to the numerical techniques which should be used at each step of our solver. We now proceed to describe the implementation of our model in more detail.
3 Our Solver 3.1 Setup
Our implementation handles both the motion of fluids and the prop- agation by the fluid of any number of substances like mass-density, temperature or texture coordinates. Each quantity is defined on ei- ther a two-dimensional (NDIM=2) or three-dimensional (NDIM=3) grid, depending on the application. The grid is defined by its phys- ical dimensions: origin O[NDIM] and length L[NDIM] of each side, and by its number of cells N[NDIM] in each coordinate. This in turn determines the size of each voxel D[i]=L[i]/N[i]. The definition of the grid is an input to our program which is speci- fied by the animator. The velocity field is defined at the center of each cell as shown in Figure 3. Notice that previous researchers, e.g., [7], defined the velocity at the boundaries of the cells. We prefer the cell-centered grid since it is more straightforward to im- plement. We allocate two grids for each component of the velocity: U0[NDIM] and U1[NDIM]. At each time step of our simulation one grid corresponds to the solution obtained in the previous step. 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 field corresponding to a substance transported by the flow. Although our implementation can handle any number of substances, for the sake of clarity we present only the algorithm for one field in this section. This scalar quantity is stored in the grids S0 and S1. The speed of interactivity is controlled by a single time step dt, which can be as large as the animator wishes, since our algorithm is stable.
In the Fourier domain the gradient operator “
the multiplication by , where . Consequently, both the diffusion step and the projection step are much simpler to solve. Indeed the diffusion operator and the projection operators in the Fourier domain are
FourierStep( add force:
advect: transform: diffuse: project: transform:
, ,
):
. The operator
projects the vector
onto the
Since the Fourier transform is of complexity
method is theoretically slightly more expensive than a method of solution relying on multi-grid solvers. However, this method is very easy to implement. We have used this algorithm to generate the “liquid textures” of Section 4.
2.4 Moving Substances through the Fluid
A non-reactive substance which is injected into the fluid will be ad- vected by it while diffusing at the same time. Common examples of this phenomenon include the patterns created by milk stirred in cof- fee or the smoke rising from a cigarette. Let be any scalar quantity which is moved through the fluid. Examples of this quantity include the density of dust, smoke or cloud droplets, the temperature of a
” is equivalent to
where
plane which is normal to the wave number . The Fourier transform of the velocity of a divergent free field is therefore always perpen- dicular to its wavenumbers. The diffusion can be interpreted as a low pass filter whose decay is proportional to both the time step and the viscosity. These simple results demonstrate the power of 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 tracer and a fast Fourier transform (FFT).
, this
?I6&61!6&2)-+1!
)1.)+)*' )0./!>&$)-+,%&)*'(!>&D&
)0.
? I
&
%#$” !&
7
< 66
&
/+7/+7 *+7< /+7
% #$! %+7%76%7.* +7 '61!77FFL7% 7
7 % 7 7 < $*7><
% # $! $
%
#
:
L
$
I 6
/+7 7
The physical properties of the fluid are a function of its viscosity visc alone. By varying the viscosity, an animator can simulate a wide range of substances ranging from glue-like matter to highly turbulent flows. The properties of the substance are modeled by a diffusion constant kS and a dissipation rate aS. Along with these parameters, the animator also must specify the values of these fields on the boundary of the grid. There are basically two types: peri- odic or fixed. The boundary conditions can be of a different type for each coordinate. When periodic boundary conditions are cho- sen, the fluid wraps around. This means that a piece of fluid which leaves the grid on one side reenters the grid on the opposite side. In the case of fixed boundaries, the value of each physical quantity must be specified at the boundary of the grid. The simplest method is to set the field 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 section we describe the boundary conditions chosen for each an- imation. For the special case when the boundary conditions are periodic in each coordinate, a very elegant solver based on the fast Fourier transform can be employed. This algorithm is described in 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 fluid is set into motion by applying external forces to it. We have written an animation system in which an animator with a mouse can apply directional forces to the fluid. The forces can also be a function of other substances in the fluid. For example, a temperature field moving through the fluid can produce buoyant and turbulent forces. In our system we allow the user to create all sorts of dependencies between the various fields, some of which are described in the results section of this paper. We do not describe our animation system in great detail since its functionality should be evident from the examples of the next section. Instead we focus on our simulator, which takes the forces and parameters set by the animator as an input.
3.2 The Simulator
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 Section 2 are in no way gratuitous but are immensely useful in cod- ing compact solvers. In particular, casting the problem into a math- ematical setting has allowed us to take advantage of the large body of work done in the numerical analysis of partial differential equa- tions. We have written the solver as a separate library of routines that are called by the interactive animation system. The entire li- brary consists of only roughly 500 lines of C code. The two main routines of this library update either the velocity field Vstep or a scalar field Sstep over a given time step. We assume that the external force is given by an array of vectors F[NDIM] and that the source is given by an array Ssource for the scalar field. The general structure of our simulator looks like
while ( simulating )
/* handle display and user interaction */
/* get forces F and sources Ssource from the UI */ Swap(U1,U0); Swap(S1,S0);
Vstep ( U1, U0, visc, F, dt );
Sstep ( S1, S0, kS, aS, U1, Ssource, dt );
The velocity solver is composed of four steps: the forces are added to the field, the field is advected by itself, the field diffuses due to viscous friction within the fluid, and in the final step the velocity is forced to conserve mass. The general structure of this routine is:
for(i=0;i