Fix edge cases for saturated PureFluid states / illustrate saturated water properties

include/cantera/tpx/Sub.h

@@ -146,8 +146,10 @@ class Substance
     //! is returned if v > Vcrit.
     double x();
 
-    //! Returns 1 if the current state is a liquid/vapor mixture, 0 otherwise
-    int TwoPhase();
+    //! Returns 1 if the current state is a liquid/vapor mixture, 0 otherwise.
+    //! By default, saturated vapor and saturated liquid are included; setting
+    //! the flag *strict* to true will exclude the boundaries.
+    int TwoPhase(bool strict=false);
     //! @}
 
     virtual double Pp()=0;

interfaces/cython/cantera/examples/thermo/vapordome.py

@@ -0,0 +1,94 @@
+"""
+This example generates a saturated steam table and plots the vapor dome. The
+steam table corresponds to data typically found in thermodynamic text books
+and uses the same customary units.
+
+Requires: Cantera >= 2.5.0, matplotlib >= 2.0, pandas >= 1.1.0, numpy >= 1.12
+"""
+
+import cantera as ct
+import pandas as pd
+import numpy as np
+from matplotlib import pyplot as plt
+
+w = ct.Water()
+
+# create colums
+columns = ['T', 'P',
+           'vf', 'vfg', 'vg',
+           'uf', 'ufg', 'ug',
+           'hf', 'hfg', 'hg',
+           'sf', 'sfg', 'sg']
+
+# temperatures correspond to Engineering Thermodynamics, Moran et al. (9th ed),
+# Table A-2; additional data points are added close to the critical point;
+# w.min_temp is equal to the triple point temperature
+degc = np.hstack([np.array([w.min_temp - 273.15, 4, 5, 6, 8]),
+                  np.arange(10, 37), np.array([38]),
+                  np.arange(40, 100, 5), np.arange(100, 300, 10),
+                  np.arange(300, 380, 20), np.arange(370, 374),
+                  np.array([w.critical_temperature - 273.15])])
+
+df = pd.DataFrame(0, index=np.arange(len(degc)), columns=columns)
+df.T = degc
+
+arr = ct.SolutionArray(w, len(degc))
+
+# saturated vapor data
+arr.TQ = degc + 273.15, 1
+df.P = arr.P_sat / 1.e5
+df.vg = arr.v
+df.ug = arr.int_energy_mass / 1.e3
+df.hg = arr.enthalpy_mass / 1.e3
+df.sg = arr.entropy_mass / 1.e3
+
+# saturated liquid data
+arr.TQ = degc + 273.15, 0
+df.vf = arr.v
+df.uf = arr.int_energy_mass / 1.e3
+df.hf = arr.enthalpy_mass / 1.e3
+df.sf = arr.entropy_mass / 1.e3
+
+# delta values
+df.vfg = df.vg - df.vf
+df.ufg = df.ug - df.uf
+df.hfg = df.hg - df.hf
+df.sfg = df.sg - df.sf
+
+# reference state (triple point; liquid state)
+w.TQ = w.min_temp, 0
+uf0 = w.int_energy_mass / 1.e3
+hf0 = w.enthalpy_mass / 1.e3
+sf0 = w.entropy_mass / 1.e3
+pv0 = w.P * w.v / 1.e3
+
+# change reference state
+df.ug -= uf0
+df.uf -= uf0
+df.hg -= hf0 - pv0
+df.hf -= hf0 - pv0
+df.sg -= sf0
+df.sf -= sf0
+
+# print and write saturated steam table to csv file
+print(df)
+df.to_csv('saturated_steam_T.csv', index=False)
+
+# illustrate the vapor dome in a P-v diagram
+plt.semilogx(df.vf.values, df.P.values, label='Saturated liquid')
+plt.semilogx(df.vg.values, df.P.values, label='Saturated vapor')
+plt.semilogx(df.vg.values[-1], df.P.values[-1], 'o', label='Critical point')
+plt.xlabel(r'Specific volume - $v$ ($\mathrm{m^3/kg}$)')
+plt.ylabel(r'Presssure - $P$ (bar)')
+plt.legend()
+
+# illustrate the vapor dome in a T-s diagram
+plt.figure()
+plt.plot(df.sf.values, df['T'].values, label='Saturated liquid')
+plt.plot(df.sg.values, df['T'].values, label='Saturated vapor')
+plt.plot(df.sg.values[-1], df['T'].values[-1], 'o', label='Critical point')
+plt.xlabel(r'Specific entropy - $s$ ($\mathrm{kJ/kg-K}$)')
+plt.ylabel(r'Temperature - $T$ (${}^\circ C$)')
+plt.legend()
+
+plt.show()

interfaces/cython/cantera/test/test_purefluid.py

@@ -34,6 +34,8 @@ def test_substance_set(self):
         self.water.TV = 400, 1.45
         self.assertNear(self.water.T, 400)
         self.assertNear(self.water.v, 1.45)
+        with self.assertRaisesRegex(ct.CanteraError, 'Negative specific volume'):
+            self.water.TV = 300, -1.
 
         self.water.PV = 101325, 1.45
         self.assertNear(self.water.P, 101325)
@@ -116,10 +118,10 @@ def check_fd_properties(self, T1, P1, T2, P2, tol):
         self.assertNear(k1, k2, tol)
         self.assertNear(alpha1, alpha2, tol)
 
-        # calculating these finite difference properties should not perturbe the
-        # state of the object
-        self.assertEqual(h1a, h1b)
-        self.assertEqual(h2a, h2b)
+        # calculating these finite difference properties should not perturb the
+        # state of the object (except for checks on edge cases)
+        self.assertNear(h1a, h1b, 1e-9)
+        self.assertNear(h2a, h2b, 1e-9)
 
     def test_properties_near_min(self):
         self.check_fd_properties(self.water.min_temp*(1+1e-5), 101325,
@@ -162,12 +164,121 @@ def test_thermal_expansion_coeff_TD(self):
             self.water.TD = T, 0.1
             self.assertNear(T * self.water.thermal_expansion_coeff, 1.0, 1e-2)
 
-    def test_fd_properties_twophase(self):
-        self.water.TQ = 400, 0.1
+    def test_pq_setter_triple_check(self):
+        self.water.PQ = 101325, .2
+        T = self.water.T
+        # change T such that it would result in a Psat larger than P
+        self.water.TP = 400, 101325
+        # ensure that correct triple point pressure is recalculated
+        # (necessary as this value is not stored by the C++ base class)
+        self.water.PQ = 101325, .2
+        self.assertNear(T, self.water.T, 1e-9)
+        with self.assertRaisesRegex(ct.CanteraError, 'below triple point'):
+            # min_temp is triple point temperature
+            self.water.TP = self.water.min_temp, 101325
+            P = self.water.P_sat # triple-point pressure
+            self.water.PQ = .999*P, .2
+
+    def test_quality_exceptions(self):
+        # Critical point
+        self.water.TP = 300, ct.one_atm
+        self.water.TQ = self.water.critical_temperature, .5
+        self.assertNear(self.water.P, self.water.critical_pressure)
+        self.water.TP = 300, ct.one_atm
+        self.water.PQ = self.water.critical_pressure, .5
+        self.assertNear(self.water.T, self.water.critical_temperature)
+
+        # Supercritical
+        with self.assertRaisesRegex(ct.CanteraError, 'supercritical'):
+            self.water.TQ = 1.001 * self.water.critical_temperature, 0.
+        with self.assertRaisesRegex(ct.CanteraError, 'supercritical'):
+            self.water.PQ = 1.001 * self.water.critical_pressure, 0.
+
+        # Q negative
+        with self.assertRaisesRegex(ct.CanteraError, 'Invalid vapor fraction'):
+            self.water.TQ = 373.15, -.001
+        with self.assertRaisesRegex(ct.CanteraError, 'Invalid vapor fraction'):
+            self.water.PQ = ct.one_atm, -.001
+
+        # Q larger than one
+        with self.assertRaisesRegex(ct.CanteraError, 'Invalid vapor fraction'):
+            self.water.TQ = 373.15, 1.001
+        with self.assertRaisesRegex(ct.CanteraError, 'Invalid vapor fraction'):
+            self.water.PQ = ct.one_atm, 1.001
+
+    def test_saturated_mixture(self):
+        self.water.TP = 300, ct.one_atm
+        with self.assertRaisesRegex(ct.CanteraError, 'Saturated mixture detected'):
+            self.water.TP = 300, self.water.P_sat
+
+        w = ct.Water()
+
+        # Saturated vapor
+        self.water.TQ = 373.15, 1.
+        self.assertEqual(self.water.phase_of_matter, 'liquid-gas-mix')
+        w.TP = self.water.T, .999 * self.water.P_sat
+        self.assertNear(self.water.cp, w.cp, 1.e-3)
+        self.assertNear(self.water.cv, w.cv, 1.e-3)
+        self.assertNear(self.water.thermal_expansion_coeff, w.thermal_expansion_coeff, 1.e-3)
+        self.assertNear(self.water.isothermal_compressibility, w.isothermal_compressibility, 1.e-3)
+
+        # Saturated mixture
+        self.water.TQ = 373.15, .5
+        self.assertEqual(self.water.phase_of_matter, 'liquid-gas-mix')
         self.assertEqual(self.water.cp, np.inf)
+        self.assertTrue(np.isnan(self.water.cv))
         self.assertEqual(self.water.isothermal_compressibility, np.inf)
         self.assertEqual(self.water.thermal_expansion_coeff, np.inf)
 
+        # Saturated liquid
+        self.water.TQ = 373.15, 0.
+        self.assertEqual(self.water.phase_of_matter, 'liquid-gas-mix')
+        w.TP = self.water.T, 1.001 * self.water.P_sat
+        self.assertNear(self.water.cp, w.cp, 1.e-3)
+        self.assertNear(self.water.cv, w.cv, 1.e-3)
+        self.assertNear(self.water.thermal_expansion_coeff, w.thermal_expansion_coeff, 1.e-3)
+        self.assertNear(self.water.isothermal_compressibility, w.isothermal_compressibility, 1.e-3)
+
+    def test_saturation_near_limits(self):
+        # Low temperature limit (triple point)
+        self.water.TP = 300, ct.one_atm
+        self.water.P_sat # ensure that solver buffers sufficiently different values
+        self.water.TP = self.water.min_temp, ct.one_atm
+        psat = self.water.P_sat
+        self.water.TP = 300, ct.one_atm
+        self.water.P_sat # ensure that solver buffers sufficiently different values
+        self.water.TP = 300, psat
+        self.assertNear(self.water.T_sat, self.water.min_temp)
+
+        # High temperature limit (critical point) - saturation temperature
+        self.water.TP = 300, ct.one_atm
+        self.water.P_sat # ensure that solver buffers sufficiently different values
+        self.water.TP = self.water.critical_temperature, self.water.critical_pressure
+        self.assertNear(self.water.T_sat, self.water.critical_temperature)
+
+        # High temperature limit (critical point) - saturation pressure
+        self.water.TP = 300, ct.one_atm
+        self.water.P_sat # ensure that solver buffers sufficiently different values
+        self.water.TP = self.water.critical_temperature, self.water.critical_pressure
+        self.assertNear(self.water.P_sat, self.water.critical_pressure)
+
+        # Supercricital
+        with self.assertRaisesRegex(ct.CanteraError, 'Illegal temperature value'):
+            self.water.TP = 1.001 * self.water.critical_temperature, self.water.critical_pressure
+            self.water.P_sat
+        with self.assertRaisesRegex(ct.CanteraError, 'Illegal pressure value'):
+            self.water.TP = self.water.critical_temperature, 1.001 * self.water.critical_pressure
+            self.water.T_sat
+
+        # Below triple point
+        with self.assertRaisesRegex(ct.CanteraError, 'Illegal temperature'):
+            self.water.TP = .999 * self.water.min_temp, ct.one_atm
+            self.water.P_sat
+        # @TODO: test disabled pending fix of GitHub issue #605
+        # with self.assertRaisesRegex(ct.CanteraError, 'Illegal pressure value'):
+        #     self.water.TP = 300, .999 * psat
+        #     self.water.T_sat
+
     def test_TPQ(self):
         self.water.TQ = 400, 0.8
         T, P, Q = self.water.TPQ