diff --git a/include/cantera/tpx/Sub.h b/include/cantera/tpx/Sub.h index e5e3841b3ef..c445a56e8f0 100644 --- a/include/cantera/tpx/Sub.h +++ b/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; diff --git a/interfaces/cython/cantera/examples/thermo/vapordome.py b/interfaces/cython/cantera/examples/thermo/vapordome.py new file mode 100644 index 00000000000..e759905b3c3 --- /dev/null +++ b/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() diff --git a/interfaces/cython/cantera/test/test_purefluid.py b/interfaces/cython/cantera/test/test_purefluid.py index 73501d2b27b..dc910f12ceb 100644 --- a/interfaces/cython/cantera/test/test_purefluid.py +++ b/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 diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp index 718e02b2511..b5d88373f0a 100644 --- a/src/tpx/Sub.cpp +++ b/src/tpx/Sub.cpp @@ -54,7 +54,14 @@ const double DeltaT = 0.000001; double Substance::cv() { - double Tsave = T, dt = 1.e-4*T; + if (TwoPhase(true)) { + // While cv can be calculated for the two-phase region (the state can + // always be continuously varied along an isochor on a T-v diagram), + // this calculation is currently not implemented + return std::numeric_limits::quiet_NaN(); + } + + double Tsave = T, dt = 1.e-4 * T; double x0 = x(); double T1 = std::max(Tmin(), Tsave - dt); double T2 = std::min(Tmax(), Tsave + dt); @@ -85,136 +92,184 @@ double Substance::cv() double Substance::cp() { - double Tsave = T, dt = 1.e-4*T; - double RhoSave = Rho; - double T1 = std::max(Tmin(), Tsave - dt); - double T2 = std::min(Tmax(), Tsave + dt); - double p0 = P(); - double x0 = x(); - if (TwoPhase()) { + if (TwoPhase(true)) { // In the two-phase region, cp is infinite return std::numeric_limits::infinity(); } - Set(PropertyPair::TP, T1, p0); - double x1 = x(); - if (x1 != x0) { - // If the initial state was pure liquid or pure vapor, and the state at - // T-dT is not, just take a one-sided difference - T1 = Tsave; - Set(PropertyPair::TP, T1, p0); - } - double s1 = s(); + double Tsave = T, dt = 1.e-4 * T; + double RhoSave = Rho; + double T1, T2, s1, s2; + double p0 = P(); - Set(PropertyPair::TP, T2, p0); - double x2 = x(); - if (x2 != x0) { - // If the initial state was pure liquid or pure vapor, and the state at - // T+dT is not, just take a one-sided difference - T2 = Tsave; + if (Rho >= Rhf) { + // initial state is pure liquid + T1 = std::max(Tmin(), Tsave - dt); + Set(PropertyPair::TP, T1, p0); + s1 = s(); + try { + T2 = std::min(Tsat(p0), Tsave + dt); + } catch (CanteraError&) { + // Tsat does not exist beyond two-phase region + T2 = Tsave + dt; + } + if (T2 < Tsave + dt) { + Set(PropertyPair::TX, T2, 0.); + } else { + Set(PropertyPair::TP, T2, p0); + } + s2 = s(); + } else { + // initial state is pure vapor + try { + T1 = std::max(Tsat(p0), Tsave - dt); + } catch (CanteraError&) { + // Tsat does not exist beyond two-phase region + T1 = Tsave - dt; + } + if (T1 > Tsave - dt) { + Set(PropertyPair::TX, T1, 1.); + } else { + Set(PropertyPair::TP, T1, p0); + } + s1 = s(); + T2 = std::min(Tmax(), Tsave + dt); Set(PropertyPair::TP, T2, p0); + s2 = s(); } - double s2 = s(); Set(PropertyPair::TV, Tsave, 1.0 / RhoSave); - return T*(s2 - s1)/(T2-T1); + return T * (s2 - s1) / (T2 - T1); } double Substance::thermalExpansionCoeff() { - double Tsave = T, dt = 1.e-4*T; - double RhoSave = Rho; - double T1 = std::max(Tmin(), Tsave - dt); - double T2 = std::min(Tmax(), Tsave + dt); - double p0 = P(); - double x0 = x(); - - if (TwoPhase()) { + if (TwoPhase(true)) { // In the two-phase region, the thermal expansion coefficient is // infinite return std::numeric_limits::infinity(); } - Set(PropertyPair::TP, T1, p0); - double x1 = x(); - if (x1 != x0) { - // If the initial state was pure liquid or pure vapor, and the state at - // T-dT is not, just take a one-sided difference - T1 = Tsave; - Set(PropertyPair::TP, T1, p0); - } - double v1 = v(); + double Tsave = T, dt = 1.e-4 * T; + double RhoSave = Rho; + double T1, T2, v1, v2; + double p0 = P(); - Set(PropertyPair::TP, T2, p0); - double x2 = x(); - if (x2 != x0) { - // If the initial state was pure liquid or pure vapor, and the state at - // T+dT is not, just take a one-sided difference - T2 = Tsave; + if (Rho >= Rhf) { + // initial state is pure liquid + T1 = std::max(Tmin(), Tsave - dt); + Set(PropertyPair::TP, T1, p0); + v1 = v(); + try { + T2 = std::min(Tsat(p0), Tsave + dt); + } catch (CanteraError&) { + // Tsat does not exist beyond two-phase region + T2 = Tsave + dt; + } + if (T2 < Tsave + dt) { + Set(PropertyPair::TX, T2, 0.); + } else { + Set(PropertyPair::TP, T2, p0); + } + v2 = v(); + } else { + // initial state is pure vapor + try { + T1 = std::max(Tsat(p0), Tsave - dt); + } catch (CanteraError&) { + // Tsat does not exist beyond two-phase region + T1 = Tsave - dt; + } + if (T1 > Tsave - dt) { + Set(PropertyPair::TX, T1, 1.); + } else { + Set(PropertyPair::TP, T1, p0); + } + v1 = v(); + T2 = std::min(Tmax(), Tsave + dt); Set(PropertyPair::TP, T2, p0); + v2 = v(); } - double v2 = v(); Set(PropertyPair::TV, Tsave, 1.0 / RhoSave); - return 2.0*(v2 - v1)/((v2 + v1)*(T2-T1)); + return 2.0 * (v2 - v1) / ((v2 + v1) * (T2 - T1)); } double Substance::isothermalCompressibility() { - double Psave = P(), dp = 1.e-4*Psave; - double RhoSave = Rho; - double x0 = x(); - - if (TwoPhase()) { + if (TwoPhase(true)) { // In the two-phase region, the isothermal compressibility is infinite return std::numeric_limits::infinity(); } + double Psave = P(), dp = 1.e-4 * Psave; + double RhoSave = Rho; + double P1, P2, v1, v2; double v0 = v(); - double P1 = Psave - dp; - double P2 = Psave + dp; - - Set(PropertyPair::TP, T, P1); - double x1 = x(); - if (x1 != x0) { - // If the initial state was pure liquid or pure vapor, and the state at - // P-dP is not, just take a one-sided difference - P1 = Psave; - Set(PropertyPair::TP, T, P1); - } - double v1 = v(); - Set(PropertyPair::TP, T, P2); - double x2 = x(); - if (x2 != x0) { - // If the initial state was pure liquid or pure vapor, and the state at - // P+dP is not, just take a one-sided difference - P2 = Psave; + if (Rho >= Rhf) { + // initial state is pure liquid + P1 = Psave - dp; + if (T > Tmin() && T <= Tcrit()) { + P1 = std::max(Ps(), P1); + } + if (P1 > Psave - dp) { + Set(PropertyPair::PX, P1, 0.); + } else { + Set(PropertyPair::TP, T, P1); + } + v1 = v(); + P2 = Psave + dp; Set(PropertyPair::TP, T, P2); + v2 = v(); + } else { + // initial state is pure vapor + P1 = std::max(1.e-7, Psave - dp); + Set(PropertyPair::TP, T, P1); + v1 = v(); + P2 = Psave + dp; + if (T > Tmin() && T <= Tcrit()) { + P2 = std::min(Ps(), P2); + } + if (P2 < Psave + dp) { + Set(PropertyPair::PX, P2, 1.); + } else { + Set(PropertyPair::TP, T, P2); + } + v2 = v(); } - double v2 = v(); Set(PropertyPair::TV, T, 1.0 / RhoSave); - return -(v2 - v1)/(v0*(P2-P1)); + return -(v2 - v1) / (v0 * (P2 - P1)); } double Substance::dPsdT() { double tsave = T; double ps1 = Ps(); - T += DeltaT; - double dpdt = (Ps() - ps1)/DeltaT; + double dpdt; + if (T + DeltaT < Tcrit()) { + T += DeltaT; + dpdt = (Ps() - ps1)/DeltaT; + } else { + T -= DeltaT; + // Ps() will fail for illegal temperature + dpdt = (ps1 - Ps())/DeltaT; + } T = tsave; return dpdt; } -int Substance::TwoPhase() +int Substance::TwoPhase(bool strict) { if (T >= Tcrit()) { return 0; } update_sat(); - return ((Rho < Rhf) && (Rho > Rhv) ? 1 : 0); + if (strict) { + return ((Rho > Rhv) && (Rho < Rhf) ? 1 : 0); + } + return ((Rho >= Rhv) && (Rho <= Rhf) ? 1 : 0); } double Substance::x() @@ -237,26 +292,29 @@ double Substance::x() double Substance::Tsat(double p) { - if (p <= 0.0 || p > Pcrit()) { - throw CanteraError("Substance::Tsat", "illegal pressure value"); + double Tsave = T; + if (p <= 0. || p > Pcrit()) { + // @TODO: this check should also fail below the triple point pressure + // (calculated as Ps() for Tmin() as it is not available as a numerical + // parameter). A bug causing low-temperature TP setters to fail for + // Heptane (GitHub issue #605) prevents this at the moment. + // Without this check, a calculation attempt fails with the less + // meaningful error "No convergence" below. + throw CanteraError("Substance::Tsat", + "Illegal pressure value: {}", p); } + T = Tsave; + Tsave = T; + int LoopCount = 0; double tol = 1.e-6*p; - double Tsave = T; - if (T < Tmin()) { - T = 0.5*(Tcrit() - Tmin()); - } - if (T >= Tcrit()) { + if (T < Tmin() || T > Tcrit()) { T = 0.5*(Tcrit() - Tmin()); } double dp = 10*tol; while (fabs(dp) > tol) { - if (T > Tcrit()) { - T = Tcrit() - 0.001; - } - if (T < Tmin()) { - T = Tmin() + 0.001; - } + T = std::min(T, Tcrit()); + T = std::max(T, Tmin()); dp = p - Ps(); double dt = dp/dPsdT(); double dta = fabs(dt); @@ -268,7 +326,7 @@ double Substance::Tsat(double p) LoopCount++; if (LoopCount > 100) { T = Tsave; - throw CanteraError("Substance::Tsat", "No convergence"); + throw CanteraError("Substance::Tsat", "No convergence: p = {}", p); } } double tsat = T; @@ -322,15 +380,17 @@ void Substance::Set(PropertyPair::type XY, double x0, double y0) x0, y0, TolAbsP, TolAbsV, TolRel, TolRel); break; case PropertyPair::TP: + set_T(x0); if (x0 < Tcrit()) { - set_T(x0); - if (y0 < Ps()) { + if (fabs(y0 - Ps()) / y0 < TolRel) { + throw CanteraError("Substance::Set", + "Saturated mixture detected: use vapor " + "fraction to specify state instead"); + } else if (y0 < Ps()) { Set(PropertyPair::TX, x0, 1.0); } else { Set(PropertyPair::TX, x0, 0.0); } - } else { - set_T(x0); } set_xy(propertyFlag::T, propertyFlag::P, x0, y0, TolAbsT, TolAbsP, TolRel, TolRel); @@ -370,31 +430,35 @@ void Substance::Set(PropertyPair::type XY, double x0, double y0) x0, y0, TolAbsS, TolAbsH, TolRel, TolRel); break; case PropertyPair::PX: - temp = Tsat(x0); + temp = Tmin(); + set_T(temp); if (y0 > 1.0 || y0 < 0.0) { throw CanteraError("Substance::Set", "Invalid vapor fraction, {}", y0); - } else if (temp >= Tcrit()) { + } + if (x0 < Ps() || x0 > Pcrit()) { throw CanteraError("Substance::Set", - "Can't set vapor fraction above the critical point"); - } else { - set_T(temp); - update_sat(); - Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv); + "Illegal pressure value: {} (supercritical " + "or below triple point)", x0); } + temp = Tsat(x0); + set_T(temp); + update_sat(); + Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv); break; case PropertyPair::TX: if (y0 > 1.0 || y0 < 0.0) { throw CanteraError("Substance::Set", "Invalid vapor fraction, {}", y0); - } else if (x0 >= Tcrit()) { + } + if (x0 < Tmin() || x0 > Tcrit()) { throw CanteraError("Substance::Set", - "Can't set vapor fraction above the critical point"); - } else { - set_T(x0); - update_sat(); - Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv); + "Illegal temperature value: {} " + "(supercritical or below triple point)", x0); } + set_T(x0); + update_sat(); + Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv); break; default: throw CanteraError("Substance::Set", "Invalid input."); @@ -417,7 +481,7 @@ void Substance::set_T(double t0) if ((t0 >= Tmin()) && (t0 <= Tmax())) { T = t0; } else { - throw CanteraError("Substance::set_T", "illegal temperature: {}", t0); + throw CanteraError("Substance::set_T", "Illegal temperature: {}", t0); } } @@ -427,7 +491,7 @@ void Substance::set_v(double v0) Rho = 1.0/v0; } else { throw CanteraError("Substance::set_v", - "negative specific volume: {}", v0); + "Negative specific volume: {}", v0); } } @@ -435,7 +499,7 @@ double Substance::Ps() { if (T < Tmin() || T > Tcrit()) { throw CanteraError("Substance::Ps", - "illegal temperature value {}", T); + "Illegal temperature value: {}", T); } update_sat(); return Pst; @@ -443,9 +507,9 @@ double Substance::Ps() void Substance::update_sat() { - if ((T != Tslast) && (T < Tcrit())) { + if (T != Tslast) { double Rho_save = Rho; - double pp = Psat(); + double pp = Psat(); // illegal temperatures are caught by Psat() double lps = log(pp); // trial value = Psat from correlation int i; @@ -474,7 +538,7 @@ void Substance::update_sat() dg = - dg; } - if (fabs(dg) < 0.001 && Rhf > Rhv) { + if (fabs(dg) < 0.001) { break; } double dp = dg/(1.0/Rhv - 1.0/Rhf); @@ -494,13 +558,9 @@ void Substance::update_sat() lps = log(pp); } } - if (Rhf <= Rhv) { - throw CanteraError("Substance::update_sat", - "wrong root found for sat. liquid or vapor at P = {}", pp); - } if (i >= 20) { - throw CanteraError("Substance::update_sat", "no convergence"); + throw CanteraError("Substance::update_sat", "No convergence"); } else { Pst = pp; Tslast = T; @@ -523,7 +583,7 @@ double Substance::vprop(propertyFlag::type ijob) case propertyFlag::P: return Pp(); default: - throw CanteraError("Substance::vprop", "invalid job index"); + throw CanteraError("Substance::vprop", "Invalid job index"); } } @@ -552,7 +612,7 @@ int Substance::Lever(int itp, double sat, double val, propertyFlag::type ifunc) return 0; } } else { - throw CanteraError("Substance::Lever","general error"); + throw CanteraError("Substance::Lever", "General error"); } Set(PropertyPair::TX, T, 1.0); double Valg = vprop(ifunc); @@ -761,7 +821,7 @@ void Substance::set_TPp(double Temp, double Pressure) Pmax = P_here; } if (Vmin >= Vmax) { - throw CanteraError("Substance::set_TPp","Vmin >= Vmax"); + throw CanteraError("Substance::set_TPp", "Vmin >= Vmax"); } else if ((Vmin > 0.0) && (Vmax < Big)) { kbr = 1; } @@ -793,14 +853,15 @@ void Substance::set_TPp(double Temp, double Pressure) } v_here += dv; if (dv == 0.0) { - throw CanteraError("Substance::set_TPp","dv = 0 and no convergence"); + throw CanteraError("Substance::set_TPp", "dv = 0 and no convergence"); } Set(PropertyPair::TV, Temp, v_here); LoopCount++; - if (LoopCount > 100) { + if (LoopCount > 200) { Set(PropertyPair::TV, Temp, v_save); throw CanteraError("Substance::set_TPp", - "no convergence for P* = {}, V* = {}", + "No convergence for P = {}, T = {}\n" + "(P* = {}, V* = {})", Pressure, Temp, Pressure/Pcrit(), v_save/Vcrit()); } } diff --git a/src/tpx/Water.cpp b/src/tpx/Water.cpp index 9079ce5f46d..8b2abf569e2 100644 --- a/src/tpx/Water.cpp +++ b/src/tpx/Water.cpp @@ -176,7 +176,7 @@ double water::ldens() { double sum=0; int i; - if ((T < Tmn) || (T >= Tc)) { + if ((T < Tmn) || (T > Tc)) { throw CanteraError("water::ldens", "Temperature out of range. T = {}", T); }