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

Update parsing of the fluid JSON files for the new state architecture

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-06-06 12:05:55 +02:00
parent 7896845a5d
commit 6e287e9454
5 changed files with 196 additions and 189 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -52,7 +52,7 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
rhoLsat = HEOS.solver_rho_Tp(HEOS._T, psatLanc, rhoLanc);
rhoVsat = HEOS.solver_rho_Tp(HEOS._T, psatVanc, rhoLanc);
if (!ValidNumber(rhoLsat) || !ValidNumber(rhoVsat) ||
if (!ValidNumber(rhoLsat) || !ValidNumber(rhoVsat) ||
fabs(rhoLsat/rhoLanc-1) > 0.1 || fabs(rhoVanc/rhoVsat-1) > 0.1)
{
throw ValueError("pseudo-pure failed");
@@ -79,7 +79,7 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
// Use Wilson iteration to obtain updated guess for pressure
pguess = SaturationSolvers::saturation_Wilson(&HEOS, HEOS._Q, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions, pguess);
// Actually call the successive substitution solver
SaturationSolvers::successive_substitution(&HEOS, HEOS._Q, HEOS._T, pguess, HEOS.mole_fractions, HEOS.K, options);
}
@@ -109,7 +109,7 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
options.use_logdelta = false;
// Actually call the solver
SaturationSolvers::saturation_PHSU_pure(&HEOS, HEOS._p, options);
// Load the outputs
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());
@@ -128,10 +128,10 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
// Use Wilson iteration to obtain updated guess for temperature
Tguess = SaturationSolvers::saturation_Wilson(&HEOS, HEOS._Q, HEOS._p, SaturationSolvers::imposed_p, HEOS.mole_fractions, Tguess);
// Actually call the successive substitution solver
SaturationSolvers::successive_substitution(&HEOS, HEOS._Q, Tguess, HEOS._p, HEOS.mole_fractions, HEOS.K, io);
PhaseEnvelope::PhaseEnvelope_GV ENV_GV;
ENV_GV.build(&HEOS, HEOS.mole_fractions, HEOS.K, io);
}
@@ -143,7 +143,7 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
class solver_resid : public FuncWrapper1D
{
public:
HelmholtzEOSMixtureBackend *HEOS;
long double r, eos, rhomolar, value, T;
int other;
@@ -164,15 +164,15 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
default:
throw ValueError(format("Input not supported"));
}
r = eos - value;
return r;
};
};
std::string errstring;
if (HEOS.imposed_phase_index > -1)
if (HEOS.imposed_phase_index > -1)
{
// Use the phase defined by the imposed phase
HEOS._phase = HEOS.imposed_phase_index;
@@ -182,17 +182,17 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
if (HEOS.is_pure_or_pseudopure)
{
CoolPropFluid * component = HEOS.components[0];
shared_ptr<HelmholtzEOSMixtureBackend> Sat;
long double rhoLtriple = component->pEOS->rhoLtriple;
long double rhoVtriple = component->pEOS->rhoVtriple;
long double rhoLtriple = component->triple_liquid.rhomolar;
long double rhoVtriple = component->triple_vapor.rhomolar;
// Check if in the "normal" region
if (HEOS._rhomolar >= rhoVtriple && HEOS._rhomolar <= rhoLtriple)
{
long double yL, yV, value, y_solid;
long double TLtriple = component->pEOS->Ttriple; //TODO: separate TL and TV for ppure
long double TVtriple = component->pEOS->Ttriple; //TODO: separate TL and TV for ppure
long double TLtriple = component->triple_liquid.T; ///TODO: separate TL and TV for ppure
long double TVtriple = component->triple_vapor.T;
// First check if solid (below the line connecting the triple point values) - this is an error for now
switch (other)
{
@@ -246,7 +246,7 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
}
// Check if vapor/solid region below triple point vapor density
else if (HEOS._rhomolar < component->pEOS->rhoVtriple)
else if (HEOS._rhomolar < component->triple_vapor.rhomolar)
{
long double y, value;
long double TVtriple = component->pEOS->Ttriple; //TODO: separate TL and TV for ppure
@@ -280,7 +280,7 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
}
// Check in the liquid/solid region above the triple point density
else
else
{
long double y, value;
long double TLtriple = component->pEOS->Ttriple; //TODO: separate TL and TV for ppure
@@ -313,7 +313,7 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
}
}
}
else
else
throw NotImplementedError("PHSU_D_flash not ready for mixtures");
}
}
@@ -323,7 +323,7 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
}
void FlashRoutines::DHSU_T_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
{
if (HEOS.imposed_phase_index > -1)
if (HEOS.imposed_phase_index > -1)
{
// Use the phase defined by the imposed phase
HEOS._phase = HEOS.imposed_phase_index;
@@ -370,9 +370,9 @@ void FlashRoutines::DHSU_T_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
default:
break;
}
HEOS.calc_pressure();
HEOS.calc_pressure();
HEOS._Q = -1;
}
}
} /* namespace CoolProp */
} /* namespace CoolProp */

194
src/Backends/Helmholtz/Fluids/FluidLibrary.h
View File

@@ -16,8 +16,8 @@ extern int get_debug_level();
/// A container for the fluid parameters for the CoolProp fluids
/**
This container holds copies of all of the fluid instances for the fluids that are loaded in CoolProp.
New fluids can be added by passing in a rapidjson::Value instance to the add_one function, or
This container holds copies of all of the fluid instances for the fluids that are loaded in CoolProp.
New fluids can be added by passing in a rapidjson::Value instance to the add_one function, or
a rapidjson array of fluids to the add_many function.
*/
class JSONFluidLibrary
@@ -31,9 +31,9 @@ protected:
/// Parse the contributions to the residual Helmholtz energy
void parse_alphar(rapidjson::Value &alphar, EquationOfState &EOS)
{
{
for (rapidjson::Value::ValueIterator itr = alphar.Begin(); itr != alphar.End(); ++itr)
{
{
// A reference for code cleanness
rapidjson::Value &contribution = *itr;
@@ -140,7 +140,7 @@ protected:
{
if (!alpha0.IsArray()){throw ValueError();}
for (rapidjson::Value::ValueIterator itr = alpha0.Begin(); itr != alpha0.End(); ++itr)
{
{
// A reference for code cleanness
rapidjson::Value &contribution = *itr;
@@ -172,7 +172,7 @@ protected:
// Retrieve the values
std::vector<long double> n = cpjson::get_long_double_array(contribution["n"]);
std::vector<long double> t = cpjson::get_long_double_array(contribution["t"]);
std::vector<long double> c = cpjson::get_long_double_array(contribution["c"]);
std::vector<long double> d = cpjson::get_long_double_array(contribution["d"]);
if (EOS.alpha0.PlanckEinstein.is_enabled() == true){
@@ -217,7 +217,7 @@ protected:
}
else if (!type.compare("IdealGasHelmholtzCP0AlyLee"))
{
std::vector<long double> constants = cpjson::get_long_double_array(contribution["c"]);
long double Tc = cpjson::get_double(contribution, "Tc");
long double T0 = cpjson::get_double(contribution, "T0");
@@ -249,7 +249,7 @@ protected:
c.push_back(1);
d.push_back(1);
}
if (EOS.alpha0.PlanckEinstein.is_enabled() == true){
EOS.alpha0.PlanckEinstein.extend(n, t, c, d);
}
@@ -297,38 +297,41 @@ protected:
EOS.R_u = cpjson::get_double(EOS_json,"gas_constant");
EOS.molar_mass = cpjson::get_double(EOS_json,"molar_mass");
EOS.accentric = cpjson::get_double(EOS_json,"accentric");
EOS.Ttriple = cpjson::get_double(EOS_json, "Ttriple");
EOS.ptriple = cpjson::get_double(EOS_json, "ptriple");
EOS.rhoLtriple = cpjson::get_double(EOS_json, "rhoLtriple");
EOS.rhoVtriple = cpjson::get_double(EOS_json, "rhoVtriple");
EOS.pseudo_pure = cpjson::get_bool(EOS_json, "pseudo_pure");
EOS.limits.Tmax = cpjson::get_double(EOS_json, "T_max");
EOS.limits.pmax = cpjson::get_double(EOS_json, "p_max");
rapidjson::Value &reducing_state = EOS_json["reducing_state"];
rapidjson::Value &reducing_state = EOS_json["STATES"]["reducing"];
rapidjson::Value &satminL_state = EOS_json["STATES"]["sat_min_liquid"];
rapidjson::Value &satminV_state = EOS_json["STATES"]["sat_min_vapor"];
// Reducing state
EOS.reduce.T = cpjson::get_double(reducing_state,"T");
EOS.reduce.rhomolar = cpjson::get_double(reducing_state,"rhomolar");
EOS.reduce.p = cpjson::get_double(reducing_state,"p");
/// todo: define limits of EOS better
EOS.limits.Tmin = cpjson::get_double(satminL_state, "T");
EOS.Ttriple = EOS.limits.Tmin;
// BibTex keys
EOS.BibTeX_EOS = cpjson::get_string(EOS_json,"BibTeX_EOS");
EOS.BibTeX_CP0 = cpjson::get_string(EOS_json,"BibTeX_CP0");
parse_alphar(EOS_json["alphar"], EOS);
parse_alpha0(EOS_json["alpha0"], EOS);
// Validate the equation of state that was just created
EOS.validate();
}
/// Parse the list of possible equations of state
void parse_EOS_listing(rapidjson::Value &EOS_array, CoolPropFluid & fluid)
{
for (rapidjson::Value::ValueIterator itr = EOS_array.Begin(); itr != EOS_array.End(); ++itr)
{
{
parse_EOS(*itr,fluid);
}
@@ -436,7 +439,7 @@ protected:
if (!type.compare("modified_Batschinski_Hildebrand")){
// Get a reference to the entry in the fluid instance to simplify the code that follows
CoolProp::ViscosityModifiedBatschinskiHildebrandData &BH = fluid.transport.viscosity_higher_order.modified_Batschinski_Hildebrand;
// Set the flag for the type of this model
fluid.transport.viscosity_higher_order.type = CoolProp::ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_BATSCHINKI_HILDEBRAND;
@@ -467,7 +470,7 @@ protected:
else if (!type.compare("friction_theory")){
// Get a reference to the entry in the fluid instance to simplify the code that follows
CoolProp::ViscosityFrictionTheoryData &F = fluid.transport.viscosity_higher_order.friction_theory;
// Set the flag for the type of this model
fluid.transport.viscosity_higher_order.type = CoolProp::ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_FRICTION_THEORY;
@@ -476,7 +479,7 @@ protected:
F.Aa = cpjson::get_long_double_array(higher["Aa"]);
F.Aaa = cpjson::get_long_double_array(higher["Aaa"]);
F.Ar = cpjson::get_long_double_array(higher["Ar"]);
F.Na = cpjson::get_integer(higher,"Na");
F.Naa = cpjson::get_integer(higher,"Naa");
@@ -487,7 +490,7 @@ protected:
assert(F.Aa.size() == 3);
assert(F.Aaa.size() == 3);
assert(F.Ar.size() == 3);
F.T_reduce = cpjson::get_double(higher,"T_reduce");
if (higher.HasMember("Arr") && !higher.HasMember("Adrdr")){
@@ -581,7 +584,7 @@ protected:
}
}
// Load dilute viscosity term
if (viscosity.HasMember("dilute")){
@@ -596,7 +599,7 @@ protected:
parse_higher_order_viscosity(viscosity["higher_order"], fluid);
}
};
/// Parse the transport properties
void parse_dilute_conductivity(rapidjson::Value &dilute, CoolPropFluid & fluid)
{
@@ -639,7 +642,7 @@ protected:
// Load up the values
data.A = cpjson::get_long_double_array(dilute["A"]);
data.t = cpjson::get_long_double_array(dilute["t"]);
data.t = cpjson::get_long_double_array(dilute["t"]);
}
else{
throw ValueError(format("type [%s] is not understood for fluid %s",type.c_str(),fluid.name.c_str()));
@@ -732,7 +735,7 @@ protected:
throw ValueError(format("type [%s] is not understood for fluid %s",type.c_str(),fluid.name.c_str()));
}
};
/// Parse the thermal conductivity data
void parse_thermal_conductivity(rapidjson::Value &conductivity, CoolPropFluid & fluid)
{
@@ -784,7 +787,7 @@ protected:
parse_viscosity(transport["viscosity"],fluid);
}
// Parse thermal conductivity
// Parse thermal conductivity
if (transport.HasMember("conductivity")){
parse_thermal_conductivity(transport["conductivity"],fluid);
}
@@ -794,7 +797,7 @@ protected:
{
// Use the method of Chung to approximate the values for epsilon_over_k and sigma_eta
// Chung, T.-H.; Ajlan, M.; Lee, L. L.; Starling, K. E. Generalized Multiparameter Correlation for Nonpolar and Polar Fluid Transport Properties. Ind. Eng. Chem. Res. 1988, 27, 671-679.
// rhoc needs to be in mol/L to yield a sigma in nm,
// rhoc needs to be in mol/L to yield a sigma in nm,
long double rho_crit_molar = fluid.pEOS->reduce.rhomolar/1000.0;// [mol/m3 to mol/L]
long double Tc = fluid.pEOS->reduce.T;
fluid.transport.sigma_eta = 0.809/pow(rho_crit_molar, static_cast<long double>(1.0/3.0))/1e9; // 1e9 is to convert from nm to m
@@ -833,7 +836,7 @@ protected:
assert(fluid.name.length() > 0);
}
public:
// Default constructor;
JSONFluidLibrary(){
_is_empty = true;
@@ -844,14 +847,14 @@ public:
void add_many(rapidjson::Value &listing)
{
for (rapidjson::Value::ValueIterator itr = listing.Begin(); itr != listing.End(); ++itr)
{
{
add_one(*itr);
}
};
void add_one(rapidjson::Value &fluid_json)
{
_is_empty = false;
// Get the next index for this fluid
std::size_t index = fluid_map.size();
@@ -863,73 +866,84 @@ public:
// Fluid name
fluid.name = fluid_json["NAME"].GetString(); name_vector.push_back(fluid.name);
// CAS number
fluid.CAS = fluid_json["CAS"].GetString();
// REFPROP alias
fluid.REFPROPname = fluid_json["REFPROP_NAME"].GetString();
// Critical state
parse_crit_state(fluid_json["CRITICAL"], fluid);
if (get_debug_level() > 5){
std::cout << format("Loading fluid %s with CAS %s; %d fluids loaded\n", fluid.name.c_str(), fluid.CAS.c_str(), index);
}
try{
// CAS number
if (!fluid_json.HasMember("CAS")){ throw ValueError(format("fluid [%s] does not have \"CAS\" member",fluid.name.c_str())); }
fluid.CAS = fluid_json["CAS"].GetString();
// REFPROP alias
if (!fluid_json.HasMember("REFPROP_NAME")){ throw ValueError(format("fluid [%s] does not have \"REFPROP_NAME\" member",fluid.name.c_str())); }
fluid.REFPROPname = fluid_json["REFPROP_NAME"].GetString();
// Critical state
if (!fluid_json.HasMember("STATES")){ throw ValueError(format("fluid [%s] does not have \"STATES\" member",fluid.name.c_str())); }
if (!fluid_json["STATES"].HasMember("critical")){ throw ValueError(format("fluid[\"STATES\"] [%s] does not have \"critical\" member",fluid.name.c_str())); }
parse_crit_state(fluid_json["STATES"]["critical"], fluid);
// Aliases
fluid.aliases = cpjson::get_string_array(fluid_json["ALIASES"]);
// EOS
parse_EOS_listing(fluid_json["EOS"], fluid);
// Validate the fluid
validate(fluid);
// Ancillaries for saturation
if (!fluid_json.HasMember("ANCILLARIES")){throw ValueError(format("Ancillary curves are missing for fluid [%s]",fluid.name.c_str()));};
parse_ancillaries(fluid_json["ANCILLARIES"],fluid);
// Surface tension
if (!(fluid_json["ANCILLARIES"].HasMember("surface_tension"))){
if (get_debug_level() > 0){
std::cout << format("Surface tension curves are missing for fluid [%s]\n", fluid.name.c_str()) ;
if (get_debug_level() > 5){
std::cout << format("Loading fluid %s with CAS %s; %d fluids loaded\n", fluid.name.c_str(), fluid.CAS.c_str(), index);
}
}
else{
parse_surface_tension(fluid_json["ANCILLARIES"]["surface_tension"], fluid);
}
// Parse the environmental parameters
if (!(fluid_json.HasMember("ENVIRONMENTAL"))){
if (get_debug_level() > 0){
std::cout << format("Environmental data are missing for fluid [%s]\n", fluid.name.c_str()) ;
// Aliases
fluid.aliases = cpjson::get_string_array(fluid_json["ALIASES"]);
// EOS
parse_EOS_listing(fluid_json["EOS"], fluid);
// Validate the fluid
validate(fluid);
// Ancillaries for saturation
if (!fluid_json.HasMember("ANCILLARIES")){throw ValueError(format("Ancillary curves are missing for fluid [%s]",fluid.name.c_str()));};
parse_ancillaries(fluid_json["ANCILLARIES"],fluid);
// Surface tension
if (!(fluid_json["ANCILLARIES"].HasMember("surface_tension"))){
if (get_debug_level() > 0){
std::cout << format("Surface tension curves are missing for fluid [%s]\n", fluid.name.c_str()) ;
}
}
else{
parse_surface_tension(fluid_json["ANCILLARIES"]["surface_tension"], fluid);
}
}
else{
parse_environmental(fluid_json["ENVIRONMENTAL"], fluid);
}
// Parse the environmental parameters
if (!(fluid_json.HasMember("TRANSPORT"))){
default_transport(fluid);
}
else{
parse_transport(fluid_json["TRANSPORT"], fluid);
}
// If the fluid is ok...
// Parse the environmental parameters
if (!(fluid_json.HasMember("ENVIRONMENTAL"))){
if (get_debug_level() > 0){
std::cout << format("Environmental data are missing for fluid [%s]\n", fluid.name.c_str()) ;
}
}
else{
parse_environmental(fluid_json["ENVIRONMENTAL"], fluid);
}
// Add CAS->index mapping
string_to_index_map[fluid.CAS] = index;
// Parse the environmental parameters
if (!(fluid_json.HasMember("TRANSPORT"))){
default_transport(fluid);
}
else{
parse_transport(fluid_json["TRANSPORT"], fluid);
}
// If the fluid is ok...
// Add CAS->index mapping
string_to_index_map[fluid.CAS] = index;
// Add name->index mapping
string_to_index_map[fluid.name] = index;
// Add the aliases
for (std::size_t i = 0; i < fluid.aliases.size(); ++i)
{
string_to_index_map[fluid.aliases[i]] = index;
}
if (get_debug_level() > 5){ std::cout << format("Loaded.\n"); }
// Add name->index mapping
string_to_index_map[fluid.name] = index;
// Add the aliases
for (std::size_t i = 0; i < fluid.aliases.size(); ++i)
{
string_to_index_map[fluid.aliases[i]] = index;
}
if (get_debug_level() > 5){ std::cout << format("Loaded.\n"); }
catch (const std::exception &e){
throw ValueError(format("Unable to load fluid [%s] due to error: %s",fluid.name.c_str(),e.what()));
}
};
/// Get a CoolPropFluid instance stored in this library
/**
@@ -982,4 +996,4 @@ std::string get_fluid_list(void);
CoolPropFluid& get_fluid(std::string fluid_string);
} /* namespace CoolProp */
#endif
#endif

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

@@ -1174,8 +1174,8 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do
// Supercritical temperature
if (_T > _crit.T)
{
long double rhomelt = components[0]->pEOS->rhoLtriple;
long double rhoc = components[0]->pEOS->reduce.rhomolar;
long double rhomelt = components[0]->triple_liquid.rhomolar;
long double rhoc = components[0]->crit.rhomolar;
long double rhomin = 1e-10;
switch(other)
@@ -1228,7 +1228,7 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do
else if (_phase == iphase_liquid)
{
long double ymelt, yL, y;
long double rhomelt = components[0]->pEOS->rhoLtriple;
long double rhomelt = components[0]->triple_liquid.rhomolar;
long double rhoL = static_cast<double>(_rhoLanc);
switch(other)

70
src/Backends/Helmholtz/TransportRoutines.cpp
View File

@@ -11,13 +11,13 @@ long double TransportRoutines::viscosity_dilute_kinetic_theory(HelmholtzEOSMixtu
long double Tstar = HEOS.T()/HEOS.components[0]->transport.epsilon_over_k;
long double sigma_nm = HEOS.components[0]->transport.sigma_eta*1e9; // 1e9 to convert from m to nm
long double molar_mass_kgkmol = HEOS.molar_mass()*1000; // 1000 to convert from kg/mol to kg/kmol
// The nondimensional empirical collision integral from Neufeld
// Neufeld, P. D.; Janzen, A. R.; Aziz, R. A. Empirical Equations to Calculate 16 of the Transport Collision Integrals (l,s)*
// The nondimensional empirical collision integral from Neufeld
// Neufeld, P. D.; Janzen, A. R.; Aziz, R. A. Empirical Equations to Calculate 16 of the Transport Collision Integrals (l,s)*
// for the Lennard-Jones (12-6) Potential. J. Chem. Phys. 1972, 57, 1100-1102
long double OMEGA22 = 1.16145*pow(Tstar, static_cast<long double>(-0.14874))+0.52487*exp(-0.77320*Tstar)+2.16178*exp(-2.43787*Tstar);
// The dilute gas component -
// The dilute gas component -
return 26.692e-9*sqrt(molar_mass_kgkmol*HEOS.T())/(pow(sigma_nm, 2)*OMEGA22); // Pa-s
}
else{
@@ -40,8 +40,8 @@ long double TransportRoutines::viscosity_dilute_collision_integral(HelmholtzEOSM
const long double sigma_nm = HEOS.components[0]->transport.sigma_eta*1e9; // 1e9 to convert from m to nm
const long double molar_mass_kgkmol = molar_mass*1000; // 1000 to convert from kg/mol to kg/kmol
/// Both the collision integral \f$\mathfrak{S}^*\f$ and effective cross section \f$\Omega^{(2,2)}\f$ have the same form,
/// in general we don't care which is used. The are related through \f$\Omega^{(2,2)} = (5/4)\mathfrak{S}^*\f$
/// Both the collision integral \f$\mathfrak{S}^*\f$ and effective cross section \f$\Omega^{(2,2)}\f$ have the same form,
/// in general we don't care which is used. The are related through \f$\Omega^{(2,2)} = (5/4)\mathfrak{S}^*\f$
/// see Vesovic(JPCRD, 1990) for CO\f$_2\f$ for further information
long double summer = 0, lnTstar = log(Tstar);
for (std::size_t i = 0; i < a.size(); ++i)
@@ -49,7 +49,7 @@ long double TransportRoutines::viscosity_dilute_collision_integral(HelmholtzEOSM
summer += a[i]*pow(lnTstar,t[i]);
}
S = exp(summer);
// The dilute gas component
return C*sqrt(molar_mass_kgkmol*HEOS.T())/(pow(sigma_nm, 2)*S); // Pa-s
}
@@ -148,7 +148,7 @@ long double TransportRoutines::viscosity_initial_density_dependence_Rainwater_Fr
long double B_eta, B_eta_star;
long double Tstar = HEOS.T()/HEOS.components[0]->transport.epsilon_over_k; // [no units]
long double sigma = HEOS.components[0]->transport.sigma_eta; // [m]
long double summer = 0;
for (unsigned int i = 0; i < b.size(); ++i){
summer += b[i]*pow(Tstar, t[i]);
@@ -225,7 +225,7 @@ long double TransportRoutines::viscosity_water_hardcoded(HelmholtzEOSMixtureBack
double x_mu=0.068,qc=1/1.9,qd=1/1.1,nu=0.630,gamma=1.239,zeta_0=0.13,LAMBDA_0=0.06,Tbar_R=1.5, pstar, Tstar, rhostar;
double delta,tau,mubar_0,mubar_1,mubar_2,drhodp,drhodp_R,DeltaChibar,zeta,w,L,Y,psi_D,Tbar,rhobar;
double drhobar_dpbar,drhobar_dpbar_R,R_Water;
pstar = 22.064e6; // [Pa]
Tstar = 647.096; // [K]
rhostar = 322; // [kg/m^3]
@@ -246,7 +246,7 @@ long double TransportRoutines::viscosity_water_hardcoded(HelmholtzEOSMixtureBack
tau=1/Tbar_R;
drhodp_R=1/(R_Water*Tbar_R*Tstar*(1+2*rhobar*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,tau,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,tau, delta)));
drhobar_dpbar_R = pstar/rhostar*drhodp_R;
DeltaChibar=rhobar*(drhobar_dpbar-drhobar_dpbar_R*Tbar_R/Tbar);
if (DeltaChibar<0)
DeltaChibar=0;
@@ -330,7 +330,7 @@ long double TransportRoutines::viscosity_higher_order_friction_theory(HelmholtzE
double deltapr = pr - pid;
double eta_f = ka*pa + kr*deltapr + ki*pid + kaa*pa*pa + kdrdr*deltapr*deltapr + krr*pr*pr + kii*pid*pid + krrr*pr*pr*pr + kaaa*pa*pa*pa;
return eta_f; //[Pa-s]
}
else{
@@ -345,10 +345,10 @@ long double TransportRoutines::viscosity_helium_hardcoded(HelmholtzEOSMixtureBac
double eta_0,eta_0_slash, eta_E_slash, B,C,D,ln_eta,x;
//
// Arp, V.D., McCarty, R.D., and Friend, D.G.,
// "Thermophysical Properties of Helium-4 from 0.8 to 1500 K with Pressures to 2000 MPa",
// "Thermophysical Properties of Helium-4 from 0.8 to 1500 K with Pressures to 2000 MPa",
// NIST Technical Note 1334 (revised), 1998.
//
// Using Arp NIST report
//
// Using Arp NIST report
// Report is not clear on viscosity, referring to REFPROP source code for clarity
// Correlation wants density in g/cm^3; kg/m^3 --> g/cm^3, divide by 1000
@@ -391,9 +391,6 @@ long double TransportRoutines::viscosity_R23_hardcoded(HelmholtzEOSMixtureBacken
rhocbar = 7.5114,
Tc = 299.2793,
DELTAeta_max = 3.967,
k = 1.380658e-23,
N_A = 6.022137e23, // 1/mol
pi = 3.141592654, //
Ru = 8.31451,
molar_mass = 70.014;
@@ -418,7 +415,7 @@ long double TransportRoutines::viscosity_R23_hardcoded(HelmholtzEOSMixtureBacken
long double TransportRoutines::viscosity_dilute_ethane(HelmholtzEOSMixtureBackend &HEOS)
{
double C[] = {0, -3.0328138281, 16.918880086, -37.189364917, 41.288861858, -24.615921140, 8.9488430959, -1.8739245042, 0.20966101390, -9.6570437074e-3};
double OMEGA_2_2 = 0, e_k = 245, sigma = 0.43682, Tstar;
double OMEGA_2_2 = 0, e_k = 245, Tstar;
Tstar = HEOS.T()/e_k;
for (int i = 1; i<= 9; i++)
@@ -572,7 +569,7 @@ long double TransportRoutines::conductivity_critical_hardcoded_R123(HelmholtzEOS
long double TransportRoutines::conductivity_critical_hardcoded_CO2_ScalabrinJPCRD2006(HelmholtzEOSMixtureBackend &HEOS){
long double nc = 0.775547504e-3*4.81384, Tr = HEOS.T()/304.1282, alpha, rhor = HEOS.keyed_output(iDmass)/467.6;
static long double a[] = {0.0, 3.0, 6.70697, 0.94604, 0.30, 0.30, 0.39751, 0.33791, 0.77963, 0.79857, 0.90, 0.02, 0.20};
// Equation 6 from Scalabrin
alpha = 1-a[10]*acosh(1+a[11]*pow(pow(1-Tr,2),a[12]));
@@ -598,8 +595,8 @@ long double TransportRoutines::conductivity_dilute_hardcoded_CO2(HelmholtzEOSMix
//Vesovic Eq. 12 [no units]
double r = sqrt(2.0/5.0*cint_k);
// According to REFPROP, 1+r^2 = cp-2.5R. This is unclear to me but seems to suggest that cint/k is the difference
// between the ideal gas specific heat and a monatomic specific heat of 5/2*R. Using the form of cint/k from Vesovic
// According to REFPROP, 1+r^2 = cp-2.5R. This is unclear to me but seems to suggest that cint/k is the difference
// between the ideal gas specific heat and a monatomic specific heat of 5/2*R. Using the form of cint/k from Vesovic
// does not yield exactly the correct values
Tstar = HEOS.T()/e_k;
@@ -670,15 +667,14 @@ long double TransportRoutines::conductivity_hardcoded_water(HelmholtzEOSMixtureB
// Finite density contribution
lambdabar_1=exp(rhobar*sum);
double nu=0.630,GAMMA =177.8514,gamma=1.239,xi_0=0.13,Lambda_0=0.06,Tr_bar=1.5,tau_ref = Tr_bar*Tstar/HEOS.T(),
double nu=0.630,GAMMA =177.8514,gamma=1.239,xi_0=0.13,Lambda_0=0.06,Tr_bar=1.5,
qd_bar=1/0.4,pi=3.141592654, delta = HEOS.delta(), R=461.51805;//J/kg/K
tau=1/Tbar;
double drhodp = 1/(R*HEOS.T()*(1+2*rhobar*HEOS.dalphar_dDelta()+rhobar*rhobar*HEOS.d2alphar_dDelta2()));
double drhobar_dpbar = pstar/rhostar*drhodp;
double drhodp_Trbar = 1/(R*Tr_bar*Tstar*(1+2*rhobar*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,1/Tr_bar,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,1/Tr_bar,delta)));
double drhobar_dpbar_Trbar = pstar/rhostar*drhodp_Trbar;
double cpmol = HEOS.cpmolar(); // [J/mol/K]
double cp = HEOS.cpmass(); // [J/kg/K]
double cv = HEOS.cvmass(); // [J/kg/K]
double cpbar = cp/R; //[-]
@@ -731,17 +727,17 @@ long double TransportRoutines::conductivity_hardcoded_R23(HelmholtzEOSMixtureBac
long double TransportRoutines::conductivity_critical_hardcoded_ammonia(HelmholtzEOSMixtureBackend &HEOS){
/*
From "Thermal Conductivity of Ammonia in a Large
/*
From "Thermal Conductivity of Ammonia in a Large
Temperature and Pressure Range Including the Critical Region"
by R. Tufeu, D.Y. Ivanov, Y. Garrabos, B. Le Neindre,
by R. Tufeu, D.Y. Ivanov, Y. Garrabos, B. Le Neindre,
Bereicht der Bunsengesellschaft Phys. Chem. 88 (1984) 422-427
*/
double T = HEOS.T(), Tc = 405.4, rhoc = 235, rho;
double LAMBDA=1.2, nu=0.63, gamma =1.24, DELTA=0.50,t,zeta_0_plus=1.34e-10,a_zeta=1,GAMMA_0_plus=0.423e-8;
double pi=3.141592654,a_chi,k_B=1.3806504e-23,X_T,DELTA_lambda,dPdT,eta_B,DELTA_lambda_id,DELTA_lambda_i;
rho = HEOS.keyed_output(CoolProp::iDmass);
t = fabs((T-Tc)/Tc);
a_chi = a_zeta/0.7;
@@ -785,14 +781,14 @@ long double TransportRoutines::conductivity_hardcoded_helium(HelmholtzEOSMixture
lambda_e = (c[0]+c[1]*T+c[2]*pow(T,1/3.0)+c[3]*pow(T,2.0/3.0))*rho
+(c[4]+c[5]*pow(T,1.0/3.0)+c[6]*pow(T,2.0/3.0))*rho*rho*rho
+(c[7]+c[8]*pow(T,1.0/3.0)+c[9]*pow(T,2.0/3.0)+c[10]/T)*rho*rho*log(rho/rhoc);
// Critical component
lambda_c = 0.0;
if (3.5 < T && T < 12)
{
double x0 = 0.392, E1 = 2.8461, E2 = 0.27156, beta = 0.3554, gamma = 1.1743, delta = 4.304, rhoc_crit = 69.158,
Tc = 5.18992, pc = 2.2746e5, R = 4.633e-10, m = 6.6455255e-27, k = 1.38066e-23, pi = M_PI;
double x0 = 0.392, E1 = 2.8461, E2 = 0.27156, beta = 0.3554, gamma = 1.1743, delta = 4.304, rhoc_crit = 69.158,
Tc = 5.18992, pc = 2.2746e5;
double DeltaT = fabs(1-T/Tc), DeltaRho = fabs(1-rho/rhoc_crit);
double eta = HEOS.viscosity(); // [Pa-s]
@@ -810,7 +806,7 @@ long double TransportRoutines::conductivity_hardcoded_helium(HelmholtzEOSMixture
double x = pow(DeltaT/DeltaRho,1/beta);
double h = E1*(1 + x/x0)*pow(1 + E2*pow(1 + x/x0, 2/beta), (gamma-1)/(2*beta));
/**
/**
dh/dx derived using sympy:
E1,x,x0,E2,beta,gamma = symbols('E1,x,x0,E2,beta,gamma')
@@ -824,7 +820,7 @@ long double TransportRoutines::conductivity_hardcoded_helium(HelmholtzEOSMixture
K_Tbar = W*K_T + (1-W)*K_Tprime;
}
// 3.4685233d-17 and 3.726229668d0 are "magical" coefficients that are present in the REFPROP source to yield the right values. Not clear why these values are needed.
// 3.4685233d-17 and 3.726229668d0 are "magical" coefficients that are present in the REFPROP source to yield the right values. Not clear why these values are needed.
// Also, the form of the critical term in REFPROP does not agree with Hands paper. EL and MH from NIST are not sure where these coefficients come from.
lambda_c = 3.4685233e-17*3.726229668*sqrt(K_Tbar)*pow(T,2)/rho/eta*pow(dpdT,2)*exp(-18.66*pow(DeltaT,2)-4.25*pow(DeltaRho,4));
}
@@ -843,7 +839,7 @@ long double TransportRoutines::viscosity_ECS(HelmholtzEOSMixtureBackend &HEOS, H
// Get a reference to the ECS data
CoolProp::ViscosityECSVariables &ECS = HEOS.components[0]->transport.viscosity_ecs;
// The correction polynomial psi_eta
double psi = 0;
for (std::size_t i=0; i < ECS.psi_a.size(); i++){
@@ -897,7 +893,7 @@ long double TransportRoutines::conductivity_ECS(HelmholtzEOSMixtureBackend &HEOS
// Get a reference to the ECS data
CoolProp::ConductivityECSVariables &ECS = HEOS.components[0]->transport.conductivity_ecs;
// The correction polynomial psi_eta in rho/rho_red
double psi = 0;
for (std::size_t i=0; i < ECS.psi_a.size(); ++i){
@@ -954,4 +950,4 @@ long double TransportRoutines::conductivity_ECS(HelmholtzEOSMixtureBackend &HEOS
}
}; /* namespace CoolProp */
}; /* namespace CoolProp */