diff --git a/Cantera/src/kinetics/InterfaceKinetics.cpp b/Cantera/src/kinetics/InterfaceKinetics.cpp index 090b8585e..2f5cd0cd7 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.cpp +++ b/Cantera/src/kinetics/InterfaceKinetics.cpp @@ -25,9 +25,8 @@ using namespace std; namespace Cantera { - ////////////////////////////////////////////////////////////////// - - /** + //==================================================================================================================== + /* * Construct an empty InterfaceKinetics reaction mechanism. * @param thermo This is an optional parameter that may be * used to initialize the inherited Kinetics class with @@ -35,8 +34,7 @@ namespace Cantera { * useful for initialization of homogeneous kinetics * mechanisms. */ - InterfaceKinetics:: - InterfaceKinetics(thermo_t* thermo) : + InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) : Kinetics(), m_kk(0), m_redo_rates(false), @@ -44,27 +42,41 @@ namespace Cantera { m_nrev(0), m_surf(0), m_integrator(0), + m_beta(0), + m_ctrxn(0), + m_ctrxn_ecdf(0), + m_StandardConc(0), + m_deltaG0(0), + m_ProdStanConcReac(0), m_finalized(false), m_has_coverage_dependence(false), m_has_electrochem_rxns(false), + m_has_exchange_current_density_formulation(false), + m_phaseExistsCheck(false), + m_phaseExists(0), + m_rxnPhaseIsReactant(0), + m_rxnPhaseIsProduct(0), m_ioFlag(0) { if (thermo != 0) addPhase(*thermo); m_kdata = new InterfaceKineticsData; m_kdata->m_temp = 0.0; } - - /** + //==================================================================================================================== + /* * Destructor */ - InterfaceKinetics:: - ~InterfaceKinetics(){ + InterfaceKinetics::~InterfaceKinetics(){ delete m_kdata; if (m_integrator) { delete m_integrator; } + for (int i = 0; i < m_ii; i++) { + delete m_rxnPhaseIsReactant[i]; + delete m_rxnPhaseIsProduct[i]; + } } - + //==================================================================================================================== // Copy Constructor for the %InterfaceKinetics object. /* * Currently, this is not fully implemented. If called it will @@ -78,9 +90,20 @@ namespace Cantera { m_nrev(0), m_surf(0), m_integrator(0), + m_beta(0), + m_ctrxn(0), + m_ctrxn_ecdf(0), + m_StandardConc(0), + m_deltaG0(0), + m_ProdStanConcReac(0), m_finalized(false), m_has_coverage_dependence(false), m_has_electrochem_rxns(false), + m_has_exchange_current_density_formulation(false), + m_phaseExistsCheck(false), + m_phaseExists(0), + m_rxnPhaseIsReactant(0), + m_rxnPhaseIsProduct(0), m_ioFlag(0) { m_kdata = new InterfaceKineticsData; @@ -90,7 +113,7 @@ namespace Cantera { */ *this = operator=(right); } - + //==================================================================================================================== // Assignment operator /* * This is NOT a virtual function. @@ -100,13 +123,21 @@ namespace Cantera { */ InterfaceKinetics& InterfaceKinetics:: operator=(const InterfaceKinetics &right) { + int i; /* * Check for self assignment. */ if (this == &right) return *this; + for (i = 0; i < m_ii; i++) { + delete (m_rxnPhaseIsReactant[i]); + delete (m_rxnPhaseIsProduct[i]); + } + Kinetics::operator=(right); + + m_kk = right.m_kk; m_revindex = right.m_revindex; m_rates = right.m_rates; @@ -129,17 +160,58 @@ namespace Cantera { m_integrator = right.m_integrator; //DANGER - shallow copy m_beta = right.m_beta; m_ctrxn = right.m_ctrxn; + m_ctrxn_ecdf = right.m_ctrxn_ecdf; + m_StandardConc = right.m_StandardConc; + m_deltaG0 = right.m_deltaG0; + m_ProdStanConcReac = right.m_ProdStanConcReac; m_finalized = right.m_finalized; m_has_coverage_dependence = right.m_has_coverage_dependence; m_has_electrochem_rxns = right.m_has_electrochem_rxns; + m_has_exchange_current_density_formulation = right.m_has_exchange_current_density_formulation; + m_phaseExistsCheck = right.m_phaseExistsCheck; + m_phaseExists = right.m_phaseExists; + + + m_rxnPhaseIsReactant.resize(m_ii, 0); + m_rxnPhaseIsProduct.resize(m_ii, 0); + int np = nPhases(); + for (i = 0; i < m_ii; i++) { + m_rxnPhaseIsReactant[i] = new bool[np]; + m_rxnPhaseIsProduct[i] = new bool[np]; + for (int p = 0; p < np; p++) { + m_rxnPhaseIsReactant[i][p] = right.m_rxnPhaseIsReactant[i][p]; + m_rxnPhaseIsProduct[i][p] = right.m_rxnPhaseIsProduct[i][p]; + } + } + + m_rxnPhaseIsProduct = right.m_rxnPhaseIsProduct; + m_ioFlag = right.m_ioFlag; return *this; } + //==================================================================================================================== + // Return the ID of the kinetics object + int InterfaceKinetics::ID() const { + return cInterfaceKinetics; + } + //==================================================================================================================== + int InterfaceKinetics::type() const { + return cInterfaceKinetics; + } + //==================================================================================================================== + // Set the electric potential in the nth phase + /* + * @param n phase Index in this kinetics object. + * @param V Electric potential (volts) + */ + void InterfaceKinetics::setElectricPotential(int n, doublereal V) { + thermo(n).setElectricPotential(V); + m_redo_rates = true; + } - - // Duplication routine for objects which inherit from - // Kinetics + //==================================================================================================================== + // Duplication routine for objects which inherit from Kinetics /* * This virtual routine can be used to duplicate %Kinetics objects * inherited from %Kinetics even if the application only has @@ -152,14 +224,18 @@ namespace Cantera { InterfaceKinetics* tp = new InterfaceKinetics(*this); return dynamic_cast(tp); } - - - /** - * Update properties that depend on temperature + //==================================================================================================================== + // Update properties that depend on temperature + /* + * This is called to update all of the properties that depend on temperature * - */ - void InterfaceKinetics:: - _update_rates_T() { + * Current objects that this function updates + * m_kdata->m_logtemp + * m_kdata->m_rfn + * m_rates. + * updateKc(); + */ + void InterfaceKinetics::_update_rates_T() { _update_rates_phi(); if (m_has_coverage_dependence) { m_surf->getCoverages(DATA_PTR(m_conc)); @@ -167,20 +243,24 @@ namespace Cantera { m_redo_rates = true; } doublereal T = thermo(surfacePhaseIndex()).temperature(); + m_redo_rates = true; if (T != m_kdata->m_temp || m_redo_rates) { m_kdata->m_logtemp = log(T); m_rates.update(T, m_kdata->m_logtemp, DATA_PTR(m_kdata->m_rfn)); - if (m_has_electrochem_rxns) + if (m_has_exchange_current_density_formulation) { + applyExchangeCurrentDensityFormulation(DATA_PTR(m_kdata->m_rfn)); + } + if (m_has_electrochem_rxns) { applyButlerVolmerCorrection(DATA_PTR(m_kdata->m_rfn)); + } m_kdata->m_temp = T; updateKc(); m_kdata->m_ROP_ok = false; m_redo_rates = false; } } - - void InterfaceKinetics:: - _update_rates_phi() { + //==================================================================================================================== + void InterfaceKinetics::_update_rates_phi() { int np = nPhases(); for (int n = 0; n < np; n++) { if (thermo(n).electricPotential() != m_phi[n]) { @@ -189,6 +269,7 @@ namespace Cantera { } } } + //==================================================================================================================== /** @@ -198,8 +279,7 @@ namespace Cantera { * representing phases should overload to return the appropriate * quantities. */ - void InterfaceKinetics:: - _update_rates_C() { + void InterfaceKinetics::_update_rates_C() { int n; int np = nPhases(); @@ -263,6 +343,8 @@ namespace Cantera { } } } + //==================================================================================================================== + void InterfaceKinetics::checkPartialEquil() { @@ -338,6 +420,36 @@ namespace Cantera { } } + void InterfaceKinetics::getExchangeCurrentQuantities() { + /* + * First collect vectors of the standard Gibbs free energies of the + * species and the standard concentrations + * - m_mu0 + * - m_logStandardConc + */ + int ik = 0; + int np = nPhases(); + + for (int n = 0; n < np; n++) { + thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]); + int nsp = thermo(n).nSpecies(); + for (int k = 0; k < nsp; k++) { + m_StandardConc[ik] = thermo(n).standardConcentration(k); + ik++; + } + } + + m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), DATA_PTR(m_deltaG0)); + + + for (int i = 0; i < m_ii; i++) { + m_ProdStanConcReac[i] = 1.0; + } + + m_rxnstoich.multiplyReactants(DATA_PTR(m_StandardConc), DATA_PTR(m_ProdStanConcReac)); + + } + // Returns the Species creation rates [kmol/m^2/s]. /* * Return the species @@ -385,12 +497,17 @@ namespace Cantera { net); } - /** + //==================================================================================================================== + // Apply corrections for interfacial charge transfer reactions + /* * For reactions that transfer charge across a potential difference, * the activation energies are modified by the potential difference. * (see, for example, ...). This method applies this correction. + * + * @param kf Vector of forward reaction rate constants on which to have + * the correction applied */ - void InterfaceKinetics::applyButlerVolmerCorrection(doublereal* kf) { + void InterfaceKinetics::applyButlerVolmerCorrection(doublereal* const kf) { int i; int n, nsp, k, ik=0; @@ -410,8 +527,7 @@ namespace Cantera { // Compute the change in electrical potential energy for each // reaction. This will only be non-zero if a potential // difference is present. - m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_pot), - DATA_PTR(m_rwork)); + m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_pot), DATA_PTR(m_rwork)); // Modify the reaction rates. Only modify those with a // non-zero activation energy. Below we decrease the @@ -432,7 +548,8 @@ namespace Cantera { for (i = 0; i < nct; i++) { irxn = m_ctrxn[i]; eamod = m_beta[i]*m_rwork[irxn]; - if (eamod != 0.0 && m_E[irxn] != 0.0) { + // if (eamod != 0.0 && m_E[irxn] != 0.0) { + if (eamod != 0.0) { #ifdef DEBUG_KIN_MODE ea = GasConstant * m_E[irxn]; if (eamod + ea < 0.0) { @@ -450,18 +567,32 @@ namespace Cantera { } } } + //==================================================================================================================== + void InterfaceKinetics::applyExchangeCurrentDensityFormulation(doublereal* const kfwd) { + getExchangeCurrentQuantities(); + int nct = m_ctrxn.size(); + doublereal rt = GasConstant*thermo(0).temperature(); + doublereal rrt = 1.0/rt; + for (int i = 0; i < nct; i++) { + int irxn = m_ctrxn[i]; + int iECDFormulation = m_ctrxn_ecdf[i]; + if (iECDFormulation) { + double tmp = exp(- m_beta[i] * m_deltaG0[irxn] * rrt); + double tmp2 = m_ProdStanConcReac[irxn]; + tmp *= 1.0 / tmp2 / Faraday; + kfwd[irxn] *= tmp; + } + } - - + } + //==================================================================================================================== /** * Update the rates of progress of the reactions in the reaciton * mechanism. This routine operates on internal data. */ void InterfaceKinetics::getFwdRateConstants(doublereal* kfwd) { - // _update_rates_T(); - // _update_rates_C(); updateROP(); const vector_fp& rf = m_kdata->m_rfn; @@ -473,7 +604,7 @@ namespace Cantera { multiply_each(kfwd, kfwd + nReactions(), m_perturb.begin()); } - + //==================================================================================================================== /** * Update the rates of progress of the reactions in the reaciton @@ -493,12 +624,12 @@ namespace Cantera { multiply_each(krev, krev + nReactions(), rkc.begin()); } } - + //==================================================================================================================== void InterfaceKinetics::getActivationEnergies(doublereal *E) { copy(m_E.begin(), m_E.end(), E); } - + //==================================================================================================================== /** * Update the rates of progress of the reactions in the reaction * mechanism. This routine operates on internal data. @@ -545,7 +676,57 @@ namespace Cantera { for (int j = 0; j != m_ii; ++j) { ropnet[j] = ropf[j] - ropr[j]; - } + } + + + /* + * For reactions involving multiple phases, we must check that the phase + * being consumed actually exists. This is particularly important for + * phases that are stoichiometric phases containing one species with a unity activity + */ + if (m_phaseExistsCheck) { + for (int j = 0; j != m_ii; ++j) { + if ((ropr[j] > ropf[j]) && (ropr[j] > 0.0)) { + for (int p = 0; p < nPhases(); p++) { + if (m_rxnPhaseIsProduct[j][p]) { + if (! m_phaseExists[p]) { + ropnet[j] = 0.0; + ropr[j] = ropf[j]; + if (ropf[j] > 0.0) { + for (int rp = 0; rp < nPhases(); rp++) { + if (m_rxnPhaseIsReactant[j][rp]) { + if (! m_phaseExists[rp]) { + ropnet[j] = 0.0; + ropr[j] = ropf[j] = 0.0;; + } + } + } + } + } + } + } + } else if ((ropf[j] > ropr[j]) && (ropf[j] > 0.0)) { + for (int p = 0; p < nPhases(); p++) { + if (m_rxnPhaseIsReactant[j][p]) { + if (! m_phaseExists[p]) { + ropnet[j] = 0.0; + ropf[j] = ropr[j]; + if (ropf[j] > 0.0) { + for (int rp = 0; rp < nPhases(); rp++) { + if (m_rxnPhaseIsProduct[j][rp]) { + if (! m_phaseExists[rp]) { + ropnet[j] = 0.0; + ropf[j] = ropr[j] = 0.0; + } + } + } + } + } + } + } + } + } + } m_kdata->m_ROP_ok = true; } @@ -738,7 +919,7 @@ namespace Cantera { m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaS); } - + //==================================================================================================================== /** * Add a single reaction to the mechanism. This routine * must be called after init() and before finalize(). @@ -752,21 +933,44 @@ namespace Cantera { * There is no difference between elementary and surface * reactions. */ - void InterfaceKinetics:: - addReaction(const ReactionData& r) { + void InterfaceKinetics::addReaction(const ReactionData& r) { addElementaryReaction(r); // operations common to all reaction types - installReagents( r ); + installReagents(r); //installGroups(reactionNumber(), r.rgroups, r.pgroups); incrementRxnCount(); m_rxneqn.push_back(r.equation); - } + m_rxnPhaseIsReactant.resize(m_ii, 0); + m_rxnPhaseIsProduct.resize(m_ii, 0); + + int np = nPhases(); + int i = m_ii -1; + m_rxnPhaseIsReactant[i] = new bool[np]; + m_rxnPhaseIsProduct[i] = new bool[np]; - void InterfaceKinetics:: - addElementaryReaction(const ReactionData& r) { + for (int p = 0; p < np; p++) { + m_rxnPhaseIsReactant[i][p] = false; + m_rxnPhaseIsProduct[i][p] = false; + } + + const vector_int& vr = reactants(i); + for (int ik = 0; ik < (int) vr.size(); ik++) { + int k = vr[ik]; + int p = speciesPhaseIndex(k); + m_rxnPhaseIsReactant[i][p] = true; + } + const vector_int& vp = products(i); + for (int ik = 0; ik < (int) vp.size(); ik++) { + int k = vp[ik]; + int p = speciesPhaseIndex(k); + m_rxnPhaseIsProduct[i][p] = true; + } + } + //==================================================================================================================== + void InterfaceKinetics::addElementaryReaction(const ReactionData& r) { int iloc; // install rate coeff calculator @@ -776,12 +980,11 @@ namespace Cantera { if (ncov > 3) { m_has_coverage_dependence = true; } - for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]); - - iloc = m_rates.install( reactionNumber(), - r.rateCoeffType, rp.size(), - DATA_PTR(rp) ); - + for (int m = 0; m < ncov; m++) { + rp.push_back(r.cov[m]); + } + // iloc = m_rates.install(reactionNumber(), r.rateCoeffType, rp.size(), DATA_PTR(rp)); + iloc = m_rates.install(reactionNumber(), ARRHENIUS_REACTION_RATECOEFF_TYPE, rp.size(), DATA_PTR(rp)); // store activation energy m_E.push_back(r.rateCoeffParameters[2]); @@ -789,13 +992,19 @@ namespace Cantera { m_has_electrochem_rxns = true; m_beta.push_back(r.beta); m_ctrxn.push_back(reactionNumber()); + if (r.rateCoeffType == EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE) { + m_has_exchange_current_density_formulation = true; + m_ctrxn_ecdf.push_back(1); + } else { + m_ctrxn_ecdf.push_back(0); + } } // add constant term to rate coeff value vector m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]); registerReaction( reactionNumber(), ELEMENTARY_RXN, iloc); } - + //==================================================================================================================== void InterfaceKinetics::setIOFlag(int ioFlag) { m_ioFlag = ioFlag; @@ -912,7 +1121,7 @@ namespace Cantera { * calculates rates of species production from reaction rates of * progress. */ - m_rxnstoich.add( reactionNumber(), r); + m_rxnstoich.add(reactionNumber(), r); /* * register reaction in lists of reversible and irreversible rxns. */ @@ -920,11 +1129,16 @@ namespace Cantera { m_revindex.push_back(reactionNumber()); m_nrev++; } else { - m_irrev.push_back( reactionNumber() ); + m_irrev.push_back(reactionNumber()); m_nirrev++; } } - + //=============================================================================================== + void InterfaceKinetics::addPhase(thermo_t &thermo) { + Kinetics::addPhase(thermo); + m_phaseExists.push_back(true); + } + //================================================================================================ /** * Prepare the class for the addition of reactions. This function * must be called after instantiation of the class, but before @@ -948,7 +1162,7 @@ namespace Cantera { m_pot.resize(m_kk, 0.0); m_phi.resize(np, 0.0); } - + //================================================================================================ /** * Finish adding reactions and prepare for use. This function * must be called after all reactions are entered into the mechanism @@ -958,6 +1172,7 @@ namespace Cantera { * since we don't know this number up to now. */ void InterfaceKinetics::finalize() { + Kinetics::finalize(); m_rwork.resize(nReactions()); int ks = reactionPhaseIndex(); if (ks < 0) throw CanteraError("InterfaceKinetics::finalize", @@ -967,10 +1182,20 @@ namespace Cantera { throw CanteraError("InterfaceKinetics::finalize", "expected interface dimension = 2, but got dimension = " +int2str(m_surf->nDim())); + + + + m_StandardConc.resize(m_nTotalSpecies, 0.0); + m_deltaG0.resize(m_ii, 0.0); + m_ProdStanConcReac.resize(m_ii, 0.0); + + if (m_thermo.size() != m_phaseExists.size()) { + throw CanteraError("InterfaceKinetics::finalize", "internal error"); + } + m_finalized = true; } - doublereal InterfaceKinetics::electrochem_beta(int irxn) const{ int n = m_ctrxn.size(); for (int i = 0; i < n; i++) { @@ -985,7 +1210,7 @@ namespace Cantera { bool InterfaceKinetics::ready() const { return (m_finalized); } - + //================================================================================================ // Advance the surface coverages in time /* * @param tstep Time value to advance the surface coverages @@ -1002,7 +1227,7 @@ namespace Cantera { delete m_integrator; m_integrator = 0; } - + //================================================================================================ // Solve for the pseudo steady-state of the surface problem /* * Solve for the steady state of the surface problem. @@ -1031,7 +1256,25 @@ namespace Cantera { */ m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride); } + //================================================================================================ + void InterfaceKinetics::setPhaseExistence(const int iphase, const bool exists) { + if (iphase < 0 || iphase >= (int) m_thermo.size()) { + throw CanteraError("InterfaceKinetics:setPhaseExistence", "out of bounds"); + } + if (exists) { + if (!m_phaseExists[iphase]) { + m_phaseExistsCheck--; + m_phaseExists[iphase] = true; + } + } else { + if (m_phaseExists[iphase]) { + m_phaseExistsCheck++; + m_phaseExists[iphase] = false; + } + } + } + //================================================================================================ void EdgeKinetics::finalize() { m_rwork.resize(nReactions()); int ks = reactionPhaseIndex(); @@ -1043,7 +1286,8 @@ namespace Cantera { "expected interface dimension = 1, but got dimension = " +int2str(m_surf->nDim())); m_finalized = true; - } + } + //================================================================================================ } diff --git a/Cantera/src/kinetics/InterfaceKinetics.h b/Cantera/src/kinetics/InterfaceKinetics.h index c1b0343a8..50ca5dc59 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.h +++ b/Cantera/src/kinetics/InterfaceKinetics.h @@ -45,11 +45,16 @@ namespace Cantera { class InterfaceKineticsData { public: InterfaceKineticsData() : + m_logp0(0.0), + m_logc0(0.0), m_ROP_ok(false), - m_temp(0.0), m_logtemp(0.0) + m_temp(0.0), + m_logtemp(0.0) {} + //! Virtual destructor - virtual ~InterfaceKineticsData(){} + virtual ~InterfaceKineticsData() { + } doublereal m_logp0; doublereal m_logc0; @@ -59,7 +64,9 @@ namespace Cantera { bool m_ROP_ok; + //! Current temperature of the data doublereal m_temp; + //! Current log of the temperature doublereal m_logtemp; vector_fp m_rfn; vector_fp m_rkcn; @@ -76,9 +83,9 @@ namespace Cantera { public: - /** - * Constructor - * + + //! Constructor + /*! * @param thermo The optional parameter may be used to initialize * the object with one ThermoPhase object. * HKM Note -> Since the interface kinetics @@ -121,20 +128,18 @@ namespace Cantera { */ virtual Kinetics *duplMyselfAsKinetics() const; + //! Return the ID of the kinetics object + virtual int ID() const; - virtual int ID() const { return cInterfaceKinetics; } - virtual int type() const { return cInterfaceKinetics; } + //! Retunr the type of the kinetics object + virtual int type() const; - /** - * Set the electric potential in the nth phase - * + //! Set the electric potential in the nth phase + /*! * @param n phase Index in this kinetics object. * @param V Electric potential (volts) */ - void setElectricPotential(int n, doublereal V) { - thermo(n).setElectricPotential(V); - m_redo_rates = true; - } + void setElectricPotential(int n, doublereal V); /// @@ -172,8 +177,17 @@ namespace Cantera { std::copy(m_kdata->m_ropnet.begin(), m_kdata->m_ropnet.end(), netROP); } + + //! Get the equilibrium constants of all reactions, whether + //! the reaction is reversible or not. + /*! + * @param kc Returns the concentration equation constant for the reaction. + * Length is the number of reactions + */ virtual void getEquilibriumConstants(doublereal* kc); + void getExchangeCurrentQuantities(); + virtual void getDeltaGibbs( doublereal* deltaG); @@ -361,6 +375,19 @@ namespace Cantera { */ //@{ + //! Add a phase to the kinetics manager object. + /*! + * This must be done before the function init() is called or + * before any reactions are input. + * + * This function calls the Kinetics operator addPhase. + * It also sets the following functions + * + * m_phaseExists[] + * + * @param thermo Reference to the ThermoPhase to be added. + */ + virtual void addPhase(thermo_t& thermo); //! Prepare the class for the addition of reactions. /*! @@ -390,12 +417,26 @@ namespace Cantera { virtual bool ready() const; - + //! Internal routine that updates the Rates of Progress of the reactions + /*! + * This is actually the guts of the functionality of the object + */ void updateROP(); - + + //! Update properties that depend on temperature + /*! + * This is called to update all of the properties that depend on temperature + * + * Current objects that this function updates + * m_kdata->m_logtemp + * m_kdata->m_rfn + * m_rates. + * updateKc(); + */ void _update_rates_T(); + void _update_rates_phi(); void _update_rates_C(); @@ -445,15 +486,63 @@ namespace Cantera { void setIOFlag(int ioFlag); void checkPartialEquil(); - + + + int reactionNumber() const { return m_ii;} + + void addElementaryReaction(const ReactionData& r); + void addGlobalReaction(const ReactionData& r); + void installReagents(const ReactionData& r); + + void updateKc(); + + //! Write values into m_index + /*! + * @param rxnNumber reaction number + * @param type reaction type + * @param loc location ?? + */ + void registerReaction(int rxnNumber, int type, int loc) { + m_index[rxnNumber] = std::pair(type, loc); + } + + //! Apply corrections for interfacial charge transfer reactions + /*! + * For reactions that transfer charge across a potential difference, + * the activation energies are modified by the potential difference. + * (see, for example, ...). This method applies this correction. + * + * @param kf Vector of forward reaction rate constants on which to have + * the correction applied + */ + void applyButlerVolmerCorrection(doublereal* const kf); + + //! When an electrode reaction rate is optionally specified in terms of its + //! exchange current density, extra vectors need to be precalculated + /*! + * + */ + void applyExchangeCurrentDensityFormulation(doublereal* const kfwd); + + //! Set the existence of a phase in the reaction object + /*! + * Tell the kinetics object whether a phase in the object exists. + * This is actually an extrinsic specification that must be carried out on top of the + * intrinsic calculation of the reaction rate + * + * @param iphase Index of the phase. This is the order within the internal thermo vector object + * @param exists Boolean indicating whether the phase exists or not + */ + void setPhaseExistence(const int iphase, const bool exists); + + protected: + //! Temporary work vector of length m_kk vector_fp m_grt; - protected: - //! m_kk is the number of species in all of the phases //! that participate in this kinetics mechanism. - int m_kk; + int m_kk; //! List of reactions numbers which are reversible reactions /*! @@ -489,7 +578,7 @@ namespace Cantera { * production rates and also handles turning thermo * properties into reaction thermo properties. */ - ReactionStoichMgr m_rxnstoich; + ReactionStoichMgr m_rxnstoich; //! Number of irreversible reactions in the mechanism int m_nirrev; @@ -564,7 +653,7 @@ namespace Cantera { */ vector_fp m_mu0; - //! Vector of phase potentials + //! Vector of phase electric potentials /*! * Temporary vector containing the potential of each phase * in the kinetics object @@ -607,30 +696,36 @@ namespace Cantera { ImplicitSurfChem* m_integrator; vector_fp m_beta; + + //! Vector of reaction indexes specifying the id of the current transfer reactions + //! in the mechanism + /*! + * Vector of reaction indecices which involve current transfers. This provides + * an index into the m_beta array. + * + * irxn = m_ctrxn[i] + */ vector_int m_ctrxn; - int reactionNumber(){ return m_ii;} + //! Vector of booleans indicating whether the charge transfer reaction may be + //! described by an exchange current density expression + vector_int m_ctrxn_ecdf; - void addElementaryReaction(const ReactionData& r); - void addGlobalReaction(const ReactionData& r); - void installReagents(const ReactionData& r); + vector_fp m_StandardConc; + vector_fp m_deltaG0; + vector_fp m_ProdStanConcReac; - void updateKc(); - //! Write values into m_index - /*! - * @param rxnNumber reaction number - * @param type reaction type - * @param loc location ?? - */ - void registerReaction(int rxnNumber, int type, int loc) { - m_index[rxnNumber] = std::pair(type, loc); - } - - void applyButlerVolmerCorrection(doublereal* kf); //! boolean indicating whether mechanism has been finalized bool m_finalized; + + //! Boolean flag indicating whether any reaction in the mechanism + //! has a coverage dependent forward reaction rate + /*! + * If this is true, then the coverage dependence is multiplied into + * the forward reaction rates constant + */ bool m_has_coverage_dependence; //! Boolean flag indicating whether any reaction in the mechanism @@ -643,6 +738,37 @@ namespace Cantera { */ bool m_has_electrochem_rxns; + //! Boolean flag indicating whether any reaction in the mechanism + //! is described by an exchange current density expression + /*! + * If this is true, the standard state gibbs free energy of the reaction and + * the product of the reactant standard concentrations must be precalculated + * in order to calculate the rate constant. + */ + bool m_has_exchange_current_density_formulation; + + //! Int flag to indicate that some phases in the kinetics mechanism are + //! non-existent. + /*! + * We change the ROP vectors to make sure that non-existent phases are treated + * correctly in the kinetics operator. The value of this is equal to the number + * of phases which don't exist. + */ + int m_phaseExistsCheck; + + //! Vector of booleans indicating whether phases exist or not + /*! + * Vector of booleans indicating whether a phase exists or not. + * We use this to set the ROP's so that unphysical things don't happen + * + * length = number of phases in the object + * By default all phases exist. + */ + std::vector m_phaseExists; + + std::vector m_rxnPhaseIsReactant; + std::vector m_rxnPhaseIsProduct; + int m_ioFlag; private: diff --git a/Cantera/src/kinetics/Kinetics.cpp b/Cantera/src/kinetics/Kinetics.cpp index 7ab58f4f8..c2338e1e1 100644 --- a/Cantera/src/kinetics/Kinetics.cpp +++ b/Cantera/src/kinetics/Kinetics.cpp @@ -29,9 +29,16 @@ using namespace std; namespace Cantera { - Kinetics::Kinetics() : m_ii(0), m_thermo(0), - m_index(-1), m_surfphase(-1), m_rxnphase(-1), - m_mindim(4) {} + Kinetics::Kinetics() : + m_ii(0), + m_nTotalSpecies(0), + m_thermo(0), + m_index(-1), + m_surfphase(-1), + m_rxnphase(-1), + m_mindim(4) + { + } Kinetics::~Kinetics(){} @@ -42,7 +49,8 @@ namespace Cantera { * throw an exception. */ Kinetics::Kinetics(const Kinetics &right) : - m_ii(0), + m_ii(0), + m_nTotalSpecies(0), m_thermo(0), m_index(-1), m_surfphase(-1), @@ -70,6 +78,7 @@ namespace Cantera { if (this == &right) return *this; m_ii = right.m_ii; + m_nTotalSpecies = right.m_nTotalSpecies; m_perturb = right.m_perturb; m_reactants = right.m_reactants; m_products = right.m_products; @@ -216,7 +225,8 @@ namespace Cantera { return thermo(0); } - /** + //============================================================================================== + /* * This function takes as an argument the kineticsSpecies index * (i.e., the list index in the list of species in the kinetics * manager) and returns the index of the phase owning the @@ -286,6 +296,15 @@ namespace Cantera { m_phaseindex[m_thermo.back()->id()] = nPhases(); } + void Kinetics::finalize() { + m_nTotalSpecies = 0; + int np = nPhases(); + for (int n = 0; n < np; n++) { + int nsp = m_thermo[n]->nSpecies(); + m_nTotalSpecies += nsp; + } + } + //! Private function of the class Kinetics, indicating that a function //! inherited from the base class hasn't had a definition assigned to it diff --git a/Cantera/src/kinetics/Kinetics.h b/Cantera/src/kinetics/Kinetics.h index 25ca5305f..33b63263b 100644 --- a/Cantera/src/kinetics/Kinetics.h +++ b/Cantera/src/kinetics/Kinetics.h @@ -811,7 +811,7 @@ namespace Cantera { * any initialization (allocating arrays, etc.) that must be * done after the reactions are entered. */ - virtual void finalize() {} + virtual void finalize(); /** * Add a single reaction to the mechanism. This routine @@ -904,6 +904,9 @@ namespace Cantera { //! Number of reactions in the mechanism int m_ii; + + //! Number of species in the species vector for this kinetics operator + int m_nTotalSpecies; /// Vector of perturbation factors for each reaction's rate of /// progress vector. It is initialized to one. diff --git a/Cantera/src/kinetics/ReactionData.h b/Cantera/src/kinetics/ReactionData.h index a02d25f8e..5c4f8c466 100644 --- a/Cantera/src/kinetics/ReactionData.h +++ b/Cantera/src/kinetics/ReactionData.h @@ -25,7 +25,7 @@ namespace Cantera { number = 0; rxn_number = 0; reversible = true; - rateCoeffType = ARRHENIUS; + rateCoeffType = ARRHENIUS_REACTION_RATECOEFF_TYPE; falloffType = NONE; error = 0; equation = ""; @@ -36,7 +36,12 @@ namespace Cantera { } virtual ~ReactionData(){} + //! type of the reaction + /*! + * The valid types are listed in the file, reaction_defs.h. + */ int reactionType; + int number; int rxn_number; vector_int reactants; @@ -51,7 +56,14 @@ namespace Cantera { //! True if the current reaction is reversible. False otherwise bool reversible; + + //! type of the rate coefficient for the forward rate constant + /*! + * The valid types are listed in the file, reaction_defs.h and they + * all end in RATECOEFF_TYPE + */ int rateCoeffType; + vector_fp rateCoeffParameters; vector_fp auxRateCoeffParameters; int falloffType; diff --git a/Cantera/src/kinetics/ReactionStoichMgr.h b/Cantera/src/kinetics/ReactionStoichMgr.h index 203bfc5e4..eba12e5d5 100644 --- a/Cantera/src/kinetics/ReactionStoichMgr.h +++ b/Cantera/src/kinetics/ReactionStoichMgr.h @@ -214,6 +214,7 @@ namespace Cantera { * \f[ * R_i = R_i * \prod_k C_k^{o_{k,i}} * \f] + * * Here \f$ o_{k,i} \f$ is the reaction order of species k in reaction i. */ virtual void multiplyReactants(const doublereal* C, doublereal* R); diff --git a/Cantera/src/kinetics/RxnRates.h b/Cantera/src/kinetics/RxnRates.h index 4e3d495ad..e8b0b5284 100644 --- a/Cantera/src/kinetics/RxnRates.h +++ b/Cantera/src/kinetics/RxnRates.h @@ -18,94 +18,99 @@ #include "ctexceptions.h" namespace Cantera { - - /** - * A rate coefficient of the form - * \f[ - * A T^b \exp (-E/RT) - * \f] + + //! Arrhenius reaction rate type depends only on temperature + /** + * A reaction rate coefficient of the following form. + * + * \f[ + * k_f = A T^b \exp (-E/RT) + * \f] + * + */ + class Arrhenius { + + public: + + /// return the rate coefficient type. + static int type(){ return ARRHENIUS_REACTION_RATECOEFF_TYPE ; } + + /// Default constructor. + Arrhenius() : + m_logA(-1.0E300), + m_b (0.0), + m_E (0.0), + m_A(0.0) {} + + /// Constructor with Arrhenius parameters specified with an array. + Arrhenius(int csize, const doublereal* c) : + m_b (c[1]), + m_E (c[2]), + m_A (c[0]) + { + if (m_A <= 0.0) { + m_logA = -1.0E300; + } else { + m_logA = log(m_A); + } + } + + /// Constructor. + /// @param A pre-exponential. The unit system is + /// (kmol, m, s). The actual units depend on the reaction + /// order and the dimensionality (surface or bulk). + /// @param b Temperature exponent. Non-dimensional. + /// @param E Activation energy in temperature units. Kelvin. + Arrhenius(doublereal A, doublereal b, doublereal E) : + m_b (b), + m_E (E), + m_A (A) + { + if (m_A <= 0.0) { + m_logA = -1.0E300; + } else { + m_logA = log(m_A); + } + } + + //! Update concentration-dependent parts of the rate coefficient. + /*! + * For this class, there are no + * concentration-dependent parts, so this method does nothing. */ - class Arrhenius { - - public: + void update_C(const doublereal* c) { + } - /// return the rate coefficient type. - static int type(){ return ARRHENIUS; } - - /// Default constructor. - Arrhenius() : - m_logA(-1.0E300), - m_b (0.0), - m_E (0.0), - m_A(0.0) {} - - /// Constructor with Arrhenius parameters specified with an array. - Arrhenius(int csize, const doublereal* c) : - m_b (c[1]), - m_E (c[2]), - m_A (c[0]) - { - if (m_A <= 0.0) { - m_logA = -1.0E300; - } else { - m_logA = log(m_A); - } - } - - /// Constructor. - /// @param A pre-exponential. The unit system is - /// (kmol, m, s). The actual units depend on the reaction - /// order and the dimensionality (surface or bulk). - /// @param b Temperature exponent. Non-dimensional. - /// @param E Activation energy in temperature units. Kelvin. - Arrhenius(doublereal A, doublereal b, doublereal E) : - m_b (b), - m_E (E), - m_A (A) - { - if (m_A <= 0.0) { - m_logA = -1.0E300; - } else { - m_logA = log(m_A); - } - } - - /// Update concentration-dependent parts of the rate - /// coefficient. For this class, there are no - /// concentration-dependent parts, so this method does - /// nothing. - void update_C(const doublereal* c) {} + /** + * Update the value of the logarithm of the rate constant. + * + * Note, this function should never be called for negative A values. + * If it does then it will produce a negative overflow result, and + * a zero net forwards reaction rate, instead of a negative reaction + * rate constant that is the expected result. + */ + doublereal update(doublereal logT, doublereal recipT) const { + return m_logA + m_b*logT - m_E*recipT; + } - /** - * Update the value of the logarithm of the rate constant. - * - * Note, this function should never be called for negative A values. - * If it does then it will produce a negative overflow result, and - * a zero net forwards reaction rate, instead of a negative reaction - * rate constant that is the expected result. - */ - doublereal update(doublereal logT, doublereal recipT) const { - return m_logA + m_b*logT - m_E*recipT; - } - - /** - * Update the value the rate constant. - * - * This function returns the actual value of the rate constant. - * It can be safely called for negative values of the pre-exponential - * factor. - */ - doublereal updateRC(doublereal logT, doublereal recipT) const { - return m_A * exp(m_b*logT - m_E*recipT); - } + /** + * Update the value the rate constant. + * + * This function returns the actual value of the rate constant. + * It can be safely called for negative values of the pre-exponential + * factor. + */ + doublereal updateRC(doublereal logT, doublereal recipT) const { + return m_A * exp(m_b*logT - m_E*recipT); + } - void writeUpdateRHS(std::ostream& s) const { - s << " exp(" << m_logA; - if (m_b != 0.0) s << " + " << m_b << " * tlog"; - if (m_E != 0.0) s << " - " << m_E << " * rt"; - s << ");" << std::endl; - } + void writeUpdateRHS(std::ostream& s) const { + s << " exp(" << m_logA; + if (m_b != 0.0) s << " + " << m_b << " * tlog"; + if (m_E != 0.0) s << " - " << m_E << " * rt"; + s << ");" << std::endl; + } doublereal activationEnergy_R() const { return m_E; @@ -121,19 +126,19 @@ namespace Cantera { class ArrheniusSum { public: - static int type(){ return ARRHENIUS_SUM; } + static int type(){ return ARRHENIUS_SUM_REACTION_RATECOEFF_TYPE; } ArrheniusSum() : m_nterms(0) {} void addArrheniusTerm(doublereal A, doublereal b, doublereal E) { - if (A > 0.0) { - m_terms.push_back(Arrhenius(A, b, E)); - m_sign.push_back(1); - } - else if (A < 0.0) { - m_terms.push_back(Arrhenius(-A, b, E)); - m_sign.push_back(-1); - } - m_nterms++; + if (A > 0.0) { + m_terms.push_back(Arrhenius(A, b, E)); + m_sign.push_back(1); + } + else if (A < 0.0) { + m_terms.push_back(Arrhenius(-A, b, E)); + m_sign.push_back(-1); + } + m_nterms++; } void update_C(const doublereal* c) {} @@ -169,38 +174,38 @@ namespace Cantera { return fsum; } - void writeUpdateRHS(std::ostream& s) const { + void writeUpdateRHS(std::ostream& s) const { ; } static bool alwaysComputeRate() { return false;} protected: - std::vector m_terms; - vector_int m_sign; - int m_nterms; + std::vector m_terms; + vector_int m_sign; + int m_nterms; }; - /** - * An Arrhenius rate with coverage-dependent terms. - */ - class SurfaceArrhenius { + /** + * An Arrhenius rate with coverage-dependent terms. + */ + class SurfaceArrhenius { - public: - static int type(){ return ARRHENIUS; } - SurfaceArrhenius() : - m_logA(-1.0E300), - m_b (0.0), - m_E (0.0), - m_A(0.0), - m_acov(0.0), - m_ecov(0.0), - m_mcov(0.0), - m_ncov(0), - m_nmcov(0) - { - } + public: + static int type(){ return ARRHENIUS_REACTION_RATECOEFF_TYPE ; } + SurfaceArrhenius() : + m_logA(-1.0E300), + m_b (0.0), + m_E (0.0), + m_A(0.0), + m_acov(0.0), + m_ecov(0.0), + m_mcov(0.0), + m_ncov(0), + m_nmcov(0) + { + } SurfaceArrhenius( int csize, const doublereal* c ) : m_b (c[1]), @@ -325,7 +330,7 @@ namespace Cantera { return exp(lres); } - void writeUpdateRHS(std::ostream& s) const {} + void writeUpdateRHS(std::ostream& s) const {} protected: doublereal delta_s0, delta_e0; @@ -335,37 +340,143 @@ namespace Cantera { #endif -} -// class LandauTeller { -// public: -// static int type(){ return LANDAUTELLER; } -// LandauTeller(){} -// LandauTeller( const vector_fp& c ) : m_c(c) { m_c[0] = log(c[0]); } + //! Arrhenius reaction rate type depends only on temperature + /** + * A reaction rate coefficient of the following form. + * + * \f[ + * k_f = A T^b \exp (-E/RT) + * \f] + * + */ + class ExchangeCurrent { + + public: -// doublereal update(doublereal logT, doublereal recipT) const { -// return m_c[0] + m_c[1]*tt[1] - m_c[2]*tt[2] -// + m_c[3]*tt[3] + m_c[4]*tt[4]; -// } - -// //void writeUpdateRHS(ostream& s) const { -// // s << exp(m_logA); -// // s << " * exp("; -// // if (m_b != 0.0) s << m_b << " * tlog"; -// // if (m_E != 0.0) s << " - " << m_E << " * rt"; -// // if (m_E != 0.0) s << " - " << m_E << " * rt"; -// // s << ");" << endl; -// // } -// //} + //! return the rate coefficient type. + static int type() { + return EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE ; + } -// protected: -// doublereal m_logA, m_b, m_E; -// }; + //! Default constructor. + ExchangeCurrent() : + m_logA(-1.0E300), + m_b (0.0), + m_E (0.0), + m_A(0.0) {} + + //! Constructor with Arrhenius parameters specified with an array. + ExchangeCurrent(int csize, const doublereal* c) : + m_b (c[1]), + m_E (c[2]), + m_A (c[0]) + { + if (m_A <= 0.0) { + m_logA = -1.0E300; + } else { + m_logA = log(m_A); + } + } + + /// Constructor. + /// @param A pre-exponential. The unit system is + /// (kmol, m, s). The actual units depend on the reaction + /// order and the dimensionality (surface or bulk). + /// @param b Temperature exponent. Non-dimensional. + /// @param E Activation energy in temperature units. Kelvin. + ExchangeCurrent(doublereal A, doublereal b, doublereal E) : + m_b (b), + m_E (E), + m_A (A) + { + if (m_A <= 0.0) { + m_logA = -1.0E300; + } else { + m_logA = log(m_A); + } + } + + //! Update concentration-dependent parts of the rate coefficient. + /*! + * For this class, there are no + * concentration-dependent parts, so this method does nothing. + */ + void update_C(const doublereal* c) { + } + + /** + * Update the value of the logarithm of the rate constant. + * + * Note, this function should never be called for negative A values. + * If it does then it will produce a negative overflow result, and + * a zero net forwards reaction rate, instead of a negative reaction + * rate constant that is the expected result. + */ + doublereal update(doublereal logT, doublereal recipT) const { + return m_logA + m_b*logT - m_E*recipT; + } + + /** + * Update the value the rate constant. + * + * This function returns the actual value of the rate constant. + * It can be safely called for negative values of the pre-exponential + * factor. + */ + doublereal updateRC(doublereal logT, doublereal recipT) const { + return m_A * exp(m_b*logT - m_E*recipT); + } + + + void writeUpdateRHS(std::ostream& s) const { + s << " exp(" << m_logA; + if (m_b != 0.0) s << " + " << m_b << " * tlog"; + if (m_E != 0.0) s << " - " << m_E << " * rt"; + s << ");" << std::endl; + } + + doublereal activationEnergy_R() const { + return m_E; + } + + static bool alwaysComputeRate() { return false;} + + protected: + doublereal m_logA, m_b, m_E, m_A; + }; + + + // class LandauTeller { + + // public: + // static int type(){ return LANDAUTELLER; } + // LandauTeller(){} + // LandauTeller( const vector_fp& c ) : m_c(c) { m_c[0] = log(c[0]); } + + // doublereal update(doublereal logT, doublereal recipT) const { + // return m_c[0] + m_c[1]*tt[1] - m_c[2]*tt[2] + // + m_c[3]*tt[3] + m_c[4]*tt[4]; + // } + + // //void writeUpdateRHS(ostream& s) const { + // // s << exp(m_logA); + // // s << " * exp("; + // // if (m_b != 0.0) s << m_b << " * tlog"; + // // if (m_E != 0.0) s << " - " << m_E << " * rt"; + // // if (m_E != 0.0) s << " - " << m_E << " * rt"; + // // s << ");" << endl; + // // } + // //} + + // protected: + // doublereal m_logA, m_b, m_E; + // }; -//} - + //} +} #endif diff --git a/Cantera/src/kinetics/importKinetics.cpp b/Cantera/src/kinetics/importKinetics.cpp index e40a1df84..267c1b76c 100644 --- a/Cantera/src/kinetics/importKinetics.cpp +++ b/Cantera/src/kinetics/importKinetics.cpp @@ -521,9 +521,24 @@ namespace Cantera { * kf should point to a XML element named "rateCoeff". * rdata is the partially filled ReactionData object for the reaction. * This function will fill in more fields in the ReactionData object. + * + * @param kf Reference to the XML Node named rateCoeff */ void getRateCoefficient(const node_t& kf, kinetics_t& kin, ReactionData& rdata, int negA) { + string type = kf.attrib("type"); + if (type == "") { + type = "Arrhenius"; + rdata.rateCoeffType = ARRHENIUS_REACTION_RATECOEFF_TYPE; + } + if (type == "ExchangeCurrentDensity") { + rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE; + } else if (type == "Arrhenius") { + + } else { + throw CanteraError("getRateCoefficient", + "Unknown type: " + type); + } int nc = kf.nChildren(); nodeset_t& kf_children = kf.children(); @@ -557,6 +572,12 @@ namespace Cantera { "negative or zero A coefficient for reaction "+int2str(rdata.number)); } } + else if (nm == "Arrhenius_ExchangeCurrentDensity") { + vector_fp coeff(3); + getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]); + chigh = coeff; + rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE; + } else if (nm == "falloff") { getFalloff(c, rdata); } @@ -580,6 +601,9 @@ namespace Cantera { rdata.auxRateCoeffParameters = clow; else if (rdata.reactionType == CHEMACT_RXN) rdata.auxRateCoeffParameters = chigh; + + + } diff --git a/Cantera/src/kinetics/reaction_defs.h b/Cantera/src/kinetics/reaction_defs.h index eb36dc544..f709d36a7 100644 --- a/Cantera/src/kinetics/reaction_defs.h +++ b/Cantera/src/kinetics/reaction_defs.h @@ -81,11 +81,12 @@ namespace Cantera { */ //@{ - const int ARRHENIUS = 1; - const int LANDAUTELLER = 2; - const int TSTRATE = 3; - const int SURF_ARRHENIUS = 4; - const int ARRHENIUS_SUM = 5; + const int ARRHENIUS_REACTION_RATECOEFF_TYPE = 1; + const int LANDAUTELLER_REACTION_RATECOEFF_TYPE = 2; + const int TSTRATE_REACTION_RATECOEFF_TYPE = 3; + const int SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE = 4; + const int ARRHENIUS_SUM_REACTION_RATECOEFF_TYPE = 5; + const int EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE = 6; //@}