address degenerate case of a=0 and b=0 for cubic EoS and added ideal gas fallback methods when degenerate case happens

wandadars · wandadars · commit 432eb0a03a14 · 2026-02-22T23:25:24.000-05:00
diff --git a/include/cantera/thermo/MixtureFugacityTP.h b/include/cantera/thermo/MixtureFugacityTP.h
@@ -311,6 +311,42 @@ class MixtureFugacityTP : public ThermoPhase
      */
     double z() const;
 
+    //! @name Ideal-Limit Utilities for Cubic EoS
+    //!
+    //! These helpers provide shared ideal-gas-limit expressions used by
+    //! cubic-EoS implementations when mixture non-ideal terms vanish.
+    //! @{
+
+    //! Molar heat capacity at constant pressure in the ideal-gas limit.
+    //! Units: J/kmol/K.
+    double cp_mole_ideal() const;
+
+    //! Molar heat capacity at constant volume in the ideal-gas limit.
+    //! Units: J/kmol/K.
+    double cv_mole_ideal() const;
+
+    //! Standard concentration in the ideal-gas limit for species @p k.
+    //! Units: kmol/m^3.
+    double standardConcentration_ideal(size_t k) const;
+
+    //! Ideal-gas-limit activity coefficients for all species.
+    //! Sets all entries of @p ac to 1.
+    void getActivityCoefficients_ideal(double* ac) const;
+
+    //! Ideal-gas-limit chemical potentials for all species.
+    //! @param mu Output vector of chemical potentials [J/kmol].
+    void getChemPotentials_ideal(double* mu) const;
+
+    //! Ideal-gas-limit partial molar enthalpies for all species.
+    //! @param hbar Output vector of partial molar enthalpies [J/kmol].
+    void getPartialMolarEnthalpies_ideal(double* hbar) const;
+
+    //! Ideal-gas-limit partial molar entropies for all species.
+    //! @param sbar Output vector of partial molar entropies [J/kmol/K].
+    void getPartialMolarEntropies_ideal(double* sbar) const;
+
+    //! @}
+
     //! Calculate the deviation terms for the total entropy of the mixture from
     //! the ideal gas mixture
     /**
diff --git a/src/thermo/MixtureFugacityTP.cpp b/src/thermo/MixtureFugacityTP.cpp
@@ -13,6 +13,8 @@
 #include "cantera/base/utilities.h"
 #include "cantera/base/global.h"
 
+#include <algorithm>
+
 using namespace std;
 
 namespace Cantera
@@ -260,6 +262,56 @@ double MixtureFugacityTP::z() const
     return pressure() * meanMolecularWeight() / (density() * RT());
 }
 
+double MixtureFugacityTP::cp_mole_ideal() const
+{
+    return GasConstant * mean_X(m_cp0_R);
+}
+
+double MixtureFugacityTP::cv_mole_ideal() const
+{
+    return GasConstant * (mean_X(m_cp0_R) - 1.0);
+}
+
+double MixtureFugacityTP::standardConcentration_ideal(size_t k) const
+{
+    getStandardVolumes(m_workS.data());
+    return 1.0 / m_workS[k];
+}
+
+void MixtureFugacityTP::getActivityCoefficients_ideal(double* ac) const
+{
+    for (size_t k = 0; k < m_kk; k++) {
+        ac[k] = 1.0;
+    }
+}
+
+void MixtureFugacityTP::getChemPotentials_ideal(double* mu) const
+{
+    double RT_ = RT();
+    getStandardChemPotentials(mu);
+    for (size_t k = 0; k < m_kk; k++) {
+        double xx = std::max(SmallNumber, moleFraction(k));
+        mu[k] += RT_ * log(xx);
+    }
+}
+
+void MixtureFugacityTP::getPartialMolarEnthalpies_ideal(double* hbar) const
+{
+    getEnthalpy_RT_ref(hbar);
+    scale(hbar, hbar + m_kk, hbar, RT());
+}
+
+void MixtureFugacityTP::getPartialMolarEntropies_ideal(double* sbar) const
+{
+    getEntropy_R_ref(sbar);
+    scale(sbar, sbar + m_kk, sbar, GasConstant);
+    double logPres = log(pressure() / refPressure());
+    for (size_t k = 0; k < m_kk; k++) {
+        double xx = std::max(SmallNumber, moleFraction(k));
+        sbar[k] -= GasConstant * (log(xx) + logPres);
+    }
+}
+
 double MixtureFugacityTP::sresid() const
 {
     throw NotImplementedError("MixtureFugacityTP::sresid");
diff --git a/src/thermo/PengRobinson.cpp b/src/thermo/PengRobinson.cpp
@@ -15,6 +15,14 @@
 namespace Cantera
 {
 
+namespace
+{
+bool isIdealCubicLimit(double aMix, double bMix)
+{
+    return aMix == 0.0 && bMix == 0.0;
+}
+}
+
 const double PengRobinson::omega_a = 4.5723552892138218E-01;
 const double PengRobinson::omega_b = 7.77960739038885E-02;
 const double PengRobinson::omega_vc = 3.07401308698703833E-01;
@@ -84,6 +92,9 @@ void PengRobinson::setBinaryCoeffs(const string& species_i,
 double PengRobinson::cp_mole() const
 {
     _updateReferenceStateThermo();
+    if (isIdealCubicLimit(m_aAlpha_mix, m_b)) {
+        return cp_mole_ideal();
+    }
     double T = temperature();
     double mv = molarVolume();
     double vpb = mv + (1 + Sqrt2) * m_b;
@@ -98,6 +109,9 @@ double PengRobinson::cp_mole() const
 double PengRobinson::cv_mole() const
 {
     _updateReferenceStateThermo();
+    if (isIdealCubicLimit(m_aAlpha_mix, m_b)) {
+        return cv_mole_ideal();
+    }
     double T = temperature();
     calculatePressureDerivatives();
     return (cp_mole() + T * m_dpdT * m_dpdT / m_dpdV);
@@ -115,12 +129,15 @@ double PengRobinson::pressure() const
 
 double PengRobinson::standardConcentration(size_t k) const
 {
-    getStandardVolumes(m_workS.data());
-    return 1.0 / m_workS[k];
+    return standardConcentration_ideal(k);
 }
 
 void PengRobinson::getActivityCoefficients(double* ac) const
 {
+    if (isIdealCubicLimit(m_aAlpha_mix, m_b)) {
+        getActivityCoefficients_ideal(ac);
+        return;
+    }
     double mv = molarVolume();
     double vpb2 = mv + (1 + Sqrt2) * m_b;
     double vmb2 = mv + (1 - Sqrt2) * m_b;
@@ -154,8 +171,13 @@ void PengRobinson::getActivityCoefficients(double* ac) const
 
 void PengRobinson::getChemPotentials(double* mu) const
 {
-    getGibbs_ref(mu);
+    if (isIdealCubicLimit(m_aAlpha_mix, m_b)) {
+        getChemPotentials_ideal(mu);
+        return;
+    }
     double RT_ = RT();
+
+    getGibbs_ref(mu);
     for (size_t k = 0; k < m_kk; k++) {
         double xx = std::max(SmallNumber, moleFraction(k));
         mu[k] += RT_ * (log(xx));
@@ -189,6 +211,10 @@ void PengRobinson::getChemPotentials(double* mu) const
 
 void PengRobinson::getPartialMolarEnthalpies(double* hbar) const
 {
+    if (isIdealCubicLimit(m_aAlpha_mix, m_b)) {
+        getPartialMolarEnthalpies_ideal(hbar);
+        return;
+    }
     // First we get the reference state contributions
     getEnthalpy_RT_ref(hbar);
     scale(hbar, hbar+m_kk, hbar, RT());
@@ -239,6 +265,11 @@ void PengRobinson::getPartialMolarEnthalpies(double* hbar) const
 
 void PengRobinson::getPartialMolarEntropies(double* sbar) const
 {
+    if (isIdealCubicLimit(m_aAlpha_mix, m_b)) {
+        getPartialMolarEntropies_ideal(sbar);
+        return;
+    }
+
     // Using the identity : (hk - T*sk) = gk
     double T = temperature();
     getPartialMolarEnthalpies(sbar);
@@ -443,6 +474,9 @@ void PengRobinson::getSpeciesParameters(const string& name, AnyMap& speciesNode)
 
 double PengRobinson::sresid() const
 {
+    if (isIdealCubicLimit(m_aAlpha_mix, m_b)) {
+        return 0.0;
+    }
     double molarV = molarVolume();
     double hh = m_b / molarV;
     double zz = z();
@@ -456,6 +490,9 @@ double PengRobinson::sresid() const
 
 double PengRobinson::hresid() const
 {
+    if (isIdealCubicLimit(m_aAlpha_mix, m_b)) {
+        return 0.0;
+    }
     double molarV = molarVolume();
     double zz = z();
     double aAlpha_1 = daAlpha_dT();
diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp
@@ -15,6 +15,14 @@
 namespace Cantera
 {
 
+namespace
+{
+bool isIdealCubicLimit(double aMix, double bMix)
+{
+    return aMix == 0.0 && bMix == 0.0;
+}
+}
+
 const double RedlichKwongMFTP::omega_a = 4.27480233540E-01;
 const double RedlichKwongMFTP::omega_b = 8.66403499650E-02;
 const double RedlichKwongMFTP::omega_vc = 3.33333333333333E-01;
@@ -86,6 +94,9 @@ void RedlichKwongMFTP::setBinaryCoeffs(const string& species_i, const string& sp
 double RedlichKwongMFTP::cp_mole() const
 {
     _updateReferenceStateThermo();
+    if (isIdealCubicLimit(m_a_current, m_b_current)) {
+        return cp_mole_ideal();
+    }
     double TKelvin = temperature();
     double sqt = sqrt(TKelvin);
     double mv = molarVolume();
@@ -102,6 +113,9 @@ double RedlichKwongMFTP::cp_mole() const
 double RedlichKwongMFTP::cv_mole() const
 {
     _updateReferenceStateThermo();
+    if (isIdealCubicLimit(m_a_current, m_b_current)) {
+        return cv_mole_ideal();
+    }
     double TKelvin = temperature();
     double sqt = sqrt(TKelvin);
     double mv = molarVolume();
@@ -126,12 +140,15 @@ double RedlichKwongMFTP::pressure() const
 
 double RedlichKwongMFTP::standardConcentration(size_t k) const
 {
-    getStandardVolumes(m_workS.data());
-    return 1.0 / m_workS[k];
+    return standardConcentration_ideal(k);
 }
 
 void RedlichKwongMFTP::getActivityCoefficients(double* ac) const
 {
+    if (isIdealCubicLimit(m_a_current, m_b_current)) {
+        getActivityCoefficients_ideal(ac);
+        return;
+    }
     double mv = molarVolume();
     double sqt = sqrt(temperature());
     double vpb = mv + m_b_current;
@@ -164,6 +181,11 @@ void RedlichKwongMFTP::getActivityCoefficients(double* ac) const
 
 void RedlichKwongMFTP::getChemPotentials(double* mu) const
 {
+    if (isIdealCubicLimit(m_a_current, m_b_current)) {
+        getChemPotentials_ideal(mu);
+        return;
+    }
+
     getGibbs_ref(mu);
     for (size_t k = 0; k < m_kk; k++) {
         double xx = std::max(SmallNumber, moleFraction(k));
@@ -198,6 +220,10 @@ void RedlichKwongMFTP::getChemPotentials(double* mu) const
 
 void RedlichKwongMFTP::getPartialMolarEnthalpies(double* hbar) const
 {
+    if (isIdealCubicLimit(m_a_current, m_b_current)) {
+        getPartialMolarEnthalpies_ideal(hbar);
+        return;
+    }
     // First we get the reference state contributions
     getEnthalpy_RT_ref(hbar);
     scale(hbar, hbar+m_kk, hbar, RT());
@@ -243,6 +269,11 @@ void RedlichKwongMFTP::getPartialMolarEnthalpies(double* hbar) const
 
 void RedlichKwongMFTP::getPartialMolarEntropies(double* sbar) const
 {
+    if (isIdealCubicLimit(m_a_current, m_b_current)) {
+        getPartialMolarEntropies_ideal(sbar);
+        return;
+    }
+
     getEntropy_R_ref(sbar);
     scale(sbar, sbar+m_kk, sbar, GasConstant);
     double TKelvin = temperature();
@@ -504,6 +535,9 @@ void RedlichKwongMFTP::getSpeciesParameters(const string& name,
 
 double RedlichKwongMFTP::sresid() const
 {
+    if (isIdealCubicLimit(m_a_current, m_b_current)) {
+        return 0.0;
+    }
     // note this agrees with tpx
     double rho = density();
     double mmw = meanMolecularWeight();
@@ -520,6 +554,9 @@ double RedlichKwongMFTP::sresid() const
 
 double RedlichKwongMFTP::hresid() const
 {
+    if (isIdealCubicLimit(m_a_current, m_b_current)) {
+        return 0.0;
+    }
     // note this agrees with tpx
     double rho = density();
     double mmw = meanMolecularWeight();
diff --git a/test/python/test_thermo.py b/test/python/test_thermo.py
