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

[Thermo] Fix calculation of derivative properties of pure fluids

Browse source 
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
This commit is contained in:
Ray Speth 2016-10-19 17:32:26 -04:00
parent b1a1ce6e23
commit a92245f683
2 changed files with 90 additions and 3 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -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')

81
src/tpx/Sub.cpp
View file

@ -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()