CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Programming & Development (
-   -   computation of K, cp, h and s ?? (

megacrout September 19, 2011 12:27

computation of K, cp, h and s ??
Hi everyone,

I am currently looking at the code behind combustion solvers, especially reactingFoam. I am interested in the way reversible reactions are implemented. I figured out that the reaction rate of the backward reaction is computed using the equilibrium constant K which is computed as

K = exp( - nMoles * g(T) / (RR * T) )

where the gibbs enthalpy is in turn computed as

g(T) = h(T) - T*s(T).

This means

K = exp ( - nMoles * (h(T) - T*s(T)) / (RR*T) )

where RR represents the reaction rate. This does not seem very close to

K = exp ( - (h(T) - T*s(T)) / (R*T) )
where R is the gas constant, which is supposed to be the correct expression for the equilibrium constant (right ?).

Futhermore, I just CANT find in the online C++ documentation where cp(T), h(T) and s(T) are computed. I know that it is based on the Janaf parameters from which the heat capacity is also derived ( But I cannot find where this computation is being done, where the equation is implemented! And since all other thermophysical parameters are then derived from those three (as mentioned in the header of specieThermo.H for example), I really want to know where to find those equations in OF.

This is not the first time I spend a few hours looking for a specific equation so thatīd be a really great help if you can give me an answer. Even just a hint ahah!!



Chris Lucas September 21, 2011 05:18


I'm not sure what you want to know?

The specific values for cp ... are calculated in specieThermo using the molar values.

The molar values are calculated (depending in the model you use) e.g. in hConst or hPolynomial or janaf ...

Best Regards,

megacrout September 22, 2011 03:58

Hi Chris,

Thanks for your answer.
Iīm trying to know how cp(T), h(T) and s(T) are computed (files AND formulas)? I canīt see what you mentioned: specieThermo contains the formulas for computing all thermophysicals parameters BUT the three Iīm looking for. Or I did miss something.
Where is - in the case of cp - the function hConst/hPolynomial/janaf being called? I did not find it in specieThermo (I must admit Iīm a C++ newbie, though).

As I mentioned, it is not the first time I get lost in the C++ documentation of OpenFOAM and in specieThermo is another good example: where is the parameter W() computed? It is used in ThermoI.H but I canīt find where its formula lies or even where it is being called.

Best regards,


adhiraj September 22, 2011 10:47

I think it was janafthermo.

Chris Lucas September 23, 2011 02:39


for W() have a look at the files in


for cp, h ... look at




Best Regards,

megacrout September 26, 2011 11:53

Hey Christian,

Thanks A LOT for the hints. I say "hints" because Iīve had to look for more details quite a long time before I could get what you mean.

The missing link were the files hsCombustionThermos.C and newhsCombustionThermo.C (which seems to have been renamed hsCombustionThermoNew.C in OF 2.0.1). Without the former file, there is no link between solver reactingFoam and the files specie or janaf/hConst/etc...

Now, I have a few more questions for you, I would greatly appreciate if you can further help me:
1) reactingFoam relies on the psiChemistrymodel (s. createFields.H), can one of the other models (rho, ODE, basic) be used as an option (i.e. without modifying the solverīs code)?
2) same question about hsCombustionThermo (s. createFields.H): could hCombustionModel, hhuCombustionModel or mixtureThermos be used instead, without altering reactingFoamīs structure, i.e. as an option?
3) Because reactingFoam relies on psiChemistrymodel and hsCombustionThermo, only 2 options can be given as thermoType in thermophysicalProperties:

hsPsiMixtureThermo<multiComponentMixture<gasThermo Physics>>
hsPsiMixtureThermo<reactingMixture<gasThermoPhysic s>>

and only 2 options can be given as psiChemistryModel in chemistryProperties:


Is that right? If not, how can I implement another thermoType? Any attempt gave me error messages such as

"Inconsistent thermo package selected..."
"Unknown hsCombustionThermo type...""
4) If what I just wrote is right, then I still donīt see where hConst/specieThermo/janaf is called.

Finally, regarding my first post (still unanswered concerning reactions): I figured out in the meanwhile that RR is used for both reacting rate and gas constant depending on the file it is used in. However, I still do not get why nMoles appears in the mentioned equation. Can you help me?



megacrout October 11, 2011 06:02

"Pretty please with sugar on top." (Harvey Keitel as Winston Wolf in Pulp Fiction)

All times are GMT -4. The time now is 07:11.