sandbox/M1EMN/Exemples/granular_sandglass_muw.c

    Flow in a Sandglass / Hourglass / Silo with a lateral orifice in a Hele-Shaw configuration

    We propose an implementation of the Jop Pouliquen Forterre µ(I) rheology for the flow in a silo with a lateral orifice of hight D_w. Like the case with an orifice at the bottom, we find that the flow follows the Beverloo-Hagen discharge law in a pure 2D case. For the 3D case, we use the ideas of Hele-Shaw cell, integrate across the width W_{domain}. This adds the influence of a moderate friction at the front and back wall to simulate a shallow 3D case by a 2D set of equations with an extra source term.

    Equations are incompressibility: \nabla \cdot u =0

    and momentum: \frac{du}{dt} = - \nabla p + \nabla \cdot (2 \eta D) + \rho g - \frac{2 \mu_w p\; u}{W_{domain} |u|} with the viscosity \eta = \frac{\mu(I)p}{\sqrt{2}D_2} (D, rate of strain tensor, D_2 his second invariant)

    If we use D_w the vertical size of the hole as length scale (do not confuse D_w, with D_{grain} the grain diameter and D the rate of strain tensor), the time scale is \sqrt{D_w/g}, the velocity scale is \sqrt{g D_w}, the pressure scale is \rho gD_w.

    without dimension momentum is:

    \displaystyle \frac{du}{dt} = - \nabla p + \nabla \cdot (\sqrt{2} \mu(I) \frac{D}{D_2}) - e_z - \frac{2 \mu_w D_w p\; u}{W_{domain} |u|}

    Grain diameter is in the \mu(I) function and \frac{2 \mu_w D_w }{W_{domain} } is the Hele-Shaw friction coefficient which contains the geometry and the wall friction. This is the parameter that we will change.

    Code

    Includes and definitions

    #include "grid/multigrid.h"
//#include "grid/quadtree.h"
////#include "grid/octree.h"
#include "navier-stokes/centered.h"
#include "vof.h"
// Domain extent
#define LDOMAIN 4.
// heap definition
double  H0,R0,Dgrain,DW,tmax,Q,muwall,sintheta,WDOMAIN;

    passive fluid small density to preserve 0 pressure and small viscocity

    #define RHOF 1e-4
#define mug  1e-5
// Maximum refinement
#define LEVEL 7
char s[80];
FILE * fpf,*fwq;
scalar f[];
scalar * interfaces = {f}; 
face vector alphav[];
face vector muv[];
scalar rhov[];

    Boundary conditions for granular flow, pressure must be zero at the surface. The pressure is zero in the hole x=1 and h_w < y< D_w+h_w, but the lithostatic gradient is given elswhere on the right wall. No slip boundary conditions on the other walls.

    p[top] = dirichlet(0);
u.n[top] = neumann(0);
u.t[bottom] =  dirichlet(0);
u.n[bottom] =  dirichlet(0);
u.n[left] = dirichlet(0);
u.t[left] = dirichlet(0);
f[left]= neumann(0);
double hw=0.09375;  // a given hight
u.t[right] = (fabs(y)>= hw && fabs(y)<= (DW+hw)) ? neumann(0):  dirichlet(0);
u.n[right] = (fabs(y)>= hw && fabs(y)<= (DW+hw)) ? neumann(0):  dirichlet(0);
p[right]   = (fabs(y)>= hw && fabs(y)<= (DW+hw)) ? dirichlet(0): neumann(0);

    main

    int main(int argc, char ** argv) {

  disable_fpe (FE_DIVBYZERO|FE_INVALID);

  L0 = LDOMAIN;
  // number of grid points
  N = 1 << LEVEL;
  // coefficient of friction of wall
  muwall=.1;
  TOLERANCE = 1e-3; 
  // Initial conditions a=.5
  H0=3.9;
  R0=20.000;
  // Grain size
  Dgrain=1./90.;
  
  
  const face vector g[] = {0.,-1.,0};
  a = g;
  alpha = alphav;
  mu = muv;
  rho=rhov;

  
  fwq = fopen ("outWQW", "w");
  fclose(fwq); 
    
  int iWW[11]={ 12 , 14 , 18, 20, 28, 76, 140 , 268 , 1024, 2048, 5000 }; 
  
  for (int iW=0;iW<11;iW++){ 
   Q = 0; 
   // maximum timestep DT = 0.002 for L7  DT =0.001 for L 8
   DT = 0.004;    
   tmax = 4.;  
   DW=0.4375; // size of the hole  
   DW=1; 
   WDOMAIN = (iWW[iW])*Dgrain; // width of the cell in grain size
   run(); 
   fprintf (stdout,"\n");
   fwq = fopen ("outWQW", "a");
   fprintf(fwq," %lf %lf %lf %lf %d  \n", DW, Q , WDOMAIN, sintheta, iW);
   fclose (fwq);
   }
}

    initial heap, a rectangle

    event init (t = 0) {
  // mask (x < 3.*L0/4. ? left : none);
  scalar phi[];
  foreach_vertex()
    phi[] = min(H0 - y, R0 - x);
  fractions (phi, f);

    lithostatic pressure, with a zero pressure near the hole to help

       foreach()
     p[] = (fabs(y-(DW/2.+ 4.*LDOMAIN/pow(2,LEVEL) ))<= DW/2. && fabs(x-LDOMAIN)<= .1) ?  0 : max(H0 - y,0) ;   
}

    total density

    #define rho(f) ((f) + RHOF*(1. - (f)))

    Viscosity computing D_2=D_{ij}D_{ji};

    In the pure shear flow D_{11}=D_{22}=0 et D_{12}=D_{21}=(1/2)(\partial u/\partial y), so that
    D_2=\sqrt{D_{ij}D_{ij}} =\sqrt{ 2 D_{12}^2} = \frac{\partial u}{ \sqrt{2}\partial y}. In a pure shear flow, \partial u/\partial y= \sqrt{2} D_2.

    The inertial number I is D \sqrt{2} D_2/\sqrt(p) and \mu = \mu_s+ \frac{\Delta \mu}{1+I/I_0} the viscosity is \eta = \frac{\mu(I)p}{\sqrt{2}D_2}:

    note that if \eta is too small an artificial small viscosity \rho D \sqrt{gD} is taken see Lagrée et al. 11 § 2.3.1

    event properties (i++) {
  trash ({alphav});
  scalar eta[];
  foreach() {
    eta[] = mug;
    if (p[] > 0.) {
      double D2 = 0.;
      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);  // this is D2     
  double sD2 = sqrt(2.)*D2; // this sD2 is (sqrt(2) D2)
  double In = sD2*Dgrain/sqrt(p[]);
  double muI = .4 + .28*In/(.4 + In);
  double etamin = sqrt(Dgrain*Dgrain*Dgrain);
  eta[] = max((muI*p[])/sD2, etamin);
  eta[] = min(eta[],100);      }
    }
  }
  boundary ({eta});
  scalar fa[];
  foreach()
    fa[] = (4.*f[] + 
      2.*(f[-1,0] + f[1,0] + f[0,-1] + f[0,1]) +
      f[1,1] + f[-1,1] + f[1,-1] + f[-1,-1])/16.;
  boundary ({fa});
  foreach_face() {
    double fm = (fa[] + fa[-1,0])/2.;
    muv.x[] = (fm*(eta[] + eta[-1,0])/2. + (1. - fm)*mug);
    // mu.x[] = 1./(2.*fm/(eta[] + eta[-1,0]) + (1. - fm)/mug);
    alphav.x[] = 1./rho(fm);
  }
foreach()
    rhov[] = rho(fa[]); 
 boundary ({muv,alphav,rhov});
}

    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,  
       exp (log (mg.resb/mg.resa)/mg.i));  
}

    convergence stats

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

    wall friction is introduced as a source term in the 2D equations (Hele-Shaw). We split the influence of the wall friction: \displaystyle \frac{du}{dt} = \frac{-2 \mu_w p\; u}{W_{domain} |u|} If we write velocity as u = |u| T, where the tangential vector is T=\frac{u}{|u|}. We introduce the Frenet base with T and N.

    So \frac{du}{dt} = \frac{-2 \mu_w p\; }{W_{domain} } T, but \frac{du}{dt} = \frac{d|u|}{dt} T + |u|^2 \frac{N}{R}, as \frac{dT}{ds}=\frac{N}{R}, hence there is no curvature, and the velocity remains in the same direction \displaystyle \frac{d|u|}{dt} = \frac{-2 \mu_w p }{W_{domain} } therefore the exact integration gives: \displaystyle u_i^{n+1}= u_i^n - \Delta t \frac{2 \mu_w p u_i^{n}}{W_{domain}|u^n|} as velocity decreases, it can not be less than 0, hence |u|^{n+1}=max(|u|^{n+1},0).

    event friction (i++) {
    foreach()
    {
        double m = 2.*muwall*dt*p[]/WDOMAIN;
        double U = norm(u);
        foreach_dimension()
        u.x[]= U > 0 ? max(U-m,0)*u.x[]/U : 0 ; 
    }
}

    note that for \mu_{wall}=0.1, and W_{domain} is measured in grain diameter (D=1/90) for comparison with experiments, we have 2 \mu_{wall}/W_{domain}= 18/n_{grains}, this coefficient is one for 18 grains. This explains that the computation breaks for about 13 grains.

    Rate of flowing materials across the hole

    event debit (t += 0.05; t <= tmax ) {
  static double Vold,V=1,Qinst=0;
  Vold=V; 
  V=0;
  foreach()
    V = V + f[]* Delta * Delta;
  Qinst = -(V-Vold)/0.05;  
  if(Qinst > Q) Q = Qinst;
   double ux=interpolate (u.x, LDOMAIN-0.01, hw+DW/2);
   double uy=interpolate (u.y, LDOMAIN-0.01, hw+DW/2);; 
   double U;
   U=sqrt(sq(ux)+sq(uy));
   sintheta = (U>0 ? fabs(ux/U) : 0);
  if(t>=.1) fprintf (stdout,"%lf %lf %lf %lf %lf \n",t,V/L0/H0,DW,Q,sintheta); 
  fflush (stdout);
}

    Run

    to run

    qcc  -g -O2 -o granular_sandglass_muw granular_sandglass_muw.c -lm 
./granular_sandglass_muw > out

    Plot of the curve showing the influence of D_w/W_{domain}

    set logscale
set xlabel "D/W"
set ylabel "Q/W^{5/2}"
set key top left  
p [:10][:100] 'outWQW' u ($1/$3):($2/($3**(1.5))) t'calc'   w lp, .76*x**1.5*sqrt(1/(1+0.26*x)) t'fit article .76*x**1.5*sqrt(1/(1+0.26*x))', 1.49*x linec -1 , .76*x**1.5 linec -1
    extended Beverloo (script)

    extended Beverloo (script)

    Related examples

    Bibliography

    • Y. Zhou Ejection de gaz et de grains suite à la rupture d’un crayon de combustible nucléaire : modélisation de la dynamique, thèse soutenue le 2 Novembre 2016

    • L. Staron, P.-Y. Lagrée, & S. Popinet (2014) “Continuum simulation of the discharge of the granular silo, A validation test for the μ(I) visco-plastic flow law” Eur. Phys. J. E (2014) 37: 5 DOI 10.1140/epje/i2014-14005-6

    • L. Staron, P.-Y. Lagrée & S. Popinet (2012) “The granular silo as a continuum plastic flow: the hour-glass vs the clepsydra” Phys. Fluids 24, 103301 (2012); doi: 10.1063/1.4757390

    • Y. Zhou, P.-Y. Lagrée, S. Popinet, P. Ruyer, P. Aussillous (2017) “Experiments on, and discrete and continuum simulations of, the discharge of granular media from silos with a lateral orifice”, Journal of Fluid Mechanic https://doi.org/10.1017/jfm.2017.543