Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (2008)

M. Gross and D. James (Editors)

Evolving Sub-Grid Turbulence for Smoke Animation

H. Schechter†1 and R. Bridson‡1

1 The University of British Columbia, Canada

Abstract
We introduce a simple turbulence model for smoke animation, qualitatively capturing the transport, diffusion, and
spectral cascade of turbulent energy unresolved on a typical simulation grid. We track the mean kinetic energy
per octave of turbulence in each grid cell, and a novel “net rotation” variable for modeling the self-advection of
turbulent eddies. These additions to a standard fluid solver drive a procedural post-process, layering plausible
dynamically evolving turbulent details on top of the large-scale simulated motion. Finally, to make the most of the
simulation grid before jumping to procedural sub-grid models, we propose a new multistep predictor to alleviate
the nonphysical dissipation of angular momentum in standard graphics fluid solvers.
Categories and Subject Descriptors (according to ACM CCS): I.3.7 [Computer Graphics]: Three-Dimensional
Graphics and Realism

1. Contributions
Animating turbulent fluid velocity fields, from the delicate
swirls of milk stirred into coffee to the violent roiling of a
volcanic eruption, poses serious challenges. While the look
of the large-scale components of motion are well captured by
direct simulation of the fluid equations, increasing the grid
resolution to capture the smallest turbulent scales hits a se-
vere scalability problem. The usual solution of augmenting
a coarse simulation with procedurally synthesized turbulent
detail generally is limited by the visual implausibility of this
detail: while some statistics or invariants may be matched,
visually important aspects of the time evolution of turbu-
lence are still missing.
This paper attempts to bridge the gap by introducing a tur-
bulence model that simulates the transfer of energy in sub-
grid scales, then providing a procedural method for instanti-
ating a high resolution turbulent velocity field from this data,
ferred to procedural models when feasible, we address one
which can be evaluated on the fly for a marker particle pass
of the chief remaining sources of nonphysical energy dissi-
when rendering. For each octave of sub-grid detail, in each
pation in graphics fluid solvers, adding a predictor step to the
grid cell, we evolve both a kinetic energy density E and a net
usual time splitting of the incompressible Euler equations.
rotation θ for generating the turbulence.
This helps make the most of a limited grid size.
In addition, since full fluid simulation is generally pre-
2. Background
† e-mail: [email protected] The field of turbulence modeling is vast; for an overview we
‡ e-mail: [email protected] refer readers to the recent text by Pope [Pop05] or the earlier

c The Eurographics Association 2008.

 H. Schechter & R. Bridson / Evolving Turbulence

sub-grid turbulence models, i.e. directly simulating the large

scales while estimate the effect of unresolved smaller scales.
These methods separate the velocity field into a sum of
the “mean flow” (the large-scale components) and a turbu-
lent fluctuation. For the Reynolds Averaged Navier-Stokes
(RANS) approach, the mean-flow averaging is taken over
appropriate time scales, and for the Large Eddy Simulation
(LES) approach the averaging is done with a spatial ker-
Figure 1: Frames from an animation with the new turbu-
nel such as a Gaussian. Averaging the Navier-Stokes equa-
lence model.
tions gives rise to very similar equations for the mean flow,
but with the addition of a “Reynolds stress” giving the ef-
fect of turbulent mixing on the mean flow—at its simplest
classic by Tennekes and Lumley [TL72]. For an introduc- modeled as a nonlinear viscosity. For this paper, where we
tion to basic fluid simulation techniques within graphics, see try to capture the qualitative look but don’t provide quan-
Bridson and Müller-Fischer’s course notes [BMF07]. titatively accurate results, we assume numerical dissipation
There is a huge disparity in length scales for turbulent in the usual first-order accurate fluid solvers dominates this
flow: e.g. in the atmosphere large-scale features may be mea- term and thus ignore it.
sured in kilometres, with the smallest features in millimetres.
Some turbulence methods go further by tracking informa-
A fine enough grid to capture everything (the idea behind
tion about the turbulent fluctuations to produce better esti-
the Direct Numerical Simulation approach in science) may
mates of their effect on the mean flow. We focus in particu-
be enormous: n3 grid cells with n > 1000. Unfortunately,
lar on the popular k − ε model, originating in work by Har-
turbulence fills the entire volume with fine-scale features,
low and Nakayama [HN67]. Here two more fields are added
eliminating the benefit of adaptive grid methods, and to re-
to the simulation, k being the kinetic energy density of the
solve the small scales time steps must be proportional to the
turbulent fluctations (energy per unit volume) and ε repre-
grid spacing, resulting in at least O(n4 ) work. This severely
senting the rate of energy dissipation in the viscous scales.
limits the scalability of straight simulation for capturing tur-
Unlike our new model, this does not explicitly track the cas-
bulence.
cade of energy from large scales to small scales, but like our
However, turbulence research has developed higher level model simulates the spatial transport, diffusion, and ultimate
models which promise efficient simulation, based on the ob- dissipation of turbulent energy.
servation that smaller scales in turbulence quickly become
isotropic and more or less statistically independent. This al- Turning to computer graphics, Shinya and
lows useful notions of average effect and evolution without Fournier [SF92] and Stam and Fiume [SF93] used the
need of resolving all details, with the chief visible actions Kolmogorov 5/3 law and the inverse Fourier transform
being: to generate incompressible velocity fields with plausible
statistics, possibly layered on top of a large-scale flow.
• Enhanced mixing, which when averaged over appropriate This was later also used by Rasmussen et al. [RNGF03]
length and time scales can be interpreted as a diffusion to break up smoothness in the interpolation of 3D velocity
process, with a so-called “eddy viscosity”. This is why fields from simulated 2D slices. Kniss and Hart [KH04]
stirring milk in coffee is much more effective than letting demonstrated that by taking the curl of vector valued noise
it sit still. functions, plausible incompressible velocity fields can be
• Transfer of energy from large-scale structures to smaller- created without need for Fourier transforms; Bridson et
scale eddies and ultimately to the smallest scales where al. [BHN07] extended the curl-noise approach to respect
molecular viscosity takes over to dissipate kinetic energy solid boundaries and conveniently construct larger-scale
into heat. This is the cause of such a large range of length patterns, and illustrated the attraction of spatial modulating
scales. turbulence. Most recently Kim et al. [KTJG08] combined
The simplest model that captures these, the Kolmogorov 5/3 curl-noise with a wavelet form of the 5/3 law, extending
power law, has long been used in graphics for synthesizing the energy density measured in the highest octave of a
velocity fields with plausible statistics. Assuming uniform, simulation into a turbulent detail layer and further advecting
steady, self-preserving turbulence, where the energy contin- the noise texture with the flow to capture spatial transport of
uously injected at the largest scales matches the energy dis- turbulent flow features, with marker particles for rendering.
sipated by viscosity at the smallest scales and an equilib-
Neyret’s work on advected textures [Ney03] is closest in
rium in the intermedate scales, this provides a simple for-
spirit to this paper. He introduces a single representative
mula for the amount of energy and rate of fluctuation in each
vorticity vector per simulation grid cell and per octave of
frequency band.
turbulence, advected with the flow. The vectors are grad-
Many researchers have augmented fluid simulations with ually blended from coarse octaves to fine octaves, and are

c The Eurographics Association 2008.

 H. Schechter & R. Bridson / Evolving Turbulence

used to rotate the gradients in flow noise [PN01]. The flow 3. Characterizing Sub-Grid Turbulence
noise is then used to distort a texture map (e.g. a smoke- Our goal is to simulate the average behavior of the sub-grid
like hypertexture), with smooth regeneration of texture de- turbulence, leaving the details of instantiating a plausible ve-
tail via blending of the original when total deformation has locity field to a post-simulation phase. In particular, we want
increased past a limit. We instead track (and conserve in all to describe the turbulence without recourse to higher reso-
but the finest octave) kinetic energy, include the spatial dif- lution grids. Note that we assume that the combination of
fusion of kinetic energy, and properly decouple the rotation FLIP with our new multistep time integration will capture
of nearby gradient vectors in flow noise (so not all of the gra- all grid-level turbulence directly in simulation, thus we only
dient vectors in one simulation cell are rotated with the same focus on the missing sub-grid part.
vorticity); we further combine it with curl-noise to produce
velocity fields for marker particle advection in rendering. We extend the k − ε approach of tracking mean kinetic en-
ergy density: instead of a single total kinetic energy density
Our predictor method for reducing nonphysical rotational per grid cell, we break it up into octaves corresponding to
dissipation in the large-scale fluid simulation also has an- spatial frequency bands (similar to the usual notion of “tur-
tecedents in CFD and graphics. The chief culprits are in bulence” textures in graphics). For example, with three oc-
the treatment of the advection term in the momentum equa- taves we store three kinetic energy densities in each grid cell
(1) (2) (3)
tion. Fedkiw et al. [FSJ01] observed that some of the strong (i, j, k): Ei jk , Ei jk , and Ei jk , using E instead of the tradi-
numerical dissipation of the trilinear interpolation semi- tional k to avoid confusion with grid cell indices. The first,
Lagrangian method introduced by Stam [Sta99] can be re- E (1) , corresponds to the components with wavelengths ap-
duced by using a sharper Catmull-Rom-based interpolant in- proximately between ∆x and 12 ∆x, the E (2) value for wave-
stead; many other improvements to pure advection followed lengths between 12 ∆x and 41 ∆x, etc. This allows us to track
(e.g. Kim et al.’s BFECC approach [KLLR05]) culminating the transfer of energy in spectrum as well as space—opening
in Zhu and Bridson’s FLIP method [ZB05] which essentially the door to handling the transition from laminar to turbulent
eliminates all numerical dissipation from advection via use flow, or the unsteady evolution and decay of unsustained tur-
of particles. However, this is not the only source of dissipa- bulence. This still has attractive scalability: to increase the
tion: the first order time splitting used in graphics, where ve- apparent resolution of a simulation by a factor of 2k our
locities are separately advected and then projected to be in- memory use only increases by O(n3 k). We cannot distin-
compressible, also introduces large errors. In particular, an- guish variations in turbulent energy at sub-grid resolution,
gular momentum is not conserved even with a perfect advec- but since the turbulent energy undergoes a diffusion process,
tion step, as can readily be seen if a rigidly rotating fluid is sub-grid variations shouldn’t be visually significant.
advected by 90◦ in one time step: the angular velocity would
be entirely transferred into a divergent field, which pressure We also address the issue of temporal coherence in the
projection would subsequently zero out. In the scientific lit- turbulence field. While the energy density is enough to pro-
erature, predictor-corrector, multistep or Runge-Kutta meth- cedurally synthesize a plausible velocity field for a single
ods are used in conjunction with high resolution discretiza- frame, we want to reflect the correlation from one frame to
tions of the advection term and small time steps to avoid the next. Inspired by flow noise [PN01] we add an additional
(b)
this problem. However, these aren’t directly applicable to scalar variable θi jk in each grid cell for each octave, an es-
the large time steps taken with semi-Lagrangian advection timate of the average amount of rotation in that region of
or FLIP in graphics. space up to the current time caused by the turbulent compo-
nents, i.e. the time integral of the magnitude of vorticity at a
Vorticity confinement [FSJ01] helps recover some of the fixed location in space. We use this to directly control flow
lost angular momentum, by adding artificial forces to en- noise when instantiating velocities.
hance spin around local maxima of the existing vorticity.
Selle et al. [SRF05] extended this with spin particles, to 4. Evolving Sub-Grid Turbulence
boost the smallest-scale vorticity present on the grid. How-
ever neither technique is applicable to the simplest case, re- We break up the problem into plausibly evolving the kinetic
(b)
ducing the dissipation in rigid rotation, and thus we propose energies Ei jk and separately updating the θ angles.
a more general solution in this paper.
4.1. Energy Evolution
We also note that alternative approaches to fluid simu-
lation, in particular vortex particle methods (e.g. [YUM86, Note that the large-scale flow on the grid transports with it
Gam95, AN05, PK05]), do not suffer from this dissipation. small-scale turbulent components. Therefore, just as in the
However vorticity methods have their own share of prob- k − ε model we begin with advecting the kinetic energy den-
lems, particular in 3D, such as in handling free surface sities with the usual fluid variables in the simulation.
and solid boundary conditions, thus most work has centred In addition, the enhanced mixing of turbulent flow is usu-
around velocity/pressure methods. ally modeled as a nonlinear spatial diffusion term, spreading

c The Eurographics Association 2008.

 H. Schechter & R. Bridson / Evolving Turbulence

turbulent energy across space. We make a crude simplifica-

tion to constant-coefficient linear diffusion instead; from a
scientific viewpoint this is probably unjustifiable, yet we ar-
gue it visually captures the qualitative effect of eddy viscos-
ity while avoiding numerical complications caused by more
accurate nonlinear models.
Finally, the nonlinear advection term in the Navier-Stokes
momentum equation causes transfer of kinetic energy be-
tween different frequency bands. While this happens in both Figure 2: Our turbulence model shows an initial distur-
directions, it is widely modeled that the dominant effect in bance in the lowest octave (left) diffusing in space and cas-
turbulence is the “spectral cascade” of energy from low fre- cading energy to higher octaves later in time (right).
quencies to higher frequencies, with the majority of the en-
ergy transfer being between nearby wavelengths. We there-
fore include a simple transfer term from the kinetic energy
in one octave to the next octave along, similar to Neyret’s found, different initial conditions and parameters for turbu-
vorticity transfer. Again, the true energy transfer is a nonlin- lence evolution can be tried out very efficiently.
ear process and our constant-coefficient linear simplification
We ran examples of the new turbulence model coupled
is scientifically invalid, but is useful as the simplest possible
with a FLIP simulation of smoke [ZB05], with the standard
model that captures the qualitative visual effect of the spec-
addition of heat and buoyancy forces to the fluid solver. We
tral cascade.
advected the turbulent kinetic energy variables with the FLIP
We discretize this for a single time step on a grid as fol- particles, transferring them to the grid, applying diffusion
lows, with given coefficients α ≥ 0 and β ≥ 0 controlling and spectral cascade on the grid, and transferring the result
spatial diffusion and spectral cascade respectively: back to the particles in the usual FLIP way.
• Advect the per-octave kinetic energies with the large-scale
(b)
flow to get intermediate energies Ēi jk . 4.2. Net Rotation Evolution
• Apply spatial diffusion and scale-to-scale transfer to get (b)
the energies at the next time step: The net rotation θi jk is an estimate of the time integral of
 (b) the magnitude of vorticity stemming from the b’th octave of
(b) 
Ēi+1, jk + Ēi−1, jk turbulence in grid cell (i, j, k). This is not a quantity to be ad-
 +Ēi(b) (b) vected: it’s an Eulerian history variable. We also emphasize
j+1,k + Ēi j−1,k 
 
(b) 
(b) Ēi jk + α∆t 

(b) (b) this is just a scalar, unlike the vector-valued vorticity, since
Ei jk =  +Ēi jk+1 + Ēi jk−1 
 (1) we are using this to estimate average rotation over a region
(b)
−6Ēi jk of space (containing many differently-oriented vorticities),
(b−1)
+β∆t(Ēi jk
(b)
− Ēi jk ) analagous to our use of scalar kinetic energy density rather
than mean velocity.
There is a time step restriction of ∆t < (6α + β)−1 , which
At every time step of the turbulence evolution simulation,
in our examples was never particularly restrictive. Strictly (b)
speaking the α and β coefficients should take into account we increment θi jk by an estimate of the average magnitude
length scales and the energies themselves, but we choose of the vorticities in the b’th octave, which we assume are
to keep them as simple tunable constants. In our evolution proportional to the mean speed (available from the kinetic
energy density) divided by the wave length:
equation the first octave Ē (1) refers to a nonexistent Ē (0) :
this ideally should be energy pulled from the large-scale sim-
q
(b)
(b)
2Ei jk /ρ
ulation itself, but we instead simply seeded it with a spatial kωki jk ∼ (2)
texture, e.g. introducing turbulent energy just in the source 2−b ∆x
of a smoke plume. At the last octave, we optionally add an- This is dimensionally correct and consistent with the obser-
other term γ∆t Ē (b) with 0 ≤ γ ≤ β to partially cancel out vation that turbulence typically remains isotropic. Introduc-
the loss of energy to untracked octaves. See figure 2 for an ing a constant of proportionality δ we use the following up-
example of the effect of both diffusion and spectral cascade. date:
q
We underscore that this evolution doesn’t influence the (b)
(b) (b)
2Ei jk /ρ
large-scale motion at all, under the assumption that errors θi jk ← θi jk + δ∆t −b (3)
2 ∆x
in the large-scale fluid simulation will dominate the ef-
fect of sub-grid turbulence. This decoupling has the attrac- For the full effect we should also introduce a θ(0) variable
tive side-effect that turbulence evolution can be done as a that is updated by the magnitude of vorticity in the large-
post-process: once a suitable large-scale simulation has been scale simulation, but we have not yet experimented with this.

c The Eurographics Association 2008.

 H. Schechter & R. Bridson / Evolving Turbulence

This has in fact proven useful as a reasonably efficient state-

less pseudo-random number generator in other work.
To evaluate a given octave of this noise at a point, we look
up the rotated gradients at the corners of the lattice cell con-
taining the point. This lattice, for octave b, has a spacing of
2−b+1 ∆x compared to the ∆x spacing of the simulation grid.
Figure 3: On the left is the large-scale velocity field inter- At each lattice point we use the hash to look up a pair of
polated from the displayed grid (unrealistically smooth for orthogonal unit vectors, p̂ and q̂, trilinearly interpolate the
illustration purposes). On the right we add sub-grid turbu- net rotation θ̂ in this octave from the simulation, and finally
lent scales to get the total velocity field. form a rotated gradient vector as:
~g = cos θ̂ p̂ + sin θ̂q̂ (5)
One of the most expensive parts of this calculation is the
5. Instantiating Velocity Fields trigonometic function evaluation, but after all the modeling
Once we have completed a large-scale fluid simulation and errors in the system it’s not crucial that these be accurate; we
then turbulence evolution, we generate a plausible total ve- substituted the following cubic spline
locity field by interpolating from the large-scale velocity grid √

1

and adding to it the curl of a vector potential derived from the sin θ ≈ 12 3t t − (t − 1) (6)
2
turbulence quantities, as in the curl-noise approach of Kniss
and Hart [KH04] and Bridson et al. [BHN07]: see figure 3. where t is the fractional part of θ/(2π), and similarly ap-
proximated cos(θ) = sin(θ + 12 π). We found this was suffi-
Perlin and Neyret’s flow noise [PN01] was originally
ciently accurate and much faster. (Note that this spline has
designed to provide fluid-like textures which animate in
zeros at the correct roots of sin and attains a max and min
time with a pseudo-advection quality (unlike regular time-
of ±1 as expected.) Finally, we evaluate three independent
animated noise which smoothly changes without coherent
noise functions to get a vector noise value ~N (b) (~x,t).
flow structure): the gradient vectors underling Perlin noise
are rotated in time as if by fluid vorticity. Neyret further We also interpolate the kinetic energy density in an octave
used this for deforming other textures in his advected tex- at any point in space from the values on the simulation grid,
ture work [Ney03], and Bridson et al. [BHN07] suggested it using quadratic B-splines, to estimate an appropriate ampli-
is well suited for use in curl-noise to generate velocity fields tude for modulating the noise:
featuring vortices (corresponding to peaks in the noise func- q
tions) which naturally swirl around and interact with other A(b) (~x,t) = C2−b ∆x 2E (b) (~x,t)/ρ (7)
vortices rather than simply smoothly appearing and disap-
The user-defined constant C should be approximately 1, so
pearing. We therefore adopt flow noise to build the vec-
that the actual kinetic energy density of the velocity field be-
tor potential ~ψ for the turbulent velocity contribution; how-
low will match up to the stored E (b) . Adding up the octaves
ever, since vorticities in turbulence should be uniformly dis-
gives the vector potential ψ:
tributed due to isotropy, we modify Neyret’s scheme which
only used a single vorticity per grid cell per octave. ~ψ(~x,t) = ∑ A(b) (~x,t)~N (b) (~x,t) (8)
b
We begin with an auxiliary array of 128 pairs of uniformly
randomly selected orthogonal unit three-dimensional vec- Further modifications to respect solid boundaries can be
tors, p̂ and q̂, which will be used to parameterize a plane of made per Bridson et al. [BHN07]. The turbulent part of the
rotation at each point in the lattice Instead of Perlin’s hash total velocity field is the curl of ~ψ.
function to map lattice points to this auxiliary array, which
For rendering we typically run many smoke marker parti-
induces severe axis-aligned artifacts in frequency space (see
cles through the total velocity field (seeded in the appropriate
figure 8(a) and (b) from Cook and DeRose’s work [CD05])
places for the given simulation), evaluating velocity wher-
we use the following:
ever and whenever needed. With several octaves of noise, all
hash(i, j, k) = H(i xor H( j xor H(k))) mod 128 (4) the interpolation and gradient evaluations can be rather more
expensive than simple interpolation of velocity from a grid,
where H(s) provides a good quality pseudo-random permu-
even with optimizations such as our trigonometic approxi-
tation of the 32-bit integers:
mation. However, we do highlight that it is embarrasingly
• H(s): parallel—each marker particle can be advected without ref-
• s ← (s xor 2747636419) · 2654435769 mod 232 erence to any of the others—and the amount of data involved
• s ← (s xor bs/216 c) · 2654435769 mod 232 is fairly small relative to the apparent resolution of the tur-
• s ← (s xor bs/216 c) · 2654435769 mod 232 bulent velocity field. Further optimization, however, may be
• return s possible by using more memory: evaluating at least some of

c The Eurographics Association 2008.

 H. Schechter & R. Bridson / Evolving Turbulence

the octaves on higher resolution grids once and then inter-

polating from there, or evaluating on smaller tiles that are
cached for reuse by multiple particles.

6. Reducing Angular Dissipation in Simulation

We have so far looked at layering in sub-grid turbulence; it
is of course also imperative that as much simulation detail
as possible is retained in the large-scale simulation. We pro-
pose a simple multistep predictor to reduce the afore men-
tioned loss of angular momentum caused by the first order
time splitting of advection.
To motivate our method, we point out that the incompress-
ibility constraint can be arrived at as the infinite limit of the
second viscosity coefficient λ:
d~u λ
= ∇(∇ ·~u) (9)
dt ρ
(Note this is not the physically correct limit in which com-
pressible flow becomes asymptotically incompressible, but
rather a useful thought experiment which avoids the need
of considering variables other than velocity and equations
other than momentum.) As λ → ∞, any divergent motions
are damped out to nothing, giving back incompressible flow;
other motions, such as shearing, are unaffected.
For a stiff problem such as equation 9, an implicit integra-
tor with stiff decay is recommended. For example, if Back-
Figure 4: Top row: large-scale motion only. Middle row: tur-
wards Euler is used for the stiff viscosity term, and we then
bulence added, but without evolution. Bottom row: our full
take the limit of the solution as λ → ∞, we find we are back
turbulence evolution model.
to the usual pressure projection splitting method: ~un+1 is the
divergence-free part of the advected ~un .
To get higher accuracy we consider better implicit meth- 7. Results
ods which still posess stiff decay, such as BDF2. This is a
multistep method which, after again assuming the advection We ran several demonstrations of the sub-grid turbulence
term is handled separately, would discretize equation 9 as: model in 3D. The time added to the simulation for track-
ing the extra turbulence quantities was negligible, giving
4 1 2 λ roughly two seconds per frame for a 30 × 120 × 30 grid on
~un+1 = ~un − ~un−1 + ∆t ∇(∇ ·~un+1 ) (10)
3 3 3 ρ a 2.2Ghz Core Duo. Running the marker particles through
If we solve this, then take the limit as λ → ∞, we get pres- the total velocity field with three octaves of noise scaled lin-
sure projection at time n + 1 of 43 ~un − 31 ~un−1 , which can be early with particle count, at about .3 seconds per frame for
seen as a simple predictor for the new velocity. (We empha- 1000 particles; this became a bottleneck for higher quality
size it is the final time n + 1 velocity that is made divergence- renders with hundreds of thousands of particles, but as we
free: there is no divergence error in the velocity field in noted earlier should allow for trivial parallelization not to
which we advect particles.) mention other optimization possibilities.
Figure 1 shows frames from a simulation of a smoke
This analysis isn’t fully worked out, but our preliminary
plume generated by a continuous source of hot smoke, with
experiments with FLIP, where we transfer the predicted new
the full turbulence model in place. Marker particles were
particle velocities, based on their current and past values,
seeded in each frame, then rendered with self-shadowing.
to the grid for pressure projection instead of the current ve-
locity, are very promising: the method remains apparently Figure 4 demonstrates the effect of layering the turbulence
unconditionally stable, with significantly improved conser- model on top of existing large-scale motion, with or with-
vation of angular momentum yet without introducing noise. out the evolution (diffusion and spectral cascade). Without
Visually flows remain a lot more lively and areas of rotation turbulence the motion is flat and boring, at least in these
are damped much more slowly, without blowing up. early frames. With turbulence but without evolution, we

c The Eurographics Association 2008.

 H. Schechter & R. Bridson / Evolving Turbulence

were forced to include energy in all octaves right from the [KH04] K NISS J., H ART D.: Volume effects: model-
start; there thus appears to be too much fine-scale detail ini- ing smoke, fire, and clouds, 2004. Section from ACM
tially, but not enough mid-scale turbulent mixing. The full SIGGRAPH 2004 courses, Real-Time Volume Graphics,
evolution model instead shows the spectral cascade, with http://www.cs.unm.edu/ jmk/sig04_modeling.ppt.
the smallest-scale vortices only appearing naturally as en- [KLLR05] K IM B., L IU Y., L LAMA I., ROSSIGNAC J.:
ergy reaches them, and the spatial diffusion lets the turbulent FlowFixer: using BFECC for fluid simulation. In Proc.
mixing extend further into the flow breaking up symmetries Eurographics Workshop on Natural Phenomena (2005).
better.
[KTJG08] K IM T., T HÜREY N., JAMES D., G ROSS M.:
Wavelet turbulence for fluid simulation. ACM Trans.
8. Conclusions Graph. (Proc. SIGGRAPH) (2008).
We presented both a new sub-grid turbulence evolution [Ney03] N EYRET F.: Advected textures. In Proc. ACM
model and a new predictor step for fluid simulation, with SIGGRAPH/Eurographics Symp. Comp. Anim. (2003),
the goal of getting closer to high quality smoke simulation. pp. 147–153.
At present our model remains too simplistic for scientific va-
lidity, but we argue it captures much of the qualitative look [PK05] PARK S. I., K IM M. J.: Vortex fluid for gaseous
of real turbulence. phenomena. In Proc. ACM SIGGRAPH/Eurographics
Symp. Comp. Anim. (2005).
There are many directions for future work, beyond simply
[PN01] P ERLIN K., N EYRET F.: Flow noise.
improving the accuracy of our model: in particular, we high-
In ACM SIGGRAPH Technical Sketches and
light the question of how to automatically pull energy from
Applications (2001), p. 187. http://www-
the large-scale simulation into the first turbulent octave. Ide-
evasion.imag.fr/Publications/2001/PN01/.
ally a solution to this would not perturb stable laminar re-
gions of flow, but naturally cause a transition to turbulence [Pop05] P OPE S. B.: Turbulent Flows. Cambridge Uni-
when an instability is detected. versity Press, 2005.
[RNGF03] R ASMUSSEN N., N GUYEN D., G EIGER W.,
Acknowledgements F EDKIW R.: Smoke simulation for large scale phenom-
ena. ACM Trans. Graph. (Proc. SIGGRAPH) 22 (2003),
This work was supported in part by a grant from the Natural 703–707.
Sciences and Engineering Research Council of Canada, and
by a scholarship from BC Innovation Council and Precarn [SF92] S HINYA M., F OURNIER A.: Stochastic motion:
Incorporated. Motion under the influence of wind. In Proc. Eurograph-
ics (1992), pp. 119–128.

References [SF93] S TAM J., F IUME E.: Turbulent wind fields for
gaseous phenomena. In Proc. ACM SIGGRAPH (1993),
[AN05] A NGELIDIS A., N EYRET F.: Simulation of pp. 369–376.
smoke based on vortex filament primitives. In Proc. ACM
SIGGRAPH/Eurographics Symp. Comp. Anim. (2005). [SRF05] S ELLE A., R ASMUSSEN N., F EDKIW R.: A
vortex particle method for smoke, water and explosions.
[BHN07] B RIDSON R., H OURIHAN J., N ORDENSTAM ACM Trans. Graph. (Proc. SIGGRAPH) (2005), 910–
M.: Curl-noise for procedural fluid flow. ACM Trans. 914.
Graph. (Proc. SIGGRAPH) 26, 3 (2007).
[Sta99] S TAM J.: Stable fluids. In Proc. SIGGRAPH
[BMF07] B RIDSON R., M ÜLLER -F ISCHER M.: Fluid (1999), pp. 121–128.
simulation, 2007. ACM SIGGRAPH 2007 courses.
[TL72] T ENNEKES H., L UMLEY J. L.: A first course in
[CD05] C OOK R. L., D E R OSE T.: Wavelet noise. ACM turbulence. MIT Press, 1972.
Trans. Graph. (Proc. SIGGRAPH) 24, 3 (2005), 803–811.
[YUM86] YAEGER L., U PSON C., M YERS R.: Combin-
[FSJ01] F EDKIW R., S TAM J., J ENSEN H. W.: Vi- ing physical and visual simulation—creation of the planet
sual simulation of smoke. In Proc. SIGGRAPH (2001), jupiter for the film 2010. In Proc. ACM SIGGRAPH
pp. 15–22. (1986), pp. 85–93.
[Gam95] G AMITO M. N.: Two dimensional Simulation
[ZB05] Z HU Y., B RIDSON R.: Animating sand as a fluid.
of Gaseous Phenomena Using Vortex Particles. In Proc.
ACM Trans. Graph. (Proc. SIGGRAPH) (2005), 965–
of the 6th Eurographics Workshop on Comput. Anim. and
972.
Sim. (1995), Springer-Verlag, pp. 3–15.
[HN67] H ARLOW F. H., NAKAYAMA P. I.: Turbulence
transport equations. Phys. Fluids 10, 11 (1967), 2323–
2332.

c The Eurographics Association 2008.