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

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

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 corresponding one-dimensional update is \displaystyle \partial_tt_j = -\nabla_x\cdot(\mathbf{u}_f t_j) + t_j\nabla_x\cdot\mathbf{u}_f The compressive second-term can be removed by defining the NO_1D_COMPRESSION macro.

```
#if !EMBED
foreach() {
c[] += dt*(flux[] - flux[1] + cc[]*(uf.x[1] - uf.x[]))/(cm[]*Delta);
#if NO_1D_COMPRESSION
scalar t, tflux;
for (t, tflux in tracers, tfluxl)
t[] += dt*(tflux[] - tflux[1])/(cm[]*Delta);
#else // !NO_1D_COMPRESSION
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);
#endif // !NO_1D_COMPRESSION
}
#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;
#if NO_1D_COMPRESSION
for (t, tflux in tracers, tfluxl)
t[] += dt*(tflux[] - tflux[1])/Delta;
#else // !NO_1D_COMPRESSION
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 // !NO_1D_COMPRESSION
}
#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) {
#if !NO_1D_COMPRESSION
scalar tc = new scalar;
tcl = list_append (tcl, tc);
#endif // !NO_1D_COMPRESSION
#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);
#if !NO_1D_COMPRESSION
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.;
}
#endif // !NO_1D_COMPRESSION
}
```

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