thermo: reject degenerate nonphysical single cubic roots

Author: wandadars

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
@@ -841,6 +843,11 @@ int MixtureFugacityTP::solveCubic(double T, double pres, double a, double b,
     }
 
     double tmp;
+    auto physicalRoot = [b](double v) {
+        // For cubic EoS, physically admissible roots satisfy V > b and V > 0.
+        double vmin = std::max(0.0, b * (1.0 + 1e-12));
+        return std::isfinite(v) && v > vmin;
+    };
     // One real root -> have to determine whether gas or liquid is the root
     if (disc > 0.0) {
         double tmpD = sqrt(disc);
@@ -921,6 +928,11 @@ int MixtureFugacityTP::solveCubic(double T, double pres, double a, double b,
     // Find an accurate root, since there might be a heavy amount of roundoff error due to bad conditioning in this solver.
     double res, dresdV = 0.0;
     for (int i = 0; i < nSolnValues; i++) {
+        if (!physicalRoot(Vroot[i])) {
+            // Non-physical roots (for example, negative molar volume) are not
+            // used for state selection and should not trigger convergence errors.
+            continue;
+        }
         for (int n = 0; n < 20; n++) {
             res = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
             if (fabs(res) < 1.0E-14) { // accurate root is obtained
@@ -946,6 +958,11 @@ int MixtureFugacityTP::solveCubic(double T, double pres, double a, double b,
                                "V = {}", T, pres, Vroot[i]);
         }
     }
+    if (nSolnValues == 1 && !physicalRoot(Vroot[0])) {
+        throw CanteraError("MixtureFugacityTP::solveCubic",
+                           "single real root is non-physical for T = {}, p = {} "
+                           "(V = {}, b = {})", T, pres, Vroot[0], b);
+    }
 
     if (nSolnValues == 1) {
         // Determine the phase of the single root.

test/python/test_thermo.py

@@ -1071,6 +1071,28 @@ def test_partial_lookup_pr_state_setting(self, temperature, pressure, compositio
         assert gas.density > 0.0
         assert gas.P == approx(pressure)
 
+    @pytest.mark.parametrize("p", [2.0e7, 4.0e7, 6.0e7])
+    def test_partial_lookup_pr_high_pressure_equilibrate(self, p):
+        gas = self._partial_lookup_peng_robinson_phase()
+
+        fuel_temp = 255.0
+        oxidizer_temp = 142.0
+
+        gas.TPY = fuel_temp, p, "CH4:1.0"
+        yin_f = gas.Y
+        gas.TPY = oxidizer_temp, p, "O2:1.0"
+        yin_o = gas.Y
+        zst = 1.0 / (1.0 + gas.stoich_air_fuel_ratio(yin_f, yin_o, "mass"))
+        yst = zst * yin_f + (1.0 - zst) * yin_o
+
+        tbar = 0.5 * (fuel_temp + oxidizer_temp)
+        gas.TPY = tbar, p, yst
+        gas.equilibrate("HP")
+
+        assert np.isfinite(gas.T)
+        assert gas.T > tbar
+        assert gas.P == approx(p)
+
     def test_partial_lookup_pr_zero_ab_properties_finite(self):
         gas = self._partial_lookup_peng_robinson_phase()
         gas.TPX = 500.0, 1.0e6, "CO2:0.333, H2O:0.666"

test/thermo/cubicSolver_Test.cpp

@@ -115,5 +115,16 @@ TEST_F(cubicSolver_Test, solve_cubic)
     p = peng_robinson_phase->pressure();
     EXPECT_NEAR(p, pCrit, 1);
 }
+
+TEST_F(cubicSolver_Test, solve_cubic_nonphysical_single_root_when_b_zero)
+{
+    auto* peng_robinson_phase = dynamic_cast<PengRobinson*>(test_phase.get());
+    ASSERT_TRUE(peng_robinson_phase != nullptr);
+
+    double Vroot[3];
+    EXPECT_THROW(
+        peng_robinson_phase->solveCubic(300.0, 1.0e5, 2.0e7, 0.0, 2.0e7, Vroot),
+        CanteraError);
+}
 #endif
 };