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).

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

    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 ((!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;
        scalar * tracers = c.tracers;
        for (scalar t in tracers) {
          t.restriction = restriction_volume_average;
          t.refine = t.prolongation = vof_concentration_refine;
          t.c = c;
        }
      }
    }
    #endif // TREE

    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);
        }
        boundary (gfl);
      }

    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 (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

    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);

    On tree grids, we need to make sure that the fluxes match at fine/coarse cell boundaries i.e. we need to restrict the fluxes from fine cells to coarse cells. This is what is usually done, for all dimensions, by the boundary_flux() function. Here, we only need to do it for a single dimension (x).

    #if TREE
      scalar * fluxl = list_concat (NULL, tfluxl);
      fluxl = list_append (fluxl, flux);
      for (int l = depth() - 1; l >= 0; l--)
        foreach_halo (prolongation, l) {
    #if dimension == 1
          if (is_refined (neighbor(-1)))
    	for (scalar fl in fluxl)
    	  fl[] = fine(fl);
          if (is_refined (neighbor(1)))
    	for (scalar fl in fluxl)
    	  fl[1] = fine(fl,2);
    #elif dimension == 2
          if (is_refined (neighbor(-1)))
    	for (scalar fl in fluxl)
    	  fl[] = (fine(fl,0,0) + fine(fl,0,1))/2.;
          if (is_refined (neighbor(1)))
    	for (scalar fl in fluxl)
    	  fl[1] = (fine(fl,2,0) + fine(fl,2,1))/2.;
    #else // dimension == 3
          if (is_refined (neighbor(-1)))
    	for (scalar fl in fluxl)
    	  fl[] = (fine(fl,0,0,0) + fine(fl,0,1,0) +
    		  fine(fl,0,0,1) + fine(fl,0,1,1))/4.;
          if (is_refined (neighbor(1)))
    	for (scalar fl in fluxl)
    	  fl[1] = (fine(fl,2,0,0) + fine(fl,2,1,0) +
    		   fine(fl,2,0,1) + fine(fl,2,1,1))/4.;
    #endif
        }
      free (fluxl);
    #endif

    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 simply \displaystyle \partial_tt_j = -\nabla_x\cdot(\mathbf{u}_f t_j)

      foreach() {
        c[] += dt*(flux[] - flux[1] + cc[]*(uf.x[1] - uf.x[]))/(cm[]*Delta + SEPS);
        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 + SEPS);
      }
      boundary ({c});
      boundary (tracers);
    
      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      
          t.restriction = restriction_volume_average;
          t.refine = t.prolongation = vof_concentration_refine;
          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[]/(1. - c[]) : 0.;
    	else
    	  tc[] = c[] > 0.5 ? 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

    [lopez2015]

    J. M. Lopez-Herrera, A. M. Ganan-Calvo, S. Popinet, and M. A. Herrada. A vof numerical study on the electrokinetic effects in the breakup of electrified jets. International Journal of Multiphase Flows, 71:14–22, 2015. [ DOI | .pdf ]

    Usage

    Tests