/[ascend]/trunk/models/johnpye/fprops/python/solve_ph_array.py
ViewVC logotype

Contents of /trunk/models/johnpye/fprops/python/solve_ph_array.py

Parent Directory Parent Directory | Revision Log Revision Log


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

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22