Skip to content
This repository was archived by the owner on Feb 20, 2021. It is now read-only.

Commit 15449e3

Browse files
raldridge11raldridge11
committed
Closes #148, adds my_allregions_pt and tests.
1 parent 1952025 commit 15449e3

File tree

3 files changed

+171
-110
lines changed

3 files changed

+171
-110
lines changed

‎XSteamPython/Viscosity.py

+126
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,126 @@
1+
# -*- coding: utf-8 -*-
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+
import math
11+
12+
import numpy as np
13+
14+
import Constants
15+
import Regions
16+
import Region1
17+
import Region2
18+
import Region3
19+
import Region5
20+
21+
h0 = np.array([0.5132047, 0.3205656, 0.0, 0.0, -0.7782567, 0.1885447])
22+
h1 = np.array([0.2151778, 0.7317883, 1.241044, 1.476783, 0.0, 0.0])
23+
h2 = np.array([-0.2818107, -1.070786, -1.263184, 0.0, 0.0, 0.0])
24+
h3 = np.array([0.1778064, 0.460504, 0.2340379, -0.4924179, 0.0, 0.0])
25+
h4 = np.array([-0.0417661, 0.0, 0.0, 0.1600435, 0.0, 0.0])
26+
h5 = np.array([0.0, -0.01578386, 0.0, 0.0, 0.0, 0.0])
27+
h6 = np.array([0.0, 0.0, 0.0, -0.003629481, 0.0, 0.0])
28+
29+
def my_allregions_pT(pressure, temperature):
30+
'''Viscosity (IAPWS formulation 1985, Revised 2003)'''
31+
32+
# Check valid area
33+
if temperature > 900.0 + 273.15 or \
34+
(temperature > 600.0 + 273.15 and pressure > 300.0) or \
35+
(temperature > 150.0 + 273.15 and pressure > 350.0) or \
36+
pressure > 500.0:
37+
return Constants._errorValue
38+
39+
region = Regions.region_pt(pressure, temperature)
40+
if region is None: return Constants._errorValue
41+
42+
density = Constants._errorValue
43+
if region == 1:
44+
density = 1.0/Region1.v1_pt(pressure, temperature)
45+
elif region == 2:
46+
density = 1.0/Region2.v2_pt(pressure, temperature)
47+
elif region == 3:
48+
density = 1.0/Region3.v3_ph(pressure, Region3.h3_pt(pressure, temperature))
49+
elif region == 4:
50+
density = Constants._errorValue
51+
elif region == 5:
52+
density = 1.0/Region5.v5_pt(pressure, temperature)
53+
54+
rhos = density/317.63
55+
ts = temperature/647.226
56+
ps = pressure/22.115
57+
58+
my_0 = ts**0.5/(1.0 + 0.978197/ts + 0.579829/(ts**2.0) - 0.202354/(ts**3.0))
59+
total = 0.0
60+
a, b = 1.0/ts - 1.0, rhos - 1.0
61+
for i in range(6):
62+
total += h0[i]*a**i + h1[i]*a**i*b + h2[i]*a**i*b**2 + h3[i]*a**i*b**3 + h4[i]*a**i*b**4 + h5[i]*a**i*b**5 + h6[i]*a**i*b**6
63+
64+
my_1 = math.exp(rhos*total)
65+
return my_0*my_1*0.000055071
66+
67+
#Rem
68+
#Rem Private Function my_AllRegions_ph(ByVal p As Double, ByVal h As Double) As Double
69+
#Rem Dim h0, h1, h2, h3, h4, h5, h6 As Variant, rho, T, Ts, ps, my0, sum, my1, rhos, v4V, v4L, xs As Double, i As Integer
70+
#Rem h0 = Array(0.5132047, 0.3205656, 0, 0, -0.7782567, 0.1885447)
71+
#Rem h1 = Array(0.2151778, 0.7317883, 1.241044, 1.476783, 0, 0)
72+
#Rem h2 = Array(-0.2818107, -1.070786, -1.263184, 0, 0, 0)
73+
#Rem h3 = Array(0.1778064, 0.460504, 0.2340379, -0.4924179, 0, 0)
74+
#Rem h4 = Array(-0.0417661, 0, 0, 0.1600435, 0, 0)
75+
#Rem h5 = Array(0, -0.01578386, 0, 0, 0, 0)
76+
#Rem h6 = Array(0, 0, 0, -0.003629481, 0, 0)
77+
#Rem
78+
#Rem 'Calcualte density.
79+
#Rem Select Case region_ph(p, h)
80+
#Rem Case 1
81+
#Rem Ts = T1_ph(p, h)
82+
#Rem T = Ts
83+
#Rem rho = 1 / v1_pT(p, Ts)
84+
#Rem Case 2
85+
#Rem Ts = T2_ph(p, h)
86+
#Rem T = Ts
87+
#Rem rho = 1 / v2_pT(p, Ts)
88+
#Rem Case 3
89+
#Rem rho = 1 / v3_ph(p, h)
90+
#Rem T = T3_ph(p, h)
91+
#Rem Case 4
92+
#Rem xs = x4_ph(p, h)
93+
#Rem If p < 16.529 Then
94+
#Rem v4V = v2_pT(p, T4_p(p))
95+
#Rem v4L = v1_pT(p, T4_p(p))
96+
#Rem Else
97+
#Rem v4V = v3_ph(p, h4V_p(p))
98+
#Rem v4L = v3_ph(p, h4L_p(p))
99+
#Rem End If
100+
#Rem rho = 1 / (xs * v4V + (1 - xs) * v4L)
101+
#Rem T = T4_p(p)
102+
#Rem Case 5
103+
#Rem Ts = T5_ph(p, h)
104+
#Rem T = Ts
105+
#Rem rho = 1 / v5_pT(p, Ts)
106+
#Rem Case Else
107+
#Rem my_AllRegions_ph = CVErr(xlErrValue)
108+
#Rem Exit Function
109+
#Rem End Select
110+
#Rem rhos = rho / 317.763
111+
#Rem Ts = T / 647.226
112+
#Rem ps = p / 22.115
113+
#Rem 'Check valid area
114+
#Rem If T > 900 + 273.15 Or (T > 600 + 273.15 And p > 300) Or (T > 150 + 273.15 And p > 350) Or p > 500 Then
115+
#Rem my_AllRegions_ph = CVErr(xlErrValue)
116+
#Rem Exit Function
117+
#Rem End If
118+
#Rem my0 = Ts ^ 0.5 / (1 + 0.978197 / Ts + 0.579829 / (Ts ^ 2) - 0.202354 / (Ts ^ 3))
119+
#Rem
120+
#Rem sum = 0
121+
#Rem For i = 0 To 5
122+
#Rem sum = sum + h0(i) * (1 / Ts - 1) ^ i + h1(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 1 + h2(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 2 + h3(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 3 + h4(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 4 + h5(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 5 + h6(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 6
123+
#Rem Next i
124+
#Rem my1 = Exp(rhos * sum)
125+
#Rem my_AllRegions_ph = my0 * my1 * 0.000055071
126+
#Rem End Function

‎XSteamPython/XSteamPython.py

+2-110
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,8 @@
99
'''
1010
import math
1111

12+
import numpy as np
13+
1214
import Constants
1315
import Convert
1416
import Region1
@@ -1441,116 +1443,6 @@ def vx_ps(pressure, entropy):
14411443
else:
14421444
return Constants._errorValue
14431445

1444-
#'***********************************************************************************************************
1445-
#Rem '*5 Transport properties
1446-
#Rem '***********************************************************************************************************
1447-
#Rem '*5.1 Viscosity (IAPWS formulation 1985, Revised 2003)
1448-
#Rem '***********************************************************************************************************
1449-
#Rem Private Function my_AllRegions_pT(ByVal p As Double, ByVal T As Double) As Double
1450-
#Rem Dim h0, h1, h2, h3, h4, h5, h6 As Variant, rho, Ts, ps, my0, sum, my1, rhos As Double, i As Integer
1451-
#Rem h0 = Array(0.5132047, 0.3205656, 0, 0, -0.7782567, 0.1885447)
1452-
#Rem h1 = Array(0.2151778, 0.7317883, 1.241044, 1.476783, 0, 0)
1453-
#Rem h2 = Array(-0.2818107, -1.070786, -1.263184, 0, 0, 0)
1454-
#Rem h3 = Array(0.1778064, 0.460504, 0.2340379, -0.4924179, 0, 0)
1455-
#Rem h4 = Array(-0.0417661, 0, 0, 0.1600435, 0, 0)
1456-
#Rem h5 = Array(0, -0.01578386, 0, 0, 0, 0)
1457-
#Rem h6 = Array(0, 0, 0, -0.003629481, 0, 0)
1458-
#Rem
1459-
#Rem 'Calcualte density.
1460-
#Rem Select Case region_pT(p, T)
1461-
#Rem Case 1
1462-
#Rem rho = 1 / v1_pT(p, T)
1463-
#Rem Case 2
1464-
#Rem rho = 1 / v2_pT(p, T)
1465-
#Rem Case 3
1466-
#Rem rho = 1 / v3_ph(p, h3_pT(p, T))
1467-
#Rem Case 4
1468-
#Rem rho = CVErr(xlErrValue)
1469-
#Rem Case 5
1470-
#Rem rho = 1 / v5_pT(p, T)
1471-
#Rem Case Else
1472-
#Rem my_AllRegions_pT = CVErr(xlErrValue)
1473-
#Rem Exit Function
1474-
#Rem End Select
1475-
#Rem
1476-
#Rem rhos = rho / 317.763
1477-
#Rem Ts = T / 647.226
1478-
#Rem ps = p / 22.115
1479-
#Rem
1480-
#Rem 'Check valid area
1481-
#Rem If T > 900 + 273.15 Or (T > 600 + 273.15 And p > 300) Or (T > 150 + 273.15 And p > 350) Or p > 500 Then
1482-
#Rem my_AllRegions_pT = CVErr(xlErrValue)
1483-
#Rem Exit Function
1484-
#Rem End If
1485-
#Rem my0 = Ts ^ 0.5 / (1 + 0.978197 / Ts + 0.579829 / (Ts ^ 2) - 0.202354 / (Ts ^ 3))
1486-
#Rem sum = 0
1487-
#Rem For i = 0 To 5
1488-
#Rem sum = sum + h0(i) * (1 / Ts - 1) ^ i + h1(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 1 + h2(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 2 + h3(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 3 + h4(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 4 + h5(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 5 + h6(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 6
1489-
#Rem Next i
1490-
#Rem my1 = Exp(rhos * sum)
1491-
#Rem my_AllRegions_pT = my0 * my1 * 0.000055071
1492-
#Rem End Function
1493-
#Rem
1494-
#Rem Private Function my_AllRegions_ph(ByVal p As Double, ByVal h As Double) As Double
1495-
#Rem Dim h0, h1, h2, h3, h4, h5, h6 As Variant, rho, T, Ts, ps, my0, sum, my1, rhos, v4V, v4L, xs As Double, i As Integer
1496-
#Rem h0 = Array(0.5132047, 0.3205656, 0, 0, -0.7782567, 0.1885447)
1497-
#Rem h1 = Array(0.2151778, 0.7317883, 1.241044, 1.476783, 0, 0)
1498-
#Rem h2 = Array(-0.2818107, -1.070786, -1.263184, 0, 0, 0)
1499-
#Rem h3 = Array(0.1778064, 0.460504, 0.2340379, -0.4924179, 0, 0)
1500-
#Rem h4 = Array(-0.0417661, 0, 0, 0.1600435, 0, 0)
1501-
#Rem h5 = Array(0, -0.01578386, 0, 0, 0, 0)
1502-
#Rem h6 = Array(0, 0, 0, -0.003629481, 0, 0)
1503-
#Rem
1504-
#Rem 'Calcualte density.
1505-
#Rem Select Case region_ph(p, h)
1506-
#Rem Case 1
1507-
#Rem Ts = T1_ph(p, h)
1508-
#Rem T = Ts
1509-
#Rem rho = 1 / v1_pT(p, Ts)
1510-
#Rem Case 2
1511-
#Rem Ts = T2_ph(p, h)
1512-
#Rem T = Ts
1513-
#Rem rho = 1 / v2_pT(p, Ts)
1514-
#Rem Case 3
1515-
#Rem rho = 1 / v3_ph(p, h)
1516-
#Rem T = T3_ph(p, h)
1517-
#Rem Case 4
1518-
#Rem xs = x4_ph(p, h)
1519-
#Rem If p < 16.529 Then
1520-
#Rem v4V = v2_pT(p, T4_p(p))
1521-
#Rem v4L = v1_pT(p, T4_p(p))
1522-
#Rem Else
1523-
#Rem v4V = v3_ph(p, h4V_p(p))
1524-
#Rem v4L = v3_ph(p, h4L_p(p))
1525-
#Rem End If
1526-
#Rem rho = 1 / (xs * v4V + (1 - xs) * v4L)
1527-
#Rem T = T4_p(p)
1528-
#Rem Case 5
1529-
#Rem Ts = T5_ph(p, h)
1530-
#Rem T = Ts
1531-
#Rem rho = 1 / v5_pT(p, Ts)
1532-
#Rem Case Else
1533-
#Rem my_AllRegions_ph = CVErr(xlErrValue)
1534-
#Rem Exit Function
1535-
#Rem End Select
1536-
#Rem rhos = rho / 317.763
1537-
#Rem Ts = T / 647.226
1538-
#Rem ps = p / 22.115
1539-
#Rem 'Check valid area
1540-
#Rem If T > 900 + 273.15 Or (T > 600 + 273.15 And p > 300) Or (T > 150 + 273.15 And p > 350) Or p > 500 Then
1541-
#Rem my_AllRegions_ph = CVErr(xlErrValue)
1542-
#Rem Exit Function
1543-
#Rem End If
1544-
#Rem my0 = Ts ^ 0.5 / (1 + 0.978197 / Ts + 0.579829 / (Ts ^ 2) - 0.202354 / (Ts ^ 3))
1545-
#Rem
1546-
#Rem sum = 0
1547-
#Rem For i = 0 To 5
1548-
#Rem sum = sum + h0(i) * (1 / Ts - 1) ^ i + h1(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 1 + h2(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 2 + h3(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 3 + h4(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 4 + h5(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 5 + h6(i) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 6
1549-
#Rem Next i
1550-
#Rem my1 = Exp(rhos * sum)
1551-
#Rem my_AllRegions_ph = my0 * my1 * 0.000055071
1552-
#Rem End Function
1553-
15541446
def tc_pTrho(pressure, temperature, density):
15551447
'''Revised release on the IAPS Formulation 1985 for the Thermal Conductivity of ordinary water IAPWS September 1998 Page 8'''
15561448
if temperature < 0.0 or pressure < Constants._pressureMin \

‎tests/Viscosity_Tests.py

+43
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,43 @@
1+
# -*- coding: utf-8 -*-
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+
import os
11+
import sys
12+
import unittest
13+
14+
file_directory = os.path.dirname(__file__)
15+
src_path = os.path.join(os.path.abspath(file_directory), "..", "XSteamPython")
16+
sys.path.append(src_path)
17+
import Viscosity
18+
19+
class Test_my_AllRegions_pT(unittest.TestCase):
20+
21+
def test_my_pT_region1(self):
22+
self.assertAlmostEqual(Viscosity.my_allregions_pT(1.0, 300.0),0.0008535, places=7)
23+
24+
def test_my_pT_region2(self):
25+
self.assertAlmostEqual(Viscosity.my_allregions_pT(1.0, 500.0), 1.71e-5, places=7)
26+
27+
def test_my_pT_region3(self):
28+
self.assertAlmostEqual(Viscosity.my_allregions_pT(19.0, 640.0), 2.54e-5, places=7)
29+
30+
def test_my_pT_region4(self):
31+
self.assertAlmostEqual(Viscosity.my_allregions_pT(2.63890, 500.0), 0.0, places=7)
32+
33+
def test_my_pT_region5(self):
34+
self.assertAlmostEqual(Viscosity.my_allregions_pT(9.0, 1100.0), 4.18e-5, places=7)
35+
36+
def test_my_pT_invalid_area(self):
37+
self.assertEqual(Viscosity.my_allregions_pT(9.0, 2000.0), 2015.0)
38+
self.assertEqual(Viscosity.my_allregions_pT(350.0, 900.0), 2015.0)
39+
self.assertEqual(Viscosity.my_allregions_pT(400.0, 500.0), 2015.0)
40+
self.assertEqual(Viscosity.my_allregions_pT(600.0, 300.0), 2015.0)
41+
42+
if __name__ == '__main__':
43+
unittest.main()

0 commit comments

Comments
 (0)