src/compressible/tension.h
Surface tension effects for the compressible solver
This module incorporates surface tension effects in the compressible solver. Unlike momentum, which is defined for the averaged mixture, the energy is defined separately for both phases. For that reason, we need to reconstruct the pressure of each phase from the averaged pressure obtained from the Helmholtz–Poisson equation.
To reconstruct p_1 and p_2 from the averaged pressure p, we interpret the averaged pressure as \displaystyle p = f p_1 + (1-f) p_2 Using the Laplace equation \displaystyle p_2 - p_1 = -\sigma \kappa + [[\mu n \cdot \tau \cdot n]] we can solve the system of the two equations for p_1 and p_2 \displaystyle p_1 = p + (1 - f) \sigma \kappa - (1-f) [[\mu n \cdot \tau \cdot n]] \displaystyle p_2 = p - f \sigma \kappa + f [[\mu n \cdot \tau \cdot n]] The flux terms depending on pressure entering into the conservative total energy equation are \displaystyle \nabla \cdot (u p_1) = \nabla \cdot (u p) + \nabla \cdot (u (1-f) \sigma \kappa) - \nabla \cdot (u (1-f) [[\mu n \cdot \tau \cdot n]]) \displaystyle \nabla \cdot (u p_2) = \nabla \cdot (u p) - \nabla \cdot (u f \sigma \kappa) + \nabla \cdot (u f [[\mu n \cdot \tau \cdot n]])
In this module we only introduce the correction due to surface tension.
event end_timestep (i++)
{
if (f.sigma > 0.) {
scalar skappa[];
curvature (f, skappa, f.sigma);Here we just introduce the correction due to the surface tension contribution to the pressure jump.
face vector upf[];
for (scalar fE in {fE1, fE2}) {
foreach_face() {
double fr = f[1], fl = f[];
if (!fE.inverse)
fr = 1. - fr, fl = 1. - fl;
if (fr + fl > 0. && fr + fl < 2.) {
if (fr > fl)
upf.x[] = uf.x[]*fl*skappa[];
else
upf.x[] = uf.x[]*fr*skappa[1];
}
else
upf.x[] = 0.;
}
foreach() {
double div = 0.;
foreach_dimension()
div += upf.x[1] - upf.x[];
fE[] -= (fE.inverse ? 1. - f[] : f[])*div/Delta*dt/cm[];
}
}
}
}