# Bagnold free surface flow with segregation: Brazil Nuts effect

An example of 2D granular flow of a binary mixture of small and large particules along an inclined plate (angle -\alpha) with a free surface is presented here. The configuration is periodic, what is injected to the left comes from the right. On the bottom solid wall there is a no slip boundary condition, on the top fixed free surface a slip condition is imposed (and zero pressure). We use the centered Navier-Stokes solver with regularization for viscosity. This flow is a simple model for of mixture of pebbles and gravels avalanche along a slope.

## SeGrayGation

Following Gray and Thornton we incorporate advection due to mean flow (\frac{\partial c_i}{\partial t}+ u \frac{\partial c_i}{\partial x}+ v\frac{\partial c_i}{\partial y}), percolation-driven segregation \frac{\partial }{\partial y}(c_i v_{pi}) and diffusion due to random particle collisions \frac{\partial }{\partial x}(D \frac{\partial }{\partial x}c_i)+ \frac{\partial }{\partial y}(D \frac{\partial }{\partial y}c_i). This gives a continuum transport equation for the volume concentration for a binary mixture of species i=l,s (i = l and i = s represent large particles and small particles, respectively). The percolation velocity may be approximated by

v_{pl}=S_r(\frac{\partial u}{\partial y})(1-c_l) and v_{ps}=-S_r(\frac{\partial u}{\partial y})(1-c_s) where S_r is parameter depending on the granular media. D may be approximated by D_2 \sqrt{2}d_g^2, here it is just a constant. The final equation:

\displaystyle \frac{\partial c_i}{\partial t}+ u \frac{\partial c_i}{\partial x}+ v\frac{\partial c_i}{\partial y}+ S_r \frac{\partial }{\partial y}\left( (\frac{\partial u}{\partial y})c_i(1-c_i) \right)= \frac{\partial }{\partial x}(D \frac{\partial }{\partial x}c_i)+ \frac{\partial }{\partial y}(D \frac{\partial }{\partial y}c_i)

Boundary condition, equilibrium at the top \displaystyle S_r (\frac{\partial u}{\partial y})|_h c_i(h)(1-c_i(h)) = - D \frac{\partial }{\partial y}c_i|_h and equilibrium at the wall \displaystyle S_r (\frac{\partial u}{\partial y})|_0 c_i(0)(1-c_i(0)) = - D \frac{\partial }{\partial y}c_i\_0

## Equations for a granular fluid: \mu(I) rheology

We have to solve the velocity field which carries the mixture. By GDR MiDi, in a scalar analysis \tau = \mu(I ) p, by Jop et al hypothesis, this is generalized : tangential stress is linked to the shear rate tensor by \displaystyle \tau_{ij} = \sqrt{2}\mu(I) p \frac{D_{ij}}{D_2} where the friction law compatible with the experiments is: \displaystyle \mu(I)= \mu_0 + \frac{\Delta \mu}{\frac{I_0}{I}+1} the coefficients depend on the nature of the granular media \mu_0=0.38 \Delta \mu = 0.26 and I_0=0.3

## Equivalent rheology

From this, one can exhibit an equivalent viscosity if we write: ${ij} = 2 {eq} D_{ij}$ we the define \mu_{eq}= \frac{\mu(I) p}{\sqrt{2} D_2 } This equivalent viscosity will be coded next (with a regularisation to avoid infinite viscosity at rest).

## Exact solution in the proposed case

This problem admits an analytical solution called “Bagnold” profile in a steady case. \displaystyle u = \frac{2}{3} \frac{I_\alpha}{d_g} \sqrt{g h^3 cos(\alpha)} (1 - (1 - \frac{y}{h})^{3/2})

just before, let us include the NS code, define the precision 2^N, define the angle 24.64 degrees

#include "navier-stokes/centered.h"
#include "tracer.h"
#include "diffusion.h"
scalar c[],q0[],q[];
scalar * tracers = {c};
mgstats mgf;

#define LEVEL 6
double mumax;
double dg;
double Sr=.5;
# define alpha 0.43
scalar mu_eq[],foo[];

the exact Bagnold shear velocity is:

double dUb( double y){
double U0=(sqrt(cos(alpha))*(-0.114 + 0.3*tan(alpha))/(dg*(0.64 - tan(alpha))));
return ((pow(1. - y, .5))*U0) ;
}

This gives the Bagnold’s solution (note the numerical value for alpha=0.43, dg=0.04 U0 = 2.06631)

double Ub( double y){
double U0=(sqrt(cos(alpha))*(-0.114 + 0.3*tan(alpha))/(dg*(0.64 - tan(alpha))));
return ((1. - pow(1. - y, 1.5))*2./3.*U0) ;
}

double zz( double t){
return 10 ;
}

The domain is one unit long. 0<x<1 0<y<1

int main() {
L0 = 1.;
origin (0., 0);

Values of grain size

    //  alpha = 0.43;
dg = 0.04;

the regularisation value of viscosity

    mumax=1000;

Boundary conditions are periodic

    periodic (right);

slip at the top

    u.t[top] = neumann(0);
//  uf.t[top] = neumann(0);
u.n[top] = dirichlet(0);

no slip at the bottom

    u.n[bottom] = dirichlet(0);
uf.n[bottom] = dirichlet(0);
u.t[bottom] = dirichlet(0);

note the pressure

    p[top] = dirichlet(0);
pf[top] = dirichlet(0);

flux note the signs (check!)

    q0[bottom] = neumann(0);
q0[top] = neumann(0);

c[bottom] = neumann(q0[]);
c[top] = neumann(-q0[]);

the \Delta t_{max} should be enough small

    DT = 0.01;
run();
}
face vector muv[];

event init (t = 0) {

prepare viscosity

    mu = muv;

pressure gradient, gravity acceleartion mdpdx \displaystyle -\frac{\partial p}{\partial x} = sin(alpha) \displaystyle -\frac{\partial p}{\partial y} = -cos(alpha)

    const face vector mdpdx[] = {sin(alpha),-cos(alpha)};

a = mdpdx;

Initialy at rest

    foreach() {
u.x[] = Ub(y);
u.y[] = 0;
p[]=cos(alpha)*(1-y);
c[] = ( (y<.75) && ( y >.25) ? 1 : 0);
c[] = ( (y<.75 && x<2./3) ? 1 : 0);
}
}

We check the number of iterations of the Poisson and viscous problems.

//event logfile (i++)
// fprintf (stderr, "%d %g %d %d\n", i, t, mgp.i, mgu.i);

old value of the velocity and concentration is saved

scalar un[],cn[];
event init_un (i = 0) {
foreach(){
un[] = u.x[];
cn[] = c[];}
}

so that when it does not more change we are converged, monitor mass conservation as well

event conv (t += 1; i <= 100000) {
double du = change (u.x, un);
double dc = change (c, cn);
double sc = statsf(c).sum;

fprintf(stdout,"t= %g u/Ub=%g mass = %g \n",t,interpolate (u.x, L0/2, .999)/Ub(1),sc);
if (i > 0 && du < 1e-5 && dc < 1e-5)
return 1; /* stop */
}

## Implementation of the Bagnold viscosity

event bagnold(i++) {

Compute \displaystyle D_2^2= D_{ij}D_{ij}= D_{11}D_{11} + D_{12}D_{21} + D_{21}D_{12} + D_{22}D_{22} with I = d_g \sqrt{2} D_2 /\sqrt{p/\rho}, and as \displaystyle \mu(I)= \mu_0 + (\Delta \mu) I/(I_0 + I) the equivalent viscosity is \displaystyle \mu_{eq}= \frac{\mu(I) p}{\sqrt{2} D_2 } the final one is the min of of \mu_{eq} and a large \mu_{max}, then the fluid flows always, it is not a solid, but a very viscous fluid.

    foreach() {
mu_eq[] =    mumax;
if (p[] > 0.) {
double D2 = 0.,In,muI;
foreach_dimension() {
double dxx = u.x[1,0] - u.x[-1,0];
double dxy = (u.x[0,1] - u.x[0,-1] + u.y[1,0] - u.y[-1,0])/2.;
D2 += sq(dxx) + sq(dxy);
}
if (D2 > 0.) {
D2 = sqrt(D2)/(2.*Delta) ;
In = sqrt(2.)*dg*D2/sqrt(fabs(p[]));
muI = .38 + (.26)*In/(.3 + In);
foo[] = D2;
mu_eq[] =  min(muI*fabs(p[])/(sqrt(2.)*D2) , mumax );}
else{
mu_eq[] =  mumax;
}
}
}
boundary ({mu_eq});
foreach_face() {
muv.x[] = (mu_eq[] + mu_eq[-1,0])/2.;
}
boundary ((scalar *){muv});
}

## Mixing

The Gray equation is solved here

event mixing(i++){
const face vector D[] = {.01 , .01 };
scalar r[],logis[],shear[];

Compute here the “logistic map source term” from percolation in computing a field \displaystyle q_0 = - \frac{S_r}{D} (\frac{\partial u}{\partial y})| c_i(1-c_i) hence at the top and bottom \displaystyle \frac{\partial }{\partial y}c_i = -q_0(x,y) so that at the top \displaystyle \frac{\partial }{\partial y}c_i|_h = -q_0(x,h) so that at the top, as y=h+n \displaystyle \frac{\partial }{\partial n}c_i|_h = -q_0(x,h) which was coded before as c[top] = neumann(-q0[]);. So that at the bottom, as y=-n \displaystyle \frac{\partial }{\partial n}c_i|_0 = q_0(x,0) mind the minus sign do to the normal in the BC at the wall !

    foreach() {
shear[] = (u.x[0,1]-u.x[0,-1])/(2*Delta);
logis[] = c[]*(1-c[])*(1-.0*c[])*shear[];
q0[] =  -c[]*(1-c[])*Sr/.01*shear[];
q[] = (c[0,1]-c[0,-1])/(2*Delta);}
boundary ({logis,shear,q0,q});

Compute the source term r for the diffusion equation \displaystyle r=-S_r \frac{\partial }{\partial y}\left( (\frac{\partial u}{\partial y})c_i(1-c_i) \right)

    foreach() {
r[] = -Sr*(logis[0,1]-logis[0,-1])/(2*Delta);}
boundary ({r});

Diffusion equation

\displaystyle \frac{\partial c_i}{\partial t}+ u \frac{\partial c_i}{\partial x}+ v\frac{\partial c_i}{\partial y}= \frac{\partial }{\partial x}(D \frac{\partial }{\partial x}c_i)+ \frac{\partial }{\partial y}(D \frac{\partial }{\partial y}c_i)+r

	mgf= diffusion (c, dt, D, r);
}

Save profiles + film

event movies (i += 10 ) {
static FILE * fp = popen ("ppm2mpeg > vort.mpg", "w");
scalar vorticity[];
foreach()
vorticity[] = (u.x[0,1] - u.x[0,-1] - u.y[1,0] + u.y[-1,0])/(2.*Delta);
boundary ({vorticity});
output_ppm (vorticity, fp, box = {{0,0},{1,1}},
min = 0, max = 2, linear = true);
static FILE * fp1 = popen ("ppm2mpeg > c.mpg", "w");
output_ppm (c, fp1, box = {{0,0},{1,1}},
linear = true, min = 0, max = 1);
}

event profiles (t += 1)
{
FILE * fp = fopen("xprof", "w");
scalar shear[];
foreach()
shear[] = (u.x[0,1] - u.x[0,-1])/(2.*Delta);
boundary ({shear});
for (double y = 0.; y < 1.0; y += 1./pow(2.,LEVEL))
fprintf (fp, "%g %g %g %g %g %g %g %g %g %g %g\n", y, interpolate (u.x, L0/2, y), interpolate (shear, L0/2, y),
Ub(y), cos(alpha)*(1-y),interpolate (p, L0/2, y),
interpolate (mu_eq, L0/2, y), interpolate (foo, L0/2, y),dUb(y), interpolate (c, L0/2, y),interpolate (q, L0/2, y));
fclose (fp);
}

We adapt according to the error on the velocity field.

event adapt (i++) {
// adapt_wavelet ({u}, (double[]){3e-3,3e-3}, 8, 6);
}

## Results and plots

To run the program

 qcc -g -O3 -o bagnold_periodic_segregation bagnold_periodic_segregation.c -lm
./bagnold_periodic_segregation

lldb bagnold_periodic

Plots of the velocity, \tau and p and concentration and flux in the middle of the domain.

 set xlabel "y"
set xlabel "U, tau, p, c"
set key left
p'xprof' u 1:2 t'U computed' ,''u 1:($7*$3) t'tau comp.','' u 1:6 t'p',''u 1:10 t'c(L0/2,y) large part.' w lp,''u 1:(\$11/15) t'flux dc/dy' w lp Bagnold flow and concentration of large particules (script)

See an animation of the concentration, red corrsponds to the large particules, blue the small

Animation of concentration field.

Gentilly Avril 2015