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

Potentially fixed critical spline code for linear and quadratic "splines". Closes https://github.com/CoolProp/CoolProp/issues/111

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-08-26 19:31:25 +02:00
parent 5fea309764
commit 05f4c64049
2 changed files with 51 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

23
include/CoolPropFluid.h
View File

@@ -46,29 +46,32 @@ struct CriticalRegionSplines{
bool enabled;
CriticalRegionSplines(){enabled = false;};
void get_densities(double T, double rho_min, double rho_crit, double rho_max, double &rhoL, double &rhoV){
bool ok1 = false, ok2 = false, ok3 = false;
int Nsoln = -1, Ngood = 0;
double rho1, rho2, rho3;
double rho1 =0, rho2 = 0, rho3 = 0;
// -----------
// Liquid part
// -----------
Ngood = 0;
long double v1 = cL[0]*pow(rho_crit,3) + cL[1]*pow(rho_crit,2) + cL[2]*pow(rho_crit,1) + cL[3];
long double v2 = cL[0]*pow(rho_max,3) + cL[1]*pow(rho_max,2) + cL[2]*pow(rho_max,1) + cL[3];
solve_cubic(cL[0], cL[1], cL[2], cL[3]-T, Nsoln, rho1, rho2, rho3);
if (rho1 < rho_max && rho1 > rho_crit){ ok1 = true; Ngood++; rhoL = rho1; }
if (rho2 < rho_max && rho2 > rho_crit){ ok2 = true; Ngood++; rhoL = rho2; }
if (rho3 < rho_max && rho3 > rho_crit){ ok3 = true; Ngood++; rhoL = rho3; }
if (Ngood != 1){ throw ValueError(format("more than one liquid solution found for critical spline for T=%g",T));};
if (rho1 < rho_max && rho1 > rho_crit){ Ngood++; rhoL = rho1; }
if (rho2 < rho_max && rho2 > rho_crit){ Ngood++; rhoL = rho2; }
if (rho3 < rho_max && rho3 > rho_crit){ Ngood++; rhoL = rho3; }
if (Ngood > 1){ throw ValueError(format("More than one liquid solution found for critical spline for T=%g",T));};
if (Ngood < 1){ throw ValueError(format("No liquid solution found for critical spline for T=%g",T));};
// ----------
// Vapor part
// ----------
Ngood = 0;
solve_cubic(cV[0], cV[1], cV[2], cV[3]-T, Nsoln, rho1, rho2, rho3);
if (rho1 > rho_min && rho1 < rho_crit){ ok1 = true; Ngood++; rhoV = rho1; }
if (rho2 > rho_min && rho2 < rho_crit){ ok2 = true; Ngood++; rhoV = rho2; }
if (rho3 > rho_min && rho3 < rho_crit){ ok3 = true; Ngood++; rhoV = rho3; }
if (Ngood != 1){ throw ValueError(format("more than one vapor solution found for critical spline for T=%g",T));};
if (rho1 > rho_min && rho1 < rho_crit){ Ngood++; rhoV = rho1; }
if (rho2 > rho_min && rho2 < rho_crit){ Ngood++; rhoV = rho2; }
if (rho3 > rho_min && rho3 < rho_crit){ Ngood++; rhoV = rho3; }
if (Ngood > 1){ throw ValueError(format("More than one vapor solution found for critical spline for T=%g",T));};
if (Ngood < 1){ throw ValueError(format("No vapor solution found for critical spline for T=%g",T));};
};
};

38
src/CoolPropTools.cpp
View File

@@ -161,9 +161,47 @@ std::string get_file_contents(const char *filename)
void solve_cubic(double a, double b, double c, double d, int &N, double &x0, double &x1, double &x2)
{
// 0 = ax^3 + b*x^2 + c*x + d
if (std::abs(a) < 10*DBL_EPSILON){
if (std::abs(b) < 10*DBL_EPSILON){
// Linear solution if a = 0 and b = 0
x0 = -d/c;
N = 1;
return;
}
else{
// Quadratic solution(s) if a = 0 and b != 0
x0 = (-c+sqrt(c*c-4*b*d))/(2*b);
x1 = (-c-sqrt(c*c-4*b*d))/(2*b);
N = 2;
return;
}
}
// Ok, it is really a cubic
// Discriminant
double DELTA = 18*a*b*c*d-4*b*b*b*d+b*b*c*c-4*a*c*c*c-27*a*a*d*d;
double DELTA0 = b*b*b-2*a*c;
// Deal with special cases
if (std::abs(DELTA) < 10*DBL_EPSILON){
if (std::abs(DELTA0)>0){
x0 = (9*a*d-b*c)/(2*DELTA0);
x1 = x0;
x2 = (4*a*b*c - 9*a*a*d - b*b*b)/(a*DELTA0);
N = 3;
return;
}
else{
x0 = -b/(3*a);
x1 = x0;
x2 = x0;
N = 3;
return;
}
}
// Coefficients for the depressed cubic t^3+p*t+q = 0
double p = (3*a*c-b*b)/(3*a*a);