Cantera  2.3.0
Reactor.cpp
Go to the documentation of this file.
1 //! @file Reactor.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"
12 
13 #include <boost/math/tools/roots.hpp>
14 
15 using namespace std;
16 namespace bmt = boost::math::tools;
17 
18 namespace Cantera
19 {
20 Reactor::Reactor() :
21  m_kin(0),
22  m_vdot(0.0),
23  m_Q(0.0),
24  m_mass(0.0),
25  m_chem(false),
26  m_energy(true),
27  m_nv(0)
28 {}
29 
30 void Reactor::setKineticsMgr(Kinetics& kin)
31 {
32  m_kin = &kin;
33  if (m_kin->nReactions() == 0) {
34  setChemistry(false);
35  } else {
36  setChemistry(true);
37  }
38 }
39 
40 void Reactor::getInitialConditions(double t0, size_t leny, double* y)
41 {
42  warn_deprecated("Reactor::getInitialConditions",
43  "Use getState instead. To be removed after Cantera 2.3.");
44  getState(y);
45 }
46 
47 void Reactor::getState(double* y)
48 {
49  if (m_thermo == 0) {
50  throw CanteraError("getState",
51  "Error: reactor is empty.");
52  }
53  m_thermo->restoreState(m_state);
54 
55  // set the first component to the total mass
56  m_mass = m_thermo->density() * m_vol;
57  y[0] = m_mass;
58 
59  // set the second component to the total volume
60  y[1] = m_vol;
61 
62  // set the third component to the total internal energy
63  y[2] = m_thermo->intEnergy_mass() * m_mass;
64 
65  // set components y+3 ... y+K+2 to the mass fractions of each species
66  m_thermo->getMassFractions(y+3);
67 
68  // set the remaining components to the surface species
69  // coverages on the walls
71 }
72 
74 {
75  size_t loc = 0;
76  for (auto& S : m_surfaces) {
77  S->getCoverages(y + loc);
78  loc += S->thermo()->nSpecies();
79  }
80 }
81 
82 void Reactor::initialize(doublereal t0)
83 {
84  if (!m_thermo || (m_chem && !m_kin)) {
85  throw CanteraError("Reactor::initialize", "Reactor contents not set"
86  " for reactor '" + m_name + "'.");
87  }
88  m_thermo->restoreState(m_state);
89  m_sdot.resize(m_nsp, 0.0);
90  m_wdot.resize(m_nsp, 0.0);
91 
92  m_enthalpy = m_thermo->enthalpy_mass();
93  m_pressure = m_thermo->pressure();
94  m_intEnergy = m_thermo->intEnergy_mass();
95 
96  for (size_t n = 0; n < m_wall.size(); n++) {
97  Wall* W = m_wall[n];
98  W->initialize();
99  if (W->kinetics(m_lr[n])) {
100  addSurface(W->reactorSurface(m_lr[n]));
101  }
102  }
103 
104  m_nv = m_nsp + 3;
105  size_t maxnt = 0;
106  for (auto& S : m_surfaces) {
107  m_nv += S->thermo()->nSpecies();
108  size_t nt = S->kinetics()->nTotalSpecies();
109  maxnt = std::max(maxnt, nt);
110  if (m_chem && &m_kin->thermo(0) != &S->kinetics()->thermo(0)) {
111  throw CanteraError("Reactor::initialize",
112  "First phase of all kinetics managers must be the gas.");
113  }
114  }
115  m_work.resize(maxnt);
116 }
117 
119 {
120  size_t ns = m_sensParams.size();
121  for (auto& S : m_surfaces) {
122  ns += S->nSensParams();
123  }
124  return ns;
125 }
126 
128 {
130  m_mass = m_thermo->density() * m_vol;
131 }
132 
133 void Reactor::updateState(doublereal* y)
134 {
135  // The components of y are [0] the total mass, [1] the total volume,
136  // [2] the total internal energy, [3...K+3] are the mass fractions of each
137  // species, and [K+3...] are the coverages of surface species on each wall.
138  m_mass = y[0];
139  m_vol = y[1];
140  m_thermo->setMassFractions_NoNorm(y+3);
141 
142  if (m_energy) {
143  double U = y[2];
144  // Residual function: error in internal energy as a function of T
145  auto u_err = [this, U](double T) {
146  m_thermo->setState_TR(T, m_mass / m_vol);
147  return m_thermo->intEnergy_mass() * m_mass - U;
148  };
149 
150  double T = m_thermo->temperature();
151  boost::uintmax_t maxiter = 100;
152  std::pair<double, double> TT;
153  try {
154  TT = bmt::bracket_and_solve_root(
155  u_err, T, 1.2, true, bmt::eps_tolerance<double>(48), maxiter);
156  } catch (std::exception& err) {
157  // Set m_thermo back to a reasonable state if root finding fails
158  m_thermo->setState_TR(T, m_mass / m_vol);
159  throw CanteraError("Reactor::updateState", err.what());
160  }
161  if (fabs(TT.first - TT.second) > 1e-7*TT.first) {
162  throw CanteraError("Reactor::updateState",
163  "bracket_and_solve_root failed");
164  }
165  m_thermo->setState_TR(TT.second, m_mass / m_vol);
166  } else {
167  m_thermo->setDensity(m_mass/m_vol);
168  }
169 
170  updateSurfaceState(y + m_nsp + 3);
171 
172  // save parameters needed by other connected reactors
173  m_enthalpy = m_thermo->enthalpy_mass();
174  m_pressure = m_thermo->pressure();
175  m_intEnergy = m_thermo->intEnergy_mass();
176  m_thermo->saveState(m_state);
177 }
178 
180 {
181  size_t loc = 0;
182  for (auto& S : m_surfaces) {
183  S->setCoverages(y+loc);
184  loc += S->thermo()->nSpecies();
185  }
186 }
187 
188 void Reactor::evalEqs(doublereal time, doublereal* y,
189  doublereal* ydot, doublereal* params)
190 {
191  double dmdt = 0.0; // dm/dt (gas phase)
192  double* dYdt = ydot + 3;
193 
194  m_thermo->restoreState(m_state);
195  applySensitivity(params);
196  evalWalls(time);
197  double mdot_surf = evalSurfaces(time, ydot + m_nsp + 3);
198  dmdt += mdot_surf; // mass added to gas phase from surface reactions
199 
200  // volume equation
201  ydot[1] = m_vdot;
202 
203  const vector_fp& mw = m_thermo->molecularWeights();
204  const doublereal* Y = m_thermo->massFractions();
205 
206  if (m_chem) {
207  m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
208  }
209 
210  for (size_t k = 0; k < m_nsp; k++) {
211  // production in gas phase and from surfaces
212  dYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k] / m_mass;
213  // dilution by net surface mass flux
214  dYdt[k] -= Y[k] * mdot_surf / m_mass;
215  }
216 
217  // Energy equation.
218  // \f[
219  // \dot U = -P\dot V + A \dot q + \dot m_{in} h_{in} - \dot m_{out} h.
220  // \f]
221  if (m_energy) {
222  ydot[2] = - m_thermo->pressure() * m_vdot - m_Q;
223  } else {
224  ydot[2] = 0.0;
225  }
226 
227  // add terms for outlets
228  for (size_t i = 0; i < m_outlet.size(); i++) {
229  double mdot_out = m_outlet[i]->massFlowRate(time);
230  dmdt -= mdot_out; // mass flow out of system
231  if (m_energy) {
232  ydot[2] -= mdot_out * m_enthalpy;
233  }
234  }
235 
236  // add terms for inlets
237  for (size_t i = 0; i < m_inlet.size(); i++) {
238  double mdot_in = m_inlet[i]->massFlowRate(time);
239  dmdt += mdot_in; // mass flow into system
240  for (size_t n = 0; n < m_nsp; n++) {
241  double mdot_spec = m_inlet[i]->outletSpeciesMassFlowRate(n);
242  // flow of species into system and dilution by other species
243  dYdt[n] += (mdot_spec - mdot_in * Y[n]) / m_mass;
244  }
245  if (m_energy) {
246  ydot[2] += mdot_in * m_inlet[i]->enthalpy_mass();
247  }
248  }
249 
250  ydot[0] = dmdt;
251  resetSensitivity(params);
252 }
253 
254 void Reactor::evalWalls(double t)
255 {
256  m_vdot = 0.0;
257  m_Q = 0.0;
258  for (size_t i = 0; i < m_wall.size(); i++) {
259  int lr = 1 - 2*m_lr[i];
260  m_vdot += lr*m_wall[i]->vdot(t);
261  m_Q += lr*m_wall[i]->Q(t);
262  }
263 }
264 
265 double Reactor::evalSurfaces(double t, double* ydot)
266 {
267  const vector_fp& mw = m_thermo->molecularWeights();
268  fill(m_sdot.begin(), m_sdot.end(), 0.0);
269  size_t loc = 0; // offset into ydot
270  double mdot_surf = 0.0; // net mass flux from surface
271 
272  for (auto S : m_surfaces) {
273  Kinetics* kin = S->kinetics();
274  SurfPhase* surf = S->thermo();
275 
276  double rs0 = 1.0/surf->siteDensity();
277  size_t nk = surf->nSpecies();
278  double sum = 0.0;
279  surf->setTemperature(m_state[0]);
280  S->syncCoverages();
281  kin->getNetProductionRates(&m_work[0]);
282  size_t ns = kin->surfacePhaseIndex();
283  size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
284  for (size_t k = 1; k < nk; k++) {
285  ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
286  sum -= ydot[loc + k];
287  }
288  ydot[loc] = sum;
289  loc += nk;
290 
291  double wallarea = S->area();
292  for (size_t k = 0; k < m_nsp; k++) {
293  m_sdot[k] += m_work[k]*wallarea;
294  mdot_surf += m_sdot[k] * mw[k];
295  }
296  }
297  return mdot_surf;
298 }
299 
301 {
302  if (!m_chem || rxn >= m_kin->nReactions()) {
303  throw CanteraError("Reactor::addSensitivityReaction",
304  "Reaction number out of range ({})", rxn);
305  }
306 
308  name()+": "+m_kin->reactionString(rxn), 1.0, 1.0);
309  m_sensParams.emplace_back(
310  SensitivityParameter{rxn, p, 1.0, SensParameterType::reaction});
311 }
312 
314 {
315  if (k >= m_thermo->nSpecies()) {
316  throw CanteraError("Reactor::addSensitivitySpeciesEnthalpy",
317  "Species index out of range ({})", k);
318  }
319 
321  name() + ": " + m_thermo->speciesName(k) + " enthalpy",
322  0.0, GasConstant * 298.15);
323  m_sensParams.emplace_back(
324  SensitivityParameter{k, p, m_thermo->Hf298SS(k),
325  SensParameterType::enthalpy});
326 }
327 
328 size_t Reactor::speciesIndex(const string& nm) const
329 {
330  // check for a gas species name
331  size_t k = m_thermo->speciesIndex(nm);
332  if (k != npos) {
333  return k;
334  }
335 
336  // check for a wall species
337  size_t offset = m_nsp;
338  for (auto& S : m_surfaces) {
339  ThermoPhase* th = S->thermo();
340  k = th->speciesIndex(nm);
341  if (k != npos) {
342  return k + offset;
343  } else {
344  offset += th->nSpecies();
345  }
346  }
347  return npos;
348 }
349 
350 size_t Reactor::componentIndex(const string& nm) const
351 {
352  size_t k = speciesIndex(nm);
353  if (k != npos) {
354  return k + 3;
355  } else if (nm == "m" || nm == "mass") {
356  if (nm == "m") {
357  warn_deprecated("Reactor::componentIndex(\"m\")",
358  "Using the name 'm' for mass is deprecated, and will be "
359  "disabled after Cantera 2.3. Use 'mass' instead.");
360  }
361  return 0;
362  } else if (nm == "V" || nm == "volume") {
363  if (nm == "V") {
364  warn_deprecated("Reactor::componentIndex(\"V\")",
365  "Using the name 'V' for volume is deprecated, and will be "
366  "disabled after Cantera 2.3. Use 'volume' instead.");
367  }
368  return 1;
369  } else if (nm == "U" || nm == "int_energy") {
370  if (nm == "U") {
371  warn_deprecated("Reactor::componentIndex(\"U\")",
372  "Using the name 'U' for internal energy is deprecated, and "
373  "will be disabled after Cantera 2.3. Use 'int_energy' instead.");
374  }
375  return 2;
376  } else {
377  return npos;
378  }
379 }
380 
381 std::string Reactor::componentName(size_t k) {
382  if (k == 0) {
383  return "mass";
384  } else if (k == 1) {
385  return "volume";
386  } else if (k == 2) {
387  return "int_energy";
388  } else if (k >= 3 && k < neq()) {
389  k -= 3;
390  if (k < m_thermo->nSpecies()) {
391  return m_thermo->speciesName(k);
392  } else {
393  k -= m_thermo->nSpecies();
394  }
395  for (auto& S : m_surfaces) {
396  ThermoPhase* th = S->thermo();
397  if (k < th->nSpecies()) {
398  return th->speciesName(k);
399  } else {
400  k -= th->nSpecies();
401  }
402  }
403  }
404  throw CanteraError("Reactor::componentName", "Index is out of bounds.");
405 }
406 
407 void Reactor::applySensitivity(double* params)
408 {
409  if (!params) {
410  return;
411  }
412  for (auto& p : m_sensParams) {
413  if (p.type == SensParameterType::reaction) {
414  p.value = m_kin->multiplier(p.local);
415  m_kin->setMultiplier(p.local, p.value*params[p.global]);
416  } else if (p.type == SensParameterType::enthalpy) {
417  m_thermo->modifyOneHf298SS(p.local, p.value + params[p.global]);
418  }
419  }
420  for (auto& S : m_surfaces) {
421  S->setSensitivityParameters(params);
422  }
423  m_thermo->invalidateCache();
424  if (m_kin) {
425  m_kin->invalidateCache();
426  }
427 }
428 
429 void Reactor::resetSensitivity(double* params)
430 {
431  if (!params) {
432  return;
433  }
434  for (auto& p : m_sensParams) {
435  if (p.type == SensParameterType::reaction) {
436  m_kin->setMultiplier(p.local, p.value);
437  } else if (p.type == SensParameterType::enthalpy) {
438  m_thermo->resetHf298(p.local);
439  }
440  }
441  for (auto& S : m_surfaces) {
442  S->resetSensitivityParameters();
443  }
444  m_thermo->invalidateCache();
445  if (m_kin) {
446  m_kin->invalidateCache();
447  }
448 }
449 
450 }
const vector_fp & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition: Phase.cpp:513
void setChemistry(bool cflag=true)
Enable or disable changes in reactor composition due to chemical reactions.
Definition: Reactor.h:75
virtual void getInitialConditions(doublereal t0, size_t leny, doublereal *y)
Called by ReactorNet to get the initial conditions.
Definition: Reactor.cpp:40
vector_fp m_sdot
Production rates of gas phase species on surfaces [kmol/s].
Definition: Reactor.h:203
void restoreState(const vector_fp &state)
Restore a state saved on a previous call to saveState.
Definition: Phase.cpp:309
doublereal m_mass
total mass
Definition: Reactor.h:199
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase, assuming an ideal solution model (see Thermodynamic Properties and class SurfPhase).
virtual size_t nSensParams()
Number of sensitivity parameters associated with this reactor (including walls)
Definition: Reactor.cpp:118
void getMassFractions(doublereal *const y) const
Get the species mass fractions.
Definition: Phase.cpp:584
Header file for class Wall.
vector_fp m_wdot
Species net molar production rates.
Definition: Reactor.h:205
doublereal temperature() const
Temperature (K).
Definition: Phase.h:601
const doublereal * massFractions() const
Return a const pointer to the mass fraction array.
Definition: Phase.h:478
size_t speciesIndex(const std::string &name) const
Returns the index of a species named &#39;name&#39; within the Phase object.
Definition: Phase.cpp:251
thermo_t & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism...
Definition: Kinetics.h:276
void saveState(vector_fp &state) const
Save the current internal state of the phase.
Definition: Phase.cpp:297
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
virtual void addSensitivityReaction(size_t rxn)
Add a sensitivity parameter associated with the reaction number rxn (in the homogeneous phase)...
Definition: Reactor.cpp:300
ReactorNet & network()
The ReactorNet that this reactor belongs to.
Definition: ReactorBase.cpp:85
virtual void evalWalls(double t)
Evaluate terms related to Walls.
Definition: Reactor.cpp:254
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:54
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:262
virtual void resetSensitivity(double *params)
Reset the reaction rate multipliers.
Definition: Reactor.cpp:429
STL namespace.
virtual doublereal density() const
Density (kg/m^3).
Definition: Phase.h:607
doublereal enthalpy_mass() const
Specific enthalpy. Units: J/kg.
Definition: ThermoPhase.h:764
virtual size_t speciesIndex(const std::string &nm) const
Return the index in the solution vector for this reactor of the species named nm, in either the homog...
Definition: Reactor.cpp:328
virtual std::string componentName(size_t k)
Return the name of the solution component with index i.
Definition: Reactor.cpp:381
virtual void initialize(doublereal t0=0.0)
Initialize the reactor.
Definition: Reactor.cpp:82
size_t registerSensitivityParameter(const std::string &name, double value, double scale)
Used by Reactor and Wall objects to register the addition of sensitivity parameters so that the React...
Definition: ReactorNet.cpp:252
virtual void getNetProductionRates(doublereal *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition: Kinetics.cpp:473
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:93
virtual void updateState(doublereal *y)
Set the state of the reactor to correspond to the state vector y.
Definition: Reactor.cpp:133
doublereal Hf298SS(const size_t k) const
Report the 298 K Heat of Formation of the standard state of one species (J kmol-1) ...
Definition: ThermoPhase.h:179
virtual void invalidateCache()
Invalidate any cached values which are normally updated only when a change in state is detected...
Kinetics * m_kin
Pointer to the homogeneous Kinetics object that handles the reactions.
Definition: Reactor.h:195
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition: SurfPhase.h:143
size_t surfacePhaseIndex()
This returns the integer index of the phase which has ThermoPhase type cSurf.
Definition: Kinetics.h:251
ReactorSurface * reactorSurface(int leftright)
Definition: Wall.h:162
virtual void applySensitivity(double *params)
Set reaction rate multipliers based on the sensitivity variables in params.
Definition: Reactor.cpp:407
doublereal multiplier(size_t i) const
The current value of the multiplier for reaction i.
Definition: Kinetics.h:821
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:267
Represents a wall between between two ReactorBase objects.
Definition: Wall.h:25
doublereal m_Q
net heat transfer through walls [W]
Definition: Reactor.h:198
virtual void syncState()
Set the state of the reactor to correspond to the state of the associated ThermoPhase object...
Definition: ReactorBase.cpp:36
virtual double evalSurfaces(double t, double *ydot)
Evaluate terms related to surface reactions.
Definition: Reactor.cpp:265
Public interface for kinetics managers.
Definition: Kinetics.h:111
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition: Kinetics.h:310
virtual void getState(doublereal *y)
Get the the current state of the reactor.
Definition: Reactor.cpp:47
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
size_t nReactions() const
Number of reactions in the reaction mechanism.
Definition: Kinetics.h:184
std::string reactionString(size_t i) const
Return a string representing the reaction.
Definition: Kinetics.h:671
virtual void initialize()
Called just before the start of integration.
Definition: Wall.h:127
virtual void resetHf298(const size_t k=npos)
Restore the original heat of formation of one or more species.
Header file for class ReactorSurface.
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
Definition: ThermoPhase.h:278
Kinetics * kinetics(int leftright)
Return a pointer to the surface kinetics object for the left (leftright=0) or right (leftright=1) wal...
Definition: Wall.h:170
virtual void setMassFractions_NoNorm(const doublereal *const y)
Set the mass fractions to the specified values without normalizing.
Definition: Phase.cpp:401
virtual void syncState()
Set the state of the reactor to correspond to the state of the associated ThermoPhase object...
Definition: Reactor.cpp:127
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 size_t componentIndex(const std::string &nm) const
Return the index in the solution vector for this reactor of the component named nm.
Definition: Reactor.cpp:350
virtual void setTemperature(const doublereal temp)
Set the internally stored temperature of the phase (K).
Definition: Phase.h:637
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
virtual void modifyOneHf298SS(const size_t k, const doublereal Hf298New)
Modify the value of the 298 K Heat of Formation of one species in the phase (J kmol-1) ...
Definition: ThermoPhase.h:194
size_t m_nsp
Number of homogeneous species in the mixture.
Definition: ReactorBase.h:231
virtual size_t neq()
Number of equations (state variables) for this reactor.
Definition: Reactor.h:99
doublereal size(size_t k) const
This routine returns the size of species k.
Definition: Phase.h:414
virtual void setMultiplier(size_t i, doublereal f)
Set the multiplier for reaction i to f.
Definition: Kinetics.h:830
virtual void evalEqs(doublereal t, doublereal *y, doublereal *ydot, doublereal *params)
Definition: Reactor.cpp:188
doublereal m_vdot
net rate of volume change from moving walls [m^3/s]
Definition: Reactor.h:197
Namespace for the Cantera kernel.
Definition: application.cpp:29
virtual void addSensitivitySpeciesEnthalpy(size_t k)
Add a sensitivity parameter associated with the enthalpy formation of species k (in the homogeneous p...
Definition: Reactor.cpp:313
doublereal intEnergy_mass() const
Specific internal energy. Units: J/kg.
Definition: ThermoPhase.h:769
void setState_TR(doublereal t, doublereal rho)
Set the internally stored temperature (K) and density (kg/m^3)
Definition: Phase.cpp:466
virtual void updateSurfaceState(double *y)
Update the state of SurfPhase objects attached to this reactor.
Definition: Reactor.cpp:179
virtual void setDensity(const doublereal density_)
Set the internally stored density (kg/m^3) of the phase.
Definition: Phase.h:622
virtual void getSurfaceInitialConditions(double *y)
Get initial conditions for SurfPhase objects attached to this reactor.
Definition: Reactor.cpp:73
doublereal siteDensity()
Returns the site density.
Definition: SurfPhase.h:327
std::string name() const
Return the name of this reactor.
Definition: ReactorBase.h:56