| 1 |
from fprops import * |
| 2 |
from pylab import * |
| 3 |
import sys |
| 4 |
|
| 5 |
#P = fluid('water','helmholtz'); |
| 6 |
#P = fluid('ammonia','pengrob'); |
| 7 |
P = fluid('isohexane','pengrob'); |
| 8 |
|
| 9 |
print "SOLVING TRIPLE POINT..." |
| 10 |
|
| 11 |
print "Fluid: %s\nData source: %s" %(P.name, P.source) |
| 12 |
|
| 13 |
try: |
| 14 |
p_t, rhof_t, rhog_t = P.triple_point() |
| 15 |
except RuntimeError,e: |
| 16 |
print "failed to solve triple point" |
| 17 |
sys.exit(1) |
| 18 |
|
| 19 |
pmax = 100e6 |
| 20 |
|
| 21 |
Tmin = P.T_t |
| 22 |
Tmax = 2 * P.T_c |
| 23 |
vmin = 1./rhof_t |
| 24 |
vmax = 2./rhog_t |
| 25 |
TT = linspace(Tmin, Tmax, 100); |
| 26 |
vv = logspace(log10(vmin),log10(vmax), 100); |
| 27 |
|
| 28 |
goodT = [] |
| 29 |
goodv = [] |
| 30 |
badT = [] |
| 31 |
badv = [] |
| 32 |
|
| 33 |
for T in TT: |
| 34 |
sys.stderr.write("+++ T = %f\r" % (T)) |
| 35 |
for v in vv: |
| 36 |
rho = 1./v |
| 37 |
S = P.set_Trho(T,rho) |
| 38 |
p = S.p |
| 39 |
if p > pmax: |
| 40 |
continue |
| 41 |
h = S.h |
| 42 |
#print " p = %f bar, h = %f kJ/kg" % (p/1e5,h/1e3) |
| 43 |
if(h > 8000e3): |
| 44 |
continue |
| 45 |
|
| 46 |
try: |
| 47 |
S = P.set_ph(p,h) |
| 48 |
T1 = S.T |
| 49 |
rho1 = S.rho |
| 50 |
except ValueError,e: |
| 51 |
print "ERROR %s at p = %f, h = %f (T = %.12e, rho = %.12e)" % (str(e),p, h,T,rho) |
| 52 |
badT.append(T); badv.append(v) |
| 53 |
continue |
| 54 |
if isnan(T1) or isnan(rho1): |
| 55 |
print "ERROR at T1 = %f, rho1 = %f (T = %.12e, rho = %.12e)" % (T1, rho1,T,rho) |
| 56 |
badT.append(T); badv.append(v) |
| 57 |
else: |
| 58 |
goodT.append(T); goodv.append(v) |
| 59 |
#print " +++ GOOD RESULT T1 = %f, rho1 = %f" % (T1, rho1) |
| 60 |
|
| 61 |
figure() |
| 62 |
|
| 63 |
print "i \tbad T \tbad v" |
| 64 |
for i in range(len(badT)): |
| 65 |
print "%d\t%e\t%e" % (i,badT[i], badv[i]) |
| 66 |
|
| 67 |
print "TOTAL %d BAD POINTS" % (len(badT)) |
| 68 |
|
| 69 |
print "AXIS =",axis() |
| 70 |
semilogx(badv, badT, 'rx') |
| 71 |
axis([vmin,vmax,Tmin,Tmax]) |
| 72 |
print "AXIS =",axis() |
| 73 |
hold(1) |
| 74 |
semilogx(goodv, goodT, 'g.') |
| 75 |
|
| 76 |
# plot saturation curves |
| 77 |
TTs = linspace(P.T_t, P.T_c, 300) |
| 78 |
TT1 = [] |
| 79 |
vf1 = [] |
| 80 |
vg1 = [] |
| 81 |
for T in TTs: |
| 82 |
try: |
| 83 |
S = P.set_Tx(T,0) |
| 84 |
p = S.p |
| 85 |
rhof = S.rho |
| 86 |
S = P.set_Tx(T,1) |
| 87 |
rhog = S.rho |
| 88 |
except: |
| 89 |
continue; |
| 90 |
TT1.append(T) |
| 91 |
vf1.append(1./rhof) |
| 92 |
vg1.append(1./rhog) |
| 93 |
|
| 94 |
semilogx(vf1,TT1,"b-") |
| 95 |
semilogx(vg1,TT1,"b-") |
| 96 |
axis([vmin,vmax,Tmin,Tmax]) |
| 97 |
title("convergence of (p,h) solver for %s" % P.name) |
| 98 |
xlabel("specific volume") |
| 99 |
ylabel("temperature") |
| 100 |
|
| 101 |
show() |
| 102 |
ion() |
| 103 |
|
| 104 |
|
Parent Directory
| 
Revision Log