Skip to content

Commit

Permalink
Test the fixed CPA where radial distribution function is specified
Browse files Browse the repository at this point in the history 
Closes #81
  • Loading branch information
@ianhbell
ianhbell committed Dec 15, 2023
1 parent 921277e commit 869f066
Show file tree
Hide file tree
Showing 2 changed files with 46 additions and 7 deletions.
22 changes: 16 additions & 6 deletions include/teqp/models/CPA.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,15 @@ inline auto get_association_classes(const std::string& s) {

enum class radial_dist { CS, KG, OT };

inline auto get_radial_dist(const std::string& s) {
if (s == "CS") { return radial_dist::CS; }
else if (s == "KG") { return radial_dist::KG; }
else if (s == "OT") { return radial_dist::OT; }
else {
throw std::invalid_argument("bad association flag:" + s);
}
}

/// Function that calculates the association binding strength between site A of molecule i and site B on molecule j
template<typename BType, typename TType, typename RhoType, typename VecType>
inline auto get_DeltaAB_pure(radial_dist dist, double epsABi, double betaABi, BType b_cubic, TType RT, RhoType rhomolar, const VecType& molefrac) {
Expand Down Expand Up @@ -69,7 +78,7 @@ inline auto get_DeltaAB_pure(radial_dist dist, double epsABi, double betaABi, BT
///

template<typename BType, typename TType, typename RhoType, typename VecType>
inline auto XA_calc_pure(int N_sites, association_classes scheme, double epsABi, double betaABi, const BType b_cubic, const TType RT, const RhoType rhomolar, const VecType& molefrac) {
inline auto XA_calc_pure(int N_sites, association_classes scheme, radial_dist dist, double epsABi, double betaABi, const BType b_cubic, const TType RT, const RhoType rhomolar, const VecType& molefrac) {

// Matrix XA(A, j) that contains all of the fractions of sites A not bonded to other active sites for each molecule i
// Start values for the iteration(set all sites to non - bonded, = 1)
Expand All @@ -79,7 +88,6 @@ inline auto XA_calc_pure(int N_sites, association_classes scheme, double epsABi,
XA.setOnes();

// Get the association strength between the associating sites
auto dist = radial_dist::KG; // TODO: pass this in
auto DeltaAiBj = get_DeltaAB_pure(dist, epsABi, betaABi, b_cubic, RT, rhomolar, molefrac);

if (scheme == association_classes::a1A) { // Acids
Expand Down Expand Up @@ -185,6 +193,7 @@ class CPAAssociation {
private:
const Cubic cubic;
const std::vector<association_classes> classes;
const radial_dist dist;
const std::valarray<double> epsABi, betaABi;
const std::vector<int> N_sites;
const double R_gas;
Expand All @@ -207,8 +216,8 @@ class CPAAssociation {
}

public:
CPAAssociation(const Cubic &&cubic, const std::vector<association_classes>& classes, const std::valarray<double> &epsABi, const std::valarray<double> &betaABi, double R_gas)
: cubic(cubic), classes(classes), epsABi(epsABi), betaABi(betaABi), N_sites(get_N_sites(classes)), R_gas(R_gas) {};
CPAAssociation(const Cubic &&cubic, const std::vector<association_classes>& classes, const radial_dist dist, const std::valarray<double> &epsABi, const std::valarray<double> &betaABi, double R_gas)
: cubic(cubic), classes(classes), dist(dist), epsABi(epsABi), betaABi(betaABi), N_sites(get_N_sites(classes)), R_gas(R_gas) {};

template<typename TType, typename RhoType, typename VecType>
auto alphar(const TType& T, const RhoType& rhomolar, const VecType& molefrac) const {
Expand All @@ -217,7 +226,7 @@ class CPAAssociation {

// Calculate the fraction of sites not bonded with other active sites
auto RT = forceeval(R_gas * T); // R times T
auto XA = XA_calc_pure(N_sites[0], classes[0], epsABi[0], betaABi[0], b_cubic, RT, rhomolar, molefrac);
auto XA = XA_calc_pure(N_sites[0], classes[0], dist, epsABi[0], betaABi[0], b_cubic, RT, rhomolar, molefrac);

using return_type = std::common_type_t<decltype(T), decltype(rhomolar), decltype(molefrac[0])>;
return_type alpha_r_asso = 0.0;
Expand Down Expand Up @@ -280,6 +289,7 @@ inline auto CPAfactory(const nlohmann::json &j){
auto build_assoc = [](const auto &&cubic, const auto& j) {
auto N = j["pures"].size();
std::vector<association_classes> classes;
radial_dist dist = get_radial_dist(j.at("radial_dist"));
std::valarray<double> epsABi(N), betaABi(N);
std::size_t i = 0;
for (auto p : j["pures"]) {
Expand All @@ -288,7 +298,7 @@ inline auto CPAfactory(const nlohmann::json &j){
classes.push_back(get_association_classes(p["class"]));
i++;
}
return CPAAssociation(std::move(cubic), classes, epsABi, betaABi, j["R_gas / J/mol/K"]);
return CPAAssociation(std::move(cubic), classes, dist, epsABi, betaABi, j["R_gas / J/mol/K"]);
};
return CPAEOS(build_cubic(j), build_assoc(build_cubic(j), j));
}
Expand Down
31 changes: 30 additions & 1 deletion src/tests/catch_tests.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -424,6 +424,35 @@ TEST_CASE("Test pure VLE with non-unity R0/Rr", "") {
auto rr = 0;
}

TEST_CASE("Test water Clapeyron.jl", "[CPA]") {
std::valarray<double> a0 = {0.12277}, bi = {0.0000145}, c1 = {0.6736}, Tc = {647.13},
molefrac = {1.0};
auto R = 8.31446261815324;
CPA::CPACubic cub(CPA::cubic_flag::SRK, a0, bi, c1, Tc, R);
double T = 303.15, rhomolar = 1/1.7915123921401366e-5;

auto z = (Eigen::ArrayXd(1) << 1).finished();

using tdx = TDXDerivatives<decltype(cub)>;
auto alphar = cub.alphar(T, rhomolar, molefrac);
CHECK(alphar == Approx(-1.2713135319123854)); // Value from Clapeyron.jl
double p_noassoc = T*rhomolar*R*(1+tdx::get_Ar01(cub, T, rhomolar, z));
CAPTURE(p_noassoc);

std::vector<CPA::association_classes> schemes = { CPA::association_classes::a4C };
std::valarray<double> epsAB = { 16655 }, betaAB = { 0.0692 };
CPA::radial_dist dist = CPA::radial_dist::CS;
CPA::CPAAssociation cpaa(std::move(cub), schemes, dist, epsAB, betaAB, R);

CPA::CPAEOS cpa(std::move(cub), std::move(cpaa));
using tdc = TDXDerivatives<decltype(cpa)>;
auto Ar01 = tdc::get_Ar01(cpa, T, rhomolar, z);
double p_withassoc = T*rhomolar*R*(1 + Ar01);
CAPTURE(p_withassoc);

REQUIRE(p_withassoc == Approx(100000.000));
}

TEST_CASE("Test water", "[CPA]") {
std::valarray<double> a0 = {0.12277}, bi = {0.000014515}, c1 = {0.67359}, Tc = {647.096},
molefrac = {1.0};
Expand All @@ -441,7 +470,7 @@ TEST_CASE("Test water", "[CPA]") {

std::vector<CPA::association_classes> schemes = { CPA::association_classes::a4C };
std::valarray<double> epsAB = { 16655 }, betaAB = { 0.0692 };
CPA::CPAAssociation cpaa(std::move(cub), schemes, epsAB, betaAB, R);
CPA::CPAAssociation cpaa(std::move(cub), schemes, CPA::radial_dist::KG, epsAB, betaAB, R);

CPA::CPAEOS cpa(std::move(cub), std::move(cpaa));
using tdc = TDXDerivatives<decltype(cpa)>;
Expand Down

0 comments on commit 869f066

Please sign in to comment.