Graphics Reference
In-Depth Information
the old value at that grid point. Here we're including ghost values of
q
in
the boundaries for some of those averages. Therefore, the maximum new
value of
q
has to be less than or equal to the maximum old value of
q
,and
similarly the minimum new value of
q
has to be greater than or equal to
the minimum old value. So unstable growth is impossible. A more detailed
analysis further can prove that spurious oscillations (the unphysical vibra-
tion we could hit before) are ruled out as well. This is all true no matter
how large Δ
t
is taken: it's unconditionally stable and monotone!
7
For the full 3D problem, with constant viscosity, things work out much
the same. For each component of velocity independently we need to solve a
system with a matrix
I
+(Δ
tη
)
/
(
ρ
Δ
x
2
)
A
where
A
is nearly of the same form
as the Laplacian matrix in the pressure solve: 6 down the main diagonal
and
1 for neighboring grid cells. The one wrinkle we need to address is
handling the different boundary conditions. As we said before, at a no-slip
solid wall, we can just substitute in the known solid velocity, eliminating
that fluid velocity from the equations. Free surfaces and no-stick solid walls
are somewhat more tricky. The simplest approach is to use ghost velocity
values equal to the neighboring fluid-velocity values. For example, when
setting up the (
i
+1
/
2
,j,k
)th equation for
u
in the liquid, where
u
i
+3
/
2
,j,k
is in the air, use the ghost value
−
u
i
+3
/
2
,j,k
=
u
i
+1
/
2
,j,k
instead. Unfortunately, not only does this have the rigid rotation issue that
arose in explicit treatment, it's not even convergent for non-grid-aligned
boundaries—much the same way the classic voxelized pressure solve fairs
poorly on non-grid-aligned boundaries. We leave a solution to this prob-
lem for further research. Finally, for no-stick inviscid solid walls, we can
substitute in either the solid velocity (if it's the component normal to the
grid cell face) or use the ghost-fluid velocity (if it's a component tangential
to the grid cell face)—or again, some weighted average of the two for a
partially viscous boundary.
Finally we turn to the variable viscosity case. This can also be done
implicitly by making the components of
τ
on the staggered grid depend
on the new velocities.
8
We don't want to actually include the
τ
values as
7
This of course doesn't mean we're necessarily getting the correct answer if Δ
t
is very
large—there is still an approximation error. However, a Fourier analysis can show that
for large Δ
t
the only “problem” is that the diffusion, the exponential decay of all the
deformation modes, is effectively slower than it should be, though it still takes place.
8
For non-Newtonian fluids where
η
is a function of the strain rate, defining that
implicitly based on the new strain rate makes the equations non-linear; it simplifies life
