In [1]:
#%matplotlib inline
%matplotlib notebook
%load_ext autoreload
%autoreload 2
from __future__ import unicode_literals
import numpy as np
import matplotlib.pyplot as plt
plt.rc('font',family='serif')
#import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import CoolProp
from ammonia_props import AmmoniaProps, massFractionToMolar
import tabulate
In [2]:
throwaway = CoolProp.AbstractState("REFPROP","water")
In [3]:
# Patek and Klomfar, my implementation
# TODO
# Ibrahim and Klein, 1995
myammonia = AmmoniaProps()
def wraph(**kwargs):
try:
return 1000*myammonia.h(**kwargs)
except:
return np.nan
hfun1 = np.vectorize(wraph)
# Tillner-Roth, 1998
#amm = lambda(x):'REFPROP::water[{}]&ammonia[{}]'.format(1-x,x)
CPRP = CoolProp.AbstractState("REFPROP","water&ammonia")
def hsat(T,Q,w,update=True):
#h = CP.PropsSI('H','T',T,'Q',Q,amm(x))
if update:
x = massFractionToMolar(w)
CPRP.set_mole_fractions([1-x,x])
try:
CPRP.update(CoolProp.QT_INPUTS,Q,T)
return CPRP.hmass()
except:
return np.nan
hsatv = np.vectorize(hsat)
In [4]:
T_ref = 273.15
h_amm_ref_rp = CoolProp.CoolProp.PropsSI("H","T",T_ref,"Q",0,"REFPROP::ammonia")
h_amm_ref_ees = wraph(T=T_ref,x=1,Qu=0)
h_offset_amm = h_amm_ref_rp - h_amm_ref_ees
print("{}\n- {}\n=======\n{}".format(h_amm_ref_rp, h_amm_ref_ees, h_offset_amm))
In [ ]:
In [5]:
help(AmmoniaProps)
In [9]:
x_w = np.logspace(-4,0)
x = 1-x_w
h_offset = x * h_offset_amm
fig1 = plt.figure(1)
plt.gca().ticklabel_format(axis='y', style = 'sci', scilimits=(-2,2), useOffset=False)
plt.xlabel('Ammonia composition $x$ [kg/kg]')
plt.ylabel('Enthalpy $h$ [J/kg]')
plt.annotate('Saturated vapor',[0.02,2e6])
plt.annotate('Saturated liquid',[0.02,1e6])
fig2 = plt.figure(2)
plt.gca().ticklabel_format(axis='y', style = 'sci', scilimits=(-2,2), useOffset=False)
plt.xlabel('Water composition $x$ [kg/kg]')
plt.ylabel('Enthalpy $h$ [J/kg]')
first = True
cc = 'b g r c m y'.split()
TT = np.linspace(273.15,573.15,6)
ees_hl_by_x = {}
ees_hv_by_x = {}
rp_hl_by_x = {}
rp_hv_by_x = {}
for T,c in zip(TT,cc):
#for T in [373]:
print(T)
ees_hl_by_x[T] = hfun1(T=T,Qu=0,x=x) + h_offset
ees_hv_by_x[T] = hfun1(T=T,Qu=1,x=x) + h_offset
rp_hl_by_x[T] = hsatv(T,0,x)
rp_hv_by_x[T] = hsatv(T,1,x)
plt.show()
In [10]:
for T,c in zip(TT,cc):
#for T in [373]:
plt.figure(1)
plt.plot(x,ees_hl_by_x[T],'--',color=c,label="EES+offset" if first else None)
plt.plot(x,ees_hv_by_x[T],'--',color=c)
plt.plot(x,rp_hl_by_x[T],'-' ,color=c,label="REFPROP" if first else None)
plt.plot(x,rp_hv_by_x[T],'-' ,color=c)
xtext = 0.4 + (T-273.15) * 0.001
htext = hfun1(T=T+5,Qu=0,x=xtext)
plt.annotate('{} $^\circ$C'.format(T-273.15),[xtext,htext])
plt.figure(2)
plt.semilogx(x_w,ees_hl_by_x[T],'--',color=c)
plt.semilogx(x_w,ees_hv_by_x[T],'--',color=c)
plt.semilogx(x_w,rp_hl_by_x[T],'-' ,color=c)
plt.semilogx(x_w,rp_hv_by_x[T],'-' ,color=c)
plt.annotate('{} $^\circ$C'.format(T-273.15),[1-xtext,htext])
first = False
plt.figure(1)
plt.legend(loc='best')
plt.show()
In [12]:
# Saturation curve of pure ammonia, because REFPROP is not good near there
rp_amm = CoolProp.AbstractState("HEOS","ammonia")
rp_amm.build_phase_envelope("")
def h_amm_rp(T,Q):
try:
rp_amm.update(CoolProp.QT_INPUTS,Q,T)
return rp_amm.hmass()
except:
return np.nan
hamham = np.vectorize(h_amm_rp)
# Pure water, same
rp_water = CoolProp.AbstractState("HEOS","water")
rp_water.build_phase_envelope("")
def h_water_rp(T,Q):
try:
rp_water.update(CoolProp.QT_INPUTS,Q,T)
return rp_water.hmass()
except:
return np.nan
humhum = np.vectorize(h_water_rp)
T_c = {}
P_c = {}
h_c = {}
T_c[1] = rp_amm.T_critical()
P_c[1] = rp_amm.p_critical()
rp_amm.specify_phase(CoolProp.iphase_critical_point)
rp_amm.update(CoolProp.PT_INPUTS,P_c[1],T_c[1])
h_c[1] = rp_amm.hmass()
rp_amm.specify_phase(CoolProp.iphase_unknown)
T_c[0] = rp_water.T_critical()
P_c[0] = rp_water.p_critical()
rp_water.specify_phase(CoolProp.iphase_critical_point)
rp_water.update(CoolProp.PT_INPUTS,P_c[0],T_c[0])
h_c[0] = rp_water.hmass()
rp_water.specify_phase(CoolProp.iphase_unknown)
ees_hl_by_T = {}
ees_hv_by_T = {}
rp_hl_by_T = {}
rp_hv_by_T = {}
x = 1
Tee = np.linspace(273,T_c[1],100)
xx=x*np.ones_like(Tee)
rp_hl_by_T[x] = hamham(Tee,0)
rp_hv_by_T[x] = hamham(Tee,1)
for x in [0.2,0.4,0.6,0.8,0.9,0.99,0.999]:
xx = x*np.ones_like(Tee)
hsatv(300,0,x)
T_c[x] = CPRP.T_critical()
P_c[x] = CPRP.p_critical()
#CPRP.specify_phase(CoolProp.constants.iphase_critical_point)
CPRP.update(CoolProp.PT_INPUTS,P_c[x],T_c[x])
h_c[x] = CPRP.hmass()
#CPRP.specify_phase(CoolProp.constants.iphase_twophase)
print("Tc_computed = {}".format(T_c[x]))
Tee = np.linspace(273,T_c[x],100)
rp_hl_by_T[x] = hsatv(Tee,0,x,False)
rp_hv_by_T[x] = hsatv(Tee,1,x,False)
x = 0
Tee = np.linspace(273, T_c[0],100)
xx=x*np.ones_like(Tee)
rp_hl_by_T[x] = humhum(Tee,0)
rp_hv_by_T[x] = humhum(Tee,1)
In [15]:
keys = list(P_c.keys())
keys.sort()
critical_data = np.array([(k,P_c[k],T_c[k],h_c[k]) for k in keys])
print(tabulate.tabulate(critical_data,"x P T h".split()))
critical_data[:,0]
Out[15]:
In [16]:
import ammonia1
In [17]:
chiller = ammonia1.AmmoniaChiller()
chiller.update()
display(chiller)
cst=chiller.stateTable()
cst['T']
Out[17]:
In [19]:
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.ticklabel_format(axis='x', style = 'sci', scilimits=(-2,2), useOffset=False)
plt.xlabel('Mass specific enthalpy [J/kg]')
plt.ylabel('Ammonia mass fraction')
ax.set_zlabel('Temperature [$^\circ$C]')
x = 1
Tee = np.linspace(273,T_c[x],100)
xx=x*np.ones_like(Tee)
plt.plot(rp_hl_by_T[x],xx,Tee-273.15,'b')
plt.plot(rp_hv_by_T[x],xx,Tee-273.15,'r')
#for x in [0.2,0.4,0.6,0.8,0.9,0.99,0.999]:
for x in [0.4,0.6,0.8,0.9,0.99,0.999]:
Tee = np.linspace(273,T_c[x],100)
xx = x*np.ones_like(Tee)
indices = Tee < 273+200
plt.plot(rp_hl_by_T[x][indices],xx[indices],Tee[indices]-273.15,'b')
plt.plot(rp_hv_by_T[x][indices],xx[indices],Tee[indices]-273.15,'r')
x = 0
Tee = np.linspace(273, T_c[x],100)
xx=x*np.ones_like(Tee)
#plt.plot(rp_hl_by_T[x],xx,Tee-273.15,'b')
#plt.plot(rp_hv_by_T[x],xx,Tee-273.15,'r')
x_w = np.logspace(-4,0)
x = 1-x_w
#for T,c in zip(TT,cc)[:3]:
for T,c in zip(TT,cc):
#for T in [373]:
indices = x >= 0.3
plt.plot(ees_hl_by_x[T][indices],x[indices],T-273.15,'b--')
plt.plot(ees_hv_by_x[T][indices],x[indices],T-273.15,'r--')
plt.plot(rp_hl_by_x[T][indices],x[indices],T-273.15,'b-' )
plt.plot(rp_hv_by_x[T][indices],x[indices],T-273.15,'r-' )
indices = critical_data[:,2] < 500
plt.plot(critical_data[indices,3],critical_data[indices,0],critical_data[indices,2]-273.15,
'k')
h = cst['h']*1e3 + cst['x'] * h_offset_amm
plt.plot(h,cst['x'],cst['T']-273.15,'o')
for p,x,y,z in zip(chiller.points,h,cst['x'],cst['T']-273.15):
#ax.text(x,y,z, p, size=20, zorder=1, color='k')
pass
from Arrow3D import Arrow3D
for (i,j) in chiller.getPaths():
#plt.plot(h[[i,j]],cst['x'][[i,j]],cst['T'][[i,j]]-273.15)
a=Arrow3D(h[[i,j]],cst['x'][[i,j]],cst['T'][[i,j]]-273.15,mutation_scale=20, lw=1, arrowstyle="-|>", color="k")
ax.add_artist(a)
ax.margins(0)
ax.autoscale_view('tight')
ax.view_init(20,-36)
fig.tight_layout()
plt.savefig('../img/view1.svg')
ax.view_init(20,-73)
fig.tight_layout()
plt.savefig('../img/view2.svg')
ax.view_init(20,-109)
fig.tight_layout()
plt.savefig('../img/view3.svg')
ax.set_ylim([0.4,1])
ax.set_zlim([0,200])
plt.show()
In [ ]:
In [20]:
from Arrow3D import Arrow3D
fig = plt.figure()
ax = fig.gca(projection='3d')
#ax.quiver3D(6,6,6,1,1,1)
a=Arrow3D((5,6),(5,6),(5,6),mutation_scale=20, lw=1, arrowstyle="-|>", color="k")
ax.add_artist(a)
ax.set_xlim([0,10])
ax.set_ylim([0,10])
ax.set_zlim([0,10])
plt.show()
In [22]:
print(CoolProp.CoolProp.PropsSI("H","T",373,"Q",0,"REFPROP::ammonia"))
print(CoolProp.CoolProp.PropsSI("H","T",195.495,"Q",0,"REFPROP::ammonia[1]&water[0]"))
print(CoolProp.CoolProp.PropsSI("H","T",195.495,"Q",0,"REFPROP::ammonia"))
print(CoolProp.CoolProp.PropsSI("H","T",273.16,"Q",0,"REFPROP::ammonia[0]&water[1]"))
print(CoolProp.CoolProp.PropsSI("H","T",273.16,"Q",0,"REFPROP::water"))
In [23]:
CoolProp.CoolProp.set_reference_state('ammonia','DEF')
CoolProp.CoolProp.PropsSI("H","T",373,"Q",0,"HEOS::ammonia")
Out[23]:
In [24]:
CPRP.set_mole_fractions([0.001,0.999])
CPRP.update(CoolProp.QT_INPUTS,0,373)
CPRP.mole_fractions_liquid()
Out[24]:
In [25]:
myammonia.h(T=273.15,x=1,Qu=0)
Out[25]:
In [28]:
print(CoolProp.CoolProp.PropsSI("H","T",273.16,"Q",0,"REFPROP::ammonia[0]&water[1]"))
print(CoolProp.CoolProp.PropsSI("H","T",273.16,"Q",0,"REFPROP::water"))
In [29]:
CoolProp.CoolProp.PropsSI("T","H",2e6,"Q",0,"REFPROP::ammonia[1]&water[0]")
In [30]:
reload(ammonia1)
In [ ]: