CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > ANSYS > FLUENT > Fluent UDF and Scheme Programming

A UDF of DEFINE_VR_RATE

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By vinerm

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 17, 2020, 22:25
Unhappy A UDF of DEFINE_VR_RATE
  #1
WZX
New Member
 
zaixing
Join Date: Apr 2020
Posts: 9
Rep Power: 2
WZX is on a distinguished road
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

Last edited by WZX; May 18, 2020 at 05:12.
WZX is offline   Reply With Quote

Old   May 17, 2020, 22:34
Default For the convenience of reading, the complicated formula is simplified
  #2
WZX
New Member
 
zaixing
Join Date: Apr 2020
Posts: 9
Rep Power: 2
WZX is on a distinguished road
#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) ));
}
WZX is offline   Reply With Quote

Old   May 18, 2020, 04:05
Default Partial Pressure
  #3
Senior Member
 
vinerm's Avatar
 
Vinerm
Join Date: Jun 2009
Location: Nederland
Posts: 2,546
Blog Entries: 1
Rep Power: 28
vinerm will become famous soon enough
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.
__________________
Regards,
Vinerm

PM to be used if and only if you do not want something to be shared on the Forum
vinerm is offline   Reply With Quote

Old   May 18, 2020, 05:11
Default
  #4
WZX
New Member
 
zaixing
Join Date: Apr 2020
Posts: 9
Rep Power: 2
WZX is on a distinguished road
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
WZX is offline   Reply With Quote

Old   May 18, 2020, 05:29
Default Ideal Gas Law
  #5
Senior Member
 
vinerm's Avatar
 
Vinerm
Join Date: Jun 2009
Location: Nederland
Posts: 2,546
Blog Entries: 1
Rep Power: 28
vinerm will become famous soon enough
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.
__________________
Regards,
Vinerm

PM to be used if and only if you do not want something to be shared on the Forum
vinerm is offline   Reply With Quote

Old   May 19, 2020, 10:19
Default
  #6
New Member
 
Vidushi
Join Date: Sep 2018
Location: Cape Town
Posts: 20
Rep Power: 3
Vidushi is on a distinguished road
Quote:
Originally Posted by WZX View Post
#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?
Vidushi is offline   Reply With Quote

Old   May 19, 2020, 10:46
Default
  #7
WZX
New Member
 
zaixing
Join Date: Apr 2020
Posts: 9
Rep Power: 2
WZX is on a distinguished road
It depends on your FLUENT /materials specieces. Do you also need to use chemical reaction UDF code? We can discuss it together.
WZX is offline   Reply With Quote

Old   May 19, 2020, 13:49
Default
  #8
New Member
 
Vidushi
Join Date: Sep 2018
Location: Cape Town
Posts: 20
Rep Power: 3
Vidushi is on a distinguished road
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 is offline   Reply With Quote

Old   May 19, 2020, 14:06
Default
  #9
New Member
 
Vidushi
Join Date: Sep 2018
Location: Cape Town
Posts: 20
Rep Power: 3
Vidushi is on a distinguished road
Quote:
Originally Posted by WZX View Post
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
....
Attached Images
File Type: png pic_1-materials panel.PNG (68.3 KB, 2 views)
File Type: png Pic_2-Species Panel in Mixture Properties.PNG (25.1 KB, 5 views)
Vidushi is offline   Reply With Quote

Old   May 19, 2020, 20:56
Default
  #10
WZX
New Member
 
zaixing
Join Date: Apr 2020
Posts: 9
Rep Power: 2
WZX is on a distinguished road
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.
WZX is offline   Reply With Quote

Old   May 20, 2020, 07:19
Default
  #11
New Member
 
Vidushi
Join Date: Sep 2018
Location: Cape Town
Posts: 20
Rep Power: 3
Vidushi is on a distinguished road
Quote:
Originally Posted by WZX View Post
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.
Vidushi is offline   Reply With Quote

Old   May 20, 2020, 07:24
Default C_P
  #12
Senior Member
 
vinerm's Avatar
 
Vinerm
Join Date: Jun 2009
Location: Nederland
Posts: 2,546
Blog Entries: 1
Rep Power: 28
vinerm will become famous soon enough
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 likes this.
__________________
Regards,
Vinerm

PM to be used if and only if you do not want something to be shared on the Forum
vinerm is offline   Reply With Quote

Old   May 20, 2020, 08:41
Default
  #13
New Member
 
Vidushi
Join Date: Sep 2018
Location: Cape Town
Posts: 20
Rep Power: 3
Vidushi is on a distinguished road
Quote:
Originally Posted by vinerm View Post
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.
Vidushi is offline   Reply With Quote

Old   May 20, 2020, 08:47
Default Backward Reaction
  #14
Senior Member
 
vinerm's Avatar
 
Vinerm
Join Date: Jun 2009
Location: Nederland
Posts: 2,546
Blog Entries: 1
Rep Power: 28
vinerm will become famous soon enough
Just attach rev_ to the variable, such as, rev_A, rev_b, rev_E, etc.
__________________
Regards,
Vinerm

PM to be used if and only if you do not want something to be shared on the Forum
vinerm is offline   Reply With Quote

Old   May 20, 2020, 09:18
Default
  #15
WZX
New Member
 
zaixing
Join Date: Apr 2020
Posts: 9
Rep Power: 2
WZX is on a distinguished road
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?
WZX is offline   Reply With Quote

Old   May 20, 2020, 09:43
Default Partial Pressure
  #16
Senior Member
 
vinerm's Avatar
 
Vinerm
Join Date: Jun 2009
Location: Nederland
Posts: 2,546
Blog Entries: 1
Rep Power: 28
vinerm will become famous soon enough
Partial Pressure is never a function of total pressure but of static absolute pressure. So, you should multiply it by C_P + Operating Pressure.
__________________
Regards,
Vinerm

PM to be used if and only if you do not want something to be shared on the Forum
vinerm is offline   Reply With Quote

Old   May 20, 2020, 14:54
Default
  #17
New Member
 
Vidushi
Join Date: Sep 2018
Location: Cape Town
Posts: 20
Rep Power: 3
Vidushi is on a distinguished road
Quote:
Originally Posted by vinerm View Post
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?
Vidushi is offline   Reply With Quote

Old   May 20, 2020, 15:00
Default Backward Reaction
  #18
Senior Member
 
vinerm's Avatar
 
Vinerm
Join Date: Jun 2009
Location: Nederland
Posts: 2,546
Blog Entries: 1
Rep Power: 28
vinerm will become famous soon enough
Yes, that's correct.
__________________
Regards,
Vinerm

PM to be used if and only if you do not want something to be shared on the Forum
vinerm is offline   Reply With Quote

Old   May 20, 2020, 22:39
Default
  #19
WZX
New Member
 
zaixing
Join Date: Apr 2020
Posts: 9
Rep Power: 2
WZX is on a distinguished road
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).
WZX is offline   Reply With Quote

Old   May 25, 2020, 10:39
Default Temperature
  #20
Senior Member
 
vinerm's Avatar
 
Vinerm
Join Date: Jun 2009
Location: Nederland
Posts: 2,546
Blog Entries: 1
Rep Power: 28
vinerm will become famous soon enough
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.
__________________
Regards,
Vinerm

PM to be used if and only if you do not want something to be shared on the Forum
vinerm is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
udf for one dimensional linear motion based on force maccheese Fluent UDF and Scheme Programming 2 September 1, 2019 02:18
Save output of udf in another udf! JuanJoMex FLUENT 0 February 8, 2018 12:43
UDF Compilation Error - Loading Library - COMMON Problem! Help! robtheslob Fluent UDF and Scheme Programming 8 July 24, 2015 00:53
UDF parallel error: chip-exec: function not found????? shankara.2 Fluent UDF and Scheme Programming 1 January 16, 2012 22:14
UDF, UDF, UDF, UDF Luc SEMINEL Main CFD Forum 0 November 25, 2002 04:01


All times are GMT -4. The time now is 06:37.