github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-01-24 05:18:04 -05:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1 from jjfPCSI1/python_plots

Browse Source 
Python plots
This commit is contained in:
Jean-Julien Joseph Fleck
2014-08-27 00:13:31 +02:00
parent 6d3330bda3 058bfd79c4
commit c8656bcf12
20 changed files with 508 additions and 238 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
CMakeLists.txt
View File

@@ -159,8 +159,11 @@ if (COOLPROP_64BIT_SHARED_LIBRARY_MODULE OR COOLPROP_64BIT_SHARED_LIBRARY)
set_target_properties(${app_name} PROPERTIES COMPILE_FLAGS "-m64" LINK_FLAGS "-m64")
endif()
add_dependencies (${app_name} generate_headers)
set_target_properties(${app_name} PROPERTIES PREFIX "")
if (MSVC)
# No lib prefix for the shared library
set_target_properties(${app_name} PROPERTIES PREFIX "")
add_custom_command(TARGET ${app_name}
POST_BUILD
COMMAND dumpbin /EXPORTS $<TARGET_FILE:CoolProp> > ${CMAKE_CURRENT_BINARY_DIR}/exports.txt)
@@ -549,6 +552,9 @@ if (COOLPROP_MATLAB_MODULE)
add_library(HAPropsSI SHARED ${APP_SOURCES} ${CMAKE_SOURCE_DIR}/wrappers/MATLAB/Matlabdef.def ${CMAKE_SOURCE_DIR}/wrappers/MATLAB/HAProps.cpp)
target_link_libraries(PropsSI ${MATLAB_LIBRARIES})
target_link_libraries(HAPropsSI ${MATLAB_LIBRARIES})
set_target_properties(PropsSI PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) #Put .lib in this directory so it won't get installed
set_target_properties(HAPropsSI PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) #Put .lib in this directory so it won't get installed
if(WIN32) # 32-bit or 64-bit mex
if (CMAKE_CL_64)

12
CoolPropBibTeXLibrary.bib
View File

@@ -1645,6 +1645,18 @@
timestamp = {2013.04.09}
}
@ARTICLE{Xiang-JPCRD-2006,
author = {Hong Wei Xiang and Arno Laesecke and Marcia L. Huber},
title = {{A New Reference Correlation for the Viscosity of Methanol}},
journal = {J. Phys. Chem. Ref. Data},
year = {2006},
volume = {35},
pages = {1597-1:24},
number = {4},
owner = {Belli},
timestamp = {2014.08.26}
}
@ARTICLE{Younglove-JPCRD-1994,
author = {Ben A. Younglove},
title = {{An International Standard Equation of State for the Thermodynamic

4
TTD.txt
View File

@@ -127,8 +127,6 @@ R113
:::::::::::: QUEUE :::::::::::::
ECS params for new fluids
Methanol
visc: Xiang
Ethylene
visc & cond: Holland 1983
Methane
@@ -149,4 +147,4 @@ No Helmholtz:
C4F10 (ECS)
C5F12 (ECS)
R245CA (ECS)
"""
"""

2
Web/_templates/layout.html
View File

@@ -6,7 +6,7 @@
<li><a href="{{ pathto('online/index') }}">CoolProp Online</a>|</li>
<li><a href="{{ pathto('coolprop/examples/examples') }}">Examples</a>|</li>
<li><a href="{{ pathto('citation') }}">Citation</a>|</li>
<li><a href="{{ pathto('_static/doxygen/html/index') }}">Code</a>|</li>
<li><a href="{{ pathto('develop/code') }}">Code</a>|</li>
<li><a href="{{ pathto('coolprop/wrappers/wrappers') }}">Downloads</a></li>
{% endblock %}

9
Web/develop/code.rst Normal file
View File

@@ -0,0 +1,9 @@
Source Code
===========
* The source code of CoolProp is stored in a github repository at https://github.com/CoolProp/CoolProp
* To see nightly builds of the `doxygen <http://www.stack.nl/~dimitri/doxygen/>`_ formatted HTML outputs for the development code, go to http://www.coolprop.dreamhosters.com/doc/CoolProp5/
* Other documentation about the source code and developer notes can be found at :ref:`develop`.

10
Web/develop/documentation.rst
View File

@@ -7,21 +7,21 @@ Documentation
Build Sphinx documentation
--------------------------
1. Check out the sources in the CoolProp/Web folder
1. Check out the sources in the CoolProp/Web folder::
git clone https://github.com/CoolProp/CoolProp --recursive
2. Make a virtualenv - helps to keep sane !
2. Make a virtualenv::
virtualenv ~/env/py27
source ~/env/py27/bin/activate # Turn on this virtual env, should see (py27) in your command shell next to the prompt to tell you that environment is active
3. Then install prerequisites into this virtualenv:
3. Then install prerequisites into this virtualenv::
pip install Cython
pip install CoolProp sphinx sphinxcontrib-doxylink sphinxcontrib-napoleon cloud-sptheme ipython matplotlib numpy scipy
4. To build the documentation, go into the CoolProp/Web folder and run:
4. To build the documentation, go into the CoolProp/Web folder and run::
make html
@@ -34,7 +34,7 @@ All the configuration is done in the ``Doxyfile`` file.
1. If you don't have doxygen, install it from `doxygen downloads <http://www.stack.nl/~dimitri/doxygen/download.html>`. Or on linux a ``sudo apt-get install doxygen`` should do it.
2. Simply fire up a shell in the root of the repo, and type
2. Simply fire up a shell in the root of the repo, and type::
doxygen Doxyfile

2
Web/develop/index.rst
View File

@@ -1,3 +1,4 @@
.. _develop:
**************************
Information for Developers
@@ -5,6 +6,7 @@ Information for Developers
.. toctree::
code.rst
cmake.rst
buildbot.rst
documentation.rst

6
dev/fluids/Methanol.json
View File

@@ -455,5 +455,11 @@
"rhomolar": 0.00012793439653979163,
"rhomolar_units": "mol/m^3"
}
},
"TRANSPORT": {
"viscosity": {
"BibTeX": "Xiang-JPCRD-2006",
"hardcoded": "Methanol"
}
}
}

160
dev/fluids/R11.json
View File

@@ -42,108 +42,108 @@
"type": "rational_polynomial"
},
"pL": {
"Tmin": 162.68000000000001,
"T_r": 471.06,
"Tmax": 471.05,
"Tmin": 162.68,
"description": "p' = pc*exp(Tc/T*sum(n_i*theta^t_i))",
"Tmax": 471.05000000000001,
"using_tau_r": true,
"max_abserror_percentage": 0.035592288515440273,
"t": [
0.003,
0.029,
1.696,
2.279,
5.742,
8.448
"max_abserror_percentage": 0.025840508063401657,
"n": [
0.002488999538674866,
0.17252290950983326,
-6.293222945340061,
-0.07387863714475532,
-3.286105292292795,
2.274584153170513
],
"reducing_value": 4394000.0,
"T_r": 471.06,
"n": [
-7.879607134124751,
0.9890113036573663,
2.609469327659576,
-2.9284269362809985,
-4.154329007797212,
2.854121787641506
"t": [
0.002,
0.901,
0.973,
2.808,
4.144,
16.732
],
"type": "pL"
"type": "pL",
"using_tau_r": true
},
"pV": {
"Tmin": 162.68000000000001,
"T_r": 471.06,
"Tmax": 471.05,
"Tmin": 162.68,
"description": "p'' = pc*exp(Tc/T*sum(n_i*theta^t_i))",
"Tmax": 471.05000000000001,
"using_tau_r": true,
"max_abserror_percentage": 0.10227628818365586,
"t": [
0.014,
0.06,
0.22,
1.174,
3.682,
3.901
"max_abserror_percentage": 0.2290674443088614,
"n": [
-11.120479554886844,
281.2061736332651,
-650.2529791068788,
385.6318003272789,
-18.198375607518727,
4.101426759626743
],
"reducing_value": 4394000.0,
"T_r": 471.06,
"n": [
-10.566092482011701,
4.291071312457424,
-0.7131411992620907,
0.9063307505531512,
1.6297078047932856,
-4.932183821910964
"t": [
1.103,
2.059,
2.227,
2.373,
3.903,
6.255
],
"type": "pV"
"type": "pV",
"using_tau_r": true
},
"rhoL": {
"Tmin": 162.68000000000001,
"T_r": 471.06,
"Tmax": 471.05,
"Tmin": 162.68,
"description": "rho' = rhoc*(1+sum(n_i*theta^t_i))",
"Tmax": 471.05000000000001,
"using_tau_r": false,
"max_abserror_percentage": 0.16765507326250706,
"t": [
0.263,
0.55,
1.629,
2.573,
3.395,
19.679
"max_abserror_percentage": 0.032682517991688975,
"n": [
0.1369999075358185,
2.5657629003804248,
-0.5933415750742217,
0.7576689326270665,
-0.21901132344678031,
0.21836410566455441
],
"reducing_value": 4113.0394269407725,
"T_r": 471.06,
"n": [
0.761043427480207,
1.852550152931646,
-0.5766105549418219,
1.2213889632373067,
-0.5943875756199383,
8.044315782092438
"t": [
0.128,
0.465,
0.948,
1.743,
2.364,
7.362
],
"type": "rhoLnoexp"
"type": "rhoLnoexp",
"using_tau_r": false
},
"rhoV": {
"Tmin": 162.68000000000001,
"T_r": 471.06,
"Tmax": 471.05,
"Tmin": 162.68,
"description": "rho'' = rhoc*exp(Tc/T*sum(n_i*theta^t_i))",
"Tmax": 471.05000000000001,
"using_tau_r": true,
"max_abserror_percentage": 4.95412803941917,
"t": [
0.42,
2.574,
5.527,
9.81,
14.598,
16.056
"max_abserror_percentage": 0.1730324085996715,
"n": [
-0.19338479907103712,
20.681546570974632,
-24.441568955282317,
-1.761137580982086,
-9.939134168011334,
7.889769508867816
],
"reducing_value": 4113.0394269407725,
"T_r": 471.06,
"n": [
-3.584040196306734,
-9.671738689741934,
62.56987026226389,
-858.2036665098312,
14388.231325320496,
-18630.024609966502
"t": [
0.078,
0.408,
0.422,
1.519,
5.592,
7.237
],
"type": "rhoV"
"type": "rhoV",
"using_tau_r": true
},
"sL": {
"A": [

184
dev/fluids/R245fa.json
View File

@@ -43,109 +43,109 @@
"max_abs_error_units": "J/mol",
"type": "rational_polynomial"
},
"pL": {
"T_r": 427.16,
"Tmax": 427.15999999999934,
"Tmin": 171.05,
"description": "p' = pc*exp(Tc/T*sum(n_i*theta^t_i))",
"max_abserror_percentage": 0.04099399491688249,
"n": [
0.0004869832757370176,
-6.919145483936883,
-1.296447452384627,
-109.63729861885537,
107.66538794686012,
-3.4622351305818055
],
"reducing_value": 3651000.0,
"t": [
0.01,
0.975,
2.763,
4.376,
4.413,
6.994
],
"type": "pL",
"using_tau_r": true
},
"pV": {
"T_r": 427.16,
"Tmax": 427.15999999999934,
"Tmin": 171.05,
"Tmin": 171.05000000000001,
"description": "p'' = pc*exp(Tc/T*sum(n_i*theta^t_i))",
"max_abserror_percentage": 0.018077333401533835,
"n": [
-125.12984314778919,
152.9911248662765,
-35.27317598226143,
-4.102406386050718,
0.25362281678701887,
-1.994325281027591
"Tmax": 427.0,
"using_tau_r": true,
"max_abserror_percentage": 0.021575129298945228,
"t": [
0.775,
0.971,
2.155,
3.232,
9.097,
11.456
],
"reducing_value": 3651000.0,
"t": [
1.141,
1.179,
1.287,
3.64,
5.737,
7.077
],
"type": "pV",
"using_tau_r": true
},
"rhoL": {
"T_r": 427.16,
"Tmax": 427.15999999999934,
"Tmin": 171.05,
"description": "rho' = rhoc*(1+sum(n_i*theta^t_i))",
"max_abserror_percentage": 0.6726827141749103,
"T_r": 427.01,
"n": [
5.25176887594682,
-4.096294176783149,
10.999081251417799,
-11.640705049730522,
4.081973016432867,
-133.73894447755137
0.21292235127871662,
-7.22997740109537,
0.4023556144042939,
-4.366079807014841,
-9.155639983057442,
14.533193893397405
],
"reducing_value": 3849.9999999999995,
"t": [
0.42,
0.565,
1.624,
1.967,
4.012,
15.818
],
"type": "rhoLnoexp",
"using_tau_r": false
"type": "pV"
},
"rhoV": {
"T_r": 427.16,
"Tmax": 427.15999999999934,
"Tmin": 171.05,
"Tmin": 171.05000000000001,
"description": "rho'' = rhoc*exp(Tc/T*sum(n_i*theta^t_i))",
"max_abserror_percentage": 0.6375071010457756,
"n": [
-2.215307322989694,
-1.1625687686537725,
-7.315572965034845,
5.341825675363525,
-7.602957384362923,
17.028040285012285
],
"reducing_value": 3849.9999999999995,
"Tmax": 427.0,
"using_tau_r": true,
"max_abserror_percentage": 0.26467661268476661,
"t": [
0.383,
0.46,
1.442,
1.848,
3.855,
12.963
0.075,
0.466,
1.018,
2.817,
6.14,
8.5
],
"type": "rhoV",
"using_tau_r": true
"reducing_value": 3875.0,
"T_r": 427.01,
"n": [
-0.12316786769188867,
-3.5467162730203325,
-2.100626576983384,
-3.02259697521286,
-10.604906302496444,
12.900898794193413
],
"type": "rhoV"
},
"pL": {
"Tmin": 171.05000000000001,
"description": "p' = pc*exp(Tc/T*sum(n_i*theta^t_i))",
"Tmax": 427.0,
"using_tau_r": true,
"max_abserror_percentage": 0.013022195594358799,
"t": [
7.194,
1.03,
1.145,
2.84,
4.613,
18.798
],
"reducing_value": 3651000.0,
"T_r": 427.01,
"n": [
-0.9670903356725931,
-11.756404923945222,
5.087651519671258,
-2.835805078027625,
-2.217939014066299,
28.8905424182831
],
"type": "pL"
},
"rhoL": {
"Tmin": 171.05000000000001,
"description": "rho' = rhoc*(1+sum(n_i*theta^t_i))",
"Tmax": 427.0,
"using_tau_r": false,
"max_abserror_percentage": 0.040900248845132658,
"t": [
0.144,
0.165,
0.53,
1.993,
2.646,
3.593
],
"reducing_value": 3875.0,
"T_r": 427.01,
"n": [
-0.645916749226158,
1.2435143264412414,
2.1374254146018647,
-1.4105539340976712,
2.697825229930464,
-1.236775453941409
],
"type": "rhoLnoexp"
},
"sL": {
"A": [

22
include/CoolPropFluid.h
View File

@@ -46,29 +46,30 @@ struct CriticalRegionSplines{
bool enabled;
CriticalRegionSplines(){enabled = false;};
void get_densities(double T, double rho_min, double rho_crit, double rho_max, double &rhoL, double &rhoV){
bool ok1 = false, ok2 = false, ok3 = false;
int Nsoln = -1, Ngood = 0;
double rho1, rho2, rho3;
double rho1 =0, rho2 = 0, rho3 = 0;
// -----------
// Liquid part
// -----------
Ngood = 0;
solve_cubic(cL[0], cL[1], cL[2], cL[3]-T, Nsoln, rho1, rho2, rho3);
if (rho1 < rho_max && rho1 > rho_crit){ ok1 = true; Ngood++; rhoL = rho1; }
if (rho2 < rho_max && rho2 > rho_crit){ ok2 = true; Ngood++; rhoL = rho2; }
if (rho3 < rho_max && rho3 > rho_crit){ ok3 = true; Ngood++; rhoL = rho3; }
if (Ngood != 1){ throw ValueError(format("more than one liquid solution found for critical spline for T=%g",T));};
if (rho1 < rho_max && rho1 > rho_crit){ Ngood++; rhoL = rho1; }
if (rho2 < rho_max && rho2 > rho_crit){ Ngood++; rhoL = rho2; }
if (rho3 < rho_max && rho3 > rho_crit){ Ngood++; rhoL = rho3; }
if (Ngood > 1){ throw ValueError(format("More than one liquid solution found for critical spline for T=%g",T));};
if (Ngood < 1){ throw ValueError(format("No liquid solution found for critical spline for T=%g",T));};
// ----------
// Vapor part
// ----------
Ngood = 0;
solve_cubic(cV[0], cV[1], cV[2], cV[3]-T, Nsoln, rho1, rho2, rho3);
if (rho1 > rho_min && rho1 < rho_crit){ ok1 = true; Ngood++; rhoV = rho1; }
if (rho2 > rho_min && rho2 < rho_crit){ ok2 = true; Ngood++; rhoV = rho2; }
if (rho3 > rho_min && rho3 < rho_crit){ ok3 = true; Ngood++; rhoV = rho3; }
if (Ngood != 1){ throw ValueError(format("more than one vapor solution found for critical spline for T=%g",T));};
if (rho1 > rho_min && rho1 < rho_crit){ Ngood++; rhoV = rho1; }
if (rho2 > rho_min && rho2 < rho_crit){ Ngood++; rhoV = rho2; }
if (rho3 > rho_min && rho3 < rho_crit){ Ngood++; rhoV = rho3; }
if (Ngood > 1){ throw ValueError(format("More than one vapor solution found for critical spline for T=%g",T));};
if (Ngood < 1){ throw ValueError(format("No vapor solution found for critical spline for T=%g",T));};
};
};
@@ -259,6 +260,7 @@ public:
enum ViscosityHardcodedEnum {VISCOSITY_HARDCODED_WATER, ///< Use \ref TransportRoutines::viscosity_water_hardcoded
VISCOSITY_HARDCODED_HELIUM, ///< Use \ref TransportRoutines::viscosity_helium_hardcoded
VISCOSITY_HARDCODED_R23, ///< Use \ref TransportRoutines::viscosity_R23_hardcoded
VISCOSITY_HARDCODED_METHANOL, ///< Use \ref TransportRoutines::viscosity_methanol_hardcoded
VISCOSITY_NOT_HARDCODED
};
enum ConductivityHardcodedEnum {

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

@@ -640,6 +640,9 @@ protected:
else if (!target.compare("R23")){
fluid.transport.hardcoded_viscosity = CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_R23; return;
}
else if (!target.compare("Methanol")){
fluid.transport.hardcoded_viscosity = CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_METHANOL; return;
}
else{
throw ValueError(format("hardcoded viscosity [%s] is not understood for fluid %s",target.c_str(), fluid.name.c_str()));
}

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

@@ -356,6 +356,8 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity(void)
return TransportRoutines::viscosity_helium_hardcoded(*this);
case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_R23:
return TransportRoutines::viscosity_R23_hardcoded(*this);
case CoolProp::TransportPropertyData::VISCOSITY_HARDCODED_METHANOL:
return TransportRoutines::viscosity_methanol_hardcoded(*this);
default:
throw ValueError(format("hardcoded viscosity type [%d] is invalid for fluid %s", component.transport.hardcoded_viscosity, name().c_str()));
}

134
src/Backends/Helmholtz/MixtureDerivatives.cpp
View File

@@ -210,7 +210,6 @@ long double MixtureDerivatives::d_ndalphardni_dDelta(HelmholtzEOSMixtureBackend
}
return term1 + term2 + term3;
}
long double MixtureDerivatives::d_ndalphardni_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i)
{
// The first line
@@ -228,5 +227,136 @@ long double MixtureDerivatives::d_ndalphardni_dTau(HelmholtzEOSMixtureBackend &H
return term1 + term2 + term3;
}
} /* namespace CoolProp */
#ifdef ENABLE_CATCH
#include "catch.hpp"
using namespace CoolProp;
TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
{
/* Set up a test class for the mixture tests */
std::vector<std::string> names(2);
names[0] = "Methane"; names[1] = "Ethane";
std::vector<long double> z(2);
z[0] = 0.5; z[1] = 1-z[0];
shared_ptr<HelmholtzEOSMixtureBackend> HEOS(new HelmholtzEOSMixtureBackend(names));
HelmholtzEOSMixtureBackend &rHEOS = *(HEOS.get());
rHEOS.set_mole_fractions(z);
// These ones only require the i index
for (std::size_t i = 0; i< z.size();++i)
{
std::ostringstream ss1;
ss1 << "dln_fugacity_coefficient_dT__constp_n, i=" << i;
SECTION(ss1.str(), "")
{
double T1 = 300, dT = 1e-3;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(PT_INPUTS, 101325, T1);
double analytic = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rHEOS, i);
rHEOS.update(PT_INPUTS, 101325, T1 + dT);
double v1 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i);
rHEOS.update(PT_INPUTS, 101325, T1 - dT);
double v2 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i);
double numeric = (v1 - v2)/(2*dT);
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-6);
}
std::ostringstream ss2;
ss2 << "dln_fugacity_coefficient_dp__constT_n, i=" << i;
SECTION(ss2.str(), "")
{
double p1 = 101325, dP = 1;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(PT_INPUTS, p1, 300);
double analytic = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(rHEOS, i);
rHEOS.update(PT_INPUTS, p1 + dP, 300);
double v1 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i);
rHEOS.update(PT_INPUTS, p1 - dP, 300);
double v2 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i);
double numeric = (v1 - v2)/(2*dP);
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-6);
}
std::ostringstream ss3;
ss3 << "d_ndalphardni_dDelta, i=" << i;
SECTION(ss3.str(), "")
{
double p1 = 101325, dP = 1e-1;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(PT_INPUTS, p1, 300);
double analytic = MixtureDerivatives::d_ndalphardni_dDelta(rHEOS, i);
rHEOS.update(PT_INPUTS, p1 + dP, 300);
double v1 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i), delta1 = rHEOS.delta();
rHEOS.update(PT_INPUTS, p1 - dP, 300);
double v2 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i), delta2 = rHEOS.delta();
double numeric = (v1 - v2)/(delta1 - delta2);
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-6);
}
std::ostringstream ss4;
ss4 << "d_ndalphardni_dTau, i=" << i;
SECTION(ss4.str(), "")
{
double p1 = 101325, dT = 1e-2;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(PT_INPUTS, 101325, 300);
double rho1 = rHEOS.rhomolar();
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double analytic = MixtureDerivatives::d_ndalphardni_dTau(rHEOS, i);
rHEOS.update(DmolarT_INPUTS, rho1, 300 + dT);
double v1 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i), tau1 = rHEOS.tau();
rHEOS.update(DmolarT_INPUTS, rho1, 300 - dT);
double v2 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i), tau2 = rHEOS.tau();
double numeric = (v1 - v2)/(tau1 - tau2);
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-6);
}
// These derivatives depend on both the i and j indices
for (std::size_t j = 0; j < z.size(); ++j){
std::ostringstream ss1;
ss1 << "dln_fugacity_coefficient_dxj__constT_p_xi, i=" << i << ", j=" << j;
SECTION(ss1.str(), "")
{
double dz = 1e-6;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(PT_INPUTS, 101325, 300);
double rho1 = rHEOS.rhomolar();
double analytic = MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rHEOS, i, j);
std::vector<long double> zp = z; /// Copy base composition
zp[j] += dz;
rHEOS.set_mole_fractions(zp);
rHEOS.update(PT_INPUTS, 101325, 300);
double v1 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i);
std::vector<long double> zm = z; /// Copy base composition
zm[j] -= dz;
rHEOS.set_mole_fractions(zm);
rHEOS.update(PT_INPUTS, 101325, 300);
double v2 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i);
double numeric = (v1 - v2)/(2*dz);
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-6);
}
}
}
}
#endif

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

@@ -372,8 +372,6 @@ long double TransportRoutines::viscosity_higher_order_friction_theory(HelmholtzE
else{
throw NotImplementedError("TransportRoutines::viscosity_higher_order_friction_theory is only for pure and pseudo-pure");
}
}
long double TransportRoutines::viscosity_helium_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
@@ -418,6 +416,67 @@ long double TransportRoutines::viscosity_helium_hardcoded(HelmholtzEOSMixtureBac
}
}
long double TransportRoutines::viscosity_methanol_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
{
long double B_eta, C_eta,
epsilon_over_k = 577.87, /* [K]*/
sigma0 = 0.3408e-9, /* [m] */
delta = 0.4575, /* NOT the reduced density, that is rhor here*/
N_A = 6.02214129e23,
M = 32.04216, /* kg/kmol */
T = HEOS.T();
long double rhomolar = HEOS.rhomolar();
long double B_eta_star, C_eta_star;
long double Tstar = T/epsilon_over_k; // [no units]
long double rhor = HEOS.rhomass()/273;
long double Tr = T/512.6;
// Rainwater-Friend initial density terms
{ // Scoped here so that we can re-use the b variable
long double b[9] = {-19.572881, 219.73999, -1015.3226, 2471.01251, -3375.1717, 2491.6597, -787.26086, 14.085455, -0.34664158};
long double t[9] = {0, -0.25, -0.5, -0.75, -1.0, -1.25, -1.5, -2.5, -5.5};
long double summer = 0;
for (unsigned int i = 0; i < 9; ++i){
summer += b[i]*pow(Tstar, t[i]);
}
B_eta_star = summer; // [no units]
B_eta = N_A*pow(sigma0, 3)*B_eta_star; // [m^3/mol]
long double c[2] = {1.86222085e-3, 9.990338};
C_eta_star = c[0]*pow(Tstar,3)*exp(c[1]*pow(Tstar,-0.5)); // [no units]
C_eta = pow(N_A*pow(sigma0, 3), 2)*C_eta_star; // [m^6/mol^2]
}
long double eta_g = 1 + B_eta*rhomolar + C_eta*rhomolar*rhomolar;
long double a[13] = {1.16145, -0.14874, 0.52487, -0.77320, 2.16178, -2.43787, 0.95976e-3, 0.10225, -0.97346, 0.10657, -0.34528, -0.44557, -2.58055};
long double d[7] = {-1.181909, 0.5031030, -0.6268461, 0.5169312, -0.2351349, 5.3980235e-2, -4.9069617e-3};
long double e[10] = {0, 4.018368, -4.239180, 2.245110, -0.5750698, 2.3021026e-2, 2.5696775e-2, -6.8372749e-3, 7.2707189e-4, -2.9255711e-5};
long double OMEGA_22_star_LJ = a[0]*pow(Tstar,a[1])+a[2]*exp(a[3]*Tstar)+a[4]*exp(a[5]*Tstar);
long double OMEGA_22_star_delta = a[7]*pow(Tstar,a[8]) + a[9]*exp(a[10]*Tstar) + a[11]*exp(a[12]*Tstar);
long double OMEGA_22_star_SM = OMEGA_22_star_LJ*(1+pow(delta,2)/(1+a[6]*pow(delta,6))*OMEGA_22_star_delta);
long double eta_0 = 2.66957e-26*sqrt(M*T)/(pow(sigma0,2)*OMEGA_22_star_SM);
long double summerd = 0;
for (unsigned int i = 0; i < 7; ++i){
summerd += d[i]/pow(Tr, i);
}
for (unsigned int j = 1; j < 10; ++j){
summerd += e[j]*pow(rhor, j);
}
long double sigmac = 0.7193422e-9; // [m]
long double sigma_HS = summerd*sigmac; // [m]
long double b = 2*M_PI*N_A*pow(sigma_HS,3)/3; // [m^3/mol]
long double zeta = b*rhomolar/4; // [-]
long double g_sigma_HS = (1 - 0.5*zeta)/pow(1 - zeta, 3); // [-]
long double eta_E = 1/g_sigma_HS + 0.8*b*rhomolar + 0.761*g_sigma_HS*pow(b*rhomolar,2); // [-]
long double f = 1/(1+exp(5*(rhor-1)));
return eta_0*(f*eta_g + (1-f)*eta_E);
}
long double TransportRoutines::viscosity_R23_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
{
double C1 = 1.3163, //

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

@@ -96,6 +96,12 @@ public:
static long double viscosity_dilute_ethane(HelmholtzEOSMixtureBackend &HEOS);
static long double viscosity_dilute_cyclohexane(HelmholtzEOSMixtureBackend &HEOS);
/** \brief Viscosity hardcoded for Methanol
*
* From Xiang et al., A New Reference Correlation for the Viscosity of Methanol, J. Phys. Chem. Ref. Data, Vol. 35, No. 4, 2006
*/
static long double viscosity_methanol_hardcoded(HelmholtzEOSMixtureBackend &HEOS);
static long double viscosity_water_hardcoded(HelmholtzEOSMixtureBackend &HEOS);
static long double viscosity_helium_hardcoded(HelmholtzEOSMixtureBackend &HEOS);
static long double viscosity_R23_hardcoded(HelmholtzEOSMixtureBackend &HEOS);

38
src/CoolPropTools.cpp
View File

@@ -161,9 +161,47 @@ std::string get_file_contents(const char *filename)
void solve_cubic(double a, double b, double c, double d, int &N, double &x0, double &x1, double &x2)
{
// 0 = ax^3 + b*x^2 + c*x + d
if (std::abs(a) < 10*DBL_EPSILON){
if (std::abs(b) < 10*DBL_EPSILON){
// Linear solution if a = 0 and b = 0
x0 = -d/c;
N = 1;
return;
}
else{
// Quadratic solution(s) if a = 0 and b != 0
x0 = (-c+sqrt(c*c-4*b*d))/(2*b);
x1 = (-c-sqrt(c*c-4*b*d))/(2*b);
N = 2;
return;
}
}
// Ok, it is really a cubic
// Discriminant
double DELTA = 18*a*b*c*d-4*b*b*b*d+b*b*c*c-4*a*c*c*c-27*a*a*d*d;
double DELTA0 = b*b*b-2*a*c;
// Deal with special cases
if (std::abs(DELTA) < 10*DBL_EPSILON){
if (std::abs(DELTA0)>0){
x0 = (9*a*d-b*c)/(2*DELTA0);
x1 = x0;
x2 = (4*a*b*c - 9*a*a*d - b*b*b)/(a*DELTA0);
N = 3;
return;
}
else{
x0 = -b/(3*a);
x1 = x0;
x2 = x0;
N = 3;
return;
}
}
// Coefficients for the depressed cubic t^3+p*t+q = 0
double p = (3*a*c-b*b)/(3*a*a);

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

@@ -208,19 +208,10 @@ vel("Ethane", "T", 430, "Dmolar", 12780, "V", 58.70e-6, 1e-2),
vel("Ethane", "T", 500, "Dmolar", 11210, "V", 48.34e-6, 1e-2),
// From Xiang, JPCRD, 2006
vel("Methanol", "T", 600, "Dmass", 800.23, "V", 0.1888e-3, 1e-4),
vel("Methanol", "T", 600, "Dmass", 833.20, "V", 0.2092e-3, 1e-4),
vel("Methanol", "T", 600, "Dmass", 861.37, "V", 0.2279e-3, 1e-4),
vel("Methanol", "T", 600, "Dmass", 908.33, "V", 0.2634e-3, 1e-4),
vel("Methanol", "T", 620, "Dmass", 788.58, "V", 0.1779e-3, 1e-4),
vel("Methanol", "T", 620, "Dmass", 822.14, "V", 0.1972e-3, 1e-4),
vel("Methanol", "T", 620, "Dmass", 850.77, "V", 0.2148e-3, 1e-4),
vel("Methanol", "T", 620, "Dmass", 898.48, "V", 0.2477e-3, 1e-4),
vel("Methanol", "T", 630, "Dmass", 782.76, "V", 0.1729e-3, 1e-4),
vel("Methanol", "T", 630, "Dmass", 811.06, "V", 0.1917e-3, 1e-4),
vel("Methanol", "T", 630, "Dmass", 840.11, "V", 0.2088e-3, 1e-4),
vel("Methanol", "T", 630, "Dmass", 888.50, "V", 0.2405e-3, 1e-4),
vel("Methanol", "T", 300, "Dmass", 0.12955, "V", 0.009696e-3, 1e-3),
vel("Methanol", "T", 300, "Dmass", 788.41, "V", 0.5422e-3, 1e-3),
vel("Methanol", "T", 630, "Dmass", 0.061183, "V", 0.02081e-3, 1e-3),
vel("Methanol", "T", 630, "Dmass", 888.50, "V", 0.2405e-3, 1e-1), // They use a different EOS in the high pressure region
// From REFPROP 9.1 since no data provided
vel("n-Butane", "T", 150, "Q", 0, "V", 0.0013697657668, 1e-4),
@@ -275,7 +266,7 @@ TEST_CASE_METHOD(TransportValidationFixture, "Compare viscosities against publis
{
vel el = viscosity_validation_data[i];
CHECK_NOTHROW(set_backend("HEOS", el.fluid));
CAPTURE(el.fluid);
CAPTURE(el.in1);
CAPTURE(el.v1);
@@ -1124,7 +1115,6 @@ TEST_CASE("Ancillary functions", "[ancillary]")
{
double T = f*Tc + (1-f)*Tt;
std::ostringstream ss1;
ss1 << "Pressure error < 2% for fluid " << fluids[i] << " at " << T << " K";
SECTION(ss1.str(), "")
@@ -1186,17 +1176,11 @@ TEST_CASE("Triple point checks", "[triple_point]")
double p_sat_min_vapor = HEOS->get_components()[0]->pEOS->sat_min_vapor.p;
double err_sat_min_liquid = std::abs(p_EOS-p_sat_min_liquid)/p_sat_min_liquid;
double err_sat_min_vapor = std::abs(p_EOS-p_sat_min_vapor)/p_sat_min_vapor;
double p_triple_liquid = HEOS->get_components()[0]->triple_liquid.p;
double p_triple_vapor = HEOS->get_components()[0]->triple_liquid.p;
double err_triple_liquid = std::abs(p_EOS-p_triple_liquid)/p_triple_liquid;
double err_triple_vapor = std::abs(p_EOS-p_triple_vapor)/p_triple_vapor;
CAPTURE(err_sat_min_liquid);
CAPTURE(err_sat_min_vapor);
CHECK(err_sat_min_liquid < 1e-3);
CHECK(err_sat_min_vapor < 1e-3);
CHECK(err_triple_liquid < 1e-3);
CHECK(err_triple_vapor < 1e-3);
}
// std::ostringstream ss2;

28
wrappers/Python/CoolProp/Plots/Common.py
View File

@@ -11,12 +11,20 @@ SMALL = 1E-5
class BasePlot(object):
#TODO: Simplify / Consolidate dictionary maps
AXIS_LABELS = {'T': ["Temperature", r"[K]"],
'P': ["Pressure", r"[kPa]"],
'S': ["Entropy", r"[kJ/kg/K]"],
'H': ["Enthalpy", r"[kJ/kg]"],
'U': ["Internal Energy", r"[kJ/kg]"],
'D': ["Density", r"[kg/m$^3$]"]
AXIS_LABELS = {'KSI': {'T': ["Temperature", r"[K]"],
'P': ["Pressure", r"[kPa]"],
'S': ["Entropy", r"[kJ/kg/K]"],
'H': ["Enthalpy", r"[kJ/kg]"],
'U': ["Internal Energy", r"[kJ/kg]"],
'D': ["Density", r"[kg/m$^3$]"]
},
'SI': {'T': ["Temperature", r"[K]"],
'P': ["Pressure", r"[Pa]"],
'S': ["Entropy", r"[J/kg/K]"],
'H': ["Enthalpy", r"[J/kg]"],
'U': ["Internal Energy", r"[J/kg]"],
'D': ["Density", r"[kg/m$^3$]"]
}
}
COLOR_MAP = {'T': 'Darkred',
@@ -236,13 +244,13 @@ class BasePlot(object):
tl_str = "%s - %s Graph for %s"
if not self.axis.get_title():
self.axis.set_title(tl_str % (self.AXIS_LABELS[y_axis_id][0],
self.AXIS_LABELS[x_axis_id][0],
self.axis.set_title(tl_str % (self.AXIS_LABELS[self.unit_system][y_axis_id][0],
self.AXIS_LABELS[self.unit_system][x_axis_id][0],
filter_fluid_ref(self.fluid_ref)))
if not self.axis.get_xlabel():
self.axis.set_xlabel(' '.join(self.AXIS_LABELS[x_axis_id]))
self.axis.set_xlabel(' '.join(self.AXIS_LABELS[self.unit_system][x_axis_id]))
if not self.axis.get_ylabel():
self.axis.set_ylabel(' '.join(self.AXIS_LABELS[y_axis_id]))
self.axis.set_ylabel(' '.join(self.AXIS_LABELS[self.unit_system][y_axis_id]))
def _draw_graph(self):
return

27
wrappers/Python/CoolProp/Plots/Plots.py
View File

@@ -183,8 +183,8 @@ def drawLines(Ref,lines,axis,plt_kwargs=None):
class IsoLines(BasePlot):
def __init__(self, fluid_ref, graph_type, iso_type, **kwargs):
BasePlot.__init__(self, fluid_ref, graph_type, **kwargs)
def __init__(self, fluid_ref, graph_type, iso_type, unit_system='SI', **kwargs):
BasePlot.__init__(self, fluid_ref, graph_type, unit_system=unit_system,**kwargs)
if not isinstance(iso_type, str):
raise TypeError("Invalid iso_type input, expected a string")
@@ -247,6 +247,11 @@ class IsoLines(BasePlot):
limits[0][1] = min([limits[0][1], max(self.axis.get_xlim())])
limits[1][0] = max([limits[1][0], min(self.axis.get_ylim())])
limits[1][1] = min([limits[1][1], max(self.axis.get_ylim())])
# Limits correction in case of KSI unit_system
if self.unit_system == 'KSI':
limits[0] = [l*self.KSI_SCALE_FACTOR[self.graph_type[1]] for l in limits[0]]
limits[1] = [l*self.KSI_SCALE_FACTOR[self.graph_type[0]] for l in limits[1]]
self.axis.set_xlim(limits[0])
self.axis.set_ylim(limits[1])
@@ -465,7 +470,7 @@ class PropsPlot(BasePlot):
See the online documentation for a list of the available fluids and
graph types
"""
BasePlot.__init__(self, fluid_name, graph_type, **kwargs)
BasePlot.__init__(self, fluid_name, graph_type, unit_system=units, **kwargs)
self.smin = kwargs.get('smin', None)
self.smax = kwargs.get('smax', None)
@@ -485,7 +490,7 @@ class PropsPlot(BasePlot):
def draw_isolines(self, iso_type, iso_range, num=10, rounding=False):
iso_lines = IsoLines(self.fluid_ref,
self.graph_type,
iso_type,
iso_type, unit_system = self.unit_system,
axis=self.axis)
iso_lines.draw_isolines(iso_range, num, rounding)
@@ -526,7 +531,7 @@ def Ts(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
>>> ax = fig.gca()
>>> Ts('R290', show=True, axis=ax)
"""
plt = PropsPlot(Ref, 'Ts', smin=Tmin, smax=Tmax, axis=axis)
plt = PropsPlot(Ref, 'Ts', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
plt._draw_graph()
if show:
plt.show()
@@ -571,7 +576,7 @@ def Ph(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
>>> ax = fig.gca()
>>> Ph('R290', show=True, axis=ax)
"""
plt = PropsPlot(Ref, 'Ph', smin=Tmin, smax=Tmax, axis=axis)
plt = PropsPlot(Ref, 'Ph', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
@@ -615,7 +620,7 @@ def Ps(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
>>> ax = fig.gca()
>>> Ps('R290', show=True, axis=ax)
"""
plt = PropsPlot(Ref, 'Ps', smin=Tmin, smax=Tmax, axis=axis)
plt = PropsPlot(Ref, 'Ps', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
@@ -658,7 +663,7 @@ def PT(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
>>> ax = fig.gca()
>>> PT('R290', show=True, axis=ax)
"""
plt = PropsPlot(Ref, 'PT', smin=Tmin, smax=Tmax, axis=axis)
plt = PropsPlot(Ref, 'PT', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
@@ -701,7 +706,7 @@ def Prho(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
>>> ax = fig.gca()
>>> Prho('R290', show=True, axis=ax)
"""
plt = PropsPlot(Ref, 'PD', smin=Tmin, smax=Tmax, axis=axis)
plt = PropsPlot(Ref, 'PD', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
@@ -744,7 +749,7 @@ def Trho(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
>>> ax = fig.gca()
>>> Trho('R290', show=True, axis=ax)
"""
plt = PropsPlot(Ref, 'TD', smin=Tmin, smax=Tmax, axis=axis)
plt = PropsPlot(Ref, 'TD', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
@@ -787,7 +792,7 @@ def hs(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
>>> ax = fig.gca()
>>> hs('R290', show=True, axis=ax)
"""
plt = PropsPlot(Ref, 'hs', smin=Tmin, smax=Tmax, axis=axis)
plt = PropsPlot(Ref, 'hs', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else: