sandbox/M1EMN/Exemples/column_SCC.c
Collapse of a Bingham flow
Physical problem
First, we consider simple viscous collapse of a Nextonian fluid and find the t^{1/5} selfsimilar solution. To be compared with all the 1D models…
Then wew consider collapse of a Non-Nextonian fluid : the Bingham fluid.
Application to mud flows,debris flows etc application to wet concrete flow: collapse of columns (Abrahams cone test).
equations
We propose an implementation of the Bingham rheology. For those flows, when stress is larger than a yield value, the media flows. For non Newtonian fluids, the “viscosity” is at least a function of the second principal invariant of the shear strain rate tensor that we define here as (other definitions proportional to this one are possible): \displaystyle D_2=\sqrt{\sum_{i,j}D_{ij}D_{ij}} In littérature two norms are used, the Euclidian: \displaystyle ||D||=\sqrt{\frac{1}{2}\sum_{i,j}D_{ij}D_{ij}} and the Frobenius one: \displaystyle |D|=\sqrt{\sum_{i,j}D_{ij}D_{ij}} Obviously \displaystyle ||D||=\sqrt{\frac{1}{2}} |D|= \frac{D_2}{\sqrt{2}} The strain rate tensor is D_{ij}=(\partial_iu_j+\partial_ju_i)/2, it has unit of \text{s}^{-1}. the components in 2D:
D_{11}=\frac{\partial u}{\partial x}, D_{12} =\frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}), D_{21} =D_{12} =\frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}), D_{22}=\frac{\partial v}{\partial y}
And where the second invariant is D_2=\sqrt{D_{ij}D_{ij}} \displaystyle D_2^2= D_{ij}D_{ij}= D_{11}D_{11} + D_{12}D_{21} + D_{21}D_{12} + D_{22}D_{22}
hence, as D_{12}D_{21} = D_{21}D_{12}: \displaystyle D_2^2= D_{ij}D_{ij}= ( \frac{\partial u}{\partial x})^2 + (\frac{\partial v}{\partial y})^2 + \frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x})^2
We have by definition of the Bingham rheology \displaystyle \tau_{ij} = 2 \mu_1 D_{ij} + \tau_y \frac{D_{ij}}{||D||} (hence if ||\tau_{ij}|| > \tau_y there is a flow). In order to identify this as an effective viscosity, \tau_{ij} = 2 \eta_{eft} D_{ij},
\displaystyle \tau_{ij} = 2 \mu D_{ij} + \tau_y \frac{D}{||D||} = 2 ( \mu + \frac{\tau_y}{2||D||})D_{ij} then the equivalent, or effective, of apparent viscosity is with D_2= \sqrt{2}||D|| \displaystyle \mu_{eq}= \mu_1 + \frac{\tau_y}{\sqrt{2} D_2 }
We can use a general Herschel-Bulkley formulation, with here N=1, \displaystyle \tau = \tau_y + \mu_N \dot{\gamma}^{N} which is coded with an effective viscosity of the form:
\displaystyle \mu(||D||)= \mu_N(2||D||)^{N-1} + {\frac{\tau_y}{2||D||}}.
\displaystyle (|| {\dot \gamma} ||^2)^{(N-1)/2} = (\sqrt{2} D_2)^{(N-1)} = = (({2} ||D||)^{(N-1)} )
Where \tau_y is the yield stress and \mu_N the generalized viscosity.
check with bingham simple, it should be the same.
Adimensionalisation
The Navier Stokes equations with Herschel-Bulkley/ Bingham are:
\displaystyle \rho _0 \frac{d u }{d t } = - \nabla p - \rho_0 g + \nabla \cdot \tau
we use the effective viscosity \tau_{ij} = 2\mu_{eff} D_{ij} with \displaystyle 2\mu_{eff} D_{ij} = (\frac{\tau_y}{ ||D||}+ 2 \mu_N||D||^{N-1}) D_{ij} so \displaystyle \tau_{ij} = 2 (\frac{\tau_y}{2||D||}+\mu_N||D||^{N-1}) D_{ij} Using x=L \bar x, y=L \bar y, and (u,v)=U_0(\bar u,\bar v) we have p= \rho_0 U_0^2 \bar p let us write the viscous part:
\displaystyle \tau = \rho_0 U_0^2 \left[ (\dfrac{\tau_y}{ 2 \rho_0 U_0^2 || \bar D||}+\frac{\mu_N }{\rho_0 U_0L} (\frac{U_0}{L})^{N-1}||\bar D||^{N-1}) 2 {\bar D} \right] or finally \displaystyle \tau = \rho_0 U_0^2 \left[ (\dfrac{\bar \tau_y}{ 2 ||\bar D||}+\frac{1 }{Re_N}||\bar D||^{N-1})( 2 {\bar D} ) \right] with the yield stress without dimension \bar \tau_y= \dfrac{\tau_y}{ \rho_0 U_0^2} and a kind of Reynolds number Re_N=\frac{\rho_0 U_0L} {\mu_N }(\frac{U_0}{L})^{1-N}.
With these two parameters the Navier Stokes Herschell-Bulkley is: \displaystyle \frac{d \bar u }{d \bar t } = - {\bar \nabla} p - \frac{gL}{U_0^2} e_y + {\bar \nabla} \left[ (\dfrac{\bar \tau_y}{ 2 ||\bar D||}+\frac{1 }{Re_N}||\bar D||^{N-1})(2 {\bar D} ) \right]
In case of Bingham flow N=1, the Reynolds is Re_1=\frac{\rho_0 U_0L} {\mu_1 } and we define the Bingham number Bi=(\dfrac{\tau_y L}{ \mu_1 U_0}) with these two parameters the Navier Stokes Bingham is: \displaystyle \frac{d \bar u }{d \bar t } = - {\bar \nabla} p - \frac{gL}{U_0^2} e_y + \frac{1 }{Re_1} {\bar \nabla} \left[ (\dfrac{Bi}{ 2 ||\bar D||}+1)(2 {\bar D} ) \right]
This is a collapse so that it is natural to put one in front of the gravity: \frac{gL}{U_0^2}=1. This means that the scales of length and time L/T^2=g, then T=\sqrt{L/g}. The velocity scale is U_0=L/T which is \sqrt{gL}, the gravitational term is 1, the viscous term: \displaystyle \frac{1}{Re_1}= \frac{\nu}{\sqrt{g L^3}}
Code
Includes and definitions
//#include "grid/cartesian.h"
//#include "grid/multigrid.h"
#include "navier-stokes/centered.h"
#include "vof.h"
#include "heights.h"
// Domain extent
#define LDOMAIN 4.000001
//passive fluid
#define RHOF 1e-3
#define mug 1e-3
double mub,Bi;
double tmax;
double etamx;
// Maximum refinement 4./2**8 = 0.015625
// Minimum refinement 2./2**4 = 0.125
#define LEVEL 7 //7 // 9 OK
#define LEVELmin 2
scalar f[];
scalar * interfaces = {f};
scalar eta[];
Boundary conditions
u.t[bottom] = dirichlet(0);
Boundary conditions
int main() {
system("mkdir SIM");
L0 = LDOMAIN;
// number of grid points
N = 1 << LEVEL;
// maximum timestep
DT = 0.025/4 ;
TOLERANCE = 1e-3;
Numerical values
The problem is a heap of viscoplastic material of height = length released on a flat plate. We may define a cone for Abrams slup test (concrete)
// Initial conditions
#define Ltas 0.9999999999
#define Ls 0.9999999999
#define Htas 0.99999999999
The fluid has a density \rho_0, time is T=\sqrt{L/g}, We define a quantity, \mu_b which is in fact the inverse of Re_1=\frac{\rho_0 U_0L}{\mu_N } and wich is here 1. and we define the Bingham number Bi=(\frac{\tau_y L}{\mu_N U_0}), we change it.
mub= 1;
etamx=1000.;
//tmax = 199.99;
tmax = 100;
DT = 0.025;
Bi = .0;
run();
system("cp velo.mp4 veloo.mp4");
// system("cp img.png imgo.png");
mub= 1;
tmax=50;
etamx=1000.;
DT = 0.01 ;
Bi = .05*2/sqrt(2);
run();
system("cp velo.mp4 veloi.mp4");
system("cp mu.mp4 mui.mp4");
// system("cp img.png imgi.png");
mub= 1;
tmax=50;
etamx=1000.;
DT = 0.001 ;
Bi = .14*2/sqrt(2);
run();
system("cp velo.mp4 veloii.mp4");
system("cp mu.mp4 muii.mp4");
// system("cp img.png imgii.png");
}
face vector alphav[];
face vector muv[];
scalar rhov[];
event init (t = 0) {
const face vector g[] = {0.,-1.};
a = g;
alpha = alphav;
mu = muv;
rho = rhov;
Defining the initial trapozeoidal shape of the Abrams cone.
fraction (f, min(Htas - y, Ltas - (fabs(x)+ y * (Ltas-Ls))));
}
event stop (t = tmax);
total density
#define rho(f) ((f) + RHOF*(1. - (f)))
Viscosity computing: D_{ij}=(u_{i,j}+u_{j,i})/2, incompressibility D_{ii}=0. The second invariant D_2 written D2
in Gerris and Basilisk D_2=\sqrt{D_{ij}D_{ij}} in 2D: D_{ij} is D_{11}=\frac{\partial u}{\partial x}, D_{12} =\frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}), D_{21} =D_{12} =\frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}), D_{22}=\frac{\partial v}{\partial y} Then \displaystyle D_{ij}D_{ij}=
D_{11}D_{11}
+D_{12}D_{21} +
D_{21}D_{12} + D_{22} D_{22}
or: \displaystyle D_{ij}D_{ij}= ( \frac{\partial u}{\partial x})^2 + (\frac{\partial v}{\partial y})^2 + 2(\frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}))^2 or: \displaystyle D_{ij}D_{ij}=( ( \frac{\partial u}{\partial x})^2 +
(\frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}))^2)
+(
(\frac{\partial v}{\partial y})^2 + (\frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}))^2) as it is coded:
dxx=u.x[1,0] - u.x[-1,0] is (\frac{\partial u}{\partial x})(2 \Delta x)
dxy is \frac{1}{2}(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial y})(2 \Delta x) then add the second dimension.
This gives the equivalent final viscosity \displaystyle \eta_{eq}=(Bi/(\sqrt{2} D_2) + 1) \mu_b
note that we do a regularisation such as for too small D_2, the viscosity is the viscosity of a very viscous newtonian fluid \displaystyle \eta_{eq}= \eta_{max}
there are several ways to do that: \displaystyle \eta_{eq} = min ( (Bi/(\sqrt{2} D_2) + 1) \mu_b ,\eta_{max}) or \displaystyle \eta_{eq} = ( \frac{Bi}{(\sqrt{2} D_2) + \mu_1 Bi/\eta_{max} } + 1) \mu_1
event properties (i++) {
trash ({alpha});
foreach() {
eta[] = etamx;
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);
eta[] = (Bi/(sqrt(2)*D2 + mub*Bi/etamx) + 1)*mub ;
// eta[] = min(eta[],etamx);
}
}
boundary ({eta});
scalar fa[];
// filtering density twice
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);
// muv.x[] = 1./(2.*fm/(eta[] + eta[-1,0]) + (1. - fm)/mug);
alphav.x[] = 1./rho(fm);
}
foreach()
rhov[] = rho(fa[]);
boundary ({muv,alphav,rhov});
}
convergence outputs
void mg_print (mgstats mg)
{
#if 0
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));
#else
if (mg.i > 0 && mg.resa > 0.)
fprintf (stderr, "# %d %g %g %g\n", mg.i, mg.resb, mg.resa,
mg.resb > 0 ? exp (log (mg.resb/mg.resa)/mg.i) : 0);
#endif
}
convergence outputs
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);
}
some outputes for comaprison
event pij0 (t = 0.) {
char s[80];
sprintf (s, "out0-%g",Bi);
FILE * fp = fopen (s, "w");
output_facets (f, fp);
fclose (fp);
}
event pij2 (t = 2.5001234) {
char s[80];
sprintf (s, "out1-%g",Bi);
FILE * fp = fopen (s, "w");
output_facets (f, fp);
fclose (fp);
}
event pij5 (t = 4.9876) {
char s[80];
sprintf (s, "out2-%g",Bi);
FILE * fp = fopen (s, "w");
output_facets (f, fp);
fclose (fp);
}
event pij10 (t = 9.9867212){
char s[80];
sprintf (s, "out3-%g",Bi);
FILE * fp = fopen (s, "w");
output_facets (f, fp);
fclose (fp);
}
event pij40 (t = 39.9867212){
char s[80];
sprintf (s, "out4-%g",Bi);
FILE * fp = fopen (s, "w");
output_facets (f, fp);
fclose (fp);
}
event pij50 (t = 49.9867212){
char s[80];
sprintf (s, "out5-%g",Bi);
FILE * fp = fopen (s, "w");
output_facets (f, fp);
fclose (fp);
}
event pijtmax (t = tmax){
char s[80];
sprintf (s, "out6-%g",Bi);
FILE * fp = fopen (s, "w");
output_facets (f, fp);
fclose (fp);
}
#if 1
event interface (t +=1 ) {
char s[80];
sprintf (s, "SIM/field-%g", t);
FILE * fp = fopen (s, "w");
output_field ({f,p,u,uf,pf,eta}, fp, linear = true);
fclose (fp);
#endif
}
#if gfsv
event gfsview (i += 10){
static FILE * fp = popen ("gfsview2D -s column_SCC.gfv", "w");
output_gfs (fp);
}
#endif
Saving the top position and runout positionfor slump measurements
vector h[];
event timeseries (t += 0.1 ) {
heights (f, h);
double maxy = - HUGE,maxx = - HUGE;;
foreach()
if ((h.y[] != nodata) && (h.x[] != nodata)) {
double yi = y + height(h.y[])*Delta;
double xi = x + height(h.x[])*Delta;
if (yi > maxy)
maxy = yi;
if (xi > maxx)
maxx = xi;
}
char s[80];
sprintf (s, "hmax-%g",Bi);
static FILE * fp0 = fopen (s, "w");
fprintf (fp0, "%g %g %g\n", t, maxx, maxy);
fflush (fp0);
}
Movies
#if 1
event pictures (t = {0, 1. , 2., 3., 4.} ) {
scalar l[];
foreach()
l[] = f[]*p[];
boundary ({l});
if(t==0.)
output_ppm (f, file = "f0.png", min=-1, max = 2, spread = 2, n = 256, linear = true,
box = {{0,0},{4,2}});
if(t==1.)
output_ppm (f, file = "f1.png", min=-1, max = 2, spread = 2, n = 256, linear = true,
box = {{0,0},{4,2}});
if(t==2.)
output_ppm (f, file = "f2.png", min=-1, max = 2, spread = 2, n = 256, linear = true,
box = {{0,0},{4,2}});
if(t==3.)
output_ppm (f, file = "f3.png", min=0, max = 2, spread = 2, n = 256, linear = true,
box = {{0,0},{4,2}});
foreach()
l[] = level;
boundary ({l});
output_ppm (f, file = "l3.png", min=0, max = 8, spread = 2, n = 256, linear = true,
box = {{0,0},{4,2}});
}
#endif
#if 0
event movie (t += 0.05) {
scalar l[];
static FILE * fp1 = popen ("ppm2mpeg > velo.mpg", "w");
foreach()
l[] = f[]*(0.025+sqrt(sq(u.x[]) + sq(u.y[])));
output_ppm (l, fp1, min = 0, max = 0.1,
n = 1024, box = {{0,-1},{LDOMAIN-.5,1.5}});
fflush (fp1);
static FILE * fp2 = popen ("ppm2mpeg > mu.mpg", "w");
foreach()
l[] = f[]*muv.x[];
boundary ({l});
output_ppm (l, fp2, min = 0, max = 5, linear = true,
n = 1024, box = {{0,-1},{3,1.5}});
fflush (fp2);
static FILE * fp3 = popen ("ppm2mpeg > f.mpg", "w");
foreach()
l[] = f[];
output_ppm (l, fp3, min = 0, max = 2, linear = true,
n = 1024, box = {{0,-1},{LDOMAIN,1.5}});
fflush (fp3);
}
#else
event movie (t += 0.05) {
scalar l[];
foreach()
l[] = level;
boundary ({l});
output_ppm (l, file = "level.mp4", min = 0, max = LEVEL,
n = 1024, box = {{0,-1},{4,3.5}});
foreach()
l[] = f[]*(1+sqrt(sq(u.x[]) + sq(u.y[])));
boundary ({l});
output_ppm (l, file = "velo.mp4", min = 0, max = 1.5, linear = true,
n = 1024, box = {{0,-1},{4,1.5}});
foreach()
l[] = f[]*muv.x[];
boundary ({l});
output_ppm (l, file = "muv.mp4", min = 0, max = 5, linear = true,
n = 1024, box = {{0,-1},{4,1.5}});
foreach()
l[] = f[]*p[];
boundary ({l});
output_ppm (l, file = "f.mp4", min = 0, max =1, linear = true,
n = 2048, box = {{0,-1},{4,1.5}});
}
#endif
#ifdef gnuX
event output (t += .1 ) {
fprintf (stdout, "p[0:4][0:2] '-' u 1:2 not w l \n");
output_facets (f, stdout);
fprintf (stdout, "e\n\n");
}
#endif
#if QUADTREE
// if no #include "grid/multigrid.h" then adapt
event adapt (i++) {
scalar K1[],K2[];
foreach()
K1[]=(f[0,1]+f[0,-1]+f[1,0]+f[-1,0])/4*noise();
boundary({K1});
for(int k=1;k<8;k++)
{
foreach()
K2[]=(K1[0,1]+K1[0,-1]+K1[1,0]+K1[-1,0])/4;
boundary({K2});
foreach()
K1[]=K2[]*noise();
}
adapt_wavelet({K1,f},(double[]){0.001,0.01}, maxlevel = LEVEL, minlevel = LEVELmin);
//adapt_wavelet ({f,u,muv.x}, (double[]){5e-3,0.02,0.02,0.01}, LEVEL, LEVELmin,list = {p,u,pf,uf,g,f});
}
#endif
Run
to run
qcc -g -O2 -DTRASH=1 -Wall -DgnuX=1 -o column_SCC column_SCC.c -lm
qcc -g -O2 -Wall -DgnuX=1 -o column_SCC column_SCC.c -lm
./column_SCC | gnuplot
qcc -g -O2 -Wall -Dgfsv=1 -o column_SCC column_SCC.c -lm
./column_SCC
make column_SCC.tst;make column_SCC/plots
make column_SCC.c.html ; open column_SCC.c.html
Results
Exemples of viscous collapse
we first revover the Huppert collapse for newtonian fluid, Bi=0, in h \sim t^{-1/5} and x \sim t^{1/5}
reset
set logscale
set xlabel "t"
set ylabel "h_m(t),x_m(t)"
p'hmax-0' u 1:2 w l t 'xm','' u 1:3 w l t 'hm',x**(1./5) t't^{1/5}',x**(-1./5) t'1/t^{1/5}'
reset
p'out2-0' u ($1/5**(1./5)):($2*5**(1./5)) t 't=2.5' w l,\
'out3-0' u ($1/10**(1./5)):($2*10**(1./5))t 't=10' w l,\
'out4-0' u ($1/40**(1./5)):($2*40**(1./5)) t 't=40'w l,\
'out5-0' u ($1/50**(1./5)):($2*50**(1./5)) t 't=50'w l,\
'out6-0' u ($1/200**(1./5)):($2*100**(1./5)) t 't=100' w l
Exemples of Bingham collapses
We next look at Bingham cases. First the height and the run out as function of time
reset
set xlabel "t"
set ylabel "h_m(t),x_m(t)"
p'hmax-0.0707107' u 1:2 w l t 'Bi=0.07','' u 1:3 w l t 'Bi=0.07',\
'hmax-0.19799' u 1:2 w l t 'Bi=0.19799','' u 1:3 w l t 'Bi=0.19799'
Second compared to Liu et al, note that the definition of Bi contains a \sqrt{2}.
reset
set size ratio -1
unset tics
p [1:1494][1:632]'./Img/liu16.png' binary filetype=png with rgbimage not,\
'out0-0.0707107' u ($1*(376-121)+121):($2*(400-166)+166) not w l linec 1,\
'out1-0.0707107' u ($1*(376-121)+121):($2*(400-166)+166) not w l linec 1,\
'out2-0.0707107' u ($1*(376-121)+121):($2*(400-166)+166) not w l linec 1,\
'out3-0.0707107' u ($1*(376-121)+121):($2*(400-166)+166) not w l linec 1,\
'out4-0.0707107' u ($1*(376-121)+121):($2*(400-166)+166) not w l linec 1,\
'out3-0.19799' u ($1*(1082-829)+829):($2*(400-166)+166) not w l linec 1,\
'out4-0.19799' u ($1*(1082-829)+829):($2*(400-166)+166) not w l linec 1,\
'out5-0.19799' u ($1*(1082-829)+829):($2*(400-166)+166) not w l linec 1
Some Films
velocity (click on image for animation)
velocity (click on image for animation)
velocity (click on image for animation)
velocity (click on image for animation)
viscosity (click on image for animation)
level (click on image for animation)
Links
- 1D viscous collapase
- Multilayer viscous collapse
- 1D Bingham collapse
- Multilayer Bingham collapse
- granular collapse
Bibliography
Lagrée, Staron, Popinet “Scaling Laws for the Slumping of a Bingham Plastic Fluid” J. Rheol. 57, 1265 (2013); doi: 10.1122/1.4802052
Dufour and Pijaudier-Cabotz “Numerical modelling of concrete flow: homogeneous approach” Int. J. Numer. Anal. Meth. Geomech., 2005; 29:395–416
Y. Liu, N.J. Balmforth, S. Hormozi, D.R. Hewitt, Two–dimensional viscoplastic dambreaks, Journal of Non-Newtonian Fluid Mechanics, Volume 238, December 2016, Pages 65-79,
Guillaume Vinay Anthony Wachs, Jean-François Agassant, Numerical simulation of non-isothermal viscoplastic waxy crude oil flows J. Non-Newtonian Fluid Mech. 128 (2005) 144–162