InterfaceKinetics rewrite -> addition of BV and Affinity reactions, next interation.
This commit is contained in:
Harry Moffat
parent
6a07722675
commit
021c8f9880
6 changed files with 717 additions and 46 deletions
|
|
@ -45,7 +45,7 @@ public:
|
|||
return cEdgeKinetics;
|
||||
}
|
||||
|
||||
// defined in InterfaceKinetics.cpp
|
||||
|
||||
virtual void finalize();
|
||||
};
|
||||
}
|
||||
|
|
|
|||
|
|
@ -74,26 +74,50 @@ public:
|
|||
*/
|
||||
virtual void updateROP();
|
||||
|
||||
virtual void determineFwdOrdersBV(ReactionData& rdata, std::vector<doublereal>& fwdFullorders);
|
||||
|
||||
//void addGlobalReaction(ReactionData& r);
|
||||
|
||||
double calcForwardROP_BV(size_t irxn, size_t iBeta);
|
||||
double calcForwardROP_BV(size_t irxn, size_t iBeta, double ioc, double nStoich, double nu, doublereal ioNet);
|
||||
|
||||
double calcForwardROP_BV_NoAct(size_t irxn, size_t iBeta, double ioc, double nStoich, double nu, doublereal ioNet);
|
||||
|
||||
bool getExchangeCurrentDensityFormulation(size_t irxn, doublereal& nStoich, doublereal& OCV, doublereal& io,
|
||||
doublereal& overPotential, doublereal& beta, doublereal& resistance);
|
||||
|
||||
//! Calculate the open circuit voltage of a given reaction
|
||||
/*!
|
||||
* If the reaction has no electron transport, then return 0.0
|
||||
*
|
||||
* @param irxn Reaction id
|
||||
*/
|
||||
double openCircuitVoltage(size_t irxn);
|
||||
|
||||
double calcCurrentDensity(double nu, double nStoich, double io, double beta, double temp, doublereal resistivity = 0.0) const;
|
||||
|
||||
double solveCurrentRes(doublereal nu, doublereal nStoich, doublereal ioc, doublereal beta, doublereal temp,
|
||||
doublereal resistivity = 0.0, int iprob = 0) const;
|
||||
|
||||
//! Prepare the class for the addition of reactions
|
||||
/*!
|
||||
* (virtual from Kinetics)
|
||||
* We determine the metal phase and solution phase here
|
||||
*/
|
||||
virtual void init();
|
||||
|
||||
|
||||
protected:
|
||||
//! Index of the metal phase in the list of phases for this kinetics object
|
||||
//! Index of the metal phase in the list of phases for this kinetics object. This is the electron phase.
|
||||
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
|
||||
size_t solnPhaseRS_;
|
||||
|
||||
//! Index of the electrons species in the list of species for this surface kinetics, if none set it to -1
|
||||
size_t kElectronRS_;
|
||||
|
||||
|
||||
|
||||
};
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -8,8 +8,11 @@
|
|||
#ifndef CT_IFACEKINETICS_H
|
||||
#define CT_IFACEKINETICS_H
|
||||
|
||||
#include "cantera/thermo/mix_defs.h"
|
||||
#include "Kinetics.h"
|
||||
#include "cantera/kinetics/RxnMolChange.h"
|
||||
|
||||
#include "cantera/base/utilities.h"
|
||||
#include "RateCoeffMgr.h"
|
||||
#include "ReactionStoichMgr.h"
|
||||
|
||||
|
|
@ -445,7 +448,7 @@ public:
|
|||
int phaseStability(const size_t iphase) const;
|
||||
|
||||
|
||||
void determineFwdOrdersBV(ReactionData& rdata, std::vector<doublereal>& fwdFullorders);
|
||||
virtual void determineFwdOrdersBV(ReactionData& rdata, std::vector<doublereal>& fwdFullorders);
|
||||
|
||||
protected:
|
||||
//! Temporary work vector of length m_kk
|
||||
|
|
@ -713,6 +716,8 @@ protected:
|
|||
*/
|
||||
std::vector<RxnOrders*> m_ctrxn_FwdOrdersList_;
|
||||
|
||||
std::vector<doublereal> m_ctrxn_resistivity_;
|
||||
|
||||
//! Vector of standard concentrations
|
||||
/*!
|
||||
* Length number of kinetic species
|
||||
|
|
|
|||
|
|
@ -434,6 +434,10 @@ public:
|
|||
* units in array kc, which must be dimensioned at least as large as the
|
||||
* total number of reactions.
|
||||
*
|
||||
* \f[
|
||||
* Kc_i = exp [ \Delta G_{ss,i} ] prod(Cs_k) exp(\sum_k \nu_{k,i} F \phi_n) ]
|
||||
* \f]
|
||||
*
|
||||
* @param kc Output vector containing the equilibrium constants.
|
||||
* Length: m_ii.
|
||||
*/
|
||||
|
|
@ -446,7 +450,7 @@ public:
|
|||
* property values \f$ z_k, k = 1, \dots, K \f$, return the
|
||||
* array of reaction values
|
||||
* \f[
|
||||
* \Delta Z_i = \sum_k \nu_{k,i} z_k, i = 1, \dots, I.
|
||||
* \Delta Z_i = \sum_k \nu_{k,i} z_k, i = 1, \dots, I.
|
||||
* \f]
|
||||
* For example, if this method is called with the array of
|
||||
* standard-state molar Gibbs free energies for the species,
|
||||
|
|
|
|||
|
|
@ -6,6 +6,8 @@
|
|||
#include "cantera/thermo/SurfPhase.h"
|
||||
#include "cantera/base/utilities.h"
|
||||
|
||||
#include <cstdio>
|
||||
|
||||
using namespace std;
|
||||
|
||||
namespace Cantera
|
||||
|
|
@ -14,7 +16,6 @@ namespace Cantera
|
|||
ElectrodeKinetics::ElectrodeKinetics(thermo_t* thermo) :
|
||||
InterfaceKinetics(thermo),
|
||||
metalPhaseRS_(npos),
|
||||
electronPhaseRS_(npos),
|
||||
solnPhaseRS_(npos),
|
||||
kElectronRS_(npos)
|
||||
{
|
||||
|
|
@ -48,7 +49,6 @@ ElectrodeKinetics& ElectrodeKinetics::operator=(const ElectrodeKinetics& right)
|
|||
InterfaceKinetics::operator=(right);
|
||||
|
||||
metalPhaseRS_ = right.metalPhaseRS_;
|
||||
electronPhaseRS_ = right.electronPhaseRS_;
|
||||
solnPhaseRS_ = right.solnPhaseRS_;
|
||||
kElectronRS_ = right.kElectronRS_;
|
||||
|
||||
|
|
@ -67,12 +67,12 @@ Kinetics* ElectrodeKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> &
|
|||
return iK;
|
||||
}
|
||||
//============================================================================================================================
|
||||
//====================================================================================================================
|
||||
// Identify the metal phase and the electrons species
|
||||
// Identify the metal phase and the electron species
|
||||
void ElectrodeKinetics::identifyMetalPhase()
|
||||
{
|
||||
metalPhaseRS_ = npos;
|
||||
kElectronRS_ = -1;
|
||||
kElectronRS_ = npos;
|
||||
solnPhaseRS_ = npos;
|
||||
size_t np = nPhases();
|
||||
//
|
||||
// Identify the metal phase as the phase with the electron species (element index of 1 for element E
|
||||
|
|
@ -105,6 +105,8 @@ void ElectrodeKinetics::identifyMetalPhase()
|
|||
// Identify the solution phase as a 3D phase, with nonzero phase charge change
|
||||
// in at least one reaction
|
||||
//
|
||||
/*
|
||||
* Haven't filled in reactions yet when this is called, unlike previous treatment.
|
||||
if (iph != metalPhaseRS_) {
|
||||
for (size_t i = 0; i < m_ii; i++) {
|
||||
RxnMolChange* rmc = rmcVector[i];
|
||||
|
|
@ -116,8 +118,38 @@ void ElectrodeKinetics::identifyMetalPhase()
|
|||
}
|
||||
}
|
||||
}
|
||||
*/
|
||||
//
|
||||
// New method is to find the first multispecies 3D phase with charged species as the solution phase
|
||||
//
|
||||
if (iph != metalPhaseRS_) {
|
||||
ThermoPhase& tp =*( m_thermo[iph]);
|
||||
size_t nsp = tp.nSpecies();
|
||||
size_t nd = tp.nDim();
|
||||
if (nd == 3 && nsp > 1) {
|
||||
for (size_t k = 0; k < nsp; k++) {
|
||||
if (tp.charge(k) != 0.0) {
|
||||
solnPhaseRS_ = iph;
|
||||
string ss = tp.name();
|
||||
// cout << "solution phase = "<< ss << endl;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
//
|
||||
// Right now, if we don't find an electron phase, let's error exit to avoid complications
|
||||
//
|
||||
if (metalPhaseRS_ == npos) {
|
||||
throw CanteraError("ElectrodeKinetics::identifyMetalPhase()",
|
||||
"Can't find electron phase -> treating this as an error right now");
|
||||
}
|
||||
if (solnPhaseRS_ == npos) {
|
||||
throw CanteraError("ElectrodeKinetics::identifyMetalPhase()",
|
||||
"Can't find solution phase -> treating this as an error right now");
|
||||
}
|
||||
}
|
||||
//============================================================================================================================
|
||||
// virtual from InterfaceKinetics
|
||||
|
|
@ -148,7 +180,6 @@ void ElectrodeKinetics::updateROP()
|
|||
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
|
||||
|
|
@ -182,7 +213,11 @@ void ElectrodeKinetics::updateROP()
|
|||
// the rxnstoich calculation contained the dependence of the current density on the activity concentrations
|
||||
// We finish up with the ROP calculation
|
||||
//
|
||||
|
||||
int iECDFormulation = m_ctrxn_ecdf[iBeta];
|
||||
if (iECDFormulation == 0) {
|
||||
throw CanteraError(" ElectrodeKinetics::updateROP()",
|
||||
"Straight kfwrd with BUTLERVOLMER_RXN not handled yet");
|
||||
}
|
||||
//
|
||||
// Get the phase mole change structure
|
||||
//
|
||||
|
|
@ -209,22 +244,151 @@ void ElectrodeKinetics::updateROP()
|
|||
// Calculate the overpotential
|
||||
//
|
||||
double nu = voltage - OCV;
|
||||
|
||||
//
|
||||
// Find the product of the standard concentrations for ROP orders that we used above
|
||||
//
|
||||
const RxnOrders* ro_rop = m_ctrxn_ROPOrdersList_[iBeta];
|
||||
if (ro_rop == 0) {
|
||||
throw CanteraError("ElectrodeKinetics::", "ROP orders pointer is zero ?!?");
|
||||
}
|
||||
double tmp2 = 1.0;
|
||||
const std::vector<size_t>& kinSpeciesIDs = ro_rop->kinSpeciesIDs_;
|
||||
const std::vector<doublereal>& kinSpeciesOrders = ro_rop->kinSpeciesOrders_;
|
||||
for (size_t j = 0; j < kinSpeciesIDs.size(); j++) {
|
||||
size_t k = kinSpeciesIDs[j];
|
||||
double oo = kinSpeciesOrders[j];
|
||||
tmp2 *= pow(m_StandardConc[k], oo);
|
||||
}
|
||||
//
|
||||
// Now have to divide this to get rid of standard concentrations. We should
|
||||
// have used just the activities in the m_rxnstoich.multiplyReactants(DATA_PTR(m_actConc), DATA_PTR(m_ropf));
|
||||
// calculation above!
|
||||
// That is because the exchange current density rate constants have the correct units in the first place.
|
||||
//
|
||||
m_ropf[irxn] /= tmp2;
|
||||
//
|
||||
// Calculate the exchange current density
|
||||
// m_ropf contains the exchange current reaction rate
|
||||
//
|
||||
double io = m_ropf[irxn] * nStoichElectrons;
|
||||
|
||||
double ioc = m_ropf[irxn] * nStoichElectrons;
|
||||
//
|
||||
// Add in the film resistance here
|
||||
//
|
||||
double resist = m_ctrxn_resistivity_[iBeta];
|
||||
double exp1 = nu * nStoichElectrons * beta / rtdf;
|
||||
double exp2 = - nu * nStoichElectrons * (1.0 - beta) / (rtdf);
|
||||
m_ropnet[irxn] = io * (exp(exp1) - exp(exp2));
|
||||
double io = ioc * (exp(exp1) - exp(exp2));
|
||||
|
||||
// Need to resurrect the forwards rate constant.
|
||||
//m_ropf[irxn] = ;
|
||||
m_ropr[irxn] = m_ropnet[irxn] - m_ropf[irxn];
|
||||
if (resist != 0.0) {
|
||||
io = solveCurrentRes(nu, nStoichElectrons, ioc, beta, TT, resist, 0);
|
||||
}
|
||||
|
||||
m_ropnet[irxn] = io / (Faraday * nStoichElectrons);
|
||||
//
|
||||
// Need to resurrect the forwards rate of progress -> there is some need to
|
||||
// calculate each direction individually
|
||||
//
|
||||
m_ropf[irxn] = calcForwardROP_BV(irxn, iBeta, ioc, nStoichElectrons, nu, io);
|
||||
//
|
||||
// Calculate the reverse rate of progress from the difference
|
||||
//
|
||||
m_ropr[irxn] = m_ropf[irxn] - m_ropnet[irxn];
|
||||
|
||||
}
|
||||
} else if (reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_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
|
||||
//
|
||||
int iECDFormulation = m_ctrxn_ecdf[iBeta];
|
||||
if (iECDFormulation == 0) {
|
||||
throw CanteraError("ElectrodeKinetics::updateROP()",
|
||||
"Straight kfwrd with BUTLERVOLMER_NOACTIVITYCOEFFS_RXN not handled yet");
|
||||
}
|
||||
//
|
||||
// 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;
|
||||
//
|
||||
// Unfortunately, we really need to recalculate everything from almost scratch
|
||||
// for this case, since it widely diverges from the thermo norm.
|
||||
//
|
||||
// Start with the exchange current reaction rate constant, which should
|
||||
// be located in m_rfn[].
|
||||
//
|
||||
double ioc = m_rfn[irxn] * nStoichElectrons * m_perturb[irxn];
|
||||
//
|
||||
// Now we need th mole fraction vector and we need the RxnOrders vector.
|
||||
//
|
||||
const RxnOrders* ro_fwd = m_ctrxn_ROPOrdersList_[iBeta];
|
||||
if (ro_fwd == 0) {
|
||||
throw CanteraError("ElectrodeKinetics::calcForwardROP_BV()", "forward orders pointer is zero ?!?");
|
||||
}
|
||||
double tmp = 1.0;
|
||||
double mfS = 0.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 ks = kinSpeciesIDs[j];
|
||||
thermo_t& th = speciesPhase(ks);
|
||||
size_t n = speciesPhaseIndex(ks);
|
||||
size_t klocal = ks - m_start[n];
|
||||
mfS = th.moleFraction(klocal);
|
||||
|
||||
double oo = kinSpeciesOrders[j];
|
||||
tmp *= pow(mfS, oo);
|
||||
}
|
||||
ioc *= tmp;
|
||||
//
|
||||
// Add in the film resistance here, later
|
||||
//
|
||||
double resist = m_ctrxn_resistivity_[iBeta];
|
||||
double exp1 = nu * nStoichElectrons * beta / rtdf;
|
||||
double exp2 = - nu * nStoichElectrons * (1.0 - beta) / (rtdf);
|
||||
double io = ioc * (exp(exp1) - exp(exp2));
|
||||
if (resist != 0.0) {
|
||||
io = solveCurrentRes(nu, nStoichElectrons, ioc, beta, TT, resist, 0);
|
||||
}
|
||||
|
||||
m_ropnet[irxn] = io / (Faraday * nStoichElectrons);
|
||||
//
|
||||
// Need to resurrect the forwards rate of progress -> there is some need to
|
||||
// calculate each direction individually
|
||||
//
|
||||
m_ropf[irxn] = calcForwardROP_BV_NoAct(irxn, iBeta, ioc, nStoichElectrons, nu, io);
|
||||
//
|
||||
// Calculate the reverse rate of progress from the difference
|
||||
//
|
||||
m_ropr[irxn] = m_ropf[irxn] - m_ropnet[irxn];
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
|
@ -299,13 +463,158 @@ void ElectrodeKinetics::updateROP()
|
|||
}
|
||||
//==================================================================================================================
|
||||
//
|
||||
// This version of takes the electrons out of the reaction rate expression
|
||||
// (note: with proper specification of the phase, this shouldn't make a numerical difference (power of 1).
|
||||
// But it certainly is a complication and unneeded work)
|
||||
// (TODO: probably can take stoichiometric solids out of the reaction order expression as well.
|
||||
// They all contribute powers of 1 as well)
|
||||
//
|
||||
void ElectrodeKinetics::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[j];
|
||||
if (kkin != kElectronRS_) {
|
||||
fwdFullorders[kkin] += betaf * oo;
|
||||
if (abs(fwdFullorders[kkin]) < 0.00001) {
|
||||
fwdFullorders[kkin] = 0.0;
|
||||
}
|
||||
} else {
|
||||
fwdFullorders[kkin] = 0.0;
|
||||
}
|
||||
}
|
||||
for (size_t j = 0; j < rdata.products.size(); j++) {
|
||||
size_t kkin = rdata.products[j];
|
||||
double oo = rdata.pstoich[j];
|
||||
if (kkin != kElectronRS_) {
|
||||
fwdFullorders[kkin] -= betaf * oo;
|
||||
if (abs(fwdFullorders[kkin]) < 0.00001) {
|
||||
fwdFullorders[kkin] = 0.0;
|
||||
}
|
||||
} else {
|
||||
fwdFullorders[kkin] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
//==================================================================================================================
|
||||
//
|
||||
// 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)
|
||||
{
|
||||
double ElectrodeKinetics::calcForwardROP_BV(size_t irxn, size_t iBeta, double ioc, double nStoich, double nu, doublereal ioNet)
|
||||
{
|
||||
double ropf;
|
||||
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];
|
||||
doublereal beta = m_beta[iBeta];
|
||||
|
||||
#ifdef DEBUG_MODE
|
||||
//
|
||||
// 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
|
||||
//
|
||||
const RxnOrders* ro_fwd = m_ctrxn_FwdOrdersList_[iBeta];
|
||||
if (ro_fwd == 0) {
|
||||
throw CanteraError("ElectrodeKinetics::calcForwardROP_BV()", "forward orders pointer is zero ?!?");
|
||||
}
|
||||
double tmp = exp(- m_beta[iBeta] * m_deltaG0[irxn] * rrt);
|
||||
double tmp2 = 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];
|
||||
tmp2 *= pow(m_StandardConc[k], oo);
|
||||
}
|
||||
|
||||
//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
|
||||
//
|
||||
tmp = 1.0;
|
||||
|
||||
for (size_t j = 0; j < kinSpeciesIDs.size(); j++) {
|
||||
size_t k = kinSpeciesIDs[j];
|
||||
double oo = kinSpeciesOrders[j];
|
||||
tmp *= pow(m_actConc[k], oo);
|
||||
}
|
||||
ropf = kf * tmp;
|
||||
#endif
|
||||
//
|
||||
// Now calculate ropf in a separate but equivalent way.
|
||||
// totally equivalent way if resistivity is zero, should be equal (HKM -> Proved exactly in one case)
|
||||
//
|
||||
double iof = ioc;
|
||||
double resistivity = m_ctrxn_resistivity_[iBeta];
|
||||
if (fabs(resistivity * ioNet) > fabs(nu)) {
|
||||
ioNet = nu / resistivity;
|
||||
}
|
||||
if (nStoich > 0.0) {
|
||||
double exp1 = nStoich * Faraday * beta * (nu - resistivity * ioNet)/ (rt);
|
||||
iof *= exp(exp1);
|
||||
} else {
|
||||
#ifdef DEBUG_MODE
|
||||
if (ioc > 0) {
|
||||
throw CanteraError(" ", "ioc should be less than zero here");
|
||||
}
|
||||
#endif
|
||||
double exp2 = -nu * nStoich * Faraday * (1.0 - beta) / (rt);
|
||||
iof = ioc * ( - exp(exp2));
|
||||
}
|
||||
ropf = iof / ( Faraday * nStoich);
|
||||
|
||||
return ropf;
|
||||
}
|
||||
//==================================================================================================================
|
||||
//
|
||||
// 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_NoAct(size_t irxn, size_t iBeta, double ioc, double nStoich, double nu,
|
||||
doublereal ioNet)
|
||||
{
|
||||
doublereal TT = thermo(0).temperature();
|
||||
doublereal rt = GasConstant * TT;
|
||||
//doublereal rrt = 1.0/rt;
|
||||
doublereal beta = m_beta[iBeta];
|
||||
|
||||
/*
|
||||
//
|
||||
// Calculate gather the exchange current reaction rate constant (where does n_s appear?)
|
||||
//
|
||||
|
|
@ -320,9 +629,10 @@ double ElectrodeKinetics::calcForwardROP_BV(size_t irxn, size_t iBeta)
|
|||
}
|
||||
//
|
||||
// Calculate the forward chemical and modify the forward reaction rate coefficient
|
||||
// (we don't use standard concentrations at all here);
|
||||
//
|
||||
double tmp = exp(- m_beta[iBeta] * m_deltaG0[irxn] * rrt);
|
||||
double tmp2 = m_ProdStanConcReac[irxn];
|
||||
double tmp2 = 1.0;
|
||||
tmp *= 1.0 / tmp2 / Faraday;
|
||||
//
|
||||
// Calculate the chemical reaction rate constant
|
||||
|
|
@ -345,16 +655,273 @@ double ElectrodeKinetics::calcForwardROP_BV(size_t irxn, size_t iBeta)
|
|||
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];
|
||||
|
||||
size_t ks = kinSpeciesIDs[j];
|
||||
thermo_t& th = speciesPhase(ks);
|
||||
size_t n = speciesPhaseIndex(ks);
|
||||
size_t klocal = ks - m_start[n];
|
||||
double mfS = th.moleFraction(klocal);
|
||||
double oo = kinSpeciesOrders[j];
|
||||
tmp *= pow(m_actConc[k], oo);
|
||||
tmp *= pow(mfS, oo);
|
||||
}
|
||||
double ropf = kf * tmp;
|
||||
|
||||
*/
|
||||
/*
|
||||
if (nStoich > 0) {
|
||||
double ropf = ioc / ( Faraday * nStoich);
|
||||
double exp1 = nu * nStoich * Faraday * beta / (rt);
|
||||
ropf *= exp(exp1);
|
||||
} else {
|
||||
double ropf = ioc / ( Faraday * nStoich);
|
||||
double exp1 = nu * nStoich * Faraday * beta / (rt);
|
||||
ropf *= exp(exp1);
|
||||
}
|
||||
*/
|
||||
//
|
||||
// With all of the thermo issues, I'm thinking this is the best we can do
|
||||
// (it certainly maintains the forward and reverse rates of progress as being positive)
|
||||
//
|
||||
double iof = ioc;
|
||||
double resistivity = m_ctrxn_resistivity_[iBeta];
|
||||
if (fabs(resistivity * ioNet) > fabs(nu)) {
|
||||
ioNet = nu / resistivity;
|
||||
}
|
||||
if (nStoich > 0) {
|
||||
double exp1 = nStoich * Faraday * beta * (nu - resistivity * ioNet)/ (rt);
|
||||
iof *= exp(exp1);
|
||||
} else {
|
||||
#ifdef DEBUG_MODE
|
||||
if (ioc > 0) {
|
||||
throw CanteraError(" ", "ioc should be less than zero here");
|
||||
}
|
||||
#endif
|
||||
double exp2 = -nu * nStoich * Faraday * (1.0 - beta) / (rt);
|
||||
iof = ioc * ( - exp(exp2));
|
||||
}
|
||||
double ropf = iof / ( Faraday * nStoich);
|
||||
return ropf;
|
||||
}
|
||||
//==================================================================================================================
|
||||
|
||||
|
||||
|
||||
double ElectrodeKinetics::openCircuitVoltage(size_t irxn)
|
||||
{
|
||||
//
|
||||
// Calculate deltaG for all reactions
|
||||
//
|
||||
getDeltaGibbs(0);
|
||||
//
|
||||
// Look up the net number of electrons that are products.
|
||||
//
|
||||
RxnMolChange* rmc = rmcVector[irxn];
|
||||
double nStoichElectrons = - rmc->m_phaseChargeChange[metalPhaseRS_];
|
||||
double OCV = 0.0;
|
||||
if (nStoichElectrons != 0.0) {
|
||||
OCV = m_deltaG[irxn] / Faraday / nStoichElectrons;
|
||||
}
|
||||
return OCV;
|
||||
}
|
||||
//==================================================================================================================
|
||||
//
|
||||
// Returns the local exchange current density formulation parameters
|
||||
//
|
||||
bool ElectrodeKinetics::
|
||||
getExchangeCurrentDensityFormulation(size_t irxn,
|
||||
doublereal& nStoichElectrons, doublereal& OCV, doublereal& io,
|
||||
doublereal& overPotential, doublereal& beta,
|
||||
doublereal& resistivity)
|
||||
{
|
||||
size_t iBeta = npos;
|
||||
beta = 0.0;
|
||||
//
|
||||
// Add logic to handle other reaction types -> return 0 if formulation isn't compatible
|
||||
//
|
||||
|
||||
// evaluate rate constants and equilibrium constants at temperature and phi (electric potential)
|
||||
_update_rates_T();
|
||||
// get updated activities (rates updated below)
|
||||
_update_rates_C();
|
||||
|
||||
updateExchangeCurrentQuantities();
|
||||
|
||||
RxnMolChange* rmc = rmcVector[irxn];
|
||||
// could also get this from reactant and product stoichiometry, maybe
|
||||
nStoichElectrons = - rmc->m_phaseChargeChange[metalPhaseRS_];
|
||||
|
||||
OCV = 0.0;
|
||||
|
||||
getDeltaGibbs(0);
|
||||
|
||||
if (nStoichElectrons != 0.0) {
|
||||
OCV = m_deltaG[irxn] / Faraday / nStoichElectrons;
|
||||
}
|
||||
|
||||
for (size_t i = 0; i < m_ctrxn.size(); i++) {
|
||||
if (m_ctrxn[i] == irxn) {
|
||||
iBeta = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
beta = m_beta[iBeta];
|
||||
|
||||
doublereal rt = GasConstant*thermo(0).temperature();
|
||||
|
||||
|
||||
double mG0 = m_deltaG0[irxn];
|
||||
int reactionType = reactionTypes_[irxn];
|
||||
|
||||
//
|
||||
// Start with the forward reaction rate
|
||||
//
|
||||
double iO = m_rfn[irxn] * m_perturb[irxn];
|
||||
int iECDFormulation = m_ctrxn_ecdf[iBeta];
|
||||
if (! iECDFormulation) {
|
||||
iO = m_rfn[irxn] * Faraday * nStoichElectrons;
|
||||
if (beta > 0.0) {
|
||||
double fac = exp(mG0 / (rt));
|
||||
iO *= pow(fac, beta);
|
||||
// Need this step because m_rfn includes the inverse of this term, while the formulas
|
||||
// only use the chemical reaction rate constant.
|
||||
fac = exp( beta * deltaElectricEnergy_[irxn] / (rt));
|
||||
iO *= fac;
|
||||
}
|
||||
} else {
|
||||
iO *= nStoichElectrons;
|
||||
}
|
||||
|
||||
double omb = 1.0 - beta;
|
||||
if (reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN) {
|
||||
const RxnOrders* ro_fwd = m_ctrxn_ROPOrdersList_[iBeta];
|
||||
if (ro_fwd == 0) {
|
||||
throw CanteraError("ElectrodeKinetics::calcForwardROP_BV()", "forward orders pointer is zero ?!?");
|
||||
}
|
||||
double 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 ks = kinSpeciesIDs[j];
|
||||
thermo_t& th = speciesPhase(ks);
|
||||
size_t n = speciesPhaseIndex(ks);
|
||||
size_t klocal = ks - m_start[n];
|
||||
double mfS = th.moleFraction(klocal);
|
||||
|
||||
double oo = kinSpeciesOrders[j];
|
||||
tmp *= pow(mfS, oo);
|
||||
}
|
||||
iO *= tmp;
|
||||
} else {
|
||||
for (size_t k = 0; k < m_kk; k++) {
|
||||
doublereal reactCoeff = reactantStoichCoeff(k, irxn);
|
||||
doublereal prodCoeff = productStoichCoeff(k, irxn);
|
||||
|
||||
if (reactCoeff != 0.0) {
|
||||
iO *= pow(m_actConc[k], reactCoeff*omb);
|
||||
iO *= pow(m_StandardConc[k], reactCoeff*beta);
|
||||
}
|
||||
if (prodCoeff != 0.0) {
|
||||
iO *= pow(m_actConc[k], prodCoeff*beta);
|
||||
iO /= pow(m_StandardConc[k], prodCoeff*omb);
|
||||
}
|
||||
}
|
||||
}
|
||||
io = iO;
|
||||
resistivity = m_ctrxn_resistivity_[iBeta];
|
||||
|
||||
double phiMetal = m_thermo[metalPhaseRS_]->electricPotential();
|
||||
double phiSoln = m_thermo[solnPhaseRS_]->electricPotential();
|
||||
double E = phiMetal - phiSoln;
|
||||
overPotential = E - OCV;
|
||||
|
||||
return true;
|
||||
}
|
||||
//====================================================================================================================
|
||||
double ElectrodeKinetics::calcCurrentDensity(double nu, double nStoich, double ioc, double beta, double temp,
|
||||
doublereal resistivity) const
|
||||
{
|
||||
double exp1 = nu * nStoich * Faraday * beta / (GasConstant * temp);
|
||||
double exp2 = -nu * nStoich * Faraday * (1.0 - beta) / (GasConstant * temp);
|
||||
double val = ioc * (exp(exp1) - exp(exp2));
|
||||
if (resistivity > 0.0) {
|
||||
val = solveCurrentRes(nu, nStoich, ioc, beta, temp, resistivity, 0);
|
||||
}
|
||||
return val;
|
||||
}
|
||||
//==================================================================================================================
|
||||
void ElectrodeKinetics::init()
|
||||
{
|
||||
InterfaceKinetics::init();
|
||||
identifyMetalPhase();
|
||||
}
|
||||
//==================================================================================================================
|
||||
|
||||
double ElectrodeKinetics::solveCurrentRes(double nu, double nStoich, doublereal ioc, doublereal beta, doublereal temp,
|
||||
doublereal resistivity, int iprob) const
|
||||
{
|
||||
// int nits = 0;
|
||||
doublereal f, dfdi, deltai, eexp1, eexp2, exp1, exp2, icurr, deltai_damp;
|
||||
doublereal nFRT = nStoich * Faraday / (GasConstant * temp);
|
||||
if (iprob == 0) {
|
||||
eexp1 = exp(nu * nFRT * beta);
|
||||
eexp2 = exp(-nu * nFRT * (1.0 - beta)) ;
|
||||
|
||||
} else {
|
||||
eexp1 = exp(nu * nFRT * beta);
|
||||
eexp2 = 0.0;
|
||||
}
|
||||
icurr = ioc * (eexp1 - eexp2);
|
||||
double icurrDamp = icurr;
|
||||
if (fabs(resistivity * icurr) > 0.9 * fabs(nu)) {
|
||||
icurrDamp = 0.9 * nu / resistivity;
|
||||
}
|
||||
if (iprob == 0) {
|
||||
eexp1 = exp( nFRT * beta * (nu - resistivity * icurrDamp));
|
||||
eexp2 = exp(- nFRT * (1.0 - beta) * (nu - resistivity * icurrDamp));
|
||||
} else {
|
||||
eexp1 = exp( nFRT * beta * (nu - resistivity * icurrDamp));
|
||||
eexp2 = 0.0;
|
||||
}
|
||||
icurr = ioc * (eexp1 - eexp2);
|
||||
if (fabs(resistivity * icurr) > 0.99 * fabs(nu)) {
|
||||
icurr = 0.99 * nu / resistivity;
|
||||
}
|
||||
|
||||
do {
|
||||
// nits++;
|
||||
if (iprob == 0) {
|
||||
exp1 = nFRT * beta * (nu - resistivity * icurr);
|
||||
exp2 = - nFRT * (1.0 - beta) * (nu - resistivity * icurr);
|
||||
eexp1 = exp(exp1);
|
||||
eexp2 = exp(exp2);
|
||||
f = icurr - ioc * (eexp1 - eexp2);
|
||||
dfdi = 1.0 - ioc * eexp1 * ( - beta * nFRT * resistivity ) +
|
||||
ioc * eexp2 * ( (1.0 - beta) * nFRT * resistivity );
|
||||
} else {
|
||||
exp1 = nFRT * beta * (nu - resistivity * icurr);
|
||||
eexp1 = exp(exp1);
|
||||
f = icurr - ioc * (eexp1);
|
||||
dfdi = 1.0 - ioc * eexp1 * ( - beta * nFRT * resistivity );
|
||||
}
|
||||
deltai = - f / dfdi;
|
||||
if (fabs(deltai) > 0.1 * fabs(icurr)) {
|
||||
deltai_damp = 0.1 * deltai;
|
||||
if (fabs(deltai_damp) > 0.1 * fabs(icurr)) {
|
||||
deltai_damp = 0.1 * icurr * (deltai_damp / fabs(deltai_damp));
|
||||
}
|
||||
} else if (fabs(deltai) > 0.01 * fabs(icurr)) {
|
||||
deltai_damp = 0.3 * deltai;
|
||||
} else if (fabs(deltai) > 0.001 * fabs(icurr)) {
|
||||
deltai_damp = 0.5 * deltai;
|
||||
} else {
|
||||
deltai_damp = deltai;
|
||||
}
|
||||
icurr += deltai_damp;
|
||||
if (fabs(resistivity * icurr) > fabs(nu)) {
|
||||
icurr = 0.999 * nu / resistivity;
|
||||
}
|
||||
|
||||
} while((fabs(deltai/icurr)> 1.0E-14) && (fabs(deltai) > 1.0E-20));
|
||||
|
||||
// printf(" its = %d\n", nits);
|
||||
|
||||
return icurr;
|
||||
}
|
||||
//==================================================================================================================
|
||||
}
|
||||
|
|
|
|||
|
|
@ -6,9 +6,10 @@
|
|||
|
||||
#include "cantera/kinetics/InterfaceKinetics.h"
|
||||
#include "cantera/kinetics/EdgeKinetics.h"
|
||||
#include "cantera/kinetics/ReactionData.h"
|
||||
#include "cantera/kinetics/RateCoeffMgr.h"
|
||||
#include "cantera/kinetics/ImplicitSurfChem.h"
|
||||
#include "cantera/thermo/SurfPhase.h"
|
||||
#include "cantera/base/vec_functions.h"
|
||||
|
||||
#include <cstdio>
|
||||
|
||||
|
|
@ -221,6 +222,9 @@ void InterfaceKinetics::setElectricPotential(int n, doublereal V)
|
|||
//============================================================================================================================
|
||||
void InterfaceKinetics::_update_rates_T()
|
||||
{
|
||||
//
|
||||
// First task is update the electrical potentials from the Phases
|
||||
//
|
||||
_update_rates_phi();
|
||||
if (m_has_coverage_dependence) {
|
||||
m_surf->getCoverages(DATA_PTR(m_actConc));
|
||||
|
|
@ -269,6 +273,9 @@ void InterfaceKinetics::_update_rates_phi()
|
|||
}
|
||||
}
|
||||
//============================================================================================================================
|
||||
//
|
||||
// Updates the internal variables m_actConc and m_conc
|
||||
//
|
||||
void InterfaceKinetics::_update_rates_C()
|
||||
{
|
||||
for (size_t n = 0; n < nPhases(); n++) {
|
||||
|
|
@ -327,6 +334,12 @@ void InterfaceKinetics::updateKc()
|
|||
//============================================================================================================================
|
||||
void InterfaceKinetics::updateMu0()
|
||||
{
|
||||
//
|
||||
// First task is update the electrical potentials from the Phases
|
||||
//
|
||||
_update_rates_phi();
|
||||
|
||||
updateExchangeCurrentQuantities();
|
||||
/*
|
||||
* Get the vector of standard state electrochemical potentials for species in the Interfacial
|
||||
* kinetics object and store it in m_mu0[] and in m_mu0_Kc[]
|
||||
|
|
@ -347,6 +360,11 @@ void InterfaceKinetics::updateMu0()
|
|||
//============================================================================================================================
|
||||
void InterfaceKinetics::checkPartialEquil()
|
||||
{
|
||||
//
|
||||
// First task is update the electrical potentials from the Phases
|
||||
//
|
||||
_update_rates_phi();
|
||||
|
||||
vector_fp dmu(nTotalSpecies(), 0.0);
|
||||
vector_fp rmu(std::max<size_t>(nReactions(), 1), 0.0);
|
||||
if (m_nrev > 0) {
|
||||
|
|
@ -424,7 +442,7 @@ void InterfaceKinetics::updateExchangeCurrentQuantities()
|
|||
* First collect vectors of the standard Gibbs free energies of the
|
||||
* species and the standard concentrations
|
||||
* - m_mu0
|
||||
* - m_logStandardConc
|
||||
* - m_StandardConc
|
||||
*/
|
||||
size_t ik = 0;
|
||||
|
||||
|
|
@ -1026,6 +1044,11 @@ void InterfaceKinetics::addElementaryReaction(ReactionData& rdata)
|
|||
}
|
||||
m_ctrxn_ROPOrdersList_.push_back(0);
|
||||
m_ctrxn_FwdOrdersList_.push_back(0);
|
||||
if (rdata.filmResistivity > 0.0) {
|
||||
throw CanteraError("InterfaceKinetics::addElementaryReaction()",
|
||||
"film resistivity set for elementary reaction");
|
||||
}
|
||||
m_ctrxn_resistivity_.push_back(rdata.filmResistivity);
|
||||
}
|
||||
|
||||
// add constant term to rate coeff value vector
|
||||
|
|
@ -1085,16 +1108,29 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata)
|
|||
if (rdata.beta > 0.0 || 1) {
|
||||
m_has_electrochem_rxns = true;
|
||||
m_beta.push_back(rdata.beta);
|
||||
// Push back the id of the reaction
|
||||
// Push back the id of the reaction
|
||||
m_ctrxn.push_back(m_ii);
|
||||
// set the default to be the normal forward / reverse calculation method
|
||||
m_ctrxn_BVform.push_back(0);
|
||||
//
|
||||
// Specify alternative forms of the electrochemical reaction
|
||||
//
|
||||
if (rdata.reactionType == BUTLERVOLMER_RXN) {
|
||||
m_ctrxn_BVform.push_back(1);
|
||||
} else if (rdata.reactionType == BUTLERVOLMER_NOACTIVITYCOEFFS_RXN) {
|
||||
m_ctrxn_BVform.push_back(2);
|
||||
} else {
|
||||
// set the default to be the normal forward / reverse calculation method
|
||||
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);
|
||||
}
|
||||
//
|
||||
// Store the film resistivity
|
||||
//
|
||||
m_ctrxn_resistivity_.push_back(rdata.filmResistivity);
|
||||
|
||||
if (rdata.forwardFullOrder_.size() > 0) {
|
||||
RxnOrders* ro = new RxnOrders();
|
||||
|
|
@ -1109,7 +1145,7 @@ void InterfaceKinetics::addGlobalReaction(ReactionData& rdata)
|
|||
std::vector<double> fwdFullorders(m_kk, 0.0);
|
||||
determineFwdOrdersBV(rdata, fwdFullorders);
|
||||
RxnOrders* ro = new RxnOrders();
|
||||
ro->fill(rdata.forwardFullOrder_);
|
||||
ro->fill(fwdFullorders);
|
||||
m_ctrxn_FwdOrdersList_[m_ii] = ro;
|
||||
}
|
||||
} else {
|
||||
|
|
@ -1258,12 +1294,15 @@ void InterfaceKinetics::finalize()
|
|||
size_t ks = reactionPhaseIndex();
|
||||
if (ks == npos) throw CanteraError("InterfaceKinetics::finalize",
|
||||
"no surface phase is present.");
|
||||
//
|
||||
// Check to see that the interface routine has a dimension of 2
|
||||
//
|
||||
m_surf = (SurfPhase*)&thermo(ks);
|
||||
if (m_surf->nDim() != 2)
|
||||
if (m_surf->nDim() != 2) {
|
||||
throw CanteraError("InterfaceKinetics::finalize",
|
||||
"expected interface dimension = 2, but got dimension = "
|
||||
+int2str(m_surf->nDim()));
|
||||
|
||||
}
|
||||
m_StandardConc.resize(m_kk, 0.0);
|
||||
m_deltaG0.resize(safe_reaction_size, 0.0);
|
||||
m_deltaG.resize(safe_reaction_size, 0.0);
|
||||
|
|
@ -1283,6 +1322,13 @@ void InterfaceKinetics::finalize()
|
|||
m_ropnet.resize(1, 0.0);
|
||||
m_rkcn.resize(1, 0.0);
|
||||
}
|
||||
//
|
||||
// Malloc and calculate all of the quantities that go into the extra description of reactions
|
||||
//
|
||||
rmcVector.resize(m_ii, 0);
|
||||
for (size_t i = 0; i < m_ii; i++) {
|
||||
rmcVector[i] = new RxnMolChange(this, i);
|
||||
}
|
||||
|
||||
m_finalized = true;
|
||||
}
|
||||
|
|
@ -1425,20 +1471,25 @@ void InterfaceKinetics::determineFwdOrdersBV(ReactionData& rdata, std::vector<do
|
|||
double betaf = rdata.beta;
|
||||
double betar = 1.0 - betaf;
|
||||
//
|
||||
// Loop over the reactants doing away the BV terms.
|
||||
// Loop over the reactants doing away with 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];
|
||||
double oo = rdata.rstoich[j];
|
||||
fwdFullorders[kkin] += betaf * oo;
|
||||
// just to make sure roundoff doesn't leave a term that should be zero (haven't checked this out yet)
|
||||
if (abs(fwdFullorders[kkin]) < 0.00001) {
|
||||
fwdFullorders[kkin] = 0.0;
|
||||
}
|
||||
}
|
||||
//
|
||||
// Loop over the products doing away with the BV terms.
|
||||
// This should leave the reactant terms only, even if they are non-mass action.
|
||||
//
|
||||
for (size_t j = 0; j < rdata.products.size(); j++) {
|
||||
size_t kkin = rdata.products[j];
|
||||
double oo = rdata.pstoich[kkin];
|
||||
double oo = rdata.pstoich[j];
|
||||
fwdFullorders[kkin] -= betaf * oo;
|
||||
if (abs(fwdFullorders[kkin]) < 0.00001) {
|
||||
fwdFullorders[kkin] = 0.0;
|
||||
|
|
@ -1448,15 +1499,36 @@ void InterfaceKinetics::determineFwdOrdersBV(ReactionData& rdata, std::vector<do
|
|||
//==================================================================================================================
|
||||
void EdgeKinetics::finalize()
|
||||
{
|
||||
deltaElectricEnergy_.resize(std::max<size_t>(m_ii, 1));
|
||||
//
|
||||
// Note we can't call the Interface::finalize() routine because we need to check for a dimension of 1 below.
|
||||
// Therefore, we have to malloc room in arrays that would normally be
|
||||
// handled by the InterfaceKinetics::finalize() call.
|
||||
//
|
||||
Kinetics::finalize();
|
||||
|
||||
size_t safe_reaction_size = std::max<size_t>(m_ii, 1);
|
||||
deltaElectricEnergy_.resize(safe_reaction_size);
|
||||
size_t ks = reactionPhaseIndex();
|
||||
if (ks == npos) throw CanteraError("EdgeKinetics::finalize",
|
||||
"no edge phase is present.");
|
||||
"no surface phase is present.");
|
||||
//
|
||||
// Check to see edge phase has a dimension of 1
|
||||
//
|
||||
m_surf = (SurfPhase*)&thermo(ks);
|
||||
if (m_surf->nDim() != 1)
|
||||
if (m_surf->nDim() != 1) {
|
||||
throw CanteraError("EdgeKinetics::finalize",
|
||||
"expected interface dimension = 1, but got dimension = "
|
||||
+int2str(m_surf->nDim()));
|
||||
}
|
||||
m_StandardConc.resize(m_kk, 0.0);
|
||||
m_deltaG0.resize(safe_reaction_size, 0.0);
|
||||
m_deltaG.resize(safe_reaction_size, 0.0);
|
||||
|
||||
m_ProdStanConcReac.resize(safe_reaction_size, 0.0);
|
||||
|
||||
if (m_thermo.size() != m_phaseExists.size()) {
|
||||
throw CanteraError("InterfaceKinetics::finalize", "internal error");
|
||||
}
|
||||
|
||||
// Guarantee that these arrays can be converted to double* even in the
|
||||
// special case where there are no reactions defined.
|
||||
|
|
@ -1467,7 +1539,6 @@ void EdgeKinetics::finalize()
|
|||
m_ropnet.resize(1, 0.0);
|
||||
m_rkcn.resize(1, 0.0);
|
||||
}
|
||||
|
||||
//
|
||||
// Malloc and calculate all of the quantities that go into the extra description of reactions
|
||||
//
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue