from __future__ import division from numpy import array,arange,sum,ones,sqrt # Taken from # The International Association for the Properties of Water and # Steam. Lucerne, Switzerland. August 2007. Revised Release on the # IAPWS Industrial Formulation 1997 for the Thermodynamic Properties # of Water and Steam. ## Functions implemented # Saturation line # h(p,T) region1, region2 and no warnings yet # constants R=0.461526 #kJ/(kg K) Tc=647.096 #Triple point temperature (K) pc=22.064 #Triple point pressure (MPa) rc=322 #Triple point density kg/m3 def psat(T): """ Returns the saturation pressure of water at a given temperature. Remember that temperature must be between 273.15K (triple point) and 647.096K (critical point) Temperatures in K, Pressures in MPa >>> psat(300) 0.0035365894130130106 >>> psat(130) 'Error: No saturation pressure for this temperature' >>> psat(700) 'Error: No saturation pressure for this temperature' """ if T < 273.15 or T > 647.096: return 'Error: No saturation pressure for this temperature' n=array([0.11670521452767E4, - 0.72421316703206E6, - 0.17073846940092E2, 0.12020824702470E5 , - 0.32325550322333E7, 0.14915108613530E2 , - 0.48232657361591E4, 0.40511340542057E6 , - 0.23855557567849 , 0.65017534844798E3],'d') v=T+n[8]/(T-n[9]) A=v**2+n[0]*v+n[1] B=n[2]*v**2+n[3]*v+n[4] C=n[5]*v**2+n[6]*v+n[7] return ((2*C)/(-B+sqrt(B**2-4*A*C)))**4 def tsat(p): """ Returns the saturation temperature of water at a given pressure. Remember that pressure must be between 0.000611213MPa (triple point) and 22.064MPa (critical point) Temperatures in K, pressures in MPa >>> tsat(0.1) 372.75591861133762 >>> tsat(0.00001) 'Error: No saturation temperature for this pressure' >>> tsat(100) 'Error: No saturation temperature for this pressure' """ if p < 0.000611213 or p > 22.064: return 'Error: No saturation temperature for this pressure' n=array([0.11670521452767E4, - 0.72421316703206E6, - 0.17073846940092E2, 0.12020824702470E5 , - 0.32325550322333E7, 0.14915108613530E2 , - 0.48232657361591E4, 0.40511340542057E6 , - 0.23855557567849 , 0.65017534844798E3],'d') beta=p**0.25 E=beta**2+n[2]*beta+n[5] F=n[0]*beta**2+n[3]*beta+n[6] G=n[1]*beta**2+n[4]*beta+n[7] D=(2*G)/(-F-sqrt(F**2-4*E*G)) return 0.5*(n[9]+D-sqrt((n[9]+D)**2-4*(n[8]+n[9]*D))) def h_pT(p,T): """ Returns specific enthalpy (kJ/kg) for a given pressure (MPa) and Temperature (K) >>> h_pT(3,300) 115.33127302143839 >>> h_pT(0.0035,300) 2549.9114508400203 Results checked against the reference. """ # region 1 implementation if p >= psat(T) : #Liquid water (pressure over saturatio pressure) pi=p/16.53 tau=1386/T raw_data=array([ [0 ,-2 , 0.14632971213167 ], [0 ,-1 ,-0.84548187169114 ], [0 , 0 ,-0.37563603672040E1 ], [0 , 1 , 0.33855169168385E1 ], [0 , 2 ,-0.95791963387872 ], [0 , 3 , 0.15772038513228 ], [0 , 4 ,-0.16616417199501E-1 ], [0 , 5 , 0.81214629983568E-3 ], [1 ,-9 , 0.28319080123804E-3 ], [1 ,-7 ,-0.60706301565874E-3 ], [1 ,-1 ,-0.18990068218419E-1 ], [1 , 0 ,-0.32529748770505E-1 ], [1 , 1 ,-0.21841717175414E-1 ], [1 , 3 ,-0.52838357969930E-4 ], [2 ,-3 ,-0.47184321073267E-3 ], [2 , 0 ,-0.30001780793026E-3 ], [2 , 1 , 0.47661393906987E-4 ], [2 , 3 ,-0.44141845330846E-5 ], [2 , 17,-0.72694996297594E-15], [3 ,-4 ,-0.31679644845054E-4 ], [3 , 0 ,-0.28270797985312E-5 ], [3 , 6 ,-0.85205128120103E-9 ], [4 ,-5 ,-0.22425281908000E-5 ], [4 ,-2 ,-0.65171222895601E-6 ], [4 , 10,-0.14341729937924E-12], [5 ,-8 ,-0.40516996860117E-6 ], [8 ,-11,-0.12734301741641E-8 ], [8 ,-6 ,-0.17424871230634E-9 ], [21,-29,-0.68762131295531E-18], [23,-31, 0.14478307828521E-19], [29,-38, 0.26335781662795E-22], [30,-39,-0.11947622640071E-22], [31,-40, 0.18228094581404E-23], [32,-41,-0.93537087292458E-25]],'d') I=raw_data[:,0] J=raw_data[:,1] n=raw_data[:,2] return R*T*tau*sum((n*(7.1-pi)**I)*J*((tau-1.222)**(J-1))) if p < psat(T): # steam, pressure under saturation pressure pi=p tau=540/T raw_data0=array([ [0 , -0.96927686500217E1 ], [1 , 0.10086655968018E2 ], [-5, -0.56087911283020E-2], [-4, 0.71452738081455E-1 ], [-3, -0.40710498223928 ], [-2, 0.14240819171444E1 ], [-1, -0.43839511319450E1 ], [2 , -0.28408632460772 ], [3 , 0.21268463753307E-1 ]],'d') J0=raw_data0[:,0] n0=raw_data0[:,1] raw_datar=array([ [1 ,0 ,-0.17731742473213E-2], [1 ,1 ,-0.17834862292358E-1], [1 ,2 ,-0.45996013696365E-1], [1 ,3 ,-0.57581259083432E-1], [1 ,6 ,-0.50325278727930E-1], [2 ,1 ,-0.33032641670203E-4], [2 ,2 ,-0.18948987516315E-3], [2 ,4 ,-0.39392777243355E-2], [2 ,7 ,-0.43797295650573E-1], [2 ,36,-0.26674547914087E-4], [3 ,0 , 0.20481737692309E-7], [3 ,1 , 0.43870667284435E-6], [3 ,3 ,-0.32277677238570E-4], [3 ,6 ,-0.15033924542148E-2], [3 ,35,-0.40668253562649E-1], [4 ,1 ,-0.78847309559367E-9], [4 ,2 , 0.12790717852285E-7], [4 ,3 , 0.48225372718507E-6], [5 ,7 , 0.22922076337661E-5], [6 ,3 ,-0.16714766451061E-1], [6 ,16,-0.21171472321355E-2], [6 ,35,-0.23895741934104E2 ], [7 ,0 ,-0.59059564324270E-1], [7 ,11,-0.12621808899101E-5], [7 ,25,-0.38946842435739E-1], [8 ,8 , 0.11256211360459E-1], [8 ,36,-0.82311340897998E1 ], [9 ,13, 0.19809712802088E-7], [10,4 , 0.10406965210174E-1], [10,10,-0.10234747095929E-1], [10,14,-0.10018179379511E-8], [16,29,-0.80882908646985E-1], [16,50, 0.10693031879409 ], [18,57,-0.33662250574171 ], [20,20, 0.89185845355421E-2], [20,35, 0.30629316876232E-1], [20,48,-0.42002467698208E-5], [21,21,-0.59056029685639E-2], [22,53, 0.37826947613457E-5], [23,39,-0.12768608934681E-1], [24,26, 0.73087610595061E-2], [24,40, 0.55414715350778E-1], [24,58,-0.94369707241210E-6]],'d') I=raw_datar[:,0] J=raw_datar[:,1] n=raw_datar[:,2] g0_tau=sum(n0*J0*tau**(J0-1)) gr_tau=sum(n*pi**I*J*(tau-0.5)**(J-1)) return R*T*tau*(g0_tau+gr_tau) def _test(): import doctest doctest.testmod() if __name__ == '__main__': _test()