From 57584a7254091d714e72ff8f5d5da74cf3e87165 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Sun, 21 Aug 2016 12:34:55 -0600 Subject: [PATCH] Merge VTPR (#1195) * Added first working version of VTPR; * Get VTPR building hopefully * Remove more constant iterators * Make VTPR accessible through factory * One last const iterator * Fix a_alpha bug and make sqrt(2) into sqrt(2.0) * Added analytic first derivative for VTPR * Fix another set of sqrt(2) -> sqrt(2.0) * Add some info on the derivatives for VTPR Looks less hopeless than I had feared * gE/RT needs to be multiplied by RT; closes #1161 * Added first working version of VTPR; * Get VTPR building hopefully * Remove more constant iterators * Make VTPR accessible through factory * One last const iterator * Fix a_alpha bug and make sqrt(2) into sqrt(2.0) * Added analytic first derivative for VTPR * Fix another set of sqrt(2) -> sqrt(2.0) * Add some info on the derivatives for VTPR Looks less hopeless than I had feared * gE/RT needs to be multiplied by RT; closes #1161 * Add VTPR code from @JonWel (#1194) * 1rst draft to implement a simple volume translation to cubics * A bit more of VT * Derivatives for volume translation * Better cm initialisation * Solves the cubic equation with volume translation * Correct the volume translation analytic development Looks good now * Update VTPR to be able to use volume translation * Unprotect cm_term This allows it to be used from the VTPR backend * Update CoolPropLib.def * Better derrivative of PI_12 The expression is simpler this way * Solves #1176 Thanks @ibell * Change the way the volume translation parrameter is set * Start the bm derivatives for VTPR * Correct one derivative * Small bug * Better bm derivatives for VTPR * Add am and bm component derivatives for VTPR @ibell I did not check yet the component derivatives of this commit, bu I checked the other ones with your code. I'll have to addapt your code to also check these ones. I separate the `am_term` and `bm_term` as the `am_bm_term` function was called twice. This reduce the call to the am_term part as this part ends up being called only once, and this helped writing the component derivatives. The tau derivative is done numerically untill we find time to develop the analytical one. The `am_bm_term` function started with a `set_temperature()`. I did not checked yet why this is needed and put this set temperature at the beginning of each of the `am_term` component derivatives. I'll try to addapt the checking code tomorow. * tab to spaces * Re-writing of cubic coefficients Introducing 3 intermediary varriables that simplify the cubic's coefficient with the volume translation. * 1rst draft to implement a simple volume translation to cubics * A bit more of VT * Derivatives for volume translation * Better cm initialisation * Solves the cubic equation with volume translation * Correct the volume translation analytic development Looks good now * Update VTPR to be able to use volume translation * Unprotect cm_term This allows it to be used from the VTPR backend * Update CoolPropLib.def * Better derrivative of PI_12 The expression is simpler this way * Solves #1176 Thanks @ibell * Change the way the volume translation parrameter is set * Start the bm derivatives for VTPR * Correct one derivative * Small bug * Better bm derivatives for VTPR * Add am and bm component derivatives for VTPR @ibell I did not check yet the component derivatives of this commit, bu I checked the other ones with your code. I'll have to addapt your code to also check these ones. I separate the `am_term` and `bm_term` as the `am_bm_term` function was called twice. This reduce the call to the am_term part as this part ends up being called only once, and this helped writing the component derivatives. The tau derivative is done numerically untill we find time to develop the analytical one. The `am_bm_term` function started with a `set_temperature()`. I did not checked yet why this is needed and put this set temperature at the beginning of each of the `am_term` component derivatives. I'll try to addapt the checking code tomorow. * tab to spaces * Re-writing of cubic coefficients Introducing 3 intermediary varriables that simplify the cubic's coefficient with the volume translation. --- ...ympy derivatives of alphar from VTPR.ipynb | 69 ++++++ include/AbstractState.h | 2 + include/Configuration.h | 1 + include/CoolPropLib.h | 13 + include/DataStructures.h | 6 +- src/AbstractState.cpp | 5 + src/Backends/Cubics/CubicBackend.cpp | 55 ++--- src/Backends/Cubics/CubicBackend.h | 3 + src/Backends/Cubics/GeneralizedCubic.cpp | 77 +++--- src/Backends/Cubics/GeneralizedCubic.h | 35 ++- src/Backends/Cubics/UNIFAQ.cpp | 231 ++++++++++++++++++ src/Backends/Cubics/UNIFAQ.h | 94 +++++++ src/Backends/Cubics/UNIFAQLibrary.cpp | 118 +++++++++ src/Backends/Cubics/UNIFAQLibrary.h | 108 ++++++++ src/Backends/Cubics/VTPRBackend.cpp | 79 ++++++ src/Backends/Cubics/VTPRBackend.h | 79 ++++++ src/Backends/Cubics/VTPRCubic.h | 166 +++++++++++++ .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 1 + .../Helmholtz/HelmholtzEOSMixtureBackend.h | 3 + src/CoolPropLib.cpp | 30 +++ src/CoolPropLib.def | 2 + src/DataStructures.cpp | 6 +- wrappers/Julia/CoolProp.jl | 21 +- 23 files changed, 1123 insertions(+), 81 deletions(-) create mode 100644 doc/notebooks/Sympy derivatives of alphar from VTPR.ipynb create mode 100644 src/Backends/Cubics/UNIFAQ.cpp create mode 100644 src/Backends/Cubics/UNIFAQ.h create mode 100644 src/Backends/Cubics/UNIFAQLibrary.cpp create mode 100644 src/Backends/Cubics/UNIFAQLibrary.h create mode 100644 src/Backends/Cubics/VTPRBackend.cpp create mode 100644 src/Backends/Cubics/VTPRBackend.h create mode 100644 src/Backends/Cubics/VTPRCubic.h 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