From d751a20162adb7342e90ad6462114fd5051f5894 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Sat, 21 May 2022 17:49:16 -0400
Subject: [PATCH] Enable setting kij for cubics
---
include/teqp/models/cubics.hpp | 32 +++++++++++++++++++++++---------
interface/cubics.cpp | 4 ++--
src/tests/catch_test_cubics.cxx | 20 ++++++++++++++++++++
3 files changed, 45 insertions(+), 11 deletions(-)
diff --git a/include/teqp/models/cubics.hpp b/include/teqp/models/cubics.hpp
index e8a4371..aa9f0c7 100644
--- a/include/teqp/models/cubics.hpp
+++ b/include/teqp/models/cubics.hpp
@@ -14,6 +14,8 @@ Implementations of the canonical cubic equations of state
#include "nlohmann/json.hpp"
+#include <Eigen/Dense>
+
namespace teqp {
/**
@@ -40,10 +42,10 @@ template <typename NumType, typename AlphaFunctions>
class GenericCubic {
protected:
std::valarray<NumType> ai, bi;
- std::valarray<std::valarray<NumType>> k;
const NumType Delta1, Delta2, OmegaA, OmegaB;
int superanc_index;
const AlphaFunctions alphas;
+ Eigen::ArrayXXd kmat;
nlohmann::json meta;
@@ -53,9 +55,21 @@ protected:
template<typename TType, typename IndexType>
auto get_bi(TType T, IndexType i) const { return bi[i]; }
+ void check_kmat(int N) {
+ if (kmat.cols() != kmat.rows()) {
+ throw teqp::InvalidArgument("kmat rows [" + std::to_string(kmat.rows()) + "] and columns [" + std::to_string(kmat.cols()) + "] are not identical");
+ }
+ if (kmat.cols() == 0) {
+ kmat.resize(N, N); kmat.setZero();
+ }
+ else if (kmat.cols() != N) {
+ throw teqp::InvalidArgument("kmat needs to be a square matrix the same size as the number of components [" + std::to_string(N) + "]");
+ }
+ };
+
public:
- GenericCubic(NumType Delta1, NumType Delta2, NumType OmegaA, NumType OmegaB, int superanc_index, const std::valarray<NumType>& Tc_K, const std::valarray<NumType>& pc_Pa, const AlphaFunctions& alphas)
- : Delta1(Delta1), Delta2(Delta2), OmegaA(OmegaA), OmegaB(OmegaB), superanc_index(superanc_index), alphas(alphas)
+ GenericCubic(NumType Delta1, NumType Delta2, NumType OmegaA, NumType OmegaB, int superanc_index, const std::valarray<NumType>& Tc_K, const std::valarray<NumType>& pc_Pa, const AlphaFunctions& alphas, const Eigen::ArrayXXd& kmat)
+ : Delta1(Delta1), Delta2(Delta2), OmegaA(OmegaA), OmegaB(OmegaB), superanc_index(superanc_index), alphas(alphas), kmat(kmat)
{
ai.resize(Tc_K.size());
bi.resize(Tc_K.size());
@@ -63,7 +77,7 @@ public:
ai[i] = OmegaA * pow2(Ru * Tc_K[i]) / pc_Pa[i];
bi[i] = OmegaB * Ru * Tc_K[i] / pc_Pa[i];
}
- k = std::valarray<std::valarray<NumType>>(std::valarray<NumType>(0.0, Tc_K.size()), Tc_K.size());
+ check_kmat(ai.size());
};
void set_meta(const nlohmann::json& j) { meta = j; }
@@ -101,7 +115,7 @@ public:
for (auto j = 0; j < molefracs.size(); ++j) {
auto alphaj = forceeval(std::visit([&](auto& t) { return t(T); }, alphas[j]));
auto aj_ = ai[j] * alphaj;
- auto aij = forceeval((1 - k[i][j]) * sqrt(ai_ * aj_));
+ auto aij = forceeval((1 - kmat(i,j)) * sqrt(ai_ * aj_));
a_ = a_ + molefracs[i] * molefracs[j] * aij;
}
}
@@ -134,7 +148,7 @@ public:
};
template <typename TCType, typename PCType, typename AcentricType>
-auto canonical_SRK(TCType Tc_K, PCType pc_K, AcentricType acentric) {
+auto canonical_SRK(TCType Tc_K, PCType pc_K, AcentricType acentric, const Eigen::ArrayXXd& kmat = {}) {
double Delta1 = 1;
double Delta2 = 0;
auto m = 0.48 + 1.574 * acentric - 0.176 * acentric * acentric;
@@ -156,13 +170,13 @@ auto canonical_SRK(TCType Tc_K, PCType pc_K, AcentricType acentric) {
{"kind", "Soave-Redlich-Kwong"}
};
- auto cub = GenericCubic(Delta1, Delta2, OmegaA, OmegaB, CubicSuperAncillary::SRK_CODE, Tc_K, pc_K, std::move(alphas));
+ auto cub = GenericCubic(Delta1, Delta2, OmegaA, OmegaB, CubicSuperAncillary::SRK_CODE, Tc_K, pc_K, std::move(alphas), kmat);
cub.set_meta(meta);
return cub;
}
template <typename TCType, typename PCType, typename AcentricType>
-auto canonical_PR(TCType Tc_K, PCType pc_K, AcentricType acentric) {
+auto canonical_PR(TCType Tc_K, PCType pc_K, AcentricType acentric, const Eigen::ArrayXXd& kmat = {}) {
double Delta1 = 1+sqrt(2);
double Delta2 = 1-sqrt(2);
AcentricType m = acentric*0.0;
@@ -189,7 +203,7 @@ auto canonical_PR(TCType Tc_K, PCType pc_K, AcentricType acentric) {
{"kind", "Peng-Robinson"}
};
- auto cub = GenericCubic(Delta1, Delta2, OmegaA, OmegaB, CubicSuperAncillary::PR_CODE, Tc_K, pc_K, std::move(alphas));
+ auto cub = GenericCubic(Delta1, Delta2, OmegaA, OmegaB, CubicSuperAncillary::PR_CODE, Tc_K, pc_K, std::move(alphas), kmat);
cub.set_meta(meta);
return cub;
}
diff --git a/interface/cubics.cpp b/interface/cubics.cpp
index 8f26cff..9598b80 100644
--- a/interface/cubics.cpp
+++ b/interface/cubics.cpp
@@ -6,8 +6,8 @@ void add_cubics(py::module& m) {
using va = std::valarray<double>;
- m.def("canonical_PR", &canonical_PR<va,va,va>, py::arg("Tc_K"), py::arg("pc_Pa"), py::arg("acentric"));
- m.def("canonical_SRK", &canonical_SRK<va, va, va>, py::arg("Tc_K"), py::arg("pc_Pa"), py::arg("acentric"));
+ m.def("canonical_PR", &canonical_PR<va,va,va>, py::arg("Tc_K"), py::arg("pc_Pa"), py::arg("acentric"), py::arg_v("kmat", Eigen::ArrayXXd(0, 0), "None"));
+ m.def("canonical_SRK", &canonical_SRK<va, va, va>, py::arg("Tc_K"), py::arg("pc_Pa"), py::arg("acentric"), py::arg_v("kmat", Eigen::ArrayXXd(0, 0), "None"));
using cub = decltype(canonical_PR(va{}, va{}, va{}));
auto wcub = py::class_<cub>(m, "GenericCubic")
diff --git a/src/tests/catch_test_cubics.cxx b/src/tests/catch_test_cubics.cxx
index 03800c5..ef98d14 100644
--- a/src/tests/catch_test_cubics.cxx
+++ b/src/tests/catch_test_cubics.cxx
@@ -35,6 +35,26 @@ TEST_CASE("Test construction of cubic", "[cubic]")
int rr = 0;
}
+TEST_CASE("Check SRK with kij setting", "[cubic]")
+{
+ // Values taken from http://dx.doi.org/10.6028/jres.121.011
+ std::valarray<double> Tc_K = { 190.564, 154.581, 150.687 },
+ pc_Pa = { 4599200, 5042800, 4863000 },
+ acentric = { 0.011, 0.022, -0.002 };
+ Eigen::ArrayXXd kij_right(3, 3); kij_right.setZero();
+ Eigen::ArrayXXd kij_bad(2, 20); kij_bad.setZero();
+
+ SECTION("No kij") {
+ CHECK_NOTHROW(canonical_SRK(Tc_K, pc_Pa, acentric));
+ }
+ SECTION("Correctly shaped kij matrix") {
+ CHECK_NOTHROW(canonical_SRK(Tc_K, pc_Pa, acentric, kij_right));
+ }
+ SECTION("Incorrectly shaped kij matrix") {
+ CHECK_THROWS(canonical_SRK(Tc_K, pc_Pa, acentric, kij_bad));
+ }
+}
+
TEST_CASE("Check calling superancillary curves", "[cubic][superanc]")
{
std::valarray<double> Tc_K = { 150.687 };
--
GitLab