diff --git a/include/cantera/kinetics/ElectrodeKinetics.h b/include/cantera/kinetics/ElectrodeKinetics.h index 8052a036b..d9159bf6f 100644 --- a/include/cantera/kinetics/ElectrodeKinetics.h +++ b/include/cantera/kinetics/ElectrodeKinetics.h @@ -68,21 +68,24 @@ public: */ void identifyMetalPhase(); - - + //! Internal routine that updates the Rates of Progress of the reactions + /*! + * This is actually the guts of the functionality of the object + */ + virtual void updateROP(); //void addGlobalReaction(ReactionData& r); - + double calcForwardROP_BV(size_t irxn, size_t iBeta); protected: - - //! index of the metal phase in the list of phases for this surface + //! Index of the metal phase in the list of phases for this kinetics object size_t metalPhaseRS_; + //! Index of the electron phase in the list of phases for this kinetics object size_t electronPhaseRS_; //! Index of the solution phase in the list of phases for this surface @@ -91,9 +94,6 @@ protected: //! Index of the electrons species in the list of species for this surface kinetics, if none set it to -1 size_t kElectronRS_; - - - }; } diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index 2165612f3..0c14bc8bc 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -24,6 +24,32 @@ class SurfPhase; class ImplicitSurfChem; class RxnMolChange; +//! forward orders +class RxnOrders { + + public: + //! constructors + RxnOrders(); + + RxnOrders(const RxnOrders &right); + + ~RxnOrders(); + + RxnOrders& operator=(const RxnOrders &right); + + //! Fill in the structure with the array. + /*! + * @param[in] Size of length kinetic species. The entries the values of the orders + */ + int fill(const std::vector& fullForwardOrders); + + //! ID's of the kinetic species + std::vector kinSpeciesIDs_; + + //! Orders of the kinetic species + std::vector kinSpeciesOrders_; +}; + //! A kinetics manager for heterogeneous reaction mechanisms. The //! reactions are assumed to occur at a 2D interface between two 3D phases. /*! @@ -139,28 +165,35 @@ public: //! @name Reaction Mechanism Informational Query Routines //! @{ - virtual doublereal reactantStoichCoeff(size_t k, size_t i) const { - return m_rrxn[k][i]; - } + //! Provide a reactant stoichiometric coefficient + /*! + * @param[in] kSpecKin Species index within the kinetics object + * @param[in] irxn Reaction index + * + * @return Returns the reactant stoichiometic coefficient within the reaction + */ + virtual doublereal reactantStoichCoeff(size_t kSpecKin, size_t irxn) const; - virtual doublereal productStoichCoeff(size_t k, size_t i) const { - return m_prxn[k][i]; - } + //! Provide a product stoichiometric coefficient + /*! + * @param[in] kSpecKin Species index within the kinetics object + * @param[in] irxn Reaction index + * + * @return Returns the product stoichiometic coefficient within the reaction + */ + virtual doublereal productStoichCoeff(size_t kSpecKin, size_t irxn) const; - //! return the reaction type of the reaction i + //! return the reaction type of the reaction irxn /*! * @param[in] Reaction index * * @return Returns the reaction type of the reaction. */ - virtual int reactionType(size_t i) const; - //virtual int reactionType(size_t i) const { - // return m_index[i].first; - //} - + virtual int reactionType(size_t irxn) const; virtual void getActivityConcentrations(doublereal* const conc); + //! Return the charge transfer rxn Beta parameter for the ith reaction /*! * Returns the beta parameter for a charge transfer reaction. This @@ -222,7 +255,7 @@ public: /*! * This is actually the guts of the functionality of the object */ - void updateROP(); + virtual void updateROP(); //! Update properties that depend on temperature /*! @@ -414,6 +447,9 @@ public: */ int phaseStability(const size_t iphase) const; + + void determineFwdOrdersBV(ReactionData& rdata, std::vector& fwdFullorders); + protected: //! Temporary work vector of length m_kk vector_fp m_grt; @@ -663,6 +699,23 @@ protected: */ vector_int m_ctrxn_ecdf; + //! Vector of booleans indicating whether the charge transfer reaction rate constant + //! is described by an exchange current density rate constant expression + /*! + * Length is equal to the number of reactions with charge transfer coefficients, m_ctrxn[] + * + * Some reactions have zero in this list, those that don't need special treatment. + */ + std::vector m_ctrxn_ROPOrdersList_; + + //! Reaction Orders for the case where the forwards rate of progress is being calculated. + /*! + * Length is equal to the number of reactions with charge transfer coefficients, m_ctrxn[] + * + * Some reactions have zero in this list, indicating that the calculation isn't necessary. + */ + std::vector m_ctrxn_FwdOrdersList_; + //! Vector of standard concentrations /*! * Length number of kinetic species diff --git a/include/cantera/kinetics/ReactionData.h b/include/cantera/kinetics/ReactionData.h index dd767cca1..5f6a6f23d 100644 --- a/include/cantera/kinetics/ReactionData.h +++ b/include/cantera/kinetics/ReactionData.h @@ -24,6 +24,7 @@ public: validate(false), number(0), rxn_number(0), + filmResistivity(0.0), reversible(true), duplicate(false), rateCoeffType(ARRHENIUS_REACTION_RATECOEFF_TYPE), @@ -102,6 +103,13 @@ public: //! products. std::map net_stoich; + //! Film Resistivity value + /*! + * Only valid for Butler-Volmer formulations. + * Units are in ohms m2. + */ + double filmResistivity; + //! True if the current reaction is reversible. False otherwise bool reversible; diff --git a/include/cantera/kinetics/importKinetics.h b/include/cantera/kinetics/importKinetics.h index 6ba288c09..828a1bfc6 100644 --- a/include/cantera/kinetics/importKinetics.h +++ b/include/cantera/kinetics/importKinetics.h @@ -87,9 +87,20 @@ bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp, std::string default std::vector& spnum, vector_fp& stoich, vector_fp& order, const ReactionRules& rules); +//! Install Butler Volmer Orders into the forward orders array. +/*! + * Install the BV order coefficients into the fullForwardsOrders vector. + * + * @param[in] rxnNode XML node pointing to the reaction element in the xml tree. + * @param[in] kin Reference to the kinetics object to install the information into. + * @param[in] rdata Reaction Data Object containing the information about one reaction + * @param[out] fullForwardsOrders Vectors of the orders of reaction. + * + */ +void installButlerVolmerOrders(const XML_Node& rxnNode, const Kinetics& kin, const ReactionData& rdata, + std::vector& fullForwardsOrders); + //! Get non-mass-action orders for a reaction - - extern bool getOrders(const XML_Node& rxnNode, Kinetics& kin, std::string default_phase, const ReactionData& rdata, vector_fp& order, vector_fp& fullForwardsOrders, diff --git a/src/kinetics/ElectrodeKinetics.cpp b/src/kinetics/ElectrodeKinetics.cpp index 4cca1fa77..2a7161118 100644 --- a/src/kinetics/ElectrodeKinetics.cpp +++ b/src/kinetics/ElectrodeKinetics.cpp @@ -3,6 +3,7 @@ */ #include "cantera/kinetics/ElectrodeKinetics.h" +#include "cantera/thermo/SurfPhase.h" using namespace std; @@ -117,6 +118,242 @@ void ElectrodeKinetics::identifyMetalPhase() } } +//============================================================================================================================ +// virtual from InterfaceKinetics +void ElectrodeKinetics::updateROP() +{ + // evaluate rate constants and equilibrium constants at temperature and phi (electric potential) + _update_rates_T(); + // get updated activities (rates updated below) + _update_rates_C(); + + double TT = m_surf->temperature(); + double rtdf = GasConstant * TT / Faraday; + + if (m_ROP_ok) { + return; + } + // + // Copy the reaction rate coefficients, m_rfn, into m_ropf + // + copy(m_rfn.begin(), m_rfn.end(), m_ropf.begin()); + // + // Multiply by the perturbation factor + // + multiply_each(m_ropf.begin(), m_ropf.end(), m_perturb.begin()); + // + // Copy the forward rate constants to the reverse rate constants + // + copy(m_ropf.begin(), m_ropf.end(), m_ropr.begin()); + + + + // + // For reverse rates computed from thermochemistry, multiply + // the forward rates copied into m_ropr by the reciprocals of + // the equilibrium constants + // + multiply_each(m_ropr.begin(), m_ropr.end(), m_rkcn.begin()); + // + // multiply ropf by the actyivity concentration reaction orders to obtain + // the forward rates of progress. + // + m_rxnstoich.multiplyReactants(DATA_PTR(m_actConc), DATA_PTR(m_ropf)); + // + // For reversible reactions, multiply ropr by the activity concentration products + // + m_rxnstoich.multiplyRevProducts(DATA_PTR(m_actConc), DATA_PTR(m_ropr)); + // + // Fix up these calculations for cases where the above formalism doesn't hold + // + double OCV = 0.0; + for (size_t iBeta = 0; iBeta < m_beta.size(); iBeta++) { + size_t irxn = m_ctrxn[iBeta]; + + int reactionType = reactionTypes_[irxn]; + if (reactionType == BUTLERVOLMER_RXN) { + // + // Get the beta value + // + double beta = m_beta[iBeta]; + // + // OK, the reaction rate constant contains the current density rate constant calculation + // the rxnstoich calculation contained the dependence of the current density on the activity concentrations + // We finish up with the ROP calculation + // + + // + // Get the phase mole change structure + // + RxnMolChange* rmc = rmcVector[irxn]; + // + // Calculate the stoichiometric eletrons for the reaction + // This is the number of electrons that are the net products of the reaction + // + double nStoichElectrons = - rmc->m_phaseChargeChange[metalPhaseRS_]; + // + // Calculate the open circuit voltage of the reaction + // + getDeltaGibbs(0); + if (nStoichElectrons != 0.0) { + OCV = m_deltaG[irxn]/Faraday/ nStoichElectrons; + } else { + OCV = 0.0; + } + // + // Calculate the voltage of the electrode. + // + double voltage = m_phi[metalPhaseRS_] - m_phi[solnPhaseRS_]; + // + // Calculate the overpotential + // + double nu = voltage - OCV; + // + // Calculate the exchange current density + // m_ropf contains the exchange current reaction rate + // + double io = m_ropf[irxn] * nStoichElectrons; + + double exp1 = nu * nStoichElectrons * beta / rtdf; + double exp2 = - nu * nStoichElectrons * (1.0 - beta) / (rtdf); + m_ropnet[irxn] = io * (exp(exp1) - exp(exp2)); + + // Need to resurrect the forwards rate constant. + //m_ropf[irxn] = ; + m_ropr[irxn] = m_ropnet[irxn] - m_ropf[irxn]; + + } + + + } + + + + for (size_t j = 0; j != m_ii; ++j) { + m_ropnet[j] = m_ropf[j] - m_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 (size_t j = 0; j != m_ii; ++j) { + if ((m_ropr[j] > m_ropf[j]) && (m_ropr[j] > 0.0)) { + for (size_t p = 0; p < nPhases(); p++) { + if (m_rxnPhaseIsProduct[j][p]) { + if (! m_phaseExists[p]) { + m_ropnet[j] = 0.0; + m_ropr[j] = m_ropf[j]; + if (m_ropf[j] > 0.0) { + for (size_t rp = 0; rp < nPhases(); rp++) { + if (m_rxnPhaseIsReactant[j][rp]) { + if (! m_phaseExists[rp]) { + m_ropnet[j] = 0.0; + m_ropr[j] = m_ropf[j] = 0.0; + } + } + } + } + } + } + if (m_rxnPhaseIsReactant[j][p]) { + if (! m_phaseIsStable[p]) { + m_ropnet[j] = 0.0; + m_ropr[j] = m_ropf[j]; + } + } + } + } else if ((m_ropf[j] > m_ropr[j]) && (m_ropf[j] > 0.0)) { + for (size_t p = 0; p < nPhases(); p++) { + if (m_rxnPhaseIsReactant[j][p]) { + if (! m_phaseExists[p]) { + m_ropnet[j] = 0.0; + m_ropf[j] = m_ropr[j]; + if (m_ropf[j] > 0.0) { + for (size_t rp = 0; rp < nPhases(); rp++) { + if (m_rxnPhaseIsProduct[j][rp]) { + if (! m_phaseExists[rp]) { + m_ropnet[j] = 0.0; + m_ropf[j] = m_ropr[j] = 0.0; + } + } + } + } + } + } + if (m_rxnPhaseIsProduct[j][p]) { + if (! m_phaseIsStable[p]) { + m_ropnet[j] = 0.0; + m_ropf[j] = m_ropr[j]; + } + } + } + } + } + } + + m_ROP_ok = true; +} +//================================================================================================================== +// +// When the BV form is used we still need to go backwards to calculate the forward rate of progress. +// This routine does that +// +double ElectrodeKinetics::calcForwardROP_BV(size_t irxn, size_t iBeta) +{ + doublereal rt = GasConstant * thermo(0).temperature(); + doublereal rrt = 1.0/rt; + // + // Calculate gather the exchange current reaction rate constant (where does n_s appear?) + // + double iorc = m_rfn[irxn] * m_perturb[irxn]; + // + // Determine whether the reaction rate constant is in an exchange current density formulation format. + // + int iECDFormulation = m_ctrxn_ecdf[iBeta]; + + if (!iECDFormulation) { + throw CanteraError("", "not handled yet"); + } + // + // Calculate the forward chemical and modify the forward reaction rate coefficient + // + double tmp = exp(- m_beta[iBeta] * m_deltaG0[irxn] * rrt); + double tmp2 = m_ProdStanConcReac[irxn]; + tmp *= 1.0 / tmp2 / Faraday; + // + // Calculate the chemical reaction rate constant + // + double kf = iorc * tmp; + // + // Calculate the electrochemical factor + // + double eamod = m_beta[iBeta] * deltaElectricEnergy_[irxn]; + kf *= exp(- eamod * rrt); + // + // Calculate the forward rate of progress + // -> get the pointer for the orders + // + const RxnOrders* ro_fwd = m_ctrxn_FwdOrdersList_[iBeta]; + if (ro_fwd == 0) { + throw CanteraError("ElectrodeKinetics::calcForwardROP_BV()", "forward orders pointer is zero ?!?"); + } + tmp = 1.0; + const std::vector& kinSpeciesIDs = ro_fwd->kinSpeciesIDs_; + const std::vector& kinSpeciesOrders = ro_fwd->kinSpeciesOrders_; + for (size_t j = 0; j < kinSpeciesIDs.size(); j++) { + size_t k = kinSpeciesIDs[j]; + double oo = kinSpeciesOrders[j]; + tmp *= pow(m_actConc[k], oo); + } + double ropf = kf * tmp; + + return ropf; +} +//================================================================================================================== + + -//================================================================================================================== } diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index e25bcf24e..e763e689f 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -61,7 +61,12 @@ InterfaceKinetics::~InterfaceKinetics() for (size_t i = 0; i < rmcVector.size(); i++) { delete rmcVector[i]; } - + for (size_t i = 0; i < m_ctrxn_ROPOrdersList_.size(); i++) { + delete m_ctrxn_ROPOrdersList_[i]; + } + for (size_t i = 0; i < m_ctrxn_FwdOrdersList_.size(); i++) { + delete m_ctrxn_FwdOrdersList_[i]; + } } //============================================================================================================================ InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics& right) : @@ -175,6 +180,24 @@ InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right) } } + for (size_t i = 0; i < m_ctrxn_ROPOrdersList_.size(); i++) { + delete m_ctrxn_ROPOrdersList_[i]; + } + m_ctrxn_ROPOrdersList_ = right.m_ctrxn_ROPOrdersList_; + for (size_t i = 0; i < m_ctrxn_ROPOrdersList_.size(); i++) { + RxnOrders* ro = right.m_ctrxn_ROPOrdersList_[i]; + m_ctrxn_ROPOrdersList_[i] = new RxnOrders(*ro); + } + + for (size_t i = 0; i < m_ctrxn_FwdOrdersList_.size(); i++) { + delete m_ctrxn_FwdOrdersList_[i]; + } + m_ctrxn_FwdOrdersList_ = right.m_ctrxn_FwdOrdersList_; + for (size_t i = 0; i < m_ctrxn_FwdOrdersList_.size(); i++) { + RxnOrders* ro = right.m_ctrxn_FwdOrdersList_[i]; + m_ctrxn_FwdOrdersList_[i] = new RxnOrders(*ro); + } + return *this; } @@ -390,6 +413,11 @@ void InterfaceKinetics::getEquilibriumConstants(doublereal* kc) //=========================================================================================================== /* * values needed to convert from exchange current density to surface reaction rate. + * Calculate: + * - m_StandardConc[] + * - m_ProdStandConcReac[] + * - m_deltaG0[] + * - m_mu0[] */ void InterfaceKinetics::updateExchangeCurrentQuantities() { @@ -419,7 +447,6 @@ void InterfaceKinetics::updateExchangeCurrentQuantities() m_ProdStanConcReac[i] = 1.0; } m_rxnstoich.multiplyReactants(DATA_PTR(m_StandardConc), DATA_PTR(m_ProdStanConcReac)); - } //=========================================================================================================== void InterfaceKinetics::getCreationRates(doublereal* cdot) @@ -476,25 +503,31 @@ void InterfaceKinetics::applyVoltageKfwdCorrection(doublereal* const kf) #endif for (size_t i = 0; i < m_beta.size(); i++) { size_t irxn = m_ctrxn[i]; - eamod = m_beta[i] * deltaElectricEnergy_[irxn]; - if (eamod != 0.0) { + // + // If we calculate the BV form directly, we don't add the voltage correction to the + // forward reaction rate constants. + // + if (m_ctrxn_BVform[i] == 0) { + eamod = m_beta[i] * deltaElectricEnergy_[irxn]; + if (eamod != 0.0) { #ifdef DEBUG_KIN_MODE - ea = GasConstant * m_E[irxn]; - if (eamod + ea < 0.0) { - writelog("Warning: act energy mod too large!\n"); - writelog(" Delta phi = "+fp2str(deltaElectricEnergy_[irxn]/Faraday)+"\n"); - writelog(" Delta Ea = "+fp2str(eamod)+"\n"); - writelog(" Ea = "+fp2str(ea)+"\n"); - for (n = 0; n < np; n++) { - writelog("Phase "+int2str(n)+": phi = " - +fp2str(m_phi[n])+"\n"); - } - } + ea = GasConstant * m_E[irxn]; + if (eamod + ea < 0.0) { + writelog("Warning: act energy mod too large!\n"); + writelog(" Delta phi = "+fp2str(deltaElectricEnergy_[irxn]/Faraday)+"\n"); + writelog(" Delta Ea = "+fp2str(eamod)+"\n"); + writelog(" Ea = "+fp2str(ea)+"\n"); + for (n = 0; n < np; n++) { + writelog("Phase "+int2str(n)+": phi = " + +fp2str(m_phi[n])+"\n"); + } + } #endif - doublereal rt = GasConstant*thermo(0).temperature(); - doublereal rrt = 1.0/rt; - kf[irxn] *= exp(-eamod*rrt); - } + doublereal rt = GasConstant*thermo(0).temperature(); + doublereal rrt = 1.0/rt; + kf[irxn] *= exp(-eamod*rrt); + } + } } } //================================================================================================================== @@ -525,8 +558,10 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con int iECDFormulation = m_ctrxn_ecdf[i]; if (iECDFormulation) { // - // If the BV form is to be converted into the normal form then we go through this process - // If it isn't to be converted, then we don't go through this process + // If the BV form is to be converted into the normal form then we go through this process. + // If it isn't to be converted, then we don't go through this process. + // + // We need to have the straight chemical reaction rate constant to come out of this calculation. // if (m_ctrxn_BVform[i] == 0) { // @@ -537,7 +572,12 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con tmp *= 1.0 / tmp2 / Faraday; kfwd[irxn] *= tmp; } + // + // If BVform is nonzero we don't need to do anything. + // } else { + // + // kfwd[] is the chemical reaction rate constant // // If we are to calculate the BV form directly, then we will do the reverse. // We will calculate the exchange current density formulation here and @@ -546,7 +586,7 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con if (m_ctrxn_BVform[i] != 0) { // // Calculate the term and modify the forward reaction rate constant so that - // it's in exchange current density formulation format + // it's in the exchange current density formulation format // double tmp = exp(m_beta[i] * m_deltaG0[irxn] * rrt); double tmp2 = m_ProdStanConcReac[irxn]; @@ -930,16 +970,20 @@ void InterfaceKinetics::addReaction(ReactionData& r) void InterfaceKinetics::addElementaryReaction(ReactionData& rdata) { - // install rate coeff calculator + // + // install rate coefficient calculator + // vector_fp& rp = rdata.rateCoeffParameters; size_t ncov = rdata.cov.size(); + // + // Turn on the global flag indicating surface coverage dependence + // if (ncov > 3) { m_has_coverage_dependence = true; } for (size_t m = 0; m < ncov; m++) { rp.push_back(rdata.cov[m]); } - // // Find out the reaction type // @@ -981,6 +1025,8 @@ void InterfaceKinetics::addElementaryReaction(ReactionData& rdata) } else { m_ctrxn_ecdf.push_back(0); } + m_ctrxn_ROPOrdersList_.push_back(0); + m_ctrxn_FwdOrdersList_.push_back(0); } // add constant term to rate coeff value vector @@ -1029,14 +1075,15 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata) * Change the reaction rate coefficient type back to its original value */ rdata.rateCoeffType = reactionRateCoeffType_orig; - - // store activation energy + // + // Store activation energy + // m_E.push_back(rdata.rateCoeffParameters[2]); // // Add the reaction into the list of electrochemical extras // - if (rdata.beta > 0.0) { + if (rdata.beta > 0.0 || 1) { m_has_electrochem_rxns = true; m_beta.push_back(rdata.beta); // Push back the id of the reaction @@ -1049,6 +1096,27 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata) } else { m_ctrxn_ecdf.push_back(0); } + + if (rdata.forwardFullOrder_.size() > 0) { + RxnOrders* ro = new RxnOrders(); + ro->fill(rdata.forwardFullOrder_); + m_ctrxn_ROPOrdersList_.push_back(ro); + m_ctrxn_FwdOrdersList_.push_back(0); + // + // + // Fill in the Fwd Orders dependence here for B-V reactions + // + if (rdata.reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN || rdata.reactionType == BUTLERVOLMER_RXN) { + std::vector fwdFullorders(m_kk, 0.0); + determineFwdOrdersBV(rdata, fwdFullorders); + RxnOrders* ro = new RxnOrders(); + ro->fill(rdata.forwardFullOrder_); + m_ctrxn_FwdOrdersList_[m_ii] = ro; + } + } else { + m_ctrxn_ROPOrdersList_.push_back(0); + m_ctrxn_FwdOrdersList_.push_back(0); + } } // add constant term to rate coeff value vector @@ -1303,9 +1371,19 @@ int InterfaceKinetics::phaseStability(const size_t iphase) const return m_phaseIsStable[iphase]; } //================================================================================================================== -int InterfaceKinetics::reactionType(size_t i) const +doublereal InterfaceKinetics::reactantStoichCoeff(size_t kSpecKin, size_t irxn) const { - return reactionType_[i]; + return m_rrxn[kSpecKin][irxn]; +} +//================================================================================================================== +doublereal InterfaceKinetics::productStoichCoeff(size_t kSpecKin, size_t irxn) const +{ + return m_prxn[kSpecKin][irxn]; +} +//================================================================================================================== +int InterfaceKinetics::reactionType(size_t irxn) const +{ + return reactionType_[irxn]; } //================================================================================================================== void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStable) @@ -1334,6 +1412,41 @@ void InterfaceKinetics::registerReaction(size_t rxnNumber, int type, size_t loc) m_index[rxnNumber] = std::pair(type, loc); } //================================================================================================================== +// +void InterfaceKinetics::determineFwdOrdersBV(ReactionData& rdata, std::vector& fwdFullorders) +{ + // + // Start out with the full ROP orders vector. + // This vector will have the BV exchange current density orders in it. + // + fwdFullorders = rdata.forwardFullOrder_; + // + // forward and reverse beta values + // + double betaf = rdata.beta; + double betar = 1.0 - betaf; + // + // Loop over the reactants doing away the BV terms. + // This should leave the reactant terms only, even if they are non-mass action. + // + for (size_t j = 0; j < rdata.reactants.size(); j++) { + size_t kkin = rdata.reactants[j]; + double oo = rdata.rstoich[kkin]; + fwdFullorders[kkin] += betaf * oo; + if (abs(fwdFullorders[kkin]) < 0.00001) { + fwdFullorders[kkin] = 0.0; + } + } + for (size_t j = 0; j < rdata.products.size(); j++) { + size_t kkin = rdata.products[j]; + double oo = rdata.pstoich[kkin]; + fwdFullorders[kkin] -= betaf * oo; + if (abs(fwdFullorders[kkin]) < 0.00001) { + fwdFullorders[kkin] = 0.0; + } + } +} +//================================================================================================================== void EdgeKinetics::finalize() { deltaElectricEnergy_.resize(std::max(m_ii, 1)); @@ -1366,6 +1479,44 @@ void EdgeKinetics::finalize() m_finalized = true; } - +//================================================================================================================== +RxnOrders::RxnOrders() +{ +} +//================================================================================================================== +RxnOrders::~RxnOrders() +{ +} +//================================================================================================================== +RxnOrders::RxnOrders(const RxnOrders& right) : + kinSpeciesIDs_(right.kinSpeciesIDs_), + kinSpeciesOrders_(right.kinSpeciesOrders_) +{ +} +//================================================================================================================== +RxnOrders& RxnOrders::operator=(const RxnOrders& right) +{ + if (this == &right) { + return *this; + } + kinSpeciesIDs_ = right.kinSpeciesIDs_; + kinSpeciesOrders_ = right.kinSpeciesOrders_; + return *this; +} +//================================================================================================================== +int RxnOrders::fill(const std::vector& fullForwardOrders) +{ + int nzeroes = 0; + kinSpeciesIDs_.clear(); + kinSpeciesOrders_.clear(); + for (size_t k = 0; k < fullForwardOrders.size(); ++k) { + if (fullForwardOrders[k] != 0.0) { + kinSpeciesIDs_.push_back(k); + kinSpeciesOrders_.push_back(fullForwardOrders[k]); + ++nzeroes; + } + } + return nzeroes; +} //================================================================================================================== } diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index d8e74ef0c..6c91119c4 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -243,6 +243,56 @@ bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp, return true; } //==================================================================================================================== +// +// Install the BV order coefficients into the fullForwardsOrders vector. +// +void installButlerVolmerOrders(const XML_Node& rxnNode, const Kinetics& kin, const ReactionData& rdata, + std::vector& fullForwardsOrders) +{ + const std::vector& reactants = rdata.reactants; + const std::vector& products = rdata.products; + const std::vector& rstoich = rdata.rstoich; + const std::vector& pstoich = rdata.pstoich; + // + // Gather the number of species in the kinetics object and resize fullForwardsOrders + // + size_t nsp = kin.nTotalSpecies(); + fullForwardsOrders.resize(nsp, 0.0); + // + // Ok first thing to do is get the electrochemical transfer coefficient + // since the order depend on the value. + // Also, if we don't find one, then it's an error. Zero is an acceptable value. + // Beta below 0 or greater than 1 are probably not good. + // + double beta = -10.0; + if (rxnNode.hasChild("rateCoeff")) { + XML_Node& rc = rxnNode.child("rateCoeff"); + if (rc.hasChild("electrochem")) { + XML_Node& eb = rc.child("electrochem"); + string sbeta = eb["beta"]; + beta = fpValueCheck(sbeta); + } + } + if (beta == -10.0) { + throw CanteraError("installButlerVolmerOrders()", + "ButlerVolmerOrders model requested but no electrochem beta input"); + } + double betar = 1.0 - beta; + for (size_t k = 0; k < nsp; k++) { + fullForwardsOrders[k] = 0.0; + } + for (size_t n = 0; n < reactants.size(); n++) { + size_t k = reactants[n]; + double fac = rstoich[n]; + fullForwardsOrders[k] += fac * betar; + } + for (size_t n = 0; n < products.size(); n++) { + size_t k = products[n]; + double fac = pstoich[n]; + fullForwardsOrders[k] += fac * beta; + } +} +//==================================================================================================================== // Fill in the fullForwardsOrders array for a specific reaction /* * rxnNode XML node for the reaction @@ -825,12 +875,13 @@ bool rxninfo::installReaction(int iRxn, const XML_Node& rxnNode, Kinetics& kin, } } // - // get the reactant and their stoichiometries + // Get the reactant and their stoichiometries // bool ok = getReagents(rxnNode, kin, 1, default_phase, rdata.reactants, rdata.rstoich, rdata.rorder, rules); - + // // Get the products. We store the id of products in rdata.products + // ok = ok && getReagents(rxnNode, kin, -1, default_phase, rdata.products, rdata.pstoich, rdata.porder, rules); @@ -863,9 +914,26 @@ bool rxninfo::installReaction(int iRxn, const XML_Node& rxnNode, Kinetics& kin, } rdata.global = true; } - // - // Fill in the forwardFullOrder_ array + // For Butler Volmer reactions, we'll install the orders for the exchange current into the + // forwardFullOrders array. It may be altered by the getOrders function below. + // + if (rdata.reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN || rdata.reactionType == BUTLERVOLMER_RXN) { + if (! rdata.reversible) { + throw CanteraError("installReaction()", "a Butler-Volmer rxn must be reversible"); + } + installButlerVolmerOrders(rxnNode, kin, rdata, rdata.forwardFullOrder_); + // + // For Butler Volmer reactions, a common addition to the formulation is to add an electrical resistance + // to the formulation. The resistance modifies the electrical current flow in both directions + // + if (rxnNode.hasChild("filmResistivity")) { + XML_Node& fNode = rxnNode.child("filmResistivity"); + rdata.filmResistivity = fpValueCheck( fNode() ); + } + } + // + // Fill in the forwardFullOrder_ array // if (rxnNode.hasChild("orders")) { ok = getOrders(rxnNode, kin, default_phase, rdata,