sandbox/M1EMN/Exemples/granular_column.c

    Collapse of granular columns

    Problem

    Geophysics, moutains…

    Equations

    We propose an implementation of the Jop Pouliquen Forterre µ(I) rheology. This is the counterpart in Basilisk of the gerris case.

    It consists to solve for three aspect ratios the collapse of three columns. Each collapse is compared to the contact dynamics solution.

    Code

    Includes and definitions, note granular.h, and here multigrid

    #include "grid/multigrid.h"
#include "granular.h"
// Domain extent
#define LDOMAIN 5.
// heap definition
double  H0,R0,D;
//film size
double  Lz;
// Maximum refinement
#define LEVEL 7
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 = 1e-2;
  TOLERANCE = 1e-3;
  
// Initial conditions a=.5
  H0=1;
  R0=2.0001;

  
  const face vector g[] = {0.,-1.};
  a = g;
    
  fpf = fopen ("out05", "w");
  Lz = LDOMAIN;
  run();
  fclose (fpf);
  system("gnuplot comp05.gnu");

// Initial conditions a=1.42
  H0=1;
  R0=1./1.42; 
  D=1./65;
  fpf = fopen ("out142", "w");
  Lz = 3;
  run();
  fclose (fpf);
  system("gnuplot comp142.gnu");
  
// Initial conditions a=6.26
  H0=1;
  R0=1./6.26; 
  D=1./65;
  fpf = fopen ("out626", "w");
  Lz = 2.5;
  run();
  fclose (fpf); 
  system("gnuplot comp626.gnu");  
  
 
}

    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);
}

    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);
}

    convergence outputs

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

    film output

    event movie (t += 0.05) {
  //static FILE * fp1 = popen ("ppm2mpeg > level.mpg", "w");
  scalar l[];
  foreach()
    l[] = level;
  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});
//  static FILE * fp2 = popen ("ppm2mpeg > velo.mpg", "w");
  output_ppm (l, file = "velo.mp4", min = 0, max = 1.5, linear = true,
	      n = 1024, box = {{0,-1},{5,1.5}});

 // static FILE * fp3 = popen ("ppm2mpeg > f.mpg", "w");
  foreach()
    l[] = f[]*p[];
  output_ppm (l, file = "f.mp4", min = 0, linear = true,
	      n = 2048, box = {{0,-1},{Lz,Lz}});
}
 event pictures (t=3) {
  output_ppm (f, file = "f.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.tst
make granular_column/plots
make granular_column.c.html;

    Results

    Exemples of comparisons Discrete Contacts Dynamics vs Continuum Model for 3 aspect ratio

    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,\
'../granular_column/out05' w l t 'h'
    collapse Contacts Dynamics vs Continuum Model a =0.5 (script)

    collapse Contacts Dynamics vs Continuum Model a =0.5 (script)

    set xlabel "x"
set ylabel "h(x,t)"
d=0.005
h0=0.324
p[0:4][0:1.5]'../Img/ShapeTime.A-04.dat' u (($1*d)/h0):(($2*d)/h0) t'DCM'w l,\
'../granular_column/out142' w l t 'h'
    collapse Contacts Dynamics vs Continuum Model a =1.42 (script)

    collapse Contacts Dynamics vs Continuum Model a =1.42 (script)

    set xlabel "x"
set ylabel "h(x,t)"
d=0.005
h0=0.67
p[0:4][0:1.5]'../Img/ShapeTime.A-03.dat' u (($1*d)/h0):(($2*d)/h0) t'DCM'w l,\
'../granular_column/out626' w l t 'h'
set term png
set output 'a.png'
 replot
    collapse Contacts Dynamics vs Continuum Model a =6.26 (script)

    collapse Contacts Dynamics vs Continuum Model a =6.26 (script)

    velocity (click on image for animation)

    Related examples

    Bibliography

    Montpellier 07/14