From 2c170d8d9de571a69d820041d44fdb606bbea9ff Mon Sep 17 00:00:00 2001 From: Bang-Shiuh Chen Date: Sat, 30 Sep 2023 14:28:15 -0400 Subject: [PATCH] [reaction, interface] add ElectronCollisionPlasmaRate --- .../kinetics/ElectronCollisionPlasmaRate.h | 121 ++++++++++++++ include/cantera/thermo/PlasmaPhase.h | 5 + interfaces/cython/cantera/reaction.pxd | 9 + interfaces/cython/cantera/reaction.pyx | 63 ++++++- interfaces/cython/cantera/thermo.pxd | 1 + interfaces/cython/cantera/thermo.pyx | 9 +- src/kinetics/ElectronCollisionPlasmaRate.cpp | 155 ++++++++++++++++++ src/kinetics/ReactionRateFactory.cpp | 6 + 8 files changed, 367 insertions(+), 2 deletions(-) create mode 100644 include/cantera/kinetics/ElectronCollisionPlasmaRate.h create mode 100644 src/kinetics/ElectronCollisionPlasmaRate.cpp diff --git a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h new file mode 100644 index 000000000..c4cac709a --- /dev/null +++ b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h @@ -0,0 +1,121 @@ +//! @file ElectronCollisionPlasmaRate.h Header for plasma reaction rates parameterized +//! by electron collision cross section and electron energy distribution. + +// 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_ELECTRONCOLLISIONPLASMARATE_H +#define CT_ELECTRONCOLLISIONPLASMARATE_H + +#include "cantera/thermo/PlasmaPhase.h" +#include "cantera/kinetics/ReactionData.h" +#include "ReactionRate.h" +#include "MultiRate.h" + +namespace Cantera +{ + +//! Data container holding shared data specific to ElectronCollisionPlasmaRate +/** + * The data container `ElectronCollisionPlasmaData` holds precalculated data common to + * all `ElectronCollisionPlasmaRate` objects. + */ +struct ElectronCollisionPlasmaData : public ReactionData +{ + ElectronCollisionPlasmaData(); + + virtual bool update(const ThermoPhase& phase, const Kinetics& kin) override; + using ReactionData::update; + + virtual void invalidateCache() override { + ReactionData::invalidateCache(); + energyLevels.resize(0); + distribution.resize(0); + } + + vector energyLevels; //!< electron energy levels + vector distribution; //!< electron energy distribution +}; + + +//! Electron collision plasma reaction rate type +/*! + * The electron collision plasma reaction rate uses the electron collision + * data and the electron energy distribution to calculate the reaction rate. + * + * @f[ + * k_f = \gamma \int_0^{\infty} \epsilon \sigma F_0 d\epsilon + * @f] + * + * where @f$ \gamma = (2e/m)^{1/2} @f$ is a constant, @f$ \epsilon @f$ is + * the electron energy, @f$ \sigma @f$ is the reaction collision cross + * section, and @f$ F_0 @f$ is the normalized electron energy distribution + * function. + * + * References: + * + * [1] Hagelaar and Pitchford @cite hagelaar2005. + * + * @ingroup plasma + */ +class ElectronCollisionPlasmaRate : public ReactionRate +{ +public: + ElectronCollisionPlasmaRate() = default; + + ElectronCollisionPlasmaRate(const AnyMap& node, + const UnitStack& rate_units={}) + : ElectronCollisionPlasmaRate() + { + setParameters(node, rate_units); + } + + virtual void setParameters(const AnyMap& node, const UnitStack& units) override; + + virtual void getParameters(AnyMap& node) const override; + + unique_ptr newMultiRate() const override { + return unique_ptr( + new MultiRate); + } + + virtual const std::string type() const override { + return "electron-collision-plasma"; + } + + virtual void setContext(const Reaction& rxn, const Kinetics& kin) override; + + //! Evaluate reaction rate + /*! + * @param shared_data data shared by all reactions of a given type + */ + double evalFromStruct(const ElectronCollisionPlasmaData& shared_data) const; + + //! Evaluate derivative of reaction rate with respect to temperature + //! divided by reaction rate + //! @param shared_data data shared by all reactions of a given type + double ddTScaledFromStruct(const ElectronCollisionPlasmaData& shared_data) const { + throw NotImplementedError("ElectronCollisionPlasmaRate::ddTScaledFromStruct"); + } + + //! The value of #m_energyLevels [eV] + const vector& energyLevels() const { + return m_energyLevels; + } + + //! The value of #m_crossSections [m2] + const vector& crossSections() const { + return m_crossSections; + } + +private: + //! electron energy levels [eV] + vector m_energyLevels; + + //! collision cross sections [m2] at #m_energyLevels + vector m_crossSections; +}; + +} + +#endif diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 239b794ea..b5020fb1e 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -249,6 +249,11 @@ public: void setParameters(const AnyMap& phaseNode, const AnyMap& rootNode=AnyMap()) override; + //! Electron species name + string electronSpeciesName() const { + return speciesName(m_electronSpeciesIndex); + } + protected: void updateThermo() const override; diff --git a/interfaces/cython/cantera/reaction.pxd b/interfaces/cython/cantera/reaction.pxd index 2e53e2390..d4180034b 100644 --- a/interfaces/cython/cantera/reaction.pxd +++ b/interfaces/cython/cantera/reaction.pxd @@ -73,6 +73,11 @@ cdef extern from "cantera/kinetics/TwoTempPlasmaRate.h" namespace "Cantera": CxxTwoTempPlasmaRate(double, double, double, double) double activationElectronEnergy() +cdef extern from "cantera/kinetics/ElectronCollisionPlasmaRate.h" namespace "Cantera": + cdef cppclass CxxElectronCollisionPlasmaRate "Cantera::ElectronCollisionPlasmaRate" (CxxReactionRate): + CxxElectronCollisionPlasmaRate(CxxAnyMap) except +translate_exception + vector[double]& energyLevels() + vector[double]& crossSections() cdef extern from "cantera/base/Array.h" namespace "Cantera": cdef cppclass CxxArray2D "Cantera::Array2D": @@ -251,6 +256,10 @@ cdef class ReactionRate: cdef class ArrheniusRateBase(ReactionRate): cdef CxxArrheniusBase* base +cdef class ElectronCollisionPlasmaRate(ReactionRate): + cdef CxxElectronCollisionPlasmaRate* base + cdef set_cxx_object(self) + cdef class FalloffRate(ReactionRate): cdef CxxFalloffRate* falloff cdef set_cxx_object(self) diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index b6a0af837..c49ee5a64 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -370,6 +370,67 @@ cdef class TwoTempPlasmaRate(ArrheniusRateBase): return self.cxx_object().activationElectronEnergy() +cdef class ElectronCollisionPlasmaRate(ReactionRate): + r""" + A reaction rate coefficient which depends on electron collision cross section + and electron energy distribution + + """ + _reaction_rate_type = "electron-collision-plasma" + + def __cinit__(self, energy_levels=None, cross_sections=None, + input_data=None, init=True): + if init: + if isinstance(input_data, dict): + self._from_dict(input_data) + elif energy_levels is not None and cross_sections is not None: + self._from_parameters(energy_levels, cross_sections) + elif input_data is None: + self._from_dict({}) + else: + raise TypeError("Invalid input parameters") + self.set_cxx_object() + + def _from_dict(self, dict input_data): + self._rate.reset( + new CxxElectronCollisionPlasmaRate(py_to_anymap(input_data, hyphenize=True)) + ) + + def _from_parameters(self, energy_levels, cross_sections): + # check length + if len(energy_levels) != len(cross_sections): + raise ValueError('Length of energy levels and ' + 'cross sections are different') + cdef np.ndarray[np.double_t, ndim=1] data_energy_levels = \ + np.ascontiguousarray(energy_levels, dtype=np.double) + cdef np.ndarray[np.double_t, ndim=1] data_cross_sections = \ + np.ascontiguousarray(cross_sections, dtype=np.double) + input_data = {'type': 'electron-collision-plasma', + 'energy-levels': data_energy_levels, + 'cross-sections': data_cross_sections} + self._from_dict(input_data) + + cdef set_cxx_object(self): + self.rate = self._rate.get() + self.base = self.rate + + property energy_levels: + """ + The energy levels [eV]. Each level corresponds to a cross section + of `cross_sections`. + """ + def __get__(self): + return np.fromiter(self.base.energyLevels(), np.double) + + property cross_sections: + """ + The cross sections [m2]. Each cross section corresponds to a energy + level of `energy_levels`. + """ + def __get__(self): + return np.fromiter(self.base.crossSections(), np.double) + + cdef class FalloffRate(ReactionRate): """ Base class for parameterizations used to describe the fall-off in reaction rates @@ -1463,7 +1524,7 @@ cdef class Reaction: Directories on Cantera's input file path will be searched for the specified file. """ - root = AnyMapFromYamlFile(stringify(filename)) + root = AnyMapFromYamlFile(stringify(str(filename))) cxx_reactions = CxxGetReactions(root[stringify(section)], deref(kinetics.kinetics)) return [Reaction.wrap(r) for r in cxx_reactions] diff --git a/interfaces/cython/cantera/thermo.pxd b/interfaces/cython/cantera/thermo.pxd index a6e4291b5..fe1db7e54 100644 --- a/interfaces/cython/cantera/thermo.pxd +++ b/interfaces/cython/cantera/thermo.pxd @@ -209,6 +209,7 @@ cdef extern from "cantera/thermo/PlasmaPhase.h": double meanElectronEnergy() size_t nElectronEnergyLevels() double electronPressure() + string electronSpeciesName() cdef extern from "cantera/cython/thermo_utils.h": diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx index 577dfb87c..3f9179192 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/interfaces/cython/cantera/thermo.pyx @@ -1865,7 +1865,7 @@ cdef class ThermoPhase(_SolutionBase): def __set__(self, distribution_type): if not self._enable_plasma: raise ThermoModelMethodError(self.thermo_model) - self.plasma.setElectronEnergyDistributionType(distribution_type) + self.plasma.setElectronEnergyDistributionType(stringify(distribution_type)) property mean_electron_energy: """ Mean electron energy [eV] """ @@ -1900,6 +1900,13 @@ cdef class ThermoPhase(_SolutionBase): raise ThermoModelMethodError(self.thermo_model) self.plasma.enableNormalizeElectronEnergyDist(enable) + property electron_species_name: + """ Electron species name """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return pystr(self.plasma.electronSpeciesName()) + cdef class InterfacePhase(ThermoPhase): """ A class representing a surface, edge phase """ diff --git a/src/kinetics/ElectronCollisionPlasmaRate.cpp b/src/kinetics/ElectronCollisionPlasmaRate.cpp new file mode 100644 index 000000000..1a6813dca --- /dev/null +++ b/src/kinetics/ElectronCollisionPlasmaRate.cpp @@ -0,0 +1,155 @@ +//! @file ElectronCollisionPlasmaRate.cpp + +// 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/kinetics/ElectronCollisionPlasmaRate.h" +#include "cantera/kinetics/Reaction.h" +#include "cantera/kinetics/Kinetics.h" +#include "cantera/numerics/funcs.h" + +namespace Cantera +{ + +ElectronCollisionPlasmaData::ElectronCollisionPlasmaData() +{ + energyLevels.assign(1, 0.0); + distribution.assign(1, 0.0); +} + +bool ElectronCollisionPlasmaData::update(const ThermoPhase& phase, + const Kinetics& kin) +{ + const PlasmaPhase& pp = dynamic_cast(phase); + + vector levels(pp.nElectronEnergyLevels()); + pp.getElectronEnergyLevels(levels.data()); + + vector dist(pp.nElectronEnergyLevels()); + pp.getElectronEnergyDistribution(dist.data()); + + bool changed = false; + + if (levels != energyLevels) { + energyLevels = std::move(levels); + changed = true; + } + + if (dist != distribution) { + distribution = std::move(dist); + changed = true; + } + + return changed; +} + +void ElectronCollisionPlasmaRate::setParameters(const AnyMap& node, const UnitStack& rate_units) +{ + ReactionRate::setParameters(node, rate_units); + if (!node.hasKey("energy-levels") && !node.hasKey("cross-sections")) { + return; + } + if (node.hasKey("energy-levels")) { + m_energyLevels = node["energy-levels"].asVector(); + } + if (node.hasKey("cross-sections")) { + m_crossSections = node["cross-sections"].asVector(); + } + if (m_energyLevels.size() != m_crossSections.size()) { + throw CanteraError("ElectronCollisionPlasmaRate::setParameters", + "Energy levels and cross section must have the same length."); + } +} + +void ElectronCollisionPlasmaRate::getParameters(AnyMap& node) const { + node["type"] = type(); + node["energy-levels"] = m_energyLevels; + node["cross-sections"] = m_crossSections; +} + +double ElectronCollisionPlasmaRate::evalFromStruct(const ElectronCollisionPlasmaData& shared_data) const +{ + // Interpolate cross-sections data to the energy levels of + // the electron energy distribution function + vector cross_sections; + for (auto i = 0; i < shared_data.energyLevels.size(); i ++ ) { + cross_sections.push_back(linearInterp(shared_data.energyLevels[i], + m_energyLevels, m_crossSections)); + } + + // Map cross sections to Eigen::ArrayXd + auto cs_array = Eigen::Map( + cross_sections.data(), cross_sections.size() + ); + + // Map energyLevels in Eigen::ArrayXd + auto eps = Eigen::Map( + shared_data.energyLevels.data(), shared_data.energyLevels.size() + ); + + // Map energyLevels in Eigen::ArrayXd + auto distribution = Eigen::Map( + shared_data.distribution.data(), shared_data.distribution.size() + ); + + // unit in kmol/m3/s + return 0.5 * pow(2.0 * ElectronCharge / ElectronMass, 0.5) * Avogadro * + simpson(distribution.cwiseProduct(cs_array), eps.pow(2.0)); +} + +void ElectronCollisionPlasmaRate::setContext(const Reaction& rxn, const Kinetics& kin) +{ + // ElectronCollisionPlasmaReaction is for a non-equilibrium plasma, and the reverse rate + // cannot be calculated from the conventional thermochemistry. + // @todo implement the reversible rate for non-equilibrium plasma + if (rxn.reversible) { + throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input, + "ElectronCollisionPlasmaRate does not support reversible reactions"); + } + + // get electron species name + string electronName; + if (kin.thermo().type() == "plasma") { + electronName = dynamic_cast(kin.thermo()).electronSpeciesName(); + } else { + throw CanteraError("ElectronCollisionPlasmaRate::setContext", + "ElectronCollisionPlasmaRate requires plasma phase"); + } + + // Check rxn.reactantString (rxn.reactants gives reduced reactants) + std::istringstream iss(rxn.reactantString()); + string token; + int nElectron = 0; + int nReactants = 0; + + // Count number of electron and reactants + // Since the reactants are one electron and one molecule, + // no token can be a number. + while (iss >> token) { + if (isdigit(token[0])) { + throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input, + "ElectronCollisionPlasmaRate requires one electron and one molecule as reactants"); + } + if (token == electronName) { + nElectron++; + } + if (token != "+") { + nReactants++; + } + } + + // Number of reactants needs to be two + if (nReactants != 2) { + throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input, + "ElectronCollisionPlasmaRate requires exactly two reactants"); + } + + // Must have only one electron + // @todo add electron-electron collision rate + if (nElectron != 1) { + throw InputFileError("ElectronCollisionPlasmaRate::setContext", rxn.input, + "ElectronCollisionPlasmaRate requires one and only one electron"); + } +} + +} diff --git a/src/kinetics/ReactionRateFactory.cpp b/src/kinetics/ReactionRateFactory.cpp index 55dc46381..3b6b54432 100644 --- a/src/kinetics/ReactionRateFactory.cpp +++ b/src/kinetics/ReactionRateFactory.cpp @@ -12,6 +12,7 @@ #include "cantera/kinetics/Arrhenius.h" #include "cantera/kinetics/ChebyshevRate.h" #include "cantera/kinetics/Custom.h" +#include "cantera/kinetics/ElectronCollisionPlasmaRate.h" #include "cantera/kinetics/Falloff.h" #include "cantera/kinetics/InterfaceRate.h" #include "cantera/kinetics/PlogRate.h" @@ -38,6 +39,11 @@ ReactionRateFactory::ReactionRateFactory() return new TwoTempPlasmaRate(node, rate_units); }); + // ElectronCollisionPlasmaRate evaluator + reg("electron-collision-plasma", [](const AnyMap& node, const UnitStack& rate_units) { + return new ElectronCollisionPlasmaRate(node, rate_units); + }); + // BlowersMaselRate evaluator reg("Blowers-Masel", [](const AnyMap& node, const UnitStack& rate_units) { return new BlowersMaselRate(node, rate_units);