Hybrid Particle Methods - Fluid Simulation for Computer Graphics - page 139

This is practically begging for a particle treatment! In this section, we'll
go through the basics of the vortex-particle method of Chorin [Chorin 73],
introduced to graphics by Gamito et al. [Gamito et al. 95].
We'll begin with an initial velocity field, say on a MAC grid, and esti-
mate its vorticity at the grid cell corners:

v i +1 ,j +1 / 2 − v i,j +1 / 2 u i +1 / 2 ,j +1 − u i +1 / 2 ,j

ω i +1 / 2 ,j +1 / 2 = − .

Δ x Δ x

We then seed particles throughout the entire fluid region (a small number,
maybe one or four, per cell) and interpolate the grid vorticity at each. Each
particle's vorticity ω p will stay fixed for the duration of the simulation: all
we need to do is advect the particles in the velocity field. Velocity itself
needn't be advected, pressure needn't be solved, etc.
Of course, if we don't explicitly simulate velocity, we are left with the
tricky question of how to determine the velocity field from just the parti-
cles. We begin by transferring the vorticity of the particles to the grid cell
corners ω i +1 / 2 ,j +1 / 2 ,k +1 / 2 in essentially the same way we transferred smoke
concentration in the previous section. Then we solve for the incompressible
velocity field that possesses this vorticity.
Recalling an identity we used at the end of Chapter 9, any incompress-
ible velocity field can be represented as the curl of a potential: u =

ψ .
In two dimensions ψ is a scalar, and the curl operator simply rotates its
gradient by 90 degrees:

∇×

ψ = ∂ψ

.

∂ψ
∂x

u = ∇×

∂y ,

−

In fact, since the gradient of ψ is then normal to the velocity field, the
isocontours of ψ are just the streamlines of the flow, and we can call ψ the
streamfunction . Storing the ψ values at the grid cell corners, we can easily
get the discrete velocity field from

ψ i +1 / 2 ,j +1 / 2 − ψ i +1 / 2 ,j− 1 / 2

Δ x

u i +1 / 2 ,j = ,

ψ i +1 / 2 ,j +1 / 2 − ψ i− 1 / 2 ,j +1 / 2

Δ x

v i,j +1 / 2 = − .

Finally, to determine the streamfunction we solve for one whose velocity
field has the desired vorticity:

∇× ( ∇× ψ )= ω.