LateralSandglassAngle40.c

Discharge of a silo through a lateral orifice with a bottom inclination
Code
Run
Bibliography

Discharge of a silo through a lateral orifice with a bottom inclination

We propose an implementation of the Jop Pouliquen Forterre µ(I) rheology for the flow in a hourglass. We propose to identify the relative role of the acceleration and the friction on the discharge flow by varying the silo geometry (bottom inclination angle), adding the influence of a moderate friction on the front and back wall of the silo, and non-slip boundary condition is chosen for the silo bottom.

Code

Includes and definitions

#include "grid/quadtree.h"
#include "navier-stokes/centered.h"
#include "vof.h"
// Domain extent
#define LDOMAIN 4.
// heap definition
double  H0, R0, D, W, tmax, Q, Wmin, muwall, WDOMAIN, Hmask, Lz, angle;

passive fluid small density to preserve 0 pressure and small viscocity

#define RHOF 1e-4
#define mug  1e-5
// Maximum refinement
#define LEVEL 8
//#define MINLEVEL 7
char s[80];
FILE * fpf, *fwq, *fqout;
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=0 and 0 < y< W, but the lithostatic gradient is given elsewhere 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[right] = dirichlet(0);
u.t[right] = dirichlet(0);
f[right] = neumann(0);
u.t[left] = (fabs(y) >= 0. &&
	fabs(y) <= (W + 0.)) ? neumann(0) : dirichlet(0);
u.n[left] = (fabs(y) >= 0.&&
	fabs(y) <= (W + 0.)) ? neumann(0) : dirichlet(0);
p[left] = (fabs(y) >= 0. &&
	fabs(y) <= (W + 0.)) ? dirichlet(0) : neumann(0);


int main(int argc, char ** argv) {
	L0 = LDOMAIN;
	// number of grid points
	N = 1 << LEVEL;
	// maximum timestep
	DT = 0.005;
	// coefficient of friction of wall
	muwall = 0.1;
	TOLERANCE = 1e-3;
	//silo height
	H0 = 3.84;
	R0 = 20.000;
	// Grain size
	D = 1. / 90.;
	// size of the hole
	W = 0.5;
	fwq = fopen("outWQ", "w");
	fclose(fwq);
	//film size
	Lz = LDOMAIN;
	//Silo thickness
	WDOMAIN = 2.0;
	const face vector g[] = { 0., -1., 0 };
	a = g;
	alpha = alphav;
	mu = muv;
	rho = rhov;
	//inclinaison angle
	angle = 40.;
	Q = 0.;
	//maximum simulation time
	tmax = 20.*0.625 / W;
	fpf = fopen("interface.txt", "w");
	fqout = fopen("qout.txt", "w");
	run();
	fclose(fqout);
	fclose(fpf);
	fprintf(stdout, "\n");
	fwq = fopen("outWQ", "a");
	fprintf(fwq, " %lf %lf \n", W, Q);
	fclose(fwq);
}

initial heap, a rectangle

// note the v
event init(t = 0) {
	mask(x > 1.*L0 / 4. ? right : none);

	scalar phi[];
	foreach_vertex()
		phi[] = min(H0 - y, R0 - x);
	fractions(phi, f);

lithostatic pressure, with a zero pressure near the hole

	foreach()
		p[] = (fabs(y - (W / 2. + 0.)) <= W / 2. &&
		fabs(x) <= .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 = \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() {
		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(2.*D2) / (2.*Delta);
				double In = D2*D / sqrt(p[]);
				double muI = .4 + .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[] = (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]) / 2.;
		muv.x[] = (fm*(eta[] + eta[-1]) / 2. + (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 += 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);
}

wall friction \displaystyle \frac{du}{dt} = \frac{-2 \mu_w p u}{W ||u||} equation with discretization: \displaystyle \frac{u^{n+1}-u^n}{\delta t} = \frac{-2 \mu_w p u^{n+1}}{W ||u^n||}

event friction(i++) {
	vector uold[], um[];
	int it = 0;
	double errmuw;

	foreach() {
		foreach_dimension(){
			uold.x[] = u.x[];
		}
	}

	do{
		it++;
		foreach() {
			double a = norm(u) == 0 ? HUGE : 1. + 2.*muwall*dt*p[] / (norm(u)*WDOMAIN);
			foreach_dimension(){
				um.x[] = u.x[];
				u.x[] = uold.x[] / a;
			}
		}
		boundary((scalar *){ u });
		errmuw = 0;
		foreach() {
			errmuw += sqrt(sq(u.x[] - um.x[]) + sq(u.y[] - um.y[]));
		}
		errmuw /= pow(2, 2 * LEVEL);
	}

convergence if |u^{m+1}-u^{m}| < 1e-6 or more 25 iterations

	while (it < 25 && errmuw>1e-6);
	fprintf(stderr, "errmuw = %g ----  %d   \n", errmuw, it);
}

non-slip boundary condition in the inclination bottom area

event unul(i++) {
	foreach() {
		foreach_dimension(){
			u.x[] = u.x[] * (y >tan(angle / 180.0*pi)*(x)+0.);
		}
	}

}

save some interfaces

event interface (t = 0; t += 1.; t <= tmax) {
#if dimension == 2
	output_facets(f, fpf);
#endif
	char s[80];
	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

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

film output

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

	foreach()
		l[] = f[] * (1 + sqrt(sq(u.x[]) + sq(u.y[])))*(y >tan(angle / 180.0*pi)*(x)+0.);
	boundary({ l });
	static FILE * fp2 = popen("ppm2mpeg > velo.mpg", "w");
	output_ppm(l, fp2, min = 0, max = 2., linear = true,
		n = 512, 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 = 512, box = { { 0, 0 }, { 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 } });
}
#endif

Run

to run

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

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, P. Ruyer, and P. Aussillous, Discharge flow of a bidisperse granular media from a silo: Discrete particle simulations Phys. Rev. E 92, 062204 (2015).

Zhou, Y., Lagree, P. Y., Popinet, S., Ruyer, P., and Aussillous, P. (2017). Experiments on, and ´ discrete and continuum simulations of, the discharge of granular media from silos with a lateral orifice. Journal of Fluid Mechanics, 829:459–485.