Skip to content

Commit

Permalink
Finalized python interfacing and example
Browse files Browse the repository at this point in the history
  • Loading branch information
ianhbell committed Oct 16, 2023
1 parent 08d751c commit 1dfaa27
Show file tree
Hide file tree
Showing 6 changed files with 301 additions and 30 deletions.
179 changes: 179 additions & 0 deletions doc/source/models/AdvancedCubicMixing.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,179 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "92b3c17a",
"metadata": {},
"source": [
"# Advanced cubic mixing rules\n",
"\n",
"\n",
"In the advanced cubic mixing rules, the conventional cubic EOS is taken as the basis for the method (usually Peng-Robinson), but different rules are used for the attractive term $a_m$. The formulation reads:\n",
"\n",
"$$\n",
"\\frac{a_m}{b_m} = \\sum_i z_i\\frac{a_i}{b_i} + \\frac{a^{E,\\gamma}_{\\rm res}}{CEoS}\n",
"$$\n",
"\n",
"where $CEoS$ is a scaling parameter that is in principle linked with the EOS coefficients, but can also be allowed to be an adjustable parameter. The $a_i$ and $b_i$ are the pure fluid values of component $i$. The $z_i$ are mole fractions. The mixture covolume is given by \n",
"$$\n",
"b_m = \\sum_i\\sum_j z_iz_jb_{ij}\n",
"$$\n",
"with \n",
"$$\n",
"b_{ij} = \\left(\\frac{b_i^{1/s} + b_j^{1/s}}{2}\\right)^s\n",
"$$\n",
"\n",
"The heart of the method is the definition of $a^{E,\\gamma}_{\\rm res}$, the residual contribution (not in the conventional thermodynamic sense) to the excess Helmholtz energy. There are many possible models here, but one that seems to work well is that of Wilson, for which the expression reads:\n",
"\n",
"$$\n",
"\\frac{a^{E,\\gamma}_{\\rm res}}{RT} = -\\sum_i z_i\\ln\\left(\\sum_j z_j \\Omega_{ji}(T)\\right) - \\sum_iz_i\\ln\\left(\\frac{\\phi_i}{z_i}\\right)\n",
"$$\n",
"with \n",
"$$\n",
"\\Omega_{ji}=\\frac{v_j}{v_i}\\exp(-A_{ij}/T)\n",
"$$\n",
"and \n",
"$$\n",
"\\frac{\\phi_i}{z_i} = \\frac{v_i}{\\sum_kz_kv_k}\n",
"$$\n",
"with the $v_i=b_i$. The parameter $A_{ij}\\neq A_{ji}$ in general, and is also given temperature dependence, which is also not supposed to be present according to the derivation. Thus, the models for $A_{ij}$ read something like this here:\n",
"$$\n",
"A_{ij} = m_{ij}T + n_{ij}\n",
"$$\n",
"so $m$ is non-dimensional and $n$ has units of temperature."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "de8cc3c3",
"metadata": {},
"outputs": [],
"source": [
"import numpy, matplotlib.pyplot as plt, numpy as np, pandas\n",
"import teqp\n",
"teqp.__version__"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7e7e2f96",
"metadata": {},
"outputs": [],
"source": [
"# Four isotherms of experimental data from doi: 10.1016/j.fluid.2016.05.015\n",
"import io, pandas\n",
"dat = pandas.read_csv(io.StringIO(\"\"\"PointID y1 uy1 x1 ux1 p/bar up T/K\n",
"1 0.0274 0.0007 0.0068 0.0002 59.830 0.053 293.10 \n",
"2 0.0664 0.0014 0.0183 0.0004 64.864 0.080 293.10 \n",
"3 0.0978 0.0020 0.0298 0.0007 69.772 0.080 293.10 \n",
"4 0.1199 0.0024 0.0424 0.0009 74.737 0.080 293.10 \n",
"5 0.1219 0.0028 0.1132 0.0023 89.869 0.080 293.10 \n",
"6 0.1339 0.0024 0.0995 0.0022 89.198 0.080 293.10 \n",
"7 0.1399 0.0026 0.0943 0.0020 88.853 0.080 293.10 \n",
"8 0.1461 0.0027 0.0823 0.0019 86.962 0.080 293.10 \n",
"9 0.1466 0.0028 0.0778 0.0017 85.942 0.080 293.10 \n",
"10 0.1466 0.0028 0.0772 0.0016 85.868 0.080 293.10 \n",
"1 0.1378 0.0027 0.0159 0.0004 42.667 0.051 273.08 \n",
"2 0.2143 0.0038 0.0297 0.0007 49.547 0.051 273.08 \n",
"3 0.2612 0.0043 0.0411 0.0009 55.238 0.051 273.08 \n",
"4 0.3209 0.0049 0.0609 0.0013 65.069 0.088 273.08 \n",
"5 0.3554 0.0051 0.0786 0.0016 73.395 0.088 273.08 \n",
"6 0.3758 0.0052 0.0978 0.0019 81.061 0.088 273.08 \n",
"7 0.3903 0.0053 0.1190 0.0023 90.706 0.088 273.08 \n",
"8 0.3914 0.0053 0.1477 0.0028 100.966 0.088 273.08 \n",
"9 0.3879 0.0053 0.1614 0.0030 104.806 0.088 273.08 \n",
"10 0.3724 0.0052 0.1875 0.0033 110.846 0.088 273.08\n",
"11 0.3550 0.0051 0.2068 0.0036 114.105 0.088 273.08\n",
"12 0.2727 0.0044 0.2531 0.0041 118.020 0.088 273.08\n",
"13 0.3343 0.0049 0.2268 0.0038 116.295 0.088 273.08\n",
"1 0.2048 0.0038 0.0106 0.0003 25.754 0.050 253.05 \n",
"2 0.3019 0.0049 0.0217 0.0005 30.479 0.050 253.05 \n",
"3 0.4638 0.0056 0.0436 0.0010 45.352 0.050 253.05 \n",
"4 0.5319 0.0056 0.0647 0.0014 58.188 0.050 253.05 \n",
"5 0.5854 0.0054 0.1077 0.0021 78.315 0.084 253.05 \n",
"6 0.5979 0.0054 0.1497 0.0028 98.276 0.084 253.05\n",
"7 0.5898 0.0054 0.1801 0.0032 109.241 0.084 253.05\n",
"8 0.5042 0.0057 0.0570 0.0012 51.343 0.084 253.05\n",
"9 0.5644 0.0055 0.0861 0.0017 67.594 0.084 253.05\n",
"10 0.5949 0.0054 0.1267 0.0024 86.883 0.084 253.05\n",
"11 0.5826 0.0054 0.2015 0.0035 116.614 0.084 253.05\n",
"12 0.5537 0.0055 0.2431 0.0040 129.873 0.084 253.05\n",
"13 0.4973 0.0055 0.2971 0.0046 139.161 0.084 253.05\n",
"14 0.4971 0.0055 0.2972 0.0046 139.261 0.084 253.05\n",
"1 0.7076 0.0050 0.0257 0.0006 27.983 0.056 223.10\n",
"2 0.7774 0.0041 0.0522 0.0011 44.918 0.056 223.10\n",
"3 0.8077 0.0036 0.0930 0.0019 64.906 0.081 223.10\n",
"4 0.8131 0.0035 0.1261 0.0024 84.799 0.081 223.10\n",
"5 0.8057 0.0035 0.1584 0.0029 104.410 0.081 223.10\n",
"6 0.7843 0.0038 0.1982 0.0035 125.782 0.081 223.10\n",
"7 0.7533 0.0041 0.2380 0.0040 144.287 0.081 223.10\n",
"8 0.7150 0.0045 0.2813 0.0044 159.015 0.081 223.10\n",
"9 0.6942 0.0047 0.3064 0.0047 165.347 0.081 223.10\n",
"\"\"\"), sep='\\s+', engine='python')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "49c01826",
"metadata": {},
"outputs": [],
"source": [
"# Model from Lasala, FPE, 2016: https://doi.org/10.1016/j.fluid.2016.05.015\n",
"j = {\n",
" \"kind\": \"advancedPRaEres\",\n",
" \"model\": {\n",
" \"Tcrit / K\": [304.21, 126.19],\n",
" \"pcrit / Pa\": [7.383e6, 3395800.0],\n",
" \"alphas\": [{\"type\": \"PR78\", \"acentric\": 0.22394}, {\"type\": \"PR78\", \"acentric\": 0.0372}],\n",
" \"aresmodel\": {\"type\": \"Wilson\", \"m\": [[0.0, -3.4768], [3.5332, 0.0]], \"n\": [[0.0, 825], [-585, 0.0]]},\n",
" \"options\": {\"s\": 2.0, \"brule\": \"Quadratic\", \"CEoS\": -0.52398}\n",
" }\n",
"}\n",
"\n",
"model = teqp.make_model(j)\n",
"for T in [223.15, 253.05, 273.08, 293.1]:\n",
" ipure = 0\n",
"\n",
" [rhoL0, rhoV0] = model.superanc_rhoLV(T, ipure)\n",
"\n",
" rhovecL0 = np.array([0.0, 0.0]); rhovecL0[ipure] = rhoL0\n",
" rhovecV0 = np.array([0.0, 0.0]); rhovecV0[ipure] = rhoV0\n",
"\n",
" J = model.trace_VLE_isotherm_binary(T, rhovecL0, rhovecV0, opt)\n",
" df = pandas.DataFrame(J)\n",
" plt.plot(df['xL_0 / mole frac.'], df['pL / Pa']/1e6,'k')\n",
" plt.plot(df['xV_0 / mole frac.'], df['pV / Pa']/1e6,'k')\n",
"\n",
"plt.plot(1-dat['x1'], dat['p/bar']/10, 'o')\n",
"plt.plot(1-dat['y1'], dat['p/bar']/10, '^')\n",
"\n",
"plt.gca().set(xlabel='$x_1$ / mole frac.', ylabel='$p$ / MPa', ylim=(0, 25))\n",
"plt.show()"
]
}
],
"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
}
1 change: 1 addition & 0 deletions doc/source/models/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ Models
make_model
cubics
QuantumPR
AdvancedCubicMixing
model_potentials
CPA
PCSAFT
Expand Down
126 changes: 96 additions & 30 deletions include/teqp/models/cubics.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,44 @@ class MathiasCopemanAlphaFunction {

using AlphaFunctionOptions = std::variant<BasicAlphaFunction<double>, TwuAlphaFunction<double>, MathiasCopemanAlphaFunction<double>>;

template<typename TC>
auto build_alpha_functions(const TC& Tc_K, const nlohmann::json& jalphas){
std::vector<AlphaFunctionOptions> alphas;
std::size_t i = 0;
if (jalphas.size() != Tc_K.size()){
throw teqp::InvalidArgument("alpha must be the same length as components");
}
for (auto alpha : jalphas){
std::string type = alpha.at("type");
if (type == "Twu"){
std::valarray<double> c_ = alpha.at("c");
Eigen::Array3d c = Eigen::Map<Eigen::Array3d>(&(c_[0]), c_.size());
alphas.emplace_back(TwuAlphaFunction(Tc_K[i], c));
}
else if (type == "Mathias-Copeman"){
std::valarray<double> c_ = alpha.at("c");
Eigen::Array3d c = Eigen::Map<Eigen::Array3d>(&(c_[0]), c_.size());
alphas.emplace_back(MathiasCopemanAlphaFunction(Tc_K[i], c));
}
else if (type == "PR78"){
double acentric = alpha.at("acentric");
double mi = 0.0;
if (acentric < 0.491) {
mi = 0.37464 + 1.54226*acentric - 0.26992*pow2(acentric);
}
else {
mi = 0.379642 + 1.48503*acentric -0.164423*pow2(acentric) + 0.016666*pow3(acentric);
}
alphas.emplace_back( BasicAlphaFunction(Tc_K[i], mi));
}
else{
throw teqp::InvalidArgument("alpha type is not understood: "+type);
}
i++;
}
return alphas;
};

template <typename NumType, typename AlphaFunctions>
class GenericCubic {
protected:
Expand Down Expand Up @@ -321,28 +359,6 @@ inline auto make_generalizedcubic(const nlohmann::json& spec){
double Delta1, Delta2, OmegaA, OmegaB;
std::string kind = "custom";

auto add_alphas = [&](const nlohmann::json& jalphas){
std::size_t i = 0;
if (jalphas.size() != Tc_K.size()){
throw teqp::InvalidArgument("alpha must be the same length as components");
}
for (auto alpha : jalphas){
std::string type = alpha.at("type");
std::valarray<double> c_ = alpha.at("c");
Eigen::Array3d c = Eigen::Map<Eigen::Array3d>(&(c_[0]), c_.size());
if (type == "Twu"){
alphas.emplace_back(TwuAlphaFunction(Tc_K[i], c));
}
else if (type == "Mathias-Copeman"){
alphas.emplace_back(MathiasCopemanAlphaFunction(Tc_K[i], c));
}
else{
throw teqp::InvalidArgument("alpha type is not understood: "+type);
}
i++;
}
};

if (spec.at("type") == "PR" ){
Delta1 = 1+sqrt(2.0);
Delta2 = 1-sqrt(2.0);
Expand All @@ -368,7 +384,7 @@ inline auto make_generalizedcubic(const nlohmann::json& spec){
if (!spec["alpha"].is_array()){
throw teqp::InvalidArgument("alpha must be array of objects");
}
add_alphas(spec.at("alpha"));
alphas = build_alpha_functions(Tc_K, spec.at("alpha"));
}
}
else if (spec.at("type") == "SRK"){
Expand All @@ -384,7 +400,7 @@ inline auto make_generalizedcubic(const nlohmann::json& spec){
if (!spec["alpha"].is_array()){
throw teqp::InvalidArgument("alpha must be array of objects");
}
add_alphas(spec.at("alpha"));
alphas = build_alpha_functions(Tc_K, spec.at("alpha"));
}
// See https://doi.org/10.1021/acs.iecr.1c00847
OmegaA = 1.0 / (9.0 * (cbrt(2) - 1));
Expand Down Expand Up @@ -416,12 +432,23 @@ inline auto make_generalizedcubic(const nlohmann::json& spec){

enum class AdvancedPRaEMixingRules {knotspecified, kLinear, kQuadratic};

NLOHMANN_JSON_SERIALIZE_ENUM( AdvancedPRaEMixingRules, {
{AdvancedPRaEMixingRules::knotspecified, nullptr},
{AdvancedPRaEMixingRules::kLinear, "Linear"},
{AdvancedPRaEMixingRules::kQuadratic, "Quadratic"},
})

struct AdvancedPRaEOptions{
AdvancedPRaEMixingRules brule = AdvancedPRaEMixingRules::kQuadratic;
double s = 2.0;
double CEoS = -sqrt(2.0)/2.0*log(1.0 + sqrt(2.0));
};

inline void from_json(const json& j, AdvancedPRaEOptions& o) {
j.at("brule").get_to(o.brule);
j.at("s").get_to(o.s);
j.at("CEoS").get_to(o.CEoS);
}

/**
A residual Helmholtz term that returns nothing (the empty term)
Expand Down Expand Up @@ -522,7 +549,7 @@ class WilsonResidualHelmholtzOverRT {
}
};

using ResidualHelmholtzOverRTOptions = std::variant<NullResidualHelmholtzOverRT<double>, WilsonResidualHelmholtzOverRT<double>>;
using ResidualHelmholtzOverRTVariant = std::variant<NullResidualHelmholtzOverRT<double>, WilsonResidualHelmholtzOverRT<double>>;

/**
Cubic EOS with advanced mixing rules, the EoS/aE method of Jaubert and co-workers
Expand All @@ -538,17 +565,21 @@ class AdvancedPRaEres {
const NumType OmegaA = 0.45723552892138218938;
const NumType OmegaB = 0.077796073903888455972;
const int superanc_code = CubicSuperAncillary::PR_CODE;
const double CEoS;


protected:

std::valarray<NumType> Tc_K, pc_Pa;

std::valarray<NumType> ai, bi;

const AlphaFunctions alphas;
const ResidualHelmholtzOverRTOptions ares;
const ResidualHelmholtzOverRTVariant ares;
Eigen::ArrayXXd lmat;
const double s;

const AdvancedPRaEMixingRules brule;
const double s;
const double CEoS;

nlohmann::json meta;

Expand All @@ -575,8 +606,8 @@ class AdvancedPRaEres {
};

public:
AdvancedPRaEres(const std::valarray<NumType>& Tc_K, const std::valarray<NumType>& pc_Pa, const AlphaFunctions& alphas, const ResidualHelmholtzOverRTOptions& ares, const Eigen::ArrayXXd& lmat, const AdvancedPRaEOptions& options = {})
: alphas(alphas), ares(ares), lmat(lmat), s(options.s), brule(options.brule), CEoS(options.CEoS)
AdvancedPRaEres(const std::valarray<NumType>& Tc_K, const std::valarray<NumType>& pc_Pa, const AlphaFunctions& alphas, const ResidualHelmholtzOverRTVariant& ares, const Eigen::ArrayXXd& lmat, const AdvancedPRaEOptions& options = {})
: Tc_K(Tc_K), pc_Pa(pc_Pa), alphas(alphas), ares(ares), lmat(lmat), brule(options.brule), s(options.s), CEoS(options.CEoS)
{
ai.resize(Tc_K.size());
bi.resize(Tc_K.size());
Expand All @@ -590,6 +621,9 @@ class AdvancedPRaEres {
void set_meta(const nlohmann::json& j) { meta = j; }
auto get_meta() const { return meta; }
auto get_lmat() const { return lmat; }
auto get_Tc_K() const { return Tc_K; }
auto get_pc_Pa() const { return pc_Pa; }

static double get_bi(double Tc_K, double pc_Pa){
const NumType OmegaB = 0.077796073903888455972;
const NumType R = 8.31446261815324;
Expand Down Expand Up @@ -691,6 +725,38 @@ class AdvancedPRaEres {
}
};

inline auto make_AdvancedPRaEres(const nlohmann::json& j){

std::valarray<double> Tc_K = j.at("Tcrit / K");
std::valarray<double> pc_Pa = j.at("pcrit / Pa");

std::vector<AlphaFunctionOptions> alphas = build_alpha_functions(Tc_K, j.at("alphas"));

auto get_ares_model = [&](const nlohmann::json& armodel) -> ResidualHelmholtzOverRTVariant {

std::string type = armodel.at("type");
if (type == "Wilson"){
std::vector<double> b;
for (auto i = 0; i < Tc_K.size(); ++i){
b.push_back(teqp::AdvancedPRaEres<double>::get_bi(Tc_K[i], pc_Pa[i]));
}
auto mWilson = build_square_matrix(armodel.at("m"));
auto nWilson = build_square_matrix(armodel.at("n"));
return WilsonResidualHelmholtzOverRT<double>(b, mWilson, nWilson);
}
else{
throw teqp::InvalidArgument("bad type of ares model: " + type);
}
};
auto aresmodel = get_ares_model(j.at("aresmodel"));

AdvancedPRaEOptions options = j.at("options");
auto model = teqp::AdvancedPRaEres<double>(Tc_K, pc_Pa, alphas, aresmodel, Eigen::ArrayXXd::Zero(2, 2), options);
return model;
}

using advancedPRaEres_t = decltype(make_AdvancedPRaEres({}));

/**
The quantum corrected Peng-Robinson model as developed in
Expand Down
Loading

0 comments on commit 1dfaa27

Please sign in to comment.