[mgc]

Hi !

I think I have found a bug in the "calculate()" function of a multicomponent mixture.
Let me go step by step:

1.- For multiComponentMixtures OpenFOAM calculates the cellMixture as:

cellMixture = ForAll (i)
{
cellMixture += Y[i] / species[i].W() * speciesData[i]
};

2.- Enthalpy in hCombustionThermo is mass-based defined, meanwhile in
specieThermo h is molar-based and H is mass-based.

3.- In hMixtureThermo::h (mass-based)

h = ForAll(i)
{
+= Y[i] / species[i].W() * speciesData[i].H()
};

where, as far as I understand, this H[i] is the mass-based enthalpy calculated in specieThermo.

4.- Having a look at the units of the hMixtureThermo::h
=> [m^2 s^(-2)] != [-][Kmol_i Kg_i^(-1)]*[m^2 s^(-2)]

Which IS DIMENSIONALLY INCORRECT!! ...

It does the same for calculating alpha/mu of the mixture, which are
defined also in a MASS basis. So:


I guess it should be fixed replacing the current cellMixture-definition of multiComponentMixture by:

cellMixture = ForAll (i)
{
cellMixture += Y[i] * speciesData[i]
};


Is there really a bug..or I am missing something??



Regards!

[niklas]

Originally Posted by mgc

1.- For multiComponentMixtures OpenFOAM calculates the cellMixture as:

cellMixture = ForAll (i)
{
cellMixture += Y[i] / species[i].W() * speciesData[i]
};

true, speciesData contain the polynomial coefficients for the individual components and cellMixture will after this contain the polynomial coefficients for the multicomponent mixture.

Originally Posted by mgc

2.- Enthalpy in hCombustionThermo is mass-based defined, meanwhile in
specieThermo h is molar-based and H is mass-based.

h in hCombustionThermo is indeed per mass, whereas it is per mole in specieThermo.
This is an inconsistency in the naming convention.

Originally Posted by mgc

3.- In hMixtureThermo::h (mass-based)

h = ForAll(i)
{
+= Y[i] / species[i].W() * H[i]
};

where, as far as I understand, this H[i] is the mass-based enthalpy calculated in specieThermo.

from where did you get this formula?
I cannot find it anywhere.
It looks like you have replaced speciesData in Q1 with H, which is incorrect.

All I can find in there that is fairly similar is this (which is correct):

136 template<class MixtureType>
137 Foam::tmp<Foam::scalarField> Foam::hMixtureThermo<MixtureType>::h
138 (
139     const scalarField& T,
140     const labelList& cells
141 ) const
142 {
143     tmp<scalarField> th(new scalarField(T.size()));
144     scalarField& h = th();
145
146     forAll(T, celli)
147     {
148         h[celli] = this->cellMixture(cells[celli]).H(T[celli]);
149     }
150
151     return th;
152 }

Originally Posted by mgc

4.- Having a look at the units of the hMixtureThermo::h
=> [m^2 s^(-2)] != [-][Kmol_i Kg_i^(-1)]*[m^2 s^(-2)]

Which IS DIMENSIONALLY INCORRECT!! ...

It does the same for calculating alpha/mu of the mixture, which are
defined also in a MASS basis. So:


I guess it should be fixed replacing the current cellMixture-definition of multiComponentMixture by:

cellMixture = ForAll (i)
{
cellMixture += Y[i] * speciesData[i]
};


Is there really a bug..or I am missing something??

I think there is an inconsistency in the naming convention which could be improved,
but all the code is correct.

I get the feeling you have started in the wrong end.
Looking at how nasa polynomials are defined, which is the speciesData.
(or it could just as well be janaf)
http://www.me.berkeley.edu/gri_mech/data/nasa_plnm.html

H/RT = a1 + a2 T /2 + a3 T^2 /3 + a4 T^3 /4 + a5 T^4 /5 + a6/T

unit for R is J/mol*K and the right hand side is dimensionless.
so H is in J/mol.
and if you want to calculate the enthalpy for the mixture you just add these
coefficients together on a molar fraction basis.
Thus the a1 coefficient for the mixture is sum of all its components = sum (X_i*a1_i)
That is what the for-loop does in Q1.

I hope all is clearer now

[mgc]

Finally, I got it.

For those who could have the same doubt, I copy the key-point (in my opinion) to understand how the "cellMixture" works:

//--------------------------------------------------------------------------------------------------------------

template<class ThermoType>
const ThermoType& multiComponentMixture<ThermoType>::cellMixture
(
const label celli
) const
{
mixture_ = Y_[0][celli]/speciesData_[0].W()*speciesData_[0];

for (label n=1; n<Y_.size(); n++)
{
mixture_ *+=* Y_[n][celli]/speciesData_[n].W()*speciesData_[n];
}

return mixture_;
}

//--------------------------------------------------------------------------------------------------------------

template<class equationOfState>
inline void janafThermo<equationOfState>::*operator+=*
(
const janafThermo<equationOfState>& jt
)
{
scalar molr1 = this->nMoles();

equationOfState:: operator+=(jt);

molr1 /= this->nMoles();
scalar molr2 = jt.nMoles()/this->nMoles();

Tlow_ = max(Tlow_, jt.Tlow_);
Thigh_ = min(Thigh_, jt.Thigh_);
Tcommon_ = molr1*Tcommon_ + molr2*jt.Tcommon_;

for
(
register label coefLabel=0;
coefLabel<janafThermo<equationOfState>::nCoeffs_;
coefLabel++
)
{
highCpCoeffs_[coefLabel] =
molr1*highCpCoeffs_[coefLabel]
+ molr2*jt.highCpCoeffs_[coefLabel];

lowCpCoeffs_[coefLabel] =
molr1*lowCpCoeffs_[coefLabel]
+ molr2*jt.lowCpCoeffs_[coefLabel];
}
}


Regards!

María

[hg2lf]

So, we could only set janafThermo? And what about the hconstThermo?