From 1a9f0a4e321cc7024d4ccbcab4c7374d3149e669 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Tue, 9 Dec 2014 18:55:52 -0500 Subject: [PATCH] Derivatives are working, seems like a bug in REFPROP PHI0dll function is the only remaining problem Signed-off-by: Ian Bell --- include/AbstractState.h | 53 ++++++++-- src/AbstractState.cpp | 96 +++++++++++++++++++ .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 30 ++++-- .../Helmholtz/HelmholtzEOSMixtureBackend.h | 41 ++++---- .../REFPROP/REFPROPMixtureBackend.cpp | 65 ++++++++++--- src/Backends/REFPROP/REFPROPMixtureBackend.h | 43 ++++++++- src/Backends/REFPROP/REFPROP_lib.h | 28 ++++++ 7 files changed, 305 insertions(+), 51 deletions(-) diff --git a/include/AbstractState.h b/include/AbstractState.h index a635af6c..71404122 100644 --- a/include/AbstractState.h +++ b/include/AbstractState.h @@ -147,7 +147,15 @@ protected: virtual long double calc_d2alphar_dDelta_dTau(void){throw NotImplementedError("calc_d2alphar_dDelta_dTau is not implemented for this backend");}; /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\tau\tau}\f$ (dimensionless) virtual long double calc_d2alphar_dTau2(void){throw NotImplementedError("calc_d2alphar_dTau2 is not implemented for this backend");}; - + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta\delta\delta}\f$ (dimensionless) + virtual long double calc_d3alphar_dDelta3(void){throw NotImplementedError("calc_d3alpha0_dDelta3 is not implemented for this backend");}; + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta\delta\tau}\f$ (dimensionless) + virtual long double calc_d3alphar_dDelta2_dTau(void){throw NotImplementedError("calc_d3alpha0_dDelta2_dTau is not implemented for this backend");}; + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta\tau\tau}\f$ (dimensionless) + virtual long double calc_d3alphar_dDelta_dTau2(void){throw NotImplementedError("calc_d3alpha0_dDelta_dTau2 is not implemented for this backend");}; + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\tau\tau\tau}\f$ (dimensionless) + virtual long double calc_d3alphar_dTau3(void){throw NotImplementedError("calc_d3alpha0_dTau3 is not implemented for this backend");}; + // Derivatives of ideal-gas helmholtz energy /// Using this backend, calculate the ideal-gas Helmholtz energy term \f$\alpha^0\f$ (dimensionless) virtual long double calc_alpha0(void){throw NotImplementedError("calc_alpha0 is not implemented for this backend");}; @@ -195,7 +203,7 @@ protected: virtual long double calc_physical_hazard(void){throw NotImplementedError("calc_physical_hazard is not implemented for this backend");}; /// Calculate the first partial derivative for the desired derivative - virtual long double calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant){throw NotImplementedError("calc_first_partial_deriv is not implemented for this backend");}; + virtual long double calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant); /// Calculate the second partial derivative using the given backend virtual long double calc_second_partial_deriv(parameters Of1, parameters Wrt1, parameters Constant1, parameters Of2, parameters Constant2){throw NotImplementedError("calc_second_partial_deriv is not implemented for this backend");}; @@ -225,10 +233,14 @@ protected: /// Using this backend, get the critical point temperature in K virtual long double calc_T_critical(void){throw NotImplementedError("calc_T_critical is not implemented for this backend");}; + /// Using this backend, get the reducing point temperature in K + virtual long double calc_T_reducing(void){throw NotImplementedError("calc_T_reducing is not implemented for this backend");}; /// Using this backend, get the critical point pressure in Pa virtual long double calc_p_critical(void){throw NotImplementedError("calc_p_critical is not implemented for this backend");}; /// Using this backend, get the critical point molar density in mol/m^3 virtual long double calc_rhomolar_critical(void){throw NotImplementedError("calc_rhomolar_critical is not implemented for this backend");}; + /// Using this backend, get the reducing point molar density in mol/m^3 + virtual long double calc_rhomolar_reducing(void){throw NotImplementedError("calc_rhomolar_reducing is not implemented for this backend");}; /// Using this backend, construct the phase envelope, the variable type describes the type of phase envelope to be built. virtual void calc_phase_envelope(const std::string &type){throw NotImplementedError("calc_phase_envelope is not implemented for this backend");}; @@ -341,6 +353,8 @@ protected: For mixtures, it is the exact critical point temperature calculated by the methods of Michelsen( \todo fill in reference) */ double T_critical(void); + /// Return the reducing point temperature in K + double T_reducing(void); /// Return the critical pressure in Pa /** For pure fluids, this is the critical point pressure as defined by the author of the EOS @@ -360,6 +374,18 @@ protected: * For mixtures, it is the exact critical point molar density calculated by the methods of Michelsen( \todo fill in reference) */ double rhomass_critical(void); + /** \brief Return the molar density at the reducing point in mol/m^3 + * + * For pure fluids, this is the critical point molar density + * For mixtures, it is the exact critical point molar density calculated by the methods of Michelsen( \todo fill in reference) + */ + double rhomolar_reducing(void); + /** \brief Return the mass density at the reducing point in kg/m^3 + * + * For pure fluids, this is the critical point molar density + * For mixtures, it is the exact critical point molar density calculated by the methods of Michelsen( \todo fill in reference) + */ + double rhomass_reducing(void); /// Return the triple point pressure double p_triple(void); @@ -617,12 +643,25 @@ protected: if (!_d2alphar_dTau2) _d2alphar_dTau2 = calc_d2alphar_dTau2(); return _d2alphar_dTau2; }; + long double d3alphar_dDelta3(void){ + if (!_d3alphar_dDelta3) _d3alphar_dDelta3 = calc_d3alphar_dDelta3(); + return _d3alphar_dDelta3; + }; + long double d3alphar_dDelta2_dTau(void){ + if (!_d3alphar_dDelta2_dTau) _d3alphar_dDelta2_dTau = calc_d3alphar_dDelta2_dTau(); + return _d3alphar_dDelta2_dTau; + }; + long double d3alphar_dDelta_dTau2(void){ + if (!_d3alphar_dDelta_dTau2) _d3alphar_dDelta_dTau2 = d3alphar_dDelta_dTau2(); + return _d3alphar_dDelta_dTau2; + }; + long double d3alphar_dTau3(void){ + if (!_d3alphar_dTau3) _d3alphar_dTau3 = calc_d3alphar_dTau3(); + return _d3alphar_dTau3; + }; + + /* - virtual double d3alphar_dDelta3(void) = 0; - virtual double d3alphar_dDelta2_dTau(void) = 0; - virtual double d3alphar_dDelta_dTau2(void) = 0; - virtual double d3alphar_dTau3(void) = 0; - virtual double dalphar_dDelta_lim(void) = 0; virtual double d2alphar_dDelta2_lim(void) = 0; virtual double d2alphar_dDelta_dTau_lim(void) = 0; diff --git a/src/AbstractState.cpp b/src/AbstractState.cpp index 40876ddb..3892573d 100644 --- a/src/AbstractState.cpp +++ b/src/AbstractState.cpp @@ -313,6 +313,9 @@ double AbstractState::pmax(void){ double AbstractState::T_critical(void){ return calc_T_critical(); } +double AbstractState::T_reducing(void){ + return calc_T_reducing(); +} double AbstractState::p_critical(void){ return calc_p_critical(); } @@ -325,6 +328,12 @@ double AbstractState::rhomolar_critical(void){ double AbstractState::rhomass_critical(void){ return calc_rhomolar_critical()*molar_mass(); } +double AbstractState::rhomolar_reducing(void){ + return calc_rhomolar_reducing(); +} +double AbstractState::rhomass_reducing(void){ + return calc_rhomolar_reducing()*molar_mass(); +} double AbstractState::hmolar(void){ if (!_hmolar) _hmolar = calc_hmolar(); return _hmolar; @@ -398,7 +407,94 @@ double AbstractState::dBvirial_dT(void){ return calc_dBvirial_dT(); } double AbstractState::dCvirial_dT(void){ return calc_dCvirial_dT(); } double AbstractState::compressibility_factor(void){ return calc_compressibility_factor(); } +// Get the derivatives of the parameters in the partial derivative with respect to T and rho +void get_dT_drho(AbstractState &AS, parameters index, long double &dT, long double &drho) +{ + long double T = AS.T(), + rho = AS.rhomolar(), + rhor = AS.rhomolar_reducing(), + Tr = AS.T_reducing(), + dT_dtau = -pow(T, 2)/Tr, + R = AS.gas_constant(), + delta = rho/rhor, + tau = Tr/T; + + switch (index) + { + case iT: + dT = 1; drho = 0; break; + case iDmolar: + dT = 0; drho = 1; break; + case iDmass: + dT = 0; drho = AS.molar_mass(); break; + case iP: + { + // dp/drho|T + drho = R*T*(1+2*delta*AS.dalphar_dDelta()+pow(delta, 2)*AS.d2alphar_dDelta2()); + // dp/dT|rho + dT = rho*R*(1+delta*AS.dalphar_dDelta() - tau*delta*AS.d2alphar_dDelta_dTau()); + break; + } + case iHmass: + case iHmolar: + { + // dh/dT|rho + dT = R*(-pow(tau,2)*(AS.d2alpha0_dTau2()+AS.d2alphar_dTau2()) + (1+delta*AS.dalphar_dDelta()-tau*delta*AS.d2alphar_dDelta_dTau())); + // dh/drhomolar|T + drho = T*R/rho*(tau*delta*AS.d2alphar_dDelta_dTau()+delta*AS.dalphar_dDelta()+pow(delta,2)*AS.d2alphar_dDelta2()); + if (index == iHmass){ + // dhmolar/drhomolar|T * dhmass/dhmolar where dhmass/dhmolar = 1/mole mass + drho /= AS.molar_mass(); + dT /= AS.molar_mass(); + } + break; + } + case iSmass: + case iSmolar: + { + // ds/dT|rho + dT = R/T*(-pow(tau,2)*(AS.d2alpha0_dTau2()+AS.d2alphar_dTau2())); + // ds/drho|T + drho = R/rho*(-(1+delta*AS.dalphar_dDelta()-tau*delta*AS.d2alphar_dDelta_dTau())); + if (index == iSmass){ + // ds/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass + drho /= AS.molar_mass(); + dT /= AS.molar_mass(); + } + break; + } + case iUmass: + case iUmolar: + { + // du/dT|rho + dT = R*(-pow(tau,2)*(AS.d2alpha0_dTau2()+AS.d2alphar_dTau2())); + // du/drho|T + drho = AS.T()*R/rho*(tau*delta*AS.d2alphar_dDelta_dTau()); + if (index == iUmass){ + // du/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass + drho /= AS.molar_mass(); + dT /= AS.molar_mass(); + } + break; + } + case iTau: + dT = 1/dT_dtau; drho = 0; break; + case iDelta: + dT = 0; drho = 1/rhor; break; + default: + throw ValueError(format("input to get_dT_drho[%s] is invalid",get_parameter_information(index,"short").c_str())); + } +} +long double AbstractState::calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant) +{ + long double dOf_dT, dOf_drho, dWrt_dT, dWrt_drho, dConstant_dT, dConstant_drho; + get_dT_drho(*this, Of, dOf_dT, dOf_drho); + get_dT_drho(*this, Wrt, dWrt_dT, dWrt_drho); + get_dT_drho(*this, Constant, dConstant_dT, dConstant_drho); + + return (dOf_dT*dConstant_drho-dOf_drho*dConstant_dT)/(dWrt_dT*dConstant_drho-dWrt_drho*dConstant_dT); +} // // ---------------------------------------- // // Smoothing functions for density // // ---------------------------------------- diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 3c829c30..0f6688e0 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -1685,16 +1685,6 @@ void get_dT_drho_second_derivatives(HelmholtzEOSMixtureBackend *HEOS, int index, throw ValueError(format("input to get_dT_drho_second_derivatives[%s] is invalid", get_parameter_information(index,"short").c_str())); } } -long double HelmholtzEOSMixtureBackend::calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant) -{ - long double dOf_dT, dOf_drho, dWrt_dT, dWrt_drho, dConstant_dT, dConstant_drho; - - get_dT_drho(this, Of, dOf_dT, dOf_drho); - get_dT_drho(this, Wrt, dWrt_dT, dWrt_drho); - get_dT_drho(this, Constant, dConstant_dT, dConstant_drho); - - return (dOf_dT*dConstant_drho-dOf_drho*dConstant_dT)/(dWrt_dT*dConstant_drho-dWrt_drho*dConstant_dT); -} long double HelmholtzEOSMixtureBackend::calc_second_partial_deriv(parameters Of1, parameters Wrt1, parameters Constant1, parameters Wrt2, parameters Constant2) { @@ -2611,6 +2601,26 @@ long double HelmholtzEOSMixtureBackend::calc_d2alphar_dDelta2(void) calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta); return static_cast(_d2alphar_dDelta2); } +long double HelmholtzEOSMixtureBackend::calc_d3alphar_dDelta3(void) +{ + calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta); + return static_cast(_d3alphar_dDelta3); +} +long double HelmholtzEOSMixtureBackend::calc_d3alphar_dDelta2_dTau(void) +{ + calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta); + return static_cast(_d3alphar_dDelta2_dTau); +} +long double HelmholtzEOSMixtureBackend::calc_d3alphar_dDelta_dTau2(void) +{ + calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta); + return static_cast(_d3alphar_dDelta_dTau2); +} +long double HelmholtzEOSMixtureBackend::calc_d3alphar_dTau3(void) +{ + calc_all_alphar_deriv_cache(mole_fractions, _tau, _delta); + return static_cast(_d3alphar_dTau3); +} long double HelmholtzEOSMixtureBackend::calc_alpha0(void) { const int nTau = 0, nDelta = 0; diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index d0374d92..68e98aa4 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -176,16 +176,26 @@ public: long double calc_fugacity_coefficient(int i); long double calc_phase_identification_parameter(void); + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r\f$ (dimensionless) long double calc_alphar(void); - long double calc_dalphar_dDelta(void); - long double calc_dalphar_dTau(void); - long double calc_d2alphar_dDelta2(void); - long double calc_d2alphar_dDelta_dTau(void); - long double calc_d2alphar_dTau2(void); - long double calc_d3alphar_dDelta3(void); - long double calc_d3alphar_dDelta2_dTau(void); - long double calc_d3alphar_dDelta_dTau2(void); - long double calc_d3alphar_dTau3(void); + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta}\f$ (dimensionless) + long double calc_dalphar_dDelta(void); + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\tau}\f$ (dimensionless) + long double calc_dalphar_dTau(void); + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta\delta}\f$ (dimensionless) + long double calc_d2alphar_dDelta2(void); + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta\tau}\f$ (dimensionless) + long double calc_d2alphar_dDelta_dTau(void); + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\tau\tau}\f$ (dimensionless) + long double calc_d2alphar_dTau2(void); + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta\delta\delta}\f$ (dimensionless) + long double calc_d3alphar_dDelta3(void); + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta\delta\tau}\f$ (dimensionless) + long double calc_d3alphar_dDelta2_dTau(void); + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta\tau\tau}\f$ (dimensionless) + long double calc_d3alphar_dDelta_dTau2(void); + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\tau\tau\tau}\f$ (dimensionless) + long double calc_d3alphar_dTau3(void); long double calc_alpha0(void); long double calc_dalpha0_dDelta(void); @@ -216,10 +226,13 @@ public: long double calc_Tmax_sat(void); void calc_Tmin_sat(long double &Tmin_satL, long double &Tmin_satV); void calc_pmin_sat(long double &pmin_satL, long double &pmin_satV); - + long double calc_T_critical(void); long double calc_p_critical(void); long double calc_rhomolar_critical(void); + + long double calc_T_reducing(void){return get_reducing_state().T;}; + long double calc_rhomolar_reducing(void){return get_reducing_state().rhomolar;}; std::string calc_name(void); @@ -257,14 +270,6 @@ public: const CoolProp::SimpleState & get_reducing_state(){calc_reducing_state(); return _reducing;}; - /** - In a general way we can calculate any first partial derivative based on calculating derivatives with respect to the fundamental variables of the EOS, \f$\tau\f$ and \f$\delta\f$ - \f[ - \left(\frac{\partial A}{\partial B}\right)_C = \frac{\left(\frac{\partial A}{\partial \tau}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_\tau-\left(\frac{\partial A}{\partial \delta}\right)_\tau\left(\frac{\partial C}{\partial \tau}\right)_\delta}{\left(\frac{\partial B}{\partial \tau}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_\tau-\left(\frac{\partial B}{\partial \delta}\right)_\tau\left(\frac{\partial C}{\partial \tau}\right)_\delta} - \f] - */ - long double calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant); - long double calc_second_partial_deriv(parameters Of1, parameters Wrt1, parameters Constant1, parameters Wrt2, parameters Constant2); void update_states(); diff --git a/src/Backends/REFPROP/REFPROPMixtureBackend.cpp b/src/Backends/REFPROP/REFPROPMixtureBackend.cpp index ab3615fc..da605015 100644 --- a/src/Backends/REFPROP/REFPROPMixtureBackend.cpp +++ b/src/Backends/REFPROP/REFPROPMixtureBackend.cpp @@ -169,11 +169,14 @@ static char default_reference_state[] = "DEF"; MELTTdll_POINTER MELTTdll; MLTH2Odll_POINTER MLTH2Odll; NAMEdll_POINTER NAMEdll; + PASSCMNdll_POINTER PASSCMNdll; PDFL1dll_POINTER PDFL1dll; PDFLSHdll_POINTER PDFLSHdll; PEFLSHdll_POINTER PEFLSHdll; PHFL1dll_POINTER PHFL1dll; PHFLSHdll_POINTER PHFLSHdll; + PHIXdll_POINTER PHIXdll; + PHI0dll_POINTER PHI0dll; PQFLSHdll_POINTER PQFLSHdll; PREOSdll_POINTER PREOSdll; PRESSdll_POINTER PRESSdll; @@ -183,6 +186,7 @@ static char default_reference_state[] = "DEF"; QMASSdll_POINTER QMASSdll; QMOLEdll_POINTER QMOLEdll; RESIDUALdll_POINTER RESIDUALdll; + RMIX2dll_POINTER RMIX2dll; SATDdll_POINTER SATDdll; SATEdll_POINTER SATEdll; SATHdll_POINTER SATHdll; @@ -296,11 +300,14 @@ double setFunctionPointers() MELTTdll = (MELTTdll_POINTER) getFunctionPointer(MELTTdll_NAME); MLTH2Odll = (MLTH2Odll_POINTER) getFunctionPointer(MLTH2Odll_NAME); NAMEdll = (NAMEdll_POINTER) getFunctionPointer(NAMEdll_NAME); - PDFL1dll = (PDFL1dll_POINTER) getFunctionPointer(PDFL1dll_NAME); + PASSCMNdll = (PASSCMNdll_POINTER) getFunctionPointer(PASSCMNdll_NAME); + PDFL1dll = (PDFL1dll_POINTER) getFunctionPointer(PDFL1dll_NAME); PDFLSHdll = (PDFLSHdll_POINTER) getFunctionPointer(PDFLSHdll_NAME); PEFLSHdll = (PEFLSHdll_POINTER) getFunctionPointer(PEFLSHdll_NAME); PHFL1dll = (PHFL1dll_POINTER) getFunctionPointer(PHFL1dll_NAME); PHFLSHdll = (PHFLSHdll_POINTER) getFunctionPointer(PHFLSHdll_NAME); + PHIXdll = (PHIXdll_POINTER) getFunctionPointer(PHIXdll_NAME); + PHI0dll = (PHI0dll_POINTER) getFunctionPointer(PHI0dll_NAME); PQFLSHdll = (PQFLSHdll_POINTER) getFunctionPointer(PQFLSHdll_NAME); PREOSdll = (PREOSdll_POINTER) getFunctionPointer(PREOSdll_NAME); PRESSdll = (PRESSdll_POINTER) getFunctionPointer(PRESSdll_NAME); @@ -308,6 +315,7 @@ double setFunctionPointers() PSFLSHdll = (PSFLSHdll_POINTER) getFunctionPointer(PSFLSHdll_NAME); PUREFLDdll = (PUREFLDdll_POINTER) getFunctionPointer(PUREFLDdll_NAME); RESIDUALdll = (RESIDUALdll_POINTER) getFunctionPointer(RESIDUALdll_NAME); + RMIX2dll = (RMIX2dll_POINTER) getFunctionPointer(RMIX2dll_NAME); QMASSdll = (QMASSdll_POINTER) getFunctionPointer(QMASSdll_NAME); QMOLEdll = (QMOLEdll_POINTER) getFunctionPointer(QMOLEdll_NAME); SATDdll = (SATDdll_POINTER) getFunctionPointer(SATDdll_NAME); @@ -708,6 +716,26 @@ long double REFPROPMixtureBackend::calc_rhomolar_critical(){ CRITPdll(&(mole_fractions[0]),&Tcrit,&pcrit_kPa,&dcrit_mol_L,&ierr,herr,255); if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());} return static_cast(dcrit_mol_L*1000); }; +long double REFPROPMixtureBackend::calc_T_reducing(){ + if (mole_fractions.size() != 1){throw ValueError("calc_T_reducing only valid for one component");}; + long ierr = 0, i, i1 = 0, i2 = 1, i3 = 0; + char herr[255], h[255], input[] = "tz"; + double Tz; + std::vector z(255); + PASSCMNdll(input, &i1, &i2, &i3, h,&i,&Tz,&(z[0]),&ierr,herr,255,255,255); + if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());} + return static_cast(Tz); +}; +long double REFPROPMixtureBackend::calc_rhomolar_reducing(){ + if (mole_fractions.size() != 1){throw ValueError("calc_rhomolar_reducing only valid for one component");}; + long ierr = 0, i = 0, i1 = 0, i2 = 1, i3 = 0; + char herr[255] = "", h[255] = "", input[255] = "rhoz"; + double rhored_mol_L; + std::vector z(20); + PASSCMNdll(input, &i1, &i2, &i3, h, &i, &rhored_mol_L,&(z[0]),&ierr,herr,255,255,255); + if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());} + return static_cast(rhored_mol_L*1000); +}; long double REFPROPMixtureBackend::calc_Ttriple(){ if (mole_fractions.size() != 1){throw ValueError("calc_Ttriple cannot be evaluated for mixtures");} long icomp = 0; @@ -715,6 +743,11 @@ long double REFPROPMixtureBackend::calc_Ttriple(){ INFOdll(&icomp, &wmm, &ttrp, &tnbpt, &tc, &pc, &Dc, &Zc, &acf, &dip, &Rgas); return static_cast(ttrp); }; +long double REFPROPMixtureBackend::calc_gas_constant(){ + double Rmix = 0; + RMIX2dll(&(mole_fractions[0]), &Rmix); + return static_cast(Rmix); +}; long double REFPROPMixtureBackend::calc_molar_mass(void) { double wmm_kg_kmol; @@ -833,19 +866,6 @@ long double REFPROPMixtureBackend::calc_cpmolar_idealgas(void) THERM0dll(&_T,&rho_mol_L,&(mole_fractions[0]),&p0,&e0,&h0,&s0,&cv0,&cp0,&w0,&A0,&G0); return static_cast(cp0); } -long double REFPROPMixtureBackend::calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant) -{ - if (Of == iP && Wrt == iT && (Constant == iDmolar || Constant == iDmass)) - { - double rho_mol_L = 0.001*_rhomolar; - double dpt; - DPDTdll(&_T, &rho_mol_L, &(mole_fractions[0]), &dpt); - return static_cast(dpt*1000); - } - else{ - throw ValueError(format("These derivative terms are not supported")); - } -} void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double value1, double value2) { @@ -1378,6 +1398,23 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu _cvmolar = cvmol; _cpmolar = cpmol; _speed_sound = w; + _tau = calc_T_critical()/_T; + _delta = _rhomolar/calc_rhomolar_critical(); +} +long double REFPROPMixtureBackend::call_phixdll(long itau, long idel) +{ + double val = 0, tau = _tau, delta = _delta; + if (PHIXdll == NULL){throw ValueError("PHIXdll function is not available in your version of REFPROP. Please upgrade");} + PHIXdll(&itau, &idel, &tau, &delta, &(mole_fractions[0]), &val); + return static_cast(val)/pow(static_cast(_delta),idel)/pow(static_cast(_tau),itau); +} +long double REFPROPMixtureBackend::call_phi0dll(long itau, long idel) +{ + double val = 0, tau = _tau, delta = _delta, __T = T(), __rho = rhomolar()/1000; + double dT_dtau; + if (PHI0dll == NULL){throw ValueError("PHI0dll function is not available in your version of REFPROP. Please upgrade");} + PHI0dll(&itau, &idel, &__T, &__rho, &(mole_fractions[0]), &val); + return static_cast(val)/pow(delta,idel)/pow(tau,itau); } } /* namespace CoolProp */ diff --git a/src/Backends/REFPROP/REFPROPMixtureBackend.h b/src/Backends/REFPROP/REFPROPMixtureBackend.h index 602451c6..4b6e3e7e 100644 --- a/src/Backends/REFPROP/REFPROPMixtureBackend.h +++ b/src/Backends/REFPROP/REFPROPMixtureBackend.h @@ -22,6 +22,12 @@ protected: static bool _REFPROP_supported; std::vector mole_fractions, mass_fractions; std::vector mole_fractions_liq, mole_fractions_vap; + + /// Call the PHIXdll function in the dll + long double call_phixdll(long itau, long idelta); + /// Call the PHI0dll function in the dll + long double call_phi0dll(long itau, long idelta); + public: REFPROPMixtureBackend(){}; @@ -58,8 +64,6 @@ public: long double calc_cpmolar_idealgas(void); - long double calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant); - /// Set the fluids in REFPROP DLL by calling the SETUPdll function /** @param fluid_names The vector of strings of the fluid components, without file ending @@ -98,9 +102,12 @@ public: bool has_melting_line(){return true;}; double calc_melt_Tmax(); long double calc_T_critical(void); + long double calc_T_reducing(void); long double calc_p_critical(void); long double calc_rhomolar_critical(void); + long double calc_rhomolar_reducing(void); long double calc_Ttriple(void); + long double calc_gas_constant(void); /// A wrapper function to calculate the limits for the EOS void limits(double &Tmin, double &Tmax, double &rhomolarmax, double &pmax); @@ -108,6 +115,38 @@ public: long double calc_pmax(void); /// Calculate the maximum temperature long double calc_Tmax(void); + + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r\f$ (dimensionless) + long double calc_alphar(void){return call_phixdll(0,0);}; + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta}\f$ (dimensionless) + long double calc_dalphar_dDelta(void){ return call_phixdll(0,1); }; + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\tau}\f$ (dimensionless) + long double calc_dalphar_dTau(void){ return call_phixdll(1,0); }; + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta\delta}\f$ (dimensionless) + long double calc_d2alphar_dDelta2(void){ return call_phixdll(0,2); }; + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta\tau}\f$ (dimensionless) + long double calc_d2alphar_dDelta_dTau(void){ return call_phixdll(1,1); }; + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\tau\tau}\f$ (dimensionless) + long double calc_d2alphar_dTau2(void){ return call_phixdll(2,0); }; + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta\delta\delta}\f$ (dimensionless) + long double calc_d3alphar_dDelta3(void){ return call_phixdll(0,3); }; + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta\delta\tau}\f$ (dimensionless) + long double calc_d3alphar_dDelta2_dTau(void){ return call_phixdll(1,2); }; + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\delta\tau\tau}\f$ (dimensionless) + long double calc_d3alphar_dDelta_dTau2(void){ return call_phixdll(2,1); }; + /// Using this backend, calculate the residual Helmholtz energy term \f$\alpha^r_{\tau\tau\tau}\f$ (dimensionless) + long double calc_d3alphar_dTau3(void){ return call_phixdll(3,0); }; + + long double calc_alpha0(void){ return call_phi0dll(0,0); }; + long double calc_dalpha0_dDelta(void){ return call_phi0dll(0,1); }; + long double calc_dalpha0_dTau(void){ return call_phi0dll(1,0); }; + long double calc_d2alpha0_dDelta2(void){ return call_phi0dll(0,2); }; + long double calc_d2alpha0_dDelta_dTau(void){ return call_phi0dll(1,1); }; + long double calc_d2alpha0_dTau2(void){ return call_phi0dll(2,0); }; + long double calc_d3alpha0_dDelta3(void){ return call_phi0dll(0,3); }; + long double calc_d3alpha0_dDelta2_dTau(void){ return call_phi0dll(1,2); }; + long double calc_d3alpha0_dDelta_dTau2(void){ return call_phi0dll(2,1); }; + long double calc_d3alpha0_dTau3(void){ return call_phi0dll(3,0); }; }; } /* namespace CoolProp */ diff --git a/src/Backends/REFPROP/REFPROP_lib.h b/src/Backends/REFPROP/REFPROP_lib.h index 372f896a..b420ba0d 100644 --- a/src/Backends/REFPROP/REFPROP_lib.h +++ b/src/Backends/REFPROP/REFPROP_lib.h @@ -75,11 +75,14 @@ # define MELTTdll MELTTdll # define MLTH2Odll MLTH2Odll # define NAMEdll NAMEdll +# define PASSCMNdll PASSCMN # define PDFL1dll PDFL1dll # define PDFLSHdll PDFLSHdll # define PEFLSHdll PEFLSHdll # define PHFL1dll PHFL1dll # define PHFLSHdll PHFLSHdll +# define PHIXdll PHIXdll +# define PHI0dll PHI0dll # define PQFLSHdll PQFLSHdll # define PREOSdll PREOSdll # define PRESSdll PRESSdll @@ -89,6 +92,7 @@ # define QMASSdll QMASSdll # define QMOLEdll QMOLEdll # define RESIDUALdll RESIDUALdll +# define RMIX2dll RMIX2dll # define SATDdll SATDdll # define SATEdll SATEdll # define SATHdll SATHdll @@ -191,11 +195,14 @@ # define MELTTdll melttdll_ # define MLTH2Odll mlth2odll_ # define NAMEdll namedll_ +# define PASSCMNdll passcmn_ # define PDFL1dll pdfl1dll_ # define PDFLSHdll pdflshdll_ # define PEFLSHdll peflshdll_ # define PHFL1dll phfl1dll_ # define PHFLSHdll phflshdll_ +# define PHIXdll phixdll_ +# define PHI0dll phi0dll_ # define PQFLSHdll pqflshdll_ # define PREOSdll preosdll_ # define PRESSdll pressdll_ @@ -205,6 +212,7 @@ # define QMASSdll qmassdll_ # define QMOLEdll qmoledll_ # define RESIDUALdll residualdll_ +# define RMIX2dll rmix2dll_ # define SATDdll satddll_ # define SATEdll satedll_ # define SATHdll sathdll_ @@ -298,11 +306,14 @@ # define MELTTdll melttdll_ # define MLTH2Odll mlth2odll_ # define NAMEdll namedll_ +# define PASSCMNdll passcmn_ # define PDFL1dll pdfl1dll_ # define PDFLSHdll pdflshdll_ # define PEFLSHdll peflshdll_ # define PHFL1dll phfl1dll_ # define PHFLSHdll phflshdll_ +# define PHIXdll phixdll_ +# define PHI0dll phi0dll_ # define PQFLSHdll pqflshdll_ # define PREOSdll preosdll_ # define PRESSdll pressdll_ @@ -312,6 +323,7 @@ # define QMASSdll qmassdll_ # define QMOLEdll qmoledll_ # define RESIDUALdll residualdll_ +# define RMIX2dll rmix2dll_ # define SATDdll satddll_ # define SATEdll satedll_ # define SATHdll sathdll_ @@ -403,11 +415,14 @@ # define MELTTdll melttdll # define MLTH2Odll mlth2odll # define NAMEdll namedll +# define PASSCMNdll passcmn # define PDFL1dll pdfl1dll # define PDFLSHdll pdflshdll # define PEFLSHdll peflshdll # define PHFL1dll phfl1dll # define PHFLSHdll phflshdll +# define PHIXdll phixdll +# define PHI0dll phi0dll # define PQFLSHdll pqflshdll # define PREOSdll preosdll # define PRESSdll pressdll @@ -417,6 +432,7 @@ # define QMASSdll qmassdll # define QMOLEdll qmoledll # define RESIDUALdll residualdll +# define RMIX2dll rmix2dll # define SATDdll satddll # define SATEdll satedll # define SATHdll sathdll @@ -520,11 +536,14 @@ #define MELTTdll_NAME FUNCTION_NAME(MELTTdll) #define MLTH2Odll_NAME FUNCTION_NAME(MLTH2Odll) #define NAMEdll_NAME FUNCTION_NAME(NAMEdll) +#define PASSCMNdll_NAME FUNCTION_NAME(PASSCMNdll) #define PDFL1dll_NAME FUNCTION_NAME(PDFL1dll) #define PDFLSHdll_NAME FUNCTION_NAME(PDFLSHdll) #define PEFLSHdll_NAME FUNCTION_NAME(PEFLSHdll) #define PHFL1dll_NAME FUNCTION_NAME(PHFL1dll) #define PHFLSHdll_NAME FUNCTION_NAME(PHFLSHdll) +#define PHIXdll_NAME FUNCTION_NAME(PHIXdll) +#define PHI0dll_NAME FUNCTION_NAME(PHI0dll) #define PQFLSHdll_NAME FUNCTION_NAME(PQFLSHdll) #define PREOSdll_NAME FUNCTION_NAME(PREOSdll) #define PRESSdll_NAME FUNCTION_NAME(PRESSdll) @@ -534,6 +553,7 @@ #define QMASSdll_NAME FUNCTION_NAME(QMASSdll) #define QMOLEdll_NAME FUNCTION_NAME(QMOLEdll) #define RESIDUALdll_NAME FUNCTION_NAME(RESIDUALdll) +#define RMIX2dll_NAME FUNCTION_NAME(RMIX2dll) #define SATDdll_NAME FUNCTION_NAME(SATDdll) #define SATEdll_NAME FUNCTION_NAME(SATEdll) #define SATHdll_NAME FUNCTION_NAME(SATHdll) @@ -634,11 +654,14 @@ extern "C" { typedef void (CALLCONV MELTTdll_TYPE)(double *,double *,double *,long *,char*,long ); typedef void (CALLCONV MLTH2Odll_TYPE)(double *,double *,double *); typedef void (CALLCONV NAMEdll_TYPE)(long *,char*,char*,char*,long ,long ,long ); + typedef void (CALLCONV PASSCMNdll_TYPE)(char *,long *,long *,long *,char *,long*,double *, double *, long*, char*, long, long, long); typedef void (CALLCONV PDFL1dll_TYPE)(double *,double *,double *,double *,long *,char*,long ); typedef void (CALLCONV PDFLSHdll_TYPE)(double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,long *,char*,long ); typedef void (CALLCONV PEFLSHdll_TYPE)(double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,long *,char*,long ); typedef void (CALLCONV PHFL1dll_TYPE)(double *,double *,double *,long *,double *,double *,long *,char*,long ); typedef void (CALLCONV PHFLSHdll_TYPE)(double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,long *,char*,long ); + typedef void (CALLCONV PHIXdll_TYPE)(long *,long *,double *,double *,double *, double *); + typedef void (CALLCONV PHI0dll_TYPE)(long *,long *,double *,double *,double *, double *); typedef void (CALLCONV PQFLSHdll_TYPE)(double *,double *,double *,long *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,long *,char*,long ); typedef void (CALLCONV PREOSdll_TYPE)(long *); typedef void (CALLCONV PRESSdll_TYPE)(double *,double *,double *,double *); @@ -648,6 +671,7 @@ extern "C" { typedef void (CALLCONV QMASSdll_TYPE)(double *,double *,double *,double *,double *,double *,double *,double *,long *,char*,long ); typedef void (CALLCONV QMOLEdll_TYPE)(double *,double *,double *,double *,double *,double *,double *,double *,long *,char*,long ); typedef void (CALLCONV RESIDUALdll_TYPE)(double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *); + typedef void (CALLCONV RMIX2dll_TYPE)(double *,double *); typedef void (CALLCONV SATDdll_TYPE)(double *,double *,long *,long *,double *,double *,double *,double *,double *,double *,long *,char*,long ); typedef void (CALLCONV SATEdll_TYPE)(double *,double *,long *,long *,long *,double *,double *,double *,long *,double *,double *,double *,long *,char*,long ); typedef void (CALLCONV SATHdll_TYPE)(double *,double *,long *,long *,long *,double *,double *,double *,long *,double *,double *,double *,long *,char*,long ); @@ -843,11 +867,14 @@ extern "C" { typedef MELTTdll_TYPE * MELTTdll_POINTER; typedef MLTH2Odll_TYPE * MLTH2Odll_POINTER; typedef NAMEdll_TYPE * NAMEdll_POINTER; + typedef PASSCMNdll_TYPE * PASSCMNdll_POINTER; typedef PDFL1dll_TYPE * PDFL1dll_POINTER; typedef PDFLSHdll_TYPE * PDFLSHdll_POINTER; typedef PEFLSHdll_TYPE * PEFLSHdll_POINTER; typedef PHFL1dll_TYPE * PHFL1dll_POINTER; typedef PHFLSHdll_TYPE * PHFLSHdll_POINTER; + typedef PHIXdll_TYPE * PHIXdll_POINTER; + typedef PHI0dll_TYPE * PHI0dll_POINTER; typedef PQFLSHdll_TYPE * PQFLSHdll_POINTER; typedef PREOSdll_TYPE * PREOSdll_POINTER; typedef PRESSdll_TYPE * PRESSdll_POINTER; @@ -857,6 +884,7 @@ extern "C" { typedef QMASSdll_TYPE * QMASSdll_POINTER; typedef QMOLEdll_TYPE * QMOLEdll_POINTER; typedef RESIDUALdll_TYPE * RESIDUALdll_POINTER; + typedef RMIX2dll_TYPE * RMIX2dll_POINTER; typedef SATDdll_TYPE * SATDdll_POINTER; typedef SATEdll_TYPE * SATEdll_POINTER; typedef SATHdll_TYPE * SATHdll_POINTER;