In [1]:
from __future__ import division
import numpy as np
In [19]:
from math import exp, fsum
b = np.zeros((33,), dtype = np.double)
fluid = 'R152a'
if fluid == 'R123':
# Start with the coefficients from Younglove and McLinden as written
b = [0,-6.57453133659e-3,2.93479845842,-9.89140469845e1,2.01029776013e4,-3.83566527886e6,2.27587641969e-3,-9.08726819450,4.34181417995e3,3.54116464954e6,-6.35394849670e-4,3.20786715274,-1.31276484299e3,-1.16360713718e-1,-1.13354409016e1,-5.37543457327e3,2.58112416120,
-1.06148632128e-1,5.00026133667e1,-2.04326706346,-2.49438345685e6,-4.63962781113e8,-2.84903429588e5,9.74392239902e9,-6.37314379308e3,3.14121189813e5,-1.45747968225e2,-8.43830261449e6,-2.41138441593,1.08508031257e3,-1.06653193965e-2,-1.21343571084e1,-2.57510383240e2]
Tc = 456.831
rhoc = 3.596417
R = 0.08314510
M = 152.930
gamma = 1
elif fluid == 'R152a':
b = [0,-0.101623317e-1,+0.215677130e1,-0.648581254e2,+0.122535596e5,-0.206805988e7,-0.379836507e-3,-0.441333233e0,+0.158248875e3,+0.564062216e6,-0.124115350e-3,+0.494972179e0,-0.208058040e3,-0.131403187e-1,+0.212083849e0,-0.151263785e3,+0.311108025e-1,-0.115280980e-2,+0.437040026e0,-0.965596535e-2,-0.242705525e6,-0.518042519e8,-0.119070546e5,+0.459333195e9,-0.719317287e2,-0.840102861e4,-0.102910957e1,-0.325913881e5,-0.412362182e-2,+0.175102808e1,-0.198636625e-4,-0.421363036e-2,-0.198696761e1]
Tc = 386.411
rhoc = 5.57145
R = 0.08314471
M = 66.051
gamma = 1
else:
raise ValueError(fluid)
# Coefficients for Equation 3.28 from Span, 2000
n_2i = b[:]
t_2i = np.array([0, 1,0.5,0,-1,-2, 1,0,-1,-2, 1,0,-1, 0, -1,-2, -1, -1,-2, -2, -2,-3, -2,-4, -2,-3, -2,-4, -2,-3, -2,-3,-4],dtype = np.longdouble)
d_2i = np.array([0, 2,2,2,2,2, 3,3,3,3, 4,4,4, 5, 6,6, 7, 8,8, 9, 3,3, 5,5, 7,7, 9,9, 11,11, 13,13,13],dtype = np.longdouble)
### Translate to the common Equation 3.26 formulation
n_3_26 = [0] + [n_2i[i]*rhoc**(d_2i[i]-1)*Tc**(t_2i[i]-1)/R for i in range(1,33)]
t_3_26 = [0] + [-t_2i[i]+1 for i in range(1,33)]
d_3_26 = [0] + [d_2i[i]-1 for i in range(1,33)]
# First polynomial part is easy (see Span, 2000, page 26)
n = [0] + [n_3_26[i]/d_3_26[i] for i in range(1,20)]
t = [0] + t_3_26[1:20]
d = [0] + d_3_26[1:20]
# Next part is harder
n += [0]*21
d += [0,0,0,0,0,0,2,2,2,4,4,4,6,6,6,8,8,8,10,10,10]
t += [3,4,5,3,4,5,3,4,5,3,4,5,3,4,5,3,4,5,3,4,5]
nn = n_3_26[:]
n[20] = nn[20]/(2.0*gamma)+nn[22]/(2.0*gamma**2)+nn[24]/(gamma**3)+3*nn[26]/(gamma**4)+12*nn[28]/(gamma**5)+60*nn[30]/(gamma**6)
n[21] = nn[21]/(2*gamma**2)+nn[25]/(gamma**3)+12*nn[29]/(gamma**5)+60*nn[31]/(gamma**6)
n[22] = nn[23]/(2*gamma**2)+3*nn[27]/(gamma**4)+60*nn[32]/(gamma**6)
n[23] = -n[20]
n[24] = -n[21]
n[25] = -n[22]
n[26] = -nn[22]/(2*gamma)-nn[24]/(gamma**2)-3*nn[26]/(gamma**3)-12*nn[28]/(gamma**4)-60*nn[30]/(gamma**5)
n[27] = -nn[25]/(gamma**2)-12*nn[29]/(gamma**4)-60*nn[31]/(gamma**5)
n[28] = -nn[23]/(2*gamma)-3*nn[27]/(gamma**3)-60*nn[32]/(gamma**5)
n[29] = -nn[24]/(2*gamma)-3*nn[26]/(2*gamma**2)-6*nn[28]/(gamma**3)-30*nn[30]/(gamma**4)
n[30] = -nn[25]/(2*gamma)-6*nn[29]/(gamma**3)-30*nn[31]/(gamma**4)
n[31] = -3*nn[27]/(2*gamma**2)-30*nn[32]/(gamma**4)
n[32] = -nn[26]/(2*gamma)-2*nn[28]/(gamma**2)-10*nn[30]/(gamma**3)
n[33] = -2*nn[29]/(gamma**2)-10*nn[31]/(gamma**3)
n[34] = -nn[27]/(2*gamma)-10*nn[32]/(gamma**3)
n[35] = -nn[28]/(2*gamma)-5*nn[30]/(2*gamma**2)
n[36] = -nn[29]/(2*gamma)-5*nn[31]/(2*gamma**2)
n[37] = -5*nn[32]/(2.0*gamma**2)
n[38] = -nn[30]/(2.0*gamma)
n[39] = -nn[31]/(2.0*gamma)
n[40] = -nn[32]/(2.0*gamma)
l = [0] + [0 for i in range(1,23)] + [2 for i in range(23,41)]
# Reorder the terms in REFPROP ordering
n = [0] + n[20:23] + n[1:20] + n[23::]
d = [0] + d[20:23] + d[1:20] + d[23::]
t = [0] + t[20:23] + t[1:20] + t[23::]
l = [0] + l[20:23] + l[1:20] + l[23::]
print 'n:', n[1:41]
print 't:', t[1:41]
print 'd:', d[1:41]
print 'l:', l[1:41]
for i in range(1,41):
print n[i],t[i],d[i],l[i]
if fluid =='R123':
def p_328(T_K, rho_molL):
T = T_K
rho = rho_molL
s = n_2i*T**t_2i*rho**d_2i
sum1 = fsum(s[1:20])
sum2 = fsum(s[20:32])
# Form of equation 3.28 from Span, 2000
return rho*R*T + sum1 + exp(-(rho/rhoc)**2)*sum2
print p_328(146+273.15, 129.92/M), 'should yield', 19.5614, 'to within rounding error'
def p_326(T_K, rho_molL):
T = T_K
rho = rho_molL
tau = Tc/T
delta = rho/rhoc
sum1 = fsum([n_3_26[i]*tau**t_3_26[i]*delta**d_3_26[i] for i in range(1,20)])
sum2 = fsum([n_3_26[i]*tau**t_3_26[i]*delta**d_3_26[i] for i in range(20,33)])
# Form of equation 3.26 from Span, 2000
return rho*R*T*(1+sum1 + exp(-(rho/rhoc)**2)*sum2)
print p_326(146+273.15, 129.92/M), 'should yield', 19.5614, 'to within rounding error'
def p_325(T_K, rho_molL):
T = T_K
rho = rho_molL
tau = Tc/T
delta = rho/rhoc
sum1 = fsum([n[i]*d[i]*tau**t[i]*delta**(d[i]-1) for i in range(1,23)])
sum2 = fsum([n[i]*delta**(d[i]-1)*(d[i]-2*delta**2)*exp(-delta**2)*tau**t[i] for i in range(23,41)])
daddelta = sum1 + sum2
print tau, delta, daddelta
# Form of equation 3.26 from Span, 2000
return rho*R*T*(1+delta*daddelta)
print p_325(146+273.15, 129.92/M), 'should yield', 19.5614, 'to within rounding error'
In [ ]: