# Multicomponent Phase Change Model with Variable Properties

This is the variable-properties extension of the multicomponent.h phase change model. Eventually, the two models will be merged in a single header.

This phase change model describes the evaporation of non-isothermal systems with multiple chemical species. The interface jump condition is computed coupling the species and temperature equation in every interfacial cell, assuming that thermodynamic equilibrium conditions occur at the gas-liquid interface. The following balance is solved:

\displaystyle \begin{cases} \dot{m}_i = \dot{m}\hat \omega_{i,l} - \rho_l \mathcal{D}_{i,l} \left.\dfrac{\partial \omega_{i,l}}{\partial\mathbf{n}_\Gamma}\right\vert_l = \dot{m}\hat\omega_{i,g} - \rho_g \mathcal{D}_{i,g} \left.\dfrac{\partial \omega_{i,g}}{\partial\mathbf{n}_\Gamma}\right\vert_g \\ \sum_{j=0}^{NS} \dot{m}_i \Delta h_{ev,i} = \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 \\ \hat x_{i,g} = k_{eq,i}(\hat T) \hat x_{i,l} \\ \end{cases}

obtaining the vaporization rate \dot{m}_i for each chemical species, the liquid interface mass fractions \hat{\omega}_{i,l}, the gas interface mass fractions \hat{\omega}_{i,g}, and the interface temperature \hat{T}. model

#define PCM_MULTICOMPONENT

We include the definitions of the interface gradients for the jump condition, the face fraction for the diffusion equation, and thermodynamic functions for the equilibrium at the gas-liquid interface.

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

## Model Options

The following compiler directives can be defined to modify some options of this code:

• NGS: number of chemical species in gas-phase
• NLS: number of chemical species in liquid-phase
• SOLVE_TEMPERATURE: non-isothermal system: solve the temperature equation
• USE_CLAPEYRON: Clausius-Clapeyron equation for the thermodynamic equilibrium
• USE_ANTOINE: Antoine equation for the thermodynamic equilibrium (see thermodynamics.h)

If USE_GSL is defined, we can solve the jump condition as a non-linear system of equations. The calculation of the interface temperature from the vaporization rate also relies of this keyword:

• undefined: temperature not solved, only decoupled solution of jump conditions
• USE_GSL = 0: decoupled solution of jump condition for the chemical species. The interface temperature is obtained from the non-linear algebraic equations solver. Required when SOLVE_TEMPERATURE is defined.
• USE_GSL = 1: decoupled solution of jump condition and interface temperature to obtain fist guess values which are then refined by the non-linear algebraic equations solver.
#ifdef USE_GSL
# include "ijc-coupled-gsl.h"
#endif

## User Data

The following data should be provided by the user:

• gas_species: list of gas phase species
• liq_species: list of liquid phase species
• inert_species: name of the inert species
• gas_start: initial mass fractions of gas species
• liq_start: initial mass fractions of liquid species
• inDmix1: diffusivities of species in liquid phase
• inDmix2: diffusivities of species in gas phase
• inKeq: thermodynamic equilibrium constant
• Tboil: boiling temperature of liquid species (only if USE_CLAPEYRON)
extern char* gas_species[NGS];
extern char* liq_species[NLS];
extern char* inert_species[1];

extern double gas_start[NGS];
extern double liq_start[NLS];

extern double inDmix1[NLS];
extern double inDmix2[NGS];

extern double inKeq[NLS];
#ifdef USE_CLAPEYRON
extern double Tboil[NLS];
#endif

By default, the molecular weights are set to 1 for each species. This can be modified in an init event. The default thermodynamic reference pressure of the system is 1atm.

double inMW[NGS];
double Pref = 101325.;

If SOLVE_TEMPERATURE is defined, additional data must be defined:

• lambda1: thermal conductivity in liquid phase
• lambda2: thermal conductivity in gas phase
• dhev: enthalpy of evaporation
• cp1: specific heat for the liquid phase
• cp2: specific heat for the gas phase
• TG0, TL0: inital temperatures in gas and liquid phase
#ifdef SOLVE_TEMPERATURE
extern double lambda1, lambda2, dhev, cp1, cp2;
extern double TL0, TG0;

scalar T[], TInt[];
scalar TL, TG;
scalar sgT[], slT[], sgTimp[], slTimp[];
face vector lambda1f[], lambda2f[];
#endif

## Fields

Since we donâ€™t know a-priori how many species we want to simulate, we initialize empty lists containing the vaporization rate for every species, the mass fractions, the interface properties and phase change source terms for the transport equations.

scalar * mEvapList = NULL;    // [NGS]
scalar * YList     = NULL;    // [NGS]
scalar * YLList    = NULL;    // [NLS]
scalar * YGList    = NULL;    // [NGS]
scalar * YLIntList = NULL;    // [NLS]
scalar * YGIntList = NULL;    // [NGS]
scalar * slexpList = NULL;    // [NLS]
scalar * slimpList = NULL;    // [NLS]
scalar * sgexpList = NULL;    // [NGS]
scalar * sgimpList = NULL;    // [NGS]
#ifdef VARPROP
scalar * Dmix1List = NULL;    // [NLS]
scalar * Dmix2List = NULL;    // [NGS]
scalar * dhevList  = NULL;    // [NLS]
scalar * Cp1List   = NULL;    // [NLS]
scalar * Cp2List   = NULL;    // [NGS]
#endif
#ifdef MOLAR_DIFFUSION
scalar * XLList    = NULL;    // [NLS]
scalar * XGList    = NULL;    // [NGS]
scalar * XLIntList = NULL;    // [NLS]
scalar * XGIntList = NULL;    // [NGS]
#endif

We declare useful fields used for loops over chemical species:

• LSI: vector with Liquid Species Indices [NLS]
• GOSI: vector with Gas-Only Species Indices [NGOS]
• NGOS: number of Gas-Only Species Indices
• inertIndex: index/position of the inert species.
int * LSI = NULL;
int * GOSI = NULL;
int NGOS;
int inertIndex;

We initilize other useful fields.

bool success;
bool init_fields = true;

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

scalar divu[], fold[];

Variable properties stuff.

#ifdef VARPROP
scalar dummy[];
scalar rho1v0[], rho2v0[];
#endif

## Defaults

In the defaults event, the lists containing the species are initialized from the lists of species names in gas and liquid phase. The tracers are assigned in this event.

event defaults (i = 0)
{

Fill Liquid Species Indices LSI array.

  Array * arrLSI = array_new();
for (int ii=0; ii<NGS; ii++) {
for (int jj=0; jj<NLS; jj++) {
if (strcmp(gas_species[ii], liq_species[jj]) == 0) {
int idx = ii; array_append (arrLSI, &idx, sizeof(int));
}
}
}
LSI = (int *) array_shrink (arrLSI);

Fill Gas Only Species Indices GOSI array.

  Array * arrGOSI = array_new();
for (int ii=0; ii<NGS; ii++) {
bool thisSpeciesIsAlsoLiquid = false;
for (int jj=0; jj<NLS; jj++) {
if (strcmp(liq_species[jj], gas_species[ii]) == 0) {
thisSpeciesIsAlsoLiquid = true;
}
}
if (!thisSpeciesIsAlsoLiquid) {
int idx = ii; array_append (arrGOSI, &idx, sizeof(int));
}
}
NGOS = arrGOSI->len/sizeof(int);
GOSI = (int *) array_shrink (arrGOSI);
assert (NGOS == (NGS - NLS));

Fill index of inert species among the gas species.

  inertIndex = -1;
for (int ii=0; ii<NGS; ii++) {
if (strcmp(gas_species[ii], inert_species[0]) == 0)
inertIndex = ii;
}

We create species fields from the list of species names.

  YList     = NULL;
YLList    = NULL;
YGList    = NULL;
YGIntList = NULL;
YLIntList = NULL;
mEvapList = NULL;
slexpList = NULL;
slimpList = NULL;
sgexpList = NULL;
sgimpList = NULL;
#ifdef VARPROP
Dmix1List = NULL;
Dmix2List = NULL;
dhevList  = NULL;
Cp1List   = NULL;
Cp2List   = NULL;
#endif
#ifdef MOLAR_DIFFUSION
XLList    = NULL;
XGList    = NULL;
XLIntList = NULL;
XGIntList = NULL;
#endif

for (int jj=0; jj<NLS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "%s_L", liq_species[jj]);
a.name = strdup (name);
YLList = list_append (YLList, a);
}
reset (YLList, 0.);

for (int jj=0; jj<NGS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "%s_G", gas_species[jj]);
a.name = strdup (name);
YGList = list_append (YGList, a);
}
reset (YGList, 0.);

for (int jj=0; jj<NLS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "%s_LInt", liq_species[jj]);
a.name = strdup (name);
YLIntList = list_append (YLIntList, a);
}
reset (YLIntList, 0.);

for (int jj=0; jj<NGS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "%s_GInt", gas_species[jj]);
a.name = strdup (name);
YGIntList = list_append (YGIntList, a);
}
reset (YGIntList, 0.);

for (int jj=0; jj<NGS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "%s", gas_species[jj]);
a.name = strdup (name);
YList = list_append (YList, a);
}
reset (YList, 0.);

for (int jj=0; jj<NGS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "mEvap_%s", gas_species[jj]);
a.name = strdup (name);
mEvapList = list_append (mEvapList, a);
}
reset (mEvapList, 0.);

for (int jj=0; jj<NLS; jj++) {
scalar a = new scalar;
scalar b = new scalar;
free (a.name);
free (b.name);
char aname[20];
char bname[20];
sprintf (aname, "slexp_%s", liq_species[jj]);
sprintf (bname, "slimp_%s", liq_species[jj]);
a.name = strdup (aname);
b.name = strdup (bname);
a.nodump = true;
b.nodump = true;
slexpList = list_append (slexpList, a);
slimpList = list_append (slimpList, b);
}
reset (slexpList, 0.);
reset (slimpList, 0.);

for (int jj=0; jj<NGS; jj++) {
scalar a = new scalar;
scalar b = new scalar;
free (a.name);
free (b.name);
char aname[20];
char bname[20];
sprintf (aname, "sgexp_%s", gas_species[jj]);
sprintf (bname, "sgimp_%s", gas_species[jj]);
a.name = strdup (aname);
b.name = strdup (bname);
a.nodump = true;
b.nodump = true;
sgexpList = list_append (sgexpList, a);
sgimpList = list_append (sgimpList, b);
}
reset (sgexpList, 0.);
reset (sgimpList, 0.);

#ifdef VARPROP
for (int jj=0; jj<NLS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "Dmix1_%s", liq_species[jj]);
a.name = strdup (name);
//a.nodump = true;
Dmix1List = list_append (Dmix1List, a);
}
for (int jj=0; jj<NGS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "Dmix2_%s", gas_species[jj]);
a.name = strdup (name);
//a.nodump = true;
Dmix2List = list_append (Dmix2List, a);
}
reset (Dmix1List, 0.);
reset (Dmix2List, 0.);

for (int jj=0; jj<NLS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "dhev_%s", liq_species[jj]);
a.name = strdup (name);
//a.nodump = true;
dhevList = list_append (dhevList, a);
}
reset (dhevList, 0.);

for (int jj=0; jj<NLS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "Cp1_%s", liq_species[jj]);
a.name = strdup (name);
a.nodump = true;
Cp1List = list_append (Cp1List, a);
}
for (int jj=0; jj<NGS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "Cp2_%s", gas_species[jj]);
a.name = strdup (name);
a.nodump = true;
Cp2List = list_append (Cp2List, a);
}
reset (Cp1List, 0.);
reset (Cp2List, 0.);
#endif
#ifdef MOLAR_DIFFUSION
for (int jj=0; jj<NLS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "XL_%s", liq_species[jj]);
a.name = strdup (name);
a.nodump = true;
a.dirty = false;
XLList = list_append (XLList, a);
}
reset (XLList, 0.);

for (int jj=0; jj<NGS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "XG_%s", gas_species[jj]);
a.name = strdup (name);
a.nodump = true;
a.dirty = false;
XGList = list_append (XGList, a);
}
reset (XGList, 0.);

for (int jj=0; jj<NLS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "XLInt_%s", liq_species[jj]);
a.name = strdup (name);
a.nodump = true;
XLIntList = list_append (XLIntList, a);
}
reset (XLIntList, 0.);

for (int jj=0; jj<NGS; jj++) {
scalar a = new scalar;
free (a.name);
char name[20];
sprintf (name, "XGInt_%s", gas_species[jj]);
a.name = strdup (name);
a.nodump = true;
XGIntList = list_append (XGIntList, a);
}
reset (XGIntList, 0.);
#endif

fL.nodump = true;
fG.nodump = true;

for (scalar s in YGList)
s.inverse = true;

for (scalar s in YLList)
s.inverse = false;

#ifdef MOLAR_DIFFUSION
for (scalar s in XGList)
s.inverse = true;

for (scalar s in XLList)
s.inverse = false;
#endif

We restore to NULL the tracers associated with fu and fuext in order to re-run simulations correctly.

  fuext.tracers = NULL;
fu.tracers    = NULL;

#ifdef CONSISTENTPHASE1
fuext.tracers = list_concat (fuext.tracers, YLList);
#else
fu.tracers = list_concat (fu.tracers, YLList);
#endif
#ifdef CONSISTENTPHASE2
fuext.tracers = list_concat (fuext.tracers, YGList);
#else
fu.tracers = list_concat (fu.tracers, YGList);
#endif

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

#if TREE
for (scalar s in YGList) {
#if EMBED
s.refine = s.prolongation = refine_embed_linear;
#else
s.refine  = refine_linear;
#endif
s.restriction = restriction_volume_average;
s.dirty = true; // boundary conditions need to be updated
}
#endif

#if TREE
for (scalar s in YLList) {
#if EMBED
s.refine = s.prolongation = refine_embed_linear;
#else
s.refine  = refine_linear;
#endif
s.restriction = restriction_volume_average;
s.dirty = true; // boundary conditions need to be updated
}
#endif

#ifdef SOLVE_TEMPERATURE

TL = new scalar;
TG = new scalar;

TL.inverse = false;
TG.inverse = true;

sgT.nodump    = true;
slT.nodump    = true;
sgTimp.nodump = true;
slTimp.nodump = 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

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

Set default multicomponent properties.

  for (int jj=0; jj<NGS; jj++)
inMW[jj] = 1.;

#ifdef VARPROP
if (is_constant (rho2v)) {
scalar * l = list_copy (all);
rho2v = new scalar;
free (all);
all = list_concat ({rho2v}, l);
free (l);
}
if (is_constant (mu2v)) {
scalar * l = list_copy (all);
mu2v = new scalar;
free (all);
all = list_concat ({mu2v}, l);
free (l);
}
if (is_constant (cp2v)) {
scalar * l = list_copy (all);
cp2v = new scalar;
free (all);
all = list_concat ({cp2v}, l);
free (l);
}
if (is_constant (lambda2v)) {
scalar * l = list_copy (all);
lambda2v = new scalar;
free (all);
all = list_concat ({lambda2v}, l);
free (l);
}
if (is_constant (rho1v)) {
scalar * l = list_copy (all);
rho1v = new scalar;
free (all);
all = list_concat ({rho1v}, l);
free (l);
}
if (is_constant (mu1v)) {
scalar * l = list_copy (all);
mu1v = new scalar;
free (all);
all = list_concat ({mu1v}, l);
free (l);
}
if (is_constant (cp1v)) {
scalar * l = list_copy (all);
cp1v = new scalar;
free (all);
all = list_concat ({cp1v}, l);
free (l);
}
if (is_constant (lambda1v)) {
scalar * l = list_copy (all);
lambda1v = new scalar;
free (all);
all = list_concat ({lambda1v}, l);
free (l);
}
#endif
}

## Init

In the init event, We initialize the temperature and chemical species and set to zero additional fields.

event init (i = 0)
{
foreach() {
for (int jj=0; jj<NLS; jj++) {
scalar s = YLList[jj];
if (init_fields)
s[] = f[]*liq_start[jj];
else
s[] = 0.;
}
for (int jj=0; jj<NGS; jj++) {
scalar s = YGList[jj];
if (init_fields)
s[] = (1. - f[])*gas_start[jj];
else
s[] = 0.;
}
for (int jj=0; jj<NLS; jj++) {
scalar s = YLIntList[jj];
s[] = 0.;
}
for (int jj=0; jj<NGS; jj++) {
scalar s = YGIntList[jj];
s[] = 0.;
}
for (int jj=0; jj<NGS; jj++) {
scalar s = mEvapList[jj];
s[] = 0.;
}
for (int jj=0; jj<NGOS; jj++) {
scalar s  = YList[GOSI[jj]];
scalar sg = YGList[GOSI[jj]];
s[] = sg[];
}
for (int jj=0; jj<NLS; jj++) {
scalar s  = YList[LSI[jj]];
scalar sg = YGList[LSI[jj]];
scalar sl = YLList[jj];
s[] = sg[] + sl[];
}
for (int jj=0; jj<NLS; jj++) {
scalar slexp = slexpList[jj];
scalar slimp = slimpList[jj];
slexp[] = 0.;
slimp[] = 0.;
}
for (int jj=0; jj<NGS; jj++) {
scalar sgexp = sgexpList[jj];
scalar sgimp = sgimpList[jj];
sgexp[] = 0.;
sgimp[] = 0.;
}
for (int jj=0; jj<NGS; jj++) {
scalar mEvap = mEvapList[jj];
mEvap[] = 0.;
}
}

#ifdef SOLVE_TEMPERATURE
if (init_fields) {
foreach() {
TL[] = TL0*f[];
TG[] = TG0*(1. - f[]);
T[]  = TL[] + TG[];
}
}
#endif
}

## Cleanup

We deallocate the various lists from the memory.

event cleanup (t = end)
{
delete (mEvapList), free (mEvapList), mEvapList = NULL;
delete (YList), free (YList), YList = NULL;
delete (YLList), free (YLList), YLList = NULL;
delete (YGList), free (YGList), YGList = NULL;
delete (YLIntList), free (YLIntList), YLIntList = NULL;
delete (YGIntList), free (YGIntList), YGIntList = NULL;
delete (slexpList), free (slexpList), slexpList = NULL;
delete (slimpList), free (slimpList), slimpList = NULL;
delete (sgexpList), free (sgexpList), sgexpList = NULL;
delete (sgimpList), free (sgimpList), sgimpList = NULL;
#ifdef VARPROP
delete (Dmix1List), free (Dmix1List), Dmix1List = NULL;
delete (Dmix2List), free (Dmix2List), Dmix2List = NULL;
delete (dhevList), free (dhevList), dhevList = NULL;
delete (Cp1List), free (Cp1List), Cp1List = NULL;
delete (Cp2List), free (Cp2List), Cp2List = NULL;
#endif
#ifdef MOLAR_DIFFUSION
delete (XLList), free (XLList), XLList = NULL;
delete (XGList), free (XGList), XGList = NULL;
delete (XLIntList), free (XLIntList), XLIntList = NULL;
delete (XGIntList), free (XGIntList), XGIntList = NULL;
#endif
delete (fu.tracers), free (fu.tracers), fu.tracers = NULL;
delete (fuext.tracers), free (fuext.tracers), fuext.tracers = NULL;
free (LSI);
free (GOSI);
#ifdef SOLVE_TEMPERATURE
delete ({TL,TG});
#endif
}

## 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() {
#ifdef SOLVE_TEMPERATURE
sgT[] = 0.;
slT[] = 0.;
sgTimp[] = 0.;
slTimp[] = 0.;
#endif
for (int jj=0; jj<NLS; jj++) {
scalar slexp = slexpList[jj];
scalar slimp = slimpList[jj];
slexp[] = 0.;
slimp[] = 0.;
}
for (int jj=0; jj<NGS; jj++) {
scalar sgexp = sgexpList[jj];
scalar sgimp = sgimpList[jj];
sgexp[] = 0.;
sgimp[] = 0.;
}
}
}

## Helper functions

We define a function which is used to update the mixture molecular weight and the mole fractions (if needed).

#include "multicomponent-properties.h"

void update_mw_moles (void) {
double MW1[NLS], MW2[NGS];
foreach_elem (YLList, jj)
MW1[jj] = inMW[LSI[jj]];
foreach_elem (YGList, jj)
MW2[jj] = inMW[jj];

#ifdef MOLAR_DIFFUSION
for (scalar s in XLList)
s.dirty = false;
for (scalar s in XGList)
s.dirty = false;
#endif

foreach() {
double x1[NLS], y1[NLS];
for (int jj=0; jj<NLS; jj++) {
scalar YL = YLList[jj];
y1[jj] = (NLS == 1.) ? 1. : YL[];
}
correctfrac (y1, NLS);
mass2molefrac (x1, y1, MW1, NLS);
MW1mix[] = mass2mw (y1, MW1, NLS);
#ifdef MOLAR_DIFFUSION
for (int jj=0; jj<NLS; jj++) {
scalar XL = XLList[jj];
XL[] = x1[jj];
}
#endif

double x2[NGS], y2[NGS];
for (int jj=0; jj<NGS; jj++) {
scalar YG = YGList[jj];
y2[jj] = (NGS == 1.) ? 1. : YG[];
}
correctfrac (y2, NGS);
mass2molefrac (x2, y2, MW2, NGS);
MW2mix[] = mass2mw (y2, MW2, NGS);
#ifdef MOLAR_DIFFUSION
for (int jj=0; jj<NGS; jj++) {
scalar XG = XGList[jj];
XG[] = x2[jj];
}
#endif
}

double MW_TOL = T_PROP;

foreach() {
if (f[] <= MW_TOL) {
double MW1mixvgh = 0.;

int counter = 0;
foreach_neighbor(1) {
if (f[] > MW_TOL) {
counter++;
MW1mixvgh += MW1mix[];
}
}
MW1mix[] = (counter != 0.) ? MW1mixvgh/counter : 0.;
}
}

foreach() {
if ((1. - f[]) <= MW_TOL) {
double MW2mixvgh = 0.;

int counter = 0;
foreach_neighbor(1) {
if ((1. - f[]) > MW_TOL) {
counter++;
MW2mixvgh += MW2mix[];
}
}
MW2mix[] = (counter != 0.) ? MW2mixvgh/counter : 0.;
}
}

foreach() {
if (f[] > F_ERR) {
double x1[NLS], y1[NLS];
for (int jj=0; jj<NLS; jj++) {
scalar YL = YLList[jj];
y1[jj] = (NLS == 1) ? 1. : YL[];
}
correctfrac (y1, NLS);
mass2molefrac (x1, y1, MW1, NLS);
MW1mix[] = mass2mw (y1, MW1, NLS);
#ifdef MOLAR_DIFFUSION
for (int jj=0; jj<NLS; jj++) {
scalar XL = XLList[jj];
XL[] = x1[jj];
}
#endif
}

if ((1. - f[]) > F_ERR) {
double x2[NGS], y2[NGS];
for (int jj=0; jj<NGS; jj++) {
scalar YG = YGList[jj];
y2[jj] = YG[];
}
correctfrac (y2, NGS);
mass2molefrac (x2, y2, MW2, NGS);
MW2mix[] = mass2mw (y2, MW2, NGS);
#ifdef MOLAR_DIFFUSION
for (int jj=0; jj<NGS; jj++) {
scalar XG = XGList[jj];
XG[] = x2[jj];
}
#endif
}
}

#ifdef MOLAR_DIFFUSION
boundary (XGList);
boundary (XLList);
#endif
boundary ({MW1mix,MW2mix});

for (int b = 0; b < nboundary; b++) {
foreach_boundary (b) {

double x1[NLS], y1[NLS];
for (int jj=0; jj<NLS; jj++) {
scalar YL = YLList[jj];
y1[jj] = 0.5*(YL[] + get_ghost (point, YL, b));
}
correctfrac (y1, NLS);
mass2molefrac (x1, y1, MW1, NLS);
double MW1face = mass2mw (y1, MW1, NLS);
set_ghost (point, MW1mix, b, MW1face);
#ifdef MOLAR_DIFFUSION
for (int jj=0; jj<NLS; jj++) {
scalar XL = XLList[jj];
set_ghost (point, XL, b, x1[jj]);
}
#endif

double x2[NGS], y2[NGS];
for (int jj=0; jj<NGS; jj++) {
scalar YG = YGList[jj];
y2[jj] = 0.5*(YG[] + get_ghost (point, YG, b));
}
correctfrac (y2, NGS);
mass2molefrac (x2, y2, MW2, NGS);
double MW2face = mass2mw (y2, MW2, NGS);
set_ghost (point, MW2mix, b, MW2face);
#ifdef MOLAR_DIFFUSION
for (int jj=0; jj<NGS; jj++) {
scalar XG = XGList[jj];
set_ghost (point, XG, b, x2[jj]);
}
#endif
}
}
}

## 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[];

#ifdef SOLVE_TEMPERATURE
TL[] = f[] > F_ERR ? TL[]/f[] : 0.;
TG[] = ((1. - f[]) > F_ERR) ? TG[]/(1. - f[]) : 0.;
#endif

We recover value of mass fractions from their tracer form.

    for (int jj=0; jj<NLS; jj++) {
scalar YL = YLList[jj];
YL[] = f[] > F_ERR ? YL[]/f[] : 0.;
}

for (int jj=0; jj<NGS; jj++) {
scalar YG = YGList[jj];
YG[] = ((1. - f[]) > F_ERR) ? YG[]/(1. - f[]) : 0.;
}
}

We compute the mole fraction fields if the diffusivity is molar-based.

#ifdef MOLAR_DIFFUSION
update_mw_moles();
#endif

The thermodynamic and transport properties are updated at the beginning of the each time-step.

#ifdef VARPROP
update_properties();
#endif

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 assign the interface temperature value.

#ifdef SOLVE_TEMPERATURE
foreach() {
TInt[] = 0.;
if (f[] > F_ERR && f[] < 1.-F_ERR)
#ifdef FIXED_INTERFACE_TEMPERATURE
TInt[] = FIXED_INTERFACE_TEMPERATURE;
#else
TInt[] = avg_neighbor (point, TL, f);
#endif
}
#endif

We compute the molecular weight of the gas-only species mixture.

  scalar MWGmix[];
foreach() {
MWGmix[] = 0.;
if (f[] > F_ERR && f[] < 1.-F_ERR) {
double yGinorm[NGOS];
double sumYGi = 0.;
for (int jj=0; jj<NGOS; jj++) {
scalar YG = YGList[GOSI[jj]];
yGinorm[jj] = YG[];
sumYGi += yGinorm[jj];
}
for (int jj=0; jj<NGOS; jj++) {
yGinorm[jj] /= (sumYGi + 1.e-10);
}
for (int jj=0; jj<NGOS; jj++) {
MWGmix[] += yGinorm[jj] / inMW[GOSI[jj]];
}
MWGmix[] = 1. / (MWGmix[] + 1.e-10);
}
}

We compute total vaporization flowrate.

  foreach() {
mEvapTot[] = 0.;

We reset to zero mEvap for every species, and we set to zero the interface mass fraction fields.

    for (int jj=0; jj<NGS; jj++) {
scalar mEvap = mEvapList[jj];
scalar YGInt = YGIntList[jj];
mEvap[] = 0.;
YGInt[] = 0.;
}
for (int jj=0; jj<NLS; jj++) {
scalar YLInt = YLIntList[jj];
YLInt[] = 0.;
}

if (f[] > F_ERR && f[] < 1.-F_ERR) {

We create fields to store the local mole fractions, the conversion is required by the thermodynamic equilibrium.

      double XGIntConv[NLS+1], YGIntConv[NLS+1];
double XLIntConv[NLS], YLIntConv[NLS];
double inMWG[NLS+1], inMWL[NLS];

for (int jj=0; jj<NLS; jj++) {
inMWL[jj] = inMW[LSI[jj]];
inMWG[jj] = inMW[LSI[jj]];
XGIntConv[jj] = 0.;
YGIntConv[jj] = 0.;
XLIntConv[jj] = 0.;
YLIntConv[jj] = 0.;
}

We convert YLInt to mole fractions XLInt.

      for (int jj=0; jj<NLS; jj++) {
scalar YL    = YLList[jj];
scalar YLInt = YLIntList[jj];
YLIntConv[jj] = avg_neighbor (point, YL, f);
YLInt[] = YLIntConv[jj];
}
mass2molefrac (XLIntConv, YLIntConv, inMWL, NLS);
#ifdef MOLAR_DIFFUSION
foreach_elem (XLIntList, jj) {
scalar XLInt = XLIntList[jj];
XLInt[] = XLIntConv[jj];
}
#endif

We compute XGInt from the thermodynamic equilibrium. Different equilibrium options are available: constant thermodynamic equilibrium constant, Clausius-Clapeyron relation, and Antoine equation.

      double sumXGi = 0.;
for (int jj=0; jj<NLS; jj++) {
XGIntConv[jj] = inKeq[jj]*XLIntConv[jj];
#ifdef USE_CLAPEYRON
double dhevvh = dhev;
# ifdef VARPROP
scalar dhevjj = dhevList[jj];
dhevvh = dhevjj[];
# endif
XGIntConv[jj] = clapeyron (min (TInt[], Tboil[jj]-1.),
Tboil[jj], dhevvh, inMW[jj])*XLIntConv[jj];
#endif
#ifdef USE_ANTOINE
scalar YL = YLList[jj];
XGIntConv[jj] = min (YL.antoine (TInt[], Pref), 0.98)*XLIntConv[jj];
#endif
#ifdef USE_ANTOINE_OPENSMOKE
XGIntConv[jj] = min (opensmoke_antoine (TInt[], Pref, jj), 0.98)*XLIntConv[jj];
#endif
sumXGi += XGIntConv[jj];
}
XGIntConv[NLS] = 1. - sumXGi;
inMWG[NLS] = MWGmix[];
mole2massfrac (YGIntConv, XGIntConv, inMWG, NLS+1);

We set the gas phase interface mass fraction values using the converted fractions.

      for (int jj=0; jj<NLS; jj++) {
scalar YGInt = YGIntList[LSI[jj]];
YGInt[] = YGIntConv[jj];
#ifdef MOLAR_DIFFUSION
scalar XGInt = XGIntList[LSI[jj]];
XGInt[] = XGIntConv[jj];
#endif
}

We adjust the interface mass fractions of the gas-only species in the system.

      double yGinorm[NGOS];
double sumYGi = 0., sumYGl = 0.;
if (f[] > F_ERR && f[] < 1.-F_ERR) {
for (int jj=0; jj<NLS; jj++) {
sumYGl += YGIntConv[jj];
}
for (int jj=0; jj<NGOS; jj++) {
scalar YG = YGList[GOSI[jj]];
yGinorm[jj] = YG[];
sumYGi += yGinorm[jj];
}
for (int jj=0; jj<NGOS; jj++) {
scalar YGInt = YGIntList[GOSI[jj]];
YGInt[] = (1. - sumYGl)*yGinorm[jj]/(sumYGi + 1.e-10);
}
}

We repeat the same operations to close the molar fractions of the gas-only species in the system.

#ifdef MOLAR_DIFFUSION
double xGinorm[NGOS];
double sumXGigo = 0., sumXGl = 0.;
if (f[] > F_ERR && f[] < 1.-F_ERR) {
for (int jj=0; jj<NLS; jj++) {
sumXGl += XGIntConv[jj];
}
for (int jj=0; jj<NGOS; jj++) {
scalar XG = XGList[GOSI[jj]];
xGinorm[jj] = XG[];
sumXGigo += xGinorm[jj];
}
for (int jj=0; jj<NGOS; jj++) {
scalar XGInt = XGIntList[GOSI[jj]];
XGInt[] = (1. - sumXGl)*xGinorm[jj]/(sumXGigo + 1.e-10);
}
}
#endif

We compute the sum of the diffusive fluxes in gas phase, to be used for the Fick corrected approach.

      double jGtot = 0.;
#ifdef FICK_CORRECTED
for (int jj=0; jj<NGS; jj++) {
double rho2vh = rho2;
double inDmix2vh = inDmix2[jj];
# ifdef VARPROP
rho2vh = rho2v[];
scalar Dmix2v = Dmix2List[jj];
inDmix2vh = Dmix2v[];
# endif
# ifdef MOLAR_DIFFUSION
scalar XGInt = XGIntList[jj];
scalar XG    = XGList[jj];
double gtrgrad = ebmgrad (point, XG, fL, fG, fsL, fsG, true, XGInt[], &success);
gtrgrad *= (MW2mix[] > 0.) ? inMW[jj]/MW2mix[] : 0.;
# else
scalar YGInt = YGIntList[jj];
scalar YG    = YGList[jj];
double gtrgrad = ebmgrad (point, YG, fL, fG, fsL, fsG, true, YGInt[], &success);
# endif
}
#endif

We compute the total vaporization rate from the sum of the interface mass balances over all the chemical species in liquid phase.

      double sum_jG = 0., sum_YGInt = 0.;

for (int jj=0; jj<NLS; jj++) {
//scalar YGInt = YGIntList[LSI[jj]];
//scalar YLInt = YLIntList[jj];
//scalar YL    = YLList[jj];
//scalar YG    = YGList[LSI[jj]];

//YLInt[] = avg_neighbor (point, YL, f);
//YGInt[] = YGIntConv[jj];

double rho2vh = rho2;
double inDmix2vh = inDmix2[LSI[jj]];
#ifdef VARPROP
rho2vh = rho2v[];
scalar Dmix2v = Dmix2List[LSI[jj]];
inDmix2vh = Dmix2v[];
#endif

scalar YGInt = YGIntList[LSI[jj]];
#ifdef MOLAR_DIFFUSION
scalar XGInt = XGIntList[LSI[jj]];
scalar XG    = XGList[LSI[jj]];
double gtrgrad = ebmgrad (point, XG, fL, fG, fsL, fsG, true, XGInt[], &success);
gtrgrad *= (MW2mix[] > 0.) ? inMW[LSI[jj]]/MW2mix[] : 0.;
#else
scalar YG    = YGList[LSI[jj]];
double gtrgrad = ebmgrad (point, YG, fL, fG, fsL, fsG, true, YGInt[], &success);
#endif
sum_YGInt += YGInt[];
}
#ifdef DIFFUSIVE
mEvapTot[] = sum_jG;
#else
mEvapTot[] = sum_jG/min(1. - sum_YGInt, 0.99);
#endif
}
}

From the knowledge of the total vaporization rate, we compute the vaporization rate for each species.

  foreach() {
if (f[] > F_ERR && f[] < 1.-F_ERR) {

double jGtot = 0.;
#ifdef FICK_CORRECTED
for (int jj=0; jj<NGS; jj++) {
double rho2vh = rho2;
double inDmix2vh = inDmix2[jj];
# ifdef VARPROP
rho2vh = rho2v[];
scalar Dmix2v = Dmix2List[jj];
inDmix2vh = Dmix2v[];
# endif

# ifdef MOLAR_DIFFUSION
scalar XGInt = XGIntList[jj];
scalar XG    = XGList[jj];
double gtrgrad = ebmgrad (point, XG, fL, fG, fsL, fsG, true, XGInt[], &success);
gtrgrad *= (MW2mix[] > 0.) ? inMW[jj]/MW2mix[] : 0.;
# else
scalar YGInt = YGIntList[jj];
scalar YG    = YGList[jj];
double gtrgrad = ebmgrad (point, YG, fL, fG, fsL, fsG, true, YGInt[], &success);
# endif
}
#endif

for (int jj=0; jj<NLS; jj++) {
scalar mEvap = mEvapList[LSI[jj]];
scalar YGInt = YGIntList[LSI[jj]];

double rho2vh = rho2;
double inDmix2vh = inDmix2[LSI[jj]];
#ifdef VARPROP
rho2vh = rho2v[];
scalar Dmix2v = Dmix2List[LSI[jj]];
inDmix2vh = Dmix2v[];
#endif

#ifdef MOLAR_DIFFUSION
scalar XGInt = XGIntList[LSI[jj]];
scalar XG    = XGList[LSI[jj]];
double gtrgrad = ebmgrad (point, XG, fL, fG, fsL, fsG, true, XGInt[], &success);
gtrgrad *= (MW2mix[] > 0.) ? inMW[LSI[jj]]/MW2mix[] : 0.;
#else
scalar YG    = YGList[LSI[jj]];
double gtrgrad = ebmgrad (point, YG, fL, fG, fsL, fsG, true, YGInt[], &success);
#endif

#ifdef DIFFUSIVE
mEvap[] = - rho2vh*inDmix2vh*gtrgrad - YGInt[]*jGtot;
#else
mEvap[] = mEvapTot[]*YGInt[] - rho2vh*inDmix2vh*gtrgrad - YGInt[]*jGtot;
#endif
}
}
}

If SOLVE_TEMPERATURE, the value of interface temperature is obtained from the total vaporization rate just computed and the interface temperature gradients. This result in an algebraic equation which is numerically zero-ed.

#ifdef USE_GSL

# ifdef SOLVE_TEMPERATURE
ijc_CoupledTemperature();
# endif

Finally, we can refine the first guess values from the decoupled solution from the fully-coupled interface jump condition. Solved as a non-linear system of equations.

# if USE_GSL > 0
ijc_CoupledNls();
# endif

#endif

The source terms for the diffusion equation of the species mass fractions in gas an liquid phase are computed here.

  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);

for (int jj=0; jj<NGS; jj++) {
scalar sgexp = sgexpList[jj];
scalar sgimp = sgimpList[jj];
scalar mEvap = mEvapList[jj];

#ifdef AXI
sgexp[] += -mEvap[]*area*(y + p.y*Delta)/(Delta*y)*cm[];
sgimp[] += +mEvapTot[]*area*(y + p.y*Delta)/(Delta*y)*cm[];
#else
sgexp[] += -mEvap[]*area/Delta*cm[];
sgimp[] += +mEvapTot[]*area/Delta*cm[];
#endif
}

for (int jj=0; jj<NLS; jj++) {
scalar slexp = slexpList[jj];
scalar slimp = slimpList[jj];
scalar mEvap = mEvapList[LSI[jj]];

#ifdef AXI
slexp[] += +mEvap[]*area*(y + p.y*Delta)/(Delta*y)*cm[];
slimp[] += -mEvapTot[]*area*(y + p.y*Delta)/(Delta*y)*cm[];
#else
slexp[] += +mEvap[]*area/Delta*cm[];
slimp[] += -mEvapTot[]*area/Delta*cm[];
#endif
}

#ifdef SOLVE_TEMPERATURE
double bc = TInt[];
double gtrgrad = ebmgrad (point, TG, fL, fG, fsL, fsG, true, bc, &success);
double ltrgrad = ebmgrad (point, TL, fL, fG, fsL, fsG, false, bc, &success);

double lambda1vh = lambda1;
double lambda2vh = lambda2;
# ifdef VARPROP
lambda1vh = lambda1v[];
lambda2vh = lambda2v[];
# endif

# ifndef VARPROP  // Without VARPROP this form converges better
lheatflux /= (rho1*cp1);
gheatflux /= (rho2*cp2);
# endif

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

We include the velocity correction for the diffusive fluxes. Using the FICK_CORRECTED approach this method enforces mass conservation in multicomponent diffusion. Using MOLAR_DIFFUSION, this procedure includes a missing term in the implicit diffusion step, which is a function of the mixture molecular weight.

#ifdef MASS_DIFFUSION_ENTHALPY
foreach() {
double mde1 = 0., mde2 = 0.;
coord gTL = {0., 0., 0.,}, gTG = {0., 0., 0.};
coord gYLj = {0., 0., 0.}, gYGj = {0., 0., 0.};
coord gYLsum = {0., 0., 0.}, gYGsum = {0., 0., 0.};

foreach_dimension() {
gTG.x = (TG[1] - TG[-1])/(2.*Delta);
gTL.x = (TL[1] - TL[-1])/(2.*Delta);
}

foreach_dimension() {
foreach_elem (YGList, jj) {
scalar Dmix2v = Dmix2List[jj];
# ifdef MOLAR_DIFFUSION
scalar XG = XGList[jj];
gYGsum.x -= (MW2mix[] > 0.) ?
rho2v[]*Dmix2v[]*inMW[jj]/MW2mix[]*(XG[1] - XG[-1])/(2.*Delta) : 0.;
# else
scalar YG = YGList[jj];
gYGsum.x -= rho2v[]*Dmix2v[]*(YG[1] - YG[-1])/(2.*Delta);
# endif
}

foreach_elem (YLList, jj) {
scalar Dmix1v = Dmix1List[jj];
# ifdef MOLAR_DIFFUSION
scalar XL = XLList[jj];
gYLsum.x -= (MW1mix[] > 0.) ?
rho1v[]*Dmix1v[]*inMW[LSI[jj]]/MW1mix[]*(XL[1] - XL[-1])/(2.*Delta) : 0.;
# else
scalar YL = YLList[jj];
gYLsum.x -= rho1v[]*Dmix1v[]*(YL[1] - YL[-1])/(2.*Delta);
# endif
}

foreach_elem (YGList, jj) {
scalar YG = YGList[jj];
scalar Cp2v = Cp2List[jj];
scalar Dmix2v = Dmix2List[jj];
# ifdef MOLAR_DIFFUSION
scalar XG = XGList[jj];
gYGj.x = (MW2mix[] > 0.) ?
-rho2v[]*Dmix2v[]*inMW[jj]/MW2mix[]*(XG[1] - XG[-1])/(2.*Delta) : 0.;
# else
gYGj.x = -rho2v[]*Dmix2v[]*(YG[1] - YG[-1])/(2.*Delta);
# endif
mde2 += Cp2v[]*(gYGj.x - YG[]*gYGsum.x)*gTG.x;
}

foreach_elem (YLList, jj) {
scalar YL = YLList[jj];
scalar Cp1v = Cp1List[jj];
scalar Dmix1v = Dmix1List[jj];
# ifdef MOLAR_DIFFUSION
scalar XL = XLList[jj];
gYLj.x = (MW1mix[] > 0.) ?
-rho1v[]*Dmix1v[]*inMW[LSI[jj]]/MW1mix[]*(XL[1] - XL[-1])/(2.*Delta) : 0.;
# else
gYLj.x = -rho1v[]*Dmix1v[]*(YL[1] - YL[-1])/(2.*Delta);
# endif
mde1 += Cp1v[]*(gYLj.x - YL[]*gYLsum.x)*gTL.x;
}
}
slT[] -= mde1*cm[]*(f[] > (1. - F_ERR));
sgT[] -= mde2*cm[]*(f[] < F_ERR);
}
#endif

The divergence of the velocity field is computed after the calculation of source terms for species and temperature.

#ifdef VARPROP
update_divergence();
#endif

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

  foreach() {
for (scalar YL in YLList)
YL[] *= f[]*(f[] > F_ERR);
for (scalar YG in YGList)
YG[] *= (1. - f[])*((1. - f[]) > F_ERR);
#ifdef SOLVE_TEMPERATURE
TL[] *= f[]*(f[] > F_ERR);
TG[] *= (1. - f[])*((1. - f[]) > F_ERR);
#endif
}
}

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 species and temperature 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.;

for (scalar YL in YLList)
YL[] = (fuext[] > F_ERR) ? YL[]/fuext[] : 0.;
for (scalar YG in YGList)
YG[] = ((1. - fu[]) > F_ERR) ? YG[]/(1. - fu[]) : 0.;

fL[] = f[]; fG[] = 1. - f[];

#ifdef SOLVE_TEMPERATURE
TL[] = (fuext[] > F_ERR) ? TL[]/fuext[] : 0.;
TG[] = ((1. - fu[]) > F_ERR) ? TG[]/(1. - fu[]) : 0.;
#endif
}

#ifdef MOLAR_DIFFUSION
update_mw_moles();
#endif

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 include the velocity correction for the diffusive fluxes. Using the FICK_CORRECTED approach this method enforces mass conservation in multicomponent diffusion. Using MOLAR_DIFFUSION, this procedure includes a missing term in the implicit diffusion step, which is a function of the mixture molecular weight.

  face vector phicGtot[];
foreach_face() {
phicGtot.x[] = 0.;
double rho2f = 0.5*(rho2v[] + rho2v[-1]);
for (int jj=0; jj<NGS; jj++) {
#ifdef VARPROP
scalar Dmix2 = Dmix2List[jj];
double Dmix2f = 0.5*(Dmix2[] + Dmix2[-1]);
#else
double Dmix2f = inDmix2[jj];
#endif

#ifdef FICK_CORRECTED
# ifdef MOLAR_DIFFUSION
scalar XG = XGList[jj];
double MW2mixf = 0.5*(MW2mix[] + MW2mix[-1]);
phicGtot.x[] += (MW2mixf > 0.) ?
# else
scalar YG = YGList[jj];
# endif // MOLAR_DIFFUSION
#else
phicGtot.x[] = 0.;
#endif  // FICK_CORRECTED
}
phicGtot.x[] *= fsG.x[]*fm.x[];
}

face vector phicLtot[];
foreach_face() {
phicLtot.x[] = 0.;
double rho1f = 0.5*(rho1v[] + rho1v[-1]);
for (int jj=0; jj<NLS; jj++) {
#ifdef VARPROP
scalar Dmix1 = Dmix1List[jj];
double Dmix1f = 0.5*(Dmix1[] + Dmix1[-1]);
#else
double Dmix1f = inDmix1[jj];
#endif

#ifdef FICK_CORRECTED
# ifdef MOLAR_DIFFUSION
scalar XL = XLList[jj];
double MW1mixf = 0.5*(MW1mix[] + MW1mix[-1]);
phicLtot.x[] += (MW1mixf > 0.) ?
# else
scalar YL = YLList[jj];
# endif // MOLAR_DIFFUSION
#else
phicLtot.x[] = 0.;
#endif  // FICK_CORRECTED
}
phicLtot.x[] *= fsL.x[]*fm.x[];
}

for (int jj=0; jj<NGS; jj++) {
face vector phicjj[];
foreach_face() {
phicjj.x[] = phicGtot.x[];
#ifdef MOLAR_DIFFUSION
scalar Dmix2 = Dmix2List[jj];
double Dmix2f = 0.5*(Dmix2[] + Dmix2[-1]);
double rho2f = 0.5*(rho2v[] + rho2v[-1]);
double MW2mixf = 0.5*(MW2mix[] + MW2mix[-1]);

phicjj.x[] -= (MW2mixf > 0.) ?
#endif  // MOLAR_DIFFUSION
}
scalar YG = YGList[jj];
face vector flux[];
tracer_fluxes (YG, phicjj, flux, dt, zeroc);

foreach()
foreach_dimension()
YG[] += (rho2v[] > 0.) ? dt/(rho2v[])*(flux.x[] - flux.x[1])/(Delta*cm[]) : 0.;
}

for (int jj=0; jj<NLS; jj++) {
face vector phicjj[];
foreach_face() {
phicjj.x[] = phicLtot.x[];
#ifdef MOLAR_DIFFUSION
scalar Dmix1 = Dmix1List[jj];
double Dmix1f = 0.5*(Dmix1[] + Dmix1[-1]);
double rho1f = 0.5*(rho1v[] + rho1v[-1]);
double MW1mixf = 0.5*(MW1mix[] + MW1mix[-1]);

phicjj.x[] -= (MW1mixf > 0.) ?
#endif  // MOLAR_DIFFUSION
}
scalar YL = YLList[jj];
face vector flux[];
tracer_fluxes (YL, phicjj, flux, dt, zeroc);

foreach()
foreach_dimension()
YL[] += (rho1v[] > 0.) ? dt/(rho1v[])*(flux.x[] - flux.x[1])/(Delta*cm[]) : 0.;
}

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

  scalar theta1[], theta2[];

#if TREE
theta1.refine = theta1.prolongation = fraction_refine;
theta2.refine = theta2.prolongation = fraction_refine;
theta1.dirty = true;
theta2.dirty = true;
#endif

for (int jj=0; jj<NLS; jj++) {
face vector Dmix1f[];
foreach_face() {
#ifdef VARPROP
scalar Dmix1jj = Dmix1List[jj];
double Dmix1vh = 0.5*(Dmix1jj[] + Dmix1jj[-1]);
double rho1vh = 0.5*(rho1v[] + rho1v[-1]);
Dmix1f.x[] = rho1vh*Dmix1vh*fsL.x[]*fm.x[];
#else
Dmix1f.x[] = rho1*inDmix1[jj]*fsL.x[]*fm.x[];
//Dmix1f.x[] = inDmix1[jj]*fsL.x[]*fm.x[];
#endif
}

foreach()
#ifdef VARPROP
theta1[] = cm[]*max(rho1v[]*fL[], F_ERR);
#else
theta1[] = cm[]*max(rho1*fL[], F_ERR);
#endif

scalar YL = YLList[jj];
scalar slexp = slexpList[jj];
scalar slimp = slimpList[jj];

//foreach() {
//  slexp[] = (f[] > F_ERR) ? slexp[] : 0.;
//  slimp[] = (f[] > F_ERR) ? slimp[] : 0.;
//}

diffusion (YL, dt, D=Dmix1f, r=slexp, beta=slimp, theta=theta1);
}

for (int jj=0; jj<NGS; jj++) {

face vector Dmix2f[];
foreach_face() {
#ifdef VARPROP
scalar Dmix2v = Dmix2List[jj];
double Dmix2vh = 0.5*(Dmix2v[] + Dmix2v[-1]);
double rho2vh = 0.5*(rho2v[] + rho2v[-1]);
Dmix2f.x[] = rho2vh*Dmix2vh*fsG.x[]*fm.x[];
#else
Dmix2f.x[] = rho2*inDmix2[jj]*fsG.x[]*fm.x[];
#endif
}

foreach()
#ifdef VARPROP
theta2[] = cm[]*max(rho2v[]*fG[], F_ERR);
#else
theta2[] = cm[]*max(rho2*fG[], F_ERR);
#endif

scalar YG = YGList[jj];
scalar sgexp = sgexpList[jj];
scalar sgimp = sgimpList[jj];

//foreach() {
//  sgexp[] = (f[] > F_ERR) ? sgexp[] : 0.;
//  sgimp[] = (f[] > F_ERR) ? sgimp[] : 0.;
//}

diffusion (YG, dt, D=Dmix2f, r=sgexp, beta=sgimp, theta=theta2);
}

#ifdef SOLVE_TEMPERATURE

foreach_face() {
#ifdef VARPROP
//double alpha1l = (rho1v[] != 0.) ? lambda1v[]/rho1v[]/cp1v[] : 0.;
//double alpha2l = (rho2v[] != 0.) ? lambda2v[]/rho2v[]/cp2v[] : 0.;
//double alpha1r = (rho1v[-1] != 0.) ? lambda1v[-1]/rho1v[-1]/cp1v[-1] : 0.;
//double alpha2r = (rho2v[-1] != 0.) ? lambda2v[-1]/rho2v[-1]/cp2v[-1] : 0.;
//lambda1f.x[] = 0.5*(alpha1r + alpha1l)*fsL.x[]*fm.x[];
//lambda2f.x[] = 0.5*(alpha2r + alpha2l)*fsG.x[]*fm.x[];
lambda1f.x[] = 0.5*(lambda1v[] + lambda1v[-1])*fsL.x[]*fm.x[];
lambda2f.x[] = 0.5*(lambda2v[] + lambda2v[-1])*fsG.x[]*fm.x[];
#else
lambda1f.x[] = lambda1/rho1/cp1*fsL.x[]*fm.x[];
lambda2f.x[] = lambda2/rho2/cp2*fsG.x[]*fm.x[];
#endif
}

foreach() {
#ifdef VARPROP
theta1[] = cm[]*max(fL[]*rho1v[]*cp1v[], F_ERR);
theta2[] = cm[]*max(fG[]*rho2v[]*cp2v[], F_ERR);
#else
theta1[] = cm[]*max(fL[], F_ERR);
theta2[] = cm[]*max(fG[], F_ERR);
#endif
}

Solve diffusion equations for temperature.

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

#endif

foreach() {

We reconstruct the mass fractions from the total mass in every cell.

    double totmassliq = 0.;
for (scalar YL in YLList)
totmassliq += YL[];
for (scalar YL in YLList)
YL[] = (totmassliq > 0.) ? YL[]/totmassliq : 0.;

double totmassgas = 0.;
for (scalar YG in YGList)
totmassgas += YG[];
for (scalar YG in YGList)
YG[] = (totmassgas > 0.) ? YG[]/totmassgas : 0.;

for (scalar YL in YLList)
YL[] *= f[];
for (scalar YG in YGList)
YG[] *= (1. - f[]);
#ifdef SOLVE_TEMPERATURE
TL[] *= f[];
TG[] *= (1. - f[]);
#endif
}

We reconstruct the volume-averaged mass fractions and the one-field mass fraction field.

  foreach() {
for (scalar YL in YLList) {
YL[] = clamp (YL[], 0., 1.);
YL[] = (YL[] > F_ERR) ? YL[] : 0.;
}
for (scalar YG in YGList) {
YG[] = clamp (YG[], 0., 1.);
YG[] = (YG[] > F_ERR) ? YG[] : 0.;
}
for (int jj=0; jj<NGOS; jj++) {
scalar Y  = YList[GOSI[jj]];
scalar YG = YGList[GOSI[jj]];

Y[] = YG[];
}
for (int jj=0; jj<NLS; jj++) {
scalar Y  = YList[LSI[jj]];
scalar YG = YGList[LSI[jj]];
scalar YL = YLList[jj];

Y[] = YL[] + YG[];
}
}

#ifdef SOLVE_TEMPERATURE
foreach()
T[] = TL[] + TG[];
#endif
}