sandbox/Antoonvh/ash.c

    Turbudity currents are driven by particles. Image via Octavio E. Sequeiros et al. (2010)

    Turbudity currents are driven by particles. Image via Octavio E. Sequeiros et al. (2010)

    Vulcanic ash settling in water

    We use 8000 spherical sand particles to model the descend of ganules in water.

    You can “see” the two-way coupling in action

    The movie displays a surprisingly pleasing bview bug as is renders only part of the grid slice…

    #include "grid/octree.h"
    #include "navier-stokes/centered.h"
    #define TWO_WAY 1
    #include "stokes-particles.h"
    #include "view.h"
    #include "scatter2.h"
    
    u.t[bottom] = dirichlet (0.);
    u.r[bottom] = dirichlet (0.);
    
    Particles sand;
    
    #define MUZ 1e-3
    
    int main() {
      const face vector muc[] = {MUZ, MUZ, MUZ}; //mu Water
      periodic (left);
      mu = muc;
      G.y = -9.81;
      L0 = 0.1;
      X0 = Y0 = Z0 = -L0/2;
      mu = muc;
      N = 64;
      run();
    }

    Initially, the flow is steady and seeded with particles near the top of the domain.

    event init (t = 0) {
      const scalar rhof[] = 1000;
      rho = rhof;
      int pni = 20;
      sand = new_inertial_particles (cube(pni));
      int n = 0;
      for (int i = 0; i < pni; i++) 
        for (int j = 0; j < pni; j++)  
          for (int k = 0; k < pni; k++) { 
    	pl[sand][n].x = (double)i/(100*pni);
    	pl[sand][n].y = (double)j/(100*pni) + L0/4;
    	pl[sand][n].z = (double)k/(100*pni);
    	pl[sand][n].u2.x = 2000;                      // approx. rho Silicon
    	pl[sand][n].u2.y = 0.05e-3 + 0.01e-3*noise(); // fine sand-grain size
    	n++;
          }
    }

    Bottom boundary

    Particles bounce at the bottom boundary. They should be able to rest on there.

    double damp = 2;
    event bottom_boundary (i++) {
      foreach_particle_in(sand) {
        if (y < Y0) {
          p().y = Y0 + Y0 - y;
          foreach_dimension()
    	p().u.y /= damp;
          p().u.y *= -1;  
        }
      }
    }
    
    event adapt (i++) 
      adapt_wavelet ((scalar*){u}, (double[]){0.01, 0.01, 0.01}, 7, 4);
    
    event mov (t += .02) {
      view (theta = t/10. - 0.5);
      scatter (sand, pc = {1, 1, 1});
      box();
      translate (z = -L0/2)
        cells();
      save ("locs.mp4");
    }
    
    event end (t = 10);