Cantera  3.1.0a1
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(0));
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
216void FlowReactor::setArea(double area) {
217 double mdot = m_rho * m_u * m_area;
218 m_area = area;
219 setMassFlowRate(mdot);
220}
221
223 if (m_sa_to_vol > 0) {
224 return m_sa_to_vol;
225 }
226
227 // assuming a cylinder, volume = Pi * r^2 * L, and perimeter = 2 * Pi * r * L
228 // where L is the length, and r is the radius of the reactor
229 // hence, perimeter / area = 2 * Pi * r * L / (Pi * L * r^2) = 2 / r
230 return 2.0 / sqrt(m_area / Pi);
231}
232
234{
235 size_t loc = 0;
236 for (auto& S : m_surfaces) {
237 // set the coverages without normalization to avoid over-constraining
238 // the system.
239 // note: the ReactorSurface class doesn't normalize when calling setCoverages
240 S->setCoverages(y+loc);
241 S->syncState();
242 loc += S->thermo()->nSpecies();
243 }
244}
245
246void FlowReactor::evalDae(double time, double* y, double* ydot, double* residual)
247{
248 m_thermo->restoreState(m_state);
249
250 evalSurfaces(ydot + m_nsp + 4, m_sdot.data());
251 const vector<double>& mw = m_thermo->molecularWeights();
252 double sk_wk = 0;
253 for (size_t i = 0; i < m_nsp; ++i) {
254 sk_wk = m_sdot[i] * mw[i];
255 }
256 m_thermo->getPartialMolarEnthalpies(m_hk.data());
257 // get net production
258 if (m_chem) {
260 }
261
262 // set dphi/dz variables
263 double drhodz = ydot[0];
264 double dudz = ydot[1];
265 double dPdz = ydot[2];
266 double dTdz = ydot[3];
267
268 // use equation of state for density residual
269 residual[0] = m_rho - m_thermo->density();
270
271 //! use mass continuity for velocity residual
272 //! Kee.'s Chemically Reacting Flow, Eq. 16.48
273 double hydraulic = surfaceAreaToVolumeRatio();
274 residual[1] = m_u * drhodz + m_rho * dudz - sk_wk * hydraulic;
275
276 //! Use conservation of momentum for pressure residual
277 //! Kee.'s Chemically Reacting Flow, Eq. 16.62
278 //! [NOTE: friction is neglected]
279 residual[2] = m_rho * m_u * dudz + m_u * hydraulic * sk_wk + dPdz;
280
281 //! use conservation of energy for temperature residual
282 //! Kee.'s Chemically Reacting Flow, Eq. 16.58
283 //! [NOTE: adiabatic]
284 // Also, as above:
285 // h_mass = h_mole / Wk
286 // hence:
287 // h_mass * Wk = h_mol
288 if (m_energy) {
289 double cp_mass = m_thermo->cp_mass();
290 residual[3] = m_rho * m_u * cp_mass * dTdz;
291 for (size_t i = 0; i < m_nsp; ++i) {
292 residual[3] += m_hk[i] * (m_wdot[i] + hydraulic * m_sdot[i]);
293 }
294 } else {
295 residual[3] = dTdz;
296 }
297
298 //! species conservation equations
299 //! Kee.'s Chemically Reacting Flow, Eq. 16.51
300 for (size_t i = 0; i < m_nsp; ++i) {
301 residual[i + m_offset_Y] = m_rho * m_u * ydot[i + m_offset_Y] +
302 y[i + m_offset_Y] * hydraulic * sk_wk -
303 mw[i] * (m_wdot[i] + hydraulic * m_sdot[i]);
304 }
305
306 // surface algebraic constraints
307 if (m_surfaces.size()) {
308 size_t loc = m_offset_Y + m_nsp; // offset into residual vector
309 for (auto &m_surf : m_surfaces) {
310 Kinetics* kin = m_surf->kinetics();
311 size_t nk = m_surf->thermo()->nSpecies();
313 double sum = y[loc];
314 for (size_t i = 1; i < nk; ++i) {
315 //! net surface production rate residuals
316 //! Kee.'s Chemically Reacting Flow, Eq. 16.63
317 residual[loc + i] = m_sdot_temp[i];
318 //! next, evaluate the algebraic constraint to explicitly conserve
319 //! surface site fractions.
320 //! Kee.'s Chemically Reacting Flow, Eq. 16.64
321 sum += y[loc + i];
322 }
323 // set residual of the first species sum - 1 to ensure
324 // surface site conservation.
325 // Note: the first species is selected to be consistent with
326 // Reactor::evalSurfaces
327 residual[loc] = sum - 1.0;
328 loc += nk;
329 }
330 }
331}
332
333void FlowReactor::getConstraints(double* constraints) {
334 // mark all variables differential equations unless otherwise specified
335 std::fill(constraints, constraints + m_nv, 1.0);
336 // the species coverages are algebraic constraints
337 std::fill(constraints + m_offset_Y + m_nsp, constraints + m_nv, 0.0);
338}
339
340size_t FlowReactor::componentIndex(const string& nm) const
341{
342 size_t k = speciesIndex(nm);
343 if (k != npos) {
344 return k + m_offset_Y;
345 } else if (nm == "density") {
346 return 0;
347 } else if (nm == "speed") {
348 return 1;
349 } else if (nm == "pressure") {
350 return 2;
351 } else if (nm == "temperature") {
352 return 3;
353 } else {
354 return npos;
355 }
356}
357
359{
360 if (k == 0) {
361 return "density";
362 } else if (k == 1) {
363 return "speed";
364 } else if (k == 2) {
365 return "pressure";
366 } else if (k == 3) {
367 return "temperature";
368 } else if (k >= 4 && k < neq()) {
369 k -= 4;
370 if (k < m_thermo->nSpecies()) {
371 return m_thermo->speciesName(k);
372 } else {
373 k -= m_thermo->nSpecies();
374 }
375 for (auto& S : m_surfaces) {
376 ThermoPhase* th = S->thermo();
377 if (k < th->nSpecies()) {
378 return th->speciesName(k);
379 } else {
380 k -= th->nSpecies();
381 }
382 }
383 }
384 throw CanteraError("FlowReactor::componentName", "Index {} is out of bounds.", k);
385}
386
387}
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.
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:125
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
Definition Kinetics.h:242
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:363
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
virtual void setMassFractions_NoNorm(const double *const y)
Set the mass fractions to the specified values without normalizing.
Definition Phase.cpp:355
double temperature() const
Temperature (K).
Definition Phase.h:562
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:655
void saveState(vector< double > &state) const
Save the current internal state of the phase.
Definition Phase.cpp:236
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
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
void getMassFractions(double *const y) const
Get the species mass fractions.
Definition Phase.cpp:471
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition Phase.h:580
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.
virtual void evalSurfaces(double *LHS, double *RHS, double *sdot)
Evaluate terms related to surface reactions.
Definition Reactor.cpp:290
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
size_t neq()
Number of equations (state variables) for this reactor.
Definition Reactor.h:109
vector< double > m_sdot
Production rates of gas phase species on surfaces [kmol/s].
Definition Reactor.h:293
virtual void getSurfaceInitialConditions(double *y)
Get initial conditions for SurfPhase objects attached to this reactor.
Definition Reactor.cpp:78
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
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:180
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.