Cantera  3.1.0a1
Loading...
Searching...
No Matches
IdealGasMoleReactor.cpp
Go to the documentation of this file.
1//! @file IdealGasMoleReactor.cpp A constant volume zero-dimensional
2//! 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("IdealGasMoleReactor::setThermoMgr",
24 "Incompatible phase type provided");
25 }
27}
28
30{
31 if (m_thermo == 0) {
32 throw CanteraError("IdealGasMoleReactor::getState",
33 "Error: reactor is empty.");
34 }
35 m_thermo->restoreState(m_state);
36
37 // get mass for calculations
38 m_mass = m_thermo->density() * m_vol;
39
40 // set the first component to the temperature
41 y[0] = m_thermo->temperature();
42
43 // set the second component to the volume
44 y[1] = m_vol;
45
46 // get moles of species in remaining state
47 getMoles(y + m_sidx);
48 // set the remaining components to the surface species moles on
49 // the walls
51}
52
53size_t IdealGasMoleReactor::componentIndex(const string& nm) const
54{
55 size_t k = speciesIndex(nm);
56 if (k != npos) {
57 return k + m_sidx;
58 } else if (nm == "temperature") {
59 return 0;
60 } else if (nm == "volume") {
61 return 1;
62 } else {
63 return npos;
64 }
65}
66
68{
69 if (k == 0) {
70 return "temperature";
71 } else {
73 }
74}
75
77{
79 m_uk.resize(m_nsp, 0.0);
80}
81
83{
84 // the components of y are: [0] the temperature, [1] the volume, [2...K+1) are the
85 // moles of each species, and [K+1...] are the moles of surface
86 // species on each wall.
88 m_vol = y[1];
89 // set state
90 m_thermo->setMolesNoTruncate(y + m_sidx);
91 m_thermo->setState_TD(y[0], m_mass / m_vol);
92 updateConnected(true);
94}
95
96void IdealGasMoleReactor::eval(double time, double* LHS, double* RHS)
97{
98 double& mcvdTdt = RHS[0]; // m * c_v * dT/dt
99 double* dndt = RHS + m_sidx; // kmol per s
100
101 evalWalls(time);
102
103 m_thermo->restoreState(m_state);
104
105 m_thermo->getPartialMolarIntEnergies(&m_uk[0]);
106 const vector<double>& imw = m_thermo->inverseMolecularWeights();
107
108 if (m_chem) {
109 m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
110 }
111
112 // evaluate surfaces
113 evalSurfaces(LHS + m_nsp + m_sidx, RHS + m_nsp + m_sidx, m_sdot.data());
114
115 // external heat transfer and compression work
116 mcvdTdt += - m_pressure * m_vdot + m_Qdot;
117
118 for (size_t n = 0; n < m_nsp; n++) {
119 // heat release from gas phase and surface reactions
120 mcvdTdt -= m_wdot[n] * m_uk[n] * m_vol;
121 mcvdTdt -= m_sdot[n] * m_uk[n];
122 // production in gas phase and from surfaces
123 dndt[n] = (m_wdot[n] * m_vol + m_sdot[n]);
124 }
125
126 // add terms for outlets
127 for (auto outlet : m_outlet) {
128 for (size_t n = 0; n < m_nsp; n++) {
129 // flow of species into system and dilution by other species
130 dndt[n] -= outlet->outletSpeciesMassFlowRate(n) * imw[n];
131 }
132 double mdot = outlet->massFlowRate();
133 mcvdTdt -= mdot * m_pressure * m_vol / m_mass; // flow work
134 }
135
136 // add terms for inlets
137 for (auto inlet : m_inlet) {
138 double mdot = inlet->massFlowRate();
139 mcvdTdt += inlet->enthalpy_mass() * mdot;
140 for (size_t n = 0; n < m_nsp; n++) {
141 double mdot_spec = inlet->outletSpeciesMassFlowRate(n);
142 // flow of species into system and dilution by other species
143 dndt[n] += inlet->outletSpeciesMassFlowRate(n) * imw[n];
144 mcvdTdt -= m_uk[n] * imw[n] * mdot_spec;
145 }
146 }
147
148 RHS[1] = m_vdot;
149 if (m_energy) {
150 LHS[0] = m_mass * m_thermo->cv_mass();
151 } else {
152 RHS[0] = 0;
153 }
154}
155
156Eigen::SparseMatrix<double> IdealGasMoleReactor::jacobian()
157{
158 if (m_nv == 0) {
159 throw CanteraError("IdealGasMoleReactor::jacobian",
160 "Reactor must be initialized first.");
161 }
162 // clear former jacobian elements
163 m_jac_trips.clear();
164 // dnk_dnj represents d(dot(n_k)) / d (n_j) but is first assigned as
165 // d (dot(omega)) / d c_j, it is later transformed appropriately.
166 Eigen::SparseMatrix<double> dnk_dnj = m_kin->netProductionRates_ddCi();
167 // species size that accounts for surface species
168 size_t ssize = m_nv - m_sidx;
169 // map derivatives from the surface chemistry jacobian
170 // to the reactor jacobian
171 if (!m_surfaces.empty()) {
172 vector<Eigen::Triplet<double>> species_trips;
173 for (int k = 0; k < dnk_dnj.outerSize(); k++) {
174 for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
175 species_trips.emplace_back(static_cast<int>(it.row()),
176 static_cast<int>(it.col()), it.value());
177 }
178 }
179 addSurfaceJacobian(species_trips);
180 dnk_dnj.resize(ssize, ssize);
181 dnk_dnj.setFromTriplets(species_trips.begin(), species_trips.end());
182 }
183 // add species to species derivatives elements to the jacobian
184 // calculate ROP derivatives, excluding the terms -n_i / (V * N) dc_i/dn_j
185 // as it substantially reduces matrix sparsity
186 for (int k = 0; k < dnk_dnj.outerSize(); k++) {
187 for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
188 m_jac_trips.emplace_back(static_cast<int>(it.row() + m_sidx),
189 static_cast<int>(it.col() + m_sidx), it.value());
190 }
191 }
192 // Temperature Derivatives
193 if (m_energy) {
194 // getting perturbed state for finite difference
195 double deltaTemp = m_thermo->temperature()
196 * std::sqrt(std::numeric_limits<double>::epsilon());
197 // finite difference temperature derivatives
198 vector<double> lhsPerturbed(m_nv, 1.0), lhsCurrent(m_nv, 1.0);
199 vector<double> rhsPerturbed(m_nv, 0.0), rhsCurrent(m_nv, 0.0);
200 vector<double> yCurrent(m_nv);
201 getState(yCurrent.data());
202 vector<double> yPerturbed = yCurrent;
203 // perturb temperature
204 yPerturbed[0] += deltaTemp;
205 // getting perturbed state
206 updateState(yPerturbed.data());
207 double time = (m_net != nullptr) ? m_net->time() : 0.0;
208 eval(time, lhsPerturbed.data(), rhsPerturbed.data());
209 // reset and get original state
210 updateState(yCurrent.data());
211 eval(time, lhsCurrent.data(), rhsCurrent.data());
212 // d ydot_j/dT
213 for (size_t j = 0; j < m_nv; j++) {
214 double ydotPerturbed = rhsPerturbed[j] / lhsPerturbed[j];
215 double ydotCurrent = rhsCurrent[j] / lhsCurrent[j];
216 m_jac_trips.emplace_back(static_cast<int>(j), 0,
217 (ydotPerturbed - ydotCurrent) / deltaTemp);
218 }
219 // d T_dot/dnj
220 Eigen::VectorXd netProductionRates = Eigen::VectorXd::Zero(ssize);
221 Eigen::VectorXd internal_energy = Eigen::VectorXd::Zero(ssize);
222 Eigen::VectorXd specificHeat = Eigen::VectorXd::Zero(ssize);
223 // getting species data
224 m_thermo->getPartialMolarIntEnergies(internal_energy.data());
225 m_kin->getNetProductionRates(netProductionRates.data());
226 m_thermo->getPartialMolarCp(specificHeat.data());
227 // convert Cp to Cv for ideal gas as Cp - Cv = R
228 for (size_t i = 0; i < m_nsp; i++) {
229 specificHeat[i] -= GasConstant;
230 netProductionRates[i] *= m_vol;
231 }
232 // scale net production rates by volume to get molar rate
233 double qdot = internal_energy.dot(netProductionRates);
234 // find the sum of n_i and cp_i
235 double NCv = 0.0;
236 double* moles = yCurrent.data() + m_sidx;
237 for (size_t i = 0; i < ssize; i++) {
238 NCv += moles[i] * specificHeat[i];
239 }
240 // make denominator beforehand
241 double denom = 1 / (NCv * NCv);
242 Eigen::VectorXd uk_dnkdnj_sums = dnk_dnj.transpose() * internal_energy;
243 // add derivatives to jacobian
244 for (size_t j = 0; j < ssize; j++) {
245 m_jac_trips.emplace_back(0, static_cast<int>(j + m_sidx),
246 (specificHeat[j] * qdot - NCv * uk_dnkdnj_sums[j]) * denom);
247 }
248 }
249 // convert triplets to sparse matrix
250 Eigen::SparseMatrix<double> jac(m_nv, m_nv);
251 jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
252 return jac;
253}
254
255}
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 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.
vector< double > m_uk
Species molar internal energies.
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.
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.
const size_t m_sidx
const value for the species start index
Definition MoleReactor.h:62
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...
string componentName(size_t k) override
Return the name of the solution component with index i.
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.
void restoreState(const vector< double > &state)
Restore a state saved on a previous call to saveState.
Definition Phase.cpp:260
void setState_TD(double t, double rho)
Set the internally stored temperature (K) and density (kg/m^3)
Definition Phase.cpp:377
double temperature() const
Temperature (K).
Definition Phase.h:562
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
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
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
double m_vdot
net rate of volume change from moving walls [m^3/s]
Definition Reactor.h:285
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 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 void getPartialMolarIntEnergies(double *ubar) const
Return an array of partial molar internal energies for the species in the mixture.
double cv_mass() const
Specific heat at constant volume. 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
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition ct_defs.h:120
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...