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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2301 - (hide annotations) (download) (as text)
Sat Aug 21 13:28:35 2010 UTC (13 years, 10 months ago) by jpye
File MIME type: text/x-python
File size: 1674 byte(s)
Regen toluene model working, next water.
1 jpye 2259 from fprops import *
2     from pylab import *
3 jpye 2267 import sys
4 jpye 2259
5 jpye 2297 D = fprops_fluid('carbondioxide');
6 jpye 2259
7 jpye 2287 print "SOLVING TRIPLE POINT..."
8    
9 jpye 2272 res, p_t, rhof_t, rhog_t = fprops_triple_point(D)
10     if res:
11     print "failed to solve triple point"
12     sys.exit(1)
13 jpye 2267
14     pmax = 100e6
15    
16 jpye 2270 Tmin = D.T_t
17     Tmax = 2 * D.T_c
18 jpye 2268 vmin = 1./rhof_t
19 jpye 2270 vmax = 2./rhog_t
20 jpye 2269 TT = linspace(Tmin, Tmax, 100);
21     vv = logspace(log10(vmin),log10(vmax), 100);
22 jpye 2259
23     goodT = []
24 jpye 2260 goodv = []
25 jpye 2259 badT = []
26 jpye 2260 badv = []
27 jpye 2259
28     for T in TT:
29 jpye 2269 sys.stderr.write("+++ T = %f\r" % (T))
30 jpye 2268 for v in vv:
31     rho = 1./v
32 jpye 2259 p = helmholtz_p(T,rho,D)
33 jpye 2267 if p > pmax:
34     continue
35 jpye 2259 h = helmholtz_h(T,rho,D)
36 jpye 2267 #print " p = %f bar, h = %f kJ/kg" % (p/1e5,h/1e3)
37 jpye 2259 if(h > 8000e3):
38     continue
39    
40     res, T1, rho1 = fprops_solve_ph(p,h,0,D);
41 jpye 2268 if res or isnan(T1) or isnan(rho1):
42 jpye 2269 print "Error at T1 = %f, rho1 = %f (T = %.12e, rho = %.12e)" % (T1, rho1,T,rho)
43 jpye 2268 badT.append(T); badv.append(v)
44 jpye 2259 else:
45 jpye 2268 goodT.append(T); goodv.append(v)
46 jpye 2267 #print " +++ GOOD RESULT T1 = %f, rho1 = %f" % (T1, rho1)
47 jpye 2259
48     figure()
49 jpye 2267
50     print "i \tbad T \tbad v"
51     for i in range(len(badT)):
52     print "%d\t%e\t%e" % (i,badT[i], badv[i])
53    
54     print "TOTAL %d BAD POINTS" % (len(badT))
55    
56 jpye 2268 print "AXIS =",axis()
57 jpye 2260 semilogx(badv, badT, 'rx')
58 jpye 2268 axis([vmin,vmax,Tmin,Tmax])
59     print "AXIS =",axis()
60 jpye 2259 hold(1)
61 jpye 2260 semilogx(goodv, goodT, 'g.')
62    
63     # plot saturation curves
64 jpye 2268 TTs = linspace(D.T_t, D.T_c, 300)
65 jpye 2260 TT1 = []
66     vf1 = []
67     vg1 = []
68 jpye 2268 for T in TTs:
69 jpye 2260 res, p, rhof, rhog = fprops_sat_T(T,D)
70     if not res:
71     TT1.append(T)
72     vf1.append(1./rhof)
73     vg1.append(1./rhog)
74    
75     semilogx(vf1,TT1,"b-")
76     semilogx(vg1,TT1,"b-")
77 jpye 2268 axis([vmin,vmax,Tmin,Tmax])
78 jpye 2301 title("convergence of (p,h) solver for %s" % D.name)
79 jpye 2266 xlabel("specific volume")
80     ylabel("temperature")
81 jpye 2260
82 jpye 2259 show()
83    

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