evaporation.h

Main Evaporation Module

Main Evaporation Module

Include this file every time you need to perform simulations with phase change. It contains the declaration of fields that are useful to every phase change model, it computes the total evaporation rate and the phase change velocity, it performs the stefan flow shifting, and it manages the advection of volume fraction and tracers. This file must be used in combination with a phase change model (i.e. multicomponent.h, temperature-gradient.h, species-gradient.h, fixedflux.h).

#include "common-evaporation.h"
#include "fracface.h"
#include "diffusion.h"

Setup tracers advection

CONSISTENTPHASE1 (default): the liquid phase tracers are advected using the velocity \mathbf{u}^E, associated with the list fuext.tracers, while the gas phase tracers are advected using \mathbf{u} associated with the list fu.tracers.
CONSISTENTPHASE2: the gas phase tracers are advected using the velocity \mathbf{u}^E, associated with the list fuext.tracers, while the liquid phase tracers are advected using \mathbf{u} using the list fu.tracers.
DIFFUSIVE: the Stefan flow is neglected and no interface velocity jump is present. Therefore, both the liquid and the gas phase tracers are transported using \mathbf{u}^E which will be equal to \mathbf{u}, using the list fuext.tracers.

#ifndef DIFFUSIVE
# ifndef CONSISTENTPHASE2
#   ifndef NOCONSISTENTPHASE
#     define CONSISTENTPHASE1
#   endif
# endif
#else
# define CONSISTENTPHASE1
# define CONSISTENTPHASE2
#endif

Setup Stefan flow

SHIFTING (default): the volume expansion term is shifted.
NOSHIFTING: do not apply shifting procedure.
SHIFT_TO_LIQ (default): shifting toward the liquid phase (i.e. droplet evaporation).
SHIFT_TO_GAS: shifting toward the gas phase (i.e. boiling problems).

#define SHIFTING
#ifndef SHIFT_TO_GAS
# define SHIFT_TO_LIQ
#endif

#ifdef NOSHIFTING
# undef SHIFTING
#endif

Fields

We initialize fields useful to any phase change model:

mEvapTot: total vaporization rate per unit of interface surface [kg/m2/s]
stefanflow: expansion term used in the projection step of navier-stokes/centered-evaporation.h
vpc: interface regression velocity [m/s]
fu: dummy volume fraction used for the advection of scalar fields using the directionally-split procedure, with the velocity \mathbf{u}
fuext: dummy volume fraction used for the advection of scalar fields using the directionally-split procedure, with the velocity \mathbf{u}^E
uf_save: helper field used to store and recover a velocity before and afer calling vof_advection

scalar mEvapTot[];
scalar stefanflow[];
face vector vpc[];
scalar fu[], fuext[];
face vector uf_save[];

Extern fields

We initialize fields that must be defined externally. In particular, mEvapList should be declared in the specific phase change model, and it must contains a list of every contribution to the phase change (for example, the vaporization rate of each chemical species).

extern scalar * mEvapList;
extern face vector ufext;

We add a bool that allows the droplet volume changes to be disabled.

bool is_shrinking;
double smoothing_expansion = 0.;

Set Variable-Properties fields for compatibility with constant-properties

simulations.

#ifndef VARPROP
//(const) scalar rho1v = unity, rho2v = unity;
//(const) scalar mu1v = zeroc, mu2v = zeroc;

// The (const) scalar approach gives problems when calling those variables
// as external by modules included before evaporation.h
scalar rho1v[], rho2v[];
scalar mu1v[], mu2v[];
#endif

Init event

We activate the shrinking and we set the dummy volume fractions to the value of the vof field.

event init (i = 0)
{
  is_shrinking = true;

  foreach() {
    fu[] = f[];
    fuext[] = f[];
  }

  scalar * interfaces2 = {fu,fuext};
#if TREE
  for (scalar c in interfaces2) {
    c.refine = c.prolongation = fraction_refine;
    c.dirty = true;
    scalar * tracers = c.tracers;
    for (scalar t in tracers) {
      t.restriction = restriction_volume_average;
      t.refine = t.prolongation = vof_concentration_refine;
      t.dirty = true;
      t.c = c;
    }
  }
#endif
  for (scalar c in interfaces2) {
    scalar * tracers = c.tracers;
    for (scalar t in tracers)
      t.depends = list_add (t.depends, c);
  }

#ifndef VARPROP
  //const scalar rho1vv[] = rho1, rho2vv[] = rho2;
  //rho1v = rho1vv, rho2v = rho2vv;

  //const scalar mu1vv[] = mu1, mu2vv[] = mu2;
  //mu1v = mu1vv, mu2v = mu2vv;

  foreach() {
    rho1v[] = rho1;
    rho2v[] = rho2;
    mu1v[] = mu1;
    mu2v[] = mu2;
  }
#endif
}

We add an event that can be overloaded by the user in order to change some initialization parameters before starting the phase change event (for example, initialize chemical species and temperature to a specific initial solution).

event end_init (i = 0);

Reset Source Terms

We impose the position of the even which resets the source terms for the governing equations solved in the specific phase change models.

event reset_sources (i++);

Chemical Reaction Step

If combustion chemistry is present we solve it here.

event chemistry (i++);

Phase Change event

We computes the total vaporization rate per unit of surface by summing all the contributions in mEvapList. We interpolate the cell-centered normals and the mEvap field in order to compute a face change velocity vpc, defined at the edges.

event phasechange (i++)
{

We restore the value of the fu and fuext fields at the beginning of the time step.

  foreach() {
    fu[] = f[];
    fuext[] = f[];
  }

We compute the total vaporization mass-flowrate mEvapTot.

  foreach() {
    mEvapTot[] = 0.;
    for (scalar mEvap in mEvapList)
      mEvapTot[] += mEvap[];
  }

We update phase-change regression velocity vpc.

  foreach_face() {
    vpc.x[] = 0.;
    coord nf;
    foreach_dimension()
      nf.x = 0.;
    double mEvapf = 0.;
    if (f[] > 0. && f[] < 1.) {
      nf = normal (point, f);
      mEvapf = mEvapTot[];
      if (f[-1] > 0. && f[-1] < 1.) {
        coord np = normal (neighborp(-1), f);
        foreach_dimension()
          nf.x = (nf.x + np.x)/2.;
        mEvapf = (mEvapf + mEvapTot[-1])/2.;
      }
    }
    else if (f[-1] > 0. && f[-1] < 1.) {
      mEvapf = mEvapTot[-1];
      nf = normal (neighborp(-1), f);
    }
#ifdef BYRHOGAS
    vpc.x[] = (rho2v[] > 0.) ? fm.x[]*(-mEvapf/rho2v[])*nf.x : 0.;
#else
    vpc.x[] = (rho1v[] > 0.) ? fm.x[]*(-mEvapf/rho1v[])*nf.x : 0.;
#endif
  }

We compute the expansion source term for the continuity equation stefanflow.

  foreach() {
    stefanflow[] = 0.;
#ifndef DIFFUSIVE
    double alpha = 0.;
    double segment = 0.;
    if (f[] > F_ERR && f[] < 1.-F_ERR) {
      coord m = mycs (point, f);
      alpha = plane_alpha (f[], m);
      coord prel;
      segment = plane_area_center (m, alpha, &prel);
#ifdef AXI
      stefanflow[] = (rho1v[] > 0. && rho2v[] > 0.) ?
        cm[]*mEvapTot[]*segment*(y + prel.y*Delta)*(1./rho2v[] - 1./rho1v[])/(Delta*y)
        : 0.;
#else
      stefanflow[] = (rho1v[] > 0. && rho2v[] > 0.) ?
        mEvapTot[]*segment*(1./rho2v[] - 1./rho1v[])/Delta*cm[]
        : 0.;
#endif
    }
#endif
  }
}

If SHIFTING is defined, we shift the expansion source term.

event shifting (i++) {

We perform shifting operations of the expansion term, in order to remove it from the interfacial cells, but still keeping its effect on the velocity field.

#ifdef SHIFTING
# ifdef SHIFT_TO_LIQ
  shift_field (stefanflow, f, 1);
# else
  shift_field (stefanflow, f, 0);
# endif
#endif

To avoid an expansion term localized on a single layer of cells, we solve a diffusion equation which smoothen the expansion, and we solve the Projection step using this modified term. The volume integral of the expansion during this step must remain constant.

  if (smoothing_expansion) {
    face vector eps[];
    foreach_face()
      eps.x[] = smoothing_expansion*fm.x[];

    scalar thetacorr[];
    foreach()
      thetacorr[] = cm[];

    diffusion (stefanflow, dt, D=eps);
  }
}

VOF event

We overload the vof event in order to include the phase-change velocity. Therefore, we perform the advection of the volume fraction and we restore the original velocity field.

event vof (i++)
{
  foreach_face() {
    uf_save.x[] = uf.x[];
#ifndef INT_USE_UF
    uf.x[] = ufext.x[];
#endif
#ifndef CONSTANT_VOLUME
    if (is_shrinking)
      uf.x[] -= vpc.x[];
#endif
  }
}

Tracer Advection event

We use the directionally-split advection procedure in order to resolve the convective transport of the scalar fields, through the use of the dummy tracers fu and fuext which exploit the velocity \mathbf{u}_f and \mathbf{u}_f^E respectively.

event tracer_advection (i++)
{

We restore the value of the \mathbf{uf} velocity field.

  foreach_face()
    uf.x[] = uf_save.x[];

If the phase tracers are not advected with the same velocity of the volume fraction field, the advection is performed here exploiting the geometric advection procedure but using a dummy volume fraction field fu and fuext.

  foreach_face() {
    uf_save.x[] = uf.x[];
    uf.x[] = ufext.x[];
  }

We call the vof_advection function to transport the volume fraction fuext and associated tracers.

  vof_advection ({fuext}, i);

#ifdef VELOCITY_JUMP
  foreach_face()
# ifdef BOILING_SETUP
    uf.x[] = uf1.x[];
# else
    uf.x[] = uf2.x[];
# endif
#else
  foreach_face()
    uf.x[] = uf_save.x[];
#endif

We call the vof_advection function to transport the volume fraction fu and associated tracers.

  vof_advection ({fu}, i);

#ifdef VELOCITY_JUMP
  foreach_face()
    uf.x[] = uf_save.x[];
#endif

}