Add isothermal_compressibility and thermal_expansion_coefficient for Redlich-Kwong and Peng-Robinson EOS

AUTHORS

@@ -29,6 +29,7 @@ John Hewson (@jchewson), Sandia National Laboratory
 Trevor Hickey (@luxe)
 Yuanjie Jiang
 Benjamin Kee (@lionkey)
+Cory Kinney (@corykinney)
 Gandhali Kogekar (@gkogekar)
 Daniel Korff (@korffdm), Colorado School of Mines
 Corey R. Randall (@c-randall), Colorado School of Mines

include/cantera/thermo/PengRobinson.h

@@ -205,6 +205,9 @@ class PengRobinson : public MixtureFugacityTP
 
 public:
 
+    virtual double isothermalCompressibility() const;
+    virtual double thermalExpansionCoeff() const;
+
     //! Calculate \f$dp/dV\f$ and \f$dp/dT\f$ at the current conditions
     /*!
      *  These are stored internally.
@@ -248,7 +251,7 @@ class PengRobinson : public MixtureFugacityTP
      */
     double m_a;
 
-    //! Value of \f$\alpha\f$ in the equation of state
+    //! Value of \f$a \alpha\f$ in the equation of state
     /*!
      *  `m_aAlpha_mix` is a function of the temperature and the mole fractions.
      */

include/cantera/thermo/RedlichKwongMFTP.h

@@ -173,6 +173,9 @@ class RedlichKwongMFTP : public MixtureFugacityTP
     virtual doublereal densSpinodalGas() const;
     virtual doublereal dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal& presCalc) const;
 
+    virtual double isothermalCompressibility() const;
+    virtual double thermalExpansionCoeff() const;
+
     //! Calculate dpdV and dpdT at the current conditions
     /*!
      *  These are stored internally.

src/thermo/PengRobinson.cpp

@@ -630,6 +630,18 @@ double PengRobinson::dpdVCalc(double T, double molarVol, double& presCalc) const
     return -GasConstant * T / (vmb * vmb) + 2 * m_aAlpha_mix * vpb / (denom*denom);
 }
 
+double PengRobinson::isothermalCompressibility() const
+{
+    calculatePressureDerivatives();
+    return -1 / (molarVolume() * m_dpdV);
+}
+
+double PengRobinson::thermalExpansionCoeff() const
+{
+    calculatePressureDerivatives();
+    return -m_dpdT / (molarVolume() * m_dpdV);
+}
+
 void PengRobinson::calculatePressureDerivatives() const
 {
     double T = temperature();

src/thermo/RedlichKwongMFTP.cpp

@@ -707,6 +707,18 @@ doublereal RedlichKwongMFTP::dpdVCalc(doublereal TKelvin, doublereal molarVol, d
     return dpdv;
 }
 
+double RedlichKwongMFTP::isothermalCompressibility() const
+{
+    pressureDerivatives();
+    return -1 / (molarVolume() * dpdV_);
+}
+
+double RedlichKwongMFTP::thermalExpansionCoeff() const
+{
+    pressureDerivatives();
+    return -dpdT_ / (molarVolume() * dpdV_);
+}
+
 void RedlichKwongMFTP::pressureDerivatives() const
 {
     doublereal TKelvin = temperature();

test/thermo/consistency.cpp

@@ -74,6 +74,8 @@ class TestConsistency : public testing::TestWithParam<std::tuple<AnyMap, AnyMap>
         atol = setup.getDouble("atol", 1e-5);
         rtol_fd = setup.getDouble("rtol_fd", 1e-6);
         atol_v = setup.getDouble("atol_v", 1e-11);
+        atol_c = setup.getDouble("atol_c", 1e-14);
+        atol_e = setup.getDouble("atol_e", 1e-18);
 
         phase = cache[key];
         phase->setState(state);
@@ -104,6 +106,8 @@ class TestConsistency : public testing::TestWithParam<std::tuple<AnyMap, AnyMap>
     double T, p, RT;
     double atol; // absolute tolerance for molar energy comparisons
     double atol_v; // absolute tolerance for molar volume comparisons
+    double atol_c; // absolute tolerance for compressibility comparison
+    double atol_e; // absolute tolerance for expansion comparison
     double rtol_fd; // relative tolerance for finite difference comparisons
 };
 
@@ -345,6 +349,58 @@ TEST_P(TestConsistency, dSdv_const_T_eq_dPdT_const_V) {
     }
 }
 
+TEST_P(TestConsistency, betaT_eq_minus_dmv_dP_const_T_div_mv)
+{
+    double betaT1;
+    try {
+        betaT1 = phase->isothermalCompressibility();
+    } catch (NotImplementedError& err) {
+        GTEST_SKIP() << err.getMethod() << " threw NotImplementedError";
+    }
+
+    double T = phase->temperature();
+    double P1 = phase->pressure();
+    double mv1 = phase->molarVolume();
+
+    double P2 = P1 * (1 + 1e-6);
+    phase->setState_TP(T, P2);
+    double betaT2 = phase->isothermalCompressibility();
+    double mv2 = phase->molarVolume();
+
+    double betaT_mid = 0.5 * (betaT1 + betaT2);
+    double mv_mid = 0.5 * (mv1 + mv2);
+    double betaT_fd = -1 / mv_mid * (mv2 - mv1) / (P2 - P1);
+
+    EXPECT_NEAR(betaT_fd, betaT_mid, 
+                max({rtol_fd * betaT_mid, rtol_fd * betaT_fd, atol_c}));
+}
+
+TEST_P(TestConsistency, alphaV_eq_dmv_dT_const_P_div_mv)
+{
+    double alphaV1;
+    try {
+        alphaV1 = phase->thermalExpansionCoeff();
+    } catch (NotImplementedError& err) {
+        GTEST_SKIP() << err.getMethod() << " threw NotImplementedError";
+    }
+
+    double P = phase->pressure();
+    double T1 = phase->temperature();
+    double mv1 = phase->molarVolume();
+
+    double T2 = T1 * (1 + 1e-6);
+    phase->setState_TP(T2, P);
+    double alphaV2 = phase->thermalExpansionCoeff();
+    double mv2 = phase->molarVolume();
+
+    double alphaV_mid = 0.5 * (alphaV1 + alphaV2);
+    double mv_mid = 0.5 * (mv1 + mv2);
+    double alphaV_fd = 1 / mv_mid * (mv2 - mv1) / (T2 - T1);
+
+    EXPECT_NEAR(alphaV_fd, alphaV_mid, 
+                max({rtol_fd * alphaV_mid, rtol_fd * alphaV_fd, atol_e}));
+}
+
 // ---------- Tests for consistency of standard state properties ---------------
 
 TEST_P(TestConsistency, hk0_eq_uk0_plus_p_vk0)