src/test/poiseuille45.c

Poiseuille flow in a periodic channel inclined at 45 degrees

We test the embedded boundaries by solving the viscous flow driven by gravity in an inclined periodic channel.

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

int main()
{
  origin (-0.5, -0.5);
  periodic (right);
  periodic (top);
  
  stokes = true;
  TOLERANCE = 1e-7;
  
  for (N = 16; N <= 64; N *= 2)
    run();
}

scalar un[];

#define WIDTH 0.5
#define EPS 1e-14

event init (t = 0) {

The gravity vector is aligned with the channel and viscosity is unity.

  const face vector g[] = {1.,1.};
  a = g;
  μ = fm;

The channel geometry is defined using the levelset function ϕ.

  vertex scalar φ[];
  foreach_vertex()
    φ[] = difference (union (y - x - EPS, x - y - 0.5 + EPS),
			y - x - 0.5 + EPS);
  fractions (φ, cs, fs);

The boundary condition is zero velocity on the embedded boundaries.

  u.n[embed] = dirichlet_embed(0);
  u.t[embed] = dirichlet_embed(0);

We initialize the reference velocity.

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

We check for a stationary solution.

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

We compute the error and display the solution using bview.

event profile (t = end) {
  printf ("\n");
  foreach()
    fprintf (stdout, "%g %g %g %g %g\n", x, y, u.x[], u.y[], p[]);
  scalar e[];
  foreach() {
    double x1 = y - x;
    e[] = u.x[] - 0.25*(sq(WIDTH/2.) - sq(x1 >= 0. ? x1 - 0.25 : x1 + 0.75));
  }
  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);
  
  draw_vof ("cs", "fs");
  squares ("u.x", linear = true, spread = -1);
  save ("u.x.png");
}
Velocity field

Velocity field

The method is almost exact.

Error as function of resolution

Error as function of resolution