[Reactor] Add optimized IdealGasReactor class
This formulation of the reactor governing equations, with temperature as a state variable, works better for ideal gas mixtures. This way, most of the Jacobian components are derivatives at constant temperature, eliminating the need to recompute the temperature-dependent part of the rate expressions when computing these entries.
This commit is contained in:
Ray Speth
parent
c89b7f1c93
commit
785d4f058e
5 changed files with 367 additions and 16 deletions
50
include/cantera/zeroD/IdealGasReactor.h
Normal file
50
include/cantera/zeroD/IdealGasReactor.h
Normal file
|
|
@ -0,0 +1,50 @@
|
|||
/**
|
||||
* @file IdealGasReactor.h
|
||||
*/
|
||||
|
||||
// Copyright 2001 California Institute of Technology
|
||||
|
||||
#ifndef CT_IDEALGASREACTOR_H
|
||||
#define CT_IDEALGASREACTOR_H
|
||||
|
||||
#include "Reactor.h"
|
||||
#include "cantera/kinetics/Kinetics.h"
|
||||
|
||||
namespace Cantera
|
||||
{
|
||||
|
||||
/**
|
||||
* Class IdealGasReactor is a class for stirred reactors that is specifically
|
||||
* optimized for ideal gases. In this formulation, temperature replaces the
|
||||
* total internal energy as a state variable.
|
||||
*/
|
||||
class IdealGasReactor : public Reactor
|
||||
{
|
||||
public:
|
||||
IdealGasReactor() {}
|
||||
|
||||
virtual int type() const {
|
||||
return IdealGasReactorType;
|
||||
}
|
||||
|
||||
virtual void setThermoMgr(ThermoPhase& thermo);
|
||||
|
||||
virtual void getInitialConditions(doublereal t0, size_t leny,
|
||||
doublereal* y);
|
||||
|
||||
virtual void initialize(doublereal t0 = 0.0);
|
||||
|
||||
virtual void evalEqs(doublereal t, doublereal* y,
|
||||
doublereal* ydot, doublereal* params);
|
||||
|
||||
virtual void updateState(doublereal* y);
|
||||
|
||||
virtual size_t componentIndex(const std::string& nm) const;
|
||||
|
||||
protected:
|
||||
vector_fp m_uk; //!< Species molar internal energies
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
|
|
@ -64,7 +64,7 @@ public:
|
|||
* a pointer to this substance is stored, and as the integration
|
||||
* proceeds, the state of the substance is modified.
|
||||
*/
|
||||
void setThermoMgr(thermo_t& thermo);
|
||||
virtual void setThermoMgr(thermo_t& thermo);
|
||||
|
||||
//! Connect an inlet FlowDevice to this reactor
|
||||
void addInlet(FlowDevice& inlet);
|
||||
|
|
|
|||
|
|
@ -228,6 +228,11 @@ cdef class ConstPressureReactor(Reactor):
|
|||
reactor_type = "ConstPressureReactor"
|
||||
|
||||
|
||||
cdef class IdealGasReactor(Reactor):
|
||||
""" A constant volume, zero-dimensional reactor for ideal gas mixtures. """
|
||||
reactor_type = "IdealGasReactor"
|
||||
|
||||
|
||||
cdef class FlowReactor(Reactor):
|
||||
"""
|
||||
A steady-state plug flow reactor with constant cross sectional area.
|
||||
|
|
|
|||
|
|
@ -7,6 +7,8 @@ from . import utilities
|
|||
|
||||
|
||||
class TestReactor(utilities.CanteraTest):
|
||||
reactorClass = ct.Reactor
|
||||
|
||||
def make_reactors(self, independent=True, n_reactors=2,
|
||||
T1=300, P1=101325, X1='O2:1.0',
|
||||
T2=300, P2=101325, X2='O2:1.0'):
|
||||
|
|
@ -15,7 +17,7 @@ class TestReactor(utilities.CanteraTest):
|
|||
|
||||
self.gas1 = ct.Solution('h2o2.xml')
|
||||
self.gas1.TPX = T1, P1, X1
|
||||
self.r1 = ct.Reactor(self.gas1)
|
||||
self.r1 = self.reactorClass(self.gas1)
|
||||
self.net.add_reactor(self.r1)
|
||||
|
||||
if independent:
|
||||
|
|
@ -25,7 +27,7 @@ class TestReactor(utilities.CanteraTest):
|
|||
|
||||
if n_reactors >= 2:
|
||||
self.gas2.TPX = T2, P2, X2
|
||||
self.r2 = ct.Reactor(self.gas2)
|
||||
self.r2 = self.reactorClass(self.gas2)
|
||||
self.net.add_reactor(self.r2)
|
||||
|
||||
def add_wall(self, **kwargs):
|
||||
|
|
@ -33,7 +35,7 @@ class TestReactor(utilities.CanteraTest):
|
|||
return self.w
|
||||
|
||||
def test_insert(self):
|
||||
R = ct.Reactor()
|
||||
R = self.reactorClass()
|
||||
f1 = lambda r: r.T
|
||||
f2 = lambda r: r.kinetics.net_production_rates
|
||||
self.assertRaises(Exception, f1, R)
|
||||
|
|
@ -470,6 +472,10 @@ class TestReactor(utilities.CanteraTest):
|
|||
self.assertNear(p1a, p1b)
|
||||
self.assertNear(p2a, p2b)
|
||||
|
||||
class TestIdealGasReactor(TestReactor):
|
||||
reactorClass = ct.IdealGasReactor
|
||||
|
||||
|
||||
class TestWellStirredReactorIgnition(utilities.CanteraTest):
|
||||
""" Ignition (or not) of a well-stirred reactor """
|
||||
def setup(self, T0, P0, mdot_fuel, mdot_ox):
|
||||
|
|
@ -486,7 +492,7 @@ class TestWellStirredReactorIgnition(utilities.CanteraTest):
|
|||
|
||||
# reactor, initially filled with N2
|
||||
self.gas.TPX = T0, P0, "N2:1.0"
|
||||
self.combustor = ct.Reactor(self.gas)
|
||||
self.combustor = ct.IdealGasReactor(self.gas)
|
||||
self.combustor.volume = 1.0
|
||||
|
||||
# outlet
|
||||
|
|
@ -587,7 +593,7 @@ class TestConstPressureReactor(utilities.CanteraTest):
|
|||
self.gas1.TPX = T0, P0, X0
|
||||
self.gas2.TPX = T0, P0, X0
|
||||
|
||||
self.r1 = ct.Reactor(self.gas1)
|
||||
self.r1 = ct.IdealGasReactor(self.gas1)
|
||||
self.r2 = ct.ConstPressureReactor(self.gas2)
|
||||
|
||||
self.r1.volume = 0.2
|
||||
|
|
@ -710,11 +716,11 @@ class TestWallKinetics(utilities.CanteraTest):
|
|||
self.solid = ct.Solution('diamond.xml', 'diamond')
|
||||
self.interface = ct.Interface('diamond.xml', 'diamond_100',
|
||||
(self.gas, self.solid))
|
||||
self.r1 = ct.Reactor(self.gas)
|
||||
self.r1 = ct.IdealGasReactor(self.gas)
|
||||
self.r1.volume = 0.01
|
||||
self.net.add_reactor(self.r1)
|
||||
|
||||
self.r2 = ct.Reactor(self.gas)
|
||||
self.r2 = ct.IdealGasReactor(self.gas)
|
||||
self.r2.volume = 0.01
|
||||
self.net.add_reactor(self.r2)
|
||||
|
||||
|
|
@ -809,7 +815,7 @@ class TestReactorSensitivities(utilities.CanteraTest):
|
|||
net = ct.ReactorNet()
|
||||
gas = ct.Solution('gri30.xml')
|
||||
gas.TPX = 1300, 20*101325, 'CO:1.0, H2:0.1, CH4:0.1, H2O:0.5'
|
||||
r1 = ct.Reactor(gas)
|
||||
r1 = ct.IdealGasReactor(gas)
|
||||
net.add_reactor(r1)
|
||||
|
||||
self.assertEqual(net.n_sensitivity_params, 0)
|
||||
|
|
@ -831,14 +837,14 @@ class TestReactorSensitivities(utilities.CanteraTest):
|
|||
solid = ct.Solution('diamond.xml', 'diamond')
|
||||
interface = ct.Interface('diamond.xml', 'diamond_100',
|
||||
(gas1, solid))
|
||||
r1 = ct.Reactor(gas1)
|
||||
r1 = ct.IdealGasReactor(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'
|
||||
r2 = ct.Reactor(gas2)
|
||||
r2 = ct.IdealGasReactor(gas2)
|
||||
net.add_reactor(r2)
|
||||
|
||||
w = ct.Wall(r1, r2)
|
||||
|
|
@ -881,7 +887,7 @@ class TestReactorSensitivities(utilities.CanteraTest):
|
|||
net = ct.ReactorNet()
|
||||
gas.TPX = 900, 101325, 'H2:0.1, OH:1e-7, O2:0.1, AR:1e-5'
|
||||
|
||||
r = ct.Reactor(gas)
|
||||
r = ct.IdealGasReactor(gas)
|
||||
net.add_reactor(r)
|
||||
return r, net
|
||||
|
||||
|
|
@ -920,11 +926,11 @@ class TestReactorSensitivities(utilities.CanteraTest):
|
|||
net = ct.ReactorNet()
|
||||
gas1 = ct.Solution('h2o2.xml')
|
||||
gas1.TPX = 900, 101325, 'H2:0.1, OH:1e-7, O2:0.1, AR:1e-5'
|
||||
rA = ct.Reactor(gas1)
|
||||
rA = ct.IdealGasReactor(gas1)
|
||||
|
||||
gas2 = ct.Solution('h2o2.xml')
|
||||
gas2.TPX = 920, 101325, 'H2:0.1, OH:1e-7, O2:0.1, AR:0.5'
|
||||
rB = ct.Reactor(gas2)
|
||||
rB = ct.IdealGasReactor(gas2)
|
||||
if reverse:
|
||||
net.add_reactor(rB)
|
||||
net.add_reactor(rA)
|
||||
|
|
@ -986,8 +992,8 @@ class TestReactorSensitivities(utilities.CanteraTest):
|
|||
gas1.TPX = 1200, 1e3, 'H:0.002, H2:1, CH4:0.01, CH3:0.0002'
|
||||
gas2.TPX = 900, 101325, 'H2:0.1, OH:1e-7, O2:0.1, AR:1e-5'
|
||||
net = ct.ReactorNet()
|
||||
rA = ct.Reactor(gas1)
|
||||
rB = ct.Reactor(gas2)
|
||||
rA = ct.IdealGasReactor(gas1)
|
||||
rB = ct.IdealGasReactor(gas2)
|
||||
|
||||
if order % 2 == 0:
|
||||
wA = ct.Wall(rA, rB)
|
||||
|
|
|
|||
290
src/zeroD/IdealGasReactor.cpp
Normal file
290
src/zeroD/IdealGasReactor.cpp
Normal file
|
|
@ -0,0 +1,290 @@
|
|||
/**
|
||||
* @file IdealGasReactor.cpp A zero-dimensional reactor
|
||||
*/
|
||||
|
||||
#include "cantera/zeroD/IdealGasReactor.h"
|
||||
#include "cantera/zeroD/FlowDevice.h"
|
||||
#include "cantera/zeroD/Wall.h"
|
||||
#include "cantera/kinetics/InterfaceKinetics.h"
|
||||
#include "cantera/thermo/SurfPhase.h"
|
||||
#include "cantera/zeroD/ReactorNet.h"
|
||||
|
||||
#include <cfloat>
|
||||
|
||||
using namespace std;
|
||||
|
||||
namespace Cantera
|
||||
{
|
||||
|
||||
void IdealGasReactor::setThermoMgr(ThermoPhase& thermo)
|
||||
{
|
||||
//! @TODO: Add a method to ThermoPhase that indicates whether a given
|
||||
//! subclass is compatible with this reactor model
|
||||
if (thermo.eosType() != cIdealGas) {
|
||||
throw CanteraError("IdealGasReactor::setThermoMgr",
|
||||
"Incompatible phase type provided");
|
||||
}
|
||||
Reactor::setThermoMgr(thermo);
|
||||
}
|
||||
|
||||
void IdealGasReactor::getInitialConditions(double t0, size_t leny, double* y)
|
||||
{
|
||||
m_init = true;
|
||||
if (m_thermo == 0) {
|
||||
cout << "Error: reactor is empty." << endl;
|
||||
return;
|
||||
}
|
||||
m_thermo->restoreState(m_state);
|
||||
|
||||
// set the first component to the total mass
|
||||
m_mass = m_thermo->density() * m_vol;
|
||||
y[0] = m_mass;
|
||||
|
||||
// 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 components y+3 ... y+K+2 to the mass fractions of each species
|
||||
m_thermo->getMassFractions(y+3);
|
||||
|
||||
// set the remaining components to the surface species
|
||||
// coverages on the walls
|
||||
size_t loc = m_nsp + 3;
|
||||
SurfPhase* surf;
|
||||
for (size_t m = 0; m < m_nwalls; m++) {
|
||||
surf = m_wall[m]->surface(m_lr[m]);
|
||||
if (surf) {
|
||||
m_wall[m]->getCoverages(m_lr[m], y + loc);
|
||||
loc += surf->nSpecies();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void IdealGasReactor::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])) {
|
||||
m_nv += m_wall[w]->surface(m_lr[w])->nSpecies();
|
||||
}
|
||||
|
||||
m_enthalpy = m_thermo->enthalpy_mass();
|
||||
m_pressure = m_thermo->pressure();
|
||||
m_intEnergy = m_thermo->intEnergy_mass();
|
||||
|
||||
size_t nt = 0, maxnt = 0;
|
||||
for (size_t m = 0; m < m_nwalls; m++) {
|
||||
m_wall[m]->initialize();
|
||||
if (m_wall[m]->kinetics(m_lr[m])) {
|
||||
nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies();
|
||||
if (nt > maxnt) {
|
||||
maxnt = nt;
|
||||
}
|
||||
if (m_wall[m]->kinetics(m_lr[m])) {
|
||||
if (&m_kin->thermo(0) !=
|
||||
&m_wall[m]->kinetics(m_lr[m])->thermo(0)) {
|
||||
throw CanteraError("IdealGasReactor::initialize",
|
||||
"First phase of all kinetics managers must be"
|
||||
" the gas.");
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
m_work.resize(maxnt);
|
||||
std::sort(m_pnum.begin(), m_pnum.end());
|
||||
m_init = true;
|
||||
}
|
||||
|
||||
void IdealGasReactor::updateState(doublereal* y)
|
||||
{
|
||||
for (size_t i = 0; i < m_nv; i++) {
|
||||
AssertFinite(y[i], "IdealGasReactor::updateState",
|
||||
"y[" + int2str(i) + "] is not finite");
|
||||
}
|
||||
|
||||
// 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.
|
||||
m_mass = y[0];
|
||||
m_vol = y[1];
|
||||
|
||||
m_thermo->setMassFractions_NoNorm(y+3);
|
||||
m_thermo->setState_TR(y[2], m_mass / m_vol);
|
||||
|
||||
size_t loc = m_nsp + 3;
|
||||
SurfPhase* surf;
|
||||
for (size_t m = 0; m < m_nwalls; m++) {
|
||||
surf = m_wall[m]->surface(m_lr[m]);
|
||||
if (surf) {
|
||||
m_wall[m]->setCoverages(m_lr[m], y+loc);
|
||||
loc += surf->nSpecies();
|
||||
}
|
||||
}
|
||||
|
||||
// save parameters needed by other connected reactors
|
||||
m_enthalpy = m_thermo->enthalpy_mass();
|
||||
m_pressure = m_thermo->pressure();
|
||||
m_intEnergy = m_thermo->intEnergy_mass();
|
||||
m_thermo->saveState(m_state);
|
||||
}
|
||||
|
||||
void IdealGasReactor::evalEqs(doublereal time, doublereal* y,
|
||||
doublereal* ydot, doublereal* params)
|
||||
{
|
||||
m_thermo->restoreState(m_state);
|
||||
|
||||
// process sensitivity parameters
|
||||
if (params) {
|
||||
size_t npar = m_pnum.size();
|
||||
for (size_t n = 0; n < npar; n++) {
|
||||
double mult = m_kin->multiplier(m_pnum[n]);
|
||||
m_kin->setMultiplier(m_pnum[n], mult*params[n]);
|
||||
}
|
||||
size_t ploc = npar;
|
||||
for (size_t m = 0; m < m_nwalls; m++) {
|
||||
if (m_nsens_wall[m] > 0) {
|
||||
m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc);
|
||||
ploc += m_nsens_wall[m];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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);
|
||||
for (size_t i = 0; i < m_nwalls; i++) {
|
||||
int lr = 1 - 2*m_lr[i];
|
||||
double vdot = lr*m_wall[i]->vdot(time);
|
||||
m_vdot += vdot;
|
||||
m_Q += lr*m_wall[i]->Q(time);
|
||||
Kinetics* kin = m_wall[i]->kinetics(m_lr[i]);
|
||||
SurfPhase* surf = m_wall[i]->surface(m_lr[i]);
|
||||
if (surf && kin) {
|
||||
double rs0 = 1.0/surf->siteDensity();
|
||||
size_t nk = surf->nSpecies();
|
||||
double sum = 0.0;
|
||||
surf->setTemperature(m_state[0]);
|
||||
m_wall[i]->syncCoverages(m_lr[i]);
|
||||
kin->getNetProductionRates(DATA_PTR(m_work));
|
||||
size_t ns = kin->surfacePhaseIndex();
|
||||
size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
|
||||
for (size_t k = 1; k < nk; k++) {
|
||||
ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
|
||||
sum -= ydot[loc + k];
|
||||
}
|
||||
ydot[loc] = sum;
|
||||
loc += nk;
|
||||
|
||||
double wallarea = m_wall[i]->area();
|
||||
for (size_t k = 0; k < m_nsp; k++) {
|
||||
m_sdot[k] += m_work[k]*wallarea;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
const vector_fp& mw = m_thermo->molecularWeights();
|
||||
const doublereal* Y = m_thermo->massFractions();
|
||||
|
||||
if (m_chem) {
|
||||
m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
|
||||
}
|
||||
|
||||
double mdot_surf = 0.0; // net mass flux from surfaces
|
||||
for (size_t k = 0; k < m_nsp; k++) {
|
||||
// production in gas phase and from surfaces
|
||||
dYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k] / m_mass;
|
||||
mdot_surf += m_sdot[k] * mw[k];
|
||||
}
|
||||
dmdt += mdot_surf;
|
||||
|
||||
// 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];
|
||||
// dilution by net surface mass flux
|
||||
dYdt[n] -= Y[n] * mdot_surf / m_mass;
|
||||
}
|
||||
|
||||
// add terms for open system
|
||||
if (m_open) {
|
||||
// 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
|
||||
}
|
||||
|
||||
// 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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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], "IdealGasReactor::evalEqs",
|
||||
"ydot[" + int2str(i) + "] is not finite");
|
||||
}
|
||||
|
||||
// reset sensitivity parameters
|
||||
if (params) {
|
||||
size_t npar = m_pnum.size();
|
||||
for (size_t n = 0; n < npar; n++) {
|
||||
double mult = m_kin->multiplier(m_pnum[n]);
|
||||
m_kin->setMultiplier(m_pnum[n], mult/params[n]);
|
||||
}
|
||||
size_t ploc = npar;
|
||||
for (size_t m = 0; m < m_nwalls; m++) {
|
||||
if (m_nsens_wall[m] > 0) {
|
||||
m_wall[m]->resetSensitivityParameters(m_lr[m]);
|
||||
ploc += m_nsens_wall[m];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
size_t IdealGasReactor::componentIndex(const string& nm) const
|
||||
{
|
||||
if (nm == "T") {
|
||||
return 2;
|
||||
} else {
|
||||
return Reactor::componentIndex(nm);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
Loading…
Add table
Reference in a new issue