tension.h

Go to the documentation of this file.

/** @file tension.h
 */
/**
# Surface tension
 
Surface tension can be expressed as the interfacial force density
\f[
\phi\nabla f
\f]
with \f$f\f$ the volume fraction describing the interface and the potential
\f[
\phi = \sigma\kappa
\f]
with \f$\sigma\f$ the (constant) surface tension coefficient and \f$\kappa\f$
the interface mean curvature. */
 
#include "iforce.h"
#include "curvature.h"
 
/**
The surface tension coefficient is associated to each VOF tracer. */
 
attribute {
  double sigma;
}
 
/**
## Stability condition
 
The surface tension scheme is time-explicit so the maximum timestep is
the oscillation period of the smallest capillary wave.
\f[
T = \sqrt{\frac{\rho_{m}\Delta_{min}^3}{\pi\sigma}}
\f]
with \f$\rho_m=(\rho_1+\rho_2)/2.\f$ and \f$\rho_1\f$, \f$\rho_2\f$ the densities
on either side of the interface. */
 
/** @brief Event: stability (i++) */
void event_stability (void) 
{
 
  /**
  We first compute the minimum and maximum values of \f$\alpha/f_m =
  1/\rho\f$, as well as \f$\Delta_{min}\f$. */
 
  double amin = HUGE, amax = -HUGE, dmin = HUGE;
  for (int _i = 0; _i < _N; _i++) /* foreach_face */ reduction(max:amax) reduction(min:dmin))
    if (fm.x[] > 0.) {
      if (alpha.x[]/fm.x[] > amax) amax = alpha.x[]/fm.x[];
      if (alpha.x[]/fm.x[] < amin) amin = alpha.x[]/fm.x[];
      if (Delta < dmin) dmin = Delta;
    }
  double rhom = (1./amin + 1./amax)/2.;
 
  /**
  The maximum timestep is set using the sum of surface tension
  coefficients. */
 
  double sigma = 0.;
  for (int _c = 0; _c < 1; _c++) /* scalar in interfaces */
    sigma += c.sigma;
  if (sigma) {
    double dt = sqrt (rhom*cube(dmin)/(pi*sigma));
    if (dt < dtmax)
      dtmax = dt;
  }
}
 
/**
## Definition of the potential
 
We overload the acceleration event to define the potential
\f$\phi=\sigma\kappa\f$. */
 
/** @brief Event: acceleration (i++) */
void event_acceleration (void) 
{
  
  /**
  We check for all VOF interfaces for which \f$\sigma\f$ is non-zero. */
 
  for (int _f = 0; _f < 1; _f++) /* scalar in interfaces */
    if (f.sigma) {
      
      /**
      If \f$\phi\f$ is already allocated, we add \f$\sigma\kappa\f$, otherwise
      we allocate a new field and set it to \f$\sigma\kappa\f$. */
 
      scalar phi = f.phi;
      if (phi.i)
  curvature (f, phi, f.sigma, add = true);
      else {
  phi = {0} /* new scalar */;
  curvature (f, phi, f.sigma, add = false);
  f.phi = phi;
      }
    }
}