# Starting flow around a cylinder

This is a canonical case of complex boundary layer separation, inspired by the experiments of Bouard & Coutanceau, 1980. Notable early numerical simulations include the results of Koumoutsakos and Leonard, 1995, hereafter K & L, which will be used in the comparisons below.

We will solve the Navier–Stokes equations and add the cylinder using an embedded boundary.

The “double projection” method is necessary to avoid noise in the pressure field which would pollute the drag force results.

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

High-resolution is needed to resolve the boundary layers properly. Mohaghegh et al., 2017 propose to use a maximum resolution of order D/10/\sqrt{Re}, with Re the Reynolds number and D the cylinder diameter. The cylinder diameter will be set to unity, and the domain size to 18, so that the corresponding levels of refinement are approximately 12 and 16 for Re=1000 and Re=9500, respectively.

int maxlevel = 12;  // 15/16 for Re = 9500, 12 for Re = 1000
double Re = 1000;   // or 9500
double cmax = 3e-3; // 1e-3 for Re = 9500, 3e-3 for Re = 1000

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

face vector muv[];

We set a unit velocity inflow on the left and an outflow on the right.

u.n[left] = dirichlet(1);

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

Command line arguments can be used to change the default parameters.

int main (int argc, char * argv[])
{
if (argc > 1)
maxlevel = atoi (argv);
if (argc > 2)
Re = atof (argv);
if (argc > 3)
cmax = atof (argv);

The domain is 18\times 18 and only half the cylinder is modelled.

  size (18);
origin (- L0/2.);

We set the viscosity field and tune the Poisson solver.

  mu = muv;
TOLERANCE = 1e-4;
NITERMIN = 2;

run();
}

The viscosity field needs to take the embedded boundary metric into account. Given that the diameter and velocity are set to one, the viscosity is just 1/Re.

event properties (i++) {
foreach_face()
muv.x[] = fm.x[]/Re;
}

## Initial conditions

event init (t = 0)
{

We can restart from a dump file.

  if (!restore ("restart")) {

Otherwise, we first create a mesh initially refined only around the cylinder.

    refine (level <= maxlevel*(1. - sqrt(fabs(sq(x) + sq(y) - sq(0.5)))/2.));

Then initialize the embedded boundary with the unit-diameter cylinder.

    solid (cs, fs, sq(x) + sq(y) - sq(0.5));

foreach()
u.x[] = cs[]; // fixme: with 1 this results in sub-optimal adaptation
}
else { // restart

When we restart, we still need to restore the face fraction field fs, since it is not dumped.

    solid (cs, fs, sq(x) + sq(y) - sq(0.5));

}

The boundary condition is zero velocity on the embedded boundary.

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

## Positions of the separation points

We would like to track the positions with time of the separation points on the surface of the cylinder, as done by K & L, 1995.

We first define a function which computes the vorticity at the surface of the cylinder and returns an array of (\theta,\omega) pairs, with \theta the angular coordinate and \omega the corresponding value of vorticity.

typedef struct {
double theta, omega;
} ThetaOmega;

int compar_theta (const void * a, const void * b)
{
const ThetaOmega * p1 = a, * p2 = b;
return p1->theta > p2->theta ? 1 : -1;
}

ThetaOmega * theta_omega()
{
Array * a = array_new();
foreach (serial)
if (cs[] > 0. && cs[] < 1.) {
coord n, b;
embed_geometry (point, &b, &n);
x += b.x*Delta, y += b.y*Delta;

ThetaOmega t;
t.omega = embed_vorticity (point, u, b, n);
t.theta = atan2(y, x);
array_append (a, &t, sizeof (ThetaOmega));
}
qsort (a->p, a->len/sizeof(ThetaOmega), sizeof(ThetaOmega), compar_theta);
ThetaOmega t = {nodata, nodata};
array_append (a, &t, sizeof (ThetaOmega));
ThetaOmega * p = a->p;
free (a);
return p;
}

The zeros of the function approximated by the (\theta,\omega) array are then recorded, together with the corresponding time, in the file pointed to by fp.

void omega_zero (FILE * fp)
{
// fixme: this function will not work with MPI
ThetaOmega * a = theta_omega();
for (ThetaOmega * o = a; (o + 1)->theta != nodata; o++) {
ThetaOmega * o1 = o + 1;
if (o1->omega*o->omega < 0.)
fprintf (fp, "%g %g\n", t,
o->theta + o->omega*(o1->theta - o->theta)/
(o->omega - o1->omega));
}
free (a);
fflush (fp);
}

## Log files

The log file contains the time, timestep, pressure and viscous forces components exerted by the fluid on the cylinder.

The “omega” file contains the locations where the surface vorticity vanishes.

event logfile (i += 10)
{
coord Fp, Fmu;
embed_force (p, u, mu, &Fp, &Fmu);
fprintf (stderr, "%d %g %g %.6f %.6f %.6f %.6f\n",
i, t, dt, Fp.x, Fp.y, Fmu.x, Fmu.y);

static FILE * fp = fopen ("omega", "w");
omega_zero (fp);
}

## Images and animations

We display the vorticity field and the corresponding adaptive mesh.

void display_omega (int width, int height)
{
view (fov = 1.68, tx = -0.0251023,
ty = 1e-12, // fixme: this is necessary to re-center the view
width = width, height = height);
draw_vof ("cs", filled = -1, fc = {1,1,1});
draw_vof ("cs", lw = 6);
squares ("omega", min = -12, max = 12, linear = true, map = blue_white_red);
mirror ({0,1}) {
draw_vof ("cs", lw = 6);
squares (color = "level", min = 6, max = maxlevel);
}
}

Still images are created at times matching those in various papers.

event snapshots (t = 1.; t <= 3.; t += 1.)
{
scalar omega[];
vorticity (u, omega);

if (t == 1. || t == 3.) {
// Fig. 4.
view (fov = 0.470916, tx = -0.0197934, ty = -0.0258659,
width = 524, height = 480);
draw_vof ("cs", filled = -1, fc = {1,1,1});
draw_vof ("cs", filled = 0, lw = 1);
double max = Re == 9500 ? 40 : 12;
squares ("omega", min = -max, max = +max, linear = 1, map = blue_white_red);
char name;
sprintf (name, "zoom-%g.png", t);
save (name);

draw_vof ("cs", lw = 2);
squares ("level");
sprintf (name, "cells-%g.png", t);
save (name);
}

if (t == 3.) {
// Fig. 3.
display_omega (640, 480);
save ("omega-3.png");
}

p.nodump = false;
char name ;
sprintf (name, "dump-%g", t);
dump (name);
}

An animation is created. The iteration interval is adjusted depending on the maximum spatial resolution.

event movie (i += 10*(1 << (maxlevel - 12)))
{
scalar omega[];
vorticity (u, omega);
display_omega (1280, 960);
save ("omega.mp4");
}

## Surface vorticity profiles

We are also interested in the details on the surface of the cylinder, in particular surface vorticity.

void cpout (FILE * fp)
{
foreach (serial)
if (cs[] > 0. && cs[] < 1.) {
coord b, n;
double area = embed_geometry (point, &b, &n);
x += b.x*Delta, y += b.y*Delta;

fprintf (fp, "%g %g %g %g %g %g\n",
x, // 1
y, // 2
atan2(y, x), // 3
embed_interpolate (point, p, b), // 4
area*Delta, // 5
embed_vorticity (point, u, b, n)); // 6
}
}

event surface_profiles (t = {0.5,0.9,1.5,2.5})
{
char name;
sprintf (name, "cp-%g", t);
FILE * fp = fopen (name, "w");
cpout (fp);
fclose (fp);
}

#if 0
event snapshot (i += 10) {
scalar omega[];
vorticity (u, omega);
p.nodump = false;
dump();

FILE * fp = fopen ("cp", "w");
cpout (fp);
fclose (fp);
}
#endif

The mesh is adapted according to the embedded boundary and velocity field.

event adapt (i++) {
}

## Results

### Re = 1000

Animation of the vorticity field and adaptive mesh, Re = 1000.

The final state at tU/D = 3 can be compared with figure 3 (top row) of Mohahegh et al. 2017.

set xlabel 'tU/D'
set ylabel 'C_D'
set grid
set pointsize 0.5
plot [][0:2]\
'fig1a.SIM' u ($1/2.):2 w l lw 2 t 'SIM (Mohahegh et al 2017)', \ 'fig1a.KL' u ($1/2.):2 pt 7 t 'K and L. 1995', \
'fig19.f' u ($1/2.):2 pt 9 t 'friction, K and L. 1995', \ 'fig19.p' u ($1/2.):2 pt 11 t 'pressure, K and L. 1995', \
'log' u 2:(4.*($4+$6)) w l lw 2 t 'Basilisk (12 levels)', \
'level-13/log' u 2:(4.*($4+$6)) w l lw 2 t 'Basilisk (13 levels)', \
'' u 2:(4.*$6) w l lw 2 t 'friction (13 levels)', \ '' u 2:(4.*$4) w l lw 2 t 'pressure (13 levels)' Time history of drag coefficient. Re = 1000. See Fig. 1a of Mohahegh et al. 2017. (script)

Note that the points of Figure 4 of K. & L. 1995 do not seem to match the data in Fig. 5a and 5b of the same paper, which explains the disagreement in the figure below. This agreement should be much better as can be seen on the more detailed surface vorticity figures below.

reset
set xlabel 'tU/D'
set ylabel 'Angle/pi'
set grid
set key bottom right
set pointsize 0.5
plot 'fig4.1000.KL' u ($1/2.):2 pt 7 t 'K and L. 1995', \ 'level-13/omega' u 1:($2/pi) pt 5 t 'Basilisk (13 levels)', \
'omega' u 1:($2/pi) pt 5 t 'Basilisk (12 levels)' Location of the points of zero surface vorticity. Re = 1000. See Fig. 4 of K. & L. 1995. (script) reset set xlabel 'theta/pi' set ylabel 'omega_sD/U' set grid plot [0:1]'fig5a.SIM' every 20 pt 7 t 'SIM (Mohahegh et al 2017)', \ 'fig5a.KL' w l lw 2 t 'K and L 1995', \ '< sort -k3,4 level-13/cp-0.5 | awk -f surface.awk' w l lw 2 \ t 'Basilisk (13 levels)', \ '< sort -k3,4 cp-0.5 | awk -f surface.awk' w l lw 2 \ t 'Basilisk (12 levels)' Surface vorticity at tU/D=0.5. Re = 1000. See Fig. 5a of Mohahegh et al. 2017. (script) plot [0:1]'fig5b.SIM' every 20 pt 7 t 'SIM (Mohahegh et al 2017)', \ 'fig5b.KL' w l lw 2 t 'K and L 1995', \ '< sort -k3,4 level-13/cp-1.5 | awk -f surface.awk' w l lw 2 \ t 'Basilisk (13 levels)', \ '< sort -k3,4 cp-1.5 | awk -f surface.awk' w l lw 2 \ t 'Basilisk (12 levels)' Surface vorticity at tU/D=1.5. Re = 1000. See Fig. 5b of Mohahegh et al. 2017. (script) ### Re = 9500 Animation of the vorticity field and adaptive mesh, Re = 9500, 16 levels. The final state at tU/D = 3 can be compared with figure 3 (bottom row) of Mohahegh et al. 2017. set xlabel 'tU/D' set ylabel 'C_D' set grid plot [][0:2.5] \ 'fig1b.SIM' u ($1/2.):2 w l lw 2 t 'SIM (Mohahegh et al 2017)', \
'fig1b.KL' u ($1/2.):2 pt 7 t 'K and L. 1995', \ 'level-15/log' u 2:(4.*($4+$6)) w l lw 2 t 'Basilisk (15 levels)', \ 'level-16/log' u 2:(4.*($4+$6)) w l lw 2 t 'Basilisk (16 levels)' Time history of drag coefficient. Re = 9500. See Fig. 1b of Mohahegh et al. 2017. (script) reset set term pngcairo enhanced font ",10" set output 'loc9500.png' set xlabel 'tU/D' set ylabel 'Angle/pi' set grid set key bottom right plot 'level-15/omega' u 1:($2/pi) pt 7 ps 0.25 t 'Basilisk (15 levels)', \
'level-16/omega' u 1:(\$2/pi) pt 5 ps 0.25 t 'Basilisk (16 levels)' Location of the points of zero surface vorticity. Re = 9500. See also Fig. 4 of K. & L. 1995. (script)

reset
set term @SVG
set xlabel 'theta/pi'
set ylabel 'omega_sD/U'
set grid
plot [0:1]'fig6a.SIM' every 20 pt 7 t 'SIM (Mohahegh et al 2017)', \
'fig6a.KL' w l lw 2 t 'K and L 1995', \
'< sort -k3,4 level-15/cp-0.9 | awk -f surface.awk' w l lw 2 \
t 'Basilisk (15 levels)', \
'< sort -k3,4 level-16/cp-0.9 | awk -f surface.awk' w l lw 2 \
t 'Basilisk (16 levels)' Surface vorticity at tU/D=0.9. Re = 9500. See Fig. 6a of Mohahegh et al. 2017. (script)

plot [0:1]'fig6b.SIM' every 20 pt 7 t 'SIM (Mohahegh et al 2017)', \
'fig6b.KL' w l lw 2 t 'K and L 1995', \
'< sort -k3,4 level-15/cp-2.5 | awk -f surface.awk' w l lw 2 \
t 'Basilisk (15 levels)', \
'< sort -k3,4 level-16/cp-2.5 | awk -f surface.awk' w l lw 2 \
t 'Basilisk (16 levels)' Surface vorticity at tU/D=2.5. Re = 9500. See Fig. 6b of Mohahegh et al. 2017. (script)

 [mohaghegh2017] Fazlolah Mohaghegh and HS Udaykumar. Comparison of sharp and smoothed interface methods for simulation of particulate flows ii: Inertial and added mass effects. Computers & Fluids, 143:103–119, 2017. [koumoutsakos1995] Petros Koumoutsakos and A Leonard. High-resolution simulations of the flow around an impulsively started cylinder using vortex methods. Journal of Fluid Mechanics, 296:1–38, 1995. [bouard1980] Roger Bouard and Madeleine Coutanceau. The early stage of development of the wake behind an impulsively started cylinder for 40 < Re < 104. Journal of Fluid Mechanics, 101(3):583–607, 1980.