Cantera  3.1.0a1
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 
17 #include "cantera/base/utilities.h"
18 
19 #include <boost/algorithm/string.hpp>
20 
21 namespace ba = boost::algorithm;
22 
23 namespace Cantera
24 {
25 
26 double 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 
42 double 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 {
65  _updateThermo();
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 {
75  _updateThermo();
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 {
85  _updateThermo();
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 {
95  _updateThermo();
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 {
105  _updateThermo();
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 {
120  _updateThermo();
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 
164 void 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 
179 void 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 {
209  _updateThermo();
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 {
220  _updateThermo();
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 {
231  _updateThermo();
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 
240 void LatticeSolidPhase::getPartialMolarCp(double* cpbar) const
241 {
242  _updateThermo();
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 {
253  _updateThermo();
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 {
264  _updateThermo();
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 
272 void LatticeSolidPhase::getGibbs_RT_ref(double* grt) const
273 {
274  _updateThermo();
275  for (size_t n = 0; n < m_lattice.size(); n++) {
276  m_lattice[n]->getGibbs_RT_ref(grt + lkstart_[n]);
277  }
278 }
279 
280 void LatticeSolidPhase::getGibbs_ref(double* g) const
281 {
282  getGibbs_RT_ref(g);
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 {
312  ThermoPhase::getParameters(phaseNode);
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 
343 bool LatticeSolidPhase::addSpecies(shared_ptr<Species> spec)
344 {
345  // Species are added from component phases in addLattice()
346  return false;
347 }
348 
349 void 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 
427 void 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  }
436  invalidateCache();
437  _updateThermo();
438 }
439 
440 void LatticeSolidPhase::resetHf298(const size_t k)
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  }
454  invalidateCache();
455  _updateThermo();
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.
Definition: ctexceptions.h:66
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.
Definition: ThermoPhase.h:1062
double m_tlast
last value of the temperature processed by reference state
Definition: ThermoPhase.h:1985
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.
Definition: ThermoPhase.h:1966
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 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
const double BigNumber
largest number to compare to inf.
Definition: ct_defs.h:160
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...