sandbox/M1EMN/Exemples/bagnold_periodic_segregation.c

    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)

    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

    Bibliography