[mAlletto]

Hello all,


I'm trying to understand how the reverse rate is calculated. Following the manual of chemkin (https://www3.nd.edu/~powers/ame.60636/chemkin2000.pdf) the reverse rate of the ith reaction is calculated as follows:





the equilibrium constant is calculated as follows:





being the sum of the stoichiometric coefficient of the ith reaction. Furthermore






being are the enthalpy and entropy of the species k.


the source code for the calculation of k_ri can be found in $FOAM_SRC/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C:

template
<
    template<class> class ReactionType,
    class ReactionThermo,
    class ReactionRate
>
Foam::scalar Foam::ReversibleReaction
<
    ReactionType,
    ReactionThermo,
    ReactionRate
>::kr
(
    const scalar kfwd,
    const scalar p,
    const scalar T,
    const scalarField& c
) const
{
    return kfwd/max(this->Kc(p, T), 1e-6);
}

and the calculation of K_ci can be found in $FOAM_SRC/thermophysicalModels/specie/thermo/thermo/thermoI.H:

template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::K(const scalar p, const scalar T) const
{
    scalar arg = -this->Y()*this->G(Pstd, T)/(RR*T);

    if (arg < 600)
    {
        return exp(arg);
    }
    else
    {
        return VGREAT;
    }
}


template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::Kp(const scalar p, const scalar T) const
{
    return K(p, T);
}


template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::Kc(const scalar p, const scalar T) const
{
    const scalar nm = this->Y()/this->W();

    if (equal(nm, SMALL))
    {
        return Kp(p, T);
    }
    else
    {
        return Kp(p, T)*pow(Pstd/(RR*T), nm);
    }
}

So what I do not understand is how the stoichiometric coefficients enter into the source code in openfoam.





Best,


Michael

[adhiraj]

As far as I can tell, OpenFOAM handles multiple species like this:

• A pseudo-species is created by taking a weighted average of all the species present in the system.
• For a mixture, the weights for each species is proportional to the actual fraction (mass or mole, I don't recall exactly) of that species in the mixture.
• For a reaction, the difference in stoichiometric coefficients for each species is used.
• For all calculations the pseudo-species is used, so stoichiometric coefficients do not show up directly in the functions you see here.

[mAlletto]

Thank you Adhiraj for pointing me in the right direction. Looking at the constructor of the Reaction class ($FOAM_SRC/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C)


we see the function setThermo

template<class ReactionThermo>
Foam::Reaction<ReactionThermo>::Reaction
(
    const speciesTable& species,
    const List<specieCoeffs>& lhs,
    const List<specieCoeffs>& rhs,
    const HashPtrTable<ReactionThermo>& thermoDatabase,
    bool initReactionThermo
)
:
    ReactionThermo::thermoType(*thermoDatabase[species[0]]),
    name_("un-named-reaction-" + Foam::name(getNewReactionID())),
    species_(species),
    lhs_(lhs),
    rhs_(rhs)
{
    if (initReactionThermo)
    {
        setThermo(thermoDatabase);
    }
}

in this function we see that the stoichiometric coefficients and the molar weight of the species enter in the construction of the thermodynamic properties:

template<class ReactionThermo>
void Foam::Reaction<ReactionThermo>::setThermo
(
    const HashPtrTable<ReactionThermo>& thermoDatabase
)
{

    typename ReactionThermo::thermoType rhsThermo
    (
        rhs_[0].stoichCoeff
        *(*thermoDatabase[species_[rhs_[0].index]]).W()
        *(*thermoDatabase[species_[rhs_[0].index]])
    );

    for (label i=1; i<rhs_.size(); ++i)
    {
        rhsThermo +=
            rhs_[i].stoichCoeff
        *(*thermoDatabase[species_[rhs_[i].index]]).W()
        *(*thermoDatabase[species_[rhs_[i].index]]);
    }

    typename ReactionThermo::thermoType lhsThermo
    (
        lhs_[0].stoichCoeff
       *(*thermoDatabase[species_[lhs_[0].index]]).W()
       *(*thermoDatabase[species_[lhs_[0].index]])
    );

    for (label i=1; i<lhs_.size(); ++i)
    {
        lhsThermo +=
            lhs_[i].stoichCoeff
           *(*thermoDatabase[species_[lhs_[i].index]]).W()
           *(*thermoDatabase[species_[lhs_[i].index]]);
    }

    ReactionThermo::thermoType::operator=(lhsThermo == rhsThermo);
}

[mAlletto]

I had a look how the exponent in is calculated.



Key feature is the function setThermo in the file Reaction.C:

   



template<class ReactionThermo>
void Foam::Reaction<ReactionThermo>::setThermo
(
    const HashPtrTable<ReactionThermo>& thermoDatabase
)
{

    typename ReactionThermo::thermoType rhsThermo
    (
        rhs_[0].stoichCoeff
        *(*thermoDatabase[species_[rhs_[0].index]]).W()
        *(*thermoDatabase[species_[rhs_[0].index]])
    );

    for (label i=1; i<rhs_.size(); ++i)
    {
        rhsThermo +=
            rhs_[i].stoichCoeff
        *(*thermoDatabase[species_[rhs_[i].index]]).W()
        *(*thermoDatabase[species_[rhs_[i].index]]);
    }

    typename ReactionThermo::thermoType lhsThermo
    (
        lhs_[0].stoichCoeff
       *(*thermoDatabase[species_[lhs_[0].index]]).W()
       *(*thermoDatabase[species_[lhs_[0].index]])
    );

    for (label i=1; i<lhs_.size(); ++i)
    {
        lhsThermo +=
            lhs_[i].stoichCoeff
           *(*thermoDatabase[species_[lhs_[i].index]]).W()
           *(*thermoDatabase[species_[lhs_[i].index]]);
    }

    ReactionThermo::thermoType::operator=(lhsThermo == rhsThermo);
}

What is done it is constructed a thermodynamic object by multiplying the thermodynamic object of each species on the right hand side with the stoichiometric coefficient and the molar weight and summing them up. The same is done for the left hand side. Afterward the == operator is applied.



To calculate the exponent we need to find out how this operation affect the Y_ and molWeight_ data of the of the thermodynamic object is affected by this operation.



Since the thermodynamic object inherits from the equation of state object and this again inherits from the species object and the latter carries the data Y_ and molWeight_ we'll have a look how the above operators are defined in $FOAM_SRC/thermophysicalModels/specie/specie/specieI.H:

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


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;
}


inline specie operator==(const specie& st1, const specie& st2)
{
    scalar diffY = st2.Y_ - st1.Y_;
    if (mag(diffY) < SMALL)
    {
        diffY = SMALL;
    }

    const scalar diffRW = st2.Y_/st2.molWeight_ - st1.Y_/st1.molWeight_;

    #ifdef __clang__
    // Using intermediate volatile bool to prevent compiler optimising out the
    // if block (above) - CLANG 3.7.1
    volatile const bool valid = (mag(diffRW) > SMALL);
    const scalar molWeight = valid ? diffY/diffRW : GREAT;
    #else
    scalar molWeight = GREAT;
    if (mag(diffRW) > SMALL)
    {
         molWeight = diffY/diffRW;
    }
    #endif

    return specie(diffY, molWeight);
}

For the multiplication of a species object with a scalar with see that only Y_ is multiplied by the scalar. For the operator += the two Y_ objects are summed. The molWeight_ object is calculated as follows:

molWeight_ = sumY/(Y_/molWeight_ + st.Y_/st.molWeight_);

So after the first two multiplications the data are





and






after the summation we obtain hence





and





for the left hand side (the forward side of the reaction). An analogous formulation is obtained also for the right hand side.



After the application of the == operator we obtain:




and





We see that the exponent nm is computed as follows:

 const scalar nm = this->Y()/this->W();

Be dividing Y_ by molWeight_ we get the same formulation as in the chemkin manual:

[mAlletto]

If someone is interest there is a description of how the reverse rate constant is calculated also here: https://openfoamwiki.net/index.php/ChemFoam

[labyrinth01]

Hi mAlletto, then how should we explain the calculation of arg = -this->Y()*this->G(Pstd, T)/(RR*T) in OpenFoam instead of -this->G(Pstd, T)/(RR*T)?

Thanks.

[labyrinth01]

It is still not clear why Y=v*W