cylinder-confined.c

Flow past a fixed confined cylinder at Re = 20

Flow past a fixed confined cylinder at Re = 20

We reproduce here a benchmark solution proposed in Schafer and Turek. In this benchmark case, we study the flow around a cylinder positioned slightly off mid-stream. Instead of using a square domain, we use embedded boundaries to define a rectangular domain.

We solve here the Navier-Stokes equations and add the cylinder using an embedded boundary.

#include "grid/quadtree.h"
#include "../myembed.h"
#include "../mycentered.h"
#include "view.h"

Reference solution

#define d    (0.1) // Cylinder characteristic length
#if RE // Re = 100
#define uref (2./3.*1.5) // Reference velocity, uref
#define Re   (100.)
#else // Re = 20
#define uref (2./3.*0.3) // Reference velocity, uref
#define Re   (20.)
#endif // RE
#define tref ((d)/(uref)) // Reference time, tref=d/u
#define um   (1.5*(uref)) // Maximum velocity

We also define the shape of the domain.

#define ht  (2.1*(d)) // Top-half width of the channel
#define hb  (2.*(d)) // Bottom-half width of the channel
#define EPS (1.e-12)

#define cylinder(x,y) (sq ((x)) + sq ((y)) - sq ((d)/2.))
#define wall(y,w)     ((y) - (w)) // + over, - under

coord p_p;

void p_shape (scalar c, face vector f, coord p)
{
  vertex scalar phi[];
  foreach_vertex()
    phi[] = intersection (
			  (cylinder ((x - p.x), (y - p.y))),
			  intersection (
					-(wall (y,  0.5*((hb) + (ht)))),
					 (wall (y, -0.5*((hb) + (ht))))
					)
			  );
  boundary ({phi});
  fractions (phi, c, f);
  fractions_cleanup (c, f,
		     smin = 1.e-14, cmin = 1.e-14);
}

Setup

We need a field for viscosity so that the embedded boundary metric can be taken into account.

face vector muv[];

We also define a reference velocity field.

scalar un[];

Since both the walls of the channel and the cylinder are described using the same embedded volume fraction cs, we declare a color field p_col that we use to color the cylinder only.

scalar p_col[];

void p_shape_col (scalar c, coord p)
{
  vertex scalar phi[];
  foreach_vertex()
    phi[] = (cylinder ((x - p.x), (y - p.y)));
  boundary ({phi});
  fractions (phi, c);
}

We define the mesh adaptation parameters.

#define lmin (6)  // Min mesh refinement level (l=6 is 3pt/d)
#define lmax (11) // Max mesh refinement level (l=11 is 92pt/d)
#define cmax (5.e-3*(uref)) // Absolute refinement criteria for the velocity field

int main ()
{

The domain is 22 d\times 22 d.

  L0 = 22.*(d);
  size (L0);
  origin (0., -L0/2.);

We set the maximum timestep.

  DT = 1.e-2*(tref);

We set the tolerance of the Poisson solver.

  TOLERANCE    = 1.e-6;
  TOLERANCE_MU = 1.e-6*(uref);

We initialize the grid.

  N = 1 << (lmin);
  init_grid (N);
  
  run();
}

Boundary conditions

We use inlet boundary conditions with a parabolic velocity profile.

#define profile(y) ((y) <= 0.5*((hb) + (ht)) && (y) >= -0.5*((hb) + (ht)) ? \
		    4.*((y) + 0.5*((hb) + (ht)))/((hb) + (ht))*(1. - ((y) + 0.5*((hb) + (ht)))/((hb) + (ht))) : \
		    0.)

u.n[left] = dirichlet ((um)*(profile (y)));
u.t[left] = dirichlet (0);
p[left]   = neumann (0);
pf[left]  = neumann (0);

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

We give boundary conditions for the face velocity to “potentially” improve the convergence of the multigrid Poisson solver.

uf.n[left]   = (um)*(profile (y));
uf.n[bottom] = 0;
uf.n[top]    = 0;

Properties

event properties (i++)
{
  foreach_face()
    muv.x[] = (uref)*(d)/(Re)*fm.x[];
  boundary ((scalar *) {muv});
}

Initial conditions

event init (i = 0)
{

We set the viscosity field in the event properties.

  mu = muv;

We use “third-order” face flux interpolation.

#if ORDER2
  for (scalar s in {u, p, pf})
    s.third = false;
#else
  for (scalar s in {u, p, pf})
    s.third = true;
#endif // ORDER2

We use a slope-limiter to reduce the errors made in small-cells.

#if SLOPELIMITER
  for (scalar s in {u}) {
    s.gradient = minmod2;
  }
#endif // SLOPELIMITER
  
#if TREE

When using TREE and in the presence of embedded boundaries, we should also define the gradient of u at the cell center of cut-cells.

#endif // TREE

We initialize the embedded boundary.

We first define the cylinder’s position.

  p_p.x = 2.*(d) + (EPS);
  p_p.y = (hb) - 0.5*((hb) + (ht)) + (EPS);
  
#if TREE

When using TREE, we refine the mesh around the embedded boundary.

  astats ss;
  int ic = 0;
  do {
    ic++;
    p_shape (cs, fs, p_p);
    ss = adapt_wavelet ({cs}, (double[]) {1.e-30},
			maxlevel = (lmax), minlevel = (1));
  } while ((ss.nf || ss.nc) && ic < 100);
#endif // TREE
  
  p_shape (cs, fs, p_p);

We also define the volume fraction at the previous timestep csm1=cs.

  csm1 = cs;

We define the no-slip boundary conditions for the velocity.

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

  uf.n[embed] = dirichlet (0);
  uf.t[embed] = dirichlet (0);
  pf[embed]   = neumann (0);

We initialize the velocity.

  foreach() 
    u.x[] = cs[]*(um)*(profile (y));
  boundary ((scalar *) {u});

We finally initialize the reference velocity field.

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

Embedded boundaries

#if TREE
event adapt (i++)
{
  adapt_wavelet ({cs,u}, (double[]) {1.e-2,(cmax),(cmax)},
  		 maxlevel = (lmax), minlevel = (1));

We also reset the embedded fractions to avoid interpolation errors on the geometry. This would affect in particular the computation of the pressure contribution to the hydrodynamic forces.

  p_shape (cs, fs, p_p);
}
#endif // TREE

Outputs

double sdt = 0., CDm1 = 0.;
double CDavg = 0., CDmin = HUGE, CDmax = -HUGE;
double CLavg = 0., CLmin = HUGE, CLmax = -HUGE;
double La = 0., Lamin = HUGE, Lamax = -HUGE;
double DP = 0., DPmin = HUGE, DPmax = -HUGE;

event logfile (i++; t <= 200.*(tref))
{
  double du = change (u.x, un);
  
  coord Fp, Fmu;
  p_shape_col (p_col, p_p);
  boundary ({p_col});
  embed_color_force (p, u, mu, p_col, &Fp, &Fmu);

  double CD = (Fp.x + Fmu.x)/(0.5*sq ((uref))*(d));
  double CL = (Fp.y + Fmu.y)/(0.5*sq ((uref))*(d));

  double dF = fabs (CDm1 - CD);
  CDm1 = CD;

  fprintf (stderr, "%g %d %g %g %d %d %d %d %d %d %g %g %g %g %g %g %g %g\n",
	   (Re),
	   i, t/(tref), dt/(tref),
	   mgp.i, mgp.nrelax, mgp.minlevel,
	   mgu.i, mgu.nrelax, mgu.minlevel,
	   mgp.resb, mgp.resa,
	   mgu.resb, mgu.resa,
	   CD, CL,
	   du, dF);
  fflush (stderr);
  
  if ((du < 1.e-5 && dF < 1.e-4 && t > 100.*(tref)) || t > 170.*(tref)) {

We compute the average, minimun and maximum drag and lift coefficients after an initial transitory state.

    CDavg += CD*dt;
    CLavg += CL*dt;
    sdt   += dt;

    if (CD > CDmax)
      CDmax = CD;
    if (CD < CDmin)
      CDmin = CD;
    if (CL > CLmax)
      CLmax = CL;
    if (CL < CLmin)
      CLmin = CL;

We compute here the pressure drop between the front and back of the cylinder end evaluate the lenght of the recirculation area.

    double step = ((L0)/(1 << (lmax)));

Pressure drop.

    double pa = interpolate (p, ((p_p.x) - (d)/2) - step/2., (p_p.y));
    double pe = interpolate (p, ((p_p.x) + (d)/2) + step/2., (p_p.y));
    DP = pa - pe;

    if (DP > DPmax)
      DPmax = DP;
    if (DP < DPmin)
      DPmin = DP;

Recirculation length.

    La = -((p_p.x) + (d)/2);
    for (double x = ((p_p.x) + (d)/2) + step/2.; x <= (L0) - step/2.; x += step) {
      double o = interpolate (u.x, x, (p_p.y));
      if (o > 0.) {
	La += x;
	break;
      }
    }

    if (La > Lamax)
      Lamax = La;
    if (La < Lamin)
      Lamin = La;
  }
}

event coeffs (t = end)
{
  char name1[80];
  sprintf (name1, "avg-min-max.dat");
  static FILE * fp = fopen (name1, "w");
  fprintf (fp, "%g %d %g %g %g %g %g %g %g %g %g %g %g %g\n",
	   (Re), (1 << (lmax)),
	   CDavg, CDmin, CDmax,
	   CLavg, CLmin, CLmax,
	   La, Lamin, Lamax,
	   DP, DPmin, DPmax);
  fflush (fp);
  fclose (fp);
}

Snapshots

event snapshot (t = end)
{
  scalar omega[];
  vorticity (u, omega);
  
  char name2[80];

We first plot the entire domain.

  view (fov = 20, camera = "front",
	tx = -(L0)/2./(L0), ty = 0.,
	bg = {1,1,1},
	width = 800, height = 800);

  draw_vof ("cs", "fs", lw = 5);
  cells ();
  sprintf (name2, "mesh.png");
  save (name2);
    
  draw_vof ("cs", "fs", filled = -1, lw = 5);
  squares ("u.x", map = cool_warm);
  sprintf (name2, "ux.png");
  save (name2);

  draw_vof ("cs", "fs", filled = -1, lw = 5);
  squares ("u.y", map = cool_warm);
  sprintf (name2, "uy.png");
  save (name2);

  draw_vof ("cs", "fs", filled = -1, lw = 5);
  squares ("p", map = cool_warm);
  sprintf (name2, "p.png");
  save (name2);

  draw_vof ("cs", "fs", filled = -1, lw = 5);
  squares ("omega", map = cool_warm);
  sprintf (name2, "omega.png");
  save (name2);

We then zoom on the cylinder.

  view (fov = 3, camera = "front",
	tx = -(p_p.x + 4.*(d))/(L0), ty = 0.,
	bg = {1,1,1},
	width = 800, height = 250);

  draw_vof ("cs", "fs", lw = 5);
  cells ();
  sprintf (name2, "mesh-zoom.png");
  save (name2);

  draw_vof ("cs", "fs", lw = 5);
  squares ("u.x", map = cool_warm);
  sprintf (name2, "ux-zoom.png");
  save (name2);

  draw_vof ("cs", "fs", lw = 5);
  squares ("u.y", map = cool_warm);
  sprintf (name2, "uy-zoom.png");
  save (name2);

  draw_vof ("cs", "fs", lw = 5);
  squares ("p", map = cool_warm);
  sprintf (name2, "p-zoom.png");
  save (name2);

  draw_vof ("cs", "fs", lw = 5);
  squares ("omega", map = cool_warm);
  sprintf (name2, "omega-zoom.png");
  save (name2);
}

Results

Snapshots for level 11

drawing drawing
drawing drawing
drawing drawing
drawing drawing
drawing drawing

Drag and lift coefficients

We compare the steady values of the drag and lift coefficients C_D and C_L with those obtained in the benchmark study of Schafer and Turek.

reset
set terminal svg font ",16"
set key top right spacing 1.1
set grid ytics
set xtics 0,20,200
set xlabel 't/(d/u)'
set ylabel 'C_{D}'
set xrange [0:200]
set yrange [0:11]
plot 5.57 w p pt 1 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=20, C_{D,min}", \
     5.59 w p pt 2 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=20, C_{D,max}", \
     3.22 w p pt 3 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=100, C_{D,min}", \
     3.24 w p pt 4 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=100, C_{D,max}", \
     'log' u 3:15 w l lc rgb "blue" t "Basilisk, l=11"

Time evolution of the drag coefficient C_D (script)

set ylabel 'C_{L}'
set yrange [-2:4]
plot 0.0104 w p pt 1 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=20, C_{L,min}", \
     0.0110 w p pt 2 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=20, C_{L,max}", \
     0.99   w p pt 3 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=100, C_{L,min}", \
     1.01   w p pt 4 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=100, C_{L,max}", \
     'log' u 3:16 w l lc rgb "blue" t "Basilisk, l=11"

Time evolution of the lift coefficient C_L (script)

set xtics 64,2,4096
set xlabel 'N'
set ylabel 'C_{D}'
set xrange [64:4096]
set yrange [1:9]
set logscale x
plot 5.57 w lp pt 1 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=20, C_{D,min}", \
     5.59 w lp pt 2 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=20, C_{D,max}", \
     3.22 w lp pt 3 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=100, C_{D,min}", \
     3.24 w lp pt 4 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=100, C_{D,max}", \
     'avg-min-max.dat' u 2:3 w p pt 7  lc rgb "black" t "Basilisk, C_{D,avg}", \
     ''                u 2:4 w p pt 5  lc rgb "blue"  t "Basilisk, C_{D,min}", \
     ''                u 2:5 w p pt 2  lc rgb "red"   t "Basilisk, C_{D,max}"

Minimum and maximum drag coeffficient C_D (script)

set ylabel 'C_{L}'
set yrange [-2:4]
plot 0.0104 w lp pt 1 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=20, C_{L,min}", \
     0.0110 w lp pt 2 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=20, C_{L,max}", \
     0.99   w lp pt 3 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=100, C_{L,min}", \
     1.01   w lp pt 4 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=100, C_{L,max}", \
     'avg-min-max.dat' u 2:6 w p pt 7  lc rgb "black" t "Basilisk, C_{L,avg}", \
     ''                u 2:7 w p pt 5  lc rgb "blue"  t "Basilisk, C_{L,min}", \
     ''                u 2:8 w p pt 2  lc rgb "red"   t "Basilisk, C_{L,max}"

Minimum and maximum lift coeffficient C_L (script)

Pressure drop and recirculation area

We do the same for the pressure drop \Delta p and the length of the recirculation area L_a.

set ylabel 'La'
set yrange [0:0.25]
plot 0.0842 w lp pt 1 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=20, La_{min}", \
     0.0852 w lp pt 2 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=20, La_{max}", \
     'avg-min-max.dat' u 2:9  w p pt 7  lc rgb "black" t "Basilisk, La_{avg}", \
     ''                u 2:10 w p pt 5  lc rgb "blue"  t "Basilisk, La_{min}", \
     ''                u 2:11 w p pt 2  lc rgb "red"   t "Basilisk, La_{max}"

Recirculation area (script)

set ylabel '{D}p'
set yrange [-1:5]
plot 0.1172 w lp pt 1 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=20, {D}p_{min}", \
     0.1176 w lp pt 2 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=20, {D}p_{max}", \
     2.46   w lp pt 3 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=100, {D}p_{min}", \
     2.5    w lp pt 4 ps 0.6 lc rgb "black" t "Schafer et al., 1996, Re=100, {D}p_{max}", \
     'avg-min-max.dat' u 2:12 w p pt 7  lc rgb "black" t "Basilisk, DP_{avg}", \
     ''                u 2:13 w p pt 5  lc rgb "blue"  t "Basilisk, DP_{min}", \
     ''                u 2:14 w p pt 2  lc rgb "red"   t "Basilisk, DP_{max}"

Pressure drop (script)

References

[schafer1996]

M. Schafer and S. Turek. Benchmark computations of laminar flow around a cylinder. Flow Simulation with High-Performance Computers II, 1996.