/[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 2680 - (show annotations) (download) (as text)
Mon Jan 28 06:30:25 2013 UTC (9 years, 4 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 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

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