Add inelastic power loss to PlasmaPhase

include/cantera/kinetics/BulkKinetics.h

@@ -81,6 +81,9 @@ class BulkKinetics : public Kinetics
     Eigen::SparseMatrix<double> netRatesOfProgress_ddCi() override;
     //! @}
 
+    virtual vector<double> fwdRateConstantsByIndices(
+        const vector<size_t>& indices) override;
+
     //! @name Rate calculation intermediate methods
     //! @{
 

include/cantera/kinetics/Kinetics.h

@@ -1244,6 +1244,11 @@
         throw NotImplementedError("Kinetics::getRevRateConstants");
     }
 
+    //! Forward rate constant by reaction index
+    virtual vector<double> fwdRateConstantsByIndices(const vector<size_t>& indices) {
+        throw NotImplementedError("Kinetics::fwdRateConstant");
+    }
+
     //! @}
     //! @name Reaction Mechanism Construction
     //! @{

include/cantera/thermo/PlasmaPhase.h

@@ -16,6 +16,7 @@ namespace Cantera
 {
 
 class Reaction;
+class Kinetics;
 class ElectronCollisionPlasmaRate;
 
 //! Base class for handling plasma properties, specifically focusing on the
@@ -311,6 +312,16 @@ class PlasmaPhase: public IdealGasPhase
      */
     double elasticPowerLoss();
 
+    /**
+     * The inelastic power loss (J/s/m³)
+     *   @f[
+     *     P_i = N_A e \sum_i U_i ROP_i,
+     *   @f]
+     * where @f$ U_i @f$ is the threshold energy (eV) and @f$ ROP_i @f$ is the
+     * net rate of progress of the collision (kmol/s/m³).
+     */
+    double inelasticPowerLoss();
+
 protected:
     void updateThermo() const override;
 
@@ -424,6 +435,9 @@ class PlasmaPhase: public IdealGasPhase
     */
     void updateElasticElectronEnergyLossCoefficients();
 
+    //! Update collision rate constants from #m_kin
+    void updateCollisionRateConstants();
+
 private:
     //! Electron energy distribution change variable. Whenever
     //! #m_electronEnergyDist changes, this int is incremented.
@@ -433,6 +447,12 @@ class PlasmaPhase: public IdealGasPhase
     //! #m_electronEnergyLevels changes, this int is incremented.
     int m_levelNum = -1;
 
+    //! Forward rate constants of collisions
+    vector<double> m_fwd_collision_rate_constants;
+
+    //! Kinetic object
+    shared_ptr<Kinetics> m_kin;
+
     //! The list of shared pointers of plasma collision reactions
     vector<shared_ptr<Reaction>> m_collisions;
 
@@ -449,12 +469,17 @@ class PlasmaPhase: public IdealGasPhase
     //! The list of whether the interpolated cross sections is ready
     vector<bool> m_interp_cs_ready;
 
+    //! Indices of electron collision plasma reactions
+    //! This is used for getting the rates of progress.
+    vector<size_t> m_electronCollisionReactionIndices;
+
     //! 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<Reaction> collision);
+    //! @param j index of the corresponding reaction for the collision
+    void addCollision(size_t j);
 
 };
 

interfaces/cython/cantera/thermo.pxd

@@ -211,6 +211,7 @@ cdef extern from "cantera/thermo/PlasmaPhase.h":
         double electronPressure()
         string electronSpeciesName()
         double elasticPowerLoss() except +translate_exception
+        double inelasticPowerLoss() except +translate_exception
 
 
 cdef extern from "cantera/cython/thermo_utils.h":

interfaces/cython/cantera/thermo.pyx

@@ -1914,6 +1914,16 @@ cdef class ThermoPhase(_SolutionBase):
                 raise ThermoModelMethodError(self.thermo_model)
             return self.plasma.elasticPowerLoss()
 
+    property inelastic_power_loss:
+        """
+        Inelastic power loss (J/s/m3)
+        .. versionadded:: 3.2
+        """
+        def __get__(self):
+            if not self._enable_plasma:
+                raise ThermoModelMethodError(self.thermo_model)
+            return self.plasma.inelasticPowerLoss()
+
 cdef class InterfacePhase(ThermoPhase):
     """ A class representing a surface, edge phase """
     def __cinit__(self, *args, **kwargs):

src/kinetics/BulkKinetics.cpp

@@ -454,6 +454,16 @@ Eigen::SparseMatrix<double> BulkKinetics::netRatesOfProgress_ddCi()
     return jac - calculateCompositionDerivatives(m_revProductStoich, rop_rates, false);
 }
 
+vector<double> BulkKinetics::fwdRateConstantsByIndices(const vector<size_t>& indices)
+{
+    updateROP();
+    vector<double> rates(indices.size());
+    for (size_t i = 0; i < indices.size(); i++) {
+        rates[i] = m_kf0[indices[i]];
+    }
+    return rates;
+}
+
 void BulkKinetics::updateROP()
 {
     static const int cacheId = m_cache.getId();

src/thermo/PlasmaPhase.cpp

@@ -315,35 +315,37 @@ void PlasmaPhase::setCollisions()
     m_collisions.clear();
     m_collisionRates.clear();
     m_targetSpeciesIndices.clear();
+    m_electronCollisionReactionIndices.clear();
 
     if (shared_ptr<Solution> soln = m_soln.lock()) {
-        shared_ptr<Kinetics> kin = soln->kinetics();
-        if (!kin) {
+        m_kin = soln->kinetics();
+        if (!m_kin) {
             return;
         }
 
         // add collision from the initial list of reactions
-        for (size_t i = 0; i < kin->nReactions(); i++) {
-            std::shared_ptr<Reaction> R = kin->reaction(i);
+        for (size_t j = 0; j < m_kin->nReactions(); j++) {
+            std::shared_ptr<Reaction> R = m_kin->reaction(j);
             if (R->rate()->type() != "electron-collision-plasma") {
                 continue;
             }
-            addCollision(R);
+            addCollision(j);
         }
 
         // 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));
+        m_kin->registerReactionAddedCallback(this, [this]() {
+            size_t j = m_kin->nReactions() - 1;
+            if (m_kin->reaction(j)->type() == "electron-collision-plasma") {
+                addCollision(j);
             }
         });
     }
 }
 
-void PlasmaPhase::addCollision(std::shared_ptr<Reaction> collision)
+void PlasmaPhase::addCollision(size_t j)
 {
+    std::shared_ptr<Reaction> collision = m_kin->reaction(j);
     size_t i = nCollisions();
 
     // setup callback to signal updating the cross-section-related
@@ -368,8 +370,11 @@ void PlasmaPhase::addCollision(std::shared_ptr<Reaction> collision)
         std::dynamic_pointer_cast<ElectronCollisionPlasmaRate>(collision->rate()));
     m_interp_cs_ready.emplace_back(false);
 
+    m_electronCollisionReactionIndices.push_back(j);
+
     // resize parameters
     m_elasticElectronEnergyLossCoefficients.resize(nCollisions());
+    m_fwd_collision_rate_constants.resize(nCollisions());
 }
 
 bool PlasmaPhase::updateInterpolatedCrossSection(size_t i)
@@ -464,6 +469,35 @@ void PlasmaPhase::updateElasticElectronEnergyLossCoefficient(size_t i)
             m_electronEnergyLevels.pow(3.0));
 }
 
+void PlasmaPhase::updateCollisionRateConstants()
+{
+    static const int cacheId = m_cache.getId();
+    CachedScalar last = m_cache.getScalar(cacheId);
+
+    // combine the distribution and energy level number
+    int stateNum = m_distNum + m_levelNum;
+
+    // update only if the reaction rates coefficients have changed
+    // which depends on the energy distribution, and energy levels
+    // which depends on the energy distribution, and energy levels
+    if (!last.validate(stateNum)) {
+        m_fwd_collision_rate_constants =
+            m_kin->fwdRateConstantsByIndices(m_electronCollisionReactionIndices);
+    }
+}
+
+double PlasmaPhase::inelasticPowerLoss()
+{
+    updateCollisionRateConstants();
+    double rate = 0.0;
+    for (size_t i = 0; i < nCollisions(); i++) {
+        rate += concentration(m_targetSpeciesIndices[i]) *
+                m_collisionRates[i]->energyLevels()[0] *
+                m_fwd_collision_rate_constants[i];
+    }
+    return Avogadro * ElectronCharge * concentration(m_electronSpeciesIndex) * rate;
+}
+
 double PlasmaPhase::elasticPowerLoss()
 {
     updateElasticElectronEnergyLossCoefficients();

test/data/oxygen-plasma.yaml

@@ -8,7 +8,7 @@ phases:
   thermo: plasma
   elements: [O, E]
   species:
-  - species: [E]
+  - species: [E, O2(a1dg)]
   - nasa_gas.yaml/species: [O2, O2-]
 
   kinetics: gas
@@ -46,6 +46,17 @@ species:
     s0: 12.679 J/mol/K
     cp0: 20.786 J/mol/K
 
+- name: O2(a1dg)
+  composition: {O: 2}
+  thermo:
+    model: NASA7
+    temperature-ranges: [200.0, 1000.0, 6000.0]
+    data:
+    - [3.78535371, -3.2192854e-03, 1.12323443e-05, -1.17254068e-08, 4.17659585e-12,
+      1.02922572e+04, 3.27320239]
+    - [3.45852381, 1.04045351e-03, -2.79664041e-07, 3.11439672e-11, -8.55656058e-16,
+      1.02229063e+04, 4.15264119]
+
 reactions:
 - equation: E + O2 + O2 => O2- + O2
   type: two-temperature-plasma

test/python/test_thermo.py

@@ -1200,6 +1200,17 @@ def collision_data(self):
                                5.49e-20, 5.64e-20, 5.77e-20, 5.87e-20, 5.97e-20]
         }
 
+    @property
+    def inelastic_collision_data(self):
+        return {
+            "equation": "O2 + E => E + O2(a1dg)",
+            "type": "electron-collision-plasma",
+            "note": "This is a electron collision process of plasma",
+            "energy-levels": [0.977, 1.5, 2.0, 3.0, 3.5, 4.0, 5.0],
+            "cross-sections": [0.0, 5.8000e-23, 1.5300e-22, 3.8000e-22, 4.9000e-22,
+                               5.7000e-22, 7.4000e-22]
+        }
+
     def test_converting_electron_energy_to_temperature(self, phase):
         phase.mean_electron_energy = 1.0
         Te = 2.0 / 3.0 * ct.electron_charge / ct.boltzmann
@@ -1297,6 +1308,16 @@ def test_elastic_power_loss_change_shape_factor(self, phase):
         phase.isotropic_shape_factor = 1.1
         assert phase.elastic_power_loss == approx(7408711810)
 
+    def test_inelastic_power_loss(self, phase):
+        phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5"
+        phase.add_reaction(ct.Reaction.from_dict(self.inelastic_collision_data, phase))
+        value = phase.inelastic_power_loss
+        assert value == approx(1470414640)
+        # change of gas temperature does not change inelastic power loss
+        # when the concentrations hold constant
+        phase.TDX = 2000, phase.density, "O2:1, E:1e-5"
+        assert phase.inelastic_power_loss == approx(value)
+
 class TestImport:
     """
     Tests the various ways of creating a Solution object