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

Improved convergence when solving for pure fluid saturated states

Browse source
This commit is contained in:
Ray Speth 2012-05-11 15:27:33 +00:00
parent 0901c8b25f
commit c7af9ddf81
2 changed files with 87 additions and 9 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

29
src/tpx/Sub.cpp
View file

@ -374,19 +374,32 @@ void Substance::update_sat()
int i;
for (i = 0; i<20; i++) {
Rho = ldens(); // trial value = liquid density
if (i==0) {
Rho = ldens(); // trial value = liquid density
} else {
Rho = Rhf;
}
set_TPp(T,pp);
Rhf = Rho; // sat liquid density
gf = hp() - T*sp();
Rho = pp*MolWt()/(8314.0*T); // trial value = ideal gas
if (i==0) {
Rho = pp*MolWt()/(8314.0*T); // trial value = ideal gas
} else {
Rho = Rhv;
}
set_TPp(T,pp);
Rhv = Rho; // sat vapor density
gv = hp() - T*sp();
dg = gv - gf;
if (fabs(dg) < 0.001) {
if (Rhv > Rhf) {
std::swap(Rhv, Rhf);
dg = - dg;
}
if (fabs(dg) < 0.001 && Rhf > Rhv) {
break;
}
dp = dg/(1.0/Rhv - 1.0/Rhf);
@ -398,14 +411,14 @@ void Substance::update_sat()
pp = exp(lps);
} else {
pp -= dp;
lps = log(pp); // added 10/5/04
lps = log(pp);
}
if (pp > Pcrit()) {
pp = psold + 0.5*(Pcrit() - psold);
lps = log(pp); // added 10/5/04
lps = log(pp);
} else if (pp < 0.0) {
pp = psold/2.0;
lps = log(pp); // added 10/5/04
lps = log(pp);
}
}
if (Rhf <= Rhv) {
@ -673,8 +686,8 @@ void Substance::set_TPp(double Temp, double Pressure)
v_here = vp();
// loop
while (P_here = Pp(), fabs(Pressure - P_here) >= ErrP*Pressure) {
while (P_here = Pp(),
fabs(Pressure - P_here) >= ErrP*Pressure || LoopCount == 0) {
if (P_here < 0.0) {
BracketSlope(Pressure);
} else {

67
test/python/testPureFluid.py
View file

@ -1,11 +1,12 @@
from __future__ import division
import itertools
import utilities
import numpy as np
import Cantera as ct
import Cantera.liquidvapor as lv
# To minimize data when transcribing tabulated data, the input units here are:
# To minimize errors when transcribing tabulated data, the input units here are:
# T: K, P: MPa, rho: kg/m3, v: m3/kg, (u,h): kJ/kg, s: kJ/kg-K
# Which are then converted to SI
class StateData(object):
@ -294,3 +295,67 @@ class Water(PureFluidTestCases, utilities.CanteraTest):
tols = Tolerances(2e-3, 2e-3, 2e-3)
PureFluidTestCases.__init__(self, 'water', refState, tols)
utilities.CanteraTest.__init__(self, *args, **kwargs)
class Convergence(utilities.CanteraTest):
def setUp(self):
self.fluid = lv.Water()
def test_TP(self):
# Focus on the region near the critical point
TT = [273.161, 300.0, 350.0, 400.0, 500.0,
600.0, 640.0, 645.0, 646.0, 647.0,
647.1, 647.2, 647.22, 647.23, 647.25,
647.26, 647.27, 647.28, 647.282, 647.284,
647.285, 647.286, 647.287, 650.0, 800.0]
PP = [1234.0, 101325.0, 5e5, 22.0e6, 22.08e6, 22.09e6, 10001000.0]
errors = ''
nErrors = 0
for T,P in itertools.product(TT,PP):
try:
self.fluid.set(T=T, P=P)
self.assertNear(self.fluid.temperature(), T, 1e-6)
self.assertNear(self.fluid.pressure(), P, 1e-6)
except Exception as e:
errors += 'Error at T=%r, P=%r:\n%s\n\n' % (T,P,e)
nErrors += 1
if errors:
errors += 'Total error count:%s\n' % nErrors
raise AssertionError(errors)
def test_UV(self):
u0 = -1.58581e7
UU = np.array([0, 100, 200, 500, 1000, 1500, 2000]) * 1000 + u0
VV = [0.001, 0.002, 0.005, 0.010, 0.10, 0.5, 1.0, 1.5, 2.0]
errors = ''
nErrors = 0
for u,v in itertools.product(UU,VV):
try:
self.fluid.set(U=u, V=v)
self.assertNear(self.fluid.intEnergy_mass(), u, 1e-6)
self.assertNear(self.fluid.volume_mass(), v, 1e-6)
except Exception as e:
errors += 'Error at u=%r, v=%r:\n%s\n\n' % (u,v,e)
nErrors += 1
if errors:
errors += 'Total error count:%s\n' % nErrors
raise AssertionError(errors)
def test_HP(self):
h0 = -1.58581e7
HH = np.array([0, 100, 200, 500, 1000, 1500, 2000]) * 1000 + h0
PP = [1234.0, 101325.0, 5e5, 22.0e6, 22.08e6, 22.09e6, 10001000.0]
errors = ''
nErrors = 0
for h,P in itertools.product(HH,PP):
try:
self.fluid.set(H=h, P=P)
self.assertNear(self.fluid.enthalpy_mass(), h, 1e-6)
self.assertNear(self.fluid.pressure(), P, 1e-6)
except Exception as e:
errors += 'Error at h=%r, P=%r:\n%s\n\n' % (h,P,e)
nErrors += 1
if errors:
errors += 'Total error count:%s\n' % nErrors
raise AssertionError(errors)