diff --git a/dev/codelite/coolprop.project b/dev/codelite/coolprop.project index 3e1d88a0..28c86cc6 100644 --- a/dev/codelite/coolprop.project +++ b/dev/codelite/coolprop.project @@ -1,10 +1,10 @@ - + - + + diff --git a/include/CoolPropTools.h b/include/CoolPropTools.h index 2763e5d8..5cd4ed2c 100644 --- a/include/CoolPropTools.h +++ b/include/CoolPropTools.h @@ -334,4 +334,22 @@ return _HUGE; \ } + /// A spline is a curve given by the form y = ax^3 + bx^2 + c*x + d + /// As there are 4 constants, 4 constraints are needed to create the spline. These constraints could be the derivative or value at a point + /// Often, the value and derivative of the value are known at two points. + class SplineClass + { + protected: + int Nconstraints; + std::vector > A; + std::vector B; + public: + double a,b,c,d; + SplineClass(); + bool build(void); + bool add_value_constraint(double x, double y); + void add_4value_constraints(double x1, double x2, double x3, double x4, double y1, double y2, double y3, double y4); + bool add_derivative_constraint(double x, double dydx); + double evaluate(double x); + }; #endif diff --git a/include/PhaseEnvelope.h b/include/PhaseEnvelope.h index 58d73adc..e7f5bb68 100644 --- a/include/PhaseEnvelope.h +++ b/include/PhaseEnvelope.h @@ -5,14 +5,22 @@ namespace CoolProp{ +/** \brief A data structure to hold the data for a phase envelope + * + */ class PhaseEnvelopeData { public: + bool TypeI; ///< True if it is a Type-I mixture that has a phase envelope that looks like a pure fluid more or less bool built; ///< True if the phase envelope has been constructed + std::size_t iTsat_max, ///< The index of the point corresponding to the maximum temperature for Type-I mixtures + ipsat_max, ///< The index of the point corresponding to the maximum pressure for Type-I mixtures + icrit; ///< The index of the point corresponding to the critical point + std::vector< std::vector > K, lnK, x, y; std::vector T, p, lnT, lnp, rhomolar_liq, rhomolar_vap, lnrhomolar_liq, lnrhomolar_vap, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap; - PhaseEnvelopeData(){ built = false; }; + PhaseEnvelopeData(){ built = false; TypeI = false; }; void resize(std::size_t N) { @@ -26,6 +34,40 @@ public: lnrhomolar_liq.clear(); lnrhomolar_vap.clear(); hmolar_liq.clear(); hmolar_vap.clear(); smolar_liq.clear(); smolar_vap.clear(); K.clear(); lnK.clear(); x.clear(); y.clear(); } + void insert_variables(const long double T, + const long double p, + const long double rhomolar_liq, + const long double rhomolar_vap, + const long double hmolar_liq, + const long double hmolar_vap, + const long double smolar_liq, + const long double smolar_vap, + const std::vector & x, + const std::vector & y, + std::size_t i) + { + std::size_t N = K.size(); + if (N==0){throw CoolProp::ValueError("Cannot insert variables in phase envelope since resize() function has not been called");} + this->p.insert(this->p.begin() + i, p); + this->T.insert(this->T.begin() + i, T); + this->lnT.insert(this->lnT.begin() + i, log(T)); + this->lnp.insert(this->lnp.begin() + i, log(p)); + this->rhomolar_liq.insert(this->rhomolar_liq.begin() + i, rhomolar_liq); + this->rhomolar_vap.insert(this->rhomolar_vap.begin() + i, rhomolar_vap); + this->hmolar_liq.insert(this->hmolar_liq.begin() + i, hmolar_liq); + this->hmolar_vap.insert(this->hmolar_vap.begin() + i, hmolar_vap); + this->smolar_liq.insert(this->smolar_liq.begin() + i, smolar_liq); + this->smolar_vap.insert(this->smolar_vap.begin() + i, smolar_vap); + this->lnrhomolar_liq.insert(this->lnrhomolar_liq.begin() + i, log(rhomolar_liq)); + this->lnrhomolar_vap.insert(this->lnrhomolar_vap.begin() + i, log(rhomolar_vap)); + for (unsigned int j = 0; j < N; j++) + { + this->K[j].insert(this->K[j].begin() + i, y[j]/x[j]); + this->lnK[j].insert(this->lnK[j].begin() + i, log(y[j]/x[j])); + this->x[j].insert(this->x[j].begin() + i, x[j]); + this->y[j].insert(this->y[j].begin() + i, y[j]); + } + }; void store_variables(const long double T, const long double p, const long double rhomolar_liq, diff --git a/src/Backends/Helmholtz/Fluids/FluidLibrary.h b/src/Backends/Helmholtz/Fluids/FluidLibrary.h index 0be7ed2a..ab075272 100644 --- a/src/Backends/Helmholtz/Fluids/FluidLibrary.h +++ b/src/Backends/Helmholtz/Fluids/FluidLibrary.h @@ -50,6 +50,8 @@ protected: assert(n.size() == d.size()); assert(n.size() == t.size()); assert(n.size() == l.size()); + + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){for(std::size_t i = 0; i < n.size(); ++i){n[i] *= EOS.R_u/R_u_CODATA;}} EOS.alphar.GenExp.add_Power(n,d,t,l); } else if (!type.compare("ResidualHelmholtzGaussian")) @@ -67,6 +69,7 @@ protected: assert(n.size() == epsilon.size()); assert(n.size() == beta.size()); assert(n.size() == gamma.size()); + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){for(std::size_t i = 0; i < n.size(); ++i){n[i] *= EOS.R_u/R_u_CODATA;}} EOS.alphar.GenExp.add_Gaussian(n,d,t,eta,epsilon,beta,gamma); } else if (!type.compare("ResidualHelmholtzNonAnalytic")) @@ -87,6 +90,7 @@ protected: assert(n.size() == B.size()); assert(n.size() == C.size()); assert(n.size() == D.size()); + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){for(std::size_t i = 0; i < n.size(); ++i){n[i] *= EOS.R_u/R_u_CODATA;}} EOS.alphar.NonAnalytic = ResidualHelmholtzNonAnalytic(n,a,b,beta,A,B,C,D); } else if (!type.compare("ResidualHelmholtzLemmon2005")) @@ -100,6 +104,7 @@ protected: assert(n.size() == t.size()); assert(n.size() == l.size()); assert(n.size() == m.size()); + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){for(std::size_t i = 0; i < n.size(); ++i){n[i] *= EOS.R_u/R_u_CODATA;}} EOS.alphar.GenExp.add_Lemmon2005(n,d,t,l,m); } else if (!type.compare("ResidualHelmholtzExponential")) @@ -113,6 +118,7 @@ protected: assert(n.size() == t.size()); assert(n.size() == g.size()); assert(n.size() == l.size()); + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){for(std::size_t i = 0; i < n.size(); ++i){n[i] *= EOS.R_u/R_u_CODATA;}} EOS.alphar.GenExp.add_Exponential(n,d,t,g,l); } else if (!type.compare("ResidualHelmholtzAssociating")) @@ -368,6 +374,10 @@ protected: // Validate the equation of state that was just created EOS.validate(); + + if (get_config_bool(NORMALIZE_GAS_CONSTANTS)){ + EOS.R_u = R_u_CODATA; + } } /// Parse the list of possible equations of state diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 3abc8eee..913cf616 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -120,6 +120,8 @@ void HelmholtzEOSMixtureBackend::calc_phase_envelope(const std::string &type) PhaseEnvelope = PhaseEnvelopeData(); // Build the phase envelope PhaseEnvelopeRoutines::build(*this); + // Finalize the phase envelope + PhaseEnvelopeRoutines::finalize(*this); }; void HelmholtzEOSMixtureBackend::set_mixture_parameters() { diff --git a/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp b/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp index 65af99c4..3cfada43 100644 --- a/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp +++ b/src/Backends/Helmholtz/PhaseEnvelopeRoutines.cpp @@ -5,6 +5,7 @@ #include "VLERoutines.h" #include "PhaseEnvelopeRoutines.h" #include "PhaseEnvelope.h" +#include "CoolPropTools.h" namespace CoolProp{ @@ -57,7 +58,8 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS) if (failure_count > 5){ // Stop since we are stuck at a bad point - throw SolutionError("stuck"); + //throw SolutionError("stuck"); + return; } if (iter - iter0 > 0){ IO.rhomolar_vap *= factor;} @@ -123,6 +125,7 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS) NR.call(HEOS, IO.y, IO.x, IO); } catch(std::exception &e){ + std::cout << e.what() << std::endl; // Try again, but with a smaller step IO.rhomolar_vap /= factor; factor = 1 + (factor-1)/2; @@ -181,6 +184,131 @@ void PhaseEnvelopeRoutines::build(HelmholtzEOSMixtureBackend &HEOS) } } +void PhaseEnvelopeRoutines::finalize(HelmholtzEOSMixtureBackend &HEOS) +{ + enum maxima_points {PMAX_SAT = 0, TMAX_SAT = 1}; + std::size_t imax; // Index of the maximal temperature or pressure + + PhaseEnvelopeData &env = HEOS.PhaseEnvelope; + + // Find the index of the point with the highest temperature + std::size_t iTmax = std::distance(env.T.begin(), std::max_element(env.T.begin(), env.T.end())); + + // Find the index of the point with the highest pressure + std::size_t ipmax = std::distance(env.p.begin(), std::max_element(env.p.begin(), env.p.end())); + + // Determine if the phase envelope corresponds to a Type I mixture + // For now we consider a mixture to be Type I if the pressure at the + // end of the envelope is lower than max pressure pressure + env.TypeI = env.p[env.p.size()-1] < env.p[ipmax]; + + // Approximate solutions for the maxima of the phase envelope + // See method in Gernert. We use our spline class to find the coefficients + if (env.TypeI){ + for (int imaxima = 0; imaxima <= 1; ++imaxima){ + maxima_points maxima; + if (imaxima == PMAX_SAT){ + maxima = PMAX_SAT; + } + else if (imaxima == TMAX_SAT){ + maxima = TMAX_SAT; + } + + // Spline using the points around it + SplineClass spline; + if (maxima == TMAX_SAT){ + imax = iTmax; + spline.add_4value_constraints(env.rhomolar_vap[iTmax-1], env.rhomolar_vap[iTmax], env.rhomolar_vap[iTmax+1], env.rhomolar_vap[iTmax+2], + env.T[iTmax-1], env.T[iTmax], env.T[iTmax+1], env.T[iTmax+2] ); + } + else{ + imax = ipmax; + spline.add_4value_constraints(env.rhomolar_vap[ipmax-1], env.rhomolar_vap[ipmax], env.rhomolar_vap[ipmax+1], env.rhomolar_vap[ipmax+2], + env.p[ipmax-1], env.p[ipmax], env.p[ipmax+1], env.p[ipmax+2] ); + } + spline.build(); // y = a*rho^3 + b*rho^2 + c*rho + d + + // Take derivative + // dy/drho = 3*a*rho^2 + 2*b*rho + c + // Solve quadratic for derivative to find rho + int Nsoln = 0; double rho0 = _HUGE, rho1 = _HUGE, rho2 = _HUGE; + solve_cubic(0, 3*spline.a, 2*spline.b, spline.c, Nsoln, rho0, rho1, rho2); + + SaturationSolvers::newton_raphson_saturation_options IO; + IO.rhomolar_vap = _HUGE; + // Find the correct solution + if (Nsoln == 1){ + IO.rhomolar_vap = rho0; + } + else if (Nsoln == 2){ + if (is_in_closed_range(env.rhomolar_vap[imax-1], env.rhomolar_vap[imax+1], (long double)rho0)){ IO.rhomolar_vap = rho0; } + if (is_in_closed_range(env.rhomolar_vap[imax-1], env.rhomolar_vap[imax+1], (long double)rho1)){ IO.rhomolar_vap = rho1; } + } + else{ + throw ValueError("More than 2 solutions found"); + } + + class solver_resid : public FuncWrapper1D + { + public: + std::size_t imax; + maxima_points maxima; + HelmholtzEOSMixtureBackend *HEOS; + SaturationSolvers::newton_raphson_saturation NR; + SaturationSolvers::newton_raphson_saturation_options IO; + solver_resid(HelmholtzEOSMixtureBackend &HEOS, std::size_t imax, maxima_points maxima) + { + this->HEOS = &HEOS, this->imax = imax; this->maxima = maxima; + }; + double call(double rhomolar_vap){ + PhaseEnvelopeData &env = HEOS->PhaseEnvelope; + IO.bubble_point = false; + IO.rhomolar_vap = rhomolar_vap; + IO.y = HEOS->get_mole_fractions(); + IO.x = IO.y; // Just to give it good size + IO.T = CubicInterp(env.rhomolar_vap, env.T, imax-1, imax, imax+1, imax+2, IO.rhomolar_vap); + IO.rhomolar_liq = CubicInterp(env.rhomolar_vap, env.rhomolar_liq, imax-1, imax, imax+1, imax+2, IO.rhomolar_vap); + for (std::size_t i = 0; i < IO.x.size()-1; ++i) // First N-1 elements + { + IO.x[i] = CubicInterp(env.rhomolar_vap, env.x[i], imax-1, imax, imax+1, imax+2, IO.rhomolar_vap); + } + IO.x[IO.x.size()-1] = 1 - std::accumulate(IO.x.begin(), IO.x.end()-1, 0.0); + NR.call(*HEOS, IO.y, IO.x, IO); + if (maxima == TMAX_SAT){ + return NR.dTsat_dPsat; + } + else{ + return NR.dPsat_dTsat; + } + }; + }; + + solver_resid resid(HEOS, imax, maxima); + std::string errstr; + double rho = Brent(resid, IO.rhomolar_vap*0.95, IO.rhomolar_vap*1.05, DBL_EPSILON, 1e-12, 100, errstr); + + // If maxima point is greater than density at point from the phase envelope, increase index by 1 so that the + // insertion will happen *after* the point in the envelope since density is monotonically increasing. + if (rho > env.rhomolar_vap[imax]){ imax++; } + + env.insert_variables(resid.IO.T, resid.IO.p, resid.IO.rhomolar_liq, resid.IO.rhomolar_vap, resid.IO.hmolar_liq, + resid.IO.hmolar_vap, resid.IO.smolar_liq, resid.IO.smolar_vap, resid.IO.x, resid.IO.y, imax); + + if (maxima == TMAX_SAT){ + env.iTsat_max = imax; + if (imax <= env.ipsat_max){ + // Psat_max goes first; an insertion at a smaller index bumps up the index + // for ipsat_max, so we bump the index to keep pace + env.ipsat_max++; + } + } + else{ + env.ipsat_max = imax; + } + } + } +} + } /* namespace CoolProp */ #endif \ No newline at end of file diff --git a/src/Backends/Helmholtz/PhaseEnvelopeRoutines.h b/src/Backends/Helmholtz/PhaseEnvelopeRoutines.h index 8925b8f1..8a1cc4b0 100644 --- a/src/Backends/Helmholtz/PhaseEnvelopeRoutines.h +++ b/src/Backends/Helmholtz/PhaseEnvelopeRoutines.h @@ -5,6 +5,10 @@ namespace CoolProp{ class PhaseEnvelopeRoutines{ public: static void build(HelmholtzEOSMixtureBackend &HEOS); + + /** \brief Finalize the phase envelope and calculate maxima values, critical point, etc. + */ + static void finalize(HelmholtzEOSMixtureBackend &HEOS); }; } /* namespace CoolProp */ \ No newline at end of file diff --git a/src/Backends/Helmholtz/VLERoutines.cpp b/src/Backends/Helmholtz/VLERoutines.cpp index ac028622..142ae08f 100644 --- a/src/Backends/Helmholtz/VLERoutines.cpp +++ b/src/Backends/Helmholtz/VLERoutines.cpp @@ -1041,6 +1041,16 @@ void SaturationSolvers::newton_raphson_saturation::build_arrays() error_rms += r[i]*r[i]; // Sum the squares } error_rms = sqrt(error_rms); // Square-root (The R in RMS) + + // Calculate derivatives along phase boundary; + long double dQ_dPsat = 0, dQ_dTsat = 0; + for (std::size_t i = 0; i < N; ++i) + { + dQ_dPsat += x[i]*(MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(rSatV, i, xN_flag)); + dQ_dTsat += x[i]*(MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatV, i, xN_flag)); + } + dTsat_dPsat = -dQ_dPsat/dQ_dTsat; + dPsat_dTsat = -dQ_dTsat/dQ_dPsat; } } /* namespace CoolProp*/ diff --git a/src/Backends/Helmholtz/VLERoutines.h b/src/Backends/Helmholtz/VLERoutines.h index 1147cbf8..10868d32 100644 --- a/src/Backends/Helmholtz/VLERoutines.h +++ b/src/Backends/Helmholtz/VLERoutines.h @@ -225,6 +225,7 @@ namespace SaturationSolvers bool logging; bool bubble_point; int Nsteps; + long double dTsat_dPsat, dPsat_dTsat; STLMatrix J; HelmholtzEOSMixtureBackend *HEOS; std::vector K, x, y, phi_ij_liq, phi_ij_vap, dlnphi_drho_liq, dlnphi_drho_vap, r, negative_r, dXdS, neg_dFdS; diff --git a/src/CoolPropTools.cpp b/src/CoolPropTools.cpp index f69087a5..22fbba72 100644 --- a/src/CoolPropTools.cpp +++ b/src/CoolPropTools.cpp @@ -9,6 +9,8 @@ #include "float.h" #include #include "CoolPropTools.h" +#include "MatrixMath.h" +#include "Exceptions.h" double root_sum_square(std::vector x) { @@ -250,3 +252,64 @@ std::string strjoin(std::vector strings, std::string delim) } return output; } + + +SplineClass::SplineClass() +{ + Nconstraints = 0; + A.resize(4, std::vector(4, 0)); + B.resize(4,0); +} +bool SplineClass::build() +{ + if (Nconstraints == 4) + { + std::vector abcd = CoolProp::linsolve(A,B); + a = abcd[0]; + b = abcd[1]; + c = abcd[2]; + d = abcd[3]; + return true; + } + else + { + throw CoolProp::ValueError(format("Number of constraints[%d] is not equal to 4", Nconstraints)); + } +} +bool SplineClass::add_value_constraint(double x, double y) +{ + int i = Nconstraints; + if (i == 4) + return false; + A[i][0] = x*x*x; + A[i][1] = x*x; + A[i][2] = x; + A[i][3] = 1; + B[i] = y; + Nconstraints++; + return true; +} +void SplineClass::add_4value_constraints(double x1, double x2, double x3, double x4, double y1, double y2, double y3, double y4) +{ + add_value_constraint(x1, y1); + add_value_constraint(x2, y2); + add_value_constraint(x3, y3); + add_value_constraint(x4, y4); +} +bool SplineClass::add_derivative_constraint(double x, double dydx) +{ + int i = Nconstraints; + if (i == 4) + return false; + A[i][0] = 3*x*x; + A[i][1] = 2*x; + A[i][2] = 1; + A[i][3] = 0; + B[i] = dydx; + Nconstraints++; + return true; +} +double SplineClass::evaluate(double x) +{ + return a*x*x*x+b*x*x+c*x+d; +} \ No newline at end of file diff --git a/wrappers/Python/CoolProp/AbstractState.pxd b/wrappers/Python/CoolProp/AbstractState.pxd index 806e678a..3218a4f5 100644 --- a/wrappers/Python/CoolProp/AbstractState.pxd +++ b/wrappers/Python/CoolProp/AbstractState.pxd @@ -7,6 +7,8 @@ cimport cAbstractState cimport constants_header cdef class PyPhaseEnvelopeData: + cpdef public bool TypeI + cpdef public size_t iTsat_max, ipsat_max, icrit cpdef public list T, p, lnT, lnp, rhomolar_liq, rhomolar_vap, lnrhomolar_liq, lnrhomolar_vap, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap cdef class AbstractState: diff --git a/wrappers/Python/CoolProp/AbstractState.pyx b/wrappers/Python/CoolProp/AbstractState.pyx index 46c67a54..17990f4c 100644 --- a/wrappers/Python/CoolProp/AbstractState.pyx +++ b/wrappers/Python/CoolProp/AbstractState.pyx @@ -133,6 +133,9 @@ cdef class AbstractState: pe_out.rhomolar_vap = pe_data.rhomolar_vap pe_out.hmolar_liq = pe_data.hmolar_liq pe_out.hmolar_vap = pe_data.hmolar_vap - + print pe_data.iTsat_max, pe_data.ipsat_max + pe_out.iTsat_max = pe_data.iTsat_max + pe_out.ipsat_max = pe_data.ipsat_max + pe_out.TypeI = pe_data.TypeI return pe_out \ No newline at end of file diff --git a/wrappers/Python/CoolProp/cAbstractState.pxd b/wrappers/Python/CoolProp/cAbstractState.pxd index 0588084a..40d2a813 100644 --- a/wrappers/Python/CoolProp/cAbstractState.pxd +++ b/wrappers/Python/CoolProp/cAbstractState.pxd @@ -6,6 +6,8 @@ cimport constants_header cdef extern from "PhaseEnvelope.h" namespace "CoolProp": cdef cppclass PhaseEnvelopeData: + bool TypeI + size_t iTsat_max, ipsat_max, icrit vector[long double] T, p, lnT, lnp, rhomolar_liq, rhomolar_vap, lnrhomolar_liq, lnrhomolar_vap, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap cdef extern from "AbstractState.h" namespace "CoolProp":