[thermo] fix plasma thermo properties

include/cantera/thermo/PlasmaPhase.h

@@ -187,11 +187,78 @@ class PlasmaPhase: public IdealGasPhase
         return m_electronTemp;
     }
 
+    //! Return the Gas Constant multiplied by the current electron temperature
+    /*!
+     *  The units are Joules kmol-1
+     */
+    double RTe() const {
+        return electronTemperature() * GasConstant;
+    }
+
+    //! Pressure
+    //! Units: Pa. For an ideal gas mixture with additional electrons,
+    //! \f[
+    //!        P = \sum_{k \neq k_e} n_k R T.
+    //! \f]
+    virtual double pressure() const {
+        double sum = 0.0;
+        for (size_t k = 0; k < m_kk; k++) {
+            if (k != m_electronSpeciesIndex) {
+                sum += GasConstant * concentration(k) * temperature();
+            }
+        }
+        return sum;
+    }
+
+    /**
+     * Electron pressure. Units: Pa.
+     * \f[P = n_{k_e} R T_e\f]
+     */
+    virtual double electronPressure() const {
+        return GasConstant * concentration(m_electronSpeciesIndex) *
+               electronTemperature();
+    }
+
     //! Number of electron levels
     size_t nElectronEnergyLevels() const {
         return m_nPoints;
     }
 
+    //! Electron Species Index
+    size_t electronSpeciesIndex() const {
+        return m_electronSpeciesIndex;
+    }
+
+    //! Return the Molar enthalpy. Units: J/kmol.
+    /*!
+     * For an ideal gas mixture with additional 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),
+     * \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
+     * thermodynamic property manager.
+     *
+     * \see MultiSpeciesThermo
+     */
+    virtual double enthalpy_mole() const;
+
+    virtual double entropy_mole() const;
+
+    virtual double gibbs_mole() const;
+
+    virtual void getGibbs_ref(double* g) const;
+
+    virtual void getStandardVolumes_ref(double* vol) const;
+
+    virtual void getStandardChemPotentials(double* mu) const;
+
+    virtual void getChemPotentials(double* mu) const;
+
+    virtual void getPartialMolarEnthalpies(double* hbar) const;
+
+    virtual void getPartialMolarIntEnergies(double* ubar) const;
+
     virtual void getParameters(AnyMap& phaseNode) const;
 
     virtual void setParameters(const AnyMap& phaseNode,

src/thermo/IdealGasPhase.cpp

@@ -54,8 +54,7 @@ void IdealGasPhase::getActivityCoefficients(doublereal* ac) const
 
 void IdealGasPhase::getStandardChemPotentials(doublereal* muStar) const
 {
-    const vector_fp& gibbsrt = gibbs_RT_ref();
-    scale(gibbsrt.begin(), gibbsrt.end(), muStar, RT());
+    getGibbs_ref(muStar);
     double tmp = log(pressure() / refPressure()) * RT();
     for (size_t k = 0; k < m_kk; k++) {
         muStar[k] += tmp; // add RT*ln(P/P_0)
@@ -152,10 +151,7 @@ void IdealGasPhase::getPureGibbs(doublereal* gpure) const
 
 void IdealGasPhase::getIntEnergy_RT(doublereal* urt) const
 {
-    const vector_fp& _h = enthalpy_RT_ref();
-    for (size_t k = 0; k < m_kk; k++) {
-        urt[k] = _h[k] - 1.0;
-    }
+    getIntEnergy_RT_ref(urt);
 }
 
 void IdealGasPhase::getCp_R(doublereal* cpr) const

src/thermo/PlasmaPhase.cpp

@@ -268,4 +268,70 @@ void PlasmaPhase::updateThermo() const
     m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
 }
 
+double PlasmaPhase::enthalpy_mole() const {
+    double value = IdealGasPhase::enthalpy_mole();
+    value += GasConstant * (electronTemperature() - temperature()) *
+             moleFraction(m_electronSpeciesIndex) *
+             m_h0_RT[m_electronSpeciesIndex];
+    return value;
+}
+
+double PlasmaPhase::entropy_mole() const {
+    warn_user("PlasmaPhase::entropy_mole",
+              "Use the same equation of IdealGasPhase::entropy_mole "
+              "which is not correct for plasma.");
+    return IdealGasPhase::entropy_mole();
+}
+
+double PlasmaPhase::gibbs_mole() const {
+    warn_user("PlasmaPhase::gibbs_mole",
+              "Use the same equation of IdealGasPhase::gibbs_mole "
+              "which is not correct for plasma.");
+    return IdealGasPhase::gibbs_mole();
+}
+
+void PlasmaPhase::getGibbs_ref(double* g) const
+{
+    IdealGasPhase::getGibbs_ref(g);
+    g[m_electronSpeciesIndex] *= electronTemperature() / temperature();
+}
+
+void PlasmaPhase::getStandardVolumes_ref(double* vol) const
+{
+    IdealGasPhase::getStandardVolumes_ref(vol);
+    vol[m_electronSpeciesIndex] *= electronTemperature() / temperature();
+}
+
+void PlasmaPhase::getStandardChemPotentials(double* muStar) const
+{
+    warn_user("PlasmaPhase::getStandardChemPotentials",
+              "Use the same equation of IdealGasPhase::getStandardChemPotentials "
+              "which is not correct for plasma.");
+    IdealGasPhase::getStandardChemPotentials(muStar);
+}
+
+void PlasmaPhase::getChemPotentials(double* mu) const
+{
+    warn_user("PlasmaPhase::getChemPotentials",
+              "Use the same equation of IdealGasPhase::getChemPotentials "
+              "which is not correct for plasma.");
+    IdealGasPhase::getChemPotentials(mu);
+}
+
+void PlasmaPhase::getPartialMolarEnthalpies(double* hbar) const
+{
+    IdealGasPhase::getPartialMolarEnthalpies(hbar);
+    hbar[m_electronSpeciesIndex] *= electronTemperature() / temperature();
+}
+
+void PlasmaPhase::getPartialMolarIntEnergies(double* ubar) const
+{
+    const vector_fp& _h = enthalpy_RT_ref();
+    for (size_t k = 0; k < m_kk; k++) {
+        ubar[k] = RT() * (_h[k] - 1.0);
+    }
+    size_t k = m_electronSpeciesIndex;
+    ubar[k] = RTe() * (_h[k] - 1.0);
+}
+
 }