 ## Fluids

Topics relate to Fluent, CFX, Turbogrid and more

• Haotian Wu
Subscriber

I'm working on an evaporating law for my DPM simulation, but I'm currently stuck on how to call fuel mass fraction on droplet surface with macros. I did not find it in User Manual. Can somebody please help me?

• Rob
Ansys Employee

There isn't a droplet surface as such: it's all "droplet". The droplet mass fraction is what you're after, so look in the DPM macros.

• Haotian Wu
Subscriber

Hi Rob,

Thanks for replying. The droplet mass fraction is unity since the droplet has only one species. Let me explain the problem more clearly. The model I want to use included the Spalding mass transfer number which is calculated by the mass fraction of droplet species in the gas-vapor mixture surrounding the droplet. Is there a macro for that value?

Regards,

Haotian

• Rob
Ansys Employee

Yes, it's the species mass fraction, https://ansyshelp.ansys.com/account/Secured?returnurl=/Views/Secured/corp/v231/en/flu_udf/flu_udf_DataAccessMacros.html%23x1-4610004.2.3   Click on Help in the solver & paste the link into that browser.

The harder part is getting the cell that the particle is in to do that maths. That's covered in the DPM part, https://ansyshelp.ansys.com/account/Secured?returnurl=/Views/Secured/corp/v231/en/flu_udf/flu_udf_DPMDEFINE.html

• Haotian Wu
Subscriber

I have found the macro but the UDF code was causing a bad termination.

Could you please give some instructions on the code?

#include "udf.h"

DEFINE_DPM_LAW(Boiling_Law, p, coupled)
{
Material* m = P_MATERIAL(p);
cell_t c = P_CELL(p);
const cphase_state_t* c0 = &(p->cphase);

int i = P_EVAP_SPECIES_INDEX(p);

real Tp = P_T(p);

real P_sat = DPM_VAPOR_PRESSURE(p, m, Tp);

real T_ambi = C_T(c, t0);
real Tb = DPM_BOILING_TEMPERATURE(p, m);//SATURATED_TEMPERATURE
real T_spe = Tb + (T_ambi - Tb) / 3;

real rho = C_R(c, t0) + (C_R(c, t0) - DPM_RHO(p, m, Tb)) / 3;

real Ap = M_PI * P_DIAM0(p) * P_DIAM0(p);//DPM_AREA(Dp)
real LTb = DPM_LATENT_HEAT(p, m);
real lamda = C_K_L(c, t0); //thermal conductivity of gas
real cpp = DPM_SPECIFIC_HEAT(p, T_spe);
real cpg = C_CP(c, t0); /*specific heat of gas*/
real r0 = 0.5 * P_DIAM0(p);
real h_infi = T_ambi * C_CP(c, t0);//enthaphy of gas
real h_b = cpp * Tp;//enthaphy of particle
real visg = DPM_Liquid_Property(p, m, T_spe); /*turbulent viscosity of gas*/

real B_T = 2.0, B_M, F_M, F_P;//spalding_heat_transfer_number;
real epsilon = DPM_EMISSIVITY(p, m);
real sigma = 5.670e-8;//STEFAN-BOLTZMANN_CONSTANT

real vel = sqrt(pow((P_VEL(p) - C_U(c, t0)), 2.0) + pow((P_VEL(p) - C_V(c, t0)), 2.0) + pow((P_VEL(p) - C_W(c, t0)), 2.0));
real Re = ((vel * P_DIAM0(p) * C_R(c, t0)) / visg);//Re num
real Pr = (cpp * visg) / lamda;//Prandtl num Pr = (c0->sHeat)*(c0->mu)/(c0->tCond);
real Sc = visg / (P_RHO(p) * DPM_DIFFUSION_COEFF(p, T_spe));//Sc = (c0->mu)/((c0->rho)*(diffusion_coeff));
real Nu, Sh;
real mp_dot;

/*compute m_flash*/
if (0 <= (P_T(p) - Tb) && (P_T(p) - Tb) <= 5) {
alphaf = 760 * pow((P_T(p) - Tb), 0.26);
}
else if ((P_T(p) - Tb) > 5 && (P_T(p) - Tb) < 25) {
alphaf = 27 * pow((P_T(p) - Tb), 2.33);
}
else if ((P_T(p) - Tb) >= 25) {
alphaf = 13800 * pow((P_T(p) - Tb), 0.39);
}

m_flash = (alphaf * Ap * (P_T(p) - Tb)) / LTb;

/* compute surface fuel mass fraction */
real Y_fuel_surf = C_DPMS_SURF_YI(c, t0, i);// P_sat/C_P(c,t0)
real Y_fuel_ambi = C_YI(c, t0, i);
/* compute initial values of Sherwood and Nusselt numbers using Frossling correlations */
Sh = 2.0 + 0.552 * pow(Re, 0.5) * pow(Sc, 1.0 / 3.0);
Nu = 2.0 + 0.552 * pow(Re, 0.5) * pow(Pr, 1.0 / 3.0);

/* calculate B_M and F_M */
B_M = (Y_fuel_surf - Y_fuel_ambi) / (1.0 - Y_fuel_surf);
F_M = pow(1.0 + B_M, 0.7) * log(1.0 + B_M) / B_M;

real Le = (lamda / (rho * cpp * DPM_DIFFUSION_COEFF(p, T_spe)));

/*spalding heat transfer number*/
real F_T, B_T_prev = 0.0, phi;
while (fabs(B_T - B_T_prev) > 1e-6) {
/* update the previous value of B_T */
B_T_prev = B_T;

/* calculate the current values of Sh and Nu */
F_T = pow(1.0 + B_T, 0.7) * log(1.0 + B_T) / B_T;
Sh = 2.0 + (Sh - 2.0) / F_M;
Nu = 2.0 + (Nu - 2.0) / F_T;

/* calculate phi and update B_T */
phi = (cpp / cpg) * (Sh / Nu) / Le;
B_T = pow(1.0 + B_M, phi) - 1.0;
}

real m_heat_old = 0.0;
real remaining_mass = 0.0;

BoilingLaw(p);
real dt = P_DT(p);
while (fabs(m_heat - m_heat_old) > 1e-4) {
m_heat_old = (P_MASS0(p) - P_MASS(p)) / dt;
m_heat = 2 * M_PI * lamda / cpp * r0 * Nu / (1 + m_flash / m_heat_old)
* log(1 + (1 + m_flash / m_heat_old) * (h_infi - h_b) / LTb);
m_heat_old = m_heat;
}

phi_radi = epsilon * sigma * Ap * (pow(T_ambi, 4.0) - pow(Tb, 4.0));

mp_dot = m_heat + m_flash + m_radi;
if (mp_dot < 0.)
mp_dot = 0.0;

UpdateTemperature(p, -mp_dot * LTb, P_MASS(p), 1.);
remaining_mass = P_MASS(p) - mp_dot * dt;
/* if the particle is very small, or the particle temperature
exceeds the boiling temperature dump the remaining mass */
if ((remaining_mass < TP_DPM_MINIMUM_MASS(p)) || (P_DIAM(p) <= DPM_LOWEST_DIAMETER))
{
P_DIAM(p) = 0.0;
P_T(p) = Tp;
P_MASS(p) = 0.;
}
else
{
P_MASS(p) -= mp_dot * dt;
P_DIAM(p) = DPM_DIAM_FROM_VOL(P_MASS(p) / P_RHO(p));
P_T(p) -= LTb / P_MASS(p) * DPM_SPECIFIC_HEAT(p, P_T(p));
}
}

• Rob
Ansys Employee

Other than recommending indents and comments to identify & record what each section is for, no. If you are seeing a compilation failure read the error carefully as they're fairly useful.

• Haotian Wu
Subscriber

The problem is not during compiling, it was after the calculation started. The only message was:

RANK 5 PID 47124
EXIT STATUS: -1

RANK 5 PID 112408
EXIT STATUS: -1

After these messages, the fluent was terminated.

• Rob
Ansys Employee

Debug on a smaller model. Then turn bits off until it works. Turn bits on until it doesn't. That's pretty much how I debug UDFs in the event I can't get a colleague to write them for me.

• Haotian Wu
Subscriber

Thanks for the tips. 