# Surface tension

Surface tension can be expressed as the interfacial force density $\varphi \nabla f$ with $f$ the volume fraction describing the interface and the potential $\varphi =\sigma \kappa$ with $\sigma$ the (constant) surface tension coefficient and $\kappa$ 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=\sqrt{\frac{{\rho }_{m}{\Delta }_{min}^{3}}{\pi \sigma }}$ with ${\rho }_{m}=\left({\rho }_{1}+{\rho }_{2}\right)/2.$ and ${\rho }_{1}$, ${\rho }_{2}$ the densities on either side of the interface.

event stability (i++) {

We first compute the minimum and maximum values of $\alpha /{f}_{m}=1/\rho$, as well as ${\Delta }_{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.;

The maximum timestep is set using the sum of surface tension coefficients.

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

## Definition of the potential

We overload the acceleration event to define the potential $\varphi =\sigma \kappa$.

event acceleration (i++)
{

We check for all VOF interfaces for which $\sigma$ is non-zero.

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

If $\varphi$ is already allocated, we add $\sigma \kappa$, otherwise we allocate a new field and set it to $\sigma \kappa$.

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