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.

    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.;
        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;
        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);

    We try each position function in turn.

      double pos = nodata;
        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;
      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.;
              hp += sq (o.x + Delta*c.x - sfm.p.x);
            hp = sqrt (hp);
          if (add)
            pos[] += 1./hp*sfm.eps*(rho2v[] - rho1v[]);
            pos[] = 1./hp*sfm.eps*(rho2v[] - rho1v[]);
          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.

    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;


    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?



    Abd Essamade Saufi, Alessio Frassoldati, T Faravelli, and A Cuoci. Dropletsmoke++: a comprehensive multiphase cfd framework for the evaporation of multidimensional fuel droplets. International Journal of Heat and Mass Transfer, 131:836–853, 2019.


    Stéphane Popinet. Numerical Models of Surface Tension. Annual Review of Fluid Mechanics, 50:49 – 75, January 2018. [ DOI | http | .pdf ]