Merge branch 'usnistgov:main' into main

include/teqp/models/multifluid.hpp

@@ -513,7 +513,7 @@ inline auto get_EOS_terms(const nlohmann::json& j)
     }
 
     // First check whether term type is allowed
-    const std::vector<std::string> allowed_types = { "ResidualHelmholtzPower", "ResidualHelmholtzGaussian", "ResidualHelmholtzNonAnalytic","ResidualHelmholtzGaoB", "ResidualHelmholtzLemmon2005", "ResidualHelmholtzExponential", "ResidualHelmholtzDoubleExponential","ResidualHelmholtzGenericCubic" };
+    const std::vector<std::string> allowed_types = { "ResidualHelmholtzPower", "ResidualHelmholtzGaussian", "ResidualHelmholtzNonAnalytic","ResidualHelmholtzGaoB", "ResidualHelmholtzLemmon2005", "ResidualHelmholtzExponential", "ResidualHelmholtzDoubleExponential","ResidualHelmholtzGenericCubic","ResidualHelmholtzPCSAFTGrossSadowski2001" };
 
     auto isallowed = [&](const auto& conventional_types, const std::string& name) {
         for (auto& a : conventional_types) { if (name == a) { return true; }; } return false;
@@ -742,6 +742,9 @@ inline auto get_EOS_terms(const nlohmann::json& j)
         else if (type == "ResidualHelmholtzGenericCubic") {
             container.add_term(GenericCubicTerm(term));
         }
+        else if (type == "ResidualHelmholtzPCSAFTGrossSadowski2001") {
+            container.add_term(PCSAFTGrossSadowski2001Term(term));
+        }
         else {
             throw std::invalid_argument("Bad term type: "+type);
         }

include/teqp/models/multifluid_eosterms.hpp

@@ -3,6 +3,7 @@
 #include "teqp/types.hpp"
 #include "teqp/exceptions.hpp"
 #include "teqp/models/cubics.hpp"
+#include "teqp/models/pcsaft.hpp"
 
 namespace teqp {
 
@@ -381,6 +382,28 @@ class GenericCubicTerm {
     }
 };
 
+/**
+ This implementation is for PC-SAFT for a pure fluid as taken from Gross & Sadowski, I&ECR, 2001
+ */
+class PCSAFTGrossSadowski2001Term {
+public:
+    const double Tred_K, rhored_molm3;
+    const PCSAFT::PCSAFTPureGrossSadowski2001 pcsaft;
+
+    PCSAFTGrossSadowski2001Term(const nlohmann::json& spec) :
+        Tred_K(spec.at("Tred / K")),
+        rhored_molm3(spec.at("rhored / mol/m^3")),
+        pcsaft(spec) // The remaining arguments will be consumed by the constructor
+    {}
+
+    template<typename TauType, typename DeltaType>
+    auto alphar(const TauType& tau, const DeltaType& delta) const {
+        auto T = forceeval(Tred_K/tau);
+        auto rhomolar = forceeval(delta*rhored_molm3);
+        return forceeval(pcsaft.alphar(T, rhomolar, Eigen::Array<double,1,1>{}));
+    }
+};
+
 template<typename... Args>
 class EOSTermContainer {  
 private:
@@ -406,7 +429,7 @@ class EOSTermContainer {
     }
 };
 
-using EOSTerms = EOSTermContainer<JustPowerEOSTerm, PowerEOSTerm, GaussianEOSTerm, NonAnalyticEOSTerm, Lemmon2005EOSTerm, GaoBEOSTerm, ExponentialEOSTerm, DoubleExponentialEOSTerm, GenericCubicTerm>;
+using EOSTerms = EOSTermContainer<JustPowerEOSTerm, PowerEOSTerm, GaussianEOSTerm, NonAnalyticEOSTerm, Lemmon2005EOSTerm, GaoBEOSTerm, ExponentialEOSTerm, DoubleExponentialEOSTerm, GenericCubicTerm, PCSAFTGrossSadowski2001Term>;
 
 using DepartureTerms = EOSTermContainer<JustPowerEOSTerm, PowerEOSTerm, GaussianEOSTerm, GERG2004EOSTerm, NullEOSTerm, DoubleExponentialEOSTerm,Chebyshev2DEOSTerm>;
 

include/teqp/models/pcsaft.hpp

@@ -516,5 +516,80 @@ inline auto PCSAFTfactory(const nlohmann::json& spec) {
     }
 }
 
+/**
+ The model of Gross & Sadowski, simplified down to the case of pure fluids
+ */
+class PCSAFTPureGrossSadowski2001{
+private:
+    Eigen::Array<double, 7, 1> aim, bim;
+public:
+    const double pi = 3.141592653589793238462643383279502884197;
+    const Eigen::Array<double, 7, 6> coeff;
+    const double m, sigma_A, eps_k;
+    double kappa1, kappa2;
+    PCSAFTPureGrossSadowski2001(const nlohmann::json&j) : coeff((Eigen::Array<double, 7, 6>() << 0.9105631445,-0.3084016918,-0.0906148351,0.7240946941,-0.5755498075,0.0976883116  ,
+                 0.6361281449,0.1860531159,0.4527842806,2.2382791861,0.6995095521,-0.2557574982    ,
+                 2.6861347891,-2.5030047259,0.5962700728,-4.0025849485,3.8925673390,-9.1558561530  ,
+                 -26.547362491,21.419793629,-1.7241829131,-21.003576815,-17.215471648,20.642075974 ,
+                 97.759208784,-65.255885330,-4.1302112531,26.855641363,192.67226447,-38.804430052  ,
+                 -159.59154087,83.318680481,13.776631870,206.55133841,-161.82646165,93.626774077   ,
+                 91.297774084,-33.746922930,-8.6728470368,-355.60235612,-165.20769346,-29.666905585).finished()),
+    m(j.at("m")), sigma_A(j.at("sigma / A")), eps_k(j.at("epsilon_over_k")) {
+        auto mfac1 = (m-1.0)/m;
+        auto mfac2 = (m-2.0)/m*mfac1;
+        aim = coeff.col(0) + coeff.col(1)*mfac1 + coeff.col(2)*mfac2;
+        bim = coeff.col(3) + coeff.col(4)*mfac1 + coeff.col(5)*mfac2;
+        kappa1 = (2.0*pi*eps_k*pow(m, 2)*pow(sigma_A, 3));
+        kappa2 = (pi*pow(eps_k, 2)*pow(m, 3)*pow(sigma_A, 3));
+    }
+
+    template<class VecType>
+    auto R(const VecType& molefrac) const {
+        return get_R_gas<decltype(molefrac[0])>();
+    }
+
+    template<typename TTYPE, typename RhoType, typename VecType>
+    auto alphar(const TTYPE& T, const RhoType& rhomolar, const VecType& /*mole_fractions*/) const {
+
+        auto rhoN_A3 = forceeval(rhomolar*N_A/1e30); // [A^3]
+
+        auto d = forceeval(sigma_A*(1.0-0.12*exp(-3.0*eps_k/T)));
+        Eigen::Array<decltype(d), 4, 1> dpowers; dpowers(0) = 1.0; for (auto i = 1U; i <= 3; ++i){ dpowers(i) = d*dpowers(i-1); }
+        auto zeta = pi/6.0*rhoN_A3*m*dpowers;
+
+        auto zeta2_to2 = zeta[2]*zeta[2];
+        auto zeta2_to3 = zeta2_to2*zeta[2];
+        auto zeta3_to2 = zeta[3]*zeta[3];
+        auto onemineta = forceeval(1.0-zeta[3]);
+        auto onemineta_to2 = onemineta*onemineta;
+        auto onemineta_to3 = onemineta*onemineta_to2;
+        auto onemineta_to4 = onemineta*onemineta_to3;
+
+        auto alpha_hs = (3.0*zeta[1]*zeta[2]/onemineta
+         + zeta2_to3/(zeta[3]*onemineta_to2)
+         + (zeta2_to3/zeta3_to2-zeta[0])*log(1.0-zeta[3]))/zeta[0];
+
+        auto fac_g_hs = d/2.0; // d*d/(2*d)
+        auto gii = (1.0/onemineta
+            + fac_g_hs*3.0*zeta[2]/onemineta_to2
+            + (fac_g_hs*fac_g_hs)*2.0*zeta2_to2/onemineta_to3);
+        auto alpha_hc = m*alpha_hs - (m-1)*log(gii);
+
+        auto eta = zeta[3];
+        auto eta2 = eta*eta;
+        auto eta3 = eta2*eta;
+        auto eta4 = eta2*eta2;
+        auto C1 = 1.0+m*(8.0*eta-2.0*eta2)/onemineta_to4+(1.0-m)*(20.0*eta-27.0*eta2+12.0*eta3-2.0*eta4)/onemineta_to2/((2.0-eta)*(2.0-eta));
+
+        Eigen::Array<decltype(eta), 7, 1> etapowers; etapowers(0) = 1.0; for (auto i = 1U; i <= 6; ++i){ etapowers(i) = eta*etapowers(i-1); }
+        auto I1 = (aim.array().template cast<decltype(eta)>()*etapowers).sum();
+        auto I2 = (bim.array().template cast<decltype(eta)>()*etapowers).sum();
+
+        auto alpha_disp = -kappa1*rhoN_A3*I1/T - kappa2*rhoN_A3*I2/C1/(T*T);
+
+        return forceeval(alpha_hc + alpha_disp);
+    }
+};
+
 } /* namespace PCSAFT */
 }; // namespace teqp

interface/CPP/model_factory_pcsaft.cpp

@@ -7,5 +7,8 @@ namespace teqp{
         std::unique_ptr<teqp::cppinterface::AbstractModel> make_PCSAFT(const nlohmann::json &spec){
             return teqp::cppinterface::adapter::make_owned(PCSAFT::PCSAFTfactory(spec));
         }
+        std::unique_ptr<teqp::cppinterface::AbstractModel> make_PCSAFTPureGrossSadowski2001(const nlohmann::json &spec){
+            return teqp::cppinterface::adapter::make_owned(PCSAFT::PCSAFTPureGrossSadowski2001(spec));
+        }
     }
 }

interface/CPP/teqp_impl_factory.cpp

@@ -13,6 +13,8 @@ namespace teqp {
 
         std::unique_ptr<teqp::cppinterface::AbstractModel> make_SAFTVRMie(const nlohmann::json &);
         std::unique_ptr<teqp::cppinterface::AbstractModel> make_PCSAFT(const nlohmann::json &);
+        std::unique_ptr<teqp::cppinterface::AbstractModel> make_PCSAFTPureGrossSadowski2001(const nlohmann::json &);
+
         std::unique_ptr<teqp::cppinterface::AbstractModel> make_GERG2004resid(const nlohmann::json &);
         std::unique_ptr<teqp::cppinterface::AbstractModel> make_GERG2008resid(const nlohmann::json &);
         std::unique_ptr<teqp::cppinterface::AbstractModel> make_GERG2004idealgas(const nlohmann::json &);
@@ -58,6 +60,7 @@ namespace teqp {
             {"SAFT-VR-Mie", [](const nlohmann::json& spec){ return make_SAFTVRMie(spec); }},
 
             {"PCSAFT", [](const nlohmann::json& spec){ return make_PCSAFT(spec); }},
+            {"PCSAFTPureGrossSadowski2001", [](const nlohmann::json& spec){ return make_PCSAFTPureGrossSadowski2001(spec); }},
 
             {"GERG2004resid", [](const nlohmann::json& spec){ return make_GERG2004resid(spec);}},
             {"GERG2008resid", [](const nlohmann::json& spec){ return make_GERG2008resid(spec);}},

src/tests/catch_test_PCSAFT.cxx

@@ -23,6 +23,32 @@ TEST_CASE("Single alphar check value", "[PCSAFT]")
     auto z = (Eigen::ArrayXd(1) << 1.0).finished();
     using tdx = teqp::TDXDerivatives<decltype(model), double>;
     CHECK(tdx::get_Ar00(model, T, Dmolar, z) == Approx(-0.032400020930842724));
+
+    nlohmann::json j = {
+        {"m", model.get_m()[0]},
+        {"sigma / A", model.get_sigma_Angstrom()[0]},
+        {"epsilon_over_k", model.get_epsilon_over_k_K()[0]}
+    };
+    PCSAFTPureGrossSadowski2001 pure(j);
+
+    auto valpure = teqp::TDXDerivatives<decltype(pure), double>::get_Ar00(pure, T, Dmolar, z);
+    CAPTURE(valpure);
+    CHECK(valpure == Approx(-0.032400020930842724));
+}
+
+TEST_CASE("Pure with neon", "[PCSAFT]")
+{
+    PCSAFTPureGrossSadowski2001 pure(R"({"m": 1.593, "sigma / A": 3.445, "epsilon_over_k": 176.47})"_json);
+    auto z = (Eigen::ArrayXd(1) << 1.0).finished();
+    auto valpure = teqp::TDXDerivatives<decltype(pure), double>::get_Ar00(pure, 450.0, 10000.0, z);
+    CAPTURE(valpure);
+    CHECK(valpure == Approx(-3.00381333e-01));
+    auto j = R"(
+    {"kind": "PCSAFTPureGrossSadowski2001", "model": {"m": 1.593, "sigma / A": 3.445, "epsilon_over_k": 176.47}}
+    )"_json;
+
+    auto model = make_model(j);
+    CHECK(model->get_Ar00(450.0, 10000.0, z) == Approx(-3.00381333e-01));
 }
 
 TEST_CASE("Check get_names and get_BibTeXKeys", "[PCSAFT]")

src/tests/catch_test_SAFTpolar.cxx

@@ -280,7 +280,7 @@ TEST_CASE("Benchmark CO2 with polar PC-SAFT model", "[CO2bench]"){
     };
 }
 
-TEST_CASE("Benchmark methane with PC-SAFT model", "[CO2bench]"){
+TEST_CASE("Benchmark methane with PC-SAFT model", "[methanebench]"){
     auto contents = R"(
     {
       "kind": "PCSAFT",
@@ -312,6 +312,30 @@ TEST_CASE("Benchmark methane with PC-SAFT model", "[CO2bench]"){
     };
 }
 
+TEST_CASE("Benchmark methane with pure PC-SAFT model", "[methanebench]"){
+    auto contents = R"(
+    {
+      "kind": "PCSAFTPureGrossSadowski2001",
+      "model": {"m": 1.593, "sigma / A": 3.445, "epsilon_over_k": 176.47}
+    }
+    )"_json;
+    auto model = teqp::cppinterface::make_model(contents);
+    auto z = (Eigen::ArrayXd(1) << 1.0).finished();
+
+    BENCHMARK("alphar"){
+        return model->get_Ar00(300, 10000, z);
+    };
+    BENCHMARK("Ar11"){
+        return model->get_Ar11(300, 10000, z);
+    };
+    BENCHMARK("Ar02"){
+        return model->get_Ar02(300, 10000, z);
+    };
+    BENCHMARK("Ar20"){
+        return model->get_Ar20(300, 10000, z);
+    };
+}
+
 
 TEST_CASE("Benchmark CO2+methane with polar PC-SAFT model", "[CO2bench]"){
     auto contents = R"(

src/tests/catch_test_multifluid.cxx

@@ -612,47 +612,82 @@ TEST_CASE("Check composition derivatives for ternary with all one component", "[
 }
 
 
-TEST_CASE("Check adding cubic EOS as pure fluid contribution in multifluid approach", "[multifluidpurecubic]") {
+TEST_CASE("Check adding cubic EOS as pure fluid contribution in multifluid approach", "[multifluidpuremodels]") {
     std::string root = FLUIDDATAPATH;
 
-    // Load some existing data from the JSON structure to avoid repeating ourselves
-    auto f = load_a_JSON_file(root+"/dev/fluids/CarbonDioxide.json");
-
-    // Flush existing contributions
-    f["EOS"][0]["alphar"].clear();
-
-    // Overwrite the residual portion with a cubic EOS
-    // In this case SRK which defines the values for
-    // OmegaA, OmegaB, Delta1 and Delta2, all other values taken
-    // from the Span&Wagner EOS
-    auto reducing = f["EOS"][0]["STATES"]["reducing"];
-    double Tc_K = reducing.at("T");
-    double pc_Pa = reducing.at("p");
-    double rhoc_molm3 = reducing.at("rhomolar");
-    f["EOS"][0]["alphar"].push_back({
-        {"R / J/mol/K", 8.31446261815324},
-        {"OmegaA", 0.42748023354034140439},
-        {"OmegaB", 0.086640349964957721589},
-        {"Delta1", 1.0},
-        {"Delta2", 0.0},
-        {"Tcrit / K", Tc_K},
-        {"pcrit / Pa", pc_Pa},
-        // Reducing state variables are taken from critical point
-        {"Tred / K", Tc_K},
-        {"rhored / mol/m^3", rhoc_molm3},
-        {"alpha", {
-            {{"type", "Twu"}, {"c", {1, 2, 3}}} // Dummy values for coefficients
-        }},
-        {"type", "ResidualHelmholtzGenericCubic"}
-    });
-//    std::cout << f["EOS"][0]["alphar"].dump(1) << std::endl;
-
-    nlohmann::json model = {
-        {"components", {f}}
-    };
-    nlohmann::json j = {
-        {"kind", "multifluid"},
-        {"model", model}
-    };
-    auto model_ = cppinterface::make_model(j);
+    SECTION("cubic EOS"){
+
+        // Load some existing data from the JSON structure to avoid repeating ourselves
+        auto f = load_a_JSON_file(root+"/dev/fluids/CarbonDioxide.json");
+
+        // Flush existing contributions
+        f["EOS"][0]["alphar"].clear();
+
+        // Overwrite the residual portion with a cubic EOS
+        // In this case SRK which defines the values for
+        // OmegaA, OmegaB, Delta1 and Delta2, all other values taken
+        // from the Span&Wagner EOS
+        auto reducing = f["EOS"][0]["STATES"]["reducing"];
+        double Tc_K = reducing.at("T");
+        double pc_Pa = reducing.at("p");
+        double rhoc_molm3 = reducing.at("rhomolar");
+        f["EOS"][0]["alphar"].push_back({
+            {"R / J/mol/K", 8.31446261815324},
+            {"OmegaA", 0.42748023354034140439},
+            {"OmegaB", 0.086640349964957721589},
+            {"Delta1", 1.0},
+            {"Delta2", 0.0},
+            {"Tcrit / K", Tc_K},
+            {"pcrit / Pa", pc_Pa},
+            // Reducing state variables are taken from critical point
+            {"Tred / K", Tc_K},
+            {"rhored / mol/m^3", rhoc_molm3},
+            {"alpha", {
+                {{"type", "Twu"}, {"c", {1, 2, 3}}} // Dummy values for coefficients
+            }},
+            {"type", "ResidualHelmholtzGenericCubic"}
+        });
+        //    std::cout << f["EOS"][0]["alphar"].dump(1) << std::endl;
+
+        nlohmann::json model = {
+            {"components", {f}}
+        };
+        nlohmann::json j = {
+            {"kind", "multifluid"},
+            {"model", model}
+        };
+        auto model_ = cppinterface::make_model(j);
+    }
+    SECTION("PC-SAFT"){
+        // Load some existing data from the JSON structure to avoid repeating ourselves
+        auto f = load_a_JSON_file(root+"/dev/fluids/CarbonDioxide.json");
+
+        // Flush existing contributions
+        f["EOS"][0]["alphar"].clear();
+
+        // Overwrite the residual portion with the PC-SAFT EOS
+        auto reducing = f["EOS"][0]["STATES"]["reducing"];
+        double Tc_K = reducing.at("T");
+        double pc_Pa = reducing.at("p");
+        double rhoc_molm3 = reducing.at("rhomolar");
+        f["EOS"][0]["alphar"].push_back({
+            // Reducing state variables are taken from critical point
+            {"Tred / K", Tc_K},
+            {"rhored / mol/m^3", rhoc_molm3},
+            {"m", 1.593}, // placeholder value for testing
+            {"sigma / A", 3.445}, // placeholder value for testing
+            {"epsilon_over_k", 176.47}, // placeholder value for testing
+            {"type", "ResidualHelmholtzPCSAFTGrossSadowski2001"}
+        });
+        //    std::cout << f["EOS"][0]["alphar"].dump(1) << std::endl;
+
+        nlohmann::json model = {
+            {"components", {f}}
+        };
+        nlohmann::json j = {
+            {"kind", "multifluid"},
+            {"model", model}
+        };
+        auto model_ = cppinterface::make_model(j);
+    }
 }