/** # Vortex Ejection from a mode 3 instability According to [Kizner et al. (2013)](https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/instabilities-of-the-flow-around-a-cylinder-and-emission-of-vortex-dipoles/2327C2CFA76059D27461A2BC6A09F146), a flow around a no-slip cylinder with radius $R$ maybe unstable and eject three dipolar vortex pairs. We study the flow using embedded boundaries and the Navier-Stokes solver. Furthermore, we will view our results. */ #include "embed.h" #include "navier-stokes/centered.h" #include "view.h" /** The maximum resolution is set to $\Delta_{min}=R/100$. This allows to run on the sandbox server. */ int maxlevel = 12; double Re = 30000., ue = 0.003; face vector muc[]; int main(){ init_grid(64); L0 = 40; mu = muc; X0 = Y0 = -L0/2; run(); } event properties (i++){ foreach_face() muc.x[] = fm.x[]/Re; boundary ((scalar*){muc}); } /** The cylinder is defined and the flow field is initialized c.f. Kizner et al. with an $m=3$ perturbation. */ #define RAD (pow(sq(x) + sq(y), 0.5)) #define THETA(M) (M*asin(x/RAD)) #define RADP(P, M) ((1 + P*sin(THETA(M)))/(pow(1 + 0.5*sq(P), 0.5))) double a1 = 1.5, b1 = 2.25; double P = 0.005, m = 3; event init (t = 0){ double gamma = (sq(a1) - 1.)/(sq(b1) - sq(a1)); refine (RAD < b1 && level < (maxlevel - 1)); refine (RAD < 1.05 && RAD > 0.95 && level < maxlevel); vertex scalar phi[]; foreach_vertex() phi[] = RAD - 1; fractions (phi, cs, fs); foreach(){ double r = RAD; double r1 = RADP(P,m)*r; double vr; if (r1 > 0.9 && r1 < a1) vr = r1 - 1./r1; else if (r1 >= a1 && r1 <= b1) vr = -gamma*r1 + ((1 + gamma)*sq(a1) - 1.)/r1; else // (0.9 > r || r > b) vr = 0; u.x[] = cm[]*0.5*vr*r*-y/(sq(x) + sq(y)); u.y[] = cm[]*0.5*vr*r* x/(sq(x) + sq(y)); } /** The boundary conditions on the embedded boundary are set: */ u.t[embed] = dirichlet (0.); u.n[embed] = dirichlet (0.); boundary (all); } /** The grid is adaptedively refined and coarsened to properly represent the boundary and the evolution of the flow field. We set $\zeta_{u,v}\approx U_{max}/150$. */ event adapt (i++) adapt_wavelet ({cs, u.x, u.y}, (double[]){0.01, ue, ue}, maxlevel); /** ## Output Movies are generated displaying the vorticity field and the adaptive-grid-cell structure. ![](kizner2/kizner.mp4) ![](kizner2/kizner_cells.mp4) */ event movie (t += 0.5){ scalar omega[]; vorticity (u, omega); boundary ({omega}); view (fov = 8, width = 600, height = 600); clear(); draw_vof ("cs", filled = -1, fc = {1., 1., 1.}); squares ("omega", min = -0.75, max = 0.75, map = cool_warm, linear = true); draw_vof ("cs", "fs", lw = 2); save ("kizner.mp4"); clear(); view (fov = 5); cells(); save ("kizner_cells.mp4"); clear(); view (fov = 2, width = 1080, height = 1080); draw_vof ("cs", filled = -1, fc = {0.6, 0.6, 0.6}); squares ("omega", min = -0.75, max = 0.75, map = cool_warm); cells(); save ("kizner_cells_zoom.mp4"); } /** Also there is this movie: ![A zoom](kizner2/kizner_cells_zoom.mp4)(width="50%") We log some solver data: */ event logger (t += 0.1){ fprintf (stderr, "%d %g %d %d %d %g %g\n", i, t, mgp.i, mgu.i, mgpf.i, DT, dt); int cn = 0; foreach(reduction(+:cn)) cn++; if (pid() == 0){ static FILE * fp = fopen ("data", "w"); fprintf (fp, "%g\t%d\t%d\t%g\t%g\n", t, i, cn, dt, DT); fflush (fp); } } /** Including the number of grid cells: ~~~gnuplot These numbers may be compared against the millions of cells that Kizner et al. (2013) employed. set yr [0 : 25000] set xlabel 'time' set ylabel 'Cells' set key off plot 'data' u 1:3 w l lw 2 ~~~ */ event the_last_event (t = 100); /** ## Flow direction partitioning We can study the radially-averaged absolute velocity in the azimutal and radial direction as a function of time and the distance from the centre. Therefore, we render movies of the radial profiles using profile6.h, gnuplot and ffmpeg. ![A more quantative perspective on the instability. Mind the logaritmic vertical axis](kizner2/mov.mp4) */ #include "profile6.h" FILE * gnuplotPipe; event init (t = 0){ gnuplotPipe = popen ("gnuplot", "w"); fprintf(gnuplotPipe, "set term pngcairo\n" "set xr [1 : 6]\n" "set yr [0.001 : 1]\n" "set logscale y\n" "set key top right\n" "set grid\n" "set xlabel 'r/R'\n" "set ylabel 'E'\n" "set size ratio 0.5\n"); } int frame = 0; event profs(t += 0.5){ scalar ur[], ua[] ,* list; list = {ur, ua}; foreach(){ double r = sqrt(sq(x) + sq(y)); ur[] = fabs(u.x[]*x/r + u.y[]*y/r); ua[] = fabs((u.y[]*x/r) - (u.x[]*y/r)); } profile (list, sqrt(sq(x) + sq(y)), "prof.dat"); fprintf (gnuplotPipe, "set output 'plot%d.png'\n", frame); fprintf (gnuplotPipe, "set title 't = %d' font ',25'\n", (int)(t + 0.9)); fprintf (gnuplotPipe, "plot 'prof.dat' u 1:2 w l lw 5 lt rgb'#11BB11' t 'Radial'," "'prof.dat' u 1:3 w l lw 5 lt rgb '#BB11BB' t 'Azimutal'\n"); fflush (gnuplotPipe); frame++; } event stop (t = end){ system("ffmpeg -r 25 -f image2 -i plot%d.png -c:v libx264 -vf format=yuv420p -y mov.mp4 && " "rm -f plot*"); return 1; } /** ## Reference Kizner, Z., Makarov, V., Kamp, L., & Van Heijst, G. (2013). *Instabilities of the flow around a cylinder and emission of vortex dipoles*. Journal of Fluid Mechanics, 730, 419-441. doi:10.1017/jfm.2013.345 */