[vismayie]

hello,

I am simulating a tubular reactor in Fluent with reaction as the wall boundary condition. The simulations didnt show any problem till temperature 775 K but suddenly at that temperature there is something like the reversed flow at the pressure outlet. and the outlet concentration of the reactant which is expected to decrease is increasing.

When i ask Fluent to compute the weighted average of the rate of reaction at the wall it shows the following : -1.#IND

So to my best knowledge something is wrong with the reaction rate at this temperature. Is it because of a very fine mesh or something ? I am using a 2D mesh.

the following is my UDF:


#include "udf.h"
#define R 8.314
#define p 101325
#define pi 3.14159


DEFINE_SR_RATE(my_rate,f,t,r,mw,yi,rr)
{
Thread *t0 = t -> t0;
cell_t c0=F_C0(f,t);


real n2 = yi[0];
real n2o = yi[1];
real co = yi[2];
real co2 = yi[3];
real he = yi[4];

real rho_w = 1.0, T_w, sd=2.68*1e-5;
real xno = 1.0, xco = 1.0, xn2o = 1.0, xn2 = 1.0, xco2 = 1.0;

real pco2, pco, pn2o, pn2;
real k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13, k14, k15, k16;
real K1, K2, K3, K4, K5, K6, K7, K8;
real zor3, kor3, denom4new, denom5new, denom8new, R4, R5, R8, rate;


T_w = F_T(f,t);
rho_w = C_R(c0,t0)*C_T(c0,t0)/T_w;


co *= rho_w/mw[2];
n2o *= rho_w/mw[1];
co2 *= rho_w/mw[3];
n2 *= rho_w/mw[0];
he *= rho_w/mw[4];



xco = co/(co+n2o+co2+n2+he);
xn2o = n2o/(co+n2o+co2+n2+he);
xn2 = n2/(co+n2o+co2+n2+he);
xco2 = co2/(co+n2o+co2+n2+he);

pn2 = xn2*p;
pco = xco*p;
pco2 = xco2*p;
pn2o= xn2o*p;



k1 = 0.67*(sqrt(1000/(R*T_w*2*pi*mw[3])));
k2 = (1e13*sd)*exp(-26*4185.8/(R*T_w));
k3 = (1e11*sd)*exp(-6.7*4185.8/(R*T_w));
k4 = (1e11*sd)*exp(-47.6*4185.8/(R*T_w));
k5 = (1e11*sd)*exp(-27*4185.8/(R*T_w));
k6 = (0.001)*(sqrt(1000/(R*T_w*2*pi*mw[1])));
k7 = (1e14*sd)*exp(-24.3*4185.8/(R*T_w));
k8 = (1e14*sd)*exp(-2.1*4185.8/(R*T_w));
k9 = (1e13*sd)*exp(-10*4185.8/(R*T_w));
k10 = 0.001*(sqrt(1000/(R*T_w*2*pi*mw[2])));
k11 = 0.67*(sqrt(1000/(R*T_w*2*pi*mw[3])));
k12 = (1e13*sd)*exp(-32*4185.8/(R*T_w));
k13 = (1e11*sd)*exp(-24.35*4185.8/(R*T_w));
k14 = (1e11*sd)*exp(-11.3*4185.8/(R*T_w));
k15 = (1e13*sd)*exp(-5.17*4185.8/(R*T_w));
k16 = 0.005*(sqrt(1000/(R*T_w*2*pi*mw[4])));

K1 = k1/k2;
K2 = k3/k4;
K3 = k5/k6;
K4 = k7/k8;
K5 = k9/k10;
K6 = k11/k12;
K7 = k13/k14;
K8 = k15/k16;

denom4new=1+((pow((pn2/K3), 0.5))+(K6*pco))+((((pow((pn2/K3), 0.5))*(pco2/(K7*K8*K6*pco)))/K2));
R4=(pow(denom4new,2))/(k8*(pn2o/K5));

denom5new=1+((pow((pn2/K3), 0.5))+(K6*pco))+(((pow((pn2/K3), 0.5))*(pco2/(K7*K8*K6*pco)))/K2);
R5=(denom5new)/(k10*(pn2o));


denom8new=1+((pow((pn2/K3), 0.5))+(K6*pco))+((pn2o/K5)/(K4*(pow((pn2/K3), 0.5))));
R8=(denom8new)/((k15)*(((K6*pco)/(K7))/(((K2*pn2o)/K5)/(K4*((pn2/K3))))));

kor3=((K2*K3*K6*K7*K8)/(K4*K5));

zor3=(((pco2*pn2))/(kor3*(pn2o*pco)));

rate = (((1-zor3)))/(R4+R5+R8);

if (STREQ(r->name, "reaction-1"))
{
*rr = 0.0013*rate;
C_UDMI(c0,t0,0) = *rr;
}


}


Please suggest me a way out. I have my project submission very soon and i am not able to figure out anything.

Thanking you,

Vismayie
IIT Madras