Cantera  2.0
ConstPressureReactor.cpp
1 /**
2  * @file Reactor.cpp
3  *
4  * A zero-dimensional reactor
5  */
6 
7 // Copyright 2001 California Institute of Technology
8 
9 #include "cantera/zeroD/ConstPressureReactor.h"
11 #include "cantera/zeroD/Wall.h"
14 
15 using namespace std;
16 
17 namespace Cantera
18 {
19 
20 ConstPressureReactor::ConstPressureReactor() : Reactor() {}
21 
22 void ConstPressureReactor::
23 getInitialConditions(double t0, size_t leny, double* y)
24 {
25  m_init = true;
26  if (m_thermo == 0) {
27  throw CanteraError("getInitialConditions",
28  "Error: reactor is empty.");
29  }
30  m_time = t0;
31  m_thermo->restoreState(m_state);
32 
33  // total mass
34  doublereal mass = m_thermo->density() * m_vol;
35 
36  // set components y + 2 ... y + K + 1 to the
37  // mass M_k of each species
38  m_thermo->getMassFractions(y+2);
39  scale(y + 2, y + m_nsp + 2, y + 2, mass);
40 
41  // set the first component to the total enthalpy
42  y[0] = m_thermo->enthalpy_mass() * mass;
43 
44  // set the second component to the total volume
45  y[1] = m_vol;
46 
47  // set the remaining components to the surface species
48  // coverages on the walls
49  size_t loc = m_nsp + 2;
50  SurfPhase* surf;
51  for (size_t m = 0; m < m_nwalls; m++) {
52  surf = m_wall[m]->surface(m_lr[m]);
53  if (surf) {
54  m_wall[m]->getCoverages(m_lr[m], y + loc);
55  loc += surf->nSpecies();
56  }
57  }
58 }
59 
61 {
62  m_thermo->restoreState(m_state);
63  m_sdot.resize(m_nsp, 0.0);
64  m_nv = m_nsp + 2;
65  for (size_t w = 0; w < m_nwalls; w++)
66  if (m_wall[w]->surface(m_lr[w])) {
67  m_nv += m_wall[w]->surface(m_lr[w])->nSpecies();
68  }
69  m_enthalpy = m_thermo->enthalpy_mass();
70  m_pressure = m_thermo->pressure();
71  m_intEnergy = m_thermo->intEnergy_mass();
72 
73  size_t nt = 0, maxnt = 0;
74  for (size_t m = 0; m < m_nwalls; m++) {
75  if (m_wall[m]->kinetics(m_lr[m])) {
76  nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies();
77  if (nt > maxnt) {
78  maxnt = nt;
79  }
80  if (m_wall[m]->kinetics(m_lr[m])) {
81  if (&m_kin->thermo(0) !=
82  &m_wall[m]->kinetics(m_lr[m])->thermo(0)) {
83  throw CanteraError("ConstPressureReactor::initialize",
84  "First phase of all kinetics managers must be"
85  " the gas.");
86  }
87  }
88  }
89  }
90  m_work.resize(maxnt);
91  m_init = true;
92 }
93 
95 {
96 
97  // The components of y are the total enthalpy,
98  // the total volume, and the mass of each species.
99  doublereal h = y[0];
100  doublereal* mss = y + 2;
101  doublereal mass = accumulate(y+2, y+2+m_nsp, 0.0);
102  m_thermo->setMassFractions(mss);
103 
104  if (m_energy) {
105  m_thermo->setState_HP(h/mass, m_pressure, 1.0e-4);
106  } else {
107  m_thermo->setPressure(m_pressure);
108  }
109  m_vol = mass / m_thermo->density();
110 
111  size_t loc = m_nsp + 2;
112  SurfPhase* surf;
113  for (size_t m = 0; m < m_nwalls; m++) {
114  surf = m_wall[m]->surface(m_lr[m]);
115  if (surf) {
116  m_wall[m]->setCoverages(m_lr[m], y+loc);
117  loc += surf->nSpecies();
118  }
119  }
120 
121  // save parameters needed by other connected reactors
122  m_enthalpy = m_thermo->enthalpy_mass();
123  m_intEnergy = m_thermo->intEnergy_mass();
124  m_thermo->saveState(m_state);
125 }
126 
127 
128 /*
129  * Called by the integrator to evaluate ydot given y at time 'time'.
130  */
131 void ConstPressureReactor::evalEqs(doublereal time, doublereal* y,
132  doublereal* ydot, doublereal* params)
133 {
134  size_t nk;
135  m_time = time;
136  m_thermo->restoreState(m_state);
137 
138  Kinetics* kin;
139  size_t npar, ploc;
140  double mult;
141 
142  // process sensitivity parameters
143  if (params) {
144 
145  npar = m_pnum.size();
146  for (size_t n = 0; n < npar; n++) {
147  mult = m_kin->multiplier(m_pnum[n]);
148  m_kin->setMultiplier(m_pnum[n], mult*params[n]);
149  }
150  ploc = npar;
151  for (size_t m = 0; m < m_nwalls; m++) {
152  if (m_nsens_wall[m] > 0) {
153  m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc);
154  ploc += m_nsens_wall[m];
155  }
156  }
157  }
158 
159  m_vdot = 0.0;
160  m_Q = 0.0;
161 
162  // compute wall terms
163  doublereal rs0, sum, wallarea;
164 
165  SurfPhase* surf;
166  size_t lr, ns, loc = m_nsp+2, surfloc;
167  fill(m_sdot.begin(), m_sdot.end(), 0.0);
168  for (size_t i = 0; i < m_nwalls; i++) {
169  lr = 1 - 2*m_lr[i];
170  m_Q += lr*m_wall[i]->Q(time);
171  kin = m_wall[i]->kinetics(m_lr[i]);
172  surf = m_wall[i]->surface(m_lr[i]);
173  if (surf && kin) {
174  rs0 = 1.0/surf->siteDensity();
175  nk = surf->nSpecies();
176  sum = 0.0;
177  surf->setTemperature(m_state[0]);
178  m_wall[i]->syncCoverages(m_lr[i]);
179  kin->getNetProductionRates(DATA_PTR(m_work));
180  ns = kin->surfacePhaseIndex();
181  surfloc = kin->kineticsSpeciesIndex(0,ns);
182  for (size_t k = 1; k < nk; k++) {
183  ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
184  sum -= ydot[loc + k];
185  }
186  ydot[loc] = sum;
187  loc += nk;
188 
189  wallarea = m_wall[i]->area();
190  for (size_t k = 0; k < m_nsp; k++) {
191  m_sdot[k] += m_work[k]*wallarea;
192  }
193  }
194  }
195 
196  // dummy equation
197  ydot[1] = 0.0;
198 
199  /* species equations
200  * Equation is:
201  * \dot M_k = \hat W_k \dot\omega_k + \dot m_{in} Y_{k,in}
202  * - \dot m_{out} Y_{k} + A \dot s_k.
203  */
204  const doublereal* mw = DATA_PTR(m_thermo->molecularWeights());
205  if (m_chem) {
206  m_kin->getNetProductionRates(ydot+2); // "omega dot"
207  } else {
208  fill(ydot + 2, ydot + 2 + m_nsp, 0.0);
209  }
210  for (size_t n = 0; n < m_nsp; n++) {
211  ydot[n+2] *= m_vol; // moles/s/m^3 -> moles/s
212  ydot[n+2] += m_sdot[n];
213  ydot[n+2] *= mw[n];
214  }
215 
216 
217  /*
218  * Energy equation.
219  * \f[
220  * \dot U = -P\dot V + A \dot q + \dot m_{in} h_{in}
221  * - \dot m_{out} h.
222  * \f]
223  */
224  if (m_energy) {
225  ydot[0] = - m_Q;
226  } else {
227  ydot[0] = 0.0;
228  }
229 
230  // add terms for open system
231  if (m_open) {
232 
233  const doublereal* mf = m_thermo->massFractions();
234  doublereal enthalpy = m_thermo->enthalpy_mass();
235 
236  // outlets
237 
238  doublereal mdot_out;
239  for (size_t i = 0; i < m_nOutlets; i++) {
240  mdot_out = m_outlet[i]->massFlowRate(time);
241  for (size_t n = 0; n < m_nsp; n++) {
242  ydot[2+n] -= mdot_out * mf[n];
243  }
244  if (m_energy) {
245  ydot[0] -= mdot_out * enthalpy;
246  }
247  }
248 
249 
250  // inlets
251 
252  doublereal mdot_in;
253  for (size_t i = 0; i < m_nInlets; i++) {
254  mdot_in = m_inlet[i]->massFlowRate(time);
255  for (size_t n = 0; n < m_nsp; n++) {
256  ydot[2+n] += m_inlet[i]->outletSpeciesMassFlowRate(n);
257  }
258  if (m_energy) {
259  ydot[0] += mdot_in * m_inlet[i]->enthalpy_mass();
260  }
261  }
262  }
263 
264  // reset sensitivity parameters
265  if (params) {
266  npar = m_pnum.size();
267  for (size_t n = 0; n < npar; n++) {
268  mult = m_kin->multiplier(m_pnum[n]);
269  m_kin->setMultiplier(m_pnum[n], mult/params[n]);
270  }
271  ploc = npar;
272  for (size_t m = 0; m < m_nwalls; m++) {
273  if (m_nsens_wall[m] > 0) {
274  m_wall[m]->resetSensitivityParameters(m_lr[m]);
275  ploc += m_nsens_wall[m];
276  }
277  }
278  }
279 }
280 
281 size_t ConstPressureReactor::componentIndex(string nm) const
282 {
283  if (nm == "H") {
284  return 0;
285  }
286  if (nm == "V") {
287  return 1;
288  }
289  // check for a gas species name
290  size_t k = m_thermo->speciesIndex(nm);
291  if (k != npos) {
292  return k + 2;
293  }
294 
295  // check for a wall species
296  size_t walloffset = 0, kp = 0;
297  thermo_t* th;
298  for (size_t m = 0; m < m_nwalls; m++) {
299  if (m_wall[m]->kinetics(m_lr[m])) {
300  kp = m_wall[m]->kinetics(m_lr[m])->reactionPhaseIndex();
301  th = &m_wall[m]->kinetics(m_lr[m])->thermo(kp);
302  k = th->speciesIndex(nm);
303  if (k != npos) {
304  return k + 2 + m_nsp + walloffset;
305  } else {
306  walloffset += th->nSpecies();
307  }
308  }
309  }
310  return npos;
311 }
312 
313 }