# sandbox/Antoonvh/ash.c 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 end (t = 10);