Brief review of the models used

The 2D incompressible Navier-Stokes equations are written into the non dimensional streamfunction (s) and vorticty (v) formulation:

$\begin{array}{*{20}{c}}\frac{\partial^2s}{\partial x^2}+ \frac{\partial^2s}{\partial y^2}=-v\\ \frac{1}{Re}\frac{\partial^2v}{\partial x^2}+\frac{1}{Re} \frac{\partial^2v}{\partial y^2}= \frac{\partial s}{\partial y}\frac{\partial v}{\partial x}-\frac{\partial s}{\partial x}\frac{\partial v}{\partial y}\end{array}$

where $Re$ is the Reynolds number and $x,y$ are the Cartesian coordinates of the inertial frame of reference. Applying central differencing the Laplacian operator of the streamfunction can be discretized as following:

$\frac{s_{i+1,j}-2s_{i,j}+s_{i-1,j}}{\Delta h^2}+\frac{s_{i,j}-2s_{i,j}+s_{i,j}}{\Delta h^2} =-v_{i,j}$

where $\Delta h=\Delta x=\Delta y$ is the uniform space spacing. Re-ordering LHS it is obtained an approximation of the stream function:

$s_{i,j}=\frac{s_{i+1,j}+s_{i-1,j}+s_{i,j+1}+s_{i,j-1}+\Delta h^2v_{i,j}}{4}$

The streamfunction equation is solved by means of Successive Over Relaxation (SOR) method to obtain the steady state solution:

$s_{i,j}=\beta\frac{s_{i+1,j}+s_{i-1,j}+s_{i,j+1}+s_{i,j-1}+\Delta h^2v_{i,j}}{4}+\left(1-\beta\right)s_{i,j}$

where $\beta$ is the over-relaxation parameter. After the (iterative) value of the streamfunction has been computed by means of the above SOR approach, the vorticity is computed similarly:

$v_{i,j}=\beta\left[\frac{v_{i+1,j}+s_{i-1,j}+s_{i,j+1}+s_{i,j-1}}{4}+ Re\frac{\left(s_{i,j+1}-s_{i,j-1}\right)\left(v_{i+1,j}-v_{i-1,j}\right)- \left(s_{i+1,j}-s_{i-1,j}\right)\left(v_{i,j+1}-v_{i,j-1}\right)}{16}\right]+ \left(1-\beta\right)v_{i,j}$

where central differencing has been used also for first order partial derivative.

The above SOR equation for $s$ and $v$ are coupled with suitable boundary conditions, see T. Stortkuhl, C. Zenger and S. Zimmer. The iteration on each cell $(i,j)$ ends when the residuals become lower than a user-defined threshold value.