sandbox/M1EMN/Exemples/cohesive_muI_collapse_ML.c

    collapse of a rectangular column of dry granular (with cohesion) on a slope,

    Problem

    A heap of dry granular (with cohesion) fluid is released along a constant slope.

    RNSP

    We solve the boundary layer (RNSP) equations corresponding to a thin layer approximation of the Navier Stokes equations: \displaystyle \left\{ \begin{aligned} \dfrac{\partial u}{\partial x} + \dfrac{\partial v}{\partial y} &= 0 \\ \rho (\dfrac{\partial u}{\partial t} + \dfrac{\partial u^2}{\partial x} + \dfrac{\partial u v}{\partial y} )&= - \rho g Z'_b - \dfrac{\partial p}{\partial x} + \dfrac{\partial \tau_{xy}}{\partial y} \\ 0 &= -\rho g -\dfrac{\partial p}{\partial y}. \end{aligned} \right. following \mu(I) rheology, with I=\dfrac{d}{\sqrt{p/\rho}}\dfrac{\partial u}{\partial y} this gives \tau_{xy}= \tau_Y + (\mu( \dfrac{d}{\sqrt{p/\rho}}\dfrac{\partial u}{\partial y})p), there is cohesion with additional \tau_Y (Yield stress) to the classical (\mu(I)p).

    We solve this with the Multilayer Solver technique from Audusse & Sainte Marie. We tune the solver to consider a non newtonian viscosity function of the shear \displaystyle \nu_{eq}= \frac{(\tau_Y+ \mu(I)p )/\rho}{|\frac{\partial u}{\partial y}|} the multilayer.h is changed in saint-venantNN.h (all is in saint-venantB.h).

    If we integrate over the depth of flow h, the pressure is hydrostatic p=\rho g (h-y). Integrating over the depth of flow h the incompressiblility equation, we obtain \displaystyle \frac{\partial h}{\partial t} + \frac{\partial Q}{\partial x}=0 where Q=\int_0^h udy. For momentum it gives: \displaystyle \frac{\partial Q}{\partial t} + \frac{\partial }{\partial x}(\frac{Q^2}{h} +g \frac{h^2}{2})= - g h Z'_b - \mu g h - \tau_Y

    this is not simple to put \tau_Y, Note that we cannot neglect inertia in the momentum equation. Indeed to solve for Q in \displaystyle 0+g \frac{\partial h }{\partial x}= - g Z'_b - \mu g with \mu constant is impossible.

    But, if \mu is function of a mean I_m=(2/5)(d/h)Q/\sqrt{g h } obtained for a Bagnold profile, then we can solve for Q in \displaystyle 0+g \frac{h^2}{2}= - g h Z'_b - \mu(I_m) g h and find a kind of 1D kinetic wave with Q = \frac{2 }{5} (\frac{-Z_b' - \mu_s -\partial_x h}{\Delta \mu}) \sqrt{g h} \frac{h^2}{d}. This flux is substituted in mass conservation; \displaystyle \frac{\partial h}{\partial t} + \frac{\partial Q}{\partial x}=0

    as in Huppert problem or Balmforth problem. This is done for this Bagnold case here

    In this file we use the Multilayer Solver to see the influence of the inertia neglected in the 1D model.

    Code

    #define pourtoutes(item, array) \
for(int keep = 1, \
count = 0,\
size = sizeof (array) / sizeof *(array); \
keep && count != size; \
keep = !keep, count++) \
for(item = (array) + count; keep; keep = !keep)


#include "grid/cartesian1D.h"
double tauY;
#include "saint-venantB.h"

    Non Newtonian viscosity

    The definition of viscosity as a function of shear (z is transverse if we use x,y in the plane): I= \dfrac{d}{\sqrt{p/\rho}}\dfrac{\partial u}{\partial z}, there is cohesion with additional \tau_Y to the classical (\mu(I)p). the non newtonian viscosity function of the shear and pressure is \displaystyle \nu_{eq}= \frac{(\tau_Y+ \mu(I)p )/\rho}{|\frac{\partial u}{\partial z}|}

    \tau_Y without dimension is refered as a “Bingham”, this is a Yield stress

    double nu_eq(double shear,double press){
    double nu_eq,In;
    In=(press>0?  1./32*sqrt(sq(shear))/sqrt((press)) : 1e30);
    nu_eq = min( ( tauY + (.4 + .28*In/(In+.4)) *  press)/sqrt(sq(shear)+dry),1e30 ) ;
    return nu_eq;}

    main

    double alpha,xfront;
char s[80];

int main() {
  X0 = 0;
  L0 = 8;
  G  = 1.;
  alpha = 0.;
  N  = 512*2;
  nl = 32*2;
  nu = 1;
  DT =0.01;
  FILE * f = fopen ("tauYxmax.txt", "w");

    a loop for the three values of Bingham parameter

       double values[] = { 0, 0.01, 0.02, 0.05,  0.1 , 0.2 , 0.3, 0.5 , 1 };
   pourtoutes(double *v, values) {
   tauY=*v;
   printf("# tauY %lf\n",tauY);
  
  sprintf (s, "xML-%.3f.txt", tauY);
  FILE * fp = fopen (s, "w"); 
  fclose(fp);
  sprintf (s, "shapeML-%.3f.txt", tauY);
  FILE * fs = fopen (s, "w");  
  fclose(fs);
        
  run();
   fprintf (f, "%g %g  \n", tauY, xfront );
   }
    fclose(f);
}

    Initialization

    event init (i = 0) {

    We impose boundary condition for h and \eta.

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

    We set a zero velocity at the inlet and a free outlet.

      for (vector u in ul) {
    u.n[left] = 0;
    u.n[right] = neumann(0.);
  }

    We initialize h.

      foreach(){
    zb[] =  (x<2?-(x-2)*alpha:0);
    h[] =   (fabs(x)<2)-zb[];
    }
}

    Output

    We print the elevation

    event outputfront (t += .1 ; t<=5) {
    double  xf=0,xe=0;

    tracking the front and the end of the heap

      foreach(){
   xf = h[] > 1e-10 ?  max(xf,x) :  xf ;
   xe = h[] > 1e-10 ?  min(xe,x) :  xe ;
  }

    save them

      xfront = xf;
  sprintf (s, "frontML-%.3f.txt", tauY);
  FILE * f = fopen (s, "w");  
  foreach()
    fprintf (f, "%g %g %g \n", fmin((x-xf),0), h[], xe-xf);   
  fclose(f);
  fprintf (stdout, "%g %g %g \n", t, xf, xe);
  sprintf (s, "xML-%.3f.txt", tauY);
  FILE * fp = fopen (s, "a");   
   fprintf (fp, "%g %g  \n", t, xf); 
  fclose(fp);
}

    save the hight the flux and the surface as a function of time

      event output  (t+=1){
    sprintf (s, "shapeML-%.3f.txt", tauY);
    FILE * fp = fopen (s, "a"); 
    foreach()
      fprintf (fp, "%g %g %g %g \n",x,h[],zb[],t);
      fprintf (fp, "\n");
    fclose(fp);
  }

    Results

    Run

    To run

    qcc -O2 -o bingham_muI_collapse_ML bingham_muI_collapse_ML.c
  ./bingham_muI_collapse_ML  2>log

    or with make

    make cohesive_muI_collapse_ML.tst; make cohesive_muI_collapse_ML/plots
make cohesive_muI_collapse_ML.c.html ;
   
source ../c2html.sh cohesive_muI_collapse_ML

    Comparisons RNSP - 1D

    Compare with 1D result and variable \mu

     set xlabel "x"
 set ylabel "h(x,t)"
 p[:8][:1.6]'../savagestaron/log' t '1D' w l,'shapeML-0.000.txt' u 1:($2+$3) t'2D' w l,'' u 1:3 not w l linec -1
    (script)

    (script)

    Compare Kzerswell Balmforth results \tau_Y=0 \mu=cst.

     set term png ;   set output 'bk.png';
   
 Xc=520
 L1=907-Xc
 H1=640-185
 Yc=185
 unset tics
 set key right top
 plot [0:1400]  '../Img/Balmforth_Kerswell05.png' binary filetype=png with rgbimage not,\
 '../savagestaron/log' u (Xc+($1-2)*L1*.4):(Yc+$2*(H1)) w l not,\
 'shapeML-0.00.txt' u (Xc+($1-2)*L1*.4):(Yc+($2)*(H1)) t'2D' w l
    cmp B&K 05 (script)

    cmp B&K 05 (script)

    with cohesion

    Compare with and without

     reset
 set multiplot layout 2,2
 set xlabel "x"
 set ylabel "h(x,t)"
 p[:8][:1.6]'../savagestaron/log' t '1D' w l,\
 'shapeML-0.000.txt' u 1:($2+$3) t'2D 0.0' w l,\
 'shapeML-0.100.txt' u 1:($2+$3) t'2D 0.1' w l,\
 '' u 1:3 not w l linec -1
   
 p[:8][:1.6]'../savagestaron/log' t '1D' w l,\
 'shapeML-0.000.txt' u 1:($2+$3) t'2D 0.0' w l,\
 'shapeML-0.200.txt' u 1:($2+$3) t'2D 0.2' w l,\
 '' u 1:3 not w l linec -1
   
 p[:8][:1.6]'../savagestaron/log' t '1D' w l,\
 'shapeML-0.000.txt' u 1:($2+$3) t'2D 0.0' w l,\
 'shapeML-1.000.txt' u 1:($2+$3) t'2D 1.0' w l,\
 '' u 1:3 not w l linec -1
   
 p[:8][:1.6]'../savagestaron/log' t '1D' w l,\
 'shapeML-0.000.txt' u 1:($2+$3) t'2D 0.0' w l,\
 'shapeML-0.500.txt' u 1:($2+$3) t'2D 0.5' w l,\
 '' u 1:3 not w l linec -1
   unset multiplot
    (script)

    (script)

    Compare with Balmforth results for position of front as function time

     reset
 mu=.4
 p[:5]'xML-0.000.txt' t'B=0.000 2D' w l,\
      'xML-0.001.txt' t'B=0.001 2D' w l,\
      'xML-0.010.txt' t'B=0.010 2D' w l,\
      'xML-0.020.txt' t'B=0.020 2D' w l,\
      'xML-0.050.txt' t'B=0.050 2D' w l,\
      'xML-0.100.txt' t'B=0.100 2D' w l,\
      'xML-0.200.txt' t'B=0.200 2D' w l,\
      'xML-0.300.txt' t'B=0.300 2D' w l,\
      'xML-1.000.txt' t'B=1.000 2D' w l,\
       (2+x*(2-x*mu/2))*(x<2/mu?1:NaN) t'K','../savagestaron/x.txt' t'num SH'
    compare front position (script)

    compare front position (script)

    Value of the final position (x_{front}-2) as function of cohesion. Plain line: comparison wit full 1D constant \mu http://basilisk.fr/sandbox/M1EMN/Exemples/cohesive_savagehutter.c

     set xlabel "tauY"
 set ylabel "extend"
 p'tauYxmax.txt' t '2D','../cohesive_savagehutter/tauYxmax.txt' t'1D' w l
    (script)

    (script)

    Bibliography