CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Bugs (http://www.cfd-online.com/Forums/openfoam-bugs/)
-   -   cellMixture: multi-component mixture (http://www.cfd-online.com/Forums/openfoam-bugs/65905-cellmixture-multi-component-mixture.html)

 mgc June 30, 2009 09:53

cellMixture: multi-component mixture

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 June 30, 2009 12:18

Quote:
 Originally Posted by mgc (Post 220966) 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.

Quote:
 Originally Posted by mgc (Post 220966) 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.

Quote:
 Originally Posted by mgc (Post 220966) 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):
Code:

```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 }```
Quote:
 Originally Posted by mgc (Post 220966) 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 June 30, 2009 13:38

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 May 12, 2013 23:06

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

 All times are GMT -4. The time now is 17:31.