958 lines
32 KiB
C++
958 lines
32 KiB
C++
/**
|
|
* @file ElectrodeKinetics.cpp
|
|
*/
|
|
|
|
#include "cantera/kinetics/ElectrodeKinetics.h"
|
|
#include "cantera/thermo/SurfPhase.h"
|
|
#include "cantera/base/utilities.h"
|
|
|
|
#include <cstdio>
|
|
|
|
using namespace std;
|
|
|
|
namespace Cantera
|
|
{
|
|
//============================================================================================================================
|
|
ElectrodeKinetics::ElectrodeKinetics(thermo_t* thermo) :
|
|
InterfaceKinetics(thermo),
|
|
metalPhaseRS_(npos),
|
|
solnPhaseRS_(npos),
|
|
kElectronRS_(npos)
|
|
{
|
|
|
|
}
|
|
//============================================================================================================================
|
|
ElectrodeKinetics::~ElectrodeKinetics()
|
|
{
|
|
|
|
}
|
|
//============================================================================================================================
|
|
ElectrodeKinetics::ElectrodeKinetics(const ElectrodeKinetics& right) :
|
|
InterfaceKinetics()
|
|
|
|
{
|
|
/*
|
|
* Call the assignment operator
|
|
*/
|
|
ElectrodeKinetics::operator=(right);
|
|
}
|
|
//============================================================================================================================
|
|
ElectrodeKinetics& ElectrodeKinetics::operator=(const ElectrodeKinetics& right)
|
|
{
|
|
/*
|
|
* Check for self assignment.
|
|
*/
|
|
if (this == &right) {
|
|
return *this;
|
|
}
|
|
|
|
InterfaceKinetics::operator=(right);
|
|
|
|
metalPhaseRS_ = right.metalPhaseRS_;
|
|
solnPhaseRS_ = right.solnPhaseRS_;
|
|
kElectronRS_ = right.kElectronRS_;
|
|
|
|
return *this;
|
|
}
|
|
//============================================================================================================================
|
|
int ElectrodeKinetics::type() const
|
|
{
|
|
return cInterfaceKinetics;
|
|
}
|
|
//============================================================================================================================
|
|
Kinetics* ElectrodeKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const
|
|
{
|
|
ElectrodeKinetics* iK = new ElectrodeKinetics(*this);
|
|
iK->assignShallowPointers(tpVector);
|
|
return iK;
|
|
}
|
|
//============================================================================================================================
|
|
// Identify the metal phase and the electron species
|
|
void ElectrodeKinetics::identifyMetalPhase()
|
|
{
|
|
metalPhaseRS_ = npos;
|
|
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
|
|
// Should probably also stipulate a charge of -1.
|
|
//
|
|
for (size_t iph = 0; iph < np; iph++) {
|
|
ThermoPhase* tp = m_thermo[iph];
|
|
size_t nSpecies = tp->nSpecies();
|
|
size_t nElements = tp->nElements();
|
|
size_t eElectron = tp->elementIndex("E");
|
|
if (eElectron != npos) {
|
|
for (size_t k = 0; k < nSpecies; k++) {
|
|
if (tp->nAtoms(k,eElectron) == 1) {
|
|
int ifound = 1;
|
|
for (size_t e = 0; e < nElements; e++) {
|
|
if (tp->nAtoms(k,e) != 0.0) {
|
|
if (e != eElectron) {
|
|
ifound = 0;
|
|
}
|
|
}
|
|
}
|
|
if (ifound == 1) {
|
|
metalPhaseRS_ = iph;
|
|
kElectronRS_ = m_start[iph] + k;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
//
|
|
// 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];
|
|
if (rmc->m_phaseChargeChange[iph] != 0) {
|
|
if (rmc->m_phaseDims[iph] == 3) {
|
|
solnPhaseRS_ = iph;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
*/
|
|
//
|
|
// 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, we will not error exit. Some functions will
|
|
// be turned off and the object will behave as an InterfaceKinetics object. This is needed
|
|
// because downstream electrode objects have internal reaction surfaces that don't have
|
|
// electrons.
|
|
//
|
|
/*
|
|
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
|
|
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
|
|
//
|
|
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
|
|
//
|
|
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
|
|
//
|
|
AssertThrow(metalPhaseRS_ != npos, "ElectrodeKinetics::updateROP()");
|
|
|
|
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;
|
|
|
|
//
|
|
// 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 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);
|
|
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(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];
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
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;
|
|
}
|
|
//==================================================================================================================
|
|
//
|
|
// 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 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?)
|
|
//
|
|
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
|
|
// (we don't use standard concentrations at all here);
|
|
//
|
|
double tmp = exp(- m_beta[iBeta] * m_deltaG0[irxn] * rrt);
|
|
double tmp2 = 1.0;
|
|
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 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);
|
|
}
|
|
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
|
|
if (metalPhaseRS_ == npos) {
|
|
nStoichElectrons = 0;
|
|
OCV = 0.0;
|
|
return false;
|
|
} else {
|
|
nStoichElectrons = - rmc->m_phaseChargeChange[metalPhaseRS_];
|
|
}
|
|
|
|
|
|
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 if (reactionType == BUTLERVOLMER_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((m_actConc[ks]/m_StandardConc[ks]), 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;
|
|
}
|
|
//==================================================================================================================
|
|
}
|