CFD Online Discussion Forums - hConstThermo: calculating Cp in addition operator

Hello all,

It 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:

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;



    return(0);

}

Run it on the following thermo.dict file:

Code:

specie1{

   specie {

      nMoles 1;

      molWeight 1;

   }

   thermodynamics {

      Cp 1;

      Hf 0;

   }

   transport {

      mu 1;

      Pr 1;

   }

}

specie2{

   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):

Code:

ct1.nMoles()/eofs.nMoles()*ct1.Cp_

      + 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

Code:

         (ct1.nMoles()/eofs.nMoles()*ct1.Cp_*ct1.W()

      + ct2.nMoles()/eofs.nMoles()*ct2.Cp_*ct2.W())

        /eofs.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:

Code:

         highCpCoeffs[coefLabel] =

            molr1*jt1.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.

Is my understanding correct and the behavior is really wrong, or is there any detail that I did not catch?

(In fact, I was trying to simulate chemistry with hConstThermo instead of janafThermo, but found that thermodynamical properties of the mixture are calculated wrongly.)

P.S. I have tried to submit this as a bug report at www.openfoam.org/bugs/ . I have registered there, filled the "Enter Report Details" form, clicked "Submit report", and got the following error:

Quote:

Forbidden
You don't have permission to access /mantisbt/bug_report.php on this server. Additionally, a 404 Not Found error was encountered while trying to use an ErrorDocument to handle the request.

I have then tried to find any way to contact the bug report system admin, or another way to contact the developers or the site admins, but found none. Can anybody suggest a way to do this?