diff --git a/include/cantera/kinetics/InterfaceRate.h b/include/cantera/kinetics/InterfaceRate.h index 4700a4e30c8..ff106ce069a 100644 --- a/include/cantera/kinetics/InterfaceRate.h +++ b/include/cantera/kinetics/InterfaceRate.h @@ -57,6 +57,7 @@ struct InterfaceData : public BlowersMaselData electricPotentials.resize(nPhases, 0.); standardChemPotentials.resize(nSpecies, 0.); standardConcentrations.resize(nSpecies, 0.); + chemPotentials.resize(nSpecies, 0.); ready = true; } @@ -67,6 +68,7 @@ struct InterfaceData : public BlowersMaselData vector_fp electricPotentials; //!< electric potentials of phases vector_fp standardChemPotentials; //!< standard state chemical potentials vector_fp standardConcentrations; //!< standard state concentrations + vector_fp chemPotentials; }; @@ -123,6 +125,11 @@ class InterfaceRateBase return m_exchangeCurrentDensityFormulation; } + // String storing the input electrochemical kinetics form used if non-standard + std::string eChemForm() { + return m_eChemForm; + } + //! Build rate-specific parameters based on Reaction and Kinetics context //! @param rxn Reaction associated with rate parameterization //! @param kin Kinetics object associated with rate parameterization @@ -161,23 +168,34 @@ class InterfaceRateBase // Calculate reaction rate correction. Only modify those with a non-zero // activation energy. double correction = 1.; + double beta_tmp = m_beta; if (m_deltaPotential_RT != 0.) { // Comments preserved from previous implementation: // Below we decrease the activation energy below zero. // NOTE, there is some discussion about this point. Should we decrease the // activation energy below zero? I don't think this has been decided in any // definitive way. The treatment below is numerically more stable, however. - correction = exp(-m_beta * m_deltaPotential_RT); + if (m_eChemForm == "Marcus") { + // If Marcus kinetics is being used, adjust beta by the reorganization + // energy m_lambdaMarcus + beta_tmp += (m_etaF / 4. / m_lambdaMarcus); //m_etaF / 4 / m_lambdaMarcus; + } + correction = exp(-beta_tmp * m_deltaPotential_RT); } // Update correction if exchange current density formulation format is used. - if (m_exchangeCurrentDensityFormulation) { + if (m_eChemForm == "Butler-Volmer" || m_exchangeCurrentDensityFormulation) { // Comment preserved from previous implementation: // We need to have the straight chemical reaction rate constant to // come out of this calculation. double tmp = exp(-m_beta * m_deltaGibbs0_RT); tmp /= m_prodStandardConcentrations * Faraday; correction *= tmp; + } else if (m_eChemForm == "Marcus") { + // If Marcus kinetics is being used, adjust the correction value. + double tmp = exp(-m_lambda_RT / 4.) * exp(-beta_tmp * m_deltaGibbs_RT); + tmp /= m_prodStandardConcentrations * Faraday; + correction *= tmp; } return correction; } @@ -204,6 +222,31 @@ class InterfaceRateBase return NAN; } + //! Return the Marcus theory reorganization energy parameter + double lambdaMarcus() const { + if (m_chargeTransfer) { + return m_lambdaMarcus; + } + return NAN; + } + + //! Return the overpotential times Faraday's constant parameter + double etaF() { + if (m_chargeTransfer) { + return m_etaF; + } + return NAN; + } + + //! Return Marcus theory reorganization energy parameter + //! Divided by the Gas Constant and temperature + double lambda_RT() { + if (m_chargeTransfer) { + return m_lambda_RT; + } + return NAN; + } + //! Return site density [kmol/m^2] /*! * @warning This method is an experimental part of the %Cantera API and @@ -233,6 +276,11 @@ class InterfaceRateBase double m_mcov; //!< Coverage term in reaction rate bool m_chargeTransfer; //!< Boolean indicating use of electrochemistry bool m_exchangeCurrentDensityFormulation; //! Electrochemistry only + std::string m_eChemForm; //! String storing echem kinetics form + double m_lambdaMarcus; //! Marcus theory reorganization energy + double m_etaF; //! Overpotential times Faraday + double m_lambda_RT; //! Marcus theory reorganization energy over RT + double m_deltaGibbs_RT; //! Normalized Gibbs free energy change double m_beta; //!< Forward value of apparent electrochemical transfer coefficient double m_deltaPotential_RT; //!< Normalized electric potential energy change double m_deltaGibbs0_RT; //!< Normalized standard state Gibbs free energy change diff --git a/src/kinetics/InterfaceRate.cpp b/src/kinetics/InterfaceRate.cpp index cd830685dc6..aa01d74a5b9 100644 --- a/src/kinetics/InterfaceRate.cpp +++ b/src/kinetics/InterfaceRate.cpp @@ -77,6 +77,7 @@ bool InterfaceData::update(const ThermoPhase& phase, const Kinetics& kin) electricPotentials[n] = ph.electricPotential(); ph.getPartialMolarEnthalpies(partialMolarEnthalpies.data() + start); ph.getStandardChemPotentials(standardChemPotentials.data() + start); + ph.getChemPotentials(chemPotentials.data() + start); size_t nsp = ph.nSpecies(); for (size_t k = 0; k < nsp; k++) { // only used for exchange current density formulation @@ -102,6 +103,10 @@ InterfaceRateBase::InterfaceRateBase() , m_chargeTransfer(false) , m_exchangeCurrentDensityFormulation(false) , m_beta(0.5) + , m_lambdaMarcus(0.) + , m_etaF(NAN) + , m_lambda_RT(NAN) + , m_deltaGibbs_RT(NAN) , m_deltaPotential_RT(NAN) , m_deltaGibbs0_RT(NAN) , m_prodStandardConcentrations(NAN) @@ -117,6 +122,12 @@ void InterfaceRateBase::setParameters(const AnyMap& node) if (node.hasKey("beta")) { m_beta = node["beta"].asDouble(); } + if (node.hasKey("echem-kinetics-form")) { + m_eChemForm = node["echem-kinetics-form"].asString(); + if (m_eChemForm == "Marcus" && node.hasKey("lambda")) { + m_lambdaMarcus = node["lambda"].asDouble(); + } + } m_exchangeCurrentDensityFormulation = node.getBool( "exchange-current-density-formulation", false); } @@ -132,6 +143,12 @@ void InterfaceRateBase::getParameters(AnyMap& node) const if (m_beta != 0.5) { node["beta"] = m_beta; } + if (m_eChemForm != "") { + node["echem-kinetics-form"] = m_eChemForm; + if (m_eChemForm == "Marcus") { + node["lambda"] = m_lambdaMarcus; + } + } if (m_exchangeCurrentDensityFormulation) { node["exchange-current-density-formulation"] = true; } @@ -238,26 +255,45 @@ void InterfaceRateBase::updateFromStruct(const InterfaceData& shared_data) { // Update change in electrical potential energy if (m_chargeTransfer) { m_deltaPotential_RT = 0.; + m_etaF = 0.; for (const auto& ch : m_netCharges) { m_deltaPotential_RT += shared_data.electricPotentials[ch.first] * ch.second; } + m_etaF += m_deltaPotential_RT; m_deltaPotential_RT /= GasConstant * shared_data.temperature; } // Update quantities used for exchange current density formulation - if (m_exchangeCurrentDensityFormulation) { - m_deltaGibbs0_RT = 0.; + if (m_eChemForm != "") { m_prodStandardConcentrations = 1.; - for (const auto& item : m_stoichCoeffs) { - m_deltaGibbs0_RT += - shared_data.standardChemPotentials[item.first] * item.second; - if (item.second > 0.) { - m_prodStandardConcentrations *= - shared_data.standardConcentrations[item.first]; + if (m_eChemForm == "Butler-Volmer" || m_exchangeCurrentDensityFormulation) { + m_deltaGibbs0_RT = 0.; + for (const auto& item : m_stoichCoeffs) { + m_deltaGibbs0_RT += + shared_data.standardChemPotentials[item.first] * item.second; + if (item.second > 0.) { + m_prodStandardConcentrations *= + shared_data.standardConcentrations[item.first]; + } + } + m_deltaGibbs0_RT /= GasConstant * shared_data.temperature; + } + else if (m_eChemForm == "Marcus") { + m_deltaGibbs_RT = 0.; + m_lambda_RT = 0.; + for (const auto& item : m_stoichCoeffs) { + m_deltaGibbs_RT += + shared_data.chemPotentials[item.first] * item.second; + if (item.second > 0.) { + m_prodStandardConcentrations *= + shared_data.standardConcentrations[item.first]; + } } + m_etaF += m_deltaGibbs_RT; + m_deltaGibbs_RT /= GasConstant * shared_data.temperature; + m_lambda_RT = m_lambdaMarcus / GasConstant / shared_data.temperature; } - m_deltaGibbs0_RT /= GasConstant * shared_data.temperature; } }