A streamfunction–vorticity solver for the Navier–Stokes equations

In two dimensions the incompressible, constant-density Navier–Stokes equations can be written ${\partial }_{t}\omega +\mathbf{\text{u}}\cdot \nabla \omega =\nu {\nabla }^{2}\omega$ ${\nabla }^{2}\psi =\omega$ with $\nu$ the viscosity coefficient. The vorticity $\omega$ and streamfunction $\psi$ are defined as $\omega ={\partial }_{x}{u}_{y}-{\partial }_{y}{u}_{x}$ ${u}_{x}=-{\partial }_{y}\psi$ ${u}_{y}={\partial }_{x}\psi$ The equation for the vorticity is an advection–diffusion equation which can be solved using the flux–based advection scheme in `advection.h`. The equation for the streamfunction is a Poisson equation.

``````#include "advection.h"
#include "poisson.h"``````

We allocate the vorticity field $\omega$, the streamfunction field $\psi$ and a structure to store the statistics on the convergence of the Poisson solver. The fields advected by the advection solver are listed in `tracers`.

``````scalar ω[], ψ[];
mgstats mgpsi;
scalar * tracers = {ω};``````

Here we set the default boundary conditions for the streamfunction. The default convention in Basilisk is no-flow through the boundaries of the domain, i.e. they are a streamline i.e. $\psi =$constant on the boundary.

``````ψ[right]  = dirichlet(0);
ψ[left]   = dirichlet(0);
ψ[top]    = dirichlet(0);
ψ[bottom] = dirichlet(0);``````

We set the default value for the `CFL` (the default in `utils.h` is 0.5). This is done once at the beginning of the simulation.

``````event defaults (i = 0)
CFL = 0.8;``````

At every timestep we update the streamfunction field $\psi$ by solving a Poisson equation with the updated vorticity field $\omega$ (which has just been advected/diffused).

``````event velocity (i++)
{
mgpsi = poisson (ψ, ω);``````

Using the new streamfunction, we can then update the components of the velocity field. Since they are defined on faces we need to average the gradients, which gives the discrete expression below (for the horizontal velocity component). The expression for the vertical velocity component is obtained by automatic permutation of the indices but requires a change of sign: this is done through the pseudo-vector `f`.

``````  trash ({u});
struct { double x, y; } f = {-1.,1.};
foreach_face()
u.x[] = f.x*(ψ[0,1] + ψ[-1,1] - ψ[0,-1] - ψ[-1,-1])/(4.*Δ);
boundary ((scalar *){u});
}``````