sandbox/M1EMN/Exemples/granular_column_cohesif.c

    Collapse of cohesive granular columns

    Problem

    We look at a 2D column of cohesive grains collapsing on a flat surface.

    Snapshot of the heap. t=0 Snapshot of the heap. t=1 Snapshot of the heap. t=2

    Equations

    We propose an implementation of the Jop Pouliquen Forterre µ(I) rheology plus cohesion. It is just an extra term: \displaystyle \tau = (Co + \mu(I) P ) \frac{D}{|D|} the equivalent viscosity \displaystyle \eta =\frac{Co + \mu(I) P}{\sqrt{2} D_2} is included in Navier Stokes in http://basilisk.fr/sandbox/M1EMN/Exemples/granular.h

    Code

    Includes and definitions, note granular.h

    #include "granular.h"
#include "heights.h"
// Domain extent
#define LDOMAIN 5.
// heap definition
double  H0,R0,xfront;
//film size
double  Lz;
// Maximum refinement
#define LEVEL 8
// 8 is OK
#define LEVELmin 2
char s[80];
FILE * fpf;

    Boundary conditions for granular flow, pressure must be zero at the surface

    p[top]      = dirichlet(-RHOF*LDOMAIN);
pf[top]     = dirichlet(-RHOF*LDOMAIN);
u.n[top]    = neumann(0);
u.t[bottom] = dirichlet(0);

    the three cases in the main

    int main() {
  L0 = LDOMAIN;
  // number of grid points
  N = 1 << LEVEL;
  // maximum timestep
  DT = 25.e-3;
  TOLERANCE = 1e-3;

    Initial conditions aspect ratio a=.5 or .25 by symmetry Grain size D, value of the \mu_s \Delta \mu and I_0 of \mu(I) law, and gravity \overrightarrow{g}=-g \overrightarrow{e}_y

      H0=1;
  R0=2.0001;
  D=1./32;
  //mus=.32;
  mus=.4;
  dmu=.28;
  I0=0.4;

  const face vector gravity[] = {0.,-1.};
  a = gravity;
 
  FILE * f = fopen ("tauYxmax.txt", "w");

  fpf = fopen ("out000", "w");
  Lz = LDOMAIN;
  Co=0.00;
  run();
  fprintf (f, "%g %g  \n", Co, xfront );
  fclose (fpf); 
  system("cp velo.mpg velo0.00.mpg");

#define loc
    
#ifdef loc
  fpf = fopen ("out005", "w");
  Co=0.05;
  run();
  fprintf (f, "%g %g  \n", Co, xfront );
  fclose (fpf); 
  system("cp velo.mpg velo0.05.mpg");
    

  fpf = fopen ("out010", "w");
  Co=0.10;
  run();
  fprintf (f, "%g %g  \n", Co, xfront );
  fclose (fpf); 
  system("cp velo.mpg velo0.10.mpg"); 

  Co=0.15;
  fpf = fopen ("out015", "w");
  run();
  fprintf (f, "%g %g  \n", Co, xfront );
  system("cp velo.mpg velo0.15.mpg"); 
  Co=0.10; 
#endif

}

    initial heap, a rectangle

    event init (t = 0) {
  scalar phi[];
  foreach_vertex()
    phi[] = min(H0 - y, R0 - x);
  fractions (phi, f);
  /*
initialisation of hydrostatic pressure for granular phase  
*/
  foreach()
    p[] = f[]*(H0-y);
}

    outputs

    convergence outputs

    void mg_print (mgstats mg)
{
  if (mg.i > 0 && mg.resa > 0.)
    fprintf (stderr, "#   %d %g %g %g\n", mg.i, mg.resb, mg.resa,
             mg.resb > 0 ? exp (log (mg.resb/mg.resa)/mg.i) : 0);
}

    convergence outputs

    event logfile (i++) {
  stats s = statsf (f);
  fprintf (stderr, "%g %d %g %g %g %g\n", t, i, dt, s.sum, s.min, s.max - 1.);
  mg_print (mgp);
  mg_print (mgpf);
  mg_print (mgu);
  fflush (stderr);
}

    surface output

    event interface (t = {0, 1. , 2., 3., 4., 5.}) {
  output_facets (f, fpf); 
  char s[80];
  sprintf (s, "field-%g-%g.txt", t, Co);
  FILE * fp = fopen (s, "w");
  output_field ({f,p,u,uf,pf}, fp, linear = true);
  fclose (fp);
}

    if QUADTREE

    if no #include “grid/multigrid.h” then adapt

    event adapt(i++){

 scalar K1[],K2[];
 foreach()
   K1[]=(f[0,1]+f[0,-1]+f[1,0]+f[-1,0])/4*noise();
 boundary({K1});
 
 for(int k=1;k<8;k++)
 {
 foreach()
   K2[]=(K1[0,1]+K1[0,-1]+K1[1,0]+K1[-1,0])/4;
 boundary({K2});
 foreach()
   K1[]=K2[]*noise();
}
 adapt_wavelet({K1,f},(double[]){0.001,0.01}, maxlevel = LEVEL, minlevel = LEVELmin);

}
    “/*" within comment [-Wcomment]

    film output

    #if 0
event movie (t += 0.05) {
  static FILE * fp1 = popen ("ppm2mpeg > level.mpg", "w");
  scalar l[];
  foreach()
    l[] = level*(1-f[]);
    boundary ({l});
  output_ppm (l, fp1, min = 0, max = LEVEL,
	      n = 2048, box = {{0,-1},{Lz,Lz}});

  foreach()
    l[] = f[]*(1+sqrt(sq(u.x[]) + sq(u.y[])));
  boundary ({l});
  static FILE * fp2 = popen ("ppm2mpeg > velo.mpg", "w");
  output_ppm (l, fp2, min = 0, max = 1.5, linear = true, 
              n = 1024, box = {{0,-1},{4,1.5}});

  static FILE * fp3 = popen ("ppm2mpeg > f.mpg", "w");
  foreach()
    l[] = f[]*p[];
    boundary ({l});
  output_ppm (l, fp3, min = 0, max =1, linear = true,
              n = 2048, box = {{0,0},{Lz,Lz}});
}
#else
event movie (t += 0.05) {
    scalar l[];
    foreach()
    l[] = level*(1-f[]);
    boundary ({l});
    output_ppm (l, file = "level.mp4", min = 0, max = LEVEL,
                n = 2048, box = {{0,-1},{Lz,Lz}});
    
    foreach()
    l[] = f[]*(1+sqrt(sq(u.x[]) + sq(u.y[])));
    boundary ({l});
    output_ppm (l, file = "velo.mp4", min = 0, max = 1.5, linear = true,
                n = 1024, box = {{0,-1},{4,1.5}});
    
    foreach()
    l[] = f[]*p[];
    boundary ({l});
    output_ppm (l, file = "f.mp4", min = 0, max =1, linear = true,
                n = 2048, box = {{0,0},{Lz,Lz}});
}

#endif

vector h[];
event timeseries (t += 0.1 ) {
    heights (f, h);
    double maxy = - HUGE,maxx = - HUGE;;
    foreach()
    if ((h.y[] != nodata) && (h.x[] != nodata)) {
        double yi = y + height(h.y[])*Delta;
        double xi = x + height(h.x[])*Delta;
        if (yi > maxy)
            maxy = yi;
        if (xi > maxx)
            maxx = xi;
    }
    char s[80];
    sprintf (s, "hmax-%g",Co);
    static FILE * fp0 = fopen (s, "w");
    fprintf (fp0, "%g %g %g\n", t, maxx, maxy);
    fflush (fp0);
    xfront=maxx;
}

event pictures (t = {0, 1. , 2., 3., 4.} ) {
    if(Co<0.0000001){
    scalar l[];
    foreach()
      l[] = f[]*p[];
    boundary ({l});
    if(t==0.)
    output_ppm (f, file = "f0.png", min=-1, max = 2,  spread = 2, n = 256, linear = true,
                box = {{0,0},{4,2}});
    if(t==1.)
    output_ppm (f, file = "f1.png", min=-1, max = 2,  spread = 2, n = 256, linear = true,
                box = {{0,0},{4,2}});
    if(t==2.)
    output_ppm (f, file = "f2.png", min=-1, max = 2,  spread = 2, n = 256, linear = true,
                box = {{0,0},{4,2}});
    if(t==3.)
    output_ppm (f, file = "f3.png", min=0, max = 2,  spread = 2, n = 256, linear = true,
                box = {{0,0},{4,2}});}

}

    Run

    to run without make

    qcc -g -O2 -DTRASH=1 -Wall  -o granular_column granular_column.c -lm
./granular_column > out

    to run with make

    make granular_column_cohesif.tst
make granular_column_cohesif/plots
make granular_column_cohesif.c.html;

    Results

    First we check that we have exactly the same results than in the pure dry column test (small differences are maybe due to mesh adaptation) and with superposition of contact dynamics method DCM:

     d=0.005
 h0=0.149
 set xlabel "x"
 set ylabel "h(x,t)"
 
 p[0:4][0:1.5]'../../granular_column/ShapeTime.A-01.dat' u (($1*d)/h0):(($2*d)/h0) t'DCM'w l,\
 'out000' w l t 'this code','../REFCASES/out05_granular_column.dat' t'for control' w l
    (script)

    (script)

    Comparison for several values of Co of the change of chape at time 1,2,3 and 4

    set xlabel "x"
set ylabel "h(x,t)"
p[0:4][0:1.5]'out000' w l t 'Co=0.0',\
             'out005' w l t 'Co=0.05',\
             'out010' w l t 'Co=0.10','out015' w l t 'Co=0.15'
    collapse Continuum Model various Co (script)

    collapse Continuum Model various Co (script)

    Iso plot of velocity and pressure at t=2

    set pm3d map
set palette rgbformulae 22,13,-31;
unset colorbox
set multiplot layout 2,1
set xlabel "x  iso u"
set ylabel "y"
 splot [:4][:1] 'field-3-0.15.txt' u 1:2:($3>.999? $5*$3 :NaN)   not
set xlabel "x  iso p"
 splot [:4][:1] 'field-3-0.15.txt' u 1:2:($3>.999? $4*$3 :NaN)   not
unset multiplot
reset
    velocity and pressure (script)

    velocity and pressure (script)

    a zoom on the solid corner… at time 1,2,3 and 4

    set view 0,0
set contour;unset surface
splot [0:3][0:1.5] 'field-1-0.1.txt' u 1:2:(sqrt($5*$5+$6*$6)*$3)  w l not, \
                  'field-2-0.1.txt' u 1:2:(sqrt($5*$5+$6*$6)*$3)   w l not,\
                  'field-3-0.1.txt' u 1:2:(sqrt($5*$5+$6*$6)*$3)   w l not,\
                  'field-4-0.1.txt' u 1:2:(sqrt($5*$5+$6*$6)*$3)   w l not
    (script)

    (script)

    Looking at the rigid corner

    set view 0,0
set contour;unset surface
splot [0:3][0:1.5] 'field-2-0.15.txt' u 1:2:(sqrt($5*$5+$6*$6)*$3)  w l not, \
 'field-2-0.15.txt' u 1:2:(sqrt($5*$5+$6*$6)*$3)   w l not,\
 'field-2-0.15.txt' u 1:2:(sqrt($5*$5+$6*$6)*$3)   w l not,\
 'field-2-0.15.txt' u 1:2:(sqrt($5*$5+$6*$6)*$3)   w l not
    (script)

    (script)

    Comparaison of runout for different cohesions as function of time

     reset
 set key left
 set xlabel "t"
 set ylabel "h_m(t),x_m(t)"
 p'hmax-0' u 1:2 w l t 'Co=0.00','' u 1:3 w l t 'Co=0.00',\
 'hmax-0.05' u 1:2 w l t 'Co=0.05','' u 1:3 w l t 'Co=0.05',\
 'hmax-0.1' u 1:2 w l t 'Co=0.10','' u 1:3 w l t 'Co=0.10',\
 'hmax-0.15' u 1:2 w l t 'Co=0.15','' u 1:3 w l t 'Co=0.15'
    (script)

    (script)

    Comparaison of runout for 1D Savage Hutter, 2D Multilayer (RNSP), and 2D NS as function of cohesion

     reset
 set xlabel "tauY"
 set ylabel "(x max- x init)"
 p'tauYxmax.txt' u ($1):($2-2)  w lp  t'NS',\
  '../cohesive_muI_collapse_ML/tauYxmax.txt'u ($1):($2-2)  t '2D RNSP' w l,\
  '../cohesive_savagehutter/tauYxmax.txt' u ($1):($2-2) t'1D SH' w l
    (script)

    (script)

    Some films, Velocity

    velocity (click on image for animation)

    level

    velocity (click on image for animation)

    collapse

    velocity (click on image for animation)

    Bibliography

    CISM 2019