From 869f06671e1c732feb33bc8cc585af22740dc08b Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Fri, 15 Dec 2023 18:03:57 -0500 Subject: [PATCH] Test the fixed CPA where radial distribution function is specified Closes #81 --- include/teqp/models/CPA.hpp | 22 ++++++++++++++++------ src/tests/catch_tests.cxx | 31 ++++++++++++++++++++++++++++++- 2 files changed, 46 insertions(+), 7 deletions(-) diff --git a/include/teqp/models/CPA.hpp b/include/teqp/models/CPA.hpp index 2e3cf8a2..21199e0a 100644 --- a/include/teqp/models/CPA.hpp +++ b/include/teqp/models/CPA.hpp @@ -26,6 +26,15 @@ inline auto get_association_classes(const std::string& s) { enum class radial_dist { CS, KG, OT }; +inline auto get_radial_dist(const std::string& s) { + if (s == "CS") { return radial_dist::CS; } + else if (s == "KG") { return radial_dist::KG; } + else if (s == "OT") { return radial_dist::OT; } + else { + throw std::invalid_argument("bad association flag:" + s); + } +} + /// Function that calculates the association binding strength between site A of molecule i and site B on molecule j template inline auto get_DeltaAB_pure(radial_dist dist, double epsABi, double betaABi, BType b_cubic, TType RT, RhoType rhomolar, const VecType& molefrac) { @@ -69,7 +78,7 @@ inline auto get_DeltaAB_pure(radial_dist dist, double epsABi, double betaABi, BT /// template -inline auto XA_calc_pure(int N_sites, association_classes scheme, double epsABi, double betaABi, const BType b_cubic, const TType RT, const RhoType rhomolar, const VecType& molefrac) { +inline auto XA_calc_pure(int N_sites, association_classes scheme, radial_dist dist, double epsABi, double betaABi, const BType b_cubic, const TType RT, const RhoType rhomolar, const VecType& molefrac) { // Matrix XA(A, j) that contains all of the fractions of sites A not bonded to other active sites for each molecule i // Start values for the iteration(set all sites to non - bonded, = 1) @@ -79,7 +88,6 @@ inline auto XA_calc_pure(int N_sites, association_classes scheme, double epsABi, XA.setOnes(); // Get the association strength between the associating sites - auto dist = radial_dist::KG; // TODO: pass this in auto DeltaAiBj = get_DeltaAB_pure(dist, epsABi, betaABi, b_cubic, RT, rhomolar, molefrac); if (scheme == association_classes::a1A) { // Acids @@ -185,6 +193,7 @@ class CPAAssociation { private: const Cubic cubic; const std::vector classes; + const radial_dist dist; const std::valarray epsABi, betaABi; const std::vector N_sites; const double R_gas; @@ -207,8 +216,8 @@ class CPAAssociation { } public: - CPAAssociation(const Cubic &&cubic, const std::vector& classes, const std::valarray &epsABi, const std::valarray &betaABi, double R_gas) - : cubic(cubic), classes(classes), epsABi(epsABi), betaABi(betaABi), N_sites(get_N_sites(classes)), R_gas(R_gas) {}; + CPAAssociation(const Cubic &&cubic, const std::vector& classes, const radial_dist dist, const std::valarray &epsABi, const std::valarray &betaABi, double R_gas) + : cubic(cubic), classes(classes), dist(dist), epsABi(epsABi), betaABi(betaABi), N_sites(get_N_sites(classes)), R_gas(R_gas) {}; template auto alphar(const TType& T, const RhoType& rhomolar, const VecType& molefrac) const { @@ -217,7 +226,7 @@ class CPAAssociation { // Calculate the fraction of sites not bonded with other active sites auto RT = forceeval(R_gas * T); // R times T - auto XA = XA_calc_pure(N_sites[0], classes[0], epsABi[0], betaABi[0], b_cubic, RT, rhomolar, molefrac); + auto XA = XA_calc_pure(N_sites[0], classes[0], dist, epsABi[0], betaABi[0], b_cubic, RT, rhomolar, molefrac); using return_type = std::common_type_t; return_type alpha_r_asso = 0.0; @@ -280,6 +289,7 @@ inline auto CPAfactory(const nlohmann::json &j){ auto build_assoc = [](const auto &&cubic, const auto& j) { auto N = j["pures"].size(); std::vector classes; + radial_dist dist = get_radial_dist(j.at("radial_dist")); std::valarray epsABi(N), betaABi(N); std::size_t i = 0; for (auto p : j["pures"]) { @@ -288,7 +298,7 @@ inline auto CPAfactory(const nlohmann::json &j){ classes.push_back(get_association_classes(p["class"])); i++; } - return CPAAssociation(std::move(cubic), classes, epsABi, betaABi, j["R_gas / J/mol/K"]); + return CPAAssociation(std::move(cubic), classes, dist, epsABi, betaABi, j["R_gas / J/mol/K"]); }; return CPAEOS(build_cubic(j), build_assoc(build_cubic(j), j)); } diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx index 701a4911..ba19bbda 100644 --- a/src/tests/catch_tests.cxx +++ b/src/tests/catch_tests.cxx @@ -424,6 +424,35 @@ TEST_CASE("Test pure VLE with non-unity R0/Rr", "") { auto rr = 0; } +TEST_CASE("Test water Clapeyron.jl", "[CPA]") { + std::valarray a0 = {0.12277}, bi = {0.0000145}, c1 = {0.6736}, Tc = {647.13}, + molefrac = {1.0}; + auto R = 8.31446261815324; + CPA::CPACubic cub(CPA::cubic_flag::SRK, a0, bi, c1, Tc, R); + double T = 303.15, rhomolar = 1/1.7915123921401366e-5; + + auto z = (Eigen::ArrayXd(1) << 1).finished(); + + using tdx = TDXDerivatives; + auto alphar = cub.alphar(T, rhomolar, molefrac); + CHECK(alphar == Approx(-1.2713135319123854)); // Value from Clapeyron.jl + double p_noassoc = T*rhomolar*R*(1+tdx::get_Ar01(cub, T, rhomolar, z)); + CAPTURE(p_noassoc); + + std::vector schemes = { CPA::association_classes::a4C }; + std::valarray epsAB = { 16655 }, betaAB = { 0.0692 }; + CPA::radial_dist dist = CPA::radial_dist::CS; + CPA::CPAAssociation cpaa(std::move(cub), schemes, dist, epsAB, betaAB, R); + + CPA::CPAEOS cpa(std::move(cub), std::move(cpaa)); + using tdc = TDXDerivatives; + auto Ar01 = tdc::get_Ar01(cpa, T, rhomolar, z); + double p_withassoc = T*rhomolar*R*(1 + Ar01); + CAPTURE(p_withassoc); + + REQUIRE(p_withassoc == Approx(100000.000)); +} + TEST_CASE("Test water", "[CPA]") { std::valarray a0 = {0.12277}, bi = {0.000014515}, c1 = {0.67359}, Tc = {647.096}, molefrac = {1.0}; @@ -441,7 +470,7 @@ TEST_CASE("Test water", "[CPA]") { std::vector schemes = { CPA::association_classes::a4C }; std::valarray epsAB = { 16655 }, betaAB = { 0.0692 }; - CPA::CPAAssociation cpaa(std::move(cub), schemes, epsAB, betaAB, R); + CPA::CPAAssociation cpaa(std::move(cub), schemes, CPA::radial_dist::KG, epsAB, betaAB, R); CPA::CPAEOS cpa(std::move(cub), std::move(cpaa)); using tdc = TDXDerivatives;