lid-bingham.c

Lid-driven cavity of a Bingham fluid at Re=0
- Results

Lid-driven cavity of a Bingham fluid at Re=0

We use the multigrid implementation (rather than the default tree implementation) and the yield-stress fluid solver.

#include "grid/multigrid.h"
#include "navier-stokes/centered.h"
#include "yield.h"

double tau_y = 0;

static void run_tauy()
{
  const face vector tauc[] = {tau_y,tau_y};
  tau_0 = tauc;
  run();  
}

int main()
{ 
  init_grid (128);
  // viscosity
  double mu0 = 1.;
  const face vector muc[] = {mu0,mu0};
  mu = muc;
  // Stokes flow
  stokes = true;
  // the "r" parameter from the augmented Lagrangian formulation
  r = 10.*mu0;

We vary the yield-stress and run a simulation for each value.

  tau_y = 0; run_tauy();
  tau_y = 1./sqrt(2); run_tauy();
  tau_y = 10./sqrt(2); run_tauy();
  tau_y = 100./sqrt(2); run_tauy();
  tau_y = 1000./sqrt(2); run_tauy();
    
#if 0 // uncomment this to get more points for max(psi) vs tau_0
  for (tau_y = 0.1; tau_y <= 1000; tau_y *= 1.5)
    run_tauy ();
#endif
}

The default boundary conditions are symmetry (i.e. slip walls). We need no-slip on three boundaries and u=1 on the top boundary i.e.

u.t[top] = dirichlet(1);

For the other no-slip boundaries this gives

u.t[bottom] = dirichlet(0);
u.t[left]   = dirichlet(0);
u.t[right]  = dirichlet(0);

For the colocated solver, imposing boundary conditions for the normal components of the (face-centered) advection velocity improves the results. Ideally, this should be done automatically by the solver.

uf.n[left]   = 0;
uf.n[right]  = 0;
uf.n[top]    = 0;
uf.n[bottom] = 0;

We want the simulation to stop when we are close to steady state. To do this we store the u.x field of the previous timestep in an auxilliary variable un.

scalar un[];

Every 0.1 time units we check whether u has changed by more than 10^-5. If it has not, the event returns 1 which stops the simulation. We also output the evolution of the difference and multigrid iterations counts on standard output.

event logfile (t += 0.1; i <= 10000) {
  double du = change (u.x, un);
  if (i > 0 && du < 1e-5)
    return 1; /* stop */
  fprintf (stdout, "du %f %d %g %d %d\n", t, i, du, mgp.i, mgu.i);
  fflush (stdout);
}

This event will happen after completion of the simulation. We compute the streamfunction \psi and the unyielded field and output that to a file.

event streamfunction (t = end)
{
  scalar psi[], omega[];

  psi[top] = dirichlet(0);
  psi[bottom] = dirichlet(0);
  psi[left]   = dirichlet(0);
  psi[right]  = dirichlet(0);
  
  foreach() {
    omega[] = (u.y[1] - u.y[-1] - u.x[0,1] + u.x[0,-1])/(2.*Delta);
    psi[] = 0.;
  }
  boundary ({omega,psi});

  poisson (psi, omega);

  scalar unyielded[];
  foreach()
    unyielded[] = (d.x.x[] == 0. && d.x.x[1] == 0. &&
		   d.y.y[] == 0. && d.y.y[0,1] == 0.);

  char name[80];
  sprintf (name, "fields-%g", tau_y);
  FILE * fp = fopen (name, "w");
  output_field ({psi, unyielded}, fp);
  fclose (fp);

We also output a vertical cross-section of the velocity components.

  sprintf (name, "xprof-%g", tau_y);
  fp = fopen (name, "w");
  for (double y = 0; y <= 1; y += 0.01)
    fprintf (fp, "%g %g %g\n", y,
	     interpolate (u.x, 0.5, y), interpolate (u.y, 0.5, y));
  fclose (fp);

And finally we output the maximum of the streamfunction on standard error.

  fprintf (stderr, "%g %d %g %g\n", t, i, tau_y, statsf(psi).max);
}

Results

set xlabel 'tau_0/sqrt(2) (Vola) or tau_0'
set ylabel 'max(psi)'
set logscale x
set xrange [0.1:1000]
set grid
plot '../vola2003-2a' u ($1/sqrt(2)):2 t 'Vola et al, 2003, fig 2a', \
     'log' u 3:4  w p t 'Basilisk'

set xlabel 'tau_0'
set ylabel ''
plot 'log' u 3:2  w p t ''

reset
set contour base
set cntrparam levels 20
unset surface
set table 'cont.dat'
splot 'fields-0' u 1:2:3 w l
unset table

set size ratio -1
unset key
unset xtics
unset ytics
set palette gray
unset colorbox
set xrange [0:1]
set yrange [0:1]
set cbrange [-1:1]
plot 'fields-0' u 1:2:(1.-$4) w image, 'cont.dat' w l lt -1

Streamlines and unyielded areas: \tau_y=0 (script)

set table 'cont.dat'
splot 'fields-0.707107' u 1:2:3 w l
unset table
plot 'fields-0.707107' u 1:2:(1.-$4) w image, 'cont.dat' w l lt -1

Streamlines and unyielded areas: $_y=1/(script) {2}

set table 'cont.dat'
splot 'fields-7.07107' u 1:2:3 w l
unset table
plot 'fields-7.07107' u 1:2:(1.-$4) w image, 'cont.dat' w l lt -1

Streamlines and unyielded areas: $_y=10/(script) {2}

set table 'cont.dat'
splot 'fields-70.7107' u 1:2:3 w l
unset table
plot 'fields-70.7107' u 1:2:(1.-$4) w image, 'cont.dat' w l lt -1

Streamlines and unyielded areas: $_y=100/(script) {2}

set table 'cont.dat'
splot 'fields-707.107' u 1:2:3 w l
unset table
plot 'fields-707.107' u 1:2:(1.-$4) w image, 'cont.dat' w l lt -1

Streamlines and unyielded areas: $_y=1000/(script) {2}

See Vola et al, 2003, Fig. 3.

reset
set xlabel 'y'
set ylabel 'u.x'
set grid
set key top left
set xrange [0:1]
plot 'xprof-0' w l t 'tau_0=0', \
     'xprof-0.707107' w l t 'tau_0=1/sqrt(2)', \
     'xprof-7.07107' w l t 'tau_0=10/sqrt(2)', \
     'xprof-70.7107' w l t 'tau_0=100/sqrt(2)', \
     'xprof-707.107' w l t 'tau_0=1000/sqrt(2)' lt -1, \
     '../vola2003-3' every :::0::0 lt 1 t 'Vola et al, 2003', \
     '../vola2003-3' every :::1::1 t '', \
     '../vola2003-3' every :::2::2 t '', \
     '../vola2003-3' every :::3::3 t '', \
     '../vola2003-3' every :::4::4 t ''

Section of horizontal velocity along the vertical mid-plane. (script)