diff --git a/include/teqp/models/cubics.hpp b/include/teqp/models/cubics.hpp
index e8a437135568c73bab51bf1a6abc2c61061cb381..aa9f0c7d4342d5fee2d38addcd39bbd5cb77de21 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 8f26cffc7716e0dc5b10b9f3c808e1cf9cde95c1..9598b805b19f022015e4c59c5813c4ca6ab83713 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 03800c530c636ce4d81af3af3a2b0a67bcbcb21d..ef98d14bb65623d16a1ec002a726261b8fa1764a 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 };