From 8fc6bff76e75f02b2ad2a8fdf29ba62fb11c8409 Mon Sep 17 00:00:00 2001 From: Bang-Shiuh Chen Date: Tue, 20 May 2025 11:38:06 -0400 Subject: [PATCH] [plasma] add inelastic power loss --- include/cantera/kinetics/Kinetics.h | 6 ++++ include/cantera/thermo/PlasmaPhase.h | 27 ++++++++++++++- interfaces/cython/cantera/thermo.pxd | 1 + interfaces/cython/cantera/thermo.pyx | 10 ++++++ src/kinetics/Kinetics.cpp | 12 +++++++ src/thermo/PlasmaPhase.cpp | 52 ++++++++++++++++++++++------ test/data/oxygen-plasma.yaml | 13 ++++++- test/python/test_thermo.py | 21 +++++++++++ 8 files changed, 130 insertions(+), 12 deletions(-) diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index aa96ec43ffe..7d3748b8d4d 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -1245,6 +1245,12 @@ class Kinetics throw NotImplementedError("Kinetics::getRevRateConstants"); } + //! Net rates of progress by reaction indices + //! @param indices The indices of the reactions for which to return the net rates + //! of progress. + //! @return A vector of the net rates of progress for the specified reactions. + virtual vector netRatesOfProgressByIndices(const vector& indices); + //! @} //! @name Reaction Mechanism Construction //! @{ diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 266fbfe6466..21d9678c2ce 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/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 rates of progress from #m_kin + void updateCollisionRatesOfProgress(); + 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; + //! Net rates of progress of collisions + vector m_netROPCollisions; + + //! Kinetic object + shared_ptr m_kin; + //! The list of shared pointers of plasma collision reactions vector> m_collisions; @@ -449,12 +469,17 @@ class PlasmaPhase: public IdealGasPhase //! The list of whether the interpolated cross sections is ready vector m_interp_cs_ready; + //! Indices of electron collision plasma reactions + //! This is used for getting the rates of progress. + vector 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 collision); + //! @param j index of the corresponding reaction for the collision + void addCollision(size_t j); }; diff --git a/interfaces/cython/cantera/thermo.pxd b/interfaces/cython/cantera/thermo.pxd index a22beca57f0..eec52444eee 100644 --- a/interfaces/cython/cantera/thermo.pxd +++ b/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": diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx index b091f708829..32022e7ab3f 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/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): diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 1f96cf4ee2b..f54bbcea30b 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -413,6 +413,18 @@ void Kinetics::getNetProductionRates(double* net) m_reactantStoich.decrementSpecies(m_ropnet.data(), net); } +vector Kinetics::netRatesOfProgressByIndices(const vector& indices) +{ + updateROP(); + + vector rates(indices.size()); + for (size_t i = 0; i < indices.size(); i++) { + rates[i] = m_ropnet[indices[i]]; + } + + return rates; +} + void Kinetics::getCreationRates_ddT(double* dwdot) { Eigen::Map out(dwdot, m_kk); diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 7bc60e04f9d..7b4a0316fc9 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/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 soln = m_soln.lock()) { - shared_ptr 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 R = kin->reaction(i); + for (size_t j = 0; j < m_kin->nReactions(); j++) { + std::shared_ptr 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 collision) +void PlasmaPhase::addCollision(size_t j) { + std::shared_ptr 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 collision) std::dynamic_pointer_cast(collision->rate())); m_interp_cs_ready.emplace_back(false); + m_electronCollisionReactionIndices.push_back(j); + // resize parameters m_elasticElectronEnergyLossCoefficients.resize(nCollisions()); + m_netROPCollisions.resize(nCollisions()); } bool PlasmaPhase::updateInterpolatedCrossSection(size_t i) @@ -464,6 +469,33 @@ void PlasmaPhase::updateElasticElectronEnergyLossCoefficient(size_t i) m_electronEnergyLevels.pow(3.0)); } +void PlasmaPhase::updateCollisionRatesOfProgress() +{ + 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 + if (!last.validate(stateNum)) { + m_netROPCollisions = + m_kin->netRatesOfProgressByIndices(m_electronCollisionReactionIndices); + } +} + +double PlasmaPhase::inelasticPowerLoss() +{ + updateCollisionRatesOfProgress(); + double powerLoss = 0.0; + for (size_t i = 0; i < nCollisions(); i++) { + powerLoss += m_netROPCollisions[i] * m_collisionRates[i]->energyLevels()[0]; + } + powerLoss *= Avogadro * ElectronCharge; + return powerLoss; +} + double PlasmaPhase::elasticPowerLoss() { updateElasticElectronEnergyLossCoefficients(); diff --git a/test/data/oxygen-plasma.yaml b/test/data/oxygen-plasma.yaml index 9c172c208bb..ff2530433f9 100644 --- a/test/data/oxygen-plasma.yaml +++ b/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 diff --git a/test/python/test_thermo.py b/test/python/test_thermo.py index 1d880285612..e2747859bb5 100644 --- a/test/python/test_thermo.py +++ b/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(285091336) + # 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