diff --git a/externals/REFPROP-headers b/externals/REFPROP-headers index 6f5e90be..4c9636ea 160000 --- a/externals/REFPROP-headers +++ b/externals/REFPROP-headers @@ -1 +1 @@ -Subproject commit 6f5e90be7ab5ebceb9f7465e0b1b81e66de0c15d +Subproject commit 4c9636ea56475da750b136b7d15f793ff9d5daea diff --git a/include/AbstractState.h b/include/AbstractState.h index ef021b4c..8d020fcb 100644 --- a/include/AbstractState.h +++ b/include/AbstractState.h @@ -90,6 +90,9 @@ protected: CachedElement _hmolar, _smolar, _umolar, _logp, _logrhomolar, _cpmolar, _cp0molar, _cvmolar, _speed_sound, _gibbsmolar, _helmholtzmolar; + /// Excess properties + CachedElement _hmolar_excess, _smolar_excess, _gibbsmolar_excess, _umolar_excess, _volumemolar_excess, _helmholtzmolar_excess; + /// Ancillary values CachedElement _rhoLanc, _rhoVanc, _pLanc, _pVanc, _TLanc, _TVanc; @@ -161,6 +164,10 @@ protected: /// Using this backend, calculate the phase identification parameter (PIP) virtual CoolPropDbl calc_PIP(void){ throw NotImplementedError("calc_PIP is not implemented for this backend"); }; + // Excess properties + /// Using this backend, calculate and cache the excess properties + virtual void calc_excess_properties(void) { throw NotImplementedError("calc_excess_properties is not implemented for this backend"); }; + // Derivatives of residual helmholtz energy /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r\f$ (dimensionless) virtual CoolPropDbl calc_alphar(void){ throw NotImplementedError("calc_alphar is not implemented for this backend"); }; @@ -288,13 +295,19 @@ protected: /// virtual CoolPropDbl calc_rhomass(void){ return rhomolar()*molar_mass(); } virtual CoolPropDbl calc_hmass(void){ return hmolar() / molar_mass(); } + virtual CoolPropDbl calc_hmass_excess(void) { return hmolar_excess() / molar_mass(); } virtual CoolPropDbl calc_smass(void){ return smolar() / molar_mass(); } + virtual CoolPropDbl calc_smass_excess(void) { return smolar_excess() / molar_mass(); } virtual CoolPropDbl calc_cpmass(void){ return cpmolar() / molar_mass(); } virtual CoolPropDbl calc_cp0mass(void){ return cp0molar() / molar_mass(); } virtual CoolPropDbl calc_cvmass(void){ return cvmolar() / molar_mass(); } virtual CoolPropDbl calc_umass(void){ return umolar() / molar_mass(); } + virtual CoolPropDbl calc_umass_excess(void) { return umolar_excess() / molar_mass(); } virtual CoolPropDbl calc_gibbsmass(void){ return gibbsmolar() / molar_mass(); } + virtual CoolPropDbl calc_gibbsmass_excess(void) { return gibbsmolar_excess() / molar_mass(); } virtual CoolPropDbl calc_helmholtzmass(void){ return helmholtzmolar() / molar_mass(); } + virtual CoolPropDbl calc_helmholtzmass_excess(void) { return helmholtzmolar_excess() / molar_mass(); } + virtual CoolPropDbl calc_volumemass_excess(void) { return volumemolar_excess() / molar_mass(); } /// Update the states after having changed the reference state for enthalpy and entropy virtual void update_states(void){ throw NotImplementedError("This backend does not implement update_states function"); }; @@ -619,14 +632,26 @@ public: double hmolar(void); /// Return the mass enthalpy in J/kg double hmass(void){ return calc_hmass(); }; + /// Return the excess molar enthalpy in J/mol + double hmolar_excess(void); + /// Return the excess mass enthalpy in J/kg + double hmass_excess(void) { return calc_hmass_excess(); }; /// Return the molar entropy in J/mol/K double smolar(void); /// Return the molar entropy in J/kg/K double smass(void){ return calc_smass(); }; + /// Return the molar entropy in J/mol/K + double smolar_excess(void); + /// Return the molar entropy in J/kg/K + double smass_excess(void) { return calc_smass_excess(); }; /// Return the molar internal energy in J/mol double umolar(void); /// Return the mass internal energy in J/kg double umass(void){ return calc_umass(); }; + /// Return the excess internal energy in J/mol + double umolar_excess(void); + /// Return the excess internal energy in J/kg + double umass_excess(void) { return calc_umass_excess(); }; /// Return the molar constant pressure specific heat in J/mol/K double cpmolar(void); /// Return the mass constant pressure specific heat in J/kg/K @@ -639,14 +664,26 @@ public: double cvmolar(void); /// Return the mass constant volume specific heat in J/kg/K double cvmass(void){ return calc_cvmass(); }; - /// Return the Gibbs function in J/mol + /// Return the Gibbs energy in J/mol double gibbsmolar(void); - /// Return the Gibbs function in J/kg + /// Return the Gibbs energy in J/kg double gibbsmass(void){ return calc_gibbsmass(); }; + /// Return the excess Gibbs energy in J/mol + double gibbsmolar_excess(void); + /// Return the excess Gibbs energy in J/kg + double gibbsmass_excess(void) { return calc_gibbsmass_excess(); }; /// Return the Helmholtz energy in J/mol double helmholtzmolar(void); /// Return the Helmholtz energy in J/kg double helmholtzmass(void){ return calc_helmholtzmass(); }; + /// Return the excess Helmholtz energy in J/mol + double helmholtzmolar_excess(void); + /// Return the excess Helmholtz energy in J/kg + double helmholtzmass_excess(void) { return calc_helmholtzmass_excess(); }; + /// Return the excess volume in m^3/mol + double volumemolar_excess(void); + /// Return the excess volume in m^3/kg + double volumemass_excess(void) { return calc_volumemass_excess(); }; /// Return the speed of sound in m/s double speed_sound(void); /// Return the isothermal compressibility \f$ \kappa = -\frac{1}{v}\left.\frac{\partial v}{\partial p}\right|_T=\frac{1}{\rho}\left.\frac{\partial \rho}{\partial p}\right|_T\f$ in 1/Pa diff --git a/src/AbstractState.cpp b/src/AbstractState.cpp index 51f48462..0669f50c 100644 --- a/src/AbstractState.cpp +++ b/src/AbstractState.cpp @@ -152,6 +152,13 @@ bool AbstractState::clear() { this->_logp.clear(); this->_logrhomolar.clear(); + this->_hmolar_excess.clear(); + this->_smolar_excess.clear(); + this->_gibbsmolar_excess.clear(); + this->_volumemolar_excess.clear(); + this->_umolar_excess.clear(); + this->_helmholtzmolar_excess.clear(); + ///// Smoothing values //this->rhospline = -_HUGE; //this->dsplinedp = -_HUGE; @@ -463,22 +470,46 @@ double AbstractState::hmolar(void){ if (!_hmolar) _hmolar = calc_hmolar(); return _hmolar; } +double AbstractState::hmolar_excess(void) { + if (!_hmolar_excess) calc_excess_properties(); + return _hmolar_excess; +} double AbstractState::smolar(void){ if (!_smolar) _smolar = calc_smolar(); return _smolar; } +double AbstractState::smolar_excess(void) { + if (!_smolar_excess) calc_excess_properties(); + return _smolar_excess; +} double AbstractState::umolar(void){ if (!_umolar) _umolar = calc_umolar(); return _umolar; } +double AbstractState::umolar_excess(void) { + if (!_umolar_excess) calc_excess_properties(); + return _umolar_excess; +} double AbstractState::gibbsmolar(void){ if (!_gibbsmolar) _gibbsmolar = calc_gibbsmolar(); return _gibbsmolar; } +double AbstractState::gibbsmolar_excess(void) { + if (!_gibbsmolar_excess) calc_excess_properties(); + return _gibbsmolar_excess; +} double AbstractState::helmholtzmolar(void){ if (!_helmholtzmolar) _helmholtzmolar = calc_helmholtzmolar(); return _helmholtzmolar; } +double AbstractState::helmholtzmolar_excess(void) { + if (!_helmholtzmolar_excess) calc_excess_properties(); + return _helmholtzmolar_excess; +} +double AbstractState::volumemolar_excess(void) { + if (!_volumemolar_excess) calc_excess_properties(); + return _volumemolar_excess; +} double AbstractState::cpmolar(void){ if (!_cpmolar) _cpmolar = calc_cpmolar(); return _cpmolar; diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index b5ec3220..60f0f5e3 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -2518,6 +2518,28 @@ CoolPropDbl HelmholtzEOSMixtureBackend::calc_gibbsmolar(void) throw ValueError(format("phase is invalid in calc_gibbsmolar")); } } +void HelmholtzEOSMixtureBackend::calc_excess_properties(void) +{ + _gibbsmolar_excess = this->gibbsmolar(), + _smolar_excess = this->smolar(), + _hmolar_excess = this->hmolar(); + _umolar_excess = this->umolar(); + _volumemolar_excess = 1/this->rhomolar(); + for (std::size_t i = 0; i < components.size(); ++i) + { + transient_pure_state.reset(new HelmholtzEOSBackend(components[i].name)); + transient_pure_state->update(PT_INPUTS, p(), T()); + double x_i = mole_fractions[i]; + double R = gas_constant(); + _gibbsmolar_excess = static_cast(_gibbsmolar_excess) - x_i*(transient_pure_state->gibbsmolar() + R*T()*log(x_i)); + _hmolar_excess = static_cast(_hmolar_excess) - x_i*transient_pure_state->hmolar(); + _umolar_excess = static_cast(_umolar_excess) - x_i*transient_pure_state->umolar(); + _smolar_excess = static_cast(_smolar_excess) - x_i*(transient_pure_state->smolar() - R*log(x_i)); + _volumemolar_excess = static_cast(_volumemolar_excess) - x_i/transient_pure_state->rhomolar(); + } + _helmholtzmolar_excess = static_cast(_umolar_excess) - _T*static_cast(_smolar_excess); +} + CoolPropDbl HelmholtzEOSMixtureBackend::calc_helmholtzmolar(void) { if (isTwoPhase()) diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index f4424762..1dd9961f 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -24,6 +24,7 @@ protected: void pre_update(CoolProp::input_pairs &input_pair, CoolPropDbl &value1, CoolPropDbl &value2 ); void post_update(bool optional_checks = true); std::vector > linked_states; ///< States that are linked to this one, and should be updated (BIP, reference state, etc.) + shared_ptr transient_pure_state; ///< A temporary state used for calculations of pure fluid properties shared_ptr TPD_state; ///< A temporary state used for calculations of the tangent-plane-distance shared_ptr critical_state; ///< A temporary state used for calculations of the critical point(s) /// Update the state class used to calculate the tangent-plane-distance @@ -253,6 +254,7 @@ public: CoolPropDbl calc_cvmolar(void); CoolPropDbl calc_cpmolar(void); CoolPropDbl calc_gibbsmolar(void); + CoolPropDbl calc_helmholtzmolar(void); CoolPropDbl calc_cpmolar_idealgas(void); CoolPropDbl calc_pressure_nocache(CoolPropDbl T, CoolPropDbl rhomolar); @@ -260,11 +262,14 @@ public: CoolPropDbl calc_smolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar); CoolPropDbl calc_hmolar(void); - CoolPropDbl calc_hmolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar); + CoolPropDbl calc_hmolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar); CoolPropDbl calc_umolar_nocache(CoolPropDbl T, CoolPropDbl rhomolar); CoolPropDbl calc_umolar(void); CoolPropDbl calc_speed_sound(void); + + void calc_excess_properties(); + /** \brief The phase identification parameter of Venkatarathnam et al., FPE, 2011 * * diff --git a/src/Backends/REFPROP/REFPROPMixtureBackend.cpp b/src/Backends/REFPROP/REFPROPMixtureBackend.cpp index 81ec3a72..1d5005a4 100644 --- a/src/Backends/REFPROP/REFPROPMixtureBackend.cpp +++ b/src/Backends/REFPROP/REFPROPMixtureBackend.cpp @@ -1730,6 +1730,48 @@ CoolPropDbl REFPROPMixtureBackend::call_phi0dll(long itau, long idel) PHI0dll(&itau, &idel, &__T, &__rho, &(mole_fractions[0]), &val); return static_cast(val)/pow(tau,itau); // Not multplied by delta^idel } +/// Calculate excess properties +void REFPROPMixtureBackend::calc_excess_properties(){ + this->check_loaded_fluid(); + long ierr = 0; + char herr[255]; + double T_K = _T, p_kPa = _p/1000.0, rho=1, vE = -1, eE = -1, hE = -1, sE = -1, aE = -1, gE = -1; + long kph = 1; + + // subroutine EXCESS(t, p, x, kph, rho, vE, eE, hE, sE, aE, gE, ierr, herr) + // c + // c compute excess properties as a function of temperature, pressure, + // c and composition. + // c + // c inputs : + //c t--temperature[K] + // c p--pressure[kPa] + // c x--composition[array of mol frac] + // c kph--phase flag : 1 = liquid + // c 2 = vapor + // c 0 = stable phase + // c outputs : + //c rho--molar density[mol / L](if input less than 0, used as initial guess) + // c vE--excess volume[L / mol] + // c eE--excess energy[J / mol] + // c hE--excess enthalpy[J / mol] + // c sE--excess entropy[J / mol - K] + // c aE--excess Helmholtz energy[J / mol] + // c gE--excess Gibbs energy[J / mol] + // c ierr--error flag : 0 = successful + // c 55 = T, p inputs in different phase for the pure fluids + // c herr--error string(character * 255 variable if ierr<>0) + EXCESSdll(&T_K, &p_kPa, &(mole_fractions[0]), &kph, + &rho, &vE, &eE, &hE, &sE, &aE, &gE, + &ierr, herr, errormessagelength); // Error message + if (static_cast(ierr) > 0) { throw ValueError(format("EXCESSdll: %s", herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());} + _volumemolar_excess = vE; + _umolar_excess = eE; + _hmolar_excess = hE; + _smolar_excess = sE; + _helmholtzmolar_excess = aE; + _gibbsmolar_excess = gE; +} void REFPROP_SETREF(char hrf[3], long ixflag, double x0[1], double &h0, double &s0, double &T0, double &p0, long &ierr, char herr[255], long l1, long l2){ std::string err; diff --git a/src/Backends/REFPROP/REFPROPMixtureBackend.h b/src/Backends/REFPROP/REFPROPMixtureBackend.h index 48402885..dd449dc6 100644 --- a/src/Backends/REFPROP/REFPROPMixtureBackend.h +++ b/src/Backends/REFPROP/REFPROPMixtureBackend.h @@ -101,6 +101,8 @@ public: void check_loaded_fluid(void); + void calc_excess_properties(); + /// Returns true if REFPROP is supported on this platform static bool REFPROP_supported(void);