research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Changes in InterfaceKinetics to allow for downstream hooks to add in the capability to model

Browse source 
arbitrarily fit open circuit potentials for intercalating electrodes.
This commit is contained in:
Harry Moffat 2014-06-05 22:43:41 +00:00
parent f2130b5e3f
commit 199cd4ceba
2 changed files with 104 additions and 58 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

75
include/cantera/kinetics/InterfaceKinetics.h
View file

@ -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.
*/

87
src/kinetics/InterfaceKinetics.cpp
View file

@ -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);