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

ECS conductivity work - nearly works, just need to sort out the conformal state solver.

Browse Source 
Goodbye to AbstractStateWrapper - can use std::tr1::shared_ptr, much nicer

Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-05-28 19:15:34 +02:00
parent 118d7e9ad2
commit f93aa76210
11 changed files with 287 additions and 159 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -437,6 +437,21 @@ protected:
}
};
void parse_ECS_conductivity(rapidjson::Value &conductivity, CoolPropFluid & fluid)
{
fluid.transport.conductivity_ecs.reference_fluid = cpjson::get_string(conductivity,"reference_fluid");
// Parameters for correction polynomials
fluid.transport.conductivity_ecs.psi_a = cpjson::get_long_double_array(conductivity["psi"]["a"]);
fluid.transport.conductivity_ecs.psi_t = cpjson::get_long_double_array(conductivity["psi"]["t"]);
fluid.transport.conductivity_ecs.psi_rhomolar_reducing = cpjson::get_double(conductivity["psi"],"rhomolar_reducing");
fluid.transport.conductivity_ecs.f_int_a = cpjson::get_long_double_array(conductivity["f_int"]["a"]);
fluid.transport.conductivity_ecs.f_int_t = cpjson::get_long_double_array(conductivity["f_int"]["t"]);
fluid.transport.conductivity_ecs.f_int_T_reducing = cpjson::get_double(conductivity["f_int"],"T_reducing");
fluid.transport.conductivity_using_ECS = true;
}
void parse_ECS_viscosity(rapidjson::Value &viscosity, CoolPropFluid & fluid)
{
fluid.transport.viscosity_ecs.reference_fluid = cpjson::get_string(viscosity,"reference_fluid");
@@ -446,7 +461,7 @@ protected:
fluid.transport.viscosity_ecs.psi_t = cpjson::get_long_double_array(viscosity["psi"]["t"]);
fluid.transport.viscosity_ecs.psi_rhomolar_reducing = cpjson::get_double(viscosity["psi"],"rhomolar_reducing");
fluid.transport.using_ECS = true;
fluid.transport.viscosity_using_ECS = true;
}
/// Parse the transport properties
@@ -635,6 +650,12 @@ protected:
/// Parse the thermal conductivity data
void parse_thermal_conductivity(rapidjson::Value &conductivity, CoolPropFluid & fluid)
{
// If it is using ECS, set ECS parameters and quit
if (conductivity.HasMember("type") && !cpjson::get_string(conductivity, "type").compare("ECS")){
parse_ECS_conductivity(conductivity, fluid);
return;
}
if (conductivity.HasMember("hardcoded")){
std::string target = cpjson::get_string(conductivity, "hardcoded");
if (!target.compare("Water")){

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

@@ -199,9 +199,14 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity(void)
CoolPropFluid &component = *(components[0]);
// Check if using ECS
if (component.transport.using_ECS)
if (component.transport.viscosity_using_ECS)
{
return TransportRoutines::viscosity_ECS(*this, *this);
// Get reference fluid name
std::string fluid_name = component.transport.viscosity_ecs.reference_fluid;
// Get a managed pointer to the reference fluid for ECS
std::tr1::shared_ptr<HelmholtzEOSMixtureBackend> ref_fluid(new HelmholtzEOSMixtureBackend(std::vector<std::string>(1, fluid_name)));
// Get the viscosity using ECS
return TransportRoutines::viscosity_ECS(*this, *(ref_fluid.get()));
}
if (component.transport.hardcoded_viscosity != CoolProp::TransportPropertyData::VISCOSITY_NOT_HARDCODED)
@@ -234,13 +239,42 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity(void)
throw NotImplementedError(format("viscosity not implemented for mixtures"));
}
}
long double HelmholtzEOSMixtureBackend::calc_conductivity_background(void)
{
// Residual part
long double lambda_residual = _HUGE;
switch(components[0]->transport.conductivity_residual.type)
{
case ConductivityResidualVariables::CONDUCTIVITY_RESIDUAL_POLYNOMIAL:
lambda_residual = TransportRoutines::conductivity_residual_polynomial(*this); break;
case ConductivityResidualVariables::CONDUCTIVITY_RESIDUAL_POLYNOMIAL_AND_EXPONENTIAL:
lambda_residual = TransportRoutines::conductivity_residual_polynomial_and_exponential(*this); break;
default:
throw ValueError(format("residual conductivity type [%d] is invalid for fluid %s", components[0]->transport.conductivity_residual.type, name().c_str()));
}
return lambda_residual;
}
long double HelmholtzEOSMixtureBackend::calc_conductivity(void)
{
if (is_pure_or_pseudopure)
{
if (components[0]->transport.hardcoded_conductivity != CoolProp::TransportPropertyData::CONDUCTIVITY_NOT_HARDCODED)
// Get a reference for code cleanness
CoolPropFluid &component = *(components[0]);
// Check if using ECS
if (component.transport.conductivity_using_ECS)
{
switch(components[0]->transport.hardcoded_conductivity)
// Get reference fluid name
std::string fluid_name = component.transport.conductivity_ecs.reference_fluid;
// Get a managed pointer to the reference fluid for ECS
std::tr1::shared_ptr<HelmholtzEOSMixtureBackend> ref_fluid(new HelmholtzEOSMixtureBackend(std::vector<std::string>(1, fluid_name)));
// Get the viscosity using ECS
return TransportRoutines::conductivity_ECS(*this, *(ref_fluid.get()));
}
if (component.transport.hardcoded_conductivity != CoolProp::TransportPropertyData::CONDUCTIVITY_NOT_HARDCODED)
{
switch(component.transport.hardcoded_conductivity)
{
case CoolProp::TransportPropertyData::CONDUCTIVITY_HARDCODED_WATER:
return TransportRoutines::conductivity_hardcoded_water(*this);
@@ -255,7 +289,7 @@ long double HelmholtzEOSMixtureBackend::calc_conductivity(void)
// Dilute part
long double lambda_dilute = _HUGE;
switch(components[0]->transport.conductivity_dilute.type)
switch(component.transport.conductivity_dilute.type)
{
case ConductivityDiluteVariables::CONDUCTIVITY_DILUTE_RATIO_POLYNOMIALS:
lambda_dilute = TransportRoutines::conductivity_dilute_ratio_polynomials(*this); break;
@@ -269,21 +303,11 @@ long double HelmholtzEOSMixtureBackend::calc_conductivity(void)
throw ValueError(format("dilute conductivity type [%d] is invalid for fluid %s", components[0]->transport.conductivity_dilute.type, name().c_str()));
}
// Residual part
long double lambda_residual = _HUGE;
switch(components[0]->transport.conductivity_residual.type)
{
case ConductivityResidualVariables::CONDUCTIVITY_RESIDUAL_POLYNOMIAL:
lambda_residual = TransportRoutines::conductivity_residual_polynomial(*this); break;
case ConductivityResidualVariables::CONDUCTIVITY_RESIDUAL_POLYNOMIAL_AND_EXPONENTIAL:
lambda_residual = TransportRoutines::conductivity_residual_polynomial_and_exponential(*this); break;
default:
throw ValueError(format("residual conductivity type [%d] is invalid for fluid %s", components[0]->transport.conductivity_residual.type, name().c_str()));
}
long double lambda_residual = calc_conductivity_background();
// Critical part
long double lambda_critical = _HUGE;
switch(components[0]->transport.conductivity_critical.type)
switch(component.transport.conductivity_critical.type)
{
case ConductivityCriticalVariables::CONDUCTIVITY_CRITICAL_SIMPLIFIED_OLCHOWY_SENGERS:
lambda_critical = TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers(*this); break;
@@ -1305,7 +1329,7 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double
if (!ValidNumber(rhomolar)){throw ValueError();}
return rhomolar;
}
catch(std::exception &e)
catch(std::exception &)
{
try{
// Next we try with Secant method
@@ -1313,7 +1337,7 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double
if (!ValidNumber(rhomolar)){throw ValueError();}
return rhomolar;
}
catch(std::exception &e)
catch(std::exception &)
{
Secant(resid, rhomolar_guess, 0.0001*rhomolar_guess, 1e-8, 100, errstring);
return _HUGE;
@@ -1538,6 +1562,19 @@ long double HelmholtzEOSMixtureBackend::calc_cpmolar(void)
return static_cast<double>(_cpmolar);
}
long double HelmholtzEOSMixtureBackend::calc_cpmolar_idealgas(void)
{
// Calculate the reducing parameters
_delta = _rhomolar/_reducing.rhomolar;
_tau = _reducing.T/_T;
// Calculate derivatives if needed, or just use cached values
long double d2a0_dTau2 = d2alpha0_dTau2();
long double R_u = static_cast<double>(_gas_constant);
// Get cp of the ideal gas
return R_u*-pow(_tau.pt(),2)*d2a0_dTau2;
}
long double HelmholtzEOSMixtureBackend::calc_speed_sound(void)
{
// Calculate the reducing parameters

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

@@ -93,6 +93,7 @@ public:
long double calc_pressure(void);
long double calc_cvmolar(void);
long double calc_cpmolar(void);
long double calc_cpmolar_idealgas(void);
long double calc_hmolar(void);
long double calc_smolar(void);
long double calc_pressure_nocache(long double T, long double rhomolar);
@@ -123,6 +124,7 @@ public:
long double calc_viscosity_background(void);
long double calc_viscosity_background(long double eta_dilute);
long double calc_conductivity(void);
long double calc_conductivity_background(void);
long double calc_Tmax(void);
long double calc_pmax(void);

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

@@ -542,7 +542,7 @@ long double TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers(
double Xref = Pcrit/pow(rhoc, 2)*HEOS.rhomolar()/dp_drho_ref*Tref/HEOS.T();
num = X - Xref;
// no critical enhancement if numerator is negative, zero, or just a tiny bit positive due to roundoff
// No critical enhancement if numerator is negative, zero, or just a tiny bit positive due to roundoff
// See also Lemmon, IJT, 2004, page 27
if (num < DBL_EPSILON*10)
return 0.0;
@@ -827,21 +827,18 @@ long double TransportRoutines::viscosity_ECS(HelmholtzEOSMixtureBackend &HEOS, H
rhocmolar = HEOS.rhomolar_critical(),
rhocmolar0 = HEOS_Reference.rhomolar_critical();
// Get a reference to the ECS data
CoolProp::ViscosityECSVariables &ECS = HEOS.components[0]->transport.viscosity_ecs;
std::vector<long double> &a = ECS.psi_a, &t = ECS.psi_t;
// The correction polynomial psi_eta
double psi = 0;
for (std::size_t i=0; i < a.size(); i++)
psi += a[i]*pow(HEOS.rhomolar()/ECS.psi_rhomolar_reducing, t[i]);
for (std::size_t i=0; i < ECS.psi_a.size(); i++){
psi += ECS.psi_a[i]*pow(HEOS.rhomolar()/ECS.psi_rhomolar_reducing, ECS.psi_t[i]);
}
// The dilute gas portion for the fluid of interest [Pa-s]
long double eta_dilute = viscosity_dilute_kinetic_theory(HEOS);
// Calculate the correction polynomial
/// \todo To be solved for...
// TODO: To be solved for...
long double theta = 1;
@@ -869,4 +866,76 @@ long double TransportRoutines::viscosity_ECS(HelmholtzEOSMixtureBackend &HEOS, H
return eta;
}
long double TransportRoutines::conductivity_ECS(HelmholtzEOSMixtureBackend &HEOS, HelmholtzEOSMixtureBackend &HEOS_Reference)
{
// Collect some parameters
long double M = HEOS.molar_mass(),
M_kmol = M*1000,
M0 = HEOS_Reference.molar_mass(),
Tc = HEOS.T_critical(),
Tc0 = HEOS_Reference.T_critical(),
rhocmolar = HEOS.rhomolar_critical(),
rhocmolar0 = HEOS_Reference.rhomolar_critical(),
R_u = HEOS.gas_constant(),
R = HEOS.gas_constant()/HEOS.molar_mass(), //[J/kg/K]
R_kJkgK = R_u/M_kmol;
// 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){
psi += ECS.psi_a[i]*pow(HEOS.rhomolar()/ECS.psi_rhomolar_reducing, ECS.psi_t[i]);
}
// The correction polynomial f_int in T/T_red
double fint = 0;
for (std::size_t i=0; i < ECS.f_int_a.size(); ++i){
fint += ECS.f_int_a[i]*pow(HEOS.T()/ECS.f_int_T_reducing, ECS.f_int_t[i]);
}
// The dilute gas density for the fluid of interest [uPa-s]
long double eta_dilute = viscosity_dilute_kinetic_theory(HEOS)*1e6;
// The mass specific ideal gas constant-pressure specific heat [J/kg/K]
long double cp0 = HEOS.calc_cpmolar_idealgas()/HEOS.molar_mass();
// The internal contribution to the thermal conductivity [W/m/K]
long double lambda_int = fint*eta_dilute*(cp0 - 2.5*R)/1e3;
// The dilute gas contribution to the thermal conductivity [W/m/K]
long double lambda_dilute = 15.0e-3/4.0*R_kJkgK*eta_dilute;
/// \todo To be solved for...
// TODO: To be solved for...
long double theta = 1;
long double phi = 1;
// The equivalent substance reducing ratios
long double f = Tc/Tc0*theta;
long double h = rhocmolar0/rhocmolar*phi; // Must be the ratio of MOLAR densities!!
// To be solved for
long double T0 = HEOS.T()/f;
long double rhomolar0 = HEOS.rhomolar()*h;
// Update the reference fluid with the conformal state
HEOS_Reference.update(DmolarT_INPUTS, rhomolar0*psi, T0);
// The reference fluid's contribution to the conductivity [W/m/K]
long double lambda_resid = HEOS_Reference.calc_conductivity_background();
// The F factor
long double F_lambda = sqrt(f)*pow(h, static_cast<long double>(-2.0/3.0))*sqrt(M0/M);
// The critical contribution from the fluid of interest [W/m/K]
long double lambda_critical = conductivity_critical_simplified_Olchowy_Sengers(HEOS);
// The total thermal conductivity considering the contributions of the fluid of interest and the reference fluid [Pa-s]
long double lambda = lambda_int + lambda_dilute + lambda_resid*F_lambda + lambda_critical;
return lambda;
}
}; /* namespace CoolProp */

2
src/Backends/Helmholtz/TransportRoutines.h
View File

@@ -193,6 +193,8 @@ public:
*/
static long double viscosity_ECS(HelmholtzEOSMixtureBackend &HEOS, HelmholtzEOSMixtureBackend &HEOS_Reference);
static long double conductivity_ECS(HelmholtzEOSMixtureBackend &HEOS, HelmholtzEOSMixtureBackend &HEOS_Reference);
}; /* class TransportRoutines */