Understanding reaction rates

vbchris · February 22, 2018, 15:02

Hello all,

I'm trying to implement a new reaction rate type in OF 5.0. I am having some trouble understanding the code and if anyone can provide clarity it would be appreciated.

In reactingFoam the species conversion is implemented in the YEqn.h file:

Code:

            fvScalarMatrix YiEqn
            (
                fvm::ddt(rho, Yi)
              + mvConvection->fvmDiv(phi, Yi)
              - fvm::laplacian(turbulence->muEff(), Yi)
             ==
                reaction->R(Yi)
              + fvOptions(rho, Yi)
            );

The function R I believe references the chemistryModel member RR which is calculated

Code:

 // Calculates the reaction rate
  // RR_[i] can then be called by chemistry.rr() in the YEqn.h file by the solver
  
  
  
  644 template<class CompType, class ThermoType>
  645 void Foam::chemistryModel<CompType, ThermoType>::calculate()
  646 {
  647     if (!this->chemistry_)
  648     {
  649         return;
  650     }
  651 
  
  // Gets the necessary property fields from thermo().  This is where all thermo properties are handled.
  
  652     tmp<volScalarField> trho(this->thermo().rho());
  653     const scalarField& rho = trho();
  654 
  655     const scalarField& T = this->thermo().T();
  656     const scalarField& p = this->thermo().p();
  657 
  
  // Main loop for calculating the reaction rates
  // The index celli refeers to a specific cell in the mesh
  
  658     forAll(rho, celli)
  659     {  
  660         const scalar rhoi = rho[celli];
  661         const scalar Ti = T[celli];
  662         const scalar pi = p[celli];
  663 
  
  // The index i referrers to a species int 
  
  664         for (label i=0; i<nSpecie_; i++)
  665         {
  666             const scalar Yi = Y_[i][celli];
  667             c_[i] = rhoi*Yi/specieThermo_[i].W();
  668         }
  669 
  
  // Omega modifies dcdt_, Ti and pi are copies and c_ isn't used in the following loop.
  
  670         omega(c_, Ti, pi, dcdt_);
  671 
  672         for (label i=0; i<nSpecie_; i++)
  673         {
  674             RR_[i][celli] = dcdt_[i]*specieThermo_[i].W();
  675         }
  676     }
  677 }

So RR = dcdt*MolecularWeight

dcdt is set by the omega function which is virtual

Code:

template<class CompType, class ThermoType>
   97 void Foam::chemistryModel<CompType, ThermoType>::omega
   98 (
   99     const scalarField& c,
  100     const scalar T,
  101     const scalar p,
  102     scalarField& dcdt
  103 ) const
  104 {
  105     scalar pf, cf, pr, cr;
  106     label lRef, rRef;
  107 
  
  //initializes dcdt to 0
  
  108     dcdt = Zero;
  109 
  
  // For each rxn
  
  110     forAll(reactions_, i)
  111     {
  112         const Reaction<ThermoType>& R = reactions_[i];
  113 
  114         scalar omegai = omega
  115         (
  116             R, c, T, p, pf, cf, lRef, pr, cr, rRef
  117         );
  118 
  119         forAll(R.lhs(), s)
  120         {
  121             const label si = R.lhs()[s].index;
  122             const scalar sl = R.lhs()[s].stoichCoeff;
  
  // dcdt = stoichCoeff*
  
  123             dcdt[si] -= sl*omegai;
  124         }
  125 
  126         forAll(R.rhs(), s)
  127         {
  128             const label si = R.rhs()[s].index;
  129             const scalar sr = R.rhs()[s].stoichCoeff;
  130             dcdt[si] += sr*omegai;
  131         }
  132     }
  133 }
  

template<class CompType, class ThermoType>
   97 void Foam::chemistryModel<CompType, ThermoType>::omega
   98 (
   99     const scalarField& c,
  100     const scalar T,
  101     const scalar p,
  102     scalarField& dcdt
  103 ) const
  104 {
  105     scalar pf, cf, pr, cr;
  106     label lRef, rRef;
  107 
  
  //initializes dcdt to 0
  
  108     dcdt = Zero;
  109 
  
  // For each rxn
  
  110     forAll(reactions_, i)
  111     {
  112         const Reaction<ThermoType>& R = reactions_[i];
  113 
  114         scalar omegai = omega
  115         (
  116             R, c, T, p, pf, cf, lRef, pr, cr, rRef
  117         );
  118 
  119         forAll(R.lhs(), s)
  120         {
  121             const label si = R.lhs()[s].index;
  122             const scalar sl = R.lhs()[s].stoichCoeff;
  
  // dcdt = stoichCoeff*
  
  123             dcdt[si] -= sl*omegai;
  124         }
  125 
  126         forAll(R.rhs(), s)
  127         {
  128             const label si = R.rhs()[s].index;
  129             const scalar sr = R.rhs()[s].stoichCoeff;
  130             dcdt[si] += sr*omegai;
  131         }
  132     }
  133 }

Here dcdt = sigma(+-stoichCoeff*omegai) where the lhs/rhs of the reaction defines the sign of the stoichCoeff

However I'm struggleing to understand the code for omegai

Code:

157 template<class CompType, class ThermoType>
  158 Foam::scalar Foam::chemistryModel<CompType, ThermoType>::omega
  159 (
  160     const Reaction<ThermoType>& R,
  161     const scalarField& c,
  162     const scalar T,
  163     const scalar p,
  164     scalar& pf,
  165     scalar& cf,
  166     label& lRef,
  167     scalar& pr,
  168     scalar& cr,
  169     label& rRef
  170 ) const
  171 {
  172     const scalar kf = R.kf(p, T, c);
  173     const scalar kr = R.kr(kf, p, T, c);
  174 
  175     pf = 1.0;
  176     pr = 1.0;
  177 
  178     const label Nl = R.lhs().size();
  179     const label Nr = R.rhs().size();
  180 
  181     label slRef = 0;
  182     lRef = R.lhs()[slRef].index;
  183 
  184     pf = kf;
  185     for (label s = 1; s < Nl; s++)
  186     {
  187         const label si = R.lhs()[s].index;
  188 
  189         if (c[si] < c[lRef])
  190         {
  191             const scalar exp = R.lhs()[slRef].exponent;
  192             pf *= pow(max(0.0, c[lRef]), exp);
  193             lRef = si;
  194             slRef = s;
  195         }
  196         else
  197         {
  198             const scalar exp = R.lhs()[s].exponent;
  199             pf *= pow(max(0.0, c[si]), exp);
  200         }
  201     }
  202     cf = max(0.0, c[lRef]);
  203 
  204     {
  205         const scalar exp = R.lhs()[slRef].exponent;
  206         if (exp < 1.0)
  207         {
  208             if (cf > SMALL)
  209             {
  210                 pf *= pow(cf, exp - 1.0);
  211             }
  212             else
  213             {
  214                 pf = 0.0;
  215             }
  216         }
  217         else
  218         {
  219             pf *= pow(cf, exp - 1.0);
  220         }
  221     }
  222 
  223     label srRef = 0;
  224     rRef = R.rhs()[srRef].index;
  225 
  226     // Find the matrix element and element position for the rhs
  227     pr = kr;
  228     for (label s = 1; s < Nr; s++)
  229     {
  230         const label si = R.rhs()[s].index;
  231         if (c[si] < c[rRef])
  232         {
  233             const scalar exp = R.rhs()[srRef].exponent;
  234             pr *= pow(max(0.0, c[rRef]), exp);
  235             rRef = si;
  236             srRef = s;
  237         }
  238         else
  239         {
  240             const scalar exp = R.rhs()[s].exponent;
  241             pr *= pow(max(0.0, c[si]), exp);
  242         }
  243     }
  244     cr = max(0.0, c[rRef]);
  245 
  246     {
  247         const scalar exp = R.rhs()[srRef].exponent;
  248         if (exp < 1.0)
  249         {
  250             if (cr>SMALL)
  251             {
  252                 pr *= pow(cr, exp - 1.0);
  253             }
  254             else
  255             {
  256                 pr = 0.0;
  257             }
  258         }
  259         else
  260         {
  261             pr *= pow(cr, exp - 1.0);
  262         }
  263     }
  264   
  265     return pf*cf - pr*cr;
  266 }

Can anyone help explain how the reaction rate is computed from Omegai?

Also, I see that the Arrhenius reaction rate is calculated in ArrheniusReactionRateI.H but I dont see how this is implemented in omegai?

Code:

inline Foam::scalar Foam::ArrheniusReactionRate::operator()
(
    const scalar p,
    const scalar T,
    const scalarField&
) const
{
    scalar ak = A_;

    if (mag(beta_) > VSMALL)
    {
        ak *= pow(T, beta_);
    }

    if (mag(Ta_) > VSMALL)
    {
        ak *= exp(-Ta_/T);
    }

    return ak;
}

Thanks!

nerv383 · February 26, 2018, 01:52

Hi,

I am having the same question. I suppose the code for omegai is to calculate the progress rate for one reaction as the expression below

$\Omega$ = $k_{f}$ * $\prod^{N}_{k=1}$ $(c_{kf})^{v_{kf}}$ - $k_{r}$ * $\prod^{N}_{k=1}$ $(c_{kr})^{v_{kr}}$

However, if that is the case, I cannot understand how the code achieve this calculation, especially for the snippet

Code:

        if (exp < 1.0)
        {
            if (cf > SMALL)
            {
                pf *= pow(cf, exp - 1.0);
            }
            else
            {
                pf = 0.0;
            }
        }
        else
        {
            pf *= pow(cf, exp - 1.0);
        }
    }

It would be very appreciated if anyone can share any ideas about it.

From my understanding, the Arrhenius law is applied when kf and kr is calculated.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Modelling species transport and implementing reaction rates	res	FLUENT	0	May 11, 2017 10:29
Segmentation fault in running alternateSteadyReactingFoam,why?	NewKid	OpenFOAM	18	January 20, 2011 16:55
combustion - reaction rates - urgent help needed	siri	Main CFD Forum	2	March 3, 2007 12:25
chemical reaction - decompostition	La S. Hyuck	CFX	1	May 23, 2001 00:07
Reaction rates	Szasz Robert	FLUENT	1	May 10, 2000 06:06

February 26, 2018, 01:52		#2
nerv383 New Member zhijie huo Join Date: Feb 2017 Posts: 1 Rep Power: 0	Hi, I am having the same question. I suppose the code for omegai is to calculate the progress rate for one reaction as the expression below $\Omega$ = $k_{f}$ * $\prod^{N}_{k=1}$ $(c_{kf})^{v_{kf}}$ - $k_{r}$ * $\prod^{N}_{k=1}$ $(c_{kr})^{v_{kr}}$ However, if that is the case, I cannot understand how the code achieve this calculation, especially for the snippet Code: if (exp < 1.0) { if (cf > SMALL) { pf = pow(cf, exp - 1.0); } else { pf = 0.0; } } else { pf = pow(cf, exp - 1.0); } } It would be very appreciated if anyone can share any ideas about it. From my understanding, the Arrhenius law is applied when kf and kr is calculated.