Cantera 2.6.0
RedlichKwongMFTP.cpp
Go to the documentation of this file.
1//! @file RedlichKwongMFTP.cpp
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
10#include "cantera/base/ctml.h"
11
12#include <boost/math/tools/roots.hpp>
13
14#include <algorithm>
15
16using namespace std;
17namespace bmt = boost::math::tools;
18
19namespace Cantera
20{
21
22const doublereal RedlichKwongMFTP::omega_a = 4.27480233540E-01;
23const doublereal RedlichKwongMFTP::omega_b = 8.66403499650E-02;
24const doublereal RedlichKwongMFTP::omega_vc = 3.33333333333333E-01;
25
26RedlichKwongMFTP::RedlichKwongMFTP(const std::string& infile, const std::string& id_) :
27 m_formTempParam(0),
28 m_b_current(0.0),
29 m_a_current(0.0),
30 NSolns_(0),
31 dpdV_(0.0),
32 dpdT_(0.0)
33{
34 fill_n(Vroot_, 3, 0.0);
35 initThermoFile(infile, id_);
36}
37
38RedlichKwongMFTP::RedlichKwongMFTP(XML_Node& phaseRefRoot, const std::string& id_) :
39 m_formTempParam(0),
40 m_b_current(0.0),
41 m_a_current(0.0),
42 NSolns_(0),
43 dpdV_(0.0),
44 dpdT_(0.0)
45{
46 fill_n(Vroot_, 3, 0.0);
47 importPhase(phaseRefRoot, this);
48}
49
50void RedlichKwongMFTP::setSpeciesCoeffs(const std::string& species,
51 double a0, double a1, double b)
52{
53 size_t k = speciesIndex(species);
54 if (k == npos) {
55 throw CanteraError("RedlichKwongMFTP::setSpeciesCoeffs",
56 "Unknown species '{}'.", species);
57 }
58
59 if (a1 != 0.0) {
60 m_formTempParam = 1; // expression is temperature-dependent
61 }
62
63 size_t counter = k + m_kk * k;
64 a_coeff_vec(0, counter) = a0;
65 a_coeff_vec(1, counter) = a1;
66
67 // standard mixing rule for cross-species interaction term
68 for (size_t j = 0; j < m_kk; j++) {
69 if (k == j) {
70 continue;
71 }
72
73 // a_coeff_vec(0) is initialized to NaN to mark uninitialized species
74 if (isnan(a_coeff_vec(0, j + m_kk * j))) {
75 // The diagonal element of the jth species has not yet been defined.
76 continue;
77 } else if (isnan(a_coeff_vec(0, j + m_kk * k))) {
78 // Only use the mixing rules if the off-diagonal element has not already been defined by a
79 // user-specified crossFluidParameters entry:
80 double a0kj = sqrt(a_coeff_vec(0, j + m_kk * j) * a0);
81 double a1kj = sqrt(a_coeff_vec(1, j + m_kk * j) * a1);
82 a_coeff_vec(0, j + m_kk * k) = a0kj;
83 a_coeff_vec(1, j + m_kk * k) = a1kj;
84 a_coeff_vec(0, k + m_kk * j) = a0kj;
85 a_coeff_vec(1, k + m_kk * j) = a1kj;
86 }
87 }
88 a_coeff_vec.getRow(0, a_vec_Curr_.data());
89 b_vec_Curr_[k] = b;
90}
91
92void RedlichKwongMFTP::setBinaryCoeffs(const std::string& species_i,
93 const std::string& species_j, double a0, double a1)
94{
95 size_t ki = speciesIndex(species_i);
96 if (ki == npos) {
97 throw CanteraError("RedlichKwongMFTP::setBinaryCoeffs",
98 "Unknown species '{}'.", species_i);
99 }
100 size_t kj = speciesIndex(species_j);
101 if (kj == npos) {
102 throw CanteraError("RedlichKwongMFTP::setBinaryCoeffs",
103 "Unknown species '{}'.", species_j);
104 }
105
106 if (a1 != 0.0) {
107 m_formTempParam = 1; // expression is temperature-dependent
108 }
109
110 m_binaryParameters[species_i][species_j] = {a0, a1};
111 m_binaryParameters[species_j][species_i] = {a0, a1};
112 size_t counter1 = ki + m_kk * kj;
113 size_t counter2 = kj + m_kk * ki;
114 a_coeff_vec(0, counter1) = a_coeff_vec(0, counter2) = a0;
115 a_coeff_vec(1, counter1) = a_coeff_vec(1, counter2) = a1;
116 a_vec_Curr_[counter1] = a_vec_Curr_[counter2] = a0;
117}
118
119// ------------Molar Thermodynamic Properties -------------------------
120
122{
124 doublereal TKelvin = temperature();
125 doublereal sqt = sqrt(TKelvin);
126 doublereal mv = molarVolume();
127 doublereal vpb = mv + m_b_current;
129 doublereal cpref = GasConstant * mean_X(m_cp0_R);
130 doublereal dadt = da_dt();
131 doublereal fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
132 doublereal dHdT_V = (cpref + mv * dpdT_ - GasConstant - 1.0 / (2.0 * m_b_current * TKelvin * sqt) * log(vpb/mv) * fac
133 +1.0/(m_b_current * sqt) * log(vpb/mv) * (-0.5 * dadt));
134 return dHdT_V - (mv + TKelvin * dpdT_ / dpdV_) * dpdT_;
135}
136
138{
140 doublereal TKelvin = temperature();
141 doublereal sqt = sqrt(TKelvin);
142 doublereal mv = molarVolume();
143 doublereal vpb = mv + m_b_current;
144 doublereal cvref = GasConstant * (mean_X(m_cp0_R) - 1.0);
145 doublereal dadt = da_dt();
146 doublereal fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
147 return (cvref - 1.0/(2.0 * m_b_current * TKelvin * sqt) * log(vpb/mv)*fac
148 +1.0/(m_b_current * sqt) * log(vpb/mv)*(-0.5*dadt));
149}
150
152{
154
155 // Get a copy of the private variables stored in the State object
156 doublereal T = temperature();
157 double molarV = meanMolecularWeight() / density();
158 double pp = GasConstant * T/(molarV - m_b_current) - m_a_current/(sqrt(T) * molarV * (molarV + m_b_current));
159 return pp;
160}
161
163{
165 return 1.0 / m_tmpV[k];
166}
167
169{
170 doublereal mv = molarVolume();
171 doublereal sqt = sqrt(temperature());
172 doublereal vpb = mv + m_b_current;
173 doublereal vmb = mv - m_b_current;
174
175 for (size_t k = 0; k < m_kk; k++) {
176 m_pp[k] = 0.0;
177 for (size_t i = 0; i < m_kk; i++) {
178 size_t counter = k + m_kk*i;
179 m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
180 }
181 }
182 doublereal pres = pressure();
183
184 for (size_t k = 0; k < m_kk; k++) {
185 ac[k] = (- RT() * log(pres * mv / RT())
186 + RT() * log(mv / vmb)
187 + RT() * b_vec_Curr_[k] / vmb
188 - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
189 + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
190 - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
191 );
192 }
193 for (size_t k = 0; k < m_kk; k++) {
194 ac[k] = exp(ac[k]/RT());
195 }
196}
197
198// ---- Partial Molar Properties of the Solution -----------------
199
200void RedlichKwongMFTP::getChemPotentials_RT(doublereal* muRT) const
201{
202 getChemPotentials(muRT);
203 for (size_t k = 0; k < m_kk; k++) {
204 muRT[k] *= 1.0 / RT();
205 }
206}
207
208void RedlichKwongMFTP::getChemPotentials(doublereal* mu) const
209{
210 getGibbs_ref(mu);
211 for (size_t k = 0; k < m_kk; k++) {
212 double xx = std::max(SmallNumber, moleFraction(k));
213 mu[k] += RT()*(log(xx));
214 }
215
216 doublereal mv = molarVolume();
217 doublereal sqt = sqrt(temperature());
218 doublereal vpb = mv + m_b_current;
219 doublereal vmb = mv - m_b_current;
220
221 for (size_t k = 0; k < m_kk; k++) {
222 m_pp[k] = 0.0;
223 for (size_t i = 0; i < m_kk; i++) {
224 size_t counter = k + m_kk*i;
225 m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
226 }
227 }
228 doublereal pres = pressure();
229 doublereal refP = refPressure();
230
231 for (size_t k = 0; k < m_kk; k++) {
232 mu[k] += (RT() * log(pres/refP) - RT() * log(pres * mv / RT())
233 + RT() * log(mv / vmb)
234 + RT() * b_vec_Curr_[k] / vmb
235 - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
236 + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
237 - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
238 );
239 }
240}
241
243{
244 // First we get the reference state contributions
245 getEnthalpy_RT_ref(hbar);
246 scale(hbar, hbar+m_kk, hbar, RT());
247
248 // We calculate dpdni_
249 doublereal TKelvin = temperature();
250 doublereal mv = molarVolume();
251 doublereal sqt = sqrt(TKelvin);
252 doublereal vpb = mv + m_b_current;
253 doublereal vmb = mv - m_b_current;
254 for (size_t k = 0; k < m_kk; k++) {
255 m_pp[k] = 0.0;
256 for (size_t i = 0; i < m_kk; i++) {
257 size_t counter = k + m_kk*i;
258 m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
259 }
260 }
261 for (size_t k = 0; k < m_kk; k++) {
262 dpdni_[k] = RT()/vmb + RT() * b_vec_Curr_[k] / (vmb * vmb) - 2.0 * m_pp[k] / (sqt * mv * vpb)
263 + m_a_current * b_vec_Curr_[k]/(sqt * mv * vpb * vpb);
264 }
265 doublereal dadt = da_dt();
266 doublereal fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
267
268 for (size_t k = 0; k < m_kk; k++) {
269 m_tmpV[k] = 0.0;
270 for (size_t i = 0; i < m_kk; i++) {
271 size_t counter = k + m_kk*i;
272 m_tmpV[k] += 2.0 * moleFractions_[i] * TKelvin * a_coeff_vec(1,counter) - 3.0 * moleFractions_[i] * a_vec_Curr_[counter];
273 }
274 }
275
277 doublereal fac2 = mv + TKelvin * dpdT_ / dpdV_;
278 for (size_t k = 0; k < m_kk; k++) {
279 double hE_v = (mv * dpdni_[k] - RT() - b_vec_Curr_[k]/ (m_b_current * m_b_current * sqt) * log(vpb/mv)*fac
280 + 1.0 / (m_b_current * sqt) * log(vpb/mv) * m_tmpV[k]
281 + b_vec_Curr_[k] / vpb / (m_b_current * sqt) * fac);
282 hbar[k] = hbar[k] + hE_v;
283 hbar[k] -= fac2 * dpdni_[k];
284 }
285}
286
288{
289 getEntropy_R_ref(sbar);
290 scale(sbar, sbar+m_kk, sbar, GasConstant);
291 doublereal TKelvin = temperature();
292 doublereal sqt = sqrt(TKelvin);
293 doublereal mv = molarVolume();
294 doublereal refP = refPressure();
295
296 for (size_t k = 0; k < m_kk; k++) {
297 doublereal xx = std::max(SmallNumber, moleFraction(k));
298 sbar[k] += GasConstant * (- log(xx));
299 }
300 for (size_t k = 0; k < m_kk; k++) {
301 m_pp[k] = 0.0;
302 for (size_t i = 0; i < m_kk; i++) {
303 size_t counter = k + m_kk*i;
304 m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
305 }
306 }
307 for (size_t k = 0; k < m_kk; k++) {
308 m_tmpV[k] = 0.0;
309 for (size_t i = 0; i < m_kk; i++) {
310 size_t counter = k + m_kk*i;
311 m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
312 }
313 }
314
315 doublereal dadt = da_dt();
316 doublereal fac = dadt - m_a_current / (2.0 * TKelvin);
317 doublereal vmb = mv - m_b_current;
318 doublereal vpb = mv + m_b_current;
319 for (size_t k = 0; k < m_kk; k++) {
320 sbar[k] -=(GasConstant * log(GasConstant * TKelvin / (refP * mv))
322 + GasConstant * log(mv/vmb)
323 + GasConstant * b_vec_Curr_[k]/vmb
324 + m_pp[k]/(m_b_current * TKelvin * sqt) * log(vpb/mv)
325 - 2.0 * m_tmpV[k]/(m_b_current * sqt) * log(vpb/mv)
326 + b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) * fac
327 - 1.0 / (m_b_current * sqt) * b_vec_Curr_[k] / vpb * fac
328 );
329 }
330
332 getPartialMolarVolumes(m_partialMolarVolumes.data());
333 for (size_t k = 0; k < m_kk; k++) {
334 sbar[k] -= -m_partialMolarVolumes[k] * dpdT_;
335 }
336}
337
339{
340 // u_k = h_k - P * v_k
341 getPartialMolarVolumes(m_partialMolarVolumes.data());
343 double p = pressure();
344 for (size_t k = 0; k < nSpecies(); k++) {
345 ubar[k] -= p * m_partialMolarVolumes[k];
346 }
347}
348
350{
351 for (size_t k = 0; k < m_kk; k++) {
352 m_pp[k] = 0.0;
353 for (size_t i = 0; i < m_kk; i++) {
354 size_t counter = k + m_kk*i;
355 m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
356 }
357 }
358 for (size_t k = 0; k < m_kk; k++) {
359 m_tmpV[k] = 0.0;
360 for (size_t i = 0; i < m_kk; i++) {
361 size_t counter = k + m_kk*i;
362 m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
363 }
364 }
365
366 doublereal sqt = sqrt(temperature());
367 doublereal mv = molarVolume();
368 doublereal vmb = mv - m_b_current;
369 doublereal vpb = mv + m_b_current;
370 for (size_t k = 0; k < m_kk; k++) {
371 doublereal num = (RT() + RT() * m_b_current/ vmb + RT() * b_vec_Curr_[k] / vmb
372 + RT() * m_b_current * b_vec_Curr_[k] /(vmb * vmb)
373 - 2.0 * m_pp[k] / (sqt * vpb)
374 + m_a_current * b_vec_Curr_[k] / (sqt * vpb * vpb)
375 );
376 doublereal denom = (pressure() + RT() * m_b_current/(vmb * vmb) - m_a_current / (sqt * vpb * vpb)
377 );
378 vbar[k] = num / denom;
379 }
380}
381
382bool RedlichKwongMFTP::addSpecies(shared_ptr<Species> spec)
383{
384 bool added = MixtureFugacityTP::addSpecies(spec);
385 if (added) {
386 a_vec_Curr_.resize(m_kk * m_kk, 0.0);
387
388 // Initialize a_vec and b_vec to NaN, to screen for species with
389 // pureFluidParameters which are undefined in the input file:
390 b_vec_Curr_.push_back(NAN);
391 a_coeff_vec.resize(2, m_kk * m_kk, NAN);
392
393 m_pp.push_back(0.0);
394 m_coeffSource.push_back(CoeffSource::EoS);
395 m_partialMolarVolumes.push_back(0.0);
396 dpdni_.push_back(0.0);
397 }
398 return added;
399}
400
401void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
402{
403 if (phaseNode.hasChild("thermo")) {
404 XML_Node& thermoNode = phaseNode.child("thermo");
405 std::string model = thermoNode["model"];
406 if (model != "RedlichKwong" && model != "RedlichKwongMFTP") {
407 throw CanteraError("RedlichKwongMFTP::initThermoXML",
408 "Unknown thermo model : " + model);
409 }
410
411 // Reset any coefficients which may have been set using values from
412 // 'critical-properties.yaml' as part of non-XML initialization, so that
413 // off-diagonal elements can be correctly initialized
414 a_coeff_vec.data().assign(a_coeff_vec.data().size(), NAN);
415
416 // Go get all of the coefficients and factors in the
417 // activityCoefficients XML block
418 if (thermoNode.hasChild("activityCoefficients")) {
419 XML_Node& acNode = thermoNode.child("activityCoefficients");
420
421 // Loop through the children and read out fluid parameters. Process
422 // all the pureFluidParameters, first:
423 // Loop back through the "activityCoefficients" children and process the
424 // crossFluidParameters in the XML tree:
425 for (size_t i = 0; i < acNode.nChildren(); i++) {
426 XML_Node& xmlACChild = acNode.child(i);
427 if (caseInsensitiveEquals(xmlACChild.name(), "purefluidparameters")) {
428 readXMLPureFluid(xmlACChild);
429 } else if (caseInsensitiveEquals(xmlACChild.name(), "crossfluidparameters")) {
430 readXMLCrossFluid(xmlACChild);
431 }
432 }
433 }
434 // If any species exist which have undefined pureFluidParameters,
435 // search the database in 'critical-properties.yaml' to find critical
436 // temperature and pressure to calculate a and b.
437
438 // Loop through all species in the CTI file
439 size_t iSpecies = 0;
440
441 for (size_t i = 0; i < m_kk; i++) {
442 string iName = speciesName(i);
443
444 // Get the index of the species
445 iSpecies = speciesIndex(iName);
446
447 // Check if a and b are already populated (only the diagonal elements of a).
448 size_t counter = iSpecies + m_kk * iSpecies;
449
450 // If not, then search the database:
451 if (isnan(a_coeff_vec(0, counter))) {
452
453 vector<double> coeffArray;
454
455 // Search the database for the species name and calculate
456 // coefficients a and b, from critical properties:
457 // coeffArray[0] = a0, coeffArray[1] = b;
458 coeffArray = getCoeff(iName);
459
460 // Check if species was found in the database of critical properties,
461 // and assign the results
462 if (!isnan(coeffArray[0])) {
463 //Assuming no temperature dependence (i,e a1 = 0)
464 setSpeciesCoeffs(iName, coeffArray[0], 0.0, coeffArray[1]);
465 m_coeffSource[i] = CoeffSource::EoS;
466 }
467 }
468 }
469 }
470
472}
473
475{
476 // Contents of 'critical-properties.yaml', loaded later if needed
477 AnyMap critPropsDb;
478 std::unordered_map<std::string, AnyMap*> dbSpecies;
479
480 for (auto& item : m_species) {
481 auto& data = item.second->input;
482 size_t k = speciesIndex(item.first);
483 if (!isnan(a_coeff_vec(0, k + m_kk * k))) {
484 continue;
485 }
486 bool foundCoeffs = false;
487
488 if (data.hasKey("equation-of-state") &&
489 data["equation-of-state"].hasMapWhere("model", "Redlich-Kwong"))
490 {
491 // Read a and b coefficients from species 'input' information (that is, as
492 // specified in a YAML input file) specific to the Redlich-Kwong model
493 auto eos = data["equation-of-state"].getMapWhere(
494 "model", "Redlich-Kwong");
495
496 if (eos.hasKey("a") && eos.hasKey("b")) {
497 double a0 = 0, a1 = 0;
498 if (eos["a"].isScalar()) {
499 a0 = eos.convert("a", "Pa*m^6/kmol^2*K^0.5");
500 } else {
501 auto avec = eos["a"].asVector<AnyValue>(2);
502 a0 = eos.units().convert(avec[0], "Pa*m^6/kmol^2*K^0.5");
503 a1 = eos.units().convert(avec[1], "Pa*m^6/kmol^2/K^0.5");
504 }
505 double b = eos.convert("b", "m^3/kmol");
506 foundCoeffs = true;
507 setSpeciesCoeffs(item.first, a0, a1, b);
508 m_coeffSource[k] = CoeffSource::EoS;
509 }
510
511 if (eos.hasKey("binary-a")) {
512 AnyMap& binary_a = eos["binary-a"].as<AnyMap>();
513 const UnitSystem& units = binary_a.units();
514 for (auto& item2 : binary_a) {
515 double a0 = 0, a1 = 0;
516 if (item2.second.isScalar()) {
517 a0 = units.convert(item2.second, "Pa*m^6/kmol^2*K^0.5");
518 } else {
519 auto avec = item2.second.asVector<AnyValue>(2);
520 a0 = units.convert(avec[0], "Pa*m^6/kmol^2*K^0.5");
521 a1 = units.convert(avec[1], "Pa*m^6/kmol^2/K^0.5");
522 }
523 setBinaryCoeffs(item.first, item2.first, a0, a1);
524 }
525 }
526
527 if (foundCoeffs) {
528 continue;
529 }
530 }
531
532 // Coefficients have not been populated from model-specific input
533 double Tc = NAN, Pc = NAN;
534 if (data.hasKey("critical-parameters")) {
535 // Use critical state information stored in the species entry to
536 // calculate a and b
537 auto& critProps = data["critical-parameters"].as<AnyMap>();
538 Tc = critProps.convert("critical-temperature", "K");
539 Pc = critProps.convert("critical-pressure", "Pa");
540 m_coeffSource[k] = CoeffSource::CritProps;
541 } else {
542 // Search 'crit-properties.yaml' to find Tc and Pc. Load data if needed
543 if (critPropsDb.empty()) {
544 critPropsDb = AnyMap::fromYamlFile("critical-properties.yaml");
545 dbSpecies = critPropsDb["species"].asMap("name");
546 }
547
548 // All names in critical-properties.yaml are upper case
549 auto ucName = boost::algorithm::to_upper_copy(item.first);
550 if (dbSpecies.count(ucName)) {
551 auto& spec = *dbSpecies.at(ucName);
552 auto& critProps = spec["critical-parameters"].as<AnyMap>();
553 Tc = critProps.convert("critical-temperature", "K");
554 Pc = critProps.convert("critical-pressure", "Pa");
555 m_coeffSource[k] = CoeffSource::Database;
556 }
557 }
558
559 // Check if critical properties were found in either location
560 if (!isnan(Tc)) {
561 // Assuming no temperature dependence (that is, a1 = 0)
562 double a = omega_a * pow(GasConstant, 2) * pow(Tc, 2.5) / Pc;
563 double b = omega_b * GasConstant * Tc / Pc;
564 setSpeciesCoeffs(item.first, a, 0.0, b);
565 }
566 // @todo: After XML is removed, this can throw an InputFileError if neither
567 // R-K parameters or critical properties were found.
568 }
569}
570
571void RedlichKwongMFTP::getSpeciesParameters(const std::string& name,
572 AnyMap& speciesNode) const
573{
575 size_t k = speciesIndex(name);
577 if (m_coeffSource[k] == CoeffSource::EoS) {
578 auto& eosNode = speciesNode["equation-of-state"].getMapWhere(
579 "model", "Redlich-Kwong", true);
580
581 size_t counter = k + m_kk * k;
582 if (a_coeff_vec(1, counter) != 0.0) {
583 vector<AnyValue> coeffs(2);
584 coeffs[0].setQuantity(a_coeff_vec(0, counter), "Pa*m^6/kmol^2*K^0.5");
585 coeffs[1].setQuantity(a_coeff_vec(1, counter), "Pa*m^6/kmol^2/K^0.5");
586 eosNode["a"] = std::move(coeffs);
587 } else {
588 eosNode["a"].setQuantity(a_coeff_vec(0, counter),
589 "Pa*m^6/kmol^2*K^0.5");
590 }
591 eosNode["b"].setQuantity(b_vec_Curr_[k], "m^3/kmol");
592 } else if (m_coeffSource[k] == CoeffSource::CritProps) {
593 auto& critProps = speciesNode["critical-parameters"];
594 double a = a_coeff_vec(0, k + m_kk * k);
595 double b = b_vec_Curr_[k];
596 double Tc = pow(a * omega_b / (b * omega_a * GasConstant), 2.0/3.0);
597 double Pc = omega_b * GasConstant * Tc / b;
598 critProps["critical-temperature"].setQuantity(Tc, "K");
599 critProps["critical-pressure"].setQuantity(Pc, "Pa");
600 }
601
602 if (m_binaryParameters.count(name)) {
603 auto& eosNode = speciesNode["equation-of-state"].getMapWhere(
604 "model", "Redlich-Kwong", true);
605 AnyMap bin_a;
606 for (const auto& item : m_binaryParameters.at(name)) {
607 if (item.second.second == 0) {
608 bin_a[item.first].setQuantity(item.second.first, "Pa*m^6/kmol^2*K^0.5");
609 } else {
610 vector<AnyValue> coeffs(2);
611 coeffs[0].setQuantity(item.second.first, "Pa*m^6/kmol^2*K^0.5");
612 coeffs[1].setQuantity(item.second.second, "Pa*m^6/kmol^2/K^0.5");
613 bin_a[item.first] = std::move(coeffs);
614 }
615 }
616 eosNode["binary-a"] = std::move(bin_a);
617 }
618}
619
620vector<double> RedlichKwongMFTP::getCoeff(const std::string& iName)
621{
622 warn_deprecated("RedlichKwongMFTP::getCoeff", "To be removed after Cantera 2.6. "
623 "Use of critical-properties.yaml is integrated into initThermo() "
624 "for YAML input files.");
625 vector_fp spCoeff{NAN, NAN};
626 AnyMap data = AnyMap::fromYamlFile("critical-properties.yaml");
627 const auto& species = data["species"].asMap("name");
628
629 // All names in critical-properties.yaml are upper case
630 auto ucName = boost::algorithm::to_upper_copy(iName);
631 if (species.count(ucName)) {
632 auto& critProps = species.at(ucName)->at("critical-parameters").as<AnyMap>();
633 double Tc = critProps.convert("critical-temperature", "K");
634 double Pc = critProps.convert("critical-pressure", "Pa");
635
636 // calculate pure fluid parameters a_k and b_k based on T_c and P_c assuming
637 // no temperature dependence
638 spCoeff[0] = omega_a * pow(GasConstant, 2) * pow(Tc, 2.5) / Pc; //coeff a
639 spCoeff[1] = omega_b * GasConstant * Tc / Pc; // coeff b
640 }
641 return spCoeff;
642}
643
645{
646 string xname = pureFluidParam.name();
647 if (xname != "pureFluidParameters") {
648 throw CanteraError("RedlichKwongMFTP::readXMLPureFluid",
649 "Incorrect name for processing this routine: " + xname);
650 }
651
652 double a0 = 0.0;
653 double a1 = 0.0;
654 double b = 0.0;
655 for (size_t iChild = 0; iChild < pureFluidParam.nChildren(); iChild++) {
656 XML_Node& xmlChild = pureFluidParam.child(iChild);
657 string nodeName = toLowerCopy(xmlChild.name());
658
659 if (nodeName == "a_coeff") {
660 vector_fp vParams;
661 string iModel = toLowerCopy(xmlChild.attrib("model"));
662 getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff");
663
664 if (iModel == "constant" && vParams.size() == 1) {
665 a0 = vParams[0];
666 a1 = 0;
667 } else if (iModel == "linear_a" && vParams.size() == 2) {
668 a0 = vParams[0];
669 a1 = vParams[1];
670 } else {
671 throw CanteraError("RedlichKwongMFTP::readXMLPureFluid",
672 "unknown model or incorrect number of parameters");
673 }
674
675 } else if (nodeName == "b_coeff") {
676 b = getFloatCurrent(xmlChild, "toSI");
677 }
678 }
679 setSpeciesCoeffs(pureFluidParam.attrib("species"), a0, a1, b);
680 m_coeffSource[speciesIndex(pureFluidParam.attrib("species"))] = CoeffSource::EoS;
681}
682
684{
685 string xname = CrossFluidParam.name();
686 if (xname != "crossFluidParameters") {
687 throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid",
688 "Incorrect name for processing this routine: " + xname);
689 }
690
691 string iName = CrossFluidParam.attrib("species1");
692 string jName = CrossFluidParam.attrib("species2");
693
694 size_t num = CrossFluidParam.nChildren();
695 for (size_t iChild = 0; iChild < num; iChild++) {
696 XML_Node& xmlChild = CrossFluidParam.child(iChild);
697 string nodeName = toLowerCopy(xmlChild.name());
698
699 if (nodeName == "a_coeff") {
700 vector_fp vParams;
701 getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff");
702 string iModel = toLowerCopy(xmlChild.attrib("model"));
703 if (iModel == "constant" && vParams.size() == 1) {
704 setBinaryCoeffs(iName, jName, vParams[0], 0.0);
705 } else if (iModel == "linear_a") {
706 setBinaryCoeffs(iName, jName, vParams[0], vParams[1]);
707 } else {
708 throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid",
709 "unknown model ({}) or wrong number of parameters ({})",
710 iModel, vParams.size());
711 }
712 }
713 }
714}
715
717{
719 std::string model = thermoNode["model"];
720}
721
722doublereal RedlichKwongMFTP::sresid() const
723{
724 // note this agrees with tpx
725 doublereal rho = density();
726 doublereal mmw = meanMolecularWeight();
727 doublereal molarV = mmw / rho;
728 double hh = m_b_current / molarV;
729 doublereal zz = z();
730 doublereal dadt = da_dt();
731 doublereal T = temperature();
732 doublereal sqT = sqrt(T);
733 doublereal fac = dadt - m_a_current / (2.0 * T);
734 double sresid_mol_R = log(zz*(1.0 - hh)) + log(1.0 + hh) * fac / (sqT * GasConstant * m_b_current);
735 return GasConstant * sresid_mol_R;
736}
737
738doublereal RedlichKwongMFTP::hresid() const
739{
740 // note this agrees with tpx
741 doublereal rho = density();
742 doublereal mmw = meanMolecularWeight();
743 doublereal molarV = mmw / rho;
744 double hh = m_b_current / molarV;
745 doublereal zz = z();
746 doublereal dadt = da_dt();
747 doublereal T = temperature();
748 doublereal sqT = sqrt(T);
749 doublereal fac = T * dadt - 3.0 *m_a_current / (2.0);
750 return GasConstant * T * (zz - 1.0) + fac * log(1.0 + hh) / (sqT * m_b_current);
751}
752
753doublereal RedlichKwongMFTP::liquidVolEst(doublereal TKelvin, doublereal& presGuess) const
754{
755 double v = m_b_current * 1.1;
756 double atmp;
757 double btmp;
758 calculateAB(TKelvin, atmp, btmp);
759 doublereal pres = std::max(psatEst(TKelvin), presGuess);
760 double Vroot[3];
761 bool foundLiq = false;
762 int m = 0;
763 while (m < 100 && !foundLiq) {
764 int nsol = solveCubic(TKelvin, pres, atmp, btmp, Vroot);
765 if (nsol == 1 || nsol == 2) {
766 double pc = critPressure();
767 if (pres > pc) {
768 foundLiq = true;
769 }
770 pres *= 1.04;
771 } else {
772 foundLiq = true;
773 }
774 }
775
776 if (foundLiq) {
777 v = Vroot[0];
778 presGuess = pres;
779 } else {
780 v = -1.0;
781 }
782 return v;
783}
784
785doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa, int phaseRequested, doublereal rhoguess)
786{
787 // It's necessary to set the temperature so that m_a_current is set correctly.
788 setTemperature(TKelvin);
789 double tcrit = critTemperature();
790 double mmw = meanMolecularWeight();
791 if (rhoguess == -1.0) {
792 if (phaseRequested >= FLUID_LIQUID_0) {
793 double lqvol = liquidVolEst(TKelvin, presPa);
794 rhoguess = mmw / lqvol;
795 }
796 } else {
797 // Assume the Gas phase initial guess, if nothing is specified to the routine
798 rhoguess = presPa * mmw / (GasConstant * TKelvin);
799 }
800
801
802 doublereal volguess = mmw / rhoguess;
803 NSolns_ = solveCubic(TKelvin, presPa, m_a_current, m_b_current, Vroot_);
804
805 doublereal molarVolLast = Vroot_[0];
806 if (NSolns_ >= 2) {
807 if (phaseRequested >= FLUID_LIQUID_0) {
808 molarVolLast = Vroot_[0];
809 } else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
810 molarVolLast = Vroot_[2];
811 } else {
812 if (volguess > Vroot_[1]) {
813 molarVolLast = Vroot_[2];
814 } else {
815 molarVolLast = Vroot_[0];
816 }
817 }
818 } else if (NSolns_ == 1) {
819 if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT || phaseRequested == FLUID_UNDEFINED) {
820 molarVolLast = Vroot_[0];
821 } else {
822 return -2.0;
823 }
824 } else if (NSolns_ == -1) {
825 if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) {
826 molarVolLast = Vroot_[0];
827 } else if (TKelvin > tcrit) {
828 molarVolLast = Vroot_[0];
829 } else {
830 return -2.0;
831 }
832 } else {
833 molarVolLast = Vroot_[0];
834 return -1.0;
835 }
836 return mmw / molarVolLast;
837}
838
840{
841 double Vroot[3];
842 double T = temperature();
843 int nsol = solveCubic(T, pressure(), m_a_current, m_b_current, Vroot);
844 if (nsol != 3) {
845 return critDensity();
846 }
847
848 auto resid = [this, T](double v) {
849 double pp;
850 return dpdVCalc(T, v, pp);
851 };
852
853 boost::uintmax_t maxiter = 100;
854 std::pair<double, double> vv = bmt::toms748_solve(
855 resid, Vroot[0], Vroot[1], bmt::eps_tolerance<double>(48), maxiter);
856
857 doublereal mmw = meanMolecularWeight();
858 return mmw / (0.5 * (vv.first + vv.second));
859}
860
862{
863 double Vroot[3];
864 double T = temperature();
865 int nsol = solveCubic(T, pressure(), m_a_current, m_b_current, Vroot);
866 if (nsol != 3) {
867 return critDensity();
868 }
869
870 auto resid = [this, T](double v) {
871 double pp;
872 return dpdVCalc(T, v, pp);
873 };
874
875 boost::uintmax_t maxiter = 100;
876 std::pair<double, double> vv = bmt::toms748_solve(
877 resid, Vroot[1], Vroot[2], bmt::eps_tolerance<double>(48), maxiter);
878
879 doublereal mmw = meanMolecularWeight();
880 return mmw / (0.5 * (vv.first + vv.second));
881}
882
883doublereal RedlichKwongMFTP::dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal& presCalc) const
884{
885 doublereal sqt = sqrt(TKelvin);
886 presCalc = GasConstant * TKelvin / (molarVol - m_b_current)
887 - m_a_current / (sqt * molarVol * (molarVol + m_b_current));
888
889 doublereal vpb = molarVol + m_b_current;
890 doublereal vmb = molarVol - m_b_current;
891 doublereal dpdv = (- GasConstant * TKelvin / (vmb * vmb)
892 + m_a_current * (2 * molarVol + m_b_current) / (sqt * molarVol * molarVol * vpb * vpb));
893 return dpdv;
894}
895
897{
898 doublereal TKelvin = temperature();
899 doublereal mv = molarVolume();
900 doublereal pres;
901
902 dpdV_ = dpdVCalc(TKelvin, mv, pres);
903 doublereal sqt = sqrt(TKelvin);
904 doublereal vpb = mv + m_b_current;
905 doublereal vmb = mv - m_b_current;
906 doublereal dadt = da_dt();
907 doublereal fac = dadt - m_a_current/(2.0 * TKelvin);
908 dpdT_ = (GasConstant / vmb - fac / (sqt * mv * vpb));
909}
910
912{
913 double temp = temperature();
914 if (m_formTempParam == 1) {
915 for (size_t i = 0; i < m_kk; i++) {
916 for (size_t j = 0; j < m_kk; j++) {
917 size_t counter = i * m_kk + j;
918 a_vec_Curr_[counter] = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
919 }
920 }
921 }
922
923 m_b_current = 0.0;
924 m_a_current = 0.0;
925 for (size_t i = 0; i < m_kk; i++) {
926 m_b_current += moleFractions_[i] * b_vec_Curr_[i];
927 for (size_t j = 0; j < m_kk; j++) {
928 m_a_current += a_vec_Curr_[i * m_kk + j] * moleFractions_[i] * moleFractions_[j];
929 }
930 }
931 if (isnan(m_b_current)) {
932 // One or more species do not have specified coefficients.
933 fmt::memory_buffer b;
934 for (size_t k = 0; k < m_kk; k++) {
935 if (isnan(b_vec_Curr_[k])) {
936 if (b.size() > 0) {
937 fmt_append(b, ", {}", speciesName(k));
938 } else {
939 fmt_append(b, "{}", speciesName(k));
940 }
941 }
942 }
943 throw CanteraError("RedlichKwongMFTP::updateMixingExpressions",
944 "Missing Redlich-Kwong coefficients for species: {}", to_string(b));
945 }
946}
947
948void RedlichKwongMFTP::calculateAB(doublereal temp, doublereal& aCalc, doublereal& bCalc) const
949{
950 bCalc = 0.0;
951 aCalc = 0.0;
952 if (m_formTempParam == 1) {
953 for (size_t i = 0; i < m_kk; i++) {
954 bCalc += moleFractions_[i] * b_vec_Curr_[i];
955 for (size_t j = 0; j < m_kk; j++) {
956 size_t counter = i * m_kk + j;
957 doublereal a_vec_Curr = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
958 aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
959 }
960 }
961 } else {
962 for (size_t i = 0; i < m_kk; i++) {
963 bCalc += moleFractions_[i] * b_vec_Curr_[i];
964 for (size_t j = 0; j < m_kk; j++) {
965 size_t counter = i * m_kk + j;
966 doublereal a_vec_Curr = a_coeff_vec(0,counter);
967 aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
968 }
969 }
970 }
971}
972
973doublereal RedlichKwongMFTP::da_dt() const
974{
975 doublereal dadT = 0.0;
976 if (m_formTempParam == 1) {
977 for (size_t i = 0; i < m_kk; i++) {
978 for (size_t j = 0; j < m_kk; j++) {
979 size_t counter = i * m_kk + j;
980 dadT+= a_coeff_vec(1,counter) * moleFractions_[i] * moleFractions_[j];
981 }
982 }
983 }
984 return dadT;
985}
986
987void RedlichKwongMFTP::calcCriticalConditions(doublereal& pc, doublereal& tc, doublereal& vc) const
988{
989 double a0 = 0.0;
990 double aT = 0.0;
991 for (size_t i = 0; i < m_kk; i++) {
992 for (size_t j = 0; j <m_kk; j++) {
993 size_t counter = i + m_kk * j;
994 a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
995 aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
996 }
997 }
998 double a = m_a_current;
999 double b = m_b_current;
1000 if (m_formTempParam != 0) {
1001 a = a0;
1002 }
1003 if (b <= 0.0) {
1004 tc = 1000000.;
1005 pc = 1.0E13;
1006 vc = omega_vc * GasConstant * tc / pc;
1007 return;
1008 }
1009 if (a <= 0.0) {
1010 tc = 0.0;
1011 pc = 0.0;
1012 vc = 2.0 * b;
1013 return;
1014 }
1015 double tmp = a * omega_b / (b * omega_a * GasConstant);
1016 double pp = 2./3.;
1017 doublereal sqrttc, f, dfdt, deltatc;
1018
1019 if (m_formTempParam == 0) {
1020 tc = pow(tmp, pp);
1021 } else {
1022 tc = pow(tmp, pp);
1023 for (int j = 0; j < 10; j++) {
1024 sqrttc = sqrt(tc);
1025 f = omega_a * b * GasConstant * tc * sqrttc / omega_b - aT * tc - a0;
1026 dfdt = 1.5 * omega_a * b * GasConstant * sqrttc / omega_b - aT;
1027 deltatc = - f / dfdt;
1028 tc += deltatc;
1029 }
1030 if (deltatc > 0.1) {
1031 throw CanteraError("RedlichKwongMFTP::calcCriticalConditions",
1032 "didn't converge");
1033 }
1034 }
1035
1036 pc = omega_b * GasConstant * tc / b;
1037 vc = omega_vc * GasConstant * tc / pc;
1038}
1039
1040int RedlichKwongMFTP::solveCubic(double T, double pres, double a, double b, double Vroot[3]) const
1041{
1042
1043 // Derive the coefficients of the cubic polynomial to solve.
1044 doublereal an = 1.0;
1045 doublereal bn = - GasConstant * T / pres;
1046 doublereal sqt = sqrt(T);
1047 doublereal cn = - (GasConstant * T * b / pres - a/(pres * sqt) + b * b);
1048 doublereal dn = - (a * b / (pres * sqt));
1049
1050 double tmp = a * omega_b / (b * omega_a * GasConstant);
1051 double pp = 2./3.;
1052 double tc = pow(tmp, pp);
1053 double pc = omega_b * GasConstant * tc / b;
1054 double vc = omega_vc * GasConstant * tc / pc;
1055
1056 int nSolnValues = MixtureFugacityTP::solveCubic(T, pres, a, b, a, Vroot, an, bn, cn, dn, tc, vc);
1057
1058 return nSolnValues;
1059}
1060
1061}
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:399
double convert(const std::string &key, const std::string &units) const
Convert the item stored by the given key to the units specified in units.
Definition: AnyMap.cpp:1508
const UnitSystem & units() const
Return the default units that should be used to convert stored values.
Definition: AnyMap.h:598
bool empty() const
Return boolean indicating whether AnyMap is empty.
Definition: AnyMap.cpp:1401
static AnyMap fromYamlFile(const std::string &name, const std::string &parent_name="")
Create an AnyMap from a YAML file.
Definition: AnyMap.cpp:1743
A wrapper for a variable whose type is determined at runtime.
Definition: AnyMap.h:84
void getRow(size_t n, double *const rw)
Get the nth row and return it in a vector.
Definition: Array.cpp:84
vector_fp & data()
Return a reference to the data vector.
Definition: Array.h:218
void resize(size_t n, size_t m, double v=0.0)
Resize the array, and fill the new entries with 'v'.
Definition: Array.cpp:52
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:61
vector_fp m_cp0_R
Temporary storage for dimensionless reference state heat capacities.
virtual bool addSpecies(shared_ptr< Species > spec)
doublereal z() const
Calculate the value of z.
virtual void _updateReferenceStateThermo() const
Updates the reference state thermodynamic functions at the current T of the solution.
virtual void getStandardVolumes(doublereal *vol) const
Get the molar volumes of each species in their standard states at the current T and P of the solution...
int solveCubic(double T, double pres, double a, double b, double aAlpha, double Vroot[3], double an, double bn, double cn, double dn, double tc, double vc) const
Solve the cubic equation of state.
virtual doublereal psatEst(doublereal TKelvin) const
Estimate for the saturation pressure.
virtual void getEntropy_R_ref(doublereal *er) const
Returns the vector of nondimensional entropies of the reference state at the current temperature of t...
virtual double critDensity() const
Critical density (kg/m3).
virtual void setTemperature(const doublereal temp)
Set the temperature of the phase.
vector_fp moleFractions_
Storage for the current values of the mole fractions of the species.
virtual double critTemperature() const
Critical temperature (K).
virtual void getGibbs_ref(doublereal *g) const
Returns the vector of the Gibbs function of the reference state at the current temperature of the sol...
vector_fp m_tmpV
Temporary storage - length = m_kk.
virtual double critPressure() const
Critical pressure (Pa).
virtual void getEnthalpy_RT_ref(doublereal *hrt) const
doublereal mean_X(const doublereal *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition: Phase.cpp:717
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range.
Definition: Phase.cpp:211
std::string name() const
Return the name of the phase.
Definition: Phase.cpp:70
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:273
size_t m_kk
Number of species in the phase.
Definition: Phase.h:943
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:200
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:751
double moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition: Phase.cpp:548
virtual double density() const
Density (kg/m^3).
Definition: Phase.h:679
doublereal temperature() const
Temperature (K).
Definition: Phase.h:654
size_t speciesIndex(const std::string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition: Phase.cpp:187
double molarVolume() const
Molar volume (m^3/kmol).
Definition: Phase.cpp:682
shared_ptr< Species > species(const std::string &name) const
Return the Species object for the named species.
Definition: Phase.cpp:950
virtual bool addSpecies(shared_ptr< Species > spec)
virtual void getActivityCoefficients(doublereal *ac) const
Get the array of non-dimensional activity coefficients at the current solution temperature,...
static const doublereal omega_b
Omega constant for b.
int m_formTempParam
Form of the temperature parameterization.
virtual doublereal hresid() const
Calculate the deviation terms for the total enthalpy of the mixture from the ideal gas mixture.
doublereal m_b_current
Value of b in the equation of state.
doublereal dpdT_
The derivative of the pressure wrt the temperature.
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
std::map< std::string, std::map< std::string, std::pair< double, double > > > m_binaryParameters
Explicitly-specified binary interaction parameters.
virtual void getPartialMolarIntEnergies(doublereal *ubar) const
Return an array of partial molar internal energies for the species in the mixture.
virtual doublereal cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
virtual std::vector< double > getCoeff(const std::string &iName)
Retrieve a and b coefficients by looking up tabulated critical parameters.
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
virtual void getPartialMolarEntropies(doublereal *sbar) const
Returns an array of partial molar entropies of the species in the solution.
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Import and initialize a ThermoPhase object using an XML tree.
void readXMLCrossFluid(XML_Node &pureFluidParam)
Read the cross species RedlichKwong input parameters.
virtual void getPartialMolarVolumes(doublereal *vbar) const
Return an array of partial molar volumes for the species in the mixture.
vector_fp m_pp
Temporary storage - length = m_kk.
RedlichKwongMFTP(const std::string &infile="", const std::string &id="")
Construct a RedlichKwongMFTP object from an input file.
virtual doublereal cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
void readXMLPureFluid(XML_Node &pureFluidParam)
Read the pure species RedlichKwong input parameters.
virtual doublereal dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal &presCalc) const
Calculate the pressure and the pressure derivative given the temperature and the molar volume.
std::vector< CoeffSource > m_coeffSource
For each species, specifies the source of the a and b coefficients.
virtual doublereal sresid() const
Calculate the deviation terms for the total entropy of the mixture from the ideal gas mixture.
virtual void updateMixingExpressions()
Update the a and b parameters.
void setBinaryCoeffs(const std::string &species_i, const std::string &species_j, double a0, double a1)
Set values for the interaction parameter between two species.
void pressureDerivatives() const
Calculate dpdV and dpdT at the current conditions.
static const doublereal omega_vc
Omega constant for the critical molar volume.
virtual doublereal liquidVolEst(doublereal TKelvin, doublereal &pres) const
Estimate for the molar volume of the liquid.
virtual void initThermo()
Initialize the ThermoPhase object after all species have been set up.
virtual doublereal densityCalc(doublereal T, doublereal pressure, int phase, doublereal rhoguess)
Calculates the density given the temperature and the pressure and a guess at the density.
doublereal m_a_current
Value of a in the equation of state.
virtual doublereal densSpinodalGas() const
Return the value of the density at the gas spinodal point (on the gas side) for the current temperatu...
void setSpeciesCoeffs(const std::string &species, double a0, double a1, double b)
Set the pure fluid interaction parameters for a species.
virtual void getSpeciesParameters(const std::string &name, AnyMap &speciesNode) const
Get phase-specific parameters of a Species object such that an identical one could be reconstructed a...
void calculateAB(doublereal temp, doublereal &aCalc, doublereal &bCalc) const
Calculate the a and the b parameters given the temperature.
virtual void getChemPotentials_RT(doublereal *mu) const
Get the array of non-dimensional species chemical potentials.
doublereal dpdV_
The derivative of the pressure wrt the volume.
virtual doublereal densSpinodalLiquid() const
Return the value of the density at the liquid spinodal point (on the liquid side) for the current tem...
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
static const doublereal omega_a
Omega constant for a -> value of a in terms of critical properties.
vector_fp dpdni_
Vector of derivatives of pressure wrt mole number.
virtual void setParametersFromXML(const XML_Node &thermoNode)
Set equation of state parameter values from XML entries.
int solveCubic(double T, double pres, double a, double b, double Vroot[3]) const
Prepare variables and call the function to solve the cubic equation of state.
virtual doublereal standardConcentration(size_t k=0) const
Returns the standard concentration , which is used to normalize the generalized concentration.
doublereal RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:782
void initThermoFile(const std::string &inputFile, const std::string &id)
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Import and initialize a ThermoPhase object using an XML tree.
virtual doublereal refPressure() const
Returns the reference pressure in Pa.
Definition: ThermoPhase.h:150
virtual void setParametersFromXML(const XML_Node &eosdata)
Set equation of state parameter values from XML entries.
Definition: ThermoPhase.h:1747
virtual void getSpeciesParameters(const std::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:1726
Unit conversion utility.
Definition: Units.h:161
double convert(double value, const std::string &src, const std::string &dest) const
Convert value from the units of src to the units of dest.
Definition: Units.cpp:558
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:103
std::string attrib(const std::string &attr) const
Function returns the value of an attribute.
Definition: xml.cpp:493
std::string name() const
Returns the name of the XML node.
Definition: xml.h:371
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
Definition: xml.cpp:529
size_t nChildren(bool discardComments=false) const
Return the number of children.
Definition: xml.cpp:557
XML_Node & child(const size_t n) const
Return a changeable reference to the n'th child of the current node.
Definition: xml.cpp:547
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data.
void importPhase(XML_Node &phase, ThermoPhase *th)
Import a phase information into an empty ThermoPhase object.
Namespace for the Cantera kernel.
Definition: AnyMap.h:29
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:192
doublereal getFloatCurrent(const XML_Node &currXML, const std::string &type="")
Get a floating-point value from the current XML element.
Definition: ctml.cpp:179
bool caseInsensitiveEquals(const std::string &input, const std::string &test)
Case insensitive equality predicate.
void warn_deprecated(const std::string &source, const AnyBase &node, const std::string &message)
A deprecation warning for syntax in an input file.
Definition: AnyMap.cpp:1901
const double SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:153
std::string toLowerCopy(const std::string &input)
Convert to lower case.
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:100
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:184
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition: ct_defs.h:113
size_t getFloatArray(const XML_Node &node, vector_fp &v, const bool convert=true, const std::string &unitsString="", const std::string &nodeName="floatArray")
This function reads the current node or a child node of the current node with the default name,...
Definition: ctml.cpp:258
Contains declarations for string manipulation functions within Cantera.