sandbox/cselcuk/starting-dlmfd.c
Starting flow around a cylinder with DLMFD
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 we will be using the fictititous-domain method.
Note that we could also use the “double projection” method that is apparently necessary for the cut-cell method but we dont.
#include "DLMFD_reverse_Uzawa.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
double deltau;
scalar un[];
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[1]);
if (argc > 2)
Re = atof (argv[2]);
if (argc > 3)
cmax = atof (argv[3]);
The domain is 18\times 18 and we model the full cylinder.
size (18);
We tune the Poisson solver.
TOLERANCE = 1e-4;
NITERMIN = 2;
run();
}
Initial conditions
event init (t = 0)
{
/* Set dynamic viscosity */
const face vector muc[] = {1/Re, 1/Re};
mu = muc;
/* set density of the flow */
const scalar rhoc[] = 1.;
rho = rhoc;
/* initial condition: particles position */
particle * pp = particles;
GeomParameter gp = {0};
pp[0].iscube = 0;
pp[0].iswall = 0;
/* write the position with respect to 0 */
gp.center.x = L0/2;
gp.center.y = L0/2;
gp.radius = 0.5;
pp[0].g = gp;
/* particle id */
pp[0].pnum = 0;
We can restart from a dump file.
if (!restore ("restart")) {
init_file_pointers(pdata, fdata, 0); // 0 - initial simulation / 1 - restarting simulation
Otherwise, we first create a mesh initially refined only around the cylinder.
refine (sqrt (sq(x-gp.center.x) + sq(y-gp.center.y)) < 1.2*gp.radius && sqrt (sq(x-gp.center.x) + sq(y-gp.center.y)) > 0.8*gp.radius && level < maxlevel);
}
else { // restart
init_file_pointers(pdata, fdata, 1);
}
foreach()
un[] = u.x[];
}
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 (and interpolates) 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;
}
void theta_omega_dlmfd (const particle * pp, ThetaOmega ** to)
{
scalar vortz[];
vorticity (u, vortz);
Array * a = array_new();
coord center = (pp->g).center;
for (int ii = 0; ii < pp->s.m; ii++) {
ThetaOmega t;
t.omega = interpolate (vortz, pp->s.x[ii], pp->s.y[ii]);
t.theta = atan2 (pp->s.y[ii] - center.y, pp->s.x[ii] - center.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);
*to = 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_dlmfd (FILE * fp, const particle * pp)
{
/* fixme: this function will not work with MPI because of the
interpolations which are thread-dependent */
ThetaOmega * a = NULL;
theta_omega_dlmfd (pp, &a);
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);
}
We log the number of iterations of the multigrid solver for pressure and viscosity.
event logfile (i++) {
deltau = change (u.x, un);
fprintf (stderr, "log output %d %g %d %d %g %g %g %ld\n", i, t, mgp.i, mgu.i, mgp.resa, mgu.resa, deltau, grid->tn);
}
Images and animations
We display the vorticity field and the corresponding adaptive mesh.
void display_omega (int width, int height)
{
view (fov = 2.16967, quat = {-0.00418115,0.0542593,-0.000125816,0.998518}, tx = -0.53112, ty = -0.499909, bg = {1,1,1}, width = width, height = height, samples = 1);
squares ("omega", min = -12, max = 12, linear = 1, map = cool_warm);
}
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 = 2.16967, quat = {-0.00418115,0.0542593,-0.000125816,0.998518}, tx = -0.53112, ty = -0.499909, bg = {1,1,1}, width = 966, height = 775, samples = 1);
double max = Re == 9500 ? 40 : 12;
squares ("omega", min = -max, max = +max, linear = 1, map = cool_warm);
char name[80];
sprintf (name, "zoom-%g.png", t);
save (name);
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 [80];
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");
}
event movie1 (i += 10*(1 << (maxlevel - 12)))
{
scalar omega[];
vorticity (u, omega);
display_omega (1280, 960);
cells();
save ("movie1.mp4");
}
Surface vorticity profiles
We are also interested in the details on the surface of the cylinder, in particular surface vorticity.
void cpout_dlmfd (FILE * fp, particle * pp)
{
scalar omega[];
vorticity(u, omega);
for (int ii = 0; ii < pp->s.m; ii++) {
GeomParameter gc = pp->g;
coord c = gc.center;
double w = atan2 (pp->s.y[ii] - c.y, pp->s.x[ii] - c.x);
if (w < 0) w += 2*pi;
fprintf (fp, "%g %g %g %g \n",
pp->s.x[ii], // 1
pp->s.y[ii], // 2
w, // 3
interpolate (omega, pp->s.x[ii], pp->s.y[ii])
);
}
}
Adaptive mesh refinement
The mesh is adapted according to the embedded boundary and velocity field.
event adapt (i++) {
particle * pp = particles;
double rhoval = 1.;
DLMFD_subproblem (pp, i, rhoval);
/* Save forces acting on particles before the adapting the mesh */
sumLambda (pp, fdata, t, dt, flagfield, DLM_lambda, index_lambda, rhoval);
static FILE * fp = fopen ("omega", "w");
omega_zero_dlmfd (fp, pp);
if (t == 0.5|| t == 0.9|| t ==1.5 || t == 2.5) {
printf("inside cp block\n");
char name[80];
sprintf (name, "cp-%g", t);
FILE * fp2 = fopen (name, "w");
cpout_dlmfd (fp2, pp);
fclose (fp2);
}
/* Free particle structures (we dont need them anymore) */
free_particles (pp, NPARTICLES);
/* Save particles trajectories */
particle_data(pp, t, i, pdata);
/* char name[80] = "test-init-starting"; */
/* scalar * list = {p, index_lambda.x, flagfield}; */
/* vector * vlist = {u}; */
/* save_data(list, vlist, i, t, name); */
/* return 1; */
adapt_wavelet ({flagfield,u}, (double[]){1e-4,cmax,cmax}, maxlevel, 5);
}
event surface_profiles (t = {0.5,0.9,1.5,2.5})
{
printf("forcing output data at t= {0.5,0.9,1.5,2.5}\n")
}
Results
Re = 1000
Animation of the vorticity field and adaptive mesh, Re = 1000.
Animation of the 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]\
"`echo $BASILISK`/test/starting/level-13/log" u 2:(4.*($4+$6)) w l lw 2 t 'Basilisk-cut-cells (13 levels)',\
"`echo $BASILISK`/test/starting/fig1a.SIM" u ($1/2.):2 w l lw 2 t 'SIM (Mohahegh et al 2017)', \
"`echo $BASILISK`/test/starting/fig1a.KL" u ($1/2.):2 pt 7 t 'K and L. 1995', \
"`echo $BASILISK`/test/starting/fig19.f" u ($1/2.):2 pt 9 t 'friction, K and L. 1995', \
"`echo $BASILISK`/test/starting/fig19.p" u ($1/2.):2 pt 11 t 'pressure, K and L. 1995', \
'sum_lambda-0' every 10 u 1:($2*2) w l lw 2 t 'Basilisk-dlmfd (12 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 "`echo $BASILISK`/test/starting/level-13/omega" u 1:($2/pi) pt 5 t 'Basilisk-cut-cell (13 levels)', \
"`echo $BASILISK`/test/starting/fig4.1000.KL" u ($1/2.):2 pt 7 ps 0.9 t 'K and L. 1995', \
'omega' every 10 u 1:($2/pi) pt 5 t 'Basilisk-dlmfd (12 levels)'
Location of the points of zero surface vorticity. Re = 1000. See Fig. 4 of K. & L. 1995. (script)
Note that the vorticity field is not constrained within the ficititous domain problem and an additional interpolation is needed to get the vorticity values at the surface of the cylinder on the location of the Lagrange multipliers. This results in a loss of accuracy in the post-process.
reset
set xlabel 'theta/pi'
set ylabel 'omega_sD/U'
set grid
plot [0:2] '< sort -k3,4 "$BASILISK/test/starting/level-13/cp-0.5" |awk -f "$BASILISK/test/starting/surface.awk"' w l lw 2 t 'Basilisk-cut-cells (13 levels)',\
"`echo $BASILISK`/test/starting/fig5a.SIM" every 20 pt 7 t 'SIM (Mohahegh et al 2017)',\
"`echo $BASILISK`/test/starting/fig5a.KL" w l lw 2 t 'K and L 1995', \
'cp-0.5' every 10 u ($3/pi):($4*2) w l lw 2 \
t 'Basilisk-dlmfd (12 levels)'
Surface vorticity at tU/D=0.5. Re = 1000. See Fig. 5a of Mohahegh et al. 2017. (script)
plot [0:2] '< sort -k3,4 "$BASILISK/test/starting/level-13/cp-1.5" |awk -f "$BASILISK/test/starting/surface.awk"' w l lw 2 t 'Basilisk-cut-cells (13 levels)',\
"`echo $BASILISK`/test/starting/fig5b.SIM" every 20 pt 7 t 'SIM (Mohahegh et al 2017)', \
"`echo $BASILISK`/test/starting/fig5b.KL" w l lw 2 t 'K and L 1995', \
'cp-1.5' every 10 u ($3/pi):($4*2) w l lw 2 \
t 'Basilisk-dlmfd (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, 15 levels.
Animation of the adaptive mesh, Re = 9500.
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] \
"`echo $BASILISK`/test/starting/fig1b.SIM" u ($1/2.):2 w l lw 2 t 'SIM (Mohahegh et al 2017)', \
"`echo $BASILISK`/test/starting/fig1b.KL" u ($1/2.):2 pt 7 t 'K and L. 1995', \
"`echo $BASILISK`/test/starting/level-15/log" u 2:(4.*($4+$6)) w l lw 2 t 'Basilisk (15 levels)', \
"`echo $BASILISK`/test/starting/level-16/log" u 2:(4.*($4+$6)) w l lw 2 t 'Basilisk (16 levels)',\
'level-15/sum_lambda-0' every 10 u 1:($2*2) w l lw 2 t 'Basilisk-dlmfd (15 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 "`echo $BASILISK`/test/starting/level-15/omega" u 1:($2/pi) pt 7 ps 0.25 t 'Basilisk (15 levels)', \
"`echo $BASILISK`/test/starting/level-16/omega" u 1:($2/pi) pt 5 ps 0.25 t 'Basilisk (16 levels)',\
'level-15/omega' every 10 u 1:($2/pi) pt 5 t 'Basilisk-dlmfd (15 levels)'
![Location of the points of zero surface vorticity. Re = 9500. See also Fig. 4 of K. & L. 1995. (script)](starting-dlmfd/loc9500.png?1710574983)
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:2] "`echo $BASILISK`/test/starting/fig6a.SIM" every 20 pt 7 t 'SIM (Mohahegh et al 2017)', \
"`echo $BASILISK`/test/starting/fig6a.KL" w l lw 2 t 'K and L 1995', \
'< sort -k3,4 "$BASILISK/test/starting/level-15/cp-0.9" | awk -f "$BASILISK/test/starting/surface.awk"' w l lw 2 \
t 'Basilisk (15 levels)', \
'< sort -k3,4 "$BASILISK/test/starting/level-16/cp-0.9" | awk -f "$BASILISK/test/starting/surface.awk"' w l lw 2 \
t 'Basilisk (16 levels)',\
'level-15/cp-0.9' every 10 u ($3/pi):($4*2) w l lw 2 \
t 'Basilisk-dlmfd (15 levels)'
Surface vorticity at tU/D=0.9. Re = 9500. See Fig. 6a of Mohahegh et al. 2017. (script)
plot [0:2] "`echo $BASILISK`/test/starting/fig6b.SIM" every 20 pt 7 t 'SIM (Mohahegh et al 2017)', \
"`echo $BASILISK`/test/starting/fig6b.KL" w l lw 2 t 'K and L 1995', \
'< sort -k3,4 "$BASILISK/test/starting/level-15/cp-2.5" | awk -f "$BASILISK/test/starting/surface.awk"' w l lw 2 \
t 'Basilisk (15 levels)', \
'< sort -k3,4 "$BASILISK/test/starting/level-16/cp-2.5" | awk -f "$BASILISK/test/starting/surface.awk"' w l lw 2 \
t 'Basilisk (16 levels)',\
'level-15/cp-2.5' every 10 u ($3/pi):($4*2) w l lw 2 \
t 'Basilisk-dlmfd (15 levels)'
Surface vorticity at tU/D=2.5. Re = 9500. See Fig. 6b of Mohahegh et al. 2017. (script)
References
[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. |