From 799422fe4028e662af85edba84c3f4c975e7125c Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Sun, 24 Aug 2014 17:26:46 +0200 Subject: [PATCH] Implemented generalized second derivative, seems to work, implemented at PropsSI level too Signed-off-by: Ian Bell --- include/AbstractState.h | 46 ++- include/DataStructures.h | 2 + .../Helmholtz/HelmholtzEOSMixtureBackend.cpp | 379 ++++++++---------- .../Helmholtz/HelmholtzEOSMixtureBackend.h | 18 +- src/CoolProp.cpp | 14 +- src/DataStructures.cpp | 45 ++- src/Tests/CoolProp-Tests.cpp | 35 +- 7 files changed, 296 insertions(+), 243 deletions(-) diff --git a/include/AbstractState.h b/include/AbstractState.h index feb51bba..522ba989 100644 --- a/include/AbstractState.h +++ b/include/AbstractState.h @@ -158,6 +158,14 @@ protected: virtual long double calc_d2alpha0_dDelta2(void){throw NotImplementedError("calc_d2alpha0_dDelta2 is not implemented for this backend");}; /// Using this backend, calculate the ideal-gas Helmholtz energy term \f$\alpha^0_{\tau\tau}\f$ (dimensionless) virtual long double calc_d2alpha0_dTau2(void){throw NotImplementedError("calc_d2alpha0_dTau2 is not implemented for this backend");}; + /// Using this backend, calculate the ideal-gas Helmholtz energy term \f$\alpha^0_{\delta\delta\delta}\f$ (dimensionless) + virtual long double calc_d3alpha0_dDelta3(void){throw NotImplementedError("calc_d3alpha0_dDelta3 is not implemented for this backend");}; + /// Using this backend, calculate the ideal-gas Helmholtz energy term \f$\alpha^0_{\delta\delta\tau}\f$ (dimensionless) + virtual long double calc_d3alpha0_dDelta2_dTau(void){throw NotImplementedError("calc_d3alpha0_dDelta2_dTau is not implemented for this backend");}; + /// Using this backend, calculate the ideal-gas Helmholtz energy term \f$\alpha^0_{\delta\tau\tau}\f$ (dimensionless) + virtual long double calc_d3alpha0_dDelta_dTau2(void){throw NotImplementedError("calc_d3alpha0_dDelta_dTau2 is not implemented for this backend");}; + /// Using this backend, calculate the ideal-gas Helmholtz energy term \f$\alpha^0_{\tau\tau\tau}\f$ (dimensionless) + virtual long double calc_d3alpha0_dTau3(void){throw NotImplementedError("calc_d3alpha0_dTau3 is not implemented for this backend");}; virtual void calc_reducing_state(void){throw NotImplementedError("calc_reducing_state is not implemented for this backend");}; @@ -308,22 +316,22 @@ public: * * The first partial derivative (\ref CoolProp::AbstractState::first_partial_deriv) can be expressed as * - * \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} = \frac{N}{D}\f] + * \f[ \left(\frac{\partial A}{\partial B}\right)_C = \frac{\left(\frac{\partial A}{\partial T}\right)_\rho\left(\frac{\partial C}{\partial \rho}\right)_T-\left(\frac{\partial A}{\partial \rho}\right)_T\left(\frac{\partial C}{\partial T}\right)_\rho}{\left(\frac{\partial B}{\partial T}\right)_\rho\left(\frac{\partial C}{\partial \rho}\right)_T-\left(\frac{\partial B}{\partial \rho}\right)_T\left(\frac{\partial C}{\partial T}\right)_\rho} = \frac{N}{D}\f] * * and the second derivative can be expressed as * * \f[ - * \frac{\partial}{\partial D}\left(\left(\frac{\partial A}{\partial B}\right)_C\right)_E = \frac{\frac{\partial}{\partial \tau}\left( \left(\frac{\partial A}{\partial B}\right)_C \right)_\delta\left(\frac{\partial E}{\partial \delta}\right)_\tau-\frac{\partial}{\partial \delta}\left(\left(\frac{\partial A}{\partial B}\right)_C\right)_\tau\left(\frac{\partial E}{\partial \tau}\right)_\delta}{\left(\frac{\partial D}{\partial \tau}\right)_\delta\left(\frac{\partial E}{\partial \delta}\right)_\tau-\left(\frac{\partial D}{\partial \delta}\right)_\tau\left(\frac{\partial E}{\partial \tau}\right)_\delta} + * \frac{\partial}{\partial D}\left(\left(\frac{\partial A}{\partial B}\right)_C\right)_E = \frac{\frac{\partial}{\partial T}\left( \left(\frac{\partial A}{\partial B}\right)_C \right)_\rho\left(\frac{\partial E}{\partial \rho}\right)_T-\frac{\partial}{\partial \rho}\left(\left(\frac{\partial A}{\partial B}\right)_C\right)_T\left(\frac{\partial E}{\partial T}\right)_\rho}{\left(\frac{\partial D}{\partial T}\right)_\rho\left(\frac{\partial E}{\partial \rho}\right)_T-\left(\frac{\partial D}{\partial \rho}\right)_T\left(\frac{\partial E}{\partial T}\right)_\rho} * \f] * * which can be expressed in parts as * - * \f[\left(\frac{\partial N}{\partial \delta}\right)_{\tau} = \left(\frac{\partial A}{\partial \tau}\right)_\delta\left(\frac{\partial^2 C}{\partial \delta^2}\right)_{\tau}+\left(\frac{\partial^2 A}{\partial \tau\partial\delta}\right)\left(\frac{\partial C}{\partial \delta}\right)_{\tau}-\left(\frac{\partial A}{\partial \delta}\right)_\tau\left(\frac{\partial^2 C}{\partial \tau\partial\delta}\right)-\left(\frac{\partial^2 A}{\partial \delta^2}\right)_{\tau}\left(\frac{\partial C}{\partial \tau}\right)_\delta\f] - * \f[\left(\frac{\partial D}{\partial \delta}\right)_{\tau} = \left(\frac{\partial B}{\partial \tau}\right)_\delta\left(\frac{\partial^2 C}{\partial \delta^2}\right)_{\tau}+\left(\frac{\partial^2 B}{\partial \tau\partial\delta}\right)\left(\frac{\partial C}{\partial \delta}\right)_{\tau}-\left(\frac{\partial B}{\partial \delta}\right)_\tau\left(\frac{\partial^2 C}{\partial \tau\partial\delta}\right)-\left(\frac{\partial^2 B}{\partial \delta^2}\right)_{\tau}\left(\frac{\partial C}{\partial \tau}\right)_\delta\f] - * \f[\left(\frac{\partial N}{\partial \tau}\right)_{\delta} = \left(\frac{\partial A}{\partial \tau}\right)_\delta\left(\frac{\partial^2 C}{\partial \delta\partial\tau}\right)+\left(\frac{\partial^2 A}{\partial \tau^2}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_{\tau}-\left(\frac{\partial A}{\partial \delta}\right)_\tau\left(\frac{\partial^2 C}{\partial \tau^2}\right)_\delta-\left(\frac{\partial^2 A}{\partial \delta\partial\tau}\right)\left(\frac{\partial C}{\partial \tau}\right)_\delta\f] - * \f[\left(\frac{\partial D}{\partial \tau}\right)_{\delta} = \left(\frac{\partial B}{\partial \tau}\right)_\delta\left(\frac{\partial^2 C}{\partial \delta\partial\tau}\right)+\left(\frac{\partial^2 B}{\partial \tau^2}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_{\tau}-\left(\frac{\partial B}{\partial \delta}\right)_\tau\left(\frac{\partial^2 C}{\partial \tau^2}\right)_\delta-\left(\frac{\partial^2 B}{\partial \delta\partial\tau}\right)\left(\frac{\partial C}{\partial \tau}\right)_\delta\f] - * \f[\frac{\partial}{\partial \tau}\left( \left(\frac{\partial A}{\partial B}\right)_C \right)_\delta = \frac{D\left(\frac{\partial N}{\partial \tau}\right)_{\delta}-N\left(\frac{\partial D}{\partial \tau}\right)_{\delta}}{D^2}\f] - * \f[\frac{\partial}{\partial \delta}\left( \left(\frac{\partial A}{\partial B}\right)_C \right)_\tau = \frac{D\left(\frac{\partial N}{\partial \delta}\right)_{\tau}-N\left(\frac{\partial D}{\partial \delta}\right)_{\tau}}{D^2}\f] + * \f[\left(\frac{\partial N}{\partial \rho}\right)_{T} = \left(\frac{\partial A}{\partial T}\right)_\rho\left(\frac{\partial^2 C}{\partial \rho^2}\right)_{T}+\left(\frac{\partial^2 A}{\partial T\partial\rho}\right)\left(\frac{\partial C}{\partial \rho}\right)_{T}-\left(\frac{\partial A}{\partial \rho}\right)_T\left(\frac{\partial^2 C}{\partial T\partial\rho}\right)-\left(\frac{\partial^2 A}{\partial \rho^2}\right)_{T}\left(\frac{\partial C}{\partial T}\right)_\rho\f] + * \f[\left(\frac{\partial D}{\partial \rho}\right)_{T} = \left(\frac{\partial B}{\partial T}\right)_\rho\left(\frac{\partial^2 C}{\partial \rho^2}\right)_{T}+\left(\frac{\partial^2 B}{\partial T\partial\rho}\right)\left(\frac{\partial C}{\partial \rho}\right)_{T}-\left(\frac{\partial B}{\partial \rho}\right)_T\left(\frac{\partial^2 C}{\partial T\partial\rho}\right)-\left(\frac{\partial^2 B}{\partial \rho^2}\right)_{T}\left(\frac{\partial C}{\partial T}\right)_\rho\f] + * \f[\left(\frac{\partial N}{\partial T}\right)_{\rho} = \left(\frac{\partial A}{\partial T}\right)_\rho\left(\frac{\partial^2 C}{\partial \rho\partial T}\right)+\left(\frac{\partial^2 A}{\partial T^2}\right)_\rho\left(\frac{\partial C}{\partial \rho}\right)_{T}-\left(\frac{\partial A}{\partial \rho}\right)_T\left(\frac{\partial^2 C}{\partial T^2}\right)_\rho-\left(\frac{\partial^2 A}{\partial \rho\partial T}\right)\left(\frac{\partial C}{\partial T}\right)_\rho\f] + * \f[\left(\frac{\partial D}{\partial T}\right)_{\rho} = \left(\frac{\partial B}{\partial T}\right)_\rho\left(\frac{\partial^2 C}{\partial \rho\partial T}\right)+\left(\frac{\partial^2 B}{\partial T^2}\right)_\rho\left(\frac{\partial C}{\partial \rho}\right)_{T}-\left(\frac{\partial B}{\partial \rho}\right)_T\left(\frac{\partial^2 C}{\partial T^2}\right)_\rho-\left(\frac{\partial^2 B}{\partial \rho\partial T}\right)\left(\frac{\partial C}{\partial T}\right)_\rho\f] + * \f[\frac{\partial}{\partial \rho}\left( \left(\frac{\partial A}{\partial B}\right)_C \right)_T = \frac{D\left(\frac{\partial N}{\partial \rho}\right)_{T}-N\left(\frac{\partial D}{\partial \rho}\right)_{\tau}}{D^2}\f] + * \f[\frac{\partial}{\partial T}\left( \left(\frac{\partial A}{\partial B}\right)_C \right)_\rho = \frac{D\left(\frac{\partial N}{\partial T}\right)_{\rho}-N\left(\frac{\partial D}{\partial T}\right)_{\rho}}{D^2}\f] * * The terms \f$ N \f$ and \f$ D \f$ are the numerator and denominator from \ref CoolProp::AbstractState::first_partial_deriv respectively */ @@ -476,12 +484,22 @@ public: if (!_d2alpha0_dTau2) _d2alpha0_dTau2 = calc_d2alpha0_dTau2(); return _d2alpha0_dTau2; }; - /* - virtual double d3alpha0_dDelta3(void) = 0; - virtual double d3alpha0_dDelta2_dTau(void) = 0; - virtual double d3alpha0_dDelta_dTau2(void) = 0; - virtual double d3alpha0_dTau3(void) = 0; - */ + long double d3alpha0_dTau3(void){ + if (!_d3alpha0_dTau3) _d3alpha0_dTau3 = calc_d3alpha0_dTau3(); + return _d3alpha0_dTau3; + }; + long double d3alpha0_dDelta_dTau2(void){ + if (!_d3alpha0_dDelta_dTau2) _d3alpha0_dDelta_dTau2 = calc_d3alpha0_dDelta_dTau2(); + return _d3alpha0_dDelta_dTau2; + }; + long double d3alpha0_dDelta2_dTau(void){ + if (!_d3alpha0_dDelta2_dTau) _d3alpha0_dDelta2_dTau = calc_d3alpha0_dDelta2_dTau(); + return _d3alpha0_dDelta2_dTau; + }; + long double d3alpha0_dDelta3(void){ + if (!_d3alpha0_dDelta3) _d3alpha0_dDelta3 = calc_d3alpha0_dDelta3(); + return _d3alpha0_dDelta3; + }; long double alphar(void){ if (!_alphar) _alphar = calc_alphar(); diff --git a/include/DataStructures.h b/include/DataStructures.h index 5325c645..bb590c16 100644 --- a/include/DataStructures.h +++ b/include/DataStructures.h @@ -89,6 +89,8 @@ bool is_valid_parameter(const std::string & name, parameters & iOutput); bool is_valid_first_derivative(const std::string & name, parameters &iOf, parameters &iWrt, parameters &iConstant); +bool is_valid_second_derivative(const std::string & name, parameters &iOf1, parameters &iWrt1, parameters &iConstant1, parameters &iWrt2, parameters &iConstant2); + std::string get_csv_parameter_list(); /// These are constants for the compositions diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 5ea8c2db..9607dafe 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -1243,120 +1243,11 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot throw ValueError(format("For now, we don't support T [%g K] below Ttriple [%g K]", _T, components[0]->pEOS->Ttriple)); } } -//void HelmholtzEOSMixtureBackend::DmolarT_phase_determination_pure_or_pseudopure() -//{ -// if (_T < _crit.T) -// { -// // Start to think about the saturation stuff -// // First try to use the ancillary equations if you are far enough away -// // You know how accurate the ancillary equations are thanks to using CoolProp code to refit them -// if (_rhomolar < 0.95*components[0]->ancillaries.rhoV.evaluate(_T)){ -// this->_phase = iphase_gas; return; -// } -// else if (_rhomolar > 1.05*components[0]->ancillaries.rhoL.evaluate(_T)){ -// this->_phase = iphase_liquid; return; -// } -// else{ -// // Actually have to use saturation information sadly -// // For the given temperature, find the saturation state -// // Run the saturation routines to determine the saturation densities and pressures -// HelmholtzEOSMixtureBackend HEOS(components); -// SaturationSolvers::saturation_T_pure_options options; -// SaturationSolvers::saturation_T_pure(&HEOS, _T, options); -// -// long double Q = (1/_rhomolar-1/HEOS.SatL->rhomolar())/(1/HEOS.SatV->rhomolar()-1/HEOS.SatL->rhomolar()); -// if (Q < -100*DBL_EPSILON){ -// this->_phase = iphase_liquid; -// } -// else if (Q > 1+100*DBL_EPSILON){ -// this->_phase = iphase_gas; -// } -// else{ -// this->_phase = iphase_twophase; -// } -// _Q = Q; -// // Load the outputs -// _p = _Q*HEOS.SatV->p() + (1-_Q)*HEOS.SatL->p(); -// _rhomolar = 1/(_Q/HEOS.SatV->rhomolar() + (1-_Q)/HEOS.SatL->rhomolar()); -// return; -// } -// } -// // Now check the states above the critical temperature. -// -// // Calculate the pressure if it is not already cached. -// calc_pressure(); -// -// if (_T > _crit.T && _p > _crit.p){ -// this->_phase = iphase_supercritical; return; -// } -// else if (_T > _crit.T && _p < _crit.p){ -// this->_phase = iphase_gas; return; -// } -// else if (_T < _crit.T && _p > _crit.p){ -// this->_phase = iphase_liquid; return; -// } -// /*else if (p < params.ptriple){ -// return iphase_gas; -// }*/ -// else{ -// throw ValueError(format("phase cannot be determined")); -// } -//} - -//void HelmholtzEOSMixtureBackend::PT_phase_determination() -//{ -// if (_T < _crit.T) -// { -// // Start to think about the saturation stuff -// // First try to use the ancillary equations if you are far enough away -// // Ancillary equations are good to within 1% in pressure in general -// // Some industrial fluids might not be within 3% -// if (_p > 1.05*components[0]->ancillaries.pL.evaluate(_T)){ -// this->_phase = iphase_liquid; return; -// } -// else if (_p < 0.95*components[0]->ancillaries.pV.evaluate(_T)){ -// this->_phase = iphase_gas; return; -// } -// else{ -// throw NotImplementedError("potentially two phase inputs not possible yet"); -// //// Actually have to use saturation information sadly -// //// For the given temperature, find the saturation state -// //// Run the saturation routines to determine the saturation densities and pressures -// //// Use the passed in variables to save calls to the saturation routine so the values can be re-used again -// //saturation_T(T, enabled_TTSE_LUT, pL, pV, rhoL, rhoV); -// //double Q = (1/rho-1/rhoL)/(1/rhoV-1/rhoL); -// //if (Q < -100*DBL_EPSILON){ -// // this->_phase = iphase_liquid; return; -// //} -// //else if (Q > 1+100*DBL_EPSILON){ -// // this->_phase = iphase_gas; return; -// //} -// //else{ -// // this->_phase = iphase_twophase; return; -// //} -// } -// } -// // Now check the states above the critical temperature. -// if (_T > _crit.T && _p > _crit.p){ -// this->_phase = iphase_supercritical; return; -// } -// else if (_T > _crit.T && _p < _crit.p){ -// this->_phase = iphase_gas; return; -// } -// else if (_T < _crit.T && _p > _crit.p){ -// this->_phase = iphase_liquid; return; -// } -// /*else if (p < params.ptriple){ -// return iphase_gas; -// }*/ -// else{ -// throw ValueError(format("phase cannot be determined")); -// } -//} - -void get_dtau_ddelta(HelmholtzEOSMixtureBackend *HEOS, long double T, long double rho, parameters index, long double &dtau, long double &ddelta) +void get_dT_drho(HelmholtzEOSMixtureBackend *HEOS, parameters index, long double &dT, long double &drho) { - long double rhor = HEOS->get_reducing_state().rhomolar, + long double T = HEOS->T(), + rho = HEOS->rhomolar(), + rhor = HEOS->get_reducing_state().rhomolar, Tr = HEOS->get_reducing_state().T, dT_dtau = -pow(T, 2)/Tr, R = HEOS->gas_constant(), @@ -1366,120 +1257,170 @@ void get_dtau_ddelta(HelmholtzEOSMixtureBackend *HEOS, long double T, long doubl switch (index) { case iT: - dtau = dT_dtau; ddelta = 0; break; + dT = 1; drho = 0; break; case iDmolar: - dtau = 0; ddelta = rhor; break; + dT = 0; drho = 1; break; case iP: - { - // dp/ddelta|tau - ddelta = rhor*R*T*(1+2*delta*HEOS->dalphar_dDelta()+pow(delta, 2)*HEOS->d2alphar_dDelta2()); - // dp/dtau|delta - dtau = dT_dtau*rho*R*(1+delta*HEOS->dalphar_dDelta() - tau*delta*HEOS->d2alphar_dDelta_dTau()); + { + // dp/drho|T + drho = R*T*(1+2*delta*HEOS->dalphar_dDelta()+pow(delta, 2)*HEOS->d2alphar_dDelta2()); + // dp/dT|rho + dT = rho*R*(1+delta*HEOS->dalphar_dDelta() - tau*delta*HEOS->d2alphar_dDelta_dTau()); break; - } + } case iHmolar: - // dh/dtau|delta - dtau = dT_dtau*R*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2()) + (1+delta*HEOS->dalphar_dDelta()-tau*delta*HEOS->d2alphar_dDelta_dTau())); - // dh/ddelta|tau - ddelta = rhor*T*R/rho*(tau*delta*HEOS->d2alphar_dDelta_dTau()+delta*HEOS->dalphar_dDelta()+pow(delta,2)*HEOS->d2alphar_dDelta2()); + // dh/dT|rho + dT = R*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2()) + (1+delta*HEOS->dalphar_dDelta()-tau*delta*HEOS->d2alphar_dDelta_dTau())); + // dh/drho|T + drho = T*R/rho*(tau*delta*HEOS->d2alphar_dDelta_dTau()+delta*HEOS->dalphar_dDelta()+pow(delta,2)*HEOS->d2alphar_dDelta2()); break; case iSmolar: - // ds/dtau|delta - dtau = dT_dtau*R/T*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2())); - // ds/ddelta|tau - ddelta = rhor*R/rho*(-(1+delta*HEOS->dalphar_dDelta()-tau*delta*HEOS->d2alphar_dDelta_dTau())); + // ds/dT|rho + dT = R/T*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2())); + // ds/drho|T + drho = R/rho*(-(1+delta*HEOS->dalphar_dDelta()-tau*delta*HEOS->d2alphar_dDelta_dTau())); break; case iUmolar: - // du/dtau|delta - dtau = dT_dtau*R*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2())); - // du/ddelta|tau - ddelta = rhor*HEOS->T()*R/rho*(tau*delta*HEOS->d2alphar_dDelta_dTau()); + // du/dT|rho + dT = R*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2())); + // du/drho|T + drho = HEOS->T()*R/rho*(tau*delta*HEOS->d2alphar_dDelta_dTau()); break; case iTau: - dtau = 1; ddelta = 0; break; + dT = 1/dT_dtau; drho = 0; break; case iDelta: - dtau = 0; ddelta = 1; break; + dT = 0; drho = 1/rhor; break; default: throw ValueError(format("input to get_dtau_ddelta[%s] is invalid",get_parameter_information(index,"short").c_str())); } } -void get_dtau_ddelta_second_derivatives(HelmholtzEOSMixtureBackend *HEOS, long double T, long double rho, int index, long double &dtau2, long double &ddelta_dtau, long double &ddelta2) +void get_dT_drho_second_derivatives(HelmholtzEOSMixtureBackend *HEOS, int index, long double &dT2, long double &drho_dT, long double &drho2) { -// long double rhor = HEOS->get_reducing_state().rhomolar, -// Tr = HEOS->get_reducing_state().T, -// dT_dtau = -pow(T, 2)/Tr, -// R = HEOS->gas_constant(), -// delta = rho/rhor, -// tau = Tr/T; -// -// HelmholtzDerivatives derivs; -// components[0]->pEOS->alphar.all(tau, delta, derivs); -// -// switch (index) -// { -// case iT: -// dtau2 = d2T_dtau2; // d2T_dtau2 -// ddelta_dtau = 0; // d2T_ddelta_dtau -// ddelta2 = 0; // d2T_ddelta2 -// break; -// case iDmolar: -// dtau2 = 0; // d2rhomolar_dtau2 -// ddelta2 = rhor; -// break; -// case iTau: -// dtau2 = 0; ddelta_dtau = 0; ddelta2 = 0; break; -// case iDelta: -// dtau2 = 0; ddelta_dtau = 0; ddelta2 = 0; break; -// case iP: -// { -// // dp/ddelta|tau -// d2pdrho2__T = R*T/rhomolar*(1 + 2*delta*dalphar_dDelta + pow(delta, 2)*d2alphar_dDelta2); -// ddelta2 = -// long double dalphar_dDelta = HEOS->calc_alphar_deriv_nocache(0, 1, HEOS->get_mole_fractions(), tau, delta); -// long double d2alphar_dDelta2 = HEOS->calc_alphar_deriv_nocache(0, 2, HEOS->get_mole_fractions(), tau, delta); -// long double d2alphar_dDelta_dTau = HEOS->calc_alphar_deriv_nocache(1, 1, HEOS->get_mole_fractions(), tau, delta); -// -// // dp/dtau|delta -// dtau = dT_dtau*rho*R*(1+delta*dalphar_dDelta-tau*delta*d2alphar_dDelta_dTau); -// break; -// } -// case iHmolar: -// // dh/dtau|delta -// dtau = dT_dtau*R*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2()) + (1+delta*HEOS->dalphar_dDelta()-tau*delta*HEOS->d2alphar_dDelta_dTau())); -// // dh/ddelta|tau -// ddelta = rhor*T*R/rho*(tau*delta*HEOS->d2alphar_dDelta_dTau()+delta*HEOS->dalphar_dDelta()+pow(delta,2)*HEOS->d2alphar_dDelta2()); -// break; -// case iSmolar: -// // ds/dtau|delta -// dtau = dT_dtau*R/T*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2())); -// // ds/ddelta|tau -// ddelta = rhor*R/rho*(-(1+delta*HEOS->dalphar_dDelta()-tau*delta*HEOS->d2alphar_dDelta_dTau())); -// break; -// case iUmolar: -// // du/dtau|delta -// dtau = dT_dtau*R*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2())); -// // du/ddelta|tau -// ddelta = rhor*HEOS->T()*R/rho*(tau*delta*HEOS->d2alphar_dDelta_dTau()); -// break; -// -// default: -// throw ValueError(format("input to get_dtau_ddelta[%s] is invalid",get_parameter_information(index,"short").c_str())); -// } -} -long double HelmholtzEOSMixtureBackend::calc_first_partial_deriv_nocache(long double T, long double rhomolar, parameters Of, parameters Wrt, parameters Constant) -{ - long double dOf_dtau, dOf_ddelta, dWrt_dtau, dWrt_ddelta, dConstant_dtau, dConstant_ddelta; + long double T = HEOS->T(), + rho = HEOS->rhomolar(), + rhor = HEOS->get_reducing_state().rhomolar, + Tr = HEOS->get_reducing_state().T, + dT_dtau = -pow(T, 2)/Tr, + R = HEOS->gas_constant(), + delta = rho/rhor, + tau = Tr/T; - get_dtau_ddelta(this, T, rhomolar, Of, dOf_dtau, dOf_ddelta); - get_dtau_ddelta(this, T, rhomolar, Wrt, dWrt_dtau, dWrt_ddelta); - get_dtau_ddelta(this, T, rhomolar, Constant, dConstant_dtau, dConstant_ddelta); - - //std::cout << dOf_dtau << " " << dOf_ddelta << " " << dWrt_dtau << " " << dWrt_ddelta << " " << dConstant_dtau << " " << dConstant_ddelta << std::endl; - return (dOf_dtau*dConstant_ddelta-dOf_ddelta*dConstant_dtau)/(dWrt_dtau*dConstant_ddelta-dWrt_ddelta*dConstant_dtau); + // Here we use T and rho as independent variables since derivations are already done by Thorade, 2013, + // Partial derivatives of thermodynamic state propertiesfor dynamic simulation, DOI 10.1007/s12665-013-2394-z + + HelmholtzDerivatives derivs = HEOS->get_components()[0]->pEOS->alphar.all(tau, delta); + + switch (index) + { + case iT: + dT2 = 0; // d2T_dT2 + drho_dT = 0; // d2T_drho_dT + drho2 = 0; + break; + case iDmolar: + dT2 = 0; // d2rhomolar_dtau2 + drho2 = 0; + drho_dT = 0; + break; + case iTau: + dT2 = 0; drho_dT = 0; drho2 = 0; break; + case iDelta: + dT2 = 0; drho_dT = 0; drho2 = 0; break; + case iP: + { + drho2 = T*R/rho*(2*delta*derivs.dalphar_ddelta+4*pow(delta,2)*derivs.d2alphar_ddelta2+pow(delta,3)*derivs.d3alphar_ddelta3); + dT2 = rho*R/T*(pow(tau,2)*delta*derivs.d3alphar_ddelta_dtau2); + drho_dT = R*(1+2*delta*derivs.dalphar_ddelta +pow(delta,2)*derivs.d2alphar_ddelta2 - 2*delta*tau*derivs.d2alphar_ddelta_dtau - tau*pow(delta, 2)*derivs.d3alphar_ddelta2_dtau); + break; + } + case iHmolar: + { + // d2h/drho2|T + drho2 = R*T/pow(rho, 2)*pow(delta,2)*(tau*HEOS->d3alpha0_dDelta2_dTau() + 2*derivs.d2alphar_ddelta2 + delta*derivs.d2alphar_ddelta2); + // d2h/dT2|rho + dT2 = R/T*pow(tau, 2)*(tau*(HEOS->d3alpha0_dTau3()+derivs.d3alphar_dtau3) + 2*(HEOS->d2alpha0_dTau2()+derivs.d2alphar_dtau2) + delta*derivs.d3alphar_ddelta_dtau2); + // d2h/drho/dT + drho_dT = R/rho*delta*(delta*derivs.d2alphar_ddelta2 - pow(tau,2)*derivs.d3alphar_ddelta_dtau2 + derivs.dalphar_ddelta - tau*delta*derivs.d3alphar_ddelta2_dtau - tau*derivs.d2alphar_ddelta_dtau); + break; + } + case iSmolar: + { + // d2s/rho2|T + drho2 = R/pow(rho,2)*(1-pow(delta,2)*derivs.d2alphar_ddelta2 + tau*pow(delta,2)*derivs.d3alphar_ddelta2_dtau); + // d2s/dT2|rho + dT2 = R*pow(tau/T, 2)*(tau*(HEOS->d3alpha0_dTau3()+derivs.d3alphar_dtau3)+3*(HEOS->d2alpha0_dTau2()+derivs.d2alphar_dtau2)); + // d2s/drho/dT + drho_dT = R/(T*rho)*(-pow(tau,2)*delta*derivs.d3alphar_ddelta_dtau2); + break; + } + case iUmolar: + { + // d2u/rho2|T + drho2 = R*T*tau*pow(delta/rho,2)*derivs.d3alphar_ddelta2_dtau; + // d2u/dT2|rho + dT2 = R/T*pow(tau, 2)*(tau*(HEOS->d3alpha0_dTau3()+derivs.d3alphar_dtau3)+2*(HEOS->d2alpha0_dTau2()+derivs.d2alphar_dtau2)); + // d2u/drho/dT + drho_dT = R/rho*(-pow(tau,2)*delta*derivs.d3alphar_ddelta_dtau2); + break; + } + default: + 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) { - return calc_first_partial_deriv_nocache(_T, _rhomolar, Of, Wrt, 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) +{ + long double dOf1_dT, dOf1_drho, dWrt1_dT, dWrt1_drho, dConstant1_dT, dConstant1_drho, d2Of1_dT2, d2Of1_drhodT, + d2Of1_drho2, d2Wrt1_dT2, d2Wrt1_drhodT, d2Wrt1_drho2, d2Constant1_dT2, d2Constant1_drhodT, d2Constant1_drho2, + dWrt2_dT, dWrt2_drho, dConstant2_dT, dConstant2_drho, N, D, dNdrho__T, dDdrho__T, dNdT__rho, dDdT__rho, + dderiv1_drho, dderiv1_dT, second; + + // First and second partials needed for terms involved in first derivative + get_dT_drho(this, Of1, dOf1_dT, dOf1_drho); + get_dT_drho(this, Wrt1, dWrt1_dT, dWrt1_drho); + get_dT_drho(this, Constant1, dConstant1_dT, dConstant1_drho); + get_dT_drho_second_derivatives(this, Of1, d2Of1_dT2, d2Of1_drhodT, d2Of1_drho2); + get_dT_drho_second_derivatives(this, Wrt1, d2Wrt1_dT2, d2Wrt1_drhodT, d2Wrt1_drho2); + get_dT_drho_second_derivatives(this, Constant1, d2Constant1_dT2, d2Constant1_drhodT, d2Constant1_drho2); + + // First derivatives of terms involved in the second derivative + get_dT_drho(this, Wrt2, dWrt2_dT, dWrt2_drho); + get_dT_drho(this, Constant2, dConstant2_dT, dConstant2_drho); + + // Numerator and denominator of first partial derivative term + N = dOf1_dT*dConstant1_drho - dOf1_drho*dConstant1_dT; + D = dWrt1_dT*dConstant1_drho - dWrt1_drho*dConstant1_dT; + + // Derivatives of the numerator and denominator of the first partial derivative term with respect to rho, T held constant + // They are of similar form, with Of1 and Wrt1 swapped + dNdrho__T = dOf1_dT*d2Constant1_drho2 + d2Of1_drhodT*dConstant1_drho - dOf1_drho*d2Constant1_drhodT - d2Of1_drho2*dConstant1_dT; + dDdrho__T = dWrt1_dT*d2Constant1_drho2 + d2Wrt1_drhodT*dConstant1_drho - dWrt1_drho*d2Constant1_drhodT - d2Wrt1_drho2*dConstant1_dT; + + // Derivatives of the numerator and denominator of the first partial derivative term with respect to T, rho held constant + // They are of similar form, with Of1 and Wrt1 swapped + dNdT__rho = dOf1_dT*d2Constant1_drhodT + d2Of1_dT2*dConstant1_drho - dOf1_drho*d2Constant1_dT2 - d2Of1_drhodT*dConstant1_dT; + dDdT__rho = dWrt1_dT*d2Constant1_drhodT + d2Wrt1_dT2*dConstant1_drho - dWrt1_drho*d2Constant1_dT2 - d2Wrt1_drhodT*dConstant1_dT; + + // First partial of first derivative term with respect to T + dderiv1_drho = (D*dNdrho__T - N*dDdrho__T)/pow(D, 2); + + // First partial of first derivative term with respect to rho + dderiv1_dT = (D*dNdT__rho - N*dDdT__rho)/pow(D, 2); + + // Complete second derivative + second = (dderiv1_dT*dConstant2_drho - dderiv1_drho*dConstant2_dT)/(dWrt2_dT*dConstant2_drho - dWrt2_drho*dConstant2_dT); + + return second; } long double HelmholtzEOSMixtureBackend::calc_pressure_nocache(long double T, long double rhomolar) @@ -1495,10 +1436,6 @@ long double HelmholtzEOSMixtureBackend::calc_pressure_nocache(long double T, lon // Get pressure return rhomolar*gas_constant()*T*(1+delta*dalphar_dDelta); } -//long double HelmholtzEOSMixtureBackend::solver_for_T_given_rho_oneof_PHSU(long double T, long double value, int other, int rhomin, int rhomax) -//{ -// -//} long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long double T, long double value, int other) { long double ymelt, yc, ymin, y; @@ -2332,6 +2269,26 @@ long double HelmholtzEOSMixtureBackend::calc_d2alpha0_dTau2(void) const int nTau = 2, nDelta = 0; return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar); } +long double HelmholtzEOSMixtureBackend::calc_d3alpha0_dDelta3(void) +{ + const int nTau = 0, nDelta = 3; + return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar); +} +long double HelmholtzEOSMixtureBackend::calc_d3alpha0_dDelta2_dTau(void) +{ + const int nTau = 1, nDelta = 2; + return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar); +} +long double HelmholtzEOSMixtureBackend::calc_d3alpha0_dDelta_dTau2(void) +{ + const int nTau = 2, nDelta = 1; + return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar); +} +long double HelmholtzEOSMixtureBackend::calc_d3alpha0_dTau3(void) +{ + const int nTau = 3, nDelta = 0; + return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar); +} long double HelmholtzEOSMixtureBackend::mixderiv_dalphar_dxi(int i) diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index 720a1a3c..94772302 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -134,7 +134,11 @@ public: 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_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); long double calc_alpha0(void); long double calc_dalpha0_dDelta(void); @@ -142,6 +146,10 @@ public: long double calc_d2alpha0_dDelta2(void); long double calc_d2alpha0_dDelta_dTau(void); long double calc_d2alpha0_dTau2(void); + long double calc_d3alpha0_dDelta3(void); + long double calc_d3alpha0_dDelta2_dTau(void); + long double calc_d3alpha0_dDelta_dTau2(void); + long double calc_d3alpha0_dTau3(void); //long double calc_surface_tension(void); long double calc_viscosity(void); @@ -206,12 +214,8 @@ public: \f] */ long double calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant); - - /** - This version doesn't use any cached values - \sa calc_first_partial_deriv - */ - long double calc_first_partial_deriv_nocache(long double T, long double rhomolar, 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/CoolProp.cpp b/src/CoolProp.cpp index e35d9b4c..ee3b2a0d 100644 --- a/src/CoolProp.cpp +++ b/src/CoolProp.cpp @@ -405,10 +405,16 @@ double _PropsSI(const std::string &Output, const std::string &Name1, double Prop // Return the value return val; } - else{// if is_valid_second_derivative(Output, iOf1, iWrt1, iConstant1, iWrt2, iConstant2){ - return _HUGE; - }; - + else if (is_valid_second_derivative(Output, iOf1, iWrt1, iConstant1, iWrt2, iConstant2)){ + // Return the desired output + double val = State->second_partial_deriv(iOf1, iWrt1, iConstant1, iWrt2, iConstant2); + + // Return the value + return val; + } + else{ + throw ValueError(format("Output [%s] is not a parameter or a string representation of a derivative",Output.c_str()).c_str()); + } } double PropsSI(const std::string &Output, const std::string &Name1, double Prop1, const std::string &Name2, double Prop2, const std::string &Ref, const std::vector &z) { diff --git a/src/DataStructures.cpp b/src/DataStructures.cpp index 3f92165a..b84ad823 100644 --- a/src/DataStructures.cpp +++ b/src/DataStructures.cpp @@ -77,13 +77,13 @@ parameter_info parameter_info_list[] = { parameter_info(iZ, "Z","O","-","Compressibility factor",false), parameter_info(ifundamental_derivative_of_gas_dynamics, "fundamental_derivative_of_gas_dynamics","O","-","Fundamental_derivative_of_gas_dynamics",false), - parameter_info(ialphar, "","O","-","Residual Helmholtz energy",false), - parameter_info(idalphar_dtau_constdelta, "","O","-","Derivative of residual Helmholtz energy with tau",false), - parameter_info(idalphar_ddelta_consttau, "","O","-","Derivative of residual Helmholtz energy with delta",false), + parameter_info(ialphar, "alphar","O","-","Residual Helmholtz energy",false), + parameter_info(idalphar_dtau_constdelta, "dalphar_dtau_constdelta","O","-","Derivative of residual Helmholtz energy with tau",false), + parameter_info(idalphar_ddelta_consttau, "dalphar_ddelta_consttau","O","-","Derivative of residual Helmholtz energy with delta",false), - parameter_info(ialpha0, "","O","-","Ideal Helmholtz energy",false), - parameter_info(idalpha0_dtau_constdelta, "","O","-","Derivative of ideal Helmholtz energy with tau",false), - parameter_info(idalpha0_ddelta_consttau, "","O","-","Derivative of ideal Helmholtz energy with delta",false), + parameter_info(ialpha0, "alpha0","O","-","Ideal Helmholtz energy",false), + parameter_info(idalpha0_dtau_constdelta, "dalpha0_dtau_constdelta","O","-","Derivative of ideal Helmholtz energy with tau",false), + parameter_info(idalpha0_ddelta_consttau, "dalpha0_ddelta_consttau","O","-","Derivative of ideal Helmholtz energy with delta",false), }; @@ -250,6 +250,39 @@ bool is_valid_first_derivative(const std::string & name, parameters &iOf, parame return false; } } + +bool is_valid_second_derivative(const std::string & name, parameters &iOf1, parameters &iWrt1, parameters &iConstant1, parameters &iWrt2, parameters &iConstant2) +{ + if (get_debug_level() > 5){std::cout << format("is_valid_second_derivative(%s)",name.c_str());} + + // Suppose we start with "d(d(P)/d(Dmolar)|T)/d(Dmolar)|T" + std::size_t i_bar = name.rfind('|'); + if (i_bar == std::string::npos){return false;} + std::string constant2 = name.substr(i_bar+1); // "T" + if (!is_valid_parameter(constant2, iConstant2)){return false;}; + std::string left_of_bar = name.substr(0, i_bar); // "d(d(P)/d(Dmolar)|T)/d(Dmolar)" + + std::size_t i_slash = left_of_bar.rfind('/'); + if (i_slash == std::string::npos){return false;} + + std::string left_of_slash = left_of_bar.substr(0, i_slash); // "d(d(P)/d(Dmolar)|T)" + std::size_t iN0 = left_of_slash.find("("); + std::size_t iN1 = left_of_slash.rfind(")"); + if (!(iN0 > 0 && iN1 > 0 && iN1 > iN0)){return false;} + std::string num = left_of_slash.substr(iN0+1, iN1-2); // "d(P)/d(Dmolar)|T" + if (!is_valid_first_derivative(num, iOf1, iWrt1, iConstant1)){return false;} + + std::string right_of_slash = left_of_bar.substr(i_slash+1); // "d(Dmolar)" + std::size_t iD0 = right_of_slash.find("("); + std::size_t iD1 = right_of_slash.rfind(")"); + if (!(iD0 > 0 && iD1 > 0 && iD1 > iD0)){return false;} + std::string den = right_of_slash.substr(iD0+1, iD1-2); // "Dmolar" + if (!is_valid_parameter(den, iWrt2)){return false;} + + // If we haven't quit yet, all is well + return true; +} + int get_parameter_index(const std::string ¶m_name) { parameters iOutput; diff --git a/src/Tests/CoolProp-Tests.cpp b/src/Tests/CoolProp-Tests.cpp index a2f46007..78a01c09 100644 --- a/src/Tests/CoolProp-Tests.cpp +++ b/src/Tests/CoolProp-Tests.cpp @@ -916,7 +916,7 @@ TEST_CASE("Multiple calls to state class are consistent", "[flashdups],[flash],[ } } -TEST_CASE("Test partial derivatives using PropsSI", "[derivatives]") +TEST_CASE("Test first partial derivatives using PropsSI", "[derivatives]") { double hmolar, hmass, T = 300; SECTION("Check drhodp|T 3 ways","") @@ -985,6 +985,39 @@ TEST_CASE("Test partial derivatives using PropsSI", "[derivatives]") } } +TEST_CASE("Test second partial derivatives using PropsSI", "[derivatives]") +{ + double hmolar, hmass, T = 300; + SECTION("Check d2pdrho2|T 3 ways","") + { + shared_ptr AS(CoolProp::AbstractState::factory("HEOS", "Water")); + double rhomolar = 60000; + AS->update(CoolProp::DmolarT_INPUTS, rhomolar, T); + double p = AS->p(); + + double d2pdrhomolar2__T_AbstractState = AS->second_partial_deriv(CoolProp::iP, CoolProp::iDmolar, CoolProp::iT, CoolProp::iDmolar, CoolProp::iT); + // Centered second derivative + double del = 1e0; + double d2pdrhomolar2__T_PropsSI_num = (PropsSI("P","T",T,"Dmolar",rhomolar+del,"Water") - 2*PropsSI("P","T",T,"Dmolar",rhomolar,"Water") + PropsSI("P","T",T,"Dmolar",rhomolar-del,"Water"))/pow(del, 2); + double d2pdrhomolar2__T_PropsSI = PropsSI("d(d(P)/d(Dmolar)|T)/d(Dmolar)|T","T",T,"Dmolar",rhomolar,"Water"); + + CAPTURE(d2pdrhomolar2__T_AbstractState); + CAPTURE(d2pdrhomolar2__T_PropsSI_num); + double rel_err_exact = std::abs((d2pdrhomolar2__T_AbstractState-d2pdrhomolar2__T_PropsSI)/d2pdrhomolar2__T_PropsSI); + double rel_err_approx = std::abs((d2pdrhomolar2__T_PropsSI_num-d2pdrhomolar2__T_AbstractState)/d2pdrhomolar2__T_AbstractState); + CHECK(rel_err_exact < 1e-5); + CHECK(rel_err_approx < 1e-5); + } + SECTION("Valid second partial derivatives","") + { + CHECK(ValidNumber(PropsSI("d(d(Hmolar)/d(P)|T)/d(T)|Dmolar","T",300,"P",101325,"n-Propane"))); + } + SECTION("Invalid second partial derivatives","") + { + CHECK(!ValidNumber(PropsSI("d(d()/d(P)|T)/d()|","T",300,"P",101325,"n-Propane"))); + CHECK(!ValidNumber(PropsSI("dd(Dmolar)/d()|T)|T","T",300,"P",101325,"n-Propane"))); + } +} TEST_CASE("Ancillary functions", "[ancillary]") {