Cantera  3.1.0a1
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 (thermo.type() != "ideal-gas") {
23 throw CanteraError("IdealGasConstPressureMoleReactor::setThermoMgr",
24 "Incompatible phase type provided");
25 }
27}
28
30{
31 if (m_thermo == 0) {
32 throw CanteraError("IdealGasConstPressureMoleReactor::getState",
33 "Error: reactor is empty.");
34 }
35 m_thermo->restoreState(m_state);
36 // get mass for calculations
37 m_mass = m_thermo->density() * m_vol;
38 // set the first component to the temperature
39 y[0] = m_thermo->temperature();
40 // get moles of species in remaining state
41 getMoles(y + m_sidx);
42 // set the remaining components to the surface species moles on the walls
44}
45
47{
49 m_hk.resize(m_nsp, 0.0);
50}
51
53{
54 // the components of y are: [0] the temperature, [1...K+1) are the
55 // moles of each species, and [K+1...] are the moles of surface
56 // species on each wall.
57 setMassFromMoles(y + m_sidx);
58 m_thermo->setMolesNoTruncate(y + m_sidx);
59 m_thermo->setState_TP(y[0], m_pressure);
60 m_vol = m_mass / m_thermo->density();
61 updateConnected(false);
62 updateSurfaceState(y + m_nsp + m_sidx);
63}
64
65void IdealGasConstPressureMoleReactor::eval(double time, double* LHS, double* RHS)
66{
67 double& mcpdTdt = RHS[0]; // m * c_p * dT/dt
68 double* dndt = RHS + m_sidx; // kmol per s
69
70 evalWalls(time);
71
72 m_thermo->restoreState(m_state);
73
74 m_thermo->getPartialMolarEnthalpies(&m_hk[0]);
75 const vector<double>& imw = m_thermo->inverseMolecularWeights();
76
77 if (m_chem) {
78 m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
79 }
80
81 // evaluate reactor surfaces
82 evalSurfaces(LHS + m_nsp + m_sidx, RHS + m_nsp + m_sidx, m_sdot.data());
83
84 // external heat transfer
85 mcpdTdt += m_Qdot;
86
87 for (size_t n = 0; n < m_nsp; n++) {
88 // heat release from gas phase and surface reactions
89 mcpdTdt -= m_wdot[n] * m_hk[n] * m_vol;
90 mcpdTdt -= m_sdot[n] * m_hk[n];
91 // production in gas phase and from surfaces
92 dndt[n] = m_wdot[n] * m_vol + m_sdot[n];
93 }
94
95 // add terms for outlets
96 for (auto outlet : m_outlet) {
97 for (size_t n = 0; n < m_nsp; n++) {
98 // flow of species into system and dilution by other species
99 dndt[n] -= outlet->outletSpeciesMassFlowRate(n) * imw[n];
100 }
101 }
102
103 // add terms for inlets
104 for (auto inlet : m_inlet) {
105 double mdot = inlet->massFlowRate();
106 mcpdTdt += inlet->enthalpy_mass() * mdot;
107 for (size_t n = 0; n < m_nsp; n++) {
108 double mdot_spec = inlet->outletSpeciesMassFlowRate(n);
109 // flow of species into system and dilution by other species
110 dndt[n] += inlet->outletSpeciesMassFlowRate(n) * imw[n];
111 mcpdTdt -= m_hk[n] * imw[n] * mdot_spec;
112 }
113 }
114
115 if (m_energy) {
116 LHS[0] = m_mass * m_thermo->cp_mass();
117 } else {
118 RHS[0] = 0.0;
119 }
120}
121
123{
124 if (m_nv == 0) {
125 throw CanteraError("IdealGasConstPressureMoleReactor::jacobian",
126 "Reactor must be initialized first.");
127 }
128 // clear former jacobian elements
129 m_jac_trips.clear();
130 // dnk_dnj represents d(dot(n_k)) / d (n_j) but is first assigned as
131 // d (dot(omega)) / d c_j, it is later transformed appropriately.
132 Eigen::SparseMatrix<double> dnk_dnj = m_kin->netProductionRates_ddCi();
133 // species size that accounts for surface species
134 size_t ssize = m_nv - m_sidx;
135 // map derivatives from the surface chemistry jacobian
136 // to the reactor jacobian
137 if (!m_surfaces.empty()) {
138 vector<Eigen::Triplet<double>> species_trips(dnk_dnj.nonZeros());
139 for (int k = 0; k < dnk_dnj.outerSize(); k++) {
140 for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
141 species_trips.emplace_back(static_cast<int>(it.row()),
142 static_cast<int>(it.col()), it.value());
143 }
144 }
145 addSurfaceJacobian(species_trips);
146 dnk_dnj.resize(ssize, ssize);
147 dnk_dnj.setFromTriplets(species_trips.begin(), species_trips.end());
148 }
149 // get net production rates
150 Eigen::VectorXd netProductionRates = Eigen::VectorXd::Zero(ssize);
151 // gas phase net production rates
152 m_kin->getNetProductionRates(netProductionRates.data());
153 // surface phase net production rates mapped to reactor gas phase
154 for (auto &S: m_surfaces) {
155 auto curr_kin = S->kinetics();
156 vector<double> prod_rates(curr_kin->nTotalSpecies());
157 curr_kin->getNetProductionRates(prod_rates.data());
158 for (size_t i = 0; i < curr_kin->nTotalSpecies(); i++) {
159 size_t row = speciesIndex(curr_kin->kineticsSpeciesName(i));
160 if (row != npos) {
161 netProductionRates[row] += prod_rates[i];
162 }
163 }
164 }
165 double molarVol = m_thermo->molarVolume();
166 // add species to species derivatives elements to the jacobian
167 // calculate ROP derivatives, excluding the terms -n_i / (V * N) dc_i/dn_j
168 // as it substantially reduces matrix sparsity
169 for (int k = 0; k < dnk_dnj.outerSize(); k++) {
170 for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
171 // gas phase species need the addition of V / N * omega_dot
172 if (static_cast<size_t>(it.row()) < m_nsp) {
173 it.valueRef() = it.value() + netProductionRates[it.row()] * molarVol;
174 }
175 m_jac_trips.emplace_back(static_cast<int>(it.row() + m_sidx),
176 static_cast<int>(it.col() + m_sidx), it.value());
177 }
178 }
179 // Temperature Derivatives
180 if (m_energy) {
181 // getting perturbed state for finite difference
182 double deltaTemp = m_thermo->temperature()
183 * std::sqrt(std::numeric_limits<double>::epsilon());
184 // get current state
185 vector<double> yCurrent(m_nv);
186 getState(yCurrent.data());
187 // finite difference temperature derivatives
188 vector<double> lhsPerturbed(m_nv, 1.0), lhsCurrent(m_nv, 1.0);
189 vector<double> rhsPerturbed(m_nv, 0.0), rhsCurrent(m_nv, 0.0);
190 vector<double> yPerturbed = yCurrent;
191 // perturb temperature
192 yPerturbed[0] += deltaTemp;
193 // getting perturbed state
194 updateState(yPerturbed.data());
195 double time = (m_net != nullptr) ? m_net->time() : 0.0;
196 eval(time, lhsPerturbed.data(), rhsPerturbed.data());
197 // reset and get original state
198 updateState(yCurrent.data());
199 eval(time, lhsCurrent.data(), rhsCurrent.data());
200 // d ydot_j/dT
201 for (size_t j = 0; j < m_nv; j++) {
202 double ydotPerturbed = rhsPerturbed[j] / lhsPerturbed[j];
203 double ydotCurrent = rhsCurrent[j] / lhsCurrent[j];
204 m_jac_trips.emplace_back(static_cast<int>(j), 0,
205 (ydotPerturbed - ydotCurrent) / deltaTemp);
206 }
207 // d T_dot/dnj
208 // allocate vectors for whole system
209 Eigen::VectorXd enthalpy = Eigen::VectorXd::Zero(ssize);
210 Eigen::VectorXd specificHeat = Eigen::VectorXd::Zero(ssize);
211 // gas phase
212 m_thermo->getPartialMolarCp(specificHeat.data());
213 m_thermo->getPartialMolarEnthalpies(enthalpy.data());
214 // scale production rates by the volume for gas species
215 for (size_t i = 0; i < m_nsp; i++) {
216 netProductionRates[i] *= m_vol;
217 }
218 double qdot = enthalpy.dot(netProductionRates);
219 // find denominator ahead of time
220 double NCp = 0.0;
221 double* moles = yCurrent.data() + m_sidx;
222 for (size_t i = 0; i < ssize; i++) {
223 NCp += moles[i] * specificHeat[i];
224 }
225 double denom = 1 / (NCp * NCp);
226 Eigen::VectorXd hk_dnkdnj_sums = dnk_dnj.transpose() * enthalpy;
227 // Add derivatives to jac by spanning columns
228 for (size_t j = 0; j < ssize; j++) {
229 m_jac_trips.emplace_back(0, static_cast<int>(j + m_sidx),
230 (specificHeat[j] * qdot - NCp * hk_dnkdnj_sums[j]) * denom);
231 }
232 }
233 // convert triplets to sparse matrix
234 Eigen::SparseMatrix<double> jac(m_nv, m_nv);
235 jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
236 return jac;
237}
238
240{
241 size_t k = speciesIndex(nm);
242 if (k != npos) {
243 return k + m_sidx;
244 } else if (nm == "temperature") {
245 return 0;
246 } else {
247 return npos;
248 }
249}
250
252 if (k == 0) {
253 return "temperature";
254 } else if (k >= m_sidx && k < neq()) {
255 k -= m_sidx;
256 if (k < m_thermo->nSpecies()) {
257 return m_thermo->speciesName(k);
258 } else {
259 k -= m_thermo->nSpecies();
260 }
261 for (auto& S : m_surfaces) {
262 ThermoPhase* th = S->thermo();
263 if (k < th->nSpecies()) {
264 return th->speciesName(k);
265 } else {
266 k -= th->nSpecies();
267 }
268 }
269 }
270 throw CanteraError("IdealGasConstPressureMoleReactor::componentName",
271 "Index is out of bounds.");
272}
273
274}
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.
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:39
void setThermoMgr(ThermoPhase &thermo) override
Specify the mixture contained in the reactor.
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.
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:363
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:260
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:231
double temperature() const
Temperature (K).
Definition Phase.h:562
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:400
virtual double density() const
Density (kg/m^3).
Definition Phase.h:587
virtual void setMolesNoTruncate(const double *const N)
Set the state of the object with moles in [kmol].
Definition Phase.cpp:528
virtual double molarVolume() const
Molar volume (m^3/kmol).
Definition Phase.cpp:581
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].
size_t m_nsp
Number of homogeneous species in the mixture.
virtual void setThermoMgr(ThermoPhase &thermo)
Specify the mixture contained in the reactor.
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:283
vector< double > m_wdot
Species net molar production rates.
Definition Reactor.h:295
virtual void evalWalls(double t)
Evaluate terms related to Walls.
Definition Reactor.cpp:278
size_t neq()
Number of equations (state variables) for this reactor.
Definition Reactor.h:109
double m_Qdot
net heat transfer into the reactor, through walls [W]
Definition Reactor.h:287
double m_mass
total mass
Definition Reactor.h:289
vector< Eigen::Triplet< double > > m_jac_trips
Vector of triplets representing the jacobian.
Definition Reactor.h:308
vector< double > m_sdot
Production rates of gas phase species on surfaces [kmol/s].
Definition Reactor.h:293
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
virtual void updateConnected(bool updatePressure)
Update the state information needed by connected reactors, flow devices, and reactor walls.
Definition Reactor.cpp:187
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 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.
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:515
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
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...