sandbox/popinet/movingcylinder.c

An example of moving solid

This is a fairly rough way of approximating boundary conditions on a moving (or stationary) solid.

A cylinder moves through a liquid with a velocity corresponding to a Reynolds number of 160 i.e. this is equivalent in principle to the von Karman vortex street example but in a different reference frame.

#include "navier-stokes/centered.h"
#include "fractions.h"

int main(){
  L0 = 8.;
  origin (-0.5, -L0/2.);
  N = 512;
  const face vector muc[] = {0.00078125,0.00078125};
  mu = muc;
  run();
}

We just make an (empty) long channel.

event init (t = 0) {
  mask(y > 0.5 ? top: y < -0.5 ? bottom : none);
}

The moving cylinder is created by first initialising a volume fraction field at each timestep corresponding to the position and shape of the moving cylinder. Note that cylinder should really be a temporary field inside moving_cylinder(). We just keep it as a global field to be able to display it in the movie below.

scalar cylinder[];

event moving_cylinder (i++) {
  coord vc = {1.,0.}; // the velocity of the cylinder
  fraction (cylinder, - (sq(x - vc.x*t) + sq(y - vc.y*t) - sq(0.0625)));

We then use this (solid) volume fraction field to impose the corresponding velocity field inside the solid (as a volume-weighted average).

  foreach()
    foreach_dimension()
      u.x[] = cylinder[]*vc.x + (1. - cylinder[])*u.x[];
  boundary ((scalar *){u});
}

We output some convergence statistics…

event logfile (i++)
  fprintf (stderr, "%d %g %d %d %d %d\n", i, t,
	   mgp.i, mgp.nrelax, mgu.i, mgu.nrelax);

… and a movie of vorticity.

Animation of the vorticity field.

event images (t += 0.1; t <= 10.0) {
  static FILE * fp = popen ("ppm2gif > vort.gif", "w");
  scalar omega[];
  vorticity (u, omega);

Cells for which m is negative will be black in the movie.

  scalar m[];
  foreach()
    m[] = 0.5 - cylinder[];
  boundary ({m});

  output_ppm (omega, fp, box = {{-0.5,-0.5},{7.5,0.5}}, mask = m,
	      min=-10, max=10, linear=true);
}

We adapt according to the error on the velocity field as usual.

event adapt (i++) {
  adapt_wavelet ((scalar *){u}, (double[]){3e-2,3e-2}, 9, 4);
}

See also

This method was used in Gerris in particular for these papers: