# Flow in a Sandglass / Hourglass

## Problem

A silo is a structure for storing bulk materials like grains, coal, cement, gravels, woodchips, food products… An hourglass or sandglass (marine term), or sand timer, or sand clock, or even egg timer, is a device used to measure the passage of time.

In both cases, a granular material (sand, grains…) is going out from a reservoir through a small aperture.

The fact that flow rate is constant has been used for time measurement for a long time. Hagen (the same Hagen from Hagen-Poiseuille) documented this and found the scaling law with aperture size,this was rediscoverd by Beverloo in the 60s.

 set arrow 1 from 4,0 to 4,1.1 front
set arrow 2 from 4,1 to 4,0 front
set arrow 3 from 0.1,.2 to 4.9,0.2 front
set arrow 4 from 4.9,.2 to 0.1,0.2 front
set arrow 5 from 2,.1 to 3.,0.1 front
set arrow 6 from 3,.1 to 2.,0.1 front
set label 1 "L0" at 3,.25 front
set label 2 "2 W" at 2.5,.05 front
set label 3 "height \nof grains" at 4.1,.5 front
set label 4 "grains" at 1.2,.9 front
set label 5 "air" at 1.2,1.2 front
set xlabel "x"
set ylabel "y"
p [0:5]0 not,1+.05*(x-2.5)**2 w filledcurves x1 linec 3 t'free surface'
unset arrow 1; unset arrow 2; unset arrow 3
unset arrow 4; unset arrow 5; unset arrow 6
unset label 1; unset label 2; unset label 3
unset label 4; unset label 5 silo (script)

## Equations

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

We find that the flow through the orifice follows the Beverloo-Hagen discharge law.

## Code

Includes and definitions

#include "grid/multigrid.h"
#include "granular.h"
// Domain extent
#define LDOMAIN 5.
// heap definition 2W is the size of the gate of the silo
double  H0,R0,D,W,tmax,Q,Wmin,DW;
//film size
double  Lz;
// Maximum refinement
#define LEVEL 6
///6 here to speed up for web site but 7 better
char s;
FILE * fpf,*fwq;

Boundary conditions for granular flow, pressure must be zero at the surface. The pressure is zero in the hole (of width 2W), but the lithostatic gradient is given elswhere on the bottom. No slip boundary conditions on the sides.

p[top] = dirichlet(0);
u.n[top] = neumann(0);
u.t[bottom] =  fabs(x-LDOMAIN/2)<= W ? neumann(0):  dirichlet(0);
u.n[bottom] =  fabs(x-LDOMAIN/2)<= W ? neumann(0):  dirichlet(0);
p[bottom]   =  fabs(x-LDOMAIN/2)<= W ? dirichlet(0): neumann(0);
u.n[right] = dirichlet(0);
u.n[left] = dirichlet(0);
u.t[right] = dirichlet(0);
u.t[left] = dirichlet(0);
f[left]= neumann(0);

the three cases in the main

int main() {
L0 = LDOMAIN;
// number of grid points
N = 1 << LEVEL;
// maximum timestep
DT = 0.01;
TOLERANCE = 1e-3;
// Initial conditions a=.5
H0=4.5;
R0=20.000;
// size of the hole
Wmin = 0.157;
DW = LDOMAIN/N;

const face vector g[] = {0.,-1.};
a = g;

fwq = fopen ("outWQ", "w");
fclose(fwq);
Lz = LDOMAIN;
tmax = 10.;

#if 1
for(W = Wmin;  W <= .6 ; W+= 2*DW ){
Q = 0;
tmax = 4.;
fpf = fopen ("interface.txt", "w");
run();
fclose (fpf);
fprintf (stdout,"\n");
fwq = fopen ("outWQ", "a");
fprintf(fwq," %lf %lf \n", W, Q);
fclose (fwq);
}
#endif
Q = 0;
tmax = 20;
W = .5;
fpf = fopen ("interface.txt", "w");
run();
fclose (fpf);
}

initial heap, a rectangle

event init (t = 0) {
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(x-LDOMAIN/2)<= W && fabs(y)<= .1) ?  0 : max(H0 - y,0) ;
}

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

save some interfaces

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

Rate of flowing materials across the hole \displaystyle -\frac{dV}{dt} \text{ with } V = \int f dv

event debit (t += 0.05 ) {
static double Vold,V=1,Qinst=0;
Vold=V;
V=0;
foreach()
V = V + f[]* Delta * Delta;
Qinst = -(V-Vold)/.05;
if(Qinst > Q) Q = Qinst;
if(t>=.1) fprintf (stdout,"%lf %lf %lf %lf \n",t,V/L0/H0,W,Q);
fflush (stdout);
}

film output

#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 = 2048, 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 = 2048, box = {{0,0},{Lz,Lz}});

static FILE * fp3 = popen ("ppm2mpeg > pressure.mpg", "w");
foreach()
l[] = f[]*p[];
output_ppm (l, fp3, min = 0, linear = true,
n = 2048, box = {{0,0},{Lz,Lz}});
}
#else

event pictures (t=3) {
scalar l[];
foreach()
l[] = f[]*(1+sqrt(sq(u.x[]) + sq(u.y[])));
output_ppm (l, file = "f.png", min=0, max = 2,  spread = 2, n = 256, linear = true,
box = {{0,0},{Lz,Lz}});
}

event movie (t += 0.05) {
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});
output_ppm (l, file = "velo.mp4", min = 0, max = 2., linear = true,
n = 2048, box = {{0,-1},{Lz,Lz}});
foreach()
l[] = f[]*p[];
output_ppm (l, file = "pressure.mp4", min = 0,max = 2., 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 = 512, linear = true, box = {{0,0},{2,2}}); output_ppm (p, file = “p.png”, min=0, max = 1, spread = 2, n = 512, linear = true, box = {{0,0},{2,2}}); }

#endif

If gfsview is installed on your system you can use this to visualise the simulation as it runs. If you manage to install Bviewit is better !

#if 0
event gfsview (i += 10) {
static FILE * fp = popen ("gfsview2D -s  ", "w");
output_gfs (fp, t = t);
}
#endif

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

#if 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 without make

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

to run with make

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

## Results

Volume as function of time for different orifices

 set xlabel "t"
set ylabel "V(t)"
set key right
p[0:]'out' t 'V' w l Beverloo (script)

Fitting the Beverloo law, 2W is the size of the hole of the silo

 set xlabel "W"
set ylabel "Debit"
f(x) = a*(x +b)**1.5
fit f(x) 'outWQ' via a,b
set key left
p[0:]'outWQ' t'debit' w lp,4.254*x**(1.5),f(x) t'fit' ,'../REFCASES/outWQ_granular_sandglass.dat' w p t'for control' Beverloo (script)

note that the '../REFCASES/outWQ_granular_sandglass.dat' was computed with #define RHOF 1e-4 #define mug 1e-5 and with LEVEL 7

Plot of the interface every at time
 reset
p[0:][0:]'interface.txt' not w l interface of grain during the flow t=0, 1, 2, 3 … (script)

Plot of pressure at time 4

reset
set pm3d; set palette rgbformulae 22,13,-31;unset surface;
set ticslevel 0;
unset border;
unset xtics;
unset ytics;
unset ztics;
unset colorbox;
#set xrange[-3:3];set yrange[-3:3];
set view 0,0
sp'field-4.txt' u 1:2:($3>.9 ?$4 :0) not pressure (script)

Plot of velocity at time 4

reset
set pm3d; set palette rgbformulae 22,13,-31;unset surface;
set ticslevel 0;
unset border;
unset xtics;
unset ytics;
unset ztics;
unset colorbox;
#set xrange[-3:3];set yrange[-3:3];
set view 0,0
sp'field-1.txt' u 1:2:($3>.9 ? sqrt($7*$7+$6*$6) :0) not velocity (script) Film of pressure velocity (click on image for animation) Film of velocity velocity (click on image for animation) if quadtree level of mesh (click on image for animation) reset set border 4095 lt -1 lw 1.000 set format cb "%.01t*10^{%T}" set samples 31, 31 set isosamples 31, 31 unset surface set ticslevel 0 set xlabel "X" set xrange [ : ] noreverse nowriteback set ylabel "Y" set yrange [ : ] noreverse nowriteback set cblabel "the colour gradient" set pm3d set palette positive nops_allcF maxcolors 0 gamma 1.5 gray set view 0,0 sp'field-3.txt' u 1:2:($3>.9 ? sqrt($7*$7+$6*$6) :0) not velocity (script)

## Bibliography

• 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

• 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

• Y. Zhou P.-Y. Lagrée S. Popinet, P. Ruyer and P. Aussillous “Experiments on, and discrete and continuum simulations of, the discharge of granular media from silos with a lateral orifice” Journal of Fluid Mechanics Volume 829 25 October 2017 , pp. 459-485 http://doi.org/10.1017/jfm.2017.543

• Y. Zhou, P.-Y. Lagrée, S. Popinet, P. Ruyer, and P. Aussillous “Gas-assisted discharge flow of granular media from silos” Phys. Rev. Fluids 4, 124305 – Published 18 December 2019 DOI: 10.1103/PhysRevFluids.4.124305

• Luke Fullard, Daniel J. Holland, Petrik Galvosas, Clive Davies, P.-Y. Lagrée, and Stéphane Popinet (2019) “Quantifying silo flow using MRI velocimetry for testing granular flow models” Phys. Rev. Fluids 4, 074302, DOI: 10.1103/PhysRevFluids.4.074302