Thermophysical models and chemistry models
 

Dear Foamers,

I am modifying rhoReactingFoam solver (OpenFOAM 1.6). I want implement Eddy break-up model. I have no problem with the reaction rate. My main problem is the energy equation source term. I need to access the specieThermo to get Hc() (I think it is the enthalpy of formation).
I cannot access it. I think, I am not understand how to access the thermophysical and chemistry models.
How are the thermophysical and chemistry models implemented in the solve? How can I access it?

Hi hk318i,

I'm not sure if this answers your question because I'm not sure where you need the thermophysical properties (e.g. the solver or the chemistry model template)

In the solver (as an expample I use rhoSimpleFoam, where hPsiThemro is used), the thermophsyical model is loaded in the "createFields.H" file. The imported lines are:

Info<< "Reading thermophysical properties\n" << endl;

autoPtr<basicPsiThermo> pThermo
(
basicPsiThermo::New(mesh)
);
basicPsiThermo& thermo = pThermo();

now lets say you want to get the enthalpy so the function is: "your variable = thermo.h();"

This function can be found in the file (OpenFoam/src/thermophysicalModels/basic/psiThemro/hPsiThemro.H). The function is

virtual const volScalarField& h() const
{
return h_;
}

These functions could also be in the following files: "basicThemro.H" and "psiBasicThermo.H".

So, if you need to get hc in the solver, you first have to make sure that a function like the one above (of course for hc) exist in your themrophysical model.

Regards,
Christian

hi chris

i've a question to your reply:

is it: yourvariable = thermo.h(); or: yourvariable = thermo.h_(); (with underscore like below?

and must i include:
#include basicPsiThermo
at the beginning?

thanks for your help
regards marc

Hi Marc

the correct term is "volScalarField& h = thermo.h();" (Have a look at the source code of rhoSimpleFoam). The other possibility that might work is:
"volScalarField& h = thermo.h_;". However, I think the variable h_ is protected, so it won't work.

h_ is the variable(field) used in the thermophysicalProperties and the function h() simply returns h_.


You have to include " #include "basicPsiThermo.H" ". In the case of rhoSimpleFoam, you can find the line in the file "rhoSimpleFoam.C"


Regards,
Christian

Dear Chirs,

rhoReactingFoam uses chemistry model not a thermophysical model as the following;

Info<< nl << "Reading thermophysicalProperties" << endl;
autoPtr<rhoChemistryModel> pChemistry
(
rhoChemistryModel::New(mesh)
);
rhoChemistryModel& chemistry = pChemistry();

hsReactionThermo& thermo = chemistry.thermo();

the variable which I am looking for is specieThermo.Hc()

my the Thermophysicalproperties is
thermoType hsRhoMixtureThermo<reactingMixture<gasThermoPhysic s>>;

my chemistryproperties is
rhoChemistryModel ODEChemistryModel<gasThermoPhysics>;

I have no idea what is gasThermoPhysics?

Thanks

hi hassan,

as i understood chemistry is a part of thermophysicalModels, because you need alway a "...Thermo"-file to get T,p, h, etc..
And in the file rhoReactingFoam.C there is
#include "hReactionThermo.H"
used, which then actually returns hc()

00134 //- Chemical enthalpy [J/kg]
00135 virtual tmp<volScalarField> hc() const = 0;

so you can use this

regards marc

Dear Marc,
my hReactionThermo.H member functions as the following;

// Member functions

//- Return the composition of the multi-component mixture
virtual basicMultiComponentMixture& composition() = 0;

//- Return the composition of the multi-component mixture
virtual const basicMultiComponentMixture& composition() const = 0;


// Access to thermodynamic state variables

//- Enthalpy [J/kg]
// Non-const access allowed for transport equations
virtual volScalarField& h()
{
return h_;
}

//- Enthalpy [J/kg]
virtual const volScalarField& h() const
{
return h_;
}


//- Update properties
virtual void correct() = 0;
};

I cannot find hc(). Actually all I need to access the enthalpy of formation to make the energy equ source term which is a function of reaction rate (which I calculated based on my model) and specie enthalpy of formation (like Hc() in specieThermo)

Regards Hassan

Just a quick note:
In reactingFoam, a transport equation for "total enthalpy " is solved. Total enthalpy is the sensible enthalpy plus the chemical enthalpy. Thus, the transport equation for this variable does not contain any chemical reaction source term, it is conserved in an adiabatic flame. The chemical reaction source terms are provided for the species mass fractions. Once you get the species mass fractions, the "thermo->correct();" method will calculate the heat of formation and desctruction of different species.

I am using OpenFOAM 1.6 which reactingFoam and rhoReactingFoam solve transport equation for sensible enthalpy not the total enthalpy.
Do you think I should use transport equation for total enthalpy instead of sensible enthalpy?

This is the hEqn.h in reactingFoam in Openfoam1.6:

{
fvScalarMatrix hEqn
(
fvm::ddt(rho, h)
+ mvConvection->fvmDiv(phi, h)
- fvm::laplacian(turbulence->alphaEff(), h)
==
DpDt
);

hEqn.relax();
hEqn.solve();

thermo.correct();
}

I do not see any source term, I do not think that it solves for sensible enthalpy, are you sure?

My version uses hs with source term as the following;

{
fvScalarMatrix hsEqn
(
fvm::ddt(rho, hs)
+ mvConvection->fvmDiv(phi, hs)
- fvm::laplacian(turbulence->alphaEff(), hs)
==
DpDt
+ chemistrySh
);

hsEqn.relax();
hsEqn.solve();

thermo.correct();

Info<< "T gas min/max = " << min(T).value() << ", "
<< max(T).value() << endl;
}


I am using OpenFoam1.6.x which I think it is like OpenFoam 1.7 use sensible enthalpy instead of total enthalpy :(

yes, I downloaded openfoam1.7 and I see now ! it's great, they finally make this. I needed this ! Thanks

Actually I am little bit confused about using chemistry and thermophysical models. How can I access the data? How can I change the model inside the solver? and what is the "Selector" ? :-O

For "selecting" you may read this:
http://openfoamwiki.net/index.php/Op...tion_mechanism

You cannot have access to everything, unless you change the source code for thermo and chemistry objects. However, you have access to reaction rates, heat release, ... .

so I cannot access to something like enthalpy of formation of each specie.

Hi everyone,
I am trying to simulate a incompressible, laminar fluid flow model with adding temperature field. I used air to be the inlet gas. Assuming that "basicPsiThemo" is only used for compressible, turbulence flow, hence I want to use "basicThermo" model only. And I was writing into the "createFields.H" like that:

line 1 Info<< "Reading thermophysical properties\n" << endl;
line 2
line 3 autoPtr<basicThermo> pThermo
line 4 (
line 5 basicThermo::New(mesh)
line 6 );
line 7 basicThermo& thermo = pThermo();
line 8 volScalarField& h = thermo.h(); ...

And then compile the code. An error is shown:
createFields.H: In function ‘int main(int, char**)’:
createFields.H:5: error: no matching function for call to ‘Foam::dictionary::New(Foam::fvMesh&)’
/home/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/dictionary.H:218: note: candidates are: static Foam::autoPtr<Foam::dictionary> Foam::dictionary::New(Foam::Istream&)

So what is the right one I have to declare? What does it mean about the line 5? And how can I check the declaration of basicThermo mode?

Please give me a hint?

Dear Anh,
I am new OpenFoam user and I am not fully aware about this problem.
As I know you cannot basicThermo like as you mentioned because basicThermo class does not have a selector function like basicpsiThermo.

Thanks. I got that. However, what is basicThermo function to use? Can basicpsiThermo quantities be used for incompressible, laminar flow? If not, what is the right one?

you can check the openFoamwiki below likw which descripe how to add temperature to icoFoam solver.
http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam

Dear All,

Anyone knows what is the CompType in ODEChemistryModel and other chemistry models?
I am still looking how to access the enthalpy of formation of each species. any ideas?

Regards,
Hassan

So back to my original question...
 

Sounds like a lot of good discussion came from this thread, which I'm glad about.

However, I'm still needing to calculate the reaction rates without solving the ODE. Any hints? Pretty please?

Which model do you want to use to calculate the reaction rates?

I'm using chemistry:

psiChemistryModel ODEChemistryModel<gasThermoPhysics>;

and thermo:

thermoType hPsiMixtureThermo<reactingMixture<gasThermoPhysics >>

I can access the chemistry.solve() function, which acts very much like chemistry.calculate (for ODE chemistry). I was temporarily using a really ugly fix by just replacing the contents of the solve() function with calculate(), and it works. But I have to remember to make changes when I use other solvers that require chemistry.solve(). I'm sure it has something to do with the way things are templatized, but I'm at a loss as to how to deal with it.

Thanks for the help

I find chemistry models and thermophysical very complicated specially I am trying to understand it myself. I attended openfoam foundation course but these models were not in the scope of the course.
I feel that everyone know something I totally missed. I hope someone will explain soon on wiki or here how these models work.
I am implementing EDM model and it starts to work good but I still need to access the enthalpy of formation of each specie to calculate energy equation source term. Can anyone tell me what is the chemical enthalpy hc()?

Ditto brother. I'm sure its just an inheritance thing. I was able to sort of gain access by ensuring that a pure virtual function existed for calculate (I think it was in the psiChemistry.H file), which allowed the scope to be resolved, but then complained that ther function wasn't implemented...

I guess I'll use the ugly hack for the time being...

Diffusion problem.
 

I have checked some books, when you contain Hc ( enthalpy of formation) inside h( total enthalpy) you have the reaction reaction source term there. So there is no extra term for reaction. Actually the chemical energy is released when a high energy material like CH4 is decreasing in Y( mass ratio) and a low energy material like CO2 is increasing.

But I have another problem. In this hEqn code, there should be one more term that should compensate enthalpy diffusion due to mass diffusion. I openFoam 1.7 this problem is fixed, because h is hs (sensible enthalpy ) which doesn't need that term. Compare this to equation 1.75 in "Turbulence Combustion" book by Nurbert Peters.
Actually I did one pure diffusion case, that proves reactingFoam in 1.6 has a problem.

you are right, there are no source term for reaction in total enthalpy equation, But I am using the new version which solve the sensibly enthalpy equation which contains a source term for reaction. I do not know how to calculate this term because it needs the enthalpy of formation of each specie which I cannot calculate in OpenFOAM. Do you have any idea how can I get it?
I checked equation 1.75 in Peter's book, but I am not sure if alphaEff() conceder the mass diffusion or not. As I noticed the hs and h equations are identical except the source term.

Trying to figure them out.
 

Actually I was trying to see if I could find how I can have access to Hc like other people. Since I want to use cantera in OpenFoam and that works with OF1.5 ( thanks to Markus Rheim) I had to change the code of reactingFoam to get rid of that problem. So I`m still trying to find out about different parts of this code.

About AlphaEff, I don't know. That alphaEff has to do with turbulence, but my cases are mostly laminar, so I don't know if that works for turbulent properly or not, but in my cases its giving me errors. Also you should notice that this aplha is not temperature diffusivity, it is alpha*cp or k/rho.
I have one more question about YEqn File, but that doesn't go with thread.
If YEqn.H we have muEff which acts simmilar to this alphaEff, the code is:
solve
(
fvm::ddt(rho, Yi)
+ mvConvection->fvmDiv(phi, Yi)
- fvm::laplacian(turbulence->muEff(), Yi)
==
kappa*chemistry.RR(i),
mesh.solver("Yi")
);

I don't know why it is muEff(), I mean it should be effective mass concentration coefficient D or DEff. Do you anything about this?

I have some cases, If you would like to see, I can send them to you, just send me your email.

sure my email is hassan.kassem@gmail.com

http://jeacfm.cse.polyu.edu.hk/downl..._MarzoukOA.pdf
I think this paper is useful. They used reactingFoam. you can check the governing equations.

Hi,
can you please help me to understand what "thermo->correct()" do in hEqn.H ?

thanks
Silvano

Hi

I am just starting to explore openfoam. I have a quick question in the following.

autoPtr<basicThermo> pThermo
I understand that in above line pThermo is an object of the template class autoptr.
basicThermo& thermo = pThermo();
But I don't understand is that how the object pThermo is used as pThermo(). I understand that thermo is a reference variable of type basicThermo. Please explain the use of pThermo().

Regards

Dear All

I am going to run a case includes some liquid species that react with each other. I am going to use the reactingFoam as a solver.
All species have constant proprties and I want to use the following thermo type:
thermoType hsPsiThermo<reactingMixture<constTransport<specieT hermo<eConstThermo<icoPolynomial>>>>>>

but when I set it in the constat/thermophysicalProperties file and run the case, I get this error:

------------------------------------------------------------------------------------
--> FOAM FATAL ERROR:
Inconsistent thermo package selected:

hsPsiThermo<reactingMixture<constTransport<specieT hermo<eConstThermo<icoPolynomial>>>>>>

Please select a thermo package based on gasThermoPhysics. Valid options include:

3
(
hsPsiMixtureThermo<singleStepReactingMixture<gasTh ermoPhysics>>
hsPsiMixtureThermo<multiComponentMixture<gasThermo Physics>>
hsPsiMixtureThermo<reactingMixture<gasThermoPhysic s>>
)



From function autoPtr<hsCombustionThermo> hsCombustionThermo::NewType(const fvMesh&, const word&)
in file combustionThermo/hsCombustionThermo/hsCombustionThermoNew.C at line 116.

FOAM exiting
-------------------------------------------------------------------------------------

What should I do to use the reactingFoam for constant propertie species?

Thanks in advance for your help.