diff --git a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h index c0d9206662..31d9e41d3b 100644 --- a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h +++ b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h @@ -31,6 +31,8 @@ struct ElectronCollisionPlasmaData : public ReactionData ReactionData::invalidateCache(); energyLevels.resize(0); distribution.resize(0); + m_dist_number = -1; + m_level_number = -1; } vector energyLevels; //!< electron energy levels @@ -89,6 +91,7 @@ struct ElectronCollisionPlasmaData : public ReactionData class ElectronCollisionPlasmaRate : public ReactionRate { public: + ElectronCollisionPlasmaRate() = default; ElectronCollisionPlasmaRate(const AnyMap& node, @@ -135,6 +138,14 @@ class ElectronCollisionPlasmaRate : public ReactionRate return m_crossSections; } + //! The value of #m_crossSectionsInterpolated [m2] + const vector& crossSectionInterpolated() const { + return m_crossSectionsInterpolated; + } + + //! Update the value of #m_crossSectionsInterpolated [m2] + void updateInterpolatedCrossSection(const vector&); + private: //! electron energy levels [eV] vector m_energyLevels; diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index 7d81acb342..8eae70e538 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -1407,6 +1407,26 @@ class Kinetics return m_root.lock(); } + //! Register a function to be called if reaction is added. + //! @param id A unique ID corresponding to the object affected by the callback. + //! Typically, this is a pointer to an object that also holds a reference to the + //! Kinetics object. + //! @param callback The callback function to be called after any reaction is added. + //! When the callback becomes invalid (for example, the corresponding object is + //! being deleted, the removeReactionAddedCallback() method must be invoked. + //! @since New in %Cantera 3.1 + void registerReactionAddedCallback(void* id, const function& callback) + { + m_reactionAddedCallbacks[id] = callback; + } + + //! Remove the reaction-changed callback function associated with the specified object. + //! @since New in %Cantera 3.1 + void removeReactionAddedCallback(void* id) + { + m_reactionAddedCallbacks.erase(id); + } + protected: //! Cache for saved calculations within each Kinetics object. ValueCache m_cache; @@ -1530,6 +1550,9 @@ class Kinetics //! reference to Solution std::weak_ptr m_root; + + //! Callback functions that are invoked when the reaction is added. + map> m_reactionAddedCallbacks; }; } diff --git a/include/cantera/kinetics/Reaction.h b/include/cantera/kinetics/Reaction.h index bd148dea15..f64fe458a8 100644 --- a/include/cantera/kinetics/Reaction.h +++ b/include/cantera/kinetics/Reaction.h @@ -162,6 +162,18 @@ class Reaction return bool(m_third_body); } + //! Register a function to be called if setRate is used. + //! @param id A unique ID corresponding to the object affected by the callback. + //! @param callback The callback function to be called after setRate is called. + //! When the callback becomes invalid (for example, the corresponding object is + //! being deleted), the removeSetRateCallback() method must be invoked. + //! @since New in %Cantera 3.2 + void registerSetRateCallback(void* id, const function& callback); + + //! Remove the setRate callback function associated with the specified object. + //! @since New in %Cantera 3.2 + void removeSetRateCallback(void* id); + protected: //! Store the parameters of a Reaction needed to reconstruct an identical //! object using the newReaction(AnyMap&, Kinetics&) function. Does not @@ -183,6 +195,9 @@ class Reaction //! Relative efficiencies of third-body species in enhancing the reaction rate //! (if applicable) shared_ptr m_third_body; + + //! Callback functions that are invoked when the rate is replaced. + map> m_setRateCallbacks; }; diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h index 116edb4be4..3d4b844ace 100644 --- a/include/cantera/thermo/Phase.h +++ b/include/cantera/thermo/Phase.h @@ -18,6 +18,7 @@ namespace Cantera class Solution; class Species; +class Kinetics; //! Class Phase is the base class for phases of matter, managing the species and //! elements in a phase, as well as the independent variables of temperature, diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 8debb8e7a1..aa6075a6a5 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -15,38 +15,63 @@ namespace Cantera { -/** - * Base class for a phase with plasma properties. This class manages the - * plasma properties such as electron energy distribution function (EEDF). - * There are two ways to define the electron distribution and electron - * temperature. The first method uses setElectronTemperature() to set - * the electron temperature which is used to calculate the electron energy - * distribution with isotropic-velocity model. The generalized electron - * energy distribution for isotropic-velocity distribution can be - * expressed as [1,2], +class Reaction; +class ElectronCollisionPlasmaRate; + +//! Base class for handling plasma properties, specifically focusing on the +//! electron energy distribution. +/*! + * This class provides functionality to manage the the electron energy distribution + * using two primary methods for defining the electron distribution and electron + * temperature. + * + * The first method utilizes setElectronTemperature(), which sets the electron + * temperature and calculates the electron energy distribution assuming an + * isotropic-velocity model. Note that all units in PlasmaPhase are in SI, except + * for electron energy, which is measured in volts. + * + * The generalized electron energy distribution for an isotropic-velocity + * distribution (as described by Gudmundsson @cite gudmundsson2001 and Khalilpour + * and Foroutan @cite khalilpour2020) + * is given by: * @f[ * f(\epsilon) = c_1 \frac{\sqrt{\epsilon}}{\epsilon_m^{3/2}} - * \exp(-c_2 (\frac{\epsilon}{\epsilon_m})^x), + * \exp \left(-c_2 \left(\frac{\epsilon}{\epsilon_m}\right)^x \right), * @f] - * where @f$ x = 1 @f$ and @f$ x = 2 @f$ correspond to the Maxwellian and - * Druyvesteyn (default) electron energy distribution, respectively. - * @f$ \epsilon_m = 3/2 T_e @f$ [eV] (mean electron energy). The second - * method uses setDiscretizedElectronEnergyDist() to manually set electron - * energy distribution and calculate electron temperature from mean electron - * energy, which is calculated as [3], + * where @f$ x = 1 @f$ corresponds to a Maxwellian distribution and + * @f$ x = 2 @f$ corresponds to a Druyvesteyn distribution (which is the + * default). Here, @f$ \epsilon_m = \frac{3}{2} T_e @f$ [V] represents the + * mean electron energy. + * + * The total probability distribution integrates to one: * @f[ - * \epsilon_m = \int_0^{\infty} \epsilon^{3/2} f(\epsilon) d\epsilon, + * \int_0^{\infty} f(\epsilon) d\epsilon = 1. * @f] - * which can be calculated using trapezoidal rule, + * According to Hagelaar and Pitchford @cite hagelaar2005, the electron energy + * probability function can be defined as + * @f$ F(\epsilon) = \frac{f(\epsilon)}{\sqrt{\epsilon}} @f$ with units of + * [V@f$^{-3/2}@f$]. The generalized form of the electron energy probability + * function for isotropic-velocity distributions is: * @f[ - * \epsilon_m = \sum_i (\epsilon^{5/2}_{i+1} - \epsilon^{5/2}_i) - * (f(\epsilon_{i+1}) + f(\epsilon_i)) / 2, + * F(\epsilon) = c_1 \frac{1}{\epsilon_m^{3/2}} + * \exp\left(-c_2 \left(\frac{\epsilon}{\epsilon_m}\right)^x\right), * @f] - * where @f$ i @f$ is the index of energy levels. + * and this form is used to model the isotropic electron energy distribution + * in PlasmaPhase. * - * For references, see Gudmundsson @cite gudmundsson2001; Khalilpour and Foroutan - * @cite khalilpour2020; Hagelaar and Pitchford @cite hagelaar2005, and BOLOS - * @cite BOLOS. + * The second method allows for manual definition of the electron energy + * distribution using setDiscretizedElectronEnergyDist(). In this approach, + * the electron temperature is derived from the mean electron energy, + * @f$ \epsilon_m @f$, which can be calculated as follows @cite hagelaar2005 : + * @f[ + * \epsilon_m = \int_0^{\infty} \epsilon^{3/2} F(\epsilon) d\epsilon. + * @f] + * This integral can be approximated using the trapezoidal rule, + * @f[ + * \epsilon_m = \sum_i \left(\epsilon_{i+1}^{5/2} - \epsilon_i^{5/2}\right) + * \frac{F(\epsilon_{i+1}) + F(\epsilon_i)}{2}, + * @f] + * where @f$ i @f$ is the index of discrete energy levels, or Simpson's rule. * * @warning This class is an experimental part of %Cantera and may be * changed or removed without notice. @@ -71,6 +96,8 @@ class PlasmaPhase: public IdealGasPhase */ explicit PlasmaPhase(const string& inputFile="", const string& id=""); + ~PlasmaPhase(); + string type() const override { return "plasma"; } @@ -81,7 +108,9 @@ class PlasmaPhase: public IdealGasPhase //! @param levels The vector of electron energy levels (eV). //! Length: #m_nPoints. //! @param length The length of the @c levels. - void setElectronEnergyLevels(const double* levels, size_t length); + //! @param updateEnergyDist update electron energy distribution + void setElectronEnergyLevels(const double* levels, size_t length, + bool updateEnergyDist=true); //! Get electron energy levels. //! @param levels The vector of electron energy levels (eV). Length: #m_nPoints @@ -99,13 +128,6 @@ class PlasmaPhase: public IdealGasPhase const double* distrb, size_t length); - //! Set discretized electron energy distribution. - //! @param distrb The vector of electron energy distribution. - //! Length: #m_nPoints. - //! @param length The length of the vectors, which equals #m_nPoints. - void setDiscretizedElectronEnergyDist(const double* distrb, - size_t length); - //! Get electron energy distribution. //! @param distrb The vector of electron energy distribution. //! Length: #m_nPoints. @@ -198,6 +220,11 @@ class PlasmaPhase: public IdealGasPhase return m_nPoints; } + //! Number of collisions + size_t nCollisions() const { + return m_collisions.size(); + } + //! Electron Species Index size_t electronSpeciesIndex() const { return m_electronSpeciesIndex; @@ -271,6 +298,19 @@ class PlasmaPhase: public IdealGasPhase return m_levelNum; } + virtual void setSolution(std::weak_ptr soln) override; + + /** + * The elastic power loss (J/s/m³) + * @f[ + * P_k = N_A N_A C_e e \sum_k C_k K_k, + * @f] + * where @f$ C_k @f$ and @f$ C_e @f$ are the concentration (kmol/m³) of the + * target species and electrons, respectively. @f$ K_k @f$ is the elastic + * electron energy loss coefficient (eV-m³/s). + */ + double elasticPowerLoss(); + protected: void updateThermo() const override; @@ -317,6 +357,12 @@ class PlasmaPhase: public IdealGasPhase //! Electron energy distribution norm void normalizeElectronEnergyDistribution(); + //! Update interpolated cross section of a collision + bool updateInterpolatedCrossSection(size_t k); + + //! Update electron energy distribution difference + void updateElectronEnergyDistDifference(); + // Electron energy order in the exponential term double m_isotropicShapeFactor = 2.0; @@ -345,6 +391,39 @@ class PlasmaPhase: public IdealGasPhase //! Flag of normalizing electron energy distribution bool m_do_normalizeElectronEnergyDist = true; + //! Data for initiate reaction + AnyMap m_root; + + //! Electron energy distribution Difference dF/dε (V^-5/2) + Eigen::ArrayXd m_electronEnergyDistDiff; + + //! Elastic electron energy loss coefficients (eV m3/s) + /*! The elastic electron energy loss coefficient for species k is, + * @f[ + * K_k = \frac{2 m_e}{m_k} \sqrt{\frac{2 e}{m_e}} \int_0^{\infty} \sigma_k + * \epsilon^2 \left( F_0 + \frac{k_B T}{e} + * \frac{\partial F_0}{\partial \epsilon} \right) d \epsilon, + * @f] + * where @f$ m_e @f$ [kg] is the electron mass, @f$ \epsilon @f$ [V] is the + * electron energy, @f$ \sigma_k @f$ [m2] is the reaction collision cross section, + * @f$ F_0 @f$ [V^(-3/2)] is the normalized electron energy distribution function. + */ + vector m_elasticElectronEnergyLossCoefficients; + + //! Updates the elastic electron energy loss coefficient for collision index i + /*! Calculates the elastic energy loss coefficient using the current electron + energy distribution and cross sections. + */ + void updateElasticElectronEnergyLossCoefficient(size_t i); + + //! Update elastic electron energy loss coefficients + /*! Used by elasticPowerLoss() and other plasma property calculations that + depends on #m_elasticElectronEnergyLossCoefficients. This function calls + updateInterpolatedCrossSection() before calling + updateElasticElectronEnergyLossCoefficient() + */ + void updateElasticElectronEnergyLossCoefficients(); + private: //! Electron energy distribution change variable. Whenever //! #m_electronEnergyDist changes, this int is incremented. @@ -353,6 +432,30 @@ class PlasmaPhase: public IdealGasPhase //! Electron energy level change variable. Whenever //! #m_electronEnergyLevels changes, this int is incremented. int m_levelNum = -1; + + //! The list of shared pointers of plasma collision reactions + vector> m_collisions; + + //! The list of shared pointers of collision rates + vector> m_collisionRates; + + //! The collision-target species indices of #m_collisions + vector m_targetSpeciesIndices; + + //! Interpolated cross sections. This is used for storing + //! interpolated cross sections temporarily. + vector m_interp_cs; + + //! The list of whether the interpolated cross sections is ready + vector m_interp_cs_ready; + + //! Set collisions. This function sets the list of collisions and + //! the list of target species using #addCollision. + void setCollisions(); + + //! Add a collision and record the target species + void addCollision(std::shared_ptr collision); + }; } diff --git a/include/cantera/thermo/ThermoPhase.h b/include/cantera/thermo/ThermoPhase.h index 635015b57a..2baec705db 100644 --- a/include/cantera/thermo/ThermoPhase.h +++ b/include/cantera/thermo/ThermoPhase.h @@ -15,6 +15,7 @@ #include "MultiSpeciesThermo.h" #include "cantera/base/Units.h" #include "cantera/base/AnyMap.h" +#include "cantera/base/Solution.h" namespace Cantera { @@ -1957,6 +1958,12 @@ class ThermoPhase : public Phase //! @} + //! Set the link to the Solution object that owns this ThermoPhase + //! @param soln Weak pointer to the parent Solution object + virtual void setSolution(std::weak_ptr soln) { + m_soln = soln; + } + protected: //! Store the parameters of a ThermoPhase object such that an identical //! one could be reconstructed using the newThermo(AnyMap&) function. This @@ -1993,6 +2000,9 @@ class ThermoPhase : public Phase //! last value of the temperature processed by reference state mutable double m_tlast = 0.0; + + //! reference to Solution + std::weak_ptr m_soln; }; } diff --git a/interfaces/cython/cantera/thermo.pxd b/interfaces/cython/cantera/thermo.pxd index fe1db7e54d..a22beca57f 100644 --- a/interfaces/cython/cantera/thermo.pxd +++ b/interfaces/cython/cantera/thermo.pxd @@ -210,6 +210,7 @@ cdef extern from "cantera/thermo/PlasmaPhase.h": size_t nElectronEnergyLevels() double electronPressure() string electronSpeciesName() + double elasticPowerLoss() except +translate_exception cdef extern from "cantera/cython/thermo_utils.h": diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx index 20c2507b47..bfcc416e15 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/interfaces/cython/cantera/thermo.pyx @@ -1904,6 +1904,15 @@ cdef class ThermoPhase(_SolutionBase): raise ThermoModelMethodError(self.thermo_model) return pystr(self.plasma.electronSpeciesName()) + property elastic_power_loss: + """ + Elastic power loss (J/s/m3) + .. versionadded:: 3.2 + """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.elasticPowerLoss() cdef class InterfacePhase(ThermoPhase): """ A class representing a surface, edge phase """ diff --git a/src/base/Solution.cpp b/src/base/Solution.cpp index 2b929259fd..eecc5a5b38 100644 --- a/src/base/Solution.cpp +++ b/src/base/Solution.cpp @@ -42,6 +42,7 @@ void Solution::setName(const string& name) { void Solution::setThermo(shared_ptr thermo) { m_thermo = thermo; + m_thermo->setSolution(weak_from_this()); for (const auto& [id, callback] : m_changeCallbacks) { callback(); } diff --git a/src/kinetics/ElectronCollisionPlasmaRate.cpp b/src/kinetics/ElectronCollisionPlasmaRate.cpp index afe8a0f9a3..abe9901383 100644 --- a/src/kinetics/ElectronCollisionPlasmaRate.cpp +++ b/src/kinetics/ElectronCollisionPlasmaRate.cpp @@ -70,17 +70,23 @@ void ElectronCollisionPlasmaRate::getParameters(AnyMap& node) const { node["cross-sections"] = m_crossSections; } +void ElectronCollisionPlasmaRate::updateInterpolatedCrossSection( + const vector& sharedLevels) { + m_crossSectionsInterpolated.clear(); + m_crossSectionsInterpolated.reserve(sharedLevels.size()); + for (double level : sharedLevels) { + m_crossSectionsInterpolated.emplace_back(linearInterp(level, + m_energyLevels, m_crossSections)); + } +} + double ElectronCollisionPlasmaRate::evalFromStruct( const ElectronCollisionPlasmaData& shared_data) { // Interpolate cross-sections data to the energy levels of // the electron energy distribution function if (shared_data.levelChanged) { - m_crossSectionsInterpolated.resize(0); - for (double level : shared_data.energyLevels) { - m_crossSectionsInterpolated.push_back(linearInterp(level, - m_energyLevels, m_crossSections)); - } + updateInterpolatedCrossSection(shared_data.energyLevels); } // Map cross sections to Eigen::ArrayXd diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 6298f8ff25..2f2819c134 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -703,6 +703,10 @@ bool Kinetics::addReaction(shared_ptr r, bool resize) m_ready = false; } + for (const auto& [id, callback] : m_reactionAddedCallbacks) { + callback(); + } + return true; } @@ -710,6 +714,13 @@ void Kinetics::modifyReaction(size_t i, shared_ptr rNew) { checkReactionIndex(i); shared_ptr& rOld = m_reactions[i]; + + if (rNew->rate()->type() == "electron-collision-plasma") { + throw CanteraError("Kinetics::modifyReaction", + "Type electron-collision-plasma is not supported. " + "Use the rate object of the reaction to modify the data."); + } + if (rNew->type() != rOld->type()) { throw CanteraError("Kinetics::modifyReaction", "Reaction types are different: {} != {}.", diff --git a/src/kinetics/Reaction.cpp b/src/kinetics/Reaction.cpp index 038f304bde..aa88c393aa 100644 --- a/src/kinetics/Reaction.cpp +++ b/src/kinetics/Reaction.cpp @@ -316,6 +316,19 @@ void Reaction::setRate(shared_ptr rate) "Reaction rate for reaction '{}' must not be empty.", equation()); } m_rate = rate; + for (const auto& [id, callback] : m_setRateCallbacks) { + callback(); + } +} + +void Reaction::registerSetRateCallback(void* id, const function& callback) +{ + m_setRateCallbacks[id] = callback; +} + +void Reaction::removeSetRateCallback(void* id) +{ + m_setRateCallbacks.erase(id); } string Reaction::reactantString() const diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 19a618324b..ce6d57a0c7 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -8,9 +8,16 @@ #include "cantera/thermo/Species.h" #include "cantera/base/global.h" #include "cantera/numerics/funcs.h" +#include "cantera/kinetics/Kinetics.h" +#include "cantera/kinetics/Reaction.h" +#include "cantera/kinetics/ElectronCollisionPlasmaRate.h" namespace Cantera { +namespace { + const double gamma = sqrt(2 * ElectronCharge / ElectronMass); +} + PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_) { initThermoFile(inputFile, id_); @@ -19,7 +26,22 @@ PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_) m_electronEnergyLevels = Eigen::ArrayXd::LinSpaced(m_nPoints, 0.0, 1.0); // initial electron temperature - setElectronTemperature(temperature()); + m_electronTemp = temperature(); + + // resize vectors + m_interp_cs.resize(m_nPoints); +} + +PlasmaPhase::~PlasmaPhase() +{ + if (shared_ptr soln = m_soln.lock()) { + soln->removeChangedCallback(this); + soln->kinetics()->removeReactionAddedCallback(this); + } + for (size_t k = 0; k < nCollisions(); k++) { + // remove callback + m_collisions[k]->removeSetRateCallback(this); + } } void PlasmaPhase::updateElectronEnergyDistribution() @@ -30,6 +52,7 @@ void PlasmaPhase::updateElectronEnergyDistribution() } else if (m_distributionType == "isotropic") { setIsotropicElectronEnergyDistribution(); } + updateElectronEnergyDistDifference(); electronEnergyDistributionChanged(); } @@ -65,8 +88,7 @@ void PlasmaPhase::setIsotropicElectronEnergyDistribution() double c1 = x * std::pow(gamma2, 1.5) / std::pow(gamma1, 2.5); double c2 = x * std::pow(gamma2 / gamma1, x); m_electronEnergyDist = - c1 * m_electronEnergyLevels.sqrt() / - std::pow(meanElectronEnergy(), 1.5) * + c1 / std::pow(meanElectronEnergy(), 1.5) * (-c2 * (m_electronEnergyLevels / meanElectronEnergy()).pow(x)).exp(); checkElectronEnergyDistribution(); @@ -82,13 +104,22 @@ void PlasmaPhase::setMeanElectronEnergy(double energy) { updateElectronEnergyDistribution(); } -void PlasmaPhase::setElectronEnergyLevels(const double* levels, size_t length) +void PlasmaPhase::setElectronEnergyLevels(const double* levels, size_t length, + bool updateEnergyDist) { m_nPoints = length; m_electronEnergyLevels = Eigen::Map(levels, length); checkElectronEnergyLevels(); electronEnergyLevelChanged(); - updateElectronEnergyDistribution(); + if (updateEnergyDist) updateElectronEnergyDistribution(); + m_interp_cs.resize(m_nPoints); + // The cross sections are interpolated on the energy levels + if (nCollisions() > 0) { + for (size_t i = 0; i < m_collisions.size(); i++) { + m_interp_cs_ready[i] = false; + updateInterpolatedCrossSection(i); + } + } } void PlasmaPhase::electronEnergyDistributionChanged() @@ -134,9 +165,7 @@ void PlasmaPhase::setDiscretizedElectronEnergyDist(const double* levels, size_t length) { m_distributionType = "discretized"; - m_nPoints = length; - m_electronEnergyLevels = - Eigen::Map(levels, length); + setElectronEnergyLevels(levels, length, false); m_electronEnergyDist = Eigen::Map(dist, length); checkElectronEnergyLevels(); @@ -144,39 +173,35 @@ void PlasmaPhase::setDiscretizedElectronEnergyDist(const double* levels, normalizeElectronEnergyDistribution(); } checkElectronEnergyDistribution(); + updateElectronEnergyDistDifference(); updateElectronTemperatureFromEnergyDist(); electronEnergyLevelChanged(); electronEnergyDistributionChanged(); } -void PlasmaPhase::setDiscretizedElectronEnergyDist(const double* dist, - size_t length) -{ - m_distributionType = "discretized"; - m_nPoints = length; - m_electronEnergyDist = - Eigen::Map(dist, length); - checkElectronEnergyLevels(); - if (m_do_normalizeElectronEnergyDist) { - normalizeElectronEnergyDistribution(); - } - checkElectronEnergyDistribution(); - updateElectronTemperatureFromEnergyDist(); - electronEnergyDistributionChanged(); -} - void PlasmaPhase::updateElectronTemperatureFromEnergyDist() { // calculate mean electron energy and electron temperature Eigen::ArrayXd eps52 = m_electronEnergyLevels.pow(5./2.); double epsilon_m = 2.0 / 5.0 * numericalQuadrature(m_quadratureMethod, m_electronEnergyDist, eps52); + if (epsilon_m < 0.0 && m_quadratureMethod == "simpson") { + // try trapezoidal method + epsilon_m = 2.0 / 5.0 * numericalQuadrature( + "trapezoidal", m_electronEnergyDist, eps52); + } + + if (epsilon_m < 0.0) { + throw CanteraError("PlasmaPhase::updateElectronTemperatureFromEnergyDist", + "The electron energy distribution produces negative electron temperature."); + } + m_electronTemp = 2.0 / 3.0 * epsilon_m * ElectronCharge / Boltzmann; } void PlasmaPhase::setIsotropicShapeFactor(double x) { m_isotropicShapeFactor = x; - setIsotropicElectronEnergyDistribution(); + updateElectronEnergyDistribution(); } void PlasmaPhase::getParameters(AnyMap& phaseNode) const @@ -202,28 +227,30 @@ void PlasmaPhase::getParameters(AnyMap& phaseNode) const void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) { IdealGasPhase::setParameters(phaseNode, rootNode); + m_root = rootNode; if (phaseNode.hasKey("electron-energy-distribution")) { const AnyMap eedf = phaseNode["electron-energy-distribution"].as(); m_distributionType = eedf["type"].asString(); if (m_distributionType == "isotropic") { if (eedf.hasKey("shape-factor")) { - setIsotropicShapeFactor(eedf["shape-factor"].asDouble()); + m_isotropicShapeFactor = eedf["shape-factor"].asDouble(); } else { throw CanteraError("PlasmaPhase::setParameters", "isotropic type requires shape-factor key."); } + if (eedf.hasKey("energy-levels")) { + setElectronEnergyLevels(eedf["energy-levels"].asVector().data(), + eedf["energy-levels"].asVector().size(), + false); + } if (eedf.hasKey("mean-electron-energy")) { double energy = eedf.convert("mean-electron-energy", "eV"); + // setMeanElectronEnergy() calls updateElectronEnergyDistribution() setMeanElectronEnergy(energy); } else { throw CanteraError("PlasmaPhase::setParameters", "isotropic type requires electron-temperature key."); } - if (eedf.hasKey("energy-levels")) { - setElectronEnergyLevels(eedf["energy-levels"].asVector().data(), - eedf["energy-levels"].asVector().size()); - } - setIsotropicElectronEnergyDistribution(); } else if (m_distributionType == "discretized") { if (!eedf.hasKey("energy-levels")) { throw CanteraError("PlasmaPhase::setParameters", @@ -272,6 +299,186 @@ void PlasmaPhase::initThermo() } } +void PlasmaPhase::setSolution(std::weak_ptr soln) { + ThermoPhase::setSolution(soln); + // register callback function to be executed + // when the thermo or kinetics object changed + if (shared_ptr soln = m_soln.lock()) { + soln->registerChangedCallback(this, [&]() { + setCollisions(); + }); + } +} + +void PlasmaPhase::setCollisions() +{ + m_collisions.clear(); + m_collisionRates.clear(); + m_targetSpeciesIndices.clear(); + + if (shared_ptr soln = m_soln.lock()) { + shared_ptr kin = soln->kinetics(); + if (!kin) { + return; + } + + // add collision from the initial list of reactions + for (size_t i = 0; i < kin->nReactions(); i++) { + std::shared_ptr R = kin->reaction(i); + if (R->rate()->type() != "electron-collision-plasma") { + continue; + } + addCollision(R); + } + + // register callback when reaction is added later + // Modifying collision reactions is not supported + kin->registerReactionAddedCallback(this, [this, kin]() { + size_t i = kin->nReactions() - 1; + if (kin->reaction(i)->type() == "electron-collision-plasma") { + addCollision(kin->reaction(i)); + } + }); + } +} + +void PlasmaPhase::addCollision(std::shared_ptr collision) +{ + size_t i = nCollisions(); + + // setup callback to signal updating the cross-section-related + // parameters + collision->registerSetRateCallback(this, [this, i, collision]() { + m_interp_cs_ready[i] = false; + m_collisionRates[i] = + std::dynamic_pointer_cast(collision->rate()); + }); + + // Identify target species for electron-collision reactions + for (const auto& [name, _] : collision->reactants) { + // Reactants are expected to be electrons and the target species + if (name != electronSpeciesName()) { + m_targetSpeciesIndices.emplace_back(speciesIndex(name)); + break; + } + } + + m_collisions.emplace_back(collision); + m_collisionRates.emplace_back( + std::dynamic_pointer_cast(collision->rate())); + m_interp_cs_ready.emplace_back(false); + + // resize parameters + m_elasticElectronEnergyLossCoefficients.resize(nCollisions()); +} + +bool PlasmaPhase::updateInterpolatedCrossSection(size_t i) +{ + if (m_interp_cs_ready[i]) { + return false; + } + vector levels(m_nPoints); + Eigen::Map(levels.data(), m_nPoints) = m_electronEnergyLevels; + m_collisionRates[i]->updateInterpolatedCrossSection(levels); + m_interp_cs_ready[i] = true; + return true; +} + +void PlasmaPhase::updateElectronEnergyDistDifference() +{ + m_electronEnergyDistDiff.resize(nElectronEnergyLevels()); + // Forward difference for the first point + m_electronEnergyDistDiff[0] = + (m_electronEnergyDist[1] - m_electronEnergyDist[0]) / + (m_electronEnergyLevels[1] - m_electronEnergyLevels[0]); + + // Central difference for the middle points + for (size_t i = 1; i < m_nPoints - 1; i++) { + double h1 = m_electronEnergyLevels[i+1] - m_electronEnergyLevels[i]; + double h0 = m_electronEnergyLevels[i] - m_electronEnergyLevels[i-1]; + m_electronEnergyDistDiff[i] = (h0 * h0 * m_electronEnergyDist[i+1] + + (h1 * h1 - h0 * h0) * m_electronEnergyDist[i] - + h1 * h1 * m_electronEnergyDist[i-1]) / + (h1 * h0) / (h1 + h0); + } + + // Backward difference for the last point + m_electronEnergyDistDiff[m_nPoints-1] = + (m_electronEnergyDist[m_nPoints-1] - + m_electronEnergyDist[m_nPoints-2]) / + (m_electronEnergyLevels[m_nPoints-1] - + m_electronEnergyLevels[m_nPoints-2]); +} + +void PlasmaPhase::updateElasticElectronEnergyLossCoefficients() +{ + // cache of cross section plus distribution plus energy-level number + static const int cacheId = m_cache.getId(); + CachedScalar last_stateNum = m_cache.getScalar(cacheId); + + // combine the distribution and energy level number + int stateNum = m_distNum + m_levelNum; + + vector interpChanged(m_collisions.size()); + for (size_t i = 0; i < m_collisions.size(); i++) { + interpChanged[i] = updateInterpolatedCrossSection(i); + } + + if (last_stateNum.validate(temperature(), stateNum)) { + // check each cross section, and only update coefficients that + // the interpolated cross sections change + for (size_t i = 0; i < m_collisions.size(); i++) { + if (interpChanged[i]) { + updateElasticElectronEnergyLossCoefficient(i); + } + } + } else { + // update every coefficient if distribution, temperature, + // or energy levels change. + for (size_t i = 0; i < m_collisions.size(); i++) { + updateElasticElectronEnergyLossCoefficient(i); + } + } +} + +void PlasmaPhase::updateElasticElectronEnergyLossCoefficient(size_t i) +{ + // @todo exclude attachment collisions + size_t k = m_targetSpeciesIndices[i]; + + // Map cross sections to Eigen::ArrayXd + auto cs_array = Eigen::Map( + m_collisionRates[i]->crossSectionInterpolated().data(), + m_collisionRates[i]->crossSectionInterpolated().size() + ); + + // Mass ratio calculation + double mass_ratio = ElectronMass / molecularWeight(k) * Avogadro; + + // Calculate the rate using Simpson's rule or trapezoidal rule + Eigen::ArrayXd f0_plus = m_electronEnergyDist + Boltzmann * temperature() / + ElectronCharge * m_electronEnergyDistDiff; + m_elasticElectronEnergyLossCoefficients[i] = 2.0 * mass_ratio * gamma * + numericalQuadrature( + m_quadratureMethod, 1.0 / 3.0 * f0_plus.cwiseProduct(cs_array), + m_electronEnergyLevels.pow(3.0)); +} + +double PlasmaPhase::elasticPowerLoss() +{ + updateElasticElectronEnergyLossCoefficients(); + // The elastic power loss includes the contributions from inelastic + // collisions (inelastic recoil effects). + double rate = 0.0; + for (size_t i = 0; i < nCollisions(); i++) { + rate += concentration(m_targetSpeciesIndices[i]) * + m_elasticElectronEnergyLossCoefficients[i]; + } + + return Avogadro * Avogadro * ElectronCharge * + concentration(m_electronSpeciesIndex) * rate; +} + void PlasmaPhase::updateThermo() const { IdealGasPhase::updateThermo(); diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index dfcd4ea788..994d3c1939 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -534,7 +534,7 @@ TEST(Reaction, TwoTempPlasmaFromYaml) TEST(Reaction, ElectronCollisionPlasmaFromYaml) { - auto sol = newSolution("oxygen-plasma.yaml", "", "none"); + auto sol = newSolution("oxygen-plasma.yaml", "discretized-electron-energy-plasma", "none"); AnyMap rxn = AnyMap::fromYamlString( "{equation: O2 + E => E + O2," " type: electron-collision-plasma," diff --git a/test/python/test_thermo.py b/test/python/test_thermo.py index be9bc5b969..60c4f5eac1 100644 --- a/test/python/test_thermo.py +++ b/test/python/test_thermo.py @@ -1183,11 +1183,22 @@ def test_multi_site_unnormalized(self, surf): class TestPlasmaPhase: - @pytest.fixture(scope='function') def phase(self): - return ct.Solution('oxygen-plasma.yaml', 'isotropic-electron-energy-plasma', + phase = ct.Solution('oxygen-plasma.yaml', 'isotropic-electron-energy-plasma', transport_model=None) + phase.isotropic_shape_factor = 1.0 + return phase + + @property + def collision_data(self): + return { + "equation": "O2 + E => E + O2", + "type": "electron-collision-plasma", + "energy-levels": [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0], + "cross-sections": [0.0, 3.83e-20, 4.47e-20, 4.79e-20, 5.07e-20, 5.31e-20, + 5.49e-20, 5.64e-20, 5.77e-20, 5.87e-20, 5.97e-20] + } def test_converting_electron_energy_to_temperature(self, phase): phase.mean_electron_energy = 1.0 @@ -1239,6 +1250,52 @@ def test_add_multiple_electron_species(self, phase): match='Only one electron species is allowed'): phase.add_species(electron) + def test_elastic_power_loss_low_T(self, phase): + phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + assert phase.elastic_power_loss == approx(6846332332) + + def test_elastic_power_loss_high_T(self, phase): + # when T is as high as Te the energy loss rate becomes small + phase.TPX = 4000, ct.one_atm, "O2:1, E:1e-5" + assert phase.elastic_power_loss == approx(2865540) + + def test_elastic_power_loss_replace_rate(self, phase): + phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + rate = ct.ReactionRate.from_dict(self.collision_data) + phase.reaction(1).rate = rate + assert phase.elastic_power_loss == approx(11765800095) + + def test_elastic_power_loss_add_reaction(self, phase): + phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + phase.add_reaction(ct.Reaction.from_dict(self.collision_data, phase)) + assert phase.elastic_power_loss == approx(18612132428) + + def test_elastic_power_loss_change_levels(self, phase): + phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + phase.electron_energy_levels = np.linspace(0,10,101) + assert phase.elastic_power_loss == approx(113058853) + + def test_elastic_power_loss_change_dist(self, phase): + phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + levels = np.array([0.0, 0.1, 1.0, 9.0, 10.0]) + dist = np.array([0.0, 0.2, 0.7, 0.01, 0.01]) + # set the electron energy levels first to test if + # set_discretized_electron_energy_distribution triggers + # updating the interpolated cross sections + phase.electron_energy_levels = np.array([0.0, 0.1, 1.0, 7.0, 10.0]) + phase.normalize_electron_energy_distribution_enabled = False + phase.set_discretized_electron_energy_distribution(levels, dist) + assert phase.elastic_power_loss == approx(7568518396) + + def test_elastic_power_loss_change_mean_electron_energy(self, phase): + phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + phase.mean_electron_energy = 2.0 + assert phase.elastic_power_loss == approx(5826212349) + + def test_elastic_power_loss_change_shape_factor(self, phase): + phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + phase.isotropic_shape_factor = 1.1 + assert phase.elastic_power_loss == approx(4273351243) class TestImport: """