Model (Navier–Stokes)

For each grid cell $D$, $$\int_D \frac{\partial\vec{u}}{\partial t} \mathrm{d}V= \int_D \left(-\nabla\cdot(\vec{u}\vec{u}) + \nu\nabla^2\vec{u}-\nabla p\right) \mathrm{d}V$$ where pressure $p$ is also unknown.

Applying the Leibniz integral rule and the divergence theorem, $$\frac{\mathrm{d}}{\mathrm{d}t}\int_D \vec{u} \mathrm{d}V =\int_{\partial D} \left(-\vec{u}\vec{u} + \nu\nabla\vec{u} \right)\mathrm{d}S - \int_{D}\nabla p\ \mathrm{d}V$$

Algorithm (Finite Volume Method)

In each iteration,

1. $\displaystyle\vec{v}^{(n+1)} = \vec{u}^{(n)} + \Delta t \frac{1}{\Delta x \Delta y}\int_{\partial D} \left(-\vec{u}^{(n)}\vec{u}^{(n)} + \nu\nabla\vec{u}^{(n)} \right)\mathrm{d}S$

2. Solve $\displaystyle \nabla^2 p^{(n+1)} = \frac{\nabla\cdot\vec{v}^{(n+1)}}{\Delta t}$ for $p$ using Gauss–Seidel.

3. $\displaystyle\vec{u}^{(n+1)} = \vec{v}^{(n+1)} - \Delta t \nabla p^{(n+1)}$