# Interfacial forces

We assume that the interfacial acceleration can be expressed as $\varphi \mathbf{\text{n}}{\delta }_{s}/\rho$ with $\mathbf{\text{n}}$ the interface normal, ${\delta }_{s}$ the interface Dirac function, $\rho$ the density and $\varphi$ a generic scalar field. Using a CSF/Peskin-like approximation, this can be expressed as $\varphi \nabla f/\rho$ with $f$ the volume fraction field describing the interface.

The interfacial force potential $\varphi$ is associated to each VOF tracer. This is done easily by adding the following field attributes.

``````attribute {
scalar φ;
}``````

Interfacial forces are 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.

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

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 $\varphi$ is allocated. The corresponding volume fraction fields will be stored in list.

``````  scalar * list = NULL;
for (scalar f in interfaces)
if (f.φ.i) {

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

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

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 interfacial force term. To do so, we apply the same prolongation to the volume fraction field as applied to the pressure field.

``````#if TREE
for (scalar f in list)
f.prolongation = p.prolongation;
boundary (list);
#endif``````

Finally, for each interface for which $\varphi$ is allocated, we compute the interfacial force acceleration $\varphi \mathbf{\text{n}}{\delta }_{s}/\rho \approx \alpha \varphi \nabla f$

``````  face vector ia = a;
foreach_face()
for (scalar f in list)
if (f[] != f[-1]) {``````

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

``````	scalar φ = f.φ;
double phif =
(φ[] < nodata && φ[-1] < nodata) ?
(φ[] + φ[-1])/2. :
φ[] < nodata ? φ[] :
φ[-1] < nodata ? φ[-1] :
0.;

ia.x[] += α.x[]/fm.x[]*phif*(f[] - f[-1])/Δ;
}``````

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

``````#if TREE
for (scalar f in list)
f.prolongation = fraction_refine;
boundary (list);
#endif``````

Finally we free the potential fields and the list of volume fractions.

``````  for (scalar f in list) {
scalar φ = f.φ;
delete ({φ});
f.φ.i = 0;
}
free (list);
}``````

## References

See Section 3, pages 8-9 of:

 [popinet2018] S. Popinet. Numerical models of surface tension. Annual Review of Fluid Mechanics, 50:1-28, 2018. [ DOI | http ]