Cantera  3.1.0a1
Loading...
Searching...
No Matches
LatticeSolidPhase.cpp
Go to the documentation of this file.
1/**
2 * @file LatticeSolidPhase.cpp
3 * Definitions for a simple thermodynamics model of a bulk solid phase
4 * derived from ThermoPhase,
5 * assuming an ideal solution model based on a lattice of solid atoms
6 * (see @ref thermoprops and class @link Cantera::LatticeSolidPhase LatticeSolidPhase@endlink).
7 */
8
9// This file is part of Cantera. See License.txt in the top-level directory or
10// at https://cantera.org/license.txt for license and copyright information.
11
18
19#include <boost/algorithm/string.hpp>
20
21namespace ba = boost::algorithm;
22
23namespace Cantera
24{
25
26double LatticeSolidPhase::minTemp(size_t k) const
27{
28 if (k != npos) {
29 for (size_t n = 0; n < m_lattice.size(); n++) {
30 if (lkstart_[n+1] < k) {
31 return m_lattice[n]->minTemp(k-lkstart_[n]);
32 }
33 }
34 }
35 double mm = 0;
36 for (auto& lattice : m_lattice) {
37 mm = std::max(mm, lattice->minTemp());
38 }
39 return mm;
40}
41
42double LatticeSolidPhase::maxTemp(size_t k) const
43{
44 if (k != npos) {
45 for (size_t n = 0; n < m_lattice.size(); n++) {
46 if (lkstart_[n+1] < k) {
47 return (m_lattice[n])->maxTemp(k - lkstart_[n]);
48 }
49 }
50 }
51 double mm = BigNumber;
52 for (auto& lattice : m_lattice) {
53 mm = std::min(mm, lattice->maxTemp());
54 }
55 return mm;
56}
57
59{
60 return m_lattice[0]->refPressure();
61}
62
64{
66 double sum = 0.0;
67 for (size_t n = 0; n < m_lattice.size(); n++) {
68 sum += theta_[n] * m_lattice[n]->enthalpy_mole();
69 }
70 return sum;
71}
72
74{
76 double sum = 0.0;
77 for (size_t n = 0; n < m_lattice.size(); n++) {
78 sum += theta_[n] * m_lattice[n]->intEnergy_mole();
79 }
80 return sum;
81}
82
84{
86 double sum = 0.0;
87 for (size_t n = 0; n < m_lattice.size(); n++) {
88 sum += theta_[n] * m_lattice[n]->entropy_mole();
89 }
90 return sum;
91}
92
94{
96 double sum = 0.0;
97 for (size_t n = 0; n < m_lattice.size(); n++) {
98 sum += theta_[n] * m_lattice[n]->gibbs_mole();
99 }
100 return sum;
101}
102
104{
106 double sum = 0.0;
107 for (size_t n = 0; n < m_lattice.size(); n++) {
108 sum += theta_[n] * m_lattice[n]->cp_mole();
109 }
110 return sum;
111}
112
114{
115 return Units(1.0);
116}
117
119{
121 size_t strt = 0;
122 for (size_t n = 0; n < m_lattice.size(); n++) {
123 m_lattice[n]->getMoleFractions(c+strt);
124 strt += m_lattice[n]->nSpecies();
125 }
126}
127
129{
130 for (size_t k = 0; k < m_kk; k++) {
131 ac[k] = 1.0;
132 }
133}
134
136{
137 return 1.0;
138}
139
141{
142 return 0.0;
143}
144
146{
147 m_press = p;
148 for (size_t n = 0; n < m_lattice.size(); n++) {
149 m_lattice[n]->setPressure(m_press);
150 }
151 calcDensity();
152}
153
155{
156 double sum = 0.0;
157 for (size_t n = 0; n < m_lattice.size(); n++) {
158 sum += theta_[n] * m_lattice[n]->density();
159 }
161 return sum;
162}
163
164void LatticeSolidPhase::setMoleFractions(const double* const x)
165{
166 size_t strt = 0;
167 for (size_t n = 0; n < m_lattice.size(); n++) {
168 size_t nsp = m_lattice[n]->nSpecies();
169 m_lattice[n]->setMoleFractions(x + strt);
170 strt += nsp;
171 }
172 for (size_t k = 0; k < strt; k++) {
173 m_x[k] = x[k] / m_lattice.size();
174 }
176 calcDensity();
177}
178
179void LatticeSolidPhase::getMoleFractions(double* const x) const
180{
181 size_t strt = 0;
182 // the ifdef block should be the way we calculate this.!!!!!
184 for (size_t n = 0; n < m_lattice.size(); n++) {
185 size_t nsp = m_lattice[n]->nSpecies();
186 double sum = 0.0;
187 for (size_t k = 0; k < nsp; k++) {
188 sum += (x + strt)[k];
189 }
190 for (size_t k = 0; k < nsp; k++) {
191 (x + strt)[k] /= sum;
192 }
193
194 // At this point we can check against the mole fraction vector of the
195 // underlying LatticePhase objects and get the same answer.
196 m_lattice[n]->getMoleFractions(&m_x[strt]);
197 for (size_t k = 0; k < nsp; k++) {
198 if (fabs((x + strt)[k] - m_x[strt+k]) > 1.0E-14) {
199 throw CanteraError("LatticeSolidPhase::getMoleFractions",
200 "internal error");
201 }
202 }
203 strt += nsp;
204 }
205}
206
208{
210 size_t strt = 0;
211 for (size_t n = 0; n < m_lattice.size(); n++) {
212 size_t nlsp = m_lattice[n]->nSpecies();
213 m_lattice[n]->getChemPotentials(mu+strt);
214 strt += nlsp;
215 }
216}
217
219{
221 size_t strt = 0;
222 for (size_t n = 0; n < m_lattice.size(); n++) {
223 size_t nlsp = m_lattice[n]->nSpecies();
224 m_lattice[n]->getPartialMolarEnthalpies(hbar + strt);
225 strt += nlsp;
226 }
227}
228
230{
232 size_t strt = 0;
233 for (size_t n = 0; n < m_lattice.size(); n++) {
234 size_t nlsp = m_lattice[n]->nSpecies();
235 m_lattice[n]->getPartialMolarEntropies(sbar + strt);
236 strt += nlsp;
237 }
238}
239
241{
243 size_t strt = 0;
244 for (size_t n = 0; n < m_lattice.size(); n++) {
245 size_t nlsp = m_lattice[n]->nSpecies();
246 m_lattice[n]->getPartialMolarCp(cpbar + strt);
247 strt += nlsp;
248 }
249}
250
252{
254 size_t strt = 0;
255 for (size_t n = 0; n < m_lattice.size(); n++) {
256 size_t nlsp = m_lattice[n]->nSpecies();
257 m_lattice[n]->getPartialMolarVolumes(vbar + strt);
258 strt += nlsp;
259 }
260}
261
263{
265 size_t strt = 0;
266 for (size_t n = 0; n < m_lattice.size(); n++) {
267 m_lattice[n]->getStandardChemPotentials(mu0+strt);
268 strt += m_lattice[n]->nSpecies();
269 }
270}
271
273{
275 for (size_t n = 0; n < m_lattice.size(); n++) {
276 m_lattice[n]->getGibbs_RT_ref(grt + lkstart_[n]);
277 }
278}
279
281{
283 for (size_t k = 0; k < m_kk; k++) {
284 g[k] *= RT();
285 }
286}
287
289 const AnyMap& rootNode)
290{
291 ThermoPhase::setParameters(phaseNode, rootNode);
292 m_rootNode = rootNode;
293}
294
296{
297 if (m_input.hasKey("composition")) {
298 Composition composition = m_input["composition"].asMap<double>();
299 for (auto& [name, stoich] : composition) {
300 AnyMap& node = m_rootNode["phases"].getMapWhere("name", name);
302 }
303 setLatticeStoichiometry(composition);
304 }
305
306 setMoleFractions(m_x.data());
308}
309
311{
313 AnyMap composition;
314 for (size_t i = 0; i < m_lattice.size(); i++) {
315 composition[m_lattice[i]->name()] = theta_[i];
316 }
317 phaseNode["composition"] = std::move(composition);
318
319 // Remove fields not used in this model
320 phaseNode.erase("species");
321 vector<string> elements;
322 for (auto& el : phaseNode["elements"].asVector<string>()) {
323 if (!ba::starts_with(el, "LC_")) {
324 elements.push_back(el);
325 }
326 }
327 phaseNode["elements"] = elements;
328}
329
331 AnyMap& speciesNode) const
332{
333 // Use child lattice phases to determine species parameters so that these
334 // are set consistently
335 for (const auto& phase : m_lattice) {
336 if (phase->speciesIndex(name) != npos) {
337 phase->getSpeciesParameters(name, speciesNode);
338 break;
339 }
340 }
341}
342
343bool LatticeSolidPhase::addSpecies(shared_ptr<Species> spec)
344{
345 // Species are added from component phases in addLattice()
346 return false;
347}
348
349void LatticeSolidPhase::addLattice(shared_ptr<ThermoPhase> lattice)
350{
351 m_lattice.push_back(lattice);
352 if (lkstart_.empty()) {
353 lkstart_.push_back(0);
354 }
355 lkstart_.push_back(lkstart_.back() + lattice->nSpecies());
356
357 if (theta_.size() == 0) {
358 theta_.push_back(1.0);
359 } else {
360 theta_.push_back(0.0);
361 }
362
363 for (size_t k = 0; k < lattice->nSpecies(); k++) {
364 ThermoPhase::addSpecies(lattice->species(k));
365 vector<double> constArr(lattice->nElements());
366 const vector<double>& aws = lattice->atomicWeights();
367 for (size_t es = 0; es < lattice->nElements(); es++) {
368 addElement(lattice->elementName(es), aws[es], lattice->atomicNumber(es),
369 lattice->entropyElement298(es), lattice->elementType(es));
370 }
371 m_x.push_back(lattice->moleFraction(k));
372 tmpV_.push_back(0.0);
373 }
374}
375
377{
378 for (size_t i = 0; i < m_lattice.size(); i++) {
379 theta_[i] = getValue(comp, m_lattice[i]->name(), 0.0);
380 }
381 // Add in the lattice stoichiometry constraint
382 for (size_t i = 1; i < m_lattice.size(); i++) {
383 string econ = fmt::format("LC_{}_{}", i, name());
384 size_t m = addElement(econ, 0.0, 0, 0.0, CT_ELEM_TYPE_LATTICERATIO);
385 size_t mm = nElements();
386 for (size_t k = 0; k < m_lattice[0]->nSpecies(); k++) {
387 m_speciesComp[k * mm + m] = -theta_[0];
388 }
389 for (size_t k = 0; k < m_lattice[i]->nSpecies(); k++) {
390 size_t ks = lkstart_[i] + k;
391 m_speciesComp[ks * mm + m] = theta_[i];
392 }
393 }
394}
395
397{
398 double tnow = temperature();
399 if (m_tlast != tnow) {
400 getMoleFractions(m_x.data());
401 size_t strt = 0;
402 for (size_t n = 0; n < m_lattice.size(); n++) {
403 m_lattice[n]->setTemperature(tnow);
404 m_lattice[n]->setMoleFractions(&m_x[strt]);
405 m_lattice[n]->setPressure(m_press);
406 strt += m_lattice[n]->nSpecies();
407 }
408 m_tlast = tnow;
409 }
410}
411
413{
414 m_lattice[nn]->setMoleFractionsByName(x);
415 size_t loc = 0;
416 for (size_t n = 0; n < m_lattice.size(); n++) {
417 size_t nsp = m_lattice[n]->nSpecies();
418 double ndens = m_lattice[n]->molarDensity();
419 for (size_t k = 0; k < nsp; k++) {
420 m_x[loc] = ndens * m_lattice[n]->moleFraction(k);
421 loc++;
422 }
423 }
424 setMoleFractions(m_x.data());
425}
426
427void LatticeSolidPhase::modifyOneHf298SS(const size_t k, const double Hf298New)
428{
429 for (size_t n = 0; n < m_lattice.size(); n++) {
430 if (lkstart_[n+1] < k) {
431 size_t kk = k-lkstart_[n];
432 MultiSpeciesThermo& l_spthermo = m_lattice[n]->speciesThermo();
433 l_spthermo.modifyOneHf298(kk, Hf298New);
434 }
435 }
438}
439
441{
442 if (k != npos) {
443 for (size_t n = 0; n < m_lattice.size(); n++) {
444 if (lkstart_[n+1] < k) {
445 size_t kk = k-lkstart_[n];
446 m_lattice[n]->speciesThermo().resetHf298(kk);
447 }
448 }
449 } else {
450 for (size_t n = 0; n < m_lattice.size(); n++) {
451 m_lattice[n]->speciesThermo().resetHf298(npos);
452 }
453 }
456}
457
458} // End namespace Cantera
#define CT_ELEM_TYPE_LATTICERATIO
Constraint associated with maintaining a fixed lattice stoichiometry in a solid.
Definition Elements.h:54
Header for a simple thermodynamics model of a bulk solid phase derived from ThermoPhase,...
Header for a general species thermodynamic property manager for a phase (see MultiSpeciesThermo).
Header for factory functions to build instances of classes that manage the standard-state thermodynam...
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:427
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition AnyMap.cpp:1423
void erase(const string &key)
Erase the value held by key.
Definition AnyMap.cpp:1428
Base class for exceptions thrown by Cantera classes.
vector< shared_ptr< ThermoPhase > > m_lattice
Vector of sublattice ThermoPhase objects.
void setLatticeStoichiometry(const Composition &comp)
Set the lattice stoichiometric coefficients, .
void getStandardChemPotentials(double *mu0) const override
Get the array of standard state chemical potentials at unit activity for the species at their standar...
double enthalpy_mole() const override
Return the Molar Enthalpy. Units: J/kmol.
double logStandardConc(size_t k=0) const override
Natural logarithm of the standard concentration of the kth species.
AnyMap m_rootNode
Root node of the AnyMap which contains this phase definition.
void getPartialMolarEnthalpies(double *hbar) const override
Returns an array of partial molar enthalpies for the species in the mixture.
void getChemPotentials(double *mu) const override
Get the species chemical potentials. Units: J/kmol.
vector< double > tmpV_
Temporary vector.
vector< double > m_x
Vector of mole fractions.
void getSpeciesParameters(const string &name, AnyMap &speciesNode) const override
Get phase-specific parameters of a Species object such that an identical one could be reconstructed a...
void getGibbs_ref(double *g) const override
Returns the vector of the Gibbs function of the reference state at the current temperature of the sol...
void modifyOneHf298SS(const size_t k, const double Hf298New) override
Modify the value of the 298 K Heat of Formation of one species in the phase (J kmol-1)
vector< double > theta_
Lattice stoichiometric coefficients.
void getParameters(AnyMap &phaseNode) const override
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
void setMoleFractions(const double *const x) override
Set the mole fractions to the specified values, and then normalize them so that they sum to 1....
void getActivityConcentrations(double *c) const override
This method returns an array of generalized concentrations.
void setPressure(double p) override
Set the pressure at constant temperature. Units: Pa.
void getPartialMolarVolumes(double *vbar) const override
returns an array of partial molar volumes of the species in the solution.
void setLatticeMoleFractionsByName(int n, const string &x)
Set the Lattice mole fractions using a string.
double m_press
Current value of the pressure.
double refPressure() const override
Returns the reference pressure in Pa.
double minTemp(size_t k=npos) const override
Minimum temperature for which the thermodynamic data for the species or phase are valid.
double intEnergy_mole() const override
Return the Molar Internal Energy. Units: J/kmol.
double entropy_mole() const override
Return the Molar Entropy. Units: J/kmol/K.
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
double cp_mole() const override
Return the constant pressure heat capacity. Units: J/kmol/K.
Units standardConcentrationUnits() const override
Returns the units of the "standard concentration" for this phase.
void getPartialMolarCp(double *cpbar) const override
Returns an array of partial molar Heat Capacities at constant pressure of the species in the solution...
double gibbs_mole() const override
Return the Molar Gibbs energy. Units: J/kmol.
double standardConcentration(size_t k=0) const override
Return the standard concentration for the kth species.
void addLattice(shared_ptr< ThermoPhase > lattice)
Add a lattice to this phase.
double calcDensity()
Calculate the density of the solid mixture.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
void _updateThermo() const
Update the reference thermodynamic functions.
void getGibbs_RT_ref(double *grt) const override
Returns the vector of nondimensional Gibbs Free Energies of the reference state at the current temper...
void getActivityCoefficients(double *ac) const override
Get the array of non-dimensional molar-based activity coefficients at the current solution temperatur...
void setParameters(const AnyMap &phaseNode, const AnyMap &rootNode=AnyMap()) override
Set equation of state parameters from an AnyMap phase description.
void resetHf298(const size_t k=npos) override
Restore the original heat of formation of one or more species.
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies of the species in the solution.
double maxTemp(size_t k=npos) const override
Maximum temperature for which the thermodynamic data for the species are valid.
A species thermodynamic property manager for a phase.
virtual void modifyOneHf298(const size_t k, const double Hf298New)
Modify the value of the 298 K Heat of Formation of the standard state of one species in the phase (J ...
void assignDensity(const double density_)
Set the internally stored constant density (kg/m^3) of the phase.
Definition Phase.cpp:597
virtual void setMoleFractions(const double *const x)
Set the mole fractions to the specified values.
Definition Phase.cpp:289
vector< double > m_speciesComp
Atomic composition of the species.
Definition Phase.h:851
size_t m_kk
Number of species in the phase.
Definition Phase.h:842
double temperature() const
Temperature (K).
Definition Phase.h:562
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition Phase.cpp:434
size_t nElements() const
Number of elements.
Definition Phase.cpp:30
size_t addElement(const string &symbol, double weight=-12345.0, int atomicNumber=0, double entropy298=ENTROPY298_UNKNOWN, int elem_type=CT_ELEM_TYPE_ABSPOS)
Add an element.
Definition Phase.cpp:635
string name() const
Return the name of the phase.
Definition Phase.cpp:20
virtual void setParameters(const AnyMap &phaseNode, const AnyMap &rootNode=AnyMap())
Set equation of state parameters from an AnyMap phase description.
virtual void getParameters(AnyMap &phaseNode) const
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
double RT() const
Return the Gas Constant multiplied by the current temperature.
double m_tlast
last value of the temperature processed by reference state
virtual void initThermo()
Initialize the ThermoPhase object after all species have been set up.
void invalidateCache() override
Invalidate any cached values which are normally updated only when a change in state is detected.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
AnyMap m_input
Data supplied via setParameters.
A representation of the units associated with a dimensional quantity.
Definition Units.h:35
shared_ptr< ThermoPhase > newThermo(const AnyMap &phaseNode, const AnyMap &rootNode)
Create a new ThermoPhase object and initialize it.
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
const double BigNumber
largest number to compare to inf.
Definition ct_defs.h:160
const U & getValue(const map< T, U > &m, const T &key, const U &default_val)
Const accessor for a value in a map.
Definition utilities.h:190
map< string, double > Composition
Map from string names to doubles.
Definition ct_defs.h:177
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...