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

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