src/test/couette.c

Couette flow between rotating cylinders

We test embedded boundaries by solving the (Stokes) Couette flow between two rotating cylinders.

#include "grid/multigrid.h"
#include "embed.h"
#include "navier-stokes/centered.h"
#include "view.h"

int main()
{
  origin (-L0/2., -L0/2.);
  
  stokes = true;
  TOLERANCE = 1e-5;
  for (N = 16; N <= 256; N *= 2)
    run();
}

scalar un[];

#define WIDTH 0.5

event init (t = 0) {

The viscosity is unity.

  μ = fm;

The geometry of the embedded boundary is defined as two cylinders of radii 0.5 and 0.25.

  vertex scalar φ[];
  foreach_vertex()
    φ[] = difference (sq(0.5) - sq(x) - sq(y),
			sq(0.25) - sq(x) - sq(y));
  fractions (φ, cs, fs);

The outer cylinder is fixed and the inner cylinder is rotating with an angular velocity unity.

  u.n[embed] = dirichlet_embed (x*x + y*y > 0.14 ? 0. : - y);
  u.t[embed] = dirichlet_embed (x*x + y*y > 0.14 ? 0. :   x);

We initialize the reference velocity field.

  foreach()
    un[] = u.y[];
}

We look for a stationary solution.

event logfile (t += 0.01; i <= 1000) {
  double du = change (u.y, un);
  if (i > 0 && du < 1e-6)
    return 1; /* stop */
}

We compute error norms and display the angular velocity, pressure and error fields using bview.

#define powerlaw(r,N) (r*(pow(0.5/r, 2./N) - 1.)/(pow(0.5/0.25, 2./N) - 1.))

event profile (t = end)
{
  scalar utheta[], e[];
  foreach() {
    double θ = atan2(y, x), r = sqrt(x*x + y*y);
    if (cs[] > 0.) {
      utheta[] = - sin(θ)*u.x[] + cos(θ)*u.y[];
      e[] = utheta[] - powerlaw (r, 1.);
    }
    else
      e[] = p[] = utheta[] = nodata;
  }

  norm n = normf (e);
  fprintf (ferr, "%d %.3g %.3g %.3g %d %d %d %d %d\n",
	   N, n.avg, n.rms, n.max, i, mgp.i, mgp.nrelax, mgu.i, mgu.nrelax);
  dump();
  
  draw_vof ("cs", "fs", filled = -1, fc = {1,1,1});
  squares ("utheta", spread = -1);
  save ("utheta.png");

  draw_vof ("cs", "fs", filled = -1, fc = {1,1,1});
  squares ("p", spread = -1);
  save ("p.png");

  draw_vof ("cs", "fs", filled = -1, fc = {1,1,1});
  squares ("e", spread = -1);
  save ("e.png");

  if (N == 32)
    foreach() {
      double θ = atan2(y, x), r = sqrt(x*x + y*y);
      fprintf (stdout, "%g %g %g %g %g %g %g\n",
	       r, θ, u.x[], u.y[], p[], utheta[], e[]);
    }
}

Results

Angular velocity

Angular velocity

Pressure field

Pressure field

Error field

Error field

Velocity profile (N = 32)

Velocity profile (N = 32)

Convergence is close to second-order.

Error convergence

Error convergence

See also