Surface tension

We will need to compute the curvature of the interface, using its Volume-Of-Fluid description.

#include "curvature.h"

The surface tension σ and interface curvature κ will be associated to each VOF tracer. This is done easily by adding the following field attributes.

attribute {
  double σ;
  scalar κ;

event defaults (i = 0) {

Surface tension is a source term in the right-hand-side of the evolution equation for the velocity of the centered Navier–Stokes solver i.e. it is an acceleration. If necessary, we allocate a new vector field to store it.

  if (is_constant(a.x)) {
    a = new face vector;
      a.x[] = 0.;
    boundary ((scalar *){a});

Each interface for which σ is not zero needs a new field to store the curvature.

  for (scalar c in interfaces)
    if (c.σ && !c.κ.i) {
      scalar κ = new_scalar ("kappa");
	κ[] = 0.;
      boundary ({κ});
      c.κ = κ;

Stability condition

The surface tension scheme is time-explicit so the maximum timestep is the oscillation period of the smallest capillary wave. T=ρmΔmin3πσ with ρm=(ρ1+ρ2)/2. and ρ1, ρ2 the densities on either side of the interface.

event stability (i++) {

We first compute the minimum and maximum values of α/fm=1/ρ, as well as Δmin.

  double amin = HUGE, amax = -HUGE, dmin = HUGE;
  foreach_face (reduction(min:amin) reduction(max:amax) reduction(min:dmin)) {
    if (α.x[]/fm.x[] > amax) amax = α.x[]/fm.x[];
    if (α.x[]/fm.x[] < amin) amin = α.x[]/fm.x[];
    if (Δ < dmin) dmin = Δ;
  double rhom = (1./amin + 1./amax)/2.;

We then consider each VOF interface with an associated value of σ different from zero and set the maximum timestep.

  for (scalar c in interfaces)
    if (c.σ) {
      double dt = sqrt (rhom*cube(dmin)/(π*c.σ));
      if (dt < dtmax)
	dtmax = dt;

Surface tension term

The calculation of the acceleration is done by this event, overloaded from its definition in the centered Navier–Stokes solver.

event acceleration (i++)

We check for all VOF interfaces for which σ is non-zero. The corresponding volume fraction fields will be stored in list.

  scalar * list = NULL;
  for (scalar c in interfaces)
    if (c.σ) {
      list = list_add (list, c);

To avoid undeterminations due to round-off errors, we remove values of the volume fraction larger than one or smaller than zero.

	c[] = clamp (c[], 0, 1);
      boundary ({c});

We update the values of κ using the height-function curvature calculation.

      assert (c.κ.i);
      curvature (c, c.κ);

On trees we need to make sure that the volume fraction gradient is computed exactly like the pressure gradient. This is necessary to ensure well-balancing of the pressure gradient and surface tension term. To do so, we apply the same prolongation to the volume fraction field as applied to the pressure field.

#if TREE
  for (scalar c in list)
    c.prolongation = p.prolongation;
  boundary (list);

Finally, for each interface for which σ is non-zero, we compute the surface tension acceleration σκnδs/ρασκc

  face vector st = a;
    for (scalar c in list)
      if (c[] != c[-1]) {
	scalar κ = c.κ;

We need to compute the curvature kf on the face, using its values at the center of the cell. If both curvatures are defined, we take the average, otherwise we take a single value. If all fails we set the curvature to zero: this should happen only because of very pathological cases e.g. weird boundary conditions for the volume fraction.

	double kf = 
	  (κ[] < nodata && κ[-1] < nodata) ?
	     (κ[] + κ[-1])/2. :
	  κ[] < nodata ? κ[] :
	  κ[-1] < nodata ? κ[-1] :
	st.x[] += α.x[]/fm.x[]*c.σ*kf*(c[] - c[-1])/Δ;

On trees, we need to restore the prolongation values for the volume fraction field.

#if TREE
  for (scalar c in list)
    c.prolongation = fraction_refine;
  boundary (list);

Finally we free the list of interfacial volume fractions.

  free (list);