diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp
index de62203304207df2c941006b6f21dc3eb458ec06..3a17d69c6e90fc2ba4f71db5ddabd146c651442f 100644
--- a/include/teqp/algorithms/VLE.hpp
+++ b/include/teqp/algorithms/VLE.hpp
@@ -168,6 +168,36 @@ auto pure_VLE_T(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV, int maxi
return do_pure_VLE_T(res, rhoL, rhoV, maxiter);
}
+/***
+ * \brief Calculate the derivative of vapor pressure with respect to temperature
+ * \param model The model to operate on
+ * \param T Temperature
+ * \param rhoL Liquid density
+ * \param rhoV Vapor density
+ *
+ * Based upon
+ * \f[
+ * \frac{dp_{\sigma}}{dT} = \frac{h''-h'}{T(v''-v')} = \frac{s''-s'}{v''-v'}
+ * \f]
+ * where the \f$h''-h'\f$ is given by the difference in residual enthalpy \f$h''-h' = h^r''-h^r'\f$ because the ideal-gas parts cancel
+ */
+template<typename Model, typename Scalar, ADBackends backend = ADBackends::autodiff>
+auto dpsatdT_pure(const Model& model, Scalar T, Scalar rhoL, Scalar rhoV) {
+
+ auto molefrac = (Eigen::ArrayXd(1) << 1.0).finished();
+
+ using tdx = TDXDerivatives<decltype(model), double, decltype(molefrac)>;
+ using iso = IsochoricDerivatives<decltype(model), double, decltype(molefrac)>;
+
+ auto R = model.R(molefrac);
+
+ auto hrVLERTV = tdx::get_Ar01(model, T, rhoV, molefrac) + tdx::get_Ar10(model, T, rhoV, molefrac);
+ auto hrVLERTL = tdx::get_Ar01(model, T, rhoL, molefrac) + tdx::get_Ar10(model, T, rhoL, molefrac);
+ auto deltahr_over_T = R*(hrVLERTV-hrVLERTL);
+ auto dpsatdT = deltahr_over_T/(1/rhoV-1/rhoL); // From Clausius-Clapeyron; dp/dT = Deltas/Deltav = Deltah/(T*Deltav); Delta=V-L
+ return dpsatdT;
+}
+
/***
* \brief Do a vapor-liquid phase equilibrium problem for a mixture (binary only for now) with mole fractions specified in the liquid phase
* \param model The model to operate on
diff --git a/include/teqp/cpp/teqpcpp.hpp b/include/teqp/cpp/teqpcpp.hpp
index ca0cac427e7da0ba4e6ba6c395697dac3bf4388d..b1981235274304f7613177a8c260c7a99061499e 100644
--- a/include/teqp/cpp/teqpcpp.hpp
+++ b/include/teqp/cpp/teqpcpp.hpp
@@ -125,6 +125,7 @@ namespace teqp {
virtual std::tuple<double, double> extrapolate_from_critical(const double Tc, const double rhoc, const double Tgiven) const = 0;
virtual std::tuple<EArrayd, EMatrixd> get_pure_critical_conditions_Jacobian(const double T, const double rho, int alternative_pure_index=-1, int alternative_length=2) const = 0;
virtual EArray2 pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const = 0;
+ virtual double dpsatdT_pure(const double T, const double rhoL, const double rhoV) const = 0;
virtual std::tuple<EArrayd, EArrayd> get_drhovecdp_Tsat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const = 0;
virtual std::tuple<EArrayd, EArrayd> get_drhovecdT_psat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const = 0;
diff --git a/interface/CPP/teqp_impl_VLE.cpp b/interface/CPP/teqp_impl_VLE.cpp
index 017d76127ddab417e610097bb497f70e98ec03be..91be45db837653c24104a13b6c0af5adf376606e 100644
--- a/interface/CPP/teqp_impl_VLE.cpp
+++ b/interface/CPP/teqp_impl_VLE.cpp
@@ -50,3 +50,9 @@ EArray2 MI::pure_VLE_T(const double T, const double rhoL, const double rhoV, int
return teqp::pure_VLE_T(model, T, rhoL, rhoV, maxiter);
}, m_model);
}
+
+double MI::dpsatdT_pure(const double T, const double rhoL, const double rhoV) const {
+ return std::visit([&](const auto& model) {
+ return teqp::dpsatdT_pure(model, T, rhoL, rhoV);
+ }, m_model);
+}
diff --git a/interface/CPP/teqpcpp.cpp b/interface/CPP/teqpcpp.cpp
index d23eaa9ef9a7d5c9d21fc53111c78d489c8db34a..11e268b7a293a6fa7b5f7846f9518beab035baaa 100644
--- a/interface/CPP/teqpcpp.cpp
+++ b/interface/CPP/teqpcpp.cpp
@@ -84,6 +84,7 @@ namespace teqp {
double get_neff(const double T, const double rho, const EArrayd& molefracs) const override;
EArray2 pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const override;
+ double dpsatdT_pure(const double T, const double rhoL, const double rhoV) const override;
std::tuple<EArrayd, EArrayd> get_drhovecdp_Tsat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const override;
std::tuple<EArrayd, EArrayd> get_drhovecdT_psat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const override;
double get_dpsat_dTsat_isopleth(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const override;
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index 6054024223273dfd9aaf907a1a61486a5a4649fc..88ac1f5fc8332dfb1e0695f0bbf8637631522e62 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -347,6 +347,7 @@ void init_teqp(py::module& m) {
.def("get_dp_dT_crit", &am::get_dp_dT_crit, "T"_a, "rhovec"_a.noconvert())
.def("pure_VLE_T", &am::pure_VLE_T, "T"_a, "rhoL"_a, "rhoV"_a, "max_iter"_a)
+ .def("dpsatdT_pure", &am::dpsatdT_pure, "T"_a, "rhoL"_a, "rhoV"_a)
.def("get_drhovecdp_Tsat", &am::get_drhovecdp_Tsat, "T"_a, "rhovecL"_a.noconvert(), "rhovecV"_a.noconvert())
.def("get_drhovecdT_psat", &am::get_drhovecdT_psat, "T"_a, "rhovecL"_a.noconvert(), "rhovecV"_a.noconvert())
diff --git a/src/tests/catch_test_cubics.cxx b/src/tests/catch_test_cubics.cxx
index 6141d2745b1b309d6f5d53e663b5bda0b136d9a4..531cca9dd3c9e1299db15f72fd595fb6c78a4b4a 100644
--- a/src/tests/catch_test_cubics.cxx
+++ b/src/tests/catch_test_cubics.cxx
@@ -80,6 +80,46 @@ TEST_CASE("Check calling superancillary curves", "[cubic][superanc]")
}
}
+TEST_CASE("Check orthobaric density derivatives for pure fluid", "[cubic][superanc]")
+{
+ std::valarray<double> Tc_K = { 150.687 };
+ std::valarray<double> pc_Pa = { 4863000.0};
+ std::valarray<double> acentric = { 0.0};
+
+ double T = 130.0, dT = 0.001;
+ auto molefrac = (Eigen::ArrayXd(1) << 1.0).finished();
+
+ auto model = canonical_PR(Tc_K, pc_Pa, acentric);
+ using tdx = TDXDerivatives<decltype(model)>;
+ using iso = IsochoricDerivatives<decltype(model)>;
+
+ auto R = model.R(molefrac);
+ auto [rhoL, rhoV] = model.superanc_rhoLV(T);
+ CHECK(rhoL > rhoV);
+
+ // Finite difference test
+ auto [rhoLp, rhoVp] = model.superanc_rhoLV(T+dT);
+ auto [rhoLm, rhoVm] = model.superanc_rhoLV(T-dT);
+ auto pLp = rhoLp*R*(T+dT) + iso::get_pr(model, T+dT, rhoLp*molefrac);
+ auto pLm = rhoLm*R*(T-dT) + iso::get_pr(model, T-dT, rhoLm*molefrac);
+
+ // Exact solution for density derivative
+ // Change in enthalpy (Deltah) is equal to change in residual enthalpy (Deltahr) because ideal parts cancel
+ auto hrVLERTV = tdx::get_Ar01(model, T, rhoV, molefrac) + tdx::get_Ar10(model, T, rhoV, molefrac);
+ auto hrVLERTL = tdx::get_Ar01(model, T, rhoL, molefrac) + tdx::get_Ar10(model, T, rhoL, molefrac);
+ auto deltahr_over_T = R*(hrVLERTV-hrVLERTL);
+ auto dpsatdT = deltahr_over_T/(1/rhoV-1/rhoL); // From Clausius-Clapeyron; dp/dT = Deltas/Deltav = Deltah/(T*Deltav); Delta=V-L
+
+ auto dpsatdT_routine = dpsatdT_pure(model, T, rhoL, rhoV);
+
+ CHECK(dpsatdT == Approx((pLp - pLm)/(2*dT)));
+ CHECK(dpsatdT_routine == Approx((pLp - pLm)/(2*dT)));
+
+// CHECK(drhovecdTL(0) == Approx((rhoLp-rhoLm)/(2*dT)));
+// auto drhoLdT =
+
+}
+
TEST_CASE("Check manual integration of subcritical VLE isotherm for binary mixture", "[cubic][isochoric][traceisotherm]")
{
using namespace boost::numeric::odeint;
@@ -389,4 +429,4 @@ TEST_CASE("Check manual integration of subcritical VLE isobar for binary mixture
std::ofstream file("isoP.json"); file << J;
}
-}
\ No newline at end of file
+}