src/navier-stokes/conserving.h

Momentum-conserving advection of velocity

This file implements momentum-conserving VOF advection of the velocity components for the two-phase Navier–Stokes solver.

On trees, we first define refinement and restriction functions which guarantee conservation of each component of the total momentum. Note that these functions do not guarantee conservation of momentum for each phase.

#if TREE
static void momentum_refine (Point point, scalar u) {
  refine_bilinear (point, u);
  double rhou = 0.;
  foreach_child()
    rhou += cm[]*ρ(f[])*u[];
  double du = u[] - rhou/((1 << dimension)*cm[]*ρ(f[]));
  foreach_child()
    u[] += du;
}

static void momentum_restriction (Point point, scalar u)
{
  double rhou = 0.;
  foreach_child()
    rhou += cm[]*ρ(f[])*u[];
  u[] = rhou/((1 << dimension)*cm[]*ρ(f[]));
}
#endif // TREE

We switch-off the default advection scheme of the centered solver.

event defaults (i = 0) {
  stokes = true;

#if TREE

On trees, the refinement and restriction functions above rely on the volume fraction field f being refined/restricted before the components of velocity. To ensure this, we move f to the front of the field list (all).

  int i = 0;
  while (all[i].i != f.i) i++;
  while (all[i].i)
    all[i] = all[i-1], i--;
  all[i] = f;

We then set the refinement and restriction functions for the components of the velocity field.

  foreach_dimension() {
    u.x.refine = u.x.prolongation = momentum_refine;
    u.x.restriction = momentum_restriction;
  }
#endif
}

We need to overload the stability event so that the CFL is taken into account (because we set stokes to true).

event stability (i++)
  dtmax = timestep (uf, dtmax);

We will transport the two components of the momentum, q1=fρ1u and q2=(1f)ρ2u. We will need to impose boundary conditions which match this definition. This is done using the functions below.

foreach_dimension()
static double boundary_q1_x (Point neighbor, Point point, scalar q1)
{
  return clamp(f[],0.,1.)*rho1*u.x[];
}

foreach_dimension()
static double boundary_q2_x (Point neighbor, Point point, scalar q2)
{
  return (1. - clamp(f[],0.,1.))*rho2*u.x[];
}

Similarly, on trees we need prolongation functions which also follow this definition.

#if TREE
foreach_dimension()
static void prolongation_q1_x (Point point, scalar q1) {
  foreach_child()
    q1[] = clamp(f[],0.,1.)*rho1*u.x[];
}

foreach_dimension()
static void prolongation_q2_x (Point point, scalar q2) {
  foreach_child()
    q2[] = (1. - clamp(f[],0.,1.))*rho2*u.x[];
}
#endif

We overload the vof() event to transport consistently the volume fraction and the momentum of each phase.

static scalar * interfaces1 = NULL;

event vof (i++) {

We allocate two temporary vector fields to store the two components of the momentum and set the boundary conditions and prolongation functions.

  vector q1[], q2[];
  for (scalar s in {q1,q2})
    foreach_dimension()
      s.v.x.i = -1; // not a vector
  for (int i = 0; i < nboundary; i++)
    foreach_dimension() {
      q1.x.boundary[i] = boundary_q1_x;
      q2.x.boundary[i] = boundary_q2_x;
    }
#if TREE
  foreach_dimension() {
    q1.x.prolongation = prolongation_q1_x;
    q2.x.prolongation = prolongation_q2_x;
  }
#endif

We split the total momentum q into its two components q1 and q2 associated with f and 1f respectively.

  foreach()
    foreach_dimension() {
      double fc = clamp(f[],0,1);
      q1.x[] = fc*rho1*u.x[];
      q2.x[] = (1. - fc)*rho2*u.x[];
    }
  boundary ((scalar *){q1,q2});

Momentum q2 is associated with 1f, so we set the inverse attribute to true. We use the same slope-limiting as for the velocity field.

  foreach_dimension() {
    q2.x.inverse = true;
    q1.x.gradient = q2.x.gradient = u.x.gradient;
  }

We associate the transport of q1 and q2 with f and transport all fields consistently using the VOF scheme.

  f.tracers = (scalar *){q1,q2};
  vof_advection ({f}, i);

We recover the advected velocity field using the total momentum and the density

  foreach()
    foreach_dimension()
      u.x[] = (q1.x[] + q2.x[])/ρ(f[]);
  boundary ((scalar *){u});

We set the list of interfaces to NULL so that the default vof() event does nothing (otherwise we would transport f twice).

  interfaces1 = interfaces, interfaces = NULL;
}

We set the list of interfaces back to its default value.

event tracer_advection (i++) {
  interfaces = interfaces1;
}

Usage

Examples