From 278abc13225409d1d140e7900cc854eebb58e155 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Thu, 19 May 2022 15:33:03 -0400
Subject: [PATCH] Add ability to set kij for PC-SAFT
---
include/teqp/models/pcsaft.hpp | 26 ++++++++++++++++++++------
interface/PCSAFT.cpp | 2 +-
src/tests/catch_test_PCSAFT.cxx | 25 +++++++++++++++++++++++++
3 files changed, 46 insertions(+), 7 deletions(-)
create mode 100644 src/tests/catch_test_PCSAFT.cxx
diff --git a/include/teqp/models/pcsaft.hpp b/include/teqp/models/pcsaft.hpp
index 9fd6d3f..d1d5677 100644
--- a/include/teqp/models/pcsaft.hpp
+++ b/include/teqp/models/pcsaft.hpp
@@ -1,6 +1,7 @@
#pragma once
#include "nlohmann/json.hpp"
+#include "teqp/exceptions.hpp"
namespace teqp {
namespace PCSAFT {
@@ -177,10 +178,23 @@ private:
sigma_Angstrom, ///<
epsilon_over_k; ///< depth of pair potential divided by Boltzman constant
std::vector<std::string> names;
- double k_ij; ///< binary interaction parameter
+ Eigen::ArrayXXd kmat; ///< binary interaction parameter matrix
+
+ void check_kmat(int N) {
+ if (kmat.cols() != kmat.rows()) {
+ throw teqp::InvalidArgument("kmat rows and columns 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");
+ }
+ };
public:
- PCSAFTMixture(const std::vector<std::string> &names) : names(names)
+ PCSAFTMixture(const std::vector<std::string> &names, const Eigen::ArrayXXd& kmat = {}) : names(names), kmat(kmat)
{
+ check_kmat(names.size());
m.resize(names.size());
mminus1.resize(names.size());
sigma_Angstrom.resize(names.size());
@@ -195,10 +209,11 @@ public:
epsilon_over_k[i] = coeff.epsilon_over_k;
i++;
}
- k_ij = 0;
};
- PCSAFTMixture(const std::vector<SAFTCoeffs> &coeffs)
+ PCSAFTMixture(const std::vector<SAFTCoeffs> &coeffs, const Eigen::ArrayXXd &kmat = {}) : kmat(kmat)
{
+ check_kmat(coeffs.size());
+
m.resize(coeffs.size());
mminus1.resize(coeffs.size());
sigma_Angstrom.resize(coeffs.size());
@@ -213,7 +228,6 @@ public:
names[i] = coeff.name;
i++;
}
- k_ij = 0;
};
auto get_m(){ return m; }
auto get_sigma_Angstrom() { return sigma_Angstrom; }
@@ -259,7 +273,7 @@ public:
for (std::size_t j = 0; j < N; ++j) {
// Eq. A.5
auto sigma_ij = 0.5 * sigma_Angstrom[i] + 0.5 * sigma_Angstrom[j];
- auto eij_over_k = sqrt(epsilon_over_k[i] * epsilon_over_k[j]) * (1.0 - k_ij);
+ auto eij_over_k = sqrt(epsilon_over_k[i] * epsilon_over_k[j]) * (1.0 - kmat(i,j));
c.m2_epsilon_sigma3_bar = c.m2_epsilon_sigma3_bar + mole_fractions[i] * mole_fractions[j] * m[i] * m[j] * eij_over_k / T * pow(sigma_ij, 3);
c.m2_epsilon2_sigma3_bar = c.m2_epsilon2_sigma3_bar + mole_fractions[i] * mole_fractions[j] * m[i] * m[j] * pow(eij_over_k / T, 2) * pow(sigma_ij, 3);
}
diff --git a/interface/PCSAFT.cpp b/interface/PCSAFT.cpp
index 578b830..5796f57 100644
--- a/interface/PCSAFT.cpp
+++ b/interface/PCSAFT.cpp
@@ -17,7 +17,7 @@ void add_PCSAFT(py::module& m) {
;
auto wPCSAFT = py::class_<PCSAFTMixture>(m, "PCSAFTEOS")
- .def(py::init<const std::vector<std::string>&>(), py::arg("names"))
+ .def(py::init<const std::vector<std::string>&, const Eigen::ArrayXXd&>(), py::arg("names"), py::arg_v("kmat", Eigen::ArrayXXd(2,2).setZero(), "None"))
.def(py::init<const std::vector<SAFTCoeffs>&>(), py::arg("coeffs"))
.def("print_info", &PCSAFTMixture::print_info)
.def("max_rhoN", &PCSAFTMixture::max_rhoN<Eigen::ArrayXd>)
diff --git a/src/tests/catch_test_PCSAFT.cxx b/src/tests/catch_test_PCSAFT.cxx
new file mode 100644
index 0000000..8ee5144
--- /dev/null
+++ b/src/tests/catch_test_PCSAFT.cxx
@@ -0,0 +1,25 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/catch_approx.hpp>
+
+using Catch::Approx;
+
+#include "teqp/core.hpp"
+#include "teqp/models/pcsaft.hpp"
+using namespace teqp::PCSAFT;
+
+TEST_CASE("Check PCSAFT with kij", "[PCSAFT]")
+{
+ std::vector<std::string> names = { "Methane", "Ethane" };
+ Eigen::ArrayXXd kij_right(2, 2); kij_right.setZero();
+ Eigen::ArrayXXd kij_bad(2, 20); kij_bad.setZero();
+
+ SECTION("No kij") {
+ CHECK_NOTHROW(PCSAFTMixture(names));
+ }
+ SECTION("Correctly shaped kij matrix") {
+ CHECK_NOTHROW(PCSAFTMixture(names, kij_right));
+ }
+ SECTION("Incorrectly shaped kij matrix") {
+ CHECK_THROWS(PCSAFTMixture(names, kij_bad));
+ }
+}
\ No newline at end of file
--
GitLab