Cantera  3.0.0
Loading...
Searching...
No Matches
FlowReactor.cpp
Go to the documentation of this file.
1//! @file FlowReactor.cpp A steady-state, ideal-gas, adiabatic,
2//! constant-area (cylindrical), frictionless plug flow reactor
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
16
17namespace Cantera
18{
19
20void FlowReactor::getStateDae(double* y, double* ydot)
21{
22 if (m_thermo == nullptr) {
23 throw CanteraError("FlowReactor::getStateDae", "Error: reactor is empty.");
24 }
25 m_thermo->restoreState(m_state);
26 m_thermo->getMassFractions(y+m_offset_Y);
27 const vector<double>& mw = m_thermo->molecularWeights();
28
29 // set the first component to the initial density
30 y[0] = m_rho;
31
32 if (m_u < 0) {
33 throw CanteraError("FlowReactor::getStateDae",
34 "Set mass flow rate before initializing reactor");
35 }
36
37 // set the second component to the initial speed
38 y[1] = m_u;
39
40 // set the third component to the initial pressure
41 y[2] = m_P;
42
43 // set the fourth component to the initial temperature
44 y[3] = m_T;
45
46 if (m_chem) {
47 m_kin->getNetProductionRates(m_wdot.data()); // "omega dot"
48 }
49
50 // need to advance the reactor surface to steady state to get the initial
51 // coverages
52 for (auto m_surf : m_surfaces) {
53 m_surf->syncState();
54 auto kin = static_cast<InterfaceKinetics*>(m_surf->kinetics());
57 auto& surf = dynamic_cast<SurfPhase&>(kin->thermo(kin->reactionPhaseIndex()));
58 vector<double> cov(surf.nSpecies());
59 surf.getCoverages(cov.data());
60 m_surf->setCoverages(cov.data());
61 }
62
63 // set the initial coverages
65
66 // reset ydot vector
67 std::fill(ydot, ydot + m_nv, 0.0);
68 // calculate initial d(coverage)/dt values
69 evalSurfaces(ydot + m_nsp + 4, m_sdot.data());
70
71 // next, we must solve for the initial derivative values
72 // a . ydot = b
73 // ydot -> [rho', u', P', T', Yk']
74 // a -> coefficients of [rho', u', P', T', Yk'] in the governing eq's, given
75 // the initial state
76 // b -> rhs constant of each conservation equation
77 //
78 // note that the species coverages are included in the algebraic constraints,
79 // hence are not included here
81
82 // first row is the ideal gas equation, differentiated with respect to z
83 a(0, 0) = - GasConstant * m_T / m_thermo->meanMolecularWeight();
84 a(0, 2) = 1;
85 a(0, 3) = - m_rho * GasConstant / m_thermo->meanMolecularWeight();
86 for (size_t i = 0; i < m_nsp; ++i) {
87 a(0, m_offset_Y + i) = - m_rho * m_T / mw[i];
88 }
89
90 // initialize the second row from mass conservation equation (Kee 16.48)
91 a(1, 0) = m_u; // rho' * u
92 a(1, 1) = m_rho; // u' * rho
93
94 // initialize the third row from momentum conservation equation (Kee 16.62),
95 // rho * u * u' + P' = -u * (P / A_c) * sum(sdot_k * W_k)
96 a(2, 1) = m_rho * m_u; // rho * u * u'
97 a(2, 2) = 1; // 1 * P'
98
99 // initialize the fourth row from conservation of energy (Kee 16.58), adiabatic
100 double cp_mass = m_thermo->cp_mass();
101 a(3, 3) = m_rho * m_u * cp_mass; // rho * u * cp * T'
102
103 // initialize the next rows from the mass-fraction equations (Kee 16.51)
104 for (size_t i = 0; i < m_nsp; ++i) {
105 a(m_offset_Y + i, m_offset_Y + i) = m_rho * m_u; // rho * u * Yk'
106 }
107
108 // now set the RHS vector
109
110 // get (perim / Ac) * sum(sk' * wk), used multiple places
111 double h_sk_wk = 0;
112 double hydraulic = surfaceAreaToVolumeRatio();
113 for (size_t i = 0; i < m_nsp; ++i) {
114 h_sk_wk += hydraulic * m_sdot[i] * mw[i];
115 }
116
117 // RHS of ideal gas eq. is zero
118 ydot[0] = 0;
119
120 // mass continuity, Kee 16.48
121 // (perim / Ac) * sum(sk' * wk)
122 ydot[1] = h_sk_wk;
123
124 // momentum conservation, Kee 16.62, no friction
125 // - u * (perim / Ac) * sum(sk' * wk)
126 ydot[2] = -m_u * h_sk_wk;
127
128 if (m_energy) {
129 // energy conservation, Kee 16.58, adiabatic
130 // -sum(wk' * Wk * hk) - (perim / Ac) * sum(sk' * Wk * hk)
131 // Note: the partial molar enthalpies are in molar units, while Kee uses mass
132 // units, where:
133 // h_mass = h_mole / Wk
134 // hence:
135 // h_mass * Wk = h_mol
136 m_thermo->getPartialMolarEnthalpies(m_hk.data());
137 for (size_t i = 0; i < m_nsp; ++i) {
138 ydot[3] -= m_hk[i] * (m_wdot[i] + hydraulic * m_sdot[i]);
139 }
140 } else {
141 ydot[3] = 0;
142 }
143
144 // mass-fraction equations Kee 16.51
145 // - Yk * (perim / Ac) * sum(sk' * Wk) + wk' * Wk + (perim / Ac) * sk' * Wk
146 for (size_t i = 0; i < m_nsp; ++i) {
147 ydot[m_offset_Y + i] = -y[m_offset_Y + i] * h_sk_wk +
148 mw[i] * (m_wdot[i] + hydraulic * m_sdot[i]);
149 }
150
151 // and solve
152 solve(a, ydot, 1, 0);
153}
154
156{
158 m_thermo->restoreState(m_state);
159 // initialize state
160 m_T = m_thermo->temperature();
161 m_rho = m_thermo->density();
162 m_P = m_thermo->pressure();
163 m_T = m_thermo->temperature();
164 // resize temporary arrays
165 m_wdot.resize(m_nsp);
166 m_hk.resize(m_nsp);
167 // set number of variables to the number of non-species equations
168 // i.e., density, velocity, pressure and temperature
169 // plus the number of species in the gas phase
170 m_nv = m_offset_Y + m_nsp;
171 if (m_surfaces.size()) {
172 size_t n_surf_species = 0;
173 size_t n_total_species = 0;
174 for (auto const &m_surf : m_surfaces) {
175 // finally, add the number of species the surface for the site-fraction
176 // equations
177 n_surf_species += m_surf->thermo()->nSpecies();
178 n_total_species += m_surf->kinetics()->nTotalSpecies();
179 }
180 m_nv += n_surf_species;
181 m_sdot_temp.resize(n_total_species);
182 }
183}
184
186{
188 m_rho = m_thermo->density();
189 m_P = m_thermo->pressure();
190 m_T = m_thermo->temperature();
191}
192
194{
195 // Set the mass fractions and density of the mixture.
196 m_rho = y[0];
197 m_u = y[1];
198 m_P = y[2];
199 m_T = y[3];
200 double* mss = y + m_offset_Y;
201 m_thermo->setMassFractions_NoNorm(mss);
202 m_thermo->setState_TP(m_T, m_P);
203
204 // update surface
206
207 m_thermo->saveState(m_state);
208}
209
211{
212 m_rho = m_thermo->density();
213 m_u = mdot/(m_rho * m_area);
214}
215
217{
218 warn_deprecated("FlowReactor::distance", "To be removed after Cantera 3.0."
219 "Access distance through the ReactorNet object.");
220 if (m_net != nullptr) {
221 return m_net->distance();
222 } else {
223 return 0.0;
224 }
225}
226
227void FlowReactor::setArea(double area) {
228 double mdot = m_rho * m_u * m_area;
229 m_area = area;
230 setMassFlowRate(mdot);
231}
232
234 if (m_sa_to_vol > 0) {
235 return m_sa_to_vol;
236 }
237
238 // assuming a cylinder, volume = Pi * r^2 * L, and perimeter = 2 * Pi * r * L
239 // where L is the length, and r is the radius of the reactor
240 // hence, perimeter / area = 2 * Pi * r * L / (Pi * L * r^2) = 2 / r
241 return 2.0 / sqrt(m_area / Pi);
242}
243
245{
246 size_t loc = 0;
247 for (auto& S : m_surfaces) {
248 // set the coverages without normalization to avoid over-constraining
249 // the system.
250 // note: the ReactorSurface class doesn't normalize when calling setCoverages
251 S->setCoverages(y+loc);
252 S->syncState();
253 loc += S->thermo()->nSpecies();
254 }
255}
256
257void FlowReactor::evalDae(double time, double* y, double* ydot, double* residual)
258{
259 m_thermo->restoreState(m_state);
260
261 evalSurfaces(ydot + m_nsp + 4, m_sdot.data());
262 const vector<double>& mw = m_thermo->molecularWeights();
263 double sk_wk = 0;
264 for (size_t i = 0; i < m_nsp; ++i) {
265 sk_wk = m_sdot[i] * mw[i];
266 }
267 m_thermo->getPartialMolarEnthalpies(m_hk.data());
268 // get net production
269 if (m_chem) {
271 }
272
273 // set dphi/dz variables
274 double drhodz = ydot[0];
275 double dudz = ydot[1];
276 double dPdz = ydot[2];
277 double dTdz = ydot[3];
278
279 // use equation of state for density residual
280 residual[0] = m_rho - m_thermo->density();
281
282 //! use mass continuity for velocity residual
283 //! Kee.'s Chemically Reacting Flow, Eq. 16.48
284 double hydraulic = surfaceAreaToVolumeRatio();
285 residual[1] = m_u * drhodz + m_rho * dudz - sk_wk * hydraulic;
286
287 //! Use conservation of momentum for pressure residual
288 //! Kee.'s Chemically Reacting Flow, Eq. 16.62
289 //! [NOTE: friction is neglected]
290 residual[2] = m_rho * m_u * dudz + m_u * hydraulic * sk_wk + dPdz;
291
292 //! use conservation of energy for temperature residual
293 //! Kee.'s Chemically Reacting Flow, Eq. 16.58
294 //! [NOTE: adiabatic]
295 // Also, as above:
296 // h_mass = h_mole / Wk
297 // hence:
298 // h_mass * Wk = h_mol
299 if (m_energy) {
300 double cp_mass = m_thermo->cp_mass();
301 residual[3] = m_rho * m_u * cp_mass * dTdz;
302 for (size_t i = 0; i < m_nsp; ++i) {
303 residual[3] += m_hk[i] * (m_wdot[i] + hydraulic * m_sdot[i]);
304 }
305 } else {
306 residual[3] = dTdz;
307 }
308
309 //! species conservation equations
310 //! Kee.'s Chemically Reacting Flow, Eq. 16.51
311 for (size_t i = 0; i < m_nsp; ++i) {
312 residual[i + m_offset_Y] = m_rho * m_u * ydot[i + m_offset_Y] +
313 y[i + m_offset_Y] * hydraulic * sk_wk -
314 mw[i] * (m_wdot[i] + hydraulic * m_sdot[i]);
315 }
316
317 // surface algebraic constraints
318 if (m_surfaces.size()) {
319 size_t loc = m_offset_Y + m_nsp; // offset into residual vector
320 for (auto &m_surf : m_surfaces) {
321 Kinetics* kin = m_surf->kinetics();
322 size_t nk = m_surf->thermo()->nSpecies();
324 size_t ns = kin->reactionPhaseIndex();
325 size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
326 double sum = y[loc];
327 for (size_t i = 1; i < nk; ++i) {
328 //! net surface production rate residuals
329 //! Kee.'s Chemically Reacting Flow, Eq. 16.63
330 residual[loc + i] = m_sdot_temp[surfloc + i];
331 //! next, evaluate the algebraic constraint to explicitly conserve
332 //! surface site fractions.
333 //! Kee.'s Chemically Reacting Flow, Eq. 16.64
334 sum += y[loc + i];
335 }
336 // set residual of the first species sum - 1 to ensure
337 // surface site conservation.
338 // Note: the first species is selected to be consistent with
339 // Reactor::evalSurfaces
340 residual[loc] = sum - 1.0;
341 loc += nk;
342 }
343 }
344}
345
346void FlowReactor::getConstraints(double* constraints) {
347 // mark all variables differential equations unless otherwise specified
348 std::fill(constraints, constraints + m_nv, 1.0);
349 // the species coverages are algebraic constraints
350 std::fill(constraints + m_offset_Y + m_nsp, constraints + m_nv, 0.0);
351}
352
353size_t FlowReactor::componentIndex(const string& nm) const
354{
355 size_t k = speciesIndex(nm);
356 if (k != npos) {
357 return k + m_offset_Y;
358 } else if (nm == "density") {
359 return 0;
360 } else if (nm == "speed") {
361 return 1;
362 } else if (nm == "pressure") {
363 return 2;
364 } else if (nm == "temperature") {
365 return 3;
366 } else {
367 return npos;
368 }
369}
370
372{
373 if (k == 0) {
374 return "density";
375 } else if (k == 1) {
376 return "speed";
377 } else if (k == 2) {
378 return "pressure";
379 } else if (k == 3) {
380 return "temperature";
381 } else if (k >= 4 && k < neq()) {
382 k -= 4;
383 if (k < m_thermo->nSpecies()) {
384 return m_thermo->speciesName(k);
385 } else {
386 k -= m_thermo->nSpecies();
387 }
388 for (auto& S : m_surfaces) {
389 ThermoPhase* th = S->thermo();
390 if (k < th->nSpecies()) {
391 return th->speciesName(k);
392 } else {
393 k -= th->nSpecies();
394 }
395 }
396 }
397 throw CanteraError("FlowReactor::componentName", "Index {} is out of bounds.", k);
398}
399
400}
Headers for the DenseMatrix object, which deals with dense rectangular matrices and description of th...
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,...
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
Base class for exceptions thrown by Cantera classes.
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition DenseMatrix.h:55
double m_area
reactor area [m^2]
vector< double > m_sdot_temp
temporary storage for surface species production rates
const size_t m_offset_Y
offset to the species equations
double area() const
The cross-sectional area of the reactor [m^2].
Definition FlowReactor.h:61
double surfaceAreaToVolumeRatio() const
The ratio of the reactor's surface area to volume ratio [m^-1].
void setMassFlowRate(double mdot)
Set the mass flow rate through the reactor [kg/s].
double m_ss_rtol
steady-state relative tolerance, used to determine initial surface coverages
int m_max_ss_error_fails
maximum number of steady-state integrator error test failures
void getConstraints(double *constraints) override
Given a vector of length neq(), mark which variables should be considered algebraic constraints.
double m_u
Axial velocity [m/s]. Second component of the state vector.
size_t componentIndex(const string &nm) const override
Return the index in the solution vector for this reactor of the component named nm.
double distance() const
void getStateDae(double *y, double *ydot) override
Get the current state and derivative vector of the reactor for a DAE solver.
int m_max_ss_steps
maximum number of steady-state coverage integrator-steps
double m_rho
Density [kg/m^3]. First component of the state vector.
void evalDae(double t, double *y, double *ydot, double *residual) override
Evaluate the reactor governing equations.
string componentName(size_t k) override
Return the name of the solution component with index i.
void syncState() override
Set the state of the reactor to correspond to the state of the associated ThermoPhase object.
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.
double m_P
Pressure [Pa]. Third component of the state vector.
double m_T
Temperature [K]. Fourth component of the state vector.
double m_sa_to_vol
reactor surface area to volume ratio [m^-1]
double m_ss_atol
steady-state absolute tolerance, used to determine initial surface coverages
vector< double > m_hk
temporary storage for species partial molar enthalpies
void setArea(double area)
Sets the area of the reactor [m^2].
A kinetics manager for heterogeneous reaction mechanisms.
void advanceCoverages(double tstep, double rtol=1.e-7, double atol=1.e-14, double maxStepSize=0, size_t maxSteps=20000, size_t maxErrTestFails=7)
Advance the surface coverages in time.
Public interface for kinetics managers.
Definition Kinetics.h:126
size_t reactionPhaseIndex() const
Phase where the reactions occur.
Definition Kinetics.h:235
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
Definition Kinetics.h:254
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition Kinetics.h:288
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:435
void restoreState(const vector< double > &state)
Restore a state saved on a previous call to saveState.
Definition Phase.cpp:275
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:245
virtual void setMassFractions_NoNorm(const double *const y)
Set the mass fractions to the specified values without normalizing.
Definition Phase.cpp:370
double temperature() const
Temperature (K).
Definition Phase.h:662
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:760
void saveState(vector< double > &state) const
Save the current internal state of the phase.
Definition Phase.cpp:251
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:151
const vector< double > & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition Phase.cpp:501
virtual double density() const
Density (kg/m^3).
Definition Phase.h:687
void getMassFractions(double *const y) const
Get the species mass fractions.
Definition Phase.cpp:584
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition Phase.h:680
ReactorNet * m_net
The ReactorNet that this reactor is part of.
virtual void syncState()
Set the state of the reactor to correspond to the state of the associated ThermoPhase object.
size_t m_nsp
Number of homogeneous species in the mixture.
double distance()
Current position [m] along the length of the reactor network, for reactors that are solved as a funct...
virtual void evalSurfaces(double *LHS, double *RHS, double *sdot)
Evaluate terms related to surface reactions.
Definition Reactor.cpp:287
Kinetics * m_kin
Pointer to the homogeneous Kinetics object that handles the reactions.
Definition Reactor.h:277
vector< double > m_wdot
Species net molar production rates.
Definition Reactor.h:289
size_t neq()
Number of equations (state variables) for this reactor.
Definition Reactor.h:103
vector< double > m_sdot
Production rates of gas phase species on surfaces [kmol/s].
Definition Reactor.h:287
virtual void getSurfaceInitialConditions(double *y)
Get initial conditions for SurfPhase objects attached to this reactor.
Definition Reactor.cpp:75
void initialize(double t0=0.0) override
Initialize the reactor.
Definition Reactor.cpp:84
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:430
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition SurfPhase.h:98
Base class for a phase with thermodynamic properties.
virtual void getPartialMolarEnthalpies(double *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
virtual void setState_TP(double t, double p)
Set the temperature (K) and pressure (Pa)
double cp_mass() const
Specific heat at constant pressure. Units: J/kg/K.
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition ct_defs.h:120
const double Pi
Pi.
Definition ct_defs.h:68
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:195
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
Definition AnyMap.cpp:1926