github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-02-09 21:35:28 -05:00
Code Issues Packages Projects Releases Wiki Activity

Added crossplatform_shared_ptr to allow shared_ptr cleanly on all architectures without needing variable namespace definition

Browse Source
This commit is contained in:
Ian Bell
2014-06-03 20:58:24 +02:00
parent a7363b067c
commit cb497b443c
12 changed files with 280 additions and 254 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/Backends/Helmholtz/ExcessHEFunction.h
View File

@@ -4,6 +4,7 @@
#include <memory>
#include <vector>
#include "CoolPropFluid.h"
#include "crossplatform_shared_ptr.h"
namespace CoolProp{

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

@@ -164,7 +164,7 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity_background(long double et
// Residual part
long double B_eta_initial = TransportRoutines::viscosity_initial_density_dependence_Rainwater_Friend(*this);
long double rho = rhomolar();
long double initial_part = eta_dilute*B_eta_initial*rhomolar();
long double initial_part = eta_dilute*B_eta_initial*rho;
// Higher order terms
long double delta_eta_h;

104
src/Backends/REFPROP/REFPROPMixtureBackend.cpp
View File

@@ -51,7 +51,7 @@ surface tension N/m
#include <stdio.h>
#include <iostream>
#include <assert.h>
#include "crossplatform_shared_ptr.h"
#if defined(_MSC_VER)
@@ -63,24 +63,24 @@ surface tension N/m
#include <sys/stat.h>
#endif
// Some constants for REFPROP... defined by macros for ease of use
// Some constants for REFPROP... defined by macros for ease of use
#define refpropcharlength 255
#define filepathlength 255
#define lengthofreference 3
#define errormessagelength 255
#define ncmax 20 // Note: ncmax is the max number of components
#define numparams 72
#define numparams 72
#define maxcoefs 50
std::string LoadedREFPROPRef;
// Some constants for REFPROP... defined by macros for ease of use
// Some constants for REFPROP... defined by macros for ease of use
#define refpropcharlength 255
#define filepathlength 255
#define lengthofreference 3
#define errormessagelength 255
#define ncmax 20 // Note: ncmax is the max number of components
#define numparams 72
#define numparams 72
#define maxcoefs 50
// Check windows
@@ -371,7 +371,7 @@ bool load_REFPROP()
// INCREASE ROBUSTNESS. ALWAYS THROW AN ERROR ON THE ELSE.
#error "Must define either ENV32BIT or ENV64BIT"
#endif
#elif defined(__ISLINUX__)
RefpropdllInstance = dlopen ("librefprop.so", RTLD_LAZY);
#elif defined(__ISAPPLE__)
@@ -401,8 +401,8 @@ bool load_REFPROP()
}
#if defined(__ISWINDOWS__)
// Get data associated with path using the windows libraries,
// Get data associated with path using the windows libraries,
// and if you can (result == 0), the path exists
#ifdef __MINGW32__
struct stat buf;
@@ -434,7 +434,7 @@ REFPROPMixtureBackend::REFPROPMixtureBackend(const std::vector<std::string>& flu
// Do the REFPROP instantiation for this fluid
_mole_fractions_set = false;
// Try to add this fluid to REFPROP - might want to think about making array of
// Try to add this fluid to REFPROP - might want to think about making array of
// components and setting mole fractions if they change a lot.
this->set_REFPROP_fluids(fluid_names);
@@ -521,8 +521,8 @@ void REFPROPMixtureBackend::set_REFPROP_fluids(const std::vector<std::string> &f
// Load REFPROP if it isn't loaded yet
load_REFPROP();
// If the name of the refrigerant doesn't match
// If the name of the refrigerant doesn't match
// that of the currently loaded refrigerant
if (LoadedREFPROPRef.compare(components_joined))
{
@@ -593,7 +593,7 @@ double REFPROPMixtureBackend::calc_melt_Tmax()
MELTPdll(&pmax_kPa, &(mole_fractions[0]),
&Tmax_melt,
&ierr,herr,errormessagelength); // Error message
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
//else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return Tmax_melt;
}
@@ -611,7 +611,7 @@ long double REFPROPMixtureBackend::calc_melt_p_T(long double T)
MELTTdll(&_T, &(mole_fractions[0]),
&p_kPa,
&ierr,herr,errormessagelength); // Error message
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
//else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return p_kPa*1000;
}
@@ -632,7 +632,7 @@ long double REFPROPMixtureBackend::calc_viscosity(void)
TRNPRPdll(&_T,&rhomol_L,&(mole_fractions[0]), // Inputs
&eta,&tcx, // Outputs
&ierr,herr,errormessagelength); // Error message
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
//else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
_viscosity = 1e-6*eta;
_conductivity = tcx;
@@ -652,7 +652,7 @@ long double REFPROPMixtureBackend::calc_surface_tension(void)
SURFTdll(&_T, &rho_mol_L, &(mole_fractions[0]), // Inputs
&sigma, // Outputs
&ierr, herr, errormessagelength); // Error message
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
//else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
_surface_tension = sigma;
return static_cast<double>(_surface_tension);
@@ -667,11 +667,11 @@ long double REFPROPMixtureBackend::calc_fugacity_coefficient(int i)
FUGCOFdll(&_T, &rho_mol_L, &(mole_fractions[0]), // Inputs
&(fug_cof[0]), // Outputs
&ierr, herr, errormessagelength); // Error message
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }
//else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return static_cast<long double>(fug_cof[i]);
}
void REFPROPMixtureBackend::update(long input_pair, double value1, double value2)
{
double rho_mol_L=_HUGE, rhoLmol_L=_HUGE, rhoVmol_L=_HUGE,
@@ -681,7 +681,7 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
char herr[255];
clear();
// Check that mole fractions have been set, etc.
check_status();
@@ -689,8 +689,8 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
WMOLdll(&(mole_fractions[0]), &mm); // returns mole mass in kg/kmol
_molar_mass = 0.001*mm; // [kg/mol]
switch(input_pair)
{
case PT_INPUTS:
@@ -724,9 +724,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
TDFLSHdll(&_T,&rho_mol_L,&(mole_fractions[0]),&p_kPa,
&rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms
&q,&emol,&hmol,&smol,&cvmol,&cpmol,&w,
&ierr,herr,errormessagelength);
&ierr,herr,errormessagelength);
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = p_kPa*1000;
if (0)
@@ -757,7 +757,7 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_rhomolar = value1;
_rhomolar = value1;
_p = value2;
if (0)
{
@@ -783,9 +783,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
DHFLSHdll(&rho_mol_L,&hmol,&(mole_fractions[0]),&_T,&p_kPa,
&rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms
&q,&emol,&smol,&cvmol,&cpmol,&w,
&ierr,herr,errormessagelength);
&ierr,herr,errormessagelength);
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = p_kPa*1000;
if (0)
@@ -813,9 +813,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
DSFLSHdll(&rho_mol_L,&smol,&(mole_fractions[0]),&_T,&p_kPa,
&rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms
&q,&emol,&hmol,&cvmol,&cpmol,&w,
&ierr,herr,errormessagelength);
&ierr,herr,errormessagelength);
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = p_kPa*1000;
if (0)
@@ -843,9 +843,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
DEFLSHdll(&rho_mol_L,&emol,&(mole_fractions[0]),&_T,&p_kPa,
&rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms
&q,&hmol,&hmol,&cvmol,&cpmol,&w,
&ierr,herr,errormessagelength);
&ierr,herr,errormessagelength);
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); }// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
// Set all cache values that can be set with unit conversion to SI
_p = p_kPa*1000;
if (0)
@@ -889,7 +889,7 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
{
// Call again, but this time with molar units
// H: [J/kg] * [kg/mol] -> [J/mol]
update(HmolarP_INPUTS, value1 * (double)_molar_mass, value2);
update(HmolarP_INPUTS, value1 * (double)_molar_mass, value2);
return;
}
case PSmolar_INPUTS:
@@ -919,7 +919,7 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
{
// Call again, but this time with molar units
// S: [J/kg/K] * [kg/mol] -> [J/mol/K]
update(PSmolar_INPUTS, value1, value2*(double)_molar_mass);
update(PSmolar_INPUTS, value1, value2*(double)_molar_mass);
return;
}
case PUmolar_INPUTS:
@@ -948,9 +948,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
}
case PUmass_INPUTS:
{
// Call again, but this time with molar units
// Call again, but this time with molar units
// U: [J/kg] * [kg/mol] -> [J/mol]
update(PUmolar_INPUTS, value1, value2*(double)_molar_mass);
update(PUmolar_INPUTS, value1, value2*(double)_molar_mass);
return;
}
case HmolarSmolar_INPUTS:
@@ -977,16 +977,16 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
}
case HmassSmass_INPUTS:
{
// Call again, but this time with molar units
// Call again, but this time with molar units
// H: [J/kg] * [kg/mol] -> [J/mol/K]
// S: [J/kg/K] * [kg/mol] -> [J/mol/K]
update(HmolarSmolar_INPUTS, value1 * (double)_molar_mass, value2 * (double)_molar_mass);
update(HmolarSmolar_INPUTS, value1 * (double)_molar_mass, value2 * (double)_molar_mass);
return;
}
case SmolarUmolar_INPUTS:
{
// Unit conversion for REFPROP
smol = value1; emol = value2;
smol = value1; emol = value2;
// from REFPROP: subroutine ESFLSH (e,s,z,t,p,D,Dl,Dv,x,y,q,h,cv,cp,w,ierr,herr)
ESFLSHdll(&emol,&smol,&(mole_fractions[0]),&_T,&p_kPa,&rho_mol_L,
@@ -1008,10 +1008,10 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
}
case SmassUmass_INPUTS:
{
// Call again, but this time with molar units
// S: [J/kg/K] * [kg/mol] -> [J/mol/K],
// Call again, but this time with molar units
// S: [J/kg/K] * [kg/mol] -> [J/mol/K],
// U: [J/kg] * [kg/mol] -> [J/mol]
update(SmolarUmolar_INPUTS, value1 * (double)_molar_mass, value2 * (double)_molar_mass);
update(SmolarUmolar_INPUTS, value1 * (double)_molar_mass, value2 * (double)_molar_mass);
return;
}
case SmolarT_INPUTS:
@@ -1047,9 +1047,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
}
case SmassT_INPUTS:
{
// Call again, but this time with molar units
// Call again, but this time with molar units
// S: [J/kg/K] * [kg/mol] -> [J/mol/K]
update(SmolarT_INPUTS, value1 * (double)_molar_mass, value2 );
update(SmolarT_INPUTS, value1 * (double)_molar_mass, value2 );
return;
}
case HmolarT_INPUTS:
@@ -1085,9 +1085,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
}
case HmassT_INPUTS:
{
// Call again, but this time with molar units
// Call again, but this time with molar units
// H: [J/kg] * [kg/mol] -> [J/mol]
update(HmolarT_INPUTS, value1 * (double)_molar_mass, value2 );
update(HmolarT_INPUTS, value1 * (double)_molar_mass, value2 );
return;
}
case TUmolar_INPUTS:
@@ -1123,9 +1123,9 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
}
case TUmass_INPUTS:
{
// Call again, but this time with molar units
// Call again, but this time with molar units
// U: [J/kg] * [kg/mol] -> [J/mol]
update(TUmolar_INPUTS, value1, value2 * (double)_molar_mass);
update(TUmolar_INPUTS, value1, value2 * (double)_molar_mass);
return;
}
case PQ_INPUTS:
@@ -1194,7 +1194,7 @@ void REFPROPMixtureBackend::update(long input_pair, double value1, double value2
{
throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
}
};
// Set these common variables that are used in every flash calculation
_hmolar = hmol;
@@ -1222,7 +1222,7 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]")
{
std::string Name = (*it);
std::string RPName = CoolProp::get_fluid_param_string((*it),"REFPROP_name");
// Skip fluids not in REFPROP
if (RPName.find("N/A") == 0){continue;}
@@ -1239,7 +1239,7 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]")
CAPTURE(RPName);
CAPTURE(rho_CP);
CAPTURE(rho_RP);
double DH = (rho_RP-rho_CP)/rho_RP;
CHECK(fabs(DH) < 0.005);
}
@@ -1252,7 +1252,7 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]")
{
std::string Name = (*it);
std::string RPName = CoolProp::get_fluid_param_string((*it),"REFPROP_name");
// Skip fluids not in REFPROP
if (RPName.find("N/A") == 0){continue;}
@@ -1270,7 +1270,7 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]")
CAPTURE(cp_CP);
CAPTURE(cp_RP);
CAPTURE(0.9*Tr);
double Dcp = (cp_RP-cp_CP)/cp_RP;
CHECK(fabs(Dcp) < 0.005);
}
@@ -1283,7 +1283,7 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]")
{
std::string Name = (*it);
std::string RPName = CoolProp::get_fluid_param_string((*it),"REFPROP_name");
// Skip fluids not in REFPROP
if (RPName.find("N/A") == 0){continue;}
@@ -1312,4 +1312,4 @@ TEST_CASE("Check REFPROP H,S reference states equal to CoolProp","[REFPROP]")
}
}
#endif
#endif

48
src/CoolProp.cpp
View File

@@ -45,8 +45,8 @@ int get_debug_level(void){return debug_level;}
//int global_Phase = -1;
void set_warning_string(std::string warning){
warning_string = warning;
void set_warning_string(std::string warning){
warning_string = warning;
}
void set_error_string(std::string error){
error_string = error;
@@ -77,7 +77,7 @@ void set_error_string(std::string error){
// // First check whether it is one of the Brines that does
// // not have a pure-fluid equivalent in CoolProp
// if (
// strcmp(Ref,"HC-10")==0 ||
// strcmp(Ref,"HC-10")==0 ||
// strncmp(Ref,"PG-",3)==0 ||
// strncmp(Ref,"EG-",3)==0 ||
// strncmp(Ref,"EA-",3)==0 ||
@@ -237,7 +237,7 @@ std::string extract_fractions(const std::string &fluid_string, std::vector<doubl
std::string all_pairs, backend_string;
// Start at the "::" if it is found to chop off the delimiter between backend and fluid
int i = fluid_string.find("::");
// If no backend/fluid delimiter
if (i < 0){
// Just use the full string
@@ -257,9 +257,9 @@ std::string extract_fractions(const std::string &fluid_string, std::vector<doubl
for (std::size_t i = 0; i < pairs.size(); ++i)
{
std::string fluid = pairs[i];
// Must end with ']'
if (fluid[fluid.size()-1] != ']')
if (fluid[fluid.size()-1] != ']')
throw ValueError(format("Fluid entry [%s] must end with ']' character",pairs[i].c_str()));
// Split at '[', but first remove the ']' from the end by taking a substring
@@ -277,7 +277,7 @@ std::string extract_fractions(const std::string &fluid_string, std::vector<doubl
// And add to vector
fractions.push_back(f);
// Add name
names.push_back(name);
}
@@ -325,7 +325,7 @@ double _PropsSI(const std::string &Output, const std::string &Name1, double Prop
// First check if it is a trivial input (critical/max parameters for instance)
// TODO: check for trivial inputs that do not require the use of the eos
/*if (State->is_trivial_output(iOutput))
{
{
double val = State->trivial_keyed_output(iOutput);
return val;
};*/
@@ -339,7 +339,7 @@ double _PropsSI(const std::string &Output, const std::string &Name1, double Prop
// Return the value
return val;
}
catch(...){
catch(...){
throw;
}
}
@@ -355,7 +355,7 @@ double PropsSI(const char *Output, const char *Name1, double Prop1, const char *
double PropsSI(const std::string &Output, const std::string &Name1, double Prop1, const std::string &Name2, double Prop2, const std::string &Ref)
{
// In this function the error catching happens;
try{
try{
// Extract fractions if they are encoded in the fluid string
if (has_fractions_in_string(Ref))
{
@@ -392,7 +392,7 @@ double Props1SI(std::string FluidName,std::string Output)
return PropsSI(Output,"",0,"",0,FluidName);
}
//EXPORT_CODE double CONVENTION IProps(long iOutput, long iName1, double Prop1, long iName2, double Prop2, long iFluid)
//{
// Prop1 = convert_from_unit_system_to_SI(iName1, Prop1, get_standard_unit_system());
@@ -408,10 +408,10 @@ double Props1SI(std::string FluidName,std::string Output)
// if (get_debug_level()>3){
// std::cout << format("%s:%d: _CoolProp_Fluid_PropsSI(%d,%d,%g,%d,%g,%s)\n",__FILE__,__LINE__,iOutput,iName1, Prop1, iName2, Prop2, pFluid->get_name().c_str()).c_str();
// }
// if (iName1 == iT){
// T = Prop1;}
// if (iName1 == iT){
// T = Prop1;}
// else if (iName2 == iT){
// T = Prop2;}
// T = Prop2;}
//
// // Generate a State instance wrapped around the Fluid instance
// CoolPropStateClassSI CPS(pFluid);
@@ -567,7 +567,7 @@ double Props1SI(std::string FluidName,std::string Output)
//
// // Update the class
// CPS.update(iT,T,iD,rho);
//
//
// // Get the output value
// val = CPS.keyed_output(iTerm);
// break;
@@ -590,7 +590,7 @@ double Props1SI(std::string FluidName,std::string Output)
// val = CPS.keyed_output(iTerm);
// break;
// }
//
//
// default:
// throw ValueError(format("Sorry DerivTerms is a work in progress, your derivative term [%d] is not available!",iTerm));
// }
@@ -642,7 +642,8 @@ double Props1SI(std::string FluidName,std::string Output)
void set_reference_stateS(std::string Ref, std::string reference_state)
{
std::tr1::shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> HEOS;
HEOS.reset(new CoolProp::HelmholtzEOSMixtureBackend(std::vector<std::string>(1, Ref)));
std::vector<std::string> _comps(1, Ref);
HEOS.reset(new CoolProp::HelmholtzEOSMixtureBackend(_comps));
if (!reference_state.compare("IIR"))
{
@@ -687,7 +688,7 @@ void set_reference_stateS(std::string Ref, std::string reference_state)
HEOS->get_components()[0]->pEOS->alpha0.EnthalpyEntropyOffset.set(0, 0, "");
}
else
{
{
throw ValueError(format("reference state string is invalid: [%s]",reference_state.c_str()));
}
}
@@ -714,10 +715,10 @@ void set_reference_stateS(std::string Ref, std::string reference_state)
//}
std::string get_BibTeXKey(std::string Ref, std::string key)
{
{
std::vector<std::string> names(1, Ref);
HelmholtzEOSMixtureBackend HEOS(names);
if (!key.compare("EOS")){ return HEOS.get_components()[0]->pEOS->BibTeX_EOS; }
else if (!key.compare("CP0")){ return HEOS.get_components()[0]->pEOS->BibTeX_CP0; }
else if (!key.compare("VISCOSITY")){ return HEOS.get_components()[0]->transport.BibTeX_viscosity; }
@@ -756,8 +757,9 @@ std::string get_global_param_string(std::string ParamName)
std::string get_fluid_param_string(std::string FluidName, std::string ParamName)
{
try{
std::tr1::shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> HEOS(new CoolProp::HelmholtzEOSMixtureBackend(std::vector<std::string>(1,FluidName)));
std::vector<std::string> comps(1,FluidName);
std::tr1::shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> HEOS(new CoolProp::HelmholtzEOSMixtureBackend(comps));
CoolProp::CoolPropFluid *fluid = HEOS->get_components()[0];
if (!ParamName.compare("aliases"))
@@ -790,4 +792,4 @@ std::string get_fluid_param_string(std::string FluidName, std::string ParamName)
}
}
} /* namespace CoolProp */
} /* namespace CoolProp */

113
src/Helmholtz.cpp
View File

@@ -375,7 +375,7 @@ long double ResidualHelmholtzExponential::dTau3(const long double &tau, const lo
void ResidualHelmholtzGaussian::to_json(rapidjson::Value &el, rapidjson::Document &doc)
{
el.AddMember("type","ResidualHelmholtzGaussian",doc.GetAllocator());
rapidjson::Value _n(rapidjson::kArrayType), _d(rapidjson::kArrayType), _t(rapidjson::kArrayType),
rapidjson::Value _n(rapidjson::kArrayType), _d(rapidjson::kArrayType), _t(rapidjson::kArrayType),
_eta(rapidjson::kArrayType), _epsilon(rapidjson::kArrayType), _beta(rapidjson::kArrayType), _gamma(rapidjson::kArrayType);
for (unsigned int i=0; i<=N; ++i)
{
@@ -515,8 +515,8 @@ long double ResidualHelmholtzGaussian::dTau3(const long double &tau, const long
void ResidualHelmholtzGERG2008Gaussian::to_json(rapidjson::Value &el, rapidjson::Document &doc)
{
el.AddMember("type","ResidualHelmholtzGERG2008Gaussian",doc.GetAllocator());
rapidjson::Value _n(rapidjson::kArrayType), _d(rapidjson::kArrayType), _t(rapidjson::kArrayType),
_eta(rapidjson::kArrayType), _epsilon(rapidjson::kArrayType),
rapidjson::Value _n(rapidjson::kArrayType), _d(rapidjson::kArrayType), _t(rapidjson::kArrayType),
_eta(rapidjson::kArrayType), _epsilon(rapidjson::kArrayType),
_beta(rapidjson::kArrayType), _gamma(rapidjson::kArrayType);
for (unsigned int i=0; i<=N; ++i)
{
@@ -581,7 +581,7 @@ long double ResidualHelmholtzGERG2008Gaussian::dDelta3(const long double &tau, c
{
/**
Term derived in sympy using
from sympy import *
n_i, d_i, t_i, tau, delta, eta, gamma, beta, epsilon = symbols('n_i, d_i, t_i, tau, delta, eta, gamma, beta, epsilon')
psi = exp(-eta*(delta-epsilon)**2-beta*(delta-gamma))
@@ -653,8 +653,8 @@ long double ResidualHelmholtzGERG2008Gaussian::dTau3(const long double &tau, con
void ResidualHelmholtzLemmon2005::to_json(rapidjson::Value &el, rapidjson::Document &doc)
{
el.AddMember("type","ResidualHelmholtzLemmon2005",doc.GetAllocator());
rapidjson::Value _n(rapidjson::kArrayType), _d(rapidjson::kArrayType),
_t(rapidjson::kArrayType), _l(rapidjson::kArrayType),
rapidjson::Value _n(rapidjson::kArrayType), _d(rapidjson::kArrayType),
_t(rapidjson::kArrayType), _l(rapidjson::kArrayType),
_m(rapidjson::kArrayType);
for (unsigned int i=0;i<=N; ++i)
{
@@ -743,7 +743,7 @@ long double ResidualHelmholtzLemmon2005::dDelta2(const long double &tau, const l
ResidualHelmholtzLemmon2005Element &el = elements[i];
long double ni = el.n, ti = el.t, di = el.d, md = el.md;
int li = el.l, mi = el.m;
if (li != 0 && mi != 0){
if (li != 0 && mi != 0){
long double pow_delta_li = pow(delta, li);
long double pow_tau_mi = pow(tau, md);
s[i] = ni*((di-li*pow_delta_li)*(di-1.0-li*pow_delta_li) - li*li*pow_delta_li)*exp(ti*log_tau+(di-2)*log_delta-pow_delta_li-pow_tau_mi);
@@ -933,8 +933,8 @@ void ResidualHelmholtzNonAnalytic::to_json(rapidjson::Value &el, rapidjson::Docu
{
el.AddMember("type","ResidualHelmholtzNonAnalytic",doc.GetAllocator());
rapidjson::Value _n(rapidjson::kArrayType), _a(rapidjson::kArrayType), _b(rapidjson::kArrayType),
_beta(rapidjson::kArrayType), _A(rapidjson::kArrayType), _B(rapidjson::kArrayType),
rapidjson::Value _n(rapidjson::kArrayType), _a(rapidjson::kArrayType), _b(rapidjson::kArrayType),
_beta(rapidjson::kArrayType), _A(rapidjson::kArrayType), _B(rapidjson::kArrayType),
_C(rapidjson::kArrayType), _D(rapidjson::kArrayType);
for (unsigned int i=0; i<=N; ++i)
{
@@ -969,7 +969,7 @@ long double ResidualHelmholtzNonAnalytic::base(const long double &tau, const lon
long double theta=(1.0-tau)+Ai*pow(pow(delta-1.0,2),1.0/(2.0*betai));
long double DELTA=pow(theta,2)+Bi*pow(pow(delta-1.0,2),ai);
long double PSI=exp(-Ci*pow(delta-1.0,2)-Di*pow(tau-1.0,2));
s[i] = ni*pow(DELTA, bi)*delta*PSI;
}
return std::accumulate(s.begin(), s.end(), 0.0);
@@ -988,7 +988,7 @@ long double ResidualHelmholtzNonAnalytic::dDelta(const long double &tau, const l
long double PSI=exp(-Ci*pow(delta-1.0,2)-Di*pow(tau-1.0,2));
long double dPSI_dDelta=-2.0*Ci*(delta-1.0)*PSI;
long double dDELTA_dDelta=(delta-1.0)*(Ai*theta*2.0/betai*pow(pow(delta-1.0,2),1.0/(2.0*betai)-1.0)+2.0*Bi*ai*pow(pow(delta-1.0,2),ai-1.0));
// At critical point, DELTA is 0, and 1/0^n is undefined
if (fabs(DELTA) < 10*DBL_EPSILON)
{
@@ -1016,7 +1016,7 @@ long double ResidualHelmholtzNonAnalytic::dTau(const long double &tau, const lon
long double dPSI_dTau=-2.0*Di*(tau-1.0)*PSI;
long double dDELTA_dDelta=(delta-1.0)*(Ai*theta*2.0/betai*pow(pow(delta-1.0,2),1.0/(2.0*betai)-1.0)+2.0*Bi*ai*pow(pow(delta-1.0,2),ai-1.0));
long double dDELTAbi_dTau=-2.0*theta*bi*pow(DELTA,bi-1.0);
// At critical point, DELTA is 0, and 1/0^n is undefined
if (fabs(DELTA) < 10*DBL_EPSILON)
{
@@ -1048,7 +1048,7 @@ long double ResidualHelmholtzNonAnalytic::dDelta2(const long double &tau, const
long double dPSI2_dDelta2=(2.0*Ci*pow(delta-1.0,2)-1.0)*2.0*Ci*PSI;
long double dDELTA2_dDelta2=1.0/(delta-1.0)*dDELTA_dDelta+pow(delta-1.0,2)*(4.0*Bi*ai*(ai-1.0)*pow(pow(delta-1.0,2),ai-2.0)+2.0*pow(Ai/betai,2)*pow(pow(pow(delta-1.0,2),1.0/(2.0*betai)-1.0),2)+Ai*theta*4.0/betai*(1.0/(2.0*betai)-1.0)*pow(pow(delta-1.0,2),1.0/(2.0*betai)-2.0));
long double dDELTAbi2_dDelta2=bi*(pow(DELTA,bi-1.0)*dDELTA2_dDelta2+(bi-1.0)*pow(DELTA,bi-2.0)*pow(dDELTA_dDelta,2));
// At critical point, DELTA is 0, and 1/0^n is undefined
if (fabs(DELTA) < 10*DBL_EPSILON)
{
@@ -1079,7 +1079,7 @@ long double ResidualHelmholtzNonAnalytic::dDelta_dTau(const long double &tau, co
long double dPSI2_dDelta_dTau=4.0*Ci*Di*(delta-1.0)*(tau-1.0)*PSI;
long double dDELTAbi2_dDelta_dTau=-Ai*bi*2.0/betai*pow(DELTA,bi-1.0)*(delta-1.0)*pow(pow(delta-1.0,2),1.0/(2.0*betai)-1.0)-2.0*theta*bi*(bi-1.0)*pow(DELTA,bi-2.0)*dDELTA_dDelta;
long double dPSI_dTau=-2.0*Di*(tau-1.0)*PSI;
long double dDELTAbi_dTau=-2.0*theta*bi*pow(DELTA,bi-1.0);
@@ -1135,12 +1135,12 @@ long double ResidualHelmholtzNonAnalytic::dDelta3(const long double &tau, const
long double dPSI2_dDelta2=(2.0*Ci*pow(delta-1.0,2)-1.0)*2.0*Ci*PSI;
long double dDELTA2_dDelta2=1.0/(delta-1.0)*dDELTA_dDelta+pow(delta-1.0,2)*(4.0*Bi*ai*(ai-1.0)*pow(pow(delta-1.0,2),ai-2.0)+2.0*pow(Ai/betai,2)*pow(pow(pow(delta-1.0,2),1.0/(2.0*betai)-1.0),2)+Ai*theta*4.0/betai*(1.0/(2.0*betai)-1.0)*pow(pow(delta-1.0,2),1.0/(2.0*betai)-2.0));
long double dDELTAbi2_dDelta2=bi*(pow(DELTA,bi-1.0)*dDELTA2_dDelta2+(bi-1.0)*pow(DELTA,bi-2.0)*pow(dDELTA_dDelta,2));
long double dtheta_dDelta = Ai/(2*betai)*pow(pow(delta-1,2),1/(2*betai)-1)*2*(delta-1);
long double dPSI3_dDelta3=2.0*Ci*PSI*(-4*Ci*Ci*pow(delta-1.0,3)+6*Ci*(delta-1));
long double PI = 4*Bi*ai*(ai-1)*pow(pow(delta-1,2),ai-2)+2*pow(Ai/betai,2)*pow(pow(delta-1,2),1/betai-2)+4*Ai*theta/betai*(1/(2*betai)-1)*pow(pow(delta-1,2),1/(2*betai)-2);
long double dPI_dDelta = -8*Bi*ai*(ai-1)*(ai-2)*pow(pow(delta-1,2),ai-5.0/2.0)-8*pow(Ai/betai,2)*(1/(2*betai)-1)*pow(pow(delta-1,2),1/betai-5.0/2.0)-(8*Ai*theta)/betai*(1/(2*betai)-1)*(1/(2*betai)-2)*pow(pow(delta-1,2),1/(2*betai)-5.0/2.0)+4*Ai/betai*(1/(2*betai)-1)*pow(pow(delta-1,2),1/(2*betai)-2)*dtheta_dDelta;
long double dDELTA3_dDelta3 = 1/(delta-1)*dDELTA2_dDelta2-1/pow(delta-1,2)*dDELTA_dDelta+pow(delta-1,2)*dPI_dDelta+2*(delta-1)*PI;
long double dDELTA3_dDelta3 = 1/(delta-1)*dDELTA2_dDelta2-1/pow(delta-1,2)*dDELTA_dDelta+pow(delta-1,2)*dPI_dDelta+2*(delta-1)*PI;
long double dDELTAbi3_dDelta3 = bi*(pow(DELTA,bi-1)*dDELTA3_dDelta3+dDELTA2_dDelta2*(bi-1)*pow(DELTA,bi-2)*dDELTA_dDelta+(bi-1)*(pow(DELTA,bi-2)*2*dDELTA_dDelta*dDELTA2_dDelta2+pow(dDELTA_dDelta,2)*(bi-2)*pow(DELTA,bi-3)*dDELTA_dDelta));
// At critical point, DELTA is 0, and 1/0^n is undefined
@@ -1174,7 +1174,7 @@ long double ResidualHelmholtzNonAnalytic::dDelta_dTau2(const long double &tau, c
long double dDELTAbi_dTau=-2.0*theta*bi*pow(DELTA,bi-1.0);
long double dPSI_dTau=-2.0*Di*(tau-1.0)*PSI;
long double dtheta_dDelta = Ai/(2*betai)*pow(pow(delta-1,2),1/(2*betai)-1)*2*(delta-1);
long double dPSI2_dDelta_dTau=4.0*Ci*Di*(delta-1.0)*(tau-1.0)*PSI;
@@ -1239,7 +1239,7 @@ long double ResidualHelmholtzNonAnalytic::dDelta2_dTau(const long double &tau, c
long double dDELTA_dTau = -2*theta;
long double dDELTA2_dDelta_dTau = 2.0*Ai/(betai)*pow(pow(delta-1,2),1.0/(2.0*betai)-0.5);
long double dDELTA3_dDelta2_dTau = 2.0*Ai*(betai-1)/(betai*betai)*pow(pow(delta-1,2),1/(2*betai)-1.0);
long double dDELTAbim1_dTau = (bi-1)*pow(DELTA,bi-2)*dDELTA_dTau;
long double dDELTAbim2_dTau = (bi-2)*pow(DELTA,bi-3)*dDELTA_dTau;
long double Line11 = dDELTAbim1_dTau*dDELTA2_dDelta2 + pow(DELTA,bi-1)*dDELTA3_dDelta2_dTau;
@@ -1288,7 +1288,7 @@ void ResidualHelmholtzSAFTAssociating::to_json(rapidjson::Value &el, rapidjson::
long double ResidualHelmholtzSAFTAssociating::Deltabar(const long double &tau, const long double &delta)
{
return this->g(this->eta(delta))*(exp(this->epsilonbar*tau)-1)*this->kappabar;
}
}
long double ResidualHelmholtzSAFTAssociating::dDeltabar_ddelta__consttau(const long double &tau, const long double &delta)
{
return this->dg_deta(this->eta(delta))*(exp(this->epsilonbar*tau)-1)*this->kappabar*this->vbarn;
@@ -1382,12 +1382,12 @@ long double ResidualHelmholtzSAFTAssociating::d2X_ddelta2(const long double &tau
long double X = this->X(delta, Deltabar);
long double alpha = this->dDeltabar_ddelta__consttau(tau, delta);
long double dalpha_ddelta = this->d2Deltabar_ddelta2__consttau(tau, delta);
long double dX_ddelta_constall = X*X*(2*Deltabar*Deltabar*X-alpha)/pow(2*Deltabar*delta*X+1,2);
long double d_dXddelta_dX = -(Deltabar+delta*alpha)*2*(Deltabar*delta*X*X+X)/pow(2*Deltabar*delta*X+1,2);
long double d_dXddelta_dDeltabar = X*X*(2*delta*delta*X*alpha-1)/pow(2*Deltabar*delta*X+1,2);
long double d_dXddelta_dalpha = -delta*X*X/(2*Deltabar*delta*X+1);
long double dX_dDeltabar = this->dX_dDeltabar__constdelta(delta, Deltabar);
long double dX_ddelta = this->dX_ddelta__constDeltabar(delta, Deltabar);
@@ -1397,7 +1397,7 @@ long double ResidualHelmholtzSAFTAssociating::d2X_ddelta2(const long double &tau
+ d_dXddelta_dDeltabar*alpha
+ d_dXddelta_dalpha*dalpha_ddelta);
return val;
}
}
long double ResidualHelmholtzSAFTAssociating::d3X_dtau3(const long double &tau, const long double &delta)
{
long double Delta = this->Deltabar(tau, delta);
@@ -1475,7 +1475,7 @@ long double ResidualHelmholtzSAFTAssociating::d3X_ddelta3(const long double &tau
long double ResidualHelmholtzSAFTAssociating::g(const long double &eta)
{
return 0.5*(2-eta)/pow(1-eta,(int)3);
}
}
long double ResidualHelmholtzSAFTAssociating::dg_deta(const long double &eta)
{
return 0.5*(5-2*eta)/pow(1-eta,(int)4);
@@ -1483,11 +1483,11 @@ long double ResidualHelmholtzSAFTAssociating::dg_deta(const long double &eta)
long double ResidualHelmholtzSAFTAssociating::d2g_deta2(const long double &eta)
{
return 3*(3-eta)/pow(1-eta,(int)5);
}
}
long double ResidualHelmholtzSAFTAssociating::d3g_deta3(const long double &eta)
{
return 6*(7-2*eta)/pow(1-eta,(int)6);
}
}
long double ResidualHelmholtzSAFTAssociating::eta(const long double &delta){
return this->vbarn*delta;
}
@@ -1592,7 +1592,7 @@ void IdealHelmholtzCP0PolyT::to_json(rapidjson::Value &el, rapidjson::Document &
el.AddMember("T0", static_cast<double>(T0), doc.GetAllocator());
}
long double IdealHelmholtzCP0PolyT::base(const long double &tau, const long double &delta) throw()
{
{
double sum=0;
for (std::size_t i = 0; i < N; ++i){
if (fabs(t[i])<10*DBL_EPSILON)
@@ -1680,7 +1680,7 @@ void IdealHelmholtzCP0AlyLee::to_json(rapidjson::Value &el, rapidjson::Document
el.AddMember("T0",static_cast<double>(T0),doc.GetAllocator());
}
long double IdealHelmholtzCP0AlyLee::base(const long double &tau, const long double &delta) throw()
{
{
if (!enabled){ return 0.0;}
return -tau*(anti_deriv_cp0_tau2(tau)-anti_deriv_cp0_tau2(tau0))+(anti_deriv_cp0_tau(tau)-anti_deriv_cp0_tau(tau0));
}
@@ -1722,13 +1722,14 @@ IdealHelmholtzEnthalpyEntropyOffset EnthalpyEntropyOffset;
#ifdef ENABLE_CATCH
#include <math.h>
#include "catch.hpp"
#include "crossplatform_shared_ptr.h"
class HelmholtzConsistencyFixture
{
public:
long double numerical, analytic;
std::tr1::shared_ptr<CoolProp::BaseHelmholtzTerm> Lead, LogTau, IGPower, PlanckEinstein,
std::tr1::shared_ptr<CoolProp::BaseHelmholtzTerm> Lead, LogTau, IGPower, PlanckEinstein,
CP0Constant, CP0PolyT, Gaussian, Lemmon2005, Power, SAFT, NonAnalytic, Exponential, GERG2008;
HelmholtzConsistencyFixture(){
@@ -1750,13 +1751,13 @@ public:
}
CP0Constant.reset(new CoolProp::IdealHelmholtzCP0Constant(4/8.314472,300,250));
{
long double beta[] = {1.24, 0.821, 15.45, 2.21, 437, 0.743},
d[] = {1, 1, 2, 2, 3, 3},
epsilon[] = {0.6734, 0.9239, 0.8636, 1.0507, 0.8482, 0.7522},
eta[] = {0.9667, 1.5154, 1.0591, 1.6642, 12.4856, 0.9662},
long double beta[] = {1.24, 0.821, 15.45, 2.21, 437, 0.743},
d[] = {1, 1, 2, 2, 3, 3},
epsilon[] = {0.6734, 0.9239, 0.8636, 1.0507, 0.8482, 0.7522},
eta[] = {0.9667, 1.5154, 1.0591, 1.6642, 12.4856, 0.9662},
gamma[] = {1.2827, 0.4317, 1.1217, 1.1871, 1.1243, 0.4203},
n[] = {1.2198, -0.4883, -0.0033293, -0.0035387, -0.51172, -0.16882},
t[] = {1, 2.124, 0.4, 3.5, 0.5, 2.7};
t[] = {1, 2.124, 0.4, 3.5, 0.5, 2.7};
Gaussian.reset(new CoolProp::ResidualHelmholtzGaussian(std::vector<long double>(n,n+sizeof(n)/sizeof(n[0])),
std::vector<long double>(d,d+sizeof(d)/sizeof(d[0])),
std::vector<long double>(t,t+sizeof(t)/sizeof(t[0])),
@@ -1767,8 +1768,8 @@ public:
));
}
{
long double d[] = {1, 1, 1, 2, 4, 1, 1, 2, 2, 3, 4, 5, 1, 5, 1, 2, 3, 5},
l[] = {0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 2, 3, 3},
long double d[] = {1, 1, 1, 2, 4, 1, 1, 2, 2, 3, 4, 5, 1, 5, 1, 2, 3, 5},
l[] = {0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 2, 3, 3},
m[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.7, 7, 6},
n[] = {5.28076, -8.67658, 0.7501127, 0.7590023, 0.01451899, 4.777189, -3.330988, 3.775673, -2.290919, 0.8888268, -0.6234864, -0.04127263, -0.08455389, -0.1308752, 0.008344962, -1.532005, -0.05883649, 0.02296658},
t[]= {0.669, 1.05,2.75, 0.956, 1, 2, 2.75, 2.38, 3.37, 3.47, 2.63, 3.45, 0.72, 4.23, 0.2, 4.5, 29, 24};
@@ -1780,9 +1781,9 @@ public:
));
}
{
long double d[] = {1, 1, 1, 3, 7, 1, 2, 5, 1, 1, 4, 2},
l[] = {0, 0, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3},
n[] = {1.0038, -2.7662, 0.42921, 0.081363, 0.00024174, 0.48246, 0.75542, -0.00743, -0.4146, -0.016558, -0.10644, -0.021704},
long double d[] = {1, 1, 1, 3, 7, 1, 2, 5, 1, 1, 4, 2},
l[] = {0, 0, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3},
n[] = {1.0038, -2.7662, 0.42921, 0.081363, 0.00024174, 0.48246, 0.75542, -0.00743, -0.4146, -0.016558, -0.10644, -0.021704},
t[] = {0.25, 1.25, 1.5, 0.25, 0.875, 2.375, 2, 2.125, 3.5, 6.5, 4.75, 12.5};
Power.reset(new CoolProp::ResidualHelmholtzPower(std::vector<long double>(n, n+sizeof(n)/sizeof(n[0])),
std::vector<long double>(d, d+sizeof(d)/sizeof(d[0])),
@@ -1791,13 +1792,13 @@ public:
));
}
{
long double a = 1, epsilonbar = 12.2735737, kappabar = 1.09117041e-05, m = 1.01871348, vbarn = 0.0444215309;
SAFT.reset(new CoolProp::ResidualHelmholtzSAFTAssociating(a,m,epsilonbar,vbarn,kappabar));
}
{
long double n[] = {-0.666422765408, 0.726086323499, 0.0550686686128},
A[] = {0.7, 0.7, 0.7}, B[] = {0.3, 0.3, 1}, C[] = {10, 10, 12.5}, D[] = {275, 275, 275},
A[] = {0.7, 0.7, 0.7}, B[] = {0.3, 0.3, 1}, C[] = {10, 10, 12.5}, D[] = {275, 275, 275},
a[] = {3.5, 3.5, 3}, b[] = {0.875, 0.925, 0.875}, beta[] = {0.3, 0.3, 0.3};
NonAnalytic.reset(new CoolProp::ResidualHelmholtzNonAnalytic(std::vector<long double>(n, n+sizeof(n)/sizeof(n[0])),
std::vector<long double>(a, a+sizeof(a)/sizeof(a[0])),
@@ -1810,10 +1811,10 @@ public:
));
}
{
long double d[] = {2, 2, 2, 0, 0, 0},
g[] = {1.65533788, 1.65533788, 1.65533788, 1.65533788, 1.65533788, 1.65533788},
l[] = {2, 2, 2, 2, 2, 2},
n[] = {-3.821884669859, 8.30345065618981, -4.4832307260286, -1.02590136933231, 2.20786016506394, -1.07889905203761},
long double d[] = {2, 2, 2, 0, 0, 0},
g[] = {1.65533788, 1.65533788, 1.65533788, 1.65533788, 1.65533788, 1.65533788},
l[] = {2, 2, 2, 2, 2, 2},
n[] = {-3.821884669859, 8.30345065618981, -4.4832307260286, -1.02590136933231, 2.20786016506394, -1.07889905203761},
t[] = {3, 4, 5, 3, 4, 5};
Exponential.reset(new CoolProp::ResidualHelmholtzExponential(std::vector<long double>(n, n+sizeof(n)/sizeof(n[0])),
std::vector<long double>(d, d+sizeof(d)/sizeof(n[0])),
@@ -1841,7 +1842,7 @@ public:
}
}
void call(std::string d, std::shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double ddelta)
void call(std::string d, shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double ddelta)
{
if (!d.compare("dTau")){return dTau(term,tau,delta,ddelta);}
else if (!d.compare("dTau2")){return dTau2(term,tau,delta,ddelta);}
@@ -1856,7 +1857,7 @@ public:
throw CoolProp::ValueError("don't understand deriv type");
}
}
std::shared_ptr<CoolProp::BaseHelmholtzTerm> get(std::string t)
shared_ptr<CoolProp::BaseHelmholtzTerm> get(std::string t)
{
if (!t.compare("Lead")){return Lead;}
else if (!t.compare("LogTau")){return LogTau;}
@@ -1876,55 +1877,55 @@ public:
throw CoolProp::ValueError(format("don't understand helmholtz type: %s",t.c_str()));
}
}
void dTau(std::shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double dtau){
void dTau(shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double dtau){
long double term_plus = term->base(tau + dtau, delta);
long double term_minus = term->base(tau - dtau, delta);
numerical = (term_plus - term_minus)/(2*dtau);
analytic = term->dTau(tau, delta);
};
void dTau2(std::shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double dtau){
void dTau2(shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double dtau){
long double term_plus = term->dTau(tau + dtau, delta);
long double term_minus = term->dTau(tau - dtau, delta);
numerical = (term_plus - term_minus)/(2*dtau);
analytic = term->dTau2(tau, delta);
};
void dTau3(std::shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double dtau){
void dTau3(shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double dtau){
long double term_plus = term->dTau2(tau + dtau, delta);
long double term_minus = term->dTau2(tau - dtau, delta);
numerical = (term_plus - term_minus)/(2*dtau);
analytic = term->dTau3(tau, delta);
};
void dDelta(std::shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double ddelta){
void dDelta(shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double ddelta){
long double term_plus = term->base(tau, delta + ddelta);
long double term_minus = term->base(tau, delta - ddelta);
numerical = (term_plus - term_minus)/(2*ddelta);
analytic = term->dDelta(tau, delta);
};
void dDelta2(std::shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double ddelta){
void dDelta2(shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double ddelta){
long double term_plus = term->dDelta(tau, delta + ddelta);
long double term_minus = term->dDelta(tau, delta - ddelta);
numerical = (term_plus - term_minus)/(2*ddelta);
analytic = term->dDelta2(tau, delta);
};
void dDelta3(std::shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double ddelta){
void dDelta3(shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double ddelta){
long double term_plus = term->dDelta2(tau, delta + ddelta);
long double term_minus = term->dDelta2(tau, delta - ddelta);
numerical = (term_plus - term_minus)/(2*ddelta);
analytic = term->dDelta3(tau, delta);
};
void dDelta_dTau(std::shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double ddelta){
void dDelta_dTau(shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double ddelta){
long double term_plus = term->dTau(tau, delta + ddelta);
long double term_minus = term->dTau(tau, delta - ddelta);
numerical = (term_plus - term_minus)/(2*ddelta);
analytic = term->dDelta_dTau(tau, delta);
};
void dDelta_dTau2(std::shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double ddelta){
void dDelta_dTau2(shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double ddelta){
long double term_plus = term->dTau2(tau, delta + ddelta);
long double term_minus = term->dTau2(tau, delta - ddelta);
numerical = (term_plus - term_minus)/(2*ddelta);
analytic = term->dDelta_dTau2(tau, delta);
};
void dDelta2_dTau(std::shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double ddelta){
void dDelta2_dTau(shared_ptr<CoolProp::BaseHelmholtzTerm> term, long double tau, long double delta, long double ddelta){
long double term_plus = term->dDelta_dTau(tau, delta + ddelta);
long double term_minus = term->dDelta_dTau(tau, delta - ddelta);
numerical = (term_plus - term_minus)/(2*ddelta);

200
src/HumidAirProp.cpp
View File

@@ -10,8 +10,7 @@
#include "CoolPropTools.h"
#include "Ice.h"
#include "CoolProp.h"
#include "crossplatform_shared_ptr.h"
#include <stdlib.h>
#include "math.h"
@@ -20,6 +19,7 @@
#include <string.h>
#include <iostream>
std::tr1::shared_ptr<CoolProp::AbstractState> Water, Air;
namespace HumidAir
@@ -105,19 +105,19 @@ void UseVirialCorrelations(int flag)
if (flag==0 || flag==1)
{
FlagUseVirialCorrelations=flag;
}
}
else
{
printf("UseVirialCorrelations takes an integer, either 0 (no) or 1 (yes)\n");
}
}
void UseIsothermCompressCorrelation(int flag)
{
if (flag==0 || flag==1)
{
FlagUseIsothermCompressCorrelation=flag;
}
}
else
{
printf("UseIsothermCompressCorrelation takes an integer, either 0 (no) or 1 (yes)\n");
@@ -128,7 +128,7 @@ void UseIdealGasEnthalpyCorrelations(int flag)
if (flag==0 || flag==1)
{
FlagUseIdealGasEnthalpyCorrelations=flag;
}
}
else
{
printf("UseIdealGasEnthalpyCorrelations takes an integer, either 0 (no) or 1 (yes)\n");
@@ -154,7 +154,7 @@ static double Brent_HAProps_T(char *OutputName, char *Input1Name, double Input1,
this->TargetVal = TargetVal;
};
~BrentSolverResids(){};
double call(double T){
return HAPropsSI(OutputName,(char *)"T",T,Input1Name,Input1,Input2Name,Input2)-TargetVal;
}
@@ -231,7 +231,7 @@ static double _B_aw(double T)
static double _dB_aw_dT(double T)
{
check_fluid_instantiation();
// Returns value in m^3/mol
// Returns value in m^3/mol
double a[]={0,0.665687e2,-0.238834e3,-0.176755e3};
double b[]={0,-0.237,-1.048,-3.183};
double rhobarstar=1000,Tstar=100;
@@ -301,11 +301,11 @@ static double B_m(double T, double psi_w)
double B_aa,B_ww,B_aw;
if (FlagUseVirialCorrelations==1)
{
B_aa=-0.000721183853646 +1.142682674467e-05*T -8.838228412173e-08*pow(T,2)
+4.104150642775e-10*pow(T,3) -1.192780880645e-12*pow(T,4) +2.134201312070e-15*pow(T,5)
B_aa=-0.000721183853646 +1.142682674467e-05*T -8.838228412173e-08*pow(T,2)
+4.104150642775e-10*pow(T,3) -1.192780880645e-12*pow(T,4) +2.134201312070e-15*pow(T,5)
-2.157430412913e-18*pow(T,6) +9.453830907795e-22*pow(T,7);
B_ww=-10.8963128394 +2.439761625859e-01*T -2.353884845100e-03*pow(T,2)
+1.265864734412e-05*pow(T,3) -4.092175700300e-08*pow(T,4) +7.943925411344e-11*pow(T,5)
B_ww=-10.8963128394 +2.439761625859e-01*T -2.353884845100e-03*pow(T,2)
+1.265864734412e-05*pow(T,3) -4.092175700300e-08*pow(T,4) +7.943925411344e-11*pow(T,5)
-8.567808759123e-14*pow(T,6) +3.958203548563e-17*pow(T,7);
}
else
@@ -313,7 +313,7 @@ static double B_m(double T, double psi_w)
B_aa = B_Air(T); // [m^3/mol]
B_ww = B_Water(T); // [m^3/mol]
}
B_aw=_B_aw(T); // [m^3/mol]
return pow(1-psi_w,2)*B_aa+2*(1-psi_w)*psi_w*B_aw+psi_w*psi_w*B_ww;
}
@@ -342,11 +342,11 @@ static double C_m(double T, double psi_w)
double C_aaa,C_www,C_aww,C_aaw;
if (FlagUseVirialCorrelations)
{
C_aaa=1.29192158975e-08 -1.776054020409e-10*T +1.359641176409e-12*pow(T,2)
-6.234878717893e-15*pow(T,3) +1.791668730770e-17*pow(T,4) -3.175283581294e-20*pow(T,5)
C_aaa=1.29192158975e-08 -1.776054020409e-10*T +1.359641176409e-12*pow(T,2)
-6.234878717893e-15*pow(T,3) +1.791668730770e-17*pow(T,4) -3.175283581294e-20*pow(T,5)
+3.184306136120e-23*pow(T,6) -1.386043640106e-26*pow(T,7);
C_www=-0.580595811134 +1.365952762696e-02*T -1.375986293288e-04*pow(T,2)
+7.687692259692e-07*pow(T,3) -2.571440816920e-09*pow(T,4) +5.147432221082e-12*pow(T,5)
C_www=-0.580595811134 +1.365952762696e-02*T -1.375986293288e-04*pow(T,2)
+7.687692259692e-07*pow(T,3) -2.571440816920e-09*pow(T,4) +5.147432221082e-12*pow(T,5)
-5.708156494894e-15*pow(T,6) +2.704605721778e-18*pow(T,7);
}
else
@@ -362,7 +362,7 @@ static double C_m(double T, double psi_w)
static double dC_m_dT(double T, double psi_w)
{
// dCm_dT has units of m^6/mol^2/K
double Tj,tau_Air,tau_Water,dC_dT_aaa,dC_dT_www,dC_dT_aww,dC_dT_aaw;
// NDG for fluid EOS for virial terms
Tj=132.6312;
@@ -391,7 +391,7 @@ static double HenryConstant(double T)
{
// Result has units of 1/Pa
double p_ws,beta_N2,beta_O2,beta_Ar,beta_a,tau,Tr,Tc=647.096;
Tr=T/Tc;
Tr=T/Tc;
tau=1-Tr;
Water->update(CoolProp::QT_INPUTS, 1.0, T);
p_ws = Water->keyed_output(CoolProp::iP); //[Pa]
@@ -404,7 +404,7 @@ static double HenryConstant(double T)
double isothermal_compressibility(double T, double p)
{
double k_T;
if (T> 273.16)
{
if (FlagUseIsothermCompressCorrelation)
@@ -421,7 +421,7 @@ double isothermal_compressibility(double T, double p)
else
{
k_T = IsothermCompress_Ice(T,p); //[1/Pa]
}
}
return k_T;
}
double f_factor(double T, double p)
@@ -432,7 +432,7 @@ double f_factor(double T, double p)
double p_ws,tau_Air,tau_Water,B_aa,B_aw,B_ww,C_aaa,C_aaw,C_aww,C_www,
line1,line2,line3,line4,line5,line6,line7,line8,k_T,beta_H,LHS,RHS,psi_ws,
vbar_ws;
// Saturation pressure [Pa]
if (T>273.16)
{
@@ -449,33 +449,33 @@ double f_factor(double T, double p)
beta_H = 0;
vbar_ws = dg_dp_Ice(T,p)*MM_Water(); //[m^3/mol]
}
k_T = isothermal_compressibility(T,p); //[1/Pa]
// Hermann: In the iteration process of the enhancement factor in Eq. (3.25), k_T is set to zero for pw,s (T) > p.
if (p_ws>p)
{
k_T=0;
beta_H=0;
}
// NDG for fluid EOS for virial terms
Tj=132.6312;
tau_Air=Tj/T;
tau_Water=Water->keyed_output(CoolProp::iT_reducing)/T;
if (FlagUseVirialCorrelations)
{
B_aa=-0.000721183853646 +1.142682674467e-05*T -8.838228412173e-08*pow(T,2)
+4.104150642775e-10*pow(T,3) -1.192780880645e-12*pow(T,4) +2.134201312070e-15*pow(T,5)
B_aa=-0.000721183853646 +1.142682674467e-05*T -8.838228412173e-08*pow(T,2)
+4.104150642775e-10*pow(T,3) -1.192780880645e-12*pow(T,4) +2.134201312070e-15*pow(T,5)
-2.157430412913e-18*pow(T,6) +9.453830907795e-22*pow(T,7);
B_ww=-10.8963128394 +2.439761625859e-01*T -2.353884845100e-03*pow(T,2)
+1.265864734412e-05*pow(T,3) -4.092175700300e-08*pow(T,4) +7.943925411344e-11*pow(T,5)
B_ww=-10.8963128394 +2.439761625859e-01*T -2.353884845100e-03*pow(T,2)
+1.265864734412e-05*pow(T,3) -4.092175700300e-08*pow(T,4) +7.943925411344e-11*pow(T,5)
-8.567808759123e-14*pow(T,6) +3.958203548563e-17*pow(T,7);
C_aaa=1.29192158975e-08 -1.776054020409e-10*T +1.359641176409e-12*pow(T,2)
-6.234878717893e-15*pow(T,3) +1.791668730770e-17*pow(T,4) -3.175283581294e-20*pow(T,5)
C_aaa=1.29192158975e-08 -1.776054020409e-10*T +1.359641176409e-12*pow(T,2)
-6.234878717893e-15*pow(T,3) +1.791668730770e-17*pow(T,4) -3.175283581294e-20*pow(T,5)
+3.184306136120e-23*pow(T,6) -1.386043640106e-26*pow(T,7);
C_www=-0.580595811134 +1.365952762696e-02*T -1.375986293288e-04*pow(T,2)
+7.687692259692e-07*pow(T,3) -2.571440816920e-09*pow(T,4) +5.147432221082e-12*pow(T,5)
C_www=-0.580595811134 +1.365952762696e-02*T -1.375986293288e-04*pow(T,2)
+7.687692259692e-07*pow(T,3) -2.571440816920e-09*pow(T,4) +5.147432221082e-12*pow(T,5)
-5.708156494894e-15*pow(T,6) +2.704605721778e-18*pow(T,7);
}
else
@@ -488,7 +488,7 @@ double f_factor(double T, double p)
B_aw = _B_aw(T); //[m^3/mol]
C_aaw = _C_aaw(T); //[m^6/mol^2]
C_aww = _C_aww(T); //[m^6/mol^2]
// Use a little secant loop to find f iteratively
// Start out with a guess value of 1 for f
while ((iter<=3 || change>eps) && iter<100)
@@ -496,12 +496,12 @@ double f_factor(double T, double p)
if (iter==1){x1=1.00; f=x1;}
if (iter==2){x2=1.00+0.000001; f=x2;}
if (iter>2) {f=x2;}
// Left-hand-side of Equation 3.25
LHS=log(f);
// Eqn 3.24
psi_ws=f*p_ws/p;
// All the terms forming the RHS of Eqn 3.25
line1=((1+k_T*p_ws)*(p-p_ws)-k_T*(p*p-p_ws*p_ws)/2.0)/(Rbar*T)*vbar_ws+log(1-beta_H*(1-psi_ws)*p);
line2=pow(1-psi_ws,2)*p/(Rbar*T)*B_aa-2*pow(1-psi_ws,2)*p/(Rbar*T)*B_aw-(p-p_ws-pow(1-psi_ws,2)*p)/(Rbar*T)*B_ww;
@@ -512,7 +512,7 @@ double f_factor(double T, double p)
line7=(6*pow(1-psi_ws,2)*psi_ws*psi_ws*p*p)/pow(Rbar*T,2)*B_ww*B_aw-(3*pow(1-psi_ws,4)*p*p)/(2*pow(Rbar*T,2))*B_aa*B_aa;
line8=-(2*pow(1-psi_ws,2)*psi_ws*(-2+3*psi_ws)*p*p)/pow(Rbar*T,2)*B_aw*B_aw-(p_ws*p_ws-(4-3*psi_ws)*pow(psi_ws,3)*p*p)/(2*pow(Rbar*T,2))*B_ww*B_ww;
RHS=line1+line2+line3+line4+line5+line6+line7+line8;
if (iter==1){y1=LHS-RHS;}
if (iter>1)
{
@@ -554,9 +554,9 @@ double Viscosity(double T, double p, double psi_w)
{
/*
Using the method of:
P.T. Tsilingiris, 2009, Thermophysical and transport properties of humid air at temperature range between 0 and 100 oC, Energy Conversion and Management, 49, 1098-1010
but using the detailed measurements for pure fluid from IAPWS formulations
*/
double mu_a,mu_w,Phi_av,Phi_va,Ma,Mw;
@@ -576,9 +576,9 @@ double Conductivity(double T, double p, double psi_w)
{
/*
Using the method of:
P.T. Tsilingiris, 2009, Thermophysical and transport properties of humid air at temperature range between 0 and 100 oC, Energy Conversion and Management, 49, 1098-1010
but using the detailed measurements for pure fluid from IAPWS formulations
*/
double mu_a,mu_w,k_a,k_w,Phi_av,Phi_va,Ma,Mw;
@@ -602,14 +602,14 @@ double MolarVolume(double T, double p, double psi_w)
// Output in m^3/mol
int iter;
double v_bar0, v_bar=0, R_bar=8.314472,x1=0,x2=0,x3,y1=0,y2,resid,eps,Bm,Cm;
// -----------------------------
// Iteratively find molar volume
// -----------------------------
// Start by assuming it is an ideal gas to get initial guess
v_bar0=R_bar*T/p;
//Bring outside the loop since not a function of v_bar
Bm=B_m(T,psi_w);
Cm=C_m(T,psi_w);
@@ -620,10 +620,10 @@ double MolarVolume(double T, double p, double psi_w)
if (iter==1){x1=v_bar0; v_bar=x1;}
if (iter==2){x2=v_bar0+0.000001; v_bar=x2;}
if (iter>2) {v_bar=x2;}
// want v_bar in m^3/mol and R_bar in J/mol-K
resid = (p-(R_bar)*T/v_bar*(1+Bm/v_bar+Cm/(v_bar*v_bar)))/p;
if (iter==1){y1=resid;}
if (iter>1)
{
@@ -653,7 +653,7 @@ double IdealGasMolarEntropy_Water(double T, double p)
tau = Water->keyed_output(CoolProp::iT_reducing)/T;
Water->update(CoolProp::DmolarT_INPUTS,p/(R_bar*T),T);
sbar_w = R_bar*(tau*Water->keyed_output(CoolProp::idalpha0_dtau_constdelta)-Water->keyed_output(CoolProp::ialpha0)); //[kJ/kmol/K]
return sbar_w;
return sbar_w;
}
double IdealGasMolarEnthalpy_Air(double T, double vmolar)
{
@@ -672,15 +672,15 @@ double IdealGasMolarEnthalpy_Air(double T, double vmolar)
double IdealGasMolarEntropy_Air(double T, double vmolar_a)
{
double sbar_0_Lem, tau, sbar_a, R_bar_Lemmon = 8.314510, T0=273.15, p0=101325, vmolar_a_0;
// Ideal-Gas contribution to entropy of air
sbar_0_Lem=-196.1375815; //[J/mol/K]
// Tj and rhoj are given by 132.6312 and 302.5507652 respectively
tau = 132.6312/T; //[no units]
vmolar_a_0 = R_bar_Lemmon*T0/p0; //[m^3/mol]
Air->update(CoolProp::DmolarT_INPUTS,1/vmolar_a_0,T);
sbar_a=sbar_0_Lem+R_bar_Lemmon*(tau*Air->keyed_output(CoolProp::idalpha0_dtau_constdelta)-Air->keyed_output(CoolProp::ialpha0))+R_bar_Lemmon*log(vmolar_a/vmolar_a_0); //[J/mol/K]
@@ -691,9 +691,9 @@ double IdealGasMolarEntropy_Air(double T, double vmolar_a)
double MolarEnthalpy(double T, double p, double psi_w, double vmolar)
{
// In units of kJ/kmol
// vbar (molar volume) in m^3/kg
double hbar_0, hbar_a, hbar_w, hbar, R_bar=8.314472;
// ----------------------------------------
// Enthalpy
@@ -701,7 +701,7 @@ double MolarEnthalpy(double T, double p, double psi_w, double vmolar)
// Constant for enthalpy
// Not clear why getting rid of this term yields the correct values in the table, but enthalpies are equal to an additive constant, so not a big deal
hbar_0=0.0;//2.924425468; //[kJ/kmol]
if (FlagUseIdealGasEnthalpyCorrelations){
hbar_w = 2.7030251618E-03*T*T + 3.1994361015E+01*T + 3.6123174929E+04;
hbar_a = 9.2486716590E-04*T*T + 2.8557221776E+01*T - 7.8616129429E+03;
@@ -710,7 +710,7 @@ double MolarEnthalpy(double T, double p, double psi_w, double vmolar)
hbar_w = IdealGasMolarEnthalpy_Water(T, vmolar);
hbar_a = IdealGasMolarEnthalpy_Air(T, vmolar);
}
hbar = hbar_0+(1-psi_w)*hbar_a+psi_w*hbar_w+R_bar*T*((B_m(T,psi_w)-T*dB_m_dT(T,psi_w))/vmolar+(C_m(T,psi_w)-T/2.0*dC_m_dT(T,psi_w))/(vmolar*vmolar));
return hbar; //[J/mol]
}
@@ -730,7 +730,7 @@ double MolarEntropy(double T, double p, double psi_w, double v_bar)
// Serious typo in RP-1485 - should use total pressure rather than
// reference pressure in density calculation for water vapor molar entropy
// vbar (molar volume) in m^3/mol
double x1=0,x2=0,x3=0,y1=0,y2=0,eps=1e-8,f=999,R_bar_Lem=8.314510;
int iter=1;
@@ -741,7 +741,7 @@ double MolarEntropy(double T, double p, double psi_w, double v_bar)
// Calculate vbar_a, the molar volume of dry air
// B_m, C_m, etc. functions take care of the units
Baa = B_m(T,0);
Baa = B_m(T,0);
B = B_m(T,psi_w);
dBdT = dB_m_dT(T,psi_w);
Caaa = C_m(T,0);
@@ -749,7 +749,7 @@ double MolarEntropy(double T, double p, double psi_w, double v_bar)
dCdT = dC_m_dT(T,psi_w);
vbar_a_guess = R_bar_Lem*T/p; //[m^3/mol] since p in [Pa]
while ((iter<=3 || fabs(f)>eps) && iter<100)
{
if (iter==1){x1=vbar_a_guess; vbar_a=x1;}
@@ -766,7 +766,7 @@ double MolarEntropy(double T, double p, double psi_w, double v_bar)
iter=iter+1;
if (iter>100){ return _HUGE; }
}
if (FlagUseIdealGasEnthalpyCorrelations){
std::cout << "Not implemented" << std::endl;
}
@@ -788,7 +788,7 @@ double DewpointTemperature(double T, double p, double psi_w)
int iter;
double p_w,eps,resid,Tdp=0,x1=0,x2=0,x3,y1=0,y2,T0;
double p_ws_dp,f_dp;
// Make sure it isn't dry air, return an impossible temperature otherwise
if ((1-psi_w)<1e-16)
{
@@ -798,9 +798,9 @@ double DewpointTemperature(double T, double p, double psi_w)
// Iteratively find the dewpoint temperature
// ------------------------------------------
// The highest dewpoint temperature possible is the dry-bulb temperature.
// The highest dewpoint temperature possible is the dry-bulb temperature.
// When they are equal, the air is saturated (R=1)
p_w = psi_w*p;
// 0.61165... is the triple point pressure of water in kPa
@@ -813,14 +813,14 @@ double DewpointTemperature(double T, double p, double psi_w)
}
// A good guess for Tdp is that enhancement factor is unity, which yields
// p_w_s = p_w, and get guess for T from saturation temperature
iter=1; eps=1e-8; resid=999;
while ((iter<=3 || fabs(resid)>eps) && iter<100)
{
if (iter==1){x1 = T0; Tdp=x1;}
if (iter==2){x2 = x1 + 0.1; Tdp=x2;}
if (iter>2) {Tdp=x2;}
if (Tdp >= 273.16)
{
// Saturation pressure at dewpoint [kPa]
@@ -836,7 +836,7 @@ double DewpointTemperature(double T, double p, double psi_w)
f_dp=f_factor(Tdp,p);
// Error between target and actual pressure [kPa]
resid=p_w-p_ws_dp*f_dp;
if (iter==1){y1=resid;}
if (iter>1)
{
@@ -875,7 +875,7 @@ public:
if (Twb > 273.16)
{
// Saturation pressure at wetbulb temperature [Pa]
Water->update(CoolProp::QT_INPUTS,0,Twb);
Water->update(CoolProp::QT_INPUTS,0,Twb);
p_ws_wb= Water->keyed_output(CoolProp::iP);
}
else
@@ -883,7 +883,7 @@ public:
// Sublimation pressure at wetbulb temperature [kPa]
p_ws_wb=psub_Ice(Twb);
}
// Vapor pressure
p_s_wb = f_wb*p_ws_wb;
// wetbulb humidity ratio
@@ -902,7 +902,7 @@ public:
h_w=h_Ice(Twb,_p);
}
// Mole masses of wetbulb and humid air
M_ha_wb=MM_Water()*psi_wb+(1-psi_wb)*28.966;
v_bar_wb=MolarVolume(Twb,_p,psi_wb);
double RHS = (MolarEnthalpy(Twb,_p,psi_wb,v_bar_wb)*(1+W_s_wb)/M_ha_wb+(_W-W_s_wb)*h_w);
@@ -934,7 +934,7 @@ double WetbulbTemperature(double T, double p, double psi_w)
// ------------------------------------------
// Iteratively find the wetbulb temperature
// ------------------------------------------
//
//
// If the temperature is less than the saturation temperature of water
// for the given atmospheric pressure, the highest wetbulb temperature that is possible is the dry bulb
// temperature
@@ -953,17 +953,17 @@ double WetbulbTemperature(double T, double p, double psi_w)
WetBulbSolver WBS(T, p, psi_w);
std::string errstr;
double return_val;
try{
return_val = Secant(WBS,Tmax,0.0001,1e-8,50,errstr);
// Solution obtained is out of range (T>Tmax)
if (return_val > Tmax) {throw CoolProp::ValueError();}
}
catch(std::exception &)
{
// The lowest wetbulb temperature that is possible for a given dry bulb temperature
// The lowest wetbulb temperature that is possible for a given dry bulb temperature
// is the saturated air temperature which yields the enthalpy of dry air at dry bulb temperature
try{
@@ -980,7 +980,7 @@ double WetbulbTemperature(double T, double p, double psi_w)
return_val = _HUGE;
}
}
return return_val;
return return_val;
}
static int Name2Type(const char *Name)
{
@@ -1023,7 +1023,7 @@ int TypeMatch(int TypeCode, const char *Input1Name, const char *Input2Name, cons
double MoleFractionWater(double T, double p, int HumInput, double InVal)
{
double p_ws,f,W,epsilon=0.621945,Tdp,p_ws_dp,f_dp,p_w_dp,p_s,RH;
if (HumInput==GIVEN_HUMRAT) //(2)
{
W=InVal;
@@ -1041,7 +1041,7 @@ double MoleFractionWater(double T, double p, int HumInput, double InVal)
{
// Sublimation pressure [Pa]
p_ws=psub_Ice(T);
}
// Enhancement Factor [-]
f=f_factor(T,p);
@@ -1107,21 +1107,21 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1,
{
// Add a check to make sure that Air and Water fluid states have been properly instantiated
check_fluid_instantiation();
int In1Type, In2Type, In3Type,iT,iW,iTdp,iRH,ip,Type1,Type2;
double vals[3],p,T,RH,W,Tdp,psi_w,M_ha,v_bar,h_bar,s_bar,MainInputValue,SecondaryInputValue,T_guess;
double Value1,Value2,W_guess;
char MainInputName[100], SecondaryInputName[100],Name1[100],Name2[100];
vals[0]=Input1;
vals[1]=Input2;
vals[2]=Input3;
// First figure out what kind of inputs you have, convert names to Macro expansions
In1Type=Name2Type(Input1Name);
In2Type=Name2Type(Input2Name);
In3Type=Name2Type(Input3Name);
// Pressure must be included
ip=TypeMatch(GIVEN_P,Input1Name,Input2Name,Input3Name);
if (ip>0)
@@ -1132,13 +1132,13 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1,
// -----------------------------------------------------------------------------------------------------
// Check whether the remaining values give explicit solution for mole fraction of water - nice and fast
// -----------------------------------------------------------------------------------------------------
// Find the codes if they are there
iT= TypeMatch(GIVEN_T,Input1Name,Input2Name,Input3Name);
iRH= TypeMatch(GIVEN_RH,Input1Name,Input2Name,Input3Name);
iW= TypeMatch(GIVEN_HUMRAT,Input1Name,Input2Name,Input3Name);
iTdp= TypeMatch(GIVEN_TDP,Input1Name,Input2Name,Input3Name);
if (iT>0) // Found T (or alias) as an input
{
T=vals[iT-1];
@@ -1189,9 +1189,9 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1,
else
{
// Need to iterate to find dry bulb temperature since temperature is not provided
// Pick one input, and alter T to match the other input
// Get the variables and their values that are NOT pressure for simplicity
// because you know you need pressure as an input and you already have
// its value in variable p
@@ -1219,17 +1219,17 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1,
else{
return _HUGE;
}
// Get the integer type codes
Type1=Name2Type(Name1);
Type2=Name2Type(Name2);
// First, if one of the inputs is something that can potentially yield
// First, if one of the inputs is something that can potentially yield
// an explicit solution at a given iteration of the solver, use it
if (Type1==GIVEN_RH || Type1==GIVEN_HUMRAT || Type1==GIVEN_TDP)
{
// First input variable is a "nice" one
// MainInput is the one that you are using in the call to HAProps
MainInputValue=Value1;
strcpy(MainInputName,Name1);
@@ -1240,7 +1240,7 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1,
else if (Type2==GIVEN_RH || Type2==GIVEN_HUMRAT || Type2==GIVEN_TDP)
{
// Second input variable is a "nice" one
// MainInput is the one that you are using in the call to HAProps
MainInputValue=Value2;
strcpy(MainInputName,Name2);
@@ -1266,8 +1266,8 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1,
T = Secant_HAProps_T(SecondaryInputName,(char *)"P",p,MainInputName,MainInputValue,SecondaryInputValue,T_guess);
double val = HAPropsSI(SecondaryInputName,(char *)"T",T,(char *)"P",p,MainInputName,MainInputValue);
if (!ValidNumber(T) || !ValidNumber(val) || !(T_min < T && T < T_max) || fabs(val-SecondaryInputValue)>1e-6)
{
throw CoolProp::ValueError();
{
throw CoolProp::ValueError();
}
else
{
@@ -1276,13 +1276,13 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1,
}
catch (std::exception &){};
}
if (T < 0) // No solution found using secant
{
// Use the Brent's method solver to find T
T = Brent_HAProps_T(SecondaryInputName,(char *)"P",p,MainInputName,MainInputValue,SecondaryInputValue,T_min,T_max);
}
// If you want the temperature, return it
if (Name2Type(OutputName)==GIVEN_T)
return T;
@@ -1294,7 +1294,7 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1,
}
}
M_ha=(1-psi_w)*0.028966+MM_Water()*psi_w; //[kg_ha/mol_ha]
// -----------------------------------------------------------------
// Calculate and return the desired value for known set of T,p,psi_w
// -----------------------------------------------------------------
@@ -1395,16 +1395,16 @@ double HAPropsSI(const char *OutputName, const char *Input1Name, double Input1,
double HAProps_Aux(const char* Name,double T, double p, double W, char *units)
{
// This function provides some things that are not usually needed, but could be interesting for debug purposes.
// Requires W since it is nice and fast and always defined. Put a dummy value if you want something that doesn't use humidity
// Takes temperature, pressure, and humidity ratio W as inputs;
double psi_w,Tj,tau_Water,tau_Air,B_aa,C_aaa,B_ww,C_www,B_aw,C_aaw,C_aww,v_bar;
Tj=132.6312;
tau_Air=Tj/T;
tau_Water=Water->keyed_output(CoolProp::iT_critical)/T;
try{
if (!strcmp(Name,"Baa"))
{
@@ -1589,8 +1589,8 @@ double HAProps_Aux(const char* Name,double T, double p, double W, char *units)
}
double cair_sat(double T)
{
// Air saturation specific heat
// Based on a correlation from EES, good from 250K to 300K.
// Air saturation specific heat
// Based on a correlation from EES, good from 250K to 300K.
// No error bound checking is carried out
// T: [K]
// cair_s: [kJ/kg-K]
@@ -1674,7 +1674,7 @@ TEST_CASE("Check HA Virials from Table A.2.1","[RP1485]")
CHECK(fabs(HumidAir::_C_aww(0+273.15)/(-224234/1e12)-1) < 1e-3);
CHECK(fabs(HumidAir::_C_aww(200+273.15)/(-8436.5/1e12)-1) < 1e-3);
CHECK(fabs(HumidAir::_C_aww(350+273.15)/(-2486.9/1e12)-1) < 1e-3);
}
}
}
TEST_CASE("Enhancement factor from Table A.3","[RP1485]")
{
@@ -1800,7 +1800,7 @@ TEST_CASE_METHOD(HAPropsConsistencyFixture, "ASHRAE RP1485 Tables", "[RP1485]")
CHECK(fabs(actual/expected-1) < 0.01);
}
}
}
#endif /* CATCH_ENABLED */

1
src/PolyMath.cpp
View File

@@ -998,7 +998,6 @@ double BaseExponential::expval(const std::vector< std::vector<double> > &coeffic
#include <math.h>
#include "catch.hpp"
class PolynomialConsistencyFixture {
public:
CoolProp::BasePolynomial poly;

1
src/SpeedTest.cpp
View File

@@ -2,6 +2,7 @@
#include "SpeedTest.h"
#include "AbstractState.h"
#include "DataStructures.h"
#include "crossplatform_shared_ptr.h"
#include <time.h>

32
src/Tests/CoolProp-Tests.cpp
View File

@@ -9,6 +9,7 @@
#if defined(ENABLE_CATCH)
#include "crossplatform_shared_ptr.h"
#include "catch.hpp"
namespace TransportValidation{
@@ -22,7 +23,7 @@ public:
vel(std::string fluid, std::string in1, double v1, std::string in2, double v2, std::string out, double expected, double tol)
{
this->in1 = in1; this->in2 = in2; this->fluid = fluid;
this->v1 = v1; this->v2 = v2; this->expected = expected;
this->v1 = v1; this->v2 = v2; this->expected = expected;
this->tol = tol;
};
};
@@ -218,7 +219,7 @@ public:
void set_backend(std::string backend, std::string fluid_name){
pState.reset(CoolProp::AbstractState::factory(backend, fluid_name));
}
void set_pair(std::string &in1, double v1, std::string &in2, double v2){
void set_pair(std::string &in1, double v1, std::string &in2, double v2){
double o1, o2;
long iin1 = CoolProp::get_parameter_index(in1);
long iin2 = CoolProp::get_parameter_index(in2);
@@ -470,24 +471,24 @@ TEST_CASE_METHOD(TransportValidationFixture, "Compare thermal conductivities aga
static int inputs[] = {
CoolProp::DmolarT_INPUTS,
//CoolProp::SmolarT_INPUTS,
//CoolProp::HmolarT_INPUTS,
//CoolProp::HmolarT_INPUTS,
//CoolProp::TUmolar_INPUTS,
CoolProp::DmolarP_INPUTS,
CoolProp::DmolarHmolar_INPUTS,
CoolProp::DmolarSmolar_INPUTS,
CoolProp::DmolarP_INPUTS,
CoolProp::DmolarHmolar_INPUTS,
CoolProp::DmolarSmolar_INPUTS,
CoolProp::DmolarUmolar_INPUTS,
/*
CoolProp::HmolarP_INPUTS,
CoolProp::PSmolar_INPUTS,
CoolProp::PUmolar_INPUTS,
CoolProp::PUmolar_INPUTS,
*/
/*
CoolProp::HmolarSmolar_INPUTS,
CoolProp::HmolarUmolar_INPUTS,
CoolProp::SmolarUmolar_INPUTS
CoolProp::HmolarSmolar_INPUTS,
CoolProp::HmolarUmolar_INPUTS,
CoolProp::SmolarUmolar_INPUTS
*/
};
@@ -503,7 +504,7 @@ public:
void set_backend(std::string backend, std::string fluid_name){
pState.reset(CoolProp::AbstractState::factory(backend, fluid_name));
}
void set_pair(int pair){
void set_pair(int pair){
this->pair = pair;
}
void set_TP(long double T, long double p)
@@ -513,14 +514,13 @@ public:
// Start with T,P as inputs, cycle through all the other pairs that are supported
State.update(CoolProp::PT_INPUTS, p, T);
// Set the other state variables
rhomolar = State.rhomolar(); hmolar = State.hmolar(); smolar = State.smolar(); umolar = State.umolar();
}
void get_variables()
{
CoolProp::AbstractState &State = *pState;
switch (pair)
{
/// In this group, T is one of the known inputs, iterate for the other one (easy)
@@ -571,7 +571,7 @@ public:
TEST_CASE_METHOD(ConsistencyFixture, "Test all input pairs for CO2 using all valid backends", "[]")
{
CHECK_NOTHROW(set_backend("HEOS", "CO2"));
int inputsN = sizeof(inputs)/sizeof(inputs[0]);
for (double p = 600000; p < 800000000.0; p *= 5)
{

6
src/Tests/Tests.cpp
View File

@@ -1,15 +1,15 @@
/**
This file includes some testing functions that will get built
This file includes some testing functions that will get built
into the program. Otherwise CTest can be used by removing this file from
the build to avoid double declaration of the main function and
Catch clashing
*/
#include "Tests.h"
#if defined ENABLE_CATCH
#if defined ENABLE_CATCH
#define CATCH_CONFIG_RUNNER
#include "catch.hpp"
static Catch::Session session; // There must be exactly one instance
#endif // ENABLE_CATCH

6
src/main.cxx
View File

@@ -66,7 +66,8 @@ int main()
{
try{
//set_reference_stateS(NBP_refs[i],"RESET");
HelmholtzEOSMixtureBackend HEOS(std::vector<std::string>(1,NBP_refs[i]));
std::vector<std::string> comps(1,NBP_refs[i]);
HelmholtzEOSMixtureBackend HEOS(comps);
HEOS.update(PQ_INPUTS, 101325, 0);
double delta_a1 = HEOS.smass()/(HEOS.gas_constant()/HEOS.molar_mass());
double delta_a2 = -HEOS.hmass()/(HEOS.gas_constant()/HEOS.molar_mass()*HEOS.get_reducing().T);
@@ -81,7 +82,8 @@ int main()
{
try{
//set_reference_stateS(IIR_refs[i],"RESET");
HelmholtzEOSMixtureBackend HEOS(std::vector<std::string>(1,IIR_refs[i]));
std::vector<std::string> comps(1,IIR_refs[i]);
HelmholtzEOSMixtureBackend HEOS(comps);
HEOS.update(QT_INPUTS, 0, 273.15);
double delta_a1 = (HEOS.smass()-1000)/(HEOS.gas_constant()/HEOS.molar_mass());
double delta_a2 = -(HEOS.hmass()-200000)/(HEOS.gas_constant()/HEOS.molar_mass()*HEOS.get_reducing().T);