Viscous Fluids - Fluid Simulation for Computer Graphics - page 118

unknowns in the linear system, however; we eliminate them by plugging in
their definition into the velocity equations. This gives three sets of coupled
equations, of the form:

Δ t
ρ Δ x 2

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

⎛ ⎞

2 η i +1 ,j,k ( u i +3 / 2 ,j,k − u i +1 / 2 ,j,k )

⎝ ⎠

2 η i,j,k ( u i +1 / 2 ,j,k − u i− 1 / 2 ,j,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 )

× ,

η i +1 / 2 ,j− 1 / 2 ,k ( u i +1 / 2 ,j,k − u i +1 / 2 ,j− 1 ,k + v i +1 ,j− 1 / 2 ,k − v i,j− 1 / 2 ,k ) −

+ η 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 )

η i +1 / 2 ,j,k− 1 / 2 ( u i +1 / 2 ,j,k − u i +1 / 2 ,j,k− 1 + w i +1 ,j,k− 1 / 2 − w i,j,k− 1 / 2 ) −

Δ t
ρ Δ x 2

v i,j +1 / 2 ,k = v i,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 )

η i− 1 / 2 ,j +1 / 2 ,k ( u i− 1 / 2 ,j +1 ,k − u i− 1 / 2 ,j,k + v i,j +1 / 2 ,k − v i− 1 ,j +1 / 2 ,k ) −

+2 η i,j +1 ,k ( v i,j +3 / 2 ,k − v i,j +1 / 2 ,k )

× ,

−

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

v i,j− 1 / 2 ,k )

+ η 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 )

η i,j +1 / 2 ,k− 1 / 2 ( v i,j +1 / 2 ,k − v i,j +1 / 2 ,k− 1 + w i,j +1 ,k− 1 / 2 − w i,j,k− 1 / 2 ) −

w i,j,k +1 / 2 = w i,j,k +1 / 2 + Δ t
ρ

⎛
⎝

⎞
⎠

η 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 )

η i− 1 / 2 ,j,k +1 / 2 ( u i− 1 / 2 ,j,k +1 − u i− 1 / 2 ,j,k + w i,j,k +1 / 2 − w i− 1 ,j,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 )

× .

η i,j− 1 / 2 ,k +1 / 2 ( v i,j− 1 / 2 ,k +1 − v i,j− 1 / 2 ,k + w i,j,k +1 / 2 − w i,j− 1 ,k +1 / 2 ) −

+2 η i,j,k +1 ( w i,j,k +3 / 2 − w i,j,k +1 / 2 )

−

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

w i,j,k− 1 / 2 )

Yes, that's them in all their gory magnificence. It's not too hard to ver-
ify that, thanks to our use of the staggered grid, we still get a symmetric
positive definite matrix to solve (for all three components of velocity si-

alottoinstead lag the value of η by basing it on the old strain rate, giving linear
equations which nonetheless should give you stability.