Surface tension

Surface tension can be expressed as the interfacial force density ϕf with f the volume fraction describing the interface and the potential ϕ=σκ with σ the (constant) surface tension coefficient and κ the interface mean curvature.

#include "iforce.h"
#include "curvature.h"

The surface tension coefficient is associated to each VOF tracer.

attribute {
  double σ;

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;

Definition of the potential

We overload the acceleration event to define the potential ϕ=σκ.

event acceleration (i++)

We check for all VOF interfaces for which σ is non-zero.

  for (scalar f in interfaces)
    if (f.σ) {

If ϕ is already allocated, we add σκ, otherwise we allocate a new field and set it to σκ.

      scalar φ = f.φ;
      if (φ.i)
	curvature (f, φ, f.σ, add = true);
      else {
	φ = new scalar;
	curvature (f, φ, f.σ, add = false);
	f.φ = φ;