# src/tension.h

# Surface tension

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

`#include "curvature.h"`

The surface tension $\sigma $ and interface curvature $\kappa $ 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;
foreach_face()
a.x[] = 0.;
boundary ((scalar *){a});
}
```

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

```
for (scalar c in interfaces)
if (c.σ && !c.κ.i) {
scalar κ = new_scalar ("kappa");
foreach()
κ[] = 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=\sqrt{\frac{{\rho}_{m}{\Delta}_{min}^{3}}{\pi \sigma}}$$ with ${\rho}_{m}=({\rho}_{1}+{\rho}_{2})/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.;
```

We then consider each VOF interface with an associated value of $\sigma $ 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 $\sigma $ 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.

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

We update the values of $\kappa $ 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);
#endif
```

Finally, for each interface for which $\sigma $ is non-zero, we compute the surface tension acceleration $$\sigma \kappa \mathbf{\text{n}}{\delta}_{s}/\rho \approx \alpha \sigma \kappa \nabla c$$

```
face vector st = a;
foreach_face()
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] :
0.;
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);
#endif
```

Finally we free the list of interfacial volume fractions.

```
free (list);
}
```