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

Added 2D function to JSON readers, fixed 2D polynomial fitting

Browse Source
This commit is contained in:
jowr
2014-06-29 18:03:37 +02:00
parent fb5c56a937
commit 00e9ddf291
12 changed files with 1324 additions and 374 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)