sandbox/ecipriano/src/icentripetal.h

Interfacial Suspending Force

We re-express the suspending force in centripetal.h as an interfacial force. The suspending force is expressed as a function of the gradient of a potential: \displaystyle \mathbf{f}_m = \rho \alpha \nabla \xi As explained by Popinet 2018, a source term expressed as a gradient of a potential can be discretized using a well-balanced scheme, similarly to the surface tension force. Setting \mathbf{u}=0 in the Navier-Stokes equations: \displaystyle -\nabla p + \rho\nabla\xi = -\nabla p + \nabla(\rho\xi) - \xi\nabla\rho = -\nabla p' - \xi\nabla\rho with \rho which jumps from a constant value \rho_l in liquid phase, to the value \rho_g in gas phase according to: \displaystyle \rho(H) = (\rho_l - \rho_g)H +\rho_l = [\rho]H + \rho_l whose gradient reads: \displaystyle \nabla\rho = [\rho]\nabla H Therefore, the body force \mathbf{f}_m is replaced by an interfacial force, using an approach similar to the one used in tension.h and reduced.h: \displaystyle -\nabla p' - \xi\nabla\rho = -\nabla p' - \phi\nabla H where the term \phi = \xi[\rho] is computed here and added to the surface forces by iforce.h. The potential is computed as: \displaystyle \xi = \dfrac{\xi_0}{|\mathbf{x}_\Gamma - \mathbf{x}_c|} where \xi_0 is an arbitrary constant that regulates the intensity of the force; \mathbf{x}_\Gamma is the center of the interface element; \mathbf{x}_c is the point where the center of the droplet should remain fixed.

#define CENTRIPETAL

We need the interfacial force module as well as some functions to compute the position of the interface.

#include "iforce.h"
#include "curvature.h"

User Data

We need to define the coordinate of the point where the center of the droplet should remain p, and a parameter that controls the intensity of the suspending force eps.

struct SuspendingForceModel {
  coord p;
  double eps;
  double sigma;
};

struct SuspendingForceModel sfm = {
  .p = {0., 0., 0.},
  .eps = 1.25e-4,
  .sigma = 0.,
};

Position of an interface

This is similar to reduced but instead of computing the position of the interface we calculate the potential: \displaystyle \xi = \dfrac{\xi_0}{|\mathbf{x}_\Gamma - \mathbf{x}_c|}

This is defined only in interfacial cells. In all the other cells it takes the value nodata.

We first need a function to compute the magnitute of the distance between the interface and the suspending center |\mathbf{x}_\Gamma - \mathbf{x}_c| of an interface. For accuracy, we first try to use height functions.

foreach_dimension()
static double pos_centripetal_x (Point point, vector h)
{
  if (fabs(height(h.x[])) > 1.)
    return nodata;
  coord o = {x, y, z};
  o.x += height(h.x[])*Delta;
  double pos = 0.;
  foreach_dimension()
    pos += sq (o.x - sfm.p.x);
  return sqrt (pos);
}

We now need to choose one of the x, y or z height functions to compute the position. This is done by the function below which returns the HF position given a volume fraction field f and a height function field h.

static double height_position_centripetal (Point point, scalar f, vector h)
{

We first define pairs of normal coordinates n (computed by simple differencing of f) and corresponding HF position function pos (defined above).

  typedef struct {
    double n;
    double (* pos) (Point, vector);
  } NormPos;
  struct { NormPos x, y, z; } n;
  foreach_dimension()
    n.x.n = f[1] - f[-1], n.x.pos = pos_centripetal_x;

We sort these pairs in decreasing order of |n|.

  if (fabs(n.x.n) < fabs(n.y.n))
    swap (NormPos, n.x, n.y);
#if dimension == 3
  if (fabs(n.x.n) < fabs(n.z.n))
    swap (NormPos, n.x, n.z);
  if (fabs(n.y.n) < fabs(n.z.n))
    swap (NormPos, n.y, n.z);
#endif

We try each position function in turn.

  double pos = nodata;
  foreach_dimension()
    if (pos == nodata)
      pos = n.x.pos (point, h);

  return pos;
}

void position_centripetal (scalar f, scalar pos, bool add = false)
{

On trees we set the prolongation and restriction functions for the position.

#if TREE
  pos.refine = pos.prolongation = curvature_prolongation;
  pos.restriction = curvature_restriction;
#endif

  vector fh = f.height, h = automatic (fh);
  if (!fh.x.i)
    heights (f, h);
  foreach() {
    if (interfacial (point, f)) {
      double hp = height_position_centripetal (point, f, h);
      if (hp == nodata) {
        coord n = mycs (point, f), o = {x,y,z}, c;
        double alpha = plane_alpha (f[], n);
        plane_area_center (n, alpha, &c);

        hp = 0.;
        foreach_dimension()
          hp += sq (o.x + Delta*c.x - sfm.p.x);
        hp = sqrt (hp);
      }
      if (add)
        pos[] += 1./hp*sfm.eps*(rho2 - rho1);
      else
        pos[] = 1./hp*sfm.eps*(rho2 - rho1);
    }
    else
      pos[] = nodata;
  }
}

Stability condition

We add the possibility to compute the stability conditions based on the time step required for the time-explicit discretization of the surface tension force. This allows to have a reliable time-step when using this module without surface tension. The coefficient \sigma is set as sfm.sigma.

event stability (i++)
{

We first compute the minimum and maximum values of \alpha/f_m = 1/\rho, as well as \Delta_{min}.

  double amin = HUGE, amax = -HUGE, dmin = HUGE;
  foreach_face (reduction(min:amin) 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 = sfm.sigma;
  if (sigma) {
    double dt = sqrt (rhom*cube(dmin)/(pi*sigma));
    if (dt < dtmax)
      dtmax = dt;
  }
}

We overload the acceleration() event to add the contribution of gravity to the interfacial potential \phi.

If \phi is already allocated, we add the suspending force potential, otherwise we allocate a new field and set it to the contribution of the suspending force potential.

event acceleration (i++)
{
  scalar phi = f.phi;
  if (phi.i)
    position_centripetal (f, phi, add = true);
  else {
    phi = new scalar;
    position_centripetal (f, phi, add = false);
    f.phi = phi;
  }
}

Notes

I noticed that this model does not work properly with AMR, except if the adapt_wavelet_leave_interface function is used in order to avoid coarsening of interface cells. Maybe different refinement and coarsening functions must be used for pos?

References