Skip to content

Commit

Permalink
Merge pull request #90 from usnistgov/LKP
Browse files Browse the repository at this point in the history 
LKP model
  • Loading branch information
@ianhbell
ianhbell authored Feb 27, 2024
2 parents f0afee5 + 838054c commit 56eab31
Show file tree
Hide file tree
Showing 4 changed files with 231 additions and 0 deletions.
93 changes: 93 additions & 0 deletions include/teqp/models/LKP.hpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,93 @@
#pragma once // Only include header once in a compilation unit if it is included multiple times

#include <set>

#include "teqp/constants.hpp" // used for R
#include "teqp/types.hpp" // needed for forceeval
#include "teqp/exceptions.hpp" // to return teqp error messages

#include "nlohmann/json.hpp"

namespace teqp{
namespace LKP{

struct LKPFluidParameters {
std::vector<double> b, c, d;
double beta, gamma_, omega;
};

class LKPMix{
public:
const LKPFluidParameters simple{{0.0, 0.1181193, 0.265728, 0.154790, 0.303230e-1}, // b, with pad to match indices
{0.0, 0.236744e-1, 0.186984e-1, 0, 0.427240e-1}, // c, with pad to match indices
{0, 0.155428e-4, 0.623689e-4}, // d, with pad to match indices
0.653920, 0.601670e-1, 0.0}, // beta, gamma, omega
ref{{0, 0.2026579, 0.331511, 0.276550e-1, 0.203488}, // b, with pad to match indices
{0, 0.313385e-1, 0.503618e-1, 0.169010e-1,0.41577e-1}, // c, with pad to match indices
{0, 0.487360e-4, 0.740336e-5}, // d, with pad to match indices
1.226, 0.03754, 0.3978}; // beta, gamma, omega
const std::vector<double> Tcrit, pcrit, acentric;
const double m_R; ///< molar gas constant to be used in this model, in J/mol/K
const std::vector<std::vector<double>> kmat;

LKPMix(const std::vector<double>& Tcrit, const std::vector<double>& pcrit, const std::vector<double>& acentric, double R, const std::vector<std::vector<double>>& kmat) : Tcrit(Tcrit), pcrit(pcrit), acentric(acentric), m_R(R), kmat(kmat){
if (std::set<std::size_t>{Tcrit.size(), pcrit.size(), acentric.size()}.size() > 1){
throw teqp::InvalidArgument("The arrays should all be the same size.");
}
}

template<class VecType>
auto R(const VecType& /*molefrac*/) const { return m_R; }

/// Calculate the contribution for one of the fluids, depending on the parameter set passed in
template<typename TTYPE, typename RhoType, typename ZcType>
auto alphar_func(const TTYPE& tau, const RhoType& delta, const ZcType& Zc, const LKPFluidParameters& params) const {
auto B = params.b[1] - params.b[2]*tau - params.b[3]*powi(tau, 2) - params.b[4]*powi(tau, 3);
auto C = params.c[1] - params.c[2]*tau + params.c[3]*powi(tau, 3);
auto D = params.d[1] + params.d[2]*tau;
auto deltaZc = forceeval(delta/Zc);
return forceeval(B/Zc*delta + 1.0/2.0*C*powi(deltaZc, 2) + 1.0/5.0*D*powi(deltaZc, 5) - params.c[4]*powi(tau, 3)/(2*params.gamma_)*(params.gamma_*powi(deltaZc, 2)+params.beta+1.0)*exp(-params.gamma_*powi(deltaZc, 2)) + params.c[4]*powi(tau,3)/(2*params.gamma_)*(params.beta+1.0));
}

template<typename TTYPE, typename RhoType, typename VecType>
auto alphar(const TTYPE& T, const RhoType& rhomolar, const VecType& mole_fractions) const {

if (static_cast<std::size_t>(mole_fractions.size()) != acentric.size()){
throw teqp::InvalidArgument("The mole fractions should be of of size "+ std::to_string(acentric.size()));
}

const VecType& x = mole_fractions; // just an alias to save typing, no copy is invoked
std::decay_t<decltype(mole_fractions[0])> summer_omega = 0.0, summer_vcmix = 0.0, summer_Tcmix = 0.0;
double Ru = m_R;

for (auto i = 0; i < mole_fractions.size(); ++i){
summer_omega += mole_fractions[i]*acentric[i];
auto v_ci = (0.2905-0.085*acentric[i])*Ru*Tcrit[i]/pcrit[i];
for (auto j = 0; j < mole_fractions.size(); ++j){
auto v_cj = (0.2905-0.085*acentric[j])*Ru*Tcrit[j]/pcrit[j];
auto v_c_ij = 1.0/8.0*powi(cbrt(v_ci) + cbrt(v_cj), 3);
auto T_c_ij = kmat[i][j]*sqrt(Tcrit[i]*Tcrit[j]);
summer_vcmix += x[i]*x[j]*v_c_ij;
summer_Tcmix += x[i]*x[j]*pow(v_c_ij, 0.25)*T_c_ij;
}
}
auto omega_mix = summer_omega;
auto vc_mix = summer_vcmix;
auto Tc_mix = 1.0/pow(summer_vcmix, 0.25)*summer_Tcmix;
// auto pc_mix = (0.2905-0.085*omega_mix)*Ru*Tc_mix/vc_mix;
auto Zc = forceeval(0.2905-0.085*omega_mix);
auto tau = forceeval(Tc_mix/T);
auto delta = forceeval(vc_mix*rhomolar);

auto retval = (1.0-omega_mix/ref.omega)*alphar_func(tau, delta, Zc, simple) + (omega_mix/ref.omega)*alphar_func(tau, delta, Zc, ref);
return forceeval(retval);
}
};

// Factory function takes in JSON data and returns an instance of the LKPMix class
inline auto make_LKPMix(const nlohmann::json& j){
return LKPMix(j.at("Tcrit / K"), j.at("pcrit / Pa"), j.at("acentric"), j.at("R / J/mol/K"), j.at("kmat"));
}

}
}
1 change: 1 addition & 0 deletions include/teqp/models/fwd.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
#include "teqp/models/mie/lennardjones.hpp"
#include "teqp/models/mie/mie.hpp"
#include "teqp/models/GERG/GERG.hpp"
#include "teqp/models/LKP.hpp"

namespace teqp {

Expand Down
2 changes: 2 additions & 0 deletions interface/CPP/teqp_impl_factory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,8 @@ namespace teqp {
{"CPA", [](const nlohmann::json& spec){ return make_owned(CPA::CPAfactory(spec));}},
{"PCSAFT", [](const nlohmann::json& spec){ return make_owned(PCSAFT::PCSAFTfactory(spec));}},

{"LKP", [](const nlohmann::json& spec){ return make_owned(LKP::make_LKPMix(spec));}},

{"multifluid", [](const nlohmann::json& spec){ return make_owned(multifluidfactory(spec));}},
{"multifluid-ECS-HuberEly1994", [](const nlohmann::json& spec){ return make_owned(ECSHuberEly::ECSHuberEly1994(spec));}},
{"SW_EspindolaHeredia2009", [](const nlohmann::json& spec){ return make_owned(squarewell::EspindolaHeredia2009(spec.at("lambda")));}},
Expand Down
135 changes: 135 additions & 0 deletions src/tests/catch_test_LKP.cxx
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,135 @@
#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_approx.hpp>
using Catch::Approx;

#include "teqp/derivs.hpp"
#include "teqp/constants.hpp"
#include "teqp/models/LKP.hpp"
#include "teqp/cpp/teqpcpp.hpp"

TEST_CASE("Check LKP", "[LKP]"){
using namespace teqp::LKP;

SECTION("pure"){
// methane, check values from TREND
std::vector<double> Tc_K = {190.564};
std::vector<double> pc_Pa = {4.5992e6};
std::vector<double> acentric = {0.011};
std::vector<std::vector<double>> kmat{{1.0}};
double R = 8.3144598;
auto model = LKPMix(Tc_K, pc_Pa, acentric, R, kmat);
auto z = (Eigen::ArrayXd(1) << 1.0).finished();

// teqp::TDXDerivatives<decltype(model)>::get_Ar00(model, 300.0, 8000.0, z);

nlohmann::json modelspec{
{"Tcrit / K", Tc_K},
{"pcrit / Pa", pc_Pa},
{"acentric", acentric},
{"R / J/mol/K", R},
{"kmat", kmat}
};
// std::cout << spec.dump(2) << std::endl;

CHECK_NOTHROW(make_LKPMix(modelspec));
nlohmann::json badspec = modelspec;
badspec["kmat"] = 4.7;
CHECK_THROWS(make_LKPMix(badspec));

nlohmann::json spec{
{"kind", "LKP"},
{"model", modelspec}
};
CHECK_NOTHROW(teqp::cppinterface::make_model(spec));
auto ptr = teqp::cppinterface::make_model(spec);

struct Point{ double T, rhomolar, alphar_expected; };
for (auto& pt : std::vector<Point>{
{2.5555555555556E+02,1.7778555555556E+03,-1.0483673775469E-01},
{2.5555555555556E+02,2.6667333333333E+03,-1.5420416099894E-01},
{2.5555555555556E+02,3.5556111111111E+03,-2.0148143793072E-01},
{2.5555555555556E+02,4.4444888888889E+03,-2.4660690326310E-01},
{2.5555555555556E+02,5.3333666666667E+03,-2.8951201335592E-01},
{2.5555555555556E+02,6.2222444444444E+03,-3.3012829826090E-01},
{2.5555555555556E+02,7.1111222222222E+03,-3.6839324530079E-01},
{2.5555555555556E+02,8.0000000000000E+03,-4.0425381376141E-01},
{2.7777777777778E+02,1.0000000000000E-01,-5.0166505188877E-06},
{2.7777777777778E+02,8.8897777777778E+02,-4.3666195153684E-02},
{2.7777777777778E+02,1.7778555555556E+03,-8.5458261689186E-02},
{2.7777777777778E+02,2.6667333333333E+03,-1.2535703634585E-01},
{2.7777777777778E+02,3.5556111111111E+03,-1.6332285849984E-01},
{2.7777777777778E+02,4.4444888888889E+03,-1.9930364880452E-01},
{2.7777777777778E+02,5.3333666666667E+03,-2.3323937266306E-01},
{2.7777777777778E+02,6.2222444444444E+03,-2.6506680393110E-01},
{2.7777777777778E+02,7.1111222222222E+03,-2.9472345661764E-01},
{2.7777777777778E+02,8.0000000000000E+03,-3.2214967196215E-01},
{3.0000000000000E+02,1.0000000000000E-01,-4.1178433147697E-06},
{3.0000000000000E+02,8.8897777777778E+02,-3.5720133653505E-02},
{3.0000000000000E+02,1.7778555555556E+03,-6.9656790541858E-02},
{3.0000000000000E+02,2.6667333333333E+03,-1.0179417406048E-01},
{3.0000000000000E+02,3.5556111111111E+03,-1.3209892931558E-01},
{3.0000000000000E+02,4.4444888888889E+03,-1.6052611602765E-01},
{3.0000000000000E+02,5.3333666666667E+03,-1.8702217613376E-01},
{3.0000000000000E+02,6.2222444444444E+03,-2.1152813972213E-01},
{3.0000000000000E+02,7.1111222222222E+03,-2.3398217049700E-01},
{3.0000000000000E+02,8.0000000000000E+03,-2.5432064692151E-01}
}){
auto alphar_actual = teqp::TDXDerivatives<decltype(model)>::get_Ar00(model, pt.T, pt.rhomolar, z);
CHECK(alphar_actual == Approx(pt.alphar_expected).margin(1e-12));

auto alphar_actual_ptr = ptr->get_Ar00(pt.T, pt.rhomolar, z);
CHECK(alphar_actual_ptr == Approx(pt.alphar_expected).margin(1e-12));
}
}

SECTION("methane + nitrogen mix"){
// methane + nitrogen, check values from TREND
std::vector<double> Tc_K = {190.564, 126.192};
std::vector<double> pc_Pa = {4.5992e6, 3.3958e6};
std::vector<double> acentric = {0.011, 0.037};
std::vector<std::vector<double>> kmat{{1.0, 0.977},{0.977, 1.0}};
double R = 8.3144598;
auto model = LKPMix(Tc_K, pc_Pa, acentric, R, kmat);
auto z = (Eigen::ArrayXd(2) << 0.8, 0.2).finished();

auto zbad = (Eigen::ArrayXd(3) << 0.3, 0.3, 0.4).finished();
CHECK_THROWS(model.alphar(300.0, 8000.0, zbad));

struct Point{ double T, rhomolar, alphar_expected; };
for (auto& pt : std::vector<Point>{
{2.55555555555555E+02, 1.00000000000000E-01, -4.91536626760729E-06},
{2.55555555555555E+02, 8.88977777777777E+02, -4.28036717751216E-02},
{2.55555555555555E+02, 1.77785555555555E+03, -8.38097840422122E-02},
{2.55555555555555E+02, 2.66673333333333E+03, -1.23001289158271E-01},
{2.55555555555555E+02, 3.55561111111111E+03, -1.60342053977241E-01},
{2.55555555555555E+02, 4.44448888888888E+03, -1.95784417210887E-01},
{2.55555555555555E+02, 5.33336666666666E+03, -2.29273119192599E-01},
{2.55555555555555E+02, 6.22224444444444E+03, -2.60749579654642E-01},
{2.55555555555555E+02, 7.11112222222222E+03, -2.90155546674414E-01},
{2.55555555555555E+02, 8.00000000000000E+03, -3.17435227813972E-01},
{2.77777777777777E+02, 1.00000000000000E-01, -3.96079440378770E-06},
{2.77777777777777E+02, 8.88977777777777E+02, -3.43639521841312E-02},
{2.77777777777777E+02, 1.77785555555555E+03, -6.70251681046093E-02},
{2.77777777777777E+02, 2.66673333333333E+03, -9.79698096970028E-02},
{2.77777777777777E+02, 3.55561111111111E+03, -1.27167929631439E-01},
{2.77777777777777E+02, 4.44448888888888E+03, -1.54578932250068E-01},
{2.77777777777777E+02, 5.33336666666666E+03, -1.80154076480535E-01},
{2.77777777777777E+02, 6.22224444444444E+03, -2.03839249743354E-01},
{2.77777777777777E+02, 7.11112222222222E+03, -2.25577252196130E-01},
{2.77777777777777E+02, 8.00000000000000E+03, -2.45308899079892E-01},
{3.00000000000000E+02, 1.00000000000000E-01, -3.18316290879055E-06},
{3.00000000000000E+02, 8.88977777777777E+02, -2.74782561342433E-02},
{3.00000000000000E+02, 1.77785555555555E+03, -5.33101587343266E-02},
{3.00000000000000E+02, 2.66673333333333E+03, -7.74841431899544E-02},
{3.00000000000000E+02, 3.55561111111111E+03, -9.99748284135246E-02},
{3.00000000000000E+02, 4.44448888888888E+03, -1.20746896159584E-01},
{3.00000000000000E+02, 5.33336666666666E+03, -1.39756583212869E-01},
{3.00000000000000E+02, 6.22224444444444E+03, -1.56953388700025E-01},
{3.00000000000000E+02, 7.11112222222222E+03, -1.72281392690949E-01},
{3.00000000000000E+02, 8.00000000000000E+03, -1.85679636571464E-01}
}){
auto alphar_actual = teqp::TDXDerivatives<decltype(model)>::get_Ar00(model, pt.T, pt.rhomolar, z);
CHECK(alphar_actual == Approx(pt.alphar_expected).margin(1e-12));
}
}
}

0 comments on commit 56eab31

Please sign in to comment.