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

Commit 44cab5e

Browse files
raldridge11raldridge11
committed
Closes #149, adds my_allregions_ph and tests.
1 parent 15449e3 commit 44cab5e

File tree

2 files changed

+80
-67
lines changed

2 files changed

+80
-67
lines changed

‎XSteamPython/Viscosity.py

+51-61
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@
1616
import Region1
1717
import Region2
1818
import Region3
19+
import Region4
1920
import Region5
2021

2122
h0 = np.array([0.5132047, 0.3205656, 0.0, 0.0, -0.7782567, 0.1885447])
@@ -53,7 +54,6 @@ def my_allregions_pT(pressure, temperature):
5354

5455
rhos = density/317.63
5556
ts = temperature/647.226
56-
ps = pressure/22.115
5757

5858
my_0 = ts**0.5/(1.0 + 0.978197/ts + 0.579829/(ts**2.0) - 0.202354/(ts**3.0))
5959
total = 0.0
@@ -64,63 +64,53 @@ def my_allregions_pT(pressure, temperature):
6464
my_1 = math.exp(rhos*total)
6565
return my_0*my_1*0.000055071
6666

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
67+
def my_allregions_ph(pressure, enthalpy):
68+
'''Viscosity (IAPWS formulation 1985, Revised 2003)'''
69+
region = Regions.region_ph(pressure, enthalpy)
70+
if region is None: return Constants._errorValue
71+
72+
temperature = Constants._errorValue
73+
density = Constants._errorValue
74+
75+
if region == 1:
76+
temperature = Region1.t1_ph(pressure, enthalpy)
77+
density = 1.0/Region1.v1_pt(pressure, temperature)
78+
elif region == 2:
79+
temperature = Region2.t2_ph(pressure, enthalpy)
80+
density = 1.0/Region2.v2_pt(pressure, temperature)
81+
elif region == 3:
82+
temperature = Region3.t3_ph(pressure, enthalpy)
83+
density = 1.0/Region3.v3_ph(pressure, enthalpy)
84+
elif region == 4:
85+
quality = Region4.x4_ph(pressure, enthalpy)
86+
specificVolumeVapor, specificVolumeLiquid = Constants._errorValue, Constants._errorValue
87+
temperature = Region4.t4_p(pressure)
88+
if pressure < Constants._pressureSubDomain:
89+
specificVolumeVapor = Region2.v2_pt(pressure, temperature)
90+
specificVolumeLiquid = Region1.v1_pt(pressure, temperature)
91+
else:
92+
specificVolumeVapor = Region3.v3_ph(pressure, Region4.h4_p(pressure, 'vap'))
93+
specificVolumeLiquid = Region3.v3_ph(pressure, Region4.h4_p(pressure, 'liq'))
94+
density = 1.0/(quality*specificVolumeVapor + (1.0 - quality)*specificVolumeLiquid)
95+
elif region == 5:
96+
temperature = Region5.t5_ph(pressure, enthalpy)
97+
density = 1.0/Region5.v5_pt(pressure, temperature)
98+
99+
# Check valid area
100+
if temperature > 900.0 + 273.15 or \
101+
(temperature > 600.0 + 273.15 and pressure > 300.0) or \
102+
(temperature > 150.0 + 273.15 and pressure > 350.0) or \
103+
pressure > 500.0:
104+
return Constants._errorValue
105+
106+
rhos = density/317.63
107+
ts = temperature/647.226
108+
109+
my_0 = ts**0.5/(1.0 + 0.978197/ts + 0.579829/(ts**2.0) - 0.202354/(ts**3.0))
110+
total = 0.0
111+
a, b = 1.0/ts - 1.0, rhos - 1.0
112+
for i in range(6):
113+
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
114+
115+
my_1 = math.exp(rhos*total)
116+
return my_0*my_1*0.000055071

‎tests/Viscosity_Tests.py

+29-6
Original file line numberDiff line numberDiff line change
@@ -18,26 +18,49 @@
1818

1919
class Test_my_AllRegions_pT(unittest.TestCase):
2020

21-
def test_my_pT_region1(self):
21+
def test_my_allregions_pT_region1(self):
2222
self.assertAlmostEqual(Viscosity.my_allregions_pT(1.0, 300.0),0.0008535, places=7)
2323

24-
def test_my_pT_region2(self):
24+
def test_my_allregions_pT_region2(self):
2525
self.assertAlmostEqual(Viscosity.my_allregions_pT(1.0, 500.0), 1.71e-5, places=7)
2626

27-
def test_my_pT_region3(self):
27+
def test_my_allregions_pT_region3(self):
2828
self.assertAlmostEqual(Viscosity.my_allregions_pT(19.0, 640.0), 2.54e-5, places=7)
2929

30-
def test_my_pT_region4(self):
30+
def test_my_allregions_pT_region4(self):
3131
self.assertAlmostEqual(Viscosity.my_allregions_pT(2.63890, 500.0), 0.0, places=7)
3232

33-
def test_my_pT_region5(self):
33+
def test_my_allregions_pT_region5(self):
3434
self.assertAlmostEqual(Viscosity.my_allregions_pT(9.0, 1100.0), 4.18e-5, places=7)
3535

36-
def test_my_pT_invalid_area(self):
36+
def test_my_allregions_pT_invalid_area(self):
3737
self.assertEqual(Viscosity.my_allregions_pT(9.0, 2000.0), 2015.0)
3838
self.assertEqual(Viscosity.my_allregions_pT(350.0, 900.0), 2015.0)
3939
self.assertEqual(Viscosity.my_allregions_pT(400.0, 500.0), 2015.0)
4040
self.assertEqual(Viscosity.my_allregions_pT(600.0, 300.0), 2015.0)
4141

42+
class Test_my_AllRegions_ph(unittest.TestCase):
43+
44+
def test_my_allregions_ph_region1(self):
45+
self.assertAlmostEqual(Viscosity.my_allregions_ph(1.0, 100.0), 0.000918, places=6)
46+
47+
def test_my_allregions_ph_region2(self):
48+
self.assertAlmostEqual(Viscosity.my_allregions_ph(1.0, 4000.0), 3.79e-5, places=7)
49+
50+
def test_my_allregions_ph_region3(self):
51+
self.assertAlmostEqual(Viscosity.my_allregions_ph(19.0, 2500.0), 2.58e-5, places=7)
52+
53+
def test_my_allregions_ph_region4(self):
54+
self.assertAlmostEqual(Viscosity.my_allregions_ph(1.0, 1000.0), 1.33e-5, places=7)
55+
56+
def test_my_allregions_ph_region5(self):
57+
self.assertAlmostEqual(Viscosity.my_allregions_ph(1.0, 4250.0), 4.19e-5, places=7)
58+
59+
def test_my_ph_invalid_area(self):
60+
self.assertEqual(Viscosity.my_allregions_ph(9.0, 6590.0), 2015.0)
61+
self.assertEqual(Viscosity.my_allregions_ph(350.0, 3490.0), 2015.0)
62+
self.assertEqual(Viscosity.my_allregions_ph(400.0, 990.0), 2015.0)
63+
self.assertEqual(Viscosity.my_allregions_ph(600.0, 170.0), 2015.0)
64+
4265
if __name__ == '__main__':
4366
unittest.main()

0 commit comments

Comments
 (0)