Cantera  3.2.0a4
Loading...
Searching...
No Matches
IdealGasConstPressureMoleReactor.cpp
Go to the documentation of this file.
1//! @file IdealGasConstPressureMoleReactor.cpp A constant pressure
2//! zero-dimensional reactor with moles as the 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
15#include <limits>
16
17namespace Cantera
18{
19
21{
22 if (m_thermo == 0) {
23 throw CanteraError("IdealGasConstPressureMoleReactor::getState",
24 "Error: reactor is empty.");
25 }
26 m_thermo->restoreState(m_state);
27 // get mass for calculations
28 m_mass = m_thermo->density() * m_vol;
29 // set the first component to the temperature
30 y[0] = m_thermo->temperature();
31 // get moles of species in remaining state
32 getMoles(y + m_sidx);
33 // set the remaining components to the surface species moles on the walls
35}
36
38{
39 if (m_thermo->type() != "ideal-gas") {
40 throw CanteraError("IdealGasConstPressureMoleReactor::initialize",
41 "Incompatible phase type '{}' provided", m_thermo->type());
42 }
44 m_hk.resize(m_nsp, 0.0);
45}
46
48{
49 // the components of y are: [0] the temperature, [1...K+1) are the
50 // moles of each species, and [K+1...] are the moles of surface
51 // species on each wall.
52 setMassFromMoles(y + m_sidx);
53 m_thermo->setMolesNoTruncate(y + m_sidx);
54 m_thermo->setState_TP(y[0], m_pressure);
55 m_vol = m_mass / m_thermo->density();
56 updateConnected(false);
57 updateSurfaceState(y + m_nsp + m_sidx);
58}
59
60void IdealGasConstPressureMoleReactor::eval(double time, double* LHS, double* RHS)
61{
62 double& mcpdTdt = RHS[0]; // m * c_p * dT/dt
63 double* dndt = RHS + m_sidx; // kmol per s
64
65 evalWalls(time);
66
67 m_thermo->restoreState(m_state);
68
69 m_thermo->getPartialMolarEnthalpies(&m_hk[0]);
70 const vector<double>& imw = m_thermo->inverseMolecularWeights();
71
72 if (m_chem) {
73 m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
74 }
75
76 // evaluate reactor surfaces
77 evalSurfaces(LHS + m_nsp + m_sidx, RHS + m_nsp + m_sidx, m_sdot.data());
78
79 // external heat transfer
80 mcpdTdt += m_Qdot;
81
82 for (size_t n = 0; n < m_nsp; n++) {
83 // heat release from gas phase and surface reactions
84 mcpdTdt -= m_wdot[n] * m_hk[n] * m_vol;
85 mcpdTdt -= m_sdot[n] * m_hk[n];
86 // production in gas phase and from surfaces
87 dndt[n] = m_wdot[n] * m_vol + m_sdot[n];
88 }
89
90 // add terms for outlets
91 for (auto outlet : m_outlet) {
92 for (size_t n = 0; n < m_nsp; n++) {
93 // flow of species into system and dilution by other species
94 dndt[n] -= outlet->outletSpeciesMassFlowRate(n) * imw[n];
95 }
96 }
97
98 // add terms for inlets
99 for (auto inlet : m_inlet) {
100 double mdot = inlet->massFlowRate();
101 mcpdTdt += inlet->enthalpy_mass() * mdot;
102 for (size_t n = 0; n < m_nsp; n++) {
103 double mdot_spec = inlet->outletSpeciesMassFlowRate(n);
104 // flow of species into system and dilution by other species
105 dndt[n] += inlet->outletSpeciesMassFlowRate(n) * imw[n];
106 mcpdTdt -= m_hk[n] * imw[n] * mdot_spec;
107 }
108 }
109
110 if (m_energy) {
111 LHS[0] = m_mass * m_thermo->cp_mass();
112 } else {
113 RHS[0] = 0.0;
114 }
115}
116
118{
119 if (m_nv == 0) {
120 throw CanteraError("IdealGasConstPressureMoleReactor::jacobian",
121 "Reactor must be initialized first.");
122 }
123 // clear former jacobian elements
124 m_jac_trips.clear();
125 // dnk_dnj represents d(dot(n_k)) / d (n_j) but is first assigned as
126 // d (dot(omega)) / d c_j, it is later transformed appropriately.
127 Eigen::SparseMatrix<double> dnk_dnj = m_kin->netProductionRates_ddCi();
128 // species size that accounts for surface species
129 size_t ssize = m_nv - m_sidx;
130 // map derivatives from the surface chemistry jacobian
131 // to the reactor jacobian
132 if (!m_surfaces.empty()) {
133 vector<Eigen::Triplet<double>> species_trips(dnk_dnj.nonZeros());
134 for (int k = 0; k < dnk_dnj.outerSize(); k++) {
135 for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
136 species_trips.emplace_back(static_cast<int>(it.row()),
137 static_cast<int>(it.col()), it.value());
138 }
139 }
140 addSurfaceJacobian(species_trips);
141 dnk_dnj.resize(ssize, ssize);
142 dnk_dnj.setFromTriplets(species_trips.begin(), species_trips.end());
143 }
144 // get net production rates
145 Eigen::VectorXd netProductionRates = Eigen::VectorXd::Zero(ssize);
146 // gas phase net production rates
147 m_kin->getNetProductionRates(netProductionRates.data());
148 // surface phase net production rates mapped to reactor gas phase
149 for (auto &S: m_surfaces) {
150 auto curr_kin = S->kinetics();
151 vector<double> prod_rates(curr_kin->nTotalSpecies());
152 curr_kin->getNetProductionRates(prod_rates.data());
153 for (size_t i = 0; i < curr_kin->nTotalSpecies(); i++) {
154 size_t row = speciesIndex(curr_kin->kineticsSpeciesName(i));
155 if (row != npos) {
156 netProductionRates[row] += prod_rates[i];
157 }
158 }
159 }
160 double molarVol = m_thermo->molarVolume();
161 // add species to species derivatives elements to the jacobian
162 // calculate ROP derivatives, excluding the terms -n_i / (V * N) dc_i/dn_j
163 // as it substantially reduces matrix sparsity
164 for (int k = 0; k < dnk_dnj.outerSize(); k++) {
165 for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
166 // gas phase species need the addition of V / N * omega_dot
167 if (static_cast<size_t>(it.row()) < m_nsp) {
168 it.valueRef() = it.value() + netProductionRates[it.row()] * molarVol;
169 }
170 m_jac_trips.emplace_back(static_cast<int>(it.row() + m_sidx),
171 static_cast<int>(it.col() + m_sidx), it.value());
172 }
173 }
174 // Temperature Derivatives
175 if (m_energy) {
176 // getting perturbed state for finite difference
177 double deltaTemp = m_thermo->temperature()
178 * std::sqrt(std::numeric_limits<double>::epsilon());
179 // get current state
180 vector<double> yCurrent(m_nv);
181 getState(yCurrent.data());
182 // finite difference temperature derivatives
183 vector<double> lhsPerturbed(m_nv, 1.0), lhsCurrent(m_nv, 1.0);
184 vector<double> rhsPerturbed(m_nv, 0.0), rhsCurrent(m_nv, 0.0);
185 vector<double> yPerturbed = yCurrent;
186 // perturb temperature
187 yPerturbed[0] += deltaTemp;
188 // getting perturbed state
189 updateState(yPerturbed.data());
190 double time = (m_net != nullptr) ? m_net->time() : 0.0;
191 eval(time, lhsPerturbed.data(), rhsPerturbed.data());
192 // reset and get original state
193 updateState(yCurrent.data());
194 eval(time, lhsCurrent.data(), rhsCurrent.data());
195 // d ydot_j/dT
196 for (size_t j = 0; j < m_nv; j++) {
197 double ydotPerturbed = rhsPerturbed[j] / lhsPerturbed[j];
198 double ydotCurrent = rhsCurrent[j] / lhsCurrent[j];
199 m_jac_trips.emplace_back(static_cast<int>(j), 0,
200 (ydotPerturbed - ydotCurrent) / deltaTemp);
201 }
202 // d T_dot/dnj
203 // allocate vectors for whole system
204 Eigen::VectorXd enthalpy = Eigen::VectorXd::Zero(ssize);
205 Eigen::VectorXd specificHeat = Eigen::VectorXd::Zero(ssize);
206 // gas phase
207 m_thermo->getPartialMolarCp(specificHeat.data());
208 m_thermo->getPartialMolarEnthalpies(enthalpy.data());
209 // scale production rates by the volume for gas species
210 for (size_t i = 0; i < m_nsp; i++) {
211 netProductionRates[i] *= m_vol;
212 }
213 double qdot = enthalpy.dot(netProductionRates);
214 // find denominator ahead of time
215 double NCp = 0.0;
216 double* moles = yCurrent.data() + m_sidx;
217 for (size_t i = 0; i < ssize; i++) {
218 NCp += moles[i] * specificHeat[i];
219 }
220 double denom = 1 / (NCp * NCp);
221 Eigen::VectorXd hk_dnkdnj_sums = dnk_dnj.transpose() * enthalpy;
222 // Add derivatives to jac by spanning columns
223 for (size_t j = 0; j < ssize; j++) {
224 m_jac_trips.emplace_back(0, static_cast<int>(j + m_sidx),
225 (specificHeat[j] * qdot - NCp * hk_dnkdnj_sums[j]) * denom);
226 }
227 }
228 // convert triplets to sparse matrix
229 Eigen::SparseMatrix<double> jac(m_nv, m_nv);
230 jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
231 return jac;
232}
233
235{
236 size_t k = speciesIndex(nm);
237 if (k != npos) {
238 return k + m_sidx;
239 } else if (nm == "temperature") {
240 return 0;
241 } else {
242 return npos;
243 }
244}
245
247 if (k == 0) {
248 return "temperature";
249 } else if (k >= m_sidx && k < neq()) {
250 k -= m_sidx;
251 if (k < m_thermo->nSpecies()) {
252 return m_thermo->speciesName(k);
253 } else {
254 k -= m_thermo->nSpecies();
255 }
256 for (auto& S : m_surfaces) {
257 ThermoPhase* th = S->thermo();
258 if (k < th->nSpecies()) {
259 return th->speciesName(k);
260 } else {
261 k -= th->nSpecies();
262 }
263 }
264 }
265 throw CanteraError("IdealGasConstPressureMoleReactor::componentName",
266 "Index is out of bounds.");
267}
268
270 if (k == 0) {
271 //@todo: Revise pending resolution of https://github.com/Cantera/enhancements/issues/229
272 return 1.5 * m_thermo->maxTemp();
273 } else {
274 return BigNumber; // moles of a bulk or surface species
275 }
276}
277
279 if (k == 0) {
280 //@todo: Revise pending resolution of https://github.com/Cantera/enhancements/issues/229
281 return 0.5 * m_thermo->minTemp();
282 } else {
284 }
285}
286
287}
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.
double lowerBound(size_t k) const override
Get the lower bound on the k-th component of the local state vector.
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
double upperBound(size_t k) const override
Get the upper bound on the k-th component of the local state vector.
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.
Eigen::SparseMatrix< double > jacobian() override
Calculate an approximate Jacobian to accelerate preconditioned solvers.
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.
vector< double > m_hk
Species molar enthalpies.
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:411
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 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 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:286
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:232
double temperature() const
Temperature (K).
Definition Phase.h:587
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
const vector< double > & inverseMolecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition Phase.cpp:428
virtual double density() const
Density (kg/m^3).
Definition Phase.h:612
virtual void setMolesNoTruncate(const double *const N)
Set the state of the object with moles in [kmol].
Definition Phase.cpp:556
virtual double molarVolume() const
Molar volume (m^3/kmol).
Definition Phase.cpp:609
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].
ReactorNet * m_net
The ReactorNet that this reactor is part of.
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 time()
Current value of the simulation time [s], for reactor networks that are solved in the time domain.
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< Eigen::Triplet< double > > m_jac_trips
Vector of triplets representing the jacobian.
Definition Reactor.h:311
vector< double > m_sdot
Production rates of gas phase species on surfaces [kmol/s].
Definition Reactor.h:299
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 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)
virtual double minTemp(size_t k=npos) const
Minimum temperature for which the thermodynamic data for the species or phase are valid.
virtual void getPartialMolarCp(double *cpbar) const
Return an array of partial molar heat capacities for the species in the mixture.
string type() const override
String indicating the thermodynamic model implemented.
virtual double maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
double cp_mass() const
Specific heat at constant pressure. Units: J/kg/K.
Eigen::SparseMatrix< double > netProductionRates_ddCi()
Calculate derivatives for species net production rates with respect to species concentration at const...
Definition Kinetics.cpp:563
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 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...