Cantera  3.1.0b1
Loading...
Searching...
No Matches
MolalityVPSSTP.cpp
Go to the documentation of this file.
1/**
2 * @file MolalityVPSSTP.cpp
3 * Definitions for intermediate ThermoPhase object for phases which
4 * employ molality based activity coefficient formulations
5 * (see @ref thermoprops
6 * and class @link Cantera::MolalityVPSSTP MolalityVPSSTP@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
15
16#include <fstream>
17
18namespace Cantera
19{
20
22{
23 // Change the default to be that charge neutrality in the phase is necessary
24 // condition for the proper specification of thermodynamic functions within
25 // the phase
27}
28
29// -------------- Utilities -------------------------------
30
31void MolalityVPSSTP::setpHScale(const int pHscaleType)
32{
33 m_pHScalingType = pHscaleType;
34 if (pHscaleType != PHSCALE_PITZER && pHscaleType != PHSCALE_NBS) {
35 throw CanteraError("MolalityVPSSTP::setpHScale",
36 "Unknown scale type: {}", pHscaleType);
37 }
38}
39
41{
42 return m_pHScalingType;
43}
44
45void MolalityVPSSTP::setMoleFSolventMin(double xmolSolventMIN)
46{
47 if (xmolSolventMIN <= 0.0) {
48 throw CanteraError("MolalityVPSSTP::setMoleFSolventMin ", "trouble");
49 } else if (xmolSolventMIN > 0.9) {
50 throw CanteraError("MolalityVPSSTP::setMoleFSolventMin ", "trouble");
51 }
52 m_xmolSolventMIN = xmolSolventMIN;
53}
54
56{
57 return m_xmolSolventMIN;
58}
59
61{
63 double xmolSolvent = std::max(m_molalities[0], m_xmolSolventMIN);
64 double denomInv = 1.0/ (m_Mnaught * xmolSolvent);
65 for (size_t k = 0; k < m_kk; k++) {
66 m_molalities[k] *= denomInv;
67 }
68}
69
70void MolalityVPSSTP::getMolalities(double* const molal) const
71{
73 for (size_t k = 0; k < m_kk; k++) {
74 molal[k] = m_molalities[k];
75 }
76}
77
78void MolalityVPSSTP::setMolalities(const double* const molal)
79{
80 double Lsum = 1.0 / m_Mnaught;
81 for (size_t k = 1; k < m_kk; k++) {
82 m_molalities[k] = molal[k];
83 Lsum += molal[k];
84 }
85 double tmp = 1.0 / Lsum;
86 m_molalities[0] = tmp / m_Mnaught;
87 double sum = m_molalities[0];
88 for (size_t k = 1; k < m_kk; k++) {
89 m_molalities[k] = tmp * molal[k];
90 sum += m_molalities[k];
91 }
92 if (sum != 1.0) {
93 tmp = 1.0 / sum;
94 for (size_t k = 0; k < m_kk; k++) {
95 m_molalities[k] *= tmp;
96 }
97 }
99
100 // Essentially we don't trust the input: We calculate the molalities from
101 // the mole fractions that we just obtained.
103}
104
106{
107 // HKM -> Might need to be more complicated here, setting neutrals so that
108 // the existing mole fractions are preserved.
109
110 // Get a vector of mole fractions
111 vector<double> mf(m_kk, 0.0);
112 getMoleFractions(mf.data());
113 double xmolSmin = std::max(mf[0], m_xmolSolventMIN);
114 for (size_t k = 0; k < m_kk; k++) {
115 double mol_k = getValue(mMap, speciesName(k), 0.0);
116 if (mol_k > 0) {
117 mf[k] = mol_k * m_Mnaught * xmolSmin;
118 }
119 }
120
121 // check charge neutrality
122 size_t largePos = npos;
123 double cPos = 0.0;
124 size_t largeNeg = npos;
125 double cNeg = 0.0;
126 double sum = 0.0;
127 for (size_t k = 0; k < m_kk; k++) {
128 double ch = charge(k);
129 if (mf[k] > 0.0) {
130 if (ch > 0.0 && ch * mf[k] > cPos) {
131 largePos = k;
132 cPos = ch * mf[k];
133 }
134 if (ch < 0.0 && fabs(ch) * mf[k] > cNeg) {
135 largeNeg = k;
136 cNeg = fabs(ch) * mf[k];
137 }
138 }
139 sum += mf[k] * ch;
140 }
141 if (sum != 0.0) {
142 if (sum > 0.0) {
143 if (cPos > sum) {
144 mf[largePos] -= sum / charge(largePos);
145 } else {
146 throw CanteraError("MolalityVPSSTP:setMolalitiesbyName",
147 "unbalanced charges");
148 }
149 } else {
150 if (cNeg > (-sum)) {
151 mf[largeNeg] -= (-sum) / fabs(charge(largeNeg));
152 } else {
153 throw CanteraError("MolalityVPSSTP:setMolalitiesbyName",
154 "unbalanced charges");
155 }
156 }
157 }
158 sum = 0.0;
159 for (size_t k = 0; k < m_kk; k++) {
160 sum += mf[k];
161 }
162 sum = 1.0/sum;
163 for (size_t k = 0; k < m_kk; k++) {
164 mf[k] *= sum;
165 }
166 setMoleFractions(mf.data());
167
168 // After we formally set the mole fractions, we calculate the molalities
169 // again and store it in this object.
171}
172
174{
177}
178
179// - Activities, Standard States, Activity Concentrations -----------
180
182{
184}
185
187{
188 throw NotImplementedError("MolalityVPSSTP::getActivityConcentrations");
189}
190
192{
193 throw NotImplementedError("MolalityVPSSTP::standardConcentration");
194}
195
196void MolalityVPSSTP::getActivities(double* ac) const
197{
198 throw NotImplementedError("MolalityVPSSTP::getActivities");
199}
200
202{
204 double xmolSolvent = std::max(moleFraction(0), m_xmolSolventMIN);
205 for (size_t k = 1; k < m_kk; k++) {
206 ac[k] /= xmolSolvent;
207 }
208}
209
211{
213 applyphScale(acMolality);
214}
215
217{
218 // First, we calculate the activities all over again
219 vector<double> act(m_kk);
220 getActivities(act.data());
221
222 // Then, we calculate the sum of the solvent molalities
223 double sum = 0;
224 for (size_t k = 1; k < m_kk; k++) {
225 sum += std::max(m_molalities[k], 0.0);
226 }
227 double oc = 1.0;
228 if (sum > 1.0E-200) {
229 oc = - log(act[0]) / (m_Mnaught * sum);
230 }
231 return oc;
232}
233
234void MolalityVPSSTP::setState_TPM(double t, double p, const double* const molalities)
235{
236 setMolalities(molalities);
237 setState_TP(t, p);
238}
239
240void MolalityVPSSTP::setState_TPM(double t, double p, const Composition& m)
241{
243 setState_TP(t, p);
244}
245
246void MolalityVPSSTP::setState_TPM(double t, double p, const string& m)
247{
249 setState_TP(t, p);
250}
251
253 AnyValue molalities;
254 if (state.hasKey("molalities")) {
255 molalities = state["molalities"];
256 } else if (state.hasKey("M")) {
257 molalities = state["M"];
258 }
259
260 if (molalities.is<string>()) {
261 setMolalitiesByName(molalities.asString());
262 } else if (molalities.is<AnyMap>()) {
263 setMolalitiesByName(molalities.asMap<double>());
264 }
265
267}
268
270{
272
273 // Find the Cl- species
275}
276
278{
279 throw NotImplementedError("MolalityVPSSTP::getUnscaledMolalityActivityCoefficients");
280}
281
282void MolalityVPSSTP::applyphScale(double* acMolality) const
283{
284 throw NotImplementedError("MolalityVPSSTP::applyphScale");
285}
286
288{
289 size_t indexCLM = npos;
290 size_t eCl = npos;
291 size_t eE = npos;
292 size_t ne = nElements();
293 for (size_t e = 0; e < ne; e++) {
294 string sn = elementName(e);
295 if (sn == "Cl" || sn == "CL") {
296 eCl = e;
297 break;
298 }
299 }
300 // We have failed if we can't find the Cl element index
301 if (eCl == npos) {
302 return npos;
303 }
304 for (size_t e = 0; e < ne; e++) {
305 string sn = elementName(e);
306 if (sn == "E" || sn == "e") {
307 eE = e;
308 break;
309 }
310 }
311 // We have failed if we can't find the E element index
312 if (eE == npos) {
313 return npos;
314 }
315 for (size_t k = 1; k < m_kk; k++) {
316 double nCl = nAtoms(k, eCl);
317 if (nCl != 1.0) {
318 continue;
319 }
320 double nE = nAtoms(k, eE);
321 if (nE != 1.0) {
322 continue;
323 }
324 for (size_t e = 0; e < ne; e++) {
325 if (e != eE && e != eCl) {
326 double nA = nAtoms(k, e);
327 if (nA != 0.0) {
328 continue;
329 }
330 }
331 }
332 string sn = speciesName(k);
333 if (sn != "Cl-" && sn != "CL-") {
334 continue;
335 }
336
337 indexCLM = k;
338 break;
339 }
340 return indexCLM;
341}
342
343bool MolalityVPSSTP::addSpecies(shared_ptr<Species> spec)
344{
345 bool added = VPStandardStateTP::addSpecies(spec);
346 if (added) {
347 if (m_kk == 1) {
348 // The solvent defaults to species 0
350 m_Mnaught = m_weightSolvent / 1000.;
351 }
352 m_molalities.push_back(0.0);
353 }
354 return added;
355}
356
357string MolalityVPSSTP::report(bool show_thermo, double threshold) const
358{
359 fmt::memory_buffer b;
360 try {
361 if (name() != "") {
362 fmt_append(b, "\n {}:\n", name());
363 }
364 fmt_append(b, "\n");
365 fmt_append(b, " temperature {:12.6g} K\n", temperature());
366 fmt_append(b, " pressure {:12.6g} Pa\n", pressure());
367 fmt_append(b, " density {:12.6g} kg/m^3\n", density());
368 fmt_append(b, " mean mol. weight {:12.6g} amu\n", meanMolecularWeight());
369
370 double phi = electricPotential();
371 fmt_append(b, " potential {:12.6g} V\n", phi);
372
373 vector<double> x(m_kk);
374 vector<double> molal(m_kk);
375 vector<double> mu(m_kk);
376 vector<double> muss(m_kk);
377 vector<double> acMolal(m_kk);
378 vector<double> actMolal(m_kk);
379 getMoleFractions(&x[0]);
380 getMolalities(&molal[0]);
381 getChemPotentials(&mu[0]);
384 getActivities(&actMolal[0]);
385
386 size_t iHp = speciesIndex("H+");
387 if (iHp != npos) {
388 double pH = -log(actMolal[iHp]) / log(10.0);
389 fmt_append(b,
390 " pH {:12.4g}\n", pH);
391 }
392
393 if (show_thermo) {
394 fmt_append(b, "\n");
395 fmt_append(b, " 1 kg 1 kmol\n");
396 fmt_append(b, " ----------- ------------\n");
397 fmt_append(b, " enthalpy {:12.6g} {:12.4g} J\n",
399 fmt_append(b, " internal energy {:12.6g} {:12.4g} J\n",
401 fmt_append(b, " entropy {:12.6g} {:12.4g} J/K\n",
403 fmt_append(b, " Gibbs function {:12.6g} {:12.4g} J\n",
405 fmt_append(b, " heat capacity c_p {:12.6g} {:12.4g} J/K\n",
406 cp_mass(), cp_mole());
407 try {
408 fmt_append(b, " heat capacity c_v {:12.6g} {:12.4g} J/K\n",
409 cv_mass(), cv_mole());
410 } catch (NotImplementedError&) {
411 fmt_append(b, " heat capacity c_v <not implemented>\n");
412 }
413 }
414
415 fmt_append(b, "\n");
416 int nMinor = 0;
417 double xMinor = 0.0;
418 if (show_thermo) {
419 fmt_append(b, " X "
420 " Molalities Chem.Pot. ChemPotSS ActCoeffMolal\n");
421 fmt_append(b, " "
422 " (J/kmol) (J/kmol)\n");
423 fmt_append(b, " ------------- "
424 " ------------ ------------ ------------ ------------\n");
425 for (size_t k = 0; k < m_kk; k++) {
426 if (x[k] > threshold) {
427 if (x[k] > SmallNumber) {
428 fmt_append(b,
429 "{:>18s} {:12.6g} {:12.6g} {:12.6g} "
430 "{:12.6g} {:12.6g}\n", speciesName(k),
431 x[k], molal[k], mu[k], muss[k], acMolal[k]);
432 } else {
433 fmt_append(b,
434 "{:>18s} {:12.6g} {:12.6g} N/A "
435 "{:12.6g} {:12.6g}\n", speciesName(k),
436 x[k], molal[k], muss[k], acMolal[k]);
437 }
438 } else {
439 nMinor++;
440 xMinor += x[k];
441 }
442 }
443 } else {
444 fmt_append(b, " X Molalities\n");
445 fmt_append(b, " ------------- ------------\n");
446 for (size_t k = 0; k < m_kk; k++) {
447 if (x[k] > threshold) {
448 fmt_append(b, "{:>18s} {:12.6g} {:12.6g}\n",
449 speciesName(k), x[k], molal[k]);
450 } else {
451 nMinor++;
452 xMinor += x[k];
453 }
454 }
455 }
456 if (nMinor) {
457 fmt_append(b, " [{:+5d} minor] {:12.6g}\n", nMinor, xMinor);
458 }
459 } catch (CanteraError& err) {
460 return to_string(b) + err.what();
461 }
462 return to_string(b);
463}
464
465}
Header for intermediate ThermoPhase object for phases which employ molality based activity coefficien...
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:431
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition AnyMap.cpp:1477
A wrapper for a variable whose type is determined at runtime.
Definition AnyMap.h:87
const string & asString() const
Return the held value, if it is a string.
Definition AnyMap.cpp:782
map< string, T > asMap() const
Return the held AnyMap as a map where all of the values have the specified type.
Definition AnyMap.inl.h:162
bool is() const
Returns true if the held value is of the specified type.
Definition AnyMap.inl.h:68
Base class for exceptions thrown by Cantera classes.
const char * what() const override
Get a description of the error.
int activityConvention() const override
We set the convention to molality here.
void setState(const AnyMap &state) override
Set the state using an AnyMap containing any combination of properties supported by the thermodynamic...
void setMolalitiesByName(const Composition &xMap)
Set the molalities of a phase.
double m_Mnaught
This is the multiplication factor that goes inside log expressions involving the molalities of specie...
size_t m_indexCLM
Index of the phScale species.
size_t findCLMIndex() const
Returns the index of the Cl- species.
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
void getActivityConcentrations(double *c) const override
This method returns an array of generalized concentrations.
int pHScale() const
Reports the pH scale, which determines the scale for single-ion activity coefficients.
virtual void getMolalityActivityCoefficients(double *acMolality) const
Get the array of non-dimensional molality based activity coefficients at the current solution tempera...
void setpHScale(const int pHscaleType)
Set the pH scale, which determines the scale for single-ion activity coefficients.
double m_xmolSolventMIN
In any molality implementation, it makes sense to have a minimum solvent mole fraction requirement,...
vector< double > m_molalities
Current value of the molalities of the species in the phase.
virtual void applyphScale(double *acMolality) const
Apply the current phScale to a set of activity Coefficients or activities.
int m_pHScalingType
Scaling to be used for output of single-ion species activity coefficients.
void setMolalities(const double *const molal)
Set the molalities of the solutes in a phase.
void getMolalities(double *const molal) const
This function will return the molalities of the species.
void setMoleFSolventMin(double xmolSolventMIN)
Sets the minimum mole fraction in the molality formulation.
virtual double osmoticCoefficient() const
Calculate the osmotic coefficient.
MolalityVPSSTP()
Default Constructor.
double moleFSolventMin() const
Returns the minimum mole fraction in the molality formulation.
void setState_TPM(double t, double p, const double *const molalities)
Set the temperature (K), pressure (Pa), and molalities (gmol kg-1) of the solutes.
void getActivities(double *ac) const override
Get the array of non-dimensional activities (molality based for this class and classes that derive fr...
virtual void getUnscaledMolalityActivityCoefficients(double *acMolality) const
Get the array of unscaled non-dimensional molality based activity coefficients at the current solutio...
void calcMolalities() const
Calculates the molality of all species and stores the result internally.
double standardConcentration(size_t k=0) const override
Return the standard concentration for the kth species.
double m_weightSolvent
Molecular weight of the Solvent.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
void getActivityCoefficients(double *ac) const override
Get the array of non-dimensional activity coefficients at the current solution temperature,...
string report(bool show_thermo=true, double threshold=1e-14) const override
returns a summary of the state of the phase as a string
An error indicating that an unimplemented function has been called.
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:854
double temperature() const
Temperature (K).
Definition Phase.h:562
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:655
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition Phase.cpp:434
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 double density() const
Density (kg/m^3).
Definition Phase.h:587
double nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
Definition Phase.cpp:103
size_t nElements() const
Number of elements.
Definition Phase.cpp:30
const vector< string > & speciesNames() const
Return a const reference to the vector of species names.
Definition Phase.cpp:148
double molecularWeight(size_t k) const
Molecular weight of species k.
Definition Phase.cpp:383
string elementName(size_t m) const
Name of the element with index m.
Definition Phase.cpp:49
double charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition Phase.h:538
string name() const
Return the name of the phase.
Definition Phase.cpp:20
double electricPotential() const
Returns the electric potential of this phase (V).
virtual double cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
virtual double enthalpy_mole() const
Molar enthalpy. Units: J/kmol.
virtual void setState(const AnyMap &state)
Set the state using an AnyMap containing any combination of properties supported by the thermodynamic...
double gibbs_mass() const
Specific Gibbs function. Units: J/kg.
bool m_chargeNeutralityNecessary
Boolean indicating whether a charge neutrality condition is a necessity.
virtual double entropy_mole() const
Molar entropy. Units: J/kmol/K.
double cv_mass() const
Specific heat at constant volume. Units: J/kg/K.
double entropy_mass() const
Specific entropy. Units: J/kg/K.
virtual void getChemPotentials(double *mu) const
Get the species chemical potentials. Units: J/kmol.
double cp_mass() const
Specific heat at constant pressure. Units: J/kg/K.
double intEnergy_mass() const
Specific internal energy. Units: J/kg.
virtual double cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
virtual double intEnergy_mole() const
Molar internal energy. Units: J/kmol.
virtual double gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
double enthalpy_mass() const
Specific enthalpy. Units: J/kg.
double pressure() const override
Returns the current pressure of the phase.
void getStandardChemPotentials(double *mu) const override
Get the array of chemical potentials at unit activity for the species at their standard states at the...
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
void setState_TP(double T, double pres) override
Set the temperature and pressure at the same time.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
void fmt_append(fmt::memory_buffer &b, const std::string &tmpl, Args... args)
Versions 6.2.0 and 6.2.1 of fmtlib do not include this define before they include windows....
Definition fmt.h:29
Composition parseCompString(const string &ss, const vector< string > &names)
Parse a composition string into a map consisting of individual key:composition pairs.
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
const int PHSCALE_PITZER
Scale to be used for the output of single-ion activity coefficients is that used by Pitzer.
const double SmallNumber
smallest number to compare to zero.
Definition ct_defs.h:158
const int PHSCALE_NBS
Scale to be used for evaluation of single-ion activity coefficients is that used by the NBS standard ...
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
const int cAC_CONVENTION_MOLALITY
Standard state uses the molality convention.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...