sandbox/M1EMN/Exemples/forchheimer.c

    Flow in a Sandglass / Hourglass with a bottom orifice

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

    Coupled trough Jackson method with a gaz with Darcy Forchheimer

    Code

    Includes and definitions

    #include "navier-stokes/centered.h"
#include "vof.h"
// Domain extent
#define LDOMAIN 4.
// heap definition
double  H0,R0,D,W,tmax,Q,Qair,Wmin,muwall,WDOMAIN,Hmask,maskr,P1,position,epai;
//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 7
// ratio of mask section
#define maskr 0.75
//#define MINLEVEL 7
char s[80];
FILE * fpf,*fwq;
scalar f[];
face vector  ug[];
scalar In[];
scalar phi[];

scalar * interfaces = {f};
face vector alphav[];
face vector muv[];
face vector av[];

//velocity of granular materials
double ugranular;

    Darcy fluid

    double betam = 100;
scalar pDarcy[],sourceDarcy[];
face vector beta[];
face vector gradpDarcy[];
mgstats mgpDarcy;

    Boundary conditions for air flow, impose pressure at the top P_1

    pDarcy[left]   = neumann(0);
pDarcy[right]  = neumann(0);
pDarcy[top]    = dirichlet(P1); 
pDarcy[bottom] = dirichlet(0);
pDarcy[right]= neumann(0);
pDarcy[left]=neumann(0);

ug.n[top]    = neumann(0);
ug.n[bottom] = neumann(0);
ug.n[right] = dirichlet(0);
ug.n[left] = dirichlet(0);

    Boundary conditions for granular flow, pressure must be zero at the surface. The pressure is zero in the hole |x|<=W/2, but the lithostatic gradient is given elsewhere on the bottom wall. No slip boundary conditions on the other walls.

    #define test 1

#if test
p[top]      = dirichlet(0);
u.n[top]    = neumann(0);
u.t[top]    = dirichlet(0);
u.t[bottom] = dirichlet(0);
u.n[bottom] = dirichlet(0);
u.n[right]  = dirichlet(0);
u.n[left]   = dirichlet(0);

/*
p[bottom]   = neumann(0);
p[right]= neumann(0);
p[left]=neumann(0);
*/

 

In[left] = neumann(0.);
In[top] = neumann(0.);
In[bottom] = neumann(0.);
In[right] = neumann(0.);

phi[left] = neumann(0.);
phi[top] = neumann(0.);
phi[bottom] = neumann(0.);
phi[right] = neumann(0.);

#else
p[top]      = dirichlet(0);
u.n[top]    = neumann(0);
u.t[bottom] = neumann(0);
u.n[bottom] = neumann(0);
p[bottom]   = dirichlet(0);
p[right]=(fabs(y)>= 0. && fabs(y)<= (LDOMAIN*position-epai)) ? dirichlet(0):  neumann(0);
p[left]=(fabs(y)>= 0. && fabs(y)<= (LDOMAIN*position-epai)) ? dirichlet(0):  neumann(0);


u.n[right] = (fabs(y)>= 0. && fabs(y)<= (LDOMAIN*position-epai)) ? neumann(0):  dirichlet(0);
u.t[right] = (fabs(y)>= 0. && fabs(y)<= (LDOMAIN*position-epai)) ? neumann(0):  dirichlet(0);
u.n[left] = (fabs(y)>= 0. && fabs(y)<= (LDOMAIN*position-epai)) ? neumann(0):  dirichlet(0);
u.t[left] = (fabs(y)>= 0. && fabs(y)<= (LDOMAIN*position-epai)) ? neumann(0):  dirichlet(0);
bid plaque;
u.t[plaque] = dirichlet(0.);
u.n[plaque] = dirichlet(0.);

p[plaque] = neumann(0.);

In[left] = neumann(0.);
In[top] = neumann(0.);
In[bottom] = neumann(0.);
In[right] = neumann(0.);

phi[left] = neumann(0.);
phi[top] = neumann(0.);
phi[bottom] = neumann(0.);
phi[right] = neumann(0.);

#endif

int main(int argc, char ** argv) {
    L0 = LDOMAIN;
    // number of grid points
    N = 1 << LEVEL;
    // maximum timestep
    DT = 0.01;
    // coefficient of friction of wall
    muwall = 0.1;
    TOLERANCE = 1e-3;
    H0 = 2.5;
    R0 = 20.000;
    // Grain size
    D = 1. / 90.;
    fwq = fopen("outWQ", "w"); // ????
    fclose(fwq);               // ????
    system("rm uprof.txt");
    Lz = LDOMAIN;
    // size of the hole
    W = 0.25;
    //Qair = 1.0;
    P1=5.;
    WDOMAIN = 2.;
    a = av;
    alpha = alphav;
    mu = muv;
    epai=0.1;
    position=1./3.;
    Q = 0;
    tmax = 2.;
    fpf = fopen("interface.txt", "w");
    run();
    fclose(fpf);
    fprintf(stdout, "\n");
    fwq = fopen("outWQ", "a");
    fprintf(fwq, " %lf %lf \n", W, Q);
    fclose(fwq);
}

#if 0
int adapt() {
#if TREE
    astats s = adapt_wavelet ({f}, (double[]){5e-3}, LEVEL, MINLEVEL);
    return s.nf;
#else // Cartesian
    return 0;
#endif
}
#endif

    initial heap, a rectangle

    // note the v
event init (t = 0) {
      
    scalar phiinit[];
    foreach_vertex()
     {
        phiinit[] = min(H0 - y, R0 - x);
    }
    fractions (phiinit, f);

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

        foreach()
     p[] = max(H0- y,0) ;
    // the boundary conditions for the pressure need to be handled by
    // the Navier--Stokes solver
    foreach()
    In[]=0;
    
    foreach()
    phi[]=0;

    initial for pressure darcy

        foreach(){
        pDarcy[] = P1;
        sourceDarcy[] = 0.;
        ug.x[]=0.;
        ug.y[]=0.;
    }
    boundary ({pDarcy, sourceDarcy});
}

    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/(\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(2.*D2)/(2.*Delta);
                In[] = D2*D/sqrt(p[]);  // this D2 is sqrt(2) D2
                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./( fm * (1. - 0. * In[]) + RHOF*(1. - fm));
    }
    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);
}

    Computation of Darcy Forchheimer pressure at each time step

    event pressureDEF (i++)
{
    double B_l=1.0;
    double B_i=1.0;
    double R_f=1.0;//4.0*1.2/2500;
    scalar un[];
    vector uold[];
    double err=1;

    Darcy Forchheimer with penalization to have a constant pressure in the gas (1-f). grains correspond to f=1.

    With $R_f= {_f} (1 + C )* (f + (1-f)) $, B_l=[ \eta \beta_l (1-\phi)] f + [\frac{1}{\beta_m}](1-f) and B_i= [\eta \beta_l (1-\phi)] f,

    and \beta_m \gg 1 the penalization parameter. We have to solve by projection \displaystyle R_f\frac{\partial \overrightarrow{u}}{\partial t}= - \overrightarrow{\nabla} p - B_l \overrightarrow{u} -B_i |\overrightarrow{u}|\overrightarrow{u} \mbox{ avec } \overrightarrow{\nabla} \cdot \overrightarrow{u} = 0

    semi implicit formulation

    \displaystyle R_f \frac{\overrightarrow{u}^{n+1} - \overrightarrow{u}^{n}} {\Delta t}= - \overrightarrow{\nabla} p^{n+1} - B_l \overrightarrow{u}^{n+1} - B_i |\overrightarrow{u}^{n}|\overrightarrow{u}^{n+1} the new velocity is \displaystyle \overrightarrow{u}^{n+1} = \frac{ R_f \overrightarrow{u}^{n}- \Delta t \overrightarrow{\nabla} p^{n+1}}{ R_f + \Delta t ( B_l + B_i |\overrightarrow{u}^{n}| )} with \displaystyle \overrightarrow{\nabla} \cdot \overrightarrow{u}^{n+1} = 0 the problem reduces to the computation of

    \displaystyle \overrightarrow{\nabla} \cdot (\beta p^{n+1} ) = S_{darcy} with \displaystyle \beta = \frac{ 1}{ R_f+ \Delta t ( B_l + B_i |\overrightarrow{u}^{n}| ) }f + (\frac{\beta_m}{(R_f+ \Delta t)})(1 - f); 

        foreach_face(){
        uold.x[] = ug.x[];
        un[] = norm(ug);
        beta.x[] =1.0/(R_f+dt*(B_l+B_i*un[]))*f[]  + (betam/(R_f+dt))*(1. - f[]);
    }
    boundary ((scalar *){beta});

    computation of source term \displaystyle \overrightarrow{u}_s = ( \frac{R_f \overrightarrow{u}^{n} }{ R_f + \Delta t ( B_l + B_i |\overrightarrow{u}^{n}| )})f+ (\frac{R_f \overrightarrow{u}^{n} }{ R_f + \Delta t })(1-f)

        vector us[];
    foreach_face(){
        us.x[]=R_f*ug.x[]/(R_f+dt*(B_l+B_i*un[]))*f[] +   R_f*ug.x[]/(R_f+dt)*(1-f[]);
    }
    boundary ((scalar *){us});

    \displaystyle S_{darcy}= \overrightarrow{\nabla} \cdot \overrightarrow{u}_s/\Delta t

        scalar div[];
    foreach(){
        div[]=0.0;
        foreach_dimension()
           div[] +=us.x[1]-us.x[];
        div[] /=dt*Delta;
        }     
    sourceDarcy=div;

    solve \nabla \cdot (\beta \nabla p_{darcy} )= S_{darcy} with Poisson solver

        mgpDarcy = poisson (pDarcy, sourceDarcy, beta);

    update \displaystyle \overrightarrow{u}^{n+1} = \frac{ R_f \overrightarrow{u}^{n}- \Delta t \overrightarrow{\nabla} p^{n+1}}{ R_f + \Delta t ( B_l + B_i |\overrightarrow{u}^{n}| )}

        foreach() {
        un[] = norm(ug);
        foreach_dimension(){
           gradpDarcy.x[] =(pDarcy[0,0] - pDarcy[-1,0])/Delta;
           ug.x[]= (R_f*uold.x[] - gradpDarcy.x[]*dt)/(R_f +dt*( B_i * un[]+B_l))*f[]+
                    (R_f/betam*uold.x[] -gradpDarcy.x[]*dt) / ((R_f +dt)/betam)*(1-f[]) ;
        }
    }
    boundary ((scalar *){ug,gradpDarcy});

    monitoring unsteadiness

       double max=0;
   foreach() {
       double du = fabs(un[] -norm(ug));
       if(du > max) max = du;
   }
        err=max;
  fprintf (stderr, "%g   du= %g\n", t, err);
  
}

event acceleration (i++)
{
    double eps = 1.0;

    Coupling gravity : Grad(p) added to gravity

        foreach_face(y)
    av.y[] += - 1. - eps*(pDarcy[0,0] - pDarcy[0,-1])/Delta;
    foreach_face(x)
    av.x[] += - eps*(pDarcy[0,0] - pDarcy[-1,0])/Delta ;
    boundary ((scalar *){av});
}

    save velocity for analytical solution

    event sauf (t += 0.1) {

  FILE * fp = fopen("uprof.txt","a");
  fprintf (fp,"%g %g   \n", t,interpolate(ug.y, 3.5, 2));
}      


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,pDarcy,gradpDarcy,In,ug}, fp, linear = true);
    fclose (fp);
}

    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;
    if (t >= 0.) fprintf (stdout,"%lf %lf %lf\n",t,V,ugranular);
    fflush (stdout);
}

    Calculate velocity of flowing materials in the bulk and add it to Qair at the top to simulate the counter flow

    event velocitybulk (i++) {
    static double VVold, VV = 3.9,Qinst = 0;
    VVold = VV;
    VV = 0;
    foreach()
    VV = VV + f[]* Delta * Delta;
    Qinst = -(VV-VVold)/dt;
    ugranular = Qinst/(LDOMAIN*(1-maskr));
}

    film output

    #if 1
event movie (t += 0.05) {
    static FILE * fp1 = popen ("ppm2mpeg > level.mpg", "w");
    scalar l[];
    foreach()
    l[] = level;
    boundary ({l});
    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[])));
    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}});
    
    foreach()
      l[] =  sqrt(sq(ug.x[]) + sq(ug.y[]));
    boundary ({l});
    static FILE * fp4 = popen ("ppm2mpeg > velo_f.mpg", "w");
    output_ppm (l, fp4, min = 0, max = 2., linear = true,
                n = 512, box = {{0,0},{Lz,Lz}});
    
    static FILE * fp3 = popen ("ppm2mpeg > pDarcy.mpg", "w");
    foreach()
    l[] = pDarcy[];
    boundary ({l});
    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


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

    Run

    to run

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

    Results

    set pm3d map
set palette rgbformulae 22,13,-31;
unset colorbox
set xlabel "x  iso p"
splot [][:] 'field-1.txt' u 1:2:10   not
reset
    pressure (script)

    pressure (script)

      set ylabel "p_darcy"
set xlabel "y"
 p 'field-1.txt' u 2:10
    pressure (script)

    pressure (script)

    Analytical solution of the problem, $ = 2 $

    \displaystyle \frac{\partial u}{\partial t}= -2 - u + u^2 \displaystyle u(t)= \frac{2 (exp(3 t )-1)}{(1+2 exp(3 t))}

    set xlabel "t"
set ylabel "u(t)"
 p'uprof.txt' t'num.' ,-2*(exp(3*x)-1)/(1+2*exp(3*x)) t'analytic'
    (script)

    (script)

    Bibliography