Viscous Fluids - Fluid Simulation for Computer Graphics - page 113

8.5.1 Explicit Treatment

The simplest thing of all is to just use central differences on the given
velocity field. From the definition τ = η (

u T )weget ∇ u + ∇

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

Δ x

τ 11

i,j,k =2 η i,j,k

τ 12

i +1 / 2 ,j +1 / 2 ,k =
η i +1 / 2 ,j +1 / 2 ,k u i +1 / 2 ,j +1 ,k −

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

+

Δ x Δ x

τ 13

i +1 / 2 ,j,k +1 / 2 =
η i +1 / 2 ,j,k +1 / 2 u i +1 / 2 ,j,k +1 −

u i +1 / 2 ,j,k w i +1 ,j,k +1 / 2 − w i,j,k +1 / 2

+

Δ x Δ x

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

Δ x

τ 22

i,j,k =2 η i,j,k

τ 23

i,j +1 / 2 ,k +1 / 2 =
η i,j +1 / 2 ,k +1 / 2 v i,j +1 / 2 ,k +1 −

v i,j +1 / 2 ,k w i,j +1 ,k +1 / 2 − w i,j,k +1 / 2

+

Δ x Δ x

w i,j,k +1 / 2 − w i,j,k− 1 / 2

Δ x

τ 33

i,j,k =2 η i,j,k

This can be simplified in the obvious way for 2D.
In the case where η is constant and u n is discretely divergence-free,
it's not too hard to verify that just as in the differential equations, the
coupling terms between velocity components will drop out, simplifying the
whole update to a discrete Laplacian for each velocity component. For

brevity, here is just the u -component in 2D:

u i +1 / 2 ,j = u i +1 / 2 ,j
+ Δ tη
ρ

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

4 u i +1 / 2 ,j

Δ x 2

which should look familiar from the discretization of the Laplacian in the
pressure problem.