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

Working on problem with solve_ph. Could be that one of the deriv routines is wrong in the saturation region?

 1 jpye 2121 from fprops import * 2 3 jpye 2680 D = fluid("isohexane") 4 jpye 2121 5 from pylab import * 6 hold(1) 7 8 jpye 2228 T_min = D.T_t 9 jpye 2678 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 jpye 2124 TT = linspace(T_min, D.T_c, 1000) 15 jpye 2121 16 jpye 2654 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 jpye 2121 21 jpye 2227 rhof2 = [] 22 rhog2 = [] 23 psat2 = [] 24 TT2 = [] 25 jpye 2124 26 jpye 2266 TT_src = linspace(T_min, D.T_c, 4000) 27 jpye 2124 TT1 = [] 28 jpye 2227 failcount = 0 29 for T in TT_src: 30 jpye 2654 try: 31 p1, rf1, rg1 = D.sat_T(T) 32 except Exception,e: 33 jpye 2669 print "error '%s' in %s saturation function T = %0.10e " % (str(e),D.name,T) 34 jpye 2227 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 jpye 2124 42 jpye 2227 print "failcount =",failcount 43 44 jpye 2225 TT = array(TT) 45 TT1 = array(TT1) 46 jpye 2228 psat2 = array(psat2) 47 jpye 2225 psat = array(psat) 48 jpye 2230 psata = array(psata) 49 jpye 2225 50 jpye 2121 plot(rhog,TT,label="vapour (Chouaieb)") 51 plot(rhof,TT,label="liquid (Rackett)") 52 53 jpye 2669 #plot(rhog1,TT1,'rx',label="vapour (unconverged)") 54 #plot(rhof1,TT1,'bx',label="liquid (unconverged)") 55 jpye 2227 plot(rhog2,TT2,'r.',label="vapour (OK, converged)") 56 plot(rhof2,TT2,'b.',label="liquid (OK, converged)") 57 jpye 2124 58 59 jpye 2121 legend(loc=8) 60 jpye 2680 xlabel('Density / [kg/m3]') 61 ylabel("Temperature / K"); 62 jpye 2121 #legend(L) 63 #axis([10,1200,1.,1e6 * D.p_c]) 64 #axis([0,1000,0,100e6]) 65 jpye 2124 66 figure() 67 hold(1) 68 69 jpye 2225 plot(TT,psat/1e5,label="Xiang") 70 jpye 2230 plot(TT,psata/1e5,label="Acentric") 71 jpye 2669 #plot(TT1,psat1/1e5,'rx',label="Maxwell (error)") 72 plot(TT2,psat2/1e5,'g.',label="FPROPS (OK)") 73 jpye 2230 legend(loc=2) 74 jpye 2680 xlabel("Temperature / K"); 75 ylabel(r"$p_\mathrm{sat}(T)$ / [bar]"); 76 jpye 2678 axis([T_min, D.T_c, 0, D.p_c/1e5]) 77 jpye 2121 show() 78