# collapse of a rectangular Bingham column on a slope,

## Problem

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. We do not solve NS but we solve RNSP equations.

## RNSP

We solve the boundary layer (RNSP) equations corresponding to a thin layer approximation of the Navier Stokes equations: \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 \tau_{xy}= (\tau_y + \mu \dfrac{\partial u}{\partial z}) with the Multilayer Solver technique from Audusse Sainte Marie. We tune the solver to consider a non newtonian viscosity function of the shear \displaystyle \nu_{eq}=\nu + \frac{\tau_y/\rho}{|\frac{\partial u}{\partial z}|} the multilayer.h is changed in saint-venantNN.h (all is in saint-venantB.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 Multilayer Solver to see the influence of the inertia neglected in the 1D model.

# Code

#include "grid/cartesian1D.h"
double BBingham;
#include "saint-venantB.h"

### Non Newtonian viscosity

The definition of viscosity as a function of shear:

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

int main() {
X0 = -.5;
L0 = 1.5;
G  = 1.;
alpha = 1;
N  = 512*4;
nl = 128/2;
nu = 1.;

a loop for the three values of Bingham parameter

 for (B = 0.5 ; B <= 2 ; B += 0.75){ //  B = 1.25;  //.5 1.25 2
BBingham = B;
sprintf (s, "xML-%.2f.txt", B);
FILE * fp = fopen (s, "w");
fclose(fp);
sprintf (s, "shapeML-%.2f.txt", B);
FILE * fs = fopen (s, "w");
fclose(fs);
run();
}
}

## Initialization

event init (i = 0) {

We impose boundary condition for h and \eta.

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

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

  for (vector u in ul) {
u.n[left] = 0;
u.n[right] = neumann(0.);
}

We initialize h.

  foreach(){
h[] =   (fabs(x)<.25);
zb[] = -(x-X0)*alpha;}
}

## Output

We print the elevation

event outputfront (t += .1 , t< 5 ) {
double  xf=0,xe=0;

tracking the front and the end of the heap

  foreach(){
xf = h[] > 1e-4 ?  max(xf,x) :  xf ;
xe = h[] > 1e-4 ?  min(xe,x) :  xe ;
}

save them

  sprintf (s, "frontML-%.2f.txt", B);
FILE * f = fopen (s, "w");
foreach()
fprintf (f, "%g %g %g \n", fmin((x-xf),0), h[], xe-xf);
fclose(f);
fprintf (stdout, "%g %g %g \n", t, xf, xe);
sprintf (s, "xML-%.2f.txt", B);
FILE * fp = fopen (s, "a");
fprintf (fp, "%g %g  \n", t, xf);
fclose(fp);
}

save the hight the flux and the yield surface as a function of time

  event output  (t = {0, 0.0625, 0.25, 1, 4}){
sprintf (s, "shapeML-%.2f.txt", B);
FILE * fp = fopen (s, "a");
foreach()
fprintf (fp, "%g %g %g\n", x, h[],t );
fprintf (fp, "\n");
fclose(fp);
}

# Results

## Run

To run

qcc -O2 -o bingham_collapse_ML bingham_collapse_ML.c
./bingham_collapse_ML  2>log

or with make

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

## Comparisons RNSP - 1D

Compare with Balmforth results B=0.5 (last t=100 not plotted in 2D)

 set xlabel "x"
set ylabel "h(x,t)"
p[-0.4:1.][:1.6]'../bingham_collapse_noSV/shape-0.50.txt' t '1D' w l,'shapeML-0.50.txt' t'2D' w lp B=0.5 (script)

Compare with Balmforth results B=1.25 (last t=100 not plotted in 2D)

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

Compare with Balmforth results B=2.0 (last t=100 not plotted in 2D)

 set xlabel "x"
set ylabel "h(x,t)"
p[-0.4:1.][:1.6]'../bingham_collapse_noSV/shape-2.00.txt' t '1D' w l,'shapeML-2.00.txt' t '2D'w lp B=2 (script)

Compare with Balmforth results for position of front as function time

p[:10]'xML-0.50.txt't'B=0.5 2D' w l,'../bingham_collapse_noSV/x-0.50.txt't'B=0.5 1D' w l,'xML-1.25.txt't'B=1.25 2D' w l,'../bingham_collapse_noSV/x-1.25.txt't'B=1.25 1D' w l,'xML-2.00.txt't'B=2.0 2D' w l,'../bingham_collapse_noSV/x-2.00.txt't'B=1.25 2D' w l compare front position (script)

Montpellier 07/17