From 9caaaef6b974203794360c8802aa32b5f1ba6006 Mon Sep 17 00:00:00 2001 From: Ian Bell <ian.bell@nist.gov> Date: Fri, 14 Oct 2022 12:11:55 -0400 Subject: [PATCH] Conversion routines from CoolProp formatted ideal-gas (#19) Add conversion routines for ideal-gas Helmholtz terms from CoolProp format --- include/teqp/ideal_eosterms.hpp | 326 +++++++++++++++++++++++++--- include/teqp/json_tools.hpp | 3 +- interface/pybind11_wrapper.cpp | 3 +- src/tests/catch_test_multifluid.cxx | 25 ++- 4 files changed, 327 insertions(+), 30 deletions(-) diff --git a/include/teqp/ideal_eosterms.hpp b/include/teqp/ideal_eosterms.hpp index 4ef0086..93fc897 100644 --- a/include/teqp/ideal_eosterms.hpp +++ b/include/teqp/ideal_eosterms.hpp @@ -1,8 +1,10 @@ #pragma once #include <variant> +#include <filesystem> #include "teqp/types.hpp" #include "teqp/exceptions.hpp" +#include "teqp/json_tools.hpp" namespace teqp { @@ -11,13 +13,13 @@ namespace teqp { */ class IdealHelmholtzConstant { public: - const double a; - IdealHelmholtzConstant(double a) : a(a) {}; + const double a, R; + IdealHelmholtzConstant(double a, double R) : a(a), R(R) {}; template<typename TType, typename RhoType> auto alphaig(const TType& T, const RhoType& rho) const { using otype = std::common_type_t <TType, RhoType>; - return static_cast<otype>(a); + return forceeval(static_cast<otype>(a)); } }; @@ -32,8 +34,8 @@ namespace teqp { */ class IdealHelmholtzLogT { public: - const double a; - IdealHelmholtzLogT(double a) : a(a) {}; + const double a, R; + IdealHelmholtzLogT(double a, double R) : a(a), R(R) {}; template<typename TType, typename RhoType> auto alphaig(const TType& T, const RhoType& rho) const { @@ -55,8 +57,8 @@ namespace teqp { */ class IdealHelmholtzLead { public: - const double a_1, a_2; - IdealHelmholtzLead(double a_1, double a_2) : a_1(a_1), a_2(a_2) {}; + const double a_1, a_2, R; + IdealHelmholtzLead(double a_1, double a_2, double R) : a_1(a_1), a_2(a_2), R(R) {}; template<typename TType, typename RhoType> auto alphaig(const TType& T, const RhoType& rho) const { @@ -65,13 +67,31 @@ namespace teqp { } }; + /** + \f$ \alpha^{\rm ig}= \sum_k n_kT^{t_k} \f$ + */ + class IdealHelmholtzPowerT { + public: + const std::valarray<double> n, t, R; + IdealHelmholtzPowerT(const std::valarray<double>& n, const std::valarray<double>& t, double R) : n(n), t(t), R(R) {}; + + template<typename TType, typename RhoType> + auto alphaig(const TType& T, const RhoType& rho) const { + std::common_type_t <TType, RhoType> summer = 0.0; + for (auto i = 0; i < n.size(); ++i) { + summer = summer + n[i] * pow(T, t[i]); + } + return forceeval(summer); + } + }; + /** \f$ \alpha^{\rm ig}= \sum_k n_k\ln(1-\exp(-\theta_k/T)) \f$ */ class IdealHelmholtzPlanckEinstein { public: - const std::valarray<double> n, theta; - IdealHelmholtzPlanckEinstein(const std::valarray<double>& n, const std::valarray<double>& theta) : n(n), theta(theta) {}; + const std::valarray<double> n, theta, R; + IdealHelmholtzPlanckEinstein(const std::valarray<double>& n, const std::valarray<double>& theta, double R) : n(n), theta(theta), R(R) {}; template<typename TType, typename RhoType> auto alphaig(const TType& T, const RhoType& rho) const { @@ -83,6 +103,30 @@ namespace teqp { } }; + /** + \f$ \alpha^{\rm ig}= \sum_k n_k\ln(c_k+d_k\exp(\theta_k/T)) \f$ + */ + class IdealHelmholtzPlanckEinsteinGeneralized { + public: + const std::valarray<double> n, c, d, theta, R; + IdealHelmholtzPlanckEinsteinGeneralized( + const std::valarray<double>& n, + const std::valarray<double>& c, + const std::valarray<double>& d, + const std::valarray<double>& theta, + double R + ) : n(n), c(c), d(d), theta(theta), R(R) {}; + + template<typename TType, typename RhoType> + auto alphaig(const TType& T, const RhoType& rho) const { + std::common_type_t <TType, RhoType> summer = 0.0; + for (auto i = 0; i < n.size(); ++i) { + summer = summer + n[i] * log(c[i] + d[i]*exp(theta[i] / T)); + } + return forceeval(summer); + } + }; + /** \f$ \alpha^{\rm ig}= \sum_k n_k ln(|cosh(theta_k/T)|) \f$ @@ -90,8 +134,8 @@ namespace teqp { */ class IdealHelmholtzGERG2004Cosh { public: - const std::valarray<double> n, theta; - IdealHelmholtzGERG2004Cosh(const std::valarray<double>& n, const std::valarray<double>& theta) : n(n), theta(theta) {}; + const std::valarray<double> n, theta, R; + IdealHelmholtzGERG2004Cosh(const std::valarray<double>& n, const std::valarray<double>& theta, double R) : n(n), theta(theta), R(R) {}; template<typename TType, typename RhoType> auto alphaig(const TType& T, const RhoType& rho) const { @@ -110,8 +154,8 @@ namespace teqp { */ class IdealHelmholtzGERG2004Sinh { public: - const std::valarray<double> n, theta; - IdealHelmholtzGERG2004Sinh(const std::valarray<double>& n, const std::valarray<double>& theta) : n(n), theta(theta) {}; + const std::valarray<double> n, theta, R; + IdealHelmholtzGERG2004Sinh(const std::valarray<double>& n, const std::valarray<double>& theta, double R) : n(n), theta(theta), R(R) {}; template<typename TType, typename RhoType> auto alphaig(const TType& T, const RhoType& rho) const { @@ -123,14 +167,60 @@ namespace teqp { } }; + /** + \f$ \alpha^{\rm ig}= c*\left( \frac{T-T0}{T}-\ln\left(T/T0\right)\right) \f$ + + from a term that is like + + \f$ \frac{c_{p0}}{R}= c \f$ + */ + class IdealHelmholtzCp0Constant { + public: + const double c, T_0, R; + IdealHelmholtzCp0Constant( + const double c, const double T_0, const double R + ) : c(c), T_0(T_0), R(R) {}; + + template<typename TType, typename RhoType> + auto alphaig(const TType& T, const RhoType& rho) const { + using otype = std::common_type_t <TType, RhoType>; + return forceeval(static_cast<otype>(c*((T-T_0)/T-log(T/T_0)))); + } + }; + + /** + \f$ \alpha^{\rm ig}= c\left[T^{t}\left(\frac{1}{t+1}-\frac{1}{t}\right)-\frac{T_0^{t+1}}{T(t+1)}+\frac{T_0^t}{t}\right] \f$ + + from a term that is like + + \f$ \frac{c_{p0}}{R}= cT^t, t \neq 0 \f$ + */ + class IdealHelmholtzCp0PowerT { + public: + const double c, t, T_0, R; + IdealHelmholtzCp0PowerT( + const double c, const double t, const double T_0, const double R + ) : c(c), t(t), T_0(T_0), R(R) {}; + + template<typename TType, typename RhoType> + auto alphaig(const TType& T, const RhoType& rho) const { + using otype = std::common_type_t <TType, RhoType>; + return forceeval(static_cast<otype>(c*(pow(T,t)*(1/(t+1)-1/t)-pow(T_0,t+1)/(T*(t+1))+pow(T_0,t)/t))); + } + }; + // The collection of possible terms that could be part of the summation using IdealHelmholtzTerms = std::variant < IdealHelmholtzConstant, IdealHelmholtzLead, - IdealHelmholtzLogT, + IdealHelmholtzLogT, + IdealHelmholtzPowerT, IdealHelmholtzPlanckEinstein, + IdealHelmholtzPlanckEinsteinGeneralized, IdealHelmholtzGERG2004Cosh, - IdealHelmholtzGERG2004Sinh + IdealHelmholtzGERG2004Sinh, + IdealHelmholtzCp0Constant, + IdealHelmholtzCp0PowerT > ; class PureIdealHelmholtz { @@ -138,37 +228,51 @@ namespace teqp { std::vector<IdealHelmholtzTerms> contributions; PureIdealHelmholtz(const nlohmann::json& jpure) { //std::string s = jpure.dump(1); - if (!jpure.is_array()) { - throw teqp::InvalidArgument("JSON data passed to PureIdealHelmholtz must be an array"); + if (jpure.is_array()) { + throw teqp::InvalidArgument("JSON data passed to PureIdealHelmholtz must be an object and contain the fields \"R\" and \"terms\""); } - for (auto& term : jpure) { + double R = jpure.at("R"); + for (auto& term : jpure.at("terms")) { if (!term.is_object()) { - throw teqp::InvalidArgument("JSON data for pure fluid must be an array"); + throw teqp::InvalidArgument("JSON data for pure fluid must be an object"); } //std::string s = term.dump(1); if (term.at("type") == "Constant") { // a - contributions.emplace_back(IdealHelmholtzConstant(term.at("a"))); + contributions.emplace_back(IdealHelmholtzConstant(term.at("a"), R)); } else if (term.at("type") == "Lead") { // ln(rho) + a_1 + a_2/T - contributions.emplace_back(IdealHelmholtzLead(term.at("a_1"), term.at("a_2"))); + contributions.emplace_back(IdealHelmholtzLead(term.at("a_1"), term.at("a_2"), R)); } else if (term.at("type") == "LogT") { // a*ln(T) - contributions.emplace_back(IdealHelmholtzLogT(term.at("a"))); + contributions.emplace_back(IdealHelmholtzLogT(term.at("a"), R)); + } + else if (term.at("type") == "PowerT") { // sum_i n_i * T^i + contributions.emplace_back(IdealHelmholtzPowerT(term.at("n"), term.at("t"), R)); } else if (term.at("type") == "PlanckEinstein") { - contributions.emplace_back(IdealHelmholtzPlanckEinstein(term.at("n"), term.at("theta"))); + contributions.emplace_back(IdealHelmholtzPlanckEinstein(term.at("n"), term.at("theta"), R)); + } + else if (term.at("type") == "PlanckEinsteinGeneralized") { + contributions.emplace_back(IdealHelmholtzPlanckEinsteinGeneralized(term.at("n"), term.at("c"), term.at("d"), term.at("theta"), R)); } else if (term.at("type") == "GERG2004Cosh") { - contributions.emplace_back(IdealHelmholtzGERG2004Cosh(term.at("n"), term.at("theta"))); + contributions.emplace_back(IdealHelmholtzGERG2004Cosh(term.at("n"), term.at("theta"), R)); } else if (term.at("type") == "GERG2004Sinh") { - contributions.emplace_back(IdealHelmholtzGERG2004Sinh(term.at("n"), term.at("theta"))); + contributions.emplace_back(IdealHelmholtzGERG2004Sinh(term.at("n"), term.at("theta"), R)); + } + else if (term.at("type") == "Cp0Constant") { + contributions.emplace_back(IdealHelmholtzCp0Constant(term.at("c"), term.at("T_0"), R)); + } + else if (term.at("type") == "Cp0PowerT") { + contributions.emplace_back(IdealHelmholtzCp0PowerT(term.at("c"), term.at("t"), term.at("T_0"), R)); } else { throw InvalidArgument("Don't understand this type: " + term.at("type").get<std::string>()); } } } + template<typename TType, typename RhoType> auto alphaig(const TType& T, const RhoType &rho) const{ std::common_type_t <TType, RhoType> ig = 0.0; @@ -185,7 +289,7 @@ namespace teqp { * * \f[ \alpha^{\rm ig} = \sum_i x_i[\alpha^{\rm ig}_{oi}(T,\rho) + x_i] \f] * - * where x_i are mole fractions + * where \f$x_i\f$ are mole fractions * */ class IdealHelmholtz { @@ -222,4 +326,172 @@ namespace teqp { return ig; } }; -} \ No newline at end of file + + inline nlohmann::json CoolProp2teqp_alphaig_term_reformatter(const nlohmann::json &term, double Tri, double rhori, double R){ + //std::string s = term.dump(1); + + if (term.at("type") == "IdealGasHelmholtzLead") { + // Was ln(delta) + a_1 + a_2*tau ==> ln(rho) + a_1 + a_2/T + // new a_1 is old a_1 - ln(rho_ri) + // new a_2 is old a_2 * Tri + return {{{"type", "Lead"}, {"a_1", term.at("a1").get<double>() - log(rhori)}, {"a_2", term.at("a2").get<double>() * Tri}, {"R", R}}}; + } + else if (term.at("type") == "IdealGasHelmholtzEnthalpyEntropyOffset") { + // Was a_1 + a_2*tau ==> a_1 + a_2/T + std::valarray<double> n = {term.at("a1").get<double>(), term.at("a2").get<double>()*Tri}; + std::valarray<double> t = {0, -1}; + return {{{"type", "PowerT"}, {"n", n}, {"t", t}, {"R", R}}}; + } + else if (term.at("type") == "IdealGasHelmholtzLogTau") { // a*ln(tau) + // Was a*ln(tau) = a*ln(Tri) - a*ln(T) ==> a*ln(T) + // Breaks into two pieces, first is a constant term a*ln(Tri), next is a*ln(T) piece + double a = term.at("a").get<double>(); + nlohmann::json term1 = {{"type", "Constant"}, {"a", a*log(Tri)}, {"R", R}}; + nlohmann::json term2 = {{"type", "LogT"}, {"a", -a}, {"R", R}}; + return {term1, term2}; + } + else if (term.at("type") == "IdealGasHelmholtzPlanckEinstein") { + // Was + std::valarray<double> n = term.at("n"); + std::valarray<double> theta = term.at("t").get<std::valarray<double>>()*Tri; + return {{{"type", "PlanckEinstein"}, {"n", n}, {"theta", theta}, {"R", R}}}; + } + else if (term.at("type") == "IdealGasHelmholtzPlanckEinsteinFunctionT") { + // Was + std::valarray<double> n = term.at("n"); + std::valarray<double> theta = term.at("v"); + return {{{"type", "PlanckEinstein"}, {"n", n}, {"theta", theta}, {"R", R}}}; + } + else if (term.at("type") == "IdealGasHelmholtzPlanckEinsteinGeneralized") { + // Was + std::valarray<double> n = term.at("n"); + std::valarray<double> c = term.at("c"); + std::valarray<double> d = term.at("d"); + std::valarray<double> theta = term.at("t").get<std::valarray<double>>()*Tri; +// std::cout << term.dump() << std::endl; + return {{{"type", "PlanckEinsteinGeneralized"}, {"n", n}, {"c", c}, {"d", d}, {"theta", theta}, {"R", R}}}; + } + else if (term.at("type") == "IdealGasHelmholtzPower") { + // Was + std::valarray<double> n = term.at("n").get<std::valarray<double>>(); + std::valarray<double> t = term.at("t").get<std::valarray<double>>(); + for (auto i = 0; i < n.size(); ++i){ + n[i] *= pow(Tri, t[i]); + t[i] *= -1; // T is in the denominator in CoolProp terms, in the numerator in teqp + } + return {{{"type", "PowerT"}, {"n", n}, {"t", t}, {"R", R}}}; + } + else if (term.at("type") == "IdealGasHelmholtzCP0PolyT") { + // Was + nlohmann::json newterms = nlohmann::json::array(); +// std::cout << term.dump() << std::endl; + std::valarray<double> t = term.at("t"), c = term.at("c"); + double T_0 = term.at("T0"); + for (auto i = 0; i < t.size(); ++i){ + if (t[i] == 0){ + newterms.push_back({{"type", "Cp0Constant"}, {"c", c[i]}, {"T_0", T_0}, {"R", R}}); + } + else{ + newterms.push_back({{"type", "Cp0PowerT"}, {"c", c[i]}, {"t", t[i]}, {"T_0", T_0}, {"R", R}}); + } + } + return newterms; + } + else if (term.at("type") == "IdealGasHelmholtzCP0Constant") { +// std::cout << term.dump() << std::endl; + double T_0 = term.at("T0"); + return {{{"type", "Cp0Constant"}, {"c", term.at("cp_over_R")}, {"T_0", T_0}, {"R", R}}}; + } + else if (term.at("type") == "IdealGasHelmholtzCP0AlyLee") { + // Was + nlohmann::json newterms = nlohmann::json::array(); +// std::cout << term.dump() << std::endl; + std::valarray<double> constants = term.at("c"); + double T_0 = term.at("T0"); + + // Take the constant term if nonzero + if (std::abs(constants[0]) > 1e-14) { + newterms.push_back({{"type", "Cp0Constant"}, {"c", constants[0]}, {"T_0", T_0}, {"R", R}}); + } + + std::vector<double> n, c, d, t; + if (std::abs(constants[1]) > 1e-14) { + // sinh term can be converted by setting a_k = C, b_k = 2*D, c_k = -1, d_k = 1 + n.push_back(constants[1]); + t.push_back(-2 * constants[2]); + c.push_back(1); + d.push_back(-1); + } + if (std::abs(constants[3]) > 1e-14) { + // cosh term can be converted by setting a_k = C, b_k = 2*D, c_k = 1, d_k = 1 + n.push_back(-constants[3]); + t.push_back(-2 * constants[4]); + c.push_back(1); + d.push_back(1); + } + newterms.push_back( + {{"type", "PlanckEinsteinGeneralized"}, {"n", n}, {"c", c}, {"d", d}, {"theta", t}, {"R", R}} + ); + return newterms; + } +// else if (term.at("type") == "GERG2004Cosh") { +// //contributions.emplace_back(IdealHelmholtzGERG2004Cosh(term.at("n"), term.at("theta"))); +// } +// else if (term.at("type") == "GERG2004Sinh") { +// //contributions.emplace_back(IdealHelmholtzGERG2004Sinh(term.at("n"), term.at("theta"))); +// } + else { + throw InvalidArgument("Don't understand this type of CoolProp ideal-gas Helmholtz energy term: " + term.at("type").get<std::string>()); + } + } + + /** + \brief Convert the ideal-gas term for a term from CoolProp-formatted JSON structure + + \param path A string, pointing to a filesystem file, or the JSON contents to be parsed + \returns j The JSON + + The key difference in the approach in CoolProp and teqp is that the contributions in teqp + are based on temperature and density as the independent variables, whereas the + implementation in CoolProp uses the pure fluid reciprocal reduced temperature and reduced + density as independent variables + */ + inline nlohmann::json convert_CoolProp_idealgas(const std::string &s, int index){ + + nlohmann::json j; + + // Get the JSON structure to be parsed + try{ + // First assume that the input argument is a path + std::filesystem::path p = s; + j = load_a_JSON_file(s); + } + catch(std::exception &){ + j = nlohmann::json::parse(s); + } + + // Extract the things from the CoolProp-formatted data structure + auto jEOS = j.at("EOS")[index]; + auto jCP = jEOS.at("alpha0"); + double Tri = jEOS.at("STATES").at("reducing").at("T"); + double rhori = jEOS.at("STATES").at("reducing").at("rhomolar"); + double R = jEOS.at("gas_constant"); + + // Now we transform the inputs to teqp-formatted terms + nlohmann::json newterms = nlohmann::json::array(); + for (const auto& term : jCP){ + // Converted can be a two-element array, so all terms are returned as array + // and then put into newterms + auto converted = CoolProp2teqp_alphaig_term_reformatter(term, Tri, rhori, R); + for (auto newterm : converted){ + newterms.push_back(newterm); + if (!newterms.back().is_object()){ + std::cout << newterm.dump() << std::endl; + } + } + } + + // And return our new data structure for this fluid + return {{"terms", newterms}, {"R", R}}; + } +} diff --git a/include/teqp/json_tools.hpp b/include/teqp/json_tools.hpp index e429482..2974fea 100644 --- a/include/teqp/json_tools.hpp +++ b/include/teqp/json_tools.hpp @@ -1,3 +1,4 @@ +#pragma once #include "nlohmann/json.hpp" #include <set> @@ -28,4 +29,4 @@ namespace teqp{ for (auto k : ks) { lengths.insert(j.at(k).size()); } return lengths.size() == 1; } -} \ No newline at end of file +} diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp index 992308a..6ca67ee 100644 --- a/interface/pybind11_wrapper.cpp +++ b/interface/pybind11_wrapper.cpp @@ -162,6 +162,7 @@ void init_teqp(py::module& m) { // The ideal gas Helmholtz energy class auto alphaig = py::class_<IdealHelmholtz>(m, "IdealHelmholtz").def(py::init<const nlohmann::json&>()); + alphaig.def_static("convert_CoolProp_format", [](const std::string &path, int index){return convert_CoolProp_idealgas(path, index);}); add_ig_derivatives<IdealHelmholtz>(m, alphaig); add_vdW(m); @@ -225,4 +226,4 @@ PYBIND11_MODULE(teqp, m) { m.doc() = "TEQP: Templated Equation of State Package"; m.attr("__version__") = TEQPVERSION; init_teqp(m); -} \ No newline at end of file +} diff --git a/src/tests/catch_test_multifluid.cxx b/src/tests/catch_test_multifluid.cxx index 2982a4f..95cf73a 100644 --- a/src/tests/catch_test_multifluid.cxx +++ b/src/tests/catch_test_multifluid.cxx @@ -1,5 +1,8 @@ #include <catch2/catch_test_macros.hpp> #include <catch2/catch_approx.hpp> +#include <catch2/generators/catch_generators.hpp> +#include <catch2/generators/catch_generators_adapters.hpp> +#include <catch2/generators/catch_generators_range.hpp> using Catch::Approx; @@ -8,6 +11,7 @@ using Catch::Approx; #include "teqp/algorithms/critical_tracing.hpp" #include "teqp/algorithms/VLE.hpp" #include "teqp/filesystem.hpp" +#include "teqp/ideal_eosterms.hpp" using namespace teqp; @@ -299,4 +303,23 @@ TEST_CASE("Calculate partial molar volume for a CO2 containing mixture", "[parti for (auto i = 0; i < expected.size(); ++i){ CHECK(expected[i] == Approx(der[i])); } -} \ No newline at end of file +} + +TEST_CASE("Check that all pure fluid ideal-gas terms can be converted", "[multifluid],[all],[ideal]") { + std::string root = "../mycp"; + auto paths = get_files_in_folder(root + "/dev/fluids", ".json"); + auto p = GENERATE_REF(from_range(paths)); + CHECK(std::filesystem::is_regular_file(p)); + CAPTURE(p); + // Check can be loaded from both path and string contents + auto jig = convert_CoolProp_idealgas(p.string(), 0 /* index of EOS */); + auto jig2 = convert_CoolProp_idealgas(load_a_JSON_file(p.string()).dump(), 0 /* index of EOS */); + // Convert to json array + nlohmann::json jaig = nlohmann::json::array(); jaig.push_back(jig); + CHECK(jaig.is_array()); + +// std::cout << jaig.dump() << std::endl; + + // Check that converted structures can be loaded + auto aig = IdealHelmholtz(jaig); +} -- GitLab