sandbox/ecipriano/src/temperature-gradient.h

    Temperature Gradient Phase Change Model

    This phase change model is suitable for boiling conditions. The vaporization rate is computed from an energy balance at the interface (in every interfacial cell), assuming that the interface is at saturation temperature and that the heat conduction to the interface is used for the phase change phenomena. Therefore, the vaporization rate is computed from the following balance:

    \displaystyle \dot{m}_{Evap} \Delta h_{ev} = \mathbf{\dot{q}}_l \cdot \mathbf{n_\Gamma} + \mathbf{\dot{q}}_g \cdot \mathbf{n_\Gamma} = - \lambda_l \left.\dfrac{\partial T_l}{\partial \mathbf{n}_\Gamma}\right\vert_l - \lambda_g \left.\dfrac{\partial T_g}{\partial \mathbf{n}_\Gamma}\right\vert_g

    where \lambda is the thermal conductivity, \Delta h_{ev} is the enthalpy of evaporation and \mathbf{n_\Gamma} is the normal to the interface. After the calculation of the vaporization rate, this model solves the transport equation for the temperature field T using a two-field formulation and including the presence of the phase change, both in the diffusion equation and in the transport of the temperature field which keeps into account the Stefan convection.

    #include "intgrad.h"
#include "fracface.h"
#include "diffusion.h"

    Memory Allocations

    This phase change model defines a list of vaporization rates mEvapList which is populated by a single scalar mEvap because a single contribution to the vaporization rate is considered.

    #ifndef PHASECHANGE
scalar mEvap[];
scalar * mEvapList = {mEvap};

    fL and fG store the value of volume fractions for the calculation of the interface gradients, while the vectors fsL and fsG contain the face fraction fields computed using the fracface.h module.

    scalar fL[], fG[], f0[];
face vector fsL[], fsG[];

face vector s[];
# define PHASECHANGE
#endif

    Field Allocations

    • T one-field temperature
    • TL liquid-phase temperature field
    • TG gas-phase temperature field
    • TInt value of temperature at the interface
    scalar T[], TL[], TG[], TInt[];

    User Data

    Using this phase change model, the user should define the following variables (SI units):

    • lambda1 Thermal conductivity in liquid phase
    • lambda2 Thermal conductivity in gas phase
    • dhev Latent heat of vaporization
    • cp1 Specific heat capacity in liquid phase
    • cp2 Specific heat capacity in gas phase
    • TIntVal Interface temperature value
    extern double lambda1, lambda2, dhev, cp1, cp2;
extern double TIntVal;

    Allocating diffusivity fields lambdaf, volume correction theta and explicit source terms sgT and sTS for the diffusion equation.

    face vector lambda1f[], lambda2f[];
scalar thetacorr1[], thetacorr2[];
scalar sgT[], slT[], sgTimp[], slTimp[];

    Defaults

    In the defaults event we setup the tracer lists for the advection of the temperature fields.

    event defaults (i = 0)
{
  TL.inverse = false;
  TG.inverse = true;

#ifdef CONSISTENTPHASE1
  fuext.tracers = list_append (fuext.tracers, TL);
#else
  fu.tracers = list_append (fu.tracers, TL);
#endif
#ifdef CONSISTENTPHASE2
  fuext.tracers = list_append (fuext.tracers, TG);
#else
  fu.tracers = list_append (fu.tracers, TG);
#endif

    On adaptive meshes, tracers need to use linear interpolation (rather than the default bilinear interpolation) to ensure conservation when refining cells.

    #if TREE
#if EMBED
      TL.refine = TL.prolongation = refine_embed_linear;
      TG.refine = TG.prolongation = refine_embed_linear;
#else
      TL.refine  = refine_linear;
      TG.refine  = refine_linear;
#endif
      TL.restriction = restriction_volume_average;
      TL.dirty = true; // boundary conditions need to be updated
      TG.restriction = restriction_volume_average;
      TG.dirty = true; // boundary conditions need to be updated
#endif
}

    Init

    In the init event, we avoid dumping all the fields that we don’t need to visualize.

    event init (i = 0)
{
  sgT.nodump = true;
  slT.nodump = true;
  fL.nodump  = true;
  fG.nodump  = true;
  thetacorr1.nodump = true;
  thetacorr2.nodump = true;
}

    Finalise

    We deallocate the various lists from the memory.

    event cleanup (t = end)
{
  free (fu.tracers), fu.tracers = NULL;
  free (fuext.tracers), fuext.tracers = NULL;
}

    Reset Source Terms

    We set to zero phase change source terms, in order to add other source term contributions from outside this module.

    event reset_sources (i++)
{
  foreach() {
    sgT[] = 0.;
    slT[] = 0.;
    sgTimp[] = 0.;
    slTimp[] = 0.;
  }
}

    Phase Change

    In the phasechange event, the vaporization rate is computed and the diffusion step for the mass fraction field (in liquid and gas phase) is solved.

    event phasechange (i++)
{

    First we compute the value of the non volume-averaged temperature fields. This procedure allows a better calculation of the gradients close to the interface.

      foreach() {
    f[] = clamp (f[], 0., 1.);
    f[] = (f[] > F_ERR) ? f[] : 0.;
    f0[] = f[];
    fL[] = f[]; fG[] = 1. - f[];
    TL[] = f[] > F_ERR ? TL[]/f[] : 0.;
    TG[] = ((1. - f[]) > F_ERR) ? TG[]/(1. - f[]) : 0.;
  }

    We compute the value of volume fraction f on the cell-faces using a geometric approach (necessary for interface gradients and diffusion equations).

      face_fraction (fL, fsL);
  face_fraction (fG, fsG);

    We compute the vaporization rate from the interface jump condition, obtaining the vaporization rate per unit of interface surface, stored in mEvap.

      foreach() {
    mEvap[] = 0.; TInt[] = 0.;
    if (f[] > F_ERR && f[] < 1.-F_ERR) {
      TInt[] = TIntVal;

      double ltrgrad = ebmgrad (point, TL, fL, fG, fsL, fsG, false, TIntVal, false);
      double gtrgrad = ebmgrad (point, TG, fL, fG, fsL, fsG, true,  TIntVal, false);

      mEvap[] += lambda1*ltrgrad/dhev;
      mEvap[] += lambda2*gtrgrad/dhev;

#ifdef SOLVE_LIQONLY
      mEvap[] = lambda1*ltrgrad/dhev;
#endif

#ifdef SOLVE_GASONLY
      mEvap[] = lambda2*gtrgrad/dhev;
#endif
    }
  }

    The calculation of the interface gradients is used also for the calculation of the source terms for the diffusion equation of the temperature fields.

      foreach() {
    if (f[] > F_ERR && f[] < 1.-F_ERR) {
      coord n = facet_normal (point, fL, fsL), p;
      double alpha = plane_alpha (fL[], n);
      double area = plane_area_center (n, alpha, &p);
      normalize (&n);

      double ltrgrad = ebmgrad (point, TL, fL, fG, fsL, fsG, false, TIntVal, false);
      double gtrgrad = ebmgrad (point, TG, fL, fG, fsL, fsG, true, TIntVal, false);

      double lheatflux = lambda1*ltrgrad;
      double gheatflux = lambda2*gtrgrad;

#ifdef AXI
      slT[] += lheatflux/rho1/cp1*area*(y + p.y*Delta)/(Delta*y)*cm[];
      sgT[] += gheatflux/rho2/cp2*area*(y + p.y*Delta)/(Delta*y)*cm[];
#else
      slT[] += lheatflux/rho1/cp1*area/Delta*cm[];
      sgT[] += gheatflux/rho2/cp2*area/Delta*cm[];
#endif
    }
  }

    We restore the tracer form of the liquid and gas-phase temperature fields.

      foreach() {
    TL[] *= f[]*(f[] > F_ERR);
    TG[] *= (1. - f[])*((1. - f[]) > F_ERR);
    T[]  = TL[] + TG[];
  }
}

    Tracer Advection

    We let the volume fractions fu and fuext to advect the fields YL and YG, as implemented in the tracer_advection event of evaporation.h

    event tracer_advection (i++);

    Tracer Diffusion

    We solve the diffusion equations for the temperature fields accounting for the phase change contributions.

    event tracer_diffusion (i++)
{

    We remove the fractions of f and mass fractions lower than F_ERR and we reconstruct the non-volume averaged form of the mass fraction fields, in order to improve the discretization of the face gradients in the diffusion equation.

      foreach() {
    f[] = clamp (f[], 0., 1.);
    f[] = (f[] > F_ERR) ? f[] : 0.;
    //fL[] = f[]; fG[] = 1. - f[];
    //TL[] = fL[] > F_ERR ? TL[]/fL[] : 0.;
    //TG[] = ((1. - fL[]) > F_ERR) ? TG[]/(1. - fL[]) : 0.;

#ifdef CONSISTENTPHASE1
    TL[] = fuext[] > F_ERR ? TL[]/fuext[] : 0.;
#else
    TL[] = fu[] > F_ERR ? TL[]/fu[] : 0.;
#endif
#ifdef CONSISTENTPHASE2
    TG[] = ((1. - fuext[]) > F_ERR) ? TG[]/(1. - fuext[]) : 0.;
#else
    TG[] = ((1. - fu[]) > F_ERR) ? TG[]/(1. - fu[]) : 0.;
#endif
    fL[] = f[]; fG[] = 1. - f[];
  }

    We compute the value of volume fraction f on the cell-faces using a geometric approach (necessary for interface gradients and diffusion equations).

      face_fraction (fL, fsL);
  face_fraction (fG, fsG);

    We solve the diffusion equations, confined by means of the face fraction fields fsL and fsG.

      foreach_face() {
    lambda1f.x[] = lambda1/rho1/cp1*fsL.x[]*fm.x[];
    lambda2f.x[] = lambda2/rho2/cp2*fsG.x[]*fm.x[];
  }

  foreach() {
    thetacorr1[] = cm[]*max(fL[], F_ERR);
    thetacorr2[] = cm[]*max(fG[], F_ERR);
    //slT[] = (f[] > F_ERR) ? slT[] : 0.;
    //sgT[] = (f[] > F_ERR) ? sgT[] : 0.;
  }

  diffusion (TG, dt, D=lambda2f, r=sgT, theta=thetacorr2);
  diffusion (TL, dt, D=lambda1f, r=slT, theta=thetacorr1);

  foreach() {
    TL[] *= fL[];
    TG[] *= (1. - fL[]);
  }

    We reconstruct the one-field temperature field summing the two fields T_L and T_G in tracer form.

      foreach() {
    TL[] = (fL[] > F_ERR) ? TL[] : 0.;
    TG[] = (fG[] > F_ERR) ? TG[] : 0.;
    T[] = TL[] + TG[];
  }
}