[Reactor] Fix using pure substances near temperature limits

Fixes #475
This commit is contained in:
Ray Speth 2017-12-02 16:49:34 -05:00
parent 2ac8a0ef08
commit deeaaed03f
2 changed files with 38 additions and 5 deletions

View file

@ -1468,6 +1468,32 @@ class PureFluidReactorTest(utilities.CanteraTest):
self.assertEqual(states.X[-1], 1)
self.assertNear(states.X[30], 0.54806, 1e-4)
def test_Reactor_2(self):
phase = ct.PureFluid('liquidvapor.xml', 'carbondioxide')
air = ct.Solution('air.xml')
phase.TP = 218, 5e6
r1 = ct.Reactor(phase)
r1.volume = 0.1
air.TP = 500, 5e6
r2 = ct.Reactor(air)
r2.volume = 10.0
w1 = ct.Wall(r1, r2, U=10000, A=1)
w1.expansion_rate_coeff = 1e-3
net = ct.ReactorNet([r1,r2])
states = ct.SolutionArray(phase, extra='t')
for t in np.arange(0.0, 60.0, 1):
net.advance(t)
states.append(TD=r1.thermo.TD, t=net.time)
self.assertEqual(states.X[0], 0)
self.assertEqual(states.X[-1], 1)
self.assertNear(states.X[20], 0.644865, 1e-4)
def test_ConstPressureReactor(self):
phase = ct.Nitrogen()
air = ct.Solution('air.xml')

View file

@ -144,13 +144,20 @@ void Reactor::updateState(doublereal* y)
TT = bmt::bracket_and_solve_root(
u_err, T, 1.2, true, bmt::eps_tolerance<double>(48), maxiter);
} catch (std::exception& err) {
// Set m_thermo back to a reasonable state if root finding fails
m_thermo->setState_TR(T, m_mass / m_vol);
throw CanteraError("Reactor::updateState", err.what());
// Try full-range bisection if bracketing fails (e.g. near
// temperature limits for the phase's equation of state)
try {
TT = bmt::bisect(u_err, m_thermo->minTemp(), m_thermo->maxTemp(),
bmt::eps_tolerance<double>(48), maxiter);
} catch (std::exception& err2) {
// Set m_thermo back to a reasonable state if root finding fails
m_thermo->setState_TR(T, m_mass / m_vol);
throw CanteraError("Reactor::updateState",
"{}\nat U = {}, rho = {}", err2.what(), U, m_mass / m_vol);
}
}
if (fabs(TT.first - TT.second) > 1e-7*TT.first) {
throw CanteraError("Reactor::updateState",
"bracket_and_solve_root failed");
throw CanteraError("Reactor::updateState", "root finding failed");
}
m_thermo->setState_TR(TT.second, m_mass / m_vol);
} else {