Peng-Robinson EoS: added setState_DP

include/cantera/thermo/IdealGasPhase.h

@@ -366,7 +366,7 @@ class IdealGasPhase: public ThermoPhase
      * @param p Pressure (Pa)
      * @since New in %Cantera 3.0.
      */
-    void setState_DP(double rho, double p) override {
+    void setState_DP(double rho, double p,double tol=1e-9, double Tguess=300) override {
         if (p <= 0) {
             throw CanteraError("IdealGasPhase::setState_DP",
                                "pressure must be positive");

include/cantera/thermo/MixtureFugacityTP.h

@@ -88,7 +88,7 @@ class MixtureFugacityTP : public ThermoPhase
      *       -1 means gas. The value -2 means unrestricted. Values of zero or
      *       greater refer to species dominated condensed phases.
      */
-    void setForcedSolutionBranch(int solnBranch);
+    void setForcedSolutionBranch(int solnBranch) override;
 
     //! Report the solution branch which the solution is restricted to
     /*!

include/cantera/thermo/PengRobinson.h

@@ -203,6 +203,7 @@ class PengRobinson : public MixtureFugacityTP
      *  These are stored internally.
      */
     double daAlpha_dT() const;
+    double daAlpha_dT(double& temp) const;
 
     //! Calculate second derivative @f$ d^2(a \alpha)/dT^2 @f$
     /*!
@@ -211,6 +212,18 @@ class PengRobinson : public MixtureFugacityTP
     double d2aAlpha_dT2() const;
 
 public:
+    //! Set the density (kg/m**3) and pressure (Pa) at constant composition
+    /*!
+     * Within %Cantera, the independent variable of the PR-EoS is the density and Temperature. Therefore, this function solves for
+     * the temperature that will yield the desired input pressure at the provided density.
+     * The composition is held constant during this process.
+     *
+     * @param rho Density (kg/m^3)
+     * @param p   Pressure (Pa)
+     * @param tol   tolerance for the iterative solver 
+     * @param TGuess   guess for the temperature (K) used as start for the iterative solver 
+     */
+    void setState_DP(double rho, double p, double tol=1e-9, double Tguess=300) override;
 
     double isothermalCompressibility() const override;
     double thermalExpansionCoeff() const override;
@@ -229,7 +242,7 @@ class PengRobinson : public MixtureFugacityTP
      *  the internal numbers based on the state of the object.
      */
     void updateMixingExpressions() override;
-
+    void updateMixingExpressions(double& temp,double& aCalc, double& bCalc, double& aAlphaCalc);
     //! Calculate the @f$ a @f$, @f$ b @f$, and @f$ \alpha @f$ parameters given the temperature
     /*!
      * This function doesn't change the internal state of the object, so it is a

include/cantera/thermo/ThermoPhase.h

@@ -1399,12 +1399,22 @@ class ThermoPhase : public Phase
      *
      * @param rho Density (kg/m^3)
      * @param p   Pressure (Pa)
+     * @param tol   tolerance for the iterative solver for cubic EoS
+     * @param TGuess   guess for the temperature (K) used as start for the iterative solver for cubic EoS
      * @since New in %Cantera 3.0.
      */
-    virtual void setState_DP(double rho, double p) {
+    virtual void setState_DP(double rho, double p, double tol=1e-9, double Tguess=300) {
         throw NotImplementedError("ThermoPhase::setState_DP",
                                   "Not implemented for phase type '{}'", type());
     }
+    //! Set the solution branch to force the ThermoPhase to exist on one branch
+    //! or another
+    /*!
+     *  @param solnBranch  Branch that the solution is restricted to. the value
+     *       -1 means gas. The value -2 means unrestricted. Values of zero or
+     *       greater refer to species dominated condensed phases.
+     */
+    virtual void setForcedSolutionBranch(int solnBranch){};
 
     //! Set the state using an AnyMap containing any combination of properties
     //! supported by the thermodynamic model

interfaces/cython/cantera/thermo.pyx

@@ -1486,32 +1486,48 @@ cdef class ThermoPhase(_SolutionBase):
         def __get__(self):
             return self.density, self.P
         def __set__(self, values):
-            assert len(values) == 2, 'incorrect number of values'
+            assert ((len(values) == 2) || (len(values) == 4)), 'incorrect number of values either 2 (DP) or 4 (D, P, tolerance, TGuess)'
             D = values[0] if values[0] is not None else self.density
             P = values[1] if values[1] is not None else self.P
-            self.thermo.setState_DP(D*self._mass_factor(), P)
+            if (len(values) == 2):
+                self.thermo.setState_DP(D*self._mass_factor(), P)
+            else:
+                tolerance = values[2] 
+                TGuess    = values[3]
+                self.thermo.setState_DP(D*self._mass_factor(), P,tolerance,TGuess)
 
     property DPX:
         """Get/Set density [kg/m^3], pressure [Pa], and mole fractions."""
         def __get__(self):
             return self.density, self.P, self.X
         def __set__(self, values):
-            assert len(values) == 3, 'incorrect number of values'
+            assert ((len(values) == 3) || (len(values) == 5)), 'incorrect number of values either 3 (DPX) or 5 (D, P, X, tolerance, TGuess)'
+            if 
             D = values[0] if values[0] is not None else self.density
             P = values[1] if values[1] is not None else self.P
             self.X = values[2]
-            self.thermo.setState_DP(D*self._mass_factor(), P)
+            if (len(values) == 3):
+                self.thermo.setState_DP(D*self._mass_factor(), P)
+            else:
+                tolerance = values[3] 
+                TGuess    = values[4]
+                self.thermo.setState_DP(D*self._mass_factor(), P,tolerance,TGuess)
 
     property DPY:
         """Get/Set density [kg/m^3], pressure [Pa], and mass fractions."""
         def __get__(self):
             return self.density, self.P, self.Y
         def __set__(self, values):
-            assert len(values) == 3, 'incorrect number of values'
+            assert ((len(values) == 3) || (len(values) == 5)), 'incorrect number of values either 3 (DPY) or 5 (D, P, Y, tolerance, TGuess)'
             D = values[0] if values[0] is not None else self.density
             P = values[1] if values[1] is not None else self.P
             self.Y = values[2]
-            self.thermo.setState_DP(D*self._mass_factor(), P)
+            if (len(values) == 3):
+                self.thermo.setState_DP(D*self._mass_factor(), P)
+            else:
+                tolerance = values[3] 
+                TGuess    = values[4]
+                self.thermo.setState_DP(D*self._mass_factor(), P,tolerance,TGuess)
 
     property HP:
         """Get/Set enthalpy [J/kg or J/kmol] and pressure [Pa]."""

src/thermo/MixtureFugacityTP.cpp

@@ -835,7 +835,8 @@ int MixtureFugacityTP::solveCubic(double T, double pres, double a, double b,
     if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) {
         if (disc > 1e-10) {
             throw CanteraError("MixtureFugacityTP::solveCubic",
-                "value of yN and h are too high, unrealistic roots may be obtained");
+                "value of yN and h are too high, unrealistic roots may be obtained"
+                );
         }
         disc = 0.0;
     }

src/thermo/PengRobinson.cpp

@@ -673,6 +673,28 @@ void PengRobinson::updateMixingExpressions()
     calculateAB(m_a,m_b,m_aAlpha_mix);
 }
 
+
+void PengRobinson::updateMixingExpressions(double& temp,double& aCalc, double& bCalc, double& aAlphaCalc)
+{
+    //double temp = temperature();
+
+    // Update individual alpha
+    for (size_t j = 0; j < m_kk; j++) {
+        double critTemp_j = speciesCritTemperature(m_a_coeffs(j,j), m_b_coeffs[j]);
+        double sqt_alpha = 1 + m_kappa[j] * (1 - sqrt(temp / critTemp_j));
+        m_alpha[j] = sqt_alpha*sqt_alpha;
+    }
+
+    // Update aAlpha_i, j
+    for (size_t i = 0; i < m_kk; i++) {
+        for (size_t j = 0; j < m_kk; j++) {
+            m_aAlpha_binary(i, j) = sqrt(m_alpha[i] * m_alpha[j]) * m_a_coeffs(i,j);
+        }
+    }
+    calculateAB(aCalc,bCalc,aAlphaCalc);
+}
+
+
 void PengRobinson::calculateAB(double& aCalc, double& bCalc, double& aAlphaCalc) const
 {
     bCalc = 0.0;
@@ -710,6 +732,30 @@ double PengRobinson::daAlpha_dT() const
     return daAlphadT;
 }
 
+double PengRobinson::daAlpha_dT(double& temp) const
+{
+    double daAlphadT = 0.0, k, Tc, sqtTr, coeff1, coeff2;
+    for (size_t i = 0; i < m_kk; i++) {
+        // Calculate first derivative of alpha for individual species
+        Tc = speciesCritTemperature(m_a_coeffs(i,i), m_b_coeffs[i]);
+        sqtTr = sqrt(temp / Tc);
+        coeff1 = 1 / (Tc*sqtTr);
+        coeff2 = sqtTr - 1;
+        k = m_kappa[i];
+        m_dalphadT[i] = coeff1 * (k*k*coeff2 - k);
+    }
+    // Calculate mixture derivative
+    for (size_t i = 0; i < m_kk; i++) {
+        for (size_t j = 0; j < m_kk; j++) {
+            daAlphadT += moleFractions_[i] * moleFractions_[j] * 0.5
+                         * m_aAlpha_binary(i, j)
+                         * (m_dalphadT[i] / m_alpha[i] + m_dalphadT[j] / m_alpha[j]);
+        }
+    }
+    return daAlphadT;
+}
+
+
 double PengRobinson::d2aAlpha_dT2() const
 {
     for (size_t i = 0; i < m_kk; i++) {
@@ -792,4 +838,45 @@ AnyMap PengRobinson::getAuxiliaryData()
     return parameters;
 }
 
+void PengRobinson::setState_DP(double rho, double p,double tol, double Tguess)
+{
+    if (p <= 0 || rho<=0 ) {
+        throw CanteraError("PengRobinson::setState_DP",
+                            "pressure and density must be positive");
+    }
+    setDensity(rho);
+    double mV =  meanMolecularWeight() / rho;
+    // Update individual alpha
+    double aCalc, bCalc, aAlpha;
+    int iteration = 0;
+    double T_new = Tguess; //(p + term2) * denum / GasConstant;
+    double p_iter=p;
+    while (iteration < 300) {
+            // Update a, b, and alpha based on the current guess for T
+        updateMixingExpressions(T_new, aCalc, bCalc, aAlpha);
+        // Calculate p using PR EoS
+        double denom = mV * mV + 2 * mV * bCalc - bCalc * bCalc;
+        double mV_b = mV - bCalc;
+        double term2 = (aAlpha) / denom;
+        p_iter = (GasConstant*T_new)/mV_b - term2;
+        // Check if the result is within the specified tolerance
+        double error =(p_iter-p)/p;
+        if (fabs(error)<tol) {
+            break;
+        }
+        m_dpdT = (GasConstant / mV_b - daAlpha_dT(T_new) / denom);
+        T_new = T_new - (error/(m_dpdT/p));
+        iteration++;
+    }
+    if (iteration==300)
+    {
+    throw CanteraError("PengRobinson::DP did not converge at",
+        "negative temperature T = {} p= {} rho={} X1={} X2={}", T_new, p, rho,moleFractions_[0],moleFractions_[1]);
+    }
+    setDensity(rho);
+    setTemperature(T_new);
+    setState_TP(T_new,p);
+    updateMixingExpressions();
+}
+
 }

test/thermo/PengRobinson_Test.cpp

@@ -179,6 +179,36 @@ TEST_F(PengRobinson_Test, CoolPropValidate)
     }
 }
 
+TEST_F(PengRobinson_Test, CoolPropValidateDP)
+{
+    // Validate the P-R EoS in Cantera with P-R EoS from CoolProp
+
+    const double rhoCoolProp[10] = {
+        9.067928191884574,
+        8.318322900591179,
+        7.6883521740498155,
+        7.150504298001246,
+        6.685330199667018,
+        6.278630757480957,
+        5.919763108091383,
+        5.600572727499541,
+        5.314694056926007,
+        5.057077678380463
+    };
+
+    double p = 5e5;
+
+    // Calculate Temperature using Peng-Robinson EoS from Cantera
+    for(int i=0; i<10; i++)
+    {
+        const double tempReference = 300 + i*25;
+        set_r(1.0);
+        test_phase->setState_DP(rhoCoolProp[i], p,1e-10,300);
+
+        EXPECT_NEAR(test_phase->temperature(),tempReference,1.e-3);
+    }
+}
+
 TEST(PengRobinson, lookupSpeciesProperties)
 {
     AnyMap phase_def = AnyMap::fromYamlString(