initial import
This commit is contained in:
Dave Goodwin
parent
2e6c5ad2a1
commit
d86a6bf097
3 changed files with 942 additions and 0 deletions
506
Cantera/src/EdgeKinetics.cpp
Normal file
506
Cantera/src/EdgeKinetics.cpp
Normal file
|
|
@ -0,0 +1,506 @@
|
|||
/**
|
||||
* @file EdgeKinetics.cpp
|
||||
*
|
||||
*/
|
||||
|
||||
// Copyright 2002 California Institute of Technology
|
||||
|
||||
|
||||
// turn off warnings under Windows
|
||||
#ifdef WIN32
|
||||
#pragma warning(disable:4786)
|
||||
#pragma warning(disable:4503)
|
||||
#endif
|
||||
|
||||
#include "EdgeKinetics.h"
|
||||
#include "SurfPhase.h"
|
||||
|
||||
#include "ReactionData.h"
|
||||
#include "StoichManager.h"
|
||||
#include "RateCoeffMgr.h"
|
||||
|
||||
#include <iostream>
|
||||
using namespace std;
|
||||
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
//////////////////////////////////////////////////////////////////
|
||||
|
||||
/**
|
||||
* Construct an empty EdgeKinetics reaction mechanism.
|
||||
* @param thermo This is an optional parameter that may be
|
||||
* used to initialize the inherited Kinetics class with
|
||||
* one ThermoPhase class object -> in other words it's
|
||||
* useful for initialization of homogeneous kinetics
|
||||
* mechanisms.
|
||||
*/
|
||||
EdgeKinetics::
|
||||
EdgeKinetics(thermo_t* thermo) :
|
||||
Kinetics(thermo),
|
||||
m_kk(0),
|
||||
m_redo_rates(false),
|
||||
m_nirrev(0),
|
||||
m_nrev(0),
|
||||
m_finalized(false)
|
||||
{
|
||||
m_kdata = new EdgeKineticsData;
|
||||
m_kdata->m_temp = 0.0;
|
||||
}
|
||||
|
||||
/**
|
||||
* Destructor
|
||||
*/
|
||||
EdgeKinetics::
|
||||
~EdgeKinetics(){
|
||||
delete m_kdata;
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Update properties that depend on temperature
|
||||
*
|
||||
*/
|
||||
void EdgeKinetics::
|
||||
_update_rates_T() {
|
||||
_update_rates_phi();
|
||||
doublereal T = thermo(surfacePhaseIndex()).temperature();
|
||||
if (T != m_kdata->m_temp || m_redo_rates) {
|
||||
m_kdata->m_logtemp = log(T);
|
||||
m_rates.update(T, m_kdata->m_logtemp, m_kdata->m_rfn.begin());
|
||||
applyButlerVolmerCorrection(m_kdata->m_rfn.begin());
|
||||
m_kdata->m_temp = T;
|
||||
updateKc();
|
||||
m_kdata->m_ROP_ok = false;
|
||||
m_redo_rates = false;
|
||||
}
|
||||
}
|
||||
|
||||
void EdgeKinetics::
|
||||
_update_rates_phi() {
|
||||
int np = nPhases();
|
||||
for (int n = 0; n < np; n++) {
|
||||
if (thermo(n).electricPotential() != m_phi[n]) {
|
||||
m_phi[n] = thermo(n).electricPotential();
|
||||
m_redo_rates = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Update properties that depend on concentrations. This method
|
||||
* fills out the array of generalized concentrations by calling
|
||||
* method getActivityConcentrations for each phase, which classes
|
||||
* representing phases should overload to return the appropriate
|
||||
* quantities.
|
||||
*/
|
||||
void EdgeKinetics::
|
||||
_update_rates_C() {
|
||||
int n;
|
||||
|
||||
//m_rates.update(m_kdata->m_temp,
|
||||
// m_kdata->m_logtemp, m_kdata->m_rfn.begin());
|
||||
|
||||
int np = nPhases();
|
||||
for (n = 0; n < np; n++) {
|
||||
thermo(n).getActivityConcentrations(m_conc.begin() + m_start[n]);
|
||||
}
|
||||
m_kdata->m_ROP_ok = false;
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Update the equilibrium constants in molar units for all
|
||||
* reversible reactions. Irreversible reactions have their
|
||||
* equilibrium constant set to zero.
|
||||
*/
|
||||
void EdgeKinetics::updateKc() {
|
||||
int i, irxn;
|
||||
vector_fp& m_rkc = m_kdata->m_rkcn;
|
||||
fill(m_rkc.begin(), m_rkc.end(), 0.0);
|
||||
|
||||
if (m_nrev > 0) {
|
||||
|
||||
int n, nsp, k, ik=0;
|
||||
doublereal rt = GasConstant*thermo(0).temperature();
|
||||
doublereal rrt = 1.0/rt;
|
||||
int np = nPhases();
|
||||
for (n = 0; n < np; n++) {
|
||||
thermo(n).getStandardChemPotentials(m_mu0.begin() + m_start[n]);
|
||||
nsp = thermo(n).nSpecies();
|
||||
for (k = 0; k < nsp; k++) {
|
||||
m_mu0[ik] -= rt*thermo(n).logStandardConc(k);
|
||||
m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k);
|
||||
ik++;
|
||||
}
|
||||
}
|
||||
|
||||
// compute Delta mu^0 for all reversible reactions
|
||||
m_reactantStoich.decrementReactions(m_mu0.begin(), m_rkc.begin());
|
||||
m_revProductStoich.incrementReactions(m_mu0.begin(), m_rkc.begin());
|
||||
|
||||
for (i = 0; i < m_nrev; i++) {
|
||||
irxn = m_revindex[i];
|
||||
m_rkc[irxn] = exp(m_rkc[irxn]*rrt);
|
||||
}
|
||||
for (i = 0; i != m_nirrev; ++i) {
|
||||
m_rkc[ m_irrev[i] ] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
void EdgeKinetics::checkPartialEquil() {
|
||||
int i, irxn;
|
||||
vector_fp dmu(nTotalSpecies(), 0.0);
|
||||
vector_fp rmu(nReactions(), 0.0);
|
||||
if (m_nrev > 0) {
|
||||
|
||||
int n, nsp, k, ik=0;
|
||||
doublereal rt = GasConstant*thermo(0).temperature();
|
||||
doublereal rrt = 1.0/rt;
|
||||
int np = nPhases();
|
||||
for (n = 0; n < np; n++) {
|
||||
thermo(n).getChemPotentials(dmu.begin() + m_start[n]);
|
||||
nsp = thermo(n).nSpecies();
|
||||
for (k = 0; k < nsp; k++) {
|
||||
dmu[ik] += Faraday * m_phi[n] * thermo(n).charge(k);
|
||||
cout << thermo(n).speciesName(k) << " " << dmu[ik] << endl;
|
||||
ik++;
|
||||
}
|
||||
}
|
||||
|
||||
// compute Delta mu^ for all reversible reactions
|
||||
m_reactantStoich.decrementReactions(dmu.begin(), rmu.begin());
|
||||
m_revProductStoich.incrementReactions(dmu.begin(), rmu.begin());
|
||||
|
||||
for (i = 0; i < m_nrev; i++) {
|
||||
irxn = m_revindex[i];
|
||||
cout << "Reaction " << irxn << " " << exp(rmu[irxn]*rrt) << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Get the equilibrium constants of all reactions, whether
|
||||
* reversible or not.
|
||||
*/
|
||||
void EdgeKinetics::getEquilibriumConstants(doublereal* kc) {
|
||||
int i;
|
||||
|
||||
int n, nsp, k, ik=0;
|
||||
doublereal rt = GasConstant*thermo(0).temperature();
|
||||
doublereal rrt = 1.0/rt;
|
||||
int np = nPhases();
|
||||
for (n = 0; n < np; n++) {
|
||||
thermo(n).getStandardChemPotentials(m_mu0.begin() + m_start[n]);
|
||||
nsp = thermo(n).nSpecies();
|
||||
for (k = 0; k < nsp; k++) {
|
||||
m_mu0[ik] -= rt*thermo(n).logStandardConc(k);
|
||||
m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k);
|
||||
ik++;
|
||||
}
|
||||
}
|
||||
|
||||
fill(kc, kc + m_ii, 0.0);
|
||||
|
||||
m_reactantStoich.decrementReactions(m_mu0.begin(), kc);
|
||||
m_revProductStoich.incrementReactions(m_mu0.begin(), kc);
|
||||
m_irrevProductStoich.incrementReactions(m_mu0.begin(), kc);
|
||||
|
||||
for (i = 0; i < m_ii; i++) {
|
||||
kc[i] = exp(-kc[i]*rrt);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* 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.
|
||||
*/
|
||||
void EdgeKinetics::applyButlerVolmerCorrection(doublereal* kf) {
|
||||
int i;
|
||||
|
||||
int n, nsp, k, ik=0;
|
||||
doublereal rt = GasConstant*thermo(0).temperature();
|
||||
doublereal rrt = 1.0/rt;
|
||||
int np = nPhases();
|
||||
|
||||
// compute the electrical potential energy of each species
|
||||
for (n = 0; n < np; n++) {
|
||||
nsp = thermo(n).nSpecies();
|
||||
for (k = 0; k < nsp; k++) {
|
||||
m_pot[ik] = Faraday*thermo(n).charge(k)*m_phi[n];
|
||||
ik++;
|
||||
}
|
||||
}
|
||||
|
||||
// compute the change in electrical potential energy for each
|
||||
// reaction. This will only be non-zero if a potential
|
||||
// difference is present.
|
||||
fill(m_rwork.begin(), m_rwork.begin() + m_ii, 0.0);
|
||||
m_reactantStoich.decrementReactions(m_pot.begin(), m_rwork.begin());
|
||||
m_revProductStoich.incrementReactions(m_pot.begin(), m_rwork.begin());
|
||||
m_irrevProductStoich.incrementReactions(m_pot.begin(), m_rwork.begin());
|
||||
|
||||
// modify the reaction rates. Only modify those with a
|
||||
// non-zero activation energy, and do not decrease the
|
||||
// activation energy below zero.
|
||||
doublereal ea, eamod;
|
||||
|
||||
for (i = 0; i < m_ii; i++) {
|
||||
eamod = 0.5*m_rwork[i];
|
||||
if (eamod != 0.0 && m_E[i] != 0.0) {
|
||||
ea = GasConstant * m_E[i];
|
||||
if (eamod + ea < 0.0) eamod = -ea;
|
||||
kf[i] *= exp(-eamod*rrt);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Update the rates of progress of the reactions in the reaciton
|
||||
* mechanism. This routine operates on internal data.
|
||||
*/
|
||||
void EdgeKinetics::updateROP() {
|
||||
|
||||
_update_rates_T();
|
||||
_update_rates_C();
|
||||
|
||||
if (m_kdata->m_ROP_ok) return;
|
||||
|
||||
const vector_fp& rf = m_kdata->m_rfn;
|
||||
const vector_fp& m_rkc = m_kdata->m_rkcn;
|
||||
array_fp& ropf = m_kdata->m_ropf;
|
||||
array_fp& ropr = m_kdata->m_ropr;
|
||||
array_fp& ropnet = m_kdata->m_ropnet;
|
||||
|
||||
// copy rate coefficients into ropf
|
||||
copy(rf.begin(), rf.end(), ropf.begin());
|
||||
|
||||
// multiply by perturbation factor
|
||||
multiply_each(ropf.begin(), ropf.end(), m_perturb.begin());
|
||||
|
||||
// copy the forward rates to the reverse rates
|
||||
copy(ropf.begin(), ropf.end(), 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(ropr.begin(), ropr.end(), m_rkc.begin());
|
||||
|
||||
// multiply ropf by concentration products
|
||||
m_reactantStoich.multiply(m_conc.begin(), ropf.begin());
|
||||
|
||||
// for reversible reactions, multiply ropr by concentration
|
||||
// products
|
||||
m_revProductStoich.multiply(m_conc.begin(), ropr.begin());
|
||||
|
||||
// do global reactions
|
||||
m_globalReactantStoich.power(m_conc.begin(), ropf.begin());
|
||||
|
||||
for (int j = 0; j != m_ii; ++j) {
|
||||
ropnet[j] = ropf[j] - ropr[j];
|
||||
}
|
||||
|
||||
m_kdata->m_ROP_ok = true;
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Add a single reaction to the mechanism. This routine
|
||||
* must be called after init() and before finalize().
|
||||
* This function branches on the types of reactions allowed
|
||||
* by the interfaceKinetics manager in order to install
|
||||
* the reaction correctly in the manager.
|
||||
* The manager allows the following reaction types
|
||||
* Elementary
|
||||
* Surface
|
||||
* Global
|
||||
* There is no difference between elementary and surface
|
||||
* reactions.
|
||||
*/
|
||||
void EdgeKinetics::
|
||||
addReaction(const ReactionData& r) {
|
||||
|
||||
int nr = r.reactants.size();
|
||||
|
||||
// a global reaction is idnetified as one with
|
||||
// a reactant stoichiometric coefficient not equal
|
||||
// to the molecularity for some reactant
|
||||
bool isglobal = false;
|
||||
for (int n = 0; n < nr; n++) {
|
||||
if (r.rstoich[n] != int(r.order[n])) {
|
||||
isglobal = true; break;
|
||||
}
|
||||
}
|
||||
if (isglobal)
|
||||
addGlobalReaction(r);
|
||||
else
|
||||
addElementaryReaction(r);
|
||||
|
||||
installReagents( r );
|
||||
installGroups(reactionNumber(), r.rgroups, r.pgroups);
|
||||
incrementRxnCount();
|
||||
m_rxneqn.push_back(r.equation);
|
||||
}
|
||||
|
||||
|
||||
void EdgeKinetics::
|
||||
addElementaryReaction(const ReactionData& r) {
|
||||
int iloc;
|
||||
// install rate coeff calculator
|
||||
vector_fp rp = r.rateCoeffParameters;
|
||||
int ncov = r.cov.size();
|
||||
for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]);
|
||||
iloc = m_rates.install( reactionNumber(),
|
||||
r.rateCoeffType, rp.size(),
|
||||
rp.begin() );
|
||||
// store activation energy
|
||||
m_E.push_back(r.rateCoeffParameters[2]);
|
||||
// add constant term to rate coeff value vector
|
||||
m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]);
|
||||
registerReaction( reactionNumber(), ELEMENTARY_RXN, iloc);
|
||||
}
|
||||
|
||||
|
||||
void EdgeKinetics::
|
||||
addGlobalReaction(const ReactionData& r) {
|
||||
|
||||
int iloc;
|
||||
// install rate coeff calculator
|
||||
vector_fp rp = r.rateCoeffParameters;
|
||||
int ncov = r.cov.size();
|
||||
for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]);
|
||||
iloc = m_rates.install( reactionNumber(),
|
||||
r.rateCoeffType, rp.size(),
|
||||
rp.begin() );
|
||||
|
||||
// add constant term to rate coeff value vector
|
||||
m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]);
|
||||
|
||||
int nr = r.order.size();
|
||||
vector_fp ordr(nr);
|
||||
for (int n = 0; n < nr; n++) {
|
||||
ordr[n] = r.order[n] - r.rstoich[n];
|
||||
}
|
||||
m_globalReactantStoich.add( reactionNumber(),
|
||||
r.reactants, ordr);
|
||||
|
||||
registerReaction( reactionNumber(), GLOBAL_RXN, iloc);
|
||||
}
|
||||
|
||||
|
||||
void EdgeKinetics::installReagents(const ReactionData& r) {
|
||||
|
||||
m_kdata->m_ropf.push_back(0.0); // extend by one for new rxn
|
||||
m_kdata->m_ropr.push_back(0.0);
|
||||
m_kdata->m_ropnet.push_back(0.0);
|
||||
int n, ns, m;
|
||||
|
||||
int rnum = reactionNumber();
|
||||
|
||||
vector_int rk;
|
||||
int nr = r.reactants.size();
|
||||
for (n = 0; n < nr; n++) {
|
||||
ns = r.rstoich[n];
|
||||
m_rrxn[r.reactants[n]][rnum] = ns;
|
||||
for (m = 0; m < ns; m++) {
|
||||
rk.push_back(r.reactants[n]);
|
||||
}
|
||||
}
|
||||
m_reactants.push_back(rk);
|
||||
|
||||
vector_int pk;
|
||||
int np = r.products.size();
|
||||
for (n = 0; n < np; n++) {
|
||||
ns = r.pstoich[n];
|
||||
m_prxn[r.products[n]][rnum] = ns;
|
||||
for (m = 0; m < ns; m++) {
|
||||
pk.push_back(r.products[n]);
|
||||
}
|
||||
}
|
||||
m_products.push_back(pk);
|
||||
|
||||
m_kdata->m_rkcn.push_back(0.0);
|
||||
|
||||
m_reactantStoich.add( reactionNumber(), rk);
|
||||
|
||||
if (r.reversible) {
|
||||
m_revProductStoich.add(reactionNumber(), pk);
|
||||
//m_dn.push_back(pk.size() - rk.size());
|
||||
m_revindex.push_back(reactionNumber());
|
||||
m_nrev++;
|
||||
}
|
||||
else {
|
||||
m_irrevProductStoich.add(reactionNumber(), pk);
|
||||
//m_dn.push_back(pk.size() - rk.size());
|
||||
m_irrev.push_back( reactionNumber() );
|
||||
m_nirrev++;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void EdgeKinetics::installGroups(int irxn,
|
||||
const vector<grouplist_t>& r, const vector<grouplist_t>& p) {
|
||||
if (!r.empty()) {
|
||||
m_rgroups[reactionNumber()] = r;
|
||||
m_pgroups[reactionNumber()] = p;
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Prepare the class for the addition of reactions. This function
|
||||
* must be called after instantiation of the class, but before
|
||||
* any reactions are actually added to the mechanism.
|
||||
* This function calculates m_kk the number of species in all
|
||||
* phases participating in the reaction mechanism. We don't know
|
||||
* m_kk previously, before all phases have been added.
|
||||
*/
|
||||
void EdgeKinetics::init() {
|
||||
int n;
|
||||
m_kk = 0;
|
||||
int np = nPhases();
|
||||
for (n = 0; n < np; n++) {
|
||||
m_kk += thermo(n).nSpecies();
|
||||
}
|
||||
m_rrxn.resize(m_kk);
|
||||
m_prxn.resize(m_kk);
|
||||
m_conc.resize(m_kk);
|
||||
m_mu0.resize(m_kk);
|
||||
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
|
||||
* and before the mechanism is used to calculate reaction rates.
|
||||
*
|
||||
* Here, we resize work arrays based on the number of reactions,
|
||||
* since we don't know this number up to now.
|
||||
*/
|
||||
void EdgeKinetics::finalize() {
|
||||
m_rwork.resize(nReactions());
|
||||
m_finalized = true;
|
||||
}
|
||||
|
||||
|
||||
bool EdgeKinetics::ready() const {
|
||||
return (m_finalized);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
396
Cantera/src/EdgeKinetics.h
Normal file
396
Cantera/src/EdgeKinetics.h
Normal file
|
|
@ -0,0 +1,396 @@
|
|||
/**
|
||||
* @file EdgeKinetics.h
|
||||
*
|
||||
* $Author$
|
||||
* $Revision$
|
||||
* $Date$
|
||||
*/
|
||||
|
||||
// Copyright 2001 California Institute of Technology
|
||||
|
||||
|
||||
#ifndef CT_EDGEKINETICS_H
|
||||
#define CT_EDGEKINETICS_H
|
||||
|
||||
#include <fstream>
|
||||
#include <math.h>
|
||||
#include <map>
|
||||
#include <stdlib.h>
|
||||
|
||||
#include "mix_defs.h"
|
||||
#include "Kinetics.h"
|
||||
|
||||
#include "utilities.h"
|
||||
#include "RateCoeffMgr.h"
|
||||
#include "StoichManager.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
// forward references
|
||||
|
||||
class ReactionData;
|
||||
class EdgeKineticsData;
|
||||
class ThermoPhase;
|
||||
class SurfPhase;
|
||||
class ImplicitSurfChem;
|
||||
|
||||
/**
|
||||
* Holds mechanism-specific data.
|
||||
*/
|
||||
class EdgeKineticsData {
|
||||
public:
|
||||
EdgeKineticsData() :
|
||||
m_ROP_ok(false),
|
||||
m_temp(0.0), m_logtemp(0.0)
|
||||
{}
|
||||
virtual ~EdgeKineticsData(){}
|
||||
|
||||
doublereal m_logp0, m_logc0;
|
||||
array_fp m_ropf, m_ropr, m_ropnet;
|
||||
array_fp m_rfn_low, m_rfn_high;
|
||||
bool m_ROP_ok;
|
||||
|
||||
doublereal m_temp, m_logtemp;
|
||||
vector_fp m_rfn;
|
||||
vector_fp m_rkcn;
|
||||
};
|
||||
|
||||
|
||||
class EdgeKinetics : public Kinetics {
|
||||
|
||||
public:
|
||||
|
||||
/**
|
||||
* Constructor
|
||||
*
|
||||
* @param thermo The optional parameter may be used to initialize
|
||||
* the object with one ThermoPhase object.
|
||||
* HKM Note -> Since the interface kinetics
|
||||
* object will probably require multiple thermophase
|
||||
* objects, this is probably not a good idea
|
||||
* to have this parameter.
|
||||
*/
|
||||
EdgeKinetics(thermo_t* thermo = 0);
|
||||
|
||||
/// Destructor.
|
||||
virtual ~EdgeKinetics();
|
||||
|
||||
/**
|
||||
* Identifies the subclass of the Kinetics manager type.
|
||||
* These are listed in mix_defs.h.
|
||||
*/
|
||||
virtual int ID() { return cEdgeKinetics; }
|
||||
|
||||
/**
|
||||
* Identifies the subclass of the Kinetics manager type.
|
||||
* These are listed in mix_defs.h.
|
||||
*/
|
||||
virtual int type() { return cEdgeKinetics; }
|
||||
|
||||
/**
|
||||
* 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;
|
||||
}
|
||||
|
||||
//@}
|
||||
/**
|
||||
* @name Reaction Rates Of Progress
|
||||
*/
|
||||
//@{
|
||||
|
||||
/**
|
||||
* Forward rates of progress.
|
||||
* Return the forward rates of progress in array fwdROP, which
|
||||
* must be dimensioned at least as large as the total number
|
||||
* of reactions.
|
||||
* Units are kmol/m2/s
|
||||
*/
|
||||
virtual void getFwdRatesOfProgress(doublereal* fwdROP) {
|
||||
updateROP();
|
||||
copy(m_kdata->m_ropf.begin(), m_kdata->m_ropf.end(), fwdROP);
|
||||
}
|
||||
|
||||
/**
|
||||
* Reverse rates of progress.
|
||||
* Return the reverse rates of progress in array revROP, which
|
||||
* must be dimensioned at least as large as the total number
|
||||
* of reactions.
|
||||
* Units are kmol/m2/s
|
||||
*/
|
||||
virtual void getRevRatesOfProgress(doublereal* revROP) {
|
||||
updateROP();
|
||||
copy(m_kdata->m_ropr.begin(), m_kdata->m_ropr.end(), revROP);
|
||||
}
|
||||
|
||||
/**
|
||||
* Net rates of progress. Return the net (forward - reverse)
|
||||
* rates of progress in array netROP, which must be
|
||||
* dimensioned at least as large as the total number of
|
||||
* reactions.
|
||||
* Units are kmol/m2/s
|
||||
*/
|
||||
virtual void getNetRatesOfProgress(doublereal* netROP) {
|
||||
updateROP();
|
||||
copy(m_kdata->m_ropnet.begin(), m_kdata->m_ropnet.end(), netROP);
|
||||
}
|
||||
|
||||
/**
|
||||
* Equilibrium constants. Return the equilibrium constants of
|
||||
* the reactions in concentration units in array kc, which
|
||||
* must be dimensioned at least as large as the total number
|
||||
* of reactions.
|
||||
*/
|
||||
virtual void getEquilibriumConstants(doublereal* kc);
|
||||
|
||||
|
||||
//@}
|
||||
/**
|
||||
* @name Species Production Rates
|
||||
*/
|
||||
//@{
|
||||
|
||||
/**
|
||||
* Species creation rates [kmol/m^2/s]. Return the species
|
||||
* creation rates in array cdot, which must be
|
||||
* dimensioned at least as large as the total number of
|
||||
* species in all phases of the kinetics
|
||||
* model
|
||||
*
|
||||
*/
|
||||
virtual void getCreationRates(doublereal* cdot) {
|
||||
updateROP();
|
||||
fill(cdot, cdot + m_kk, 0.0);
|
||||
m_revProductStoich.incrementSpecies(
|
||||
m_kdata->m_ropf.begin(), cdot);
|
||||
m_irrevProductStoich.incrementSpecies(
|
||||
m_kdata->m_ropf.begin(), cdot);
|
||||
m_reactantStoich.incrementSpecies(
|
||||
m_kdata->m_ropr.begin(), cdot);
|
||||
}
|
||||
|
||||
/**
|
||||
* Species destruction rates [kmol/m^2/s]. Return the species
|
||||
* destruction rates in array ddot, which must be
|
||||
* dimensioned at least as large as the total number of
|
||||
* species in all phases of the kinetics
|
||||
* model
|
||||
*
|
||||
*/
|
||||
virtual void getDestructionRates(doublereal* ddot) {
|
||||
updateROP();
|
||||
fill(ddot, ddot + m_kk, 0.0);
|
||||
m_revProductStoich.incrementSpecies(
|
||||
m_kdata->m_ropr.begin(), ddot);
|
||||
m_reactantStoich.incrementSpecies(
|
||||
m_kdata->m_ropf.begin(), ddot);
|
||||
}
|
||||
|
||||
/**
|
||||
* Species net production rates [kmol/m^2/s]. Return the species
|
||||
* net production rates (creation - destruction) in array
|
||||
* wdot, which must be dimensioned at least as large as the
|
||||
* total number of species in all phases of the kinetics
|
||||
* model
|
||||
*/
|
||||
virtual void getNetProductionRates(doublereal* net) {
|
||||
updateROP();
|
||||
fill(net, net + m_kk, 0.0);
|
||||
m_revProductStoich.incrementSpecies(
|
||||
m_kdata->m_ropnet.begin(), net);
|
||||
m_irrevProductStoich.incrementSpecies(
|
||||
m_kdata->m_ropnet.begin(), net);
|
||||
m_reactantStoich.decrementSpecies(
|
||||
m_kdata->m_ropnet.begin(), net);
|
||||
}
|
||||
|
||||
//@}
|
||||
/**
|
||||
* @name Reaction Mechanism Informational Query Routines
|
||||
*/
|
||||
//@{
|
||||
|
||||
/**
|
||||
* Stoichiometric coefficient of species k as a reactant in
|
||||
* reaction i.
|
||||
*/
|
||||
virtual doublereal reactantStoichCoeff(int k, int i) const {
|
||||
return m_rrxn[k][i];
|
||||
}
|
||||
|
||||
/**
|
||||
* Stoichiometric coefficient of species k as a product in
|
||||
* reaction i.
|
||||
*/
|
||||
virtual doublereal productStoichCoeff(int k, int i) const {
|
||||
return m_prxn[k][i];
|
||||
}
|
||||
|
||||
/**
|
||||
* Flag specifying the type of reaction. The legal values and
|
||||
* their meaning are specific to the particular kinetics
|
||||
* manager.
|
||||
*/
|
||||
virtual int reactionType(int i) const {
|
||||
return m_index[i].first;
|
||||
}
|
||||
|
||||
/**
|
||||
* True if reaction i has been declared to be reversible. If
|
||||
* isReversible(i) is false, then the reverse rate of progress
|
||||
* for reaction i is always zero.
|
||||
*/
|
||||
virtual bool isReversible(int i) {
|
||||
if (find(m_revindex.begin(), m_revindex.end(), i)
|
||||
< m_revindex.end()) return true;
|
||||
else return false;
|
||||
}
|
||||
|
||||
/**
|
||||
* Return a string representing the reaction.
|
||||
*/
|
||||
virtual string reactionString(int i) const {
|
||||
return m_rxneqn[i];
|
||||
}
|
||||
|
||||
//@}
|
||||
/**
|
||||
* @name Reaction Mechanism Construction
|
||||
*/
|
||||
//@{
|
||||
|
||||
/**
|
||||
* Prepare the class for the addition of reactions. This function
|
||||
* must be called after instantiation of the class, but before
|
||||
* any reactions are actually added to the mechanism.
|
||||
* This function calculates m_kk the number of species in all
|
||||
* phases participating in the reaction mechanism. We don't know
|
||||
* m_kk previously, before all phases have been added.
|
||||
*/
|
||||
virtual void init();
|
||||
|
||||
/**
|
||||
* Add a single reaction to the mechanism.
|
||||
*/
|
||||
virtual void addReaction(const ReactionData& r);
|
||||
|
||||
/**
|
||||
* Finish adding reactions and prepare for use. This function
|
||||
* must be called after all reactions are entered into the mechanism
|
||||
* and before the mechanism is used to calculate reaction rates.
|
||||
*/
|
||||
virtual void finalize();
|
||||
virtual bool ready() const;
|
||||
|
||||
|
||||
void updateROP();
|
||||
|
||||
|
||||
const vector<grouplist_t>& reactantGroups(int i)
|
||||
{ return m_rgroups[i]; }
|
||||
const vector<grouplist_t>& productGroups(int i)
|
||||
{ return m_pgroups[i]; }
|
||||
|
||||
void _update_rates_T();
|
||||
void _update_rates_phi();
|
||||
void _update_rates_C();
|
||||
void checkPartialEquil();
|
||||
|
||||
protected:
|
||||
/**
|
||||
* m_kk here is the number of species in all of the phases
|
||||
* that participate in the kinetics mechanism.
|
||||
*/
|
||||
int m_kk;
|
||||
|
||||
Rate1<SurfaceArrhenius> m_rates;
|
||||
//Rate1<Arrhenius> m_rates;
|
||||
bool m_redo_rates;
|
||||
|
||||
/**
|
||||
* Vector of information about reactions in the
|
||||
* mechanism.
|
||||
* The key is the reaction index (0 < i < m_ii).
|
||||
* The first pair is the reactionType of the reaction.
|
||||
* The second pair is ...
|
||||
*/
|
||||
mutable map<int, pair<int, int> > m_index;
|
||||
|
||||
vector<int> m_irrev;
|
||||
|
||||
StoichManagerN m_reactantStoich;
|
||||
StoichManagerN m_revProductStoich;
|
||||
StoichManagerN m_irrevProductStoich;
|
||||
|
||||
StoichManagerN m_globalReactantStoich;
|
||||
|
||||
int m_nirrev;
|
||||
|
||||
/**
|
||||
* Number of reversible reactions in the mechanism
|
||||
*/
|
||||
int m_nrev;
|
||||
|
||||
map<int, vector<grouplist_t> > m_rgroups;
|
||||
map<int, vector<grouplist_t> > m_pgroups;
|
||||
|
||||
vector<int> m_rxntype;
|
||||
|
||||
mutable vector<map<int, doublereal> > m_rrxn;
|
||||
mutable vector<map<int, doublereal> > m_prxn;
|
||||
|
||||
vector_int m_revindex;
|
||||
vector<string> m_rxneqn;
|
||||
|
||||
/**
|
||||
* Temporary data storage used in calculating the rates of
|
||||
* of reactions.
|
||||
*/
|
||||
EdgeKineticsData* m_kdata;
|
||||
|
||||
/**
|
||||
* An array of generalized concentrations
|
||||
* \f$ C_k \f$ that are defined such that \f$ a_k = C_k /
|
||||
* C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration/
|
||||
* These generalized concentrations are used
|
||||
* by this kinetics manager class to compute the forward and
|
||||
* reverse rates of elementary reactions. The "units" for the
|
||||
* concentrations of each phase depend upon the implementation
|
||||
* of kinetics within that phase.
|
||||
* The order of the species within the vector is based on
|
||||
* the order of listed ThermoPhase objects in the class, and the
|
||||
* order of the species within each ThermoPhase class.
|
||||
*/
|
||||
vector_fp m_conc;
|
||||
|
||||
vector_fp m_mu0;
|
||||
vector_fp m_phi;
|
||||
vector_fp m_pot;
|
||||
vector_fp m_rwork;
|
||||
vector_fp m_E;
|
||||
|
||||
private:
|
||||
|
||||
int reactionNumber(){ return m_ii;}
|
||||
void addElementaryReaction(const ReactionData& r);
|
||||
void addGlobalReaction(const ReactionData& r);
|
||||
void installReagents(const ReactionData& r);
|
||||
|
||||
void installGroups(int irxn, const vector<grouplist_t>& r,
|
||||
const vector<grouplist_t>& p);
|
||||
void updateKc();
|
||||
|
||||
void registerReaction(int rxnNumber, int type, int loc) {
|
||||
m_index[rxnNumber] = pair<int, int>(type, loc);
|
||||
}
|
||||
void applyButlerVolmerCorrection(doublereal* kf);
|
||||
bool m_finalized;
|
||||
};
|
||||
}
|
||||
|
||||
#endif
|
||||
40
Cantera/src/EdgePhase.h
Normal file
40
Cantera/src/EdgePhase.h
Normal file
|
|
@ -0,0 +1,40 @@
|
|||
/**
|
||||
*
|
||||
* @file EdgePhase.h
|
||||
*
|
||||
*/
|
||||
|
||||
/* $Author$
|
||||
* $Date$
|
||||
* $Revision$
|
||||
*
|
||||
* Copyright 2002 California Institute of Technology
|
||||
*
|
||||
*/
|
||||
|
||||
|
||||
#ifndef CT_EDGEPHASE_H
|
||||
#define CT_EDGEPHASE_H
|
||||
|
||||
#include "mix_defs.h"
|
||||
#include "ThermoPhase.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
class EdgePhase : public SurfPhase {
|
||||
|
||||
public:
|
||||
|
||||
EdgePhase(doublereal n0 = 0.0);
|
||||
virtual ~EdgePhase() {}
|
||||
virtual int eosType() const { return cEdge; }
|
||||
|
||||
};
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Loading…
Add table
Reference in a new issue