Cantera  3.2.0a2
Loading...
Searching...
No Matches
ConstPressureReactor.cpp
Go to the documentation of this file.
1//! @file ConstPressureReactor.cpp A constant pressure zero-dimensional reactor
2
3// This file is part of Cantera. See License.txt in the top-level directory or
4// at https://cantera.org/license.txt for license and copyright information.
5
13
14namespace Cantera
15{
16
18{
19 if (m_thermo == 0) {
20 throw CanteraError("ConstPressureReactor::getState",
21 "Error: reactor is empty.");
22 }
23 m_thermo->restoreState(m_state);
24
25 // set the first component to the total mass
26 y[0] = m_thermo->density() * m_vol;
27
28 // set the second component to the total enthalpy
29 y[1] = m_thermo->enthalpy_mass() * m_thermo->density() * m_vol;
30
31 // set components y+2 ... y+K+1 to the mass fractions Y_k of each species
32 m_thermo->getMassFractions(y+2);
33
34 // set the remaining components to the surface species
35 // coverages on the walls
37}
38
40{
42 m_nv -= 1; // Constant pressure reactor has one fewer state variable
43}
44
46{
47 // The components of y are [0] the total mass, [1] the total enthalpy,
48 // [2...K+2) are the mass fractions of each species, and [K+2...] are the
49 // coverages of surface species on each wall.
50 m_mass = y[0];
51 m_thermo->setMassFractions_NoNorm(y+2);
52 if (m_energy) {
53 m_thermo->setState_HP(y[1]/m_mass, m_pressure);
54 } else {
55 m_thermo->setPressure(m_pressure);
56 }
57 m_vol = m_mass / m_thermo->density();
58 updateConnected(false);
60}
61
62void ConstPressureReactor::eval(double time, double* LHS, double* RHS)
63{
64 double& dmdt = RHS[0];
65 double* mdYdt = RHS + 2; // mass * dY/dt
66
67 dmdt = 0.0;
68
69 evalWalls(time);
70
71 m_thermo->restoreState(m_state);
72 const vector<double>& mw = m_thermo->molecularWeights();
73 const double* Y = m_thermo->massFractions();
74
75 evalSurfaces(LHS + m_nsp + 2, RHS + m_nsp + 2, m_sdot.data());
76 double mdot_surf = dot(m_sdot.begin(), m_sdot.end(), mw.begin());
77 dmdt += mdot_surf;
78
79 if (m_chem) {
80 m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
81 }
82
83 for (size_t k = 0; k < m_nsp; k++) {
84 // production in gas phase and from surfaces
85 mdYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k];
86 // dilution by net surface mass flux
87 mdYdt[k] -= Y[k] * mdot_surf;
88 //Assign left-hand side of dYdt ODE as total mass
89 LHS[k+2] = m_mass;
90 }
91
92 // external heat transfer
93 double dHdt = m_Qdot;
94
95 // add terms for outlets
96 for (auto outlet : m_outlet) {
97 double mdot = outlet->massFlowRate();
98 dmdt -= mdot;
99 dHdt -= mdot * m_enthalpy;
100 }
101
102 // add terms for inlets
103 for (auto inlet : m_inlet) {
104 double mdot = inlet->massFlowRate();
105 dmdt += mdot; // mass flow into system
106 for (size_t n = 0; n < m_nsp; n++) {
107 double mdot_spec = inlet->outletSpeciesMassFlowRate(n);
108 // flow of species into system and dilution by other species
109 mdYdt[n] += mdot_spec - mdot * Y[n];
110 }
111 dHdt += mdot * inlet->enthalpy_mass();
112 }
113
114 if (m_energy) {
115 RHS[1] = dHdt;
116 } else {
117 RHS[1] = 0.0;
118 }
119}
120
122 if (nSurfs() != 0) {
123 throw CanteraError("ConstPressureReactor::steadyConstraints",
124 "Steady state solver cannot currently be used with ConstPressureReactor"
125 " when reactor surfaces are present.\n"
126 "See https://github.com/Cantera/enhancements/issues/234");
127 }
128 return {0}; // mass
129}
130
131size_t ConstPressureReactor::componentIndex(const string& nm) const
132{
133 size_t k = speciesIndex(nm);
134 if (k != npos) {
135 return k + 2;
136 } else if (nm == "mass") {
137 return 0;
138 } else if (nm == "enthalpy") {
139 return 1;
140 } else {
141 return npos;
142 }
143}
144
146 if (k == 0) {
147 return "mass";
148 } else if (k == 1) {
149 return "enthalpy";
150 } else if (k >= 2 && k < neq()) {
151 k -= 2;
152 if (k < m_thermo->nSpecies()) {
153 return m_thermo->speciesName(k);
154 } else {
155 k -= m_thermo->nSpecies();
156 }
157 for (auto& S : m_surfaces) {
158 ThermoPhase* th = S->thermo();
159 if (k < th->nSpecies()) {
160 return th->speciesName(k);
161 } else {
162 k -= th->nSpecies();
163 }
164 }
165 }
166 throw CanteraError("ConstPressureReactor::componentName",
167 "Index is out of bounds.");
168}
169
170double ConstPressureReactor::upperBound(size_t k) const {
171 if (k == 0) {
172 return BigNumber; // mass
173 } else if (k == 1) {
174 return BigNumber; // enthalpy
175 } else if (k >= 2 && k < m_nv) {
176 return 1.0; // species mass fraction or surface coverage
177 } else {
178 throw CanteraError("ConstPressureReactor::upperBound",
179 "Index {} is out of bounds.", k);
180 }
181}
182
183double ConstPressureReactor::lowerBound(size_t k) const {
184 if (k == 0) {
185 return 0; // mass
186 } else if (k == 1) {
187 return -BigNumber; // enthalpy
188 } else if (k >= 2 && k < m_nv) {
189 return -Tiny; // species mass fraction or surface coverage
190 } else {
191 throw CanteraError("ConstPressureReactor::lowerBound",
192 "Index {} is out of bounds.", k);
193 }
194}
195
197 for (size_t k = 2; k < m_nv; k++) {
198 y[k] = std::max(y[k], 0.0);
199 }
200}
201
202}
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 base class WallBase.
Base class for exceptions thrown by Cantera classes.
double upperBound(size_t k) const override
Get the upper bound on the k-th component of the local state vector.
void resetBadValues(double *y) override
Reset physically or mathematically problematic values, such as negative species concentrations.
vector< size_t > steadyConstraints() const override
Get the indices of equations that are algebraic constraints when solving the steady-state problem.
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.
void getState(double *y) override
Get the the current state of the reactor.
double lowerBound(size_t k) const override
Get the lower bound on the k-th component of the local state vector.
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.
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:36
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:411
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:232
virtual void setMassFractions_NoNorm(const double *const y)
Set the mass fractions to the specified values without normalizing.
Definition Phase.cpp:355
virtual void setPressure(double p)
Set the internally stored pressure (Pa) at constant temperature and composition.
Definition Phase.h:617
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
const double * massFractions() const
Return a const pointer to the mass fraction array.
Definition Phase.h:443
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:588
void getMassFractions(double *const y) const
Get the species mass fractions.
Definition Phase.cpp:471
FlowDevice & outlet(size_t n=0)
Return a reference to the n-th outlet FlowDevice connected to this reactor.
double m_pressure
Current pressure in the reactor [Pa].
virtual size_t nSurfs() const
Return the number of surfaces in a 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].
double m_mass
Current mass of the reactor [kg].
size_t m_nsp
Number of homogeneous species in the mixture.
double m_enthalpy
Current specific enthalpy of the reactor [J/kg].
virtual void evalSurfaces(double *LHS, double *RHS, double *sdot)
Evaluate terms related to surface reactions.
Definition Reactor.cpp:295
virtual void updateSurfaceState(double *y)
Update the state of SurfPhase objects attached to this reactor.
Definition Reactor.cpp:179
Kinetics * m_kin
Pointer to the homogeneous Kinetics object that handles the reactions.
Definition Reactor.h:290
vector< double > m_wdot
Species net molar production rates.
Definition Reactor.h:301
virtual void evalWalls(double t)
Evaluate terms related to Walls.
Definition Reactor.cpp:283
size_t neq()
Number of equations (state variables) for this reactor.
Definition Reactor.h:97
double m_Qdot
net heat transfer into the reactor, through walls [W]
Definition Reactor.h:294
vector< double > m_sdot
Production rates of gas phase species on surfaces [kmol/s].
Definition Reactor.h:299
virtual void getSurfaceInitialConditions(double *y)
Get initial conditions for SurfPhase objects attached to this reactor.
Definition Reactor.cpp:85
void initialize(double t0=0.0) override
Initialize the reactor.
Definition Reactor.cpp:94
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:450
virtual void updateConnected(bool updatePressure)
Update the state information needed by connected reactors, flow devices, and reactor walls.
Definition Reactor.cpp:188
Base class for a phase with thermodynamic properties.
virtual void setState_HP(double h, double p, double tol=1e-9)
Set the internally stored specific enthalpy (J/kg) and pressure (Pa) of the phase.
double enthalpy_mass() const
Specific enthalpy. Units: J/kg.
double dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
Definition utilities.h:82
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
const double Tiny
Small number to compare differences of mole fractions against.
Definition ct_defs.h:173
const double BigNumber
largest number to compare to inf.
Definition ct_defs.h:160
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...