sandbox/ecipriano/src/evaporation.h

    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.

    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.

    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
    
    }