From e10fa97af025eb6371f575a067388cc512152e71 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Wed, 15 Dec 2010 02:54:08 +0000 Subject: [PATCH] Added the phaseIsStable idea --- Cantera/src/kinetics/InterfaceKinetics.cpp | 30 ++++++++ Cantera/src/kinetics/InterfaceKinetics.h | 85 ++++++++++++++++++++-- 2 files changed, 108 insertions(+), 7 deletions(-) diff --git a/Cantera/src/kinetics/InterfaceKinetics.cpp b/Cantera/src/kinetics/InterfaceKinetics.cpp index 6b6f60c9d..3ad664d24 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.cpp +++ b/Cantera/src/kinetics/InterfaceKinetics.cpp @@ -92,6 +92,7 @@ namespace Cantera { m_has_exchange_current_density_formulation(false), m_phaseExistsCheck(false), m_phaseExists(0), + m_phaseIsStable(0), m_rxnPhaseIsReactant(0), m_rxnPhaseIsProduct(0), m_ioFlag(0) @@ -140,6 +141,7 @@ namespace Cantera { m_has_exchange_current_density_formulation(false), m_phaseExistsCheck(false), m_phaseExists(0), + m_phaseIsStable(0), m_rxnPhaseIsReactant(0), m_rxnPhaseIsProduct(0), m_ioFlag(0) @@ -208,6 +210,7 @@ namespace Cantera { m_has_exchange_current_density_formulation = right.m_has_exchange_current_density_formulation; m_phaseExistsCheck = right.m_phaseExistsCheck; m_phaseExists = right.m_phaseExists; + m_phaseIsStable = right.m_phaseIsStable; m_rxnPhaseIsReactant.resize(m_ii, 0); m_rxnPhaseIsProduct.resize(m_ii, 0); @@ -753,6 +756,12 @@ namespace Cantera { } } } + if (m_rxnPhaseIsReactant[j][p]) { + if (! m_phaseIsStable[p]) { + ropnet[j] = 0.0; + ropr[j] = ropf[j]; + } + } } } else if ((ropf[j] > ropr[j]) && (ropf[j] > 0.0)) { for (int p = 0; p < nPhases(); p++) { @@ -772,6 +781,12 @@ namespace Cantera { } } } + if (m_rxnPhaseIsProduct[j][p]) { + if (! m_phaseIsStable[p]) { + ropnet[j] = 0.0; + ropf[j] = ropr[j]; + } + } } } } @@ -1229,6 +1244,7 @@ namespace Cantera { void InterfaceKinetics::addPhase(thermo_t &thermo) { Kinetics::addPhase(thermo); m_phaseExists.push_back(true); + m_phaseIsStable.push_back(true); } //================================================================================================ /** @@ -1358,14 +1374,28 @@ namespace Cantera { if (!m_phaseExists[iphase]) { m_phaseExistsCheck--; m_phaseExists[iphase] = true; + m_phaseIsStable[iphase] = true; } } else { if (m_phaseExists[iphase]) { m_phaseExistsCheck++; m_phaseExists[iphase] = false; + m_phaseIsStable[iphase] = false; } } } + //================================================================================================ + + void InterfaceKinetics::setPhaseStability(const int iphase, const bool isStable) { + if (iphase < 0 || iphase >= (int) m_thermo.size()) { + throw CanteraError("InterfaceKinetics:setPhaseStability", "out of bounds"); + } + if (isStable) { + m_phaseIsStable[iphase] = true; + } else { + m_phaseIsStable[iphase] = false; + } + } //================================================================================================ void EdgeKinetics::finalize() { m_rwork.resize(nReactions()); diff --git a/Cantera/src/kinetics/InterfaceKinetics.h b/Cantera/src/kinetics/InterfaceKinetics.h index 9c8c4ac5c..82f30bd4c 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.h +++ b/Cantera/src/kinetics/InterfaceKinetics.h @@ -73,12 +73,38 @@ namespace Cantera { }; - /// - /// A kinetics manager for heterogeneous reaction mechanisms. The - /// reactions are assumed to occur at a 2D interface between two - /// 3D phases. - /// - /// @ingroup chemkinetics + + //! A kinetics manager for heterogeneous reaction mechanisms. The + //! reactions are assumed to occur at a 2D interface between two 3D phases. + /*! + * + * There are some important additions to the behavior of the kinetics class due to the + * presence of multiple phases and a heterogeneous interface. If a reactant phase + * doesn't exists, i.e., has a mole number of zero, a heterogeneous reaction can not + * proceed from reactants to products. Note it could perhaps proceed from products to + * reactants if all of the product phases exist. + * + * In order to make the determination of whether a phase exists or not actually involves + * the specification of additional information to the kinetics object., which heretofore + * has only had access to intrinsic field information about the phases (i.e., temperature + * pressure, and mole fraction). + * + * The extrinsic specification of whether a phase exists or not must be specified on top of the + * intrinsic calculation of the reaction rate. This routine carries a set of + * booleans indicating whether a phase in the heterogeneous mechanism exists or not. + * + * Additionally, the routine carries a set of booleans around indicating whether a product + * phase is stable or not. If a phase is not thermodynamically stable, it may be the case that + * a particular reaction in a heterogeneous mechanism will create a product species in the + * unstable phase. However, other reactions in the mechanism will destruct that species. + * This may cause oscillations in the formation of the unstable phase from time step to time + * step within a ODE solver, in practice. In order to avoid this situation, a set of + * booleans is tracked which sets the stability of a phase. If a phase is deemed to be unstable, + * then species in that phase will not be allowed to be birthed by the kinetics operator. + * Nonexistent phases are deemed to be unstable by default, but this can be changed. + * + * @ingroup chemkinetics + */ class InterfaceKinetics : public Kinetics { public: @@ -557,13 +583,34 @@ namespace Cantera { /*! * Tell the kinetics object whether a phase in the object exists. * This is actually an extrinsic specification that must be carried out on top of the - * intrinsic calculation of the reaction rate + * intrinsic calculation of the reaction rate. + * The routine will also flip the IsStable boolean within the kinetics object as well. * * @param iphase Index of the phase. This is the order within the internal thermo vector object * @param exists Boolean indicating whether the phase exists or not */ void setPhaseExistence(const int iphase, const bool exists); + + //! Set the stability of a phase in the reaction object + /*! + * Tell the kinetics object whether a phase in the object is stable. Species in an unstable phase + * will not be allowed to have a positive rate of formation from this kinetics object. + * This is actually an extrinsic specification that must be carried out on top of the + * intrinsic calculation of the reaction rate. + * + * While conceptually not needed since kinetics is consistent with thermo when taken as a whole, + * in practice it has found to be very useful to turn off the creation of phases which shouldn't + * be forming. Typically thais can reduce the oscillations in phase formation and destruction + * which are observed. + * + * @param iphase Index of the phase. This is the order within the internal thermo vector object + * @param exists Boolean indicating whether the phase exists or not + */ + void setPhaseStability(const int iphase, const bool isStable); + + + protected: //! Temporary work vector of length m_kk @@ -795,7 +842,31 @@ namespace Cantera { */ std::vector m_phaseExists; + //! Vector of booleans indicating whether phases are stable or not + /*! + * Vector of booleans indicating whether a phase is stable or not + * under the current conditions. + * We use this to set the ROP's so that unphysical things don't happen + * + * length = number of phases in the object + * By default all phases are stable + */ + std::vector m_phaseIsStable; + + //! Vector of vector of booleans indicating whether a phase participates in a + //! reaction as a reactant + /*! + * m_rxnPhaseIsReactant[j][p] indicates whether a species in phase p + * participates in reaction j as a reactant. + */ std::vector m_rxnPhaseIsReactant; + + //! Vector of vector of booleans indicating whether a phase participates in a + //! reaction as a product + /*! + * m_rxnPhaseIsReactant[j][p] indicates whether a species in phase p + * participates in reaction j as a product. + */ std::vector m_rxnPhaseIsProduct; int m_ioFlag;