diff --git a/doc/notebooks/Sympy derivatives of alphar from VTPR.ipynb b/doc/notebooks/Sympy derivatives of alphar from VTPR.ipynb new file mode 100644 index 00000000..1f7eb70f --- /dev/null +++ b/doc/notebooks/Sympy derivatives of alphar from VTPR.ipynb @@ -0,0 +1,69 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "from __future__ import division\n", + "from sympy import *\n", + "init_session(quiet=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "T,delta,rho_r,b_m,c_m,a_m,R_u = symbols('T,delta,rho_r,b_m,c_m,a_m,R_u')\n", + "W = symbols('W', cls=Function)(delta)\n", + "\n", + "alphar = -log(1-delta*rho_r*(b_m-c_m)) - sqrt(2)*a_m/(4*R_u*T*b_m)*log(W);\n", + "display(alphar)\n", + "for ndelta in range(1,5):\n", + " ss = simplify(diff(alphar, delta, ndelta))\n", + " display(ss)\n", + "\n", + "W =(1+delta*rho_r*(b_m*(1+sqrt(2)+c_m))) / (1+delta*rho_r*(b_m*(1-sqrt(2)+c_m)))\n", + "for ndelta in range(1,5):\n", + " display(diff(W,delta,ndelta))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.11" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/include/AbstractState.h b/include/AbstractState.h index d7b16c56..b3f90ff5 100644 --- a/include/AbstractState.h +++ b/include/AbstractState.h @@ -503,6 +503,8 @@ public: virtual std::string get_binary_interaction_string(const std::string &CAS1, const std::string &CAS2, const std::string ¶meter){ throw NotImplementedError("get_binary_interaction_string is not implemented for this backend"); }; /// Apply a simple mixing rule (EXPERT USE ONLY!!!) virtual void apply_simple_mixing_rule(std::size_t i, std::size_t j, const std::string &model) { throw NotImplementedError("apply_simple_mixing_rule is not implemented for this backend"); }; + // Set fluid parameter (currently the volume translation parameter for cubic) + virtual void set_fluid_parameter_double(const size_t i, const std::string parameter, const double value) { throw ValueError("set_fluid_parameter_double only defined for cubic backends"); }; /// Clear all the cached values virtual bool clear(); diff --git a/include/Configuration.h b/include/Configuration.h index bf867549..24bcc1ac 100644 --- a/include/Configuration.h +++ b/include/Configuration.h @@ -33,6 +33,7 @@ X(HENRYS_LAW_TO_GENERATE_VLE_GUESSES, "HENRYS_LAW_TO_GENERATE_VLE_GUESSES", false, "If true, when doing water-based mixture dewpoint calculations, use Henry's Law to generate guesses for liquid-phase composition") \ X(PHASE_ENVELOPE_STARTING_PRESSURE_PA, "PHASE_ENVELOPE_STARTING_PRESSURE_PA", 100.0, "Starting pressure [Pa] for phase envelope construction") \ X(R_U_CODATA, "R_U_CODATA", 8.3144598, "The value for the ideal gas constant in J/mol/K according to CODATA 2014. This value is used to harmonize all the ideal gas constants. This is especially important in the critical region.") \ + X(VTPR_UNIFAQ_PATH, "VTPR_UNIFAQ_PATH", "", "The path to the directory containing the UNIFAQ JSON files. Should be slash terminated") \ // Use preprocessor to create the Enum enum configuration_keys{ diff --git a/include/CoolPropLib.h b/include/CoolPropLib.h index 3ebef291..10b8a899 100644 --- a/include/CoolPropLib.h +++ b/include/CoolPropLib.h @@ -412,6 +412,19 @@ */ EXPORT_CODE void CONVENTION AbstractState_set_binary_interaction_double(const long handle, const long i, const long j, const char* parameter, const double value, long *errcode, char *message_buffer, const long buffer_length); + /** + * @brief Set some fluid parameter (ie volume translation for cubic) + * @param handle The integer handle for the state class stored in memory + * @param i indice of the fluid the parramter should be applied too (for mixtures) + * @param parameter the string specifying the parameter to use, ex "cm" for volume translation + * @param value the value of the parameter + * @param errcode The errorcode that is returned (0 = no error, !0 = error) + * @param message_buffer A buffer for the error code + * @param buffer_length The length of the buffer for the error code + * @return + */ + EXPORT_CODE void CONVENTION AbstractState_set_fluid_parameter_double(const long handle, const size_t i, const char* parameter, const double value, long *errcode, char *message_buffer, const long buffer_length); + // ************************************************************************************* // ************************************************************************************* // ***************************** DEPRECATED ******************************************* diff --git a/include/DataStructures.h b/include/DataStructures.h index 79599737..e6abfb08 100644 --- a/include/DataStructures.h +++ b/include/DataStructures.h @@ -406,7 +406,8 @@ enum backend_families { TTSE_BACKEND_FAMILY, BICUBIC_BACKEND_FAMILY, SRK_BACKEND_FAMILY, - PR_BACKEND_FAMILY + PR_BACKEND_FAMILY, + VTPR_BACKEND_FAMILY }; enum backends { INVALID_BACKEND = 0, @@ -420,7 +421,8 @@ enum backends { TTSE_BACKEND, BICUBIC_BACKEND, SRK_BACKEND, - PR_BACKEND + PR_BACKEND, + VTPR_BACKEND }; /// Convert a string into the enum values diff --git a/src/AbstractState.cpp b/src/AbstractState.cpp index 0669f50c..81b3787d 100644 --- a/src/AbstractState.cpp +++ b/src/AbstractState.cpp @@ -19,6 +19,7 @@ #include "Backends/Helmholtz/Fluids/FluidLibrary.h" #include "Backends/IF97/IF97Backend.h" #include "Backends/Cubics/CubicBackend.h" +#include "Backends/Cubics/VTPRBackend.h" #if !defined(NO_TABULAR_BACKENDS) #include "Backends/Tabular/TTSEBackend.h" @@ -89,6 +90,10 @@ AbstractState * AbstractState::factory(const std::string &backend, const std::ve { return new PengRobinsonBackend(fluid_names, get_config_double(R_U_CODATA)); } + else if (f1==VTPR_BACKEND_FAMILY) + { + return new VTPRBackend(fluid_names, get_config_double(R_U_CODATA)); + } else if (!backend.compare("?") || backend.empty()) { std::size_t idel = fluid_names[0].find("::"); diff --git a/src/Backends/Cubics/CubicBackend.cpp b/src/Backends/Cubics/CubicBackend.cpp index a3cfb59a..6ade26f1 100644 --- a/src/Backends/Cubics/CubicBackend.cpp +++ b/src/Backends/Cubics/CubicBackend.cpp @@ -210,43 +210,25 @@ void CoolProp::AbstractCubicBackend::update(CoolProp::input_pairs input_pair, do void CoolProp::AbstractCubicBackend::rho_Tp_cubic(CoolPropDbl T, CoolPropDbl p, int &Nsolns, double &rho0, double &rho1, double &rho2){ AbstractCubic *cubic = get_cubic().get(); double R = cubic->get_R_u(); - double Delta_1 = cubic->get_Delta_1(); - double Delta_2 = cubic->get_Delta_2(); double am = cubic->am_term(cubic->T_r/T, mole_fractions_double, 0); double bm = cubic->bm_term(mole_fractions); + double cm = cubic->cm_term(); + + // Introducing new variables to simplify the equation: + double d1 = cm - bm; + double d2 = cm + cubic->get_Delta_1()*bm; + double d3 = cm + cubic->get_Delta_2()*bm; - // ** SYMPY CODE *** - // R,T,v,b,a,Delta_1,Delta_2,p,Z,rho = symbols('R,T,v,b,a,Delta_1,Delta_2,p,Z,rho') - // eqn = (R*T/(v-b)-a/(v+Delta_1*b)/(v+Delta_2*b)-p).subs(v,1/rho) - // eqn2 = eqn*(-b + 1/rho)*(Delta_1*b + 1/rho)*(Delta_2*b + 1/rho) - // display(simplify(factor(expand(eqn2),rho))) - // - // yields: - // (b*rho**3*(Delta_1*Delta_2*R*T*b + Delta_1*Delta_2*b**2*p + a) - p + rho**2*(-Delta_1*Delta_2*b**2*p + Delta_1*R*T*b + Delta_1*b**2*p + Delta_2*R*T*b + Delta_2*b**2*p - a) - rho*(Delta_1*b*p + Delta_2*b*p - R*T - b*p))/rho**3 - + // Cubic coefficients: double crho0 = -p; - double crho1 = -1*((Delta_1+Delta_2-1)*bm*p - R*T); - double crho2 = -Delta_1*Delta_2*bm*bm*p + (Delta_1+Delta_2)*(R*T*bm + bm*bm*p) - am; - double crho3 = bm*(Delta_1*Delta_2*(R*T*bm + bm*bm*p) + am); + double crho1 = R*T - p*(d1+d2+d3); + double crho2 = R*T*(d2 + d3) - p*(d1*(d2 + d3) + d2*d3) - am; + double crho3 = R*T*d2*d3 - p*d1*d2*d3 - d1*am; + + // Solving the cubic: solve_cubic(crho3, crho2, crho1, crho0, Nsolns, rho0, rho1, rho2); sort3(rho0, rho1, rho2); return; - -// double A = cubic->am_term(cubic->T_r/T, mole_fractions_double, 0)*p/(POW2(R*T)); -// double B = cubic->bm_term(mole_fractions)*p/(R*T); -// double Z0=0, Z1=0, Z2=0; -// solve_cubic(1, -// B*(Delta_1+Delta_2-1)-1, -// A + B*B*(Delta_1*Delta_2-Delta_1-Delta_2) - B*(Delta_1+Delta_2), -// -A*B-Delta_1*Delta_2*(POW2(B)+POW3(B)), -// Nsolns, Z0, Z1, Z2); -// if (Nsolns == 1){ rho0 = p/(Z0*R*T); } -// else if (Nsolns == 3){ -// rho0 = p/(Z0*R*T); -// rho1 = p/(Z1*R*T); -// rho2 = p/(Z2*R*T); -// sort3(rho0, rho1, rho2); -// } } class SaturationResidual : public CoolProp::FuncWrapper1D{ @@ -507,3 +489,16 @@ void CoolProp::AbstractCubicBackend::set_C_Twu(double L, double M, double N){ ACB->set_C_Twu(L,M,N); } } + +void CoolProp::AbstractCubicBackend::set_fluid_parameter_double(const size_t i, const std::string parameter, const double value) +{ + // Set the volume translation parrameter, currently applied to the whole fluid, not to components. + if (parameter == "c" || parameter == "cm" || parameter == "c_m") { + get_cubic()->set_cm(value); + if (this->SatL.get() != NULL) { this->SatL->set_fluid_parameter_double(i, "cm", value); } + if (this->SatV.get() != NULL) { this->SatV->set_fluid_parameter_double(i, "cm", value); } + } + else { + throw ValueError(format("I don't know what to do with parameter [%s]", parameter.c_str())); + } +} diff --git a/src/Backends/Cubics/CubicBackend.h b/src/Backends/Cubics/CubicBackend.h index 8d4ddf8f..3485ac9a 100644 --- a/src/Backends/Cubics/CubicBackend.h +++ b/src/Backends/Cubics/CubicBackend.h @@ -186,6 +186,9 @@ public: // Set the Twu constants L,M,N for a pure fluid void set_C_Twu(double L, double M, double N); + + // Set fluid parameter (currently the volume translation parameter) + void set_fluid_parameter_double(const size_t i, const std::string parameter, const double value); }; diff --git a/src/Backends/Cubics/GeneralizedCubic.cpp b/src/Backends/Cubics/GeneralizedCubic.cpp index 5967a55c..28a4e3f8 100644 --- a/src/Backends/Cubics/GeneralizedCubic.cpp +++ b/src/Backends/Cubics/GeneralizedCubic.cpp @@ -81,6 +81,11 @@ double AbstractCubic::d3_bm_term_dxidxjdxk(const std::vector &x, std::si return 0; } +double AbstractCubic::cm_term() +{ + return cm; +} + double AbstractCubic::aii_term(double tau, std::size_t i, std::size_t itau) { double Tr_over_Tci = T_r/Tc[i]; @@ -247,20 +252,20 @@ double AbstractCubic::aij_term(double tau, std::size_t i, std::size_t j, std::si double AbstractCubic::psi_minus(double delta, const std::vector &x, std::size_t itau, std::size_t idelta) { if (itau > 0) return 0.0; - double b = bm_term(x); - double bracket = 1-b*delta*rho_r; + double bmc = bm_term(x)-cm_term(); // appears only in the form (b-c) in the equations + double bracket = 1-bmc*delta*rho_r; switch(idelta){ case 0: return -log(bracket); case 1: - return b*rho_r/bracket; + return bmc*rho_r/bracket; case 2: - return pow(b*rho_r/bracket, 2); + return pow(bmc*rho_r/bracket, 2); case 3: - return 2*pow(b*rho_r/bracket, 3); + return 2*pow(bmc*rho_r/bracket, 3); case 4: - return 6*pow(b*rho_r/bracket, 4); + return 6*pow(bmc*rho_r/bracket, 4); default: throw -1; } @@ -268,9 +273,9 @@ double AbstractCubic::psi_minus(double delta, const std::vector &x, std: double AbstractCubic::d_psi_minus_dxi(double delta, const std::vector &x, std::size_t itau, std::size_t idelta, std::size_t i, bool xN_independent) { if (itau > 0) return 0.0; - double b = bm_term(x); + double bmc = bm_term(x) - cm_term(); // appears only in the form (b-c) in the equations double db_dxi = d_bm_term_dxi(x, i, xN_independent); - double bracket = 1-b*delta*rho_r; + double bracket = 1-bmc*delta*rho_r; switch(idelta){ case 0: @@ -278,11 +283,11 @@ double AbstractCubic::d_psi_minus_dxi(double delta, const std::vector &x case 1: return rho_r*db_dxi/pow(bracket, 2); case 2: - return 2*pow(rho_r,2)*b*db_dxi/pow(bracket, 3); + return 2*pow(rho_r,2)*bmc*db_dxi/pow(bracket, 3); case 3: - return 6*pow(rho_r,3)*pow(b, 2)*db_dxi/pow(bracket, 4); + return 6*pow(rho_r,3)*pow(bmc, 2)*db_dxi/pow(bracket, 4); case 4: - return 24*pow(rho_r,4)*pow(b, 3)*db_dxi/pow(bracket, 5); + return 24*pow(rho_r,4)*pow(bmc, 3)*db_dxi/pow(bracket, 5); default: throw -1; } @@ -290,11 +295,11 @@ double AbstractCubic::d_psi_minus_dxi(double delta, const std::vector &x double AbstractCubic::d2_psi_minus_dxidxj(double delta, const std::vector &x, std::size_t itau, std::size_t idelta, std::size_t i, std::size_t j, bool xN_independent) { if (itau > 0) return 0.0; - double b = bm_term(x); + double bmc = bm_term(x) - cm_term(); // appears only in the form (b-c) in the equations double db_dxi = d_bm_term_dxi(x, i, xN_independent), db_dxj = d_bm_term_dxi(x, j, xN_independent), d2b_dxidxj = d2_bm_term_dxidxj(x, i, j, xN_independent); - double bracket = 1-b*delta*rho_r; + double bracket = 1-bmc*delta*rho_r; switch(idelta){ case 0: @@ -302,11 +307,11 @@ double AbstractCubic::d2_psi_minus_dxidxj(double delta, const std::vector &x, std::size_t itau, std::size_t idelta, std::size_t i, std::size_t j, std::size_t k, bool xN_independent) { if (itau > 0) return 0.0; - double b = bm_term(x); + double bmc = bm_term(x) - cm_term(); // appears only in the form (b-c) in the equations double db_dxi = d_bm_term_dxi(x, i, xN_independent), db_dxj = d_bm_term_dxi(x, j, xN_independent), db_dxk = d_bm_term_dxi(x, k, xN_independent), @@ -322,7 +327,7 @@ double AbstractCubic::d3_psi_minus_dxidxjdxk(double delta, const std::vector &x, std::size_t idelta) { - double b = bm_term(x); + double bm = bm_term(x); + double cm = cm_term(); switch(idelta){ case 0: - return (1+Delta_1*b*rho_r*delta)*(1+Delta_2*b*rho_r*delta); + return (1+(Delta_1*bm+cm)*rho_r*delta)*(1+(Delta_2*bm+cm)*rho_r*delta); case 1: - return b*rho_r*(2*Delta_1*Delta_2*b*delta*rho_r+Delta_1+Delta_2); + return rho_r*(2*(bm*Delta_1 + cm)*(bm*Delta_2 + cm)*delta*rho_r + (Delta_1 + Delta_2)*bm + 2*cm); case 2: - return 2*Delta_1*Delta_2*pow(b*rho_r, 2); + return 2*(Delta_1*bm+cm)*(Delta_2*bm + cm)*pow(rho_r, 2); case 3: return 0; case 4: @@ -361,15 +367,16 @@ double AbstractCubic::PI_12(double delta, const std::vector &x, std::siz } double AbstractCubic::d_PI_12_dxi(double delta, const std::vector &x, std::size_t idelta, std::size_t i, bool xN_independent) { - double b = bm_term(x); + double bm = bm_term(x); + double cm = cm_term(); double db_dxi = d_bm_term_dxi(x, i, xN_independent); switch(idelta){ case 0: - return delta*rho_r*db_dxi*(2*Delta_1*Delta_2*b*delta*rho_r+Delta_1+Delta_2); + return delta*rho_r*db_dxi*(2*Delta_1*Delta_2*bm*delta*rho_r+(Delta_1+Delta_2)*(1+cm*delta*rho_r)); case 1: - return rho_r*db_dxi*(4*Delta_1*Delta_2*b*delta*rho_r+Delta_1+Delta_2); + return rho_r*db_dxi*(4*Delta_1*Delta_2*bm*delta*rho_r+(Delta_1 + Delta_2)*(1+2*cm*delta*rho_r)); case 2: - return 4*Delta_1*Delta_2*pow(rho_r, 2)*b*db_dxi; + return 2*pow(rho_r, 2)*(2*Delta_1*Delta_2*bm+ (Delta_1+Delta_2)*cm)*db_dxi; case 3: return 0; case 4: @@ -380,17 +387,18 @@ double AbstractCubic::d_PI_12_dxi(double delta, const std::vector &x, st } double AbstractCubic::d2_PI_12_dxidxj(double delta, const std::vector &x, std::size_t idelta, std::size_t i, std::size_t j, bool xN_independent) { - double b = bm_term(x); + double bm = bm_term(x); + double cm = cm_term(); double db_dxi = d_bm_term_dxi(x, i, xN_independent), db_dxj = d_bm_term_dxi(x, j, xN_independent), d2b_dxidxj = d2_bm_term_dxidxj(x, i, j, xN_independent); switch(idelta){ case 0: - return delta*rho_r*(2*Delta_1*Delta_2*delta*rho_r*db_dxi*db_dxj + (2*Delta_1*Delta_2*delta*rho_r*b+Delta_1+Delta_2)*d2b_dxidxj); + return delta*rho_r*(2*Delta_1*Delta_2*delta*rho_r*db_dxi*db_dxj + (2* Delta_1*Delta_2*bm*delta*rho_r + (Delta_1 + Delta_2)*(1 + cm*delta*rho_r))*d2b_dxidxj); case 1: - return rho_r*(4*Delta_1*Delta_2*delta*rho_r*db_dxi*db_dxj + (4*Delta_1*Delta_2*delta*rho_r*b+Delta_1+Delta_2)*d2b_dxidxj); + return rho_r*(4*Delta_1*Delta_2*delta*rho_r*db_dxi*db_dxj + (4*Delta_1*Delta_2*bm*delta*rho_r + (Delta_1 + Delta_2)*(1 + 2*cm*delta*rho_r))*d2b_dxidxj); case 2: - return 4*Delta_1*Delta_2*pow(rho_r,2)*(db_dxi*db_dxj + b*d2b_dxidxj); + return 2*pow(rho_r,2)*(2 * Delta_1*Delta_2*db_dxi*db_dxj + (2*Delta_1*Delta_2*bm + (Delta_1 + Delta_2)*cm)*d2b_dxidxj); case 3: return 0; case 4: @@ -401,7 +409,8 @@ double AbstractCubic::d2_PI_12_dxidxj(double delta, const std::vector &x } double AbstractCubic::d3_PI_12_dxidxjdxk(double delta, const std::vector &x, std::size_t idelta, std::size_t i, std::size_t j, std::size_t k, bool xN_independent) { - double b = bm_term(x); + double bm = bm_term(x); + double cm = cm_term(); double db_dxi = d_bm_term_dxi(x, i, xN_independent), db_dxj = d_bm_term_dxi(x, j, xN_independent), db_dxk = d_bm_term_dxi(x, k, xN_independent), @@ -411,14 +420,14 @@ double AbstractCubic::d3_PI_12_dxidxjdxk(double delta, const std::vector d3b_dxidxjdxk = d3_bm_term_dxidxjdxk(x, i, j, k, xN_independent); switch(idelta){ case 0: - return delta*rho_r*((2*Delta_1*Delta_2*delta*rho_r*b+Delta_1+Delta_2)*d3b_dxidxjdxk + return delta*rho_r*((2 * Delta_1*Delta_2*bm*delta*rho_r + (Delta_1 + Delta_2)*(1 + cm*delta*rho_r))*d3b_dxidxjdxk + 2*Delta_1*Delta_2*delta*rho_r*(db_dxi*d2b_dxjdxk +db_dxj*d2b_dxidxk +db_dxk*d2b_dxidxj ) ); case 1: - return rho_r*((4*Delta_1*Delta_2*delta*rho_r*b+Delta_1+Delta_2)*d3b_dxidxjdxk + return rho_r*((4.*Delta_1*Delta_2*bm*delta*rho_r + (Delta_1 + Delta_2)*(1 + 2*cm*delta*rho_r))*d3b_dxidxjdxk + 4*Delta_1*Delta_2*delta*rho_r*(db_dxi*d2b_dxjdxk + db_dxj*d2b_dxidxk + db_dxk*d2b_dxidxj @@ -692,4 +701,4 @@ double PengRobinson::m_ii(std::size_t i) double omega = acentric[i]; double m = 0.37464 + 1.54226*omega - 0.26992*omega*omega; return m; -} \ No newline at end of file +} diff --git a/src/Backends/Cubics/GeneralizedCubic.h b/src/Backends/Cubics/GeneralizedCubic.h index d3dd97b9..6e12ae34 100644 --- a/src/Backends/Cubics/GeneralizedCubic.h +++ b/src/Backends/Cubics/GeneralizedCubic.h @@ -30,6 +30,7 @@ protected: std::vector C1, C2, C3; ///< The Mathias-Copeman coefficients for a_ii std::vector L_Twu, M_Twu, N_Twu; ///< The Twu coefficients for a_ii AII_Model aii_model; ///< Enumeration for the aii model in use + double cm; ///< The volume translation parameter public: static const double rho_r, T_r; /** @@ -55,6 +56,7 @@ public: { N = static_cast(Tc.size()); k.resize(N, std::vector(N, 0)); + cm = 0.; /// If no Mathias-Copeman coefficients are passed in (all empty vectors), use the predictive scheme for m_ii simple_aii = (C1.empty() && C2.empty() && C3.empty() && L_Twu.empty() && M_Twu.empty() && N_Twu.empty()); }; @@ -118,13 +120,13 @@ public: virtual double m_ii(std::size_t i) = 0; /// The residual non-dimensionalized Helmholtz energy \f$\alpha^r\f$ - double alphar(double tau, double delta, const std::vector &x, std::size_t itau, std::size_t idelta); + virtual double alphar(double tau, double delta, const std::vector &x, std::size_t itau, std::size_t idelta); /// The first composition derivative of \f$\alpha^r\f$ as well as derivatives with respect to \f$\tau\f$ and \f$\delta\f$ - double d_alphar_dxi(double tau, double delta, const std::vector &x, std::size_t itau, std::size_t idelta, std::size_t i, bool xN_independent); + virtual double d_alphar_dxi(double tau, double delta, const std::vector &x, std::size_t itau, std::size_t idelta, std::size_t i, bool xN_independent); /// The second composition derivative of \f$\alpha^r\f$ as well as derivatives with respect to \f$\tau\f$ and \f$\delta\f$ - double d2_alphar_dxidxj(double tau, double delta, const std::vector &x, std::size_t itau, std::size_t idelta, std::size_t i, std::size_t j, bool xN_independent); + virtual double d2_alphar_dxidxj(double tau, double delta, const std::vector &x, std::size_t itau, std::size_t idelta, std::size_t i, std::size_t j, bool xN_independent); /// The third composition derivative of \f$\alpha^r\f$ as well as derivatives with respect to \f$\tau\f$ and \f$\delta\f$ - double d3_alphar_dxidxjdxk(double tau, double delta, const std::vector &x, std::size_t itau, std::size_t idelta, std::size_t i, std::size_t j, std::size_t k, bool xN_independent); + virtual double d3_alphar_dxidxjdxk(double tau, double delta, const std::vector &x, std::size_t itau, std::size_t idelta, std::size_t i, std::size_t j, std::size_t k, bool xN_independent); /** * \brief The n-th derivative of \f$a_m\f$ with respect to \f$\tau\f$ @@ -190,7 +192,15 @@ public: * \param xN_independent True if \f$x_N\f$ is an independent variable, false otherwise (dependent on other \f$N-1\f$ mole fractions) */ virtual double d3_bm_term_dxidxjdxk(const std::vector &x, std::size_t i, std::size_t j, std::size_t k, bool xN_independent); - + /** + * \brief The term \f$c_{\rm m}\f$ (volume translation) + */ + virtual double cm_term(); + /// Set the volume translation parameter + void set_cm(double val) { cm = val; } + +protected: + /** * \brief The n-th \f$\tau\f$ derivative of \f$a_{ij}(\tau)\f$ * \param tau The reciprocal reduced temperature \f$\tau=T_r/T\f$ @@ -401,8 +411,8 @@ public: * \param xN_independent True if \f$x_N\f$ is an independent variable, false otherwise (dependent on other \f$N-1\f$ mole fractions) */ double d2_c_term_dxidxj(const std::vector &x, std::size_t i, std::size_t j, bool xN_independent){ - double b = bm_term(x); - return (2*d_bm_term_dxi(x, i, xN_independent)*d_bm_term_dxi(x, j, xN_independent) - b*d2_bm_term_dxidxj(x, i,j,xN_independent))/pow(b, 3); + double bm = bm_term(x); + return (2*d_bm_term_dxi(x, i, xN_independent)*d_bm_term_dxi(x, j, xN_independent) - bm*d2_bm_term_dxidxj(x, i,j,xN_independent))/pow(bm, 3); }; /** * \brief The third composition derivative of the term \f$c\f$ used in the pure composition partial derivatives of \f$\psi^{(+)}\f$ @@ -413,12 +423,12 @@ public: * \param xN_independent True if \f$x_N\f$ is an independent variable, false otherwise (dependent on other \f$N-1\f$ mole fractions) */ double d3_c_term_dxidxjdxk(const std::vector &x, std::size_t i, std::size_t j, std::size_t k, bool xN_independent){ - double b = bm_term(x); - return 1/pow(b,4)*(2*b*(d_bm_term_dxi(x, i, xN_independent)*d2_bm_term_dxidxj(x, j, k, xN_independent) + double bm = bm_term(x); + return 1/pow(bm,4)*(2*bm*(d_bm_term_dxi(x, i, xN_independent)*d2_bm_term_dxidxj(x, j, k, xN_independent) +d_bm_term_dxi(x, j, xN_independent)*d2_bm_term_dxidxj(x, i, k, xN_independent) +d_bm_term_dxi(x, k, xN_independent)*d2_bm_term_dxidxj(x, i, j, xN_independent) ) - - pow(b,2)*d3_bm_term_dxidxjdxk(x, i,j,k,xN_independent) + - pow(bm,2)*d3_bm_term_dxidxjdxk(x, i,j,k,xN_independent) -6*d_bm_term_dxi(x, i, xN_independent)*d_bm_term_dxi(x, j, xN_independent)*d_bm_term_dxi(x, k, xN_independent) ); }; @@ -434,8 +444,9 @@ public: * \param x The vector of mole fractions */ double A_term(double delta, const std::vector &x){ - double b = bm_term(x); - return log((Delta_1*delta*rho_r*b+1)/(Delta_2*delta*rho_r*b+1)); + double bm = bm_term(x); + double cm = cm_term(); + return log((delta*rho_r*(Delta_1*bm+cm)+1)/(delta*rho_r*(Delta_2*bm + cm) +1)); }; /** * \brief The first composition derivative of the term \f$A\f$ used in the pure composition partial derivatives of \f$\psi^{(+)}\f$ diff --git a/src/Backends/Cubics/UNIFAQ.cpp b/src/Backends/Cubics/UNIFAQ.cpp new file mode 100644 index 00000000..e5835441 --- /dev/null +++ b/src/Backends/Cubics/UNIFAQ.cpp @@ -0,0 +1,231 @@ +#include "UNIFAQ.h" + +void UNIFAQ::UNIFAQMixture::set_interaction_parameters() { + for (int i = 0; i < unique_groups.size(); ++i) { + for (int j = i + 1; j < unique_groups.size(); ++j) { + int mgi1 = unique_groups[i].mgi, mgi2 = unique_groups[j].mgi; + std::pair< std::pair, UNIFAQLibrary::InteractionParameters> m_pair(std::pair(mgi1, mgi2), library.get_interaction_parameters(mgi1, mgi2)); + interaction.insert(m_pair); + } + } +} + +/// Set the mole fractions of the components in the mixtures (not the groups) +void UNIFAQ::UNIFAQMixture::set_mole_fractions(const std::vector &z) { + pure_data.clear(); + this->mole_fractions = z; + std::size_t N = z.size(); + std::vector &r = m_r, &q = m_q, &l = m_l, &phi = m_phi, &theta = m_theta, &ln_Gamma_C = m_ln_Gamma_C; + r.resize(N); q.resize(N); l.resize(N); phi.resize(N); theta.resize(N); ln_Gamma_C.resize(N); + double summerzr = 0, summerzq = 0, summerzl = 0; + for (std::size_t i = 0; i < z.size(); ++i) { + double summerr = 0, summerq = 0; + const UNIFAQLibrary::Component &c = components[i]; + for (std::size_t j = 0; j < c.groups.size(); ++j) { + const UNIFAQLibrary::ComponentGroup &cg = c.groups[j]; + summerr += cg.count*cg.group.R_k; + summerq += cg.count*cg.group.Q_k; + } + r[i] = summerr; + q[i] = summerq; + summerzr += z[i] * r[i]; + summerzq += z[i] * q[i]; + } + for (std::size_t i = 0; i < z.size(); ++i) { + phi[i] = z[i] * r[i] / summerzr; + theta[i] = z[i] * q[i] / summerzq; + l[i] = 10.0 / 2.0*(r[i] - q[i]) - (r[i] - 1); + summerzl += z[i] * l[i]; + } + for (std::size_t i = 0; i < z.size(); ++i) { + ln_Gamma_C[i] = log(phi[i]/z[i]) + 10.0/2.0*q[i] * log(theta[i]/phi[i]) + l[i] - phi[i]/z[i]*summerzl; + } + + /// Calculate the parameters X and theta for the pure components, which does not depend on temperature + for (std::size_t i = 0; i < z.size(); ++i){ + int totalgroups = 0; + const UNIFAQLibrary::Component &c = components[i]; + ComponentData cd; + double summerxq = 0; + cd.group_count = 0; + for (std::size_t j = 0; j < c.groups.size(); ++j) { + const UNIFAQLibrary::ComponentGroup &cg = c.groups[j]; + double x = static_cast(cg.count); + double theta = static_cast(cg.count*cg.group.Q_k); + cd.X.insert( std::pair(cg.group.sgi, x) ); + cd.theta.insert(std::pair(cg.group.sgi, theta)); + cd.group_count += cg.count; + totalgroups += cg.count; + summerxq += x*cg.group.Q_k; + } + /// Now come back through and divide by the total # groups for this fluid + for (std::map::iterator it = cd.X.begin(); it != cd.X.end(); ++it) { + it->second /= totalgroups; + //printf("X^(%d)_{%d}: %g\n", static_cast(i + 1), static_cast(it->first), it->second); + } + /// Now come back through and divide by the sum(X*Q) for this fluid + for (std::map::iterator it = cd.theta.begin(); it != cd.theta.end(); ++it){ + it->second /= summerxq; + //printf("theta^(%d)_{%d}: %g\n", static_cast(i+1), static_cast(it->first), it->second); + } + pure_data.push_back(cd); + } + for (std::size_t i = 0; i < z.size(); ++i) { + //printf("%g %g %g %g %g %g\n", l[i], phi[i], q[i], r[i], theta[i], ln_Gamma_C[i]); + } +} + +double UNIFAQ::UNIFAQMixture::Psi(std::size_t sgi1, std::size_t sgi2) const { + std::size_t mgi1 = m_sgi_to_mgi.find(sgi1)->second; + std::size_t mgi2 = m_sgi_to_mgi.find(sgi2)->second; + std::map, UNIFAQLibrary::InteractionParameters>::const_iterator it = this->interaction.find(std::pair(mgi1,mgi2)); + if (it != this->interaction.end()){ + return exp(-(it->second.a_ij + it->second.b_ij*this->m_T + it->second.c_ij*this->m_T*this->m_T)/this->m_T); + } + else{ + return 1; + } +} + +std::size_t UNIFAQ::UNIFAQMixture::group_count(std::size_t i, std::size_t sgi) const { + const UNIFAQLibrary::Component &c = components[i]; + for (std::vector::const_iterator it = c.groups.begin(); it != c.groups.end(); ++it){ + if (it->group.sgi == sgi){ return it->count; } + } + return 0; +} + +double UNIFAQ::UNIFAQMixture::theta_pure(std::size_t i, std::size_t sgi) const { + return pure_data[i].theta.find(sgi)->second; +} + +void UNIFAQ::UNIFAQMixture::set_temperature(const double T){ + this->m_T = T; + for (std::size_t i = 0; i < this->mole_fractions.size(); ++i) { + const UNIFAQLibrary::Component &c = components[i]; + for (std::size_t k = 0; k < c.groups.size(); ++k) { + double Q = c.groups[k].group.Q_k; + int sgik = c.groups[k].group.sgi; + double sum1 = 0; + for (std::size_t m = 0; m < c.groups.size(); ++m) { + int sgim = c.groups[m].group.sgi; + sum1 += theta_pure(i, sgim)*Psi(sgim, sgik); + } + double s = 1 - log(sum1); + for (std::size_t m = 0; m < c.groups.size(); ++m) { + int sgim = c.groups[m].group.sgi; + double sum2 = 0; + for (std::size_t n = 0; n < c.groups.size(); ++n) { + int sgin = c.groups[n].group.sgi; + sum2 += theta_pure(i, sgin)*Psi(sgin, sgim); + } + s -= theta_pure(i, sgim)*Psi(sgik, sgim)/sum2; + } + ComponentData &cd = pure_data[i]; + cd.lnGamma.insert(std::pair(sgik, Q*s)); + //printf("ln(Gamma)^(%d)_{%d}: %g\n", static_cast(i + 1), sgik, Q*s); + } + } + + if (this->mole_fractions.empty()){ + throw CoolProp::ValueError("mole fractions must be set before calling set_temperature"); + } + + std::map &Xg = m_Xg, &thetag = m_thetag, &lnGammag = m_lnGammag; + Xg.clear(); thetag.clear(); lnGammag.clear(); + + // Iterate over the fluids + double X_summer = 0; + for (std::size_t i = 0; i < this->mole_fractions.size(); ++i) { + X_summer += this->mole_fractions[i] * pure_data[i].group_count; + } + /// Calculations for each group in the total mixture + for (std::vector::iterator it = unique_groups.begin(); it != unique_groups.end(); ++it){ + double X = 0; + // Iterate over the fluids + for (std::size_t i = 0; i < this->mole_fractions.size(); ++i) { + X += this->mole_fractions[i]*group_count(i, it->sgi); + } + Xg.insert(std::pair(it->sgi, X)); + } + /// Now come back through and divide by the sum(z_i*count) for this fluid + for (std::map::iterator it = Xg.begin(); it != Xg.end(); ++it) { + it->second /= X_summer; + //printf("X_{%d}: %g\n", it->first, it->second); + } + double theta_summer = 0; + for (std::vector::iterator it = unique_groups.begin(); it != unique_groups.end(); ++it) { + double cont = Xg.find(it->sgi)->second * it->Q_k; + theta_summer += cont; + thetag.insert(std::pair(it->sgi, cont)); + } + /// Now come back through and divide by the sum(X*Q) for this fluid + for (std::map::iterator it = thetag.begin(); it != thetag.end(); ++it) { + it->second /= theta_summer; + //printf("theta_{%d}: %g\n", it->first, it->second); + } + + for (std::vector::iterator itk = unique_groups.begin(); itk != unique_groups.end(); ++itk) { + double sum1 = 0; + for (std::vector::iterator itm = unique_groups.begin(); itm != unique_groups.end(); ++itm) { + sum1 += thetag.find(itm->sgi)->second*this->Psi(itm->sgi, itk->sgi); + } + double s = 1-log(sum1); + for (std::vector::iterator itm = unique_groups.begin(); itm != unique_groups.end(); ++itm) { + double sum3 = 0; + for (std::vector::iterator itn = unique_groups.begin(); itn != unique_groups.end(); ++itn) { + sum3 += thetag.find(itn->sgi)->second*this->Psi(itn->sgi, itm->sgi); + } + s -= thetag.find(itm->sgi)->second*Psi(itk->sgi, itm->sgi)/sum3; + } + lnGammag.insert(std::pair(itk->sgi, itk->Q_k*s)); + //printf("log(Gamma)_{%d}: %g\n", itk->sgi, itk->Q_k*s); + } +} +double UNIFAQ::UNIFAQMixture::ln_gamma_R(std::size_t i) const{ + double summer = 0; + for (std::vector::const_iterator it = unique_groups.begin(); it != unique_groups.end(); ++it) { + std::size_t k = it->sgi; + std::size_t count = group_count(i, k); + if (count > 0){ + summer += count*(m_lnGammag.find(k)->second - pure_data[i].lnGamma.find(k)->second); + } + } + //printf("log(gamma)_{%d}: %g\n", i+1, summer); + return summer; +} +double UNIFAQ::UNIFAQMixture::activity_coefficient(std::size_t i) const { + return exp(ln_gamma_R(i) + m_ln_Gamma_C[i]); +} + +/// Add a component with the defined groups defined by (count, sgi) pairs +void UNIFAQ::UNIFAQMixture::add_component(const UNIFAQLibrary::Component &comp) { + components.push_back(comp); + // Check if you also need to add group into list of unique groups + for (std::vector::const_iterator it = comp.groups.begin(); it != comp.groups.end(); ++it) { + bool insert_into_unique = true; + // if already in unique_groups, don't save it, go to next one + for (std::vector::const_iterator it2 = unique_groups.begin(); it2 != unique_groups.end(); ++it2) { + if (it2->sgi == it->group.sgi) { insert_into_unique = false; break; } + } + if (insert_into_unique) { + unique_groups.push_back(it->group); + m_sgi_to_mgi.insert(std::pair(it->group.sgi, it->group.mgi)); + } + } +} + +void UNIFAQ::UNIFAQMixture::set_components(const std::string &identifier_type, std::vector identifiers) { + components.clear(); + if (identifier_type == "name") { + // Iterate over the provided names + for (std::vector::const_iterator it = identifiers.begin(); it != identifiers.end(); ++it) { + // Get and add the component + UNIFAQLibrary::Component c = library.get_component("name", *it); + add_component(c); + } + } + else { + throw CoolProp::ValueError("Cannot understand identifier_type"); + } +} \ No newline at end of file diff --git a/src/Backends/Cubics/UNIFAQ.h b/src/Backends/Cubics/UNIFAQ.h new file mode 100644 index 00000000..a9cbed0a --- /dev/null +++ b/src/Backends/Cubics/UNIFAQ.h @@ -0,0 +1,94 @@ +#ifndef UNIFAQ_H_ +#define UNIFAQ_H_ + +#include + +#include "UNIFAQLibrary.h" + +/// Structure containing data for the pure fluid in the mixture +struct ComponentData { + std::map X, theta, lnGamma; + int group_count; ///< The total number of groups in the pure fluid +}; + +namespace UNIFAQ +{ + class UNIFAQMixture + { + private: + double m_T; ///< The temperature in K + + std::vector m_r, + m_q, + m_l, + m_phi, + m_theta, + m_ln_Gamma_C; + + std::map m_Xg, ///< Map from sgi to mole fraction of group in the mixture + m_thetag, ///< Map from sgi to theta for the group in the mixture + m_lnGammag; ///< Map from sgi to ln(Gamma) for the group in the mixture + + /// A const reference to the library of group and interaction parameters + const UNIFAQLibrary::UNIFAQParameterLibrary &library; + + /// A map from (i, j) indices for subgroup, subgroup indices to the interaction parameters for this pair + std::map, UNIFAQLibrary::InteractionParameters> interaction; + + /// A map from SGI to MGI + std::map m_sgi_to_mgi; + + /// A vector of unique groups in this mixture + std::vector unique_groups; + + std::vector mole_fractions; + + std::vector components; + + std::vector pure_data; + + public: + + UNIFAQMixture(const UNIFAQLibrary::UNIFAQParameterLibrary &library) : library(library) {}; + + /** + * \brief Set all the interaction parameters between groups + * + * \param subgroups A vector of the set of the unique Group forming the mixture - these + * permutations represent the set of posisble binary interactions + */ + void set_interaction_parameters(); + + /// Set the mole fractions of the components in the mixtures (not the groups) + void set_mole_fractions(const std::vector &z); + + /// Get the mole fractions of the components in the mixtures (not the groups) + const std::vector & get_mole_fractions() { return mole_fractions; } + + /// Set the mole fractions of the components in the mixtures (not the groups) + void set_temperature(const double T); + + /// Get the temperature + double get_temperature() const { return m_T; } + + double Psi(std::size_t sgi1, std::size_t sgi2) const; + + double theta_pure(std::size_t i, std::size_t sgi) const; + + double activity_coefficient(std::size_t i) const; + + double ln_gamma_R(std::size_t i) const; + + std::size_t group_count(std::size_t i, std::size_t sgi) const; + + /// Add a component with the defined groups defined by (count, sgi) pairs + void add_component(const UNIFAQLibrary::Component &comp); + + void set_components(const std::string &identifier_type, std::vector identifiers); + + const std::vector & get_components() { return components; }; + }; + +} /* namespace UNIFAQ */ + +#endif \ No newline at end of file diff --git a/src/Backends/Cubics/UNIFAQLibrary.cpp b/src/Backends/Cubics/UNIFAQLibrary.cpp new file mode 100644 index 00000000..82966c04 --- /dev/null +++ b/src/Backends/Cubics/UNIFAQLibrary.cpp @@ -0,0 +1,118 @@ +#include "UNIFAQLibrary.h" + +namespace UNIFAQLibrary{ + + void UNIFAQParameterLibrary::jsonize(std::string &s, rapidjson::Document &d) + { + d.Parse<0>(s.c_str()); + if (d.HasParseError()) { + throw -1; + } + else { + return; + } + } + void UNIFAQParameterLibrary::populate(rapidjson::Value &group_data, rapidjson::Value &interaction_data, rapidjson::Value &comp_data) + { + // Schema should have been used to validate the data already, so by this point we are can safely consume the data without checking ... + for (rapidjson::Value::ValueIterator itr = group_data.Begin(); itr != group_data.End(); ++itr) + { + Group g; + g.sgi = (*itr)["sgi"].GetInt(); + g.mgi = (*itr)["mgi"].GetInt(); + g.R_k = (*itr)["R_k"].GetDouble(); + g.Q_k = (*itr)["Q_k"].GetDouble(); + groups.push_back(g); + } + for (rapidjson::Value::ValueIterator itr = interaction_data.Begin(); itr != interaction_data.End(); ++itr) + { + InteractionParameters ip; + ip.mgi1 = (*itr)["mgi1"].GetInt(); + ip.mgi2 = (*itr)["mgi2"].GetInt(); + ip.a_ij = (*itr)["a_ij"].GetDouble(); + ip.a_ji = (*itr)["a_ji"].GetDouble(); + ip.b_ij = (*itr)["b_ij"].GetDouble(); + ip.b_ji = (*itr)["b_ji"].GetDouble(); + ip.c_ij = (*itr)["c_ij"].GetDouble(); + ip.c_ji = (*itr)["c_ji"].GetDouble(); + interaction_parameters.push_back(ip); + } + for (rapidjson::Value::ValueIterator itr = comp_data.Begin(); itr != comp_data.End(); ++itr) + { + Component c; + c.inchikey = (*itr)["inchikey"].GetString(); + c.registry_number = (*itr)["registry_number"].GetString(); + c.name = (*itr)["name"].GetString(); + c.Tc = (*itr)["Tc"].GetDouble(); + c.pc = (*itr)["pc"].GetDouble(); + c.acentric = (*itr)["acentric"].GetDouble(); + // userid is an optional user identifier + if ((*itr).HasMember("userid")){ + c.userid = (*itr)["userid"].GetString(); + } + rapidjson::Value &groups = (*itr)["groups"]; + for (rapidjson::Value::ValueIterator itrg = groups.Begin(); itrg != groups.End(); ++itrg) + { + int count = (*itrg)["count"].GetInt(); + int sgi = (*itrg)["sgi"].GetInt(); + if (has_group(sgi)){ + ComponentGroup cg(count, get_group(sgi)); + c.groups.push_back(cg); + } + } + components.push_back(c); + } + } + void UNIFAQParameterLibrary::populate(std::string &group_data, std::string &interaction_data, std::string &decomp_data) + { + rapidjson::Document group_JSON; jsonize(group_data, group_JSON); + rapidjson::Document interaction_JSON; jsonize(interaction_data, interaction_JSON); + rapidjson::Document decomp_JSON; jsonize(decomp_data, decomp_JSON); + populate(group_JSON, interaction_JSON, decomp_JSON); + } + Group UNIFAQParameterLibrary::get_group(int sgi) const { + for (std::vector::const_iterator it = groups.begin(); it != groups.end(); ++it) { + if (it->sgi == sgi) { return *it; } + } + throw CoolProp::ValueError("Could not find group"); + } + bool UNIFAQParameterLibrary::has_group(int sgi) const { + for (std::vector::const_iterator it = groups.begin(); it != groups.end(); ++it) { + if (it->sgi == sgi) { return true; } + } + return false; + } + + InteractionParameters UNIFAQParameterLibrary::get_interaction_parameters(int mgi1, int mgi2) const { + + // If both mgi are the same, yield all zeros for the interaction parameters + if (mgi1 == mgi2){ + InteractionParameters ip; ip.mgi1 = mgi1; ip.mgi2 = mgi2; + ip.zero_out(); + return ip; + } + for (std::vector::const_iterator it = interaction_parameters.begin(); it != interaction_parameters.end(); ++it) { + if (it->mgi1 == mgi1 && it->mgi2 == mgi2) { + // Correct order, return it + return *it; + } + if (it->mgi2 == mgi1 && it->mgi1 == mgi2) { + // Backwards, swap the parameters + InteractionParameters ip = *it; + ip.swap(); + return ip; + } + } + throw CoolProp::ValueError("Could not find interaction pair"); + } + + Component UNIFAQParameterLibrary::get_component(const std::string &identifier, const std::string &value) const { + if (identifier == "name"){ + for (std::vector::const_iterator it = components.begin(); it != components.end(); ++it ){ + if (it->name == value ){ return *it; } + } + } + throw CoolProp::ValueError("Could not find component"); + } + +}; /* namespace UNIFAQLibrary */ \ No newline at end of file diff --git a/src/Backends/Cubics/UNIFAQLibrary.h b/src/Backends/Cubics/UNIFAQLibrary.h new file mode 100644 index 00000000..c59cb0a3 --- /dev/null +++ b/src/Backends/Cubics/UNIFAQLibrary.h @@ -0,0 +1,108 @@ +#ifndef UNIFAQ_LIBRARY_H +#define UNIFAQ_LIBRARY_H + +#include +#include + +#include "rapidjson_include.h" + +namespace UNIFAQLibrary{ + + /// A structure containing references for a single group (its multiplicity, main group index, etc.) + struct Group{ + int sgi, ///< Sub group index + mgi; ///< Main group index + double R_k, ///< R_k + Q_k; ///< Q_k + }; + + /// A structure containing the parameters for a given mgi-mgi pair + struct InteractionParameters{ + int mgi1, ///< The first main group index + mgi2; ///< The second main group index + double a_ij, ///< + a_ji, ///< + b_ij, ///< + b_ji, ///< + c_ij, ///< + c_ji; ///< + /// Swap a_ij with a_ji, b_ij with b_ji, etc. + void swap() { + std::swap(a_ij, a_ji); + std::swap(b_ij, b_ji); + std::swap(c_ij, c_ji); + } + /// Set all the values to 0 + void zero_out() { + a_ij = 0; a_ji = 0; + b_ij = 0; b_ji = 0; + c_ij = 0; c_ji = 0; + } + }; + + /// A structure containing a group (its count, index, etc.) for a subgroup forming a part of a component + struct ComponentGroup { + int count; + UNIFAQLibrary::Group group; + ComponentGroup(const int count, const UNIFAQLibrary::Group group) : count(count), group(group) {}; + }; + + /// A structure containing the groups and additional information for a component + struct Component{ + std::string name, ///< A user-readable name (not guaranteed unique) + inchikey, ///< The InChI key for the component + registry_number, ///< The registry number for the component in xxxxxxxxxx-xx-x format + userid; ///< A user-specified string identifier + double Tc, ///< The critical temperature in K + pc, ///< The critical pressure in Pa + acentric; ///< The acentric factor + std::vector groups; + }; + + /** + * \brief A container for the parameters for a given UNIFAQ model + * + * This container is intended to be sufficiently generic to allow the user to populate it with UNIFAQ parameters from + * any of the standard UNIFAQ models + * + * Input of parameters (population) is done using JSON-formatted strings, and the class can be interrogated to return + * the desired group information and/or interaction parameters + */ + struct UNIFAQParameterLibrary{ + private: + bool m_populated; ///< True if the library has been populated + std::vector groups; ///< The collection of groups forming the component from the group decomposition + std::vector interaction_parameters; ///< The collection of interaction parameters between main groups in the library + std::vector components; ///< The collection of components that are included in this library + + /// Convert string to JSON document + void jsonize(std::string &s, rapidjson::Document &doc); + + /// Populate internal data structures based on rapidjson Documents + void populate(rapidjson::Value &group_data, rapidjson::Value &interaction_data, rapidjson::Value &decomp_data); + + public: + UNIFAQParameterLibrary() : m_populated(false) {}; + + /// Return true if library has been populated + bool is_populated(){ return m_populated; }; + + /// Populate internal data structures based on JSON-formatted strings + void populate(std::string &group_data, std::string &interaction_data, std::string &decomp_data); + + /// Get the data for group with given sub group index + Group get_group(int sgi) const; + + /// Check if the sub group index can be retrieved + bool has_group(int sgi) const; + + /// Get the group decomposition for a given component + Component get_component(const std::string &identifier, const std::string &value) const; + + /// Get the interaction parameters for given mgi-mgi pair + InteractionParameters get_interaction_parameters(int mgi1, int mgi2) const; + }; + +}; /* namespace UNIFAQLibrary*/ + +#endif \ No newline at end of file diff --git a/src/Backends/Cubics/VTPRBackend.cpp b/src/Backends/Cubics/VTPRBackend.cpp new file mode 100644 index 00000000..ea454037 --- /dev/null +++ b/src/Backends/Cubics/VTPRBackend.cpp @@ -0,0 +1,79 @@ + +#include +#include +#include +#include +#include + +#include "VTPRBackend.h" +#include "Configuration.h" +#include "Exceptions.h" + +static UNIFAQLibrary::UNIFAQParameterLibrary lib; + +void CoolProp::VTPRBackend::setup(const std::vector &names, bool generate_SatL_and_SatV){ + + R = get_config_double(R_U_CODATA); + + // Reset the residual Helmholtz energy class + residual_helmholtz.reset(new CubicResidualHelmholtz(this)); + + // If pure, set the mole fractions to be unity + if (is_pure_or_pseudopure){ + mole_fractions = std::vector(1, 1.0); + mole_fractions_double = std::vector(1, 1.0); + } + + // Now set the reducing function for the mixture + Reducing.reset(new ConstantReducingFunction(cubic->T_r, cubic->rho_r)); + + VTPRCubic * _cubic= static_cast(cubic.get()); + _cubic->get_unifaq().set_components("name", names); + _cubic->get_unifaq().set_interaction_parameters(); + + // Store the fluid names + m_fluid_names = names; + + // Resize the vectors + resize(names.size()); + + // Top-level class can hold copies of the base saturation classes, + // saturation classes cannot hold copies of the saturation classes + if (generate_SatL_and_SatV) + { + bool SatLSatV = false; + SatL.reset(this->get_copy(SatLSatV)); + SatL->specify_phase(iphase_liquid); + linked_states.push_back(SatL); + SatV.reset(this->get_copy(SatLSatV)); + SatV->specify_phase(iphase_gas); + linked_states.push_back(SatV); + + if (is_pure_or_pseudopure) { + std::vector z(1, 1.0); + set_mole_fractions(z); + SatL->set_mole_fractions(z); + SatV->set_mole_fractions(z); + } + } +} + +const UNIFAQLibrary::UNIFAQParameterLibrary & CoolProp::VTPRBackend::LoadLibrary(){ + if (!lib.is_populated()){ + std::string UNIFAQ_path = get_config_string(VTPR_UNIFAQ_PATH); + if (UNIFAQ_path.empty()){ + throw ValueError("You must provide the path to the UNIFAQ library files as VTPR_UNIFAQ_PATH"); + } + if (!(UNIFAQ_path[UNIFAQ_path.size()-1] == '\\' || UNIFAQ_path[UNIFAQ_path.size()-1] == '/')){ + throw ValueError("VTPR_UNIFAQ_PATH must end with / or \\ character"); + } + std::string group_path = UNIFAQ_path + "/group_data.json"; + std::string groups = get_file_contents(group_path.c_str()); + std::string interaction_path = UNIFAQ_path + "/interaction_parameters.json"; + std::string interaction = get_file_contents(interaction_path.c_str()); + std::string decomps_path = UNIFAQ_path + "/decompositions.json"; + std::string decomps = get_file_contents(decomps_path.c_str()); + lib.populate(groups, interaction, decomps); + } + return lib; +} diff --git a/src/Backends/Cubics/VTPRBackend.h b/src/Backends/Cubics/VTPRBackend.h new file mode 100644 index 00000000..69170134 --- /dev/null +++ b/src/Backends/Cubics/VTPRBackend.h @@ -0,0 +1,79 @@ +// +// VTPRBackend.h +// CoolProp +// +// Created by Ian on 7/17/16. +// +// + +#ifndef VTPRBackend_h +#define VTPRBackend_h + +#include "CoolPropTools.h" +#include "DataStructures.h" +#include "GeneralizedCubic.h" +#include "AbstractState.h" +#include "Backends/Helmholtz/HelmholtzEOSMixtureBackend.h" +#include "Exceptions.h" +#include +#include "CubicBackend.h" +#include "UNIFAQLibrary.h" +#include "UNIFAQ.h" +#include "VTPRCubic.h" + +namespace CoolProp { + + +class VTPRBackend : public PengRobinsonBackend { + +private: + + std::vector Tc, pc, omega, m_ii; + double R; + std::vector m_fluid_names; +public: + + VTPRBackend(const std::vector fluid_identifiers, + const double R_u, + bool generate_SatL_and_SatV = true) + : PengRobinsonBackend(fluid_identifiers, R_u) + { + std::vector Tc, pc, acentric; + N = fluid_identifiers.size(); + components.resize(N); + for (std::size_t i = 0; i < fluid_identifiers.size(); ++i){ + components[i] = CubicLibrary::get_cubic_values(fluid_identifiers[i]); + Tc.push_back(components[i].Tc); + pc.push_back(components[i].pc); + acentric.push_back(components[i].acentric); + } + cubic.reset(new VTPRCubic(Tc, pc, acentric, R_u, LoadLibrary())); + setup(fluid_identifiers, generate_SatL_and_SatV); + }; + AbstractCubicBackend * get_copy(bool generate_SatL_and_SatV){ + AbstractCubicBackend * ACB = new VTPRBackend(fluid_names(), cubic->get_R_u(), generate_SatL_and_SatV); + //ACB->copy_k(this); + return ACB; + } + + /// Return the fluid names + std::vector calc_fluid_names(void) { return m_fluid_names; } + + /// Set the pointer to the residual helmholtz class, etc. + void setup(const std::vector &names, bool generate_SatL_and_SatV = true); + + /// Load the UNIFAQ library if needed and get const reference to it + const UNIFAQLibrary::UNIFAQParameterLibrary &LoadLibrary(); + + void set_mole_fractions(const std::vector &z){ + mole_fractions = z; + mole_fractions_double = z; + VTPRCubic * _cubic= static_cast(cubic.get()); + _cubic->get_unifaq().set_mole_fractions(z); + }; + +}; + +}; /* namespace CoolProp */ + +#endif /* VTPRBackend_h */ diff --git a/src/Backends/Cubics/VTPRCubic.h b/src/Backends/Cubics/VTPRCubic.h new file mode 100644 index 00000000..46a8f283 --- /dev/null +++ b/src/Backends/Cubics/VTPRCubic.h @@ -0,0 +1,166 @@ +// +// VTPRCubic.h +// CoolProp +// +// Created by Ian on 7/17/16. +// +// + +#include "GeneralizedCubic.h" + +#ifndef VTPRCubic_h +#define VTPRCubic_h + +class VTPRCubic : public PengRobinson +{ +private: + UNIFAQ::UNIFAQMixture unifaq; +public: + VTPRCubic(std::vector Tc, + std::vector pc, + std::vector acentric, + double R_u, + const UNIFAQLibrary::UNIFAQParameterLibrary & lib + ) + : PengRobinson(Tc, pc, acentric, R_u), unifaq(lib) {}; + + VTPRCubic(double Tc, + double pc, + double acentric, + double R_u, + const UNIFAQLibrary::UNIFAQParameterLibrary & lib) + : PengRobinson(std::vector(1, Tc), std::vector(1, pc), std::vector(1, acentric), R_u), unifaq(lib) {}; + + /// Get a reference to the managed UNIFAQ instance + UNIFAQ::UNIFAQMixture &get_unifaq() { return unifaq; } + + /// The attractive part in cubic EOS + double a_alpha(double T, std::size_t i) { + return pow(1 + m_ii(i)*(1 - sqrt(T / Tc[i])), 2); + } + /// Calculate the non-dimensionalized gE/RT term + double gE_R_RT() { + const std::vector &z = unifaq.get_mole_fractions(); + double summer = 0; + for (std::size_t i = 0; i < z.size(); ++i) { + summer += z[i] * unifaq.ln_gamma_R(i); + } + return summer; + } + double d_gE_R_RT_dxi(const std::vector &x, std::size_t i, bool xN_independent) { + if (xN_independent) + { + return unifaq.ln_gamma_R(i); + } + else { + return unifaq.ln_gamma_R(i) - unifaq.ln_gamma_R(N - 1); + } + } + /// The co-volume for the i-th pure component + double b_ii(std::size_t i) { + return 0.0778*R_u*Tc[i] / pc[i]; + } + /// The attractive parameter for the i-th pure component + double a_ii(std::size_t i) { + double a0 = 0.45724*pow(R_u*Tc[i], 2) / pc[i]; + double alpha = a_alpha(unifaq.get_temperature(), i); + return a0*alpha; + } + double am_term(double tau, const std::vector &x, std::size_t itau) { + if (itau == 0) { + set_temperature(T_r / tau); + return bm_term(x)*(sum_xi_aii_bii(x) + R_u*unifaq.get_temperature()*gE_R_RT() / (-0.53087)); + } + else { + double dtau = 0.01*tau; + return (am_term(tau + dtau, x, itau - 1) - am_term(tau - dtau, x, itau - 1)) / (2 * dtau); + } + } + double sum_xi_aii_bii(const std::vector &x) { + double summeram = 0; + for (std::size_t i = 0; i < N; ++i) { + summeram += x[i] * a_ii(i) / b_ii(i); + } + return summeram; + } + double d_sum_xi_aii_bii_dxi(const std::vector &x, std::size_t i, bool xN_independent) { + if (xN_independent) + { + return x[i] * a_ii(i) / b_ii(i); + } + else { + return x[i] * a_ii(i) / b_ii(i) - x[N - 1] * a_ii(N - 1) / b_ii(N - 1); + } + } + double d_am_term_dxi(double tau, const std::vector &x, std::size_t itau, std::size_t i, bool xN_independent) + { + set_temperature(T_r / tau); + return d_bm_term_dxi(x, i, xN_independent)*(sum_xi_aii_bii(x) + R_u*unifaq.get_temperature()*gE_R_RT() / (-0.53087)) + + bm_term(x)*(d_sum_xi_aii_bii_dxi(x, i, xN_independent) + R_u*unifaq.get_temperature()*d_gE_R_RT_dxi(x, i, xN_independent) / (-0.53087)); + } + double d2_am_term_dxidxj(double tau, const std::vector &x, std::size_t itau, std::size_t i, std::size_t j, bool xN_independent) + { + set_temperature(T_r / tau); + return d2_bm_term_dxidxj(x, i, j, xN_independent)*(sum_xi_aii_bii(x) + R_u*unifaq.get_temperature()*gE_R_RT() / (-0.53087)) + + d_bm_term_dxi(x, i, xN_independent)*(d_sum_xi_aii_bii_dxi(x, j, xN_independent) + R_u*unifaq.get_temperature()*d_gE_R_RT_dxi(x, j, xN_independent) / (-0.53087)) + + d_bm_term_dxi(x, j, xN_independent)*(d_sum_xi_aii_bii_dxi(x, i, xN_independent) + R_u*unifaq.get_temperature()*d_gE_R_RT_dxi(x, i, xN_independent) / (-0.53087)); + } + double d3_am_term_dxidxjdxk(double tau, const std::vector &x, std::size_t itau, std::size_t i, std::size_t j, std::size_t k, bool xN_independent) + { + set_temperature(T_r / tau); + return d3_bm_term_dxidxjdxk(x, i, j, k, xN_independent)*(sum_xi_aii_bii(x) + R_u*unifaq.get_temperature()*gE_R_RT() / (-0.53087)) + + d2_bm_term_dxidxj(x, i, k, xN_independent)*(d_sum_xi_aii_bii_dxi(x, j, xN_independent) + R_u*unifaq.get_temperature()*d_gE_R_RT_dxi(x, j, xN_independent) / (-0.53087)) + + d2_bm_term_dxidxj(x, j, k, xN_independent)*(d_sum_xi_aii_bii_dxi(x, i, xN_independent) + R_u*unifaq.get_temperature()*d_gE_R_RT_dxi(x, i, xN_independent) / (-0.53087)); + } + + double bm_term(const std::vector &x) { + double summerbm = 0; + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < N; ++j) { + summerbm += x[i] * x[j] * bij_term(i, j); + } + } + return summerbm; + } + double bij_term(std::size_t i, std::size_t j) + { + return pow((pow(b_ii(i), 0.75) + pow(b_ii(j), 0.75)) / 2.0, 4.0 / 3.0); + } + double d_bm_term_dxi(const std::vector &x, std::size_t i, bool xN_independent) + { + double summer = 0; + if (xN_independent) + { + for (int j = N - 1; j >= 0; --j) + { + summer += x[j] * bij_term(i, j); + } + return 2 * summer; + } + else { + for (int k = N - 2; k >= 0; --k) + { + summer += x[k] * (bij_term(i, k) - bij_term(k, N - 1)); + } + return 2 * (summer + x[N - 1] * (bij_term(N - 1, i) - bij_term(N - 1, N - 1))); + } + } + double d2_bm_term_dxidxj(const std::vector &x, std::size_t i, std::size_t j, bool xN_independent) + { + if (xN_independent) + { + return 2 * bij_term(i, j); + } + else { + return 2 * (bij_term(i, j) - bij_term(j, N - 1) - bij_term(N - 1, i) + bij_term(N - 1, N - 1)); + } + } + double d3_bm_term_dxidxjdxk(const std::vector &x, std::size_t i, std::size_t j, std::size_t k, bool xN_independent) + { + return 0; + } + + void set_temperature(const double T) { unifaq.set_temperature(T); } +}; + +#endif /* VTPRCubic_h */ diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index b2089ef3..c219686d 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -149,6 +149,7 @@ void HelmholtzEOSMixtureBackend::set_mass_fractions(const std::vectormole_fractions.resize(N); + this->mole_fractions_double.resize(N); this->K.resize(N); this->lnK.resize(N); for (std::vector >::iterator it = linked_states.begin(); it != linked_states.end(); ++it) { diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index 6553a0f5..0f206317 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -100,6 +100,9 @@ public: /// Apply a simple mixing rule void apply_simple_mixing_rule(std::size_t i, std::size_t j, const std::string &model); + // Set fluid parameter (currently the volume translation parameter for cubic) + virtual void set_fluid_parameter_double(const size_t i, const std::string parameter, const double value) { throw ValueError("set_fluid_parameter_double only defined for cubic backends"); }; + phases calc_phase(void){return _phase;}; /** \brief Specify the phase - this phase will always be used in calculations diff --git a/src/CoolPropLib.cpp b/src/CoolPropLib.cpp index c4c4001b..700c9feb 100644 --- a/src/CoolPropLib.cpp +++ b/src/CoolPropLib.cpp @@ -781,6 +781,36 @@ EXPORT_CODE void CONVENTION AbstractState_set_binary_interaction_double(const lo } } +EXPORT_CODE void CONVENTION AbstractState_set_fluid_parameter_double(const long handle, const size_t i, const char* parameter, const double value , long *errcode, char *message_buffer, const long buffer_length) { + *errcode = 0; + try { + shared_ptr &AS = handle_manager.get(handle); + AS->set_fluid_parameter_double(static_cast(i), parameter, value); + } + catch (CoolProp::HandleError &e) { + std::string errmsg = std::string("HandleError: ") + e.what(); + if (errmsg.size() < static_cast(buffer_length)) { + *errcode = 1; + strcpy(message_buffer, errmsg.c_str()); + } + else { + *errcode = 2; + } + } + catch (CoolProp::CoolPropBaseError &e) { + std::string errmsg = std::string("Error: ") + e.what(); + if (errmsg.size() < static_cast(buffer_length)) { + *errcode = 1; + strcpy(message_buffer, errmsg.c_str()); + } + else { + *errcode = 2; + } + } + catch (...) { + *errcode = 3; + } +} /// ********************************************************************************* /// ********************************************************************************* diff --git a/src/CoolPropLib.def b/src/CoolPropLib.def index f3bd523c..89725b42 100644 --- a/src/CoolPropLib.def +++ b/src/CoolPropLib.def @@ -6,6 +6,7 @@ EXPORTS AbstractState_keyed_output = _AbstractState_keyed_output@20 AbstractState_set_binary_interaction_double = AbstractState_set_binary_interaction_double@36 AbstractState_set_fractions = _AbstractState_set_fractions@24 + AbstractState_set_volume_translation = _AbstractState_set_volume_translation@24 AbstractState_update = _AbstractState_update@36 AbstractState_first_saturation_deriv = _AbstractState_first_saturation_deriv@24 AbstractState_first_partial_deriv = _AbstractState_first_partial_deriv@28 @@ -36,6 +37,7 @@ EXPORTS propssi_ = _propssi_@28 redirect_stdout = _redirect_stdout@4 saturation_ancillary = _saturation_ancillary@24 + set_config_string = _set_config_string@8 set_debug_level = _set_debug_level@4 set_reference_stateD = _set_reference_stateD@36 set_reference_stateS = _set_reference_stateS@8 diff --git a/src/DataStructures.cpp b/src/DataStructures.cpp index b92a60ec..8a31eb5d 100644 --- a/src/DataStructures.cpp +++ b/src/DataStructures.cpp @@ -550,7 +550,8 @@ const backend_family_info backend_family_list[] = { { TTSE_BACKEND_FAMILY, "TTSE" }, { BICUBIC_BACKEND_FAMILY, "BICUBIC" }, { SRK_BACKEND_FAMILY, "SRK" }, - { PR_BACKEND_FAMILY, "PR" } + { PR_BACKEND_FAMILY, "PR" }, + { VTPR_BACKEND_FAMILY, "VTPR" } }; const backend_info backend_list[] = { @@ -564,7 +565,8 @@ const backend_info backend_list[] = { { TTSE_BACKEND, "TTSEBackend", TTSE_BACKEND_FAMILY }, { BICUBIC_BACKEND, "BicubicBackend", BICUBIC_BACKEND_FAMILY }, { SRK_BACKEND, "SRKBackend", SRK_BACKEND_FAMILY }, - { PR_BACKEND, "PengRobinsonBackend", PR_BACKEND_FAMILY } + { PR_BACKEND, "PengRobinsonBackend", PR_BACKEND_FAMILY }, + { VTPR_BACKEND, "VTPRBackend", VTPR_BACKEND_FAMILY } }; class BackendInformation { diff --git a/wrappers/Julia/CoolProp.jl b/wrappers/Julia/CoolProp.jl index f968921a..7f53f9fd 100644 --- a/wrappers/Julia/CoolProp.jl +++ b/wrappers/Julia/CoolProp.jl @@ -1,7 +1,7 @@ VERSION < v"0.4.0" && __precompile__() module CoolProp -export PropsSI, PhaseSI, get_global_param_string, get_parameter_information_string,get_fluid_param_string,set_reference_stateS, get_param_index, get_input_pair_index, set_config_string, F2K, K2F, HAPropsSI, AbstractState_factory, AbstractState_free, AbstractState_set_fractions, AbstractState_update, AbstractState_specify_phase, AbstractState_unspecify_phase, AbstractState_keyed_output, AbstractState_output, AbstractState_update_and_common_out, AbstractState_update_and_1_out, AbstractState_update_and_5_out, AbstractState_set_binary_interaction_double +export PropsSI, PhaseSI, get_global_param_string, get_parameter_information_string,get_fluid_param_string,set_reference_stateS, get_param_index, get_input_pair_index, set_config_string, F2K, K2F, HAPropsSI, AbstractState_factory, AbstractState_free, AbstractState_set_fractions, AbstractState_update, AbstractState_specify_phase, AbstractState_unspecify_phase, AbstractState_keyed_output, AbstractState_output, AbstractState_update_and_common_out, AbstractState_update_and_1_out, AbstractState_update_and_5_out, AbstractState_set_binary_interaction_double, AbstractState_set_fluid_parameter_double # Check the current Julia version to make this Julia 0.4 code compatible with older version if VERSION <= VersionNumber(0,4) @@ -401,4 +401,23 @@ function AbstractState_set_binary_interaction_double(handle::Clong,i::Int, j::In return nothing end +# Set some fluid parameter (ie volume translation for cubic) +# handle The integer handle for the state class stored in memory +# i indice of the fluid the parramter should be applied too (for mixtures) +# parameter string wit the name of the parameter +# value the value of the parameter +function AbstractState_set_fluid_parameter_double(handle::Clong, i::Int, parameter::AbstractString, value::Cdouble) + ccall( (:AbstractState_set_fluid_parameter_double, "CoolProp"), Void, (Clong,Clong,Ptr{UInt8},Cdouble,Ref{Clong},Ptr{UInt8},Clong), handle,i,parameter,value,errcode,message_buffer::Array{UInt8,1},buffer_length) + if errcode[] != 0 + if errcode[] == 1 + error("CoolProp: ", bytestring(convert(Ptr{UInt8}, pointer(message_buffer)))) + elseif errcode[] == 2 + error("CoolProp: message buffer too small") + else # == 3 + error("CoolProp: unknown error") + end + end + return nothing +end + end #module