[S.Colucci]

Hi,
I'm looking into thermophysicalModels to understand how mixture properties are calculated in multiphase solvers (in particular in reactingTwoPhaseEulerFoam). I suppose that the mixture thermophysical properties are calculated in multicomponentMixture.C:

--------------------------------------------------------------------------------------------------------------------
template<class ThermoType>
const ThermoType& Foam::multiComponentMixture<ThermoType>::cellMixtu re
(
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_;
}
--------------------------------------------------------------------------------------------------------------------

There are few things that I do not understand:

1) if you print speciesData (as well as mixture) you find an array made by 7 elements, the 1st is the component's (or phase for mixture) name, 2nd and 3rd, respectively, the number of moles and the molecular weight (nMoles and MolWeight in thermophysicalProperties.*), what about the other ones?

2) I do not understand why the mass (or volume?) fraction of each component (Y_[n]) is divided by the molecular weight (mixture_ += Y_[n][celli]/speciesData_[n].W()*speciesData_[n]). How can it be dimensionally consistent?

Thanks in advance

Simone

[marupio]

Hi Simone,

OpenFOAM is designed to be easy to develop for on the higher levels, e.g. the application, or top level libraries. But when you delve deep into it, you will find areas where things get really complex to try to make the top level more intuitive.

speciesData_ is a list of "ThermoType"s. mixture_ and mixuterVol_ are also ThermoTypes. The nature of a ThermoType varies depending on the thermo package selected at run time (there are hundreds). So it will not be consistent, what data it holds. It is implemented as a curiously recurring template pattern, which is designed to be fast and efficient... but it makes it difficult to tell what anything does.

A ThermoType contains data about a specie of fluid (or solid also), and also contains functions that calculate different values, given different inputs, such as pressure and density (equation of state), viscosity, and so on. ThermoTypes are only relevant for a single cell. That means, you will find somewhere on a higher level, a forAll(cells) when calculating thermodynamic properties (see hePsiThermo, for instance).

The speciesData_ is just a holder for a bunch of ThermoTypes. Depending on the mixture that inherits this (e.g. homogeneousMixture, reactingMixture), these species could be different things. species_[0] could be the raw fuel, species_[1] could be the air. They do not necessarily mean that this mixture is made up of these species. In reactingMixture, they are the constituent components of the mixture.

mixture_ and mixtureVol_ are temporary memory space holding the current mixture properties. multiComponentMixture (or the derived type, eg homogeneousMixture) will fill these, and then the top level thermo package (eg hePsiThermo) will use the total mixture to produce relevant fluid fields. mixture_ mixes by molar fraction, and mixtureVol_ mixes by volume fraction, I believe.

Now to your question, you are looking at the mixture code.

mixture_ += Y_[n][celli]/speciesData_[n].W()*speciesData_[n];

These are all overloaded operators, designed to make it look like we're simply mixing them together. Working on the ThermoType we have +=, and *. The / only operates on scalars. The net effect is:

ThermoType += scalar * ThermoType

You have to look through the various components in the specie directory to see how these operator overloads work. But the multiplication only works on the moles, (specie:perator*) I believe. And the += applies a before-and-after molar fraction calculation (e.g. rhoConst:perator+=).

Best of luck!

[Neka]

Dear OF users,

I have exactly the same question as Simone had.

David explained it allready, and gave a lot of very useful information, especially about the basic concept of how ThermoTypes work. However, due to my poor C++ and OF knowledge I still have problems with units multiComponentMixture.C. And again, this is about this code part:

mixture_ += Y_[n][celli]/speciesData_[n].W()*speciesData_[n]

I also found another thread about the same problem:
cellMixture: multi-component mixture

As far as I know Y is the mass fraction ([kg/kg]) and is dimensionless. speciesData_[n].W() gives me the molecular weight of the species[n] in a cell. The units of the molecular weight W are [kg/kmol]. If I need a mixture for say energy E, which is given in [J/kg], I get the following units for my mixture: cellMixture += [dimless] / [kg/kmol] * [J/kg]. I get [(kmol J)/(kg^2)].

So, obviously, I understand something wrong here. If I get my energy/enthalpy from the janafThermo, the units of cp are in [J/(kmol K)] and for e/h in [J/kmol]. This is because RR(), defined in specie.C, is in [J/(kmol K)]. In specieThermo e is converted to E by e/W (W beeng molecular weight here), so in units [J/kmol]/[kg/kmol] = [J/kg].

So, my mixture must have [J/kg] as well.

David says, that

mixture_ += Y_[n][celli]/speciesData_[n].W()*speciesData_[n]

is like :
ThermoType += scalar * ThermoType

But my scalar has wrong units, namely the units of Y/W.

In the thread, in the link above, Maria said that the key-point of understanding this problem lies in the definition of the operator "+=" in janafThermoI.H.
As far as I understand David said the similar thing.

My OF knowledge is not jet sufficient to understand the code part where this "+=" operator is defined:

template<class equationOfState>
inline void Foam::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];
    }
}

Could anybody explain me a bit more precisely, what I dont see in:

mixture_ += Y_[n][celli]/speciesData_[n].W()*speciesData_[n]

Thanks in advance,

Alex

[Zhiheng Wang]

Hi,
I want to calculate the species diffusivity D (fick law) on basis of binary diffusivity I am knowing the formulation ,
My problem is how to extract the data.

As in creatFields.H spicies defined as
basicMultiComponentMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y();

every property kappa, mu, alpha,or rho for example is called as composition.rho(i,1e5,T)

How can I implement D (molecular diffusion in same form).

Or if any one suggest how to create PtrList<volScalrField>& D = some variable which can calculate this.

[Neka]

Hello Zhiheng,

I think your question is wrong in this thread. I've seen you opend allready a thread for your question yesterday:
Calcolation of spicies diffusivity

So, why do you post it here again.

[Neka]

Dear OF users,

after spending some time with delving into OF operator definitions and asking for some help of an advanced OF user, I would like to post, what I’ve found out concerning my question above. Now (I hope) I understand much better the answer of David Gaden in post #2. BTW, I use foam-extend and in my concrete case I use the thermoType, which includes multiComponentMixture, janafThermo and perfectGas classes. But I think It should be very similar for OpenFOAM. Let’s say I would like to have an enthalpy for a species mixture.
In foam-extend in e.g. in hPsiMixtureThermo.C that would be:
In multiComponentMixture we can find the definition of the cellMixture:

 
template<class ThermoType>
const ThermoType& Foam::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_;
}

What was bothering me is the += operator in this line:

 
        mixture_ += Y_[n][celli]/speciesData_[n].W()*speciesData_[n];

So we need to find out where this operator is defined. In foam-extend the enthalpy H [J/kg] = h/W (h in [J/kmol] and W in [kg/kmol]), is in terms of mass and is a member function of specie/thermo/specieThermo. So, when I look at the definition of the += member operator in specieThermo I see, that it is taken from the += member operator of janafThermo:

 
template<class equationOfState>
inline void Foam::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];
    }
}

The line

 
    const janafThermo<equationOfState>& jt

In the janafThermo += member operator definition comes from the += member operator of specieI.H:

 
inline void specie::operator+=(const specie& st)
{
    scalar sumNmoles_ = max(nMoles_ + st.nMoles_, SMALL);
 
    molWeight_ =
        nMoles_/sumNmoles_*molWeight_
      + st.nMoles_/sumNmoles_*st.molWeight_;
 
    nMoles_ = sumNmoles_
}

When I use janafThermo, for my member function h/W I need to mix (ore to treat in some way):
- W (molecular weight),
- Tlow_,
- Thigh_,
- Tcommon_,
- highCpCoeffs_
- lowCpCoeffs_
This is all I need, not more. And I need to do it on the mole fraction basis, meaning, the mole fraction X_i = nMoles_i/sumNmoles. So, the calculation of the h_mix happens via coefficients! And W_mix (for calculating H_mix) is calculated on the mole fraction basis in the += operator of specie class (molWeight_ = nMoles/sumNmoles_*….). David Gaden said, that the net effect is:

 
thermoType += scalar*thermoType

where this scalar is scalar = Y_[n][celli]/speciesData_[n].W().
And because the friend operator * of specie only works on nMoles_, I think this “scalar” does not participate in the enthalpy calculation for a mixture! Maybe it is used in some other thermoTypes, but not in mine (janaf…, idealGas… and so on).

I am really not sure about it!
So please, if someone of the advanced users could confirm my vague “theory” – that would be very kind.

Thanks in advance,

Alex

[Neka]

Small correction to the post above:
…
In foam-extend in e.g. in hPsiMixtureThermo.C that would be:

 
    forAll(hCells, celli)
    {
        hCells[celli] = this->cellMixture(celli).H(TCells[celli]);
    }

…

[Daniel_P]

Dear Alex,

I am looking into the same issue right now and also have some troubles in understanding how the mixtures are calculated exactly.

I think you are right with your statement that the calculation is done with coefficients that are averaged based on the components' mole fractions.

As far as I got it, in the code this happens here:

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 not be the same.

Is this corrected at a later stage within the mixture property calculation?

Thanks already!

[openfoammaofnepo]

Hello,

It seems that I am using the same solver as Simone's (reactingTwoPhaseEulerFoam) and thinking about the same questions. My case is gas-particle flows. gas phase is the continuous phase, while particle phase is dispersed phase. For the particle phase, it contains three differerent components of species, so I need to use multiComponentMixture model. My specification in the thermophysicalProperties.particles is as follows:

  thermoType
  {
     type            heRhoThermo;
      mixture         multiComponentMixture;
      transport       const;
      thermo          hRefConst;
      equationOfState rhoConst;
      specie          specie;
      energy          sensibleInternalEnergy;
 }

  species
  (
     char
      volatile
      ash

  );

  inertSpecie ash;

When I run this three-component particle phase, the solver always crashed. However, if I set the components as ash and char, or ash and volatile, the solver can run properly. [These two cases can run, at least indicating that all the necessary thermoType parameters for them are reasonable] Do you have the similar problems before? I do not understand why there is error (floating point exeception) when I changed to three species.

Thank you so much

[Zhiheng Wang]

Hi,
As you say it seems more Lagrangian-Eulerian flow. why dont you try reactingParcelFilmFoam or reactingParcelFoam with reactingCloud option active , As per my knowledge for char burning or oxidation you can also use FireFoam check in
$WM_PROJECT_DIR\solvers\combustion and $WM_PROJECT_DIR\solvers\Lagrangian\

[Neka]

Hello Daniel,



here is what I've found out:



1) Everything which is defined in the specie directory is dimensionless (scalar) and defined per cell.
2) The averaged thermodynamic values (Cp, h, s…) are definitely calculated on the mole weighted basis (if you calculate a gas mixture).

Actually, David Gaden’s post heats the spot.

3) In each class you have the operators defined, which are used by the members of this class. That means you can define whatever you want under e.g. plus-operator. In the janaf class under += operator the whole law of mixtures on the mole weighted basis is defined and not just += operation.
4) Moreover, the += operator in janaf includes the += operator of the equationOfState-class e.g.:

    equationOfState::operator+=(jt);

which on its part uses the definition of the += operator from the specie-class (see e.g. perfectGas class):

    specie::operator+=(pg);

So, basically you always have to take a look at operators defined in the class that you use or modify.

I’m far away of beeng a professional OpenFOAM user and much less a programmer, but here is how I see it:
Let’s stay with this code snippet:

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

, where Y[i] / species[i].W() is a scalar. The operator * in specie.H gives me as output the following thing:

s[i]*st.nMoles[i] = nMoles[i]

where s[i] = Y[i]/W[i].
Basically I get a new nMoles[i].
The operator * in janafThermo gives me as output:

Tlow[i], Thigh[i], Tcommon[i], lowCpCoeffs[i], highCpCoeffs[i]

By calling the += operator all these data will be processed. You will get molar fractions molr1/2. These fractions (variables) are used to enlarge your mixture at each for-loop by a new species. Like: specie_a = mixture, then mixture + specie_b = mixture … (all mole weighted!).

Then, this whole information goes into e.g. the H-Function(enthalpy) in e.g. Janaf, like h[mix] = RR()*(NASA-polynomial with mole averaged coefficients).

Correct me if I’m wrong.

But one thing still bothers me. The scalar is as I said s = Y[i] / species[i].W(). That means during mixing I already divide by W to get my enthalpy in terms of mass. So actually, if I calculate H[mix] I don’t need the operation H = h/W anymore. But this operation still happens in specieThermo-class (I use foam-extend). Why?

I hope I could help.

Regards,

Alex

[Daniel_P]

Hello Alex,

thank you very much for your reply and excuse my really late answer. I have been out of office for some time and did not find the time to reply earlier... :-/
Your explanations really clarified some things for me. However, I am still not sure on my initial question.

Let's start again at the code snippet:

mixture_ += Y_[n][celli]/speciesData_[n].W()*speciesData_[n];

From your post I understand that the dimensions are neglected in the specie directory.

Y[i] is the massfraction and speciesData_[i].W() the molar mass of component i. From general mixing theory, we know that the ratio of both equals molar fraction phi divided by molar mass of the mixture, i.e.: Y[i] / M[i] = phi[i] / M_mix .

Thus, we have a molar-based averaging on the right-hand side of the code snippet. However, this would mean that the variable "mixture" corresponds to a specific value whereas the variable "speciesData" is molar-based. Is this correct?

Looking for example at the specific heat capacity cp_mix of a binary mixture (components '1' and '2') this would yield:
cp_mix = 1 / M_mix * [phi_1 * cp_1 + phi_2 * cp_2] ,
with the dimensions cp_mix [J/(kg K)], M_mix [kg/mol], phi [-] and cp_1/2 [J/(mol K)].


Moreover, as far as I got it (and as you state), one has to be careful with the operator "+=" as it is defined differently for the various properties in the corresponding classes.

As an example, I looked into the files for the Sutherland coefficients As and Ts. There we have a "mixing rule" which is given by

template<class Thermo>
inline void Foam::sutherlandTransport<Thermo>::operator+=
(
    const sutherlandTransport<Thermo>& st
)
{
    scalar molr1 = this->nMoles();

    Thermo::operator+=(st);

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

    As_ = molr1*As_ + molr2*st.As_;
    Ts_ = molr1*Ts_ + molr2*st.Ts_;
}

This "mixing rule" makes use of the nMoles variable which is used frequently. However, its meaning is not completely clear to me.
Are you sure that nMoles describes the actual number of moles in a cell, i.e. "nMoles_species[i] / sumNmoles" which would be equal to the molar fraction phi?
As far as I understood it up to now, the usage of nMoles is limited to combustion models and for simple mixing processes without reactions we have nMoles=1 for all species. Is this not the case?

With respect to the operator '*' I find the following code in specieI.h (OF 2.3.x):

inline specie operator*(const scalar s, const specie& st)
{
    return specie
    (
        s*st.nMoles_,
        st.molWeight_
    );
}

I am not sure how this translates to a mathematical expression but are you sure that for the JANAF class this is not just the standard mathematical multiplication?

Where did you find the definition of the scalar s which you gave as s[i]]=Y[i]/W[i] ?


You see I am a bit lost in understanding the code.
For the simulations I am doing, the mixture property calculation is working the way it should, so I am fine from that end. Still, I would really like to understand better what is done in the code...

Thank you one more time for taking the time to write your previous post. I hope you can help me out again!

Regards,
Daniel

[Neka]

Hello Daniel,

Sorry for the delay, I was out of office and have currently a lot of things to do.

In the specie directory all things (thermodynamic properties …) are defined as scalars (thus dimensionless). You can see it if you open e.g. specieThermo.H. The dimensioning happens later, when fields are created.

By the way, I use foam-extend, so it could be that some class names differ a bit from OpenFOAM.

As far as I know the /-operator works only on scalars, thus Y[i]/W[i] is a scalar and I replaced it in my explanation with s.

Thus:

mixture_ += Y_[n][celli]/speciesData_[n].W()*speciesData_[n];

can be seen as:

mixture_ += scalar*speciesData_[n];

where mixture_ and speciesData_[n] are place holders for the operands (members).


In the second post of this thread David says, that you can see it as:

ThermoType += scalar * ThermoType

As I already mentioned I’m not a professional either and it would be nice to continue this discussion, because it can be, that I’m completely wrong. I’ll try to summarize:

We have two operators in the expression:

mixture_ += Y_[n][celli]/speciesData_[n].W()*speciesData_[n];

namely * and += (forget about /-operator and assume that Y_[n][celli]/speciesData_[n].W() = scalar “s”).

In OpenFoam the thermoType consists of several classes. An example taken from: /OpenFOAM/OpenFOAM-2.3.0/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/thermophysicalProperties

shows this:

thermoType
{
type hePsiThermo;
mixture reactingMixture;
transport sutherland;
thermo janaf;
energy sensibleEnthalpy;
equationOfState perfectGas;
specie specie;
}

That means, that in the classes hePsiThermo, reactingMixture, sutherland, janaf, sensibleEnthalpy, perfectGas and specie the operators (* and += for lhs = thermoType and rhs = scalar and thermoType) have to be defined. In other words, it is defined what must happened if those operators are called.

As you remember, I told you, that I think, that the definition of those operators is spread over several classes. What I mean is that if I calculate h_mix, than the janafThermo-class is actually responsible for the * and += operator definitions, because h is a member of janafThermo-class. But when I take a look at how they are defined, I find for the e.g. * operator the following line in janafThermo:

s*static_cast<const equationOfState>(jt),

which says me that a part of the definition of the * operator is situated in the e.g. perfectGas-class (or inherited? – I don’t know the correct expression for this process yet). And when I take a look there I find the following line:

return perfectGas(s*static_cast<const specie&>(pg));

which means that a part of the definition of the * operator is situated in the specie-class. So, at the end, we can say that the * operator for a member of the janafThermo-class is actually defined in specie-class plus some additional stuff in janafThermo-class itself.

Now let’s take a look at specie-class. The member functions of the specie-class are: W(), nMoles() and R() (see specie.H). They can be read on different ways e.g. via chemkinLexer. The mechanism for reading is at the end of specie.C file. So, we have the information about name, nMoles and molWeight for each species.
Now let’s take a look at * operator in the specieI.H:

 
inline specie operator*(const scalar s, const specie& st)
{
    return specie
    (
        s*st.nMoles_,
        st.molWeight_
    );
}

As INPUT we have here the information about the species itself (“st” has: name, nMoles and molWeight) and the scalar “s”. As OUTPUT we get a new species if we call * operator, namely s*nMoles and the molWeight of this specie. Here the private member nMoles is initialized with s*st.nMoles_ and “s” is in our case Y/W.
That means to me that if I have nMoles of specie1, than after calling * operator first time I’ll get back nMoles_specie1*(Y[1]/W[1]) and the units would be (later) [kmol/kg].

The rest of the *operator in janafThermo loads Tlow, Thigh, Tcommon, lowCpCoeffs, highCpCoeffs of this specie. So, this is how we get the necessary information of species for making a mixture. The mixture is made by the += operator in the loop. First you get specie1, than you get a mixture of specie1 and specie2 = mixture1_2. In the next loop you load specie3 and make mixture1_2 + specie3 = mixture1_2_3 and so on…

Am I right, or is there a mistake?

Let’s continue:
The definition of the += operator in janafThermo starts with:

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

So, here we first initialize molr1 (mole fraction of species1) with nMoles[species1]. Then the += operator of specie-class is done, which is:

 
inline void specie::operator+=(const specie& st)
{
    scalar sumNmoles_ = max(nMoles_ + st.nMoles_, SMALL);
 
    molWeight_ =
        nMoles_/sumNmoles_*molWeight_
      + st.nMoles_/sumNmoles_*st.molWeight_;
 
    nMoles_ = sumNmoles_;
}

Here you first calculate the sumNmoles, which is nMoles[species1] + nMoles[species2]. Or if your mixture consists of say 6 species and you are already in the middle of mixing (e.g. already looping for the 4th species) it would be sumNmoles = nMoles[mixture1_2_3] + nMoles[species4]. After calculation the sumNmoles the W_mix is calculated. And now comes an important step:

nMoles_ = sumNmoles_;

So, if we stay with our example of 6 species and have just admixed the 4th species this would be: nMoles = mixture1_2_3_4 = sumNmoles. In the next loop it will be nMoles = mixture1_2_3_4 and st.nMoles = species5 and the sumNmoles = mixture1_2_3_4_5. Am I correct?

Let’s continue with the operator +=:
…

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

Remember that nMoles is sumNmoles now! So, what we do here is we calculate mole fractions molr1 and molr2. If we would stay with our example of 6 species and would admix the 5th species molr1 would be the molar fraction of mixture1_2_3_4 and molr2 would be the molar fraction of species5. What is important is since we have the mole fractions we can calculate a molar weighted enthalpy of our mixture. And this is what is done in the remaining definition of the += operator in the janafThermo:

 
…
…
    Tlow_ = max(Tlow_, jt.Tlow_);
    Thigh_ = min(Thigh_, jt.Thigh_);
    TcommonL_ = molr1*TcommonL_ + molr2*jt.TcommonL_;
    TcommonH_ = molr1*TcommonH_ + molr2*jt.TcommonH_;
 
    for
    (
        register label coefLabel=0;
        coefLabel<janafThermo<equationOfState>::nCoeffs_;
        coefLabel++
    )
    {
        highCpCoeffs_[coefLabel] =
            molr1*highCpCoeffs_[coefLabel]
          + molr2*jt.highCpCoeffs_[coefLabel];
 
        middleCpCoeffs_[coefLabel] =
            molr1*middleCpCoeffs_[coefLabel]
          + molr2*jt.middleCpCoeffs_[coefLabel];
 
        lowCpCoeffs_[coefLabel] =
            molr1*lowCpCoeffs_[coefLabel]
          + molr2*jt.lowCpCoeffs_[coefLabel];
    }
}

If we would stay with our example of 6 species and be in the fourth loop than we would mix the janaf coefficients of mixture1_2_3_4 with the coefficients of species5. And the next loop would be the last one to get the whole mixture of 6 species.

I know, I didn’t answer all your questions, but let’s first clarify this.

This is approximately how I see it.

Please correct me when I’m wrong. I myself still couldn’t answer all questions (especially concerning units consistency). So, please keep me posted, so we can carry out this discussion till all questions are clarified.

Regards,

Alex

[Neka]

Hello Daniel,

just an addendum to my last post.

I more and more think, that the scalar „s“ is actually nMoles_. I think, that nMoles is by default = 1, and what happens during calling * operator in specie-class is nMoles = s*nMoles, thus nMoles = s*1. The dimension of s = Y/W would be [ kmol/kg], which is the number of moles of this component in one kilogram of the mixture. What do you think? That would clarify everything for me and I wouldn’t see any unit inconsistency in the mixing procedure.

Regards,

Alex

[Daniel_P]

Hello Alex,

thanks for your reply and your explanations.

I think you are right with respect to the operators, i.e. " * " or " += ". We need to look at the respective class (for example janafThermo) and then might have to go one level deeper (specie).

Moreover, I agree with your last posts. I am not entirely sure though where the assignment " nMoles=s*nMoles=s*1 " is happening.


Summing up and repeating what I hope to understand, the code is performing 4 steps starting from the following code:

mixture_ += Y_[n][celli]/speciesData_[n].W()*speciesData_[n];

### Step 1:

Operator * is taken from specieI.H:

inline specie operator*(const scalar s, const specie& st)
{
    return specie
    (
        s*st.nMoles_,
        st.molWeight_
    );
}

Here the scalar s is taken from the equation above, i.e. s=Y_[n][celli]/speciesData_[n].W() or short s=Yi/Wi.
nMoles is equal to 1 per default.

This could be where the new nMoles is assigned and the code returns nMoles=Yi/Wi. Do you agree?


### Step 2:

The operator += from specie class determines the sum of molar fractions of species1 and species2 for the first step in the mixing calculation.

As you said, first a mixture of species1 and species2 is done, then a mixture of mixedSpecies1_2 and species_3 is done. This is continued until all species are included in the mixture.

inline void specie::operator+=(const specie& st)
{
    scalar sumNmoles = max(nMoles_ + st.nMoles_, SMALL);

    molWeight_ =
        nMoles_/sumNmoles*molWeight_
      + st.nMoles_/sumNmoles*st.molWeight_;

    nMoles_ = sumNmoles;
}

Here, sumNmoles is the sum of the molar fractions of species1 and species2.
As an example, we have a mixture of gas_A, gas_B and gas_C in our cell and each has the molar fraction of 1/3 and the same molar weight W.
After the first step of mixing gas_A and gas_B, we would have something like sumNmoles=1/3+1/3=2/3 (here the molar weight is taken as 1 for all gases).


### Step 3:

The += operator in janafThermoI.H includes:

    scalar molr1 = this->nMoles();

    EquationOfState::operator+=(jt);

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

Here molr1 = Y[gas_A] / ( Y[gas_A]+Y[gas_B] ) and molr2 = Y[gas_B] / ( Y[gas_A]+Y[gas_B] ) which gives us the molar fraction of the mixture of gas_A and gas_B (not the full mixture yet).
Note that I took the molar weight equal to 1 again to simplify the writing.
In my example, molr1=molr2=1/2.


### Step 4:

Then the actual property calculation is done by molar averaging of gas_A and gas_B:

    {
        highCpCoeffs_[coefLabel] =
            molr1*highCpCoeffs_[coefLabel]
          + molr2*jt.highCpCoeffs_[coefLabel];

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

###

The process is done until all species have contributed to the mixture. In my example Step 1 to 4 would be performed twice.
1) Mixing of gas_A and gas_B --> mixture_AB
2) Mixing of mixture_AB and gas_C --> mixture_ABC


I am really not sure about these explanations but this is how far I understood it up to now from the code and your posts. Can you confirm this?


Kind regards,
Daniel

[Neka]

Hello Daniel,

sounds about right to me.

It would be nice though if someone of advanced users could confirm this.

[bbita]

Dear Foamers,

I added YEqn.H to compressibleMultiphaseInterFoam and relate that to basicMultiComponentMixture. However by changing Y, I am not observing any changes in viscosity and density.

Thanks.

[Zhiheng Wang]

Can u send detail how did u do it or atleast some tar or zip folder of solver ???

[bbita]

It is very similar to adding CEqn.H to interFoam based on Haroun et al. (2010) model.

[HosamAlrefaie]

Hi openfoammaofnepo,

can you please explain how or where to add the mass or mole fractions for these three species.

I am trying to model a gas mixture using multiComponentMixture & buoyantSimpleFoam, but I don't know how to specify the fractions.