From 996107a4e301d6375248646cf63b2968837cc363 Mon Sep 17 00:00:00 2001 From: pag1pag <66677348+pag1pag@users.noreply.github.com> Date: Tue, 26 May 2026 12:00:35 +0200 Subject: [PATCH 01/11] [Plasma] Reorder code to match method grouping of IdealGasPhase --- include/cantera/thermo/PlasmaPhase.h | 297 ++++++++++-------- src/thermo/PlasmaPhase.cpp | 436 +++++++++++++++------------ 2 files changed, 412 insertions(+), 321 deletions(-) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 61bb85c9338..67fd11adc8a 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -105,6 +105,37 @@ class PlasmaPhase: public IdealGasPhase void initThermo() override; + // ================================================================= // + // ================================================================= // + //! @name Overridden from IdealGasPhase or ThermoPhase + //! @{ + bool addSpecies(shared_ptr spec) override; + virtual void setSolution(std::weak_ptr soln) override; + void getParameters(AnyMap& phaseNode) const override; + void setParameters(const AnyMap& phaseNode, + const AnyMap& rootNode=AnyMap()) override; + //! @} + // ================================================================= // + // ================================================================= // + //! @name Electron Species Information + //! @{ + + //! Electron Species Index + size_t electronSpeciesIndex() const { + return m_electronSpeciesIndex; + } + + //! Electron species name + string electronSpeciesName() const { + return speciesName(m_electronSpeciesIndex); + } + + //! @} + // ================================================================= // + // ================================================================= // + //! @name Electron Energy Distribution Functions + //! @{ + //! Set electron energy levels. //! @param levels The vector of electron energy levels (eV). //! Length: #m_nPoints. @@ -186,31 +217,12 @@ class PlasmaPhase: public IdealGasPhase return m_do_normalizeElectronEnergyDist; } - bool addSpecies(shared_ptr spec) override; - //! Electron Temperature (K) //! @return The electron temperature of the phase double electronTemperature() const override { 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; - } - - /** - * 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; @@ -234,139 +246,50 @@ class PlasmaPhase: public IdealGasPhase return m_collisionRates[i]; } - //! 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 - */ - double enthalpy_mole() const override; - - //! Return the molar entropy. Units: J/kmol/K. - /*! - * For an ideal gas mixture with an additional electron species, - * @f[ - * \hat s(T,T_e,p,X) - * = \sum_{k\neq k_e} X_k\left[\hat s_k^0(T) - R\ln\left(\frac{X_k p}{p^0}\right)\right] - * + X_{k_e}\left[\hat s_{k_e}^0(T_e) - R\ln\left(\frac{X_{k_e} p_e}{p^0}\right)\right], - * @f] - * where heavy-species properties are evaluated at @f$T@f$, electron properties at - * @f$T_e@f$, and the electron mixing term uses the electron pressure - * @f$p_e = n_{k_e} R T_e@f$. - * - * @see MultiSpeciesThermo - */ - double entropy_mole() const override; - - //! Return the molar Gibbs free energy. Units: J/kmol. - /*! - * For an ideal gas mixture with an additional electron species, - * @f[ - * \hat g(T, T_e, p, X) = \sum_k X_k \mu_k(T_k, p_k, X), - * @f] - * where heavy species use the gas temperature @f$T@f$ and bulk pressure, - * while the electron chemical potential uses @f$T_e@f$ and - * @f$p_e = n_{k_e} R T_e@f$. - * - * @see MultiSpeciesThermo - */ - double gibbs_mole() const override; - - //! Return the molar internal energy. Units: J/kmol. - /*! - * For an ideal gas mixture with an additional electron species, - * @f[ - * \hat u(T,T_e) = \sum_{k \neq k_e} X_k \hat u^0_k(T) + X_{k_e} \hat u^0_{k_e}(T_e), - * @f] - * where @f$\hat u^0_k(T) = \hat h^0_k(T) - R T@f$ for heavy species and - * @f$\hat u^0_{k_e}(T_e) = \hat h^0_{k_e}(T_e) - R T_e@f$ for electrons. - * - * @see MultiSpeciesThermo - */ - double intEnergy_mole() const override; - - void getEntropy_R(span sr) const override; - - void getGibbs_RT(span grt) const override; - - void getGibbs_ref(span g) const override; - - void getStandardVolumes_ref(span vol) const override; - - void getChemPotentials(span mu) const override; - - void getStandardChemPotentials(span muStar) const override; - - void getPartialMolarEnthalpies(span hbar) const override; - - void getPartialMolarEntropies(span sbar) const override; - - void getPartialMolarIntEnergies(span ubar) const override; - - void getParameters(AnyMap& phaseNode) const override; - - void setParameters(const AnyMap& phaseNode, - const AnyMap& rootNode=AnyMap()) override; //! Update the electron energy distribution. void updateElectronEnergyDistribution(); - //! Electron species name - string electronSpeciesName() const { - return speciesName(m_electronSpeciesIndex); - } - - //! Return the distribution Number #m_distNum + //! Return the distribution number #m_distNum. int distributionNumber() const { return m_distNum; } - //! Return the electron energy level Number #m_levelNum + //! Return the electron energy level number #m_levelNum. int levelNumber() const { return m_levelNum; } - //! Get the indicies for inelastic electron collisions + //! Get the indicies for inelastic electron collisions. //! @since New in %Cantera 3.2. const vector& kInelastic() const { return m_kInelastic; } - //! Get the indices for elastic electron collisions + //! Get the indices for elastic electron collisions. //! @since New in %Cantera 3.2. const vector& kElastic() const { return m_kElastic; } - //! target of a specific process + //! Return the target of a specific process. //! @since New in %Cantera 3.2. size_t targetIndex(size_t i) const { return m_targetSpeciesIndices[i]; } - //! Get the frequency of the applied electric field [Hz] + //! Get the frequency of the applied electric field [Hz]. //! @since New in %Cantera 3.2. double electricFieldFrequency() const { return m_electricFieldFrequency; } - //! Get the applied electric field strength [V/m] + //! Get the applied electric field strength [V/m]. double electricField() const { return m_electricField; } - //! Set the absolute electric field strength [V/m] + //! Set the absolute electric field strength [V/m]. void setElectricField(double E) { m_electricField = E; } @@ -378,20 +301,18 @@ class PlasmaPhase: public IdealGasPhase // return ne / n_total; //} - //! Get the reduced electric field strength [V·m²] + //! Get the reduced electric field strength [V·m²]. double reducedElectricField() const { return m_electricField / (molarDensity() * Avogadro); } - //! Set reduced electric field given in [V·m²] + //! Set reduced electric field given in [V·m²]. void setReducedElectricField(double EN) { m_electricField = EN * molarDensity() * Avogadro; // [V/m] } - virtual void setSolution(std::weak_ptr soln) override; - /** - * The elastic power loss (J/s/m³) + * The elastic power loss [J/s/m³] * @f[ * P_k = N_A N_A C_e e \sum_k C_k K_k, * @f] @@ -429,6 +350,140 @@ class PlasmaPhase: public IdealGasPhase void endEquilibrate() override; double intrinsicHeating() override; + //! @} + // ================================================================= // + // ================================================================= // + //! @name Molar Thermodynamic Properties of the Solution + //! @{ + + //! 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 + */ + double enthalpy_mole() const override; + + //! Return the molar entropy. Units: J/kmol/K. + /*! + * For an ideal gas mixture with an additional electron species, + * @f[ + * \hat s(T,T_e,p,X) + * = \sum_{k\neq k_e} X_k\left[\hat s_k^0(T) - R\ln\left(\frac{X_k p}{p^0}\right)\right] + * + X_{k_e}\left[\hat s_{k_e}^0(T_e) - R\ln\left(\frac{X_{k_e} p_e}{p^0}\right)\right], + * @f] + * where heavy-species properties are evaluated at @f$T@f$, electron properties at + * @f$T_e@f$, and the electron mixing term uses the electron pressure + * @f$p_e = n_{k_e} R T_e@f$. + * + * @see MultiSpeciesThermo + */ + double entropy_mole() const override; + + //! Return the molar Gibbs free energy. Units: J/kmol. + /*! + * For an ideal gas mixture with an additional electron species, + * @f[ + * \hat g(T, T_e, p, X) = \sum_k X_k \mu_k(T_k, p_k, X), + * @f] + * where heavy species use the gas temperature @f$T@f$ and bulk pressure, + * while the electron chemical potential uses @f$T_e@f$ and + * @f$p_e = n_{k_e} R T_e@f$. + * + * @see MultiSpeciesThermo + */ + double gibbs_mole() const override; + + //! Return the molar internal energy. Units: J/kmol. + /*! + * For an ideal gas mixture with an additional electron species, + * @f[ + * \hat u(T,T_e) = \sum_{k \neq k_e} X_k \hat u^0_k(T) + X_{k_e} \hat u^0_{k_e}(T_e), + * @f] + * where @f$\hat u^0_k(T) = \hat h^0_k(T) - R T@f$ for heavy species and + * @f$\hat u^0_{k_e}(T_e) = \hat h^0_{k_e}(T_e) - R T_e@f$ for electrons. + * + * @see MultiSpeciesThermo + */ + double intEnergy_mole() const override; + + double cp_mole() const override; + // double cp_mass() const; // Already defined in ThermoPhase + // double cv_mole() const; // Already defined in IdealGasPhase + + //! @} + // ================================================================= // + // ================================================================= // + //! @name Mechanical Equation of State + //! @{ + + //! Return the Gas Constant multiplied by the current electron temperature + /*! + * The units are Joules kmol-1 + */ + double RTe() const { + return electronTemperature() * GasConstant; + } + + /** + * Electron pressure. Units: Pa. + * @f[P = n_{k_e} R T_e @f] + */ + virtual double electronPressure() const { + return GasConstant * concentration(m_electronSpeciesIndex) * + electronTemperature(); + } + + //! @} + // ================================================================= // + // ================================================================= // + //! @name Partial Molar Properties of the Solution + //! @{ + + void getChemPotentials(span mu) const override; + void getPartialMolarEnthalpies(span hbar) const override; + void getPartialMolarEntropies(span sbar) const override; + void getPartialMolarIntEnergies(span ubar) const override; + // void getPartialMolarCp(span cpbar) const override; + // void getPartialMolarVolumes(span vbar) const override; + + //! @} + // ================================================================= // + // ================================================================= // + //! @name Properties of the Standard State of the Species in the Solution + //! @{ + + void getStandardChemPotentials(span muStar) const override; + // void getEnthalpy_RT(span hrt) const override; + void getEntropy_R(span sr) const override; + void getGibbs_RT(span grt) const override; + // void getIntEnergy_RT(span urt) const override; + // void getCp_R(span cpr) const override; + // void getStandardVolumes(span vol) const override; + + //! @} + // ================================================================= // + // ================================================================= // + //! @name Thermodynamic Values for the Species Reference States + //! @{ + + // void getEnthalpy_RT_ref(span hrt) const override; + // void getGibbs_RT_ref(span grt) const override; + void getGibbs_ref(span g) const override; + // void getEntropy_R_ref(span er) const override; + // void getIntEnergy_RT_ref(span urt) const override; + // void getCp_R_ref(span cprt) const override; + void getStandardVolumes_ref(span vol) const override; + + //! @} + + protected: void updateThermo() const override; diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 26f83aeb469..50f0d23f069 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -52,6 +52,167 @@ PlasmaPhase::~PlasmaPhase() } } +void PlasmaPhase::initThermo() +{ + IdealGasPhase::initThermo(); + + // Check if there is an electron species in the phase. + if (m_electronSpeciesIndex == npos) { + throw CanteraError("PlasmaPhase::initThermo", + "No electron species found."); + } +} + +void PlasmaPhase::updateThermo() const +{ + // Update the heavy species thermodynamic properties + // before updating the electron species properties. + IdealGasPhase::updateThermo(); + static const int cacheId = m_cache.getId(); + CachedScalar cached = m_cache.getScalar(cacheId); + double tempNow = temperature(); + double electronTempNow = electronTemperature(); + size_t k = m_electronSpeciesIndex; + // If the electron temperature has changed since the last time these + // properties were computed, recompute them. + if (cached.state1 != tempNow || cached.state2 != electronTempNow) { + // Evaluate the electron species thermodynamic properties + // at the electron temperature. + m_spthermo.update_single(k, electronTemperature(), + m_cp0_R[k], m_h0_RT[k], m_s0_R[k]); + cached.state1 = tempNow; + cached.state2 = electronTempNow; + + // Update the electron Gibbs functions, with the electron temperature. + m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; + } +} + +// ================================================================= // +// Overridden from IdealGasPhase or ThermoPhase // +// ================================================================= // + +bool PlasmaPhase::addSpecies(shared_ptr spec) +{ + bool added = IdealGasPhase::addSpecies(spec); + size_t k = m_kk - 1; + + if ((spec->name == "e" || spec->name == "Electron") || + (spec->composition.find("E") != spec->composition.end() && + spec->composition.size() == 1 && + spec->composition["E"] == 1)) { + if (m_electronSpeciesIndex == npos) { + m_electronSpeciesIndex = k; + } else { + throw CanteraError("PlasmaPhase::addSpecies", + "Cannot add species, {}. " + "Only one electron species is allowed.", spec->name); + } + } + return added; +} + +void PlasmaPhase::setSolution(std::weak_ptr soln) { + ThermoPhase::setSolution(soln); + // Register callback function to be executed + // when the thermo or kinetics object changed. + if (shared_ptr soln = m_soln.lock()) { + soln->registerChangedCallback(this, [&]() { + setCollisions(); + }); + } +} + +void PlasmaPhase::getParameters(AnyMap& phaseNode) const +{ + IdealGasPhase::getParameters(phaseNode); + AnyMap eedf; + eedf["type"] = m_distributionType; + vector levels(m_nPoints); + Eigen::Map(levels.data(), m_nPoints) = m_electronEnergyLevels; + eedf["energy-levels"] = levels; + if (m_distributionType == "isotropic") { + eedf["shape-factor"] = m_isotropicShapeFactor; + eedf["mean-electron-energy"].setQuantity(meanElectronEnergy(), "eV"); + } else if (m_distributionType == "discretized") { + vector dist(m_nPoints); + Eigen::Map(dist.data(), m_nPoints) = m_electronEnergyDist; + eedf["distribution"] = dist; + eedf["normalize"] = m_do_normalizeElectronEnergyDist; + } + phaseNode["electron-energy-distribution"] = std::move(eedf); +} + +void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) +{ + IdealGasPhase::setParameters(phaseNode, rootNode); + if (phaseNode.hasKey("electron-energy-distribution")) { + const AnyMap eedf = phaseNode["electron-energy-distribution"].as(); + m_distributionType = eedf["type"].asString(); + if (m_distributionType == "isotropic") { + if (eedf.hasKey("shape-factor")) { + setIsotropicShapeFactor(eedf["shape-factor"].asDouble()); + } else { + throw CanteraError("PlasmaPhase::setParameters", + "isotropic type requires shape-factor key."); + } + if (eedf.hasKey("mean-electron-energy")) { + double energy = eedf.convert("mean-electron-energy", "eV"); + setMeanElectronEnergy(energy); + } else { + throw CanteraError("PlasmaPhase::setParameters", + "isotropic type requires electron-temperature key."); + } + if (eedf.hasKey("energy-levels")) { + auto levels = eedf["energy-levels"].asVector(); + setElectronEnergyLevels(levels); + } + setIsotropicElectronEnergyDistribution(); + } else if (m_distributionType == "discretized") { + if (!eedf.hasKey("energy-levels")) { + throw CanteraError("PlasmaPhase::setParameters", + "Cannot find key energy-levels."); + } + if (!eedf.hasKey("distribution")) { + throw CanteraError("PlasmaPhase::setParameters", + "Cannot find key distribution."); + } + if (eedf.hasKey("normalize")) { + enableNormalizeElectronEnergyDist(eedf["normalize"].asBool()); + } + auto levels = eedf["energy-levels"].asVector(); + auto distribution = eedf["distribution"].asVector(levels.size()); + setDiscretizedElectronEnergyDist(levels, distribution); + } + } + + if (rootNode.hasKey("electron-collisions")) { + for (const auto& item : rootNode["electron-collisions"].asVector()) { + auto rate = make_shared(item); + Composition reactants, products; + reactants[item["target"].asString()] = 1; + reactants[electronSpeciesName()] = 1; + if (item.hasKey("product")) { + products[item["product"].asString()] = 1; + } else { + products[item["target"].asString()] = 1; + } + products[electronSpeciesName()] = 1; + if (rate->kind() == "ionization") { + products[electronSpeciesName()] += 1; + } else if (rate->kind() == "attachment") { + products[electronSpeciesName()] -= 1; + } + auto R = make_shared(reactants, products, rate); + addCollision(R); + } + } +} + +// ================================================================= // +// Electron Energy Distribution Functions // +// ================================================================= // + void PlasmaPhase::updateElectronEnergyDistribution() { if (m_distributionType == "discretized") { @@ -177,12 +338,14 @@ void PlasmaPhase::electronEnergyDistributionChanged() void PlasmaPhase::electronEnergyLevelChanged() { m_levelNum++; - // Cross sections are interpolated on the energy levels + // Cross sections are interpolated on the energy levels. if (m_collisions.size() > 0) { + vector energyLevels(m_nPoints); + MappedVector(energyLevels.data(), m_nPoints) = m_electronEnergyLevels; for (shared_ptr collision : m_collisions) { const auto& rate = boost::polymorphic_pointer_downcast (collision->rate()); - rate->updateInterpolatedCrossSection(asSpan(m_electronEnergyLevels)); + rate->updateInterpolatedCrossSection(energyLevels); } } } @@ -220,8 +383,10 @@ void PlasmaPhase::setDiscretizedElectronEnergyDist(span levels, { m_distributionType = "discretized"; m_nPoints = levels.size(); - m_electronEnergyLevels = asVectorXd(levels); - m_electronEnergyDist = asVectorXd(dist); + m_electronEnergyLevels = + Eigen::Map(levels.data(), m_nPoints); + m_electronEnergyDist = + Eigen::Map(dist.data(), m_nPoints); checkElectronEnergyLevels(); if (m_do_normalizeElectronEnergyDist) { normalizeElectronEnergyDistribution(); @@ -258,134 +423,6 @@ void PlasmaPhase::setIsotropicShapeFactor(double x) { updateElectronEnergyDistribution(); } -void PlasmaPhase::getParameters(AnyMap& phaseNode) const -{ - IdealGasPhase::getParameters(phaseNode); - AnyMap eedf; - eedf["type"] = m_distributionType; - vector levels(m_nPoints); - Eigen::Map(levels.data(), m_nPoints) = m_electronEnergyLevels; - eedf["energy-levels"] = levels; - if (m_distributionType == "isotropic") { - eedf["shape-factor"] = m_isotropicShapeFactor; - eedf["mean-electron-energy"].setQuantity(meanElectronEnergy(), "eV"); - } else if (m_distributionType == "discretized") { - vector dist(m_nPoints); - Eigen::Map(dist.data(), m_nPoints) = m_electronEnergyDist; - eedf["distribution"] = dist; - eedf["normalize"] = m_do_normalizeElectronEnergyDist; - } - phaseNode["electron-energy-distribution"] = std::move(eedf); -} - -void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) -{ - IdealGasPhase::setParameters(phaseNode, rootNode); - if (phaseNode.hasKey("electron-energy-distribution")) { - const AnyMap eedf = phaseNode["electron-energy-distribution"].as(); - m_distributionType = eedf["type"].asString(); - if (m_distributionType == "isotropic") { - if (eedf.hasKey("shape-factor")) { - setIsotropicShapeFactor(eedf["shape-factor"].asDouble()); - } else { - throw CanteraError("PlasmaPhase::setParameters", - "isotropic type requires shape-factor key."); - } - if (eedf.hasKey("mean-electron-energy")) { - double energy = eedf.convert("mean-electron-energy", "eV"); - setMeanElectronEnergy(energy); - } else { - throw CanteraError("PlasmaPhase::setParameters", - "isotropic type requires electron-temperature key."); - } - if (eedf.hasKey("energy-levels")) { - auto levels = eedf["energy-levels"].asVector(); - setElectronEnergyLevels(levels); - } - setIsotropicElectronEnergyDistribution(); - } else if (m_distributionType == "discretized") { - if (!eedf.hasKey("energy-levels")) { - throw CanteraError("PlasmaPhase::setParameters", - "Cannot find key energy-levels."); - } - if (!eedf.hasKey("distribution")) { - throw CanteraError("PlasmaPhase::setParameters", - "Cannot find key distribution."); - } - if (eedf.hasKey("normalize")) { - enableNormalizeElectronEnergyDist(eedf["normalize"].asBool()); - } - auto levels = eedf["energy-levels"].asVector(); - auto distribution = eedf["distribution"].asVector(levels.size()); - setDiscretizedElectronEnergyDist(levels, distribution); - } - } - - if (rootNode.hasKey("electron-collisions")) { - for (const auto& item : rootNode["electron-collisions"].asVector()) { - auto rate = make_shared(item); - Composition reactants, products; - reactants[item["target"].asString()] = 1; - reactants[electronSpeciesName()] = 1; - if (item.hasKey("product")) { - products[item["product"].asString()] = 1; - } else { - products[item["target"].asString()] = 1; - } - products[electronSpeciesName()] = 1; - if (rate->kind() == "ionization") { - products[electronSpeciesName()] += 1; - } else if (rate->kind() == "attachment") { - products[electronSpeciesName()] -= 1; - } - auto R = make_shared(reactants, products, rate); - addCollision(R); - } - } -} - -bool PlasmaPhase::addSpecies(shared_ptr spec) -{ - bool added = IdealGasPhase::addSpecies(spec); - size_t k = m_kk - 1; - - if ((spec->name == "e" || spec->name == "Electron") || - (spec->composition.find("E") != spec->composition.end() && - spec->composition.size() == 1 && - spec->composition["E"] == 1)) { - if (m_electronSpeciesIndex == npos) { - m_electronSpeciesIndex = k; - } else { - throw CanteraError("PlasmaPhase::addSpecies", - "Cannot add species, {}. " - "Only one electron species is allowed.", spec->name); - } - } - return added; -} - -void PlasmaPhase::initThermo() -{ - IdealGasPhase::initThermo(); - - // Check electron species - if (m_electronSpeciesIndex == npos) { - throw CanteraError("PlasmaPhase::initThermo", - "No electron species found."); - } -} - -void PlasmaPhase::setSolution(std::weak_ptr soln) { - ThermoPhase::setSolution(soln); - // register callback function to be executed - // when the thermo or kinetics object changed - if (shared_ptr soln = m_soln.lock()) { - soln->registerChangedCallback(this, [&]() { - setCollisions(); - }); - } -} - void PlasmaPhase::setCollisions() { if (shared_ptr soln = m_soln.lock()) { @@ -601,26 +638,55 @@ double PlasmaPhase::elasticPowerLoss() return q_elastic; } -void PlasmaPhase::updateThermo() const +double PlasmaPhase::electronMobility() const { - IdealGasPhase::updateThermo(); - static const int cacheId = m_cache.getId(); - CachedScalar cached = m_cache.getScalar(cacheId); - double tempNow = temperature(); - double electronTempNow = electronTemperature(); - size_t k = m_electronSpeciesIndex; - // If the electron temperature has changed since the last time these - // properties were computed, recompute them. - if (cached.state1 != tempNow || cached.state2 != electronTempNow) { - m_spthermo.update_single(k, electronTemperature(), - m_cp0_R[k], m_h0_RT[k], m_s0_R[k]); - cached.state1 = tempNow; - cached.state2 = electronTempNow; + // Only implemented when using the Boltzmann two-term EEDF + if (m_distributionType == "Boltzmann-two-term") { + return m_eedfSolver->getElectronMobility(); + } else { + throw NotImplementedError("PlasmaPhase::electronMobility", + "Electron mobility is only available for 'Boltzmann-two-term' " + "electron energy distributions."); } - // update the species Gibbs functions - m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; } + +double PlasmaPhase::jouleHeatingPower() const +{ + // sigma = e * n_e * mu_e [S/m]; q_J = sigma * E^2 [W/m^3] + const double mu_e = electronMobility(); // m^2 / (V·s) + if (mu_e <= 0.0) { + return 0.0; + } + const double ne = concentration(m_electronSpeciesIndex) * Avogadro; // m^-3 + if (ne <= 0.0) { + return 0.0; + } + const double E = electricField(); // V/m + if (E <= 0.0) { + return 0.0; + } + const double sigma = ElectronCharge * ne * mu_e; // S/m + return sigma * E * E; // W/m^3 +} + +double PlasmaPhase::intrinsicHeating() +{ + // Joule heating: sigma * E^2 [W/m^3] + const double qJ = jouleHeatingPower(); + checkFinite(qJ); + + // Elastic + inelastic recoil power loss [W/m^3] + double qElastic = elasticPowerLoss(); + checkFinite(qElastic); + + return qJ + qElastic; +} + +// ================================================================= // +// Molar Thermodynamic Properties of the Solution // +// ================================================================= // + double PlasmaPhase::enthalpy_mole() const { double value = IdealGasPhase::enthalpy_mole(); value += GasConstant * (electronTemperature() - temperature()) * @@ -662,16 +728,16 @@ double PlasmaPhase::gibbs_mole() const return g; } -void PlasmaPhase::getGibbs_ref(span g) const -{ - IdealGasPhase::getGibbs_ref(g); - g[m_electronSpeciesIndex] *= electronTemperature() / temperature(); -} +// ================================================================= // +// Partial Molar Properties of the Solution // +// ================================================================= // -void PlasmaPhase::getStandardVolumes_ref(span vol) const +void PlasmaPhase::getChemPotentials(span mu) const { - IdealGasPhase::getStandardVolumes_ref(vol); - vol[m_electronSpeciesIndex] *= electronTemperature() / temperature(); + IdealGasPhase::getChemPotentials(mu); + size_t k = m_electronSpeciesIndex; + double xx = std::max(SmallNumber, moleFraction(k)); + mu[k] += (RTe() - RT()) * log(xx); } void PlasmaPhase::getPartialMolarEnthalpies(span hbar) const @@ -695,17 +761,14 @@ void PlasmaPhase::getPartialMolarIntEnergies(span ubar) const for (size_t k = 0; k < m_kk; k++) { ubar[k] = RT() * (_h[k] - 1.0); } + // Redefine it for the electron species. size_t k = m_electronSpeciesIndex; ubar[k] = RTe() * (_h[k] - 1.0); } -void PlasmaPhase::getChemPotentials(span mu) const -{ - IdealGasPhase::getChemPotentials(mu); - size_t k = m_electronSpeciesIndex; - double xx = std::max(SmallNumber, moleFraction(k)); - mu[k] += (RTe() - RT()) * log(xx); -} +// ================================================================= // +// Properties of the Standard State of the Species in the Solution // +// ================================================================= // void PlasmaPhase::getStandardChemPotentials(span muStar) const { @@ -745,49 +808,22 @@ void PlasmaPhase::getGibbs_RT(span grt) const } } -double PlasmaPhase::electronMobility() const -{ - // Only implemented when using the Boltzmann two-term EEDF - if (m_distributionType == "Boltzmann-two-term") { - return m_eedfSolver->getElectronMobility(); - } else { - throw NotImplementedError("PlasmaPhase::electronMobility", - "Electron mobility is only available for 'Boltzmann-two-term' " - "electron energy distributions."); - } -} -double PlasmaPhase::jouleHeatingPower() const +// ================================================================= // +// Thermodynamic Values for the Species Reference States // +// ================================================================= // + +void PlasmaPhase::getGibbs_ref(span g) const { - // sigma = e * n_e * mu_e [S/m]; q_J = sigma * E^2 [W/m^3] - const double mu_e = electronMobility(); // m^2 / (V·s) - if (mu_e <= 0.0) { - return 0.0; - } - const double ne = concentration(m_electronSpeciesIndex) * Avogadro; // m^-3 - if (ne <= 0.0) { - return 0.0; - } - const double E = electricField(); // V/m - if (E <= 0.0) { - return 0.0; - } - const double sigma = ElectronCharge * ne * mu_e; // S/m - return sigma * E * E; // W/m^3 + IdealGasPhase::getGibbs_ref(g); + g[m_electronSpeciesIndex] *= electronTemperature() / temperature(); } -double PlasmaPhase::intrinsicHeating() +void PlasmaPhase::getStandardVolumes_ref(span vol) const { - // Joule heating: sigma * E^2 [W/m^3] - const double qJ = jouleHeatingPower(); - checkFinite(qJ); - - // Elastic + inelastic recoil power loss [W/m^3] - double qElastic = elasticPowerLoss(); - checkFinite(qElastic); - - return qJ + qElastic; + IdealGasPhase::getStandardVolumes_ref(vol); + vol[m_electronSpeciesIndex] *= electronTemperature() / temperature(); } From ce860be9ead86bd1cc28e61709fadda231b4613e Mon Sep 17 00:00:00 2001 From: pag1pag <66677348+pag1pag@users.noreply.github.com> Date: Tue, 26 May 2026 11:21:48 +0200 Subject: [PATCH 02/11] [Plasma] Update `enthalpy_mole` for consistency with other methods --- src/thermo/PlasmaPhase.cpp | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 50f0d23f069..aca59481480 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -687,12 +687,15 @@ double PlasmaPhase::intrinsicHeating() // Molar Thermodynamic Properties of the Solution // // ================================================================= // -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::enthalpy_mole() const +{ + m_work.resize(m_kk); + getPartialMolarEnthalpies(m_work); + double h = 0.0; + for (size_t k = 0; k < m_kk; ++k) { + h += moleFraction(k) * m_work[k]; + } + return h; } double PlasmaPhase::intEnergy_mole() const From 164b1777e081f4cac4259973dc8fb6ef693f57b6 Mon Sep 17 00:00:00 2001 From: pag1pag <66677348+pag1pag@users.noreply.github.com> Date: Tue, 26 May 2026 09:48:28 +0200 Subject: [PATCH 03/11] [Plasma] Add `pressure` and `meanTemperature` --- include/cantera/thermo/PlasmaPhase.h | 56 ++++++++++++++++++++++++++++ src/thermo/PlasmaPhase.cpp | 16 ++++++++ 2 files changed, 72 insertions(+) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 67fd11adc8a..e07d1edc6f3 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -423,6 +423,62 @@ class PlasmaPhase: public IdealGasPhase //! @name Mechanical Equation of State //! @{ + //! Return the mean temperature of the plasma phase [K]. + /*! + * In a plasma phase, the electron temperature can differ from the + * heavy-species (gas) temperature. Therefore, the mean temperature is + * defined as a mole-fraction-weighted average of the electron and + * heavy-species temperatures: + * @f[ + * \bar{T} = \sum_{k \neq k_e} X_k T_g + X_{k_e} T_e + * = (1 - X_{k_e}) T_g + X_{k_e} T_e + * = T_g + X_{k_e} (T_e - T_g) + * @f] + * See the `pressure()` method for usage of the mean temperature in + * calculating the pressure of the plasma phase. + */ + double meanTemperature() const; + + //! Return the pressure of the plasma phase [Pa]. + /*! + * The pressure of the plasma phase is calculated using the mean temperature, + * which is a mole-fraction-weighted average of the electron and heavy-species + * temperatures. + * @f[ + * P = \sum_k n_k k_B T_k + * = \sum_{k \neq e} n_k k_B T_g + n_{e} k_B T_e + * = (n_{total} - n_{e}) k_B T_g + n_{e} k_B T_e + * = n_{total} (1 - X_e) k_B T_g + n_{total} X_e k_B T_e + * = n_{total} k_B \bar{T} + * @f] + * where @f$ \bar{T} @f$ is the mean temperature of the plasma phase, + * defined in the `meanTemperature()` method. + * Here, @f$ n_k @f$ is the number density of species @f$ k @f$, + * @f$ n_{total} @f$ is the total number density of the phase, + * @f$ X_e @f$ is the mole fraction of electrons, @f$ T_g @f$ is the + * heavy-species (gas) temperature, @f$ T_e @f$ is the electron temperature, + * and @f$ k_B @f$ is the Boltzmann constant. + * + * The number density times Boltzmann constant can be expressed as + * @f$ n_{total} k_B = C_{total} R @f$, where @f$ C_{total} @f$ is the + * total molar concentration of the phase and @f$ R @f$ is the gas constant. + */ + double pressure() const override; + + //! Set the pressure at constant temperature and composition. + /*! + * Units: Pa. + * This method is implemented by setting the mass density to + * @f[ + * \rho = \frac{P \overline{W}}{\hat R \overline{T} }. + * @f] + * + * @param p Pressure (Pa) + */ + void setPressure(double p) override { + setDensity(p * meanMolecularWeight() / (GasConstant * meanTemperature())); + } + //! Return the Gas Constant multiplied by the current electron temperature /*! * The units are Joules kmol-1 diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index aca59481480..6cc55518ea6 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -731,6 +731,22 @@ double PlasmaPhase::gibbs_mole() const return g; } +// ================================================================= // +// Mechanical Equation of State // +// ================================================================= // + +double PlasmaPhase::meanTemperature() const +{ + double T_g = temperature(); + double T_e = electronTemperature(); + double X_e = moleFraction(m_electronSpeciesIndex); + return T_g + X_e * (T_e - T_g); +} + +double PlasmaPhase::pressure() const { + return GasConstant * meanTemperature() * density() / meanMolecularWeight(); +} + // ================================================================= // // Partial Molar Properties of the Solution // // ================================================================= // From 8e1cdbda6e21c021cc322e5b069f3d83aa643645 Mon Sep 17 00:00:00 2001 From: pag1pag <66677348+pag1pag@users.noreply.github.com> Date: Tue, 26 May 2026 10:18:02 +0200 Subject: [PATCH 04/11] [Plasma] Make `mean_temperature` a Python method --- interfaces/cython/cantera/composite.pyi | 2 ++ interfaces/cython/cantera/thermo.pxd | 1 + interfaces/cython/cantera/thermo.pyi | 2 ++ interfaces/cython/cantera/thermo.pyx | 9 +++++++++ 4 files changed, 14 insertions(+) diff --git a/interfaces/cython/cantera/composite.pyi b/interfaces/cython/cantera/composite.pyi index 1599e6cc414..1b2ac51825c 100644 --- a/interfaces/cython/cantera/composite.pyi +++ b/interfaces/cython/cantera/composite.pyi @@ -747,6 +747,8 @@ class Quantity: def electron_species_name(self) -> str: ... @property def elastic_power_loss(self) -> float: ... + @property + def mean_temperature(self) -> float: ... class SolutionArray(SolutionArrayBase, Generic[_P]): _phase: _P | None diff --git a/interfaces/cython/cantera/thermo.pxd b/interfaces/cython/cantera/thermo.pxd index ef8be3b98e1..9ed25467f21 100644 --- a/interfaces/cython/cantera/thermo.pxd +++ b/interfaces/cython/cantera/thermo.pxd @@ -215,6 +215,7 @@ cdef extern from "cantera/thermo/PlasmaPhase.h": double electricField() void updateElectronEnergyDistribution() except +translate_exception double elasticPowerLoss() except +translate_exception + double meanTemperature() cdef extern from "cantera/cython/thermo_utils.h": diff --git a/interfaces/cython/cantera/thermo.pyi b/interfaces/cython/cantera/thermo.pyi index 799cc4155e4..84f2803d3c4 100644 --- a/interfaces/cython/cantera/thermo.pyi +++ b/interfaces/cython/cantera/thermo.pyi @@ -525,6 +525,8 @@ class ThermoPhase(_SolutionBase): def electron_species_name(self) -> str: ... @property def elastic_power_loss(self) -> float: ... + @property + def mean_temperature(self) -> float: ... class InterfacePhase(ThermoPhase): @property diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx index d504e44f303..c4418c84bfb 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/interfaces/cython/cantera/thermo.pyx @@ -1889,6 +1889,15 @@ cdef class ThermoPhase(_SolutionBase): span[double](&data_levels[0], data_levels.size), span[double](&data_dist[0], data_dist.size)) + property mean_temperature: + """Get the mean temperature of the plasma [K]. + This is a mole-weighted average of the electron and heavy species temperatures. + """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.meanTemperature() + def update_electron_energy_distribution(self): """ Update the electron energy distribution function to account for changes in From f0c0ae3699fb0d83ebce2e1cf35f1250feb83313 Mon Sep 17 00:00:00 2001 From: pag1pag <66677348+pag1pag@users.noreply.github.com> Date: Tue, 26 May 2026 10:37:32 +0200 Subject: [PATCH 05/11] [Plasma] Update `getStandardChemPotentials` and comment `getEntropy_R`, `getGibbs_RT`, and `getPartialMolarEntropies` --- include/cantera/thermo/PlasmaPhase.h | 27 +++++++++++++-- src/thermo/PlasmaPhase.cpp | 49 +++++++--------------------- 2 files changed, 35 insertions(+), 41 deletions(-) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index e07d1edc6f3..2945e895c73 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -504,7 +504,7 @@ class PlasmaPhase: public IdealGasPhase void getChemPotentials(span mu) const override; void getPartialMolarEnthalpies(span hbar) const override; - void getPartialMolarEntropies(span sbar) const override; + // void getPartialMolarEntropies(span sbar) const override; void getPartialMolarIntEnergies(span ubar) const override; // void getPartialMolarCp(span cpbar) const override; // void getPartialMolarVolumes(span vbar) const override; @@ -515,10 +515,20 @@ class PlasmaPhase: public IdealGasPhase //! @name Properties of the Standard State of the Species in the Solution //! @{ + //! Return the standard chemical potentials of the species [J/kmol]. + /*! + * For heavy species, this is identical to the IdealGasPhase + * implementation. For electrons, the standard chemical potential + * is evaluated at the electron temperature: + * @f[ + * \mu^0_{e}(T_e) = g_k^0(T_e) + RT_e \ln \left(\frac{P}{P^0}\right). + * @f] + */ void getStandardChemPotentials(span muStar) const override; + // void getEnthalpy_RT(span hrt) const override; - void getEntropy_R(span sr) const override; - void getGibbs_RT(span grt) const override; + // void getEntropy_R(span sr) const override; + // void getGibbs_RT(span grt) const override; // void getIntEnergy_RT(span urt) const override; // void getCp_R(span cpr) const override; // void getStandardVolumes(span vol) const override; @@ -531,7 +541,18 @@ class PlasmaPhase: public IdealGasPhase // void getEnthalpy_RT_ref(span hrt) const override; // void getGibbs_RT_ref(span grt) const override; + + //! Return the reference-state Gibbs free energys of the species [J/kmol]. + /*! + * For heavy species, this is identical to the IdealGasPhase + * implementation. For electrons, the reference-state Gibbs free energy + * is evaluated at the electron temperature: + * @f[ + * \hat{g}^0_{e}(T_e) = \hat{h}^0_{e}(T_e) - T_e \hat{s}^0_{e}(T_e). + * @f] + */ void getGibbs_ref(span g) const override; + // void getEntropy_R_ref(span er) const override; // void getIntEnergy_RT_ref(span urt) const override; // void getCp_R_ref(span cprt) const override; diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 6cc55518ea6..583c41f770e 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -765,13 +765,6 @@ void PlasmaPhase::getPartialMolarEnthalpies(span hbar) const hbar[m_electronSpeciesIndex] *= electronTemperature() / temperature(); } -void PlasmaPhase::getPartialMolarEntropies(span sbar) const -{ - IdealGasPhase::getPartialMolarEntropies(sbar); - double logp = log(pressure()); - double logpe = log(electronPressure()); - sbar[m_electronSpeciesIndex] += GasConstant * (logp - logpe); -} void PlasmaPhase::getPartialMolarIntEnergies(span ubar) const { @@ -791,40 +784,20 @@ void PlasmaPhase::getPartialMolarIntEnergies(span ubar) const void PlasmaPhase::getStandardChemPotentials(span muStar) const { - IdealGasPhase::getStandardChemPotentials(muStar); - size_t k = m_electronSpeciesIndex; - muStar[k] -= log(pressure() / refPressure()) * RT(); - muStar[k] += log(electronPressure() / refPressure()) * RTe(); -} + // After calling PlasmaPhase::getGibbs_ref, muStar = g_k(T_k). + // g_k is evaluated at T for heavy species and at Te for electrons. + getGibbs_ref(muStar); -void PlasmaPhase::getEntropy_R(span sr) const -{ - checkArraySize("PlasmaPhase::getEntropy_R", sr.size(), m_kk); - auto _s = entropy_R_ref(); - copy(_s.begin(), _s.end(), sr.begin()); - double tmp = log(pressure() / refPressure()); + // Then, we need to add R*T_k*ln(P/Pref) to g_k. + // .. For heavy species, mu_star = g_k(T) + R*T*ln(P/Pref) + double tmp = log(pressure() / refPressure()) * RT(); for (size_t k = 0; k < m_kk; k++) { - if (k != m_electronSpeciesIndex) { - sr[k] -= tmp; - } else { - sr[k] -= log(electronPressure() / refPressure()); - } - } -} - -void PlasmaPhase::getGibbs_RT(span grt) const -{ - checkArraySize("PlasmaPhase::getGibbs_RT", grt.size(), m_kk); - auto gibbsrt = gibbs_RT_ref(); - copy(gibbsrt.begin(), gibbsrt.end(), grt.begin()); - double tmp = log(pressure() / refPressure()); - for (size_t k = 0; k < m_kk; k++) { - if (k != m_electronSpeciesIndex) { - grt[k] += tmp; - } else { - grt[k] += log(electronPressure() / refPressure()); - } + muStar[k] += tmp; } + // .. For electrons, mu_star = g_k(Te) + R*T_e*ln(P/Pref) + size_t k = m_electronSpeciesIndex; + muStar[k] -= log(pressure() / refPressure()) * RT(); + muStar[k] += log(pressure() / refPressure()) * RTe(); } From a03f3d2fec55a11c794bd87b49047a3dc9bd4eb0 Mon Sep 17 00:00:00 2001 From: pag1pag <66677348+pag1pag@users.noreply.github.com> Date: Tue, 26 May 2026 11:16:01 +0200 Subject: [PATCH 06/11] [Plasma] Add `standardConcentration`, `getPartialMolarVolumes` and `getStandardVolumes` --- include/cantera/thermo/PlasmaPhase.h | 45 ++++++++++++++++++++++++++-- src/thermo/PlasmaPhase.cpp | 29 +++++++++++++++++- 2 files changed, 71 insertions(+), 3 deletions(-) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 2945e895c73..e5e4c710762 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -496,6 +496,39 @@ class PlasmaPhase: public IdealGasPhase electronTemperature(); } + double thermalExpansionCoeff() const override { + // Which temperature to use here? + throw NotImplementedError("PlasmaPhase::thermalExpansionCoeff"); + } + + double soundSpeed() const override { + // Which temperature to use here? + throw NotImplementedError("PlasmaPhase::soundSpeed"); + } + + //! @} + // ================================================================= // + // ================================================================= // + //! @name Chemical Potentials and Activities + //! @{ + + //! Returns the standard concentration @f$ C^0_k @f$, which is used to + //! normalize the generalized concentration. + /*! + * This is defined as the concentration by which the generalized + * concentration is normalized to produce the activity. Since the activity + * for an ideal gas mixture is simply the mole fraction, for an ideal gas + * @f$ C^0_k = P/\hat R T @f$. + * For a multi-temperature plasma phase, the mean temperature is used instead: + * @f[$ C^0_k = P/\hat R \overline{T} @f$. + * + * @param k Parameter indicating the species. The default + * is to assume this refers to species 0. + * @return + * Returns the standard Concentration in units of m3 kmol-1. + */ + double standardConcentration(size_t k=0) const override; + //! @} // ================================================================= // // ================================================================= // @@ -507,7 +540,11 @@ class PlasmaPhase: public IdealGasPhase // void getPartialMolarEntropies(span sbar) const override; void getPartialMolarIntEnergies(span ubar) const override; // void getPartialMolarCp(span cpbar) const override; - // void getPartialMolarVolumes(span vbar) const override; + + // Whenever a temperature can be defined, the following relation holds: + // h_k = u_k + R * T_k = u_k + p * v_k + // Therefore, v_k = R * T_k / p + void getPartialMolarVolumes(span vbar) const override; //! @} // ================================================================= // @@ -531,7 +568,11 @@ class PlasmaPhase: public IdealGasPhase // void getGibbs_RT(span grt) const override; // void getIntEnergy_RT(span urt) const override; // void getCp_R(span cpr) const override; - // void getStandardVolumes(span vol) const override; + + // Whenever a temperature can be defined, the following relation holds: + // h_k = u_k + R * T_k = u_k + p * v_k + // Therefore, v_k = R * T_k / p + void getStandardVolumes(span vol) const override; //! @} // ================================================================= // diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 583c41f770e..17526877171 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -747,6 +747,17 @@ double PlasmaPhase::pressure() const { return GasConstant * meanTemperature() * density() / meanMolecularWeight(); } + +// ================================================================= // +// Chemical Potentials and Activities // +// ================================================================= // + +double PlasmaPhase::standardConcentration(size_t k) const +{ + return pressure() / (GasConstant * meanTemperature()); +} + + // ================================================================= // // Partial Molar Properties of the Solution // // ================================================================= // @@ -778,6 +789,15 @@ void PlasmaPhase::getPartialMolarIntEnergies(span ubar) const ubar[k] = RTe() * (_h[k] - 1.0); } +void PlasmaPhase::getPartialMolarVolumes(span vbar) const +{ + double vol = RT() / pressure(); + for (size_t k = 0; k < m_kk; k++) { + vbar[k] = vol; + } + vbar[m_electronSpeciesIndex] = RTe() / pressure(); +} + // ================================================================= // // Properties of the Standard State of the Species in the Solution // // ================================================================= // @@ -800,7 +820,14 @@ void PlasmaPhase::getStandardChemPotentials(span muStar) const muStar[k] += log(pressure() / refPressure()) * RTe(); } - +void PlasmaPhase::getStandardVolumes(span vol) const +{ + double tmp = RT() / pressure(); + for (size_t k = 0; k < m_kk; k++) { + vol[k] = tmp; + } + vol[m_electronSpeciesIndex] = RTe() / pressure(); +} // ================================================================= // // Thermodynamic Values for the Species Reference States // From 31ef6fa27017b2a3db9f5c4e0539ff0ec17f777a Mon Sep 17 00:00:00 2001 From: pag1pag <66677348+pag1pag@users.noreply.github.com> Date: Mon, 15 Dec 2025 16:03:53 +0100 Subject: [PATCH 07/11] [Phase] Make `setState_TD` virtual in Phase class Changed the setState_TD method in the Phase class to be virtual, allowing derived classes to override this method for custom behavior. --- include/cantera/thermo/Phase.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h index 7e32209e2d6..d5fb54f27cc 100644 --- a/include/cantera/thermo/Phase.h +++ b/include/cantera/thermo/Phase.h @@ -378,7 +378,7 @@ class Phase //! @param t Temperature in kelvin //! @param rho Density (kg/m^3) //! @since New in %Cantera 3.0. - void setState_TD(double t, double rho); + virtual void setState_TD(double t, double rho); //! @} end group set thermo state From 3c7294cfaa072f532a6252b5ce56ecc04830fea2 Mon Sep 17 00:00:00 2001 From: pag1pag <66677348+pag1pag@users.noreply.github.com> Date: Tue, 26 May 2026 11:16:58 +0200 Subject: [PATCH 08/11] [Plasma] Add `setState` and friends --- include/cantera/thermo/PlasmaPhase.h | 56 +++++++++ src/thermo/PlasmaPhase.cpp | 169 +++++++++++++++++++++++++++ 2 files changed, 225 insertions(+) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index e5e4c710762..f6e33defa36 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -600,7 +600,63 @@ class PlasmaPhase: public IdealGasPhase void getStandardVolumes_ref(span vol) const override; //! @} + // ================================================================= // + // ================================================================= // + //! @name Setting the State + //! + //! For a plasma phase, setting the state requires specifying both + //! the heavy-species (gas) temperature and the electron temperature. + //! @{ + + //! Set the state using an AnyMap containing any combination of properties + //! supported by the thermodynamic model + /*! + * Accepted keys are: + * * `X` (mole fractions) + * * `Y` (mass fractions) + * * `T` or `Tg` or `gas-temperature` [K] + * * `Te` or `electron-temperature` [K] + * * `P` or `pressure` [Pa] + * * `D` or `density` [kg/m^3] + * + * Composition can be specified as either an AnyMap of species names to + * values or as a composition string. All other values can be given as + * floating point values in Cantera's default units, or as strings with the + * units specified, which will be converted using the Units class. + */ + void setState(const AnyMap& state) override; + + //! Set the gas and electron temperature [K] and pressure [Pa]. + /*! + * @param t Temperature [K] + * @param p Pressure [Pa] + */ + void setState_TP(double t, double p) override; + //! Set the gas temperature [K], electron temperature [K], and pressure [Pa]. + /*! + * @param Tg Gas (heavy-species) temperature [K] + * @param Te Electron temperature [K] + * @param p Pressure [Pa] + */ + virtual void setState_TgTeP(double Tg, double Te, double p); + + //! Set the gas and electron temperature [K] and mass density [kg/m^3]. + /*! + * @param t Temperature [K] + * @param rho Mass density [kg/m^3] + */ + void setState_TD(double t, double rho) override; + + //! Set the gas temperature [K], electron temperature [K], and mass density [kg/m^3]. + /*! + * @param Tg Gas (heavy-species) temperature [K] + * @param Te Electron temperature [K] + * @param rho Mass density [kg/m^3] + */ + virtual void setState_TgTeD(double Tg, double Te, double rho); + + //! @} protected: diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 17526877171..245a3a4b4f0 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -846,4 +846,173 @@ void PlasmaPhase::getStandardVolumes_ref(span vol) const } +// ================================================================= // +// Setting the State // +// ================================================================= // + + +void PlasmaPhase::setState(const AnyMap& input_state) +{ + AnyMap state = input_state; + + // Remap allowable synonyms. + if (state.hasKey("mass-fractions")) { + state["Y"] = state["mass-fractions"]; + state.erase("mass-fractions"); + } + if (state.hasKey("mole-fractions")) { + state["X"] = state["mole-fractions"]; + state.erase("mole-fractions"); + } + if (state.hasKey("gas-temperature")) { + state["T"] = state["gas-temperature"]; + } + if (state.hasKey("Tg")) { + state["T"] = state["Tg"]; + } + if (state.hasKey("electron-temperature")) { + state["Te"] = state["electron-temperature"]; + } + if (state.hasKey("pressure")) { + state["P"] = state["pressure"]; + } + if (state.hasKey("density")) { + state["D"] = state["density"]; + } + + // Set composition. + if (state.hasKey("X")) { + if (state["X"].is()) { + setMoleFractionsByName(state["X"].asString()); + } else { + setMoleFractionsByName(state["X"].asMap()); + } + state.erase("X"); + } else if (state.hasKey("Y")) { + if (state["Y"].is()) { + setMassFractionsByName(state["Y"].asString()); + } else { + setMassFractionsByName(state["Y"].asMap()); + } + state.erase("Y"); + } + + // Set thermodynamic state using whichever property set is found + if (state.hasKey("T") && state.hasKey("Te") && state.hasKey("P")) { + setState_TgTeP( + state.convert("T", "K"), + state.convert("Te", "K"), + state.convert("P", "Pa") + ); + } else if (state.hasKey("T") && state.hasKey("Te") && state.hasKey("D")) { + setState_TgTeD( + state.convert("T", "K"), + state.convert("Te", "K"), + state.convert("D", "kg/m^3") + ); + } else if (state.hasKey("T") && state.hasKey("P")) { + setState_TP( + state.convert("T", "K"), + state.convert("P", "Pa") + ); + } else if (state.hasKey("T") && state.hasKey("D")) { + setState_TD( + state.convert("T", "K"), + state.convert("D", "kg/m^3") + ); + } else { + throw CanteraError("PlasmaPhase::setState", + "'state' did not specify a recognized set of properties.\n" + "Keys provided were: {}", input_state.keys_str()); + } +} + +void PlasmaPhase::setState_TP(double t, double p) +{ + vector state(partialStateSize()); + savePartialState(state); + try { + setTemperature(t); + if (m_distributionType == "isotropic") { + setElectronTemperature(t); + } else { + // Warning for other distribution types. + warn_user("PlasmaPhase::setState_TP", + "The electron temperature is not set equal to the gas temperature " + "when using '{}' electron energy distribution. " + "This may not be intended.", m_distributionType); + } + setPressure(p); + } catch (std::exception&) { + restorePartialState(state); + throw; + } +} + +void PlasmaPhase::setState_TgTeP(double Tg, double Te, double p) +{ + vector state(partialStateSize()); + savePartialState(state); + try { + setTemperature(Tg); + if (m_distributionType == "isotropic") { + setElectronTemperature(Te); + } else { + // Warning for other distribution types. + warn_user("PlasmaPhase::setState_TgTeP", + "The electron temperature cannot be set " + "when using '{}' electron energy distribution. " + "This may not be intended.", m_distributionType); + } + setPressure(p); + } catch (std::exception&) { + restorePartialState(state); + throw; + } +} + +void PlasmaPhase::setState_TD(double t, double rho) +{ + vector state(partialStateSize()); + savePartialState(state); + try { + setTemperature(t); + if (m_distributionType == "isotropic") { + setElectronTemperature(t); + } else { + // Warning for other distribution types. + warn_user("PlasmaPhase::setState_TD", + "The electron temperature is not set equal to the gas temperature " + "when using '{}' electron energy distribution. " + "This may not be intended.", m_distributionType); + } + setDensity(rho); + } catch (std::exception&) { + restorePartialState(state); + throw; + } +} + +void PlasmaPhase::setState_TgTeD(double Tg, double Te, double rho) +{ + vector state(partialStateSize()); + savePartialState(state); + try { + setTemperature(Tg); + if (m_distributionType == "isotropic") { + setElectronTemperature(Te); + } else { + // Warning for other distribution types. + warn_user("PlasmaPhase::setState_TgTeD", + "The electron temperature cannot be set " + "when using '{}' electron energy distribution. " + "This may not be intended.", m_distributionType); + } + setDensity(rho); + } catch (std::exception&) { + restorePartialState(state); + throw; + } +} + } From fa1156f7e7331a0f15e7583f9069a00addb23b87 Mon Sep 17 00:00:00 2001 From: pag1pag <66677348+pag1pag@users.noreply.github.com> Date: Thu, 19 Mar 2026 18:26:09 +0100 Subject: [PATCH 09/11] [Plasma] Check electron temperature is positive --- src/thermo/PlasmaPhase.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 245a3a4b4f0..0fc61543599 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -289,6 +289,11 @@ void PlasmaPhase::setIsotropicElectronEnergyDistribution() } void PlasmaPhase::setElectronTemperature(const double Te) { + // If the electron temperature is negative, throw an error. + if (Te < 0.0) { + throw CanteraError("PlasmaPhase::setElectronTemperature", + "Electron temperature cannot be negative."); + } m_electronTemp = Te; updateElectronEnergyDistribution(); } From 47dd9e129f38bd002c7eef9bf70447bb87075d2c Mon Sep 17 00:00:00 2001 From: pag1pag <66677348+pag1pag@users.noreply.github.com> Date: Tue, 26 May 2026 11:20:34 +0200 Subject: [PATCH 10/11] [Plasma] Update docstrings and comments --- include/cantera/thermo/PlasmaPhase.h | 432 ++++++++++++++++++--------- src/thermo/PlasmaPhase.cpp | 84 ++++-- 2 files changed, 339 insertions(+), 177 deletions(-) diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index f6e33defa36..3e456a6b238 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -83,11 +83,11 @@ class ElectronCollisionPlasmaRate; class PlasmaPhase: public IdealGasPhase { public: - //! Construct and initialize a PlasmaPhase object - //! directly from an input file. The constructor initializes the electron - //! energy distribution to be Druyvesteyn distribution (m_x = 2.0). The initial - //! electron energy grid is set to a linear space which starts at 0.01 eV and ends - //! at 1 eV with 1000 points. + //! Construct and initialize a PlasmaPhase object directly from an input file. + //! The constructor initializes the electron energy distribution to be a + //! Maxwellian distribution (m_isotropicShapeFactor = 1.0). + //! The initial electron energy grid is set to a linear space which starts + //! at 0.01 eV and ends at 1 eV with 1000 points. /*! * @param inputFile Name of the input file containing the phase definition * to set up the object. If blank, an empty phase will be @@ -137,18 +137,18 @@ class PlasmaPhase: public IdealGasPhase //! @{ //! Set electron energy levels. - //! @param levels The vector of electron energy levels (eV). + //! @param levels The vector of electron energy levels [eV]. //! Length: #m_nPoints. void setElectronEnergyLevels(span levels); //! Get electron energy levels. - //! @param levels The vector of electron energy levels (eV). Length: #m_nPoints + //! @param levels The vector of electron energy levels [eV]. Length: #m_nPoints void getElectronEnergyLevels(span levels) const { Eigen::Map(levels.data(), levels.size()) = m_electronEnergyLevels; } //! Set discretized electron energy distribution. - //! @param levels The vector of electron energy levels (eV). + //! @param levels The vector of electron energy levels [eV]. //! Length: #m_nPoints. //! @param distrb The vector of electron energy distribution. //! Length: #m_nPoints. @@ -168,25 +168,35 @@ class PlasmaPhase: public IdealGasPhase //! @param x The shape factor void setIsotropicShapeFactor(double x); - //! The shape factor of isotropic electron energy distribution + //! The shape factor of isotropic electron energy distribution. double isotropicShapeFactor() const { return m_isotropicShapeFactor; } - //! Set the internally stored electron temperature of the phase (K). - //! @param Te Electron temperature in Kelvin + //! Electron Temperature [K]. + //! @return The electron temperature of the phase. + double electronTemperature() const override { + return m_electronTemp; + } + //! Set the internally stored electron temperature of the phase [K]. + //! @param Te Electron temperature in Kelvin. void setElectronTemperature(double Te) override; - //! Set mean electron energy [eV]. This method also sets electron temperature - //! accordingly. + //! Mean electron energy [eV]. + double meanElectronEnergy() const { + return 3.0 / 2.0 * electronTemperature() * Boltzmann / ElectronCharge; + } + + //! Set mean electron energy [eV]. + //! This method also sets electron temperature accordingly. void setMeanElectronEnergy(double energy); - //! Get electron energy distribution type + //! Get electron energy distribution type. string electronEnergyDistributionType() const { return m_distributionType; } - //! Set electron energy distribution type + //! Set electron energy distribution type. void setElectronEnergyDistributionType(const string& type); //! Numerical quadrature method. Method: #m_quadratureMethod @@ -200,12 +210,7 @@ class PlasmaPhase: public IdealGasPhase m_quadratureMethod = method; } - //! Mean electron energy [eV] - double meanElectronEnergy() const { - return 3.0 / 2.0 * electronTemperature() * Boltzmann / ElectronCharge; - } - - //! Set flag of automatically normalize electron energy distribution + //! Set flag of automatically normalize electron energy distribution. //! Flag: #m_do_normalizeElectronEnergyDist void enableNormalizeElectronEnergyDist(bool enable) { m_do_normalizeElectronEnergyDist = enable; @@ -217,18 +222,12 @@ class PlasmaPhase: public IdealGasPhase return m_do_normalizeElectronEnergyDist; } - //! Electron Temperature (K) - //! @return The electron temperature of the phase - double electronTemperature() const override { - return m_electronTemp; - } - - //! Number of electron levels + //! Number of electron levels. size_t nElectronEnergyLevels() const { return m_nPoints; } - //! Number of electron collision cross sections + //! Number of electron collision cross sections. size_t nCollisions() const { return m_collisions.size(); } @@ -246,7 +245,6 @@ class PlasmaPhase: public IdealGasPhase return m_collisionRates[i]; } - //! Update the electron energy distribution. void updateElectronEnergyDistribution(); @@ -311,6 +309,16 @@ class PlasmaPhase: public IdealGasPhase m_electricField = EN * molarDensity() * Avogadro; // [V/m] } + /** + * The electron mobility (m²/V/s) + * @f[ + * \mu = \nu_d / E, + * @f] + * where @f$ \nu_d @f$ is the drift velocity (m²/s), and @f$ E @f$ is the electric + * field strength (V/m). + */ + double electronMobility() const; + /** * The elastic power loss [J/s/m³] * @f[ @@ -322,16 +330,6 @@ class PlasmaPhase: public IdealGasPhase */ double elasticPowerLoss(); - /** - * The electron mobility (m²/V/s) - * @f[ - * \mu = \nu_d / E, - * @f] - * where @f$ \nu_d @f$ is the drift velocity (m²/s), and @f$ E @f$ is the electric - * field strength (V/m). - */ - double electronMobility() const; - /** * The joule heating power (W/m³) * @f[ @@ -350,6 +348,7 @@ class PlasmaPhase: public IdealGasPhase void endEquilibrate() override; double intrinsicHeating() override; + //! @} // ================================================================= // // ================================================================= // @@ -358,151 +357,191 @@ class PlasmaPhase: public IdealGasPhase //! Return the Molar enthalpy. Units: J/kmol. /*! - * For an ideal gas mixture with additional electron, + * For an ideal gas mixture with electrons at a different temperature + * than the heavy species, the molar enthalpy is calculated as: * @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), + * \begin{align} + * \hat h(T, T_\text{e}, P) + * &= \sum_{k} X_k \tilde{h}_k(T_k) + * = \sum_{k} X_k h^\text{ref}_k(T_k) \\ + * &= \sum_{k \neq e} X_k h^\text{ref}_k(T) + * + X_\text{e} h_\text{e}^\text{ref}(T_\text{e}) + * \end{align} * @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 + * where heavy-species properties are evaluated at @f$ T @f$, electron properties + * at @f$ T_\text{e} @f$, and @f$ P @f$ is the total pressure of the mixture. + * For an ideal gas, enthalpy is indepedant of pressure. + * The standard-state pure-species enthalpies @f$ h^\text{ref}_k(T_k) @f$ are + * computed by the species thermodynamic property manager. */ double enthalpy_mole() const override; //! Return the molar entropy. Units: J/kmol/K. /*! - * For an ideal gas mixture with an additional electron species, - * @f[ - * \hat s(T,T_e,p,X) - * = \sum_{k\neq k_e} X_k\left[\hat s_k^0(T) - R\ln\left(\frac{X_k p}{p^0}\right)\right] - * + X_{k_e}\left[\hat s_{k_e}^0(T_e) - R\ln\left(\frac{X_{k_e} p_e}{p^0}\right)\right], - * @f] - * where heavy-species properties are evaluated at @f$T@f$, electron properties at - * @f$T_e@f$, and the electron mixing term uses the electron pressure - * @f$p_e = n_{k_e} R T_e@f$. - * - * @see MultiSpeciesThermo - */ + * For an ideal gas mixture with electrons at a different temperature + * than the heavy species, the molar entropy is calculated as: + * @f[ + * \begin{align} + * \hat s(T,T_\text{e},P) + * &= \sum_{k} X_k \tilde{s}_k(T_k, P) \\ + * &= \sum_{k} X_k \left[s^\text{ref}_k(T_k) + * - R \ln \frac{X_k P}{P^\text{ref}} \right] \\ + * \end{align} + * @f] + * where heavy-species properties are evaluated at @f$ T @f$, electron properties + * at @f$ T_\text{e} @f$, and @f$ P @f$ is the total pressure of the mixture. + */ double entropy_mole() const override; //! Return the molar Gibbs free energy. Units: J/kmol. /*! - * For an ideal gas mixture with an additional electron species, - * @f[ - * \hat g(T, T_e, p, X) = \sum_k X_k \mu_k(T_k, p_k, X), - * @f] - * where heavy species use the gas temperature @f$T@f$ and bulk pressure, - * while the electron chemical potential uses @f$T_e@f$ and - * @f$p_e = n_{k_e} R T_e@f$. - * - * @see MultiSpeciesThermo - */ + * For an ideal gas mixture with electrons at a different temperature + * than the heavy species, the molar Gibbs free energy is calculated as: + * @f[ + * \begin{align} + * \hat g(T,T_\text{e},P) + * &= \sum_{k} X_k \tilde{g}_k(T_k, P) + * = \sum_{k} X_k \mu_k(T_k, P) \\ + * &= \sum_{k} X_k \left[\mu^\text{ref}_k(T_k) + * + R T_k \ln \left( \frac{X_k P}{P^\text{ref}} \right) \right] \\ + * \end{align} + * @f] + * where heavy-species properties are evaluated at @f$ T @f$, electron properties + * at @f$ T_\text{e} @f$, and @f$ P @f$ is the total pressure of the mixture. + */ double gibbs_mole() const override; //! Return the molar internal energy. Units: J/kmol. /*! - * For an ideal gas mixture with an additional electron species, - * @f[ - * \hat u(T,T_e) = \sum_{k \neq k_e} X_k \hat u^0_k(T) + X_{k_e} \hat u^0_{k_e}(T_e), - * @f] - * where @f$\hat u^0_k(T) = \hat h^0_k(T) - R T@f$ for heavy species and - * @f$\hat u^0_{k_e}(T_e) = \hat h^0_{k_e}(T_e) - R T_e@f$ for electrons. - * - * @see MultiSpeciesThermo - */ + * For an ideal gas mixture with electrons at a different temperature + * than the heavy species, the molar internal energy is calculated as: + * @f[ + * \hat u(T,T_\text{e},P) + * = \sum_{k} X_k \tilde{u}_k(T_k) + * = \sum_{k} X_k (\tilde{h}_k(T_k) - R T_k) + * @f] + * where heavy-species properties are evaluated at @f$ T @f$, electron properties + * at @f$ T_\text{e} @f$, and @f$ P @f$ is the total pressure of the mixture. + * For an ideal gas, enthalpy is indepedant of pressure. + */ double intEnergy_mole() const override; - double cp_mole() const override; + // No need to override IdealGasPhase:cp_mole, as the computation + // does not imply temperature. + // double cp_mole() const override; // double cp_mass() const; // Already defined in ThermoPhase // double cv_mole() const; // Already defined in IdealGasPhase + //! @} // ================================================================= // // ================================================================= // //! @name Mechanical Equation of State //! @{ - //! Return the mean temperature of the plasma phase [K]. + //! Return the mean temperature of the plasma phase. Units: K. /*! * In a plasma phase, the electron temperature can differ from the * heavy-species (gas) temperature. Therefore, the mean temperature is * defined as a mole-fraction-weighted average of the electron and * heavy-species temperatures: * @f[ - * \bar{T} = \sum_{k \neq k_e} X_k T_g + X_{k_e} T_e - * = (1 - X_{k_e}) T_g + X_{k_e} T_e - * = T_g + X_{k_e} (T_e - T_g) + * \begin{align} + * \overline{T} &= \sum_{k \neq e} X_k T + X_\text{e} T_\text{e} \\ + * &= (1 - X_\text{e}) T + X_\text{e} T_\text{e} \\ + * &= T + X_\text{e} (T_\text{e} - T) + * \end{align} * @f] - * See the `pressure()` method for usage of the mean temperature in + * See the #pressure() method for usage of the mean temperature in * calculating the pressure of the plasma phase. */ double meanTemperature() const; - //! Return the pressure of the plasma phase [Pa]. + //! Return the pressure of the plasma phase. Units: Pa. /*! * The pressure of the plasma phase is calculated using the mean temperature, * which is a mole-fraction-weighted average of the electron and heavy-species * temperatures. * @f[ - * P = \sum_k n_k k_B T_k - * = \sum_{k \neq e} n_k k_B T_g + n_{e} k_B T_e - * = (n_{total} - n_{e}) k_B T_g + n_{e} k_B T_e - * = n_{total} (1 - X_e) k_B T_g + n_{total} X_e k_B T_e - * = n_{total} k_B \bar{T} + * \begin{align} + * P &= \sum_k n_k k_\text{B} T_k \\ + * &= \sum_{k \neq e} n_k k_\text{B} T + * + n_\text{e} k_\text{B} T_\text{e} \\ + * &= (n_\text{total} - n_\text{e}) k_\text{B} T + * + n_\text{e} k_\text{B} T_\text{e} \\ + * &= n_\text{total} (1 - X_\text{e}) k_\text{B} T + * + n_\text{total} X_\text{e} k_\text{B} T_\text{e} \\ + * &= n_\text{total} k_\text{B} \overline{T} \\ + * \end{align} * @f] - * where @f$ \bar{T} @f$ is the mean temperature of the plasma phase, - * defined in the `meanTemperature()` method. - * Here, @f$ n_k @f$ is the number density of species @f$ k @f$, - * @f$ n_{total} @f$ is the total number density of the phase, - * @f$ X_e @f$ is the mole fraction of electrons, @f$ T_g @f$ is the - * heavy-species (gas) temperature, @f$ T_e @f$ is the electron temperature, - * and @f$ k_B @f$ is the Boltzmann constant. + * where @f$ \overline{T} @f$ is the mean temperature of the plasma phase, + * defined in the #meanTemperature() method. + * Here, @f$ n_k @f$ is the number density of species *k* [in 1/m³], + * @f$ n_\text{total} @f$ is the total number density of the phase [in 1/m³], + * @f$ X_\text{e} @f$ is the mole fraction of electrons, + * @f$ T @f$ is the heavy-species (gas) temperature [K], + * @f$ T_\text{e} @f$ is the electron temperature [K], + * and @f$ k_\text{B} @f$ is the Boltzmann constant [J/K]. * * The number density times Boltzmann constant can be expressed as - * @f$ n_{total} k_B = C_{total} R @f$, where @f$ C_{total} @f$ is the - * total molar concentration of the phase and @f$ R @f$ is the gas constant. + * @f$ n_\text{total} k_\text{B} = C_\text{total} R @f$, where + * @f$ C_\text{total} @f$ is the total molar concentration of the phase [kmol/m³] + * and @f$ R @f$ is the gas constant [J/kmol/K], so that: + * @f[ + * P = C_\text{total} R \overline{T}. + * @f] + * + * Here, the pressure is set via the density, via: + * @f[ + * P = R \overline{T} \frac{\rho}{\overline{W}} + * @f] + * where @f$ \overline{W} @f$ is the mean molecular weight. */ double pressure() const override; - //! Set the pressure at constant temperature and composition. + //! Set the pressure at constant temperature and composition. Units: Pa. /*! - * Units: Pa. * This method is implemented by setting the mass density to * @f[ * \rho = \frac{P \overline{W}}{\hat R \overline{T} }. * @f] * - * @param p Pressure (Pa) + * where @f$ \overline{W} @f$ is the mean molecular weight, + * and @f$ \overline{T} @f$ is the mean temperature (see #meanTemperature()). + * + * @param p Pressure [Pa]. */ void setPressure(double p) override { setDensity(p * meanMolecularWeight() / (GasConstant * meanTemperature())); } - //! Return the Gas Constant multiplied by the current electron temperature - /*! - * The units are Joules kmol-1 - */ + //! Return the Gas Constant multiplied by the current electron temperature [J/kmol]. double RTe() const { return electronTemperature() * GasConstant; } - /** - * Electron pressure. Units: Pa. - * @f[P = n_{k_e} R T_e @f] + //! Return the electron pressure. Units: Pa. + /*! + * The partial pressure of electrons is calculated using the electron temperature: + * @f[ + * P_\text{e} = n_\text{e} k_\text{B} T_\text{e} = C_\text{e} R T_\text{e}, + * @f] + * where @f$ C_\text{e} @f$ is the molar concentration of electrons [in kmol/m³], + * @f$ R @f$ is the gas constant [J/kmol/K], + * and @f$ T_\text{e} @f$ is the electron temperature [K]. */ virtual double electronPressure() const { return GasConstant * concentration(m_electronSpeciesIndex) * electronTemperature(); } + //! Raise NotImplementedError, as there is an ambiguity on the temperature to use. double thermalExpansionCoeff() const override { - // Which temperature to use here? throw NotImplementedError("PlasmaPhase::thermalExpansionCoeff"); } + //! Raise NotImplementedError, as there is an ambiguity on the temperature to use. double soundSpeed() const override { - // Which temperature to use here? throw NotImplementedError("PlasmaPhase::soundSpeed"); } @@ -513,19 +552,26 @@ class PlasmaPhase: public IdealGasPhase //! @{ //! Returns the standard concentration @f$ C^0_k @f$, which is used to - //! normalize the generalized concentration. + //! normalize the generalized concentration. Units: m³/kmol. /*! * This is defined as the concentration by which the generalized * concentration is normalized to produce the activity. Since the activity * for an ideal gas mixture is simply the mole fraction, for an ideal gas * @f$ C^0_k = P/\hat R T @f$. - * For a multi-temperature plasma phase, the mean temperature is used instead: - * @f[$ C^0_k = P/\hat R \overline{T} @f$. + * For a multi-temperature system, this translates to: + * @f[ + * \begin{align} + * C^0_k &= \frac{C_k}{a_k} + * &= \frac{X_k \frac{P}{R \overline{T}}}{X_k} + * &= \frac{P}{R \overline{T}} + * \end{align} + * @f] + * where @f$C_k@f$ is the molar concentration of species @f$k@f$ [in kmol/m³], + * @f$ a_k @f$ is the activity of species *k*, which is equal to the + * mole fraction @f$ X_k @f$. * * @param k Parameter indicating the species. The default * is to assume this refers to species 0. - * @return - * Returns the standard Concentration in units of m3 kmol-1. */ double standardConcentration(size_t k=0) const override; @@ -535,15 +581,75 @@ class PlasmaPhase: public IdealGasPhase //! @name Partial Molar Properties of the Solution //! @{ + //! Return the chemical potentials of the species in the solution. Units: J/kmol. + /*! + * The chemical potential of species *k* is calculated as: + * @f[ + * \begin{align} + * \mu_k(T_k, X_k, P) + * &= \mu_k^o(T_k, P) + R T_k \ln(X_k) \\ + * &= \mu^\text{ref}_k(T_k)(T_k) + * + R T_k \ln\left(\frac{P}{P^\text{ref}}\right) + * + R T_k \ln(X_k) \\ + * &= h^\text{ref}_k(T_k) + * - T_k s^\text{ref}_k(T_k) + * + R T_k \ln\left(\frac{P X_k}{P^\text{ref}}\right) \\ + * \end{align} + * @f] + * Note that it is also equal to the partial molar Gibbs free energy. + */ void getChemPotentials(span mu) const override; + + //! Return the partial molar enthalpies of the species in the solution. Units: J/kmol. + /*! + * The partial molar enthalpy of species *k* is + * @f$ \tilde{h}_k(T_k,P) = h^o_k(T,P) = h^{ref}_k(T_k) @f$, + * where heavy-species properties are evaluated at the gas temperature @f$ T @f$, + * and electron properties are evaluated at the electron temperature @f$ T_e @f$. + */ void getPartialMolarEnthalpies(span hbar) const override; - // void getPartialMolarEntropies(span sbar) const override; + + //! Return the partial molar entropies of the species in the solution. Units: J/kmol/K. + /*! + * The partial molar enthalpy of species *k* is: + * @f[ + * \tilde{s}_k(T_k, P) = s^\text{ref}_k(T_k) - R \ln \frac{X_k P}{P^\text{ref}} + * @f] + * where heavy-species properties are evaluated at the gas temperature, + * and electron properties are evaluated at the electron temperature. + * Here, @f$ P @f$ is the total pressure of the mixture, as defined in #pressure(). + */ + void getPartialMolarEntropies(span sbar) const override; + + //! Return the partial molar internal energies of the species in the solution. Units: J/kmol. + /*! + * The partial molar internal energy of species *k* is calculated as: + * @f[ + * \tilde{u}_k(T_k,P) = \tilde{h}_k(T_k,P) - R T_k = h^{ref}_k(T_k) - R T_k + * @f] + * where @f$ \tilde{h}_k @f$ is the partial molar enthalpy of species *k*, and + * @f$ T_k @f$ is the temperature at which the partial molar enthalpy is evaluated. + * For heavy species, the partial molar enthalpy is evaluated at the gas temperature, + * while for electrons, it is evaluated at the electron temperature. + */ void getPartialMolarIntEnergies(span ubar) const override; + + //! Return the partial molar heat capacities of the species in the solution. Units: J/kmol/K. + /*! + * Since the computation of the partial molar enthalpy does not depend on the temperature, + * there is no need to override the method (since #updateThermo has been overriden). + */ // void getPartialMolarCp(span cpbar) const override; - // Whenever a temperature can be defined, the following relation holds: - // h_k = u_k + R * T_k = u_k + p * v_k - // Therefore, v_k = R * T_k / p + //! Return the partial molar volumes of the species in the solution. Units: m³/kmol. + /*! + * For a multitemperature system, + * @f[ + * v_k = \frac{R T_k}{P} + * @f] + * where @f$ T_k @f$ is the temperature at which the + * partial molar enthalpy is evaluated. + */ void getPartialMolarVolumes(span vbar) const override; //! @} @@ -552,26 +658,45 @@ class PlasmaPhase: public IdealGasPhase //! @name Properties of the Standard State of the Species in the Solution //! @{ - //! Return the standard chemical potentials of the species [J/kmol]. + //! Return the standard chemical potentials of the species. Units: J/kmol. /*! - * For heavy species, this is identical to the IdealGasPhase - * implementation. For electrons, the standard chemical potential - * is evaluated at the electron temperature: + * The standard chemical potentials (or standard state Gibbs free energy) + * of species *k* is calculated as: * @f[ - * \mu^0_{e}(T_e) = g_k^0(T_e) + RT_e \ln \left(\frac{P}{P^0}\right). + * \begin{align} + * \mu^0_k(T_k, X_k, P) + * &= mu^\text{ref}_k(T_k)(T_k) + * + R T_k \ln\left(\frac{P}{P^\text{ref}}\right) \\ + * &= h^\text{ref}_k(T_k) - T_k s^\text{ref}_k(T_k) + * + R T_k \ln\left(\frac{P X_k}{P^\text{ref}}\right) \\ + * \end{align} * @f] + * Here, @f$ P @f$ is the total pressure of the mixture, as defined in #pressure(). */ void getStandardChemPotentials(span muStar) const override; + // Below are 5 methods already defined in IdealGasPhase, that do not need + // to be overridden since `updateThermo` already updates the standard-state + // properties for electrons at the electron temperature. + // And for entropy (and thus Gibbs free energy), for all species, the total + // pressure is used in the logarithmic term (not the partial pressure). + // { // void getEnthalpy_RT(span hrt) const override; // void getEntropy_R(span sr) const override; // void getGibbs_RT(span grt) const override; // void getIntEnergy_RT(span urt) const override; // void getCp_R(span cpr) const override; + // } - // Whenever a temperature can be defined, the following relation holds: - // h_k = u_k + R * T_k = u_k + p * v_k - // Therefore, v_k = R * T_k / p + //! Return the standard molar volumes of the species. Units: m³/kmol. + /*! + * For a multitemperature system, + * @f[ + * v_k = \frac{R T_k}{P}, + * @f] + * where @f$ T_k @f$ is the temperature at which + * the standard molar volume is evaluated. + */ void getStandardVolumes(span vol) const override; //! @} @@ -580,23 +705,38 @@ class PlasmaPhase: public IdealGasPhase //! @name Thermodynamic Values for the Species Reference States //! @{ + // Below are 5 methods already defined in IdealGasPhase, that do not need + // to be overridden, since they do not depend on the temperature. + // { // void getEnthalpy_RT_ref(span hrt) const override; // void getGibbs_RT_ref(span grt) const override; + // void getEntropy_R_ref(span er) const override; + // void getIntEnergy_RT_ref(span urt) const override; + // void getCp_R_ref(span cprt) const override; + // } - //! Return the reference-state Gibbs free energys of the species [J/kmol]. + //! Return the reference chemical potentials of the species. Units: J/kmol. /*! - * For heavy species, this is identical to the IdealGasPhase - * implementation. For electrons, the reference-state Gibbs free energy - * is evaluated at the electron temperature: + * The reference chemical potentials (or reference state Gibbs free energy) + * of species *k* is calculated as: * @f[ - * \hat{g}^0_{e}(T_e) = \hat{h}^0_{e}(T_e) - T_e \hat{s}^0_{e}(T_e). + * \mu^\text{ref}_k(T_k) = h^\text{ref}_k(T_k) - T_k s^\text{ref}_k(T_k) * @f] */ void getGibbs_ref(span g) const override; - // void getEntropy_R_ref(span er) const override; - // void getIntEnergy_RT_ref(span urt) const override; - // void getCp_R_ref(span cprt) const override; + //! Return the molar volumes of the species reference states. Units: m³/kmol. + /*! + * The molar volumes of the species reference states of species *k* + * is calculated as: + * @f[ + * v^o_k(T_k) = \frac{R T_k}{P^\text{ref}} + * @f] + * + * @param vol Output vector containing the standard state volumes. + * Length: m_kk. + + */ void getStandardVolumes_ref(span vol) const override; //! @} @@ -721,31 +861,31 @@ class PlasmaPhase: public IdealGasPhase //! Length: #m_nPoints Eigen::ArrayXd m_electronEnergyDist; - //! Index of electron species + //! Index of electron species. size_t m_electronSpeciesIndex = npos; - //! Electron temperature [K] + //! Electron temperature [K]. double m_electronTemp; - //! Electron energy distribution type + //! Electron energy distribution type. Can be "isotropic", "discretized" or "Boltzmann-two-term". string m_distributionType = "isotropic"; - //! Numerical quadrature method for electron energy distribution + //! Numerical quadrature method for electron energy distribution. string m_quadratureMethod = "simpson"; - //! Flag of normalizing electron energy distribution + //! Flag of normalizing electron energy distribution. bool m_do_normalizeElectronEnergyDist = true; - //! Indices of inelastic collisions in m_crossSections + //! Indices of inelastic collisions in m_crossSections. vector m_kInelastic; - //! Indices of elastic collisions in m_crossSections + //! Indices of elastic collisions in m_crossSections. vector m_kElastic; - //! electric field [V/m] + //! electric field [V/m]. double m_electricField = 0.0; - //! electric field freq [Hz] + //! electric field freq [Hz]. double m_electricFieldFrequency = 0.0; //! Cross section data. m_crossSections[i][j], where i is the specific process, @@ -766,7 +906,7 @@ class PlasmaPhase: public IdealGasPhase /*! The elastic electron energy loss coefficient for species k is, * @f[ * K_k = \frac{2 m_e}{m_k} \sqrt{\frac{2 e}{m_e}} \int_0^{\infty} \sigma_k - * \epsilon^2 \left( F_0 + \frac{k_B T}{e} + * \epsilon^2 \left( F_0 + \frac{k_\text{B} T}{e} * \frac{\partial F_0}{\partial \epsilon} \right) d \epsilon, * @f] * where @f$ m_e @f$ [kg] is the electron mass, @f$ \epsilon @f$ [V] is the diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 0fc61543599..a2b6f371586 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -25,13 +25,13 @@ PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_) { initThermoFile(inputFile, id_); - // initial electron temperature + // Initial electron temperature. m_electronTemp = temperature(); - // Initialize the Boltzmann Solver + // Initialize the Boltzmann Solver. m_eedfSolver = make_unique(this); - // Set Energy Grid (Hardcoded Defaults for Now) + // Set Energy Grid (Hardcoded Defaults for Now). double kTe_max = 60; size_t nGridCells = 301; m_nPoints = nGridCells + 1; @@ -47,7 +47,7 @@ PlasmaPhase::~PlasmaPhase() soln->kinetics()->removeReactionAddedCallback(this); } for (size_t k = 0; k < nCollisions(); k++) { - // remove callback + // Remove callback. m_collisions[k]->removeSetRateCallback(this); } } @@ -405,12 +405,12 @@ void PlasmaPhase::setDiscretizedElectronEnergyDist(span levels, void PlasmaPhase::updateElectronTemperatureFromEnergyDist() { - // calculate mean electron energy and electron temperature + // Calculate mean electron energy and electron temperature. Eigen::ArrayXd eps52 = m_electronEnergyLevels.pow(5./2.); double epsilon_m = 2.0 / 5.0 * numericalQuadrature(m_quadratureMethod, m_electronEnergyDist, eps52); if (epsilon_m < 0.0 && m_quadratureMethod == "simpson") { - // try trapezoidal method + // Try trapezoidal method. epsilon_m = 2.0 / 5.0 * numericalQuadrature( "trapezoidal", m_electronEnergyDist, eps52); } @@ -436,8 +436,8 @@ void PlasmaPhase::setCollisions() return; } - // add collision from the initial list of reactions. Only add reactions we - // haven't seen before + // Add collision from the initial list of reactions. + // Only add reactions we haven't seen before. set existing; for (auto& R : m_collisions) { existing.insert(R.get()); @@ -451,8 +451,8 @@ void PlasmaPhase::setCollisions() addCollision(R); } - // register callback when reaction is added later - // Modifying collision reactions is not supported + // 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") { @@ -466,18 +466,18 @@ void PlasmaPhase::addCollision(shared_ptr collision) { size_t i = nCollisions(); - // setup callback to signal updating the cross-section-related - // parameters + // Setup callback to signal updating the cross-section-related + // parameters. collision->registerSetRateCallback(this, [this, i, collision]() { m_interp_cs_ready[i] = false; m_collisionRates[i] = std::dynamic_pointer_cast(collision->rate()); }); - // Identify target species for electron-collision reactions + // Identify target species for electron-collision reactions. string target; for (const auto& [name, _] : collision->reactants) { - // Reactants are expected to be electrons and the target species + // Reactants are expected to be electrons and the target species. if (name != electronSpeciesName()) { m_targetSpeciesIndices.emplace_back(speciesIndex(name, true)); target = name; @@ -494,11 +494,11 @@ void PlasmaPhase::addCollision(shared_ptr collision) std::dynamic_pointer_cast(collision->rate())); m_interp_cs_ready.emplace_back(false); - // resize parameters + // Resize parameters. m_elasticElectronEnergyLossCoefficients.resize(nCollisions()); updateInterpolatedCrossSection(i); - // Set up data used by Boltzmann solver + // Set up data used by Boltzmann solver. auto& rate = *m_collisionRates.back(); string kind = m_collisionRates.back()->kind(); @@ -539,12 +539,12 @@ bool PlasmaPhase::updateInterpolatedCrossSection(size_t i) void PlasmaPhase::updateElectronEnergyDistDifference() { m_electronEnergyDistDiff.resize(nElectronEnergyLevels()); - // Forward difference for the first point + // Forward difference for the first point. m_electronEnergyDistDiff[0] = (m_electronEnergyDist[1] - m_electronEnergyDist[0]) / (m_electronEnergyLevels[1] - m_electronEnergyLevels[0]); - // Central difference for the middle points + // Central difference for the middle points. for (size_t i = 1; i < m_nPoints - 1; i++) { double h1 = m_electronEnergyLevels[i+1] - m_electronEnergyLevels[i]; double h0 = m_electronEnergyLevels[i] - m_electronEnergyLevels[i-1]; @@ -554,7 +554,7 @@ void PlasmaPhase::updateElectronEnergyDistDifference() (h1 * h0) / (h1 + h0); } - // Backward difference for the last point + // Backward difference for the last point. m_electronEnergyDistDiff[m_nPoints-1] = (m_electronEnergyDist[m_nPoints-1] - m_electronEnergyDist[m_nPoints-2]) / @@ -564,11 +564,11 @@ void PlasmaPhase::updateElectronEnergyDistDifference() void PlasmaPhase::updateElasticElectronEnergyLossCoefficients() { - // cache of cross section plus distribution plus energy-level number + // Cache of cross section plus distribution plus energy-level number. static const int cacheId = m_cache.getId(); CachedScalar last_stateNum = m_cache.getScalar(cacheId); - // combine the distribution and energy level number + // Combine the distribution and energy level number. int stateNum = m_distNum + m_levelNum; vector interpChanged(m_collisions.size()); @@ -577,15 +577,15 @@ void PlasmaPhase::updateElasticElectronEnergyLossCoefficients() } if (last_stateNum.validate(temperature(), stateNum)) { - // check each cross section, and only update coefficients that - // the interpolated cross sections change + // Check each cross section, and only update coefficients that + // the interpolated cross sections change. for (size_t i = 0; i < m_collisions.size(); i++) { if (interpChanged[i]) { updateElasticElectronEnergyLossCoefficient(i); } } } else { - // update every coefficient if distribution, temperature, + // Update every coefficient if distribution, temperature, // or energy levels change. for (size_t i = 0; i < m_collisions.size(); i++) { updateElasticElectronEnergyLossCoefficient(i); @@ -604,10 +604,10 @@ void PlasmaPhase::updateElasticElectronEnergyLossCoefficient(size_t i) m_collisionRates[i]->crossSectionInterpolated().size() ); - // Mass ratio calculation + // Mass ratio calculation. double mass_ratio = ElectronMass / molecularWeight(k) * Avogadro; - // Calculate the rate using Simpson's rule or trapezoidal rule + // Calculate the rate using Simpson's rule or trapezoidal rule. Eigen::ArrayXd f0_plus = m_electronEnergyDist + Boltzmann * temperature() / ElectronCharge * m_electronEnergyDistDiff; m_elasticElectronEnergyLossCoefficients[i] = 2.0 * mass_ratio * gamma * @@ -777,10 +777,25 @@ void PlasmaPhase::getChemPotentials(span mu) const void PlasmaPhase::getPartialMolarEnthalpies(span hbar) const { + // Since the `updateThermo` is overriden in `PlasmaPhase`, + // `enthalpy_RT_ref` returns \tilde{h}_k(T_k) / (R * T_k). + // When calling `IdealGasPhase::getPartialMolarEnthalpies(hbar)`, + // the `hbar` array is equal to \tilde{h}_k(T_k) * (R * T) / (R * T_k). + // For all heavy species, T_k == T, so we get \tilde{h}_k(T). + // For electrons, we need to multiply by T_e/T to get \tilde{h}_k(T_e). IdealGasPhase::getPartialMolarEnthalpies(hbar); hbar[m_electronSpeciesIndex] *= electronTemperature() / temperature(); } +void PlasmaPhase::getPartialMolarEntropies(span sbar) const +{ + // Since the `updateThermo` is overriden in `PlasmaPhase`, + // `entropy_R_ref` returns s^\text{ref}_k(T_k)/R. + // When calling `IdealGasPhase::getPartialMolarEntropies(hbar)`, + // the `sbar` array is equal to s^\text{ref}_k(T_k)*R/R - R ln(X_k P/P^ref). + // Therefore, there is no need to correct for temperature. + IdealGasPhase::getPartialMolarEntropies(sbar); +} void PlasmaPhase::getPartialMolarIntEnergies(span ubar) const { @@ -809,17 +824,17 @@ void PlasmaPhase::getPartialMolarVolumes(span vbar) const void PlasmaPhase::getStandardChemPotentials(span muStar) const { - // After calling PlasmaPhase::getGibbs_ref, muStar = g_k(T_k). - // g_k is evaluated at T for heavy species and at Te for electrons. + // After calling PlasmaPhase::getGibbs_ref, muStar = mu^\text{ref}_k(T_k)(T_k). + // mu^\text{ref} is evaluated at T for heavy species and at Te for electrons. getGibbs_ref(muStar); - // Then, we need to add R*T_k*ln(P/Pref) to g_k. - // .. For heavy species, mu_star = g_k(T) + R*T*ln(P/Pref) + // Then, we need to add R*T_k*ln(P/Pref) to mu^\text{ref}. + // .. For heavy species, mu_star = mu^\text{ref}(T) + R*T*ln(P/Pref) double tmp = log(pressure() / refPressure()) * RT(); for (size_t k = 0; k < m_kk; k++) { muStar[k] += tmp; } - // .. For electrons, mu_star = g_k(Te) + R*T_e*ln(P/Pref) + // .. For electrons, mu_star = mu^\text{ref}(Te) + R*T_e*ln(P/Pref) size_t k = m_electronSpeciesIndex; muStar[k] -= log(pressure() / refPressure()) * RT(); muStar[k] += log(pressure() / refPressure()) * RTe(); @@ -834,12 +849,19 @@ void PlasmaPhase::getStandardVolumes(span vol) const vol[m_electronSpeciesIndex] = RTe() / pressure(); } + // ================================================================= // // Thermodynamic Values for the Species Reference States // // ================================================================= // void PlasmaPhase::getGibbs_ref(span g) const { + // Since the `updateThermo` is overriden in `PlasmaPhase`, + // `gibbs_RT_ref` returns \mu^\text{ref}_k(T_k) / (R * T_k). + // When calling `IdealGasPhase::getGibbs_ref(g)`, + // the `g` array is equal to \mu^\text{ref}_k(T_k) * (R * T) / (R * T_k). + // For all heavy species, T_k == T, so we get \mu^\text{ref}_k(T). + // For electrons, we need to multiply by T_e/T to get \mu^\text{ref}_k(T_e). IdealGasPhase::getGibbs_ref(g); g[m_electronSpeciesIndex] *= electronTemperature() / temperature(); } From c62696335278ede4857faade97cd765e57d96a66 Mon Sep 17 00:00:00 2001 From: pag1pag <66677348+pag1pag@users.noreply.github.com> Date: Tue, 26 May 2026 11:42:59 +0200 Subject: [PATCH 11/11] [Plasma] Update tests --- test/data/consistency-cases.yaml | 30 +++--- test/data/thermo_plasma_simple.yaml | 49 +++++++++ test/python/test_thermo.py | 78 ++++++++++++++ test/thermo/PlasmaPhase_Test.cpp | 136 ++++++++++++++++++++++++ test/thermo/thermoToYaml.cpp | 4 +- test/thermo_consistency/consistency.cpp | 86 +++++++++++---- 6 files changed, 349 insertions(+), 34 deletions(-) create mode 100644 test/data/thermo_plasma_simple.yaml create mode 100644 test/thermo/PlasmaPhase_Test.cpp diff --git a/test/data/consistency-cases.yaml b/test/data/consistency-cases.yaml index 3fcc34794f9..ea2fd9e4fc1 100644 --- a/test/data/consistency-cases.yaml +++ b/test/data/consistency-cases.yaml @@ -125,23 +125,29 @@ nitrogen-purefluid: plasma: setup: file: oxygen-plasma.yaml - phase: discretized-electron-energy-plasma + phase: isotropic-electron-energy-plasma known-failures: - g_eq_h_minus_Ts/.+: Test does not account for distinct electron temperature - u_eq_sum_uk_Xk/.+: Test does not account for distinct electron temperature - s_eq_sum_sk_Xk/.+: Test does not account for distinct electron temperature - cp_eq_dhdT/(0|2): Test does not account for distinct electron temperature - cv_eq_dudT/(0|2): Test does not account for distinct electron temperature - c_eq_sqrt_dP_drho_const_s/1: Test does not account for distinct electron temperature - cp_eq_.+/1: Test does not account for distinct electron temperature - cv_eq_.+/1: Test does not account for distinct electron temperature - c._eq_dsdT_const_._times_T: Test does not account for distinct electron temperature - cpk0_eq_dhk0dT: Test does not account for distinct electron temperature - cpRef_eq_dhRefdT: Test does not account for distinct electron temperature + c_eq_sqrt_dP_drho_const_s/.+: For now, `setState_SV` is not supported in PlasmaPhase. + cp_eq_dhdT/[345]: The differential methods are not yet implemented for two-temperature systems. + cv_eq_dudT/[345]: The differential methods are not yet implemented for two-temperature systems. + cp_eq_dsdT_const_p_times_T/[345]: + The differential methods are not yet implemented for two-temperature systems. + cv_eq_dsdT_const_v_times_T/[345]: + The differential methods are not yet implemented for two-temperature systems. + dSdv_const_T_eq_dPdT_const_V/[345]: + The differential methods are not yet implemented for two-temperature systems. + betaT_eq_minus_dmv_dP_const_T_div_mv/[345]: + The differential methods are not yet implemented for two-temperature systems. states: + # In the test cases below, it is implicitly assumed that the electron temperature + # equals the gas temperature. - {T: 300, P: 1 atm, X: {O2: 1.0, O2-: 1e-5, E: 1e-5}} - {T: 300, P: 1 atm, X: {E: 1.0}} - {T: 3500, P: 10 atm, X: {O2: 1.0, O2-: 2e-5, E: 2e-5}} + # In the test case below, set Te != T + - {T: 300, P: 1 atm, X: {O2: 1.0, O2-: 1e-5, E: 1e-5}, Te: 10000} + - {T: 300, P: 1 atm, X: {E: 1.0}, Te: 10000} + - {T: 3500, P: 10 atm, X: {O2: 1.0, O2-: 2e-5, E: 2e-5}, Te: 10000} debye-huckel-dilute: setup: diff --git a/test/data/thermo_plasma_simple.yaml b/test/data/thermo_plasma_simple.yaml new file mode 100644 index 00000000000..a77f046d15c --- /dev/null +++ b/test/data/thermo_plasma_simple.yaml @@ -0,0 +1,49 @@ +description: |- + This is a simple test of the plasma thermodynamic model, with a single plasma phase and a single gas phase. + The plasma phase has an isotropic electron energy distribution, and the gas phase is an ideal gas. + +phases: + - name: plasma + thermo: plasma + elements: [O, E] + species: [O, O+, e-] + kinetics: gas + state: { T: 300.0, P: 1 atm } + electron-energy-distribution: + type: isotropic + shape-factor: 1.0 + energy-levels: [0.0, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0] + mean-electron-energy: 1.0 eV + reactions: none + +species: + - name: O + composition: { O: 1 } + thermo: + model: NASA9 + temperature-ranges: [200.0, 1000.0, 6000.0, 2.0e+04] + data: + - [0.0, 0.0, 3, 0.0, 0.0, 0.0, 0.0, 3e+04, 8.5] + - [0.0, 0.0, 3, 0.0, 0.0, 0.0, 0.0, 3e+04, 8.5] + - [0.0, 0.0, 3, 0.0, 0.0, 0.0, 0.0, 3e+04, 8.5] + note: False, but for testing purposes. + - name: O+ + composition: { O: 1, E: -1 } + thermo: + model: NASA9 + temperature-ranges: [298.15, 1000.0, 6000.0, 2.0e+04] + data: + - [0.0, 0.0, 8, 0.0, 0.0, 0.0, 0.0, 2e+05, 10.0] + - [0.0, 0.0, 8, 0.0, 0.0, 0.0, 0.0, 2e+05, 10.0] + - [0.0, 0.0, 8, 0.0, 0.0, 0.0, 0.0, 2e+05, 10.0] + note: False, but for testing purposes. + - name: e- + composition: { E: 1 } + thermo: + model: NASA9 + temperature-ranges: [298.15, 1000.0, 6000.0, 2.0e+04] + data: + - [0.0, 0.0, 2.5, 0.0, 0.0, 0.0, 0.0, -750, -11.7] + - [0.0, 0.0, 2.5, 0.0, 0.0, 0.0, 0.0, -750, -11.7] + - [0.0, 0.0, 2.5, 0.0, 0.0, 0.0, 0.0, -750, -11.7] + note: False, but for testing purposes. diff --git a/test/python/test_thermo.py b/test/python/test_thermo.py index 8e71ce310ef..b608dc81f94 100644 --- a/test/python/test_thermo.py +++ b/test/python/test_thermo.py @@ -1439,15 +1439,18 @@ def test_add_multiple_electron_species(self, phase): def test_elastic_power_loss_low_T(self, phase): phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + phase.mean_electron_energy = 1.0 assert phase.elastic_power_loss == approx(6846332332) def test_elastic_power_loss_high_T(self, phase): # when T is as high as Te the energy loss rate becomes small phase.TPX = 4000, ct.one_atm, "O2:1, E:1e-5" + phase.mean_electron_energy = 1.0 assert phase.elastic_power_loss == approx(2865540) def test_elastic_power_loss_replace_rate(self, phase): phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + phase.mean_electron_energy = 1.0 rate = ct.ReactionRate.from_dict(self.collision_data) phase.reaction(1).rate = rate assert phase.elastic_power_loss == approx(11765800095) @@ -1456,11 +1459,13 @@ def test_elastic_power_loss_add_reaction(self, phase): phase2 = ct.Solution(thermo="plasma", kinetics="bulk", species=phase.species(), reactions=[]) phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + phase.mean_electron_energy = 1.0 phase.add_reaction(ct.Reaction.from_dict(self.collision_data, phase)) assert phase.elastic_power_loss == approx(18612132428) def test_elastic_power_loss_change_levels(self, phase): phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + phase.mean_electron_energy = 1.0 phase.electron_energy_levels = np.linspace(0,10,101) assert phase.elastic_power_loss == approx(113058853) @@ -1483,6 +1488,7 @@ def test_elastic_power_loss_change_mean_electron_energy(self, phase): def test_elastic_power_loss_change_shape_factor(self, phase): phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + phase.mean_electron_energy = 1.0 phase.isotropic_shape_factor = 1.1 assert phase.elastic_power_loss == approx(7408711810) @@ -1521,6 +1527,78 @@ def test_eedf_solver(self): l2_norm = np.linalg.norm(interp - reference_eedf) assert l2_norm < 1e-3 + def _make_simple_plasma_phase(self, Te, set_te_first=False): + phase = ct.Solution("thermo_plasma_simple.yaml", "plasma", transport_model=None) + phase.isotropic_shape_factor = 1.0 + + T = 300.0 # K + P = ct.one_atm + composition = "O: 0.8, O+: 0.1, e-: 0.1" + + if set_te_first: + phase.Te = Te + phase.TPX = T, P, composition + if not set_te_first: + phase.Te = Te + + return phase + + @staticmethod + def _expected_cp_mole(X): + cp_o = 3.0 * ct.gas_constant + cp_o_plus = 8.0 * ct.gas_constant + cp_e = 2.5 * ct.gas_constant + return X[0] * cp_o + X[1] * cp_o_plus + X[2] * cp_e + + @staticmethod + def _expected_h_mole(X, T, Te): + h_o = 3.0 * ct.gas_constant * T + 3e4 * ct.gas_constant + h_o_plus = 8.0 * ct.gas_constant * T + 2e5 * ct.gas_constant + h_e = 2.5 * ct.gas_constant * Te - 750.0 * ct.gas_constant + return X[0] * h_o + X[1] * h_o_plus + X[2] * h_e + + @pytest.mark.parametrize( + "Te,set_te_first", + [ + (300.0, True), + (300.0, False), + (30_000.0, True), + (30_000.0, False), + ], + ) + def test_rho_two_temperature_modes(self, Te, set_te_first): + phase = self._make_simple_plasma_phase(Te, set_te_first=set_te_first) + T_ref = phase.mean_temperature + expected_rho = phase.P * phase.mean_molecular_weight / (ct.gas_constant * T_ref) + assert phase.density == approx(expected_rho) + + @pytest.mark.parametrize( + "Te,set_te_first", + [ + (300.0, True), + (300.0, False), + (30_000.0, True), + (30_000.0, False), + ], + ) + def test_cp_mole_two_temperature_modes(self, Te, set_te_first): + phase = self._make_simple_plasma_phase(Te, set_te_first=set_te_first) + assert phase.cp_mole == approx(self._expected_cp_mole(phase.X)) + + @pytest.mark.parametrize( + "Te,set_te_first", + [ + (300.0, True), + (300.0, False), + (30_000.0, True), + (30_000.0, False), + ], + ) + def test_enthalpy_mole_two_temperature_modes(self, Te, set_te_first): + phase = self._make_simple_plasma_phase(Te, set_te_first=set_te_first) + expected_h = self._expected_h_mole(phase.X, phase.T, phase.Te) + assert phase.enthalpy_mole == approx(expected_h) + class TestImport: """ diff --git a/test/thermo/PlasmaPhase_Test.cpp b/test/thermo/PlasmaPhase_Test.cpp new file mode 100644 index 00000000000..a5a4ff2db30 --- /dev/null +++ b/test/thermo/PlasmaPhase_Test.cpp @@ -0,0 +1,136 @@ +#include "gtest/gtest.h" +#include "cantera/thermo/PlasmaPhase.h" +#include "cantera/thermo/ThermoFactory.h" +#include "cantera/thermo/Species.h" + +#include + + +namespace Cantera +{ + +class PlasmaPhase_Test : public testing::Test +{ +public: + PlasmaPhase_Test() { + auto thermo = newThermo("example_data/oxygen-plasma-itikawa.yaml"); + test_phase = std::dynamic_pointer_cast(thermo); + } + + shared_ptr test_phase; +}; + +TEST_F(PlasmaPhase_Test, setState_TP) +{ + test_phase->setState_TP(1000, 1e5); + EXPECT_DOUBLE_EQ(test_phase->temperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->pressure(), 1e5); + + test_phase->setState_TP(2000, 1e5); + EXPECT_DOUBLE_EQ(test_phase->temperature(), 2000); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 2000); + EXPECT_DOUBLE_EQ(test_phase->pressure(), 1e5); + + test_phase->setState_TP(2000, 2e5); + EXPECT_DOUBLE_EQ(test_phase->temperature(), 2000); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 2000); + EXPECT_DOUBLE_EQ(test_phase->pressure(), 2e5); +} + +TEST_F(PlasmaPhase_Test, setState_TgTeP) +{ + test_phase->setState_TgTeP(1000, 1000, 1e5); + EXPECT_DOUBLE_EQ(test_phase->temperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->pressure(), 1e5); + + test_phase->setState_TgTeP(2000, 1000, 1e5); + EXPECT_DOUBLE_EQ(test_phase->temperature(), 2000); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->pressure(), 1e5); + + test_phase->setState_TgTeP(1000, 2000, 1e5); + EXPECT_DOUBLE_EQ(test_phase->temperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 2000); + EXPECT_DOUBLE_EQ(test_phase->pressure(), 1e5); + + test_phase->setState_TgTeP(1000, 1000, 2e5); + EXPECT_DOUBLE_EQ(test_phase->temperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->pressure(), 2e5); +} + +TEST_F(PlasmaPhase_Test, setState_TD) +{ + test_phase->setState_TD(1000, 1e-2); + EXPECT_DOUBLE_EQ(test_phase->temperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->density(), 1e-2); + + test_phase->setState_TD(2000, 1e-2); + EXPECT_DOUBLE_EQ(test_phase->temperature(), 2000); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 2000); + EXPECT_DOUBLE_EQ(test_phase->density(), 1e-2); + + test_phase->setState_TD(2000, 2e-2); + EXPECT_DOUBLE_EQ(test_phase->temperature(), 2000); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 2000); + EXPECT_DOUBLE_EQ(test_phase->density(), 2e-2); +} + +TEST_F(PlasmaPhase_Test, setState_TgTeD) +{ + test_phase->setState_TgTeD(1000, 1000, 1e-2); + EXPECT_DOUBLE_EQ(test_phase->temperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->density(), 1e-2); + + test_phase->setState_TgTeD(2000, 1000, 1e-2); + EXPECT_DOUBLE_EQ(test_phase->temperature(), 2000); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->density(), 1e-2); + + test_phase->setState_TgTeD(1000, 2000, 1e-2); + EXPECT_DOUBLE_EQ(test_phase->temperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 2000); + EXPECT_DOUBLE_EQ(test_phase->density(), 1e-2); + + test_phase->setState_TgTeD(1000, 1000, 2e-2); + EXPECT_DOUBLE_EQ(test_phase->temperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 1000); + EXPECT_DOUBLE_EQ(test_phase->density(), 2e-2); +} + +TEST_F(PlasmaPhase_Test, setState) +{ + AnyMap state; + state["gas-temperature"] = 1234; + state["P"] = "5 bar"; + state["X"] = "O2: 0.8, O2+: 0.1, E: 0.1"; + test_phase->setState(state); + + EXPECT_DOUBLE_EQ(test_phase->temperature(), 1234); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 1234); + EXPECT_DOUBLE_EQ(test_phase->pressure(), 5e5); + EXPECT_DOUBLE_EQ(test_phase->moleFraction("O2"), 0.8); + EXPECT_DOUBLE_EQ(test_phase->moleFraction("O2+"), 0.1); + EXPECT_DOUBLE_EQ(test_phase->moleFraction("E"), 0.1); + + + AnyMap state2; + state2["T"] = 1234; + state2["Te"] = 4321; + state2["D"] = 1e-2; + state2["Y"] = "O2: 0.8, O2+: 0.1, E: 0.1"; + test_phase->setState(state2); + + EXPECT_DOUBLE_EQ(test_phase->temperature(), 1234); + EXPECT_DOUBLE_EQ(test_phase->electronTemperature(), 4321); + EXPECT_DOUBLE_EQ(test_phase->density(), 1e-2); + EXPECT_DOUBLE_EQ(test_phase->massFraction("O2"), 0.8); + EXPECT_DOUBLE_EQ(test_phase->massFraction("O2+"), 0.1); + EXPECT_DOUBLE_EQ(test_phase->massFraction("E"), 0.1); +} + +}; diff --git a/test/thermo/thermoToYaml.cpp b/test/thermo/thermoToYaml.cpp index bac65a4a1f5..d31346885b2 100644 --- a/test/thermo/thermoToYaml.cpp +++ b/test/thermo/thermoToYaml.cpp @@ -572,8 +572,8 @@ TEST_F(ThermoYamlRoundTrip, CoverageDependentSurface) TEST_F(ThermoYamlRoundTrip, IsotropicElectronEnergyPlasma) { roundtrip("oxygen-plasma.yaml", "isotropic-electron-energy-plasma"); - skip_cp = true; // Not implemented for PlasmaPhase - skip_entropy = true; // Not implemented for PlasmaPhase + skip_cp = false; // Not implemented for PlasmaPhase + skip_entropy = false; // Not implemented for PlasmaPhase compareThermo(800, 2*OneAtm); auto origPlasma = std::dynamic_pointer_cast(original); auto duplPlasma = std::dynamic_pointer_cast(duplicate); diff --git a/test/thermo_consistency/consistency.cpp b/test/thermo_consistency/consistency.cpp index 7ac47d1115f..bca5adab3f0 100644 --- a/test/thermo_consistency/consistency.cpp +++ b/test/thermo_consistency/consistency.cpp @@ -147,17 +147,7 @@ TEST_P(TestConsistency, h_eq_u_plus_Pv) { } catch (NotImplementedError& err) { GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; } - if (phase->type() == "plasma" && ke != npos) { - // Two-temperature identity for PlasmaPhase: - // h = u + p*v + X_e * R * (Te - T) - vector X(nsp); - phase->getMoleFractions(X); - double Xe = X[ke]; - EXPECT_NEAR(h, u + p * v + Xe * (RTe - RT), atol); - } else { - // Standard single-temperature identity - EXPECT_NEAR(h, u + p * v, atol); - } + EXPECT_NEAR(h, u + p * v, atol); } TEST_P(TestConsistency, g_eq_h_minus_Ts) { @@ -169,7 +159,32 @@ TEST_P(TestConsistency, g_eq_h_minus_Ts) { } catch (NotImplementedError& err) { GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; } - EXPECT_NEAR(g, h - T * s, atol); + if (phase->type() == "plasma") { + // Whenever a temperature can be defined, the following relation holds: + // g_k = h_k - T_k * s_k + // When multiplying by mole fraction and summing over all species, we get: + // sum_k X_k * g_k = sum_k X_k * h_k - sum_k X_k * T_k * s_k + // The left side is simply g. The first term on the right side is h. + // The second term on the right side can be separated into electron and + // non-electron contributions (for a 2 temperature plasma model): + // sum_k X_k * T_k * s_k = T * sum_{k!=e} X_k * s_k + Te * X_e * s_e + // = T * (s - X_e * s_e) + Te * X_e * s_e + // = T * s + X_e * s_e * (Te - T) + // Rearranging gives: + // g = h - T * s - X_e * s_e * (Te - T) + + // Get partial molar entropy of electron species. + vector sk(nsp); + phase->getPartialMolarEntropies(sk); + double se = sk[ke]; + // Get mole fraction of electron species. + double Xe = phase->moleFraction(ke); + + EXPECT_NEAR(g, h - T * s - Xe * se * (Te - T), atol); + } + else { + EXPECT_NEAR(g, h - T * s, atol); + } } TEST_P(TestConsistency, hk_eq_uk_plus_P_vk) @@ -183,9 +198,7 @@ TEST_P(TestConsistency, hk_eq_uk_plus_P_vk) GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; } for (size_t k = 0; k < nsp; k++) { - if (k != ke) { - EXPECT_NEAR(hk[k], uk[k] + p * vk[k], atol) << "k = " << k; - } // not applicable for electron + EXPECT_NEAR(hk[k], uk[k] + p * vk[k], atol) << "k = " << k; } } @@ -256,7 +269,10 @@ TEST_P(TestConsistency, gk_eq_hk_minus_T_sk) for (size_t k = 0; k < nsp; k++) { if (k != ke) { EXPECT_NEAR(gk[k], hk[k] - T * sk[k], atol) << "k = " << k; - } // not applicable for electron + } + else { + EXPECT_NEAR(gk[k], hk[k] - Te * sk[k], atol) << "k = " << k; + } } } @@ -580,7 +596,11 @@ TEST_P(TestConsistency, hk0_eq_uk0_plus_p_vk0) GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; } for (size_t k = 0; k < nsp; k++) { - EXPECT_NEAR(h0[k] * RT, u0[k] * RT + p * v0[k], atol) << "k = " << k; + if (k != ke) { + EXPECT_NEAR(h0[k] * RT, u0[k] * RT + p * v0[k], atol) << "k = " << k; + } else { + EXPECT_NEAR(h0[k] * RTe, u0[k] * RTe + p * v0[k], atol) << "k = " << k; + } } } @@ -609,14 +629,27 @@ TEST_P(TestConsistency, cpk0_eq_dhk0dT) } catch (NotImplementedError& err) { GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; } + // Perturb temperature. double T1 = phase->temperature(); double dT = 1e-5 * phase->temperature(); + // For plasma phases, also perturb the electron temperature. + double Te = phase->electronTemperature(); + double dTe = 1e-5 * Te; + phase->setState_TP(T1 + dT, phase->pressure()); + if (phase->type() == "plasma") { + phase->setElectronTemperature(Te + dTe); + } phase->getEnthalpy_RT(h2); phase->getCp_R(cp2); + for (size_t k = 0; k < nsp; k++) { + // Determine effective temperature and perturbation for species k. + double T_eff = (k == ke) ? Te : T1; + double dT_eff = (k == ke) ? dTe : dT; + double cp_mid = 0.5 * (cp1[k] + cp2[k]) * GasConstant; - double cp_fd = (h2[k] * (T1 + dT) - h1[k] * T1) / dT * GasConstant; + double cp_fd = (h2[k] * (T_eff + dT_eff) - h1[k] * T_eff) / dT_eff * GasConstant; double tol = max({rtol_fd * std::abs(cp_mid), rtol_fd * std::abs(cp_fd), atol}); EXPECT_NEAR(cp_fd, cp_mid, tol) << "k = " << k; } @@ -765,16 +798,29 @@ TEST_P(TestConsistency, cpRef_eq_dhRefdT) } catch (NotImplementedError& err) { GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; } + // Perturb temperature. double T1 = phase->temperature(); double dT = 1e-5 * phase->temperature(); + // For plasma phases, also perturb the electron temperature. + double Te = phase->electronTemperature(); + double dTe = 1e-5 * Te; + phase->setState_TP(T1 + dT, phase->pressure()); + if (phase->type() == "plasma") { + phase->setElectronTemperature(Te + dTe); + } phase->getEnthalpy_RT_ref(h2); phase->getCp_R_ref(cp2); + for (size_t k = 0; k < nsp; k++) { + // Determine effective temperature and perturbation for species k. + double T_eff = (k == ke) ? Te : T1; + double dT_eff = (k == ke) ? dTe : dT; + double cp_mid = 0.5 * (cp1[k] + cp2[k]) * GasConstant; - double cp_fd = (h2[k] * (T1 + dT) - h1[k] * T1) / dT * GasConstant; + double cp_fd = (h2[k] * (T_eff + dT_eff) - h1[k] * T_eff) / dT_eff * GasConstant; double tol = max({rtol_fd * std::abs(cp_mid), rtol_fd * std::abs(cp_fd), atol}); - EXPECT_NEAR(cp_fd, cp_mid, tol) << "k = " << k; + EXPECT_NEAR(cp_fd, cp_mid, tol) << "k = " << k << "(ke = " << ke << ")"; } }