?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Correction d'une faute d'orthographe
[memoirecycle.git] / XSteam.py
blobabc97fae271f92e9d1b5abd63515f15e5c07e283
1 """
2 %***********************************************************************************************************
3 %* Water and steam properties according to IAPWS IF-97 *
4 %* By Magnus Holmgren, www.x-eng.com *
5 %* The steam tables are free and provided as is. *
6 %* We take no responsibilities for any errors in the code or damage thereby. *
7 %* You are free to use, modify and distribute the code as long as authorship is properly acknowledged. *
8 %* Please notify me at magnus@x-eng.com if the code is used in commercial applications *
9 %***********************************************************************************************************
10 %Translated to Python and check with check() by Gauthier Zeimes
11 %
12 % XSteam provides accurate steam and water properties from 0 - 1000 bar and from 0 - 2000 deg C according to
13 % the standard IAPWS IF-97. For accuracy of the functions in different regions see IF-97 (www.iapws.org)
14 %
15 % *** Using XSteam *****************************************************************************************
16 %XSteam take 2 or 3 arguments. The first argument must always be the steam table function you want to use.
17 %The other arguments are the inputs to that function.
18 %Example: XSteam('h_pt',1,20) Returns the enthalpy of water at 1 bar and 20 degC
19 %Example: XSteam('TSat_p',1) Returns the saturation temperature of water at 1 bar.
20 %For a list of valid Steam Table functions se bellow or the XSteam macros for MS Excel.
21 %
22 %*** Nomenclature ******************************************************************************************
23 % First the wanted property then a _ then the wanted input properties.
24 % Example. T_ph is temperature as a function of pressure and enthalpy.
25 % For a list of valid functions se bellow or XSteam for MS Excel.
26 % T Temperature (deg C)
27 % p Pressure (bar)
28 % h Enthalpy (kJ/kg)
29 % v Specific volume (m3/kg)
30 % rho Density
31 % s Specific entropy
32 % u Specific internal energy
33 % Cp Specific isobaric heat capacity
34 % Cv Specific isochoric heat capacity
35 % w Speed of sound
36 % my Viscosity
37 % tc Thermal Conductivity
38 % st Surface Tension
39 % x Vapour fraction
40 % vx Vapour Volume Fraction
41 %
42 %*** Valid Steam table functions. ****************************************************************************
43 %
44 %Temperature
45 %Tsat_p Saturation temperature
46 %T_ph Temperture as a function of pressure and enthalpy
47 %T_ps Temperture as a function of pressure and entropy
48 %T_hs Temperture as a function of enthalpy and entropy
49 %
50 %Pressure
51 %psat_T Saturation pressure
52 %p_hs Pressure as a function of h and s.
53 %p_hrho Pressure as a function of h and rho. Very unaccurate for solid water region since it's almost incompressible!
54 %
55 %Enthalpy
56 %hV_p Saturated vapour enthalpy
57 %hL_p Saturated liquid enthalpy
58 %hV_T Saturated vapour enthalpy
59 %hL_T Saturated liquid enthalpy
60 %h_pT Entalpy as a function of pressure and temperature.
61 %h_ps Entalpy as a function of pressure and entropy.
62 %h_px Entalpy as a function of pressure and vapour fraction
63 %h_prho Entalpy as a function of pressure and density. Observe for low temperatures (liquid) this equation has 2 solutions.
64 %h_Tx Entalpy as a function of temperature and vapour fraction
65 %
66 %Specific volume
67 %vV_p Saturated vapour volume
68 %vL_p Saturated liquid volume
69 %vV_T Saturated vapour volume
70 %vL_T Saturated liquid volume
71 %v_pT Specific volume as a function of pressure and temperature.
72 %v_ph Specific volume as a function of pressure and enthalpy
73 %v_ps Specific volume as a function of pressure and entropy.
74 %
75 %Density
76 %rhoV_p Saturated vapour density
77 %rhoL_p Saturated liquid density
78 %rhoV_T Saturated vapour density
79 %rhoL_T Saturated liquid density
80 %rho_pT Density as a function of pressure and temperature.
81 %rho_ph Density as a function of pressure and enthalpy
82 %rho_ps Density as a function of pressure and entropy.
83 %
84 %Specific entropy
85 %sV_p Saturated vapour entropy
86 %sL_p Saturated liquid entropy
87 %sV_T Saturated vapour entropy
88 %sL_T Saturated liquid entropy
89 %s_pT Specific entropy as a function of pressure and temperature (Returns saturated vapour entalpy if mixture.)
90 %s_ph Specific entropy as a function of pressure and enthalpy
91 %
92 %Specific internal energy
93 %uV_p Saturated vapour internal energy
94 %uL_p Saturated liquid internal energy
95 %uV_T Saturated vapour internal energy
96 %uL_T Saturated liquid internal energy
97 %u_pT Specific internal energy as a function of pressure and temperature.
98 %u_ph Specific internal energy as a function of pressure and enthalpy
99 %u_ps Specific internal energy as a function of pressure and entropy.
100 %
101 %Specific isobaric heat capacity
102 %CpV_p Saturated vapour heat capacity
103 %CpL_p Saturated liquid heat capacity
104 %CpV_T Saturated vapour heat capacity
105 %CpL_T Saturated liquid heat capacity
106 %Cp_pT Specific isobaric heat capacity as a function of pressure and temperature.
107 %Cp_ph Specific isobaric heat capacity as a function of pressure and enthalpy
108 %Cp_ps Specific isobaric heat capacity as a function of pressure and entropy.
109 %
110 %Specific isochoric heat capacity
111 %CvV_p Saturated vapour isochoric heat capacity
112 %CvL_p Saturated liquid isochoric heat capacity
113 %CvV_T Saturated vapour isochoric heat capacity
114 %CvL_T Saturated liquid isochoric heat capacity
115 %Cv_pT Specific isochoric heat capacity as a function of pressure and temperature.
116 %Cv_ph Specific isochoric heat capacity as a function of pressure and enthalpy
117 %Cv_ps Specific isochoric heat capacity as a function of pressure and entropy.
118 %
119 %Speed of sound
120 %wV_p Saturated vapour speed of sound
121 %wL_p Saturated liquid speed of sound
122 %wV_T Saturated vapour speed of sound
123 %wL_T Saturated liquid speed of sound
124 %w_pT Speed of sound as a function of pressure and temperature.
125 %w_ph Speed of sound as a function of pressure and enthalpy
126 %w_ps Speed of sound as a function of pressure and entropy.
127 %
128 %Viscosity
129 %Viscosity is not part of IAPWS Steam IF97. Equations from
130 %"Revised Release on the IAPWS Formulation 1985 for the Viscosity of Ordinary Water Substance", 2003 are used.
131 %Viscosity in the mixed region (4) is interpolated according to the density. This is not true since it will be two fases.
132 %my_pT Viscosity as a function of pressure and temperature.
133 %my_ph Viscosity as a function of pressure and enthalpy
134 %my_ps Viscosity as a function of pressure and entropy.
135 %
136 %Thermal Conductivity
137 %Revised release on the IAPS Formulation 1985 for the Thermal Conductivity of ordinary water substance (IAPWS 1998)
138 %tcL_p Saturated vapour thermal conductivity
139 %tcV_p Saturated liquid thermal conductivity
140 %tcL_T Saturated vapour thermal conductivity
141 %tcV_T Saturated liquid thermal conductivity
142 %tc_pT Thermal conductivity as a function of pressure and temperature.
143 %tc_ph Thermal conductivity as a function of pressure and enthalpy
144 %tc_hs Thermal conductivity as a function of enthalpy and entropy
145 %
146 %Surface tension
147 %st_T Surface tension for two phase water/steam as a function of T
148 %st_p Surface tension for two phase water/steam as a function of T
149 %Vapour fraction
150 %x_ph Vapour fraction as a function of pressure and enthalpy
151 %x_ps Vapour fraction as a function of pressure and entropy.
152 %
153 %Vapour volume fraction
154 %vx_ph Vapour volume fraction as a function of pressure and enthalpy
155 %vx_ps Vapour volume fraction as a function of pressure and entropy.
156 """
157
158 """
159 function Out=XSteam(fun,In1,In2)
160 %*Contents.
161 %*1 Calling functions
162 %*1.1
163 %*1.2 Temperature (T)
164 %*1.3 Pressure (p)
165 %*1.4 Enthalpy (h)
166 %*1.5 Specific Volume (v)
167 %*1.7 Specific entropy (s)
168 %*1.17 Vapour fraction
169
170 NOT TRANSLATED:
171 %*1.6 Density (rho)
172 %*1.8 Specific internal energy (u)
173 %*1.9 Specific isobaric heat capacity (Cp)
174 %*1.10 Specific isochoric heat capacity (Cv)
175 %*1.11 Speed of sound
176 %*1.12 Viscosity
177 %*1.13 Prandtl
178 %*1.14 Kappa
179 %*1.15 Surface tension
180 %*1.16 Heat conductivity
181 %*1.18 Vapour Volume Fraction
182
183 %
184 %*2 IAPWS IF 97 Calling functions
185 %*2.1 Functions for region 1
186 %*2.2 Functions for region 2
187 %*2.3 Functions for region 3
188 %*2.4 Functions for region 4
189 %*2.5 Functions for region 5
190 %
191 %*3 Region Selection
192 %*3.1 Regions as a function of pT
193 %*3.2 Regions as a function of ph
194 %*3.3 Regions as a function of ps
195 %*3.4 Regions as a function of hs
196 %
197 %4 Region Borders
198 %4.1 Boundary between region 1 and 3.
199 %4.2 Region 3. pSat_h and pSat_s
200 %4.3 Region boundary 1to3 and 3to2 as a functions of s
201 %
202 %5 Transport properties
203 %5.1 Viscosity (IAPWS formulation 1985)
204 %5.2 Thermal Conductivity (IAPWS formulation 1985)
205 %5.3 Surface Tension
206 %
207 %6 Units
208 %
209 %7 Verification
210 %7.1 Verifiy region 1
211 %7.2 Verifiy region 2
212 %7.3 Verifiy region 3
213 %7.4 Verifiy region 4
214 %7.5 Verifiy region 5
215 NOT DONE 7.6
216 """
217
218
219 import numpy
220 import math
221
222 def log(x):
223 return math.log(abs(x))
224
225 '''
226 ***********************************************************************************************************
227 *1 Calling functions *
228 ***********************************************************************************************************
229 '''
230
231 def function(fun,In1,In2):
232 fun = fun.lower()
233 #***************************************************************************
234 #1.2 Temperature
235 #***************************************************************************
236 if fun == 'tsat_p':
237 p = toSIunit_p(In1)
238 if p > 0.000611657 and p < 22.06395 :
239 return fromSIunit_T(T4_p(p))
240 else:
241 return float('NaN')
242
243 elif fun == 'tsat_s':
244 s = toSIunit_s(In1)
245 if s > -0.0001545495919 and s < 9.155759395 :
246 ps = p4_s(s)
247 return fromSIunit_T(T4_p(ps))
248 else:
249 return float('NaN')
250
251 elif fun == 't_ph':
252 p = toSIunit_p(In1)
253 h = toSIunit_h(In2)
254 Region = region_ph(p, h)
255 if Region == 1:
256 return fromSIunit_T(T1_ph(p, h))
257 elif Region == 2:
258 return fromSIunit_T(T2_ph(p, h))
259 elif Region == 3:
260 return fromSIunit_T(T3_ph(p, h))
261 elif Region == 4:
262 return fromSIunit_T(T4_p(p))
263 elif Region == 5:
264 return fromSIunit_T(T5_ph(p, h))
265 else:
266 return float('NaN')
267
268 elif fun == 't_ps':
269 p = toSIunit_p(In1)
270 s = toSIunit_s(In2)
271 Region = region_ps(p, s)
272 if Region == 1:
273 return fromSIunit_T(T1_ps(p, s))
274 elif Region == 2:
275 return fromSIunit_T(T2_ps(p, s))
276 elif Region == 3:
277 return fromSIunit_T(T3_ps(p, s))
278 elif Region == 4:
279 return fromSIunit_T(T4_p(p))
280 elif Region == 5:
281 return fromSIunit_T(T5_ps(p, s))
282 else:
283 return float('NaN')
284
285 elif fun == 't_hs':
286 h = toSIunit_h(In1)
287 s = toSIunit_s(In2)
288 Region = region_hs(h, s)
289 if Region == 1:
290 p1 = p1_hs(h, s)
291 return fromSIunit_T(T1_ph(p1, h))
292 elif Region == 2:
293 p2 = p2_hs(h, s);
294 return fromSIunit_T(T2_ph(p2, h))
295 elif Region == 3:
296 p3 = p3_hs(h, s);
297 return fromSIunit_T(T3_ph(p3, h))
298 elif Region == 4:
299 return fromSIunit_T(T4_hs(h, s))
300 elif Region == 5:
301 print 'functions of hs is not avlaible in region 5'
302 else:
303 return float('NaN')
304
305 #***************************************************************************
306 #1.3 Pressure (p)
307 #***************************************************************************
308 elif fun == 'psat_t':
309 T = toSIunit_T(In1)
310 if T < 647.096 and T > 273.15 :
311 return fromSIunit_p(p4_T(T))
312 else:
313 return float('NaN')
314 elif fun == 'psat_s':
315 s = toSIunit_s(In1)
316 if s > -0.0001545495919 and s < 9.155759395 :
317 return fromSIunit_p(p4_s(s))
318 else:
319 return float('NaN')
320
321 elif fun == 'p_hs':
322 h = toSIunit_h(In1)
323 s = toSIunit_s(In2)
324 Region = region_hs(h,s)
325 if Region == 1:
326 return fromSIunit_p(p1_hs(h, s))
327 elif Region == 2:
328 return fromSIunit_p(p2_hs(h, s))
329 elif Region == 3:
330 return fromSIunit_p(p3_hs(h, s))
331 elif Region == 4:
332 tSat = T4_hs(h, s)
333 return fromSIunit_p(p4_T(tSat))
334 elif Region == 5:
335 print 'functions of hs is not avlaible in region 5'
336 else:
337 return float('NaN')
338
339 elif fun == 'p_hrho':
340 h=In1
341 rho=In2
342 print 'Function p_hrho not implemented'
343 '''
344 %Not valid for water or sumpercritical since water rho does not change very much with p.
345 %Uses iteration to find p.
346 High_Bound = fromSIunit_p(100);
347 Low_Bound = fromSIunit_p(0.000611657);
348 ps = fromSIunit_p(10);
349 rhos = 1 / XSteam('v_ph',ps, h);
350 while abs(rho - rhos) > 0.0000001
351 rhos = 1 / XSteam('v_ph',ps, h);
352 if rhos >= rho
353 High_Bound = ps;
354 else
355 Low_Bound = ps;
356 end
357 ps = (Low_Bound + High_Bound) / 2.0;
358 end
359 Out = ps;
360 '''
361
362
363 #***************************************************************************
364 #1.4 Enthalpy (h)
365 #***************************************************************************
366
367 elif fun == 'hv_p':
368 p = toSIunit_p(In1)
369 if p > 0.000611657 and p < 22.06395 :
370 return fromSIunit_h(h4V_p(p))
371 else:
372 return float('NaN')
373
374 elif fun == 'hl_p':
375 p = toSIunit_p(In1)
376 if p > 0.000611657 and p < 22.06395 :
377 return fromSIunit_h(h4L_p(p))
378 else:
379 return float('NaN')
380
381 elif fun == 'hv_t':
382 T = toSIunit_T(In1)
383 if T > 273.15 and T < 647.096 :
384 p = p4_T(T)
385 return fromSIunit_h(h4V_p(p))
386 else:
387 return float('NaN')
388
389 elif fun == 'hl_t':
390 T = toSIunit_T(In1)
391 if T > 273.15 and T < 647.096 :
392 p = p4_T(T)
393 return fromSIunit_h(h4L_p(p))
394 else:
395 return float('NaN')
396
397 elif fun == 'h_pt':
398 p = toSIunit_p(In1)
399 T = toSIunit_T(In2)
400 Region = region_pT(p, T)
401 if Region == 1 :
402 return fromSIunit_h(h1_pT(p, T))
403 elif Region == 2 :
404 return fromSIunit_h(h2_pT(p, T))
405 elif Region == 3 :
406 return fromSIunit_h(h3_pT(p, T))
407 elif Region == 4 :
408 return float('NaN')
409 elif Region == 5 :
410 return fromSIunit_h(h5_pT(p, T))
411 else:
412 return float('NaN')
413
414 elif fun == 'h_ps':
415 p = toSIunit_p(In1)
416 s = toSIunit_s(In2)
417 Region = region_ps(p, s)
418 if Region == 1 :
419 return fromSIunit_h(h1_pT(p, T1_ps(p, s)))
420 elif Region == 2 :
421 return fromSIunit_h(h2_pT(p, T2_ps(p, s)))
422 elif Region == 3 :
423 return fromSIunit_h(h3_rhoT(1 / v3_ps(p, s), T3_ps(p, s)))
424 elif Region == 4 :
425 xs = x4_ps(p, s)
426 return fromSIunit_h(xs * h4V_p(p) + (1 - xs) * h4L_p(p))
427 elif Region == 5 :
428 return fromSIunit_h(h5_pT(p, T5_ps(p, s)))
429 else:
430 return float('NaN')
431
432 elif fun == 'h_px':
433 p = toSIunit_p(In1)
434 x = toSIunit_x(In2)
435 if x > 1 or x < 0 or p >= 22.064 :
436 return float('NaN')
437 else:
438 hL = h4L_p(p)
439 hV = h4V_p(p)
440 return hL + x * (hV - hL)
441
442 elif fun == 'h_prho':
443 p = toSIunit_p(In1)
444 rho = 1.0 / toSIunit_v(1.0 / In2)
445 Region = Region_prho(p, rho)
446 if Region == 1:
447 return fromSIunit_h(h1_pT(p, T1_prho(p, rho)))
448 elif Region == 2:
449 return fromSIunit_h(h2_pT(p, T2_prho(p, rho)))
450 elif Region == 3:
451 return fromSIunit_h(h3_rhoT(rho, T3_prho(p, rho)))
452 elif Region == 4:
453 if p < 16.529:
454 vV = v2_pT(p, T4_p(p))
455 vL = v1_pT(p, T4_p(p))
456 else:
457 vV = v3_ph(p, h4V_p(p))
458 vL = v3_ph(p, h4L_p(p))
459 hV = h4V_p(p)
460 hL = h4L_p(p)
461 x = (1 / rho - vL) / (vV - vL)
462 return fromSIunit_h((1 - x) * hL + x * hV)
463 elif Region == 5:
464 return fromSIunit_h(h5_pT(p, T5_prho(p, rho)))
465 else :
466 return float('NaN')
467
468
469 elif fun == 'h_tx':
470 T = toSIunit_T(In1)
471 x = toSIunit_x(In2)
472 if x > 1 or x < 0 or T >= 647.096 :
473 return float('NaN')
474 else :
475 p = p4_T(T)
476 hL = h4L_p(p)
477 hV = h4V_p(p)
478 return hL + x * (hV - hL)
479
480 #***************************************************************************
481 #1.7 Specific volume (v)
482 #***************************************************************************
483 elif fun == 'vv_p':
484 p = toSIunit_p(In1)
485 if p > 0.000611657 and p < 22.06395 :
486 if p < 16.529 :
487 return fromSIunit_v(v2_pT(p, T4_p(p)))
488 else :
489 return fromSIunit_v(v3_ph(p, h4V_p(p)))
490 else:
491 return float('NaN')
492
493 elif fun == 'vl_p':
494 p = toSIunit_p(In1)
495 if p > 0.000611657 and p < 22.06395 :
496 if p < 16.529 :
497 return fromSIunit_v(v1_pT(p, T4_p(p)))
498 else:
499 return fromSIunit_v(v3_ph(p, h4L_p(p)))
500 else:
501 return float('NaN')
502
503 elif fun == 'vv_t':
504 T = toSIunit_T(In1)
505 if T > 273.15 and T < 647.096 :
506 if T <= 623.15 :
507 return fromSIunit_v(v2_pT(p4_T(T), T))
508 else:
509 return fromSIunit_v(v3_ph(p4_T(T), h4V_p(p4_T(T))))
510 else:
511 return None
512
513 elif fun == 'vl_t':
514 T = toSIunit_T(In1)
515 if T > 273.15 and T < 647.096 :
516 if T <= 623.15 :
517 return fromSIunit_v(v1_pT(p4_T(T), T))
518 else :
519 return fromSIunit_v(v3_ph(p4_T(T), h4L_p(p4_T(T))))
520 else :
521 return float('NaN')
522
523 elif fun == 'v_pt':
524 p = toSIunit_p(In1)
525 T = toSIunit_T(In2)
526 Region = region_pT(p, T)
527 if Region == 1:
528 return fromSIunit_v(v1_pT(p, T))
529 elif Region == 2:
530 return fromSIunit_v(v2_pT(p, T))
531 elif Region == 3:
532 return fromSIunit_v(v3_ph(p, h3_pT(p, T)))
533 elif Region == 4:
534 return None
535 elif Region == 5:
536 return fromSIunit_v(v5_pT(p, T))
537 else:
538 return float('NaN')
539
540 elif fun == 'v_ph':
541 p = toSIunit_p(In1)
542 h = toSIunit_h(In2)
543 Region = region_ph(p, h)
544 if Region == 1:
545 return fromSIunit_v(v1_pT(p, T1_ph(p, h)))
546 elif Region == 2:
547 return fromSIunit_v(v2_pT(p, T2_ph(p, h)))
548 elif Region == 3:
549 return fromSIunit_v(v3_ph(p, h))
550 elif Region == 4:
551 xs = x4_ph(p,h)
552 if p < 16.529:
553 v4v = v2_pT(p, T4_p(p))
554 v4L = v1_pT(p, T4_p(p))
555 else:
556 v4v = v3_ph(p, h4V_p(p))
557 v4L = v3_ph(p, h4L_p(p))
558 return fromSIunit_v((xs * v4v + (1 - xs) * v4L))
559 elif Region == 5:
560 Ts = T5_ph(p, h)
561 return fromSIunit_v(v5_pT(p, Ts))
562 else:
563 return float('NaN')
564
565
566 elif fun == 'v_ps':
567 p = toSIunit_p(In1)
568 s = toSIunit_s(In2)
569 Region = region_ps(p, s)
570 if Region == 1:
571 return fromSIunit_v(v1_pT(p, T1_ps(p, s)))
572 elif Region == 2:
573 return fromSIunit_v(v2_pT(p, T2_ps(p, s)))
574 elif Region == 3:
575 return fromSIunit_v(v3_ps(p, s))
576 elif Region == 4:
577 xs = x4_ps(p,s)
578 if p < 16.529 :
579 v4v = v2_pT(p, T4_p(p))
580 v4L = v1_pT(p, T4_p(p))
581 else:
582 v4v = v3_ph(p, h4V_p(p))
583 v4L = v3_ph(p, h4L_p(p))
584 return fromSIunit_v((xs * v4v + (1 - xs) * v4L))
585 elif Region == 5:
586 Ts = T5_ps(p, s)
587 return fromSIunit_v(v5_pT(p, Ts))
588 else:
589 return float('NaN')
590
591
592 #***************************************************************************
593 #1.7 Specific entropy (s)
594 #***************************************************************************
595
596 elif fun == 'sv_p':
597 p = toSIunit_p(In1)
598 if p > 0.000611657 and p < 22.06395 :
599 if p < 16.529 :
600 return fromSIunit_s(s2_pT(p, T4_p(p)))
601 else :
602 return fromSIunit_s(s3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)))
603 else :
604 return float('NaN')
605
606 elif fun == 'sl_p':
607 p = toSIunit_p(In1)
608 if p > 0.000611657 and p < 22.06395 :
609 if p < 16.529:
610 return fromSIunit_s(s1_pT(p, T4_p(p)))
611 else:
612 return fromSIunit_s(s3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)))
613 else:
614 return float('NaN')
615
616 elif fun == 'sv_t':
617 T = toSIunit_T(In1)
618 if T > 273.15 and T < 647.096 :
619 if T <= 623.15 :
620 return fromSIunit_s(s2_pT(p4_T(T), T))
621 else :
622 return fromSIunit_s(s3_rhoT(1 / (v3_ph(p4_T(T), h4V_p(p4_T(T)))), T))
623 else :
624 return float('NaN')
625
626 elif fun == 'sl_t':
627 T = toSIunit_T(In1)
628 if T > 273.15 and T < 647.096 :
629 if T <= 623.15 :
630 return fromSIunit_s(s1_pT(p4_T(T), T))
631 else:
632 return fromSIunit_s(s3_rhoT(1 / (v3_ph(p4_T(T), h4L_p(p4_T(T)))), T))
633 else:
634 return float('NaN')
635
636 elif fun == 's_pt':
637 p = toSIunit_p(In1)
638 T = toSIunit_T(In2)
639 Region = region_pT(p, T)
640 if Region == 1:
641 return fromSIunit_s(s1_pT(p, T))
642 elif Region == 2:
643 return fromSIunit_s(s2_pT(p, T))
644 elif Region == 3:
645 hs = h3_pT(p, T)
646 rhos = 1 / v3_ph(p, hs)
647 return fromSIunit_s(s3_rhoT(rhos, T))
648 elif Region == 4:
649 return float('NaN')
650 elif Region == 5:
651 return fromSIunit_s(s5_pT(p, T))
652 else:
653 return float('NaN')
654
655 elif fun == 's_ph':
656 p = toSIunit_p(In1)
657 h = toSIunit_h(In2)
658 Region = region_ph(p, h)
659 if Region == 1:
660 T = T1_ph(p, h)
661 return fromSIunit_s(s1_pT(p, T))
662 elif Region == 2:
663 T = T2_ph(p, h)
664 return fromSIunit_s(s2_pT(p, T))
665 elif Region == 3:
666 rhos = 1.0 / v3_ph(p, h)
667 Ts = T3_ph(p, h)
668 return fromSIunit_s(s3_rhoT(rhos, Ts))
669 elif Region == 4:
670 Ts = T4_p(p)
671 xs = x4_ph(p, h)
672 if p < 16.529 :
673 s4v = s2_pT(p, Ts)
674 s4L = s1_pT(p, Ts)
675 else:
676 v4v = v3_ph(p, h4V_p(p))
677 s4v = s3_rhoT(1.0 / v4v, Ts)
678 v4L = v3_ph(p, h4L_p(p))
679 s4L = s3_rhoT(1.0 / v4L, Ts)
680 return fromSIunit_s((xs * s4v + (1 - xs) * s4L))
681 elif Region == 5:
682 T = T5_ph(p, h)
683 return fromSIunit_s(s5_pT(p, T))
684 else:
685 return float('NaN')
686
687
688
689 #***************************************************************************
690 #*1.17 Vapour fraction
691 #***************************************************************************
692
693 elif fun == 'x_ph':
694 p = toSIunit_p(In1)
695 h = toSIunit_h(In2)
696 if p > 0.000611657 and p < 22.06395 :
697 return fromSIunit_x(x4_ph(p, h))
698 else:
699 return float('NaN')
700
701 elif fun == 'x_ps':
702 p = toSIunit_p(In1)
703 s = toSIunit_s(In2)
704 if p > 0.000611657 and p < 22.06395 :
705 return fromSIunit_x(x4_ps(p, s))
706 else:
707 return float('NaN')
708
709 else:
710 print 'Unknown calling function to XSteam, ',fun, ' See help XSteam for valid calling functions'
711
712
713 '''
714 %***********************************************************************************************************
715 %*2 IAPWS IF 97 Calling functions *
716 %***********************************************************************************************************
717 '''
718
719 #***********************************************************************************************************
720 #*2.1 Functions for region 1
721 #***********************************************************************************************************
722 def v1_pT(p,T):
723 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
724 #%5 Equations for Region 1, Section. 5.1 Basic Equation
725 #%Eqution 7, Table 3, Page 6
726 I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32]
727 J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41]
728 n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26]
729 R = 0.461526 #kJ/(kg K)
730 Pi = p / 16.53
731 tau = 1386.0 / float(T)
732 gamma_der_pi = 0.0
733 for i in range(0,34):
734 gamma_der_pi = gamma_der_pi - n1[i] * I1[i] * (7.1 - Pi) ** (I1[i] - 1) * (tau - 1.222) ** J1[i]
735 return R * T / p * Pi * gamma_der_pi / 1000.0
736
737
738 def h1_pT(p, T):
739 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
740 #%5 Equations for Region 1, Section. 5.1 Basic Equation
741 #%Eqution 7, Table 3, Page 6
742 I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32]
743 J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41]
744 n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26]
745 R = 0.461526 #%kJ/(kg K)
746 Pi = p / 16.53
747 tau = 1386.0 / float(T)
748 gamma_der_tau = 0.0
749 for i in range(0,34):
750 gamma_der_tau = gamma_der_tau + (n1[i] * (7.1 - Pi) ** I1[i] * J1[i] * (tau - 1.222) ** (J1[i] - 1))
751 return R * T * tau * gamma_der_tau
752
753
754 def u1_pT(p, T):
755 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
756 #%5 Equations for Region 1, Section. 5.1 Basic Equation
757 #%Eqution 7, Table 3, Page 6
758 I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32]
759 J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41]
760 n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26]
761 R = 0.461526 #%kJ/(kg K)
762 Pi = p / 16.53
763 tau = 1386.0 / float(T)
764 gamma_der_tau = 0.0
765 gamma_der_pi = 0.0
766 for i in range(0,34):
767 gamma_der_pi = gamma_der_pi - n1[i] * I1[i] * (7.1 - Pi) ** (I1[i] - 1) * (tau - 1.222) ** J1[i]
768 gamma_der_tau = gamma_der_tau + (n1[i] * (7.1 - Pi) ** I1[i] * J1[i] * (tau - 1.222) ** (J1[i] - 1))
769 return R * T * (tau * gamma_der_tau - Pi * gamma_der_pi)
770
771
772 def s1_pT(p, T):
773 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
774 #%5 Equations for Region 1, Section. 5.1 Basic Equation
775 #%Eqution 7, Table 3, Page 6
776 I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32]
777 J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41]
778 n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26]
779 R = 0.461526 # %kJ/(kg K)
780 Pi = p / 16.53
781 tau = 1386.0 / float(T)
782 gamma = 0.0
783 gamma_der_tau = 0.0
784 for i in range(0,34):
785 gamma_der_tau = float(gamma_der_tau) + (float(n1[i]) * (7.1 - float(Pi)) ** float(I1[i]) * float(J1[i]) * (float(tau) - 1.222) ** (float(J1[i]) - 1))
786 gamma = float(gamma) + float(n1[i]) * (7.1 - float(Pi)) ** float(I1[i]) * (float(tau) - 1.222) ** float(J1[i])
787 return R * tau * gamma_der_tau - R * gamma
788
789
790 def Cp1_pT(p, T):
791 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
792 #%5 Equations for Region 1, Section. 5.1 Basic Equation
793 #%Eqution 7, Table 3, Page 6
794 I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32]
795 J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41]
796 n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26]
797 R = 0.461526 # %kJ/(kg K)
798 Pi = p / 16.53
799 tau = 1386.0 / float(T)
800 gamma_der_tautau = 0.0
801 for i in range(0,34):
802 gamma_der_tautau = gamma_der_tautau + (n1[i] * (7.1 - Pi) ** I1[i] * J1[i] * (J1[i] - 1) * (tau - 1.222) ** (J1[i] - 2))
803 return -R * tau ** 2 * gamma_der_tautau
804
805
806 def Cv1_pT(p, T):
807 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
808 #%5 Equations for Region 1, Section. 5.1 Basic Equation
809 #%Eqution 7, Table 3, Page 6
810 I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32]
811 J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41]
812 n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26]
813 R = 0.461526 # %kJ/(kg K)
814 Pi = p / 16.53
815 tau = 1386.0 / float(T)
816 gamma_der_pi = 0.0
817 gamma_der_pipi = 0.0
818 gamma_der_pitau = 0.0
819 gamma_der_tautau = 0.0
820 for i in range(0,34):
821 gamma_der_pi = gamma_der_pi - n1[i] * I1[i] * (7.1 - Pi) ** (I1[i] - 1) * (tau - 1.222) ** J1[i]
822 gamma_der_pipi = gamma_der_pipi + n1[i] * I1[i] * (I1[i] - 1) * (7.1 - Pi) ** (I1[i] - 2) * (tau - 1.222) ** J1[i]
823 gamma_der_pitau = gamma_der_pitau - n1[i] * I1[i] * (7.1 - Pi) ** (I1[i] - 1) * J1[i] * (tau - 1.222) ** (J1[i] - 1)
824 gamma_der_tautau = gamma_der_tautau + n1[i] * (7.1 - Pi) ** I1[i] * J1[i] * (J1[i] - 1) * (tau - 1.222) ** (J1[i] - 2)
825 return R * (-tau ** 2 * gamma_der_tautau + (gamma_der_pi - tau * gamma_der_pitau) ** 2 / gamma_der_pipi)
826
827
828 def w1_pT(p, T):
829 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
830 #%5 Equations for Region 1, Section. 5.1 Basic Equation
831 #%Eqution 7, Table 3, Page 6
832 I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32]
833 J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41]
834 n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26]
835 R = 0.461526 # %kJ/(kg K)
836 Pi = p / 16.53
837 tau = 1386.0 / float(T)
838 gamma_der_pi = 0.0
839 gamma_der_pipi = 0.0
840 gamma_der_pitau = 0.0
841 gamma_der_tautau = 0.0
842 for i in range(0,34):
843 gamma_der_pi = gamma_der_pi - n1[i] * I1[i] * (7.1 - Pi) ** (I1[i] - 1) * (tau - 1.222) ** J1[i]
844 gamma_der_pipi = gamma_der_pipi + n1[i] * I1[i] * (I1[i] - 1) * (7.1 - Pi) ** (I1[i] - 2) * (tau - 1.222) ** J1[i]
845 gamma_der_pitau = gamma_der_pitau - n1[i] * I1[i] * (7.1 - Pi) ** (I1[i] - 1) * J1[i] * (tau - 1.222) ** (J1[i] - 1)
846 gamma_der_tautau = gamma_der_tautau + n1[i] * (7.1 - Pi) ** I1[i] * J1[i] * (J1[i] - 1) * (tau - 1.222) ** (J1[i] - 2)
847 return (1000.0 * R * T * gamma_der_pi ** 2 / ((gamma_der_pi - tau * gamma_der_pitau) ** 2 / (tau ** 2 * gamma_der_tautau) - gamma_der_pipi)) ** 0.5
848
849
850 def T1_ph(p, h):
851 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
852 #%5 Equations for Region 1, Section. 5.1 Basic Equation, 5.2.1 The Backward Equation T ( p,h )
853 #%Eqution 11, Table 6, Page 10
854 I1 = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 4, 5, 6]
855 J1 = [0, 1, 2, 6, 22, 32, 0, 1, 2, 3, 4, 10, 32, 10, 32, 10, 32, 32, 32, 32]
856 n1 = [-238.72489924521, 404.21188637945, 113.49746881718, -5.8457616048039, -1.528548241314E-04, -1.0866707695377E-06, -13.391744872602, 43.211039183559, -54.010067170506, 30.535892203916, -6.5964749423638, 9.3965400878363E-03, 1.157364750534E-07, -2.5858641282073E-05, -4.0644363084799E-09, 6.6456186191635E-08, 8.0670734103027E-11, -9.3477771213947E-13, 5.8265442020601E-15, -1.5020185953503E-17]
857 Pi = p / 1.0
858 eta = h / 2500.0
859 T = 0.0
860 for i in range(0,20):
861 T = T + n1[i] * Pi ** I1[i] * (eta + 1) ** J1[i]
862 return T
863
864
865 def T1_ps(p, s):
866 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
867 #%5 Equations for Region 1, Section. 5.1 Basic Equation, 5.2.2 The Backward Equation T ( p, s )
868 #%Eqution 13, Table 8, Page 11
869 I1 = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 4]
870 J1 = [0, 1, 2, 3, 11, 31, 0, 1, 2, 3, 12, 31, 0, 1, 2, 9, 31, 10, 32, 32]
871 n1 = [174.78268058307, 34.806930892873, 6.5292584978455, 0.33039981775489, -1.9281382923196E-07, -2.4909197244573E-23, -0.26107636489332, 0.22592965981586, -0.064256463395226, 7.8876289270526E-03, 3.5672110607366E-10, 1.7332496994895E-24, 5.6608900654837E-04, -3.2635483139717E-04, 4.4778286690632E-05, -5.1322156908507E-10, -4.2522657042207E-26, 2.6400441360689E-13, 7.8124600459723E-29, -3.0732199903668E-31]
872 Pi = p / 1.0
873 Sigma = s / 1.0
874 T = 0.0
875 for i in range(0,20):
876 T = T + n1[i] * Pi ** I1[i] * (Sigma + 2) ** J1[i]
877 return T
878
879
880 def p1_hs(h, s):
881 #%Supplementary Release on Backward Equations for Pressure as a Function of Enthalpy and Entropy p(h,s) to the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
882 #%5 Backward Equation p(h,s) for Region 1
883 #%Eqution 1, Table 2, Page 5
884 I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 5]
885 J1 = [0, 1, 2, 4, 5, 6, 8, 14, 0, 1, 4, 6, 0, 1, 10, 4, 1, 4, 0]
886 n1 = [-0.691997014660582, -18.361254878756, -9.28332409297335, 65.9639569909906, -16.2060388912024, 450.620017338667, 854.68067822417, 6075.23214001162, 32.6487682621856, -26.9408844582931, -319.9478483343, -928.35430704332, 30.3634537455249, -65.0540422444146, -4309.9131651613, -747.512324096068, 730.000345529245, 1142.84032569021, -436.407041874559]
887 eta = h / 3400.0
888 Sigma = s / 7.6
889 p = 0.0
890 for i in range(0,19):
891 p = p + n1[i] * (eta + 0.05) ** I1[i] * (Sigma + 0.05) ** J1[i]
892 return p * 100.0
893
894 def T1_prho(p ,rho):
895 # %Solve by iteration. Observe that for low temperatures this equation has 2 solutions.
896 # %Solve with half interval method
897 '''
898 Low_Bound = 273.15;
899 High_Bound = T4_p(p);
900 rhos=-1000;
901 while abs(rho - rhos) > 0.00001
902 Ts = (Low_Bound + High_Bound) / 2;
903 rhos = 1 / v1_pT(p, Ts);
904 if rhos < rho
905 High_Bound = Ts;
906 else
907 Low_Bound = Ts;
908 end
909 end
910 T1_prho = Ts;
911 '''
912 print 'T1_prho Not yet translated, because of iteration'
913
914
915
916 #***********************************************************************************************************
917 #*2.2 Functions for region 2
918 #***********************************************************************************************************
919 def v2_pT(p, T):
920 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
921 #%6 Equations for Region 2, Section. 6.1 Basic Equation
922 #%Table 11 and 12, Page 14 and 15
923 J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3]
924 n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307]
925 Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24]
926 Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58]
927 nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07]
928 R = 0.461526 # %kJ/(kg K)
929 Pi = p
930 tau = 540.0 / float(T)
931 g0_pi = 1.0 / float(Pi)
932 gr_pi = 0.0
933 for i in range(0,43):
934 gr_pi = gr_pi + nr[i] * Ir[i] * Pi ** (Ir[i] - 1) * (tau - 0.5) ** Jr[i]
935 return R * T / p * Pi * (g0_pi + gr_pi) / 1000.0
936
937
938 def h2_pT(p, T):
939 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
940 #%6 Equations for Region 2, Section. 6.1 Basic Equation
941 #%Table 11 and 12, Page 14 and 15
942 J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3]
943 n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307]
944 Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24]
945 Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58]
946 nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07]
947 R = 0.461526 # %kJ/(kg K)
948 Pi = p
949 tau = 540.0 / float(T)
950 g0_tau = 0.0
951 for i in range(0,9):
952 g0_tau = g0_tau + n0[i] * J0[i] * tau ** (J0[i] - 1)
953 gr_tau = 0.0
954 for i in range(0,43):
955 gr_tau = gr_tau + nr[i] * Pi ** Ir[i] * Jr[i] * (tau - 0.5) ** (Jr[i] - 1)
956 return R * T * tau * (g0_tau + gr_tau)
957
958
959 def u2_pT(p, T):
960 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
961 #%6 Equations for Region 2, Section. 6.1 Basic Equation
962 #%Table 11 and 12, Page 14 and 15
963 J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3]
964 n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307]
965 Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24]
966 Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58]
967 nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07]
968 R = 0.461526 # %kJ/(kg K)
969 Pi = p
970 tau = 540 / float(T)
971 g0_pi = 1 / float(Pi)
972 g0_tau = 0.0
973 for i in range(0,9):
974 g0_tau = g0_tau + n0[i] * J0[i] * tau ** (J0[i] - 1)
975 gr_pi = 0.0
976 gr_tau = 0.0
977 for i in range(0,43):
978 gr_pi = gr_pi + nr[i] * Ir[i] * Pi ** (Ir[i] - 1) * (tau - 0.5) ** Jr[i]
979 gr_tau = gr_tau + nr[i] * Pi ** Ir[i] * Jr[i] * (tau - 0.5) ** (Jr[i] - 1)
980 u2_pT = R * T * (tau * (g0_tau + gr_tau) - Pi * (g0_pi + gr_pi))
981
982
983 def s2_pT(p, T):
984 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
985 #%6 Equations for Region 2, Section. 6.1 Basic Equation
986 #%Table 11 and 12, Page 14 and 15
987 J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3]
988 n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307]
989 Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24]
990 Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58]
991 nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07]
992 R = 0.461526 # %kJ/(kg K)
993 Pi = float(p)
994 tau = 540 / float(T)
995 g0 = log(Pi)
996 g0_tau = 0.0
997 for i in range(0,9):
998 g0 = g0 + n0[i] * tau ** J0[i]
999 g0_tau = g0_tau + n0[i] * J0[i] * tau ** (J0[i] - 1)
1000
1001 gr = 0.0
1002 gr_tau = 0.0
1003 for i in range(0,43):
1004 gr = gr + nr[i] * Pi ** Ir[i] * (tau - 0.5) ** Jr[i]
1005 gr_tau = gr_tau + nr[i] * Pi ** Ir[i] * Jr[i] * (tau - 0.5) ** (Jr[i] - 1)
1006 return R * (tau * (g0_tau + gr_tau) - (g0 + gr))
1007
1008
1009
1010 def Cp2_pT(p, T):
1011 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1012 #%6 Equations for Region 2, Section. 6.1 Basic Equation
1013 #%Table 11 and 12, Page 14 and 15
1014 J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3]
1015 n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307]
1016 Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24]
1017 Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58]
1018 nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07]
1019 R = 0.461526 # %kJ/(kg K)
1020 Pi = float(p)
1021 tau = 540 / float(T)
1022 g0_tautau = 0.0
1023 for i in range(01,9):
1024 g0_tautau = g0_tautau + n0[i] * J0[i] * (J0[i] - 1) * tau ** (J0[i] - 2)
1025 gr_tautau = 0.0
1026 for i in range(0,43):
1027 gr_tautau = gr_tautau + nr[i] * Pi ** Ir[i] * Jr[i] * (Jr[i] - 1) * (tau - 0.5) ** (Jr[i] - 2)
1028 Cp2_pT = -R * tau ** 2 * (g0_tautau + gr_tautau)
1029
1030
1031
1032 def Cv2_pT(p, T):
1033 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1034 #%6 Equations for Region 2, Section. 6.1 Basic Equation
1035 #%Table 11 and 12, Page 14 and 15
1036 J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3]
1037 n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307]
1038 Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24]
1039 Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58]
1040 nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07]
1041 R = 0.461526 #%kJ/(kg K)
1042 Pi = float(p)
1043 tau = 540 / float(T)
1044 g0_tautau = 0.0
1045 for i in range(0,9):
1046 g0_tautau = g0_tautau + n0[i] * J0[i] * (J0[i] - 1) * tau ** (J0[i] - 2)
1047 gr_pi = 0.0
1048 gr_pitau = 0.0
1049 gr_pipi = 0.0
1050 gr_tautau = 0.0
1051 for i in range(0,43):
1052 gr_pi = gr_pi + nr[i] * Ir[i] * Pi ** (Ir[i] - 1) * (tau - 0.5) ** Jr[i]
1053 gr_pipi = gr_pipi + nr[i] * Ir[i] * (Ir[i] - 1) * Pi ** (Ir[i] - 2) * (tau - 0.5) ** Jr[i]
1054 gr_pitau = gr_pitau + nr[i] * Ir[i] * Pi ** (Ir[i] - 1) * Jr[i] * (tau - 0.5) ** (Jr[i] - 1)
1055 gr_tautau = gr_tautau + nr[i] * Pi ** Ir[i] * Jr[i] * (Jr[i] - 1) * (tau - 0.5) ** (Jr[i] - 2)
1056 return R * (-tau ** 2 * (g0_tautau + gr_tautau) - (1 + Pi * gr_pi - tau * Pi * gr_pitau) ** 2 / (1 - Pi ** 2 * gr_pipi))
1057
1058
1059 def w2_pT(p, T):
1060 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1061 #%6 Equations for Region 2, Section. 6.1 Basic Equation
1062 #%Table 11 and 12, Page 14 and 15
1063 J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3]
1064 n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307]
1065 Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24]
1066 Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58]
1067 nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07]
1068 R = 0.461526 #%kJ/(kg K)
1069 Pi = float(p)
1070 tau = 540 / float(T)
1071 g0_tautau = 0.0
1072 for i in range(0,9):
1073 g0_tautau = g0_tautau + n0[i] * J0[i] * (J0[i] - 1) * tau ** (J0[i] - 2)
1074 gr_pi = 0.0
1075 gr_pitau = 0.0
1076 gr_pipi = 0.0
1077 gr_tautau = 0.0
1078 for i in range(0,43):
1079 gr_pi = gr_pi + nr[i] * Ir[i] * Pi ** (Ir[i] - 1) * (tau - 0.5) ** Jr[i]
1080 gr_pipi = gr_pipi + nr[i] * Ir[i] * (Ir[i] - 1) * Pi ** (Ir[i] - 2) * (tau - 0.5) ** Jr[i]
1081 gr_pitau = gr_pitau + nr[i] * Ir[i] * Pi ** (Ir[i] - 1) * Jr[i] * (tau - 0.5) ** (Jr[i] - 1)
1082 gr_tautau = gr_tautau + nr[i] * Pi ** Ir[i] * Jr[i] * (Jr[i] - 1) * (tau - 0.5) ** (Jr[i] - 2)
1083 return (1000.0 * R * T * (1 + 2 * Pi * gr_pi + Pi ** 2 * gr_pi ** 2) / ((1 - Pi ** 2 * gr_pipi) + (1 + Pi * gr_pi - tau * Pi * gr_pitau) ** 2 / (tau ** 2 * (g0_tautau + gr_tautau)))) ** 0.5
1084
1085
1086 def T2_ph(p, h):
1087 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1088 #%6 Equations for Region 2,6.3.1 The Backward Equations T( p, h ) for Subregions 2a, 2b, and 2c
1089 if p < 4 :
1090 sub_reg = 1
1091 else:
1092 if p < (905.84278514723 - 0.67955786399241 * h + 1.2809002730136E-04 * h ** 2) :
1093 sub_reg = 2
1094 else:
1095 sub_reg = 3
1096 if sub_reg == 1:
1097 #%Subregion A
1098 #%Table 20, Eq 22, page 22
1099 Ji = [0, 1, 2, 3, 7, 20, 0, 1, 2, 3, 7, 9, 11, 18, 44, 0, 2, 7, 36, 38, 40, 42, 44, 24, 44, 12, 32, 44, 32, 36, 42, 34, 44, 28]
1100 Ii = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 7]
1101 ni = [1089.8952318288, 849.51654495535, -107.81748091826, 33.153654801263, -7.4232016790248, 11.765048724356, 1.844574935579, -4.1792700549624, 6.2478196935812, -17.344563108114, -200.58176862096, 271.96065473796, -455.11318285818, 3091.9688604755, 252266.40357872, -6.1707422868339E-03, -0.31078046629583, 11.670873077107, 128127984.04046, -985549096.23276, 2822454697.3002, -3594897141.0703, 1722734991.3197, -13551.334240775, 12848734.66465, 1.3865724283226, 235988.32556514, -13105236.545054, 7399.9835474766, -551966.9703006, 3715408.5996233, 19127.72923966, -415351.64835634, -62.459855192507]
1102 Ts = 0.0
1103 hs = h / 2000.0
1104 for i in range(0,34):
1105 Ts = Ts + ni[i] * p ** (Ii[i]) * (hs - 2.1) ** Ji[i]
1106 return Ts
1107 elif sub_reg == 2:
1108 #%Subregion B
1109 #%Table 21, Eq 23, page 23
1110 Ji = [0, 1, 2, 12, 18, 24, 28, 40, 0, 2, 6, 12, 18, 24, 28, 40, 2, 8, 18, 40, 1, 2, 12, 24, 2, 12, 18, 24, 28, 40, 18, 24, 40, 28, 2, 28, 1, 40]
1111 Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 7, 7, 9, 9]
1112 ni = [1489.5041079516, 743.07798314034, -97.708318797837, 2.4742464705674, -0.63281320016026, 1.1385952129658, -0.47811863648625, 8.5208123431544E-03, 0.93747147377932, 3.3593118604916, 3.3809355601454, 0.16844539671904, 0.73875745236695, -0.47128737436186, 0.15020273139707, -0.002176411421975, -0.021810755324761, -0.10829784403677, -0.046333324635812, 7.1280351959551E-05, 1.1032831789999E-04, 1.8955248387902E-04, 3.0891541160537E-03, 1.3555504554949E-03, 2.8640237477456E-07, -1.0779857357512E-05, -7.6462712454814E-05, 1.4052392818316E-05, -3.1083814331434E-05, -1.0302738212103E-06, 2.821728163504E-07, 1.2704902271945E-06, 7.3803353468292E-08, -1.1030139238909E-08, -8.1456365207833E-14, -2.5180545682962E-11, -1.7565233969407E-18, 8.6934156344163E-15]
1113 Ts = 0.0
1114 hs = h / 2000.0
1115 for i in range(0,38):
1116 Ts = Ts + ni[i] * (p - 2) ** (Ii[i]) * (hs - 2.6) ** Ji[i]
1117 return Ts
1118 else:
1119 #%Subregion C
1120 #%Table 22, Eq 24, page 24
1121 Ji = [0, 4, 0, 2, 0, 2, 0, 1, 0, 2, 0, 1, 4, 8, 4, 0, 1, 4, 10, 12, 16, 20, 22]
1122 Ii = [-7, -7, -6, -6, -5, -5, -2, -2, -1, -1, 0, 0, 1, 1, 2, 6, 6, 6, 6, 6, 6, 6, 6]
1123 ni = [-3236839855524.2, 7326335090218.1, 358250899454.47, -583401318515.9, -10783068217.47, 20825544563.171, 610747.83564516, 859777.2253558, -25745.72360417, 31081.088422714, 1208.2315865936, 482.19755109255, 3.7966001272486, -10.842984880077, -0.04536417267666, 1.4559115658698E-13, 1.126159740723E-12, -1.7804982240686E-11, 1.2324579690832E-07, -1.1606921130984E-06, 2.7846367088554E-05, -5.9270038474176E-04, 1.2918582991878E-03]
1124 Ts = 0.0
1125 hs = h / 2000.0
1126 for i in range(0,23):
1127 Ts = Ts + ni[i] * (p + 25) ** (Ii[i]) * (hs - 1.8) ** Ji[i]
1128 return Ts
1129
1130
1131 def T2_ps(p, s):
1132 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1133 #%6 Equations for Region 2,6.3.2 The Backward Equations T( p, s ) for Subregions 2a, 2b, and 2c
1134 #%Page 26
1135 if p < 4 :
1136 sub_reg = 1
1137 else :
1138 if s < 5.85 :
1139 sub_reg = 3
1140 else:
1141 sub_reg = 2
1142
1143 if sub_reg == 1:
1144 #%Subregion A
1145 #%Table 25, Eq 25, page 26
1146 Ii = [-1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.25, -1.25, -1.25, -1, -1, -1, -1, -1, -1, -0.75, -0.75, -0.5, -0.5, -0.5, -0.5, -0.25, -0.25, -0.25, -0.25, 0.25, 0.25, 0.25, 0.25, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.75, 0.75, 0.75, 0.75, 1, 1, 1.25, 1.25, 1.5, 1.5]
1147 Ji = [-24, -23, -19, -13, -11, -10, -19, -15, -6, -26, -21, -17, -16, -9, -8, -15, -14, -26, -13, -9, -7, -27, -25, -11, -6, 1, 4, 8, 11, 0, 1, 5, 6, 10, 14, 16, 0, 4, 9, 17, 7, 18, 3, 15, 5, 18]
1148 ni = [-392359.83861984, 515265.7382727, 40482.443161048, -321.93790923902, 96.961424218694, -22.867846371773, -449429.14124357, -5011.8336020166, 0.35684463560015, 44235.33584819, -13673.388811708, 421632.60207864, 22516.925837475, 474.42144865646, -149.31130797647, -197811.26320452, -23554.39947076, -19070.616302076, 55375.669883164, 3829.3691437363, -603.91860580567, 1936.3102620331, 4266.064369861, -5978.0638872718, -704.01463926862, 338.36784107553, 20.862786635187, 0.033834172656196, -4.3124428414893E-05, 166.53791356412, -139.86292055898, -0.78849547999872, 0.072132411753872, -5.9754839398283E-03, -1.2141358953904E-05, 2.3227096733871E-07, -10.538463566194, 2.0718925496502, -0.072193155260427, 2.074988708112E-07, -0.018340657911379, 2.9036272348696E-07, 0.21037527893619, 2.5681239729999E-04, -0.012799002933781, -8.2198102652018E-06]
1149 Pi = float(p)
1150 Sigma = s / 2.0
1151 teta = 0.0
1152 for i in range(0,46):
1153 teta = teta + ni[i] * Pi ** Ii[i] * (Sigma - 2) ** Ji[i]
1154 return teta
1155 elif sub_reg == 2:
1156 #%Subregion B
1157 #%Table 26, Eq 26, page 27
1158 Ii = [-6, -6, -5, -5, -4, -4, -4, -3, -3, -3, -3, -2, -2, -2, -2, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 5]
1159 Ji = [0, 11, 0, 11, 0, 1, 11, 0, 1, 11, 12, 0, 1, 6, 10, 0, 1, 5, 8, 9, 0, 1, 2, 4, 5, 6, 9, 0, 1, 2, 3, 7, 8, 0, 1, 5, 0, 1, 3, 0, 1, 0, 1, 2]
1160 ni = [316876.65083497, 20.864175881858, -398593.99803599, -21.816058518877, 223697.85194242, -2784.1703445817, 9.920743607148, -75197.512299157, 2970.8605951158, -3.4406878548526, 0.38815564249115, 17511.29508575, -1423.7112854449, 1.0943803364167, 0.89971619308495, -3375.9740098958, 471.62885818355, -1.9188241993679, 0.41078580492196, -0.33465378172097, 1387.0034777505, -406.63326195838, 41.72734715961, 2.1932549434532, -1.0320050009077, 0.35882943516703, 5.2511453726066E-03, 12.838916450705, -2.8642437219381, 0.56912683664855, -0.099962954584931, -3.2632037778459E-03, 2.3320922576723E-04, -0.1533480985745, 0.029072288239902, 3.7534702741167E-04, 1.7296691702411E-03, -3.8556050844504E-04, -3.5017712292608E-05, -1.4566393631492E-05, 5.6420857267269E-06, 4.1286150074605E-08, -2.0684671118824E-08, 1.6409393674725E-09]
1161 Pi = float(p)
1162 Sigma = s / 0.7853
1163 teta = 0.0
1164 for i in range(0,44):
1165 teta = teta + ni[i] * Pi ** Ii[i] * (10 - Sigma) ** Ji[i]
1166 return teta
1167 else:
1168 #%Subregion C
1169 #%Table 27, Eq 27, page 28
1170 Ii = [-2, -2, -1, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 7, 7, 7, 7, 7]
1171 Ji = [0, 1, 0, 0, 1, 2, 3, 0, 1, 3, 4, 0, 1, 2, 0, 1, 5, 0, 1, 4, 0, 1, 2, 0, 1, 0, 1, 3, 4, 5]
1172 ni = [909.68501005365, 2404.566708842, -591.6232638713, 541.45404128074, -270.98308411192, 979.76525097926, -469.66772959435, 14.399274604723, -19.104204230429, 5.3299167111971, -21.252975375934, -0.3114733441376, 0.60334840894623, -0.042764839702509, 5.8185597255259E-03, -0.014597008284753, 5.6631175631027E-03, -7.6155864584577E-05, 2.2440342919332E-04, -1.2561095013413E-05, 6.3323132660934E-07, -2.0541989675375E-06, 3.6405370390082E-08, -2.9759897789215E-09, 1.0136618529763E-08, 5.9925719692351E-12, -2.0677870105164E-11, -2.0874278181886E-11, 1.0162166825089E-10, -1.6429828281347E-10]
1173 Pi = float(p)
1174 Sigma = s / 2.9251
1175 teta = 0.0
1176 for i in range(0,30):
1177 teta = teta + ni[i] * Pi ** Ii[i] * (2 - Sigma) ** Ji[i]
1178 return teta
1179
1180
1181 def p2_hs(h, s):
1182 #%Supplementary Release on Backward Equations for Pressure as a function of Enthalpy and Entropy p(h,s) to the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
1183 #%Chapter 6:Backward Equations p(h,s) for Region 2
1184 if h < (-3498.98083432139 + 2575.60716905876 * s - 421.073558227969 * s ** 2 + 27.6349063799944 * s ** 3):
1185 sub_reg = 1
1186 else:
1187 if s < 5.85 :
1188 sub_reg = 3
1189 else:
1190 sub_reg = 2
1191
1192 if sub_reg == 1:
1193 #%Subregion A
1194 #%Table 6, Eq 3, page 8
1195 Ii = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 4, 5, 5, 6, 7]
1196 Ji = [1, 3, 6, 16, 20, 22, 0, 1, 2, 3, 5, 6, 10, 16, 20, 22, 3, 16, 20, 0, 2, 3, 6, 16, 16, 3, 16, 3, 1]
1197 ni = [-1.82575361923032E-02, -0.125229548799536, 0.592290437320145, 6.04769706185122, 238.624965444474, -298.639090222922, 0.051225081304075, -0.437266515606486, 0.413336902999504, -5.16468254574773, -5.57014838445711, 12.8555037824478, 11.414410895329, -119.504225652714, -2847.7798596156, 4317.57846408006, 1.1289404080265, 1974.09186206319, 1516.12444706087, 1.41324451421235E-02, 0.585501282219601, -2.97258075863012, 5.94567314847319, -6236.56565798905, 9659.86235133332, 6.81500934948134, -6332.07286824489, -5.5891922446576, 4.00645798472063E-02]
1198 eta = h / 4200.0
1199 Sigma = s / 12.0
1200 Pi = 0.0
1201 for i in range(0,29):
1202 Pi = Pi + ni[i] * (eta - 0.5) ** Ii[i] * (Sigma - 1.2) ** Ji[i]
1203 return Pi ** 4.0 * 4.0
1204 elif sub_reg == 2:
1205 #%Subregion B
1206 #%Table 7, Eq 4, page 9
1207 Ii = [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 8, 8, 12, 14]
1208 Ji = [0, 1, 2, 4, 8, 0, 1, 2, 3, 5, 12, 1, 6, 18, 0, 1, 7, 12, 1, 16, 1, 12, 1, 8, 18, 1, 16, 1, 3, 14, 18, 10, 16]
1209 ni = [8.01496989929495E-02, -0.543862807146111, 0.337455597421283, 8.9055545115745, 313.840736431485, 0.797367065977789, -1.2161697355624, 8.72803386937477, -16.9769781757602, -186.552827328416, 95115.9274344237, -18.9168510120494, -4334.0703719484, 543212633.012715, 0.144793408386013, 128.024559637516, -67230.9534071268, 33697238.0095287, -586.63419676272, -22140322476.9889, 1716.06668708389, -570817595.806302, -3121.09693178482, -2078413.8463301, 3056059461577.86, 3221.57004314333, 326810259797.295, -1441.04158934487, 410.694867802691, 109077066873.024, -24796465425889.3, 1888019068.65134, -123651009018773]
1210 eta = h / 4100.0
1211 Sigma = s / 7.9
1212 Pi = 0.0
1213 for i in range(0,33):
1214 Pi = Pi + ni[i] * (eta - 0.6) ** Ii[i] * (Sigma - 1.01) ** Ji[i]
1215 return Pi ** 4 * 100.0
1216 else:
1217 #%Subregion C
1218 #%Table 8, Eq 5, page 10
1219 Ii = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 5, 5, 5, 5, 6, 6, 10, 12, 16]
1220 Ji = [0, 1, 2, 3, 4, 8, 0, 2, 5, 8, 14, 2, 3, 7, 10, 18, 0, 5, 8, 16, 18, 18, 1, 4, 6, 14, 8, 18, 7, 7, 10]
1221 ni = [0.112225607199012, -3.39005953606712, -32.0503911730094, -197.5973051049, -407.693861553446, 13294.3775222331, 1.70846839774007, 37.3694198142245, 3581.44365815434, 423014.446424664, -751071025.760063, 52.3446127607898, -228.351290812417, -960652.417056937, -80705929.2526074, 1626980172256.69, 0.772465073604171, 46392.9973837746, -13731788.5134128, 1704703926305.12, -25110462818730.8, 31774883083552, 53.8685623675312, -55308.9094625169, -1028615.22421405, 2042494187562.34, 273918446.626977, -2.63963146312685E+15, -1078908541.08088, -29649262098.0124, -1.11754907323424E+15]
1222 eta = h / 3500.0
1223 Sigma = s / 5.9
1224 Pi = 0.0
1225 for i in range(0,31):
1226 Pi = Pi + ni[i] * (eta - 0.7) ** Ii[i] * (Sigma - 1.1) ** Ji[i]
1227 return Pi ** 4.0 * 100.0
1228
1229
1230 def T2_prho(p,rho):
1231 # %Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
1232 # %Solve with half interval method
1233 '''
1234 if p < 16.5292
1235 Low_Bound = T4_p(p);
1236 else
1237 Low_Bound = B23T_p(p);
1238 end
1239 High_Bound = 1073.15;
1240 rhos=-1000;
1241 while abs(rho - rhos) > 0.000001
1242 Ts = (Low_Bound + High_Bound) / 2;
1243 rhos = 1 / v2_pT(p, Ts);
1244 if rhos < rho
1245 High_Bound = Ts;
1246 else
1247 Low_Bound = Ts;
1248 end
1249 end
1250 T2_prho = Ts;
1251 '''
1252 print 'T2_prho Not yet translated, because of iteration'
1253
1254
1255
1256
1257 #***********************************************************************************************************
1258 #*2.3 Functions for region 3
1259 #***********************************************************************************************************
1260 def p3_rhoT(rho, T):
1261 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1262 #%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
1263 #%Table 30 and 31, Page 30 and 31
1264 Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11]
1265 Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26]
1266 ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05]
1267 R = 0.461526 #%kJ/(KgK)
1268 tc = 647.096 #%K
1269 pc = 22.064 #%MPa
1270 rhoc = 322.0 #%kg/m3
1271 delta = rho / float(rhoc)
1272 tau = tc / float(T)
1273 fidelta = 0.0
1274 for i in range(1,40):
1275 fidelta = fidelta + ni[i] * Ii[i] * delta ** (Ii[i] - 1) * tau ** Ji[i]
1276 fidelta = fidelta + ni[0] / delta
1277 return rho * R * T * delta * fidelta / 1000.0
1278
1279
1280 def u3_rhoT(rho, T):
1281 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1282 #%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
1283 #%Table 30 and 31, Page 30 and 31
1284 Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11]
1285 Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26]
1286 ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05]
1287 R = 0.461526 #%kJ/(KgK)
1288 tc = 647.096 #%K
1289 pc = 22.064 #%MPa
1290 rhoc = 322 #%kg/m3
1291 delta = rho / float(rhoc)
1292 tau = tc / float(T)
1293 fitau = 0.0
1294 for i in range(1,40):
1295 fitau = fitau + ni[i] * delta ** Ii[i] * Ji[i] * tau ** (Ji[i] - 1)
1296 return R * T * (tau * fitau)
1297
1298
1299 def h3_rhoT(rho, T):
1300 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1301 #%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
1302 #%Table 30 and 31, Page 30 and 31
1303 Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11]
1304 Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26]
1305 ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05]
1306 R = 0.461526 # %kJ/(KgK)
1307 tc = 647.096 # %K
1308 pc = 22.064 # %MPa
1309 rhoc = 322 # %kg/m3
1310 delta = rho / float(rhoc)
1311 tau = tc / float(T)
1312 fidelta = 0.0
1313 fitau = 0.0
1314 for i in range(1,40):
1315 fidelta = fidelta + ni[i] * Ii[i] * delta ** (Ii[i] - 1) * tau ** Ji[i]
1316 fitau = fitau + ni[i] * delta ** Ii[i] * Ji[i] * tau ** (Ji[i] - 1)
1317 fidelta = fidelta + ni[0] / delta
1318 return R * T * (tau * fitau + delta * fidelta)
1319
1320
1321 def s3_rhoT(rho, T):
1322 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1323 #%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
1324 #%Table 30 and 31, Page 30 and 31
1325 Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11]
1326 Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26]
1327 ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05]
1328 R = 0.461526 #%kJ/(KgK)
1329 tc = 647.096 #%K
1330 pc = 22.064 #%MPa
1331 rhoc = 322 #%kg/m3
1332 delta = rho / float(rhoc)
1333 tau = tc / float(T)
1334 fi = 0.0
1335 fitau = 0.0
1336 for i in range(1,40):
1337 fi = fi + ni[i] * delta ** Ii[i] * tau ** Ji[i]
1338 fitau = fitau + ni[i] * delta ** Ii[i] * Ji[i] * tau ** (Ji[i] - 1)
1339 fi = fi + ni[0] * log(delta)
1340 return R * (tau * fitau - fi)
1341
1342
1343 def Cp3_rhoT(rho, T):
1344 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1345 #%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
1346 #%Table 30 and 31, Page 30 and 31
1347 Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11]
1348 Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26]
1349 ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05]
1350 R = 0.461526 #%kJ/(KgK)
1351 tc = 647.096 #%K
1352 pc = 22.064 #%MPa
1353 rhoc = 322.0 #%kg/m3
1354 delta = rho / float(rhoc)
1355 tau = tc / float(T)
1356 fitautau = 0.0
1357 fidelta = 0.0
1358 fideltatau = 0.0
1359 fideltadelta = 0.0
1360 for i in range(1,40):
1361 fitautau = fitautau + ni[i] * delta ** Ii[i] * Ji[i] * (Ji[i] - 1) * tau ** (Ji[i] - 2)
1362 fidelta = fidelta + ni[i] * Ii[i] * delta ** (Ii[i] - 1) * tau ** Ji[i]
1363 fideltatau = fideltatau + ni[i] * Ii[i] * delta ** (Ii[i] - 1) * Ji[i] * tau ** (Ji[i] - 1)
1364 fideltadelta = fideltadelta + ni[i] * Ii[i] * (Ii[i] - 1) * delta ** (Ii[i] - 2) * tau ** Ji[i]
1365 fidelta = fidelta + ni[0] / delta
1366 fideltadelta = fideltadelta - ni[0] / float((delta ** 2))
1367 return R * (-tau ** 2 * fitautau + (delta * fidelta - delta * tau * fideltatau) ** 2 / (2 * delta * fidelta + delta ** 2 * fideltadelta))
1368
1369
1370 def Cv3_rhoT(rho, T):
1371 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1372 #%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
1373 #%Table 30 and 31, Page 30 and 31
1374 Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11]
1375 Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26]
1376 ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05]
1377 R = 0.461526 #%kJ/(KgK)
1378 tc = 647.096 #%K
1379 pc = 22.064 #%MPa
1380 rhoc = 322 #%kg/m3
1381 delta = rho / float(rhoc)
1382 tau = tc / float(T)
1383 fitautau = 0.0
1384 for i in range(0,40):
1385 fitautau = fitautau + ni[i] * delta ** Ii[i] * Ji[i] * (Ji[i] - 1) * tau ** (Ji[i] - 2)
1386 return R * -(tau * tau * fitautau)
1387
1388
1389 def w3_rhoT(rho, T):
1390 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1391 #%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
1392 #%Table 30 and 31, Page 30 and 31
1393 Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11]
1394 Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26]
1395 ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05]
1396 R = 0.461526 #%kJ/(KgK)
1397 tc = 647.096 #%K
1398 pc = 22.064 #%MPa
1399 rhoc = 322 #%kg/m3
1400 delta = rho / float(rhoc)
1401 tau = tc / float(T)
1402 fitautau = 0.0
1403 fidelta = 0.0
1404 fideltatau = 0.0
1405 fideltadelta = 0.0
1406 for i in range(1,40):
1407 fitautau = fitautau + ni[i] * delta ** Ii[i] * Ji[i] * (Ji[i] - 1) * tau ** (Ji[i] - 2)
1408 fidelta = fidelta + ni[i] * Ii[i] * delta ** (Ii[i] - 1) * tau ** Ji[i]
1409 fideltatau = fideltatau + ni[i] * Ii[i] * delta ** (Ii[i] - 1) * Ji[i] * tau ** (Ji[i] - 1)
1410 fideltadelta = fideltadelta + ni[i] * Ii[i] * (Ii[i] - 1) * delta ** (Ii[i] - 2) * tau ** Ji[i]
1411 fidelta = fidelta + ni[0] / delta
1412 fideltadelta = fideltadelta - ni[0] / (delta ** 2)
1413 return (1000.0 * R * T * (2 * delta * fidelta + delta ** 2 * fideltadelta - (delta * fidelta - delta * tau * fideltatau) ** 2 / (tau ** 2 * fitautau))) ** 0.5
1414
1415
1416 def T3_ph(p, h):
1417 #%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
1418 #%2004
1419 #%Section 3.3 Backward Equations T(p,h) and v(p,h) for Subregions 3a and 3b
1420 #%Boundary equation, Eq 1 Page 5
1421 h3ab = 2014.64004206875 + 3.74696550136983 * p - 2.19921901054187E-02 * p ** 2 + 8.7513168600995E-05 * p ** 3
1422 if h < h3ab :
1423 #%Subregion 3a
1424 #%Eq 2, Table 3, Page 7
1425 Ii = [-12, -12, -12, -12, -12, -12, -12, -12, -10, -10, -10, -8, -8, -8, -8, -5, -3, -2, -2, -2, -1, -1, 0, 0, 1, 3, 3, 4, 4, 10, 12]
1426 Ji = [0, 1, 2, 6, 14, 16, 20, 22, 1, 5, 12, 0, 2, 4, 10, 2, 0, 1, 3, 4, 0, 2, 0, 1, 1, 0, 1, 0, 3, 4, 5]
1427 ni = [-1.33645667811215E-07, 4.55912656802978E-06, -1.46294640700979E-05, 6.3934131297008E-03, 372.783927268847, -7186.54377460447, 573494.7521034, -2675693.29111439, -3.34066283302614E-05, -2.45479214069597E-02, 47.8087847764996, 7.64664131818904E-06, 1.28350627676972E-03, 1.71219081377331E-02, -8.51007304583213, -1.36513461629781E-02, -3.84460997596657E-06, 3.37423807911655E-03, -0.551624873066791, 0.72920227710747, -9.92522757376041E-03, -0.119308831407288, 0.793929190615421, 0.454270731799386, 0.20999859125991, -6.42109823904738E-03, -0.023515586860454, 2.52233108341612E-03, -7.64885133368119E-03, 1.36176427574291E-02, -1.33027883575669E-02]
1428 ps = p / 100.0
1429 hs = h / 2300.0
1430 Ts = 0.0
1431 for i in range(0,31):
1432 Ts = Ts + ni[i] * (ps + 0.24) ** Ii[i] * (hs - 0.615) ** Ji[i]
1433 return Ts * 760.0
1434 else :
1435 #%Subregion 3b
1436 #%Eq 3, Table 4, Page 7,8
1437 Ii = [-12, -12, -10, -10, -10, -10, -10, -8, -8, -8, -8, -8, -6, -6, -6, -4, -4, -3, -2, -2, -1, -1, -1, -1, -1, -1, 0, 0, 1, 3, 5, 6, 8]
1438 Ji = [0, 1, 0, 1, 5, 10, 12, 0, 1, 2, 4, 10, 0, 1, 2, 0, 1, 5, 0, 4, 2, 4, 6, 10, 14, 16, 0, 2, 1, 1, 1, 1, 1]
1439 ni = [3.2325457364492E-05, -1.27575556587181E-04, -4.75851877356068E-04, 1.56183014181602E-03, 0.105724860113781, -85.8514221132534, 724.140095480911, 2.96475810273257E-03, -5.92721983365988E-03, -1.26305422818666E-02, -0.115716196364853, 84.9000969739595, -1.08602260086615E-02, 1.54304475328851E-02, 7.50455441524466E-02, 2.52520973612982E-02, -6.02507901232996E-02, -3.07622221350501, -5.74011959864879E-02, 5.03471360939849, -0.925081888584834, 3.91733882917546, -77.314600713019, 9493.08762098587, -1410437.19679409, 8491662.30819026, 0.861095729446704, 0.32334644281172, 0.873281936020439, -0.436653048526683, 0.286596714529479, -0.131778331276228, 6.76682064330275E-03]
1440 hs = h / 2800.0
1441 ps = p / 100.0
1442 Ts = 0.0
1443 for i in range(0,33):
1444 Ts = Ts + ni[i] * (ps + 0.298) ** Ii[i] * (hs - 0.72) ** Ji[i]
1445 return Ts * 860.0
1446
1447
1448 def v3_ph(p, h):
1449 #%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
1450 #%2004
1451 #%Section 3.3 Backward Equations T(p,h) and v(p,h) for Subregions 3a and 3b
1452 #%Boundary equation, Eq 1 Page 5
1453 h3ab = 2014.64004206875 + 3.74696550136983 * p - 2.19921901054187E-02 * p ** 2 + 8.7513168600995E-05 * p ** 3
1454
1455 if h < h3ab :
1456 #%Subregion 3a
1457 #%Eq 4, Table 6, Page 9
1458 Ii = [-12, -12, -12, -12, -10, -10, -10, -8, -8, -6, -6, -6, -4, -4, -3, -2, -2, -1, -1, -1, -1, 0, 0, 1, 1, 1, 2, 2, 3, 4, 5, 8]
1459 Ji = [6, 8, 12, 18, 4, 7, 10, 5, 12, 3, 4, 22, 2, 3, 7, 3, 16, 0, 1, 2, 3, 0, 1, 0, 1, 2, 0, 2, 0, 2, 2, 2]
1460 ni = [5.29944062966028E-03, -0.170099690234461, 11.1323814312927, -2178.98123145125, -5.06061827980875E-04, 0.556495239685324, -9.43672726094016, -0.297856807561527, 93.9353943717186, 1.92944939465981E-02, 0.421740664704763, -3689141.2628233, -7.37566847600639E-03, -0.354753242424366, -1.99768169338727, 1.15456297059049, 5683.6687581596, 8.08169540124668E-03, 0.172416341519307, 1.04270175292927, -0.297691372792847, 0.560394465163593, 0.275234661176914, -0.148347894866012, -6.51142513478515E-02, -2.92468715386302, 6.64876096952665E-02, 3.52335014263844, -1.46340792313332E-02, -2.24503486668184, 1.10533464706142, -4.08757344495612E-02]
1461 ps = p / 100.0
1462 hs = h / 2100.0
1463 vs = 0.0
1464 for i in range(0,32):
1465 vs = vs + ni[i] * (ps + 0.128) ** Ii[i] * (hs - 0.727) ** Ji[i]
1466 return vs * 0.0028
1467 else:
1468 #%Subregion 3b
1469 #%Eq 5, Table 7, Page 9
1470 Ii = [-12, -12, -8, -8, -8, -8, -8, -8, -6, -6, -6, -6, -6, -6, -4, -4, -4, -3, -3, -2, -2, -1, -1, -1, -1, 0, 1, 1, 2, 2]
1471 Ji = [0, 1, 0, 1, 3, 6, 7, 8, 0, 1, 2, 5, 6, 10, 3, 6, 10, 0, 2, 1, 2, 0, 1, 4, 5, 0, 0, 1, 2, 6]
1472 ni = [-2.25196934336318E-09, 1.40674363313486E-08, 2.3378408528056E-06, -3.31833715229001E-05, 1.07956778514318E-03, -0.271382067378863, 1.07202262490333, -0.853821329075382, -2.15214194340526E-05, 7.6965608822273E-04, -4.31136580433864E-03, 0.453342167309331, -0.507749535873652, -100.475154528389, -0.219201924648793, -3.21087965668917, 607.567815637771, 5.57686450685932E-04, 0.18749904002955, 9.05368030448107E-03, 0.285417173048685, 3.29924030996098E-02, 0.239897419685483, 4.82754995951394, -11.8035753702231, 0.169490044091791, -1.79967222507787E-02, 3.71810116332674E-02, -5.36288335065096E-02, 1.6069710109252]
1473 ps = p / 100.0
1474 hs = h / 2800.0
1475 vs = 0.0
1476 for i in range(0,30):
1477 vs = vs + ni[i] * (ps + 0.0661) ** Ii[i] * (hs - 0.72) ** Ji[i]
1478 return vs * 0.0088
1479
1480
1481 def T3_ps(p, s):
1482 #%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
1483 #%2004
1484 #%3.4 Backward Equations T(p,s) and v(p,s) for Subregions 3a and 3b
1485 #%Boundary equation, Eq 6 Page 11
1486 if s <= 4.41202148223476 :
1487 #%Subregion 3a
1488 #%Eq 6, Table 10, Page 11
1489 Ii = [-12, -12, -10, -10, -10, -10, -8, -8, -8, -8, -6, -6, -6, -5, -5, -5, -4, -4, -4, -2, -2, -1, -1, 0, 0, 0, 1, 2, 2, 3, 8, 8, 10]
1490 Ji = [28, 32, 4, 10, 12, 14, 5, 7, 8, 28, 2, 6, 32, 0, 14, 32, 6, 10, 36, 1, 4, 1, 6, 0, 1, 4, 0, 0, 3, 2, 0, 1, 2]
1491 ni = [1500420082.63875, -159397258480.424, 5.02181140217975E-04, -67.2057767855466, 1450.58545404456, -8238.8953488889, -0.154852214233853, 11.2305046746695, -29.7000213482822, 43856513263.5495, 1.37837838635464E-03, -2.97478527157462, 9717779473494.13, -5.71527767052398E-05, 28830.794977842, -74442828926270.3, 12.8017324848921, -368.275545889071, 6.64768904779177E+15, 0.044935925195888, -4.22897836099655, -0.240614376434179, -4.74341365254924, 0.72409399912611, 0.923874349695897, 3.99043655281015, 3.84066651868009E-02, -3.59344365571848E-03, -0.735196448821653, 0.188367048396131, 1.41064266818704E-04, -2.57418501496337E-03, 1.23220024851555E-03]
1492 Sigma = s / 4.4
1493 Pi = p / 100.0
1494 teta = 0.0
1495 for i in range(0,33):
1496 teta = teta + ni[i] * (Pi + 0.24) ** Ii[i] * (Sigma - 0.703) ** Ji[i]
1497 return teta * 760.0
1498 else:
1499 #%Subregion 3b
1500 #%Eq 7, Table 11, Page 11
1501 Ii = [-12, -12, -12, -12, -8, -8, -8, -6, -6, -6, -5, -5, -5, -5, -5, -4, -3, -3, -2, 0, 2, 3, 4, 5, 6, 8, 12, 14]
1502 Ji = [1, 3, 4, 7, 0, 1, 3, 0, 2, 4, 0, 1, 2, 4, 6, 12, 1, 6, 2, 0, 1, 1, 0, 24, 0, 3, 1, 2]
1503 ni = [0.52711170160166, -40.1317830052742, 153.020073134484, -2247.99398218827, -0.193993484669048, -1.40467557893768, 42.6799878114024, 0.752810643416743, 22.6657238616417, -622.873556909932, -0.660823667935396, 0.841267087271658, -25.3717501764397, 485.708963532948, 880.531517490555, 2650155.92794626, -0.359287150025783, -656.991567673753, 2.41768149185367, 0.856873461222588, 0.655143675313458, -0.213535213206406, 5.62974957606348E-03, -316955725450471, -6.99997000152457E-04, 1.19845803210767E-02, 1.93848122022095E-05, -2.15095749182309E-05]
1504 Sigma = s / 5.3
1505 Pi = p / 100.0
1506 teta = 0.0
1507 for i in range(0,28):
1508 teta = teta + ni[i] * (Pi + 0.76) ** Ii[i] * (Sigma - 0.818) ** Ji[i]
1509 return teta * 860.0
1510
1511
1512 def v3_ps(p, s):
1513 #%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
1514 #%2004
1515 #%3.4 Backward Equations T(p,s) and v(p,s) for Subregions 3a and 3b
1516 #%Boundary equation, Eq 6 Page 11
1517 if s <= 4.41202148223476 :
1518 #%Subregion 3a
1519 #%Eq 8, Table 13, Page 14
1520 Ii = [-12, -12, -12, -10, -10, -10, -10, -8, -8, -8, -8, -6, -5, -4, -3, -3, -2, -2, -1, -1, 0, 0, 0, 1, 2, 4, 5, 6]
1521 Ji = [10, 12, 14, 4, 8, 10, 20, 5, 6, 14, 16, 28, 1, 5, 2, 4, 3, 8, 1, 2, 0, 1, 3, 0, 0, 2, 2, 0]
1522 ni = [79.5544074093975, -2382.6124298459, 17681.3100617787, -1.10524727080379E-03, -15.3213833655326, 297.544599376982, -35031520.6871242, 0.277513761062119, -0.523964271036888, -148011.182995403, 1600148.99374266, 1708023226634.27, 2.46866996006494E-04, 1.6532608479798, -0.118008384666987, 2.537986423559, 0.965127704669424, -28.2172420532826, 0.203224612353823, 1.10648186063513, 0.52612794845128, 0.277000018736321, 1.08153340501132, -7.44127885357893E-02, 1.64094443541384E-02, -6.80468275301065E-02, 0.025798857610164, -1.45749861944416E-04]
1523 Pi = p / 100.0
1524 Sigma = s / 4.4
1525 omega = 0.0
1526 for i in range(0,28):
1527 omega = omega + ni[i] * (Pi + 0.187) ** Ii[i] * (Sigma - 0.755) ** Ji[i]
1528 return omega * 0.0028
1529 else:
1530 #%Subregion 3b
1531 #%Eq 9, Table 14, Page 14
1532 Ii = [-12, -12, -12, -12, -12, -12, -10, -10, -10, -10, -8, -5, -5, -5, -4, -4, -4, -4, -3, -2, -2, -2, -2, -2, -2, 0, 0, 0, 1, 1, 2]
1533 Ji = [0, 1, 2, 3, 5, 6, 0, 1, 2, 4, 0, 1, 2, 3, 0, 1, 2, 3, 1, 0, 1, 2, 3, 4, 12, 0, 1, 2, 0, 2, 2]
1534 ni = [5.91599780322238E-05, -1.85465997137856E-03, 1.04190510480013E-02, 5.9864730203859E-03, -0.771391189901699, 1.72549765557036, -4.67076079846526E-04, 1.34533823384439E-02, -8.08094336805495E-02, 0.508139374365767, 1.28584643361683E-03, -1.63899353915435, 5.86938199318063, -2.92466667918613, -6.14076301499537E-03, 5.76199014049172, -12.1613320606788, 1.67637540957944, -7.44135838773463, 3.78168091437659E-02, 4.01432203027688, 16.0279837479185, 3.17848779347728, -3.58362310304853, -1159952.60446827, 0.199256573577909, -0.122270624794624, -19.1449143716586, -1.50448002905284E-02, 14.6407900162154, -3.2747778718823]
1535 Pi = p / 100.0
1536 Sigma = s / 5.3
1537 omega = 0.0
1538 for i in range(0,31):
1539 omega = omega + ni[i] * (Pi + 0.298) ** Ii[i] * (Sigma - 0.816) ** Ji[i]
1540 return omega * 0.0088
1541
1542
1543 def p3_hs(h, s):
1544 #%Supplementary Release on Backward Equations ( ) , p h s for Region 3,
1545 #%Equations as a function of h and s for the Region Boundaries, and an Equation
1546 #%( ) sat , T hs for Region 4 of the IAPWS Industrial formulation 1997 for the
1547 #%Thermodynamic Properties of Water and Steam
1548 #%2004
1549 #%Section 3 Backward functions p(h,s), T(h,s), and v(h,s) for Region 3
1550 if s < 4.41202148223476 :
1551 #%Subregion 3a
1552 #%Eq 1, Table 3, Page 8
1553 Ii = [0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 6, 7, 8, 10, 10, 14, 18, 20, 22, 22, 24, 28, 28, 32, 32]
1554 Ji = [0, 1, 5, 0, 3, 4, 8, 14, 6, 16, 0, 2, 3, 0, 1, 4, 5, 28, 28, 24, 1, 32, 36, 22, 28, 36, 16, 28, 36, 16, 36, 10, 28]
1555 ni = [7.70889828326934, -26.0835009128688, 267.416218930389, 17.2221089496844, -293.54233214597, 614.135601882478, -61056.2757725674, -65127225.1118219, 73591.9313521937, -11664650591.4191, 35.5267086434461, -596.144543825955, -475.842430145708, 69.6781965359503, 335.674250377312, 25052.6809130882, 146997.380630766, 5.38069315091534E+19, 1.43619827291346E+21, 3.64985866165994E+19, -2547.41561156775, 2.40120197096563E+27, -3.93847464679496E+29, 1.47073407024852E+24, -4.26391250432059E+31, 1.94509340621077E+38, 6.66212132114896E+23, 7.06777016552858E+33, 1.75563621975576E+41, 1.08408607429124E+28, 7.30872705175151E+43, 1.5914584739887E+24, 3.77121605943324E+40]
1556 Sigma = s / 4.4
1557 eta = h / 2300.0
1558 Pi = 0.0
1559 for i in range(0,33):
1560 Pi = Pi + ni[i] * (eta - 1.01) ** Ii[i] * (Sigma - 0.75) ** Ji[i]
1561 return Pi * 99.0
1562 else:
1563 #%Subregion 3b
1564 #%Eq 2, Table 4, Page 8
1565 Ii = [-12, -12, -12, -12, -12, -10, -10, -10, -10, -8, -8, -6, -6, -6, -6, -5, -4, -4, -4, -3, -3, -3, -3, -2, -2, -1, 0, 2, 2, 5, 6, 8, 10, 14, 14]
1566 Ji = [2, 10, 12, 14, 20, 2, 10, 14, 18, 2, 8, 2, 6, 7, 8, 10, 4, 5, 8, 1, 3, 5, 6, 0, 1, 0, 3, 0, 1, 0, 1, 1, 1, 3, 7]
1567 ni = [1.25244360717979E-13, -1.26599322553713E-02, 5.06878030140626, 31.7847171154202, -391041.161399932, -9.75733406392044E-11, -18.6312419488279, 510.973543414101, 373847.005822362, 2.99804024666572E-08, 20.0544393820342, -4.98030487662829E-06, -10.230180636003, 55.2819126990325, -206.211367510878, -7940.12232324823, 7.82248472028153, -58.6544326902468, 3550.73647696481, -1.15303107290162E-04, -1.75092403171802, 257.98168774816, -727.048374179467, 1.21644822609198E-04, 3.93137871762692E-02, 7.04181005909296E-03, -82.910820069811, -0.26517881813125, 13.7531682453991, -52.2394090753046, 2405.56298941048, -22736.1631268929, 89074.6343932567, -23923456.5822486, 5687958081.29714]
1568 Sigma = s / 5.3
1569 eta = h / 2800.0
1570 Pi = 0.0
1571 for i in range(0,35):
1572 Pi = Pi + ni[i] * (eta - 0.681) ** Ii[i] * (Sigma - 0.792) ** Ji[i]
1573 return 16.6 / float(Pi)
1574
1575
1576 def h3_pT(p, T):
1577 #%Not avalible with if 97
1578 #%Solve function T3_ph-T=0 with half interval method.
1579 #%ver2.6 Start corrected bug
1580 if p < 22.06395: #%Bellow tripple point
1581 Ts = T4_p(p) #%Saturation temperature
1582 if T <= Ts: #%Liquid side
1583 High_Bound = h4L_p(p) #%Max h ar liauid h.
1584 Low_Bound = h1_pT(p, 623.15)
1585 else:
1586 Low_Bound = h4V_p(p) #%Min h ar Vapour h.
1587 High_Bound = h2_pT(p, B23T_p(p))
1588 else : #%Above tripple point. R3 from R2 till R3.
1589 Low_Bound = h1_pT(p, 623.15)
1590 High_Bound = h2_pT(p, B23T_p(p))
1591
1592 #%ver2.6 End corrected bug
1593 Ts = T+1;
1594 while abs(T - Ts) > 0.00001 :
1595 hs = (Low_Bound + High_Bound) / 2.0
1596 Ts = T3_ph(p, hs)
1597 if Ts > T :
1598 High_Bound = hs
1599 else:
1600 Low_Bound = hs
1601 return hs
1602
1603
1604 def T3_prho(p, rho):
1605 #%Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
1606 #%Solve with half interval method
1607 Low_Bound = 623.15
1608 High_Bound = 1073.15
1609 ps=-1000.0
1610 while abs(p - ps) > 0.00000001:
1611 Ts = (Low_Bound + High_Bound) / 2.0
1612 ps = p3_rhoT(rho, Ts)
1613 if ps > p :
1614 High_Bound = Ts
1615 else:
1616 Low_Bound = Ts
1617 return Ts
1618
1619
1620
1621 #***********************************************************************************************************
1622 #*2.4 Functions for region 4
1623 #***********************************************************************************************************
1624 def p4_T(T):
1625 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1626 #%Section 8.1 The Saturation-Pressure Equation
1627 #%Eq 30, Page 33
1628 teta = T - 0.23855557567849 / (T - 650.17534844798)
1629 a = teta ** 2 + 1167.0521452767 * teta - 724213.16703206
1630 B = -17.073846940092 * teta ** 2 + 12020.82470247 * teta - 3232555.0322333
1631 C = 14.91510861353 * teta ** 2 - 4823.2657361591 * teta + 405113.40542057
1632 return (2 * C / (-B + (B ** 2.0 - 4.0 * a * C) ** 0.5)) ** 4.0
1633
1634
1635 def T4_p(p):
1636 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1637 #%Section 8.2 The Saturation-Temperature Equation
1638 #%Eq 31, Page 34
1639 beta = p**0.25
1640 E = beta ** 2 - 17.073846940092 * beta + 14.91510861353
1641 f = 1167.0521452767 * beta ** 2 + 12020.82470247 * beta - 4823.2657361591
1642 G = -724213.16703206 * beta ** 2 - 3232555.0322333 * beta + 405113.40542057
1643 D = 2 * G / (-f - (f ** 2 - 4 * E * G) ** 0.5)
1644 return (650.17534844798 + D - ((650.17534844798 + D) ** 2 - 4 * (-0.23855557567849 + 650.17534844798 * D)) ** 0.5) / 2
1645
1646
1647
1648 def h4_s(s):
1649 #%Supplementary Release on Backward Equations ( ) , p h s for Region 3,Equations as a function of h and s for the Region Boundaries, and an Equation( ) sat , T hs for Region 4 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
1650 #%4 Equations for Region Boundaries Given Enthalpy and Entropy
1651 #% Se picture page 14
1652 if s > -0.0001545495919 and s <= 3.77828134 :
1653 #%hL1_s
1654 #%Eq 3,Table 9,Page 16
1655 Ii = [0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 7, 8, 12, 12, 14, 14, 16, 20, 20, 22, 24, 28, 32, 32]
1656 Ji = [14, 36, 3, 16, 0, 5, 4, 36, 4, 16, 24, 18, 24, 1, 4, 2, 4, 1, 22, 10, 12, 28, 8, 3, 0, 6, 8]
1657 ni = [0.332171191705237, 6.11217706323496E-04, -8.82092478906822, -0.45562819254325, -2.63483840850452E-05, -22.3949661148062, -4.28398660164013, -0.616679338856916, -14.682303110404, 284.523138727299, -113.398503195444, 1156.71380760859, 395.551267359325, -1.54891257229285, 19.4486637751291, -3.57915139457043, -3.35369414148819, -0.66442679633246, 32332.1885383934, 3317.66744667084, -22350.1257931087, 5739538.75852936, 173.226193407919, -3.63968822121321E-02, 8.34596332878346E-07, 5.03611916682674, 65.5444787064505]
1658 Sigma = s / 3.8
1659 eta = 0.0
1660 for i in range(0,27):
1661 eta = eta + ni[i] * (Sigma - 1.09) ** Ii[i] * (Sigma + 0.0000366) ** Ji[i]
1662 return eta * 1700.0
1663 elif s > 3.77828134 and s <= 4.41202148223476:
1664 #%hL3_s
1665 #%Eq 4,Table 10,Page 16
1666 Ii = [0, 0, 0, 0, 2, 3, 4, 4, 5, 5, 6, 7, 7, 7, 10, 10, 10, 32, 32]
1667 Ji = [1, 4, 10, 16, 1, 36, 3, 16, 20, 36, 4, 2, 28, 32, 14, 32, 36, 0, 6]
1668 ni = [0.822673364673336, 0.181977213534479, -0.011200026031362, -7.46778287048033E-04, -0.179046263257381, 4.24220110836657E-02, -0.341355823438768, -2.09881740853565, -8.22477343323596, -4.99684082076008, 0.191413958471069, 5.81062241093136E-02, -1655.05498701029, 1588.70443421201, -85.0623535172818, -31771.4386511207, -94589.0406632871, -1.3927384708869E-06, 0.63105253224098]
1669 Sigma = s / 3.8
1670 eta = 0.0
1671 for i in range(0,19):
1672 eta = eta + ni[i] * (Sigma - 1.09) ** Ii[i] * (Sigma + 0.0000366) ** Ji[i]
1673 return eta * 1700.0
1674 elif s > 4.41202148223476 and s <= 5.85 :
1675 #%Section 4.4 Equations ( ) 2ab " h s and ( ) 2c3b "h s for the Saturated Vapor Line
1676 #%Page 19, Eq 5
1677 #%hV2c3b_s(s)
1678 Ii = [0, 0, 0, 1, 1, 5, 6, 7, 8, 8, 12, 16, 22, 22, 24, 36]
1679 Ji = [0, 3, 4, 0, 12, 36, 12, 16, 2, 20, 32, 36, 2, 32, 7, 20]
1680 ni = [1.04351280732769, -2.27807912708513, 1.80535256723202, 0.420440834792042, -105721.24483466, 4.36911607493884E+24, -328032702839.753, -6.7868676080427E+15, 7439.57464645363, -3.56896445355761E+19, 1.67590585186801E+31, -3.55028625419105E+37, 396611982166.538, -4.14716268484468E+40, 3.59080103867382E+18, -1.16994334851995E+40]
1681 Sigma = s / 5.9
1682 eta = 0.0
1683 for i in range(0,16):
1684 eta = eta + ni[i] * (Sigma - 1.02) ** Ii[i] * (Sigma - 0.726) ** Ji[i]
1685 return eta ** 4 * 2800.0
1686 elif s > 5.85 and s < 9.155759395:
1687 #%Section 4.4 Equations ( ) 2ab " h s and ( ) 2c3b "h s for the Saturated Vapor Line
1688 #%Page 20, Eq 6
1689 Ii = [1, 1, 2, 2, 4, 4, 7, 8, 8, 10, 12, 12, 18, 20, 24, 28, 28, 28, 28, 28, 32, 32, 32, 32, 32, 36, 36, 36, 36, 36]
1690 Ji = [8, 24, 4, 32, 1, 2, 7, 5, 12, 1, 0, 7, 10, 12, 32, 8, 12, 20, 22, 24, 2, 7, 12, 14, 24, 10, 12, 20, 22, 28]
1691 ni = [-524.581170928788, -9269472.18142218, -237.385107491666, 21077015581.2776, -23.9494562010986, 221.802480294197, -5104725.33393438, 1249813.96109147, 2000084369.96201, -815.158509791035, -157.612685637523, -11420042233.2791, 6.62364680776872E+15, -2.27622818296144E+18, -1.71048081348406E+31, 6.60788766938091E+15, 1.66320055886021E+22, -2.18003784381501E+29, -7.87276140295618E+29, 1.51062329700346E+31, 7957321.70300541, 1.31957647355347E+15, -3.2509706829914E+23, -4.18600611419248E+25, 2.97478906557467E+34, -9.53588761745473E+19, 1.66957699620939E+24, -1.75407764869978E+32, 3.47581490626396E+34, -7.10971318427851E+38]
1692 Sigma1 = s / 5.21
1693 Sigma2 = s / 9.2
1694 eta = 0.0
1695 for i in range(0,30):
1696 eta = eta + ni[i] * (1 / Sigma1 - 0.513) ** Ii[i] * (Sigma2 - 0.524) ** Ji[i]
1697 return math.exp(eta) * 2800.0
1698 else:
1699 return -99999.9
1700
1701
1702 def p4_s(s):
1703 #%Uses h4_s and p_hs for the diffrent regions to determine p4_s
1704 h_sat = h4_s(s)
1705 if s > -0.0001545495919 and s <= 3.77828134 :
1706 return p1_hs(h_sat, s)
1707 elif s > 3.77828134 and s <= 5.210887663 :
1708 return p3_hs(h_sat, s)
1709 elif s > 5.210887663 and s < 9.155759395 :
1710 return p2_hs(h_sat, s)
1711 else:
1712 return -99999.9
1713
1714
1715 def h4L_p(p):
1716 if p > 0.000611657 and p < 22.06395 :
1717 Ts = T4_p(p)
1718 if p < 16.529 :
1719 return h1_pT(p, Ts)
1720 else:
1721 #%Iterate to find the the backward solution of p3sat_h
1722 Low_Bound = 1670.858218
1723 High_Bound = 2087.23500164864
1724 ps=-1000.0
1725 while abs(p - ps) > 0.00001 :
1726 hs = (Low_Bound + High_Bound) / 2.0
1727 ps = p3sat_h(hs)
1728 if ps > p :
1729 High_Bound = hs
1730 else:
1731 Low_Bound = hs
1732 return hs
1733 else:
1734 return -99999.9
1735
1736
1737 def h4V_p(p):
1738 if p > 0.000611657 and p < 22.06395:
1739 Ts = T4_p(p)
1740 if p < 16.529:
1741 return h2_pT(p, Ts)
1742 else:
1743 #%Iterate to find the the backward solution of p3sat_h
1744 Low_Bound = 2087.23500164864
1745 High_Bound = 2563.592004+5
1746 ps=-1000.0
1747 while abs(p - ps) > 0.000001 :
1748 hs = (Low_Bound + High_Bound) / 2.0
1749 ps = p3sat_h(hs)
1750 if ps < p :
1751 High_Bound = hs
1752 else :
1753 Low_Bound = hs
1754 return hs
1755 else :
1756 return -99999.9
1757
1758
1759 def x4_ph(p, h):
1760 #%Calculate vapour fraction from hL and hV for given p
1761 h4v = h4V_p(p)
1762 h4L = h4L_p(p)
1763 if h > h4v :
1764 return float('NaN')
1765 #return 1
1766 elif h < h4L :
1767 return float('NaN')
1768 #return 0
1769 else :
1770 return (h - h4L) / (h4v - h4L)
1771
1772
1773 def x4_ps(p, s):
1774 if p < 16.529 :
1775 ssv = s2_pT(p, T4_p(p))
1776 ssL = s1_pT(p, T4_p(p))
1777 else:
1778 ssv = s3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p))
1779 ssL = s3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p))
1780 if s < ssL :
1781 return float('NaN')
1782 #return 0
1783 elif s > ssv :
1784 return float('NaN')
1785 #return 1
1786 else:
1787 return (s - ssL) / (ssv - ssL)
1788
1789
1790 def T4_hs(h, s):
1791 #%Supplementary Release on Backward Equations ( ) , p h s for Region 3,
1792 #%Chapter 5.3 page 30.
1793 #%The if 97 function is only valid for part of region4. Use iteration outsida.
1794 Ii = [0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 5, 5, 5, 5, 6, 6, 6, 8, 10, 10, 12, 14, 14, 16, 16, 18, 18, 18, 20, 28]
1795 Ji = [0, 3, 12, 0, 1, 2, 5, 0, 5, 8, 0, 2, 3, 4, 0, 1, 1, 2, 4, 16, 6, 8, 22, 1, 20, 36, 24, 1, 28, 12, 32, 14, 22, 36, 24, 36]
1796 ni = [0.179882673606601, -0.267507455199603, 1.162767226126, 0.147545428713616, -0.512871635973248, 0.421333567697984, 0.56374952218987, 0.429274443819153, -3.3570455214214, 10.8890916499278, -0.248483390456012, 0.30415322190639, -0.494819763939905, 1.07551674933261, 7.33888415457688E-02, 1.40170545411085E-02, -0.106110975998808, 1.68324361811875E-02, 1.25028363714877, 1013.16840309509, -1.51791558000712, 52.4277865990866, 23049.5545563912, 2.49459806365456E-02, 2107964.67412137, 366836848.613065, -144814105.365163, -1.7927637300359E-03, 4899556021.00459, 471.262212070518, -82929439019.8652, -1715.45662263191, 3557776.82973575, 586062760258.436, -12988763.5078195, 31724744937.1057]
1797 if s > 5.210887825 and s < 9.15546555571324:
1798 Sigma = s / 9.2
1799 eta = h / 2800.0
1800 teta = 0.0
1801 for i in range(0,36):
1802 teta = teta + ni[i] * (eta - 0.119) ** Ii[i] * (Sigma - 1.07) ** Ji[i]
1803 return teta * 550.0
1804 else:
1805 #%function psat_h
1806 if s > -0.0001545495919 and s <= 3.77828134:
1807 Low_Bound = 0.000611
1808 High_Bound = 165.291642526045
1809 hL=-1000.0
1810 while abs(hL - h) > 0.00001 and abs(High_Bound - Low_Bound) > 0.0001:
1811 PL = (Low_Bound + High_Bound) / 2.0
1812 Ts = T4_p(PL)
1813 hL = h1_pT(PL, Ts)
1814 if hL > h :
1815 High_Bound = PL
1816 else :
1817 Low_Bound = PL
1818 elif s > 3.77828134 and s <= 4.41202148223476 :
1819 PL = p3sat_h(h)
1820 elif s > 4.41202148223476 and s <= 5.210887663 :
1821 PL = p3sat_h(h)
1822 Low_Bound = 0.000611
1823 High_Bound = float(PL)
1824 sss=-1000.0
1825 while abs(s - sss) > 0.000001 and abs(High_Bound - Low_Bound) > 0.0000001:
1826 p = float((Low_Bound + High_Bound) / 2.0)
1827 #%Calculate s4_ph
1828 Ts = T4_p(p)
1829 xs = x4_ph(p, h)
1830 if p < 16.529 :
1831 s4v = s2_pT(p, Ts)
1832 s4L = s1_pT(p, Ts)
1833
1834 else:
1835 v4v = v3_ph(p, h4V_p(p))
1836 s4v = s3_rhoT(1.0 / v4v, Ts)
1837 v4L = v3_ph(p, h4L_p(p))
1838 s4L = s3_rhoT(1.0 / v4L, Ts)
1839
1840 sss = (xs * s4v + (1 - xs) * s4L)
1841 if sss < s :
1842 High_Bound = p
1843 else :
1844 Low_Bound = p
1845 return T4_p(p)
1846
1847
1848
1849 #***********************************************************************************************************
1850 #*2.5 Functions for region 5
1851 #***********************************************************************************************************
1852 def h5_pT(p, T):
1853 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1854 #%Basic Equation for Region 5
1855 #%Eq 32,33, Page 36, Tables 37-41
1856 Ji0 = [0, 1, -3, -2, -1, 2]
1857 ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917]
1858 Iir = [1, 1, 1, 2, 3]
1859 Jir = [0, 1, 3, 9, 3]
1860 nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07]
1861 R = 0.461526 #%kJ/(kg K)
1862 tau = 1000 / float(T)
1863 Pi = float(p)
1864 gamma0_tau = 0.0
1865 for i in range(0,6):
1866 gamma0_tau = gamma0_tau + ni0[i] * Ji0[i] * tau ** (Ji0[i] - 1)
1867 gammar_tau = 0.0
1868 for i in range(0,5):
1869 gammar_tau = gammar_tau + nir[i] * Pi ** Iir[i] * Jir[i] * tau ** (Jir[i] - 1)
1870 return R * T * tau * (gamma0_tau + gammar_tau)
1871
1872
1873 def v5_pT(p, T):
1874 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1875 #%Basic Equation for Region 5
1876 #%Eq 32,33, Page 36, Tables 37-41
1877 Ji0 = [0, 1, -3, -2, -1, 2]
1878 ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917]
1879 Iir = [1, 1, 1, 2, 3]
1880 Jir = [0, 1, 3, 9, 3]
1881 nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07]
1882 R = 0.461526 #%kJ/(kg K)
1883 tau = 1000 / float(T)
1884 Pi = float(p)
1885 gamma0_pi = 1 / float(Pi)
1886 gammar_pi = 0.0
1887 for i in range(0,5):
1888 gammar_pi = gammar_pi + nir[i] * Iir[i] * Pi ** (Iir[i] - 1) * tau ** Jir[i]
1889 return R * T / p * Pi * (gamma0_pi + gammar_pi) / 1000.0
1890
1891
1892 def u5_pT(p, T):
1893 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1894 #%Basic Equation for Region 5
1895 #%Eq 32,33, Page 36, Tables 37-41
1896 Ji0 = [0, 1, -3, -2, -1, 2]
1897 ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917]
1898 Iir = [1, 1, 1, 2, 3]
1899 Jir = [0, 1, 3, 9, 3]
1900 nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07]
1901 R = 0.461526 #%kJ/(kg K)
1902 tau = 1000 / float(T)
1903 Pi = float(p)
1904 gamma0_pi = 1 / float(Pi)
1905 gamma0_tau = 0.0
1906 for i in range(0,6):
1907 gamma0_tau = gamma0_tau + ni0[i] * Ji0[i] * tau ** (Ji0[i] - 1)
1908 gammar_pi = 0.0
1909 gammar_tau = 0.0
1910 for i in range(0,5):
1911 gammar_pi = gammar_pi + nir[i] * Iir[i] * Pi ** (Iir[i] - 1) * tau ** Jir[i]
1912 gammar_tau = gammar_tau + nir[i] * Pi ** Iir[i] * Jir[i] * tau ** (Jir[i] - 1)
1913 return R * T * (tau * (gamma0_tau + gammar_tau) - Pi * (gamma0_pi + gammar_pi))
1914
1915
1916 def Cp5_pT(p, T):
1917 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1918 #%Basic Equation for Region 5
1919 #%Eq 32,33, Page 36, Tables 37-41
1920 Ji0 = [0, 1, -3, -2, -1, 2]
1921 ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917]
1922 Iir = [1, 1, 1, 2, 3]
1923 Jir = [0, 1, 3, 9, 3]
1924 nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07]
1925 R = 0.461526 #%kJ/(kg K)
1926 tau = 1000 / float(T)
1927 Pi = float(p)
1928 gamma0_tautau = 0.0
1929 for i in range(0,6):
1930 gamma0_tautau = gamma0_tautau + ni0[i] * Ji0[i] * (Ji0[i] - 1) * tau ** (Ji0[i] - 2)
1931 gammar_tautau = 0.0
1932 for i in range(0,5):
1933 gammar_tautau = gammar_tautau + nir[i] * Pi ** Iir[i] * Jir[i] * (Jir[i] - 1) * tau ** (Jir[i] - 2)
1934 return -R * tau ** 2.0 * (gamma0_tautau + gammar_tautau)
1935
1936
1937 def s5_pT(p, T):
1938 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1939 #%Basic Equation for Region 5
1940 #%Eq 32,33, Page 36, Tables 37-41
1941 Ji0 = [0, 1, -3, -2, -1, 2]
1942 ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917]
1943 Iir = [1, 1, 1, 2, 3]
1944 Jir = [0, 1, 3, 9, 3]
1945 nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07]
1946 R = 0.461526 #%kJ/(kg K)
1947 tau = 1000 / float(T)
1948 Pi = float(p)
1949 gamma0 = log(Pi)
1950 gamma0_tau = 0.0
1951 for i in range(0,6):
1952 gamma0_tau = gamma0_tau + ni0[i] * Ji0[i] * tau ** (Ji0[i] - 1)
1953 gamma0 = gamma0 + ni0[i] * tau ** Ji0[i]
1954 gammar = 0.0
1955 gammar_tau = 0.0
1956 for i in range(0,5):
1957 gammar = gammar + nir[i] * Pi ** Iir[i] * tau ** Jir[i]
1958 gammar_tau = gammar_tau + nir[i] * Pi ** Iir[i] * Jir[i] * tau ** (Jir[i] - 1)
1959 return R * (tau * (gamma0_tau + gammar_tau) - (gamma0 + gammar))
1960
1961
1962 def Cv5_pT(p, T):
1963 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1964 #%Basic Equation for Region 5
1965 #%Eq 32,33, Page 36, Tables 37-41
1966 Ji0 = [0, 1, -3, -2, -1, 2]
1967 ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917]
1968 Iir = [1, 1, 1, 2, 3]
1969 Jir = [0, 1, 3, 9, 3]
1970 nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07]
1971 R = 0.461526 #%kJ/(kg K)
1972 tau = 1000 / float(T)
1973 Pi = float(p)
1974 gamma0_tautau = 0.0
1975 for i in range(0,6):
1976 gamma0_tautau = gamma0_tautau + ni0[i] * (Ji0[i] - 1) * Ji0[i] * tau ** (Ji0[i] - 2)
1977 gammar_pi = 0.0
1978 gammar_pitau = 0.0
1979 gammar_pipi = 0.0
1980 gammar_tautau = 0.0
1981 for i in range(0,5):
1982 gammar_pi = gammar_pi + nir[i] * Iir[i] * Pi ** (Iir[i] - 1) * tau ** Jir[i]
1983 gammar_pitau = gammar_pitau + nir[i] * Iir[i] * Pi ** (Iir[i] - 1) * Jir[i] * tau ** (Jir[i] - 1)
1984 gammar_pipi = gammar_pipi + nir[i] * Iir[i] * (Iir[i] - 1) * Pi ** (Iir[i] - 2) * tau ** Jir[i]
1985 gammar_tautau = gammar_tautau + nir[i] * Pi ** Iir[i] * Jir[i] * (Jir[i] - 1) * tau ** (Jir[i] - 2)
1986 return R * (-(tau ** 2 *(gamma0_tautau + gammar_tautau)) - (1 + Pi * gammar_pi - tau * Pi * gammar_pitau)**2 / (1 - Pi ** 2 * gammar_pipi))
1987
1988
1989 def w5_pT(p, T):
1990 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
1991 #%Basic Equation for Region 5
1992 #%Eq 32,33, Page 36, Tables 37-41
1993 Ji0 = [0, 1, -3, -2, -1, 2]
1994 ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917]
1995 Iir = [1, 1, 1, 2, 3]
1996 Jir = [0, 1, 3, 9, 3]
1997 nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07]
1998 R = 0.461526 #%kJ/(kg K)
1999 tau = 1000 / float(T)
2000 Pi = float(p)
2001 gamma0_tautau = 0.0
2002 for i in range(0,6):
2003 gamma0_tautau = gamma0_tautau + ni0[i] * (Ji0[i] - 1) * Ji0[i] * tau ** (Ji0[i] - 2)
2004 gammar_pi = 0.0
2005 gammar_pitau = 0.0
2006 gammar_pipi = 0.0
2007 gammar_tautau = 0.0
2008 for i in range(0,5):
2009 gammar_pi = gammar_pi + nir[i] * Iir[i] * Pi ** (Iir[i] - 1) * tau ** Jir[i]
2010 gammar_pitau = gammar_pitau + nir[i] * Iir[i] * Pi ** (Iir[i] - 1) * Jir[i] * tau ** (Jir[i] - 1)
2011 gammar_pipi = gammar_pipi + nir[i] * Iir[i] * (Iir[i] - 1) * Pi ** (Iir[i] - 2) * tau ** Jir[i]
2012 gammar_tautau = gammar_tautau + nir[i] * Pi ** Iir[i] * Jir[i] * (Jir[i] - 1) * tau ** (Jir[i] - 2)
2013 return (1000.0 * R * T * (1 + 2 * Pi * gammar_pi + Pi ** 2 * gammar_pi ** 2) / ((1 - Pi ** 2 * gammar_pipi) + (1 + Pi * gammar_pi - tau * Pi * gammar_pitau) ** 2 / (tau ** 2 * (gamma0_tautau + gammar_tautau)))) ** 0.5
2014
2015
2016 def T5_ph(p, h):
2017 #%Solve with half interval method
2018 Low_Bound = 1073.15
2019 High_Bound = 2273.15
2020 hs=h-1
2021 while abs(h - hs) > 0.00001 :
2022 Ts = (Low_Bound + High_Bound) / 2.0
2023 hs = h5_pT(p, Ts)
2024 if hs > h :
2025 High_Bound = Ts
2026 else:
2027 Low_Bound = Ts
2028 return Ts
2029
2030
2031 def T5_ps(p, s):
2032 #%Solve with half interval method
2033 Low_Bound = 1073.15
2034 High_Bound = 2273.15
2035 ss=s-1
2036 while abs(s - ss) > 0.00001:
2037 Ts = (Low_Bound + High_Bound) / 2.0
2038 ss = s5_pT(p, Ts)
2039 if ss > s :
2040 High_Bound = Ts
2041 else:
2042 Low_Bound = Ts
2043 return Ts
2044
2045
2046 def T5_prho(p,rho):
2047 #%Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
2048 #%Solve with half interval method
2049
2050 Low_Bound = 1073.15
2051 High_Bound = 2073.15
2052 rhos=-1000.0
2053 while abs(rho - rhos) > 0.000001:
2054 Ts = (Low_Bound + High_Bound) / 2.0
2055 rhos = 1 / v2_pT(p, Ts)
2056 if rhos < rho:
2057 High_Bound = Ts
2058 else:
2059 Low_Bound = Ts
2060 return Ts
2061
2062
2063
2064
2065 '''
2066 %***********************************************************************************************************
2067 %*3 Region Selection
2068 %***********************************************************************************************************
2069 '''
2070 #***********************************************************************************************************
2071 #*3.1 Regions as a function of pT
2072 #***********************************************************************************************************
2073 def region_pT(p,T):
2074 if T > 1073.15 and p < 10 and T < 2273.15 and p > 0.000611 :
2075 return 5
2076 elif T <= 1073.15 and T > 273.15 and p <= 100 and p > 0.000611:
2077 if T > 623.15:
2078 if p > B23p_T(T):
2079 if T < 647.096:
2080 ps = p4_T(T)
2081 if abs(p - ps) < 0.00001:
2082 return 4
2083 else :
2084 return 3
2085 else:
2086 return 2
2087 else :
2088 ps = p4_T(T)
2089 if abs(p-ps)<0.00001:
2090 return 4
2091 elif p>ps:
2092 return 1
2093 else:
2094 return 2
2095
2096 else:
2097 print 'Error, Outside valid area'
2098
2099
2100 '''
2101 MATLAB VERSION
2102 %*3.1 Regions as a function of pT
2103 function region_pT = region_pT(p, T)
2104 if T > 1073.15 && p < 10 && T < 2273.15 && p > 0.000611
2105 region_pT = 5;
2106 elseif T <= 1073.15 && T > 273.15 && p <= 100 && p > 0.000611
2107 if T > 623.15
2108 if p > B23p_T(T)
2109 region_pT = 3;
2110 if T < 647.096
2111 ps = p4_T(T);
2112 if abs(p - ps) < 0.00001
2113 region_pT = 4;
2114 end
2115 end
2116 else
2117 region_pT = 2;
2118 end
2119 else
2120 ps = p4_T(T);
2121 if abs(p - ps) < 0.00001
2122 region_pT = 4;
2123 elseif p > ps
2124 region_pT = 1;
2125 else
2126 region_pT = 2;
2127 end
2128 end
2129 else
2130 region_pT = 0; %**Error, Outside valid area
2131 end
2132 '''
2133
2134
2135 #***********************************************************************************************************
2136 #*3.2 Regions as a function of ph
2137 #***********************************************************************************************************
2138 def region_ph(p,h):
2139 #%Check if outside pressure limits
2140 if p < 0.000611657 or p > 100 :
2141 return 0
2142
2143 #%Check if outside low h.
2144 elif h < 0.963 * p + 2.2 : # %Linear adaption to h1_pt()+2 to speed up calcualations.
2145 if h < h1_pT(p, 273.15) :
2146 return 0
2147
2148 elif p < 16.5292 : # %Bellow region 3,Check region 1,4,2,5
2149 #%Check Region 1
2150 Ts = T4_p(p)
2151 hL = 109.6635 * log(p) + 40.3481 * p + 734.58 #%Approximate function for hL_p
2152 if abs(h - hL) < 100 : #%if approximate is not god enough use real function
2153 hL = h1_pT(p, Ts)
2154 if h <= hL :
2155 return 1
2156 #%Check Region 4
2157 hV = 45.1768 * log(p) - 20.158 * p + 2804.4 #%Approximate function for hV_p
2158 if abs(h - hV) < 50: #%if approximate is not god enough use real function
2159 hV = h2_pT(p, Ts)
2160 if h < hV :
2161 return 4
2162 #%Check upper limit of region 2 Quick Test
2163 if h < 4000.0 :
2164 return 2
2165 #%Check region 2 (Real value)
2166 h_45 = h2_pT(p, 1073.15)
2167 if h <= h_45 :
2168 return 2
2169 #%Check region 5
2170 h_5u = h5_pT(p, 2273.15)
2171 if p > 10.0 :
2172 return 0
2173 if h < h_5u :
2174 return 5
2175 else:
2176 return 0
2177
2178 else: # %for p>16.5292
2179 #%Check if in region1
2180 if h < h1_pT(p, 623.15) :
2181 return 1
2182 #%Check if in region 3 or 4 (Bellow Reg 2)
2183 elif h < h2_pT(p, B23T_p(p)) :
2184 #%Region 3 or 4
2185 if p > p3sat_h(h) :
2186 return 3
2187 else :
2188 return 4
2189 #%Check if region 2
2190 elif h < h2_pT(p, 1073.15) :
2191 return 2
2192 else :
2193 return 0
2194
2195
2196 #***********************************************************************************************************
2197 #*3.3 Regions as a function of ps
2198 #***********************************************************************************************************
2199 def region_ps(p,s):
2200 if p < 0.000611657 or p > 100.0 or s < 0 or (s > s5_pT(p, 2273.15) and s5_pT(p, 2273.15) != float('NaN')):
2201 return 0
2202
2203 #%Check region 5
2204 if s > s2_pT(p, 1073.15) and s2_pT(p,1073.15) != float('NaN') :
2205 if p <= 10.0 :
2206 return 5
2207 else:
2208 return 0
2209
2210 #%Check region 2
2211 if p > 16.529 :
2212 ss = s2_pT(p, B23T_p(p)) #%Between 5.047 & 5.261. Use to speed up!
2213 else:
2214 ss = s2_pT(p, T4_p(p))
2215 if ss != float('NaN') and s > ss :
2216 return 2
2217
2218 #%Check region 3
2219 ss = s1_pT(p, 623.15)
2220 if p > 16.529 and s > ss and ss != float('NaN') :
2221 if p > p3sat_s(s) and p3sat_s(s) != float('NaN') :
2222 return 3
2223 else:
2224 return 4
2225
2226 #%Check region 4 (Not inside region 3)
2227 if p < 16.529 and s > s1_pT(p, T4_p(p)) and s1_pT(p, T4_p(p)) != float('NaN'):
2228 return 4
2229
2230 #%Check region 1
2231 if p > 0.000611657 and s > s1_pT(p, 273.15) and s1_pT(p, 273.15) != float('NaN') :
2232 return 1
2233
2234 return 1
2235
2236
2237 #***********************************************************************************************************
2238 #*3.4 Regions as a function of hs
2239 #***********************************************************************************************************
2240 def region_hs(h,s):
2241 if s < -0.0001545495919:
2242 return 0
2243
2244 #%Check linear adaption to p=0.000611. if bellow region 4.
2245 hMin = (((-0.0415878 - 2500.89262) / (-0.00015455 - 9.155759)) * s)
2246 if s < 9.155759395 and h < hMin :
2247 return 0
2248
2249 #%******Kolla 1 eller 4. (+liten bit over B13)
2250 if s >= -0.0001545495919 and s <= 3.77828134 :
2251 if h < h4_s(s) :
2252 return 4
2253 elif s < 3.397782955 : #%100MPa line is limiting
2254 TMax = T1_ps(100.0, s)
2255 hMax = h1_pT(100.0, TMax)
2256 if h < hMax :
2257 return 1
2258
2259 else:
2260 return 0
2261 else: #%The point is either in region 4,1,3. Check B23
2262 hB = hB13_s(s)
2263 if h < hB :
2264 return 1
2265 TMax = T3_ps(100.0, s)
2266 vmax = v3_ps(100.0, s)
2267 hMax = h3_rhoT(1.0 / vmax, TMax)
2268 if h < hMax :
2269 return 3
2270 else:
2271 return 0
2272
2273
2274 #%******Kolla region 2 eller 4. (Ovre delen av omrade b23-> max)
2275 if s >= 5.260578707 and s <= 11.9212156897728 :
2276 if s > 9.155759395 : #%Above region 4
2277 Tmin = T2_ps(0.000611, s)
2278 hMin = h2_pT(0.000611, Tmin)
2279 #%function adapted to h(1073.15,s)
2280 hMax = -0.07554022 * s ** 4 + 3.341571 * s ** 3 - 55.42151 * s ** 2 + 408.515 * s + 3031.338
2281 if h > hMin and h < hMax :
2282 return 2
2283 else:
2284 return 0
2285 hV = h4_s(s)
2286 if h < hV : #%Region 4. Under region 3.
2287 return 4
2288 if s < 6.04048367171238 :
2289 TMax = T2_ps(100.0, s)
2290 hMax = h2_pT(100.0, TMax)
2291 else:
2292 #%function adapted to h(1073.15,s)
2293 hMax = -2.988734 * s ** 4 + 121.4015 * s ** 3 - 1805.15 * s ** 2 + 11720.16 * s - 23998.33
2294 if h < hMax: #%Region 2. Over region 4.
2295 return 2
2296 else :
2297 return 0
2298
2299 #%Kolla region 3 eller 4. Under kritiska punkten.
2300 if s >= 3.77828134 and s <= 4.41202148223476 :
2301 hL = h4_s(s)
2302 if h < hL :
2303 return 4
2304 TMax = T3_ps(100.0, s)
2305 vmax = v3_ps(100.0, s)
2306 hMax = h3_rhoT(1.0 / vmax, TMax)
2307 if h < hMax :
2308 return 3
2309 else:
2310 return 0
2311
2312
2313 #%Kolla region 3 eller 4 fran kritiska punkten till ovre delen av b23
2314 if s >= 4.41202148223476 and s <= 5.260578707 :
2315 hV = h4_s(s)
2316 if h < hV :
2317 return 4
2318 #%Kolla om vi ar under b23 giltighetsomrade.
2319 if s <= 5.048096828 :
2320 TMax = T3_ps(100.0, s)
2321 vmax = v3_ps(100.0, s)
2322 hMax = h3_rhoT(1.0 / vmax, TMax)
2323 if h < hMax :
2324 return 3
2325 else:
2326 return 0
2327
2328 else : #%Inom omradet for B23 i s led.
2329 if (h > 2812.942061): #%Ovanfor B23 i h_led
2330 if s > 5.09796573397125 :
2331 TMax = T2_ps(100.0, s)
2332 hMax = h2_pT(100.0, TMax)
2333 if h < hMax :
2334 return 2
2335 else:
2336 return 0
2337 else:
2338 return 0
2339
2340 if (h < 2563.592004) : #%Nedanfor B23 i h_led men vi har redan kollat ovanfor hV2c3b
2341 return 3
2342
2343 #%Vi ar inom b23 omradet i bade s och h led.
2344 Tact = TB23_hs(h, s)
2345 pact = p2_hs(h, s)
2346 pBound = B23p_T(Tact)
2347 if pact > pBound :
2348 return 3
2349 else:
2350 return 2
2351 return 0
2352
2353
2354
2355
2356
2357
2358 #***********************************************************************************************************
2359 #*3.5 Regions as a function of p and rho
2360 #***********************************************************************************************************
2361 def Region_prho(p,rho):
2362 v = 1 / rho
2363 if p < 0.000611657 or p > 100.0:
2364 return 0
2365
2366 if p < 16.5292 : #%Bellow region 3, Check region 1,4,2
2367 if v < v1_pT(p, 273.15): #%Observe that this is not actually min of v. Not valid Water of 4 C is ligther.
2368 return 0
2369
2370 if v <= v1_pT(p, T4_p(p)):
2371 return 1
2372
2373 if v < v2_pT(p, T4_p(p)) :
2374 return 4
2375
2376 if v <= v2_pT(p, 1073.15):
2377 return 2
2378
2379 if p > 10.0: #%Above region 5
2380 return 0
2381
2382 if v <= v5_pT(p, 2073.15):
2383 return 5
2384
2385 else: #%Check region 1,3,4,3,2 (Above the lowest point of region 3.)
2386 if v < v1_pT(p, 273.15): #%Observe that this is not actually min of v. Not valid Water of 4C is ligther.
2387 return 0
2388
2389 if v < v1_pT(p, 623.15):
2390 return 1
2391
2392 #%Check if in region 3 or 4 (Bellow Reg 2)
2393 if v < v2_pT(p, B23T_p(p)):
2394 #%Region 3 or 4
2395 if p > 22.064: #%Above region 4
2396 return 3
2397
2398 if v < v3_ph(p, h4L_p(p)) or (v > v3_ph(p, h4V_p(p)) and v3_ph(p, h4V_p(p)) != float('NaN')): #%Uses iteration!!
2399 return 3
2400 else:
2401 return 4
2402
2403 #%Check if region 2
2404 if v < v2_pT(p, 1073.15):
2405 return 2
2406 return 0
2407
2408
2409
2410
2411 '''
2412 %***********************************************************************************************************
2413 %*4 Region Borders
2414 %***********************************************************************************************************
2415 '''
2416
2417 #***********************************************************************************************************
2418 #*4.1 Boundary between region 2 and 3.
2419 #***********************************************************************************************************
2420 def B23p_T(T):
2421 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
2422 #%1997
2423 #%Section 4 Auxiliary Equation for the Boundary between Regions 2 and 3
2424 #%Eq 5, Page 5
2425 return 348.05185628969 - 1.1671859879975 * T + 1.0192970039326E-03 * T ** 2.0
2426
2427
2428 def B23T_p(p):
2429 #%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
2430 #%1997
2431 #%Section 4 Auxiliary Equation for the Boundary between Regions 2 and 3
2432 #%Eq 6, Page 6
2433 return 572.54459862746 + ((p - 13.91883977887) / 1.0192970039326E-03) ** 0.5
2434
2435
2436
2437 #***********************************************************************************************************
2438 #*4.2 Region 3. pSat_h & pSat_s
2439 #***********************************************************************************************************
2440 def p3sat_h(h):
2441 #%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h) & T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water & Steam
2442 #%2004
2443 #%Section 4 Boundary Equations psat(h) & psat(s) for the Saturation Lines of Region 3
2444 #%Se pictures Page 17, Eq 10, Table 17, Page 18
2445 Ii = [0, 1, 1, 1, 1, 5, 7, 8, 14, 20, 22, 24, 28, 36]
2446 Ji = [0, 1, 3, 4, 36, 3, 0, 24, 16, 16, 3, 18, 8, 24]
2447 ni = [0.600073641753024, -9.36203654849857, 24.6590798594147, -107.014222858224, -91582131580576.8, -8623.32011700662, -23.5837344740032, 2.52304969384128E+17, -3.89718771997719E+18, -3.33775713645296E+22, 35649946963.6328, -1.48547544720641E+26, 3.30611514838798E+18, 8.13641294467829E+37]
2448 hs = h/2600.0
2449 ps = 0.0
2450 for i in range(0,14):
2451 ps = ps + ni[i] * (pow((hs - 1.02), Ii[i])) * (pow((hs - 0.608), Ji[i]))
2452 return ps * 22.0
2453
2454
2455
2456 def p3sat_s(s):
2457 Ii = [0, 1, 1, 4, 12, 12, 16, 24, 28, 32]
2458 Ji = [0, 1, 32, 7, 4, 14, 36, 10, 0, 18]
2459 ni = [0.639767553612785, -12.9727445396014, -2.24595125848403E+15, 1774667.41801846, 7170793495.71538, -3.78829107169011E+17, -9.55586736431328E+34, 1.87269814676188E+23, 119254746466.473, 1.10649277244882E+36]
2460 Sigma = s / 5.2
2461 Pi = 0.0
2462 for i in range(0,10):
2463 Pi = Pi + ni[i] * (Sigma - 1.03) ** Ii[i] * (Sigma - 0.699) ** Ji[i]
2464 return Pi * 22.0
2465
2466
2467 #***********************************************************************************************************
2468 #*4.3 Region boundary 1to3 & 3to2 as a functions of s
2469 #***********************************************************************************************************
2470 def hB13_s(s):
2471 #%Supplementary Release on Backward Equations ( ) , p h s for Region 3,
2472 #%Chapter 4.5 page 23.
2473 Ii = [0, 1, 1, 3, 5, 6]
2474 Ji = [0, -2, 2, -12, -4, -3]
2475 ni = [0.913965547600543, -4.30944856041991E-05, 60.3235694765419, 1.17518273082168E-18, 0.220000904781292, -69.0815545851641]
2476 Sigma = s / 3.8
2477 eta = 0.0
2478 for i in range(0,6):
2479 eta = eta + ni[i] * (Sigma - 0.884) ** Ii[i] * (Sigma - 0.864) ** Ji[i]
2480 return eta * 1700.0
2481
2482
2483 def TB23_hs(h, s):
2484 #%Supplementary Release on Backward Equations ( ) , p h s for Region 3,
2485 #%Chapter 4.6 page 25.
2486 Ii = [-12, -10, -8, -4, -3, -2, -2, -2, -2, 0, 1, 1, 1, 3, 3, 5, 6, 6, 8, 8, 8, 12, 12, 14, 14]
2487 Ji = [10, 8, 3, 4, 3, -6, 2, 3, 4, 0, -3, -2, 10, -2, -1, -5, -6, -3, -8, -2, -1, -12, -1, -12, 1]
2488 ni = [6.2909626082981E-04, -8.23453502583165E-04, 5.15446951519474E-08, -1.17565945784945, 3.48519684726192, -5.07837382408313E-12, -2.84637670005479, -2.36092263939673, 6.01492324973779, 1.48039650824546, 3.60075182221907E-04, -1.26700045009952E-02, -1221843.32521413, 0.149276502463272, 0.698733471798484, -2.52207040114321E-02, 1.47151930985213E-02, -1.08618917681849, -9.36875039816322E-04, 81.9877897570217, -182.041861521835, 2.61907376402688E-06, -29162.6417025961, 1.40660774926165E-05, 7832370.62349385]
2489 Sigma = s / 5.3
2490 eta = h / 3000.0
2491 teta = 0.0
2492 for i in range(0,25):
2493 teta = teta + ni[i] * (eta - 0.727) ** Ii[i] * (Sigma - 0.864) ** Ji[i]
2494 return teta * 900.0
2495
2496
2497 '''
2498 %***********************************************************************************************************
2499 %*5 Transport properties
2500 %***********************************************************************************************************
2501 '''
2502
2503 #***********************************************************************************************************
2504 #*5.1 Viscosity (IAPWS formulation 1985, Revised 2003)
2505 #***********************************************************************************************************
2506 def my_AllRegions_pT(p,T):
2507 h0 = [0.5132047, 0.3205656, 0, 0, -0.7782567, 0.1885447]
2508 h1 = [0.2151778, 0.7317883, 1.241044, 1.476783, 0, 0]
2509 h2 = [-0.2818107, -1.070786, -1.263184, 0, 0, 0]
2510 h3 = [0.1778064, 0.460504, 0.2340379, -0.4924179, 0, 0]
2511 h4 = [-0.0417661, 0, 0, 0.1600435, 0, 0]
2512 h5 = [0, -0.01578386, 0, 0, 0, 0]
2513 h6 = [0, 0, 0, -0.003629481, 0, 0]
2514
2515 #%Calcualte density.
2516 if region_pT(p,T) == 1:
2517 rho = 1.0/ v1_pT(p, T)
2518 elif region_pT(p,T) == 2:
2519 rho = 1.0 / v2_pT(p, T)
2520 elif region_pT(p,T) == 3:
2521 hs = h3_pT(p, T)
2522 rho = 1.0 / v3_ph(p, hs)
2523 elif region_pT(p,T) == 4:
2524 rho = float('NaN')
2525 elif region_pT(p,T) == 5:
2526 rho = 1.0 / v5_pT(p, T)
2527 else :
2528 return float('NaN')
2529
2530 rhos = rho / 317.763
2531 Ts = T / 647.226
2532 ps = p / 22.115
2533
2534 #%Check valid area
2535 if T > 900.0 + 273.15 or (T > 600.0 + 273.15 and p > 300.0) or (T > 150.0 + 273.15 and p > 350.0) or p > 500.0 :
2536 return float('NaN')
2537 my0 = Ts ** 0.5 / (1 + 0.978197 / Ts + 0.579829 / (Ts ** 2) - 0.202354 / (Ts ** 3))
2538 Sum = 0.0
2539 for i in range(0,6): #in Matlab for i = 0 : 5
2540 Sum = Sum + h0(i+1) * (1 / Ts - 1) ** i + h1(i+1) * (1 / Ts - 1) ** i * (rhos - 1) ** 1 + h2(i+1) * (1 / Ts - 1) ** i * (rhos - 1) ** 2 + h3(i+1) * (1 / Ts - 1) ** i * (rhos - 1) ** 3 + h4(i+1) * (1 / Ts - 1) ** i * (rhos - 1) ** 4 + h5(i+1) * (1 / Ts - 1) ** i * (rhos - 1) ** 5 + h6(i+1) * (1 / Ts - 1) ** i * (rhos - 1) ** 6.0
2541 my1 = math.exp(rhos * Sum)
2542 mys = my0 * my1
2543 return mys * 0.000055071
2544
2545
2546 def my_AllRegions_ph(p,h):
2547 h0 = [0.5132047, 0.3205656, 0, 0, -0.7782567, 0.1885447]
2548 h1 = [0.2151778, 0.7317883, 1.241044, 1.476783, 0, 0]
2549 h2 = [-0.2818107, -1.070786, -1.263184, 0, 0, 0]
2550 h3 = [0.1778064, 0.460504, 0.2340379, -0.4924179, 0, 0]
2551 h4 = [-0.0417661, 0, 0, 0.1600435, 0, 0]
2552 h5 = [0, -0.01578386, 0, 0, 0, 0]
2553 h6 = [0, 0, 0, -0.003629481, 0, 0]
2554
2555 #%Calcualte density.
2556 if region_ph(p,h) == 1:
2557 Ts = T1_ph(p, h)
2558 T = Ts
2559 rho = 1.0 / v1_pT(p, Ts)
2560 elif region_ph(p,h) == 2:
2561 Ts = T2_ph(p, h)
2562 T = Ts
2563 rho = 1.0 / v2_pT(p, Ts)
2564 elif region_ph(p,h) == 3:
2565 rho = 1.0 / v3_ph(p, h)
2566 T = T3_ph(p, h)
2567 elif region_ph(p,h) == 4:
2568 xs = x4_ph(p, h)
2569 if p < 16.529 :
2570 v4v = v2_pT(p, T4_p(p))
2571 v4L = v1_pT(p, T4_p(p))
2572 else :
2573 v4v = v3_ph(p, h4V_p(p))
2574 v4L = v3_ph(p, h4L_p(p))
2575
2576 rho = 1.0 / (xs * v4v + (1 - xs) * v4L)
2577 T = T4_p(p)
2578 elif region_ph(p,h) == 5:
2579 Ts = T5_ph(p, h)
2580 T = Ts
2581 rho = 1.0 / v5_pT(p, Ts)
2582 else:
2583 return float('NaN')
2584 rhos = rho / 317.763
2585 Ts = T / 647.226
2586 ps = p / 22.115
2587 #%Check valid area
2588 if T > 900 + 273.15 or (T > 600 + 273.15 and p > 300) or (T > 150 + 273.15 and p > 350) or p > 500 :
2589 return float('NaN')
2590 my0 = Ts ** 0.5 / (1 + 0.978197 / Ts + 0.579829 / (Ts ** 2) - 0.202354 / (Ts ** 3))
2591
2592 Sum = 0.0
2593 for i in range(0,6):
2594 Sum = Sum + h0(i+1) * (1 / Ts - 1) ** i + h1(i+1) * (1 / Ts - 1) ** i * (rhos - 1) ** 1 + h2(i+1) * (1 / Ts - 1) ** i * (rhos - 1) ** 2 + h3(i+1) * (1 / Ts - 1) ** i * (rhos - 1) ** 3 + h4(i+1) * (1 / Ts - 1) ** i * (rhos - 1) ** 4 + h5(i+1) * (1 / Ts - 1) ** i * (rhos - 1) ** 5 + h6(i+1) * (1 / Ts - 1) ** i * (rhos - 1) ** 6.0
2595 my1 = math.exp(rhos * Sum)
2596 mys = my0 * my1
2597 return mys * 0.000055071
2598
2599
2600 #***********************************************************************************************************
2601 #*5.2 Thermal Conductivity (IAPWS formulation 1985)
2602 #***********************************************************************************************************
2603 def tc_ptrho( p, T, rho):
2604 #%Revised release on the IAPWS formulation 1985 for the Thermal Conductivity of ordinary water
2605 #%IAPWS September 1998
2606 #%Page 8
2607 #%ver2.6 Start corrected bug
2608 if T < 273.15:
2609 return float('NaN') # %Out of range of validity (para. B4)
2610 elif T < 500 + 273.15:
2611 if p > 100:
2612 return float('NaN') # %Out of range of validity (para. B4)
2613
2614 elif T <= 650 + 273.15:
2615 if p > 70:
2616 return float('NaN') #%Out of range of validity (para. B4)
2617
2618 elif T <= 800 + 273.15:
2619 if p > 40:
2620 return float('NaN') # %Out of range of validity (para. B4)
2621
2622 #%ver2.6 End corrected bug
2623 T = T / 647.26
2624 rho = rho / 317.7
2625 tc0 = T ** 0.5 * (0.0102811 + 0.0299621 * T + 0.0156146 * T ** 2 - 0.00422464 * T ** 3)
2626 tc1 = -0.39707 + 0.400302 * rho + 1.06 * math.exp(-0.171587 * (rho + 2.39219) ** 2)
2627 dT = abs(T - 1) + 0.00308976
2628 Q = 2 + 0.0822994 / dT ** (3 / 5)
2629 if T >= 1 :
2630 s = 1.0 / dT
2631 else:
2632 s = 10.0932 / dT ** (3 / 5)
2633 tc2 = (0.0701309 / T ** 10 + 0.011852) * rho ** (9 / 5) * math.exp(0.642857 * (1 - rho ** (14 / 5))) + 0.00169937 * s * rho ** Q * math.exp((Q / (1 + Q)) * (1 - rho ** (1 + Q))) - 1.02 * math.exp(-4.11717 * T ** (3 / 2) - 6.17937 / rho ** 5)
2634 return tc0 + tc1 + tc2
2635
2636
2637
2638
2639 #***********************************************************************************************************
2640 #*5.3 Surface Tension
2641 #***********************************************************************************************************
2642 def Surface_Tension_T(T):
2643 #%IAPWS Release on Surface Tension of Ordinary Water Substance,
2644 #%September 1994
2645 tc = 647.096 # %K
2646 B = 0.2358 # %N/m
2647 bb = -0.625
2648 my = 1.256
2649 if T < 0.01 or T > tc :
2650 return float('NaN') # %"Out of valid region"
2651 tau = 1 - T / tc
2652 return B * tau ** my * (1.0 + bb * tau)
2653
2654 '''
2655 %***********************************************************************************************************
2656 %*6 Units *
2657 %***********************************************************************************************************
2658 '''
2659 def toSIunit_p( Ins ):
2660 #%Translate bar to MPa
2661 return Ins / 10.0
2662 def fromSIunit_p( Ins ):
2663 #%Translate bar to MPa
2664 return Ins * 10.0
2665 def toSIunit_T( Ins ):
2666 #%Translate degC to Kelvon
2667 return Ins + 273.15
2668 def fromSIunit_T( Ins ):
2669 #%Translate Kelvin to degC
2670 return Ins - 273.15
2671 def toSIunit_h( Ins ):
2672 return Ins
2673 def fromSIunit_h( Ins ):
2674 return Ins
2675 def toSIunit_v( Ins ):
2676 return Ins
2677 def fromSIunit_v( Ins ):
2678 return Ins
2679 def toSIunit_s( Ins ):
2680 return Ins
2681 def fromSIunit_s( Ins ):
2682 return Ins
2683 def toSIunit_u( Ins ):
2684 return Ins
2685 def fromSIunit_u( Ins ):
2686 return Ins
2687 def toSIunit_Cp( Ins ):
2688 return Ins
2689 def fromSIunit_Cp( Ins ):
2690 return Ins
2691 def toSIunit_Cv( Ins ):
2692 return Ins
2693 def fromSIunit_Cv( Ins ):
2694 return Ins
2695 def toSIunit_w( Ins ):
2696 return Ins
2697 def fromSIunit_w( Ins ):
2698 return Ins
2699 def toSIunit_tc( Ins ):
2700 return Ins
2701 def fromSIunit_tc( Ins ):
2702 return Ins
2703 def toSIunit_st( Ins ):
2704 return Ins
2705 def fromSIunit_st( Ins ):
2706 return Ins
2707 def toSIunit_x( Ins ):
2708 return Ins
2709 def fromSIunit_x( Ins ):
2710 return Ins
2711 def toSIunit_vx( Ins ):
2712 return Ins
2713 def fromSIunit_vx( Ins ):
2714 return Ins
2715 def toSIunit_my( Ins ):
2716 return Ins
2717 def fromSIunit_my( Ins ):
2718 return Ins
2719
2720 '''
2721 %***********************************************************************************************************
2722 %*7 Verification *
2723 %***********************************************************************************************************
2724 '''
2725 def check():
2726 err=0.0
2727 print 'Region 1'
2728 #***********************************************************************************************************
2729 #* 7.1 Verifiy region 1
2730 #***********************************************************************************************************
2731 #%IF-97 Table 5, Page 9
2732 p=[30/10.0,800/10.0,30/10.0]
2733 T=[300.0,300.0,500.0]
2734 #Fun={'v1_pT','h1_pT','u1_pT','s1_pT','Cp1_pT','w1_pT'}
2735 Fun=[v1_pT,h1_pT,u1_pT,s1_pT,Cp1_pT,w1_pT]
2736 IF97=numpy.matrix('0.00100215168 0.000971180894 0.001202418; \
2737 115.331273 184.142828 975.542239;\
2738 112.324818 106.448356 971.934985;\
2739 0.392294792 0.368563852 2.58041912;\
2740 4.17301218 4.01008987 4.65580682;\
2741 1507.73921 1634.69054 1240.71337')
2742
2743 R1=numpy.zeros((6,3),float)
2744 for i in range(0,3):
2745 for j in range(0,6):
2746 #R1[j,i]=eval([char(Fun[j]),'(',num2str(p[i]),',',num2str(T[i]),');'])
2747 R1[j,i]=Fun[j](p[i],T[i])
2748
2749 Region1_error=abs(numpy.divide((R1-IF97),IF97))
2750 #err=err+sum(sum(Region1_error>1E-8))
2751 if numpy.sum(Region1_error)>1E-8 :
2752 err=err+numpy.sum(numpy.sum(Region1_error))
2753
2754 print err #OK 10e-8
2755
2756
2757 #%IF-97 Table 7, Page 11
2758 p=[30/10.0,800/10.0,800/10.0]
2759 h=[500.0,500.0,1500.0]
2760 IF97=[391.798509,378.108626,611.041229]
2761 R1=numpy.zeros((3,1),float)
2762 for i in range(0,3):
2763 R1[i]=T1_ph(p[i],h[i])
2764 T1_ph_error=abs(numpy.divide((numpy.transpose(R1)-IF97),IF97))
2765 if numpy.sum(T1_ph_error) > 1E-8 :
2766 err=err+numpy.sum(numpy.sum(T1_ph_error))
2767 print err
2768
2769 #%IF-97 Table 9, Page 12
2770 p=[30/10.0,800/10.0,800/10.0]
2771 s=[0.5,0.5,3.0]
2772 IF97=[307.842258,309.979785,565.899909]
2773 R1=numpy.zeros((3,1),float)
2774 for i in range(0,3):
2775 R1[i]=T1_ps(p[i],s[i])
2776 T1_ps_error=abs(numpy.divide((numpy.transpose(R1)-IF97),IF97))
2777 if numpy.sum(T1_ps_error)>1E-8:
2778 err=err+numpy.sum(numpy.sum(T1_ps_error))
2779
2780 print err
2781 #%Supplementary Release on Backward Equations
2782 #%for Pressure as a Function of Enthalpy and Entropy p(h,s)
2783 #%Table 3, Page 6
2784 h=[0.001,90.0,1500.0]
2785 s=[0,0,3.4];
2786 IF97=[0.0009800980612,91.929547272,58.68294423]
2787 R1=numpy.zeros((3,1),float)
2788 for i in range(0,3):
2789 R1[i]=p1_hs(h[i],s[i])
2790
2791 p1_hs_error=abs(numpy.divide((numpy.transpose(R1)-IF97),IF97))
2792 if numpy.sum(p1_hs_error)>1E8:
2793 err=err+numpy.sum(numpy.sum(p1_hs_error))
2794 print err
2795 #***********************************************************************************************************
2796 #* 7.2 Verifiy region 2
2797 #***********************************************************************************************************
2798 print 'Region 2'
2799 #% IF-97 Table 15, Page 17
2800 p=[0.035/10.0,0.035/10.0,300/10.0]
2801 T=[300.0,700.0,700.0]
2802 Fun=[v2_pT,h2_pT,u2_pT,s2_pT,Cp2_pT,w2_pT]
2803 IF97=numpy.matrix('39.4913866 92.3015898 0.00542946619;\
2804 2549.91145 3335.68375 2631.49474;\
2805 2411.6916 3012.62819 2468.61076;\
2806 8.52238967 10.1749996 5.17540298;\
2807 1.91300162 2.08141274 10.3505092;\
2808 427.920172 644.289068 480.386523')
2809 R2=numpy.zeros((6,3),float)
2810 for i in range(0,3):
2811 for j in range (0,6):
2812 #R2(j,i)=eval([char(Fun(j)),'(',num2str(p[i]),',',num2str(T[i]),');']);
2813 R2[j,i]=Fun[j](p[i],T[i])
2814
2815 Region2_error=abs(numpy.divide((R2-IF97),IF97))
2816 if numpy.sum(Region2_error)>1E-8:
2817 err=err+numpy.sum(numpy.sum(Region2_error))
2818 print err
2819 #%IF-97 Table 24, Page 25
2820 p=[0.01/10.0,30/10.0,30/10.0,50/10.0,50/10.0,250/10.0,400/10.0,600/10.0,600/10.0]
2821 h=[3000.0,3000.0,4000.0,3500.0,4000.0,3500.0,2700.0,2700.0,3200.0]
2822 IF97=[534.433241,575.37337,1010.77577,801.299102,1015.31583,875.279054,743.056411,791.137067,882.75686]
2823 R2=numpy.zeros((9,1),float)
2824 for i in range(0,9):
2825 R2[i]=T2_ph(p[i],h[i])
2826 T2_ph_error=abs(numpy.divide((numpy.transpose(R2)-IF97),IF97))
2827 if numpy.sum(T2_ph_error) > 1E-8:
2828 err=err+numpy.sum(numpy.sum(T2_ph_error))
2829 print err
2830 #%IF-97 Table 29, Page 29
2831 p=[1/10.0,1/10.0,25/10.0,80/10.0,80/10.0,900/10.0,200/10.0,800/10.0,800/10.0]
2832 s=[7.5,8.0,8.0,6.0,7.5,6.0,5.75,5.25,5.75]
2833 IF97=[399.517097,514.127081,1039.84917,600.48404,1064.95556,1038.01126,697.992849,854.011484,949.017998]
2834 R2=numpy.zeros((9,1),float)
2835 for i in range(0,9):
2836 R2[i]=T2_ps(p[i],s[i])
2837 T2_ps_error=abs(numpy.divide((numpy.transpose(R2)-IF97),IF97))
2838 if numpy.sum(T2_ps_error) >1E-8:
2839 err=err+numpy.sum(numpy.sum(T2_ps_error))
2840 print err
2841 #%Supplementary Release on Backward Equations for Pressure as a Function of Enthalpy and Entropy p(h,s)
2842 #%Table 3, Page 6
2843 h=[2800.0,2800.0,4100.0,2800.0,3600.0,3600.0,2800.0,2800,3400.0]
2844 s=[6.5,9.5,9.5,6,6,7,5.1,5.8,5.8]
2845 IF97=[1.371012767,0.001879743844,0.1024788997,4.793911442,83.95519209,7.527161441,94.3920206,8.414574124,83.76903879]
2846 R2=numpy.zeros((9,1),float)
2847 for i in range(0,9):
2848 R2[i]=p2_hs(h[i],s[i])
2849 p2_hs_error=abs(numpy.divide((numpy.transpose(R2)-IF97),IF97))
2850 if numpy.sum(p2_hs_error)>1E-8:
2851 err=err+numpy.sum(numpy.sum(p2_hs_error))
2852
2853 print err
2854 #***********************************************************************************************************
2855 #* 7.3 Verifiy region 3
2856 #***********************************************************************************************************
2857 print 'Region 3'
2858 #% IF-97 Table 33, Page 32
2859 T=[650.0,650.0,750.0]
2860 rho=[500.0,200.0,500.0]
2861 Fun=[p3_rhoT,h3_rhoT,u3_rhoT,s3_rhoT,Cp3_rhoT,w3_rhoT]
2862 IF97=numpy.matrix('25.5837018 22.2930643 78.3095639;\
2863 1863.43019 2375.12401 2258.68845;\
2864 1812.26279 2263.65868 2102.06932;\
2865 4.05427273 4.85438792 4.46971906;\
2866 13.8935717 44.6579342 6.34165359;\
2867 502.005554 383.444594 760.696041')
2868 R3=numpy.zeros((6,3),float)
2869 for i in range(0,3):
2870 for j in range(0,6):
2871 #R3(j,i)=eval([char(Fun(j)),'(',num2str(rho[i]),',',num2str(T[i]),');']);
2872 R3[j,i]=Fun[j](rho[i],T[i])
2873 Region3_error=abs(numpy.divide((R3-IF97),IF97))
2874 if numpy.sum(Region3_error)>1E-8:
2875 err=err+numpy.sum(numpy.sum(Region3_error))
2876 print err
2877 #%T3_ph
2878 p=[200/10.0,500/10.0,1000/10.0,200/10.0,500/10.0,1000/10.0]
2879 h=[1700.0,2000.0,2100.0,2500.0,2400.0,2700.0]
2880 IF97=[629.3083892,690.5718338,733.6163014,641.8418053,735.1848618,842.0460876]
2881 R3=numpy.zeros((6,1),float)
2882 for i in range(0,6):
2883 R3[i]=T3_ph(p[i],h[i])
2884 T3_ph_error=abs(numpy.divide((numpy.transpose(R3)-IF97),IF97))
2885 if numpy.sum(T3_ph_error)>1E-8 :
2886 err=err+numpy.sum(numpy.sum(T3_ph_error))
2887 print err
2888 #%v3_ph
2889 p=[200/10.0,500/10.0,1000/10.0,200/10.0,500/10.0,1000/10.0]
2890 h=[1700.0,2000.0,2100.0,2500.0,2400.0,2700.0]
2891 IF97=[0.001749903962,0.001908139035,0.001676229776,0.006670547043,0.0028012445,0.002404234998]
2892 R3=numpy.zeros((6,1),float)
2893 for i in range(0,6):
2894 R3[i]=v3_ph(p[i],h[i])
2895 v3_ph_error=abs(numpy.divide((numpy.transpose(R3)-IF97),IF97))
2896 if numpy.sum(v3_ph_error)>1E-7 :
2897 err=err+numpy.sum(numpy.sum(v3_ph_error))
2898 print err
2899 #%T3_ps
2900 p=[200/10.0,500/10.0,1000/10.0,200/10.0,500/10.0,1000/10.0]
2901 s=[3.7,3.5,4,5,4.5,5]
2902 IF97=[620.8841563,618.1549029,705.6880237,640.1176443,716.3687517,847.4332825]
2903 R3=numpy.zeros((6,1),float)
2904 for i in range(0,6):
2905 R3[i]=T3_ps(p[i],s[i])
2906 T3_ps_error=abs(numpy.divide((numpy.transpose(R3)-IF97),IF97))
2907 if numpy.sum(T3_ps_error) >1E-8:
2908 err=err+numpy.sum(numpy.sum(T3_ps_error))
2909 print err
2910 #%v3_ps
2911 p=[200/10.0,500/10.0,1000/10.0,200/10.0,500/10.0,1000/10.0]
2912 s=[3.7,3.5,4,5,4.5,5]
2913 IF97=[0.001639890984,0.001423030205,0.001555893131,0.006262101987,0.002332634294,0.002449610757]
2914 R3=numpy.zeros((6,1),float)
2915 for i in range(0,6):
2916 R3[i]=v3_ps(p[i],s[i])
2917 v3_ps_error=abs(numpy.divide((numpy.transpose(R3)-IF97),IF97))
2918 if numpy.sum(v3_ps_error) >1E-8:
2919 err=err+numpy.sum(numpy.sum(v3_ps_error))
2920 print err
2921 #%p3_hs
2922 h=[1700,2000,2100,2500,2400,2700]
2923 s=[3.8,4.2,4.3,5.1,4.7,5]
2924 IF97=[25.55703246,45.40873468,60.7812334,17.20612413,63.63924887,88.39043281]
2925 R3=numpy.zeros((6,1),float)
2926 for i in range(0,6):
2927 R3[i]=p3_hs(h[i],s[i])
2928 p3_hs_error=abs(numpy.divide((numpy.transpose(R3)-IF97),IF97))
2929 if numpy.sum(p3_hs_error)>1E-8:
2930 err=err+numpy.sum(numpy.sum(p3_hs_error))
2931 print err
2932 #%h3_pT (Iteration)
2933 p=[255.83702/10.0,222.93064/10.0,783.09564/10.0]
2934 T=[650.0,650.0,750.0]
2935 IF97=[1863.271389,2375.696155,2258.626582]
2936 R3=numpy.zeros((3,1),float)
2937 for i in range(0,3):
2938 R3[i]=h3_pT(p[i],T[i])
2939 h3_pT_error=abs(numpy.divide((numpy.transpose(R3)-IF97),IF97))
2940 if numpy.sum(h3_pT_error)>1E-6:
2941 err=err+numpy.sum(numpy.sum(h3_pT_error>1E-6)) #%Decimals in IF97
2942 print err
2943
2944 #***********************************************************************************************************
2945 #* 7.4 Verifiy region 4
2946 #***********************************************************************************************************
2947 print 'Region 4'
2948 #%Saturation pressure, If97, Table 35, Page 34
2949 T=[300.0,500.0,600.0]
2950 IF97=[0.0353658941/10.0, 26.3889776/10.0, 123.443146/10.0]
2951 R4=numpy.zeros((3,1),float)
2952 for i in range(0,3):
2953 R4[i]=p4_T(T[i])
2954 p4_t_error=abs(numpy.divide((numpy.transpose(R4)-IF97),IF97))
2955 if numpy.sum(p4_t_error)>1E-7:
2956 err=err+numpy.sum(numpy.sum( p4_t_error))
2957 print err
2958
2959 T=[1.0/10.0,10.0/10.0,100.0/10.0]
2960 IF97=[372.755919,453.035632,584.149488]
2961 R4=numpy.zeros((3,1),float)
2962 for i in range(0,3):
2963 R4[i]=T4_p(T[i])
2964 T4_p_error=abs(numpy.divide((numpy.transpose(R4)-IF97),IF97))
2965 if numpy.sum(T4_p_error)>1E-7:
2966 err=err+numpy.sum(numpy.sum( T4_p_error))
2967 print err
2968 s=[1.0,2.0,3.0,3.8,4.0,4.2,7.0,8.0,9.0,5.5,5,4.5]
2969 IF97=[308.5509647,700.6304472,1198.359754,1685.025565,1816.891476,1949.352563,2723.729985,2599.04721,2511.861477,2687.69385,2451.623609,2144.360448]
2970 R4=numpy.zeros((12,1),float)
2971 for i in range(0,12):
2972 R4[i]=h4_s(s[i])
2973 h4_s_error=abs(numpy.divide((numpy.transpose(R4)-IF97),IF97))
2974 if numpy.sum( h4_s_error)>1E-7:
2975 err=err+numpy.sum(numpy.sum( h4_s_error))
2976 print err
2977 #***********************************************************************************************************
2978 #* 7.5 Verifiy region 5
2979 #***********************************************************************************************************
2980 print 'Region 5'
2981 #%* 7.5 Verifiy region 5
2982 #% IF-97 Table 42, Page 39
2983 T=[1500.0,1500.0,2000.0]
2984 p=[5.0/10.0,80.0/10.0,80.0/10.0]
2985 Fun=[v5_pT,h5_pT,u5_pT,s5_pT,Cp5_pT,w5_pT]
2986 IF97=numpy.matrix('1.38455354 0.0865156616 0.115743146;\
2987 5219.76332 5206.09634 6583.80291;\
2988 4527.48654 4513.97105 5657.85774;\
2989 9.65408431 8.36546724 9.15671044;\
2990 2.61610228 2.64453866 2.8530675;\
2991 917.071933 919.708859 1054.35806')
2992 R5=numpy.zeros((6,3),float)
2993 for i in range(0,3):
2994 for j in range(0,6):
2995 #R5(j,i)=eval([char(Fun(j)),'(',num2str(p[i]),',',num2str(T[i]),');']);
2996 R5[j,i]= Fun[j](p[i],T[i])
2997 Region5_error=abs(numpy.divide((R5-IF97),IF97))
2998 if numpy.sum(Region5_error)>1E-8:
2999 err=err+numpy.sum(numpy.sum(Region5_error))
3000 print err
3001 #%T5_ph (Iteration)
3002 p=[0.5,8.0,8.0]
3003 h=[5219.76331549428,5206.09634477373,6583.80290533381]
3004 IF97=[1500.0,1500.0,2000.0]
3005 R5=numpy.zeros((3,1),float)
3006 for i in range(0,3):
3007 R5[i]=T5_ph(p[i],h[i])
3008 T5_ph_error=abs(numpy.divide((numpy.transpose(R5)-IF97),IF97))
3009 if numpy.sum(T5_ph_error) >1E-7:
3010 err=err+numpy.sum(numpy.sum(T5_ph_error>1E-7)) #%Decimals in IF97
3011 print err
3012 #%T5_ps (Iteration)
3013 p=[0.5,8.0,8.0]
3014 s=[9.65408430982588,8.36546724495503,9.15671044273249]
3015 IF97=[1500.0,1500.0,2000.0]
3016 R5=numpy.zeros((3,1),float)
3017 for i in range(0,3):
3018 R5[i]=T5_ps(p[i],s[i])
3019 T5_ps_error=abs(numpy.divide((numpy.transpose(R5)-IF97),IF97))
3020 if numpy.sum(T5_ps_error)>1E-4:
3021 err=err+numpy.sum(numpy.sum(T5_ps_error)) #%Decimals in IF97
3022 print err
memoire cycle