Graphics Reference
In-Depth Information
and cosines, which implies the Fourier coecients
u
ijk
may be complex even
though
u
is real-valued: this helps to simplify some of the other notation,
and more to the point, matches the API of most Fast Fourier Transform
packages. Note also that the Fourier coecients
u
ijk
are 3D vectors of
complex numbers, not just scalars: you can think of this really being the
Fourier transform of
u
, the separate Fourier transform of
v
, and the further
separate Fourier transform of
w
all wrapped up into one equation.
In practice, of course, we'll use just a discrete Fourier transform, over
an
m
m
array of Fourier coecients. The length
L
should be chosen
large enough that the periodic tiling isn't too conspicuous, but not too large
relative to the simulation grid spacing Δ
x
—after all, we won't be able to
afford to take
m
too large, and we want the Fourier grid spacing
L/m
(the
smallest details we'll procedurally add) to be a lot smaller than Δ
x
.
Shinya and Fournier [Shinya and Fournier 92] and Stam and Fiume
[Stam and Fiume 93] introduced to graphics perhaps the simplest physically
reasonable turbulence model, the Kolmogorov “5/3-law”. This states that
for fully developed steady-state turbulence, the kinetic energy contained
in the Fourier modes of spatial frequency
k
should scale like
k
−
5
/
3
.This
means that the (
i, j, k
) Fourier coecient should have magnitude on the
order of
×
m
×
(
i
2
+
j
2
+
k
2
)
−
11
/
12
u
ijk
∼
for
i, j, k
large enough (the low spatial frequencies are assumed to not
belong to the isotropic turbulence regime.) We can take it, in fact, to be a
random vector uniformly sampled from a ball of radius
C
(
i
2
+
j
2
+
k
2
)
−
11
/
12
,
for some user-tunable constant
C
. The cut-off frequency, below which we
keep
u
ijk
zero, should be on the order of
L/
Δ
x
or more, so that we don't
add procedural details at scales that were captured in the simulation.
The divergence of velocity becomes a simple algebraic operation on the
Fourier series:
√
−
u
(
x
)=
∞
12
π
L
u
ijk
]
e
√
−
12
π
(
ix
+
jy
+
kz
)
/L
.
∇
[(
i, j, k
)
·
i,j,k
=
−∞
Therefore,
u
= 0 is equivalent to requiring that each Fourier coe-
cient
u
ijk
is perpendicular to its
wave vector
(
i, j, k
). Makingavelocity
field divergence-free then simplifies to just fixing each coecient individu-
ally, subtracting off their component in the radial direction—a simple 3D
projection.
Finally, once the Fourier coecients have been determined, an inverse
FFTcanbeappliedoneachofthe
u
-,
v
∇·
−
,and
w
-components to get a
