mirror of
https://github.com/CoolProp/CoolProp.git
synced 2026-04-23 03:00:17 -04:00
19a4875879
* autopep8 rule-groups E101,W1,W2,W3 * autopep8 with rule group E3 (blank lines) autopep8 --in-place --recursive --max-line-length=200 --exclude="externals" --select="E101,E3,W1,W2,W3" . * tabs and space W191 * autopep8 aggressive
558 lines
23 KiB
Python
558 lines
23 KiB
Python
import numpy as np
|
|
from CoolProp.CoolProp import Props
|
|
|
|
import matplotlib.pyplot as plt
|
|
import matplotlib.mlab as mlab
|
|
|
|
import scipy.optimize
|
|
import scipy.stats
|
|
import random
|
|
import h5py
|
|
from templates import *
|
|
|
|
indices = []
|
|
|
|
|
|
class TermLibrary():
|
|
"""
|
|
Build a term library using the coefficients from Wagner and Pruss (IAPWS95)
|
|
"""
|
|
|
|
def __init__(self):
|
|
L,D,T = [],[],[]
|
|
|
|
for i in range(1,6):
|
|
for j in range(-4,9):
|
|
T.append(float(j)/8.0)
|
|
D.append(float(i))
|
|
L.append(float(0))
|
|
for i in range(1,16):
|
|
for j in range(1,16):
|
|
T.append(float(j))
|
|
D.append(float(i))
|
|
L.append(float(1))
|
|
for i in range(1,13):
|
|
for j in range(1,11):
|
|
T.append(float(j))
|
|
D.append(float(i))
|
|
L.append(float(2))
|
|
for i in range(1,6):
|
|
for j in range(10,24):
|
|
T.append(float(j))
|
|
D.append(float(i))
|
|
L.append(float(3))
|
|
for i in range(1,10):
|
|
for j in range(10,21):
|
|
T.append(float(j))
|
|
D.append(float(i)*2)
|
|
L.append(float(4))
|
|
|
|
self.T = T
|
|
self.D = D
|
|
self.L = L
|
|
|
|
|
|
from Helmholtz import helmholtz
|
|
|
|
|
|
def rsquared(x, y):
|
|
""" Return R^2 where x and y are array-like."""
|
|
|
|
slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(x, y)
|
|
return r_value**2
|
|
|
|
|
|
def get_fluid_constants(Ref):
|
|
if Ref == 'R407F':
|
|
RefString = 'REFPROP-MIX:R32[0.47319469]&R125[0.2051091]&R134a[0.32169621]'
|
|
elif Ref == 'R410A':
|
|
RefString = 'REFPROP-MIX:R32[0.6976147]&R125[0.3023853]'
|
|
|
|
LIBRARY = TermLibrary()
|
|
|
|
# Coefficients for HFC blends
|
|
LIBRARY.N = np.array([9.87252E-01, -1.03017E+00, 1.17666E+00, 6.10984E+00, -7.79453E+00, 1.83377E-02, 1.05880E+00, -1.12018E+00, 6.29064E-01, 6.24982E+00, -8.07855E+00, 2.64843E-02, -2.53639E+00, 8.50922E-01, -5.20084E-01, -4.64225E-02, -1.75725E+00, 1.38469E+00, -9.22473E-01, -5.03562E-02, 6.79757E-01, -6.52431E-01, 2.33779E-01, -2.91790E-01, -1.38991E-01, 2.62270E-01, -3.51688E-03, -3.51953E-01, 2.86215E-01, -5.07076E-03, -1.96680E+00, 6.21190E-01, -1.95505E-01, -1.12009E+00, 2.77353E-02, 8.22098E-01, -2.77727E-01, -7.58262E-02, -8.15653E-02, 2.00760E-02, -1.39761E-02, 6.89437E-02, -4.42552E-03, 7.55927E-02, -8.30480E-01, 3.36159E-01, 8.98881E-01, -1.17591E+00, 3.58172E-01, -2.21041E-02, -2.33323E-02, -5.07525E-02, -5.42007E-02, 1.16181E-02, 1.09552E-02, -3.76062E-02, -1.26426E-02, 5.53849E-02, -7.10970E-02, 3.10441E-02, 1.32798E-02, 1.54776E-02, -3.14579E-02, 3.52952E-02, 1.59566E-02, -1.85110E-02, -1.01254E-02, 3.02373E-03, 4.55978E-03, 1.72477E-01, -2.61116E-01, -7.45473E-02, 8.18591E-02, -7.94097E-02, -1.04047E-05, 1.71939E-02, 1.61382E-02, 9.15953E-03, 1.70451E-02, 1.05992E-03, 1.16124E-03, -4.82049E-03, -3.61575E-03, -6.36579E-03, -6.07010E-03, -8.75332E-04])
|
|
LIBRARY.T = np.array([0.44, 1.2, 2.97, 0.67, 0.91, 5.96, 0.241, 0.69, 2.58, 0.692, 0.943, 5.8, 1.93, 1.7, 3.3, 7, 2.15, 2, 3, 7, 2.1, 4.3, 3.3, 4.7, 2.95, 0.7, 6, 1.15, 0.77, 5.84, 1.78, 2.05, 4.3, 2.43, 5.3, 2.2, 4.3, 12, 12, 13, 16, 13, 16.2, 13, 3, 2.7, 0.76, 1.48, 2.7, 6, 6, 17, 17, 0.3, 0.24, 1.8, 1.2, 0.25, 7, 8.7, 11.6, 0.45, 8.4, 8.5, 11.5, 25, 26, 0.2, 0.248, 0.2, 0.74, 3, 0.24, 2.86, 8, 17, 16, 16, 16.2, 0.7, 0.69, 7.4, 8.7, 1.25, 1.23, 4.7])
|
|
LIBRARY.D = np.array([1.0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 9])
|
|
LIBRARY.L = np.array([0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 2, 2, 3, 3, 0, 0, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 0, 0, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 0, 0, 3, 3, 1, 1, 2])
|
|
|
|
global indices
|
|
indices = set()
|
|
while len(indices) < 23:
|
|
indices.add(random.randint(0, len(LIBRARY.T)-1))
|
|
print indices, len(LIBRARY.T)
|
|
|
|
T0 = np.array([LIBRARY.T[i] for i in indices])
|
|
D0 = np.array([LIBRARY.D[i] for i in indices])
|
|
L0 = np.array([LIBRARY.L[i] for i in indices])
|
|
N0 = np.array([LIBRARY.N[i] for i in indices])
|
|
|
|
# Values from Span short(2003) (polar)
|
|
# D0 = np.array([0, 1.0, 1.0, 1.0, 3.0, 7.0, 1.0, 2.0, 5.0, 1.0, 1.0, 4.0, 2.0])
|
|
# L0 = np.array([0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 3.0])
|
|
# T0 = np.array([0, 0.25, 1.25, 1.5, 0.25, 0.875, 2.375, 2.0, 2.125, 3.5, 6.5, 4.75, 12.5])
|
|
# N0 = 0.5*np.ones_like(D0)
|
|
|
|
# values from R410A
|
|
N0 = np.array([0.0, 0.987252, -1.03017, 1.17666, -0.138991, 0.00302373, -2.53639, -1.96680, -0.830480, 0.172477, -0.261116, -0.0745473, 0.679757, -0.652431, 0.0553849, -0.0710970, -0.000875332, 0.0200760, -0.0139761, -0.0185110, 0.0171939, -0.00482049])
|
|
T0 = np.array([0.0,0.44,1.2,2.97,2.95,0.2,1.93,1.78,3.0,0.2,0.74,3.0,2.1,4.3,0.25,7.0,4.7,13.0,16.0,25.0,17.0,7.4])
|
|
D0 = np.array([0,1.0,1,1,2,5,1,2,3,5,5,5,1,1,4,4,9,2,2,4,5,6])
|
|
L0 = np.array([0,0.0,0,0,0,0,1,1,1,1,1,1,2,2,2,2,2,3,3,3,3,3])
|
|
|
|
indices = set()
|
|
while len(indices) < 5:
|
|
indices.add(random.randint(0, len(LIBRARY.T)-1))
|
|
print indices, len(LIBRARY.T)
|
|
|
|
T0 = np.append(T0, [LIBRARY.T[i] for i in indices])
|
|
D0 = np.append(D0, [LIBRARY.D[i] for i in indices])
|
|
L0 = np.append(L0, [LIBRARY.L[i] for i in indices])
|
|
N0 = np.append(N0, [LIBRARY.N[i] for i in indices])
|
|
|
|
# values from R407C
|
|
# N0 = np.array([0.0, 1.0588,-1.12018, 0.629064,-0.351953, 0.00455978,-1.75725,-1.12009, 0.0277353, 0.898881,-1.17591, 0.0818591,-0.0794097,-0.0000104047, 0.233779,-0.291790, 0.0154776,-0.0314579,-0.00442552,-0.0101254, 0.00915953,-0.003615])
|
|
# T0 = np.array([0.0,0.241,0.69,2.58,1.15,0.248,2.15,2.43,5.3,0.76,1.48,0.24,2.86,8.0,3.3,4.7,0.45,8.4,16.2,26.0,16.0,8.7])
|
|
# D0 = np.array([0.0,1,1,1,2,5,1,2,2,3,3,5,5,5,1,1,4,4,2,4,5,6])
|
|
# L0 = np.array([0.0,0,0,0,0,0,1,1,1,1,1,1,1,1,2,2,2,2,3,3,3,3])
|
|
|
|
return RefString, N0, T0, D0, L0
|
|
|
|
|
|
class IdealPartFitter(object):
|
|
def __init__(self, Ref):
|
|
self.Ref = Ref
|
|
self.RefString, N0, T0, D0, L0 = get_fluid_constants(Ref)
|
|
self.molemass = Props(self.RefString,'molemass')
|
|
self.Tc = Props(self.RefString, 'Tcrit')
|
|
self.rhoc = Props(self.RefString, 'rhocrit')
|
|
self.pc = Props(self.RefString, 'pcrit')
|
|
self.T = np.linspace(100, 450, 200)
|
|
self.tau = self.Tc/self.T
|
|
self.C = Props('C', 'T', self.T, 'D', 1e-15, self.RefString)
|
|
R = 8.314472/self.molemass
|
|
self.cp0_R = self.C/R
|
|
|
|
def cp0_R_from_fit(self, a_e):
|
|
a = a_e[0:len(a_e)//2]
|
|
e = a_e[len(a_e)//2::]
|
|
u1 = e[1]/self.T
|
|
u2 = e[2]/self.T
|
|
u3 = e[3]/self.T
|
|
return a[0]*self.T**e[0]+a[1]*u1**2*np.exp(u1)/(np.exp(u1)-1)**2+a[2]*u2**2*np.exp(u2)/(np.exp(u2)-1)**2+a[3]*u3**2*np.exp(u3)/(np.exp(u3)-1)**2
|
|
|
|
def OBJECTIVE_cp0_R(self, a_e):
|
|
cp0_R_fit = self.cp0_R_from_fit(a_e)
|
|
RMS = np.sqrt(np.mean(np.power((self.cp0_R-cp0_R_fit)/self.cp0_R, 2)))
|
|
return RMS
|
|
|
|
def fit(self):
|
|
a_e = [2.8749, 2.0623, 5.9751, 1.5612, 0.1, 697.0, 1723.0, 3875.0]
|
|
|
|
a_e = scipy.optimize.minimize(self.OBJECTIVE_cp0_R, a_e).x
|
|
|
|
self.a = a_e[0:len(a_e)//2]
|
|
self.e = a_e[len(a_e)//2::]
|
|
|
|
cp0_over_R_check = 1-self.tau**2*self.d2phi0_dTau2(self.tau)
|
|
|
|
plt.plot(self.T, (self.cp0_R_from_fit(a_e)/self.cp0_R-1)*100, '-', self.T, (cp0_over_R_check/self.cp0_R-1)*100, '^')
|
|
plt.xlabel('Temperature [K]')
|
|
plt.ylabel('($c_{p0}/R$ (fit) / $c_{p0}/R$ (REFPROP) -1)*100 [%]')
|
|
plt.savefig('cp0.pdf')
|
|
plt.close()
|
|
|
|
def d2phi0_dTau2(self, tau):
|
|
d = []
|
|
for _tau in tau:
|
|
#lead term is killed
|
|
d.append(helmholtz.phi0_logtau(-1.0).dTau2(_tau, _tau)
|
|
+ helmholtz.phi0_cp0_poly(self.a[0],self.e[0],self.Tc,298.15).dTau2(_tau, _tau)
|
|
+ helmholtz.phi0_Planck_Einstein(self.a,self.e/self.Tc,1,len(self.a)-1).dTau2(_tau, _tau)
|
|
)
|
|
return np.array(d)
|
|
|
|
|
|
class ResidualPartFitter(object):
|
|
|
|
def __init__(self, Ref, IPF):
|
|
self.Ref = Ref
|
|
self.IPF = IPF
|
|
self.RefString, self.N0, self.T0, self.D0, self.L0 = get_fluid_constants(Ref)
|
|
self.Tc = Props(self.RefString,'Tcrit')
|
|
self.rhoc = Props(self.RefString,'rhocrit')
|
|
molemass = Props(self.RefString,'molemass')
|
|
self.R = 8.314472/ molemass
|
|
|
|
def termwise_Rsquared(self):
|
|
|
|
keepers = []
|
|
values = []
|
|
print len(self.N0), 'terms at start'
|
|
for i in range(len(self.N0)):
|
|
|
|
n = helmholtz.vectord([float(1)])
|
|
d = helmholtz.vectord([self.D0[i]])
|
|
t = helmholtz.vectord([self.T0[i]])
|
|
l = helmholtz.vectord([self.L0[i]])
|
|
|
|
self.phir = helmholtz.phir_power(n, d, t, l, 0, 0)
|
|
|
|
PPF = self.evaluate_EOS(np.array(list(n)))
|
|
|
|
R2 = rsquared(PPF.p,self.phir.dDeltaV(self.tauV,self.deltaV))
|
|
|
|
values.append((R2,i))
|
|
if R2 > 0.9:
|
|
keepers.append(i)
|
|
|
|
values,indices = zip(*reversed(sorted(values)))
|
|
|
|
keepers = list(indices[0:30])
|
|
|
|
self.N0 = self.N0[keepers]
|
|
self.T0 = self.T0[keepers]
|
|
self.D0 = self.D0[keepers]
|
|
self.L0 = self.L0[keepers]
|
|
|
|
print len(self.N0), 'terms at end'
|
|
|
|
def generate_1phase_data(self):
|
|
|
|
Tc = Props(self.RefString, 'Tcrit')
|
|
rhoc = Props(self.RefString, 'rhocrit')
|
|
TTT, RHO, PPP, CPP, CVV, AAA = [], [], [], [], [], []
|
|
|
|
for _T in np.linspace(220, 450, 100):
|
|
print _T
|
|
for _rho in np.logspace(np.log10(1e-2), np.log10(rhoc), 100):
|
|
try:
|
|
if _T > Tc:
|
|
p = Props('P', 'T', _T, 'D', _rho, self.RefString)
|
|
cp = Props('C', 'T', _T, 'D', _rho, self.RefString)
|
|
cv = Props('O', 'T', _T, 'D', _rho, self.RefString)
|
|
a = Props('A', 'T', _T, 'D', _rho, self.RefString)
|
|
else:
|
|
DL = Props('D', 'T', _T, 'Q', 0, self.RefString)
|
|
DV = Props('D', 'T', _T, 'Q', 1, self.RefString)
|
|
if _rho < DV or _rho > DL:
|
|
p = Props('P', 'T', _T, 'D', _rho, self.RefString)
|
|
cp = Props('C', 'T', _T, 'D', _rho, self.RefString)
|
|
cv = Props('O', 'T', _T, 'D', _rho, self.RefString)
|
|
a = Props('A', 'T', _T, 'D', _rho, self.RefString)
|
|
else:
|
|
p = None
|
|
if p is not None:
|
|
TTT.append(_T)
|
|
RHO.append(_rho)
|
|
PPP.append(p)
|
|
CPP.append(cp)
|
|
CVV.append(cv)
|
|
AAA.append(a)
|
|
|
|
except ValueError as VE:
|
|
print VE
|
|
pass
|
|
|
|
for _rho in np.linspace(rhoc, 3.36*rhoc, 50):
|
|
try:
|
|
if _T > Tc:
|
|
p = Props('P', 'T', _T, 'D', _rho, self.RefString)
|
|
cp = Props('C', 'T', _T, 'D', _rho, self.RefString)
|
|
cv = Props('O', 'T', _T, 'D', _rho, self.RefString)
|
|
a = Props('A', 'T', _T, 'D', _rho, self.RefString)
|
|
else:
|
|
DL = Props('D', 'T', _T, 'Q', 0, self.RefString)
|
|
DV = Props('D', 'T', _T, 'Q', 1, self.RefString)
|
|
if _rho < DV or _rho > DL:
|
|
p = Props('P', 'T', _T, 'D', _rho, self.RefString)
|
|
cp = Props('C', 'T', _T, 'D', _rho, self.RefString)
|
|
cv = Props('O', 'T', _T, 'D', _rho, self.RefString)
|
|
a = Props('A', 'T', _T, 'D', _rho, self.RefString)
|
|
else:
|
|
p = None
|
|
if p is not None:
|
|
TTT.append(_T)
|
|
RHO.append(_rho)
|
|
PPP.append(p)
|
|
CPP.append(cp)
|
|
CVV.append(cv)
|
|
AAA.append(a)
|
|
|
|
except ValueError as VE:
|
|
print VE
|
|
pass
|
|
|
|
h = h5py.File('T_rho_p.h5','w')
|
|
grp = h.create_group(self.Ref)
|
|
grp.create_dataset("T",data = np.array(TTT),compression = "gzip")
|
|
grp.create_dataset("rho", data = np.array(RHO),compression = "gzip")
|
|
grp.create_dataset("p", data = np.array(PPP),compression = "gzip")
|
|
grp.create_dataset("cp", data = np.array(CPP),compression = "gzip")
|
|
grp.create_dataset("cv", data = np.array(CVV),compression = "gzip")
|
|
grp.create_dataset("speed_sound", data = np.array(AAA),compression = "gzip")
|
|
h.close()
|
|
|
|
def load_data(self):
|
|
h = h5py.File('T_rho_p.h5','r')
|
|
self.T = h.get(self.Ref + '/T').value
|
|
self.rho = h.get(self.Ref + '/rho').value
|
|
self.p = h.get(self.Ref + '/p').value
|
|
self.cp = h.get(self.Ref + '/cp').value
|
|
self.cv = h.get(self.Ref + '/cv').value
|
|
self.speed_sound = h.get(self.Ref + '/speed_sound').value
|
|
|
|
self.tau = self.Tc/self.T
|
|
self.delta = self.rho/self.rhoc
|
|
self.tauV = helmholtz.vectord(self.tau)
|
|
self.deltaV = helmholtz.vectord(self.delta)
|
|
|
|
# Get the derivative d2phi0_dTau2 from the ideal part fitter
|
|
self.d2phi0_dTau2 = self.IPF.d2phi0_dTau2(self.tau)
|
|
|
|
def evaluate_EOS(self, N):
|
|
|
|
self.phir.n = helmholtz.vectord(N)
|
|
|
|
dDelta = self.phir.dDeltaV(self.tauV, self.deltaV)
|
|
dTau2 = self.phir.dTau2V(self.tauV, self.deltaV)
|
|
dDelta2 = self.phir.dDelta2V(self.tauV, self.deltaV)
|
|
dDelta_dTau = self.phir.dDelta_dTauV(self.tauV, self.deltaV)
|
|
|
|
# Evaluate the pressure
|
|
p = (self.rho*self.R*self.T)*(1 + self.delta*dDelta)
|
|
# Evaluate the specific heat at constant volume
|
|
cv_over_R = -self.tau**2*(self.d2phi0_dTau2 + dTau2)
|
|
cv = cv_over_R*self.R
|
|
# Evaluate the specific heat at constant pressure
|
|
cp_over_R = cv_over_R+(1.0+self.delta*dDelta-self.delta*self.tau*dDelta_dTau)**2/(1+2*self.delta*dDelta+self.delta**2*dDelta2)
|
|
cp = cp_over_R*self.R
|
|
# Evaluate the speed of sound
|
|
w = np.sqrt(1000*self.R*self.T*cp_over_R/cv_over_R*(1+2*self.delta*dDelta+self.delta**2*dDelta2))
|
|
|
|
class stub: pass
|
|
PPF = stub()
|
|
PPF.p = np.array(p, ndmin = 1).T
|
|
PPF.cp = np.array(cp, ndmin = 1).T
|
|
PPF.cv = np.array(cv, ndmin = 1).T
|
|
PPF.w = np.array(w, ndmin = 1).T
|
|
|
|
return PPF
|
|
|
|
def OBJECTIVE(self, N):
|
|
PPF = self.evaluate_EOS(N)
|
|
|
|
## plt.plot(PPF.p, self.p); plt.show()
|
|
## plt.plot(PPF.cp, self.cp); plt.show()
|
|
## plt.plot(PPF.cv, self.cv); plt.show()
|
|
## plt.plot(PPF.w, self.speed_sound); plt.show()
|
|
w_p = 1.0
|
|
w_cv = 1.0
|
|
w_w = 1.0
|
|
w_cp = 1.0
|
|
w_total = (w_p+w_cv+w_w+w_cp)/4
|
|
|
|
w_p_norm = w_p/w_total
|
|
w_cv_norm = w_cv/w_total
|
|
w_cp_norm = w_cp/w_total
|
|
w_w_norm = w_w/w_total
|
|
residuals = np.r_[(PPF.p/self.p-1),(PPF.cv/self.cv-1),(PPF.cp/self.cp-1)]#,(PPF.w**2/self.speed_sound**2-1)]
|
|
RMS = np.sqrt(np.mean(np.power(residuals, 2)))
|
|
|
|
print 'RMS:',RMS*100, '% Max',np.max(np.abs(residuals))*100,'%'
|
|
self.RMS = RMS
|
|
self.MaxError = np.max(np.abs(residuals))
|
|
return RMS
|
|
|
|
def fit(self):
|
|
|
|
# Kill off some not as good terms
|
|
#self.termwise_Rsquared()
|
|
|
|
# Load up the residual Helmholtz term with parameters
|
|
n = helmholtz.vectord(self.N0)
|
|
d = helmholtz.vectord(self.D0)
|
|
t = helmholtz.vectord(self.T0)
|
|
l = helmholtz.vectord(self.L0)
|
|
self.phir = helmholtz.phir_power(n, d, t, l, 1, len(self.N0)-1)
|
|
|
|
# Solve for the coefficients
|
|
Nbounds = [(-10,10) for _ in range(len(self.N0))]
|
|
tbounds = [(-1,30) for _ in range(len(self.T0))]
|
|
print self.OBJECTIVE(np.array(list(self.N0)))
|
|
#self.N = self.N0
|
|
#self.N = scipy.optimize.minimize(self.OBJECTIVE, np.array(list(self.N0)), bounds = Nbounds, options = dict(maxiter = 5)).x
|
|
self.N = scipy.optimize.minimize(self.OBJECTIVE, np.array(list(self.N0)), method = 'L-BFGS-B', bounds = Nbounds, options = dict(maxiter = 100)).x
|
|
|
|
# Write the coefficients to HDF5 file
|
|
h = h5py.File('fit_coeffs.h5','w')
|
|
grp = h.create_group(self.Ref)
|
|
grp.create_dataset("n", data = np.array(self.N), compression = "gzip")
|
|
print self.N
|
|
#grp.create_dataset("t", data = np.array(self.N[len(self.N)//2::]), compression = "gzip")
|
|
h.close()
|
|
|
|
def evaluate_REFPROP(self, Ref, T, rho):
|
|
|
|
p,cp,cv,w = [],[],[],[]
|
|
R = 8.314472/Props(Ref,'molemass')
|
|
for _T,_rho in zip(T, rho):
|
|
p.append(Props("P",'T',_T,'D',_rho,Ref))
|
|
cp.append(Props("C",'T',_T,'D',_rho,Ref))
|
|
cv.append(Props("O",'T',_T,'D',_rho,Ref))
|
|
w.append(Props("A",'T',_T,'D',_rho,Ref))
|
|
|
|
class stub: pass
|
|
PPF = stub()
|
|
PPF.p = np.array(p, ndmin = 1).T
|
|
PPF.cp = np.array(cp, ndmin = 1).T
|
|
PPF.cv = np.array(cv, ndmin = 1).T
|
|
PPF.w = np.array(w, ndmin = 1).T
|
|
|
|
return PPF
|
|
|
|
def check(self):
|
|
# Load the coefficients from file
|
|
h = h5py.File('fit_coeffs.h5','r')
|
|
grp = h.get(self.Ref)
|
|
n = grp.get('n').value
|
|
h.close()
|
|
|
|
print n
|
|
|
|
import matplotlib.colors as colors
|
|
cNorm = colors.LogNorm(vmin=1e-3, vmax=50)
|
|
PPF = self.evaluate_EOS(np.array(list(n)))
|
|
self.OBJECTIVE(np.array(list(n)))
|
|
|
|
print 'max error (p)',np.max(np.abs(PPF.p/self.p-1)*100),'%'
|
|
SC1 = plt.scatter(self.rho, self.T, s = 8, c = np.abs(PPF.p/self.p-1)*100, edgecolors = 'none', cmap = plt.get_cmap('jet'), norm = cNorm)
|
|
plt.gca().set_xscale('log')
|
|
cb = plt.colorbar()
|
|
cb.set_label('abs(PPF.p/self.p-1)*100')
|
|
plt.savefig('pressure.png')
|
|
plt.show()
|
|
|
|
print 'max error (cp)',np.max(np.abs(PPF.cp/self.cp-1)*100),'%'
|
|
SC1 = plt.scatter(self.rho, self.T, s = 8, c = np.abs(PPF.cp/self.cp-1)*100, edgecolors = 'none', cmap = plt.get_cmap('jet'), norm = cNorm)
|
|
plt.gca().set_xscale('log')
|
|
cb = plt.colorbar()
|
|
cb.set_label('abs(PPF.cp/self.cp-1)*100')
|
|
plt.savefig('cp.png')
|
|
plt.show()
|
|
|
|
## plt.plot(self.T,PPF.p/self.p,'.'); plt.show()
|
|
## plt.plot(self.T,PPF.cp/self.cp,'.'); plt.show()
|
|
## plt.plot(self.T,PPF.cv/self.cv,'.'); plt.show()
|
|
## plt.plot(self.T,PPF.w/self.speed_sound,'.'); plt.show()
|
|
|
|
|
|
class PPFFitterClass(object):
|
|
|
|
def __init__(self, Ref, regenerate_data = True, fit = True):
|
|
|
|
self.Ref = Ref
|
|
|
|
self.IPF = IdealPartFitter(Ref)
|
|
self.IPF.fit()
|
|
for i in range(1):
|
|
|
|
self.RPF = ResidualPartFitter(Ref, IPF = self.IPF)
|
|
if regenerate_data:
|
|
self.RPF.generate_1phase_data()
|
|
|
|
self.RPF.load_data()
|
|
|
|
if fit:
|
|
self.RPF.fit()
|
|
|
|
f = open('results.txt','a+')
|
|
print >> f, indices, self.RPF.RMS, self.RPF.MaxError
|
|
f.close()
|
|
|
|
self.RPF.check()
|
|
quit()
|
|
self.output_files()
|
|
|
|
def contour_plot(values):
|
|
"""
|
|
Parameters
|
|
----------
|
|
values : iterable, same size as T and rho
|
|
"""
|
|
|
|
plt.semilogx(self.RPF.rho,self.RPF.T,'o')
|
|
plt.show()
|
|
|
|
# Generate a regular grid to interpolate the data.
|
|
xi = np.linspace(min(self.RPF.T), max(self.RPF.T), 100)
|
|
yi = np.linspace(min(self.RPF.rho), max(self.RPF.rho), 100)
|
|
xi, yi = np.meshgrid(xi, yi)
|
|
# Interpolate using delaunay triangularization
|
|
zi = mlab.griddata(np.array(self.RPF.T),np.array(self.RPF.rho),np.array(values),xi,yi)
|
|
cont = plt.contourf(yi,xi,zi,30)
|
|
plt.colorbar()
|
|
plt.show()
|
|
|
|
def output_files(self):
|
|
h = h5py.File('fit_coeffs.h5','r')
|
|
n = h.get(self.Ref+'/n').value
|
|
#t = h.get(self.Ref+'/t').value
|
|
|
|
# Output the header file
|
|
header = PPF_h_template.format(Ref = self.Ref, RefUpper = self.Ref.upper())
|
|
|
|
acoeffs = '0, '+', '.join(['{a:0.6f}'.format(a=_) for _ in self.IPF.a])
|
|
# First one doesn't get divided by critical temperature, later ones do
|
|
bcoeffs = '0, '
|
|
bcoeffs += str(self.IPF.e[0])+', '
|
|
bcoeffs += ', '.join(['{b:0.4f}/{Tcrit:g}'.format(b=_,Tcrit = self.IPF.Tc) for _ in self.IPF.e[1::]])
|
|
|
|
ncoeffs = ', '.join(['{a:0.6g}'.format(a=_) for _ in n])
|
|
tcoeffs = ', '.join(['{a:0.6g}'.format(a=_) for _ in self.RPF.T0])
|
|
dcoeffs = ', '.join(['{a:0.6g}'.format(a=_) for _ in self.RPF.D0])
|
|
lcoeffs = ', '.join(['{a:0.6g}'.format(a=_) for _ in self.RPF.L0])
|
|
|
|
import sys
|
|
sys.path.append('..')
|
|
from fit_ancillary_ODRPACK import saturation_pressure, saturation_density
|
|
pL = saturation_pressure(self.IPF.RefString, self.IPF.Ref, LV = 'L')
|
|
pV = saturation_pressure(self.IPF.RefString, self.IPF.Ref, LV = 'V')
|
|
rhoL = saturation_density(self.IPF.RefString, self.IPF.Ref, form='A', LV='L', add_critical = False)
|
|
rhoV = saturation_density(self.IPF.RefString, self.IPF.Ref, form='B', LV='V', add_critical = False)
|
|
|
|
code = PPF_cpp_template.format(Ref = self.Ref,
|
|
RefUpper = self.Ref.upper(),
|
|
acoeffs = acoeffs,
|
|
bcoeffs = bcoeffs,
|
|
Ncoeffs = ncoeffs,
|
|
tcoeffs = tcoeffs,
|
|
dcoeffs = dcoeffs,
|
|
Lcoeffs = lcoeffs,
|
|
N_phir = len(n),
|
|
N_cp0 = len(self.IPF.a),
|
|
molemass = self.IPF.molemass,
|
|
Ttriple = 200,
|
|
accentric = 0.7,
|
|
pcrit = self.IPF.pc,
|
|
Tcrit = self.IPF.Tc,
|
|
rhocrit = self.IPF.rhoc,
|
|
pL = pL,
|
|
pV = pV,
|
|
rhoL = rhoL,
|
|
rhoV = rhoV
|
|
)
|
|
|
|
f = open(self.IPF.Ref+'.h','w')
|
|
f.write(header)
|
|
f.close()
|
|
|
|
f = open(self.IPF.Ref+'.cpp','w')
|
|
f.write(code)
|
|
f.close()
|
|
|
|
|
|
if __name__=='__main__':
|
|
Ref = 'R407F'
|
|
|
|
PPFFitterClass(Ref)
|