two-phase.h

Go to the documentation of this file.

/** @file two-phase.h
 */
/**
# A compressible two-phase flow solver
 
This solves the two-phase compressible Navier-Stokes equations
including the total energy equation.
\f[
\frac{\partial (f \rho_i)}{\partial t } + 
\nabla \cdot (f \rho_i \mathbf{u}) = 0 
\f]
\f[
\frac{\partial (\rho_i \mathbf{u})}{\partial t } + 
\nabla \cdot ( \rho_i  \mathbf{u} \mathbf{u}) = 
-\nabla p + \nabla \cdot \tau_i' 
\f]
\f[
\frac{\partial [\rho_i (e_i + \mathbf{u}^2/2)]}{\partial t } 
+ \nabla \cdot [ \rho_i \mathbf{u} (e_i  + \mathbf{u}^2/2)] =
-\nabla \cdot (\mathbf{u} p_i) + 
\nabla \cdot \left( \tau'_i \mathbf{u} \right)
\f]
an advection equation for the volume fraction \f$f\f$
\f[
\frac{\partial f}{\partial t} + \mathbf{u} \cdot \nabla f = 0
\f]
and an Equation Of State (EOS) that needs to be defined. 
 
The method is described in detail in [Fuster & Popinet,
2018](#fuster2018) and relies on a time splitting technique were we
first solve the advection step using the VOF method for \f$f\f$, \f$f_i
\rho_i\f$, \f$f_i \rho_i \mathbf{u}\f$ \f$f_i \rho_i e_{tot,i}\f$ and then
perform the projection step using the all-Mach solver. */
 
#include "all-mach.h"
 
/**
We transport VOF tracers without the one-dimensional compressive term. */
 
#define NO_1D_COMPRESSION 1
#include "vof.h"
 
/**
The primary fields are:
\f[
\begin{aligned}
\text{frho1} & = f \rho_1, \\
\text{frho2} & = (1-f) \rho_2, \\
\text{q} & = f \rho_1 \mathbf{u} + (1-f) \rho_2 \mathbf{u}, \\
\text{fE1} & = f \rho_1 (e_1 + \mathbf{u}^2/2), \\
\text{fE2} & = (1-f) \rho_2 (e_2 + \mathbf{u}^2/2)
\end{aligned}
\f]
*/
 
scalar f[], * interfaces = {f};
scalar frho1[], frho2[], fE1[], fE2[];
 
/**
The timestep can be limited by a CFL based on the speed of sound. */
 
double CFLac = HUGE;
 
/**
The dynamic viscosities for each phase, as well as the volumetric
viscosity coefficients. */
 
double mu1 = 0., mu2 = 0.;
double lambdav1 = 0., lambdav2 = 0. ;
 
/**
These functions are provided by the Equation Of State. See for example
[the Mie--Gruneisen Equation of State](Mie-Gruneisen.h). */
 
extern double sound_speed          (Point point);
extern double average_pressure     (Point point);
extern double bulk_compressibility (Point point);
extern double internal_energy      (Point point, double fc);
 
/**
By default the Harmonic mean is used to compute the phase-averaged
dynamic viscosity. */
 
#ifndef mu
# define mu(f) (mu1*mu2/(clamp(f,0,1)*(mu2 - mu1) + mu1))
#endif
 
/**
The volumetric viscosity uses arithmetic average by default. */
 
// fixme: Include the term depending on \f$\nabla \cdot u\f$. It is tricky because this
// quantity can be very different upon the phase
#ifndef lambdav
# define lambdav(f)  (clamp(f,0,1)*(lambdav1 - lambdav2) + lambdav2)
// # define lambdav(f)  (clamp(f,0.,1.)*(lambdav1 - lambdav2 - 2./3*(mu1 - mu2)) + lambdav2 - 2./3*mu2)
#endif
 
/**
## Auxilliary fields
 
Auxilliary fields need to be allocated. The quantity \f$\rho c^2\f$, the
average density `rhov` and its inverse `alphav`. */
 
scalar rhoc2v[];
vector alphav[];
scalar rhov[];
const vector lambdav0[] = {0,0,0};
const vector lambdav = lambdav0;
 
/**
## Time step restriction based on the speed of sound */
 
/** @brief Event: stability (i++) */
void event_stability (void) 
{
  if (CFLac < HUGE)
    foreach (reduction (min:dtmax)) {
      double c = sound_speed (point);
      if (CFLac*Delta < c*dtmax)
  dtmax = CFLac*Delta/c;
    }
}
 
#if TREE
/**
## Energy refinement function
 
The energy is refined from the refined pressures, momentum and
densities, using the equation of state. */
void fE_refine (Point point, scalar fE)
{
  for (int _c = 0; _c < 4; _c++) /* foreach_child */ {
    double Ek = 0.;
    for (int _d = 0; _d < dimension; _d++)
      Ek += sq(q.x[]);
    Ek /= 2.*(frho1[] + frho2[]);
    fE[] = fE.inverse ?
      (internal_energy (point, 0.) + Ek)*(1. - f[]) :
      (internal_energy (point, 1.) + Ek)*f[];
  }
}
#endif // TREE
 
/**
## Initialisation
 
We set the default values. */
 
/** @brief Event: defaults (i = 0) */
void event_defaults (void) 
{
  alpha = alphav;
  rho = rhov;
 
  /**
  If the viscosity is non-zero, we need to allocate the face-centered
  viscosity field. */
  
  if (mu1 || mu2) {
    mu = {0} /* new vector */;
    lambdav = {0} /* new vector */;
  }
 
  rhoc2 = rhoc2v;  
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    frho1[] = 1., fE1[] = 1.;
    frho2[] = 0., fE2[] = 0.;
    p[] = 1.;
    f[] = 1.;
  }
  
  /** 
  We also initialize the list of tracers to be advected with the VOF
  function \f$f\f$ (or its complementary function). */
  
  f.tracers = list_copy ({frho1, frho2, fE1, fE2});
  for (int _s = 0; _s < 1; _s++) /* scalar in {frho2, fE2} */
    s.inverse = true;
 
  /**
  We set limiting. */
  
  for (int _s = 0; _s < 1; _s++) /* scalar in {frho1, frho2, fE1, fE2, q} */ {
    s.gradient = minmod2;
#if TREE
    /**
    On trees, we ensure that limiting is also applied to prolongation
    and refinement. */
 
    s.prolongation = s.refine = refine_linear;
#endif
  }
  
  /**
  We add the interface and the density to the default display. */
 
  display ("draw_vof (c = 'f');");
  display ("squares (color = 'rhov', spread = -1);");
}
 
/** @brief Event: init (i = 0) */
void event_init (void) 
{
  /* trash: uf */;
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    uf.x[] = fm.x[]*(q.x[]/(frho1[] + frho2[]) + q.x[-1]/(frho1[-1] + frho2[-1]))/2.;
 
  /**
  We update the fluid properties. */
 
  event ("properties");
 
  /**
  We set the initial timestep (this is useful only when restoring from
  a previous run). */
 
  dtmax = DT;
  event ("stability");
 
  /** 
  For the associated tracers we use the gradient defined by
  f.gradient. */
 
  if (f.gradient)
    for (int _s = 0; _s < 1; _s++) /* scalar in {frho1, frho2, fE1, fE2, q} */
      s.gradient = f.gradient;  
}
 
/**
## VOF advection of momentum 
 
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[]/(frho1[] + frho2[]);
      q1.x[] = frho1[]*u;
      q2.x[] = frho2[]*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[];
    
    /**
    We avoid negative densities and energies which may have been
    caused by round-off during VOF advection. */
    
    if (f[] <= 0.){
      frho1[] = 0.;
      fE1[] = 0.;
    }
    if (f[] >= 1.){
      frho2[] = 0.;
      fE2[] = 0.;
    }
  }
  
#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);
  }
 
#if 0
  /**
  This is switched off by default for now as the standard refinement
  in [/src/vof.h#vof_concentration_refine]() seems to work fine. */
  
  for (int _s = 0; _s < 1; _s++) /* scalar in {fE1,fE2} */ {
    s.refine = fE_refine;
    set_prolongation (s, fE_refine);
  }
#endif
#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;
}
 
/**
## Pressure and density
 
During the projection step we compute the provisional pressure from
the EOS using the values after advection. */
 
/** @brief Event: properties (i++) */
void event_properties (void) 
{
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    rhov[] = frho1[] + frho2[];
    ps[] = average_pressure (point);
 
    /** 
    We also compute \f$\rho c^2\f$. */
    
    rhoc2v[] = bulk_compressibility (point);
  }
  
  for (int _i = 0; _i < _N; _i++) /* foreach_face */ {
    
    /** 
    If viscosity is present we obtain the averaged viscosity at the cell faces. */
  
    if (mu1 || mu2) {
      vector muv = mu, lambdavv = lambdav;
      double ff = (f[] + f[-1])/2.;
      muv.x[] = fm.x[]*mu(ff);
      lambdavv.x[] = lambdav(ff);
    }
    
    /** 
    We also compute the averaged density at the cell faces. */
    
    alphav.x[] = 2.*fm.x[]/(rhov[] + rhov[-1]);
 
  }
 
  /**
  The all-Mach solver needs rho*cm. */
  
  for (int _i = 0; _i < _N; _i++) /* foreach */
    rhov[] *= cm[];
  
  /* I still miss a source term in the divergence related to viscous
   * dissipation. Usually it is small and a little bit complicated to
   * calculate. In the current form of allmach it cannot be added
   * because it does not allow user-defined divergence sources. */
}
 
/**
## Update of the total energy
 
After projection we update the values of the total energy adding the
term missing from the projection step.
 
For a Newtonian fluid, the stress tensor comprises the pressure term,
the viscous uncompressible term and the viscous compressible term,
 
\f[ \tau = -p \mathbf{I} + \mu (\nabla \mathbf{u} + \nabla \mathbf{u}^T)
+ \lambda_v  (\nabla \cdot \mathbf{u}) \mathbf{I}
\f]
 
where \f$\mathbf{I}\f$ is the unity tensor and \f$\lambda_v = \mu_v - 2/3
\mu\f$. The volumetric viscosity \f$\mu_v\f$ is negligible in most cases. */
 
/** @brief Event: end_timestep (i++) */
void event_end_timestep (void) 
{
#if 0
  /**
  We first compute the divergence of the velocity at cell centers
  using the cell face velocity. Note that the vector *uf*
  incorporates the metric factor. */
 
  scalar divU[];
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    divU[] = 0.;
    for (int _d = 0; _d < dimension; _d++)
      divU[] += (uf.x[1] - uf.x[])/(Delta*cm[]);
  }
 
  /**
  We compute explicitly the contribution of the compressible viscous
  term to the momentum equation and the total energy equation,
  \f[
  \partial_t \mathbf{q} = \nabla \cdot [\lambda_v(\nabla \cdot \mathbf{u}) \mathbf{I}] 
  \f]
  \f[
  \partial_t [\rho (e + \mathbf{u}^2/2)] = 
  \nabla \cdot [\lambda_v (\nabla \cdot \mathbf{u}) \mathbf{u}] 
  \f]
 
  The axysimmetric case allows some simplifications. Since
  \f[
  \mathbf{S} = \lambda_v(\nabla \cdot \mathbf{u}) \mathbf{I} =
  \left( 
  \begin{array}{ccc}
  S_{xx} & 0 &  0 \       \
  0 &  S_{yy} &  0 \        \
  0 & 0 & S_{\theta \theta}  
  \end{array}
  \right)
  \f]
  with \f$S = S_{xx} = S_{yy} = S_{\theta \theta}= \lambda_v \nabla
  \cdot \mathbf{u}\f$. Hence, the radial component of \f$\nabla \cdot
  \mathbf{S}\f$ is simply
  \f[
  (\nabla \cdot \mathbf{S}) \cdot \mathbf{e}_y  = 
  \frac{1}{y} \frac{\partial(y S_{yy}) }{\partial y} 
  - \frac{S_{\theta \theta}}{y} 
  = \frac{\partial S }{\partial y}
  \f]
  Also \f$u_\theta = 0\f$. */ 
 
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    double momentum = 0., energy = 0.;
    for (int _d = 0; _d < dimension; _d++) {
      double right = lambdav.x[1]*(divU[1]  + divU[])/2.;
      double left = lambdav.x[]*(divU[-1] + divU[])/2.;
      momentum += right - left;
      energy += uf.x[1]*right - uf.x[]*left;
    }
    for (int _d = 0; _d < dimension; _d++)
      q.x[] += dt*momentum/Delta;
    energy *= dt/(Delta*cm[]);
    double fc = clamp(f[],0,1);
    fE1[] += energy*fc;
    fE2[] += energy*(1. - fc);
  }
#endif
  
  /**
  The contribution of the pressure to the energy of each phase is
  lacking. */
 
  {
    vector upf[];
    for (int _i = 0; _i < _N; _i++) /* foreach_face */
      upf.x[] = - uf.x[]*(p[] + p[-1])/2.;
 
    for (int _i = 0; _i < _N; _i++) /* foreach */ {
      double energy = 0.; 
      for (int _d = 0; _d < dimension; _d++)
  energy += upf.x[1] - upf.x[];
      energy *= dt/(Delta*cm[]);
      double fc = clamp(f[],0,1);
      fE1[] += energy*fc;
      fE2[] += energy*(1. - fc);
    }
  }
  
  /**
  This is the contribution of the incompressible viscous term. */
 
  if (mu1 || mu2) {
    
    /**
    The velocity \f$\mathbf{u}\f$ is first obtained from the updated
    momentum \f$\mathbf{q}\f$. */
    
    vector u = q;
    for (int _i = 0; _i < _N; _i++) /* foreach */
      for (int _d = 0; _d < dimension; _d++) 
        u.x[] = q.x[]/(frho1[] + frho2[]);
 
    /**
    We add the incompressible contribution of the \f$\nabla \cdot (u_i
    \tau'_i)\f$ term. Note that the compressible contribution, which is
    typically small, is neglected. */
 
    // fixme: add the compressible contribution (careful, \f$\nabla \cdot
    // u\f$ can be very different upon the phase and also very different
    // from the value defined for the mixture with an averaged velocity
    vector ueijk[];
    for (int _d = 0; _d < dimension; _d++) {
      for (int _i = 0; _i < _N; _i++) /* foreach_face */
  ueijk.x[] = 2.*uf.x[]*(u.x[] - u.x[-1])/Delta;
#if dimension > 1
      for (int _i = 0; _i < _N; _i++) /* foreach_face */
  ueijk.y[] = uf.y[]*(u.x[] - u.x[0,-1] + 
          (u.y[1,-1] + u.y[1,0])/4. -
          (u.y[-1,-1] + u.y[-1,0])/4.)/Delta;
#endif
#if dimension > 2
      for (int _i = 0; _i < _N; _i++) /* foreach_face */
  ueijk.z[] = uf.z[]*(u.x[] - u.x[0,0,-1] + 
          (u.z[1,0,-1] + u.z[1,0,0])/4. -
          (u.z[-1,0,-1] + u.z[-1,0,0])/4.)/Delta;
#endif
      for (int _i = 0; _i < _N; _i++) /* foreach */ {
  double energy = 0.; 
  for (int _d = 0; _d < dimension; _d++)
    energy += ueijk.x[1] - ueijk.x[];
  energy *= dt/(Delta*cm[]);
  double fc = clamp(f[],0,1);
  fE1[] += mu1*energy*fc;
  fE2[] += mu2*energy*(1. - fc);
      }   
 
      // fixme: Formally, we need to include a correction term due to
      // the normal viscous stress jump (to be done)
    }
 
    /**
    Finally we recover momentum.*/
 
    for (int _i = 0; _i < _N; _i++) /* foreach */
      for (int _d = 0; _d < dimension; _d++) 
        q.x[] *= frho1[] + frho2[];
  }
}
 
/**
At the end of the simulation we clean the tracer list. */
 
/** @brief Event: cleanup (i = end) */
void event_cleanup (void) {
  free (f.tracers);
  f.tracers = NULL;
}
 
/**
## References
 
~~~bib
@hal{fuster2018, hal-01845218}
~~~
 
## See also
 
* [The Mie--Gruneisen Equation of State](Mie-Gruneisen.h)
*/