From 021c8f9880115a77da729089ca171b6fab8ffc9f Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Wed, 3 Sep 2014 04:55:27 +0000 Subject: [PATCH] InterfaceKinetics rewrite -> addition of BV and Affinity reactions, next interation. --- include/cantera/kinetics/EdgeKinetics.h | 2 +- include/cantera/kinetics/ElectrodeKinetics.h | 34 +- include/cantera/kinetics/InterfaceKinetics.h | 7 +- include/cantera/kinetics/Kinetics.h | 6 +- src/kinetics/ElectrodeKinetics.cpp | 613 ++++++++++++++++++- src/kinetics/InterfaceKinetics.cpp | 101 ++- 6 files changed, 717 insertions(+), 46 deletions(-) diff --git a/include/cantera/kinetics/EdgeKinetics.h b/include/cantera/kinetics/EdgeKinetics.h index a3b93fd9a..a70cef951 100644 --- a/include/cantera/kinetics/EdgeKinetics.h +++ b/include/cantera/kinetics/EdgeKinetics.h @@ -45,7 +45,7 @@ public: return cEdgeKinetics; } - // defined in InterfaceKinetics.cpp + virtual void finalize(); }; } diff --git a/include/cantera/kinetics/ElectrodeKinetics.h b/include/cantera/kinetics/ElectrodeKinetics.h index d9159bf6f..86afdd8fd 100644 --- a/include/cantera/kinetics/ElectrodeKinetics.h +++ b/include/cantera/kinetics/ElectrodeKinetics.h @@ -74,26 +74,50 @@ public: */ virtual void updateROP(); + virtual void determineFwdOrdersBV(ReactionData& rdata, std::vector& fwdFullorders); //void addGlobalReaction(ReactionData& r); - double calcForwardROP_BV(size_t irxn, size_t iBeta); + double calcForwardROP_BV(size_t irxn, size_t iBeta, double ioc, double nStoich, double nu, doublereal ioNet); + double calcForwardROP_BV_NoAct(size_t irxn, size_t iBeta, double ioc, double nStoich, double nu, doublereal ioNet); + + bool getExchangeCurrentDensityFormulation(size_t irxn, doublereal& nStoich, doublereal& OCV, doublereal& io, + doublereal& overPotential, doublereal& beta, doublereal& resistance); + + //! Calculate the open circuit voltage of a given reaction + /*! + * If the reaction has no electron transport, then return 0.0 + * + * @param irxn Reaction id + */ + double openCircuitVoltage(size_t irxn); + + double calcCurrentDensity(double nu, double nStoich, double io, double beta, double temp, doublereal resistivity = 0.0) const; + + double solveCurrentRes(doublereal nu, doublereal nStoich, doublereal ioc, doublereal beta, doublereal temp, + doublereal resistivity = 0.0, int iprob = 0) const; + + //! Prepare the class for the addition of reactions + /*! + * (virtual from Kinetics) + * We determine the metal phase and solution phase here + */ + virtual void init(); protected: - //! Index of the metal phase in the list of phases for this kinetics object + //! Index of the metal phase in the list of phases for this kinetics object. This is the electron phase. 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 size_t solnPhaseRS_; //! 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 a6a4aed82..2c0e1b752 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -8,8 +8,11 @@ #ifndef CT_IFACEKINETICS_H #define CT_IFACEKINETICS_H +#include "cantera/thermo/mix_defs.h" #include "Kinetics.h" #include "cantera/kinetics/RxnMolChange.h" + +#include "cantera/base/utilities.h" #include "RateCoeffMgr.h" #include "ReactionStoichMgr.h" @@ -445,7 +448,7 @@ public: int phaseStability(const size_t iphase) const; - void determineFwdOrdersBV(ReactionData& rdata, std::vector& fwdFullorders); + virtual void determineFwdOrdersBV(ReactionData& rdata, std::vector& fwdFullorders); protected: //! Temporary work vector of length m_kk @@ -713,6 +716,8 @@ protected: */ std::vector m_ctrxn_FwdOrdersList_; + std::vector m_ctrxn_resistivity_; + //! Vector of standard concentrations /*! * Length number of kinetic species diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index c06e7af0b..13a1a78b7 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -434,6 +434,10 @@ public: * units in array kc, which must be dimensioned at least as large as the * total number of reactions. * + * \f[ + * Kc_i = exp [ \Delta G_{ss,i} ] prod(Cs_k) exp(\sum_k \nu_{k,i} F \phi_n) ] + * \f] + * * @param kc Output vector containing the equilibrium constants. * Length: m_ii. */ @@ -446,7 +450,7 @@ public: * property values \f$ z_k, k = 1, \dots, K \f$, return the * array of reaction values * \f[ - * \Delta Z_i = \sum_k \nu_{k,i} z_k, i = 1, \dots, I. + * \Delta Z_i = \sum_k \nu_{k,i} z_k, i = 1, \dots, I. * \f] * For example, if this method is called with the array of * standard-state molar Gibbs free energies for the species, diff --git a/src/kinetics/ElectrodeKinetics.cpp b/src/kinetics/ElectrodeKinetics.cpp index 9a80439f0..7ecdba98b 100644 --- a/src/kinetics/ElectrodeKinetics.cpp +++ b/src/kinetics/ElectrodeKinetics.cpp @@ -6,6 +6,8 @@ #include "cantera/thermo/SurfPhase.h" #include "cantera/base/utilities.h" +#include + using namespace std; namespace Cantera @@ -14,7 +16,6 @@ namespace Cantera ElectrodeKinetics::ElectrodeKinetics(thermo_t* thermo) : InterfaceKinetics(thermo), metalPhaseRS_(npos), - electronPhaseRS_(npos), solnPhaseRS_(npos), kElectronRS_(npos) { @@ -48,7 +49,6 @@ ElectrodeKinetics& ElectrodeKinetics::operator=(const ElectrodeKinetics& right) InterfaceKinetics::operator=(right); metalPhaseRS_ = right.metalPhaseRS_; - electronPhaseRS_ = right.electronPhaseRS_; solnPhaseRS_ = right.solnPhaseRS_; kElectronRS_ = right.kElectronRS_; @@ -67,12 +67,12 @@ Kinetics* ElectrodeKinetics::duplMyselfAsKinetics(const std::vector & return iK; } //============================================================================================================================ -//==================================================================================================================== -// Identify the metal phase and the electrons species +// Identify the metal phase and the electron species void ElectrodeKinetics::identifyMetalPhase() { metalPhaseRS_ = npos; - kElectronRS_ = -1; + kElectronRS_ = npos; + solnPhaseRS_ = npos; size_t np = nPhases(); // // Identify the metal phase as the phase with the electron species (element index of 1 for element E @@ -105,6 +105,8 @@ void ElectrodeKinetics::identifyMetalPhase() // Identify the solution phase as a 3D phase, with nonzero phase charge change // in at least one reaction // + /* + * Haven't filled in reactions yet when this is called, unlike previous treatment. if (iph != metalPhaseRS_) { for (size_t i = 0; i < m_ii; i++) { RxnMolChange* rmc = rmcVector[i]; @@ -116,8 +118,38 @@ void ElectrodeKinetics::identifyMetalPhase() } } } + */ + // + // New method is to find the first multispecies 3D phase with charged species as the solution phase + // + if (iph != metalPhaseRS_) { + ThermoPhase& tp =*( m_thermo[iph]); + size_t nsp = tp.nSpecies(); + size_t nd = tp.nDim(); + if (nd == 3 && nsp > 1) { + for (size_t k = 0; k < nsp; k++) { + if (tp.charge(k) != 0.0) { + solnPhaseRS_ = iph; + string ss = tp.name(); + // cout << "solution phase = "<< ss << endl; + break; + } + } + } + } } + // + // Right now, if we don't find an electron phase, let's error exit to avoid complications + // + if (metalPhaseRS_ == npos) { + throw CanteraError("ElectrodeKinetics::identifyMetalPhase()", + "Can't find electron phase -> treating this as an error right now"); + } + if (solnPhaseRS_ == npos) { + throw CanteraError("ElectrodeKinetics::identifyMetalPhase()", + "Can't find solution phase -> treating this as an error right now"); + } } //============================================================================================================================ // virtual from InterfaceKinetics @@ -148,7 +180,6 @@ void ElectrodeKinetics::updateROP() 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 @@ -182,7 +213,11 @@ void ElectrodeKinetics::updateROP() // the rxnstoich calculation contained the dependence of the current density on the activity concentrations // We finish up with the ROP calculation // - + int iECDFormulation = m_ctrxn_ecdf[iBeta]; + if (iECDFormulation == 0) { + throw CanteraError(" ElectrodeKinetics::updateROP()", + "Straight kfwrd with BUTLERVOLMER_RXN not handled yet"); + } // // Get the phase mole change structure // @@ -209,22 +244,151 @@ void ElectrodeKinetics::updateROP() // Calculate the overpotential // double nu = voltage - OCV; + + // + // Find the product of the standard concentrations for ROP orders that we used above + // + const RxnOrders* ro_rop = m_ctrxn_ROPOrdersList_[iBeta]; + if (ro_rop == 0) { + throw CanteraError("ElectrodeKinetics::", "ROP orders pointer is zero ?!?"); + } + double tmp2 = 1.0; + const std::vector& kinSpeciesIDs = ro_rop->kinSpeciesIDs_; + const std::vector& kinSpeciesOrders = ro_rop->kinSpeciesOrders_; + for (size_t j = 0; j < kinSpeciesIDs.size(); j++) { + size_t k = kinSpeciesIDs[j]; + double oo = kinSpeciesOrders[j]; + tmp2 *= pow(m_StandardConc[k], oo); + } + // + // Now have to divide this to get rid of standard concentrations. We should + // have used just the activities in the m_rxnstoich.multiplyReactants(DATA_PTR(m_actConc), DATA_PTR(m_ropf)); + // calculation above! + // That is because the exchange current density rate constants have the correct units in the first place. + // + m_ropf[irxn] /= tmp2; // // Calculate the exchange current density // m_ropf contains the exchange current reaction rate // - double io = m_ropf[irxn] * nStoichElectrons; - + double ioc = m_ropf[irxn] * nStoichElectrons; + // + // Add in the film resistance here + // + double resist = m_ctrxn_resistivity_[iBeta]; double exp1 = nu * nStoichElectrons * beta / rtdf; double exp2 = - nu * nStoichElectrons * (1.0 - beta) / (rtdf); - m_ropnet[irxn] = io * (exp(exp1) - exp(exp2)); + double io = ioc * (exp(exp1) - exp(exp2)); - // Need to resurrect the forwards rate constant. - //m_ropf[irxn] = ; - m_ropr[irxn] = m_ropnet[irxn] - m_ropf[irxn]; + if (resist != 0.0) { + io = solveCurrentRes(nu, nStoichElectrons, ioc, beta, TT, resist, 0); + } + + m_ropnet[irxn] = io / (Faraday * nStoichElectrons); + // + // Need to resurrect the forwards rate of progress -> there is some need to + // calculate each direction individually + // + m_ropf[irxn] = calcForwardROP_BV(irxn, iBeta, ioc, nStoichElectrons, nu, io); + // + // Calculate the reverse rate of progress from the difference + // + m_ropr[irxn] = m_ropf[irxn] - m_ropnet[irxn]; - } + } else if (reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_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 + // + int iECDFormulation = m_ctrxn_ecdf[iBeta]; + if (iECDFormulation == 0) { + throw CanteraError("ElectrodeKinetics::updateROP()", + "Straight kfwrd with BUTLERVOLMER_NOACTIVITYCOEFFS_RXN not handled yet"); + } + // + // 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; + // + // Unfortunately, we really need to recalculate everything from almost scratch + // for this case, since it widely diverges from the thermo norm. + // + // Start with the exchange current reaction rate constant, which should + // be located in m_rfn[]. + // + double ioc = m_rfn[irxn] * nStoichElectrons * m_perturb[irxn]; + // + // Now we need th mole fraction vector and we need the RxnOrders vector. + // + const RxnOrders* ro_fwd = m_ctrxn_ROPOrdersList_[iBeta]; + if (ro_fwd == 0) { + throw CanteraError("ElectrodeKinetics::calcForwardROP_BV()", "forward orders pointer is zero ?!?"); + } + double tmp = 1.0; + double mfS = 0.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 ks = kinSpeciesIDs[j]; + thermo_t& th = speciesPhase(ks); + size_t n = speciesPhaseIndex(ks); + size_t klocal = ks - m_start[n]; + mfS = th.moleFraction(klocal); + + double oo = kinSpeciesOrders[j]; + tmp *= pow(mfS, oo); + } + ioc *= tmp; + // + // Add in the film resistance here, later + // + double resist = m_ctrxn_resistivity_[iBeta]; + double exp1 = nu * nStoichElectrons * beta / rtdf; + double exp2 = - nu * nStoichElectrons * (1.0 - beta) / (rtdf); + double io = ioc * (exp(exp1) - exp(exp2)); + if (resist != 0.0) { + io = solveCurrentRes(nu, nStoichElectrons, ioc, beta, TT, resist, 0); + } + + m_ropnet[irxn] = io / (Faraday * nStoichElectrons); + // + // Need to resurrect the forwards rate of progress -> there is some need to + // calculate each direction individually + // + m_ropf[irxn] = calcForwardROP_BV_NoAct(irxn, iBeta, ioc, nStoichElectrons, nu, io); + // + // Calculate the reverse rate of progress from the difference + // + m_ropr[irxn] = m_ropf[irxn] - m_ropnet[irxn]; + } } @@ -299,13 +463,158 @@ void ElectrodeKinetics::updateROP() } //================================================================================================================== // +// This version of takes the electrons out of the reaction rate expression +// (note: with proper specification of the phase, this shouldn't make a numerical difference (power of 1). +// But it certainly is a complication and unneeded work) +// (TODO: probably can take stoichiometric solids out of the reaction order expression as well. +// They all contribute powers of 1 as well) +// +void ElectrodeKinetics::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[j]; + if (kkin != kElectronRS_) { + fwdFullorders[kkin] += betaf * oo; + if (abs(fwdFullorders[kkin]) < 0.00001) { + fwdFullorders[kkin] = 0.0; + } + } else { + fwdFullorders[kkin] = 0.0; + } + } + for (size_t j = 0; j < rdata.products.size(); j++) { + size_t kkin = rdata.products[j]; + double oo = rdata.pstoich[j]; + if (kkin != kElectronRS_) { + fwdFullorders[kkin] -= betaf * oo; + if (abs(fwdFullorders[kkin]) < 0.00001) { + fwdFullorders[kkin] = 0.0; + } + } else { + fwdFullorders[kkin] = 0.0; + } + } +} +//================================================================================================================== +// // 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) -{ +double ElectrodeKinetics::calcForwardROP_BV(size_t irxn, size_t iBeta, double ioc, double nStoich, double nu, doublereal ioNet) +{ + double ropf; 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]; + doublereal beta = m_beta[iBeta]; + +#ifdef DEBUG_MODE + // + // 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 + // + const RxnOrders* ro_fwd = m_ctrxn_FwdOrdersList_[iBeta]; + if (ro_fwd == 0) { + throw CanteraError("ElectrodeKinetics::calcForwardROP_BV()", "forward orders pointer is zero ?!?"); + } + double tmp = exp(- m_beta[iBeta] * m_deltaG0[irxn] * rrt); + double tmp2 = 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]; + tmp2 *= pow(m_StandardConc[k], oo); + } + + //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 + // + tmp = 1.0; + + for (size_t j = 0; j < kinSpeciesIDs.size(); j++) { + size_t k = kinSpeciesIDs[j]; + double oo = kinSpeciesOrders[j]; + tmp *= pow(m_actConc[k], oo); + } + ropf = kf * tmp; +#endif + // + // Now calculate ropf in a separate but equivalent way. + // totally equivalent way if resistivity is zero, should be equal (HKM -> Proved exactly in one case) + // + double iof = ioc; + double resistivity = m_ctrxn_resistivity_[iBeta]; + if (fabs(resistivity * ioNet) > fabs(nu)) { + ioNet = nu / resistivity; + } + if (nStoich > 0.0) { + double exp1 = nStoich * Faraday * beta * (nu - resistivity * ioNet)/ (rt); + iof *= exp(exp1); + } else { +#ifdef DEBUG_MODE + if (ioc > 0) { + throw CanteraError(" ", "ioc should be less than zero here"); + } +#endif + double exp2 = -nu * nStoich * Faraday * (1.0 - beta) / (rt); + iof = ioc * ( - exp(exp2)); + } + ropf = iof / ( Faraday * nStoich); + + return ropf; +} +//================================================================================================================== +// +// 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_NoAct(size_t irxn, size_t iBeta, double ioc, double nStoich, double nu, + doublereal ioNet) +{ + doublereal TT = thermo(0).temperature(); + doublereal rt = GasConstant * TT; + //doublereal rrt = 1.0/rt; + doublereal beta = m_beta[iBeta]; + +/* // // Calculate gather the exchange current reaction rate constant (where does n_s appear?) // @@ -320,9 +629,10 @@ double ElectrodeKinetics::calcForwardROP_BV(size_t irxn, size_t iBeta) } // // Calculate the forward chemical and modify the forward reaction rate coefficient + // (we don't use standard concentrations at all here); // double tmp = exp(- m_beta[iBeta] * m_deltaG0[irxn] * rrt); - double tmp2 = m_ProdStanConcReac[irxn]; + double tmp2 = 1.0; tmp *= 1.0 / tmp2 / Faraday; // // Calculate the chemical reaction rate constant @@ -345,16 +655,273 @@ double ElectrodeKinetics::calcForwardROP_BV(size_t irxn, size_t iBeta) 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]; + + size_t ks = kinSpeciesIDs[j]; + thermo_t& th = speciesPhase(ks); + size_t n = speciesPhaseIndex(ks); + size_t klocal = ks - m_start[n]; + double mfS = th.moleFraction(klocal); double oo = kinSpeciesOrders[j]; - tmp *= pow(m_actConc[k], oo); + tmp *= pow(mfS, oo); } double ropf = kf * tmp; - +*/ +/* + if (nStoich > 0) { + double ropf = ioc / ( Faraday * nStoich); + double exp1 = nu * nStoich * Faraday * beta / (rt); + ropf *= exp(exp1); + } else { + double ropf = ioc / ( Faraday * nStoich); + double exp1 = nu * nStoich * Faraday * beta / (rt); + ropf *= exp(exp1); + } +*/ + // + // With all of the thermo issues, I'm thinking this is the best we can do + // (it certainly maintains the forward and reverse rates of progress as being positive) + // + double iof = ioc; + double resistivity = m_ctrxn_resistivity_[iBeta]; + if (fabs(resistivity * ioNet) > fabs(nu)) { + ioNet = nu / resistivity; + } + if (nStoich > 0) { + double exp1 = nStoich * Faraday * beta * (nu - resistivity * ioNet)/ (rt); + iof *= exp(exp1); + } else { +#ifdef DEBUG_MODE + if (ioc > 0) { + throw CanteraError(" ", "ioc should be less than zero here"); + } +#endif + double exp2 = -nu * nStoich * Faraday * (1.0 - beta) / (rt); + iof = ioc * ( - exp(exp2)); + } + double ropf = iof / ( Faraday * nStoich); return ropf; } //================================================================================================================== - - - +double ElectrodeKinetics::openCircuitVoltage(size_t irxn) +{ + // + // Calculate deltaG for all reactions + // + getDeltaGibbs(0); + // + // Look up the net number of electrons that are products. + // + RxnMolChange* rmc = rmcVector[irxn]; + double nStoichElectrons = - rmc->m_phaseChargeChange[metalPhaseRS_]; + double OCV = 0.0; + if (nStoichElectrons != 0.0) { + OCV = m_deltaG[irxn] / Faraday / nStoichElectrons; + } + return OCV; +} +//================================================================================================================== +// +// Returns the local exchange current density formulation parameters +// +bool ElectrodeKinetics:: +getExchangeCurrentDensityFormulation(size_t irxn, + doublereal& nStoichElectrons, doublereal& OCV, doublereal& io, + doublereal& overPotential, doublereal& beta, + doublereal& resistivity) +{ + size_t iBeta = npos; + beta = 0.0; + // + // Add logic to handle other reaction types -> return 0 if formulation isn't compatible + // + + // evaluate rate constants and equilibrium constants at temperature and phi (electric potential) + _update_rates_T(); + // get updated activities (rates updated below) + _update_rates_C(); + + updateExchangeCurrentQuantities(); + + RxnMolChange* rmc = rmcVector[irxn]; + // could also get this from reactant and product stoichiometry, maybe + nStoichElectrons = - rmc->m_phaseChargeChange[metalPhaseRS_]; + + OCV = 0.0; + + getDeltaGibbs(0); + + if (nStoichElectrons != 0.0) { + OCV = m_deltaG[irxn] / Faraday / nStoichElectrons; + } + + for (size_t i = 0; i < m_ctrxn.size(); i++) { + if (m_ctrxn[i] == irxn) { + iBeta = i; + break; + } + } + beta = m_beta[iBeta]; + + doublereal rt = GasConstant*thermo(0).temperature(); + + + double mG0 = m_deltaG0[irxn]; + int reactionType = reactionTypes_[irxn]; + + // + // Start with the forward reaction rate + // + double iO = m_rfn[irxn] * m_perturb[irxn]; + int iECDFormulation = m_ctrxn_ecdf[iBeta]; + if (! iECDFormulation) { + iO = m_rfn[irxn] * Faraday * nStoichElectrons; + if (beta > 0.0) { + double fac = exp(mG0 / (rt)); + iO *= pow(fac, beta); + // Need this step because m_rfn includes the inverse of this term, while the formulas + // only use the chemical reaction rate constant. + fac = exp( beta * deltaElectricEnergy_[irxn] / (rt)); + iO *= fac; + } + } else { + iO *= nStoichElectrons; + } + + double omb = 1.0 - beta; + if (reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN) { + const RxnOrders* ro_fwd = m_ctrxn_ROPOrdersList_[iBeta]; + if (ro_fwd == 0) { + throw CanteraError("ElectrodeKinetics::calcForwardROP_BV()", "forward orders pointer is zero ?!?"); + } + double 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 ks = kinSpeciesIDs[j]; + thermo_t& th = speciesPhase(ks); + size_t n = speciesPhaseIndex(ks); + size_t klocal = ks - m_start[n]; + double mfS = th.moleFraction(klocal); + + double oo = kinSpeciesOrders[j]; + tmp *= pow(mfS, oo); + } + iO *= tmp; + } else { + for (size_t k = 0; k < m_kk; k++) { + doublereal reactCoeff = reactantStoichCoeff(k, irxn); + doublereal prodCoeff = productStoichCoeff(k, irxn); + + if (reactCoeff != 0.0) { + iO *= pow(m_actConc[k], reactCoeff*omb); + iO *= pow(m_StandardConc[k], reactCoeff*beta); + } + if (prodCoeff != 0.0) { + iO *= pow(m_actConc[k], prodCoeff*beta); + iO /= pow(m_StandardConc[k], prodCoeff*omb); + } + } + } + io = iO; + resistivity = m_ctrxn_resistivity_[iBeta]; + + double phiMetal = m_thermo[metalPhaseRS_]->electricPotential(); + double phiSoln = m_thermo[solnPhaseRS_]->electricPotential(); + double E = phiMetal - phiSoln; + overPotential = E - OCV; + + return true; +} +//==================================================================================================================== +double ElectrodeKinetics::calcCurrentDensity(double nu, double nStoich, double ioc, double beta, double temp, + doublereal resistivity) const +{ + double exp1 = nu * nStoich * Faraday * beta / (GasConstant * temp); + double exp2 = -nu * nStoich * Faraday * (1.0 - beta) / (GasConstant * temp); + double val = ioc * (exp(exp1) - exp(exp2)); + if (resistivity > 0.0) { + val = solveCurrentRes(nu, nStoich, ioc, beta, temp, resistivity, 0); + } + return val; +} +//================================================================================================================== +void ElectrodeKinetics::init() +{ + InterfaceKinetics::init(); + identifyMetalPhase(); +} +//================================================================================================================== + +double ElectrodeKinetics::solveCurrentRes(double nu, double nStoich, doublereal ioc, doublereal beta, doublereal temp, + doublereal resistivity, int iprob) const +{ + // int nits = 0; + doublereal f, dfdi, deltai, eexp1, eexp2, exp1, exp2, icurr, deltai_damp; + doublereal nFRT = nStoich * Faraday / (GasConstant * temp); + if (iprob == 0) { + eexp1 = exp(nu * nFRT * beta); + eexp2 = exp(-nu * nFRT * (1.0 - beta)) ; + + } else { + eexp1 = exp(nu * nFRT * beta); + eexp2 = 0.0; + } + icurr = ioc * (eexp1 - eexp2); + double icurrDamp = icurr; + if (fabs(resistivity * icurr) > 0.9 * fabs(nu)) { + icurrDamp = 0.9 * nu / resistivity; + } + if (iprob == 0) { + eexp1 = exp( nFRT * beta * (nu - resistivity * icurrDamp)); + eexp2 = exp(- nFRT * (1.0 - beta) * (nu - resistivity * icurrDamp)); + } else { + eexp1 = exp( nFRT * beta * (nu - resistivity * icurrDamp)); + eexp2 = 0.0; + } + icurr = ioc * (eexp1 - eexp2); + if (fabs(resistivity * icurr) > 0.99 * fabs(nu)) { + icurr = 0.99 * nu / resistivity; + } + + do { + // nits++; + if (iprob == 0) { + exp1 = nFRT * beta * (nu - resistivity * icurr); + exp2 = - nFRT * (1.0 - beta) * (nu - resistivity * icurr); + eexp1 = exp(exp1); + eexp2 = exp(exp2); + f = icurr - ioc * (eexp1 - eexp2); + dfdi = 1.0 - ioc * eexp1 * ( - beta * nFRT * resistivity ) + + ioc * eexp2 * ( (1.0 - beta) * nFRT * resistivity ); + } else { + exp1 = nFRT * beta * (nu - resistivity * icurr); + eexp1 = exp(exp1); + f = icurr - ioc * (eexp1); + dfdi = 1.0 - ioc * eexp1 * ( - beta * nFRT * resistivity ); + } + deltai = - f / dfdi; + if (fabs(deltai) > 0.1 * fabs(icurr)) { + deltai_damp = 0.1 * deltai; + if (fabs(deltai_damp) > 0.1 * fabs(icurr)) { + deltai_damp = 0.1 * icurr * (deltai_damp / fabs(deltai_damp)); + } + } else if (fabs(deltai) > 0.01 * fabs(icurr)) { + deltai_damp = 0.3 * deltai; + } else if (fabs(deltai) > 0.001 * fabs(icurr)) { + deltai_damp = 0.5 * deltai; + } else { + deltai_damp = deltai; + } + icurr += deltai_damp; + if (fabs(resistivity * icurr) > fabs(nu)) { + icurr = 0.999 * nu / resistivity; + } + + } while((fabs(deltai/icurr)> 1.0E-14) && (fabs(deltai) > 1.0E-20)); + + // printf(" its = %d\n", nits); + + return icurr; +} +//================================================================================================================== } diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 2de81d73a..ae4539b01 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -6,9 +6,10 @@ #include "cantera/kinetics/InterfaceKinetics.h" #include "cantera/kinetics/EdgeKinetics.h" +#include "cantera/kinetics/ReactionData.h" +#include "cantera/kinetics/RateCoeffMgr.h" #include "cantera/kinetics/ImplicitSurfChem.h" #include "cantera/thermo/SurfPhase.h" -#include "cantera/base/vec_functions.h" #include @@ -221,6 +222,9 @@ void InterfaceKinetics::setElectricPotential(int n, doublereal V) //============================================================================================================================ void InterfaceKinetics::_update_rates_T() { + // + // First task is update the electrical potentials from the Phases + // _update_rates_phi(); if (m_has_coverage_dependence) { m_surf->getCoverages(DATA_PTR(m_actConc)); @@ -269,6 +273,9 @@ void InterfaceKinetics::_update_rates_phi() } } //============================================================================================================================ +// +// Updates the internal variables m_actConc and m_conc +// void InterfaceKinetics::_update_rates_C() { for (size_t n = 0; n < nPhases(); n++) { @@ -327,6 +334,12 @@ void InterfaceKinetics::updateKc() //============================================================================================================================ void InterfaceKinetics::updateMu0() { + // + // First task is update the electrical potentials from the Phases + // + _update_rates_phi(); + + updateExchangeCurrentQuantities(); /* * Get the vector of standard state electrochemical potentials for species in the Interfacial * kinetics object and store it in m_mu0[] and in m_mu0_Kc[] @@ -347,6 +360,11 @@ void InterfaceKinetics::updateMu0() //============================================================================================================================ void InterfaceKinetics::checkPartialEquil() { + // + // First task is update the electrical potentials from the Phases + // + _update_rates_phi(); + vector_fp dmu(nTotalSpecies(), 0.0); vector_fp rmu(std::max(nReactions(), 1), 0.0); if (m_nrev > 0) { @@ -424,7 +442,7 @@ void InterfaceKinetics::updateExchangeCurrentQuantities() * First collect vectors of the standard Gibbs free energies of the * species and the standard concentrations * - m_mu0 - * - m_logStandardConc + * - m_StandardConc */ size_t ik = 0; @@ -1026,6 +1044,11 @@ void InterfaceKinetics::addElementaryReaction(ReactionData& rdata) } m_ctrxn_ROPOrdersList_.push_back(0); m_ctrxn_FwdOrdersList_.push_back(0); + if (rdata.filmResistivity > 0.0) { + throw CanteraError("InterfaceKinetics::addElementaryReaction()", + "film resistivity set for elementary reaction"); + } + m_ctrxn_resistivity_.push_back(rdata.filmResistivity); } // add constant term to rate coeff value vector @@ -1085,16 +1108,29 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata) if (rdata.beta > 0.0 || 1) { m_has_electrochem_rxns = true; m_beta.push_back(rdata.beta); - // Push back the id of the reaction + // Push back the id of the reaction m_ctrxn.push_back(m_ii); - // set the default to be the normal forward / reverse calculation method - m_ctrxn_BVform.push_back(0); + // + // Specify alternative forms of the electrochemical reaction + // + if (rdata.reactionType == BUTLERVOLMER_RXN) { + m_ctrxn_BVform.push_back(1); + } else if (rdata.reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN) { + m_ctrxn_BVform.push_back(2); + } else { + // set the default to be the normal forward / reverse calculation method + m_ctrxn_BVform.push_back(0); + } if (rdata.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); } + // + // Store the film resistivity + // + m_ctrxn_resistivity_.push_back(rdata.filmResistivity); if (rdata.forwardFullOrder_.size() > 0) { RxnOrders* ro = new RxnOrders(); @@ -1109,7 +1145,7 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata) std::vector fwdFullorders(m_kk, 0.0); determineFwdOrdersBV(rdata, fwdFullorders); RxnOrders* ro = new RxnOrders(); - ro->fill(rdata.forwardFullOrder_); + ro->fill(fwdFullorders); m_ctrxn_FwdOrdersList_[m_ii] = ro; } } else { @@ -1258,12 +1294,15 @@ void InterfaceKinetics::finalize() size_t ks = reactionPhaseIndex(); if (ks == npos) throw CanteraError("InterfaceKinetics::finalize", "no surface phase is present."); + // + // Check to see that the interface routine has a dimension of 2 + // m_surf = (SurfPhase*)&thermo(ks); - if (m_surf->nDim() != 2) + if (m_surf->nDim() != 2) { throw CanteraError("InterfaceKinetics::finalize", "expected interface dimension = 2, but got dimension = " +int2str(m_surf->nDim())); - + } m_StandardConc.resize(m_kk, 0.0); m_deltaG0.resize(safe_reaction_size, 0.0); m_deltaG.resize(safe_reaction_size, 0.0); @@ -1283,6 +1322,13 @@ void InterfaceKinetics::finalize() m_ropnet.resize(1, 0.0); m_rkcn.resize(1, 0.0); } + // + // Malloc and calculate all of the quantities that go into the extra description of reactions + // + rmcVector.resize(m_ii, 0); + for (size_t i = 0; i < m_ii; i++) { + rmcVector[i] = new RxnMolChange(this, i); + } m_finalized = true; } @@ -1425,20 +1471,25 @@ void InterfaceKinetics::determineFwdOrdersBV(ReactionData& rdata, std::vector(m_ii, 1)); + // + // Note we can't call the Interface::finalize() routine because we need to check for a dimension of 1 below. + // Therefore, we have to malloc room in arrays that would normally be + // handled by the InterfaceKinetics::finalize() call. + // + Kinetics::finalize(); + + size_t safe_reaction_size = std::max(m_ii, 1); + deltaElectricEnergy_.resize(safe_reaction_size); size_t ks = reactionPhaseIndex(); if (ks == npos) throw CanteraError("EdgeKinetics::finalize", - "no edge phase is present."); + "no surface phase is present."); + // + // Check to see edge phase has a dimension of 1 + // m_surf = (SurfPhase*)&thermo(ks); - if (m_surf->nDim() != 1) + if (m_surf->nDim() != 1) { throw CanteraError("EdgeKinetics::finalize", "expected interface dimension = 1, but got dimension = " +int2str(m_surf->nDim())); + } + m_StandardConc.resize(m_kk, 0.0); + m_deltaG0.resize(safe_reaction_size, 0.0); + m_deltaG.resize(safe_reaction_size, 0.0); + + m_ProdStanConcReac.resize(safe_reaction_size, 0.0); + + if (m_thermo.size() != m_phaseExists.size()) { + throw CanteraError("InterfaceKinetics::finalize", "internal error"); + } // Guarantee that these arrays can be converted to double* even in the // special case where there are no reactions defined. @@ -1467,7 +1539,6 @@ void EdgeKinetics::finalize() m_ropnet.resize(1, 0.0); m_rkcn.resize(1, 0.0); } - // // Malloc and calculate all of the quantities that go into the extra description of reactions //