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

Two major capabilities have been merged into the main branch

Browse source 
from liquidTransportDevelop branch

The Exchange_current_reaction_ratecoeff_type method for
formulating reaction rate constants for exchange current
reactions has been ported into the main branch.

Fixed the InterfaceKinetics object so that we can
specify that a phase doesn't exist. If a phase doesn't exist,
species in that phase can't be reactants in the phase. Before,
this fix, they could be !?!!
This was because InterfaceKinetics objects use ThermoPhase
objects, which are intrinsic only. We need an extrinsic notion 
that is in addition to ThermoPhase info in order to get this right.
This commit is contained in:
Harry Moffat 2010-05-08 03:15:22 +00:00
parent c292e98c8a
commit eaa2875d10
9 changed files with 794 additions and 253 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

368
Cantera/src/kinetics/InterfaceKinetics.cpp
View file

@ -25,9 +25,8 @@ using namespace std;
namespace Cantera {
//////////////////////////////////////////////////////////////////
/**
//====================================================================================================================
/*
* Construct an empty InterfaceKinetics reaction mechanism.
* @param thermo This is an optional parameter that may be
* used to initialize the inherited Kinetics class with
@ -35,8 +34,7 @@ namespace Cantera {
* useful for initialization of homogeneous kinetics
* mechanisms.
*/
InterfaceKinetics::
InterfaceKinetics(thermo_t* thermo) :
InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) :
Kinetics(),
m_kk(0),
m_redo_rates(false),
@ -44,27 +42,41 @@ namespace Cantera {
m_nrev(0),
m_surf(0),
m_integrator(0),
m_beta(0),
m_ctrxn(0),
m_ctrxn_ecdf(0),
m_StandardConc(0),
m_deltaG0(0),
m_ProdStanConcReac(0),
m_finalized(false),
m_has_coverage_dependence(false),
m_has_electrochem_rxns(false),
m_has_exchange_current_density_formulation(false),
m_phaseExistsCheck(false),
m_phaseExists(0),
m_rxnPhaseIsReactant(0),
m_rxnPhaseIsProduct(0),
m_ioFlag(0)
{
if (thermo != 0) addPhase(*thermo);
m_kdata = new InterfaceKineticsData;
m_kdata->m_temp = 0.0;
}
/**
//====================================================================================================================
/*
* Destructor
*/
InterfaceKinetics::
~InterfaceKinetics(){
InterfaceKinetics::~InterfaceKinetics(){
delete m_kdata;
if (m_integrator) {
delete m_integrator;
}
for (int i = 0; i < m_ii; i++) {
delete m_rxnPhaseIsReactant[i];
delete m_rxnPhaseIsProduct[i];
}
}
//====================================================================================================================
// Copy Constructor for the %InterfaceKinetics object.
/*
* Currently, this is not fully implemented. If called it will
@ -78,9 +90,20 @@ namespace Cantera {
m_nrev(0),
m_surf(0),
m_integrator(0),
m_beta(0),
m_ctrxn(0),
m_ctrxn_ecdf(0),
m_StandardConc(0),
m_deltaG0(0),
m_ProdStanConcReac(0),
m_finalized(false),
m_has_coverage_dependence(false),
m_has_electrochem_rxns(false),
m_has_exchange_current_density_formulation(false),
m_phaseExistsCheck(false),
m_phaseExists(0),
m_rxnPhaseIsReactant(0),
m_rxnPhaseIsProduct(0),
m_ioFlag(0)
{
m_kdata = new InterfaceKineticsData;
@ -90,7 +113,7 @@ namespace Cantera {
*/
*this = operator=(right);
}
//====================================================================================================================
// Assignment operator
/*
* This is NOT a virtual function.
@ -100,13 +123,21 @@ namespace Cantera {
*/
InterfaceKinetics& InterfaceKinetics::
operator=(const InterfaceKinetics &right) {
int i;
/*
* Check for self assignment.
*/
if (this == &right) return *this;
for (i = 0; i < m_ii; i++) {
delete (m_rxnPhaseIsReactant[i]);
delete (m_rxnPhaseIsProduct[i]);
}
Kinetics::operator=(right);
m_kk = right.m_kk;
m_revindex = right.m_revindex;
m_rates = right.m_rates;
@ -129,17 +160,58 @@ namespace Cantera {
m_integrator = right.m_integrator; //DANGER - shallow copy
m_beta = right.m_beta;
m_ctrxn = right.m_ctrxn;
m_ctrxn_ecdf = right.m_ctrxn_ecdf;
m_StandardConc = right.m_StandardConc;
m_deltaG0 = right.m_deltaG0;
m_ProdStanConcReac = right.m_ProdStanConcReac;
m_finalized = right.m_finalized;
m_has_coverage_dependence = right.m_has_coverage_dependence;
m_has_electrochem_rxns = right.m_has_electrochem_rxns;
m_has_exchange_current_density_formulation = right.m_has_exchange_current_density_formulation;
m_phaseExistsCheck = right.m_phaseExistsCheck;
m_phaseExists = right.m_phaseExists;
m_rxnPhaseIsReactant.resize(m_ii, 0);
m_rxnPhaseIsProduct.resize(m_ii, 0);
int np = nPhases();
for (i = 0; i < m_ii; i++) {
m_rxnPhaseIsReactant[i] = new bool[np];
m_rxnPhaseIsProduct[i] = new bool[np];
for (int p = 0; p < np; p++) {
m_rxnPhaseIsReactant[i][p] = right.m_rxnPhaseIsReactant[i][p];
m_rxnPhaseIsProduct[i][p] = right.m_rxnPhaseIsProduct[i][p];
}
}
m_rxnPhaseIsProduct = right.m_rxnPhaseIsProduct;
m_ioFlag = right.m_ioFlag;
return *this;
}
//====================================================================================================================
// Return the ID of the kinetics object
int InterfaceKinetics::ID() const {
return cInterfaceKinetics;
}
//====================================================================================================================
int InterfaceKinetics::type() const {
return cInterfaceKinetics;
}
//====================================================================================================================
// Set the electric potential in the nth phase
/*
* @param n phase Index in this kinetics object.
* @param V Electric potential (volts)
*/
void InterfaceKinetics::setElectricPotential(int n, doublereal V) {
thermo(n).setElectricPotential(V);
m_redo_rates = true;
}
// Duplication routine for objects which inherit from
// Kinetics
//====================================================================================================================
// Duplication routine for objects which inherit from Kinetics
/*
* This virtual routine can be used to duplicate %Kinetics objects
* inherited from %Kinetics even if the application only has
@ -152,14 +224,18 @@ namespace Cantera {
InterfaceKinetics* tp = new InterfaceKinetics(*this);
return dynamic_cast<Kinetics *>(tp);
}
/**
* Update properties that depend on temperature
//====================================================================================================================
// Update properties that depend on temperature
/*
* This is called to update all of the properties that depend on temperature
*
*/
void InterfaceKinetics::
_update_rates_T() {
* Current objects that this function updates
* m_kdata->m_logtemp
* m_kdata->m_rfn
* m_rates.
* updateKc();
*/
void InterfaceKinetics::_update_rates_T() {
_update_rates_phi();
if (m_has_coverage_dependence) {
m_surf->getCoverages(DATA_PTR(m_conc));
@ -167,20 +243,24 @@ namespace Cantera {
m_redo_rates = true;
}
doublereal T = thermo(surfacePhaseIndex()).temperature();
m_redo_rates = true;
if (T != m_kdata->m_temp || m_redo_rates) {
m_kdata->m_logtemp = log(T);
m_rates.update(T, m_kdata->m_logtemp, DATA_PTR(m_kdata->m_rfn));
if (m_has_electrochem_rxns)
if (m_has_exchange_current_density_formulation) {
applyExchangeCurrentDensityFormulation(DATA_PTR(m_kdata->m_rfn));
}
if (m_has_electrochem_rxns) {
applyButlerVolmerCorrection(DATA_PTR(m_kdata->m_rfn));
}
m_kdata->m_temp = T;
updateKc();
m_kdata->m_ROP_ok = false;
m_redo_rates = false;
}
}
void InterfaceKinetics::
_update_rates_phi() {
//====================================================================================================================
void InterfaceKinetics::_update_rates_phi() {
int np = nPhases();
for (int n = 0; n < np; n++) {
if (thermo(n).electricPotential() != m_phi[n]) {
@ -189,6 +269,7 @@ namespace Cantera {
}
}
}
//====================================================================================================================
/**
@ -198,8 +279,7 @@ namespace Cantera {
* representing phases should overload to return the appropriate
* quantities.
*/
void InterfaceKinetics::
_update_rates_C() {
void InterfaceKinetics::_update_rates_C() {
int n;
int np = nPhases();
@ -263,6 +343,8 @@ namespace Cantera {
}
}
}
//====================================================================================================================
void InterfaceKinetics::checkPartialEquil() {
@ -338,6 +420,36 @@ namespace Cantera {
}
}
void InterfaceKinetics::getExchangeCurrentQuantities() {
/*
* First collect vectors of the standard Gibbs free energies of the
* species and the standard concentrations
* - m_mu0
* - m_logStandardConc
*/
int ik = 0;
int np = nPhases();
for (int n = 0; n < np; n++) {
thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]);
int nsp = thermo(n).nSpecies();
for (int k = 0; k < nsp; k++) {
m_StandardConc[ik] = thermo(n).standardConcentration(k);
ik++;
}
}
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), DATA_PTR(m_deltaG0));
for (int i = 0; i < m_ii; i++) {
m_ProdStanConcReac[i] = 1.0;
}
m_rxnstoich.multiplyReactants(DATA_PTR(m_StandardConc), DATA_PTR(m_ProdStanConcReac));
}
// Returns the Species creation rates [kmol/m^2/s].
/*
* Return the species
@ -385,12 +497,17 @@ namespace Cantera {
net);
}
/**
//====================================================================================================================
// Apply corrections for interfacial charge transfer reactions
/*
* For reactions that transfer charge across a potential difference,
* the activation energies are modified by the potential difference.
* (see, for example, ...). This method applies this correction.
*
* @param kf Vector of forward reaction rate constants on which to have
* the correction applied
*/
void InterfaceKinetics::applyButlerVolmerCorrection(doublereal* kf) {
void InterfaceKinetics::applyButlerVolmerCorrection(doublereal* const kf) {
int i;
int n, nsp, k, ik=0;
@ -410,8 +527,7 @@ namespace Cantera {
// Compute the change in electrical potential energy for each
// reaction. This will only be non-zero if a potential
// 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(m_rwork));
// Modify the reaction rates. Only modify those with a
// non-zero activation energy. Below we decrease the
@ -432,7 +548,8 @@ namespace Cantera {
for (i = 0; i < nct; i++) {
irxn = m_ctrxn[i];
eamod = m_beta[i]*m_rwork[irxn];
if (eamod != 0.0 && m_E[irxn] != 0.0) {
// if (eamod != 0.0 && m_E[irxn] != 0.0) {
if (eamod != 0.0) {
#ifdef DEBUG_KIN_MODE
ea = GasConstant * m_E[irxn];
if (eamod + ea < 0.0) {
@ -450,18 +567,32 @@ namespace Cantera {
}
}
}
//====================================================================================================================
void InterfaceKinetics::applyExchangeCurrentDensityFormulation(doublereal* const kfwd) {
getExchangeCurrentQuantities();
int nct = m_ctrxn.size();
doublereal rt = GasConstant*thermo(0).temperature();
doublereal rrt = 1.0/rt;
for (int i = 0; i < nct; i++) {
int irxn = m_ctrxn[i];
int iECDFormulation = m_ctrxn_ecdf[i];
if (iECDFormulation) {
double tmp = exp(- m_beta[i] * m_deltaG0[irxn] * rrt);
double tmp2 = m_ProdStanConcReac[irxn];
tmp *= 1.0 / tmp2 / Faraday;
kfwd[irxn] *= tmp;
}
}
}
//====================================================================================================================
/**
* Update the rates of progress of the reactions in the reaciton
* mechanism. This routine operates on internal data.
*/
void InterfaceKinetics::getFwdRateConstants(doublereal* kfwd) {
// _update_rates_T();
// _update_rates_C();
updateROP();
const vector_fp& rf = m_kdata->m_rfn;
@ -473,7 +604,7 @@ namespace Cantera {
multiply_each(kfwd, kfwd + nReactions(), m_perturb.begin());
}
//====================================================================================================================
/**
* Update the rates of progress of the reactions in the reaciton
@ -493,12 +624,12 @@ namespace Cantera {
multiply_each(krev, krev + nReactions(), rkc.begin());
}
}
//====================================================================================================================
void InterfaceKinetics::getActivationEnergies(doublereal *E) {
copy(m_E.begin(), m_E.end(), E);
}
//====================================================================================================================
/**
* Update the rates of progress of the reactions in the reaction
* mechanism. This routine operates on internal data.
@ -545,7 +676,57 @@ namespace Cantera {
for (int j = 0; j != m_ii; ++j) {
ropnet[j] = ropf[j] - 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 (int j = 0; j != m_ii; ++j) {
if ((ropr[j] > ropf[j]) && (ropr[j] > 0.0)) {
for (int p = 0; p < nPhases(); p++) {
if (m_rxnPhaseIsProduct[j][p]) {
if (! m_phaseExists[p]) {
ropnet[j] = 0.0;
ropr[j] = ropf[j];
if (ropf[j] > 0.0) {
for (int rp = 0; rp < nPhases(); rp++) {
if (m_rxnPhaseIsReactant[j][rp]) {
if (! m_phaseExists[rp]) {
ropnet[j] = 0.0;
ropr[j] = ropf[j] = 0.0;;
}
}
}
}
}
}
}
} else if ((ropf[j] > ropr[j]) && (ropf[j] > 0.0)) {
for (int p = 0; p < nPhases(); p++) {
if (m_rxnPhaseIsReactant[j][p]) {
if (! m_phaseExists[p]) {
ropnet[j] = 0.0;
ropf[j] = ropr[j];
if (ropf[j] > 0.0) {
for (int rp = 0; rp < nPhases(); rp++) {
if (m_rxnPhaseIsProduct[j][rp]) {
if (! m_phaseExists[rp]) {
ropnet[j] = 0.0;
ropf[j] = ropr[j] = 0.0;
}
}
}
}
}
}
}
}
}
}
m_kdata->m_ROP_ok = true;
}
@ -738,7 +919,7 @@ namespace Cantera {
m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaS);
}
//====================================================================================================================
/**
* Add a single reaction to the mechanism. This routine
* must be called after init() and before finalize().
@ -752,21 +933,44 @@ namespace Cantera {
* There is no difference between elementary and surface
* reactions.
*/
void InterfaceKinetics::
addReaction(const ReactionData& r) {
void InterfaceKinetics::addReaction(const ReactionData& r) {
addElementaryReaction(r);
// operations common to all reaction types
installReagents( r );
installReagents(r);
//installGroups(reactionNumber(), r.rgroups, r.pgroups);
incrementRxnCount();
m_rxneqn.push_back(r.equation);
}
m_rxnPhaseIsReactant.resize(m_ii, 0);
m_rxnPhaseIsProduct.resize(m_ii, 0);
int np = nPhases();
int i = m_ii -1;
m_rxnPhaseIsReactant[i] = new bool[np];
m_rxnPhaseIsProduct[i] = new bool[np];
void InterfaceKinetics::
addElementaryReaction(const ReactionData& r) {
for (int p = 0; p < np; p++) {
m_rxnPhaseIsReactant[i][p] = false;
m_rxnPhaseIsProduct[i][p] = false;
}
const vector_int& vr = reactants(i);
for (int ik = 0; ik < (int) vr.size(); ik++) {
int k = vr[ik];
int p = speciesPhaseIndex(k);
m_rxnPhaseIsReactant[i][p] = true;
}
const vector_int& vp = products(i);
for (int ik = 0; ik < (int) vp.size(); ik++) {
int k = vp[ik];
int p = speciesPhaseIndex(k);
m_rxnPhaseIsProduct[i][p] = true;
}
}
//====================================================================================================================
void InterfaceKinetics::addElementaryReaction(const ReactionData& r) {
int iloc;
// install rate coeff calculator
@ -776,12 +980,11 @@ namespace Cantera {
if (ncov > 3) {
m_has_coverage_dependence = true;
}
for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]);
iloc = m_rates.install( reactionNumber(),
r.rateCoeffType, rp.size(),
DATA_PTR(rp) );
for (int m = 0; m < ncov; m++) {
rp.push_back(r.cov[m]);
}
// iloc = m_rates.install(reactionNumber(), r.rateCoeffType, rp.size(), DATA_PTR(rp));
iloc = m_rates.install(reactionNumber(), ARRHENIUS_REACTION_RATECOEFF_TYPE, rp.size(), DATA_PTR(rp));
// store activation energy
m_E.push_back(r.rateCoeffParameters[2]);
@ -789,13 +992,19 @@ namespace Cantera {
m_has_electrochem_rxns = true;
m_beta.push_back(r.beta);
m_ctrxn.push_back(reactionNumber());
if (r.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_kdata->m_rfn.push_back(r.rateCoeffParameters[0]);
registerReaction( reactionNumber(), ELEMENTARY_RXN, iloc);
}
//====================================================================================================================
void InterfaceKinetics::setIOFlag(int ioFlag) {
m_ioFlag = ioFlag;
@ -912,7 +1121,7 @@ namespace Cantera {
* calculates rates of species production from reaction rates of
* progress.
*/
m_rxnstoich.add( reactionNumber(), r);
m_rxnstoich.add(reactionNumber(), r);
/*
* register reaction in lists of reversible and irreversible rxns.
*/
@ -920,11 +1129,16 @@ namespace Cantera {
m_revindex.push_back(reactionNumber());
m_nrev++;
} else {
m_irrev.push_back( reactionNumber() );
m_irrev.push_back(reactionNumber());
m_nirrev++;
}
}
//===============================================================================================
void InterfaceKinetics::addPhase(thermo_t &thermo) {
Kinetics::addPhase(thermo);
m_phaseExists.push_back(true);
}
//================================================================================================
/**
* Prepare the class for the addition of reactions. This function
* must be called after instantiation of the class, but before
@ -948,7 +1162,7 @@ namespace Cantera {
m_pot.resize(m_kk, 0.0);
m_phi.resize(np, 0.0);
}
//================================================================================================
/**
* Finish adding reactions and prepare for use. This function
* must be called after all reactions are entered into the mechanism
@ -958,6 +1172,7 @@ namespace Cantera {
* since we don't know this number up to now.
*/
void InterfaceKinetics::finalize() {
Kinetics::finalize();
m_rwork.resize(nReactions());
int ks = reactionPhaseIndex();
if (ks < 0) throw CanteraError("InterfaceKinetics::finalize",
@ -967,10 +1182,20 @@ namespace Cantera {
throw CanteraError("InterfaceKinetics::finalize",
"expected interface dimension = 2, but got dimension = "
+int2str(m_surf->nDim()));
m_StandardConc.resize(m_nTotalSpecies, 0.0);
m_deltaG0.resize(m_ii, 0.0);
m_ProdStanConcReac.resize(m_ii, 0.0);
if (m_thermo.size() != m_phaseExists.size()) {
throw CanteraError("InterfaceKinetics::finalize", "internal error");
}
m_finalized = true;
}
doublereal InterfaceKinetics::electrochem_beta(int irxn) const{
int n = m_ctrxn.size();
for (int i = 0; i < n; i++) {
@ -985,7 +1210,7 @@ namespace Cantera {
bool InterfaceKinetics::ready() const {
return (m_finalized);
}
//================================================================================================
// Advance the surface coverages in time
/*
* @param tstep Time value to advance the surface coverages
@ -1002,7 +1227,7 @@ namespace Cantera {
delete m_integrator;
m_integrator = 0;
}
//================================================================================================
// Solve for the pseudo steady-state of the surface problem
/*
* Solve for the steady state of the surface problem.
@ -1031,7 +1256,25 @@ namespace Cantera {
*/
m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride);
}
//================================================================================================
void InterfaceKinetics::setPhaseExistence(const int iphase, const bool exists) {
if (iphase < 0 || iphase >= (int) m_thermo.size()) {
throw CanteraError("InterfaceKinetics:setPhaseExistence", "out of bounds");
}
if (exists) {
if (!m_phaseExists[iphase]) {
m_phaseExistsCheck--;
m_phaseExists[iphase] = true;
}
} else {
if (m_phaseExists[iphase]) {
m_phaseExistsCheck++;
m_phaseExists[iphase] = false;
}
}
}
//================================================================================================
void EdgeKinetics::finalize() {
m_rwork.resize(nReactions());
int ks = reactionPhaseIndex();
@ -1043,7 +1286,8 @@ namespace Cantera {
"expected interface dimension = 1, but got dimension = "
+int2str(m_surf->nDim()));
m_finalized = true;
}
}
//================================================================================================
}

202
Cantera/src/kinetics/InterfaceKinetics.h
View file

@ -45,11 +45,16 @@ namespace Cantera {
class InterfaceKineticsData {
public:
InterfaceKineticsData() :
m_logp0(0.0),
m_logc0(0.0),
m_ROP_ok(false),
m_temp(0.0), m_logtemp(0.0)
m_temp(0.0),
m_logtemp(0.0)
{}
//! Virtual destructor
virtual ~InterfaceKineticsData(){}
virtual ~InterfaceKineticsData() {
}
doublereal m_logp0;
doublereal m_logc0;
@ -59,7 +64,9 @@ namespace Cantera {
bool m_ROP_ok;
//! Current temperature of the data
doublereal m_temp;
//! Current log of the temperature
doublereal m_logtemp;
vector_fp m_rfn;
vector_fp m_rkcn;
@ -76,9 +83,9 @@ namespace Cantera {
public:
/**
* Constructor
*
//! Constructor
/*!
* @param thermo The optional parameter may be used to initialize
* the object with one ThermoPhase object.
* HKM Note -> Since the interface kinetics
@ -121,20 +128,18 @@ namespace Cantera {
*/
virtual Kinetics *duplMyselfAsKinetics() const;
//! Return the ID of the kinetics object
virtual int ID() const;
virtual int ID() const { return cInterfaceKinetics; }
virtual int type() const { return cInterfaceKinetics; }
//! Retunr the type of the kinetics object
virtual int type() const;
/**
* Set the electric potential in the nth phase
*
//! Set the electric potential in the nth phase
/*!
* @param n phase Index in this kinetics object.
* @param V Electric potential (volts)
*/
void setElectricPotential(int n, doublereal V) {
thermo(n).setElectricPotential(V);
m_redo_rates = true;
}
void setElectricPotential(int n, doublereal V);
///
@ -172,8 +177,17 @@ namespace Cantera {
std::copy(m_kdata->m_ropnet.begin(), m_kdata->m_ropnet.end(), netROP);
}
//! Get the equilibrium constants of all reactions, whether
//! the reaction is reversible or not.
/*!
* @param kc Returns the concentration equation constant for the reaction.
* Length is the number of reactions
*/
virtual void getEquilibriumConstants(doublereal* kc);
void getExchangeCurrentQuantities();
virtual void getDeltaGibbs( doublereal* deltaG);
@ -361,6 +375,19 @@ namespace Cantera {
*/
//@{
//! Add a phase to the kinetics manager object.
/*!
* This must be done before the function init() is called or
* before any reactions are input.
*
* This function calls the Kinetics operator addPhase.
* It also sets the following functions
*
* m_phaseExists[]
*
* @param thermo Reference to the ThermoPhase to be added.
*/
virtual void addPhase(thermo_t& thermo);
//! Prepare the class for the addition of reactions.
/*!
@ -390,12 +417,26 @@ namespace Cantera {
virtual bool ready() const;
//! Internal routine that updates the Rates of Progress of the reactions
/*!
* This is actually the guts of the functionality of the object
*/
void updateROP();
//! Update properties that depend on temperature
/*!
* This is called to update all of the properties that depend on temperature
*
* Current objects that this function updates
* m_kdata->m_logtemp
* m_kdata->m_rfn
* m_rates.
* updateKc();
*/
void _update_rates_T();
void _update_rates_phi();
void _update_rates_C();
@ -445,15 +486,63 @@ namespace Cantera {
void setIOFlag(int ioFlag);
void checkPartialEquil();
int reactionNumber() const { return m_ii;}
void addElementaryReaction(const ReactionData& r);
void addGlobalReaction(const ReactionData& r);
void installReagents(const ReactionData& r);
void updateKc();
//! Write values into m_index
/*!
* @param rxnNumber reaction number
* @param type reaction type
* @param loc location ??
*/
void registerReaction(int rxnNumber, int type, int loc) {
m_index[rxnNumber] = std::pair<int, int>(type, loc);
}
//! Apply corrections for interfacial charge transfer reactions
/*!
* For reactions that transfer charge across a potential difference,
* the activation energies are modified by the potential difference.
* (see, for example, ...). This method applies this correction.
*
* @param kf Vector of forward reaction rate constants on which to have
* the correction applied
*/
void applyButlerVolmerCorrection(doublereal* const kf);
//! When an electrode reaction rate is optionally specified in terms of its
//! exchange current density, extra vectors need to be precalculated
/*!
*
*/
void applyExchangeCurrentDensityFormulation(doublereal* const kfwd);
//! Set the existence of a phase in the reaction object
/*!
* Tell the kinetics object whether a phase in the object exists.
* This is actually an extrinsic specification that must be carried out on top of the
* intrinsic calculation of the reaction rate
*
* @param iphase Index of the phase. This is the order within the internal thermo vector object
* @param exists Boolean indicating whether the phase exists or not
*/
void setPhaseExistence(const int iphase, const bool exists);
protected:
//! Temporary work vector of length m_kk
vector_fp m_grt;
protected:
//! m_kk is the number of species in all of the phases
//! that participate in this kinetics mechanism.
int m_kk;
int m_kk;
//! List of reactions numbers which are reversible reactions
/*!
@ -489,7 +578,7 @@ namespace Cantera {
* production rates and also handles turning thermo
* properties into reaction thermo properties.
*/
ReactionStoichMgr m_rxnstoich;
ReactionStoichMgr m_rxnstoich;
//! Number of irreversible reactions in the mechanism
int m_nirrev;
@ -564,7 +653,7 @@ namespace Cantera {
*/
vector_fp m_mu0;
//! Vector of phase potentials
//! Vector of phase electric potentials
/*!
* Temporary vector containing the potential of each phase
* in the kinetics object
@ -607,30 +696,36 @@ namespace Cantera {
ImplicitSurfChem* m_integrator;
vector_fp m_beta;
//! Vector of reaction indexes specifying the id of the current transfer reactions
//! in the mechanism
/*!
* Vector of reaction indecices which involve current transfers. This provides
* an index into the m_beta array.
*
* irxn = m_ctrxn[i]
*/
vector_int m_ctrxn;
int reactionNumber(){ return m_ii;}
//! Vector of booleans indicating whether the charge transfer reaction may be
//! described by an exchange current density expression
vector_int m_ctrxn_ecdf;
void addElementaryReaction(const ReactionData& r);
void addGlobalReaction(const ReactionData& r);
void installReagents(const ReactionData& r);
vector_fp m_StandardConc;
vector_fp m_deltaG0;
vector_fp m_ProdStanConcReac;
void updateKc();
//! Write values into m_index
/*!
* @param rxnNumber reaction number
* @param type reaction type
* @param loc location ??
*/
void registerReaction(int rxnNumber, int type, int loc) {
m_index[rxnNumber] = std::pair<int, int>(type, loc);
}
void applyButlerVolmerCorrection(doublereal* kf);
//! boolean indicating whether mechanism has been finalized
bool m_finalized;
//! Boolean flag indicating whether any reaction in the mechanism
//! has a coverage dependent forward reaction rate
/*!
* If this is true, then the coverage dependence is multiplied into
* the forward reaction rates constant
*/
bool m_has_coverage_dependence;
//! Boolean flag indicating whether any reaction in the mechanism
@ -643,6 +738,37 @@ namespace Cantera {
*/
bool m_has_electrochem_rxns;
//! Boolean flag indicating whether any reaction in the mechanism
//! is described by an exchange current density expression
/*!
* If this is true, the standard state gibbs free energy of the reaction and
* the product of the reactant standard concentrations must be precalculated
* in order to calculate the rate constant.
*/
bool m_has_exchange_current_density_formulation;
//! Int flag to indicate that some phases in the kinetics mechanism are
//! non-existent.
/*!
* We change the ROP vectors to make sure that non-existent phases are treated
* correctly in the kinetics operator. The value of this is equal to the number
* of phases which don't exist.
*/
int m_phaseExistsCheck;
//! 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
*
* length = number of phases in the object
* By default all phases exist.
*/
std::vector<bool> m_phaseExists;
std::vector<bool *> m_rxnPhaseIsReactant;
std::vector<bool *> m_rxnPhaseIsProduct;
int m_ioFlag;
private:

29
Cantera/src/kinetics/Kinetics.cpp
View file

@ -29,9 +29,16 @@ using namespace std;
namespace Cantera {
Kinetics::Kinetics() : m_ii(0), m_thermo(0),
m_index(-1), m_surfphase(-1), m_rxnphase(-1),
m_mindim(4) {}
Kinetics::Kinetics() :
m_ii(0),
m_nTotalSpecies(0),
m_thermo(0),
m_index(-1),
m_surfphase(-1),
m_rxnphase(-1),
m_mindim(4)
{
}
Kinetics::~Kinetics(){}
@ -42,7 +49,8 @@ namespace Cantera {
* throw an exception.
*/
Kinetics::Kinetics(const Kinetics &right) :
m_ii(0),
m_ii(0),
m_nTotalSpecies(0),
m_thermo(0),
m_index(-1),
m_surfphase(-1),
@ -70,6 +78,7 @@ namespace Cantera {
if (this == &right) return *this;
m_ii = right.m_ii;
m_nTotalSpecies = right.m_nTotalSpecies;
m_perturb = right.m_perturb;
m_reactants = right.m_reactants;
m_products = right.m_products;
@ -216,7 +225,8 @@ namespace Cantera {
return thermo(0);
}
/**
//==============================================================================================
/*
* This function takes as an argument the kineticsSpecies index
* (i.e., the list index in the list of species in the kinetics
* manager) and returns the index of the phase owning the
@ -286,6 +296,15 @@ namespace Cantera {
m_phaseindex[m_thermo.back()->id()] = nPhases();
}
void Kinetics::finalize() {
m_nTotalSpecies = 0;
int np = nPhases();
for (int n = 0; n < np; n++) {
int nsp = m_thermo[n]->nSpecies();
m_nTotalSpecies += nsp;
}
}
//! Private function of the class Kinetics, indicating that a function
//! inherited from the base class hasn't had a definition assigned to it

5
Cantera/src/kinetics/Kinetics.h
View file

@ -811,7 +811,7 @@ namespace Cantera {
* any initialization (allocating arrays, etc.) that must be
* done after the reactions are entered.
*/
virtual void finalize() {}
virtual void finalize();
/**
* Add a single reaction to the mechanism. This routine
@ -904,6 +904,9 @@ namespace Cantera {
//! Number of reactions in the mechanism
int m_ii;
//! Number of species in the species vector for this kinetics operator
int m_nTotalSpecies;
/// Vector of perturbation factors for each reaction's rate of
/// progress vector. It is initialized to one.

14
Cantera/src/kinetics/ReactionData.h
View file

@ -25,7 +25,7 @@ namespace Cantera {
number = 0;
rxn_number = 0;
reversible = true;
rateCoeffType = ARRHENIUS;
rateCoeffType = ARRHENIUS_REACTION_RATECOEFF_TYPE;
falloffType = NONE;
error = 0;
equation = "";
@ -36,7 +36,12 @@ namespace Cantera {
}
virtual ~ReactionData(){}
//! type of the reaction
/*!
* The valid types are listed in the file, reaction_defs.h.
*/
int reactionType;
int number;
int rxn_number;
vector_int reactants;
@ -51,7 +56,14 @@ namespace Cantera {
//! True if the current reaction is reversible. False otherwise
bool reversible;
//! type of the rate coefficient for the forward rate constant
/*!
* The valid types are listed in the file, reaction_defs.h and they
* all end in RATECOEFF_TYPE
*/
int rateCoeffType;
vector_fp rateCoeffParameters;
vector_fp auxRateCoeffParameters;
int falloffType;

1
Cantera/src/kinetics/ReactionStoichMgr.h
View file

@ -214,6 +214,7 @@ namespace Cantera {
* \f[
* R_i = R_i * \prod_k C_k^{o_{k,i}}
* \f]
*
* Here \f$ o_{k,i} \f$ is the reaction order of species k in reaction i.
*/
virtual void multiplyReactants(const doublereal* C, doublereal* R);

393
Cantera/src/kinetics/RxnRates.h
View file

@ -18,94 +18,99 @@
#include "ctexceptions.h"
namespace Cantera {
/**
* A rate coefficient of the form
* \f[
* A T^b \exp (-E/RT)
* \f]
//! Arrhenius reaction rate type depends only on temperature
/**
* A reaction rate coefficient of the following form.
*
* \f[
* k_f = A T^b \exp (-E/RT)
* \f]
*
*/
class Arrhenius {
public:
/// return the rate coefficient type.
static int type(){ return ARRHENIUS_REACTION_RATECOEFF_TYPE ; }
/// Default constructor.
Arrhenius() :
m_logA(-1.0E300),
m_b (0.0),
m_E (0.0),
m_A(0.0) {}
/// Constructor with Arrhenius parameters specified with an array.
Arrhenius(int csize, const doublereal* c) :
m_b (c[1]),
m_E (c[2]),
m_A (c[0])
{
if (m_A <= 0.0) {
m_logA = -1.0E300;
} else {
m_logA = log(m_A);
}
}
/// Constructor.
/// @param A pre-exponential. The unit system is
/// (kmol, m, s). The actual units depend on the reaction
/// order and the dimensionality (surface or bulk).
/// @param b Temperature exponent. Non-dimensional.
/// @param E Activation energy in temperature units. Kelvin.
Arrhenius(doublereal A, doublereal b, doublereal E) :
m_b (b),
m_E (E),
m_A (A)
{
if (m_A <= 0.0) {
m_logA = -1.0E300;
} else {
m_logA = log(m_A);
}
}
//! Update concentration-dependent parts of the rate coefficient.
/*!
* For this class, there are no
* concentration-dependent parts, so this method does nothing.
*/
class Arrhenius {
public:
void update_C(const doublereal* c) {
}
/// return the rate coefficient type.
static int type(){ return ARRHENIUS; }
/// Default constructor.
Arrhenius() :
m_logA(-1.0E300),
m_b (0.0),
m_E (0.0),
m_A(0.0) {}
/// Constructor with Arrhenius parameters specified with an array.
Arrhenius(int csize, const doublereal* c) :
m_b (c[1]),
m_E (c[2]),
m_A (c[0])
{
if (m_A <= 0.0) {
m_logA = -1.0E300;
} else {
m_logA = log(m_A);
}
}
/// Constructor.
/// @param A pre-exponential. The unit system is
/// (kmol, m, s). The actual units depend on the reaction
/// order and the dimensionality (surface or bulk).
/// @param b Temperature exponent. Non-dimensional.
/// @param E Activation energy in temperature units. Kelvin.
Arrhenius(doublereal A, doublereal b, doublereal E) :
m_b (b),
m_E (E),
m_A (A)
{
if (m_A <= 0.0) {
m_logA = -1.0E300;
} else {
m_logA = log(m_A);
}
}
/// Update concentration-dependent parts of the rate
/// coefficient. For this class, there are no
/// concentration-dependent parts, so this method does
/// nothing.
void update_C(const doublereal* c) {}
/**
* Update the value of the logarithm of the rate constant.
*
* Note, this function should never be called for negative A values.
* If it does then it will produce a negative overflow result, and
* a zero net forwards reaction rate, instead of a negative reaction
* rate constant that is the expected result.
*/
doublereal update(doublereal logT, doublereal recipT) const {
return m_logA + m_b*logT - m_E*recipT;
}
/**
* Update the value of the logarithm of the rate constant.
*
* Note, this function should never be called for negative A values.
* If it does then it will produce a negative overflow result, and
* a zero net forwards reaction rate, instead of a negative reaction
* rate constant that is the expected result.
*/
doublereal update(doublereal logT, doublereal recipT) const {
return m_logA + m_b*logT - m_E*recipT;
}
/**
* Update the value the rate constant.
*
* This function returns the actual value of the rate constant.
* It can be safely called for negative values of the pre-exponential
* factor.
*/
doublereal updateRC(doublereal logT, doublereal recipT) const {
return m_A * exp(m_b*logT - m_E*recipT);
}
/**
* Update the value the rate constant.
*
* This function returns the actual value of the rate constant.
* It can be safely called for negative values of the pre-exponential
* factor.
*/
doublereal updateRC(doublereal logT, doublereal recipT) const {
return m_A * exp(m_b*logT - m_E*recipT);
}
void writeUpdateRHS(std::ostream& s) const {
s << " exp(" << m_logA;
if (m_b != 0.0) s << " + " << m_b << " * tlog";
if (m_E != 0.0) s << " - " << m_E << " * rt";
s << ");" << std::endl;
}
void writeUpdateRHS(std::ostream& s) const {
s << " exp(" << m_logA;
if (m_b != 0.0) s << " + " << m_b << " * tlog";
if (m_E != 0.0) s << " - " << m_E << " * rt";
s << ");" << std::endl;
}
doublereal activationEnergy_R() const {
return m_E;
@ -121,19 +126,19 @@ namespace Cantera {
class ArrheniusSum {
public:
static int type(){ return ARRHENIUS_SUM; }
static int type(){ return ARRHENIUS_SUM_REACTION_RATECOEFF_TYPE; }
ArrheniusSum() : m_nterms(0) {}
void addArrheniusTerm(doublereal A, doublereal b, doublereal E) {
if (A > 0.0) {
m_terms.push_back(Arrhenius(A, b, E));
m_sign.push_back(1);
}
else if (A < 0.0) {
m_terms.push_back(Arrhenius(-A, b, E));
m_sign.push_back(-1);
}
m_nterms++;
if (A > 0.0) {
m_terms.push_back(Arrhenius(A, b, E));
m_sign.push_back(1);
}
else if (A < 0.0) {
m_terms.push_back(Arrhenius(-A, b, E));
m_sign.push_back(-1);
}
m_nterms++;
}
void update_C(const doublereal* c) {}
@ -169,38 +174,38 @@ namespace Cantera {
return fsum;
}
void writeUpdateRHS(std::ostream& s) const {
void writeUpdateRHS(std::ostream& s) const {
;
}
static bool alwaysComputeRate() { return false;}
protected:
std::vector<Arrhenius> m_terms;
vector_int m_sign;
int m_nterms;
std::vector<Arrhenius> m_terms;
vector_int m_sign;
int m_nterms;
};
/**
* An Arrhenius rate with coverage-dependent terms.
*/
class SurfaceArrhenius {
/**
* An Arrhenius rate with coverage-dependent terms.
*/
class SurfaceArrhenius {
public:
static int type(){ return ARRHENIUS; }
SurfaceArrhenius() :
m_logA(-1.0E300),
m_b (0.0),
m_E (0.0),
m_A(0.0),
m_acov(0.0),
m_ecov(0.0),
m_mcov(0.0),
m_ncov(0),
m_nmcov(0)
{
}
public:
static int type(){ return ARRHENIUS_REACTION_RATECOEFF_TYPE ; }
SurfaceArrhenius() :
m_logA(-1.0E300),
m_b (0.0),
m_E (0.0),
m_A(0.0),
m_acov(0.0),
m_ecov(0.0),
m_mcov(0.0),
m_ncov(0),
m_nmcov(0)
{
}
SurfaceArrhenius( int csize, const doublereal* c ) :
m_b (c[1]),
@ -325,7 +330,7 @@ namespace Cantera {
return exp(lres);
}
void writeUpdateRHS(std::ostream& s) const {}
void writeUpdateRHS(std::ostream& s) const {}
protected:
doublereal delta_s0, delta_e0;
@ -335,37 +340,143 @@ namespace Cantera {
#endif
}
// class LandauTeller {
// public:
// static int type(){ return LANDAUTELLER; }
// LandauTeller(){}
// LandauTeller( const vector_fp& c ) : m_c(c) { m_c[0] = log(c[0]); }
//! Arrhenius reaction rate type depends only on temperature
/**
* A reaction rate coefficient of the following form.
*
* \f[
* k_f = A T^b \exp (-E/RT)
* \f]
*
*/
class ExchangeCurrent {
public:
// doublereal update(doublereal logT, doublereal recipT) const {
// return m_c[0] + m_c[1]*tt[1] - m_c[2]*tt[2]
// + m_c[3]*tt[3] + m_c[4]*tt[4];
// }
// //void writeUpdateRHS(ostream& s) const {
// // s << exp(m_logA);
// // s << " * exp(";
// // if (m_b != 0.0) s << m_b << " * tlog";
// // if (m_E != 0.0) s << " - " << m_E << " * rt";
// // if (m_E != 0.0) s << " - " << m_E << " * rt";
// // s << ");" << endl;
// // }
// //}
//! return the rate coefficient type.
static int type() {
return EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE ;
}
// protected:
// doublereal m_logA, m_b, m_E;
// };
//! Default constructor.
ExchangeCurrent() :
m_logA(-1.0E300),
m_b (0.0),
m_E (0.0),
m_A(0.0) {}
//! Constructor with Arrhenius parameters specified with an array.
ExchangeCurrent(int csize, const doublereal* c) :
m_b (c[1]),
m_E (c[2]),
m_A (c[0])
{
if (m_A <= 0.0) {
m_logA = -1.0E300;
} else {
m_logA = log(m_A);
}
}
/// Constructor.
/// @param A pre-exponential. The unit system is
/// (kmol, m, s). The actual units depend on the reaction
/// order and the dimensionality (surface or bulk).
/// @param b Temperature exponent. Non-dimensional.
/// @param E Activation energy in temperature units. Kelvin.
ExchangeCurrent(doublereal A, doublereal b, doublereal E) :
m_b (b),
m_E (E),
m_A (A)
{
if (m_A <= 0.0) {
m_logA = -1.0E300;
} else {
m_logA = log(m_A);
}
}
//! Update concentration-dependent parts of the rate coefficient.
/*!
* For this class, there are no
* concentration-dependent parts, so this method does nothing.
*/
void update_C(const doublereal* c) {
}
/**
* Update the value of the logarithm of the rate constant.
*
* Note, this function should never be called for negative A values.
* If it does then it will produce a negative overflow result, and
* a zero net forwards reaction rate, instead of a negative reaction
* rate constant that is the expected result.
*/
doublereal update(doublereal logT, doublereal recipT) const {
return m_logA + m_b*logT - m_E*recipT;
}
/**
* Update the value the rate constant.
*
* This function returns the actual value of the rate constant.
* It can be safely called for negative values of the pre-exponential
* factor.
*/
doublereal updateRC(doublereal logT, doublereal recipT) const {
return m_A * exp(m_b*logT - m_E*recipT);
}
void writeUpdateRHS(std::ostream& s) const {
s << " exp(" << m_logA;
if (m_b != 0.0) s << " + " << m_b << " * tlog";
if (m_E != 0.0) s << " - " << m_E << " * rt";
s << ");" << std::endl;
}
doublereal activationEnergy_R() const {
return m_E;
}
static bool alwaysComputeRate() { return false;}
protected:
doublereal m_logA, m_b, m_E, m_A;
};
// class LandauTeller {
// public:
// static int type(){ return LANDAUTELLER; }
// LandauTeller(){}
// LandauTeller( const vector_fp& c ) : m_c(c) { m_c[0] = log(c[0]); }
// doublereal update(doublereal logT, doublereal recipT) const {
// return m_c[0] + m_c[1]*tt[1] - m_c[2]*tt[2]
// + m_c[3]*tt[3] + m_c[4]*tt[4];
// }
// //void writeUpdateRHS(ostream& s) const {
// // s << exp(m_logA);
// // s << " * exp(";
// // if (m_b != 0.0) s << m_b << " * tlog";
// // if (m_E != 0.0) s << " - " << m_E << " * rt";
// // if (m_E != 0.0) s << " - " << m_E << " * rt";
// // s << ");" << endl;
// // }
// //}
// protected:
// doublereal m_logA, m_b, m_E;
// };
//}
//}
}
#endif

24
Cantera/src/kinetics/importKinetics.cpp
View file

@ -521,9 +521,24 @@ namespace Cantera {
* kf should point to a XML element named "rateCoeff".
* rdata is the partially filled ReactionData object for the reaction.
* This function will fill in more fields in the ReactionData object.
*
* @param kf Reference to the XML Node named rateCoeff
*/
void getRateCoefficient(const node_t& kf, kinetics_t& kin,
ReactionData& rdata, int negA) {
string type = kf.attrib("type");
if (type == "") {
type = "Arrhenius";
rdata.rateCoeffType = ARRHENIUS_REACTION_RATECOEFF_TYPE;
}
if (type == "ExchangeCurrentDensity") {
rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE;
} else if (type == "Arrhenius") {
} else {
throw CanteraError("getRateCoefficient",
"Unknown type: " + type);
}
int nc = kf.nChildren();
nodeset_t& kf_children = kf.children();
@ -557,6 +572,12 @@ namespace Cantera {
"negative or zero A coefficient for reaction "+int2str(rdata.number));
}
}
else if (nm == "Arrhenius_ExchangeCurrentDensity") {
vector_fp coeff(3);
getArrhenius(c, highlow, coeff[0], coeff[1], coeff[2]);
chigh = coeff;
rdata.rateCoeffType = EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE;
}
else if (nm == "falloff") {
getFalloff(c, rdata);
}
@ -580,6 +601,9 @@ namespace Cantera {
rdata.auxRateCoeffParameters = clow;
else if (rdata.reactionType == CHEMACT_RXN)
rdata.auxRateCoeffParameters = chigh;
}

11
Cantera/src/kinetics/reaction_defs.h
View file

@ -81,11 +81,12 @@ namespace Cantera {
*/
//@{
const int ARRHENIUS = 1;
const int LANDAUTELLER = 2;
const int TSTRATE = 3;
const int SURF_ARRHENIUS = 4;
const int ARRHENIUS_SUM = 5;
const int ARRHENIUS_REACTION_RATECOEFF_TYPE = 1;
const int LANDAUTELLER_REACTION_RATECOEFF_TYPE = 2;
const int TSTRATE_REACTION_RATECOEFF_TYPE = 3;
const int SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE = 4;
const int ARRHENIUS_SUM_REACTION_RATECOEFF_TYPE = 5;
const int EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE = 6;
//@}