From 65526d5512f80de76b5ecff76913974b92247aff Mon Sep 17 00:00:00 2001 From: Bang-Shiuh Chen Date: Sat, 10 May 2025 16:59:28 -0400 Subject: [PATCH 1/2] [plasma] update plasma thermo definition --- include/cantera/thermo/PlasmaPhase.h | 46 ++++++++++++++++++++++++++-- interfaces/cython/cantera/thermo.pxd | 3 ++ interfaces/cython/cantera/thermo.pyx | 27 ++++++++++++++++ src/thermo/PlasmaPhase.cpp | 2 +- 4 files changed, 74 insertions(+), 4 deletions(-) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 266fbfe6466..7138dea87c5 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -206,6 +206,15 @@ class PlasmaPhase: public IdealGasPhase return electronTemperature() * GasConstant; } + /** + * Pressure. Units: Pa. + * For an ideal gas mixture, + * @f[ P = n \hat R T. @f] + */ + double pressure() const override { + return GasConstant * gasMolarDensity() * temperature(); + } + /** * Electron pressure. Units: Pa. * @f[P = n_{k_e} R T_e @f] @@ -230,11 +239,24 @@ class PlasmaPhase: public IdealGasPhase return m_electronSpeciesIndex; } + //! Return the molar density of heavy species (non-electron species) + /*! + * Calculates the molar density of the gas excluding electrons by multiplying + * the total molar density by the fraction of non-electron species. + * + * Units: kmol/m³ + * + * @returns The molar density of heavy species (all species except electrons) + */ + double gasMolarDensity() const { + return molarDensity() * (1.0 - moleFraction(m_electronSpeciesIndex)); + } + //! Return the Molar enthalpy. Units: J/kmol. /*! - * For an ideal gas mixture with additional electron, + * For an ideal gas mixture excluding electron, * @f[ - * \hat h(T) = \sum_{k \neq k_e} X_k \hat h^0_k(T) + X_{k_e} \hat h^0_{k_e}(T_e), + * \hat h(T) = \sum_{k \neq k_e} X_k \hat h^0_k(T), * @f] * and is a function only of temperature. The standard-state pure-species * enthalpies @f$ \hat h^0_k(T) @f$ are computed by the species @@ -244,6 +266,19 @@ class PlasmaPhase: public IdealGasPhase */ double enthalpy_mole() const override; + //! Return the Molar enthalpy of electron. Units: J/kmol. + /*! + * @f[ + * \hat h(T) = X_{k_e} \hat h^0_{k_e}(T_e), + * @f] + * and is a function only of electron temperature. + */ + double electronEnthalpy_mole() const { + return GasConstant * electronTemperature() * + moleFraction(m_electronSpeciesIndex) * + enthalpy_RT_ref()[m_electronSpeciesIndex]; + } + double cp_mole() const override { throw NotImplementedError("PlasmaPhase::cp_mole"); } @@ -257,7 +292,12 @@ class PlasmaPhase: public IdealGasPhase } double intEnergy_mole() const override { - throw NotImplementedError("PlasmaPhase::intEnergy_mole"); + return enthalpy_mole() - pressure() / gasMolarDensity(); + } + + double electronIntEnergy_mole() const { + return electronEnthalpy_mole() - electronPressure() / + concentration(m_electronSpeciesIndex); } void getEntropy_R(double* sr) const override; diff --git a/interfaces/cython/cantera/thermo.pxd b/interfaces/cython/cantera/thermo.pxd index a22beca57f0..72438f1916e 100644 --- a/interfaces/cython/cantera/thermo.pxd +++ b/interfaces/cython/cantera/thermo.pxd @@ -211,6 +211,9 @@ cdef extern from "cantera/thermo/PlasmaPhase.h": double electronPressure() string electronSpeciesName() double elasticPowerLoss() except +translate_exception + double electronEnthalpy_mole() except +translate_exception + double electronIntEnergy_mole() except +translate_exception + double gasMolarDensity() 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..a1796e5f22e 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/interfaces/cython/cantera/thermo.pyx @@ -1914,6 +1914,33 @@ cdef class ThermoPhase(_SolutionBase): raise ThermoModelMethodError(self.thermo_model) return self.plasma.elasticPowerLoss() + property electron_enthalpy_mole: + """ + Molar enthalpy of electron (J/kmol) + """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.electronEnthalpy_mole() + + property electron_int_energy_mole: + """ + Molar internal energy of electron (J/kmol) + """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.electronIntEnergy_mole() + + property gas_molar_density: + """ + Molar density of heavy species (non-electron species) in kmol/m³. + """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.gasMolarDensity() + cdef class InterfacePhase(ThermoPhase): """ A class representing a surface, edge phase """ def __cinit__(self, *args, **kwargs): diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 7bc60e04f9d..a2cec63eb51 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -501,7 +501,7 @@ void PlasmaPhase::updateThermo() const double PlasmaPhase::enthalpy_mole() const { double value = IdealGasPhase::enthalpy_mole(); - value += GasConstant * (electronTemperature() - temperature()) * + value -= GasConstant * temperature() * moleFraction(m_electronSpeciesIndex) * m_h0_RT[m_electronSpeciesIndex]; return value; From 8c607014f6b22397afb0dfbb60b8a34bd3d22257 Mon Sep 17 00:00:00 2001 From: Bang-Shiuh Chen Date: Sun, 11 May 2025 22:10:05 -0400 Subject: [PATCH 2/2] [test] add test for plasma thermo properties --- test/data/oxygen-electron.yaml | 38 ++++++++++++++++++++++++++++++++ test/python/test_thermo.py | 40 ++++++++++++++++++++++++++++++++++ 2 files changed, 78 insertions(+) create mode 100644 test/data/oxygen-electron.yaml diff --git a/test/data/oxygen-electron.yaml b/test/data/oxygen-electron.yaml new file mode 100644 index 00000000000..20fee544cbf --- /dev/null +++ b/test/data/oxygen-electron.yaml @@ -0,0 +1,38 @@ +description: |- + This file is for testing purpose. We use constant Cp for electron with no Tmax. + +units: {length: cm, quantity: molec, activation-energy: K} + +phases: +- name: oxygen-electron + thermo: ideal-gas + elements: [O, E] + species: + - species: [E] + - nasa_gas.yaml/species: [O2] + + kinetics: gas + reactions: none + transport: none + +species: +- name: E + composition: {E: 1} + thermo: + model: constant-cp + T0: 200 K + h0: -2.04 kJ/mol + s0: 12.679 J/mol/K + cp0: 20.786 J/mol/K + +# reactions: +# - equation: E + O2 + O2 => O2- + O2 +# type: two-temperature-plasma +# rate-constant: {A: 4.2e-27, b: -1.0, Ea-gas: 600, Ea-electron: 700} + +# - equation: O2 + E => E + O2 +# type: electron-collision-plasma +# note: This is a electron collision process of plasma +# energy-levels: [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0] +# cross-sections: [0.0, 5.97e-20, 6.45e-20, 6.74e-20, 6.93e-20, 7.2e-20, +# 7.52e-20, 7.86e-20, 8.21e-20, 8.49e-20, 8.8e-20] diff --git a/test/python/test_thermo.py b/test/python/test_thermo.py index 1d880285612..9276cfa47f4 100644 --- a/test/python/test_thermo.py +++ b/test/python/test_thermo.py @@ -1190,6 +1190,11 @@ def phase(self): phase.isotropic_shape_factor = 1.0 return phase + @pytest.fixture(scope='function') + def gas(self): + gas = ct.Solution('oxygen-electron.yaml') + return gas + @property def collision_data(self): return { @@ -1297,6 +1302,41 @@ 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_enthalpy_mole(self, phase, gas): + phase.TPX = 1000, ct.one_atm, "O2:1, E:1" + gas.TPX = 1000, ct.one_atm, "O2:1" + assert (2.0 * phase.enthalpy_mole) == approx(gas.enthalpy_mole) + + def test_electron_enthalpy_mole(self, phase, gas): + phase.TPX = 1000, ct.one_atm, "E:1" + phase.Te = 1000 + gas.TPX = 1000, ct.one_atm, "E:1" + assert (phase.electron_enthalpy_mole) == approx(gas.enthalpy_mole) + + def test_pressure(self, phase, gas): + phase.TPX = 1000, ct.one_atm, "O2:1, E:1" + gas.TPX = 1000, ct.one_atm, "O2:1" + assert (2.0 * phase.P) == approx(gas.P) + + def test_electron_pressure(self, phase, gas): + phase.TPX = 1000, ct.one_atm, "E:1" + phase.Te = 2000 + gas.TPX = 1000, 2.0 * ct.one_atm, "E:1" + assert phase.Pe == approx(gas.P) + + def test_internal_energy_mole(self, phase): + phase.TPX = 1000, ct.one_atm, "O2:1, E:1" + phase.Te = 1000 + u_mole = phase.enthalpy_mole - phase.P / phase.gas_molar_density + assert phase.int_energy_mole == approx(u_mole) + + def test_electron_internal_energy_mole(self, phase): + phase.TPX = 1000, ct.one_atm, "O2:1, E:1" + phase.Te = 1000 + u_mole = (phase.electron_enthalpy_mole - phase.P / + phase.concentrations[phase.species_index('E')]) + assert phase.electron_int_energy_mole == approx(u_mole) + class TestImport: """ Tests the various ways of creating a Solution object