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

Merge branch 'eigenPolynomials'

Browse Source
This commit is contained in:
jowr
2014-06-29 19:15:58 +02:00
parent 134297644a 3b02f98879
commit 5778067a39
19 changed files with 5999 additions and 2534 deletions
Download Patch File Download Diff File Expand all files Collapse all files

306
dev/IncompressibleLiquids/LinearAlgebra.ipynb Normal file
View File

File diff suppressed because one or more lines are too long

57
dev/IncompressibleLiquids/PureExample.json Normal file
View File

@@ -0,0 +1,57 @@
{
"T_freeze": {
"coeffs": "null",
"type": 0
},
"Tbase": 0.0,
"Tmax": 533.15,
"Tmin": 188.14999999999998,
"TminPsat": 313.15,
"conductivity": {
"coeffs": "null",
"type": 0
},
"density": {
"coeffs": [
[
-5.439005439004754e-06
],
[
0.005552564102563395
],
[
-2.6422007556329925
],
[
1197.8367716325445
]
],
"type": "polynomial"
},
"description": "Heat transfer fluid TherminolD12 by Solutia",
"mass2mole": {
"coeffs": "null",
"type": 0
},
"name": "PureExample",
"reference": "Solutia data sheet",
"saturation_pressure": {
"coeffs": "null",
"type": 0
},
"specific_heat": {
"coeffs": "null",
"type": 0
},
"viscosity": {
"coeffs": "null",
"type": 0
},
"volume2mass": {
"coeffs": "null",
"type": 0
},
"xbase": 0.0,
"xmax": null,
"xmin": null
}

182
dev/IncompressibleLiquids/SecCoolExample.json Normal file
View File

@@ -0,0 +1,182 @@
{
"T_freeze": {
"coeffs": [
27.7555556,
-22.9732217,
-1.10405072,
-0.0120762281,
-9.343458e-05
],
"type": "polynomial"
},
"Tbase": 268.66999999999996,
"Tmax": 293.15,
"Tmin": 223.14999999999998,
"TminPsat": 293.15,
"conductivity": {
"coeffs": [
[
0.40820667,
-0.003981687,
1.583368e-05,
-3.552049e-07,
-9.884176e-10,
4.46e-10
],
[
0.0006629321,
-2.686475e-05,
9.03915e-07,
-2.128257e-08,
-5.562e-10,
0.0
],
[
3.685975e-07,
7.188416e-08,
-1.041773e-08,
2.278001e-10,
0.0,
0.0
],
[
4.703395e-08,
7.612361e-11,
-2.734e-10,
0.0,
0.0,
0.0
]
],
"type": "polynomial"
},
"density": {
"coeffs": [
[
960.246658,
1.29038391,
-0.016104252,
-0.0001969888,
1.131559e-05,
9.181999e-08
],
[
-0.402034827,
-0.0162463989,
0.0001623301,
4.367343e-06,
1.199e-08,
0.0
],
[
-0.0025204776,
0.0001101514,
-2.320217e-07,
7.794999e-08,
0.0,
0.0
],
[
9.937483e-06,
-1.346886e-06,
4.141999e-08,
0.0,
0.0,
0.0
]
],
"type": "polynomial"
},
"description": "Methanol solution",
"mass2mole": {
"coeffs": "null",
"type": 0
},
"name": "SecCoolExample",
"reference": "SecCool software",
"saturation_pressure": {
"coeffs": "null",
"type": 0
},
"specific_heat": {
"coeffs": [
[
3822.97123,
-23.1224095,
0.0678775826,
0.0022413893,
-0.0003045332,
-4.758e-06
],
[
2.35014495,
0.178883941,
0.0006828,
0.0002101166,
-9.812e-06,
0.0
],
[
-0.0004724176,
-0.0003317949,
0.0001002032,
-5.306e-06,
0.0,
0.0
],
[
4.242194e-05,
2.34719e-05,
-1.894e-06,
0.0,
0.0,
0.0
]
],
"type": "polynomial"
},
"viscosity": {
"coeffs": [
[
1.47255255,
0.0022218998,
-0.0004406139,
6.047984e-06,
-1.95473e-07,
-2.372e-09
],
[
-0.0411841566,
0.0001784479,
-3.564413e-06,
4.064671e-08,
1.915e-08,
0.0
],
[
0.0002572862,
-9.226343e-07,
-2.178577e-08,
-9.529999e-10,
0.0,
0.0
],
[
-1.699844e-06,
-1.023552e-07,
4.482e-09,
0.0,
0.0,
0.0
]
],
"type": "polynomial"
},
"volume2mass": {
"coeffs": "null",
"type": 0
},
"xbase": 0.3157,
"xmax": 0.5,
"xmin": 0.0
}

63
dev/IncompressibleLiquids/SolutionExample.json Normal file
View File

@@ -0,0 +1,63 @@
{
"T_freeze": {
"coeffs": "null",
"type": 0
},
"Tbase": 0.0,
"Tmax": null,
"Tmin": null,
"TminPsat": null,
"conductivity": {
"coeffs": "null",
"type": 0
},
"density": {
"coeffs": [
[
-1.1876190476190525,
-1.0978571428571455,
-1.008571428571432,
-0.926190476190484,
-0.8433333333333355,
-0.7626190476190485,
-0.6845238095238149
],
[
1338.9886190476204,
1308.851107142858,
1278.9555714285727,
1250.893690476193,
1222.8398333333341,
1195.3998690476192,
1168.7407738095253
]
],
"type": "polynomial"
},
"description": "Ethanol ice slurry",
"mass2mole": {
"coeffs": "null",
"type": 0
},
"name": "SolutionExample",
"reference": "SecCool software",
"saturation_pressure": {
"coeffs": "null",
"type": 0
},
"specific_heat": {
"coeffs": "null",
"type": 0
},
"viscosity": {
"coeffs": "null",
"type": 0
},
"volume2mass": {
"coeffs": "null",
"type": 0
},
"xbase": 0.0,
"xmax": null,
"xmin": null
}

418
dev/IncompressibleLiquids/data_incompressible.py
View File

@@ -1,4 +1,4 @@
import numpy
import numpy,itertools,scipy.interpolate
class LiquidData(object):
@@ -532,29 +532,405 @@ class ZS55(LiquidData):
#dataObj = Therminol72()
#for i in range(len(dataObj.T)):
# print str(dataObj.T[i])+", "+str(numpy.log(dataObj.mu_dyn[i]))
class IncompressibleData(object):
def __init__(self):
self.INCOMPRESSIBLE_NOT_SET = 0
self.INCOMPRESSIBLE_POLYNOMIAL = 'polynomial'
self.INCOMPRESSIBLE_EXPONENTIAL = 'exponential'
self.INCOMPRESSIBLE_EXPPOLYNOMIAL = 'exppolynomial'
self.type = self.INCOMPRESSIBLE_NOT_SET
self.coeffs = None #numpy.zeros((4,4))
self.data = None #numpy.zeros((10,10))
def fit(self, T, x=0.0, Tbase=0.0, xbase=0.0):
(cr,cc) = self.coeffs.shape
(dr,dc) = self.data.shape
(Tr,Tc) = (len(T),1) #T.shape #(len(T),1)
if (Tc!=1): raise ValueError("Temperature has to be a 2D array with one column.")
if (Tr!=dr): raise ValueError("Temperature and fitting data have to have the same number of rows.")
if self.type==self.INCOMPRESSIBLE_POLYNOMIAL:
if (numpy.max(x)==0.0): # x not given
x = numpy.array([0.0])
if (cc==1): # requires 1D coefficients
if (dc==1):
#z = numpy.polyfit(T, self.data[:,0], cr-1)
#self.coeffs = numpy.zeros((cr,1))
#self.coeffs[:,0] = z[::-1]
#print self.coeffs
self.coeffs = self.getCoeffs2d(T-Tbase, x-xbase, self.data, cr-1, cc-1)
#print self.coeffs
else:
raise ValueError("Cannot use 2D data with 1D references")
else:
raise ValueError("Cannot use 2D coefficients without concentration input")
else: # Assume 2D input
(xr,xc) = (1,len(x))#x.shape#(1,len(x))
if (xr!=1): raise ValueError("Concentration has to be a 2D array with one row.")
if (xc!=dc): raise ValueError("Concentration and fitting data have to have the same number of columns.")
#if (cr!=cc): raise ValueError("For now, you have to have the same number of exponents for x and y.")
self.coeffs = self.getCoeffs2d(T-Tbase, x-xbase, self.data, cr-1, cc-1)
#raise ValueError("Cannot use 2D coefficients yet")
print self.coeffs
else:
raise ValueError("Cannot fit that function.")
class SolutionData(LiquidData):
def getCoeffs1d(self, x,z,order):
if (len(x)<order+1):
raise ValueError("You have only {0} elements and try to fit {1} coefficients, please reduce the order.".format(len(x),order+1))
A = numpy.vander(x,order+1)[:,::-1]
#Anew = numpy.dot(A.T,A)
#znew = numpy.dot(A.T,z)
#coeffs = numpy.linalg.solve(Anew, znew)
coeffs = numpy.linalg.lstsq(A, z)[0]
return coeffs
def getCoeffs2d(self, x_in,y_in,z_in,x_order,y_order):
x_order += 1
y_order += 1
#To avoid overfitting, we only use the upper left triangle of the coefficient matrix
x_exp = range(x_order)
y_exp = range(y_order)
limit = max(x_order,y_order)
xy_exp = []
# Construct the upper left triangle of coefficients
for i in x_exp:
for j in y_exp:
if(i+j<limit): xy_exp.append((i,j))
# Construct input pairs
xx, yy = numpy.meshgrid(x_in,y_in,indexing='ij')
xx = numpy.array(xx.flat)
yy = numpy.array(yy.flat)
zz = numpy.array(z_in.flat)
x_num = len(x_in)
y_num = len(y_in)
cols = len(xy_exp)
eqns = x_num * y_num
#if (eqns<cols):
# raise ValueError("You have only {0} equations and try to fit {1} coefficients, please reduce the order.".format(eqns,cols))
if (x_num<x_order):
raise ValueError("You have only {0} x-entries and try to fit {1} x-coefficients, please reduce the x_order.".format(x_num,x_order))
if (y_num<y_order):
raise ValueError("You have only {0} y-entries and try to fit {1} y-coefficients, please reduce the y_order.".format(y_num,y_order))
# Build the functional matrix
A = numpy.zeros((eqns,cols))
for i in range(eqns): # row loop
for j, (xj,yj) in enumerate(xy_exp): # makes columns
A[i][j] = xx[i]**xj * yy[i]**yj
coeffs = numpy.linalg.lstsq(A, zz)[0]
#Rearrange coefficients to a matrix shape
C = numpy.zeros((x_order,y_order))
for i, (xi,yi) in enumerate(xy_exp): # makes columns
C[xi][yi] = coeffs[i]
return C
def toJSON(self):
j = {}
try:
j['coeffs'] = self.coeffs.tolist()
except:
j['coeffs'] = 'null'
j['type'] = self.type
return j
class SolutionData(object):
"""
A base class that defines all the variables needed
in order to make a proper fit for solution data.
in order to make a proper fit. You can copy this code
put in your data and add some documentation for where the
information came from.
"""
Name = None # Name of the current fluid
Desc = None # Name of the current fluid
Tmin = None # Minimum temperature in K
TminPsat = None # Minimum saturation temperature in K
Tmax = None # Maximum temperature in K
xmin = None # Minimum concentration in weight fraction
xmax = None # Minimum concentration in weight fraction
Tbase = None # Base temperature for data fit in K
xbase = None # Base concentration for data fit in weight fraction
def __init__(self):
self.name = None # Name of the current fluid
self.description = None # Description of the current fluid
self.reference = None # Reference data for the current fluid
self.Tmax = None # Maximum temperature in K
self.Tmin = None # Minimum temperature in K
self.xmax = None # Maximum concentration
self.xmin = None # Minimum concentration
self.TminPsat = None # Minimum saturation temperature in K
self.Tbase = 0.0 # Base value for temperature fits
self.xbase = 0.0 # Base value for concentration fits
# Data points
x = None # Concentration data points in weight fraction
T = None # Temperature for data points in K
rho = None # Density in kg/m3
c_p = None # Heat capacity in J/(kg.K)
lam = None # Thermal conductivity in W/(m.K)
mu_dyn = None # Dynamic viscosity in Pa.s
psat = None # Saturation pressure in Pa
Tfreeze = None # Freezing temperature in K
self.temperature = IncompressibleData() # Temperature for data points in K
self.concentration = IncompressibleData() # Concentration data points in weight fraction
self.density = IncompressibleData() # Density in kg/m3
self.specific_heat = IncompressibleData() # Heat capacity in J/(kg.K)
self.viscosity = IncompressibleData() # Dynamic viscosity in Pa.s
self.conductivity = IncompressibleData() # Thermal conductivity in W/(m.K)
self.saturation_pressure = IncompressibleData() # Saturation pressure in Pa
self.T_freeze = IncompressibleData() # Freezing temperature in K
self.volume2mass = IncompressibleData() # dd
self.mass2mole = IncompressibleData() # dd
self.xref = None
self.Tref = None
self.pref = None
self.href = None
self.sref = None
self.uref = None
self.rhoref = None
def rho (T, p, x=0.0, c=None):
if c==None: return numpy.polynomial.polynomial.polyval2d(T-self.Tbase, x-self.xbase, self.density.coeffs)
else: return numpy.polynomial.polynomial.polyval2d(T-self.Tbase, x-self.xbase, c)
def c (T, p, x=0.0, c=None):
if c==None: return numpy.polynomial.polynomial.polyval2d(T-self.Tbase, x-self.xbase, data.specific_heat.coeffs)
else: return numpy.polynomial.polynomial.polyval2d(T-self.Tbase, x-self.xbase, c)
def cp (T, p, x=0.0, c=None):
return self.c(T,p,x,c)
def cv (T, p, x=0.0, c=None):
return self.c(T,p,x,c)
#def s (T, p, x);
#def u (T, p, x);
def h (T, p, x=0.0):
return h_u(T,p,x)
#def visc(T, p, x=0.0){return baseFunction(viscosity, T, x);};
#def cond(T, p, x=0.0){return baseFunction(conductivity, T, x);};
#def psat(T, x=0.0){return baseFunction(p_sat, T, x);};
#def Tfreeze( p, x);
#def V2M (T, y);
#def M2M (T, x);
def h_u(T, p, x):
return u(T,p,x)+p/rho(T,p,x)-self.href
def u_h(T, p, x):
return h(T,p,x)-p/rho(T,p,x)+self.href
def set_reference_state(T0, p0, x0=0.0, h0=0.0, s0=0.0):
self.rhoref = rho(T0,p0,x0)
self.pref = p0
self.uref = h0 - p0/rhoref
self.uref = u(T0,p0,x0)
self.href = h0
self.sref = s0
self.href = h(T0,p0,x0)
self.sref = s(T0,p0,x0)
class SecCoolData(SolutionData):
"""
Ethanol-Water mixture according to Melinder book
Source: SecCool Software
"""
def __init__(self):
SolutionData.__init__(self)
self.reference = "SecCool software"
def convertSecCoolArray(self, array):
if len(array)!=18:
raise ValueError("The lenght is not equal to 18!")
array = numpy.array(array)
tmp = numpy.zeros((4,6))
tmp[0][0] = array[0]
tmp[0][1] = array[1]
tmp[0][2] = array[2]
tmp[0][3] = array[3]
tmp[0][4] = array[4]
tmp[0][5] = array[5]
tmp[1][0] = array[6]
tmp[1][1] = array[7]
tmp[1][2] = array[8]
tmp[1][3] = array[9]
tmp[1][4] = array[10]
#tmp[1][5] = array[11]
tmp[2][0] = array[11]
tmp[2][1] = array[12]
tmp[2][2] = array[13]
tmp[2][3] = array[14]
#tmp[2][4] = array[4]
#tmp[2][5] = array[5]
tmp[3][0] = array[15]
tmp[3][1] = array[16]
tmp[3][2] = array[17]
#tmp[3][3] = array[3]
#tmp[3][4] = array[4]
#tmp[3][5] = array[5]
return tmp
class SecCoolExample(SecCoolData):
"""
Ethanol-Water mixture according to Melinder book
Source: SecCool Software
"""
def __init__(self):
SecCoolData.__init__(self)
self.name = "SecCoolExample"
self.description = "Methanol solution"
#self.reference = "SecCool software"
self.Tmax = 20 + 273.15
self.Tmin = -50 + 273.15
self.xmax = 0.5
self.xmin = 0.0
self.TminPsat = 20 + 273.15
self.Tbase = -4.48 + 273.15
self.xbase = 31.57 / 100.0
self.density.type = self.density.INCOMPRESSIBLE_POLYNOMIAL
self.density.coeffs = self.convertSecCoolArray(numpy.array([
960.24665800,
1.2903839100,
-0.0161042520,
-0.0001969888,
1.131559E-05,
9.181999E-08,
-0.4020348270,
-0.0162463989,
0.0001623301,
4.367343E-06,
1.199000E-08,
-0.0025204776,
0.0001101514,
-2.320217E-07,
7.794999E-08,
9.937483E-06,
-1.346886E-06,
4.141999E-08]))
self.specific_heat.type = self.density.INCOMPRESSIBLE_POLYNOMIAL
self.specific_heat.coeffs = self.convertSecCoolArray(numpy.array([
3822.9712300,
-23.122409500,
0.0678775826,
0.0022413893,
-0.0003045332,
-4.758000E-06,
2.3501449500,
0.1788839410,
0.0006828000,
0.0002101166,
-9.812000E-06,
-0.0004724176,
-0.0003317949,
0.0001002032,
-5.306000E-06,
4.242194E-05,
2.347190E-05,
-1.894000E-06]))
self.conductivity.type = self.density.INCOMPRESSIBLE_POLYNOMIAL
self.conductivity.coeffs = self.convertSecCoolArray(numpy.array([
0.4082066700,
-0.0039816870,
1.583368E-05,
-3.552049E-07,
-9.884176E-10,
4.460000E-10,
0.0006629321,
-2.686475E-05,
9.039150E-07,
-2.128257E-08,
-5.562000E-10,
3.685975E-07,
7.188416E-08,
-1.041773E-08,
2.278001E-10,
4.703395E-08,
7.612361E-11,
-2.734000E-10]))
self.viscosity.type = self.density.INCOMPRESSIBLE_POLYNOMIAL
self.viscosity.coeffs = self.convertSecCoolArray(numpy.array([
1.4725525500,
0.0022218998,
-0.0004406139,
6.047984E-06,
-1.954730E-07,
-2.372000E-09,
-0.0411841566,
0.0001784479,
-3.564413E-06,
4.064671E-08,
1.915000E-08,
0.0002572862,
-9.226343E-07,
-2.178577E-08,
-9.529999E-10,
-1.699844E-06,
-1.023552E-07,
4.482000E-09]))
self.T_freeze.type = self.density.INCOMPRESSIBLE_POLYNOMIAL
self.T_freeze.coeffs = numpy.array([
27.755555600,
-22.973221700,
-1.1040507200,
-0.0120762281,
-9.343458E-05])
class PureExample(SolutionData):
def __init__(self):
SolutionData.__init__(self)
self.name = "PureExample"
self.description = "Heat transfer fluid TherminolD12 by Solutia"
self.reference = "Solutia data sheet"
self.Tmax = 260 + 273.15
self.Tmin = -85 + 273.15
self.TminPsat = 40 + 273.15
self.temperature.data = numpy.array([ 50 , 60 , 70 , 80 , 90 , 100 , 110 , 120 , 130 , 140 , 150 ])+273.15 # Kelvin
self.density.data = numpy.array([[ 740],[ 733],[ 726],[ 717],[ 710],[ 702],[ 695],[ 687],[ 679],[ 670],[ 662]]) # kg/m3
self.specific_heat.data = numpy.array([[ 2235],[ 2280],[ 2326],[ 2361],[ 2406],[ 2445],[ 2485],[ 2528],[ 2571],[ 2607],[ 2645]]) # J/kg-K
self.viscosity.data = numpy.array([[0.804],[ 0.704],[ 0.623],[ 0.556],[ 0.498],[ 0.451],[ 0.410],[ 0.374],[ 0.346],[ 0.317],[ 0.289]]) # Pa-s
self.conductivity.data = numpy.array([[0.105],[ 0.104],[ 0.102],[ 0.100],[ 0.098],[ 0.096],[ 0.095],[ 0.093],[ 0.091],[ 0.089],[ 0.087]]) # W/m-K
self.saturation_pressure.data = numpy.array([[ 0.5],[ 0.9],[ 1.4],[ 2.3],[ 3.9],[ 6.0],[ 8.7],[ 12.4],[ 17.6],[ 24.4],[ 33.2]]) # Pa
class SolutionExample(SolutionData):
def __init__(self):
SolutionData.__init__(self)
self.name = "SolutionExample"
self.description = "Ethanol ice slurry"
self.reference = "SecCool software"
#self.Tmax = 260 + 273.15
#self.Tmin = -85 + 273.15
#self.TminPsat = 40 + 273.15
self.temperature.data = numpy.array([ -45 , -40 , -35 , -30 , -25 , -20 , -15 , -10])+273.15 # Kelvin
self.concentration.data = numpy.array([ 5 , 10 , 15 , 20 , 25 , 30 , 35 ])/100.0 # mass fraction
self.density.data = numpy.array([
[1064.0, 1054.6, 1045.3, 1036.3, 1027.4, 1018.6, 1010.0],
[1061.3, 1052.1, 1043.1, 1034.3, 1025.6, 1017.0, 1008.6],
[1057.6, 1048.8, 1040.1, 1031.5, 1023.1, 1014.8, 1006.7],
[1053.1, 1044.6, 1036.2, 1028.0, 1019.9, 1012.0, 1004.1],
[1047.5, 1039.4, 1031.5, 1023.7, 1016.0, 1008.4, 1000.9],
[1040.7, 1033.2, 1025.7, 1018.4, 1011.2, 1004.0, 997.0],
[1032.3, 1025.3, 1018.5, 1011.7, 1005.1, 998.5, 992.0],
[1021.5, 1015.3, 1009.2, 1003.1, 997.1, 991.2, 985.4]]) # kg/m3

122
dev/IncompressibleLiquids/fit_incompressible.py
View File

@@ -6,6 +6,128 @@ from matplotlib.ticker import MaxNLocator
import os
import numpy as np
INCOMPRESSIBLE_NOT_SET = 0
INCOMPRESSIBLE_POLYNOMIAL = 1
INCOMPRESSIBLE_EXPONENTIAL = 2
INCOMPRESSIBLE_EXPPOLYNOMIAL = 3
class IncompressibleData(object):
def __init__(self):
self.type = INCOMPRESSIBLE_NOT_SET
self.coeffs = numpy.zeros(4,4)
def polyfit2d(x, y, z, order=3, linear=False):
"""Two-dimensional polynomial fit. Based uppon code provided by
Joe Kington.
References:
http://stackoverflow.com/questions/7997152/
python-3d-polynomial-surface-fit-order-dependent/7997925#7997925
"""
ncols = (order + 1)**2
G = numpy.zeros((x.size, ncols))
ij = itertools.product(range(order+1), range(order+1))
for k, (i,j) in enumerate(ij):
G[:,k] = x**i * y**j
if linear & (i != 0.) & (j != 0.):
G[:, k] = 0
m, _, _, _ = numpy.linalg.lstsq(G, z)
return m
def polyval2d(x, y, m):
"""Values to two-dimensional polynomial fit. Based uppon code
provided by Joe Kington.
"""
order = int(numpy.sqrt(len(m))) - 1
ij = itertools.product(range(order+1), range(order+1))
z = numpy.zeros_like(x)
for a, (i,j) in zip(m, ij):
z += a * x**i * y**j
return z
class IncompSolutionFit(object):
"""
A class for fitting data sheet data to predefined functions.
Some functions are only used during the fitting procedure.
Note that the order in which you fit the different properties
might impact the coefficients. Usually, the fitting order should be:
1) Density
2) Heat capacity
3) Thermal conductivity
4) Viscosity
5) Vapour pressure
"""
def __init__(self):
self.DEBUG = False
self._name= None
self._description= None
self._reference= None
self._Tmin= None
self._Tmax= None
self._xmin= None
self._xmax= None
self._TminPsat= None
self._xref= None
self._Tref= None
self._pref= None
self._href= None
self._sref= None
self._ure= None
self._rhoref= None
self._xbase= None
self._Tbase= None
# parameters for the different fits
self._density = numpy.ones(4,4) # Typically 4 parameters
self._cHeatCapacity = numpy.ones(4,4) # Typically 4 parameters
self._cTConductivity = numpy.ones(3,3) # Typically 3 parameters
self._cViscosity = numpy.ones(3,3) # Typically 3 parameters
self._cPsat = numpy.ones(3,3) # Typically 3 parameters
self._density;
self._specific_heat;
self._viscosity;
self._conductivity;
self._p_sat;
self._T_freeze;
self._volToMass;
self._massToMole;
def set_reference_state(T0, p0, x0=0.0, h0=0.0, s0=0.0):
this._rhoref = rho(T0,p0,x0);
this._pref = p0;
this._uref = h0 - p0/rhoref;
this._uref = u(T0,p0,x0);
this._href = h0; // set new reference value
this._sref = s0; // set new reference value
this._href = h(T0,p0,x0); // adjust offset to fit to equations
this._sref = s(T0,p0,x0); // adjust offset to fit to equations
}
class IncompLiquidFit(object):
"""
A class for fitting data sheet data to predefined functions.

333
dev/IncompressibleLiquids/fit_incompressible_simple.py
View File

@@ -1,333 +0,0 @@
import matplotlib.pyplot
import numpy
from scipy.optimize._minimize import minimize
class IncompLiquidFit_simple(object):
"""
A class for fitting data sheet data to predefined functions.
Some functions are only used during the fitting procedure.
Note that the order in which you fit the different properties
might impact the coefficients. Usually, the fitting order should be:
1) Density
2) Heat capacity
3) Thermal conductivity
4) Viscosity
5) Vapour pressure
"""
def __init__(self):
self.DEBUG = False
# parameters for the different fits
self._cDensity = numpy.ones(2) # Typically 2 parameters
self._cHeatCapacity = numpy.ones(2) # Typically 2 parameters
self._cTConductivity = numpy.ones(2) # Typically 2 parameters
self._cViscosity = numpy.ones(3) # Typically 3 parameters
self._cPsat = numpy.ones(3) # Typically 3 parameters
# bounds for fit
self._Tmin = None
self._TminPsat = None
self._Tmax = None
# some flags to set
self._TinC = False # Temperature in Celsius
self._DynVisc = True # Data for dynamic viscosity
self._minPoints = 5
def setParams(self,fluid):
if fluid=='init':
# initial parameters for the different fits
self._cDensity = [+9.2e+2, -0.5e+0]
self._cHeatCapacity = [+1.0e+0, +2.9e-3]
self._cTConductivity = [+1.1e-4, +7.8e-7]
self._cViscosity = [+7.1e+0, -2.3e-2, +3.4e-5]
self._cPsat = [-5.3e+3, -3.2e+1, -1.6e+1]
else:
raise (ValueError("No coefficients available for "+str(fluid)))
def _checkT(self,T=0):
Tmin = self.Props('Tmin')
Tmax = self.Props('Tmax')
if Tmin is None:
raise (ValueError("Please specify the minimum temperature."))
if Tmax is None:
raise (ValueError("Please specify the maximum temperature."))
if not (Tmin<=T<=Tmax):
raise (ValueError("Temperature out of range: "+str(T)+" not in "+str(Tmin)+"-"+str(Tmax)+". "))
def _checkP(self,T=0,P=0):
Psat = self.Props('Psat',T=T)
if P<Psat:
raise (ValueError("Equations are valid for liquid phase only: "+str(P)+" < "+str(Psat)+". "))
def _checkTP(self,T=0,P=0):
self._checkT(T=T)
self._checkP(T=T, P=P)
def _basePolynomial(self,coefficients,x):
""" Base function to produce polynomials of
order len(coefficients) with the coefficients
"""
result = 0.
for i in range(len(coefficients)):
result += coefficients[i] * x**i
return result
def _baseExponential(self,coefficients,x,num=3):
""" Base function to produce exponential
with defined coefficients
"""
if len(coefficients)==num:
if num==3:
return numpy.exp(coefficients[0]/(x+coefficients[1]) + coefficients[2])
else:
print "Error!"
else:
print "Error!"
def Props(self,out,T=0,P=0):
if out=='D':
self._checkTP(T=T,P=P)
return self._basePolynomial(self._cDensity,T)
elif out=='C':
self._checkTP(T=T,P=P)
return self._basePolynomial(self._cHeatCapacity,T)
elif out=='L':
self._checkTP(T=T,P=P)
return self._basePolynomial(self._cTConductivity,T)
elif out=='V':
self._checkTP(T=T,P=P)
return numpy.exp(self._basePolynomial(self._cViscosity,T))
elif out=='Psat':
self._checkT(T=T)
if T<self._TminPsat:
return 1e-14
return self._baseExponential(self._cPsat,T)
elif out=='Tmin':
return self._Tmin
elif out=='Tmax':
return self._Tmax
else:
raise (ValueError("Error: You used an unknown output qualifier."))
def _PropsFit(self,coefficients,inVal,T=0):
"""
Calculates a property from a given set of
coefficents for a certain temperature. Is used
to obtain data to feed to the optimisation
procedures.
"""
if inVal=='D':
self._checkT(T=T)
return self._basePolynomial(coefficients,T)
elif inVal=='C':
self._checkT(T=T)
return self._basePolynomial(coefficients,T)
elif inVal=='L':
self._checkT(T=T)
return self._basePolynomial(coefficients,T)
elif inVal=='V':
self._checkT(T=T)
return numpy.exp(self._basePolynomial(coefficients,T))
elif inVal=='Psat':
self._checkT(T=T)
if T<self._TminPsat:
return 1e-14
return self._baseExponential(coefficients,T)
else:
raise (ValueError("Error: You used an unknown property qualifier."))
def getCoefficients(self,inVal):
"""
Get the array with coefficients.
"""
if inVal=='D':
return self._cDensity
elif inVal=='C':
return self._cHeatCapacity
elif inVal=='L':
return self._cTConductivity
elif inVal=='V':
return self._cViscosity
elif inVal=='Psat':
return self._cPsat
else:
raise (ValueError("Error: You used an unknown property qualifier."))
def setCoefficients(self,inVal,coeffs):
"""
Set the array of coefficients.
"""
if inVal=='D':
self._cDensity = coeffs
elif inVal=='C':
self._cHeatCapacity = coeffs
elif inVal=='L':
self._cTConductivity = coeffs
elif inVal=='V':
self._cViscosity = coeffs
elif inVal=='Psat':
self._cPsat = coeffs
else:
raise (ValueError("Error: You used an unknown property qualifier."))
def setTmin(self,T):
self._Tmin = T
def setTmax(self,T):
self._Tmax = T
def setTminPsat(self,T):
self._TminPsat = T
def fitCoefficients(self,xName,T=[],xData=[]):
if (len(T)!=len(xData)):
raise (ValueError("Error: There has to be the same number of temperature and data points."))
if len(T)<self._minPoints:
raise (ValueError("Error: You should use at least "+str(self._minPoints)+" points."))
def fun(coefficients,xName,T,xData):
# Values for conductivity are very small,
# algorithms prefer larger values
if xName=='L':
calculated = numpy.array([self._PropsFit(coefficients,xName,T=Ti) for Ti in T])*1e6
data = numpy.array(xData)*1e6
# Fit logarithms for viscosity and saturation pressure
elif xName=='V' or xName=='Psat':
calculated = numpy.log([self._PropsFit(coefficients,xName,T=Ti) for Ti in T])
data = numpy.log(xData)
else:
calculated = numpy.array([self._PropsFit(coefficients,xName,T=Ti) for Ti in T])
data = numpy.array(xData)
res = numpy.sum((calculated-data)**2.)
return res
initValues = self.getCoefficients(xName)[:]
arguments = (xName,T,xData)
options = {'maxiter': 1e4, 'maxfev': 1e6}
res = minimize(fun, initValues, method='Powell', args=arguments, tol=1e-15, options=options)
if res.success:
return res.x
else:
print res
return False
### Load the data
from data_incompressible import TherminolD12 as DataContainer
data = DataContainer()
### Some test case
liqObj = IncompLiquidFit_simple()
liqObj.setParams("init")
liqObj.setTmin(data.Tmin)
liqObj.setTminPsat(data.TminPsat)
liqObj.setTmax(data.Tmax)
numpy.set_printoptions(formatter={'float': lambda x: format(x, '+1.10E')})
print "minimum T: "+str(data.Tmin)
print "maximum T: "+str(data.Tmax)
print "min T pSat:"+str(data.TminPsat)
print
# row and column sharing for test plots
f, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = matplotlib.pyplot.subplots(3, 2, sharex='col')
f.set_size_inches(matplotlib.pyplot.figaspect(1.2)*1.5)
### This is the actual fitting
tData = data.T
cData = tData - 273.15
Pin = 10e5 # Dummy pressure
inVal = 'D'
xData = data.rho
oldCoeffs = liqObj.getCoefficients(inVal)
newCoeffs = liqObj.fitCoefficients(inVal,T=tData,xData=xData)
print "Density, old: "+str(oldCoeffs)
print "Density, new: "+str(newCoeffs)
print
liqObj.setCoefficients(inVal,newCoeffs)
fData = numpy.array([liqObj.Props(inVal, T=Tin, P=Pin) for Tin in tData])
ax1.plot(cData, fData, label="CoolProp")
ax1.plot(cData, xData, label="Data Sheet")
ax1.set_ylabel(r'$\mathregular{Density\/(kg\/m^{-3})}$')
inVal = 'C'
xData = data.c_p
oldCoeffs = liqObj.getCoefficients(inVal)
newCoeffs = liqObj.fitCoefficients(inVal,T=tData,xData=xData)
print "Heat c., old: "+str(oldCoeffs)
print "Heat c., new: "+str(newCoeffs)
print
liqObj.setCoefficients(inVal,newCoeffs)
fData = numpy.array([liqObj.Props(inVal, T=Tin, P=Pin) for Tin in tData])
ax2.plot(cData, fData, label="CoolProp")
ax2.plot(cData, xData, label="Data Sheet")
ax2.set_ylabel(r'$\mathregular{Heat\/Cap.\/(kJ\/kg^{-1}\/K^{-1})}$')
inVal = 'L'
xData = data.lam
oldCoeffs = liqObj.getCoefficients(inVal)
newCoeffs = liqObj.fitCoefficients(inVal,T=tData,xData=xData)
print "Th. Co., old: "+str(oldCoeffs)
print "Th. Co., new: "+str(newCoeffs)
print
liqObj.setCoefficients(inVal,newCoeffs)
fData = numpy.array([liqObj.Props(inVal, T=Tin, P=Pin) for Tin in tData])
ax3.plot(cData, fData*1e6 , label="CoolProp")
ax3.plot(cData, xData*1e6 , label="Data Sheet")
ax3.set_ylabel(r'$\mathregular{Th.\/Cond.\/(mW\/m^{-1}\/K^{-1})}$')
inVal = 'V'
xData = data.mu_dyn
oldCoeffs = liqObj.getCoefficients(inVal)
newCoeffs = liqObj.fitCoefficients(inVal,T=tData,xData=xData)
print "Viscos., old: "+str(oldCoeffs)
print "Viscos., new: "+str(newCoeffs)
print
liqObj.setCoefficients(inVal,newCoeffs)
fData = numpy.array([liqObj.Props(inVal, T=Tin, P=Pin) for Tin in tData])
ax4.plot(cData, fData*1e3, label="CoolProp")
ax4.plot(cData, xData*1e3, label="Data Sheet")
ax4.set_ylabel(r'$\mathregular{Dyn.\/Viscosity\/(mPa\/s)}$')
ax4.set_yscale('log')
inVal = 'Psat'
xData = data.psat
oldCoeffs = liqObj.getCoefficients(inVal)
newCoeffs = liqObj.fitCoefficients(inVal,T=tData,xData=xData)
print "P sat. , old: "+str(oldCoeffs)
print "P sat. , new: "+str(newCoeffs)
print
liqObj.setCoefficients(inVal,newCoeffs)
fData = numpy.array([liqObj.Props(inVal, T=Tin, P=Pin) for Tin in tData])
ax5.plot(cData, fData, label="CoolProp")
ax5.plot(cData, xData, label="Data Sheet")
ax5.set_ylabel(r'$\mathregular{Vap.\/Pressure\/(kPa)}$')
ax5.set_yscale('log')
ax5.set_xlabel(ur'$\mathregular{Temperature\/(\u00B0C)}$')
ax6.set_xlabel(ur'$\mathregular{Temperature\/(\u00B0C)}$')
ax1.legend(loc=3)
matplotlib.pyplot.tight_layout()
matplotlib.pyplot.savefig(data.Name+"_simple.pdf")

104
dev/IncompressibleLiquids/json_incompressible.py Normal file
View File

@@ -0,0 +1,104 @@
import numpy
from data_incompressible import SolutionData
class SolutionDataWriter(object):
"""
A base class that defines all the variables needed
in order to make a proper fit. You can copy this code
put in your data and add some documentation for where the
information came from.
"""
#def __init__(self):
# SolutionData.__init__(self)
def fitAll(self):
def toJSON(self,data):
jobj = {}
jobj['name'] = data.name # Name of the current fluid
jobj['description'] = data.description # Description of the current fluid
jobj['reference'] = data.reference # Reference data for the current fluid
jobj['Tmax'] = data.Tmax # Maximum temperature in K
jobj['Tmin'] = data.Tmin # Minimum temperature in K
jobj['xmax'] = data.xmax # Maximum concentration
jobj['xmin'] = data.xmin # Minimum concentration
jobj['TminPsat'] = data.TminPsat # Minimum saturation temperature in K
jobj['Tbase'] = data.Tbase # Base value for temperature fits
jobj['xbase'] = data.xbase # Base value for concentration fits
#data.temperature # Temperature for data points in K
#data.concentration # Concentration data points in weight fraction
jobj['density'] = data.density.toJSON() # Density in kg/m3
jobj['specific_heat'] = data.specific_heat.toJSON() # Heat capacity in J/(kg.K)
jobj['viscosity'] = data.viscosity.toJSON() # Dynamic viscosity in Pa.s
jobj['conductivity'] = data.conductivity.toJSON() # Thermal conductivity in W/(m.K)
jobj['saturation_pressure'] = data.saturation_pressure.toJSON() # Saturation pressure in Pa
jobj['T_freeze'] = data.T_freeze.toJSON() # Freezing temperature in K
jobj['volume2mass'] = data.volume2mass.toJSON() # dd
jobj['mass2mole'] = data.mass2mole.toJSON() # dd
import json
dump = json.dumps(jobj, indent = 2, sort_keys = True)
print dump
fp = open(jobj['name']+'.json', 'w')
fp.write(dump)
fp.close()
if __name__ == '__main__':
writer = SolutionDataWriter()
# data = SecCoolExample()
# writer.toJSON(data)
data = PureExample()
data.density.coeffs = numpy.zeros((4,1))
data.density.type = data.density.INCOMPRESSIBLE_POLYNOMIAL
data.density.data = data.density.data[0:4]
data.temperature.data = data.temperature.data[0:4]
data.density.fit(data.temperature.data)
#writer.toJSON(data)
print data.density.data[0][0]
print numpy.polynomial.polynomial.polyval2d(data.temperature.data[0], 0, data.density.coeffs)
print data.density.data[1][0]
print numpy.polynomial.polynomial.polyval2d(data.temperature.data[1], 0, data.density.coeffs)
data = SolutionExample()
data.density.coeffs = numpy.zeros((3,2))
data.density.type = data.density.INCOMPRESSIBLE_POLYNOMIAL
data.density.data = data.density.data[0:4][:,0:2]
data.temperature.data = data.temperature.data[0:4]
data.density.fit(data.temperature.data,data.concentration.data[0:2])
#writer.toJSON(data)
print data.density.data[0][0]
print numpy.polynomial.polynomial.polyval2d(data.temperature.data[0], data.concentration.data[0], data.density.coeffs)
print data.density.data[1][1]
print numpy.polynomial.polynomial.polyval2d(data.temperature.data[1], data.concentration.data[1], data.density.coeffs)

4
include/CoolPropTools.h
View File

@@ -244,8 +244,8 @@
max = min;
min = fabs(A);
}
if (max>D) return ( ( 1.0-min/max*1e0 ) < D );
else throw CoolProp::ValueError(format("Too small numbers: %f cannot be tested with an accepted error of %f. ",max,D));
if (max>DBL_EPSILON*1e3) return ( ( 1.0-min/max*1e0 ) < D );
else throw CoolProp::ValueError(format("Too small numbers: %f cannot be tested with an accepted error of %f for a machine precision of %f. ",max,D,DBL_EPSILON));
};
bool inline check_abs(double A, double B){
return check_abs(A,B,1e5*DBL_EPSILON);

234
include/IncompressibleFluid.h
View File

@@ -19,70 +19,200 @@
#include <assert.h>
#include <iterator>
#include <Eigen/Core>
#include "PolyMath.h"
#include "MatrixMath.h"
namespace CoolProp {
struct IncompressiblePolynomialData{
std::vector<double> coeffs;
};
struct IncompressibleViscosityVariables{
enum IncompressibleViscosityEnum{
INCOMPRESSIBLE_VISCOSITY_POLYNOMIAL,
INCOMPRESSIBLE_VISCOSITY_NOT_SET
};
int type;
IncompressiblePolynomialData poly;
IncompressibleViscosityVariables(){type = INCOMPRESSIBLE_VISCOSITY_NOT_SET;};
};
struct IncompressibleConductivityVariables{
enum IncompressibleConductivityEnum{
INCOMPRESSIBLE_CONDUCTIVITY_POLYNOMIAL,
INCOMPRESSIBLE_CONDUCTIVITY_NOT_SET
};
int type;
IncompressiblePolynomialData poly;
IncompressibleConductivityVariables(){type = INCOMPRESSIBLE_CONDUCTIVITY_NOT_SET;};
};
struct IncompressibleDensityVariables{
enum IncompressibleDensityEnum{
INCOMPRESSIBLE_DENSITY_POLYNOMIAL,
INCOMPRESSIBLE_DENSITY_NOT_SET
};
int type;
IncompressiblePolynomialData poly;
IncompressibleDensityVariables(){type = INCOMPRESSIBLE_DENSITY_NOT_SET;};
};
struct IncompressibleSpecificHeatVariables{
enum IncompressibleSpecificHeatEnum{
INCOMPRESSIBLE_SPECIFIC_HEAT_POLYNOMIAL,
INCOMPRESSIBLE_SPECIFIC_HEAT_NOT_SET
};
int type;
IncompressiblePolynomialData poly;
IncompressibleSpecificHeatVariables(){type = INCOMPRESSIBLE_SPECIFIC_HEAT_NOT_SET;};
struct IncompressibleData {
int type;
enum IncompressibleTypeEnum {
INCOMPRESSIBLE_NOT_SET,
INCOMPRESSIBLE_POLYNOMIAL,
INCOMPRESSIBLE_EXPONENTIAL,
INCOMPRESSIBLE_EXPPOLYNOMIAL
};
Eigen::MatrixXd coeffs; //TODO: Can we store the Eigen::Matrix objects more efficiently?
//std::vector<std::vector<double> > coeffs;
IncompressibleData() {
type = INCOMPRESSIBLE_NOT_SET;
};
};
/// A thermophysical property provider for critical and reducing values as well as derivatives of Helmholtz energy
/**
This fluid instance is populated using an entry from a JSON file
*/
struct IncompressibleFluid {
class IncompressibleFluid{
std::string name;
std::string reference;
protected:
std::string name;
std::string description;
std::string reference;
double Tmin, Tmax;
double Tmin, Tmax;
double xmin, xmax;
IncompressibleViscosityVariables viscosity;
IncompressibleConductivityVariables conductivity;
IncompressibleSpecificHeatVariables specific_heat;
IncompressibleDensityVariables density;
double TminPsat;
double xref, Tref, pref;
double href, sref;
double uref, rhoref;
double xbase, Tbase;
IncompressibleData density;
IncompressibleData specific_heat;
IncompressibleData viscosity;
IncompressibleData conductivity;
IncompressibleData p_sat;
IncompressibleData T_freeze;
IncompressibleData volToMass;
IncompressibleData massToMole;
Polynomial2DFrac poly;
// Forward declaration of the some internal functions
//double h_u(double T, double p, double x);
//double u_h(double T, double p, double x);
public:
IncompressibleFluid(){};
virtual ~IncompressibleFluid(){};
std::string getName() const {return name;}
std::string get_name() const {return getName();}// For backwards-compatibility.
std::string getDescription() const {return description;}
std::string getReference() const {return reference;}
double getTmax() const {return Tmax;}
double getTmin() const {return Tmin;}
double getxmax() const {return xmax;}
double getxmin() const {return xmin;}
double getTminPsat() const {return TminPsat;}
double getTref() const {return Tref;}
double getpref() const {return pref;}
double getxref() const {return xref;}
double gethref() const {return href;}
double getsref() const {return sref;}
double getTbase() const {return Tbase;}
double getxbase() const {return xbase;}
void setName(std::string name) {this->name = name;}
void setDescription(std::string description) {this->description = description;}
void setReference(std::string reference) {this->reference = reference;}
void setTmax(double Tmax) {this->Tmax = Tmax;}
void setTmin(double Tmin) {this->Tmin = Tmin;}
void setxmax(double xmax) {this->xmax = xmax;}
void setxmin(double xmin) {this->xmin = xmin;}
void setTminPsat(double TminPsat) {this->TminPsat = TminPsat;}
//void setTref(double Tref) {this->Tref = Tref;}
//void setpref(double pref) {this->pref = pref;}
//void setxref(double xref) {this->xref = xref;}
void set_reference_state(double T0, double p0, double x0, double h0, double s0);
void setTbase(double Tbase) {this->Tbase = Tbase;}
void setxbase(double xbase) {this->xbase = xbase;}
/// Setters for the coefficients
void setDensity(IncompressibleData density){this->density = density;}
void setSpecificHeat(IncompressibleData specific_heat){this->specific_heat = specific_heat;}
void setViscosity(IncompressibleData viscosity){this->viscosity = viscosity;}
void setConductivity(IncompressibleData conductivity){this->conductivity = conductivity;}
void setPsat(IncompressibleData p_sat){this->p_sat = p_sat;}
void setTfreeze(IncompressibleData T_freeze){this->T_freeze = T_freeze;}
void setVolToMass(IncompressibleData volToMass){this->volToMass = volToMass;}
void setMassToMole(IncompressibleData massToMole){this->massToMole = massToMole;}
/// A function to check coefficients and equation types.
void validate();
protected:
/// Base function that handles the custom data type, just a place holder to show the structure.
double baseFunction(IncompressibleData data, double x_in, double y_in);
public:
/* All functions need T and p as input. Might not
* be necessary, but gives a clearer structure.
*/
/// Density as a function of temperature, pressure and composition.
double rho (double T, double p, double x=0.0){return baseFunction(density, T, x);};
/// Heat capacities as a function of temperature, pressure and composition.
double c (double T, double p, double x=0.0){return baseFunction(specific_heat, T, x);};
double cp (double T, double p, double x=0.0){return c(T,p,x);};
double cv (double T, double p, double x=0.0){return c(T,p,x);};
/// Entropy as a function of temperature, pressure and composition.
double s (double T, double p, double x);
/// Internal energy as a function of temperature, pressure and composition.
double u (double T, double p, double x);
/// Enthalpy as a function of temperature, pressure and composition.
double h (double T, double p, double x=0.0){return h_u(T,p,x);};
/// Viscosity as a function of temperature, pressure and composition.
double visc(double T, double p, double x=0.0){return baseFunction(viscosity, T, x);};
/// Thermal conductivity as a function of temperature, pressure and composition.
double cond(double T, double p, double x=0.0){return baseFunction(conductivity, T, x);};
/// Saturation pressure as a function of temperature and composition.
double psat(double T, double x=0.0){return baseFunction(p_sat, T, x);};
/// Freezing temperature as a function of pressure and composition.
double Tfreeze( double p, double x);
/// Conversion from volume-based to mass-based composition.
double V2M (double T, double y);
/// Conversion from mass-based to mole-based composition.
double M2M (double T, double x);
protected:
/* Define internal energy and enthalpy as functions of the
* other properties to provide data in case there are no
* coefficients.
*/
/// Enthalpy from u, p and rho.
/** Calculate enthalpy as a function of temperature and
* pressure employing functions for internal energy and
* density. Provides consistent formulations. */
double h_u(double T, double p, double x) {
return u(T,p,x)+p/rho(T,p,x)-href;
};
/// Internal energy from h, p and rho.
/** Calculate internal energy as a function of temperature
* and pressure employing functions for enthalpy and
* density. Provides consistent formulations. */
double u_h(double T, double p, double x) {
return h(T,p,x)-p/rho(T,p,x)+href;
};
/*
* Some more functions to provide a single implementation
* of important routines.
* We start with the check functions that can validate input
* in terms of pressure p, temperature T and composition x.
*/
/// Check validity of temperature input.
/** Compares the given temperature T to the result of a
* freezing point calculation. This is not necessarily
* defined for all fluids, default values do not cause errors. */
bool checkT(double T, double p, double x);
/// Check validity of pressure input.
/** Compares the given pressure p to the saturation pressure at
* temperature T and throws and exception if p is lower than
* the saturation conditions.
* The default value for psat is -1 yielding true if psat
* is not redefined in the subclass.
* */
bool checkP(double T, double p, double x);
/// Check validity of composition input.
/** Compares the given composition x to a stored minimum and
* maximum value. Enforces the redefinition of xmin and
* xmax since the default values cause an error. */
bool checkX(double x);
/// Check validity of temperature, pressure and composition input.
bool checkTPX(double T, double p, double x);
};
} /* namespace CoolProp */
#endif /* COOLPROPFLUID_H_ */
#endif /* INCOMPRESSIBLEFLUID_H_ */

383
include/MatrixMath.h
View File

@@ -10,8 +10,328 @@
#include <sstream>
#include "float.h"
#include <Eigen/Core>
/// A wrapper around std::vector
/** This wrapper makes the standard vector multi-dimensional.
* A useful thing even though we might not need it that
* much. However, it makes the code look better and the
* polynomial class really is a mess...
* Source: http://stackoverflow.com/questions/13105514/n-dimensional-vector
*/
template<size_t dimcount, typename T> struct VectorNd {
typedef std::vector< typename VectorNd<dimcount-1, T>::type > type;
};
template<typename T> struct VectorNd<0,T> {
typedef T type;
};
namespace CoolProp{
/// Convert vectors and matrices
/** Conversion functions for the different kinds of object-like
* parameters. This might be obsolete at a later stage, but now
* it is still needed.
*/
/// @param coefficients matrix containing the ordered coefficients
template <typename T> std::vector<T> eigen_to_vec1D(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &coefficients, int axis = 0){
std::vector<T> result;
size_t r = coefficients.rows(), c = coefficients.cols();
if (axis==0) {
if (c!=1) throw ValueError(format("Your matrix has the wrong dimensions: %d,%d",r,c));
result.resize(r);
for (size_t i = 0; i < r; ++i) {
result[i] = coefficients(i,0);
}
} else if (axis==1) {
if (r!=1) throw ValueError(format("Your matrix has the wrong dimensions: %d,%d",r,c));
result.resize(c);
for (size_t i = 0; i < c; ++i) {
result[i] = coefficients(0,i);
}
} else {
throw ValueError(format("You have to provide axis information: %d is not valid. ",axis));
}
return result;
}
/// @param coefficients matrix containing the ordered coefficients
template <class T> std::vector<std::vector<T> > eigen_to_vec(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &coefficients){
// Eigen uses columns as major axis, this might be faster than the row iteration.
// However, the 2D vector stores things differently, no idea what is faster...
std::vector<std::vector<T> > result;
size_t r = coefficients.rows(), c = coefficients.cols();
result.resize(r, std::vector<T>(c, 0)); // extends vector if necessary
for (size_t i = 0; i < r; ++i) {
result[i].resize(c, 0);
for (size_t j = 0; j < c; ++j) {
result[i][j] = coefficients(i,j);
}
}
return result;
}
/// @param coefficients matrix containing the ordered coefficients
template <class T> Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> vec_to_eigen(const std::vector<std::vector<T> > &coefficients){
size_t nRows = num_rows(coefficients), nCols = num_cols(coefficients);
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> result(nRows,nCols);
for (size_t i = 0; i < nCols; ++i) {
for (size_t j = 0; j < nRows; ++j) {
result(j,i) = coefficients[j][i];
}
}
return result;
}
/// @param coefficients matrix containing the ordered coefficients
template <class T> Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> vec_to_eigen(const std::vector<T> &coefficients, int axis = 0){
size_t nRows = num_rows(coefficients);
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> result;
if (axis==0) result.resize(nRows,1);
else if (axis==1) result.resize(1,nRows);
else throw ValueError(format("You have to provide axis information: %d is not valid. ",axis));
for (size_t i = 0; i < nRows; ++i) {
if (axis==0) result(i,0) = coefficients[i];
if (axis==1) result(0,i) = coefficients[i];
}
return result;
}
/// @param coefficient
template <class T> Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> vec_to_eigen(const T &coefficient){
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> result = Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>(1,1);
result(0,0) = coefficient;
return result;
}
/// Remove rows and columns from matrices
/** A set of convenience functions inspired by http://stackoverflow.com/questions/13290395/how-to-remove-a-certain-row-or-column-while-using-eigen-library-c
* but altered to respect templates.
*/
template< class T> void removeRow(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &matrix, std::size_t rowToRemove){
//template<class T> void removeRow(Eigen::MatrixXd& matrix, std::size_t rowToRemove){
//void removeRow(Eigen::MatrixXd& matrix, unsigned int rowToRemove){
//template <typename Derived> void removeRow(Eigen::MatrixBase<Derived> &matrix, std::size_t rowToRemove){
std::size_t numRows = matrix.rows()-1;
std::size_t numCols = matrix.cols();
if( rowToRemove < numRows ){
matrix.block(rowToRemove,0,numRows-rowToRemove,numCols) = matrix.block(rowToRemove+1,0,numRows-rowToRemove,numCols);
} else {
if( rowToRemove > numRows ){
throw ValueError(format("Your matrix does not have enough rows, %d is not greater or equal to %d.",numRows,rowToRemove));
}
// Do nothing, resize removes the last row
}
matrix.conservativeResize(numRows,numCols);
}
template <class T> void removeColumn(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &matrix, std::size_t colToRemove){
//template<class T> void removeColumn(Eigen::MatrixXd& matrix, std::size_t colToRemove){
//void removeColumn(Eigen::MatrixXd& matrix, unsigned int colToRemove){
//template <typename Derived> void removeColumn(Eigen::MatrixBase<Derived> &matrix, std::size_t colToRemove){
std::size_t numRows = matrix.rows();
std::size_t numCols = matrix.cols()-1;
if( colToRemove < numCols ) {
matrix.block(0,colToRemove,numRows,numCols-colToRemove) = matrix.block(0,colToRemove+1,numRows,numCols-colToRemove);
} else {
if( colToRemove > numCols ) {
throw ValueError(format("Your matrix does not have enough columns, %d is not greater or equal to %d.",numCols,colToRemove));
}
// Do nothing, resize removes the last column
}
matrix.conservativeResize(numRows,numCols);
}
///// @param coefficients matrix containing the ordered coefficients
//template <class T> Eigen::Matrix<T, Eigen::Dynamic,Eigen::Dynamic> convert(const std::vector<std::vector<T> > &coefficients){
// size_t nRows = num_rows(coefficients), nCols = num_cols(coefficients);
// Eigen::MatrixBase<T> result(nRows,nCols);
// for (size_t i = 0; i < nCols; ++i) {
// for (size_t j = 0; j < nRows; ++j) {
// result(j,i) = coefficients[j][i];
// }
// }
// return result;
//}
//
///// @param coefficients matrix containing the ordered coefficients
//template <class T, int R, int C> void convert(const std::vector<std::vector<T> > &coefficients, Eigen::Matrix<T,R,C> &result){
// size_t nRows = num_rows(coefficients), nCols = num_cols(coefficients);
// //Eigen::MatrixBase<T> result(nRows,nCols);
// for (size_t i = 0; i < nCols; ++i) {
// for (size_t j = 0; j < nRows; ++j) {
// result(j,i) = coefficients[j][i];
// }
// }
// //return result;
//}
//
//template <class T> void convert(const std::vector<std::vector<T> > &coefficients, Eigen::MatrixBase<T> &result){
// size_t nRows = num_rows(coefficients), nCols = num_cols(coefficients);
// //Eigen::MatrixBase<T> result;
// //if ((R!=nRows) || (C!=nCols))
// result.resize(nRows,nCols);
// for (size_t i = 0; i < nCols; ++i) {
// for (size_t j = 0; j < nRows; ++j) {
// result(j,i) = coefficients[j][i];
// }
// }
// //return result;
//}
//template <class Derived>
//inline void func1(MatrixBase<Derived> &out_mat ){
// // Do something then return a matrix
// out_mat = ...
//}
//template <class Derived>
//Eigen::Matrix<class Derived::Scalar, Derived::RowsAtCompileTime, Derived::ColsAtCompileTime>
//Multiply(const Eigen::MatrixBase<DerivedA>& p1,
// const Eigen::MatrixBase<DerivedB>& p2)
//{
// return (p1 * p2);
//}
//
//
//template <typename DerivedA, typename DerivedB>
//Eigen::Matrix<typename DerivedA::Scalar, DerivedA::RowsAtCompileTime, DerivedB::ColsAtCompileTime>
//Multiply(const Eigen::MatrixBase<DerivedA>& p1,
// const Eigen::MatrixBase<DerivedB>& p2)
//{
// return (p1 * p2);
//}
/// Templates for printing numbers, vectors and matrices
static const char* stdFmt = "%8.3f";
///Templates for turning vectors (1D-matrices) into strings
template<class T> std::string vec_to_string(const std::vector<T> &a, const char *fmt) {
if (a.size()<1) return std::string("");
std::stringstream out;
out << "[ " << format(fmt,a[0]);
for (size_t j = 1; j < a.size(); j++) {
out << ", " << format(fmt, a[j]);
}
out << " ]";
return out.str();
};
template<class T> std::string vec_to_string(const std::vector<T> &a) {
return vec_to_string(a, stdFmt);
};
/// Templates for turning numbers (0D-matrices) into strings
template<class T> std::string vec_to_string(const T &a, const char *fmt) {
std::vector<T> vec;
vec.push_back(a);
return vec_to_string(vec, fmt);
};
template<class T> std::string vec_to_string(const T &a) {
return vec_to_string(a, stdFmt);
};
///Templates for turning 2D-matrices into strings
template<class T> std::string vec_to_string(const std::vector<std::vector<T> > &A, const char *fmt) {
if (A.size()<1) return std::string("");
std::stringstream out;
out << "[ " << vec_to_string(A[0], fmt);
for (size_t j = 1; j < A.size(); j++) {
out << ", " << std::endl << " " << vec_to_string(A[j], fmt);
}
out << " ]";
return out.str();
};
template<class T> std::string vec_to_string(const std::vector<std::vector<T> > &A) {
return vec_to_string(A, stdFmt);
};
///Templates for turning Eigen matrices into strings
template <class T> std::string mat_to_string(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &A, const char *fmt) {
//std::string mat_to_string(const Eigen::MatrixXd &A, const char *fmt) {
std::size_t r = A.rows();
std::size_t c = A.cols();
if ((r<1)||(c<1)) return std::string("");
std::stringstream out;
out << "[ ";
if (r==1) {
out << format(fmt, A(0,0));
for (size_t j = 1; j < c; j++) {
out << ", " << format(fmt, A(0,j));
}
} else {
out << mat_to_string(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>(A.row(0)), fmt);
for (size_t i = 1; i < r; i++) {
out << ", " << std::endl << " " << mat_to_string(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>(A.row(i)), fmt);
}
}
out << " ]";
return out.str();
};
template <class T> std::string mat_to_string(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &A) {
//std::string vec_to_string(const Eigen::MatrixXd &A) {
return mat_to_string(A, stdFmt);
};
/// Template class for turning numbers (0D-matrices) into strings
//template<class T> std::string vec_to_string(const T &a){
// return vec_to_string(a, stdFmt);
// std::stringstream out;
// out << format("[ %7.3f ]",a);
// return out.str();
//};
//template<class T> std::string vec_to_string(const VectorNd<0, T> &a){
// return vec_to_string(a, stdFmt);
//};
//template<class T> std::string vec_to_string(const VectorNd<0, T> &a, const char *fmt) {
// VectorNd<1, T> vec;
// vec.push_back(a);
// return vec_to_string(vec, fmt);
//};
//
/////Template classes for turning vectors (1D-matrices) into strings
//template<class T> std::string vec_to_string(const VectorNd<1, T> &a) {
// return vec_to_string(a, stdFmt);
//};
//template<class T> std::string vec_to_string(const VectorNd<1, T> &a, const char *fmt) {
// if (a.size()<1) {
// return std::string("");
// } else {
// std::stringstream out;
// out << "[ ";
// out << format(fmt,a[0]);
// for (size_t j = 1; j < a.size(); j++) {
// out << ", ";
// out << format(fmt,a[j]);
// }
// out << " ]";
// return out.str();
// }
//};
//
/////Template classes for turning 2D-matrices into strings
//template<class T> std::string vec_to_string(const VectorNd<2, T> &A) {
// return vec_to_string(A, stdFmt);
//}
//template<class T> std::string vec_to_string(const VectorNd<2, T> &A, const char *fmt) {
// if (A.size()<1) return std::string("");
// std::stringstream out;
// out << "[ " << format(fmt,A[0]);
// for (size_t j = 1; j < A.size(); j++) {
// out << ", " << std::endl << " " << vec_to_string(A[j], fmt);
// }
// out << " ]";
// return out.str();
//}
///// Publish the linear algebra solver
//template<class T> std::vector<T> linsolve(std::vector<std::vector<T> > const& A, std::vector<T> const& b);
//template<class T> std::vector<std::vector<T> > linsolve(std::vector<std::vector<T> > const& A, std::vector<std::vector<T> > const& B);
@@ -43,10 +363,20 @@ namespace CoolProp{
//template<class T> std::string vec_to_string( std::vector<T> const& a, const char *fmt);
//template<class T> std::string vec_to_string(std::vector<std::vector<T> > const& A, const char *fmt);
// Forward definitions
template<class T> std::size_t num_rows (std::vector<std::vector<T> > const& in);
template<class T> std::size_t max_cols (std::vector<std::vector<T> > const& in);
/// Some shortcuts and regularly needed operations
template<class T> std::size_t num_rows ( std::vector<T> const& in){ return in.size(); }
template<class T> std::size_t num_rows (std::vector<std::vector<T> > const& in){ return in.size(); }
template<class T> std::size_t max_cols (std::vector<std::vector<T> > const& in){
std::size_t cols = 0;
std::size_t col = 0;
for (std::size_t i = 0; i < in.size(); i++) {
col = in[i].size();
if (cols<col) {cols = col;}
}
return cols;
};
template<class T> bool is_squared(std::vector<std::vector<T> > const& in){
std::size_t cols = max_cols(in);
if (cols!=num_rows(in)) { return false;}
@@ -57,17 +387,8 @@ template<class T> bool is_squared(std::vector<std::vector<T> > co
}
return true;
};
template<class T> std::size_t max_cols (std::vector<std::vector<T> > const& in){
std::size_t cols = 0;
std::size_t col = 0;
for (std::size_t i = 0; i < in.size(); i++) {
col = in[i].size();
if (cols<col) {cols = col;}
}
return cols;
};
/// Some shortcuts and regularly needed operations
template<class T> std::size_t num_rows (std::vector<std::vector<T> > const& in){ return in.size(); }
template<class T> std::size_t num_cols ( std::vector<T> const& in){ return 1; }
template<class T> std::size_t num_cols (std::vector<std::vector<T> > const& in){
if (num_rows(in)>0) {
if (is_squared(in)) {
@@ -398,42 +719,8 @@ template<class T> std::vector< std::vector<T> > invert(std::vector<std::vecto
return linsolve(in,identity);
};
template<class T> std::string vec_to_string( T const& a){
std::stringstream out;
out << format("[ %7.3f ]",a);
return out.str();
};
template<class T> std::string vec_to_string( std::vector<T> const& a) {
return vec_to_string(a,"%7.3g");
};
template<class T> std::string vec_to_string( std::vector<T> const& a, const char *fmt) {
if (a.size()<1) {
return std::string("");
} else {
std::stringstream out;
out << format("[ ");
out << format(fmt,a[0]);
for (size_t j = 1; j < a.size(); j++) {
out << ", ";
out << format(fmt,a[j]);
}
out << " ]";
return out.str();
}
};
template<class T> std::string vec_to_string(std::vector<std::vector<T> > const& A) {
return vec_to_string(A, "%7.3g");
}
template<class T> std::string vec_to_string(std::vector<std::vector<T> > const& A, const char *fmt) {
std::stringstream out;
for (size_t j = 0; j < A.size(); j++) {
out << vec_to_string(A[j], fmt);
}
return out.str();
}
}; /* namespace CoolProp */
#endif

1251
include/PolyMath.h
View File

File diff suppressed because it is too large Load Diff

68
include/rapidjson/rapidjson_include.h
View File

@@ -100,6 +100,46 @@ namespace cpjson
return out;
};
/// A convenience function to get a 2D double array compactly
inline std::vector< std::vector<double> > get_double_array2D(rapidjson::Value &v)
{
std::vector< std::vector<double> > out;
std::vector<double> tmp;
if (!v.IsArray()) { throw CoolProp::ValueError("input is not an array"); }
for (rapidjson::Value::ValueIterator itr = v.Begin(); itr != v.End(); ++itr)
{
if (!itr->IsArray()) { throw CoolProp::ValueError("input is not a 2D array"); }
tmp.clear();
for (rapidjson::Value::ValueIterator i = itr->Begin(); i != itr->End(); ++i)
{
if (!i->IsNumber()){throw CoolProp::ValueError("input is not a number");}
tmp.push_back(i->GetDouble());
}
out.push_back(tmp);
}
return out;
};
/// A convenience function to get a 2D long double array compactly
inline std::vector< std::vector<long double> > get_long_double_array2D(rapidjson::Value &v)
{
std::vector< std::vector<long double> > out;
std::vector<long double> tmp;
if (!v.IsArray()) { throw CoolProp::ValueError("input is not an array"); }
for (rapidjson::Value::ValueIterator itr = v.Begin(); itr != v.End(); ++itr)
{
if (!itr->IsArray()) { throw CoolProp::ValueError("input is not a 2D array"); }
tmp.clear();
for (rapidjson::Value::ValueIterator i = itr->Begin(); i != itr->End(); ++i)
{
if (!i->IsNumber()){throw CoolProp::ValueError("input is not a number");}
tmp.push_back(i->GetDouble());
}
out.push_back(tmp);
}
return out;
};
/// A convenience function to get a long double array compactly
inline std::vector<long double> get_long_double_array(rapidjson::Value &v, std::string name)
{
@@ -135,6 +175,34 @@ namespace cpjson
return buffer.GetString();
};
// /// A convenience function to set an array compactly
// template<typename T>
// inline void set_array(const char *key, const std::vector<T> &vec, rapidjson::Value &value, rapidjson::Document &doc)
// {
// rapidjson::Value _v(rapidjson::kArrayType);
// for (unsigned int i = 0; i < vec.size(); ++i)
// {
// _v.PushBack(vec[i],doc.GetAllocator());
// }
// value.AddMember(key, _v, doc.GetAllocator());
// };
//
// /// A convenience function to set an array compactly
// template<typename T>
// inline void set_array2D(const char *key, const std::vector< std::vector<T> > &vec, rapidjson::Value &value, rapidjson::Document &doc)
// {
// rapidjson::Value _i(rapidjson::kArrayType);
// rapidjson::Value _j;
// for (unsigned int i = 0; i < vec.size(); ++i) {
// _j = rapidjson::Value(rapidjson::kArrayType);
// for (unsigned int j = 0; j < vec[i].size(); ++j) {
// _j.PushBack(vec[j],doc.GetAllocator());
// }
// _i.PushBack(_j,doc.GetAllocator());
// }
// value.AddMember(key, _i, doc.GetAllocator());
// };
/// A convenience function to set a double array compactly
inline void set_double_array(const char *key, const std::vector<double> &vec, rapidjson::Value &value, rapidjson::Document &doc)
{

130
src/Backends/Helmholtz/Fluids/Incompressible.cpp
View File

@@ -1,130 +0,0 @@
/*
* Incompressible.cpp
*
* Created on: 20 Dec 2013
* Author: jowr
*/
#include "Incompressible.h"
namespace CoolProp {
std::vector<std::vector<double> > Incompressible::changeAxis(const std::vector<double> &input){
std::vector<std::vector<double> > output;
std::vector<double> tmp;
size_t sizeX = input.size();
for (size_t i = 0; i < sizeX; ++i){
tmp.clear();
tmp.push_back(input[i]);
output.push_back(tmp);
}
return output;
}
/* All functions need T and p as input. Might not
* be necessary, but gives a clearer structure.
*/
/// Density as a function of temperature, pressure and composition.
double Incompressible::rho(double T_K, double p) {
return poly.polyval(cRho, getxInput(x), getTInput(T_K));
}
/// Heat capacities as a function of temperature, pressure and composition.
double Incompressible::c(double T_K, double p) {
return poly.polyval(cHeat, getxInput(x), getTInput(T_K));
}
/// Enthalpy as a function of temperature, pressure and composition.
double Incompressible::h(double T_K, double p) {
return h_u(T_K, p);
}
/// Entropy as a function of temperature, pressure and composition.
double Incompressible::s(double T_K, double p) {
return poly.polyfracintcentral(cHeat, getxInput(x), T_K, Tbase)
- poly.polyfracintcentral(cHeat, getxInput(x), Tref, Tbase);
}
/// Viscosity as a function of temperature, pressure and composition.
double Incompressible::visc(double T_K, double p) {
return expo.expval(cVisc, getxInput(x), getTInput(T_K), 2) / 1e3;
}
/// Thermal conductivity as a function of temperature, pressure and composition.
double Incompressible::cond(double T_K, double p) {
return poly.polyval(cCond, getxInput(x), getTInput(T_K));
}
/// Internal energy as a function of temperature, pressure and composition.
double Incompressible::u(double T_K, double p) {
return poly.polyint(cHeat, getxInput(x), getTInput(T_K))
- poly.polyint(cHeat, getxInput(x), getTInput(Tref));
}
/// Saturation pressure as a function of temperature and composition.
double Incompressible::psat(double T_K ){throw NotImplementedError("Psat is not available");};
/// Freezing temperature as a function of pressure and composition.
double Incompressible::Tfreeze( double p){throw NotImplementedError("Tfreeze is not available");};
/*
* Some more functions to provide a single implementation
* of important routines.
* We start with the check functions that can validate input
* in terms of pressure p, temperature T and composition x.
*/
/// Check validity of temperature input.
/** Compares the given temperature T to the result of a
* freezing point calculation. This is not necessarily
* defined for all fluids, default values do not
* cause errors. */
bool Incompressible::checkT(double T_K, double p){
if( Tmin < 0. ) {
throw ValueError("Please specify the minimum temperature.");
} else if( Tmax < 0.) {
throw ValueError("Please specify the maximum temperature.");
} else if ( (Tmin>T_K) || (T_K>Tmax) ) {
throw ValueError(format("Your temperature %f is not between %f and %f.",T_K,Tmin,Tmax));
} else if (T_K < Tfreeze(p)) {
throw ValueError(format("Your temperature %f is below the freezing point of %f.",T_K,Tfreeze(p)));
} else {
return true;
}
return false;
}
/// Check validity of pressure input.
/** Compares the given pressure p to the saturation pressure at
* temperature T and throws and exception if p is lower than
* the saturation conditions.
* The default value for psat is -1 yielding true if psat
* is not redefined in the subclass.
* */
bool Incompressible::checkP(double T_K, double p) {
double ps = psat(T_K);
if (p<ps) {
throw ValueError(format("Equations are valid for solution phase only: %f < %f. ",p,ps));
} else {
return true;
}
}
/// Check validity of composition input.
/** Compares the given composition x to a stored minimum and
* maximum value. Enforces the redefinition of xmin and
* xmax since the default values cause an error. */
bool Incompressible::checkX(double x){
if( xmin < 0. ) {
throw ValueError("Please specify the minimum concentration.");
} else if( xmax < 0.) {
throw ValueError("Please specify the maximum concentration.");
} else if ( (xmin>x) || (x>xmax) ) {
throw ValueError(format("Your composition %f is not between %f and %f.",x,xmin,xmax));
} else {
return true;
}
return false;
}
/// Check validity of temperature, pressure and composition input.
bool Incompressible::checkTPX(double T, double p, double x) {
return (checkT(T,p) && checkP(T,p) && checkX(x));
}
} /* namespace CoolProp */

177
src/Backends/Helmholtz/Fluids/Incompressible.h
View File

@@ -1,177 +0,0 @@
/*
* Incompressible.h
*
* Created on: 20 Dec 2013
* Author: jowr
*/
#ifndef INCOMPRESSIBLE_H_
#define INCOMPRESSIBLE_H_
#include "PolyMath.h"
namespace CoolProp {
class Incompressible{
protected:
std::string name;
std::string description;
std::string reference;
double Tmin, Tmax;
double xmin, xmax;
double TminPsat;
double xref, Tref;
double xbase, Tbase;
double x;
std::vector< std::vector<double> > cRho;
std::vector< std::vector<double> > cHeat;
std::vector< std::vector<double> > cVisc;
std::vector< std::vector<double> > cCond;
std::vector< std::vector<double> > cPsat;
std::vector<double> cTfreeze;
std::vector<std::vector<double> > changeAxis(const std::vector<double> &input);
BasePolynomial poly;
PolynomialSolver solver;
BaseExponential expo;
public:
Incompressible();
virtual ~Incompressible();
std::string getName() const {return name;}
std::string get_name() const {return getName();}// For backwards-compatibility.
std::string getDescription() const {return description;}
std::string getReference() const {return reference;}
double getTmax() const {return Tmax;}
double getTmin() const {return Tmin;}
double getxmax() const {return xmax;}
double getxmin() const {return xmin;}
double getTminPsat() const {return TminPsat;}
double getTref() const {return Tref;}
double getxref() const {return xref;}
double getTbase() const {return Tbase;}
double getxbase() const {return xbase;}
double getx() const {return x;}
void setName(std::string name) {this->name = name;}
void setDescription(std::string description) {this->description = description;}
void setReference(std::string reference) {this->reference = reference;}
void setTmax(double Tmax) {this->Tmax = Tmax;}
void setTmin(double Tmin) {this->Tmin = Tmin;}
void setxmax(double xmax) {this->xmax = xmax;}
void setxmin(double xmin) {this->xmin = xmin;}
void setTminPsat(double TminPsat) {this->TminPsat = TminPsat;}
void setTref(double Tref) {this->Tref = Tref;}
void setxref(double xref) {this->xref = xref;}
void setTbase(double Tbase) {this->Tbase = Tbase;}
void setxbase(double xbase) {this->xbase = xbase;}
void setx(double x) {this->x = x;}
// Setters for the coefficients
void setcRho(std::vector<std::vector<double> > cRho){this->cRho = cRho;}
void setcHeat(std::vector<std::vector<double> > cHeat){this->cHeat = cHeat;}
void setcVisc(std::vector<std::vector<double> > cVisc){this->cVisc = cVisc;}
void setcCond(std::vector<std::vector<double> > cCond){this->cCond = cCond;}
void setcPsat(std::vector<std::vector<double> > cPsat){this->cPsat = cPsat;}
void setcTfreeze(std::vector<double> cTfreeze){this->cTfreeze = cTfreeze;}
void setcRho(std::vector<double> cRho){this->cRho = changeAxis(cRho);}
void setcHeat(std::vector<double> cHeat){this->cHeat = changeAxis(cHeat);}
void setcVisc(std::vector<double> cVisc){this->cVisc = changeAxis(cVisc);}
void setcCond(std::vector<double> cCond){this->cCond = changeAxis(cCond);}
void setcPsat(std::vector<double> cPsat){this->cPsat = changeAxis(cPsat);}
double getTInput(double curTValue){return curTValue-Tbase;}
double getxInput(double curxValue){return (curxValue-xbase)*100.0;}
/* All functions need T and p as input. Might not
* be necessary, but gives a clearer structure.
*/
/// Density as a function of temperature, pressure and composition.
virtual double rho (double T_K, double p);
/// Heat capacities as a function of temperature, pressure and composition.
virtual double c (double T_K, double p);
virtual double cp (double T_K, double p){return c(T_K,p);};
virtual double cv (double T_K, double p){return c(T_K,p);};
/// Entropy as a function of temperature, pressure and composition.
virtual double s (double T_K, double p);
/// Internal energy as a function of temperature, pressure and composition.
virtual double u (double T_K, double p);
/// Enthalpy as a function of temperature, pressure and composition.
virtual double h (double T_K, double p);
/// Viscosity as a function of temperature, pressure and composition.
virtual double visc(double T_K, double p);
/// Thermal conductivity as a function of temperature, pressure and composition.
virtual double cond(double T_K, double p);
/// Saturation pressure as a function of temperature and composition.
virtual double psat(double T_K );
/// Freezing temperature as a function of pressure and composition.
virtual double Tfreeze( double p);
protected:
/* Define internal energy and enthalpy as functions of the
* other properties to provide data in case there are no
* coefficients.
*/
/// Enthalpy from u, p and rho.
/** Calculate enthalpy as a function of temperature and
* pressure employing functions for internal energy and
* density. Provides consistent formulations. */
double h_u(double T_K, double p) {
return u(T_K,p)+p/rho(T_K,p);
};
/// Internal energy from h, p and rho.
/** Calculate internal energy as a function of temperature
* and pressure employing functions for enthalpy and
* density. Provides consistent formulations. */
double u_h(double T_K, double p) {
return h(T_K,p)-p/rho(T_K,p);
};
/*
* Some more functions to provide a single implementation
* of important routines.
* We start with the check functions that can validate input
* in terms of pressure p, temperature T and composition x.
*/
/// Check validity of temperature input.
/** Compares the given temperature T to the result of a
* freezing point calculation. This is not necessarily
* defined for all fluids, default values do not cause errors. */
bool checkT(double T_K, double p);
/// Check validity of pressure input.
/** Compares the given pressure p to the saturation pressure at
* temperature T and throws and exception if p is lower than
* the saturation conditions.
* The default value for psat is -1 yielding true if psat
* is not redefined in the subclass.
* */
bool checkP(double T_K, double p);
/// Check validity of composition input.
/** Compares the given composition x to a stored minimum and
* maximum value. Enforces the redefinition of xmin and
* xmax since the default values cause an error. */
bool checkX(double x);
/// Check validity of temperature, pressure and composition input.
bool checkTPX(double T, double p, double x);
};
} /* namespace CoolProp */
#endif /* INCOMPRESSIBLE_H_ */

179
src/Backends/Incompressible/IncompressibleLibrary.cpp
View File

@@ -1,7 +1,186 @@
#include "IncompressibleLibrary.h"
#include "MatrixMath.h"
#include "rapidjson/rapidjson_include.h"
#include "all_incompressibles_JSON.h" // Makes a std::string variable called all_fluids_JSON
namespace CoolProp{
/// A general function to parse the json files that hold the coefficient matrices
IncompressibleData JSONIncompressibleLibrary::parse_coefficients(rapidjson::Value &obj, std::string id, bool vital){
IncompressibleData fluidData;
if (obj.HasMember(id.c_str())) {
//rapidjson::Value value = obj[id.c_str()];
if (obj[id.c_str()].HasMember("type")){
if (obj[id.c_str()].HasMember("coeffs")){
std::string type = cpjson::get_string(obj[id.c_str()], "type");
if (!type.compare("polynomial")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array2D(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else if (!type.compare("exponential")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else if (!type.compare("exppolynomial")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array2D(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else{
throw ValueError(format("The type [%s] is not understood for [%s] of incompressible fluids. Please check your JSON file.", type.c_str(), id.c_str()));
}
}
else{
throw ValueError(format("Your file does not have an entry for \"coeffs\" in [%s], which is vital for this function.", id.c_str()));
}
}
else{
throw ValueError(format("Your file does not have an entry for \"type\" in [%s], which is vital for this function.", id.c_str()));
}
}
else{
if (vital) {
throw ValueError(format("Your file does not have information for [%s], which is vital for an incompressible fluid.", id.c_str()));
}
}
return fluidData;
}
/// Get a double from the JSON storage if it is defined, otherwise return def
double JSONIncompressibleLibrary::parse_value(rapidjson::Value &obj, std::string id, bool vital, double def=0.0){
if (obj.HasMember(id.c_str())) {return cpjson::get_double(obj, id);}
else{
if (vital) {
throw ValueError(format("Your file does not have information for [%s], which is vital for an incompressible fluid.", id.c_str()));
}
else{
return def;
}
}
}
/// Add all the fluid entries in the rapidjson::Value instance passed in
void JSONIncompressibleLibrary::add_many(rapidjson::Value &listing) {
for (rapidjson::Value::ValueIterator itr = listing.Begin();
itr != listing.End(); ++itr) {
add_one(*itr);
}
};
void JSONIncompressibleLibrary::add_one(rapidjson::Value &fluid_json) {
_is_empty = false;
// Get the next index for this fluid
std::size_t index = fluid_map.size();
// Add index->fluid mapping
fluid_map[index] = IncompressibleFluid();
// Create an instance of the fluid
IncompressibleFluid &fluid = fluid_map[index];
fluid.setName(cpjson::get_string(fluid_json, "name"));
fluid.setDescription(cpjson::get_string(fluid_json, "description"));
fluid.setReference(cpjson::get_string(fluid_json, "reference"));
fluid.setTmax(parse_value(fluid_json, "Tmax", true, 0.0));
fluid.setTmin(parse_value(fluid_json, "Tmin", true, 0.0));
fluid.setxmax(parse_value(fluid_json, "xmax", false, 1.0));
fluid.setxmin(parse_value(fluid_json, "xmin", false, 0.0));
fluid.setTminPsat(parse_value(fluid_json, "TminPsat", false, 0.0));
fluid.setTbase(parse_value(fluid_json, "Tbase", false, 0.0));
fluid.setxbase(parse_value(fluid_json, "xbase", false, 0.0));
/// Setters for the coefficients
fluid.setDensity(parse_coefficients(fluid_json, "density", true));
fluid.setSpecificHeat(parse_coefficients(fluid_json, "specific_heat", true));
fluid.setViscosity(parse_coefficients(fluid_json, "viscosity", false));
fluid.setConductivity(parse_coefficients(fluid_json, "conductivity", false));
fluid.setPsat(parse_coefficients(fluid_json, "saturation_pressure", false));
fluid.setTfreeze(parse_coefficients(fluid_json, "T_freeze", false));
fluid.setVolToMass(parse_coefficients(fluid_json, "volume2mass", false));
fluid.setMassToMole(parse_coefficients(fluid_json, "mass2mole", false));
fluid.set_reference_state(
parse_value(fluid_json, "Tref", false, 25+273.15) ,
parse_value(fluid_json, "pref", false, 1.01325e5) ,
parse_value(fluid_json, "xref", false, 0.0) ,
parse_value(fluid_json, "href", false, 0.0) ,
parse_value(fluid_json, "sref", false, 0.0)
);
/// A function to check coefficients and equation types.
fluid.validate();
// Add name->index mapping
string_to_index_map[fluid.getName()] = index;
};
/// Get an IncompressibleFluid instance stored in this library
/**
@param name Name of the fluid
*/
IncompressibleFluid& JSONIncompressibleLibrary::get(std::string key) {
std::map<std::string, std::size_t>::iterator it;
// Try to find it
it = string_to_index_map.find(key);
// If it is found
if (it != string_to_index_map.end()) {
return get(it->second);
} else {
throw ValueError(
format(
"key [%s] was not found in string_to_index_map in JSONIncompressibleLibrary",
key.c_str()
)
);
}
};
/// Get a CoolPropFluid instance stored in this library
/**
@param key The index of the fluid in the map
*/
IncompressibleFluid& JSONIncompressibleLibrary::get(std::size_t key) {
std::map<std::size_t, IncompressibleFluid>::iterator it;
// Try to find it
it = fluid_map.find(key);
// If it is found
if (it != fluid_map.end()) {
return it->second;
} else {
throw ValueError(
format("key [%d] was not found in JSONIncompressibleLibrary",key));
}
};
static JSONIncompressibleLibrary library;

162
src/Backends/Incompressible/IncompressibleLibrary.h
View File

@@ -22,172 +22,44 @@ a rapidjson array of fluids to the add_many function.
*/
class JSONIncompressibleLibrary
{
/// Map from CAS code to JSON instance. For pseudo-pure fluids, use name in place of CAS code since no CASE number is defined for mixtures
/// Map from CAS code to JSON instance.
/** This is not practical for the incomressibles, the CAS may not be
* defined for blends of heat transfer fluids and solutions.
*/
std::map<std::size_t, IncompressibleFluid> fluid_map;
std::vector<std::string> name_vector;
std::map<std::string, std::size_t> string_to_index_map;
bool _is_empty;
protected:
/// A general function to parse the json files that hold the coefficient matrices
IncompressibleData parse_coefficients(rapidjson::Value &obj, std::string id, bool vital);
double parse_value(rapidjson::Value &obj, std::string id, bool vital, double def);
/// Parse the viscosity
void parse_viscosity(rapidjson::Value &viscosity, IncompressibleFluid & fluid)
{
if (viscosity.HasMember("type")){
std::string type = cpjson::get_string(viscosity, "type");
if (!type.compare("polynomial")){
fluid.viscosity.type = CoolProp::IncompressibleViscosityVariables::INCOMPRESSIBLE_VISCOSITY_POLYNOMIAL;
fluid.viscosity.poly.coeffs = cpjson::get_double_array(viscosity["coeffs"]);
return;
}
else{
throw ValueError(format("viscosity type [%s] is not understood for fluid %s", type.c_str(), fluid.name.c_str()));
}
}
else{
throw ValueError(format("viscosity does not have \"type\" for fluid %s", fluid.name.c_str()));
}
};
/// Parse the conductivity
void parse_conductivity(rapidjson::Value &conductivity, IncompressibleFluid & fluid)
{
if (conductivity.HasMember("type")){
std::string type = cpjson::get_string(conductivity, "type");
if (!type.compare("polynomial")){
fluid.conductivity.type = CoolProp::IncompressibleConductivityVariables::INCOMPRESSIBLE_CONDUCTIVITY_POLYNOMIAL;
fluid.conductivity.poly.coeffs = cpjson::get_double_array(conductivity["coeffs"]);
return;
}
else{
throw ValueError(format("conductivity type [%s] is not understood for fluid %s", type.c_str(), fluid.name.c_str()));
}
}
else{
throw ValueError(format("conductivity does not have \"type\" for fluid %s", fluid.name.c_str()));
}
};
/// Parse the specific_heat
void parse_specific_heat(rapidjson::Value &specific_heat, IncompressibleFluid & fluid)
{
if (specific_heat.HasMember("type")){
std::string type = cpjson::get_string(specific_heat, "type");
if (!type.compare("polynomial")){
fluid.specific_heat.type = CoolProp::IncompressibleSpecificHeatVariables::INCOMPRESSIBLE_SPECIFIC_HEAT_POLYNOMIAL; return;
fluid.specific_heat.poly.coeffs = cpjson::get_double_array(specific_heat["coeffs"]);
}
else{
throw ValueError(format("specific_heat type [%s] is not understood for fluid %s", type.c_str(), fluid.name.c_str()));
}
}
else{
throw ValueError(format("specific_heat does not have \"type\" for fluid %s", fluid.name.c_str()));
}
};
/// Parse the density
void parse_density(rapidjson::Value &density, IncompressibleFluid & fluid)
{
if (density.HasMember("type")){
std::string type = cpjson::get_string(density, "type");
if (!type.compare("polynomial")){
fluid.density.type = CoolProp::IncompressibleDensityVariables::INCOMPRESSIBLE_DENSITY_POLYNOMIAL; return;
fluid.density.poly.coeffs = cpjson::get_double_array(density["coeffs"]);
}
else{
throw ValueError(format("density type [%s] is not understood for fluid %s", type.c_str(), fluid.name.c_str()));
}
}
else{
throw ValueError(format("density does not have \"type\" for fluid %s", fluid.name.c_str()));
}
};
/// Validate the fluid file that was just constructed
void validate(IncompressibleFluid & fluid)
{
}
public:
// Default constructor;
JSONIncompressibleLibrary(){
_is_empty = true;
};
JSONIncompressibleLibrary(){ _is_empty = true;};
bool is_empty(void){ return _is_empty;};
/// Add all the fluid entries in the rapidjson::Value instance passed in
void add_many(rapidjson::Value &listing)
{
for (rapidjson::Value::ValueIterator itr = listing.Begin(); itr != listing.End(); ++itr)
{
add_one(*itr);
}
};
void add_one(rapidjson::Value &fluid_json)
{
_is_empty = false;
void add_many(rapidjson::Value &listing);
void add_one(rapidjson::Value &fluid_json);
// Get the next index for this fluid
std::size_t index = fluid_map.size();
// Add index->fluid mapping
fluid_map[index] = IncompressibleFluid();
// Create an instance of the fluid
IncompressibleFluid &fluid = fluid_map[index];
fluid.name = cpjson::get_string(fluid_json, "name");
fluid.Tmin = cpjson::get_double(fluid_json, "Tmin");
fluid.Tmax = cpjson::get_double(fluid_json, "Tmax");
parse_conductivity(fluid_json["conductivity"], fluid);
parse_density(fluid_json["density"], fluid);
parse_viscosity(fluid_json["viscosity"], fluid);
parse_specific_heat(fluid_json["specific_heat"], fluid);
// Add name->index mapping
string_to_index_map[fluid.name] = index;
};
/// Get an IncompressibleFluid instance stored in this library
/**
@param name Name of the fluid
*/
IncompressibleFluid& get(std::string key)
{
std::map<std::string, std::size_t>::iterator it;
// Try to find it
it = string_to_index_map.find(key);
// If it is found
if (it != string_to_index_map.end()){
return get(it->second);
}
else{
throw ValueError(format("key [%s] was not found in string_to_index_map in JSONIncompressibleLibrary",key.c_str()));
}
};
IncompressibleFluid& get(std::string key);
/// Get a CoolPropFluid instance stored in this library
/**
@param key The index of the fluid in the map
*/
IncompressibleFluid& get(std::size_t key)
{
std::map<std::size_t, IncompressibleFluid>::iterator it;
// Try to find it
it = fluid_map.find(key);
// If it is found
if (it != fluid_map.end()){
return it->second;
}
else{
throw ValueError(format("key [%d] was not found in JSONIncompressibleLibrary",key));
}
};
IncompressibleFluid& get(std::size_t key);
/// Return a comma-separated list of fluid names
std::string get_fluid_list(void)
{
return strjoin(name_vector, ",");
};
std::string get_fluid_list(void){ return strjoin(name_vector, ",");};
};
/// Get a reference to the library instance

693
src/IncompressibleFluid.cpp Normal file
View File

@@ -0,0 +1,693 @@
#include "IncompressibleFluid.h"
#include "math.h"
#include "MatrixMath.h"
#include "PolyMath.h"
namespace CoolProp {
/// A thermophysical property provider for critical and reducing values as well as derivatives of Helmholtz energy
/**
This fluid instance is populated using an entry from a JSON file
*/
//IncompressibleFluid::IncompressibleFluid();
void IncompressibleFluid::set_reference_state(double T0, double p0, double x0=0.0, double h0=0.0, double s0=0.0){
this->rhoref = rho(T0,p0,x0);
this->pref = p0;
this->uref = h0 - p0/rhoref;
this->uref = u(T0,p0,x0);
this->href = h0; // set new reference value
this->sref = s0; // set new reference value
this->href = h(T0,p0,x0); // adjust offset to fit to equations
this->sref = s(T0,p0,x0); // adjust offset to fit to equations
}
void IncompressibleFluid::validate(){throw NotImplementedError("TODO");}
/// Base function that handles the custom data type
double IncompressibleFluid::baseFunction(IncompressibleData data, double T_in, double x_in=0.0){
switch (data.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
//throw NotImplementedError("Here you should implement the polynomial.");
return poly.evaluate(data.coeffs, T_in, x_in, 0, 0, Tbase, xbase);
break;
case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
//throw NotImplementedError("Here you should implement the exponential.");
poly.checkCoefficients(data.coeffs, 3,1);
return exp( data.coeffs(0,0) / ( (T_in-Tbase)+data.coeffs(1,0) ) - data.coeffs(2,0) );
break;
case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
//throw NotImplementedError("Here you should implement the exponential polynomial.");
return exp(poly.evaluate(data.coeffs, T_in, x_in, 0, 0, Tbase, xbase));
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,data.type));
break;
default:
throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.",__FILE__,__LINE__,data.type));
break;
}
return -_HUGE;
}
/// Entropy as a function of temperature, pressure and composition.
double IncompressibleFluid::s (double T, double p, double x=0.0){
IncompressibleData data = specific_heat;
switch (data.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
//throw NotImplementedError("Here you should implement the polynomial.");
return poly.integral(data.coeffs, T, x, 0, -1, 0, Tbase, xbase) - sref;
break;
case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
throw NotImplementedError(format("%s (%d): There is no automatic integration of the exponential function.",__FILE__,__LINE__));
break;
case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
throw NotImplementedError(format("%s (%d): There is no automatic integration of the exponential polynomial function.",__FILE__,__LINE__));
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,data.type));
break;
default:
throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.",__FILE__,__LINE__,data.type));
break;
}
return -_HUGE;
}
/// Internal energy as a function of temperature, pressure and composition.
double IncompressibleFluid::u (double T, double p, double x=0.0){
IncompressibleData data = specific_heat;
switch (data.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
//throw NotImplementedError("Here you should implement the polynomial.");
return poly.integral(data.coeffs, T, x, 0, 0, 0, Tbase, xbase) - uref;
break;
case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
throw NotImplementedError(format("%s (%d): There is no automatic integration of the exponential function.",__FILE__,__LINE__));
break;
case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
throw NotImplementedError(format("%s (%d): There is no automatic integration of the exponential polynomial function.",__FILE__,__LINE__));
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,data.type));
break;
default:
throw ValueError(format("%s (%d): Your function type \"[%d]\" is unknown.",__FILE__,__LINE__,data.type));
break;
}
return -_HUGE;
}
/// Freezing temperature as a function of pressure and composition.
double IncompressibleFluid::Tfreeze( double p, double x){
//IncompressibleClass::checkCoefficients(cTfreeze,5);
//std::vector<double> tmpVector(cTfreeze.begin()+1,cTfreeze.end());
//return polyval(tmpVector, x*100.0-cTfreeze[0])+273.15;
//Eigen::MatrixXd tmp = T_freeze.coeffs;
//removeRow(tmp, 0);
//double x_in = x-T_freeze.coeffs(0,0);
return poly.evaluate(Eigen::MatrixXd(T_freeze.coeffs.block(1,0,T_freeze.coeffs.rows()-1,1)).transpose(), x, 0, T_freeze.coeffs(0,0));
}
/// Define freezing point calculations
//double Tfreeze(double p, double x){
// IncompressibleClass::checkCoefficients(cTfreeze,5);
// std::vector<double> tmpVector(cTfreeze.begin()+1,cTfreeze.end());
// return polyval(tmpVector, x*100.0-cTfreeze[0])+273.15;
//}
/// Conversion from volume-based to mass-based composition.
double V2M (double T, double y){throw NotImplementedError("TODO");}
/// Conversion from mass-based to mole-based composition.
double M2M (double T, double x){throw NotImplementedError("TODO");}
/*
* Some more functions to provide a single implementation
* of important routines.
* We start with the check functions that can validate input
* in terms of pressure p, temperature T and composition x.
*/
/// Check validity of temperature input.
/** Compares the given temperature T to the result of a
* freezing point calculation. This is not necessarily
* defined for all fluids, default values do not cause errors. */
bool IncompressibleFluid::checkT(double T, double p, double x=0.0){throw NotImplementedError("TODO");}
/// Check validity of pressure input.
/** Compares the given pressure p to the saturation pressure at
* temperature T and throws and exception if p is lower than
* the saturation conditions.
* The default value for psat is -1 yielding true if psat
* is not redefined in the subclass.
* */
bool IncompressibleFluid::checkP(double T, double p, double x=0.0){throw NotImplementedError("TODO");}
/// Check validity of composition input.
/** Compares the given composition x to a stored minimum and
* maximum value. Enforces the redefinition of xmin and
* xmax since the default values cause an error. */
bool IncompressibleFluid::checkX(double x){throw NotImplementedError("TODO");}
/// Check validity of temperature, pressure and composition input.
bool IncompressibleFluid::checkTPX(double T, double p, double x=0.0){throw NotImplementedError("TODO");}
} /* namespace CoolProp */
// Testing still needs to be enhanced.
/* Below, I try to carry out some basic tests for both 2D and 1D
* polynomials as well as the exponential functions for vapour
* pressure etc.
*/
#ifdef ENABLE_CATCH
#include <math.h>
#include <iostream>
#include "catch.hpp"
Eigen::MatrixXd makeMatrix(const std::vector<double> &coefficients){
//IncompressibleClass::checkCoefficients(coefficients,18);
std::vector< std::vector<double> > matrix;
std::vector<double> tmpVector;
tmpVector.clear();
tmpVector.push_back(coefficients[0]);
tmpVector.push_back(coefficients[6]);
tmpVector.push_back(coefficients[11]);
tmpVector.push_back(coefficients[15]);
matrix.push_back(tmpVector);
tmpVector.clear();
tmpVector.push_back(coefficients[1]);
tmpVector.push_back(coefficients[7]);
tmpVector.push_back(coefficients[12]);
tmpVector.push_back(coefficients[16]);
matrix.push_back(tmpVector);
tmpVector.clear();
tmpVector.push_back(coefficients[2]);
tmpVector.push_back(coefficients[8]);
tmpVector.push_back(coefficients[13]);
tmpVector.push_back(coefficients[17]);
matrix.push_back(tmpVector);
tmpVector.clear();
tmpVector.push_back(coefficients[3]);
tmpVector.push_back(coefficients[9]);
tmpVector.push_back(coefficients[14]);
tmpVector.push_back(0.0);
matrix.push_back(tmpVector);
tmpVector.clear();
tmpVector.push_back(coefficients[4]);
tmpVector.push_back(coefficients[10]);
tmpVector.push_back(0.0);
tmpVector.push_back(0.0);
matrix.push_back(tmpVector);
tmpVector.clear();
tmpVector.push_back(coefficients[5]);
tmpVector.push_back(0.0);
tmpVector.push_back(0.0);
tmpVector.push_back(0.0);
matrix.push_back(tmpVector);
tmpVector.clear();
return CoolProp::vec_to_eigen(matrix).transpose();
}
TEST_CASE("Internal consistency checks and example use cases for the incompressible fluids","[IncompressibleFluids]")
{
bool PRINT = false;
std::string tmpStr;
std::vector<double> tmpVector;
std::vector< std::vector<double> > tmpMatrix;
SECTION("Test case for \"SylthermXLT\" by Dow Chemicals") {
std::vector<double> cRho;
cRho.push_back(+1.1563685145E+03);
cRho.push_back(-1.0269048032E+00);
cRho.push_back(-9.3506079577E-07);
cRho.push_back(+1.0368116627E-09);
CoolProp::IncompressibleData density;
density.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
density.coeffs = CoolProp::vec_to_eigen(cRho);
std::vector<double> cHeat;
cHeat.push_back(+1.1562261074E+03);
cHeat.push_back(+2.0994549103E+00);
cHeat.push_back(+7.7175381057E-07);
cHeat.push_back(-3.7008444051E-20);
CoolProp::IncompressibleData specific_heat;
specific_heat.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
specific_heat.coeffs = CoolProp::vec_to_eigen(cHeat);
std::vector<double> cCond;
cCond.push_back(+1.6121957379E-01);
cCond.push_back(-1.3023781944E-04);
cCond.push_back(-1.4395238766E-07);
CoolProp::IncompressibleData conductivity;
conductivity.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
conductivity.coeffs = CoolProp::vec_to_eigen(cCond);
std::vector<double> cVisc;
cVisc.push_back(+1.0337654989E+03);
cVisc.push_back(-4.3322764383E+01);
cVisc.push_back(+1.0715062356E+01);
CoolProp::IncompressibleData viscosity;
viscosity.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL;
viscosity.coeffs = CoolProp::vec_to_eigen(cVisc);
CoolProp::IncompressibleFluid XLT;
XLT.setName("XLT");
XLT.setDescription("SylthermXLT");
XLT.setReference("Dow Chemicals data sheet");
XLT.setTmax(533.15);
XLT.setTmin(173.15);
XLT.setxmax(0.0);
XLT.setxmin(0.0);
XLT.setTminPsat(533.15);
XLT.setTbase(0.0);
XLT.setxbase(0.0);
/// Setters for the coefficients
XLT.setDensity(density);
XLT.setSpecificHeat(specific_heat);
XLT.setViscosity(viscosity);
XLT.setConductivity(conductivity);
//XLT.setPsat(parse_coefficients(fluid_json, "saturation_pressure", false));
//XLT.setTfreeze(parse_coefficients(fluid_json, "T_freeze", false));
//XLT.setVolToMass(parse_coefficients(fluid_json, "volume2mass", false));
//XLT.setMassToMole(parse_coefficients(fluid_json, "mass2mole", false));
//XLT.set_reference_state(25+273.15, 1.01325e5, 0.0, 0.0, 0.0);
double Tref = 25+273.15;
double pref = 0.0;
double xref = 0.0;
double href = 0.0;
double sref = 0.0;
XLT.set_reference_state(Tref, pref, xref, href, sref);
/// A function to check coefficients and equation types.
//XLT.validate();
// Prepare the results and compare them to the calculated values
double acc = 0.0001;
double T = 273.15+50;
double p = 10e5;
double val = 0;
double res = 0;
// Compare density
val = 824.4615702148608;
res = XLT.rho(T,p);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
// Compare cp
val = 1834.7455527670554;
res = XLT.c(T,p);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
// Compare s
val = 145.59157247249246;
res = XLT.s(T,p);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
val = 0.0;
res = XLT.s(Tref,pref);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( val==res );
}
// Compare u
val = 45212.407309106304;
res = XLT.u(T,p);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
val = href - pref/XLT.rho(Tref,pref);
res = XLT.u(Tref,pref);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( val==res );
}
// Compare h
val = 46425.32011926845;
res = XLT.h(T,p);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
val = 0.0;
res = XLT.h(Tref,pref);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( val==res );
}
// Compare v
val = 0.0008931435169681835;
res = XLT.visc(T,p);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
// Compare l
val = 0.10410086156049088;
res = XLT.cond(T,p);
{
CAPTURE(T);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
}
SECTION("Test case for Methanol from SecCool") {
tmpVector.clear();
tmpVector.push_back( 960.24665800);
tmpVector.push_back(-1.2903839100);
tmpVector.push_back(-0.0161042520);
tmpVector.push_back(-0.0001969888);
tmpVector.push_back( 1.131559E-05);
tmpVector.push_back( 9.181999E-08);
tmpVector.push_back(-0.4020348270);
tmpVector.push_back(-0.0162463989);
tmpVector.push_back( 0.0001623301);
tmpVector.push_back( 4.367343E-06);
tmpVector.push_back( 1.199000E-08);
tmpVector.push_back(-0.0025204776);
tmpVector.push_back( 0.0001101514);
tmpVector.push_back(-2.320217E-07);
tmpVector.push_back( 7.794999E-08);
tmpVector.push_back( 9.937483E-06);
tmpVector.push_back(-1.346886E-06);
tmpVector.push_back( 4.141999E-08);
CoolProp::IncompressibleData density;
density.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
density.coeffs = makeMatrix(tmpVector);
tmpVector.clear();
tmpVector.push_back( 3822.9712300);
tmpVector.push_back(-23.122409500);
tmpVector.push_back( 0.0678775826);
tmpVector.push_back( 0.0022413893);
tmpVector.push_back(-0.0003045332);
tmpVector.push_back(-4.758000E-06);
tmpVector.push_back( 2.3501449500);
tmpVector.push_back( 0.1788839410);
tmpVector.push_back( 0.0006828000);
tmpVector.push_back( 0.0002101166);
tmpVector.push_back(-9.812000E-06);
tmpVector.push_back(-0.0004724176);
tmpVector.push_back(-0.0003317949);
tmpVector.push_back( 0.0001002032);
tmpVector.push_back(-5.306000E-06);
tmpVector.push_back( 4.242194E-05);
tmpVector.push_back( 2.347190E-05);
tmpVector.push_back(-1.894000E-06);
CoolProp::IncompressibleData specific_heat;
specific_heat.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
specific_heat.coeffs = makeMatrix(tmpVector);
tmpVector.clear();
tmpVector.push_back( 0.4082066700);
tmpVector.push_back(-0.0039816870);
tmpVector.push_back( 1.583368E-05);
tmpVector.push_back(-3.552049E-07);
tmpVector.push_back(-9.884176E-10);
tmpVector.push_back( 4.460000E-10);
tmpVector.push_back( 0.0006629321);
tmpVector.push_back(-2.686475E-05);
tmpVector.push_back( 9.039150E-07);
tmpVector.push_back(-2.128257E-08);
tmpVector.push_back(-5.562000E-10);
tmpVector.push_back( 3.685975E-07);
tmpVector.push_back( 7.188416E-08);
tmpVector.push_back(-1.041773E-08);
tmpVector.push_back( 2.278001E-10);
tmpVector.push_back( 4.703395E-08);
tmpVector.push_back( 7.612361E-11);
tmpVector.push_back(-2.734000E-10);
CoolProp::IncompressibleData conductivity;
conductivity.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
conductivity.coeffs = makeMatrix(tmpVector);
tmpVector.clear();
tmpVector.push_back( 1.4725525500);
tmpVector.push_back( 0.0022218998);
tmpVector.push_back(-0.0004406139);
tmpVector.push_back( 6.047984E-06);
tmpVector.push_back(-1.954730E-07);
tmpVector.push_back(-2.372000E-09);
tmpVector.push_back(-0.0411841566);
tmpVector.push_back( 0.0001784479);
tmpVector.push_back(-3.564413E-06);
tmpVector.push_back( 4.064671E-08);
tmpVector.push_back( 1.915000E-08);
tmpVector.push_back( 0.0002572862);
tmpVector.push_back(-9.226343E-07);
tmpVector.push_back(-2.178577E-08);
tmpVector.push_back(-9.529999E-10);
tmpVector.push_back(-1.699844E-06);
tmpVector.push_back(-1.023552E-07);
tmpVector.push_back( 4.482000E-09);
CoolProp::IncompressibleData viscosity;
viscosity.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL;
viscosity.coeffs = makeMatrix(tmpVector);
tmpVector.clear();
tmpVector.push_back( 27.755555600); // reference concentration in per cent
tmpVector.push_back(-22.973221700);
tmpVector.push_back(-1.1040507200);
tmpVector.push_back(-0.0120762281);
tmpVector.push_back(-9.343458E-05);
CoolProp::IncompressibleData T_freeze;
T_freeze.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_NOT_SET;
T_freeze.coeffs = CoolProp::vec_to_eigen(tmpVector);
// After preparing the coefficients, we have to create the objects
CoolProp::IncompressibleFluid CH3OH;
CH3OH.setName("CH3OH");
CH3OH.setDescription("Methanol solution");
CH3OH.setReference("SecCool software");
CH3OH.setTmax( 20 + 273.15);
CH3OH.setTmin(-50 + 273.15);
CH3OH.setxmax(0.5*100);
CH3OH.setxmin(0.0);
CH3OH.setTminPsat( 20 + 273.15);
CH3OH.setTbase(-4.48 + 273.15);
CH3OH.setxbase(31.57 / 100.0 * 100);
/// Setters for the coefficients
CH3OH.setDensity(density);
CH3OH.setSpecificHeat(specific_heat);
CH3OH.setViscosity(viscosity);
CH3OH.setConductivity(conductivity);
//CH3OH.setPsat(saturation_pressure);
CH3OH.setTfreeze(T_freeze);
//CH3OH.setVolToMass(volume2mass);
//CH3OH.setMassToMole(mass2mole);
//XLT.set_reference_state(25+273.15, 1.01325e5, 0.0, 0.0, 0.0);
double Tref = 25+273.15;
double pref = 0.0;
double xref = 0.25*100;
double href = 0.0;
double sref = 0.0;
CH3OH.set_reference_state(Tref, pref, xref, href, sref);
/// A function to check coefficients and equation types.
//CH3OH.validate();
// Prepare the results and compare them to the calculated values
double acc = 0.0001;
double T = 273.15+10;
double p = 10e5;
double x = 0.25*100;
double val = 0;
double res = 0;
// Compare density
val = 963.2886528091547;
res = CH3OH.rho(T,p,x);
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
// Compare cp
val = 3993.9748117022423;
res = CH3OH.c(T,p,x);
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
// Compare s
val = -206.62646783739274;
res = CH3OH.s(T,p,x);
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
val = 0.0;
res = CH3OH.s(Tref,pref,xref);
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(val);
CAPTURE(res);
CHECK( val==res );
}
// Compare u
val = -60043.78429641827;
res = CH3OH.u(T,p,x);
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
val = href - pref/CH3OH.rho(Tref,pref,xref);
res = CH3OH.u(Tref,pref,xref);
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(val);
CAPTURE(res);
CHECK( val==res );
}
// Compare h
val = -59005.67386390795;
res = CH3OH.h(T,p,x);
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
val = 0.0;
res = CH3OH.h(Tref,pref,xref);
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(val);
CAPTURE(res);
CHECK( val==res );
}
// Compare v
val = 0.0023970245009602097;
res = CH3OH.visc(T,p,x)/1e3;
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
// Compare l
val = 0.44791148414693727;
res = CH3OH.cond(T,p,x);
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
// Compare Tfreeze
val = 253.1293105454671;
res = CH3OH.Tfreeze(p,x)+273.15;
{
CAPTURE(T);
CAPTURE(p);
CAPTURE(x);
CAPTURE(val);
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
}
}
#endif /* ENABLE_CATCH */

3667
src/PolyMath.cpp
View File

File diff suppressed because it is too large Load Diff