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

[Reactor] Refactor wall / surface chemistry evaluation

Browse source 
This creates a single implementation of the calculation of the contributions of
walls and surface chemistry to the governing equations for all reactor types.
This commit is contained in:
Ray Speth 2014-06-10 16:21:56 +00:00
parent 54c14c3ac1
commit 0da063cdef
5 changed files with 88 additions and 186 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

12
include/cantera/zeroD/Reactor.h
View file

@ -157,6 +157,18 @@ protected:
//! specific reactor implementations.
virtual size_t speciesIndex(const std::string& nm) const;
//! Evaluate terms related to Walls
//! Calculates #m_vdot and #m_Q based on wall movement and heat transfer
//! @param t the current time
virtual void evalWalls(double t);
//! Evaluate terms related to surface reactions
//! Calculates #m_sdot and rate of change in surface species coverages
//! @param t the current time
//! @param[out] ydot array of d(coverage)/dt for surface species
//! @returns Net mass flux from surfaces
virtual double evalSurfaces(double t, double* ydot);
//! Pointer to the homogeneous Kinetics object that handles the reactions
Kinetics* m_kin;

55
src/zeroD/ConstPressureReactor.cpp
View file

@ -117,49 +117,14 @@ void ConstPressureReactor::updateState(doublereal* y)
void ConstPressureReactor::evalEqs(doublereal time, doublereal* y,
doublereal* ydot, doublereal* params)
{
size_t nk;
m_thermo->restoreState(m_state);
Kinetics* kin;
applySensitivity(params);
m_Q = 0.0;
// compute wall terms
doublereal rs0, sum, wallarea;
double dmdt = 0.0; // dm/dt (gas phase)
double* dYdt = ydot + 2;
SurfPhase* surf;
size_t lr, ns, loc = m_nsp+2, surfloc;
fill(m_sdot.begin(), m_sdot.end(), 0.0);
for (size_t i = 0; i < m_nwalls; i++) {
lr = 1 - 2*m_lr[i];
m_Q += lr*m_wall[i]->Q(time);
kin = m_wall[i]->kinetics(m_lr[i]);
surf = m_wall[i]->surface(m_lr[i]);
if (surf && kin) {
rs0 = 1.0/surf->siteDensity();
nk = surf->nSpecies();
sum = 0.0;
surf->setTemperature(m_state[0]);
m_wall[i]->syncCoverages(m_lr[i]);
kin->getNetProductionRates(DATA_PTR(m_work));
ns = kin->surfacePhaseIndex();
surfloc = kin->kineticsSpeciesIndex(0,ns);
for (size_t k = 1; k < nk; k++) {
ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
sum -= ydot[loc + k];
}
ydot[loc] = sum;
loc += nk;
wallarea = m_wall[i]->area();
for (size_t k = 0; k < m_nsp; k++) {
m_sdot[k] += m_work[k]*wallarea;
}
}
}
m_thermo->restoreState(m_state);
applySensitivity(params);
evalWalls(time);
double mdot_surf = evalSurfaces(time, ydot + m_nsp + 2);
dmdt += mdot_surf;
const vector_fp& mw = m_thermo->molecularWeights();
const doublereal* Y = m_thermo->massFractions();
@ -168,22 +133,16 @@ void ConstPressureReactor::evalEqs(doublereal time, doublereal* y,
m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
}
double mdot_surf = 0.0; // net mass flux from surface
for (size_t k = 0; k < m_nsp; k++) {
// production in gas phase and from surfaces
dYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k] / m_mass;
mdot_surf += m_sdot[k] * mw[k];
// dilution by net surface mass flux
dYdt[k] -= Y[k] * mdot_surf / m_mass;
}
dmdt += mdot_surf;
// external heat transfer
double dHdt = - m_Q;
for (size_t n = 0; n < m_nsp; n++) {
// dilution by net surface mass flux
dYdt[n] -= Y[n] * mdot_surf / m_mass;
}
// add terms for open system
if (m_open) {
double enthalpy = m_thermo->enthalpy_mass();

58
src/zeroD/IdealGasConstPressureReactor.cpp
View file

@ -125,51 +125,17 @@ void IdealGasConstPressureReactor::updateState(doublereal* y)
void IdealGasConstPressureReactor::evalEqs(doublereal time, doublereal* y,
doublereal* ydot, doublereal* params)
{
size_t nk;
double dmdt = 0.0; // dm/dt (gas phase)
double mcpdTdt = 0.0; // m * c_p * dT/dt
double* dYdt = ydot + 2;
m_thermo->restoreState(m_state);
applySensitivity(params);
evalWalls(time);
double mdot_surf = evalSurfaces(time, ydot + m_nsp + 2);
dmdt += mdot_surf;
Kinetics* kin;
m_Q = 0.0;
// compute wall terms
doublereal rs0, sum, wallarea;
double mcpdTdt = 0.0; // m * c_p * dT/dt
double dmdt = 0.0; // dm/dt (gas phase)
double* dYdt = ydot + 2;
m_thermo->getPartialMolarEnthalpies(&m_hk[0]);
SurfPhase* surf;
size_t lr, ns, loc = m_nsp+2, surfloc;
fill(m_sdot.begin(), m_sdot.end(), 0.0);
for (size_t i = 0; i < m_nwalls; i++) {
lr = 1 - 2*m_lr[i];
m_Q += lr*m_wall[i]->Q(time);
kin = m_wall[i]->kinetics(m_lr[i]);
surf = m_wall[i]->surface(m_lr[i]);
if (surf && kin) {
rs0 = 1.0/surf->siteDensity();
nk = surf->nSpecies();
sum = 0.0;
surf->setTemperature(m_state[0]);
m_wall[i]->syncCoverages(m_lr[i]);
kin->getNetProductionRates(DATA_PTR(m_work));
ns = kin->surfacePhaseIndex();
surfloc = kin->kineticsSpeciesIndex(0,ns);
for (size_t k = 1; k < nk; k++) {
ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
sum -= ydot[loc + k];
}
ydot[loc] = sum;
loc += nk;
wallarea = m_wall[i]->area();
for (size_t k = 0; k < m_nsp; k++) {
m_sdot[k] += m_work[k]*wallarea;
}
}
}
const vector_fp& mw = m_thermo->molecularWeights();
const doublereal* Y = m_thermo->massFractions();
@ -177,14 +143,6 @@ void IdealGasConstPressureReactor::evalEqs(doublereal time, doublereal* y,
m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
}
double mdot_surf = 0.0; // net mass flux from surface
for (size_t k = 0; k < m_nsp; k++) {
// production in gas phase and from surfaces
dYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k] / m_mass;
mdot_surf += m_sdot[k] * mw[k];
}
dmdt += mdot_surf;
// external heat transfer
mcpdTdt -= m_Q;
@ -192,6 +150,8 @@ void IdealGasConstPressureReactor::evalEqs(doublereal time, doublereal* y,
// heat release from gas phase and surface reations
mcpdTdt -= m_wdot[n] * m_hk[n] * m_vol;
mcpdTdt -= m_sdot[n] * m_hk[n];
// production in gas phase and from surfaces
dYdt[n] = (m_wdot[n] * m_vol + m_sdot[n]) * mw[n] / m_mass;
// dilution by net surface mass flux
dYdt[n] -= Y[n] * mdot_surf / m_mass;
}

53
src/zeroD/IdealGasReactor.cpp
View file

@ -135,50 +135,13 @@ void IdealGasReactor::updateState(doublereal* y)
void IdealGasReactor::evalEqs(doublereal time, doublereal* y,
doublereal* ydot, doublereal* params)
{
m_thermo->restoreState(m_state);
applySensitivity(params);
m_vdot = 0.0;
m_Q = 0.0;
double mcvdTdt = 0.0; // m * c_v * dT/dt
double dmdt = 0.0; // dm/dt (gas phase)
double mcvdTdt = 0.0; // m * c_v * dT/dt
double* dYdt = ydot + 3;
m_thermo->restoreState(m_state);
applySensitivity(params);
m_thermo->getPartialMolarIntEnergies(&m_uk[0]);
// compute wall terms
size_t loc = m_nsp+3;
fill(m_sdot.begin(), m_sdot.end(), 0.0);
for (size_t i = 0; i < m_nwalls; i++) {
int lr = 1 - 2*m_lr[i];
double vdot = lr*m_wall[i]->vdot(time);
m_vdot += vdot;
m_Q += lr*m_wall[i]->Q(time);
Kinetics* kin = m_wall[i]->kinetics(m_lr[i]);
SurfPhase* surf = m_wall[i]->surface(m_lr[i]);
if (surf && kin) {
double rs0 = 1.0/surf->siteDensity();
size_t nk = surf->nSpecies();
double sum = 0.0;
surf->setTemperature(m_state[0]);
m_wall[i]->syncCoverages(m_lr[i]);
kin->getNetProductionRates(DATA_PTR(m_work));
size_t ns = kin->surfacePhaseIndex();
size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
for (size_t k = 1; k < nk; k++) {
ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
sum -= ydot[loc + k];
}
ydot[loc] = sum;
loc += nk;
double wallarea = m_wall[i]->area();
for (size_t k = 0; k < m_nsp; k++) {
m_sdot[k] += m_work[k]*wallarea;
}
}
}
const vector_fp& mw = m_thermo->molecularWeights();
const doublereal* Y = m_thermo->massFractions();
@ -186,12 +149,8 @@ void IdealGasReactor::evalEqs(doublereal time, doublereal* y,
m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
}
double mdot_surf = 0.0; // net mass flux from surfaces
for (size_t k = 0; k < m_nsp; k++) {
// production in gas phase and from surfaces
dYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k] / m_mass;
mdot_surf += m_sdot[k] * mw[k];
}
evalWalls(time);
double mdot_surf = evalSurfaces(time, ydot + m_nsp + 3);
dmdt += mdot_surf;
// compression work and external heat transfer
@ -201,6 +160,8 @@ void IdealGasReactor::evalEqs(doublereal time, doublereal* y,
// heat release from gas phase and surface reations
mcvdTdt -= m_wdot[n] * m_uk[n] * m_vol;
mcvdTdt -= m_sdot[n] * m_uk[n];
// production in gas phase and from surfaces
dYdt[n] = (m_wdot[n] * m_vol + m_sdot[n]) * mw[n] / m_mass;
// dilution by net surface mass flux
dYdt[n] -= Y[n] * mdot_surf / m_mass;
}

96
src/zeroD/Reactor.cpp
View file

@ -177,46 +177,14 @@ void Reactor::updateState(doublereal* y)
void Reactor::evalEqs(doublereal time, doublereal* y,
doublereal* ydot, doublereal* params)
{
m_thermo->restoreState(m_state);
applySensitivity(params);
m_vdot = 0.0;
m_Q = 0.0;
double dmdt = 0.0; // dm/dt (gas phase)
double* dYdt = ydot + 3;
// compute wall terms
size_t loc = m_nsp+3;
fill(m_sdot.begin(), m_sdot.end(), 0.0);
for (size_t i = 0; i < m_nwalls; i++) {
int lr = 1 - 2*m_lr[i];
double vdot = lr*m_wall[i]->vdot(time);
m_vdot += vdot;
m_Q += lr*m_wall[i]->Q(time);
Kinetics* kin = m_wall[i]->kinetics(m_lr[i]);
SurfPhase* surf = m_wall[i]->surface(m_lr[i]);
if (surf && kin) {
double rs0 = 1.0/surf->siteDensity();
size_t nk = surf->nSpecies();
double sum = 0.0;
surf->setTemperature(m_state[0]);
m_wall[i]->syncCoverages(m_lr[i]);
kin->getNetProductionRates(DATA_PTR(m_work));
size_t ns = kin->surfacePhaseIndex();
size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
for (size_t k = 1; k < nk; k++) {
ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
sum -= ydot[loc + k];
}
ydot[loc] = sum;
loc += nk;
double wallarea = m_wall[i]->area();
for (size_t k = 0; k < m_nsp; k++) {
m_sdot[k] += m_work[k]*wallarea;
}
}
}
m_thermo->restoreState(m_state);
applySensitivity(params);
evalWalls(time);
double mdot_surf = evalSurfaces(time, ydot + m_nsp + 3);
dmdt += mdot_surf; // mass added to gas phase from surface reations
// volume equation
ydot[1] = m_vdot;
@ -228,15 +196,9 @@ void Reactor::evalEqs(doublereal time, doublereal* y,
m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
}
double mdot_surf = 0.0; // net mass flux from surfaces
for (size_t k = 0; k < m_nsp; k++) {
// production in gas phase and from surfaces
dYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k] / m_mass;
mdot_surf += m_sdot[k] * mw[k];
}
dmdt += mdot_surf; // mass added to gas phase from surface reations
for (size_t k = 0; k < m_nsp; k++) {
// dilution by net surface mass flux
dYdt[k] -= Y[k] * mdot_surf / m_mass;
}
@ -292,6 +254,54 @@ void Reactor::evalEqs(doublereal time, doublereal* y,
resetSensitivity(params);
}
void Reactor::evalWalls(double t)
{
m_vdot = 0.0;
m_Q = 0.0;
for (size_t i = 0; i < m_nwalls; i++) {
int lr = 1 - 2*m_lr[i];
m_vdot += lr*m_wall[i]->vdot(t);
m_Q += lr*m_wall[i]->Q(t);
}
}
double Reactor::evalSurfaces(double t, double* ydot)
{
const vector_fp& mw = m_thermo->molecularWeights();
fill(m_sdot.begin(), m_sdot.end(), 0.0);
size_t loc = 0; // offset into ydot
double mdot_surf = 0.0; // net mass flux from surface
for (size_t i = 0; i < m_nwalls; i++) {
Kinetics* kin = m_wall[i]->kinetics(m_lr[i]);
SurfPhase* surf = m_wall[i]->surface(m_lr[i]);
if (surf && kin) {
double rs0 = 1.0/surf->siteDensity();
size_t nk = surf->nSpecies();
double sum = 0.0;
surf->setTemperature(m_state[0]);
m_wall[i]->syncCoverages(m_lr[i]);
kin->getNetProductionRates(&m_work[0]);
size_t ns = kin->surfacePhaseIndex();
size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
for (size_t k = 1; k < nk; k++) {
ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
sum -= ydot[loc + k];
}
ydot[loc] = sum;
loc += nk;
double wallarea = m_wall[i]->area();
for (size_t k = 0; k < m_nsp; k++) {
m_sdot[k] += m_work[k]*wallarea;
mdot_surf += m_sdot[k] * mw[k];
}
}
}
return mdot_surf;
}
void Reactor::addSensitivityReaction(size_t rxn)
{
if (rxn >= m_kin->nReactions())