diff --git a/include/AbstractState.h b/include/AbstractState.h index 1964e222..58970cec 100644 --- a/include/AbstractState.h +++ b/include/AbstractState.h @@ -327,6 +327,9 @@ protected: /// Using this backend, get the molar density in mol/m^3 virtual long double calc_rhomolar(void){return _rhomolar;} + /// Using this backend, return true critical point where dp/drho|T = 0 and d2p/drho^2|T = 0 + virtual void calc_true_critical_point(double &T, double &rho){throw NotImplementedError("calc_true_critical_point is not implemented for this backend");}; + public: @@ -540,6 +543,9 @@ public: //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" double PIP(){return calc_PIP();}; + + /// Find the "true" critical point where dpdrho|T and d2p/drho2|T are equal to zero + void true_critical_point(double &T, double &rho){ calc_true_critical_point(T, rho); } std::vector mole_fractions_liquid(void){return calc_mole_fractions_liquid();}; std::vector mole_fractions_vapor(void){return calc_mole_fractions_vapor();}; diff --git a/include/Solvers.h b/include/Solvers.h index ab8e0b7b..770fe735 100644 --- a/include/Solvers.h +++ b/include/Solvers.h @@ -61,7 +61,7 @@ inline double Halley(FuncWrapper1DWithTwoDerivs &f, double x0, double ftol, int } // Multi-Dimensional solvers -std::vector NDNewtonRaphson_Jacobian(FuncWrapperND *f, const std::vector &x0, double tol, int maxiter, std::string *errstring); +std::vector NDNewtonRaphson_Jacobian(FuncWrapperND *f, std::vector &x0, double tol, int maxiter, std::string *errstring); }; /*namespace CoolProp*/ #endif diff --git a/src/Backends/REFPROP/REFPROPMixtureBackend.cpp b/src/Backends/REFPROP/REFPROPMixtureBackend.cpp index dba67c86..306d1e73 100644 --- a/src/Backends/REFPROP/REFPROPMixtureBackend.cpp +++ b/src/Backends/REFPROP/REFPROPMixtureBackend.cpp @@ -30,6 +30,7 @@ surface tension N/m #include "Exceptions.h" #include "Configuration.h" #include "CoolProp.h" +#include "Solvers.h" #include #include @@ -1265,6 +1266,34 @@ void REFPROP_SETREF(char hrf[3], long ixflag, double x0[1], double &h0, double & SETREFdll(hrf, &ixflag, x0, &h0, &s0, &T0, &p0, &ierr, herr, l1, l2); } +void REFPROPMixtureBackend::calc_true_critical_point(double &T, double &rho) +{ + + class wrapper : public FuncWrapperND + { + public: + const std::vector z; + wrapper(const std::vector &z) : z(z) {}; + std::vector call(const std::vector& x){ + std::vector r(2); + double dpdrho__constT = _HUGE, d2pdrho2__constT = _HUGE; + DPDDdll(const_cast(&(x[0])), const_cast(&(x[1])), const_cast(&(z[0])), &dpdrho__constT); + DPDD2dll(const_cast(&(x[0])), const_cast(&(x[1])), const_cast(&(z[0])), &d2pdrho2__constT); + r[0] = dpdrho__constT; + r[1] = d2pdrho2__constT; + return r; + }; + }; + wrapper resid(mole_fractions); + + T = calc_T_critical(); + double rho_moldm3 = calc_rhomolar_critical()/1000.0; + std::vector x(2,T); x[1] = rho_moldm3; + std::string errstr; + std::vector xfinal = NDNewtonRaphson_Jacobian(&resid, x, 1e-9, 30, &errstr); + T = xfinal[0]; rho = xfinal[1]*1000.0; +} + } /* namespace CoolProp */ diff --git a/src/Backends/REFPROP/REFPROPMixtureBackend.h b/src/Backends/REFPROP/REFPROPMixtureBackend.h index f45204e0..a144ac06 100644 --- a/src/Backends/REFPROP/REFPROPMixtureBackend.h +++ b/src/Backends/REFPROP/REFPROPMixtureBackend.h @@ -130,6 +130,9 @@ public: CoolPropDbl calc_Ttriple(void); CoolPropDbl calc_gas_constant(void); + /// Calculate the "true" critical point where dp/drho|T and d2p/drho2|T are zero + void calc_true_critical_point(double &T, double &rho); + /// A wrapper function to calculate the limits for the EOS void limits(double &Tmin, double &Tmax, double &rhomolarmax, double &pmax); /// Calculate the maximum pressure diff --git a/src/Solvers.cpp b/src/Solvers.cpp index 0558f03e..81aaf433 100644 --- a/src/Solvers.cpp +++ b/src/Solvers.cpp @@ -4,6 +4,7 @@ #include "MatrixMath.h" #include #include "CoolPropTools.h" +#include namespace CoolProp{ @@ -48,24 +49,32 @@ functions, each of which take the vector x. The data is managed using std::vecto @param errstring A string with the returned error. If the length of errstring is zero, no errors were found @returns If no errors are found, the solution. Otherwise, _HUGE, the value for infinity */ -std::vector NDNewtonRaphson_Jacobian(FuncWrapperND *f, std::vector x0, double tol, int maxiter, std::string *errstring) +std::vector NDNewtonRaphson_Jacobian(FuncWrapperND *f, std::vector &x0, double tol, int maxiter, std::string *errstring) { int iter=0; *errstring=std::string(""); - std::vector f0,v,negative_f0; - std::vector > J; + std::vector f0,v; + std::vector > JJ; + Eigen::VectorXd r(x0.size()); + Eigen::Matrix2d J(x0.size(), x0.size()); double error = 999; while (iter==0 || std::abs(error)>tol){ f0 = f->call(x0); - J = f->Jacobian(x0); + JJ = f->Jacobian(x0); + + for (std::size_t i = 0; i < x0.size(); ++i) + { + r(i) = f0[i]; + for (std::size_t j = 0; j < x0.size(); ++j) + { + J(i,j) = JJ[i][j]; + } + } + + Eigen::Vector2d v = J.colPivHouseholderQr().solve(-r); - // Negate f0 - negative_f0 = f0; - for (std::size_t i = 0; imaxiter){ *errstring=std::string("reached maximum number of iterations"); diff --git a/wrappers/Python/CoolProp/AbstractState.pxd b/wrappers/Python/CoolProp/AbstractState.pxd index 26f76803..2cf693cd 100644 --- a/wrappers/Python/CoolProp/AbstractState.pxd +++ b/wrappers/Python/CoolProp/AbstractState.pxd @@ -71,6 +71,7 @@ cdef class AbstractState: cpdef double molar_mass(self) except * cpdef double acentric_factor(self) except* + cpdef tuple true_critical_point(self) cpdef double keyed_output(self, constants_header.parameters) except * cpdef double trivial_keyed_output(self, constants_header.parameters) except * cpdef double saturated_liquid_keyed_output(self, constants_header.parameters) except * diff --git a/wrappers/Python/CoolProp/AbstractState.pyx b/wrappers/Python/CoolProp/AbstractState.pyx index 2790e2c3..d381f060 100644 --- a/wrappers/Python/CoolProp/AbstractState.pyx +++ b/wrappers/Python/CoolProp/AbstractState.pyx @@ -188,6 +188,12 @@ cdef class AbstractState: cpdef mole_fractions_vapor(self): """ Get the mole fractions of the vapor phase - wrapper of c++ function :cpapi:`CoolProp::AbstractState::mole_fractions_vapor(void)` """ return self.thisptr.mole_fractions_vapor() + + cpdef tuple true_critical_point(self): + """ Get the "true" critical point where dp/drho|T = 0 & d2p/drho^2|T = 0 - wrapper of c++ function :cpapi:`CoolProp::AbstractState::true_critical_point(void)` """ + cdef double T = 1e99, rho = 1e99 + self.thisptr.true_critical_point(T, rho) + return T, rho ## ---------------------------------------- ## Derivatives diff --git a/wrappers/Python/CoolProp/cAbstractState.pxd b/wrappers/Python/CoolProp/cAbstractState.pxd index cc8ca8d3..82d76d91 100644 --- a/wrappers/Python/CoolProp/cAbstractState.pxd +++ b/wrappers/Python/CoolProp/cAbstractState.pxd @@ -100,6 +100,7 @@ cdef extern from "AbstractState.h" namespace "CoolProp": double first_two_phase_deriv(constants_header.parameters Of, constants_header.parameters Wrt, constants_header.parameters Constant) except+ValueError double second_two_phase_deriv(constants_header.parameters Of, constants_header.parameters Wrt1, constants_header.parameters Constant1, constants_header.parameters Wrt2, constants_header.parameters Constant2) except+ValueError double first_two_phase_deriv_splined(constants_header.parameters Of, constants_header.parameters Wrt, constants_header.parameters Constant, double x_end) except+ValueError + void true_critical_point(double &T, double &rho) except +ValueError double melting_line(int,int,double) except+ValueError bool has_melting_line() except+ValueError