Graphics Reference
In-Depth Information
scale = dt / (rho*dx);
loop over i,j,k where cell(i,j,k)==FLUID:
u(i,j,k) -= scale * p(i,j,k);
u(i+1,j,k) += scale * p(i,j,k);
v(i,j,k) -= scale * p(i,j,k);
v(i,j+1,k) += scale * p(i,j,k);
w(i,j,k) -= scale * p(i,j,k);
w(i,j,k+1) += scale * p(i,j,k);
loop over i,j,k where cell(i,j,k)==SOLID:
u(i,j,k) = usolid(i,j,k);
u(i+1,j,k) = usolid(i+1,j,k);
v(i,j,k) = vsolid(i,j,k);
v(i,j+1,k) = vsolid(i,j+1,k);
w(i,j,k) = wsolid(i,j,k);
w(i,j,k+1) = wsolid(i,j,k+1);
Figure 4.1.
Pseudocode for the pressure update.
there could be multiple values for it) but instead just use formulas like
Equation (4.6) whenever we need to know the pressure difference across a
boundary face.
Using the same convention for storage from Chapter 2, Equations (2.14)-
(2.17), the pressure update can be translated into code similar to Figure
4.1. Note that instead of looping over velocity locations and updating them
with pressure differences, we are looping over pressure values and updating
the velocities they affect. Also, we use the trick of directly setting solid wall
velocities instead of working it out as a pressure update. The terms such as
usolid(i,j,k)
may well be replaced with simple expressions rather than
actually stored in an array.
Boundary conditions can be complicated and are the usual culprit when
bugs show up. It could be worth your time going over this section slowly,
with a drawing of the MAC grid (like Figure 2.1) in front of you, looking at
different configurations of solid, fluid, and air cells until you feel confident
about all this.
4.2 The Discrete Divergence
Now for the easy part of the chapter! In the continuum case, we want our
fluid to be incompressible:
∇·
u
= 0. On the grid, we will approximate this
