CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Implement code to calculate the reaction rate in reactingFoam

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

Like Tree3Likes
  • 3 Post By mturcios777

Reply
 
LinkBack Thread Tools Display Modes
Old   July 26, 2012, 09:50
Default Implement code to calculate the reaction rate in reactingFoam
  #1
New Member
 
Join Date: Jul 2012
Posts: 11
Rep Power: 5
matari is on a distinguished road
Hello,

I'm new in programming in OpenFoam.
I use the reactingFoam solver in OpenFoam-2.1.x and want to calculate the reaction rate with a special equation which is not implemented in OpenFoam and which is similar to the arrhenius equation. The arrhenius reaction rate can be calculated in OpenFoam, so I would like to modify this calculation. How can I modify the arrhenius equation and compile the new reaction type to use it in reactingFoam?

Thanks in advance for any hint,

Matari

Last edited by matari; July 27, 2012 at 05:11.
matari is offline   Reply With Quote

Old   July 29, 2012, 04:16
Default
  #2
Senior Member
 
Elvis
Join Date: Mar 2009
Location: Sindelfingen, Germany
Posts: 577
Blog Entries: 5
Rep Power: 13
elvis is on a distinguished road
Hi Matari,
there is a Howto that shows how you should alter an existing solver (icoFoam)
http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam

this is a good guideline for your reactingFoam
elvis is offline   Reply With Quote

Old   August 1, 2012, 04:04
Default
  #3
New Member
 
Join Date: Jul 2012
Posts: 11
Rep Power: 5
matari is on a distinguished road
Hello Elvis,

thanks for your reply. I've done this tutorial but I have problems to adapt it to the reaction calculation in reactingFoam. In reactingFoam I can find the reaction type in my CASE/constant/reactions file. For example:

propaneReaction
{
type irreversibleArrheniusReaction;
reaction "CH4 + 2O2 = CO2 + 2H20";
A 9.5e+11;
beta 0;
Ta 23650;
}

I found the calculation for the Arrhenius reaction rate in $FOAM_SRC/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H
In this folder are the two files ArrheniusReactionRate.H and ArrheniusReactionRateI.H.
My problems now are:
1. How can I compile changes in ArrheniusReactionRateI.H.
2. The Arrhenius equation also includes the concentration of the elements. In ArrheniusReactionRateI.H it seems that there is only defined the first part of the equation without the concentrations:
{
scalar ak = A_;

if (mag(beta_) > VSMALL)
{
ak *= pow(T, beta_);
}

if (mag(Ta_) > VSMALL)
{
ak *= exp(-Ta_/T);
}

return ak;
}
3. How can I define a new type to use in the CASE/constant/reactions file instead of irreversibleArrheniusReaction. I don't find a file where irreversibleArrheniusReaction is defined.

Does anyone have a suggestion to solve these problems?

Best regards,

Matari
matari is offline   Reply With Quote

Old   August 1, 2012, 10:21
Default
  #4
New Member
 
Join Date: Jul 2012
Posts: 11
Rep Power: 5
matari is on a distinguished road
Hi,

ok, maybe I should describe my problem a bit more precisely.
The Arrhenius reaction rate for the reaction
CH4^0.7 + 2O2^0.8 = CO2 + 2H20
has to be calculated using the formula R = A*T^b*exp(-Ta/T)*c(CH4)^0.7 * c(O2)^0.8
It seems like the first part A*T^b*exp(-Ta/T) is defined in ArrheniusReactionRateI.H
I'm not sure, but it looks like the concentrations are calculated in the Reaction.C and/or the ReactionI.H file. Can anyone tell me where I find the calculation of the second part "c(CH4)^0.7 * c(O2)^0.8"? This is the part I need to change for my calculations.
matari is offline   Reply With Quote

Old   August 1, 2012, 11:16
Default
  #5
Member
 
MSR CHANDRA MURTHY
Join Date: Mar 2009
Posts: 32
Rep Power: 8
chandramurthy is on a distinguished road
You many look at $FOAM_SRC/combustionModels folder, several models are implemented to calculated rate of production of species (w-dot) which appear in species conservation equation. You may get some idea while closely looking at various combustion models.

Quote:
Originally Posted by matari View Post
Hi,

ok, maybe I should describe my problem a bit more precisely.
The Arrhenius reaction rate for the reaction
CH4^0.7 + 2O2^0.8 = CO2 + 2H20
has to be calculated using the formula R = A*T^b*exp(-Ta/T)*c(CH4)^0.7 * c(O2)^0.8
It seems like the first part A*T^b*exp(-Ta/T) is defined in ArrheniusReactionRateI.H
I'm not sure, but it looks like the concentrations are calculated in the Reaction.C and/or the ReactionI.H file. Can anyone tell me where I find the calculation of the second part "c(CH4)^0.7 * c(O2)^0.8"? This is the part I need to change for my calculations.
chandramurthy is offline   Reply With Quote

Old   August 1, 2012, 12:54
Default
  #6
Senior Member
 
mturcios777's Avatar
 
Marco A. Turcios
Join Date: Mar 2009
Location: Vancouver, BC, Canada
Posts: 725
Rep Power: 18
mturcios777 will become famous soon enough
The easiest long term way to do this (though it requires some extra work up front) is to add a new reaction rate (call it irreversibleArrheniusReactionRate) and link it to your case/solver so that you can have access to it.

This is very similar to creating your own boundary condition. This tutorial is a little old and sparse on detail, but should give you a start.

http://physics.ucsd.edu/students/cou.../Tutorial2.pdf

Basic steps are:

1) Copy over a class that almost does what you want
2) Rename/edit the source files and compile options so that the class you have compiles under the new name
3) Make changes to your new class to do what you want

You can either compile the new class as part of the OpenFOAM library the original class or to create a library with only your new class that you can link (statically or dynamically) to your case/solver. I prefer the latter as it keeps the original installation clean and you don't break anything.
olivier78, Hanzo and neiht like this.
mturcios777 is offline   Reply With Quote

Old   August 2, 2012, 08:59
Default
  #7
New Member
 
Join Date: Jul 2012
Posts: 11
Rep Power: 5
matari is on a distinguished road
Hello,

thank you for your reply. The tutorial is really good and helps a lot!
matari is offline   Reply With Quote

Old   August 2, 2012, 11:30
Default
  #8
New Member
 
JOMAA Ghassan
Join Date: Jan 2012
Location: Paris
Posts: 5
Rep Power: 5
ghas is on a distinguished road
Hi matari,

Quote:
Originally Posted by matari View Post
Hi,

ok, maybe I should describe my problem a bit more precisely.
The Arrhenius reaction rate for the reaction
CH4^0.7 + 2O2^0.8 = CO2 + 2H20
has to be calculated using the formula R = A*T^b*exp(-Ta/T)*c(CH4)^0.7 * c(O2)^0.8
It seems like the first part A*T^b*exp(-Ta/T) is defined in ArrheniusReactionRateI.H
I'm not sure, but it looks like the concentrations are calculated in the Reaction.C and/or the ReactionI.H file. Can anyone tell me where I find the calculation of the second part "c(CH4)^0.7 * c(O2)^0.8"? This is the part I need to change for my calculations.
You have to see $FOAM_SRC/thermophysicalModels/chemisrtymodel. The calculation of all reaction rates is done there using the complete formula of reaction rate ( R = A*T^b*exp(-Ta/T)*c(CH4)^0.7 in your case). The reactionRate classes are used just to calculate the kinetic constant at given temperature for different kinetic models. In Reactions classes (in $FOAM_SRC/thermophysicalModels/specie/reaction/reaction ) the other kinetic parameters are stored such as the reaction equations (CH4 + 2O2 = CO2 + 2H20 in your case), their stoichiometric coefficients, their partial kinetic orders (0.7 and 0.8 in your case) ... In the chemisrtymodel classes, the former classes are used to calculate the reaction rate( R = A*T^b*exp(-Ta/T)*c(CH4)^0.7 * c(O2)^0.8) which done in omega(...).Also, depend on which chemistry solver you use, there are many other calculation will be done.

Best regards


Ghassan
ghas is offline   Reply With Quote

Old   August 2, 2012, 12:29
Default
  #9
Senior Member
 
mturcios777's Avatar
 
Marco A. Turcios
Join Date: Mar 2009
Location: Vancouver, BC, Canada
Posts: 725
Rep Power: 18
mturcios777 will become famous soon enough
Hey Matari,

The forum software says you posted another message, and I read it in my email but don't see it on the thread. Did you sort out the issue with building your new library? From the error messages you posted it seems like you have a typo in the EXE_LIBS entry.

Let me know if you fixed it!
mturcios777 is offline   Reply With Quote

Old   August 3, 2012, 09:19
Default
  #10
New Member
 
Join Date: Jul 2012
Posts: 11
Rep Power: 5
matari is on a distinguished road
Hi,

@mturcios777: Thank you for the reply. I fixed the problem with the library a few minutes after I posted the message. There were some files missing.

In reactingFoam I use the ODEChemistryModel. From Ghassans message I understand that the reactionRate class calculates the first part of the equation. In my example: ak = A*T^b*exp(-Ta/T).
The chemistryModel class then calculates the ReactionRate R = ak * c(CH4)^0.7 * c(O2)^0.8.
My goal is to calculate the reactionRate with the equation R = rA * c(CH4) * c(O2) where rA depends on the temperature, pressure and concentration of CH4, O2, CO2 and H2O. In the ArrheniusReactionRate class ak only depends on T. So I modified the LangmuirHinshelwoodReactionRate instead and called it myLangmuirHinshelwoodReactionRate.
Here is the myLangmuirHinshelwoodReactionRateI.H file:

// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

double rA;

inline Foam::myLangmuirHinshelwoodReactionRate::myLangmui rHinshelwoodReactionRate
(
const label ch4,
const label o2,
const label co2,
const label h2o
)
:
ch4_(ch4),
o2_(o2),
co2_(co2),
h2o_(h2o)
{}
inline Foam::myLangmuirHinshelwoodReactionRate::myLangmui rHinshelwoodReactionRate
(
const speciesTable& st,
Istream& is
)
:
ch4_(st["CH4"]),
o2_(st["O2"]),
co2_(st["CO2"]),
h2o_(st["H2O"])
{}
inline Foam::myLangmuirHinshelwoodReactionRate::myLangmui rHinshelwoodReactionRate
(
const speciesTable& st,
const dictionary& dict
)
:
ch4_(st["CH4"]),
o2_(st["O2"]),
co2_(st["CO2"]),
h2o_(st["H2O"])
{}

// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //

inline Foam::scalar Foam::myLangmuirHinshelwoodReactionRate:perator( )
(
const scalar T,
const scalar p,
const scalarField& c
) const
{
rA = ... ; //calculation of rA that depends on T, p, c[ch4_], c[co2_], c[h2o_], c[h2o_]
return rA;
}
inline void Foam::myLangmuirHinshelwoodReactionRate::write(Ost ream& os) const
{}


inline Foam::Ostream& Foam:perator<<
(
Ostream& os,
const myLangmuirHinshelwoodReactionRate& lhrr
)
{
return os;
}

The compilation is successful but the calculated reaction rate is too little compared to an ANSYS Fluent simulation with an UDF that calculates the reaction rate using the equation R = rA * c(CH4) * c(O2). Is the reaction rate in my case with ODEChemistryModel really calculated using the equation R = rA * c(CH4) * c(O2) or did I misunderstood that?
matari is offline   Reply With Quote

Old   August 3, 2012, 12:18
Default
  #11
Senior Member
 
mturcios777's Avatar
 
Marco A. Turcios
Join Date: Mar 2009
Location: Vancouver, BC, Canada
Posts: 725
Rep Power: 18
mturcios777 will become famous soon enough
You have two options:

1) Compile a debug version of OF and try trace to make sure that your function is being called as it should be

2) Add your own debug messages to make sure the proper reaction rate is being calculated.

A related question: how are you linking your own code to the solver?
mturcios777 is offline   Reply With Quote

Old   August 6, 2012, 05:40
Default
  #12
New Member
 
Join Date: Jul 2012
Posts: 11
Rep Power: 5
matari is on a distinguished road
Hello Marco,

in my Make directory I generated to files:
1. files:

makeLangmuirHinshelwoodReactions.C

LIB = $(FOAM_USER_LIBBIN)/libMyLangmuirHinshelwood


2. options:

EXE_INC = \
-I$(SIB_SRC)/thermophysicalModels/specie/lnInclude

EXE_LIBS = \
-1thermophysicalModels/specie


To link my code to the solver I add libs ("libMyLangmuirHinshelwood.so") to my system/controlDict.

Best regards,

Matari
matari is offline   Reply With Quote

Old   August 7, 2012, 14:43
Default
  #13
Senior Member
 
mturcios777's Avatar
 
Marco A. Turcios
Join Date: Mar 2009
Location: Vancouver, BC, Canada
Posts: 725
Rep Power: 18
mturcios777 will become famous soon enough
Sorry its taken a while to respond, I was on vacation. Sounds like everything is linked and working well. From this point in, you will need to monitor the values that are passed to the operator() function and those that the function returns. It may be that your function is doing exactly what you intend it to do based on the inputs it receives, but what OF does with that after is no doubt different than what Fluent does, and may need adjustments down the line.

Good luck!
mturcios777 is offline   Reply With Quote

Old   August 13, 2012, 04:24
Default
  #14
New Member
 
Join Date: Jul 2012
Posts: 11
Rep Power: 5
matari is on a distinguished road
Hi,

the operator() function does what I intend it to do, so I would like to monitor the omega and dQ() that are calculated in the ODEChemistryModel.C file. Therefore I renamed ODEChemistryModel* as myChemistryModel*. I renamed all necessary files and replaced in each file "ODE" with "my". The files in my Make directory is:

makeChemistrySolvers.C
LIB = $(FOAM_USER_LIBBIN)/libmyChemistryModel


and options:

EXE_INC = \
-I$(SIB_SRC)/thermophysicalModels/chemistryModel/lnInclude

EXE_LIBS = \
-1thermophysicalModels/chemistryModel


I don't get any errors during the compilation:

options:4:12: warning: backslash-newline at end of file [enabled by default]
wmakeLnInclude: linking include files to ./lnInclude
Making dependency list for source file makeChemistrySolvers.C
SOURCE=makeChemistrySolvers.C ; g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -O3 -DNoRepository -ftemplate-depth-100 -I/thermophysicalModels/chemistryModel/lnInclude -IlnInclude -I. -I/home/fluent/OpenFOAM/OpenFOAM-2.1.x/src/OpenFOAM/lnInclude -I/home/fluent/OpenFOAM/OpenFOAM-2.1.x/src/OSspecific/POSIX/lnInclude -fPIC -c $SOURCE -o Make/linux64Gcc46DPOpt/makeChemistrySolvers.o
'/home/fluent/OpenFOAM/fluent-2.1.x/platforms/linux64Gcc46DPOpt/lib/libmyChemistryModel.so' is up to date.


I included libmyChemistryModel.so in my system/controlDict. In the constant/chemistryProperties file I have:

psiChemistryModel myChemistryModel<gasThermoPhysics>;

chemistry on;

chemistrySolver ode;

initialChemicalTimeStep 1e-07;

sequentialCoeffs
{
cTauChem 0.001;
}

EulerImplicitCoeffs
{
cTauChem 0.05;
equilibriumRateLimiter off;
}

odeCoeffs
{
solver SIBS;
eps 0.05;
scale 1;
}

When I try to run the calculation with reactingFoam it seems like there is something wrong with the library:

Create time

--> FOAM Warning :
From function dlOpen(const fileName&, const bool)
in file POSIX.C at line 1175
dlopen error : /home/fluent/OpenFOAM/fluent-2.1.x/platforms/linux64Gcc46DPOpt/lib/libmyChemistryModel.so: undefined symbol: _ZN4Foam16myChemistryModelINS_17psiChemistryModelE NS_19sutherlandTransportINS_12specieThermoINS_11ja nafThermoINS_18isobaricPerfectGasEEEEEEEE8typeName E
--> FOAM Warning :
From function dlLibraryTable:pen(const fileName&, const bool)
in file db/dynamicLibrary/dlLibraryTable/dlLibraryTable.C at line 96
could not load "libmyChemistryModel.so"

Create mesh for time = 0


Reading g
Creating combustion model

Selecting combustion model PaSR<psiChemistryCombustionModel>
Selecting psiChemistryModel myChemistryModel<gasThermoPhysics>


--> FOAM FATAL ERROR:
Unknown psiChemistryModel type myChemistryModel<gasThermoPhysics>

Valid psiChemistryModel types are:
5
(
ODEChemistryModel<constGasThermoPhysics>
ODEChemistryModel<constIsobaricGasThermoPhysics>
ODEChemistryModel<icoPoly8ThermoPhysics>
ODEChemistryModel<gasThermoPhysics>
ODEChemistryModel<isobaricGasThermoPhysics>
)

From function psiChemistryModel::New(const mesh&)
in file chemistryModel/psiChemistryModel/psiChemistryModelNew.C at line 117.

FOAM exiting


What did I miss?
matari is offline   Reply With Quote

Old   November 6, 2012, 14:44
Default
  #15
Senior Member
 
mturcios777's Avatar
 
Marco A. Turcios
Join Date: Mar 2009
Location: Vancouver, BC, Canada
Posts: 725
Rep Power: 18
mturcios777 will become famous soon enough
Sorry I hadn't seen this. Hopefully you got it resolved, but it seems like the linking hasn't been completed.
mturcios777 is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
reactingFoam Enthalpy of Reaction Problem Cyberholmes OpenFOAM 1 August 23, 2011 05:23
surface reaction rate with udf yellow-stuff FLUENT 4 January 29, 2010 13:53
help! Creating iso-surface of reaction rate James Willie FLUENT 1 September 28, 2005 21:56
Combustion reaction rate Karthick FLUENT 2 April 22, 2004 05:23
Design Integration with CFD? John C. Chien Main CFD Forum 19 May 17, 2001 15:56


All times are GMT -4. The time now is 09:14.