Reply To: UDF for surface reaction rate〈Please help a beginner〉

[Yuya.S]

Here is the modified macro. 

#include "udf.h"
DEFINE_SR_RATE(my_rate,f,t,r,mw,yi,rr)
{
 Thread *t0=THREAD_T0(t);
 cell_t c0=F_C0(f,t);
/*mass fraction of species i at the wall*/
real y_h2=yi[0];
real y_co2=yi[1];
real y_ch4=yi[2];
real y_h2o=yi[3];
real Nsum, R, A1, A2, E1, E2, k1, k2, Keq, T_w, TP, r1;

/*calculate species i in the unit of kmol i/kg mix*/
 y_h2 *= 1/mw[0];
 y_co2 *= 1/mw[1];
 y_ch4 *= 1/mw[2];
 y_h2o *= 1/mw[3];
/*total mole number per kg mix */
 Nsum = y_h2o + y_ch4 + y_h2 + y_co2;

/*calculate mole fraction of species i in the unit of kgmol i/kgmol mix*/
 y_h2o *= 1/Nsum;
 y_ch4 *= 1/Nsum;
 y_h2 *= 1/Nsum;
 y_co2 *= 1/Nsum;
/*gas constant, J/molK */
 R=8.314;
/*equilibrium constant calculated by the empirical formula*/
 T_w=F_T(f,t);

Keq=137*pow(T_w,-3.998)/exp(158.7e+3/R/T_w);
/*reaction rate constant, need to be guessed, trial-and-error*/
 A1=3.47e+4, A2=6.80e+7;
 E1=56.0e+3, E2=89.0e+3;          /*J/mol*/
 k1=A1*exp(-E1/R/T_w);
 k2=A2*exp(-E2/R/T_w);

/*total pressure in the cell near wall, bar*/
 TP=C_P(c0,t0)/1.0e+5;

r1 = (k1*pow(TP*y_h2, 0.49)*pow(TP*y_co2, 0.42) / (1 + k2*TP*y_h2o*pow(TP*y_h2, -0.5)))*(1 - (TP*y_ch4*pow(TP*y_h2o, 2)) / pow(TP*y_h2, 4)*TP*y_co2 * Keq);

if(FLUID_THREAD_P(t) && THREAD_VAR(t).fluid.porous)

 *rr = r1;  
}

At least the error is resolved and the computation itself works correctly after hooking the UDF. However, the key products, methane and water vapor, are not produced at all. I would be happy to take advice from anyone.