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