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

Modifications to saturation routines (work great now)

Browse Source 
Including:
adding critical pressure parameter
Adding relaxation parameter for NR sat_T
Only evaluate correction step if not near critical temp
Fixed bug with fluids like "n-Propane" that was getting interpreted as an incompressible fluid

Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-08-03 23:36:22 +02:00
parent e9ac0bcc9d
commit 17458261b5
7 changed files with 39 additions and 25 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
include/CoolPropFluid.h
View File

@@ -395,7 +395,12 @@ public:
solver_resid resid(this, value);
std::string errstring;
return Brent(resid,Tmin,Tmax,DBL_EPSILON,1e-12,100,errstring);
try{
return Brent(resid,Tmin,Tmax,DBL_EPSILON,1e-12,100,errstring);
}
catch(std::exception &e){
return Secant(resid,Tmax, -0.01, 1e-12, 100, errstring);
}
}
};

2
include/DataStructures.h
View File

@@ -31,7 +31,7 @@ enum parameters{
// General parameters
imolar_mass, irhomolar_reducing, irhomolar_critical, iT_reducing, iT_critical,
irhomass_reducing, irhomass_critical,
irhomass_reducing, irhomass_critical, iP_critical,
// Bulk properties
iT, iP, iQ, iTau, iDelta,

6
src/AbstractState.cpp
View File

@@ -167,6 +167,8 @@ double AbstractState::trivial_keyed_output(int key)
return get_reducing().T;
case irhomolar_reducing:
return get_reducing().rhomolar;
case iP_critical:
return this->p_critical();
case iT_critical:
return this->T_critical();
case irhomolar_critical:
@@ -220,10 +222,6 @@ double AbstractState::keyed_output(int key)
return get_reducing().rhomolar;
case ispeed_sound:
return speed_sound();
//case iT_critical:
// return get_critical().T;
//case irhomolar_critical:
// return get_critical().rhomolar; // TODO
case ialpha0:
return alpha0();
case idalpha0_ddelta_consttau:

3
src/Backends/Helmholtz/FlashRoutines.cpp
View File

@@ -28,7 +28,7 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
{
// Set some imput options
SaturationSolvers::saturation_T_pure_options options;
options.omega = 1.0;
options.omega = 0.1;
options.use_guesses = false;
// Actually call the solver
SaturationSolvers::saturation_T_pure(&HEOS, HEOS._T, options);
@@ -114,6 +114,7 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
options.specified_variable = SaturationSolvers::saturation_PHSU_pure_options::IMPOSED_PL;
// Use logarithm of delta as independent variables
options.use_logdelta = false;
options.omega = 0.3;
// Actually call the solver
SaturationSolvers::saturation_PHSU_pure(&HEOS, HEOS._p, options);

39
src/Backends/Helmholtz/VLERoutines.cpp
View File

@@ -35,12 +35,19 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l
{
// Invert liquid density ancillary to get temperature
// TODO: fit inverse ancillaries too
T = HEOS->get_components()[0]->ancillaries.pL.invert(specified_value);
try{
T = HEOS->get_components()[0]->ancillaries.pL.invert(specified_value);
}
catch(std::exception &e)
{
throw ValueError("Unable to invert ancillary equation");
}
}
else
{
throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid",options.specified_variable));
}
if (T > HEOS->T_critical()-1){ T -= 1; }
// Evaluate densities from the ancillary equations
rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T);
@@ -391,7 +398,7 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE
long double rhoL,rhoV,JL,JV,KL,KV,dJL,dJV,dKL,dKV;
long double DELTA, deltaL=0, deltaV=0, tau, error, PL, PV, stepL, stepV;
int iter=0;
// Use the density ancillary function as the starting point for the solver
try
{
if (options.use_guesses)
@@ -402,25 +409,27 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE
}
else
{
// Use the density ancillary function as the starting point for the solver
// If very close to the critical temp, evaluate the ancillaries for a slightly lower temperature
if (T > 0.99*HEOS->get_reducing().T){
rhoL = HEOS->get_components()[0]->ancillaries.rhoL.evaluate(T-1);
rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T-1);
rhoL = HEOS->get_components()[0]->ancillaries.rhoL.evaluate(T-0.1);
rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T-0.1);
}
else{
rhoL = HEOS->get_components()[0]->ancillaries.rhoL.evaluate(T);
rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T);
// Apply a single step of Newton's method to improve guess value for liquid
// based on the error between the gas pressure (which is usually very close already)
// and the liquid pressure, which can sometimes (especially at low pressure),
// be way off, and often times negative
SatL->update(DmolarT_INPUTS, rhoL, T);
SatV->update(DmolarT_INPUTS, rhoV, T);
rhoL += -(SatL->p()-SatV->p())/SatL->first_partial_deriv(iP, iDmolar, iT);
}
}
// Apply a single step of Newton's method to improve guess value for liquid
// based on the error between the gas pressure (which is usually very close already)
// and the liquid pressure, which can sometimes (especially at low pressure),
// be way off, and often times negative
SatL->update(DmolarT_INPUTS, rhoL, T);
SatV->update(DmolarT_INPUTS, rhoV, T);
rhoL += -(SatL->p()-SatV->p())/SatL->first_partial_deriv(iP, iDmolar, iT);
deltaL = rhoL/reduce.rhomolar;
deltaV = rhoV/reduce.rhomolar;
tau = reduce.T/T;
@@ -474,7 +483,7 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE
DELTA = dJV*dKL-dJL*dKV;
error = fabs(KL-KV)+fabs(JL-JV);
error = sqrt((KL-KV)*(KL-KV)+(JL-JV)*(JL-JV));
// Get the predicted step
stepL = options.omega/DELTA*( (KV-KL)*dJV-(JV-JL)*dKV);
@@ -485,14 +494,14 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE
rhoL = deltaL*reduce.rhomolar; rhoV = deltaV*reduce.rhomolar;
}
else{
throw ValueError(format("rhosatPure_Akasaka failed"));
throw ValueError(format("rhosatPure_Akasaka densities are crossed"));
}
iter++;
if (iter > 100){
throw SolutionError(format("Akasaka solver did not converge after 100 iterations"));
}
}
while (error > 1e-10 && fabs(stepL) > 10*DBL_EPSILON*fabs(stepL) && fabs(stepV) > 10*DBL_EPSILON*fabs(stepV));
while (error > 1e-12 && fabs(stepL) > 10*DBL_EPSILON*fabs(stepL) && fabs(stepV) > 10*DBL_EPSILON*fabs(stepV));
double p_error_limit = 1e-3;
double p_error = (PL - PV)/PL;

3
src/CoolProp.cpp
View File

@@ -295,7 +295,7 @@ bool has_solution_concentration(const std::string &fluid_string)
// Start at the "::" if it is found to chop off the delimiter between backend and fluid
std::size_t i = fluid_string.find("::");
// If can find "-", expect mass fractions encoded as string
if (fluid_string.find('-',i+2)!=std::string::npos)
if (fluid_string.find('-',i+2)!=std::string::npos && fluid_string.find("INCOMP")==0)
return true;
return false;
}
@@ -419,7 +419,6 @@ double PropsSI(const std::string &Output, const std::string &Name1, double Prop1
// Extract the fractions and reformulate the list of fluids INCOMP::EG-0.2 -> INCOMP::EG and [0.2]
std::string Ref2 = extract_concentrations(Ref, fractions);
return _PropsSI(Output, Name1, Prop1, Name2, Prop2, Ref2, fractions);
}
else
{

4
src/DataStructures.cpp
View File

@@ -49,9 +49,10 @@ parameter_info parameter_info_list[] = {
parameter_info(imolar_mass, "molar_mass","O","kg/mol","Molar mass",true),
parameter_info(irhomolar_reducing, "rhomolar_reducing","O","mol/m^3","Molar density at reducing point",true),
parameter_info(irhomolar_critical, "rhomolar_critical","O","mol/m^3","Molar density at critical point",true),
parameter_info(irhomass_critical, "rhomass_critical","O","kg/m^3","Molar density at critical point",true),
parameter_info(irhomass_critical, "rhomass_critical","O","kg/m^3","Mass density at critical point",true),
parameter_info(iT_reducing, "T_reducing","O","K","Temperature at the reducing point",true),
parameter_info(iT_critical, "T_critical","O","K","Temperature at the critical point",true),
parameter_info(iP_critical, "p_critical","O","Pa","Pressure at the critical point",true),
parameter_info(iisothermal_compressibility, "isothermal_compressibility","O","1/Pa","Isothermal compressibility",false),
parameter_info(ispeed_sound, "speed_of_sound","O","m/s","Speed of sound",false),
parameter_info(iviscosity, "viscosity","O","Pa-s","Viscosity",false),
@@ -86,6 +87,7 @@ public:
index_map.insert(std::pair<std::string, int>("O", iCvmass));
index_map.insert(std::pair<std::string, int>("V", iviscosity));
index_map.insert(std::pair<std::string, int>("L", iconductivity));
index_map.insert(std::pair<std::string, int>("pcrit", iP_critical));
index_map.insert(std::pair<std::string, int>("Tcrit", iT_critical));
index_map.insert(std::pair<std::string, int>("rhocrit", irhomass_critical));
}