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

Rewrite of InterfaceKinetics - next iteration.

Browse source
This commit is contained in:
Harry Moffat 2014-08-14 00:00:00 +00:00
parent 6741d8f7c6
commit 81556ea8aa
7 changed files with 1867 additions and 99 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
bin/make_doxygen
View file

@ -4,6 +4,7 @@
CUR=`pwd` CUR=`pwd`
CANTERA_SRC_ROOT=${CANTERA_SRC_ROOT:="$CUR"} CANTERA_SRC_ROOT=${CANTERA_SRC_ROOT:="$CUR"}
# /bin/rm -rf build/docs/doxygen
cd $CANTERA_SRC_ROOT cd $CANTERA_SRC_ROOT
# doxygen doc/doxygen/Doxyfile #doxygen doc/doxygen/Doxyfile
doxygen doc/doxygen/Doxyfile.tested doxygen doc/doxygen/Doxyfile.tested

1551
doc/doxygen/Doxyfile.tested Normal file
View file

File diff suppressed because it is too large Load diff

17
include/cantera/base/global.h
View file

@ -85,9 +85,19 @@ void popError();
* *
* There are two different types of input files within %Cantera: * There are two different types of input files within %Cantera:
* - ctml: This is an xml file laid out in such a way that %Cantera can * - ctml: This is an xml file laid out in such a way that %Cantera can
* interpret the contents. * interpret the contents. This is the essential input file within
* - cti: A human-readable ascii format for information that %Cantera * Cantera, and contains all elements that are involved with simulation,
* will read. * error propagation, data support, and versioning.
*
*
* - cti: A Chemkin-like input file that %Cantera will read. This file
* may not contain all of the information storred in the ctml file,
* nor all of options and equations of state that %Cantera can read.
* The file supports backwards compatibility with gas-phase mechanisms
* written for Chemkin.
*
* %Cantera takes its input from the ctml file. Given a file in cti format,
* %Cantera will perform a translation from the cti file into a ctml file.
* *
* %Cantera can take its input from both types of files. However, given a file * %Cantera can take its input from both types of files. However, given a file
* in cti format, the initial operation that %Cantera will perform is to * in cti format, the initial operation that %Cantera will perform is to
@ -99,6 +109,7 @@ void popError();
* *
* Other input routines in other modules: * Other input routines in other modules:
* @see importKinetics() * @see importKinetics()
*
* @{ * @{
*/ */

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

@ -273,12 +273,23 @@ public:
*/ */
virtual void updateMu0(); virtual void updateMu0();
//! Number of reactions in the mechanism
/*!
* @deprecated This is a duplicate of Kinetics::nReactions()
*/
size_t reactionNumber() const { size_t reactionNumber() const {
return m_ii; return m_ii;
} }
void addElementaryReaction(ReactionData& r); //! Add a single elementary reaction to the list of reactions for the object
//void addGlobalReaction(const ReactionData& r); /*!
* @param rdata
*/
void addElementaryReaction(ReactionData& rdata);
void addGlobalReaction(ReactionData& r);
void installReagents(const ReactionData& r); void installReagents(const ReactionData& r);
@ -301,23 +312,31 @@ public:
m_index[rxnNumber] = std::pair<int, size_t>(type, loc); m_index[rxnNumber] = std::pair<int, size_t>(type, loc);
} }
//! Apply corrections for interfacial charge transfer reactions //! Apply modifications for the fowward reaction rate for interfacial charge transfer reactions
/*! /*!
* For reactions that transfer charge across a potential difference, * For reactions that transfer charge across a potential difference,
* the activation energies are modified by the potential difference. * the activation energies are modified by the potential difference.
* (see, for example, ...). This method applies this correction. * (see, for example, ...). This method applies this correction.
* *
* @param kf Vector of forward reaction rate constants on which to have * @param kfwd Vector of forward reaction rate constants on which to have
* the correction applied * the voltage correction applied
*/ */
void applyButlerVolmerCorrection(doublereal* const kf); void applyVoltageKfwdCorrection(doublereal* const kfwd);
//! When an electrode reaction rate is optionally specified in terms of its //! When an electrode reaction rate is optionally specified in terms of its
//! exchange current density, adjust kfwd to the standard reaction rate constant form and units. //! exchange current density, adjust kfwd to the standard reaction rate constant form and units.
//! When the BV reaction types are used, keep the exchange current density form.
/*! /*!
* For a reaction rate constant that was given in units of Amps/m2 (exchange current * For a reaction rate constant that was given in units of Amps/m2 (exchange current
* density formulation with iECDFormulation == true), convert the rate to * density formulation with iECDFormulation == true), convert the rate to
* kmoles/m2/s. * kmoles/m2/s.
*
* For a reaction rate constant that was given in units of kmol/m2/sec when the
* reaction type is a butler-volmer form, convert it to exchange current density
* form (amps/m2).
*
* @param kfwd Vector of forward reaction rate constants, given in either
* normal form or in exchange current density form.
*/ */
void convertExchangeCurrentDensityFormulation(doublereal* const kfwd); void convertExchangeCurrentDensityFormulation(doublereal* const kfwd);
@ -444,6 +463,13 @@ protected:
*/ */
mutable std::vector<std::map<size_t, doublereal> > m_prxn; mutable std::vector<std::map<size_t, doublereal> > m_prxn;
//! Vector of reactionType for the reactions defined within this object
/*!
* Length = number of reactions, m_ii
* contains the type of reaction.
*/
vector_int reactionType_;
//! String expression for each rxn //! String expression for each rxn
/*! /*!
* Vector of strings of length m_ii, the number of * Vector of strings of length m_ii, the number of
@ -452,7 +478,7 @@ protected:
*/ */
std::vector<std::string> m_rxneqn; std::vector<std::string> m_rxneqn;
//! an array of generalized concentrations for each species //! Array of concentrations for each species in the kinetics mechanism
/*! /*!
* An array of generalized concentrations \f$ C_k \f$ that are defined * An array of generalized concentrations \f$ C_k \f$ that are defined
* such that \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$ is a standard * such that \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$ is a standard
@ -466,6 +492,20 @@ protected:
*/ */
vector_fp m_conc; vector_fp m_conc;
//! Array of activity concentrations for each species in the kinetics object
/*!
* An array of activity concentrations \f$ Ca_k \f$ that are defined
* such that \f$ a_k = Ca_k / C^0_k, \f$ where \f$ C^0_k \f$ is a standard
* concentration. These activity concentrations are used by this
* kinetics manager class to compute the forward and reverse rates of
* elementary reactions. The "units" for the concentrations of each phase
* depend upon the implementation of kinetics within that phase. The order
* of the species within the vector is based on the order of listed
* ThermoPhase objects in the class, and the order of the species within
* each ThermoPhase class.
*/
vector_fp m_actConc;
//! Vector of standard state chemical potentials for all species //! Vector of standard state chemical potentials for all species
/*! /*!
* This vector contains a temporary vector of standard state chemical * This vector contains a temporary vector of standard state chemical
@ -507,12 +547,14 @@ protected:
*/ */
vector_fp m_pot; vector_fp m_pot;
//! Vector temporary //! Storage for the net electric energy change due to reaction.
/*! /*!
* Length is number of reactions. It's used to store the * Length is number of reactions. It's used to store the
* voltage contribution to the activation energy. * net electric potential energy change due to the reaction.
*
* deltaElectricEnergy_[jrxn] = sum_i ( F V_i z_i nu_ij)
*/ */
vector_fp m_rwork; vector_fp deltaElectricEnergy_;
//! Vector of raw activation energies for the reactions //! Vector of raw activation energies for the reactions
/*! /*!
@ -549,7 +591,7 @@ protected:
//! reactions in the mechanism //! reactions in the mechanism
/*! /*!
* Vector of reaction indices which involve current transfers. This provides * Vector of reaction indices which involve current transfers. This provides
* an index into the m_beta, ctrxn_BVform array. * an index into the m_beta and m_ctrxn_BVform array.
* *
* irxn = m_ctrxn[i] * irxn = m_ctrxn[i]
*/ */
@ -573,6 +615,11 @@ protected:
//! is described by an exchange current density rate constant expression //! is described by an exchange current density rate constant expression
/*! /*!
* Length is equal to the number of reactions with charge transfer coefficients, m_ctrxn[] * Length is equal to the number of reactions with charge transfer coefficients, m_ctrxn[]
*
* m_ctrxn_ecdf[irxn] = 0 This means that the rate coefficient calculator will calculate
* the rate constant as a chemical forward rate constant, a standard format.
* m_ctrxn_ecdf[irxn] = 1 this means that the rate coefficient calculator will calculate
* the rate constant as an exchange current density rate constant expression.
*/ */
vector_int m_ctrxn_ecdf; vector_int m_ctrxn_ecdf;

10
include/cantera/kinetics/reaction_defs.h
View file

@ -19,14 +19,18 @@ const int NONE = 0;
//@{ //@{
/** //! A reaction with a rate coefficient that depends only on temperature and voltage
* A reaction with a rate coefficient that depends only on //! that also obeys mass-action kinetics.
/*!
* Here mass-action kinetics is defined as the reaction orders being equal to
* the reaction's stoichiometry.
*
* temperature. Example: O + OH <-> O2 + H * temperature. Example: O + OH <-> O2 + H
*/ */
const int ELEMENTARY_RXN = 1; const int ELEMENTARY_RXN = 1;
/** /**
* A reaction that requires a third-body collision partner. Example: * A gas-phase reaction that requires a third-body collision partner. Example:
* O2 + M <-> O + O + M * O2 + M <-> O + O + M
*/ */
const int THREE_BODY_RXN = 2; const int THREE_BODY_RXN = 2;

297
src/kinetics/InterfaceKinetics.cpp
View file

@ -17,6 +17,7 @@ using namespace std;
namespace Cantera namespace Cantera
{ {
//============================================================================================================================
InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) : InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) :
Kinetics(), Kinetics(),
m_redo_rates(false), m_redo_rates(false),
@ -26,6 +27,7 @@ InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) :
m_integrator(0), m_integrator(0),
m_beta(0), m_beta(0),
m_ctrxn(0), m_ctrxn(0),
m_ctrxn_BVform(0),
m_ctrxn_ecdf(0), m_ctrxn_ecdf(0),
m_StandardConc(0), m_StandardConc(0),
m_deltaG0(0), m_deltaG0(0),
@ -50,12 +52,12 @@ InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) :
addPhase(*thermo); addPhase(*thermo);
} }
} }
//============================================================================================================================
InterfaceKinetics::~InterfaceKinetics() InterfaceKinetics::~InterfaceKinetics()
{ {
delete m_integrator; delete m_integrator;
} }
//============================================================================================================================
InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics& right) : InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics& right) :
Kinetics(), Kinetics(),
m_redo_rates(false), m_redo_rates(false),
@ -65,6 +67,7 @@ InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics& right) :
m_integrator(0), m_integrator(0),
m_beta(0), m_beta(0),
m_ctrxn(0), m_ctrxn(0),
m_ctrxn_BVform(0),
m_ctrxn_ecdf(0), m_ctrxn_ecdf(0),
m_StandardConc(0), m_StandardConc(0),
m_deltaG0(0), m_deltaG0(0),
@ -90,7 +93,7 @@ InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics& right) :
*/ */
operator=(right); operator=(right);
} }
//============================================================================================================================
InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right) InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right)
{ {
/* /*
@ -113,18 +116,21 @@ InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right)
m_nrev = right.m_nrev; m_nrev = right.m_nrev;
m_rrxn = right.m_rrxn; m_rrxn = right.m_rrxn;
m_prxn = right.m_prxn; m_prxn = right.m_prxn;
reactionType_ = right.reactionType_;
m_rxneqn = right.m_rxneqn; m_rxneqn = right.m_rxneqn;
m_conc = right.m_conc; m_conc = right.m_conc;
m_actConc = right.m_actConc;
m_mu0 = right.m_mu0; m_mu0 = right.m_mu0;
m_mu0_Kc = right.m_mu0_Kc; m_mu0_Kc = right.m_mu0_Kc;
m_phi = right.m_phi; m_phi = right.m_phi;
m_pot = right.m_pot; m_pot = right.m_pot;
m_rwork = right.m_rwork; deltaElectricEnergy_ = right.deltaElectricEnergy_;
m_E = right.m_E; m_E = right.m_E;
m_surf = right.m_surf; //DANGER - shallow copy m_surf = right.m_surf; //DANGER - shallow copy
m_integrator = right.m_integrator; //DANGER - shallow copy m_integrator = right.m_integrator; //DANGER - shallow copy
m_beta = right.m_beta; m_beta = right.m_beta;
m_ctrxn = right.m_ctrxn; m_ctrxn = right.m_ctrxn;
m_ctrxn_BVform = right.m_ctrxn_BVform;
m_ctrxn_ecdf = right.m_ctrxn_ecdf; m_ctrxn_ecdf = right.m_ctrxn_ecdf;
m_StandardConc = right.m_StandardConc; m_StandardConc = right.m_StandardConc;
m_deltaG0 = right.m_deltaG0; m_deltaG0 = right.m_deltaG0;
@ -152,12 +158,12 @@ InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right)
return *this; return *this;
} }
//============================================================================================================================
int InterfaceKinetics::type() const int InterfaceKinetics::type() const
{ {
return cInterfaceKinetics; return cInterfaceKinetics;
} }
//============================================================================================================================
Kinetics* InterfaceKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const Kinetics* InterfaceKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const
{ {
InterfaceKinetics* iK = new InterfaceKinetics(*this); InterfaceKinetics* iK = new InterfaceKinetics(*this);
@ -170,25 +176,36 @@ void InterfaceKinetics::setElectricPotential(int n, doublereal V)
thermo(n).setElectricPotential(V); thermo(n).setElectricPotential(V);
m_redo_rates = true; m_redo_rates = true;
} }
//============================================================================================================================
void InterfaceKinetics::_update_rates_T() void InterfaceKinetics::_update_rates_T()
{ {
_update_rates_phi(); _update_rates_phi();
if (m_has_coverage_dependence) { if (m_has_coverage_dependence) {
m_surf->getCoverages(DATA_PTR(m_conc)); m_surf->getCoverages(DATA_PTR(m_actConc));
m_rates.update_C(DATA_PTR(m_conc)); m_rates.update_C(DATA_PTR(m_actConc));
m_redo_rates = true; m_redo_rates = true;
} }
//
// Go find the temperature from the surface
//
doublereal T = thermo(surfacePhaseIndex()).temperature(); doublereal T = thermo(surfacePhaseIndex()).temperature();
m_redo_rates = true; m_redo_rates = true;
if (T != m_temp || m_redo_rates) { if (T != m_temp || m_redo_rates) {
m_logtemp = log(T); m_logtemp = log(T);
//
// Calculate the forward rate constant by calling m_rates and store it in m_rfn[]
//
m_rates.update(T, m_logtemp, DATA_PTR(m_rfn)); m_rates.update(T, m_logtemp, DATA_PTR(m_rfn));
//
// If we need to do conversions between exchange current density formulation and regular formulation
// (either way) do it here.
//
if (m_has_exchange_current_density_formulation) { if (m_has_exchange_current_density_formulation) {
convertExchangeCurrentDensityFormulation(DATA_PTR(m_rfn)); convertExchangeCurrentDensityFormulation(DATA_PTR(m_rfn));
} }
if (m_has_electrochem_rxns) { if (m_has_electrochem_rxns) {
applyButlerVolmerCorrection(DATA_PTR(m_rfn)); applyVoltageKfwdCorrection(DATA_PTR(m_rfn));
} }
m_temp = T; m_temp = T;
updateKc(); updateKc();
@ -196,7 +213,7 @@ void InterfaceKinetics::_update_rates_T()
m_redo_rates = false; m_redo_rates = false;
} }
} }
//============================================================================================================================
void InterfaceKinetics::_update_rates_phi() void InterfaceKinetics::_update_rates_phi()
{ {
// //
@ -209,10 +226,11 @@ void InterfaceKinetics::_update_rates_phi()
} }
} }
} }
//============================================================================================================================
void InterfaceKinetics::_update_rates_C() void InterfaceKinetics::_update_rates_C()
{ {
for (size_t n = 0; n < nPhases(); n++) { for (size_t n = 0; n < nPhases(); n++) {
const ThermoPhase* tp = m_thermo[n];
/* /*
* We call the getActivityConcentrations function of each * We call the getActivityConcentrations function of each
* ThermoPhase class that makes up this kinetics object to * ThermoPhase class that makes up this kinetics object to
@ -221,7 +239,11 @@ void InterfaceKinetics::_update_rates_C()
* are integer indices for that vector denoting the start of the * are integer indices for that vector denoting the start of the
* species for each phase. * species for each phase.
*/ */
thermo(n).getActivityConcentrations(DATA_PTR(m_conc) + m_start[n]); tp->getActivityConcentrations(DATA_PTR(m_actConc) + m_start[n]);
//
// Get regular concentrations too
//
tp->getConcentrations(DATA_PTR(m_conc) + m_start[n]);
} }
m_ROP_ok = false; m_ROP_ok = false;
} }
@ -229,7 +251,7 @@ void InterfaceKinetics::_update_rates_C()
void InterfaceKinetics::getActivityConcentrations(doublereal* const conc) void InterfaceKinetics::getActivityConcentrations(doublereal* const conc)
{ {
_update_rates_C(); _update_rates_C();
copy(m_conc.begin(), m_conc.end(), conc); copy(m_actConc.begin(), m_actConc.end(), conc);
} }
//============================================================================================================================ //============================================================================================================================
void InterfaceKinetics::updateKc() void InterfaceKinetics::updateKc()
@ -394,10 +416,12 @@ void InterfaceKinetics::getNetProductionRates(doublereal* net)
updateROP(); updateROP();
m_rxnstoich.getNetProductionRates(m_kk, &m_ropnet[0], net); m_rxnstoich.getNetProductionRates(m_kk, &m_ropnet[0], net);
} }
//===========================================================================================================
void InterfaceKinetics::applyButlerVolmerCorrection(doublereal* const kf) void InterfaceKinetics::applyVoltageKfwdCorrection(doublereal* const kf)
{ {
// compute the electrical potential energy of each species //
// Compute the electrical potential energy of each species
//
size_t ik = 0; size_t ik = 0;
for (size_t n = 0; n < nPhases(); n++) { for (size_t n = 0; n < nPhases(); n++) {
size_t nsp = thermo(n).nSpecies(); size_t nsp = thermo(n).nSpecies();
@ -406,11 +430,12 @@ void InterfaceKinetics::applyButlerVolmerCorrection(doublereal* const kf)
ik++; ik++;
} }
} }
//
// Compute the change in electrical potential energy for each // Compute the change in electrical potential energy for each
// reaction. This will only be non-zero if a potential // reaction. This will only be non-zero if a potential
// difference is present. // difference is present.
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_pot), DATA_PTR(m_rwork)); //
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_pot), DATA_PTR(deltaElectricEnergy_));
// Modify the reaction rates. Only modify those with a // Modify the reaction rates. Only modify those with a
// non-zero activation energy. Below we decrease the // non-zero activation energy. Below we decrease the
@ -428,13 +453,13 @@ void InterfaceKinetics::applyButlerVolmerCorrection(doublereal* const kf)
#endif #endif
for (size_t i = 0; i < m_beta.size(); i++) { for (size_t i = 0; i < m_beta.size(); i++) {
size_t irxn = m_ctrxn[i]; size_t irxn = m_ctrxn[i];
eamod = m_beta[i]*m_rwork[irxn]; eamod = m_beta[i] * deltaElectricEnergy_[irxn];
if (eamod != 0.0) { if (eamod != 0.0) {
#ifdef DEBUG_KIN_MODE #ifdef DEBUG_KIN_MODE
ea = GasConstant * m_E[irxn]; ea = GasConstant * m_E[irxn];
if (eamod + ea < 0.0) { if (eamod + ea < 0.0) {
writelog("Warning: act energy mod too large!\n"); writelog("Warning: act energy mod too large!\n");
writelog(" Delta phi = "+fp2str(m_rwork[irxn]/Faraday)+"\n"); writelog(" Delta phi = "+fp2str(deltaElectricEnergy_[irxn]/Faraday)+"\n");
writelog(" Delta Ea = "+fp2str(eamod)+"\n"); writelog(" Delta Ea = "+fp2str(eamod)+"\n");
writelog(" Ea = "+fp2str(ea)+"\n"); writelog(" Ea = "+fp2str(ea)+"\n");
for (n = 0; n < np; n++) { for (n = 0; n < np; n++) {
@ -451,25 +476,61 @@ void InterfaceKinetics::applyButlerVolmerCorrection(doublereal* const kf)
} }
//================================================================================================================== //==================================================================================================================
/* /*
* For a reaction rate that was given in units of Amps/m2 (exchange current * For a reaction rate constant that was given in units of Amps/m2 (exchange current
* density formulation with iECDFormulation == true), convert the rate to * density formulation with iECDFormulation == true), convert the rate to
* kmoles/m2/s. * kmoles/m2/s.
* RENAMED THIS METHOD from "apply" to "convert" *
* For a reaction rate constant that was given in units of kmol/m2/sec when the
* reaction type is a butler-volmer form, convert it to exchange current density
* form (amps/m2).
*
*/ */
void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* const kfwd) void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* const kfwd)
{ {
updateExchangeCurrentQuantities(); updateExchangeCurrentQuantities();
doublereal rt = GasConstant*thermo(0).temperature(); doublereal rt = GasConstant * thermo(0).temperature();
doublereal rrt = 1.0/rt; doublereal rrt = 1.0/rt;
//
// Loop over all reactions which are defined to have a voltage transfer coefficient that
// affects the activity energy for the reaction
//
for (size_t i = 0; i < m_ctrxn.size(); i++) { for (size_t i = 0; i < m_ctrxn.size(); i++) {
size_t irxn = m_ctrxn[i]; size_t irxn = m_ctrxn[i];
int iECDFormulation = m_ctrxn_ecdf[i]; //
if (iECDFormulation) { // Determine whether the reaction rate constant is in an exchange current density formulation format.
double tmp = exp(- m_beta[i] * m_deltaG0[irxn] * rrt); //
double tmp2 = m_ProdStanConcReac[irxn]; int iECDFormulation = m_ctrxn_ecdf[i];
tmp *= 1.0 / tmp2 / Faraday; if (iECDFormulation) {
kfwd[irxn] *= tmp; //
} // 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 (m_ctrxn_BVform[i] == 0) {
//
// Calculate the term and modify the forward reaction
//
double tmp = exp(- m_beta[i] * m_deltaG0[irxn] * rrt);
double tmp2 = m_ProdStanConcReac[irxn];
tmp *= 1.0 / tmp2 / Faraday;
kfwd[irxn] *= tmp;
}
} else {
//
// 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
// substitute it.
//
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
//
double tmp = exp(m_beta[i] * m_deltaG0[irxn] * rrt);
double tmp2 = m_ProdStanConcReac[irxn];
tmp *= Faraday * tmp2;
kfwd[irxn] *= tmp;
}
}
} }
} }
//================================================================================================================== //==================================================================================================================
@ -498,10 +559,10 @@ void InterfaceKinetics::getRevRateConstants(doublereal* krev, bool doIrreversibl
multiply_each(krev, krev + nReactions(), m_rkcn.begin()); multiply_each(krev, krev + nReactions(), m_rkcn.begin());
} }
} }
//============================================================================================================================
void InterfaceKinetics::updateROP() void InterfaceKinetics::updateROP()
{ {
// evaluate rate and equilibrium constants at temperature and phi (electric potential) // evaluate rate constants and equilibrium constants at temperature and phi (electric potential)
_update_rates_T(); _update_rates_T();
// get updated activities (rates updated below) // get updated activities (rates updated below)
_update_rates_C(); _update_rates_C();
@ -527,11 +588,11 @@ void InterfaceKinetics::updateROP()
// multiply ropf by the actyivity concentration reaction orders to obtain // multiply ropf by the actyivity concentration reaction orders to obtain
// the forward rates of progress. // the forward rates of progress.
// //
m_rxnstoich.multiplyReactants(DATA_PTR(m_conc), DATA_PTR(m_ropf)); m_rxnstoich.multiplyReactants(DATA_PTR(m_actConc), DATA_PTR(m_ropf));
// for reversible reactions, multiply ropr by concentration // for reversible reactions, multiply ropr by the activity concentration
// products // products
m_rxnstoich.multiplyRevProducts(DATA_PTR(m_conc), m_rxnstoich.multiplyRevProducts(DATA_PTR(m_actConc),
DATA_PTR(m_ropr)); DATA_PTR(m_ropr));
for (size_t j = 0; j != m_ii; ++j) { for (size_t j = 0; j != m_ii; ++j) {
@ -601,7 +662,7 @@ void InterfaceKinetics::updateROP()
m_ROP_ok = true; m_ROP_ok = true;
} }
//==================================================================================================================
void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG) void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG)
{ {
/* /*
@ -617,7 +678,7 @@ void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG)
*/ */
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaG); m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaG);
} }
//==================================================================================================================
void InterfaceKinetics::getDeltaElectrochemPotentials(doublereal* deltaM) void InterfaceKinetics::getDeltaElectrochemPotentials(doublereal* deltaM)
{ {
/* /*
@ -633,7 +694,7 @@ void InterfaceKinetics::getDeltaElectrochemPotentials(doublereal* deltaM)
*/ */
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaM); m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaM);
} }
//==================================================================================================================
void InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH) void InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH)
{ {
/* /*
@ -728,11 +789,25 @@ void InterfaceKinetics::getDeltaSSEntropy(doublereal* deltaS)
//============================================================================================================================ //============================================================================================================================
void InterfaceKinetics::addReaction(ReactionData& r) void InterfaceKinetics::addReaction(ReactionData& r)
{ {
/* int reactionType = r.reactionType;
* Install the rate coefficient for the current reaction
* in the appropriate data structure. reactionType_.push_back(reactionType);
*/
addElementaryReaction(r); if ((reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN ) ||
(reactionType == BUTLERVOLMER_RXN ) ||
(reactionType == SURFACEAFFINITY_RXN) ||
(reactionType == GLOBAL_RXN)) {
//
// Add global reactions
//
addGlobalReaction(r);
} else {
/*
* Install the rate coefficient for the current reaction
* in the appropriate data structure.
*/
addElementaryReaction(r);
}
/* /*
* Add the reactants and products for m_ropnet;the current reaction * Add the reactants and products for m_ropnet;the current reaction
* to the various stoichiometric coefficient arrays. * to the various stoichiometric coefficient arrays.
@ -772,47 +847,48 @@ void InterfaceKinetics::addReaction(ReactionData& r)
} }
//============================================================================================================================ //============================================================================================================================
void InterfaceKinetics::addElementaryReaction(ReactionData& r) void InterfaceKinetics::addElementaryReaction(ReactionData& rdata)
{ {
// install rate coeff calculator // install rate coeff calculator
vector_fp& rp = r.rateCoeffParameters; vector_fp& rp = rdata.rateCoeffParameters;
size_t ncov = r.cov.size(); size_t ncov = rdata.cov.size();
if (ncov > 3) { if (ncov > 3) {
m_has_coverage_dependence = true; m_has_coverage_dependence = true;
} }
for (size_t m = 0; m < ncov; m++) { for (size_t m = 0; m < ncov; m++) {
rp.push_back(r.cov[m]); rp.push_back(rdata.cov[m]);
} }
/* /*
* Temporarily change the reaction rate coefficient type to surface arrhenius. * Temporarily change the reaction rate coefficient type to surface arrhenius.
* This is what is expected. We'll handle exchange current types below by hand. * This is what is expected. We'll handle exchange current types below by hand.
*/ */
int reactionRateCoeffType_orig = r.rateCoeffType; int reactionRateCoeffType_orig = rdata.rateCoeffType;
if (r.rateCoeffType == EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE) { if (rdata.rateCoeffType == EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE) {
r.rateCoeffType = SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE; rdata.rateCoeffType = SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE;
} }
if (r.rateCoeffType == ARRHENIUS_REACTION_RATECOEFF_TYPE) { if (rdata.rateCoeffType == ARRHENIUS_REACTION_RATECOEFF_TYPE) {
r.rateCoeffType = SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE; rdata.rateCoeffType = SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE;
} }
/* /*
* Install the reaction rate into the vector of reactions handled by this class * Install the reaction rate into the vector of reactions handled by this class
*/ */
size_t iloc = m_rates.install(m_ii, r); size_t iloc = m_rates.install(m_ii, rdata);
/* /*
* Change the reaction rate coefficient type back to its original value * Change the reaction rate coefficient type back to its original value
*/ */
r.rateCoeffType = reactionRateCoeffType_orig; rdata.rateCoeffType = reactionRateCoeffType_orig;
// store activation energy // store activation energy
m_E.push_back(r.rateCoeffParameters[2]); m_E.push_back(rdata.rateCoeffParameters[2]);
if (r.beta > 0.0) { if (rdata.beta > 0.0) {
m_has_electrochem_rxns = true; m_has_electrochem_rxns = true;
m_beta.push_back(r.beta); m_beta.push_back(rdata.beta);
m_ctrxn.push_back(reactionNumber()); m_ctrxn.push_back(m_ii);
if (r.rateCoeffType == EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE) { m_ctrxn_BVform.push_back(0);
if (rdata.rateCoeffType == EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE) {
m_has_exchange_current_density_formulation = true; m_has_exchange_current_density_formulation = true;
m_ctrxn_ecdf.push_back(1); m_ctrxn_ecdf.push_back(1);
} else { } else {
@ -821,10 +897,72 @@ void InterfaceKinetics::addElementaryReaction(ReactionData& r)
} }
// add constant term to rate coeff value vector // add constant term to rate coeff value vector
m_rfn.push_back(r.rateCoeffParameters[0]); m_rfn.push_back(rdata.rateCoeffParameters[0]);
registerReaction(reactionNumber(), ELEMENTARY_RXN, iloc); registerReaction(reactionNumber(), ELEMENTARY_RXN, iloc);
} }
//============================================================================================================================ //============================================================================================================================
void InterfaceKinetics::addGlobalReaction(ReactionData& rdata)
{
//
// Install rate coeff calculator
// This is done no matter what the type of reaction it is
//
vector_fp& rp = rdata.rateCoeffParameters;
size_t ncov = rdata.cov.size();
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
//
int reactionType = rdata.reactionType;
/*
* Temporarily change the reaction rate coefficient type to surface arrhenius.
* This is what is expected. We'll handle exchange current types below by hand.
*/
int reactionRateCoeffType_orig = rdata.rateCoeffType;
if (rdata.rateCoeffType == EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE) {
rdata.rateCoeffType = SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE;
}
if (rdata.rateCoeffType == ARRHENIUS_REACTION_RATECOEFF_TYPE) {
rdata.rateCoeffType = SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE;
}
/*
* Install the reaction rate into the vector of reactions handled by this class
*/
size_t iloc = m_rates.install(m_ii, rdata);
/*
* Change the reaction rate coefficient type back to its original value
*/
rdata.rateCoeffType = reactionRateCoeffType_orig;
// store activation energy
m_E.push_back(rdata.rateCoeffParameters[2]);
if (rdata.beta > 0.0) {
m_has_electrochem_rxns = true;
m_beta.push_back(rdata.beta);
m_ctrxn.push_back(m_ii);
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);
}
}
// add constant term to rate coeff value vector
m_rfn.push_back(rdata.rateCoeffParameters[0]);
registerReaction(m_ii, ELEMENTARY_RXN, iloc);
}
//==================================================================================================================
void InterfaceKinetics::setIOFlag(int ioFlag) void InterfaceKinetics::setIOFlag(int ioFlag)
{ {
m_ioFlag = ioFlag; m_ioFlag = ioFlag;
@ -832,7 +970,7 @@ void InterfaceKinetics::setIOFlag(int ioFlag)
m_integrator->setIOFlag(ioFlag); m_integrator->setIOFlag(ioFlag);
} }
} }
//==================================================================================================================
void InterfaceKinetics::installReagents(const ReactionData& r) void InterfaceKinetics::installReagents(const ReactionData& r)
{ {
@ -926,14 +1064,14 @@ void InterfaceKinetics::installReagents(const ReactionData& r)
m_nirrev++; m_nirrev++;
} }
} }
//==================================================================================================================
void InterfaceKinetics::addPhase(thermo_t& thermo) void InterfaceKinetics::addPhase(thermo_t& thermo)
{ {
Kinetics::addPhase(thermo); Kinetics::addPhase(thermo);
m_phaseExists.push_back(true); m_phaseExists.push_back(true);
m_phaseIsStable.push_back(true); m_phaseIsStable.push_back(true);
} }
//==================================================================================================================
void InterfaceKinetics::init() void InterfaceKinetics::init()
{ {
m_kk = 0; m_kk = 0;
@ -942,6 +1080,7 @@ void InterfaceKinetics::init()
} }
m_rrxn.resize(m_kk); m_rrxn.resize(m_kk);
m_prxn.resize(m_kk); m_prxn.resize(m_kk);
m_actConc.resize(m_kk);
m_conc.resize(m_kk); m_conc.resize(m_kk);
m_mu0.resize(m_kk); m_mu0.resize(m_kk);
m_mu0_Kc.resize(m_kk); m_mu0_Kc.resize(m_kk);
@ -949,12 +1088,12 @@ void InterfaceKinetics::init()
m_pot.resize(m_kk, 0.0); m_pot.resize(m_kk, 0.0);
m_phi.resize(nPhases(), 0.0); m_phi.resize(nPhases(), 0.0);
} }
//==================================================================================================================
void InterfaceKinetics::finalize() void InterfaceKinetics::finalize()
{ {
Kinetics::finalize(); Kinetics::finalize();
size_t safe_reaction_size = std::max<size_t>(nReactions(), 1); size_t safe_reaction_size = std::max<size_t>(m_ii, 1);
m_rwork.resize(safe_reaction_size); deltaElectricEnergy_.resize(safe_reaction_size);
size_t ks = reactionPhaseIndex(); size_t ks = reactionPhaseIndex();
if (ks == npos) throw CanteraError("InterfaceKinetics::finalize", if (ks == npos) throw CanteraError("InterfaceKinetics::finalize",
"no surface phase is present."); "no surface phase is present.");
@ -984,7 +1123,7 @@ void InterfaceKinetics::finalize()
m_finalized = true; m_finalized = true;
} }
//==================================================================================================================
doublereal InterfaceKinetics::electrochem_beta(size_t irxn) const doublereal InterfaceKinetics::electrochem_beta(size_t irxn) const
{ {
for (size_t i = 0; i < m_ctrxn.size(); i++) { for (size_t i = 0; i < m_ctrxn.size(); i++) {
@ -994,12 +1133,12 @@ doublereal InterfaceKinetics::electrochem_beta(size_t irxn) const
} }
return 0.0; return 0.0;
} }
//==================================================================================================================
bool InterfaceKinetics::ready() const bool InterfaceKinetics::ready() const
{ {
return m_finalized; return m_finalized;
} }
//==================================================================================================================
void InterfaceKinetics::advanceCoverages(doublereal tstep) void InterfaceKinetics::advanceCoverages(doublereal tstep)
{ {
if (m_integrator == 0) { if (m_integrator == 0) {
@ -1012,7 +1151,7 @@ void InterfaceKinetics::advanceCoverages(doublereal tstep)
delete m_integrator; delete m_integrator;
m_integrator = 0; m_integrator = 0;
} }
//==================================================================================================================
void InterfaceKinetics::solvePseudoSteadyStateProblem( void InterfaceKinetics::solvePseudoSteadyStateProblem(
int ifuncOverride, doublereal timeScaleOverride) int ifuncOverride, doublereal timeScaleOverride)
{ {
@ -1029,7 +1168,7 @@ void InterfaceKinetics::solvePseudoSteadyStateProblem(
*/ */
m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride); m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride);
} }
//==================================================================================================================
void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists) void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists)
{ {
if (iphase >= m_thermo.size()) { if (iphase >= m_thermo.size()) {
@ -1051,7 +1190,7 @@ void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists)
} }
} }
//==================================================================================================================
int InterfaceKinetics::phaseExistence(const size_t iphase) const int InterfaceKinetics::phaseExistence(const size_t iphase) const
{ {
if (iphase >= m_thermo.size()) { if (iphase >= m_thermo.size()) {
@ -1059,7 +1198,7 @@ int InterfaceKinetics::phaseExistence(const size_t iphase) const
} }
return m_phaseExists[iphase]; return m_phaseExists[iphase];
} }
//==================================================================================================================
int InterfaceKinetics::phaseStability(const size_t iphase) const int InterfaceKinetics::phaseStability(const size_t iphase) const
{ {
if (iphase >= m_thermo.size()) { if (iphase >= m_thermo.size()) {
@ -1067,7 +1206,7 @@ int InterfaceKinetics::phaseStability(const size_t iphase) const
} }
return m_phaseIsStable[iphase]; return m_phaseIsStable[iphase];
} }
//==================================================================================================================
void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStable) void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStable)
{ {
if (iphase >= m_thermo.size()) { if (iphase >= m_thermo.size()) {
@ -1079,10 +1218,10 @@ void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStabl
m_phaseIsStable[iphase] = false; m_phaseIsStable[iphase] = false;
} }
} }
//==================================================================================================================
void EdgeKinetics::finalize() void EdgeKinetics::finalize()
{ {
m_rwork.resize(std::max<size_t>(nReactions(), 1)); deltaElectricEnergy_.resize(std::max<size_t>(m_ii, 1));
size_t ks = reactionPhaseIndex(); size_t ks = reactionPhaseIndex();
if (ks == npos) throw CanteraError("EdgeKinetics::finalize", if (ks == npos) throw CanteraError("EdgeKinetics::finalize",
"no edge phase is present."); "no edge phase is present.");
@ -1104,5 +1243,5 @@ void EdgeKinetics::finalize()
m_finalized = true; m_finalized = true;
} }
//==================================================================================================================
} }

19
src/kinetics/importKinetics.cpp
View file

@ -485,7 +485,22 @@ static void getStick(const XML_Node& node, Kinetics& kin,
E = getFloat(node, "E", "actEnergy"); E = getFloat(node, "E", "actEnergy");
E /= GasConstant; E /= GasConstant;
} }
//=====================================================================================================
//! Read the XML data concerning the coverage dependence of an interfacial reaction
/*!
* @param node XML node with name reaction containing the reaction information
* @param surfphase Surface phase
* @param rdata Reaction data for the reaction.
*
* Example:
* @verbatim
<coverage species="CH3*">
<a> 1.0E-5 </a>
<m> 0.0 </m>
<actEnergy> 0.0 </actEnergy>
</coverage>
@endverbatim
*/
static void getCoverageDependence(const XML_Node& node, static void getCoverageDependence(const XML_Node& node,
thermo_t& surfphase, ReactionData& rdata) thermo_t& surfphase, ReactionData& rdata)
{ {
@ -507,7 +522,7 @@ static void getCoverageDependence(const XML_Node& node,
} }
} }
} }
//=====================================================================================================
//! Get falloff parameters for a reaction. //! Get falloff parameters for a reaction.
/*! /*!
* This routine reads the falloff XML node and extracts parameters into a * This routine reads the falloff XML node and extracts parameters into a