CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Bugs

cellMixture: multi-component mixture

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree6Likes
  • 2 Post By mgc
  • 1 Post By niklas
  • 1 Post By mgc
  • 1 Post By hg2lf
  • 1 Post By raunakbardia

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 30, 2009, 09:53
Default cellMixture: multi-component mixture
  #1
mgc
New Member
 
Maria
Join Date: Apr 2009
Posts: 12
Rep Power: 13
mgc is on a distinguished road
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!
Zhiheng Wang and Kummi like this.

Last edited by mgc; July 1, 2009 at 04:14.
mgc is offline   Reply With Quote

Old   June 30, 2009, 12:18
Default
  #2
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 694
Rep Power: 25
niklas will become famous soon enoughniklas will become famous soon enough
Quote:
Originally Posted by mgc View Post
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 View Post
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 View Post
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 View Post
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
mickbatti likes this.
niklas is offline   Reply With Quote

Old   June 30, 2009, 13:38
Smile
  #3
mgc
New Member
 
Maria
Join Date: Apr 2009
Posts: 12
Rep Power: 13
mgc is on a distinguished road
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
Zhiheng Wang likes this.

Last edited by mgc; July 15, 2009 at 05:57.
mgc is offline   Reply With Quote

Old   May 12, 2013, 23:06
Default
  #4
New Member
 
何刚
Join Date: Jan 2011
Posts: 25
Rep Power: 11
hg2lf is on a distinguished road
So, we could only set janafThermo? And what about the hconstThermo?
Zhiheng Wang likes this.
hg2lf is offline   Reply With Quote

Old   April 27, 2017, 11:18
Default
  #5
New Member
 
Join Date: Dec 2013
Posts: 11
Rep Power: 8
Daniel_P is on a distinguished road
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!
Daniel_P is offline   Reply With Quote

Old   August 9, 2018, 17:59
Default Understanding cellMixture: OpenFOAM incorrect density calculation in rhoReactingTherm
  #6
Member
 
Raunak Bardia
Join Date: Jan 2015
Posts: 32
Rep Power: 7
raunakbardia is on a distinguished road
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:
{
...

    EEqn.relax();

    fvOptions.constrain(EEqn);

    EEqn.solve();

    fvOptions.correct(he);

    thermo.correct(); // properties are probably calculated in this

    Info<< "min/max(T) = "
        << min(T).value() << ", " << max(T).value() << endl;
}
2) thermo is a variable type of rhoReactionThermo ($FOAM_SRC/thermophysicalModels/reactionThermo/rhoReactionThermo). This class does not implement the correct() function we are looking for but inherits the function from rhoThermo.

rhoReactionThermo.H
Code:
class rhoReactionThermo
:
public rhoThermo
{
...
3) Based on my definition of thermophysicalProperties dictionary
(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>
void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::calculate()
{
    const scalarField& hCells = this->he();
    const scalarField& pCells = this->p_;

    scalarField& TCells = this->T_.primitiveFieldRef();
    scalarField& psiCells = this->psi_.primitiveFieldRef();
    scalarField& rhoCells = this->rho_.primitiveFieldRef();
    scalarField& muCells = this->mu_.primitiveFieldRef();
    scalarField& alphaCells = this->alpha_.primitiveFieldRef();

    forAll(TCells, celli)
    {
        const typename MixtureType::thermoType& mixture_ =
            this->cellMixture(celli);

        TCells[celli] = mixture_.THE
        (
            hCells[celli],
            pCells[celli],
            TCells[celli]
        );

        psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);
        rhoCells[celli] = mixture_.rho(pCells[celli], TCells[celli]);

        muCells[celli] = mixture_.mu(pCells[celli], TCells[celli]);
        alphaCells[celli] = mixture_.alphah(pCells[celli], TCells[celli]);
    }

 ..... Handle Boundary Patches
4) Turns out that this calculate function also does not really compute the thermophysical properties. It simply creates a thermotype variable from the mixture class using the "cellMixture" function (Above description should e helpful for the new FOAMers). Again looking at the thermophysicalProperties dictionary I use in my case, I have defined my mixture as a reactingMixture ($FOAM_SRC/thermophysicalModels/reactionThermo/mixtures/reactingMixture/).
This class does not implement the cellMixture function but instead inherits it from the multiComponentMixture.
Code:
template<class ThermoType>
class reactingMixture
:
    public speciesTable,
    public autoPtr<chemistryReader<ThermoType>>,
    public multiComponentMixture<ThermoType>,
    public PtrList<Reaction<ThermoType>>
{
...
5) cellMixture function is defined in the multicomponentMixture class ($FOAM_SRC/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/)

multiComponentMixture.C
Code:
template<class ThermoType>
const ThermoType& Foam::multiComponentMixture<ThermoType>::cellMixture
(
    const label celli
) const
{
    mixture_ = Y_[0][celli]*speciesData_[0];   
/* mixture is of datatype ThermoType and is initialized with speciesdata of
1st specie multiplied by its mass fraction*/

    for (label n=1; n<Y_.size(); n++)
    {
        mixture_ += Y_[n][celli]*speciesData_[n];
/* All the weighted species data of remaining species is added to the mixture
variable. Now, any property needed for the mixture can directly be accessed from the mixture itself. */
    }

    return mixture_;
}
The interesting part follows. The operator "+=" does not exactly do what we think it does. For eg. Usually,
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>
inline void Foam::hPolynomialThermo<EquationOfState, PolySize>::operator+=
(
    const hPolynomialThermo<EquationOfState, PolySize>& pt
)
{
    scalar Y1 = this->Y();

    EquationOfState::operator+=(pt);

    if (mag(this->Y()) > SMALL)
    {
        Y1 /= this->Y();
        const scalar Y2 = pt.Y()/this->Y();

        Hf_ = Y1*Hf_ + Y2*pt.Hf_;
        Sf_ = Y1*Sf_ + Y2*pt.Sf_;
        CpCoeffs_ = Y1*CpCoeffs_ + Y2*pt.CpCoeffs_;
        hCoeffs_ = Y1*hCoeffs_ + Y2*pt.hCoeffs_;
        sCoeffs_ = Y1*sCoeffs_ + Y2*pt.sCoeffs_;
    }
}
While reading this operator, let's go back to our a+b example above. In this function "this" pointer is a, and "pt" is b. This is how you read any of the other operators.

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>
inline void Foam::icoPolynomial<Specie, PolySize>::operator+=
(
    const icoPolynomial<Specie, PolySize>& ip
)
{
    scalar Y1 = this->Y();
    Specie::operator+=(ip);

    if (mag(this->Y()) > SMALL)
    {
        Y1 /= this->Y();
        const scalar Y2 = ip.Y()/this->Y();

        rhoCoeffs_ = Y1*rhoCoeffs_ + Y2*ip.rhoCoeffs_;
    }
}
8) A new call to another specie operator leads to the END OF THIS TRAIL .

Code:
inline void specie::operator+=(const specie& st)
{
    const scalar sumY = Y_ + st.Y_;
    if (mag(sumY) > SMALL)
    {
        molWeight_ = sumY/(Y_/molWeight_ + st.Y_/st.molWeight_);
    }

    Y_ = sumY;
}
So, now that we have all the pieces to figure out how OpenFOAM calculates any thermophysical property for a mixture. I took an example.
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:
rho_mix = 1/(summation_over_all_species(Yi/rhoi))
,
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:
rho_mix_incorrect = summation_over_all_species(rhoi * Yi)
. If you run through the math of the operators we described in points 6, 7, and 8, you can see that we do run into this incorrect calculation.


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>
const ThermoType& Foam::multiComponentMixture<ThermoType>::cellVolMixture
(
    const scalar p,
    const scalar T,
    const label celli
) const
{
    scalar rhoInv = 0.0;
    forAll(speciesData_, i)
    {
        rhoInv += Y_[i][celli]/speciesData_[i].rho(p, T);
    }

    mixtureVol_ =
        Y_[0][celli]/speciesData_[0].rho(p, T)/rhoInv*speciesData_[0];

    for (label n=1; n<Y_.size(); n++)
    {
        mixtureVol_ +=
            Y_[n][celli]/speciesData_[n].rho(p, T)/rhoInv*speciesData_[n];
    }

    return mixtureVol_;
}
Those who have worked with this code or this problem before, your inputs will be very valuable. Density calculation seems to be incorrect but it is hard to believe that statement. Please correct me if I am wrong.

## THE END ##
Attached Files
File Type: txt thermophysicalProperties.txt (1.2 KB, 11 views)
File Type: docx OpenFOAM_log_singlecellcase.docx (14.8 KB, 14 views)
mickbatti likes this.
raunakbardia is offline   Reply With Quote

Old   August 12, 2018, 22:28
Default
  #7
Senior Member
 
zhangyan's Avatar
 
Yan Zhang
Join Date: May 2014
Posts: 108
Rep Power: 8
zhangyan is on a distinguished road
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>>
or it can be icoPoly8HThermoPhysics, i.e.
Code:
polynomialTransport<species::thermo<hPolynomialThermo<icoPolynomial<specie, 8>,8>,sensibleEnthalpy>,8>
Ref: https://github.com/OpenFOAM/OpenFOAM...PhysicsTypes.H
__________________
https://openfoam.top/en
zhangyan is offline   Reply With Quote

Old   July 12, 2019, 17:28
Default
  #8
New Member
 
Join Date: Jul 2012
Posts: 7
Rep Power: 10
mickbatti is on a distinguished road
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
mickbatti is offline   Reply With Quote

Old   July 12, 2019, 17:51
Default
  #9
Member
 
Raunak Bardia
Join Date: Jan 2015
Posts: 32
Rep Power: 7
raunakbardia is on a distinguished road
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 thermo-physics 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.
raunakbardia is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
cellMixture: multi-component mixture mgc OpenFOAM 2 July 15, 2009 05:59
Improvement to the update of T_ et al in hMixtureThermoC and hhuMixtureThermoC richpaj OpenFOAM Bugs 2 October 9, 2008 22:14
multi component ranap Main CFD Forum 1 October 8, 2008 03:00
multi component diffusivity mahut FLUENT 0 October 10, 2007 04:05
doubt about background fluid in multico mixture... Giosuč Siemens 0 June 11, 2004 05:18


All times are GMT -4. The time now is 03:55.