diff --git a/doc/sphinx/cython/zerodim.rst b/doc/sphinx/cython/zerodim.rst index 301041160a..48e2cac000 100644 --- a/doc/sphinx/cython/zerodim.rst +++ b/doc/sphinx/cython/zerodim.rst @@ -42,6 +42,10 @@ Reactor ^^^^^^^ .. autoclass:: Reactor +MoleReactor +^^^^^^^^^^^ +.. autoclass:: MoleReactor(contents=None, *, name=None, energy='on') + IdealGasReactor ^^^^^^^^^^^^^^^ .. autoclass:: IdealGasReactor(contents=None, *, name=None, energy='on') @@ -54,6 +58,10 @@ ConstPressureReactor ^^^^^^^^^^^^^^^^^^^^ .. autoclass:: ConstPressureReactor(contents=None, *, name=None, energy='on') +ConstPressureMoleReactor +^^^^^^^^^^^^^^^^^^^^^^^^ +.. autoclass:: ConstPressureMoleReactor(contents=None, *, name=None, energy='on') + IdealGasConstPressureReactor ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. autoclass:: IdealGasConstPressureReactor(contents=None, *, name=None, energy='on') @@ -82,6 +90,22 @@ ExtensibleIdealGasConstPressureReactor ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. autoclass:: ExtensibleIdealGasConstPressureReactor(contents=None, *, name=None, energy='on') +ExtensibleMoleReactor +^^^^^^^^^^^^^^^^^^^^^ +.. autoclass:: ExtensibleMoleReactor(contents=None, *, name=None, energy='on') + +ExtensibleIdealGasMoleReactor +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +.. autoclass:: ExtensibleIdealGasMoleReactor(contents=None, *, name=None, energy='on') + +ExtensibleConstPressureMoleReactor +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +.. autoclass:: ExtensibleConstPressureMoleReactor(contents=None, *, name=None, energy='on') + +ExtensibleIdealGasConstPressureMoleReactor +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +.. autoclass:: ExtensibleIdealGasConstPressureMoleReactor(contents=None, *, name=None, energy='on') + Walls ----- diff --git a/include/cantera/zeroD/ConstPressureMoleReactor.h b/include/cantera/zeroD/ConstPressureMoleReactor.h new file mode 100644 index 0000000000..e7fb04730a --- /dev/null +++ b/include/cantera/zeroD/ConstPressureMoleReactor.h @@ -0,0 +1,45 @@ +//! @file ConstPressureMoleReactor.h + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#ifndef CT_CONSTPRESSMOLE_REACTOR_H +#define CT_CONSTPRESSMOLE_REACTOR_H + +#include "cantera/zeroD/MoleReactor.h" + +namespace Cantera +{ + +/*! + * ConstPressureMoleReactor is a class for constant-pressure reactors + * which use a state of moles. + */ +class ConstPressureMoleReactor : public MoleReactor +{ +public: + ConstPressureMoleReactor() {} + + virtual std::string type() const { + return "ConstPressureMoleReactor"; + }; + + virtual size_t componentIndex(const std::string& nm) const; + + virtual std::string componentName(size_t k); + + virtual void getState(double* y); + + virtual void initialize(double t0 = 0.0); + + virtual void eval(double t, double* LHS, double* RHS); + + virtual void updateState(double* y); + +protected: + const size_t m_sidx = 1; +}; + +} + +#endif diff --git a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h index 190f4c1b02..945e0cd8e3 100644 --- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h +++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h @@ -6,7 +6,7 @@ #ifndef CT_IDEALGASCONSTPRESSMOLE_REACTOR_H #define CT_IDEALGASCONSTPRESSMOLE_REACTOR_H -#include "cantera/zeroD/MoleReactor.h" +#include "cantera/zeroD/ConstPressureMoleReactor.h" namespace Cantera { @@ -15,7 +15,7 @@ namespace Cantera * IdealGasConstPressureMoleReactor is a class for ideal gas constant-pressure reactors * which use a state of moles. */ -class IdealGasConstPressureMoleReactor : public MoleReactor +class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor { public: IdealGasConstPressureMoleReactor() {} @@ -47,8 +47,6 @@ class IdealGasConstPressureMoleReactor : public MoleReactor protected: vector_fp m_hk; //!< Species molar enthalpies - - const size_t m_sidx = 1; }; } diff --git a/include/cantera/zeroD/MoleReactor.h b/include/cantera/zeroD/MoleReactor.h index e76e489f8c..4f33a495f1 100644 --- a/include/cantera/zeroD/MoleReactor.h +++ b/include/cantera/zeroD/MoleReactor.h @@ -13,8 +13,7 @@ namespace Cantera /*! * MoleReactor is meant to serve the same purpose as the reactor class but with a state - * vector composed of moles. It is currently not functional and only serves as a base - * class for other reactors. + * vector composed of moles. It also serves as the base class for other mole reactors. */ class MoleReactor : public Reactor { @@ -27,19 +26,25 @@ class MoleReactor : public Reactor virtual void initialize(double t0 = 0.0); - virtual void eval(double t, double* LHS, double* RHS) { - throw NotImplementedError("MoleReactor::eval()"); - } + virtual void getState(double* y); - size_t componentIndex(const std::string& nm) const { - throw NotImplementedError("MoleReactor::componentIndex"); - } + virtual void updateState(double* y); - std::string componentName(size_t k) { - throw NotImplementedError("MoleReactor::componentName"); - } + virtual void eval(double t, double* LHS, double* RHS); + + size_t componentIndex(const std::string& nm) const; + + std::string componentName(size_t k); protected: + //! Get moles of the system from mass fractions stored by thermo object + //! @param y vector for moles to be put into + virtual void getMoles(double* y); + + //! Set internal mass variable based on moles given + //! @param y vector of moles of the system + virtual void setMassFromMoles(double* y); + virtual void evalSurfaces(double* LHS, double* RHS, double* sdot); virtual void updateSurfaceState(double* y); diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd index 782e80d06d..2d0f6325ea 100644 --- a/interfaces/cython/cantera/reactor.pxd +++ b/interfaces/cython/cantera/reactor.pxd @@ -193,12 +193,18 @@ cdef class Reactor(ReactorBase): cdef CxxReactor* reactor cdef object _kinetics +cdef class MoleReactor(Reactor): + pass + cdef class Reservoir(ReactorBase): pass cdef class ConstPressureReactor(Reactor): pass +cdef class ConstPressureMoleReactor(Reactor): + pass + cdef class IdealGasReactor(Reactor): pass diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index 45772f5c93..1bf4e1f88a 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -369,6 +369,14 @@ cdef class Reactor(ReactorBase): limit = -1. self.reactor.setAdvanceLimit(stringify(name), limit) +cdef class MoleReactor(Reactor): + """ + A homogeneous zero-dimensional reactor with a mole based state vector. By default, + they are closed (no inlets or outlets), have fixed volume, and have adiabatic, + chemically-inert walls. These properties may all be changed by adding + appropriate components such as `Wall`, `MassFlowController` and `Valve`. + """ + reactor_type = "MoleReactor" cdef class Reservoir(ReactorBase): """ @@ -386,6 +394,13 @@ cdef class ConstPressureReactor(Reactor): """ reactor_type = "ConstPressureReactor" +cdef class ConstPressureMoleReactor(Reactor): + """A homogeneous, constant pressure, zero-dimensional reactor with a mole based + state vector. The volume of the reactor changes as a function of time in order to + keep the pressure constant. + """ + reactor_type = "ConstPressureMoleReactor" + cdef class IdealGasReactor(Reactor): """ A constant volume, zero-dimensional reactor for ideal gas mixtures. """ @@ -613,6 +628,38 @@ cdef class ExtensibleIdealGasConstPressureReactor(ExtensibleReactor): reactor_type = "ExtensibleIdealGasConstPressureReactor" +cdef class ExtensibleMoleReactor(ExtensibleReactor): + """ + A variant of `ExtensibleReactor` where the base behavior corresponds to the + `MoleReactor` class. + """ + reactor_type = "ExtensibleMoleReactor" + + +cdef class ExtensibleIdealGasMoleReactor(ExtensibleReactor): + """ + A variant of `ExtensibleReactor` where the base behavior corresponds to the + `IdealGasMoleReactor` class. + """ + reactor_type = "ExtensibleIdealGasMoleReactor" + + +cdef class ExtensibleConstPressureMoleReactor(ExtensibleReactor): + """ + A variant of `ExtensibleReactor` where the base behavior corresponds to the + `ConstPressureMoleReactor` class. + """ + reactor_type = "ExtensibleConstPressureMoleReactor" + + +cdef class ExtensibleIdealGasConstPressureMoleReactor(ExtensibleReactor): + """ + A variant of `ExtensibleReactor` where the base behavior corresponds to the + `IdealGasConstPressureMoleReactor` class. + """ + reactor_type = "ExtensibleIdealGasConstPressureMoleReactor" + + cdef class ReactorSurface: """ Represents a surface in contact with the contents of a reactor. diff --git a/src/zeroD/ConstPressureMoleReactor.cpp b/src/zeroD/ConstPressureMoleReactor.cpp new file mode 100644 index 0000000000..917d8b4d32 --- /dev/null +++ b/src/zeroD/ConstPressureMoleReactor.cpp @@ -0,0 +1,146 @@ +//! @file ConstPressureMoleReactor.cpp A constant pressure +//! zero-dimensional reactor with moles as the state + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#include "cantera/zeroD/Wall.h" +#include "cantera/zeroD/FlowDevice.h" +#include "cantera/zeroD/ConstPressureMoleReactor.h" +#include "cantera/base/utilities.h" +#include "cantera/thermo/SurfPhase.h" +#include "cantera/kinetics/Kinetics.h" + +using namespace std; + +namespace Cantera +{ + +void ConstPressureMoleReactor::getState(double* y) +{ + if (m_thermo == 0) { + throw CanteraError("ConstPressureMoleReactor::getState", + "Error: reactor is empty."); + } + m_thermo->restoreState(m_state); + // set mass to be used in getMoles function + m_mass = m_thermo->density() * m_vol; + // set the first array element to enthalpy + y[0] = m_thermo->enthalpy_mass() * m_thermo->density() * m_vol; + // get moles of species in remaining state + getMoles(y + m_sidx); + // set the remaining components to the surface species moles on the walls + getSurfaceInitialConditions(y+m_nsp+m_sidx); +} + +void ConstPressureMoleReactor::initialize(double t0) +{ + MoleReactor::initialize(t0); + m_nv -= 1; // const pressure system loses 1 more variable from MoleReactor +} + +void ConstPressureMoleReactor::updateState(double* y) +{ + // the components of y are: [0] the enthalpy, [1...K+1) are the + // moles of each species, and [K+1...] are the moles of surface + // species on each wall. + setMassFromMoles(y + m_sidx); + m_thermo->setMolesNoTruncate(y + m_sidx); + if (m_energy) { + m_thermo->setState_HP(y[0] / m_mass, m_pressure); + } else { + m_thermo->setPressure(m_pressure); + } + m_vol = m_mass / m_thermo->density(); + updateConnected(false); + updateSurfaceState(y + m_nsp + m_sidx); +} + +void ConstPressureMoleReactor::eval(double time, double* LHS, double* RHS) +{ + double* dndt = RHS + m_sidx; // kmol per s + + evalWalls(time); + + m_thermo->restoreState(m_state); + + const vector_fp& imw = m_thermo->inverseMolecularWeights(); + + if (m_chem) { + m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" + } + + // evaluate reactor surfaces + evalSurfaces(LHS + m_nsp + m_sidx, RHS + m_nsp + m_sidx, m_sdot.data()); + + // external heat transfer + double dHdt = m_Qdot; + + for (size_t n = 0; n < m_nsp; n++) { + // production in gas phase and from surfaces + dndt[n] = m_wdot[n] * m_vol + m_sdot[n]; + } + + // add terms for outlets + for (auto outlet : m_outlet) { + // determine enthalpy contribution + dHdt -= outlet->massFlowRate() * m_enthalpy; + // flow of species into system and dilution by other species + for (size_t n = 0; n < m_nsp; n++) { + dndt[n] -= outlet->outletSpeciesMassFlowRate(n) * imw[n]; + } + } + + // add terms for inlets + for (auto inlet : m_inlet) { + // enthalpy contribution from inlets + dHdt += inlet->enthalpy_mass() * inlet->massFlowRate(); + // flow of species into system and dilution by other species + for (size_t n = 0; n < m_nsp; n++) { + dndt[n] += inlet->outletSpeciesMassFlowRate(n) * imw[n]; + } + } + + if (m_energy) { + RHS[0] = dHdt; + } else { + RHS[0] = 0.0; + } +} + +size_t ConstPressureMoleReactor::componentIndex(const string& nm) const +{ + size_t k = speciesIndex(nm); + if (k != npos) { + return k + m_sidx; + } else if (nm == "enthalpy") { + return 0; + } else { + return npos; + } +} + +std::string ConstPressureMoleReactor::componentName(size_t k) { + if (k == 0) { + return "enthalpy"; + } else if (k >= m_sidx && k < neq()) { + k -= m_sidx; + if (k < m_thermo->nSpecies()) { + return m_thermo->speciesName(k); + } else { + k -= m_thermo->nSpecies(); + } + for (auto& S : m_surfaces) { + ThermoPhase* th = S->thermo(); + if (k < th->nSpecies()) { + return th->speciesName(k); + } else { + k -= th->nSpecies(); + } + } + } + throw CanteraError("ConstPressureMoleReactor::componentName", + "Index is out of bounds."); +} + +} diff --git a/src/zeroD/IdealGasConstPressureMoleReactor.cpp b/src/zeroD/IdealGasConstPressureMoleReactor.cpp index 172dc7fffe..e44e8529ee 100644 --- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp +++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp @@ -25,7 +25,7 @@ void IdealGasConstPressureMoleReactor::setThermoMgr(ThermoPhase& thermo) throw CanteraError("IdealGasConstPressureMoleReactor::setThermoMgr", "Incompatible phase type provided"); } - MoleReactor::setThermoMgr(thermo); + ConstPressureMoleReactor::setThermoMgr(thermo); } void IdealGasConstPressureMoleReactor::getState(double* y) @@ -39,21 +39,15 @@ void IdealGasConstPressureMoleReactor::getState(double* y) m_mass = m_thermo->density() * m_vol; // set the first component to the temperature y[0] = m_thermo->temperature(); - // Use inverse molecular weights - const double* Y = m_thermo->massFractions(); - const vector_fp& imw = m_thermo->inverseMolecularWeights(); - double *ys = y + m_sidx; - for (size_t i = 0; i < m_nsp; i++) { - ys[i] = m_mass * imw[i] * Y[i]; - } + // get moles of species in remaining state + getMoles(y + m_sidx); // set the remaining components to the surface species moles on the walls getSurfaceInitialConditions(y + m_nsp + m_sidx); } void IdealGasConstPressureMoleReactor::initialize(double t0) { - MoleReactor::initialize(t0); - m_nv -= 1; // const pressure system loses 1 more variable from MoleReactor + ConstPressureMoleReactor::initialize(t0); m_hk.resize(m_nsp, 0.0); } @@ -62,14 +56,9 @@ void IdealGasConstPressureMoleReactor::updateState(double* y) // the components of y are: [0] the temperature, [1...K+1) are the // moles of each species, and [K+1...] are the moles of surface // species on each wall. + setMassFromMoles(y + m_sidx); m_thermo->setMolesNoTruncate(y + m_sidx); m_thermo->setState_TP(y[0], m_pressure); - // calculate mass from moles - const vector_fp& mw = m_thermo->molecularWeights(); - m_mass = 0; - for (size_t i = 0; i < m_nsp; i++) { - m_mass += y[i + m_sidx] * mw[i]; - } m_vol = m_mass / m_thermo->density(); updateConnected(false); updateSurfaceState(y + m_nsp + m_sidx); diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp index ce2f65e977..a553f5aedb 100644 --- a/src/zeroD/IdealGasMoleReactor.cpp +++ b/src/zeroD/IdealGasMoleReactor.cpp @@ -45,13 +45,8 @@ void IdealGasMoleReactor::getState(double* y) // set the second component to the volume y[1] = m_vol; - // use inverse molecular weights - const double* Y = m_thermo->massFractions(); - const vector_fp& imw = m_thermo->inverseMolecularWeights(); - double* ys = y + m_sidx; - for (size_t i = 0; i < m_nsp; i++) { - ys[i] = m_mass * imw[i] * Y[i]; - } + // get moles of species in remaining state + getMoles(y + m_sidx); // set the remaining components to the surface species moles on // the walls getSurfaceInitialConditions(y + m_nsp + m_sidx); @@ -75,25 +70,9 @@ std::string IdealGasMoleReactor::componentName(size_t k) { if (k == 0) { return "temperature"; - } else if (k == 1) { - return "volume"; - } else if (k >= m_sidx && k < neq()) { - k -= m_sidx; - if (k < m_thermo->nSpecies()) { - return m_thermo->speciesName(k); - } else { - k -= m_thermo->nSpecies(); - } - for (auto& S : m_surfaces) { - ThermoPhase* th = S->thermo(); - if (k < th->nSpecies()) { - return th->speciesName(k); - } else { - k -= th->nSpecies(); - } - } + } else { + return MoleReactor::componentName(k); } - throw IndexError("IdealGasMoleReactor::componentName", "components", k, neq() - 1); } void IdealGasMoleReactor::initialize(double t0) @@ -107,12 +86,7 @@ void IdealGasMoleReactor::updateState(double* y) // the components of y are: [0] the temperature, [1] the volume, [2...K+1) are the // moles of each species, and [K+1...] are the moles of surface // species on each wall. - const vector_fp& mw = m_thermo->molecularWeights(); - // calculate mass from moles - m_mass = 0; - for (size_t i = 0; i < m_nv - m_sidx; i++) { - m_mass += y[i + m_sidx] * mw[i]; - } + setMassFromMoles(y + m_sidx); m_vol = y[1]; // set state m_thermo->setMolesNoTruncate(y + m_sidx); diff --git a/src/zeroD/MoleReactor.cpp b/src/zeroD/MoleReactor.cpp index a06fe25946..9d4bedacd3 100644 --- a/src/zeroD/MoleReactor.cpp +++ b/src/zeroD/MoleReactor.cpp @@ -5,11 +5,15 @@ // at https://cantera.org/license.txt for license and copyright information. #include "cantera/zeroD/MoleReactor.h" -#include "cantera/thermo/SurfPhase.h" +#include "cantera/zeroD/FlowDevice.h" #include "cantera/zeroD/ReactorSurface.h" +#include "cantera/thermo/SurfPhase.h" #include "cantera/kinetics/Kinetics.h" +#include "cantera/base/utilities.h" +#include using namespace std; +namespace bmt = boost::math::tools; namespace Cantera { @@ -80,4 +84,186 @@ void MoleReactor::evalSurfaces(double* LHS, double* RHS, double* sdot) } } +void MoleReactor::getMoles(double* y) +{ + // Use inverse molecular weights to convert to moles + const double* Y = m_thermo->massFractions(); + const vector_fp& imw = m_thermo->inverseMolecularWeights(); + for (size_t i = 0; i < m_nsp; i++) { + y[i] = m_mass * imw[i] * Y[i]; + } +} + +void MoleReactor::setMassFromMoles(double* y) +{ + const vector_fp& mw = m_thermo->molecularWeights(); + // calculate mass from moles + m_mass = 0; + for (size_t i = 0; i < m_nsp; i++) { + m_mass += y[i] * mw[i]; + } +} + +void MoleReactor::getState(double* y) +{ + if (m_thermo == 0) { + throw CanteraError("MoleReactor::getState", + "Error: reactor is empty."); + } + m_thermo->restoreState(m_state); + // set the first component to the internal energy + m_mass = m_thermo->density() * m_vol; + y[0] = m_thermo->intEnergy_mass() * m_mass; + // set the second component to the total volume + y[1] = m_vol; + // set components y+2 ... y+K+2 to the moles of each species + getMoles(y + m_sidx); + // set the remaining components to the surface species + // moles on walls + getSurfaceInitialConditions(y+m_nsp+m_sidx); +} + + +void MoleReactor::updateState(double* y) +{ + // The components of y are [0] total internal energy, [1] the total volume, and + // [2...K+3] are the moles of each species, and [K+3...] are the moles + // of surface species on each wall. + setMassFromMoles(y + m_sidx); + m_vol = y[1]; + m_thermo->setMolesNoTruncate(y + m_sidx); + if (m_energy) { + double U = y[0]; + // Residual function: error in internal energy as a function of T + auto u_err = [this, U](double T) { + m_thermo->setState_TR(T, m_mass / m_vol); + return m_thermo->intEnergy_mass() * m_mass - U; + }; + double T = m_thermo->temperature(); + boost::uintmax_t maxiter = 100; + std::pair TT; + try { + TT = bmt::bracket_and_solve_root( + u_err, T, 1.2, true, bmt::eps_tolerance(48), maxiter); + } catch (std::exception&) { + // Try full-range bisection if bracketing fails (for example, near + // temperature limits for the phase's equation of state) + try { + TT = bmt::bisect(u_err, m_thermo->minTemp(), m_thermo->maxTemp(), + bmt::eps_tolerance(48), maxiter); + } catch (std::exception& err2) { + // Set m_thermo back to a reasonable state if root finding fails + m_thermo->setState_TR(T, m_mass / m_vol); + throw CanteraError("MoleReactor::updateState", + "{}\nat U = {}, rho = {}", err2.what(), U, m_mass / m_vol); + } + } + if (fabs(TT.first - TT.second) > 1e-7*TT.first) { + throw CanteraError("MoleReactor::updateState", "root finding failed"); + } + m_thermo->setState_TR(TT.second, m_mass / m_vol); + } else { + m_thermo->setDensity(m_mass / m_vol); + } + updateConnected(true); + updateSurfaceState(y + m_nsp + m_sidx); +} + +void MoleReactor::eval(double time, double* LHS, double* RHS) +{ + double* dndt = RHS + m_sidx; // moles per time + + evalWalls(time); + m_thermo->restoreState(m_state); + + evalSurfaces(LHS + m_nsp + m_sidx, RHS + m_nsp + m_sidx, m_sdot.data()); + // inverse molecular weights for conversion + const vector_fp& imw = m_thermo->inverseMolecularWeights(); + // volume equation + RHS[1] = m_vdot; + + if (m_chem) { + m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" + } + + // Energy equation. + // \f[ + // \dot U = -P\dot V + A \dot q + \dot m_{in} h_{in} - \dot m_{out} h. + // \f] + if (m_energy) { + RHS[0] = - m_thermo->pressure() * m_vdot + m_Qdot; + } else { + RHS[0] = 0.0; + } + + for (size_t k = 0; k < m_nsp; k++) { + // production in gas phase and from surfaces + dndt[k] = m_wdot[k] * m_vol + m_sdot[k]; + } + + // add terms for outlets + for (auto outlet : m_outlet) { + // flow of species into system and dilution by other species + for (size_t n = 0; n < m_nsp; n++) { + dndt[n] -= outlet->outletSpeciesMassFlowRate(n) * imw[n]; + } + // energy update based on mass flow + double mdot = outlet->massFlowRate(); + if (m_energy) { + RHS[0] -= mdot * m_enthalpy; + } + } + + // add terms for inlets + for (auto inlet : m_inlet) { + double mdot = inlet->massFlowRate(); + for (size_t n = 0; n < m_nsp; n++) { + // flow of species into system and dilution by other species + dndt[n] += inlet->outletSpeciesMassFlowRate(n) * imw[n]; + } + if (m_energy) { + RHS[0] += mdot * inlet->enthalpy_mass(); + } + } +} + + +size_t MoleReactor::componentIndex(const string& nm) const +{ + size_t k = speciesIndex(nm); + if (k != npos) { + return k + m_sidx; + } else if (nm == "int_energy") { + return 0; + } else if (nm == "volume") { + return 1; + } else { + return npos; + } +} + +std::string MoleReactor::componentName(size_t k) { + if (k == 0) { + return "int_energy"; + } else if (k == 1) { + return "volume"; + } else if (k >= m_sidx && k < neq()) { + k -= m_sidx; + if (k < m_thermo->nSpecies()) { + return m_thermo->speciesName(k); + } else { + k -= m_thermo->nSpecies(); + } + for (auto& S : m_surfaces) { + ThermoPhase* th = S->thermo(); + if (k < th->nSpecies()) { + return th->speciesName(k); + } else { + k -= th->nSpecies(); + } + } + } + throw CanteraError("MoleReactor::componentName", "Index is out of bounds."); +} + } diff --git a/src/zeroD/ReactorFactory.cpp b/src/zeroD/ReactorFactory.cpp index 29ae1baa1b..a618203756 100644 --- a/src/zeroD/ReactorFactory.cpp +++ b/src/zeroD/ReactorFactory.cpp @@ -6,8 +6,10 @@ #include "cantera/zeroD/ReactorFactory.h" #include "cantera/zeroD/Reservoir.h" #include "cantera/zeroD/Reactor.h" +#include "cantera/zeroD/MoleReactor.h" #include "cantera/zeroD/FlowReactor.h" #include "cantera/zeroD/ConstPressureReactor.h" +#include "cantera/zeroD/ConstPressureMoleReactor.h" #include "cantera/zeroD/IdealGasReactor.h" #include "cantera/zeroD/IdealGasMoleReactor.h" #include "cantera/zeroD/IdealGasConstPressureReactor.h" @@ -38,8 +40,19 @@ ReactorFactory::ReactorFactory() []() { return new ReactorDelegator(); }); reg("ExtensibleIdealGasConstPressureReactor", []() { return new ReactorDelegator(); }); - reg("IdealGasConstPressureMoleReactor", []() { return new IdealGasConstPressureMoleReactor(); }); + reg("ExtensibleMoleReactor", + []() { return new ReactorDelegator(); }); + reg("ExtensibleConstPressureMoleReactor", + []() { return new ReactorDelegator(); }); + reg("ExtensibleIdealGasMoleReactor", + []() { return new ReactorDelegator(); }); + reg("ExtensibleIdealGasConstPressureMoleReactor", + []() { return new ReactorDelegator(); }); + reg("IdealGasConstPressureMoleReactor", []() { return new + IdealGasConstPressureMoleReactor(); }); reg("IdealGasMoleReactor", []() { return new IdealGasMoleReactor(); }); + reg("ConstPressureMoleReactor", []() { return new ConstPressureMoleReactor(); }); + reg("MoleReactor", []() { return new MoleReactor(); }); } ReactorBase* ReactorFactory::newReactor(const std::string& reactorType) diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index c409c2645a..25feeb899c 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -89,6 +89,7 @@ def test_component_index(self): def test_component_names(self): self.make_reactors(n_reactors=2) + self.net.initialize() N = self.net.n_vars // 2 for i in range(N): self.assertEqual(self.r1.component_index(self.r1.component_name(i)), i) @@ -258,10 +259,9 @@ def integrate(atol, rtol): n_baseline = integrate(1e-10, 1e-20) n_rtol = integrate(5e-7, 1e-20) - n_atol = integrate(1e-10, 1e-6) - - self.assertTrue(n_baseline > n_rtol) - self.assertTrue(n_baseline > n_atol) + n_atol = integrate(1e-10, 1e-5) + assert n_baseline > n_rtol + assert n_baseline > n_atol def test_advance_limits(self): P0 = 10 * ct.one_atm @@ -823,6 +823,56 @@ def test_bad_kwarg(self): with self.assertRaises(TypeError): r1 = self.reactorClass(foobar=3.14) +class TestMoleReactor(TestReactor): + reactorClass = ct.MoleReactor + + def test_mole_reactor_surface_chem(self): + model = "ptcombust.yaml" + # initial conditions + T0 = 1500 + P0 = ct.one_atm + equiv_ratio = 1 + surf_area = 0.527 + fuel = "CH4" + air = "O2:1.0, N2:3.773" + # reactor 1 + gas1 = ct.Solution(model, transport_model=None) + gas1.TP = T0, P0 + gas1.set_equivalence_ratio(equiv_ratio, fuel, air) + r1 = self.reactorClass(gas1) + # comparison reactor + gas2 = ct.Solution(model, transport_model=None) + gas2.TP = T0, P0 + gas2.set_equivalence_ratio(equiv_ratio, fuel, air) + if "ConstPressure" in r1.type: + r2 = ct.ConstPressureReactor(gas2) + else: + r2 = ct.Reactor(gas2) + # surf 1 + surf1 = ct.Interface(model, "Pt_surf", [gas1]) + surf1.TP = T0, P0 + surf1.coverages = {"PT(S)":1} + rsurf1 = ct.ReactorSurface(surf1, r1, A=surf_area) + # surf 2 + surf2 = ct.Interface(model, "Pt_surf", [gas2]) + surf2.TP = T0, P0 + surf2.coverages = {"PT(S)":1} + rsurf2 = ct.ReactorSurface(surf2, r2, A=surf_area) + # reactor network setup + net1 = ct.ReactorNet([r1,]) + net2 = ct.ReactorNet([r2,]) + net1.rtol = net2.rtol = 1e-9 + net1.atol = net2.atol = 1e-18 + # steady state occurs at ~0.002 seconds + for i in np.linspace(0, 0.0025, 50)[1:]: + net1.advance(i) + net2.advance(i) + self.assertArrayNear(r1.thermo.Y, r2.thermo.Y, rtol=5e-4, atol=1e-6) + self.assertNear(r1.T, r2.T, rtol=5e-5) + self.assertNear(r1.thermo.P, r2.thermo.P, rtol=1e-6) + self.assertArrayNear(rsurf1.coverages, rsurf2.coverages, rtol=1e-4, + atol=1e-8) + class TestIdealGasReactor(TestReactor): reactorClass = ct.IdealGasReactor @@ -1070,12 +1120,21 @@ def test_with_surface_reactions(self): self.net1.rtol = self.net2.rtol = 1e-9 self.integrate(surf=True) +class TestConstPressureMoleReactor(TestConstPressureReactor): + """ + The constant pressure reactor should give the same results as + as a regular "Reactor" with a wall with a very high expansion rate + coefficient. + """ + reactorClass = ct.ConstPressureMoleReactor + test_mole_reactor_surface_chem = TestMoleReactor.test_mole_reactor_surface_chem + class TestIdealGasConstPressureReactor(TestConstPressureReactor): reactorClass = ct.IdealGasConstPressureReactor -class TestIdealGasConstPressureMoleReactor(TestIdealGasConstPressureReactor): +class TestIdealGasConstPressureMoleReactor(TestConstPressureMoleReactor): reactorClass = ct.IdealGasConstPressureMoleReactor def create_reactors(self, **kwargs): @@ -1097,7 +1156,7 @@ def test_component_index(self): with pytest.raises(ct.CanteraError): super().test_component_index() -class TestIdealGasMoleReactor(TestReactor): +class TestIdealGasMoleReactor(TestMoleReactor): reactorClass = ct.IdealGasMoleReactor def test_adaptive_precon_integration(self):