github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-01-22 04:18:02 -05:00
Code Issues Packages Projects Releases Wiki Activity

Fix some edge cases in critical point calculation

Browse Source
This commit is contained in:
Ian Bell
2016-02-24 21:44:41 -07:00
parent 250feab6fb
commit b2506e404b
3 changed files with 21 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/Backends/Cubics/CubicBackend.cpp
View File

@@ -38,8 +38,8 @@ void CoolProp::AbstractCubicBackend::get_critical_point_search_radii(double &R_d
// Now we scale them to get the appropriate search radii
double Tr_GERGlike, rhor_GERGlike;
get_linear_reducing_parameters(rhor_GERGlike, Tr_GERGlike);
R_delta *= rhor_GERGlike/rhomolar_reducing();
R_tau *= T_reducing()/Tr_GERGlike;
R_delta *= rhor_GERGlike/rhomolar_reducing()*5;
R_tau *= T_reducing()/Tr_GERGlike*5;
}
bool CoolProp::AbstractCubicBackend::get_critical_is_terminated(double &delta, double &tau)

16
src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
View File

@@ -3117,7 +3117,7 @@ CoolProp::CriticalState HelmholtzEOSMixtureBackend::calc_critical_point(double r
tau_delta[0] = T_reducing()/T0;
tau_delta[1] = rho0/rhomolar_reducing();
std::string errstr;
x = NDNewtonRaphson_Jacobian(&resid, tau_delta, 1e-11, 100, &errstr);
x = NDNewtonRaphson_Jacobian(&resid, tau_delta, 1e-10, 100, &errstr);
_critical.T = T_reducing()/x[0];
_critical.rhomolar = x[1]*rhomolar_reducing();
_critical.p = calc_pressure_nocache(_critical.T, _critical.rhomolar);
@@ -3247,11 +3247,21 @@ public:
// have an analytic solution for the derivative
R_tau = R_tau_tracer; R_delta = R_delta_tracer;
theta = Newton(this, theta_last, 1e-10, 100, errstr);
// If the solver takes a U-turn, going in the opposite direction of travel
// this is not a good thing, and force a Brent's method solver call to make sure we keep
// tracing in the same direction
if (std::abs(angle_difference(theta, theta_last)) > M_PI/2.0){
theta = Brent(this, theta_last-M_PI/3.5, theta_last+M_PI/3.5, DBL_EPSILON, 1e-10, 100, errstr);
// We have found at least one critical point and we are now well above the density
// associated with the first critical point
if (N_critical_points > 0 && delta > 2.5*critical_points[0].rhomolar/HEOS.rhomolar_reducing()){
// Stopping the search - probably we have a kink in the L1*=0 contour
// caused by poor low-temperature behavior
continue;
}
else{
theta = Brent(this, theta_last-M_PI/3.5, theta_last+M_PI/3.5, DBL_EPSILON, 1e-10, 100, errstr);
}
}
}
@@ -3272,7 +3282,7 @@ public:
// you have bracketed a critical point, run the full critical point solver to
// find the critical point and store it
if (i > 0 && M1*M1_last < 0){
double rhomolar = HEOS.rhomolar_reducing()*delta_new, T = HEOS.T_reducing()/tau_new;
double rhomolar = HEOS.rhomolar_reducing()*(delta+delta_new)/2.0, T = HEOS.T_reducing()/((tau+tau_new)/2.0);
CoolProp::CriticalState crit = HEOS.calc_critical_point(rhomolar, T);
critical_points.push_back(crit);
N_critical_points++;

8
src/Solvers.cpp
View File

@@ -75,6 +75,11 @@ std::vector<double> NDNewtonRaphson_Jacobian(FuncWrapperND *f, std::vector<doubl
// Update the guess
for (std::size_t i = 0; i<x0.size(); i++){ x0[i] += v(i);}
// Stop if the solution is not changing by more than numerical precision
if (v.cwiseAbs().maxCoeff() < DBL_EPSILON*100){
return x0;
}
error = root_sum_square(f0);
if (iter>maxiter){
*errstring=std::string("reached maximum number of iterations");
@@ -112,8 +117,7 @@ double Newton(FuncWrapper1DWithDeriv* f, double x0, double ftol, int maxiter, st
x += dx;
if (std::abs(dx/x) < 10*DBL_EPSILON)
{
if (std::abs(dx/x) < 1e-11){
return x;
}