Cantera  3.0.0
Loading...
Searching...
No Matches
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
11
12#include <boost/algorithm/string.hpp>
13#include <boost/math/tools/roots.hpp>
14
15namespace bmt = boost::math::tools;
16
17namespace Cantera
18{
19
20const double RedlichKwongMFTP::omega_a = 4.27480233540E-01;
21const double RedlichKwongMFTP::omega_b = 8.66403499650E-02;
22const double RedlichKwongMFTP::omega_vc = 3.33333333333333E-01;
23
24RedlichKwongMFTP::RedlichKwongMFTP(const string& infile, const string& id_)
25{
26 initThermoFile(infile, id_);
27}
28
29void 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
71void 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{
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 warn_deprecated("RedlichKwongMFTP::getChemPotentials_RT",
182 "To be removed after Cantera 3.0. Use getChemPotentials instead.");
183 getChemPotentials(muRT);
184 for (size_t k = 0; k < m_kk; k++) {
185 muRT[k] *= 1.0 / RT();
186 }
187}
188
190{
191 getGibbs_ref(mu);
192 for (size_t k = 0; k < m_kk; k++) {
193 double xx = std::max(SmallNumber, moleFraction(k));
194 mu[k] += RT()*(log(xx));
195 }
196
197 double mv = molarVolume();
198 double sqt = sqrt(temperature());
199 double vpb = mv + m_b_current;
200 double vmb = mv - m_b_current;
201
202 for (size_t k = 0; k < m_kk; k++) {
203 m_pp[k] = 0.0;
204 for (size_t i = 0; i < m_kk; i++) {
205 size_t counter = k + m_kk*i;
206 m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
207 }
208 }
209 double pres = pressure();
210 double refP = refPressure();
211
212 for (size_t k = 0; k < m_kk; k++) {
213 mu[k] += (RT() * log(pres/refP) - RT() * log(pres * mv / RT())
214 + RT() * log(mv / vmb)
215 + RT() * b_vec_Curr_[k] / vmb
216 - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
217 + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
218 - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
219 );
220 }
221}
222
224{
225 // First we get the reference state contributions
226 getEnthalpy_RT_ref(hbar);
227 scale(hbar, hbar+m_kk, hbar, RT());
228
229 // We calculate dpdni_
230 double TKelvin = temperature();
231 double mv = molarVolume();
232 double sqt = sqrt(TKelvin);
233 double vpb = mv + m_b_current;
234 double vmb = mv - m_b_current;
235 for (size_t k = 0; k < m_kk; k++) {
236 m_pp[k] = 0.0;
237 for (size_t i = 0; i < m_kk; i++) {
238 size_t counter = k + m_kk*i;
239 m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
240 }
241 }
242 for (size_t k = 0; k < m_kk; k++) {
243 dpdni_[k] = RT()/vmb + RT() * b_vec_Curr_[k] / (vmb * vmb) - 2.0 * m_pp[k] / (sqt * mv * vpb)
244 + m_a_current * b_vec_Curr_[k]/(sqt * mv * vpb * vpb);
245 }
246 double dadt = da_dt();
247 double fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
248
249 for (size_t k = 0; k < m_kk; k++) {
250 m_tmpV[k] = 0.0;
251 for (size_t i = 0; i < m_kk; i++) {
252 size_t counter = k + m_kk*i;
253 m_tmpV[k] += 2.0 * moleFractions_[i] * TKelvin * a_coeff_vec(1,counter) - 3.0 * moleFractions_[i] * a_vec_Curr_[counter];
254 }
255 }
256
258 double fac2 = mv + TKelvin * dpdT_ / dpdV_;
259 for (size_t k = 0; k < m_kk; k++) {
260 double hE_v = (mv * dpdni_[k] - RT() - b_vec_Curr_[k]/ (m_b_current * m_b_current * sqt) * log(vpb/mv)*fac
261 + 1.0 / (m_b_current * sqt) * log(vpb/mv) * m_tmpV[k]
262 + b_vec_Curr_[k] / vpb / (m_b_current * sqt) * fac);
263 hbar[k] = hbar[k] + hE_v;
264 hbar[k] -= fac2 * dpdni_[k];
265 }
266}
267
269{
270 getEntropy_R_ref(sbar);
271 scale(sbar, sbar+m_kk, sbar, GasConstant);
272 double TKelvin = temperature();
273 double sqt = sqrt(TKelvin);
274 double mv = molarVolume();
275 double refP = refPressure();
276
277 for (size_t k = 0; k < m_kk; k++) {
278 double xx = std::max(SmallNumber, moleFraction(k));
279 sbar[k] += GasConstant * (- log(xx));
280 }
281 for (size_t k = 0; k < m_kk; k++) {
282 m_pp[k] = 0.0;
283 for (size_t i = 0; i < m_kk; i++) {
284 size_t counter = k + m_kk*i;
285 m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
286 }
287 }
288 for (size_t k = 0; k < m_kk; k++) {
289 m_tmpV[k] = 0.0;
290 for (size_t i = 0; i < m_kk; i++) {
291 size_t counter = k + m_kk*i;
292 m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
293 }
294 }
295
296 double dadt = da_dt();
297 double fac = dadt - m_a_current / (2.0 * TKelvin);
298 double vmb = mv - m_b_current;
299 double vpb = mv + m_b_current;
300 for (size_t k = 0; k < m_kk; k++) {
301 sbar[k] -=(GasConstant * log(GasConstant * TKelvin / (refP * mv))
303 + GasConstant * log(mv/vmb)
304 + GasConstant * b_vec_Curr_[k]/vmb
305 + m_pp[k]/(m_b_current * TKelvin * sqt) * log(vpb/mv)
306 - 2.0 * m_tmpV[k]/(m_b_current * sqt) * log(vpb/mv)
307 + b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) * fac
308 - 1.0 / (m_b_current * sqt) * b_vec_Curr_[k] / vpb * fac
309 );
310 }
311
313 getPartialMolarVolumes(m_partialMolarVolumes.data());
314 for (size_t k = 0; k < m_kk; k++) {
315 sbar[k] -= -m_partialMolarVolumes[k] * dpdT_;
316 }
317}
318
320{
321 // u_k = h_k - P * v_k
322 getPartialMolarVolumes(m_partialMolarVolumes.data());
324 double p = pressure();
325 for (size_t k = 0; k < nSpecies(); k++) {
326 ubar[k] -= p * m_partialMolarVolumes[k];
327 }
328}
329
331{
332 for (size_t k = 0; k < m_kk; k++) {
333 m_pp[k] = 0.0;
334 for (size_t i = 0; i < m_kk; i++) {
335 size_t counter = k + m_kk*i;
336 m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
337 }
338 }
339 for (size_t k = 0; k < m_kk; k++) {
340 m_tmpV[k] = 0.0;
341 for (size_t i = 0; i < m_kk; i++) {
342 size_t counter = k + m_kk*i;
343 m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
344 }
345 }
346
347 double sqt = sqrt(temperature());
348 double mv = molarVolume();
349 double vmb = mv - m_b_current;
350 double vpb = mv + m_b_current;
351 for (size_t k = 0; k < m_kk; k++) {
352 double num = (RT() + RT() * m_b_current/ vmb + RT() * b_vec_Curr_[k] / vmb
353 + RT() * m_b_current * b_vec_Curr_[k] /(vmb * vmb)
354 - 2.0 * m_pp[k] / (sqt * vpb)
355 + m_a_current * b_vec_Curr_[k] / (sqt * vpb * vpb)
356 );
357 double denom = (pressure() + RT() * m_b_current/(vmb * vmb) - m_a_current / (sqt * vpb * vpb)
358 );
359 vbar[k] = num / denom;
360 }
361}
362
363bool RedlichKwongMFTP::addSpecies(shared_ptr<Species> spec)
364{
365 bool added = MixtureFugacityTP::addSpecies(spec);
366 if (added) {
367 a_vec_Curr_.resize(m_kk * m_kk, 0.0);
368
369 // Initialize a_vec and b_vec to NaN, to screen for species with
370 // pureFluidParameters which are undefined in the input file:
371 b_vec_Curr_.push_back(NAN);
372 a_coeff_vec.resize(2, m_kk * m_kk, NAN);
373
374 m_pp.push_back(0.0);
375 m_coeffSource.push_back(CoeffSource::EoS);
376 m_partialMolarVolumes.push_back(0.0);
377 dpdni_.push_back(0.0);
378 }
379 return added;
380}
381
383{
384 // Contents of 'critical-properties.yaml', loaded later if needed
385 AnyMap critPropsDb;
386 std::unordered_map<string, AnyMap*> dbSpecies;
387
388 for (auto& [name, species] : m_species) {
389 auto& data = species->input;
390 size_t k = speciesIndex(name);
391 if (!isnan(a_coeff_vec(0, k + m_kk * k))) {
392 continue;
393 }
394 bool foundCoeffs = false;
395
396 if (data.hasKey("equation-of-state") &&
397 data["equation-of-state"].hasMapWhere("model", "Redlich-Kwong"))
398 {
399 // Read a and b coefficients from species 'input' information (that is, as
400 // specified in a YAML input file) specific to the Redlich-Kwong model
401 auto eos = data["equation-of-state"].getMapWhere(
402 "model", "Redlich-Kwong");
403
404 if (eos.hasKey("a") && eos.hasKey("b")) {
405 double a0 = 0, a1 = 0;
406 if (eos["a"].isScalar()) {
407 a0 = eos.convert("a", "Pa*m^6/kmol^2*K^0.5");
408 } else {
409 auto avec = eos["a"].asVector<AnyValue>(2);
410 a0 = eos.units().convert(avec[0], "Pa*m^6/kmol^2*K^0.5");
411 a1 = eos.units().convert(avec[1], "Pa*m^6/kmol^2/K^0.5");
412 }
413 double b = eos.convert("b", "m^3/kmol");
414 foundCoeffs = true;
415 setSpeciesCoeffs(name, a0, a1, b);
416 m_coeffSource[k] = CoeffSource::EoS;
417 }
418
419 if (eos.hasKey("binary-a")) {
420 AnyMap& binary_a = eos["binary-a"].as<AnyMap>();
421 const UnitSystem& units = binary_a.units();
422 for (auto& [name2, coeff] : binary_a) {
423 double a0 = 0, a1 = 0;
424 if (coeff.isScalar()) {
425 a0 = units.convert(coeff, "Pa*m^6/kmol^2*K^0.5");
426 } else {
427 auto avec = coeff.asVector<AnyValue>(2);
428 a0 = units.convert(avec[0], "Pa*m^6/kmol^2*K^0.5");
429 a1 = units.convert(avec[1], "Pa*m^6/kmol^2/K^0.5");
430 }
431 setBinaryCoeffs(name, name2, a0, a1);
432 }
433 }
434
435 if (foundCoeffs) {
436 continue;
437 }
438 }
439
440 // Coefficients have not been populated from model-specific input
441 double Tc = NAN, Pc = NAN;
442 if (data.hasKey("critical-parameters")) {
443 // Use critical state information stored in the species entry to
444 // calculate a and b
445 auto& critProps = data["critical-parameters"].as<AnyMap>();
446 Tc = critProps.convert("critical-temperature", "K");
447 Pc = critProps.convert("critical-pressure", "Pa");
448 m_coeffSource[k] = CoeffSource::CritProps;
449 } else {
450 // Search 'crit-properties.yaml' to find Tc and Pc. Load data if needed
451 if (critPropsDb.empty()) {
452 critPropsDb = AnyMap::fromYamlFile("critical-properties.yaml");
453 dbSpecies = critPropsDb["species"].asMap("name");
454 }
455
456 // All names in critical-properties.yaml are upper case
457 auto ucName = boost::algorithm::to_upper_copy(name);
458 if (dbSpecies.count(ucName)) {
459 auto& spec = *dbSpecies.at(ucName);
460 auto& critProps = spec["critical-parameters"].as<AnyMap>();
461 Tc = critProps.convert("critical-temperature", "K");
462 Pc = critProps.convert("critical-pressure", "Pa");
463 m_coeffSource[k] = CoeffSource::Database;
464 }
465 }
466
467 // Check if critical properties were found in either location
468 if (!isnan(Tc)) {
469 // Assuming no temperature dependence (that is, a1 = 0)
470 double a = omega_a * pow(GasConstant, 2) * pow(Tc, 2.5) / Pc;
471 double b = omega_b * GasConstant * Tc / Pc;
472 setSpeciesCoeffs(name, a, 0.0, b);
473 } else {
474 throw InputFileError("RedlichKwongMFTP::initThermo", data,
475 "No critical property or Redlich-Kwong parameters found "
476 "for species {}.", name);
477 }
478 }
479}
480
482 AnyMap& speciesNode) const
483{
485 size_t k = speciesIndex(name);
487 if (m_coeffSource[k] == CoeffSource::EoS) {
488 auto& eosNode = speciesNode["equation-of-state"].getMapWhere(
489 "model", "Redlich-Kwong", true);
490
491 size_t counter = k + m_kk * k;
492 if (a_coeff_vec(1, counter) != 0.0) {
493 vector<AnyValue> coeffs(2);
494 coeffs[0].setQuantity(a_coeff_vec(0, counter), "Pa*m^6/kmol^2*K^0.5");
495 coeffs[1].setQuantity(a_coeff_vec(1, counter), "Pa*m^6/kmol^2/K^0.5");
496 eosNode["a"] = std::move(coeffs);
497 } else {
498 eosNode["a"].setQuantity(a_coeff_vec(0, counter),
499 "Pa*m^6/kmol^2*K^0.5");
500 }
501 eosNode["b"].setQuantity(b_vec_Curr_[k], "m^3/kmol");
502 } else if (m_coeffSource[k] == CoeffSource::CritProps) {
503 auto& critProps = speciesNode["critical-parameters"];
504 double a = a_coeff_vec(0, k + m_kk * k);
505 double b = b_vec_Curr_[k];
506 double Tc = pow(a * omega_b / (b * omega_a * GasConstant), 2.0/3.0);
507 double Pc = omega_b * GasConstant * Tc / b;
508 critProps["critical-temperature"].setQuantity(Tc, "K");
509 critProps["critical-pressure"].setQuantity(Pc, "Pa");
510 }
511
512 if (m_binaryParameters.count(name)) {
513 auto& eosNode = speciesNode["equation-of-state"].getMapWhere(
514 "model", "Redlich-Kwong", true);
515 AnyMap bin_a;
516 for (const auto& [name2, coeffs] : m_binaryParameters.at(name)) {
517 if (coeffs.second == 0) {
518 bin_a[name2].setQuantity(coeffs.first, "Pa*m^6/kmol^2*K^0.5");
519 } else {
520 vector<AnyValue> C(2);
521 C[0].setQuantity(coeffs.first, "Pa*m^6/kmol^2*K^0.5");
522 C[1].setQuantity(coeffs.second, "Pa*m^6/kmol^2/K^0.5");
523 bin_a[name2] = std::move(C);
524 }
525 }
526 eosNode["binary-a"] = std::move(bin_a);
527 }
528}
529
531{
532 // note this agrees with tpx
533 double rho = density();
534 double mmw = meanMolecularWeight();
535 double molarV = mmw / rho;
536 double hh = m_b_current / molarV;
537 double zz = z();
538 double dadt = da_dt();
539 double T = temperature();
540 double sqT = sqrt(T);
541 double fac = dadt - m_a_current / (2.0 * T);
542 double sresid_mol_R = log(zz*(1.0 - hh)) + log(1.0 + hh) * fac / (sqT * GasConstant * m_b_current);
543 return GasConstant * sresid_mol_R;
544}
545
547{
548 // note this agrees with tpx
549 double rho = density();
550 double mmw = meanMolecularWeight();
551 double molarV = mmw / rho;
552 double hh = m_b_current / molarV;
553 double zz = z();
554 double dadt = da_dt();
555 double T = temperature();
556 double sqT = sqrt(T);
557 double fac = T * dadt - 3.0 *m_a_current / (2.0);
558 return GasConstant * T * (zz - 1.0) + fac * log(1.0 + hh) / (sqT * m_b_current);
559}
560
561double RedlichKwongMFTP::liquidVolEst(double TKelvin, double& presGuess) const
562{
563 double v = m_b_current * 1.1;
564 double atmp;
565 double btmp;
566 calculateAB(TKelvin, atmp, btmp);
567 double pres = std::max(psatEst(TKelvin), presGuess);
568 double Vroot[3];
569 bool foundLiq = false;
570 int m = 0;
571 while (m < 100 && !foundLiq) {
572 int nsol = solveCubic(TKelvin, pres, atmp, btmp, Vroot);
573 if (nsol == 1 || nsol == 2) {
574 double pc = critPressure();
575 if (pres > pc) {
576 foundLiq = true;
577 }
578 pres *= 1.04;
579 } else {
580 foundLiq = true;
581 }
582 }
583
584 if (foundLiq) {
585 v = Vroot[0];
586 presGuess = pres;
587 } else {
588 v = -1.0;
589 }
590 return v;
591}
592
593double RedlichKwongMFTP::densityCalc(double TKelvin, double presPa, int phaseRequested,
594 double rhoguess)
595{
596 // It's necessary to set the temperature so that m_a_current is set correctly.
597 setTemperature(TKelvin);
598 double tcrit = critTemperature();
599 double mmw = meanMolecularWeight();
600 if (rhoguess == -1.0) {
601 if (phaseRequested >= FLUID_LIQUID_0) {
602 double lqvol = liquidVolEst(TKelvin, presPa);
603 rhoguess = mmw / lqvol;
604 }
605 } else {
606 // Assume the Gas phase initial guess, if nothing is specified to the routine
607 rhoguess = presPa * mmw / (GasConstant * TKelvin);
608 }
609
610
611 double volguess = mmw / rhoguess;
612 NSolns_ = solveCubic(TKelvin, presPa, m_a_current, m_b_current, Vroot_);
613
614 double molarVolLast = Vroot_[0];
615 if (NSolns_ >= 2) {
616 if (phaseRequested >= FLUID_LIQUID_0) {
617 molarVolLast = Vroot_[0];
618 } else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
619 molarVolLast = Vroot_[2];
620 } else {
621 if (volguess > Vroot_[1]) {
622 molarVolLast = Vroot_[2];
623 } else {
624 molarVolLast = Vroot_[0];
625 }
626 }
627 } else if (NSolns_ == 1) {
628 if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT || phaseRequested == FLUID_UNDEFINED) {
629 molarVolLast = Vroot_[0];
630 } else {
631 return -2.0;
632 }
633 } else if (NSolns_ == -1) {
634 if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) {
635 molarVolLast = Vroot_[0];
636 } else if (TKelvin > tcrit) {
637 molarVolLast = Vroot_[0];
638 } else {
639 return -2.0;
640 }
641 } else {
642 molarVolLast = Vroot_[0];
643 return -1.0;
644 }
645 return mmw / molarVolLast;
646}
647
649{
650 double Vroot[3];
651 double T = temperature();
652 int nsol = solveCubic(T, pressure(), m_a_current, m_b_current, Vroot);
653 if (nsol != 3) {
654 return critDensity();
655 }
656
657 auto resid = [this, T](double v) {
658 double pp;
659 return dpdVCalc(T, v, pp);
660 };
661
662 boost::uintmax_t maxiter = 100;
663 auto [lower, upper] = bmt::toms748_solve(
664 resid, Vroot[0], Vroot[1], bmt::eps_tolerance<double>(48), maxiter);
665
666 double mmw = meanMolecularWeight();
667 return mmw / (0.5 * (lower + upper));
668}
669
671{
672 double Vroot[3];
673 double T = temperature();
674 int nsol = solveCubic(T, pressure(), m_a_current, m_b_current, Vroot);
675 if (nsol != 3) {
676 return critDensity();
677 }
678
679 auto resid = [this, T](double v) {
680 double pp;
681 return dpdVCalc(T, v, pp);
682 };
683
684 boost::uintmax_t maxiter = 100;
685 auto [lower, upper] = bmt::toms748_solve(
686 resid, Vroot[1], Vroot[2], bmt::eps_tolerance<double>(48), maxiter);
687
688 double mmw = meanMolecularWeight();
689 return mmw / (0.5 * (lower + upper));
690}
691
692double RedlichKwongMFTP::dpdVCalc(double TKelvin, double molarVol, double& presCalc) const
693{
694 double sqt = sqrt(TKelvin);
695 presCalc = GasConstant * TKelvin / (molarVol - m_b_current)
696 - m_a_current / (sqt * molarVol * (molarVol + m_b_current));
697
698 double vpb = molarVol + m_b_current;
699 double vmb = molarVol - m_b_current;
700 double dpdv = (- GasConstant * TKelvin / (vmb * vmb)
701 + m_a_current * (2 * molarVol + m_b_current) / (sqt * molarVol * molarVol * vpb * vpb));
702 return dpdv;
703}
704
706{
708 return -1 / (molarVolume() * dpdV_);
709}
710
712{
714 return -dpdT_ / (molarVolume() * dpdV_);
715}
716
718{
720 return molarVolume() * sqrt(-cp_mole() / cv_mole() * dpdV_ / meanMolecularWeight());
721}
722
724{
725 double TKelvin = temperature();
726 double mv = molarVolume();
727 double pres;
728
729 dpdV_ = dpdVCalc(TKelvin, mv, pres);
730 double sqt = sqrt(TKelvin);
731 double vpb = mv + m_b_current;
732 double vmb = mv - m_b_current;
733 double dadt = da_dt();
734 double fac = dadt - m_a_current/(2.0 * TKelvin);
735 dpdT_ = (GasConstant / vmb - fac / (sqt * mv * vpb));
736}
737
739{
740 double temp = temperature();
741 if (m_formTempParam == 1) {
742 for (size_t i = 0; i < m_kk; i++) {
743 for (size_t j = 0; j < m_kk; j++) {
744 size_t counter = i * m_kk + j;
745 a_vec_Curr_[counter] = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
746 }
747 }
748 }
749
750 m_b_current = 0.0;
751 m_a_current = 0.0;
752 for (size_t i = 0; i < m_kk; i++) {
753 m_b_current += moleFractions_[i] * b_vec_Curr_[i];
754 for (size_t j = 0; j < m_kk; j++) {
755 m_a_current += a_vec_Curr_[i * m_kk + j] * moleFractions_[i] * moleFractions_[j];
756 }
757 }
758 if (isnan(m_b_current)) {
759 // One or more species do not have specified coefficients.
760 fmt::memory_buffer b;
761 for (size_t k = 0; k < m_kk; k++) {
762 if (isnan(b_vec_Curr_[k])) {
763 if (b.size() > 0) {
764 fmt_append(b, ", {}", speciesName(k));
765 } else {
766 fmt_append(b, "{}", speciesName(k));
767 }
768 }
769 }
770 throw CanteraError("RedlichKwongMFTP::updateMixingExpressions",
771 "Missing Redlich-Kwong coefficients for species: {}", to_string(b));
772 }
773}
774
775void RedlichKwongMFTP::calculateAB(double temp, double& aCalc, double& bCalc) const
776{
777 bCalc = 0.0;
778 aCalc = 0.0;
779 if (m_formTempParam == 1) {
780 for (size_t i = 0; i < m_kk; i++) {
781 bCalc += moleFractions_[i] * b_vec_Curr_[i];
782 for (size_t j = 0; j < m_kk; j++) {
783 size_t counter = i * m_kk + j;
784 double a_vec_Curr = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
785 aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
786 }
787 }
788 } else {
789 for (size_t i = 0; i < m_kk; i++) {
790 bCalc += moleFractions_[i] * b_vec_Curr_[i];
791 for (size_t j = 0; j < m_kk; j++) {
792 size_t counter = i * m_kk + j;
793 double a_vec_Curr = a_coeff_vec(0,counter);
794 aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
795 }
796 }
797 }
798}
799
800double RedlichKwongMFTP::da_dt() const
801{
802 double dadT = 0.0;
803 if (m_formTempParam == 1) {
804 for (size_t i = 0; i < m_kk; i++) {
805 for (size_t j = 0; j < m_kk; j++) {
806 size_t counter = i * m_kk + j;
807 dadT+= a_coeff_vec(1,counter) * moleFractions_[i] * moleFractions_[j];
808 }
809 }
810 }
811 return dadT;
812}
813
814void RedlichKwongMFTP::calcCriticalConditions(double& pc, double& tc, double& vc) const
815{
816 double a0 = 0.0;
817 double aT = 0.0;
818 for (size_t i = 0; i < m_kk; i++) {
819 for (size_t j = 0; j <m_kk; j++) {
820 size_t counter = i + m_kk * j;
821 a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
822 aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
823 }
824 }
825 double a = m_a_current;
826 double b = m_b_current;
827 if (m_formTempParam != 0) {
828 a = a0;
829 }
830 if (b <= 0.0) {
831 tc = 1000000.;
832 pc = 1.0E13;
833 vc = omega_vc * GasConstant * tc / pc;
834 return;
835 }
836 if (a <= 0.0) {
837 tc = 0.0;
838 pc = 0.0;
839 vc = 2.0 * b;
840 return;
841 }
842 double tmp = a * omega_b / (b * omega_a * GasConstant);
843 double pp = 2./3.;
844 double sqrttc, f, dfdt, deltatc;
845
846 if (m_formTempParam == 0) {
847 tc = pow(tmp, pp);
848 } else {
849 tc = pow(tmp, pp);
850 for (int j = 0; j < 10; j++) {
851 sqrttc = sqrt(tc);
852 f = omega_a * b * GasConstant * tc * sqrttc / omega_b - aT * tc - a0;
853 dfdt = 1.5 * omega_a * b * GasConstant * sqrttc / omega_b - aT;
854 deltatc = - f / dfdt;
855 tc += deltatc;
856 }
857 if (deltatc > 0.1) {
858 throw CanteraError("RedlichKwongMFTP::calcCriticalConditions",
859 "didn't converge");
860 }
861 }
862
863 pc = omega_b * GasConstant * tc / b;
864 vc = omega_vc * GasConstant * tc / pc;
865}
866
867int RedlichKwongMFTP::solveCubic(double T, double pres, double a, double b, double Vroot[3]) const
868{
869
870 // Derive the coefficients of the cubic polynomial to solve.
871 double an = 1.0;
872 double bn = - GasConstant * T / pres;
873 double sqt = sqrt(T);
874 double cn = - (GasConstant * T * b / pres - a/(pres * sqt) + b * b);
875 double dn = - (a * b / (pres * sqt));
876
877 double tmp = a * omega_b / (b * omega_a * GasConstant);
878 double pp = 2./3.;
879 double tc = pow(tmp, pp);
880 double pc = omega_b * GasConstant * tc / b;
881 double vc = omega_vc * GasConstant * tc / pc;
882
883 int nSolnValues = MixtureFugacityTP::solveCubic(T, pres, a, b, a, Vroot, an, bn, cn, dn, tc, vc);
884
885 return nSolnValues;
886}
887
888}
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
const UnitSystem & units() const
Return the default units that should be used to convert stored values.
Definition AnyMap.h:630
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
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.
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:162
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:245
size_t m_kk
Number of species in the phase.
Definition Phase.h:947
double temperature() const
Temperature (K).
Definition Phase.h:662
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:760
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:151
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition Phase.cpp:138
double moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition Phase.cpp:545
virtual double density() const
Density (kg/m^3).
Definition Phase.h:687
double mean_X(const double *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition Phase.cpp:737
shared_ptr< Species > species(const string &name) const
Return the Species object for the named species.
Definition Phase.cpp:977
virtual double molarVolume() const
Molar volume (m^3/kmol).
Definition Phase.cpp:702
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 getChemPotentials_RT(double *mu) const override
Get the array of non-dimensional species chemical potentials.
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.
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...
virtual double refPressure() const
Returns the reference pressure in Pa.
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:195
const double SmallNumber
smallest number to compare to zero.
Definition ct_defs.h:158
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
Definition AnyMap.cpp:1926
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...