/[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 2680 - (hide annotations) (download) (as text)
Mon Jan 28 06:30:25 2013 UTC (9 years, 6 months ago) by jpye
File MIME type: text/x-python
File size: 2014 byte(s)
Working on problem with solve_ph. Could be that one of the deriv routines is wrong in the saturation region?
1 jpye 2259 from fprops import *
2     from pylab import *
3 jpye 2267 import sys
4 jpye 2259
5 jpye 2662 #P = fluid('water','helmholtz');
6 jpye 2668 #P = fluid('ammonia','pengrob');
7 jpye 2680 P = fluid('isohexane','pengrob');
8 jpye 2259
9 jpye 2287 print "SOLVING TRIPLE POINT..."
10    
11 jpye 2662 print "Fluid: %s\nData source: %s" %(P.name, P.source)
12    
13 jpye 2654 try:
14 jpye 2662 p_t, rhof_t, rhog_t = P.triple_point()
15 jpye 2654 except RuntimeError,e:
16 jpye 2272 print "failed to solve triple point"
17     sys.exit(1)
18 jpye 2267
19     pmax = 100e6
20    
21 jpye 2662 Tmin = P.T_t
22     Tmax = 2 * P.T_c
23 jpye 2268 vmin = 1./rhof_t
24 jpye 2270 vmax = 2./rhog_t
25 jpye 2269 TT = linspace(Tmin, Tmax, 100);
26     vv = logspace(log10(vmin),log10(vmax), 100);
27 jpye 2259
28     goodT = []
29 jpye 2260 goodv = []
30 jpye 2259 badT = []
31 jpye 2260 badv = []
32 jpye 2259
33     for T in TT:
34 jpye 2269 sys.stderr.write("+++ T = %f\r" % (T))
35 jpye 2268 for v in vv:
36     rho = 1./v
37 jpye 2662 S = P.set_Trho(T,rho)
38     p = S.p
39 jpye 2267 if p > pmax:
40     continue
41 jpye 2662 h = S.h
42 jpye 2267 #print " p = %f bar, h = %f kJ/kg" % (p/1e5,h/1e3)
43 jpye 2259 if(h > 8000e3):
44     continue
45    
46 jpye 2654 try:
47 jpye 2662 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 jpye 2654 badT.append(T); badv.append(v)
53     continue
54     if isnan(T1) or isnan(rho1):
55 jpye 2662 print "ERROR at T1 = %f, rho1 = %f (T = %.12e, rho = %.12e)" % (T1, rho1,T,rho)
56 jpye 2268 badT.append(T); badv.append(v)
57 jpye 2259 else:
58 jpye 2268 goodT.append(T); goodv.append(v)
59 jpye 2267 #print " +++ GOOD RESULT T1 = %f, rho1 = %f" % (T1, rho1)
60 jpye 2259
61     figure()
62 jpye 2267
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 jpye 2268 print "AXIS =",axis()
70 jpye 2260 semilogx(badv, badT, 'rx')
71 jpye 2268 axis([vmin,vmax,Tmin,Tmax])
72     print "AXIS =",axis()
73 jpye 2259 hold(1)
74 jpye 2260 semilogx(goodv, goodT, 'g.')
75    
76     # plot saturation curves
77 jpye 2662 TTs = linspace(P.T_t, P.T_c, 300)
78 jpye 2260 TT1 = []
79     vf1 = []
80     vg1 = []
81 jpye 2268 for T in TTs:
82 jpye 2654 try:
83 jpye 2662 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 jpye 2654 except:
89     continue;
90     TT1.append(T)
91     vf1.append(1./rhof)
92     vg1.append(1./rhog)
93 jpye 2260
94     semilogx(vf1,TT1,"b-")
95     semilogx(vg1,TT1,"b-")
96 jpye 2268 axis([vmin,vmax,Tmin,Tmax])
97 jpye 2662 title("convergence of (p,h) solver for %s" % P.name)
98 jpye 2266 xlabel("specific volume")
99     ylabel("temperature")
100 jpye 2260
101 jpye 2259 show()
102 jpye 2680 ion()
103 jpye 2259
104 jpye 2662

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