Julio Guzman
Subscriber

Hi Rob,

The UDF I have written is to calculate saturation temperature given partial pressure of water and so does not involve any phase transfer values. I have diagnosed each line of the code using the Message() function (during initialisation) and everything outputted seems sensible! I haven't tracked any of them with a UDM yet!

#include "udf.h"

#define MOLAR_MASS_WATER 18.01534 //g/mol
#define MOLAR_MASS_AIR 28.97 // g/mol
#define RHO_WV 0.5542 // kg/m3
#define RHO_AIR 1.225 // kg/m3

DEFINE_PROPERTY(saturation_temp, c, ct)
{
// c: cell variable

// Cell volume
real vol = C_VOLUME(c, ct);

// Get the volume fraction of both phases
real vf_sp = C_VOF(c, sp_ct);
real vf_pp = 1 - vf_sp;

// Get the pressure of the mixture
real p_mix = C_P(c, ct);

// Get the operating pressure
real p_op = RP_Get_Real("operating-pressure");

// Primary phase density
real rho_pp = C_R(c, pp_ct);

// Get mass fractions in primary phase
real mf[2]; // Array to store mass fractions

Material *spe = NULL;
int i; // Species index - 0 for water vapor and 1 for air

mixture_species_loop(m, spe, i)
{
mf[i] = C_YI(c, pp_ct, i);
}

real p_wv; // h20 pressure for cell
real t_sat; // saturation temperature

// If secondary phase (liquid water) only
if (vf_sp == 1)
{
p_wv = p_mix + p_op; // Gauge pressure of water + operating pressure

// Inverted Tetens equation (see wikipedia)
t_sat = (4947.46 - 82.46 * log10(p_wv)) / (23.69 - 2.3 * log10(p_wv));
}

// If primary phase (moist) or mixture of phases
else
{
// Find the partial pressure of water vapour
// partial pressure = cell pressure * water mole fraction

// mass of primary phase in cell
real m_pp = rho_pp * vol * vf_pp;

// mass of water vapour and air in cell
real m_wv = mf[0] * m_pp;
real m_air = m_pp - m_wv;

// No of moles in water vapour and air
real N_wv = m_wv / MOLAR_MASS_WATER;
real N_air = m_air / MOLAR_MASS_AIR;

real N_total = N_wv + N_air;

// Water vapour partial pressure
p_wv = (C_P(c, ct) + p_op) * (N_wv / N_total);

if (mf[0] < 0.01)
{
// Match Tetens to linear interpolation from NIST data to account for low partial water pressure behaviour
t_sat = 8.64e-3 * p_wv + 273.15;
}

else
{
// Inverted Tetens equation (see wikipedia)
t_sat = (4947.46 - 82.46 * log10(p_wv)) / (23.69 - 2.3 * log10(p_wv));
}

}

return t_sat;
}

Regarding the high temperatures, I don't need to operate eventually at such high temperatures: I just wanted to check that, given solver settings etc, the simulation runs appropriately for the t_sat=const model (which it does), implying that nothing is wrong on the Solution Controls side of things.

There are two fluid domains: one for the water (so that I could patch), and one for the humid air mixture. The settings for each are:
In grey are the wall boundaries and in green the interal boundary between the two fluid domains:

Does this all look sensible to you?

Thanks for all your help so far!

Best,
Nick