cellMixture: multicomponent 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 massbased defined, meanwhile in specieThermo h is molarbased and H is massbased. 3. In hMixtureThermo::h (massbased) h = ForAll(i) { += Y[i] / species[i].W() * speciesData[i].H() }; where, as far as I understand, this H[i] is the massbased 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 cellMixturedefinition of multiComponentMixture by: cellMixture = ForAll (i) { cellMixture += Y[i] * speciesData[i] }; Is there really a bug..or I am missing something?? Regards! 
Quote:
Quote:
This is an inconsistency in the naming convention. Quote:
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> Quote:
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 forloop does in Q1. I hope all is clearer now 
Finally, I got it.
For those who could have the same doubt, I copy the keypoint (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 
So, we could only set janafThermo? And what about the hconstThermo?

Dear Maria and Niklas,
sorry for reviving this old thread but I am currently looking into the same issue and haven't fully understood how the mixtures are calculated in OpenFOAM, yet. I work with release 2.3.x and am looking at the equation given below: forAll (i) { cellMixture += Y[i] / species[i].W() * speciesData[i] } As I understand, speciesData[i] and cellMixture may for example be the coefficients of the JANAF polynomials: speciesData[i] for component i and cellMixture for the mixture. However, I feel like the molar mass of the mixture is missing in the equation. As it is now, the value of the coefficient related to the mixture should be scaled by the molar mass of the mixture or do I forget something? Moreover, the dimensions of "speciesData" and "cellMixture" should be different. Is this corrected at a later stage within the mixture property calculation? I hope you can help me find my error in reasoning! Thanks already! 
Understanding cellMixture: OpenFOAM incorrect density calculation in rhoReactingTherm
2 Attachment(s)
Hi everyone,
LONG POST, PLEASE READ! At present, I am working with rhoReactingFoam in OpenFOAM v 5.0. On doing some digging into the solver, I was confused about the cellMixture calculations. Here are my observations, and I would request the senior members to please take some time and help the new Foamers understand this method. Disclaimer: The calculations seem to have been modified from previous versions of OpenFOAM so I will be posting everything afresh. 1) Where does rhoReactingFoam update thermophysical properties? The only viable candidate was in EEqn.H Code:
{ rhoReactionThermo.H Code:
class rhoReactionThermo (attached with this post) of my case in the constant directory, I use the heRhoThermo ($FOAM_SRC/thermophysicalModels/basic/rhoThermo/heRhoThermo.*) type of rhoThermo. The correct function has been defined in this class. It simply calls another function in this class called calculate(), which is shown below: heRhoThermo.C Code:
template<class BasicPsiThermo, class MixtureType> This class does not implement the cellMixture function but instead inherits it from the multiComponentMixture. Code:
template<class ThermoType> multiComponentMixture.C Code:
template<class ThermoType> a = 1; b = 2; a += b; => a = 3; But, in OpenFOAM, they have redefined the way this and other operators behave for the thermoType class (mixture is a thermoType variable). 6) In my case again, thermophysicalProperties defines thermoType thermo as hPolynomial. So, I go to $FOAM_SRC/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H, where this operator has been defined. hPolynomialThermoI.H Code:
template<class EquationOfState, int PolySize> 7) The bold part calls another operator for equationofstate (in my case icoPolyomial) $FOAM_SRC/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H Code:
template<class Specie, int PolySize> Code:
inline void specie::operator+=(const specie& st) Single Cell: Details: Specie 1: Molecular Weight = 36.46, Mass Fraction = 0.025, Density = (458.6 0 0 0 0 0 0 0) // icoPolynomial definition Specie 2: Molecular Weight = 112.56, Mass Fraction = 0.60, Density = 969.8 Specie 3: Molecular Weight = 98.92, Mass Fraction = 0.375, Density = 1089.4 Quick density calculation: Quote:
where Yi is mass fraction of each specie, and rhoi is density of each specie. For our example it comes out to be 982.874 But OpenFOAM calculates it to be 1001.87 (Output Log attached with this post). This number matches exactly if we calculate, Quote:
Interestingly, multiComponentMixture class has another function cellVolMixture, which seems to calculate the density correctly. But this is not the function called to define the mixture type in rhoReactionFoam. Code:
template<class ThermoType> ## THE END ## 
Hi Raunak,
Actually, the operator += in the calculation of mixture_ is defined in the class sutherlandTransport (or other transport class). You can have a look at this post: https://openfoam.top/en/thermodynamicLIB/. Sorry its in Chinese, but I think the UML will help. As you can see in the UML, multiComponentMixture<ThermoType> needs a template, it can be gasHThermoPhysics, i.e. Code:
sutherlandTransport<species::thermo<janafThermo<perfectGas<specie>>,sensibleEnthalpy>> Code:
polynomialTransport<species::thermo<hPolynomialThermo<icoPolynomial<specie, 8>,8>,sensibleEnthalpy>,8> 
Hi Raunak,
can you share the code you used to probe the density calculation in a singlecell? I think I'm noticing the same issue, and would like to doublecheck Thank you very much 
Hey,
I dealt with this issue about a year ago and can't really recall the exact solution I found for this. It was proprietary work so I don't even have my code to go back and look. If I recall, I played around a bit with this for a while but then skipped trying to understand the existing thermophysics model and defined one for my own case. Added it as another run time selectable model and used it in my simulations. Sorry, I am not much help on this one. Hope you will figure it out. 
All times are GMT 4. The time now is 17:58. 