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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2680 - (show annotations) (download) (as text)
Mon Jan 28 06:30:25 2013 UTC (11 years, 5 months ago) by jpye
File MIME type: text/x-python
File size: 1641 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
3 D = fluid("isohexane")
4
5 from pylab import *
6 hold(1)
7
8 T_min = D.T_t
9
10 if T_min == 0:
11 print "WARNING: no triple point specified, using 0.4*T_c"
12 T_min = 0.4 * D.T_c
13
14 TT = linspace(T_min, D.T_c, 1000)
15
16 rhog = array([D.rhog_T_chouaieb(T) for T in TT])
17 rhof = array([D.rhof_T_rackett(T) for T in TT])
18 psat = array([D.psat_T_xiang(T) for T in TT])
19 psata = array([D.psat_T_acentric(T) for T in TT])
20
21 rhof2 = []
22 rhog2 = []
23 psat2 = []
24 TT2 = []
25
26 TT_src = linspace(T_min, D.T_c, 4000)
27 TT1 = []
28 failcount = 0
29 for T in TT_src:
30 try:
31 p1, rf1, rg1 = D.sat_T(T)
32 except Exception,e:
33 print "error '%s' in %s saturation function T = %0.10e " % (str(e),D.name,T)
34 failcount += 1
35 TT1.append(T)
36 continue
37 rhof2.append(rf1)
38 rhog2.append(rg1)
39 psat2.append(p1)
40 TT2.append(T)
41
42 print "failcount =",failcount
43
44 TT = array(TT)
45 TT1 = array(TT1)
46 psat2 = array(psat2)
47 psat = array(psat)
48 psata = array(psata)
49
50 plot(rhog,TT,label="vapour (Chouaieb)")
51 plot(rhof,TT,label="liquid (Rackett)")
52
53 #plot(rhog1,TT1,'rx',label="vapour (unconverged)")
54 #plot(rhof1,TT1,'bx',label="liquid (unconverged)")
55 plot(rhog2,TT2,'r.',label="vapour (OK, converged)")
56 plot(rhof2,TT2,'b.',label="liquid (OK, converged)")
57
58
59 legend(loc=8)
60 xlabel('Density / [kg/m3]')
61 ylabel("Temperature / K");
62 #legend(L)
63 #axis([10,1200,1.,1e6 * D.p_c])
64 #axis([0,1000,0,100e6])
65
66 figure()
67 hold(1)
68
69 plot(TT,psat/1e5,label="Xiang")
70 plot(TT,psata/1e5,label="Acentric")
71 #plot(TT1,psat1/1e5,'rx',label="Maxwell (error)")
72 plot(TT2,psat2/1e5,'g.',label="FPROPS (OK)")
73 legend(loc=2)
74 xlabel("Temperature / K");
75 ylabel(r"$p_\mathrm{sat}(T)$ / [bar]");
76 axis([T_min, D.T_c, 0, D.p_c/1e5])
77 show()
78

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