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

InterfaceKinetics rewrite -> implementation of BV and affinity reactions

Browse source 
next iteration
This commit is contained in:
Harry Moffat 2014-08-26 20:55:38 +00:00
parent 9e05e85b66
commit 979bdffe6e
7 changed files with 586 additions and 58 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

16
include/cantera/kinetics/ElectrodeKinetics.h
View file

@ -68,21 +68,24 @@ public:
*/
void identifyMetalPhase();
//! Internal routine that updates the Rates of Progress of the reactions
/*!
* This is actually the guts of the functionality of the object
*/
virtual void updateROP();
//void addGlobalReaction(ReactionData& r);
double calcForwardROP_BV(size_t irxn, size_t iBeta);
protected:
//! index of the metal phase in the list of phases for this surface
//! Index of the metal phase in the list of phases for this kinetics object
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
@ -91,9 +94,6 @@ protected:
//! Index of the electrons species in the list of species for this surface kinetics, if none set it to -1
size_t kElectronRS_;
};
}

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

@ -24,6 +24,32 @@ class SurfPhase;
class ImplicitSurfChem;
class RxnMolChange;
//! forward orders
class RxnOrders {
public:
//! constructors
RxnOrders();
RxnOrders(const RxnOrders &right);
~RxnOrders();
RxnOrders& operator=(const RxnOrders &right);
//! Fill in the structure with the array.
/*!
* @param[in] Size of length kinetic species. The entries the values of the orders
*/
int fill(const std::vector<doublereal>& fullForwardOrders);
//! ID's of the kinetic species
std::vector<size_t> kinSpeciesIDs_;
//! Orders of the kinetic species
std::vector<doublereal> kinSpeciesOrders_;
};
//! A kinetics manager for heterogeneous reaction mechanisms. The
//! reactions are assumed to occur at a 2D interface between two 3D phases.
/*!
@ -139,28 +165,35 @@ public:
//! @name Reaction Mechanism Informational Query Routines
//! @{
virtual doublereal reactantStoichCoeff(size_t k, size_t i) const {
return m_rrxn[k][i];
}
//! Provide a reactant stoichiometric coefficient
/*!
* @param[in] kSpecKin Species index within the kinetics object
* @param[in] irxn Reaction index
*
* @return Returns the reactant stoichiometic coefficient within the reaction
*/
virtual doublereal reactantStoichCoeff(size_t kSpecKin, size_t irxn) const;
virtual doublereal productStoichCoeff(size_t k, size_t i) const {
return m_prxn[k][i];
}
//! Provide a product stoichiometric coefficient
/*!
* @param[in] kSpecKin Species index within the kinetics object
* @param[in] irxn Reaction index
*
* @return Returns the product stoichiometic coefficient within the reaction
*/
virtual doublereal productStoichCoeff(size_t kSpecKin, size_t irxn) const;
//! return the reaction type of the reaction i
//! return the reaction type of the reaction irxn
/*!
* @param[in] Reaction index
*
* @return Returns the reaction type of the reaction.
*/
virtual int reactionType(size_t i) const;
//virtual int reactionType(size_t i) const {
// return m_index[i].first;
//}
virtual int reactionType(size_t irxn) const;
virtual void getActivityConcentrations(doublereal* const conc);
//! Return the charge transfer rxn Beta parameter for the ith reaction
/*!
* Returns the beta parameter for a charge transfer reaction. This
@ -222,7 +255,7 @@ public:
/*!
* This is actually the guts of the functionality of the object
*/
void updateROP();
virtual void updateROP();
//! Update properties that depend on temperature
/*!
@ -414,6 +447,9 @@ public:
*/
int phaseStability(const size_t iphase) const;
void determineFwdOrdersBV(ReactionData& rdata, std::vector<doublereal>& fwdFullorders);
protected:
//! Temporary work vector of length m_kk
vector_fp m_grt;
@ -663,6 +699,23 @@ protected:
*/
vector_int m_ctrxn_ecdf;
//! Vector of booleans indicating whether the charge transfer reaction rate constant
//! is described by an exchange current density rate constant expression
/*!
* Length is equal to the number of reactions with charge transfer coefficients, m_ctrxn[]
*
* Some reactions have zero in this list, those that don't need special treatment.
*/
std::vector<RxnOrders*> m_ctrxn_ROPOrdersList_;
//! Reaction Orders for the case where the forwards rate of progress is being calculated.
/*!
* Length is equal to the number of reactions with charge transfer coefficients, m_ctrxn[]
*
* Some reactions have zero in this list, indicating that the calculation isn't necessary.
*/
std::vector<RxnOrders*> m_ctrxn_FwdOrdersList_;
//! Vector of standard concentrations
/*!
* Length number of kinetic species

8
include/cantera/kinetics/ReactionData.h
View file

@ -24,6 +24,7 @@ public:
validate(false),
number(0),
rxn_number(0),
filmResistivity(0.0),
reversible(true),
duplicate(false),
rateCoeffType(ARRHENIUS_REACTION_RATECOEFF_TYPE),
@ -102,6 +103,13 @@ public:
//! products.
std::map<int, doublereal> net_stoich;
//! Film Resistivity value
/*!
* Only valid for Butler-Volmer formulations.
* Units are in ohms m2.
*/
double filmResistivity;
//! True if the current reaction is reversible. False otherwise
bool reversible;

15
include/cantera/kinetics/importKinetics.h
View file

@ -87,9 +87,20 @@ bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp, std::string default
std::vector<size_t>& spnum, vector_fp& stoich,
vector_fp& order, const ReactionRules& rules);
//! Install Butler Volmer Orders into the forward orders array.
/*!
* Install the BV order coefficients into the fullForwardsOrders vector.
*
* @param[in] rxnNode XML node pointing to the reaction element in the xml tree.
* @param[in] kin Reference to the kinetics object to install the information into.
* @param[in] rdata Reaction Data Object containing the information about one reaction
* @param[out] fullForwardsOrders Vectors of the orders of reaction.
*
*/
void installButlerVolmerOrders(const XML_Node& rxnNode, const Kinetics& kin, const ReactionData& rdata,
std::vector<doublereal>& fullForwardsOrders);
//! Get non-mass-action orders for a reaction
extern bool getOrders(const XML_Node& rxnNode, Kinetics& kin,
std::string default_phase, const ReactionData& rdata,
vector_fp& order, vector_fp& fullForwardsOrders,

239
src/kinetics/ElectrodeKinetics.cpp
View file

@ -3,6 +3,7 @@
*/
#include "cantera/kinetics/ElectrodeKinetics.h"
#include "cantera/thermo/SurfPhase.h"
using namespace std;
@ -117,6 +118,242 @@ void ElectrodeKinetics::identifyMetalPhase()
}
}
//============================================================================================================================
// virtual from InterfaceKinetics
void ElectrodeKinetics::updateROP()
{
// evaluate rate constants and equilibrium constants at temperature and phi (electric potential)
_update_rates_T();
// get updated activities (rates updated below)
_update_rates_C();
double TT = m_surf->temperature();
double rtdf = GasConstant * TT / Faraday;
if (m_ROP_ok) {
return;
}
//
// Copy the reaction rate coefficients, m_rfn, into m_ropf
//
copy(m_rfn.begin(), m_rfn.end(), m_ropf.begin());
//
// Multiply by the perturbation factor
//
multiply_each(m_ropf.begin(), m_ropf.end(), m_perturb.begin());
//
// Copy the forward rate constants to the reverse rate constants
//
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
// the equilibrium constants
//
multiply_each(m_ropr.begin(), m_ropr.end(), m_rkcn.begin());
//
// multiply ropf by the actyivity concentration reaction orders to obtain
// the forward rates of progress.
//
m_rxnstoich.multiplyReactants(DATA_PTR(m_actConc), DATA_PTR(m_ropf));
//
// For reversible reactions, multiply ropr by the activity concentration products
//
m_rxnstoich.multiplyRevProducts(DATA_PTR(m_actConc), DATA_PTR(m_ropr));
//
// Fix up these calculations for cases where the above formalism doesn't hold
//
double OCV = 0.0;
for (size_t iBeta = 0; iBeta < m_beta.size(); iBeta++) {
size_t irxn = m_ctrxn[iBeta];
int reactionType = reactionTypes_[irxn];
if (reactionType == BUTLERVOLMER_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
//
//
// 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;
//
// Calculate the exchange current density
// m_ropf contains the exchange current reaction rate
//
double io = m_ropf[irxn] * nStoichElectrons;
double exp1 = nu * nStoichElectrons * beta / rtdf;
double exp2 = - nu * nStoichElectrons * (1.0 - beta) / (rtdf);
m_ropnet[irxn] = io * (exp(exp1) - exp(exp2));
// Need to resurrect the forwards rate constant.
//m_ropf[irxn] = ;
m_ropr[irxn] = m_ropnet[irxn] - m_ropf[irxn];
}
}
for (size_t j = 0; j != m_ii; ++j) {
m_ropnet[j] = m_ropf[j] - m_ropr[j];
}
/*
* For reactions involving multiple phases, we must check that the phase
* being consumed actually exists. This is particularly important for
* phases that are stoichiometric phases containing one species with a unity activity
*/
if (m_phaseExistsCheck) {
for (size_t j = 0; j != m_ii; ++j) {
if ((m_ropr[j] > m_ropf[j]) && (m_ropr[j] > 0.0)) {
for (size_t p = 0; p < nPhases(); p++) {
if (m_rxnPhaseIsProduct[j][p]) {
if (! m_phaseExists[p]) {
m_ropnet[j] = 0.0;
m_ropr[j] = m_ropf[j];
if (m_ropf[j] > 0.0) {
for (size_t rp = 0; rp < nPhases(); rp++) {
if (m_rxnPhaseIsReactant[j][rp]) {
if (! m_phaseExists[rp]) {
m_ropnet[j] = 0.0;
m_ropr[j] = m_ropf[j] = 0.0;
}
}
}
}
}
}
if (m_rxnPhaseIsReactant[j][p]) {
if (! m_phaseIsStable[p]) {
m_ropnet[j] = 0.0;
m_ropr[j] = m_ropf[j];
}
}
}
} else if ((m_ropf[j] > m_ropr[j]) && (m_ropf[j] > 0.0)) {
for (size_t p = 0; p < nPhases(); p++) {
if (m_rxnPhaseIsReactant[j][p]) {
if (! m_phaseExists[p]) {
m_ropnet[j] = 0.0;
m_ropf[j] = m_ropr[j];
if (m_ropf[j] > 0.0) {
for (size_t rp = 0; rp < nPhases(); rp++) {
if (m_rxnPhaseIsProduct[j][rp]) {
if (! m_phaseExists[rp]) {
m_ropnet[j] = 0.0;
m_ropf[j] = m_ropr[j] = 0.0;
}
}
}
}
}
}
if (m_rxnPhaseIsProduct[j][p]) {
if (! m_phaseIsStable[p]) {
m_ropnet[j] = 0.0;
m_ropf[j] = m_ropr[j];
}
}
}
}
}
}
m_ROP_ok = true;
}
//==================================================================================================================
//
// 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)
{
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];
//
// 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
//
double tmp = exp(- m_beta[iBeta] * m_deltaG0[irxn] * rrt);
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
//
const RxnOrders* ro_fwd = m_ctrxn_FwdOrdersList_[iBeta];
if (ro_fwd == 0) {
throw CanteraError("ElectrodeKinetics::calcForwardROP_BV()", "forward orders pointer is zero ?!?");
}
tmp = 1.0;
const std::vector<size_t>& kinSpeciesIDs = ro_fwd->kinSpeciesIDs_;
const std::vector<doublereal>& kinSpeciesOrders = ro_fwd->kinSpeciesOrders_;
for (size_t j = 0; j < kinSpeciesIDs.size(); j++) {
size_t k = kinSpeciesIDs[j];
double oo = kinSpeciesOrders[j];
tmp *= pow(m_actConc[k], oo);
}
double ropf = kf * tmp;
return ropf;
}
//==================================================================================================================
//==================================================================================================================
}

211
src/kinetics/InterfaceKinetics.cpp
View file

@ -61,7 +61,12 @@ InterfaceKinetics::~InterfaceKinetics()
for (size_t i = 0; i < rmcVector.size(); i++) {
delete rmcVector[i];
}
for (size_t i = 0; i < m_ctrxn_ROPOrdersList_.size(); i++) {
delete m_ctrxn_ROPOrdersList_[i];
}
for (size_t i = 0; i < m_ctrxn_FwdOrdersList_.size(); i++) {
delete m_ctrxn_FwdOrdersList_[i];
}
}
//============================================================================================================================
InterfaceKinetics::InterfaceKinetics(const InterfaceKinetics& right) :
@ -175,6 +180,24 @@ InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right)
}
}
for (size_t i = 0; i < m_ctrxn_ROPOrdersList_.size(); i++) {
delete m_ctrxn_ROPOrdersList_[i];
}
m_ctrxn_ROPOrdersList_ = right.m_ctrxn_ROPOrdersList_;
for (size_t i = 0; i < m_ctrxn_ROPOrdersList_.size(); i++) {
RxnOrders* ro = right.m_ctrxn_ROPOrdersList_[i];
m_ctrxn_ROPOrdersList_[i] = new RxnOrders(*ro);
}
for (size_t i = 0; i < m_ctrxn_FwdOrdersList_.size(); i++) {
delete m_ctrxn_FwdOrdersList_[i];
}
m_ctrxn_FwdOrdersList_ = right.m_ctrxn_FwdOrdersList_;
for (size_t i = 0; i < m_ctrxn_FwdOrdersList_.size(); i++) {
RxnOrders* ro = right.m_ctrxn_FwdOrdersList_[i];
m_ctrxn_FwdOrdersList_[i] = new RxnOrders(*ro);
}
return *this;
}
@ -390,6 +413,11 @@ void InterfaceKinetics::getEquilibriumConstants(doublereal* kc)
//===========================================================================================================
/*
* values needed to convert from exchange current density to surface reaction rate.
* Calculate:
* - m_StandardConc[]
* - m_ProdStandConcReac[]
* - m_deltaG0[]
* - m_mu0[]
*/
void InterfaceKinetics::updateExchangeCurrentQuantities()
{
@ -419,7 +447,6 @@ void InterfaceKinetics::updateExchangeCurrentQuantities()
m_ProdStanConcReac[i] = 1.0;
}
m_rxnstoich.multiplyReactants(DATA_PTR(m_StandardConc), DATA_PTR(m_ProdStanConcReac));
}
//===========================================================================================================
void InterfaceKinetics::getCreationRates(doublereal* cdot)
@ -476,25 +503,31 @@ void InterfaceKinetics::applyVoltageKfwdCorrection(doublereal* const kf)
#endif
for (size_t i = 0; i < m_beta.size(); i++) {
size_t irxn = m_ctrxn[i];
eamod = m_beta[i] * deltaElectricEnergy_[irxn];
if (eamod != 0.0) {
//
// If we calculate the BV form directly, we don't add the voltage correction to the
// forward reaction rate constants.
//
if (m_ctrxn_BVform[i] == 0) {
eamod = m_beta[i] * deltaElectricEnergy_[irxn];
if (eamod != 0.0) {
#ifdef DEBUG_KIN_MODE
ea = GasConstant * m_E[irxn];
if (eamod + ea < 0.0) {
writelog("Warning: act energy mod too large!\n");
writelog(" Delta phi = "+fp2str(deltaElectricEnergy_[irxn]/Faraday)+"\n");
writelog(" Delta Ea = "+fp2str(eamod)+"\n");
writelog(" Ea = "+fp2str(ea)+"\n");
for (n = 0; n < np; n++) {
writelog("Phase "+int2str(n)+": phi = "
+fp2str(m_phi[n])+"\n");
}
}
ea = GasConstant * m_E[irxn];
if (eamod + ea < 0.0) {
writelog("Warning: act energy mod too large!\n");
writelog(" Delta phi = "+fp2str(deltaElectricEnergy_[irxn]/Faraday)+"\n");
writelog(" Delta Ea = "+fp2str(eamod)+"\n");
writelog(" Ea = "+fp2str(ea)+"\n");
for (n = 0; n < np; n++) {
writelog("Phase "+int2str(n)+": phi = "
+fp2str(m_phi[n])+"\n");
}
}
#endif
doublereal rt = GasConstant*thermo(0).temperature();
doublereal rrt = 1.0/rt;
kf[irxn] *= exp(-eamod*rrt);
}
doublereal rt = GasConstant*thermo(0).temperature();
doublereal rrt = 1.0/rt;
kf[irxn] *= exp(-eamod*rrt);
}
}
}
}
//==================================================================================================================
@ -525,8 +558,10 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con
int iECDFormulation = m_ctrxn_ecdf[i];
if (iECDFormulation) {
//
// 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 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.
//
// We need to have the straight chemical reaction rate constant to come out of this calculation.
//
if (m_ctrxn_BVform[i] == 0) {
//
@ -537,7 +572,12 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con
tmp *= 1.0 / tmp2 / Faraday;
kfwd[irxn] *= tmp;
}
//
// If BVform is nonzero we don't need to do anything.
//
} else {
//
// kfwd[] is the chemical reaction rate constant
//
// 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
@ -546,7 +586,7 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con
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
// it's in the exchange current density formulation format
//
double tmp = exp(m_beta[i] * m_deltaG0[irxn] * rrt);
double tmp2 = m_ProdStanConcReac[irxn];
@ -930,16 +970,20 @@ void InterfaceKinetics::addReaction(ReactionData& r)
void InterfaceKinetics::addElementaryReaction(ReactionData& rdata)
{
// install rate coeff calculator
//
// install rate coefficient calculator
//
vector_fp& rp = rdata.rateCoeffParameters;
size_t ncov = rdata.cov.size();
//
// Turn on the global flag indicating surface coverage dependence
//
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
//
@ -981,6 +1025,8 @@ void InterfaceKinetics::addElementaryReaction(ReactionData& rdata)
} else {
m_ctrxn_ecdf.push_back(0);
}
m_ctrxn_ROPOrdersList_.push_back(0);
m_ctrxn_FwdOrdersList_.push_back(0);
}
// add constant term to rate coeff value vector
@ -1029,14 +1075,15 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata)
* Change the reaction rate coefficient type back to its original value
*/
rdata.rateCoeffType = reactionRateCoeffType_orig;
// store activation energy
//
// Store activation energy
//
m_E.push_back(rdata.rateCoeffParameters[2]);
//
// Add the reaction into the list of electrochemical extras
//
if (rdata.beta > 0.0) {
if (rdata.beta > 0.0 || 1) {
m_has_electrochem_rxns = true;
m_beta.push_back(rdata.beta);
// Push back the id of the reaction
@ -1049,6 +1096,27 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata)
} else {
m_ctrxn_ecdf.push_back(0);
}
if (rdata.forwardFullOrder_.size() > 0) {
RxnOrders* ro = new RxnOrders();
ro->fill(rdata.forwardFullOrder_);
m_ctrxn_ROPOrdersList_.push_back(ro);
m_ctrxn_FwdOrdersList_.push_back(0);
//
//
// Fill in the Fwd Orders dependence here for B-V reactions
//
if (rdata.reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN || rdata.reactionType == BUTLERVOLMER_RXN) {
std::vector<double> fwdFullorders(m_kk, 0.0);
determineFwdOrdersBV(rdata, fwdFullorders);
RxnOrders* ro = new RxnOrders();
ro->fill(rdata.forwardFullOrder_);
m_ctrxn_FwdOrdersList_[m_ii] = ro;
}
} else {
m_ctrxn_ROPOrdersList_.push_back(0);
m_ctrxn_FwdOrdersList_.push_back(0);
}
}
// add constant term to rate coeff value vector
@ -1303,9 +1371,19 @@ int InterfaceKinetics::phaseStability(const size_t iphase) const
return m_phaseIsStable[iphase];
}
//==================================================================================================================
int InterfaceKinetics::reactionType(size_t i) const
doublereal InterfaceKinetics::reactantStoichCoeff(size_t kSpecKin, size_t irxn) const
{
return reactionType_[i];
return m_rrxn[kSpecKin][irxn];
}
//==================================================================================================================
doublereal InterfaceKinetics::productStoichCoeff(size_t kSpecKin, size_t irxn) const
{
return m_prxn[kSpecKin][irxn];
}
//==================================================================================================================
int InterfaceKinetics::reactionType(size_t irxn) const
{
return reactionType_[irxn];
}
//==================================================================================================================
void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStable)
@ -1334,6 +1412,41 @@ void InterfaceKinetics::registerReaction(size_t rxnNumber, int type, size_t loc)
m_index[rxnNumber] = std::pair<int, size_t>(type, loc);
}
//==================================================================================================================
//
void InterfaceKinetics::determineFwdOrdersBV(ReactionData& rdata, std::vector<doublereal>& 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[kkin];
fwdFullorders[kkin] += betaf * oo;
if (abs(fwdFullorders[kkin]) < 0.00001) {
fwdFullorders[kkin] = 0.0;
}
}
for (size_t j = 0; j < rdata.products.size(); j++) {
size_t kkin = rdata.products[j];
double oo = rdata.pstoich[kkin];
fwdFullorders[kkin] -= betaf * oo;
if (abs(fwdFullorders[kkin]) < 0.00001) {
fwdFullorders[kkin] = 0.0;
}
}
}
//==================================================================================================================
void EdgeKinetics::finalize()
{
deltaElectricEnergy_.resize(std::max<size_t>(m_ii, 1));
@ -1366,6 +1479,44 @@ void EdgeKinetics::finalize()
m_finalized = true;
}
//==================================================================================================================
RxnOrders::RxnOrders()
{
}
//==================================================================================================================
RxnOrders::~RxnOrders()
{
}
//==================================================================================================================
RxnOrders::RxnOrders(const RxnOrders& right) :
kinSpeciesIDs_(right.kinSpeciesIDs_),
kinSpeciesOrders_(right.kinSpeciesOrders_)
{
}
//==================================================================================================================
RxnOrders& RxnOrders::operator=(const RxnOrders& right)
{
if (this == &right) {
return *this;
}
kinSpeciesIDs_ = right.kinSpeciesIDs_;
kinSpeciesOrders_ = right.kinSpeciesOrders_;
return *this;
}
//==================================================================================================================
int RxnOrders::fill(const std::vector<doublereal>& fullForwardOrders)
{
int nzeroes = 0;
kinSpeciesIDs_.clear();
kinSpeciesOrders_.clear();
for (size_t k = 0; k < fullForwardOrders.size(); ++k) {
if (fullForwardOrders[k] != 0.0) {
kinSpeciesIDs_.push_back(k);
kinSpeciesOrders_.push_back(fullForwardOrders[k]);
++nzeroes;
}
}
return nzeroes;
}
//==================================================================================================================
}

76
src/kinetics/importKinetics.cpp
View file

@ -243,6 +243,56 @@ bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp,
return true;
}
//====================================================================================================================
//
// Install the BV order coefficients into the fullForwardsOrders vector.
//
void installButlerVolmerOrders(const XML_Node& rxnNode, const Kinetics& kin, const ReactionData& rdata,
std::vector<doublereal>& fullForwardsOrders)
{
const std::vector<size_t>& reactants = rdata.reactants;
const std::vector<size_t>& products = rdata.products;
const std::vector<doublereal>& rstoich = rdata.rstoich;
const std::vector<doublereal>& pstoich = rdata.pstoich;
//
// Gather the number of species in the kinetics object and resize fullForwardsOrders
//
size_t nsp = kin.nTotalSpecies();
fullForwardsOrders.resize(nsp, 0.0);
//
// Ok first thing to do is get the electrochemical transfer coefficient
// since the order depend on the value.
// Also, if we don't find one, then it's an error. Zero is an acceptable value.
// Beta below 0 or greater than 1 are probably not good.
//
double beta = -10.0;
if (rxnNode.hasChild("rateCoeff")) {
XML_Node& rc = rxnNode.child("rateCoeff");
if (rc.hasChild("electrochem")) {
XML_Node& eb = rc.child("electrochem");
string sbeta = eb["beta"];
beta = fpValueCheck(sbeta);
}
}
if (beta == -10.0) {
throw CanteraError("installButlerVolmerOrders()",
"ButlerVolmerOrders model requested but no electrochem beta input");
}
double betar = 1.0 - beta;
for (size_t k = 0; k < nsp; k++) {
fullForwardsOrders[k] = 0.0;
}
for (size_t n = 0; n < reactants.size(); n++) {
size_t k = reactants[n];
double fac = rstoich[n];
fullForwardsOrders[k] += fac * betar;
}
for (size_t n = 0; n < products.size(); n++) {
size_t k = products[n];
double fac = pstoich[n];
fullForwardsOrders[k] += fac * beta;
}
}
//====================================================================================================================
// Fill in the fullForwardsOrders array for a specific reaction
/*
* rxnNode XML node for the reaction
@ -825,12 +875,13 @@ bool rxninfo::installReaction(int iRxn, const XML_Node& rxnNode, Kinetics& kin,
}
}
//
// get the reactant and their stoichiometries
// Get the reactant and their stoichiometries
//
bool ok = getReagents(rxnNode, kin, 1, default_phase, rdata.reactants,
rdata.rstoich, rdata.rorder, rules);
//
// Get the products. We store the id of products in rdata.products
//
ok = ok && getReagents(rxnNode, kin, -1, default_phase, rdata.products,
rdata.pstoich, rdata.porder, rules);
@ -863,9 +914,26 @@ bool rxninfo::installReaction(int iRxn, const XML_Node& rxnNode, Kinetics& kin,
}
rdata.global = true;
}
//
// Fill in the forwardFullOrder_ array
// For Butler Volmer reactions, we'll install the orders for the exchange current into the
// forwardFullOrders array. It may be altered by the getOrders function below.
//
if (rdata.reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN || rdata.reactionType == BUTLERVOLMER_RXN) {
if (! rdata.reversible) {
throw CanteraError("installReaction()", "a Butler-Volmer rxn must be reversible");
}
installButlerVolmerOrders(rxnNode, kin, rdata, rdata.forwardFullOrder_);
//
// For Butler Volmer reactions, a common addition to the formulation is to add an electrical resistance
// to the formulation. The resistance modifies the electrical current flow in both directions
//
if (rxnNode.hasChild("filmResistivity")) {
XML_Node& fNode = rxnNode.child("filmResistivity");
rdata.filmResistivity = fpValueCheck( fNode() );
}
}
//
// Fill in the forwardFullOrder_ array
//
if (rxnNode.hasChild("orders")) {
ok = getOrders(rxnNode, kin, default_phase, rdata,