forked from Cantera/cantera
-
Notifications
You must be signed in to change notification settings - Fork 0
/
PSRv0-3.cpp
110 lines (94 loc) · 3.75 KB
/
PSRv0-3.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
/// @file PSRv0-3.cpp
/// Perfectly Stirred Reactor Solver (v0.3)
/// Changes since v0.2.1:
/// - Created a separate module for nonlinear system solution (Cantera1D_NonLinSol)
/// - Removed nonlinear solver capability from PSR
#include "cantera/numerics/Cantera_NonLinSol.h" //a nonlinear algebraic system solver
const char *MECHANISM = "gri30.yaml";
const double PRESSURE = Cantera::OneAtm;
class PSR : public Cantera_NonLinSol //inherit nonlinear solver functionality
{
public:
// Constructor
PSR(std::shared_ptr<Cantera::Solution> inletSol, std::shared_ptr<Cantera::Solution> reactorSol, double mdot, double volume)
// ------------------------ INITIALIZE GLOBAL VARS ------------------------
: inletThermo(inletSol->thermo()),
reactorThermo(reactorSol->thermo()),
reactorKinetics(reactorSol->kinetics()),
nSpecies(inletSol->thermo()->nSpecies()),
Y_in(inletSol->thermo()->massFractions()), // mass fractions at the inlet
molecularWeight(inletSol->thermo()->molecularWeights()), // molecular weight of each species (these never change)
volume(volume),
mdot(mdot)
{
}
/// Providing implementations for the solver's virtual functions:
// Specify guesses for the initial values.
// Note: called during Sim1D initialization
doublereal ctNLS_initialValue(size_t i)
{
return reactorThermo->massFraction(i); // mass fractions of the provided intial reactorSol (aka the initial guess)
}
// Number of equations (state variables) for this reactor
size_t ctNLS_nEqs()
{
return nSpecies;
}
// Advance the PSR to steady state
void solveSteady()
{
Cantera_NonLinSol::solve();
}
// Specify the residual function for the system
// sol - iteration solution vector (input)
// rsd - residual vector (output)
void ctNLS_residFunction(double *sol, double *rsd)
{
// ----------------------- UPDATE REACTOR STATE -----------------------
reactorThermo->setMassFractions_NoNorm(sol);
reactorThermo->setState_HP(inletThermo->enthalpy_mass(), PRESSURE); // keep total enthalpy constant, allow Cantera to control reactor temperature
// ----------------------- GET REQ'D PROPERTIES -----------------------
doublereal wdot[nSpecies];
reactorKinetics->getNetProductionRates(wdot);
// ----------------------- SPECIES CONSERVATION -----------------------
for (int k = 0; k < nSpecies; k++)
{
rsd[k] = wdot[k] * molecularWeight[k] * volume + mdot * (Y_in[k] - sol[k]); //PSR residual equation
}
}
private:
std::shared_ptr<Cantera::ThermoPhase> inletThermo;
std::shared_ptr<Cantera::ThermoPhase> reactorThermo;
std::shared_ptr<Cantera::Kinetics> reactorKinetics;
int nSpecies;
const double *Y_in;
Cantera::vector_fp molecularWeight;
double volume;
double mdot;
};
int main()
{
try
{
std::shared_ptr<Cantera::Solution> mixInlet(Cantera::newSolution(MECHANISM));
std::shared_ptr<Cantera::Solution> mixOutlet(Cantera::newSolution(MECHANISM));
mixInlet->thermo()->setState_TPX(300, PRESSURE, "O2:1 CH4:0.5 N2:3.76");
mixOutlet->thermo()->setState_TPX(300, PRESSURE, "O2:1 CH4:0.5 N2:3.76");
mixOutlet->thermo()->equilibrate("HP");
double mdot;
double volume;
std::cout << "Enter mdot: ";
std::cin >> mdot;
std::cout << "Enter volume: ";
std::cin >> volume;
PSR reactor(mixInlet, mixOutlet, mdot, volume);
reactor.solveSteady();
std::cout << mixOutlet->thermo()->report() << "\n";
return 0;
}
catch (Cantera::CanteraError &err)
{
std::cerr << err.what() << std::endl;
return -1;
}
}