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

understanding reverse rate calculation for chemical reactions

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

Like Tree2Likes
  • 1 Post By mAlletto
  • 1 Post By adhiraj

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   March 10, 2020, 07:41
Default understanding reverse rate calculation for chemical reactions
  #1
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 335
Rep Power: 6
mAlletto is on a distinguished road
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 H^0_k, S^0_k 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
Tobi likes this.
mAlletto is offline   Reply With Quote

Old   March 10, 2020, 16:23
Default
  #2
Senior Member
 
Adhiraj
Join Date: Sep 2010
Location: Maryland, United States
Posts: 162
Rep Power: 11
adhiraj is on a distinguished road
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.
mAlletto likes this.
adhiraj is offline   Reply With Quote

Old   March 11, 2020, 17:25
Default
  #3
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 335
Rep Power: 6
mAlletto is on a distinguished road
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);
}
mAlletto is offline   Reply With Quote

Old   March 12, 2020, 11:02
Default
  #4
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 335
Rep Power: 6
mAlletto is on a distinguished road
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



molWeight_ =  W


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''
mAlletto is offline   Reply With Quote

Old   March 16, 2020, 11:52
Default
  #5
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Location: Bremen
Posts: 335
Rep Power: 6
mAlletto is on a distinguished road
If someone is interest there is a description of how the reverse rate constant is calculated also here: https://openfoamwiki.net/index.php/ChemFoam
mAlletto 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
calculation of mass flow rate from the inner part of a vertical cylinder in natural manoj kumar dash FLUENT 5 August 24, 2018 10:57
Understanding Surface Growth Rate alw7 STAR-CCM+ 2 August 13, 2018 19:53
PaSR + infinite reaction rate in reactingFoam --> no reactions occurring tatu OpenFOAM Running, Solving & CFD 2 November 24, 2016 18:34
Reverse rate parameters with CHEMKIN imports Kevin L FLUENT 1 June 27, 2013 04:22
Heat transfer rate calculation Lotta_questions CFX 3 February 11, 2010 18:35


All times are GMT -4. The time now is 10:30.