2016-09-25 15:11 BST

View Issue Details Jump to Notes ]
IDProjectCategoryView StatusLast Update
0000669OpenFOAM[All Projects] Bugpublic2012-10-29 12:37
Assigned Tohenry 
Product Version 
Target VersionFixed in Version 
Summary0000669: Wrong calculation of calculating Cp in addition operator in hConstThermo
DescriptionIt seems to me that operator+ in hConstThermo calculates the Cp in a wrong way (and this applies also to Hf and to other arithmetical operators).

 Consider the following code:

int main(int argc, char *argv[])
    typedef constTransport<specieThermo<hConstThermo<perfectGas> > > ThermoType;
    IFstream f("thermo.dict");
    dictionary dict(f);
    ThermoType t1(dict.subDict("specie1"));
    ThermoType t2(dict.subDict("specie2"));
    Info << "W= " << t1.W() << " " << t2.W() << " " << (t1+t2).W() << endl;
    Info << "Cp= " << t1.cp(1) << " " << t2.cp(1) << " " << (t1+t2).cp(1) << endl;

Run it on the following thermo.dict file:

   specie {
      nMoles 1;
      molWeight 1;
   thermodynamics {
      Cp 1;
      Hf 0;
   transport {
      mu 1;
      Pr 1;
   specie {
      nMoles 1;
      molWeight 0.5;
   thermodynamics {
      Cp 2;
      Hf 0;
   transport {
      mu 1;
      Pr 1;
The program will output
 W= 1 0.5 0.75
 Cp= 1 1 1.125
 However, the I expect the following output:
 W= 1 0.5 0.75
 Cp= 1 1 1

 Indeed, here two species are defined having different molar masses (1 and 0.5 kg/kmol) and different heat capacitance per kg (1 and 2), but note that they have the same heat capacitance per mole --- 1 J/(kmol K).

 The program reads both specie data, adds them, and outputs the resulting heat capacitance. As expected, cp() for both species separately is 1 (cp() returns J/(kmol K)). The heat capacitance per kilomole of mixture of these two species should therefore be 1 too, independently of the mixture composition, therefore the program should output 1, not 1.125.


 As I understand it, the problems comes from the hConstThermo addition operator that calculates the thermal capacity by the following code (hConstThermoI.H, lines 209-210):

      + ct2.nMoles()/eofs.nMoles()*ct2.Cp_
This would be correct if hConstThermo::Cp_ was the thermal capacity per mole (J/(kmol K) or J/(mol K)). However, Cp_ is the thermal capacity per kg (J/(kg K)) --- see, for example, hConstThermo::cp(), which returns J/(kmol K) and is implemented as Cp_*this->W().

 So, the correct code should be something like

      + ct2.nMoles()/eofs.nMoles()*ct2.Cp_*ct2.W())
The same applies to Hf calculations and to other arithmetical operators in hConstThermo (operator-, operator+=, operator -=).

 Note also that a similar code exists in janafThermo:
         highCpCoeffs[coefLabel] =
          + molr2*jt2.highCpCoeffs_[coefLabel];
However, here highCpCpeffs_ are per mole, not per kg (and janafThermo::cp() does not refer to W() ), so this equation is a correct one.
TagsNo tags attached.
Attached Files




henry (manager)

Thanks for the detailed bug-report.

On further consideration we have concluded that it would be better to convert the Cp_ and Hf_ to mole-based on input and back again on output in the same manner as in hPolynomialThermo which is more consistent and also resolves the issue you reported.

Resolved by commit 74eb38ab125c177165c9704f405bcbfe12cf26c8

-Issue History
Date Modified Username Field Change
2012-10-26 17:05 pkalinin New Issue
2012-10-29 12:37 henry Note Added: 0001757
2012-10-29 12:37 henry Status new => resolved
2012-10-29 12:37 henry Resolution open => fixed
2012-10-29 12:37 henry Assigned To => henry
+Issue History