Graphics Reference
In-Depth Information
The simplest and cheapest method that makes sense to use is a second-
order Runge-Kutta integrator, an example of which we introduced before.
However, it's actually slightly unstable in regions of rigid rotation: parti-
cles will steadily migrate away from the center of rotation, perturbing the
distribution. A much better choice, albeit a little more expensive choice,
is to use a “three-stage third-order Runge-Kutta scheme:” these
are
sta-
ble for rigid rotations as long as Δ
t
isn't too large. There are infinitely
many such schemes, but Ralston [Ralston 62] showed the following scheme
is probably the best when it comes to minimizing error:
k
1
=
u
(
x
n
)
,
k
2
=
u
(
x
n
+
2
Δ
tk
1
)
,
k
3
=
u
(
x
n
+
4
Δ
tk
2
)
,
x
n
+1
=
x
n
+
9
Δ
tk
1
+
9
Δ
tk
2
+
9
Δ
tk
3
.
Here I've ignored the time dependence of the velocity field, which techni-
cally reduces the accuracy back to first order in time, just like our overall
time-splitting algorithm; the error we care more about here has to do with
the variation of velocity in space.
Chapter 3 also mentions the possibility of using substeps in time in-
tegration to better control error, say constraining each particle to move
roughly Δ
x
each substep. This can easily be incorporated into Runge-
Kutta schemes by first evaluating
u
(
x
n
), which will be used regardless of
the time-step size and then setting the substep size Δ
τ
to keep
below a threshold. Of course, substep sizes must be adjusted so that when
you add them up they exactly equal the global time step Δ
t
.
Δ
τu
(
x
n
)
A fairly
effective solution is given in Figure 10.1.
There's finally also the question of boundaries: what should happen to
a particle when it leaves the fluid? For the usual solid wall boundaries, this
is presumably just due to small numerical errors, which in the limit should
go to zero. The obvious fix if a particle does end up on the other side of a
solid wall is to project it out to the closest point on the boundary, perhaps
plus a small fraction of Δ
x
back into the fluid. This projection is a trivial
operation if the solid is described with a level set and not too dicult for
other geometry. The same may be true of a liquid free surface—unless the
particles themselves are helping to define where that surface is, in which
case they should be free to move. However, for open edges on the grid
with a
p
= 0 “free surface” condition, or inflow/outflow boundaries where
we specify
u
n
to be something different than the solid's normal velocity,
it makes sense instead to simply delete particles that stray outside the
·
