research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[Thermo] Fix finite-difference properties in two-phase region

Browse source 
These properties are actually infinite in the two-phase region, but attempting
to compute them by finite difference would incorrectly give a finite result, so
they need to be treated as a special case.
This commit is contained in:
Ray Speth 2016-10-20 15:50:08 -04:00
parent a92245f683
commit 3f3b96587b
2 changed files with 29 additions and 6 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

7
interfaces/cython/cantera/test/test_purefluid.py
View file

@ -132,6 +132,13 @@ class TestPureFluid(utilities.CanteraTest):
self.water.TP = 450, 12
self.assertNear(ref.thermal_expansion_coeff,
self.water.thermal_expansion_coeff, 1e-5)
def test_fd_properties_twophase(self):
self.water.TX = 400, 0.1
self.assertEqual(self.water.cp, np.inf)
self.assertEqual(self.water.isothermal_compressibility, np.inf)
self.assertEqual(self.water.thermal_expansion_coeff, np.inf)
def test_TPX(self):
self.water.TX = 400, 0.8
T,P,X = self.water.TPX

28
src/tpx/Sub.cpp
View file

@ -90,10 +90,14 @@ double Substance::cp()
double T2 = std::min(Tmax(), Tsave + dt);
double p0 = P();
double x0 = x();
if (TwoPhase()) {
// In the two-phase region, cp is infinite
return std::numeric_limits<double>::infinity();
}
Set(PropertyPair::TP, T1, p0);
double x1 = x();
if ((x0 == 1.0 || x0 == 0.0) && x1 != x0) {
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;
@ -103,7 +107,7 @@ double Substance::cp()
Set(PropertyPair::TP, T2, p0);
double x2 = x();
if ((x0 == 1.0 || x0 == 0.0) && x2 != x0) {
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;
@ -123,9 +127,15 @@ double Substance::thermalExpansionCoeff()
double p0 = P();
double x0 = x();
if (TwoPhase()) {
// In the two-phase region, the thermal expansion coefficient is
// infinite
return std::numeric_limits<double>::infinity();
}
Set(PropertyPair::TP, T1, p0);
double x1 = x();
if ((x0 == 1.0 || x0 == 0.0) && x1 != x0) {
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;
@ -135,7 +145,7 @@ double Substance::thermalExpansionCoeff()
Set(PropertyPair::TP, T2, p0);
double x2 = x();
if ((x0 == 1.0 || x0 == 0.0) && x2 != x0) {
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;
@ -151,13 +161,19 @@ double Substance::isothermalCompressibility()
{
double Psave = P(), dp = 1.e-4*Psave;
double x0 = x();
if (TwoPhase()) {
// In the two-phase region, the isothermal compressibility is infinite
return std::numeric_limits<double>::infinity();
}
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 (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;
@ -167,7 +183,7 @@ double Substance::isothermalCompressibility()
Set(PropertyPair::TP, T, P2);
double x2 = x();
if ((x0 == 1.0 || x0 == 0.0) && x2 != x0) {
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;