[WZX]

Hello, friends！

I use CFD simulation to study the hydrogen production by Methanol-Steam Reforming Reaction (MSRR). The overall reactions of MSRR can be written as follows:

The reaction rate expressions for these overall reactions involved in the
process of methanol-steam reforming are as follows：

Parameters of the comprehensive kinetic model of the methanol-steam reforming reaction process

And these three reaction rates can be implemented by a UDF of DEFINE_VR_RATE(MSR_rate,c,t,r,mw,yi,rr,rr_t ). However, when I was editing the code, I encountered the following difficulties:
a)In the formula, the partial pressure of each gas component is used, but in DEFINE_VR_RATE, there is only the mass fraction yi, so I want to get the mole fractionXi through the mass fraction, and then multiply it by the total pressure to get the partial pressure of each component.Is this thinking right?And how to express the total pressure?C_P(C,t)？
b) In the formula, the denominator contains the partial pressure of products such as H2 and CO. however, my import condition gives the mass fraction of reactant CH3OH & H2, which leads to the fact that the denominator is 0 at the beginning and cannot operate. How can I solve this problem?
Can anyone help me

[WZX]

#include "udf.h"

#define C_S1_T 7.5e-6
#define C_S2_T 7.5e-6
#define C_S1a_T 1.5e-5
#define C_S2a_T 1.5e-5
#define S_c 5.0e+4
#define rou_b 1.3e+3

#define R_U 8.314


DEFINE_VR_RATE(MSR_rate,c,t,r,mw,yi,rr,rr_t )
{
real Y_CH3OH, Y_H2O, Y_H2, Y_CO, Y_CO2; /* Define the mass fraction of each component of the gas mixture*/ *
real X_CH3OH, X_H2O, X_H2, X_CO, X_CO2; /* mole fractaion */
real mw_CH3OH, mw_H2O, mw_H2, mw_CO, mw_CO2;
real MW_mix; /* the average molecular mass of the gas mixture*/

real kR, kD, kW;
real K_R, K_D, K_W;

real K_CH3O_1, K_CH3O_2, K_HCOO_1, K_OH_1, K_OH_2, K_H_1a, K_H_2a;

real r_MSR,r_MD,r_rWGS;

real temp= C_T(c,t);

Y_CH3OH =yi[0];
Y_H2O = yi[1];
Y_H2 = yi[2];
Y_CO = yi[3];
Y_CO2 = yi[4];

mw_CH3OH = mw[0];
mw_H2O = mw[1];
mw_H2 = mw[2];
mw_CO = mw[3];
mw_CO2 = mw[4];


/* Partial Pressure is equal to mole fraction when operating pressure is 1 atm */

MW_mix=1.0 / ((Y_CH3OH/mw_CH3OH)+(Y_H2O/mw_H2O)+(Y_CO/mw_CO)+(Y_H2/mw_H2)+(Y_CO2/mw_CO2));

/* Is that ok? or I still have to multiply the mole fraction by the total pressure*/
X_CH3OH = Y_CH3OH*MW_mix/mw_CH3OH;
X_H2O = Y_H2O*MW_mix/mw_H2O;
X_H2 = Y_H2*MW_mix/mw_H2;
X_CO2 = Y_CO2*MW_mix/mw_CO2;

/* Equilibrium Constants*/

K_CH3O_1 = exp(-41.8/R_U+20.0/(R_U*temp));

K_HCOO_1 = exp(179.2/R_U-10.0/(R_U*temp));

K_OH_1 = exp(-44.5/R_U+20.0/(R_U*temp));

K_H_1a = exp(-100.8/R_U+50.0/(R_U*temp));

/*Kinetics Constants */
kR = (7.4e+14)*exp(-102.8/(R_U*temp));

K_R = pow(10.0, 10.592-3066.0/temp);


if (!strcmp(r->name, "reaction-1"))

{
r_MSR = C_S1_T*C_S1a_T *S_c*rou_b *(kR *K_CH3O_1* (X_CH3OH /(pow(X_H2,0.5)))*(1-(X_CO2*pow(X_H2,3.0)/(K_R*X_CH3OH*X_H2O))))/((1+K_CH3O_1* (X_CH3OH /(pow(X_H2,0.5)))+K_HCOO_1* X_CO2 *(pow(X_H2,0.5))+K_OH_1*(X_H2O /(pow(X_H2,0.5))))*(1+pow(K_H_1a,0.5)*pow(X_H2,0.5) ));
}

[vinerm]

1. If you are using Ideal-gas law for mixture density, then it is straightforward. You can determine pressure of each gas using ideal-gas law. Number of moles or mass can be derived from mass fraction, volume, and density in the cell. Molecular weight is, of course, known. Do note that for last species, you should determine pressure from Dalton's law else the total pressure will not match C_P. Furthermore, C_P returns gauge pressure while ideal gas law is based on absolute pressure. So, you must ensure that the operating pressure is set as 0 in your case.

2. There are two ways to deal with it. One is to check for the partial pressure and use the equation only when partial pressure is above a certain threshold. However, if the equation is valid even for very low partial pressure values, i.e., it has realizable asymptote, then you can just add a very small value to the denominator. This will ensure that there is no division by 0.

[WZX]

Thank you for your reply.
1. I use the pressure-based model and the incompressible Ideal-gas law for mixture density in fluent.Do I need to turn incompressible into Ideal-gas law for mixture density？
2. According to the FLUENT USER'S GUIDE, the last species in the Selected Species list is the most abundant species. However, I can't determine which gas speciece is the most after the final reaction. what should I do?
3. I choose the inlet as velocity-inlet, then I use pressure-outlet boundary conditions at the outlet. The gauage pressure is set as 0 in pressure-outlet B.C and 101325 Pa for operating pressure. You said I must ensure that the operating pressure is set as 0 .If the operating pressure is set as 0, should I set the gauage pressure as 101325 pa?

__________________
Regards,
WZX

[vinerm]

If Mach number is less than 0.2, then incompressible ideal gas law is good enough. If it is higher, prefer using Ideal Gas Law.

Operating pressure need to be set to 0 only if you use ideal gas law, not otherwise. Yes, if it is set to 0, then outlet pressure must be set to 101325 Pa.

Setting most abundant as last specie is to ensure that round-off error does not cause significant error in mass fractions of species. It does not have to be really the most abundant one but the one that is abundant enough. If a specie has a low mass fraction, say less than 0.02 then keeping it last could lead to wrong prediction. But if you have two species that have higher than 0.2 or 0.3 mass fraction, then any one of those can be used as last one.

[Vidushi]

Originally Posted by WZX

#include "udf.h"

#define C_S1_T 7.5e-6
#define C_S2_T 7.5e-6
#define C_S1a_T 1.5e-5
#define C_S2a_T 1.5e-5
#define S_c 5.0e+4
#define rou_b 1.3e+3

#define R_U 8.314


DEFINE_VR_RATE(MSR_rate,c,t,r,mw,yi,rr,rr_t )
{
real Y_CH3OH, Y_H2O, Y_H2, Y_CO, Y_CO2; /* Define the mass fraction of each component of the gas mixture*/ *
real X_CH3OH, X_H2O, X_H2, X_CO, X_CO2; /* mole fractaion */
real mw_CH3OH, mw_H2O, mw_H2, mw_CO, mw_CO2;
real MW_mix; /* the average molecular mass of the gas mixture*/

real kR, kD, kW;
real K_R, K_D, K_W;

real K_CH3O_1, K_CH3O_2, K_HCOO_1, K_OH_1, K_OH_2, K_H_1a, K_H_2a;

real r_MSR,r_MD,r_rWGS;

real temp= C_T(c,t);

Y_CH3OH =yi[0];
Y_H2O = yi[1];
Y_H2 = yi[2];
Y_CO = yi[3];
Y_CO2 = yi[4];

mw_CH3OH = mw[0];
mw_H2O = mw[1];
mw_H2 = mw[2];
mw_CO = mw[3];
mw_CO2 = mw[4];


/* Partial Pressure is equal to mole fraction when operating pressure is 1 atm */

MW_mix=1.0 / ((Y_CH3OH/mw_CH3OH)+(Y_H2O/mw_H2O)+(Y_CO/mw_CO)+(Y_H2/mw_H2)+(Y_CO2/mw_CO2));

/* Is that ok? or I still have to multiply the mole fraction by the total pressure*/
X_CH3OH = Y_CH3OH*MW_mix/mw_CH3OH;
X_H2O = Y_H2O*MW_mix/mw_H2O;
X_H2 = Y_H2*MW_mix/mw_H2;
X_CO2 = Y_CO2*MW_mix/mw_CO2;

/* Equilibrium Constants*/

K_CH3O_1 = exp(-41.8/R_U+20.0/(R_U*temp));

K_HCOO_1 = exp(179.2/R_U-10.0/(R_U*temp));

K_OH_1 = exp(-44.5/R_U+20.0/(R_U*temp));

K_H_1a = exp(-100.8/R_U+50.0/(R_U*temp));

/*Kinetics Constants */
kR = (7.4e+14)*exp(-102.8/(R_U*temp));

K_R = pow(10.0, 10.592-3066.0/temp);


if (!strcmp(r->name, "reaction-1"))

{
r_MSR = C_S1_T*C_S1a_T *S_c*rou_b *(kR *K_CH3O_1* (X_CH3OH /(pow(X_H2,0.5)))*(1-(X_CO2*pow(X_H2,3.0)/(K_R*X_CH3OH*X_H2O))))/((1+K_CH3O_1* (X_CH3OH /(pow(X_H2,0.5)))+K_HCOO_1* X_CO2 *(pow(X_H2,0.5))+K_OH_1*(X_H2O /(pow(X_H2,0.5))))*(1+pow(K_H_1a,0.5)*pow(X_H2,0.5) ));
}

Sir, how did you know the order of species mass fractions stored in yi?
I mean how did you figure out that yi[0] is storing the mass fraction of CH3OH and not the mass fraction of any other specie in the reaction?

What is the order yi or wi follows?

[WZX]

It depends on your FLUENT /materials specieces. Do you also need to use chemical reaction UDF code? We can discuss it together.

[Vidushi]

Yes, please. I am doing something very similar to you. I would love to discuss it with you. I am writing a DEFINE_VR_Rate for the Fischer Tropsch Polymerisation process.

[Vidushi]

Originally Posted by WZX

It depends on your FLUENT /materials specieces. Do you also need to use chemical reaction UDF code? We can discuss it together.

Do you mean it is in the order of species in the materials panel? Or is it in the order of species in the Mixture properties? I have attached the pictures for the two cases. please let me know which is the right one. If we assume it is the pic 1, i.e. the materials panel, Is the following true:

yi[0]=ethane
yi[1]=methane
yi[2]=hydrogen
yi[3]=oxygen
....

[WZX]

THE SELECTED SPECIES
yi[0]=H2O
yi[1]=N2
yi[2]=CH4
...
By the way, do you need the partial pressure of the gas component in your formula？What did you do with it? This has puzzled me for a long time.

[Vidushi]

Originally Posted by WZX

THE SELECTED SPECIES
yi[0]=H2O
yi[1]=N2
yi[2]=CH4
...
By the way, do you need the partial pressure of the gas component in your formula？What did you do with it? This has puzzled me for a long time.

I am working with activities, which is given by

activity = mole frac x fugacity x total pressure/Ref. pressure

I guess the total pressure should be the one returned by C_P(c,t) here.

[vinerm]

C_P returns static gauge pressure and not total pressure. You have to add dynamic pressure and operating pressure to it. If operating pressure is set to 0, then C_P returns static absolute pressure. So, you have to add only dynamic pressure to it.

[Vidushi]

Originally Posted by vinerm

C_P returns static gauge pressure and not total pressure. You have to add dynamic pressure and operating pressure to it. If operating pressure is set to 0, then C_P returns static absolute pressure. So, you have to add only dynamic pressure to it.

Thank you Vinerm for the information.
I have one more query regarding DEFINE VR_RATE. I request you to please help me with it if you can.

So, in DEFINE_VR_RATE, while defining the Arrhenius rate:

*rate = r->A * exp( - r->E / (UNIVERSAL_GAS_CONSTANT * C_T(c,t))) *
pow(C_T(c,t), r->b) * prod;


r->A and r->E retrieves the pre-exponential factor and the activation energy for the current reaction that the user has inputted in the reactions panel from the data structure reaction stored in Fluent.
If I enable the backward reactions and input the backward reaction pre-exponential factor and the backward activation energy, how can I retrieve their values in my UDF.

Thank you in advance.

[vinerm]

Just attach rev_ to the variable, such as, rev_A, rev_b, rev_E, etc.

[WZX]

Dear sir,
According to your suggestions, I multiply the mole fraction of each component by the total pressure to get the partial pressure of the corresponding component.Can I use the method shown in the figure below in UDF? P is the total pressure.

The following figure shows the settings in fluent, operating pressure is 0, gauage pressuure is 101235Ppa, is it right?


And for the mole fracation, I took the following method, is it right?

[vinerm]

Partial Pressure is never a function of total pressure but of static absolute pressure. So, you should multiply it by C_P + Operating Pressure.

[Vidushi]

Originally Posted by vinerm

Just attach rev_ to the variable, such as, rev_A, rev_b, rev_E, etc.

Okay. So you mean I should write it as r->rev_A, r->b etc?

[vinerm]

Yes, that's correct.

[WZX]

Hi,friend.Why the temperature of chemical reactiuon is negative? Even if I gave a constant reaction rate. I took the function C_T(f,t).

[vinerm]

What do you mean by temperature of chemical reaction? C_T requires proper arguments, i.e., of type cell_t and Thread *. If those are of wrong type, then return value would be wrong.