sandbox/Antoonvh/two_rings.c

    drawing

    The experiment of Lim and Nickels (1992).

    The collision of two vortex rings.

    400000 dye particles (Video via antoonvanhooft.nl)

    Set-up

    #include "grid/octree.h"
    #include "navier-stokes/centered.h"
    #include "navier-stokes/perfs.h"
    #include "fractions.h"
    #include "view.h"
    #include "lambda2.h"
    #include "tracer-particles.h"
    #include "scatter2.h"
    
    #define RADIUS (sqrt(sq(y) + sq(z)))
    
    Particles red, blue;
    
    scalar f[];
    int maxlevel = 11;
    double ti = 4., ue = 0.008;
    double Re = 1750.;
    double tend = 120. + 0.1;
    int np = 2e5;

    Boundary conditions

    Initially, there are two opposing jets at the left an right wall. The top and bottom boundary can leak fluid.

    u.n[left]   = dirichlet( f[]  *(1.) * (t <= ti));
    u.n[right]  = dirichlet(-f[-1]*(1.) * (t <= ti));
    u.n[top]    = neumann (0.);
    p[top]      = dirichlet (0.);
    pf[top]     = dirichlet (0.);
    u.n[bottom] = neumann (0.);
    p[bottom]   = dirichlet (0.);
    pf[bottom]  = dirichlet (0.);
    
    int main() {
      L0 = 32.;
      X0 = Y0 = Z0 = -L0/2;
      const face vector muc[] = {1./Re, 1./Re, 1./Re};
      mu = muc;
      run();
    }

    The simulations starts when the jets are triggered. We have a guess at an initial grid.

    event init (t = 0) {
      refine (RADIUS < 2.5 && fabs(x) > 9.*L0/20.   && level < (maxlevel - 1));
      refine (RADIUS < 1.5 && fabs(x) > 19.5*L0/40. && level < (maxlevel));
      f.refine = f.prolongation = fraction_refine;
      fraction (f, 1. - RADIUS);
      boundary ({f});
      new_tracer_particles (0); //This should not be needed
    }

    During the injection phase, the orifice shape is recomputed.

    event inject (i++; t <= ti) {
      fraction (f, 1. - RADIUS);
      boundary ({f});
    }

    Particles

    Tracer particles are added as a sort of dye. This is done right after the injection. Special care is taken to make this step MPI compatible.

    event add_tp (t = ti) {
      if (pid() == 0) {
        red  = new_tracer_particles (np);
        blue = new_tracer_particles (np);
        int j = 0;
        double Rseed = 1.8; //The rings are larger than the opening
        while (j < np) {
          double yp = noise();
          double zp = noise();
          if (sq(yp) + sq(zp) < sq(1.)) {
    	pl[red][j].x  = X0 +      fabs(noise()*ti*0.35);
    	pl[red][j].y  = yp*Rseed;
    	pl[red][j].z  = zp*Rseed;
    	pl[blue][j].x = X0 + L0 - fabs(noise()*ti*0.35);
    	pl[blue][j].y = yp*Rseed;
    	pl[blue][j].z = zp*Rseed;
    	j++;
          }
        }
      } else {
        red  = new_tracer_particles (0);
        blue = new_tracer_particles (0);
      }
      particle_boundary (red);
      particle_boundary (blue);
      set_particle_attributes (red);
      set_particle_attributes (blue);
    }

    Output

    Two movies are generated, plotting an adaptive \lambda_2-isosurface value, a slice of the the cells and the the particles.

    event snapshots (t += 0.1) {
      char str[99];
      sprintf (str,"Re = %g", Re);
      double val = -0.0001;
      scalar l2[];
      lambda2 (u, l2);
      foreach()
        l2[] = l2[] < val ? l2[] : nodata;
      stats l2s = statsf(l2);
      foreach()
        l2[] = l2[] < val ? l2[] : 0;
      boundary ({l2});
      int w = 1280;
      double tzoom = 100;
      double fov = t < tzoom ? 15 - t/80 : 15 - t/80 - (t - tzoom)/2.5;
      view (fov = 15 - t/80, width = w, height = 9*w/16, bg = {0.3, 0.4, 0.6},
    	theta = 0.6 + 0.15*cos(t/15), phi = 0.6, ty = 0.05);
      isosurface ("l2", min (-l2s.stddev, val));
      translate (y = -L0/4) {
        cells (n = {0,1,0});
      }
      draw_string (str, pos = 2, lw = 3);
      save ("movl2.mp4");
      
      if (blue) {
        clear();
        view (bg = {0.9, 0.9, 0.9});
        scatter (red,  s = 7, pc = {0.8, 0.2, 0.2});
        scatter (blue, s = 7, pc = {0.2, 0.2, 0.8});
        translate (y = -L0/4) {
          cells (n = {0,1,0});
        }
        draw_string (str, pos = 2, lw = 3);
        save ("movp.mp4");
      }
    }
    
    event dumps (t += 10) {
      char fn[99];
      sprintf (fn, "dump%g", t);
      dump (fn);
    
      sprintf (fn, "red_blue%g", t);
      pdump (fn);
    }
    
    event adapt (i++)
      adapt_wavelet ((scalar*){u}, (double[]){ue, ue, ue}, maxlevel);
    
    event stop (t = tend);

    Results

    Visualization of the \lambda_2 iso-surface reveals detail:

    See also via vimeo.

    drawing

    There is also this visualization from another run with tracers:

    drawing