Reimplemented critical region splines all the way from scripts to build to C++ code

Closes https://github.com/CoolProp/CoolProp/issues/94

Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-08-22 22:14:48 +02:00
parent cee9801423
commit 69bbdc693b
73 changed files with 1747 additions and 15 deletions

View File

@@ -37,12 +37,33 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
HEOS.calc_Tmin_sat(Tmin_satL, Tmin_satV);
Tmin_sat = std::max(Tmin_satL, Tmin_satV);
// Get a reference to keep the code a bit cleaner
CriticalRegionSplines &splines = HEOS.components[0]->pEOS->critical_region_splines;
// Check limits
if (!is_in_closed_range(Tmin_sat, Tmax_sat, static_cast<long double>(HEOS._T))){
throw ValueError(format("Temperature to QT_flash [%6g K] must be in range [%8g K, %8g K]",HEOS._T, Tmin_sat, Tmax_sat));
}
if (!(HEOS.components[0]->pEOS->pseudo_pure))
//splines.enabled = false;
// If exactly at the critical temperature, liquid and vapor have the critial density
if (std::abs(T-HEOS.T_critical())< 1e-14){
HEOS.SatL->update(DmolarT_INPUTS, HEOS.rhomolar_critical(), HEOS._T);
HEOS.SatV->update(DmolarT_INPUTS, HEOS.rhomolar_critical(), HEOS._T);
HEOS._rhomolar = HEOS.rhomolar_critical();
HEOS._p = HEOS.SatL->p();
}
else if (splines.enabled && HEOS._T > splines.T_min){
double rhoL, rhoV;
// Use critical region spline if it has it and temperature is in its range
splines.get_densities(T, splines.rhomolar_min, HEOS.rhomolar_critical(), splines.rhomolar_max, rhoL, rhoV);
HEOS.SatL->update(DmolarT_INPUTS, rhoL, HEOS._T);
HEOS.SatV->update(DmolarT_INPUTS, rhoV, HEOS._T);
HEOS._p = HEOS._Q*HEOS.SatV->p() + (1- HEOS._Q)*HEOS.SatL->p();
HEOS._rhomolar = 1/(HEOS._Q/HEOS.SatV->rhomolar() + (1 - HEOS._Q)/HEOS.SatL->rhomolar());
}
else if (!(HEOS.components[0]->pEOS->pseudo_pure))
{
// Set some imput options
SaturationSolvers::saturation_T_pure_options options;
@@ -73,11 +94,11 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
SaturationSolvers::saturation_T_pure_1D_P(&HEOS, T, options);
}
catch(std::exception &){
SaturationSolvers::saturation_critical(&HEOS, iT, T);
throw;
//SaturationSolvers::saturation_critical(&HEOS, iT, T);
}
}
// Load the outputs
HEOS._phase = iphase_twophase;
HEOS._p = HEOS._Q*HEOS.SatV->p() + (1- HEOS._Q)*HEOS.SatL->p();
HEOS._rhomolar = 1/(HEOS._Q/HEOS.SatV->rhomolar() + (1 - HEOS._Q)/HEOS.SatL->rhomolar());
}
@@ -108,12 +129,13 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
rhoLsat = rhoLanc;
rhoVsat = rhoVanc;
}
HEOS._phase = iphase_twophase;
HEOS._p = HEOS._Q*psatVanc + (1-HEOS._Q)*psatLanc;
HEOS._rhomolar = 1/(HEOS._Q/rhoVsat + (1-HEOS._Q)/rhoLsat);
HEOS.SatL->update(DmolarT_INPUTS, rhoLsat, HEOS._T);
HEOS.SatV->update(DmolarT_INPUTS, rhoVsat, HEOS._T);
}
// Load the outputs
HEOS._phase = iphase_twophase;
}
else
{
@@ -221,7 +243,7 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
double Tc = HEOS.components[0]->pEOS->reduce.T;
double Tt = HEOS.components[0]->pEOS->Ttriple;
double pt = HEOS.components[0]->pEOS->ptriple;
double Tsat_guess = 1/(1/Tc-(1/Tt-1/Tc)/log(pc/pt)*log(HEOS._p/pc));
//double Tsat_guess = 1/(1/Tc-(1/Tt-1/Tc)/log(pc/pt)*log(HEOS._p/pc));
// Set some imput options
SaturationSolvers::mixture_VLE_IO io;

View File

@@ -352,6 +352,16 @@ protected:
EOS.max_sat_p.rhomolar = cpjson::get_double(s, "rhomolar");
}
if (EOS_json.HasMember("critical_region_splines")){
rapidjson::Value &spline = EOS_json["critical_region_splines"];
EOS.critical_region_splines.T_min = cpjson::get_double(spline, "T_min");
EOS.critical_region_splines.T_max = cpjson::get_double(spline, "T_max");
EOS.critical_region_splines.rhomolar_min = cpjson::get_double(spline, "rhomolar_min");
EOS.critical_region_splines.rhomolar_max = cpjson::get_double(spline, "rhomolar_max");
EOS.critical_region_splines.cL = cpjson::get_double_array(spline["cL"]);
EOS.critical_region_splines.cV = cpjson::get_double_array(spline["cV"]);
EOS.critical_region_splines.enabled = true;
}
// Validate the equation of state that was just created
EOS.validate();