sandbox/lopez/src/vof.h

Volume-Of-Fluid advection

We want to approximate the solution of the advection equations \displaystyle \partial_tc_i + \mathbf{u}_f\cdot\nabla c_i = 0 where c_i are volume fraction fields describing sharp interfaces.

This can be done using a conservative, non-diffusive geometric VOF scheme.

We also add the option to transport diffusive tracers (aka “VOF concentrations”) confined to one side of the interface i.e. solve the equations \displaystyle \partial_tt_{i,j} + \nabla\cdot(\mathbf{u}_ft_{i,j}) = 0 with t_{i,j} = c_if_j (or t_{i,j} = (1 - c_i)f_j) and f_j is a volumetric tracer concentration (see Lopez-Herrera et al, 2015).

The list of tracers associated with the volume fraction is stored in the tracers attribute. For each tracer, the “side” of the interface (i.e. either c or 1 - c) is controlled by the inverse attribute). The boolean attribute comp switch off the divergence free correction in the tracer equation if it is set to “true”.

attribute {
  scalar * tracers, c;
  bool inverse, comp;
}

We will need basic functions for volume fraction computations.

#include "fractions.h"

The list of volume fraction fields interfaces, will be provided by the user.

The face velocity field uf will be defined by a solver as well as the timestep.

extern scalar * interfaces;
extern face vector uf;
extern double dt;

The gradient of a VOF-concentration t is computed using a standard three-point scheme if we are far enough from the interface (as controlled by cmin), otherwise a two-point scheme biased away from the interface is used.

foreach_dimension()
static double vof_concentration_gradient_x (Point point, scalar c, scalar t)
{
  static const double cmin = 0.5;
  double cl = c[-1], cc = c[], cr = c[1];
  if (t.inverse)
    cl = 1. - cl, cc = 1. - cc, cr = 1. - cr;
  if (cc >= cmin && t.gradient != zero) {
    if (cr >= cmin) {
      if (cl >= cmin) {
	if (t.gradient)
	  return t.gradient (t[-1]/cl, t[]/cc, t[1]/cr)/Delta;
	else
	  return (t[1]/cr - t[-1]/cl)/(2.*Delta);
      }
      else
	return (t[1]/cr - t[]/cc)/Delta;
    }
    else if (cl >= cmin)
      return (t[]/cc - t[-1]/cl)/Delta;
  }
  return 0.;
}

On trees, VOF concentrations need to be refined properly i.e. using volume-fraction-weighted linear interpolation of the concentration.

#if TREE
static void vof_concentration_refine (Point point, scalar s)
{
  scalar f = s.c;
  if (cm[] == 0. || (!s.inverse && f[] <= 0.) || (s.inverse && f[] >= 1.))
    foreach_child()
      s[] = 0.;
  else {
    coord g;
    foreach_dimension()
      g.x = Delta*vof_concentration_gradient_x (point, f, s);
    double sc = s.inverse ? s[]/(1. - f[]) : s[]/f[], cmc = 4.*cm[];
    foreach_child() {
      s[] = sc;
      foreach_dimension()
	s[] += child.x*g.x*cm[-child.x]/cmc;
      s[] *= s.inverse ? 1. - f[] : f[];
    }
  }
}

On trees, we need to setup the appropriate prolongation and refinement functions for the volume fraction fields.

event defaults (i = 0)
{
  for (scalar c in interfaces) {
    c.refine = c.prolongation = fraction_refine;
    c.dirty = true;
    scalar * tracers = c.tracers;
    for (scalar t in tracers) {
      t.restriction = restriction_volume_average;
      t.refine = t.prolongation = vof_concentration_refine;
      t.dirty = true;
      t.c = c;
    }
  }
}
#endif // TREE

Boundary conditions for VOF-advected tracers usually depend on boundary conditions for the VOF field.

event defaults (i = 0)
{
  for (scalar c in interfaces) {
    scalar * tracers = c.tracers;
    for (scalar t in tracers)
      t.depends = list_add (t.depends, c);
  }
}

We need to make sure that the CFL is smaller than 0.5 to ensure stability of the VOF scheme.

event stability (i++) {
  if (CFL > 0.5)
    CFL = 0.5;
}

One-dimensional advection

The simplest way to implement a multi-dimensional VOF advection scheme is to use dimension-splitting i.e. advect the field along each dimension successively using a one-dimensional scheme.

We implement the one-dimensional scheme along the x-dimension and use the foreach_dimension() operator to automatically derive the corresponding functions along the other dimensions.

foreach_dimension()
static void sweep_x (scalar c, scalar cc, scalar * tcl)
{
  vector n[];
  scalar alpha[], flux[];
  double cfl = 0.;

If we are also transporting tracers associated with c, we need to compute their gradient i.e. \partial_xf_j = \partial_x(t_j/c) or \partial_xf_j = \partial_x(t_j/(1 - c)) (for higher-order upwinding) and we need to store the computed fluxes. We first allocate the corresponding lists.

  scalar * tracers = c.tracers, * gfl = NULL, * tfluxl = NULL;
  if (tracers) {
    for (scalar t in tracers) {
      scalar gf = new scalar, flux = new scalar;
      gfl = list_append (gfl, gf);
      tfluxl = list_append (tfluxl, flux);
    }

The gradient is computed using the “interface-biased” scheme above.

    foreach() {
      scalar t, gf;
      for (t,gf in tracers,gfl)
	gf[] = vof_concentration_gradient_x (point, c, t);
    }
  }

We reconstruct the interface normal \mathbf{n} and the intercept \alpha for each cell. Then we go through each (vertical) face of the grid.

  reconstruction (c, n, alpha);
  foreach_face(x, reduction (max:cfl)) {

To compute the volume fraction flux, we check the sign of the velocity component normal to the face and compute the index i of the corresponding upwind cell (either 0 or -1).

    double un = uf.x[]*dt/(Delta*fm.x[] + SEPS), s = sign(un);
    int i = -(s + 1.)/2.;

We also check that we are not violating the CFL condition.

#if EMBED
    if (cs[] >= 1.)
#endif
    if (un*fm.x[]*s/(cm[] + SEPS) > cfl)
      cfl = un*fm.x[]*s/(cm[] + SEPS);

If we assume that un is negative i.e. s is -1 and i is 0, the volume fraction flux through the face of the cell is given by the dark area in the figure below. The corresponding volume fraction can be computed using the rectangle_fraction() function.

Volume fraction flux

When the upwind cell is entirely full or empty we can avoid this computation.

    double cf = (c[i] <= 0. || c[i] >= 1.) ? c[i] :
      rectangle_fraction ((coord){-s*n.x[i], n.y[i], n.z[i]}, alpha[i],
			  (coord){-0.5, -0.5, -0.5},
			  (coord){s*un - 0.5, 0.5, 0.5});

Once we have the upwind volume fraction cf, the volume fraction flux through the face is simply:

    flux[] = cf*uf.x[];

If we are transporting tracers, we compute their flux using the upwind volume fraction cf and a tracer value upwinded using the Bell–Collela–Glaz scheme and the gradient computed above.

    scalar t, gf, tflux;
    for (t,gf,tflux in tracers,gfl,tfluxl) {
      double cf1 = cf, ci = c[i];
      if (t.inverse)
	cf1 = 1. - cf1, ci = 1. - ci;
      if (ci > 1e-10) {
	double ff = t[i]/ci + s*min(1., 1. - s*un)*gf[i]*Delta/2.;
	tflux[] = ff*cf1*uf.x[];
      }
      else
	tflux[] = 0.;
    }
  }
  delete (gfl); free (gfl);

We warn the user if the CFL condition has been violated.

  if (cfl > 0.5 + 1e-6)
    fprintf (ferr, 
	     "WARNING: CFL must be <= 0.5 for VOF (cfl - 0.5 = %g)\n", 
	     cfl - 0.5), fflush (ferr);

Once we have computed the fluxes on all faces, we can update the volume fraction field according to the one-dimensional advection equation \displaystyle \partial_tc = -\nabla_x\cdot(\mathbf{u}_f c) + c\nabla_x\cdot\mathbf{u}_f The first term is computed using the fluxes. The second term – which is non-zero for the one-dimensional velocity field – is approximated using a centered volume fraction field cc which will be defined below.

For tracers, the one-dimensional update is, following the same scheme that above, simply \displaystyle \partial_tt_j = -\nabla_x\cdot(\mathbf{u}_f t_j) + f_j \nabla_x\cdot\mathbf{u}_f The second term, a correction for non perfect divergence free velocity field, can be removed setting to true the attribute comp.

#if !EMBED
  foreach() {
    c[] += dt*(flux[] - flux[1] + cc[]*(uf.x[1] - uf.x[]))/(cm[]*Delta);
    scalar t, tc, tflux;
    for (t, tc, tflux in tracers, tcl, tfluxl)
      t[] += dt*(tflux[] - tflux[1] + tc[]*(uf.x[1] - uf.x[]))/(cm[]*Delta);
  }
#else // EMBED

When dealing with embedded boundaries, we simply ignore the fraction occupied by the solid. This is a simple approximation which has the advantage of ensuring boundedness of the volume fraction and conservation of the total tracer mass (if it is computed also ignoring the volume occupied by the solid in partial cells).

  foreach()
    if (cs[] > 0.) {
      c[] += dt*(flux[] - flux[1] + cc[]*(uf.x[1] - uf.x[]))/Delta;
      scalar t, tc, tflux;
      for (t, tc, tflux in tracers, tcl, tfluxl)
	t[] += dt*(tflux[] - tflux[1] + tc[]*(uf.x[1] - uf.x[]))/Delta;
    }
#endif // EMBED

  delete (tfluxl); free (tfluxl);
}

Multi-dimensional advection

The multi-dimensional advection is performed by the event below.

void vof_advection (scalar * interfaces, int i)
{
  for (scalar c in interfaces) {

We first define the volume fraction field used to compute the divergent term in the one-dimensional advection equation above. We follow Weymouth & Yue, 2010 and use a step function which guarantees exact mass conservation for the multi-dimensional advection scheme (provided the advection velocity field is exactly non-divergent).

    scalar cc[], * tcl = NULL, * tracers = c.tracers;    
    for (scalar t in tracers) {
      scalar tc = new scalar;
      tcl = list_append (tcl, tc);
#if TREE
      if (t.refine != vof_concentration_refine) {
	t.refine = t.prolongation = vof_concentration_refine;
	t.restriction = restriction_volume_average;
	t.dirty = true;
	t.c = c;
      }
#endif // TREE
    }
    foreach() {
      cc[] = (c[] > 0.5);
      scalar t, tc;
      for (t, tc in tracers, tcl) {
	if (t.inverse)
	  tc[] = (c[] < 0.5 && !t.comp ? t[]/(1.-c[]) : 0.);
	else
	  tc[] = (c[] > 0.5 && !t.comp ? t[]/c[] : 0.);
      }
    }

We then apply the one-dimensional advection scheme along each dimension. To try to minimise phase errors, we alternate dimensions according to the parity of the iteration index i.

    void (* sweep[dimension]) (scalar, scalar, scalar *);
    int d = 0;
    foreach_dimension()
      sweep[d++] = sweep_x;
    for (d = 0; d < dimension; d++)
      sweep[(i + d) % dimension] (c, cc, tcl);
    delete (tcl), free (tcl);
  }
}

event vof (i++)
  vof_advection (interfaces, i);

References