cyl2d.c

Flow around a cylinder in 2D

Flow around a cylinder in 2D

Following Stephane’s setups for the Von Karman-street example and the starting flow around a cylinder test, we also setup a similar scenario.

We use automatically scaling refinement near the embedded surface. This maybe useful when multiple cylinders at various radii are used.

Vorticity and tracer particles

set xlabel 't'
set ylabel 'C_d'
plot 'log10-10' u 2:3

(script)

#include "embed.h"
#include "navier-stokes/centered.h"
#include "navier-stokes/double-projection.h"
#include "view.h"
#define BVIEW 1
#include "../particles.h"
#define CYLINDER (sqrt(sq(x) + sq(y - 0.01*R)) - R)

double R = 1, U = 1, Re = 2000, c = 10., C = 1.25, ue, nu;
int maxlevel = 10;

u.n[left]  = dirichlet (U);
uf.n[left] = U;
u.t[left]  = dirichlet (0);
p[left]    = dirichlet (0.);
pf[left]   = dirichlet (0.);

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

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

FILE * fp;
face vector muc[];
int main() {
  nu = U*R/Re;
  periodic (top);
  L0 = 30;
  X0 = -L0/3.5;
  Y0 = Z0 = -L0/2.;
  mu = muc;
  NITERMIN = 2;
  TOLERANCE = 1e-4;
  char logname[99];
  ue = U/c;
  sprintf (logname, "log%d-%g", maxlevel, c);
  fp = fopen (logname, "w");
  run();
}

event init (t = 0) {
  refine (CYLINDER < 0.2*R && CYLINDER > -0.2*R && level < maxlevel);
  vertex scalar phi[];
  foreach_vertex()
    phi[] = CYLINDER;
  boundary ({phi});
  fractions (phi, cs, fs);
  init_particles_2D_square_grid (10, -5, 0.01, 2);
  foreach()
    u.x[] = cs[] > 0;
}

In order to omit issues with inconsitent boundary conditions, we use a damping layer near the in and out flow boundaries.

event damp (i++) {
  coord Uinf = {U, 0, 0};
  foreach() {
    if (fabs(x - (X0 + L0/2.)) > 4*L0/10.) 
      foreach_dimension()
	u.x[] += dt*(Uinf.x - u.x[])/2.;
  }
  boundary ((scalar*){u});
}

event properties (i++) {
  foreach_face()
    muc.x[] = fm.x[]*nu; 
  boundary ((scalar*){muc});
}

void prolongate_ratio (Point point, scalar s) {
  foreach_child() {
    if (s[] != nodata)
      s[] += s[]*Delta;
  }
}
event adapt (i++) {
  scalar res[];
  foreach() {
    res[] = nodata;
    if (cs[] > 0 && cs[] < 1) {
      res[] = U/sqrt(R*nu);
    }
  }
  res.prolongation = prolongate_ratio;
  adapt_wavelet ((scalar*){res, u}, (double[]){C, ue, ue, ue}, maxlevel, 5);
  unrefine (level > 5 && (x - X0) < L0/10. && (X0 + L0 - x) < L0/10.);
}

event logger (i += 5) {
  coord Fp, Fmu;
  embed_force (p, u, mu, &Fp, &Fmu);
  fprintf (fp, "%d %g %g %g %g %g %ld\n",
	   i, t, Fp.x, Fp.y, Fmu.x, Fmu.y, grid->n);
}

event movies (t += 0.4) {
  view (fov = 6, width = 1200, height = 400, tx = -0.25);
  scalar omega[];
  vorticity (u, omega);
  boundary ({omega});
  translate (z = 0.05) {
    draw_vof ("cs", "fs", filled = -1, fc = {1,1,1});
    draw_vof ("cs", "fs", lw = 2);
  }
  squares ("omega", min = -2, max = 2, linear = true, map = cool_warm);
  cells();
  char str[99];
  sprintf (str, "Re = %g, C = %g, ML = %d", Re, c, maxlevel);
  draw_string (str, 1, lw = 3, lc = {1, 0, 1});
  scatter (loc);
  save ("mov.mp4");
}

event stop (t = 150) 
  fclose (fp);