From 199cd4cebab590ca4f08395b0bb0db89f95eb210 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 5 Jun 2014 22:43:41 +0000 Subject: [PATCH] Changes in InterfaceKinetics to allow for downstream hooks to add in the capability to model arbitrarily fit open circuit potentials for intercalating electrodes. --- include/cantera/kinetics/InterfaceKinetics.h | 75 ++++++++++++++--- src/kinetics/InterfaceKinetics.cpp | 87 ++++++++++---------- 2 files changed, 104 insertions(+), 58 deletions(-) diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index 3b3258c19..a92d2fcad 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -269,19 +269,27 @@ public: void checkPartialEquil(); + //! Update the standard state chemical potentials and species equilibrium constant entries + /*! + * Virtual because it is overwritten when dealing with experimental open circuit voltage overrides + */ + virtual void updateMu0(); + size_t reactionNumber() const { return m_ii; } void addElementaryReaction(ReactionData& r); - void addGlobalReaction(const ReactionData& r); + //void addGlobalReaction(const ReactionData& r); void installReagents(const ReactionData& r); - /** - * Update the equilibrium constants in molar units for all reversible - * reactions. Irreversible reactions have their equilibrium constant set - * to zero. For reactions involving charged species the equilibrium - * constant is adjusted according to the electrostatic potential. + + //! Update the equilibrium constants and stored electrochemical potentials + //! in molar units for all reversible reactions and for all species. + /*! + * Irreversible reactions have their equilibrium constant set + * to zero. For reactions involving charged species the equilibrium + * constant is adjusted according to the electrostatic potential. */ void updateKc(); @@ -460,15 +468,31 @@ protected: */ vector_fp m_conc; - //! Vector of standard state chemical potentials + //! Vector of standard state chemical potentials for all species /*! * This vector contains a temporary vector of standard state chemical * potentials for all of the species in the kinetics object * - * Length = m_k. Units = J/kmol. + * Length = m_kk. Units = J/kmol. */ vector_fp m_mu0; + //! Vector of standard state electrochemical potentials modified by + //! a standard concentration term. + /*! + * This vector contains a temporary vector of standard state electrochemical + * potentials + RTln(Cs) for all of the species in the kinetics object + * + * In order to get the units correct for the concentration equilibrium + * constant, each species needs to have an + * RT ln(Cs) added to its contribution to the equilibrium constant + * Cs is the standard concentration for the species. Frequently, for + * solid species, Cs is equal to 1. However, for gases Cs is P/RT. + * + * Length = m_kk. Units = J/kmol. + */ + vector_fp m_mu0_Kc; + //! Vector of phase electric potentials /*! * Temporary vector containing the potential of each phase in the kinetics @@ -494,8 +518,8 @@ protected: //! Vector of raw activation energies for the reactions /*! - * units are in Kelvin - * Length is number of reactions. + * Units are in Kelvin. + * Length is number of reactions. */ vector_fp m_E; @@ -526,7 +550,18 @@ protected: //! be described by an exchange current density expression vector_int m_ctrxn_ecdf; + //! Vector of standard concentrations + /*! + * Length number of kinetic species + * units depend on the definition of the standard concentration within each phase + */ vector_fp m_StandardConc; + + //! Vector of delta G^0, the standard state gibbs free energies for each reaction + /*! + * Length is the number of reactions + * units are Joule kmol-1 + */ vector_fp m_deltaG0; vector_fp m_ProdStanConcReac; @@ -540,12 +575,22 @@ protected: //! Current temperature of the data doublereal m_temp; + //! Current log of the temperature doublereal m_logtemp; + vector_fp m_rfn; + + //! Equilibrium constant for all reactions including the voltage term + /*! + * Kc = exp(deltaG/RT) + * + * where deltaG is the electrochemical potential difference between + * products minus reactants. + */ vector_fp m_rkcn; - //! boolean indicating whether mechanism has been finalized + //! Boolean indicating whether mechanism has been finalized bool m_finalized; //! Boolean flag indicating whether any reaction in the mechanism @@ -587,7 +632,11 @@ protected: //! 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 + * this to set the ROP's so that unphysical things don't happen. + * For example, a reaction can't go in the forwards direction if a + * phase in which a reactant is present doesn't exist. Because InterfaceKinetics + * deals with intrinsic quantities only normally, nowhere else is this extrinsic + * concept introduced except here. * * length = number of phases in the object. By default all phases exist. */ @@ -597,7 +646,7 @@ protected: /*! * Vector of booleans indicating whether a phase is stable or not under * the current conditions. We use this to set the ROP's so that - * unphysical things don't happen + * unphysical things don't happen. * * length = number of phases in the object. By default all phases are stable. */ diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 9f44d6bc1..3a520983f 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -88,7 +88,7 @@ InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics& right) : /* * Call the assignment operator */ - *this = operator=(right); + operator=(right); } InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right) @@ -116,6 +116,7 @@ InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right) m_rxneqn = right.m_rxneqn; m_conc = right.m_conc; m_mu0 = right.m_mu0; + m_mu0_Kc = right.m_mu0_Kc; m_phi = right.m_phi; m_pot = right.m_pot; m_rwork = right.m_rwork; @@ -227,7 +228,7 @@ void InterfaceKinetics::getActivityConcentrations(doublereal* const conc) _update_rates_C(); copy(m_conc.begin(), m_conc.end(), conc); } - +//============================================================================================================================ void InterfaceKinetics::updateKc() { fill(m_rkcn.begin(), m_rkcn.end(), 0.0); @@ -235,30 +236,18 @@ void InterfaceKinetics::updateKc() if (m_nrev > 0) { /* * Get the vector of standard state electrochemical potentials for species in the Interfacial - * kinetics object and store it in m_mu0[] + * kinetics object and store it in m_mu0[] and m_mu0_Kc[] */ - size_t nsp, ik = 0; - doublereal rt = GasConstant*thermo(0).temperature(); - doublereal rrt = 1.0 / rt; - size_t np = nPhases(); - for (size_t n = 0; n < np; n++) { - thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]); - nsp = thermo(n).nSpecies(); - for (size_t k = 0; k < nsp; k++) { - m_mu0[ik] -= rt * thermo(n).logStandardConc(k); - m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k); - ik++; - } - } - + updateMu0(); + doublereal rrt = 1.0 / (GasConstant * thermo(0).temperature()); + // compute Delta mu^0 for all reversible reactions - m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(m_mu0), DATA_PTR(m_rkcn)); + m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(m_mu0_Kc), DATA_PTR(m_rkcn)); for (size_t i = 0; i < m_nrev; i++) { size_t irxn = m_revindex[i]; if (irxn == npos || irxn >= nReactions()) { - throw CanteraError("InterfaceKinetics", - "illegal value: irxn = "+int2str(irxn)); + throw CanteraError("InterfaceKinetics", "illegal value: irxn = "+int2str(irxn)); } // WARNING this may overflow HKM m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt); @@ -268,7 +257,27 @@ void InterfaceKinetics::updateKc() } } } - +//============================================================================================================================ +void InterfaceKinetics::updateMu0() +{ + /* + * 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[] + */ + size_t nsp, ik = 0; + doublereal rt = GasConstant * thermo(0).temperature(); + size_t np = nPhases(); + for (size_t n = 0; n < np; n++) { + thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]); + nsp = thermo(n).nSpecies(); + for (size_t k = 0; k < nsp; k++) { + m_mu0_Kc[ik] = m_mu0[ik] + Faraday * m_phi[n] * thermo(n).charge(k); + m_mu0_Kc[ik] -= rt * thermo(n).logStandardConc(k); + ik++; + } + } +} +//============================================================================================================================ void InterfaceKinetics::checkPartialEquil() { vector_fp dmu(nTotalSpecies(), 0.0); @@ -322,22 +331,12 @@ void InterfaceKinetics::getNetRatesOfProgress(doublereal* netROP) void InterfaceKinetics::getEquilibriumConstants(doublereal* kc) { - size_t ik=0; - doublereal rt = GasConstant*thermo(0).temperature(); - doublereal rrt = 1.0/rt; - for (size_t n = 0; n < nPhases(); n++) { - thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]); - size_t nsp = thermo(n).nSpecies(); - for (size_t k = 0; k < nsp; k++) { - m_mu0[ik] -= rt*thermo(n).logStandardConc(k); - m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k); - ik++; - } - } + updateMu0(); + doublereal rrt = 1.0 / (GasConstant * thermo(0).temperature()); - fill(kc, kc + m_ii, 0.0); + std::fill(kc, kc + m_ii, 0.0); - m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), kc); + m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0_Kc), kc); for (size_t i = 0; i < m_ii; i++) { kc[i] = exp(-kc[i]*rrt); @@ -601,13 +600,12 @@ void InterfaceKinetics::updateROP() void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG) { /* - * Get the chemical potentials of the species in the - * ideal gas solution. + * Get the chemical potentials of the species in the all of the phases used in the + * kinetics mechanism */ for (size_t n = 0; n < nPhases(); n++) { thermo(n).getChemPotentials(DATA_PTR(m_grt) + m_start[n]); } - /* * Use the stoichiometric manager to find deltaG for each * reaction. @@ -618,8 +616,7 @@ void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG) void InterfaceKinetics::getDeltaElectrochemPotentials(doublereal* deltaM) { /* - * Get the chemical potentials of the species in the - * ideal gas solution. + * Get the chemical potentials of the species */ size_t np = nPhases(); for (size_t n = 0; n < np; n++) { @@ -635,8 +632,7 @@ void InterfaceKinetics::getDeltaElectrochemPotentials(doublereal* deltaM) void InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH) { /* - * Get the partial molar enthalpy of all species in the - * ideal gas. + * Get the partial molar enthalpy of all species */ for (size_t n = 0; n < nPhases(); n++) { thermo(n).getPartialMolarEnthalpies(DATA_PTR(m_grt) + m_start[n]); @@ -664,7 +660,7 @@ void InterfaceKinetics::getDeltaEntropy(doublereal* deltaS) m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaS); } -void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaG) +void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaGSS) { /* * Get the standard state chemical potentials of the species. @@ -673,13 +669,13 @@ void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaG) * species at the temperature and pressure of the solution. */ for (size_t n = 0; n < nPhases(); n++) { - thermo(n).getStandardChemPotentials(DATA_PTR(m_grt) + m_start[n]); + thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]); } /* * Use the stoichiometric manager to find deltaG for each * reaction. */ - m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaG); + m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), deltaGSS); } void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH) @@ -941,6 +937,7 @@ void InterfaceKinetics::init() m_prxn.resize(m_kk); m_conc.resize(m_kk); m_mu0.resize(m_kk); + m_mu0_Kc.resize(m_kk); m_grt.resize(m_kk); m_pot.resize(m_kk, 0.0); m_phi.resize(nPhases(), 0.0);