sandbox/M1EMN/Exemples/bingham_collapse_NH.c

    collapse of a rectangular Bingham column on a slope (Eugène and Eugène)

    Problem of Liu and Mei/ Balmforth: collapse of Bingham fluid

    A heap of fluid following Bingham rheology is released along a constant slope. We see the front moving to the right, and the left front going slowly up hill to the left. A real Bingham flow should stop. See Balmforth 1D related examples with only mass equation and lubrication

    animation of the collapse

    animation of the collapse

    We solve here the simplified system of RNSP equations with the “vertically-Lagrangian multilayer solver for free-surface flows” (not Navier Stokes and not Multilayer).

    RNSP Model equations

    We solve the boundary layer Prandtl (RNSP) equations corresponding to a thin layer approximation of the Navier Stokes equations (hydrostatic): \displaystyle \left\{ \begin{aligned} \dfrac{\partial u}{\partial x} + \dfrac{\partial v}{\partial y} &= 0 \\ \rho (\dfrac{\partial u}{\partial t} + \dfrac{\partial u^2}{\partial x} + \dfrac{\partial u v}{\partial y} )&= - \rho g Z'_b - \dfrac{\partial p}{\partial x} + \dfrac{\partial \tau_{xy}}{\partial y} \\ 0 &= -\rho g -\dfrac{\partial p}{\partial y}. \end{aligned} \right. with Bingham rheology \tau_{xy}= \tau_y + \mu \dfrac{\partial u}{\partial z}. We define an equivalent viscosity \displaystyle \tau_{xy}= \rho \nu_{eq} \dfrac{\partial u}{\partial z} \text{ with } \nu_{eq}=\rho \nu (1 + \dfrac{\tau_y/{\rho \nu}}{|\dfrac{\partial u}{\partial z}|}) the layered/hydro.h is changed in http://basilisk.fr/sandbox/M1EMN/Exemples/hydroNN.h

    If we integrate over the depth of flow h the incompressiblility equation, we obtain \frac{\partial h}{\partial t} + \frac{\partial Q}{\partial x}=0, where Q=\int_0^h udy. If we neglect inertia in the momentum equation, we solve for u in \displaystyle 0 = - \rho g Z'_b - \dfrac{\partial p}{\partial x} + \dfrac{\partial \tau_{xy}}{\partial y} were the pressure is hydrostatic p=\rho g (h-y), and with the Bingham rheology. This gives then Q=\int_0^h udy that we put in the mass conservation, and hence we obtain the 1D kinetic wave from the paper of Balmforth. Here we use the Hydro Solver to see the influence of the inertia neglected in the 1D model.

    The example is linked to the collapse of a viscous fluid (Huppert 82 “The propagation of two-dimensional and axisymmetric viscous gravity currents over a rigid horizontal surface”), here we have a slope. It is done with only mass equation and lubrication here and with shallow water there. The same RNSP equations are solved with Multilayer shallow water (using the Multilayer Shallow Water (Saint Venant Multi Couches) strategy of Audusse Sainte-Marie et al 2011. See De Vita 2020 for details). We compare here with Popinet (2019) both resolutions of RNSP system. The hydrostatic http://basilisk.fr/src/layered/hydro.h “hydro.h” has been changed in “hydroNN.h” to include variable viscosity. “difusion.h” has been chenged in it and is now “difusionNN.h” (Non-Newtonian).

    it is possible to test the non hydrostatic http://basilisk.fr/src/layered/nh.h, it should be done next.

    Snapshot of the heap.

    Snapshot of the heap.

    Code

    Includes

    #include "grid/multigrid1D.h"
double BBingham;
#include "hydroNN.h"
//#include "layered/nh.h"

    Non Newtonian viscosity

    The definition of viscosity as a function of shear (computed in “difusion.h”):

    double nu_eq(double shear,double pipi){
    double nu_eq;
    nu_eq = nu*(1 + BBingham/sqrt(sq(shear) + sq(.1e-8)));
    return nu_eq;}

    definition of viscosity (it is a function of shear: Bingham flow) \displaystyle \tau = \nu (\frac{\partial u}{\partial z} + B) where \nu is the viscosity and \nu B is the yeld stress, we define an equivalent viscosity: \displaystyle \nu_{eq}=\nu (1 + \frac{B}{|\frac{\partial u}{\partial z}|}) Note the regularization introduced to avoid division by zero. The \nu and the B have no dimension….

    Main and BC

    double tmax,alpha;


int main() {
    X0 = -0.5;
    L0 = 1.5;
    G  = 1.;
    N  = 2048/2;
    nl = 128*2;
    nu = 1.;
    alpha=1;
    tmax = 1;
    BBingham = 1.25 ;
    run();
}

    We impose boundary condition for the position of the surface \eta (as h not defined in ‘hydro.h’, but by construction \eta= z_b+h).

    eta[left] = neumann (0);
eta[right] = neumann (0);

    Initialization

    event init (i = 0) {

    We set a zero velocity at the inlet and a free outlet.

            u.n[left] = dirichlet(0);
        u.n[right] = neumann(0.);

    We initialize h.

        foreach() {
        zb[] = -(x)*alpha;
        double h0 = (fabs(x)<.25) + dry;
        foreach_layer()
        h[] = (h0 )/nl;}
}

    Output

    We use gnuplot to visualise the profile during time in live with X11 or we generate an animatioingif`

    animation of the collapse

    animation of the collapse

    #if 0
void plot_profile (double t, FILE * fp)
{
    fprintf (fp,
             "set title 't = %.2f'\n"
             " h(x) = (0.9*(1.28338-x*x))**(1/3.) \n"
             "t = %.2f'\n"
             "p [0:2][0:1.5]'-' u ($1/(t**.2)):($2*(t**.2)) w lp lc 3 t 'num' , h(x) t'huppert' \n", t,t);
    foreach()
    fprintf (fp, "%g %g \n", x,eta[] );
    fprintf (fp, "e\n\n");
    fflush (fp);
}
#else
void plot_profile (double t, FILE * fp)
{
    fprintf (fp,
             "set title 't = %.2f'\n"
             "p [-.5:1. ][:1.25]'-' u 1:3:2 w filledcu lc 3 t 'surface', '../bingham_collapse_noSV/shape-1.25.txt' t 'B=1.25 h' w l\n", t);
    foreach()
    fprintf (fp, "%g %g %g\n", x,eta[]-zb[],0.00);
    fprintf (fp, "e\n\n");
    fflush (fp);
}
#endif

event profiles (t += .01) {
    static FILE * fp = popen ("gnuplot --persist 2> /dev/null", "w");
     if(t==0) fprintf (fp,"set term gif animate;set output 'animate.gif';set size ratio .333333\n");
    // comment to have direct X11 animation, uncomment to generate the gif
    plot_profile (t, fp);
    
    if(t==tmax){
        fprintf (fp,"! cp animate.gif ../animate.gif\n");
        fprintf (fp,"! cp animate.gif a2.gif \n");
    }
}

    generate a snapshot at t=t_{max}.

    event gnuplot (t = end) {
    FILE * fp = popen ("gnuplot", "w");
    fprintf (fp,
             "set term png enhanced size 640,200 font \",8\"\n"
             "set output 'snapshot.png'\n");
    plot_profile (t, fp);
    //foreach()
    // fprintf (stdout, "%g %g %g \n", x,eta[],t);
    //getchar();
}

    We print the elevation h=\eta -z_b at last time.

    event output (t = tmax) {
    foreach()
    fprintf (stdout, "%g %g %g \n", x, eta[]-zb[],t );
    fprintf (stdout, "\n");
}

    Compilation

    Usual make file

     make bingham_collapse_NH.tst
 make bingham_collapse_NH/plots;make bingham_collapse_NH.c.html
 source ../c2html.sh  bingham_collapse_NH

    Results

    Comparison with 1D Balmforth with B=1.25 at time t=1

     set xlabel "x"
 set ylabel "h(x,t)"
 p [-.4:.5][0:1.2]'out' w lp t'NH',\
'../bingham_collapse_noSV/shape-1.25.txt' t 'B=1.25 h' w l
    Comparison (script)

    Comparison (script)

    Comparison with 1D Balmforth and with 2D RNSP Multilayer multilayer.h with B=1.25 at time t=1

    set xlabel "x"
set ylabel "h(x,t)"
p[-0.4:.5][:1.2]'out' w lp t'2D NH',\
'../bingham_collapse_noSV/shape-1.25.txt' t '1D   ' w l,\
'../bingham_collapse_ML/shapeML-1.25.txt' t'2D ML' w l
    B=1.25 (script)

    B=1.25 (script)

    confine 03/20 (instead of using alcoholic gel, use hydro include file!)

    Bibliography