From a92245f683818d3fc70129ca1150de0e9e9f0ccd Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 19 Oct 2016 17:32:26 -0400 Subject: [PATCH] [Thermo] Fix calculation of derivative properties of pure fluids Calculating specific heat capacities (cp or cv) for states near the saturation curve would give incorrect results if the finite difference method used points within the saturation region. To avoid this, we now check the points used for computing the properties and use values outside the saturation region to compute the derivatives. Fixes #273 --- .../cython/cantera/test/test_purefluid.py | 12 +++ src/tpx/Sub.cpp | 81 ++++++++++++++++++- 2 files changed, 90 insertions(+), 3 deletions(-) diff --git a/interfaces/cython/cantera/test/test_purefluid.py b/interfaces/cython/cantera/test/test_purefluid.py index bebee5fc4..6660cae2e 100644 --- a/interfaces/cython/cantera/test/test_purefluid.py +++ b/interfaces/cython/cantera/test/test_purefluid.py @@ -105,6 +105,18 @@ class TestPureFluid(utilities.CanteraTest): self.check_fd_properties(self.water.max_temp*(1-1e-5), 101325, self.water.max_temp*(1-1e-4), 101325, 1e-2) + def test_properties_near_sat1(self): + for T in [340,390,420]: + self.water.TX = T, 0.0 + P = self.water.P + self.check_fd_properties(T, P+0.01, T, P+0.5, 1e-4) + + def test_properties_near_sat2(self): + for T in [340,390,420]: + self.water.TX = T, 0.0 + P = self.water.P + self.check_fd_properties(T, P-0.01, T, P-0.5, 1e-4) + def test_isothermal_compressibility_lowP(self): # Low-pressure limit corresponds to ideal gas ref = ct.Solution('gri30.xml') diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp index a6cdf6f73..6e6f448e2 100644 --- a/src/tpx/Sub.cpp +++ b/src/tpx/Sub.cpp @@ -55,12 +55,30 @@ const double DeltaT = 0.000001; double Substance::cv() { 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); + set_T(T1); + double x1 = x(); + if ((x0 == 1.0 || x0 == 0.0) && 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_T(T1); + } double s1 = s(); + set_T(T2); + double x2 = x(); + if ((x0 == 1.0 || x0 == 0.0) && 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; + set_T(T2); + } double s2 = s(); + set_T(Tsave); return T*(s2 - s1)/(T2-T1); } @@ -71,10 +89,28 @@ double Substance::cp() double T1 = std::max(Tmin(), Tsave - dt); double T2 = std::min(Tmax(), Tsave + dt); double p0 = P(); + double x0 = x(); + Set(PropertyPair::TP, T1, p0); + double x1 = x(); + if ((x0 == 1.0 || x0 == 0.0) && 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(); + Set(PropertyPair::TP, T2, p0); + double x2 = x(); + if ((x0 == 1.0 || x0 == 0.0) && 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; + Set(PropertyPair::TP, T2, p0); + } double s2 = s(); + Set(PropertyPair::TP, Tsave, p0); return T*(s2 - s1)/(T2-T1); } @@ -85,10 +121,28 @@ double Substance::thermalExpansionCoeff() double T1 = std::max(Tmin(), Tsave - dt); double T2 = std::min(Tmax(), Tsave + dt); double p0 = P(); + double x0 = x(); + Set(PropertyPair::TP, T1, p0); + double x1 = x(); + if ((x0 == 1.0 || x0 == 0.0) && 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(); + Set(PropertyPair::TP, T2, p0); + double x2 = x(); + if ((x0 == 1.0 || x0 == 0.0) && 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; + Set(PropertyPair::TP, T2, p0); + } double v2 = v(); + Set(PropertyPair::TP, Tsave, p0); return 2.0*(v2 - v1)/((v2 + v1)*(T2-T1)); } @@ -96,12 +150,33 @@ double Substance::thermalExpansionCoeff() double Substance::isothermalCompressibility() { double Psave = P(), dp = 1.e-4*Psave; - Set(PropertyPair::TP, T, Psave - dp); + double x0 = x(); + double v0 = v(); + double P1 = Psave - dp; + double P2 = Psave + dp; + + Set(PropertyPair::TP, T, P1); + double x1 = x(); + if ((x0 == 1.0 || x0 == 0.0) && 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, Psave + dp); + + Set(PropertyPair::TP, T, P2); + double x2 = x(); + if ((x0 == 1.0 || x0 == 0.0) && 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; + Set(PropertyPair::TP, T, P2); + } double v2 = v(); + Set(PropertyPair::TP, T, Psave); - return -(v2 - v1)/((v2 + v1)*dp); + return -(v2 - v1)/(v0*(P2-P1)); } double Substance::dPsdT()