CFD Online Discussion Forums - understanding reverse rate calculation for chemical reactions

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM (https://www.cfd-online.com/Forums/openfoam/)

- - understanding reverse rate calculation for chemical reactions (https://www.cfd-online.com/Forums/openfoam/224970-understanding-reverse-rate-calculation-chemical-reactions.html)

understanding reverse rate calculation for chemical reactions

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:

$k_{ri} = \frac{k_{fi}}{K_{ci}}$

the equilibrium constant is calculated as follows:

${K_{ci}} = K_{pi} (\frac{p_{atm}}{R T})^{\sum_k \nu_{ki}}$

being $\nu_{ki}$ the sum of the stoichiometric coefficient of the ith reaction. Furthermore

$K_{pi} = exp ( \sum_k \nu_{ki} \frac{S^0_k}{R} - \sum_k \nu_{ki} \frac{H^0_k}{RT})$

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:

Code:

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:

Code:

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 $\nu_{ki}$ enter into the source code in openfoam.

Best,

Michael

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 $\nu_{ki}$ do not show up directly in the functions you see here.

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

Code:

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:

Code:

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

}

I had a look how the exponent in $K_{ci}= K_{pi}(\frac{p_{atm}}{RT})^{\sum_k \nu_{ki} }$ is calculated.

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

Code:

   







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:

Code:

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:

Code:

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

So after the first two multiplications the data are

$Y_ = \nu W$

and

after the summation we obtain hence

$Y' = \sum_k \nu_k' W_k'$

and

$molWeight' = \frac{ \sum_k \nu_k' W_k'}{ \sum_k \nu_k' }$

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:

$Y_ = \sum_k \nu_k' W_k' - \sum_k \nu_k'' W_k''$

and

$molWeight_ = \frac{ \sum_k \nu_k' W_k' - \sum_k \nu_k'' W_k'}{ \sum_k \nu_k' - \sum_k \nu_k'' }$

We see that the exponent nm is computed as follows:

Code:

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

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

$nm = \sum_k \nu_k' - \sum_k \nu_k''$

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

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.

It is still not clear why Y=v*W