Curvature and surface tension force

Interfacial curvature is often needed to describe the surface tension force. Although not directly relevant to the volume-of-fluid method, here I briefly explain how it is computed.

Mathematically speaking, the surface tension force (only normal component is considered, and the tangential components are ignored) is given by

\[\frac{1}{We} \kappa \delta n_i,\]

which is approximated as

\[\frac{1}{We} \kappa \der{\phi}{x_i},\]

where \(\phi\) is the volume-of-fluid.

The pre-factor (inverse of \(We\)) is computed in the initialisation process:

src/interface/init.c

double We = 0.;
if(0 != config.get_double("We", &We)) return 1;
interface->tension = 1. / We;

\(\kappa\) is the local mean curvature, whose definition is

\[\kappa \equiv -\der{n_i}{x_i}.\]

In this project, this is directly evaluated using the normal vectors computed at the cell corners, i.e.:

\[\begin{split}\vat{\kappa}{\pic, \pjc} & = -\vat{\der{n_x}{x}}{\pic, \pjc} -\vat{\der{n_y}{y}}{\pic, \pjc} \\ & \approx -\frac{\vat{n_x}{\pip, \pjc} - \vat{n_x}{\pim, \pjc}}{\vat{\Delta x}{\pic}} -\frac{\vat{n_y}{\pic, \pjp} - \vat{n_y}{\pic, \pjm}}{ \Delta y }.\end{split}\]

The right-hand-side terms have already been computed in the previous step.

src/interface/curvature_tensor.c

const double dnxdx = 1. / dx * (
    - DVOF(i  , j  )[0] + DVOF(i+1, j  )[0]
    - DVOF(i  , j+1)[0] + DVOF(i+1, j+1)[0]
);
const double dnydy = 1. / dy * (
    - DVOF(i  , j  )[1] - DVOF(i+1, j  )[1]
    + DVOF(i  , j+1)[1] + DVOF(i+1, j+1)[1]
);
CURV(i, j) = 0.5 * (
  - dnxdx
  - dnydy
);

src/interface/curvature_tensor.c

const double dnxdx = 1. / dx * (
    - DVOF(i  , j  , k  )[0] + DVOF(i+1, j  , k  )[0]
    - DVOF(i  , j+1, k  )[0] + DVOF(i+1, j+1, k  )[0]
    - DVOF(i  , j  , k+1)[0] + DVOF(i+1, j  , k+1)[0]
    - DVOF(i  , j+1, k+1)[0] + DVOF(i+1, j+1, k+1)[0]
);
const double dnydy = 1. / dy * (
    - DVOF(i  , j  , k  )[1] - DVOF(i+1, j  , k  )[1]
    + DVOF(i  , j+1, k  )[1] + DVOF(i+1, j+1, k  )[1]
    - DVOF(i  , j  , k+1)[1] - DVOF(i+1, j  , k+1)[1]
    + DVOF(i  , j+1, k+1)[1] + DVOF(i+1, j+1, k+1)[1]
);
const double dnzdz = 1. / dz * (
    - DVOF(i  , j  , k  )[2] - DVOF(i+1, j  , k  )[2]
    - DVOF(i  , j+1, k  )[2] - DVOF(i+1, j+1, k  )[2]
    + DVOF(i  , j  , k+1)[2] + DVOF(i+1, j  , k+1)[2]
    + DVOF(i  , j+1, k+1)[2] + DVOF(i+1, j+1, k+1)[2]
);
CURV(i, j, k) = 0.25 * (
  - dnxdx
  - dnydy
  - dnzdz
);

As a result, the surface tension force in the \(x\) direction is described as

\[\vat{f_x^{st}}{\pip, \pjc} \approx \frac{1}{We} \frac{ \vat{\kappa}{\pipp, \pjc} + \vat{\kappa}{\pic, \pjc} }{2} \frac{ \vat{\phi}{\pipp, \pjc} - \vat{\phi}{\pic, \pjc} }{\vat{\Delta x}{\pip}}.\]

Same applies to the other directions.

\(x\) direction:

src/interface/force.c

const double dx = DXC(i  );
const double grad = 1. / dx * (
    - VOF(i-1, j  )
    + VOF(i  , j  )
);
const double kappa = 0.5 * (
    + CURV(i-1, j  )
    + CURV(i  , j  )
);
IFRCX(i, j) = tension * grad * kappa;

src/interface/force.c

const double dx = DXC(i  );
const double grad = 1. / dx * (
    - VOF(i-1, j  , k  )
    + VOF(i  , j  , k  )
);
const double kappa = 0.5 * (
    + CURV(i-1, j  , k  )
    + CURV(i  , j  , k  )
);
IFRCX(i, j, k) = tension * grad * kappa;

\(y\) direction:

src/interface/force.c

const double grad = 1. / dy * (
    - VOF(i  , j-1)
    + VOF(i  , j  )
);
const double kappa = 0.5 * (
    + CURV(i  , j-1)
    + CURV(i  , j  )
);
IFRCY(i, j) = tension * grad * kappa;

src/interface/force.c

const double grad = 1. / dy * (
    - VOF(i  , j-1, k  )
    + VOF(i  , j  , k  )
);
const double kappa = 0.5 * (
    + CURV(i  , j-1, k  )
    + CURV(i  , j  , k  )
);
IFRCY(i, j, k) = tension * grad * kappa;

\(z\) direction:

src/interface/force.c

const double grad = 1. / dz * (
    - VOF(i  , j  , k-1)
    + VOF(i  , j  , k  )
);
const double kappa = 0.5 * (
    + CURV(i  , j  , k-1)
    + CURV(i  , j  , k  )
);
IFRCZ(i, j, k) = tension * grad * kappa;

Note

Here I adopt a very simple approach, but more accurate evaluation of the surface curvature is a challenging problem and not conclusive at all. See Popinet, Annu. Rev. Fluid Mech. (50), 2018 for instance.