Consider two cells of unit area, touching along a boundary of unit length. We have the average mass and momentum of each cell, and we want to know how these will change in a unit time step. (Using units based on the cell size and time step makes things simpler. Scaling factors can be folded into constants, saving some multiplications.)

The flow of mass is density times velocity. And since density times velocity is the momentum density, the flow in unit time is

mass flow = dot(momentum, flow direction)

The flow of momentum is only slightly more tricky. It’s the mass flow times the velocity, and velocity is momentum divided by density, so

momentum flow (parallel) = mass flow * mass flow / density
momentum flow (perpendicular) = dot(momentum, perpendicular) * mass flow / density

All these values are interpolated at the midpoint of the two cells.

Each cell has several neighbours, so updating the grid is a two pass process. First, calculate all the flows. Then, apply them to the cells. To update the grid,

source density -= mass flow
dest density += mass flow
source momentum -= momentum flow
dest momentum += momentum flow

Note that the same amount is added to one cell as is subtracted from the other, so both mass and momentum are conserved. The only problem comes if one of the densities becomes negative. There are various ways of dealing with this, but just making sure that no more mass can flow than is actually in the cell is the easiest.

This is all fine, but run it and nothing actually happens because all the initial momentums are zero. Trying to start an explosion by putting a huge density in one cell doesn’t do anything either. Why is that? If each cell is full of particles and the average momentum is zero, that doesn’t mean the momentum of every particle is zero. In fact, they will be flying in all directions and it just happens that it all averages out. If a cell has more particles than its neighbour, some particles will just randomly cross over, and fewer will come back. So there is a flow proportional to the difference in density between the cells.

mass flow += (source density - dest density) * diffusion constant

If there’s mass flow, then there’s momentum flow, from the above formula.

Once things get moving they aren’t necessarily going to stop. You can end up with a sort of lava lamp effect as the mass flows around. Applying some damping can sort this out. For each cell,

momentum *= damping factor

And draining the excess mass can be useful too,

density = density * drain factor + (1 - drain factor)

Also, some friction between cells helps to smooth out the flow.

momentum flow (perpendicular) += momentum difference * friction constant

What are the neighbours of a cell? Technically, only the four touching cells. But eight-way flow looks better. Hexagons would be an alternative, but I like the cells to line up with the walls.

Boundary conditions are important. A cell can be part of the active area, or it can be behind a wall. Nothing can flow through a wall, so the flow is zeroed if a cell or its neighbour are out of the active area. It is, however possible for momentum to bounce off a wall.

So the idea is simple, but there are several subtleties to the implementation. It helps to have an extra row of wall cells all the way around the grid, and the flow calculation stops one short of the edge. For rendering, I want the cell corners, so I resample the grid, averaging the four cells around each corner.

The speed is not bad. A full screen of cells at 8 pixel resolution runs at 60 fps. Still, I wanted it faster and the algorithm seems a perfect match for vectorization. So that’s the next post.