momentum.h

Go to the documentation of this file.

/** @file momentum.h
 */
/**
# Momentum-conserving formulation for two-phase interfacial flows
 
The interface between the fluids is tracked with a Volume-Of-Fluid
method. The volume fraction in fluid 1 is \f$f=1\f$ and \f$f=0\f$ in fluid
2. The densities and dynamic viscosities for fluid 1 and 2 are *rho1*,
*mu1*, *rho2*, *mu2*, respectively. */
 
#include "all-mach.h"
#include "vof.h"
 
scalar f[], * interfaces = {f};
double rho1 = 1., mu1 = 0., rho2 = 1., mu2 = 0.;
 
/**
Auxilliary fields are necessary to define the (variable) specific
volume \f$\alpha=1/\rho\f$ and average viscosity \f$\mu\f$ (on faces) as well
as the cell-centered density. */
 
vector alphav[], muv[];
scalar rhov[];
 
/** @brief Event: defaults (i = 0) */
void event_defaults (void) {
  alpha = alphav;
  rho = rhov;
  mu = muv;
 
  /**
  We use (strict) minmod slope limiting for all components. */
 
  theta = 1.;
  for (int _d = 0; _d < dimension; _d++)
    q.x.gradient = minmod2;
}
 
/**
The density and viscosity are defined using arithmetic averages by
default. The user can overload these definitions to use other types of
averages (i.e. harmonic). */
 
#ifndef rho
# define rho(f) (clamp(f,0,1)*(rho1 - rho2) + rho2)
#endif
#ifndef mu
# define mu(f)  (clamp(f,0,1)*(mu1 - mu2) + mu2)
#endif
 
/** @brief Event: properties (i++) */
void event_properties (void) {
  for (int _i = 0; _i < _N; _i++) /* foreach */
    rhov[] = rho(f[])*cm[];
  for (int _i = 0; _i < _N; _i++) /* foreach_face */ {
    double ff = (f[] + f[-1])/2.;
    alphav.x[] = fm.x[]/rho(ff);
    muv.x[] = fm.x[]*mu(ff);
  }
}
 
/**
We overload the *vof()* event to transport consistently the volume
fraction and the momentum of each phase. */
 
static scalar * interfaces1 = NULL;
 
/** @brief Event: vof (i++) */
void event_vof (void) {
 
  /**
  We split the total momentum \f$q\f$ into its two components \f$q1\f$ and
  \f$q2\f$ associated with \f$f\f$ and \f$1 - f\f$ respectively. */
  
  vector q1 = q, q2[];
  for (int _i = 0; _i < _N; _i++) /* foreach */
    for (int _d = 0; _d < dimension; _d++) {
      double u = q.x[]/rho(f[]);
      q1.x[] = f[]*rho1*u;
      q2.x[] = (1. - f[])*rho2*u;
    }
 
  /**
  Momentum \f$q2\f$ is associated with \f$1 - f\f$, so we set the *inverse*
  attribute to *true*. We use the same limiting for q1 and q2. */
 
  for (int _d = 0; _d < dimension; _d++) {
    q2.x.inverse = true;
    q2.x.gradient = q1.x.gradient;
  }
 
#if TREE
  /**
  The refinement function is modified by *vof_advection()*. To be able
  to restore it, we store its value. */
  
  void (* refine) (Point, scalar) = q1.x.refine;
#endif
  
  /**
  We associate the transport of \f$q1\f$ and \f$q2\f$ with \f$f\f$ and transport
  all fields consistently using the VOF scheme. */
 
  scalar * tracers = f.tracers;
  f.tracers = list_concat (tracers, (scalar *){q1, q2});
  vof_advection ({f}, i);
  free (f.tracers);
  f.tracers = tracers;
  
  /**
  We recover the total momentum. */
  
  for (int _i = 0; _i < _N; _i++) /* foreach */
    for (int _d = 0; _d < dimension; _d++)
      q.x[] = q1.x[] + q2.x[];
 
#if TREE
  /**
  We restore the refinement function for the total momentum. */
  
  for (int _s = 0; _s < 1; _s++) /* scalar in {q} */ {
    s.refine = refine;
    set_prolongation (s, refine);
  }
#endif // TREE  
 
  /**
  We set the list of interfaces to NULL so that the default *vof()*
  event does nothing (otherwise we would transport \f$f\f$ twice). */
  
  interfaces1 = interfaces, interfaces = NULL;
}
 
/**
We set the list of interfaces back to its default value. */
 
/** @brief Event: tracer_advection (i++) */
void event_tracer_advection (void) {
  interfaces = interfaces1;
}
 
/**
## See also
 
* [Two-phase interfacial flows](two-phase.h)
*/