Cantera  2.4.0
IdealGasReactor.cpp
Go to the documentation of this file.
1 //! @file IdealGasReactor.cpp A zero-dimensional reactor
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at http://www.cantera.org/license.txt for license and copyright information.
5 
8 #include "cantera/zeroD/Wall.h"
9 
10 using namespace std;
11 
12 namespace Cantera
13 {
14 
15 void IdealGasReactor::setThermoMgr(ThermoPhase& thermo)
16 {
17  //! @TODO: Add a method to ThermoPhase that indicates whether a given
18  //! subclass is compatible with this reactor model
19  if (thermo.type() != "IdealGas") {
20  throw CanteraError("IdealGasReactor::setThermoMgr",
21  "Incompatible phase type provided");
22  }
23  Reactor::setThermoMgr(thermo);
24 }
25 
26 void IdealGasReactor::getState(double* y)
27 {
28  if (m_thermo == 0) {
29  throw CanteraError("getState",
30  "Error: reactor is empty.");
31  }
32  m_thermo->restoreState(m_state);
33 
34  // set the first component to the total mass
35  m_mass = m_thermo->density() * m_vol;
36  y[0] = m_mass;
37 
38  // set the second component to the total volume
39  y[1] = m_vol;
40 
41  // Set the third component to the temperature
42  y[2] = m_thermo->temperature();
43 
44  // set components y+3 ... y+K+2 to the mass fractions of each species
45  m_thermo->getMassFractions(y+3);
46 
47  // set the remaining components to the surface species
48  // coverages on the walls
49  getSurfaceInitialConditions(y + m_nsp + 3);
50 }
51 
52 void IdealGasReactor::initialize(doublereal t0)
53 {
54  Reactor::initialize(t0);
55  m_uk.resize(m_nsp, 0.0);
56 }
57 
58 void IdealGasReactor::updateState(doublereal* y)
59 {
60  // The components of y are [0] the total mass, [1] the total volume,
61  // [2] the temperature, [3...K+3] are the mass fractions of each species,
62  // and [K+3...] are the coverages of surface species on each wall.
63  m_mass = y[0];
64  m_vol = y[1];
65  m_thermo->setMassFractions_NoNorm(y+3);
66  m_thermo->setState_TR(y[2], m_mass / m_vol);
67  updateSurfaceState(y + m_nsp + 3);
68 
69  // save parameters needed by other connected reactors
70  m_enthalpy = m_thermo->enthalpy_mass();
71  m_pressure = m_thermo->pressure();
72  m_intEnergy = m_thermo->intEnergy_mass();
73  m_thermo->saveState(m_state);
74 }
75 
76 void IdealGasReactor::evalEqs(doublereal time, doublereal* y,
77  doublereal* ydot, doublereal* params)
78 {
79  double dmdt = 0.0; // dm/dt (gas phase)
80  double mcvdTdt = 0.0; // m * c_v * dT/dt
81  double* dYdt = ydot + 3;
82 
83  m_thermo->restoreState(m_state);
84  applySensitivity(params);
85  m_thermo->getPartialMolarIntEnergies(&m_uk[0]);
86  const vector_fp& mw = m_thermo->molecularWeights();
87  const doublereal* Y = m_thermo->massFractions();
88 
89  if (m_chem) {
90  m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
91  }
92 
93  evalWalls(time);
94  double mdot_surf = evalSurfaces(time, ydot + m_nsp + 3);
95  dmdt += mdot_surf;
96 
97  // compression work and external heat transfer
98  mcvdTdt += - m_pressure * m_vdot - m_Q;
99 
100  for (size_t n = 0; n < m_nsp; n++) {
101  // heat release from gas phase and surface reactions
102  mcvdTdt -= m_wdot[n] * m_uk[n] * m_vol;
103  mcvdTdt -= m_sdot[n] * m_uk[n];
104  // production in gas phase and from surfaces
105  dYdt[n] = (m_wdot[n] * m_vol + m_sdot[n]) * mw[n] / m_mass;
106  // dilution by net surface mass flux
107  dYdt[n] -= Y[n] * mdot_surf / m_mass;
108  }
109 
110  // add terms for outlets
111  for (size_t i = 0; i < m_outlet.size(); i++) {
112  double mdot_out = m_outlet[i]->massFlowRate(time);
113  dmdt -= mdot_out; // mass flow out of system
114  mcvdTdt -= mdot_out * m_pressure * m_vol / m_mass; // flow work
115  }
116 
117  // add terms for inlets
118  for (size_t i = 0; i < m_inlet.size(); i++) {
119  double mdot_in = m_inlet[i]->massFlowRate(time);
120  dmdt += mdot_in; // mass flow into system
121  mcvdTdt += m_inlet[i]->enthalpy_mass() * mdot_in;
122  for (size_t n = 0; n < m_nsp; n++) {
123  double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
124  // flow of species into system and dilution by other species
125  dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass;
126 
127  // In combination with h_in*mdot_in, flow work plus thermal
128  // energy carried with the species
129  mcvdTdt -= m_uk[n] / mw[n] * mdot_spec;
130  }
131  }
132 
133  ydot[0] = dmdt;
134  ydot[1] = m_vdot;
135  if (m_energy) {
136  ydot[2] = mcvdTdt / (m_mass * m_thermo->cv_mass());
137  } else {
138  ydot[2] = 0;
139  }
140 
141  resetSensitivity(params);
142 }
143 
144 size_t IdealGasReactor::componentIndex(const string& nm) const
145 {
146  size_t k = speciesIndex(nm);
147  if (k != npos) {
148  return k + 3;
149  } else if (nm == "mass") {
150  return 0;
151  } else if (nm == "volume") {
152  return 1;
153  } else if (nm == "temperature") {
154  return 2;
155  } else {
156  return npos;
157  }
158 }
159 
160 std::string IdealGasReactor::componentName(size_t k) {
161  if (k == 2) {
162  return "temperature";
163  } else {
164  return Reactor::componentName(k);
165  }
166 }
167 
168 
169 }
Header file for class Wall.
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
STL namespace.
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:93
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
virtual std::string type() const
String indicating the thermodynamic model implemented.
Definition: ThermoPhase.h:108
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8