From 111371df9857d4bd151da7d27ec7b4ce15192546 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Tue, 30 Jun 2015 22:04:08 -0600 Subject: [PATCH] Implement fugacity through low-level interface; closes #699 --- externals/REFPROP-headers | 2 +- include/AbstractState.h | 9 ++++++-- src/AbstractState.cpp | 6 ++++- .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 7 +++++- .../Helmholtz/HelmholtzEOSMixtureBackend.h | 5 ++++- .../REFPROP/REFPROPMixtureBackend.cpp | 22 ++++++++++++++++--- src/Backends/REFPROP/REFPROPMixtureBackend.h | 3 ++- wrappers/Python/CoolProp/AbstractState.pxd | 2 ++ wrappers/Python/CoolProp/AbstractState.pyx | 6 +++++ wrappers/Python/CoolProp/cAbstractState.pxd | 2 ++ 10 files changed, 54 insertions(+), 10 deletions(-) diff --git a/externals/REFPROP-headers b/externals/REFPROP-headers index f452de52..62a12441 160000 --- a/externals/REFPROP-headers +++ b/externals/REFPROP-headers @@ -1 +1 @@ -Subproject commit f452de5261c59f49f1166cd6133cc08efc49b46e +Subproject commit 62a12441151511e1bcf19b1d2c310c96195ec861 diff --git a/include/AbstractState.h b/include/AbstractState.h index 921bb357..9ec1bf25 100644 --- a/include/AbstractState.h +++ b/include/AbstractState.h @@ -144,7 +144,9 @@ protected: /// Using this backend, calculate the universal gas constant \f$R_u\f$ in J/mol/K virtual CoolPropDbl calc_gas_constant(void){ throw NotImplementedError("calc_gas_constant is not implemented for this backend"); }; /// Using this backend, calculate the fugacity coefficient (dimensionless) - virtual CoolPropDbl calc_fugacity_coefficient(int i){ throw NotImplementedError("calc_fugacity_coefficient is not implemented for this backend"); }; + virtual CoolPropDbl calc_fugacity_coefficient(std::size_t i){ throw NotImplementedError("calc_fugacity_coefficient is not implemented for this backend"); }; + /// Using this backend, calculate the fugacity in Pa + virtual CoolPropDbl calc_fugacity(std::size_t i){ throw NotImplementedError("calc_fugacity is not implemented for this backend"); }; /// Using this backend, calculate the phase identification parameter (PIP) virtual CoolPropDbl calc_PIP(void){ throw NotImplementedError("calc_PIP is not implemented for this backend"); }; @@ -544,7 +546,10 @@ public: double isothermal_compressibility(void); /// Return the isobaric expansion coefficient \f$ \beta = \frac{1}{v}\left.\frac{\partial v}{\partial T}\right|_p = -\frac{1}{\rho}\left.\frac{\partial \rho}{\partial T}\right|_p\f$ in 1/K double isobaric_expansion_coefficient(void); - double fugacity_coefficient(int i); + /// Return the fugacity coefficient of the i-th component of the mixture + double fugacity_coefficient(std::size_t i); + /// Return the fugacity of the i-th component of the mixture + double fugacity(std::size_t i); /// Return the fundamental derivative of gas dynamics //double fundamental_derivative_of_gas_dynamics(void){return this->second_partial_deriv(iP, iDmolar, iSmolar, iDmolar, iSmolar)/pow(speed_sound(), 2)/2/pow(this->rhomolar(),3);}; /// Return the phase identification parameter (PIP) of G. Venkatarathnam and L.R. Oellrich, "Identification of the phase of a fluid using partial derivatives of pressure, volume, and temperature without reference to saturation properties: Applications in phase equilibria calculations" diff --git a/src/AbstractState.cpp b/src/AbstractState.cpp index e2d88cb4..a88789cd 100644 --- a/src/AbstractState.cpp +++ b/src/AbstractState.cpp @@ -431,10 +431,14 @@ double AbstractState::gas_constant(void){ if (!_gas_constant) _gas_constant = calc_gas_constant(); return _gas_constant; } -double AbstractState::fugacity_coefficient(int i){ +double AbstractState::fugacity_coefficient(std::size_t i){ // TODO: Cache the fug. coeff for each component return calc_fugacity_coefficient(i); } +double AbstractState::fugacity(std::size_t i){ + // TODO: Cache the fug. coeff for each component + return calc_fugacity(i); +} void AbstractState::build_phase_envelope(const std::string &type) { calc_phase_envelope(type); diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index cddb95a2..3ac82913 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -2341,11 +2341,16 @@ CoolPropDbl HelmholtzEOSMixtureBackend::calc_gibbsmolar(void) throw ValueError(format("phase is invalid in calc_gibbsmolar")); } } -CoolPropDbl HelmholtzEOSMixtureBackend::calc_fugacity_coefficient(int i) +CoolPropDbl HelmholtzEOSMixtureBackend::calc_fugacity_coefficient(std::size_t i) { x_N_dependency_flag xN_flag = XN_DEPENDENT; return exp(MixtureDerivatives::ln_fugacity_coefficient(*this, i, xN_flag)); } +CoolPropDbl HelmholtzEOSMixtureBackend::calc_fugacity(std::size_t i) +{ + x_N_dependency_flag xN_flag = XN_DEPENDENT; + return MixtureDerivatives::fugacity_i(*this, i, xN_flag); +} CoolPropDbl HelmholtzEOSMixtureBackend::calc_phase_identification_parameter(void) { return 2 - rhomolar()*(second_partial_deriv(iP, iDmolar, iT, iT, iDmolar)/first_partial_deriv(iP, iT, iDmolar) - second_partial_deriv(iP, iDmolar, iT, iDmolar, iT)/first_partial_deriv(iP, iDmolar, iT)); diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index 885aa2b1..d7dfb032 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -192,8 +192,11 @@ public: * * */ - CoolPropDbl calc_fugacity_coefficient(int i); + + CoolPropDbl calc_phase_identification_parameter(void); + CoolPropDbl calc_fugacity(std::size_t i); + CoolPropDbl calc_fugacity_coefficient(std::size_t i); /// Using this backend, calculate the flame hazard CoolPropDbl calc_flame_hazard(void){ return components[0].environment.FH;}; diff --git a/src/Backends/REFPROP/REFPROPMixtureBackend.cpp b/src/Backends/REFPROP/REFPROPMixtureBackend.cpp index c6e9a362..0a4f1b5c 100644 --- a/src/Backends/REFPROP/REFPROPMixtureBackend.cpp +++ b/src/Backends/REFPROP/REFPROPMixtureBackend.cpp @@ -22,8 +22,10 @@ surface tension N/m #define _CRT_SECURE_NO_WARNINGS #define REFPROP_IMPLEMENTATION +#define REFPROP_CSTYLE_REFERENCES #include "REFPROP_lib.h" #undef REFPROP_IMPLEMENTATION +#undef REFPROP_CSTYLE_REFERENCES #include "CoolPropTools.h" #include "REFPROPMixtureBackend.h" @@ -584,7 +586,7 @@ CoolPropDbl REFPROPMixtureBackend::calc_surface_tension(void) _surface_tension = sigma; return static_cast(_surface_tension); } -CoolPropDbl REFPROPMixtureBackend::calc_fugacity_coefficient(int i) +CoolPropDbl REFPROPMixtureBackend::calc_fugacity_coefficient(std::size_t i) { this->check_loaded_fluid(); double rho_mol_L = 0.001*_rhomolar; @@ -593,12 +595,26 @@ CoolPropDbl REFPROPMixtureBackend::calc_fugacity_coefficient(int i) fug_cof.resize(mole_fractions.size()); char herr[255]; FUGCOFdll(&_T, &rho_mol_L, &(mole_fractions[0]), // Inputs - &(fug_cof[0]), // Outputs - &ierr, herr, errormessagelength); // Error message + &(fug_cof[0]), // Outputs + &ierr, herr, errormessagelength); // Error message if (static_cast(ierr) > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());} return static_cast(fug_cof[i]); } +CoolPropDbl REFPROPMixtureBackend::calc_fugacity(std::size_t i) +{ + this->check_loaded_fluid(); + double rho_mol_L = 0.001*_rhomolar; + long ierr = 0; + std::vector f(mole_fractions.size()); + char herr[255]; + FGCTY2dll(&_T, &rho_mol_L, &(mole_fractions[0]), // Inputs + &(f[0]), // Outputs + &ierr, herr, errormessagelength); // Error message + if (static_cast(ierr) > 0) { throw ValueError(format("%s", herr).c_str()); } + //else if (ierr < 0) {set_warning(format("%s",herr).c_str());} + return static_cast(f[i]*1000); +} void REFPROPMixtureBackend::calc_phase_envelope(const std::string &type) { diff --git a/src/Backends/REFPROP/REFPROPMixtureBackend.h b/src/Backends/REFPROP/REFPROPMixtureBackend.h index a0b5ae6c..a220cc36 100644 --- a/src/Backends/REFPROP/REFPROPMixtureBackend.h +++ b/src/Backends/REFPROP/REFPROPMixtureBackend.h @@ -130,7 +130,8 @@ public: /// Calc the C virial coefficient CoolPropDbl calc_Cvirial(void); - CoolPropDbl calc_fugacity_coefficient(int i); + CoolPropDbl calc_fugacity_coefficient(std::size_t i); + CoolPropDbl calc_fugacity(std::size_t i); CoolPropDbl calc_melting_line(int param, int given, CoolPropDbl value); bool has_melting_line(){return true;}; double calc_melt_Tmax(); diff --git a/wrappers/Python/CoolProp/AbstractState.pxd b/wrappers/Python/CoolProp/AbstractState.pxd index d17ea904..995e9093 100644 --- a/wrappers/Python/CoolProp/AbstractState.pxd +++ b/wrappers/Python/CoolProp/AbstractState.pxd @@ -68,6 +68,8 @@ cdef class AbstractState: cpdef double Bvirial(self) except * cpdef double Cvirial(self) except * cpdef double PIP(self) except * + cpdef double fugacity(self, size_t) except * + cpdef double fugacity_coefficient(self, size_t) except * cpdef double molar_mass(self) except * cpdef double acentric_factor(self) except* diff --git a/wrappers/Python/CoolProp/AbstractState.pyx b/wrappers/Python/CoolProp/AbstractState.pyx index 9c583dc3..f2760970 100644 --- a/wrappers/Python/CoolProp/AbstractState.pyx +++ b/wrappers/Python/CoolProp/AbstractState.pyx @@ -182,6 +182,12 @@ cdef class AbstractState: cpdef double PIP(self) except *: """ Get the phase identification parameter - wrapper of c++ function :cpapi:`CoolProp::AbstractState::PIP(void)` """ return self.thisptr.PIP() + cpdef double fugacity(self, size_t i) except *: + """ Get the fugacity of the i-th component - wrapper of c++ function :cpapi:`CoolProp::AbstractState::fugacity(std::size_t)` """ + return self.thisptr.fugacity(i) + cpdef double fugacity_coefficient(self, size_t i) except *: + """ Get the fugacity coefficient of the i-th component - wrapper of c++ function :cpapi:`CoolProp::AbstractState::fugacity_coefficient(std::size_t)` """ + return self.thisptr.fugacity_coefficient(i) cpdef mole_fractions_liquid(self): """ Get the mole fractions of the liquid phase - wrapper of c++ function :cpapi:`CoolProp::AbstractState::mole_fractions_liquid(void)` """ diff --git a/wrappers/Python/CoolProp/cAbstractState.pxd b/wrappers/Python/CoolProp/cAbstractState.pxd index 7ebf988d..795fbb02 100644 --- a/wrappers/Python/CoolProp/cAbstractState.pxd +++ b/wrappers/Python/CoolProp/cAbstractState.pxd @@ -85,6 +85,8 @@ cdef extern from "AbstractState.h" namespace "CoolProp": double Bvirial() except +ValueError double Cvirial() except +ValueError double PIP() except +ValueError + double fugacity(size_t) except +ValueError + double fugacity_coefficient(size_t) except +ValueError double keyed_output(constants_header.parameters) except+ValueError double trivial_keyed_output(constants_header.parameters) except+ValueError