/[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 2266 - (show annotations) (download) (as text)
Thu Aug 5 13:53:01 2010 UTC (13 years, 10 months ago) by jpye
File MIME type: text/x-python
File size: 1627 byte(s)
The (p,h) function is converging almost everywhere. The home-made Newton iteration wasn't starting well for supercritical states
when the starting guess had rho = rho_c, so changing slightly away from fixed almost all the problems.
1 from fprops import *
2
3 D = helmholtz_data_methane;
4
5 from pylab import *
6 hold(1)
7
8 T_min = D.T_t
9 TT = linspace(T_min, D.T_c, 1000)
10
11 rhog = array([fprops_rhog_T_chouaieb(T,D) for T in TT])
12 rhof = array([fprops_rhof_T_rackett(T,D) for T in TT])
13 psat = array([fprops_psat_T_xiang(T,D) for T in TT])
14 psata = array([fprops_psat_T_acentric(T,D) for T in TT])
15
16 rhof1 = []
17 rhog1 = []
18 psat1 = []
19 rhof2 = []
20 rhog2 = []
21 psat2 = []
22 TT2 = []
23
24 TT_src = linspace(T_min, D.T_c, 4000)
25 TT1 = []
26 failcount = 0
27 for T in TT_src:
28 res, p1, rf1, rg1 = fprops_sat_T(T,D)
29 #print "T=%f, psat=%f bar, rhof=%f, rhog=%f" % (T,p1/1e5,rf1,rg1)
30 if res:
31 print "error %d in %s saturation function T = %0.10e " % (res,D.name,T)
32 failcount += 1
33 rhof1.append(rf1)
34 rhog1.append(rg1)
35 psat1.append(p1)
36 TT1.append(T)
37 continue
38 rhof2.append(rf1)
39 rhog2.append(rg1)
40 psat2.append(p1)
41 TT2.append(T)
42
43 print "failcount =",failcount
44
45 TT = array(TT)
46 TT1 = array(TT1)
47 psat1 = array(psat1)
48 psat2 = array(psat2)
49 psat = array(psat)
50 psata = array(psata)
51
52 plot(rhog,TT,label="vapour (Chouaieb)")
53 plot(rhof,TT,label="liquid (Rackett)")
54
55 plot(rhog1,TT1,'rx',label="vapour (unconverged)")
56 plot(rhof1,TT1,'bx',label="liquid (unconverged)")
57 plot(rhog2,TT2,'r.',label="vapour (OK, converged)")
58 plot(rhof2,TT2,'b.',label="liquid (OK, converged)")
59
60
61 legend(loc=8)
62 xlabel('Density')
63 ylabel('Temperature')
64 #legend(L)
65 #axis([10,1200,1.,1e6 * D.p_c])
66 #axis([0,1000,0,100e6])
67
68 figure()
69 hold(1)
70
71 plot(TT,psat/1e5,label="Xiang")
72 plot(TT,psata/1e5,label="Acentric")
73 plot(TT1,psat1/1e5,'rx',label="Maxwell (error)")
74 plot(TT2,psat2/1e5,'g.',label="Maxwell (OK)")
75 legend(loc=2)
76
77 show()
78

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