lid-driven-cavity.c

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

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

#include "grid/multigrid.h"
#include "../mycentered2.h"
#include "../myviscosity-viscoplastic.h"

Reference solution

#define l    (1.) // Reference length
#define nu   (1.) // Viscosity
#define uref (1.) // Reference velocity
#define tref (sq (l)/(nu)) // Reference time

Set up

We need a field for the viscosity and the yield stress. It also seems that we need a field for the density, but this may be a bug.

face vector muv[], Tv[];
scalar rhov[];
double T1 = 0.;
int iT = 0;

We also need a reference velocity field un.

scalar un[];

Finally, we define the level of refinement.

#define lvl (7)

Code

int main()
{

The domain is 1\times 1.

  L0 = (l);
  size (L0);
  origin (0., 0.);

We set the maximum timestep.

  DT = 1.e-2*(tref);

We set the tolerance of the Poisson solver.

  stokes       = true;
  TOLERANCE    = 1.e-5;
  TOLERANCE_MU = 1.e-5*(uref);
  NITERMAX     = 500;

We initialize the grid. We start with a Newtonian fluid.

  N = 1 << (lvl);
  init_grid (N);

  T1 = 0.;
  
  run ();
}

Boundary conditions

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

We give boundary conditions for the face velocity to “potentially” improve the convergence of the multigrid Poisson solver and to make sure the compatibility condition for the Poisson equation is not violated.

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

Properties

event properties (i++)
{

We set the density and account for the metric.

  foreach()
    rhov[] = cm[];
  boundary ({rhov});

We now set the viscosity and yield stress T and account for the metric.

  foreach_face() {
    muv.x[] = fm.x[];
    Tv.x[]  = fm.x[]*(T1);
  }
  boundary ((scalar *) {muv, Tv});
}

Initial conditions

event init (i = 0)
{

We set the viscosity and density fields in the event properties.

  rho = rhov;
  mu  = muv;

We also set the yield stress and the regularization parameters.

  T = new face vector;
  Tv = T;

We finally initialize the reference velocity field.

  foreach()
    un[] = u.x[];
}

Outputs

event logfile (i++; i <= 1000)
{
  double du = change (u.x, un);
  fprintf (ferr, "%d %g %g %g %g %g %d %d\n",
	   i, t, dt, eps,
	   T1,
	   du,
	   mgp.i, mgu.i);
  fflush (ferr);

  if (i > 0 && du < 1e-5) {

We compute the streamfunction \psi and the unyielded field and output that to a file.

    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 ((scalar *) {omega, psi});

    poisson (psi, omega);
    
    yielded_region ();

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

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

    sprintf (name, "xprof-%g", T1);
    fp = fopen (name, "w");

    double step = (L0)/(N);
    for (double y = step/2.; y <= (L0) - step/2.; y += step)
      fprintf (fp, "%g %g %g\n", y,
	       interpolate (u.x, (L0)/2., y),
	       interpolate (u.y, (L0)/2., y));
    fclose (fp);

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

    fprintf (fout, "%g %d %g %g\n", t, i, T1, statsf(psi).max);
    fflush (fout);

We increase the value of the yield stress and use the previous solution as an initial guess.

    T1 = pow (10., iT)/sqrt (2);

To ensure convergence of the multigrid solver for large Bingham numbers, we reduce the time step and increase the regularization parameter. These parameters were chosen by trial and error.

    DT  = 1.e-2*(tref)/max (1, pow (10., iT));
    eps = min (5.e-2, max (1.e-4, 5.e-4*pow (10., iT))); 
    
    event ("properties");
    boundary (all);

    if (iT == 4)
      return 1; /* stop */
    iT++;
  }
}

Results

See Vola et al, 2003, Fig. 3.

set key top left spacing 1.1
set xlabel 'y'
set ylabel 'u.x'
set grid
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)', \
     '../data/Vola2003/vola2003-3' every :::0::0 lt 1 t 'Vola et al, 2003', \
     '../data/Vola2003/vola2003-3' every :::1::1 t '', \
     '../data/Vola2003/vola2003-3' every :::2::2 t '', \
     '../data/Vola2003/vola2003-3' every :::3::3 t '', \
     '../data/Vola2003/vola2003-3' every :::4::4 t ''

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

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

Vortex intensity (script)

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 for \tau_y=0 (script)$

Streamlines for \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 for $_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 for $_y=10/(script) {2}