#include "udf.h"
/Constants used in psat_h2o to calculate saturation pressure/
#define PSAT_A 0.01
#define PSAT_TP 338.15
#define C_LOOP 8
#define H2O_PC 22.089E6
#define H2O_TC 647.286
/user inputs/
#define MAX_SPE_EQNS_PRIM 2 /total number of species in primary phase/
#define index_evap_primary 0 /evaporating species index in primary phase/
#define prim_index 0 /index of primary phase/
#define P_OPER 101325 /operating pressure equal to GUI value/
/end of user inputs/
//
/ UDF for specifying an interfacial area density /
//
double psat_h2o(double tsat)
/* /
/ Computes saturation pressure of water vapor /
/ as function of temperature /
/ Equation is taken from THERMODYNAMIC PROPERTIES IN SI, /
/ by Reynolds, 1979 /
/ Returns pressure in PASCALS, given temperature in KELVIN /
{
int i;
double var1,sum1,ans1,psat;
double constants[8]={-7.4192420, 2.97221E-1, -1.155286E-1,
8.68563E-3, 1.094098E-3, -4.39993E-3, 2.520658E-3, -5.218684E-4};
/ var1 is an expression that is used in the summation loop /
var1 = PSAT_A(tsat-PSAT_TP);
/* Compute summation loop /
i = 0;
sum1 = 0.0;
while (i < C_LOOP){
sum1+=constants[i]pow(var1,i);
++i;
}
ans1 == sum1(H2O_TC/tsat-1.0);
/ compute exponential to determine result /
/ psat has units of Pascals /
psat = H2O_PCexp(ans1);
return psat;
}
DEFINE_HET_RXN_RATE(user_evap_condens_react, c, t, hr, mw, yi, rr, rr_t)
{
Thread **pt = THREAD_SUB_THREADS(t);
Thread *tp = pt[0];
Thread *ts = pt[1];
int i;
real concentration_evap_primary, accum = 0., mole_frac_evap_prim,
concentration_sat ;
real T_prim = C_T(c,tp); /primary phase (gas) temperature/
real T_sec = C_T(c,ts); /secondary phase (droplet) temperature/
real diam = C_PHASE_DIAMETER(c,ts); /secondary phase diameter/
real D_evap_prim = C_DIFF_EFF(c,tp,index_evap_primary)
0.7C_MU_T(c,tp)/C_R(c,tp);
/primary phase species turbulent diffusivity/
real Re, Sc, Nu, urel, urelx,urely,urelz=0., mass_coeff, area_density,
flux_evap ;
if(Data_Valid_P())
{
urelx = C_U(c,tp) - C_U(c,ts);
urely = C_V(c,tp) - C_V(c,ts);
#if RP_3D
urelz = C_W(c,tp) - C_W(c,ts);
#endif
urel = sqrt(urelxurelx + urelyurely + urelzurelz);
/relative velocity/
Re = urel * diam * C_R(c,tp) / C_MU_L(c,tp);
Sc = C_MU_L(c,tp) / C_R(c,tp) / D_evap_prim ;
Nu = 2. + 0.6 * pow(Re, 0.5)* pow(Sc, 0.333);
mass_coeff = Nu * D_evap_prim / diam ;
for (i=0; i < MAX_SPE_EQNS_PRIM ; i++)
{
accum = accum + C_YI(c,tp,i)/mw[i][prim_index];
}
mole_frac_evap_prim = C_YI(c,tp,index_evap_primary)
/ mw[index_evap_primary][prim_index] / accum;
concentration_evap_primary = mole_frac_evap_prim * P_OPER
/ UNIVERSAL_GAS_CONSTANT / T_prim ;
concentration_sat = psat_h2o(T_sec)/UNIVERSAL_GAS_CONSTANT/T_sec ;
area_density = 6. * C_VOF(c,ts) / diam ;
flux_evap = mass_coeff *
(concentration_sat - concentration_evap_primary) ;
*rr = area_density * flux_evap ;
}
}
