Cantera  3.1.0a1
Loading...
Searching...
No Matches
MoleReactor.cpp
Go to the documentation of this file.
1//! @file MoleReactor.cpp A zero-dimensional reactor with a moles as the
2//! state
3
4// This file is part of Cantera. See License.txt in the top-level directory or
5// at https://cantera.org/license.txt for license and copyright information.
6
13#include <boost/math/tools/roots.hpp>
14
15using namespace std;
16namespace bmt = boost::math::tools;
17
18namespace Cantera
19{
20
22{
23 size_t loc = 0;
24 for (auto& S : m_surfaces) {
25 double area = S->area();
26 auto currPhase = S->thermo();
27 size_t tempLoc = currPhase->nSpecies();
28 double surfDensity = currPhase->siteDensity();
29 S->getCoverages(y + loc);
30 // convert coverages to moles
31 for (size_t i = 0; i < tempLoc; i++) {
32 y[i + loc] = y[i + loc] * area * surfDensity / currPhase->size(i);
33 }
34 loc += tempLoc;
35 }
36}
37
39{
41 m_nv -= 1; // moles gives the state one fewer variables
42}
43
45{
46 size_t loc = 0;
47 vector<double> coverages(m_nv_surf, 0.0);
48 for (auto& S : m_surfaces) {
49 auto surf = S->thermo();
50 double invArea = 1/S->area();
51 double invSurfDensity = 1/surf->siteDensity();
52 size_t tempLoc = surf->nSpecies();
53 for (size_t i = 0; i < tempLoc; i++) {
54 coverages[i + loc] = y[i + loc] * invArea * surf->size(i) * invSurfDensity;
55 }
56 S->setCoverages(coverages.data()+loc);
57 loc += tempLoc;
58 }
59}
60
61void MoleReactor::evalSurfaces(double* LHS, double* RHS, double* sdot)
62{
63 fill(sdot, sdot + m_nsp, 0.0);
64 size_t loc = 0; // offset into ydot
65 for (auto S : m_surfaces) {
66 Kinetics* kin = S->kinetics();
67 SurfPhase* surf = S->thermo();
68 double wallarea = S->area();
69 size_t nk = surf->nSpecies();
70 S->syncState();
71 kin->getNetProductionRates(&m_work[0]);
72 for (size_t k = 0; k < nk; k++) {
73 RHS[loc + k] = m_work[k] * wallarea / surf->size(k);
74 }
75 loc += nk;
76
77 size_t bulkloc = kin->kineticsSpeciesIndex(m_thermo->speciesName(0));
78
79 for (size_t k = 0; k < m_nsp; k++) {
80 sdot[k] += m_work[bulkloc + k] * wallarea;
81 }
82 }
83}
84
85void MoleReactor::addSurfaceJacobian(vector<Eigen::Triplet<double>> &triplets)
86{
87 size_t offset = m_nsp;
88 for (auto& S : m_surfaces) {
89 S->syncState();
90 double A = S->area();
91 auto kin = S->kinetics();
92 size_t nk = S->thermo()->nSpecies();
93 // index of gas and surface phases to check if the species is in gas or surface
94 size_t spi = 0;
95 size_t gpi = kin->speciesPhaseIndex(kin->kineticsSpeciesIndex(
96 m_thermo->speciesName(0)));
97 // get surface jacobian in concentration units
98 Eigen::SparseMatrix<double> surfJac = kin->netProductionRates_ddCi();
99 // loop through surface specific jacobian and add elements to triplets vector
100 // accordingly
101 for (int k=0; k<surfJac.outerSize(); ++k) {
102 for (Eigen::SparseMatrix<double>::InnerIterator it(surfJac, k); it; ++it) {
103 size_t row = it.row();
104 size_t col = it.col();
105 auto& rowPhase = kin->speciesPhase(row);
106 auto& colPhase = kin->speciesPhase(col);
107 size_t rpi = kin->phaseIndex(rowPhase.name());
108 size_t cpi = kin->phaseIndex(colPhase.name());
109 // check if the reactor kinetics object contains both phases to avoid
110 // any solid phases which may be included then use phases to map surf
111 // kinetics indicies to reactor kinetic indices
112 if ((rpi == spi || rpi == gpi) && (cpi == spi || cpi == gpi) ) {
113 // subtract start of phase
114 row -= kin->kineticsSpeciesIndex(0, rpi);
115 col -= kin->kineticsSpeciesIndex(0, cpi);
116 // since the gas phase is the first phase in the reactor state
117 // vector add the offset only if it is a surf species index to
118 // both row and col
119 row = (rpi != spi) ? row : row + offset;
120 // determine appropriate scalar to account for dimensionality
121 // gas phase species indices will be less than m_nsp
122 // so use volume if that is the case or area otherwise
123 double scalar = A;
124 if (cpi == spi) {
125 col += offset;
126 scalar /= A;
127 } else {
128 scalar /= m_vol;
129 }
130 // push back scaled value triplet
131 triplets.emplace_back(static_cast<int>(row), static_cast<int>(col),
132 scalar * it.value());
133 }
134 }
135 }
136 // add species in this surface to the offset
137 offset += nk;
138 }
139}
140
142{
143 // Use inverse molecular weights to convert to moles
144 const double* Y = m_thermo->massFractions();
145 const vector<double>& imw = m_thermo->inverseMolecularWeights();
146 for (size_t i = 0; i < m_nsp; i++) {
147 y[i] = m_mass * imw[i] * Y[i];
148 }
149}
150
152{
153 const vector<double>& mw = m_thermo->molecularWeights();
154 // calculate mass from moles
155 m_mass = 0;
156 for (size_t i = 0; i < m_nsp; i++) {
157 m_mass += y[i] * mw[i];
158 }
159}
160
162{
163 if (m_thermo == 0) {
164 throw CanteraError("MoleReactor::getState",
165 "Error: reactor is empty.");
166 }
167 m_thermo->restoreState(m_state);
168 // set the first component to the internal energy
169 m_mass = m_thermo->density() * m_vol;
170 y[0] = m_thermo->intEnergy_mass() * m_mass;
171 // set the second component to the total volume
172 y[1] = m_vol;
173 // set components y+2 ... y+K+2 to the moles of each species
174 getMoles(y + m_sidx);
175 // set the remaining components to the surface species
176 // moles on walls
178}
179
181{
182 // The components of y are [0] total internal energy, [1] the total volume, and
183 // [2...K+3] are the moles of each species, and [K+3...] are the moles
184 // of surface species on each wall.
186 m_vol = y[1];
187 m_thermo->setMolesNoTruncate(y + m_sidx);
188 if (m_energy) {
189 double U = y[0];
190 // Residual function: error in internal energy as a function of T
191 auto u_err = [this, U](double T) {
192 m_thermo->setState_TD(T, m_mass / m_vol);
193 return m_thermo->intEnergy_mass() * m_mass - U;
194 };
195 double T = m_thermo->temperature();
196 boost::uintmax_t maxiter = 100;
197 pair<double, double> TT;
198 try {
199 TT = bmt::bracket_and_solve_root(
200 u_err, T, 1.2, true, bmt::eps_tolerance<double>(48), maxiter);
201 } catch (std::exception&) {
202 // Try full-range bisection if bracketing fails (for example, near
203 // temperature limits for the phase's equation of state)
204 try {
205 TT = bmt::bisect(u_err, m_thermo->minTemp(), m_thermo->maxTemp(),
206 bmt::eps_tolerance<double>(48), maxiter);
207 } catch (std::exception& err2) {
208 // Set m_thermo back to a reasonable state if root finding fails
209 m_thermo->setState_TD(T, m_mass / m_vol);
210 throw CanteraError("MoleReactor::updateState",
211 "{}\nat U = {}, rho = {}", err2.what(), U, m_mass / m_vol);
212 }
213 }
214 if (fabs(TT.first - TT.second) > 1e-7*TT.first) {
215 throw CanteraError("MoleReactor::updateState", "root finding failed");
216 }
217 m_thermo->setState_TD(TT.second, m_mass / m_vol);
218 } else {
219 m_thermo->setDensity(m_mass / m_vol);
220 }
221 updateConnected(true);
223}
224
225void MoleReactor::eval(double time, double* LHS, double* RHS)
226{
227 double* dndt = RHS + m_sidx; // moles per time
228
229 evalWalls(time);
230 m_thermo->restoreState(m_state);
231
232 evalSurfaces(LHS + m_nsp + m_sidx, RHS + m_nsp + m_sidx, m_sdot.data());
233 // inverse molecular weights for conversion
234 const vector<double>& imw = m_thermo->inverseMolecularWeights();
235 // volume equation
236 RHS[1] = m_vdot;
237
238 if (m_chem) {
239 m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
240 }
241
242 // Energy equation.
243 // @f[
244 // \dot U = -P\dot V + A \dot q + \dot m_{in} h_{in} - \dot m_{out} h.
245 // @f]
246 if (m_energy) {
247 RHS[0] = - m_thermo->pressure() * m_vdot + m_Qdot;
248 } else {
249 RHS[0] = 0.0;
250 }
251
252 for (size_t k = 0; k < m_nsp; k++) {
253 // production in gas phase and from surfaces
254 dndt[k] = m_wdot[k] * m_vol + m_sdot[k];
255 }
256
257 // add terms for outlets
258 for (auto outlet : m_outlet) {
259 // flow of species into system and dilution by other species
260 for (size_t n = 0; n < m_nsp; n++) {
261 dndt[n] -= outlet->outletSpeciesMassFlowRate(n) * imw[n];
262 }
263 // energy update based on mass flow
264 double mdot = outlet->massFlowRate();
265 if (m_energy) {
266 RHS[0] -= mdot * m_enthalpy;
267 }
268 }
269
270 // add terms for inlets
271 for (auto inlet : m_inlet) {
272 double mdot = inlet->massFlowRate();
273 for (size_t n = 0; n < m_nsp; n++) {
274 // flow of species into system and dilution by other species
275 dndt[n] += inlet->outletSpeciesMassFlowRate(n) * imw[n];
276 }
277 if (m_energy) {
278 RHS[0] += mdot * inlet->enthalpy_mass();
279 }
280 }
281}
282
283
284size_t MoleReactor::componentIndex(const string& nm) const
285{
286 size_t k = speciesIndex(nm);
287 if (k != npos) {
288 return k + m_sidx;
289 } else if (nm == "int_energy") {
290 return 0;
291 } else if (nm == "volume") {
292 return 1;
293 } else {
294 return npos;
295 }
296}
297
299 if (k == 0) {
300 return "int_energy";
301 } else if (k == 1) {
302 return "volume";
303 } else if (k >= m_sidx && k < neq()) {
304 k -= m_sidx;
305 if (k < m_thermo->nSpecies()) {
306 return m_thermo->speciesName(k);
307 } else {
308 k -= m_thermo->nSpecies();
309 }
310 for (auto& S : m_surfaces) {
311 ThermoPhase* th = S->thermo();
312 if (k < th->nSpecies()) {
313 return th->speciesName(k);
314 } else {
315 k -= th->nSpecies();
316 }
317 }
318 }
319 throw CanteraError("MoleReactor::componentName", "Index is out of bounds.");
320}
321
322}
Base class for kinetics managers and also contains the kineticsmgr module documentation (see Kinetics...
Header file for class ReactorSurface.
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase,...
Base class for exceptions thrown by Cantera classes.
double outletSpeciesMassFlowRate(size_t k)
Mass flow rate (kg/s) of outlet species k.
double enthalpy_mass()
specific enthalpy
double massFlowRate()
Mass flow rate (kg/s).
Definition FlowDevice.h:39
Public interface for kinetics managers.
Definition Kinetics.h:125
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition Kinetics.h:276
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:363
void evalSurfaces(double *LHS, double *RHS, double *sdot) override
Evaluate terms related to surface reactions.
void getSurfaceInitialConditions(double *y) override
Get initial conditions for SurfPhase objects attached to this reactor.
void getMoles(double *y)
Get moles of the system from mass fractions stored by thermo object.
void eval(double t, double *LHS, double *RHS) override
Evaluate the reactor governing equations.
size_t componentIndex(const string &nm) const override
Return the index in the solution vector for this reactor of the component named nm.
const size_t m_sidx
const value for the species start index
Definition MoleReactor.h:62
void setMassFromMoles(double *y)
Set internal mass variable based on moles given.
virtual void addSurfaceJacobian(vector< Eigen::Triplet< double > > &triplets)
For each surface in the reactor, update vector of triplets with all relevant surface jacobian derivat...
void getState(double *y) override
Get the the current state of the reactor.
string componentName(size_t k) override
Return the name of the solution component with index i.
void updateState(double *y) override
Set the state of the reactor to correspond to the state vector y.
void initialize(double t0=0.0) override
Initialize the reactor.
void updateSurfaceState(double *y) override
Update the state of SurfPhase objects attached to this reactor.
void restoreState(const vector< double > &state)
Restore a state saved on a previous call to saveState.
Definition Phase.cpp:260
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:231
void setState_TD(double t, double rho)
Set the internally stored temperature (K) and density (kg/m^3)
Definition Phase.cpp:377
double temperature() const
Temperature (K).
Definition Phase.h:562
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
virtual void setDensity(const double density_)
Set the internally stored density (kg/m^3) of the phase.
Definition Phase.cpp:586
const vector< double > & inverseMolecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition Phase.cpp:400
const double * massFractions() const
Return a const pointer to the mass fraction array.
Definition Phase.h:442
const vector< double > & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition Phase.cpp:395
virtual double density() const
Density (kg/m^3).
Definition Phase.h:587
virtual void setMolesNoTruncate(const double *const N)
Set the state of the object with moles in [kmol].
Definition Phase.cpp:528
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition Phase.h:580
FlowDevice & outlet(size_t n=0)
Return a reference to the n-th outlet FlowDevice connected to this reactor.
FlowDevice & inlet(size_t n=0)
Return a reference to the n-th inlet FlowDevice connected to this reactor.
double m_vol
Current volume of the reactor [m^3].
size_t m_nsp
Number of homogeneous species in the mixture.
double m_enthalpy
Current specific enthalpy of the reactor [J/kg].
Kinetics * m_kin
Pointer to the homogeneous Kinetics object that handles the reactions.
Definition Reactor.h:283
vector< double > m_wdot
Species net molar production rates.
Definition Reactor.h:295
virtual void evalWalls(double t)
Evaluate terms related to Walls.
Definition Reactor.cpp:278
size_t neq()
Number of equations (state variables) for this reactor.
Definition Reactor.h:109
double m_Qdot
net heat transfer into the reactor, through walls [W]
Definition Reactor.h:287
double m_mass
total mass
Definition Reactor.h:289
vector< double > m_sdot
Production rates of gas phase species on surfaces [kmol/s].
Definition Reactor.h:293
double m_vdot
net rate of volume change from moving walls [m^3/s]
Definition Reactor.h:285
void initialize(double t0=0.0) override
Initialize the reactor.
Definition Reactor.cpp:87
virtual size_t speciesIndex(const 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:429
virtual void updateConnected(bool updatePressure)
Update the state information needed by connected reactors, flow devices, and reactor walls.
Definition Reactor.cpp:187
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition SurfPhase.h:98
double size(size_t k) const
Returns the number of sites occupied by one molecule of species k.
Definition SurfPhase.h:221
Base class for a phase with thermodynamic properties.
virtual double minTemp(size_t k=npos) const
Minimum temperature for which the thermodynamic data for the species or phase are valid.
virtual double maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
double intEnergy_mass() const
Specific internal energy. Units: J/kg.
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
offset
Offsets of solution components in the 1D solution array.
Definition StFlow.h:24
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...