Sandglass3D.c

Flow in a Sandglass / Hourglass with an orifice at the bottom in 3D case
- Code
- Result
Bibliography

Flow in a Sandglass / Hourglass with an orifice at the bottom in 3D case

We propose an implementation of the Jop Pouliquen Forterre µ(I) rheology for the flow in a hourglass.

Code

Includes and definitions

// #include "grid/multigrid.h"
#include "grid/octree.h"
#include "navier-stokes/centered.h"
#include "vof.h"
// Domain extent
#define LDOMAIN 1.
// heap definition
double  H0,R0,D,W,tmax,Q,Wmin,DW;
//film size
double  Lz;

passive fluid small density to preserve 0 pressure and small viscocity

#define RHOF 1e-4
#define mug  1e-5
// Maximum refinement
#define LEVEL 5 //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 4\deltax < y< W+4\deltax, 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] =  (x-LDOMAIN/2)*(x-LDOMAIN/2)+ (z-LDOMAIN/2)*(z-LDOMAIN/2) <= W/2.*W/2. ? neumann(0):  dirichlet(0);
u.n[bottom] =  (x-LDOMAIN/2)*(x-LDOMAIN/2)+ (z-LDOMAIN/2)*(z-LDOMAIN/2) <= W/2.*W/2. ? neumann(0):  dirichlet(0);
//u.n[bottom] = dirichlet(0);
p[bottom]   =  (x-LDOMAIN/2)*(x-LDOMAIN/2)+ (z-LDOMAIN/2)*(z-LDOMAIN/2) <= W/2.*W/2. ? dirichlet(0): neumann(0); 
u.n[right] = dirichlet(0);
u.n[left] = dirichlet(0);
u.t[right] = dirichlet(0);
u.t[left] = dirichlet(0);
u.n[front] = dirichlet(0);
u.t[front] = dirichlet(0);
u.n[back] = dirichlet(0);
u.t[back] = dirichlet(0);
f[left]= neumann(0);

the three cases in the main

int main(int argc, char ** argv) {
  L0 = LDOMAIN;
  // number of grid points
  N = 1 << LEVEL;
  // maximum timestep
  DT = 0.001;
 
  TOLERANCE = 1e-3; 
  // Initial conditions a=.5
  H0=0.96;
  R0=20.000;
  // Grain size
  D=1./90.;
  // size of the hole
  W = 0.25;
  //Wmin = 0.15625;
  DW = LDOMAIN/N;
  fwq = fopen ("outWQ", "w");
  fclose(fwq);
  Lz = LDOMAIN;
  const face vector g[] = {0.,-1.,0};
  a = g;
  alpha = alphav;
  mu = muv;
  rho = rhov;

  
  Q = 0; 
  tmax = 3;
  fpf = fopen ("interface.txt", "w");
  run();
  fclose (fpf);
  fprintf (stdout,"\n");
  fwq = fopen ("outWQ", "a");
  fprintf(fwq," %lf %lf \n", W, Q);
  fclose (fwq);
  
}

initial heap, a paralepipede

// note the v
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[] = ((x-LDOMAIN/2)*(x-LDOMAIN/2)+ (z-LDOMAIN/2)*(z-LDOMAIN/2) <= W/2.*W/2. && fabs(y)<= .1) ?  0 : max(H0 - y,0) ;
  boundary ({p});
}

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 = \mu(I)p/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_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); }

  foreach() {
    eta[] = mug;
    if (p[] > 0.) {
      double D2 = 0.;
      double dxx = u.x[1,0,0] - u.x[-1,0,0];
      double dyy =  u.y[0,1,0] - u.y[0,-1,0];
      double dzz =  u.z[0,0,1] - u.z[0,0,-1];
      double dxy = (u.x[0,1,0] - u.x[0,-1,0] + u.y[1,0,0] - u.y[-1,0,0])/2.;
      double dxz = (u.x[0,0,1] - u.x[0,0,-1] + u.z[1,0,0] - u.z[-1,0,0])/2.;
      double dyz = (u.y[0,0,1] - u.y[0,0,-1] + u.z[0,1,0] - u.z[0,-1,0])/2.;
      D2 = sq(dxx) + sq(dyy) +sq(dzz) +2.*sq(dxy) +2.*sq(dxz) +2.*sq(dyz);
      if (D2 > 0.) {
	D2 = sqrt(2.*D2)/(2.*Delta);
	double In = D2*D/sqrt(p[]);
	double muI = .32 + .28*In/(.4 + In);
	double etamin = sqrt(D*D*D);
	eta[] = max((muI*p[])/D2, etamin);
	eta[] = min(eta[],100);
      }
    }
  }
  boundary ({eta});
  scalar fa[];
  foreach()
    fa[] = (8.*f[] + 
	    4.*(f[-1,0,0] + f[1,0,0] + f[0,-1,0] + f[0,1,0]+ f[0,0,1]+ f[0,0,-1]) +
	        2*(f[1,1,0] + f[-1,1,0] + f[1,-1,0] + f[-1,-1,0]+f[0,1,1] + f[0,-1,1] +
                   f[0,1,-1] + f[0,-1,-1]+f[1,0,1] + f[-1,0,1] + f[1,0,-1] + f[-1,0,-1])+ 
          f[1,1,1]+ f[1,1,-1]+f[-1,1,-1]+f[-1,1,1]+f[1,-1,1]+ f[1,-1,-1]+f[-1,-1,-1]+f[-1,-1,1])/64.;
  boundary ({fa});
  foreach_face() {
    double fm = (fa[] + fa[-1])/2.;
    muv.x[] = (fm*(eta[] + eta[-1])/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});

  //boundary ((scalar *){muv,alphav});
}

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 += 1) {
  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);
}

Rate of flowing materials across the hole

event debit (t = 0 ;t += 0.1 ; t<= tmax  ) {
  static double V=1; 
  V=0;
  foreach()
    V = V + f[]* Delta * Delta* Delta; 
  if(t>=0.) fprintf (stdout,"%lf %lf \n",t,V);
  fflush (stdout);
}

film output for 2D

#if 0
event movie (t += 0.05) {
  static FILE * fp1 = popen ("ppm2mpeg > level.mpg", "w");
  scalar l[];
  foreach()
    l[] = level;
  output_ppm (l, fp1, min = 0, max = LEVEL, 
	      n = 256, box = {{0,0},{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 = 2., linear = true, 
	      n = 256, box = {{0,0},{Lz,Lz}});

  static FILE * fp3 = popen ("ppm2mpeg > f.mpg", "w");
  foreach()
    l[] = f[]*p[];
  output_ppm (l, fp3, min = 0, linear = true,
	      n = 256, box = {{0,0},{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},{2,2}});
}
#endif

gfsview…

#if 0
event gfsview (i += 1) {
  static FILE * fp = popen (dimension == 2 ?
			    "gfsview2D gfsview.gfv" :
			    "gfsview3D gfsview3D.gfv",
			    "w");
  output_gfs (fp, t = t);
}
#endif

if “grid/multigrid.h” is not included, then this is the quadtree with possibility of adaptation

#if 0 // QUADTREE
event adapt (i++) {
adapt_wavelet ({f,u.x,u.y}, (double[]){5e-3,0.001,0.001}, LEVEL, LEVEL-2,
	       list = {p,u,pf,uf,g,f});
}
#endif

# Run

to run

   qcc -g -O2 -Wall -o Sandglass3D Sandglass3D.c -lm
   ./Sandglass3D > out

Result

We plot the volume as function of time. it is linear after a transition, this is Beverloo law:

p'out' w l

V(t) (script)

Bibliography

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 PhD