From a247d0f4ebc8ea8294367434bd172914dff49209 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 14 Jun 2019 17:39:13 -0400 Subject: [PATCH] [Reactor] Use correct phase state after mass flow rate evaluation A user-defined mass flow rate function can modify the ThermoPhase object used by a reactor, for example if it depends on calculating some property of a different reactor. To make sure that the reactor governing equations are evaluated correctly, the ThermoPhase state needs to be set after all user-defined functions have been called. --- include/cantera/zeroD/Reactor.h | 6 +++++ include/cantera/zeroD/ReactorBase.h | 7 +++++ .../cython/cantera/test/test_reactor.py | 24 ++++++++++++++--- src/zeroD/ConstPressureReactor.cpp | 19 ++++++------- src/zeroD/IdealGasConstPressureReactor.cpp | 13 ++++----- src/zeroD/IdealGasReactor.cpp | 19 ++++++------- src/zeroD/Reactor.cpp | 27 ++++++++++++------- src/zeroD/ReactorBase.cpp | 2 ++ 8 files changed, 81 insertions(+), 36 deletions(-) diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 54ea826c0..d02a672c8 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -164,6 +164,12 @@ protected: //! @param t the current time virtual void evalWalls(double t); + //! Evaluate inlet and outlet mass flow rates. This is called in evalEqs() + //! before setting the state of #m_thermo, since calling the mass flow rate + //! functions may modify ThermoPhase objects that are shared with other + //! reactors. + virtual void evalFlowDevices(double t); + //! Evaluate terms related to surface reactions. Calculates #m_sdot and rate //! of change in surface species coverages. //! @param t the current time diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 054079c0a..5489dd828 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -239,6 +239,13 @@ protected: doublereal m_pressure; vector_fp m_state; std::vector m_inlet, m_outlet; + + //! Temporary storage for mass flow rates from each inlet FlowDevice + vector_fp m_mdot_in; + + //! Temporary storage for mass flow rates from each outlet FlowDevice + vector_fp m_mdot_out; + std::vector m_wall; std::vector m_surfaces; vector_int m_lr; diff --git a/interfaces/cython/cantera/test/test_reactor.py b/interfaces/cython/cantera/test/test_reactor.py index b39c77c3c..670e45ad4 100644 --- a/interfaces/cython/cantera/test/test_reactor.py +++ b/interfaces/cython/cantera/test/test_reactor.py @@ -1425,9 +1425,9 @@ class CombustorTestImplementation: # create and install the mass flow controllers. Controllers # m1 and m2 provide constant mass flow rates, and m3 provides # a short Gaussian pulse only to ignite the mixture - m1 = ct.MassFlowController(fuel_in, self.combustor, mdot=fuel_mdot) - m2 = ct.MassFlowController(oxidizer_in, self.combustor, mdot=oxidizer_mdot) - m3 = ct.MassFlowController(self.igniter, self.combustor, mdot=igniter_mdot) + self.m1 = ct.MassFlowController(fuel_in, self.combustor, mdot=fuel_mdot) + self.m2 = ct.MassFlowController(oxidizer_in, self.combustor, mdot=oxidizer_mdot) + self.m3 = ct.MassFlowController(self.igniter, self.combustor, mdot=igniter_mdot) # put a valve on the exhaust line to regulate the pressure self.v = ct.Valve(self.combustor, self.exhaust, K=1.0) @@ -1463,6 +1463,24 @@ class CombustorTestImplementation: rtol=1e-6, atol=1e-12) self.assertFalse(bad, bad) + def test_invasive_mdot_function(self): + def igniter_mdot(t, t0=0.1, fwhm=0.05, amplitude=0.1): + # Querying properties of the igniter changes the state of the + # underlying ThermoPhase object, but shouldn't affect the + # integration + self.igniter.density + return amplitude * math.exp(-(t-t0)**2 * 4 * math.log(2) / fwhm**2) + self.m3.set_mass_flow_rate(igniter_mdot) + + self.data = [] + for t in np.linspace(0, 0.25, 101)[1:]: + self.sim.advance(t) + self.data.append([t, self.combustor.T] + + list(self.combustor.thermo.X)) + + bad = utilities.compareProfiles(self.referenceFile, self.data, + rtol=1e-6, atol=1e-12) + self.assertFalse(bad, bad) class WallTestImplementation: """ diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp index 0c3eb4823..75af74068 100644 --- a/src/zeroD/ConstPressureReactor.cpp +++ b/src/zeroD/ConstPressureReactor.cpp @@ -67,9 +67,12 @@ void ConstPressureReactor::evalEqs(doublereal time, doublereal* y, { double dmdt = 0.0; // dm/dt (gas phase) double* dYdt = ydot + 2; - m_thermo->restoreState(m_state); - applySensitivity(params); + + evalFlowDevices(time); evalWalls(time); + applySensitivity(params); + + m_thermo->restoreState(m_state); double mdot_surf = evalSurfaces(time, ydot + m_nsp + 2); dmdt += mdot_surf; @@ -92,21 +95,19 @@ void ConstPressureReactor::evalEqs(doublereal time, doublereal* y, // add terms for outlets for (size_t i = 0; i < m_outlet.size(); i++) { - double mdot_out = m_outlet[i]->massFlowRate(time); // mass flow out of system - dmdt -= mdot_out; - dHdt -= mdot_out * m_enthalpy; + dmdt -= m_mdot_out[i]; + dHdt -= m_mdot_out[i] * m_enthalpy; } // add terms for inlets for (size_t i = 0; i < m_inlet.size(); i++) { - double mdot_in = m_inlet[i]->massFlowRate(time); - dmdt += mdot_in; // mass flow into system + dmdt += m_mdot_in[i]; // mass flow into system 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; + dYdt[n] += (mdot_spec - m_mdot_in[i] * Y[n]) / m_mass; } - dHdt += mdot_in * m_inlet[i]->enthalpy_mass(); + dHdt += m_mdot_in[i] * m_inlet[i]->enthalpy_mass(); } ydot[0] = dmdt; diff --git a/src/zeroD/IdealGasConstPressureReactor.cpp b/src/zeroD/IdealGasConstPressureReactor.cpp index 772a6a733..44b893aca 100644 --- a/src/zeroD/IdealGasConstPressureReactor.cpp +++ b/src/zeroD/IdealGasConstPressureReactor.cpp @@ -74,9 +74,11 @@ void IdealGasConstPressureReactor::evalEqs(doublereal time, doublereal* y, double mcpdTdt = 0.0; // m * c_p * dT/dt double* dYdt = ydot + 2; - m_thermo->restoreState(m_state); + evalFlowDevices(time); applySensitivity(params); evalWalls(time); + + m_thermo->restoreState(m_state); double mdot_surf = evalSurfaces(time, ydot + m_nsp + 2); dmdt += mdot_surf; @@ -103,18 +105,17 @@ void IdealGasConstPressureReactor::evalEqs(doublereal time, doublereal* y, // add terms for outlets for (size_t i = 0; i < m_outlet.size(); i++) { - dmdt -= m_outlet[i]->massFlowRate(time); // mass flow out of system + dmdt -= m_mdot_out[i]; // mass flow out of system } // add terms for inlets for (size_t i = 0; i < m_inlet.size(); i++) { - double mdot_in = m_inlet[i]->massFlowRate(time); - dmdt += mdot_in; // mass flow into system - mcpdTdt += m_inlet[i]->enthalpy_mass() * mdot_in; + dmdt += m_mdot_in[i]; // mass flow into system + mcpdTdt += m_inlet[i]->enthalpy_mass() * m_mdot_in[i]; 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; + dYdt[n] += (mdot_spec - m_mdot_in[i] * Y[n]) / m_mass; mcpdTdt -= m_hk[n] / mw[n] * mdot_spec; } } diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp index 05b407d22..05316091b 100644 --- a/src/zeroD/IdealGasReactor.cpp +++ b/src/zeroD/IdealGasReactor.cpp @@ -80,8 +80,10 @@ void IdealGasReactor::evalEqs(doublereal time, doublereal* y, double mcvdTdt = 0.0; // m * c_v * dT/dt double* dYdt = ydot + 3; - m_thermo->restoreState(m_state); + evalFlowDevices(time); + evalWalls(time); applySensitivity(params); + m_thermo->restoreState(m_state); m_thermo->getPartialMolarIntEnergies(&m_uk[0]); const vector_fp& mw = m_thermo->molecularWeights(); const doublereal* Y = m_thermo->massFractions(); @@ -90,7 +92,6 @@ void IdealGasReactor::evalEqs(doublereal time, doublereal* y, m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" } - evalWalls(time); double mdot_surf = evalSurfaces(time, ydot + m_nsp + 3); dmdt += mdot_surf; @@ -109,20 +110,20 @@ void IdealGasReactor::evalEqs(doublereal time, doublereal* y, // add terms for outlets for (size_t i = 0; i < m_outlet.size(); 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 + // double mdot_out = m_outlet[i]->massFlowRate(time); + dmdt -= m_mdot_out[i]; // mass flow out of system + mcvdTdt -= m_mdot_out[i] * m_pressure * m_vol / m_mass; // flow work } // add terms for inlets for (size_t i = 0; i < m_inlet.size(); i++) { - double mdot_in = m_inlet[i]->massFlowRate(time); - dmdt += mdot_in; // mass flow into system - mcvdTdt += m_inlet[i]->enthalpy_mass() * mdot_in; + // double mdot_in = m_inlet[i]->massFlowRate(time); + dmdt += m_mdot_in[i]; // mass flow into system + mcvdTdt += m_inlet[i]->enthalpy_mass() * m_mdot_in[i]; 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; + dYdt[n] += (mdot_spec - m_mdot_in[i] * Y[n]) / m_mass; // In combination with h_in*mdot_in, flow work plus thermal // energy carried with the species diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 924126de2..80c284b28 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -189,9 +189,10 @@ void Reactor::evalEqs(doublereal time, doublereal* y, double dmdt = 0.0; // dm/dt (gas phase) double* dYdt = ydot + 3; - m_thermo->restoreState(m_state); - applySensitivity(params); + evalFlowDevices(time); evalWalls(time); + applySensitivity(params); + m_thermo->restoreState(m_state); double mdot_surf = evalSurfaces(time, ydot + m_nsp + 3); dmdt += mdot_surf; // mass added to gas phase from surface reactions @@ -224,24 +225,22 @@ void Reactor::evalEqs(doublereal time, doublereal* y, // add terms for outlets for (size_t i = 0; i < m_outlet.size(); i++) { - double mdot_out = m_outlet[i]->massFlowRate(time); - dmdt -= mdot_out; // mass flow out of system + dmdt -= m_mdot_out[i]; // mass flow out of system if (m_energy) { - ydot[2] -= mdot_out * m_enthalpy; + ydot[2] -= m_mdot_out[i] * m_enthalpy; } } // add terms for inlets for (size_t i = 0; i < m_inlet.size(); i++) { - double mdot_in = m_inlet[i]->massFlowRate(time); - dmdt += mdot_in; // mass flow into system + dmdt += m_mdot_in[i]; // mass flow into system 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; + dYdt[n] += (mdot_spec - m_mdot_in[i] * Y[n]) / m_mass; } if (m_energy) { - ydot[2] += mdot_in * m_inlet[i]->enthalpy_mass(); + ydot[2] += m_mdot_in[i] * m_inlet[i]->enthalpy_mass(); } } @@ -260,6 +259,16 @@ void Reactor::evalWalls(double t) } } +void Reactor::evalFlowDevices(double t) +{ + for (size_t i = 0; i < m_outlet.size(); i++) { + m_mdot_out[i] = m_outlet[i]->massFlowRate(t); + } + for (size_t i = 0; i < m_inlet.size(); i++) { + m_mdot_in[i] = m_inlet[i]->massFlowRate(t); + } +} + double Reactor::evalSurfaces(double t, double* ydot) { const vector_fp& mw = m_thermo->molecularWeights(); diff --git a/src/zeroD/ReactorBase.cpp b/src/zeroD/ReactorBase.cpp index edcb7e153..566250e45 100644 --- a/src/zeroD/ReactorBase.cpp +++ b/src/zeroD/ReactorBase.cpp @@ -47,11 +47,13 @@ void ReactorBase::syncState() void ReactorBase::addInlet(FlowDevice& inlet) { m_inlet.push_back(&inlet); + m_mdot_in.push_back(0.0); } void ReactorBase::addOutlet(FlowDevice& outlet) { m_outlet.push_back(&outlet); + m_mdot_out.push_back(0.0); } void ReactorBase::addWall(WallBase& w, int lr)