Collapse of granular heap


    This is the problem of the collapse of a initial rectangular heap of grains. It is a kind of Dam Break problem.

    animation of the collapse

    animation of the collapse


    We compare the contact dynamic solution of granular collapse with “Shallow Water” or “Savage Hutter equations” or “depth averaged equations” or “Saint Venant” simplified equations with a basal friction \mu(I)P. Here P=\rho gh and I= 5./2.d (u/h)(g h)^{-1/2} . The system is \displaystyle \left\{\begin{array}{l} \frac{\partial }{\partial t} h \; +\; \frac{\partial }{\partial x} uh=s\\ \frac{\partial }{\partial t} hu + \frac{\partial }{\partial x} \dfrac{(hu)^2}{h} + \frac{\partial }{\partial x}g\dfrac{h^2}{2} = - gh \frac{\partial }{\partial x} Z-\mu g h\frac{u}{|u|} \end{array}\right.

    Some pluviation s \ne 0 will be added as suggested in Larieu et all. (but is not yet implemented: here s = 0)

    We solve the problem by splitting, first (note that we use Q=hu) \displaystyle \frac{h^*-h^n}{\Delta t} + \frac{\partial Q^n}{\partial x} =0, \text{ and } \frac{Q^*-Q^n}{\Delta t}+ \frac{\partial }{\partial x}\frac{Q^n}{h^n} + \frac{\partial }{\partial x}g\dfrac{h^{n2}}{2} = - g \frac{\partial }{\partial x} Z

    is solved with

    Second, (h is simplified) \displaystyle \frac{u^{n+1}-u^*}{\Delta t} =- \mu g \frac{u}{|u|} is solved.


    #include "grid/multigrid1D.h"
    #include "saint-venant.h"

    The domain is 5 long, the height is the unit of length. The problem is without dimensions. The problem is solved in one dimension.

    Wall symmetry at the left Neumann conditions at the exit

    u.n[left] = dirichlet(0);
    h[left] = neumann(0);
    u.n[right] = neumann(0);
    h[right] = neumann(0);

    Main and parameters

    //int source;
    double Ltas,Htas,tmax;
    int main()
      X0 = 0.;
      L0 = 7;
      G = 1.; 
      N = 1024;
    //  source = 0;
      Htas = 1;
      Ltas = 2.; 

    The initial conditions are a given heap of length L_{tas} (the double by symetry) of height H_{tas}. to left.

    event init (i = 0){
      double h0;

    initial heap

        h0=Htas;//0.1*Htas*source + (1.-source)*Htas;
        zb[] = 0;	
        h[] = (x < Ltas) ? h0 : 0;
        u.x[]= 0;}


    We use a simple implicit scheme to implement coulomb bottom friction i.e. \frac{d\mathbf{u}}{dt}=-\mu g \frac{\mathbf{u}}{|\mathbf{u}|} with \mu fonction of a mean I (written with Q/h^{3/2}).

    Of course we have splitted, first the Rieman Problem, second the friction problem.

    • first we tried (as in viscous or turbulent cases) a semi implicit discretisation like: u(t+\Delta t) = u(t ) -\mu g (\Delta t) \frac{u(t+\Delta t) }{u(t ) } to obtain \displaystyle u(t+\Delta t) = \frac{u(t)}{1 + \mu g (\Delta t) \frac{1 }{u(t)}}

    • finaly (see granular sand glass with friction) we noticed that there is an exact solution of this splited problem, defining the norm of velocity U=|\mathbf{u}| and \overrightarrow{T}=\frac{\mathbf{u}}{|\mathbf{u}|} the equation \frac{d (U \overrightarrow{T})}{dt}=-\mu g \overrightarrow{T} is solved explicitely for t'>t, obviously, as \frac{d \overrightarrow{T}}{ds}=\frac{\overrightarrow{N}}{R} we have just to solve: \displaystyle \frac{d U}{dt}=-\mu g it is linear and the solution is U(t')=U(t)- \mu g (t'-t) down to stop, so: \displaystyle U(t+\Delta t) = max(U(t) - \mu \Delta t g,0) This gives a real stop.

    If \mu is function of a mean I across the layer, this mean I is obtained for a Bagnold profile as (2/5)(d/h)Q/\sqrt{g h } and \mu(I)=\mu_0 +\Delta \mu I/(I+I_0)

    event coulomb_friction (i++) {
      double In,nbrg=32.,eps=0.00000,mu,U,tauY=.001;
      foreach() {
        mu = (0.4 + 0.28*In/(0.4+In)) ;
      // mu = (0.4 + 0.26*exp(-0.136/In));  
       //   mu=.4;//.4; // to tes Kerswell
    // obsolete
    //    ff = U > 0 ? (1. +  dt *mu*G/U) : HUGE ;
    //    foreach_dimension()
    //      u.x[] /= ff;
    // new and better
              u.x[] = max(U -dt *( mu * G + tauY/h[]),0)*u.x[]/U;}
      boundary ({u.x});


    event outputfront (t += .1 ) {
        double  xf=0,xe=0;
        static FILE * ff = fopen("x.txt", "w");

    tracking the front and the end of the heap

            xf = h[] > 1e-20 ?  max(xf,x) :  xf ;
            xe = h[] > 1e-20 ?  min(xe,x) :  xe ;
         fprintf (ff, "%g %g %g \n", t, xf , xe);
    event output (t += .01; t < tmax) {
        static FILE * fp = popen ("gnuplot -persist 2> /dev/null", "w");
    #ifdef gnuX
        if(t==0) fprintf (fp,"set term gif animate;set output 'animate.gif';set size ratio .5\n");
        fprintf (fp,"set title ' collapse --- t= %.2lf '\n"
                 "p[0:5][-.5:2]  '-' u 1:($2) t'free surface' w lp lt 3,0 not w l linec -1\n",t);
          fprintf (fp, "%g %g \n", x, h[]);
      fprintf (fp,"e\n\n");
      event outputlog (t += 1; t < tmax) {
         fprintf (stderr, "%g %g %g \n", x, h[], t);
       fprintf (stderr, "\n");


    to run and plot with gnuplot during the run (note the -DgnuX)

    qcc -g -O2 -DTRASH=1 -Wall -DgnuX=1 -o savagestaron savagestaron.c -lm

    with make

     make savagestaron.tst;make savagestaron/plots; make savagestaron.c.html
     source ../ savagestaron

    plot with gnuplot of the results, dynamics is not very good, but deposit is not so bad.

    set xlabel 'x'
    set ylabel 'h'
    set xlabel "x"
    set ylabel "h(x,t)"
    p[0:6][0:1.5]'../../granular_column/ShapeTime.A-01.dat' u (($1*d)/h0):(($2*d)/h0) t'DCM'w l,'../savagestaron/log' w l t 'h'
    Fluid depth profile, Contact Dynamics vs Savage Hutter (script)

    Fluid depth profile, Contact Dynamics vs Savage Hutter (script)

    Comparison with Balmforth & Kerswell 05 for the runout proposed by Kerswell 05

    u = \frac{dx}{dt}=2 - \mu t so that x = 2 (t - \mu t^2/4) and t_{max} =\frac{2}{\mu} x_{max} =\frac{2}{\mu}

     set xlabel "t"
     set ylabel "xmax"
     mu =.4
     set key bottom
     p [][0:]'x.txt' u ($1):(($2-2)) t' calc',(x*(2-x*mu/2))*(x<2/mu?1:NaN) t'caract pred.'


    Comparison with Balmforth & Kerswell 05 who rescaled to have a canonical problem:

     set term png ;   set output 'bk.png';
     unset tics
     set key center top
     plot [0:1300]  '../Img/Balmforth_Kerswell05.png' binary filetype=png with rgbimage not,\
       'log' u (Xc+($1-2)*L1*mu):(Yc+$2*(H1)) w l not
    cmp B&K 05 (script)

    cmp B&K 05 (script)


    • Larrieu Staron Hinch “Raining into shallow water as a description of the collapse of a column of grains” J. Fluid Mech. (2006), vol. 554, pp. 259–270.

    • Balmforth & Kerswell 05 “Granular collapse in two dimensions” J. Fluid Mech. (2005), vol. 538, pp. 399–428

    • R. R. Kerswell Dam break with Coulomb friction: “A model for granular slumping?” PHYSICS OF FLUIDS 17, 057101 2005

    v1: Montpellier 11 juillet 15