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.