tension.h

Go to the documentation of this file.

/** @file 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 \f$p_1\f$ and \f$p_2\f$ from the averaged pressure \f$p\f$,
we interpret the averaged pressure as
\f[p = f p_1 + (1-f) p_2\f]
Using the Laplace equation
\f[p_2 - p_1 = -\sigma \kappa + [[\mu n \cdot \tau \cdot n]]\f]
we can solve the system of the two equations for \f$p_1\f$ and \f$p_2\f$
\f[  p_1  = p + (1 - f) \sigma \kappa - (1-f) [[\mu n \cdot \tau \cdot n]]\f]
\f[  p_2  = p - f \sigma \kappa + f [[\mu n \cdot \tau \cdot n]]\f]
The flux terms depending on pressure entering into the conservative total energy equation are
\f[\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]])\f]
\f[\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]])\f]
 
In this module we only introduce the correction due to surface tension. */
 
/** @brief Event: end_timestep (i++) */
void event_end_timestep (void) 
{
  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. */
    
    vector upf[];
    for (int _fE = 0; _fE < 1; _fE++) /* scalar in {fE1, fE2} */ {
      for (int _i = 0; _i < _N; _i++) /* 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.;
      }
 
      for (int _i = 0; _i < _N; _i++) /* foreach */ {
  double div = 0.;
  for (int _d = 0; _d < dimension; _d++)
    div += upf.x[1] - upf.x[];
  fE[] -= (fE.inverse ? 1. - f[] : f[])*div/Delta*dt/cm[];
      }
    }
  }
}