From f807a018a297419a1b910143febeb32274644837 Mon Sep 17 00:00:00 2001 From: Cory Kinney Date: Thu, 12 Jan 2023 15:38:49 -0500 Subject: [PATCH 1/3] [Test] Add thermo consistency tests for compressibility functions Add finite difference checks for `isothermalCompressibility` and `thermalExpansionCoeff` for thermo models. --- test/thermo/consistency.cpp | 56 +++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/test/thermo/consistency.cpp b/test/thermo/consistency.cpp index 0004d242120..18a8d10288d 100644 --- a/test/thermo/consistency.cpp +++ b/test/thermo/consistency.cpp @@ -74,6 +74,8 @@ class TestConsistency : public testing::TestWithParam 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 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) From 931de2450e34226a431846a07d7db075939cd7ec Mon Sep 17 00:00:00 2001 From: Cory Kinney Date: Sat, 31 Dec 2022 15:07:22 -0500 Subject: [PATCH 2/3] [Thermo] Add real gas compressibility functions Add `isothermalCompressibility` and `thermalExpansionCoeff` definitions to the `RedlichKwongMFTP` and `PengRobinson` models --- include/cantera/thermo/PengRobinson.h | 5 ++++- include/cantera/thermo/RedlichKwongMFTP.h | 3 +++ src/thermo/PengRobinson.cpp | 12 ++++++++++++ src/thermo/RedlichKwongMFTP.cpp | 12 ++++++++++++ 4 files changed, 31 insertions(+), 1 deletion(-) diff --git a/include/cantera/thermo/PengRobinson.h b/include/cantera/thermo/PengRobinson.h index 6680b1087e1..ace5722707b 100644 --- a/include/cantera/thermo/PengRobinson.h +++ b/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. */ diff --git a/include/cantera/thermo/RedlichKwongMFTP.h b/include/cantera/thermo/RedlichKwongMFTP.h index f6ceba56d59..4c2d42a3211 100644 --- a/include/cantera/thermo/RedlichKwongMFTP.h +++ b/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. diff --git a/src/thermo/PengRobinson.cpp b/src/thermo/PengRobinson.cpp index 913a392bac1..2fe5fffe49a 100644 --- a/src/thermo/PengRobinson.cpp +++ b/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(); diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp index 3267af53572..8c615630fe5 100644 --- a/src/thermo/RedlichKwongMFTP.cpp +++ b/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(); From 79da95aafcaa588a2e06b79443e89e2ecbaa2e00 Mon Sep 17 00:00:00 2001 From: corykinney Date: Sat, 13 May 2023 14:43:08 -0400 Subject: [PATCH 3/3] Update AUTHORS --- AUTHORS | 1 + 1 file changed, 1 insertion(+) diff --git a/AUTHORS b/AUTHORS index e7078c2ddcd..86f59b2e80e 100644 --- a/AUTHORS +++ b/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