Cantera  2.2.1
ConstPressureReactor.cpp
Go to the documentation of this file.
1 /**
2  * @file ConstPressureReactor.cpp A constant pressure zero-dimensional
3  * reactor
4  */
5
6 // Copyright 2001 California Institute of Technology
7
10
11 using namespace std;
12
13 namespace Cantera
14 {
15
16 void ConstPressureReactor::getInitialConditions(double t0, size_t leny, double* y)
17 {
18  if (m_thermo == 0) {
19  throw CanteraError("getInitialConditions",
20  "Error: reactor is empty.");
21  }
22  m_thermo->restoreState(m_state);
23
24  // set the first component to the total mass
25  y[0] = m_thermo->density() * m_vol;
26
27  // set the second component to the total enthalpy
28  y[1] = m_thermo->enthalpy_mass() * m_thermo->density() * m_vol;
29
30  // set components y+2 ... y+K+1 to the mass fractions Y_k of each species
31  m_thermo->getMassFractions(y+2);
32
33  // set the remaining components to the surface species
34  // coverages on the walls
35  getSurfaceInitialConditions(y + m_nsp + 2);
36 }
37
38 void ConstPressureReactor::initialize(doublereal t0)
39 {
40  Reactor::initialize(t0);
41  m_nv -= 1; // Constant pressure reactor has one fewer state variable
42 }
43
45 {
46  // The components of y are [0] the total mass, [1] the total enthalpy,
47  // [2...K+2) are the mass fractions of each species, and [K+2...] are the
48  // coverages of surface species on each wall.
49  m_mass = y[0];
50  m_thermo->setMassFractions_NoNorm(y+2);
51  if (m_energy) {
52  m_thermo->setState_HP(y[1]/m_mass, m_pressure, 1.0e-4);
53  } else {
54  m_thermo->setPressure(m_pressure);
55  }
56  m_vol = m_mass / m_thermo->density();
57  updateSurfaceState(y + m_nsp + 2);
58
59  // save parameters needed by other connected reactors
60  m_enthalpy = m_thermo->enthalpy_mass();
61  m_intEnergy = m_thermo->intEnergy_mass();
62  m_thermo->saveState(m_state);
63 }
64
65 void ConstPressureReactor::evalEqs(doublereal time, doublereal* y,
66  doublereal* ydot, doublereal* params)
67 {
68  double dmdt = 0.0; // dm/dt (gas phase)
69  double* dYdt = ydot + 2;
70
71  m_thermo->restoreState(m_state);
72  applySensitivity(params);
73  evalWalls(time);
74  double mdot_surf = evalSurfaces(time, ydot + m_nsp + 2);
75  dmdt += mdot_surf;
76
77  const vector_fp& mw = m_thermo->molecularWeights();
78  const doublereal* Y = m_thermo->massFractions();
79
80  if (m_chem) {
81  m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
82  }
83
84  for (size_t k = 0; k < m_nsp; k++) {
85  // production in gas phase and from surfaces
86  dYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k] / m_mass;
87  // dilution by net surface mass flux
88  dYdt[k] -= Y[k] * mdot_surf / m_mass;
89  }
90
91  // external heat transfer
92  double dHdt = - m_Q;
93
94  // add terms for outlets
95  for (size_t i = 0; i < m_outlet.size(); i++) {
96  double mdot_out = m_outlet[i]->massFlowRate(time); // mass flow out of system
97  dmdt -= mdot_out;
98  dHdt -= mdot_out * m_enthalpy;
99  }
100
101  // add terms for inlets
102  for (size_t i = 0; i < m_inlet.size(); i++) {
103  double mdot_in = m_inlet[i]->massFlowRate(time);
104  dmdt += mdot_in; // mass flow into system
105  for (size_t n = 0; n < m_nsp; n++) {
106  double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
107  // flow of species into system and dilution by other species
108  dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass;
109  }
110  dHdt += mdot_in * m_inlet[i]->enthalpy_mass();
111  }
112
113  ydot[0] = dmdt;
114  if (m_energy) {
115  ydot[1] = dHdt;
116  } else {
117  ydot[1] = 0.0;
118  }
119
120  // reset sensitivity parameters
121  resetSensitivity(params);
122 }
123
124 size_t ConstPressureReactor::componentIndex(const string& nm) const
125 {
126  size_t k = speciesIndex(nm);
127  if (k != npos) {
128  return k + 2;
129  } else if (nm == "m" || nm == "mass") {
130  return 0;
131  } else if (nm == "H" || nm == "enthalpy") {
132  return 1;
133  } else {
134  return npos;
135  }
136 }
137
138 }
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
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