Cantera  3.1.0a1
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/utilities.h"
11 
12 #include <boost/algorithm/string.hpp>
13 #include <boost/math/tools/roots.hpp>
14 
15 namespace bmt = boost::math::tools;
16 
17 namespace Cantera
18 {
19 
20 const double RedlichKwongMFTP::omega_a = 4.27480233540E-01;
21 const double RedlichKwongMFTP::omega_b = 8.66403499650E-02;
22 const double RedlichKwongMFTP::omega_vc = 3.33333333333333E-01;
23 
24 RedlichKwongMFTP::RedlichKwongMFTP(const string& infile, const string& id_)
25 {
26  initThermoFile(infile, id_);
27 }
28 
29 void RedlichKwongMFTP::setSpeciesCoeffs(const string& species,
30  double a0, double a1, double b)
31 {
32  size_t k = speciesIndex(species);
33  if (k == npos) {
34  throw CanteraError("RedlichKwongMFTP::setSpeciesCoeffs",
35  "Unknown species '{}'.", species);
36  }
37 
38  if (a1 != 0.0) {
39  m_formTempParam = 1; // expression is temperature-dependent
40  }
41 
42  size_t counter = k + m_kk * k;
43  a_coeff_vec(0, counter) = a0;
44  a_coeff_vec(1, counter) = a1;
45 
46  // standard mixing rule for cross-species interaction term
47  for (size_t j = 0; j < m_kk; j++) {
48  if (k == j) {
49  continue;
50  }
51 
52  // a_coeff_vec(0) is initialized to NaN to mark uninitialized species
53  if (isnan(a_coeff_vec(0, j + m_kk * j))) {
54  // The diagonal element of the jth species has not yet been defined.
55  continue;
56  } else if (isnan(a_coeff_vec(0, j + m_kk * k))) {
57  // Only use the mixing rules if the off-diagonal element has not already been defined by a
58  // user-specified crossFluidParameters entry:
59  double a0kj = sqrt(a_coeff_vec(0, j + m_kk * j) * a0);
60  double a1kj = sqrt(a_coeff_vec(1, j + m_kk * j) * a1);
61  a_coeff_vec(0, j + m_kk * k) = a0kj;
62  a_coeff_vec(1, j + m_kk * k) = a1kj;
63  a_coeff_vec(0, k + m_kk * j) = a0kj;
64  a_coeff_vec(1, k + m_kk * j) = a1kj;
65  }
66  }
67  a_coeff_vec.getRow(0, a_vec_Curr_.data());
68  b_vec_Curr_[k] = b;
69 }
70 
71 void RedlichKwongMFTP::setBinaryCoeffs(const string& species_i, const string& species_j,
72  double a0, double a1)
73 {
74  size_t ki = speciesIndex(species_i);
75  if (ki == npos) {
76  throw CanteraError("RedlichKwongMFTP::setBinaryCoeffs",
77  "Unknown species '{}'.", species_i);
78  }
79  size_t kj = speciesIndex(species_j);
80  if (kj == npos) {
81  throw CanteraError("RedlichKwongMFTP::setBinaryCoeffs",
82  "Unknown species '{}'.", species_j);
83  }
84 
85  if (a1 != 0.0) {
86  m_formTempParam = 1; // expression is temperature-dependent
87  }
88 
89  m_binaryParameters[species_i][species_j] = {a0, a1};
90  m_binaryParameters[species_j][species_i] = {a0, a1};
91  size_t counter1 = ki + m_kk * kj;
92  size_t counter2 = kj + m_kk * ki;
93  a_coeff_vec(0, counter1) = a_coeff_vec(0, counter2) = a0;
94  a_coeff_vec(1, counter1) = a_coeff_vec(1, counter2) = a1;
95  a_vec_Curr_[counter1] = a_vec_Curr_[counter2] = a0;
96 }
97 
98 // ------------Molar Thermodynamic Properties -------------------------
99 
101 {
103  double TKelvin = temperature();
104  double sqt = sqrt(TKelvin);
105  double mv = molarVolume();
106  double vpb = mv + m_b_current;
108  double cpref = GasConstant * mean_X(m_cp0_R);
109  double dadt = da_dt();
110  double fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
111  double dHdT_V = (cpref + mv * dpdT_ - GasConstant - 1.0 / (2.0 * m_b_current * TKelvin * sqt) * log(vpb/mv) * fac
112  +1.0/(m_b_current * sqt) * log(vpb/mv) * (-0.5 * dadt));
113  return dHdT_V - (mv + TKelvin * dpdT_ / dpdV_) * dpdT_;
114 }
115 
117 {
119  double TKelvin = temperature();
120  double sqt = sqrt(TKelvin);
121  double mv = molarVolume();
122  double vpb = mv + m_b_current;
123  double cvref = GasConstant * (mean_X(m_cp0_R) - 1.0);
124  double dadt = da_dt();
125  double fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
126  return (cvref - 1.0/(2.0 * m_b_current * TKelvin * sqt) * log(vpb/mv)*fac
127  +1.0/(m_b_current * sqt) * log(vpb/mv)*(-0.5*dadt));
128 }
129 
131 {
133 
134  // Get a copy of the private variables stored in the State object
135  double T = temperature();
136  double molarV = meanMolecularWeight() / density();
137  double pp = GasConstant * T/(molarV - m_b_current) - m_a_current/(sqrt(T) * molarV * (molarV + m_b_current));
138  return pp;
139 }
140 
142 {
143  getStandardVolumes(m_tmpV.data());
144  return 1.0 / m_tmpV[k];
145 }
146 
148 {
149  double mv = molarVolume();
150  double sqt = sqrt(temperature());
151  double vpb = mv + m_b_current;
152  double vmb = mv - m_b_current;
153 
154  for (size_t k = 0; k < m_kk; k++) {
155  m_pp[k] = 0.0;
156  for (size_t i = 0; i < m_kk; i++) {
157  size_t counter = k + m_kk*i;
158  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
159  }
160  }
161  double pres = pressure();
162 
163  for (size_t k = 0; k < m_kk; k++) {
164  ac[k] = (- RT() * log(pres * mv / RT())
165  + RT() * log(mv / vmb)
166  + RT() * b_vec_Curr_[k] / vmb
167  - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
168  + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
169  - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
170  );
171  }
172  for (size_t k = 0; k < m_kk; k++) {
173  ac[k] = exp(ac[k]/RT());
174  }
175 }
176 
177 // ---- Partial Molar Properties of the Solution -----------------
178 
180 {
181  getGibbs_ref(mu);
182  for (size_t k = 0; k < m_kk; k++) {
183  double xx = std::max(SmallNumber, moleFraction(k));
184  mu[k] += RT()*(log(xx));
185  }
186 
187  double mv = molarVolume();
188  double sqt = sqrt(temperature());
189  double vpb = mv + m_b_current;
190  double vmb = mv - m_b_current;
191 
192  for (size_t k = 0; k < m_kk; k++) {
193  m_pp[k] = 0.0;
194  for (size_t i = 0; i < m_kk; i++) {
195  size_t counter = k + m_kk*i;
196  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
197  }
198  }
199  double pres = pressure();
200  double refP = refPressure();
201 
202  for (size_t k = 0; k < m_kk; k++) {
203  mu[k] += (RT() * log(pres/refP) - RT() * log(pres * mv / RT())
204  + RT() * log(mv / vmb)
205  + RT() * b_vec_Curr_[k] / vmb
206  - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
207  + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
208  - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
209  );
210  }
211 }
212 
214 {
215  // First we get the reference state contributions
216  getEnthalpy_RT_ref(hbar);
217  scale(hbar, hbar+m_kk, hbar, RT());
218 
219  // We calculate dpdni_
220  double TKelvin = temperature();
221  double mv = molarVolume();
222  double sqt = sqrt(TKelvin);
223  double vpb = mv + m_b_current;
224  double vmb = mv - m_b_current;
225  for (size_t k = 0; k < m_kk; k++) {
226  m_pp[k] = 0.0;
227  for (size_t i = 0; i < m_kk; i++) {
228  size_t counter = k + m_kk*i;
229  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
230  }
231  }
232  for (size_t k = 0; k < m_kk; k++) {
233  dpdni_[k] = RT()/vmb + RT() * b_vec_Curr_[k] / (vmb * vmb) - 2.0 * m_pp[k] / (sqt * mv * vpb)
234  + m_a_current * b_vec_Curr_[k]/(sqt * mv * vpb * vpb);
235  }
236  double dadt = da_dt();
237  double fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
238 
239  for (size_t k = 0; k < m_kk; k++) {
240  m_tmpV[k] = 0.0;
241  for (size_t i = 0; i < m_kk; i++) {
242  size_t counter = k + m_kk*i;
243  m_tmpV[k] += 2.0 * moleFractions_[i] * TKelvin * a_coeff_vec(1,counter) - 3.0 * moleFractions_[i] * a_vec_Curr_[counter];
244  }
245  }
246 
248  double fac2 = mv + TKelvin * dpdT_ / dpdV_;
249  for (size_t k = 0; k < m_kk; k++) {
250  double hE_v = (mv * dpdni_[k] - RT() - b_vec_Curr_[k]/ (m_b_current * m_b_current * sqt) * log(vpb/mv)*fac
251  + 1.0 / (m_b_current * sqt) * log(vpb/mv) * m_tmpV[k]
252  + b_vec_Curr_[k] / vpb / (m_b_current * sqt) * fac);
253  hbar[k] = hbar[k] + hE_v;
254  hbar[k] -= fac2 * dpdni_[k];
255  }
256 }
257 
259 {
260  getEntropy_R_ref(sbar);
261  scale(sbar, sbar+m_kk, sbar, GasConstant);
262  double TKelvin = temperature();
263  double sqt = sqrt(TKelvin);
264  double mv = molarVolume();
265  double refP = refPressure();
266 
267  for (size_t k = 0; k < m_kk; k++) {
268  double xx = std::max(SmallNumber, moleFraction(k));
269  sbar[k] += GasConstant * (- log(xx));
270  }
271  for (size_t k = 0; k < m_kk; k++) {
272  m_pp[k] = 0.0;
273  for (size_t i = 0; i < m_kk; i++) {
274  size_t counter = k + m_kk*i;
275  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
276  }
277  }
278  for (size_t k = 0; k < m_kk; k++) {
279  m_tmpV[k] = 0.0;
280  for (size_t i = 0; i < m_kk; i++) {
281  size_t counter = k + m_kk*i;
282  m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
283  }
284  }
285 
286  double dadt = da_dt();
287  double fac = dadt - m_a_current / (2.0 * TKelvin);
288  double vmb = mv - m_b_current;
289  double vpb = mv + m_b_current;
290  for (size_t k = 0; k < m_kk; k++) {
291  sbar[k] -=(GasConstant * log(GasConstant * TKelvin / (refP * mv))
292  + GasConstant
293  + GasConstant * log(mv/vmb)
294  + GasConstant * b_vec_Curr_[k]/vmb
295  + m_pp[k]/(m_b_current * TKelvin * sqt) * log(vpb/mv)
296  - 2.0 * m_tmpV[k]/(m_b_current * sqt) * log(vpb/mv)
297  + b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) * fac
298  - 1.0 / (m_b_current * sqt) * b_vec_Curr_[k] / vpb * fac
299  );
300  }
301 
303  getPartialMolarVolumes(m_partialMolarVolumes.data());
304  for (size_t k = 0; k < m_kk; k++) {
305  sbar[k] -= -m_partialMolarVolumes[k] * dpdT_;
306  }
307 }
308 
310 {
311  // u_k = h_k - P * v_k
312  getPartialMolarVolumes(m_partialMolarVolumes.data());
314  double p = pressure();
315  for (size_t k = 0; k < nSpecies(); k++) {
316  ubar[k] -= p * m_partialMolarVolumes[k];
317  }
318 }
319 
321 {
322  for (size_t k = 0; k < m_kk; k++) {
323  m_pp[k] = 0.0;
324  for (size_t i = 0; i < m_kk; i++) {
325  size_t counter = k + m_kk*i;
326  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
327  }
328  }
329  for (size_t k = 0; k < m_kk; k++) {
330  m_tmpV[k] = 0.0;
331  for (size_t i = 0; i < m_kk; i++) {
332  size_t counter = k + m_kk*i;
333  m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
334  }
335  }
336 
337  double sqt = sqrt(temperature());
338  double mv = molarVolume();
339  double vmb = mv - m_b_current;
340  double vpb = mv + m_b_current;
341  for (size_t k = 0; k < m_kk; k++) {
342  double num = (RT() + RT() * m_b_current/ vmb + RT() * b_vec_Curr_[k] / vmb
343  + RT() * m_b_current * b_vec_Curr_[k] /(vmb * vmb)
344  - 2.0 * m_pp[k] / (sqt * vpb)
345  + m_a_current * b_vec_Curr_[k] / (sqt * vpb * vpb)
346  );
347  double denom = (pressure() + RT() * m_b_current/(vmb * vmb) - m_a_current / (sqt * vpb * vpb)
348  );
349  vbar[k] = num / denom;
350  }
351 }
352 
353 bool RedlichKwongMFTP::addSpecies(shared_ptr<Species> spec)
354 {
355  bool added = MixtureFugacityTP::addSpecies(spec);
356  if (added) {
357  a_vec_Curr_.resize(m_kk * m_kk, 0.0);
358 
359  // Initialize a_vec and b_vec to NaN, to screen for species with
360  // pureFluidParameters which are undefined in the input file:
361  b_vec_Curr_.push_back(NAN);
362  a_coeff_vec.resize(2, m_kk * m_kk, NAN);
363 
364  m_pp.push_back(0.0);
365  m_coeffSource.push_back(CoeffSource::EoS);
366  m_partialMolarVolumes.push_back(0.0);
367  dpdni_.push_back(0.0);
368  }
369  return added;
370 }
371 
373 {
374  // Contents of 'critical-properties.yaml', loaded later if needed
375  AnyMap critPropsDb;
376  std::unordered_map<string, AnyMap*> dbSpecies;
377 
378  for (auto& [name, species] : m_species) {
379  auto& data = species->input;
380  size_t k = speciesIndex(name);
381  if (!isnan(a_coeff_vec(0, k + m_kk * k))) {
382  continue;
383  }
384  bool foundCoeffs = false;
385 
386  if (data.hasKey("equation-of-state") &&
387  data["equation-of-state"].hasMapWhere("model", "Redlich-Kwong"))
388  {
389  // Read a and b coefficients from species 'input' information (that is, as
390  // specified in a YAML input file) specific to the Redlich-Kwong model
391  auto eos = data["equation-of-state"].getMapWhere(
392  "model", "Redlich-Kwong");
393 
394  if (eos.hasKey("a") && eos.hasKey("b")) {
395  double a0 = 0, a1 = 0;
396  if (eos["a"].isScalar()) {
397  a0 = eos.convert("a", "Pa*m^6/kmol^2*K^0.5");
398  } else {
399  auto avec = eos["a"].asVector<AnyValue>(2);
400  a0 = eos.units().convert(avec[0], "Pa*m^6/kmol^2*K^0.5");
401  a1 = eos.units().convert(avec[1], "Pa*m^6/kmol^2/K^0.5");
402  }
403  double b = eos.convert("b", "m^3/kmol");
404  foundCoeffs = true;
405  setSpeciesCoeffs(name, a0, a1, b);
406  m_coeffSource[k] = CoeffSource::EoS;
407  }
408 
409  if (eos.hasKey("binary-a")) {
410  AnyMap& binary_a = eos["binary-a"].as<AnyMap>();
411  const UnitSystem& units = binary_a.units();
412  for (auto& [name2, coeff] : binary_a) {
413  double a0 = 0, a1 = 0;
414  if (coeff.isScalar()) {
415  a0 = units.convert(coeff, "Pa*m^6/kmol^2*K^0.5");
416  } else {
417  auto avec = coeff.asVector<AnyValue>(2);
418  a0 = units.convert(avec[0], "Pa*m^6/kmol^2*K^0.5");
419  a1 = units.convert(avec[1], "Pa*m^6/kmol^2/K^0.5");
420  }
421  setBinaryCoeffs(name, name2, a0, a1);
422  }
423  }
424 
425  if (foundCoeffs) {
426  continue;
427  }
428  }
429 
430  // Coefficients have not been populated from model-specific input
431  double Tc = NAN, Pc = NAN;
432  if (data.hasKey("critical-parameters")) {
433  // Use critical state information stored in the species entry to
434  // calculate a and b
435  auto& critProps = data["critical-parameters"].as<AnyMap>();
436  Tc = critProps.convert("critical-temperature", "K");
437  Pc = critProps.convert("critical-pressure", "Pa");
438  m_coeffSource[k] = CoeffSource::CritProps;
439  } else {
440  // Search 'crit-properties.yaml' to find Tc and Pc. Load data if needed
441  if (critPropsDb.empty()) {
442  critPropsDb = AnyMap::fromYamlFile("critical-properties.yaml");
443  dbSpecies = critPropsDb["species"].asMap("name");
444  }
445 
446  // All names in critical-properties.yaml are upper case
447  auto ucName = boost::algorithm::to_upper_copy(name);
448  if (dbSpecies.count(ucName)) {
449  auto& spec = *dbSpecies.at(ucName);
450  auto& critProps = spec["critical-parameters"].as<AnyMap>();
451  Tc = critProps.convert("critical-temperature", "K");
452  Pc = critProps.convert("critical-pressure", "Pa");
453  m_coeffSource[k] = CoeffSource::Database;
454  }
455  }
456 
457  // Check if critical properties were found in either location
458  if (!isnan(Tc)) {
459  // Assuming no temperature dependence (that is, a1 = 0)
460  double a = omega_a * pow(GasConstant, 2) * pow(Tc, 2.5) / Pc;
461  double b = omega_b * GasConstant * Tc / Pc;
462  setSpeciesCoeffs(name, a, 0.0, b);
463  } else {
464  throw InputFileError("RedlichKwongMFTP::initThermo", data,
465  "No critical property or Redlich-Kwong parameters found "
466  "for species {}.", name);
467  }
468  }
469 }
470 
472  AnyMap& speciesNode) const
473 {
475  size_t k = speciesIndex(name);
477  if (m_coeffSource[k] == CoeffSource::EoS) {
478  auto& eosNode = speciesNode["equation-of-state"].getMapWhere(
479  "model", "Redlich-Kwong", true);
480 
481  size_t counter = k + m_kk * k;
482  if (a_coeff_vec(1, counter) != 0.0) {
483  vector<AnyValue> coeffs(2);
484  coeffs[0].setQuantity(a_coeff_vec(0, counter), "Pa*m^6/kmol^2*K^0.5");
485  coeffs[1].setQuantity(a_coeff_vec(1, counter), "Pa*m^6/kmol^2/K^0.5");
486  eosNode["a"] = std::move(coeffs);
487  } else {
488  eosNode["a"].setQuantity(a_coeff_vec(0, counter),
489  "Pa*m^6/kmol^2*K^0.5");
490  }
491  eosNode["b"].setQuantity(b_vec_Curr_[k], "m^3/kmol");
492  } else if (m_coeffSource[k] == CoeffSource::CritProps) {
493  auto& critProps = speciesNode["critical-parameters"];
494  double a = a_coeff_vec(0, k + m_kk * k);
495  double b = b_vec_Curr_[k];
496  double Tc = pow(a * omega_b / (b * omega_a * GasConstant), 2.0/3.0);
497  double Pc = omega_b * GasConstant * Tc / b;
498  critProps["critical-temperature"].setQuantity(Tc, "K");
499  critProps["critical-pressure"].setQuantity(Pc, "Pa");
500  }
501 
502  if (m_binaryParameters.count(name)) {
503  auto& eosNode = speciesNode["equation-of-state"].getMapWhere(
504  "model", "Redlich-Kwong", true);
505  AnyMap bin_a;
506  for (const auto& [name2, coeffs] : m_binaryParameters.at(name)) {
507  if (coeffs.second == 0) {
508  bin_a[name2].setQuantity(coeffs.first, "Pa*m^6/kmol^2*K^0.5");
509  } else {
510  vector<AnyValue> C(2);
511  C[0].setQuantity(coeffs.first, "Pa*m^6/kmol^2*K^0.5");
512  C[1].setQuantity(coeffs.second, "Pa*m^6/kmol^2/K^0.5");
513  bin_a[name2] = std::move(C);
514  }
515  }
516  eosNode["binary-a"] = std::move(bin_a);
517  }
518 }
519 
521 {
522  // note this agrees with tpx
523  double rho = density();
524  double mmw = meanMolecularWeight();
525  double molarV = mmw / rho;
526  double hh = m_b_current / molarV;
527  double zz = z();
528  double dadt = da_dt();
529  double T = temperature();
530  double sqT = sqrt(T);
531  double fac = dadt - m_a_current / (2.0 * T);
532  double sresid_mol_R = log(zz*(1.0 - hh)) + log(1.0 + hh) * fac / (sqT * GasConstant * m_b_current);
533  return GasConstant * sresid_mol_R;
534 }
535 
537 {
538  // note this agrees with tpx
539  double rho = density();
540  double mmw = meanMolecularWeight();
541  double molarV = mmw / rho;
542  double hh = m_b_current / molarV;
543  double zz = z();
544  double dadt = da_dt();
545  double T = temperature();
546  double sqT = sqrt(T);
547  double fac = T * dadt - 3.0 *m_a_current / (2.0);
548  return GasConstant * T * (zz - 1.0) + fac * log(1.0 + hh) / (sqT * m_b_current);
549 }
550 
551 double RedlichKwongMFTP::liquidVolEst(double TKelvin, double& presGuess) const
552 {
553  double v = m_b_current * 1.1;
554  double atmp;
555  double btmp;
556  calculateAB(TKelvin, atmp, btmp);
557  double pres = std::max(psatEst(TKelvin), presGuess);
558  double Vroot[3];
559  bool foundLiq = false;
560  int m = 0;
561  while (m < 100 && !foundLiq) {
562  int nsol = solveCubic(TKelvin, pres, atmp, btmp, Vroot);
563  if (nsol == 1 || nsol == 2) {
564  double pc = critPressure();
565  if (pres > pc) {
566  foundLiq = true;
567  }
568  pres *= 1.04;
569  } else {
570  foundLiq = true;
571  }
572  }
573 
574  if (foundLiq) {
575  v = Vroot[0];
576  presGuess = pres;
577  } else {
578  v = -1.0;
579  }
580  return v;
581 }
582 
583 double RedlichKwongMFTP::densityCalc(double TKelvin, double presPa, int phaseRequested,
584  double rhoguess)
585 {
586  // It's necessary to set the temperature so that m_a_current is set correctly.
587  setTemperature(TKelvin);
588  double tcrit = critTemperature();
589  double mmw = meanMolecularWeight();
590  if (rhoguess == -1.0) {
591  if (phaseRequested >= FLUID_LIQUID_0) {
592  double lqvol = liquidVolEst(TKelvin, presPa);
593  rhoguess = mmw / lqvol;
594  }
595  } else {
596  // Assume the Gas phase initial guess, if nothing is specified to the routine
597  rhoguess = presPa * mmw / (GasConstant * TKelvin);
598  }
599 
600 
601  double volguess = mmw / rhoguess;
602  NSolns_ = solveCubic(TKelvin, presPa, m_a_current, m_b_current, Vroot_);
603 
604  double molarVolLast = Vroot_[0];
605  if (NSolns_ >= 2) {
606  if (phaseRequested >= FLUID_LIQUID_0) {
607  molarVolLast = Vroot_[0];
608  } else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
609  molarVolLast = Vroot_[2];
610  } else {
611  if (volguess > Vroot_[1]) {
612  molarVolLast = Vroot_[2];
613  } else {
614  molarVolLast = Vroot_[0];
615  }
616  }
617  } else if (NSolns_ == 1) {
618  if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT || phaseRequested == FLUID_UNDEFINED) {
619  molarVolLast = Vroot_[0];
620  } else {
621  return -2.0;
622  }
623  } else if (NSolns_ == -1) {
624  if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) {
625  molarVolLast = Vroot_[0];
626  } else if (TKelvin > tcrit) {
627  molarVolLast = Vroot_[0];
628  } else {
629  return -2.0;
630  }
631  } else {
632  molarVolLast = Vroot_[0];
633  return -1.0;
634  }
635  return mmw / molarVolLast;
636 }
637 
639 {
640  double Vroot[3];
641  double T = temperature();
642  int nsol = solveCubic(T, pressure(), m_a_current, m_b_current, Vroot);
643  if (nsol != 3) {
644  return critDensity();
645  }
646 
647  auto resid = [this, T](double v) {
648  double pp;
649  return dpdVCalc(T, v, pp);
650  };
651 
652  boost::uintmax_t maxiter = 100;
653  auto [lower, upper] = bmt::toms748_solve(
654  resid, Vroot[0], Vroot[1], bmt::eps_tolerance<double>(48), maxiter);
655 
656  double mmw = meanMolecularWeight();
657  return mmw / (0.5 * (lower + upper));
658 }
659 
661 {
662  double Vroot[3];
663  double T = temperature();
664  int nsol = solveCubic(T, pressure(), m_a_current, m_b_current, Vroot);
665  if (nsol != 3) {
666  return critDensity();
667  }
668 
669  auto resid = [this, T](double v) {
670  double pp;
671  return dpdVCalc(T, v, pp);
672  };
673 
674  boost::uintmax_t maxiter = 100;
675  auto [lower, upper] = bmt::toms748_solve(
676  resid, Vroot[1], Vroot[2], bmt::eps_tolerance<double>(48), maxiter);
677 
678  double mmw = meanMolecularWeight();
679  return mmw / (0.5 * (lower + upper));
680 }
681 
682 double RedlichKwongMFTP::dpdVCalc(double TKelvin, double molarVol, double& presCalc) const
683 {
684  double sqt = sqrt(TKelvin);
685  presCalc = GasConstant * TKelvin / (molarVol - m_b_current)
686  - m_a_current / (sqt * molarVol * (molarVol + m_b_current));
687 
688  double vpb = molarVol + m_b_current;
689  double vmb = molarVol - m_b_current;
690  double dpdv = (- GasConstant * TKelvin / (vmb * vmb)
691  + m_a_current * (2 * molarVol + m_b_current) / (sqt * molarVol * molarVol * vpb * vpb));
692  return dpdv;
693 }
694 
696 {
698  return -1 / (molarVolume() * dpdV_);
699 }
700 
702 {
704  return -dpdT_ / (molarVolume() * dpdV_);
705 }
706 
708 {
710  return molarVolume() * sqrt(-cp_mole() / cv_mole() * dpdV_ / meanMolecularWeight());
711 }
712 
714 {
715  double TKelvin = temperature();
716  double mv = molarVolume();
717  double pres;
718 
719  dpdV_ = dpdVCalc(TKelvin, mv, pres);
720  double sqt = sqrt(TKelvin);
721  double vpb = mv + m_b_current;
722  double vmb = mv - m_b_current;
723  double dadt = da_dt();
724  double fac = dadt - m_a_current/(2.0 * TKelvin);
725  dpdT_ = (GasConstant / vmb - fac / (sqt * mv * vpb));
726 }
727 
729 {
730  double temp = temperature();
731  if (m_formTempParam == 1) {
732  for (size_t i = 0; i < m_kk; i++) {
733  for (size_t j = 0; j < m_kk; j++) {
734  size_t counter = i * m_kk + j;
735  a_vec_Curr_[counter] = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
736  }
737  }
738  }
739 
740  m_b_current = 0.0;
741  m_a_current = 0.0;
742  for (size_t i = 0; i < m_kk; i++) {
743  m_b_current += moleFractions_[i] * b_vec_Curr_[i];
744  for (size_t j = 0; j < m_kk; j++) {
745  m_a_current += a_vec_Curr_[i * m_kk + j] * moleFractions_[i] * moleFractions_[j];
746  }
747  }
748  if (isnan(m_b_current)) {
749  // One or more species do not have specified coefficients.
750  fmt::memory_buffer b;
751  for (size_t k = 0; k < m_kk; k++) {
752  if (isnan(b_vec_Curr_[k])) {
753  if (b.size() > 0) {
754  fmt_append(b, ", {}", speciesName(k));
755  } else {
756  fmt_append(b, "{}", speciesName(k));
757  }
758  }
759  }
760  throw CanteraError("RedlichKwongMFTP::updateMixingExpressions",
761  "Missing Redlich-Kwong coefficients for species: {}", to_string(b));
762  }
763 }
764 
765 void RedlichKwongMFTP::calculateAB(double temp, double& aCalc, double& bCalc) const
766 {
767  bCalc = 0.0;
768  aCalc = 0.0;
769  if (m_formTempParam == 1) {
770  for (size_t i = 0; i < m_kk; i++) {
771  bCalc += moleFractions_[i] * b_vec_Curr_[i];
772  for (size_t j = 0; j < m_kk; j++) {
773  size_t counter = i * m_kk + j;
774  double a_vec_Curr = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
775  aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
776  }
777  }
778  } else {
779  for (size_t i = 0; i < m_kk; i++) {
780  bCalc += moleFractions_[i] * b_vec_Curr_[i];
781  for (size_t j = 0; j < m_kk; j++) {
782  size_t counter = i * m_kk + j;
783  double a_vec_Curr = a_coeff_vec(0,counter);
784  aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
785  }
786  }
787  }
788 }
789 
790 double RedlichKwongMFTP::da_dt() const
791 {
792  double dadT = 0.0;
793  if (m_formTempParam == 1) {
794  for (size_t i = 0; i < m_kk; i++) {
795  for (size_t j = 0; j < m_kk; j++) {
796  size_t counter = i * m_kk + j;
797  dadT+= a_coeff_vec(1,counter) * moleFractions_[i] * moleFractions_[j];
798  }
799  }
800  }
801  return dadT;
802 }
803 
804 void RedlichKwongMFTP::calcCriticalConditions(double& pc, double& tc, double& vc) const
805 {
806  double a0 = 0.0;
807  double aT = 0.0;
808  for (size_t i = 0; i < m_kk; i++) {
809  for (size_t j = 0; j <m_kk; j++) {
810  size_t counter = i + m_kk * j;
811  a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
812  aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
813  }
814  }
815  double a = m_a_current;
816  double b = m_b_current;
817  if (m_formTempParam != 0) {
818  a = a0;
819  }
820  if (b <= 0.0) {
821  tc = 1000000.;
822  pc = 1.0E13;
823  vc = omega_vc * GasConstant * tc / pc;
824  return;
825  }
826  if (a <= 0.0) {
827  tc = 0.0;
828  pc = 0.0;
829  vc = 2.0 * b;
830  return;
831  }
832  double tmp = a * omega_b / (b * omega_a * GasConstant);
833  double pp = 2./3.;
834  double sqrttc, f, dfdt, deltatc;
835 
836  if (m_formTempParam == 0) {
837  tc = pow(tmp, pp);
838  } else {
839  tc = pow(tmp, pp);
840  for (int j = 0; j < 10; j++) {
841  sqrttc = sqrt(tc);
842  f = omega_a * b * GasConstant * tc * sqrttc / omega_b - aT * tc - a0;
843  dfdt = 1.5 * omega_a * b * GasConstant * sqrttc / omega_b - aT;
844  deltatc = - f / dfdt;
845  tc += deltatc;
846  }
847  if (deltatc > 0.1) {
848  throw CanteraError("RedlichKwongMFTP::calcCriticalConditions",
849  "didn't converge");
850  }
851  }
852 
853  pc = omega_b * GasConstant * tc / b;
854  vc = omega_vc * GasConstant * tc / pc;
855 }
856 
857 int RedlichKwongMFTP::solveCubic(double T, double pres, double a, double b, double Vroot[3]) const
858 {
859 
860  // Derive the coefficients of the cubic polynomial to solve.
861  double an = 1.0;
862  double bn = - GasConstant * T / pres;
863  double sqt = sqrt(T);
864  double cn = - (GasConstant * T * b / pres - a/(pres * sqt) + b * b);
865  double dn = - (a * b / (pres * sqt));
866 
867  double tmp = a * omega_b / (b * omega_a * GasConstant);
868  double pp = 2./3.;
869  double tc = pow(tmp, pp);
870  double pc = omega_b * GasConstant * tc / b;
871  double vc = omega_vc * GasConstant * tc / pc;
872 
873  int nSolnValues = MixtureFugacityTP::solveCubic(T, pres, a, b, a, Vroot, an, bn, cn, dn, tc, vc);
874 
875  return nSolnValues;
876 }
877 
878 }
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 empty() const
Return boolean indicating whether AnyMap is empty.
Definition: AnyMap.cpp:1418
double convert(const string &key, const string &units) const
Convert the item stored by the given key to the units specified in units.
Definition: AnyMap.cpp:1535
const UnitSystem & units() const
Return the default units that should be used to convert stored values.
Definition: AnyMap.h:630
static AnyMap fromYamlFile(const string &name, const string &parent_name="")
Create an AnyMap from a YAML file.
Definition: AnyMap.cpp:1771
A wrapper for a variable whose type is determined at runtime.
Definition: AnyMap.h:86
void getRow(size_t n, double *const rw)
Get the nth row and return it in a vector.
Definition: Array.cpp:79
virtual 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:47
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
Error thrown for problems processing information contained in an AnyMap or AnyValue.
Definition: AnyMap.h:738
double critPressure() const override
Critical pressure (Pa).
double critDensity() const override
Critical density (kg/m3).
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 critTemperature() const override
Critical temperature (K).
virtual void _updateReferenceStateThermo() const
Updates the reference state thermodynamic functions at the current T of the solution.
vector< double > m_tmpV
Temporary storage - length = m_kk.
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.
vector< double > moleFractions_
Storage for the current values of the mole fractions of the species.
void getEntropy_R_ref(double *er) const override
Returns the vector of nondimensional entropies of the reference state at the current temperature of t...
void setTemperature(const double temp) override
Set the temperature of the phase.
virtual double psatEst(double TKelvin) const
Estimate for the saturation pressure.
void getStandardVolumes(double *vol) const override
Get the molar volumes of each species in their standard states at the current T and P of the solution...
vector< double > m_cp0_R
Temporary storage for dimensionless reference state heat capacities.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
double z() const
Calculate the value of z.
void getEnthalpy_RT_ref(double *hrt) const override
Returns the vector of nondimensional enthalpies of the reference state at the current temperature of ...
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range.
Definition: Phase.cpp:153
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:231
size_t m_kk
Number of species in the phase.
Definition: Phase.h:842
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
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 mean_X(const double *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition: Phase.cpp:616
shared_ptr< Species > species(const string &name) const
Return the Species object for the named species.
Definition: Phase.cpp:856
virtual double molarVolume() const
Molar volume (m^3/kmol).
Definition: Phase.cpp:581
string name() const
Return the name of the phase.
Definition: Phase.cpp:20
double dpdT_
The derivative of the pressure wrt the temperature.
void calculateAB(double temp, double &aCalc, double &bCalc) const
Calculate the a and the b parameters given the temperature.
double thermalExpansionCoeff() const override
Return the volumetric thermal expansion coefficient. Units: 1/K.
void setBinaryCoeffs(const string &species_i, const string &species_j, double a0, double a1)
Set values for the interaction parameter between two species.
double densSpinodalLiquid() const override
Return the value of the density at the liquid spinodal point (on the liquid side) for the current tem...
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.
double sresid() const override
Calculate the deviation terms for the total entropy of the mixture from the ideal gas mixture.
double soundSpeed() const override
Return the speed of sound. Units: m/s.
double pressure() const override
Return the thermodynamic pressure (Pa).
int m_formTempParam
Form of the temperature parameterization.
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...
RedlichKwongMFTP(const string &infile="", const string &id="")
Construct a RedlichKwongMFTP object from an input file.
static const double omega_b
Omega constant for b.
vector< double > m_pp
Temporary storage - length = m_kk.
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
double densSpinodalGas() const override
Return the value of the density at the gas spinodal point (on the gas side) for the current temperatu...
void getPartialMolarVolumes(double *vbar) const override
Return an array of partial molar volumes for the species in the mixture.
double cv_mole() const override
Molar heat capacity at constant volume. Units: J/kmol/K.
void pressureDerivatives() const
Calculate dpdV and dpdT at the current conditions.
double isothermalCompressibility() const override
Returns the isothermal compressibility. Units: 1/Pa.
double hresid() const override
Calculate the deviation terms for the total enthalpy of the mixture from the ideal gas mixture.
static const double omega_a
Omega constant for a -> value of a in terms of critical properties.
double liquidVolEst(double TKelvin, double &pres) const override
Estimate for the molar volume of the liquid.
double dpdVCalc(double TKelvin, double molarVol, double &presCalc) const override
Calculate the pressure and the pressure derivative given the temperature and the molar volume.
static const double omega_vc
Omega constant for the critical molar volume.
vector< CoeffSource > m_coeffSource
For each species, specifies the source of the a and b coefficients.
void getPartialMolarIntEnergies(double *ubar) const override
Return an array of partial molar internal energies for the species in the mixture.
void updateMixingExpressions() override
Update the a and b parameters.
double cp_mole() const override
Molar heat capacity at constant pressure. Units: J/kmol/K.
double m_a_current
Value of a in the equation of state.
double standardConcentration(size_t k=0) const override
Returns the standard concentration , which is used to normalize the generalized concentration.
double dpdV_
The derivative of the pressure wrt the volume.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
double m_b_current
Value of b in the equation of state.
void getActivityCoefficients(double *ac) const override
Get the array of non-dimensional activity coefficients at the current solution temperature,...
double densityCalc(double T, double pressure, int phase, double rhoguess) override
Calculates the density given the temperature and the pressure and a guess at the density.
vector< double > dpdni_
Vector of derivatives of pressure wrt mole number.
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies of the species in the solution.
map< string, map< string, pair< double, double > > > m_binaryParameters
Explicitly-specified binary interaction parameters.
void setSpeciesCoeffs(const string &species, double a0, double a1, double b)
Set the pure fluid interaction parameters for a species.
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.
double RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:1062
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
virtual double refPressure() const
Returns the reference pressure in Pa.
Definition: ThermoPhase.h:436
Unit conversion utility.
Definition: Units.h:169
double convert(double value, const string &src, const string &dest) const
Convert value from the units of src to the units of dest.
Definition: Units.cpp:538
void fmt_append(fmt::memory_buffer &b, 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
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:104
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition: ct_defs.h:120
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:180
const double 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...