diff --git a/Cantera/src/GasKinetics.cpp b/Cantera/src/GasKinetics.cpp index 2d5e3f94e..c258774e9 100755 --- a/Cantera/src/GasKinetics.cpp +++ b/Cantera/src/GasKinetics.cpp @@ -17,7 +17,6 @@ #include "GasKinetics.h" #include "ReactionData.h" -//#include "StoichManager.h" #include "Enhanced3BConc.h" #include "ThirdBodyMgr.h" #include "RateCoeffMgr.h" @@ -29,10 +28,10 @@ using namespace std; #include "mkl_vml.h" #endif +#ifdef HWMECH void update_kc(const double* grt, double c0, double* rkc); -void update_rates(double t, double tlog, double* rf); -void mult_by_conc(const double* c, double* ropf, double* ropr); void eval_ropnet(const double* c, const double* rf, const double* rkc, double* r); +#endif namespace Cantera { @@ -130,8 +129,6 @@ namespace Cantera { // compute Delta G^0 for all reversible reactions m_rxnstoich.getRevReactionDelta(m_ii, m_grt.begin(), m_rkc.begin()); - //m_reactantStoich.decrementReactions(m_grt.begin(), m_rkc.begin()); - //m_revProductStoich.incrementReactions(m_grt.begin(), m_rkc.begin()); doublereal logc0 = m_kdata->m_logc0; doublereal rrt = 1.0/(GasConstant * thermo().temperature()); @@ -160,12 +157,6 @@ namespace Cantera { // compute Delta G^0 for all reactions m_rxnstoich.getReactionDelta(m_ii, m_grt.begin(), rkc.begin()); - - // m_reactantStoich.decrementReactions(m_grt.begin(), rkc.begin()); - //m_revProductStoich.incrementReactions(m_grt.begin(), - //rkc.begin()); - //m_irrevProductStoich.incrementReactions(m_grt.begin(), - //rkc.begin()); doublereal logc0 = m_kdata->m_logc0; doublereal rrt = 1.0/(GasConstant * thermo().temperature()); @@ -174,6 +165,162 @@ namespace Cantera { } } + /** + * + * getDeltaGibbs(): + * + * Return the vector of values for the reaction gibbs free energy + * change + * These values depend upon the concentration + * of the ideal gas. + * + * units = J kmol-1 + */ + void GasKinetics::getDeltaGibbs(doublereal* deltaG) { + /* + * Get the chemical potentials of the species in the + * ideal gas solution. + */ + thermo().getChemPotentials(m_grt.begin()); + /* + * Use the stoichiometric manager to find deltaG for each + * reaction. + */ + m_rxnstoich.getReactionDelta(m_ii, m_grt.begin(), deltaG); + } + + /** + * + * getDeltaEnthalpy(): + * + * Return the vector of values for the reactions change in + * enthalpy. + * These values depend upon the concentration + * of the solution. + * + * units = J kmol-1 + */ + void GasKinetics::getDeltaEnthalpy(doublereal* deltaH) { + /* + * Get the partial molar enthalpy of all species in the + * ideal gas. + */ + thermo().getPartialMolarEnthalpies(m_grt.begin()); + /* + * Use the stoichiometric manager to find deltaG for each + * reaction. + */ + m_rxnstoich.getReactionDelta(m_ii, m_grt.begin(), deltaH); + } + + /************************************************************************ + * + * getDeltaEntropy(): + * + * Return the vector of values for the reactions change in + * entropy. + * These values depend upon the concentration + * of the solution. + * + * units = J kmol-1 Kelvin-1 + */ + void GasKinetics::getDeltaEntropy( doublereal* deltaS) { + /* + * Get the partial molar entropy of all species in the + * solid solution. + */ + thermo().getPartialMolarEntropies(m_grt.begin()); + /* + * Use the stoichiometric manager to find deltaS for each + * reaction. + */ + m_rxnstoich.getReactionDelta(m_ii, m_grt.begin(), deltaS); + } + + /** + * + * getDeltaSSGibbs(): + * + * Return the vector of values for the reaction + * standard state gibbs free energy change. + * These values don't depend upon the concentration + * of the solution. + * + * units = J kmol-1 + */ + void GasKinetics::getDeltaSSGibbs(doublereal* deltaG) { + /* + * Get the standard state chemical potentials of the species. + * This is the array of chemical potentials at unit activity + * We define these here as the chemical potentials of the pure + * species at the temperature and pressure of the solution. + */ + thermo().getStandardChemPotentials(m_grt.begin()); + /* + * Use the stoichiometric manager to find deltaG for each + * reaction. + */ + m_rxnstoich.getReactionDelta(m_ii, m_grt.begin(), deltaG); + } + + /** + * + * getDeltaSSEnthalpy(): + * + * Return the vector of values for the change in the + * standard state enthalpies of reaction. + * These values don't depend upon the concentration + * of the solution. + * + * units = J kmol-1 + */ + void GasKinetics::getDeltaSSEnthalpy(doublereal* deltaH) { + /* + * Get the standard state enthalpies of the species. + * This is the array of chemical potentials at unit activity + * We define these here as the enthalpies of the pure + * species at the temperature and pressure of the solution. + */ + thermo().getEnthalpy_RT(m_grt.begin()); + doublereal RT = thermo().temperature() * GasConstant; + for (int k = 0; k < m_kk; k++) { + m_grt[k] *= RT; + } + /* + * Use the stoichiometric manager to find deltaG for each + * reaction. + */ + m_rxnstoich.getReactionDelta(m_ii, m_grt.begin(), deltaH); + } + + /********************************************************************* + * + * getDeltaSSEntropy(): + * + * Return the vector of values for the change in the + * standard state entropies for each reaction. + * These values don't depend upon the concentration + * of the solution. + * + * units = J kmol-1 Kelvin-1 + */ + void GasKinetics::getDeltaSSEntropy(doublereal* deltaS) { + /* + * Get the standard state entropy of the species. + * We define these here as the entropies of the pure + * species at the temperature and pressure of the solution. + */ + thermo().getEntropy_R(m_grt.begin()); + doublereal R = GasConstant; + for (int k = 0; k < m_kk; k++) { + m_grt[k] *= R; + } + /* + * Use the stoichiometric manager to find deltaS for each + * reaction. + */ + m_rxnstoich.getReactionDelta(m_ii, m_grt.begin(), deltaS); + } void GasKinetics::processFalloffReactions() { @@ -259,6 +406,77 @@ namespace Cantera { m_kdata->m_ROP_ok = true; } + /** + * + * getFwdRateConstants(): + * + * Update the rate of progress for the reactions. + * This key routine makes sure that the rate of progress vectors + * located in the solid kinetics data class are up to date. + */ + void GasKinetics:: + getFwdRateConstants(doublereal *kfwd) { + _update_rates_T(); + _update_rates_C(); + + // copy rate coefficients into ropf + const vector_fp& rf = m_kdata->m_rfn; + array_fp& ropf = m_kdata->m_ropf; + copy(rf.begin(), rf.end(), ropf.begin()); + + // multiply ropf by enhanced 3b conc for all 3b rxns + m_3b_concm.multiply(ropf.begin(), m_kdata->concm_3b_values.begin() ); + + /* + * This routine is hardcoded to replace some of the values + * of the ropf vector. + */ + processFalloffReactions(); + + // multiply by perturbation factor + multiply_each(ropf.begin(), ropf.end(), m_perturb.begin()); + + for (int i = 0; i < m_ii; i++) { + kfwd[i] = ropf[i]; + } + } + + /** + * + * getRevRateConstants(): + * + * Return a vector of the reverse reaction rate constants + * + * Length is the number of reactions. units depends + * on many issues. Note, this routine will return rate constants + * for irreversible reactions if the default for + * doIrreversible is overridden. + */ + void GasKinetics:: + getRevRateConstants(doublereal *krev, bool doIrreversible) { + /* + * go get the forward rate constants. -> note, we don't + * really care about speed or redundancy in these + * informational routines. + */ + getFwdRateConstants(krev); + + if (doIrreversible) { + doublereal *tmpKc = m_kdata->m_ropnet.begin(); + getEquilibriumConstants(tmpKc); + for (int i = 0; i < m_ii; i++) { + krev[i] /= tmpKc[i]; + } + } else { + /* + * m_rkc[] is zero for irreversibly reactions + */ + const vector_fp& m_rkc = m_kdata->m_rkcn; + for (int i = 0; i < m_ii; i++) { + krev[i] *= m_rkc[i]; + } + } + } void GasKinetics:: addReaction(const ReactionData& r) { diff --git a/Cantera/src/GasKinetics.h b/Cantera/src/GasKinetics.h index da54dd822..0bc46d839 100755 --- a/Cantera/src/GasKinetics.h +++ b/Cantera/src/GasKinetics.h @@ -75,7 +75,10 @@ namespace Cantera { class GasKinetics : public Kinetics { public: - + /** + * @name Constructors and General Information about Mechanism + */ + //@{ /// Constructor. GasKinetics(thermo_t* thermo = 0); @@ -92,21 +95,123 @@ namespace Cantera { return m_prxn[k][i]; } + //@} + /** + * @name Reaction Rates Of Progress + */ + //@{ + /** + * Forward rates of progress. + * Return the forward rates of progress in array fwdROP, which + * must be dimensioned at least as large as the total number + * of reactions. + */ virtual void getFwdRatesOfProgress(doublereal* fwdROP) { updateROP(); copy(m_kdata->m_ropf.begin(), m_kdata->m_ropf.end(), fwdROP); } - + /** + * Reverse rates of progress. + * Return the reverse rates of progress in array revROP, which + * must be dimensioned at least as large as the total number + * of reactions. + */ virtual void getRevRatesOfProgress(doublereal* revROP) { updateROP(); copy(m_kdata->m_ropr.begin(), m_kdata->m_ropr.end(), revROP); } - + /** + * Net rates of progress. Return the net (forward - reverse) + * rates of progress in array netROP, which must be + * dimensioned at least as large as the total number of + * reactions. + */ virtual void getNetRatesOfProgress(doublereal* netROP) { updateROP(); copy(m_kdata->m_ropnet.begin(), m_kdata->m_ropnet.end(), netROP); } + + /** + * Equilibrium constants. Return the equilibrium constants of + * the reactions in concentration units in array kc, which + * must be dimensioned at least as large as the total number + * of reactions. + */ + virtual void getEquilibriumConstants(doublereal* kc); + + /** + * Return the vector of values for the reaction gibbs free energy + * change. + * These values depend upon the concentration + * of the solution. + * + * units = J kmol-1 + */ + virtual void getDeltaGibbs( doublereal* deltaG); + + /** + * Return the vector of values for the reactions change in + * enthalpy. + * These values depend upon the concentration + * of the solution. + * + * units = J kmol-1 + */ + virtual void getDeltaEnthalpy( doublereal* deltaH); + + /** + * Return the vector of values for the reactions change in + * entropy. + * These values depend upon the concentration + * of the solution. + * + * units = J kmol-1 Kelvin-1 + */ + virtual void getDeltaEntropy(doublereal* deltaS); + + /** + * Return the vector of values for the reaction + * standard state gibbs free energy change. + * These values don't depend upon the concentration + * of the solution. + * + * units = J kmol-1 + */ + virtual void getDeltaSSGibbs(doublereal* deltaG); + + /** + * Return the vector of values for the change in the + * standard state enthalpies of reaction. + * These values don't depend upon the concentration + * of the solution. + * + * units = J kmol-1 + */ + virtual void getDeltaSSEnthalpy(doublereal* deltaH); + + /** + * Return the vector of values for the change in the + * standard state entropies for each reaction. + * These values don't depend upon the concentration + * of the solution. + * + * units = J kmol-1 Kelvin-1 + */ + virtual void getDeltaSSEntropy(doublereal* deltaS); + + //@} + /** + * @name Species Production Rates + */ + //@{ + + /** + * Species net production rates [kmol/m^3]. Return the species + * net production rates (creation - destruction) in array + * wdot, which must be dimensioned at least as large as the + * total number of species. + */ virtual void getNetProductionRates(doublereal* net) { updateROP(); #ifdef HWMECH @@ -123,6 +228,13 @@ namespace Cantera { #endif } + /** + * Species creation rates [kmol/m^3]. Return the species + * creation rates in array cdot, which must be + * dimensioned at least as large as the total number of + * species. + * + */ virtual void getCreationRates(doublereal* cdot) { updateROP(); m_rxnstoich.getCreationRates(m_kk, m_kdata->m_ropf.begin(), @@ -136,6 +248,13 @@ namespace Cantera { // m_kdata->m_ropr.begin(), cdot); } + /** + * Species destruction rates [kmol/m^3]. Return the species + * destruction rates in array ddot, which must be + * dimensioned at least as large as the total number of + * species. + * + */ virtual void getDestructionRates(doublereal* ddot) { updateROP(); m_rxnstoich.getDestructionRates(m_kk, m_kdata->m_ropf.begin(), @@ -147,9 +266,63 @@ namespace Cantera { // m_kdata->m_ropf.begin(), ddot); } - virtual void getEquilibriumConstants(doublereal* kc); + //@} + /** + * @name Reaction Mechanism Informational Query Routines + */ + //@{ /** + * Flag specifying the type of reaction. The legal values and + * their meaning are specific to the particular kinetics + * manager. + */ + virtual int reactionType(int i) const { + return m_index[i].first; + } + + virtual string reactionString(int i) const { + return m_rxneqn[i]; + } + + /** + * True if reaction i has been declared to be reversible. If + * isReversible(i) is false, then the reverse rate of progress + * for reaction i is always zero. + */ + virtual bool isReversible(int i) { + if (find(m_revindex.begin(), m_revindex.end(), i) + < m_revindex.end()) return true; + else return false; + } + + /** + * Return the forward rate constants + * + * length is the number of reactions. units depends + * on many issues. + */ + virtual void getFwdRateConstants(doublereal *kfwd); + + /** + * Return the reverse rate constants. + * + * length is the number of reactions. units depends + * on many issues. Note, this routine will return rate constants + * for irreversible reactions if the default for + * doIrreversible is overridden. + */ + virtual void getRevRateConstants(doublereal *krev, + bool doIrreversible = false); + + //@} + /** + * @name Reaction Mechanism Setup Routines + */ + //@{ + + + /** * Set delta T threshold for updating temperature-dependent * rates. */ @@ -170,28 +343,18 @@ namespace Cantera { void updateROP(); - virtual int reactionType(int i) const { - return m_index[i].first; - } - - virtual string reactionString(int i) const { - return m_rxneqn[i]; - } const vector& reactantGroups(int i) { return m_rgroups[i]; } const vector& productGroups(int i) { return m_pgroups[i]; } - virtual bool isReversible(int i) { - if (find(m_revindex.begin(), m_revindex.end(), i) - < m_revindex.end()) return true; - else return false; - } void _update_rates_T(); void _update_rates_C(); + //@} + protected: int m_kk, m_nfall;