Cantera  2.0
Reactor.cpp
Go to the documentation of this file.
1 /**
2  * @file Reactor.cpp
3  *
4  * A zero-dimensional reactor
5  */
6 
7 // Copyright 2001 California Institute of Technology
8 
10 //#include "../CVode.h"
12 #include "cantera/zeroD/Wall.h"
15 
16 using namespace std;
17 
18 namespace Cantera
19 {
20 Reactor::Reactor() : ReactorBase(),
21  m_kin(0),
22  m_temp_atol(1.e-11),
23  m_maxstep(0.0),
24  m_vdot(0.0),
25  m_Q(0.0),
26  m_rtol(1.e-9),
27  m_chem(true),
28  m_energy(true),
29  m_nsens(npos)
30 {}
31 
32 // overloaded method of FuncEval. Called by the integrator to
33 // get the initial conditions.
34 void Reactor::getInitialConditions(double t0, size_t leny, double* y)
35 {
36  m_init = true;
37  if (m_thermo == 0) {
38  cout << "Error: reactor is empty." << endl;
39  return;
40  }
41  m_time = t0;
42  m_thermo->restoreState(m_state);
43 
44  // total mass
45  doublereal mass = m_thermo->density() * m_vol;
46 
47  // set components y + 2 ... y + K + 1 to the
48  // mass M_k of each species
49  m_thermo->getMassFractions(y+2);
50  scale(y + 2, y + m_nsp + 2, y + 2, mass);
51 
52  // set the first component to the total internal
53  // energy
54  y[0] = m_thermo->intEnergy_mass() * mass;
55 
56  // set the second component to the total volume
57  y[1] = m_vol;
58 
59  // set the remaining components to the surface species
60  // coverages on the walls
61  size_t loc = m_nsp + 2;
62  SurfPhase* surf;
63  for (size_t m = 0; m < m_nwalls; m++) {
64  surf = m_wall[m]->surface(m_lr[m]);
65  if (surf) {
66  m_wall[m]->getCoverages(m_lr[m], y + loc);
67  //surf->getCoverages(y+loc);
68  loc += surf->nSpecies();
69  }
70  }
71 }
72 
73 /*
74  * Must be called before calling method 'advance'
75  */
76 void Reactor::initialize(doublereal t0)
77 {
78  m_thermo->restoreState(m_state);
79  m_sdot.resize(m_nsp, 0.0);
80  m_nv = m_nsp + 2;
81  for (size_t w = 0; w < m_nwalls; w++)
82  if (m_wall[w]->surface(m_lr[w])) {
83  m_nv += m_wall[w]->surface(m_lr[w])->nSpecies();
84  }
85 
86  m_enthalpy = m_thermo->enthalpy_mass();
87  m_pressure = m_thermo->pressure();
88  m_intEnergy = m_thermo->intEnergy_mass();
89 
90  size_t nt = 0, maxnt = 0;
91  for (size_t m = 0; m < m_nwalls; m++) {
92  if (m_wall[m]->kinetics(m_lr[m])) {
93  nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies();
94  if (nt > maxnt) {
95  maxnt = nt;
96  }
97  if (m_wall[m]->kinetics(m_lr[m])) {
98  if (&m_kin->thermo(0) !=
99  &m_wall[m]->kinetics(m_lr[m])->thermo(0)) {
100  throw CanteraError("Reactor::initialize",
101  "First phase of all kinetics managers must be"
102  " the gas.");
103  }
104  }
105  }
106  }
107  m_work.resize(maxnt);
108  m_init = true;
109 }
110 
111 size_t Reactor::nSensParams()
112 {
113  if (m_nsens == npos) {
114  // determine the number of sensitivity parameters
115  size_t m, ns;
116  m_nsens = m_pnum.size();
117  for (m = 0; m < m_nwalls; m++) {
118  ns = m_wall[m]->nSensParams(m_lr[m]);
119  m_nsens_wall.push_back(ns);
120  m_nsens += ns;
121  }
122  }
123  return m_nsens;
124 }
125 
126 void Reactor::updateState(doublereal* y)
127 {
128  ThermoPhase& mix = *m_thermo; // define for readability
129 
130  // The components of y are the total internal energy,
131  // the total volume, and the mass of each species.
132 
133  // Set the mass fractions and density of the mixture.
134 
135 
136  doublereal u = y[0];
137  m_vol = y[1];
138  doublereal* mss = y + 2;
139  doublereal mass = accumulate(y+2, y+2+m_nsp, 0.0);
140  m_thermo->setMassFractions(mss);
141 
142  m_thermo->setDensity(mass/m_vol);
143 
144  doublereal temp = temperature();
145  mix.setTemperature(temp);
146 
147  if (m_energy) {
148  // Decreased the tolerance on delta_T to 1.0E-7 so that T is
149  // accurate to 9 sig digits, because this is
150  // used in the numerical jacobian routines where relative values
151  // of 1.0E-7 are used in the deltas.
152  m_thermo->setState_UV(u/mass,m_vol/mass, 1.0e-7);
153  temp = mix.temperature(); //mix.setTemperature(temp);
154  }
155  //m_state[0] = temp;
156 
157  size_t loc = m_nsp + 2;
158  SurfPhase* surf;
159  for (size_t m = 0; m < m_nwalls; m++) {
160  surf = m_wall[m]->surface(m_lr[m]);
161  if (surf) {
162  // surf->setTemperature(temp);
163  //surf->setCoverages(y+loc);
164  m_wall[m]->setCoverages(m_lr[m], y+loc);
165  loc += surf->nSpecies();
166  }
167  }
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 
176 
177 /*
178  * Called by the integrator to evaluate ydot given y at time 'time'.
179  */
180 void Reactor::evalEqs(doublereal time, doublereal* y,
181  doublereal* ydot, doublereal* params)
182 {
183  m_time = time;
184  m_thermo->restoreState(m_state);
185 
186  // process sensitivity parameters
187  if (params) {
188  size_t npar = m_pnum.size();
189  for (size_t n = 0; n < npar; n++) {
190  double mult = m_kin->multiplier(m_pnum[n]);
191  m_kin->setMultiplier(m_pnum[n], mult*params[n]);
192  }
193  size_t ploc = npar;
194  for (size_t m = 0; m < m_nwalls; m++) {
195  if (m_nsens_wall[m] > 0) {
196  m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc);
197  ploc += m_nsens_wall[m];
198  }
199  }
200  }
201 
202  m_vdot = 0.0;
203  m_Q = 0.0;
204 
205  // compute wall terms
206  size_t loc = m_nsp+2;
207  fill(m_sdot.begin(), m_sdot.end(), 0.0);
208  for (size_t i = 0; i < m_nwalls; i++) {
209  int lr = 1 - 2*m_lr[i];
210  double vdot = lr*m_wall[i]->vdot(time);
211  m_vdot += vdot;
212  m_Q += lr*m_wall[i]->Q(time);
213  Kinetics* kin = m_wall[i]->kinetics(m_lr[i]);
214  SurfPhase* surf = m_wall[i]->surface(m_lr[i]);
215  if (surf && kin) {
216  double rs0 = 1.0/surf->siteDensity();
217  size_t nk = surf->nSpecies();
218  double sum = 0.0;
219  surf->setTemperature(m_state[0]);
220  m_wall[i]->syncCoverages(m_lr[i]);
221  kin->getNetProductionRates(DATA_PTR(m_work));
222  size_t ns = kin->surfacePhaseIndex();
223  size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
224  for (size_t k = 1; k < nk; k++) {
225  ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
226  sum -= ydot[loc + k];
227  }
228  ydot[loc] = sum;
229  loc += nk;
230 
231  double wallarea = m_wall[i]->area();
232  for (size_t k = 0; k < m_nsp; k++) {
233  m_sdot[k] += m_work[k]*wallarea;
234  }
235  }
236  }
237 
238  // volume equation
239  ydot[1] = m_vdot;
240 
241  /* species equations
242  * Equation is:
243  * \dot M_k = \hat W_k \dot\omega_k + \dot m_{in} Y_{k,in}
244  * - \dot m_{out} Y_{k} + A \dot s_k.
245  */
246  const vector_fp& mw = m_thermo->molecularWeights();
247  if (m_chem) {
248  m_kin->getNetProductionRates(ydot+2); // "omega dot"
249  } else {
250  fill(ydot + 2, ydot + 2 + m_nsp, 0.0);
251  }
252  for (size_t n = 0; n < m_nsp; n++) {
253  ydot[n+2] *= m_vol; // moles/s/m^3 -> moles/s
254  ydot[n+2] += m_sdot[n];
255  ydot[n+2] *= mw[n];
256  }
257 
258  /*
259  * Energy equation.
260  * \f[
261  * \dot U = -P\dot V + A \dot q + \dot m_{in} h_{in}
262  * - \dot m_{out} h.
263  * \f]
264  */
265  if (m_energy) {
266  ydot[0] = - m_thermo->pressure() * m_vdot - m_Q;
267  } else {
268  ydot[0] = 0.0;
269  }
270 
271  // add terms for open system
272  if (m_open) {
273  const doublereal* mf = m_thermo->massFractions();
274  doublereal enthalpy = m_thermo->enthalpy_mass();
275 
276  // outlets
277  for (size_t i = 0; i < m_nOutlets; i++) {
278  double mdot_out = m_outlet[i]->massFlowRate(time);
279  for (size_t n = 0; n < m_nsp; n++) {
280  ydot[2+n] -= mdot_out * mf[n];
281  }
282  if (m_energy) {
283  ydot[0] -= mdot_out * enthalpy;
284  }
285  }
286 
287  // inlets
288  for (size_t i = 0; i < m_nInlets; i++) {
289  double mdot_in = m_inlet[i]->massFlowRate(time);
290  for (size_t n = 0; n < m_nsp; n++) {
291  ydot[2+n] += m_inlet[i]->outletSpeciesMassFlowRate(n);
292  }
293  if (m_energy) {
294  ydot[0] += mdot_in * m_inlet[i]->enthalpy_mass();
295  }
296  }
297  }
298 
299  // reset sensitivity parameters
300  if (params) {
301  size_t npar = m_pnum.size();
302  for (size_t n = 0; n < npar; n++) {
303  double mult = m_kin->multiplier(m_pnum[n]);
304  m_kin->setMultiplier(m_pnum[n], mult/params[n]);
305  }
306  size_t ploc = npar;
307  for (size_t m = 0; m < m_nwalls; m++) {
308  if (m_nsens_wall[m] > 0) {
309  m_wall[m]->resetSensitivityParameters(m_lr[m]);
310  ploc += m_nsens_wall[m];
311  }
312  }
313  }
314 }
315 
316 void Reactor::addSensitivityReaction(size_t rxn)
317 {
318  if (rxn >= m_kin->nReactions())
319  throw CanteraError("Reactor::addSensitivityReaction",
320  "Reaction number out of range ("+int2str(rxn)+")");
321  m_pnum.push_back(rxn);
322  m_pname.push_back(name()+": "+m_kin->reactionString(rxn));
323  m_mult_save.push_back(1.0);
324 }
325 
326 
327 size_t Reactor::componentIndex(string nm) const
328 {
329  if (nm == "U") {
330  return 0;
331  }
332  if (nm == "V") {
333  return 1;
334  }
335  // check for a gas species name
336  size_t k = m_thermo->speciesIndex(nm);
337  if (k != npos) {
338  return k + 2;
339  }
340 
341  // check for a wall species
342  size_t walloffset = 0, kp = 0;
343  thermo_t* th;
344  for (size_t m = 0; m < m_nwalls; m++) {
345  if (m_wall[m]->kinetics(m_lr[m])) {
346  kp = m_wall[m]->kinetics(m_lr[m])->reactionPhaseIndex();
347  th = &m_wall[m]->kinetics(m_lr[m])->thermo(kp);
348  k = th->speciesIndex(nm);
349  if (k != npos) {
350  return k + 2 + m_nsp + walloffset;
351  } else {
352  walloffset += th->nSpecies();
353  }
354  }
355  }
356  return npos;
357 }
358 
359 }