src/test/swirl.c

Boundary layer on a rotating disk

Von Kármán, 1921 showed that the steady flow of an incompressible liquid of kinematic viscosity \nu induced by an infinite plane disk rotating at angular velocity \Omega can be described by a similarity solution. In effect, using \zeta=z\sqrt{\Omega/\nu} and setting the axial velocity U, radial velocity V and azimuthal velocity W as \displaystyle U=\sqrt{\nu \Omega} F(\zeta) \quad V=\Omega r H(\zeta) \quad W = \Omega r G(\zeta) the Navier-Stokes equations reduce to a couple of ODEs: \displaystyle F'''-F\, F'' +F'^2/2 +2G^2 = 0 \quad \mathrm{and} \quad G''-F\,G'+G\,F' = 0 with boundary conditions \displaystyle F(0)=F'(0)=0 \: G(0)=1. \quad \mathrm{and} \quad F'(\infty)=G(\infty)=0. where the prime denotes differentiation with respect to \zeta.

To reproduce this solution numerically, we use the axisymmetric Navier–Stokes solver with azimuthal velocity (swirl).

#include "grid/multigrid.h"
#include "axi.h"
#include "navier-stokes/centered.h"
#include "navier-stokes/swirl.h"

The left boundary is the rotating disk with \Omega = 1 and a no-slip condition for the tangential velocity i.e.

u.t[left] = dirichlet(0);
w[left]   = dirichlet(y);

We use an open (outflow) boundary condition for the right boundary.

u.n[right] = neumann(0);
p[right] = dirichlet(0);
pf[right] = dirichlet(0);

The top boundary condition is more tricky but the following seems to work.

u.n[top] = neumann(0);
p[top] = neumann(0);

We use a constant viscosity but it needs to be weighted by the (axisymmetric) metric.

face vector muv[];

event properties (i++) {
  foreach_face()
    muv.x[] = 0.2*fm.x[];
}

The computational domain is 12\times 12 and we limit the timestep.

int main()
{
  size (12);
  N = 128;
  mu = muv;
  DT = 2e-2;
  run();
}

We wait until the boundary layer is fully developed and quasi-stationary. We only consider values close to the origin to minimize the influence of boundaries (von Kármán’s solution is valid in an infinite domain).

event end (t = 20)
{
  foreach()
    if (x*x + y*y < 8)
      fprintf (stderr, "%g %g %.4g %.4g\n", x, y, u.x[], w[]);
}

set xlabel '{/Symbol z}'
set key center right
nu = 0.2
plot [0:6]'analytical' u 1:2 w l t '-F({/Symbol z})',	     \
          'log' u ($1/sqrt(nu)):(-$3/sqrt(nu)) t 'Basilisk', \
          'analytical' u 1:3 w l t 'G({/Symbol z})',	     \
          'log' u ($1/sqrt(nu)):($4/$2) t 'Basilisk'

Axial F(\zeta) and azimuthal G(\zeta) dimensionless velocity components (script)

References

See also

• Same case with Gerris