Cantera  3.1.0a1
Loading...
Searching...
No Matches
MultiPhaseEquil.h
Go to the documentation of this file.
1//! @file MultiPhaseEquil.h
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
6#ifndef CT_MULTIPHASE_EQUIL
7#define CT_MULTIPHASE_EQUIL
8
9#include "MultiPhase.h"
10
11namespace Cantera
12{
13
14/**
15 * Multiphase chemical equilibrium solver. Class MultiPhaseEquil is designed
16 * to be used to set a mixture containing one or more phases to a state of
17 * chemical equilibrium. It implements the VCS algorithm, described in Smith
18 * and Missen @cite smith1982.
19 *
20 * This class only handles chemical equilibrium at a specified temperature and
21 * pressure. To compute equilibrium holding other properties fixed, it is
22 * necessary to iterate on T and P in an "outer" loop, until the specified
23 * properties have the desired values. This is done, for example, in method
24 * equilibrate of class MultiPhase.
25 *
26 * This class is primarily meant to be used internally by the equilibrate
27 * method of class MultiPhase, although there is no reason it cannot be used
28 * directly in application programs if desired.
29 *
30 * @ingroup equilGroup
31 */
33{
34public:
35 //! Construct a multiphase equilibrium manager for a multiphase mixture.
36 //! @param mix Pointer to a multiphase mixture object.
37 //! @param start If true, the initial composition will be determined by a
38 //! linear Gibbs minimization, otherwise the initial mixture
39 //! composition will be used.
40 //! @param loglevel Desired level of debug printing. loglevel = 0 suppresses
41 //! printing. Higher values request more verbose logging output.
42 MultiPhaseEquil(MultiPhase* mix, bool start=true, int loglevel = 0);
43
44 virtual ~MultiPhaseEquil() {}
45
46 size_t constituent(size_t m) {
47 if (m < m_nel) {
48 return m_order[m];
49 } else {
50 return npos;
51 }
52 }
53
54 void getStoichVector(size_t rxn, vector<double>& nu) {
55 nu.resize(m_nsp, 0.0);
56 if (rxn > nFree()) {
57 return;
58 }
59 for (size_t k = 0; k < m_nsp; k++) {
60 nu[m_order[k]] = m_N(k, rxn);
61 }
62 }
63
64 int iterations() {
65 return m_iter;
66 }
67
68 double equilibrate(int XY, double err = 1.0e-9,
69 int maxsteps = 1000, int loglevel=-99);
70 double error();
71
72 string reactionString(size_t j) {
73 return "";
74 }
75 void setInitialMixMoles(int loglevel = 0) {
76 setInitialMoles(loglevel);
77 finish();
78 }
79
80 size_t componentIndex(size_t n) {
81 return m_species[m_order[n]];
82 }
83
84 void reportCSV(const string& reportFile);
85
86 double phaseMoles(size_t iph) const;
87
88protected:
89 //! This method finds a set of component species and a complete set of
90 //! formation reactions for the non-components in terms of the components.
91 //! In most cases, many different component sets are possible, and
92 //! therefore neither the components returned by this method nor the
93 //! formation reactions are unique. The algorithm used here is described
94 //! in Smith and Missen @cite smith1982.
95 //!
96 //! The component species are taken to be the first M species in array
97 //! 'species' that have linearly-independent compositions.
98 //!
99 //! @param order On entry, vector @e order should contain species index
100 //! numbers in the order of decreasing desirability as a component.
101 //! For example, if it is desired to choose the components from among
102 //! the major species, this array might list species index numbers in
103 //! decreasing order of mole fraction. If array 'species' does not
104 //! have length = nSpecies(), then the species will be considered as
105 //! candidates to be components in declaration order, beginning with
106 //! the first phase added.
107 void getComponents(const vector<size_t>& order);
108
109 //! Estimate the initial mole numbers. This is done by running each
110 //! reaction as far forward or backward as possible, subject to the
111 //! constraint that all mole numbers remain non-negative. Reactions for
112 //! which @f$ \Delta \mu^0 @f$ are positive are run in reverse, and ones
113 //! for which it is negative are run in the forward direction. The end
114 //! result is equivalent to solving the linear programming problem of
115 //! minimizing the linear Gibbs function subject to the element and non-
116 //! negativity constraints.
117 int setInitialMoles(int loglevel = 0);
118
119 void computeN();
120
121 //! Take one step in composition, given the gradient of G at the starting
122 //! point, and a vector of reaction steps dxi.
123 double stepComposition(int loglevel = 0);
124
125 //! Re-arrange a vector of species properties in sorted form
126 //! (components first) into unsorted, sequential form.
127 void unsort(vector<double>& x);
128
129 void step(double omega, vector<double>& deltaN, int loglevel = 0);
130
131 //! Compute the change in extent of reaction for each reaction.
132 double computeReactionSteps(vector<double>& dxi);
133
134 void updateMixMoles();
135
136 //! Clean up the composition. The solution algorithm can leave some
137 //! species in stoichiometric condensed phases with very small negative
138 //! mole numbers. This method simply sets these to zero.
139 void finish();
140
141 // moles of the species with sorted index ns
142 double moles(size_t ns) const {
143 return m_moles[m_order[ns]];
144 }
145 double& moles(size_t ns) {
146 return m_moles[m_order[ns]];
147 }
148 int solutionSpecies(size_t n) const {
149 return m_dsoln[m_order[n]];
150 }
151 bool isStoichPhase(size_t n) const {
152 return (m_dsoln[m_order[n]] == 0);
153 }
154 double mu(size_t n) const {
155 return m_mu[m_species[m_order[n]]];
156 }
157 string speciesName(size_t n) const {
158 return
159 m_mix->speciesName(m_species[m_order[n]]);
160 }
161
162 //! Number of degrees of freedom
163 size_t nFree() const {
164 return (m_nsp > m_nel) ? m_nsp - m_nel : 0;
165 }
166
167 size_t m_nel_mix, m_nsp_mix;
168 size_t m_nel = 0;
169 size_t m_nsp = 0;
170 size_t m_eloc = 1000;
171 int m_iter;
172 MultiPhase* m_mix;
173 double m_press, m_temp;
174 vector<size_t> m_order;
175 DenseMatrix m_N, m_A;
176 vector<double> m_work, m_work2, m_work3;
177 vector<double> m_moles, m_lastmoles, m_dxi;
178 vector<double> m_deltaG_RT, m_mu;
179 vector<bool> m_majorsp;
180 vector<size_t> m_sortindex;
181 vector<int> m_lastsort;
182 vector<int> m_dsoln;
183 vector<int> m_incl_element, m_incl_species;
184
185 // Vector of indices for species that are included in the calculation.
186 // This is used to exclude pure-phase species with invalid thermo data
187 vector<size_t> m_species;
188 vector<size_t> m_element;
189 vector<bool> m_solnrxn;
190 bool m_force = true;
191};
192
193}
194
195#endif
Headers for the MultiPhase object that is used to set up multiphase equilibrium problems (see Chemica...
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition DenseMatrix.h:55
Multiphase chemical equilibrium solver.
double computeReactionSteps(vector< double > &dxi)
Compute the change in extent of reaction for each reaction.
double stepComposition(int loglevel=0)
Take one step in composition, given the gradient of G at the starting point, and a vector of reaction...
void finish()
Clean up the composition.
void unsort(vector< double > &x)
Re-arrange a vector of species properties in sorted form (components first) into unsorted,...
size_t nFree() const
Number of degrees of freedom.
int setInitialMoles(int loglevel=0)
Estimate the initial mole numbers.
void getComponents(const vector< size_t > &order)
This method finds a set of component species and a complete set of formation reactions for the non-co...
A class for multiphase mixtures.
Definition MultiPhase.h:57
string speciesName(const size_t kGlob) const
Name of species with global index kGlob.
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