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

[Kinetics] Always use temperature of "reacting" phase

Browse source 
This is always the lowest-dimensional phase, e.g. surface or edge for
heterogeneous systems.
This commit is contained in:
Ray Speth 2019-02-26 18:21:41 -05:00
parent 5601ee9067
commit 6c82b61c3d
5 changed files with 20 additions and 35 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
src/kinetics/ImplicitSurfChem.cpp
View file

@ -180,7 +180,7 @@ void ImplicitSurfChem::solvePseudoSteadyStateProblem(int ifuncOverride,
// 2) Temperature and pressure
getConcSpecies(m_concSpecies.data());
InterfaceKinetics* ik = m_vecKinPtrs[0];
ThermoPhase& tp = ik->thermo(0);
ThermoPhase& tp = ik->thermo(ik->reactionPhaseIndex());
doublereal TKelvin = tp.temperature();
doublereal PGas = tp.pressure();

17
src/kinetics/InterfaceKinetics.cpp
View file

@ -129,7 +129,7 @@ void InterfaceKinetics::updateKc()
* and m_mu0_Kc[]
*/
updateMu0();
doublereal rrt = 1.0 / thermo(0).RT();
doublereal rrt = 1.0 / thermo(reactionPhaseIndex()).RT();
// compute Delta mu^0 for all reversible reactions
getRevReactionDelta(m_mu0_Kc.data(), m_rkcn.data());
@ -159,7 +159,8 @@ void InterfaceKinetics::updateMu0()
thermo(n).getStandardChemPotentials(m_mu0.data() + m_start[n]);
for (size_t k = 0; k < thermo(n).nSpecies(); k++) {
m_mu0_Kc[ik] = m_mu0[ik] + Faraday * m_phi[n] * thermo(n).charge(k);
m_mu0_Kc[ik] -= thermo(0).RT() * thermo(n).logStandardConc(k);
m_mu0_Kc[ik] -= thermo(reactionPhaseIndex()).RT()
* thermo(n).logStandardConc(k);
ik++;
}
}
@ -168,7 +169,7 @@ void InterfaceKinetics::updateMu0()
void InterfaceKinetics::getEquilibriumConstants(doublereal* kc)
{
updateMu0();
doublereal rrt = 1.0 / thermo(0).RT();
doublereal rrt = 1.0 / thermo(reactionPhaseIndex()).RT();
std::fill(kc, kc + nReactions(), 0.0);
getReactionDelta(m_mu0_Kc.data(), kc);
for (size_t i = 0; i < nReactions(); i++) {
@ -239,7 +240,7 @@ void InterfaceKinetics::applyVoltageKfwdCorrection(doublereal* const kf)
if (m_ctrxn_BVform[i] == 0) {
double eamod = m_beta[i] * deltaElectricEnergy_[irxn];
if (eamod != 0.0) {
kf[irxn] *= exp(-eamod/thermo(0).RT());
kf[irxn] *= exp(-eamod/thermo(reactionPhaseIndex()).RT());
}
}
}
@ -265,7 +266,8 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con
// come out of this calculation.
if (m_ctrxn_BVform[i] == 0) {
// Calculate the term and modify the forward reaction
double tmp = exp(- m_beta[i] * m_deltaG0[irxn] / thermo(0).RT());
double tmp = exp(- m_beta[i] * m_deltaG0[irxn]
/ thermo(reactionPhaseIndex()).RT());
tmp *= 1.0 / m_ProdStanConcReac[irxn] / Faraday;
kfwd[irxn] *= tmp;
}
@ -280,7 +282,8 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con
// Calculate the term and modify the forward reaction rate
// constant so that it's in the exchange current density
// formulation format
double tmp = exp(m_beta[i] * m_deltaG0[irxn] * thermo(0).RT());
double tmp = exp(m_beta[i] * m_deltaG0[irxn]
* thermo(reactionPhaseIndex()).RT());
tmp *= Faraday * m_ProdStanConcReac[irxn];
kfwd[irxn] *= tmp;
}
@ -468,7 +471,7 @@ void InterfaceKinetics::getDeltaSSEnthalpy(doublereal* deltaH)
thermo(n).getEnthalpy_RT(m_grt.data() + m_start[n]);
}
for (size_t k = 0; k < m_kk; k++) {
m_grt[k] *= thermo(0).RT();
m_grt[k] *= thermo(reactionPhaseIndex()).RT();
}
// Use the stoichiometric manager to find deltaH for each reaction.

4
test_problems/surfSolverTest/haca2.cti
View file

@ -15,7 +15,7 @@ units(length = 'cm', quantity = 'mol', act_energy = 'kcal/mol')
ideal_gas(name = 'gas',
elements = 'O H C N Ar',
species = 'gri30: H N2 CH3 CH4 C2H2 H2 OH H2O CO O2',
initial_state = state(temperature = 1400.0,
initial_state = state(temperature = 1100.0,
pressure = OneAtm,
mole_fractions = 'H:0.01, N2:0.8899, H2:0.04, CH4:0.01 C2H2:0.01 \
OH:0.0001 H2O:0.04 O2:0.001'))
@ -53,7 +53,7 @@ ideal_interface(name = 'soot_interface',
reactions = 'all',
phases = 'gas soot',
site_density = (3.8E-9, 'mol/cm2'),
initial_state = state(temperature= 1000.0,
initial_state = state(temperature= 1400.0,
coverages = 'Csoot-*:0.1, Csoot-H:0.9'))
# HKM -> Note, thermo from the following source:

17
test_problems/surfSolverTest/surfaceSolver.cpp
View file

@ -371,19 +371,10 @@ int main(int argc, char** argv)
/* Now Tweak the inputs and do a quick calculation */
/****************************************************************************/
/*
* Set the Gas State:
* -> note that the states are set in the XML files too
*/
/*
* Set the Gas State:
* -> note that the states are set in the XML files too
*/
pres = gasTP->pressure();
double temp = gasTP->temperature();
pres = surfPhaseTP->pressure();
double temp = surfPhaseTP->temperature();
temp += 95;
gasTP->setState_TP(temp, pres);
surfPhaseTP->setState_TP(temp, pres);
iKin_ptr->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);
@ -399,7 +390,7 @@ int main(int argc, char** argv)
/*****************************************************************************/
/* Now Don't Tweak the inputs at all */
/****************************************************************************/
gasTP->setState_TP(temp, pres);
surfPhaseTP->setState_TP(temp, pres);
iKin_ptr->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);

15
test_problems/surfSolverTest/surfaceSolver2.cpp
View file

@ -399,19 +399,10 @@ int main(int argc, char** argv)
/* Now Tweak the inputs and do a quick calculation */
/****************************************************************************/
/*
* Set the Gas State:
* -> note that the states are set in the XML files too
*/
/*
* Set the Gas State:
* -> note that the states are set in the XML files too
*/
pres = gasTP->pressure();
double temp = gasTP->temperature();
pres = surfPhaseTP->pressure();
double temp = surfPhaseTP->temperature();
temp += 95;
gasTP->setState_TP(temp, pres);
surfPhaseTP->setState_TP(temp, pres);
surfaceProb->solvePseudoSteadyStateProblem();
iKin_ptr->getNetProductionRates(src);