[matari]

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

[elvis]

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

[matari]

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]

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]

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.

Originally Posted by matari

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.

[mturcios777]

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.

[matari]

Hello,

thank you for your reply. The tutorial is really good and helps a lot!

[ghas]

Hi matari,

Originally Posted by matari

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

[mturcios777]

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!

[matari]

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?

[mturcios777]

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?

[matari]

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

[mturcios777]

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!

[matari]

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?

[mturcios777]

Sorry I hadn't seen this. Hopefully you got it resolved, but it seems like the linking hasn't been completed.

[iiigoriii]

Dear users,

I really don't understand one thing (maybe more ) and please I need help. I modified both ArrheniusReactionRate.C /.H

src\thermophysicalModels\specie\reaction\reactionRate\ArrheniusReactionRate

I changed ArrheniusReactionRate.H equation, wmake (also wmake libso) pass without errors. I can also found new libspecie.dll in lib directory. So I guess everything is OK.

But, I have a testing example and I see that nothing changed. My modified Arrhenius equation wont work. Please have you some explanation? I read some general manuals, also try search at forum... (solver modifications works without problems). Thanks, ii

[iiigoriii]

I am a little despairing, I also did Allwmake for whole OpenFOAM-dev and nothing. I have all libraries with new date (today) and nothing. Why changes in reaction rate equation do nothing?

despairing, ii

[iiigoriii]

OK, sorry for desperate last two spams. Problem was, as usual, between monitor and keyboard. I tried whole process under linux and everything works perfectly. My unsuccessful attempts were made under win version....

I created my own reaction rate, and everything works perfectly.
ii