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.

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

Vortex intensity (script)

$Streamlines for \tau_y=0 (script)$

Streamlines for $\tau_y=0$ (script)

Streamlines for $_y=1/(script) {2}

Streamlines for $_y=10/(script) {2}