diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 3d17c49e7..7112717b5 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -19,6 +19,7 @@ const int ReservoirType = 1; const int ReactorType = 2; const int FlowReactorType = 3; const int ConstPressureReactorType = 4; +const int IdealGasReactorType = 5; /** * Base class for stirred reactors. Allows using any substance model, with diff --git a/interfaces/cython/cantera/test/test_reactor.py b/interfaces/cython/cantera/test/test_reactor.py index 030f83dcf..1e40d661c 100644 --- a/interfaces/cython/cantera/test/test_reactor.py +++ b/interfaces/cython/cantera/test/test_reactor.py @@ -149,7 +149,7 @@ class TestReactor(utilities.CanteraTest): return nSteps n_baseline = integrate(1e-10, 1e-20) - n_rtol = integrate(1e-6, 1e-20) + n_rtol = integrate(5e-7, 1e-20) n_atol = integrate(1e-10, 1e-6) self.assertTrue(n_baseline > n_rtol) @@ -315,7 +315,8 @@ class TestReactor(utilities.CanteraTest): ma = self.r1.volume * self.r1.density Ya = self.r1.Y - self.net.rtol = 1e-12 + self.net.rtol = 1e-11 + self.net.set_max_time_step(0.05) self.net.advance(2.5) mb = self.r1.volume * self.r1.density @@ -832,6 +833,8 @@ class TestReactorSensitivities(utilities.CanteraTest): (gas1, solid)) r1 = ct.Reactor(gas1) net.add_reactor(r1) + net.atol_sensitivity = 1e-10 + net.rtol_sensitivity = 1e-8 gas2 = ct.Solution('h2o2.xml') gas2.TPX = 900, 101325, 'H2:0.1, OH:1e-7, O2:0.1, AR:1e-5' @@ -867,8 +870,8 @@ class TestReactorSensitivities(utilities.CanteraTest): self.assertArrayNear(S[K2+1,:], np.zeros(2)) # Sensitivity coefficients for the disjoint reactors should be zero - self.assertNear(np.linalg.norm(S[Ns:K2,1]), 0.0) - self.assertNear(np.linalg.norm(S[K2+Ns:,0]), 0.0) + self.assertNear(np.linalg.norm(S[Ns:K2,1]), 0.0, atol=1e-5) + self.assertNear(np.linalg.norm(S[K2+Ns:,0]), 0.0, atol=1e-5) def test_parameter_order1(self): # Single reactor, changing the order in which parameters are added diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 65aab0f4e..d56782ea2 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -42,8 +42,8 @@ void Reactor::getInitialConditions(double t0, size_t leny, double* y) // set the second component to the total volume y[1] = m_vol; - // Set the third component to the temperature - y[2] = m_thermo->temperature(); + // set the third component to the total internal energy + y[2] = m_thermo->intEnergy_mass() * m_mass; // set components y+3 ... y+K+2 to the mass fractions of each species m_thermo->getMassFractions(y+3); @@ -66,7 +66,6 @@ void Reactor::initialize(doublereal t0) m_thermo->restoreState(m_state); m_sdot.resize(m_nsp, 0.0); m_wdot.resize(m_nsp, 0.0); - m_uk.resize(m_nsp, 0.0); m_nv = m_nsp + 3; for (size_t w = 0; w < m_nwalls; w++) if (m_wall[w]->surface(m_lr[w])) { @@ -123,13 +122,42 @@ void Reactor::updateState(doublereal* y) } // The components of y are [0] the total mass, [1] the total volume, - // [2] the temperature, [3...K+3] are the mass fractions of each species, - // and [K+3...] are the coverages of surface species on each wall. + // [2] the total internal energy, [3...K+3] are the mass fractions of each + // species, and [K+3...] are the coverages of surface species on each wall. m_mass = y[0]; m_vol = y[1]; m_thermo->setMassFractions_NoNorm(y+3); - m_thermo->setState_TR(y[2], m_mass / m_vol); + + if (m_energy) { + // Use Newton's method to determine the mixture temperature. Tight + // tolerances are required both for Jacobian evaluation and for + // sensitivity analysis to work correctly. + + doublereal U = y[2]; + doublereal T = temperature(); + double dT = 100; + + int i = 0; + while (abs(dT / T) > 10 * DBL_EPSILON) { + m_thermo->setState_TR(T, m_mass / m_vol); + double dUdT = m_thermo->cv_mass() * m_mass; + dT = (m_thermo->intEnergy_mass() * m_mass - U) / dUdT; + dT = std::min(dT, 0.5 * T); + T -= dT; + i++; + if (i > 100) { + std::string message = "no convergence"; + message += "\nU/m = " + fp2str(U / m_mass); + message += "\nT = " + fp2str(T); + message += "\nrho = " + fp2str(m_mass / m_vol); + message += "\n"; + throw CanteraError("Reactor::updateState", message); + } + } + } else { + m_thermo->setDensity(m_mass/m_vol); + } size_t loc = m_nsp + 3; SurfPhase* surf; @@ -171,12 +199,9 @@ void Reactor::evalEqs(doublereal time, doublereal* y, m_vdot = 0.0; m_Q = 0.0; - double mcvdTdt = 0.0; // m * c_v * dT/dt double dmdt = 0.0; // dm/dt (gas phase) double* dYdt = ydot + 3; - m_thermo->getPartialMolarIntEnergies(&m_uk[0]); - // compute wall terms size_t loc = m_nsp+3; fill(m_sdot.begin(), m_sdot.end(), 0.0); @@ -210,6 +235,9 @@ void Reactor::evalEqs(doublereal time, doublereal* y, } } + // volume equation + ydot[1] = m_vdot; + const vector_fp& mw = m_thermo->molecularWeights(); const doublereal* Y = m_thermo->massFractions(); @@ -223,52 +251,55 @@ void Reactor::evalEqs(doublereal time, doublereal* y, dYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k] / m_mass; mdot_surf += m_sdot[k] * mw[k]; } - dmdt += mdot_surf; + dmdt += mdot_surf; // mass added to gas phase from surface reations - // compression work and external heat transfer - mcvdTdt += - m_pressure * m_vdot - m_Q; - - for (size_t n = 0; n < m_nsp; n++) { - // heat release from gas phase and surface reations - mcvdTdt -= m_wdot[n] * m_uk[n] * m_vol; - mcvdTdt -= m_sdot[n] * m_uk[n]; + for (size_t k = 0; k < m_nsp; k++) { // dilution by net surface mass flux - dYdt[n] -= Y[n] * mdot_surf / m_mass; + dYdt[k] -= Y[k] * mdot_surf / m_mass; + } + + /* + * Energy equation. + * \f[ + * \dot U = -P\dot V + A \dot q + \dot m_{in} h_{in} + * - \dot m_{out} h. + * \f] + */ + if (m_energy) { + ydot[2] = - m_thermo->pressure() * m_vdot - m_Q; + } else { + ydot[2] = 0.0; } // add terms for open system if (m_open) { + doublereal enthalpy = m_thermo->enthalpy_mass(); + // outlets for (size_t i = 0; i < m_nOutlets; i++) { double mdot_out = m_outlet[i]->massFlowRate(time); dmdt -= mdot_out; // mass flow out of system - mcvdTdt -= mdot_out * m_pressure * m_vol / m_mass; // flow work + if (m_energy) { + ydot[2] -= mdot_out * enthalpy; + } } // inlets for (size_t i = 0; i < m_nInlets; i++) { double mdot_in = m_inlet[i]->massFlowRate(time); dmdt += mdot_in; // mass flow into system - mcvdTdt += m_inlet[i]->enthalpy_mass() * mdot_in; for (size_t n = 0; n < m_nsp; n++) { double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n); // flow of species into system and dilution by other species dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass; - - // In combintion with h_in*mdot_in, flow work plus thermal - // energy carried with the species - mcvdTdt -= m_uk[n] / mw[n] * mdot_spec; + } + if (m_energy) { + ydot[2] += mdot_in * m_inlet[i]->enthalpy_mass(); } } } ydot[0] = dmdt; - ydot[1] = m_vdot; - if (m_energy) { - ydot[2] = mcvdTdt / (m_mass * m_thermo->cv_mass()); - } else { - ydot[2] = 0; - } for (size_t i = 0; i < m_nv; i++) { AssertFinite(ydot[i], "Reactor::evalEqs", @@ -324,7 +355,7 @@ size_t Reactor::componentIndex(const string& nm) const if (nm == "V") { return 1; } - if (nm == "T") { + if (nm == "U") { return 2; } diff --git a/src/zeroD/ReactorFactory.cpp b/src/zeroD/ReactorFactory.cpp index 071fe1241..aa8c14ce2 100644 --- a/src/zeroD/ReactorFactory.cpp +++ b/src/zeroD/ReactorFactory.cpp @@ -9,6 +9,7 @@ #include "cantera/zeroD/Reactor.h" #include "cantera/zeroD/FlowReactor.h" #include "cantera/zeroD/ConstPressureReactor.h" +#include "cantera/zeroD/IdealGasReactor.h" using namespace std; namespace Cantera @@ -17,14 +18,14 @@ namespace Cantera ReactorFactory* ReactorFactory::s_factory = 0; mutex_t ReactorFactory::reactor_mutex; -static int ntypes = 4; +static int ntypes = 5; static string _types[] = {"Reservoir", "Reactor", "ConstPressureReactor", - "FlowReactor" + "FlowReactor", "IdealGasReactor" }; // these constants are defined in ReactorBase.h static int _itypes[] = {ReservoirType, ReactorType, ConstPressureReactorType, - FlowReactorType + FlowReactorType, IdealGasReactorType }; /** @@ -56,6 +57,8 @@ ReactorBase* ReactorFactory::newReactor(int ir) return new FlowReactor(); case ConstPressureReactorType: return new ConstPressureReactor(); + case IdealGasReactorType: + return new IdealGasReactor(); default: throw Cantera::CanteraError("ReactorFactory::newReactor", "unknown reactor type!");