sandbox/yixian/granular_bottom_dynamicNGF.c
Flow in a Sandglass / Hourglass / Silo with a bottom orifice
We propose an implementation of the dynamic non-local granular fluidity(dynamic NGF) model proposed by Kamrin and Henann(2015) for the flow in a silo with a bottom orifice.
Code
Includes and definitions
#include "grid/multigrid.h"
#include "navier-stokes/centered.h"
#include "vof.h"
#include "poisson.h"
#include "diffusion.h"
#define LDOMAIN 1.
double H0,R0,D,W,tmax;
double zetaA; // nonlocal amplitude A
#define RHOF 1e-4
#define mug 1e-5
#define LEVEL 6
char s[80];
FILE * fpf,*fwq;
scalar f[];
scalar * interfaces = {f};
face vector alphav[];
face vector muv[];
scalar rhov[];
scalar duxx[],duxy[],duyx[],duyy[],dpdx[],dpdy[];
/* Nonlocal fields and params*/
scalar eta[];
scalar zeta_loc[]; //
face vector alphaf[]; //
scalar muIloc[]; /* local mu(I) */
scalar muI[];
scalar muIold[]; // μ^n
scalar gloc[]; /* local fluidity */
scalar gnonloc[]; /* nonlocal fluidity (solution) */
scalar gold[]; // g^n
double mu_s = 0.4;
double mu_2 = 0.68; //
double I0 = 0.4; //
double t0_g = 0.001; //
double rhos = 1.; // density of granular
double t_switch = 0.1;
/* Boundary conditions */
p[top] = dirichlet(0);
pf[top] = dirichlet(0);
u.n[top] = neumann(0);
u.t[bottom] = fabs(x-LDOMAIN/2)<= W/2 ? neumann(0): dirichlet(0);
u.n[bottom] = fabs(x-LDOMAIN/2)<= W/2 ? neumann(0): dirichlet(0);
p[bottom] = fabs(x-LDOMAIN/2)<= W/2 ? dirichlet(0): neumann(0);
pf[bottom] = fabs(x-LDOMAIN/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);
// ==== boundary condition for g ====
gnonloc[left] = dirichlet(0);
gnonloc[right] = dirichlet(0);
gnonloc[top] = neumann(0);
gnonloc[bottom] = fabs(x-LDOMAIN/2)<= W/2 ? neumann(0): dirichlet(0);
int main(int argc, char ** argv) {
disable_fpe (FE_DIVBYZERO|FE_INVALID);
L0 = LDOMAIN;
N = 1 << LEVEL;
TOLERANCE = 1e-3;
H0 = 0.9;
R0 = 20.000;
const face vector g[] = {0.,-1.,0};
a = g;
alpha = alphav;
mu = muv;
rho = rhov;
DT = 0.001;
tmax = 4.;
W = 0.25;
zetaA = 0.5;
D = 0.04375;
fpf = fopen ("interface.txt", "w");
run();
fclose (fpf);
fprintf (stdout,"\n");
return 0;
}
/* initial heap, a rectangle */
event init (t = 0) {
scalar phi[];
foreach_vertex()
phi[] = min(H0 - y, R0 - x);
fractions (phi, f);
foreach(){
p[] = (fabs(x-LDOMAIN/2)<= W/2 && fabs(y)<= .1) ? 0 : max(H0 - y,0) ;
}
}
#define rho(f) ((f) + RHOF*(1. - (f)))We use the full dynamic form of the NGF equation:
t_0\dot{g} = A^{2}d^{2}\nabla^{2}g - \Delta{\mu}\left(\frac{\mu_{s} - \mu}{\mu_{2} - \mu}\right)g - b\sqrt{\frac{\rho_{s}d^{2}}{p}}\mu g^{2}.
In this equation, t_0 is a relaxation time scale, and b = \Delta\mu / I_0 is a material constant. Following Dunatunga and Kamrin(2022), we employ an operator-splitting scheme via an intermediate variable g^*:
\frac{g^* - g^n}{\Delta t} = \frac{1}{t_0} A^2 d^2 \nabla^2 g^*_{i,j},
\frac{g^{n + 1} - g^*}{\Delta t} = \frac{1}{t_0}\left(-\Delta \mu \left(\frac{\mu_s - \mu^{n+1}_\mathrm{loc}}{\mu_2 - \mu^{n+1}_\mathrm{loc}}\right)g^{n + 1} - \frac{\Delta\mu}{I_0}\sqrt{\frac{\rho_s d^2}{p^{n + 1}}}\mu^n g^n g^{n + 1}\right).
In the diffusion step, g^* is solved implicitly. In the reaction step, we use the local friction coefficient \mu_{\mathrm{loc}}^{n+1}, whereas employs the global \mu^{n+1}.
event properties (i++) {
trash ({alphav});
foreach()
eta[] = mug;
// ===== 0. save g^n and μ^n =====
foreach () {
gold[] = gnonloc[];// g^n
muIold[] = muIloc[]; // μ^n
}
// ===== 1. Solve：μ^{n+1} and g_loc =====
foreach () {
gloc[] = 0.;
muIloc[] = 0.;
zeta_loc[] = 0.;
if (p[] > 0. && f[] > 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); //
double sD2 = sqrt(2.)*D2; // shear rate gammatot
double In = sD2*D/sqrt(p[]); // I
double mu_loc = mu_s + (mu_2 - mu_s)*In/(I0 + In);
muIloc[] = mu_loc;
gloc[] = sD2/(mu_loc); // g_loc = gammatot/μ_loc
double xi = zetaA * D;
zeta_loc[] = xi;
}
}
else {
gloc[] = 0.;
muIloc[] = 0.;
muI[] = 0.;
zeta_loc[] = 0.;
}
}
boundary ({gloc, muIloc, zeta_loc});
// ===== 2. t < t_switch：Local rheology =====
if (t < t_switch) {
foreach() {
if (p[] > 0. && f[] > 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);
double sD2 = sqrt(2.)*D2;
double mu_loc = muIloc[];
double etamin = sqrt(D*D*D);
double eta_loc = (mu_loc * p[]) / sD2;
eta_loc = max(eta_loc, etamin);
eta_loc = min(eta_loc, 100.);
eta[] = eta_loc;
gnonloc[] = gloc[];
muI[] = mu_loc;
} else {
eta[] = 100.;
gnonloc[] = 0.;
muI[] = 0.;
}
}
else {
eta[] = mug;
gnonloc[] = 0.;
muI[] = 0.;
}
}
boundary ({eta, gnonloc, muI});
}
// ============ 3. t >= t_switch：Dynamic NGF ============
if (t >= t_switch) {
foreach_face(x) {
double zf = 0.5*(zeta_loc[] + zeta_loc[-1,0]);
alphaf.x[] = sq(zf);
}
foreach_face(y) {
double zg = 0.5*(zeta_loc[] + zeta_loc[0,-1]);
alphaf.y[] = sq(zg);
}
boundary ((scalar *){alphaf});
foreach()
gnonloc[] = gloc[];
boundary ({gnonloc});
// 3.1 solve： g* - dt ∇·(ξ² ∇g*) = g^n
mgstats mg_g = diffusion (gnonloc, dt/t0_g, alphaf);
if (mg_g.i > 0 && mg_g.resa > 0.)
fprintf (stderr,
"# gnonloc diffusion: i=%d resa=%g resb=%g\n",
mg_g.i, mg_g.resa, mg_g.resb);
boundary ({gnonloc});
// --- 3.2 solve g^{n+1} ---
foreach() {
if (p[] > 0. && f[] > 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);
double sD2 = sqrt(2.)*D2;
double mu_new = muIloc[]; // μ^{n+1}
double mu_old = muIold[]; // μ^n
double g_old = gold[]; // g^n
double gstar = gnonloc[]; // g*
// splitting scheme α^{n+1} = Δμ (μ_s - μ^{n+1})/(μ_2 - μ^{n+1})
double alpha = (mu_2-mu_s)*(mu_s - mu_new)/
(mu_2 - mu_new);
// β^{n+1} = (Δμ/I0) sqrt(ρ_s d^2 / p^{n+1}) μ^n g^n
double beta = ((mu_2-mu_s)/I0)*
sqrt(rhos*D*D/p[])*mu_old*g_old;
double denom = 1. + (dt/t0_g)*(alpha + beta);
double gnew = gstar/denom;
gnonloc[] = gnew; // g^{n+1}
double mu_eff = sD2/(gnew);
muI[] = mu_eff;
double etamin = sqrt(D*D*D);
double eta_nonloc = (mu_eff * p[]) / sD2;
eta_nonloc = max(eta_nonloc, etamin);
eta_nonloc = min(eta_nonloc, 100.);
eta[] = eta_nonloc;
} else {
eta[] = 100.;
gnonloc[] = 0.;
muI[] = 0.;
}
}
else {
eta[] = mug;
gnonloc[] = 0.;
muI[] = 0.;
}
}
boundary ({gnonloc, eta, muI});
}
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);
alphav.x[] = 1.0 / rho(fm);
}
foreach()
rhov[] = rho(fa[]);
boundary ({muv, alphav, rhov});
foreach(){
duxy[] = (u.x[0,1] - u.x[0,-1])/2./Delta + (u.y[1,0] - u.y[-1,0])/2./Delta ;
duyx[] = (u.y[1,0] - u.y[-1,0])/2./Delta + (u.x[0,1] - u.x[0,-1])/2./Delta ;
duxx[] = 2 * (u.x[1,0] - u.x[-1,0])/2./Delta ;
duyy[] = 2 * (u.y[0,1] - u.y[0,-1])/2./Delta;
dpdx[] = (p[1,0] - p[-1,0])/2./Delta;
dpdy[] = (p[0,1] - p[0,-1])/2./Delta;
}
}
/* 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));
}
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);
}
/* output including gnonloc for diagnostics */
event interface (t = 2 ; t+=1. ; t<= 4.) {
#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,muIloc,muI,gloc,gnonloc,eta,duxx,duxy,duyx,duyy,dpdx,dpdy}, fp, linear = true);
fclose (fp);
}
/* debit: use explicit dx */
event debit (t += 0.05; t <= tmax ) {
static double Vold,V=1,Qinst=0;
Vold=V;
V=0;
double dx = L0 / N;
foreach()
V = V + f[] * dx * dx;
Qinst = -(V-Vold)/0.05;
if(t>=.1) fprintf (stdout,"%lf %lf %lf \n",t,V,Qinst);
fflush (stdout);
}to run
qcc -g -O2 -o granular_bottom_dynamicNGF granular_bottom_dynamicNGF.c -lm
./granular_bottom_dynamicNGF > volPlots
Plot of vertical velocity at time 2
reset
set pm3d map
set pm3d; set palette rgbformulae 22,13,-31;
splot 'field-2.txt' u 1:2:(abs($6)*$3) with pm3d notitle
set label "Vy" at graph 0.45, 0.85 font ",24" front
set terminal postscript color enhanced set output “t2vy.eps” replot
Plot of pressure at time 2
reset
set pm3d map
set pm3d; set palette rgbformulae 22,13,-31;
splot 'field-2.txt' u 1:2:(abs($4)**$3) with pm3d notitle
set label "pressure" at graph 0.45, 0.85 font ",24" front
set terminal postscript color enhanced
set output “t2p.eps” replot
Plot of viscosity at time 2
reset
set cbrange [0:2]
set pm3d map
set pm3d; set palette rgbformulae 22,13,-31;
set label "Viscosity" at graph 0.45, 0.85 font ",24" front
splot 'field-2.txt' u 1:2:(abs($14)**$3) with pm3d notitle
set terminal postscript color enhanced set output “t2eta.eps” replot
Links
https://basilisk.fr/sandbox/yixian/granular_bottom_NGF.c
https://basilisk.fr/sandbox/yixian/granular_bottom_constantNGF.c
https://basilisk.fr/sandbox/yixian/granular_bottom_dynamicNGF.c
https://basilisk.fr/sandbox/yixian/granular_bottom_linearisedconstantNGF.c
https://basilisk.fr/sandbox/yixian/granular_bottom_linearisedNGF.c
https://basilisk.fr/sandbox/yixian/granular_bottom_Local.c
https://basilisk.fr/sandbox/yixian/granular_bottom_muItheta.c
https://basilisk.fr/sandbox/yixian/granular_bottom_Igradient.c
Bibliography
Y. Zhou PhD https://theses.fr/2016AIXM4731
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-6L. 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.4757390Kamrin, K. and Koval, G. (2012). Nonlocal constitutive relation for steady granular flow. Phys. Rev. Lett., 108:178301.
Kamrin, K. and Henann, D. L. (2015). Nonlocal modeling of granular flows down inclines. Soft Matter, 11:179–185.
Dunatunga, S. and Kamrin, K. (2022). Modelling silo clogging with non-local granular rheology. J. Fluid Mech., 940:A14.
