/[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 2301 - (show annotations) (download) (as text)
Sat Aug 21 13:28:35 2010 UTC (13 years, 9 months ago) by jpye
File MIME type: text/x-python
File size: 1674 byte(s)
Regen toluene model working, next water.
1 from fprops import *
2 from pylab import *
3 import sys
4
5 D = fprops_fluid('carbondioxide');
6
7 print "SOLVING TRIPLE POINT..."
8
9 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
14 pmax = 100e6
15
16 Tmin = D.T_t
17 Tmax = 2 * D.T_c
18 vmin = 1./rhof_t
19 vmax = 2./rhog_t
20 TT = linspace(Tmin, Tmax, 100);
21 vv = logspace(log10(vmin),log10(vmax), 100);
22
23 goodT = []
24 goodv = []
25 badT = []
26 badv = []
27
28 for T in TT:
29 sys.stderr.write("+++ T = %f\r" % (T))
30 for v in vv:
31 rho = 1./v
32 p = helmholtz_p(T,rho,D)
33 if p > pmax:
34 continue
35 h = helmholtz_h(T,rho,D)
36 #print " p = %f bar, h = %f kJ/kg" % (p/1e5,h/1e3)
37 if(h > 8000e3):
38 continue
39
40 res, T1, rho1 = fprops_solve_ph(p,h,0,D);
41 if res or isnan(T1) or isnan(rho1):
42 print "Error at T1 = %f, rho1 = %f (T = %.12e, rho = %.12e)" % (T1, rho1,T,rho)
43 badT.append(T); badv.append(v)
44 else:
45 goodT.append(T); goodv.append(v)
46 #print " +++ GOOD RESULT T1 = %f, rho1 = %f" % (T1, rho1)
47
48 figure()
49
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 print "AXIS =",axis()
57 semilogx(badv, badT, 'rx')
58 axis([vmin,vmax,Tmin,Tmax])
59 print "AXIS =",axis()
60 hold(1)
61 semilogx(goodv, goodT, 'g.')
62
63 # plot saturation curves
64 TTs = linspace(D.T_t, D.T_c, 300)
65 TT1 = []
66 vf1 = []
67 vg1 = []
68 for T in TTs:
69 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 axis([vmin,vmax,Tmin,Tmax])
78 title("convergence of (p,h) solver for %s" % D.name)
79 xlabel("specific volume")
80 ylabel("temperature")
81
82 show()
83

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