Cantera  3.1.0a1
IdealSolidSolnPhase.cpp
Go to the documentation of this file.
1 /**
2  * @file IdealSolidSolnPhase.cpp Implementation file for an ideal solid
3  * solution model with incompressible thermodynamics (see @ref
4  * thermoprops and @link Cantera::IdealSolidSolnPhase
5  * IdealSolidSolnPhase@endlink).
6  */
7 
8 // This file is part of Cantera. See License.txt in the top-level directory or
9 // at https://cantera.org/license.txt for license and copyright information.
10 
13 #include "cantera/thermo/Species.h"
15 #include "cantera/base/utilities.h"
16 
17 namespace Cantera
18 {
19 
20 IdealSolidSolnPhase::IdealSolidSolnPhase(const string& inputFile, const string& id_)
21 {
22  initThermoFile(inputFile, id_);
23 }
24 
25 // Molar Thermodynamic Properties of the Solution
26 
28 {
29  double htp = RT() * mean_X(enthalpy_RT_ref());
30  return htp + (pressure() - m_Pref)/molarDensity();
31 }
32 
34 {
35  return GasConstant * (mean_X(entropy_R_ref()) - sum_xlogx());
36 }
37 
39 {
40  double Pv = (pressure() - m_Pref)/molarDensity();
41  return RT() * (mean_X(gibbs_RT_ref()) + sum_xlogx()) + Pv;
42 }
43 
45 {
46  return GasConstant * mean_X(cp_R_ref());
47 }
48 
49 // Mechanical Equation of State
50 
52 {
53  // Calculate the molarVolume of the solution (m**3 kmol-1)
54  double v_mol = mean_X(m_speciesMolarVolume);
55 
56  // Set the density in the parent object directly, by calling the
57  // Phase::assignDensity() function.
59 }
60 
62 {
63  m_Pcurrent = p;
64  calcDensity();
65 }
66 
68 {
70  calcDensity();
71 }
72 
73 // Chemical Potentials and Activities
74 
76 {
77  if (m_formGC == 0) {
78  return Units(1.0); // dimensionless
79  } else {
80  // kmol/m^3 for bulk phases
81  return Units(1.0, 0, -static_cast<double>(nDim()), 0, 0, 0, 1);
82  }
83 }
84 
86 {
88  switch (m_formGC) {
89  case 0:
90  break;
91  case 1:
92  for (size_t k = 0; k < m_kk; k++) {
93  c[k] /= m_speciesMolarVolume[k];
94  }
95  break;
96  case 2:
97  for (size_t k = 0; k < m_kk; k++) {
98  c[k] /= m_speciesMolarVolume[m_kk-1];
99  }
100  break;
101  }
102 }
103 
105 {
106  switch (m_formGC) {
107  case 0:
108  return 1.0;
109  case 1:
110  return 1.0 / m_speciesMolarVolume[k];
111  case 2:
112  return 1.0/m_speciesMolarVolume[m_kk-1];
113  }
114  return 0.0;
115 }
116 
118 {
119  for (size_t k = 0; k < m_kk; k++) {
120  ac[k] = 1.0;
121  }
122 }
123 
125 {
126  double delta_p = m_Pcurrent - m_Pref;
127  const vector<double>& g_RT = gibbs_RT_ref();
128  for (size_t k = 0; k < m_kk; k++) {
129  double xx = std::max(SmallNumber, moleFraction(k));
130  mu[k] = RT() * (g_RT[k] + log(xx))
131  + delta_p * m_speciesMolarVolume[k];
132  }
133 }
134 
135 // Partial Molar Properties
136 
138 {
139  const vector<double>& _h = enthalpy_RT_ref();
140  double delta_p = m_Pcurrent - m_Pref;
141  for (size_t k = 0; k < m_kk; k++) {
142  hbar[k] = _h[k]*RT() + delta_p * m_speciesMolarVolume[k];
143  }
144  // scale(_h.begin(), _h.end(), hbar, RT());
145 }
146 
148 {
149  const vector<double>& _s = entropy_R_ref();
150  for (size_t k = 0; k < m_kk; k++) {
151  double xx = std::max(SmallNumber, moleFraction(k));
152  sbar[k] = GasConstant * (_s[k] - log(xx));
153  }
154 }
155 
156 void IdealSolidSolnPhase::getPartialMolarCp(double* cpbar) const
157 {
158  getCp_R(cpbar);
159  for (size_t k = 0; k < m_kk; k++) {
160  cpbar[k] *= GasConstant;
161  }
162 }
163 
165 {
166  getStandardVolumes(vbar);
167 }
168 
169 // Properties of the Standard State of the Species in the Solution
170 
171 void IdealSolidSolnPhase::getPureGibbs(double* gpure) const
172 {
173  const vector<double>& gibbsrt = gibbs_RT_ref();
174  double delta_p = (m_Pcurrent - m_Pref);
175  for (size_t k = 0; k < m_kk; k++) {
176  gpure[k] = RT() * gibbsrt[k] + delta_p * m_speciesMolarVolume[k];
177  }
178 }
179 
180 void IdealSolidSolnPhase::getGibbs_RT(double* grt) const
181 {
182  const vector<double>& gibbsrt = gibbs_RT_ref();
183  double delta_prt = (m_Pcurrent - m_Pref)/ RT();
184  for (size_t k = 0; k < m_kk; k++) {
185  grt[k] = gibbsrt[k] + delta_prt * m_speciesMolarVolume[k];
186  }
187 }
188 
189 void IdealSolidSolnPhase::getEnthalpy_RT(double* hrt) const
190 {
191  const vector<double>& _h = enthalpy_RT_ref();
192  double delta_prt = (m_Pcurrent - m_Pref) / RT();
193  for (size_t k = 0; k < m_kk; k++) {
194  hrt[k] = _h[k] + delta_prt * m_speciesMolarVolume[k];
195  }
196 }
197 
198 void IdealSolidSolnPhase::getEntropy_R(double* sr) const
199 {
200  const vector<double>& _s = entropy_R_ref();
201  copy(_s.begin(), _s.end(), sr);
202 }
203 
205 {
206  const vector<double>& _h = enthalpy_RT_ref();
207  double prefrt = m_Pref / RT();
208  for (size_t k = 0; k < m_kk; k++) {
209  urt[k] = _h[k] - prefrt * m_speciesMolarVolume[k];
210  }
211 }
212 
213 void IdealSolidSolnPhase::getCp_R(double* cpr) const
214 {
215  const vector<double>& _cpr = cp_R_ref();
216  copy(_cpr.begin(), _cpr.end(), cpr);
217 }
218 
220 {
221  copy(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), vol);
222 }
223 
224 // Thermodynamic Values for the Species Reference States
225 
227 {
228  _updateThermo();
229  for (size_t k = 0; k != m_kk; k++) {
230  hrt[k] = m_h0_RT[k];
231  }
232 }
233 
235 {
236  _updateThermo();
237  for (size_t k = 0; k != m_kk; k++) {
238  grt[k] = m_g0_RT[k];
239  }
240 }
241 
243 {
244  _updateThermo();
245  double tmp = RT();
246  for (size_t k = 0; k != m_kk; k++) {
247  g[k] = tmp * m_g0_RT[k];
248  }
249 }
250 
252 {
253  const vector<double>& _h = enthalpy_RT_ref();
254  double prefrt = m_Pref / RT();
255  for (size_t k = 0; k < m_kk; k++) {
256  urt[k] = _h[k] - prefrt * m_speciesMolarVolume[k];
257  }
258 }
259 
261 {
262  _updateThermo();
263  for (size_t k = 0; k != m_kk; k++) {
264  er[k] = m_s0_R[k];
265  }
266 }
267 
268 void IdealSolidSolnPhase::getCp_R_ref(double* cpr) const
269 {
270  _updateThermo();
271  for (size_t k = 0; k != m_kk; k++) {
272  cpr[k] = m_cp0_R[k];
273  }
274 }
275 
276 const vector<double>& IdealSolidSolnPhase::enthalpy_RT_ref() const
277 {
278  _updateThermo();
279  return m_h0_RT;
280 }
281 
282 const vector<double>& IdealSolidSolnPhase::entropy_R_ref() const
283 {
284  _updateThermo();
285  return m_s0_R;
286 }
287 
288 // Utility Functions
289 
290 bool IdealSolidSolnPhase::addSpecies(shared_ptr<Species> spec)
291 {
292  bool added = ThermoPhase::addSpecies(spec);
293  if (added) {
294  if (m_kk == 1) {
295  // Obtain the reference pressure by calling the ThermoPhase function
296  // refPressure, which in turn calls the species thermo reference
297  // pressure function of the same name.
298  m_Pref = refPressure();
299  }
300 
301  m_h0_RT.push_back(0.0);
302  m_g0_RT.push_back(0.0);
303  m_expg0_RT.push_back(0.0);
304  m_cp0_R.push_back(0.0);
305  m_s0_R.push_back(0.0);
306  m_pp.push_back(0.0);
307  if (spec->input.hasKey("equation-of-state")) {
308  auto& eos = spec->input["equation-of-state"].getMapWhere("model", "constant-volume");
309  double mv;
310  if (eos.hasKey("density")) {
311  mv = molecularWeight(m_kk-1) / eos.convert("density", "kg/m^3");
312  } else if (eos.hasKey("molar-density")) {
313  mv = 1.0 / eos.convert("molar-density", "kmol/m^3");
314  } else if (eos.hasKey("molar-volume")) {
315  mv = eos.convert("molar-volume", "m^3/kmol");
316  } else {
317  throw CanteraError("IdealSolidSolnPhase::addSpecies",
318  "equation-of-state entry for species '{}' is missing "
319  "'density', 'molar-volume', or 'molar-density' "
320  "specification", spec->name);
321  }
322  m_speciesMolarVolume.push_back(mv);
323  } else {
324  throw CanteraError("IdealSolidSolnPhase::addSpecies",
325  "Molar volume not specified for species '{}'", spec->name);
326  }
327  if (ready()) {
328  calcDensity();
329  }
330  }
331  return added;
332 }
333 
335 {
336  if (m_input.hasKey("standard-concentration-basis")) {
337  setStandardConcentrationModel(m_input["standard-concentration-basis"].asString());
338  }
340 }
341 
343 {
344  ThermoPhase::getParameters(phaseNode);
345  if (m_formGC == 1) {
346  phaseNode["standard-concentration-basis"] = "species-molar-volume";
347  } else if (m_formGC == 2) {
348  phaseNode["standard-concentration-basis"] = "solvent-molar-volume";
349  }
350 }
351 
353  AnyMap& speciesNode) const
354 {
356  size_t k = speciesIndex(name);
357  const auto S = species(k);
358  auto& eosNode = speciesNode["equation-of-state"].getMapWhere(
359  "model", "constant-volume", true);
360  // Output volume information in a form consistent with the input
361  if (S->input.hasKey("equation-of-state")) {
362  auto& eosIn = S->input["equation-of-state"];
363  if (eosIn.hasKey("density")) {
364  eosNode["density"].setQuantity(
365  molecularWeight(k) / m_speciesMolarVolume[k], "kg/m^3");
366  } else if (eosIn.hasKey("molar-density")) {
367  eosNode["molar-density"].setQuantity(1.0 / m_speciesMolarVolume[k],
368  "kmol/m^3");
369  } else {
370  eosNode["molar-volume"].setQuantity(m_speciesMolarVolume[k],
371  "m^3/kmol");
372  }
373  } else {
374  eosNode["molar-volume"].setQuantity(m_speciesMolarVolume[k],
375  "m^3/kmol");
376  }
377 }
378 
379 void IdealSolidSolnPhase::setToEquilState(const double* mu_RT)
380 {
381  const vector<double>& grt = gibbs_RT_ref();
382 
383  // Within the method, we protect against inf results if the exponent is too
384  // high.
385  //
386  // If it is too low, we set the partial pressure to zero. This capability is
387  // needed by the elemental potential method.
388  double pres = 0.0;
389  double m_p0 = refPressure();
390  for (size_t k = 0; k < m_kk; k++) {
391  double tmp = -grt[k] + mu_RT[k];
392  if (tmp < -600.) {
393  m_pp[k] = 0.0;
394  } else if (tmp > 500.0) {
395  // Protect against inf results if the exponent is too high
396  double tmp2 = tmp / 500.;
397  tmp2 *= tmp2;
398  m_pp[k] = m_p0 * exp(500.) * tmp2;
399  } else {
400  m_pp[k] = m_p0 * exp(tmp);
401  }
402  pres += m_pp[k];
403  }
404  // set state
405  setMoleFractions(m_pp.data());
406  setPressure(pres);
407 }
408 
410 {
411  if (caseInsensitiveEquals(model, "unity")) {
412  m_formGC = 0;
413  } else if (caseInsensitiveEquals(model, "species-molar-volume")
414  || caseInsensitiveEquals(model, "molar_volume")) {
415  m_formGC = 1;
416  } else if (caseInsensitiveEquals(model, "solvent-molar-volume")
417  || caseInsensitiveEquals(model, "solvent_volume")) {
418  m_formGC = 2;
419  } else {
420  throw CanteraError("IdealSolidSolnPhase::setStandardConcentrationModel",
421  "Unknown standard concentration model '{}'", model);
422  }
423 }
424 
426 {
427  return m_speciesMolarVolume[k];
428 }
429 
431 {
432  copy(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), smv);
433 }
434 
436 {
437  double tnow = temperature();
438  if (m_tlast != tnow) {
439 
440  // Update the thermodynamic functions of the reference state.
441  m_spthermo.update(tnow, m_cp0_R.data(), m_h0_RT.data(), m_s0_R.data());
442  m_tlast = tnow;
443  for (size_t k = 0; k < m_kk; k++) {
444  m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
445  }
446  m_tlast = tnow;
447  }
448 }
449 
450 } // end namespace Cantera
Header file for an ideal solid solution model with incompressible thermodynamics (see Thermodynamic P...
Declaration for class Cantera::Species.
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
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
const vector< double > & entropy_R_ref() const
Returns a reference to the vector of nondimensional enthalpies of the reference state at the current ...
double enthalpy_mole() const override
Molar enthalpy of the solution.
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.
double pressure() const override
Pressure.
vector< double > m_g0_RT
Vector containing the species reference Gibbs functions at T = m_tlast.
const vector< double > & cp_R_ref() const
Returns a reference to the vector of nondimensional enthalpies of the reference state at the current ...
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...
const vector< double > & gibbs_RT_ref() const
Returns a reference to the vector of nondimensional enthalpies of the reference state at the current ...
void getEntropy_R(double *sr) const override
Get the nondimensional Entropies for the species standard states at the current T and P of the soluti...
vector< double > m_h0_RT
Vector containing the species reference enthalpies at T = m_tlast.
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...
double speciesMolarVolume(int k) const
Report the molar volume of species k.
vector< double > m_pp
Temporary array used in equilibrium calculations.
double m_Pref
Value of the reference pressure for all species in this phase.
void getCp_R(double *cpr) const override
Get the nondimensional heat capacity at constant pressure function for the species standard states at...
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 getActivityConcentrations(double *c) const override
This method returns the array of generalized concentrations.
void getSpeciesMolarVolumes(double *smv) const
Fill in a return vector containing the species molar volumes.
void setPressure(double p) override
Set the pressure at constant temperature.
void getPartialMolarVolumes(double *vbar) const override
returns an array of partial molar volumes of the species in the solution.
void getPureGibbs(double *gpure) const override
Get the Gibbs functions for the pure species at the current T and P of the solution.
double standardConcentration(size_t k) const override
The standard concentration used to normalize the generalized concentration.
void setStandardConcentrationModel(const string &model)
Set the form for the standard and generalized concentrations.
void getIntEnergy_RT_ref(double *urt) const override
Returns the vector of nondimensional internal Energies of the reference state at the current temperat...
void getEnthalpy_RT(double *hrt) const override
Get the array of nondimensional Enthalpy functions for the standard state species at the current T an...
void getEntropy_R_ref(double *er) const override
Returns the vector of nondimensional entropies of the reference state at the current temperature of t...
vector< double > m_s0_R
Vector containing the species reference entropies at T = m_tlast.
void getGibbs_RT(double *grt) const override
Get the nondimensional Gibbs function for the species standard states at the current T and P of the s...
double entropy_mole() const override
Molar entropy of the solution.
int m_formGC
The standard concentrations can have one of three different forms: 0 = 'unity', 1 = 'species-molar-vo...
vector< double > m_speciesMolarVolume
Vector of molar volumes for each species in the solution.
void getCp_R_ref(double *cprt) const override
Returns the vector of nondimensional constant pressure heat capacities of the reference state at the ...
void getStandardVolumes(double *vol) const override
Get the molar volumes of the species standard states at the current T and P of the solution.
double cp_mole() const override
Molar heat capacity at constant pressure of the solution.
void getIntEnergy_RT(double *urt) const override
Returns the vector of nondimensional Internal Energies of the standard state species at the current T...
Units standardConcentrationUnits() const override
Returns the units of the "standard concentration" for this phase.
IdealSolidSolnPhase(const string &infile="", const string &id="")
Construct and initialize an IdealSolidSolnPhase ThermoPhase object directly from an input file.
void getPartialMolarCp(double *cpbar) const override
Returns an array of partial molar Heat Capacities at constant pressure of the species in the solution...
void compositionChanged() override
Apply changes to the state which are needed after the composition changes.
double gibbs_mole() const override
Molar Gibbs free energy of the solution.
vector< double > m_cp0_R
Vector containing the species reference constant pressure heat capacities at T = m_tlast.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
virtual void _updateThermo() const
This function gets called for every call to functions in this class.
void setToEquilState(const double *mu_RT) override
This method is used by the ChemEquil equilibrium solver.
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 species activity coefficients.
vector< double > m_expg0_RT
Vector containing the species reference exp(-G/RT) functions at T = m_tlast.
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies of the species in the solution.
double m_Pcurrent
m_Pcurrent = The current pressure Since the density isn't a function of pressure, but only of the mol...
void getEnthalpy_RT_ref(double *hrt) const override
Returns the vector of nondimensional enthalpies of the reference state at the current temperature of ...
virtual void calcDensity()
Calculate the density of the mixture using the partial molar volumes and mole fractions as input.
const vector< double > & enthalpy_RT_ref() const
Returns a reference to the vector of nondimensional enthalpies of the reference state at the current ...
virtual void update(double T, double *cp_R, double *h_RT, double *s_R) const
Compute the reference-state properties for all species.
virtual double molarDensity() const
Molar density (kmol/m^3).
Definition: Phase.cpp:576
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
size_t m_kk
Number of species in the phase.
Definition: Phase.h:842
size_t nDim() const
Returns the number of spatial dimensions (1, 2, or 3)
Definition: Phase.h:546
double temperature() const
Temperature (K).
Definition: Phase.h:562
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:655
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:434
double sum_xlogx() const
Evaluate .
Definition: Phase.cpp:626
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition: Phase.cpp:129
double moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition: Phase.cpp:439
virtual void compositionChanged()
Apply changes to the state which are needed after the composition changes.
Definition: Phase.cpp:905
double mean_X(const double *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition: Phase.cpp:616
virtual bool ready() const
Returns a bool indicating whether the object is ready for use.
Definition: Phase.cpp:885
double molecularWeight(size_t k) const
Molecular weight of species k.
Definition: Phase.cpp:383
shared_ptr< Species > species(const string &name) const
Return the Species object for the named species.
Definition: Phase.cpp:856
string name() const
Return the name of the phase.
Definition: Phase.cpp:20
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 initThermoFile(const string &inputFile, const string &id)
Initialize a ThermoPhase object using an input file.
virtual void getSpeciesParameters(const string &name, AnyMap &speciesNode) const
Get phase-specific parameters of a Species object such that an identical one could be reconstructed a...
Definition: ThermoPhase.h:1831
MultiSpeciesThermo m_spthermo
Pointer to the calculation manager for species reference-state thermodynamic properties.
Definition: ThermoPhase.h:1962
virtual double refPressure() const
Returns the reference pressure in Pa.
Definition: ThermoPhase.h:436
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
bool caseInsensitiveEquals(const string &input, const string &test)
Case insensitive equality predicate.
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition: ct_defs.h:120
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:564
const double SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:158
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...