[Sqnderby]

Dear all,

I have looked thoroughly at forums to try and find a solution for my problem, but sadly it failed.

I am trying to simulate an evaporating ammonia/water spray. I have a long rectangular duct with 970 C hot flue gas flowing upwards (approx 5 m/s) in the Y-direction. The particles are injected in the Z direction and is supposed to evaporate into both ammonia and water. The initial conditions for the spray is specified by a particle file from previous simulation. A compressed air stream(2,5bar) directly behind the injection point accelerates the particles. Radiation is included with the DO-model and the k-e standard model is used to account for turbulence. I am using the multicomponent law for the particles, which doesnt seem to take the change of particle mass and diameter into consideration (atleast when i plot the particle trajectories after reaching convergence).

In order to calculate the evaporation rate of both species, I implemented an UDF using the DEFINE_DPM_HEAT_MASS function. I use the same code as the example from the UDF folder. In order to change the diameter of the particle as it evaporates, i implemented another UDF using the DEFINE_DPM_LAW. Both codes are compiled without issues (using Visual Studio Express).

However, after convergence I notice:

• The particle temperatures exceeds the boiling point of the components even though the 1st UDF should prevent it from doing so by invoking constant temperature during boiling
• The mass and diameter of the particles remains unchanged along the trajectories, even though I can see both ammonia and water evaporating into the continous phase

I would greatly appreciate any suggestion as to why the UDF's don't have the desired impact on my simulation.

Best Regards
Simon

The two UDF's are as follows:

/************************************************** *****************
UDF the models a custom law for evaporation of species
************************************************** *****************/

#include "udf.h"
#define REMOVE_PARTICLES FALSE

DEFINE_DPM_HEAT_MASS(multivap,p,Cp,hgas,hvap,cvap_ surf,dydt,dzdt)
{
int ns;
Material *sp;

real dens_total =0.0;
real P_total =0.0;

int nc = TP_N_COMPONENTS( p ); /* number of particle components */
cell_t c0 = RP_CELL(&(p->cCell)); /* cell and thread */
Thread *t0 = RP_THREAD( &(p->cCell) ); /* where the particle is in */
Material *gas_mix = THREAD_MATERIAL( t0 ); /* gas mixture material */
Material *cond_mix = p->injection->material;/* particle mixture material */
cphase_state_t *c = &(p->cphase); /* cell info of particle location */
real molwt[MAX_SPE_EQNS]; /* molecular weight of gas species */
real Tp = P_T(p); /* particle temperature */
real mp = P_MASS(p); /* particle mass */
real molwt_bulk = 0.; /* average molecular weight in bulk gas */
real Dp = DPM_DIAM_FROM_VOL( mp / P_RHO(p) ); /* particle diameter */
real Ap = DPM_AREA(Dp); /* particle surface */
real Pr = c->sHeat * c->mu / c->tCond; /* Prandtl number */
real Nu = 2.0 + 0.6 * sqrt( p->Re ) * pow( Pr, 1./3. ); /* Nusselt number */
real h = Nu * c->tCond / Dp; /* Heat transfer coefficient */
real dh_dt = h * ( c->temp - Tp ) * Ap; /* heat source term */
dydt[0] += dh_dt / ( mp * Cp );
dzdt->energy -= dh_dt;

mixture_species_loop(gas_mix,sp,ns)
{
molwt[ns] = MATERIAL_PROP(sp,PROP_mwi); /* molecular weight of gas species */
molwt_bulk += C_YI(c0,t0,ns) / molwt[ns]; /* average molecular weight */
}


/* prevent division by zero */
molwt_bulk = MAX(molwt_bulk,DPM_SMALL);

for( ns = 0; ns < nc; ns++ )
{
/* gas species index of vaporization */
int gas_index = TP_COMPONENT_INDEX_I(p,ns);
if( gas_index >= 0 )
{
/* condensed material */
Material * cond_c = MIXTURE_COMPONENT(cond_mix, ns );
/* vaporization temperature */
real vap_temp = MATERIAL_PROP(cond_c,PROP_vap_temp);
/* diffusion coefficient */
real D = MATERIAL_PROP_POLYNOMIAL( cond_c, PROP_binary_diffusivity, c->temp);
/* Schmidt number */
real Sc = c->mu / ( c->rho * D );
/* mass transfer coefficient */
real k = ( 2. + 0.6 * sqrt(p->Re) * pow( Sc, 1./3. ) ) * D / Dp;
/* bulk gas concentration */
real cvap_bulk = c->pressure / UNIVERSAL_GAS_CONSTANT / c->temp
* c->yi[gas_index] / molwt_bulk / solver_par.molWeight[gas_index];
/* vaporization rate */
real vap_rate = k * molwt[gas_index] * Ap * ( cvap_surf[ns] - cvap_bulk );
/* only condensation below vaporization temperature */
if( 0. < vap_rate && Tp < vap_temp )
vap_rate = 0.;

dydt[1+ns] -= vap_rate;
dzdt->species[gas_index] += vap_rate;
/* dT/dt = dh/dt / (m Cp)*/
dydt[0] -= hvap[gas_index] * vap_rate / ( mp * Cp );
/* gas enthalpy source term */
dzdt->energy += hgas[gas_index] * vap_rate;

P_total += cvap_surf[ns];
dens_total += cvap_surf[ns]*molwt[gas_index];
}
}
/*multicomponent boiling*/
P_total *= UNIVERSAL_GAS_CONSTANT * Tp;

if (P_total > c->pressure && dydt[0] > 0.)
{
real h_boil = dydt[0] *mp *Cp;
dydt[0] = 0.;
for(ns=0; ns < nc; ns++)
{
int gas_index = TP_COMPONENT_INDEX_I(p,ns);
if (gas_index >=0)
{
real boil_rate = h_boil / hvap[gas_index] * cvap_surf[ns] * molwt[gas_index] / dens_total;
dydt[1+ns] -= boil_rate;
dzdt->species[gas_index] += boil_rate;
dzdt->energy += hgas[gas_index]*boil_rate;
}
}
}

}

#include "udf.h"
#define REMOVE_PARTICLES TRUE;

DEFINE_DPM_LAW(Diam,p,ci)
{
P_CURRENT_LAW(p) = DPM_LAW_MULTICOMPONENT;
{
real dm = P_MASS0(p)-P_MASS(p);
real rho = P_RHO(p);
real vol = rho/dm;
real mp = P_MASS(p);
P_DIAM(p) = P_DIAM0(p) - DPM_DIAM_FROM_VOL(vol);
if (mp <= 0 )
REMOVE_PARTICLES TRUE;
}
}