From 0da063cdef0729c66c7ee526e814ad1e5133319b Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 10 Jun 2014 16:21:56 +0000 Subject: [PATCH] [Reactor] Refactor wall / surface chemistry evaluation This creates a single implementation of the calculation of the contributions of walls and surface chemistry to the governing equations for all reactor types. --- include/cantera/zeroD/Reactor.h | 12 +++ src/zeroD/ConstPressureReactor.cpp | 55 ++----------- src/zeroD/IdealGasConstPressureReactor.cpp | 58 ++----------- src/zeroD/IdealGasReactor.cpp | 53 ++---------- src/zeroD/Reactor.cpp | 96 ++++++++++++---------- 5 files changed, 88 insertions(+), 186 deletions(-) diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 6489c4d51..287f15e04 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -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; diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp index c02454047..03f0572d0 100644 --- a/src/zeroD/ConstPressureReactor.cpp +++ b/src/zeroD/ConstPressureReactor.cpp @@ -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(); diff --git a/src/zeroD/IdealGasConstPressureReactor.cpp b/src/zeroD/IdealGasConstPressureReactor.cpp index 7bcf3d2b3..68a83fe6d 100644 --- a/src/zeroD/IdealGasConstPressureReactor.cpp +++ b/src/zeroD/IdealGasConstPressureReactor.cpp @@ -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; } diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp index 91f59010e..6338f0a4f 100644 --- a/src/zeroD/IdealGasReactor.cpp +++ b/src/zeroD/IdealGasReactor.cpp @@ -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; } diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 5e0739327..6190f1674 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -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())