Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

RK-PR #76

Merged
merged 2 commits into from
Nov 7, 2023
Merged

RK-PR #76

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5890129
Size check for advanced cubic
ianhbell Nov 7, 2023
8319cbd
Implement the RK-PR method of Cismondi
ianhbell Nov 7, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
202 changes: 202 additions & 0 deletions doc/source/models/RK-PR.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,202 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "e86a8de0",
"metadata": {},
"source": [
"# RK-PR\n",
"\n",
"The EOS can be given as \n",
"\n",
"$$ \\alpha^{\\rm r} = \\psi^{(-)} - \\dfrac{a_m}{RT } \\psi^{(+)} $$\n",
"\n",
"$$ \\psi^{(-)} =-\\ln(1-b_m\\rho ) $$\n",
"\n",
"$$ \\psi^{(+)} = \\dfrac{\\ln\\left(\\dfrac{\\Delta_1 b_m\\rho+1}{\\Delta_2b_m\\rho+1}\\right)}{b_m(\\Delta_1-\\Delta_2)} $$\n",
"\n",
"with the EOS fixed constants of\n",
"$$\n",
"\\Delta_1 = \\sum_i x_i \\delta_{1,i}\n",
"$$\n",
"$$\n",
"\\Delta_2 = \\frac{1-\\Delta_1}{1+\\Delta_1}\n",
"$$\n",
"\n",
"The attractive term goes like\n",
"$$\n",
"a_{i} = a_{c,i}\\left(\\frac{2}{3+T/T_{c,i}}\\right)^{k_i}\n",
"$$\n",
"with quadratic mixing rules\n",
"$$\n",
"a_m = \\sum_i\\sum_jx_ix_j(1-k_{ij})\\sqrt{a_{i}(T)a_{j}(T)}\n",
"$$\n",
"And the covolume also gets quadratic mixing rules\n",
"$$\n",
"b_m = \\sum_i\\sum_jx_ix_j(1-l_{ij})(b_{i}+b_{j})/2\n",
"$$\n",
"\n",
"Thus, to implement the RK-PR model in predictive mode, the following steps are required:\n",
"\n",
"1. Obtain the critical parameters Tc, pc\n",
"2. Solve for delta_1 from the experimental critical compressibility factor, begin with the values from the correlation\n",
"3. Solve for k by fixing the pressure at the T=0.7*Tc. In the case (e.g, CO$_2$) that Tt < 0.7*Tc, use instead Tr=Tt/Tc \n",
"\n",
"It may be necessary to adjust the values of $\\delta_{1,i}$ and $k_i$ for an individual component to better match the behavior of more polar components."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e1aa4062",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np \n",
"import scipy.optimize \n",
"import matplotlib.pyplot as plt \n",
"import teqp, numpy as np\n",
"import CoolProp.CoolProp as CP \n",
"import pandas\n",
"\n",
"def delta1_correlation(Zc):\n",
" # Eq. B.4 of Cismondi FPE 2005\n",
" d1 = 0.428363\n",
" d2 = 18.496215\n",
" d3 = 0.338426\n",
" d4 = 0.660000\n",
" d5 = 789.723105\n",
" d6 = 2.512392\n",
" return d1 + d2*(d3-Zc)**d4 + d5*(d3-Zc)**d6 \n",
"\n",
"def Zc_delta1(delta1):\n",
" # Eqs. B.1 to B.3 of Cismondi FPE 2005\n",
" d1 = (1+delta1**2)/(1+delta1)\n",
" y = 1 + (2*(1+delta1))**(1/3) + (4/(1+delta1))**(1/3)\n",
" return y/(3*y + d1 - 1)\n",
"\n",
"DELTA1 = np.linspace(np.sqrt(2)-1, 25, 1000)\n",
"ZZ = Zc_delta1(DELTA1)\n",
"plt.plot(DELTA1, ZZ, label='values')\n",
"DELTA1back = delta1_correlation(ZZ)\n",
"plt.axvline(np.sqrt(2)-1)\n",
"plt.plot(DELTA1back, ZZ, dashes=[2,2], label='correlation')\n",
"plt.gca().set(ylabel='$Z_c$', xlabel='$\\delta_1$')\n",
"plt.legend(loc='best')\n",
"plt.show()\n",
"\n",
"# for Zc in np.linspace(0.2, 0.3383, 1000):\n",
"# resid = lambda x: Zc_delta1(x)-Zc\n",
"# # print(resid(delta1_correlation(Zc)))\n",
"# print(Zc, scipy.optimize.newton(resid, delta1_correlation(Zc)), delta1_correlation(Zc))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8fb5d777",
"metadata": {},
"outputs": [],
"source": [
"names = ['CO2', 'n-Decane']\n",
"\n",
"R = 8.31446261815324\n",
"Tc = np.array([CP.PropsSI(k,\"Tcrit\") for k in names])\n",
"pc = np.array([CP.PropsSI(k,\"pcrit\") for k in names])\n",
"rhoc = np.array([CP.PropsSI(k,\"rhomolar_critical\") for k in names])\n",
"Zcexp = pc/(rhoc*R*Tc)\n",
"\n",
"# Use a rescaled Zc to obtain delta_1\n",
"Zc = 1.168*Zcexp\n",
"\n",
"delta_1 = [scipy.optimize.newton(lambda x: Zc_delta1(x)-Zc_, delta1_correlation(Zc_)) for Zc_ in Zc]\n",
"\n",
"def solve_for_k(i, p_target, Tr):\n",
" \"\"\" \n",
" The value of k for the i-th component is based on getting \n",
" the right vapor pressure, so a rootfinding routing is\n",
" used to obtain these values\n",
" \"\"\"\n",
" def objective(k):\n",
" j = {\n",
" \"kind\": \"RKPRCismondi2005\",\n",
" \"model\": {\n",
" \"delta_1\": [delta_1[i]],\n",
" \"Tcrit / K\": [Tc[i]],\n",
" \"pcrit / Pa\": [pc[i]],\n",
" \"k\": [k],\n",
" \"kmat\": [[0.0]],\n",
" \"lmat\": [[0.0]],\n",
" }\n",
" }\n",
" model = teqp.make_model(j)\n",
" T = Tr*Tc[i]\n",
" z = np.array([1.0])\n",
" a, b = model.get_ab(T, z)\n",
"\n",
" anc = teqp.build_ancillaries(model, Tc[i], rhoc[i], 150)\n",
" rhoL, rhoV = model.pure_VLE_T(T, anc.rhoL(T), anc.rhoV(T), 10)\n",
" p = T*R*rhoL*(1+model.get_Ar01(T, rhoL, z))\n",
" \n",
" return p-p_target\n",
" return scipy.optimize.newton(objective, 2.1)\n",
"\n",
"Tr = 0.7\n",
"i = 1\n",
"k_C10 = solve_for_k(i, CP.PropsSI('P','T',Tr*Tc[i],'Q',0,names[i]), Tr)\n",
"\n",
"model = teqp.make_model({\n",
" \"kind\": \"RKPRCismondi2005\",\n",
" \"model\": {\n",
" \"delta_1\": delta_1,\n",
" \"Tcrit / K\": Tc.tolist(),\n",
" \"pcrit / Pa\": pc.tolist(),\n",
" \"k\": [2.23854, k_C10],\n",
" \"kmat\": [[0,0],[0,0]],\n",
" \"lmat\": [[0,0],[0,0]],\n",
" }\n",
"})\n",
"\n",
"# Start at both pures\n",
"for ipure in [0, 1]:\n",
" Tc, rhoc = model.solve_pure_critical(300, 5000, {\"alternative_pure_index\":ipure, \"alternative_length\": 2})\n",
" z = np.array([0.0, 0.0]); z[ipure] = 1.0\n",
" pc = Tc*R*rhoc*(1+model.get_Ar01(Tc, rhoc, z))\n",
" plt.plot(Tc, pc/1e5, 'o')\n",
" \n",
" opt = teqp.TCABOptions(); opt.polish=True; opt.verbosity=100; opt.integration_order=5; opt.rel_err=1e-10; opt.abs_err=1e-10\n",
" trace = model.trace_critical_arclength_binary(Tc, z*rhoc, options=opt)\n",
" df = pandas.DataFrame(trace)\n",
" plt.plot(df['T / K'], df['p / Pa']/1e5)\n",
" \n",
"# Overlay the data from Reamer and Sage, Cismondi additional data points not present in Reamer and Sage\n",
"Tc_K = [310.928, 344.261, 377.594, 410.928, 444.261, 477.594, 510.928]\n",
"pc_kPa = np.array([7997.92, 12824.25, 16492.26, 18560.69, 18836.48, 17836.74, 15333.94])\n",
"plt.plot(Tc_K, pc_kPa/1e2, 'o')\n",
"\n",
"plt.gca().set(xlabel='$T$ / K', ylabel='$p$ / bar');"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
119 changes: 119 additions & 0 deletions include/teqp/models/cubics.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -849,6 +849,9 @@ class QuantumCorrectedPR{

template<typename TType, typename RhoType, typename FractionsType>
auto alphar(const TType& T, const RhoType& rhoinit, const FractionsType& molefrac) const {
if (molefrac.size() != alphas.size()) {
throw std::invalid_argument("Sizes do not match");
}
// First shift the volume by the volume translation
auto c = get_c(T, molefrac);
auto rho = 1.0/(1.0/rhoinit + c);
Expand Down Expand Up @@ -878,4 +881,120 @@ class QuantumCorrectedPR{
}
};


/**
The RK-PR method as defined in:

Martín Cismondi, Jørgen Mollerup,
Development and application of a three-parameter RK–PR equation of state,
Fluid Phase Equilibria,
Volume 232, Issues 1–2,
2005,
Pages 74-89,
https://doi.org/10.1016/j.fluid.2005.03.020.
*/
class RKPRCismondi2005{
public:
const double Ru = get_R_gas<double>(); /// Universal gas constant, exact number

private:
const std::vector<double> delta_1, Tc_K, pc_Pa, k;
const Eigen::ArrayXXd kmat, lmat;
const std::vector<double> a_c, b_c;

/// A convenience function to save some typing
std::vector<double> get_(const nlohmann::json &j, const std::string& k) const { return j.at(k).get<std::vector<double>>(); }

/// Calculate the parameters \f$y\f$ and \f$d_1\f$
auto get_yd1(double delta_1){
return std::make_tuple(1 + cbrt(2*(1+delta_1)) + cbrt(4/(1+delta_1)), (1+delta_1*delta_1)/(1+delta_1));
}

auto build_ac(){
std::vector<double> a_c(delta_1.size());
for (auto i = 0; i < delta_1.size(); ++i){
auto [y, d_1] = get_yd1(delta_1[i]);
auto Omega_a = (3*y*y + 3*y*d_1 + d_1*d_1 + d_1 - 1.0)/pow(3.0*y + d_1 - 1.0, 2);
a_c[i] = Omega_a*pow(Ru*Tc_K[i], 2)/pc_Pa[i];
}
return a_c;
}
auto build_bc(){
std::vector<double> b_c(delta_1.size());
for (auto i = 0; i < delta_1.size(); ++i){
auto [y, d_1] = get_yd1(delta_1[i]);
auto Omega_b = 1.0/(3.0*y + d_1 - 1.0);
b_c[i] = Omega_b*Ru*Tc_K[i]/pc_Pa[i];
}
return b_c;
}
public:

RKPRCismondi2005(const nlohmann::json &j) : delta_1(get_(j, "delta_1")), Tc_K(get_(j, "Tcrit / K")), pc_Pa(get_(j, "pcrit / Pa")), k(get_(j, "k")), kmat(build_square_matrix(j.at("kmat"))), lmat(build_square_matrix(j.at("lmat"))), a_c(build_ac()), b_c(build_bc()) {}

template<class VecType>
auto R(const VecType& /*molefrac*/) const {
return Ru;
}
auto get_delta_1() const { return delta_1; }
auto get_Tc_K() const { return Tc_K; }
auto get_pc_Pa() const { return pc_Pa; }
auto get_k() const { return k; }
auto get_kmat() const { return kmat; }
auto get_lmat() const { return lmat; }

template<typename TType>
auto get_bi(std::size_t i, const TType& T) const {
return b_c[i];
}

template<typename TType>
auto get_ai(std::size_t i, const TType& T) const {
return forceeval(a_c[i]*pow(3.0/(2.0+T/Tc_K[i]), k[i]));
}

template<typename TType, typename FractionsType>
auto get_ab(const TType& T, const FractionsType& z) const{
using numtype = std::common_type_t<TType, decltype(z[0])>;
numtype b = 0.0;
numtype a = 0.0;
std::size_t N = delta_1.size();
for (auto i = 0; i < N; ++i){
auto bi = get_bi(i, T);
auto ai = get_ai(i, T);
for (auto j = 0; j < N; ++j){
auto aj = get_ai(j, T);
auto bj = get_bi(j, T);

a += z[i]*z[j]*sqrt(ai*aj)*(1.0 - kmat(i,j));
b += z[i]*z[j]*(bi + bj)/2.0*(1.0 - lmat(i,j));

}
}
return std::make_tuple(a, b);
}

template<typename TType, typename RhoType, typename FractionsType>
auto alphar(const TType& T, const RhoType& rho, const FractionsType& molefrac) const {
if (molefrac.size() != delta_1.size()) {
throw std::invalid_argument("Sizes do not match");
}

// The mixture Delta_1 is a mole-fraction-weighted average of the pure values of delta_1
const auto delta_1_view = Eigen::Map<const Eigen::ArrayXd>(&(delta_1[0]), delta_1.size());
const decltype(molefrac[0]) Delta1 = (molefrac*delta_1_view).eval().sum();

auto Delta2 = (1.0-Delta1)/(1.0+Delta1);
auto [a, b] = get_ab(T, molefrac);

auto Psiminus = -log(1.0 - b*rho);
auto Psiplus = log((Delta1 * b * rho + 1.0) / (Delta2 * b * rho + 1.0)) / (b * (Delta1 - Delta2));
auto val = Psiminus - a/(Ru*T)*Psiplus;
return forceeval(val);
}
};
using RKPRCismondi2005_t = decltype(RKPRCismondi2005({}));



}; // namespace teqp
1 change: 1 addition & 0 deletions interface/CPP/teqp_impl_factory.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ namespace teqp {
{"cubic", [](const nlohmann::json& spec){ return make_owned(make_generalizedcubic(spec));}},
{"QCPRAasen", [](const nlohmann::json& spec){ return make_owned(QuantumCorrectedPR(spec));}},
{"advancedPRaEres", [](const nlohmann::json& spec){ return make_owned(make_AdvancedPRaEres(spec));}},
{"RKPRCismondi2005", [](const nlohmann::json& spec){ return make_owned(RKPRCismondi2005(spec));}},

{"CPA", [](const nlohmann::json& spec){ return make_owned(CPA::CPAfactory(spec));}},
{"PCSAFT", [](const nlohmann::json& spec){ return make_owned(PCSAFT::PCSAFTfactory(spec));}},
Expand Down
13 changes: 13 additions & 0 deletions interface/pybind11_wrapper.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,7 @@ const std::type_index EXP6_Kataoka1992_i{std::type_index(typeid(EXP6_Kataoka1992
const std::type_index twocenterLJF_i{std::type_index(typeid(twocenterLJF_t))};
const std::type_index QuantumPR_i{std::type_index(typeid(QuantumPR_t))};
const std::type_index advancedPRaEres_i{std::type_index(typeid(advancedPRaEres_t))};
const std::type_index RKPRCismondi2005_i{std::type_index(typeid(RKPRCismondi2005_t))};

/**
At runtime we can add additional model-specific methods that only apply for a particular model. We take in a Python-wrapped
Expand Down Expand Up @@ -186,6 +187,18 @@ void attach_model_specific_methods(py::object& obj){
setattr("get_pc_Pa", MethodType(py::cpp_function([](py::object& o){ return get_typed<advancedPRaEres_t>(o).get_pc_Pa(); }), obj));
// setattr("get_meta", MethodType(py::cpp_function([](py::object& o){ return get_typed<advancedPRaEres_t>(o).get_meta(); }), obj));
// setattr("set_meta", MethodType(py::cpp_function([](py::object& o, const std::string& s){ return get_mutable_typed<advancedPRaEres_t>(o).set_meta(s); }, "self"_a, "s"_a), obj));
}
else if (index == RKPRCismondi2005_i){
setattr("get_ab", MethodType(py::cpp_function([](py::object& o, double T, REArrayd& molefrac){ return get_typed<RKPRCismondi2005_t>(o).get_ab(T, molefrac); }, "self"_a, "T"_a, "molefrac"_a), obj));
setattr("get_delta_1", MethodType(py::cpp_function([](py::object& o){ return get_typed<RKPRCismondi2005_t>(o).get_delta_1(); }), obj));
setattr("get_Tc_K", MethodType(py::cpp_function([](py::object& o){ return get_typed<RKPRCismondi2005_t>(o).get_Tc_K(); }), obj));
setattr("get_pc_Pa", MethodType(py::cpp_function([](py::object& o){ return get_typed<RKPRCismondi2005_t>(o).get_pc_Pa(); }), obj));
setattr("get_k", MethodType(py::cpp_function([](py::object& o){ return get_typed<RKPRCismondi2005_t>(o).get_k(); }), obj));

setattr("get_kmat", MethodType(py::cpp_function([](py::object& o){ return get_typed<RKPRCismondi2005_t>(o).get_kmat(); }), obj));
setattr("get_lmat", MethodType(py::cpp_function([](py::object& o){ return get_typed<RKPRCismondi2005_t>(o).get_lmat(); }), obj));
// setattr("get_meta", MethodType(py::cpp_function([](py::object& o){ return get_typed<RKPRCismondi2005_t>(o).get_meta(); }), obj));
// setattr("set_meta", MethodType(py::cpp_function([](py::object& o, const std::string& s){ return get_mutable_typed<RKPRCismondi2005_t>(o).set_meta(s); }, "self"_a, "s"_a), obj));
}
else if (index == AmmoniaWaterTillnerRoth_i){
setattr("TcNH3", get_typed<AmmoniaWaterTillnerRoth>(obj).TcNH3);
Expand Down
Loading