Plasma thermo

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;

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":

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):

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;

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]

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