ash.c

Vulcanic ash settling in water
- Bottom boundary

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

Vulcanic ash settling in water

We use 64000 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; 
  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 = 40;
  if (pid() == 0) {
    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++;
	}
  } else
    sand = new_inertial_particles (0);
  particle_boundary (sand);
}

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);