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

[Reactor] Use correct phase state after mass flow rate evaluation

Browse source 
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.
This commit is contained in:
Ray Speth 2019-06-14 17:39:13 -04:00
parent edcc9c59fd
commit a247d0f4eb
8 changed files with 81 additions and 36 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

6
include/cantera/zeroD/Reactor.h
View file

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

7
include/cantera/zeroD/ReactorBase.h
View file

@ -239,6 +239,13 @@ protected:
doublereal m_pressure;
vector_fp m_state;
std::vector<FlowDevice*> 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<WallBase*> m_wall;
std::vector<ReactorSurface*> m_surfaces;
vector_int m_lr;

24
interfaces/cython/cantera/test/test_reactor.py
View file

@ -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:
"""

19
src/zeroD/ConstPressureReactor.cpp
View file

@ -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;

13
src/zeroD/IdealGasConstPressureReactor.cpp
View file

@ -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;
}
}

19
src/zeroD/IdealGasReactor.cpp
View file

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

27
src/zeroD/Reactor.cpp
View file

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

2
src/zeroD/ReactorBase.cpp
View file

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