22 if (m_thermo ==
nullptr) {
23 throw CanteraError(
"FlowReactor::getStateDae",
"Error: reactor is empty.");
34 "Set mass flow rate before initializing reactor");
52 for (
auto m_surf : m_surfaces) {
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());
67 std::fill(ydot, ydot + m_nv, 0.0);
86 for (
size_t i = 0; i <
m_nsp; ++i) {
100 double cp_mass = m_thermo->
cp_mass();
104 for (
size_t i = 0; i <
m_nsp; ++i) {
113 for (
size_t i = 0; i <
m_nsp; ++i) {
114 h_sk_wk += hydraulic *
m_sdot[i] * mw[i];
126 ydot[2] = -
m_u * h_sk_wk;
137 for (
size_t i = 0; i <
m_nsp; ++i) {
146 for (
size_t i = 0; i <
m_nsp; ++i) {
152 solve(a, ydot, 1, 0);
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) {
177 n_surf_species += m_surf->thermo()->nSpecies();
178 n_total_species += m_surf->kinetics()->nTotalSpecies();
180 m_nv += n_surf_species;
218 warn_deprecated(
"FlowReactor::distance",
"To be removed after Cantera 3.0."
219 "Access distance through the ReactorNet object.");
220 if (
m_net !=
nullptr) {
247 for (
auto& S : m_surfaces) {
251 S->setCoverages(y+loc);
253 loc += S->thermo()->nSpecies();
264 for (
size_t i = 0; i <
m_nsp; ++i) {
265 sk_wk =
m_sdot[i] * mw[i];
274 double drhodz = ydot[0];
275 double dudz = ydot[1];
276 double dPdz = ydot[2];
277 double dTdz = ydot[3];
285 residual[1] =
m_u * drhodz +
m_rho * dudz - sk_wk * hydraulic;
290 residual[2] =
m_rho *
m_u * dudz +
m_u * hydraulic * sk_wk + dPdz;
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) {
311 for (
size_t i = 0; i <
m_nsp; ++i) {
318 if (m_surfaces.size()) {
320 for (
auto &m_surf : m_surfaces) {
327 for (
size_t i = 1; i < nk; ++i) {
340 residual[loc] = sum - 1.0;
348 std::fill(constraints, constraints + m_nv, 1.0);
358 }
else if (nm ==
"density") {
360 }
else if (nm ==
"speed") {
362 }
else if (nm ==
"pressure") {
364 }
else if (nm ==
"temperature") {
380 return "temperature";
381 }
else if (k >= 4 && k <
neq()) {
383 if (k < m_thermo->nSpecies()) {
388 for (
auto& S : m_surfaces) {
390 if (k < th->nSpecies()) {
397 throw CanteraError(
"FlowReactor::componentName",
"Index {} is out of bounds.", k);
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...
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].
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.
size_t reactionPhaseIndex() const
Phase where the reactions occur.
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
void restoreState(const vector< double > &state)
Restore a state saved on a previous call to saveState.
size_t nSpecies() const
Returns the number of species in the phase.
virtual void setMassFractions_NoNorm(const double *const y)
Set the mass fractions to the specified values without normalizing.
double temperature() const
Temperature (K).
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
void saveState(vector< double > &state) const
Save the current internal state of the phase.
string speciesName(size_t k) const
Name of the species with index k.
const vector< double > & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
virtual double density() const
Density (kg/m^3).
void getMassFractions(double *const y) const
Get the species mass fractions.
virtual double pressure() const
Return the thermodynamic pressure (Pa).
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.
Kinetics * m_kin
Pointer to the homogeneous Kinetics object that handles the reactions.
vector< double > m_wdot
Species net molar production rates.
size_t neq()
Number of equations (state variables) for this reactor.
vector< double > m_sdot
Production rates of gas phase species on surfaces [kmol/s].
virtual void getSurfaceInitialConditions(double *y)
Get initial conditions for SurfPhase objects attached to this reactor.
void initialize(double t0=0.0) override
Initialize the reactor.
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...
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
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].
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
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.