Cantera  2.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
IdealGasConstPressureReactor.cpp
1 /**
2  * @file ConstPressureReactor.cpp A constant pressure zero-dimensional
3  * reactor
4  */
5 
6 // Copyright 2001 California Institute of Technology
7 
8 #include "cantera/zeroD/IdealGasConstPressureReactor.h"
10 
11 using namespace std;
12 
13 namespace Cantera
14 {
15 
16 void IdealGasConstPressureReactor::setThermoMgr(ThermoPhase& thermo)
17 {
18  //! @TODO: Add a method to ThermoPhase that indicates whether a given
19  //! subclass is compatible with this reactor model
20  if (thermo.eosType() != cIdealGas) {
21  throw CanteraError("IdealGasReactor::setThermoMgr",
22  "Incompatible phase type provided");
23  }
24  Reactor::setThermoMgr(thermo);
25 }
26 
27 
28 void IdealGasConstPressureReactor::getInitialConditions(double t0, size_t leny,
29  double* y)
30 {
31  if (m_thermo == 0) {
32  throw CanteraError("getInitialConditions",
33  "Error: reactor is empty.");
34  }
35  m_thermo->restoreState(m_state);
36 
37  // set the first component to the total mass
38  y[0] = m_thermo->density() * m_vol;
39 
40  // set the second component to the temperature
41  y[1] = m_thermo->temperature();
42 
43  // set components y+2 ... y+K+1 to the mass fractions Y_k of each species
44  m_thermo->getMassFractions(y+2);
45 
46  // set the remaining components to the surface species
47  // coverages on the walls
48  getSurfaceInitialConditions(y + m_nsp + 2);
49 }
50 
51 void IdealGasConstPressureReactor::initialize(doublereal t0)
52 {
53  ConstPressureReactor::initialize(t0);
54  m_hk.resize(m_nsp, 0.0);
55 }
56 
57 void IdealGasConstPressureReactor::updateState(doublereal* y)
58 {
59  // The components of y are [0] the total mass, [1] the temperature,
60  // [2...K+2) are the mass fractions of each species, and [K+2...] are the
61  // coverages of surface species on each wall.
62  m_mass = y[0];
63  m_thermo->setMassFractions_NoNorm(y+2);
64  m_thermo->setState_TP(y[1], m_pressure);
65  m_vol = m_mass / m_thermo->density();
66  updateSurfaceState(y + m_nsp + 2);
67 
68  // save parameters needed by other connected reactors
69  m_enthalpy = m_thermo->enthalpy_mass();
70  m_intEnergy = m_thermo->intEnergy_mass();
71  m_thermo->saveState(m_state);
72 }
73 
74 void IdealGasConstPressureReactor::evalEqs(doublereal time, doublereal* y,
75  doublereal* ydot, doublereal* params)
76 {
77  double dmdt = 0.0; // dm/dt (gas phase)
78  double mcpdTdt = 0.0; // m * c_p * dT/dt
79  double* dYdt = ydot + 2;
80 
81  m_thermo->restoreState(m_state);
82  applySensitivity(params);
83  evalWalls(time);
84  double mdot_surf = evalSurfaces(time, ydot + m_nsp + 2);
85  dmdt += mdot_surf;
86 
87  m_thermo->getPartialMolarEnthalpies(&m_hk[0]);
88  const vector_fp& mw = m_thermo->molecularWeights();
89  const doublereal* Y = m_thermo->massFractions();
90 
91  if (m_chem) {
92  m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
93  }
94 
95  // external heat transfer
96  mcpdTdt -= m_Q;
97 
98  for (size_t n = 0; n < m_nsp; n++) {
99  // heat release from gas phase and surface reations
100  mcpdTdt -= m_wdot[n] * m_hk[n] * m_vol;
101  mcpdTdt -= m_sdot[n] * m_hk[n];
102  // production in gas phase and from surfaces
103  dYdt[n] = (m_wdot[n] * m_vol + m_sdot[n]) * mw[n] / m_mass;
104  // dilution by net surface mass flux
105  dYdt[n] -= Y[n] * mdot_surf / m_mass;
106  }
107 
108  // add terms for outlets
109  for (size_t i = 0; i < m_outlet.size(); i++) {
110  dmdt -= m_outlet[i]->massFlowRate(time); // mass flow out of system
111  }
112 
113  // add terms for inlets
114  for (size_t i = 0; i < m_inlet.size(); i++) {
115  double mdot_in = m_inlet[i]->massFlowRate(time);
116  dmdt += mdot_in; // mass flow into system
117  mcpdTdt += m_inlet[i]->enthalpy_mass() * mdot_in;
118  for (size_t n = 0; n < m_nsp; n++) {
119  double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
120  // flow of species into system and dilution by other species
121  dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass;
122  mcpdTdt -= m_hk[n] / mw[n] * mdot_spec;
123  }
124  }
125 
126  ydot[0] = dmdt;
127  if (m_energy) {
128  ydot[1] = mcpdTdt / (m_mass * m_thermo->cp_mass());
129  } else {
130  ydot[1] = 0.0;
131  }
132 
133  resetSensitivity(params);
134 }
135 
136 size_t IdealGasConstPressureReactor::componentIndex(const string& nm) const
137 {
138  size_t k = speciesIndex(nm);
139  if (k != npos) {
140  return k + 2;
141  } else if (nm == "m" || nm == "mass") {
142  return 0;
143  } else if (nm == "T" || nm == "temperature") {
144  return 1;
145  } else {
146  return npos;
147  }
148 }
149 
150 }
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:97
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:99
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
virtual int eosType() const
Equation of state type flag.
Definition: ThermoPhase.h:142
const int cIdealGas
Equation of state types:
Definition: mix_defs.h:37