iforce.h

Go to the documentation of this file.

/** @file iforce.h
 */
/**
# Interfacial forces
 
We assume that the interfacial acceleration can be expressed as
\f[
\phi\mathbf{n}\delta_s/\rho
\f]
with \f$\mathbf{n}\f$ the interface normal, \f$\delta_s\f$ the interface Dirac
function, \f$\rho\f$ the density and \f$\phi\f$ a generic scalar field. Using
a CSF/Peskin-like approximation, this can be expressed as
\f[
\phi\nabla f/\rho
\f]
with \f$f\f$ the volume fraction field describing the interface.
 
The interfacial force potential \f$\phi\f$ is associated to each VOF
tracer. This is done easily by adding the following [field
attributes](/Basilisk C#field-attributes). */
 
attribute {
  scalar phi;
}
 
/**
Interfacial forces are a source term in the right-hand-side of the
evolution equation for the velocity of the [centered Navier--Stokes
solver](navier-stokes/centered.h) i.e. it is an acceleration. If
necessary, we allocate a new vector field to store it. */
 
/** @brief Event: defaults (i = 0) */
void event_defaults (void) {
  if (is_constant(a.x)) {
    a = {0} /* new vector */;
    for (int _i = 0; _i < _N; _i++) /* foreach_face */ {
      a.x[] = 0.;
      /* dim: a.x[] == Delta/sq(DT */);
    }
  }
}
 
/**
The calculation of the acceleration is done by this event, overloaded
from [its definition](navier-stokes/centered.h#acceleration-term) in
the centered Navier--Stokes solver. */
 
/** @brief Event: acceleration (i++) */
void event_acceleration (void) 
{
  
  /**
  We check for all VOF interfaces for which \f$\phi\f$ is allocated. The
  corresponding volume fraction fields will be stored in *list*. */
 
  scalar * list = NULL;
  for (int _f = 0; _f < 1; _f++) /* scalar in interfaces */
    if (f.phi.i) {
      list = list_add (list, f);
 
      /**
      To avoid undeterminations due to round-off errors, we remove
      values of the volume fraction larger than one or smaller than
      zero. */
 
      for (int _i = 0; _i < _N; _i++) /* foreach */
  f[] = clamp (f[], 0., 1.);
    }
 
  /**
  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 (int _f = 0; _f < 1; _f++) /* scalar in list */
    set_prolongation (f, p.prolongation);
#endif
 
  /**
  Finally, for each interface for which \f$\phi\f$ is allocated, we
  compute the interfacial force acceleration
  \f[
  \phi\mathbf{n}\delta_s/\rho \approx \alpha\phi\nabla f
  \f] 
  */
 
  vector ia = a;
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    for (int _f = 0; _f < 1; _f++) /* scalar in list */
      if (f[] != f[-1] && fm.x[] > 0.) {
 
  /**
  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 phi = f.phi;
  double phif =
    (phi[] < nodata && phi[-1] < nodata) ?
    (phi[] + phi[-1])/2. :
    phi[] < nodata ? phi[] :
    phi[-1] < nodata ? phi[-1] :
    0.;
 
  ia.x[] += alpha.x[]/(fm.x[] + SEPS)*phif*(f[] - f[-1])/Delta;
      }
 
  /**
  On trees, we need to restore the prolongation values for the
  volume fraction field. */
  
#if TREE
  for (int _f = 0; _f < 1; _f++) /* scalar in list */
    set_prolongation (f, fraction_refine);
#endif
  
  /**
  Finally we free the potential fields and the list of volume
  fractions. */
 
  for (int _f = 0; _f < 1; _f++) /* scalar in list */ {
    scalar phi = f.phi;
    delete ({phi});
    f.phi.i = 0;
  }
  free (list);
}
 
/**
## References
 
See Section 3, pages 8-9 of:
 
~~~bib
@hal{popinet2018, hal-01528255}
~~~
*/